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"
12 #include "CAFAna/Vars/Vars.h"
13 using namespace ana;
14 
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;
216  const std::string selNames[kNumSels] = {"nocut", "presel", "lem", "lid"};
217 
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 Cut kNueNDFirstAnaPreselPlusLID
const Var kShwNHit
const XML_Char * name
Definition: expat.h:151
int isinf(const stan::math::var &a)
Definition: std_isinf.hpp:16
const Cut kNueNDFirstAnaPresel
const int kNumVars
Definition: vars.h:14
const Var kLEM
PID
Definition: Vars.cxx:24
Oscillation analysis framework, runs over CAF files outside of ART.
const Var kEPi0LLL
#define SHWLIDVAR(VAR)
virtual void SaveTo(TDirectory *dir) const
const Var kShwStopX
const Var kShwCalE
const Var kEGLLL
const Var kEGLLT
const Var kENLLT
const Var kShwStartY
HistAxis axis
Definition: NuePlotLists.h:13
nhit
Definition: demo1.py:25
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:37
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 kPtp([](const caf::SRProxy *sr){if(isnan(sr->sel.nuecosrej.photptp)|| isinf(sr->sel.nuecosrej.photptp)) return-1.f;return sr->sel.nuecosrej.photptp;})
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:84
const Cut sels[kNumSels]
Definition: vars.h:44
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:50
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:91
const Var kNHit
Definition: Vars.cxx:69
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:12
const Var kEMuLLL
const Var kEPiLLT
const Var kEPi0LLT
Float_t d
Definition: plot.C:236
const Var kShwGap
virtual void Go() override
Load all the registered spectra.
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
const Var kELLL
std::string name
Definition: NuePlotLists.h:12
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
const Var kVtxZ
const Var kEPLLT
const Cut kNueNDFirstAnaPreselPlusLEM
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:28
const Var kEPLLL
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void spects(osc::IOscCalculatorAdjustable *calc, IPrediction *predFHC, IPrediction *predRHC)
Definition: sensitivity.C:100
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:104
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:102
const Var kShwStartZ
const Var kShwLen
Float_t X
Definition: plot.C:38
const Var kElecId
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
sel
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.
static constexpr Double_t sr
Definition: Munits.h:164
const Var kQFrac
void SaveTo(TDirectory *dir) const
Definition: Spectrum.cxx:1029
#define ELASTICVAR(VAR)