BPFTmvaTrainer_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file BPFTmvaTrainer_module.cc
3 // \brief Module to create the TTree(s) used to train TMVA for the
4 // BPF particle PIDs and the BPF energy estimator. This module
5 // is ONLY intended to run over MC files.
6 // \version
7 // \author $Author: mbaird42 $
8 ////////////////////////////////////////////////////////////////////////
9 
10 // C/C++ includes
11 #include <cmath>
12 #include <iostream>
13 #include <vector>
14 
15 // ROOT includes
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TFile.h"
19 #include "TTree.h"
20 #include "TMVA/Factory.h"
21 #include "TMVA/Tools.h"
22 #include "TMVA/Reader.h"
23 #include "TMVA/MethodKNN.h"
24 
25 // Framework includes
36 #include "fhiclcpp/ParameterSet.h"
38 
39 // NOvASoft includes
40 #include "RecoBase/CellHit.h"
41 #include "RecoBase/RecoHit.h"
42 #include "RecoBase/Cluster.h"
43 #include "RecoBase/Prong.h"
44 #include "RecoBase/Track.h"
45 #include "RecoBase/FitSum.h"
46 
47 #include "ReMId/ReMId.h"
48 
49 #include "Utilities/AssociationUtil.h"
50 #include "Geometry/Geometry.h"
51 #include "Calibrator/Calibrator.h"
52 
55 
56 #include "MCCheater/BackTracker.h"
59 
60 #include "SummaryData/POTSum.h"
61 
62 
63 
64 // BPFTmvaTrainer header
65 namespace bpfit {
67  public:
68  explicit BPFTmvaTrainer(fhicl::ParameterSet const& pset);
69  ~BPFTmvaTrainer();
70 
71  void analyze(const art::Event& evt);
72  void reconfigure(const fhicl::ParameterSet& pset);
73  void beginJob();
74  void endJob();
75  virtual void endSubRun(const art::SubRun& sr);
76 
77  void resetVars();
78  void setTruthVars(const std::vector<cheat::NeutrinoEffPur>& nuEP);
79 
80  private:
81  // FCL parameters:
84  std::string fProngLabel; ///< label for module that made the prongs
85  std::string fProngInstance; ///< instance label for the prongs to be used
86  std::string fTrackLabel; ///< label for module that made the tracks
87  std::string fHistoFile; ///< location of dE/dx log-likelihood histogram
88  std::string fXMLFile; ///< location of the XML file for the muon PID
89  double fdxTOL; ///< lower limit for dx values to be used in dE/dx calculation
90  unsigned int fTrainMode; ///< training mode (energy estimator or muon PID?)
91  int fPIDmode; ///< make the training TTree using the muon, pion, or proton fit
92 
93  // Things used for internal processing and helper classes:
94  TH2F *fdEdxVSbg[2][2]; ///< histo of computed dE/dx for cell hits (vs. beta*gamma)
95  bpfit::dEdxCalculator fdEdxCalc; ///< helper class for computing dEdx values
96  TMVA::Reader fMuonPID; ///< TMVA reader object for the muon PID
97  murem::TrackCleanUpAlg fTrackCleanUpAlg; ///< object to determine hadronic energy on the muon track
98  double fTotalPOT; ///< total POT counted for all events in this job
99 
100  // Output TTrees:
101  TTree *fBPFTrackSignalTree; ///< TTree with track level signal info (for PID training)
102  TTree *fBPFTrackBackgroundTree; ///< TTree with track level background info (for PID training)
103  TTree *fBPFSliceTree; ///< TTree with slice level info (for energy estimator training)
104 
105 
106 
107  //
108  // Variables to go into the Track TTrees.
109  //
110 
111  // training variables:
112  float fLength; ///< track length
113  float fChi2H; ///< chi^2 hits value per DOF
114  float fChi2A; ///< chi^2 angle value per DOF
115  float fChi2T; ///< total chi^2
116  float fdEdxLL; ///< dEdx LL value
117  float fHitRatio; ///< ratio of track hits to prong hits (used as a track/shower discriminator)
118  float fLengthRatio; ///< ratio of track to longest track in the slice
119  float fDChi2Tmupi; ///< delta chi2 total muon - pion
120  float fDChi2Tmupr; ///< delta chi2 total muon - proton
121  float fDChi2Tpipr; ///< delta chi2 total pion - proton
122 
123  // track PDG code used to separate signal and background
124  float fTrackPDG; ///< PDG code for the primary track particle
125 
126 
127 
128  //
129  // Variables to go into the Slice TTree. Input variables are mapped to the
130  // target variables using the cut variables for event selection.
131  //
132  // fE and the cut variables are also used for the PID tree
133  //
134 
135  // target variables:
136  float fE; ///< true energy
137  float fEsq; ///< true energy squared (used to estiamte st. dev.)
138 
139  // input variables:
140  float fPmu; ///< reconstructed muon momentum from the track selected to be the muon
141  float fDirZmu; ///< muon track dirZ from the track selected to be the muon
142  float fN3Dprongs; ///< total number of 3D prongs
143  float fE3Dprongs; ///< sum of hit energies in 3D prongs (excluding the "muon" prong) - this does NOT account for dead material
144  float fEremain; ///< sum of remaining slice energy not in 3D prongs PLUS extra energy added back in from TrackCleanUpAlg
145  float fSumPE; ///< sum of the PE for all hits not on the muon track
146 
147  // spectator variables:
148  float fEhadMu; ///< Hadronic energy contamination on the muon track as computed by TrackCleanUpAlg
149 
150  // cut variables:
151  float fMinX; ///< minimum hit X value in the slice
152  float fMinY; ///< minimum hit Y value in the slice
153  float fMinZ; ///< minimum hit Z value in the slice
154  float fMaxX; ///< maximum hit X value in the slice
155  float fMaxY; ///< maximum hit Y value in the slice
156  float fMaxZ; ///< maximum hit Z value in the slice
157  int fNHitX; ///< number of X-view hits in the slice
158  int fNHitY; ///< number of Y-view hits in the slice
159  int fCCNC; ///< is this CC or NC (from truth)
160  int fIntType; ///< neutrino interaction type (from truth)
161  int fnuPDG; ///< PDG code for the best matched nu (from truth)
162  float fBestBPFmuPID; ///< best BPF muon PID value
163  float fBestRemid; ///< best remid value
164  int fBestBPFmuPIDpdg; ///< pdg code for the particle that contributed the most energy to the track with the best BPF muon PID
165 
166  TH1D *fPOT; ///< Total POT accumulated for all events processed in this job.
167 
168  };
169 }
170 
171 
172 
173 namespace bpfit
174 {
175  //.......................................................................
177  : EDAnalyzer(pset),
178  fSlicerToken(consumes<std::vector<rb::Cluster>>(pset.get<std::string>("SlicerLabel"))),
179  fPOTSumToken(consumes<sumdata::POTSum, art::InSubRun>("generator")),
180  fdEdxCalc(),
181  fMuonPID(),
182  fTrackCleanUpAlg(pset.get<fhicl::ParameterSet>("TrackCleanUpAlgPSet")),
183  fTotalPOT(0.0)
184  {
185  this->reconfigure(pset);
186  }
187 
188  //......................................................................
190  {
191  //======================================================================
192  // Clean up any memory allocated by your module... ...if you dare...
193  //======================================================================
194  }
195 
196  //......................................................................
198  {
199  fProngLabel = pset.get<std::string> ("ProngLabel");
200  fProngInstance = pset.get<std::string> ("ProngInstance");
201  fTrackLabel = pset.get<std::string> ("TrackLabel");
202  fHistoFile = pset.get<std::string> ("HistoFile");
203  fXMLFile = pset.get<std::string> ("XMLFile");
204  fdxTOL = pset.get<double> ("dxTOL");
205  fTrainMode = pset.get<unsigned int>("TrainMode");
206  fPIDmode = pset.get<int> ("PIDmode");
207  }
208 
209  //......................................................................
211  {
212 
213  //
214  // Check that the requested PID training mode is supported.
215  //
216  if(fTrainMode == 0) {
217  if(fPIDmode != 13 && fPIDmode != 211 && fPIDmode != 2212) {
218  std::cerr << "\n\n\nWARNING: PID training mode for PDG code = "
219  << fPIDmode
220  << " is not currently supported. ABORTING...\n\n\n";
221  abort();
222  }
223  }
224 
225  //
226  // Read in the log-likelihood histos from the input file
227  //
228  TFile infile(fHistoFile.c_str());
229 
230  fdEdxVSbg[0][0] = (TH2F*)infile.FindObjectAny("fdEdxVSbg_W0_X0");
231  fdEdxVSbg[0][1] = (TH2F*)infile.FindObjectAny("fdEdxVSbg_W0_X1");
232  fdEdxVSbg[1][0] = (TH2F*)infile.FindObjectAny("fdEdxVSbg_W1_X0");
233  fdEdxVSbg[1][1] = (TH2F*)infile.FindObjectAny("fdEdxVSbg_W1_X1");
234 
235  infile.Close();
236 
237 
238 
239  //
240  // Set up the TMVA reader for the muon PID.
241  //
242  if(fTrainMode == 1) {
243  fMuonPID.AddVariable("length",&fLength);
244  fMuonPID.AddVariable("chi2T",&fChi2T);
245  fMuonPID.AddVariable("dEdxLL",&fdEdxLL);
246  fMuonPID.AddVariable("hitRatio",&fHitRatio);
247  fMuonPID.BookMVA("KNN Method",fXMLFile);
248  }
249 
250 
251 
253 
254  //
255  // Book the PID Training TTrees
256  //
257 
258  // The signal TTree...
259  fBPFTrackSignalTree = tfs->make<TTree>("BPFTrackSignalTree","Tree for TMVA training with signal track level BPF info");
260 
261  // input variables:
262  fBPFTrackSignalTree->Branch("length",&fLength,"length/F");
263  fBPFTrackSignalTree->Branch("chi2H",&fChi2H,"chi2H/F");
264  fBPFTrackSignalTree->Branch("chi2A",&fChi2A,"chi2A/F");
265  fBPFTrackSignalTree->Branch("chi2T",&fChi2T,"chi2T/F");
266  fBPFTrackSignalTree->Branch("dEdxLL",&fdEdxLL,"dEdxLL/F");
267  fBPFTrackSignalTree->Branch("hitRatio",&fHitRatio,"hitRatio/F");
268  fBPFTrackSignalTree->Branch("lengthRatio",&fLengthRatio,"lengthRatio/F");
269  fBPFTrackSignalTree->Branch("DChi2Tmupi",&fDChi2Tmupi,"DChi2Tmupi/F");
270  fBPFTrackSignalTree->Branch("DChi2Tmupr",&fDChi2Tmupr,"DChi2Tmupr/F");
271  fBPFTrackSignalTree->Branch("DChi2Tpipr",&fDChi2Tpipr,"DChi2Tpipr/F");
272 
273  // signal/background discriminator
274  fBPFTrackSignalTree->Branch("trackPDG",&fTrackPDG,"trackPDG/F");
275 
276  // cut/weight variables:
277  fBPFTrackSignalTree->Branch("E",&fE,"E/F");
278  fBPFTrackSignalTree->Branch("MinX",&fMinX,"MinX/F");
279  fBPFTrackSignalTree->Branch("MinY",&fMinY,"MinY/F");
280  fBPFTrackSignalTree->Branch("MinZ",&fMinZ,"MinZ/F");
281  fBPFTrackSignalTree->Branch("MaxX",&fMaxX,"MaxX/F");
282  fBPFTrackSignalTree->Branch("MaxY",&fMaxY,"MaxY/F");
283  fBPFTrackSignalTree->Branch("MaxZ",&fMaxZ,"MaxZ/F");
284  fBPFTrackSignalTree->Branch("NHitX",&fNHitX,"NHitX/I");
285  fBPFTrackSignalTree->Branch("NHitY",&fNHitY,"NHitY/I");
286  fBPFTrackSignalTree->Branch("CCNC",&fCCNC,"CCNC/I");
287  fBPFTrackSignalTree->Branch("IntType",&fIntType,"IntType/I");
288  fBPFTrackSignalTree->Branch("nuPDG",&fnuPDG,"nuPDG/I");
289 
290 
291 
292  // The background TTree...
293  fBPFTrackBackgroundTree = tfs->make<TTree>("BPFTrackBackgroundTree","Tree for TMVA training with background track level BPF info");
294 
295  // input variables:
296  fBPFTrackBackgroundTree->Branch("length",&fLength,"length/F");
297  fBPFTrackBackgroundTree->Branch("chi2H",&fChi2H,"chi2H/F");
298  fBPFTrackBackgroundTree->Branch("chi2A",&fChi2A,"chi2A/F");
299  fBPFTrackBackgroundTree->Branch("chi2T",&fChi2T,"chi2T/F");
300  fBPFTrackBackgroundTree->Branch("dEdxLL",&fdEdxLL,"dEdxLL/F");
301  fBPFTrackBackgroundTree->Branch("hitRatio",&fHitRatio,"hitRatio/F");
302  fBPFTrackBackgroundTree->Branch("lengthRatio",&fLengthRatio,"lengthRatio/F");
303  fBPFTrackBackgroundTree->Branch("DChi2Tmupi",&fDChi2Tmupi,"DChi2Tmupi/F");
304  fBPFTrackBackgroundTree->Branch("DChi2Tmupr",&fDChi2Tmupr,"DChi2Tmupr/F");
305  fBPFTrackBackgroundTree->Branch("DChi2Tpipr",&fDChi2Tpipr,"DChi2Tpipr/F");
306 
307 
308  // signal/background discriminator
309  fBPFTrackBackgroundTree->Branch("trackPDG",&fTrackPDG,"trackPDG/F");
310 
311  // cut/weight variables:
312  fBPFTrackBackgroundTree->Branch("E",&fE,"E/F");
313  fBPFTrackBackgroundTree->Branch("MinX",&fMinX,"MinX/F");
314  fBPFTrackBackgroundTree->Branch("MinY",&fMinY,"MinY/F");
315  fBPFTrackBackgroundTree->Branch("MinZ",&fMinZ,"MinZ/F");
316  fBPFTrackBackgroundTree->Branch("MaxX",&fMaxX,"MaxX/F");
317  fBPFTrackBackgroundTree->Branch("MaxY",&fMaxY,"MaxY/F");
318  fBPFTrackBackgroundTree->Branch("MaxZ",&fMaxZ,"MaxZ/F");
319  fBPFTrackBackgroundTree->Branch("NHitX",&fNHitX,"NHitX/I");
320  fBPFTrackBackgroundTree->Branch("NHitY",&fNHitY,"NHitY/I");
321  fBPFTrackBackgroundTree->Branch("CCNC",&fCCNC,"CCNC/I");
322  fBPFTrackBackgroundTree->Branch("IntType",&fIntType,"IntType/I");
323  fBPFTrackBackgroundTree->Branch("nuPDG",&fnuPDG,"nuPDG/I");
324 
325 
326 
327  //
328  // Book the Energy Estimator Training TTree
329  //
330  fBPFSliceTree = tfs->make<TTree>("BPFSliceTree","Tree for TMVA training with slice level BPF info");
331 
332  // target variables:
333  fBPFSliceTree->Branch("E",&fE,"E/F");
334  fBPFSliceTree->Branch("Esq",&fEsq,"Esq/F");
335 
336  // input variables:
337  fBPFSliceTree->Branch("Pmu",&fPmu,"Pmu/F");
338  fBPFSliceTree->Branch("DirZmu",&fDirZmu,"DirZmu/F");
339  fBPFSliceTree->Branch("N3Dprongs",&fN3Dprongs,"N3Dprongs/F");
340  fBPFSliceTree->Branch("E3Dprongs",&fE3Dprongs,"E3Dprongs/F");
341  fBPFSliceTree->Branch("Eremain",&fEremain,"Eremain/F");
342  fBPFSliceTree->Branch("SumPE",&fSumPE,"SumPE/F");
343 
344  // spectator variables:
345  fBPFSliceTree->Branch("EhadMu",&fEhadMu,"EhadMu/F");
346 
347  // cut variables:
348  fBPFSliceTree->Branch("MinX",&fMinX,"MinX/F");
349  fBPFSliceTree->Branch("MinY",&fMinY,"MinY/F");
350  fBPFSliceTree->Branch("MinZ",&fMinZ,"MinZ/F");
351  fBPFSliceTree->Branch("MaxX",&fMaxX,"MaxX/F");
352  fBPFSliceTree->Branch("MaxY",&fMaxY,"MaxY/F");
353  fBPFSliceTree->Branch("MaxZ",&fMaxZ,"MaxZ/F");
354  fBPFSliceTree->Branch("NHitX",&fNHitX,"NHitX/I");
355  fBPFSliceTree->Branch("NHitY",&fNHitY,"NHitY/I");
356  fBPFSliceTree->Branch("CCNC",&fCCNC,"CCNC/I");
357  fBPFSliceTree->Branch("IntType",&fIntType,"IntType/I");
358  fBPFSliceTree->Branch("nuPDG",&fnuPDG,"nuPDG/I");
359  fBPFSliceTree->Branch("BestBPFmuPID",&fBestBPFmuPID,"BestBPFmuPID/F");
360  fBPFSliceTree->Branch("BestRemid",&fBestRemid,"BestRemid/F");
361  fBPFSliceTree->Branch("BestBPFmuPIDpdg",&fBestBPFmuPIDpdg,"BestBPFmuPIDpdg/I");
362 
363 
364 
365  //
366  // Book histograms...
367  //
368  fPOT = tfs->make<TH1D>("fPOT","Total POT",1,0,1);
369 
370  }
371 
372  //......................................................................
374  {
375  fPOT->Fill(0.5,fTotalPOT);
376  }
377 
378  //......................................................................
380  {
382  sr.getByToken(fPOTSumToken, pot);
383  if(!pot.failedToGet()) fTotalPOT += pot->totgoodpot;
384  }
385 
386  //......................................................................
388  {
389  //
390  // PID variables:
391  //
392  fLength = -5.0;
393  fChi2H = 1.0e9;
394  fChi2A = 1.0e9;
395  fChi2T = 1.0e9;
396  fdEdxLL = -1.0e9;
397  fHitRatio = -5.0;
398  fLengthRatio = -5.0;
399  fDChi2Tmupi = 1.0e9;
400  fDChi2Tmupr = 1.0e9;
401  fDChi2Tpipr = 1.0e9;
402  fTrackPDG = 0;
403 
404 
405 
406  //
407  // energy estimator variables:
408  //
409  fE = -5.0;
410  fEsq = -5.0;
411  fPmu = -5.0;
412  fDirZmu = -5.0;
413  fN3Dprongs = -5.0;
414  fE3Dprongs = -5.0;
415  fEremain = -5.0;
416  fSumPE = -5.0;
417  fEhadMu = -5.0;
418  fMinX = 1.0e9;
419  fMinY = 1.0e9;
420  fMinZ = 1.0e9;
421  fMaxX = -1.0e9;
422  fMaxY = -1.0e9;
423  fMaxZ = -1.0e9;
424  fNHitX = -5;
425  fNHitY = -5;
426  fCCNC = -5;
427  fIntType = -5;
428  fnuPDG = 0;
429  fBestBPFmuPID = -5.0;
430  fBestRemid = -5.0;
431  fBestBPFmuPIDpdg = 0;
432 
433  }
434 
435  //......................................................................
436  void BPFTmvaTrainer::setTruthVars(const std::vector<cheat::NeutrinoEffPur>& nuEP)
437  {
438  //
439  // This function sets any variables for the TTree that come from truth.
440  // Currently, the most pure neutrino is assumed to be the best
441  // match.
442  //
443 
444  float mostPure = 0.0;
445  unsigned int match = 0;
446  for(unsigned int i = 0; i < nuEP.size(); ++i) {
447  if(nuEP[i].purity > mostPure && nuEP[i].neutrinoInt->NeutrinoSet()) {
448  mostPure = nuEP[i].purity;
449  match = i;
450  }
451  }
452 
453  if(mostPure > 0 && nuEP[match].neutrinoInt->NeutrinoSet()) {
454  simb::MCNeutrino nu = nuEP[match].neutrinoInt->GetNeutrino();
455  simb::MCParticle Nu = nu.Nu();
456  fE = Nu.E();
457  fEsq = fE*fE;
458  fCCNC = nu.CCNC();
459  fIntType = nu.InteractionType();
460  fnuPDG = Nu.PdgCode();
461  }
462 
463  }
464 
465  //......................................................................
467  {
468 
469  //
470  // General note on efficiency:
471  //
472  // This module has two training modes: PID and energy estimator.
473  // Currently, all of the information for both training modes is
474  // gathered and only the info needed for the specified training mode
475  // is stored in the TTree output. This is inefficient but this module
476  // is relatively fast anyway...
477 
478  //
479  // NOTE: Currently for the muon PID training, we only use the variables
480  // fLength, fChi2T, fdEdxLL, and fHitRatio. All other variables were
481  // found to be less effective but further exploration in the future
482  // could be done to make them better.
483  //
484 
485  // get the slices
487  evt.getByToken(fSlicerToken, slices);
488 
489  // create the FindMany object for getting prongs
491 
492  // create the FindMany object for getting Kalman tracks (used later for getting remid values)
493  art::FindManyP<rb::Track> KtrackFMP(slices, evt, "kalmantrackmerge");
494 
495  // Make the required service handles...
499 
500  // Put the slices into a PtrVector and get the list neutrinos matched
501  // to each slice from backtracker.
502  art::PtrVector<rb::Cluster> ARTslices;
503  for(unsigned int i = 0; i < slices->size(); ++i) {
504  ARTslices.push_back(art::Ptr<rb::Cluster>(slices,i));
505  }
506  std::vector< std::vector< cheat::NeutrinoEffPur > > ALLnuEffPur;
507  ALLnuEffPur = BT->SlicesToMCTruthsTable(ARTslices);
508 
509 
510 
511  // loop over all slices and get the fuzzyk 3D prongs
512  for(unsigned int islice = 0; islice < slices->size(); ++islice) {
513 
514  // As usual, skip the noise slice...
515  if((*slices)[islice].IsNoise()) continue;
516 
517  // NOTE: Currently, the fundamental object of analysis is the slice.
518  // If in the future we change things so that the fundamental object
519  // of analysis becomes the vertex, then we will have to get prongs
520  // associated with the vertex instead.
521 
522  // Reset all TTree variables.
523  this->resetVars();
524 
525  // Set slice level and truth variables for the TTrees.
526  fMinX = (*slices)[islice].MinX();
527  fMinY = (*slices)[islice].MinY();
528  fMinZ = (*slices)[islice].MinZ();
529  fMaxX = (*slices)[islice].MaxX();
530  fMaxY = (*slices)[islice].MaxY();
531  fMaxZ = (*slices)[islice].MaxZ();
532  fNHitX = (*slices)[islice].NXCell();
533  fNHitY = (*slices)[islice].NYCell();
534  this->setTruthVars(ALLnuEffPur[islice]);
535 
536 
537 
538  // get the 3D prongs associated with this slice
539  std::vector< art::Ptr< rb::Prong > > prongs;
540  if(prongFMP.isValid()) prongs = prongFMP.at(islice);
541 
542  // create the FindMany object for getting tracks
543  art::FindManyP<rb::Track> trackFMP(prongs, evt, fTrackLabel);
544 
545 
546 
547  // Variables used to keep track of the longest track (for the lengthRatio variable)
548  // and for the track with the best muon PID.
549  float longest = -1.0; // Length of longest track.
550  float bestPID = -1.0; // best muonPID value
551  int prongPick = -1; // Prong picked as having the longest track.
552 
553 
554 
555  //
556  // For the LengthRatio variable, determine first the lenght of the
557  // longest track.
558  //
559  for(unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
560 
561  std::vector< art::Ptr< rb::Track > > tracks;
562  if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
563  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
564 
565  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
566 
567  std::vector< art::Ptr< rb::FitSum > > fitsums;
568  if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
569  if(fitsums.size() != 1) abort();
570 
571  float len = tracks[itrack]->TotalLength();
572  if(len > longest && fitsums[0]->PDG() == fPIDmode) {
573  longest = len;
574  }
575 
576  } // end loop over tracks
577  } // end loop over prongs
578 
579 
580 
581  //////////
582  //
583  // Step 1: Determine which prong is the most muon-like.
584  //
585  //////////
586 
587  // loop over all prongs and get the tracks
588  for(unsigned int iprong = 0; iprong < prongs.size(); ++iprong) {
589 
590  // get the tracks associated with this prong
591  std::vector< art::Ptr< rb::Track > > tracks;
592  if(trackFMP.isValid()) tracks = trackFMP.at(iprong);
593 
594  // create the FindMany object for getting FitSum objects
595  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
596 
597  // Keep track of the total chi2 values for each of the track fit assupmtions.
598  float chi2Tmu = 1.0e9;
599  float chi2Tpi = 1.0e9;
600  float chi2Tpr = 1.0e9;
601 
602  // Loop over all tracks...
603  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
604 
605  // get the FitSum object associated with this track
606  std::vector< art::Ptr< rb::FitSum > > fitsums;
607  if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
608 
609  // There can be only one!!!
610  // (something went really wrong with BPF if there's not one and
611  // only one of these per track.)
612  if(fitsums.size() != 1) abort();
613 
614 
615 
616  float LLtemp = 0.0;
617 
618  // fill in chi2 and NDOF values
619  float tempH = 0.0;
620  float tempA = 0.0;
621  float tempN = 0.0;
622  tempH = fitsums[0]->Chi2(0)+fitsums[0]->Chi2(1); // chi2 from the hits
623  tempA = fitsums[0]->Chi2(2)+fitsums[0]->Chi2(3); // chi2 from the angles
624  tempN = fitsums[0]->Ndof(0)+fitsums[0]->Ndof(1)-4; // NDOF (NOTE: currently 4 is NOT subtracted from NDOF when the fitsum object is made by BPF)
625 
626  // Keep track of the appropriate track fit info...
627  if(fitsums[0]->PDG() == fPIDmode && tempN > 0.0) {
628  fChi2H = tempH/tempN;
629  fChi2A = tempA/tempN;
630  fChi2T = fChi2H+fChi2A;
631  fLength = tracks[itrack]->TotalLength();
632  fLengthRatio = fLength/longest;
633  }
634 
635  if(fitsums[0]->PDG() == 13 && tempN > 0.0) chi2Tmu = tempH/tempN+tempA/tempN;
636  if(fitsums[0]->PDG() == 211 && tempN > 0.0) chi2Tpi = tempH/tempN+tempA/tempN;
637  if(fitsums[0]->PDG() == 2212 && tempN > 0.0) chi2Tpr = tempH/tempN+tempA/tempN;
638 
639  //
640  // Calculate dE/dx for each cell the track passes through using
641  // the helper class.
642  //
643  fdEdxCalc.computeDEDX(tracks[itrack],fitsums[0]->PDG());
644 
645  std::vector<double> dE;
646  std::vector<double> dx;
647  std::vector<double> W;
648  std::vector<double> BG;
649  std::vector<double> s; // this is a throw away variable...
650 
651  fdEdxCalc.getResults(dE,dx,W,BG,s,fdxTOL);
652 
653  //
654  // Compute the dE/dx likelihood values for all cells.
655  //
656  for(unsigned int a = 0; a < dE.size(); ++a) {
657 
658  unsigned int w = 0;
659  unsigned int x = 0;
660 
661  if(W[a] > 0.0) w = 1;
662  if(dx[a] > 8.0) x = 1;
663 
664  // For the NearDet, we will just the "close to the readout end" FarDet histos.
665  if(geom->DetId() == novadaq::cnv::kNEARDET) w = 1;
666 
667  double content = fdEdxVSbg[w][x]->GetBinContent(fdEdxVSbg[w][x]->FindBin(log(BG[a]),dE[a]/dx[a]));
668  // We can't take the log of zero, so put a minimum on it.
669  if(content < 1.0e-9) content = 1.0e-9;
670  LLtemp += log(content);
671 
672  }
673 
674  if(dE.size() > 0) {
675  // Keep track of the appropriate track fit info...
676  if(fitsums[0]->PDG() == fPIDmode) {
677  fdEdxLL = LLtemp/(float)dE.size();
678  // NOTE: For this ratio, the number of hits on the track is taken to be the
679  // number of hits for which a value of dE/dx was successfully calculated.
680  fHitRatio = (float)dE.size()/(float)prongs[iprong]->NCell();
681  }
682  }
683  // end computing dE/dx LL info...
684 
685 
686 
687  // Set the variables for the track with the best muonPID
688  if(fTrainMode == 1 && fitsums[0]->PDG() == fPIDmode) {
689  if(fdEdxLL > -1.0e9 && fChi2T < 1.0e9) {
690  // Get the PID value...
691  float pid = fMuonPID.EvaluateMVA("KNN Method");
692  if(pid > bestPID) {
693  bestPID = pid;
694  prongPick = (int)iprong;
695  fPmu = fitsums[0]->FourMom().P();
696  fDirZmu = tracks[itrack]->Dir().Z();
697  }
698  }
699  }
700 
701  } // end loop over tracks (itrack)
702 
703 
704 
705  // Set the delta-chi2 values. We will only do this if BPF succeeded in making all
706  // three types of track (muon, pion, proton) otherwise these values won't make any
707  // sense.
708  //
709  // NOTE: If in the future, BPF makes more than three types of track then this section
710  // will need to be adjusted!
711  if(tracks.size() == 3) {
712  fDChi2Tmupi = chi2Tmu - chi2Tpi;
713  fDChi2Tmupr = chi2Tmu - chi2Tpr;
714  fDChi2Tpipr = chi2Tpi - chi2Tpr;
715  }
716 
717 
718  // Determine the true particle type for this prong.
719  const std::vector< const sim::Particle * > particles = BT->HitsToParticle(prongs[iprong]->AllCells());
720  if(particles.size() > 0) fTrackPDG = particles[0]->PdgCode();
721 
722  // Fill the PID training TTree.
723  if(fLength > 0.0 && fTrainMode == 0) {
724  // WARNING to the user! We are treating particle and anti-particle the same here...
725  if(std::abs(fTrackPDG) == fPIDmode) {
726  fBPFTrackSignalTree->Fill();
727  }
728  else {
729  fBPFTrackBackgroundTree->Fill();
730  }
731  }
732 
733 
734 
735  } // end loop over prongs (iprong)
736 
737 
738 
739  //////////
740  //
741  // Step 2: Compute the values to be used for the energy estimator.
742  //
743  //////////
744 
745  //
746  // Add up the energy of the non-muon 3D prongs.
747  //
748 
749  float N3Dpr = 0.0;
750  float E3Dpr = 0.0;
751  float sumPE = 0.0;
752  rb::Cluster hits3Dprongs; // a container for all hits on 3D prongs
753 
754  for(unsigned int p = 0; p < prongs.size(); ++p) {
755 
756  // Skip 2D prongs.
757  if(prongs[p]->Is2D()) continue;
758  N3Dpr++;
759 
760  // Keep track of all hits in 3D prongs (to be used later when adding up
761  // the remaining slice energy.)
762  // NOTE: FuzzyK prongs can sometimes share hits, so we have to be careful
763  // not to double count.
764  for(unsigned int h = 0; h < prongs[p]->NCell(); ++h) {
765 
766  const rb::CellHit *pronghit = prongs[p]->Cell(h).get();
767 
768  // Check to see if this hit is already in the list.
769  bool onlist = false;
770  for(unsigned int hpr = 0; hpr < hits3Dprongs.NCell(); ++hpr) {
771  const rb::CellHit *pronghit3D = hits3Dprongs.Cell(hpr).get();
772  if((*pronghit) == (*pronghit3D)) {
773  onlist = true;
774  break;
775  }
776  } // end loop over hits in 3D prongs (hpr)
777 
778  // Add the hit if its not already on the list.
779  if(!onlist) hits3Dprongs.Add(prongs[p]->Cell(h));
780 
781  } // end loop over hits (h)
782 
783  // Skip this prong if it is the identified muon prong.
784  if((int)p == prongPick) continue;
785 
786  // Loop over all hits in this 3D, non-muon prong and add up their energies.
787  //
788  // NOTE: For the hit W position, I am currently using the prong W instead of
789  // the track W. Since I haven't determined which of the three BPF fits
790  // is the best one for this prong, this saves me from having to try
791  // to guess at which track to use. My assumption is that prong->W will
792  // be good enough since I am lumping all hits from 3D prongs together
793  // into one 3D prong energy number...
794  //
795  // WARNING: Since hits can be shared by fuzzyK prongs, there is the potential
796  // below to double (or triple) count energy on shared cell hits. For
797  // now, I am assuming that this is a small effect but something
798  // smarter could/should be put into place in the future.
799  for(unsigned int h = 0; h < prongs[p]->NCell(); ++h) {
800 
801  const rb::CellHit *chit = prongs[p]->Cell(h).get();
802  float W = prongs[p]->W(chit);
803  rb::RecoHit rhit(cal->MakeRecoHit(*chit,W));
804 
805  sumPE += chit->PE();
806  if(rhit.IsCalibrated()) E3Dpr += rhit.GeV();
807 
808  } // end loop over hits (h)
809 
810  } // end loop over prongs (p)
811 
812  // Set the prong variables.
813  fN3Dprongs = N3Dpr;
814  fE3Dprongs = E3Dpr;
815 
816 
817 
818  //
819  // Add up the remaining energy in the slice.
820  //
821 
822  rb::Cluster hitsRemain; // a container for all hits NOT on 3D prongs
823 
824  // Loop over all hits in this slice.
825  for(unsigned int hsl = 0; hsl < (*slices)[islice].NCell(); ++hsl) {
826 
827  const rb::CellHit *slicehit = (*slices)[islice].Cell(hsl).get();
828 
829  // Check to see if this hit is in the list of hits on 3D prongs.
830  bool onlist = false;
831  for(unsigned int hpr = 0; hpr < hits3Dprongs.NCell(); ++hpr) {
832  const rb::CellHit *pronghit = hits3Dprongs.Cell(hpr).get();
833  if((*pronghit) == (*slicehit)) {
834  onlist = true;
835  break;
836  }
837  } // end loop over hits in 3D prongs (hpr)
838 
839  // If this hit is not in a 3D prong, add it to the list of remaining slice hits.
840  if(!onlist) {
841  hitsRemain.Add((*slices)[islice].Cell(hsl));
842  sumPE += slicehit->PE();
843  }
844 
845  } // end loop over hits in this slice (hsl)
846 
847  // Set the remaining slice energy variables.
848  if(hitsRemain.NCell() > 0) fEremain = hitsRemain.TotalGeV();
849  else fEremain = 0.0;
850 
851  fSumPE = sumPE;
852 
853 
854 
855  //
856  // Compute the hadronic contamination on the muon track and add that back in.
857  //
858  float muonTrackHadE = 0.0;
859 
860  // Only do this if we actually picked a track to be the muon.
861  if(prongPick >= 0) {
862 
863  std::vector< art::Ptr< rb::Track > > tracks;
864  if(trackFMP.isValid()) tracks = trackFMP.at(prongPick);
865  art::FindManyP<rb::FitSum> fitsumFMP(tracks, evt, fTrackLabel);
866 
867  // Figure out which BPF track was made with the muon assumption.
868  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack) {
869 
870  std::vector< art::Ptr< rb::FitSum > > fitsums;
871  if(fitsumFMP.isValid()) fitsums = fitsumFMP.at(itrack);
872  if(fitsums.size() != 1) abort();
873 
874  // Skip if not the muon track...
875  if(fitsums[0]->PDG() != 13) continue;
876 
877  // If it is the muon track, do the extra Ehad calculation
878  muonTrackHadE = fTrackCleanUpAlg.ExtraEOnTrackInGeV(*tracks[itrack],*ARTslices[islice]);
879 
880  } // end loop over tracks (itrack)
881 
882  fEhadMu = muonTrackHadE;
883 
884  } // end if(prongPick >= 0)
885 
886  fEremain += muonTrackHadE;
887 
888 
889 
890  //
891  // Set the PID variables for this event:
892  //
893 
894  // set the best BPF muon PID value
895  fBestBPFmuPID = bestPID;
896 
897  // loop over Kalman tracks to find the best remid value
898  if(KtrackFMP.isValid()) {
899 
900  std::vector< art::Ptr<rb::Track> > KalTracks = KtrackFMP.at(islice);
901  art::FindManyP<remid::ReMId> remidFMP(KalTracks, evt, "remid");
902  float bestRemid = -5.0;
903 
904  if(!remidFMP.isValid()) continue;
905 
906  for(unsigned int ikt = 0; ikt < KalTracks.size(); ++ikt) {
907  std::vector< art::Ptr<remid::ReMId> > remids = remidFMP.at(ikt);
908  // should only be one remid value per track
909  if(remids.size() != 1) continue;
910  if(remids[0]->Value() > bestRemid) bestRemid = remids[0]->Value();
911  }
912 
913  fBestRemid = bestRemid;
914 
915  }
916 
917 
918 
919  // Record the pdg code for the particle that contributed the most energy
920  // to the track with the highest BPF muon PID value.
921  if(prongPick >= 0) {
922  const std::vector< const sim::Particle * > particles = BT->HitsToParticle(prongs[prongPick]->AllCells());
923  int pdg = 0;
924  if(particles.size() > 0) pdg = abs(particles[0]->PdgCode());
926  }
927 
928 
929 
930  // If there is a neutrino match to this slice and we have at least one
931  // 3D prong for which BPF made at least one track under the muon assumption,
932  // fill the TTree...
933  if(fE > 0.0 && fPmu > 0.0 && fTrainMode == 1) fBPFSliceTree->Fill();
934 
935  } // end loop over slices (islice)
936 
937  return;
938 
939  }
940 
941 } // end namespace bpfit
942 ////////////////////////////////////////////////////////////////////////
943 
944 
945 
947 
948 namespace bpfit
949 {
951 }
const XML_Char int len
Definition: expat.h:262
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
unsigned int fTrainMode
training mode (energy estimator or muon PID?)
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
std::string fXMLFile
location of the XML file for the muon PID
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::string fHistoFile
location of dE/dx log-likelihood histogram
double fTotalPOT
total POT counted for all events in this job
TTree * fBPFSliceTree
TTree with slice level info (for energy estimator training)
murem::TrackCleanUpAlg fTrackCleanUpAlg
object to determine hadronic energy on the muon track
float fEsq
true energy squared (used to estiamte st. dev.)
float fDChi2Tmupi
delta chi2 total muon - pion
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
int fBestBPFmuPIDpdg
pdg code for the particle that contributed the most energy to the track with the best BPF muon PID ...
int fnuPDG
PDG code for the best matched nu (from truth)
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void analyze(const art::Event &evt)
A collection of associated CellHits.
Definition: Cluster.h:47
OStream cerr
Definition: OStream.cxx:7
int fCCNC
is this CC or NC (from truth)
float fMaxY
maximum hit Y value in the slice
void abs(TH1 *hist)
void computeDEDX(art::Ptr< rb::Track > &track, int pdg=13)
Compute dE/dx for all cells along this track and the fitsum object that went with it...
float fSumPE
sum of the PE for all hits not on the muon track
int fNHitY
number of Y-view hits in the slice
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
float fDChi2Tpipr
delta chi2 total pion - proton
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
bpfit::dEdxCalculator fdEdxCalc
helper class for computing dEdx values
double dE
float abs(float number)
Definition: d0nt_math.hpp:39
int fPIDmode
make the training TTree using the muon, pion, or proton fit
float fMaxX
maximum hit X value in the slice
float fLengthRatio
ratio of track to longest track in the slice
float fN3Dprongs
total number of 3D prongs
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
TH1D * fPOT
Total POT accumulated for all events processed in this job.
const art::ProductToken< sumdata::POTSum > fPOTSumToken
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
Calculates dE/dx in cells passed through by a track.
float fChi2H
chi^2 hits value per DOF
std::string fProngInstance
instance label for the prongs to be used
unsigned short Cell() const
Definition: CellHit.h:40
float fMinZ
minimum hit Z value in the slice
const XML_Char * s
Definition: expat.h:262
TTree * fBPFTrackBackgroundTree
TTree with track level background info (for PID training)
int InteractionType() const
Definition: MCNeutrino.h:150
TMVA::Reader fMuonPID
TMVA reader object for the muon PID.
int fNHitX
number of X-view hits in the slice
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
double dx[NP][NC]
string infile
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
float fE3Dprongs
sum of hit energies in 3D prongs (excluding the "muon" prong) - this does NOT account for dead materi...
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
#define pot
float PE() const
Definition: CellHit.h:42
Near Detector in the NuMI cavern.
float fDirZmu
muon track dirZ from the track selected to be the muon
float fEremain
sum of remaining slice energy not in 3D prongs PLUS extra energy added back in from TrackCleanUpAlg ...
caf::StandardRecord * sr
float fTrackPDG
PDG code for the primary track particle.
double fdxTOL
lower limit for dx values to be used in dE/dx calculation
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Perform a "2 point" Hough transform on a collection of hits.
virtual void endSubRun(const art::SubRun &sr)
TH2F * fdEdxVSbg[2][2]
histo of computed dE/dx for cell hits (vs. beta*gamma)
void setTruthVars(const std::vector< cheat::NeutrinoEffPur > &nuEP)
std::string fProngLabel
label for module that made the prongs
std::string fTrackLabel
label for module that made the tracks
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float fdEdxLL
dEdx LL value
float fBestBPFmuPID
best BPF muon PID value
float fPmu
reconstructed muon momentum from the track selected to be the muon
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
T const * get() const
Definition: Ptr.h:321
float fHitRatio
ratio of track hits to prong hits (used as a track/shower discriminator)
void geom(int which=0)
Definition: geom.C:163
float fMaxZ
maximum hit Z value in the slice
float fMinX
minimum hit X value in the slice
float fEhadMu
Hadronic energy contamination on the muon track as computed by TrackCleanUpAlg.
TODO.
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
BPFTmvaTrainer(fhicl::ParameterSet const &pset)
const art::ProductToken< std::vector< rb::Cluster > > fSlicerToken
float fDChi2Tmupr
delta chi2 total muon - proton
int fIntType
neutrino interaction type (from truth)
void reconfigure(const fhicl::ParameterSet &pset)
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
TTree * fBPFTrackSignalTree
TTree with track level signal info (for PID training)
Float_t w
Definition: plot.C:20
float fBestRemid
best remid value
#define W(x)
float ExtraEOnTrackInGeV(const rb::Track &track, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
ProductToken< T > consumes(InputTag const &)
float fChi2A
chi^2 angle value per DOF
void getResults(std::vector< double > &dEtemp, std::vector< double > &dxtemp, std::vector< double > &Wtemp, std::vector< double > &bgtemp, std::vector< double > &stemp, double dxTol)
Return the results of the dE/dx calculation.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
float fMinY
minimum hit Y value in the slice
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
enum BeamMode string