BPFEnergyVars.cxx
Go to the documentation of this file.
1 // Helper function guts for the BPF energy estimator.
2 
4 
5 #include "CAFAna/Core/Var.h"
9 
11 
12 
13 namespace ana
14 {
15 
16  //
17  // Define helper functions for computing energies:
18  //
19 
20 
21 
22  /// getMuonEnergy
23  /// \brief: function to compute the muon energy
24  ///
25  double getMuonEnergy(const TSpline3 & spline_calE, const TSpline3 & spline_BPFE, const caf::SRProxy* sr) {
26 
27  double muE = 0.0;
28 
29  // get best muon prong index
30  unsigned int muIdx = (unsigned int)kCVNMuonIdx(sr);
31  // std::cout << muIdx << std::endl;
32  // get the muon calE and BPFE
33  // As usual, here is a hardcoded assumption about there being only one vertex...
34  double muBPFE = -1.0;
35  double muCalE = sr->vtx.elastic.fuzzyk.png[muIdx].calE; // need to subtrack overlapE in the loop below
36 
37  // Assuming we've passed the BPF quality cuts, take the energy from the muon track.
38  if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid) {
39  muBPFE = sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.energy;
40  muCalE = muCalE - sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE; // subtract overlapping E from the muon track (to be put in the garbage bin)
41  }
42 
43  // adjust both E values with spline fits
44  double newCalE = spline_calE.Eval(muCalE);
45  double newBPFE = spline_BPFE.Eval(muBPFE);
46 
47  // Finally, combine newBPFE and newCalE.
48  // If no BPF muon track existed, fall back on calE only.
49  if(muBPFE < 0.0) muE = newCalE;
50  else muE = kMuonFitWeight * newBPFE + (1.0 - kMuonFitWeight) * newCalE;
51 
52  return muE;
53 
54  } // end getMuonEnergy
55 
56 
57 
58  /// getEMShowerEnergy
59  /// \brief: function to compute the EM shower energy
60  ///
61  double getEMShowerEnergy(const TSpline3 & spline_calE, const caf::SRProxy* sr) {
62 
63  double E = 0.0;
64  double calE = 0.0;
65 
66  // get best muon prong index
67  unsigned int muIdx = (unsigned int)kCVNMuonIdx(sr);
68 
69  // loop over all prongs, skip the muon prong, take all that pass EM criteria
70  // Again with the assumption about only having one vertex
71  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
72  if(p == muIdx) continue; // skip the muon
73  if(EMScore(p,sr) > 0.5) calE += sr->vtx.elastic.fuzzyk.png[p].weightedCalE;
74  } // end loop over 3D prongs
75 
76  // adjust calE with spline fits
77  E = spline_calE.Eval(calE);
78 
79  return E;
80 
81  } // end getEMShowerEnergy
82 
83 
84 
85  /// getHadronEnergy
86  /// \brief: function to compute the hadron energy
87  ///
88  double getHadronEnergy(const TSpline3 & spline_calE, const caf::SRProxy* sr) {
89 
90  double E = 0.0;
91  double calE = 0.0;
92 
93  // get best muon prong index
94  unsigned int muIdx = (unsigned int)kCVNMuonIdx(sr);
95 
96  // loop over all prongs, skip the muon prong, take all that pass general hadron criteria
97  // Again with the assumption about only having one vertex
98  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png.size(); ++p) {
99  if(p == muIdx) continue; // skip the muon
100  if(HadScoreWithMuon(p,sr) >= 0.5) calE += sr->vtx.elastic.fuzzyk.png[p].weightedCalE;
101  } // end loop over 3D prongs
102 
103  // adjust calE with spline fits
104  E = spline_calE.Eval(calE);
105 
106  return E;
107 
108  } // end getHadronEnergy
109 
110 
111 
112 
113 
114  /// takeOutTrash
115  /// \brief: function to compute v3 of the garbage bin.
116  /// this takes out muons, EM showers, and general hadrons, everything else in the garbage bin
117  ///
118  double takeOutTrash(const TSpline3 & spline_calE, const caf::SRProxy* sr) {
119 
120  double energy = 0.0;
121  double calE = 0.0;
122 
123  unsigned int muIdx = (unsigned int)kCVNMuonIdx(sr);
124 
125  // loop over all 2D prongs, sum weighted calE
126  for(unsigned int p = 0; p < sr->vtx.elastic.fuzzyk.png2d.size(); ++p) {
127  calE += sr->vtx.elastic.fuzzyk.png2d[p].weightedCalE;
128  } // end loop over 2D prongs
129 
130  // add in the orphaned hit calE
131  calE += sr->vtx.elastic.fuzzyk.orphCalE;
132 
133  // add in the overlapping energy from the muon track
134  double overlapE = 0.0;
135 
136  // Assuming we've passed the quality cuts, take the value from the muon track.
137  if(sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.IsValid)
138  overlapE = sr->vtx.elastic.fuzzyk.png[muIdx].bpf.muon.overlapE;
139  calE += overlapE;
140 
141 
142  // convert calE with the spline
143  energy = spline_calE.Eval(calE);
144 
145  return energy;
146 
147  } // end takeOutTrash
148 
149 
150 
151  const Var kBPFMuonE([](const caf::SRProxy *sr) {
152 
153  double muE = 0.0;
154 
155  int run = sr->hdr.run;
156  caf::Det_t det = sr->hdr.det;
157 
158  // split by detector
159  if(det == caf::kNEARDET) {
160  unsigned int period = PeriodFromRunND(run);
161  muE = getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);
162  }
163 
164  if(det == caf::kFARDET) {
165  unsigned int period = PeriodFromRunFD(run);
166  muE = getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);
167  }
168 
169  return muE;
170 
171  }); // end kBPFMuonE
172 
173 
174  const Var kBPFEME([](const caf::SRProxy *sr){
175 
176  double EME = 0.0;
177 
178  int run = sr->hdr.run;
179  caf::Det_t det = sr->hdr.det;
180 
181  // split by detector
182  if(det == caf::kNEARDET) {
183  unsigned int period = PeriodFromRunND(run);
184  EME = getEMShowerEnergy(allNDEMFits[period], sr);
185  }
186 
187  if(det == caf::kFARDET) {
188  unsigned int period = PeriodFromRunFD(run);
189  EME = getEMShowerEnergy(allFDEMFits[period], sr);
190  }
191 
192  return EME;
193 
194  });
195 
196 
197  const Var kBPFHadE([](const caf::SRProxy *sr){
198 
199  double hadE = 0.0;
200 
201  int run = sr->hdr.run;
202  caf::Det_t det = sr->hdr.det;
203 
204  // split by detector
205  if(det == caf::kNEARDET) {
206  unsigned int period = PeriodFromRunND(run);
207  hadE = getHadronEnergy(allNDHadFits[period], sr);
208  }
209 
210  if(det == caf::kFARDET) {
211  unsigned int period = PeriodFromRunFD(run);
212  hadE = getHadronEnergy(allFDHadFits[period], sr);
213  }
214 
215  return hadE;
216 
217  });
218 
219  const Var kBPFGarbageE([](const caf::SRProxy *sr){
220  double garE = 0.0;
221 
222  int run = sr->hdr.run;
223  caf::Det_t det = sr->hdr.det;
224 
225  // split by detector
226  if(det == caf::kNEARDET) {
227  unsigned int period = PeriodFromRunND(run);
228  garE = takeOutTrash(allNDGarbageFits[period], sr);
229  }
230 
231  if(det == caf::kFARDET) {
232  unsigned int period = PeriodFromRunFD(run);
233  garE = takeOutTrash(allFDGarbageFits[period], sr);
234  }
235 
236  return garE;
237 
238  });
239 
240 
241  const Var kBPFSliceE([](const caf::SRProxy *sr) {
242 
243  double energy = 0.0;
244 
245  int run = sr->hdr.run;
246  caf::Det_t det = sr->hdr.det;
247 
248  // split by detector
249  if(det == caf::kNEARDET) {
250 
251  // get the period
252  unsigned int period = PeriodFromRunND(run);
253 
254  // get the muon energy
255  energy += getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);
256 
257  // get the EM shower energy
258  energy += getEMShowerEnergy(allNDEMFits[period], sr);
259 
260  // get the hadron energy
261  energy += getHadronEnergy(allNDHadFits[period], sr);
262  // get the garbage bin energy
263  energy += takeOutTrash(allNDGarbageFits[period], sr);
264 
265  }
266 
267  if(det == caf::kFARDET) {
268 
269  // get the period
270  unsigned int period = PeriodFromRunFD(run);
271 
272  // get the muon energy
273  energy += getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);
274 
275  // get the EM shower energy
276  energy += getEMShowerEnergy(allFDEMFits[period], sr);
277 
278  // get the hadron energy
279  energy += getHadronEnergy(allFDHadFits[period], sr);
280 
281  // get the garbage bin energy
282  energy += takeOutTrash(allFDGarbageFits[period], sr);
283 
284  }
285 
286  return energy;
287 
288  }); // end kBPFSliceE
289 
290 
291 } //end namespace ana
const Var kBPFSliceE([](const caf::SRProxy *sr){double energy=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);energy+=getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergy(allNDEMFits[period], sr);energy+=getHadronEnergy(allNDHadFits[period], sr);energy+=takeOutTrash(allNDGarbageFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);energy+=getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);energy+=getEMShowerEnergy(allFDEMFits[period], sr);energy+=getHadronEnergy(allFDHadFits[period], sr);energy+=takeOutTrash(allFDGarbageFits[period], sr);}return energy;})
: Main Var for reconstructed neutrino energy.
Definition: BPFEnergyVars.h:58
Near Detector underground.
Definition: SREnums.h:10
Det_t
Which NOvA detector?
Definition: SREnums.h:7
caf::Proxy< caf::SRFuzzyK > fuzzyk
Definition: SRProxy.h:2059
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kBPFHadE([](const caf::SRProxy *sr){double hadE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);hadE=getHadronEnergy(allNDHadFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);hadE=getHadronEnergy(allFDHadFits[period], sr);}return hadE;})
Definition: BPFEnergyVars.h:51
caf::Proxy< std::vector< caf::SRProng > > png2d
Definition: SRProxy.h:2044
const std::vector< TSpline3 > allFDMuonFits_BPFE
double EMScore(unsigned int ProngIdx, const caf::SRProxy *sr)
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
const char * p
Definition: xmltok.h:285
const std::vector< TSpline3 > allNDHadFits
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
double HadScoreWithMuon(unsigned int ProngIdx, const caf::SRProxy *sr)
const Var kCVNMuonIdx([](const caf::SRProxy *sr){float longest_idx=-5.0;float longest_len=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.len > longest_len){longest_len=png.len;longest_idx=png_idx;}}} if(longest_len > 500.0) return longest_idx;float best_idx=-5.0;float best_score=-5.0;if(sr->vtx.elastic.IsValid){for(unsigned int png_idx=0;png_idx< sr->vtx.elastic.fuzzyk.png.size();png_idx++){auto &png=sr->vtx.elastic.fuzzyk.png[png_idx];if(png.cvnpart.muonid > best_score){best_score=png.cvnpart.muonid;best_idx=png_idx;}}}return best_idx;})
: Prong index of best muon prong by CVN score & length
Definition: CVNProngVars.h:19
const Var kBPFEME([](const caf::SRProxy *sr){double EME=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);EME=getEMShowerEnergy(allNDEMFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);EME=getEMShowerEnergy(allFDEMFits[period], sr);}return EME;})
Definition: BPFEnergyVars.h:50
caf::Proxy< unsigned int > run
Definition: SRProxy.h:248
unsigned int PeriodFromRunFD(int run)
: Helper function to compute the period given a run number for the FD.
Definition: RunPeriods.cxx:15
const std::vector< TSpline3 > allFDMuonFits_calE
caf::Proxy< caf::SRElastic > elastic
Definition: SRProxy.h:2118
caf::Proxy< std::vector< caf::SRFuzzyKProng > > png
Definition: SRProxy.h:2043
const std::vector< TSpline3 > allNDMuonFits_calE
const std::vector< TSpline3 > allNDEMFits
const Var kBPFMuonE([](const caf::SRProxy *sr){double muE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);muE=getMuonEnergy(allNDMuonFits_calE[period], allNDMuonFits_BPFE[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);muE=getMuonEnergy(allFDMuonFits_calE[period], allFDMuonFits_BPFE[period], sr);}return muE;})
: Returns the recoE of the CVN selected muon.
Definition: BPFEnergyVars.h:48
Float_t E
Definition: plot.C:20
double getEMShowerEnergy(const TSpline3 &spline_calE, const caf::SRProxy *sr)
: function to compute the EM shower energy
static const double kMuonFitWeight
double energy
Definition: plottest35.C:25
const std::vector< TSpline3 > allFDHadFits
caf::StandardRecord * sr
const std::vector< TSpline3 > allFDGarbageFits
double takeOutTrash(const TSpline3 &spline_calE, const caf::SRProxy *sr)
: function to compute v3 of the garbage bin. this takes out muons, EM showers, and general hadrons...
Definition: run.py:1
caf::Proxy< float > orphCalE
Definition: SRProxy.h:2042
double getMuonEnergy(const TSpline3 &spline_calE, const TSpline3 &spline_BPFE, const caf::SRProxy *sr)
: function to compute the muon energy
double getHadronEnergy(const TSpline3 &spline_calE, const caf::SRProxy *sr)
: function to compute the hadron energy
const std::vector< TSpline3 > allNDGarbageFits
const Var kBPFGarbageE([](const caf::SRProxy *sr){double garE=0.0;int run=sr->hdr.run;caf::Det_t det=sr->hdr.det;if(det==caf::kNEARDET){unsigned int period=PeriodFromRunND(run);garE=takeOutTrash(allNDGarbageFits[period], sr);}if(det==caf::kFARDET){unsigned int period=PeriodFromRunFD(run);garE=takeOutTrash(allFDGarbageFits[period], sr);}return garE;})
Definition: BPFEnergyVars.h:52
const std::vector< TSpline3 > allNDMuonFits_BPFE
caf::Proxy< caf::SRVertexBranch > vtx
Definition: SRProxy.h:2146
const std::vector< TSpline3 > allFDEMFits
unsigned int PeriodFromRunND(int run)
: Helper function to compute the period given a run number for the ND.
Definition: RunPeriods.cxx:34
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:232