Binning.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Binning.h"
2 
3 #include "TDirectory.h"
4 #include "TH1.h"
5 #include "TObjString.h"
6 #include "TVectorD.h"
7 
8 namespace ana
9 {
10  int Binning::fgNextID = 0;
11 
12  //----------------------------------------------------------------------
14  : fID(-1)
15  {
16  }
17 
18  //----------------------------------------------------------------------
19  Binning Binning::SimpleHelper(int n, double lo, double hi,
20  const std::vector<std::string>& labels)
21  {
22  assert(labels.empty() || int(labels.size()) == n);
23 
24  Binning bins;
25  bins.fNBins = n;
26  bins.fMin = lo;
27  bins.fMax = hi;
28  bins.fEdges.resize(n+1);
29  for (int i = 0; i <= n; i++)
30  bins.fEdges[i] = lo + i*(hi-lo)/n;
31  bins.fLabels = labels;
32  bins.fIsSimple = true;
33 
34  return bins;
35  }
36 
37  //----------------------------------------------------------------------
38  Binning Binning::Simple(int n, double lo, double hi,
39  const std::vector<std::string>& labels)
40  {
41  Binning bins = SimpleHelper(n, lo, hi, labels);
42 
43  auto it = IDMap().find(bins);
44  if(it == IDMap().end()){
45  bins.fID = fgNextID++;
46  IDMap().emplace(bins, bins.fID);
47  }
48  else{
49  bins.fID = it->second;
50  }
51 
52  return bins;
53  }
54 
55  //----------------------------------------------------------------------
56  Binning Binning::LogUniform(int n, double lo, double hi)
57  {
58  std::vector<double> edges(n+1);
59  const double logSpacing = exp( (log(hi) - log(lo)) / n );
60  for (int i = 0; i <= n; ++i) {
61  if (i == 0) edges[i] = lo;
62  else edges[i] = logSpacing*edges[i-1];
63  }
64  return Custom(edges);
65  }
66 
67  //----------------------------------------------------------------------
68  Binning Binning::CustomHelper(const std::vector<double>& edges)
69  {
70  assert(edges.size() > 1);
71 
72  Binning bins;
73  bins.fEdges = edges;
74  bins.fNBins = edges.size()-1;
75  bins.fMin = edges.front();
76  bins.fMax = edges.back();
77  bins.fIsSimple = false;
78 
79  return bins;
80  }
81 
82  //----------------------------------------------------------------------
83  Binning Binning::Custom(const std::vector<double>& edges)
84  {
85  Binning bins = CustomHelper(edges);
86 
87  auto it = IDMap().find(bins);
88  if(it == IDMap().end()){
89  bins.fID = fgNextID++;
90  IDMap().emplace(bins, bins.fID);
91  }
92  else{
93  bins.fID = it->second;
94  }
95 
96  return bins;
97  }
98 
99  //----------------------------------------------------------------------
100  int Binning::FindBin(float x) const
101  {
102  // Treat anything outside [fMin, fMax) at Underflow / Overflow
103  if (x < fMin) return 0; // Underflow
104  if (x >= fMax) return fEdges.size(); // Overflow
105 
106  // Follow ROOT convention, first bin of histogram is bin 1
107 
108  if (this->IsSimple()){
109  double binwidth = (fMax - fMin) / fNBins;
110  int bin = (x - fMin) / binwidth + 1;
111  return bin;
112  }
113 
114  int bin =
115  std::lower_bound(fEdges.begin(), fEdges.end(), x) - fEdges.begin();
116  if (x == fEdges[bin]) bin++;
117  assert(bin >= 0 && bin < (int)fEdges.size());
118  return bin;
119  }
120 
121  //----------------------------------------------------------------------
123  {
124  // Evenly spaced binning
125  if(!ax->IsVariableBinSize()){
126  return Binning::Simple(ax->GetNbins(), ax->GetXmin(), ax->GetXmax());
127  }
128  else{
129  std::vector<double> edges(ax->GetNbins()+1);
130  ax->GetLowEdge(&edges.front());
131  edges[ax->GetNbins()] = ax->GetBinUpEdge(ax->GetNbins());
132 
133  return Binning::Custom(edges);
134  }
135  }
136 
137  // Can we give trueE bin a special ID?
138  //----------------------------------------------------------------------
140  {
141  // Osc P is ~sin^2(1/E). Difference in prob across a bin is ~dP/dE which
142  // goes as 1/E^2 times a trigonometric term depending on the parameters but
143  // bounded between -1 and +1. So bins should have width ~1/E^2. E_i~1/i
144  // achieves that.
145 
146  const int kNumTrueEnergyBins = 100;
147 
148  // N+1 bin low edges
149  std::vector<double> edges(kNumTrueEnergyBins+1);
150 
151  const double Emin = .5; // 500 MeV: there's really no events below there
152 
153  // How many edges to generate. Allow room for 0-Emin bin
154  const double N = kNumTrueEnergyBins-1;
155  const double A = N*Emin;
156 
157  edges[0] = 0;
158 
159  for(int i = 1; i <= N; ++i){
160  edges[kNumTrueEnergyBins-i] = A/i;
161  }
162 
163  edges[kNumTrueEnergyBins] = 120; // Replace the infinity that would be here
164 
165  return Binning::Custom(edges);
166  }
167 
168 
169  //----------------------------------------------------------------------
171  {
172  return Binning::LogUniform(400, 0.005, 5.);
173  }
174 
175  //----------------------------------------------------------------------
176  void Binning::SaveTo(TDirectory* dir) const
177  {
178  TDirectory* tmp = gDirectory;
179  dir->cd();
180 
181  TObjString("Binning").Write("type");
182 
183  TVectorD nminmax(3);
184 
185  nminmax[0] = fNBins;
186  nminmax[1] = fMin;
187  nminmax[2] = fMax;
188 
189  nminmax.Write("nminmax");
190 
191  TVectorD issimple(1);
192  issimple[0] = fIsSimple;
193  issimple.Write("issimple");
194 
195  TVectorD edges(fEdges.size());
196  for(unsigned int i = 0; i < fEdges.size(); ++i)
197  edges[i] = fEdges[i];
198 
199  edges.Write("edges");
200 
201  for(unsigned int i = 0; i < fLabels.size(); ++i)
202  TObjString(fLabels[i].c_str()).Write(TString::Format("label%d", i).Data());
203 
204  tmp->cd();
205  }
206 
207  //----------------------------------------------------------------------
208  std::unique_ptr<Binning> Binning::LoadFrom(TDirectory* dir)
209  {
210  TObjString* tag = (TObjString*)dir->Get("type");
211  assert(tag);
212  assert(tag->GetString() == "Binning");
213 
214  TVectorD* vMinMax = (TVectorD*)dir->Get("nminmax");
215  assert(vMinMax);
216 
217  Binning ret;
218 
219  const TVectorD* issimple = (TVectorD*)dir->Get("issimple");
220  if((*issimple)[0]){
221  ret = Binning::Simple((*vMinMax)[0],
222  (*vMinMax)[1],
223  (*vMinMax)[2]);
224  }
225  else{
226  const TVectorD* vEdges = (TVectorD*)dir->Get("edges");
227  std::vector<double> edges(vEdges->GetNrows());
228  for(int i = 0; i < vEdges->GetNrows(); ++i) edges[i] = (*vEdges)[i];
229 
230  ret = Binning::Custom(edges);
231  }
232 
233  for(unsigned int i = 0; ; ++i){
234  TObjString* s = (TObjString*)dir->Get(TString::Format("label%d", i).Data());
235  if(!s) break;
236  ret.fLabels.push_back(s->GetString().Data());
237  }
238 
239  return std::make_unique<Binning>(ret);
240  }
241 
242  //----------------------------------------------------------------------
243  bool Binning::operator==(const Binning& rhs) const
244  {
245  // NB don't look at ID here or in < because we use these in the maps below
246  // that are used to find the IDs in the first place
247  if(fIsSimple != rhs.fIsSimple) return false;
248  if(fNBins != rhs.fNBins) return false;
249 
250  if(fIsSimple){
251  return fMin == rhs.fMin && fMax == rhs.fMax;
252  }
253  else{
254  return fEdges == rhs.fEdges;
255  }
256  }
257 
258  //----------------------------------------------------------------------
259  bool Binning::operator<(const Binning& rhs) const
260  {
261  if(fIsSimple != rhs.fIsSimple) return fIsSimple < rhs.fIsSimple;
262  if(fNBins != rhs.fNBins) return fNBins < rhs.fNBins;
263 
264  if(fIsSimple){
265  return std::make_tuple(fMin, fMax) < std::make_tuple(rhs.fMin, rhs.fMax);
266  }
267  else{
268  return fEdges < rhs.fEdges;
269  }
270  }
271 
272  //----------------------------------------------------------------------
273  std::map<Binning, int>& Binning::IDMap()
274  {
275  static std::map<Binning, int> ret;
276  return ret;
277  }
278 
279 
280 
281  /// Default true-energy bin edges
283 
284  /// Default true-energy bin edges
286 
287  // These are just guesses
288  // No optimization done for these binnings.
291 
292  /// Energy binnings used in Ana01 for nus extrapolation
294 
295  /// Energy binnings used in Nus17 and planning for Nus18 for nus analysis
297 
298  /// Energy binnings for Nus18 for nus analysis
300 
301  /// \brief Binning for plotting remid attractively
302  ///
303  /// There are 81 possible values, from 0/80 to 80/80. Want them centered in
304  /// their bins.
305  const Binning kRemidBinning = Binning::Simple(81, -1/160., 1+1/160.);
306 
307  /// The energy part of the SA 2D binning
309 
310  /// \brief Nonlinear binning for LEM
311  ///
312  /// There are 3 bins. The boundary of the highest one is set at the
313  /// value of the S/Sqrt(B) cut
314 
315  const Binning kLEMNLBinning = Binning::Custom({0.47,0.65,0.80,1});
316 
317  /// \brief Nonlinear binning for LID
318  ///
319  /// There are 3 bins. The boundary of the highest one is set at the
320  /// value of the S/Sqrt(B) cut
321 
322  const Binning kLIDNLBinning = Binning::Custom({0.63,0.8,0.95,2});
323 
324  /// \brief Nonlinear binning for CVN
325  ///
326  /// There are 3 bins. The boundary of the highest one is set at the
327  /// value of the S/Sqrt(B) cut
328 
329  const Binning kCVNNLBinning = Binning::Custom({0.75,0.87,0.95,1});
330 
331  /// 10 energy bins x 3 pid bins is 30 bins in 1D
332  const Binning kNue2DBinning = Binning::Simple(30,0.0,30.0);
333 
334  /// 8 Cos theta bins. Theta is the angle between the muon and the incident
335  /// neutrino. Based on angular resolution from Numu CC Inclusive Cross
336  /// section studies
337  const Binning kNumuCCIncluXSecCosThetaBinning = Binning::Custom({0.0, 0.75, 0.8, 0.85, 0.88, 0.91, 0.94, 0.97, 1.001});
338 
339 
340  /// \brief Optimised binning for numuCCE from L. Vinton. See docdb 16332.
341  /// This was close to 'custC' in that talk, with a slight modification
342  /// agreed at the collaboration meeting
343  const Binning kNumuCCEOptimisedBinning = Binning::Custom( { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 } );
344 
345  /// NC custom binning
346  const Binning kNCNDBinning = Binning::Custom( { 0, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 10, 20 } );
347  const Binning kNCFDBinning = Binning::Custom( { 0, 1, 1.25, 1.5, 1.75, 2, 2.5, 3.25, 4, 5, 6, 8, 10, 20 } );
348 
349 }
TSpline3 lo("lo", xlo, ylo, 12,"0")
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: CutFlow_header.h:5
set< int >::iterator it
bool operator<(const Binning &rhs) const
Definition: Binning.cxx:259
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:305
Binning TrueLOverTrueEBins()
Default trueL Over true-energy bin edges.
Definition: Binning.cxx:170
const Binning kTrueLOverTrueEBins
Default true-energy bin edges.
Definition: Binning.cxx:285
const Binning kCVNNLBinning
Nonlinear binning for CVN.
Definition: Binning.cxx:329
static Binning FromTAxis(const TAxis *ax)
Definition: Binning.cxx:122
static std::map< Binning, int > & IDMap()
Definition: Binning.cxx:273
std::vector< std::string > fLabels
Definition: Binning.h:54
static Binning CustomHelper(const std::vector< double > &edges)
Definition: Binning.cxx:68
const Binning kBEN2020Binning
Definition: Binning.cxx:290
Float_t tmp
Definition: plot.C:36
const Binning kNueSAEnergyBinning
The energy part of the SA 2D binning.
Definition: Binning.cxx:308
std::vector< double > fEdges
Definition: Binning.h:53
bool fIsSimple
Definition: Binning.h:57
bool IsSimple() const
Definition: Binning.h:29
static int fgNextID
The next ID that hasn&#39;t yet been assigned.
Definition: Binning.h:61
const XML_Char * s
Definition: expat.h:262
TSpline3 hi("hi", xhi, yhi, 18,"0")
const Binning kNus18EnergyBinning
Energy binnings for Nus18 for nus analysis.
Definition: Binning.cxx:299
bool operator==(const Binning &rhs) const
Definition: Binning.cxx:243
const Binning kNumuEnergyBinning
Definition: Binning.cxx:289
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
void SaveTo(TDirectory *dir) const
Definition: Binning.cxx:176
int FindBin(float x) const
Definition: Binning.cxx:100
const Binning kTrueEnergyBins
Default true-energy bin edges.
Definition: Binning.cxx:282
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:83
float bin[41]
Definition: plottest35.C:14
const Binning kNue2DBinning
10 energy bins x 3 pid bins is 30 bins in 1D
Definition: Binning.cxx:332
const Binning bins
Definition: NumuCC_CPiBin.h:8
double fMax
Definition: Binning.h:56
static const double A
Definition: Units.h:82
const Binning kNCNDBinning
NC custom binning.
Definition: Binning.cxx:346
Binning TrueEnergyBins()
Default true-energy bin edges.
Definition: Binning.cxx:139
const Binning kLIDNLBinning
Nonlinear binning for LID.
Definition: Binning.cxx:322
static Binning LogUniform(int n, double lo, double hi)
Definition: Binning.cxx:56
TDirectory * dir
Definition: macro.C:5
const Binning kNumuCCIncluXSecCosThetaBinning
Definition: Binning.cxx:337
const Binning kLEMNLBinning
Nonlinear binning for LEM.
Definition: Binning.cxx:315
const Binning kNumuCCEOptimisedBinning
Optimised binning for numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39; in that talk...
Definition: Binning.cxx:343
const Binning kNCFDBinning
Definition: Binning.cxx:347
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
double fMin
Definition: Binning.h:56
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
static std::unique_ptr< Binning > LoadFrom(TDirectory *dir)
Definition: Binning.cxx:208
const Binning kNus17EnergyBinning
Energy binnings used in Nus17 and planning for Nus18 for nus analysis.
Definition: Binning.cxx:296
int fNBins
Definition: Binning.h:55
const Binning kNCDisappearanceEnergyBinning
Energy binnings used in Ana01 for nus extrapolation.
Definition: Binning.cxx:293
gm Write()
static Binning SimpleHelper(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:19