Spectrum.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
2 
4 #include "CAFAna/Core/MultiVar.h"
6 #include "CAFAna/Core/HistogramTypePolicy.h" // gFloatMode
7 #include "CAFAna/Core/Ratio.h"
9 
10 #include "Utilities/func/MathUtil.h"
11 
12 #include "TDirectory.h"
13 #include "TH2.h"
14 #include "TH3.h"
15 #include "THnSparse.h"
16 #include "TObjString.h"
17 #include "TRandom3.h"
18 
19 #include <cassert>
20 #include <iostream>
21 
22 namespace ana
23 {
24  //----------------------------------------------------------------------
25  // The one constructor to rule them all. (Almost) everyone forwards here
26  Spectrum::Spectrum(const std::vector<std::string>& labels,
27  const std::vector<Binning>& bins,
28  ESparse sparse)
29  : fHistD(0), fHistF(0), fHistSparse(0), fPOT(0), fLivetime(0),
30  fLabels(labels), fBins(bins)
31  {
32  ConstructHistogram(sparse);
33  }
34 
35  //----------------------------------------------------------------------
36  // Or here...
38  const Binning& bins,
39  ESparse sparse)
40  : fHistD(0), fHistF(0), fHistSparse(0), fPOT(0), fLivetime(0),
41  fLabels(1, label), fBins(1, bins)
42  {
43  ConstructHistogram(sparse);
44  }
45 
46  //----------------------------------------------------------------------
49  const Var& var,
50  const Cut& cut,
51  const SystShifts& shift,
52  const Var& wei)
53  : Spectrum(label, bins)
54  {
55  loader.AddSpectrum(*this, var, cut, shift, wei);
56  }
57 
58  //----------------------------------------------------------------------
61  const MultiVar& var,
62  const Cut& cut,
63  const SystShifts& shift,
64  const Var& wei)
65  : Spectrum(label, bins)
66  {
67  loader.AddSpectrum(*this, var, cut, shift, wei);
68  }
69 
70  //----------------------------------------------------------------------
73  const NuTruthVar& var,
74  const NuTruthCut& cut,
75  const SystShifts& shifts,
76  const NuTruthVar& wei)
77  : Spectrum(label, bins)
78  {
79  loader.AddSpectrum(*this, var, cut, shifts, wei);
80  }
81 
82  //----------------------------------------------------------------------
84  const Binning& binsx, const NuTruthVar& varx,
85  const Binning& binsy, const NuTruthVar& vary,
86  const NuTruthCut& cut,
87  const SystShifts& shifts,
88  const NuTruthVar& wei)
89  : Spectrum({label, ""}, {binsx, binsy}) // TODO why no labelY?
90  {
91  loader.AddSpectrum(*this, Var2D(varx, binsx, vary, binsy), cut, shifts, wei);
92  }
93 
94  //----------------------------------------------------------------------
96  const HistAxis& axis,
97  const Cut& cut,
98  const SystShifts& shift,
99  const Var& wei)
100  : Spectrum(axis.GetLabels(), axis.GetBinnings())
101  {
102  loader.AddSpectrum(*this, axis.GetMultiDVar(), cut, shift, wei);
103  }
104 
105  //----------------------------------------------------------------------
107  const MultiVarHistAxis& axis,
108  const Cut& cut,
109  const SystShifts& shift,
110  const Var& wei)
111  : Spectrum(axis.GetLabels(), axis.GetBinnings())
112  {
113  loader.AddSpectrum(*this, axis.GetMultiDVar(), cut, shift, wei);
114  }
115 
116  //----------------------------------------------------------------------
118  const NuTruthHistAxis& axis,
119  const NuTruthCut& cut,
120  const SystShifts& shifts,
121  const NuTruthVar& weight)
122  : Spectrum(axis.GetLabels(), axis.GetBinnings())
123  {
124  loader.AddSpectrum(*this, axis.GetMultiDVar(), cut, shifts, weight);
125  }
126 
127  //----------------------------------------------------------------------
129  const Binning& bins)
130  : Spectrum(label, bins)
131  {
132  fPOT = pot;
134  }
135 
136  //----------------------------------------------------------------------
138  const std::vector<std::string>& labels,
139  const std::vector<Binning>& bins,
140  double pot, double livetime)
141  : fHistD(0), fHistF(0), fHistSparse(0), fPOT(pot), fLivetime(livetime), fLabels(labels), fBins(bins)
142  {
143  if(!h){
144  fHistD = 0;
145  fHistF = 0;
146  return;
147  }
148 
149  DontAddDirectory guard;
150 
151  const TString className = h->ClassName();
152 
153  if(className == "TH1D"){
154  // Shortcut if types match
155  fHistD = HistCache::Copy((TH1D*)h);
156  }
157  else if(className == "TH1F"){
158  fHistF = HistCache::Copy((TH1F*)h);
159  }
160  else{
161  fHistD = HistCache::New("", h->GetXaxis());
162  fHistD->Add(h);
163  }
164  }
165 
166  //----------------------------------------------------------------------
167  Spectrum::Spectrum(TH1* h, double pot, double livetime)
168  : Spectrum(h,
169  std::vector<std::string>(1, h->GetXaxis()->GetTitle()),
170  std::vector<Binning>(1, Binning::FromTAxis(h->GetXaxis())),
171  pot, livetime)
172  {
173  }
174 
175  //----------------------------------------------------------------------
176  Spectrum::Spectrum(std::unique_ptr<TH1D> h,
177  const std::vector<std::string>& labels,
178  const std::vector<Binning>& bins,
179  double pot, double livetime)
180  : fHistD(h.release()), fHistF(0), fHistSparse(0), fPOT(pot), fLivetime(livetime), fLabels(labels), fBins(bins)
181  {
182  }
183 
184  //----------------------------------------------------------------------
185  Spectrum::Spectrum(std::unique_ptr<TH1F> h,
186  const std::vector<std::string>& labels,
187  const std::vector<Binning>& bins,
188  double pot, double livetime)
189  : fHistD(0), fHistF(h.release()), fHistSparse(0), fPOT(pot), fLivetime(livetime), fLabels(labels), fBins(bins)
190  {
191  }
192 
193  //----------------------------------------------------------------------
195  const Binning& binsx, const Var& varx,
196  const Binning& binsy, const Var& vary,
197  const Cut& cut,
198  const SystShifts& shift,
199  const Var& wei)
200  : Spectrum(label, "", loader, binsx, varx, binsy, vary, cut, shift, wei)
201  {
202  // TODO do we want this variant when there's one with a labelY just below?
203  }
204 
205  //----------------------------------------------------------------------
207  const Binning& binsx, const MultiVar& varx,
208  const Binning& binsy, const MultiVar& vary,
209  const Cut& cut,
210  const SystShifts& shift,
211  const Var& wei)
212  : Spectrum(label, "", loader, binsx, varx, binsy, vary, cut, shift, wei)
213  {
214  }
215 
216  //----------------------------------------------------------------------
218  const HistAxis& xAxis,
219  const HistAxis& yAxis,
220  const Cut& cut,
221  const SystShifts& shift,
222  const Var& wei)
223  : Spectrum(xAxis.GetLabels()[0], yAxis.GetLabels()[0],
224  loader,
225  xAxis.GetBinnings()[0], xAxis.GetVars()[0],
226  yAxis.GetBinnings()[0], yAxis.GetVars()[0],
227  cut, shift, wei)
228  {
229  // TODO - do we want to keep this variant around?
230  assert(xAxis.NDimensions() == 1);
231  assert(yAxis.NDimensions() == 1);
232  }
233 
234  //----------------------------------------------------------------------
236  const std::string& yLabel,
238  const Binning& binsx, const Var& varx,
239  const Binning& binsy, const Var& vary,
240  const Cut& cut,
241  const SystShifts& shift,
242  const Var& wei)
243  : Spectrum({xLabel, yLabel}, {binsx, binsy})
244  {
245  Var multiDVar = Var2D(varx, binsx, vary, binsy);
246 
247  loader.AddSpectrum(*this, multiDVar, cut, shift, wei);
248  }
249 
250  //----------------------------------------------------------------------
252  const std::string& yLabel,
254  const Binning& binsx, const MultiVar& varx,
255  const Binning& binsy, const MultiVar& vary,
256  const Cut& cut,
257  const SystShifts& shift,
258  const Var& wei)
259  : Spectrum({xLabel, yLabel}, {binsx, binsy})
260  {
261  MultiVar multiDVar = MultiVar2D(varx, binsx, vary, binsy);
262 
263  loader.AddSpectrum(*this, multiDVar, cut, shift, wei);
264  }
265 
266  //----------------------------------------------------------------------
268  const Binning& binsx, const Var& varx,
269  const Binning& binsy, const Var& vary,
270  const Binning& binsz, const Var& varz,
271  const Cut& cut,
272  const SystShifts& shift,
273  const Var& wei,
274  ESparse sparse)
275  : Spectrum(label, "", "", loader, binsx, varx, binsy, vary, binsz, varz, cut, shift, wei, sparse)
276  {
277  // TODO do we want this variant when there's one with a labelY and labelZ
278  // just below?
279  }
280 
281  //----------------------------------------------------------------------
283  const std::string& yLabel,
284  const std::string& zLabel,
286  const Binning& binsx, const Var& varx,
287  const Binning& binsy, const Var& vary,
288  const Binning& binsz, const Var& varz,
289  const Cut& cut,
290  const SystShifts& shift,
291  const Var& wei,
292  ESparse sparse)
293  : Spectrum({xLabel, yLabel, zLabel}, {binsx, binsy, binsz}, sparse)
294  {
295  Var multiDVar = Var3D(varx, binsx, vary, binsy, varz, binsz);
296 
297  loader.AddSpectrum(*this, multiDVar, cut, shift, wei);
298  }
299 
300  //----------------------------------------------------------------------
302  const HistAxis& xAxis,
303  const HistAxis& yAxis,
304  const HistAxis& zAxis,
305  const Cut& cut,
306  const SystShifts& shift,
307  const Var& wei,
308  ESparse sparse)
309  : Spectrum(xAxis.GetLabels()[0], loader,
310  xAxis.GetBinnings()[0], xAxis.GetVars()[0],
311  yAxis.GetBinnings()[0], yAxis.GetVars()[0],
312  zAxis.GetBinnings()[0], zAxis.GetVars()[0],
313  cut, shift, wei, sparse)
314  {
315  // TODO - do we want to keep this variant around?
316  assert(xAxis.NDimensions() == 1);
317  assert(yAxis.NDimensions() == 1);
318  assert(zAxis.NDimensions() == 1);
319  }
320 
321  //----------------------------------------------------------------------
323  {
324  if(fHistD && fHistD->GetDirectory()){
325  static bool once = true;
326  if(once){
327  once = false;
328  std::cerr << "Spectrum's fHist (" << fHistD << ") is associated with a directory (" << fHistD->GetDirectory() << ". How did that happen?" << std::endl;
329  }
330  }
331 
333  { loader->RemoveSpectrum(this); }
334 
337 
338  delete fHistSparse;
339  }
340 
341  //----------------------------------------------------------------------
343  fHistD(0),
344  fHistF(0),
345  fHistSparse(0),
346  fPOT(rhs.fPOT),
347  fLivetime(rhs.fLivetime),
348  fLabels(rhs.fLabels),
349  fBins(rhs.fBins)
350  {
351  DontAddDirectory guard;
352 
353  assert(rhs.fHistD || rhs.fHistF || rhs.fHistSparse);
354  if(rhs.fHistD)
356  if(rhs.fHistF)
358  if(rhs.fHistSparse){
359  // Doesn't exist?
360  // fHistSparse = new THnSparseD(*rhs.fHistSparse);
361  fHistSparse = (THnSparseD*)rhs.fHistSparse->Clone();
362  }
363 
364  assert( rhs.fLoaderCount.empty() ); // Copying with pending loads is unexpected
365  }
366 
367  //----------------------------------------------------------------------
369  fHistD(0),
370  fHistF(0),
371  fHistSparse(0),
372  fPOT(rhs.fPOT),
373  fLivetime(rhs.fLivetime),
374  fLabels(rhs.fLabels),
375  fBins(rhs.fBins)
376  {
377  assert(rhs.fHistD || rhs.fHistF || rhs.fHistSparse);
378 
379  if(rhs.fHistD){
380  fHistD = rhs.fHistD;
381  rhs.fHistD = 0;
382  }
383  if(rhs.fHistF){
384  fHistF = rhs.fHistF;
385  rhs.fHistF = 0;
386  }
387  if(rhs.fHistSparse){
388  fHistSparse = rhs.fHistSparse;
389  rhs.fHistSparse = 0;
390  }
391 
392  for (auto & loader : rhs.fLoaderCount)
393  {
394  loader->ReplaceSpectrum(&rhs, this);
395  fLoaderCount.insert(loader);
396  }
397  rhs.fLoaderCount.clear();
398 
399  }
400 
401  //----------------------------------------------------------------------
403  : fHistD(nullptr), fHistF(nullptr), fHistSparse(nullptr),
404  fPOT(rhs.POT()), fLivetime(rhs.Livetime()),
405  fLabels(rhs.GetLabels()), fBins(rhs.GetBinnings())
406  {
407  // for now we don't have any other use cases besides 1D...
408  assert(rhs.GetBinnings().size() == 1);
409 
411 
412  // now fill in the bin content...
413  auto binC = rhs.ToBins(rhs.POT());
414  for (int binIdx = 0; binIdx < rhs.GetBinnings()[0].NBins()+2; binIdx++) // be sure to get under- and overflow
415  {
416  if (gFloatMode)
417  fHistF->SetBinContent(binIdx, binC[binIdx].val());
418  else
419  fHistD->SetBinContent(binIdx, binC[binIdx].val());
420  }
421  }
422 
423  //----------------------------------------------------------------------
425  {
426  if(this == &rhs) return *this;
427 
428  DontAddDirectory guard;
429 
432  delete fHistSparse;
433 
434  assert(rhs.fHistD || rhs.fHistF || rhs.fHistSparse);
435 
436  if(rhs.fHistD){
438  fHistF = 0;
439  fHistSparse = 0;
440  }
441 
442  if(rhs.fHistF){
444  fHistD = 0;
445  fHistSparse = 0;
446  }
447 
448  if(rhs.fHistSparse){
449  fHistSparse = (THnSparseD*)rhs.fHistSparse->Clone();
450  fHistD = 0;
451  fHistF = 0;
452  }
453 
454  fPOT = rhs.fPOT;
455  fLivetime = rhs.fLivetime;
456  fLabels = rhs.fLabels;
457  fBins = rhs.fBins;
458 
459  assert( fLoaderCount.empty() ); // Copying with pending loads is unexpected
460 
461  return *this;
462  }
463 
464  //----------------------------------------------------------------------
466  {
467  if(this == &rhs) return *this;
468 
471  delete fHistSparse;
472 
473  assert(rhs.fHistD || rhs.fHistF || rhs.fHistSparse);
474 
475  fHistD = rhs.fHistD;
476  fHistF = rhs.fHistF;
477  fHistSparse = rhs.fHistSparse;
478 
479  fPOT = rhs.fPOT;
480  fLivetime = rhs.fLivetime;
481  fLabels = rhs.fLabels;
482  fBins = rhs.fBins;
483 
484  rhs.fHistD = 0;
485  rhs.fHistF = 0;
486  rhs.fHistSparse = 0;
487 
488  for (auto & loader : rhs.fLoaderCount)
489  {
490  loader->ReplaceSpectrum(&rhs, this);
491  fLoaderCount.insert(loader);
492  }
493  rhs.fLoaderCount.clear();
494 
495  return *this;
496  }
497 
498  //----------------------------------------------------------------------
500  {
501  // use placement new to just overwrite this object using the relevant constructor
502  this->~Spectrum();
503  new (this) Spectrum(rhs);
504  return *this;
505  }
506 
507  //----------------------------------------------------------------------
509  {
510  DontAddDirectory guard;
511 
512  assert(!fHistD && !fHistF && !fHistSparse);
513 
514  Binning bins1D = fBins[0];
515  if(fBins.size() > 1){
516  int n = 1;
517  for(const Binning& b: fBins) n *= b.NBins();
518  bins1D = Binning::Simple(n, 0, n);
519  }
520 
521  if(sparse == kSparse){
522  assert(bins1D.IsSimple());
523  const int nbins = bins1D.NBins();
524  const double xmin = bins1D.Min();
525  const double xmax = bins1D.Max();
526  fHistSparse = new THnSparseD(UniqueName().c_str(), UniqueName().c_str(),
527  1, &nbins, &xmin, &xmax);
528 
529  // Ensure errors get accumulated properly
530  fHistSparse->Sumw2();
531  }
532  else if(sparse == kDense && !gFloatMode){
533  fHistD = HistCache::New("", bins1D);
534  fHistF = 0;
535 
536  // first clear the sum-of-weights.
537  // can be pre-filled with 0s if you happen to have made several Spectra
538  // out of the same Binning and a preceding one is still in the cache.
539  fHistD->Sumw2(false);
540  fHistD->Sumw2(); // Ensure errors get accumulated properly
541  }
542  else if(sparse == kFloat || gFloatMode){
543  fHistF = HistCache::NewF("", bins1D);
544  fHistD = 0;
545 
546  // see above
547  fHistF->Sumw2(false);
548  fHistF->Sumw2(); // Ensure errors get accumulated properly
549  }
550  }
551 
552  //----------------------------------------------------------------------
553  TH1D* Spectrum::ToTH1(double exposure, Color_t col, Style_t style,
554  EExposureType expotype,
555  EBinType bintype) const
556  {
557  // Could have a file temporarily open
558  DontAddDirectory guard;
559 
560  TH1D* ret = 0;
561  if(fHistD){
562  ret = HistCache::Copy(fHistD);
563  }
564  else if(fHistF){
565  ret = HistCache::New(fHistF->GetTitle(), fHistF->GetXaxis());
566  ret->Add(fHistF);
567  }
568  else{
569  ret = fHistSparse->Projection(0);
570  }
571 
572  if(expotype == kPOT){
573  const double pot = exposure;
574  if(fPOT){
575  ret->Scale(pot/fPOT);
576  }
577  else{
578  // Allow zero POT if there are also zero events
579  if(ret->Integral() > 0){
580  std::cout << "Error: Spectrum with " << ret->Integral()
581  << " entries has zero POT, no way to scale to "
582  << exposure << " POT.";
583  if(fLivetime > 0){
584  std::cout << " Spectrum has " << fLivetime << " seconds livetime. "
585  << "Did you mean to pass kLivetime to ToTH1()?";
586  }
587  std::cout << std::endl;
588  abort();
589  }
590  }
591  }
592  if(expotype == kLivetime){
593  const double livetime = exposure;
594  if(fLivetime){
595  ret->Scale(livetime/fLivetime);
596  }
597  else{
598  // Allow zero exposure if there are also zero events
599  if(ret->Integral() > 0){
600  std::cout << "Error: Spectrum with " << ret->Integral()
601  << " entries has zero livetime, no way to scale to "
602  << livetime << " seconds.";
603  if(fPOT > 0){
604  std::cout << " Spectrum has " << fPOT << " POT. "
605  << "Did you mean to pass kPOT to ToTH1()?";
606  }
607  std::cout << std::endl;
608  abort();
609  }
610  }
611  }
612 
614  for(const std::string& l: fLabels) label += l + " and ";
615  label.resize(label.size()-5); // drop the last "and"
616  ret->GetXaxis()->SetTitle(label.c_str());
617 
618  if(fBins.size() == 1){
619  for(unsigned int i = 0; i < fBins[0].Labels().size(); ++i){
620  ret->GetXaxis()->SetBinLabel(i+1, fBins[0].Labels()[i].c_str());
621  }
622  }
623 
624  ret->GetYaxis()->SetTitle("Events");
625 
626 
627  ret->SetLineColor(col);
628  ret->SetMarkerColor(col);
629  ret->SetLineStyle(style);
630 
631  if(bintype == kBinDensity) ret->Scale(1, "width");
632 
633  // Allow GetMean() and friends to work even if this histogram never had any
634  // explicit Fill() calls made.
635  if(ret->GetEntries() == 0) ret->SetEntries(1);
636 
637  return ret;
638  }
639 
640  //----------------------------------------------------------------------
641  TH2* Spectrum::ToTH2(double exposure, EExposureType expotype, EBinType bintype) const
642  {
643  if(fBins.size() != 2){
644  std::cout << "Error: This Spectrum does not appear to be 2D." << std::endl;
645  abort();
646  }
647 
648  TH2* ret = ana::ToTH2(*this, exposure, expotype, fBins[0], fBins[1]);
649 
650  for(int iax = 0; iax < 2; ++iax){
651  TAxis* ax = ret->GetXaxis();
652  if(iax == 1) ax = ret->GetYaxis();
653 
654  ax->SetTitle(fLabels[iax].c_str());
655 
656  for(unsigned int i = 0; i < fBins[iax].Labels().size(); ++i){
657  ax->SetBinLabel(i+1, fBins[iax].Labels()[i].c_str());
658  }
659  }
660 
661  if(bintype == kBinDensity) ret->Scale(1, "width");
662 
663  // Allow GetMean() and friends to work even if this histogram never had any
664  // explicit Fill() calls made.
665  if(ret->GetEntries() == 0) ret->SetEntries(1);
666 
667  return ret;
668  }
669 
670  //----------------------------------------------------------------------
671  TH2* Spectrum::ToTH2NormX(double exposure, EExposureType expotype) const
672  {
673  TH2* xyhist = ToTH2(exposure, expotype);
674  if(!xyhist) return nullptr;
675 
676  const int nbinsx = fBins[0].NBins();
677  const int nbinsy = fBins[1].NBins();
678 
679  // Normalize 2D histogram to X-axis spectrum
680  for(int i=1; i<=nbinsx; ++i){
681  double norm = 0.0;
682  for(int j=1; j<=nbinsy; ++j){
683  norm += xyhist->GetBinContent(i, j);
684  }
685  /// If no entries in the column, skip normalization
686  if(norm < 0.0000001) continue;
687 
688  norm = 1.0 / norm;
689  for(int j=1; j<=nbinsy; ++j){
690  xyhist->SetBinContent(i,j, xyhist->GetBinContent(i, j) * norm);
691  }
692  }
693 
694  // Allow GetMean() and friends to work even if this histogram never had any
695  // explicit Fill() calls made.
696  if(xyhist->GetEntries() == 0) xyhist->SetEntries(1);
697 
698  return xyhist;
699  }
700 
701  //----------------------------------------------------------------------
702  TH3* Spectrum::ToTH3(double exposure, EExposureType expotype, EBinType bintype) const
703  {
704  if(fBins.size() != 3){
705  std::cout << "Error: This Spectrum does not appear to be 3D." << std::endl;
706  abort();
707  }
708 
709  TH3* ret = ana::ToTH3(*this, exposure, expotype,
710  fBins[0], fBins[1], fBins[2]);
711 
712  for(int iax = 0; iax < 3; ++iax){
713  TAxis* ax = ret->GetXaxis();
714  if(iax == 1) ax = ret->GetYaxis();
715  if(iax == 2) ax = ret->GetZaxis();
716 
717  ax->SetTitle(fLabels[iax].c_str());
718 
719  for(unsigned int i = 0; i < fBins[iax].Labels().size(); ++i){
720  ax->SetBinLabel(i+1, fBins[iax].Labels()[i].c_str());
721  }
722  }
723 
724  if(bintype == kBinDensity) ret->Scale(1, "width");
725 
726  // Allow GetMean() and friends to work even if this histogram never had any
727  // explicit Fill() calls made.
728  if(ret->GetEntries() == 0) ret->SetEntries(1);
729 
730  return ret;
731  }
732 
733  //----------------------------------------------------------------------
734  void Spectrum::Scale(double c)
735  {
736  if(fHistD) fHistD->Scale(c);
737  if(fHistF) fHistF->Scale(c);
738  if(fHistSparse) fHistSparse->Scale(c);
739  }
740 
741  //----------------------------------------------------------------------
742  double Spectrum::Integral(double exposure, double* err,
743  EExposureType expotype) const
744  {
745  const double ratio = (expotype == kPOT) ? exposure/fPOT : exposure/fLivetime;
746 
747  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
748 
749  if(err){
750  *err = 0;
751 
752  for(int i = 0; i < h->GetNbinsX()+2; ++i){
753  *err += util::sqr(h->GetBinError(i));
754  }
755  *err = sqrt(*err) * ratio;
756  }
757 
758  // TODO how to integrate fHistSparse?
759 
760  return h->Integral(0, -1) * ratio;
761  }
762 
763  //----------------------------------------------------------------------
764  double Spectrum::Mean() const
765  {
766  // Allow GetMean() to work even if this histogram never had any explicit
767  // Fill() calls made.
768  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
769  if(h->GetEntries() == 0) h->SetEntries(1);
770  return h->GetMean();
771  }
772 
773  //----------------------------------------------------------------------
774  void Spectrum::Fill(double x, double w)
775  {
776  assert( (fHistD || fHistF || fHistSparse) && "Somehow both fHist and fHistSparse are null in Spectrum::Fill" );
777 
778  if(fHistD)
779  fHistD->Fill(x, w);
780  else if(fHistF)
781  fHistF->Fill(x, w);
782  else if(fHistSparse)
783  fHistSparse->Fill(&x, w);
784  }
785 
786  //----------------------------------------------------------------------
787  Spectrum Spectrum::MockData(double pot, int idx) const
788  {
789  Spectrum ret = FakeData(pot);
790 
791  TRandom3 rnd(idx); // zero seeds randomly
792 
793  if(ret.fHistD || ret.fHistF){
794  TH1* h = ret.fHistD ? (TH1*)ret.fHistD : (TH1*)ret.fHistF;
795  for(int i = 0; i < h->GetNbinsX()+2; ++i){
796  h->SetBinContent(i, rnd.Poisson(h->GetBinContent(i)));
797  }
798  }
799  if(ret.fHistSparse){
800  for(int i = 0; i < ret.fHistSparse->GetNbins(); ++i)
801  ret.fHistSparse->SetBinContent(i, rnd.Poisson(ret.fHistSparse->GetBinContent(i)));
802  }
803 
804  // Drop old errors, which are based on the MC statistics, and create new
805  // ones that are based on the prediction for the data
806  if(ret.fHistD){
807  ret.fHistD->Sumw2(false);
808  ret.fHistD->Sumw2();
809  }
810  if(ret.fHistF){
811  ret.fHistF->Sumw2(false);
812  ret.fHistF->Sumw2();
813  }
814 
815  return ret;
816  }
817 
818  //----------------------------------------------------------------------
820  {
821  Spectrum ret = *this;
822  if(fPOT > 0){
823  if(ret.fHistD)
824  ret.fHistD->Scale(pot/fPOT);
825  if(ret.fHistF)
826  ret.fHistF->Scale(pot/fPOT);
827  if(ret.fHistSparse)
828  ret.fHistSparse->Scale(pot/fPOT);
829  }
830  ret.fPOT = pot;
831 
832  // Drop old errors, which are based on the MC statistics, and create new
833  // ones that are based on the prediction for the data
834  if(ret.fHistD){
835  ret.fHistD->Sumw2(false);
836  ret.fHistD->Sumw2();
837  }
838  if(ret.fHistF){
839  ret.fHistF->Sumw2(false);
840  ret.fHistF->Sumw2();
841  }
842 
843  return ret;
844  }
845 
846  //----------------------------------------------------------------------
847  Spectrum Spectrum::FakeData(double pot, double livetime) const
848  {
849  Spectrum ret = *this;
850 
851  if(fPOT > 0){
852  if(ret.fHistD)
853  ret.fHistD->Scale(pot/fPOT);
854  if(ret.fHistF)
855  ret.fHistF->Scale(pot/fPOT);
856  if(ret.fHistSparse)
857  ret.fHistSparse->Scale(pot/fPOT);
858  }
859  ret.fPOT = pot;
860 
861  // overwrite livetime for fake data
862  // PlusEqualsHelper will take over the rest
863  ret.fLivetime = livetime;
864 
865  if(ret.fHistD){
866  ret.fHistD->Sumw2(false);
867  ret.fHistD->Sumw2();
868  }
869  if(ret.fHistF){
870  ret.fHistF->Sumw2(false);
871  ret.fHistF->Sumw2();
872  }
873 
874  return ret;
875  }
876 
877  //----------------------------------------------------------------------
879  {
880  if(fHistD) fHistD->Reset();
881  if(fHistF) fHistF->Reset();
882  if(fHistSparse) fHistSparse->Reset();
883  }
884 
885  //----------------------------------------------------------------------
887  { fLoaderCount.erase(p); }
888 
889  //----------------------------------------------------------------------
891  { fLoaderCount.insert(p); }
892 
893  //----------------------------------------------------------------------
895  {
896  // In this case it would be OK to have no POT/livetime
897  if(rhs.fHistD && rhs.fHistD->Integral(0, -1) == 0) return *this;
898  if(rhs.fHistF && rhs.fHistF->Integral(0, -1) == 0) return *this;
899 
900  if((!fPOT && !fLivetime) || (!rhs.fPOT && !rhs.fLivetime)){
901  std::cout << "Error: can't sum Spectrum with no POT or livetime."
902  << std::endl;
903  abort();
904  }
905 
906  if(!fLivetime && !rhs.fPOT){
907  std::cout << "Error: can't sum Spectrum with POT ("
908  << fPOT << ") but no livetime and Spectrum with livetime ("
909  << rhs.fLivetime << " sec) but no POT." << std::endl;
910  abort();
911  }
912 
913  if(!fPOT && !rhs.fLivetime){
914  std::cout << "Error: can't sum Spectrum with livetime ("
915  << fLivetime << " sec) but no POT and Spectrum with POT ("
916  << rhs.fPOT << ") but no livetime." << std::endl;
917  abort();
918  }
919 
920  // And now there are still a bunch of good cases to consider
921 
922  if(fPOT && rhs.fPOT){
923  // Scale by POT when possible
924  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
925  TH1* rh = rhs.fHistD ? (TH1*)rhs.fHistD : (TH1*)rhs.fHistF;
926  if(rh) h->Add(rh, sign*fPOT/rhs.fPOT);
927  if(rhs.fHistSparse) fHistSparse->Add(rhs.fHistSparse, sign*fPOT/rhs.fPOT);
928 
929  if(fLivetime && rhs.fLivetime){
930  // If POT/livetime ratios match, keep regular lifetime, otherwise zero
931  // it out.
932  if(AlmostEqual(fLivetime*rhs.fPOT, rhs.fLivetime*fPOT))
933  fLivetime = 0;
934  }
935  if(!fLivetime && rhs.fLivetime){
936  // If the RHS has a livetime and we don't, copy it in (suitably scaled)
937  fLivetime = rhs.fLivetime * fPOT/rhs.fPOT;
938  }
939  // Otherwise, keep our own livetime (if any)
940 
941  return *this;
942  }
943 
944  if(fLivetime && rhs.fLivetime){
945  // Scale by livetime, the only thing in common
946  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
947  TH1* rh = rhs.fHistD ? (TH1*)rhs.fHistD : (TH1*)rhs.fHistF;
948  if(rh) h->Add(rh, sign*fLivetime/rhs.fLivetime);
949  if(rhs.fHistSparse) fHistSparse->Add(rhs.fHistSparse, sign*fLivetime/rhs.fLivetime);
950 
951  if(!fPOT && rhs.fPOT){
952  // If the RHS has a POT and we don't, copy it in (suitably scaled)
953  fPOT = rhs.fPOT * fLivetime/rhs.fLivetime;
954  }
955  // Otherwise, keep our own POT (if any)
956 
957  return *this;
958  }
959 
960  // That should have been all the cases. I definitely want to know what
961  // happened if it wasn't.
962  std::cout << "Spectrum::operator+=(). How did we get here? "
963  << fPOT << " " << fLivetime << " "
964  << rhs.fPOT << " " << rhs.fLivetime << std::endl;
965  abort();
966  }
967 
968  //----------------------------------------------------------------------
970  {
971  return PlusEqualsHelper(rhs, +1);
972  }
973 
974  //----------------------------------------------------------------------
976  {
977  Spectrum ret = *this;
978  ret += rhs;
979  return ret;
980  }
981 
982  //----------------------------------------------------------------------
984  {
985  return PlusEqualsHelper(rhs, -1);
986  }
987 
988  //----------------------------------------------------------------------
990  {
991  Spectrum ret = *this;
992  ret -= rhs;
993  return ret;
994  }
995 
996  //----------------------------------------------------------------------
998  {
999  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
1000  h->Multiply(rhs.fHist);
1001  return *this;
1002  }
1003 
1004  //----------------------------------------------------------------------
1006  {
1007  Spectrum ret = *this;
1008  ret *= rhs;
1009  return ret;
1010  }
1011 
1012  //----------------------------------------------------------------------
1014  {
1015  TH1* h = fHistD ? (TH1*)fHistD : (TH1*)fHistF;
1016  h->Divide(rhs.fHist);
1017  return *this;
1018  }
1019 
1020  //----------------------------------------------------------------------
1022  {
1023  Spectrum ret = *this;
1024  ret /= rhs;
1025  return ret;
1026  }
1027 
1028  //----------------------------------------------------------------------
1029  void Spectrum::SaveTo(TDirectory* dir) const
1030  {
1031  TDirectory* tmp = gDirectory;
1032  dir->cd();
1033 
1034  TObjString("Spectrum").Write("type");
1035 
1036  if(fHistD) fHistD->Write("hist");
1037  if(fHistF) fHistF->Write("hist");
1038  if(fHistSparse) fHistSparse->Write("hist_sparse");
1039  TH1D hPot("", "", 1, 0, 1);
1040  hPot.Fill(.5, fPOT);
1041  hPot.Write("pot");
1042  TH1D hLivetime("", "", 1, 0, 1);
1043  hLivetime.Fill(.5, fLivetime);
1044  hLivetime.Write("livetime");
1045 
1046  for(unsigned int i = 0; i < fBins.size(); ++i){
1047  TObjString(fLabels[i].c_str()).Write(TString::Format("label%d", i).Data());
1048  fBins[i].SaveTo(dir->mkdir(TString::Format("bins%d", i)));
1049  }
1050 
1051  tmp->cd();
1052  }
1053 
1054  //----------------------------------------------------------------------
1055  std::unique_ptr<Spectrum> Spectrum::LoadFrom(TDirectory* dir)
1056  {
1057  DontAddDirectory guard;
1058 
1059  TObjString* tag = (TObjString*)dir->Get("type");
1060  assert(tag);
1061  assert(tag->GetString() == "Spectrum");
1062  delete tag;
1063 
1064  TH1D* spectD = 0;
1065  TH1F* spectF = 0;
1066  dir->GetObject("hist", spectD);
1067  dir->GetObject("hist", spectF);
1068  THnSparseD* spectSparse = (THnSparseD*)dir->Get("hist_sparse");
1069  assert(spectD || spectF || spectSparse);
1070  TH1* hPot = (TH1*)dir->Get("pot");
1071  assert(hPot);
1072  TH1* hLivetime = (TH1*)dir->Get("livetime");
1073  assert(hLivetime);
1074 
1075  std::vector<std::string> labels;
1076  std::vector<Binning> bins;
1077  for(int i = 0; ; ++i){
1078  TDirectory* subdir = dir->GetDirectory(TString::Format("bins%d", i));
1079  if(!subdir) break;
1080  bins.push_back(*Binning::LoadFrom(subdir));
1081  TObjString* label = (TObjString*)dir->Get(TString::Format("label%d", i));
1082  labels.push_back(label ? label->GetString().Data() : "");
1083  delete subdir;
1084  delete label;
1085  }
1086 
1087  if(bins.empty() && labels.empty()){
1088  // Must be an old file. Make an attempt at backwards compatibility.
1089  if(spectD){
1090  bins.push_back(Binning::FromTAxis(spectD->GetXaxis()));
1091  labels.push_back(spectD->GetXaxis()->GetTitle());
1092  }
1093  else{
1094  bins.push_back(Binning::FromTAxis(spectSparse->GetAxis(0)));
1095  labels.push_back(spectSparse->GetAxis(0)->GetTitle());
1096  }
1097  }
1098 
1099  if(spectSparse){
1100  std::unique_ptr<Spectrum> ret = std::make_unique<Spectrum>((TH1*)0, labels, bins, hPot->GetBinContent(1), hLivetime->GetBinContent(1));
1101  ret->fHistSparse = spectSparse;
1102 
1103  delete spectD;
1104  delete spectF;
1105  delete hPot;
1106  delete hLivetime;
1107  return ret;
1108  }
1109  else{
1110  std::unique_ptr<Spectrum> ret;
1111  if(spectD)
1112  ret = std::make_unique<Spectrum>(std::unique_ptr<TH1D>(spectD), labels, bins, hPot->GetBinContent(1), hLivetime->GetBinContent(1));
1113  else
1114  ret = std::make_unique<Spectrum>(std::unique_ptr<TH1F>(spectF), labels, bins, hPot->GetBinContent(1), hLivetime->GetBinContent(1));
1115  delete spectSparse;
1116  delete hPot;
1117  delete hLivetime;
1118  return ret;
1119  }
1120  }
1121 }
TH1D * fHistD
Definition: Spectrum.h:314
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:641
Spectrum operator*(const Ratio &rhs) const
Definition: Spectrum.cxx:1005
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Oscillation analysis framework, runs over CAF files outside of ART.
TH2 * rh
Definition: drawXsec.C:5
unsigned int NDimensions() const
Definition: HistAxis.h:75
static TH1D * Copy(const TH1D *h)
Definition: HistCache.cxx:175
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:553
const SpillVar fPOT([](const caf::SRSpillProxy *spill){return spill->spillpot;})
const Var weight
subdir
Definition: cvnie.py:7
std::vector< stan::math::var > ToBins(double pot) const
virtual void AddSpectrum(Spectrum &spect, const Var &var, const Cut &cut, const SystShifts &shift, const Var &wei=kUnweighted)
For use by the Spectrum constructor.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Spectrum & PlusEqualsHelper(const Spectrum &rhs, int sign)
Helper for operator+= and operator-=.
Definition: Spectrum.cxx:894
static Binning FromTAxis(const TAxis *ax)
Definition: Binning.cxx:122
const char * p
Definition: xmltok.h:285
Spectrum & operator+=(const Spectrum &rhs)
Definition: Spectrum.cxx:969
Divide bin contents by bin widths.
Definition: Utilities.h:49
TH1F * fHistF
Definition: Spectrum.h:315
OStream cerr
Definition: OStream.cxx:7
var
TH2 * ToTH2NormX(double exposure, EExposureType expotype=kPOT) const
Spectrum must be 2D to obtain TH2. Normalized to X axis.
Definition: Spectrum.cxx:671
void AddLoader(SpectrumLoaderBase *)
Definition: Spectrum.cxx:890
void ConstructHistogram(ESparse sparse=kDense)
Definition: Spectrum.cxx:508
void Clear()
Definition: Spectrum.cxx:878
double fPOT
Definition: Spectrum.h:317
Float_t tmp
Definition: plot.C:36
void Fill(double x, double w=1)
Definition: Spectrum.cxx:774
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
Spectrum & operator-=(const Spectrum &rhs)
Definition: Spectrum.cxx:983
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:734
MultiVar MultiVar2D(const MultiVar &a, const Binning &binsa, const MultiVar &b, const Binning &binsb)
Definition: MultiVar.cxx:57
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const char * label
bool IsSimple() const
Definition: Binning.h:29
std::vector< std::string > fLabels
Definition: Spectrum.h:323
EBinType
Definition: Utilities.h:46
const int nbins
Definition: cellShifts.C:15
std::set< SpectrumLoaderBase * > fLoaderCount
This count is maintained by SpectrumLoader, as a sanity check.
Definition: Spectrum.h:321
double Min() const
Definition: Binning.h:26
THnSparseT< TArrayD > THnSparseD
Definition: Spectrum.h:24
Spectrum & operator/=(const Ratio &rhs)
Definition: Spectrum.cxx:1013
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
Spectrum operator/(const Ratio &rhs) const
Definition: Spectrum.cxx:1021
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
Int_t col[ntarg]
Definition: Style.C:29
EExposureType
For use as an argument to Spectrum::ToTH1.
Definition: Utilities.h:54
correl_xv GetXaxis() -> SetDecimals()
TH1D * fHist
Definition: Ratio.h:37
#define pot
Spectrum & operator=(const Spectrum &rhs)
Definition: Spectrum.cxx:424
double fLivetime
Definition: Spectrum.h:318
const double j
Definition: BetheBloch.cxx:29
bool AlmostEqual(double a, double b)
Definition: Utilities.cxx:732
loader
Definition: demo0.py:10
T GetMultiDVar() const
Definition: HistAxis.h:83
void RemoveLoader(SpectrumLoaderBase *)
Definition: Spectrum.cxx:886
const HistogramTypePolicy gFloatMode
double POT() const
Definition: Spectrum.h:263
THnSparseD * fHistSparse
Definition: Spectrum.h:316
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:374
Represent the ratio between two spectra.
Definition: Ratio.h:8
def sign(x)
Definition: canMan.py:204
static TH1F * NewF(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:46
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
double Max() const
Definition: Binning.h:27
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:787
const std::vector< VarDef > GetVars(TString opt, TString type="all")
Spectrum operator+(const Spectrum &rhs) const
Definition: Spectrum.cxx:975
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:702
Float_t norm
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
std::vector< Binning > fBins
Definition: Spectrum.h:324
int NBins() const
Definition: Binning.h:25
const hit & b
Definition: hits.cxx:21
static TH1D * New(const std::string &title, const Binning &bins)
Definition: HistCache.cxx:24
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
virtual ~Spectrum()
Definition: Spectrum.cxx:322
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1055
GenericVar< T > Var2D(const GenericVar< T > &a, const Binning &binsa, const GenericVar< T > &b, const Binning &binsb)
Variable formed from two input variables.
Definition: Var.cxx:245
double POT() const
Definition: SpectrumStan.h:47
Prevent histograms being added to the current directory.
Definition: Utilities.h:66
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Definition: Utilities.cxx:424
Float_t w
Definition: plot.C:20
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:221
Spectrum operator-(const Spectrum &rhs) const
Definition: Spectrum.cxx:989
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
GenericVar< T > Var3D(const GenericVar< T > &a, const Binning &binsa, const GenericVar< T > &b, const Binning &binsb, const GenericVar< T > &c, const Binning &binsc)
This is just like a Var2D, but useful for 3D Spectra.
Definition: Var.cxx:271
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299
Spectrum & operator*=(const Ratio &rhs)
Definition: Spectrum.cxx:997
h
Definition: demo3.py:41
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029
double Mean() const
Return mean of 1D histogram.
Definition: Spectrum.cxx:764
Spectrum(const std::string &label, const Binning &bins, SpectrumLoaderBase &loader, const Var &var, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
Definition: Spectrum.cxx:47
const std::vector< Binning > & GetBinnings() const
Definition: SpectrumStan.h:43