Public Types | Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Static Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
ana::BayesianSurface Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-25/CAFAna/Fit/BayesianSurface.h"

Inheritance diagram for ana::BayesianSurface:
ana::ISurface ana::BayesianMarginal

Public Types

enum  MarginalMode { MarginalMode::kHistogram, MarginalMode::kLLWgtdHistogram, MarginalMode::kKNN }
 How to construct marginal distributions. More...
 

Public Member Functions

 BayesianSurface (const MCMCSamples &samples, const IFitVar *xvar, int nbinsx, double xmin, double xmax, const IFitVar *yvar, int nbinsy, double ymin, double ymax, MarginalMode mode=MarginalMode::kHistogram)
 Build a surface from the output of an MCMC fit. More...
 
 BayesianSurface (const MCMCSamples &samples, const FitAxis &xax, const FitAxis &yax, MarginalMode mode=MarginalMode::kHistogram)
 
 BayesianSurface (const MCMCSamples &samples, const IFitVar *xsyst, int nbinsx, double xmin, double xmax, const ISyst *ysyst, int nbinsy, double ymin, double ymax, MarginalMode mode=MarginalMode::kHistogram)
 Same as other constructors, but with one IFitVar and one ISyst. More...
 
 BayesianSurface (const MCMCSamples &samples, const ISyst *xsyst, int nbinsx, double xmin, double xmax, const IFitVar *ysyst, int nbinsy, double ymin, double ymax, MarginalMode mode=MarginalMode::kHistogram)
 Same as other constructors, but with one ISyst and one IFitVar. More...
 
 BayesianSurface (const MCMCSamples &samples, const ISyst *xsyst, int nbinsx, double xmin, double xmax, const ISyst *ysyst, int nbinsy, double ymin, double ymax, MarginalMode mode=MarginalMode::kHistogram)
 Same as other constructors, but with two ISysts. More...
 
TH2 * QuantileSurface (Quantile quantile) const
 
TH2 * QuantileSurface (double quantile, const MCMCSamples &mcmsamples) const
 
void SaveTo (TDirectory *dir) 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)
 
double QuantileThreshold (Quantile quantile, const TH1 *pdf=nullptr) const
 
double QuantileThreshold (double quantile, const MCMCSamples &samples, const TH1 *pdf=nullptr) const
 
std::vector< const ISyst * > Systs () const
 
std::vector< const IFitVar * > Vars () const
 

Static Public Member Functions

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

Static Public Attributes

static const std::vector< std::pair< Quantile, double > > QuantileUpVals
 

Protected Member Functions

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
 
double EstimateLLFromKNN (const std::vector< float > &vals) const
 
std::unique_ptr< TH1 > ToHistogram (const std::vector< Binning > &bins) const
 

Static Protected Member Functions

static void FillSurfObj (ISurface &surf, TDirectory *dir)
 
static double ThresholdFromTH1 (const TH1 *pdf, double quantile)
 Determine the threshold for the 'histogram' mode from the given histogram with a brute force method. More...
 
static void LoadInto (BayesianMarginal *marg, TDirectory *dir)
 

Protected Attributes

double fBestLikelihood
 
double fBestFitX
 
double fBestFitY
 
TH2F * fHist
 
bool fLogX
 
bool fLogY
 
std::vector< double > fSeedValues
 
MarginalMode fMode
 What mode are we using? More...
 
const MCMCSamplesfMCMCSamples
 Ref to the MCMCSamples. Will not be persisted. More...
 
std::vector< std::stringfOrderedBrNames
 the variables/systs to be marginalized over, in axis order (x-axis, y-axis, ...) More...
 

Private Member Functions

 BayesianSurface ()=default
 
template<typename SystOrVar1 , typename SystOrVar2 >
void BuildHist (const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
 Build the TH2 from the kNN. More...
 

Detailed Description

Definition at line 15 of file BayesianSurface.h.

Member Enumeration Documentation

How to construct marginal distributions.

Enumerator
kHistogram 
kLLWgtdHistogram 
kKNN 

Definition at line 45 of file BayesianMarginal.h.

46  {
47  kHistogram,
48  kLLWgtdHistogram,
49  kKNN,
50  };

Constructor & Destructor Documentation

ana::BayesianSurface::BayesianSurface ( const MCMCSamples samples,
const IFitVar xvar,
int  nbinsx,
double  xmin,
double  xmax,
const IFitVar yvar,
int  nbinsy,
double  ymin,
double  ymax,
MarginalMode  mode = MarginalMode::kHistogram 
)

Build a surface from the output of an MCMC fit.

Parameters
samplesMCMCSamples object from the fit.
xvarOscillation parameter to place on the x axis
nbinsxNumber of bins along x axis
xminMinimum value of x axis
xmaxMaximum value of x axis
yvarOscillation parameter to place on the y axis
nbinsyNumber of bins along y axis
yminMinimum value of y axis
ymaxMaximum value of y axis

Definition at line 25 of file BayesianSurface.cxx.

References BuildHist(), ana::ISurface::fLogX, ana::ISurface::fLogY, submit_nova_art::mode, nbinsx, xmax, make_mec_shifts_plots::xmin, ymax, and ymin.

29  : BayesianMarginal(samples, {xvar, yvar}, mode)
30  {
31  fLogX = fLogY = false;
32  BuildHist(samples, xvar, nbinsx, xmin, xmax, yvar, nbinsy, ymin, ymax);
33  }
std::map< std::string, double > xmax
Double_t ymax
Definition: plot.C:25
void BuildHist(const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
Build the TH2 from the kNN.
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
Int_t nbinsx
Definition: plot.C:23
Double_t ymin
Definition: plot.C:24
ana::BayesianSurface::BayesianSurface ( const MCMCSamples samples,
const FitAxis xax,
const FitAxis yax,
MarginalMode  mode = MarginalMode::kHistogram 
)

Definition at line 13 of file BayesianSurface.cxx.

References BuildHist(), ana::ISurface::fLogX, ana::ISurface::fLogY, submit_nova_art::mode, and ana::FitAxis::var.

17  : BayesianMarginal(samples, {xax.var, yax.var}, mode)
18  {
19  fLogX = xax.islog;
20  fLogY = yax.islog;
21  BuildHist(samples, xax.var, xax.nbins, xax.min, xax.max, yax.var, yax.nbins, yax.min, yax.max);
22  }
void BuildHist(const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
Build the TH2 from the kNN.
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
ana::BayesianSurface::BayesianSurface ( const MCMCSamples samples,
const IFitVar xsyst,
int  nbinsx,
double  xmin,
double  xmax,
const ISyst ysyst,
int  nbinsy,
double  ymin,
double  ymax,
MarginalMode  mode = MarginalMode::kHistogram 
)

Same as other constructors, but with one IFitVar and one ISyst.

Definition at line 47 of file BayesianSurface.cxx.

References BuildHist(), ana::ISurface::fLogX, ana::ISurface::fLogY, submit_nova_art::mode, nbinsx, xmax, make_mec_shifts_plots::xmin, ymax, and ymin.

51  : BayesianMarginal(samples, {xvar, ysyst}, mode)
52  {
53  fLogX = fLogY = false;
54  BuildHist(samples, xvar, nbinsx, xmin, xmax, ysyst, nbinsy, ymin, ymax);
55  }
std::map< std::string, double > xmax
Double_t ymax
Definition: plot.C:25
void BuildHist(const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
Build the TH2 from the kNN.
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
Int_t nbinsx
Definition: plot.C:23
Double_t ymin
Definition: plot.C:24
ana::BayesianSurface::BayesianSurface ( const MCMCSamples samples,
const ISyst xsyst,
int  nbinsx,
double  xmin,
double  xmax,
const IFitVar ysyst,
int  nbinsy,
double  ymin,
double  ymax,
MarginalMode  mode = MarginalMode::kHistogram 
)

Same as other constructors, but with one ISyst and one IFitVar.

Definition at line 36 of file BayesianSurface.cxx.

References BuildHist(), ana::ISurface::fLogX, ana::ISurface::fLogY, submit_nova_art::mode, nbinsx, xmax, make_mec_shifts_plots::xmin, ymax, and ymin.

40  : BayesianMarginal(samples, {xsyst, yvar}, mode)
41  {
42  fLogX = fLogY = false;
43  BuildHist(samples, xsyst, nbinsx, xmin, xmax, yvar, nbinsy, ymin, ymax);
44  }
std::map< std::string, double > xmax
Double_t ymax
Definition: plot.C:25
void BuildHist(const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
Build the TH2 from the kNN.
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
Int_t nbinsx
Definition: plot.C:23
Double_t ymin
Definition: plot.C:24
ana::BayesianSurface::BayesianSurface ( const MCMCSamples samples,
const ISyst xsyst,
int  nbinsx,
double  xmin,
double  xmax,
const ISyst ysyst,
int  nbinsy,
double  ymin,
double  ymax,
MarginalMode  mode = MarginalMode::kHistogram 
)

Same as other constructors, but with two ISysts.

Definition at line 58 of file BayesianSurface.cxx.

References BuildHist(), ana::ISurface::fLogX, ana::ISurface::fLogY, submit_nova_art::mode, nbinsx, xmax, make_mec_shifts_plots::xmin, ymax, and ymin.

62  : BayesianMarginal(samples, {xsyst, ysyst}, mode)
63  {
64  fLogX = fLogY = false;
65  BuildHist(samples, xsyst, nbinsx, xmin, xmax, ysyst, nbinsy, ymin, ymax);
66  }
std::map< std::string, double > xmax
Double_t ymax
Definition: plot.C:25
void BuildHist(const MCMCSamples &samples, const SystOrVar1 *x, int nbinsx, double xmin, double xmax, const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
Build the TH2 from the kNN.
bool fLogX
Definition: ISurface.h:67
bool fLogY
Definition: ISurface.h:67
Int_t nbinsx
Definition: plot.C:23
Double_t ymin
Definition: plot.C:24
ana::BayesianSurface::BayesianSurface ( )
privatedefault

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 BuildHist(), ana::ISurface::EnsureAxes(), ana::FrequentistSurface::FillSurface(), ana::FrequentistSurface::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 BuildHist(), ana::ISurface::EnsureAxes(), ana::FrequentistSurface::FillSurface(), ana::FrequentistSurface::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
template<typename SystOrVar1 , typename SystOrVar2 >
template void ana::BayesianSurface::BuildHist ( const MCMCSamples samples,
const SystOrVar1 *  x,
int  nbinsx,
double  xmin,
double  xmax,
const SystOrVar2 *  y,
int  nbinsy,
double  ymin,
double  ymax 
)
private

Build the TH2 from the kNN.

Definition at line 70 of file BayesianSurface.cxx.

References ana::MCMCSamples::BestFitSampleIdx(), ana::ISurface::BinCenterX(), ana::ISurface::BinCenterY(), ana::BayesianMarginal::EstimateLLFromKNN(), ana::ExpandedHistogram(), ana::ISurface::fBestFitX, ana::ISurface::fBestFitY, ana::ISurface::fBestLikelihood, ana::ISurface::fHist, ana::ISurface::fLogX, ana::ISurface::fLogY, ana::BayesianMarginal::fMode, ana::BayesianMarginal::kHistogram, ana::BayesianMarginal::kKNN, ana::BayesianMarginal::kLLWgtdHistogram, nbinsx, ana::MCMCSamples::SampleLL(), ana::MCMCSamples::SampleValue(), ana::Binning::Simple(), ana::BayesianMarginal::ToHistogram(), febshutoff_auto::val, submit_syst::x, xmax, make_mec_shifts_plots::xmin, submit_syst::y, ymax, and ymin.

Referenced by BayesianSurface().

73  {
76  "BayesianSurface::BuildHist() can only be used with IFitVar or ISyst objects");
77 
78  auto bestSampleIdx = samples.BestFitSampleIdx();
79  fBestLikelihood = samples.SampleLL(bestSampleIdx);
80  fBestFitX = samples.SampleValue(x, bestSampleIdx);
81  fBestFitY = samples.SampleValue(y, bestSampleIdx);
82 
83  fHist = ExpandedHistogram(";"+x->LatexName()+";"+y->LatexName(),
84  nbinsx, xmin, xmax, fLogX,
85  nbinsy, ymin, ymax, fLogY);
86 
87  std::unique_ptr<TH1> tmpHist;
89  tmpHist = ToHistogram({Binning::Simple(nbinsx, xmin, xmax),
90  Binning::Simple(nbinsy, ymin, ymax)});
91  for (int xBin = 0; xBin < fHist->GetNbinsX() + 2; xBin++)
92  {
93  for (int yBin = 0; yBin < fHist->GetNbinsY() + 2; yBin++)
94  {
95  double val = std::numeric_limits<double>::signaling_NaN();
97  val = tmpHist->GetBinContent(xBin, yBin);
98  else if (fMode == MarginalMode::kKNN)
99  {
100  float xBinCtr = BinCenterX(xBin);
101  float yBinCtr = BinCenterY(yBin);
102  val = EstimateLLFromKNN({xBinCtr, yBinCtr});
103  }
104  fHist->SetBinContent(xBin, yBin, val);
105  }
106  }
107  }
std::map< std::string, double > xmax
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
double EstimateLLFromKNN(const std::vector< float > &vals) const
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:149
Double_t ymax
Definition: plot.C:25
const XML_Char int const XML_Char * value
Definition: expat.h:331
TH2F * fHist
Definition: ISurface.h:66
bool fLogX
Definition: ISurface.h:67
double fBestFitY
Definition: ISurface.h:65
bool fLogY
Definition: ISurface.h:67
MarginalMode fMode
What mode are we using?
Int_t nbinsx
Definition: plot.C:23
double fBestFitX
Definition: ISurface.h:64
Double_t ymin
Definition: plot.C:24
std::unique_ptr< TH1 > ToHistogram(const std::vector< Binning > &bins) const
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
double fBestLikelihood
Definition: ISurface.h:63
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
double ana::BayesianMarginal::EstimateLLFromKNN ( const std::vector< float > &  vals) const
protectedinherited

Definition at line 136 of file BayesianMarginal.cxx.

References ana::assert(), ana::BayesianMarginal::fkNN, ana::BayesianMarginal::fMode, and ana::BayesianMarginal::kKNN.

Referenced by BuildHist(), ana::Bayesian1DMarginal::QuantileRanges(), and ana::Bayesian1DMarginal::ToTH1().

137  {
139 
140  assert(vals.size() == fkNN->DataInfo().GetNVariables());
141 
142  TMVA::Event ev(vals, fkNN->DataInfo().GetNVariables());
143  return dynamic_cast<TMVA::MethodBase*>(fkNN.get())->GetRegressionValues(&ev)[0];
144  }
std::unique_ptr< TMVA::MethodKNN > fkNN
the kNN used for approximating the surface from the points when in kNN mode
MarginalMode fMode
What mode are we using?
assert(nhit_max >=nhit_nbins)
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(), LoadFrom(), and ana::FrequentistSurface::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
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::unique_ptr< BayesianSurface > ana::BayesianSurface::LoadFrom ( TDirectory *  dir,
const std::string name 
)
static

Definition at line 123 of file BayesianSurface.cxx.

References ana::assert(), dir, ana::ISurface::FillSurfObj(), ana::BayesianMarginal::LoadInto(), logx, logy, demo5::surf, and getGoodRuns4SAM::tag.

124  {
125  dir = dir->GetDirectory(name.c_str()); // switch to subdir
126  assert(dir);
127 
128  DontAddDirectory guard;
129 
130  auto tag = (TObjString *) dir->Get("type");
131  assert(tag);
132  assert(tag->GetString() == "BayesianSurface");
133 
134  // can't use make_unique<> because the one-arg constructor is protected
135  // and the std::forward() inside make_unique can't forward the visibility
136  std::unique_ptr<BayesianSurface> surf(new BayesianSurface);
138 
140 
141  const TObjString* logx = (TObjString*)dir->Get("logx");
142  const TObjString* logy = (TObjString*)dir->Get("logy");
143  // Tolerate missing log tags for backwards compatibility
144  surf->fLogX = (logx && logx->GetString() == "yes");
145  surf->fLogY = (logy && logy->GetString() == "yes");
146 
147  delete dir;
148 
149  return std::move(surf);
150  }
const XML_Char * name
Definition: expat.h:151
static void LoadInto(BayesianMarginal *marg, TDirectory *dir)
bool logy
bool logx
BayesianSurface()=default
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
static void FillSurfObj(ISurface &surf, TDirectory *dir)
Definition: ISurface.cxx:214
surf
Definition: demo5.py:35
void ana::BayesianMarginal::LoadInto ( BayesianMarginal marg,
TDirectory *  dir 
)
staticprotectedinherited

Note that the signature here is different from usual. This version expects to be passed a derived class instance created by a derived LoadFrom() as its argument...

Definition at line 147 of file BayesianMarginal.cxx.

References ana::BayesianMarginal::fdataLoader, ana::BayesianMarginal::fkNN, ana::BayesianMarginal::fMode, ana::BayesianMarginal::fQuantileSampleMap, ana::BayesianMarginal::fSysts, ana::BayesianMarginal::fVars, ana::BayesianMarginal::kKNN, submit_nova_art::mode, ana::BayesianMarginal::SetupKNN(), ana::BayesianMarginal::sKNNName, systNames, ana::BayesianMarginal::Systs(), compareCafs::tree, and ana::BayesianMarginal::Vars().

Referenced by LoadFrom().

148  {
149  auto varNames = dynamic_cast<TList*>(dir->Get("vars"));
150  for (const auto & obj : *varNames)
151  marg->fVars.push_back(Registry<IFitVar>::ShortNameToPtr(obj->GetName()));
152  auto systNames = dynamic_cast<TList*>(dir->Get("systs"));
153  for (const auto & obj : *systNames)
154  marg->fSysts.push_back(Registry<ISyst>::ShortNameToPtr(obj->GetName()));
155 
156  auto mode = dynamic_cast<TParameter<int>*>(dir->Get("mode"));
157  marg->fMode = MarginalMode(mode->GetVal());
158 
159  if (marg->fMode == MarginalMode::kKNN)
160  {
161  // again we're forced to make a temporary file (see SaveTo()). c'est la vie...
162  marg->fdataLoader = std::make_unique<TMVA::DataLoader>(sKNNName.c_str());
163  marg->SetupKNN(marg->fdataLoader.get());
164  auto tree = dynamic_cast<TTree *>(dir->Get("kNN_weights"));
165  tree->SetDirectory(nullptr);
166  TFile tmpFile("tmp_surface.root", "recreate");
167  tmpFile.cd();
168  tree->Write("knn");
169  marg->fkNN->ReadWeightsFromStream(tmpFile);
170  tmpFile.Close();
171  gSystem->Unlink("tmp_surface.root");
172  }
173 
174  auto quantMap = dynamic_cast<TMap*>(dir->Get("quantileSamples"));
175  for (const auto & quantPair : *quantMap)
176  {
177  auto pair = dynamic_cast<TPair*>(quantPair);
178  auto quant = Quantile(dynamic_cast<TParameter<int>*>(pair->Key())->GetVal());
179  marg->fQuantileSampleMap.emplace(std::piecewise_construct,
180  std::forward_as_tuple(quant),
181  std::forward_as_tuple(*dynamic_cast<TMap*>(pair->Value()),
182  marg->Vars(), marg->Systs()));
183  }
184 
185  }
MarginalMode
How to construct marginal distributions.
Quantile
Type to specify quantile values, so that we don&#39;t run into floating-point matching issues...
static const std::string sKNNName
TDirectory * dir
Definition: macro.C:5
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
static const T * ShortNameToPtr(const std::string &s, bool allowFail=false)
Definition: Registry.cxx:60
TH2 * ana::BayesianSurface::QuantileSurface ( Quantile  quantile) const

Up-value surface for a given quantile (e.g. 90%) Idea: Traditional 'highest posterior density' (HPD) estimate using the PDF of the LL we made in FillSurface. Use the "Density Quantile Approach" (R. Hyndman, Am. Statist. 50, 120): i.e., determine the 'critical' LL for significance alpha by taking the N sampled points, ordering them by their LL, and choosing the (1-alpha)*(N points)th point's LL as the cutoff. (Technically this is an estimator for the true critical LL, but it is normally distributed around the true value – see Hyndman – and can be made arbitrarily close by increasing the number of samples. Supposedly it has lots of nice properties, for instance, exact Frequentist coverage at level 1-alpha; see, e.g., http://astrostatistics.psu.edu/su11scma5/BayesComp11-CollectedSlides.pdf, pp. 270ff.) Since the surface's fHist already contains estimates of the PDF in each bin, we can just add a flat surface with the correct LL level (like is done in the Frequentist version) to get the right contours at LL = 0 in the sum.

Parameters
quantileThe quantile you want the surface made at (see Quantile)
Returns
A TH2 filled with the correct LL that can be added to fHist by GetContour() inherited from ISurface

Note that this is a method instead of a free function because it needs access to the stored LL values for each quantile. There is another version you can use if you need an abitrary quantile (though it requires you to supply the MCMCSamples to determine the LL).

Definition at line 153 of file BayesianSurface.cxx.

References ana::ISurface::fHist, ana::Flat(), and ana::BayesianMarginal::QuantileThreshold().

Referenced by plot_3flavor_withsysts(), test_stanfit_statsonly(), and test_stanfit_withsysts().

154  {
155  return Flat(QuantileThreshold(quantile, fHist), *this);
156  }
TH2F * fHist
Definition: ISurface.h:66
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
double QuantileThreshold(Quantile quantile, const TH1 *pdf=nullptr) const
TH2 * Flat(double level, const ISurface &s)
Helper function for the gaussian approximation surfaces.
Definition: ISurface.cxx:232
TH2 * ana::BayesianSurface::QuantileSurface ( double  quantile,
const MCMCSamples mcmsamples 
) const

Create quantile surface for an arbitrary quantile. Prefer the version above when the value is one of the known quantiles as those have been stored and don't require recomputation.

Parameters
quantileThe fractional quantile the surface should be made at (i.e. 90% –> 0.9).
mcmcsamplesThe MCMCSamples object originally used to make this surface (must match exactly!)

Definition at line 160 of file BayesianSurface.cxx.

References ana::ISurface::fHist, ana::Flat(), and ana::BayesianMarginal::QuantileThreshold().

161  {
162  return Flat(QuantileThreshold(quantile, mcmcsamples, fHist), *this);
163  }
TH2F * fHist
Definition: ISurface.h:66
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
double QuantileThreshold(Quantile quantile, const TH1 *pdf=nullptr) const
TH2 * Flat(double level, const ISurface &s)
Helper function for the gaussian approximation surfaces.
Definition: ISurface.cxx:232
double ana::BayesianMarginal::QuantileThreshold ( ana::Quantile  quantile,
const TH1 *  pdf = nullptr 
) const
inherited

Retrieve the value above which the remaining points comprise the given quantile.

Parameters
quantileThe quantile value desired
pdfIf fMode==MarginalMode::kHistogram, the returned value depends on the histogram that results from the marginalization (see on return val, below). Thus the relevant histogram must be supplied for that mode. For the other modes it may be left as nullptr.
Returns
Interpretation of return value depends on the MarginalMode fMode:

Definition at line 188 of file BayesianMarginal.cxx.

References ana::assert(), srt_file_template::find_if, ana::BayesianMarginal::fMode, ana::BayesianMarginal::fQuantileSampleMap, ana::BayesianMarginal::kHistogram, quantile(), ana::BayesianMarginal::QuantileUpVals, moon_position_table_new3::second, and ana::BayesianMarginal::ThresholdFromTH1().

Referenced by ana::Bayesian1DMarginal::QuantileRanges(), and QuantileSurface().

189  {
190  // in "histogram" mode, the output isn't weighted by LL;
191  // what's returned is simply the best-guess pdf
192  // approximated by the histogram itself.
193  // so instead of simply yielding the LL itself,
194  // we need to find the pdf value corresponding to the right breakpoint.
196  {
197  assert(pdf && "In 'histogram' marginal mode, must supply the histogram for which the quantile threshold is to be calculated");
199  QuantileUpVals.end(),
200  [quantile](const auto & pair) { return pair.first == quantile;})->second);
201  }
202 
203  return fQuantileSampleMap.at(quantile).LL();
204  }
std::map< Quantile, MCMCSample > fQuantileSampleMap
The samples that are the boundaries for various commonly-used quantiles.
static const std::vector< std::pair< Quantile, double > > QuantileUpVals
static double ThresholdFromTH1(const TH1 *pdf, double quantile)
Determine the threshold for the &#39;histogram&#39; mode from the given histogram with a brute force method...
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
MarginalMode fMode
What mode are we using?
assert(nhit_max >=nhit_nbins)
double ana::BayesianMarginal::QuantileThreshold ( double  quantile,
const MCMCSamples samples,
const TH1 *  pdf = nullptr 
) const
inherited

See other variant of QuantileThreshold() for discussion. This version computes the threshold for an arbitrary quantile, rather than the prestored ones, and so requires the original MCMCSamples.

Definition at line 207 of file BayesianMarginal.cxx.

References ana::assert(), ana::BayesianMarginal::fMode, ana::BayesianMarginal::kHistogram, ana::MCMCSamples::QuantileLL(), and ana::BayesianMarginal::ThresholdFromTH1().

208  {
209  if (quantile == -1)
210  quantile = 0;
211  assert(quantile >= 0 && quantile <= 1 && "QuantileSurface(): quantile must be between 0 and 1");
212 
213  auto quantilePair = samples.QuantileLL(quantile);
214 
216  return ThresholdFromTH1(pdf, quantile);
217  return quantilePair.second;
218  }
static double ThresholdFromTH1(const TH1 *pdf, double quantile)
Determine the threshold for the &#39;histogram&#39; mode from the given histogram with a brute force method...
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
MarginalMode fMode
What mode are we using?
assert(nhit_max >=nhit_nbins)
void ana::BayesianSurface::SaveTo ( TDirectory *  dir) const

Definition at line 166 of file BayesianSurface.cxx.

References ana::BayesianMarginal::SaveTo(), and ana::ISurface::SaveToHelper().

167  {
168  auto oldDir = gDirectory;
169  dir->cd();
170 
171  dir->cd();
172  TObjString("BayesianSurface").Write("type");
173 
175 
177 
178  if (oldDir)
179  oldDir->cd();
180  }
void SaveTo(TDirectory *dir) const
void SaveToHelper(TDirectory *dir) const
dir should already be the appropriate sub-directory
Definition: ISurface.cxx:189
TDirectory * dir
Definition: macro.C:5
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(), SaveTo(), and ana::FrequentistSurface::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
std::vector<const ISyst*> ana::BayesianMarginal::Systs ( ) const
inlineinherited

Definition at line 76 of file BayesianMarginal.h.

Referenced by ana::BayesianMarginal::LoadInto(), and ana::Bayesian1DMarginal::ToTH1().

76 { return fSysts; }
std::vector< const ISyst * > fSysts
Systs we want to marginalize to (see also fVars)
double ana::BayesianMarginal::ThresholdFromTH1 ( const TH1 *  pdf,
double  quantile 
)
staticprotectedinherited

Determine the threshold for the 'histogram' mode from the given histogram with a brute force method.

Definition at line 294 of file BayesianMarginal.cxx.

References std::abs(), om::cerr, allTimeWatchdog::endl, compare_h5_caf::idx, it, RunSnowGlobes::pdf, quantile(), and Scale().

Referenced by ana::BayesianMarginal::QuantileThreshold().

295  {
296  if (pdf->Integral() <= 0)
297  return std::numeric_limits<double>::signaling_NaN();
298 
299  // easy cases
300  if (quantile == 0)
301  return 1.0;
302  else if (quantile == 1)
303  return 0.0;
304  else if (quantile < 0 || quantile > 1)
305  {
306  std::cerr << "BayesianMarginal::ThresholdFromTH1(): requested quantile = " << quantile << " is not between 0 and 1. Abort" << std::endl;
307  abort();
308  }
309 
310  // check if the histogram I was passed was already unit normalized.
311  // if not, make a copy that is
312  if (std::abs(pdf->Integral() - 1) > 0.00001)
313  {
314  pdf = dynamic_cast<TH1*>(pdf->Clone());
315  const_cast<TH1*>(pdf)->Scale(1./pdf->Integral());
316  }
317 
318  // first make a list of the fractional bin contents.
319  std::vector<double> contents(pdf->GetNcells());
320  double overflowC = 0;
321  double underflowC = 0;
322  for (int binIdx = 0; binIdx <= pdf->GetNcells(); binIdx++)
323  {
324  double binC = pdf->GetBinContent(binIdx);
325  if (binC == 0)
326  continue; // don't go inserting useless 0s into the vector to make it much longer and slower to sort
327 
328  if (pdf->IsBinOverflow(binIdx))
329  overflowC += binC;
330  else if (pdf->IsBinUnderflow(binIdx))
331  underflowC += binC;
332 
333  contents.push_back(binC);
334  }
335  if ((1 - quantile) < underflowC || (1 - quantile) < overflowC)
336  {
337  std::cerr << "Warning: Histogram underflow or overflow contains more than requested quantile (" << quantile << ") of events." << std::endl;
338  std::cerr << " Bayesian credible region will be incorrect." << std::endl;
339  }
340 
341  // then do a search to find the value above which the integral corresponds to the desired quantile.
342  // begin at the specified quantile index as a first guess
343 
344  std::sort(contents.begin(), contents.end());
345  std::size_t idx = (1 - quantile) * (contents.size() - 1);
346  auto it = contents.begin();
347  decltype(it) oldIt = it;
348  std::advance(it, idx);
349  double oldSum = std::accumulate(it, contents.end(), 0.0);
350  // i.e., are we including too much, so we need to shrink the interval by moving towards the end?
351  bool forward = oldSum > quantile;
352  while (it != contents.begin() && it != contents.end())
353  {
354  // if we're exactly right, then we're good
355  if (oldSum == quantile)
356  return *it;
357 
358  std::advance(it, forward ? +1 : -1);
359  double newSum = std::accumulate(it, contents.end(), 0.0);
360 
361  // if we just crossed over the quantile boundary, we're also done.
362  if ((newSum - quantile) / (oldSum - quantile) < 0)
363  break;
364 
365  oldSum = newSum;
366  oldIt = it;
367  }
368 
369  // the average is not unbiased (it depends on the steepness of the distribution near the breakpoint)
370  // but it's probably closer than the values on either side?
371  return (*it + *oldIt) / 2;
372  }
set< int >::iterator it
OStream cerr
Definition: OStream.cxx:7
float abs(float number)
Definition: d0nt_math.hpp:39
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
while(!feof(fp))
simulatedPeak Scale(1/simulationNormalisationFactor)
std::unique_ptr< TH1 > ana::BayesianMarginal::ToHistogram ( const std::vector< Binning > &  bins) const
protectedinherited

Definition at line 375 of file BayesianMarginal.cxx.

References ana::assert(), b, om::cerr, allTimeWatchdog::endl, stan::math::exp(), Fill(), ana::BayesianMarginal::fMCMCSamples, ana::BayesianMarginal::fMode, ana::BayesianMarginal::fOrderedBrNames, ana::BayesianMarginal::fSysts, ana::BayesianMarginal::fVars, std::isnan(), ana::BayesianMarginal::kHistogram, ana::BayesianMarginal::kLLWgtdHistogram, ana::MakeTH1D(), ana::MakeTH2D(), ana::MCMCSamples::NumSamples(), runNovaSAM::ret, ana::MCMCSamples::SampleLL(), ana::MCMCSamples::SampleValue(), ana::MCMCSamples::SampleValues(), string, ana::UniqueName(), and PandAna.Demos.tute_pid_validation::var.

Referenced by BuildHist(), ana::Bayesian1DMarginal::QuantileRanges(), and ana::Bayesian1DMarginal::ToTH1().

376  {
377  // idea:
378  // * for mode==kHistogram:
379  // simply histogram the MCMC samples.
380  // * for mode==kLLWgtdHistogram:
381  // find the bin-averaged probability distribution
382  // by making the histogram weighted by exp(LL)
383  // and dividing it by the simple histogram of the events
385  assert(bins.size() == fOrderedBrNames.size());
386  assert(bins.size() > 0 && bins.size() < 3); // could port NOvA CAFAna 3D hist stuff if needed
387 
388  auto axisNameStr = std::accumulate(fOrderedBrNames.begin(), fOrderedBrNames.end(),
389  std::string(""),
390  [this](std::string s, const std::string & b)
391  {
392 
395  name = var->LatexName();
396  else if (auto syst = Registry<ISyst>::ShortNameToPtr(fOrderedBrNames[0], true))
397  name = syst->LatexName();
398  return s + ";" + name;
399  });
400  std::unique_ptr<TH1> ret;
401  if (bins.size() == 1)
402  ret.reset(MakeTH1D(UniqueName().c_str(), axisNameStr.c_str(), bins[0]));
403  else if (bins.size() == 2)
404  ret.reset(MakeTH2D(UniqueName().c_str(), axisNameStr.c_str(), bins[0], bins[1]));
405  // we haven't implemented MakeTH3D here. but won't worry about it for now
406  //else if (bins.size() == 3)
407  // ret.reset(MakeTH3D(UniqueName().c_str(), axisNameStr.c_str(), bins[0], bins[1], bins[2]));
408 
409  std::unique_ptr<TH1> denom;
410  if (bins.size() == 1)
411  denom = std::make_unique<TH1D>(*dynamic_cast<TH1D*>(ret.get()));
412  else if (bins.size() == 2)
413  denom = std::make_unique<TH2D>(*dynamic_cast<TH2D*>(ret.get()));
414  //else if (bins.size() == 3)
415  // denom = std::make_unique<TH3D>(*dynamic_cast<TH3D*>(ret.get()));
416  denom->SetName(UniqueName().c_str());
417 
418  // this could certainly be done less sloppily
419  // but I wanted to ensure we don't do lots of expensive lookups in std::maps
420  // inside the loop over samples, which could run into the millions of iterations
421  union VarOrSyst
422  {
423  const IFitVar * var;
424  const ISyst * syst;
425  };
426  std::vector<bool> isVar(fVars.size() + fSysts.size(), true);
427  std::vector<VarOrSyst> varsOrSysts(fVars.size() + fSysts.size());
428  for (std::size_t brIdx = 0; brIdx < fOrderedBrNames.size(); brIdx++)
429  {
430  const auto & brName = fOrderedBrNames[brIdx];
431  if (auto var = Registry<IFitVar>::ShortNameToPtr(brName, true))
432  varsOrSysts[brIdx].var = var;
433  else if (auto syst = Registry<ISyst>::ShortNameToPtr(brName, true))
434  {
435  isVar[brIdx] = false;
436  varsOrSysts[brIdx].syst = syst;
437  }
438  }
439  for (std::size_t sample = 0; sample < fMCMCSamples->NumSamples(); ++sample)
440  {
441  double logprob = fMCMCSamples->SampleLL(sample);
442  if (std::isnan(logprob))
443  std::cerr << "Warning: Encountered NaN log-probability in an MCMC sample. Other things will probably go wrong..." << std::endl;
444  std::map<const IFitVar*, double> varVals;
445  std::map<const ISyst*, double> systVals;
446  fMCMCSamples->SampleValues(sample, fVars, varVals, fSysts, systVals);
447  std::vector<double> vals;
448  for (std::size_t brIdx = 0; brIdx < varsOrSysts.size(); ++brIdx)
449  vals.push_back(
450  isVar[brIdx] ? fMCMCSamples->SampleValue(varsOrSysts[brIdx].var, sample)
451  : fMCMCSamples->SampleValue(varsOrSysts[brIdx].syst, sample));
452 
453  // could template this, but probably not worth the effort?
454  double numWgt = (fMode == MarginalMode::kHistogram) ? 1.0 : logprob;
455  if (bins.size() == 1)
456  ret->Fill(vals[0], numWgt);
457  else if (bins.size() == 2)
458  {
459  dynamic_cast<TH2*>(ret.get())->Fill(vals[0], vals[1], numWgt);
460  dynamic_cast<TH2*>(denom.get())->Fill(vals[0], vals[1]);
461  }
462  else if (bins.size() == 3)
463  {
464  dynamic_cast<TH3*>(ret.get())->Fill(vals[0], vals[1], vals[2], numWgt);
465  dynamic_cast<TH3*>(denom.get())->Fill(vals[0], vals[1], vals[2]);
466  }
467  }
468 
469  // now exponentiate the log-probs to get probabilities
471  {
472  ret->Divide(denom.get());
473  for (int binIdx = 0; binIdx < ret->GetNcells(); binIdx++)
474  {
475  // if it's *exactly* 0, this is a region where there were no samples altogether. just leave empty
476  if (ret->GetBinContent(binIdx) != 0.0)
477  ret->SetBinContent(binIdx, exp(ret->GetBinContent(binIdx)));
478  }
479  }
480  else if (fMode == MarginalMode::kHistogram)
481  {
482  if (ret->Integral() > 0)
483  ret->Scale(1. / ret->Integral());
484  }
485 
486  return ret;
487  }
const XML_Char * name
Definition: expat.h:151
std::vector< const ISyst * > fSysts
Systs we want to marginalize to (see also fVars)
std::size_t NumSamples() const
How many samples do we have?
Definition: MCMCSamples.h:135
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
Definition: UtilsExt.cxx:74
OStream cerr
Definition: OStream.cxx:7
double SampleValue(const IFitVar *var, std::size_t idx) const
Get the value of Var var for sample number .
Definition: MCMCSamples.h:169
std::vector< std::string > fOrderedBrNames
the variables/systs to be marginalized over, in axis order (x-axis, y-axis, ...)
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
void SampleValues(std::size_t idx, const std::vector< const ana::IFitVar * > &vars, std::map< const ana::IFitVar *, double > &varVals) const
Get the values of FitVars vars for sample number idx.
const XML_Char * s
Definition: expat.h:262
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
MarginalMode fMode
What mode are we using?
TH2D * MakeTH2D(const char *name, const char *title, const Binning &binsx, const Binning &binsy)
Utility function to avoid 4-way combinatorial explosion on the bin types.
Definition: UtilsExt.cxx:86
double SampleLL(std::size_t idx) const
Get the LL for sample number .
Definition: MCMCSamples.h:162
const Binning bins
Definition: NumuCC_CPiBin.h:8
const MCMCSamples * fMCMCSamples
Ref to the MCMCSamples. Will not be persisted.
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
std::vector< const IFitVar * > fVars
Variables we want to marginalize to (see also fSysts)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
static const T * ShortNameToPtr(const std::string &s, bool allowFail=false)
Definition: Registry.cxx:60
enum BeamMode string
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
std::vector<const IFitVar*> ana::BayesianMarginal::Vars ( ) const
inlineinherited

Definition at line 77 of file BayesianMarginal.h.

References ana::bins, dir, RunSnowGlobes::pdf, quantile(), and ana::SaveTo().

Referenced by ana::BayesianMarginal::LoadInto(), and ana::Bayesian1DMarginal::ToTH1().

77 { return fVars; }
std::vector< const IFitVar * > fVars
Variables we want to marginalize to (see also fSysts)

Member Data Documentation

double ana::ISurface::fBestFitX
protectedinherited
double ana::ISurface::fBestFitY
protectedinherited
double ana::ISurface::fBestLikelihood
protectedinherited
TH2F* ana::ISurface::fHist
protectedinherited
bool ana::ISurface::fLogX
protectedinherited
bool ana::ISurface::fLogY
protectedinherited
const MCMCSamples* ana::BayesianMarginal::fMCMCSamples
protectedinherited

Ref to the MCMCSamples. Will not be persisted.

Definition at line 99 of file BayesianMarginal.h.

Referenced by ana::BayesianMarginal::ToHistogram().

MarginalMode ana::BayesianMarginal::fMode
protectedinherited
std::vector<std::string> ana::BayesianMarginal::fOrderedBrNames
protectedinherited

the variables/systs to be marginalized over, in axis order (x-axis, y-axis, ...)

Definition at line 102 of file BayesianMarginal.h.

Referenced by ana::BayesianMarginal::BayesianMarginal(), and ana::BayesianMarginal::ToHistogram().

std::vector<double> ana::ISurface::fSeedValues
protectedinherited
const std::vector< std::pair< Quantile, double > > ana::BayesianMarginal::QuantileUpVals
staticinherited

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