ParallelHitsSeedFinderAlg.cxx
Go to the documentation of this file.
1 /**
2  * @file ParallelHitsSeedFinderAlg.cxx
3  *
4  * @brief Implementation of the Seed Finder Algorithm
5  * The intent of this algorithm is to take an input list of 3D space points and from those
6  * to find candidate track start points and directions
7  */
8 
9 // The main include
11 // Framework Includes
12 
13 // LArSoft includes
17 #include "lardata/RecoObjects/Cluster3D.h"
20 
21 // ROOT includes
22 #include "TTree.h"
23 #include "TVectorD.h"
24 #include "TMatrixD.h"
25 #include "TDecompSVD.h"
26 
27 // std includes
28 #include <string>
29 #include <functional>
30 #include <iostream>
31 #include <memory>
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
35 //------------------------------------------------------------------------------------------------------------------------------------------
36 // implementation follows
37 
38 namespace lar_cluster3d {
39 
41  m_maxNumEdgeHits(1000),
42  m_gapDistance(20.),
43  m_numSeed2DHits(80),
44  m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
45 {
46  this->reconfigure(pset);
47 
49 
50  m_geometry = &*geometry;
51  // m_detector = detectorProperties->provider();
52 }
53 
54 //------------------------------------------------------------------------------------------------------------------------------------------
55 
57 {
58 }
59 
61 {
62  m_maxNumEdgeHits = pset.get<size_t>("MaxNumEdgeHits", 1000);
63  m_gapDistance = pset.get<double>("GapDistance", 20.);
64  m_numSeed2DHits = pset.get<size_t>("NumSeed2DHits", 80);
65  m_pcaAlg.reconfigure(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
66 }
67 
69  reco::PrincipalComponents& inputPCA,
70  SeedHitPairListPairVec& seedHitPairVec) const
71 {
72  // This routine can fail...
73  bool foundGoodSeeds(false);
74 
75  // Make sure we are using the right pca
76  reco::PrincipalComponents pca = inputPCA;
77 
78  if (pca.getSvdOK())
79  {
80  // Presume CR muons will be "downward going"...
81  if (pca.getEigenVectors()[0][1] > 0.) pca.flipAxis(0);
82 
83  // This routine is typically called when there are LOTS of hits... so we are going to try
84  // to reduce the number of operations on the full list as much as possible. However, we
85  // really need the input hit list to be sorted by th input PCA so do that now
86  m_pcaAlg.PCAAnalysis_calc3DDocas(inputHitPairListPtr, pca);
87 
88  // Use this info to sort the hits along the principle axis
89  // Note that this will sort hits from the "middle" to the "outside"
90  inputHitPairListPtr.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
91 
92  // We need to determine the number of hits to drop... and this is really driven by the number of unique
93  // 2D hits we want to keep at one end of the cluster. So, make some containers and do some looping
94  reco::HitPairListPtr seedHit3DList;
95  std::set<const reco::ClusterHit2D*> hit2DSet;
96 
97  // Look for numSeed2DHits which are continuous
98  double lastArcLen = inputHitPairListPtr.front()->getArclenToPoca();
99 
100  for(const auto& hit3D : inputHitPairListPtr)
101  {
102  double arcLen = hit3D->getArclenToPoca();
103 
104  if (fabs(arcLen-lastArcLen) > m_gapDistance)
105  {
106  seedHit3DList.clear();
107  hit2DSet.clear();
108  }
109 
110  seedHit3DList.push_back(hit3D);
111 
112  for(const auto& hit2D : hit3D->getHits())
113  {
114  hit2DSet.insert(hit2D);
115  }
116 
117  if (hit2DSet.size() > m_numSeed2DHits) break;
118 
119  lastArcLen = arcLen;
120  }
121 
122  // We require a minimum number of seed hits in order to proceed
123  if (hit2DSet.size() > m_numSeed2DHits)
124  {
125  // The above derivation determines the number of hits to keep each side of the cloud
126  size_t num3DHitsToKeep = std::min(2*seedHit3DList.size(), inputHitPairListPtr.size());
127  size_t numEdgeHits = std::min(size_t(num3DHitsToKeep/2), m_maxNumEdgeHits);
128 
129  // Get an iterator for the hits
130  reco::HitPairListPtr::iterator edgeHitItr = inputHitPairListPtr.begin();
131 
132  std::advance(edgeHitItr, numEdgeHits);
133 
134  // Make a container for copying in the edge hits and size it to hold all the hits
135  reco::HitPairListPtr hit3DList;
136 
137  hit3DList.resize(2*numEdgeHits);
138 
139  // Copy the low edge hit range into the list
140  reco::HitPairListPtr::iterator nextHit3DItr = std::copy(inputHitPairListPtr.begin(), edgeHitItr, hit3DList.begin());
141 
142  // reset the "seed hits"
143  seedHit3DList.clear();
144  seedHit3DList.resize(numEdgeHits);
145 
146  std::copy(inputHitPairListPtr.begin(), edgeHitItr, seedHit3DList.begin());
147 
148  // Now advance the iterator into the main container and copy the rest of the elements
149  std::advance(edgeHitItr, inputHitPairListPtr.size() - 2 * numEdgeHits);
150 
151  std::copy(edgeHitItr, inputHitPairListPtr.end(), nextHit3DItr);
152 
153  // At this point we should now have trimmed out the bulk of the 3D hits from the middle
154  // of the input cloud of hits. Next step is to rerun the PCA on our reduced set of hits
156 
157  m_pcaAlg.PCAAnalysis_3D(hit3DList, seedPCA, false);
158 
159  if (seedPCA.getSvdOK())
160  {
161  // Still looking to point "down"
162  if (seedPCA.getEigenVectors()[0][1] > 0.) seedPCA.flipAxis(0);
163 
164  // Recompute the 3D docas and arc lengths
165  m_pcaAlg.PCAAnalysis_calc3DDocas(hit3DList, seedPCA);
166 
167  // Now sort hits in regular order
168  //hit3DList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
169  seedHit3DList.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
170 
171  // Now translate the seedCenter by the arc len to the first hit
172  double seedDir[3] = {seedPCA.getEigenVectors()[0][0], seedPCA.getEigenVectors()[0][1], seedPCA.getEigenVectors()[0][2]};
173  double seedStart[3] = {seedHit3DList.front()->getX(), seedHit3DList.front()->getY(), seedHit3DList.front()->getZ()};
174 
175  // Move from the first hit to the halfway point but from the first hit, not from the axis position
176  double halfArcLen = 0.5 * fabs(seedHit3DList.back()->getArclenToPoca() - seedHit3DList.front()->getArclenToPoca());
177 
178  seedStart[0] += halfArcLen * seedDir[0];
179  seedStart[1] += halfArcLen * seedDir[1];
180  seedStart[2] += halfArcLen * seedDir[2];
181 
182  for(const auto& hit3D : seedHit3DList) hit3D->setStatusBit(0x40000000);
183 
184  seedHitPairVec.emplace_back(std::pair<recob::Seed, reco::HitPairListPtr>(recob::Seed(seedStart, seedDir), seedHit3DList));
185 
186  inputPCA = seedPCA;
187 
188  foundGoodSeeds = true;
189  }
190  }
191  }
192 
193  return foundGoodSeeds;
194 }
195 
196 } // namespace lar_cluster3d
intermediate_table::iterator iterator
Define a comparator which will sort hits by arc length along a PCA axis.
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
bool getSvdOK() const
Definition: Cluster3D.h:244
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
double m_gapDistance
Maximum allowed distance between hits.
void flipAxis(size_t axis)
Definition: Cluster3D.cxx:208
virtual void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
size_t m_maxNumEdgeHits
Maximum number hits each end of PCA axis.
art framework interface to geometry description
ParallelHitsSeedFinderAlg(fhicl::ParameterSet const &pset)
Constructor.
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
Encapsulate the geometry of a wire.
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
Encapsulate the construction of a single detector plane.
size_t m_numSeed2DHits
Number 2D seed hits desired.
T copy(T const &v)
std::vector< SeedHitPairListPair > SeedHitPairListPairVec
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
virtual bool findTrackSeeds(reco::HitPairListPtr &hitPairListPtr, reco::PrincipalComponents &inputPCA, SeedHitPairListPairVec &seedHitMap) const
Given the list of hits this will search for candidate Seed objects and return them.
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247