NumuEAlg.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 // $Id: NumuEAlg.cxx,v 1.10 2012-11-15 20:46:24 lein Exp $
3 //
4 // Numu Energy Algorithms
5 //
6 // lein@physics.umn.edu
7 //
8 ////////////////////////////////////////////////////////////////////////
9 
10 // NOvA includes
11 #include "NumuEnergy/NumuEAlg.h"
12 #include "NumuEnergy/func/Numu2017Fits.h"
13 #include "NumuEnergy/func/Numu2018Fits.h"
14 #include "NumuEnergy/func/Numu2020Fits.h"
15 #include "ReMId/ReMId.h"
16 #include "RecoBase/RecoHit.h"
17 #include "Geometry/Geometry.h"
18 #include "GeometryObjects/CellGeo.h"
19 #include "NovaDAQConventions/DAQConventions.h"
20 #include "Utilities/func/EnvExpand.h"
21 #include "Utilities/func/MathUtil.h"
22 #include "NumuEnergy/ScatterEnergy.h"
23 #include "MCCheater/BackTracker.h"
24 
25 // Art includes
26 #include "messagefacility/MessageLogger/MessageLogger.h"
27 #include "SummaryData/SpillData.h"
28 
29 #include "cetlib/search_path.h"
30 
31 #include <cmath>
32 
33 namespace numue {
34 
35  //----------------------------------------------------------------------
36  NumuEAlg::NumuEAlg(fhicl::ParameterSet const& pset):
37  fPIDModuleLabel (pset.get< std::string>("PIDModuleLabel" )),
38  fPIDMethod (pset.get< bool >("PIDMethod" )),
39  fTMVAName (pset.get< std::string>("TMVAFileName", "UC_Eres_BDTG.weights.xml")),
40  fTMVANameSingle (pset.get< std::string>("TMVAFileSingle", "UC_MuonE_Single_BDTG.weights.xml")),
41  fTMVANameNonSingle (pset.get< std::string>("TMVAFileNonSingle", "UC_MuonE_NonSingle_BDTG.weights.xml")),
42  //fTMVApathUCMuonE (pset.get< std::string >("tmvapathUCMuonE")), //For local tests
43  fFD (false),
44  fND (false),
45  fMaxTrkLen (-1),
46  fNumiDir (-1,-1,-1),
47  fTranPlane (-1),
48  fTranZPos (-1)
49  {
50  art::ServiceHandle<geo::Geometry> geom;
51 
52  novadaq::cnv::DetId detID = geom->DetId();
53  if (detID == novadaq::cnv::kNEARDET) fND = true;
54  if (detID == novadaq::cnv::kFARDET) fFD = true;
55 
56  //We want the plane just BEFORE the muon catcher, hence the -1.
57  //We want the general center of the z plane - so I picked cell 32.
58  //This could be made more advanced / clever in the future if needed.
59  if (fND){
60  fTranPlane = (geom->FirstPlaneInMuonCatcher())-1;
61  double xyz[3];
62  geom->Plane(fTranPlane)->Cell(32)->GetCenter(xyz);
63  fTranZPos = xyz[2];
64  }
65 
66  double detHalfWidth = geom->DetHalfWidth();
67  double detHalfHeight = geom->DetHalfHeight();
68  double detLength = geom->DetLength();
69 
70  fMaxTrkLen = util::pythag(detHalfWidth*2, detHalfHeight*2, detLength);
71  fNumiDir = geom->NuMIBeamDirection();
72 
73  fReaderE = new TMVA::Reader;
74  fReaderE->AddVariable("slc.calE", &Evars[0]);
75  fReaderE->AddVariable("sel.cosrej.anglekal", &Evars[1]);
76  fReaderE->AddVariable("sel.cosrej.scatt", &Evars[2]);
77  fReaderE->AddVariable("sel.cosrej.eratio", &Evars[3]);
78  fReaderE->AddVariable("sel.cosrej.kalfwdcell", &Evars[4]);
79  fReaderE->AddVariable("trk.cosmic[0].nhit", &Evars[5]);
80  fReaderE->AddVariable("sel.cosrej.ntracks3d", &Evars[6]);
81 
82  std::string fileName;
83  cet::search_path sp("UCANA_LIB_PATH");
84  sp.find_file(fTMVAName, fileName);
85 
86  fReaderE->BookMVA("BDTG",fileName.c_str());
87 
88  //---
89  fReaderUCMuonSingle = new TMVA::Reader;
90  fReaderUCMuonSingle->AddVariable("anglekal", &TMVAvarsSingle[0]);
91  fReaderUCMuonSingle->AddVariable("tracklen", &TMVAvarsSingle[1]);
92  fReaderUCMuonSingle->AddVariable("trkEPerNHit", &TMVAvarsSingle[2]);
93  fReaderUCMuonSingle->AddVariable("scatt_D_tracklen", &TMVAvarsSingle[3]);
94  fReaderUCMuonSingle->AddVariable("hadclustcalE", &TMVAvarsSingle[4]);
95 
96  fReaderUCMuonSingle->AddSpectator("BDTA_700", &TMVAvarsSingle[5]);
97  fReaderUCMuonSingle->AddSpectator("anglekal", &TMVAvarsSingle[6]);
98  fReaderUCMuonSingle->AddSpectator("pngptp", &TMVAvarsSingle[7]);
99  fReaderUCMuonSingle->AddSpectator("slcT", &TMVAvarsSingle[8]);
100  fReaderUCMuonSingle->AddSpectator("osc", &TMVAvarsSingle[9]);
101  fReaderUCMuonSingle->AddSpectator("oscVal1", &TMVAvarsSingle[10]);
102  fReaderUCMuonSingle->AddSpectator("nhit", &TMVAvarsSingle[11]);
103  fReaderUCMuonSingle->AddSpectator("trueE", &TMVAvarsSingle[12]);
104  fReaderUCMuonSingle->AddSpectator("numuE", &TMVAvarsSingle[13]);
105  fReaderUCMuonSingle->AddSpectator("truemuonE", &TMVAvarsSingle[14]);
106  fReaderUCMuonSingle->AddSpectator("recomuonE", &TMVAvarsSingle[15]);
107  fReaderUCMuonSingle->AddSpectator("hadclustcalE", &TMVAvarsSingle[16]);
108  fReaderUCMuonSingle->AddSpectator("hadcalE", &TMVAvarsSingle[17]);
109  fReaderUCMuonSingle->AddSpectator("hadtrkE", &TMVAvarsSingle[18]);
110  fReaderUCMuonSingle->AddSpectator("ntracks3d", &TMVAvarsSingle[19]);
111  fReaderUCMuonSingle->AddSpectator("npngs3d", &TMVAvarsSingle[20]);
112  fReaderUCMuonSingle->AddSpectator("vtx", &TMVAvarsSingle[21]);
113  fReaderUCMuonSingle->AddSpectator("vty", &TMVAvarsSingle[22]);
114  fReaderUCMuonSingle->AddSpectator("vtz", &TMVAvarsSingle[23]);
115  fReaderUCMuonSingle->AddSpectator("endx", &TMVAvarsSingle[24]);
116  fReaderUCMuonSingle->AddSpectator("endy", &TMVAvarsSingle[25]);
117  fReaderUCMuonSingle->AddSpectator("endz", &TMVAvarsSingle[26]);
118  fReaderUCMuonSingle->AddSpectator("vetoangl", &TMVAvarsSingle[27]);
119  fReaderUCMuonSingle->AddSpectator("tracklen2nd", &TMVAvarsSingle[28]);
120  fReaderUCMuonSingle->AddSpectator("trkE2nd", &TMVAvarsSingle[29]);
121  fReaderUCMuonSingle->AddSpectator("angle2nd", &TMVAvarsSingle[30]);
122  fReaderUCMuonSingle->AddSpectator("run", &TMVAvarsSingle[31]);
123  fReaderUCMuonSingle->AddSpectator("subrun", &TMVAvarsSingle[32]);
124  fReaderUCMuonSingle->AddSpectator("evt", &TMVAvarsSingle[33]);
125  fReaderUCMuonSingle->AddSpectator("subevt", &TMVAvarsSingle[34]);
126 
127  // Only if path is an ENV
128  std::string fileNameSingle;
129  cet::search_path spSingle("UCANA_LIB_PATH");
130  spSingle.find_file(fTMVANameSingle, fileNameSingle);
131  fReaderUCMuonSingle->BookMVA("UC_BDTG_S",fileNameSingle.c_str());
132 
133  //std::string pidpathUCMuonE = util::EnvExpansion(fTMVApathUCMuonE);//TMP
134  //fReaderUCMuonSingle->BookMVA("UC_BDTG_S",pidpathUCMuonE+"UC_MuonE_Single_BDTG.weights.xml");//TMP
135 
136 
137  fReaderUCMuonNonSingle = new TMVA::Reader;
138  fReaderUCMuonNonSingle->AddVariable("anglekal", &TMVAvarsNonSingle[0]);
139  fReaderUCMuonNonSingle->AddVariable("tracklen", &TMVAvarsNonSingle[1]);
140  fReaderUCMuonNonSingle->AddVariable("trkEPerNHit", &TMVAvarsNonSingle[2]);
141  fReaderUCMuonNonSingle->AddVariable("scatt_D_tracklen", &TMVAvarsNonSingle[3]);
142  fReaderUCMuonNonSingle->AddVariable("hadclustcalE", &TMVAvarsNonSingle[4]);
143  fReaderUCMuonNonSingle->AddVariable("tracklen2nd", &TMVAvarsNonSingle[5]);
144  fReaderUCMuonNonSingle->AddVariable("trkE2nd", &TMVAvarsNonSingle[6]);
145  fReaderUCMuonNonSingle->AddVariable("angle2nd", &TMVAvarsNonSingle[7]);
146 
147  fReaderUCMuonNonSingle->AddSpectator("BDTA_700", &TMVAvarsNonSingle[8]);
148  fReaderUCMuonNonSingle->AddSpectator("anglekal", &TMVAvarsNonSingle[9]);
149  fReaderUCMuonNonSingle->AddSpectator("pngptp", &TMVAvarsNonSingle[10]);
150  fReaderUCMuonNonSingle->AddSpectator("slcT", &TMVAvarsNonSingle[11]);
151  fReaderUCMuonNonSingle->AddSpectator("osc", &TMVAvarsNonSingle[12]);
152  fReaderUCMuonNonSingle->AddSpectator("oscVal1", &TMVAvarsNonSingle[13]);
153  fReaderUCMuonNonSingle->AddSpectator("nhit", &TMVAvarsNonSingle[14]);
154  fReaderUCMuonNonSingle->AddSpectator("trueE", &TMVAvarsNonSingle[15]);
155  fReaderUCMuonNonSingle->AddSpectator("numuE", &TMVAvarsNonSingle[16]);
156  fReaderUCMuonNonSingle->AddSpectator("truemuonE", &TMVAvarsNonSingle[17]);
157  fReaderUCMuonNonSingle->AddSpectator("recomuonE", &TMVAvarsNonSingle[18]);
158  fReaderUCMuonNonSingle->AddSpectator("hadclustcalE", &TMVAvarsNonSingle[19]);
159  fReaderUCMuonNonSingle->AddSpectator("hadcalE", &TMVAvarsNonSingle[20]);
160  fReaderUCMuonNonSingle->AddSpectator("hadtrkE", &TMVAvarsNonSingle[21]);
161  fReaderUCMuonNonSingle->AddSpectator("ntracks3d", &TMVAvarsNonSingle[22]);
162  fReaderUCMuonNonSingle->AddSpectator("npngs3d", &TMVAvarsNonSingle[23]);
163  fReaderUCMuonNonSingle->AddSpectator("vtx", &TMVAvarsNonSingle[24]);
164  fReaderUCMuonNonSingle->AddSpectator("vty", &TMVAvarsNonSingle[25]);
165  fReaderUCMuonNonSingle->AddSpectator("vtz", &TMVAvarsNonSingle[26]);
166  fReaderUCMuonNonSingle->AddSpectator("endx", &TMVAvarsNonSingle[27]);
167  fReaderUCMuonNonSingle->AddSpectator("endy", &TMVAvarsNonSingle[28]);
168  fReaderUCMuonNonSingle->AddSpectator("endz", &TMVAvarsNonSingle[29]);
169  fReaderUCMuonNonSingle->AddSpectator("vetoangl", &TMVAvarsNonSingle[30]);
170  fReaderUCMuonNonSingle->AddSpectator("tracklen2nd", &TMVAvarsNonSingle[31]);
171  fReaderUCMuonNonSingle->AddSpectator("trkE2nd", &TMVAvarsNonSingle[32]);
172  fReaderUCMuonNonSingle->AddSpectator("angle2nd", &TMVAvarsNonSingle[33]);
173  fReaderUCMuonNonSingle->AddSpectator("run", &TMVAvarsNonSingle[34]);
174  fReaderUCMuonNonSingle->AddSpectator("subrun", &TMVAvarsNonSingle[35]);
175  fReaderUCMuonNonSingle->AddSpectator("evt", &TMVAvarsNonSingle[36]);
176  fReaderUCMuonNonSingle->AddSpectator("subevt", &TMVAvarsNonSingle[37]);
177 
178  // Only if path is an ENV
179  std::string fileNameNonSingle;
180  cet::search_path spNonSingle("UCANA_LIB_PATH");
181  spNonSingle.find_file(fTMVANameNonSingle, fileNameNonSingle);
182  fReaderUCMuonNonSingle->BookMVA("UC_BDTG_NS",fileNameNonSingle.c_str());
183 
184  //std::string pidpathUCMuonE = util::EnvExpansion(fTMVApathUCMuonE);//NO NEED, DECLARED ABOVE
185  //fReaderUCMuonNonSingle->BookMVA("UC_BDTG_NS",pidpathUCMuonE+"UC_MuonE_NonSingle_BDTG.weights.xml");//TMP
186  //---
187 
188  return;
189  }
190 
191  //----------------------------------------------------------------------
192  NumuEAlg::~NumuEAlg()
193  {
194  if (fReaderE) delete fReaderE;
195  }
196 
197  //----------------------------------------------------------------------
198  NumuE NumuEAlg::FDEnergy(std::vector< art::Ptr<rb::Track> > const sliceTracks,
199  art::PtrVector< rb::CellHit > sliceHits,
200  art::Ptr< rb::Cluster> sliceCluster,
201  art::Event const& e,
202  murem::TrackCleanUpAlg* trkCleanUpAlg,
203  bool isRHC)
204  {
205  /* This returns multiple neutrino energy estimates, some of which are better than others.
206  They are: CalCC, TrkQEE, TrkNonQEE, and TrkCCE. The default E value is set to be indentically
207  equal to TrkCCE and comes with the same caveats. At the same time, we are also
208  setting HadCalE and HadTrkE. If no 3D track is found, no energy estimates are made.
209  When using the PIDMethod (which is default), the muon track is defined as the 3D track
210  with the highest PID. Otherwise, we simply take the longest 3D track as the muon.
211 
212  CalCC is a very simple, first pass at energy estimation for the numu cc group.
213  ONE SHOULD ALMOST NEVER USE CalCC! THE OTHER ENERGIES ARE USUALLY BETTER CHOICES!
214  The weights and methods used were optimized for a contained population that
215  doesn't distinguish between nonqe and qe events. No background existed in the
216  optimized sample. It used a linear fit to find the two weights and it preformed
217  the fit simultaneously. The optimization was done only for neutrinos with true
218  energy between 1 and 3 GeV and basic activity cuts. The muon energy is calculated
219  by summing the reco GeV for every hit on the track. The vertex energy is simply
220  a summed GeV of every hit NOT on the muon track.
221 
222  The form of the neutrino energy for CalCC is:
223  Reco Nu Energy = Weight_muon * Summed GeV of Muon Track
224  + Weight_vertex * Summed GeV of all non-Muon Track hits
225 
226  The CalCC method is NOT our best energy estimator and is mostly an accident of history.
227 
228  TrkQE, TrkNonQE and TrkCC are very similar. The weights and methods used were optimized
229  for a contained population that distinguishes by truth between QE and NonQE events.
230  The muon energy is found by a spline fit
231  to muon track length and muon true energy. This function is MuonEFromTrackLength and more info
232  can be found at the top of that function. Then the hadronic energy is defined as the sum of HadCalE
233  and HadTrkE. HadCalE is the same vertex energy as used in CalCC - it is just a sum of the GeV
234  of all hits not on the muon track. HadTrkE is the hadronic energy determined to be on
235  the muon track. This is found by a function in MuonRemove. Summed together, this is the visible hadronic
236  energy. Then this energy is transformed into our reco hadronic energy
237  by a spline fit. This fit is different for the NonQE, QE and combined populations. It is a fit between
238  the summed hadronic energy and (True Neutrino Energy - Reco Muon Energy). More information can
239  be found in the functions HadEQE, HadENonQE, and HadECC. Finally, the neutrino energy is just as
240  sum of the Muon Energy and the adjusted Hadronic Energy.
241 
242  Currently, all the energy methods are only really meant for contained events; however, this function
243  and the module that calls it doesn't not attempt to determine containment. Please
244  use with discretion!
245  */
246 
247  NumuE sliceEnergy;
248 
249  int bestTrack = -1;
250  double trackLen = 0.0;
251 
252  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
253  else bestTrack = LongestTrack(sliceTracks, trackLen);
254 
255  if(bestTrack == -1) return sliceEnergy;
256  if(bestTrack == 999) return sliceEnergy;
257 
258  if (fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
259 
260  //Gets GeV of hits on track
261  double recoTrkGeV = sliceTracks[bestTrack]->TotalGeV();
262 
263  //Calculates GeV of muon from track length
264  fRun = e.run();
265  bool ismc = checkIsMC(e);
266 
267 
268  double muonGeVFromLen = FDMuonEFromTrackLength(trackLen,
269  fRun,
270  ismc,
271  isRHC);
272 
273  art::PtrVector< rb::CellHit > trackHits;
274  trackHits = sliceTracks[bestTrack]->AllCells();
275 
276  //Getting all the hits in the slice that are not on our best track
277  rb::Cluster vertexHits = MakeVertexCluster(sliceHits,trackHits);
278 
279  //Gets GeV of hadronic hits NOT on track
280  double vertexEnergy = 0.0;
281  if (vertexHits.NCell()>0){
282  vertexEnergy = vertexHits.TotalGeV();
283  }
284 
285  //Gets GeV of hadronic hits on track
286  double extraHadE = trkCleanUpAlg->ExtraEOnTrackInGeV(*sliceTracks[bestTrack],*sliceCluster);
287 
288  //Get adjusted GeV of total hadronic energy assuming QE/NonQE
289  double totHadE = vertexEnergy + extraHadE;
290  double hadQEGeV = HadEQE(totHadE);
291  double hadNonQEGeV = HadENonQE(totHadE);
292  double hadCCGeV = HadECC(totHadE, fRun, ismc, isRHC);
293 
294  //***************Fill values in our slice energy object
295 
296  double calCCMuonWeight = 1.81134; //Using Summed GeV method and Susan's fit from MN talk on 9/7/2012
297  double calCCVertWeight = 2.17923;
298  double calcc = calCCMuonWeight*recoTrkGeV + calCCVertWeight*vertexEnergy;
299  sliceEnergy.SetCalCCE(calcc);
300 
301  sliceEnergy.SetHadCalE(vertexEnergy);
302  sliceEnergy.SetHadTrkE(extraHadE);
303 
304  double qeE = muonGeVFromLen + hadQEGeV;
305  sliceEnergy.SetTrkQEE(qeE);
306 
307  double nonqeE = muonGeVFromLen + hadNonQEGeV;
308  sliceEnergy.SetTrkNonQEE(nonqeE);
309 
310  double ccE = muonGeVFromLen + hadCCGeV;
311  sliceEnergy.SetTrkCCE(ccE);
312  sliceEnergy.SetRecoMuonE(muonGeVFromLen);
313  sliceEnergy.SetRecoTrkCCHadE(hadCCGeV);
314 
315  sliceEnergy.SetE(ccE); //For now, E is just set to be identical to Trk CC E
316 
317  sliceEnergy.SetHadCluster(vertexHits);
318 
319  return sliceEnergy;
320 
321  }//End of FDEnergy function
322 
323  //----------------------------------------------------------------------
324  rb::Energy NumuEAlg::QEFormulaEnergy(std::vector< art::Ptr<rb::Track> > const sliceTracks,
325  art::PtrVector< rb::CellHit > sliceHits,
326  double &error,
327  art::Event const& e,
328  bool isRHC)
329  {
330  /* This estimates the neutrino energy using the quasi-elastic formula.
331  If no 3D track is found, it doesn't return an energy.
332  It uses the 3D Kalman track with the best ReMId PID as the "muon", or the
333  longest track (as set by fcl, best PID recommended).
334 
335  Nuclear binding energy is set to 25 MeV. This is consistant with Minerba's
336  cross-section analysis, Susan's work from Feb. 2012, and what GENIE uses
337  in the simulation.
338 
339  The NuMI beam direction is determined from a geometry function for this purpose.
340 
341  Future: one could consider returning a numuE object instead of an Energy as well
342  as a separate error value.
343  */
344 
345  int bestTrack = -1;
346  double trackLen = 0.0;
347 
348  MuonType trackType;
349 
350  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
351  else bestTrack = LongestTrack(sliceTracks, trackLen);
352 
353  if(bestTrack == -1) return 0;
354  if(bestTrack == 999) return 0;
355 
356  if (fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
357 
358  art::PtrVector< rb::CellHit > trackHits;
359 
360  TVector3 mudir = sliceTracks[bestTrack]->Dir(); // get direction of muon
361  if(mudir.Mag()!=1) mudir = mudir.Unit(); // double check is proper unit vector
362 
363  double costheta = mudir * fNumiDir; // get cosine of angle
364 
365  fRun = e.run();
366  bool ismc = checkIsMC(e);
367 
368  double muE = 0.0; // muon energy (GeV)
369  if (fFD) {
370  muE = FDMuonEFromTrackLength(trackLen, fRun, ismc, isRHC);
371  trackType = kFDMuon;
372  }
373  if (fND) {
374  muE = NDMuonEFromTrackLength(sliceTracks[bestTrack], isRHC, &trackType);
375  }
376 
377  double Mp = .938272; // proton mass (GeV)
378  double Mu = .105658; // muon mass (GeV)
379  double Mb = 0.025; // nuclear binding energy (GeV)
380  double Mn = 0.939565; // neutron mass
381  double Mnp = Mn - Mb; // neutron mass - carbon nuclear binding energy (GeV)
382  double pmu_2 = muE*muE - (Mu*Mu); // Muon momentum squared
383 
384  double recoNu = (2*Mnp*muE-(Mnp*Mnp+Mu*Mu-(Mp*Mp))) / (2*(Mnp-muE+sqrt(pmu_2)*costheta)); // QE Formula
385 
386  // calculate the error on this energy estimation
387  double errMb = 0.005;// Uncertainty in nuclear binding energy (GeV)
388 
389  double resE = 0.0;//These resolutions are the RMS of (Reco-True)/True Muon Energy plots - should be updated
390  if (trackType == kFDMuon) {resE = 0.043;}
391  else if (trackType == kNDActiveMuon) {resE = 0.068;}
392  else if (trackType == kNDCatcherMuon){resE = 0.102;}
393  else if (trackType == kNDActCatMuon) {resE = 0.057;}
394  else if (trackType == kNDPiddleMuon) {resE = 0.084;}
395  else {error = -1;}
396 
397  double errE = resE*muE;
398 
399  double rescostheta = 0.05; // resolution of track direction
400 
401  // variance of nuclear binding energy, track length, direction
402  double varMb = errMb*errMb;
403  double varE = errE*errE;
404  double varcostheta = (rescostheta*costheta)*(rescostheta*costheta);
405 
406  double x = 2*Mnp*muE - (Mnp*Mnp + Mu*Mu - (Mp*Mp));
407  double varx = 4*(Mn*Mn+Mb*Mb)*varE + 4*(Mn*Mn+muE*muE +0.5*Mb*Mb)*varMb;
408  double vary = 4*varMb + 4*varE + 4*((muE*muE*costheta*costheta*varE)/(2*pmu_2) + pmu_2*varcostheta);
409 
410  double sigmaE_2 = recoNu*recoNu*((varx+recoNu*recoNu*vary)/(x*x));
411  error = sqrt(sigmaE_2);
412 
413  if (std::isnan(error)) error = -1;
414  if (std::isnan(recoNu)) recoNu = -1;
415 
416  rb::Energy sliceEnergy(recoNu);
417  return sliceEnergy;
418 
419  }//End of QEFormulaEnergy function
420 
421 
422  //----------------------------------------------------------------------
423  NumuE NumuEAlg::NDEnergy(std::vector< art::Ptr<rb::Track> > const sliceTracks,
424  art::PtrVector< rb::CellHit > sliceHits,
425  art::Ptr< rb::Cluster> sliceCluster,
426  art::Event const& e,
427  murem::TrackCleanUpAlg* trkCleanUpAlg,
428  bool isRHC)
429  {
430  /* This returns multiple neutrino energy estimates, some of which are better than others.
431  They are: CalCC, TrkQEE, TrkNonQEE and TrkCCE. The default E value is set to be indentically
432  equal to TrkCCE and comes with the same caveats. At the same time, we are also
433  setting HadCalE and ND specific information about the distribution of energy relative
434  to the muon catcher. If no 3D track is found, no energy estimates are made.
435  When using the PIDMethod (which is default), the muon track is defined as the 3D track
436  with the highest PID. Otherwise, we simply take the longest 3D track as the muon.
437 
438  CalCC is a very simple, first pass at energy estimation for the numu cc group.
439  ONE SHOULD ALMOST NEVER USE CalCC! THE OTHER ENERGIES ARE USUALLY BETTER CHOICES!
440  The weights and methods used were optimized for a contained population that
441  doesn't distinguish between nonqe and qe events. No background existed in the
442  optimized sample. It used a linear fit to find the weights and it performed
443  the fit simultaneously. It was optimized using Kalaman tracks.
444  The optimization was done only for neutrinos with true energy between 1 and 3
445  GeV and basic activity cuts.
446 
447  The form of the neutrino energy for CalCC is:
448  Reco Nu Energy = Weight_muon_active * Summed GeV of Muon Track in active region
449  + Weight_muon_tran * Summed GeV of Muon Track in transition plane
450  + Weight_muon_steel * Summed GeV of Muon Track in steel section
451  + Weight_vertex * Summed GeV of all non-Muon Track hits
452  + Weight_offset
453 
454  The weights are listed directly following this comment block.
455 
456  The CalCC method is NOT our best energy estimator and exists as an accident of history.
457 
458  TrkQE, TrkNonQE and TrkCCE are very similar. The weights and methods used were optimized
459  for a contained population that distinguishes by truth between QE and NonQE events.
460  The muon energy is found by a spline fit
461  to muon track length and muon true energy. This relies heavily on MuonEFromTrackLength but first
462  takes into account the muon track length in the active region vs. the catcher region.
463  Then the hadronic energy is defined as the sum of HadCalE and HadTrkE. HadCalE is the same vertex
464  energy as used in CalCC - it is just a sum of the GeV of all hits not on the muon track. HadTrkE
465  is the hadronic energy determined to be on the muon track. This is found by a function in
466  MuonRemove. Summed together, this is the visible hadronic energy. This energy
467  is broken down into three regions - the active region, the transition plane, and the muon
468  catcher. For now, all fits have been done on a population that has primarily hadronic energy in
469  the active region. However, in order to return an energy for every case, we actually use
470  the total hadronic region without regards to where it was located. The slices which have
471  hadronic energy in regions other than fully active will have a worse and systematically biased energy
472  resolution and should be used with caution. The system can be updated in the future to deal with these
473  cases more carefully. The summed GeV of the hadronic energy is transformed into our reco hadronic energy
474  by a spline fit. This fit is different for the NonQE, QE and combined populations. It is a fit between
475  the summed hadronic energy and (True Neutrino Energy - Reco Muon Energy). More information can
476  be found in the functions NDHadEQE, NDHadENonQE and NDHadECC. Finally, the neutrino energy is just as
477  sum of the Muon Energy and the reco Hadronic Energy.
478 
479  The transition plane is the plane that divides the active region from the beginning of
480  the muon catcher. It is the very last plane in the active region - after that comes
481  steel and planes which are fully encased in the muon catcher.
482 
483  Currently, all the described energy methods are only really meant for contained events; however, this
484  function and the module that calls it doesn't not attempt to determine containment. Please
485  use with discretion!
486  */
487 
488  double muonWeight_act = 0.900903; //For CalCC Energy - Susan's fit from Collaboration Meeting (docdb 7738) on 7/30/2012
489  double muonWeight_tran = 3.84225;
490  double muonWeight_steel = 6.58836;
491  double vertWeight = 1.61335;
492  double offsetWeight = 0.985003;
493 
494  NumuE sliceEnergy;
495 
496  int bestTrack = -1;
497  double trackLen = 0.0;
498 
499  //Need to have valid locations of transition plane!
500  if ((fTranPlane < 0)||(fTranZPos < 0)) return sliceEnergy;
501 
502  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
503  else bestTrack = LongestTrack(sliceTracks, trackLen);
504 
505  if(bestTrack == -1) return sliceEnergy;
506  if(bestTrack == 999) return sliceEnergy;
507 
508  if (fPIDMethod) trackLen = sliceTracks[bestTrack]->TotalLength();
509 
510  //Obtaining values of track length / energy in different regions - written to CAF and can be
511  //used to reoptimize the energy later
512 
513  double recoTrkGeV_act = 0.0;
514  double recoTrkGeV_tran = 0.0;
515  double recoTrkGeV_steel = 0.0;
516  double actTrkLen = 0.0;
517  double catTrkLen = 0.0;
518  double tranX = 999.0;
519  double tranY = 999.0;
520  double hadCalGeV_act = 0.0;
521  double hadCalGeV_tran = 0.0;
522  double hadCalGeV_steel = 0.0;
523  double hadTrkGeV_act = 0.0;
524  double hadTrkGeV_tran = 0.0;
525  double hadTrkGeV_steel = 0.0;
526  double vertexEnergy = 0.0;
527  double vertexTrkEnergy = 0.0;
528 
529  //Need to sum the muon track energy in three different regions
530  NDTrackEnergySplitting(sliceTracks[bestTrack],&recoTrkGeV_act,&recoTrkGeV_tran,&recoTrkGeV_steel);
531 
532  //Need to determine how much track length is in each region
533  NDTrackLength(sliceTracks[bestTrack],&actTrkLen,&catTrkLen,&tranX,&tranY);
534 
535  art::PtrVector< rb::CellHit > trackHits;
536  trackHits = sliceTracks[bestTrack]->AllCells();
537 
538  //Getting all the hits in the slice that are not on our best track
539  rb::Cluster vertexHits = MakeVertexCluster(sliceHits,trackHits);
540 
541  //Gets GeV of hadronic hits NOT on track
542  if (vertexHits.NCell()>0){
543  vertexEnergy = vertexHits.TotalGeV();
544  }
545 
546  //Need to sum the hadron energy in three different regions
547  for (unsigned int hitHad = 0; hitHad <vertexHits.NCell(); ++hitHad){
548  if (vertexHits.RecoHit(hitHad).IsCalibrated()){
549  if (vertexHits.Cell(hitHad)->Plane() < fTranPlane){
550  hadCalGeV_act += vertexHits.RecoHit(hitHad).GeV();
551  }
552  else if (vertexHits.Cell(hitHad)->Plane() == fTranPlane){
553  hadCalGeV_tran += vertexHits.RecoHit(hitHad).GeV();
554  }
555  else {
556  hadCalGeV_steel += vertexHits.RecoHit(hitHad).GeV();
557  }
558  }//End of loop over calibrated nd hadron hits
559  }//End of loop over nd hadron hits
560 
561  //Need hadronic contamination of muon track
562  std::map<int,float> extraHadE = trkCleanUpAlg->ExtraEOnTrackPlanesInGeV(*sliceTracks[bestTrack],*sliceCluster);
563 
564  for (std::map<int,float>::iterator it = extraHadE.begin(); it!=extraHadE.end(); ++it){
565  if (it->first < fTranPlane){
566  hadTrkGeV_act += it->second;
567  }
568  else if (it->first == fTranPlane){
569  hadTrkGeV_tran += it->second;
570  }
571  else {
572  hadTrkGeV_steel += it->second;
573  }
574  }//End of loop over extra hadronic energy
575 
576  vertexTrkEnergy = hadTrkGeV_act + hadTrkGeV_tran + hadTrkGeV_steel;
577 
578 
579  //Calculates GeV of muon from track length in regions of ND
580  double muonGeVFromLen = NDMuonEFromTrackLength(actTrkLen,
581  catTrkLen,
582  recoTrkGeV_act,
583  recoTrkGeV_tran,
584  recoTrkGeV_steel,
585  isRHC);
586 
587  //Gets adjusted GeV of total hadronic energy assuming QE/NonQE
588  double totHadE = vertexEnergy + vertexTrkEnergy;
589  double hadQEGeV = NDHadEQE(totHadE);
590  double hadNonQEGeV = NDHadENonQE(totHadE);
591  double hadCCGeV = NDHadECC(totHadE, isRHC);
592 
593  //***************Fill values in our slice energy object
594 
595  double calcc = muonWeight_act*recoTrkGeV_act + muonWeight_tran*recoTrkGeV_tran + muonWeight_steel*recoTrkGeV_steel
596  + vertWeight*vertexEnergy + offsetWeight;
597  sliceEnergy.SetCalCCE(calcc);
598 
599  double qeE = muonGeVFromLen + hadQEGeV;
600  sliceEnergy.SetTrkQEE(qeE);
601 
602  double nonqeE = muonGeVFromLen + hadNonQEGeV;
603  sliceEnergy.SetTrkNonQEE(nonqeE);
604 
605  double ccE = muonGeVFromLen + hadCCGeV;
606  sliceEnergy.SetTrkCCE(ccE);
607  sliceEnergy.SetRecoMuonE(muonGeVFromLen);
608  sliceEnergy.SetRecoTrkCCHadE(hadCCGeV);
609 
610  sliceEnergy.SetE(ccE); //For now, E is just set to be identical to Trk CC E
611 
612  sliceEnergy.SetNDTrkLenAct(actTrkLen);
613  sliceEnergy.SetNDTrkLenCat(catTrkLen);
614 
615  sliceEnergy.SetNDTrkCalAct(recoTrkGeV_act);
616  sliceEnergy.SetNDTrkCalTran(recoTrkGeV_tran);
617  sliceEnergy.SetNDTrkCalCat(recoTrkGeV_steel);
618 
619  sliceEnergy.SetHadCalE(vertexEnergy);
620  sliceEnergy.SetHadTrkE(vertexTrkEnergy);
621 
622  sliceEnergy.SetNDHadCalAct(hadCalGeV_act);
623  sliceEnergy.SetNDHadCalTran(hadCalGeV_tran);
624  sliceEnergy.SetNDHadCalCat(hadCalGeV_steel);
625  sliceEnergy.SetNDHadTrkAct(hadTrkGeV_act);
626  sliceEnergy.SetNDHadTrkTran(hadTrkGeV_tran);
627  sliceEnergy.SetNDHadTrkCat(hadTrkGeV_steel);
628 
629  sliceEnergy.SetNDTrkTranX(tranX);
630  sliceEnergy.SetNDTrkTranY(tranY);
631 
632  sliceEnergy.SetHadCluster(vertexHits);
633 
634  return sliceEnergy;
635 
636  }//End of SimpleNDEnergy function
637 
638  //----------------------------------------------------------------------
639  NumuE NumuEAlg::Energy(std::vector< art::Ptr<rb::Track> > const sliceTracks,
640  art::PtrVector< rb::CellHit > sliceHits,
641  art::Ptr< rb::Cluster> sliceCluster,
642  art::Event const& e,
643  murem::TrackCleanUpAlg* trkCleanUpAlg,
644  bool isRHC)
645  {
646  if(fND) return NDEnergy(sliceTracks, sliceHits, sliceCluster, e, trkCleanUpAlg, isRHC);
647  if(fFD) return FDEnergy(sliceTracks, sliceHits, sliceCluster, e, trkCleanUpAlg, isRHC);
648 
649  std::cerr << "NumuEAlg::Energy: unknown detector" << std::endl;
650  abort();
651  }
652 
653  //----------------------------------------------------------------------
654  NumuE NumuEAlg::
655  MCTruthEnergyVariables(const std::vector<art::Ptr<rb::Track>>& sliceTracks,
656  const art::Ptr<rb::Cluster>& slice,
657  const std::vector<rb::CellHit>& allHits,
658  const art::Event& e)
659  {
660  /*These variables are useful for making new fits. By writing them to the
661  object, we can fill caf files with them. However, it is dangerous and
662  nasty to touch truth stuff around real analysis work. Please always
663  use with caution and don't do something dumb, like ever use truth
664  information in reconstructed variables. */
665 
666  NumuE truthEnergy;
667 
668  int bestTrack = -1;
669  double trackLen = -1;
670 
671  //Need to have valid locations of transition plane!
672  if (fND && ((fTranPlane < 0)||(fTranZPos < 0))) return truthEnergy;
673 
674  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
675  else bestTrack = LongestTrack(sliceTracks, trackLen);
676 
677  if(bestTrack == -1) return truthEnergy;
678  if(bestTrack == 999) return truthEnergy;
679 
680  art::ServiceHandle<cheat::BackTracker> bt;
681  if(!bt->HaveTruthInfo()){
682  return truthEnergy;
683  }
684 
685  std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slice,allHits);
686  if (funcReturn.empty()) return truthEnergy;
687 
688  std::vector<const sim::Particle*> particleList = bt->MCTruthToParticles(funcReturn[0].neutrinoInt);
689  std::set<int> neutrinoTrackIDs;
690  for (unsigned int i = 0; i < particleList.size(); ++i){
691  neutrinoTrackIDs.insert(particleList[i]->TrackId());
692  }
693 
694  std::set<int>::iterator itr = neutrinoTrackIDs.begin();
695 
696  bool goodMuon = false;
697  double trueCatcherE = -5.0;
698  double trueMuonE = -5.0;
699 
700  while( itr != neutrinoTrackIDs.end() ){
701 
702  int id_int = *itr;
703  const sim::Particle* particle = bt->TrackIDToParticle(id_int);
704  int PDG = particle->PdgCode();
705 
706  if ((PDG != 13)&&(PDG != -13)){
707  ++itr;
708  continue;
709  }//End of loop over non-muons
710 
711  //Trying to find main muon for the event
712  const sim::Particle* mother = bt->TrackIDToMotherParticle(id_int);
713  if (mother->TrackId() != id_int){
714  ++itr;
715  continue;
716  }//End of loop over muons who don't come from the neutrino
717 
718  bool passCuts = bt->PassMuonCuts(id_int,slice->AllCells());
719 
720  if (!passCuts){
721  ++itr;
722  continue;
723  }//End of loop over muons who don't pass cuts
724 
725  goodMuon = true;
726  trueMuonE = particle->E();
727 
728  if (fND){
729 
730  double closeToCenter = 50.0;
731 
732  std::vector<sim::FLSHit> flshits = bt->ParticleToFLSHit(particle->TrackId());
733 
734  for(unsigned int q = 0; q < flshits.size(); ++q){
735  if((flshits[q].GetPlaneID() == fTranPlane) && (std::abs(flshits[q].GetPDG())==13)){
736  if(fabs(flshits[q].GetEntryY()) < closeToCenter){
737  //Adding true muon mass to the energy reported by fls hit
738  trueCatcherE = flshits[q].GetEntryEnergy() + 0.105658;
739  closeToCenter = fabs(flshits[q].GetEntryY());
740  }
741  if(fabs(flshits[q].GetExitY()) < closeToCenter){
742  trueCatcherE = flshits[q].GetExitEnergy() + 0.105658;
743  closeToCenter = fabs(flshits[q].GetExitY());
744  }
745  } // End of loop over transition plane
746  } // End of loop over flshits
747  }//End of loop for ND only
748 
749  break;
750 
751  }//End of while loop over neutrino particles
752 
753  truthEnergy.SetMCTrueMuonE(trueMuonE);
754  truthEnergy.SetMCTrueMuonCatcherE(trueCatcherE);
755  truthEnergy.SetMCGoodTrueMuon(goodMuon);
756 
757  return truthEnergy;
758 
759  }
760 
761  //-------------------------------------------------------------------
762  rb::Cluster NumuEAlg::MakeVertexCluster(art::PtrVector< rb::CellHit > const& sliceHits,
763  art::PtrVector< rb::CellHit > const& trackHits) const
764 
765  {
766  rb::Cluster vertexCluster;
767 
768  for (unsigned int sliceHit = 0; sliceHit < sliceHits.size(); ++sliceHit){
769 
770  bool match = 0;
771 
772  for (unsigned int trackHit = 0; trackHit < trackHits.size(); ++trackHit){
773  match = (sliceHits[sliceHit] == trackHits[trackHit]);
774  if (match) break;
775  }//End of loop over hits in track
776 
777  if (!match) vertexCluster.Add(sliceHits[sliceHit]);
778 
779  }//End of loop over hits in the slice
780 
781  return vertexCluster;
782 
783  }//End of MakeVertexCluster
784 
785  //-------------------------------------------------------------------
786  int NumuEAlg::LongestTrack(std::vector< art::Ptr<rb::Track> > const& sliceTracks,
787  double& longestTrackLength) const
788  {
789  int bestTrack = -1;
790 
791  for(size_t c = 0; c < sliceTracks.size(); ++c){
792 
793  geo::View_t view = sliceTracks[c]->View();
794 
795  //Only use 3D tracks
796  if (view != geo::kXorY){
797  continue;
798  } //End of loop over 2d tracks
799 
800  if ((sliceTracks[c]->TotalLength() > longestTrackLength)&&(sliceTracks[c]->TotalLength() < fMaxTrkLen)) {
801  longestTrackLength = sliceTracks[c]->TotalLength();
802  bestTrack = c;
803  }
804 
805  }//End of for loop over tracks
806 
807  return bestTrack;
808 
809  }//End of longest track
810 
811  //-------------------------------------------------------------------
812  void NumuEAlg::NDTrackLength(art::Ptr<rb::Track> const& track,
813  double* actLength,
814  double* catLength,
815  double* transitionX,
816  double* transitionY) const
817  {
818  if(track->NTrajectoryPoints() == 1) {
819  std::cout<<"You have a problem! Your track only has one trajectory point!"<<std::endl;
820  return;
821  }
822 
823  bool foundActRegion = false;
824  bool foundCatRegion = false;
825 
826  double activeLen = 0.0;
827  double catcherLen = 0.0;
828 
829  TVector3 actTraj;
830  TVector3 catTraj;
831 
832  double tranX = 999.0;
833  double tranY = 999.0;
834 
835  const int N = track->NTrajectoryPoints();
836 
837  for(int n = 0; n < N-1; ++n){
838 
839  if (track->TrajectoryPoint(n+1).Z() <= fTranZPos){
840  foundActRegion = true;
841  activeLen += (track->TrajectoryPoint(n+1)-track->TrajectoryPoint(n)).Mag();
842  }
843  if ((track->TrajectoryPoint(n+1).Z() > fTranZPos)&&!foundCatRegion){
844  foundCatRegion = true;
845  if (track->TrajectoryPoint(n).Z() <= fTranZPos){
846  foundActRegion = true;
847  actTraj = track->TrajectoryPoint(n);
848  catTraj = track->TrajectoryPoint(n+1);
849  }
850  else {
851  catcherLen += (track->TrajectoryPoint(n+1)-track->TrajectoryPoint(n)).Mag();
852  }
853  continue;
854  }
855  if ((track->TrajectoryPoint(n+1).Z() > fTranZPos)&&foundCatRegion){
856  catcherLen += (track->TrajectoryPoint(n+1)-track->TrajectoryPoint(n)).Mag();
857  }
858  }
859 
860  if (foundActRegion && foundCatRegion){
861  TVector3 diff = catTraj-actTraj;
862  double addActLen = (fTranZPos-actTraj.Z())*diff.Mag()/diff.Z();
863  double addCatLen = (catTraj.Z()-fTranZPos)*diff.Mag()/diff.Z();
864  activeLen += addActLen;
865  catcherLen += addCatLen;
866 
867  //Finding x, y location when crosses transition plane
868  tranX = actTraj.X() + (fTranZPos - actTraj.Z())*diff.X()/diff.Z();
869  tranY = actTraj.Y() + (fTranZPos - actTraj.Z())*diff.Y()/diff.Z();
870  }
871 
872  if (actLength) *actLength = activeLen;
873  if (catLength) *catLength = catcherLen;
874  if (transitionX) *transitionX = tranX;
875  if (transitionY) *transitionY = tranY;
876 
877  }//End of NDTrackLength
878 
879  //-------------------------------------------------------------------
880  void NumuEAlg::NDTrackEnergySplitting(art::Ptr<rb::Track> const& track,
881  double* trkCalAct,
882  double* trkCalTran,
883  double* trkCalCat) const
884  {
885  double trkEnergyInActive = 0.0;
886  double trkEnergyInTran = 0.0;
887  double trkEnergyInCatcher = 0.0;
888 
889  //Need to sum the muon track energy in three different regions
890  for (unsigned int hitIdx = 0; hitIdx <track->NCell(); ++hitIdx){
891  if (track->RecoHit(hitIdx).IsCalibrated()){
892  if (track->Cell(hitIdx)->Plane() < fTranPlane){
893  trkEnergyInActive += track->RecoHit(hitIdx).GeV();
894  }
895  else if (track->Cell(hitIdx)->Plane() == fTranPlane){
896  trkEnergyInTran += track->RecoHit(hitIdx).GeV();
897  }
898  else {
899  trkEnergyInCatcher += track->RecoHit(hitIdx).GeV();
900  }
901  }//End of loop over calibrated nd track hits
902  }//End of loop over nd track hits
903 
904  if (trkCalAct) *trkCalAct = trkEnergyInActive;
905  if (trkCalTran) *trkCalTran = trkEnergyInTran;
906  if (trkCalCat) *trkCalCat = trkEnergyInCatcher;
907 
908  }//End of NDTrackEnergySplitting
909 
910  //-------------------------------------------------------------------
911  double NumuEAlg::MuonEFromTrackLength(double trkLen) const
912  {
913 
914  //This function returns muon energy, in GeV, given the muon track length in cm.
915  //It uses a spline fit from May 2015.
916 
917  double muonE = 0.0;
918 
919  if (trkLen <= 0.0) return muonE;
920 
921  if (trkLen < FDMstitch1){
922  muonE = FDMslope1*trkLen + FDMoffset;
923  }
924  else if (trkLen < FDMstitch2){
925  muonE = FDMslope2*trkLen + ((FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset);
926  }
927  else if (trkLen < FDMstitch3){
928  muonE = FDMslope3*trkLen + ((FDMslope2-FDMslope3)*FDMstitch2 + (FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset);
929  }
930  else {
931  muonE = FDMslope4*trkLen + ((FDMslope3-FDMslope4)*FDMstitch3 + (FDMslope2-FDMslope3)*FDMstitch2
932  +(FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset);
933  }
934 
935  return muonE;
936 
937  }//End of MuonEFromTrackLength
938 
939  //-------------------------------------------------------------------
940  double NumuEAlg::FDMuonEFromTrackLength(double trkLen,
941  int run,
942  bool ismc,
943  bool isRHC) const
944  {
945  //Using 2018 energy functions. 2017: predict_prod3_fd_muon_energy(trkLen, run, ismc, isRHC)
946 // return NumuEnergyFunc::predict_prod4_fd_muon_energy(trkLen, run, ismc, isRHC);
947  return NumuEnergyFunc::predict_prod5_fd_muon_energy(trkLen, run, ismc, isRHC);
948  }
949 
950  //-------------------------------------------------------------------
951  double NumuEAlg::ActiveTrackLengthFromMuonE(double muonE) const
952  {
953  /*
954  * TODO: as of 2017 energy estimators, this does not invert the fit in
955  * the correct way.
956  */
957 
958  //This function returns effective active-region track length, in cm, given the muon energy in GeV.
959  //It reverses a spline fit from May 2015.
960 
961  double trkLen = 0.0;
962 
963  if (muonE <= 0.0) return trkLen;
964 
965  if (muonE <= FDMoffset ) return trkLen;
966 
967  else if (muonE < (FDMslope1*FDMstitch1 + FDMoffset)){
968  trkLen = (muonE - FDMoffset)/FDMslope1;
969  }
970  else if (muonE < (FDMslope2*FDMstitch2 + ((FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))){
971  trkLen = (muonE - ((FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope2;
972  }
973  else if (muonE < (FDMslope3*FDMstitch3 + ((FDMslope2-FDMslope3)*FDMstitch2 + (FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))){
974  trkLen = (muonE - ((FDMslope2-FDMslope3)*FDMstitch2 +(FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope3;
975  }
976  else {
977  trkLen = (muonE - ((FDMslope3-FDMslope4)*FDMstitch3 + (FDMslope2-FDMslope3)*FDMstitch2 + (FDMslope1-FDMslope2)*FDMstitch1 + FDMoffset))/FDMslope4;
978  }
979 
980  return trkLen;
981 
982  }//End of ActiveTrackLengthFromMuonE
983  //-------------------------------------------------------------------
984  double NumuEAlg:: NDMuonEFromTrackLength(double actTrkLen,
985  double catTrkLen,
986  double trkCalAct,
987  double trkCalTran,
988  double trkCalCat,
989  bool isRHC,
990  MuonType* muonType) const
991  {
992 
993  //This function returns muon energy, in GeV, given the muon track length in cm.
994  //It distinguishes between the four types of muon tracks and calls the
995  //correct energy estimator for each.
996 
997  double muonE = 0.0;
998 
999  if ((catTrkLen <= 0.0) && (actTrkLen <= 0.0)) return muonE;
1000 
1001  /*
1002  * Numu 2017 Fits
1003  *
1004  * NOTE: apparently there are situations, when there is no calorimetric
1005  * energy deposited in the muon catcher, while there is nonzero
1006  * catcher track length. There are about 10% of such tracks in the
1007  * mc files. None of them pass the numu containment cut though.
1008  *
1009  * NOTE: implementations in CAFAna and here differ in how they treat
1010  * this case. NumuEnergy ignores catcher track if catcalE is 0.
1011  * CAFAna does not, since I believe it is more correct.
1012  * However, the fit was done on the sample with ignored tracks.
1013  *
1014  * NOTE: no 2017 fits were made for the muon, fully contained in the
1015  * catcher. Reason -- they do not pass the numu containment cut.
1016  */
1017 
1018  if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat == 0.0)){
1019  muonE = NDMuonEInActiveOnly(actTrkLen, isRHC);
1020  if (muonType){*muonType = kNDActiveMuon;}
1021  }//Entirely contained in active region and transition plane
1022 
1023  if ((trkCalAct == 0.0)&&(trkCalTran == 0.0)&&(trkCalCat > 0.0)){
1024  muonE = NDMuonEInCatcherOnly(catTrkLen);
1025  if (muonType){*muonType = kNDCatcherMuon;}
1026  }//Entirely contained in catcher region
1027 
1028  if (((trkCalAct > 0.0)||(trkCalTran > 0.0))&&(trkCalCat > 0.0)){
1029  muonE = NDMuonEInActiveAndCatcher(actTrkLen, catTrkLen, isRHC);
1030  if (muonType){*muonType = kNDActCatMuon;}
1031  }//Length in active and catcher region
1032 
1033  return muonE;
1034 
1035  }//End of NDMuonEFromTrackLength
1036  //-------------------------------------------------------------------
1037  double NumuEAlg:: NDMuonEFromTrackLength(art::Ptr<rb::Track> const& track,
1038  bool isRHC,
1039  MuonType* muonType) const
1040  {
1041 
1042  //This function returns muon energy, in GeV, given the muon track length in cm.
1043  //It distinguishes between the four types of muon tracks and calls the
1044  //correct energy estimator for each.
1045 
1046  double muonE = 0.0;
1047 
1048  double actTrkLen = 0.0;
1049  double catTrkLen = 0.0;
1050  double trkCalAct = 0.0;
1051  double trkCalTran = 0.0;
1052  double trkCalCat = 0.0;
1053 
1054  NDTrackLength(track,&actTrkLen,&catTrkLen);
1055  NDTrackEnergySplitting(track,&trkCalAct,&trkCalTran,&trkCalCat);
1056  muonE = NDMuonEFromTrackLength(actTrkLen,
1057  catTrkLen,
1058  trkCalAct,
1059  trkCalTran,
1060  trkCalCat,
1061  isRHC,
1062  muonType);
1063 
1064  return muonE;
1065 
1066  }//End of NDMuonEFromTrackLength
1067  //-------------------------------------------------------------------
1068  double NumuEAlg::NDMuonEInActiveOnly(double actTrkLen, bool isRHC) const
1069  {
1070  //Using 2018 energy functions. 2017: predict_prod3_nd_act_energy(actTrkLen, isRHC)
1071 // return NumuEnergyFunc::predict_prod4_nd_act_energy(actTrkLen, isRHC);
1072  return NumuEnergyFunc::predict_prod5_nd_act_energy(actTrkLen, isRHC);
1073  }
1074  //-------------------------------------------------------------------
1075  double NumuEAlg:: NDMuonEInCatcherOnly(double catTrkLen) const
1076  {
1077 
1078  //This function returns muon energy, in GeV, given the muon track length in cm.
1079  //It is for ND muons that are entirely contained in the muon catcher region of the
1080  //detector. It is a linear fit since the catcher removes our ability to see
1081  //fine structure like we can in the active region. This fit is from July 2013.
1082  //It could use a new tune, but for FA, none of these muons were used so the tuning
1083  //wasn't done.
1084 
1085  double muonE = 0.0;
1086 
1087  if (catTrkLen <= 0.0) return muonE;
1088 
1089  double slope = 0.00513341; // GeV/cm
1090  double offset = 0.230294; // GeV
1091 
1092  muonE = slope*catTrkLen + offset;
1093 
1094  return muonE;
1095 
1096  }//End of NDMuonEInCatcherOnly
1097  //-------------------------------------------------------------------
1098  double NumuEAlg:: NDMuonCatcherEForActiveAndCatcher(double catTrkLen) const
1099  {
1100 
1101  //This function returns muon energy, in GeV, given the muon track length in cm.
1102  //It is for ND muons that have energy in both the active and catcher region.
1103  //This energy is ONLY for the part of the track in the catcher region.
1104  //It is a linear fit since the catcher removes our ability to see
1105  //fine structure like we can in the active region. This fit is from March 2016.
1106  //This is used by ReMId.
1107 
1108  double muonE = 0.0;
1109 
1110  if (catTrkLen <= 0.0) return muonE;
1111 
1112  muonE = NDMCslope*catTrkLen + NDMCoffset;
1113 
1114  return muonE;
1115 
1116  }//End of NDMuonCatcherEForActiveAndCatcher
1117  //-------------------------------------------------------------------
1118  double NumuEAlg::NDMuonCatcherEForActiveAndCatcher(double catTrkLen,
1119  bool isRHC) const
1120  {
1121  //Using 2018 energy functions. 2017 it was predict_prod3_nd_cat_energy
1122 // return NumuEnergyFunc::predict_prod4_nd_cat_energy(catTrkLen, isRHC);
1123  return NumuEnergyFunc::predict_prod5_nd_cat_energy(catTrkLen, isRHC);
1124  }
1125  //-------------------------------------------------------------------
1126  double NumuEAlg::NDMuonEInActiveAndCatcher(double actTrkLen,
1127  double catTrkLen,
1128  bool isRHC) const
1129  {
1130  //Using 2018 energy functions. In 2017 they were predict_prod3_nd_act_energy and predict_prod3_nd_cat_energy
1131 // return NumuEnergyFunc::predict_prod4_nd_act_energy(actTrkLen, isRHC) +
1132 // NumuEnergyFunc::predict_prod4_nd_cat_energy(catTrkLen, isRHC);
1133  return NumuEnergyFunc::predict_prod5_nd_act_energy(actTrkLen, isRHC) + NumuEnergyFunc::predict_prod5_nd_cat_energy(catTrkLen, isRHC);
1134 
1135  }
1136  //-------------------------------------------------------------------
1137  std::vector<rb::Energy> NumuEAlg::MuonEnergies(
1138  std::vector< art::Ptr<rb::Track> > const sliceTracks,
1139  int run,
1140  bool ismc,
1141  bool isRHC
1142  ) const
1143  {
1144 
1145  //This function returns energy for each track, given the assumption
1146  //that the track is a muon.
1147 
1148  std::vector<rb::Energy> energiesOfTracks;
1149 
1150  for(size_t c = 0; c < sliceTracks.size(); ++c){
1151 
1152  rb::Energy trackEnergy(-1.0);
1153 
1154  geo::View_t view = sliceTracks[c]->View();
1155 
1156  //Only give 3D tracks valid energies
1157  if (view != geo::kXorY){
1158  energiesOfTracks.push_back(trackEnergy);
1159  continue;
1160  } //End of loop over 2d tracks
1161 
1162  double energy = -1.0;
1163 
1164  if (fFD){
1165  double trkLen = sliceTracks[c]->TotalLength();
1166  energy = FDMuonEFromTrackLength(trkLen, run, ismc, isRHC);
1167  trackEnergy.SetE(energy);
1168  energiesOfTracks.push_back(trackEnergy);
1169  }//End of FD track
1170 
1171  if (fND){
1172  energy = NDMuonEFromTrackLength(sliceTracks[c], isRHC);
1173  trackEnergy.SetE(energy);
1174  energiesOfTracks.push_back(trackEnergy);
1175  }//End of ND track
1176 
1177  }//End of loop over tracks
1178 
1179  return energiesOfTracks;
1180 
1181  }//End of MuonEnergies
1182  //-------------------------------------------------------------------
1183  int NumuEAlg:: NDTransitionPlane() const
1184  {
1185 
1186  //Function that returns the plane number of the transition plane -
1187  //the last active plane before the steel of the muon catcher.
1188 
1189  return fTranPlane;
1190 
1191  }//End of NDTransitionPlane
1192  //-------------------------------------------------------------------
1193  double NumuEAlg:: NDTransitionZPos() const
1194  {
1195 
1196  //Function that returns the z position of the middle of the transition plane -
1197  //the last active plane before the steel of the muon catcher.
1198 
1199  return fTranZPos;
1200 
1201  }//End of NDTransitionZPos
1202  //-------------------------------------------------------------------
1203  double NumuEAlg:: HadEQE(double calE) const
1204  {
1205 
1206  //This function returns weighted hadron energy, in GeV, given the sum of energy not on
1207  //the muon track and energy from hadron on muon track.
1208  //It uses a spline fit from May 2015. It was fit for a contained QE-only reco population.
1209 
1210  double hadE = 0.0;
1211 
1212  if (calE < 0.0) return hadE;
1213 
1214  double slope1 = 0.621; // Unitless
1215  double slope2 = 1.47; // Unitless
1216  double slope3 = 1.81; // Unitless
1217  double slope4 = 2.06; // Unitless
1218  double offset = 0.0515; // GeV
1219  double stitch1 = 0.0597; // GeV
1220  double stitch2 = 0.139; // GeV
1221  double stitch3 = 0.800; // GeV
1222 
1223  if (calE < stitch1){
1224  hadE = slope1*calE + offset;
1225  }
1226  else if (calE < stitch2){
1227  hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1228  }
1229  else if (calE < stitch3){
1230  hadE = slope3*calE + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset);
1231  }
1232  else {
1233  hadE = slope4*calE + ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2
1234  +(slope1-slope2)*stitch1 + offset);
1235  }
1236 
1237  return hadE;
1238 
1239  }//End of HadEQE
1240  //-------------------------------------------------------------------
1241  double NumuEAlg:: HadENonQE(double calE) const
1242  {
1243 
1244  //This function returns weighted hadron energy, in GeV, given the sum of energy not on
1245  //the muon track and energy from hadron on muon track.
1246  //It uses a spline fit from May 2015. It was fit for a contained NonQE-only reco population.
1247 
1248  double hadE = 0.0;
1249 
1250  if (calE < 0.0) return hadE;
1251 
1252  double slope1 = 0.94; // Unitless
1253  double slope2 = 1.21; // Unitless
1254  double slope3 = 1.78; // Unitless
1255  double slope4 = 2.20; // Unitless
1256  double offset = 0.241; // GeV
1257  double stitch1 = 0.132; // GeV
1258  double stitch2 = 0.223; // GeV
1259  double stitch3 = 0.964; // GeV
1260 
1261  if (calE < stitch1){
1262  hadE = slope1*calE + offset;
1263  }
1264  else if (calE < stitch2){
1265  hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1266  }
1267  else if (calE < stitch3){
1268  hadE = slope3*calE + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset);
1269  }
1270  else {
1271  hadE = slope4*calE + ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2
1272  +(slope1-slope2)*stitch1 + offset);
1273  }
1274 
1275  return hadE;
1276 
1277  }//End of HadENonQE
1278  //-------------------------------------------------------------------
1279  double NumuEAlg::HadECC(double calE, int run, bool ismc, bool isRHC) const
1280  {
1281  //Using 2018 energy functions. In 2017 it was predict_prod3_fd_had_energy
1282 // return NumuEnergyFunc::predict_prod4_fd_had_energy(calE, run, ismc, isRHC);
1283  return NumuEnergyFunc::predict_prod5_fd_had_energy(calE, run, ismc, isRHC);
1284  }//End of HadECC
1285  //-------------------------------------------------------------------
1286  double NumuEAlg:: NDHadEQE(double calE) const
1287  {
1288 
1289  //This function returns weighted hadron energy, in GeV, given the sum of energy not on
1290  //the muon track and energy from hadron on muon track.
1291  //It uses a spline fit from May 2015. It was fit for a contained QE-only reco population.
1292  //The population used for the fit had ALMOST NO hadronic activity in the transition plane or muon catcher!
1293 
1294  double hadE = 0.0;
1295 
1296  if (calE < 0.0) return hadE;
1297 
1298  double slope1 = 1.567; // Unitless
1299  double slope2 = 1.822; // Unitless
1300  double offset = 0.0450; // GeV
1301  double stitch1 = 0.0820; // GeV
1302 
1303  if (calE < stitch1){
1304  hadE = slope1*calE + offset;
1305  }
1306  else {
1307  hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1308  }
1309 
1310  return hadE;
1311 
1312  }//End of NDHadEQE
1313  //-------------------------------------------------------------------
1314  double NumuEAlg:: NDHadENonQE(double calE) const
1315  {
1316 
1317  //This function returns weighted hadron energy, in GeV, given the sum of energy not on
1318  //the muon track and energy from hadron on muon track.
1319  //It uses a spline fit from May 2015. It was fit for a contained nonQE-only reco population.
1320  //The population used for the fit had ALMOST NO hadronic activity in the transition plane or muon catcher!
1321 
1322  double hadE = 0.0;
1323 
1324  if (calE < 0.0) return hadE;
1325 
1326  double slope1 = 1.080; // Unitless
1327  double slope2 = 1.889; // Unitless
1328  double offset = 0.254; // GeV
1329  double stitch1 = 0.169; // GeV
1330 
1331  if (calE < stitch1){
1332  hadE = slope1*calE + offset;
1333  }
1334  else {
1335  hadE = slope2*calE + ((slope1-slope2)*stitch1 + offset);
1336  }
1337 
1338 
1339  return hadE;
1340 
1341  }//End of NDHadENonQE
1342  //-------------------------------------------------------------------
1343  double NumuEAlg::NDHadECC(double calE, bool isRHC) const
1344  {
1345  //Using 2018 energy functions. In 2017 it was predict_prod3_nd_had_energy
1346 // return NumuEnergyFunc::predict_prod4_nd_had_energy(calE, isRHC);
1347  return NumuEnergyFunc::predict_prod5_nd_had_energy(calE, isRHC);
1348  }
1349  //-------------------------------------------------------------------
1350  double NumuEAlg::GetUCE(art::FindManyP<cosrej::CosRejObj> CRO, size_t slicenum, float totalgev, int nhitslice)
1351  {
1352 
1353  std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1354  Evars[0] = totalgev;
1355  Evars[1] = cosmicrej[0]->AngleKal();
1356  Evars[2] = cosmicrej[0]->Scatt();
1357  Evars[3] = cosmicrej[0]->ERatio();
1358  Evars[4] = cosmicrej[0]->KalFwdCell();
1359  Evars[5] = cosmicrej[0]->CosHitRatio()*1.0*nhitslice;
1360  Evars[6] = cosmicrej[0]->NTracks3D();
1361  return fReaderE->EvaluateRegression("BDTG")[0];
1362 
1363  }//End of GetUCE
1364  //-------------------------------------------------------------------
1365  double NumuEAlg::GetUCMuonESingle(std::vector< art::Ptr<rb::Track> > const sliceTracks,
1366  art::PtrVector< rb::CellHit > sliceHits,
1367  art::Event const& e,
1368  art::FindManyP<cosrej::CosRejObj> CRO,
1369  size_t slicenum)
1370  {
1371  int bestTrack = -1;
1372  double trackLen = 0.0;
1373  double hadclustcalE = 0.0;
1374 
1375  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
1376  else bestTrack = LongestTrack(sliceTracks, trackLen);
1377 
1378  if(bestTrack == -1) return -5.0;
1379  if(bestTrack == 999) return -5.0;
1380 
1381  art::PtrVector< rb::CellHit > trackHits;
1382  trackHits = sliceTracks[bestTrack]->AllCells();
1383 
1384  //Getting all the hits in the slice that are not on our best track
1385  rb::Cluster vertexHits = MakeVertexCluster(sliceHits,trackHits);
1386 
1387  if (!(vertexHits.NCell()==0)){
1388  hadclustcalE = vertexHits.CalorimetricEnergy();
1389  }
1390  else{
1391  hadclustcalE = -5.;
1392  }
1393 
1394  std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1395  TMVAvarsSingle[0] = cosmicrej[0]->AngleKal();
1396  TMVAvarsSingle[1] = cosmicrej[0]->TrackLenKal();
1397  TMVAvarsSingle[2] = cosmicrej[0]->TrkEPerNHit();
1398  TMVAvarsSingle[3] = (cosmicrej[0]->Scatt())/(cosmicrej[0]->TrackLenKal());;
1399  TMVAvarsSingle[4] = hadclustcalE;
1400  return fReaderUCMuonSingle->EvaluateRegression("UC_BDTG_S")[0];
1401 
1402  }//End of GetUCMuonESingle
1403  //-------------------------------------------------------------------
1404  double NumuEAlg::GetUCMuonENonSingle(std::vector< art::Ptr<rb::Track> > const sliceTracks,
1405  art::PtrVector< rb::CellHit > sliceHits,
1406  art::Event const& e,
1407  art::FindManyP<cosrej::CosRejObj> CRO,
1408  size_t slicenum)
1409  {
1410  int bestTrack = -1;
1411  double trackLen = 0.0;
1412  double hadclustcalE = 0.0;
1413 
1414  if (fPIDMethod) bestTrack = remid::HighestPIDTrack(sliceTracks, fPIDModuleLabel, e);
1415  else bestTrack = LongestTrack(sliceTracks, trackLen);
1416 
1417  if(bestTrack == -1) return -5.0;
1418  if(bestTrack == 999) return -5.0;
1419 
1420  art::PtrVector< rb::CellHit > trackHits;
1421  trackHits = sliceTracks[bestTrack]->AllCells();
1422 
1423  //Getting all the hits in the slice that are not on our best track
1424  rb::Cluster vertexHits = MakeVertexCluster(sliceHits,trackHits);
1425 
1426  if (!(vertexHits.NCell()==0)){
1427  hadclustcalE = vertexHits.CalorimetricEnergy();
1428  }
1429  else{
1430  hadclustcalE = -5.;
1431  }
1432 
1433  std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1434  TMVAvarsNonSingle[0] = cosmicrej[0]->AngleKal();
1435  TMVAvarsNonSingle[1] = cosmicrej[0]->TrackLenKal();
1436  TMVAvarsNonSingle[2] = cosmicrej[0]->TrkEPerNHit();
1437  TMVAvarsNonSingle[3] = (cosmicrej[0]->Scatt())/(cosmicrej[0]->TrackLenKal());;
1438  TMVAvarsNonSingle[4] = hadclustcalE;
1439  TMVAvarsNonSingle[5] = cosmicrej[0]->TrackLenKal2nd();
1440  TMVAvarsNonSingle[6] = cosmicrej[0]->TrackE2nd();
1441  TMVAvarsNonSingle[7] = cosmicrej[0]->AngleBtwKalTrks();
1442  return fReaderUCMuonNonSingle->EvaluateRegression("UC_BDTG_NS")[0];
1443 
1444  }//End of GetUCMuonESingle
1445  //-------------------------------------------------------------------
1446  double NumuEAlg::GetFSE(art::FindManyP<cosrej::CosRejObj> CRO, size_t slicenum)
1447  {
1448  std::vector< art::Ptr<cosrej::CosRejObj> > cosmicrej = CRO.at(slicenum);
1449  double pars[4];
1450  pars[0] = cosmicrej[0]->FScattExt();
1451  pars[1] = cosmicrej[0]->FScattSig();
1452  pars[2] = cosmicrej[0]->FScattMax();
1453  pars[3] = cosmicrej[0]->FScattSum();
1454  ScatterEnergy *nnEnergy = new ScatterEnergy();
1455  double scatterMuEnergy = nnEnergy->Value(0, pars);
1456  delete nnEnergy;
1457  return scatterMuEnergy;
1458  }//End of GetFSE
1459 
1460 
1461 
1462 
1463  bool NumuEAlg::checkIsMC(art::Event const& evt) const
1464  {
1465  art::ServiceHandle<cheat::BackTracker> bt;
1466  return bt->HaveTruthInfo();
1467  }
1468 
1469 
1470 } // End of namespace
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:232
rb::Track::NTrajectoryPoints
size_t NTrajectoryPoints() const
Definition: Track.h:83
numue::NumuEAlg::fTranZPos
double fTranZPos
Z position of middle of transition plane.
Definition: NumuEAlg.h:181
numue::NumuEAlg::NDTrackEnergySplitting
void NDTrackEnergySplitting(art::Ptr< rb::Track > const &track, double *trkCalAct=NULL, double *trkCalTran=NULL, double *trkCalCat=NULL) const
Function for ND that divides track visible energy into three regions: active, transiton, and catcher.
Definition: NumuEAlg.cxx:880
numue::NumuEAlg::GetUCMuonENonSingle
double GetUCMuonENonSingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
Definition: NumuEAlg.cxx:1404
BackTracker.h
back track the reconstruction to the simulation
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:211
numue::NumuEAlg::NDMuonCatcherEForActiveAndCatcher
double NDMuonCatcherEForActiveAndCatcher(double catTrkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
Definition: NumuEAlg.cxx:1098
plotROC.fileName
fileName
Definition: plotROC.py:78
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
geo::CellGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
rb::Energy::SetE
void SetE(float energy)
Definition: Energy.cxx:21
rb::Track::TrajectoryPoint
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
numue
Energy estimators for CC events.
Definition: FillEnergies.h:7
it
set< int >::iterator it
Definition: ProduceBadRunsByFailureMode.C:101
numue::NumuEAlg::fReaderUCMuonSingle
TMVA::Reader * fReaderUCMuonSingle
Reader for TMVA to get muon uncontained energy.
Definition: NumuEAlg.h:156
numue::NumuEAlg::FDMslope3
double FDMslope3
GeV/cm; Third slope for 4 spline fit for muon trk len -> energy in FD.
Definition: NumuEAlg.h:169
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
numue::NumuEAlg::GetFSE
double GetFSE(art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Fernanda&#39;s scattering energy. Scattering variables themselves stored in CosRej.
Definition: NumuEAlg.cxx:1446
numue::NumuEAlg::FDMstitch3
double FDMstitch3
cm; Third stitch location on 4 spline fit for muon trk len -> energy in FD
Definition: NumuEAlg.h:165
numue::NumuEAlg::Evars
float Evars[8]
Inputs for uncontained energy BDT (TMVA)
Definition: NumuEAlg.h:158
numue::NumuEAlg::FDMslope2
double FDMslope2
GeV/cm; Second slope for 4 spline fit for muon trk len -> energy in FD.
Definition: NumuEAlg.h:168
cheat::BackTracker::TrackIDToMotherParticle
const sim::Particle * TrackIDToMotherParticle(int const &id) const
Definition: BackTracker_service.cc:496
numue::NumuEAlg::HadENonQE
double HadENonQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1241
cheat::BackTracker::SliceToNeutrinoInteractions
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
Definition: BackTracker_service.cc:1062
numue::NumuEAlg::NDHadECC
double NDHadECC(double calE, bool isRHC) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1343
geo::kXorY
X or Y views.
Definition: PlaneGeo.h:30
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
NumuEnergyFunc::predict_prod5_fd_muon_energy
double predict_prod5_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
Definition: Numu2020Fits.cxx:7
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
numue::NumuE::SetNDTrkCalAct
void SetNDTrkCalAct(float ndtrkcalact)
Definition: NumuE.cxx:159
make_mec_shifts_plots.isRHC
string isRHC
Definition: make_mec_shifts_plots.py:100
numue::NumuEAlg::NDHadEQE
double NDHadEQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1286
numue::NumuE::SetMCGoodTrueMuon
void SetMCGoodTrueMuon(bool mcgoodtruemuon)
Definition: NumuE.cxx:229
numue::NumuEAlg::fTMVANameNonSingle
std::string fTMVANameNonSingle
path to TMVA muon weight file
Definition: NumuEAlg.h:152
numue::NumuEAlg::Energy
NumuE Energy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for the current detector. See comments at top of function in ...
Definition: NumuEAlg.cxx:639
numue::NumuEAlg::FDMuonEFromTrackLength
double FDMuonEFromTrackLength(double trkLen, int run, bool ismc, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is done with a spline ...
Definition: NumuEAlg.cxx:940
NumuEnergyFunc::predict_prod5_nd_had_energy
double predict_prod5_nd_had_energy(double hadVisE, bool isRHC)
Definition: Numu2020Fits.cxx:105
geo::PlaneGeo::Cell
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cheat::BackTracker::ParticleToFLSHit
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Definition: BackTracker_service.cc:1514
numue::NumuEAlg::FDEnergy
NumuE FDEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for FD. See comments at top of function in .cxx for full explanati...
Definition: NumuEAlg.cxx:198
numue::NumuE::SetMCTrueMuonE
void SetMCTrueMuonE(float mctruemuone)
Definition: NumuE.cxx:219
numue::NumuE::SetTrkNonQEE
void SetTrkNonQEE(float trknonqe)
Definition: NumuE.cxx:92
track
Definition: event.h:19
rb::Cluster
A collection of associated CellHits.
Definition: Cluster.h:47
numue::NumuEAlg::LongestTrack
int LongestTrack(std::vector< art::Ptr< rb::Track > > const &sliceTracks, double &longestTrackLength) const
Function that finds the longest track that is still smaller than maxTrkLen.
Definition: NumuEAlg.cxx:786
geo::GeometryBase::DetLength
double DetLength() const
Definition: GeometryBase.cxx:481
numue::NumuEAlg::NumuEAlg
NumuEAlg(fhicl::ParameterSet const &pset)
Definition: NumuEAlg.cxx:36
om::cerr
OStream cerr
Definition: OStream.cxx:7
site_stats_from_log.get
get
Definition: site_stats_from_log.py:192
numue::NumuEAlg::NDTransitionPlane
int NDTransitionPlane() const
Function that returns the plane number of the transition plane - the last active plane before the ste...
Definition: NumuEAlg.cxx:1183
NumuEnergyFunc::predict_prod5_nd_cat_energy
double predict_prod5_nd_cat_energy(double ndtrklencat, bool isRHC)
Definition: Numu2020Fits.cxx:93
numue::NumuE::SetRecoTrkCCHadE
void SetRecoTrkCCHadE(float recotrkcchade)
Definition: NumuE.cxx:117
geo::GeometryBase::Plane
const PlaneGeo * Plane(unsigned int i) const
Definition: GeometryBase.cxx:749
numue::NumuEAlg::HadECC
double HadECC(double calE, int run, bool ismc, bool isRHC) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1279
numue::NumuEAlg::NDMCoffset
double NDMCoffset
GeV; Offset for linear fit for muon trk len -> energy in ND muon catcher.
Definition: NumuEAlg.h:172
cet::search_path::find_file
std::string find_file(std::string const &filename) const
numue::NumuEAlg::NDMuonEInActiveAndCatcher
double NDMuonEInActiveAndCatcher(double actTrkLen, double catTrkLen, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
Definition: NumuEAlg.cxx:1126
std::abs
float abs(float number)
Definition: d0nt_math.hpp:39
rb::Cluster::AllCells
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
numue::NumuEAlg::MuonEFromTrackLength
double MuonEFromTrackLength(double trkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is done with a spline ...
Definition: NumuEAlg.cxx:911
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:209
std::isnan
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
numue::NumuEAlg::NDTrackLength
void NDTrackLength(art::Ptr< rb::Track > const &track, double *actLength=NULL, double *catLength=NULL, double *transitionX=NULL, double *transitionY=NULL) const
Function for ND that divides track length into length in active material and length in muon catcher...
Definition: NumuEAlg.cxx:812
rb::Cluster::Add
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
numue::NumuEAlg::FDMoffset
double FDMoffset
GeV; Offset for 4 spline fit for muon trk len -> energy in FD.
Definition: NumuEAlg.h:166
numue::NumuE::SetNDHadCalCat
void SetNDHadCalCat(float ndhadcalcat)
Definition: NumuE.cxx:184
numue::NumuE::SetMCTrueMuonCatcherE
void SetMCTrueMuonCatcherE(float mctruemuoncatchere)
Definition: NumuE.cxx:224
murem::TrackCleanUpAlg::ExtraEOnTrackPlanesInGeV
std::map< int, float > ExtraEOnTrackPlanesInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Definition: TrackCleanUpAlg.cxx:614
numue::NumuEAlg::NDMCslope
double NDMCslope
GeV/cm; First slope for linear fit for muon trk len -> energy in ND muon catcher. ...
Definition: NumuEAlg.h:173
novadaq::cnv::kFARDET
Far Detector at Ash River, MN.
Definition: DAQConventions.h:31
numue::NumuE::SetNDHadCalAct
void SetNDHadCalAct(float ndhadcalact)
Definition: NumuE.cxx:174
submit_syst.x
x
Definition: submit_syst.py:249
numue::NumuEAlg::~NumuEAlg
virtual ~NumuEAlg()
Definition: NumuEAlg.cxx:192
rb::Cluster::CalorimetricEnergy
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
numue::NumuEAlg::fTranPlane
int fTranPlane
Last plane of scintillator before steel begins in ND.
Definition: NumuEAlg.h:180
numue::NumuE::SetNDHadTrkCat
void SetNDHadTrkCat(float ndhadtrkcat)
Definition: NumuE.cxx:199
numue::NumuEAlg::fNumiDir
TVector3 fNumiDir
Direction of NuMI Beam.
Definition: NumuEAlg.h:179
numue::MuonType
MuonType
Definition: NumuE.h:17
numue::NumuEAlg::GetUCMuonESingle
double GetUCMuonESingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
Definition: NumuEAlg.cxx:1365
numue::NumuEAlg::NDMuonEFromTrackLength
double NDMuonEFromTrackLength(double actTrkLen, double catTrkLen, double trkCalAct, double trkCalTran, double trkCalCat, bool isRHC, MuonType *muonType=NULL) const
Function that, given muon track length in cm, returns muon energy in GeV. This distinguishes between ...
Definition: NumuEAlg.cxx:984
geo::GeometryBase::DetId
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
evt
int evt
Definition: FillTrainingTrees.h:22
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
numue::NumuE::SetTrkQEE
void SetTrkQEE(float trkqe)
Definition: NumuE.cxx:87
energy
double energy
Definition: plottest35.C:25
numue::NumuE::SetNDHadTrkAct
void SetNDHadTrkAct(float ndhadtrkact)
Definition: NumuE.cxx:189
rb::Cluster::RecoHit
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
novadaq::cnv::kNEARDET
Near Detector in the NuMI cavern.
Definition: DAQConventions.h:30
numue::NumuEAlg::FDMslope1
double FDMslope1
GeV/cm; First slope for 4 spline fit for muon trk len -> energy in FD.
Definition: NumuEAlg.h:167
numue::NumuEAlg::QEFormulaEnergy
rb::Energy QEFormulaEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, double &error, art::Event const &e, bool isRHC)
Returns QE energy estimation using formula. See comments at top of function in .cxx for full explanat...
Definition: NumuEAlg.cxx:324
numue::NumuE::SetNDTrkTranX
void SetNDTrkTranX(float ndtrktranx)
Definition: NumuE.cxx:209
release_diff.diff
diff
Definition: release_diff.py:80
rb::RecoHit::IsCalibrated
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
numue::NumuEAlg::fPIDModuleLabel
std::string fPIDModuleLabel
Label for module writing ReMId objects to file.
Definition: NumuEAlg.h:147
numue::NumuEAlg::fReaderE
TMVA::Reader * fReaderE
Reader for TMVA to get uncontained energy.
Definition: NumuEAlg.h:155
numue::NumuEAlg::MCTruthEnergyVariables
NumuE MCTruthEnergyVariables(const std::vector< art::Ptr< rb::Track >> &sliceTracks, const art::Ptr< rb::Cluster > &slice, const std::vector< rb::CellHit > &allHits, const art::Event &e)
Returns various useful truth energy values, good for fitting.
Definition: NumuEAlg.cxx:655
numue::NumuE::SetTrkCCE
void SetTrkCCE(float trkcce)
Definition: NumuE.cxx:97
numue::NumuEAlg::fND
bool fND
Is detector ND?
Definition: NumuEAlg.h:176
geo::GeometryBase::DetHalfHeight
double DetHalfHeight() const
Definition: GeometryBase.cxx:474
numue::NumuEAlg::NDTransitionZPos
double NDTransitionZPos() const
Function that returns the z position of the middle of the transition plane - the last active plane be...
Definition: NumuEAlg.cxx:1193
numue::NumuEAlg::checkIsMC
bool checkIsMC(art::Event const &evt) const
Definition: NumuEAlg.cxx:1463
numue::NumuE::SetNDTrkLenCat
void SetNDTrkLenCat(float ndtrklencat)
Definition: NumuE.cxx:154
numue::NumuE::SetNDTrkCalTran
void SetNDTrkCalTran(float ndtrkcaltran)
Definition: NumuE.cxx:164
numue::NumuEAlg::fPIDMethod
bool fPIDMethod
If true, use highest PID to identify muon track instead of longest track (Only FD currently) ...
Definition: NumuEAlg.h:148
run
Definition: run.py:1
art::PtrVector::size
size_type size() const
Definition: PtrVector.h:308
numue::NumuEAlg::FDMstitch2
double FDMstitch2
cm; Second stitch location on 4 spline fit for muon trk len -> energy in FD
Definition: NumuEAlg.h:164
numue::NumuEAlg::ActiveTrackLengthFromMuonE
double ActiveTrackLengthFromMuonE(double muonE) const
Function that, given muon energy in GeV, returns an effective track length in active material in cm...
Definition: NumuEAlg.cxx:951
cheat::BackTracker::TrackIDToParticle
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
Definition: BackTracker_service.cc:481
numue::NumuEAlg::fTMVANameSingle
std::string fTMVANameSingle
path to TMVA muon weight file
Definition: NumuEAlg.h:151
om::cout
OStream cout
Definition: OStream.cxx:6
numue::NumuEAlg::GetUCE
double GetUCE(art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum, float totalgev, int nhitslice)
Get values from cosrej object and use them to get TMVA trained value of uncontained energy for FD...
Definition: NumuEAlg.cxx:1350
PandAna.reco_validation.add_data.offset
offset
Definition: add_data.py:36
numue::NumuEAlg::fReaderUCMuonNonSingle
TMVA::Reader * fReaderUCMuonNonSingle
Reader for TMVA to get muon uncontained energy.
Definition: NumuEAlg.h:157
numue::NumuE::SetNDHadTrkTran
void SetNDHadTrkTran(float ndhadtrktran)
Definition: NumuE.cxx:194
geo::GeometryBase::NuMIBeamDirection
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
Definition: GeometryBase.cxx:523
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
numue::NumuE::SetNDTrkTranY
void SetNDTrkTranY(float ndtrktrany)
Definition: NumuE.cxx:214
numue::NumuEAlg::MakeVertexCluster
rb::Cluster MakeVertexCluster(art::PtrVector< rb::CellHit > const &sliceHits, art::PtrVector< rb::CellHit > const &trackHits) const
Function that makes cluster from all hits in slice that do not belong to track.
Definition: NumuEAlg.cxx:762
rb::RecoHit::GeV
float GeV() const
Definition: RecoHit.cxx:69
numue::NumuE::SetNDTrkLenAct
void SetNDTrkLenAct(float ndtrklenact)
Definition: NumuE.cxx:149
geo::GeometryBase::DetHalfWidth
double DetHalfWidth() const
Definition: GeometryBase.cxx:467
pars
std::string pars("Th23Dmsq32")
rb::Energy
A container for energy information.
Definition: Energy.h:20
numue::NumuEAlg::fRun
int fRun
What run is it?
Definition: NumuEAlg.h:177
util::pythag
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
numue::NumuEAlg::FDMslope4
double FDMslope4
GeV/cm; Fourth slope for 4 spline fit for muon trk len -> energy in FD.
Definition: NumuEAlg.h:170
numue::NumuEAlg::fFD
bool fFD
Is detector FD?
Definition: NumuEAlg.h:175
numue::NumuE::SetHadTrkE
void SetHadTrkE(float hadtrk)
Definition: NumuE.cxx:127
numue::NumuEAlg::NDMuonEInActiveOnly
double NDMuonEInActiveOnly(double actTrkLen, bool isRHC) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
Definition: NumuEAlg.cxx:1068
numue::NumuE::SetCalCCE
void SetCalCCE(float calcc)
Definition: NumuE.cxx:82
numue::NumuE::SetNDHadCalTran
void SetNDHadCalTran(float ndhadcaltran)
Definition: NumuE.cxx:179
geom
void geom(int which=0)
Definition: geom.C:163
cheat::BackTracker::HaveTruthInfo
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
Mag
float Mag() const
Definition: SRLorentzVectorAddons.h:10
rb::Cluster::TotalGeV
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
search_path.h
demo.particleList
particleList
Definition: demo.py:41
numue::NumuEAlg::NDEnergy
NumuE NDEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for ND. See comments at top of function in .cxx for full explanati...
Definition: NumuEAlg.cxx:423
numue::NumuE::SetHadCluster
void SetHadCluster(rb::Cluster hadcluster)
Definition: NumuE.cxx:234
numue::NumuEAlg::fMaxTrkLen
double fMaxTrkLen
Maximum track length allowed in detector (diagonal length)
Definition: NumuEAlg.h:178
cheat::BackTracker::PassMuonCuts
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
Definition: BackTracker_service.cc:965
numue::NumuEAlg::HadEQE
double HadEQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1203
NumuEnergyFunc::predict_prod5_fd_had_energy
double predict_prod5_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
Definition: Numu2020Fits.cxx:46
e
Float_t e
Definition: plot.C:35
art::Event::run
RunNumber_t run() const
Definition: Event.h:77
numue::NumuE::SetRecoMuonE
void SetRecoMuonE(float recomuone)
Definition: NumuE.cxx:112
numue::NumuEAlg::TMVAvarsSingle
float TMVAvarsSingle[38]
Inputs for muon uncontained energy BDT (TMVA) for single track events.
Definition: NumuEAlg.h:160
numue::NumuE::SetHadCalE
void SetHadCalE(float hadcal)
Definition: NumuE.cxx:122
numue::NumuE::SetNDTrkCalCat
void SetNDTrkCalCat(float ndtrkcalcat)
Definition: NumuE.cxx:169
numue::NumuEAlg::TMVAvarsNonSingle
float TMVAvarsNonSingle[38]
Inputs for muon uncontained energy BDT (TMVA) for non-single track events.
Definition: NumuEAlg.h:161
NumuEnergyFunc::predict_prod5_nd_act_energy
double predict_prod5_nd_act_energy(double ndtrklenact, bool isRHC)
Definition: Numu2020Fits.cxx:81
numue::NumuEAlg::FDMstitch1
double FDMstitch1
cm; First stitch location on 4 spline fit for muon trk len -> energy in FD
Definition: NumuEAlg.h:163
murem::TrackCleanUpAlg::ExtraEOnTrackInGeV
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
Definition: TrackCleanUpAlg.cxx:571
numue::NumuEAlg::fTMVAName
std::string fTMVAName
path to TMVA weight file
Definition: NumuEAlg.h:149
numue::ScatterEnergy::Value
double Value(int index, double in0, double in1, double in2, double in3)
Definition: ScatterEnergy.cxx:17
numue::NumuEAlg::NDHadENonQE
double NDHadENonQE(double calE) const
Function that, given total hadron cal energy on and off muon track in GeV, returns weighted hadron en...
Definition: NumuEAlg.cxx:1314
cet::search_path
Definition: search_path.h:37
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
remid::HighestPIDTrack
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
numue::NumuEAlg::MuonEnergies
std::vector< rb::Energy > MuonEnergies(std::vector< art::Ptr< rb::Track > > const sliceTracks, int run, bool ismc, bool isRHC) const
Function that returns energy for each track given the assumption that it is a muon. It uses the functions above, like MuonEFromTrackLength and NDMuonEFromTrackLength.
Definition: NumuEAlg.cxx:1137
numue::NumuEAlg::NDMuonEInCatcherOnly
double NDMuonEInCatcherOnly(double catTrkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
Definition: NumuEAlg.cxx:1075
cheat::BackTracker::MCTruthToParticles
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
Definition: BackTracker_service.cc:527
geo::GeometryBase::FirstPlaneInMuonCatcher
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
Definition: GeometryBase.cxx:979
string
enum BeamMode string
Generated on Sun Feb 28 2021 14:13:24 for NOvA by   doxygen 1.8.11