Multiverse.cxx
Go to the documentation of this file.
2 
3 #include <algorithm>
4 #include <iostream>
5 
7 
8 #include "TObjString.h"
9 #include "TDirectory.h"
10 #include "Math/ProbFunc.h"
11 #include "TH1.h"
12 #include "TH2.h"
13 #include "TH3.h"
14 #include "TVectorD.h"
15 
16 namespace ana {
17  ///////////////////////////////////////////////////////////////////////
18  void
20  {
21  if(fSpectra.size() > 0 && fHists.size() == 0) InitFromSpectra();
22  for(auto i = 0u; i < fHists.size(); i++) {
23  fHists[i]->Scale(scale);
24  }
25  }
26 
27  ///////////////////////////////////////////////////////////////////////
30  const std::function<TransformFunc_t>& transform)
31  {
32  Multiverse ret(rhs);
33  ret.Transform(transform);
34  return ret;
35  }
36 
37  ///////////////////////////////////////////////////////////////////////
38  void
39  Multiverse::Transform(const std::function<TransformFunc_t>& transform)
40  {
41  // if rhs also has spectra and we want to transform those,
42  // do that
43  double POT = 0;
44  double livetime = 0;
45  std::vector<std::string> labels;
46  std::vector<Binning> binning;
47  if(fSpectra.size() != 0) {
48  POT = fSpectra[0]->POT();
49  livetime = fSpectra[0]->Livetime();
50  labels = fSpectra[0]->GetLabels();
51  binning = fSpectra[0]->GetBinnings();
52  if(fHists.size() == 0) InitFromSpectra();
53  }
54 
55  std::vector<TH1*> xformed;
56  std::vector<Spectrum*> specs;
57  for(auto i = 0u; i < fHists.size(); i++) {
58  xformed.push_back(transform(fHists[i]));
59  if(fSpectra.size() != 0) {
60  specs.push_back(ToSpectrum(xformed.back(), POT, livetime));
61  delete fSpectra[i];
62  }
63  delete fHists[i];
64  }
65 
66  fHists = xformed;
67  fSpectra = specs;
68  }
69 
70  ///////////////////////////////////////////////////////////////////////
71  void
73  const std::function<ProductFunc_t>& transform)
74  {
75  // if rhs also has spectra and we want to transform those,
76  // do that
77  double POT = 0;
78  double livetime = 0;
79  std::vector<std::string> labels;
80  std::vector<Binning> binning;
81  if(fSpectra.size() != 0) {
82  POT = rhs.fSpectra[0]->POT();
83  livetime = rhs.fSpectra[0]->Livetime();
84  labels = rhs.fSpectra[0]->GetLabels();
85  binning = rhs.fSpectra[0]->GetBinnings();
86  if(fHists.size() == 0) InitFromSpectra();
87  }
88 
89  std::vector<TH1*> xformed;
90  std::vector<Spectrum*> specs;
91  for(auto i = 0u; i < fHists.size(); i++) {
92  xformed.push_back(transform(fHists[i], rhs.fHists[i]));
93  if(fSpectra.size() != 0) {
94  specs.push_back(ToSpectrum(xformed.back(), POT, livetime));
95  delete fSpectra[i];
96  }
97  delete fHists[i];
98 
99  }
100 
101  fHists = xformed;
102  fSpectra = specs;
103  }
104 
105  ///////////////////////////////////////////////////////////////////////
106  const std::vector<TH1*>
107  Multiverse::GetUniversesHist(int line_color, int line_style)
108  {
109  if(fHists.size() == 0)
110  InitFromSpectra();
111  for(auto ihist = 0u; ihist < fHists.size(); ihist++) {
112  fHists[ihist]->SetLineColor(line_color);
113  fHists[ihist]->SetLineStyle(line_style);
114  }
115  return this->GetUniversesHist();
116  }
117 
118  ///////////////////////////////////////////////////////////////////////
119  const std::vector<TH1*>
121  {
122  if(fHists.size() == 0)
123  InitFromSpectra();
124  return fHists;
125  }
126 
127  ///////////////////////////////////////////////////////////////////////
128  void
129  Multiverse::DivisionHelper(TH1 * h1, TH1 * h2, bool binomial_errors)
130  {
131  if(binomial_errors)
132  h1->Divide(h1, h2, 1, 1, "B");
133  else
134  h1->Divide(h2);
135  }
136 
137  ///////////////////////////////////////////////////////////////////////
138  void
139  Multiverse::Divide(Multiverse & rhs, bool binomial_errors)
140  {
141  if(fHists.size() == 0)
142  this->InitFromSpectra();
143  if(rhs.fHists.size() == 0)
144  rhs.InitFromSpectra();
145  assert(fHists.size() == rhs.fHists.size() &&
146  "Multiverses are not compatible");
147  for(int i = 0; i < fNUniverses; i++) {
148  DivisionHelper(fHists[i], rhs.fHists[i], binomial_errors);
149  }
150  }
151 
152  ///////////////////////////////////////////////////////////////////////
153  Multiverse &
155  {
156  assert(fHists.size() == rhs.fHists.size() &&
157  fHists.size() != 0 &&
158  "Multiverses are not compatible");
159  for(int i = 0; i < fNUniverses; i++) {
160  fHists[i]->Multiply(rhs.fHists[i]);
161  }
162  return *this;
163  }
164 
165  ///////////////////////////////////////////////////////////////////////
166  Multiverse
168  {
169  Multiverse ret = *this;
170  assert(ret.fHists.size() == rhs.fHists.size() &&
171  ret.fHists.size() != 0 &&
172  "Multiverses are not compatible");
173  ret *= rhs;
174  return ret;
175  }
176 
177  ///////////////////////////////////////////////////////////////////////
178  Multiverse &
180  {
181  assert(fHists.size() != 0 &&
182  "Multiverse must have histograms initialized");
183  TH1 * h = ToHist(&rhs);
184  for(int i = 0; i < fNUniverses; i++) {
185  fHists[i]->Divide(h);
186  }
187  delete h;
188  return *this;
189  }
190 
191  ///////////////////////////////////////////////////////////////////////
192  Multiverse
193  Multiverse::operator/(const Spectrum & rhs) const
194  {
195  Multiverse ret = *this;
196  assert(fHists.size() != 0 &&
197  "Multiverse must have histograms initialized");
198  ret /= rhs;
199  return ret;
200  }
201 
202  ///////////////////////////////////////////////////////////////////////
203  Multiverse &
205  {
206  assert(fHists.size() != 0 &&
207  "Multiverse must have histograms initialized");
208  TH1 * h = ToHist(&rhs);
209  for(int i = 0; i < fNUniverses; i++) {
210  fHists[i]->Multiply(h);
211  }
212  delete h;
213  return *this;
214  }
215 
216  ///////////////////////////////////////////////////////////////////////
217  Multiverse
218  Multiverse::operator*(const Spectrum & rhs) const
219  {
220  Multiverse ret = *this;
221  assert(fHists.size() != 0 &&
222  "Multiverse must have histograms initialized");
223  ret *= rhs;
224  return ret;
225  }
226 
227  ///////////////////////////////////////////////////////////////////////
228  Multiverse &
229  Multiverse::operator/=(const TH1 & rhs)
230  {
231  assert(fHists.size() != 0 &&
232  "Multiverse must have histograms initialized");
233  for(int i = 0; i < fNUniverses; i++) {
234  fHists[i]->Divide(&rhs);
235  }
236  return *this;
237  }
238 
239  ///////////////////////////////////////////////////////////////////////
240  Multiverse
241  Multiverse::operator/(const TH1 & rhs) const
242  {
243  Multiverse ret = *this;
244  assert(fHists.size() != 0 &&
245  "Multiverse must have histograms initialized");
246  ret /= rhs;
247  return ret;
248  }
249 
250  ///////////////////////////////////////////////////////////////////////
251  Multiverse &
252  Multiverse::operator*=(const TH1 & rhs)
253  {
254  assert(fHists.size() != 0 &&
255  "Multiverse must have histograms initialized");
256  for(int i = 0; i < fNUniverses; i++) {
257  fHists[i]->Multiply(&rhs);
258  }
259  return *this;
260  }
261 
262  ///////////////////////////////////////////////////////////////////////
263  Multiverse
264  Multiverse::operator*(const TH1 & rhs) const
265  {
266  Multiverse ret = *this;
267  assert(fHists.size() != 0 &&
268  "Multiverse must have histograms initialized");
269  ret *= rhs;
270  return ret;
271  }
272 
273  ///////////////////////////////////////////////////////////////////////
275  const NuTruthHistAxis & histaxis,
276  const NuTruthCut & cut,
277  std::vector<SystShifts> shifts,
278  const NuTruthVar & weight)
279  {
280  fNUniverses = shifts.size();
281  for(int i = 0; i < (int) shifts.size(); i++) {
282  fSpectra.push_back(new Spectrum(loader,
283  histaxis,
284  cut,
285  shifts[i],
286  weight));
287  }
288  }
289 
290  ///////////////////////////////////////////////////////////////////////
291  Multiverse::Multiverse(std::vector<SpectrumLoaderBase*> loaders,
292  const NuTruthHistAxis & histaxis,
293  const NuTruthCut & cut,
294  std::vector<SystShifts> shifts,
295  const NuTruthVar & weight)
296  {
297  fNUniverses = shifts.size();
298  for(int i = 0; i < (int) shifts.size(); i++) {
299  fSpectra.push_back(new Spectrum(*loaders[0],
300  histaxis,
301  cut,
302  shifts[i],
303  weight));
304  for(auto iload = 1u; iload < loaders.size(); iload++) {
305  loaders[iload]->AddSpectrum(*fSpectra.back(),
306  histaxis.GetVar1D(),
307  cut,
308  shifts[i],
309  weight);
310  }
311  }
312  }
313 
314  ///////////////////////////////////////////////////////////////////////
316  const HistAxis & histaxis,
317  const Cut & cut,
318  std::vector<SystShifts> shifts,
319  const Var & weight)
320  {
321  fNUniverses = shifts.size();
322  for(int i = 0; i < (int) shifts.size(); i++) {
323  fSpectra.push_back(new Spectrum(loader,
324  histaxis,
325  cut,
326  shifts[i],
327  weight));
328  }
329  }
330 
331  ///////////////////////////////////////////////////////////////////////
332  Multiverse::Multiverse(std::vector<SpectrumLoaderBase*> loaders,
333  const HistAxis & histaxis,
334  const Cut & cut,
335  std::vector<SystShifts> shifts,
336  const Var & weight)
337  {
338  fNUniverses = shifts.size();
339  for(int i = 0; i < (int) shifts.size(); i++) {
340  fSpectra.push_back(new Spectrum(*loaders[0],
341  histaxis,
342  cut,
343  shifts[i],
344  weight));
345  for(auto iload = 1u; iload < loaders.size(); iload++) {
346  loaders[iload]->AddSpectrum(*fSpectra.back(),
347  histaxis.GetVar1D(),
348  cut,
349  shifts[i],
350  weight);
351  }
352  }
353  }
354 
355  ///////////////////////////////////////////////////////////////////////
357  const NuTruthHistAxis & histaxis,
358  const NuTruthCut & cut,
359  const SystShifts & shift,
360  std::vector<NuTruthVar> weights,
361  const NuTruthVar & common_weight)
362  {
363  fNUniverses = weights.size();
364  for(int i = 0; i < (int) weights.size(); i++) {
365  fSpectra.push_back(new Spectrum(loader,
366  histaxis,
367  cut,
368  shift,
369  weights[i] * common_weight));
370  }
371  }
372 
373  ///////////////////////////////////////////////////////////////////////
374  Multiverse::Multiverse(std::vector<SpectrumLoaderBase*> loaders,
375  const NuTruthHistAxis & histaxis,
376  const NuTruthCut & cut,
377  const SystShifts & shift,
378  std::vector<NuTruthVar> weights,
379  const NuTruthVar & common_weight)
380  {
381  fNUniverses = weights.size();
382  for(int i = 0; i < (int) weights.size(); i++) {
383  fSpectra.push_back(new Spectrum(*loaders[0],
384  histaxis,
385  cut,
386  shift,
387  weights[i] * common_weight));
388  for(auto iload = 1u; iload < loaders.size(); iload++) {
389  loaders[iload]->AddSpectrum(*fSpectra.back(),
390  histaxis.GetVar1D(),
391  cut,
392  shift,
393  weights[i] * common_weight);
394  }
395  }
396  }
397 
398  ///////////////////////////////////////////////////////////////////////
400  const HistAxis & histaxis,
401  const Cut & cut,
402  const SystShifts & shift,
403  std::vector<Var> weights,
404  const Var & common_weight)
405  {
406  fNUniverses = weights.size();
407  for(int i = 0; i < (int) weights.size(); i++) {
408  fSpectra.push_back(new Spectrum(loader,
409  histaxis,
410  cut,
411  shift,
412  weights[i] * common_weight));
413  }
414  }
415 
416  ///////////////////////////////////////////////////////////////////////
417  Multiverse::Multiverse(std::vector<SpectrumLoaderBase*> loaders,
418  const HistAxis & histaxis,
419  const Cut & cut,
420  const SystShifts & shift,
421  std::vector<Var> weights,
422  const Var & common_weight)
423  {
424  fNUniverses = weights.size();
425  for(int i = 0; i < (int) weights.size(); i++) {
426  fSpectra.push_back(new Spectrum(*loaders[0],
427  histaxis,
428  cut,
429  shift,
430  weights[i] * common_weight));
431  for(auto iload = 1u; iload < loaders.size(); iload++) {
432  loaders[iload]->AddSpectrum(*fSpectra.back(),
433  histaxis.GetVar1D(),
434  cut,
435  shift,
436  weights[i] * common_weight);
437  }
438  }
439  }
440 
441  ///////////////////////////////////////////////////////////////////////
442  Spectrum *
444  {
445  return GetNSigmaShift(nominal, 1);
446  }
447 
448  ///////////////////////////////////////////////////////////////////////
449  Spectrum *
451  {
452  return GetNSigmaShift(nominal, -1);
453  }
454 
455  ///////////////////////////////////////////////////////////////////////
456  TH1 *
457  Multiverse::GetPlusOneSigmaShift(const TH1 * nominal)
458  {
459  return GetNSigmaShift(nominal, 1);
460  }
461 
462  ///////////////////////////////////////////////////////////////////////
463  TH1 *
465  {
466  return GetNSigmaShift(nominal, -1);
467  }
468 
469  ///////////////////////////////////////////////////////////////////////
470  Spectrum *
471  Multiverse::GetNSigmaShift(const Spectrum * nominal, double nsigma)
472  {
473  // keep track of binning info to make a correct spectrum at the end
474  double POT = nominal->POT();
475  double livetime = nominal->Livetime();
476 
477 
478  // transform to 1D hists
479  TH1 * hnominal = nominal->ToTH1(fSpectra[0]->POT());
480  TH1D * shift = (TH1D*)GetNSigmaShift(hnominal, nsigma);
481  return ToSpectrum(shift, POT, livetime);
482  }
483 
484  ///////////////////////////////////////////////////////////////////////
485  TH1 *
486  Multiverse::GetNSigmaShift(const TH1 * nominal, double nsigma)
487  {
488  if(fHists.empty())
489  InitFromSpectra();
490  // transform to 1D hists
491  TH1 * hret = (TH1*) nominal->Clone();
492  for(int i = 0; i < (int) fSpectra.size(); i++) {
493  // verify histograms are compatable
494  assert(nominal->GetNbinsX() == fHists[i]->GetNbinsX());
495  assert(nominal->GetNbinsY() == fHists[i]->GetNbinsY());
496  assert(nominal->GetNbinsZ() == fHists[i]->GetNbinsZ());
497  }
498 
499  // for each bin, calculate the shift
500  for(int i = 1; i <= nominal->GetNbinsX(); i++) {
501  for(int j = 1; j <= nominal->GetNbinsY(); j++) {
502  for(int k = 1; k <= nominal->GetNbinsZ(); k++) {
503  std::vector<double> vals;
504  for(int ihist = 0; ihist < (int) fHists.size(); ihist++) {
505  vals.push_back(fHists[ihist]->GetBinContent(i, j, k));
506  }
507  double cv = hret->GetBinContent(i, j, k);
508  hret->SetBinContent(i, j, k, BinSigma(vals, nsigma, cv));
509  }
510  }
511  }
512  return hret;
513  }
514 
515 
516  ///////////////////////////////////////////////////////////////////////
517  // Calculate the universe that marks the n sigma shift
518  // copied from ana::GenieMultiverse::BinSigma()
519  double
520  Multiverse::BinSigma(std::vector<double> events,
521  double nsigma,
522  double pivot)
523  {
524  int pivotbin = 0;
525  double pivotbincenter = 0;
526  std::sort(events.begin(), events.end());
527  for(int i = 0; i < (int) events.size() - 1; i++) {
528  if(pivot >= events[i] && pivot < events[i+1]) {
529  pivotbin = i;
530  break;
531  }
532  }
533  pivotbincenter = pivotbin+0.5;
534  double count_fraction = 2.0 * (ROOT::Math::normal_cdf(nsigma) - ROOT::Math::normal_cdf(0));
535 
536  int nsideevents = 0;
537  int lastbinindex = (int) events.size() - 1;
538  if(nsigma >= 0) nsideevents = lastbinindex - pivotbin;
539  else nsideevents = pivotbin;
540  int boundIdx = pivotbincenter + count_fraction*(double)nsideevents;
541 
542  int index = 0;
543  if(nsigma >= 0) index = std::min(boundIdx, (int)events.size() - 1);
544  else index = std::max(boundIdx, 0);
545  return events.at(index);
546  }
547 
548  ///////////////////////////////////////////////////////////////////////
549  void
550  Multiverse::SaveTo(TDirectory * dir, const std::string& name) const
551  {
552  TDirectory * tmp = gDirectory;
553 
554  dir = dir->mkdir(name.c_str()); // switch to subdir
555  dir->cd();
556 
557  TObjString("Multiverse").Write("type");
558  for(int i = 0; i < fNUniverses; i++) {
559  fSpectra[i]->SaveTo(dir, TString::Format("mvspectra_%d",i).Data());
560  }
561  TVectorD nuniv(1);
562  nuniv[0] = fNUniverses;
563  nuniv.Write("nuniv");
564 
565  dir->Write();
566  delete dir;
567 
568  tmp->cd();
569  }
570 
571  ///////////////////////////////////////////////////////////////////////
572  std::unique_ptr<Multiverse>
573  Multiverse::LoadFrom(TDirectory * dir, const std::string& name)
574  {
575  dir = dir->GetDirectory(name.c_str()); // switch to subdir
576  assert(dir);
577 
578  TObjString * ptag = (TObjString*) dir->Get("type");
579  assert(ptag);
580  assert(ptag->GetString() == "Multiverse" && "Type does not match Multiverse");
581  delete ptag;
582 
583  TVectorD * vnuniv = (TVectorD*) dir->Get("nuniv");
584  int nuniv = (*vnuniv)[0];
585  std::vector<Spectrum*> specs;
586  for(int i = 0; i < nuniv; i++) {
587  specs.push_back(Spectrum::LoadFrom(dir, TString::Format("mvspectra_%d", i).Data()).release());
588  }
589 
590  delete dir;
591 
592  std::unique_ptr<Multiverse> ret(new Multiverse(specs));
593  ret->InitFromSpectra();
594 
595  return ret;
596  }
597 
598  ///////////////////////////////////////////////////////////////////////
599  Multiverse &
601  {
602  if( this == & rhs )
603  return *this;
604  this->fNUniverses = rhs.fNUniverses;
605  for(auto i = 0u; i < rhs.fSpectra.size(); i++) {
606  this->fSpectra.push_back(new Spectrum(*rhs.fSpectra[i]));
607  }
608  for(auto i = 0u; i< rhs.fHists.size(); i++) {
609  this->fHists.push_back((TH1*) rhs.fHists[i]->Clone());
610  }
611  return *this;
612  }
613 
614  ///////////////////////////////////////////////////////////////////////
615  Multiverse &
617  {
618  if( this == & rhs )
619  return *this;
620  this->fNUniverses = rhs.fNUniverses;
621  for(auto i = 0u; i < rhs.fSpectra.size(); i++) {
622  this->fSpectra.push_back(std::move(rhs.fSpectra[i]));
623  }
624  for(auto i = 0u; i< rhs.fHists.size(); i++) {
625  this->fHists.push_back(std::move(rhs.fHists[i]));
626  }
627 
628  return *this;
629  }
630 
631  ///////////////////////////////////////////////////////////////////////
633  {
634  for(int i = 0; i < (int) fSpectra.size(); i++) {
635  delete fSpectra[i];
636  }
637  for(int i = 0; i < (int) fHists.size(); i++) {
638  delete fHists[i];
639  }
640 
641  }
642 
643  ///////////////////////////////////////////////////////////////////////
645  {
646  this->fNUniverses = rhs.fNUniverses;
647  for(auto i = 0u; i < rhs.fSpectra.size(); i++) {
648  this->fSpectra.push_back(new Spectrum(*rhs.fSpectra[i]));
649  }
650  for(auto i = 0u; i< rhs.fHists.size(); i++) {
651  this->fHists.push_back((TH1*) rhs.fHists[i]->Clone());
652  }
653 
654  }
655 
656  ///////////////////////////////////////////////////////////////////////
658  {
659  this->fNUniverses = rhs.fNUniverses;
660  for(auto i = 0u; i < rhs.fSpectra.size(); i++) {
661  this->fSpectra.push_back(std::move(rhs.fSpectra[i]));
662  }
663  for(auto i = 0u; i< rhs.fHists.size(); i++) {
664  this->fHists.push_back(std::move(rhs.fHists[i]));
665  }
666 
667  }
668 
669  ///////////////////////////////////////////////////////////////////////
670  Multiverse::Multiverse(std::vector<Spectrum*> specs)
671  {
672  this->fNUniverses = specs.size();
673  assert(fNUniverses == (int) specs.size() && "fNUniverses does not match size of spectrum vector");
674  for(int i = 0; i < this->fNUniverses; i++) {
675  fSpectra.push_back(specs[i]);
676  specs[i] = NULL;
677  }
678  }
679 
680  ///////////////////////////////////////////////////////////////////////
681  Multiverse::Multiverse(std::vector<TH1 *> hists,
682  int ndims)
683  {
684  this->fNUniverses = hists.size();
685  assert(fNUniverses == (int) hists.size() && "fNUniverses does not match size of histogram vector");
686  for(int i = 0; i < this->fNUniverses; i++) {
687  fHists.push_back(hists[i]);
688  hists[i] = NULL;
689  }
690  }
691 
692  ///////////////////////////////////////////////////////////////////////
693  void
695  {
696  assert(fSpectra.size() != 0);
697  // assume universes share binning and labels
698  fBinning = fSpectra[0]->GetBinnings();
699  fLabels = fSpectra[0]->GetLabels();
700  for(int i = 0; i < (int) fSpectra.size(); i++) {
701  fHists.push_back(ToHist(fSpectra[i]));
702  }
703  }
704 
705  ///////////////////////////////////////////////////////////////////////
706  TH1 *
708  double POT)
709  {
710  int dims = spec->NDimensions();
711  if(POT < 0)
712  POT = spec->POT();
713 
714  switch(dims) {
715  case(1) : {
716  return spec->ToTH1(POT);
717  }
718  case(2) : {
719  return spec->ToTH2(POT);
720  }
721  case(3) : {
722  return spec->ToTH3(POT);
723  }
724  default : {
725  std::cout << "Dimension of multiverse spectra undetermined." << std::endl;
726  abort();
727  }
728  } // end switch case
729  }
730 
731  ///////////////////////////////////////////////////////////////////////
732  Spectrum *
733  Multiverse::ToSpectrum(TH1 * h, double POT, double livetime)
734  {
735  int dims = h->GetDimension();
736 
737  switch(dims) {
738  case(1) : {
739  std::vector<Binning> bins = {Binning::FromTAxis(h->GetXaxis())};
740  std::vector<std::string> labels = {h->GetXaxis()->GetTitle()};
741  return new Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(((TH1D*)h)->GetArray(), h->GetNbinsX()+2)),
742  HistAxis(labels, bins), POT, livetime);
743  }
744  case(2) : {
745 
746  auto NX = h->GetNbinsX();
747  auto NY = h->GetNbinsY();
748 
749  TH1D* ret = new TH1D(UniqueName().c_str(), "",
750  NX*NY, 0, NX*NY);
751  std::vector<Binning> bins = {Binning::FromTAxis(h->GetXaxis()),
752  Binning::FromTAxis(h->GetYaxis())};
753  std::vector<std::string> labels = {h->GetXaxis()->GetTitle(),
754  h->GetYaxis()->GetTitle()};
755 
756 
757  for(auto ibinx = 0; ibinx < NX; ibinx++) {
758  for(auto ibiny = 0; ibiny < NY; ibiny++) {
759  int bin = ibinx * bins[1].NBins() + ibiny;
760  auto val = h->GetBinContent(ibinx+1, ibiny+1);
761  auto err = h->GetBinError(ibinx+1, ibiny+1);
762  ret->SetBinContent(bin+1, val);
763  ret->SetBinError(bin+1, err);
764  }
765  }
766  return new Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(ret->GetArray(), ret->GetNbinsX()+2)),
767  HistAxis(labels, bins), POT, livetime);
768  }
769  case(3) : {
770 
771  auto NX = h->GetNbinsX();
772  auto NY = h->GetNbinsY();
773  auto NZ = h->GetNbinsZ();
774  TH1D* ret = new TH1D(UniqueName().c_str(), "",
775  NX*NY*NZ, 0, NX*NY*NZ);
776  std::vector<Binning> bins = {Binning::FromTAxis(h->GetXaxis()),
777  Binning::FromTAxis(h->GetYaxis()),
778  Binning::FromTAxis(h->GetZaxis())};
779  std::vector<std::string> labels = {h->GetXaxis()->GetTitle(),
780  h->GetYaxis()->GetTitle(),
781  h->GetZaxis()->GetTitle()};
782 
783 
784  for(auto ibinx = 0; ibinx < NX; ibinx++) {
785  for(auto ibiny = 0; ibiny < NY; ibiny++) {
786  for(auto ibinz = 0; ibinz < NZ; ibinz++) {
787  int bin = ibinx * NY * NZ + ibiny * NX * ibinz;
788  auto val = h->GetBinContent(ibinx+1, ibiny+1, ibinz+1);
789  auto err = h->GetBinError(ibinx+1, ibiny+1, ibinz+1);
790  ret->SetBinContent(bin+1, val);
791  ret->SetBinError(bin+1, err);
792  }
793  }
794  }
795 
796  return new Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(ret->GetArray(), ret->GetNbinsX()+2)),
797  HistAxis(labels, bins), POT, livetime);
798 
799  }
800  default : {
801  std::cout << "Dimension " << dims << " not supported." << std::endl;
802  abort();
803  }
804  }
805 
806  }
807 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Multiverse.cxx:550
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
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:148
const Var weight
Multiverse operator/(const Spectrum &) const
Definition: Multiverse.cxx:193
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
static Binning FromTAxis(const TAxis *ax)
Definition: Binning.cxx:173
std::vector< std::string > fLabels
Definition: Multiverse.h:117
TH1 * ToHist(const Spectrum *, double=-1)
Definition: Multiverse.cxx:707
Float_t tmp
Definition: plot.C:36
std::vector< Binning > fBinning
Definition: Multiverse.h:116
TString hists[nhists]
Definition: bdt_com.C:3
std::vector< Spectrum * > fSpectra
Definition: Multiverse.h:115
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
double BinSigma(std::vector< double >, double, double)
Definition: Multiverse.cxx:520
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Double_t scale
Definition: plot.C:25
Multiverse & operator/=(const Spectrum &)
Definition: Multiverse.cxx:179
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const std::string cv[Ncv]
std::vector< TH1 * > fHists
Definition: Multiverse.h:119
Multiverse & operator=(const Multiverse &)
Definition: Multiverse.cxx:600
Multiverse operator*(const Multiverse &) const
Definition: Multiverse.cxx:167
Var weights
const double j
Definition: BetheBloch.cxx:29
unsigned int NDimensions() const
Definition: Spectrum.h:262
loader
Definition: demo0.py:10
static Multiverse CopyAndTransform(const Multiverse &rhs, const std::function< TransformFunc_t > &transform)
Definition: Multiverse.cxx:29
Multiverse & operator*=(const Multiverse &)
Definition: Multiverse.cxx:154
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
unsigned int nuniv
Definition: PlotUnfolding.C:14
void InitFromSpectra()
Definition: Multiverse.cxx:694
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Multiverse.cxx:573
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
TH1F * h2
Definition: plot.C:45
Base class for the various types of spectrum loader.
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1F * h1
const Cut cut
Definition: exporter_fd.C:30
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
Definition: Spectrum.cxx:187
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void Scale(double scale)
Definition: Multiverse.cxx:19
const std::vector< TH1 * > GetUniversesHist()
Definition: Multiverse.cxx:120
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
T min(const caf::Proxy< T > &a, T b)
Template for Cut and SpillCut.
Definition: Cut.h:15
Template for Var and SpillVar.
void Transform(const Multiverse &rhs, const std::function< ProductFunc_t > &transform)
Definition: Multiverse.cxx:72
Spectrum * ToSpectrum(TH1 *h, double POT=1, double livetime=0)
Definition: Multiverse.cxx:733
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
return_type< T_y, T_loc, T_scale >::type normal_cdf(const T_y &y, const T_loc &mu, const T_scale &sigma)
Definition: normal_cdf.hpp:39
void events(int which)
Definition: Cana.C:52
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Spectrum * GetNSigmaShift(const Spectrum *, double)
Definition: Multiverse.cxx:471
void Divide(Multiverse &, bool=false)
Definition: Multiverse.cxx:139
T GetVar1D() const
Definition: HistAxis.cxx:43
Multiverse(SpectrumLoaderBase &, const HistAxis &, const Cut &, std::vector< SystShifts >, const Var &)
Definition: Multiverse.cxx:315
void DivisionHelper(TH1 *, TH1 *, bool)
Definition: Multiverse.cxx:129
enum BeamMode string