pion_multiverse.C
Go to the documentation of this file.
1 // $ setup_nova -b maxopt -r R17-03-01-prod3reco.d; cd $WORK/test_releases/R17-03-01-prod3reco.d; srt_setup -a; setup_proxy
2 // $ setup_nova; cd $WORK/test_releases/development; srt_setup -a; setup_proxy; cd -
3 // $ setup_nova -b maxopt; cd $WORK/test_releases/development; srt_setup -a; setup_proxy; cd -
4 
5 #include "limits.h"
6 
7 #ifdef __CINT__
8 void pion_multiverse(int test = 1, int numUniverses = 300)
9 {
10  std::cout << "Sorry, you must run in compiled mode" << std::endl;
11 }
12 #else
13 
14 #include "CAFAna/Analysis/Plots.h"
15 #include "CAFAna/Core/HistAxis.h"
16 #include "CAFAna/Core/Spectrum.h"
20 #include "CAFAna/Systs/BeamSysts.h"
23 #include "CAFAna/Vars/XsecTunes.h"
30 //#include "./NDXSecTrackPID.cxx"
35 
36 
37 #include "TCanvas.h"
38 #include "TROOT.h"
39 
40 using namespace ana;
41 
42 void pion_multiverse(int test = 1, int numUniverses = 100)
43 {
44 
45  //Spectrum Loaders for MC
46  std::string tmpfname;
47  if(test==0) tmpfname = "dataset_def_name_newest_snapshot nd_fhc_cvnprongrespin_caf_minus_muonid_training";
48  else if (test==1) tmpfname = "dataset_def_name_newest_snapshot nd_fhc_cvnprongrespin_caf_minus_muonid_training_minus_fakedata with limit 1";
49  else if (test==2) tmpfname = "dataset_def_name_newest_snapshot nd_fhc_cvnprongrespin_caf_minus_muonid_training_minus_fakedata";
50  else tmpfname = "dataset_def_name_newest_snapshot nd_fhc_cvnprongrespin_caf_minus_muonid_training_minus_fakedata with limit 50";
51 
52  SpectrumLoader mcloader(tmpfname);
53 
54  HistAxis kScoreAxis("Pion Score", anglepionresbinning, kBestEventPionBDTScore);
55  HistAxis kMuonIDScoreAxis("MuonID", anglepionresbinning, kMuonID);
56  HistAxis kRecoEnergyAxis("Reco Pion Energy", energy2GeVmaxbinning, kRecoPionEnergy);
57  HistAxis kTruePionEnergyAxis("True Pion Energy", energy1point5GeVmaxbinning, kTruePionE);
58  HistAxis kRecoPionEnergyAxis("Reco Pion Energy", pionkebinning, kRecoPionEnergy);
59  HistAxis kRecoMuonEnergyAxis("Reco Pion Energy", mukebinsDisplay, kRecoMuKE);
60 
61  std::vector<const ISyst*> systs = getAllXsecSysts_2018();
62  GenieMultiverseParameters genie_shifts(numUniverses,systs);
63 
64  std::vector<SystShifts> genie_univ = genie_shifts.GetSystShifts();
65 
66  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
67 
68  Spectrum* genie_pionscore_spects_sel[numUniverses];
69  Spectrum* genie_pionscore_spects_sig[numUniverses];
70  Spectrum* genie_pionscore_spects_bkg[numUniverses];
71  Spectrum* genie_pionscore_spects_tru[numUniverses];
72 
73  Spectrum* genie_recoenergy_spects_sel[numUniverses];
74  Spectrum* genie_recoenergy_spects_sig[numUniverses];
75  Spectrum* genie_recoenergy_spects_bkg[numUniverses];
76  Spectrum* genie_recoenergy_spects_tru[numUniverses];
77 
78  Spectrum* genie_recopionenergy_spects_sel[numUniverses];
79  Spectrum* genie_recopionenergy_spects_sig[numUniverses];
80  Spectrum* genie_recopionenergy_spects_bkg[numUniverses];
81  Spectrum* genie_recopionenergy_spects_tru[numUniverses];
82 
83  Spectrum* genie_recomuonenergy_spects_sel[numUniverses];
84  Spectrum* genie_recomuonenergy_spects_sig[numUniverses];
85  Spectrum* genie_recomuonenergy_spects_bkg[numUniverses];
86  Spectrum* genie_recomuonenergy_spects_tru[numUniverses];
87 
88  Spectrum* genie_truepionenergy_spects_sel[numUniverses];
89  Spectrum* genie_truepionenergy_spects_sig[numUniverses];
90  Spectrum* genie_truepionenergy_spects_tru[numUniverses];
91 
92  Spectrum* genie_truepionenergyisreco_spects_sel[numUniverses];
93  Spectrum* genie_truepionenergyisreco_spects_sig[numUniverses];
94  Spectrum* genie_truepionenergyisreco_spects_tru[numUniverses];
95 
96  Spectrum* genie_template_spects_sel[numUniverses];
97  Spectrum* genie_template_spects_sig[numUniverses];
98  Spectrum* genie_template_spects_cc0pi[numUniverses];
99  Spectrum* genie_template_spects_otherbkg[numUniverses];
100 
101  Spectrum* ppfx_template_spects_sel[numUniverses];
102 
103 
104 
105  for(int k = 0; k < numUniverses; k++){
106  /*
107  genie_pionscore_spects_sel[k] = new Spectrum(mcloader,kScoreAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
108  genie_pionscore_spects_sig[k] = new Spectrum(mcloader,kScoreAxis,kPreselection&&PionSelec&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
109  genie_pionscore_spects_bkg[k] = new Spectrum(mcloader,kScoreAxis,kPreselection&&PionSelec&&!kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
110  genie_pionscore_spects_tru[k] = new Spectrum(mcloader,kScoreAxis,kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
111 
112  genie_recoenergy_spects_sel[k] = new Spectrum(mcloader,kRecoEnergyAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
113  genie_recoenergy_spects_sig[k] = new Spectrum(mcloader,kRecoEnergyAxis,kPreselection&&PionSelec&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
114  genie_recoenergy_spects_bkg[k] = new Spectrum(mcloader,kRecoEnergyAxis,kPreselection&&PionSelec&&!kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
115  genie_recoenergy_spects_tru[k] = new Spectrum(mcloader,kRecoEnergyAxis,kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
116 
117  genie_recopionenergy_spects_sel[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
118  genie_recopionenergy_spects_sig[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kPreselection&&PionSelec&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
119  genie_recopionenergy_spects_bkg[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kPreselection&&PionSelec&&!kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
120  genie_recopionenergy_spects_tru[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
121 
122  genie_recomuonenergy_spects_sel[k] = new Spectrum(mcloader,kRecoMuonEnergyAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
123  genie_recomuonenergy_spects_sig[k] = new Spectrum(mcloader,kRecoMuonEnergyAxis,kPreselection&&PionSelec&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
124  genie_recomuonenergy_spects_bkg[k] = new Spectrum(mcloader,kRecoMuonEnergyAxis,kPreselection&&PionSelec&&!kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
125  genie_recomuonenergy_spects_tru[k] = new Spectrum(mcloader,kRecoMuonEnergyAxis,kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);//*/
126 
127  genie_template_spects_sel[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kScoreAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
128  if(k<100)ppfx_template_spects_sel[k] = new Spectrum(mcloader,kRecoPionEnergyAxis,kScoreAxis,kPreselection&&PionSelec,kNoShift,kXSecCVWgt2018*ppfx_univ[k]);
129 
130  genie_truepionenergy_spects_sel[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kPreselection&&PionSelec,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
131  genie_truepionenergy_spects_sig[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kPreselection&&PionSelec&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
132  genie_truepionenergy_spects_tru[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
133 
134  genie_truepionenergyisreco_spects_sel[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kPreselection&&PionSelec&&kIsPionReco,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
135  genie_truepionenergyisreco_spects_sig[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kPreselection&&PionSelec&&kIsPionReco&&kTrueSignalEvent,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
136  genie_truepionenergyisreco_spects_tru[k] = new Spectrum(mcloader,kTruePionEnergyAxis,kTrueSignalEvent&&kIsPionReco,genie_univ[k],kXSecCVWgt2018*kPPFXFluxCVWgt);
137  //*/
138 
139  }
140 
141  // load spectra
142  mcloader.Go();
143 
144  string ofn;
145  ofn = "ccpi_fhc_multiverse_spectra.root";
146  TFile * outf = new TFile(ofn.c_str(), "recreate");
147 
148 
149  //.SaveTo(outf->mkdir(""));
150  for(int k = 0; k < numUniverses; k++){
151  char name[50];
152  /*
153 
154  sprintf(name, "mc_geniemv_pionscore_%i_sel",k);
155  genie_pionscore_spects_sel[k]->SaveTo(outf->mkdir(name));
156  sprintf(name, "mc_geniemv_pionscore_%i_sig",k);
157  genie_pionscore_spects_sig[k]->SaveTo(outf->mkdir(name));
158  sprintf(name, "mc_geniemv_pionscore_%i_bkg",k);
159  genie_pionscore_spects_bkg[k]->SaveTo(outf->mkdir(name));
160  sprintf(name, "mc_geniemv_pionscore_%i_tru",k);
161  genie_pionscore_spects_tru[k]->SaveTo(outf->mkdir(name));
162  ///////////////////////////////////////////////////////////
163  sprintf(name, "mc_geniemv_recoenergy_%i_sel",k);
164  genie_recoenergy_spects_sel[k]->SaveTo(outf->mkdir(name));
165  sprintf(name, "mc_geniemv_recoenergy_%i_sig",k);
166  genie_recoenergy_spects_sig[k]->SaveTo(outf->mkdir(name));
167  sprintf(name, "mc_geniemv_recoenergy_%i_bkg",k);
168  genie_recoenergy_spects_bkg[k]->SaveTo(outf->mkdir(name));
169  sprintf(name, "mc_geniemv_recoenergy_%i_tru",k);
170  genie_recoenergy_spects_tru[k]->SaveTo(outf->mkdir(name));
171  ///////////////////////////////////////////////////////////
172  sprintf(name, "mc_geniemv_recopionenergy_%i_sel",k);
173  genie_recopionenergy_spects_sel[k]->SaveTo(outf->mkdir(name));
174  sprintf(name, "mc_geniemv_recopionenergy_%i_sig",k);
175  genie_recopionenergy_spects_sig[k]->SaveTo(outf->mkdir(name));
176  sprintf(name, "mc_geniemv_recopionenergy_%i_bkg",k);
177  genie_recopionenergy_spects_bkg[k]->SaveTo(outf->mkdir(name));
178  sprintf(name, "mc_geniemv_recopionenergy_%i_tru",k);
179  genie_recopionenergy_spects_tru[k]->SaveTo(outf->mkdir(name));
180  ///////////////////////////////////////////////////////////
181  sprintf(name, "mc_geniemv_recomuonenergy_%i_sel",k);
182  genie_recomuonenergy_spects_sel[k]->SaveTo(outf->mkdir(name));
183  sprintf(name, "mc_geniemv_recomuonenergy_%i_sig",k);
184  genie_recomuonenergy_spects_sig[k]->SaveTo(outf->mkdir(name));
185  sprintf(name, "mc_geniemv_recomuonenergy_%i_bkg",k);
186  genie_recomuonenergy_spects_bkg[k]->SaveTo(outf->mkdir(name));
187  sprintf(name, "mc_geniemv_recomuonenergy_%i_tru",k);
188  genie_recomuonenergy_spects_tru[k]->SaveTo(outf->mkdir(name));
189  ///////////////////////////////////////////////////////////
190  //*/
191  sprintf(name, "mc_geniemv_template_%i_sel",k);
192  genie_template_spects_sel[k]->SaveTo(outf->mkdir(name));
193  sprintf(name, "mc_ppfxmv_template_%i_sel",k);
194  if(k<100)ppfx_template_spects_sel[k]->SaveTo(outf->mkdir(name));
195 
196  sprintf(name, "mc_geniemv_truepionenergy_%i_sel",k);
197  genie_truepionenergy_spects_sel[k]->SaveTo(outf->mkdir(name));
198  sprintf(name, "mc_geniemv_truepionenergy_%i_sig",k);
199  genie_truepionenergy_spects_sig[k]->SaveTo(outf->mkdir(name));
200  sprintf(name, "mc_geniemv_truepionenergy_%i_tru",k);
201  genie_truepionenergy_spects_tru[k]->SaveTo(outf->mkdir(name));
202 
203  sprintf(name, "mc_geniemv_truepionenergyisreco_%i_sel",k);
204  genie_truepionenergyisreco_spects_sel[k]->SaveTo(outf->mkdir(name));
205  sprintf(name, "mc_geniemv_truepionenergyisreco_%i_sig",k);
206  genie_truepionenergyisreco_spects_sig[k]->SaveTo(outf->mkdir(name));
207  sprintf(name, "mc_geniemv_truepionenergyisreco_%i_tru",k);
208  genie_truepionenergyisreco_spects_tru[k]->SaveTo(outf->mkdir(name));
209 
210  }
211 
212 
213  outf->Close();
214 }
215 #endif
const XML_Char * name
Definition: expat.h:151
const Var kRecoPionEnergy([](const caf::SRProxy *sr){float pionidx=GetBestPionTrack()(sr);float pionscore=GetBestPionID()(sr);if(pionscore< -1) return-2.0f;if(pionidx< 0) return-3.0f;float len=sr->trk.kalman.tracks[pionidx].len;float recoE=0.0699059+0.00454418 *len+-2.44087e-05 *len *len+1.19292e-07 *len *len *len+-2.61589e-10 *len *len *len *len+2.06403e-13 *len *len *len *len *len;return recoE;})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kRecoMuKE
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut PionSelec
const Var kBestEventPionBDTScore([](const caf::SRProxy *sr){float pionscore=GetBestPionID()(sr);if(pionscore<-1) return-2.0f;return pionscore;})
void pion_multiverse(int test=1, int numUniverses=100)
const Binning energy1point5GeVmaxbinning
Definition: NumuCC_CPiBin.h:36
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Binning anglepionresbinning
Definition: NumuCC_CPiBin.h:54
const Var kMuonID
Definition: NDXSecMuonPID.h:32
const Binning energy2GeVmaxbinning
Definition: NumuCC_CPiBin.h:32
const Cut kPreselection
const Var kTruePionE([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return-1.0;double trupionE=-1;int pionidx=-1;for(uint i=0;i< sr->mc.nu[0].prim.size();i++){if(abs(sr->mc.nu[0].prim[i].pdg)==211 &&sr->mc.nu[0].prim[i].p.T()>trupionE){pionidx=i;trupionE=sr->mc.nu[0].prim[i].p.T();}}if(pionidx==-1) return-1.0;double gamma=sr->mc.nu[0].prim[pionidx].p.Gamma();double mass=sqrt(sr->mc.nu[0].prim[pionidx].p.T()*sr->mc.nu[0].prim[pionidx].p.T()-sr->mc.nu[0].prim[pionidx].p.Mag()*sr->mc.nu[0].prim[pionidx].p.Mag());double KE=(gamma-1)*mass;return KE;})
TFile * outf
Definition: testXsec.C:51
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
const Binning mukebinsDisplay
Definition: NumuCCIncBins.h:91
std::vector< float > Spectrum
Definition: Constants.h:759
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< SystShifts > GetSystShifts()
const Cut kIsPionReco([](const caf::SRProxy *sr){if(sr->mc.nnu< 1) return false;double trupionE=-1;int pionidx=-1;for(uint i=0;i< sr->mc.nu[0].prim.size();i++){if(abs(sr->mc.nu[0].prim[i].pdg)==211 &&sr->mc.nu[0].prim[i].p.T()>trupionE){pionidx=i;trupionE=sr->mc.nu[0].prim[i].p.T();}}if(pionidx==-1) return false;if(sr->trk.kalman.ntracks< 1) return false;for(uint i=0;i< sr->trk.kalman.ntracks;i++){if(abs(sr->trk.kalman.tracks[i].truth.pdg)==211 &&sr->trk.kalman.tracks[i].truth.p.T()==trupionE) return true;}return false;})
size_t numUniverses
const Binning pionkebinning
const Cut kTrueSignalEvent
enum BeamMode string