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