GenieMultiverseSyst.cxx
Go to the documentation of this file.
1 #include <algorithm>
2 #include <boost/property_tree/json_parser.hpp>
3 #include <boost/property_tree/ptree.hpp>
4 #include <iostream>
5 #include <string>
6 
7 #include "CAFAna/Core/Spectrum.h"
8 #include "CAFAna/Systs/Systs.h"
12 #include "Math/ProbFunc.h"
13 #include "nugen/NuReweight/ReweightLabels.h"
14 #include "TKDE.h"
15 
16 #include "TF1.h"
17 #include "TH2.h"
18 #include "TRandom3.h"
19 #include "TObjString.h"
20 
21 
22 using namespace std;
23 
24 namespace ana
25 {
26 
27 // list of enumerations from rwgt::ReweightLabels because iterating over an
28 // enumerated list is not allowed with gcc 6
29 // This should really go into nugen/NuReweight/ReweightLabels.h at some point
30 
31 const std::list<rwgt::ReweightLabel_t> kReweightLabels({rwgt::fReweightMaNCEL,
32 rwgt::fReweightEtaNCEL,
33 rwgt::fReweightNormCCQE,
34 rwgt::fReweightNormCCQEenu,
35 rwgt::fReweightMaCCQEshape,
36 rwgt::fReweightMaCCQE,
37 rwgt::fReweightVecCCQEshape,
38 rwgt::fReweightNormCCRES,
39 rwgt::fReweightMaCCRESshape,
40 rwgt::fReweightMvCCRESshape,
41 rwgt::fReweightMaCCRES,
42 rwgt::fReweightMvCCRES,
43 rwgt::fReweightNormNCRES,
44 rwgt::fReweightMaNCRESshape,
45 rwgt::fReweightMvNCRESshape,
46 rwgt::fReweightMaNCRES,
47 rwgt::fReweightMvNCRES,
48 rwgt::fReweightMaCOHpi,
49 rwgt::fReweightR0COHpi,
50 rwgt::fReweightRvpCC1pi,
51 rwgt::fReweightRvpCC2pi,
52 rwgt::fReweightRvpNC1pi,
53 rwgt::fReweightRvpNC2pi,
54 rwgt::fReweightRvnCC1pi,
55 rwgt::fReweightRvnCC2pi,
56 rwgt::fReweightRvnNC1pi,
57 rwgt::fReweightRvnNC2pi,
58 rwgt::fReweightRvbarpCC1pi,
59 rwgt::fReweightRvbarpCC2pi,
60 rwgt::fReweightRvbarpNC1pi,
61 rwgt::fReweightRvbarpNC2pi,
62 rwgt::fReweightRvbarnCC1pi,
63 rwgt::fReweightRvbarnCC2pi,
64 rwgt::fReweightRvbarnNC1pi,
65 rwgt::fReweightRvbarnNC2pi,
66 rwgt::fReweightAhtBY,
67 rwgt::fReweightBhtBY,
68 rwgt::fReweightCV1uBY,
69 rwgt::fReweightCV2uBY,
70 rwgt::fReweightAhtBYshape,
71 rwgt::fReweightBhtBYshape,
72 rwgt::fReweightCV1uBYshape,
73 rwgt::fReweightCV2uBYshape,
74 rwgt::fReweightNormDISCC,
75 rwgt::fReweightRnubarnuCC,
76 rwgt::fReweightDISNuclMod,
77 rwgt::fReweightNC,
78 rwgt::fReweightAGKY_xF1pi,
79 rwgt::fReweightAGKY_pT1pi,
80 rwgt::fReweightFormZone,
81 rwgt::fReweightMFP_pi,
82 rwgt::fReweightMFP_N,
83 rwgt::fReweightFrCEx_pi,
84 //rwgt::fReweightFrElas_pi,
85 rwgt::fReweightFrInel_pi,
86 rwgt::fReweightFrAbs_pi,
87 rwgt::fReweightFrPiProd_pi,
88 rwgt::fReweightFrCEx_N,
89 //rwgt::fReweightFrElas_N,
90 rwgt::fReweightFrInel_N,
91 rwgt::fReweightFrAbs_N,
92 rwgt::fReweightFrPiProd_N,
93 rwgt::fReweightCCQEPauliSupViaKF,
94 rwgt::fReweightCCQEMomDistroFGtoSF,
95 rwgt::fReweightBR1gamma,
96 rwgt::fReweightBR1eta,
97 rwgt::fReweightTheta_Delta2Npi,
98 rwgt::fReweightZNormCCQE,
99 rwgt::fReweightZExpA1CCQE,
100 rwgt::fReweightZExpA2CCQE,
101 rwgt::fReweightZExpA3CCQE,
102 rwgt::fReweightZExpA4CCQE,
103 rwgt::fReweightAxFFCCQEshape});
104 
105  GenieMultiverseSpectra::GenieMultiverseSpectra(
106  unsigned int nuniverses,
107  SpectrumLoaderBase& loader,
108  const HistAxis& axis,
109  const Cut& cut,
110  const SystShifts& shift,
111  const Var& wei,
112  string config_pathname): fNUniverses(nuniverses)
113  {
114  BuildKnobConfigTable(config_pathname);
115 
116  // universe with index 0 is the nominal universe
117  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
118  for(unsigned int i = 0; i < nuniverses-1; i++)
119  fSpectra.push_back(CreateUniverse(loader, axis, cut, shift, wei));
120  }
121 
123  unsigned int nuniverses,
124  SpectrumLoaderBase& loader,
125  const HistAxis& axis,
126  const Cut& cut,
127  const SystShifts& shift,
128  const Var& wei,
129  std::vector<const ISyst*> systs): fNUniverses(nuniverses)
130  {
131  // universe with index 0 is the nominal universe
132  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
133  for(unsigned int i = 0; i < nuniverses-1; i++)
134  fSpectra.push_back(CreateUniverse(loader, axis, cut, systs, wei));
135  }
136 
138  unsigned int nuniverses,
139  SpectrumLoaderBase& loader,
140  const NuTruthHistAxis& axis,
141  const NuTruthCut& cut,
142  const SystShifts& shift,
143  const NuTruthVar& wei,
144  string config_pathname): fNUniverses(nuniverses)
145  {
146  BuildKnobConfigTable(config_pathname);
147 
148  // universe with index 0 is the nominal universe
149  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
150  for(unsigned int i = 0; i < nuniverses-1; i++)
151  fSpectra.push_back(CreateUniverse(loader, axis, cut, shift, wei));
152  }
153 
155  unsigned int nuniverses,
156  SpectrumLoaderBase& loader,
157  const NuTruthHistAxis& axis,
158  const NuTruthCut& cut,
159  const SystShifts& shift,
160  const NuTruthVar& wei,
161  vector<const ISyst*> systs): fNUniverses(nuniverses)
162  {
163  // universe with index 0 is the nominal universe
164  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
165  for(unsigned int i = 0; i < nuniverses-1; i++)
166  fSpectra.push_back(CreateUniverse(loader, axis, cut, systs, wei));
167  }
168 
170  std::vector<std::unique_ptr<ana::Spectrum>>& spectra, bool fromfile): fLoadFromFile(fromfile)
171  {
173  for(auto& it: spectra)
174  fSpectra.push_back(it.release());
175  }
176 
177  GenieMultiverseSpectra::GenieMultiverseSpectra(vector<std::unique_ptr<Spectrum>>& spectra)
178  {
179  for(auto& it: spectra)
180  fSpectra.push_back(it.release());
181  fNUniverses = fSpectra.size();
182  }
183 
184  //----------------------------------------------------------------------------
185  pair<TH2*, TH2*> GenieMultiverseSpectra::GetSigmaHistograms(vector<TH2*> huniv, BandOptions opt)
186  {
187  // create up and down histograms
188  TH2* hup = (TH2*)huniv[0]->Clone("hup");
189  TH2* hdown = (TH2*)huniv[0]->Clone("hdown");
190 
191  // loop through all bins
192  int nbinsx = hup->GetNbinsX();
193  int nbinsy = hup->GetNbinsY();
194  for(int ix = 0; ix <= nbinsx; ix++){
195  for(int iy = 0; iy <= nbinsy; iy++)
196  {
197  // collect counts from all universes
198  vector<float> counts_in_one_bin;
199  for(unsigned int iu = 0; iu < huniv.size(); iu++)
200  counts_in_one_bin.push_back(huniv[iu]->GetBinContent(ix, iy));
201 
202  // determine from where the error band is drawn
203  double pivot = 0.;
204  double mean = TMath::Mean(counts_in_one_bin.begin(), counts_in_one_bin.end());
205  if(opt == kBandFromMean) pivot = mean;
206  else if(opt == kBandFromNominal || opt == kNoBand) pivot = huniv[0]->GetBinContent(ix, iy);
207  else if(opt == kBandFromMedian){
208  float* a = &counts_in_one_bin[0];
209  pivot = TMath::Median(counts_in_one_bin.size(), a);
210  }
211  else if(opt == kBandFromMPV){
212  double a[counts_in_one_bin.size()];
213  for(unsigned int i = 0; i < counts_in_one_bin.size(); i++) a[i] = counts_in_one_bin[i];
214  double min = *min_element(counts_in_one_bin.begin(), counts_in_one_bin.end());
215  double max = *max_element(counts_in_one_bin.begin(), counts_in_one_bin.end());
216  TKDE kde(counts_in_one_bin.size(), a, min, max, "kUnbinned");
217  pivot = kde.GetFunction(counts_in_one_bin.size(), min, max)->GetMaximumX();
218  }
219 
220  // draw band
221  hup->SetBinContent(ix, iy, BinSigma(counts_in_one_bin, 1, pivot));
222  hdown->SetBinContent(ix, iy, BinSigma(counts_in_one_bin, -1, pivot));
223  }
224  }
225 
226  return pair<TH2*, TH2*>(hup, hdown);
227  }
228 
230  {
231  if(!fLoadFromFile)
232  for (auto it = fSpectra.begin(); it != fSpectra.end(); ++it) delete *it;
233  delete fUpperRMS;
234  delete fLowerRMS;
235  delete fMean;
236  }
237 
239  {
241  return fMean;
242  }
243 
245  {
246  Spectrum* nominal_u = 0;
247  try {
248  nominal_u = fSpectra.at(0);
249  } catch(const std::out_of_range& e) {
250  cerr << "Nominal universe is accessed before existence." << endl;
251  }
252 
253  return nominal_u;
254  }
255 
256  TGraphAsymmErrors* GenieMultiverseSpectra::ToAreaNormalizedTH1(double pot, int col,
257  int errCol, float headroom, bool newaxis) const
258  {
259  if(col == -1){
260  col = kRed;
261  errCol = kRed-10;
262  }
263  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
264 
265  TH1* nom = Nominal()->ToTH1(pot);
266 
267  std::vector<TH1*> ups, dns;
268  ups.push_back(fUpperRMS->ToTH1(pot));
269  dns.push_back(fLowerRMS->ToTH1(pot));
270 
271  nom->SetLineColor(col);
272  nom->GetXaxis()->CenterTitle();
273  nom->GetYaxis()->CenterTitle();
274  //~ if(newaxis) nom->Draw("hist ]["); // Set the axes up
275 
276  double yMax = nom->GetBinContent(nom->GetMaximumBin());
277 
278  TGraphAsymmErrors* g = new TGraphAsymmErrors;
279 
280  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
281  const double y = nom->GetBinContent(binIdx);
282  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
283 
284  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
285 
286  double errUp = 0;
287  double errDn = 0;
288 
289  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
290  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
291  double lo = y-dns[systIdx]->GetBinContent(binIdx);
292 
293  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
294 
295  errUp += hi*hi;
296  errDn += lo*lo;
297 
298  // TODO: what happens if they're both high or both low?
299  } // end for systIdx
300 
301 
302  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
303  } // end for i
304 
305  g->SetFillColor(errCol);
306  //~ g->Draw("e2 same");
307  g->GetYaxis()->SetRangeUser(0, headroom*yMax);
308  nom->GetYaxis()->SetRangeUser(0, headroom*yMax);
309 
310  //~ nom->Draw("hist ][ same");
311  return g;
312  }
313 
315  {
316  return ToTH1(Nominal()->POT());
317  }
318 
320  {
321  return ToTH2(Nominal()->POT());
322  }
323 
324  // return a TH1* with RMS of all universes as the error band
326  {
327  TH1* hNominal = Nominal()->ToTH1(pot);
328  vector<TH1*> hUniverses;
329  map< int, vector<float> > counts_in_bins;
330  for(auto& it: fSpectra)
331  {
332  TH1* curhist = it->ToTH1(pot);
333  hUniverses.push_back(curhist);
334  int nbinsx = curhist->GetNbinsX();
335  for(int i = 1; i <= nbinsx; i++)
336  counts_in_bins[i].push_back(curhist->GetBinContent(i));
337  }
338 
339  for(auto& it: counts_in_bins)
340  {
341  hNominal->SetBinError(it.first, TMath::RMS(it.second.begin(), it.second.end()));
342  //~ cout << TMath::RMS(it.second.begin(), it.second.end()) << endl;
343  }
344  hNominal->SetLineColor(kRed);
345  hNominal->SetMarkerColor(kRed);
346 
347  return hNominal;
348  }
349 
350  // return a TH2* with RMS of all universes as the error band
352  {
353  TH2* hNominal = Nominal()->ToTH2(pot);
354  vector<TH2*> hUniverses;
355  map< std::pair<int,int>, vector<float> > counts_in_bins;
356  for(auto& it: fSpectra)
357  {
358  TH2* curhist = it->ToTH2(pot);
359  hUniverses.push_back(curhist);
360  int nbinsx = curhist->GetNbinsX();
361  int nbinsy = curhist->GetNbinsY();
362  for(int i = 1; i <= nbinsx; i++)
363  {
364  for(int j = 1; j<= nbinsy; j++)
365  {
366  counts_in_bins[std::make_pair(i,j)].push_back(curhist->GetBinContent(i,j));
367  }
368  }
369  }
370 
371  for(auto& it: counts_in_bins)
372  {
373  hNominal->SetBinError(it.first.first,it.first.second, TMath::RMS(it.second.begin(), it.second.end()));
374  //~ cout << TMath::RMS(it.second.begin(), it.second.end()) << endl;
375  }
376 
377  return hNominal;
378  }
379 
380  // return upper boundary
382  {
384  return fUpperRMS;
385  }
386 
387  // return upper boundary by sigma (34% away from the mean)
389  {
390  if(fCurBandOpt != opt) FindBandBoundaries(opt);
391  return fUpperSigma;
392  }
393 
394  // return boundaries by specified number of sigmas, accounting for different # of events on either side.
395  float GenieMultiverseSpectra::BinSigma(vector<float> events, float nsigmas, float pivot)
396  {
397  int pivotbin = 0;
398  double pivotbincenter = 0;
399  sort(events.begin(), events.end());
400  for(unsigned int i = 0; i < events.size()-1; i++)
401  if(pivot >= events[i] && pivot < events[i+1]){
402  pivotbin = i;
403  break;
404  }
405  pivotbincenter = pivotbin+0.5;
406  double count_fraction = 2.0 * (ROOT::Math::normal_cdf(nsigmas)-ROOT::Math::normal_cdf(0));
407 
408  int nsideevents = 0;
409  int lastbinindex = (int)events.size() - 1;
410  if (nsigmas >= 0) nsideevents = lastbinindex - pivotbin;
411  else nsideevents = pivotbin;
412  int boundIdx = pivotbincenter + count_fraction*(double)nsideevents;
413 
414  int index = 0;
415  if(nsigmas >= 0) index = min(boundIdx, (int)events.size()-1);
416  else index = max(boundIdx, 0);
417 
418  return events.at(index);
419  }
420 
421  // build a knob configuration table from a config file
423  {
424  // Create an empty property tree object
425  //~ using boost::property_tree::ptree;
426  //~ ptree pt;
427 
428  // check file existence
429  ifstream f(fpn.c_str());
430  if(!f.good()){
431  cerr << "ERROR: GENIE knob configuration file \033[1;31m" << fpn << "\033[0m is not found." << endl;
432  exit(0);
433  }
434 
435  // read configuration from the json file
436  // This line is where the boost bug emerges.
437  //~ boost::property_tree::json_parser::read_json(fpn, pt);
438 
439  string knobname;
440  int sample_mode;
441  int index = 0;
442  while(f >> knobname >> sample_mode)
443  fKnobConfigTable.insert(knob_sampling_mode(index++,knobname,sample_mode));
444 
445  // print out key values
446  //~ for(auto& nameIt: fKnobConfigTable.get<0>()) cout << nameIt.id << " " << nameIt.name << " " << nameIt.mode << endl;
447  }
448 
449  // function to create one universe
451  const HistAxis& axis,
452  const Cut& cut,
453  const SystShifts& shift,
454  const Var& wei)
455  {
456  const vector<const ISyst*> systs = getAllXsecNuTruthSysts_2018();
457 
458  map<const ISyst*, double> knob_weight_table;
459  for(auto& idIt: fKnobConfigTable.get<0>())
460  {
461  TRandom3 rand(fUniverseSeed++);
462  rwgt::ReweightLabel_t knoblabel = (rwgt::ReweightLabel_t)idIt.id;
463  double nsigma = 0;
464  if(idIt.mode == 0) continue;
465  else if(idIt.mode == 1) nsigma = rand.Gaus();
466  else if(idIt.mode == 2) nsigma = rand.Integer(2);
467  // For tags later than 2017/08/14, when the underlying framework
468  // underwent an overhaul, this conditional check is needed.
469 
470  // Note that in GetGenieKnobSyst function in XSecSysts.cxx,
471  // rwgt::kReweightNull is excluded. Therefore a call with argument
472  // rwgt::kReweightNull will lead to seg fault. Skip it!
473  if(knoblabel == rwgt::kReweightNull) continue;
474  if(find(systs.begin(), systs.end(), GetGenieKnobSyst(knoblabel)) != systs.end())
475  knob_weight_table[GetGenieKnobSyst(knoblabel)] = nsigma;
476 
477  }
478  //~ for(auto& it: knob_weight_table) cout << "number of sigmas: " << it.second << endl;
479 
480  auto spec = new Spectrum(loader, axis, cut, knob_weight_table, wei);
481  fShiftTable.push_back(std::move(knob_weight_table));
482  return spec;
483  }
484 
485  // Function to create one universe
487  const HistAxis& axis,
488  const Cut& cut,
489  std::vector<const ISyst*> systs,
490  const Var& wei)
491  {
492  map<const ISyst*, double> knob_weight_table;
493 
494  for (const auto & syst: systs)
495  {
496  TRandom3 rand(fUniverseSeed++);
497  knob_weight_table[syst] = rand.Gaus();
498  }
499 
500  auto spec = new Spectrum(loader, axis, cut, knob_weight_table, wei);
501  fShiftTable.push_back(std::move(knob_weight_table));
502  return spec;
503  }
504 
505  // Function to create one universe for NuTruth variables.
507  const NuTruthHistAxis& axis,
508  const NuTruthCut& cut,
509  const SystShifts& shift,
510  const NuTruthVar& wei)
511  {
512  const vector<const ISyst*> systs = getAllXsecNuTruthSysts_2017();
513 
514  map<const ISyst*, double> knob_weight_table;
515  for(auto& idIt: fKnobConfigTable.get<0>())
516  {
517  TRandom3 rand(fUniverseSeed++);
518  rwgt::ReweightLabel_t knoblabel = (rwgt::ReweightLabel_t)idIt.id;
519  double nsigma = 0;
520  if(idIt.mode == 0) continue;
521  else if(idIt.mode == 1) nsigma = rand.Gaus();
522  else if(idIt.mode == 2) nsigma = rand.Integer(2);
523  // For tags later than 2017/08/14, when the underlying framework
524  // underwent an overhaul, this conditional check is needed.
525 
526  // Note that in GetGenieKnobSyst function in XSecSysts.cxx,
527  // rwgt::kReweightNull is excluded. Therefore a call with argument
528  // rwgt::kReweightNull will lead to seg fault. Skip it!
529  if(knoblabel == rwgt::kReweightNull) continue;
530  if(find(systs.begin(), systs.end(), GetGenieKnobSyst(knoblabel)) != systs.end())
531  knob_weight_table[GetGenieKnobSyst(knoblabel)] = nsigma;
532  }
533 
534  return new Spectrum(loader, axis, cut, SystShifts(knob_weight_table), wei);
535  }
536 
537  // Function to create one universe for NuTruth variables.
539  const NuTruthHistAxis& axis,
540  const NuTruthCut& cut,
541  vector<const ISyst*> systs,
542  const NuTruthVar& wei)
543  {
544  map<const ISyst*, double> knob_weight_table;
545  for(auto knoblabel : ana::kReweightLabels)
546  {
547  TRandom3 rand(fUniverseSeed++);
548  // Note that in GetGenieKnobSyst function in XSecSysts.cxx,
549  // rwgt::kReweightNull is excluded. Therefore a call with argument
550  // rwgt::kReweightNull will lead to seg fault. Skip it!
551  if(knoblabel == rwgt::kReweightNull) continue;
552  if(find(systs.begin(), systs.end(), GetGenieKnobSyst(knoblabel)) != systs.end())
553  knob_weight_table[GetGenieKnobSyst(knoblabel)] = rand.Gaus();
554  }
555 
556  return new Spectrum(loader, axis, cut, SystShifts(knob_weight_table), wei);
557  }
558 
559  // function to fill the boudary spactra encircling the error band
561  {
562  Eigen::ArrayXd aNominal = fSpectra[0]->GetEigen(fSpectra[0]->POT());
563  map<int, vector<float> > events_in_bins;
564  for(auto& it: fSpectra)
565  {
566  Eigen::ArrayXd curhist = it->GetEigen(it->POT());
567  int nbinsx = curhist.size()-2;
568  for(int i = 1; i <= nbinsx; i++)
569  events_in_bins[i].push_back(curhist[i]);
570  }
571 
572  Eigen::ArrayXd aUpRMS = aNominal;
573  Eigen::ArrayXd aUpSigma = aNominal;
574  Eigen::ArrayXd aLowRMS = aNominal;
575  Eigen::ArrayXd aLowSigma = aNominal;
576  Eigen::ArrayXd aMean = aNominal;
577  int binIt = 1;
578  for(auto& it: events_in_bins){
579 
580  // Find the pivot point from which the boundaries are calculated.
581  double pivot = 0;
582  double mean = TMath::Mean(it.second.begin(), it.second.end());
583  if(opt == kBandFromMean)
584  pivot = mean;
585  else if(opt == kBandFromNominal)
586  pivot = aNominal[binIt];
587  else if(opt == kBandFromMedian){
588  float* a = &it.second[0];
589  pivot = TMath::Median(it.second.size(), a);
590  }
591  else if(opt == kBandFromMPV){
592  double a[it.second.size()];
593  for(unsigned int i = 0; i < it.second.size(); i++) a[i] = it.second[i];
594  double min = *min_element(it.second.begin(), it.second.end());
595  double max = *max_element(it.second.begin(), it.second.end());
596  TKDE kde(fNUniverses, a, min, max, "kUnbinned");
597  pivot = kde.GetFunction(fNUniverses, min, max)->GetMaximumX();
598  }
599 
600  // Deal with no band option
601  if(opt != kNoBand){
602  double rms = TMath::RMS(it.second.begin(), it.second.end());
603  aUpRMS[binIt] = mean+rms;
604  aUpSigma[binIt] = BinSigma(it.second, 1, pivot);
605  aLowRMS[binIt] = mean-rms;
606  aLowSigma[binIt] = BinSigma(it.second, -1, pivot);
607  }
608  else{
609  aUpRMS[binIt] = aNominal[binIt];
610  aUpSigma[binIt] = aNominal[binIt];
611  aLowRMS[binIt] = aNominal[binIt];
612  aLowSigma[binIt] = aNominal[binIt];
613  }
614  aMean[binIt] = mean;
615 
616  binIt++;
617  }
618 
619  // release memory if they are already allocated
620  if(fUpperRMS) delete fUpperRMS;
621  if(fLowerRMS) delete fLowerRMS;
622  if(fMean) delete fMean;
623  if(fUpperSigma) delete fUpperSigma;
624  if(fLowerSigma) delete fLowerSigma;
625  HistAxis ax(fSpectra[0]->GetLabels(), fSpectra[0]->GetBinnings());
626  fUpperRMS = new Spectrum(std::move(aUpRMS), ax, fSpectra[0]->POT(), fSpectra[0]->Livetime());
627  fLowerRMS = new Spectrum(std::move(aLowRMS), ax, fSpectra[0]->POT(), fSpectra[0]->Livetime());
628  fMean = new Spectrum(std::move(aMean), ax, fSpectra[0]->POT(), fSpectra[0]->Livetime());
629  fUpperSigma = new Spectrum(std::move(aUpSigma), ax, fSpectra[0]->POT(), fSpectra[0]->Livetime());
630  fLowerSigma = new Spectrum(std::move(aLowSigma), ax, fSpectra[0]->POT(), fSpectra[0]->Livetime());
631  fBoundariesFound = true;
632  fCurBandOpt = opt;
633  // area normalized spectra
634  // Based on Kirk Bay's code
635  // http://nusoft.fnal.gov/nova/novasoft/doxygen/html/sa__nd__data__mc__systs_8C_source.html
636  // function DrawMCNormSyst()
637  // Actual code removed. The class GenieMultiverseNormalizedSpectra will take care of this automatically.
638  }
639 
640  // return lower boundary
642  {
644  return fLowerRMS;
645  }
646 
647  // return lower boundary by sigma (34% away from the mean)
649  {
650  if(fCurBandOpt != opt) FindBandBoundaries(opt);
651  return fLowerSigma;
652  }
653 
654  //----------------------------------------------------------------------
655  void GenieMultiverseSpectra::SaveTo(TDirectory* dir, const std::string& name) const {
656  TDirectory *tmp = gDirectory;
657 
658  dir = dir->mkdir(name.c_str()); // switch to subdir
659  dir->cd();
660 
661  //~ TVectorD universes(1);
662  //~ universes[0] = {(Double_t)fSpectra.size()};
663  //~ universes.Write("nUniverses");
664  // Could not get the above code to store a resonable unsigned int value,
665  // which leads to failure when executing hadd_cafana.
666  // Among the hadd_cafana data types, TObjString types are also not added.
667  // Gonna try this.
668  TObjString universes(Form("%lu",fSpectra.size()));
669  universes.Write("nUniverses");
670 
671  int count = 0;
672  for(auto& it: fSpectra){
673  it->SaveTo(dir, "universe_"+std::to_string(count));
674  count++;
675  }// end loop over spectra
676 
677  dir->Write();
678  delete dir;
679 
680  tmp->cd();
681  }// end of SaveTo
682 
683 
684  std::unique_ptr<GenieMultiverseSpectra> GenieMultiverseSpectra::LoadFrom(TDirectory* dir, const std::string& name)
685  {
686  dir = dir->GetDirectory(name.c_str()); // switch to subdir
687  assert(dir);
688 
689  //~ TObjString* tag = (TObjString*)dir->Get("type");
690  //~ assert(tag);
691  //~ assert(tag->GetString() == "GenieMultiverseSpectra");
692  // The above code is departure from the Spectrum::LoadFrom implementation.
693  // In multiverse root directory the only information needed is the number
694  // of universes created.
695 
696  TObjString* universes = (TObjString*)dir->Get("nUniverses");
697  unsigned int nUniverses = universes->GetString().Atoi();
698 
699 
700  std::vector<std::unique_ptr<ana::Spectrum>> spectra;
701  for(unsigned int i = 0; i < nUniverses; i++){
702  spectra.push_back(ana::LoadFrom<Spectrum>(dir, "universe_"+std::to_string(i)));
703  }
704 
705  delete dir;
706 
707  return std::make_unique<GenieMultiverseSpectra>(nUniverses, spectra, true);
708  }// end of LoadFrom
709 
711  unsigned int nuniverses,
712  std::string config_pathname):
713  fUniverseSeed(1001)
714  {
715  // Retrieve knob configuration from file.
716  BuildKnobConfigTable(config_pathname);
717  // Generate tables of knobs for all universes.
718  // The first entry corresponds to the nominal universe.
719  for(unsigned int uIdx = 0; uIdx < nuniverses; uIdx++)
720  {
721  map<const ISyst*, double> knob_weight_table;
722  for(auto& idIt: fKnobConfigTable.get<0>())
723  {
724  rwgt::ReweightLabel_t knoblabel = (rwgt::ReweightLabel_t)idIt.id;
725  double nsigma = 0;
726  if(uIdx > 0){
727  TRandom3 rand(fUniverseSeed++);
728  if(idIt.mode == 0) continue;
729  else if(idIt.mode == 1) nsigma = rand.Gaus();
730  else if(idIt.mode == 2) nsigma = rand.Integer(2);
731  }
732  // For tags later than 2017/08/14, when the underlying framework
733  // underwent an overhaul, this conditional check is needed.
734  vector<const ISyst*> systs = getAllXsecNuTruthSysts_2017();
735  if(find(systs.begin(), systs.end(), GetGenieKnobSyst(knoblabel)) != systs.end())
736  knob_weight_table[GetGenieKnobSyst(knoblabel)] = nsigma;
737  }
738  fShiftTables.push_back(knob_weight_table);
739  fSystShifts.push_back(SystShifts(knob_weight_table));
740  }
741  }
742 
743 //------------------------------------------------------------------------------
745  unsigned int nuniverses,
746  vector<const ISyst*> sel_systs):
747  fUniverseSeed(1001)
748  {
749  // Generate tables of knobs for all universes.
750  // The first entry corresponds to the nominal universe.
751  for(unsigned int uIdx = 0; uIdx < nuniverses; uIdx++)
752  {
753  map<const ISyst*, double> knob_weight_table;
754  for(auto& knobIt: sel_systs)
755  {
756  double nsigma = 0;
757  if(uIdx > 0){
758  TRandom3 rand(fUniverseSeed++);
759  nsigma = rand.Gaus();
760  }
761  knob_weight_table[knobIt] = nsigma;
762  }
763  fShiftTables.push_back(knob_weight_table);
764  fSystShifts.push_back(SystShifts(knob_weight_table));
765  }
766  fIsNuTruth = false;
767  }
768 
770  {
771  // check file existence
772  ifstream f(fpn.c_str());
773  if(!f.good()){
774  cerr << "ERROR: GENIE knob configuration file \033[1;31m" << fpn << "\033[0m is not found." << endl;
775  exit(0);
776  }
777 
778  string knobname;
779  int sample_mode;
780  int index = 0;
781  while(f >> knobname >> sample_mode)
782  fKnobConfigTable.insert(knob_sampling_mode(index++,knobname,sample_mode));
783  }
784 
785  //
786  // Implementation of shape-only spectra.
787  //
788  // Note that the normalization only happens in the later stages, not in the
789  // constructor because they are not instanced by calling Go() yet.
791  unsigned int nuniverses,
792  SpectrumLoaderBase& loader,
793  const HistAxis& axis,
794  const Cut& cut,
795  const SystShifts& shift,
796  const Var& wei,
797  string config_pathname)
798  {
799  BuildKnobConfigTable(config_pathname);
800 
801  // universe with index 0 is the nominal universe
802  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
803  for(unsigned int i = 0; i < nuniverses-1; i++)
804  fSpectra.push_back(CreateUniverse(loader, axis, cut, shift, wei));
805  }
806 
808  unsigned int nuniverses,
809  SpectrumLoaderBase& loader,
810  const HistAxis& axis,
811  const Cut& cut,
812  const SystShifts& shift,
813  const Var& wei,
814  vector<const ISyst*> systs)
815  {
816  // universe with index 0 is the nominal universe
817  fSpectra.push_back(new Spectrum(loader, axis, cut, shift, wei));
818  for(unsigned int i = 0; i < nuniverses-1; i++)
819  fSpectra.push_back(GenieMultiverseSpectra::CreateUniverse(loader, axis, cut, systs, wei));
820  }
821 
823  SpectrumLoaderBase& loader,
824  const HistAxis& axis,
825  const Cut& cut,
826  const SystShifts& shift,
827  const Var& wei)
828  {
829  map<const ISyst*, double> knob_weight_table;
830  for(auto& idIt: fKnobConfigTable.get<0>())
831  {
832  TRandom3 rand(fUniverseSeed++);
833  rwgt::ReweightLabel_t knoblabel = (rwgt::ReweightLabel_t)idIt.id;
834  double nsigma = 0;
835  if(idIt.mode == 0) continue;
836  else if(idIt.mode == 1) nsigma = rand.Gaus();
837  else if(idIt.mode == 2) nsigma = rand.Integer(2);
838  // For tags later than 2017/08/14, when the underlying framework
839  // underwent an overhaul, this conditional check is needed.
840  vector<const ISyst*> systs = getAllXsecNuTruthSysts_2017();
841  if(find(systs.begin(), systs.end(), GetGenieKnobSyst(knoblabel)) != systs.end())
842  knob_weight_table[GetGenieKnobSyst(knoblabel)] = nsigma;
843  }
844  Spectrum* normSpec = new Spectrum(loader, axis, cut, knob_weight_table, wei);
845 
846  return normSpec;
847  }
848 
849  std::unique_ptr<GenieMultiverseNormalizedSpectra> GenieMultiverseNormalizedSpectra::LoadFrom(TDirectory* dir, const std::string& name)
850  {
851  dir = dir->GetDirectory(name.c_str()); // switch to subdir
852  assert(dir);
853 
854  TObjString* universes = (TObjString*)dir->Get("nUniverses");
855  unsigned int nUniverses = universes->GetString().Atoi();
856 
857  std::vector<std::unique_ptr<ana::Spectrum>> spectra;
858  for(unsigned int i = 0; i < nUniverses; i++){
859  spectra.push_back(ana::LoadFrom<Spectrum>(dir, "universe_"+std::to_string(i)));
860  }
861 
862  delete dir;
863 
864  return std::make_unique<GenieMultiverseNormalizedSpectra>(nUniverses, spectra, true);
865  }
866 
868  {
869  if(!fSpectraNormalized) NormalizeSpectra();
871  }
872 
874  {
875  fNominalCount = fSpectra[0]->Integral(fSpectra[0]->POT());
876  for(unsigned int i = 1; i < fSpectra.size(); i++)
877  {
878  Spectrum* normSpec = fSpectra[i];
879  // Scale this spectrum to the area of the nominal spectrum.
880  double scale = fNominalCount/normSpec->Integral(normSpec->POT());
881  normSpec->Scale(scale);
882  }
883  fSpectraNormalized = true;
884  }
885 
887  {
888  if(!fSpectraNormalized) NormalizeSpectra();
890  }
891 
893  {
894  if(!fSpectraNormalized) NormalizeSpectra();
896  }
897 }
std::vector< Spectrum * > fSpectra
static float BinSigma(std::vector< float > events, float nsigmas, float mean=-1)
std::vector< const ISyst * > getAllXsecNuTruthSysts_2018()
TSpline3 lo("lo", xlo, ylo, 12,"0")
enum BeamMode kRed
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
Spectrum * fLowerRMS
spectra for lower rms
std::vector< std::map< const ISyst *, double > > fShiftTables
Container of shift tables for all universes.
set< int >::iterator it
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const std::list< rwgt::ReweightLabel_t > kReweightLabels({rwgt::fReweightMaNCEL, rwgt::fReweightEtaNCEL, rwgt::fReweightNormCCQE, rwgt::fReweightNormCCQEenu, rwgt::fReweightMaCCQEshape, rwgt::fReweightMaCCQE, rwgt::fReweightVecCCQEshape, rwgt::fReweightNormCCRES, rwgt::fReweightMaCCRESshape, rwgt::fReweightMvCCRESshape, rwgt::fReweightMaCCRES, rwgt::fReweightMvCCRES, rwgt::fReweightNormNCRES, rwgt::fReweightMaNCRESshape, rwgt::fReweightMvNCRESshape, rwgt::fReweightMaNCRES, rwgt::fReweightMvNCRES, rwgt::fReweightMaCOHpi, rwgt::fReweightR0COHpi, rwgt::fReweightRvpCC1pi, rwgt::fReweightRvpCC2pi, rwgt::fReweightRvpNC1pi, rwgt::fReweightRvpNC2pi, rwgt::fReweightRvnCC1pi, rwgt::fReweightRvnCC2pi, rwgt::fReweightRvnNC1pi, rwgt::fReweightRvnNC2pi, rwgt::fReweightRvbarpCC1pi, rwgt::fReweightRvbarpCC2pi, rwgt::fReweightRvbarpNC1pi, rwgt::fReweightRvbarpNC2pi, rwgt::fReweightRvbarnCC1pi, rwgt::fReweightRvbarnCC2pi, rwgt::fReweightRvbarnNC1pi, rwgt::fReweightRvbarnNC2pi, rwgt::fReweightAhtBY, rwgt::fReweightBhtBY, rwgt::fReweightCV1uBY, rwgt::fReweightCV2uBY, rwgt::fReweightAhtBYshape, rwgt::fReweightBhtBYshape, rwgt::fReweightCV1uBYshape, rwgt::fReweightCV2uBYshape, rwgt::fReweightNormDISCC, rwgt::fReweightRnubarnuCC, rwgt::fReweightDISNuclMod, rwgt::fReweightNC, rwgt::fReweightAGKY_xF1pi, rwgt::fReweightAGKY_pT1pi, rwgt::fReweightFormZone, rwgt::fReweightMFP_pi, rwgt::fReweightMFP_N, rwgt::fReweightFrCEx_pi, rwgt::fReweightFrInel_pi, rwgt::fReweightFrAbs_pi, rwgt::fReweightFrPiProd_pi, rwgt::fReweightFrCEx_N, rwgt::fReweightFrInel_N, rwgt::fReweightFrAbs_N, rwgt::fReweightFrPiProd_N, rwgt::fReweightCCQEPauliSupViaKF, rwgt::fReweightCCQEMomDistroFGtoSF, rwgt::fReweightBR1gamma, rwgt::fReweightBR1eta, rwgt::fReweightTheta_Delta2Npi, rwgt::fReweightZNormCCQE, rwgt::fReweightZExpA1CCQE, rwgt::fReweightZExpA2CCQE, rwgt::fReweightZExpA3CCQE, rwgt::fReweightZExpA4CCQE, rwgt::fReweightAxFFCCQEshape})
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< std::map< const ISyst *, double > > fShiftTable
OStream cerr
Definition: OStream.cxx:7
void SaveTo(TDirectory *dir, const std::string &name)
Save the Spectrum structure of all universe to a TFile directory.
static std::unique_ptr< GenieMultiverseNormalizedSpectra > LoadFrom(TDirectory *dir, const std::string &name)
Load the Spectrum structure of all universe from a TFile directory.
Spectrum * CreateUniverse(SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
Float_t tmp
Definition: plot.C:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
declarations for knob configuration table
void SaveTo(TDirectory *dir, const std::string &name) const
Save the Spectrum structure of all universe to a TFile directory.
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
std::vector< Spectrum * > fSpectra
spectra of all universes including the central value:
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
const Spectrum * LowerRMS() const
unsigned int fNUniverses
number of universes to generate
void Cut(double x)
Definition: plot_outliers.C:1
Double_t scale
Definition: plot.C:25
TSpline3 hi("hi", xhi, yhi, 18,"0")
static std::pair< TH2 *, TH2 * > GetSigmaHistograms(std::vector< TH2 * >, BandOptions opt=kBandFromNominal)
Calculate error band given 2D histograms.
Spectrum * CreateUniverse(SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const SystShifts &shift=kNoShift, const Var &wei=kUnweighted)
void BuildKnobConfigTable(std::string config_pathname)
void BuildKnobConfigTable(std::string config_pathname)
Function for building knob configuration table.
Spectrum * fMean
spectra for mean
Int_t col[ntarg]
Definition: Style.C:29
const double a
GenieMultiverseNormalizedSpectra(unsigned int nuniverses, SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &cut, const SystShifts &shift, const Var &wei=kUnweighted, std::string config_pathname="knob_config.txt")
base_types push_back(int_type())
const Spectrum * Mean() const
#define pot
knob_sampling_mode_set fKnobConfigTable
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal)
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal)
const double j
Definition: BetheBloch.cxx:29
Spectrum * fLowerSigma
spectra for lower sigma
TGraphAsymmErrors * ToAreaNormalizedTH1(double pot, int col=-1, int errCol=-1, float headroom=1.3, bool newaxis=true) const
loader
Definition: demo0.py:10
Spectrum * fUpperRMS
spectra for upper rms
std::vector< float > Spectrum
Definition: Constants.h:610
double POT() const
Definition: Spectrum.h:227
const Spectrum * UpperRMS() const
return RMS extremes
std::vector< double > POT
void out_of_range(const char *function, int max, int index, const char *msg1="", const char *msg2="")
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Base class for the various types of spectrum loader.
GenieMultiverseParameters(unsigned int nuniverses, std::string config_pathname="knob_config.txt")
This class generates the parameters of the multi-universes. No more, no less.
std::vector< SystShifts > fSystShifts
const Cut cut
Definition: exporter_fd.C:30
knob_sampling_mode_set fKnobConfigTable
knob configuration table
Spectrum * fUpperSigma
spectra for upper sigma
Int_t nbinsx
Definition: plot.C:23
TDirectory * dir
Definition: macro.C:5
exit(0)
assert(nhit_max >=nhit_nbins)
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
std::vector< const ISyst * > getAllXsecNuTruthSysts_2017()
Get master XSec syst list for 2017 analyses (NuTruthSyst variant)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static std::unique_ptr< GenieMultiverseSpectra > LoadFrom(TDirectory *dir, const std::string &name)
Load the Spectrum structure of all universe from a TFile directory.
const int nuniverses
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
void NormalizeSpectra()
Normalize spectra. It has to be done AFTER SpectrumLoader::Go() is called.
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
Float_t w
Definition: plot.C:20
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
unsigned int fNUniverses
number of universes to generate
void FindBandBoundaries(BandOptions opt=kBandFromNominal) const
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string