LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
PCASeedFinderAlg.cxx
Go to the documentation of this file.
1 
10 // The main include
12 // Framework Includes
13 
14 // LArSoft includes
20 
21 // ROOT includes
22 #include "TTree.h"
23 #include "TCanvas.h"
24 #include "TFrame.h"
25 #include "TH2D.h"
26 #include "TVectorD.h"
27 #include "TMatrixD.h"
28 #include "TDecompSVD.h"
29 
30 // std includes
31 #include <string>
32 #include <functional>
33 #include <iostream>
34 #include <memory>
35 
36 //------------------------------------------------------------------------------------------------------------------------------------------
37 
38 //------------------------------------------------------------------------------------------------------------------------------------------
39 // implementation follows
40 
41 namespace lar_cluster3d {
42 
44  m_gapDistance(5.),
45  m_numSeed2DHits(80),
46  m_minAllowedCosAng(0.7),
47  m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
48 {
49  this->reconfigure(pset);
50 
52  // auto const* detectorProperties = lar::providerFrom<detinfo::DetectorPropertiesService>();
53 
54  m_geometry = &*geometry;
55  // m_detector = detectorProperties->provider();
56 }
57 
58 //------------------------------------------------------------------------------------------------------------------------------------------
59 
61 {
62 }
63 
65 {
66  m_gapDistance = pset.get<double>("GapDistance", 5.);
67  m_numSeed2DHits = pset.get<size_t>("NumSeed2DHits", 80);
68  m_minAllowedCosAng = pset.get<double>("MinAllowedCosAng", 0.7);
69  m_pcaAlg.reconfigure(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
70 }
71 
73  reco::PrincipalComponents& inputPCA,
74  SeedHitPairListPairVec& seedHitPairVec) const
75 {
76  bool foundGoodSeed(false);
77 
78  // Make sure we are using the right pca
79  reco::HitPairListPtr hitPairListPtr = inputHitPairListPtr;
80 
81  // Make a local copy of the input pca
82  reco::PrincipalComponents pca = inputPCA;
83 
84  // We also require that there be some spread in the data, otherwise not worth running?
85  double eigenVal0 = 3. * sqrt(pca.getEigenValues()[0]);
86  double eigenVal1 = 3. * sqrt(pca.getEigenValues()[1]);
87 
88  if (eigenVal0 > 5. && eigenVal1 > 0.001)
89  {
90  // Presume CR muons will be "downward going"...
91  if (pca.getEigenVectors()[0][1] > 0.) pca.flipAxis(0);
92 
93  // Use the following to set the 3D doca and arclength for each hit
94  m_pcaAlg.PCAAnalysis_calc3DDocas(hitPairListPtr, pca);
95 
96  // Use this info to sort the hits along the principle axis
97  hitPairListPtr.sort(SeedFinderAlgBase::Sort3DHitsByArcLen3D());
98  //hitPairListPtr.sort(PcaSort3DHitsByAbsArcLen3D());
99 
100  // Make a local copy of the 3D hits
101  reco::HitPairListPtr hit3DList;
102 
103  hit3DList.resize(hitPairListPtr.size());
104 
105  std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), hit3DList.begin());
106 
107  reco::PrincipalComponents seedPca = pca;
108 
109  if (getHitsAtEnd(hit3DList, seedPca))
110  {
111  // We can use the same routine to check hits at the opposite end to make sure
112  // we have consistency between both ends of the track.
113  // So... follow the same general program as above
114  reco::HitPairListPtr checkList;
115 
116  checkList.resize(hitPairListPtr.size());
117 
118  std::copy(hitPairListPtr.begin(), hitPairListPtr.end(), checkList.begin());
119 
120  std::reverse(checkList.begin(), checkList.end());
121 
122  reco::PrincipalComponents checkPca = pca;
123 
124  if (getHitsAtEnd(checkList, checkPca))
125  {
126  // We want our seed to be in the middle of our seed points
127  // First step to getting there is to reorder the hits along the
128  // new pca axis
129  m_pcaAlg.PCAAnalysis_calc3DDocas(hit3DList, seedPca);
130 
131  // Use this info to sort the hits along the principle axis - note by absolute value of arc length
133 
134  // Now translate the seedCenter by the arc len to the first hit
135  double seedDir[3] = {seedPca.getEigenVectors()[0][0], seedPca.getEigenVectors()[0][1], seedPca.getEigenVectors()[0][2]};
136  double seedStart[3] = {hit3DList.front()->getX(), hit3DList.front()->getY(), hit3DList.front()->getZ()};
137 
138  if (hit3DList.size() > 10)
139  {
140  TVector3 newSeedPos;
141  TVector3 newSeedDir;
142  double chiDOF;
143 
144  LineFit2DHits(hit3DList, seedStart[0], newSeedPos, newSeedDir, chiDOF);
145 
146  if (chiDOF > 0.)
147  {
148  // check angles between new/old directions
149  double cosAng = seedDir[0]*newSeedDir[0]+seedDir[1]*newSeedDir[1]+seedDir[2]*newSeedDir[2];
150 
151  if (cosAng < 0.) newSeedDir *= -1.;
152 
153  seedStart[0] = newSeedPos[0];
154  seedStart[1] = newSeedPos[1];
155  seedStart[2] = newSeedPos[2];
156  seedDir[0] = newSeedDir[0];
157  seedDir[1] = newSeedDir[1];
158  seedDir[2] = newSeedDir[2];
159  }
160  }
161 
162  for(const auto& hit3D : hit3DList) hit3D->setStatusBit(0x40000000);
163 
164  seedHitPairVec.emplace_back(std::pair<recob::Seed, reco::HitPairListPtr>(recob::Seed(seedStart, seedDir), hit3DList));
165 
166  foundGoodSeed = true;
167  }
168  }
169  }
170 
171  return foundGoodSeed;
172 }
173 
175 {
176  bool foundGoodSeedHits(false);
177 
178  // Use a set to count 2D hits by keeping only a single copy
179  std::set<const reco::ClusterHit2D*> hit2DSet;
180 
181  // Look for numSeed2DHits which are continuous
182  double lastArcLen = hit3DList.front()->getArclenToPoca();
183 
184  reco::HitPairListPtr::iterator startItr = hit3DList.begin();
185  reco::HitPairListPtr::iterator lastItr = hit3DList.begin();
186 
187  while(++lastItr != hit3DList.end())
188  {
189  const reco::ClusterHit3D* hit3D = *lastItr;
190  double arcLen = hit3D->getArclenToPoca();
191 
192  if (fabs(arcLen-lastArcLen) > m_gapDistance)
193  {
194  startItr = lastItr;
195  hit2DSet.clear();
196  }
197 
198  for(const auto& hit2D : hit3D->getHits())
199  {
200  hit2DSet.insert(hit2D);
201  }
202 
203  if (hit2DSet.size() > m_numSeed2DHits) break;
204 
205  lastArcLen = arcLen;
206  }
207 
208  if (hit2DSet.size() > m_numSeed2DHits)
209  {
210  if (startItr != hit3DList.begin()) hit3DList.erase(hit3DList.begin(), startItr);
211  if (lastItr != hit3DList.end() ) hit3DList.erase(lastItr, hit3DList.end());
212 
213  // On input, the seedPca will contain the original values so we can recover the original axis now
214  TVector3 planeVec0(seedPca.getEigenVectors()[0][0],seedPca.getEigenVectors()[0][1],seedPca.getEigenVectors()[0][2]);
215 
216  m_pcaAlg.PCAAnalysis_3D(hit3DList, seedPca, true);
217 
218  if (seedPca.getSvdOK())
219  {
220  // Still looking to point "down"
221  if (seedPca.getEigenVectors()[0][1] > 0.) seedPca.flipAxis(0);
222 
223  // Check that the seed PCA we have found is consistent with the input PCA
224  TVector3 primarySeedAxis(seedPca.getEigenVectors()[0][0],seedPca.getEigenVectors()[0][1],seedPca.getEigenVectors()[0][2]);
225 
226  double cosAng = primarySeedAxis.Dot(planeVec0);
227 
228  // If the proposed seed axis is not relatively aligned with the input PCA then
229  // we should not be using this method to return seeds. Check that here
230  if (cosAng > m_minAllowedCosAng) foundGoodSeedHits = true;
231  }
232  }
233 
234  return foundGoodSeedHits;
235 }
236 
237 //------------------------------------------------------------------------------
239  double XOrigin,
240  TVector3& Pos,
241  TVector3& Dir,
242  double& ChiDOF) const
243 {
244  // The following is lifted from Bruce Baller to try to get better
245  // initial parameters for a candidate Seed. It is slightly reworked
246  // which is why it is included here instead of used as is.
247  //
248  // Linear fit using X as the independent variable. Hits to be fitted
249  // are passed in the hits vector in a pair form (X, WireID). The
250  // fitted track position at XOrigin is returned in the Pos vector.
251  // The direction cosines are returned in the Dir vector.
252  //
253  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
254  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
255  // a matrix to calculate a track projected to a point at X, and v is
256  // a vector (Yo, Zo, dY/dX, dZ/dX).
257  //
258  // Note: The covariance matrix should also be returned
259  // B. Baller August 2014
260 
261  // assume failure
262  ChiDOF = -1;
263 
264  // Make a set of 2D hits so we consider only unique hits
265  std::set<const reco::ClusterHit2D*> hit2DSet;
266 
267  for(const auto& hit3D : hit3DList)
268  {
269  for(const auto& hit2D : hit3D->getHits()) hit2DSet.insert(hit2D);
270  }
271 
272  if(hit2DSet.size() < 4) return;
273 
274  const unsigned int nvars = 4;
275  unsigned int npts = 3 * hit3DList.size();
276 
277  TMatrixD A(npts, nvars);
278  // vector holding the Wire number
279  TVectorD w(npts);
280  unsigned short ninpl[3] = {0};
281  unsigned short nok = 0;
282  unsigned short iht(0), cstat, tpc, ipl;
283  double x, cw, sw, off;
284 
285  // Loop over the 2D hits in the above
286  for (const auto& hit : hit2DSet)
287  {
288  geo::WireID wireID = hit->getHit().WireID();
289 
290  cstat = wireID.Cryostat;
291  tpc = wireID.TPC;
292  ipl = wireID.Plane;
293 
294  // get the wire plane offset
295  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
296 
297  // get the "cosine-like" component
298  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
299 
300  // the "sine-like" component
301  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
302 
303  x = hit->getXPosition() - XOrigin;
304 
305  A[iht][0] = cw;
306  A[iht][1] = sw;
307  A[iht][2] = cw * x;
308  A[iht][3] = sw * x;
309  w[iht] = wireID.Wire - off;
310 
311  ++ninpl[ipl];
312 
313  // need at least two points in a plane
314  if(ninpl[ipl] == 2) ++nok;
315 
316  iht++;
317  }
318 
319  // need at least 2 planes with at least two points
320  if(nok < 2) return;
321 
322  TDecompSVD svd(A);
323  bool ok;
324  TVectorD tVec = svd.Solve(w, ok);
325 
326  ChiDOF = 0;
327 
328  // not enough points to calculate Chisq
329  if(npts <= 4) return;
330 
331  double ypr, zpr, diff;
332 
333  for (const auto& hit : hit2DSet)
334  {
335  geo::WireID wireID = hit->getHit().WireID();
336 
337  cstat = wireID.Cryostat;
338  tpc = wireID.TPC;
339  ipl = wireID.Plane;
340  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
341  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
342  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
343  x = hit->getXPosition() - XOrigin;
344  ypr = tVec[0] + tVec[2] * x;
345  zpr = tVec[1] + tVec[3] * x;
346  diff = ypr * cw + zpr * sw - (wireID.Wire - off);
347  ChiDOF += diff * diff;
348  }
349 
350  float werr2 = m_geometry->WirePitch() * m_geometry->WirePitch();
351  ChiDOF /= werr2;
352  ChiDOF /= (float)(npts - 4);
353 
354  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
355  Dir[0] = 1 / norm;
356  Dir[1] = tVec[2] / norm;
357  Dir[2] = tVec[3] / norm;
358 
359  Pos[0] = XOrigin;
360  Pos[1] = tVec[0];
361  Pos[2] = tVec[1];
362 
363 } // TrkLineFit()
364 
365  /*
366  //------------------------------------------------------------------------------
367  void PCASeedFinderAlg::LineFit2DHits(const reco::HitPairListPtr& hit3DList,
368  double XOrigin,
369  TVector3& Pos,
370  TVector3& Dir,
371  double& ChiDOF) const
372  {
373  // Linear fit using X as the independent variable. Hits to be fitted
374  // are passed in the hits vector in a pair form (X, WireID). The
375  // fitted track position at XOrigin is returned in the Pos vector.
376  // The direction cosines are returned in the Dir vector.
377  //
378  // SVD fit adapted from $ROOTSYS/tutorials/matrix/solveLinear.C
379  // Fit equation is w = A(X)v, where w is a vector of hit wires, A is
380  // a matrix to calculate a track projected to a point at X, and v is
381  // a vector (Yo, Zo, dY/dX, dZ/dX).
382  //
383  // Note: The covariance matrix should also be returned
384  // B. Baller August 2014
385 
386  // assume failure
387  ChiDOF = -1;
388 
389  if(hit3DList.size() < 4) return;
390 
391  const unsigned int nvars = 4;
392  unsigned int npts = 3 * hit3DList.size();
393 
394  TMatrixD A(npts, nvars);
395  // vector holding the Wire number
396  TVectorD w(npts);
397  unsigned short ninpl[3] = {0};
398  unsigned short nok = 0;
399  unsigned short iht(0), cstat, tpc, ipl;
400  double x, cw, sw, off;
401 
402  for (const auto& hit3D : hit3DList)
403  {
404  // Inner loop over the 2D hits in the above
405  for (const auto& hit : hit3D->getHits())
406  {
407  geo::WireID wireID = hit->getHit().WireID();
408 
409  cstat = wireID.Cryostat;
410  tpc = wireID.TPC;
411  ipl = wireID.Plane;
412 
413  // get the wire plane offset
414  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
415 
416  // get the "cosine-like" component
417  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
418 
419  // the "sine-like" component
420  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
421 
422  x = hit->getXPosition() - XOrigin;
423 
424  A[iht][0] = cw;
425  A[iht][1] = sw;
426  A[iht][2] = cw * x;
427  A[iht][3] = sw * x;
428  w[iht] = wireID.Wire - off;
429 
430  ++ninpl[ipl];
431 
432  // need at least two points in a plane
433  if(ninpl[ipl] == 2) ++nok;
434 
435  iht++;
436  }
437  }
438 
439  // need at least 2 planes with at least two points
440  if(nok < 2) return;
441 
442  TDecompSVD svd(A);
443  bool ok;
444  TVectorD tVec = svd.Solve(w, ok);
445 
446  ChiDOF = 0;
447 
448  // not enough points to calculate Chisq
449  if(npts <= 4) return;
450 
451  double ypr, zpr, diff;
452 
453  for(const auto& hit3D : hit3DList)
454  {
455  for (const auto& hit : hit3D->getHits())
456  {
457  geo::WireID wireID = hit->getHit().WireID();
458 
459  cstat = wireID.Cryostat;
460  tpc = wireID.TPC;
461  ipl = wireID.Plane;
462  off = m_geometry->WireCoordinate(0, 0, ipl, tpc, cstat);
463  cw = m_geometry->WireCoordinate(1, 0, ipl, tpc, cstat) - off;
464  sw = m_geometry->WireCoordinate(0, 1, ipl, tpc, cstat) - off;
465  x = hit->getXPosition() - XOrigin;
466  ypr = tVec[0] + tVec[2] * x;
467  zpr = tVec[1] + tVec[3] * x;
468  diff = ypr * cw + zpr * sw - (wireID.Wire - off);
469  ChiDOF += diff * diff;
470  }
471  }
472 
473  float werr2 = m_geometry->WirePitch() * m_geometry->WirePitch();
474  ChiDOF /= werr2;
475  ChiDOF /= (float)(npts - 4);
476 
477  double norm = sqrt(1 + tVec[2] * tVec[2] + tVec[3] * tVec[3]);
478  Dir[0] = 1 / norm;
479  Dir[1] = tVec[2] / norm;
480  Dir[2] = tVec[3] / norm;
481 
482  Pos[0] = XOrigin;
483  Pos[1] = tVec[0];
484  Pos[2] = tVec[1];
485 
486  } // TrkLineFit()
487  */
488 
489 } // namespace lar_cluster3d
Float_t x
Definition: compare.C:6
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
std::vector< SeedHitPairListPair > SeedHitPairListPairVec
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 getHitsAtEnd(reco::HitPairListPtr &hit3DList, reco::PrincipalComponents &seedPca) const
Separate function to find hits at the ends of the input hits.
bool getSvdOK() const
Definition: Cluster3D.h:224
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
void flipAxis(size_t axis)
Definition: Cluster3D.cxx:195
intermediate_table::iterator iterator
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
Declaration of signal hit object.
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.
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
PrincipalComponentsAlg m_pcaAlg
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
This is an algorithm for finding recob::Seed objects in 3D clusters.
Define a comparator which will sort hits by the absolute value of arc length so hits are ordered clos...
virtual void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
parameter set interface
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:315
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
Detector simulation of raw signals on wires.
Encapsulate the geometry of a wire.
double m_minAllowedCosAng
The minimum cos(ang) between input and seed axes.
Float_t sw
Definition: plot.C:23
virtual ~PCASeedFinderAlg()
Destructor.
PCASeedFinderAlg(fhicl::ParameterSet const &pset)
Constructor.
const float * getEigenValues() const
Definition: Cluster3D.h:226
Float_t norm
Encapsulate the construction of a single detector plane.
float getArclenToPoca() const
Definition: Cluster3D.h:155
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
void LineFit2DHits(const reco::HitPairListPtr &hitList, double XOrigin, TVector3 &Pos, TVector3 &Dir, double &ChiDOF) const
Float_t w
Definition: plot.C:23
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:156
art framework interface to geometry description
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:227