getSpectra_ForFitting.C
Go to the documentation of this file.
1 //Script to save spectra for both MC and fakedata that are necessary for background estimation study.
2 //Currently set up for the pi^{0} kinetic energy
3 
4 #include <string>
5 
6 #include "CAFAna/Core/Var.h"
7 #include "CAFAna/Core/Cut.h"
8 #include "CAFAna/Vars/Vars.h"
9 
11 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Core/ISyst.h"
13 #include "CAFAna/Core/SystShifts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
24 #include "CAFAna/Systs/XSecSysts.h"
25 
30 
31 // Names for the GENIE systs
32 #include "nugen/NuReweight/ReweightLabels.h"
33 #include "ReWeight/GSyst.h"
34 
35 
36 using namespace ana;
37 
39 {
40 
41  //Weights for Background Estimation Tests
42  const NuTruthVar weightSigUpST = kPPFXFluxCVWgtST * kXSecCVWgt2017ST * kSigUpST;
43  const Var kMySignalWgt = VarFromNuTruthVar(weightSigUpST, 1);
44 
45  const NuTruthVar weightHighSigUpST = kPPFXFluxCVWgtST * kXSecCVWgt2017ST * kSigHighWtST;
46  const Var kMySignalHighWgt = VarFromNuTruthVar(weightHighSigUpST, 1);
47 
48  const NuTruthVar weightBkgUpST = kPPFXFluxCVWgtST * kXSecCVWgt2017ST * kBkgUpST;
49  const Var kMyBkgWgt = VarFromNuTruthVar(weightBkgUpST, 1);
50 
51 
52  const std::string fFakeData = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with stride 4 with offset 0 with limit 1";// with limit 1 with limit 1";
53 
54  const std::string fGENIE = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with stride 4 with offset 1 with limit 1";// with limit 1 with limit 1";
55 
56  const std::string fCalibDwn = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-neg-offset_v1 with limit 1";// with limit 1 with limit 1"; 9710 files
57 
58  const std::string fCalibUp = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-pos-offset_v1 with limit 1";// with limit 1 with limit 1"; 8967 files
59 
60  const std::string fLightDwn = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1 with limit 1";// with limit 1 with limit 1"; 7992 files
61 
62  const std::string fLightUp = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1 with limit 1";// with limit 1 with limit 1"; 7842 files
63 
64  const std::string fCalibShape = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.j_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1 with limit 1";// with limit 1 with limit 1"; 9073 files
65 
66  const std::string fCherenkov = "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1 with limit 1";// with limit 1"; 7427 files
67 
68 
69  SpectrumLoader lNDDA(fFakeData);
70  SpectrumLoader lNDMC(fGENIE);
71 
72  SpectrumLoader lNDUPCal(fCalibUp);
73  SpectrumLoader lNDDWNCal(fCalibDwn);
74 
75  SpectrumLoader lNDUPLight(fLightUp);
76  SpectrumLoader lNDDWNLight(fLightDwn);
77 
78  SpectrumLoader lNDCalibShape(fCalibShape);
79  SpectrumLoader lNDCherenkov(fCherenkov);
80 
81 
90 
91  //PPFX Multiverse
92  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
93  int numUniverses = (int)ppfx_univ.size();
94  Spectrum* ppfx_spects[numUniverses];
95  for(int i = 0; i < numUniverses; i++){
96  ppfx_spects[i] = new Spectrum(lNDMC, analysis_vars[4].axis,
97  analysis_vars_anabins[0].axis,
98  kIsPreSel,
99  kNoShift,
101  ppfx_univ[i]);
102  }
103 
104 
105  //Nominal MC and Components
106  Spectrum* spect_cv = new Spectrum(lNDMC, analysis_vars[4].axis,
107  analysis_vars_anabins[0].axis,
108  kIsPreSel,
109  kNoShift,
111 
112 
113  Spectrum* spect_templates[knumchns];
114  for(int i = 0; i < knumchns; i++){
115  spect_templates[i] = new Spectrum(lNDMC, analysis_vars[4].axis,
116  analysis_vars_anabins[0].axis,
117  kIsPreSel && chns[i].cut,
118  kNoShift,
119  kPPFXFluxCVWgt * kXSecCVWgt2017);
120  }
121 
122 
123 
124  //Genie Multiverse
125  std::vector<const ISyst*> systs = getAllXsecSysts_2017();
126 
127  GenieMultiverseParameters genieshifts(100,
128  systs,
129  {});//ana::getAllXsecNuTruthSysts_2017());
130 
131  std::vector<SystShifts> genie_univ = genieshifts.GetSystShifts();
132 
133  Spectrum* genie_spects[100];
134  for(int i = 0; i < 100; i++){
135  genie_spects[i] = new Spectrum(lNDMC, analysis_vars[4].axis,
136  analysis_vars_anabins[0].axis,
137  kIsPreSel,
138  genie_univ[i],
139  kPPFXFluxCVWgt * kXSecCVWgt2017);
140  }
141 
142 
143  // for various other systematic samples
144 
145  //Systematic Samples
146  Spectrum *upcal_spect =
147  new Spectrum(lNDUPCal, analysis_vars[4].axis,
148  analysis_vars_anabins[0].axis,
149  kIsPreSel,
150  kNoShift,
151  kXSecCVWgt2017 * kPPFXFluxCVWgt);
152 
153 
154  Spectrum *dwncal_spect =
155  new Spectrum(lNDDWNCal,analysis_vars[4].axis,
156  analysis_vars_anabins[0].axis,
157  kIsPreSel,
158  kNoShift,
159  kXSecCVWgt2017 * kPPFXFluxCVWgt);
160 
161 
162 
163  Spectrum *lightup_spect =
164  new Spectrum(lNDUPLight, analysis_vars[4].axis,
165  analysis_vars_anabins[0].axis,
166  kIsPreSel,
167  kNoShift,
168  kXSecCVWgt2017 * kPPFXFluxCVWgt);
169 
170  Spectrum *lightdwn_spect =
171  new Spectrum(lNDDWNLight, analysis_vars[4].axis,
172  analysis_vars_anabins[0].axis,
173  kIsPreSel,
174  kNoShift,
175  kXSecCVWgt2017 * kPPFXFluxCVWgt);
176 
177 
178  Spectrum *calibshape_spect =
179  new Spectrum(lNDCalibShape, analysis_vars[4].axis,
180  analysis_vars_anabins[0].axis,
181  kIsPreSel,
182  kNoShift,
183  kXSecCVWgt2017 * kPPFXFluxCVWgt);
184 
185  Spectrum *cherenkov_spect =
186  new Spectrum(lNDCherenkov, analysis_vars[4].axis,
187  analysis_vars_anabins[0].axis,
188  kIsPreSel,
189  kNoShift,
190  kXSecCVWgt2017 * kPPFXFluxCVWgt);
191 
192 
193 
194 
195  //Data and Components (if fake data)
196 
197  Spectrum* data_cv = new Spectrum(lNDDA, analysis_vars[4].axis,
198  analysis_vars_anabins[0].axis,
199  kIsPreSel,
200  kNoShift,
201  kXSecCVWgt2017 * kPPFXFluxCVWgt);
202 
203  Spectrum* data_templates[knumchns];
204  for(int i = 0; i < knumchns; i++){
205  data_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
206  analysis_vars_anabins[0].axis,
207  kIsPreSel && chns[i].cut,
208  kNoShift,
209  kPPFXFluxCVWgt * kXSecCVWgt2017);
210  }
211 
212  Spectrum* data_sig = new Spectrum(lNDDA, analysis_vars[4].axis,
213  analysis_vars_anabins[0].axis,
214  kIsPreSel,
215  kNoShift,
216  kMySignalWgt);
217 
218  Spectrum* data_sig_templates[knumchns];
219  for(int i = 0; i < knumchns; i++){
220  data_sig_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
221  analysis_vars_anabins[0].axis,
222  kIsPreSel && chns[i].cut,
223  kNoShift,
224  kMySignalWgt);
225  }
226 
227 
228  Spectrum* data_sighigh = new Spectrum(lNDDA, analysis_vars[4].axis,
229  analysis_vars_anabins[0].axis,
230  kIsPreSel,
231  kNoShift,
232  kMySignalHighWgt);
233 
234  Spectrum* data_sighigh_templates[knumchns];
235  for(int i = 0; i < knumchns; i++){
236  data_sighigh_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
237  analysis_vars_anabins[0].axis,
238  kIsPreSel && chns[i].cut,
239  kNoShift,
240  kMySignalHighWgt);
241  }
242 
243  Spectrum* data_bkg = new Spectrum(lNDDA, analysis_vars[4].axis,
244  analysis_vars_anabins[0].axis,
245  kIsPreSel,
246  kNoShift,
247  kMyBkgWgt);
248 
249  Spectrum* data_bkg_templates[knumchns];
250  for(int i = 0; i < knumchns; i++){
251  data_bkg_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
252  analysis_vars_anabins[0].axis,
253  kIsPreSel && chns[i].cut,
254  kNoShift,
255  kMyBkgWgt);
256  }
257 
258 
259 
260  Spectrum* data_mancres_sighigh = new Spectrum(lNDDA, analysis_vars[4].axis,
261  analysis_vars_anabins[0].axis,
262  kIsPreSel,
263  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaNCRES), +1),
264  kMySignalHighWgt);
265 
266  Spectrum* data_mancres_sighigh_templates[knumchns];
267  for(int i = 0; i < knumchns; i++){
268  data_mancres_sighigh_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
269  analysis_vars_anabins[0].axis,
270  kIsPreSel && chns[i].cut,
271  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaNCRES), +1),
272  kMySignalHighWgt);
273  }
274 
275  Spectrum* data_mancres = new Spectrum(lNDDA, analysis_vars[4].axis,
276  analysis_vars_anabins[0].axis,
277  kIsPreSel,
278  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaNCRES), +1),
279  kXSecCVWgt2017 * kPPFXFluxCVWgt);
280 
281  Spectrum* data_mancres_templates[knumchns];
282  for(int i = 0; i < knumchns; i++){
283  data_mancres_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
284  analysis_vars_anabins[0].axis,
285  kIsPreSel && chns[i].cut,
286  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaNCRES), +1),
287  kXSecCVWgt2017 * kPPFXFluxCVWgt);
288  }
289 
290  Spectrum* data_maccres = new Spectrum(lNDDA, analysis_vars[4].axis,
291  analysis_vars_anabins[0].axis,
292  kIsPreSel,
293  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaCCRES), +1),
294  kXSecCVWgt2017 * kPPFXFluxCVWgt);
295 
296  Spectrum* data_maccres_templates[knumchns];
297  for(int i = 0; i < knumchns; i++){
298  data_maccres_templates[i] = new Spectrum(lNDDA, analysis_vars[4].axis,
299  analysis_vars_anabins[0].axis,
300  kIsPreSel && chns[i].cut,
301  SystShifts(GetGenieRwgtNuTruthSyst(rwgt::fReweightMaCCRES), +1),
302  kXSecCVWgt2017 * kPPFXFluxCVWgt);
303  }
304 
305 
306 
307  lNDDA.Go();
308  lNDMC.Go();
309  lNDUPCal.Go();
310  lNDDWNCal.Go();
311  lNDUPLight.Go();
312  lNDDWNLight.Go();
313  lNDCalibShape.Go();
314  lNDCherenkov.Go();
315 
316 
317  std::string outname_flux = "PPFX_universes_bdte_energy_new_bdt_bins.root";
318  TFile *flux_file = new TFile(outname_flux.c_str(), "recreate");
319 
320  std::string outname_genie = "GENIE_universes_bdte_energy_new_bdt_bins.root";
321  TFile *genie_file = new TFile(outname_genie.c_str(), "recreate");
322 
323  std::string outname_systs = "Systs_samples_bdte_energy_new_bdt_bins.root";
324  TFile *systs_file = new TFile(outname_systs.c_str(), "recreate");
325 
326  std::string outname_nominal = "Nominal_MC_bdte_energy_new_bdt_bins.root";
327  TFile *nominal_file = new TFile(outname_nominal.c_str(),"recreate");
328 
329  std::string outname_data = "Data_bdte_energy_new_bdt_bins.root";
330  TFile *data_file = new TFile(outname_data.c_str(),"recreate");
331 
332 
333  for(int i = 0; i < numUniverses; i++){
334  char name_ppfx[50];
335  sprintf(name_ppfx, "%s_%i", "flux_universe", i);
336  ppfx_spects[i]->SaveTo(flux_file->mkdir(name_ppfx));
337 
338  }
339 
340  for(int i = 0; i < 100; i++){
341  char name_genie[50];
342  sprintf(name_genie, "%s_%il", "genie_universe", i);
343  genie_spects[i]->SaveTo(genie_file->mkdir(name_genie));
344  }
345 
346  for(int i = 0; i < knumchns; i++){
347  char name[50];
348  sprintf(name, "%s_%s", "template", chns[i].name.c_str());
349  spect_templates[i]->SaveTo(nominal_file->mkdir(name));
350  }
351  spect_cv->SaveTo(nominal_file->mkdir("nominal_mc"));
352  upcal_spect->SaveTo(systs_file->mkdir("upcal_systs"));
353  dwncal_spect->SaveTo(systs_file->mkdir("dwncal_systs"));
354  lightup_spect->SaveTo(systs_file->mkdir("lightup_systs"));
355  lightdwn_spect->SaveTo(systs_file->mkdir("lightdwn_systs"));
356  calibshape_spect->SaveTo(systs_file->mkdir("calibshape_systs"));
357  cherenkov_spect->SaveTo(systs_file->mkdir("cherenkov_systs"));
358 
359 
360  for(int i = 0; i < knumchns; i++){
361  char name[50];
362  sprintf(name, "%s_%s", "data", chns[i].name.c_str());
363  data_templates[i]->SaveTo(data_file->mkdir(name));
364  }
365  data_cv->SaveTo(data_file->mkdir("nominal_data"));
366 
367 
368  for(int i = 0; i < knumchns; i++){
369  char name[50];
370  sprintf(name, "%s_%s", "data_sig", chns[i].name.c_str());
371  data_sig_templates[i]->SaveTo(data_file->mkdir(name));
372  }
373  data_sig->SaveTo(data_file->mkdir("nominal_data_sig"));
374 
375 
376  for(int i = 0; i < knumchns; i++){
377  char name[50];
378  sprintf(name, "%s_%s", "data_sighigh", chns[i].name.c_str());
379  data_sighigh_templates[i]->SaveTo(data_file->mkdir(name));
380  }
381  data_sighigh->SaveTo(data_file->mkdir("nominal_data_sighigh"));
382 
383  for(int i = 0; i < knumchns; i++){
384  char name[50];
385  sprintf(name, "%s_%s", "data_bkg", chns[i].name.c_str());
386  data_bkg_templates[i]->SaveTo(data_file->mkdir(name));
387  }
388  data_bkg->SaveTo(data_file->mkdir("nominal_data_bkg"));
389 
390 
391  for(int i = 0; i < knumchns; i++){
392  char name[50];
393  sprintf(name, "%s_%s", "data_mancres_sighigh", chns[i].name.c_str());
394  data_mancres_sighigh_templates[i]->SaveTo(data_file->mkdir(name));
395  }
396  data_mancres_sighigh->SaveTo(data_file->mkdir("nominal_data_mancres_sighigh"));
397 
398  for(int i = 0; i < knumchns; i++){
399  char name[50];
400  sprintf(name, "%s_%s", "data_mancres", chns[i].name.c_str());
401  data_mancres_templates[i]->SaveTo(data_file->mkdir(name));
402  }
403  data_mancres->SaveTo(data_file->mkdir("nominal_data_mancres"));
404 
405  for(int i = 0; i < knumchns; i++){
406  char name[50];
407  sprintf(name, "%s_%s", "data_maccres", chns[i].name.c_str());
408  data_maccres_templates[i]->SaveTo(data_file->mkdir(name));
409  }
410  data_maccres->SaveTo(data_file->mkdir("nominal_data_maccres"));
411 
412  data_file->Close();
413  nominal_file->Close();
414  flux_file->Close();
415  genie_file->Close();
416  systs_file->Close();
417 
418 }
419 
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const SelDef chns[knumchns]
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
void SetSpillCut(const SpillCut &cut)
const int knumchns
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const mHistAxisDef analysis_vars[6]
void getSpectra_ForFitting()
const NuTruthVar kSigUpST([](const caf::SRNeutrinoProxy *nu){float MassOfPi0=0.135;float kinetic=0.0f;float en=0.0f;int nbOfPi=0;int nbofprim=nu->prim.size();for(int i=0;i< nbofprim;i++){if(nu->prim[i].pdg==111){nbOfPi++;en=nu->prim[i].p.E;kinetic=en-MassOfPi0;}}if(((nu->pdg==14 &&nu->pdgorig==14)||(nu->pdg==12 &&nu->pdgorig==12))&&(nu->iscc==0)&&(kinetic >0.1)) return 1.10;else return 1.;})
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:13
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
const mHistAxisDef analysis_vars_anabins[kEnAngVars]
std::vector< float > Spectrum
Definition: Constants.h:527
const Cut kIsPreSel
Definition: ncpi0Cuts.h:1099
const NuTruthVar kSigHighWtST([](const caf::SRNeutrinoProxy *nu){float MassOfPi0=0.135;float kinetic=0.0f;float en=0.0f;int nbOfPi=0;int nbofprim=nu->prim.size();for(int i=0;i< nbofprim;i++){if(nu->prim[i].pdg==111){nbOfPi++;en=nu->prim[i].p.E;kinetic=en-MassOfPi0;}}double wt=1.1+(kinetic *0.11);if(((nu->pdg==14 &&nu->pdgorig==14)||(nu->pdg==12 &&nu->pdgorig==12))&&(nu->iscc==0)&&(kinetic >0.1)) return wt;else return 1.;})
const SystShifts kNoShift
Definition: SystShifts.h:115
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Var VarFromNuTruthVar(const NuTruthVar &stv, double _default)
Definition: Var.cxx:339
const SpillCut kTightBeamQualityCuts([](const caf::SRSpillProxy *s){if(s->ismc) return true; if(s->trigger==2) return true;if(s->spilltimesec==0 &&s->deltaspilltimensec==0 &&s->widthx==0) return bool(s->isgoodspill);if(std::abs(s->deltaspilltimensec) > 0.5e9) return false;if(s->spillpot< 2e12) return false;if(s->hornI< -202|| s->hornI >-198) return false;if(s->posx< -2.00|| s->posx >+2.00) return false;if(s->posy< -2.00|| s->posy >+2.00) return false;int run=s->run; if(s->det==1){if(run< 13471){if(s->widthx< 0.57||s->widthx > 1.58) return false;}else{if(s->widthx< 0.57|| s->widthx > 1.88) return false;}}else{if(run< 34574){if(s->widthx< 0.57|| s->widthx > 1.58) return false;}else{if(s->widthx< 0.57|| s->widthx > 1.88) return false;}}if(s->widthy< 0.57|| s->widthy > 1.58) return false;return true;})
Definition: SpillCuts.h:10
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::vector< SystShifts > GetSystShifts()
const NuTruthVar kXSecCVWgt2017ST
2017 analysis master weight
Definition: XsecTunes.h:36
size_t numUniverses
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 Var kXSecCVWgt2017
Definition: XsecTunes.h:37
const NuTruthVar kBkgUpST([](const caf::SRNeutrinoProxy *nu){float MassOfPi0=0.135;float kinetic=0.0f;float en=0.0f;int nbOfPi=0;int nbofprim=nu->prim.size();for(int i=0;i< nbofprim;i++){if(nu->prim[i].pdg==111){nbOfPi++;en=nu->prim[i].p.E;kinetic=en-MassOfPi0;}}if(!(((nu->pdg==14 &&nu->pdgorig==14)||(nu->pdg==12 &&nu->pdgorig==12))&&(nu->iscc==0)&&(kinetic >0.1))) return 1.1;else return 1.;})
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
std::vector< const ISyst * > getAllXsecSysts_2017()
Get master XSec syst list for 2017 analyses.
Template for Var and SpillVar.
Definition: Var.h:16