ProtoDUNETruthUtils.cxx
Go to the documentation of this file.
5 
11 
13 
14 #include "TVector3.h"
15 
17 
18 }
19 
21 
22 }
23 
24 std::vector< const recob::Hit* > protoana::ProtoDUNETruthUtils::FillSharedHits
25  (detinfo::DetectorClocksData const& clockData,
26  const simb::MCParticle & mcpart, const std::vector< const recob::Hit * > & hitsVec,
27  bool delta_ray, bool use_eve ) const {
28 
29  std::vector<const recob::Hit*> outVec;
30 
31  // Push back all hits that are shared with this MCParticle
33  for(const recob::Hit* hitp : hitsVec) {
34  if (use_eve) {
35  for(const sim::TrackIDE & ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
36  int trackId = ide.trackID;
37  if( !delta_ray && ( trackId == mcpart.TrackId() ) ){
38  outVec.push_back(hitp);
39  break;
40  }
41  else if( delta_ray && ( trackId == (-1*mcpart.TrackId()) ) ) {
42  outVec.push_back(hitp);
43  break;
44  }
45  }
46  }
47  else {
48  for(const int & trackId : bt_serv->HitToTrackIds(clockData, *hitp)) {
49  if( !delta_ray && ( trackId == mcpart.TrackId() ) ){
50  outVec.push_back(hitp);
51  break;
52  }
53  else if( delta_ray && ( trackId == (-1*mcpart.TrackId()) ) ) {
54  outVec.push_back(hitp);
55  break;
56  }
57  }
58  }
59  }
60 
61  return outVec;
62 }
63 
64 // Function that produces a vector of hit pointers that appear in both the given MCParticle and
65 // PFParticle.
66 std::vector<const recob::Hit*> protoana::ProtoDUNETruthUtils::GetSharedHits
67  (detinfo::DetectorClocksData const& clockData,
68  const simb::MCParticle& mcpart, const recob::PFParticle& pfpart, const art::Event& evt, std::string pfparticleModule, bool delta_ray) const {
69 
70  // Get the hits associated with this PFParticle
72  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);
73 
74  return FillSharedHits( clockData, mcpart, pfpHits, delta_ray );
75 }
76 
77 // Function that produces a vector of hit pointers that appear in both the given MCParticle and
78 // reconstructed track.
79 std::vector<const recob::Hit*> protoana::ProtoDUNETruthUtils::GetSharedHits
80  (detinfo::DetectorClocksData const& clockData,
81  const simb::MCParticle& mcpart, const recob::Track& track, const art::Event& evt, std::string trackModule, bool delta_ray) const {
82 
83  // Get the hits associated with this track
85  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);
86  return FillSharedHits( clockData, mcpart, trackHits, delta_ray );
87 }
88 
89 // Function that produces a vector of hit pointers that appear in both the given MCParticle and
90 // reconstructed shower.
91 std::vector<const recob::Hit*> protoana::ProtoDUNETruthUtils::GetSharedHits
92  (detinfo::DetectorClocksData const& clockData,
93  const simb::MCParticle& mcpart, const recob::Shower& shower, const art::Event& evt, std::string showerModule, bool delta_ray) const {
94 
95  // Get the hits associated with this shower
97  std::vector<const recob::Hit*> shHits = shUtils.GetRecoShowerHits(shower, evt, showerModule);
98  return FillSharedHits( clockData, mcpart, shHits, delta_ray );
99 }
100 
101 // Function that loops over all hits in an event and returns those that an MCParticle
102 // contributed to.
104  detinfo::DetectorClocksData const& clockData,
105  const simb::MCParticle &mcpart,
106  const art::Event &evt,
107  std::string hitModule,
108  bool use_eve) const {
109  std::vector<const recob::Hit*> outVec;
110 
111  auto hitHandle = evt.getHandle<std::vector<recob::Hit> >(hitModule);
112  if (!hitHandle) {
113  std::cout << "protoana::ProtoDUNETruthUtils::GetMCParticleHits: could not find hits in event.\n";
114  return outVec;
115  }
116 
117  // Backtrack all hits to verify whether they belong to the current MCParticle.
120  for(const recob::Hit& hit : *hitHandle) {
121  if (use_eve) {
122  for(const sim::TrackIDE & ide : bt_serv->HitToEveTrackIDEs(clockData, hit)) {
123  int trackId = ide.trackID;
124  if(pi_serv->TrackIdToParticle_P(trackId) ==
125  pi_serv->TrackIdToParticle_P(mcpart.TrackId())) {
126  outVec.push_back(&hit);
127  break;
128  }
129  }
130  }
131  else {
132  for(const int trackId : bt_serv->HitToTrackIds(clockData, hit)) {
133  if(pi_serv->TrackIdToParticle_P(trackId) ==
134  pi_serv->TrackIdToParticle_P(mcpart.TrackId())) {
135  outVec.push_back(&hit);
136  break;
137  }
138  }
139  }
140  }
141 
142  return outVec;
143 }
144 
145 // Function to return the completeness of a reconstructed object. By definition:
146 // Completeness = number of shared hits between MCParticle and reco object /
147 // number of hits in MCParticle
148 template <typename T>
149 double
151  const T &recobj,
152  const art::Event &evt,
153  std::string recoModule,
154  std::string hitModule) const
155 {
156  // Find the hits of the corresponding MCParticle.
157  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, recobj, evt, recoModule);
158  if(mcmatch == 0x0) { return 0; }
159  // std::cout << "MCmatch: " << mcmatch << '\n';
160  const std::vector<const recob::Hit*> MCHits = GetMCParticleHits(clockData, *mcmatch, evt, hitModule);
161  if(MCHits.size() == 0) { return 0; }
162  // Get shared hits between reco object and MCParticle.
163  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, recobj, evt, recoModule);
164 
165  return (double)sharedHits.size() / MCHits.size();
166 }
167 template double protoana::ProtoDUNETruthUtils::GetCompleteness<recob::PFParticle>
169 template double protoana::ProtoDUNETruthUtils::GetCompleteness<recob::Track>
170  (const detinfo::DetectorClocksData&, const recob::Track&, const art::Event&, std::string, std::string) const;
171 template double protoana::ProtoDUNETruthUtils::GetCompleteness<recob::Shower>
172  (const detinfo::DetectorClocksData&, const recob::Shower&, const art::Event&, std::string, std::string) const;
173 
174 // Function to return the purity of a PFParticle. By definition:
175 // Purity = number of shared hits between MCParticle and reco object /
176 // number of hits in reco object
177 double protoana::ProtoDUNETruthUtils::GetPurity(detinfo::DetectorClocksData const& clockData,
178  const recob::PFParticle &pfpart,
179  const art::Event &evt, std::string pfparticleModule) const {
180 
181  // Find the corresponding MCParticle and shared hits.
182  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, pfpart, evt, pfparticleModule);
183  if(mcmatch == 0x0) { return 0; }
184  // Get shared hits between PFParticle and MCParticle.
185  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, pfpart, evt, pfparticleModule);
186 
187  // Get the hits associated with this PFParticle
189  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);
190 
191  return (double)sharedHits.size() / pfpHits.size();
192 }
193 
194 // Function to return the purity of a reco track. By definition:
195 // Purity = number of shared hits between MCParticle and reco object /
196 // number of hits in reco object
197 double protoana::ProtoDUNETruthUtils::GetPurity(detinfo::DetectorClocksData const& clockData,
198  const recob::Track &track,
199  const art::Event &evt, std::string trackModule) const {
200 
201  // Find the corresponding MCParticle and shared hits.
202  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, track, evt, trackModule);
203  if(mcmatch == 0x0) { return 0; }
204  // Get shared hits between track and MCParticle.
205  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, track, evt, trackModule);
206 
207  // Get the hits associated with this track
209  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);
210 
211  return (double)sharedHits.size() / trackHits.size();
212 }
213 
214 // Function to return the purity of a reco shower. By definition:
215 // Purity = number of shared hits between MCParticle and reco object /
216 // number of hits in reco object
217 double protoana::ProtoDUNETruthUtils::GetPurity(detinfo::DetectorClocksData const& clockData,
218  const recob::Shower &shower,
219  const art::Event &evt, std::string showerModule) const {
220 
221  // Find the corresponding MCParticle and shared hits.
222  const simb::MCParticle* mcmatch = GetMCParticleFromReco(clockData, shower, evt, showerModule);
223  if(mcmatch == 0x0) { return 0; }
224  // Get shared hits between shower and MCParticle.
225  const std::vector<const recob::Hit*> sharedHits = GetSharedHits(clockData, *mcmatch, shower, evt, showerModule);
226 
227  // Get the hits associated with this shower
228  protoana::ProtoDUNEShowerUtils showerUtils;
229  std::vector<const recob::Hit*> showerHits = showerUtils.GetRecoShowerHits(shower, evt, showerModule);
230 
231  return (double)sharedHits.size() / showerHits.size();
232 }
233 
234 // Function to find a weighted and sorted vector of the true particles matched to a
235 // vector of reconstructed track hits. In case of trouble, return an empty vector.
236 std::vector<std::pair<const simb::MCParticle*, double>>
238  (detinfo::DetectorClocksData const& clockData,
239  const std::vector<const recob::Hit*>& hitVec, bool use_eve) const {
240 
241  using weightedMCPair = std::pair<const simb::MCParticle*, double>;
242  std::vector<weightedMCPair> outVec;
243 
244  // Loop over all hits in the input vector and record the contributing MCParticles.
247  std::unordered_map<const simb::MCParticle*, double> mcEMap;
248  double hitTotalE = 0;
249  for(const recob::Hit* hitp : hitVec) {
250  if (use_eve) {
251  for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
252  const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
253  mcEMap[curr_part] += ide.energy;
254  hitTotalE += ide.energy;
255  }
256  }
257  else {
258  for(const sim::TrackIDE& ide : bt_serv->HitToTrackIDEs(clockData, *hitp)) {
259  const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
260  mcEMap[curr_part] += ide.energy;
261  hitTotalE += ide.energy;
262  }
263  }
264  }
265 
266  // Fill and sort the output vector
267  for (weightedMCPair const& p : mcEMap) {
268  outVec.push_back(p);
269  }
270  std::sort(outVec.begin(), outVec.end(),
271  [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});
272 
273  // Normalise the weights by the total track energy.
274  if (hitTotalE < 1e-5) { hitTotalE = 1; } // Protect against zero division
275  for (weightedMCPair& p : outVec) {
276  p.second /= hitTotalE;
277  }
278 
279  return outVec;
280 }
281 
282 // Function to find a weighted and sorted vector of the true particles matched to a
283 // vector of reconstructed shower hits. In case of trouble, return an empty vector.
284 std::vector<std::pair<const simb::MCParticle*, double>>
286  (detinfo::DetectorClocksData const& clockData,
287  const std::vector<const recob::Hit*>& hitVec) const {
288 
289  using weightedMCPair = std::pair<const simb::MCParticle*, double>;
290  std::vector<weightedMCPair> outVec;
291 
292  // Loop over all hits in the input vector and record the contributing MCParticles.
295  std::unordered_map<const simb::MCParticle*, double> mcEMap;
296  double hitTotalE = 0;
297  for(const recob::Hit* hitp : hitVec) {
298  for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
299  const simb::MCParticle* curr_part = pi_serv->TrackIdToParticle_P(ide.trackID);
300  mcEMap[curr_part] += ide.energy;
301  hitTotalE += ide.energy;
302  }
303  }
304 
305  // Fill and sort the output vector
306  for (weightedMCPair const& p : mcEMap) {
307  outVec.push_back(p);
308  }
309  std::sort(outVec.begin(), outVec.end(),
310  [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});
311 
312  // Normalise the weights by the total track energy.
313  if (hitTotalE < 1e-5) { hitTotalE = 1; } // Protect against zero division
314  for (weightedMCPair& p : outVec) {
315  p.second /= hitTotalE;
316  }
317 
318  return outVec;
319 }
320 
321 // Function to find a weighted and sorted vector of the true particles matched to
322 // a PFParticle. Return an empty vector in case of trouble.
323 std::vector<std::pair<const simb::MCParticle*, double>>
325  (detinfo::DetectorClocksData const& clockData,
326  const recob::PFParticle &pfpart, art::Event const & evt, std::string pfparticleModule) const {
327 
328  std::vector<std::pair<const simb::MCParticle*, double>> outVec;
329 
330  // We must have MC for this module to make sense
331  if(evt.isRealData()) return outVec;
332 
333  // Get the hits associated with this PFParticle
335  std::vector<const recob::Hit*> pfpHits = pfpUtils.GetPFParticleHits(pfpart, evt, pfparticleModule);
336 
337  // Call the appropriate function to determine MC matches.
338  if(abs(pfpart.PdgCode()) == 11 || pfpart.PdgCode() == 22) {
339  outVec = GetMCParticleListFromShowerHits(clockData, pfpHits);
340  } else {
341  outVec = GetMCParticleListFromTrackHits(clockData, pfpHits);
342  }
343 
344  return outVec;
345 }
346 
347 std::vector<std::pair<const recob::PFParticle*, double>>
349  (detinfo::DetectorClocksData const& clockData,
350  const simb::MCParticle &part, art::Event const & evt, std::string pfparticleModule, bool use_eve) const {
351 
352  using weightedPFPair = std::pair<const recob::PFParticle*, double>;
353 
354  std::vector<weightedPFPair> outVec;
355 
356  // We must have MC for this module to make sense
357  if(evt.isRealData()) return outVec;
358 
359  // Get the reconstructed PFParticles
360  auto allPFParticles = evt.getValidHandle<std::vector<recob::PFParticle> >(pfparticleModule);
361 
362  // Loop over all PFParticles and save the energetic overlap with each MCParticle in a map
364  std::unordered_map<const recob::PFParticle*, double> pfpEMap;
365  double pfpTotalE = 0;
366  for(const recob::PFParticle& pfpart : *allPFParticles) {
367  for(const recob::Hit* hitp : GetSharedHits(clockData, part, pfpart, evt, pfparticleModule)) {
368  if (use_eve) {
369  for(const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, *hitp)) {
370  if(ide.trackID == part.TrackId()) {
371  pfpEMap[&pfpart] += ide.energy;
372  }
373  pfpTotalE += ide.energy;
374  }
375  }
376  else {
377  for(const sim::TrackIDE& ide : bt_serv->HitToTrackIDEs(clockData, *hitp)) {
378  if(ide.trackID == part.TrackId()) {
379  pfpEMap[&pfpart] += ide.energy;
380  }
381  pfpTotalE += ide.energy;
382  }
383  }
384  }
385  }
386 
387  // Fill and sort the output vector
388  for (weightedPFPair const& p : pfpEMap) {
389  outVec.push_back(p);
390  }
391  std::sort(outVec.begin(), outVec.end(),
392  [](weightedPFPair a, weightedPFPair b){ return a.second > b.second;});
393 
394  // Normalise the weights by the total track energy.
395  if (pfpTotalE < 1e-5) { pfpTotalE = 1; } // Protect against zero division
396  for (weightedPFPair& p : outVec) {
397  p.second /= pfpTotalE;
398  }
399 
400  return outVec;
401 }
402 
403 // Function to find a weighted and sorted vector of the true particles matched to
404 // a track. Return an empty vector in case of trouble.
405 std::vector<std::pair<const simb::MCParticle*, double>>
407  (detinfo::DetectorClocksData const& clockData,
408  const recob::Track &track, art::Event const & evt, std::string trackModule) const {
409 
410  // We must have MC for this module to make sense
411  if(evt.isRealData()) return {};
412 
413  // Get the reconstructed track hits
415  std::vector<const recob::Hit*> trackHits = trackUtils.GetRecoTrackHits(track, evt, trackModule);
416  return GetMCParticleListFromTrackHits(clockData, trackHits);
417 }
418 
419 // Function to find the best matched reconstructed particle tracks to a true
420 // particle. In case of problems, or if the particle was not reconstruced as a
421 // track, returns an empty vector.
422 std::vector<std::pair<const recob::Track*, double>>
424  (detinfo::DetectorClocksData const& clockData,
425  const simb::MCParticle &part, art::Event const & evt, std::string trackModule) const{
426 
427  using weightedTrackPair = std::pair<const recob::Track*, double>;
428 
429  std::vector<weightedTrackPair> outVec;
430 
431  // We must have MC for this module to make sense
432  if(evt.isRealData()) return outVec;
433 
434  // Get the reconstructed tracks
435  auto allRecoTracks = evt.getValidHandle<std::vector<recob::Track> >(trackModule);
436 
437  // We need the association between the tracks and the hits
438  const art::FindManyP<recob::Hit> findTrackHits(allRecoTracks, evt, trackModule);
439 
441  std::unordered_map<int, double> recoTrack;
442 
443  // Record the energy contribution to the MCParticle of every relevant reco track
444  double part_E = 0;
445  for (recob::Track const & track : *allRecoTracks)
446  {
447  // This finds all hits that are shared between a reconstructed track and MCParticle
448  for (art::Ptr<recob::Hit> const & hptr :
449  bt_serv->TrackIdToHits_Ps(clockData, part.TrackId(), findTrackHits.at(track.ID())))
450  {
451  // Loop over hit IDEs to find our particle's part in it
452  for (const sim::IDE* ide : bt_serv->HitToSimIDEs_Ps(clockData, hptr))
453  {
454  if (ide->trackID == part.TrackId())
455  {
456  recoTrack[track.ID()] += ide->energy;
457  }
458  part_E += ide->energy;
459  } // sim::IDE*
460  } // art::Ptr<recob::Hit>
461  } // const recob::Track&
462 
463  // Fill and sort the output vector
464  for (std::pair<int, double> const& p : recoTrack)
465  {
466  auto const trackIt = std::find_if(allRecoTracks->begin(), allRecoTracks->end(),
467  [&](recob::Track tr){ return tr.ID() == p.first; });
468  outVec.push_back(std::make_pair(&*trackIt, p.second));
469  }
470  std::sort(outVec.begin(), outVec.end(),
471  [](weightedTrackPair a, weightedTrackPair b){ return a.second > b.second;});
472 
473  // Normalise the vector weights
474  if (part_E < 1e-5) { part_E = 1; } // Protect against zero division
475  for (weightedTrackPair& p : outVec)
476  {
477  p.second /= part_E;
478  }
479 
480  return outVec;
481 }
482 
483 // Function to find the best matched true particle to a reconstructed particle
484 // shower. In case of problems, returns a null pointer
485 std::vector<std::pair<const simb::MCParticle*, double>>
487  (detinfo::DetectorClocksData const& clockData,
488  const recob::Shower &shower, art::Event const & evt, std::string showerModule) const{
489 
490  // We must have MC for this module to make sense
491  if(evt.isRealData()) return {};
492 
493  // Get the reconstructed shower hits
494  protoana::ProtoDUNEShowerUtils showerUtils;
495  std::vector<const recob::Hit*> showerHits = showerUtils.GetRecoShowerHits(shower, evt, showerModule);
496  return GetMCParticleListFromShowerHits(clockData, showerHits);
497 }
498 
499 // Function to find the best matched reconstructed particle tracks to a true
500 // particle. In case of problems, or if the particle was not reconstruced as a
501 // shower, returns an empty vector.
502 std::vector<std::pair<const recob::Shower*, double>>
504  (detinfo::DetectorClocksData const& clockData,
505  const simb::MCParticle &part, art::Event const & evt, std::string showerModule) const{
506 
507  using weightedShowerPair = std::pair<const recob::Shower*, double>;
508 
509  std::vector<weightedShowerPair> outVec;
510 
511  // We must have MC for this module to make sense
512  if(evt.isRealData()) return outVec;
513 
514  // Get the reconstructed showers
515  auto allRecoShowers = evt.getValidHandle<std::vector<recob::Shower> >(showerModule);
516 
517  // We need the association between the showers and the hits
518  const art::FindManyP<recob::Hit> findShowerHits(allRecoShowers, evt, showerModule);
519 
521  ProtoDUNEShowerUtils shUtils;
522  std::unordered_map<int, double> recoShower;
523 
524  // Record the energy contribution to the MCParticle of every relevant reco shower
525  double part_E = 0;
526  for (recob::Shower const & shower : *allRecoShowers)
527  {
528  const unsigned showerIndex = shUtils.GetShowerIndex(shower, evt, showerModule);
529  // Since MCParticles do not have shower hits associated with them, loop over all shower hits
530  for (art::Ptr<recob::Hit> hit : findShowerHits.at(showerIndex))
531  {
532  // Loop over hit IDEs to find our particle's part in it
533  for (const sim::TrackIDE& ide : bt_serv->HitToEveTrackIDEs(clockData, hit))
534  {
535  if (ide.trackID == part.TrackId())
536  {
537  recoShower[showerIndex] += ide.energy;
538  }
539  part_E += ide.energy;
540  } // sim::IDE*
541  } // art::Ptr<recob::Hit>
542  } // const recob::Shower&
543 
544  // Fill and sort the output vector
545  for (std::pair<int, double> const& p : recoShower)
546  {
547  auto const showerIt = std::find_if(allRecoShowers->begin(), allRecoShowers->end(),
548  [&](recob::Shower sh){ return shUtils.GetShowerIndex(sh, evt, showerModule) == p.first; });
549  outVec.push_back(std::make_pair(&*showerIt, p.second));
550  }
551  std::sort(outVec.begin(), outVec.end(),
552  [](weightedShowerPair a, weightedShowerPair b){ return a.second > b.second;});
553 
554  // Normalise the vector weights
555  if (part_E < 1e-5) { part_E = 1; } // Protect against zero division
556  std::for_each(outVec.begin(), outVec.end(),
557  [&](weightedShowerPair& p){ p.second /= part_E; });
558 
559  return outVec;
560 }
561 
562 // Function to find the best matched true particle to a reconstructed PFParticle.
563 // In case of problems, returns a null pointer
565  (detinfo::DetectorClocksData const& clockData,
566  const recob::PFParticle &pfpart, art::Event const & evt, std::string pfparticleModule) const{
567 
568  const simb::MCParticle* outPart = 0x0;
569 
570  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
571  GetMCParticleListFromPFParticle(clockData, pfpart, evt, pfparticleModule);
572 
573  if (inVec.size() != 0) outPart = inVec[0].first;
574 
575  return outPart;
576 }
577 
578 // Function to find the best matched reconstructed PFParticle to a true particle. In
579 // case of problems, or if the true particle is not a primary contributor to any
580 // PFParticle, return a null pointer
582  (detinfo::DetectorClocksData const& clockData,
583  const simb::MCParticle &part, art::Event const & evt, std::string pfparticleModule) const {
584 
585  const recob::PFParticle* outPFPart = 0x0;
586 
587  // We must have MC for this module to make sense
588  if(evt.isRealData()) return outPFPart;
589 
590  // Get the list of contributing PFParticles for the MCParticle and see if any of
591  // them have the MCParticle as a primary contributor.
592  for (std::pair<const recob::PFParticle*, double> p :
593  GetPFParticleListFromMCParticle(clockData, part, evt, pfparticleModule))
594  {
595  const recob::PFParticle* pfp = p.first;
596  if (GetMCParticleFromPFParticle(clockData, *pfp, evt, pfparticleModule)->TrackId() == part.TrackId())
597  {
598  outPFPart = pfp;
599  break;
600  }
601  }
602 
603  return outPFPart;
604 }
605 
606 // Function to find the best matched true particle to a reconstructed particle
607 // track. In case of problems, returns a null pointer
609  (detinfo::DetectorClocksData const& clockData,
610  const recob::Track &track, art::Event const & evt, std::string trackModule) const{
611 
612  const simb::MCParticle* outPart = 0x0;
613 
614  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
615  GetMCParticleListFromRecoTrack(clockData, track, evt, trackModule);
616 
617  if (inVec.size() != 0) outPart = inVec[0].first;
618 
619  return outPart;
620 }
621 
622 // Function to find the best matched reconstructed track to a true particle. In
623 // case of problems, or if the true particle is not a primary contributor to any
624 // track, return a null pointer
626  (detinfo::DetectorClocksData const& clockData,
627  const simb::MCParticle &part, art::Event const & evt, std::string trackModule) const {
628 
629  const recob::Track* outTrack = 0x0;
630 
631  // We must have MC for this module to make sense
632  if(evt.isRealData()) return outTrack;
633 
634  // Get the list of contributing tracks for the MCParticle and see if any of
635  // them have the MCParticle as a primary contributor.
636  for (std::pair<const recob::Track*, double> p :
637  GetRecoTrackListFromMCParticle(clockData, part, evt, trackModule))
638  {
639  const recob::Track* tr = p.first;
640  if (GetMCParticleFromRecoTrack(clockData, *tr, evt, trackModule)->TrackId() == part.TrackId())
641  {
642  outTrack = tr;
643  break;
644  }
645  }
646 
647  return outTrack;
648 }
649 
650 // Function to find the best matched true particle to a reconstructed particle
651 // shower. In case of problems, returns a null pointer
653  (detinfo::DetectorClocksData const& clockData,
654  const recob::Shower &shower, art::Event const & evt, std::string showerModule) const{
655 
656  const simb::MCParticle* outPart = 0x0;
657 
658  const std::vector<std::pair<const simb::MCParticle*, double>> inVec =
659  GetMCParticleListFromRecoShower(clockData, shower, evt, showerModule);
660 
661  if (inVec.size() != 0) outPart = inVec[0].first;
662 
663  return outPart;
664 }
665 
666 // Function to find the best matched reconstructed shower to a true particle. In
667 // case of problems, or if the true particle is not a primary contributor to any
668 // shower, returns a null pointer
670  (detinfo::DetectorClocksData const& clockData,
671  const simb::MCParticle &part, art::Event const & evt, std::string showerModule) const {
672 
673  const recob::Shower* outShower = 0x0;
674 
675  // We must have MC for this module to make sense
676  if(evt.isRealData()) return outShower;
677 
678  // Get the list of contributing showers for the MCParticle and see if any of
679  // them have the MCParticle as a primary contributor.
680  for (std::pair<const recob::Shower*, double> p :
681  GetRecoShowerListFromMCParticle(clockData, part, evt, showerModule))
682  {
683  const recob::Shower* tr = p.first;
684  if (GetMCParticleFromRecoShower(clockData, *tr, evt, showerModule)->TrackId() == part.TrackId())
685  {
686  outShower = tr;
687  break;
688  }
689  }
690 
691  return outShower;
692 }
693 
694 // General overloaded functions to get the best MCParticle match for reconstructed objects
696  (detinfo::DetectorClocksData const& clockData,
697  const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const {
698  return GetMCParticleFromPFParticle(clockData, pfpart, evt, pfparticleModule);
699 }
701  (detinfo::DetectorClocksData const& clockData,
702  const recob::Track &track, art::Event const &evt, std::string trackModule) const {
703  return GetMCParticleFromRecoTrack(clockData, track, evt, trackModule);
704 }
706  (detinfo::DetectorClocksData const& clockData,
707  const recob::Shower &shower, art::Event const &evt, std::string showerModule) const {
708  return GetMCParticleFromRecoShower(clockData, shower, evt, showerModule);
709 }
710 
711 // General overloaded functions to get the contributing MCParticles to a reconstructed object
712 std::vector<std::pair<const simb::MCParticle*, double>>
714  (detinfo::DetectorClocksData const& clockData,
715  const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const {
716  return GetMCParticleListFromPFParticle(clockData, pfpart, evt, pfparticleModule);
717 }
718 std::vector<std::pair<const simb::MCParticle*, double>>
720  (detinfo::DetectorClocksData const& clockData,
721  const recob::Track &track, art::Event const &evt, std::string trackModule) const {
722  return GetMCParticleListFromRecoTrack(clockData, track, evt, trackModule);
723 }
724 std::vector<std::pair<const simb::MCParticle*, double>>
726  (detinfo::DetectorClocksData const& clockData,
727  const recob::Shower &shower, art::Event const &evt, std::string showerModule) const {
728  return GetMCParticleListFromRecoShower(clockData, shower, evt, showerModule);
729 }
730 
731 // Match Reco to True by Number of Hits contributed as opposed to Energy contributed
732 //
733 // List
734 template <typename T>
735 //std::vector< std::pair< const simb::MCParticle *, size_t > >
736 std::vector< protoana::MCParticleSharedHits >
737  protoana::ProtoDUNETruthUtils::GetMCParticleListByHits(detinfo::DetectorClocksData const& clockData,
738  const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const {
739 
740  //Get the particle list. This is sorted by energy contributed, but we'll turn this into nHits
741  const std::vector< std::pair< const simb::MCParticle*, double > > mcparts = GetMCParticleListFromReco(clockData, recobj, evt, recoModule );
742 
743  //using weightedMCPair = std::pair<const simb::MCParticle*, size_t>;
744  //std::vector< weightedMCPair > results;
745  std::vector< protoana::MCParticleSharedHits > results;
746 
747  for( const std::pair< const simb::MCParticle*, double > part : mcparts ){
748  //Get the non-delta ray hits
749  const std::vector<const recob::Hit*> hits = GetSharedHits(clockData, *(part.first), recobj, evt, recoModule);
750 
751  //Get the delta ray hits
752  const std::vector<const recob::Hit*> dr_hits = GetSharedHits(clockData, *(part.first), recobj, evt, recoModule, true);
753 
755  entry.particle = part.first;
756  entry.nSharedHits = hits.size();
757  entry.nSharedDeltaRayHits = dr_hits.size();
758  results.push_back( entry );
759  //results.push_back( std::make_pair( part.first, hits.size() ) );
760  }
761 
762  //std::sort(results.begin(), results.end(),
763  // [](weightedMCPair a, weightedMCPair b){ return a.second > b.second;});
764 
765  std::sort( results.begin(), results.end(),
767  { return ( (a.nSharedHits + a.nSharedDeltaRayHits) > (b.nSharedHits + b.nSharedDeltaRayHits) ); });
768 
769  return results;
770 }
771 
772 //template std::vector< std::pair< const simb::MCParticle *, size_t > > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::PFParticle >
773 template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::PFParticle >
774  (detinfo::DetectorClocksData const&,
775  const recob::PFParticle&, const art::Event&, std::string, std::string) const;
776 //template std::vector< std::pair< const simb::MCParticle *, size_t > > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Track >
777 template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Track >
778  (detinfo::DetectorClocksData const&,
779  const recob::Track&, const art::Event&, std::string, std::string) const;
780 //template std::vector< std::pair< const simb::MCParticle *, size_t > > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Shower >
781 template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Shower >
782  (detinfo::DetectorClocksData const&,
783  const recob::Shower&, const art::Event&, std::string, std::string) const;
784 /////////////////
785 
786 // Best Match
787 template <typename T>
788 const /*simb::MCParticle **/ protoana::MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits (detinfo::DetectorClocksData const& clockData,
789  const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule ) const {
790 
791  auto list = GetMCParticleListByHits(clockData, recobj, evt, recoModule, hitModule );
792 
794  dummy.particle = 0x0;
795  dummy.nSharedHits = 0;
796  dummy.nSharedDeltaRayHits = 0;
797  //if( !list.size() ) return 0x0;
798  if( !list.size() ) return dummy;
799 
800  //return list[0].first;
801  //return list[0].particle;
802  return list[0];
803 }
804 
806 (const detinfo::DetectorClocksData&, const recob::PFParticle&, const art::Event&, std::string, std::string) const;
808  (const detinfo::DetectorClocksData&, const recob::Track&, const art::Event&, std::string, std::string) const;
810  (const detinfo::DetectorClocksData&, const recob::Shower&, const art::Event&, std::string, std::string) const;
811 /////////////////
812 
813 const simb::MCParticle* protoana::ProtoDUNETruthUtils::MatchPduneMCtoG4( const simb::MCParticle & pDunePart, const art::Event & evt )
814 { // Function that will match the protoDUNE MC particle to the Geant 4 MC particle, and return the matched particle (or a null pointer).
815 
816  // Find the energy of the procided MC Particle
817  double pDuneEnergy = pDunePart.E();
818 
819  // Get list of the g4 particles. plist should be a std::map< int, simb::MCParticle* >
821  const sim::ParticleList & plist = pi_serv->ParticleList();
822 
823  // Check if plist is empty
824  if ( !plist.size() ) {
825  std::cerr << "\n\n#####################################\n"
826  << "\nEvent " << evt.id().event() << "\n"
827  << "sim::ParticleList from cheat::ParticleInventoryService is empty\n"
828  << "A null pointer will be returned\n"
829  << "#####################################\n\n";
830  return nullptr;
831  }
832 
833  // Go through the list of G4 particles
834  for ( auto partIt = plist.begin() ; partIt != plist.end() ; partIt++ ) {
835 
836  const simb::MCParticle* pPart = partIt->second;
837  if (!pPart) {
838  std::cerr << "\n\n#####################################\n"
839  << "\nEvent " << evt.id().event() << "\n"
840  << "GEANT particle #" << partIt->first << " returned a null pointer\n"
841  << "This is not necessarily bad. It just means at least one\n"
842  << "of the G4 particles returned a null pointer. It may well\n"
843  << "have still matched a PD particle and a G4 particle.\n"
844  << "#####################################\n\n";
845  continue;
846  }
847 
848  // If the initial energy of the g4 particle is very close to the energy of the protoDUNE particle, call it a day and have a cuppa.
849  if ( (pDunePart.PdgCode() == pPart->PdgCode()) && fabs(pPart->E() - pDuneEnergy) < 0.00001 ) {
850  return pPart;
851  }
852 
853  } // G4 particle list loop end.
854 
855  std::cout << "No G4 particle was matched for Event " << evt.id().event() << ". Null pointer returned\n";
856  return nullptr;
857 
858 } // End MatchPduneMCtoG4
859 
860 const simb::MCParticle* protoana::ProtoDUNETruthUtils::GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const{
861 
862  // Get the good particle from the MCTruth
863  simb::MCParticle goodPart;
864  bool found = false;
865  for(int t = 0; t < genTruth.NParticles(); ++t){
866  simb::MCParticle part = genTruth.GetParticle(t);
867  if(part.Process() == "primary"){
868  goodPart = part;
869  found = true;
870  break;
871  }
872  }
873 
874  if(!found){
875  std::cerr << "No good particle found, returning null pointer" << std::endl;
876  return nullptr;
877  }
878 
879  // Now loop over geant particles to find the one that matches
880  // Get list of the g4 particles. plist should be a std::map< int, simb::MCParticle* >
882  const sim::ParticleList & plist = pi_serv->ParticleList();
883 
884  for(auto const part : plist){
885  if((goodPart.PdgCode() == part.second->PdgCode()) && fabs(part.second->E() - goodPart.E()) < 1e-5){
886  return part.second;
887  }
888  }
889 
890  // If we get here something has gone wrong
891  std::cerr << "No G4 version of the good particle was found, returning null pointer" << std::endl;
892  return nullptr;
893 
894 }
895 
896 // Converting times in LArSoft can be a bit of a minefield. These functions convert true times in ns
897 // to pandora times in ns
898 const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeNano(detinfo::DetectorClocksData const& clockData,
899  const simb::MCParticle &part) const{
900  return ConvertTrueTimeToPandoraTimeNano(clockData, part.T());
901 }
902 
903 const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeNano(detinfo::DetectorClocksData const& clockData,
904  const float trueTime) const{
905  return 1000. * ConvertTrueTimeToPandoraTimeMicro(clockData, trueTime);
906 }
907 
908 // Microsecond versions
909 const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeMicro(detinfo::DetectorClocksData const& clockData,
910  const simb::MCParticle &part) const{
911  return ConvertTrueTimeToPandoraTimeMicro(clockData, part.T());
912 }
913 
914 const float protoana::ProtoDUNETruthUtils::ConvertTrueTimeToPandoraTimeMicro(detinfo::DetectorClocksData const& clockData,
915  const float trueTime) const{
916 
917  return clockData.G4ToElecTime(trueTime);
918 }
919 
920 // Get process key.
922 
923  if(process.compare("primary") == 0) return 0;
924  else if(process.compare("hadElastic") == 0) return 1;
925  else if(process.compare("pi-Inelastic") == 0) return 2;
926  else if(process.compare("pi+Inelastic") == 0) return 3;
927  else if(process.compare("kaon-Inelastic") == 0) return 4;
928  else if(process.compare("kaon+Inelastic") == 0) return 5;
929  else if(process.compare("protonInelastic") == 0) return 6;
930  else if(process.compare("neutronInelastic") == 0) return 7;
931  else if(process.compare("kaon0SInelastic") == 0) return 8;
932  else if(process.compare("kaon0LInelastic") == 0) return 9;
933  else if(process.compare("lambdaInelastic") == 0) return 10;
934  else if(process.compare("omega-Inelastic") == 0) return 11;
935  else if(process.compare("sigma+Inelastic") == 0) return 12;
936  else if(process.compare("sigma-Inelastic") == 0) return 13;
937  else if(process.compare("sigma0Inelastic") == 0) return 14;
938  else if(process.compare("xi-Inelastic") == 0) return 15;
939  else if(process.compare("xi0Inelastic") == 0) return 16;
940  else if(process.compare("anti_protonInelastic") == 0) return 20;
941  else if(process.compare("anti_neutronInelastic") == 0) return 21;
942  else if(process.compare("anti_lambdaInelastic") == 0) return 22;
943  else if(process.compare("anti_omega-Inelastic") == 0) return 23;
944  else if(process.compare("anti_sigma+Inelastic") == 0) return 24;
945  else if(process.compare("anti_sigma-Inelastic") == 0) return 25;
946  else if(process.compare("anti_xi-Inelastic") == 0) return 26;
947  else if(process.compare("anti_xi0Inelastic") == 0) return 27;
948 
949  else if(process.compare("Decay") == 0) return 30;
950  else if(process.compare("FastScintillation") == 0) return 31;
951  else if(process.compare("nKiller") == 0) return 32; // Remove unwanted neutrons: neutron kinetic energy threshold (default 0) or time limit for neutron track
952  else if(process.compare("nCapture") == 0) return 33; // Neutron capture
953 
954  else if(process.compare("compt") == 0) return 40; // Compton Scattering
955  else if(process.compare("rayleigh") == 0) return 41; // Rayleigh Scattering
956  else if(process.compare("phot") == 0) return 42; // Photoelectric Effect
957  else if(process.compare("conv") == 0) return 43; // Pair production
958  else if(process.compare("CoupledTransportation") == 0) return 44; //
959 
960  else return -1;
961 }
962 
963 // Get estimated particle energy deposit. The G4 trackId must be provided
964 double protoana::ProtoDUNETruthUtils::GetDepEnergyMC(const art::Event &evt, geo::GeometryCore const * fGeom, int trackid, int whichview) const {
965 
966  double edep = 0.0;
967 
968  auto simchannelHandle = evt.getHandle< std::vector<sim::SimChannel> >("largeant");
969  if (simchannelHandle) {
970  // Loop over sim channels
971  for(auto const& simchannel : (*simchannelHandle)){
972  // Only keep channels in the selected view
973  if(fGeom->View(simchannel.Channel()) != whichview) continue;
974  // Get all time slices
975  auto const& alltimeslices = simchannel.TDCIDEMap();
976  // Loop over time slices
977  for(auto const& tslice : alltimeslices){
978  auto const& simide = tslice.second;
979  // Loop over energy deposits
980  for(auto const& eDep : simide){
981  if(eDep.trackID == trackid || eDep.trackID == -trackid)
982  edep += eDep.energy;
983  }
984  }
985  }
986  }
987 
988  return edep;
989 
990 }
991 
992 // Get first trajectory point in TPC active volume
993 int protoana::ProtoDUNETruthUtils::GetFirstTrajectoryPointInTPCActiveVolume(const simb::MCParticle& mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh){
994 
995  int firstpoint = -999;
996  for(unsigned int i = 0; i < mcpart.NumberTrajectoryPoints(); ++i) {
997  if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){
998 
999  firstpoint = i;
1000  break;
1001  }
1002  }
1003 
1004  return firstpoint;
1005 }
1006 
1007 // Get MC Particle length in TPC active volume
1008 double protoana::ProtoDUNETruthUtils::GetMCParticleLengthInTPCActiveVolume(const simb::MCParticle& mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh){
1009 
1010  double length = 0.0;
1011  int firstpoint = GetFirstTrajectoryPointInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);
1012 
1013  if(firstpoint < 0) return length;
1014 
1015  TVector3 pos = mcpart.Position(firstpoint).Vect();
1016  for(unsigned int i = firstpoint+1; i < mcpart.NumberTrajectoryPoints(); ++i) {
1017  if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){
1018 
1019  pos -= mcpart.Position(i).Vect();
1020  length += pos.Mag();
1021  pos = mcpart.Position(i).Vect();
1022  }
1023  }
1024 
1025  return length;
1026 }
1027 
1028 // Estimate last point energy loss
1029 double protoana::ProtoDUNETruthUtils::GetDepEnergyAtLastTrajPoint(const simb::MCParticle& mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh){
1030 
1031  // First trajectory point inside TPC active volume
1032  int firstpoint = GetFirstTrajectoryPointInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);
1033  if(firstpoint < 0) return 0.0;
1034 
1035  // True track length in TPC active volume
1036  double length = GetMCParticleLengthInTPCActiveVolume(mcpart, tpcactiveXlow, tpcactiveXhigh, tpcactiveYlow, tpcactiveYhigh, tpcactiveZlow, tpcactiveZhigh);
1037 
1038  // Get the true track length up to the next to the last trajectory point
1039  int ntrajpoints = mcpart.NumberTrajectoryPoints();
1040  double plength = 0.0;
1041  TVector3 pos = mcpart.Position(firstpoint).Vect();
1042  for(int i = firstpoint+1; i < ntrajpoints-1; ++i) {
1043  if(mcpart.Vx(i) >= tpcactiveXlow && mcpart.Vx(i) <= tpcactiveXhigh && mcpart.Vy(i) >= tpcactiveYlow && mcpart.Vy(i) <= tpcactiveYhigh && mcpart.Vz(i) >= tpcactiveZlow && mcpart.Vz(i) <= tpcactiveZhigh){
1044 
1045  pos -= mcpart.Position(i).Vect();
1046  plength += pos.Mag();
1047  pos = mcpart.Position(i).Vect();
1048  }
1049  }
1050  if(plength == 0.0) return 0.0;
1051 
1052  double length_laststep = length - plength;
1053  if(length_laststep <= 0.0) return 0.0;
1054 
1055  // Kinetic energy of the track
1056  double kinene_start = std::sqrt(mcpart.P(firstpoint)*mcpart.P(firstpoint) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
1057  double kinene_end = std::sqrt(mcpart.P(ntrajpoints-2)*mcpart.P(ntrajpoints-2) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
1058  double eneloss = (kinene_start - kinene_end);
1059  if(eneloss < 0.0 || kinene_start <= 0.0) return 0.0;
1060 
1061  // Average energy loss due to last trajectory point
1062  double aveneloss_laststep = length_laststep*eneloss/plength;
1063 
1064  return aveneloss_laststep;
1065 
1066 }
1067 
1068 // Get particle's kinetic energy at the interaction vertex
1069 double protoana::ProtoDUNETruthUtils::GetKinEnergyAtVertex(const simb::MCParticle& mcpart, double kinene_lastpoint){
1070 
1071  int ntrajpoints = mcpart.NumberTrajectoryPoints();
1072  double kinene_end = std::sqrt(mcpart.P(ntrajpoints-2)*mcpart.P(ntrajpoints-2) + mcpart.Mass()*mcpart.Mass()) - mcpart.Mass();
1073 
1074  return (kinene_end - kinene_lastpoint);
1075 }
1076 
1077 // Get the sim::IDEs from the MCParticle, organized by the trajectory points
1078 std::map< size_t, std::vector< const sim::IDE * > >
1080 {
1082  //art::ServiceHandle< geo::Geometry > geom;
1083 
1084  const simb::MCTrajectory & mctraj = mcpart.Trajectory();
1085 
1086  //Get all ides from the MCParticle
1087  std::vector< const sim::IDE * > ides = bt_serv->TrackIdToSimIDEs_Ps( mcpart.TrackId()/*, geom->View(2)*/ );
1088  //Sort by z position, assume traveling in beam direction.
1089  std::sort( ides.begin(), ides.end(), sort_IDEs );
1090 
1091  std::map< size_t, std::vector< const sim::IDE * > > results;
1092 
1093  size_t ide_start = 0;
1094 
1095  for( size_t i = 0; i < mctraj.size(); ++i ){
1096  results[i] = std::vector< const sim::IDE * >();
1097 
1098  for( size_t j = ide_start; j < ides.size(); ++j ){
1099 
1100  if( ( ides[j]->z > mctraj.Z(i) ) && ( ides[j]->z < mctraj.Z(i+1) ) ){
1101  results[i].push_back( ides[j] );
1102  }
1103 
1104  else if( ( ides[j]->z < mctraj.Z(i) ) && ( ides[j]->z < mctraj.Z(i+1) ) ){
1105  continue;
1106  }
1107 
1108  else{
1109  ide_start = j;
1110  break;
1111  }
1112  }
1113  }
1114 
1115  return results;
1116 }
1117 
1118 std::vector<const sim::IDE*>
1120  const TLorentzVector & p1,
1121  const TLorentzVector &p2)
1122 {
1125  std::vector<const sim::IDE *> ides =
1126  bt_serv->TrackIdToSimIDEs_Ps(mcpart.TrackId(), geom->View(2));
1127 
1128  std::vector<const sim::IDE *> results;
1129  for (const sim::IDE * theIDE : ides) {
1130  if (theIDE->z > p1.Z() && theIDE->z < p2.Z() /*&&
1131  theIDE->x > p1.X() && theIDE->x < p2.X() &&
1132  theIDE->y > p1.Y() && theIDE->y < p2.Y()*/) {
1133  results.push_back(theIDE);
1134  }
1135  }
1136 
1137  return results;
1138 }
1139 
1140 std::map< int, std::vector< int > > protoana::ProtoDUNETruthUtils::GetMapMCToPFPs_ByHits(detinfo::DetectorClocksData const& clockData,
1141  const art::Event & evt, std::string pfpTag, std::string hitTag ){
1142 
1143  std::map< int, std::vector< int > > results;
1144 
1145  auto allPFPs = evt.getValidHandle< std::vector< recob::PFParticle > >( pfpTag );
1146 
1147  for( size_t i = 0; i < allPFPs->size(); ++i ){
1148  auto thePFP = &(allPFPs->at(i));
1149  protoana::MCParticleSharedHits match = GetMCParticleByHits( clockData, *thePFP, evt, pfpTag, hitTag );
1150  if( match.particle ){
1151  results[ match.particle->TrackId() ].push_back( i );
1152  }
1153  }
1154 
1155  return results;
1156 }
1157 
1158 bool protoana::sort_IDEs( const sim::IDE * i1, const sim::IDE * i2){
1159  return( i1->z < i2->z );
1160 }
double E(const int i=0) const
Definition: MCParticle.h:233
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
double Z(const size_type i) const
Definition: MCTrajectory.h:151
const std::vector< const recob::Hit * > GetRecoShowerHits(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
Get the hits from a given reco shower.
float z
z position of ionization [cm]
Definition: SimChannel.h:117
int PdgCode() const
Definition: MCParticle.h:212
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromTrackHits(detinfo::DetectorClocksData const &clockData, const std::vector< const recob::Hit * > &hitVec, bool use_eve=true) const
QList< Entry > entry
double GetDepEnergyMC(const art::Event &evt, geo::GeometryCore const *fGeom, int trackid, int whichview) const
const simb::MCParticle * GetGeantGoodParticle(const simb::MCTruth &genTruth, const art::Event &evt) const
const recob::Shower * GetRecoShowerFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string showerModule) const
std::vector< art::Ptr< recob::Hit > > TrackIdToHits_Ps(detinfo::DetectorClocksData const &clockData, int tkId, std::vector< art::Ptr< recob::Hit >> const &hitsIn) const
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
std::string string
Definition: nybbler.cc:12
template const protoana::MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits< recob::PFParticle >(const detinfo::DetectorClocksData &, const recob::PFParticle &, const art::Event &, std::string, std::string) const
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
const simb::MCParticle * TrackIdToParticle_P(int id) const
double GetCompleteness(detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
double Mass() const
Definition: MCParticle.h:239
template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Shower >(detinfo::DetectorClocksData const &, const recob::Shower &, const art::Event &, std::string, std::string) const
const std::vector< const recob::Hit * > GetPFParticleHits(const recob::PFParticle &particle, art::Event const &evt, const std::string particleLabel) const
Get the hits associated to the PFParticle.
template const protoana::MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits< recob::Shower >(const detinfo::DetectorClocksData &, const recob::Shower &, const art::Event &, std::string, std::string) const
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
int PdgCode() const
Return the type of particle as a PDG ID.
Definition: PFParticle.h:83
const MCParticleSharedHits GetMCParticleByHits(detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
int NParticles() const
Definition: MCTruth.h:75
const float ConvertTrueTimeToPandoraTimeNano(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part) const
std::string Process() const
Definition: MCParticle.h:215
template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::Track >(detinfo::DetectorClocksData const &, const recob::Track &, const art::Event &, std::string, std::string) const
const recob::PFParticle * GetPFParticleFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string pfparticleModule) const
std::vector< int > HitToTrackIds(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
int TrackId() const
Definition: MCParticle.h:210
template std::vector< protoana::MCParticleSharedHits > protoana::ProtoDUNETruthUtils::GetMCParticleListByHits< recob::PFParticle >(detinfo::DetectorClocksData const &, const recob::PFParticle &, const art::Event &, std::string, std::string) const
double GetKinEnergyAtVertex(const simb::MCParticle &mcpart, double kinene_lastpoint=0.0)
double GetDepEnergyAtLastTrajPoint(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
double GetMCParticleLengthInTPCActiveVolume(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
std::vector< const recob::Hit * > GetMCParticleHits(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const art::Event &evt, std::string hitModule, bool use_eve=true) const
const simb::MCParticle * GetMCParticleFromReco(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
def process(f, kind)
Definition: search.py:254
bool isRealData() const
T abs(T value)
std::vector< MCParticleSharedHits > GetMCParticleListByHits(detinfo::DetectorClocksData const &clockData, const T &recobj, const art::Event &evt, std::string recoModule, std::string hitModule) const
const double e
std::vector< const recob::Hit * > FillSharedHits(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const std::vector< const recob::Hit * > &hitsVec, bool delta_ray, bool use_eve=true) const
std::map< int, std::vector< int > > GetMapMCToPFPs_ByHits(detinfo::DetectorClocksData const &clockData, const art::Event &evt, std::string pfpTag, std::string hitTag)
const recob::Track * GetRecoTrackFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string trackModule) const
double GetPurity(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule) const
int GetShowerIndex(const recob::Shower &shower, art::Event const &evt, const std::string showerModule) const
If the shower.ID() isn&#39;t filled we must find the actual shower index ourselves.
Ionization at a point of the TPC sensitive volume.
Definition: SimChannel.h:84
const double a
const simb::MCParticle * GetMCParticleFromPFParticle(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
double P(const int i=0) const
Definition: MCParticle.h:234
double T(const int i=0) const
Definition: MCParticle.h:224
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:441
p
Definition: test.py:223
const simb::MCParticle * MatchPduneMCtoG4(const simb::MCParticle &pDunePart, const art::Event &evt)
bool sort_IDEs(const sim::IDE *i1, const sim::IDE *i2)
const simb::MCParticle * GetMCParticleFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
Description of geometry of one entire detector.
std::vector< std::pair< const recob::Shower *, double > > GetRecoShowerListFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string showerModule) const
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
template const protoana::MCParticleSharedHits protoana::ProtoDUNETruthUtils::GetMCParticleByHits< recob::Track >(const detinfo::DetectorClocksData &, const recob::Track &, const art::Event &, std::string, std::string) const
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
std::map< size_t, std::vector< const sim::IDE * > > GetSimIDEs(const simb::MCParticle &mcpart)
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
std::vector< const sim::IDE * > GetSimIDEsBetweenPoints(const simb::MCParticle &mcpart, const TLorentzVector &p1, const TLorentzVector &p2)
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromRecoShower(detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
double Vx(const int i=0) const
Definition: MCParticle.h:221
int ID() const
Definition: Track.h:198
const std::vector< const recob::Hit * > GetRecoTrackHits(const recob::Track &track, art::Event const &evt, const std::string trackModule) const
Get the hits from a given reco track.
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromReco(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
Declaration of signal hit object.
const simb::MCParticle * particle
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
double G4ToElecTime(double const g4_time) const
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromPFParticle(detinfo::DetectorClocksData const &clockData, const recob::PFParticle &pfpart, art::Event const &evt, std::string pfparticleModule) const
size_type size() const
Definition: MCTrajectory.h:166
Contains all timing reference information for the detector.
const float ConvertTrueTimeToPandoraTimeMicro(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part) const
std::vector< sim::TrackIDE > HitToEveTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
double Vz(const int i=0) const
Definition: MCParticle.h:223
static bool * b
Definition: config.cpp:1043
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromShowerHits(detinfo::DetectorClocksData const &clockData, const std::vector< const recob::Hit * > &hitVec) const
cet::LibraryManager dummy("noplugin")
int GetFirstTrajectoryPointInTPCActiveVolume(const simb::MCParticle &mcpart, double tpcactiveXlow, double tpcactiveXhigh, double tpcactiveYlow, double tpcactiveYhigh, double tpcactiveZlow, double tpcactiveZhigh)
EventNumber_t event() const
Definition: EventID.h:116
std::vector< const sim::IDE * > HitToSimIDEs_Ps(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
std::vector< std::pair< const recob::PFParticle *, double > > GetPFParticleListFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string pfparticleModule, bool use_eve=true) const
std::vector< std::pair< const recob::Track *, double > > GetRecoTrackListFromMCParticle(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &part, art::Event const &evt, std::string trackModule) const
TCEvent evt
Definition: DataStructs.cxx:7
Event generator information.
Definition: MCTruth.h:32
const simb::MCParticle * GetMCParticleFromRecoShower(detinfo::DetectorClocksData const &clockData, const recob::Shower &shower, art::Event const &evt, std::string showerModule) const
std::vector< std::pair< const simb::MCParticle *, double > > GetMCParticleListFromRecoTrack(detinfo::DetectorClocksData const &clockData, const recob::Track &track, art::Event const &evt, std::string trackModule) const
Ionization energy from a Geant4 track.
Definition: SimChannel.h:25
std::vector< const recob::Hit * > GetSharedHits(detinfo::DetectorClocksData const &clockData, const simb::MCParticle &mcpart, const recob::PFParticle &pfpart, const art::Event &evt, std::string pfparticleModule, bool delta_ray=false) const
EventID id() const
Definition: Event.cc:34
double Vy(const int i=0) const
Definition: MCParticle.h:222
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
int GetProcessKey(std::string process)
QTextStream & endl(QTextStream &s)