Bayesian1DMarginal.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Binning.h"
4 #include "CAFAna/Core/IFitVar.h"
7 
8 #include "Utilities/func/Stan.h"
9 
10 #include "TH1D.h"
11 
12 #include <cmath>
13 
14 namespace ana
15 {
16  //----------------------------------------------------------------------
18  : BayesianMarginal(samples, {var}, mode)
19  {}
20 
21  //----------------------------------------------------------------------
23  : BayesianMarginal(samples, {syst}, mode)
24  {}
25 
26  //----------------------------------------------------------------------
28  {
29  if (fCachedBinning && bins == *fCachedBinning)
30  return *fCachedHist;
31 
32  fCachedBinning = std::make_unique<Binning>(bins);
34  if (Vars().size() > 0)
35  name = Vars()[0]->LatexName();
36  else if (Systs().size() > 0)
37  name = Systs()[0]->LatexName();
39  fCachedHist.reset(dynamic_cast<TH1D*>(ToHistogram({bins}).release()));
40  else if (fMode == MarginalMode::kKNN)
41  {
42  fCachedHist = std::make_unique<TH1D>(UniqueName().c_str(), (";"+name).c_str(), bins.NBins(), &bins.Edges()[0]);
43  for (int bin = 1; bin <= bins.NBins(); bin++)
44  fCachedHist->SetBinContent(bin, EstimateLLFromKNN({float(fCachedHist->GetBinCenter(bin))}));
45  }
46  return *fCachedHist;
47  }
48 
49 
50  //----------------------------------------------------------------------
51  std::vector<std::pair<double, double>> Bayesian1DMarginal::QuantileRanges(Quantile quantile, const Binning & bins) const
52  {
53 
54  assert(bins.NBins() >= 2);
55 
56  std::unique_ptr<TH1> hist;
58  hist = ToHistogram({bins});
59  double threshold = QuantileThreshold(quantile, hist.get());
60 
61  // walk through the bins.
62  // if we cross from above to below the threshold (or vice versa),
63  // take that as a boundary of a range
64  std::vector<std::pair<double, double>> ranges;
65  std::pair<double, double> pair;
66  bool above = false;
67  for (unsigned int bin = 0; bin < bins.Edges().size(); ++bin)
68  {
69 
70  double bnd, y;
72  {
73  // the uppermost bin edge isn't the right thing to use when we use bin centers, like for the hist
74  if (bin >= bins.Edges().size())
75  break;
76 
77  bnd = hist->GetBinCenter(bin);
78  y = hist->GetBinContent(bin);
79  }
80  else
81  {
82  bnd = bins.Edges()[bin];
83  y = EstimateLLFromKNN({float(bnd)});
84  }
85 
86  // only do anything if we are changing from being above threshold to below
87  // or vice versa
88  bool transition = (y >= threshold && !above) || (y < threshold && above)
89  || bin==0 || bin>=bins.Edges().size()-1;
90  if (!transition)
91  continue;
92 
93  // we're opening a new interval if we went from below to above threshold
94  bool opening = !above;
95  double &val = opening ? pair.first : pair.second;
96  val = bnd;
97 
98  // store the pair if we're closing it
99  if (!opening)
100  ranges.emplace_back(pair);
101 
102  above = y >= threshold;
103  }
104 
105  return ranges;
106  }
107 
108 }
const XML_Char * name
Definition: expat.h:151
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::unique_ptr< Binning > fCachedBinning
Bayesian1DMarginal(const MCMCSamples &samples, const IFitVar *var, MarginalMode mode=MarginalMode::kHistogram)
Build a 1D marginal distribution from the output of an MCMC fit.
std::vector< std::pair< double, double > > QuantileRanges(Quantile quantile, const Binning &bins) const
Return a list of ranges of the variable where LL is above the specified quantile. ...
std::vector< const IFitVar * > Vars() const
std::unique_ptr< TH1D > fCachedHist
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
double EstimateLLFromKNN(const std::vector< float > &vals) const
MarginalMode
How to construct marginal distributions.
Quantile
Type to specify quantile values, so that we don&#39;t run into floating-point matching issues...
float bin[41]
Definition: plottest35.C:14
std::vector< const ISyst * > Systs() const
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
const std::vector< double > & Edges() const
Definition: Binning.h:34
MarginalMode fMode
What mode are we using?
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The result of marginalizing over a set of Bayesian MCMC samples down to a few dimensions.
int NBins() const
Definition: Binning.h:29
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
double QuantileThreshold(Quantile quantile, const TH1 *pdf=nullptr) const
TH1D ToTH1(const Binning &bins) const
std::unique_ptr< TH1 > ToHistogram(const std::vector< Binning > &bins) const
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29