Public Member Functions | Static Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
ana::FrequentistSurface Class Reference

Log-likelihood scan across two parameters. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-10-20/CAFAna/Fit/FrequentistSurface.h"

Inheritance diagram for ana::FrequentistSurface:
ana::ISurface ana::SurfaceKrige

Public Member Functions

 FrequentistSurface (const IExperiment *expt, osc::IOscCalcAdjustable *calc, const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars={}, const std::vector< const ISyst * > &profSysts={}, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, bool parallel=false, MinuitFitter::FitOpts opts=MinuitFitter::kNormal)
 
 FrequentistSurface (const IExperiment *expt, osc::IOscCalcAdjustable *calc, 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={}, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, bool parallel=false, MinuitFitter::FitOpts opts=MinuitFitter::kNormal)
 
virtual ~FrequentistSurface ()
 
std::vector< TH2 * > GetProfiledHists ()
 Maps of the values taken on by the profiled parameters. More...
 
void SaveTo (TDirectory *dir, const std::string &name) const
 
double BestLikelihood () const
 
double GetBestFitX () const
 
double GetBestFitY () const
 
void Draw () const
 Draw the surface itself. More...
 
void DrawBestFit (Color_t color, Int_t marker=kFullCircle) const
 Draw the best fit point. More...
 
void DrawContour (TH2 *fc, Style_t style, Color_t color, double minchi=-1)
 
std::vector< TGraph * > GetGraphs (TH2 *fc, double minchi=-1)
 For expert use, custom painting of contours. More...
 
TH2 * ToTH2 (double minchi=-1) const
 
void SetTitle (const char *str)
 

Static Public Member Functions

static std::unique_ptr< FrequentistSurfaceLoadFrom (TDirectory *dir, const std::string &name)
 

Protected Member Functions

 FrequentistSurface ()
 
void CreateHistograms (const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts)
 
std::string ProgressBarTitle (const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts) const
 
virtual void FillSurface (const IExperiment *expt, osc::IOscCalcAdjustable *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)
 
double FillSurfacePoint (const IExperiment *expt, osc::IOscCalcAdjustable *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)
 
void FindMinimum (const IExperiment *expt, osc::IOscCalcAdjustable *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)
 
double BinCenterX (int bin) const
 
double BinCenterY (int bin) const
 
void SaveToHelper (TDirectory *dir) const
 dir should already be the appropriate sub-directory More...
 
void EnsureAxes () const
 

Static Protected Member Functions

static void FillSurfObj (ISurface &surf, TDirectory *dir)
 

Protected Attributes

bool fParallel
 
MinuitFitter::FitOpts fFitOpts
 
std::vector< TH2 * > fProfHists
 
double fBestLikelihood
 
double fBestFitX
 
double fBestFitY
 
TH2F * fHist
 
bool fLogX
 
bool fLogY
 
std::vector< double > fSeedValues
 

Detailed Description

Log-likelihood scan across two parameters.

Definition at line 27 of file FrequentistSurface.h.

Constructor & Destructor Documentation

ana::FrequentistSurface::FrequentistSurface ( const IExperiment expt,
osc::IOscCalcAdjustable calc,
const FitAxis xax,
const FitAxis yax,
const std::vector< const IFitVar * > &  profVars = {},
const std::vector< const ISyst * > &  profSysts = {},
const SeedList seedPts = SeedList(),
const std::vector< SystShifts > &  systSeedPts = {},
bool  parallel = false,
MinuitFitter::FitOpts  opts = MinuitFitter::kNormal 
)
Parameters
exptThe experiment object to draw $ \chi^2 $ values from
calcValues for oscillation parameters to be held fixed
xaxDefinition of the x-axis
yaxDefinition of the y-axis
profVarsOscillation parameters to profile over
profSystsSystematic parameters to profile over
seedPtsTry all combinations of these params as seeds
systSeedPtsTry all of these systematic combinations as seeds
parallelUse all the cores on this machine? Be careful...

Definition at line 32 of file FrequentistSurface.cxx.

References CreateHistograms(), FillSurface(), FindMinimum(), ana::ISurface::fLogX, ana::ISurface::fLogY, ana::ISurface::fSeedValues, ana::FitAxis::islog, registry_explorer::v, and ana::FitAxis::var.

42  : fParallel(parallel), fFitOpts(opts)
43  {
44  fLogX = xax.islog;
45  fLogY = yax.islog;
46 
47  CreateHistograms(xax, yax, profVars, profSysts);
48 
49  for(const IFitVar* v: profVars) fSeedValues.push_back(v->GetValue( calc));
50 
51  FillSurface(expt, calc, xax.var, yax.var, profVars, profSysts, seedPts, systSeedPts);
52 
53  FindMinimum(expt, calc, xax.var, yax.var, profVars, profSysts, seedPts, systSeedPts);
54  }
MinuitFitter::FitOpts fFitOpts
expt
Definition: demo5.py:34
void FindMinimum(const IExperiment *expt, osc::IOscCalcAdjustable *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)
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
std::vector< double > fSeedValues
Definition: ISurface.h:68
void CreateHistograms(const FitAxis &xax, const FitAxis &yax, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts)
virtual void FillSurface(const IExperiment *expt, osc::IOscCalcAdjustable *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)
ana::FrequentistSurface::FrequentistSurface ( const IExperiment expt,
osc::IOscCalcAdjustable calc,
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 = {},
const SeedList seedPts = SeedList(),
const std::vector< SystShifts > &  systSeedPts = {},
bool  parallel = false,
MinuitFitter::FitOpts  opts = MinuitFitter::kNormal 
)
inline
Parameters
exptThe experiment object to draw $ \chi^2 $ values from
calcValues for oscillation parameters to be held fixed
xvarOscillation parameter to place on the x axis
nbinsxNumber of bins along x axis
xminMinimum value of x axis
xmaxMaximum value of x axis
nbinsyNumber of bins along y axis
yminMinimum value of y axis
ymaxMaximum value of y axis
profVarsOscillation parameters to profile over
profSystsSystematic parameters to profile over
seedPtsTry all combinations of these params as seeds
systSeedPtsTry all of these systematic combinations as seeds
parallelUse all the cores on this machine? Be careful...

Definition at line 63 of file FrequentistSurface.h.

References FrequentistSurface(), ana::MinuitFitter::kNormal, plot_validation_datamc::opts, and ~FrequentistSurface().

67  {},
68  const std::vector<const ISyst*>& profSysts = {},
69  const SeedList& seedPts = SeedList(),
70  const std::vector<SystShifts>& systSeedPts = {},
71  bool parallel = false,
73  : FrequentistSurface(expt, calc,
74  FitAxis(xvar, nbinsx, xmin, xmax, false),
75  FitAxis(yvar, nbinsy, ymin, ymax, false),
76  profVars, profSysts, seedPts, systSeedPts, parallel, opts)
77  {
78  }
std::map< std::string, double > xmax
Double_t ymax
Definition: plot.C:25
default if un-specified
Definition: MinuitFitter.h:25
expt
Definition: demo5.py:34
Int_t nbinsx
Definition: plot.C:23
Double_t ymin
Definition: plot.C:24
ana::FrequentistSurface::~FrequentistSurface ( )
virtual

Definition at line 57 of file FrequentistSurface.cxx.

References CreateHistograms().

Referenced by FrequentistSurface().

58  {
59  }
ana::FrequentistSurface::FrequentistSurface ( )
inlineprotected

Member Function Documentation

double ana::ISurface::BestLikelihood ( ) const
inlineinherited

Definition at line 29 of file ISurface.h.

References ana::ISurface::fBestLikelihood.

29 {return fBestLikelihood;}
double fBestLikelihood
Definition: ISurface.h:63
double ana::ISurface::BinCenterX ( int  bin) const
protectedinherited

Definition at line 244 of file ISurface.cxx.

References visualisationForPaperMasterPlot::ax, ana::ISurface::fHist, and ana::ISurface::fLogX.

Referenced by ana::BayesianSurface::BuildHist(), ana::ISurface::EnsureAxes(), FillSurface(), FindMinimum(), and ana::ISurface::GetBestFitY().

245  {
246  const TAxis* ax = fHist->GetXaxis();
247  return fLogX ? ax->GetBinCenterLog(bin) : ax->GetBinCenter(bin);
248  }
TH2F * fHist
Definition: ISurface.h:66
bool fLogX
Definition: ISurface.h:67
float bin[41]
Definition: plottest35.C:14
double ana::ISurface::BinCenterY ( int  bin) const
protectedinherited

Definition at line 251 of file ISurface.cxx.

References visualisationForPaperMasterPlot::ax, ana::ISurface::fHist, and ana::ISurface::fLogY.

Referenced by ana::BayesianSurface::BuildHist(), ana::ISurface::EnsureAxes(), FillSurface(), FindMinimum(), and ana::ISurface::GetBestFitY().

252  {
253  const TAxis* ax = fHist->GetYaxis();
254  return fLogY ? ax->GetBinCenterLog(bin) : ax->GetBinCenter(bin);
255  }
TH2F * fHist
Definition: ISurface.h:66
float bin[41]
Definition: plottest35.C:14
bool fLogY
Definition: ISurface.h:67
void ana::FrequentistSurface::CreateHistograms ( const FitAxis xax,
const FitAxis yax,
const std::vector< const IFitVar * > &  profVars,
const std::vector< const ISyst * > &  profSysts 
)
protected

Definition at line 63 of file FrequentistSurface.cxx.

References ana::ExpandedHistogram(), ana::ISurface::fHist, fProfHists, MECModelEnuComparisons::i, ana::FitAxis::islog, ana::IFitVar::LatexName(), ana::FitAxis::max, ana::FitAxis::min, ana::FitAxis::nbins, ProgressBarTitle(), plotROC::title, and ana::FitAxis::var.

Referenced by FrequentistSurface(), and ~FrequentistSurface().

66  {
67  fHist = ExpandedHistogram(";"+xax.var->LatexName()+";"+yax.var->LatexName(),
68  xax.nbins, xax.min, xax.max, xax.islog,
69  yax.nbins, yax.min, yax.max, yax.islog);
70 
71  for(unsigned int i = 0; i < profVars.size()+profSysts.size(); ++i){
73  if(i < profVars.size())
74  title = profVars[i]->LatexName();
75  else
76  title = profSysts[i-profVars.size()]->LatexName();
77 
78  fProfHists.push_back(ExpandedHistogram(title+";"+xax.var->LatexName()+";"+yax.var->LatexName(),
79  xax.nbins, xax.min, xax.max, xax.islog,
80  yax.nbins, yax.min, yax.max, yax.islog));
81  }
82 
83  }
TH2F * ExpandedHistogram(const std::string &title, int nbinsx, double xmin, double xmax, bool xlog, int nbinsy, double ymin, double ymax, bool ylog)
Internal helper for Surface and FCSurface.
Definition: Utilities.cxx:153
std::vector< TH2 * > fProfHists
TH2F * fHist
Definition: ISurface.h:66
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void ana::ISurface::Draw ( ) const
inherited

Draw the surface itself.

Definition at line 19 of file ISurface.cxx.

References ana::ISurface::EnsureAxes(), and ana::ISurface::fHist.

Referenced by cc(), demo5(), ana::ISurface::GetBestFitY(), modularextrap_demo_nue(), modularextrap_demo_numu(), Plot2DSlices(), plot_3flavor_withsysts(), test_ana(), test_stanfit_statsonly(), and test_stanfit_withsysts().

20  {
21  EnsureAxes();
22 
23  fHist->Draw("colz same");
24 
25  // colz obliterated them
26  gPad->RedrawAxis();
27 
28  gPad->Update();
29  }
TH2F * fHist
Definition: ISurface.h:66
void EnsureAxes() const
Definition: ISurface.cxx:63
void ana::ISurface::DrawBestFit ( Color_t  color,
Int_t  marker = kFullCircle 
) const
inherited
void ana::ISurface::DrawContour ( TH2 *  fc,
Style_t  style,
Color_t  color,
double  minchi = -1 
)
inherited
Parameters
fcSurface to compare against for this significance level
styleLine style for TAttLine, solid, dotted, dashed etc
colorLine color for TAttLine
minchi$\chi^2$ of best fit to compare against. Default: best fit from this surface.

Definition at line 44 of file ISurface.cxx.

References ana::ISurface::EnsureAxes(), MECModelEnuComparisons::g, ana::ISurface::GetGraphs(), and APDGainPoints::gs.

Referenced by cc(), demo5(), Draw2DSurface(), drawSensitivity(), ana::ISurface::GetBestFitY(), getSensitivity(), make_contours(), make_nom_expt(), modularextrap_demo_nue(), modularextrap_demo_numu(), Plot2DSlice(), Plot2DSlices(), plot_3flavor_withsysts(), plotContProf(), plots(), starPlot(), syst_test(), template_basic(), template_GENIE_systs(), template_nonGENIE_systs(), test_ana(), test_stanfit_statsonly(), and test_stanfit_withsysts().

46  {
47  EnsureAxes();
48 
49  std::vector<TGraph *> gs = GetGraphs(fc, minchi);
50 
51  for (TGraph *g: gs)
52  {
53  g->SetLineWidth(3);//2);
54  g->SetLineStyle(style);
55  g->SetLineColor(color);
56  g->Draw("l");
57  }
58 
59  gPad->Update();
60  }
std::vector< TGraph * > GetGraphs(TH2 *fc, double minchi=-1)
For expert use, custom painting of contours.
Definition: ISurface.cxx:106
void EnsureAxes() const
Definition: ISurface.cxx:63
void ana::ISurface::EnsureAxes ( ) const
protectedinherited

Definition at line 63 of file ISurface.cxx.

References visualisationForPaperMasterPlot::ax, file_size_ana::axes, ana::ISurface::BinCenterX(), ana::ISurface::BinCenterY(), ana::ISurface::fHist, ana::ISurface::fLogX, ana::ISurface::fLogY, genie::utils::style::Format(), and ana::UniqueName().

Referenced by ana::ISurface::Draw(), ana::ISurface::DrawBestFit(), and ana::ISurface::DrawContour().

64  {
65  // Could have a file temporarily open
66  DontAddDirectory guard;
67 
68  // If this pad has already been drawn in, already has axes
69  if (gPad && !gPad->GetListOfPrimitives()->IsEmpty()) return;
70 
71  const TAxis *ax = fHist->GetXaxis();
72  const TAxis *ay = fHist->GetYaxis();
73  const double Nx = ax->GetNbins();
74  const double Ny = ay->GetNbins();
75 
76  // Axes with limits where the user originally requested them, which we
77  // adjusted to be the centres of the first and last bins.
78  TH2 *axes = new TH2C(UniqueName().c_str(),
79  TString::Format(";%s;%s",
80  ax->GetTitle(), ay->GetTitle()),
81  Nx - 1, BinCenterX(1), BinCenterX(Nx),
82  Ny - 1, BinCenterY(1), BinCenterY(Ny));
83  axes->Draw();
84 
85  if(fHist){
86  // "colz same" will reuse axis's min and max, so set them helpfully here
87  axes->SetMinimum(fHist->GetMinimum());
88  axes->SetMaximum(fHist->GetMaximum());
89  }
90 
91  axes->SetTitle(fHist->GetTitle());
92  axes->GetXaxis()->SetLabelSize(ax->GetLabelSize());
93  axes->GetYaxis()->SetLabelSize(ay->GetLabelSize());
94  axes->GetXaxis()->SetLabelOffset(ax->GetLabelOffset());
95  axes->GetYaxis()->SetLabelOffset(ay->GetLabelOffset());
96  axes->GetXaxis()->CenterTitle();
97  axes->GetYaxis()->CenterTitle();
98 
99  if(fLogX) gPad->SetLogx();
100  if(fLogY) gPad->SetLogy();
101 
102  gPad->Update();
103  }
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
TH2F * fHist
Definition: ISurface.h:66
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void ana::FrequentistSurface::FillSurface ( const IExperiment expt,
osc::IOscCalcAdjustable 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 
)
protectedvirtual

Reimplemented in ana::SurfaceKrige.

Definition at line 108 of file FrequentistSurface.cxx.

References ana::ThreadPool::AddTask(), bin, ana::ISurface::BinCenterX(), ana::ISurface::BinCenterY(), om::cerr, ana::IExperiment::ChiSq(), ana::Progress::Done(), e, allTimeWatchdog::endl, ana::ISurface::fHist, FillSurfacePoint(), ana::ThreadPool::Finish(), genie::utils::style::Format(), fParallel, ana::ThreadPool::NThreads(), ana::IFitVar::Penalty(), time_estimates::pool, cacheDefinitionData::prog, ProgressBarTitle(), ana::Progress::SetProgress(), ana::IFitVar::ShortName(), ana::ThreadPool::ShowProgress(), fillBadChanDBTables::step, submit_syst::x, and submit_syst::y.

Referenced by FrequentistSurface().

115  {
116  // Nothing created during surface filling belongs in a
117  // directory. Unfortunately the local guards in Spectrum etc are racey when
118  // run in parallel. But this should cover the whole lot safely.
119  DontAddDirectory guard;
120 
121  const std::string progTitle = ProgressBarTitle(xvar, yvar, profVars, profSysts);
122 
123  Progress *prog = 0;
124  // Difficult to keep a progress bar properly up to date when threaded
125  if (!fParallel) prog = new Progress(progTitle);
126  ThreadPool *pool = 0;
127 
128  if(fParallel){
129  // A hack/workaround needed for parallel running:
130  //
131  // Give all the constituents of the Prediction a chance to do lazy
132  // initialization, before they race themselves trying to do it in
133  // parallel.
134  expt->ChiSq(calc);
135 
136  pool = new ThreadPool;
137  pool->ShowProgress(progTitle + TString::Format(" using %d threads", pool->NThreads()).Data());
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 = BinCenterX(x);
160  const double yv = BinCenterY(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  ThreadPool::func_t task = [=](){
178  FillSurfacePoint(expt, calc,
179  xvar, xv, yvar, yv,
180  profVars, profSysts, seedPts, systSeedPts);
181  };
182 
183  if(fParallel){
184  pool->AddTask(task);
185  }
186  else{
187  task(); // Just do it straight away
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  pool->Finish();
198  delete pool;
199  }
200  else{
201  prog->Done();
202  delete prog;
203  }
204  }
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
double FillSurfacePoint(const IExperiment *expt, osc::IOscCalcAdjustable *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)
std::function< void(void)> func_t
The type of the user&#39;s worker functions.
Definition: ThreadPool.h:33
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
OStream cerr
Definition: OStream.cxx:7
std::string ProgressBarTitle(const IFitVar *xvar, const IFitVar *yvar, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts) const
expt
Definition: demo5.py:34
TH2F * fHist
Definition: ISurface.h:66
float bin[41]
Definition: plottest35.C:14
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Float_t e
Definition: plot.C:35
double ana::FrequentistSurface::FillSurfacePoint ( const IExperiment expt,
osc::IOscCalcAdjustable 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 
)
protected

Definition at line 207 of file FrequentistSurface.cxx.

References calc, ana::IExperiment::ChiSq(), osc::_IOscCalcAdjustable< T >::Copy(), fFitOpts, ana::ISurface::fHist, ana::IFitter::Fit(), fParallel, fProfHists, ana::ISurface::fSeedValues, ana::SystShifts::GetShift(), MECModelEnuComparisons::i, makeTrainCVSamples::int, calib::j, ana::IFitter::kQuiet, ana::MinuitFitter::SetFitOpts(), and ana::IFitVar::SetValue().

Referenced by FillSurface(), and FrequentistSurface().

215  {
216  if(fParallel){
217  // Need to take our own copy so that we don't get overwritten by someone
218  // else's changes.
219  calc = calc->Copy();
220  }
221 
222  xvar->SetValue(calc, x);
223  yvar->SetValue(calc, y);
224 
225  //Make sure that the profiled values of fitvars do not persist between steps.
226  for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
227 
228  double chi;
229  if(profVars.empty() && profSysts.empty()){
230  chi = expt->ChiSq(calc);
231  }
232  else{
233  MinuitFitter fitter(expt, profVars, profSysts);
234  fitter.SetFitOpts(fFitOpts);
235  SystShifts bestSysts;
236  chi = fitter.Fit(calc, bestSysts, seedPts, systSeedPts, MinuitFitter::kQuiet)->EvalMetricVal();
237 
238  for(unsigned int i = 0; i < profVars.size(); ++i){
239  fProfHists[i]->Fill(x, y, profVars[i]->GetValue(calc));
240  }
241  for(unsigned int j = 0; j < profSysts.size(); ++j){
242  fProfHists[j+profVars.size()]->Fill(x, y, bestSysts.GetShift(profSysts[j]));
243  }
244  }
245 
246  fHist->Fill(x, y, chi);
247 
248  if(fParallel) delete calc;
249 
250  return chi;
251  }
virtual _IOscCalcAdjustable< T > * Copy() const =0
MinuitFitter::FitOpts fFitOpts
osc::OscCalcDumb calc
std::vector< TH2 * > fProfHists
expt
Definition: demo5.py:34
TH2F * fHist
Definition: ISurface.h:66
const double j
Definition: BetheBloch.cxx:29
std::vector< double > fSeedValues
Definition: ISurface.h:68
void ana::ISurface::FillSurfObj ( ISurface surf,
TDirectory *  dir 
)
staticprotectedinherited

Definition at line 214 of file ISurface.cxx.

References ana::ISurface::fBestFitX, ana::ISurface::fBestFitY, ana::ISurface::fBestLikelihood, ana::ISurface::fHist, ana::ISurface::fSeedValues, compare_h5_caf::idx, and registry_explorer::v.

Referenced by ana::ISurface::GetBestFitY(), ana::BayesianSurface::LoadFrom(), and LoadFrom().

215  {
216 
217  const TVectorD v = *(TVectorD *) dir->Get("minValues");
218  const TVectorD s = *(TVectorD *) dir->Get("seeds");
219 
220  surf.fHist = (TH2F *) dir->Get("hist");
221  surf.fBestLikelihood = v[0];
222  surf.fBestFitX = v[1];
223  surf.fBestFitY = v[2];
224 
225  for (int idx = 0; idx < s.GetNrows(); ++idx)
226  surf.fSeedValues.push_back(s[idx]);
227 
228  }
const XML_Char * s
Definition: expat.h:262
TDirectory * dir
Definition: macro.C:5
surf
Definition: demo5.py:35
void ana::FrequentistSurface::FindMinimum ( const IExperiment expt,
osc::IOscCalcAdjustable 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 
)
protected

Definition at line 255 of file FrequentistSurface.cxx.

References ana::ISurface::BinCenterX(), ana::ISurface::BinCenterY(), ana::ISurface::fBestFitX, ana::ISurface::fBestFitY, ana::ISurface::fBestLikelihood, fFitOpts, ana::ISurface::fHist, ana::IFitter::Fit(), ana::ISurface::fSeedValues, ana::IFitVar::GetValue(), MECModelEnuComparisons::i, makeTrainCVSamples::int, ana::SystShifts::Nominal(), ana::MinuitFitter::SetFitOpts(), ana::IFitVar::SetValue(), submit_syst::x, and submit_syst::y.

Referenced by FrequentistSurface().

262  {
263  // Location of the best minimum found from filled surface
264  double minchi = 1e10;
265  int minx = fHist->GetNbinsX()/2;
266  int miny = fHist->GetNbinsY()/2;
267  for(int x = 1; x <= fHist->GetNbinsX(); ++x){
268  for(int y = 1; y <= fHist->GetNbinsY(); ++y){
269  const double chi = fHist->GetBinContent(x, y);
270  if(chi < minchi){
271  minchi = chi;
272  minx = x;
273  miny = y;
274  }
275  }
276  }
277 
278  std::vector<const IFitVar*> allVars = {xvar, yvar};
279  allVars.insert(allVars.end(), profVars.begin(), profVars.end());
280  MinuitFitter fit(expt, allVars, profSysts);
281  fit.SetFitOpts(fFitOpts);
282  // Seed from best grid point
283  xvar->SetValue(calc, BinCenterX(minx));
284  yvar->SetValue(calc, BinCenterY(miny));
285  for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
286  SystShifts systSeed = SystShifts::Nominal();
287  fBestLikelihood = fit.Fit(calc, systSeed, seedPts)->EvalMetricVal();
288  fBestFitX = xvar->GetValue(calc);
289  fBestFitY = yvar->GetValue(calc);
290 
291  for(int x = 0; x < fHist->GetNbinsX()+2; ++x){
292  for(int y = 0; y < fHist->GetNbinsY()+2; ++y){
293  fHist->SetBinContent(x, y, fHist->GetBinContent(x, y)-fBestLikelihood);
294  }
295  }
296 
297  fHist->SetMinimum(0);
298  }
MinuitFitter::FitOpts fFitOpts
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
static SystShifts Nominal()
Definition: SystShifts.h:34
expt
Definition: demo5.py:34
TH2F * fHist
Definition: ISurface.h:66
double fBestFitY
Definition: ISurface.h:65
std::vector< double > fSeedValues
Definition: ISurface.h:68
double fBestFitX
Definition: ISurface.h:64
double fBestLikelihood
Definition: ISurface.h:63
double ana::ISurface::GetBestFitX ( ) const
inlineinherited

Definition at line 30 of file ISurface.h.

References ana::ISurface::fBestFitX.

Referenced by test_stanfit_statsonly().

30 {return fBestFitX;}
double fBestFitX
Definition: ISurface.h:64
double ana::ISurface::GetBestFitY ( ) const
inlineinherited
std::vector< TGraph * > ana::ISurface::GetGraphs ( TH2 *  fc,
double  minchi = -1 
)
inherited

For expert use, custom painting of contours.

Definition at line 106 of file ISurface.cxx.

References ana::ISurface::fBestLikelihood, ana::ISurface::fHist, runNovaSAM::ret, PandAna.Demos.tute_pid_validation::specs, demo5::surf, tmp, and ana::UniqueName().

Referenced by ana::ISurface::DrawContour(), and ana::ISurface::GetBestFitY().

107  {
108  std::vector<TGraph*> ret;
109 
110  if(minchi < 0) minchi = fBestLikelihood;
111  std::unique_ptr<TH2> surf((TH2*)fHist->Clone(UniqueName().c_str()));
112  surf->Add(fc, -1);
113 
114  TVirtualPad* bak = gPad;
115 
116  const bool wasbatch = gROOT->IsBatch();
117  gROOT->SetBatch(); // User doesn't want to see our temporary canvas
118  TCanvas tmp;
119 
120  gStyle->SetOptStat(0);
121 
122  const double level = minchi-fBestLikelihood;
123  surf->SetContour(1, &level);
124  surf->Draw("cont list");
125 
126  tmp.Update();
127  tmp.Paint();
128 
129  gROOT->SetBatch(wasbatch);
130  gPad = bak;
131 
132  // The graphs we need (contained inside TLists, contained inside
133  // TObjArrays) are in the list of specials. But we need to be careful about
134  // types, because other stuff can get in here too (TDatabasePDG for
135  // example).
136  TCollection* specs = gROOT->GetListOfSpecials();
137 
138  TIter nextSpec(specs);
139  while(TObject* spec = nextSpec()){
140  if(!spec->InheritsFrom(TObjArray::Class())) continue;
141  TObjArray* conts = (TObjArray*)spec;
142 
143  if(conts->IsEmpty()) continue;
144 
145  if(!conts->At(0)->InheritsFrom(TList::Class())) continue;
146  TList* cont = (TList*)conts->At(0);
147 
148  TIter nextObj(cont);
149  // Contour could be split into multiple pieces
150  while(TObject* obj = nextObj()){
151  if(!obj->InheritsFrom(TGraph::Class())) continue;
152 
153  ret.push_back((TGraph*)obj->Clone(UniqueName().c_str()));
154  } // end for obj
155  } // end for spec
156 
157  return ret;
158  }
Float_t tmp
Definition: plot.C:36
TH2F * fHist
Definition: ISurface.h:66
double fBestLikelihood
Definition: ISurface.h:63
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
surf
Definition: demo5.py:35
std::vector<TH2*> ana::FrequentistSurface::GetProfiledHists ( )
inline

Maps of the values taken on by the profiled parameters.

Definition at line 83 of file FrequentistSurface.h.

References dir, fProfHists, LoadFrom(), and SaveTo().

Referenced by run_joint_fit_2020_contours().

83 {return fProfHists;}
std::vector< TH2 * > fProfHists
std::unique_ptr< FrequentistSurface > ana::FrequentistSurface::LoadFrom ( TDirectory *  dir,
const std::string &  name 
)
static

Definition at line 327 of file FrequentistSurface.cxx.

References ana::assert(), dir, ana::ISurface::FillSurfObj(), genie::utils::style::Format(), make_syst_table_plots::h, compare_h5_caf::idx, logx, logy, demo5::surf, and getGoodRuns4SAM::tag.

Referenced by compare_fits(), demoFitContours(), DrawExtrapSurface(), GetProfiledHists(), joint_fit_2018_contours(), joint_fit_2019_contours(), median_contours(), plot_joint_fit_2020_contours(), Plotting_CompareMultipleContours(), sensitivity2018(), and sensitivity2020().

328  {
329  dir = dir->GetDirectory(name.c_str()); // switch to subdir
330  assert(dir);
331 
332  DontAddDirectory guard;
333 
334  TObjString *tag = (TObjString *) dir->Get("type");
335  assert(tag);
336  assert(tag->GetString() == "FrequentistSurface" || tag->GetString() == "Surface");
337  delete tag;
338 
339  std::unique_ptr<FrequentistSurface> surf(new FrequentistSurface);
341 
342  for (std::size_t idx = 0; idx < surf->fSeedValues.size(); ++idx)
343  {
344  // Search for old "marg" name here too for backwards compatibility
345  TH2 *h = (TH2 *) dir->Get(TString::Format("profHists/hist%lu", idx));
346  if (h)
347  surf->fProfHists.push_back(h);
348  else
349  surf->fProfHists.push_back((TH2 *) dir->Get(TString::Format("margHists/hist%lu", idx)));
350  }
351 
352  const TObjString* logx = (TObjString*)dir->Get("logx");
353  const TObjString* logy = (TObjString*)dir->Get("logy");
354  // Tolerate missing log tags for backwards compatibility
355  surf->fLogX = (logx && logx->GetString() == "yes");
356  surf->fLogY = (logy && logy->GetString() == "yes");
357 
358  delete dir;
359 
360  return surf;
361  }
const XML_Char * name
Definition: expat.h:151
bool logy
bool logx
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
static void FillSurfObj(ISurface &surf, TDirectory *dir)
Definition: ISurface.cxx:214
surf
Definition: demo5.py:35
std::string ana::FrequentistSurface::ProgressBarTitle ( const IFitVar xvar,
const IFitVar yvar,
const std::vector< const IFitVar * > &  profVars,
const std::vector< const ISyst * > &  profSysts 
) const
protected

Definition at line 87 of file FrequentistSurface.cxx.

References ana::IFitVar::ShortName(), plotROC::title, and registry_explorer::v.

Referenced by CreateHistograms(), FillSurface(), and FrequentistSurface().

90  {
91  std::string title = "Filling surface for "+yvar->ShortName()+" vs "+xvar->ShortName();
92 
93  if(!profVars.empty() || !profSysts.empty()){
94  title += " (profiling ";
95 
96  for(const IFitVar* v: profVars) title += v->ShortName() + ", ";
97  for(const ISyst* s: profSysts) title += s->ShortName() + ", ";
98 
99  // Always have one superfluous ", " at the end
100  title.resize(title.size()-2);
101  title += ")";
102  }
103 
104  return title;
105  }
const XML_Char * s
Definition: expat.h:262
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void ana::FrequentistSurface::SaveTo ( TDirectory *  dir,
const std::string &  name 
) const

Definition at line 301 of file FrequentistSurface.cxx.

References dir, genie::utils::style::Format(), fProfHists, compare_h5_caf::idx, it, ana::ISurface::SaveToHelper(), and tmp.

Referenced by getContProf(), getContProf_Sensitivity(), GetProfiledHists(), joint_fit_future_contour_univ(), make_surfprof(), make_surfprof_sensitivity(), MakeExtrapSurface(), run_joint_fit_2020_contours(), sensitivity2018(), and sensitivity2020().

302  {
303  TDirectory *tmp = gDirectory;
304 
305  dir = dir->mkdir(name.c_str()); // switch to subdir
306  dir->cd();
307 
308  TObjString("FrequentistSurface").Write("type");
309 
311 
312  TDirectory *profDir = dir->mkdir("profHists");
313  int idx = 0;
315  {
316  profDir->cd();
317  it->Write(TString::Format("hist%d", idx++));
318  }
319 
320  dir->Write();
321  delete dir;
322 
323  tmp->cd();
324  }
const XML_Char * name
Definition: expat.h:151
set< int >::iterator it
Float_t tmp
Definition: plot.C:36
std::vector< TH2 * > fProfHists
void SaveToHelper(TDirectory *dir) const
dir should already be the appropriate sub-directory
Definition: ISurface.cxx:189
TDirectory * dir
Definition: macro.C:5
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void ana::ISurface::SaveToHelper ( TDirectory *  dir) const
protectedinherited

dir should already be the appropriate sub-directory

Definition at line 189 of file ISurface.cxx.

References ana::ISurface::fBestFitX, ana::ISurface::fBestFitY, ana::ISurface::fBestLikelihood, ana::ISurface::fHist, ana::ISurface::fLogX, ana::ISurface::fLogY, ana::ISurface::fSeedValues, and registry_explorer::v.

Referenced by ana::ISurface::GetBestFitY(), ana::BayesianSurface::SaveTo(), and SaveTo().

190  {
191  TDirectory *oldDir = gDirectory;
192 
193  dir->cd();
194 
195  TVectorD v(3);
196  v[0] = fBestLikelihood;
197  v[1] = fBestFitX;
198  v[2] = fBestFitY;
199  v.Write("minValues");
200 
201  fHist->Write("hist");
202 
203  TVectorD s(fSeedValues.size(), &fSeedValues[0]);
204  s.Write("seeds");
205 
206  TObjString(fLogX ? "yes" : "no").Write("logx");
207  TObjString(fLogY ? "yes" : "no").Write("logy");
208 
209  if (oldDir)
210  oldDir->cd();
211  }
const XML_Char * s
Definition: expat.h:262
TH2F * fHist
Definition: ISurface.h:66
bool fLogX
Definition: ISurface.h:67
double fBestFitY
Definition: ISurface.h:65
bool fLogY
Definition: ISurface.h:67
std::vector< double > fSeedValues
Definition: ISurface.h:68
TDirectory * dir
Definition: macro.C:5
double fBestFitX
Definition: ISurface.h:64
double fBestLikelihood
Definition: ISurface.h:63
void ana::ISurface::SetTitle ( const char *  str)
inherited

Definition at line 183 of file ISurface.cxx.

References ana::ISurface::fHist.

Referenced by ana::ISurface::GetBestFitY(), and plots().

184  {
185  fHist->SetTitle(str);
186  }
TH2F * fHist
Definition: ISurface.h:66
TH2 * ana::ISurface::ToTH2 ( double  minchi = -1) const
inherited

Definition at line 161 of file ISurface.cxx.

References ana::ISurface::fBestLikelihood, ana::ISurface::fHist, runNovaSAM::ret, submit_syst::x, and submit_syst::y.

Referenced by CAF_makeCAFSensitivities_for_FNEX(), demoFitContours(), ana::Flat(), ana::ISurface::GetBestFitY(), joint_fit_2017_contours(), joint_fit_2018_contours(), joint_fit_2019_contours(), joint_fit_future_contour_univ(), MakeCAFSensitivities_for_FNEX(), Plot2DSlice(), Plot2DSlices(), plot_3flavor_withsysts(), run_joint_fit_2020_contours(), sensitivity2018(), sensitivity2020(), and test_stanfit_withsysts().

162  {
163  // Could have a file temporarily open
164  DontAddDirectory guard;
165 
166  TH2 *ret = new TH2F(*fHist);
167 
168  if (minchi >= 0)
169  {
170  for (int x = 0; x < ret->GetNbinsX() + 2; ++x)
171  {
172  for (int y = 0; y < ret->GetNbinsY() + 2; ++y)
173  {
174  ret->SetBinContent(x, y, ret->GetBinContent(x, y) + fBestLikelihood - minchi);
175  }
176  }
177  }
178 
179  return ret;
180  }
TH2F * fHist
Definition: ISurface.h:66
double fBestLikelihood
Definition: ISurface.h:63

Member Data Documentation

double ana::ISurface::fBestFitX
protectedinherited
double ana::ISurface::fBestFitY
protectedinherited
double ana::ISurface::fBestLikelihood
protectedinherited
MinuitFitter::FitOpts ana::FrequentistSurface::fFitOpts
protected

Definition at line 126 of file FrequentistSurface.h.

Referenced by FillSurfacePoint(), and FindMinimum().

TH2F* ana::ISurface::fHist
protectedinherited
bool ana::ISurface::fLogX
protectedinherited
bool ana::ISurface::fLogY
protectedinherited
bool ana::FrequentistSurface::fParallel
protected

Definition at line 124 of file FrequentistSurface.h.

Referenced by FillSurface(), and FillSurfacePoint().

std::vector<TH2*> ana::FrequentistSurface::fProfHists
protected

Definition at line 129 of file FrequentistSurface.h.

Referenced by CreateHistograms(), FillSurfacePoint(), GetProfiledHists(), and SaveTo().

std::vector<double> ana::ISurface::fSeedValues
protectedinherited

The documentation for this class was generated from the following files: