caf_nue_data_mc.C
Go to the documentation of this file.
1 #ifdef __CINT__
2 void caf_nue_data_mc(){}
3 #else
4 
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 #include "CAFAna/Cuts/NueCutsFirstAna.h"
12 #include "CAFAna/Vars/Vars.h"
13 using namespace ana;
14 
15 #include "OscLib/func/IOscCalculator.h"
16 
17 #include "TCanvas.h"
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TLegend.h"
21 #include "TPad.h"
22 
23 #define ELASTICVAR(VAR) \
24  Var({"vtx.nelastic", "vtx.elastic.*"}, \
25  [](const caf::StandardRecord* sr) \
26  { \
27  if(sr->vtx.nelastic < 1) return -1000.; \
28  return double(sr->vtx.elastic[0].vtx.VAR); \
29  });
30 
31 const Var kVtxX = ELASTICVAR(X());
32 const Var kVtxY = ELASTICVAR(Y());
33 const Var kVtxZ = ELASTICVAR(Z());
34 
35 const Var kEPerHit({"slc.calE","slc.nhit"},
36  [](const caf::StandardRecord* sr)
37  {
38  return sr->slc.calE/sr->slc.nhit;
39  });
40 
41 #define SHWLIDVAR(VAR) \
42  Var({"shw.nshwlid", "shw.shwlid.*"}, \
43  [](const caf::StandardRecord* sr) \
44  { \
45  if(sr->shw.nshwlid < 1) return -1000.; \
46  return double(sr->shw.shwlid[0].VAR); \
47  });
48 
52 
53 const Var kShwStopX = SHWLIDVAR(stop.X());
54 const Var kShwStopY = SHWLIDVAR(stop.Y());
55 const Var kShwStopZ = SHWLIDVAR(stop.Z());
56 
57 const Var kShwCalE = SHWLIDVAR(calE);
58 
60 const Var kShwWidth = SHWLIDVAR(width);
62 const Var kShwGap = SHWLIDVAR(gap);
63 
64 const Var kHadCalE({"shw.nshwlid", "shw.shwlid.*","slc.calE"},
65  [](const caf::StandardRecord* sr)
66  {
67  if(sr->shw.nshwlid < 1) return sr->slc.calE;
68  return (sr->slc.calE -
69  sr->shw.shwlid[0].calE);
70  });
71 
72 #define ELECIDVAR(VAR) \
73  Var({"sel.elecid.nshwlid", "sel.elecid.shwlid.*"}, \
74  [](const caf::StandardRecord* sr) \
75  { \
76  if(sr->sel.elecid.nshwlid < 1) return -1000.; \
77  return double(sr->sel.elecid.shwlid[0].VAR); \
78  });
79 
80 const Var kCosTheta = ELECIDVAR(costheta);
81 const Var kPi0Mass = ELECIDVAR(pi0mass);
82 const Var kShwEFrac = ELECIDVAR(shwEFrac);
83 
84 const Var kMuLLL = ELECIDVAR(mulll);
85 const Var kMuLLT = ELECIDVAR(mullt);
86 const Var kELLL = ELECIDVAR(elll);
87 const Var kELLT = ELECIDVAR(ellt);
88 const Var kEGLLL = ELECIDVAR(eglll);
89 const Var kEGLLT = ELECIDVAR(egllt);
90 const Var kEPLLL = ELECIDVAR(eplll);
91 const Var kEPLLT = ELECIDVAR(epllt);
92 const Var kENLLL = ELECIDVAR(enlll);
93 const Var kENLLT = ELECIDVAR(enllt);
94 const Var kEMuLLL = ELECIDVAR(emulll);
95 const Var kEMuLLT = ELECIDVAR(emullt);
96 const Var kEPi0LLL = ELECIDVAR(epi0lll);
97 const Var kEPi0LLT = ELECIDVAR(epi0llt);
98 const Var kEPiLLL = ELECIDVAR(epilll);
99 const Var kEPiLLT = ELECIDVAR(epillt);
100 
101 
102 const Var kConfusion({"sel.lem.pid", "sel.elecid.ann"},
103  [](const caf::StandardRecord* sr)
104  {
105  if(sr->sel.lem.pid < .8 && sr->sel.elecid.ann < .95) return .5;
106  if(sr->sel.lem.pid < .8 && sr->sel.elecid.ann >= .95) return 1.5;
107  if(sr->sel.lem.pid >= .8 && sr->sel.elecid.ann < .95) return 2.5;
108  if(sr->sel.lem.pid >= .8 && sr->sel.elecid.ann >= .95) return 3.5;
109 
110  });
111 
113 {
114  SpectrumLoader loaderData("defname: prod_caf_S15-05-22a_nd_numi_goodruns with stride 100");
115  SpectrumLoader loaderDataStagger("defname: prod_caf_S15-05-22a_nd_numi_stagger_goodruns with stride 100");
116  SpectrumLoader loaderMC("defname: prod_caf_S15-05-22_nd_genie_fhc_nonswap with stride 100");
117 
118  loaderData.SetSpillCut(kStandardSpillCuts);
119  loaderDataStagger.SetSpillCut(kStandardSpillCuts);
120  loaderMC.SetSpillCut(kStandardDQCuts);
121 
122  const Var kPIDExp = SIMPLEVAR(sel.lem.pidexp);
123  const Var kMeany = SIMPLEVAR(sel.lem.meanyexp);
124  const Var kQFrac = SIMPLEVAR(sel.lem.meanqfracexp);
125  const Var kEDiff = SIMPLEVAR(sel.lem.energydiffexp);
126  const Var kEnrich = SIMPLEVAR(sel.lem.enrichfracexp);
127  const Var kNElastic = SIMPLEVAR(vtx.nelastic);
128  const Var kNShwLID = SIMPLEVAR(shw.nshwlid);
129  const Var kLongestProng = SIMPLEVAR(sel.rvp.longtr);
130 
131  // const Var kPtp = SIMPLEVAR(sel.nuecosrej.photptp);
132  const Var kPtp({"sel.nuecosrej.photptp"},
133  [](const caf::StandardRecord* sr)
134  {
135  // This happened in the MC
136  if(isnan(sr->sel.nuecosrej.photptp) ||
137  isinf(sr->sel.nuecosrej.photptp)) return -1.f;
138  return sr->sel.nuecosrej.photptp;
139  });
140 
141  const Var kHitsPerPlane = SIMPLEVAR(sel.nuecosrej.hitsperplane);
142 
143  const Var kElecId = SIMPLEVAR(sel.elecid.ann);
144 
145 
146  struct HistDef
147  {
149  HistAxis axis;
150  };
151 
152  const int kNumVars = 52;
153  const HistDef defs[kNumVars] = {
154  {"lem", {"LEM", Binning::Simple(100, 0, 1), kLEM}},
155  {"lem_zoom", {"LEM", Binning::Simple(20, .8, 1), kLEM}},
156  {"nhits", {"Number of hits", Binning::Simple(50, 0, 250), kNHit}},
157  {"lem_pidexp", {"f_{sig}", Binning::Simple(100, 0, 1), kPIDExp}},
158  {"lem_meany", {"<y>", Binning::Simple(80, 0, .8), kMeany}},
159  {"lem_qfrac", {"<f_{Q}>", Binning::Simple(80, 0, .8), kQFrac}},
160  {"lem_ediff", {"D", Binning::Simple(100, -.02, +.02), kEDiff}},
161  {"lem_enrich", {"f_{enr}", Binning::Simple(100, 0, 1), kEnrich}},
162  {"calE", {"E_{cal}", Binning::Simple(100, 0, 5), kCaloE}},
163  {"nvtx", {"Number of vertices", Binning::Simple(3, 0, 3), kNElastic}},
164  {"vtxx", {"Vertex X", Binning::Simple(50, -200, +200), kVtxX}},
165  {"vtxy", {"Vertex Y", Binning::Simple(50, -200, +200), kVtxY}},
166  {"vtxz", {"Vertex Z", Binning::Simple(40, 0, +800), kVtxZ}},
167 
168  {"nshw", {"Number of showers", Binning::Simple(10, 0, 10), kNShwLID}},
169  {"shwstartx", {"Shower X Start", Binning::Simple(50, -200, +200), kShwStartX}},
170  {"shwstarty", {"Shower Y Start", Binning::Simple(50, -200, +200), kShwStartY}},
171  {"shwstartz", {"Shower Z Start", Binning::Simple(50, 0, +1000), kShwStartZ}},
172  {"shwstopx", {"Shower X Stop", Binning::Simple(50, -200, +200), kShwStopX}},
173  {"shwstopy", {"Shower Y Stop", Binning::Simple(50, -200, +200), kShwStopY}},
174  {"shwstopz", {"Shower Z Stop", Binning::Simple(70, 0, +1400), kShwStopZ}},
175  {"longprong", {"Longest prong (cm)", Binning::Simple(60, 0, 600), kLongestProng}},
176  {"e_per_hit", {"Calorimetric energy / hit (GeV)", Binning::Simple(40, 0, .04), kEPerHit}},
177  {"shw_nhit", {"Number of hits in primary shower", Binning::Simple(50, 0, 300), kShwNHit}},
178  {"shw_len", {"Length of primary shower (cm)", Binning::Simple(50, 0, 500), kShwLen}},
179  {"shw_width", {"Width of primary shower (cm)", Binning::Simple(20, 0, 20), kShwWidth}},
180  {"shw_gap", {"Primary shower gap to vertex (cm)", Binning::Simple(25, 0, 200), kShwGap}},
181  {"shw_efrac", {"Fraction of slice energy in shower", Binning::Simple(52, -.02, 1.02), kShwEFrac}},
182 
183  {"shw_e", {"Primary shower calorimetric energy (GeV)", Binning::Simple(30, 0, 3), kShwCalE}},
184  {"had_e", {"(Slice - Shower) calorimetric energy (GeV)", Binning::Simple(30, 0, 3), kHadCalE}},
185 
186  {"ptp", {"p_{T}/p", Binning::Simple(52, -.02, 1.02), kPtp}},
187  {"costheta", {"cos#theta", Binning::Simple(25, 0, +1), kCosTheta}},
188  {"hits_per_plane", {"Hits per plane", Binning::Simple(10, 0, 10), kHitsPerPlane}},
189  {"pi0mass", {"#pi^{0} mass (GeV)", Binning::Simple(50, 0, .5), kPi0Mass}},
190 
191  {"lid", {"LID", Binning::Simple(120, -0.1, 1.1), kElecId}},
192  {"lid_zoom", {"LID", Binning::Simple(20, 0.95, 1.05), kElecId}},
193 
194  {"mulll", {"#mu LLL", Binning::Simple(50, -1.2, 3), kMuLLL}},
195  {"elll", {"e LLL", Binning::Simple(50, -1.2, 3), kELLL}},
196  {"emulll", {"e/#mu LLL", Binning::Simple(50, -1.2, 3), kEMuLLL}},
197  {"eglll", {"e/#gamma LLL", Binning::Simple(50, -1.2, 1.2), kEGLLL}},
198  {"eplll", {"e/p LLL", Binning::Simple(50, -1.2, 1.2), kEPLLL}},
199  {"enlll", {"e/n LLL", Binning::Simple(50, -1.2, 1.2), kENLLL}},
200  {"epilll", {"e/#pi LLL", Binning::Simple(50, -1.2, 1.2), kEPiLLL}},
201  {"epi0lll", {"e/#pi^{0} LLL", Binning::Simple(50, -1.2, 1.2), kEPi0LLL}},
202 
203  {"mullt", {"#mu LLT", Binning::Simple(80, 0, 8), kMuLLT}},
204  {"ellt", {"e LLT", Binning::Simple(80, 0, 8), kELLT}},
205  {"emullt", {"e/#mu LLT", Binning::Simple(30, -1, +2), kEMuLLT}},
206  {"egllt", {"e/#gamma LLT", Binning::Simple(30, -1, +2), kEGLLT}},
207  {"epllt", {"e/p LLT", Binning::Simple(30, -1, +2), kEPLLT}},
208  {"enllt", {"e/n LLT", Binning::Simple(30, -1, +2), kENLLT}},
209  {"epillt", {"e/#pi LLT", Binning::Simple(30, -1, +2), kEPiLLT}},
210  {"epi0llt", {"e/#pi^{0} LLT", Binning::Simple(30, -1, +2), kEPi0LLT}},
211  {"confuse", {"Neither, LID, LEM, Both", Binning::Simple(4, 0, 4), kConfusion}},
212  };
213 
214  const int kNumSels = 4;
215  const Cut sels[kNumSels] = {kNoCut, kNueNDFirstAnaPresel, kNueNDFirstAnaPreselPlusLEM, kNueNDFirstAnaPreselPlusLID};
216  const std::string selNames[kNumSels] = {"nocut", "presel", "lem", "lid"};
217 
218  Spectrum* spects[kNumSels][kNumVars];
219  Spectrum* spectsStagger[kNumSels][kNumVars];
220  IPrediction* preds[kNumSels][kNumVars];
221 
222  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
223  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
224  const HistAxis& axis = defs[varIdx].axis;
225  spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx]);
226  spectsStagger[selIdx][varIdx] = new Spectrum(loaderDataStagger, axis, sels[selIdx]);
227  preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
228  axis.label, axis.bins, axis.var,
229  sels[selIdx]);
230  }
231  }
232 
233  loaderMC.Go();
234  loaderData.Go();
235  loaderDataStagger.Go();
236 
237  TFile* fout = new TFile("fout.root", "RECREATE");
238  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
239  TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
240  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
241  const char* name = defs[varIdx].name.c_str();
242  spects[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_%s", name)));
243  spectsStagger[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("spect_%s_stagger", name)));
244  preds[selIdx][varIdx]->SaveTo(d->mkdir(TString::Format("pred_%s", name)));
245  }
246  }
247 }
248 
249 #endif
const XML_Char int len
Definition: expat.h:262
const Var kShwNHit
const XML_Char * name
Definition: expat.h:151
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
const int kNumVars
Definition: vars.h:14
const Var kLEM
PID
Definition: Vars.cxx:26
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kEPi0LLL
#define SHWLIDVAR(VAR)
const Var kShwStopX
const Var kShwCalE
const Var kEGLLL
float pid
The default PID value - normally use this one.
Definition: SRLem.h:24
const Var kEGLLT
const Var kENLLT
std::string name
Definition: NuePlotLists.h:12
const Var kShwStartY
nhit
Definition: demo1.py:25
virtual void SaveTo(TDirectory *dir, const std::string &name) const
void SetSpillCut(const SpillCut &cut)
const Var kEPerHit([](const caf::SRProxy *sr){if(sr->slc.nhit >0) return 1000.0 *(sr->slc.calE/sr->slc.nhit);else return-5.;})
Definition: NueVarsExtra.h:14
#define ELECIDVAR(VAR)
const Var kVtxX
const Var kShwEFrac
const Var kHadCalE([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return float(sr->slc.calE);if(sr->vtx.elastic.fuzzyk.nshwlid< 1) return float(sr->slc.calE);return((sr->slc.calE- sr->vtx.elastic.fuzzyk.png[0].shwlid.calE));})
Definition: NueVarsExtra.h:40
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const Var kMeany
Float_t Y
Definition: plot.C:38
const Var kShwWidth
const Var kEnrich
Float_t Z
Definition: plot.C:38
const Var kENLLL
const Var kHitsPerPlane
const Var kEMuLLT
const Var kShwStopY
const Var kMuLLT
const Var kShwStartX
const Var kLongestProng([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return 0.f;if(sr->vtx.elastic.fuzzyk.npng==0) return 0.f;auto idx=sr->vtx.elastic.fuzzyk.longestidx;return float(sr->vtx.elastic.fuzzyk.png[idx].len);})
Definition: Vars.h:89
const Cut sels[kNumSels]
Definition: vars.h:44
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
const Var kEPiLLL
const Var kConfusion([](const caf::SRProxy *sr){if(sr->sel.lem.pid< .47 &&sr->sel.lid.ann< .63) return.5;if(sr->sel.lem.pid< .47 &&sr->sel.lid.ann >=.63) return 1.5;if(sr->sel.lem.pid >=.47 &&sr->sel.lid.ann< .63) return 2.5;if(sr->sel.lem.pid >=.47 &&sr->sel.lid.ann >=.63) return 3.5;return-5.;})
Definition: NueVarsExtra.h:95
const Var kNHit
Definition: Vars.cxx:71
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
Definition: SRIDBranch.h:48
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
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
const Var kEMuLLL
const Var kEPiLLT
const Var kEPi0LLT
Float_t d
Definition: plot.C:236
const Var kShwGap
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
const Var kMuLLL
const Var kShwStopZ
const HistDef defs[kNumVars]
Definition: vars.h:15
const Var kELLT
const Var kEDiff
const int kNumSels
Definition: vars.h:43
std::vector< float > Spectrum
Definition: Constants.h:527
unsigned int nhit
number of hits
Definition: SRSlice.h:22
const Var kELLL
HistAxis axis
Definition: NuePlotLists.h:13
const Var kPIDExp
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:94
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Var kVtxY
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kPi0Mass
SRIDBranch sel
Selector (PID) branch.
const Var kVtxZ
const Var kEPLLT
SRSlice slc
Slice branch: nhit, extents, time, etc.
void caf_nue_data_mc()
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
SRLem lem
Output from Library Event Matching (LEM)
Definition: SRIDBranch.h:43
const Var kEPLLL
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Var kNShwLID([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return 0u;return(unsigned int) sr->vtx.elastic.fuzzyk.nshwlid;})
Definition: NueVarsExtra.h:108
Prediction that just uses FD MC, with no extrapolation.
const Var kCosTheta
const std::string selNames[kNumSels]
Definition: vars.h:46
const Var kNElastic([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid){return 0.;}else{return 1.;}})
Definition: NueVarsExtra.h:106
const Var kShwStartZ
const Var kShwLen
Float_t X
Definition: plot.C:38
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:117
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const Var kQFrac
#define ELASTICVAR(VAR)