getCrossSectionAnalysis_Spectra_systematics.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
7 #include "CAFAna/XSec/Flux.h"
9 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
16 #include "TVector3.h"
20 
21 using namespace ana;
22 
23 const Var kcvn ([](const caf::SRProxy* sr){
24  return (sr->sel.cvn.nueid);
25  });
26 
27 const Var kcvn2017([] (const caf::SRProxy* sr){
28  return (sr->sel.cvn2017.nueid);
29  });
30 
31 const Var kcvnProd3([] (const caf::SRProxy* sr){
32  return (sr->sel.cvnProd3Train.nueid);
33  });
34 
35 const Var klid([](const caf::SRProxy* sr){
36  return sr->sel.lid.ann;
37  });
38 
39 //sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.photonid
40 const Var kprongCVN([](const caf::SRProxy* sr){
41  if (sr->vtx.nelastic < 1) return -1000.f;
42  if (sr->vtx.elastic[0].fuzzyk.npng < 1) return -1000.f;
43 
44  int nProngs = sr->vtx.elastic[0].fuzzyk.npng;
45  float maxProngCVNe = -5.;
46  for(int iProng = 0; iProng < nProngs; iProng++){
47  if(sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid >
48  maxProngCVNe) {
49  maxProngCVNe =
50  sr->vtx.elastic[0].fuzzyk.png[iProng].cvnpart.electronid;
51  }
52  }
53  return maxProngCVNe;
54  });
55 
56 const int kNumPIDVars = 1;
58  //{"cvne", {"CVNe", Binning::Simple(40,0,1), kcvn}},
59  //{"cvne2017", {"CVNe 2017", pidbins, kcvn2017}},
60  //{"cvneProd3", {"CVNe Prod 3", Binning::Simple(40,0,1), kcvnProd3}},
61  //{"lid", {"LID", Binning::Simple(40,-0.1,1.1), klid}},
62  {"prongCVN", {"Electron Prong Selector", pidbins, kprongCVN}},
63  //{"twoviewCVN", {"Electron Prong Selector 2 View",
64  // Binning::Simple(40,0,1), ktwoviewCVN}},
65  //{"flatfluxCVN", {"Electron Prong Selector Flat Flux",
66  // Binning::Simple(40,0,1), kflatfluxCVN}},
67  //{"genielikeCVN", {"Electron Prong Selector GENIE Like",
68  // Binning::Simple(40,0,1), kgenielikeCVN}},
69 };
70 
71 
72 const Var kMySignalWgt([](const caf::SRProxy* sr)->float
73  {
74  if(sr->mc.nnu == 0) return 0.;
75  assert(sr->mc.nnu == 1);
76  if( (sr->mc.nu[0].iscc)
77  && (abs(sr->mc.nu[0].pdg) == 12)) return 1.06;
78  else return 1.;
79  });
80 
81 const NuTruthVar kMySignalWgtNT([](const caf::SRNeutrinoProxy* sr)->float
82  {
83  if( (sr->iscc)
84  && (abs(sr->pdg) == 12)) return 1.06;
85  else return 1.;
86  });
87 
88 //Get Flux and Efficiency Systematic Histograms
90 {
91  std::vector<std::string> dataset_labels = {"nominal","lightdown",
92  "lightup", "ckv", "calibneg",
93  "calibpos","calibshape"};
94 
95  std::vector<std::string> datasets =
96  {
104  };
105 
106  std::vector<ana::SpectrumLoaderBase*> loaders;
107  for(uint idataset = 0; idataset < datasets.size(); ++idataset){
109  loader = new SpectrumLoader(datasets[idataset]);
111  loaders.push_back(loader);
112  }//end loop over loaders
113 
114 
115  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
116  std::vector<ana::NuTruthVar> ppfx_univ_st = ana::GetkPPFXFluxUnivWgtST();
117 
118  GenieMultiverseParameters verse(100,
120  std::vector<ana::SystShifts> genie_shifts = verse.GetSystShifts();
121 
122  //Spectra
123  std::vector<ana::Spectrum*> spec_efficiency_num_1d;
124  std::vector<ana::Spectrum*> spec_efficiency_num_2d;
125  std::vector<ana::Spectrum*> spec_efficiency_denom_1d;
126  std::vector<ana::Spectrum*> spec_efficiency_denom_2d;
127 
128 
129  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_xsec_num_1d;
130  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_xsec_num_2d;
131  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_xsec_denom_1d;
132  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_xsec_denom_2d;
133 
134  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_ppfx_num_1d;
135  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_ppfx_num_2d;
136  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_ppfx_denom_1d;
137  std::vector<std::vector<ana::Spectrum*>> spec_efficiency_ppfx_denom_2d;
138 
139  std::vector<ana::Spectrum*> spec_flux;
140  std::vector<std::vector<ana::Spectrum*>> spec_flux_ppfx;
141 
142  const TVector3* fid_min = new TVector3(-130,-140,150);
143  const TVector3* fid_max = new TVector3(150,140,800);
144 
145 
146  for(uint iload = 0; iload < loaders.size(); ++iload){
147  spec_efficiency_num_1d.push_back(new Spectrum(*loaders[iload],
149  kcPresel &&
151  kNoShift,
153  kXSecCVWgt2018));
154  spec_efficiency_num_2d.push_back(new Spectrum(*loaders[iload],
156  kcPresel &&
158  kNoShift,
160  kXSecCVWgt2018));
161 
162  spec_efficiency_denom_1d.push_back(new Spectrum(*loaders[iload],
163  analysis_vars_truth[0].axis, // true nu energy
164  kIsTrueSigST,
165  kNoShift,
168  spec_efficiency_denom_2d.push_back(new Spectrum(*loaders[iload],
170  kIsTrueSigST,
171  kNoShift,
173  kXSecCVWgt2018_smallerDISScale_NT));
174  spec_flux.push_back((Spectrum*)
175  DeriveFlux(*loaders[iload],
176  enubins,
177  -12, //12,
178  fid_min,
179  fid_max,
180  kNoShift,
182 
183  if(iload == 0){
184  std::vector<ana::Spectrum*>vspec_xsec_num_1d;
185  std::vector<ana::Spectrum*>vspec_xsec_num_2d;
186  std::vector<ana::Spectrum*>vspec_xsec_denom_1d;
187  std::vector<ana::Spectrum*>vspec_xsec_denom_2d;
188  std::vector<ana::Spectrum*>vspec_ppfx_num_1d;
189  std::vector<ana::Spectrum*>vspec_ppfx_num_2d;
190  std::vector<ana::Spectrum*>vspec_ppfx_denom_1d;
191  std::vector<ana::Spectrum*>vspec_ppfx_denom_2d;
192  std::vector<ana::Spectrum*>vspec_flux_ppfx;
193 
194  for(int iuniv = 0; iuniv < 100; iuniv++){
195  vspec_xsec_num_1d.push_back(new Spectrum
196  (*loaders[iload],
197  HistAxisFromNuTruthHistAxis(analysis_vars_truth[0].axis), // true nu energy
198  kcPresel &&
200  genie_shifts[iuniv],
202  kXSecCVWgt2018));
203  vspec_xsec_num_2d.push_back(new Spectrum
204  (*loaders[iload],
206  kcPresel &&
208  genie_shifts[iuniv],
209  kPPFXFluxCVWgt*kXSecCVWgt2018));
210  vspec_ppfx_num_1d.push_back(new Spectrum
211  (*loaders[iload],
212  HistAxisFromNuTruthHistAxis(analysis_vars_truth[0].axis), // true nu energy
213  kcPresel &&
215  kNoShift,
216  ppfx_univ[iuniv] *
217  kXSecCVWgt2018));
218  vspec_ppfx_num_2d.push_back(new Spectrum
219  (*loaders[iload],
221  kcPresel &&
223  kNoShift,
224  ppfx_univ[iuniv] *
225  kXSecCVWgt2018));
226 
227  vspec_xsec_denom_1d.push_back(new Spectrum
228  (*loaders[iload],
229  analysis_vars_truth[0].axis, // true nu energy
230  kIsTrueSigST,
231  genie_shifts[iuniv],
233  kXSecCVWgt2018_smallerDISScale_NT));
234  vspec_xsec_denom_2d.push_back(new Spectrum
235  (*loaders[iload],
237  kIsTrueSigST,
238  genie_shifts[iuniv],
240  kXSecCVWgt2018_smallerDISScale_NT));
241 
242  vspec_ppfx_denom_1d.push_back(new Spectrum
243  (*loaders[iload],
244  analysis_vars_truth[0].axis, // true nu energy
245  kIsTrueSigST,
246  kNoShift,
247  ppfx_univ_st[iuniv]*
248  kXSecCVWgt2018_smallerDISScale_NT));
249  vspec_ppfx_denom_2d.push_back(new Spectrum
250  (*loaders[iload],
252  kIsTrueSigST,
253  kNoShift,
254  ppfx_univ_st[iuniv] *
255  kXSecCVWgt2018_smallerDISScale_NT));
256 
257  vspec_flux_ppfx.push_back((Spectrum*)
258  DeriveFlux(*loaders[iload],
259  enubins,
260  -12, // 12,
261  fid_min,
262  fid_max,
263  kNoShift,
264  ppfx_univ_st[iuniv]));
265  }
266 
267  spec_efficiency_xsec_num_1d.push_back(vspec_xsec_num_1d);
268  spec_efficiency_xsec_num_2d.push_back(vspec_xsec_num_2d);
269  spec_efficiency_xsec_denom_1d.push_back(vspec_xsec_denom_1d);
270  spec_efficiency_xsec_denom_2d.push_back(vspec_xsec_denom_2d);
271  spec_efficiency_ppfx_num_1d.push_back(vspec_ppfx_num_1d);
272  spec_efficiency_ppfx_num_2d.push_back(vspec_ppfx_num_2d);
273  spec_efficiency_ppfx_denom_1d.push_back(vspec_ppfx_denom_1d);
274  spec_efficiency_ppfx_denom_2d.push_back(vspec_ppfx_denom_2d);
275  spec_flux_ppfx.push_back(vspec_flux_ppfx);
276  } //if statement
277  }//loop over loaders
278 
279  for(uint iload = 0; iload < loaders.size(); ++iload){
280  loaders[iload]->Go();
281  }
282 
283  TFile* outf = new TFile("fCrossSection_EfficiencyFlux_Systs.root","recreate");
284 
285  for(uint iload = 0; iload < loaders.size(); ++iload){
286 
287  char name[50];
288  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
289  "efficiency", "num_1d");
290  spec_efficiency_num_1d[iload]->SaveTo(outf, name);
291  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
292  "efficiency", "num_2d");
293  spec_efficiency_num_2d[iload]->SaveTo(outf, name);
294  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
295  "efficiency", "denom_1d");
296  spec_efficiency_denom_1d[iload]->SaveTo(outf, name);
297  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[iload].c_str(),
298  "efficiency", "denom_2d");
299  spec_efficiency_denom_2d[iload]->SaveTo(outf, name);
300  sprintf(name, "%s_%s_%s", "mc", dataset_labels[iload].c_str(),
301  "flux");
302  spec_flux[iload]->SaveTo(outf, name);
303 
304  if(iload == 0){
305  for(int iuniv = 0; iuniv < 100; iuniv++){
306  sprintf(name, "%s_%s_%i_%s_%s", "mc", "genie", iuniv,
307  "efficiency", "num_1d");
308  spec_efficiency_xsec_num_1d[iload][iuniv]->SaveTo(outf, name);
309  sprintf(name, "%s_%s_%i_%s_%s", "mc", "genie", iuniv,
310  "efficiency", "num_2d");
311  spec_efficiency_xsec_num_2d[iload][iuniv]->SaveTo(outf, name);
312  sprintf(name, "%s_%s_%i_%s_%s", "mc", "genie", iuniv,
313  "efficiency", "denom_1d");
314  spec_efficiency_xsec_denom_1d[iload][iuniv]->SaveTo(outf, name);
315  sprintf(name, "%s_%s_%i_%s_%s", "mc", "genie", iuniv,
316  "efficiency", "denom_2d");
317  spec_efficiency_xsec_denom_2d[iload][iuniv]->SaveTo(outf, name);
318  sprintf(name, "%s_%s_%i_%s_%s", "mc", "ppfx", iuniv,
319  "efficiency", "num_1d");
320  spec_efficiency_ppfx_num_1d[iload][iuniv]->SaveTo(outf, name);
321  sprintf(name, "%s_%s_%i_%s_%s", "mc", "ppfx", iuniv,
322  "efficiency", "num_2d");
323  spec_efficiency_ppfx_num_2d[iload][iuniv]->SaveTo(outf, name);
324  sprintf(name, "%s_%s_%i_%s_%s", "mc", "ppfx", iuniv,
325  "efficiency", "denom_1d");
326  spec_efficiency_ppfx_denom_1d[iload][iuniv]->SaveTo(outf, name);
327  sprintf(name, "%s_%s_%i_%s_%s", "mc", "ppfx", iuniv,
328  "efficiency", "denom_2d");
329  spec_efficiency_ppfx_denom_2d[iload][iuniv]->SaveTo(outf, name);
330  sprintf(name, "%s_%s_%i_%s", "mc", "ppfx", iuniv,
331  "flux");
332  spec_flux_ppfx[iload][iuniv]->SaveTo(outf, name);
333  }//loop over universes
334  } //if statement
335  }//loop over loaders
336  outf->Close();
337 
338 }
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
const Var kMySignalWgt([](const caf::SRProxy *sr) ->float{if(sr->mc.nnu==0) return 0.;assert(sr->mc.nnu==1);if((sr->mc.nu[0].iscc) &&(abs(sr->mc.nu[0].pdg)==12)) return 1.06;else return 1.;})
const XML_Char * name
Definition: expat.h:151
const Var kcvn2017([](const caf::SRProxy *sr){return(sr->sel.cvn2017.nueid);})
const std::string NominalMC_entire_RHC_prod4
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 Var kcvnProd3([](const caf::SRProxy *sr){return(sr->sel.cvnProd3Train.nueid);})
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kcPresel
Definition: NueCCIncCuts.h:139
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
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;})
std::vector< Dataset > datasets
Definition: Periods.h:3
const mHistAxisSTDef analysis_vars_truth[kvAnalysisVarsST]
Definition: NueCCIncVars.h:489
void abs(TH1 *hist)
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
const Var kcvn([](const caf::SRProxy *sr){return(sr->sel.cvn.nueid);})
const std::string DWNLightFiles_full_RHC_prod4
void getCrossSectionAnalysis_Spectra_systematics()
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const ana::Binning pidbins
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
const std::string ShapeCalibFiles_full_RHC_prod4
const NuTruthVar kMySignalWgtNT([](const caf::SRNeutrinoProxy *sr) ->float{if((sr->iscc) &&(abs(sr->pdg)==12)) return 1.06;else return 1.;})
const Binning enubins
if(dump)
const mHistAxisDef pid_vars[kNumPIDVars]
const std::string UPCalibFiles_full_RHC_prod4
caf::Proxy< caf::SRCVNResult > cvn
Definition: SRProxy.h:1253
caf::Proxy< float > nueid
Definition: SRProxy.h:906
caf::Proxy< bool > iscc
Definition: SRProxy.h:538
caf::StandardRecord * sr
HistAxis HistAxisFromNuTruthHistAxis(NuTruthHistAxis ntha, double _default)
Definition: HistAxis.cxx:9
TFile * outf
Definition: testXsec.C:51
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const std::string UPLightFiles_full_RHC_prod4
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:47
loader
Definition: demo0.py:10
const HistAxis kTrueElecEVsCosStandardAxis
Definition: NueCCIncVars.h:533
std::vector< float > Spectrum
Definition: Constants.h:759
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::SRELid > lid
Definition: SRProxy.h:1261
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< float > ann
Definition: SRProxy.h:971
Spectrum * DeriveFlux(SpectrumLoaderBase &loader, const Binning &bins, int pdg, const TVector3 *min, const TVector3 *max, const SystShifts &shift, const NuTruthVar weight)
Definition: Flux.cxx:58
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
caf::Proxy< short int > pdg
Definition: SRProxy.h:552
const std::string DWNCalibFiles_full_RHC_prod4
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 CherenkovFiles_full_RHC_prod4
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const Var klid([](const caf::SRProxy *sr){return sr->sel.lid.ann;})
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
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:7
unsigned int uint