IFitter.cxx
Go to the documentation of this file.
1 #include "CAFAna/Fit/IFitter.h"
2 
3 #include "CAFAna/Core/IFitVar.h"
4 #include "CAFAna/Core/ISyst.h"
7 
8 #include "OscLib/IOscCalc.h"
9 
10 #include "Utilities/func/MathUtil.h"
12 
13 #include <cassert>
14 
15 namespace ana
16 {
17 
18  SystShifts IFitter::junkShifts = SystShifts();
19 
20  //----------------------------------------------------------------------
21  IFitter::IFitter(const std::vector<const IFitVar*>& vars,
22  const std::vector<const ISyst*>& systs)
23  : fVars(vars), fSysts(systs)
24  {
25  }
26 
27  //----------------------------------------------------------------------
29  : fVars(other.fVars),
30  fSysts(other.fSysts),
31  fShifts(std::make_unique<SystShifts>(*other.fShifts))
32  {}
33 
34  //----------------------------------------------------------------------
36  {
37  }
38 
39  //----------------------------------------------------------------------
41  {
42  }
43 
44  //----------------------------------------------------------------------
45  std::vector<IFitter::SeedPt>
47  const std::vector<SystShifts>& systSeedPts) const
48  {
49  std::vector<SeedPt> ret;
50  for(Seed seed: seedPts.GetSeeds()) ret.emplace_back(seed, SystShifts::Nominal());
51 
52  // Now duplicate as many times as required for the syst seeds
53  if(!systSeedPts.empty()){
54  std::vector<SeedPt> newret;
55  for(const SystShifts& s: systSeedPts){
56  for(SeedPt pt: ret){
57  pt.shift = s;
58  newret.push_back(pt);
59  }
60  }
61  ret = newret;
62  }
63 
64  return ret;
65  }
66 
67  //----------------------------------------------------------------------
68  std::unique_ptr<IFitter::IFitSummary>
70  SystShifts& bestSysts,
71  const SeedList& seedPts,
72  const std::vector<SystShifts> &systSeedPts,
73  Verbosity verb) const
74  {
75  fVerb = verb;
76  ValidateSeeds(seed, seedPts, systSeedPts);
77 
78  if (fVerb == kVerbose)
79  {
80  std::cout << "Finding best fit for";
81  for (const auto *v: fVars) std::cout << " " << v->ShortName();
82  for (const auto *s: fSysts) std::cout << " " << s->ShortName();
83 // if (fSupportsDerivatives) std::cout << " using analytic derivatives";
84  std::cout << "..." << std::endl;
85  }
86 
87  // Do all the actual work. This wrapper function is just so we can have
88  // better control of printing.
89  auto fitSummary = FitHelper(seed, bestSysts, seedPts, systSeedPts, verb);
90 
91  if (fVerb == kVerbose)
92  {
93  std::cout << "Best fit";
94  for (const auto *v: fVars)
95  {
96  std::cout << ", " << v->ShortName() << " = " << v->GetValue(seed);
97  }
98  for (const auto *s: fSysts)
99  {
100  std::cout << ", " << s->ShortName() << " = " << bestSysts.GetShift(s);
101  }
102  std::cout << ", eval metric = " << fitSummary->EvalMetricVal() << std::endl;
103 
104 // std::cout << " found with " << fNEval << " evaluations of the likelihood";
105 // if (fNEvalFiniteDiff > 0)
106 // std::cout << " (" << fNEvalFiniteDiff << " to evaluate gradients numerically)";
107 // if (fNEvalGrad > 0)
108 // std::cout << " and " << fNEvalGrad << " evaluations of the gradient";
109 // std::cout << std::endl;
110  }
111 
112  return fitSummary;
113  }
114 
115  //----------------------------------------------------------------------
116  std::unique_ptr<IFitter::IFitSummary>
118  SystShifts& bestSysts,
119  const SeedList& seedPts,
120  const std::vector<SystShifts>& systSeedPts,
121  Verbosity verb) const
122  {
123  // if user passed a derived kind of SystShifts, this preserves it
124  fShifts = bestSysts.Copy();
125  fShifts->ResetToNominal();
126 
127  const std::vector<SeedPt> pts = ExpandSeeds(seedPts, systSeedPts);
128 
129  std::unique_ptr<IFitSummary> bestFitSummary;
130  std::vector<double> bestFitPars, bestSystPars;
131 
132  for (const SeedPt &pt: pts)
133  {
134  osc::IOscCalcAdjustable* seed = nullptr;
135  if(initseed) seed = initseed->Copy();
136 
137  pt.fitvars.ResetCalc(seed);
138 
139  // be sure to keep any derived class stuff around
140  auto shift = fShifts->Copy();
141  *shift = pt.shift;
142  auto fitSummary = FitHelperSeeded(seed, *shift, verb);
143  if (fitSummary->IsBetterThan(bestFitSummary.get()))
144  {
145  bestFitSummary = std::move(fitSummary);
146  // Store the best fit values of all the parameters we know are being
147  // varied.
148  bestFitPars.clear();
149  bestSystPars.clear();
150  for (const auto *v: fVars)
151  bestFitPars.push_back(v->GetValue(seed));
152  for (const auto *s: fSysts)
153  bestSystPars.push_back(shift->GetShift<double>(s));
154  }
155 
156  delete seed;
157  } // end for pt
158 
159  assert(bestFitSummary);
160  // Stuff the results of the actual best fit back into the seeds
161  for (unsigned int i = 0; i < fVars.size(); ++i)
162  fVars[i]->SetValue(initseed, bestFitPars[i]);
163  for (unsigned int i = 0; i < fSysts.size(); ++i)
164  bestSysts.SetShift(fSysts[i], bestSystPars[i]);
165 
166  return std::move(bestFitSummary);
167  }
168 
169  //----------------------------------------------------------------------
171  const SeedList& seedPts,
172  const std::vector<SystShifts>& systSeedPts) const
173  {
174  if(!seed && !fVars.empty()){
175  std::cout << "ERROR: MinuitFitter::Fit() trying to fit oscillation parameters without an oscillation calculator" << std::endl;
176  abort();
177  }
178 
179  for(const IFitVar* v: seedPts.ActiveFitVars()){
180  if (std::find(fVars.begin(), fVars.end(), v) == fVars.end()){
181  std::cout << "ERROR MinuitFitter::Fit() trying to seed '"
182  << v->ShortName()
183  << "' which is not part of the fit." << std::endl;
184  abort();
185  }
186  }
187 
188  for(const SystShifts& pt: systSeedPts){
189  for(const ISyst* s: pt.ActiveSysts()){
190  if(std::find(fSysts.begin(), fSysts.end(), s) == fSysts.end()){
191  std::cout << "ERROR MinuitFitter::Fit() trying to seed '"
192  << s->ShortName()
193  << "' which is not part of the fit." << std::endl;
194  abort();
195  }
196  }
197  }
198  }
199 
200 } // namespace ana
virtual _IOscCalcAdjustable< T > * Copy() const =0
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual std::unique_ptr< SystShifts > Copy() const
Definition: SystShifts.cxx:67
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const std::vector< Seed > & GetSeeds() const
Definition: SeedList.h:47
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void ValidateSeeds(osc::IOscCalcAdjustable *seed, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) const
Check that the seeds that were specified are compatible with the vars being fitted.
Definition: IFitter.cxx:170
static SystShifts Nominal()
Definition: SystShifts.h:34
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::vector< const ISyst * > fSysts
Definition: IFitter.h:115
const XML_Char * s
Definition: expat.h:262
virtual ~IFitter()
Definition: IFitter.cxx:35
unsigned int seed
Definition: runWimpSim.h:102
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
T GetShift(const ISyst *syst) const
virtual std::unique_ptr< IFitSummary > FitHelperSeeded(osc::IOscCalcAdjustable *seed, SystShifts &systSeed, Verbosity verb) const =0
Helper for FitHelper – actually does fitting. Reimplement in derived classes.
const std::map< std::pair< std::string, std::string >, Variable > vars
Base class for fitters.
Definition: IFitter.h:16
OStream cout
Definition: OStream.cxx:6
std::vector< const IFitVar * > fVars
Definition: IFitter.h:114
static SystShifts junkShifts
Definition: IFitter.h:19
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::vector< SeedPt > ExpandSeeds(const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts) const
Definition: IFitter.cxx:46
std::unique_ptr< SystShifts > fShifts
Definition: IFitter.h:116
std::set< const IFitVar * > ActiveFitVars() const
Definition: SeedList.cxx:53
virtual std::unique_ptr< IFitSummary > FitHelper(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, Verbosity verb) const
Helper for Fit.
Definition: IFitter.cxx:117
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
Verbosity fVerb
Definition: IFitter.h:117
IFitter(const std::vector< const IFitVar * > &vars, const std::vector< const ISyst * > &systs={})
Definition: IFitter.cxx:21