getCrossSectionAnalysisSpectra.C
Go to the documentation of this file.
1 #include "CAFAna/Core/SpectrumLoader.h"
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
4 #include "CAFAna/Core/SystShifts.h"
5 #include "CAFAna/Cuts/SpillCuts.h"
6 #include "CAFAna/Cuts/TruthCuts.h"
7 #include "CAFAna/Systs/XSecSystLists.h"
8 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
9 #include "3FlavorAna/Vars/NueVars.h"
10 #include "3FlavorAna/Vars/NueVarsExtra.h"
11 #include "CAFAna/Vars/GenieWeights.h"
12 #include "CAFAna/Vars/PPFXWeights.h"
13 #include "CAFAna/XSec/GenieMultiverseSyst.h"
14 #include "CAFAna/XSec/FluxMultiverseSyst.h"
15 #include "CAFAna/Vars/XsecTunes.h"
16 
17 #include "NDAna/nuecc_inc/NueCCIncVars.h"
18 #include "NDAna/nuecc_inc/NueCCIncCuts.h"
19 #include "NDAna/nuecc_inc/NueCCIncExtra.h"
20 #include "NDAna/nuecc_inc/NueCCIncCrossSectionAnalysis.h"
21 
22 using namespace ana;
23 
24 const Var kcvn ([](const caf::SRProxy* sr){
25  return (sr->sel.cvn.nueid);
26  });
27 
28 const Var kcvn2017([] (const caf::SRProxy* sr){
29  return (sr->sel.cvn2017.nueid);
30  });
31 
32 const Var kcvnProd3([] (const caf::SRProxy* sr){
33  return (sr->sel.cvnProd3Train.nueid);
34  });
35 
36 const Var klid([](const caf::SRProxy* sr){
37  return sr->sel.lid.ann;
38  });
39 
40 //sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.photonid
41 const Var kprongCVN([](const caf::SRProxy* sr){
42  if (sr->vtx.nelastic < 1) return -1000.f;
43  if (sr->vtx.elastic[0].fuzzyk.npng < 1) return -1000.f;
44 
45  int nProngs = sr->vtx.elastic[0].fuzzyk.npng;
46  float maxProngCVNe = -5.;
47  for(int iProng = 0; iProng < nProngs; iProng++){
48  if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >
49  maxProngCVNe) {
50  maxProngCVNe =
51  sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;
52  }
53  }
54  return maxProngCVNe;
55  });
56 
57 const int kNumPIDVars = 1;
58 const mHistAxisDef pid_vars[kNumPIDVars] = {
59  {"cvne2017", {"CVNe 2017", pidbins, kcvn2017}},
60  //{"cvneProd3", {"CVNe Prod 3", pidbins, kcvnProd3}},
61  //{"lid", {"LID", pidbins, klid}},
62  //{"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
63  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
64  // pidbins, kflatfluxCVN}},
65 };
66 
67 
68 void getCrossSectionAnalysisSpectra()
69 {
70  SpectrumLoader lNDMC(NominalMC_full_prod4_CVNProng);
71  SpectrumLoader lNDLightDwn(DWNLightFiles_full_prod4_CVNProng);
72  SpectrumLoader lNDLightUp(UPLightFiles_full_prod4_CVNProng);
73  SpectrumLoader lNDCkv(CherenkovFiles_full_prod4_CVNProng);
74  SpectrumLoader lNDCalibDwn(DWNCalibFiles_full_prod4_CVNProng);
75  SpectrumLoader lNDCalibUp(UPCalibFiles_full_prod4_CVNProng);
76  SpectrumLoader lNDCalibShape(ShapeCalibFiles_full_prod4_CVNProng);
77 
78  lNDMC.SetSpillCut(kStandardSpillCuts);
79  lNDLightDwn.SetSpillCut(kStandardSpillCuts);
80  lNDLightUp.SetSpillCut(kStandardSpillCuts);
81  lNDCalibDwn.SetSpillCut(kStandardSpillCuts);
82  lNDCalibUp.SetSpillCut(kStandardSpillCuts);
83  lNDCkv.SetSpillCut(kStandardSpillCuts);
84  lNDCalibShape.SetSpillCut(kStandardSpillCuts);
85 
86 
87  std::vector<Cut> template_cuts;
88  for(int i = 0; i < kcNumChns; i++)
89  template_cuts.push_back(chns[i].cut);
90 
91 
92  //Now For Analysis Spectra
93  NueCCIncCrossSectionAnalysis nominal_ana(lNDMC,
94  analysis_vars[2].axis,
95  analysis_vars[1].axis,
96  analysis_vars[0].axis,
97  kRecoElecEVsCosStandardAxis,
98  kTrueElecEVsCosStandardAxisST,
99  analysis_vars_truth[0].axis,
100  pid_vars[0].axis,
101  enubins,
102  kcPresel,
103  kIsTrueSigST,
104  vtxmin,
105  vtxmax,
106  kNoShift,
107  kPPFXFluxCVWgtST *
108  kXSecCVWgt2018_smallerDISScale_NT,
109  kPPFXFluxCVWgt * kXSecCVWgt2018);
110 
111  NueCCIncCrossSectionTemplates nominal_templates(lNDMC,
112  analysis_vars[2].axis,
113  analysis_vars[1].axis,
114  analysis_vars[0].axis,
115  pid_vars[0].axis,
116  kcPresel,
117  template_cuts,
118  kcNumChns,
119  kNoShift,
120  kPPFXFluxCVWgt*kXSecCVWgt2018);
121 
122 
123 
124  NueCCIncCrossSectionAnalysis light_dwn_ana(lNDLightDwn,
125  analysis_vars[2].axis,
126  analysis_vars[1].axis,
127  analysis_vars[0].axis,
128  kRecoElecEVsCosStandardAxis,
129  kTrueElecEVsCosStandardAxisST,
130  analysis_vars_truth[0].axis,
131  pid_vars[0].axis,
132  enubins,
133  kcPresel,
134  kIsTrueSigST,
135  vtxmin,
136  vtxmax,
137  kNoShift,
138  kPPFXFluxCVWgtST *
139  kXSecCVWgt2018_smallerDISScale_NT,
140  kPPFXFluxCVWgt * kXSecCVWgt2018);
141 
142  NueCCIncCrossSectionAnalysis light_up_ana(lNDLightUp,
143  analysis_vars[2].axis,
144  analysis_vars[1].axis,
145  analysis_vars[0].axis,
146  kRecoElecEVsCosStandardAxis,
147  kTrueElecEVsCosStandardAxisST,
148  analysis_vars_truth[0].axis,
149  pid_vars[0].axis,
150  enubins,
151  kcPresel,
152  kIsTrueSigST,
153  vtxmin,
154  vtxmax,
155  kNoShift,
156  kPPFXFluxCVWgtST *
157  kXSecCVWgt2018_smallerDISScale_NT,
158  kPPFXFluxCVWgt * kXSecCVWgt2018);
159 
160  NueCCIncCrossSectionAnalysis calib_dwn_ana(lNDCalibDwn,
161  analysis_vars[2].axis,
162  analysis_vars[1].axis,
163  analysis_vars[0].axis,
164  kRecoElecEVsCosStandardAxis,
165  kTrueElecEVsCosStandardAxisST,
166  analysis_vars_truth[0].axis,
167  pid_vars[0].axis,
168  enubins,
169  kcPresel,
170  kIsTrueSigST,
171  vtxmin,
172  vtxmax,
173  kNoShift,
174  kPPFXFluxCVWgtST *
175  kXSecCVWgt2018_smallerDISScale_NT,
176  kPPFXFluxCVWgt * kXSecCVWgt2018);
177 
178  NueCCIncCrossSectionAnalysis calib_up_ana(lNDCalibUp,
179  analysis_vars[2].axis,
180  analysis_vars[1].axis,
181  analysis_vars[0].axis,
182  kRecoElecEVsCosStandardAxis,
183  kTrueElecEVsCosStandardAxisST,
184  analysis_vars_truth[0].axis,
185  pid_vars[0].axis,
186  enubins,
187  kcPresel,
188  kIsTrueSigST,
189  vtxmin,
190  vtxmax,
191  kNoShift,
192  kPPFXFluxCVWgtST *
193  kXSecCVWgt2018_smallerDISScale_NT,
194  kPPFXFluxCVWgt * kXSecCVWgt2018);
195 
196  NueCCIncCrossSectionAnalysis ckv_ana(lNDCkv,
197  analysis_vars[2].axis,
198  analysis_vars[1].axis,
199  analysis_vars[0].axis,
200  kRecoElecEVsCosStandardAxis,
201  kTrueElecEVsCosStandardAxisST,
202  analysis_vars_truth[0].axis,
203  pid_vars[0].axis,
204  enubins,
205  kcPresel,
206  kIsTrueSigST,
207  vtxmin,
208  vtxmax,
209  kNoShift,
210  kPPFXFluxCVWgtST *
211  kXSecCVWgt2018_smallerDISScale_NT,
212  kPPFXFluxCVWgt * kXSecCVWgt2018);
213 
214  NueCCIncCrossSectionAnalysis
215  calib_shape_ana(lNDCalibShape,
216  analysis_vars[2].axis,
217  analysis_vars[1].axis,
218  analysis_vars[0].axis,
219  kRecoElecEVsCosStandardAxis,
220  kTrueElecEVsCosStandardAxisST,
221  analysis_vars_truth[0].axis,
222  pid_vars[0].axis,
223  enubins,
224  kcPresel,
225  kIsTrueSigST,
226  vtxmin,
227  vtxmax,
228  kNoShift,
229  kPPFXFluxCVWgtST *
230  kXSecCVWgt2018_smallerDISScale_NT,
231  kPPFXFluxCVWgt * kXSecCVWgt2018);
232 
233 
234  std::vector<NueCCIncCrossSectionAnalysis*> genie_universes_ana;
235  std::vector<NueCCIncCrossSectionAnalysis*> ppfx_universes_ana;
236 
237  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
238  std::vector<ana::NuTruthVar> ppfx_univ_st = ana::GetkPPFXFluxUnivWgtST();
239  GenieMultiverseParameters verse(100,ana::getAllXsecSysts_2018());
240  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
241 
242  for(int i = 0; i < 100; i++){
243 
244  genie_universes_ana.push_back(new NueCCIncCrossSectionAnalysis
245  (lNDMC,
246  analysis_vars[2].axis,
247  analysis_vars[1].axis,
248  analysis_vars[0].axis,
249  kRecoElecEVsCosStandardAxis,
250  kTrueElecEVsCosStandardAxisST,
251  analysis_vars_truth[0].axis,
252  pid_vars[0].axis,
253  enubins,
254  kcPresel,
255  kIsTrueSigST,
256  vtxmin,
257  vtxmax,
258  genie_shifts[i],
259  kPPFXFluxCVWgtST *
260  kXSecCVWgt2018_smallerDISScale_NT,
261  kPPFXFluxCVWgt * kXSecCVWgt2018));
262 
263  ppfx_universes_ana.push_back(new NueCCIncCrossSectionAnalysis
264  (lNDMC,
265  analysis_vars[2].axis,
266  analysis_vars[1].axis,
267  analysis_vars[0].axis,
268  kRecoElecEVsCosStandardAxis,
269  kTrueElecEVsCosStandardAxisST,
270  analysis_vars_truth[0].axis,
271  pid_vars[0].axis,
272  enubins,
273  kcPresel,
274  kIsTrueSigST,
275  vtxmin,
276  vtxmax,
277  kNoShift,
278  ppfx_univ_st[i] *
279  kXSecCVWgt2018_smallerDISScale_NT,
280  ppfx_univ[i] * kXSecCVWgt2018));
281  }
282 
283 
284  lNDMC.Go();
285  lNDLightDwn.Go();
286  lNDLightUp.Go();
287  lNDCalibDwn.Go();
288  lNDCalibUp.Go();
289  lNDCkv.Go();
290  lNDCalibShape.Go();
291 
292  //Save Files
293  TFile* fout = new TFile("fAnalysisSpectra.root", "recreate");
294  fout->cd();
295 
296  TDirectory* nom_dir = fout->mkdir(Form("%s","Nominal"));
297  nom_dir->cd();
298  nominal_ana.SaveTo(nom_dir->mkdir("Analysis"));
299  nominal_templates.SaveTo(nom_dir->mkdir("Templates"));
300 
301  TDirectory* light_dwn_dir = fout->mkdir(Form("%s","LightDwn"));
302  light_dwn_dir->cd();
303  light_dwn_ana.SaveTo(light_dwn_dir->mkdir("Analysis"));
304 
305  TDirectory* light_up_dir = fout->mkdir(Form("%s","LightUp"));
306  light_up_dir->cd();
307  light_up_ana.SaveTo(light_up_dir->mkdir("Analysis"));
308 
309  TDirectory* calib_dwn_dir = fout->mkdir(Form("%s","CalibDwn"));
310  calib_dwn_dir->cd();
311  calib_dwn_ana.SaveTo(calib_dwn_dir->mkdir("Analysis"));
312 
313  TDirectory* calib_up_dir = fout->mkdir(Form("%s","CalibUp"));
314  calib_up_dir->cd();
315  calib_up_ana.SaveTo(calib_up_dir->mkdir("Analysis"));
316 
317  TDirectory* calib_shape_dir = fout->mkdir(Form("%s","CalibShape"));
318  calib_shape_dir->cd();
319  calib_shape_ana.SaveTo(calib_shape_dir->mkdir("Analysis"));
320 
321  TDirectory* ckv_dir = fout->mkdir(Form("%s","Ckv"));
322  ckv_dir->cd();
323  ckv_ana.SaveTo(ckv_dir->mkdir("Analysis"));
324 
325  TDirectory* genie_dir = fout->mkdir(Form("%s","GENIE"));
326  genie_dir->cd();
327  for(int i = 0; i < 100; i++){
328  genie_universes_ana[i]->SaveTo(genie_dir->mkdir(Form("%s_%i","Universe",i)));
329  }
330 
331  TDirectory* ppfx_dir = fout->mkdir(Form("%s","PPFX"));
332  ppfx_dir->cd();
333  for(int i = 0; i < 100; i++){
334  ppfx_universes_ana[i]->SaveTo(ppfx_dir->mkdir(Form("%s_%i","Universe",i)));
335  }
336 
337 
338  fout->Close();
339 
340 
341 
342 }
caf::Proxy< caf::SRFuzzyK >::npng
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
ana::nueccinc::UPLightFiles_full_prod4_CVNProng
const std::string UPLightFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:53
caf::Proxy< caf::SRElastic >::fuzzyk
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
ana::kIsTrueSigST
const NuTruthCut kIsTrueSigST
Definition: NumuCCIncCuts.cxx:205
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::getAllXsecSysts_2018
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
Definition: XSecSystLists.cxx:460
ana::nueccinc::DWNLightFiles_full_prod4_CVNProng
const std::string DWNLightFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:56
kcvn2017
const Var kcvn2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
klid
const Var klid([](const caf::SRProxy *sr){return sr->sel.lid.ann;})
ana::chns
const SelDef chns[knumchns]
Definition: NCPi0SemiIncCuts.h:309
ana::vtxmin
const TVector3 vtxmin(-130,-176, 225)
ana::nueccinc::kcPresel
const Cut kcPresel
Definition: NueCCIncCuts.h:139
submit_syst.fout
fout
Definition: submit_syst.py:244
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
ana::nueccinc::analysis_vars_truth
const mHistAxisSTDef analysis_vars_truth[kvAnalysisVarsST]
Definition: NueCCIncVars.h:489
ana::SpectrumLoaderBase::SetSpillCut
void SetSpillCut(const SpillCut &cut)
Definition: SpectrumLoaderBase.cxx:64
ana::nueccinc::CherenkovFiles_full_prod4_CVNProng
const std::string CherenkovFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:59
ana::GetkPPFXFluxUnivWgt
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
nuebarccinc::pidbins
const ana::Binning pidbins
Definition: NuebarCCIncBins.h:31
caf::Proxy< caf::SRVertexBranch >::elastic
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
kprongCVN
const Var kprongCVN([](const caf::SRProxy *sr){if(sr->vtx.nelastic< 1) return-1000.f;if(sr->vtx.elastic[0].fuzzyk.npng< 1) return-1000.f;int nProngs=sr->vtx.elastic[0].fuzzyk.npng;float maxProngCVNe=-5.;for(int iProng=0;iProng< nProngs;iProng++){if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >maxProngCVNe){maxProngCVNe=sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;}}return maxProngCVNe;})
ana::analysis_vars
const mHistAxisDef analysis_vars[6]
Definition: NCPi0SemiIncVars.h:449
ana::enubins
const Binning enubins
Definition: NumuCCIncBins.h:109
if
if(dump)
Definition: checkMeasured_Target.C:110
caf::Proxy< caf::SRIDBranch >::cvn
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< caf::SRCVNResult >::nueid
caf::Proxy< float > nueid
Definition: SRProxy.h:906
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::SpectrumLoader::Go
virtual void Go() override
Load all the registered spectra.
Definition: SpectrumLoader.cxx:66
ana::GetkPPFXFluxUnivWgtST
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
ana::kXSecCVWgt2018_smallerDISScale_NT
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:47
kcvnProd3
const Var kcvnProd3([](const caf::SRProxy *sr){return(sr->sel.cvnProd3Train.nueid);})
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::SRIDBranch >::lid
caf::Proxy< caf::SRELid > lid
Definition: SRProxy.h:1261
ana::nueccinc::kRecoElecEVsCosStandardAxis
const HistAxis kRecoElecEVsCosStandardAxis("Reconstructed E_{e} vs cos #theta; Elec_{e} (GeV); cos #{theta}", double_diff_bins, kRecoElecEVsCos)
caf::Proxy< caf::SRELid >::ann
caf::Proxy< float > ann
Definition: SRProxy.h:971
ana::nueccinc::ShapeCalibFiles_full_prod4_CVNProng
const std::string ShapeCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:50
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::nueccinc::kcNumChns
const int kcNumChns
Definition: NueCCIncCuts.h:293
ana::kXSecCVWgt2018
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
ana::GenieMultiverseParameters::GetSystShifts
std::vector< SystShifts > GetSystShifts()
Definition: GenieMultiverseSyst.h:248
genie_dir
string genie_dir(std::getenv("GENIE"))
ana::kStandardSpillCuts
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< caf::StandardRecord >::sel
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
ana::nueccinc::kTrueElecEVsCosStandardAxisST
const NuTruthHistAxis kTrueElecEVsCosStandardAxisST("True E_{e} vs cos #theta;cs Elec_{e} (GeV); cos #{theta}", double_diff_bins, kTrueElecEVsCosST)
ana::nueccinc::NominalMC_full_prod4_CVNProng
const std::string NominalMC_full_prod4_CVNProng
Definition: NueCCIncExtra.h:29
caf::Proxy< caf::StandardRecord >::vtx
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
ana::vtxmax
const TVector3 vtxmax(160, 160, 1000)
ana::nueccinc::DWNCalibFiles_full_prod4_CVNProng
const std::string DWNCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:44
ana::_Var< caf::SRProxy >
getCrossSectionAnalysisSpectra
void getCrossSectionAnalysisSpectra()
Definition: getCrossSectionAnalysisSpectra.C:68
kcvn
const Var kcvn([](const caf::SRProxy *sr){return(sr->sel.cvn.nueid);})
ana::nueccinc::UPCalibFiles_full_prod4_CVNProng
const std::string UPCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:47
ana::kPPFXFluxCVWgtST
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
kNumPIDVars
const int kNumPIDVars
Definition: getCrossSectionAnalysisSpectra.C:57
pid_vars
const mHistAxisDef pid_vars[kNumPIDVars]
Definition: getCrossSectionAnalysisSpectra.C:58
Generated on Tue Jan 26 2021 14:12:46 for NOvA by   doxygen 1.8.11