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:40
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:16
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:506
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:759
unsigned int nhit
number of hits
Definition: SRSlice.h:22
const Var kELLL
HistAxis axis
Definition: NuePlotLists.h:13
const Var kPIDExp
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
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
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:49
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:107
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const Var kQFrac
#define ELASTICVAR(VAR)
enum BeamMode string