NuMuAnalysisSetup.cxx
Go to the documentation of this file.
1 /*
2  * \file NuMuAnalysisSetup.cxx
3  * \brief analysis configuration for nu_mu analysis, First Analysis version
4  *
5  * Created on: \date Feb 18, 2016
6  * Original author: \author J. Wolcott <jwolcott@fnal.gov>
7  *
8  */
9 
11 
17 
18 #include "RecoBase/Cluster.h"
19 #include "RecoBase/Energy.h"
20 #include "RecoBase/Track.h"
21 #include "NumuEnergy/NumuE.h"
24 #include "ReMId/ReMId.h"
25 
26 #include "FNEX/core/VarVals.h"
27 
32 
33 namespace fnex
34 {
35  //----------------------------------------------------------------------------
37  {
38  this->Reconfigure(pset);
39  }
40 
41  //----------------------------------------------------------------------------
43  art::ValidHandle<std::vector<rb::Cluster> > const& slices,
44  art::InputTag const& tag,
45  std::vector<fnex::DataVarVals> & dataVarVals,
46  std::vector<fnex::MCVarVals> & mcVarVals,
47  std::vector<fnex::EventId> & evIds,
48  fnex::BeamType_t const& beamType ) const
49  {
50  dataVarVals.clear();
51  mcVarVals .clear();
52  evIds .clear();
53 
54  if( !slices.isValid() ){
55  LOG_WARNING("NuMuAnalysisSetup")
56  << "could not find valid handle of rb::Cluster, ignore this event";
57  }
58  else{
59 
60  LOG_DEBUG("NuMuAnalysisSetup")
61  << "Start filling variables for tag: "
62  << tag
63  << " and quantile "
65 
66  // get the FindOnes
67  art::FindOne <numue::NumuE > fone (slices, evt, tag);
68  art::FindManyP<rb::Track > fmt (slices, evt, tag);
69  art::FindOneP <simb::MCTruth > fomct(slices, evt, tag);
70  art::FindOne <simb::MCFlux > fomcf(slices, evt, tag);
71  art::FindOne <fxwgt::FluxWeights> fofw (slices, evt, tag);
72  art::FindOne <simb::GTruth > fogt (slices, evt, tag);
73 
74  std::vector< art::Ptr<rb::Track> > trackPtrs;
75 
76  size_t bestPIDTrack = std::numeric_limits<size_t>::max();
77 
78  for(size_t sl = 0; sl < slices->size(); ++sl){
79  fnex::DataVarVals dataVars;
80  fnex::MCVarVals mcVars;
81  fnex::EventId evId(evt.run(),
82  evt.subRun(),
83  evt.event(),
85  (*slices)[sl].ID());
86 
87  // check that the FindOnes are valid
88  // only need the 0th entry in the FindOnes because we have only one slice here
89  if(fone.isValid() ){
90 
91  cet::maybe_ref<numue::NumuE const> nume = fone.at(sl);
92 
93  fmt.get(sl, trackPtrs);
94  if(trackPtrs.size() < 1) continue;
95 
96  // get the best pid for the slice - this index will not be the
97  // same as in the preskimmed slice
98  bestPIDTrack = remid::HighestPIDTrack(trackPtrs, tag, evt);
99  //bestPIDTrack = this->HighestPIDTrack(trackPtrs, tag, evt); // Reimplement here because Is3D does not work
100  // and the remid default does not exclude 2d tracks for Skimmed files
101 
102  // this is not terribly efficient, as we are making the FindOneP
103  // for every slice in the event, but I don't see a good way around
104  // that while still getting the bookkeeping right
105  art::FindOne<rb::Energy> foe(trackPtrs, evt, tag);
106 
107  if( nume.isValid() && foe.isValid() )
108  this->FillRecoVars(dataVars, nume.ref(), foe.at(bestPIDTrack).ref(), *(trackPtrs[bestPIDTrack]), !evt.isRealData(), beamType);
109 
110  }// end filling recovariables
111  else{
112  LOG_WARNING("NuMuAnalysisSetup")
113  << "something is wrong with the FindOne<numue::NumuE> "
114  << " skip this slice";
115  continue;
116  }
117 
118  // now for the MC variables, if this is a MC event
119  if( !this->FillTruthVars(mcVars, evt, sl, fomct, fomcf, fofw, fogt) ) continue;
120 
121  LOG_DEBUG("NuMuAnalysisSetup")
122  << "numu truth is "
123  << mcVars;
124 
125  // check the hadronic eenrgy fraction quantile and only add the vars to the collection
126  // if we are in the desired range (expects 0 for SA
127  double lepE = dataVars.val_at(fnex::kLep_RecoE, fnex::MetaData());
128  double hadE = dataVars.val_at(fnex::kHad_RecoE, fnex::MetaData());
129 
130  int hadEFracQ = FourthAnalysisQuantile(lepE, hadE, beamType);
131  if( fHadEFracQuantile == 0 || hadEFracQ == fHadEFracQuantile ){
132  dataVarVals.push_back(dataVars);
133  mcVarVals .push_back(mcVars);
134  evIds .push_back(evId);
135  }
136  } // end loop over slices
137 
138  } // end if the handle was valid
139 
140  return;
141  } // NuMuAnalysisSetup::FillVars
142 
143  //----------------------------------------------------------------------------
145  {
146  // call base class to get Cuts & shifters registered
147  this->AnalysisSetupBase::Reconfigure(pset);
148  fNuMuEEstimator = pset.get<std::string>("NuMuEEstimator", "SA");
149 
150  // get the selection type and determine the hadronic fraction qunatile
151  auto const& type = pset.get<std::string>("AnalysisType");
152 
153  fHadEFracQuantile = 0; // defalut for selections prior to TA
154 
155  if (type.compare("NuMu_Q1") == 0){
156  fHadEFracQuantile = 1;
157  }
158  else if(type.compare("NuMu_Q2") == 0){
159  fHadEFracQuantile = 2;
160  }
161  else if(type.compare("NuMu_Q3") == 0){
162  fHadEFracQuantile = 3;
163  }
164  else if(type.compare("NuMu_Q4") == 0){
165  fHadEFracQuantile = 4;
166  }
167  }
168 
169  //----------------------------------------------------------------------------
171  numue::NumuE const& numue,
172  rb::Energy const& energy,
173  rb::Track const& track,
174  bool isMC,
175  fnex::BeamType_t const& beamType ) const
176  {
177  // lepton energy is muon energy
178  double hadE = (numue.TrkCCE() - energy.E());
179  double muonE = energy.E();
180  if(!fNuMuEEstimator.compare("SA"))
181  this->SecondAnalysisEnergy(numue, track, hadE, muonE);
182  else if(!fNuMuEEstimator.compare("TA"))
183  this->ThirdAnalysisEnergy(numue, track, hadE, muonE, isMC, beamType);
184  else if(!fNuMuEEstimator.compare("2018")){
185  LOG_DEBUG("NuMuAnalysisSetup") << "Using 2018 numu energy estimator";
186  this->Analysis2018Energy(numue, track, hadE, muonE, isMC, beamType);
187  }
188  else
189  throw cet::exception("NuMuAnalysisSetup:FillRecoVars : Need to set NuMuEEstimator to SA or TA (for 2nd and 3rd analysis respectively)");
190 
191  vars.set_val_at(fnex::kLep_RecoE, muonE);
192  vars.set_val_at(fnex::kHad_RecoE, hadE );
193 
194  // fraction of leptonic energy in muon catcher for ND
195  auto det = fDetService->DetId();
196  //auto run = fRH->RunNumber();
197 
198  if(det == novadaq::cnv::kNEARDET) {
199 
200  bool isRHC = false; // FHC mode
201  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
202 
203  double lepEAct = predict_prod3_nd_act_energy(numue.NDTrkLenAct(), isRHC);
204  double lepECat = predict_prod3_nd_cat_energy(numue.NDTrkLenCat(), isRHC);
205  double lepMCFrac = lepECat / (lepEAct + lepECat);
206  vars.set_val_at(fnex::kLep_RecoE_MCFrac, lepMCFrac);
207  }
208  else vars.set_val_at(fnex::kLep_RecoE_MCFrac, 0);
209 
210  // reconstructed Q^2
211  double M_mu_sqr = std::pow(0.1056, 2);
212  double p_mu = std::sqrt(std::pow(energy.E(), 2) - M_mu_sqr);
213 
214  TVector3 beamDir;
216  beamDir = TVector3(-8.424e-04,
217  -6.174e-02,
218  +9.981e-1).Unit();
219  else if(det == novadaq::cnv::kFARDET)
220  beamDir = TVector3(-6.833e-05,
221  +6.388e-02,
222  +9.980e-1).Unit();
223 
226 
227  double RecoQ2 = 2 * vars.val_at(fnex::kNu_RecoE, md) * (energy.E() - p_mu * track.Dir().Dot(beamDir) - M_mu_sqr);
228 
229  vars.set_val_at(fnex::kRecoQ2, RecoQ2);
230 
231  return;
232  }
233 
234  //----------------------------------------------------------------------------
236  rb::Track const& track,
237  double & hadE,
238  double & muonE) const
239  {
240  auto det = fDetService->DetId();
241  auto run = fRH->RunNumber();
242 
243  //New muon tuning values
244  double stitch1 = 342; // cm
245  double stitch2 = 520; // cm
246  double stitch3 = 1090; // cm
247  double offset = 0.1491; // GeV
248  double slope1 = 0.001951; // GeV/cm
249  double slope2 = 0.002022; // GeV/cm
250  double slope3 = 0.002067; // GeV/cm
251  double slope4 = 0.002196; // GeV/cm
252 
253  double trkLen = track.TotalLength();
254  double hadVisE = numue.HadTrkE() + numue.HadCalE();
255 
256  if(det == novadaq::cnv::kFARDET){
257  if (trkLen <= 0.0){
258  muonE = 0.;
259  hadE = 0.;
260  return;
261  }
262 
263  if (trkLen < stitch1){
264  muonE = slope1*trkLen + offset;
265  }
266  else if (trkLen < stitch2){
267  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
268  }
269  else if (trkLen <stitch3){
270  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
271  }
272  else{
273  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
274  }
275 
276  //New Tune Values with Error Rounding - Period 1
277  double stitchH1 = 0.0799; // GeV
278  double stitchH2 = 0.373; // GeV
279  double stitchH3 = 0.925; // GeV
280  double offsetH = 0.0419; // GeV
281  double slopeH1 = 0.816; // unitless
282  double slopeH2 = 2.04; // unitless
283  double slopeH3 = 1.68; // unitless
284  double slopeH4 = 2.45; // unitless
285 
286  if(run >= 17140 && run < 20753){
287  //New Tune Values with Error Rounding - Period 2
288  stitchH1 = 0.0694; // GeV
289  stitchH2 = 0.288; // GeV
290  stitchH3 = 0.930; // GeV
291  offsetH = 0.0392; // GeV
292  slopeH1 = 1.45; // unitless
293  slopeH2 = 2.13; // unitless
294  slopeH3 = 1.76; // unitless
295  slopeH4 = 2.17; // unitless
296  }
297  else if(run >= 20753){
298  //New Tune Values with Error Rounding - Period 3
299  stitchH1 = 0.0702; // GeV
300  stitchH2 = 0.321; // GeV
301  stitchH3 = 1.06; // GeV
302  offsetH = 0.0365; // GeV
303  slopeH1 = 1.34; // unitless
304  slopeH2 = 2.05; // unitless
305  slopeH3 = 1.68; // unitless
306  slopeH4 = 2.28; // unitless
307  }
308 
309  if (hadVisE < 0.0){
310  hadE = offsetH;
311  }
312  else if (hadVisE < stitchH1){
313  hadE = slopeH1*hadVisE + offsetH;
314  }
315  else if (hadVisE < stitchH2){
316  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
317  }
318  else if (hadVisE <stitchH3){
319  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
320  }
321  else{
322  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
323  }
324  }//FD energy
325 
327  double trkLenAct = numue.NDTrkLenAct();
328  double trkLenCat = numue.NDTrkLenCat();
329 
330  if(trkLenAct <= 0.0 && trkLenCat <= 0.0){
331  muonE = 0.;
332  hadE = 0.;
333  return;
334  }
335 
336  bool actOnly = ((numue.NDTrkCalAct() > 0.0 ||
337  numue.NDTrkCalTran() > 0.0 ) &&
338  numue.NDTrkCalCat() == 0.0 );
339  bool actAndCat = ((numue.NDTrkCalAct() > 0.0 ||
340  numue.NDTrkCalTran() > 0.0 ) &&
341  numue.NDTrkCalCat() > 0.0 );
342 
343  if( !actOnly && !actAndCat){
344  muonE = 0.;
345  hadE = 0.;
346  return;
347  }
348  if( actOnly && actAndCat){
349  muonE = 0.;
350  hadE = 0.;
351  }
352 
353  if(actOnly){
354  trkLen = trkLenAct;
355  }
356 
357  if(actAndCat){
358 
359  double muonEMC = 0.0;
360 
361  //New tuning values rounded for errors
362  double offsetMC = 0.139; // GeV
363  double slope1MC = 0.00537; // GeV/cm
364 
365  muonEMC = slope1MC*trkLenCat + offsetMC;
366 
367  double effActTrkLen = 0.0;
368 
369  if(muonEMC <= 0.0){
370  muonE = 0.;
371  hadE = 0.;
372  return;
373  }
374 
375  if(muonEMC <= offset) {
376  effActTrkLen = 0.0;
377  }
378  else if(muonEMC < (slope1*stitch1 + offset)){
379  effActTrkLen = (muonEMC - offset)/slope1;
380  }
381  else if(muonEMC < (slope2*stitch2 + ((slope1-slope2)*stitch1 + offset))){
382  effActTrkLen = (muonEMC - ((slope1-slope2)*stitch1 + offset))/slope2;
383  }
384  else if(muonEMC < (slope3*stitch3 + ((slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))){
385  effActTrkLen = (muonEMC - ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset))/slope3;
386  }
387  else{
388  effActTrkLen = (muonEMC - ((slope3-slope4)*stitch3 + (slope2-slope3)*stitch2 + (slope1-slope2)*stitch1 + offset))/slope4;
389  }
390 
391  trkLen = trkLenAct + effActTrkLen;
392 
393  }//End of actAndCat
394 
395  //Now get final muon energy from total "active" track length
396  if(trkLen < stitch1){
397  muonE = slope1*trkLen + offset;
398  }
399  else if (trkLen < stitch2){
400  muonE = slope2*trkLen + ((slope1-slope2)*stitch1 + offset);
401  }
402  else if (trkLen <stitch3){
403  muonE = slope3*trkLen + ((slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
404  }
405  else{
406  muonE = slope4*trkLen + ((slope3-slope4)*stitch3 +(slope2-slope3)*stitch2 +(slope1-slope2)*stitch1 + offset);
407  }
408 
409 
410  //New Tune Values with Error Rounding
411  double stitchH1 = 0.0394; // GeV
412  double stitchH2 = 0.362; // GeV
413  double stitchH3 = 0.955; // GeV
414  double offsetH = 0.0534; // GeV
415  double slopeH1 = 0.956; // unitless
416  double slopeH2 = 2.14; // unitless
417  double slopeH3 = 1.78; // unitless
418  double slopeH4 = 2.00; // unitless
419 
420  if (hadVisE < 0.0){
421  hadE = offsetH;
422  }
423  else if (hadVisE < stitchH1){
424  hadE = slopeH1*hadVisE + offsetH;
425  }
426  else if (hadVisE < stitchH2){
427  hadE = slopeH2*hadVisE + ((slopeH1-slopeH2)*stitchH1 + offsetH);
428  }
429  else if (hadVisE <stitchH3){
430  hadE = slopeH3*hadVisE + ((slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
431  }
432  else{
433  hadE = slopeH4*hadVisE + ((slopeH3-slopeH4)*stitchH3 +(slopeH2-slopeH3)*stitchH2 +(slopeH1-slopeH2)*stitchH1 + offsetH);
434  }
435 
436  }//ND energy
437 
438  // START TERRIBLE HACK FOR SYSTEMATICS
439 // if(tMuE > 0) muonE = tMuE;
440 // if(tHadE > 0) hadE = tHadE;
441 // double nuE = muonE + hadE;
442 // if(tNuEScale > 0) nuE *= tNuEScale;
443 // return nuE;
444 
445  }
446 
447  //----------------------------------------------------------------------------
449  rb::Track const& track,
450  double & hadE,
451  double & muonE,
452  bool isMC,
453  fnex::BeamType_t const& beamType ) const
454  {
455  auto det = fDetService->DetId();
456  auto run = fRH->RunNumber();
457 
458  double hadVisE = numue.HadTrkE() + numue.HadCalE();
459 
460  bool isRHC = false; // FHC mode
461  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
462 
463 
465  {
466 
467  double trkLen = track.TotalLength();
468 
469  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
470  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
471 
472  muonE = predict_prod3_fd_muon_energy(trkLen, run, isMC, isRHC);
473  hadE = predict_prod3_fd_had_energy(hadVisE, run, isMC, isRHC);
474 
475  LOG_DEBUG("NuMuAnalysisSetup") << "LepE = " << muonE << "\t HadE = " << hadE;
476  }
477 
479  {
480 
481  double trkLenAct = numue.NDTrkLenAct();
482  double trkLenCat = numue.NDTrkLenCat();
483 
484  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
485  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
486 
487  muonE = predict_prod3_nd_act_energy(trkLenAct, isRHC) + predict_prod3_nd_cat_energy(trkLenCat, isRHC);
488  hadE = predict_prod3_nd_had_energy(hadVisE, isRHC);
489  }
490 
491  }// end TA energy estimator
492 
493 
494  //----------------------------------------------------------------------------
496  double & hadE) const
497  {
498 
499  double recoE = lepE + hadE;
500  if( recoE == 0 ) {
501  LOG_DEBUG("NuMuAnalysisSetup")
502  << "cannot compute quantile: total energy is zero";
503  return -1;
504  }
505 
506  double hadEFrac = hadE / ( lepE + hadE );
507 
508  if( recoE >= 0 && recoE < 0.75 ){
509  if ( hadEFrac >= 0 && hadEFrac < 0.160864) return 1;
510  else if( hadEFrac >= 0.160864 && hadEFrac < 0.231946) return 2;
511  else if( hadEFrac >= 0.231946 && hadEFrac < 0.294404) return 3;
512  else if( hadEFrac >= 0.294404 && hadEFrac <= 1 ) return 4;
513  }
514  else if( recoE >= 0.75 && recoE < 1 ){
515  if ( hadEFrac >= 0 && hadEFrac < 0.201787) return 1;
516  else if( hadEFrac >= 0.201787 && hadEFrac < 0.290059) return 2;
517  else if( hadEFrac >= 0.290059 && hadEFrac < 0.376756) return 3;
518  else if( hadEFrac >= 0.376756 && hadEFrac <= 1 ) return 4;
519  }
520  else if( recoE >= 1 && recoE < 1.1 ){
521  if ( hadEFrac >= 0 && hadEFrac < 0.209212) return 1;
522  else if( hadEFrac >= 0.209212 && hadEFrac < 0.306823) return 2;
523  else if( hadEFrac >= 0.306823 && hadEFrac < 0.404893) return 3;
524  else if( hadEFrac >= 0.404893 && hadEFrac <= 1 ) return 4;
525  }
526  else if( recoE >= 1.1 && recoE < 1.2 ){
527  if ( hadEFrac >= 0 && hadEFrac < 0.214098) return 1;
528  else if( hadEFrac >= 0.214098 && hadEFrac < 0.314586) return 2;
529  else if( hadEFrac >= 0.314586 && hadEFrac < 0.417883) return 3;
530  else if( hadEFrac >= 0.417883 && hadEFrac <= 1 ) return 4;
531  }
532  else if( recoE >= 1.2 && recoE < 1.3 ){
533  if ( hadEFrac >= 0 && hadEFrac < 0.209203) return 1;
534  else if( hadEFrac >= 0.209203 && hadEFrac < 0.316993) return 2;
535  else if( hadEFrac >= 0.316993 && hadEFrac < 0.428731) return 3;
536  else if( hadEFrac >= 0.428731 && hadEFrac <= 1 ) return 4;
537  }
538  else if( recoE >= 1.3 && recoE < 1.4 ){
539  if ( hadEFrac >= 0 && hadEFrac < 0.203099) return 1;
540  else if( hadEFrac >= 0.203099 && hadEFrac < 0.315997) return 2;
541  else if( hadEFrac >= 0.315997 && hadEFrac < 0.432587) return 3;
542  else if( hadEFrac >= 0.432587 && hadEFrac <= 1 ) return 4;
543  }
544  else if( recoE >= 1.4 && recoE < 1.5 ){
545  if ( hadEFrac >= 0 && hadEFrac < 0.198472) return 1;
546  else if( hadEFrac >= 0.198472 && hadEFrac < 0.31283 ) return 2;
547  else if( hadEFrac >= 0.31283 && hadEFrac < 0.442388) return 3;
548  else if( hadEFrac >= 0.442388 && hadEFrac <= 1 ) return 4;
549  }
550  else if( recoE >= 1.5 && recoE < 1.6 ){
551  if ( hadEFrac >= 0 && hadEFrac < 0.188711) return 1;
552  else if( hadEFrac >= 0.188711 && hadEFrac < 0.302839) return 2;
553  else if( hadEFrac >= 0.302839 && hadEFrac < 0.442194) return 3;
554  else if( hadEFrac >= 0.442194 && hadEFrac <= 1 ) return 4;
555  }
556  else if( recoE >= 1.6 && recoE < 1.7 ){
557  if ( hadEFrac >= 0 && hadEFrac < 0.180827) return 1;
558  else if( hadEFrac >= 0.180827 && hadEFrac < 0.293246) return 2;
559  else if( hadEFrac >= 0.293246 && hadEFrac < 0.443012) return 3;
560  else if( hadEFrac >= 0.443012 && hadEFrac <= 1 ) return 4;
561  }
562  else if( recoE >= 1.7 && recoE < 1.8 ){
563  if ( hadEFrac >= 0 && hadEFrac < 0.172282) return 1;
564  else if( hadEFrac >= 0.172282 && hadEFrac < 0.283574) return 2;
565  else if( hadEFrac >= 0.283574 && hadEFrac < 0.441167) return 3;
566  else if( hadEFrac >= 0.441167 && hadEFrac <= 1 ) return 4;
567  }
568  else if( recoE >= 1.8 && recoE < 1.9 ){
569  if ( hadEFrac >= 0 && hadEFrac < 0.16746) return 1;
570  else if( hadEFrac >= 0.16746 && hadEFrac < 0.27588) return 2;
571  else if( hadEFrac >= 0.27588 && hadEFrac < 0.43916) return 3;
572  else if( hadEFrac >= 0.43916 && hadEFrac <= 1 ) return 4;
573  }
574  else if( recoE >= 1.9 && recoE < 2 ){
575  if ( hadEFrac >= 0 && hadEFrac < 0.1646 ) return 1;
576  else if( hadEFrac >= 0.1646 && hadEFrac < 0.272764) return 2;
577  else if( hadEFrac >= 0.272764 && hadEFrac < 0.444039) return 3;
578  else if( hadEFrac >= 0.444039 && hadEFrac <= 1 ) return 4;
579  }
580  else if( recoE >= 2 && recoE < 2.25 ){
581  if ( hadEFrac >= 0 && hadEFrac < 0.165357) return 1;
582  else if( hadEFrac >= 0.165357 && hadEFrac < 0.27372) return 2;
583  else if( hadEFrac >= 0.27372 && hadEFrac < 0.451054) return 3;
584  else if( hadEFrac >= 0.451054 && hadEFrac <= 1 ) return 4;
585  }
586  else if( recoE >= 2.25 && recoE < 2.5 ){
587  if ( hadEFrac >= 0 && hadEFrac < 0.169381) return 1;
588  else if( hadEFrac >= 0.169381 && hadEFrac < 0.282544) return 2;
589  else if( hadEFrac >= 0.282544 && hadEFrac < 0.460158) return 3;
590  else if( hadEFrac >= 0.460158 && hadEFrac <= 1 ) return 4;
591  }
592  else if( recoE >= 2.5 && recoE < 2.75 ){
593  if ( hadEFrac >= 0 && hadEFrac < 0.165923) return 1;
594  else if( hadEFrac >= 0.165923 && hadEFrac < 0.283755) return 2;
595  else if( hadEFrac >= 0.283755 && hadEFrac < 0.456392) return 3;
596  else if( hadEFrac >= 0.456392 && hadEFrac <= 1 ) return 4;
597  }
598  else if( recoE >= 2.75 && recoE < 3 ){
599  if ( hadEFrac >= 0 && hadEFrac < 0.158863) return 1;
600  else if( hadEFrac >= 0.158863 && hadEFrac < 0.283625) return 2;
601  else if( hadEFrac >= 0.283625 && hadEFrac < 0.45765 ) return 3;
602  else if( hadEFrac >= 0.45765 && hadEFrac <= 1 ) return 4;
603  }
604  else if( recoE >= 3 && recoE < 3.5 ){
605  if ( hadEFrac >= 0 && hadEFrac < 0.157646) return 1;
606  else if( hadEFrac >= 0.157646 && hadEFrac < 0.290466) return 2;
607  else if( hadEFrac >= 0.290466 && hadEFrac < 0.468987) return 3;
608  else if( hadEFrac >= 0.468987 && hadEFrac <= 1 ) return 4;
609  }
610  else if( recoE >= 3.5 && recoE < 4 ){
611  if ( hadEFrac >= 0 && hadEFrac < 0.153759) return 1;
612  else if( hadEFrac >= 0.153759 && hadEFrac < 0.295047) return 2;
613  else if( hadEFrac >= 0.295047 && hadEFrac < 0.495104) return 3;
614  else if( hadEFrac >= 0.495104 && hadEFrac <= 1 ) return 4;
615  }
616  else if( recoE >= 4 && recoE < 5 ){
617  if ( hadEFrac >= 0 && hadEFrac < 0.154679) return 1;
618  else if( hadEFrac >= 0.154679 && hadEFrac < 0.299486) return 2;
619  else if( hadEFrac >= 0.299486 && hadEFrac < 0.542474) return 3;
620  else if( hadEFrac >= 0.542474 && hadEFrac <= 1 ) return 4;
621  }
622 
623  LOG_DEBUG("NuMuAnalysisSetup")
624  << "cannot compute quantile: total energy outside range";
625  return -1;
626 
627  }// quantiles
628 
629  //----------------------------------------------------------------------------
630 
632  double & hadE,
633  fnex::BeamType_t const& beamType ) const
634  {
635 
636  if (beamType == fnex::kFHC) return FourthAnalysisQuantileFHC(lepE, hadE);
637  else if(beamType == fnex::kRHC) return FourthAnalysisQuantileRHC(lepE, hadE);
638 
639  LOG_DEBUG("NuMuAnalysisSetup")
640  << "uknown BeamType revert to Third Analysis Numu quantiles";
641  return ThirdAnalysisQuantile(lepE, hadE);
642 
643  } // fourth analysis quantiles
644 
645  //----------------------------------------------------------------------------
646 
648  double & hadE ) const
649  {
650 
651  double recoE = lepE + hadE;
652  if( recoE == 0 ) {
653  LOG_DEBUG("NuMuAnalysisSetup")
654  << "cannot compute quantile: total energy is zero";
655  return -1;
656  }
657 
658  double hadEFrac = hadE / ( lepE + hadE );
659 
660  if( recoE >= 0 && recoE < 0.75 ){
661  if ( hadEFrac >= 0 && hadEFrac < 0.130486) return 1;
662  else if( hadEFrac >= 0.130486 && hadEFrac < 0.182243) return 2;
663  else if( hadEFrac >= 0.182243 && hadEFrac < 0.233157) return 3;
664  else if( hadEFrac >= 0.233157 && hadEFrac < 1.0) return 4;
665  }
666  else if( recoE >= 0.75 && recoE < 1 ){
667  if ( hadEFrac >= 0 && hadEFrac < 0.161791) return 1;
668  else if( hadEFrac >= 0.161791 && hadEFrac < 0.234801) return 2;
669  else if( hadEFrac >= 0.234801 && hadEFrac < 0.319627) return 3;
670  else if( hadEFrac >= 0.319627 && hadEFrac < 1.0) return 4;
671  }
672  else if( recoE >= 1 && recoE < 1.1 ){
673  if ( hadEFrac >= 0 && hadEFrac < 0.166412) return 1;
674  else if( hadEFrac >= 0.166412 && hadEFrac < 0.262069) return 2;
675  else if( hadEFrac >= 0.262069 && hadEFrac < 0.358449) return 3;
676  else if( hadEFrac >= 0.358449 && hadEFrac < 1.0) return 4;
677  }
678  else if( recoE >= 1.1 && recoE < 1.2 ){
679  if ( hadEFrac >= 0 && hadEFrac < 0.169974) return 1;
680  else if( hadEFrac >= 0.169974 && hadEFrac < 0.272577) return 2;
681  else if( hadEFrac >= 0.272577 && hadEFrac < 0.379373) return 3;
682  else if( hadEFrac >= 0.379373 && hadEFrac < 1.0) return 4;
683  }
684  else if( recoE >= 1.2 && recoE < 1.3 ){
685  if ( hadEFrac >= 0 && hadEFrac < 0.169908) return 1;
686  else if( hadEFrac >= 0.169908 && hadEFrac < 0.279256) return 2;
687  else if( hadEFrac >= 0.279256 && hadEFrac < 0.394825) return 3;
688  else if( hadEFrac >= 0.394825 && hadEFrac < 1.0) return 4;
689  }
690  else if( recoE >= 1.3 && recoE < 1.4 ){
691  if ( hadEFrac >= 0 && hadEFrac < 0.163625) return 1;
692  else if( hadEFrac >= 0.163625 && hadEFrac < 0.280182) return 2;
693  else if( hadEFrac >= 0.280182 && hadEFrac < 0.406361) return 3;
694  else if( hadEFrac >= 0.406361 && hadEFrac < 1.0) return 4;
695  }
696  else if( recoE >= 1.4 && recoE < 1.5 ){
697  if ( hadEFrac >= 0 && hadEFrac < 0.162347) return 1;
698  else if( hadEFrac >= 0.162347 && hadEFrac < 0.278829) return 2;
699  else if( hadEFrac >= 0.278829 && hadEFrac < 0.412204) return 3;
700  else if( hadEFrac >= 0.412204 && hadEFrac < 1.0) return 4;
701  }
702  else if( recoE >= 1.5 && recoE < 1.6 ){
703  if ( hadEFrac >= 0 && hadEFrac < 0.156687) return 1;
704  else if( hadEFrac >= 0.156687 && hadEFrac < 0.276816) return 2;
705  else if( hadEFrac >= 0.276816 && hadEFrac < 0.416131) return 3;
706  else if( hadEFrac >= 0.416131 && hadEFrac < 1.0) return 4;
707  }
708  else if( recoE >= 1.6 && recoE < 1.7 ){
709  if ( hadEFrac >= 0 && hadEFrac < 0.149375) return 1;
710  else if( hadEFrac >= 0.149375 && hadEFrac < 0.267368) return 2;
711  else if( hadEFrac >= 0.267368 && hadEFrac < 0.415353) return 3;
712  else if( hadEFrac >= 0.415353 && hadEFrac < 1.0) return 4;
713  }
714  else if( recoE >= 1.7 && recoE < 1.8 ){
715  if ( hadEFrac >= 0 && hadEFrac < 0.141755) return 1;
716  else if( hadEFrac >= 0.141755 && hadEFrac < 0.260274) return 2;
717  else if( hadEFrac >= 0.260274 && hadEFrac < 0.416843) return 3;
718  else if( hadEFrac >= 0.416843 && hadEFrac < 1.0) return 4;
719  }
720  else if( recoE >= 1.8 && recoE < 1.9 ){
721  if ( hadEFrac >= 0 && hadEFrac < 0.139631) return 1;
722  else if( hadEFrac >= 0.139631 && hadEFrac < 0.254985) return 2;
723  else if( hadEFrac >= 0.254985 && hadEFrac < 0.416642) return 3;
724  else if( hadEFrac >= 0.416642 && hadEFrac < 1.0) return 4;
725  }
726  else if( recoE >= 1.9 && recoE < 2 ){
727  if ( hadEFrac >= 0 && hadEFrac < 0.137587) return 1;
728  else if( hadEFrac >= 0.137587 && hadEFrac < 0.252339) return 2;
729  else if( hadEFrac >= 0.252339 && hadEFrac < 0.417286) return 3;
730  else if( hadEFrac >= 0.417286 && hadEFrac < 1.0) return 4;
731  }
732  else if( recoE >= 2 && recoE < 2.25 ){
733  if ( hadEFrac >= 0 && hadEFrac < 0.141741) return 1;
734  else if( hadEFrac >= 0.141741 && hadEFrac < 0.25798) return 2;
735  else if( hadEFrac >= 0.25798 && hadEFrac < 0.428896) return 3;
736  else if( hadEFrac >= 0.428896 && hadEFrac < 1.0) return 4;
737  }
738  else if( recoE >= 2.25 && recoE < 2.5 ){
739  if ( hadEFrac >= 0 && hadEFrac < 0.14987) return 1;
740  else if( hadEFrac >= 0.14987 && hadEFrac < 0.271534) return 2;
741  else if( hadEFrac >= 0.271534 && hadEFrac < 0.445448) return 3;
742  else if( hadEFrac >= 0.445448 && hadEFrac < 1.0) return 4;
743  }
744  else if( recoE >= 2.5 && recoE < 2.75 ){
745  if ( hadEFrac >= 0 && hadEFrac < 0.148543) return 1;
746  else if( hadEFrac >= 0.148543 && hadEFrac < 0.27805) return 2;
747  else if( hadEFrac >= 0.27805 && hadEFrac < 0.453015) return 3;
748  else if( hadEFrac >= 0.453015 && hadEFrac < 1.0) return 4;
749  }
750  else if( recoE >= 2.75 && recoE < 3 ){
751  if ( hadEFrac >= 0 && hadEFrac < 0.14472) return 1;
752  else if( hadEFrac >= 0.14472 && hadEFrac < 0.279848) return 2;
753  else if( hadEFrac >= 0.279848 && hadEFrac < 0.460109) return 3;
754  else if( hadEFrac >= 0.460109 && hadEFrac < 1.0) return 4;
755  }
756  else if( recoE >= 3 && recoE < 3.5 ){
757  if ( hadEFrac >= 0 && hadEFrac < 0.142453) return 1;
758  else if( hadEFrac >= 0.142453 && hadEFrac < 0.283427) return 2;
759  else if( hadEFrac >= 0.283427 && hadEFrac < 0.46781) return 3;
760  else if( hadEFrac >= 0.46781 && hadEFrac < 1.0) return 4;
761  }
762  else if( recoE >= 3.5 && recoE < 4 ){
763  if ( hadEFrac >= 0 && hadEFrac < 0.141182) return 1;
764  else if( hadEFrac >= 0.141182 && hadEFrac < 0.286712) return 2;
765  else if( hadEFrac >= 0.286712 && hadEFrac < 0.488556) return 3;
766  else if( hadEFrac >= 0.488556 && hadEFrac < 1.0) return 4;
767  }
768  else if( recoE >= 4 && recoE < 5 ){
769  if ( hadEFrac >= 0 && hadEFrac < 0.14111) return 1;
770  else if( hadEFrac >= 0.14111 && hadEFrac < 0.28538) return 2;
771  else if( hadEFrac >= 0.28538 && hadEFrac < 0.513618) return 3;
772  else if( hadEFrac >= 0.513618 && hadEFrac < 1.0) return 4;
773  }
774 
775  LOG_DEBUG("NuMuAnalysisSetup")
776  << "cannot compute quantile: total energy outside range";
777  return -1;
778 
779  }// fourth analysis quantiles FHC
780 
781  //----------------------------------------------------------------------------
782 
784  double & hadE ) const
785  {
786 
787  double recoE = lepE + hadE;
788  if( recoE == 0 ) {
789  LOG_DEBUG("NuMuAnalysisSetup")
790  << "cannot compute quantile: total energy is zero";
791  return -1;
792  }
793 
794  double hadEFrac = hadE / ( lepE + hadE );
795 
796  LOG_VERBATIM("NuMuAnalysisSetup") << "FourthAnalysisQuantileRHC ";
797 
798  if( recoE >= 0 && recoE < 0.75 ){
799  if ( hadEFrac >= 0 && hadEFrac < 0.0980177) return 1;
800  else if( hadEFrac >= 0.0980177 && hadEFrac < 0.140308) return 2;
801  else if( hadEFrac >= 0.140308 && hadEFrac < 0.201423) return 3;
802  else if( hadEFrac >= 0.201423 && hadEFrac < 1.0) return 4;
803  }
804  else if( recoE >= 0.75 && recoE < 1 ){
805  if ( hadEFrac >= 0 && hadEFrac < 0.123805) return 1;
806  else if( hadEFrac >= 0.123805 && hadEFrac < 0.200619) return 2;
807  else if( hadEFrac >= 0.200619 && hadEFrac < 0.293348) return 3;
808  else if( hadEFrac >= 0.293348 && hadEFrac < 1.0) return 4;
809  }
810  else if( recoE >= 1 && recoE < 1.1 ){
811  if ( hadEFrac >= 0 && hadEFrac < 0.129862) return 1;
812  else if( hadEFrac >= 0.129862 && hadEFrac < 0.220222) return 2;
813  else if( hadEFrac >= 0.220222 && hadEFrac < 0.322866) return 3;
814  else if( hadEFrac >= 0.322866 && hadEFrac < 1.0) return 4;
815  }
816  else if( recoE >= 1.1 && recoE < 1.2 ){
817  if ( hadEFrac >= 0 && hadEFrac < 0.130338) return 1;
818  else if( hadEFrac >= 0.130338 && hadEFrac < 0.225444) return 2;
819  else if( hadEFrac >= 0.225444 && hadEFrac < 0.331177) return 3;
820  else if( hadEFrac >= 0.331177 && hadEFrac < 1.0) return 4;
821  }
822  else if( recoE >= 1.2 && recoE < 1.3 ){
823  if ( hadEFrac >= 0 && hadEFrac < 0.126173) return 1;
824  else if( hadEFrac >= 0.126173 && hadEFrac < 0.228904) return 2;
825  else if( hadEFrac >= 0.228904 && hadEFrac < 0.339258) return 3;
826  else if( hadEFrac >= 0.339258 && hadEFrac < 1.0) return 4;
827  }
828  else if( recoE >= 1.3 && recoE < 1.4 ){
829  if ( hadEFrac >= 0 && hadEFrac < 0.116904) return 1;
830  else if( hadEFrac >= 0.116904 && hadEFrac < 0.215843) return 2;
831  else if( hadEFrac >= 0.215843 && hadEFrac < 0.339619) return 3;
832  else if( hadEFrac >= 0.339619 && hadEFrac < 1.0) return 4;
833  }
834  else if( recoE >= 1.4 && recoE < 1.5 ){
835  if ( hadEFrac >= 0 && hadEFrac < 0.112387) return 1;
836  else if( hadEFrac >= 0.112387 && hadEFrac < 0.212124) return 2;
837  else if( hadEFrac >= 0.212124 && hadEFrac < 0.342516) return 3;
838  else if( hadEFrac >= 0.342516 && hadEFrac < 1.0) return 4;
839  }
840  else if( recoE >= 1.5 && recoE < 1.6 ){
841  if ( hadEFrac >= 0 && hadEFrac < 0.108383) return 1;
842  else if( hadEFrac >= 0.108383 && hadEFrac < 0.203952) return 2;
843  else if( hadEFrac >= 0.203952 && hadEFrac < 0.340037) return 3;
844  else if( hadEFrac >= 0.340037 && hadEFrac < 1.0) return 4;
845  }
846  else if( recoE >= 1.6 && recoE < 1.7 ){
847  if ( hadEFrac >= 0 && hadEFrac < 0.102762) return 1;
848  else if( hadEFrac >= 0.102762 && hadEFrac < 0.196337) return 2;
849  else if( hadEFrac >= 0.196337 && hadEFrac < 0.337452) return 3;
850  else if( hadEFrac >= 0.337452 && hadEFrac < 1.0) return 4;
851  }
852  else if( recoE >= 1.7 && recoE < 1.8 ){
853  if ( hadEFrac >= 0 && hadEFrac < 0.099041) return 1;
854  else if( hadEFrac >= 0.099041 && hadEFrac < 0.188072) return 2;
855  else if( hadEFrac >= 0.188072 && hadEFrac < 0.330353) return 3;
856  else if( hadEFrac >= 0.330353 && hadEFrac < 1.0) return 4;
857  }
858  else if( recoE >= 1.8 && recoE < 1.9 ){
859  if ( hadEFrac >= 0 && hadEFrac < 0.0965183) return 1;
860  else if( hadEFrac >= 0.0965183 && hadEFrac < 0.184834) return 2;
861  else if( hadEFrac >= 0.184834 && hadEFrac < 0.329758) return 3;
862  else if( hadEFrac >= 0.329758 && hadEFrac < 1.0) return 4;
863  }
864  else if( recoE >= 1.9 && recoE < 2 ){
865  if ( hadEFrac >= 0 && hadEFrac < 0.0938676) return 1;
866  else if( hadEFrac >= 0.0938676 && hadEFrac < 0.180883) return 2;
867  else if( hadEFrac >= 0.180883 && hadEFrac < 0.325615) return 3;
868  else if( hadEFrac >= 0.325615 && hadEFrac < 1.0) return 4;
869  }
870  else if( recoE >= 2 && recoE < 2.25 ){
871  if ( hadEFrac >= 0 && hadEFrac < 0.095024) return 1;
872  else if( hadEFrac >= 0.095024 && hadEFrac < 0.183673) return 2;
873  else if( hadEFrac >= 0.183673 && hadEFrac < 0.332639) return 3;
874  else if( hadEFrac >= 0.332639 && hadEFrac < 1.0) return 4;
875  }
876  else if( recoE >= 2.25 && recoE < 2.5 ){
877  if ( hadEFrac >= 0 && hadEFrac < 0.100057) return 1;
878  else if( hadEFrac >= 0.100057 && hadEFrac < 0.191958) return 2;
879  else if( hadEFrac >= 0.191958 && hadEFrac < 0.349162) return 3;
880  else if( hadEFrac >= 0.349162 && hadEFrac < 1.0) return 4;
881  }
882  else if( recoE >= 2.5 && recoE < 2.75 ){
883  if ( hadEFrac >= 0 && hadEFrac < 0.10519) return 1;
884  else if( hadEFrac >= 0.10519 && hadEFrac < 0.203698) return 2;
885  else if( hadEFrac >= 0.203698 && hadEFrac < 0.366594) return 3;
886  else if( hadEFrac >= 0.366594 && hadEFrac < 1.0) return 4;
887  }
888  else if( recoE >= 2.75 && recoE < 3 ){
889  if ( hadEFrac >= 0 && hadEFrac < 0.108997) return 1;
890  else if( hadEFrac >= 0.108997 && hadEFrac < 0.214948) return 2;
891  else if( hadEFrac >= 0.214948 && hadEFrac < 0.371937) return 3;
892  else if( hadEFrac >= 0.371937 && hadEFrac < 1.0) return 4;
893  }
894  else if( recoE >= 3 && recoE < 3.5 ){
895  if ( hadEFrac >= 0 && hadEFrac < 0.118638) return 1;
896  else if( hadEFrac >= 0.118638 && hadEFrac < 0.232254) return 2;
897  else if( hadEFrac >= 0.232254 && hadEFrac < 0.389001) return 3;
898  else if( hadEFrac >= 0.389001 && hadEFrac < 1.0) return 4;
899  }
900  else if( recoE >= 3.5 && recoE < 4 ){
901  if ( hadEFrac >= 0 && hadEFrac < 0.122369) return 1;
902  else if( hadEFrac >= 0.122369 && hadEFrac < 0.243635) return 2;
903  else if( hadEFrac >= 0.243635 && hadEFrac < 0.422904) return 3;
904  else if( hadEFrac >= 0.422904 && hadEFrac < 1.0) return 4;
905  }
906  else if( recoE >= 4 && recoE < 5 ){
907  if ( hadEFrac >= 0 && hadEFrac < 0.128739) return 1;
908  else if( hadEFrac >= 0.128739 && hadEFrac < 0.248248) return 2;
909  else if( hadEFrac >= 0.248248 && hadEFrac < 0.455024) return 3;
910  else if( hadEFrac >= 0.455024 && hadEFrac < 1.0) return 4;
911  }
912 
913  LOG_DEBUG("NuMuAnalysisSetup")
914  << "cannot compute quantile: total energy outside range";
915  return -1;
916 
917 
918  }// fourth analysis quantiles RHC
919 
920  //----------------------------------------------------------------------------
921 
923  rb::Track const& track,
924  double & hadE,
925  double & muonE,
926  bool isMC,
927  fnex::BeamType_t const& beamType ) const
928  {
929  auto det = fDetService->DetId();
930  auto run = fRH->RunNumber();
931 
932  double hadVisE = numue.HadTrkE() + numue.HadCalE();
933 
934  bool isRHC = false; // FHC mode
935  if(beamType == fnex::kRHC) isRHC = true ; // RHC mode
936  LOG_DEBUG("NuMuAnalysisSetup")
937  << "Inside NuMuAnalysisSetup::Analysis2018Energy function. ";
938  //if(isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "RHC file";
939  //else LOG_VERBATIM("NuMuAnalysisSetup") << "FHC file";
940 
942  {
943 
944  double trkLen = track.TotalLength();
945 
946  //if(isMC && !isRHC) LOG_DEBUG("NuMuAnalysisSetup") << "MC, not RHC";
947  //if(!isMC && !isRHC) LOG_DEBUG("NuMuAnalysisSetup") << "Not MC, not RHC";
948 
949  muonE = NumuEnergyFunc::predict_prod4_fd_muon_energy(trkLen, run, isMC, isRHC);
950  hadE = NumuEnergyFunc::predict_prod4_fd_had_energy(hadVisE, run, isMC, isRHC);
951 
952  LOG_DEBUG("NuMuAnalysisSetup")
953  << "Inside NuMuAnalysisSetup::Analysis2018Energy function. "
954  << "LepE = " << muonE << "\t HadE = " << hadE;
955  }
956 
958  {
959 
960  double trkLenAct = numue.NDTrkLenAct();
961  double trkLenCat = numue.NDTrkLenCat();
962 
963  // if(isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "MC, not RHC";
964  // if(!isMC && !isRHC) LOG_VERBATIM("NuMuAnalysisSetup") << "Not MC, not RHC";
965 
967  hadE = NumuEnergyFunc::predict_prod4_nd_had_energy(hadVisE, isRHC);
968  }
969 
970  }// end Ana2018 energy estimator
971 
972  //----------------------------------------------------------------------------
973  unsigned int NuMuAnalysisSetup::HighestPIDTrack(std::vector< art::Ptr<rb::Track> > const& sliceTracks,
974  art::InputTag const& remidTag,
975  art::Event const& e) const
976  {
977  int bestTrack = -1;
978  double highestPID = -2.0;
979  double bestTrackLength = 0.0;
980  double bestStartZ = std::numeric_limits<double>::max();
981  double bestADC = std::numeric_limits<double>::max();
982 
983  art::FindOneP<remid::ReMId> tracksToReMId(sliceTracks, e, remidTag);
984 
985  for(size_t c = 0; c < sliceTracks.size(); ++c){
986 
987  //geo::View_t view = sliceTracks[c]->View();
988  //Only use 3D tracks
989  if( sliceTracks[c]->NXCell() == 0 || sliceTracks[c]->NYCell() == 0 ){
990  continue;
991  } //End of loop over 2d tracks
992 
993  if(tracksToReMId.isValid()){
994  art::Ptr<remid::ReMId> trackReMId = tracksToReMId.at(c);
995 
996  // If the PID is lower, we can bail, nothing to update
997  if(trackReMId->Value() == highestPID)
998  { // If the new guy is shorter, we can bail
999 
1000  if (sliceTracks[c]->TotalLength() < bestTrackLength) continue;
1001 
1002  // If equal, need to fall to the next level, low start Z
1003  if (sliceTracks[c]->TotalLength() == bestTrackLength){
1004 
1005  // If start Z is higher, we can bail
1006  if (sliceTracks[c]->Start().Z() > bestStartZ) continue;
1007 
1008  // If equal, need to fall to the next level, total ADC
1009  // You might think this is stupid, but it's equivalent to flipping
1010  // a coin which always lands on the same side. This is good, this
1011  // will happen so rarely that you should just pee up a rope.
1012  if (sliceTracks[c]->Start().Z() == bestStartZ){
1013 
1014  if(sliceTracks[c]->TotalADC() > bestADC) continue;
1015 
1016  if(sliceTracks[c]->TotalADC() == bestADC)
1017  throw cet::exception("ReMId Matched Track Exception")<<std::endl
1018  << "Can't break a tie between two ReMId/tracks. "
1019  << "This should have never happend, but it did. "
1020  << "The numu group is going to have to reflect on this "
1021  << "and find a way to fix it." << std::endl;
1022  }
1023 
1024 
1025  }
1026  }
1027  highestPID = trackReMId->Value();
1028  bestTrackLength = sliceTracks[c]->TotalLength();
1029  bestStartZ = sliceTracks[c]->Start().Z();
1030  bestADC = sliceTracks[c]->TotalADC();
1031  bestTrack = c;
1032 
1033  }//End of valid ReMId object
1034  }//End of for loop over tracks
1035 
1036  if (bestTrack < 0)
1037  return 999;
1038 
1039  return (unsigned int)bestTrack;
1040 
1041  } // highestpidtrack
1042 
1043 } // namespace fnex
1044 
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
novadaq::cnv::DetId DetId() const
What detector are we in?
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
SubRunNumber_t subRun() const
Definition: Event.h:72
float NDTrkLenCat() const
Definition: NumuE.cxx:309
Energy estimators for CC events.
Definition: FillEnergies.h:7
int fHadEFracQuantile
Hadronic Energy Fraction quantile.
float E() const
Definition: Energy.cxx:27
void FillVars(art::Event const &evt, art::ValidHandle< std::vector< rb::Cluster > > const &slices, art::InputTag const &tag, std::vector< fnex::DataVarVals > &dataVarVals, std::vector< fnex::MCVarVals > &mcVarVals, std::vector< fnex::EventId > &evIds, fnex::BeamType_t const &beamType) const override
void ThirdAnalysisEnergy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:930
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
float NDTrkLenAct() const
Definition: NumuE.cxx:304
static const unsigned char kRecoQ2
Definition: VarVals.h:38
Definition: event.h:19
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
float TrkCCE() const
Definition: NumuE.cxx:254
double predict_prod3_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
Create a list of fnex::Events to be used in fits.
bool isRealData() const
Definition: Event.h:83
bool isValid() const
Definition: maybe_ref.h:60
float NDTrkCalAct() const
Definition: NumuE.cxx:314
double Value() const
Definition: PID.h:22
object containing MC flux information
void Analysis2018Energy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE, bool isMC, fnex::BeamType_t const &beamType) const
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
Far Detector at Ash River, MN.
int FourthAnalysisQuantileRHC(double &lepE, double &hadE) const
art::ServiceHandle< nova::dbi::RunHistoryService > fRH
RunHistory service.
int FourthAnalysisQuantileFHC(double &lepE, double &hadE) const
void Reconfigure(fhicl::ParameterSet const &pset) override
virtual void Reconfigure(fhicl::ParameterSet const &pset)
float HadTrkE() const
Definition: NumuE.cxx:284
static const unsigned char kLep_RecoE_MCFrac
Definition: VarVals.h:37
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
void FillRecoVars(fnex::DataVarVals &vars, numue::NumuE const &, rb::Energy const &, rb::Track const &, bool isMC, fnex::BeamType_t const &beamType) const
double energy
Definition: plottest35.C:25
Near Detector in the NuMI cavern.
double predict_prod4_nd_cat_energy(double ndtrklencat, bool isRHC)
float val_at(unsigned char const &varkey, fnex::MetaData const &md) const
Definition: VarVals.cxx:499
const std::map< std::pair< std::string, std::string >, Variable > vars
EventNumber_t event() const
Definition: Event.h:67
double predict_prod4_nd_act_energy(double ndtrklenact, bool isRHC)
int FourthAnalysisQuantile(double &lepE, double &hadE, fnex::BeamType_t const &beamType) const
unsigned int HighestPIDTrack(std::vector< art::Ptr< rb::Track > > const &sliceTracks, art::InputTag const &remidTag, art::Event const &e) const
int ThirdAnalysisQuantile(double &recoE, double &hadEFrac) const
Definition: run.py:1
#define LOG_WARNING(category)
double predict_prod4_fd_muon_energy(double trkLen, int run, bool ismc, bool isRHC)
Definition: Numu2018Fits.cxx:7
art::ServiceHandle< ds::DetectorService > fDetService
which detector?
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
float HadCalE() const
Definition: NumuE.cxx:279
fnex::SelectionType_t selectionType
Definition: Structs.h:52
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double predict_prod3_nd_cat_energy(double ndtrklencat, bool isRHC)
double predict_prod3_nd_had_energy(double hadVisE, bool isRHC)
A container for energy information.
Definition: Energy.h:20
double predict_prod4_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
enum fnex::beam_type BeamType_t
void set_val_at(std::string const &varkey, float const &val)
Definition: VarVals.cxx:525
double predict_prod4_nd_had_energy(double hadVisE, bool isRHC)
static const unsigned char kHad_RecoE
Definition: VarVals.h:34
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:469
bool FillTruthVars(fnex::MCVarVals &vars, art::Event const &evt, size_t const &sliceNum, art::FindOneP< simb::MCTruth > const &fomct, art::FindOne< simb::MCFlux > const &fomcf, art::FindOne< fxwgt::FluxWeights > const &fofw, art::FindOne< simb::GTruth > const &fogt) const
NuMuAnalysisSetup(fhicl::ParameterSet const &pset)
#define LOG_VERBATIM(category)
float NDTrkCalTran() const
Definition: NumuE.cxx:319
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static const unsigned char kLep_RecoE
Definition: VarVals.h:35
int RunNumber() const
Definition: RunHistory.h:374
Base container for the MC related Vars that constitute an event.
Definition: VarVals.h:874
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
void SecondAnalysisEnergy(numue::NumuE const &numue, rb::Track const &track, double &hadE, double &muonE) const
double predict_prod3_fd_had_energy(double hadVisE, int run, bool ismc, bool isRHC)
float NDTrkCalCat() const
Definition: NumuE.cxx:324
double predict_prod3_nd_act_energy(double ndtrklenact, bool isRHC)