BayesianSurface.cxx
Go to the documentation of this file.
1 
2 #include "TH2F.h"
3 
4 #include "CAFAna/Core/Binning.h"
5 #include "CAFAna/Core/IFitVar.h"
9 
10 namespace ana
11 {
12  //----------------------------------------------------------------------
14  const FitAxis& xax,
15  const FitAxis& yax,
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  }
23 
24  //----------------------------------------------------------------------
26  const IFitVar * xvar, int nbinsx, double xmin, double xmax,
27  const IFitVar * yvar, int nbinsy, double ymin, double ymax,
29  : BayesianMarginal(samples, {xvar, yvar}, mode)
30  {
31  fLogX = fLogY = false;
32  BuildHist(samples, xvar, nbinsx, xmin, xmax, yvar, nbinsy, ymin, ymax);
33  }
34 
35  //----------------------------------------------------------------------
37  const ISyst * xsyst, int nbinsx, double xmin, double xmax,
38  const IFitVar * yvar, int nbinsy, double ymin, double ymax,
40  : BayesianMarginal(samples, {xsyst, yvar}, mode)
41  {
42  fLogX = fLogY = false;
43  BuildHist(samples, xsyst, nbinsx, xmin, xmax, yvar, nbinsy, ymin, ymax);
44  }
45 
46  //----------------------------------------------------------------------
48  const IFitVar *xvar, int nbinsx, double xmin, double xmax,
49  const ISyst *ysyst, int nbinsy, double ymin, double ymax,
51  : BayesianMarginal(samples, {xvar, ysyst}, mode)
52  {
53  fLogX = fLogY = false;
54  BuildHist(samples, xvar, nbinsx, xmin, xmax, ysyst, nbinsy, ymin, ymax);
55  }
56 
57  //----------------------------------------------------------------------
59  const ISyst *xsyst, int nbinsx, double xmin, double xmax,
60  const ISyst *ysyst, int nbinsy, double ymin, double ymax,
62  : BayesianMarginal(samples, {xsyst, ysyst}, mode)
63  {
64  fLogX = fLogY = false;
65  BuildHist(samples, xsyst, nbinsx, xmin, xmax, ysyst, nbinsy, ymin, ymax);
66  }
67 
68  //----------------------------------------------------------------------
69  template <typename SystOrVar1, typename SystOrVar2>
71  const SystOrVar1 *x, int nbinsx, double xmin, double xmax,
72  const SystOrVar2 *y, int nbinsy, double ymin, double ymax)
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  }
108  // explicitly instantiate the ones we need
109  template void BayesianSurface::BuildHist(const MCMCSamples & samples,
110  const IFitVar * x, int nbinsx, double xmin, double xmax,
111  const IFitVar * y, int nbinsy, double ymin, double ymax);
112  template void BayesianSurface::BuildHist(const MCMCSamples & samples,
113  const IFitVar * x, int nbinsx, double xmin, double xmax,
114  const ISyst * y, int nbinsy, double ymin, double ymax);
115  template void BayesianSurface::BuildHist(const MCMCSamples & samples,
116  const ISyst * x, int nbinsx, double xmin, double xmax,
117  const IFitVar * y, int nbinsy, double ymin, double ymax);
118  template void BayesianSurface::BuildHist(const MCMCSamples & samples,
119  const ISyst * x, int nbinsx, double xmin, double xmax,
120  const ISyst * y, int nbinsy, double ymin, double ymax);
121 
122  //----------------------------------------------------------------------
123  std::unique_ptr<BayesianSurface> BayesianSurface::LoadFrom(TDirectory *dir, const std::string& name)
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);
137  ISurface::FillSurfObj(*surf, dir);
138 
139  BayesianMarginal::LoadInto(surf.get(), dir);
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  }
151 
152  //----------------------------------------------------------------------
154  {
155  return Flat(QuantileThreshold(quantile, fHist), *this);
156  }
157 
158 
159  //----------------------------------------------------------------------
160  TH2 * BayesianSurface::QuantileSurface(double quantile, const MCMCSamples& mcmcsamples) const
161  {
162  return Flat(QuantileThreshold(quantile, mcmcsamples, fHist), *this);
163  }
164 
165  //----------------------------------------------------------------------
166  void BayesianSurface::SaveTo(TDirectory *dir) const
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  }
181 
182 
183 }
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
double BinCenterX(int bin) const
Definition: ISurface.cxx:244
static void LoadInto(BayesianMarginal *marg, TDirectory *dir)
void SaveTo(TDirectory *dir) const
TH2 * QuantileSurface(Quantile quantile) const
const IFitVar * var
Definition: FitAxis.h:17
double BinCenterY(int bin) const
Definition: ISurface.cxx:251
double SampleValue(const IFitVar *var, std::size_t idx) const
Get the value of Var var for sample number .
Definition: MCMCSamples.h:169
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
double EstimateLLFromKNN(const std::vector< float > &vals) const
static std::unique_ptr< BayesianSurface > LoadFrom(TDirectory *dir, const std::string &name)
bool logy
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
MarginalMode
How to construct marginal distributions.
const XML_Char int const XML_Char * value
Definition: expat.h:331
TH2F * fHist
Definition: ISurface.h:66
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
double fBestFitY
Definition: ISurface.h:65
void SaveToHelper(TDirectory *dir) const
dir should already be the appropriate sub-directory
Definition: ISurface.cxx:189
Quantile
Type to specify quantile values, so that we don&#39;t run into floating-point matching issues...
bool fLogY
Definition: ISurface.h:67
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?
bool logx
BayesianSurface()=default
double SampleLL(std::size_t idx) const
Get the LL for sample number .
Definition: MCMCSamples.h:162
The result of marginalizing over a set of Bayesian MCMC samples down to a few dimensions.
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
double fBestFitX
Definition: ISurface.h:64
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double QuantileThreshold(Quantile quantile, const TH1 *pdf=nullptr) const
Double_t ymin
Definition: plot.C:24
static void FillSurfObj(ISurface &surf, TDirectory *dir)
Definition: ISurface.cxx:214
std::unique_ptr< TH1 > ToHistogram(const std::vector< Binning > &bins) const
TH2 * Flat(double level, const ISurface &s)
Helper function for the gaussian approximation surfaces.
Definition: ISurface.cxx:232
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
void SaveTo(TDirectory *dir) const
double fBestLikelihood
Definition: ISurface.h:63
std::size_t BestFitSampleIdx() const
Retrieve the index of the best-fit sample (lowest LL), which can then be used with the SampleLL() or ...
Collect information describing the axis of a fit variable.
Definition: FitAxis.h:10
surf
Definition: demo5.py:35
enum BeamMode string