TrivialCrossSectionAnalysis.cxx
Go to the documentation of this file.
4 #include "TObjString.h"
5 #include "CAFAna/Core/Ratio.h"
6 #include "CAFAna/XSec/Flux.h"
7 
8 #include "TH2.h"
9 #include "TH3.h"
10 #include "TVectorD.h"
11 
12 namespace ana
13 {
14  ///////////////////////////////////////////////////////////////////////
15  void
18  {
20  "TrivialCrossSectionAnalysis object not yet ready to setup systematic spectra." );
21 
22  fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
23  (syst,
24  new CrossSectionSpectra(*loader,
25  *fRecoAxis,
26  *fTruthAxis,
28  *fSignalCut,
29  new TrivialSignalEstimator(loader,
30  fRecoAxis,
32  fBkgdCuts,
33  *fSystShifts[syst],
34  *fWeights[syst]),
35  fFidMax,
36  fFidMin,
37  *fENuBinning,
38  fPDG,
39  *fWeights[syst],
40  *fNuTruthWeights[syst],
41  *fSystShifts[syst])));
42  }
43 
44  ///////////////////////////////////////////////////////////////////////
45  void
48  const Var * weight,
49  const NuTruthVar * weight_nt,
50  const SystShifts * shift)
51  {
52  SetWeight(syst, weight);
53  SetNuTruthWeight(syst, weight_nt);
54  SetSystShifts(syst, shift);
55  InitSystematicSpectra(syst, loader);
56  }
57 
58  ///////////////////////////////////////////////////////////////////////
59  void
61  std::vector<SystShifts> genie_shifts)
62  {
64  "CrossSectionAnalysis object not yet ready to setup genie multiverses." );
65  fGenieMVEffNum = new Multiverse(*loader,
68  genie_shifts,
70  fGenieMVEffDenom = new Multiverse(*loader,
71  *fTruthAxis,
72  *fSignalCut,
73  genie_shifts,
74  *fNuTruthWeights[XSecSysts::kNominal]);
75  fHasGenieMultiverses = true;
76  }
77 
78  ///////////////////////////////////////////////////////////////////////
79  void
81  std::vector<Var> ppfx_weights,
82  std::vector<NuTruthVar> ppfx_weights_nt)
83  {
85  "CrossSectionAnalysis object not yet ready to setup PPFX multiverses.");
86  std::vector<Spectrum *> mv_flux;
87  for(int i = 0; i < (int) ppfx_weights.size(); i++) {
88  mv_flux.push_back(DeriveFlux(*loader,
89  *fENuBinning,
90  fPDG,
91  &fFidMin,
92  &fFidMax,
94  ppfx_weights_nt[i]));
95  }
96  fPPFXMVFlux = new Multiverse(mv_flux);
97  fPPFXMVEffNum = new Multiverse(*loader,
101  ppfx_weights,
102  *fWeights[XSecSysts::kNominal]);
103  fPPFXMVEffDenom = new Multiverse(*loader,
104  *fTruthAxis,
105  *fSignalCut,
106  *fSystShifts[XSecSysts::kNominal],
107  ppfx_weights_nt,
108  *fNuTruthWeights[XSecSysts::kNominal]);
109 
110  // make sure we transferred ownership completely to the multiverse object
111  for(int i = 0; i < (int) ppfx_weights.size(); i++) {
112  mv_flux[i] = NULL;
113  }
114  fHasPPFXMultiverses = true;
115  }
116 
117  ///////////////////////////////////////////////////////////////////////
118  TH1 *
120  {
121  TH1 * num = ToHist(*fXSecSpectra[syst]->GetSelectedSignalRecoSpace(), fPOT);
122  TH1 * denom = ToHist(*fXSecSpectra[syst]->GetSelectedAllRecoSpace(), fPOT);
123  num->Divide(num, denom, 1, 1, "B"); // calculate binomial errors
124  return num;
125  }
126 
127  ///////////////////////////////////////////////////////////////////////
128  // basic constructor
130  const HistAxis * reco_axis,
131  const NuTruthCut & signal_cut,
132  const NuTruthHistAxis * truth_axis,
133  TVector3 fid_max,
134  TVector3 fid_min,
135  const Binning * enu_binning,
136  std::vector<Cut> bkgd_cuts,
137  int pdg,
138  bool is_differential)
139  : fBkgdCuts(bkgd_cuts)
140  {
141  this->fPDG = pdg;
142  this->fENuBinning = enu_binning;
143  this->fFidMax = fid_max;
144  this->fFidMin = fid_min;
145  this->fTarget = TargetCount(fFidMin, fFidMax, 1e6);
146  this->fSelectionCut = new Cut(selection_cut);
147  this->fSignalCut = new NuTruthCut(signal_cut);
148  this->fRecoAxis = reco_axis;
149  this->fTruthAxis = truth_axis;
150  this->fNDims = fRecoAxis->NDimensions();
151  this->fIsDifferential = is_differential;
152  }
153 
154  ///////////////////////////////////////////////////////////////////////
155  std::pair<TH1 *, TH1 *>
157  {
158  // load file systematic shifts
159  std::vector<XSecSysts::Syst_t> systs = XSecSysts::GetAllSystematics();
160  std::vector<TH1 *> shifts;
161  for(XSecSysts::Syst_t syst : systs) {
162  shifts.push_back(RelativeUncertainty(syst));
163  }
164  // genie and ppfx shifts
165  std::pair<TH1 *, TH1 *> flux_error = FluxError();
166  std::pair<TH1 *, TH1 *> genie_error = GenieError();
167  TH1 * flux_up = flux_error.first;
168  TH1 * flux_down = flux_error.second;
169  TH1 * genie_up = genie_error.first;
170  TH1 * genie_down = genie_error.second;
171 
172  TH1 * nominal = fXSec[XSecSysts::kNominal];
173  TH1 * errup = (TH1*) nominal->Clone();
174  TH1 * errdown = (TH1*) nominal->Clone();
175  double up, down, delta;
176 
177  for(int i = 1; i <= nominal->GetNbinsX(); i++) {
178  for(int j = 1; j <= nominal->GetNbinsY(); j++) {
179  for(int k = 1; k <= nominal->GetNbinsZ(); k++) {
180  double totalup = 0;
181  double totaldown = 0;
182  // process file systematics
183  for(int ishift = 0; ishift < (int) shifts.size(); ishift++) {
184  delta =
185  shifts[ishift]->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
186  if(delta > 0) {
187  up = delta;
188  down = 0;
189  }
190  else {
191  up = 0;
192  down = delta;
193  }
194  totalup += up * up;
195  totaldown += down * down;
196  }
197  // now genie and ppfx universes
198  up = genie_up->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
199  down = genie_down->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
200  totalup += up * up;
201  totaldown += down * down;
202 
203  up = flux_up->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
204  down = flux_down->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
205  totalup += up * up;
206  totaldown += down * down;
207  errup->SetBinContent(i, j, k, totalup);
208  errdown->SetBinContent(i, j, k, totaldown);
209  }
210  }
211  }
212  return std::pair<TH1 *, TH1 *> (errup, errdown);
213  }
214 
215  ///////////////////////////////////////////////////////////////////////
216  TH1 *
218  {
219  TH1 * rel = (TH1*) fXSec[XSecSysts::kNominal]->Clone();
220  double shift, cv, delta;
221  for(int i = 1; i <= rel->GetNbinsX(); i++) {
222  for(int j = 1; j <= rel->GetNbinsY(); j++) {
223  for(int k = 1; k <= rel->GetNbinsZ(); k++) {
224  cv = fXSec[XSecSysts::kNominal]->GetBinContent(i, j, k);
225  shift = fXSec[syst]->GetBinContent(i, j, k);
226  delta = fabs(shift - cv);
227  rel->SetBinContent(i, j, k, delta / cv);
228  }
229  }
230  }
231  return rel;
232  }
233 
234  ///////////////////////////////////////////////////////////////////////
235  std::pair<TH1 *, TH1 *>
237  {
238  // xsec is proportional to 1/efficiency so efficiency denom goes up in the
239  // numerator
240  Multiverse * genie_xsec = new Multiverse(*fGenieMVEffDenom);
241  genie_xsec->Divide(*fGenieMVEffNum, true);
242  *genie_xsec /= *Flux(XSecSysts::kNominal);
243  *genie_xsec *= *UnfoldedSignal(XSecSysts::kNominal);
244  return std::pair<TH1 *, TH1 *>
246  genie_xsec->GetMinusOneSigmaShift(fXSec[XSecSysts::kNominal]));
247  }
248 
249  ///////////////////////////////////////////////////////////////////////
250  std::pair<TH1 *, TH1 *>
252  {
253  // xsec is proportional to 1/efficiency so efficiency denom goes up in the
254  // numerator
255  Multiverse * ppfx_xsec = new Multiverse(*fPPFXMVEffDenom);
256  ppfx_xsec->Divide(*fPPFXMVEffNum, true);
257  ppfx_xsec->Divide(*fPPFXMVFlux);
258  *ppfx_xsec *= *UnfoldedSignal(XSecSysts::kNominal);
259 
260  return std::pair<TH1 *, TH1 *>
262  ppfx_xsec->GetMinusOneSigmaShift(fXSec[XSecSysts::kNominal]));
263  }
264 
265  ///////////////////////////////////////////////////////////////////////
266  TH1 *
268  {
269  return ToHist(*fXSecSpectra[syst]->GetFlux(), fPOT);
270  }
271 
272  ///////////////////////////////////////////////////////////////////////
273  TH1 *
275  {
276  TH1 * num = ToHist(*fXSecSpectra[syst]->GetSelectedSignalTrueSpace(), fPOT);
277  TH1 * denom = ToHist(*fXSecSpectra[syst]->GetTrueSignalTrueSpace(), fPOT);
278  num->Divide(num, denom, 1, 1, "B"); // calculate binomial errors
279  return num;
280  }
281 
282  ///////////////////////////////////////////////////////////////////////
283  std::unique_ptr<TrivialCrossSectionAnalysis>
285  {
286  dir = dir->GetDirectory(name.c_str()); // switch to subdir
287  assert(dir);
288 
289  TObjString * ptag = (TObjString*) dir->Get("type");
290  assert(ptag);
291 
292  const TString tag = ptag->GetString();
293  assert(tag == "TrivialCrossSectionAnalysis" && "Type does not match TrivialCrossSectionAnalysis");
294  delete ptag;
295 
296  Spectrum * data = Spectrum::LoadFrom(dir, "fData").release();
297  std::map<XSecSysts::Syst_t, CrossSectionSpectra *> xsec_spectra;
298  std::vector<XSecSysts::Syst_t> systs
300 
301  for(XSecSysts::Syst_t syst : systs) {
302  xsec_spectra.
303  insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
304  (syst,
306  }
307 
308  TVector3 * fidmax = (TVector3*) dir->Get("fFidMax");
309  TVector3 * fidmin = (TVector3*) dir->Get("fFidMin");
310  TVectorD * ndims = (TVectorD*) dir->Get("ndims");
311  Multiverse * genie_mv_eff_num = NULL;
312  Multiverse * genie_mv_eff_denom = NULL;
313  Multiverse * ppfx_mv_eff_num = NULL;
314  Multiverse * ppfx_mv_eff_denom= NULL;
315  Multiverse * ppfx_mv_flux = NULL;
316  bool isdifferential = false;
317  if(dir->GetDirectory("fIsDifferential"))
318  isdifferential = true;
319  if(dir->GetDirectory("fGenieMVEffNum")) {
320  genie_mv_eff_num =
321  Multiverse::LoadFrom(dir, "fGenieMVEffNum").release();
322  }
323  if(dir->GetDirectory("fGenieMVEffDenom")) {
324  genie_mv_eff_denom =
325  Multiverse::LoadFrom(dir, "fGenieMVEffDenom").release();
326  }
327  if(dir->GetDirectory("fPPFXMVFlux")) {
328  ppfx_mv_flux =
329  Multiverse::LoadFrom(dir, "fPPFXMVFlux").release();
330  }
331  if(dir->GetDirectory("fPPFXMVEffNum")) {
332  ppfx_mv_eff_num =
333  Multiverse::LoadFrom(dir, "fPPFXMVEffNum").release();
334  }
335  if(dir->GetDirectory("fPPFXMVEffDenom")) {
336  ppfx_mv_eff_denom =
337  Multiverse::LoadFrom(dir, "fPPFXMVEffDenom").release();
338  }
339 
340  std::unique_ptr<TrivialCrossSectionAnalysis> ret
341  ( new TrivialCrossSectionAnalysis(data,
342  xsec_spectra,
343  ppfx_mv_flux,
344  ppfx_mv_eff_num,
345  ppfx_mv_eff_denom,
346  genie_mv_eff_num,
347  genie_mv_eff_denom,
348  fidmax,
349  fidmin,
350  ndims,
351  isdifferential));
352  double POT = ret->fData->POT();
353  ret->SetPOT(POT);
354  for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : xsec_spectra)
355  if(xsec.second)
356  xsec.second = NULL;
357  data = NULL;
358  genie_mv_eff_num = NULL;
359  genie_mv_eff_denom = NULL;
360  ppfx_mv_flux = NULL;
361  ppfx_mv_eff_num = NULL;
362  ppfx_mv_eff_denom = NULL;
363  fidmax = NULL;
364  fidmin = NULL;
365  ndims = NULL;
366 
367  delete dir;
368 
369  return ret;
370  }
371 
372  ///////////////////////////////////////////////////////////////////////
373  void
375  {
376  TDirectory * tmp = gDirectory;
377 
378  dir = dir->mkdir(name.c_str()); // switch to subdir
379  dir->cd();
380 
381  TObjString("TrivialCrossSectionAnalysis").Write("type");
382  if(fIsDifferential)
383  TObjString("true").Write("fIsDifferential");
384  for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : fXSecSpectra) {
385  if(xsec.second)
386  xsec.second->SaveTo(dir, XSecSysts::LabelFromSystematic(xsec.first).Data());
387  }
388  fFidMax.Write("fFidMax");
389  fFidMin.Write("fFidMin");
390  TVectorD ndims(1);
391  ndims[0] = fNDims;
392  ndims.Write("ndims");
394  fGenieMVEffNum->SaveTo(dir, "fGenieMVEffNum");
395  fGenieMVEffDenom->SaveTo(dir, "fGenieMVEffDenom");
396  }
397  if(fHasPPFXMultiverses) {
398  fPPFXMVEffNum->SaveTo(dir, "fPPFXMVEffNum");
399  fPPFXMVEffDenom->SaveTo(dir, "fPPFXMVEffDenom");
400  fPPFXMVFlux->SaveTo(dir, "fPPFXMVFlux");
401  }
402  if(HasData())
403  fData->SaveTo(dir, "fData");
404 
405  dir->Write();
406  delete dir;
407 
408  tmp->cd();
409  }
410 
411  ///////////////////////////////////////////////////////////////////////
414  {
415  return new TrivialCrossSectionAnalysis(*this);
416  }
417 
418  ///////////////////////////////////////////////////////////////////////
419  // Make sure things get properly allocated during assignment
422  {
423  if ( this == &rhs )
424  return *this;
425  std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
426  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
427  {
428  this->fXSecSpectra
429  .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
430  (xsec.first,
431  new CrossSectionSpectra(*xsec.second)));
432  });
433 
434  std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
435  [this](std::pair<XSecSysts::Syst_t, const Var *> weight)
436  {
437  this->fWeights
438  .insert(std::pair<XSecSysts::Syst_t, const Var *>
439  (weight.first,
440  Copy(weight.second)));
441  });
442 
443  std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
444  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
445  {
446  this->fSystShifts
447  .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
448  (shift.first,
449  Copy(shift.second)));
450  });
451 
452  this->fTarget = rhs.fTarget;
453  this->fFidMax = rhs.fFidMax;
454  this->fFidMin = rhs.fFidMin;
455  this->fSelectionCut = Copy(rhs.fSelectionCut);
456  this->fRecoAxis = Copy(rhs.fRecoAxis);
457  this->fSignalCut = Copy(rhs.fSignalCut);
458  this->fTruthAxis = Copy(rhs.fTruthAxis);
459  this->fUnfoldMethod = rhs.fUnfoldMethod;
460  this->fUnfoldReg = rhs.fUnfoldReg;
461  this->fNDims = rhs.fNDims;
462  this->fIsDifferential = rhs.fIsDifferential;
463  if(rhs.fPPFXXSec.size() != 0) {
464  std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
465  [this](TH1 * ppfx)
466  {
467  this->fPPFXXSec.push_back((TH1*) ppfx->Clone());
468  });
469  }
470  if(rhs.fGenieXSec.size() != 0) {
471  std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
472  [this](TH1 * genie)
473  {
474  this->fGenieXSec.push_back((TH1*) genie->Clone());
475  });
476  }
477 
478  this->fGenieMVEffNum = new Multiverse(*rhs.fGenieMVEffNum);
479  this->fGenieMVEffDenom = new Multiverse(*rhs.fGenieMVEffDenom);
480  this->fPPFXMVFlux = new Multiverse(*rhs.fPPFXMVFlux);
481  this->fPPFXMVEffNum = new Multiverse(*rhs.fPPFXMVEffNum);
482  this->fPPFXMVEffDenom = new Multiverse(*rhs.fPPFXMVEffDenom);
483  this->fData = new Spectrum(*rhs.fData);
484  return *this;
485  }
486 
487  ///////////////////////////////////////////////////////////////////////
488  // Move assignment
491  {
492  if ( this == &rhs )
493  return *this;
494  std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
495  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
496  {
497  this->fXSecSpectra
498  .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
499  (xsec.first,
500  xsec.second));
501  xsec.second = NULL;
502  });
503 
504  std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
505  [this](std::pair<XSecSysts::Syst_t, const Var *> weight)
506  {
507  this->fWeights
508  .insert(std::pair<XSecSysts::Syst_t, const Var *>
509  (weight.first,
510  weight.second));
511  weight.second = NULL;
512  });
513 
514  std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
515  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
516  {
517  this->fSystShifts
518  .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
519  (shift.first,
520  shift.second));
521  shift.second = NULL;
522  });
523 
524  this->fTarget = rhs.fTarget;
525  this->fFidMax = rhs.fFidMax;
526  this->fFidMin = rhs.fFidMin;
527  this->fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
528  this->fRecoAxis = rhs.fRecoAxis; rhs.fRecoAxis = NULL;
529  this->fSignalCut = rhs.fSignalCut; rhs.fSignalCut = NULL;
530  this->fTruthAxis = rhs.fTruthAxis; rhs.fTruthAxis = NULL;
531  this->fUnfoldMethod = rhs.fUnfoldMethod;
532  this->fUnfoldReg = rhs.fUnfoldReg;
533  this->fNDims = rhs.fNDims;
534  this->fIsDifferential = rhs.fIsDifferential;
535  if(rhs.fPPFXXSec.size() != 0) {
536  std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
537  [this](TH1 * ppfx)
538  {
539  this->fPPFXXSec.push_back(ppfx);
540  ppfx = NULL;
541  });
542  }
543  if(rhs.fGenieXSec.size() != 0) {
544  std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
545  [this](TH1 * genie)
546  {
547  this->fGenieXSec.push_back(genie);
548  genie = NULL;
549  });
550  }
551 
552  this->fGenieMVEffNum = rhs.fGenieMVEffNum; rhs.fGenieMVEffNum = NULL;
553  this->fGenieMVEffDenom = rhs.fGenieMVEffDenom; rhs.fGenieMVEffDenom = NULL;
554  this->fPPFXMVFlux = rhs.fPPFXMVFlux; rhs.fPPFXMVFlux = NULL;
555  this->fPPFXMVEffNum = rhs.fPPFXMVEffNum; rhs.fPPFXMVEffNum = NULL;
556  this->fPPFXMVEffDenom = rhs.fPPFXMVEffDenom; rhs.fPPFXMVEffDenom = NULL;
557  this->fData = rhs.fData; rhs.fData = NULL;
558  return *this;
559  }
560 
561  ///////////////////////////////////////////////////////////////////////
562  // Move constructor
564  {
565  std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
566  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
567  {
568  this->fXSecSpectra
569  .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
570  (xsec.first,
571  xsec.second));
572  xsec.second = NULL;
573  });
574 
575  std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
576  [this](std::pair<XSecSysts::Syst_t, const Var *> weight)
577  {
578  this->fWeights
579  .insert(std::pair<XSecSysts::Syst_t, const Var *>
580  (weight.first,
581  weight.second));
582  weight.second = NULL;
583  });
584 
585  std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
586  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
587  {
588  this->fSystShifts
589  .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
590  (shift.first,
591  shift.second));
592  shift.second = NULL;
593  });
594 
595  this->fTarget = rhs.fTarget;
596  this->fFidMax = rhs.fFidMax;
597  this->fFidMin = rhs.fFidMin;
598  this->fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
599  this->fRecoAxis = rhs.fRecoAxis; rhs.fRecoAxis = NULL;
600  this->fSignalCut = rhs.fSignalCut; rhs.fSignalCut = NULL;
601  this->fTruthAxis = rhs.fTruthAxis; rhs.fTruthAxis = NULL;
602  this->fUnfoldMethod = rhs.fUnfoldMethod;
603  this->fUnfoldReg = rhs.fUnfoldReg;
604  this->fNDims = rhs.fNDims;
605  this->fIsDifferential = rhs.fIsDifferential;
606  if(rhs.fPPFXXSec.size() != 0) {
607  std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
608  [this](TH1 * ppfx)
609  {
610  this->fPPFXXSec.push_back(ppfx);
611  ppfx = NULL;
612  });
613  }
614  if(rhs.fGenieXSec.size() != 0) {
615  std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
616  [this](TH1 * genie)
617  {
618  this->fGenieXSec.push_back(genie);
619  genie = NULL;
620  });
621  }
622 
623  this->fGenieMVEffNum = rhs.fGenieMVEffNum; rhs.fGenieMVEffNum = NULL;
624  this->fGenieMVEffDenom = rhs.fGenieMVEffDenom; rhs.fGenieMVEffDenom = NULL;
625  this->fPPFXMVFlux = rhs.fPPFXMVFlux; rhs.fPPFXMVFlux = NULL;
626  this->fPPFXMVEffNum = rhs.fPPFXMVEffNum; rhs.fPPFXMVEffNum = NULL;
627  this->fPPFXMVEffDenom = rhs.fPPFXMVEffDenom; rhs.fPPFXMVEffDenom = NULL;
628  this->fData = rhs.fData; rhs.fData = NULL;
629  }
630 
631  ///////////////////////////////////////////////////////////////////////
632  // Copy constructor
634  {
635  std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
636  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
637  {
638  this->fXSecSpectra
639  .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
640  (xsec.first,
641  new CrossSectionSpectra(*xsec.second)));
642  });
643 
644  std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
645  [this](std::pair<XSecSysts::Syst_t, const Var *> weight)
646  {
647  this->fWeights
648  .insert(std::pair<XSecSysts::Syst_t, const Var *>
649  (weight.first,
650  Copy(weight.second)));
651  });
652 
653  std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
654  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift)
655  {
656  this->fSystShifts
657  .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
658  (shift.first,
659  Copy(shift.second)));
660  });
661 
662  this->fTarget = rhs.fTarget;
663  this->fFidMax = rhs.fFidMax;
664  this->fFidMin = rhs.fFidMin;
665  this->fSelectionCut = Copy(rhs.fSelectionCut);
666  this->fSignalCut = Copy(rhs.fSignalCut);
667  this->fRecoAxis = Copy(rhs.fRecoAxis);
668  this->fTruthAxis = Copy(rhs.fTruthAxis);
669  this->fUnfoldMethod = rhs.fUnfoldMethod;
670  this->fUnfoldReg = rhs.fUnfoldReg;
671  this->fNDims = rhs.fNDims;
672  this->fIsDifferential = rhs.fIsDifferential;
673  if(rhs.fPPFXXSec.size() != 0) {
674  std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
675  [this](TH1 * ppfx)
676  {
677  this->fPPFXXSec.push_back((TH1*) ppfx->Clone());
678  });
679  }
680  if(rhs.fGenieXSec.size() != 0) {
681  std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
682  [this](TH1 * genie)
683  {
684  this->fGenieXSec.push_back((TH1*) genie->Clone());
685  });
686  }
687  this->fGenieMVEffNum = new Multiverse(*rhs.fGenieMVEffNum);
688  this->fGenieMVEffDenom = new Multiverse(*rhs.fGenieMVEffDenom);
689  this->fPPFXMVFlux = new Multiverse(*rhs.fPPFXMVFlux);
690  this->fPPFXMVEffNum = new Multiverse(*rhs.fPPFXMVEffNum);
691  this->fPPFXMVEffDenom = new Multiverse(*rhs.fPPFXMVEffDenom);
692  this->fData = new Spectrum(*rhs.fData);
693  }
694 
695 
696  ///////////////////////////////////////////////////////////////////////
698  std::map<XSecSysts::Syst_t, CrossSectionSpectra*> _fXSecSpectra,
699  Multiverse * ppfx_mv_flux,
700  Multiverse * ppfx_mv_eff_num,
701  Multiverse * ppfx_mv_eff_denom,
702  Multiverse * genie_mv_eff_num,
703  Multiverse * genie_mv_eff_denom,
704  TVector3 * _fFidMax,
705  TVector3 * _fFidMin,
706  TVectorD * ndims,
707  bool isdifferential)
708  {
709  this->fData = data;
710  std::for_each(_fXSecSpectra.begin(), _fXSecSpectra.end(),
711  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec)
712  {
713  this->fXSecSpectra
714  .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
715  (xsec.first,
716  xsec.second));
717  });
718  this->fFidMax = *_fFidMax;
719  this->fFidMin = *_fFidMin;
720  this->fNDims = (*ndims)[0];
721  this->fIsDifferential = isdifferential;
722  this->fGenieMVEffNum = genie_mv_eff_num;
723  this->fGenieMVEffDenom = genie_mv_eff_denom;
724  this->fPPFXMVFlux = ppfx_mv_flux;
725  this->fPPFXMVEffNum = ppfx_mv_eff_num;
726  this->fPPFXMVEffDenom = ppfx_mv_eff_denom;
727 
728  for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : _fXSecSpectra)
729  if(xsec.second != NULL)
730  xsec.second = NULL;
731  genie_mv_eff_num = NULL;
732  genie_mv_eff_denom = NULL;
733  ppfx_mv_flux = NULL;
734  ppfx_mv_eff_num = NULL;
735  ppfx_mv_eff_denom = NULL;
736  }
737 
738  ///////////////////////////////////////////////////////////////////////
740  {
741  for(int ipxsec = 0; ipxsec < (int) fPPFXXSec.size(); ipxsec++) {
742  delete fPPFXXSec[ipxsec];
743  }
744  for(int igxsec = 0; igxsec < (int) fGenieXSec.size(); igxsec++) {
745  delete fGenieXSec[igxsec];
746  }
747  delete fData;
748  delete fGenieMVEffNum;
749  delete fGenieMVEffDenom;
750  delete fPPFXMVEffNum;
751  delete fPPFXMVEffDenom;
752  delete fPPFXMVFlux;
753  }
754 }
755 
static TString LabelFromSystematic(Syst_t)
Definition: XSecSysts.cxx:25
const XML_Char * name
Definition: expat.h:151
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
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
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int NDimensions() const
Definition: HistAxis.h:87
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:443
TrivialCrossSectionAnalysis & operator=(const TrivialCrossSectionAnalysis &)
double delta
Definition: runWimpSim.h:98
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
GenericCut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
TH1 * Flux(XSecSysts::Syst_t) override
Float_t tmp
Definition: plot.C:36
GenericCut< caf::SRNeutrinoProxy > NuTruthCut
Cut designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Cut.h:104
ICrossSectionAnalysis * Clone() const override
TH1 * Efficiency(XSecSysts::Syst_t) override
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
Definition: Multiverse.cxx:450
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const TH1 * UnfoldedSignal(XSecSysts::Syst_t)
TH1 * RelativeUncertainty(XSecSysts::Syst_t) override
Function that calculates relative uncertainty on the cross section result.
GFluxI * GetFlux(void)
Definition: gAtmoEvGen.cxx:441
const XML_Char const XML_Char * data
Definition: expat.h:268
std::pair< TH1 *, TH1 * > GenieError()
Function that calculates the up (first) and down (second) error on the genie universes.
static std::unique_ptr< CrossSectionSpectra > LoadFrom(TDirectory *, const std::string &)
void SetWeight(XSecSysts::Syst_t, const Var *)
std::pair< TH1 *, TH1 * > TotalErrors() override
TrivialCrossSectionAnalysis()
default constructor. Do nothing
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
std::map< XSecSysts::Syst_t, CrossSectionSpectra * > fXSecSpectra
const std::string cv[Ncv]
void SaveTo(TDirectory *, const std::string &name) const override
static std::vector< Syst_t > GetAllSystematics()
Definition: XSecSysts.cxx:9
std::pair< TH1 *, TH1 * > FluxError() override
Function that calculates the up (first) and down (second) error on the flux estimation.
void SetSystShifts(XSecSysts::Syst_t, const SystShifts *)
static std::unique_ptr< TrivialCrossSectionAnalysis > LoadFrom(TDirectory *, const std::string &name)
std::map< XSecSysts::Syst_t, const NuTruthVar * > fNuTruthWeights
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
Definition: HistAxis.cxx:95
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
const double j
Definition: BetheBloch.cxx:29
string rel
Definition: shutoffs.py:11
void InitPPFXUniverses(SpectrumLoaderBase *, std::vector< Var >, std::vector< NuTruthVar >)
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:527
void SetNuTruthWeight(XSecSysts::Syst_t, const NuTruthVar *)
Double_t xsec[nknots]
Definition: testXsec.C:47
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Multiverse.cxx:573
std::vector< double > POT
Base class for the various types of spectrum loader.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::map< XSecSysts::Syst_t, const Var * > fWeights
TH1 * Purity(XSecSysts::Syst_t) override
Function that returns the purity of the specified sample.
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:56
TDirectory * dir
Definition: macro.C:5
void InitGenieUniverses(SpectrumLoaderBase *, std::vector< SystShifts >)
int num
Definition: f2_nu.C:119
TGraph * genie
Definition: Xsec_final.C:116
std::map< XSecSysts::Syst_t, const SystShifts * > fSystShifts
assert(nhit_max >=nhit_nbins)
const NuTruthHistAxis * fTruthAxis
std::map< XSecSysts::Syst_t, TH1 * > fXSec
TH1 * ToHist(const Spectrum &, double=-1)
Count number of targets within the main part of the ND (vSuperBlockND)
Definition: TargetCount.h:10
void InitSystematicSpectra(XSecSysts::Syst_t, SpectrumLoaderBase *)
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:128
void Divide(Multiverse &, bool=false)
Definition: Multiverse.cxx:139