FrequentistSurface.cxx
Go to the documentation of this file.
1 #include <iostream>
2 #include <functional>
3 
4 #include "TCanvas.h"
5 #include "TGraph.h"
6 #include "TH2F.h"
7 #include "Minuit2/StackAllocator.h"
8 #include "TObjArray.h"
9 #include "TPad.h"
10 #include "TROOT.h"
11 #include "TStyle.h"
12 #include "TKey.h"
13 #include "TVectorD.h"
14 #include "TObjString.h"
15 #include "TCollection.h"
16 
21 #include "CAFAna/Core/IFitVar.h"
22 #include "CAFAna/Core/Progress.h"
23 #include "CAFAna/Core/ThreadPool.h"
24 #include "CAFAna/Core/Utilities.h"
25 
27 
28 namespace ana
29 {
30  //----------------------------------------------------------------------
33  const IFitVar* xvar, int nbinsx, double xmin, double xmax,
34  const IFitVar* yvar, int nbinsy, double ymin, double ymax,
35  const std::vector<const IFitVar*>& profVars,
36  const std::vector<const ISyst*>& profSysts,
37  const SeedList& seedPts,
38  const std::vector<SystShifts>& systSeedPts,
39  bool parallel,
41  : fParallel(parallel), fFitOpts(opts)
42  {
43  CreateHistograms(xvar, nbinsx, xmin, xmax,
44  yvar, nbinsy, ymin, ymax,
45  profVars, profSysts);
46 
47  for(const IFitVar* v: profVars) fSeedValues.push_back(v->GetValue( calc));
48 
49  FillSurface(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);
50 
51  FindMinimum(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);
52  }
53 
54  //---------------------------------------------------------------------
56  {
57  }
58 
59  //---------------------------------------------------------------------
61  CreateHistograms(const IFitVar* xvar, int nbinsx, double xmin, double xmax,
62  const IFitVar* yvar, int nbinsy, double ymin, double ymax,
63  const std::vector<const IFitVar*>& profVars,
64  const std::vector<const ISyst*>& profSysts)
65  {
66  fHist = ExpandedHistogram(";"+xvar->LatexName()+";"+yvar->LatexName(),
67  nbinsx, xmin, xmax,
68  nbinsy, ymin, ymax);
69 
70  for(unsigned int i = 0; i < profVars.size()+profSysts.size(); ++i){
72  if(i < profVars.size())
73  title = profVars[i]->LatexName();
74  else
75  title = profSysts[i-profVars.size()]->LatexName();
76 
77  fProfHists.push_back(ExpandedHistogram(title+";"+xvar->LatexName()+";"+yvar->LatexName(),
78  nbinsx, xmin, xmax,
79  nbinsy, ymin, ymax));
80  }
81 
82  }
83 
84  //---------------------------------------------------------------------
86  ProgressBarTitle(const IFitVar* xvar, const IFitVar* yvar,
87  const std::vector<const IFitVar*>& profVars,
88  const std::vector<const ISyst*>& profSysts) const
89  {
90  std::string title = "Filling surface for "+yvar->ShortName()+" vs "+xvar->ShortName();
91 
92  if(!profVars.empty() || !profSysts.empty()){
93  title += " (profiling ";
94 
95  for(const IFitVar* v: profVars) title += v->ShortName() + ", ";
96  for(const ISyst* s: profSysts) title += s->ShortName() + ", ";
97 
98  // Always have one superfluous ", " at the end
99  title.resize(title.size()-2);
100  title += ")";
101  }
102 
103  return title;
104  }
105 
106  //---------------------------------------------------------------------
109  const IFitVar *xvar, const IFitVar *yvar,
110  const std::vector<const IFitVar *> &profVars,
111  const std::vector<const ISyst *> &profSysts,
112  const SeedList& seedPts,
113  const std::vector<SystShifts> &systSeedPts)
114  {
115  if (fParallel && !(profVars.empty() && profSysts.empty()))
116  {
117  // Minuit calls this at some point, and it creates a static. Which
118  // doesn't like happening on multiple threads at once. So do it upfront
119  // while we're still single-threaded.
121  }
122 
123  // Nothing created during surface filling belongs in a
124  // directory. Unfortunately the local guards in Spectrum etc are racey when
125  // run in parallel. But this should cover the whole lot safely.
126  DontAddDirectory guard;
127 
128  const std::string progTitle = ProgressBarTitle(xvar, yvar, profVars, profSysts);
129 
130  Progress *prog = 0;
131  // Difficult to keep a progress bar properly up to date when threaded
132  if (!fParallel) prog = new Progress(progTitle);
133  ThreadPool *pool = 0;
134  if (fParallel)
135  {
136  pool = new ThreadPool;
137  pool->ShowProgress(progTitle);
138  }
139 
140  const int Nx = fHist->GetNbinsX();
141  const int Ny = fHist->GetNbinsY();
142 
143  // Fill bins in "random" order so that the progress bar is accurate
144  // straight away instead of possibly being misled by whatever atypical
145  // points we start with. This step is a prime which guarantees we get every
146  // cell.
147  int step = 7919;
148  // Very unlikely (Nx or Ny is a multiple of step), but just to be safe.
149  if ((Nx * Ny) % step == 0) step = 1;
150 
151  int bin = 0;
152  int neval = 0;
153 
154  do
155  {
156  const int x = bin % Nx + 1;
157  const int y = bin / Nx + 1;
158 
159  const double xv = fHist->GetXaxis()->GetBinCenter(x);
160  const double yv = fHist->GetYaxis()->GetBinCenter(y);
161 
162  if (xvar->Penalty(xv, calc) > 1e-10)
163  {
164  std::cerr << "Warning! " << xvar->ShortName() << " = " << xv
165  << " has penalty of " << xvar->Penalty(xv, calc)
166  << " that could have been applied in surface. "
167  << "This should never happen." << std::endl;
168  }
169  if (yvar->Penalty(yv, calc) > 1e-10)
170  {
171  std::cerr << "Warning! " << yvar->ShortName() << " = " << yv
172  << " has penalty of " << yvar->Penalty(yv, calc)
173  << " that could have been applied in surface. "
174  << "This should never happen." << std::endl;
175  }
176 
177  if (fParallel)
178  {
180  expt, calc,
181  xvar, xv, yvar, yv,
182  profVars, profSysts, seedPts, systSeedPts);
183  } else
184  {
185  FillSurfacePoint(expt, calc,
186  xvar, xv, yvar, yv,
187  profVars, profSysts, seedPts, systSeedPts);
188  ++neval;
189  prog->SetProgress(neval / double(Nx * Ny));
190  }
191 
192  bin = (bin + step) % (Nx * Ny);
193  } while (bin != 0);
194 
195 
196  if (fParallel)
197  {
198  pool->Finish();
199  delete pool;
200  } else
201  {
202  prog->Done();
203  delete prog;
204  }
205  }
206 
207  /// \brief Helper for FrequentistSurface::FillSurfacePoint
208  ///
209  /// The cacheing of the nominal done in PredictionInterp is not
210  /// threadsafe. This is an inelegant but pragmatic way of suppressing it.
212  {
213  public:
215 
217  {
218  std::cout << "FrequentistSurface::OscCalcNoHash not copyable." << std::endl;
219  abort();
220  }
221 
223 
224  double P(int a, int b, double E) override {return fCalc->P(a, b, E);}
225  /// Marks this calculator "unhashable" so the cacheing won't occur
226  TMD5* GetParamsHash() const override {return 0;}
227 
228  // It's a shame we have to forward all the getters and setters explicitly,
229  // but I don't see how to avoid it. Take away some drudgery with a macro.
230 #define F(v)\
231  void Set##v(double x) override {fCalc->Set##v(x);}\
232  double Get##v() const override {return fCalc->Get##v();}
233  F(L); F(Rho);
234 #undef F
235 #define F(v)\
236  void Set##v(const double& x) override {fCalc->Set##v(x);}\
237  double Get##v() const override {return fCalc->Get##v();}
238  F(Dmsq21); F(Dmsq32); F(Th12); F(Th13); F(Th23); F(dCP);
239 #undef F
240  protected:
242  };
243 
244  //----------------------------------------------------------------------
247  const IFitVar* xvar, double x,
248  const IFitVar* yvar, double y,
249  const std::vector<const IFitVar*>& profVars,
250  const std::vector<const ISyst*>& profSysts,
251  const SeedList& seedPts,
252  const std::vector<SystShifts>& systSeedPts)
253  {
254  osc::IOscCalculatorAdjustable* calcNoHash = 0; // specific to parallel mode
255 
256  if(fParallel){
257  // Need to take our own copy so that we don't get overwritten by someone
258  // else's changes.
259  calc = calc->Copy();
260  }
261 
262  xvar->SetValue(calc, x);
263  yvar->SetValue(calc, y);
264 
265  if (fParallel){
266  calcNoHash = new OscCalcNoHash(calc);
267  }
268 
269  //Make sure that the profiled values of fitvars do not persist between steps.
270  for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
271 
272  double chi;
273  if(profVars.empty() && profSysts.empty()){
274  chi = expt->ChiSq(fParallel ? calcNoHash : calc);
275  }
276  else{
277  MinuitFitter fitter(expt, profVars, profSysts);
278  fitter.SetFitOpts(fFitOpts);
279  SystShifts bestSysts;
280  chi = fitter.Fit(calc, bestSysts, seedPts, systSeedPts, MinuitFitter::kQuiet);
281 
282  for(unsigned int i = 0; i < profVars.size(); ++i){
283  fProfHists[i]->Fill(x, y, profVars[i]->GetValue(calc));
284  }
285  for(unsigned int j = 0; j < profSysts.size(); ++j){
286  fProfHists[j+profVars.size()]->Fill(x, y, bestSysts.GetShift(profSysts[j]));
287  }
288  }
289 
290  fHist->Fill(x, y, chi);
291 
292  if(fParallel){
293  delete calc;
294  delete calcNoHash;
295  }
296  return chi;
297  }
298 
299 
300  //---------------------------------------------------------------------
303  const IFitVar* xvar, const IFitVar* yvar,
304  const std::vector<const IFitVar*>& profVars,
305  const std::vector<const ISyst*>& profSysts,
306  const SeedList& seedPts,
307  const std::vector<SystShifts>& systSeedPts)
308  {
309  // Location of the best minimum found from filled surface
310  double minchi = 1e10;
311  int minx = fHist->GetNbinsX()/2;
312  int miny = fHist->GetNbinsY()/2;
313  for(int x = 1; x <= fHist->GetNbinsX(); ++x){
314  for(int y = 1; y <= fHist->GetNbinsY(); ++y){
315  const double chi = fHist->GetBinContent(x, y);
316  if(chi < minchi){
317  minchi = chi;
318  minx = x;
319  miny = y;
320  }
321  }
322  }
323 
324  std::vector<const IFitVar*> allVars = {xvar, yvar};
325  allVars.insert(allVars.end(), profVars.begin(), profVars.end());
326  MinuitFitter fit(expt, allVars, profSysts);
327  fit.SetFitOpts(fFitOpts);
328  // Seed from best grid point
329  xvar->SetValue(calc, fHist->GetXaxis()->GetBinCenter(minx));
330  yvar->SetValue(calc, fHist->GetYaxis()->GetBinCenter(miny));
331  for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
332  SystShifts systSeed = SystShifts::Nominal();
333  fBestLikelihood = fit.Fit(calc, systSeed, seedPts);
334  fBestFitX = xvar->GetValue(calc);
335  fBestFitY = yvar->GetValue(calc);
336 
337  for(int x = 0; x < fHist->GetNbinsX()+2; ++x){
338  for(int y = 0; y < fHist->GetNbinsY()+2; ++y){
339  fHist->SetBinContent(x, y, fHist->GetBinContent(x, y)-fBestLikelihood);
340  }
341  }
342 
343  fHist->SetMinimum(0);
344  }
345 
346  //----------------------------------------------------------------------
347  void FrequentistSurface::SaveTo(TDirectory *dir) const
348  {
349  TDirectory *tmp = gDirectory;
350  dir->cd();
351  TObjString("FrequentistSurface").Write("type");
352 
353  ISurface::SaveTo(dir);
354 
355  TDirectory *profDir = dir->mkdir("profHists");
356  int idx = 0;
358  {
359  profDir->cd();
360  it->Write(TString::Format("hist%d", idx++));
361  }
362 
363  tmp->cd();
364  }
365 
366 //----------------------------------------------------------------------
367  std::unique_ptr<FrequentistSurface> FrequentistSurface::LoadFrom(TDirectory *dir)
368  {
369  DontAddDirectory guard;
370 
371  TObjString *tag = (TObjString *) dir->Get("type");
372  assert(tag);
373  assert(tag->GetString() == "FrequentistSurface" || tag->GetString() == "Surface");
374 
375  std::unique_ptr<FrequentistSurface> surf(new FrequentistSurface);
376  ISurface::FillSurfObj(*surf, dir);
377 
378  for (std::size_t idx = 0; idx < surf->fSeedValues.size(); ++idx)
379  {
380  // Search for old "marg" name here too for backwards compatibility
381  TH2 *h = (TH2 *) dir->Get(TString::Format("profHists/hist%lu", idx));
382  if (h)
383  surf->fProfHists.push_back(h);
384  else
385  surf->fProfHists.push_back((TH2 *) dir->Get(TString::Format("margHists/hist%lu", idx)));
386  }
387 
388  return surf;
389  }
390 
391 
392  // See eg the statistics section of the PDG
393  TH2* Gaussian68Percent2D(const FrequentistSurface& s){return Flat(2.30, s);}
394  TH2* Gaussian90Percent2D(const FrequentistSurface& s){return Flat(4.61, s);}
395  TH2* Gaussian95Percent2D(const FrequentistSurface& s){return Flat(5.99, s);}
396  TH2* Gaussian2Sigma2D (const FrequentistSurface& s){return Flat(6.18, s);}
397  TH2* Gaussian99Percent2D(const FrequentistSurface& s){return Flat(9.21, s);}
398  TH2* Gaussian3Sigma2D (const FrequentistSurface& s){return Flat(11.83, s);}
399 
400  TH2* Gaussian68Percent1D(const FrequentistSurface& s){return Flat(1.00, s);}
401  TH2* Gaussian90Percent1D(const FrequentistSurface& s){return Flat(2.71, s);}
402  TH2* Gaussian95Percent1D(const FrequentistSurface& s){return Flat(3.84, s);}
403  TH2* Gaussian2Sigma1D (const FrequentistSurface& s){return Flat(4.00, s);}
404  TH2* Gaussian99Percent1D(const FrequentistSurface& s){return Flat(6.63, s);}
405  TH2* Gaussian3Sigma1D (const FrequentistSurface& s){return Flat(9.00, s);}
406 
407 
408 } // namespace
virtual double Penalty(double, osc::IOscCalculatorAdjustable *) const
Definition: IFitVar.h:34
void CreateHistograms(const IFitVar *xvar, int nbinsx, double xmin, double xmax, const IFitVar *yvar, int nbinsy, double ymin, double ymax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts)
void Finish()
Wait for all threads to complete before returning.
Definition: ThreadPool.cxx:44
const std::string & LatexName() const
Definition: IFitVar.h:37
virtual void FillSurface(const IExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
TH2 * Gaussian90Percent1D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 1D in gaussian approximation.
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
std::map< std::string, double > xmax
set< int >::iterator it
MinuitFitter::FitOpts fFitOpts
Helper for FrequentistSurface::FillSurfacePoint.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
OscCalcNoHash(osc::IOscCalculatorAdjustable *c)
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
void SaveTo(TDirectory *dir) const
double FillSurfacePoint(const IExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *xvar, double x, const IFitVar *yvar, double y, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
OStream cerr
Definition: OStream.cxx:7
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
const float Rho() const
Definition: HoughPeak.h:87
Float_t tmp
Definition: plot.C:36
TH2 * Gaussian2Sigma1D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 1D in gaussian approximation.
static SystShifts Nominal()
Definition: SystShifts.h:34
TMD5 * GetParamsHash() const override
Marks this calculator "unhashable" so the cacheing won&#39;t occur.
std::string ProgressBarTitle(const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts) const
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::vector< TH2 * > fProfHists
Double_t ymax
Definition: plot.C:25
osc::IOscCalculatorAdjustable * Copy() const override
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
TH2 * Gaussian95Percent2D(const FrequentistSurface &s)
Up-value surface for 95% confidence in 2D in gaussian approximation.
expt
Definition: demo5.py:34
TH2F * fHist
Definition: ISurface.h:62
static constexpr double L
Float_t E
Definition: plot.C:20
double dCP
void SaveTo(TDirectory *dir) const
Definition: ISurface.cxx:198
const double a
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
T GetShift(const ISyst *syst) const
void AddMemberTask(T *obj, M meth, A...args)
Add member function task, with arguments.
Definition: ThreadPool.h:78
TH2 * Gaussian99Percent2D(const FrequentistSurface &s)
Up-value surface for 99% confidence in 2D in gaussian approximation.
double fBestFitY
Definition: ISurface.h:61
TH2 * Gaussian3Sigma1D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 1D in gaussian approximation.
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir)
const double j
Definition: BetheBloch.cxx:29
float bin[41]
Definition: plottest35.C:14
virtual void SetValue(osc::IOscCalculatorAdjustable *osc, double val) const =0
osc::OscCalculatorDumb calc
OStream cout
Definition: OStream.cxx:6
std::vector< double > fSeedValues
Definition: ISurface.h:63
virtual double P(int flavBefore, int flavAfter, double E)=0
E in GeV; flavors as PDG codes (so, neg==>antinu)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
const std::string & ShortName() const
Definition: IFitVar.h:36
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
double fBestFitX
Definition: ISurface.h:60
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
virtual double ChiSq(osc::IOscCalculatorAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
Definition: IExperiment.h:18
virtual _IOscCalculatorAdjustable< T > * Copy() const =0
Base class defining interface for experiments.
Definition: IExperiment.h:14
const hit & b
Definition: hits.cxx:21
void FindMinimum(const IExperiment *expt, osc::IOscCalculatorAdjustable *calc, const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts)
A very simple thread pool for use by Surface.
Definition: ThreadPool.h:17
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
void ShowProgress(const std::string &title)
Definition: ThreadPool.cxx:38
Interface definition for fittable variables.
Definition: IFitVar.h:16
Double_t ymin
Definition: plot.C:24
virtual double GetValue(const osc::IOscCalculatorAdjustable *osc) const =0
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
TH2 * Gaussian99Percent1D(const FrequentistSurface &s)
Up-value surface for 99% confidence in 1D in gaussian approximation.
static void FillSurfObj(ISurface &surf, TDirectory *dir)
Definition: ISurface.cxx:219
osc::IOscCalculatorAdjustable * fCalc
#define F(v)
TH2 * Flat(double level, const ISurface &s)
Helper function for the gaussian approximation surfaces.
Definition: ISurface.cxx:237
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
Float_t e
Definition: plot.C:35
virtual double Fit(osc::IOscCalculatorAdjustable *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:56
TH2 * Gaussian95Percent1D(const FrequentistSurface &s)
Up-value surface for 95% confidence in 1D in gaussian approximation.
void Done()
Call this when action is completed.
Definition: Progress.cxx:91
double fBestLikelihood
Definition: ISurface.h:59
TH2 * Gaussian68Percent1D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 1D in gaussian approximation.
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax)
Definition: Utilities.cxx:198
double P(int a, int b, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
surf
Definition: demo5.py:35
void SetFitOpts(FitOpts opts)
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:15