LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
lar_cluster3d::PCASeedFinderAlg Class Referencefinal

PCASeedFinderAlg class. More...

#include "PCASeedFinderAlg.h"

Inheritance diagram for lar_cluster3d::PCASeedFinderAlg:
lar_cluster3d::SeedFinderAlgBase

Public Member Functions

 PCASeedFinderAlg (fhicl::ParameterSet const &pset)
 Constructor. More...
 
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. More...
 

Private Member Functions

bool getHitsAtEnd (reco::HitPairListPtr &hit3DList, reco::PrincipalComponents &seedPca) const
 Separate function to find hits at the ends of the input hits. More...
 
void LineFit2DHits (const reco::HitPairListPtr &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
 

Private Attributes

geo::Geometry const * m_geometry
 
double m_gapDistance
 
size_t m_numSeed2DHits
 
double m_minAllowedCosAng
 The minimum cos(ang) between input and seed axes. More...
 
PrincipalComponentsAlg m_pcaAlg
 

Detailed Description

PCASeedFinderAlg class.

Definition at line 33 of file PCASeedFinderAlg.h.

Constructor & Destructor Documentation

lar_cluster3d::PCASeedFinderAlg::PCASeedFinderAlg ( fhicl::ParameterSet const &  pset)

Constructor.

Parameters
pset

Definition at line 36 of file PCASeedFinderAlg.cxx.

References art::ServiceHandle< T, SCOPE >::get(), fhicl::ParameterSet::get(), m_gapDistance, m_geometry, m_minAllowedCosAng, and m_numSeed2DHits.

37  : m_gapDistance(5.)
38  , m_numSeed2DHits(80)
39  , m_minAllowedCosAng(0.7)
40  , m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
41  {
42  m_gapDistance = pset.get<double>("GapDistance", 5.);
43  m_numSeed2DHits = pset.get<size_t>("NumSeed2DHits", 80);
44  m_minAllowedCosAng = pset.get<double>("MinAllowedCosAng", 0.7);
46  }
geo::Geometry const * m_geometry
T * get() const
Definition: ServiceHandle.h:69
PrincipalComponentsAlg m_pcaAlg
double m_minAllowedCosAng
The minimum cos(ang) between input and seed axes.

Member Function Documentation

bool lar_cluster3d::PCASeedFinderAlg::findTrackSeeds ( reco::HitPairListPtr hitPairListPtr,
reco::PrincipalComponents inputPCA,
SeedHitPairListPairVec seedHitMap 
) const
overridevirtual

Given the list of hits this will search for candidate Seed objects and return them.

Implements lar_cluster3d::SeedFinderAlgBase.

Definition at line 48 of file PCASeedFinderAlg.cxx.

References reco::PrincipalComponents::flipAxis(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), getHitsAtEnd(), LineFit2DHits(), m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas().

Referenced by lar_cluster3d::Cluster3D::findTrackSeeds().

51  {
52  bool foundGoodSeed(false);
53 
54  // Make sure we are using the right pca
55  reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
56 
57  // Make a local copy of the input pca
58  reco::PrincipalComponents pca = inputPCA;
59 
60  // We also require that there be some spread in the data, otherwise not worth running?
61  double eigenVal0 = 3. * sqrt(pca.getEigenValues()[2]);
62  double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
63 
64  if (eigenVal0 > 5. && eigenVal1 > 0.001) {
65  // Presume CR muons will be "downward going"...
66  if (pca.getEigenVectors().row(2)(1) > 0.) pca.flipAxis(0);
67 
68  // Use the following to set the 3D doca and arclength for each hit
69  m_pcaAlg.PCAAnalysis_calc3DDocas(hitPairListPtr, pca);
70 
71  // Use this info to sort the hits along the principle axis
72  hitPairListPtr.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
73  //hitPairListPtr.sort(PcaSort3DHitsByAbsArcLen3D());
74 
75  // Make a local copy of the 3D hits
76  reco::HitPairListPtr hit3DList;
77 
78  hit3DList.resize(hitPairListPtr.size());
79 
80  std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), hit3DList.begin());
81 
82  reco::PrincipalComponents seedPca = pca;
83 
84  if (getHitsAtEnd(hit3DList, seedPca)) {
85  // We can use the same routine to check hits at the opposite end to make sure
86  // we have consistency between both ends of the track.
87  // So... follow the same general program as above
88  reco::HitPairListPtr checkList;
89 
90  checkList.resize(hitPairListPtr.size());
91 
92  std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), checkList.begin());
93 
94  std::reverse(checkList.begin(), checkList.end());
95 
96  reco::PrincipalComponents checkPca = pca;
97 
98  if (getHitsAtEnd(checkList, checkPca)) {
99  // We want our seed to be in the middle of our seed points
100  // First step to getting there is to reorder the hits along the
101  // new pca axis
102  m_pcaAlg.PCAAnalysis_calc3DDocas(hit3DList, seedPca);
103 
104  // Use this info to sort the hits along the principle axis - note by absolute value of arc length
105  hit3DList.sort(SeedFinderAlgBase::Sort3DHitsByAbsArcLen3D());
106 
107  // Now translate the seedCenter by the arc len to the first hit
108  double seedDir[3] = {seedPca.getEigenVectors().row(2)(0),
109  seedPca.getEigenVectors().row(2)(1),
110  seedPca.getEigenVectors().row(2)(2)};
111  double seedStart[3] = {
112  hit3DList.front()->getX(), hit3DList.front()->getY(), hit3DList.front()->getZ()};
113 
114  if (hit3DList.size() > 10) {
115  TVector3 newSeedPos;
116  TVector3 newSeedDir;
117  double chiDOF;
118 
119  LineFit2DHits(hit3DList, seedStart[0], newSeedPos, newSeedDir, chiDOF);
120 
121  if (chiDOF > 0.) {
122  // check angles between new/old directions
123  double cosAng = seedDir[0] * newSeedDir[0] + seedDir[1] * newSeedDir[1] +
124  seedDir[2] * newSeedDir[2];
125 
126  if (cosAng < 0.) newSeedDir *= -1.;
127 
128  seedStart[0] = newSeedPos[0];
129  seedStart[1] = newSeedPos[1];
130  seedStart[2] = newSeedPos[2];
131  seedDir[0] = newSeedDir[0];
132  seedDir[1] = newSeedDir[1];
133  seedDir[2] = newSeedDir[2];
134  }
135  }
136 
137  for (const auto& hit3D : hit3DList)
138  hit3D->setStatusBit(0x40000000);
139 
140  seedHitPairVec.emplace_back(std::pair<recob::Seed, reco::HitPairListPtr>(
141  recob::Seed(seedStart, seedDir), hit3DList));
142 
143  foundGoodSeed = true;
144  }
145  }
146  }
147 
148  return foundGoodSeed;
149  }
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
bool getHitsAtEnd(reco::HitPairListPtr &hit3DList, reco::PrincipalComponents &seedPca) const
Separate function to find hits at the ends of the input hits.
void flipAxis(size_t axis)
Definition: Cluster3D.cxx:232
PrincipalComponentsAlg m_pcaAlg
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:242
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
void LineFit2DHits(const reco::HitPairListPtr &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
bool lar_cluster3d::PCASeedFinderAlg::getHitsAtEnd ( reco::HitPairListPtr hit3DList,
reco::PrincipalComponents seedPca 
) const
private

Separate function to find hits at the ends of the input hits.

Definition at line 151 of file PCASeedFinderAlg.cxx.

References reco::PrincipalComponents::flipAxis(), reco::ClusterHit3D::getArclenToPoca(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterHit3D::getHits(), reco::PrincipalComponents::getSvdOK(), m_gapDistance, m_minAllowedCosAng, m_numSeed2DHits, m_pcaAlg, and lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D().

Referenced by findTrackSeeds().

153  {
154  bool foundGoodSeedHits(false);
155 
156  // Use a set to count 2D hits by keeping only a single copy
157  std::set<const reco::ClusterHit2D*> hit2DSet;
158 
159  // Look for numSeed2DHits which are continuous
160  double lastArcLen = hit3DList.front()->getArclenToPoca();
161 
162  reco::HitPairListPtr::iterator startItr = hit3DList.begin();
163  reco::HitPairListPtr::iterator lastItr = hit3DList.begin();
164 
165  while (++lastItr != hit3DList.end()) {
166  const reco::ClusterHit3D* hit3D = *lastItr;
167  double arcLen = hit3D->getArclenToPoca();
168 
169  if (fabs(arcLen - lastArcLen) > m_gapDistance) {
170  startItr = lastItr;
171  hit2DSet.clear();
172  }
173 
174  for (const auto& hit2D : hit3D->getHits()) {
175  hit2DSet.insert(hit2D);
176  }
177 
178  if (hit2DSet.size() > m_numSeed2DHits) break;
179 
180  lastArcLen = arcLen;
181  }
182 
183  if (hit2DSet.size() > m_numSeed2DHits) {
184  if (startItr != hit3DList.begin()) hit3DList.erase(hit3DList.begin(), startItr);
185  if (lastItr != hit3DList.end()) hit3DList.erase(lastItr, hit3DList.end());
186 
187  // On input, the seedPca will contain the original values so we can recover the original axis now
188  Eigen::Vector3f planeVec0(seedPca.getEigenVectors().row(2));
189 
190  m_pcaAlg.PCAAnalysis_3D(hit3DList, seedPca, true);
191 
192  if (seedPca.getSvdOK()) {
193  // Still looking to point "down"
194  if (seedPca.getEigenVectors().row(2)(1) > 0.) seedPca.flipAxis(0);
195 
196  // Check that the seed PCA we have found is consistent with the input PCA
197  Eigen::Vector3f primarySeedAxis(seedPca.getEigenVectors().row(2));
198 
199  double cosAng = primarySeedAxis.dot(planeVec0);
200 
201  // If the proposed seed axis is not relatively aligned with the input PCA then
202  // we should not be using this method to return seeds. Check that here
203  if (cosAng > m_minAllowedCosAng) foundGoodSeedHits = true;
204  }
205  }
206 
207  return foundGoodSeedHits;
208  }
intermediate_table::iterator iterator
bool getSvdOK() const
Definition: Cluster3D.h:240
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
void flipAxis(size_t axis)
Definition: Cluster3D.cxx:232
PrincipalComponentsAlg m_pcaAlg
double m_minAllowedCosAng
The minimum cos(ang) between input and seed axes.
float getArclenToPoca() const
Definition: Cluster3D.h:167
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:168
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
void lar_cluster3d::PCASeedFinderAlg::LineFit2DHits ( const reco::HitPairListPtr hitList,
double  XOrigin,
TVector3 &  Pos,
TVector3 &  Dir,
double &  ChiDOF 
) const
private

Definition at line 211 of file PCASeedFinderAlg.cxx.

References m_geometry, norm, geo::PlaneID::Plane, sw, w, geo::WireID::Wire, geo::GeometryCore::WireCoordinate(), geo::GeometryCore::WirePitch(), and x.

Referenced by findTrackSeeds().

216  {
217  // The following is lifted from Bruce Baller to try to get better
218  // initial parameters for a candidate Seed. It is slightly reworked
219  // which is why it is included here instead of used as is.
220  //
221  // Linear fit using X as the independent variable. Hits to be fitted
222  // are passed in the hits vector in a pair form (X, WireID). The
223  // fitted track position at XOrigin is returned in the Pos vector.
224  // The direction cosines are returned in the Dir vector.
225  //
226  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
227  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
228  // a matrix to calculate a track projected to a point at X, and v is
229  // a vector (Yo, Zo, dY/dX, dZ/dX).
230  //
231  // Note: The covariance matrix should also be returned
232  // B. Baller August 2014
233 
234  // assume failure
235  ChiDOF = -1;
236 
237  // Make a set of 2D hits so we consider only unique hits
238  std::set<const reco::ClusterHit2D*> hit2DSet;
239 
240  for (const auto& hit3D : hit3DList) {
241  for (const auto& hit2D : hit3D->getHits())
242  hit2DSet.insert(hit2D);
243  }
244 
245  if (hit2DSet.size() < 4) return;
246 
247  const unsigned int nvars = 4;
248  unsigned int npts = 3 * hit3DList.size();
249 
250  TMatrixD A(npts, nvars);
251  // vector holding the Wire number
252  TVectorD w(npts);
253  unsigned short ninpl[3] = {0};
254  unsigned short nok = 0;
255  unsigned short iht(0);
256 
257  // Loop over the 2D hits in the above
258  for (const auto& hit : hit2DSet) {
259  geo::WireID const& wireID = hit->WireID();
260  unsigned int ipl = wireID.Plane;
261 
262  // get the wire plane offset
263  double const off = m_geometry->WireCoordinate(geo::Point_t{0, 0, 0}, wireID);
264 
265  // get the "cosine-like" component
266  double const cw = m_geometry->WireCoordinate(geo::Point_t{0, 1, 0}, wireID) - off;
267 
268  // the "sine-like" component
269  double const sw = m_geometry->WireCoordinate(geo::Point_t{0, 0, 1}, wireID) - off;
270 
271  double const x = hit->getXPosition() - XOrigin;
272 
273  A[iht][0] = cw;
274  A[iht][1] = sw;
275  A[iht][2] = cw * x;
276  A[iht][3] = sw * x;
277  w[iht] = wireID.Wire - off;
278 
279  ++ninpl[ipl];
280 
281  // need at least two points in a plane
282  if (ninpl[ipl] == 2) ++nok;
283 
284  iht++;
285  }
286 
287  // need at least 2 planes with at least two points
288  if (nok < 2) return;
289 
290  TDecompSVD svd(A);
291  bool ok;
292  TVectorD tVec = svd.Solve(w, ok);
293 
294  ChiDOF = 0;
295 
296  // not enough points to calculate Chisq
297  if (npts <= 4) return;
298 
299  for (const auto& hit : hit2DSet) {
300  geo::WireID const& wireID = hit->WireID();
301  double const off = m_geometry->WireCoordinate(geo::Point_t{0, 0, 0}, wireID);
302  double const cw = m_geometry->WireCoordinate(geo::Point_t{0, 1, 0}, wireID) - off;
303  double const sw = m_geometry->WireCoordinate(geo::Point_t{0, 0, 1}, wireID) - off;
304  double const x = hit->getXPosition() - XOrigin;
305  double const ypr = tVec[0] + tVec[2] * x;
306  double const zpr = tVec[1] + tVec[3] * x;
307  double const diff = ypr * cw + zpr * sw - (wireID.Wire - off);
308  ChiDOF += diff * diff;
309  }
310 
311  float werr2 = m_geometry->WirePitch() * m_geometry->WirePitch();
312  ChiDOF /= werr2;
313  ChiDOF /= (float)(npts - 4);
314 
315  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
316  Dir[0] = 1 / norm;
317  Dir[1] = tVec[2] / norm;
318  Dir[2] = tVec[3] / norm;
319 
320  Pos[0] = XOrigin;
321  Pos[1] = tVec[0];
322  Pos[2] = tVec[1];
323 
324  } // TrkLineFit()
Float_t x
Definition: compare.C:6
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
geo::Geometry const * m_geometry
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Detector simulation of raw signals on wires.
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:180
Float_t sw
Definition: plot.C:20
Float_t norm
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.

Member Data Documentation

double lar_cluster3d::PCASeedFinderAlg::m_gapDistance
private

Definition at line 63 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and PCASeedFinderAlg().

geo::Geometry const* lar_cluster3d::PCASeedFinderAlg::m_geometry
private

Definition at line 61 of file PCASeedFinderAlg.h.

Referenced by LineFit2DHits(), and PCASeedFinderAlg().

double lar_cluster3d::PCASeedFinderAlg::m_minAllowedCosAng
private

The minimum cos(ang) between input and seed axes.

Definition at line 65 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and PCASeedFinderAlg().

size_t lar_cluster3d::PCASeedFinderAlg::m_numSeed2DHits
private

Definition at line 64 of file PCASeedFinderAlg.h.

Referenced by getHitsAtEnd(), and PCASeedFinderAlg().

PrincipalComponentsAlg lar_cluster3d::PCASeedFinderAlg::m_pcaAlg
private

Definition at line 67 of file PCASeedFinderAlg.h.

Referenced by findTrackSeeds(), and getHitsAtEnd().


The documentation for this class was generated from the following files: