Spectrum.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
2 
3 #include "CAFAna/Core/Ratio.h"
4 #include "CAFAna/Core/Stan.h"
5 
6 #include "TDirectory.h"
7 #include "TH2.h"
8 #include "TH3.h"
9 #include "THnSparse.h"
10 #include "TObjString.h"
11 #include "TRandom3.h"
12 
13 #include <cassert>
14 #include <iostream>
15 
16 namespace
17 {
18  namespace util{template<class T> T sqr(const T& x){return x*x;}}
19 }
20 
21 namespace ana
22 {
23  //----------------------------------------------------------------------
25  : fHist(Hist::Uninitialized()), fPOT(0), fLivetime(0), fAxis(axis)
26  {
27  const Binning bins1D = fAxis.GetBins1D();
28 
29  if(sparse == kSparse){
30  fHist = Hist::ZeroSparse(bins1D.NBins());
31  }
32  else{
33  fHist = Hist::Zero(bins1D.NBins());
34  }
35  }
36 
37  //----------------------------------------------------------------------
38  Spectrum::Spectrum(Eigen::ArrayXd&& h,
39  const LabelsAndBins& axis,
40  double pot, double livetime)
41  : fHist(Hist::Adopt(std::move(h))), fPOT(pot), fLivetime(livetime), fAxis(axis)
42  {
43  }
44 
45  //----------------------------------------------------------------------
47  const LabelsAndBins& axis,
48  double pot, double livetime)
49  : fHist(Hist::AdoptStan(std::move(h))), fPOT(pot), fLivetime(livetime), fAxis(axis)
50  {
51  }
52 
53  //----------------------------------------------------------------------
55  fHist(rhs.fHist),
56  fPOT(rhs.fPOT),
57  fLivetime(rhs.fLivetime),
58  fAxis(rhs.fAxis)
59  {
60  assert(rhs.fReferences.empty()); // Copying with pending loads is unexpected
61  }
62 
63  //----------------------------------------------------------------------
65  fHist(std::move(rhs.fHist)),
66  fPOT(rhs.fPOT),
67  fLivetime(rhs.fLivetime),
68  fAxis(rhs.fAxis)
69  {
70  std::swap(fReferences, rhs.fReferences);
71  for(Spectrum** ref: fReferences) *ref = this;
72  }
73 
74  //----------------------------------------------------------------------
76  {
77  if(this == &rhs) return *this;
78 
79  fHist = rhs.fHist;
80  fPOT = rhs.fPOT;
81  fLivetime = rhs.fLivetime;
82  fAxis = rhs.fAxis;
83 
84  assert(fReferences.empty()); // Copying with pending loads is unexpected
85 
86  return *this;
87  }
88 
89  //----------------------------------------------------------------------
91  {
92  if(this == &rhs) return *this;
93 
94  fHist = std::move(rhs.fHist);
95  fPOT = rhs.fPOT;
96  fLivetime = rhs.fLivetime;
97  fAxis = rhs.fAxis;
98 
99  std::swap(fReferences, rhs.fReferences);
100  for(Spectrum** ref: fReferences) *ref = this;
101 
102  return *this;
103  }
104 
105  //----------------------------------------------------------------------
107  {
108  // Unregister self from anything that might still want to fill us
109  for(Spectrum** ref: fReferences) *ref = 0;
110  }
111 
112  //----------------------------------------------------------------------
113  TH1D* Spectrum::ToTH1(double exposure,
114  EExposureType expotype,
115  EBinType bintype) const
116  {
117  // Could have a file temporarily open
118  DontAddDirectory guard;
119 
120  TH1D* ret = fHist.ToTH1(fAxis.GetBins1D());
121 
122  ret->GetXaxis()->SetTitle(fAxis.GetLabel1D().c_str());
123  ret->GetYaxis()->SetTitle("Events");
124 
125  if(expotype == kPOT){
126  const double pot = exposure;
127  if(fPOT){
128  ret->Scale(pot/fPOT);
129  }
130  else{
131  // Allow zero POT if there are also zero events
132  if(ret->Integral() > 0){
133  std::cout << "Error: Spectrum with " << ret->Integral()
134  << " entries has zero POT, no way to scale to "
135  << exposure << " POT.";
136  if(fLivetime > 0){
137  std::cout << " Spectrum has " << fLivetime << " seconds livetime. "
138  << "Did you mean to pass kLivetime to ToTH1()?";
139  }
140  std::cout << std::endl;
141  abort();
142  }
143  }
144  }
145  if(expotype == kLivetime){
146  const double livetime = exposure;
147  if(fLivetime){
148  ret->Scale(livetime/fLivetime);
149  }
150  else{
151  // Allow zero exposure if there are also zero events
152  if(ret->Integral() > 0){
153  std::cout << "Error: Spectrum with " << ret->Integral()
154  << " entries has zero livetime, no way to scale to "
155  << livetime << " seconds.";
156  if(fPOT > 0){
157  std::cout << " Spectrum has " << fPOT << " POT. "
158  << "Did you mean to pass kPOT to ToTH1()?";
159  }
160  std::cout << std::endl;
161  abort();
162  }
163  }
164  }
165 
166  if(bintype == kBinDensity) ret->Scale(1, "width");
167 
168  // Allow GetMean() and friends to work even if this histogram never had any
169  // explicit Fill() calls made.
170  if(ret->GetEntries() == 0) ret->SetEntries(1);
171 
172  return ret;
173  }
174 
175  //----------------------------------------------------------------------
176  TH1D* Spectrum::ToTH1(double exposure, Color_t col, Style_t style,
177  EExposureType expotype,
178  EBinType bintype) const
179  {
180  // Could have a file temporarily open
181  DontAddDirectory guard;
182 
183  TH1D* ret = ToTH1(exposure, expotype, bintype);
184 
185  ret->SetLineColor(col);
186  ret->SetMarkerColor(col);
187  ret->SetLineStyle(style);
188 
189  return ret;
190  }
191 
192  //----------------------------------------------------------------------
193  TH2* Spectrum::ToTH2(double exposure, EExposureType expotype, EBinType bintype) const
194  {
195  if(NDimensions() != 2){
196  std::cout << "Error: This Spectrum does not appear to be 2D." << std::endl;
197  abort();
198  }
199 
200  TH2* ret = ana::ToTH2(*this, exposure, expotype, GetBinnings()[0], GetBinnings()[1]);
201 
202  ret->GetXaxis()->SetTitle(GetLabels()[0].c_str());
203  ret->GetYaxis()->SetTitle(GetLabels()[1].c_str());
204 
205  if(bintype == kBinDensity) ret->Scale(1, "width");
206 
207  // Allow GetMean() and friends to work even if this histogram never had any
208  // explicit Fill() calls made.
209  if(ret->GetEntries() == 0) ret->SetEntries(1);
210 
211  return ret;
212  }
213 
214  //----------------------------------------------------------------------
215  TH3* Spectrum::ToTH3(double exposure, EExposureType expotype, EBinType bintype) const
216  {
217  if(NDimensions() != 3){
218  std::cout << "Error: This Spectrum does not appear to be 3D." << std::endl;
219  abort();
220  }
221 
222  TH3* ret = ana::ToTH3(*this, exposure, expotype,
223  GetBinnings()[0], GetBinnings()[1], GetBinnings()[2]);
224 
225  ret->GetXaxis()->SetTitle(GetLabels()[0].c_str());
226  ret->GetYaxis()->SetTitle(GetLabels()[1].c_str());
227  ret->GetZaxis()->SetTitle(GetLabels()[2].c_str());
228 
229  if(bintype == kBinDensity) ret->Scale(1, "width");
230 
231  // Allow GetMean() and friends to work even if this histogram never had any
232  // explicit Fill() calls made.
233  if(ret->GetEntries() == 0) ret->SetEntries(1);
234 
235  return ret;
236  }
237 
238  //----------------------------------------------------------------------
239  TH1* Spectrum::ToTHX(double exposure, EExposureType expotype, EBinType bintype) const
240  {
241  if(NDimensions() == 2) return ToTH2(exposure, expotype, bintype);
242  if(NDimensions() == 3) return ToTH3(exposure, expotype, bintype);
243  return ToTH1(exposure, expotype, bintype);
244  }
245 
246  //----------------------------------------------------------------------
247  Eigen::ArrayXd Spectrum::GetEigen(double exposure, EExposureType expotype) const
248  {
249  if(expotype == kPOT)
250  return (exposure/fPOT) * fHist.GetEigen();
251  else
252  return (exposure/fLivetime) * fHist.GetEigen();
253  }
254 
255  //----------------------------------------------------------------------
256  Eigen::ArrayXstan Spectrum::GetEigenStan(double exposure, EExposureType expotype) const
257  {
258  if(expotype == kPOT)
259  return (exposure/fPOT) * fHist.GetEigenStan();
260  else
261  return (exposure/fLivetime) * fHist.GetEigenStan();
262  }
263 
264  //----------------------------------------------------------------------
265  void Spectrum::Scale(double c)
266  {
267  fHist.Scale(c);
268  }
269 
270  //----------------------------------------------------------------------
272  {
273  fHist.Scale(c);
274  }
275 
276  //----------------------------------------------------------------------
277  double Spectrum::Integral(double exposure, double* err,
278  EExposureType expotype) const
279  {
280  const double ratio = (expotype == kPOT) ? exposure/fPOT : exposure/fLivetime;
281 
282  if(err){
283  *err = 0;
284 
285  for(int i = 0; i < fHist.GetNbinsX()+2; ++i){
286  *err += util::sqr(fHist.GetBinError(i));
287  }
288  *err = sqrt(*err) * ratio;
289  }
290 
291  return fHist.Integral() * ratio;
292  }
293 
294  //----------------------------------------------------------------------
295  double Spectrum::Mean() const
296  {
297  const Binning bins = fAxis.GetBins1D();
298 
299  if(fHist.GetBinContent(0) != 0){
300  std::cout << "Spectrum::Mean(): Warning ignoring underflow bin content " << fHist.GetBinContent(0) << std::endl;
301  }
302 
303  if(fHist.GetBinContent(bins.NBins()+1) != 0){
304  std::cout << "Spectrum::Mean(): Warning ignoring overflow bin content " << fHist.GetBinContent(bins.NBins()+1) << std::endl;
305  }
306 
307  double mean = 0;
308  double W = 0;
309  for(int i = 1; i <= bins.NBins(); ++i){
310  const double w = fHist.GetBinContent(i);
311  W += w;
312  const double x0 = bins.Edges()[i-1];
313  const double x1 = bins.Edges()[i];
314  mean += w * (x0+x1)/2;
315  }
316 
317  return mean/W;
318  }
319 
320  //----------------------------------------------------------------------
321  void Spectrum::Fill(double x, double w)
322  {
323  fHist.Fill(fAxis.GetBins1D(), x, w);
324  // TODO Pull binning out of Hist entirely and just update an index?
325  }
326 
327  //----------------------------------------------------------------------
328  Spectrum Spectrum::MockData(double pot, int seed) const
329  {
330  Spectrum ret = FakeData(pot);
331 
332  TRandom3 rnd(seed); // zero seeds randomly
333 
334  for(int i = 0; i < ret.fHist.GetNbinsX()+2; ++i){
335  ret.fHist.SetBinContent(i, rnd.Poisson(ret.fHist.GetBinContent(i)));
336  }
337 
338  // Drop old errors, which are based on the MC statistics, and create new
339  // ones that are based on the prediction for the data
340  ret.fHist.ResetErrors();
341 
342  return ret;
343  }
344 
345  //----------------------------------------------------------------------
347  {
348  Spectrum ret = *this;
349  if(fPOT > 0) ret.fHist.Scale(pot/fPOT);
350  ret.fPOT = pot;
351 
352  // Drop old errors, which are based on the MC statistics, and create new
353  // ones that are based on the prediction for the data
354  ret.fHist.ResetErrors();
355 
356  return ret;
357  }
358 
359  //----------------------------------------------------------------------
361  {
362  Spectrum ret = *this;
363 
364  if(fPOT > 0) ret.fHist.Scale(pot/fPOT);
365  ret.fPOT = pot;
366 
367  // overwrite livetime for fake data
368  // PlusEqualsHelper will take over the rest
369  ret.fLivetime = livetime;
370 
371  ret.fHist.ResetErrors();
372 
373  return ret;
374  }
375 
376  //----------------------------------------------------------------------
378  {
379  return AsimovData(pot);
380  }
381 
382  //----------------------------------------------------------------------
383  Spectrum Spectrum::FakeData(double pot, double livetime) const
384  {
385  return AsimovData(pot, livetime);
386  }
387 
388  //----------------------------------------------------------------------
390  {
391  fHist.Reset();
392  }
393 
394  //----------------------------------------------------------------------
396  {
397  fReferences.erase(ref);
398  }
399 
400  //----------------------------------------------------------------------
402  {
403  fReferences.insert(ref);
404  }
405 
406  //----------------------------------------------------------------------
408  {
409  // In this case it would be OK to have no POT/livetime
410  if(rhs.fHist.Initialized() && rhs.fHist.Integral() == 0) return *this;
411 
412  if((!fPOT && !fLivetime) || (!rhs.fPOT && !rhs.fLivetime)){
413  std::cout << "Error: can't sum Spectrum with no POT or livetime: "
414  << fPOT << " " << rhs.fPOT << " " << fLivetime << " " << rhs.fLivetime
415  << std::endl;
416  abort();
417  }
418 
419  if(!fLivetime && !rhs.fPOT){
420  std::cout << "Error: can't sum Spectrum with POT ("
421  << fPOT << ") but no livetime and Spectrum with livetime ("
422  << rhs.fLivetime << " sec) but no POT." << std::endl;
423  abort();
424  }
425 
426  if(!fPOT && !rhs.fLivetime){
427  std::cout << "Error: can't sum Spectrum with livetime ("
428  << fLivetime << " sec) but no POT and Spectrum with POT ("
429  << rhs.fPOT << ") but no livetime." << std::endl;
430  abort();
431  }
432 
433  // And now there are still a bunch of good cases to consider
434 
435  if(fPOT && rhs.fPOT){
436  // Scale by POT when possible
437  fHist.Add(rhs.fHist, sign*fPOT/rhs.fPOT);
438 
439  if(fLivetime && rhs.fLivetime){
440  // If POT/livetime ratios match, keep regular lifetime, otherwise zero
441  // it out.
442  if(AlmostEqual(fLivetime*rhs.fPOT, rhs.fLivetime*fPOT))
443  fLivetime = 0;
444  }
445  if(!fLivetime && rhs.fLivetime){
446  // If the RHS has a livetime and we don't, copy it in (suitably scaled)
447  fLivetime = rhs.fLivetime * fPOT/rhs.fPOT;
448  }
449  // Otherwise, keep our own livetime (if any)
450 
451  return *this;
452  }
453 
454  if(fLivetime && rhs.fLivetime){
455  // Scale by livetime, the only thing in common
456  fHist.Add(rhs.fHist, sign*fLivetime/rhs.fLivetime);
457 
458  if(!fPOT && rhs.fPOT){
459  // If the RHS has a POT and we don't, copy it in (suitably scaled)
460  fPOT = rhs.fPOT * fLivetime/rhs.fLivetime;
461  }
462  // Otherwise, keep our own POT (if any)
463 
464  return *this;
465  }
466 
467  // That should have been all the cases. I definitely want to know what
468  // happened if it wasn't.
469  std::cout << "Spectrum::operator+=(). How did we get here? "
470  << fPOT << " " << fLivetime << " "
471  << rhs.fPOT << " " << rhs.fLivetime << std::endl;
472  abort();
473  }
474 
475  //----------------------------------------------------------------------
477  {
478  return PlusEqualsHelper(rhs, +1);
479  }
480 
481  //----------------------------------------------------------------------
483  {
484  Spectrum ret = *this;
485  ret += rhs;
486  return ret;
487  }
488 
489  //----------------------------------------------------------------------
491  {
492  return PlusEqualsHelper(rhs, -1);
493  }
494 
495  //----------------------------------------------------------------------
497  {
498  Spectrum ret = *this;
499  ret -= rhs;
500  return ret;
501  }
502 
503  //----------------------------------------------------------------------
505  {
506  fHist.Multiply(rhs.fHist);
507  return *this;
508  }
509 
510  //----------------------------------------------------------------------
512  {
513  Spectrum ret = *this;
514  ret *= rhs;
515  return ret;
516  }
517 
518  //----------------------------------------------------------------------
520  {
521  fHist.Divide(rhs.fHist);
522  return *this;
523  }
524 
525  //----------------------------------------------------------------------
527  {
528  Spectrum ret = *this;
529  ret /= rhs;
530  return ret;
531  }
532 
533  //----------------------------------------------------------------------
534  void Spectrum::SaveTo(TDirectory* dir, const std::string& name) const
535  {
536  TDirectory* tmp = gDirectory;
537 
538  dir = dir->mkdir(name.c_str()); // switch to subdir
539  dir->cd();
540 
541  TObjString("Spectrum").Write("type");
542 
544  TH1D hPot("", "", 1, 0, 1);
545  hPot.Fill(.5, fPOT);
546  hPot.Write("pot");
547  TH1D hLivetime("", "", 1, 0, 1);
548  hLivetime.Fill(.5, fLivetime);
549  hLivetime.Write("livetime");
550 
551  for(unsigned int i = 0; i < NDimensions(); ++i){
552  TObjString(GetLabels()[i].c_str()).Write(TString::Format("label%d", i).Data());
553  GetBinnings()[i].SaveTo(dir, TString::Format("bins%d", i).Data());
554  }
555 
556  dir->Write();
557  delete dir;
558 
559  tmp->cd();
560  }
561 
562  //----------------------------------------------------------------------
563  std::unique_ptr<Spectrum> Spectrum::LoadFrom(TDirectory* dir, const std::string& name)
564  {
565  dir = dir->GetDirectory(name.c_str()); // switch to subdir
566  assert(dir);
567 
568  DontAddDirectory guard;
569 
570  TObjString* tag = (TObjString*)dir->Get("type");
571  assert(tag);
572  assert(tag->GetString() == "Spectrum");
573  delete tag;
574 
575  TH1* hPot = (TH1*)dir->Get("pot");
576  assert(hPot);
577  TH1* hLivetime = (TH1*)dir->Get("livetime");
578  assert(hLivetime);
579 
580  std::vector<std::string> labels;
581  std::vector<Binning> bins;
582  for(int i = 0; ; ++i){
583  const std::string subname = TString::Format("bins%d", i).Data();
584  TDirectory* subdir = dir->GetDirectory(subname.c_str());
585  if(!subdir) break;
586  delete subdir;
587  bins.push_back(*Binning::LoadFrom(dir, subname));
588  TObjString* label = (TObjString*)dir->Get(TString::Format("label%d", i));
589  labels.push_back(label ? label->GetString().Data() : "");
590  delete label;
591  }
592 
593  std::unique_ptr<Spectrum> ret = std::make_unique<Spectrum>(Eigen::ArrayXd(), LabelsAndBins(labels, bins), hPot->GetBinContent(1), hLivetime->GetBinContent(1));
594  ret->fHist = Hist::FromDirectory(dir);
595 
596  delete hPot;
597  delete hLivetime;
598 
599  delete dir;
600 
601  return ret;
602  }
603 }
static std::unique_ptr< Binning > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Binning.cxx:230
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:263
const XML_Char * name
Definition: expat.h:151
void SetBinContent(int i, double x)
Definition: Hist.cxx:363
void Write(const Binning &bins) const
Definition: Hist.cxx:583
Filter events based on their run/event numbers.
const Eigen::ArrayXd & GetEigen() const
Definition: Hist.h:45
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:193
Spectrum operator*(const Ratio &rhs) const
Definition: Spectrum.cxx:511
void Add(const Hist &rhs, double scale=1)
Definition: Hist.cxx:455
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
void AddLoader(Spectrum **)
Definition: Spectrum.cxx:401
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
LabelsAndBins fAxis
Definition: Spectrum.h:290
double Integral() const
Definition: Hist.cxx:252
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:176
const SpillVar fPOT([](const caf::SRSpillProxy *spill){return spill->spillpot;})
subdir
Definition: cvnie.py:7
void ResetErrors()
Definition: Hist.cxx:342
Spectrum & PlusEqualsHelper(const Spectrum &rhs, int sign)
Helper for operator+= and operator-=.
Definition: Spectrum.cxx:407
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
Spectrum & operator+=(const Spectrum &rhs)
Definition: Spectrum.cxx:476
void RemoveLoader(Spectrum **)
Definition: Spectrum.cxx:395
Definition: Hist.h:21
std::set< Spectrum ** > fReferences
Things that point at this Spectrum. Maintained by SpectrumLoader.
Definition: Spectrum.h:288
static Hist FromDirectory(TDirectory *dir)
Definition: Hist.cxx:158
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
int GetNbinsX() const
Definition: Hist.cxx:228
std::vector< double > Spectrum
Definition: Constants.h:746
Divide bin contents by bin widths.
Definition: UtilsExt.h:31
var_value< double > var
Definition: StanTypedefs.h:14
TH1D * ToTH1(const Binning &bins) const
Definition: Hist.cxx:207
void Clear()
Definition: Spectrum.cxx:389
double fPOT
Definition: Spectrum.h:284
Float_t tmp
Definition: plot.C:36
void Fill(double x, double w=1)
Definition: Spectrum.cxx:321
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:277
Spectrum & operator-=(const Spectrum &rhs)
Definition: Spectrum.cxx:490
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:265
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:328
TH1 * ToTHX(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Definition: Spectrum.cxx:239
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const char * label
static Hist ZeroSparse(int nbins)
Definition: Hist.cxx:51
EBinType
Definition: UtilsExt.h:28
bool Initialized() const
Definition: Hist.h:31
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:563
const std::string & GetLabel1D() const
unsigned int seed
Definition: runWimpSim.h:102
bool AlmostEqual(double a, double b, double eps)
Definition: UtilsExt.cxx:40
Spectrum & operator/=(const Ratio &rhs)
Definition: Spectrum.cxx:519
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:377
Spectrum operator/(const Ratio &rhs) const
Definition: Spectrum.cxx:526
const Binning & GetBins1D() const
Appropriate binning and labelling for that 1D Var.
Int_t col[ntarg]
Definition: Style.C:29
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: UtilsExt.h:35
Spectrum AsimovData(double pot) const
Asimov data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:346
#define pot
Spectrum & operator=(const Spectrum &rhs)
Definition: Spectrum.cxx:75
Spectrum()
Constructor for Uninitialized()
Definition: Spectrum.h:267
double fLivetime
Definition: Spectrum.h:285
void Scale(double s)
Definition: Hist.cxx:296
void Reset()
Definition: Hist.cxx:382
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
unsigned int NDimensions() const
Definition: Spectrum.h:261
Represent the ratio between two spectra.
Definition: Ratio.h:8
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
void Fill(const Binning &bins, double x, double w)
Definition: Hist.cxx:265
const Binning bins
Definition: NumuCC_CPiBin.h:8
Spectrum operator+(const Spectrum &rhs) const
Definition: Spectrum.cxx:482
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:215
void Divide(const Hist &rhs)
Definition: Hist.cxx:531
Eigen::Array< stan::math::var, Eigen::Dynamic, 1 > ArrayXstan
Definition: StanUtils.h:7
TDirectory * dir
Definition: macro.C:5
Hist fHist
Definition: Ratio.h:56
double livetime
Definition: saveFDMCHists.C:21
int NBins() const
Definition: Binning.h:29
double GetBinContent(int i) const
Definition: Hist.cxx:349
static Hist Zero(int nbins)
Definition: Hist.cxx:41
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 T
Definition: Xdiff_gwt.C:5
virtual ~Spectrum()
Definition: Spectrum.cxx:106
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:262
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
double GetBinError(int i) const
Definition: Hist.cxx:241
Float_t w
Definition: plot.C:20
#define W(x)
Spectrum operator-(const Spectrum &rhs) const
Definition: Spectrum.cxx:496
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: UtilsExt.cxx:162
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
def sign(x)
Definition: canMan.py:197
Spectrum & operator*=(const Ratio &rhs)
Definition: Spectrum.cxx:504
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Spectrum.h:189
void Multiply(const Hist &rhs)
Definition: Hist.cxx:479
const Eigen::ArrayXstan & GetEigenStan() const
Definition: Hist.h:46
double Mean() const
Return mean of 1D histogram.
Definition: Spectrum.cxx:295
enum BeamMode string