All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRUMBS_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CRUMBS
3 // Plugin Type: producer
4 // File: CRUMBS_module.cc
5 //
6 // Generated at Wed Jan 5 08:25:29 2022 by Henry Lay using cetskelgen
7 // from version .
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDProducer.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include "art/Framework/Principal/Event.h"
13 #include "art/Framework/Principal/Handle.h"
14 #include "art/Framework/Principal/Run.h"
15 #include "art/Framework/Principal/SubRun.h"
16 #include "canvas/Utilities/InputTag.h"
17 #include "fhiclcpp/ParameterSet.h"
19 #include "messagefacility/MessageLogger/MessageLogger.h"
20 
21 #include "canvas/Persistency/Common/FindManyP.h"
22 #include "canvas/Persistency/Common/FindOneP.h"
23 #include "art_root_io/TFileService.h"
24 
25 #include "nusimdata/SimulationBase/MCTruth.h"
26 #include "nusimdata/SimulationBase/MCParticle.h"
27 
35 
38 
44 
45 #include "TTree.h"
46 #include "TMVA/Reader.h"
47 
48 class CRUMBS;
49 
50 namespace sbn {
51  class CRUMBS : public art::EDProducer {
52  public:
53  explicit CRUMBS(fhicl::ParameterSet const& p);
54  // The compiler-generated destructor is fine for non-base
55  // classes without bare pointers or other resource use.
56 
57  // Plugins should not be copied or assigned.
58  CRUMBS(CRUMBS const&) = delete;
59  CRUMBS(CRUMBS&&) = delete;
60  CRUMBS& operator=(CRUMBS const&) = delete;
61  CRUMBS& operator=(CRUMBS&&) = delete;
62 
63  // Required functions.
64  void produce(art::Event& e) override;
65 
66  void InitialiseMVAReader(TMVA::Reader &mvaReader, std::string &mvaName, std::string &mvaFileName);
67 
68  void ResetVars();
69 
70  void GetMaps(art::Event const& e, std::map<int, int> &trackIDToGenMap, std::map<int, std::string> &genTypeMap,
71  std::map<int, int> &genCCNCMap, std::map<int, int> &genNuTypeMap);
72 
73  art::Ptr<recob::PFParticle> GetSlicePrimary(art::Event const& e,
74  const art::Ptr<recob::Slice> &slice,
75  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
76 
77  std::vector<art::Ptr<anab::T0> > GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice,
78  const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
79  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
80 
81  std::vector<art::Ptr<anab::T0> > GetCRTHitT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice,
82  const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
83  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
84 
85  float GetLongestTrackStoppingChi2Ratio(art::Event const& e, const art::Ptr<recob::Slice> &slice,
86  const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
87  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
88 
89  void FillCRTVars(const std::vector<art::Ptr<anab::T0> > &trackT0s, const std::vector<art::Ptr<anab::T0> > &hitT0s);
90 
91  void FillPandoraNuScoreVars(std::map<std::string, float> &propertiesMap);
92 
93  std::vector<art::Ptr<recob::Hit> > GetAllSliceHits(art::Event const& e,
94  const art::Ptr<recob::Slice> &slice,
95  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices);
96 
97  void GetTruthMatching(art::Event const& e, const std::vector<art::Ptr<recob::Hit> > &sliceHits, const std::vector<art::Ptr<recob::Hit> > &allHits,
98  std::map<int, int> &trackIDToGenMap, int &matchedID, double &purity, double &completeness);
99 
100  int SliceTruthId(std::map<int, float> &purities);
101 
102  private:
103 
104  // Bools to control training
106 
107  // Module labels
110 
111  // MVA location and type for loading
113 
114  // Parameter set to pass to the stopping chi2 alg
115  fhicl::ParameterSet fChi2FitParams;
116 
117  // Tree for storing training information
118  TTree *fSliceTree;
119 
120  // TMVA reader for calculating CRUMBS score
122 
123  // Other useful information for training tree
124  float tpc_NuScore;
126  int ccnc, nutype;
127  std::string matchedType;
129 
130  // Algorithms used for calculating variables
133 
134  // ======================== //
135  // CRUMBS INPUT VARIABLES //
136  // ======================== //
137 
138  // Pandora Cosmic Hypothesis Variables
139  float tpc_CRFracHitsInLongestTrack; // fraction of slice’s space points in longest track
140  float tpc_CRLongestTrackDeflection; // 1 - the cosine of the angle between the starting and finishing directions of the longest track
141  float tpc_CRLongestTrackDirY; // relative direction of the longest track in Y
142  float tpc_CRNHitsMax; // the number of space points in the largest pfp
143 
144  // Pandora Neutrino Hypothesis Variables
145  float tpc_NuEigenRatioInSphere; // the ratio between the first and second eigenvalues from a PCA of spacepoints within 10cm of the vertex
146  float tpc_NuNFinalStatePfos; // the number of final state pfos
147  float tpc_NuNHitsTotal; // the total number of space points
148  float tpc_NuNSpacePointsInSphere; // the total number of space points within 10cm of the vertex
149  float tpc_NuVertexY; // the vertex position in Y [cm]
150  float tpc_NuWeightedDirZ; // the Z component of the space-point weighted direction of the final state pfos
151 
152  // Other TPC Variables
153  float tpc_StoppingChi2CosmicRatio; // a ratio of chi2 values intended to find Bragg peaks in stopping muon tracks
154 
155  // SBN Simple Flash Match Variables
156  float pds_FMTotalScore; // the total score
157  float pds_FMPE; // the total number of photoelectrons in the associated flash
158  float pds_FMTime; // the time associated with the flash [us]
159 
160  // CRT Track and Hit Matching Variables
161  float crt_TrackScore; // a combination of the DCA and angle between the best matched TPC & CRT tracks
162  float crt_HitScore; // the best distance from an extrapolated TPC track to a CRT hit [cm]
163  float crt_TrackTime; // the time associated with the matched CRT track [us]
164  float crt_HitTime; // the time associated with the matched CRT hit [us]
165  };
166 
167 
168  CRUMBS::CRUMBS(fhicl::ParameterSet const& p)
169  : EDProducer{p},
170  fTrainingMode (p.get<bool>("TrainingMode",false)),
171  fEvaluateResultInTrainingMode (p.get<bool>("EvaluateResultInTrainingMode",false)),
172  fProcessNeutrinos (p.get<bool>("ProcessNeutrinos",true)),
173  fProcessCosmics (p.get<bool>("ProcessCosmics",true)),
174  fMCParticleModuleLabel (p.get<std::string>("MCParticleModuleLabel","")),
175  fGeneratorModuleLabel (p.get<std::string>("GeneratorModuleLabel","")),
176  fCosmicModuleLabel (p.get<std::string>("CosmicModuleLabel","")),
177  fPFParticleModuleLabel (p.get<std::string>("PFParticleModuleLabel")),
178  fHitModuleLabel (p.get<std::string>("HitModuleLabel")),
179  fTrackModuleLabel (p.get<std::string>("TrackModuleLabel")),
180  fSliceModuleLabel (p.get<std::string>("SliceModuleLabel")),
181  fFlashMatchModuleLabel (p.get<std::string>("FlashMatchModuleLabel")),
182  fCRTTrackMatchModuleLabel (p.get<std::string>("CRTTrackMatchModuleLabel")),
183  fCRTHitMatchModuleLabel (p.get<std::string>("CRTHitMatchModuleLabel")),
184  fCalorimetryModuleLabel (p.get<std::string>("CalorimetryModuleLabel")),
185  fMVAName (p.get<std::string>("MVAName")),
186  fMVAFileName (p.get<std::string>("MVAFileName")),
187  fCCNuMuMVAName (p.get<std::string>("CCNuMuMVAName")),
188  fCCNuMuMVAFileName (p.get<std::string>("CCNuMuMVAFileName")),
189  fCCNuEMVAName (p.get<std::string>("CCNuEMVAName")),
190  fCCNuEMVAFileName (p.get<std::string>("CCNuEMVAFileName")),
191  fNCMVAName (p.get<std::string>("NCMVAName")),
192  fNCMVAFileName (p.get<std::string>("NCMVAFileName")),
193  fChi2FitParams (p.get<fhicl::ParameterSet>("Chi2FitParams")),
194  fTrackStoppingChi2Alg(fChi2FitParams)
195  {
196  if(!fTrainingMode || fEvaluateResultInTrainingMode)
197  {
198  produces<std::vector<CRUMBSResult>>();
199  produces<art::Assns<recob::Slice, CRUMBSResult>>();
200 
201  InitialiseMVAReader(fMVAReader, fMVAName, fMVAFileName);
202  InitialiseMVAReader(fCCNuMuMVAReader, fCCNuMuMVAName, fCCNuMuMVAFileName);
203  InitialiseMVAReader(fCCNuEMVAReader, fCCNuEMVAName, fCCNuEMVAFileName);
204  InitialiseMVAReader(fNCMVAReader, fNCMVAName, fNCMVAFileName);
205  }
206 
207  art::ServiceHandle<art::TFileService> tfs;
208  if(fTrainingMode)
209  {
210  fSliceTree = tfs->make<TTree>("SliceTree","Slice data TTree");
211 
212  fSliceTree->Branch("tpc_NuScore",&tpc_NuScore);
213  fSliceTree->Branch("tpc_CRFracHitsInLongestTrack",&tpc_CRFracHitsInLongestTrack);
214  fSliceTree->Branch("tpc_CRLongestTrackDeflection",&tpc_CRLongestTrackDeflection);
215  fSliceTree->Branch("tpc_CRLongestTrackDirY",&tpc_CRLongestTrackDirY);
216  fSliceTree->Branch("tpc_CRNHitsMax",&tpc_CRNHitsMax);
217  fSliceTree->Branch("tpc_NuEigenRatioInSphere",&tpc_NuEigenRatioInSphere);
218  fSliceTree->Branch("tpc_NuNFinalStatePfos",&tpc_NuNFinalStatePfos);
219  fSliceTree->Branch("tpc_NuNHitsTotal",&tpc_NuNHitsTotal);
220  fSliceTree->Branch("tpc_NuNSpacePointsInSphere",&tpc_NuNSpacePointsInSphere);
221  fSliceTree->Branch("tpc_NuVertexY",&tpc_NuVertexY);
222  fSliceTree->Branch("tpc_NuWeightedDirZ",&tpc_NuWeightedDirZ);
223  fSliceTree->Branch("tpc_StoppingChi2CosmicRatio",&tpc_StoppingChi2CosmicRatio);
224 
225  fSliceTree->Branch("pds_FMTotalScore",&pds_FMTotalScore);
226  fSliceTree->Branch("pds_FMPE",&pds_FMPE);
227  fSliceTree->Branch("pds_FMTime",&pds_FMTime);
228 
229  fSliceTree->Branch("crt_TrackScore",&crt_TrackScore);
230  fSliceTree->Branch("crt_HitScore",&crt_HitScore);
231  fSliceTree->Branch("crt_TrackTime",&crt_TrackTime);
232  fSliceTree->Branch("crt_HitTime",&crt_HitTime);
233 
234  fSliceTree->Branch("eventID",&eventID);
235  fSliceTree->Branch("subRunID",&subRunID);
236  fSliceTree->Branch("runID",&runID);
237  fSliceTree->Branch("slicePDG",&slicePDG);
238  fSliceTree->Branch("matchedType",&matchedType);
239  fSliceTree->Branch("matchedPurity",&matchedPurity);
240  fSliceTree->Branch("matchedCompleteness",&matchedCompleteness);
241  fSliceTree->Branch("ccnc",&ccnc);
242  fSliceTree->Branch("nutype",&nutype);
243  }
244  }
245 
246  void CRUMBS::InitialiseMVAReader(TMVA::Reader &mvaReader, std::string &mvaName, std::string &mvaFileName)
247  {
248  mvaReader.AddVariable("tpc_CRFracHitsInLongestTrack",&tpc_CRFracHitsInLongestTrack);
249  mvaReader.AddVariable("tpc_CRLongestTrackDeflection",&tpc_CRLongestTrackDeflection);
250  mvaReader.AddVariable("tpc_CRLongestTrackDirY",&tpc_CRLongestTrackDirY);
251  mvaReader.AddVariable("tpc_CRNHitsMax",&tpc_CRNHitsMax);
252  mvaReader.AddVariable("tpc_NuEigenRatioInSphere",&tpc_NuEigenRatioInSphere);
253  mvaReader.AddVariable("tpc_NuNFinalStatePfos",&tpc_NuNFinalStatePfos);
254  mvaReader.AddVariable("tpc_NuNHitsTotal",&tpc_NuNHitsTotal);
255  mvaReader.AddVariable("tpc_NuNSpacePointsInSphere",&tpc_NuNSpacePointsInSphere);
256  mvaReader.AddVariable("tpc_NuVertexY",&tpc_NuVertexY);
257  mvaReader.AddVariable("tpc_NuWeightedDirZ",&tpc_NuWeightedDirZ);
258  mvaReader.AddVariable("tpc_StoppingChi2CosmicRatio",&tpc_StoppingChi2CosmicRatio);
259 
260  mvaReader.AddVariable("pds_FMTotalScore",&pds_FMTotalScore);
261  mvaReader.AddVariable("pds_FMPE",&pds_FMPE);
262  mvaReader.AddVariable("pds_FMTime",&pds_FMTime);
263 
264  mvaReader.AddVariable("crt_TrackScore",&crt_TrackScore);
265  mvaReader.AddVariable("crt_HitScore",&crt_HitScore);
266  mvaReader.AddVariable("crt_TrackTime",&crt_TrackTime);
267  mvaReader.AddVariable("crt_HitTime",&crt_HitTime);
268 
269  cet::search_path searchPath("FW_SEARCH_PATH");
270  std::string weightFileFullPath;
271  if (!searchPath.find_file(mvaFileName, weightFileFullPath))
272  throw cet::exception("CRUMBS") << "Unable to find weight file: " << mvaFileName << " in FW_SEARCH_PATH: " << searchPath.to_string();
273 
274  mvaReader.BookMVA(mvaName, weightFileFullPath);
275  }
276 
278  {
279  tpc_NuScore = -999999.; tpc_CRFracHitsInLongestTrack = -999999.; tpc_CRLongestTrackDeflection = -999999.; tpc_CRLongestTrackDirY = -999999.; tpc_CRNHitsMax = -999999.;
280  tpc_NuEigenRatioInSphere = -999999.; tpc_NuNFinalStatePfos = -999999.; tpc_NuNHitsTotal = -999999.; tpc_NuNSpacePointsInSphere = -999999.; tpc_NuVertexY = -999999.;
282 
283  pds_FMTotalScore = -999999.; pds_FMPE = -999999.; pds_FMTime = -500.;
284 
285  crt_TrackScore = -4.; crt_HitScore = -4.; crt_TrackTime = -3000; crt_HitTime = -3000;
286 
287  slicePDG = 999999;
288  matchedType = "";
289  matchedPurity = -999999.; matchedCompleteness = -999999.;
290  ccnc = 999999; nutype = 999999;
291  }
292 
293  void CRUMBS::GetMaps(art::Event const& e, std::map<int, int> &trackIDToGenMap, std::map<int, std::string> &genTypeMap,
294  std::map<int, int> &genCCNCMap, std::map<int, int> &genNuTypeMap)
295  {
296 
297  unsigned nNu(0), nCos(0);
298 
300  {
301  art::Handle<std::vector<simb::MCTruth> > handleMCTruthNu;
302  e.getByLabel(fGeneratorModuleLabel, handleMCTruthNu);
303  art::FindManyP<simb::MCParticle> truthNuMCPAssn(handleMCTruthNu,e,fMCParticleModuleLabel);
304 
305  for (unsigned int i = 0; i < handleMCTruthNu->size(); ++i){
306  const art::Ptr<simb::MCTruth> mcTruth(handleMCTruthNu, i);
307  const simb::MCParticle nu = mcTruth->GetNeutrino().Nu();
308 
309  if(!fTpcGeo.InVolume(nu))
310  genTypeMap[i] = "DirtNu";
311  else
312  genTypeMap[i] = "Nu";
313 
314  const std::vector<art::Ptr<simb::MCParticle> > particles = truthNuMCPAssn.at(mcTruth.key());
315 
316  for (auto const& particle : particles)
317  {
318  trackIDToGenMap[particle->TrackId()] = i;
319  }
320  ++nNu;
321 
322  genCCNCMap[i] = mcTruth->GetNeutrino().CCNC();
323  genNuTypeMap[i] = mcTruth->GetNeutrino().Nu().PdgCode();
324  }
325  }
326 
327  if(fProcessCosmics)
328  {
329  art::Handle<std::vector<simb::MCTruth> > handleMCTruthCosmic;
330  e.getByLabel(fCosmicModuleLabel, handleMCTruthCosmic);
331 
332  art::FindManyP<simb::MCParticle> truthCosmicMCPAssn(handleMCTruthCosmic,e,fMCParticleModuleLabel);
333 
334  for (unsigned int i = 0; i < handleMCTruthCosmic->size(); ++i){
335  const art::Ptr<simb::MCTruth> mcTruth(handleMCTruthCosmic, i);
336 
337  genTypeMap[i + nNu] = "Cosmic";
338 
339  const std::vector<art::Ptr<simb::MCParticle> > particles = truthCosmicMCPAssn.at(mcTruth.key());
340 
341  for (auto const& particle : particles)
342  {
343  trackIDToGenMap[particle->TrackId()] = i + nNu;
344  }
345  ++nCos;
346 
347  genCCNCMap[i + nNu] = -1;
348  genNuTypeMap[i + nNu] = -1;
349  }
350  }
351 
352  eventID = e.event();
353  subRunID = e.subRun();
354  runID = e.run();
355  }
356 
357  void CRUMBS::produce(art::Event& e)
358  {
359  std::map<int, int> trackIDToGenMap;
360  std::map<int, std::string> genTypeMap;
361  std::map<int, int> genCCNCMap;
362  std::map<int, int> genNuTypeMap;
363 
364  if(fTrainingMode)
365  this->GetMaps(e, trackIDToGenMap, genTypeMap, genCCNCMap, genNuTypeMap);
366 
367  auto resultsVec = std::make_unique<std::vector<CRUMBSResult>>();
368  auto sliceAssns = std::make_unique<art::Assns<recob::Slice, CRUMBSResult>>();
369 
370  auto const handleSlices(e.getValidHandle<std::vector<recob::Slice>>(fSliceModuleLabel));
371  std::vector<art::Ptr<recob::Slice>> slices;
372  art::fill_ptr_vector(slices, handleSlices);
373 
374  auto const handlePFPs(e.getValidHandle<std::vector<recob::PFParticle>>(fPFParticleModuleLabel));
375  std::vector<art::Ptr<recob::PFParticle>> pfps;
376  art::fill_ptr_vector(pfps, handlePFPs);
377 
378  auto const handleHits(e.getValidHandle<std::vector<recob::Hit>>(fHitModuleLabel));
379  std::vector<art::Ptr<recob::Hit>> allHits;
380  art::fill_ptr_vector(allHits, handleHits);
381 
382  art::FindManyP<larpandoraobj::PFParticleMetadata> pfpMetadataAssoc(handlePFPs, e, fPFParticleModuleLabel);
383  art::FindManyP<sbn::SimpleFlashMatch> pfpFMAssoc(handlePFPs, e, fFlashMatchModuleLabel);
384 
385  for(auto const &slice : slices)
386  {
387  this->ResetVars();
388 
389  auto const primary = this->GetSlicePrimary(e, slice, handleSlices);
390 
391  if(primary.isNull())
392  continue;
393 
394  if(primary->PdgCode() == 13 || primary->PdgCode() == 11)
395  continue;
396 
397  const std::vector<art::Ptr<larpandoraobj::PFParticleMetadata> > pfpMetaVec = pfpMetadataAssoc.at(primary.key());
398  const std::vector<art::Ptr<sbn::SimpleFlashMatch> > pfpFMVec = pfpFMAssoc.at(primary.key());
399  const std::vector<art::Ptr<anab::T0> > sliceCRTTrackT0s = this->GetCRTTrackT0s(e, slice, handlePFPs, handleSlices);
400  const std::vector<art::Ptr<anab::T0> > sliceCRTHitT0s = this->GetCRTHitT0s(e, slice, handlePFPs, handleSlices);
401 
402  this->FillCRTVars(sliceCRTTrackT0s, sliceCRTHitT0s);
403 
404  const art::Ptr<larpandoraobj::PFParticleMetadata> pfpMeta = pfpMetaVec.front();
405  std::map<std::string, float> propertiesMap = pfpMeta->GetPropertiesMap();
406 
407  this->FillPandoraNuScoreVars(propertiesMap);
408 
409  tpc_StoppingChi2CosmicRatio = this->GetLongestTrackStoppingChi2Ratio(e, slice, handlePFPs, handleSlices);
410 
411  const art::Ptr<sbn::SimpleFlashMatch> flashmatch = pfpFMVec.front();
412  pds_FMTotalScore = flashmatch->score.total;
413  pds_FMPE = flashmatch->light.pe;
414  pds_FMTime = std::max(flashmatch->time, -100.);
415 
417  {
418  const float score = fMVAReader.EvaluateMVA(fMVAName);
419  const float ccnumuscore = fCCNuMuMVAReader.EvaluateMVA(fCCNuMuMVAName);
420  const float ccnuescore = fCCNuEMVAReader.EvaluateMVA(fCCNuEMVAName);
421  const float ncscore = fNCMVAReader.EvaluateMVA(fNCMVAName);
422 
423  const float bestscore = (ccnumuscore > ccnuescore && ccnumuscore > ncscore) ? ccnumuscore : (ccnuescore > ncscore) ? ccnuescore : ncscore;
424  const int bestid = (ccnumuscore > ccnuescore && ccnumuscore > ncscore) ? 14 : (ccnuescore > ncscore) ? 12 : 1;
425 
426  resultsVec->emplace_back(score, ccnumuscore, ccnuescore, ncscore, bestscore, bestid, tpc_CRFracHitsInLongestTrack, tpc_CRLongestTrackDeflection,
431 
432  util::CreateAssn(*this, e, *resultsVec, slice, *sliceAssns);
433  }
434 
435  if(fTrainingMode)
436  {
437  std::vector<art::Ptr<recob::Hit> > sliceHits = this->GetAllSliceHits(e, slice, handleSlices);
438 
439  int matchedID(-1);
440  this->GetTruthMatching(e, sliceHits, allHits, trackIDToGenMap, matchedID, matchedPurity, matchedCompleteness);
441 
442  slicePDG = primary->PdgCode();
443  matchedType = genTypeMap[matchedID];
444 
445  ccnc = genCCNCMap[matchedID];
446  nutype = genNuTypeMap[matchedID];
447 
448  fSliceTree->Fill();
449  }
450  }
451 
453  {
454  e.put(std::move(resultsVec));
455  e.put(std::move(sliceAssns));
456  }
457  }
458 
459  void CRUMBS::FillCRTVars(const std::vector<art::Ptr<anab::T0> > &trackT0s, const std::vector<art::Ptr<anab::T0> > &hitT0s)
460  {
461  if (!trackT0s.empty()){
462  crt_TrackScore = std::numeric_limits<float>::max();
463  for(auto const crttrackmatcht0 : trackT0s)
464  {
465  if(crttrackmatcht0->TriggerConfidence() < crt_TrackScore)
466  {
467  crt_TrackScore = crttrackmatcht0->TriggerConfidence();
468  crt_TrackTime = crttrackmatcht0->Time() * 1e-3;
469  }
470  }
471  }
472 
473  if (!hitT0s.empty()){
474  crt_HitScore = std::numeric_limits<float>::max();
475  for(auto const crthitmatcht0 : hitT0s)
476  {
477  if(crthitmatcht0->TriggerConfidence() < crt_HitScore)
478  {
479  crt_HitScore = crthitmatcht0->TriggerConfidence();
480  crt_HitTime = crthitmatcht0->Time() * 1e-3;
481  }
482  }
483  }
484  }
485 
486  void CRUMBS::FillPandoraNuScoreVars(std::map<std::string, float> &propertiesMap)
487  {
488  auto propertiesMapIter = propertiesMap.find("NuScore");
489  if (propertiesMapIter == propertiesMap.end()){
490  std::cout << "CRUMBS_module: Error finding variable -- NuScore" << std::endl;
491  abort();
492  }
493  tpc_NuScore = propertiesMapIter->second;
494 
495  propertiesMapIter = propertiesMap.find("CRFracHitsInLongestTrack");
496  if (propertiesMapIter == propertiesMap.end()){
497  std::cout << "CRUMBS_module: Error finding variable -- CRFracHitsInLongestTrack" << std::endl;
498  abort();
499  }
500  tpc_CRFracHitsInLongestTrack = propertiesMapIter->second;
501 
502  propertiesMapIter = propertiesMap.find("CRLongestTrackDeflection");
503  if (propertiesMapIter == propertiesMap.end()){
504  std::cout << "CRUMBS_module: Error finding variable -- CRLongestTrackDeflection" << std::endl;
505  abort();
506  }
507  tpc_CRLongestTrackDeflection = propertiesMapIter->second;
508 
509  propertiesMapIter = propertiesMap.find("CRLongestTrackDirY");
510  if (propertiesMapIter == propertiesMap.end()){
511  std::cout << "CRUMBS_module: Error finding variable -- CRLongestTrackDirY" << std::endl;
512  abort();
513  }
514  tpc_CRLongestTrackDirY = propertiesMapIter->second;
515 
516  propertiesMapIter = propertiesMap.find("CRNHitsMax");
517  if (propertiesMapIter == propertiesMap.end()){
518  std::cout << "CRUMBS_module: Error finding variable -- CRNHitsMax" << std::endl;
519  abort();
520  }
521  tpc_CRNHitsMax = propertiesMapIter->second;
522 
523  propertiesMapIter = propertiesMap.find("NuEigenRatioInSphere");
524  if (propertiesMapIter == propertiesMap.end()){
525  std::cout << "CRUMBS_module: Error finding variable -- NuEigenRatioInSphere" << std::endl;
526  abort();
527  }
528  tpc_NuEigenRatioInSphere = propertiesMapIter->second;
529 
530  propertiesMapIter = propertiesMap.find("NuNFinalStatePfos");
531  if (propertiesMapIter == propertiesMap.end()){
532  std::cout << "CRUMBS_module: Error finding variable -- NuNFinalStatePfos" << std::endl;
533  abort();
534  }
535  tpc_NuNFinalStatePfos = propertiesMapIter->second;
536 
537  propertiesMapIter = propertiesMap.find("NuNHitsTotal");
538  if (propertiesMapIter == propertiesMap.end()){
539  std::cout << "CRUMBS_module: Error finding variable -- NuNHitsTotal" << std::endl;
540  abort();
541  }
542  tpc_NuNHitsTotal = propertiesMapIter->second;
543 
544  propertiesMapIter = propertiesMap.find("NuNSpacePointsInSphere");
545  if (propertiesMapIter == propertiesMap.end()){
546  std::cout << "CRUMBS_module: Error finding variable -- NuNSpacePointsInSphere" << std::endl;
547  abort();
548  }
549  tpc_NuNSpacePointsInSphere = propertiesMapIter->second;
550 
551  propertiesMapIter = propertiesMap.find("NuVertexY");
552  if (propertiesMapIter == propertiesMap.end()){
553  std::cout << "CRUMBS_module: Error finding variable -- NuVertexY" << std::endl;
554  abort();
555  }
556  tpc_NuVertexY = propertiesMapIter->second;
557 
558  propertiesMapIter = propertiesMap.find("NuWeightedDirZ");
559  if (propertiesMapIter == propertiesMap.end()){
560  std::cout << "CRUMBS_module: Error finding variable -- NuWeightedDirZ" << std::endl;
561  abort();
562  }
563  tpc_NuWeightedDirZ = propertiesMapIter->second;
564  }
565 
566  std::vector<art::Ptr<recob::Hit> > CRUMBS::GetAllSliceHits(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
567  {
568  art::FindManyP<recob::Hit> sliceHitAssn(handleSlices,e,fSliceModuleLabel);
569  return sliceHitAssn.at(slice.key());
570  }
571 
572  art::Ptr<recob::PFParticle> CRUMBS::GetSlicePrimary(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
573  {
574  art::FindManyP<recob::PFParticle> slicePfpAssn(handleSlices,e,fSliceModuleLabel);
575  std::vector<art::Ptr<recob::PFParticle> > pfps = slicePfpAssn.at(slice.key());
576 
577  for(auto const &pfp : pfps)
578  {
579  if(pfp->IsPrimary())
580  return pfp;
581  }
582 
583  art::Ptr<recob::PFParticle> nullReturn;
584  return nullReturn;
585  }
586 
587  void CRUMBS::GetTruthMatching(art::Event const& e, const std::vector<art::Ptr<recob::Hit> > &sliceHits, const std::vector<art::Ptr<recob::Hit> > &allHits,
588  std::map<int, int> &trackIDToGenMap, int &matchedID, double &purity, double &completeness)
589  {
590  std::map<int, int> sliceHitMap;
591  std::map<int, float> slicePurityMap;
592 
593  auto clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
594 
595  for (auto const& hit : sliceHits)
596  {
597  ++sliceHitMap[trackIDToGenMap[TruthMatchUtils::TrueParticleID(clockData,hit,true)]];
598  }
599 
600  for (auto const& [id, nHits] : sliceHitMap)
601  {
602  slicePurityMap[id] = (float) nHits / (float) sliceHits.size();
603  }
604 
605  for (auto const& [id, pur] : slicePurityMap)
606  {
607  if(pur > purity)
608  {
609  matchedID = id;
610  purity = pur;
611  }
612  }
613 
614  int totalTrueHits(0);
615 
616  for (auto const& hit : allHits)
617  {
618  if(trackIDToGenMap[TruthMatchUtils::TrueParticleID(clockData,hit,true)] == matchedID)
619  ++totalTrueHits;
620  }
621 
622  if(totalTrueHits == 0)
623  completeness = 0;
624  else
625  completeness = sliceHitMap[matchedID] / (float) totalTrueHits;
626  }
627 
628  std::vector<art::Ptr<anab::T0> > CRUMBS::GetCRTTrackT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
629  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
630  {
631  std::vector<art::Ptr<anab::T0> > t0Vec;
632 
633  art::Handle<std::vector<recob::Track> > handleTracks;
634  e.getByLabel(fTrackModuleLabel, handleTracks);
635 
636  art::FindManyP<recob::PFParticle> slicePFPAssn(handleSlices,e,fSliceModuleLabel);
637  art::FindManyP<recob::Track> pfpTrackAssn(handlePFPs,e,fTrackModuleLabel);
638  art::FindManyP<anab::T0> trackT0Assn(handleTracks,e,fCRTTrackMatchModuleLabel);
639 
640  const std::vector<art::Ptr<recob::PFParticle> > pfps = slicePFPAssn.at(slice.key());
641 
642  for(auto const& pfp : pfps)
643  {
644  if(pfp->PdgCode() != 13)
645  continue;
646 
647  const std::vector<art::Ptr<recob::Track> > tracks = pfpTrackAssn.at(pfp.key());
648 
649  if(tracks.size() != 1)
650  continue;
651 
652  const art::Ptr<recob::Track> track = tracks.front();
653 
654  const std::vector<art::Ptr<anab::T0> > t0s = trackT0Assn.at(track.key());
655  t0Vec.insert(t0Vec.end(), t0s.begin(), t0s.end());
656  }
657 
658  return t0Vec;
659  }
660 
661  std::vector<art::Ptr<anab::T0> > CRUMBS::GetCRTHitT0s(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
662  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
663  {
664  std::vector<art::Ptr<anab::T0> > t0Vec;
665 
666  art::Handle<std::vector<recob::Track> > handleTracks;
667  e.getByLabel(fTrackModuleLabel, handleTracks);
668 
669  art::FindManyP<recob::PFParticle> slicePFPAssn(handleSlices,e,fSliceModuleLabel);
670  art::FindManyP<recob::Track> pfpTrackAssn(handlePFPs,e,fTrackModuleLabel);
671  art::FindManyP<anab::T0> trackT0Assn(handleTracks,e,fCRTHitMatchModuleLabel);
672 
673  const std::vector<art::Ptr<recob::PFParticle> > pfps = slicePFPAssn.at(slice.key());
674 
675  for(auto const& pfp : pfps)
676  {
677  if(pfp->PdgCode() != 13)
678  continue;
679 
680  const std::vector<art::Ptr<recob::Track> > tracks = pfpTrackAssn.at(pfp.key());
681 
682  if(tracks.size() != 1)
683  continue;
684 
685  const art::Ptr<recob::Track> track = tracks.front();
686 
687  const std::vector<art::Ptr<anab::T0> > t0s = trackT0Assn.at(track.key());
688  t0Vec.insert(t0Vec.end(), t0s.begin(), t0s.end());
689  }
690 
691  return t0Vec;
692  }
693 
694  float CRUMBS::GetLongestTrackStoppingChi2Ratio(art::Event const& e, const art::Ptr<recob::Slice> &slice, const art::ValidHandle<std::vector<recob::PFParticle> > &handlePFPs,
695  const art::ValidHandle<std::vector<recob::Slice> > &handleSlices)
696  {
697  art::Ptr<anab::Calorimetry> longestTrackCalo;
698  float maxLength = -std::numeric_limits<float>::max();
699 
700  art::Handle<std::vector<recob::Track> > handleTracks;
701  e.getByLabel(fTrackModuleLabel, handleTracks);
702 
703  art::FindManyP<recob::PFParticle> slicePFPAssn(handleSlices,e,fSliceModuleLabel);
704  art::FindOneP<recob::Track> pfpTrackAssn(handlePFPs,e,fTrackModuleLabel);
705  art::FindManyP<anab::Calorimetry> trackCaloAssn(handleTracks,e,fCalorimetryModuleLabel);
706 
707  const std::vector<art::Ptr<recob::PFParticle> > pfps = slicePFPAssn.at(slice.key());
708 
709  for(auto const& pfp : pfps)
710  {
711  if(pfp->PdgCode() != 13)
712  continue;
713 
714  const art::Ptr<recob::Track> track = pfpTrackAssn.at(pfp.key());
715 
716  if(track.isNull())
717  continue;
718 
719  const std::vector<art::Ptr<anab::Calorimetry> > calos = trackCaloAssn.at(track.key());
720 
721  const unsigned int maxHits(std::max({ calos[0]->dEdx().size(), calos[1]->dEdx().size(), calos[2]->dEdx().size() }));
722  const int bestPlane((calos[2]->dEdx().size() == maxHits) ? 2 : (calos[0]->dEdx().size() == maxHits) ? 0 : (calos[1]->dEdx().size() == maxHits) ? 1 : -1);
723 
724  if (bestPlane == -1)
725  continue;
726 
727  const art::Ptr<anab::Calorimetry> calo = calos.at(bestPlane);
728 
729  if(track->Length() > maxLength)
730  {
731  maxLength = track->Length();
732  longestTrackCalo = calo;
733  }
734  }
735 
736  sbn::StoppingChi2Fit longestTrackFit = longestTrackCalo.isNull() ? StoppingChi2Fit() : fTrackStoppingChi2Alg.RunFitForCosmicID(*longestTrackCalo);
737 
738  if(longestTrackFit.pol0Chi2 < 0 || longestTrackFit.expChi2 <= 0)
739  return -4.f;
740 
741  return longestTrackFit.pol0Chi2 / longestTrackFit.expChi2;
742  }
743 }
744 
745 DEFINE_ART_MODULE(sbn::CRUMBS)
unsigned eventID
ClusterModuleLabel join with tracks
std::string matchedType
std::string fNCMVAFileName
float crt_HitScore
float tpc_NuWeightedDirZ
std::string fTrackModuleLabel
float tpc_NuNSpacePointsInSphere
void FillCRTVars(const std::vector< art::Ptr< anab::T0 > > &trackT0s, const std::vector< art::Ptr< anab::T0 > > &hitT0s)
Declaration of signal hit object.
std::string fMVAFileName
float crt_TrackTime
void InitialiseMVAReader(TMVA::Reader &mvaReader, std::string &mvaName, std::string &mvaFileName)
pdgs p
Definition: selectors.fcl:22
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
std::string fFlashMatchModuleLabel
G4ID TrueParticleID(detinfo::DetectorClocksData const &clockData, const art::Ptr< recob::Hit > &pHit, const bool rollupUnsavedIDs)
The G4 ID of the true particle which deposits the most energy in the recob::Hit.
float pds_FMTotalScore
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
float crt_TrackScore
std::string fSliceModuleLabel
process_name use argoneut_mc_hitfinder track
process_name hit
Definition: cheaterreco.fcl:51
std::string fHitModuleLabel
void FillPandoraNuScoreVars(std::map< std::string, float > &propertiesMap)
std::string fCosmicModuleLabel
std::string fPFParticleModuleLabel
sbn::TPCGeoAlg fTpcGeo
std::string fMVAName
bool fEvaluateResultInTrainingMode
void GetTruthMatching(art::Event const &e, const std::vector< art::Ptr< recob::Hit > > &sliceHits, const std::vector< art::Ptr< recob::Hit > > &allHits, std::map< int, int > &trackIDToGenMap, int &matchedID, double &purity, double &completeness)
TMVA::Reader fNCMVAReader
unsigned subRunID
process_name can override from command line with o or output calo
Definition: pid.fcl:40
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
sbn::TrackStoppingChi2Alg fTrackStoppingChi2Alg
std::string fCCNuMuMVAFileName
Utilities for matching a recob::Hit or vector of recob::Hit to the ID of the most significantly contr...
float tpc_StoppingChi2CosmicRatio
void ResetVars()
std::string fCRTTrackMatchModuleLabel
float tpc_CRFracHitsInLongestTrack
art::Ptr< recob::PFParticle > GetSlicePrimary(art::Event const &e, const art::Ptr< recob::Slice > &slice, const art::ValidHandle< std::vector< recob::Slice > > &handleSlices)
float tpc_CRLongestTrackDirY
fhicl::ParameterSet fChi2FitParams
std::string fNCMVAName
bool fProcessCosmics
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2687
Provides recob::Track data product.
std::vector< TCSlice > slices
Definition: DataStructs.cxx:13
std::vector< art::Ptr< anab::T0 > > GetCRTHitT0s(art::Event const &e, const art::Ptr< recob::Slice > &slice, const art::ValidHandle< std::vector< recob::PFParticle > > &handlePFPs, const art::ValidHandle< std::vector< recob::Slice > > &handleSlices)
float tpc_CRNHitsMax
float GetLongestTrackStoppingChi2Ratio(art::Event const &e, const art::Ptr< recob::Slice > &slice, const art::ValidHandle< std::vector< recob::PFParticle > > &handlePFPs, const art::ValidHandle< std::vector< recob::Slice > > &handleSlices)
float tpc_NuNHitsTotal
double matchedPurity
float tpc_NuEigenRatioInSphere
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
unsigned slicePDG
std::string fCCNuMuMVAName
int SliceTruthId(std::map< int, float > &purities)
TMVA::Reader fCCNuEMVAReader
float tpc_NuVertexY
std::string fCCNuEMVAName
StoppingChi2Fit RunFitForCosmicID(const anab::Calorimetry &calo) const
float tpc_CRLongestTrackDeflection
do i e
std::string fCalorimetryModuleLabel
std::string fCRTHitMatchModuleLabel
std::string fMCParticleModuleLabel
art::ServiceHandle< art::TFileService > tfs
void GetMaps(art::Event const &e, std::map< int, int > &trackIDToGenMap, std::map< int, std::string > &genTypeMap, std::map< int, int > &genCCNCMap, std::map< int, int > &genNuTypeMap)
bool fProcessNeutrinos
CRUMBS(fhicl::ParameterSet const &p)
double matchedCompleteness
float tpc_NuNFinalStatePfos
std::string fCCNuEMVAFileName
BEGIN_PROLOG SN nu
TMVA::Reader fMVAReader
std::string fGeneratorModuleLabel
TMVA::Reader fCCNuMuMVAReader
TTree * fSliceTree
std::vector< art::Ptr< anab::T0 > > GetCRTTrackT0s(art::Event const &e, const art::Ptr< recob::Slice > &slice, const art::ValidHandle< std::vector< recob::PFParticle > > &handlePFPs, const art::ValidHandle< std::vector< recob::Slice > > &handleSlices)
void produce(art::Event &e) override
BEGIN_PROLOG could also be cout
bool InVolume(const simb::MCParticle &particle)
std::vector< art::Ptr< recob::Hit > > GetAllSliceHits(art::Event const &e, const art::Ptr< recob::Slice > &slice, const art::ValidHandle< std::vector< recob::Slice > > &handleSlices)
unsigned runID
CRUMBS & operator=(CRUMBS const &)=delete