MinuitFitter.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/IFitVar.h"
4 #include "CAFAna/Core/ISyst.h"
8 
9 #include "OscLib/IOscCalc.h"
10 #include "Utilities/func/MathUtil.h"
11 
12 #include "TError.h"
13 
14 #include "Minuit2/FCNBase.h"
15 #include "Minuit2/MnMigrad.h"
16 #include "Minuit2/MnUserParameters.h"
17 #include "Minuit2/FunctionMinimum.h"
18 
19 #include <cassert>
20 #include <iostream>
21 #include <memory>
22 #include <vector>
23 
24 namespace ana
25 {
26  // Minuit calls this at some point, and it creates a static. Which doesn't
27  // like happening on multiple threads at once. So do it upfront while we're
28  // still single-threaded.
30  {
33 
34  //----------------------------------------------------------------------
35  MinuitFitter::MinuitFitSummary::MinuitFitSummary(std::unique_ptr<ROOT::Minuit2::FunctionMinimum> && minimum)
36  : fMinimum(std::move(minimum))
37  {
38  assert (fMinimum.get() && "Can't make a MinuitFitSummary from a null ROOT::Minuit2::FunctionMinimum");
39  }
40 
41  //----------------------------------------------------------------------
43  {
44  if (!other)
45  return true;
46 
47  if (const auto mnFitSummary = dynamic_cast<const MinuitFitSummary*>(other))
48  return EvalMetricVal() < mnFitSummary->EvalMetricVal();
49  else
50  {
51  assert(false && "Can't compare a MinuitFitSummary to another kind of IFitSummary");
52  return false; // prevent compiler warnings
53  }
54  }
55 
56  //----------------------------------------------------------------------
58  {
59  return fMinimum->Fval();
60  }
61 
62 
63  //----------------------------------------------------------------------
65  std::vector<const IFitVar *> vars,
66  std::vector<const ISyst *> systs,
67  FitOpts opts)
68  : IFitter(vars, systs),
69  fExpt(expt),
70  fFitOpts(opts),
72  {
73  }
74 
75  //----------------------------------------------------------------------
77  {
78  }
79 
80  //----------------------------------------------------------------------
82  {
83  return false;
84  /*
85  // Make completely opt-in for now
86  if (getenv("CAFANA_ANALYTIC_DERIVATIVES") == 0) return false;
87 
88  // No point using derivatives for FitVars only, we do finite differences,
89  // probably worse than MINUIT would.
90  if (fSysts.empty()) return false;
91 
92  // Otherwise, do the minimal trial to see if the experiment will return a
93  // gradient.
94  std::unordered_map<const ISyst *, double> dchi = {{fSysts[0], 0}};
95  osc::NoOscillations calc;
96  fExpt->Derivative(&calc, SystShifts::Nominal(), dchi);
97  return !dchi.empty();
98  */
99  }
100 
101  //----------------------------------------------------------------------
102  std::unique_ptr<IFitter::IFitSummary>
104  SystShifts &systSeed,
105  Verbosity verb) const
106  {
107  fVerb = verb;
108  if((fFitOpts & kPrefitSysts) && !fVars.empty()){
109  // Update systSeed to the best syst values, while holding the oscillation
110  // parameters fixed
111  MinuitFitter prefit(fExpt, {}, fSysts, fFitOpts);
112  prefit.Fit(seed, systSeed, verb);
113  // Then continue with the full fit as usual
114  }
115 
116  fCalc = seed;
117  *fShifts = systSeed;
118 
119  ROOT::Minuit2::MnUserParameters mnPars;
120 
121  for (const IFitVar *v: fVars)
122  {
123  const double val = v->GetValue(seed);
124  // name, value, error
125  mnPars.Add(v->ShortName(), val, val ? val / 2 : .1);
126  }
127  // One way this can go wrong is if two variables have the same ShortName
128  assert(mnPars.Params().size() == fVars.size());
129  for (const ISyst *s: fSysts)
130  {
131  const double val = systSeed.GetShift(s);
132  // name, value, error
133  mnPars.Add(s->ShortName(), val, 1);
134  }
135  // One way this can go wrong is if two variables have the same ShortName
136  assert(mnPars.Params().size() == fVars.size() + fSysts.size());
137 
138  ROOT::Minuit2::MnApplication *mnApp = 0;
139 
140  if ((fFitOpts & kPrecisionMask) == kGradDesc)
141  {
142  mnApp = new GradientDescent(*this, mnPars);
143  } else
144  {
146  {
147  mnApp = new ROOT::Minuit2::MnMigrad(*this, mnPars, GetMINUITPrecision());
148  } else
149  {
150  mnApp = new ROOT::Minuit2::MnMigrad(*((ROOT::Minuit2::FCNBase *) this), mnPars, GetMINUITPrecision());
151  }
152  }
153 
154  // Minuit2 doesn't give a good way to control verbosity...
155  const int olderr = gErrorIgnoreLevel;
156  if (fVerb <= Verbosity::kQuiet) gErrorIgnoreLevel = 1001; // Ignore warnings
157  auto minpt = std::make_unique<ROOT::Minuit2::FunctionMinimum>((*mnApp)());
158 
159  gErrorIgnoreLevel = olderr;
160 
161  if(minpt->IsValid()){
162  const std::vector<double> cov = minpt->UserCovariance().Data();
163  const unsigned int nPars = mnPars.Params().size();
164 
165  assert(cov.size() == (nPars*(nPars+1))/2);
166 
167  fCovariance.ResizeTo(nPars, nPars);
168  int pos = 0;
169  for(unsigned int iRow = 0; iRow < nPars; iRow++){
170  for(unsigned int iCol = 0; iCol < iRow+1; iCol++){
171  fCovariance(iRow,iCol) = fCovariance(iCol, iRow) = cov[pos];
172  pos++;
173  }
174  }
175  }
176  else{
177  fCovariance.ResizeTo(0, 0);
178  if(fVerb > Verbosity::kQuiet)
179  std::cout << "*** ERROR: minimum is not valid ***" << std::endl;
180  }
181 
182  const std::vector<double> minvec = minpt->UserParameters().Params();
183 
184  // Store results back to the "seed" variable
185  for (unsigned int i = 0; i < fVars.size(); ++i)
186  {
187  fVars[i]->SetValue(seed, minvec[i]);
188  }
189  // Store systematic results back into "systSeed"
190  for (unsigned int j = 0; j < fSysts.size(); ++j)
191  {
192  systSeed.SetShift(fSysts[j], minvec[fVars.size() + j]);
193  }
194 
195  delete mnApp;
196 
197  return std::make_unique<MinuitFitSummary>(std::move(minpt));
198  }
199 
200  //----------------------------------------------------------------------
202  {
203  if ((opts & kPrecisionMask) == kGradDesc &&
204  !fSysts.empty() &&
206  {
207  std::cout
208  << "Warning - not setting precision to kGradDesc, since analytic gradients are not supported by this experiment"
209  << std::endl;
210  return;
211  }
212 
213  fFitOpts = opts;
214  }
215 
216  //----------------------------------------------------------------------
217  double MinuitFitter::operator()(const std::vector<double> &pars) const
218  {
219  ++fNEval;
220 
221  assert(pars.size() == fVars.size() + fSysts.size());
222 
223  DecodePars(pars, fCalc); // Updates fCalc and fShifts
224 
225  // Have to re-fetch the FitVar values because DecodePars() will have
226  // truncated values to the physical range where necessary.
227  double penalty = 0;
228  for (unsigned int i = 0; i < fVars.size(); ++i)
229  {
230  penalty += fVars[i]->Penalty(pars[i], fCalc);
231  }
232 
233  return fExpt->ChiSq(fCalc, *fShifts) + penalty + fShifts->Penalty();
234  }
235 
236  //----------------------------------------------------------------------
237  std::vector<double> MinuitFitter::Gradient(const std::vector<double> &pars) const
238  {
239  abort();
240  /*
241  ++fNEvalGrad;
242 
243  std::vector<double> ret(pars.size());
244 
245  // TODO handling of FitVars including penalty terms
246 
247  if (!fVars.empty())
248  {
249  // Have to use finite differences to calculate these derivatives
250  const double dx = 1e-9;
251  const double nom = (*this)(pars);
252  ++fNEvalFiniteDiff;
253  std::vector<double> parsCopy = pars;
254  for (unsigned int i = 0; i < fVars.size(); ++i)
255  {
256  parsCopy[i] += dx;
257  ret[i] = ((*this)(parsCopy) - nom) / dx;
258  ++fNEvalFiniteDiff;
259  parsCopy[i] = pars[i];
260  }
261  }
262 
263  // Get systematic parts analytically
264 
265  DecodePars(pars, fCalc); // Updates fCalc and fShifts
266 
267  std::unordered_map<const ISyst *, double> dchi;
268  for (const ISyst *s: fSysts) dchi[s] = 0;
269  fExpt->Derivative(fCalc, *fShifts, dchi);
270 
271  for (unsigned int i = 0; i < fSysts.size(); ++i)
272  {
273  // Include the derivative of the penalty terms too. TODO this is only
274  // right for quadratic penalties (ie all the currently existing ones)
275  ret[fVars.size() + i] = dchi[fSysts[i]] + 2 * fShifts->GetShift(fSysts[i]);
276  }
277 
278  return ret;
279  */
280  }
281 
282  //----------------------------------------------------------------------
283  void MinuitFitter::DecodePars(const std::vector<double> &pars, osc::IOscCalcAdjustable *calc) const
284  {
285  assert(pars.size() == fVars.size()+fSysts.size());
286 
287  if (fVars.size() > 0)
288  {
289  assert(calc);
290  for(unsigned int i = 0; i < fVars.size(); ++i){
291  auto val = pars[i];
292  fVars[i]->SetValue(calc, val);
293  }
294  }
295 
296  for(unsigned int j = 0; j < fSysts.size(); ++j){
297  auto val = pars[fVars.size()+j];
298  fShifts->SetShift(fSysts[j], val);
299  }
300  }
301 }
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
Definition: IExperiment.h:18
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
struct ana::MinuitStaticInit gMinuitStaticInit
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
bool SupportsDerivatives() const
MinuitFitSummary(std::unique_ptr< ROOT::Minuit2::FunctionMinimum > &&minimum)
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
osc::IOscCalcAdjustable * fCalc
Definition: MinuitFitter.h:88
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const IExperiment * fExpt
Definition: MinuitFitter.h:89
const XML_Char * s
Definition: expat.h:262
expt
Definition: demo5.py:34
unsigned int seed
Definition: runWimpSim.h:102
std::vector< double > Gradient(const std::vector< double > &pars) const override
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
bool IsBetterThan(const IFitSummary *other) const override
Is the current fit better than other?
int GetMINUITPrecision() const
Definition: MinuitFitter.h:58
T GetShift(const ISyst *syst) const
void DecodePars(const std::vector< double > &pars, osc::IOscCalcAdjustable *calc) const
Stuff the parameters into the calculator and/or syst shifts object.
const std::map< std::pair< std::string, std::string >, Variable > vars
Base class for fitters.
Definition: IFitter.h:16
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
TMatrixDSym fCovariance
Definition: MinuitFitter.h:98
std::unique_ptr< IFitSummary > FitHelperSeeded(osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const override
Helper for FitHelper.
double EvalMetricVal() const override
The chi^2 passed to the fitter.
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
std::string pars("Th23Dmsq32")
MinuitFitter(const IExperiment *expt, std::vector< const IFitVar * > vars, std::vector< const ISyst * > systs={}, FitOpts opts=kNormal)
Base class defining interface for experiments.
Definition: IExperiment.h:14
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
virtual double operator()(const std::vector< double > &pars) const override
Evaluate the log-likelihood, as required by MINUT interface.
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
Verbosity fVerb
Definition: IFitter.h:117
std::unique_ptr< ROOT::Minuit2::FunctionMinimum > fMinimum
the minimizer instance
Definition: MinuitFitter.h:76
void SetFitOpts(FitOpts opts)
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17