Fit.cxx
Go to the documentation of this file.
1 #include "CAFAna/Fit/Fit.h"
3 
4 #include "CAFAna/Core/Progress.h"
5 #include "CAFAna/Core/IFitVar.h"
8 
9 #include "OscLib/IOscCalc.h"
10 #include "Utilities/func/MathUtil.h"
11 
12 #include "TError.h"
13 #include "TGraph.h"
14 #include "TH1.h"
15 
16 #include <cassert>
17 #include <memory>
18 #include <iostream>
19 
20 namespace ana
21 {
22  //----------------------------------------------------------------------
23  double SimpleFOM(const Spectrum& obs, const Spectrum& unosc, double pot)
24  {
25  if(pot == 0) pot = obs.POT();
26 
27  std::unique_ptr<TH1> oh(obs.ToTH1(pot));
28  std::unique_ptr<TH1> bh(unosc.ToTH1(pot));
29  assert(oh->GetNbinsX() == bh->GetNbinsX());
30 
31  double fomSq = 0;
32 
33  // Combine s/sqrt(s+b) in quadrature between bins
34  for(int i = 0; i < oh->GetNbinsX(); ++i){
35  const double o = oh->GetBinContent(i);
36  const double b = bh->GetBinContent(i);
37  const double s = o-b;
38 
39  if(s <= 0) continue;
40 
41  fomSq += util::sqr(s)/(s+b);
42  }
43 
44  return sqrt(fomSq);
45  }
46 
47  //----------------------------------------------------------------------
48  TGraph* Profile(const IExperiment* expt,
50  int nbinsx, double minx, double maxx,
51  double input_minchi,
52  const std::vector<const IFitVar*>& profVars,
53  const std::vector<const ISyst*>& profSysts,
54  const SeedList& seedPts,
55  const std::vector<SystShifts>& systSeedPts,
56  std::map<const IFitVar*, TGraph*>& profVarsMap,
57  std::map<const ISyst*, TGraph*>& profSystsMap,
59  {
60  Progress prog ("Filling profile");
61  // If we're called with the default arguments they could already have stuff
62  // in from before.
63  for(auto it: profVarsMap) delete it.second;
64  for(auto it: profSystsMap) delete it.second;
65  profVarsMap.clear();
66  profSystsMap.clear();
67 
68  // And then create the plots we'll be filling
69  for(const IFitVar* v: profVars) profVarsMap[v] = new TGraph;
70  for(const ISyst* s: profSysts) profSystsMap[s] = new TGraph;
71 
72  TGraph* ret = new TGraph;
73 
74  // Save the values of the fit vars as they were in the seed so we can put
75  // them back to that value every iteration.
76  std::vector<double> seedValues;
77  for(const IFitVar* v: profVars) seedValues.push_back(v->GetValue(calc));
78 
79  double minpos = 0;
80  double minchi = 1e10;
81 
82  MinuitFitter fit(expt, profVars, profSysts, opts);
83 
84  for(int n = 0; n <= nbinsx; ++n){
85  prog.SetProgress((double) n/(nbinsx+1));
86 
87  const double x = minx + (maxx-minx)*double(n)/nbinsx;
88  v->SetValue(calc, x);
89 
90  // Put oscillation values back to their seed position each iteration
91  for(unsigned int i = 0; i < seedValues.size(); ++i)
92  profVars[i]->SetValue( calc, seedValues[i] );
93 
94  SystShifts systshift = SystShifts::Nominal();
95  const double chi = fit.Fit(calc, systshift, seedPts, systSeedPts, MinuitFitter::kQuiet)->EvalMetricVal();
96 
97  ret->SetPoint(ret->GetN(), x, chi);
98 
99  if(chi < minchi){
100  minchi = chi;
101  minpos = x;
102  }
103  for(const IFitVar* var: profVars){
104  profVarsMap[var]->SetPoint(n, x, var->GetValue(calc));
105  }
106  for(const ISyst* s: profSysts){
107  profSystsMap[s]->SetPoint(n, x, systshift.GetShift(s));
108  }
109  }
110  prog.Done();
111  // If we weren't given an explicit minimum chisq, go find one
112  if(input_minchi == -1){
113  std::vector<const IFitVar*> allVars = {v};
114  for(unsigned int i = 0; i < seedValues.size(); ++i) {
115  profVars[i]->SetValue(calc, seedValues[i]);
116  allVars.push_back(profVars[i]);
117  }
118  MinuitFitter fit(expt, allVars, profSysts, opts);
119  // Seed from best grid point
120  v->SetValue(calc, minpos);
121  minchi = fit.Fit(calc)->EvalMetricVal(); // get a better value
122 
123  // Add in the best fit point explicitly so that the graph goes to zero
124  ret->SetPoint(ret->GetN(), v->GetValue(calc), minchi);
125  // And put it in the right place in the order (default is sort by X)
126  ret->Sort();
127  }
128  else{
129  minchi = input_minchi;
130  }
131 
132  // Zero-subtract the result graph
133  for(int i = 0; i < ret->GetN(); ++i){
134  double x, y;
135  ret->GetPoint(i, x, y);
136  ret->SetPoint(i, x, y-minchi);
137  }
138 
139  return ret;
140  }
141 
142  //----------------------------------------------------------------------
143  TGraph* SqrtProfile(const IExperiment* expt,
145  int nbinsx, double minx, double maxx, double minchi,
146  std::vector<const IFitVar*> profVars,
147  std::vector<const ISyst*> profSysts,
148  const SeedList& seedPts,
149  const std::vector<SystShifts>& systSeedPts,
150  std::map<const IFitVar*, TGraph*>& profVarsMap,
151  std::map<const ISyst*, TGraph*>& systsMap,
153  {
154  TGraph* ret = Profile(expt, calc,
155  v, nbinsx, minx, maxx,
156  minchi, profVars, profSysts, seedPts, systSeedPts,
157  profVarsMap, systsMap, opts);
158 
159  for(int i = 0; i < ret->GetN(); ++i){
160  double x, y;
161  ret->GetPoint(i, x, y);
162  ret->SetPoint(i, x, y > 0 ? sqrt(y) : 0);
163  }
164 
165  return ret;
166  }
167 
168  //----------------------------------------------------------------------
169  TGraph* Slice(const IExperiment* expt,
171  int nbinsx, double minx, double maxx,
172  double minchi,
174  {
175  return Profile(expt, calc, v, nbinsx, minx, maxx, minchi,
176  {}, {}, {}, {}, empty_vars_map, empty_syst_map, opts);
177  }
178 
179  //----------------------------------------------------------------------
180  TGraph* SqrtSlice(const IExperiment* expt,
182  int nbinsx, double minx, double maxx, double minchi,
184  {
185  TGraph* ret = Slice(expt, calc, v, nbinsx, minx, maxx, minchi, opts);
186 
187  for(int i = 0; i < ret->GetN(); ++i){
188  double x, y;
189  ret->GetPoint(i, x, y);
190  ret->SetPoint(i, x, y > 0 ? sqrt(y) : 0);
191  }
192 
193  return ret;
194  }
195 
196  //----------------------------------------------------------------------
197  TGraph* FindValley(const IExperiment* expt,
199  const IFitVar& scanVar,
200  const IFitVar& fitVar,
201  int nbinsx, double xmin, double xmax,
202  const std::vector<const IFitVar*>& profVars,
203  const std::vector<const ISyst*>& profSysts,
204  const SeedList& seedPts,
205  const std::vector<SystShifts>& systSeedPts,
206  bool transpose,
208  {
209  TGraph* g = new TGraph;
210  g->SetLineWidth(2);
211 
212  std::vector<const IFitVar*> vars = profVars;
213  vars.push_back(&fitVar);
214 
215  for(int i = 0; i <= nbinsx; ++i){
216  const double x = xmin + i*(xmax-xmin)/nbinsx;
217  scanVar.SetValue(calc, x);
218  MinuitFitter fit(expt, vars, profSysts, opts);
219  SystShifts shiftSeed = SystShifts::Nominal();
220  fit.Fit(calc, shiftSeed, seedPts, systSeedPts, MinuitFitter::kQuiet);
221  const double y = fitVar.GetValue(calc);
222  if(transpose)
223  g->SetPoint(i, y, x);
224  else
225  g->SetPoint(i, x, y);
226  }
227 
228  return g;
229  }
230 
231  //----------------------------------------------------------------------
232  std::vector<double> FindCurveCrossings(TH1* h, double critVal)
233  {
234  std::vector<double> ret;
235 
236  // Don't look at the under/overflow bins because they don't have properly
237  // defined centre positions in any case. Stop one short because we compare
238  // i and i+1 inside the loop.
239  for(int i = 1; i < h->GetNbinsX(); ++i){
240  const double x0 = h->GetXaxis()->GetBinCenter(i);
241  const double x1 = h->GetXaxis()->GetBinCenter(i+1);
242 
243  const double y0 = h->GetBinContent(i);
244  const double y1 = h->GetBinContent(i+1);
245 
246  // Passed the critical value between the two points
247  if((y0 < critVal) != (y1 < critVal)){
248  // Interpolate the actual crossing point
249  const double x = x0 + (x1-x0)*(critVal-y0)/(y1-y0);
250  ret.push_back(x);
251  }
252  } // end for i
253 
254  return ret;
255  }
256 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
set< int >::iterator it
Float_t y1[n_points_granero]
Definition: compare.C:5
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:209
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
std::vector< double > FindCurveCrossings(TH1 *h, double critVal)
Intended for use on the output of Profile.
Definition: Fit.cxx:232
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
static std::map< const ISyst *, TGraph * > empty_syst_map
Definition: Fit.h:27
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const XML_Char * s
Definition: expat.h:262
expt
Definition: demo5.py:34
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:64
#define pot
T GetShift(const ISyst *syst) const
static std::map< const IFitVar *, TGraph * > empty_vars_map
Definition: Fit.h:26
const std::map< std::pair< std::string, std::string >, Variable > vars
TGraph * Slice(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
scan in one variable, holding all others constant
Definition: Fit.cxx:169
double POT() const
Definition: Spectrum.h:231
TGraph * SqrtSlice(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
Forward to Slice but sqrt the result for a crude significance.
Definition: Fit.cxx:180
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
Int_t nbinsx
Definition: plot.C:23
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
double SimpleFOM(const Spectrum &obs, const Spectrum &unosc, double pot)
Figure-of-merit with no systematics, for binned data.
Definition: Fit.cxx:23
Base class defining interface for experiments.
Definition: IExperiment.h:14
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
Interface definition for fittable variables.
Definition: IFitVar.h:16
static float_mat transpose(const float_mat &a)
returns the transposed matrix.
void Done()
Call this when action is completed.
Definition: Progress.cxx:91
TGraph * FindValley(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar &scanVar, const IFitVar &fitVar, int nbinsx, double xmin, double xmax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, bool transpose, MinuitFitter::FitOpts opts)
Find the minimum in one variable as a function of another.
Definition: Fit.cxx:197
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.
Definition: Fit.cxx:143
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17