getResolutionSpectra.C
Go to the documentation of this file.
1 // I found that the energy estimator used by the nuebar oscillation analysis
2 // was not sufficient for this analysis due to training on a sample of
3 // neutrinos with energy less than 4.5 GeV.
4 //
5 // This script marks the first piece of training a new estimator.
6 // Following docdb-26696, we need to normalize the true energy distribution
7 // to avoid biasing the energy estimation toward our peak beam energy.
8 //
9 // Grab flat-weighted true energy, recoEME, and recoHADE energy distributions
10 // dddoyle@colostate.edu
11 
13 #include "CAFAna/Core/Spectrum.h"
14 #include "CAFAna/Core/Var.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
20 #include "CAFAna/Vars/XsecTunes.h"
22 
27 
28 using namespace ana;
29 
30 // selection used in the nuebar CC inclusive analysis
33 Cut cuts = kcPresel && kIsTrueSig; // ktIsNuebarCCST && kTrueFiducialST
34 std::string training_data = NominalMC_energy_compare_RHC_prod4;
35 
38 
40 {
41  std::vector<std::pair <const Var, std::string> > resolution_vars =
42  {{kFracEnergyRes_poly2, "poly2" },
43  {kFracEnergyRes_poly2X, "poly2x" },
44  {kFracEnergyRes_poly3, "poly3" },
45  {kFracEnergyRes_poly3X, "poly3x" },
46  {kFracEnergyRes_poly2_shifted, "poly2_shifted" },
47  {kFracEnergyRes_poly2X_shifted, "poly2x_shifted" },
48  {kFracEnergyRes_poly3_shifted, "poly3_shifted" },
49  {kFracEnergyRes_poly3X_shifted, "poly3x_shifted" },
50  {kFracEnergyRes_nue2018, "nue2018"},
51  {kFracEnergyRes_nue2018_shifted, "nue2018_shifted"}};
52 
53  std::vector<Spectrum*> frac_resolution;
54  std::vector<Spectrum*> frac_resolution_nocut_shape;
55  std::vector<Spectrum*> frac_resolution_cc_shape;
56  std::vector<Spectrum*> frac_resolution_qe_shape;
57  std::vector<Spectrum*> frac_resolution_res_shape;
58  std::vector<Spectrum*> frac_resolution_dis_shape;
59  std::vector<Spectrum*> frac_resolution_coh_shape;
60  std::vector<Spectrum*> frac_resolution_mec_shape;
61  std::vector<Spectrum*> frac_resolution_other_shape;
62  std::vector<Spectrum*> frac_resolution_nocut_flat;
63  std::vector<Spectrum*> frac_resolution_cc_flat;
64  std::vector<Spectrum*> frac_resolution_qe_flat;
65  std::vector<Spectrum*> frac_resolution_res_flat;
66  std::vector<Spectrum*> frac_resolution_dis_flat;
67  std::vector<Spectrum*> frac_resolution_coh_flat;
68  std::vector<Spectrum*> frac_resolution_mec_flat;
69  std::vector<Spectrum*> frac_resolution_other_flat;
70 
71  std::vector<Spectrum*> resolution_scan_trueE_shape;
72  std::vector<Spectrum*> resolution_scan_EME_shape;
73  std::vector<Spectrum*> resolution_scan_HADE_shape;
74  std::vector<Spectrum*> resolution_scan_trueE_flat;
75  std::vector<Spectrum*> resolution_scan_EME_flat;
76  std::vector<Spectrum*> resolution_scan_HADE_flat;
77  std::vector<Spectrum*> resolution_scan_recoE_shape;
78  std::vector<Spectrum*> resolution_scan_recoE_flat;
79 
81  for(int i = 0; i < (int) resolution_vars.size(); i++) {
82  frac_resolution.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
83  resolution_binning,
84  *loader,
85  resolution_vars[i].first,
86  cuts,
87  kNoShift,
88  weight_shape));
89 
90  frac_resolution_nocut_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
91  resolution_binning,
92  *loader,
93  resolution_vars[i].first,
94  cuts && interaction_chns[0].cut,
95  kNoShift,
96  weight_shape));
97  frac_resolution_cc_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
98  resolution_binning,
99  *loader,
100  resolution_vars[i].first,
101  cuts && interaction_chns[1].cut,
102  kNoShift,
103  weight_shape));
104  frac_resolution_qe_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
105  resolution_binning,
106  *loader,
107  resolution_vars[i].first,
108  cuts && interaction_chns[2].cut,
109  kNoShift,
110  weight_shape));
111  frac_resolution_res_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
112  resolution_binning,
113  *loader,
114  resolution_vars[i].first,
115  cuts && interaction_chns[3].cut,
116  kNoShift,
117  weight_shape));
118  frac_resolution_dis_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
119  resolution_binning,
120  *loader,
121  resolution_vars[i].first,
122  cuts && interaction_chns[4].cut,
123  kNoShift,
124  weight_shape));
125  frac_resolution_coh_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
126  resolution_binning,
127  *loader,
128  resolution_vars[i].first,
129  cuts && interaction_chns[5].cut,
130  kNoShift,
131  weight_shape));
132  frac_resolution_mec_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
133  resolution_binning,
134  *loader,
135  resolution_vars[i].first,
136  cuts && interaction_chns[6].cut,
137  kNoShift,
138  weight_shape));
139  frac_resolution_other_shape.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
140  resolution_binning,
141  *loader,
142  resolution_vars[i].first,
143  cuts && interaction_chns[7].cut,
144  kNoShift,
145  weight_shape));
146  //
147  frac_resolution_nocut_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
148  resolution_binning,
149  *loader,
150  resolution_vars[i].first,
151  cuts && interaction_chns[0].cut,
152  kNoShift,
153  weight_flat));
154  frac_resolution_cc_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
155  resolution_binning,
156  *loader,
157  resolution_vars[i].first,
158  cuts && interaction_chns[1].cut,
159  kNoShift,
160  weight_flat));
161  frac_resolution_qe_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
162  resolution_binning,
163  *loader,
164  resolution_vars[i].first,
165  cuts && interaction_chns[2].cut,
166  kNoShift,
167  weight_flat));
168  frac_resolution_res_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
169  resolution_binning,
170  *loader,
171  resolution_vars[i].first,
172  cuts && interaction_chns[3].cut,
173  kNoShift,
174  weight_flat));
175  frac_resolution_dis_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
176  resolution_binning,
177  *loader,
178  resolution_vars[i].first,
179  cuts && interaction_chns[4].cut,
180  kNoShift,
181  weight_flat));
182  frac_resolution_coh_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
183  resolution_binning,
184  *loader,
185  resolution_vars[i].first,
186  cuts && interaction_chns[5].cut,
187  kNoShift,
188  weight_flat));
189  frac_resolution_mec_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
190  resolution_binning,
191  *loader,
192  resolution_vars[i].first,
193  cuts && interaction_chns[6].cut,
194  kNoShift,
195  weight_flat));
196  frac_resolution_other_flat.push_back( new Spectrum("(RecoE - TrueE) / TrueE",
197  resolution_binning,
198  *loader,
199  resolution_vars[i].first,
200  cuts && interaction_chns[7].cut,
201  kNoShift,
202  weight_flat));
203 
204 
205 
206 
207  resolution_scan_trueE_flat.push_back( new Spectrum("TrueE (flat) [GeV]",
208  "(RecoE - TrueE) / TrueE",
209  *loader,
210  energy_binning,
211  kTrueE,
212  resolution_binning,
213  resolution_vars[i].first,
214  cuts,
215  kNoShift,
216  weight_flat));
217 
218  resolution_scan_EME_flat.push_back( new Spectrum("EM Shower Energy (flat) [GeV]",
219  "(RecoE - TrueE) / TrueE",
220  *loader,
221  energy_binning,
222  kRecoEME,
223  resolution_binning,
224  resolution_vars[i].first,
225  cuts,
226  kNoShift,
227  weight_flat));
228 
229  resolution_scan_HADE_flat.push_back( new Spectrum("Hadronic Energy (flat) [GeV]",
230  "(RecoE - TrueE) / TrueE",
231  *loader,
232  energy_binning,
233  kRecoHADE,
234  resolution_binning,
235  resolution_vars[i].first,
236  cuts,
237  kNoShift,
238  weight_flat));
239 
240  resolution_scan_recoE_flat.push_back( new Spectrum("RecoE (flat) [GeV]",
241  "(RecoE - TrueE) / TrueE",
242  *loader,
243  energy_binning,
244  resolution_vars[i].first,
245  resolution_binning,
246  resolution_vars[i].first,
247  cuts,
248  kNoShift,
249  weight_flat));
250 
251 
252  resolution_scan_trueE_shape.push_back( new Spectrum("TrueE (shape) [GeV]",
253  "(RecoE - TrueE) / TrueE",
254  *loader,
255  energy_binning,
256  kTrueE,
257  resolution_binning,
258  resolution_vars[i].first,
259  cuts,
260  kNoShift,
261  weight_shape));
262 
263  resolution_scan_EME_shape.push_back( new Spectrum("EM Shower Energy (shape) [GeV]",
264  "(RecoE - TrueE) / TrueE",
265  *loader,
266  energy_binning,
267  kRecoEME,
268  resolution_binning,
269  resolution_vars[i].first,
270  cuts,
271  kNoShift,
272  weight_shape));
273 
274  resolution_scan_HADE_shape.push_back( new Spectrum("Hadronic Energy (shape) [GeV]",
275  "(RecoE - TrueE) / TrueE",
276  *loader,
277  energy_binning,
278  kRecoHADE,
279  resolution_binning,
280  resolution_vars[i].first,
281  cuts,
282  kNoShift,
283  weight_shape));
284 
285  resolution_scan_recoE_shape.push_back( new Spectrum("RecoE (shape) [GeV]",
286  "(RecoE - TrueE) / TrueE",
287  *loader,
288  energy_binning,
289  resolution_vars[i].first,
290  resolution_binning,
291  resolution_vars[i].first,
292  cuts,
293  kNoShift,
294  weight_shape));
295  }
296  loader->Go();
297 
298  TFile* out_file = new TFile("energy_resolution_spectra.root", "recreate");
299  for(int i = 0; i < (int) resolution_vars.size(); i++) {
300  frac_resolution[i]->SaveTo(out_file, ("res_" + resolution_vars[i].second).c_str());
301 
302  frac_resolution_nocut_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[0].name + "_" + resolution_vars[i].second + "_shape").c_str());
303  frac_resolution_cc_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[1].name + "_" + resolution_vars[i].second + "_shape").c_str());
304  frac_resolution_qe_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[2].name + "_" + resolution_vars[i].second + "_shape").c_str());
305  frac_resolution_res_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[3].name + "_" + resolution_vars[i].second + "_shape").c_str());
306  frac_resolution_dis_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[4].name + "_" + resolution_vars[i].second + "_shape").c_str());
307  frac_resolution_coh_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[5].name + "_" + resolution_vars[i].second + "_shape").c_str());
308  frac_resolution_mec_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[6].name + "_" + resolution_vars[i].second + "_shape").c_str());
309  frac_resolution_other_shape[i]->SaveTo(out_file, ("res_" + interaction_chns[7].name + "_" + resolution_vars[i].second + "_shape").c_str());
310 
311  frac_resolution_nocut_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[0].name + "_" + resolution_vars[i].second + "_flat").c_str());
312  frac_resolution_cc_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[1].name + "_" + resolution_vars[i].second + "_flat").c_str());
313  frac_resolution_qe_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[2].name + "_" + resolution_vars[i].second + "_flat").c_str());
314  frac_resolution_res_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[3].name + "_" + resolution_vars[i].second + "_flat").c_str());
315  frac_resolution_dis_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[4].name + "_" + resolution_vars[i].second + "_flat").c_str());
316  frac_resolution_coh_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[5].name + "_" + resolution_vars[i].second + "_flat").c_str());
317  frac_resolution_mec_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[6].name + "_" + resolution_vars[i].second + "_flat").c_str());
318  frac_resolution_other_flat[i]->SaveTo(out_file, ("res_" + interaction_chns[7].name + "_" + resolution_vars[i].second + "_flat").c_str());
319 
320  resolution_scan_trueE_flat[i]->SaveTo(out_file, ("res_scan_trueE" + resolution_vars[i].second + "_flat").c_str());
321  resolution_scan_EME_flat[i]-> SaveTo(out_file, ("res_scan_EME" + resolution_vars[i].second + "_flat").c_str());
322  resolution_scan_HADE_flat[i]-> SaveTo(out_file, ("res_scan_HADE" + resolution_vars[i].second + "_flat").c_str());
323  resolution_scan_recoE_flat[i]->SaveTo(out_file, ("res_scan_recoE" + resolution_vars[i].second + "_flat").c_str());
324 
325  resolution_scan_trueE_shape[i]->SaveTo(out_file, ("res_scan_trueE" + resolution_vars[i].second + "_shape").c_str());
326  resolution_scan_EME_shape[i]-> SaveTo(out_file, ("res_scan_EME" + resolution_vars[i].second + "_shape").c_str());
327  resolution_scan_HADE_shape[i]-> SaveTo(out_file, ("res_scan_HADE" + resolution_vars[i].second + "_shape").c_str());
328  resolution_scan_recoE_shape[i]->SaveTo(out_file, ("res_scan_recoE" + resolution_vars[i].second + "_shape").c_str());
329  }
330 
331 }
const XML_Char * name
Definition: expat.h:151
const Var kFracEnergyRes_nue2018_shifted
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kFracEnergyRes_poly2X_shifted
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Cut kcPresel
Definition: NueCCIncCuts.h:139
std::string training_data
const Binning energy_binning
const Var kFracEnergyRes_poly3_shifted
const mHistAxisDef resolution_vars[kvResolutionVars]
Definition: NueCCIncVars.h:455
const Var kFlatWeightNDSig([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return 0.0;double trueE=sr->mc.nu[0].E;double index=trueE/0.2;if(trueE > 10) return 0.0;int i=(int) index;if(trueE/0.2==i) return flat_weight[i];return flat_weight[i+1];})
const Binning resolution_binning
Cut cuts
virtual void Go()=0
Load all the registered spectra.
out_file
Append EOF lines.
Definition: modifyFHiCL.py:113
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
const SelDef interaction_chns[kNumIntChns]
void getResolutionSpectra()
const Var kRecoEME([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;double CVNem_CalE=0.0;const caf::SRFuzzyKProngProxy &prim_png=sr->vtx.elastic.fuzzyk.png[0];for(const caf::SRFuzzyKProngProxy &png:sr->vtx.elastic.fuzzyk.png){double png_CalE=png.shwlid.calE;double emPID=((double) png.cvnpart.photonid+ (double) png.cvnpart.pizeroid+ (double) png.cvnpart.electronid);double haPID=((double) png.cvnpart.protonid+ (double) png.cvnpart.pionid+ (double) png.cvnpart.neutronid+ (double) png.cvnpart.otherid+ (double) png.cvnpart.muonid);if(emPID< 0) continue;if(emPID >=haPID ) CVNem_CalE+=png_CalE;else continue;}if(CVNem_CalE==0.0) CVNem_CalE=prim_png.shwlid.calE;return CVNem_CalE *CalibrationBugCorrectionFactor(sr->hdr);})
Definition: NueEnergy2018.h:17
const Var kFracEnergyRes_poly2
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:527
const Var kFracEnergyRes_poly2_shifted
const Var kFracEnergyRes_nue2018
const SystShifts kNoShift
Definition: SystShifts.h:115
const Var kFracEnergyRes_poly2X
Base class for the various types of spectrum loader.
const Cut kIsTrueSig
Var weight_shape
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Var weight_flat
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const Var kFracEnergyRes_poly3
const Var kFracEnergyRes_poly3X
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Var kFracEnergyRes_poly3X_shifted
const Var kRecoHADE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1.0;if(sr->vtx.elastic.fuzzyk.npng< 1) return-1.0;const double CVNha_CalE=sr->slc.calE *CalibrationBugCorrectionFactor(sr->hdr);return std::max(CVNha_CalE-kRecoEME(sr), 0.);})
Definition: NueEnergy2018.h:18