Functions
neutronE_macro.C File Reference
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "Utilities/rootlogon.C"
#include "TCanvas.h"
#include "TH2.h"
#include "TProfile.h"
#include "TSystem.h"

Go to the source code of this file.

Functions

void neutronE_macro (bool IsFHC=false, bool IsFarDet=true)
 

Function Documentation

void neutronE_macro ( bool  IsFHC = false,
bool  IsFarDet = true 
)

Definition at line 33 of file neutronE_macro.C.

References CutNames, fillBadChanDBTables::det, E, caf::Proxy< caf::SRVertexBranch >::elastic, plot_validation_datamc::fname, for, caf::Proxy< caf::SRElastic >::fuzzyk, GENIEStr, ana::SpectrumLoader::Go(), if(), inFile, IsFHC, caf::Proxy< caf::SRElastic >::IsValid, ana::kHadEFracAxis, ana::kNeutrinoEn(), ana::kNeutronEn(), ana::kNNeutrons(), ana::kNoCut, ana::kNoShift, ana::kNue2018FD, ana::kNue2018NDCVNSsb, ana::kNumuCCOptimisedAxis, ana::kNumuContainFD2017, ana::kNumuContainND2017, ana::kNumuCosmicRej2018, ana::kNumuPID2018, ana::kNumuQuality, ana::kPPFXFluxCVWgt, ana::kReMId(), ana::kRemidBinning, ana::kStandardDQCuts, ana::kXSecCVWgt2018, demo0::loader, ana::MassNeut, ana::MassProt, caf::Proxy< caf::StandardRecord >::mc, caf::Proxy< caf::SRTruthBranch >::nnu, caf::Proxy< caf::SRFuzzyK >::npng, caf::Proxy< caf::SRTruthBranch >::nu, OutFile, caf::Proxy< caf::SRRemid >::pid, ana::pnfs2xrootd(), caf::Proxy< caf::SRFuzzyK >::png, caf::Proxy< caf::SRFuzzyK >::png2d, ana::QuantileCutsFromTH2(), QuantNames, caf::Proxy< caf::SRIDBranch >::remid, ana::SaveTo(), caf::Proxy< caf::StandardRecord >::sel, ana::SpectrumLoaderBase::SetSpillCut(), ana::Binning::Simple(), sr, string, caf::Proxy< caf::StandardRecord >::vtx, and ana::weight.

34 {
35 
36  std::string sFHC_RHC = ( IsFHC == true ? "fhc" : "rhc" );
37  std::string sFD_ND = ( IsFarDet == true ? "FD" : "ND" );
38 
39  // --- Define my loader, and the file / dataset which it will load...
40  std::string fname = "";
41  if ( IsFHC ) {
42  if (IsFarDet) {
43 // fname = "prod_sumdecaf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
44  fname = "prod_caf_R17-11-14-prod4reco.d_fd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
45  }
46  else {
47 // fname = "prod_sumdecaf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1_numu2018"; // numu concats
48  fname = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1"; // full cafs
49  }
50  }
51  else {
52  if (IsFarDet) {
53 // fname = "prod_sumdecaf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
54  fname = "prod_caf_R17-11-14-prod4reco.e_fd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
55  }
56  else {
57 // fname = "prod_sumdecaf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1_numu2018"; // numu concats
58  fname = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1"; // full cafs
59  }
60  }
61 
62  SpectrumLoader loader(fname);
63  loader.SetSpillCut(kStandardDQCuts);
64 
65  // weights, shifts
66  const Var weight = kPPFXFluxCVWgt * kXSecCVWgt2018;
67  const SystShifts shift = kNoShift;
68 
69  // --- Define my output directory.
70  std::string sOutFile = "SpectrumFiles_"+sFD_ND+"_"+sFHC_RHC+".root";
71 
72  // --- Define the number of GENIE interaction types I'm looking at....
73  const std::vector<Cut> GENIECuts = { kNoCut };
74  const std::vector<std::string> GENIEStr = { "All" };
75  //const std::vector<Cut> GENIECuts = { kNoCut, kIsDytmanMEC, kIsDIS, kIsRes, kIsCoh, kIsQE };
76  //const std::vector<std::string> GENIEStr = { "All", "MEC", "DIS", "Resonance", "Coherent", "QuasiElactic" };
77  const unsigned int NumGENIE = GENIECuts.size();
78 
79  // --- Define my broad cuts
80  std::vector<std::string> CutNames; std::vector<Cut> kDetCuts;
81  if (IsFarDet) {
82  CutNames.push_back( "FD_2018-PID" ); kDetCuts.push_back( kNumuQuality && kNumuContainFD2017 && kNumuCosmicRej2018);
83  CutNames.push_back( "FD_2018-Cont" ); kDetCuts.push_back( kNumuQuality && kNumuPID2018 && kNumuCosmicRej2018 );
84  CutNames.push_back( "FD_2018!PID" ); kDetCuts.push_back( kNumuQuality && kNumuContainFD2017 && !kNumuPID2018 && kNumuCosmicRej2018 );
85  CutNames.push_back( "FD_2018!Cont" ); kDetCuts.push_back( kNumuQuality && !kNumuContainFD2017 && kNumuPID2018 && kNumuCosmicRej2018 );
86  CutNames.push_back( "FD_2018full" ); kDetCuts.push_back( kNumuQuality && kNumuContainFD2017 && kNumuPID2018 && kNumuCosmicRej2018 );
87  CutNames.push_back( "FD_nue2018full"); kDetCuts.push_back( kNue2018FD );
88  } else {
89  CutNames.push_back( "ND_2018-PID" ); kDetCuts.push_back( kNumuQuality && kNumuContainND2017);
90  CutNames.push_back( "ND_2018-Cont" ); kDetCuts.push_back( kNumuQuality && kNumuPID2018 );
91  CutNames.push_back( "ND_2018!PID" ); kDetCuts.push_back( kNumuQuality && kNumuContainND2017 && !kNumuPID2018 );
92  CutNames.push_back( "ND_2018!Cont" ); kDetCuts.push_back( kNumuQuality && !kNumuContainND2017 && kNumuPID2018 );
93  CutNames.push_back( "ND_2018full" ); kDetCuts.push_back( kNumuQuality && kNumuContainND2017 && kNumuPID2018 );
94  CutNames.push_back( "ND_nue2018full"); kDetCuts.push_back( kNue2018NDCVNSsb );
95  }
96  unsigned int nDetCuts = kDetCuts.size();
97 
98  // --- Figure out where to load my quantile boundaries from.
99  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__" + sFHC_RHC + "_full__numu2018.root";
100 
101  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile.c_str()).c_str() );
102  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
103  // How many quantiles?
104  const int NHadEFracQuantiles = 4; // defines how many divisions of hadEFrac are used
105  // Make my cuts.
106  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
107  HadEFracQuantCuts.push_back( kNoCut );
108  std::vector<std::string> QuantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };
109  const unsigned int nQuants = QuantNames.size();
110 
111  // --- Define some bin sets which I want to use....
112  const unsigned int NumEBins = 100;
113  const Binning LepBins = Binning::Simple(500, 0., 5.);
114  const Binning HadBins = Binning::Simple(100, 0., 1.);
115  const Binning FracBins = Binning::Simple(100, 0., 1.);
116  const Binning EpHBins = Binning::Simple(50, 0., 100.);
117  const Binning CVNBins = Binning::Simple(20, 0., 1.);
118 
119  static const double MassNeut = 0.939565;
120  static const double MassProt = 0.938272;
121 
122  // ----- Neutrino energy
123  const Var kNeutrinoEn([](const caf::SRProxy *sr) {
124  double ThisEn = 0.;
125  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
126  if (sr->mc.nu[0].E<=0) return -5.; // avoid nan in energy fraction
127 
128  // --- Check I have some truth information, and if do calc energy
129  ThisEn = sr->mc.nu[0].E;
130  // --- Now return our final value...
131  return ThisEn;
132  });
133 
134  // ----- Total energy from primaries
135  const Var kTotPrimE([](const caf::SRProxy* sr) {
136  double ThisEn = 0.;
137  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
138  if (sr->mc.nu[0].E<=0) return -5.; // avoid nan in energy fraction
139 
140  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
141  double E = 0;
142  if(sr->mc.nu[0].prim[PrimL].p.E<0.0) continue;
143 
144  if (sr->mc.nu[0].prim[PrimL].pdg == 2112) E = sr->mc.nu[0].prim[PrimL].p.E - MassNeut;
145  else if (sr->mc.nu[0].prim[PrimL].pdg == 2212) E = sr->mc.nu[0].prim[PrimL].p.E - MassProt;
146  else E = sr->mc.nu[0].prim[PrimL].p.E;
147 
148  ThisEn += E;
149  }
150  return ThisEn;
151  });
152 
153  // ----- Total energy from primaries and neutron KE +20%
154  const Var kTotPrimEnKEp20([](const caf::SRProxy* sr) {
155  double ThisEn = 0.;
156  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
157  if (sr->mc.nu[0].E<=0) return -5.; // avoid nan in energy fraction
158 
159  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
160  double E = 0;
161  if(sr->mc.nu[0].prim[PrimL].p.E<0.0) continue;
162 
163  if (sr->mc.nu[0].prim[PrimL].pdg == 2112) E = 1.2*(sr->mc.nu[0].prim[PrimL].p.E - MassNeut);
164  else if (sr->mc.nu[0].prim[PrimL].pdg == 2212) E = sr->mc.nu[0].prim[PrimL].p.E - MassProt;
165  else E = sr->mc.nu[0].prim[PrimL].p.E;
166 
167  ThisEn += E;
168  }
169  return ThisEn;
170  });
171 
172  // ----- Total energy from primaries and neutron KE -20%
173  const Var kTotPrimEnKEm20([](const caf::SRProxy* sr) {
174  double ThisEn = 0.;
175  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
176  if (sr->mc.nu[0].E<=0) return -5.; // avoid nan in energy fraction
177 
178 
179  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
180  double E = 0;
181  if(sr->mc.nu[0].prim[PrimL].p.E<0.0) continue;
182 
183  if (sr->mc.nu[0].prim[PrimL].pdg == 2112) E = 0.8*(sr->mc.nu[0].prim[PrimL].p.E - MassNeut);
184  else if (sr->mc.nu[0].prim[PrimL].pdg == 2212) E = sr->mc.nu[0].prim[PrimL].p.E - MassProt;
185  else E = sr->mc.nu[0].prim[PrimL].p.E;
186 
187  ThisEn += E;
188  }
189  return ThisEn;
190  });
191  // ----- Neutron multiplicity
192  const Var kNNeutrons([](const caf::SRProxy* sr){
193  if (sr->mc.nnu == 0) return -5.; //if (sr->mc.nu.size() == 0) return false;
194  return double(sr->mc.nu[0].nneutron);
195  });
196 
197  // ----- Neutron energy
198  const MultiVar kNeutronEn([](const caf::SRProxy* sr) {
199 
200  std::vector<double> E;
201  if (sr->mc.nnu == 0) return E; //if (sr->mc.nu.size() == 0) return false;
202  if (sr->mc.nu[0].E<=0) return E; // avoid nan in energy fraction
203 
204  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
205  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
206  E.push_back(sr->mc.nu[0].prim[PrimL].p.E - MassNeut);
207  }
208  return E;
209  });
210 
211  // ----- Neutron energy fraction
212  const MultiVar kNeutronEF([](const caf::SRProxy* sr) {
213 
214  std::vector<double> E;
215  if (sr->mc.nnu == 0) return E; //if (sr->mc.nu.size() == 0) return false;
216  if (sr->mc.nu[0].E<=0) return E; // avoid nan in energy fraction
217 
218  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
219  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
220  E.push_back((sr->mc.nu[0].prim[PrimL].p.E - MassNeut)/sr->mc.nu[0].E);
221  }
222  return E;
223  });
224 
225  // ----- Neutron prong calE
226  const MultiVar kNeutronProngCalE([](const caf::SRProxy* sr) {
227 
228  std::vector<double> E;
229 
230  if(!sr->vtx.elastic.IsValid) return E;
231  if(sr->vtx.elastic.fuzzyk.npng == 0) return E;
232 
233  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
234  if (png.truth.motherpdg != 2112) continue;
235  E.push_back(png.calE);
236  }
237  return E;
238  });
239 
240  // ----- Neutron prong calE per hit
241  const MultiVar kNeutronProngCalEperHit([](const caf::SRProxy* sr) {
242  std::vector<double> EperH;
243  if(!sr->vtx.elastic.IsValid) return EperH;
244  if(sr->vtx.elastic.fuzzyk.npng == 0) return EperH;
245 
246  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
247  if (png.truth.motherpdg != 2112) continue;
248  EperH.push_back((png.calE*1000)/double(png.nhit));
249  }
250  return EperH;
251  });
252 
253  // ----- Neutron prong calE with 2D prongs
254  const MultiVar kNeutronProng2dCalE([](const caf::SRProxy* sr) {
255 
256  std::vector<double> E;
257  if( !sr->vtx.elastic.IsValid ) return E;
258  if(sr->vtx.elastic.fuzzyk.npng == 0) return E;
259 
260  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
261  if (png.truth.motherpdg != 2112) continue;
262  E.push_back(png.calE);
263  }
264 
265  for(const auto& png : sr->vtx.elastic.fuzzyk.png2d){
266  if (png.truth.motherpdg != 2112) continue;
267  E.push_back(png.calE);
268  }
269 
270  return E;
271  });
272 
273  // ----- Neutron prong calE per hit with 2D prongs
274  const MultiVar kNeutronProng2dCalEperHit([](const caf::SRProxy* sr) {
275 
276  std::vector<double> EperH;
277 
278  if( !sr->vtx.elastic.IsValid) return EperH;
279  if(sr->vtx.elastic.fuzzyk.npng == 0) return EperH;
280 
281  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
282  if (png.truth.motherpdg != 2112) continue;
283  EperH.push_back((png.calE*1000)/double(png.nhit));
284  }
285 
286  for(const auto& png : sr->vtx.elastic.fuzzyk.png2d){
287  if (png.truth.motherpdg != 2112) continue;
288  EperH.push_back((png.calE*1000)/double(png.nhit));
289  }
290 
291  return EperH;
292  });
293 
294  // ----- Neutron prong weighted calE
295  const MultiVar kNeutronProngWtCalE([](const caf::SRProxy* sr) {
296 
297  std::vector<double> E;
298 
299  if( !sr->vtx.elastic.IsValid ) return E;
300  if(sr->vtx.elastic.fuzzyk.npng == 0) return E;
301 
302  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
303  if (png.truth.motherpdg != 2112) continue;
304  E.push_back(png.weightedCalE);
305  }
306  return E;
307  });
308 
309  // ----- Neutron prong weighted calE per hit
310  const MultiVar kNeutronProngWtCalEperHit([](const caf::SRProxy* sr) {
311 
312  std::vector<double> EperH;
313 
314  if( !sr->vtx.elastic.IsValid ) return EperH;
315  if(sr->vtx.elastic.fuzzyk.npng == 0) return EperH;
316 
317  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
318  if (png.truth.motherpdg != 2112) continue;
319  EperH.push_back((png.weightedCalE*1000)/double(png.nhit));
320  }
321  return EperH;
322  });
323 
324  // ----- Neutron prong weighted calE with 2D prongs
325  const MultiVar kNeutronProngWt2dCalE([](const caf::SRProxy* sr) {
326 
327 
328  std::vector<double> E;
329  if( !sr->vtx.elastic.IsValid ) return E;
330  if(sr->vtx.elastic.fuzzyk.npng == 0) return E;
331 
332  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
333  if (png.truth.motherpdg != 2112) continue;
334  E.push_back(png.weightedCalE);
335  }
336 
337  for(const auto& png : sr->vtx.elastic.fuzzyk.png2d){
338  if (png.truth.motherpdg != 2112) continue;
339  E.push_back(png.weightedCalE);
340  }
341 
342  return E;
343  });
344 
345  // ----- Neutron prong weighted calE per hit with 2D prongs
346  const MultiVar kNeutronProngWt2dCalEperHit([](const caf::SRProxy* sr) {
347 
348  std::vector<double> EperH;
349  if(!sr->vtx.elastic.IsValid) return EperH;
350  if(sr->vtx.elastic.fuzzyk.npng == 0) return EperH;
351 
352  for(const auto& png : sr->vtx.elastic.fuzzyk.png){
353  if (png.truth.motherpdg != 2112) continue;
354  EperH.push_back((png.weightedCalE*1000)/double(png.nhit));
355  }
356 
357  for(const auto& png : sr->vtx.elastic.fuzzyk.png2d){
358  if (png.truth.motherpdg != 2112) continue;
359  EperH.push_back((png.weightedCalE*1000)/double(png.nhit));
360  }
361 
362  return EperH;
363  });
364 
365 
366  // need multivars for these or I can't plot them against neutron multivars
367  // ----- ReMId
368  const MultiVar kReMId([](const caf::SRProxy* sr) {
369 
370  std::vector<double> id;
371  if (sr->mc.nnu == 0) return id; //if (sr->mc.nu.size() == 0) return false;
372  if (sr->mc.nu[0].E<=0) return id; // avoid nan in energy fraction
373 
374  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
375  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
376  id.push_back(sr->sel.remid.pid);
377  }
378  return id;
379  });
380 
381  // ----- CVNid
382  const MultiVar kCVNeid([](const caf::SRProxy* sr){
383 
384  std::vector<double> id;
385 
386  if (sr->mc.nnu == 0) return id; //if (sr->mc.nu.size() == 0) return false;
387  if (sr->mc.nu[0].E<=0) return id; // avoid nan in energy fraction
388 
389  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
390  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
391  id.push_back(sr->sel.cvnProd3Train.nueid);
392  }
393  return id;
394  });
395 
396 
397  const MultiVar kCVNmid([](const caf::SRProxy* sr){
398 
399  std::vector<double> id;
400 
401  if (sr->mc.nnu == 0) return id; //if (sr->mc.nu.size() == 0) return false;
402  if (sr->mc.nu[0].E<=0) return id; // avoid nan in energy fraction
403 
404  for (unsigned int PrimL=0; PrimL < sr->mc.nu[0].prim.size(); ++PrimL) {
405  if (sr->mc.nu[0].prim[PrimL].pdg != 2112) continue;
406  id.push_back(sr->sel.cvnProd3Train.numuid);
407  }
408  return id;
409  });
410 
411  // --- 1D Spectra
412  // Energy plots
413  Spectrum *NeutrinoE [NumGENIE][nDetCuts][nQuants];
414  Spectrum *TotPrimE [NumGENIE][nDetCuts][nQuants];
415  Spectrum *TotPrimEnKEp20 [NumGENIE][nDetCuts][nQuants];
416  Spectrum *TotPrimEnKEm20 [NumGENIE][nDetCuts][nQuants];
417  Spectrum *NeutMult [NumGENIE][nDetCuts][nQuants];
418  Spectrum *NeutEnergy [NumGENIE][nDetCuts][nQuants];
419  Spectrum *NeutEFrac [NumGENIE][nDetCuts][nQuants];
420 
421  Spectrum *NeutPngCalE [NumGENIE][nDetCuts][nQuants];
422  Spectrum *NeutPngCalEpH [NumGENIE][nDetCuts][nQuants];
423  Spectrum *NeutPng2DCalE [NumGENIE][nDetCuts][nQuants];
424  Spectrum *NeutPng2DCalEpH [NumGENIE][nDetCuts][nQuants];
425 
426  Spectrum *NeutPngWtCalE [NumGENIE][nDetCuts][nQuants];
427  Spectrum *NeutPngWtCalEpH [NumGENIE][nDetCuts][nQuants];
428  Spectrum *NeutPng2DWtCalE [NumGENIE][nDetCuts][nQuants];
429  Spectrum *NeutPng2DWtCalEpH[NumGENIE][nDetCuts][nQuants];
430 
431  // 2d spectra
432 
433  Spectrum *NEFvsCVNe [NumGENIE][nDetCuts][nQuants];
434  Spectrum *NEFvsCVNm [NumGENIE][nDetCuts][nQuants];
435  Spectrum *NEFvsReMId[NumGENIE][nDetCuts][nQuants];
436 
437  // --- Define my Spectra which are to be filled from the loader
438  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
439  for (unsigned int det=0; det<nDetCuts; ++det) {
440  for (unsigned int quant=0; quant<nQuants; ++quant) {
441  // Define my new cut.
442  Cut GENIECut = GENIECuts[gen] && kDetCuts[det] && HadEFracQuantCuts[quant];
443  std::string MyStr = " for GENIE Int. "+GENIEStr[gen]+", Cut "+CutNames[det]+", "+QuantNames[quant];
444  // --- Now fill our 1D spectra!
445  // Energy spectra
446  NeutrinoE[gen][det][quant] = new Spectrum( "Neutrino energy (GeV)" + MyStr, LepBins, loader, kNeutrinoEn, GENIECut, kNoShift, weight);
447  TotPrimE[gen][det][quant] = new Spectrum( "Total primE (GeV)" + MyStr, LepBins, loader, kTotPrimE, GENIECut, kNoShift, weight);
448  TotPrimEnKEp20[gen][det][quant] = new Spectrum( "Total primE (GeV), n KE +20%"+ MyStr, LepBins, loader, kTotPrimEnKEp20, GENIECut, kNoShift, weight);
449  TotPrimEnKEm20[gen][det][quant] = new Spectrum( "Total primE (GeV), n KE +20%"+ MyStr, LepBins, loader, kTotPrimEnKEm20, GENIECut, kNoShift, weight);
450  NeutMult [gen][det][quant] = new Spectrum( "Neutron multiplicity" +MyStr, Binning::Simple(5,0,5), loader, kNNeutrons, GENIECut, kNoShift, weight);
451  NeutEnergy [gen][det][quant] = new Spectrum( "True Neutron KE (MeV)" +MyStr, HadBins, loader, kNeutronEn, GENIECut, kNoShift, weight);
452  NeutEFrac [gen][det][quant] = new Spectrum( "True Neutron KE / True Neutrino E" +MyStr, FracBins, loader, kNeutronEF, GENIECut, kNoShift, weight);
453  NeutPngCalE [gen][det][quant] = new Spectrum( "Neutron Prong calE"+MyStr, HadBins, loader, kNeutronProngCalE, GENIECut, kNoShift, weight);
454  NeutPngCalEpH [gen][det][quant] = new Spectrum( "Neutron Prong calE per hit"+MyStr, EpHBins, loader, kNeutronProngCalEperHit, GENIECut, kNoShift, weight);
455  NeutPng2DCalE [gen][det][quant] = new Spectrum( "Neutron Prong (inc. 2d) calE" +MyStr, HadBins, loader, kNeutronProng2dCalE, GENIECut, kNoShift, weight);
456  NeutPng2DCalEpH [gen][det][quant] = new Spectrum( "Neutron Prong (inc. 2d) calE per hit" +MyStr, EpHBins, loader, kNeutronProng2dCalEperHit, GENIECut, kNoShift, weight);
457  NeutPngWtCalE [gen][det][quant] = new Spectrum( "Neutron Prong weighted calE"+MyStr, HadBins, loader, kNeutronProngWtCalE, GENIECut, kNoShift, weight);
458  NeutPngWtCalEpH [gen][det][quant] = new Spectrum( "Neutron Prong weighted calE per hit" +MyStr, EpHBins, loader, kNeutronProngWtCalEperHit, GENIECut, kNoShift, weight);
459  NeutPng2DWtCalE [gen][det][quant] = new Spectrum( "Neutron Prong (inc. 2d) weighted calE"+MyStr, HadBins, loader, kNeutronProngWt2dCalE, GENIECut, kNoShift, weight);
460  NeutPng2DWtCalEpH [gen][det][quant] = new Spectrum( "Neutron Prong (inc. 2d) weighted calE per hit"+MyStr, HadBins, loader, kNeutronProngWt2dCalEperHit, GENIECut, kNoShift, weight);
461  NEFvsCVNe [gen][det][quant] = new Spectrum( "Neutron KE/nu E vs. CVNe"+MyStr, loader, CVNBins, kCVNeid, FracBins, kNeutronEF, GENIECut, kNoShift, weight);
462  NEFvsCVNm [gen][det][quant] = new Spectrum( "Neutron KE/nu E vs. CVNm"+MyStr, loader, CVNBins, kCVNmid, FracBins, kNeutronEF, GENIECut, kNoShift, weight);
463  NEFvsReMId [gen][det][quant] = new Spectrum( "Neutron KE/nu E vs. ReMId"+MyStr, loader, kRemidBinning, kReMId, FracBins, kNeutronEF, GENIECut, kNoShift, weight);
464 
465  } // Loop through my Spectrum arrays.
466  }
467  }
468 
469  // --- Do it!
470  loader.Go();
471 
472  // --- Open my outfile so that I save my histograms!
473  TFile *OutFile = new TFile(sOutFile.c_str(),"RECREATE");
474 
475  // --- And finally, lets make my plots a little nicer and save the fuckers!
476  for (unsigned int gen=0; gen<NumGENIE; ++gen) {
477  for (unsigned int det=0; det<nDetCuts; ++det) {
478  for (unsigned int quant=0; quant<nQuants; ++quant) {
479  std::string MyStr = GENIEStr[gen]+"_"+CutNames[det]+"_"+QuantNames[quant];;
480  // --- Save my 1D histograms!
481  // Energy spectra.
482  NeutrinoE [gen][det][quant] -> SaveTo( OutFile, TString("NeutrinoE_") + TString(MyStr) ) ;
483  TotPrimE [gen][det][quant] -> SaveTo( OutFile, TString("TotPrimE_") + TString(MyStr) ) ;
484  TotPrimEnKEp20 [gen][det][quant] -> SaveTo( OutFile, TString("TotPrimEnKEp20_") + TString(MyStr) ) ;
485  TotPrimEnKEm20 [gen][det][quant] -> SaveTo( OutFile, TString("TotPrimEnKEm20_") + TString(MyStr) ) ;
486  NeutMult [gen][det][quant] -> SaveTo( OutFile, TString("NeutMult") + TString(MyStr) ) ;
487  NeutEnergy [gen][det][quant] -> SaveTo( OutFile, TString("NeutEnergy") + TString(MyStr) ) ;
488  NeutEFrac [gen][det][quant] -> SaveTo( OutFile, TString("NeutEFrac") + TString(MyStr) ) ;
489  NeutPngCalE [gen][det][quant] -> SaveTo( OutFile, TString("NeutPngCalE") + TString(MyStr) ) ;
490  NeutPngCalEpH [gen][det][quant] -> SaveTo( OutFile, TString("NeutPngCalEpH") + TString(MyStr) ) ;
491  NeutPng2DCalE [gen][det][quant] -> SaveTo( OutFile, TString("NeutPng2DCalE") + TString(MyStr) ) ;
492  NeutPng2DCalEpH [gen][det][quant] -> SaveTo( OutFile, TString("NeutPng2DCalEpH") + TString(MyStr) ) ;
493  NeutPngWtCalE [gen][det][quant] -> SaveTo( OutFile, TString("NeutPngWtCalE") + TString(MyStr) ) ;
494  NeutPngWtCalEpH [gen][det][quant] -> SaveTo( OutFile, TString("NeutPngWtCalEpH") + TString(MyStr) ) ;
495  NeutPng2DWtCalE [gen][det][quant] -> SaveTo( OutFile, TString("NeutPng2DWtCalE") + TString(MyStr) ) ;
496  NeutPng2DWtCalEpH [gen][det][quant] -> SaveTo( OutFile, TString("NeutPng2DWtCalEpH") + TString(MyStr) ) ;
497  // and 2D
498  NEFvsCVNe [gen][det][quant] -> SaveTo( OutFile, TString("NEFvsCVNe_") + TString(MyStr) ) ;
499  NEFvsCVNm [gen][det][quant] -> SaveTo( OutFile, TString("NEFvsCVNm_") + TString(MyStr) ) ;
500  NEFvsReMId [gen][det][quant] -> SaveTo( OutFile, TString("NEFvsReMId_")+ TString(MyStr) ) ;
501 
502  }
503  }
504  }
505 }
xml_schema::id
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
caf::Proxy< caf::SRFuzzyK >::npng
caf::Proxy< size_t > npng
Definition: SRProxy.h:2038
plot_validation_datamc.fname
fname
Definition: plot_validation_datamc.py:180
caf::Proxy< caf::SRElastic >::fuzzyk
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
ana::Binning
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
ana::SaveTo
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
Definition: LoadFromFile.cxx:105
ana::kRemidBinning
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
caf::Proxy< caf::SRFuzzyK >::png2d
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
ana::kPPFXFluxCVWgt
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
ana::weight
const Var weight
Definition: PIDCutTuningDefs.h:96
ana::kNumuContainND2017
const Cut kNumuContainND2017([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;++i){TVector3 start=sr->vtx.elastic.fuzzyk.png[i].shwlid.start;TVector3 stop=sr->vtx.elastic.fuzzyk.png[i].shwlid.stop;if(std::min(start.X(), stop.X())< -180.0) return false;if(std::max(start.X(), stop.X()) > 180.0) return false;if(std::min(start.Y(), stop.Y())< -180.0) return false;if(std::max(start.Y(), stop.Y()) > 180.0) return false;if(std::min(start.Z(), stop.Z())< 20.0) return false;if(std::max(start.Z(), stop.Z()) > 1525.0) return false;}if(sr->trk.kalman.ntracks< 1) return false;for(unsigned int i=0;i< sr->trk.kalman.ntracks;++i){if(i==sr->trk.kalman.idxremid) continue;else if(sr->trk.kalman.tracks[i].start.Z() > 1275||sr->trk.kalman.tracks[i].stop.Z() > 1275) return false;}return(sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1100 &&(sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&sr->sel.contain.kalfwdcellnd > 5 &&sr->sel.contain.kalbakcellnd > 10);})
Definition: NumuCuts2017.h:11
IsFHC
bool IsFHC
Definition: ND_DataMC_Comp_Systs.C:4
ana::SystShifts
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
caf::Proxy< caf::StandardRecord >
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< caf::SRTruthBranch >::nu
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
ana::kReMId
const Var kReMId([](const caf::SRProxy *sr){return sr->sel.remid.pid;})
ana::kNue2018NDCVNSsb
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
ana::kHadEFracAxis
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
caf::Proxy< caf::SRRemid >::pid
caf::Proxy< float > pid
Definition: SRProxy.h:1136
caf::Proxy< caf::SRTruthBranch >::nnu
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
ana::pnfs2xrootd
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
caf::Proxy< caf::SRVertexBranch >::elastic
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
ana::kNumuCosmicRej2018
const Cut kNumuCosmicRej2018([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.53 && sr->slc.nhit< 400 && sr->sel.nuecosrej.pngptp< 0.9 );})
Definition: NumuCuts2018.h:19
caf::Proxy< caf::SRFuzzyK >::png
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
ana::kNumuContainFD2017
const Cut kNumuContainFD2017
Definition: NumuCuts2017.h:21
if
if(dump)
Definition: checkMeasured_Target.C:110
E
Float_t E
Definition: plot.C:20
ana::kStandardDQCuts
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:16
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
ana::MassProt
static const float MassProt
Definition: EnergyCont_help.h:11
kNue2018FD
const Cut kNue2018FD
Definition: make_nue_pidtuning.C:292
caf::Proxy< caf::SRIDBranch >::remid
caf::Proxy< caf::SRRemid > remid
Definition: SRProxy.h:1269
QuantNames
const std::string QuantNames[NumQuant]
Definition: DoThePlots.C:207
demo0.loader
loader
Definition: demo0.py:10
ana::kNNeutrons
const Var kNNeutrons([](const caf::SRProxy *sr){float NPar=0.;if(DEBUGGING){std::cout<< "\nLooking at event "<< sr->hdr.evt<< ", there are "<< sr->mc.nu.size()<< ", the 0th neutrino is a "<< sr->mc.nu[0].pdg<< ", with Energy "<< sr->mc.nu[0].E<< ". There were a total of "<< sr->mc.nu[0].prim.size()<< " other primaries associated with that neutrino, "<< sr->mc.nu[0].nneutron<< " were neutrons."<< std::endl;if(kIsDytmanMEC(sr)) std::cout<< " The neutrino was a MEC"<< std::endl;else if(kIsDIS(sr)) std::cout<< " The neutrino was a DIS"<< std::endl;else if(kIsRes(sr)) std::cout<< " The neutrino was a Resonance"<< std::endl;else if(kIsCoh(sr)) std::cout<< " The neutrino was a Coherent"<< std::endl;float TotEn=0.0;for(unsigned int PrimL=0;PrimL< sr->mc.nu[0].prim.size();++PrimL){std::cout<< "\t\tLooking at Prim "<< PrimL<< " of "<< sr->mc.nu[0].prim.size()<< ", it was a "<< sr->mc.nu[0].prim[PrimL].pdg<< ", VisE "<< sr->mc.nu[0].prim[PrimL].visE<< ", VisEInslc "<< sr->mc.nu[0].prim[PrimL].visEinslc<< ", E0 "<< sr->mc.nu[0].prim[PrimL].p.E<< ", enteringE "<< sr->mc.nu[0].prim[PrimL].enteringE<< ", totEscE "<< sr->mc.nu[0].prim[PrimL].totEscE<< std::endl;TotEn+=sr->mc.nu[0].prim[PrimL].totEscE;}std::cout<< "\tThe total uncontained energy was "<< TotEn<< std::endl;NPar=sr->mc.nu[0].nneutron;}return NPar;})
cmf::Spectrum
std::vector< float > Spectrum
Definition: Constants.h:759
ana::kNumuCCOptimisedAxis
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
ana::_MultiVar< caf::SRProxy >
ana::kNoShift
const SystShifts kNoShift
Definition: SystShifts.cxx:22
caf::Proxy< caf::StandardRecord >::mc
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
caf::Proxy< caf::SRElastic >::IsValid
caf::Proxy< bool > IsValid
Definition: SRProxy.h:2058
ana::SpectrumLoader
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Definition: SpectrumLoader.h:21
ana::kNumuPID2018
const Cut kNumuPID2018([](const caf::SRProxy *sr){std::cout<< "ERROR::kNumuPID2018, cutting on both cvnProd3Train and cvn2017."<< " Neither branch exists anymore. Returning False."<< std::endl;abort();return false;})
Definition: NumuCuts2018.h:22
ana::kNeutronEn
const Var kNeutronEn([](const caf::SRProxy *sr){float ThisEn=0.;ThisEn=CalcPrimE(sr, 2112, false);return ThisEn;})
ana::kXSecCVWgt2018
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
ana::MassNeut
static const float MassNeut
Definition: EnergyCont_help.h:10
GENIEStr
const std::vector< std::string > GENIEStr
Definition: EnergyCont_plot.C:30
ana::QuantileCutsFromTH2
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
Definition: QuantileCuts.cxx:198
ana::kNeutrinoEn
const Var kNeutrinoEn([](const caf::SRProxy *sr){float ThisEn=0.;ThisEn=sr->mc.nu[0].E;return ThisEn;})
OutFile
TFile * OutFile
Definition: NuMu2020_BasicPIDPlots_Plot.C:70
caf::Proxy< caf::StandardRecord >::sel
caf::Proxy< caf::SRIDBranch > sel
Definition: SRProxy.h:2141
ana::_Cut< caf::SRProxy >
caf::Proxy< caf::StandardRecord >::vtx
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
ana::kNumuQuality
const Cut kNumuQuality
Definition: NumuCuts.h:18
ana::_Var< caf::SRProxy >
for
#define for
Definition: msvc_pragmas.h:3
CutNames
std::vector< std::string > CutNames
Definition: MakeCutFlow.C:49
string
enum BeamMode string
Generated on Sat Jan 23 2021 14:14:50 for NOvA by   doxygen 1.8.11