getCrossSectionAnalysisSpectra.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
8 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
15 #include "CAFAna/Vars/XsecTunes.h"
16 
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;
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 
69 {
77 
79  lNDLightDwn.SetSpillCut(kStandardSpillCuts);
80  lNDLightUp.SetSpillCut(kStandardSpillCuts);
81  lNDCalibDwn.SetSpillCut(kStandardSpillCuts);
82  lNDCalibUp.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,
95  analysis_vars[1].axis,
96  analysis_vars[0].axis,
99  analysis_vars_truth[0].axis,
100  pid_vars[0].axis,
101  enubins,
102  kcPresel,
103  kIsTrueSigST,
104  vtxmin,
105  vtxmax,
106  kNoShift,
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,
130  analysis_vars_truth[0].axis,
131  pid_vars[0].axis,
132  enubins,
133  kcPresel,
134  kIsTrueSigST,
135  vtxmin,
136  vtxmax,
137  kNoShift,
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,
148  analysis_vars_truth[0].axis,
149  pid_vars[0].axis,
150  enubins,
151  kcPresel,
152  kIsTrueSigST,
153  vtxmin,
154  vtxmax,
155  kNoShift,
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,
166  analysis_vars_truth[0].axis,
167  pid_vars[0].axis,
168  enubins,
169  kcPresel,
170  kIsTrueSigST,
171  vtxmin,
172  vtxmax,
173  kNoShift,
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,
184  analysis_vars_truth[0].axis,
185  pid_vars[0].axis,
186  enubins,
187  kcPresel,
188  kIsTrueSigST,
189  vtxmin,
190  vtxmax,
191  kNoShift,
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,
202  analysis_vars_truth[0].axis,
203  pid_vars[0].axis,
204  enubins,
205  kcPresel,
206  kIsTrueSigST,
207  vtxmin,
208  vtxmax,
209  kNoShift,
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,
221  analysis_vars_truth[0].axis,
222  pid_vars[0].axis,
223  enubins,
224  kcPresel,
225  kIsTrueSigST,
226  vtxmin,
227  vtxmax,
228  kNoShift,
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();
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,
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],
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,
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< size_t > npng
Definition: SRProxy.h:2038
const std::string UPLightFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:53
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
const NuTruthCut kIsTrueSigST
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
const std::string DWNLightFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:56
const Var kcvn2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var klid([](const caf::SRProxy *sr){return sr->sel.lid.ann;})
const SelDef chns[knumchns]
const TVector3 vtxmin(-130,-176, 225)
const Cut kcPresel
Definition: NueCCIncCuts.h:139
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const mHistAxisSTDef analysis_vars_truth[kvAnalysisVarsST]
Definition: NueCCIncVars.h:489
void SetSpillCut(const SpillCut &cut)
const std::string CherenkovFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:59
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const ana::Binning pidbins
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
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;})
const mHistAxisDef analysis_vars[6]
const Binning enubins
if(dump)
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< float > nueid
Definition: SRProxy.h:906
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:47
const Var kcvnProd3([](const caf::SRProxy *sr){return(sr->sel.cvnProd3Train.nueid);})
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::SRELid > lid
Definition: SRProxy.h:1261
const HistAxis kRecoElecEVsCosStandardAxis("Reconstructed E_{e} vs cos #theta; Elec_{e} (GeV); cos #{theta}", double_diff_bins, kRecoElecEVsCos)
caf::Proxy< float > ann
Definition: SRProxy.h:971
const std::string ShapeCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:50
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const int kcNumChns
Definition: NueCCIncCuts.h:293
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< SystShifts > GetSystShifts()
string genie_dir(std::getenv("GENIE"))
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
const NuTruthHistAxis kTrueElecEVsCosStandardAxisST("True E_{e} vs cos #theta;cs Elec_{e} (GeV); cos #{theta}", double_diff_bins, kTrueElecEVsCosST)
const std::string NominalMC_full_prod4_CVNProng
Definition: NueCCIncExtra.h:29
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const TVector3 vtxmax(160, 160, 1000)
const std::string DWNCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:44
void getCrossSectionAnalysisSpectra()
const Var kcvn([](const caf::SRProxy *sr){return(sr->sel.cvn.nueid);})
const std::string UPCalibFiles_full_prod4_CVNProng
Definition: NueCCIncExtra.h:47
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
const int kNumPIDVars
const mHistAxisDef pid_vars[kNumPIDVars]