LArPandoraHelper.cxx
Go to the documentation of this file.
1 /**
2  * @file larpandora/LArPandoraInterface/LArPandoraHelper.cxx
3  *
4  * @brief helper function for LArPandoraInterface producer module
5  *
6  */
7 
9 #include "canvas/Persistency/Common/FindManyP.h"
10 #include "canvas/Persistency/Common/FindOneP.h"
11 #include "cetlib_except/exception.h"
13 
30 
34 
35 #include "Objects/ParticleFlowObject.h"
36 #include "Pandora/PandoraInternal.h"
37 #include "Pandora/PdgTable.h"
38 
40 
41 #include <iostream>
42 #include <limits>
43 
44 namespace lar_pandora {
45 
46  void
48  const std::string& label,
49  WireVector& wireVector)
50  {
52  evt.getByLabel(label, theWires);
53 
54  if (!theWires.isValid()) {
55  mf::LogDebug("LArPandora") << " Failed to find wires... " << std::endl;
56  return;
57  }
58  else {
59  mf::LogDebug("LArPandora") << " Found: " << theWires->size() << " Wires " << std::endl;
60  }
61 
62  for (unsigned int i = 0; i < theWires->size(); ++i) {
63  const art::Ptr<recob::Wire> wire(theWires, i);
64  wireVector.push_back(wire);
65  }
66  }
67 
68  //------------------------------------------------------------------------------------------------------------------------------------------
69 
70  void
72  const std::string& label,
73  HitVector& hitVector)
74  {
76  evt.getByLabel(label, theHits);
77 
78  if (!theHits.isValid()) {
79  mf::LogDebug("LArPandora") << " Failed to find hits... " << std::endl;
80  return;
81  }
82  else {
83  mf::LogDebug("LArPandora") << " Found: " << theHits->size() << " Hits " << std::endl;
84  }
85 
86  for (unsigned int i = 0; i < theHits->size(); ++i) {
87  const art::Ptr<recob::Hit> hit(theHits, i);
88  hitVector.push_back(hit);
89  }
90  }
91 
92  //------------------------------------------------------------------------------------------------------------------------------------------
93 
94  void
96  const std::string& label,
97  PFParticleVector& particleVector)
98  {
100  evt.getByLabel(label, theParticles);
101 
102  if (!theParticles.isValid()) {
103  mf::LogDebug("LArPandora") << " Failed to find particles... " << std::endl;
104  return;
105  }
106  else {
107  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " PFParticles "
108  << std::endl;
109  }
110 
111  for (unsigned int i = 0; i < theParticles->size(); ++i) {
112  const art::Ptr<recob::PFParticle> particle(theParticles, i);
113  particleVector.push_back(particle);
114  }
115  }
116 
117  //------------------------------------------------------------------------------------------------------------------------------------------
118 
119  void
121  const std::string& label,
122  SpacePointVector& spacePointVector,
123  SpacePointsToHits& spacePointsToHits)
124  {
125  HitsToSpacePoints hitsToSpacePoints;
127  evt, label, spacePointVector, spacePointsToHits, hitsToSpacePoints);
128  }
129 
130  //------------------------------------------------------------------------------------------------------------------------------------------
131 
132  void
134  const std::string& label,
135  SpacePointVector& spacePointVector,
136  SpacePointsToHits& spacePointsToHits,
137  HitsToSpacePoints& hitsToSpacePoints)
138  {
140  evt.getByLabel(label, theSpacePoints);
141 
142  if (!theSpacePoints.isValid()) {
143  mf::LogDebug("LArPandora") << " Failed to find spacepoints... " << std::endl;
144  return;
145  }
146  else {
147  mf::LogDebug("LArPandora") << " Found: " << theSpacePoints->size() << " SpacePoints "
148  << std::endl;
149  }
150 
151  art::FindOneP<recob::Hit> theHitAssns(theSpacePoints, evt, label);
152  for (unsigned int i = 0; i < theSpacePoints->size(); ++i) {
153  const art::Ptr<recob::SpacePoint> spacepoint(theSpacePoints, i);
154  spacePointVector.push_back(spacepoint);
155  const art::Ptr<recob::Hit> hit = theHitAssns.at(i);
156  spacePointsToHits[spacepoint] = hit;
157  hitsToSpacePoints[hit] = spacepoint;
158  }
159  }
160 
161  //------------------------------------------------------------------------------------------------------------------------------------------
162 
163  void
165  const std::string& label,
166  ClusterVector& clusterVector,
167  ClustersToHits& clustersToHits)
168  {
170  evt.getByLabel(label, theClusters);
171 
172  if (!theClusters.isValid()) {
173  mf::LogDebug("LArPandora") << " Failed to find clusters... " << std::endl;
174  return;
175  }
176  else {
177  mf::LogDebug("LArPandora") << " Found: " << theClusters->size() << " Clusters " << std::endl;
178  }
179 
180  art::FindManyP<recob::Hit> theHitAssns(theClusters, evt, label);
181  for (unsigned int i = 0; i < theClusters->size(); ++i) {
182  const art::Ptr<recob::Cluster> cluster(theClusters, i);
183  clusterVector.push_back(cluster);
184 
185  const std::vector<art::Ptr<recob::Hit>> hits = theHitAssns.at(i);
186  for (unsigned int j = 0; j < hits.size(); ++j) {
187  const art::Ptr<recob::Hit> hit = hits.at(j);
188  clustersToHits[cluster].push_back(hit);
189  }
190  }
191  }
192 
193  //------------------------------------------------------------------------------------------------------------------------------------------
194 
195  void
197  const std::string& label,
198  PFParticleVector& particleVector,
199  PFParticlesToSpacePoints& particlesToSpacePoints)
200  {
202  evt.getByLabel(label, theParticles);
203 
204  if (!theParticles.isValid()) {
205  mf::LogDebug("LArPandora") << " Failed to find particles... " << std::endl;
206  return;
207  }
208  else {
209  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " PFParticles "
210  << std::endl;
211  }
212 
213  art::FindManyP<recob::SpacePoint> theSpacePointAssns(theParticles, evt, label);
214  for (unsigned int i = 0; i < theParticles->size(); ++i) {
215  const art::Ptr<recob::PFParticle> particle(theParticles, i);
216  particleVector.push_back(particle);
217 
218  const std::vector<art::Ptr<recob::SpacePoint>> spacepoints = theSpacePointAssns.at(i);
219  for (unsigned int j = 0; j < spacepoints.size(); ++j) {
220  const art::Ptr<recob::SpacePoint> spacepoint = spacepoints.at(j);
221  particlesToSpacePoints[particle].push_back(spacepoint);
222  }
223  }
224  }
225 
226  //------------------------------------------------------------------------------------------------------------------------------------------
227 
228  void
230  const std::string& label,
231  PFParticleVector& particleVector,
232  PFParticlesToClusters& particlesToClusters)
233  {
235  evt.getByLabel(label, theParticles);
236 
237  if (!theParticles.isValid()) {
238  mf::LogDebug("LArPandora") << " Failed to find particles... " << std::endl;
239  return;
240  }
241  else {
242  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " PFParticles "
243  << std::endl;
244  }
245 
246  art::FindManyP<recob::Cluster> theClusterAssns(theParticles, evt, label);
247  for (unsigned int i = 0; i < theParticles->size(); ++i) {
248  const art::Ptr<recob::PFParticle> particle(theParticles, i);
249  particleVector.push_back(particle);
250 
251  const std::vector<art::Ptr<recob::Cluster>> clusters = theClusterAssns.at(i);
252  for (unsigned int j = 0; j < clusters.size(); ++j) {
253  const art::Ptr<recob::Cluster> cluster = clusters.at(j);
254  particlesToClusters[particle].push_back(cluster);
255  }
256  }
257  }
258 
259  //------------------------------------------------------------------------------------------------------------------------------------------
260 
261  void
263  const std::string& label,
264  PFParticleVector& particleVector,
265  PFParticlesToMetadata& particlesToMetadata)
266  {
268  evt.getByLabel(label, theParticles);
269 
270  if (!theParticles.isValid()) {
271  mf::LogDebug("LArPandora") << " Failed to find particles... " << std::endl;
272  return;
273  }
274  else {
275  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " PFParticles "
276  << std::endl;
277  }
278 
279  art::FindManyP<larpandoraobj::PFParticleMetadata> theMetadataAssns(theParticles, evt, label);
280  for (unsigned int i = 0; i < theParticles->size(); ++i) {
281  const art::Ptr<recob::PFParticle> particle(theParticles, i);
282  particleVector.push_back(particle);
283 
284  const std::vector<art::Ptr<larpandoraobj::PFParticleMetadata>> pfParticleMetadataList =
285  theMetadataAssns.at(i);
286  for (unsigned int j = 0; j < pfParticleMetadataList.size(); ++j) {
287  const art::Ptr<larpandoraobj::PFParticleMetadata> pfParticleMetadata =
288  pfParticleMetadataList.at(j);
289  particlesToMetadata[particle].push_back(pfParticleMetadata);
290  }
291  }
292  }
293 
294  //------------------------------------------------------------------------------------------------------------------------------------------
295 
296  void
298  const std::string& label,
299  ShowerVector& showerVector,
300  PFParticlesToShowers& particlesToShowers)
301  {
303  evt.getByLabel(label, theShowers);
304 
305  if (!theShowers.isValid()) {
306  mf::LogDebug("LArPandora") << " Failed to find showers... " << std::endl;
307  return;
308  }
309  else {
310  mf::LogDebug("LArPandora") << " Found: " << theShowers->size() << " Showers " << std::endl;
311  }
312 
313  art::FindManyP<recob::PFParticle> theShowerAssns(theShowers, evt, label);
314  for (unsigned int i = 0; i < theShowers->size(); ++i) {
315  const art::Ptr<recob::Shower> shower(theShowers, i);
316  showerVector.push_back(shower);
317 
318  const std::vector<art::Ptr<recob::PFParticle>> particles = theShowerAssns.at(i);
319  for (unsigned int j = 0; j < particles.size(); ++j) {
320  const art::Ptr<recob::PFParticle> particle = particles.at(j);
321  particlesToShowers[particle].push_back(shower);
322  }
323  }
324  }
325 
326  //------------------------------------------------------------------------------------------------------------------------------------------
327 
328  void
330  const std::string& label,
331  TrackVector& trackVector,
332  PFParticlesToTracks& particlesToTracks)
333  {
335  evt.getByLabel(label, theTracks);
336 
337  if (!theTracks.isValid()) {
338  mf::LogDebug("LArPandora") << " Failed to find tracks... " << std::endl;
339  return;
340  }
341  else {
342  mf::LogDebug("LArPandora") << " Found: " << theTracks->size() << " Tracks " << std::endl;
343  }
344 
345  art::FindManyP<recob::PFParticle> theTrackAssns(theTracks, evt, label);
346  for (unsigned int i = 0; i < theTracks->size(); ++i) {
347  const art::Ptr<recob::Track> track(theTracks, i);
348  trackVector.push_back(track);
349 
350  const std::vector<art::Ptr<recob::PFParticle>> particles = theTrackAssns.at(i);
351  for (unsigned int j = 0; j < particles.size(); ++j) {
352  const art::Ptr<recob::PFParticle> particle = particles.at(j);
353  particlesToTracks[particle].push_back(track);
354  }
355  }
356  }
357 
358  //------------------------------------------------------------------------------------------------------------------------------------------
359 
360  void
362  const std::string& label,
363  TrackVector& trackVector,
364  TracksToHits& tracksToHits)
365  {
367  evt.getByLabel(label, theTracks);
368 
369  if (!theTracks.isValid()) {
370  mf::LogDebug("LArPandora") << " Failed to find tracks... " << std::endl;
371  return;
372  }
373  else {
374  mf::LogDebug("LArPandora") << " Found: " << theTracks->size() << " Tracks " << std::endl;
375  }
376 
377  art::FindManyP<recob::Hit> theHitAssns(theTracks, evt, label);
378  for (unsigned int i = 0; i < theTracks->size(); ++i) {
379  const art::Ptr<recob::Track> track(theTracks, i);
380  trackVector.push_back(track);
381 
382  const std::vector<art::Ptr<recob::Hit>> hits = theHitAssns.at(i);
383  for (unsigned int j = 0; j < hits.size(); ++j) {
384  const art::Ptr<recob::Hit> hit = hits.at(j);
385  tracksToHits[track].push_back(hit);
386  }
387  }
388  }
389 
390  //------------------------------------------------------------------------------------------------------------------------------------------
391 
392  void
394  const std::string& label,
395  ShowerVector& showerVector,
396  ShowersToHits& showersToHits)
397  {
399  evt.getByLabel(label, theShowers);
400 
401  if (!theShowers.isValid()) {
402  mf::LogDebug("LArPandora") << " Failed to find showers... " << std::endl;
403  return;
404  }
405  else {
406  mf::LogDebug("LArPandora") << " Found: " << theShowers->size() << " Showers " << std::endl;
407  }
408 
409  art::FindManyP<recob::Hit> theHitAssns(theShowers, evt, label);
410  for (unsigned int i = 0; i < theShowers->size(); ++i) {
411  const art::Ptr<recob::Shower> shower(theShowers, i);
412  showerVector.push_back(shower);
413 
414  const std::vector<art::Ptr<recob::Hit>> hits = theHitAssns.at(i);
415  for (unsigned int j = 0; j < hits.size(); ++j) {
416  const art::Ptr<recob::Hit> hit = hits.at(j);
417  showersToHits[shower].push_back(hit);
418  }
419  }
420  }
421 
422  //------------------------------------------------------------------------------------------------------------------------------------------
423 
424  void
426  const std::string& label,
427  SeedVector& seedVector,
428  PFParticlesToSeeds& particlesToSeeds)
429  {
431  evt.getByLabel(label, theSeeds);
432 
433  if (!theSeeds.isValid()) {
434  mf::LogDebug("LArPandora") << " Failed to find seeds... " << std::endl;
435  return;
436  }
437  else {
438  mf::LogDebug("LArPandora") << " Found: " << theSeeds->size() << " Seeds " << std::endl;
439  }
440 
441  art::FindManyP<recob::PFParticle> theSeedAssns(theSeeds, evt, label);
442  for (unsigned int i = 0; i < theSeeds->size(); ++i) {
443  const art::Ptr<recob::Seed> seed(theSeeds, i);
444  seedVector.push_back(seed);
445 
446  const std::vector<art::Ptr<recob::PFParticle>> particles = theSeedAssns.at(i);
447  for (unsigned int j = 0; j < particles.size(); ++j) {
448  const art::Ptr<recob::PFParticle> particle = particles.at(j);
449  particlesToSeeds[particle].push_back(seed);
450  }
451  }
452  }
453 
454  //------------------------------------------------------------------------------------------------------------------------------------------
455 
456  void
458  const std::string& label,
459  SeedVector& seedVector,
460  SeedsToHits& seedsToHits)
461  {
463  evt.getByLabel(label, theSeeds);
464 
465  if (!theSeeds.isValid()) {
466  mf::LogDebug("LArPandora") << " Failed to find seeds... " << std::endl;
467  return;
468  }
469  else {
470  mf::LogDebug("LArPandora") << " Found: " << theSeeds->size() << " Seeds " << std::endl;
471  }
472 
473  art::FindOneP<recob::Hit> theHitAssns(theSeeds, evt, label);
474 
475  if (!theHitAssns.isValid()) {
476  mf::LogDebug("LArPandora") << " Failed to find seed associations... " << std::endl;
477  return;
478  }
479 
480  for (unsigned int i = 0; i < theSeeds->size(); ++i) {
481  const art::Ptr<recob::Seed> seed(theSeeds, i);
482  seedVector.push_back(seed);
483  const art::Ptr<recob::Hit> hit = theHitAssns.at(i);
484  seedsToHits[seed] = hit;
485  }
486  }
487 
488  //------------------------------------------------------------------------------------------------------------------------------------------
489 
490  void
492  const std::string& label,
493  VertexVector& vertexVector,
494  PFParticlesToVertices& particlesToVertices)
495  {
497  evt.getByLabel(label, theVertices);
498 
499  if (!theVertices.isValid()) {
500  mf::LogDebug("LArPandora") << " Failed to find vertices... " << std::endl;
501  return;
502  }
503  else {
504  mf::LogDebug("LArPandora") << " Found: " << theVertices->size() << " Vertices " << std::endl;
505  }
506 
507  art::FindManyP<recob::PFParticle> theVerticesAssns(theVertices, evt, label);
508  for (unsigned int i = 0; i < theVertices->size(); ++i) {
509  const art::Ptr<recob::Vertex> vertex(theVertices, i);
510  vertexVector.push_back(vertex);
511 
512  const std::vector<art::Ptr<recob::PFParticle>> particles = theVerticesAssns.at(i);
513  for (unsigned int j = 0; j < particles.size(); ++j) {
514  const art::Ptr<recob::PFParticle> particle = particles.at(j);
515  particlesToVertices[particle].push_back(vertex);
516  }
517  }
518  }
519 
520  //------------------------------------------------------------------------------------------------------------------------------------------
521 
522  void
524  const PFParticlesToSpacePoints& particlesToSpacePoints,
525  const SpacePointsToHits& spacePointsToHits,
526  PFParticlesToHits& particlesToHits,
527  HitsToPFParticles& hitsToParticles,
528  const DaughterMode daughterMode)
529  {
530  // Build mapping from particle to particle ID for parent/daughter navigation
531  PFParticleMap particleMap;
532 
533  for (PFParticleVector::const_iterator iter1 = particleVector.begin(),
534  iterEnd1 = particleVector.end();
535  iter1 != iterEnd1;
536  ++iter1) {
537  const art::Ptr<recob::PFParticle> particle = *iter1;
538  particleMap[particle->Self()] = particle;
539  }
540 
541  // Loop over hits and build mapping between reconstructed final-state particles and reconstructed hits
542  for (PFParticlesToSpacePoints::const_iterator iter1 = particlesToSpacePoints.begin(),
543  iterEnd1 = particlesToSpacePoints.end();
544  iter1 != iterEnd1;
545  ++iter1) {
546  const art::Ptr<recob::PFParticle> thisParticle = iter1->first;
547  const art::Ptr<recob::PFParticle> particle(
548  (kAddDaughters == daughterMode) ?
549  LArPandoraHelper::GetFinalStatePFParticle(particleMap, thisParticle) :
550  thisParticle);
551 
552  if ((kIgnoreDaughters == daughterMode) &&
553  !LArPandoraHelper::IsFinalState(particleMap, particle))
554  continue;
555 
556  const SpacePointVector& spacePointVector = iter1->second;
557 
558  for (SpacePointVector::const_iterator iter2 = spacePointVector.begin(),
559  iterEnd2 = spacePointVector.end();
560  iter2 != iterEnd2;
561  ++iter2) {
562  const art::Ptr<recob::SpacePoint> spacepoint = *iter2;
563 
564  SpacePointsToHits::const_iterator iter3 = spacePointsToHits.find(spacepoint);
565  if (spacePointsToHits.end() == iter3)
566  throw cet::exception("LArPandora") << " PandoraCollector::BuildPFParticleHitMaps --- "
567  "Found a space point without an associated hit ";
568 
569  const art::Ptr<recob::Hit> hit = iter3->second;
570 
571  particlesToHits[particle].push_back(hit);
572  hitsToParticles[hit] = particle;
573  }
574  }
575  }
576 
577  //------------------------------------------------------------------------------------------------------------------------------------------
578 
579  void
581  const PFParticlesToClusters& particlesToClusters,
582  const ClustersToHits& clustersToHits,
583  PFParticlesToHits& particlesToHits,
584  HitsToPFParticles& hitsToParticles,
585  const DaughterMode daughterMode)
586  {
587  // Build mapping from particle to particle ID for parent/daughter navigation
588  PFParticleMap particleMap;
589 
590  for (PFParticleVector::const_iterator iter1 = particleVector.begin(),
591  iterEnd1 = particleVector.end();
592  iter1 != iterEnd1;
593  ++iter1) {
594  const art::Ptr<recob::PFParticle> particle = *iter1;
595  particleMap[particle->Self()] = particle;
596  }
597 
598  // Loop over hits and build mapping between reconstructed final-state particles and reconstructed hits
599  for (PFParticlesToClusters::const_iterator iter1 = particlesToClusters.begin(),
600  iterEnd1 = particlesToClusters.end();
601  iter1 != iterEnd1;
602  ++iter1) {
603  const art::Ptr<recob::PFParticle> thisParticle = iter1->first;
604  const art::Ptr<recob::PFParticle> particle(
605  (kAddDaughters == daughterMode) ?
606  LArPandoraHelper::GetFinalStatePFParticle(particleMap, thisParticle) :
607  thisParticle);
608 
609  if ((kIgnoreDaughters == daughterMode) &&
610  !LArPandoraHelper::IsFinalState(particleMap, particle))
611  continue;
612 
613  const ClusterVector& clusterVector = iter1->second;
614  for (ClusterVector::const_iterator iter2 = clusterVector.begin(),
615  iterEnd2 = clusterVector.end();
616  iter2 != iterEnd2;
617  ++iter2) {
618  const art::Ptr<recob::Cluster> cluster = *iter2;
619 
620  ClustersToHits::const_iterator iter3 = clustersToHits.find(cluster);
621  if (clustersToHits.end() == iter3)
622  throw cet::exception("LArPandora") << " PandoraCollector::BuildPFParticleHitMaps --- "
623  "Found a space point without an associated hit ";
624 
625  const HitVector& hitVector = iter3->second;
626  for (HitVector::const_iterator iter4 = hitVector.begin(), iterEnd4 = hitVector.end();
627  iter4 != iterEnd4;
628  ++iter4) {
629  const art::Ptr<recob::Hit> hit = *iter4;
630 
631  particlesToHits[particle].push_back(hit);
632  hitsToParticles[hit] = particle;
633  }
634  }
635  }
636  }
637 
638  //------------------------------------------------------------------------------------------------------------------------------------------
639 
640  void
642  const std::string& label,
643  PFParticlesToHits& particlesToHits,
644  HitsToPFParticles& hitsToParticles,
645  const DaughterMode daughterMode,
646  const bool useClusters)
647  {
649  evt, label, label, particlesToHits, hitsToParticles, daughterMode, useClusters);
650  }
651 
652  //------------------------------------------------------------------------------------------------------------------------------------------
653 
654  void
656  const std::string& label_pfpart,
657  const std::string& label_middle,
658  PFParticlesToHits& particlesToHits,
659  HitsToPFParticles& hitsToParticles,
660  const DaughterMode daughterMode,
661  const bool useClusters)
662  {
663  // Use intermediate clusters
664  if (useClusters) {
665  PFParticleVector particleVector;
666  PFParticlesToClusters particlesToClusters;
667 
668  ClusterVector clusterVector;
669  ClustersToHits clustersToHits;
670 
671  LArPandoraHelper::CollectPFParticles(evt, label_pfpart, particleVector, particlesToClusters);
672  LArPandoraHelper::CollectClusters(evt, label_middle, clusterVector, clustersToHits);
673 
675  particlesToClusters,
676  clustersToHits,
677  particlesToHits,
678  hitsToParticles,
679  daughterMode);
680  }
681 
682  // Use intermediate space points
683  else {
684  PFParticleVector particleVector;
685  PFParticlesToSpacePoints particlesToSpacePoints;
686 
687  SpacePointVector spacePointVector;
688  SpacePointsToHits spacePointsToHits;
689 
691  evt, label_pfpart, particleVector, particlesToSpacePoints);
692  LArPandoraHelper::CollectSpacePoints(evt, label_middle, spacePointVector, spacePointsToHits);
693 
695  particlesToSpacePoints,
696  spacePointsToHits,
697  particlesToHits,
698  hitsToParticles,
699  daughterMode);
700  }
701  }
702 
703  //------------------------------------------------------------------------------------------------------------------------------------------
704 
705  void
707  PFParticleVector& outputParticles)
708  {
709  for (PFParticleVector::const_iterator iter = inputParticles.begin(),
710  iterEnd = inputParticles.end();
711  iter != iterEnd;
712  ++iter) {
713  const art::Ptr<recob::PFParticle> particle = *iter;
714 
715  if (LArPandoraHelper::IsNeutrino(particle)) outputParticles.push_back(particle);
716  }
717  }
718 
719  //------------------------------------------------------------------------------------------------------------------------------------------
720 
721  void
723  PFParticleVector& outputParticles)
724  {
725  // Build mapping from particle to particle ID for parent/daughter navigation
726  PFParticleMap particleMap;
727 
728  for (PFParticleVector::const_iterator iter = inputParticles.begin(),
729  iterEnd = inputParticles.end();
730  iter != iterEnd;
731  ++iter) {
732  const art::Ptr<recob::PFParticle> particle = *iter;
733  particleMap[particle->Self()] = particle;
734  }
735 
736  // Select final-state particles
737  for (PFParticleVector::const_iterator iter = inputParticles.begin(),
738  iterEnd = inputParticles.end();
739  iter != iterEnd;
740  ++iter) {
741  const art::Ptr<recob::PFParticle> particle = *iter;
742 
743  if (LArPandoraHelper::IsFinalState(particleMap, particle))
744  outputParticles.push_back(particle);
745  }
746  }
747 
748  //------------------------------------------------------------------------------------------------------------------------------------------
749 
750  void
752  const std::string& label,
753  CosmicTagVector& cosmicTagVector,
754  TracksToCosmicTags& tracksToCosmicTags)
755  {
757  evt.getByLabel(label, theCosmicTags);
758 
759  if (theCosmicTags.isValid()) {
760  art::FindOneP<recob::Track> theCosmicAssns(
761  theCosmicTags, evt, label); // We assume there is one tag per algorithm
762  for (unsigned int i = 0; i < theCosmicTags->size(); ++i) {
763  const art::Ptr<anab::CosmicTag> cosmicTag(theCosmicTags, i);
764  const art::Ptr<recob::Track> track = theCosmicAssns.at(i);
765  tracksToCosmicTags[track].push_back(
766  cosmicTag); // We assume there could be multiple algorithms
767  cosmicTagVector.push_back(cosmicTag);
768  }
769  }
770  }
771 
772  //------------------------------------------------------------------------------------------------------------------------------------------
773 
774  void
776  const std::string& label,
777  T0Vector& t0Vector,
778  PFParticlesToT0s& particlesToT0s)
779  {
781  evt.getByLabel(label, theT0s);
782 
783  if (theT0s.isValid()) {
784  art::FindManyP<recob::PFParticle> theAssns(theT0s, evt, label);
785  for (unsigned int i = 0; i < theT0s->size(); ++i) {
786  const art::Ptr<anab::T0> theT0(theT0s, i);
787  t0Vector.push_back(theT0);
788 
789  const std::vector<art::Ptr<recob::PFParticle>> particles = theAssns.at(i);
790  for (unsigned int j = 0; j < particles.size(); ++j) {
791  const art::Ptr<recob::PFParticle> theParticle = particles.at(j);
792  particlesToT0s[theParticle].push_back(
793  theT0); // We assume there could be multiple T0s per PFParticle
794  }
795  }
796  }
797  }
798 
799  //------------------------------------------------------------------------------------------------------------------------------------------
800 
801  void
803  const std::string& label,
804  SimChannelVector& simChannelVector,
805  bool& areSimChannelsValid)
806  {
808  evt.getByLabel(label, theSimChannels);
809 
810  if (!theSimChannels.isValid()) {
811  mf::LogDebug("LArPandora") << " Failed to find sim channels... " << std::endl;
812  areSimChannelsValid = false;
813  return;
814  }
815  else {
816  mf::LogDebug("LArPandora") << " Found: " << theSimChannels->size() << " SimChannels "
817  << std::endl;
818  areSimChannelsValid = true;
819  }
820 
821  for (unsigned int i = 0; i < theSimChannels->size(); ++i) {
822  const art::Ptr<sim::SimChannel> channel(theSimChannels, i);
823  simChannelVector.push_back(channel);
824  }
825  }
826 
827  //------------------------------------------------------------------------------------------------------------------------------------------
828 
829  void
831  const std::string& label,
832  MCParticleVector& particleVector)
833  {
835  evt.getByLabel(label, theParticles);
836 
837  if (!theParticles.isValid()) {
838  mf::LogDebug("LArPandora") << " Failed to find MC particles... " << std::endl;
839  return;
840  }
841  else {
842  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " MC particles "
843  << std::endl;
844  }
845 
846  for (unsigned int i = 0; i < theParticles->size(); ++i) {
847  const art::Ptr<simb::MCParticle> particle(theParticles, i);
848  particleVector.push_back(particle);
849  }
850  }
851 
852  //------------------------------------------------------------------------------------------------------------------------------------------
853 
854  void
856  const std::string& label,
857  RawMCParticleVector& particleVector)
858  {
860  evt.getByLabel(label, mcTruthBlocks);
861 
862  if (!mcTruthBlocks.isValid()) {
863  mf::LogDebug("LArPandora") << " Failed to find MC truth blocks from generator... "
864  << std::endl;
865  return;
866  }
867  else {
868  mf::LogDebug("LArPandora") << " Found: " << mcTruthBlocks->size() << " MC truth blocks "
869  << std::endl;
870  }
871 
872  if (mcTruthBlocks->size() != 1)
873  throw cet::exception("LArPandora") << " PandoraCollector::CollectGeneratorMCParticles --- "
874  "Unexpected number of MC truth blocks ";
875 
876  const art::Ptr<simb::MCTruth> mcTruth(mcTruthBlocks, 0);
877 
878  for (int i = 0; i < mcTruth->NParticles(); ++i) {
879  particleVector.push_back(mcTruth->GetParticle(i));
880  }
881  }
882 
883  //------------------------------------------------------------------------------------------------------------------------------------------
884 
885  void
887  const std::string& label,
888  MCTruthToMCParticles& truthToParticles,
889  MCParticlesToMCTruth& particlesToTruth)
890  {
892  evt.getByLabel(label, theParticles);
893 
894  if (!theParticles.isValid()) {
895  mf::LogDebug("LArPandora") << " Failed to find MC particles... " << std::endl;
896  return;
897  }
898  else {
899  mf::LogDebug("LArPandora") << " Found: " << theParticles->size() << " MC particles "
900  << std::endl;
901  }
902 
903  art::FindOneP<simb::MCTruth> theTruthAssns(theParticles, evt, label);
904 
905  for (unsigned int i = 0, iEnd = theParticles->size(); i < iEnd; ++i) {
906  const art::Ptr<simb::MCParticle> particle(theParticles, i);
907  const art::Ptr<simb::MCTruth> truth(theTruthAssns.at(i));
908  truthToParticles[truth].push_back(particle);
909  particlesToTruth[particle] = truth;
910  }
911  }
912 
913  //------------------------------------------------------------------------------------------------------------------------------------------
914 
915  void
917  const HitVector& hitVector,
918  const SimChannelVector& simChannelVector,
919  HitsToTrackIDEs& hitsToTrackIDEs)
920  {
921  auto const clock_data =
923 
924  SimChannelMap simChannelMap;
925 
926  for (SimChannelVector::const_iterator iter = simChannelVector.begin(),
927  iterEnd = simChannelVector.end();
928  iter != iterEnd;
929  ++iter) {
930  const art::Ptr<sim::SimChannel> simChannel = *iter;
931  simChannelMap.insert(SimChannelMap::value_type(simChannel->Channel(), simChannel));
932  }
933 
934  for (HitVector::const_iterator iter = hitVector.begin(), iterEnd = hitVector.end();
935  iter != iterEnd;
936  ++iter) {
937  const art::Ptr<recob::Hit> hit = *iter;
938 
939  SimChannelMap::const_iterator sIter = simChannelMap.find(hit->Channel());
940  if (simChannelMap.end() == sIter) continue; // Hit has no truth information [continue]
941 
942  // ATTN: Need to convert TDCtick (integer) to TDC (unsigned integer) before passing to simChannel
943  const raw::TDCtick_t start_tick(clock_data.TPCTick2TDC(hit->PeakTimeMinusRMS()));
944  const raw::TDCtick_t end_tick(clock_data.TPCTick2TDC(hit->PeakTimePlusRMS()));
945  const unsigned int start_tdc((start_tick < 0) ? 0 : start_tick);
946  const unsigned int end_tdc(end_tick);
947 
948  if (start_tdc > end_tdc) continue; // Hit undershoots the readout window [continue]
949 
950  const art::Ptr<sim::SimChannel> simChannel = sIter->second;
951  const TrackIDEVector trackCollection(simChannel->TrackIDEs(start_tdc, end_tdc));
952 
953  if (trackCollection.empty()) continue; // Hit has no truth information [continue]
954 
955  for (unsigned int iTrack = 0, iTrackEnd = trackCollection.size(); iTrack < iTrackEnd;
956  ++iTrack) {
957  const sim::TrackIDE trackIDE = trackCollection.at(iTrack);
958  hitsToTrackIDEs[hit].push_back(trackIDE);
959  }
960  }
961  }
962 
963  //------------------------------------------------------------------------------------------------------------------------------------------
964 
965  void
967  const MCTruthToMCParticles& truthToParticles,
968  MCParticlesToHits& particlesToHits,
969  HitsToMCParticles& hitsToParticles,
970  const DaughterMode daughterMode)
971  {
972  // Build mapping between particles and track IDs for parent/daughter navigation
973  MCParticleMap particleMap;
974 
975  for (MCTruthToMCParticles::const_iterator iter1 = truthToParticles.begin(),
976  iterEnd1 = truthToParticles.end();
977  iter1 != iterEnd1;
978  ++iter1) {
979  const MCParticleVector& particleVector = iter1->second;
980  for (MCParticleVector::const_iterator iter2 = particleVector.begin(),
981  iterEnd2 = particleVector.end();
982  iter2 != iterEnd2;
983  ++iter2) {
984  const art::Ptr<simb::MCParticle> particle = *iter2;
985  particleMap[particle->TrackId()] = particle;
986  }
987  }
988 
989  // Loop over hits and build mapping between reconstructed hits and true particles
990  for (HitsToTrackIDEs::const_iterator iter1 = hitsToTrackIDEs.begin(),
991  iterEnd1 = hitsToTrackIDEs.end();
992  iter1 != iterEnd1;
993  ++iter1) {
994  const art::Ptr<recob::Hit> hit = iter1->first;
995  const TrackIDEVector& trackCollection = iter1->second;
996 
997  int bestTrackID(-1);
998  float bestEnergyFrac(0.f);
999 
1000  for (TrackIDEVector::const_iterator iter2 = trackCollection.begin(),
1001  iterEnd2 = trackCollection.end();
1002  iter2 != iterEnd2;
1003  ++iter2) {
1004  const sim::TrackIDE& trackIDE = *iter2;
1005  const int trackID(std::abs(trackIDE.trackID)); // TODO: Find out why std::abs is needed
1006  const float energyFrac(trackIDE.energyFrac);
1007 
1008  if (energyFrac > bestEnergyFrac) {
1009  bestEnergyFrac = energyFrac;
1010  bestTrackID = trackID;
1011  }
1012  }
1013 
1014  if (bestTrackID >= 0) {
1015  MCParticleMap::const_iterator iter3 = particleMap.find(bestTrackID);
1016  if (particleMap.end() == iter3)
1017  throw cet::exception("LArPandora") << " PandoraCollector::BuildMCParticleHitMaps --- "
1018  "Found a track ID without an MC Particle ";
1019 
1020  try {
1021  const art::Ptr<simb::MCParticle> thisParticle = iter3->second;
1022  const art::Ptr<simb::MCParticle> primaryParticle(
1023  LArPandoraHelper::GetFinalStateMCParticle(particleMap, thisParticle));
1024  const art::Ptr<simb::MCParticle> selectedParticle(
1025  (kAddDaughters == daughterMode) ? primaryParticle : thisParticle);
1026 
1027  if ((kIgnoreDaughters == daughterMode) && (selectedParticle != primaryParticle)) continue;
1028 
1029  if (!(LArPandoraHelper::IsVisible(selectedParticle))) continue;
1030 
1031  particlesToHits[selectedParticle].push_back(hit);
1032  hitsToParticles[hit] = selectedParticle;
1033  }
1034  catch (cet::exception& e) {
1035  }
1036  }
1037  }
1038  }
1039 
1040  //------------------------------------------------------------------------------------------------------------------------------------------
1041 
1042  void
1044  const std::string& label,
1045  const HitVector& hitVector,
1046  MCParticlesToHits& particlesToHits,
1047  HitsToMCParticles& hitsToParticles,
1048  const DaughterMode daughterMode)
1049  {
1050  SimChannelVector simChannelVector;
1051  MCTruthToMCParticles truthToParticles;
1052  MCParticlesToMCTruth particlesToTruth;
1053  HitsToTrackIDEs hitsToTrackIDEs;
1054 
1055  bool areSimChannelsValid(false);
1056  LArPandoraHelper::CollectSimChannels(evt, label, simChannelVector, areSimChannelsValid);
1057 
1058  LArPandoraHelper::CollectMCParticles(evt, label, truthToParticles, particlesToTruth);
1059  LArPandoraHelper::BuildMCParticleHitMaps(evt, hitVector, simChannelVector, hitsToTrackIDEs);
1061  hitsToTrackIDEs, truthToParticles, particlesToHits, hitsToParticles, daughterMode);
1062  }
1063 
1064  //------------------------------------------------------------------------------------------------------------------------------------------
1065 
1066  void
1068  const std::string& hitLabel,
1069  const std::string& backtrackLabel,
1070  HitsToTrackIDEs& hitsToTrackIDEs)
1071  {
1072  // Start by getting the collection of Hits
1074  evt.getByLabel(hitLabel, theHits);
1075 
1076  if (!theHits.isValid()) {
1077  mf::LogDebug("LArPandora") << " Failed to find hits... " << std::endl;
1078  return;
1079  }
1080 
1081  HitVector hitVector;
1082 
1083  for (unsigned int i = 0; i < theHits->size(); ++i) {
1084  const art::Ptr<recob::Hit> hit(theHits, i);
1085  hitVector.push_back(hit);
1086  }
1087 
1088  // Now get the associations between Hits and MCParticles
1089  std::vector<anab::BackTrackerHitMatchingData const*> backtrackerVector;
1090 
1091  MCParticleVector particleVector;
1092 
1093  art::FindManyP<simb::MCParticle, anab::BackTrackerHitMatchingData> particles_per_hit(
1094  theHits, evt, backtrackLabel);
1095 
1096  if (!particles_per_hit.isValid()) {
1097  mf::LogDebug("LArPandora") << " Failed to find reco-truth matching... " << std::endl;
1098  return;
1099  }
1100 
1101  // Now loop over the hits and build a collection of IDEs
1102  for (HitVector::const_iterator iter = hitVector.begin(), iterEnd = hitVector.end();
1103  iter != iterEnd;
1104  ++iter) {
1105  const art::Ptr<recob::Hit> hit = *iter;
1106 
1107  particleVector.clear();
1108  backtrackerVector.clear();
1109  particles_per_hit.get(hit.key(), particleVector, backtrackerVector);
1110 
1111  for (unsigned int j = 0; j < particleVector.size(); ++j) {
1112  const art::Ptr<simb::MCParticle> particle = particleVector[j];
1113 
1114  sim::TrackIDE trackIDE;
1115  trackIDE.trackID = particle->TrackId();
1116  trackIDE.energy = backtrackerVector[j]->energy;
1117  trackIDE.energyFrac = backtrackerVector[j]->ideFraction;
1118 
1119  hitsToTrackIDEs[hit].push_back(trackIDE);
1120  }
1121  }
1122  }
1123 
1124  //------------------------------------------------------------------------------------------------------------------------------------------
1125 
1126  void
1128  const std::string& truthLabel,
1129  const std::string& hitLabel,
1130  const std::string& backtrackLabel,
1131  MCParticlesToHits& particlesToHits,
1132  HitsToMCParticles& hitsToParticles,
1133  const DaughterMode daughterMode)
1134  {
1135  MCTruthToMCParticles truthToParticles;
1136  MCParticlesToMCTruth particlesToTruth;
1137  HitsToTrackIDEs hitsToTrackIDEs;
1138 
1139  LArPandoraHelper::CollectMCParticles(evt, truthLabel, truthToParticles, particlesToTruth);
1140  LArPandoraHelper::BuildMCParticleHitMaps(evt, hitLabel, backtrackLabel, hitsToTrackIDEs);
1142  hitsToTrackIDEs, truthToParticles, particlesToHits, hitsToParticles, daughterMode);
1143  }
1144 
1145  //------------------------------------------------------------------------------------------------------------------------------------------
1146 
1147  template <typename T>
1148  void
1150  const std::string& label,
1151  const std::vector<art::Ptr<T>>& inputVector,
1152  HitVector& associatedHits,
1153  const pandora::IntVector* const indexVector)
1154  {
1155 
1157  evt.getByLabel(label, handle);
1158  art::FindManyP<recob::Hit> hitAssoc(handle, evt, label);
1159 
1160  if (indexVector != nullptr) {
1161  if (inputVector.size() != indexVector->size())
1162  throw cet::exception("LArPandora") << " PandoraHelper::GetAssociatedHits --- trying to use "
1163  "an index vector not matching input vector";
1164 
1165  // If indexVector is filled, sort hits according to trajectory points order
1166  for (int index : (*indexVector)) {
1167  const art::Ptr<T>& element = inputVector.at(index);
1168  const HitVector& hits = hitAssoc.at(element.key());
1169  associatedHits.insert(associatedHits.end(), hits.begin(), hits.end());
1170  }
1171  }
1172  else {
1173  // If indexVector is empty just loop through inputSpacePoints
1174  for (const art::Ptr<T>& element : inputVector) {
1175  const HitVector& hits = hitAssoc.at(element.key());
1176  associatedHits.insert(associatedHits.end(), hits.begin(), hits.end());
1177  }
1178  }
1179  }
1180 
1181  //------------------------------------------------------------------------------------------------------------------------------------------
1182 
1183  void
1185  MCParticleMap& particleMap)
1186  {
1187  for (MCParticleVector::const_iterator iter = particleVector.begin(),
1188  iterEnd = particleVector.end();
1189  iter != iterEnd;
1190  ++iter) {
1191  const art::Ptr<simb::MCParticle> particle = *iter;
1192  particleMap[particle->TrackId()] = particle;
1193  particleMap[particle->TrackId()] = particle;
1194  }
1195  }
1196 
1197  //------------------------------------------------------------------------------------------------------------------------------------------
1198 
1199  void
1201  PFParticleMap& particleMap)
1202  {
1203  for (PFParticleVector::const_iterator iter = particleVector.begin(),
1204  iterEnd = particleVector.end();
1205  iter != iterEnd;
1206  ++iter) {
1207  const art::Ptr<recob::PFParticle> particle = *iter;
1208  particleMap[particle->Self()] = particle;
1209  }
1210  }
1211 
1212  //------------------------------------------------------------------------------------------------------------------------------------------
1213 
1216  const art::Ptr<recob::PFParticle> inputParticle)
1217  {
1218  // Navigate upward through PFO daughter/parent links - return the top-level PF Particle
1219  int primaryTrackID(inputParticle->Self());
1220 
1221  if (!inputParticle->IsPrimary()) {
1222  while (1) {
1223  PFParticleMap::const_iterator pIter1 = particleMap.find(primaryTrackID);
1224  if (particleMap.end() == pIter1)
1225  throw cet::exception("LArPandora") << " PandoraCollector::GetParentPFParticle --- Found "
1226  "a PFParticle without a particle ID ";
1227 
1228  const art::Ptr<recob::PFParticle> primaryParticle = pIter1->second;
1229  if (primaryParticle->IsPrimary()) break;
1230 
1231  primaryTrackID = primaryParticle->Parent();
1232  }
1233  }
1234 
1235  PFParticleMap::const_iterator pIter2 = particleMap.find(primaryTrackID);
1236  if (particleMap.end() == pIter2)
1237  throw cet::exception("LArPandora")
1238  << " PandoraCollector::GetParentPFParticle --- Found a PFParticle without a particle ID ";
1239 
1240  const art::Ptr<recob::PFParticle> outputParticle = pIter2->second;
1241  return outputParticle;
1242  }
1243 
1244  //------------------------------------------------------------------------------------------------------------------------------------------
1245 
1248  const art::Ptr<recob::PFParticle> inputParticle)
1249  {
1250  // Navigate upward through PFO daughter/parent links - return the top-level non-neutrino PF Particle
1251  int primaryTrackID(inputParticle->Self());
1252 
1253  if (!inputParticle->IsPrimary()) {
1254  int parentTrackID(inputParticle->Parent());
1255 
1256  while (1) {
1257  PFParticleMap::const_iterator pIter1 = particleMap.find(parentTrackID);
1258  if (particleMap.end() == pIter1)
1259  throw cet::exception("LArPandora") << " PandoraCollector::GetFinalStatePFParticle --- "
1260  "Found a PFParticle without a particle ID ";
1261 
1262  const art::Ptr<recob::PFParticle> parentParticle = pIter1->second;
1263  if (LArPandoraHelper::IsNeutrino(parentParticle)) break;
1264 
1265  primaryTrackID = parentTrackID;
1266 
1267  if (parentParticle->IsPrimary()) break;
1268 
1269  parentTrackID = parentParticle->Parent();
1270  }
1271  }
1272 
1273  PFParticleMap::const_iterator pIter2 = particleMap.find(primaryTrackID);
1274  if (particleMap.end() == pIter2)
1275  throw cet::exception("LArPandora") << " PandoraCollector::GetFinalStatePFParticle --- Found "
1276  "a PFParticle without a particle ID ";
1277 
1278  const art::Ptr<recob::PFParticle> outputParticle = pIter2->second;
1279  return outputParticle;
1280  }
1281 
1282  //------------------------------------------------------------------------------------------------------------------------------------------
1283 
1286  const art::Ptr<simb::MCParticle> inputParticle)
1287  {
1288  // Navigate upward through MC daughter/parent links - return the top-level MC particle
1289  int primaryTrackID(inputParticle->TrackId());
1290  int parentTrackID(inputParticle->Mother());
1291 
1292  while (1) {
1293  MCParticleMap::const_iterator pIter1 = particleMap.find(parentTrackID);
1294  if (particleMap.end() == pIter1) break; // Can't find MC Particle for this track ID [break]
1295 
1296  const art::Ptr<simb::MCParticle> particle = pIter1->second;
1297 
1298  primaryTrackID = parentTrackID;
1299  parentTrackID = particle->Mother();
1300  }
1301 
1302  MCParticleMap::const_iterator pIter2 = particleMap.find(primaryTrackID);
1303  if (particleMap.end() == pIter2)
1304  throw cet::exception("LArPandora")
1305  << " PandoraCollector::GetParentMCParticle --- Found a track ID without a MC particle ";
1306 
1307  const art::Ptr<simb::MCParticle> outputParticle = pIter2->second;
1308  return outputParticle;
1309  }
1310 
1311  //------------------------------------------------------------------------------------------------------------------------------------------
1312 
1315  const art::Ptr<simb::MCParticle> inputParticle)
1316  {
1317  // Navigate upward through MC daughter/parent links - collect this particle and all its parents
1318  MCParticleVector mcVector;
1319 
1320  int trackID(inputParticle->TrackId());
1321 
1322  while (1) {
1323  MCParticleMap::const_iterator pIter = particleMap.find(trackID);
1324  if (particleMap.end() == pIter) break; // Can't find MC Particle for this track ID [break]
1325 
1326  const art::Ptr<simb::MCParticle> particle = pIter->second;
1327  mcVector.push_back(particle);
1328 
1329  trackID = particle->Mother();
1330  }
1331 
1332  // Navigate downward through MC parent/daughter links - return the first long-lived charged particle
1333  for (MCParticleVector::const_reverse_iterator iter = mcVector.rbegin(),
1334  iterEnd = mcVector.rend();
1335  iter != iterEnd;
1336  ++iter) {
1337  const art::Ptr<simb::MCParticle> nextParticle = *iter;
1338 
1339  if (LArPandoraHelper::IsVisible(nextParticle)) return nextParticle;
1340  }
1341 
1342  throw cet::exception("LArPandora"); // need to catch this exception
1343  }
1344 
1345  //------------------------------------------------------------------------------------------------------------------------------------------
1346 
1349  const art::Ptr<recob::PFParticle> particle)
1350  {
1351  PFParticlesToTracks::const_iterator tIter = particlesToTracks.find(particle);
1352 
1353  if (particlesToTracks.end() == tIter || tIter->second.empty())
1354  throw cet::exception("LArPandora")
1355  << " PandoraCollector::GetPrimaryTrack --- Failed to find associated track ";
1356 
1357  if (tIter->second.size() != 1)
1358  throw cet::exception("LArPandora")
1359  << " PandoraCollector::GetPrimaryTrack --- Found more than one associated track ";
1360 
1361  const art::Ptr<recob::Track> primaryTrack = *(tIter->second.begin());
1362  return primaryTrack;
1363  }
1364 
1365  //------------------------------------------------------------------------------------------------------------------------------------------
1366 
1367  int
1369  const art::Ptr<recob::PFParticle> inputParticle)
1370  {
1371  // Navigate upward through PFO daughter/parent links - return the top-level PF Particle
1372  int nGenerations(0);
1373  int primaryTrackID(inputParticle->Self());
1374 
1375  while (1) {
1376  PFParticleMap::const_iterator pIter = particleMap.find(primaryTrackID);
1377  if (particleMap.end() == pIter)
1378  throw cet::exception("LArPandora")
1379  << " PandoraCollector::GetGeneration --- Found a PFParticle without a particle ID ";
1380 
1381  ++nGenerations;
1382 
1383  const art::Ptr<recob::PFParticle> primaryParticle = pIter->second;
1384  if (primaryParticle->IsPrimary()) break;
1385 
1386  primaryTrackID = primaryParticle->Parent();
1387  }
1388 
1389  return nGenerations;
1390  }
1391 
1392  //------------------------------------------------------------------------------------------------------------------------------------------
1393 
1394  int
1396  const art::Ptr<recob::PFParticle> daughterParticle)
1397  {
1398  art::Ptr<recob::PFParticle> parentParticle =
1399  LArPandoraHelper::GetParentPFParticle(particleMap, daughterParticle);
1400 
1401  if (LArPandoraHelper::IsNeutrino(parentParticle)) return parentParticle->PdgCode();
1402 
1403  if (parentParticle->IsPrimary()) return 0;
1404 
1405  const int parentID(parentParticle->Parent());
1406 
1407  PFParticleMap::const_iterator pIter = particleMap.find(parentID);
1408  if (particleMap.end() == pIter)
1409  throw cet::exception("LArPandora")
1410  << " PandoraCollector::GetParentNeutrino --- Found a PFParticle without a particle ID ";
1411 
1412  const art::Ptr<recob::PFParticle> neutrinoParticle = pIter->second;
1413  return neutrinoParticle->PdgCode();
1414  }
1415 
1416  //------------------------------------------------------------------------------------------------------------------------------------------
1417 
1418  bool
1420  const art::Ptr<recob::PFParticle> daughterParticle)
1421  {
1422  if (LArPandoraHelper::IsNeutrino(daughterParticle)) return false;
1423 
1424  if (daughterParticle->IsPrimary()) return true;
1425 
1426  const int parentID(daughterParticle->Parent());
1427 
1428  PFParticleMap::const_iterator pIter = particleMap.find(parentID);
1429  if (particleMap.end() == pIter)
1430  throw cet::exception("LArPandora")
1431  << " PandoraCollector::IsFinalState --- Found a PFParticle without a particle ID ";
1432 
1433  const art::Ptr<recob::PFParticle> parentParticle = pIter->second;
1434 
1435  if (LArPandoraHelper::IsNeutrino(parentParticle)) return true;
1436 
1437  return false;
1438  }
1439 
1440  //------------------------------------------------------------------------------------------------------------------------------------------
1441 
1442  bool
1444  {
1445  const int pdg(particle->PdgCode());
1446 
1447  // electron, muon, tau (use Pandora PDG tables)
1448  return ((pandora::NU_E == std::abs(pdg)) || (pandora::NU_MU == std::abs(pdg)) ||
1449  (pandora::NU_TAU == std::abs(pdg)));
1450  }
1451 
1452  //------------------------------------------------------------------------------------------------------------------------------------------
1453 
1454  bool
1456  {
1457  const int pdg(particle->PdgCode());
1458 
1459  // muon, pion, proton, kaon (use Pandora PDG tables)
1460  return ((pandora::MU_MINUS == std::abs(pdg)) || (pandora::PI_PLUS == std::abs(pdg)) ||
1461  (pandora::PROTON == std::abs(pdg)) || (pandora::K_PLUS == std::abs(pdg)));
1462  }
1463 
1464  //------------------------------------------------------------------------------------------------------------------------------------------
1465 
1466  bool
1468  {
1469  const int pdg(particle->PdgCode());
1470 
1471  // electron, photon (use Pandora PDG tables)
1472  return ((pandora::E_MINUS == std::abs(pdg)) || (pandora::PHOTON == std::abs(pdg)));
1473  }
1474 
1475  //------------------------------------------------------------------------------------------------------------------------------------------
1476 
1477  bool
1479  {
1480  // Include long-lived charged particles
1481  const int pdg(particle->PdgCode());
1482 
1483  if ((pandora::E_MINUS == std::abs(pdg)) || (pandora::MU_MINUS == std::abs(pdg)) ||
1484  (pandora::PROTON == std::abs(pdg)) || (pandora::PI_PLUS == std::abs(pdg)) ||
1485  (pandora::K_PLUS == std::abs(pdg)) || (pandora::SIGMA_MINUS == std::abs(pdg)) ||
1486  (pandora::SIGMA_PLUS == std::abs(pdg)) || (pandora::HYPERON_MINUS == std::abs(pdg)) ||
1487  (pandora::PHOTON == std::abs(pdg)) || (pandora::NEUTRON == std::abs(pdg)))
1488  return true;
1489 
1490  // TODO: What about ions, neutrons, photons? (Have included neutrons and photons for now)
1491 
1492  return false;
1493  }
1494 
1495  //------------------------------------------------------------------------------------------------------------------------------------------
1496 
1498  LArPandoraHelper::GetPFParticleMetadata(const pandora::ParticleFlowObject* const pPfo)
1499  {
1500  return larpandoraobj::PFParticleMetadata(pPfo->GetPropertiesMap());
1501  }
1502 
1503  //------------------------------------------------------------------------------------------------------------------------------------------
1504  //------------------------------------------------------------------------------------------------------------------------------------------
1505 
1506  template void LArPandoraHelper::GetAssociatedHits(const art::Event&,
1507  const std::string&,
1509  HitVector&,
1510  const pandora::IntVector* const);
1511 
1512  template void LArPandoraHelper::GetAssociatedHits(const art::Event&,
1513  const std::string&,
1515  HitVector&,
1516  const pandora::IntVector* const);
1517 
1518 } // namespace lar_pandora
static void BuildPFParticleHitMaps(const PFParticleVector &particleVector, const PFParticlesToSpacePoints &particlesToSpacePoints, const SpacePointsToHits &spacePointsToHits, PFParticlesToHits &particlesToHits, HitsToPFParticles &hitsToParticles, const DaughterMode daughterMode=kUseDaughters)
Build mapping between PFParticles and Hits using PFParticle/SpacePoint/Hit maps.
static void CollectClusters(const art::Event &evt, const std::string &label, ClusterVector &clusterVector, ClustersToHits &clustersToHits)
Collect the reconstructed Clusters and associated hits from the ART event record. ...
std::map< art::Ptr< recob::Hit >, art::Ptr< simb::MCParticle > > HitsToMCParticles
int PdgCode() const
Definition: MCParticle.h:212
static void BuildPFParticleMap(const PFParticleVector &particleVector, PFParticleMap &particleMap)
Build particle maps for reconstructed particles.
std::map< int, art::Ptr< sim::SimChannel > > SimChannelMap
void cluster(In first, In last, Out result, Pred *pred)
Definition: NNClusters.h:41
static int GetParentNeutrino(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the parent neutrino PDG code (or zero for cosmics) for a given reconstructed particle...
size_t Self() const
Returns the index of this particle.
Definition: PFParticle.h:92
std::map< art::Ptr< recob::PFParticle >, ClusterVector > PFParticlesToClusters
static bool IsNeutrino(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as a neutrino.
std::string string
Definition: nybbler.cc:12
int Mother() const
Definition: MCParticle.h:213
std::map< art::Ptr< simb::MCParticle >, art::Ptr< simb::MCTruth > > MCParticlesToMCTruth
static art::Ptr< recob::PFParticle > GetFinalStatePFParticle(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the final-state parent particle by navigating up the chain of parent/daughter associations...
static void CollectWires(const art::Event &evt, const std::string &label, WireVector &wireVector)
Collect the reconstructed wires from the ART event record.
struct vector vector
static bool IsShower(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as shower-like.
int PdgCode() const
Return the type of particle as a PDG ID.
Definition: PFParticle.h:83
std::vector< art::Ptr< recob::Shower > > ShowerVector
intermediate_table::const_iterator const_iterator
int NParticles() const
Definition: MCTruth.h:75
uint8_t channel
Definition: CRTFragment.hh:201
std::vector< int > IntVector
std::vector< art::Ptr< anab::CosmicTag > > CosmicTagVector
Cluster finding and building.
float energy
energy from the particle with this trackID [MeV]
Definition: SimChannel.h:28
Particle class.
static art::Ptr< recob::PFParticle > GetParentPFParticle(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the top-level parent particle by navigating up the chain of parent/daughter associations...
std::map< art::Ptr< recob::PFParticle >, TrackVector > PFParticlesToTracks
static void CollectGeneratorMCParticles(const art::Event &evt, const std::string &label, RawMCParticleVector &particleVector)
Collect a vector of MCParticle objects from the generator in the ART event record. ATTN: This function is needed as accessing generator (opposed to Geant4) level MCParticles requires use of MCTruth block.
art framework interface to geometry description
std::map< art::Ptr< recob::PFParticle >, VertexVector > PFParticlesToVertices
std::map< int, art::Ptr< recob::PFParticle > > PFParticleMap
int TrackId() const
Definition: MCParticle.h:210
std::map< art::Ptr< recob::Track >, CosmicTagVector > TracksToCosmicTags
int TDCtick_t
Type representing a TDC tick.
Definition: RawTypes.h:25
std::map< int, art::Ptr< simb::MCParticle > > MCParticleMap
std::map< art::Ptr< recob::Hit >, TrackIDEVector > HitsToTrackIDEs
bool isValid() const noexcept
Definition: Handle.h:191
static void GetAssociatedHits(const art::Event &evt, const std::string &label, const std::vector< art::Ptr< T >> &inputVector, HitVector &associatedHits, const pandora::IntVector *const indexVector=nullptr)
Get all hits associated with input clusters.
static void CollectSpacePoints(const art::Event &evt, const std::string &label, SpacePointVector &spacePointVector, SpacePointsToHits &spacePointsToHits)
Collect the reconstructed SpacePoints and associated hits from the ART event record.
std::vector< art::Ptr< recob::Seed > > SeedVector
static void CollectSeeds(const art::Event &evt, const std::string &label, SeedVector &seedVector, PFParticlesToSeeds &particlesToSeeds)
Collect the reconstructed PFParticles and associated Seeds from the ART event record.
T abs(T value)
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
static bool IsVisible(const art::Ptr< simb::MCParticle > particle)
Determine whether a particle is visible (i.e. long-lived charged particle)
std::map< art::Ptr< recob::PFParticle >, MetadataVector > PFParticlesToMetadata
std::map< art::Ptr< recob::PFParticle >, T0Vector > PFParticlesToT0s
const double e
std::map< art::Ptr< recob::PFParticle >, ShowerVector > PFParticlesToShowers
std::vector< simb::MCParticle > RawMCParticleVector
static art::Ptr< simb::MCParticle > GetParentMCParticle(const MCParticleMap &particleMap, const art::Ptr< simb::MCParticle > daughterParticle)
Return the top-level parent particle by navigating up the chain of parent/daughter associations...
Metadata associated to PFParticles.
std::map< art::Ptr< simb::MCParticle >, HitVector > MCParticlesToHits
std::map< art::Ptr< recob::PFParticle >, SeedVector > PFParticlesToSeeds
std::vector< art::Ptr< recob::PFParticle > > PFParticleVector
size_t Parent() const
Definition: PFParticle.h:96
std::map< art::Ptr< recob::Cluster >, HitVector > ClustersToHits
std::map< art::Ptr< recob::Shower >, HitVector > ShowersToHits
std::map< art::Ptr< simb::MCTruth >, MCParticleVector > MCTruthToMCParticles
static void CollectShowers(const art::Event &evt, const std::string &label, ShowerVector &showerVector, PFParticlesToShowers &particlesToShowers)
Collect the reconstructed PFParticles and associated Showers from the ART event record.
key_type key() const noexcept
Definition: Ptr.h:216
std::vector< art::Ptr< recob::SpacePoint > > SpacePointVector
static art::Ptr< recob::Track > GetPrimaryTrack(const PFParticlesToTracks &particlesToTracks, const art::Ptr< recob::PFParticle > particle)
Return the primary track associated with a PFParticle.
std::map< art::Ptr< recob::PFParticle >, HitVector > PFParticlesToHits
std::vector< sim::TrackIDE > TrackIDEs(TDC_t startTDC, TDC_t endTDC) const
Returns energies collected for each track within a time interval.
Definition: SimChannel.cxx:246
static void CollectVertices(const art::Event &evt, const std::string &label, VertexVector &vertexVector, PFParticlesToVertices &particlesToVertices)
Collect the reconstructed PFParticles and associated Vertices from the ART event record.
float energyFrac
fraction of hit energy from the particle with this trackID
Definition: SimChannel.h:27
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > HitsToSpacePoints
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
std::vector< sim::TrackIDE > TrackIDEVector
static void CollectSimChannels(const art::Event &evt, const std::string &label, SimChannelVector &simChannelVector, bool &areSimChannelsValid)
Collect a vector of SimChannel objects from the ART event record.
bool IsPrimary() const
Returns whether the particle is the root of the flow.
Definition: PFParticle.h:86
static void SelectFinalStatePFParticles(const PFParticleVector &inputParticles, PFParticleVector &outputParticles)
Select final-state reconstructed particles from a list of all reconstructed particles.
std::vector< art::Ptr< recob::Track > > TrackVector
static void CollectPFParticleMetadata(const art::Event &evt, const std::string &label, PFParticleVector &particleVector, PFParticlesToMetadata &particlesToMetadata)
Collect the reconstructed PFParticle Metadata from the ART event record.
std::map< art::Ptr< recob::PFParticle >, SpacePointVector > PFParticlesToSpacePoints
std::map< art::Ptr< recob::Track >, HitVector > TracksToHits
static void CollectTracks(const art::Event &evt, const std::string &label, TrackVector &trackVector, PFParticlesToTracks &particlesToTracks)
Collect the reconstructed PFParticles and associated Tracks from the ART event record.
float PeakTimeMinusRMS(float sigmas=+1.) const
Definition: Hit.h:239
static void CollectPFParticles(const art::Event &evt, const std::string &label, PFParticleVector &particleVector)
Collect the reconstructed PFParticles from the ART event record.
static void CollectHits(const art::Event &evt, const std::string &label, HitVector &hitVector)
Collect the reconstructed Hits from the ART event record.
Detector simulation of raw signals on wires.
DaughterMode
DaughterMode enumeration.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
static void SelectNeutrinoPFParticles(const PFParticleVector &inputParticles, PFParticleVector &outputParticles)
Select reconstructed neutrino particles from a list of all reconstructed particles.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:329
std::vector< art::Ptr< recob::Hit > > HitVector
Declaration of signal hit object.
int trackID
Geant4 supplied trackID.
Definition: SimChannel.h:26
std::vector< art::Ptr< sim::SimChannel > > SimChannelVector
std::vector< art::Ptr< recob::Wire > > WireVector
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::PFParticle > > HitsToPFParticles
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
static bool IsTrack(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as track-like.
Provides recob::Track data product.
std::vector< art::Ptr< recob::Vertex > > VertexVector
Declaration of basic channel signal object.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
float PeakTimePlusRMS(float sigmas=+1.) const
Returns a time sigmas RMS away from the peak time.
Definition: Hit.h:236
static int GetGeneration(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the generation of this particle (first generation if primary)
TCEvent evt
Definition: DataStructs.cxx:7
std::vector< art::Ptr< anab::T0 > > T0Vector
static void CollectT0s(const art::Event &evt, const std::string &label, T0Vector &t0Vector, PFParticlesToT0s &particlesToT0s)
Collect a vector of T0s from the ART event record.
std::map< art::Ptr< recob::Seed >, art::Ptr< recob::Hit > > SeedsToHits
static void BuildMCParticleHitMaps(const art::Event &evt, const HitVector &hitVector, const SimChannelVector &simChannelVector, HitsToTrackIDEs &hitsToTrackIDEs)
Collect the links from reconstructed hits to their true energy deposits.
static void CollectMCParticles(const art::Event &evt, const std::string &label, MCParticleVector &particleVector)
Collect a vector of MCParticle objects from the ART event record.
helper function for LArPandoraInterface producer module
Ionization energy from a Geant4 track.
Definition: SimChannel.h:25
raw::ChannelID_t Channel() const
ID of the readout channel the hit was extracted from.
Definition: Hit.h:230
static art::Ptr< simb::MCParticle > GetFinalStateMCParticle(const MCParticleMap &particleMap, const art::Ptr< simb::MCParticle > daughterParticle)
Return the final-state parent particle by navigating up the chain of parent/daughter associations...
static bool IsFinalState(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Determine whether a particle has been reconstructed as a final-state particle.
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > SpacePointsToHits
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
static larpandoraobj::PFParticleMetadata GetPFParticleMetadata(const pandora::ParticleFlowObject *const pPfo)
Get metadata associated to a PFO.
QTextStream & endl(QTextStream &s)
static void CollectCosmicTags(const art::Event &evt, const std::string &label, CosmicTagVector &cosmicTagVector, TracksToCosmicTags &tracksToCosmicTags)
Collect a vector of cosmic tags from the ART event record.
static void BuildMCParticleMap(const MCParticleVector &particleVector, MCParticleMap &particleMap)
Build particle maps for true particles.
vertex reconstruction