QuantileCuts.cxx
Go to the documentation of this file.
2 #include "TH2.h"
3 
4 #include <iostream>
5 #include <iomanip>
6 
7 namespace ana{
8  /// \brief Returns a TH2D with xAxis as independentAxis and yAxis as quantileAxis
9  /// Quantiles in yAxis can be calculated from this using \ref ana::GetQuantileBins
11  const HistAxis& independentAxis,
12  const HistAxis& quantileAxis,
13  const Cut& cut,
14  const SpillCut& spillCut,
15  const SystShifts& shift,
16  const Var& wei){
17 
18  loader.SetSpillCut(spillCut);
19  Spectrum spectrum2D(loader, independentAxis, quantileAxis, cut, shift, wei);
20  loader.Go();
21  return spectrum2D.ToTH2(spectrum2D.POT());
22 
23  }
24 
25  /// \brief Returns a 2D vector
26  /// First index is the independentAxis bin number
27  /// Second index is the high bin edge for a quantile
28  /// [i][j-1] is jth quantile's low edge for ith independentAxis bin
29  /// [i][j] is jth quantile's high edge for ith independentAxis bin
30  std::vector<std::vector<double>> GetQuantileBins(TH2* quantileHist,
31  const HistAxis& independentAxis,
32  const HistAxis& quantileAxis,
33  const unsigned int& numQuantiles,
34  const bool verbose){
35  if(independentAxis.NDimensions() != 1 ||
36  quantileAxis.NDimensions() != 1){
37  std::cout << "Error: GetQuantileBins() not equipped for multi-dimensional axes." << std::endl;
38  abort();
39  }
40 
41  const Binning independentBins = independentAxis.GetBinnings()[0];
42  const Binning quantileBins = quantileAxis.GetBinnings()[0];
43 
44  if(verbose){
45  std::cout << "GetQuantileBins - Calculating Quantile bins" << std::endl;
46  std::cout << "independentAxis has " << independentBins.NBins() << " bins and is labelled \'" << independentAxis.GetLabels()[0] << "\'" << std::endl;
47  std::cout << "quantileAxis has " << quantileBins.NBins() << " bins and is labelled \'" << quantileAxis.GetLabels()[0] << "\'" << std::endl;
48  }
49 
50  std::vector<std::vector<double>> xBinQuantilesVec;
51 
52  //Copied from Luke Vinton
53  //bin=0 is underflow
54  //independentBins.NBins()+1 is overflow
55  for(int bin = 0; bin <= independentBins.NBins()+1; bin++){
56 
57  TH1* hTempProj = quantileHist->ProjectionY(Form("quantileHistProjYBin%d", bin), bin, bin);
58  hTempProj->Scale( (double)1 / (double)hTempProj->Integral());
59  double probs[numQuantiles];
60  double q[numQuantiles];
61  for(unsigned int i = 0; i<numQuantiles; i++){
62  probs[i] = (i+1) * ((double)1)/ ((double)numQuantiles);
63  }
64 
65  hTempProj->GetQuantiles(numQuantiles, q, probs);
66  //Need to return what the low bin edge is as well, so return a vector of size numQuantiles+1
67  std::vector< double > quantilesVec(numQuantiles+1);
68  quantilesVec[0] = quantileHist->GetYaxis()->GetBinLowEdge(1);
69  for(unsigned int i = 0; i<numQuantiles; i++){
70  quantilesVec[i+1] = q[i];
71  }
72 
73  xBinQuantilesVec.push_back(quantilesVec);
74  if(hTempProj) delete hTempProj;
75  }//loop over independentAxis bins
76 
77  //print out the bin edges
78  if(verbose){
79  for(unsigned int bin=0;bin<xBinQuantilesVec.size();bin++){
80  std::cout << "independentAxis bin: " << std::setw(4) << bin << " "
81  << "lowEdge: " << std::setw(10) << quantileHist->GetXaxis()->GetBinLowEdge(bin) << " "
82  << "centre: " << std::setw(10) << quantileHist->GetXaxis()->GetBinCenter(bin) << " "
83  << "quantile bin edges:";
84  for(unsigned int edge=0;edge<xBinQuantilesVec[bin].size();edge++) std::cout << " " << std::setw(10) << xBinQuantilesVec[bin][edge];
86  }
87  }
88 
89  return xBinQuantilesVec;
90  }
91 
92  /// \brief Returns a 2D vector
93  /// First index is the independentAxis bin number
94  /// Second index is the high bin edge for a quantile
95  /// [i][j-1] is jth quantile's low edge for ith independentAxis bin
96  /// [i][j] is jth quantile's high edge for ith independentAxis bin
97  std::vector<std::vector<double>> GetQuantileBins(SpectrumLoader& loader,
98  const HistAxis& independentAxis,
99  const HistAxis& quantileAxis,
100  const Cut& cut,
101  const unsigned int& numQuantiles,
102  const SpillCut& spillCut,
103  const SystShifts& shift,
104  const Var& wei,
105  const bool verbose){
106 
107  TH2* quantileHist = MakeQuantileHistogram(loader, independentAxis, quantileAxis, cut, spillCut, shift, wei);
108 
109  std::vector<std::vector<double>> xBinQuantilesVec = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);
110 
111  //clean up after ourselves
112  if(quantileHist) delete quantileHist;
113 
114  return xBinQuantilesVec;
115  }
116 
117  /// \brief Class to Generate a Quantile \ref ana::Cut that is specialised via the constructor arguments
118  /// \param independentAxis - each bin on this axis will be split into the quantiles (i.e. a range of quantileAxis bins)
119  /// \param quantileAxis - the axis to split into quantiles. This is done separately for each independentAxis bin
120  /// \param quantile - the quantile to pass, other quantiles fail
121  /// \param quantileBins - 2D array specifying quantileAxis quantile edges. 1st index is over independentAxis, second over quantile
122 
124  const HistAxis& quantileAxis,
125  const unsigned int& quantile,
126  const std::vector<std::vector<double>>& quantileBins):
127  fIndependentAxis(independentAxis), fQuantileAxis(quantileAxis), fQuantile(quantile), fQuantileBins(quantileBins)
128  {
129  if(independentAxis.NDimensions() != 1 ||
130  quantileAxis.NDimensions() != 1){
131  std::cout << "Error: QuantileCutGenerator not equipped for multi-dimensional axes." << std::endl;
132  abort();
133  }
134 
135  assert(fQuantileBins.size() > 0 && "fQuantileBins must have more than 0 bins in the independent variable");
136  //fQuantileBins[] has the TH1 binning schema
137  //fQuantileBins[].size = number of fQuantiles + 1 as we keep both the low edge and high edge
138  assert((fQuantile > 0 && fQuantile < fQuantileBins[0].size()) && "fQuantile must be > 0 and <= number of fQuantiles");
139  }
140 
142  double varX = fIndependentAxis.GetVars()[0](sr);
143  double varY = fQuantileAxis.GetVars()[0](sr);
144  unsigned int binX = (unsigned int)(fIndependentAxis.GetBinnings()[0].FindBin(varX));
145 
146  //fQuantileBins first index is independentAxis bin number, where 0 is underflow
147  assert(binX < fQuantileBins.size() && "binX >= fQuantileBins.size()");
148  std::vector<double> thisQuantilesBins = fQuantileBins[binX];
149  assert(fQuantile > 0 && fQuantile < thisQuantilesBins.size() && "fQuantile <= 0 || fQuantile >= thisQuantilesBins.size()");
150  double lowEdge=thisQuantilesBins[fQuantile-1];
151  double highEdge=thisQuantilesBins[fQuantile];
152  return varY >= lowEdge && varY < highEdge;
153  }
154 
155 
156  /// \brief Returns a Cut that will only pass events falling into a particular quantile of quantileAxis, where quantiles have been calculated for each independentAxis bin
157  /// Quantiles have been precalculated in 'quantileBins' (number of quantiles, bin edges etc.)
158  const Cut QuantileCut(const std::vector<std::vector<double>>& quantileBins,
159  const HistAxis& independentAxis,
160  const HistAxis& quantileAxis,
161  const unsigned int& quantile){
162  if(independentAxis.NDimensions() != 1 ||
163  quantileAxis.NDimensions() != 1){
164  std::cout << "Error: QuantileCut() not equipped for multi-dimensional axes." << std::endl;
165  abort();
166  }
167 
168  assert((quantileBins.size() > 0) && "quantileBins.size() must be > 0");
169  unsigned int numQuantiles = quantileBins[0].size();
170  assert((quantile > 0 && quantile <= numQuantiles) && "quantile must be > 0 and <= number of quantiles");
171 
172  return Cut(QuantileCutGenerator(independentAxis, quantileAxis, quantile, quantileBins));
173 
174  }
175 
176  /// \brief Returns a vector of Cuts, each one for a different quantile
177  /// An individual cut will only pass events falling into a particular quantile of quantileAxis, where quantiles have been calculated for each independentAxis bin
178  std::vector<Cut> QuantileCuts(SpectrumLoader& loader,
179  const HistAxis& independentAxis,
180  const HistAxis& quantileAxis,
181  const Cut& cut,
182  const unsigned int& numQuantiles,
183  const SpillCut& spillCut,
184  const SystShifts& shift,
185  const Var& wei,
186  const bool verbose){
187 
188  std::vector<std::vector<double>> quantileBins = GetQuantileBins(loader, independentAxis, quantileAxis, cut, numQuantiles, spillCut, shift, wei, verbose);
189 
190  std::vector<Cut> cutVec;
191  for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
192  cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
193  }
194  return cutVec;
195  }
196 
197  /// \breif: Do the same as the QuantileCuts function but taking in the TH2 instead of making it.
198  std::vector<Cut> QuantileCutsFromTH2(TH2* quantileHist,
199  const HistAxis& independentAxis,
200  const HistAxis& quantileAxis,
201  const unsigned int& numQuantiles,
202  const bool verbose){
203 
204  std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);
205 
206  std::vector<Cut> cutVec;
207  for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
208  cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
209  }
210  return cutVec;
211  }
212 
213  std::vector<Cut> QuantileAndPIDCutsFromTH2(TH2* quantileHist,
214  const HistAxis& independentAxis,
215  const HistAxis& quantileAxis,
216  const unsigned int& numQuantiles,
217  const bool rhc,
218  const bool verbose){
219 
220  std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist, independentAxis, quantileAxis, numQuantiles, verbose);
221 
222  std::vector<Cut> cutVec;
223  for(unsigned int quantile=1;quantile<=numQuantiles;quantile++){
224  cutVec.push_back(QuantileCut(quantileBins, independentAxis, quantileAxis, quantile));
225  }
226  if (cutVec.size()==4 && rhc){
227  cutVec[0] = cutVec[0] && k2018PIDs(0.47, 0.40, 0.80);
228  cutVec[1] = cutVec[1] && k2018PIDs(0.48, 0.40, 0.65);
229  cutVec[2] = cutVec[2] && k2018PIDs(0.56, 0.40, 0.60);
230  cutVec[3] = cutVec[3] && k2018PIDs(0.59, 0.65, 0.50);
231  }
232  else if (cutVec.size()==4 && !rhc){
233  cutVec[0] = cutVec[0] && k2018PIDs(0.56, 0.60, 0.55);
234  cutVec[1] = cutVec[1] && k2018PIDs(0.56, 0.60, 0.55);
235  cutVec[2] = cutVec[2] && k2018PIDs(0.58, 0.60, 0.55);
236  cutVec[3] = cutVec[3] && k2018PIDs(0.59, 0.60, 0.55);
237  }
238  else for(unsigned int quant=0;quant<numQuantiles;quant++) cutVec[quant] = cutVec[quant] && k2018PIDs();
239 
240  return cutVec;
241  }
242 
243 
244 
245 
246 }//namespace
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:226
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
unsigned int NDimensions() const
Definition: HistAxis.h:87
const std::vector< std::string > & GetLabels() const
Definition: HistAxis.h:90
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut QuantileCut(const std::vector< std::vector< double >> &quantileBins, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &quantile)
Returns a Cut that will only pass events falling into a particular quantile of quantileAxis, where quantiles have been calculated for each independentAxis bin Quantiles have been precalculated in &#39;quantileBins&#39; (number of quantiles, bin edges etc.)
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
Collect information describing the x-axis of an analysis histogram.
Definition: HistAxis.h:15
GenericCut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
std::vector< std::vector< double > > GetQuantileBins(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
Returns a 2D vector First index is the independentAxis bin number Second index is the high bin edge f...
std::vector< Cut > QuantileCuts(SpectrumLoader &loader, const HistAxis &independentAxis, const HistAxis &quantileAxis, const Cut &cut, const unsigned int &numQuantiles, const SpillCut &spillCut, const SystShifts &shift, const Var &wei, const bool verbose)
Returns a vector of Cuts, each one for a different quantile An individual cut will only pass events f...
void SetSpillCut(const SpillCut &cut)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Cut k2018PIDs(double bdtCut, double remidCut, double cvnCut)
std::vector< std::vector< double > > fQuantileBins
Definition: QuantileCuts.h:70
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:231
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
TH2 * MakeQuantileHistogram(SpectrumLoader &loader, const HistAxis &independentAxis, const HistAxis &quantileAxis, const Cut &cut, const SpillCut &spillCut, const SystShifts &shift, const Var &wei)
Returns a TH2D with xAxis as independentAxis and yAxis as quantileAxis Quantiles in yAxis can be calc...
OStream cout
Definition: OStream.cxx:6
const Cut cut
Definition: exporter_fd.C:30
bool operator()(const caf::SRProxy *sr) const
Defines the Cut this Instance generates. Returns true if the slice is in the quantile set in quantile...
Template for Cut and SpillCut.
Definition: Cut.h:15
std::vector< Cut > QuantileAndPIDCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool rhc, const bool verbose)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
int NBins() const
Definition: Binning.h:25
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
assert(nhit_max >=nhit_nbins)
QuantileCutGenerator(const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &quantile, const std::vector< std::vector< double >> &quantileBins)
Class to Generate a Quantile ana::Cut that is specialised via the constructor arguments.
static constexpr Double_t sr
Definition: Munits.h:164
const std::vector< Binning > & GetBinnings() const
Definition: HistAxis.h:91
Template for Var and SpillVar.
Definition: Var.h:16