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  fExpt->Reset();
119 
120  ROOT::Minuit2::MnUserParameters mnPars;
121 
122  for (const IFitVar *v: fVars)
123  {
124  const double val = v->GetValue(seed);
125  // name, value, error
126  mnPars.Add(v->ShortName(), val, val ? val / 2 : .1);
127  }
128  // One way this can go wrong is if two variables have the same ShortName
129  assert(mnPars.Params().size() == fVars.size());
130  for (const ISyst *s: fSysts)
131  {
132  const double val = systSeed.GetShift(s);
133  // name, value, error
134  mnPars.Add(s->ShortName(), val, 1);
135  }
136  // One way this can go wrong is if two variables have the same ShortName
137  assert(mnPars.Params().size() == fVars.size() + fSysts.size());
138 
139  ROOT::Minuit2::MnApplication *mnApp = 0;
140 
141  if ((fFitOpts & kPrecisionMask) == kGradDesc)
142  {
143  mnApp = new GradientDescent(*this, mnPars);
144  } else
145  {
147  {
148  mnApp = new ROOT::Minuit2::MnMigrad(*this, mnPars, GetMINUITPrecision());
149  } else
150  {
151  mnApp = new ROOT::Minuit2::MnMigrad(*((ROOT::Minuit2::FCNBase *) this), mnPars, GetMINUITPrecision());
152  }
153  }
154 
155  // Minuit2 doesn't give a good way to control verbosity...
156  const int olderr = gErrorIgnoreLevel;
157  if (fVerb <= Verbosity::kQuiet) gErrorIgnoreLevel = 1001; // Ignore warnings
158  auto minpt = std::make_unique<ROOT::Minuit2::FunctionMinimum>((*mnApp)());
159 
160  gErrorIgnoreLevel = olderr;
161 
162  if(minpt->IsValid()){
163  const std::vector<double> cov = minpt->UserCovariance().Data();
164  const unsigned int nPars = mnPars.Params().size();
165 
166  assert(cov.size() == (nPars*(nPars+1))/2);
167 
168  fCovariance.ResizeTo(nPars, nPars);
169  int pos = 0;
170  for(unsigned int iRow = 0; iRow < nPars; iRow++){
171  for(unsigned int iCol = 0; iCol < iRow+1; iCol++){
172  fCovariance(iRow,iCol) = fCovariance(iCol, iRow) = cov[pos];
173  pos++;
174  }
175  }
176  }
177  else{
178  fCovariance.ResizeTo(0, 0);
179  if(fVerb > Verbosity::kQuiet)
180  std::cout << "*** ERROR: minimum is not valid ***" << std::endl;
181  }
182 
183  const std::vector<double> minvec = minpt->UserParameters().Params();
184 
185  // Store results back to the "seed" variable
186  for (unsigned int i = 0; i < fVars.size(); ++i)
187  {
188  fVars[i]->SetValue(seed, minvec[i]);
189  }
190  // Store systematic results back into "systSeed"
191  for (unsigned int j = 0; j < fSysts.size(); ++j)
192  {
193  systSeed.SetShift(fSysts[j], minvec[fVars.size() + j]);
194  }
195 
196  delete mnApp;
197 
198  return std::make_unique<MinuitFitSummary>(std::move(minpt));
199  }
200 
201  //----------------------------------------------------------------------
203  {
204  if ((opts & kPrecisionMask) == kGradDesc &&
205  !fSysts.empty() &&
207  {
208  std::cout
209  << "Warning - not setting precision to kGradDesc, since analytic gradients are not supported by this experiment"
210  << std::endl;
211  return;
212  }
213 
214  fFitOpts = opts;
215  }
216 
217  //----------------------------------------------------------------------
218  double MinuitFitter::operator()(const std::vector<double> &pars) const
219  {
220  ++fNEval;
221 
222  assert(pars.size() == fVars.size() + fSysts.size());
223 
224  DecodePars(pars, fCalc); // Updates fCalc and fShifts
225 
226  // Have to re-fetch the FitVar values because DecodePars() will have
227  // truncated values to the physical range where necessary.
228  double penalty = 0;
229  for (unsigned int i = 0; i < fVars.size(); ++i)
230  {
231  penalty += fVars[i]->Penalty(pars[i], fCalc);
232  }
233 
234  return fExpt->ChiSq(fCalc, *fShifts) + penalty + fShifts->Penalty();
235  }
236 
237  //----------------------------------------------------------------------
238  std::vector<double> MinuitFitter::Gradient(const std::vector<double> &pars) const
239  {
240  abort();
241  /*
242  ++fNEvalGrad;
243 
244  std::vector<double> ret(pars.size());
245 
246  // TODO handling of FitVars including penalty terms
247 
248  if (!fVars.empty())
249  {
250  // Have to use finite differences to calculate these derivatives
251  const double dx = 1e-9;
252  const double nom = (*this)(pars);
253  ++fNEvalFiniteDiff;
254  std::vector<double> parsCopy = pars;
255  for (unsigned int i = 0; i < fVars.size(); ++i)
256  {
257  parsCopy[i] += dx;
258  ret[i] = ((*this)(parsCopy) - nom) / dx;
259  ++fNEvalFiniteDiff;
260  parsCopy[i] = pars[i];
261  }
262  }
263 
264  // Get systematic parts analytically
265 
266  DecodePars(pars, fCalc); // Updates fCalc and fShifts
267 
268  std::unordered_map<const ISyst *, double> dchi;
269  for (const ISyst *s: fSysts) dchi[s] = 0;
270  fExpt->Derivative(fCalc, *fShifts, dchi);
271 
272  for (unsigned int i = 0; i < fSysts.size(); ++i)
273  {
274  // Include the derivative of the penalty terms too. TODO this is only
275  // right for quadratic penalties (ie all the currently existing ones)
276  ret[fVars.size() + i] = dchi[fSysts[i]] + 2 * fShifts->GetShift(fSysts[i]);
277  }
278 
279  return ret;
280  */
281  }
282 
283  //----------------------------------------------------------------------
284  void MinuitFitter::DecodePars(const std::vector<double> &pars, osc::IOscCalcAdjustable *calc) const
285  {
286  assert(pars.size() == fVars.size()+fSysts.size());
287 
288  if (fVars.size() > 0)
289  {
290  assert(calc);
291  for(unsigned int i = 0; i < fVars.size(); ++i){
292  auto val = pars[i];
293  fVars[i]->SetValue(calc, val);
294  }
295  }
296 
297  for(unsigned int j = 0; j < fSysts.size(); ++j){
298  auto val = pars[fVars.size()+j];
299  fShifts->SetShift(fSysts[j], val);
300  }
301  }
302 }
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:17
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.
virtual void Reset() const
Definition: IExperiment.h:32
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
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