genie_syst_make.C
Go to the documentation of this file.
1 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Systs/Systs.h"
12 
14 
15 #include "TFile.h"
16 
17 #include <iostream>
18 
19 using namespace ana;
20 
21 void genie_syst_make(const unsigned int nuniverses = 1000)
22 {
23  SpectrumLoader loaderNDFHC("prod_caf_R19-11-18-prod5reco.d.h.l_nd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1");
24  SpectrumLoader loaderFDFHC("prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_fhc_nova_v08_full_v1");
25  SpectrumLoader loaderNDRHC("prod_caf_R19-11-18-prod5reco.d_nd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1");
26  SpectrumLoader loaderFDRHC("prod_caf_R19-11-18-prod5reco.f_fd_genie_N1810j0211a_nonswap_rhc_nova_v08_full_v1");
27 
28  loaderNDFHC.SetSpillCut(kStandardDQCuts);
29  loaderNDRHC.SetSpillCut(kStandardDQCuts);
30  loaderFDFHC.SetSpillCut(kStandardDQCuts);
31  loaderFDRHC.SetSpillCut(kStandardDQCuts);
32 
33  struct AnaDef{
35  Cut anaCut;
36  bool isND;
37  };
38 
39  struct cutDef{
41  Cut cut;
42  };
43 
44  const std::vector<AnaDef> anas = {{"NumuND",kNumu2020ND,1},
45  {"NueND",kNue2020ND,1},
46  {"NumuFD",kNumu2020FD,0},
47  {"NueFD",kNue2020FD,0}};
48 
49  const std::vector<cutDef> nuCuts = {{"nues",kIsNue && !kIsNC},
50  {"numus",kIsNumu && !kIsNC},
51  {"ncs",kIsNC}};
52 
53  const std::vector<cutDef> qeCuts = {{"All",kNoCut},
54  {"QE",kIsQE},
55  {"nonQE",!kIsQE}};
56 
57  const Binning bins = Binning::Simple(20,0,5);
58  const HistAxis axis("True E (GeV)",bins,kTrueE);
59 
61 
62  //auto allsysts = getAllXsecSysts_2018();
63  //auto allsysts = getAna2018SmallXsecSystsForPCA();
64  //auto allsysts = getAllXsecSysts_2020();
65  auto allsysts = get3FlavorAna2020SmallXsecSysts();
66  const int totalsysts = allsysts.size();
67 
68  std::vector <std::pair <const ISyst*, SystMode> > systConfigs;
69  for (int idx = 0; idx < totalsysts; ++idx)
70  {
71  // default seed is 1001
72  systConfigs.push_back({allsysts[idx], kSystGaussian});
73  }
74 
75  auto multiverse_shifts = GetSystShiftsMultiverse(systConfigs, nuniverses);
76  assert(multiverse_shifts.size() == nuniverses);
77 
78  Spectrum* multiv[anas.size()][2][nuCuts.size()][qeCuts.size()][nuniverses];
79  for(unsigned int i = 0;i < anas.size();++i){
80  for(unsigned int j = 0;j < nuCuts.size();++j){
81  for(unsigned int k = 0;k < qeCuts.size();++k){
82  for(unsigned int l = 0;l < nuniverses;++l){
83  // 0 is the nominal universe
84  if(l == 0){
85  multiv[i][0][j][k][l] = new Spectrum(anas[i].isND ? loaderNDFHC : loaderFDFHC,
86  axis,anas[i].anaCut && nuCuts[j].cut && qeCuts[k].cut,
87  kNoShift,weight);
88  multiv[i][1][j][k][l] = new Spectrum(anas[i].isND ? loaderNDRHC : loaderFDRHC,
89  axis,anas[i].anaCut && nuCuts[j].cut && qeCuts[k].cut,
90  kNoShift,weight);
91  }
92  else{
93  multiv[i][0][j][k][l] = new Spectrum(anas[i].isND ? loaderNDFHC : loaderFDFHC,
94  axis,anas[i].anaCut && nuCuts[j].cut && qeCuts[k].cut,
95  multiverse_shifts[l],weight);
96  multiv[i][1][j][k][l] = new Spectrum(anas[i].isND ? loaderNDRHC : loaderFDRHC,
97  axis,anas[i].anaCut && nuCuts[j].cut && qeCuts[k].cut,
98  multiverse_shifts[l],weight);
99  }
100  }
101  }
102  }
103  }
104 
105  loaderNDFHC.Go();
106  loaderFDFHC.Go();
107  loaderNDRHC.Go();
108  loaderFDRHC.Go();
109 
110  TFile fout("make_genieverse.root","RECREATE");
111  for(unsigned int i = 0;i < anas.size();++i){
112  for(unsigned int j = 0;j < nuCuts.size();++j){
113  for(unsigned int k = 0;k < qeCuts.size();++k){
114  TDirectory* fhcdir = fout.mkdir(("genieverse_FHC_"+anas[i].name+"_"+nuCuts[j].name+"_"+qeCuts[k].name).c_str());
115  fhcdir->cd();
116  for(unsigned int l = 0;l < nuniverses;++l)
117  multiv[i][0][j][k][l]->SaveTo(fhcdir, ("Univ"+std::to_string(l)).c_str());
118 
119  TDirectory* rhcdir = fout.mkdir(("genieverse_RHC_"+anas[i].name+"_"+nuCuts[j].name+"_"+qeCuts[k].name).c_str());
120  rhcdir->cd();
121  for(unsigned int l = 0;l < nuniverses;++l)
122  multiv[i][1][j][k][l]->SaveTo(rhcdir, ("Univ"+std::to_string(l)).c_str());
123  }
124  }
125  }
126 
127 }
const Cut kIsQE
Definition: TruthCuts.cxx:104
std::vector< const ISyst * > get3FlavorAna2020SmallXsecSysts(const EAnaType2020 ana)
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(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
const Cut kIsNumu([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==14);})
Definition: TruthCuts.h:10
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
void SetSpillCut(const SpillCut &cut)
const Cut kNue2020ND
Definition: NueCuts2020.h:178
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
void genie_syst_make(const unsigned int nuniverses=1000)
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
const Cut kIsNue([](const caf::SRProxy *sr){return(sr->mc.nnu==1 &&abs(sr->mc.nu[0].pdg)==12);})
Definition: TruthCuts.h:9
virtual void Go() override
Load all the registered spectra.
const double j
Definition: BetheBloch.cxx:29
std::vector< SystShifts > GetSystShiftsMultiverse(const std::vector< std::pair< const ISyst *, SystMode > > &systConfigs, const int nUniverses, const int seed)
Definition: SystShifts.cxx:283
const Cut kNumu2020FD
Definition: NumuCuts2020.h:59
std::vector< float > Spectrum
Definition: Constants.h:759
const Cut kNue2020FD
Definition: NueCuts2020.h:65
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Binning bins
Definition: NumuCC_CPiBin.h:8
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
assert(nhit_max >=nhit_nbins)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const int nuniverses
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const Var kXSecCVWgt2020
Definition: XsecTunes.h:105
enum BeamMode string