LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
Cluster3D_module.cc
Go to the documentation of this file.
1 
46 // Framework Includes
53 #include "art_root_io/TFileService.h"
54 #include "cetlib/cpu_timer.h"
55 
56 // LArSoft includes
71 
86 
87 // ROOT includes
88 #include "TTree.h"
89 #include "TVector3.h"
90 
91 // std includes
92 #include <iostream>
93 #include <memory>
94 #include <string>
95 
96 //------------------------------------------------------------------------------------------------------------------------------------------
97 
98 namespace lar_cluster3d {
99  using Hit3DToSPPtrMap =
100  std::unordered_map<const reco::ClusterHit3D*, art::Ptr<recob::SpacePoint>>;
103 
104  //------------------------------------------------------------------------------------------------------------------------------------------
105  // Definition of the producer module here
106 
110  class Cluster3D : public art::EDProducer {
111  public:
112  explicit Cluster3D(fhicl::ParameterSet const& pset);
113 
114  private:
115  void beginJob() override;
116  void produce(art::Event& evt) override;
117 
119  public:
121  std::string& pathName,
122  std::string& vertexName,
123  std::string& extremeName)
124  : artPCAxisVector(new std::vector<recob::PCAxis>)
125  , artPFParticleVector(new std::vector<recob::PFParticle>)
126  , artClusterVector(new std::vector<recob::Cluster>)
127  , artSpacePointVector(new std::vector<recob::SpacePoint>)
128  , artPathPointVector(new std::vector<recob::SpacePoint>)
129  , artVertexPointVector(new std::vector<recob::SpacePoint>)
130  , artExtremePointVector(new std::vector<recob::SpacePoint>)
131  , artSeedVector(new std::vector<recob::Seed>)
132  , artEdgeVector(new std::vector<recob::Edge>)
135  , artClusterAssociations(new art::Assns<recob::Cluster, recob::Hit>)
136  , artPFPartAxisAssociations(new art::Assns<recob::PFParticle, recob::PCAxis>)
137  , artPFPartClusAssociations(new art::Assns<recob::PFParticle, recob::Cluster>)
138  , artPFPartSPAssociations(new art::Assns<recob::SpacePoint, recob::PFParticle>)
139  , artPFPartPPAssociations(new art::Assns<recob::SpacePoint, recob::PFParticle>)
140  , artPFPartSeedAssociations(new art::Assns<recob::PFParticle, recob::Seed>)
141  , artPFPartEdgeAssociations(new art::Assns<recob::Edge, recob::PFParticle>)
142  , artPFPartPathEdgeAssociations(new art::Assns<recob::Edge, recob::PFParticle>)
143  , artSeedHitAssociations(new art::Assns<recob::Seed, recob::Hit>)
144  , artSPHitAssociations(new art::Assns<recob::Hit, recob::SpacePoint>)
145  , artPPHitAssociations(new art::Assns<recob::Hit, recob::SpacePoint>)
146  , artEdgeSPAssociations(new art::Assns<recob::SpacePoint, recob::Edge>)
147  , artEdgePPAssociations(new art::Assns<recob::SpacePoint, recob::Edge>)
148  , fEvt(evt)
149  , fSPPtrMaker(evt)
150  , fSPPtrMakerPath(evt, pathName)
151  , fEdgePtrMaker(evt)
152  , fEdgePtrMakerPath(evt, pathName)
153  , fPathName(pathName)
154  , fVertexName(vertexName)
155  , fExtremeName(extremeName)
156  {}
157 
159  {
161  }
162 
163  void makeSpacePointHitAssns(std::vector<recob::SpacePoint>& spacePointVector,
164  RecobHitVector& recobHits,
166  const std::string& path = "")
167  {
168  for (auto& hit : recobHits)
169  util::CreateAssn(fEvt, spacePointVector, hit, spHitAssns, path);
170  }
171 
173  {
178  artPCAxisVector->size() - 2,
179  artPCAxisVector->size());
180  }
181 
182  void makePFPartSeedAssns(size_t numSeedsStart)
183  {
186  *artSeedVector,
188  numSeedsStart,
189  artSeedVector->size());
190  }
191 
192  void makePFPartClusterAssns(size_t clusterStart)
193  {
198  clusterStart,
199  artClusterVector->size());
200  }
201 
203  std::vector<recob::SpacePoint>& spacePointVector,
205  size_t spacePointStart,
206  const std::string& instance = "")
207  {
208  for (size_t idx = spacePointStart; idx < spacePointVector.size(); idx++) {
210  util::CreateAssn(fEvt, *artPFParticleVector, spacePoint, pfPartSPAssociations);
211  }
212  }
213 
214  void makePFPartEdgeAssns(std::vector<recob::Edge>& edgeVector,
215  art::Assns<recob::Edge, recob::PFParticle>& pfPartEdgeAssociations,
216  size_t edgeStart,
217  const std::string& instance = "")
218  {
219  for (size_t idx = edgeStart; idx < edgeVector.size(); idx++) {
221 
222  util::CreateAssn(fEvt, *artPFParticleVector, edge, pfPartEdgeAssociations);
223  }
224  }
225 
226  void makeEdgeSpacePointAssns(std::vector<recob::Edge>& edgeVector,
227  RecobSpacePointVector& spacePointVector,
228  art::Assns<recob::SpacePoint, recob::Edge>& edgeSPAssociations,
229  const std::string& path = "")
230  {
231  for (auto& spacePoint : spacePointVector)
232  util::CreateAssn(fEvt, edgeVector, spacePoint, edgeSPAssociations, path);
233  }
234 
236  {
237  fEvt.put(std::move(artPCAxisVector));
238  fEvt.put(std::move(artPFParticleVector));
239  fEvt.put(std::move(artClusterVector));
240  fEvt.put(std::move(artSpacePointVector));
241  fEvt.put(std::move(artSeedVector));
242  fEvt.put(std::move(artEdgeVector));
243  fEvt.put(std::move(artPFPartAxisAssociations));
244  fEvt.put(std::move(artPFPartClusAssociations));
245  fEvt.put(std::move(artClusterAssociations));
246  fEvt.put(std::move(artPFPartSPAssociations));
247  fEvt.put(std::move(artPFPartSeedAssociations));
248  fEvt.put(std::move(artPFPartEdgeAssociations));
249  fEvt.put(std::move(artSeedHitAssociations));
250  fEvt.put(std::move(artSPHitAssociations));
251  fEvt.put(std::move(artEdgeSPAssociations));
252  fEvt.put(std::move(artPathEdgeVector), fPathName);
253  fEvt.put(std::move(artPathPointVector), fPathName);
257  fEvt.put(std::move(artPPHitAssociations), fPathName);
261  }
262 
263  art::Ptr<recob::SpacePoint> makeSpacePointPtr(size_t index, const std::string& instance = "")
264  {
265  if (empty(instance)) { return fSPPtrMaker(index); }
266  return fSPPtrMakerPath(index);
267  };
268 
269  art::Ptr<recob::Edge> makeEdgePtr(size_t index, const std::string& instance = "")
270  {
271  if (empty(instance)) { return fEdgePtrMaker(index); }
272  return fEdgePtrMakerPath(index);
273  };
274 
275  std::unique_ptr<std::vector<recob::PCAxis>> artPCAxisVector;
276  std::unique_ptr<std::vector<recob::PFParticle>> artPFParticleVector;
277  std::unique_ptr<std::vector<recob::Cluster>> artClusterVector;
278  std::unique_ptr<std::vector<recob::SpacePoint>> artSpacePointVector;
279  std::unique_ptr<std::vector<recob::SpacePoint>> artPathPointVector;
280  std::unique_ptr<std::vector<recob::SpacePoint>> artVertexPointVector;
281  std::unique_ptr<std::vector<recob::SpacePoint>> artExtremePointVector;
282  std::unique_ptr<std::vector<recob::Seed>> artSeedVector;
283  std::unique_ptr<std::vector<recob::Edge>> artEdgeVector;
284  std::unique_ptr<std::vector<recob::Edge>> artPathEdgeVector;
285  std::unique_ptr<std::vector<recob::Edge>> artVertexEdgeVector;
286 
287  std::unique_ptr<art::Assns<recob::Cluster, recob::Hit>> artClusterAssociations;
288  std::unique_ptr<art::Assns<recob::PFParticle, recob::PCAxis>> artPFPartAxisAssociations;
289  std::unique_ptr<art::Assns<recob::PFParticle, recob::Cluster>> artPFPartClusAssociations;
290  std::unique_ptr<art::Assns<recob::SpacePoint, recob::PFParticle>> artPFPartSPAssociations;
291  std::unique_ptr<art::Assns<recob::SpacePoint, recob::PFParticle>> artPFPartPPAssociations;
292  std::unique_ptr<art::Assns<recob::PFParticle, recob::Seed>> artPFPartSeedAssociations;
293  std::unique_ptr<art::Assns<recob::Edge, recob::PFParticle>> artPFPartEdgeAssociations;
294  std::unique_ptr<art::Assns<recob::Edge, recob::PFParticle>> artPFPartPathEdgeAssociations;
295  std::unique_ptr<art::Assns<recob::Seed, recob::Hit>> artSeedHitAssociations;
296  std::unique_ptr<art::Assns<recob::Hit, recob::SpacePoint>> artSPHitAssociations;
297  std::unique_ptr<art::Assns<recob::Hit, recob::SpacePoint>> artPPHitAssociations;
298  std::unique_ptr<art::Assns<recob::SpacePoint, recob::Edge>> artEdgeSPAssociations;
299  std::unique_ptr<art::Assns<recob::SpacePoint, recob::Edge>> artEdgePPAssociations;
300 
301  private:
307  std::string& fPathName;
308  std::string& fVertexName;
309  std::string& fExtremeName;
310  };
311 
317  void PrepareEvent(const art::Event& evt);
318 
322  void InitializeMonitoring();
323 
333  void findTrackSeeds(art::Event& evt,
335  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
336  std::vector<recob::Seed>& seedVec,
337  art::Assns<recob::Seed, recob::Hit>& seedHitAssns) const;
338 
345  void splitClustersWithHough(reco::ClusterParameters& clusterParameters,
346  reco::ClusterParametersList& clusterParametersList) const;
347 
357  std::vector<recob::SpacePoint>& spacePointVec,
359  reco::HitPairListPtr& clusHitPairVector,
360  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
361  Hit3DToSPPtrMap& hit3DToSPPtrMap,
362  const std::string& path = "") const;
363 
371  reco::ConvexHullKinkTupleList& clusHitPairVector) const;
372 
382  dcel2d::HalfEdgeList&) const;
383 
390  using IdxToPCAMap = std::map<size_t, const reco::PrincipalComponents*>;
391 
394  IdxToPCAMap&) const;
395 
406  ArtOutputHandler& output,
407  reco::ClusterParameters& clusterParameters,
408  size_t pfParticleParent,
409  IdxToPCAMap& idxToPCAMap,
410  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
411  Hit3DToSPPtrMap& hit3DToSPPtrMap,
412  Hit3DToSPPtrMap& best3DToSPPtrMap) const;
413 
422  ArtOutputHandler& output,
423  reco::HitPairList& hitPairList,
424  reco::ClusterParametersList& clusterParametersList,
425  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap) const;
426 
435  size_t ConvertToArtOutput(util::GeometryUtilities const& gser,
436  ArtOutputHandler& output,
437  reco::ClusterParameters& clusterParameters,
438  size_t pfParticleParent,
439  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
440  Hit3DToSPPtrMap& hit3DToSPPtrMap,
441  Hit3DToSPPtrMap& best3DToSPPtrMap) const;
442 
450  {
451  return fabs(pca.getEigenVectors().row(2)(0)) > m_parallelHitsCosAng &&
452  3. * sqrt(pca.getEigenValues()(1)) > m_parallelHitsTransWid;
453  }
454 
460  size_t countUltimateDaughters(reco::ClusterParameters& clusterParameters) const;
461 
469 
473  TTree* m_pRecoTree;
474  int m_run;
475  int m_event;
476  int m_hits;
477  int m_hits3D;
478  float m_totalTime;
482  float m_dbscanTime;
485  float m_finishTime;
486  std::string m_pathInstance;
487  std::string m_vertexInstance;
488  std::string m_extremeInstance;
489 
490  // Algorithms
491  std::unique_ptr<lar_cluster3d::IHit3DBuilder>
493  std::unique_ptr<lar_cluster3d::IClusterAlg>
495  std::unique_ptr<lar_cluster3d::IClusterModAlg>
497  std::unique_ptr<lar_cluster3d::IClusterModAlg>
499  std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
506  };
507 
509 
510 } // namespace lar_cluster3d
511 
512 //------------------------------------------------------------------------------------------------------------------------------------------
513 // implementation follows
514 
515 namespace lar_cluster3d {
516 
518  : EDProducer{pset}
519  , m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
520  , m_skeletonAlg(pset.get<fhicl::ParameterSet>("SkeletonAlg"))
521  , m_seedFinderAlg(pset.get<fhicl::ParameterSet>("SeedFinderAlg"))
522  , m_pcaSeedFinderAlg(pset.get<fhicl::ParameterSet>("PCASeedFinderAlg"))
523  , m_parallelHitsAlg(pset.get<fhicl::ParameterSet>("ParallelHitsAlg"))
524  {
525  m_onlyMakSpacePoints = pset.get<bool>("MakeSpacePointsOnly", false);
526  m_enableMonitoring = pset.get<bool>("EnableMonitoring", false);
527  m_parallelHitsCosAng = pset.get<float>("ParallelHitsCosAng", 0.999);
528  m_parallelHitsTransWid = pset.get<float>("ParallelHitsTransWid", 25.0);
529  m_pathInstance = pset.get<std::string>("PathPointsName", "Path");
530  m_vertexInstance = pset.get<std::string>("VertexPointsName", "Vertex");
531  m_extremeInstance = pset.get<std::string>("ExtremePointsName", "Extreme");
532 
533  m_hit3DBuilderAlg = art::make_tool<lar_cluster3d::IHit3DBuilder>(
534  pset.get<fhicl::ParameterSet>("Hit3DBuilderAlg"));
535  m_clusterAlg =
536  art::make_tool<lar_cluster3d::IClusterAlg>(pset.get<fhicl::ParameterSet>("ClusterAlg"));
537  m_clusterMergeAlg = art::make_tool<lar_cluster3d::IClusterModAlg>(
538  pset.get<fhicl::ParameterSet>("ClusterMergeAlg"));
539  m_clusterPathAlg = art::make_tool<lar_cluster3d::IClusterModAlg>(
540  pset.get<fhicl::ParameterSet>("ClusterPathAlg"));
541  m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
542  pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
543 
544  // Handle special case of Space Point building outputting a new hit collection
545  m_hit3DBuilderAlg->produces(producesCollector());
546 
547  produces<std::vector<recob::PCAxis>>();
548  produces<std::vector<recob::PFParticle>>();
549  produces<std::vector<recob::Cluster>>();
550  produces<std::vector<recob::Seed>>();
551  produces<std::vector<recob::Edge>>();
552 
553  produces<art::Assns<recob::PFParticle, recob::PCAxis>>();
554  produces<art::Assns<recob::PFParticle, recob::Cluster>>();
555  produces<art::Assns<recob::PFParticle, recob::Seed>>();
556  produces<art::Assns<recob::Edge, recob::PFParticle>>();
557  produces<art::Assns<recob::Seed, recob::Hit>>();
558  produces<art::Assns<recob::Cluster, recob::Hit>>();
559 
560  produces<std::vector<recob::SpacePoint>>();
561  produces<art::Assns<recob::SpacePoint, recob::PFParticle>>();
562  produces<art::Assns<recob::Hit, recob::SpacePoint>>();
563  produces<art::Assns<recob::SpacePoint, recob::Edge>>();
564 
565  produces<std::vector<recob::SpacePoint>>(m_pathInstance);
566  produces<std::vector<recob::Edge>>(m_pathInstance);
567  produces<art::Assns<recob::SpacePoint, recob::PFParticle>>(m_pathInstance);
568  produces<art::Assns<recob::Edge, recob::PFParticle>>(m_pathInstance);
569  produces<art::Assns<recob::Hit, recob::SpacePoint>>(m_pathInstance);
570  produces<art::Assns<recob::SpacePoint, recob::Edge>>(m_pathInstance);
571 
572  produces<std::vector<recob::Edge>>(m_vertexInstance);
573  produces<std::vector<recob::SpacePoint>>(m_vertexInstance);
574 
575  produces<std::vector<recob::SpacePoint>>(m_extremeInstance);
576  }
577 
578  //------------------------------------------------------------------------------------------------------------------------------------------
579 
581  {
587  }
588 
589  //------------------------------------------------------------------------------------------------------------------------------------------
590 
592  {
593  mf::LogInfo("Cluster3D") << " *** Cluster3D::produce(...) [Run=" << evt.run()
594  << ", Event=" << evt.id().event() << "] Starting Now! *** "
595  << std::endl;
596 
597  // Set up for monitoring the timing... at some point this should be removed in favor of
598  // external profilers
599  cet::cpu_timer theClockTotal;
600  cet::cpu_timer theClockFinish;
601 
602  if (m_enableMonitoring) theClockTotal.start();
603 
604  // This really only does anything if we are monitoring since it clears our tree variables
605  this->PrepareEvent(evt);
606 
607  // Get instances of the primary data structures needed
608  reco::ClusterParametersList clusterParametersList;
609  IHit3DBuilder::RecobHitToPtrMap clusterHitToArtPtrMap;
610  std::unique_ptr<reco::HitPairList> hitPairList(
611  new reco::HitPairList); // Potentially lots of hits, use heap instead of stack
612 
613  // Call the algorithm that builds 3D hits and stores the hit collection
614  m_hit3DBuilderAlg->Hit3DBuilder(evt, *hitPairList, clusterHitToArtPtrMap);
615 
616  // Only do the rest if we are not in the mode of only building space points (requested by ML folks)
617  if (!m_onlyMakSpacePoints) {
618  // Call the main workhorse algorithm for building the local version of candidate 3D clusters
619  m_clusterAlg->Cluster3DHits(*hitPairList, clusterParametersList);
620 
621  // Try merging clusters
622  m_clusterMergeAlg->ModifyClusters(clusterParametersList);
623 
624  // Run the path finding
625  m_clusterPathAlg->ModifyClusters(clusterParametersList);
626  }
627 
628  if (m_enableMonitoring) theClockFinish.start();
629 
630  // Get the art ouput object
632 
633  // Call the module that does the end processing (of which there is quite a bit of work!)
634  // This goes here to insure that something is always written to the data store
635  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
636  auto const detProp =
638  util::GeometryUtilities const gser{*lar::providerFrom<geo::Geometry>(),
640  clockData,
641  detProp};
642  ProduceArtClusters(gser, output, *hitPairList, clusterParametersList, clusterHitToArtPtrMap);
643 
644  // Output to art
645  output.outputObjects();
646 
647  // If monitoring then deal with the fallout
648  if (m_enableMonitoring) {
649  theClockFinish.stop();
650  theClockTotal.stop();
651 
652  m_run = evt.run();
653  m_event = evt.id().event();
654  m_totalTime = theClockTotal.accumulated_real_time();
658  m_dbscanTime = m_clusterAlg->getTimeToExecute(IClusterAlg::RUNDBSCAN) +
659  m_clusterAlg->getTimeToExecute(IClusterAlg::BUILDCLUSTERINFO);
660  m_clusterMergeTime = m_clusterMergeAlg->getTimeToExecute();
661  m_pathFindingTime = m_clusterPathAlg->getTimeToExecute();
662  m_finishTime = theClockFinish.accumulated_real_time();
663  m_hits = static_cast<int>(clusterHitToArtPtrMap.size());
664  m_hits3D = static_cast<int>(hitPairList->size());
665  m_pRecoTree->Fill();
666 
667  mf::LogDebug("Cluster3D") << "*** Cluster3D total time: " << m_totalTime
668  << ", art: " << m_artHitsTime << ", make: " << m_makeHitsTime
669  << ", build: " << m_buildNeighborhoodTime
670  << ", clustering: " << m_dbscanTime
671  << ", merge: " << m_clusterMergeTime
672  << ", path: " << m_pathFindingTime << ", finish: " << m_finishTime
673  << std::endl;
674  }
675 
676  // Will we ever get here? ;-)
677  return;
678  }
679 
680  //------------------------------------------------------------------------------------------------------------------------------------------
681 
683  {
685  m_pRecoTree = tfs->make<TTree>("monitoring", "LAr Reco");
686  m_pRecoTree->Branch("run", &m_run, "run/I");
687  m_pRecoTree->Branch("event", &m_event, "event/I");
688  m_pRecoTree->Branch("hits", &m_hits, "hits/I");
689  m_pRecoTree->Branch("hits3D", &m_hits3D, "hits3D/I");
690  m_pRecoTree->Branch("totalTime", &m_totalTime, "time/F");
691  m_pRecoTree->Branch("artHitsTime", &m_artHitsTime, "time/F");
692  m_pRecoTree->Branch("makeHitsTime", &m_makeHitsTime, "time/F");
693  m_pRecoTree->Branch("buildneigborhoodTime", &m_buildNeighborhoodTime, "time/F");
694  m_pRecoTree->Branch("dbscanTime", &m_dbscanTime, "time/F");
695  m_pRecoTree->Branch("clusterMergeTime", &m_clusterMergeTime, "time/F");
696  m_pRecoTree->Branch("pathfindingtime", &m_pathFindingTime, "time/F");
697  m_pRecoTree->Branch("finishTime", &m_finishTime, "time/F");
698 
699  m_clusterPathAlg->initializeHistograms(*tfs.get());
700  }
701 
702  //------------------------------------------------------------------------------------------------------------------------------------------
703 
705  {
706  m_run = evt.run();
707  m_event = evt.id().event();
708  m_hits = 0;
709  m_hits3D = 0;
710  m_totalTime = 0.f;
711  m_artHitsTime = 0.f;
712  m_makeHitsTime = 0.f;
714  m_dbscanTime = 0.f;
715  m_pathFindingTime = 0.f;
716  m_finishTime = 0.f;
717  }
718 
719  //------------------------------------------------------------------------------------------------------------------------------------------
720 
723  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
724  std::vector<recob::Seed>& seedVec,
725  art::Assns<recob::Seed, recob::Hit>& seedHitAssns) const
726  {
732  // Make sure we are using the right pca
733  reco::PrincipalComponents& fullPCA = cluster.getFullPCA();
734  reco::PrincipalComponents& skeletonPCA = cluster.getSkeletonPCA();
735  reco::HitPairListPtr& hitPairListPtr = cluster.getHitPairListPtr();
736  reco::HitPairListPtr skeletonListPtr;
737 
738  // We want to work with the "skeleton" hits so first step is to call the algorithm to
739  // recover only these hits from the entire input collection
740  m_skeletonAlg.GetSkeletonHits(hitPairListPtr, skeletonListPtr);
741 
742  // Skeleton hits are nice but we can do better if we then make a pass through to "average"
743  // the skeleton hits position in the Y-Z plane
744  m_skeletonAlg.AverageSkeletonPositions(skeletonListPtr);
745 
746  SeedHitPairListPairVec seedHitPairVec;
747 
748  // Some combination of the elements below will be used to determine which seed finding algorithm
749  // to pursue below
750  float eigenVal0 = 3. * sqrt(skeletonPCA.getEigenValues()[0]);
751  float eigenVal1 = 3. * sqrt(skeletonPCA.getEigenValues()[1]);
752  float eigenVal2 = 3. * sqrt(skeletonPCA.getEigenValues()[2]);
753  float transRMS = std::hypot(eigenVal0, eigenVal1);
754 
755  bool foundGoodSeed(false);
756 
757  // Choose a method for finding the seeds based on the PCA that was run...
758  // Currently we have an ad hoc if-else block which I hope will be improved soon!
759  if (aParallelHitsCluster(fullPCA)) {
760  // In this case we have a track moving relatively parallel to the wire plane with lots of
761  // ambiguous 3D hits. Your best bet here is to use the "parallel hits" algorithm to get the
762  // best axis and seeds
763  // This algorithm does not fail (foundGoodSeed will always return true)
764  foundGoodSeed = m_parallelHitsAlg.findTrackSeeds(hitPairListPtr, skeletonPCA, seedHitPairVec);
765  }
766  else if (eigenVal2 > 40. && transRMS < 5.) {
767  // If the input cluster is relatively "straight" then chances are it is a single straight track,
768  // probably a CR muon, and we can simply use the PCA to determine the seed
769  // This algorithm will check both "ends" of the input hits and if the angles become inconsistent
770  // then it will "fail"
771  foundGoodSeed =
772  m_pcaSeedFinderAlg.findTrackSeeds(skeletonListPtr, skeletonPCA, seedHitPairVec);
773  }
774 
775  // In the event the above two methods failed then we hit it with the real seed finder
776  if (!foundGoodSeed) {
777  // If here then we have a complicated 3D cluster and we'll use the hough transform algorithm to
778  // return a list of candidate seeds and seed hits
779  m_seedFinderAlg.findTrackSeeds(skeletonListPtr, skeletonPCA, seedHitPairVec);
780  }
781 
782  // Go through the returned lists and build out the art friendly seeds and hits
783  for (const auto& seedHitPair : seedHitPairVec) {
784  seedVec.push_back(seedHitPair.first);
785 
786  // We use a set here because our 3D hits can share 2D hits
787  // The set will make sure we get unique combinations of 2D hits
788  std::set<art::Ptr<recob::Hit>> seedHitSet;
789  for (const auto& hit3D : seedHitPair.second) {
790  for (const auto& hit2D : hit3D->getHits()) {
791  if (!hit2D) continue;
792 
793  const recob::Hit* recobHit = hit2D->getHit();
794  seedHitSet.insert(hitToPtrMap[recobHit]);
795  }
796  }
797 
798  RecobHitVector seedHitVec;
799  for (const auto& hit2D : seedHitSet)
800  seedHitVec.push_back(hit2D);
801 
802  util::CreateAssn(evt, seedVec, seedHitVec, seedHitAssns);
803  }
804  }
805 
807  bool operator()(const std::pair<float, const reco::ClusterHit3D*>& left,
808  const std::pair<float, const reco::ClusterHit3D*>& right)
809  {
810  return left.first < right.first;
811  }
812  };
813 
815  public:
816  CopyIfInRange(float maxRange) : m_maxRange(maxRange) {}
817 
818  bool operator()(const reco::ClusterHit3D* hit3D) const
819  {
820  return hit3D->getDocaToAxis() < m_maxRange;
821  }
822 
823  private:
824  float m_maxRange;
825  };
826 
828  reco::ClusterParametersList& clusterParametersList) const
829  {
830  // @brief A method for splitted "crossed tracks" clusters into separate clusters
831  //
832  // If this routine is called then we believe we have a cluster which needs splitting.
833  // The specific topology we are looking for is two long straight tracks which cross at some
834  // point in close proximity so their hits were joined into a single 3D cluster. The method
835  // to split this topology is to let the hough transform algorithm find the two leading candidates
836  // and then to see if we use those to build two clusters instead of one.
837 
838  // Recover the hits we'll work on.
839  // Note that we use on the skeleton hits so will need to recover them
840  reco::HitPairListPtr& hitPairListPtr = clusterParameters.getHitPairListPtr();
841  reco::HitPairListPtr skeletonListPtr;
842 
843  // We want to work with the "skeleton" hits so first step is to call the algorithm to
844  // recover only these hits from the entire input collection
845  m_skeletonAlg.GetSkeletonHits(hitPairListPtr, skeletonListPtr);
846 
847  // Skeleton hits are nice but we can do better if we then make a pass through to "average"
848  // the skeleton hits position in the Y-Z plane
849  m_skeletonAlg.AverageSkeletonPositions(skeletonListPtr);
850 
851  // Define the container for our lists of hits
852  reco::HitPairListPtrList hitPairListPtrList;
853 
854  // Now feed this to the Hough Transform to find candidate straight lines
856  skeletonListPtr, clusterParameters.getSkeletonPCA(), hitPairListPtrList);
857 
858  // We need at least two lists or else there is nothing to do
859  if (hitPairListPtrList.size() < 2) return;
860 
861  // The game plan will be the following:
862  // 1) Take the first list of hits and run the PCA on this to get an axis
863  // - Then calculate the 3d doca for ALL hits in the cluster to this axis
864  // - Move all hits within "3 sigam" of the axis to a new list
865  // 2) run the PCA on the second list of hits to get that axis
866  // - Then calculate the 3d doca for all hits in our first list
867  // - Copy hits in the first list which are within 3 sigma of the new axis
868  // back into our original cluster - these are shared hits
869  reco::HitPairListPtrList::iterator hitPairListIter = hitPairListPtrList.begin();
870  reco::HitPairListPtr& firstHitList = *hitPairListIter++;
871  reco::PrincipalComponents firstHitListPCA;
872 
873  m_pcaAlg.PCAAnalysis_3D(firstHitList, firstHitListPCA);
874 
875  // Make sure we have a successful calculation.
876  if (firstHitListPCA.getSvdOK()) {
877  // The fill routines below will expect to see unused 2D hits so we need to clear the
878  // status bits... and I am not sure of a better way...
879  for (const auto& hit3D : hitPairListPtr) {
880  for (const auto& hit2D : hit3D->getHits())
881  if (hit2D) hit2D->clearStatusBits(0x1);
882  }
883 
884  // Calculate the 3D doca's for the hits which were used to make this PCA
885  m_pcaAlg.PCAAnalysis_calc3DDocas(firstHitList, firstHitListPCA);
886 
887  // Divine from the ether some maximum allowed range for transfering hits
888  float allowedHitRange = 6. * firstHitListPCA.getAveHitDoca();
889 
890  // Now go through and calculate the 3D doca's for ALL the hits in the original cluster
891  m_pcaAlg.PCAAnalysis_calc3DDocas(hitPairListPtr, firstHitListPCA);
892 
893  // Let's make a new cluster to hold the hits
894  clusterParametersList.push_back(reco::ClusterParameters());
895 
896  // Can we get a reference to what we just created?
897  reco::ClusterParameters& newClusterParams = clusterParametersList.back();
898 
899  reco::HitPairListPtr& newClusterHitList = newClusterParams.getHitPairListPtr();
900 
901  newClusterHitList.resize(hitPairListPtr.size());
902 
903  // Do the actual copy of the hits we want
904  reco::HitPairListPtr::iterator newListEnd = std::copy_if(hitPairListPtr.begin(),
905  hitPairListPtr.end(),
906  newClusterHitList.begin(),
907  CopyIfInRange(allowedHitRange));
908 
909  // Shrink to fit
910  newClusterHitList.resize(std::distance(newClusterHitList.begin(), newListEnd));
911 
912  // And now remove these hits from the original cluster
913  hitPairListPtr.remove_if(CopyIfInRange(allowedHitRange));
914 
915  // Get an empty hit to cluster map...
916  reco::Hit2DToClusterMap hit2DToClusterMap;
917 
918  // Now "fill" the cluster parameters but turn off the hit rejection
919  m_clusterBuilder->FillClusterParams(newClusterParams, hit2DToClusterMap, 0., 1.);
920 
921  // Set the skeleton pca to the value calculated above on input
922  clusterParameters.getSkeletonPCA() = firstHitListPCA;
923 
924  // We are done with splitting out one track. Because the two tracks cross in
925  // close proximity, this is the one case where we might consider sharing 3D hits
926  // So let's make a little detour here to try to copy some of those hits back into
927  // the main hit list
928  reco::HitPairListPtr& secondHitList = *hitPairListIter;
929  reco::PrincipalComponents secondHitListPCA;
930 
931  m_pcaAlg.PCAAnalysis_3D(secondHitList, secondHitListPCA);
932 
933  // Make sure we have a successful calculation.
934  if (secondHitListPCA.getSvdOK()) {
935  // Calculate the 3D doca's for the hits which were used to make this PCA
936  m_pcaAlg.PCAAnalysis_calc3DDocas(secondHitList, secondHitListPCA);
937 
938  // Since this is the "other" cluster, we'll be a bit more generous in adding back hits
939  float newAllowedHitRange = 6. * secondHitListPCA.getAveHitDoca();
940 
941  // Go through and calculate the 3D doca's for the hits in our new candidate cluster
942  m_pcaAlg.PCAAnalysis_calc3DDocas(newClusterHitList, secondHitListPCA);
943 
944  // Create a temporary list to fill with the hits we might want to save
945  reco::HitPairListPtr tempHitList(newClusterHitList.size());
946 
947  // Do the actual copy of the hits we want...
948  reco::HitPairListPtr::iterator tempListEnd =
949  std::copy_if(newClusterHitList.begin(),
950  newClusterHitList.end(),
951  tempHitList.begin(),
952  CopyIfInRange(newAllowedHitRange));
953 
954  hitPairListPtr.insert(hitPairListPtr.end(), tempHitList.begin(), tempListEnd);
955  }
956 
957  // Of course, now we need to modify the original cluster parameters
958  reco::ClusterParameters originalParams(hitPairListPtr);
959 
960  // Now "fill" the cluster parameters but turn off the hit rejection
961  m_clusterBuilder->FillClusterParams(originalParams, hit2DToClusterMap, 0., 1.);
962 
963  // Overwrite original cluster parameters with our new values
964  clusterParameters.getClusterParams() = originalParams.getClusterParams();
965  clusterParameters.getFullPCA() = originalParams.getFullPCA();
966  clusterParameters.getSkeletonPCA() = secondHitListPCA;
967  }
968  }
969 
971  ArtOutputHandler& output,
972  reco::HitPairList& hitPairVector,
973  reco::ClusterParametersList& clusterParametersList,
974  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap) const
975  {
980  mf::LogDebug("Cluster3D") << " *** Cluster3D::ProduceArtClusters() *** " << std::endl;
981 
982  // Make sure there is something to do here!
983  if (!clusterParametersList.empty()) {
984  // This is the loop over candidate 3D clusters
985  // Note that it might be that the list of candidate clusters is modified by splitting
986  // So we use the following construct to make sure we get all of them
987  for (auto& clusterParameters : clusterParametersList) {
988  // It should be straightforward at this point to transfer information from our vector of clusters
989  // to the larsoft objects... of course we still have some work to do first, in particular to
990  // find the candidate seeds and their seed hits
991 
992  // The chances of getting here and this condition not being true are probably zero... but check anyway
993  if (!clusterParameters.getFullPCA().getSvdOK()) {
994  mf::LogDebug("Cluster3D")
995  << "--> no feature extraction done on this cluster!!" << std::endl;
996  continue;
997  }
998 
999  // Keep track of hit 3D to SP for when we do edges
1000  Hit3DToSPPtrMap hit3DToSPPtrMap;
1001  Hit3DToSPPtrMap best3DToSPPtrMap;
1002 
1003  // Do a special output of voronoi vertices here...
1004  dcel2d::VertexList& vertexList = clusterParameters.getVertexList();
1005  dcel2d::HalfEdgeList& halfEdgeList = clusterParameters.getHalfEdgeList();
1006 
1007  MakeAndSaveVertexPoints(output, vertexList, halfEdgeList);
1008 
1009  // Special case handling... if no daughters then call standard conversion routine to make sure space points
1010  // created, etc.
1011  if (clusterParameters.daughterList().empty()) {
1012  ConvertToArtOutput(gser,
1013  output,
1014  clusterParameters,
1016  hitToPtrMap,
1017  hit3DToSPPtrMap,
1018  best3DToSPPtrMap);
1019 
1020  // Get the extreme points
1021  MakeAndSaveKinkPoints(output,
1022  clusterParameters.getConvexHull().getConvexHullKinkPoints());
1023  }
1024  // Otherwise, the cluster has daughters so we handle specially
1025  else {
1026  // Set up to keep track of parent/daughters
1027  IdxToPCAMap idxToPCAMap;
1028  size_t numTotalDaughters = countUltimateDaughters(clusterParameters);
1029  size_t pfParticleIdx(output.artPFParticleVector->size() + numTotalDaughters);
1030 
1031  FindAndStoreDaughters(gser,
1032  output,
1033  clusterParameters,
1034  pfParticleIdx,
1035  idxToPCAMap,
1036  hitToPtrMap,
1037  hit3DToSPPtrMap,
1038  best3DToSPPtrMap);
1039 
1040  // Need to make a daughter vec from our map
1041  std::vector<size_t> daughterVec;
1042 
1043  for (auto& idxToPCA : idxToPCAMap)
1044  daughterVec.emplace_back(idxToPCA.first);
1045 
1046  // Now create/handle the parent PFParticle
1047  recob::PFParticle pfParticle(
1048  13, pfParticleIdx, recob::PFParticle::kPFParticlePrimary, daughterVec);
1049  output.artPFParticleVector->push_back(pfParticle);
1050 
1051  recob::PCAxis::EigenVectors eigenVecs;
1052  double eigenVals[] = {0., 0., 0.};
1053  double avePosition[] = {0., 0., 0.};
1054 
1055  eigenVecs.resize(3);
1056 
1057  reco::PrincipalComponents& skeletonPCA = clusterParameters.getSkeletonPCA();
1058 
1059  for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
1060  avePosition[outerIdx] = skeletonPCA.getAvePosition()[outerIdx];
1061  eigenVals[outerIdx] = skeletonPCA.getEigenValues()[outerIdx];
1062 
1063  eigenVecs[outerIdx].resize(3);
1064 
1065  // Be careful here... eigen stores in column major order buy default
1066  for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
1067  eigenVecs[outerIdx][innerIdx] = skeletonPCA.getEigenVectors().row(outerIdx)[innerIdx];
1068  }
1069 
1070  recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
1071  skeletonPCA.getNumHitsUsed(),
1072  eigenVals,
1073  eigenVecs,
1074  avePosition,
1075  skeletonPCA.getAveHitDoca(),
1076  output.artPCAxisVector->size());
1077 
1078  output.artPCAxisVector->push_back(skelPcAxis);
1079 
1080  reco::PrincipalComponents& fullPCA = clusterParameters.getFullPCA();
1081 
1082  for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
1083  avePosition[outerIdx] = fullPCA.getAvePosition()[outerIdx];
1084  eigenVals[outerIdx] = fullPCA.getEigenValues()[outerIdx];
1085 
1086  for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
1087  eigenVecs[outerIdx][innerIdx] = fullPCA.getEigenVectors().row(outerIdx)(innerIdx);
1088  }
1089 
1090  recob::PCAxis fullPcAxis(fullPCA.getSvdOK(),
1091  fullPCA.getNumHitsUsed(),
1092  eigenVals,
1093  eigenVecs,
1094  avePosition,
1095  fullPCA.getAveHitDoca(),
1096  output.artPCAxisVector->size());
1097 
1098  output.artPCAxisVector->push_back(fullPcAxis);
1099 
1100  // Create associations to the PFParticle
1101  output.makePFPartPCAAssns();
1102 
1103  // Make associations to all space points for this cluster
1104  MakeAndSaveSpacePoints(output,
1105  *output.artSpacePointVector.get(),
1106  *output.artSPHitAssociations.get(),
1107  clusterParameters.getHitPairListPtr(),
1108  hitToPtrMap,
1109  hit3DToSPPtrMap);
1110 
1111  // Get the extreme points
1112  MakeAndSaveKinkPoints(output,
1113  clusterParameters.getConvexHull().getConvexHullKinkPoints());
1114 
1115  // Build the edges now
1116  size_t edgeStart(output.artEdgeVector->size());
1117 
1118  for (const auto& edge : clusterParameters.getConvexHull().getConvexHullEdgeList()) {
1119  RecobSpacePointVector spacePointVec;
1120 
1121  try {
1122  spacePointVec.push_back(hit3DToSPPtrMap.at(std::get<0>(edge)));
1123  spacePointVec.push_back(hit3DToSPPtrMap.at(std::get<1>(edge)));
1124  }
1125  catch (...) {
1126  std::cout << "Caught exception in looking up space point ptr... " << std::get<0>(edge)
1127  << ", " << std::get<1>(edge) << std::endl;
1128  continue;
1129  }
1130 
1131  output.artEdgeVector->emplace_back(std::get<2>(edge),
1132  spacePointVec[0].key(),
1133  spacePointVec[1].key(),
1134  output.artEdgeVector->size());
1135 
1136  output.makeEdgeSpacePointAssns(
1137  *output.artEdgeVector.get(), spacePointVec, *output.artEdgeSPAssociations.get());
1138  }
1139 
1140  output.makePFPartEdgeAssns(
1141  *output.artEdgeVector.get(), *output.artPFPartEdgeAssociations.get(), edgeStart);
1142  }
1143  }
1144  }
1145 
1146  // Right now error matrix is uniform...
1147  int nFreePoints(0);
1148 
1149  // Run through the HitPairVector and add any unused hit pairs to the list
1150  for (auto& hitPair : hitPairVector) {
1151  if (hitPair.bitsAreSet(reco::ClusterHit3D::MADESPACEPOINT)) continue;
1152 
1153  double spacePointPos[] = {
1154  hitPair.getPosition()[0], hitPair.getPosition()[1], hitPair.getPosition()[2]};
1155  double spacePointErr[] = {1., 0., 0., 1., 0., 1.};
1156  double chisq = hitPair.getHitChiSquare();
1157 
1158  RecobHitVector recobHits;
1159 
1160  for (const auto hit : hitPair.getHits()) {
1161  if (!hit) {
1162  chisq = -1000.;
1163  continue;
1164  }
1165 
1166  art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit->getHit()];
1167  recobHits.push_back(hitPtr);
1168  }
1169 
1170  nFreePoints++;
1171 
1172  spacePointErr[1] = hitPair.getTotalCharge();
1173  spacePointErr[3] = hitPair.getChargeAsymmetry();
1174 
1175  output.artSpacePointVector->push_back(
1176  recob::SpacePoint(spacePointPos, spacePointErr, chisq, output.artSpacePointVector->size()));
1177 
1178  if (!recobHits.empty())
1179  output.makeSpacePointHitAssns(
1180  *output.artSpacePointVector.get(), recobHits, *output.artSPHitAssociations.get());
1181  }
1182 
1183  std::cout << "++++>>>> total num hits: " << hitPairVector.size()
1184  << ", num free: " << nFreePoints << std::endl;
1185  }
1186 
1188  {
1189  size_t localCount(0);
1190 
1191  if (!clusterParameters.daughterList().empty()) {
1192  for (auto& clusterParams : clusterParameters.daughterList())
1193  localCount += countUltimateDaughters(clusterParams);
1194  }
1195  else
1196  localCount++;
1197 
1198  return localCount;
1199  }
1200 
1202  ArtOutputHandler& output,
1203  reco::ClusterParameters& clusterParameters,
1204  size_t pfParticleParent,
1205  IdxToPCAMap& idxToPCAMap,
1206  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
1207  Hit3DToSPPtrMap& hit3DToSPPtrMap,
1208  Hit3DToSPPtrMap& best3DToSPPtrMap) const
1209  {
1210  // This is a recursive routine, we keep calling ourself as long as the daughter list is non empty
1211  if (!clusterParameters.daughterList().empty()) {
1212  for (auto& clusterParams : clusterParameters.daughterList())
1213  FindAndStoreDaughters(gser,
1214  output,
1215  clusterParams,
1216  pfParticleParent,
1217  idxToPCAMap,
1218  hitToPtrMap,
1219  hit3DToSPPtrMap,
1220  best3DToSPPtrMap);
1221  }
1222  // Otherwise we want to store the information
1223  else {
1224  size_t daughterIdx = ConvertToArtOutput(gser,
1225  output,
1226  clusterParameters,
1227  pfParticleParent,
1228  hitToPtrMap,
1229  hit3DToSPPtrMap,
1230  best3DToSPPtrMap);
1231 
1232  idxToPCAMap[daughterIdx] = &clusterParameters.getFullPCA();
1233  }
1234 
1235  return idxToPCAMap.size();
1236  }
1237 
1239  ArtOutputHandler& output,
1240  reco::ClusterParameters& clusterParameters,
1241  size_t pfParticleParent,
1242  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
1243  Hit3DToSPPtrMap& hit3DToSPPtrMap,
1244  Hit3DToSPPtrMap& best3DToSPPtrMap) const
1245  {
1246  // prepare the algorithm to compute the cluster characteristics;
1247  // we use the "standard" one here, except that we override selected items
1248  // (so, thanks to metaprogramming, we finally have wrappers of wrappers);
1249  // configuration would happen here, but we are using the default
1250  // configuration for that algorithm
1251  using OverriddenClusterParamsAlg_t =
1253 
1255 
1256  // It should be straightforward at this point to transfer information from our vector of clusters
1257  // to the larsoft objects... of course we still have some work to do first, in particular to
1258  // find the candidate seeds and their seed hits
1259 
1260  // We keep track of 2 PCA axes, the first is the "full" PCA run over all the 3D hits in the
1261  // candidate cluster. The second will be that derived from just using the "skeleton" hits.
1262  // Make a copy of the full PCA to keep that, then get a reference for the skeleton PCA
1263  reco::PrincipalComponents& fullPCA = clusterParameters.getFullPCA();
1264  reco::PrincipalComponents& skeletonPCA = clusterParameters.getSkeletonPCA();
1265 
1266  size_t clusterStart = output.artClusterVector->size();
1267 
1268  // Start loop over views to build out the hit lists and the 2D cluster objects
1269  for (const auto& clusParametersPair : clusterParameters.getClusterParams()) {
1270  const reco::RecobClusterParameters& clusParams = clusParametersPair.second;
1271 
1272  // Protect against a missing view
1273  if (clusParams.m_view == geo::kUnknown) continue;
1274 
1275  // Get a vector of hit pointers to seed the cluster algorithm
1276  std::vector<const recob::Hit*> recobHitVec(clusParams.m_hitVector.size());
1277 
1278  std::transform(clusParams.m_hitVector.begin(),
1279  clusParams.m_hitVector.end(),
1280  recobHitVec.begin(),
1281  [](const auto& hit2D) { return hit2D->getHit(); });
1282 
1283  // Get the tdc/wire slope... from the unit vector...
1284  double startWire(clusParams.m_startWire);
1285  double endWire(clusParams.m_endWire);
1286  double startTime(clusParams.m_startTime);
1287  double endTime(clusParams.m_endTime);
1288 
1289  // feed the algorithm with all the cluster hits
1290  ClusterParamAlgo.ImportHits(gser, recobHitVec);
1291 
1292  // create the recob::Cluster directly in the vector
1293  cluster::ClusterCreator artCluster(gser,
1294  ClusterParamAlgo, // algo
1295  startWire, // start_wire
1296  0., // sigma_start_wire
1297  startTime, // start_tick
1298  clusParams.m_sigmaStartTime, // sigma_start_tick
1299  endWire, // end_wire
1300  0., // sigma_end_wire,
1301  endTime, // end_tick
1302  clusParams.m_sigmaEndTime, // sigma_end_tick
1303  output.artClusterVector->size(), // ID
1304  clusParams.m_view, // view
1305  clusParams.m_plane, // plane
1306  recob::Cluster::Sentry // sentry
1307  );
1308 
1309  output.artClusterVector->emplace_back(artCluster.move());
1310 
1311  // We love looping. In this case, our list of hits is comprised of "ClusterHits" and we need to get a RecobHitVector instead...
1312  RecobHitVector recobHits;
1313 
1314  for (const auto& hit2D : clusParams.m_hitVector) {
1315  if (hit2D == nullptr || hit2D->getHit() == nullptr) continue;
1316 
1317  art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit2D->getHit()];
1318  recobHits.push_back(hitPtr);
1319  }
1320 
1321  output.makeClusterHitAssns(recobHits);
1322  }
1323 
1324  // Last, let's try to get seeds for tracking..
1325  // Keep track of how many we have so far
1326  size_t numSeedsStart = output.artSeedVector->size();
1327 
1328  // Empty daughter vector for now
1329  std::vector<size_t> nullVector;
1330  size_t pfParticleIdx(output.artPFParticleVector->size());
1331 
1332  recob::PFParticle pfParticle(13, pfParticleIdx, pfParticleParent, nullVector);
1333  output.artPFParticleVector->push_back(pfParticle);
1334 
1335  // Assume that if we are a daughter particle then there is a set of "best" 3D points and the convex hull will have been calculated from them
1336  // If there is no best list then assume the convex hull is from all of the points...
1337  std::string instance("");
1338  reco::HitPairListPtr* hit3DListPtr = &clusterParameters.getHitPairListPtr();
1339  Hit3DToSPPtrMap* hit3DToPtrMap = &hit3DToSPPtrMap;
1340 
1341  auto spVector = output.artSpacePointVector.get();
1342  auto edgeVector = output.artEdgeVector.get();
1343  auto pfPartEdgeAssns = output.artPFPartEdgeAssociations.get();
1344  auto spEdgeAssns = output.artEdgeSPAssociations.get();
1345  auto spHitAssns = output.artSPHitAssociations.get();
1346  auto pfPartSPAssns = output.artPFPartSPAssociations.get();
1347 
1348  // Make associations to all space points for this cluster
1349  int spacePointStart = spVector->size();
1350 
1352  output, *spVector, *spHitAssns, *hit3DListPtr, hitToPtrMap, *hit3DToPtrMap);
1353 
1354  // And make sure to have associations to the PFParticle
1355  output.makePFPartSpacePointAssns(*spVector, *pfPartSPAssns, spacePointStart);
1356 
1357  // If they exist, make space points for the Path points
1358  if (!clusterParameters.getBestHitPairListPtr().empty()) {
1359  instance = m_pathInstance;
1360  hit3DListPtr = &clusterParameters.getBestHitPairListPtr();
1361  hit3DToPtrMap = &best3DToSPPtrMap;
1362 
1363  spVector = output.artPathPointVector.get();
1364  edgeVector = output.artPathEdgeVector.get();
1365  pfPartEdgeAssns = output.artPFPartPathEdgeAssociations.get();
1366  spEdgeAssns = output.artEdgePPAssociations.get();
1367  spHitAssns = output.artPPHitAssociations.get();
1368  pfPartSPAssns = output.artPFPartPPAssociations.get();
1369 
1370  spacePointStart = spVector->size();
1371 
1373  output, *spVector, *spHitAssns, *hit3DListPtr, hitToPtrMap, *hit3DToPtrMap, instance);
1374 
1375  // And make sure to have associations to the PFParticle
1376  output.makePFPartSpacePointAssns(*spVector, *pfPartSPAssns, spacePointStart, instance);
1377  }
1378 
1379  // Now handle the edges according to whether associated with regular or "best" space points
1380  if (!clusterParameters.getBestEdgeList().empty()) {
1381  size_t edgeStart = edgeVector->size();
1382 
1383  for (const auto& edge : clusterParameters.getBestEdgeList()) {
1384  RecobSpacePointVector spacePointVec;
1385 
1386  try {
1387  spacePointVec.push_back(hit3DToPtrMap->at(std::get<0>(edge)));
1388  spacePointVec.push_back(hit3DToPtrMap->at(std::get<1>(edge)));
1389  }
1390  catch (...) {
1391  std::cout << "Caught exception in looking up space point ptr... " << std::get<0>(edge)
1392  << ", " << std::get<1>(edge) << std::endl;
1393  continue;
1394  }
1395 
1396  edgeVector->emplace_back(
1397  std::get<2>(edge), spacePointVec[0].key(), spacePointVec[1].key(), edgeVector->size());
1398 
1399  output.makeEdgeSpacePointAssns(*edgeVector, spacePointVec, *spEdgeAssns, instance);
1400  }
1401 
1402  output.makePFPartEdgeAssns(*edgeVector, *pfPartEdgeAssns, edgeStart, instance);
1403  }
1404 
1405  // Look at making the PCAxis and associations - for both the skeleton (the first) and the full
1406  // First need some float to double conversion containers
1407  recob::PCAxis::EigenVectors eigenVecs;
1408  double eigenVals[] = {0., 0., 0.};
1409  double avePosition[] = {0., 0., 0.};
1410 
1411  eigenVecs.resize(3);
1412 
1413  for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
1414  avePosition[outerIdx] = skeletonPCA.getAvePosition()[outerIdx];
1415  eigenVals[outerIdx] = skeletonPCA.getEigenValues()[outerIdx];
1416 
1417  eigenVecs[outerIdx].resize(3);
1418 
1419  for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
1420  eigenVecs[outerIdx][innerIdx] = skeletonPCA.getEigenVectors().row(outerIdx)(innerIdx);
1421  }
1422 
1423  recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
1424  skeletonPCA.getNumHitsUsed(),
1425  eigenVals,
1426  eigenVecs,
1427  avePosition,
1428  skeletonPCA.getAveHitDoca(),
1429  output.artPCAxisVector->size());
1430 
1431  output.artPCAxisVector->push_back(skelPcAxis);
1432 
1433  for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
1434  avePosition[outerIdx] = fullPCA.getAvePosition()[outerIdx];
1435  eigenVals[outerIdx] = fullPCA.getEigenValues()[outerIdx];
1436 
1437  for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
1438  eigenVecs[outerIdx][innerIdx] = fullPCA.getEigenVectors().row(outerIdx)(innerIdx);
1439  }
1440 
1441  recob::PCAxis fullPcAxis(fullPCA.getSvdOK(),
1442  fullPCA.getNumHitsUsed(),
1443  eigenVals, //fullPCA.getEigenValues(),
1444  eigenVecs, //fullPCA.getEigenVectors(),
1445  avePosition, //fullPCA.getAvePosition(),
1446  fullPCA.getAveHitDoca(),
1447  output.artPCAxisVector->size());
1448 
1449  output.artPCAxisVector->push_back(fullPcAxis);
1450 
1451  // Create associations to the PFParticle
1452  output.makePFPartPCAAssns();
1453  output.makePFPartSeedAssns(numSeedsStart);
1454  output.makePFPartClusterAssns(clusterStart);
1455 
1456  return pfParticleIdx;
1457  }
1458 
1460  std::vector<recob::SpacePoint>& spacePointVec,
1462  reco::HitPairListPtr& clusHitPairVector,
1463  IHit3DBuilder::RecobHitToPtrMap& hitToPtrMap,
1464  Hit3DToSPPtrMap& hit3DToSPPtrMap,
1465  const std::string& instance) const
1466  {
1467  // Reserve space...
1468  spacePointVec.reserve(spacePointVec.size() + clusHitPairVector.size());
1469 
1470  // Right now error matrix is uniform...
1471  double spError[] = {1., 0., 1., 0., 0., 1.};
1472 
1473  // Copy these hits to the vector to be stored with the event
1474  for (auto& hitPair : clusHitPairVector) {
1475  // Skip those space points that have already been created
1476  if (hit3DToSPPtrMap.find(hitPair) != hit3DToSPPtrMap.end()) continue;
1477 
1478  // Don't make space point if this hit was "rejected"
1479  if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
1480 
1481  double chisq = hitPair->getHitChiSquare(); // secret handshake...
1482 
1483  // Mark this hit pair as in use
1484  hitPair->setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);
1485 
1486  // Create and store the space point
1487  size_t spacePointID = spacePointVec.size();
1488  double spacePointPos[] = {
1489  hitPair->getPosition()[0], hitPair->getPosition()[1], hitPair->getPosition()[2]};
1490 
1491  spError[1] = hitPair->getTotalCharge();
1492  spError[3] = hitPair->getChargeAsymmetry();
1493 
1494  spacePointVec.emplace_back(spacePointPos, spError, chisq, spacePointID);
1495 
1496  // Update mapping
1497  hit3DToSPPtrMap[hitPair] = output.makeSpacePointPtr(spacePointID, instance);
1498 
1499  // space point hits associations
1500  RecobHitVector recobHits;
1501 
1502  for (const auto& hit : hitPair->getHits()) {
1503  if (!hit) continue;
1504 
1505  art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit->getHit()];
1506  recobHits.push_back(hitPtr);
1507  }
1508 
1509  if (!recobHits.empty())
1510  output.makeSpacePointHitAssns(spacePointVec, recobHits, spHitAssns, instance);
1511  }
1512 
1513  return;
1514  }
1515 
1517  reco::ConvexHullKinkTupleList& kinkTupleVec) const
1518  {
1519  // Right now error matrix is uniform...
1520  double spError[] = {1., 0., 1., 0., 0., 1.};
1521 
1522  // Copy these hits to the vector to be stored with the event
1523  for (auto& kinkTuple : kinkTupleVec) {
1524  const reco::ClusterHit3D* hit = std::get<2>(std::get<0>(kinkTuple));
1525 
1526  double chisq = hit->getHitChiSquare(); // secret handshake...
1527 
1528  // Create and store the space point
1529  double spacePointPos[] = {
1530  hit->getPosition()[0], hit->getPosition()[1], hit->getPosition()[2]};
1531 
1532  spError[1] = hit->getTotalCharge();
1533  spError[3] = hit->getChargeAsymmetry();
1534 
1535  output.artExtremePointVector->push_back(
1536  recob::SpacePoint(spacePointPos, spError, chisq, output.artExtremePointVector->size()));
1537  }
1538 
1539  return;
1540  }
1541 
1543  dcel2d::VertexList& vertexList,
1544  dcel2d::HalfEdgeList& halfEdgeList) const
1545  {
1546  // We actually do two things here:
1547  // 1) Create space points to represent the vertex locations of the voronoi diagram
1548  // 2) Create the edges that link the space points together
1549 
1550  // Set up the space point creation
1551  // Right now error matrix is uniform...
1552  double spError[] = {1., 0., 1., 0., 0., 1.};
1553  double chisq = 1.;
1554 
1555  // Keep track of the vertex to space point association
1556  std::map<const dcel2d::Vertex*, size_t> vertexToSpacePointMap;
1557 
1558  // Copy these hits to the vector to be stored with the event
1559  for (auto& vertex : vertexList) {
1560  const dcel2d::Coords& coords = vertex.getCoords();
1561 
1562  // Create and store the space point
1563  double spacePointPos[] = {coords[0], coords[1], coords[2]};
1564 
1565  vertexToSpacePointMap[&vertex] = output.artVertexPointVector->size();
1566 
1567  output.artVertexPointVector->emplace_back(
1568  spacePointPos, spError, chisq, output.artVertexPointVector->size());
1569  }
1570 
1571  // Try to avoid double counting
1572  std::set<const dcel2d::HalfEdge*> halfEdgeSet;
1573 
1574  // Build the edges now
1575  for (const auto& halfEdge : halfEdgeList) {
1576  // Recover twin
1577  const dcel2d::HalfEdge* twin = halfEdge.getTwinHalfEdge();
1578 
1579  // It can happen that we have no twin... and also check that we've not been here before
1580  if (twin && halfEdgeSet.find(twin) == halfEdgeSet.end()) {
1581  // Recover the vertices
1582  const dcel2d::Vertex* fromVertex = twin->getTargetVertex();
1583  const dcel2d::Vertex* toVertex = halfEdge.getTargetVertex();
1584 
1585  // It can happen for the open edges that there is no target vertex
1586  if (!toVertex || !fromVertex) continue;
1587 
1588  if (vertexToSpacePointMap.find(fromVertex) == vertexToSpacePointMap.end() ||
1589  vertexToSpacePointMap.find(toVertex) == vertexToSpacePointMap.end())
1590  continue;
1591 
1592  // Need the distance between vertices
1593  Eigen::Vector3f distVec = toVertex->getCoords() - fromVertex->getCoords();
1594 
1595  output.artVertexEdgeVector->emplace_back(distVec.norm(),
1596  vertexToSpacePointMap[fromVertex],
1597  vertexToSpacePointMap[toVertex],
1598  output.artEdgeVector->size());
1599 
1600  halfEdgeSet.insert(&halfEdge);
1601  }
1602  }
1603 
1604  return;
1605  }
1606 
1608  const reco::PrincipalComponents& fullPCA,
1609  IdxToPCAMap& idxToPCAMap) const
1610  {
1611  // We actually do two things here:
1612  // 1) Create space points from the centroids of the PCA for each cluster
1613  // 2) Create the edges that link the space points together
1614 
1615  // The first task is to put the list of PCA's into some semblance of order... they may be
1616  // preordered by likely they are piecewise ordered so fix that here
1617 
1618  // We'll need the current PCA axis to determine doca and arclen
1619  Eigen::Vector3f avePosition(
1620  fullPCA.getAvePosition()[0], fullPCA.getAvePosition()[1], fullPCA.getAvePosition()[2]);
1621  Eigen::Vector3f axisDirVec(fullPCA.getEigenVectors().row(2));
1622 
1623  using DocaToPCAPair = std::pair<float, const reco::PrincipalComponents*>;
1624  using DocaToPCAVec = std::vector<DocaToPCAPair>;
1625 
1626  DocaToPCAVec docaToPCAVec;
1627 
1628  // Outer loop over views
1629  for (const auto& idxToPCA : idxToPCAMap) {
1630  const reco::PrincipalComponents* pca = idxToPCA.second;
1631 
1632  // Now we need to calculate the doca and poca...
1633  // Start by getting this hits position
1634  Eigen::Vector3f pcaPos(
1635  pca->getAvePosition()[0], pca->getAvePosition()[1], pca->getAvePosition()[2]);
1636 
1637  // Form a TVector from this to the cluster average position
1638  Eigen::Vector3f pcaToAveVec = pcaPos - avePosition;
1639 
1640  // With this we can get the arclength to the doca point
1641  float arclenToPoca = pcaToAveVec.dot(axisDirVec);
1642 
1643  docaToPCAVec.emplace_back(arclenToPoca, pca);
1644  }
1645 
1646  std::sort(docaToPCAVec.begin(), docaToPCAVec.end(), [](const auto& left, const auto& right) {
1647  return left.first < right.first;
1648  });
1649 
1650  // Set up the space point creation
1651  // Right now error matrix is uniform...
1652  double spError[] = {1., 0., 1., 0., 0., 1.};
1653  double chisq = 1.;
1654 
1655  const reco::PrincipalComponents* lastPCA(nullptr);
1656 
1657  // Set up to loop through the clusters
1658  for (const auto& docaToPCAPair : docaToPCAVec) {
1659  // Recover the PCA for this cluster
1660  const reco::PrincipalComponents* curPCA = docaToPCAPair.second;
1661 
1662  if (lastPCA) {
1663  double lastPointPos[] = {
1664  lastPCA->getAvePosition()[0], lastPCA->getAvePosition()[1], lastPCA->getAvePosition()[2]};
1665  size_t lastPointBin = output.artVertexPointVector->size();
1666  double curPointPos[] = {
1667  curPCA->getAvePosition()[0], curPCA->getAvePosition()[1], curPCA->getAvePosition()[2]};
1668  size_t curPointBin = lastPointBin + 1;
1669 
1670  output.artVertexPointVector->emplace_back(lastPointPos, spError, chisq, lastPointBin);
1671  output.artVertexPointVector->emplace_back(curPointPos, spError, chisq, curPointBin);
1672 
1673  Eigen::Vector3f distVec(curPointPos[0] - lastPointPos[0],
1674  curPointPos[1] - lastPointPos[1],
1675  curPointPos[2] - lastPointPos[2]);
1676 
1677  output.artVertexEdgeVector->emplace_back(
1678  distVec.norm(), lastPointBin, curPointBin, output.artEdgeVector->size());
1679  }
1680 
1681  lastPCA = curPCA;
1682  }
1683  }
1684 
1685 } // namespace lar_cluster3d
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:467
intermediate_table::iterator iterator
std::unique_ptr< std::vector< recob::Edge > > artPathEdgeVector
float m_buildNeighborhoodTime
Keeps track of time to build epsilon neighborhood.
std::unique_ptr< art::Assns< recob::SpacePoint, recob::PFParticle > > artPFPartSPAssociations
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
std::list< reco::ClusterHit3D > HitPairList
Definition: Cluster3D.h:330
std::unique_ptr< std::vector< recob::SpacePoint > > artVertexPointVector
std::unique_ptr< lar_cluster3d::IClusterModAlg > m_clusterPathAlg
Algorithm to do cluster path finding.
bool getSvdOK() const
Definition: Cluster3D.h:240
void makeEdgeSpacePointAssns(std::vector< recob::Edge > &edgeVector, RecobSpacePointVector &spacePointVector, art::Assns< recob::SpacePoint, recob::Edge > &edgeSPAssociations, const std::string &path="")
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
Utilities related to art service access.
An object to define a "edge" which is used to connect space points in a triangulation algorithm...
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
Class managing the creation of a new recob::Cluster object.
std::unique_ptr< std::vector< recob::Seed > > artSeedVector
Reconstruction base classes.
T * get() const
Definition: ServiceHandle.h:69
std::string m_vertexInstance
Special instance name for vertex points.
static constexpr size_t kPFParticlePrimary
Define index to signify primary particle.
Definition: PFParticle.h:57
std::list< HalfEdge > HalfEdgeList
Definition: DCEL.h:171
void makePFPartSpacePointAssns(std::vector< recob::SpacePoint > &spacePointVector, art::Assns< recob::SpacePoint, recob::PFParticle > &pfPartSPAssociations, size_t spacePointStart, const std::string &instance="")
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
std::unique_ptr< art::Assns< recob::Seed, recob::Hit > > artSeedHitAssociations
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
float m_parallelHitsTransWid
Cut on transverse width of cluster (PCA 2nd eigenvalue)
HalfEdge * getTwinHalfEdge() const
Definition: DCEL.h:150
Unknown view.
Definition: geo_types.h:138
Float_t x1[n_points_granero]
Definition: compare.C:5
Cluster3D class.
Definition: SkeletonAlg.h:28
float m_parallelHitsCosAng
Cut for PCA 3rd axis angle to X axis.
bool aParallelHitsCluster(const reco::PrincipalComponents &pca) const
There are several places we will want to know if a candidate cluster is a "parallel hits" type of clu...
Declaration of signal hit object.
reco::PrincipalComponents & getSkeletonPCA()
Definition: Cluster3D.h:465
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
void findTrackSeeds(art::Event &evt, reco::ClusterParameters &cluster, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, std::vector< recob::Seed > &seedVec, art::Assns< recob::Seed, recob::Hit > &seedHitAssns) const
An interface to the seed finding algorithm.
bool findTrackSeeds(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitPairVec) const override
Given the list of hits this will search for candidate Seed objects and return them.
void makeClusterHitAssns(RecobHitVector &recobHits)
SkeletonAlg m_skeletonAlg
Skeleton point finder.
void makePFPartEdgeAssns(std::vector< recob::Edge > &edgeVector, art::Assns< recob::Edge, recob::PFParticle > &pfPartEdgeAssociations, size_t edgeStart, const std::string &instance="")
std::unique_ptr< std::vector< recob::Edge > > artVertexEdgeVector
const std::string instance
int m_hits3D
Keeps track of the number of 3D hits made.
STL namespace.
float getTotalCharge() const
Definition: Cluster3D.h:159
Hit has been rejected for any reason.
Definition: Cluster3D.h:96
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
int getNumHitsUsed() const
Definition: Cluster3D.h:241
HoughSeedFinderAlg class.
boost::graph_traits< ModuleGraph >::edge_descriptor Edge
Definition: ModuleGraph.h:24
std::list< HitPairListPtr > HitPairListPtrList
Definition: Cluster3D.h:328
A utility class used in construction of 3D clusters.
Definition: Cluster3D.h:292
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
void ProduceArtClusters(util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::HitPairList &hitPairList, reco::ClusterParametersList &clusterParametersList, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap) const
Top level output routine, allows checking cluster status.
bool m_onlyMakSpacePoints
If true we don&#39;t do the full cluster 3D processing.
void makePFPartSeedAssns(size_t numSeedsStart)
Cluster finding and building.
reco::PlaneToClusterParamsMap & getClusterParams()
Definition: Cluster3D.h:461
HoughSeedFinderAlg m_seedFinderAlg
Seed finder.
std::unique_ptr< lar_cluster3d::IClusterAlg > m_clusterAlg
Algorithm to do 3D space point clustering.
std::unique_ptr< art::Assns< recob::PFParticle, recob::Seed > > artPFPartSeedAssociations
float m_makeHitsTime
Keeps track of time to build 3D hits.
std::unique_ptr< std::vector< recob::PFParticle > > artPFParticleVector
void makeSpacePointHitAssns(std::vector< recob::SpacePoint > &spacePointVector, RecobHitVector &recobHits, art::Assns< recob::Hit, recob::SpacePoint > &spHitAssns, const std::string &path="")
void GetSkeletonHits(const reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &skeletonHitList) const
Return the skeleton hits from the input list.
art::Ptr< recob::Edge > makeEdgePtr(size_t index, const std::string &instance="")
This is an algorithm for finding recob::Seed objects in 3D clusters.
float getDocaToAxis() const
Definition: Cluster3D.h:166
static const SentryArgument_t Sentry
An instance of the sentry object.
Definition: Cluster.h:174
float m_dbscanTime
Keeps track of time to run DBScan.
void ImportHits(util::GeometryUtilities const &gser, Iter begin, Iter end)
Calls SetHits() with the hits in the sequence.
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
ClusterParametersList & daughterList()
Definition: Cluster3D.h:438
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
std::unique_ptr< art::Assns< recob::Hit, recob::SpacePoint > > artPPHitAssociations
float getAveHitDoca() const
Definition: Cluster3D.h:245
Overrides another ClusterParamsAlgBase class with selected constants.
This is an algorithm for finding recob::Seed objects in 3D clusters.
bool findTrackSeeds(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitMap) const override
Given the list of hits this will search for candidate Seed objects and return them.
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:242
art::PtrMaker< recob::Edge > fEdgePtrMaker
std::string m_extremeInstance
Instance name for the extreme points.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
size_t countUltimateDaughters(reco::ClusterParameters &clusterParameters) const
Count number of end of line daughters.
float getChargeAsymmetry() const
Definition: Cluster3D.h:165
art::PtrMaker< recob::SpacePoint > fSPPtrMakerPath
This is an algorithm for finding recob::Seed objects in 3D clusters.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
std::unique_ptr< std::vector< recob::SpacePoint > > artPathPointVector
void InitializeMonitoring()
Initialize the internal monitoring.
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
Definition: Cluster3D.h:499
float m_finishTime
Keeps track of time to run output module.
Helper functions to create a cluster.
Wrapper for ClusterParamsAlgBase objects to accept diverse input.
std::unique_ptr< std::vector< recob::SpacePoint > > artSpacePointVector
float m_clusterMergeTime
Keeps track of the time to merge clusters.
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
std::unique_ptr< art::Assns< recob::PFParticle, recob::PCAxis > > artPFPartAxisAssociations
std::unique_ptr< art::Assns< recob::Cluster, recob::Hit > > artClusterAssociations
PCASeedFinderAlg m_pcaSeedFinderAlg
Use PCA axis to find seeds.
ArtOutputHandler(art::Event &evt, std::string &pathName, std::string &vertexName, std::string &extremeName)
std::list< ConvexHullKinkTuple > ConvexHullKinkTupleList
Definition: Cluster3D.h:349
void AverageSkeletonPositions(reco::HitPairListPtr &skeletonHitList) const
Modifies the position of input skeleton hits by averaging along the "best" wire direction.
Wrapper for ClusterParamsAlgBase objects to accept arbitrary input.
bool findTrackSeeds(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitMap) const override
Given the list of hits this will search for candidate Seed objects and return them.
ParallelHitsSeedFinderAlg m_parallelHitsAlg
Deal with parallel hits clusters.
std::unique_ptr< std::vector< recob::Edge > > artEdgeVector
void MakeAndSavePCAPoints(ArtOutputHandler &, const reco::PrincipalComponents &, IdxToPCAMap &) const
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
std::unique_ptr< art::Assns< recob::Hit, recob::SpacePoint > > artSPHitAssociations
This provides an art tool interface definition for 3D Cluster algorithms.
std::unique_ptr< art::Assns< recob::SpacePoint, recob::Edge > > artEdgePPAssociations
int m_hits
Keeps track of the number of hits seen.
void produce(art::Event &evt) override
Declaration of cluster object.
size_t FindAndStoreDaughters(util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, IdxToPCAMap &idxToPCAMap, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, Hit3DToSPPtrMap &best3DToSPPtrMap) const
This will produce art output for daughters noting that it needs to be done recursively.
bool empty() const
Definition: PtrVector.h:330
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
bool operator()(const std::pair< float, const reco::ClusterHit3D * > &left, const std::pair< float, const reco::ClusterHit3D * > &right)
void MakeAndSaveSpacePoints(ArtOutputHandler &output, std::vector< recob::SpacePoint > &spacePointVec, art::Assns< recob::Hit, recob::SpacePoint > &spHitAssns, reco::HitPairListPtr &clusHitPairVector, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, const std::string &path="") const
Special routine to handle creating and saving space points.
This provides an art tool interface definition for tools which construct 3D hits used in 3D clusterin...
This header file defines the interface to a principal components analysis designed to be used within ...
std::unique_ptr< art::Assns< recob::SpacePoint, recob::PFParticle > > artPFPartPPAssociations
std::unique_ptr< art::Assns< recob::PFParticle, recob::Cluster > > artPFPartClusAssociations
ProducesCollector & producesCollector() noexcept
float m_pathFindingTime
Keeps track of the path finding time.
bool findTrackHits(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, reco::HitPairListPtrList &hitPairListPtrList) const
Given the list of hits this will return the sets of hits which belong on the same line...
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
bool operator()(const reco::ClusterHit3D *hit3D) const
art::Ptr< recob::SpacePoint > makeSpacePointPtr(size_t index, const std::string &instance="")
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > m_clusterBuilder
Common cluster builder tool.
Cluster3D(fhicl::ParameterSet const &pset)
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
art::PtrMaker< recob::SpacePoint > fSPPtrMaker
std::unique_ptr< std::vector< recob::SpacePoint > > artExtremePointVector
Hierarchical representation of particle flow.
Definition: PFParticle.h:44
Utility object to perform functions of association.
void makePFPartClusterAssns(size_t clusterStart)
Header file to define the interface to the SkeletonAlg.
Hit has been made into Space Point.
Definition: Cluster3D.h:100
std::vector< SeedHitPairListPair > SeedHitPairListPairVec
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
This provides an art tool interface definition for 3D Cluster algorithms.
void MakeAndSaveVertexPoints(ArtOutputHandler &, dcel2d::VertexList &, dcel2d::HalfEdgeList &) const
Special routine to handle creating and saving space points & edges associated to voronoi diagrams...
std::unique_ptr< std::vector< recob::Cluster > > artClusterVector
std::unique_ptr< lar_cluster3d::IClusterModAlg > m_clusterMergeAlg
Algorithm to do cluster merging.
size_t ConvertToArtOutput(util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, Hit3DToSPPtrMap &best3DToSPPtrMap) const
Produces the art output from all the work done in this producer module.
std::unique_ptr< lar_cluster3d::IHit3DBuilder > m_hit3DBuilderAlg
Builds the 3D hits to operate on.
PrincipalComponentsAlg m_pcaAlg
Principal Components algorithm.
Definition: MVAAlg.h:12
EventNumber_t event() const
Definition: EventID.h:116
art::PtrMaker< recob::Edge > fEdgePtrMakerPath
ClusterHit2DVec m_hitVector
Definition: Cluster3D.h:319
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:46
Algorithm collection class computing cluster parameters.
float getHitChiSquare() const
Definition: Cluster3D.h:163
Definition of the Cluster3D class.
float m_artHitsTime
Keeps track of time to recover hits.
Interface to class computing cluster parameters.
std::unique_ptr< std::vector< recob::PCAxis > > artPCAxisVector
TCEvent evt
Definition: DataStructs.cxx:8
Vertex * getTargetVertex() const
Definition: DCEL.h:148
Eigen::Vector3f Coords
Definition: DCEL.h:44
std::unique_ptr< art::Assns< recob::SpacePoint, recob::Edge > > artEdgeSPAssociations
std::unordered_map< const reco::ClusterHit3D *, art::Ptr< recob::SpacePoint >> Hit3DToSPPtrMap
RunNumber_t run() const
Definition: Event.cc:29
void MakeAndSaveKinkPoints(ArtOutputHandler &output, reco::ConvexHullKinkTupleList &clusHitPairVector) const
Special routine to handle creating and saving space points.
PCASeedFinderAlg class.
std::vector< std::vector< double > > EigenVectors
Definition: PCAxis.h:26
std::unordered_map< const recob::Hit *, art::Ptr< recob::Hit >> RecobHitToPtrMap
Defines a structure mapping art representation to internal.
Definition: IHit3DBuilder.h:43
void PrepareEvent(const art::Event &evt)
Event Preparation.
const Coords & getCoords() const
Definition: DCEL.h:62
std::list< Vertex > VertexList
Definition: DCEL.h:169
void splitClustersWithHough(reco::ClusterParameters &clusterParameters, reco::ClusterParametersList &clusterParametersList) const
Attempt to split clusters using the output of the Hough Filter.
EventID id() const
Definition: Event.cc:23
std::string m_pathInstance
Special instance for path points.
art framework interface to geometry description
float m_totalTime
Keeps track of total execution time.
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:393
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
std::unique_ptr< art::Assns< recob::Edge, recob::PFParticle > > artPFPartEdgeAssociations
std::unique_ptr< art::Assns< recob::Edge, recob::PFParticle > > artPFPartPathEdgeAssociations
std::map< size_t, const reco::PrincipalComponents * > IdxToPCAMap
Special routine to handle creating and saving space points & edges PCA points.
vertex reconstruction
bool m_enableMonitoring
Turn on monitoring of this algorithm.