LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
PrincipalComponentsAlg.cxx
Go to the documentation of this file.
1 
9 
10 // Framework Includes
12 #include "fhiclcpp/ParameterSet.h"
14 
15 // LArSoft includes
22 
23 // std includes
24 #include <functional>
25 #include <iostream>
26 #include <numeric>
27 
28 // Eigen includes
29 #include "Eigen/Core"
30 #include "Eigen/Dense"
31 #include "Eigen/Eigenvalues"
32 #include "Eigen/Geometry"
33 #include "Eigen/Jacobi"
34 
35 //------------------------------------------------------------------------------------------------------------------------------------------
36 // implementation follows
37 
38 namespace lar_cluster3d {
39 
41  {
42  m_parallel = pset.get<float>("ParallelLines", 0.00001);
44  }
45 
48  {
49  return left->getDocaToAxis() < right->getDocaToAxis();
50  }
51  };
52 
55  {
56  return left->getArclenToPoca() < right->getArclenToPoca();
57  }
58  };
59 
62  {
63  return fabs(left->getArclenToPoca()) < fabs(right->getArclenToPoca());
64  }
65  };
66 
68  const reco::HitPairListPtr& hitPairVector,
70  float doca3DScl) const
71  {
72  // This is the controlling outside function for running
73  // a Principal Components Analysis on the hits in our
74  // input cluster.
75  // There is a bootstrap process to be followed
76  // 1) Get the initial results from all 3D hits associated
77  // with the cluster
78  // 2) Refine the axis by using only the 2D hits on wires
79 
80  // Get the first axis from all 3D hits
81  PCAAnalysis_3D(hitPairVector, pca);
82 
83  // Make sure first pass was good
84  if (!pca.getSvdOK()) return;
85 
86  // First attempt to refine it using only 2D information
87  reco::PrincipalComponents pcaLoop = pca;
88 
89  PCAAnalysis_2D(detProp, hitPairVector, pcaLoop);
90 
91  // If valid result then go to next steps
92  if (pcaLoop.getSvdOK()) {
93  // Let's check the angle between the original and the updated axis
94  float cosAngle = pcaLoop.getEigenVectors().row(0) * pca.getEigenVectors().row(0).transpose();
95 
96  // Set the scale factor for the outlier rejection
97  float sclFctr(3.);
98 
99  // If we had a significant change then let's do some outlier rejection, etc.
100  if (cosAngle < 1.0) // pretty much everyone takes a turn
101  {
102  int maxIterations(3);
103  float maxRange = 3. * sqrt(pcaLoop.getEigenValues()[1]);
104  float aveDoca = pcaLoop.getAveHitDoca(); // was 0.2
105 
106  maxRange = sclFctr * 0.5 * (maxRange + aveDoca); // was std::max(maxRange, aveDoca);
107 
108  int numRejHits = PCAAnalysis_reject2DOutliers(hitPairVector, maxRange);
109  int totalRejects(numRejHits);
110  int maxRejects(0.4 * hitPairVector.size());
111 
112  // Try looping to see if we improve things
113  while (maxIterations-- && numRejHits > 0 && totalRejects < maxRejects) {
114  // Run the PCA
115  PCAAnalysis_2D(detProp, hitPairVector, pcaLoop, true);
116 
117  maxRange =
118  sclFctr * 0.5 * (3. * sqrt(pcaLoop.getEigenValues()[1]) + pcaLoop.getAveHitDoca());
119 
120  numRejHits = PCAAnalysis_reject2DOutliers(hitPairVector, maxRange);
121  }
122  }
123 
124  // Ok at this stage copy the latest results back into the cluster
125  pca = pcaLoop;
126 
127  // Now we make one last pass through the 3D hits to reject outliers there
128  PCAAnalysis_reject3DOutliers(hitPairVector, pca, doca3DScl * pca.getAveHitDoca());
129  }
130  else
131  pca = pcaLoop;
132 
133  return;
134  }
135 
138  bool skeletonOnly) const
139  {
140  // We want to run a PCA on the input TkrVecPoints...
141  // The steps are:
142  // 1) do a mean normalization of the input vec points
143  // 2) compute the covariance matrix
144  // 3) run the SVD
145  // 4) extract the eigen vectors and values
146  // see what happens
147 
148  // Run through the HitPairList and get the mean position of all the hits
149  Eigen::Vector3d meanPos(Eigen::Vector3d::Zero());
150  double meanWeightSum(0.);
151  int numPairsInt(0);
152 
153  // const float minimumDeltaPeakSig(0.00001);
154  double minimumDeltaPeakSig(0.00001);
155 
156  // Want to use the hit "chi square" to weight the hits but we need to put a lower limit on its value
157  // to prevent a few hits being over counted.
158  // This is a bit experimental until we can evaluate the cost (time to calculate) vs the benefit
159  // (better fits)..
160  std::vector<double> hitChiSquareVec;
161 
162  hitChiSquareVec.resize(hitPairVector.size());
163 
164  std::transform(hitPairVector.begin(),
165  hitPairVector.end(),
166  hitChiSquareVec.begin(),
167  [](const auto& hit) { return hit->getHitChiSquare(); });
168  std::sort(hitChiSquareVec.begin(), hitChiSquareVec.end());
169 
170  size_t numToKeep = 0.8 * hitChiSquareVec.size();
171 
172  hitChiSquareVec.resize(numToKeep);
173 
174  double aveValue = std::accumulate(hitChiSquareVec.begin(), hitChiSquareVec.end(), double(0.)) /
175  double(hitChiSquareVec.size());
176  double rms = std::sqrt(std::inner_product(hitChiSquareVec.begin(),
177  hitChiSquareVec.end(),
178  hitChiSquareVec.begin(),
179  0.,
180  std::plus<>(),
181  [aveValue](const auto& left, const auto& right) {
182  return (left - aveValue) * (right - aveValue);
183  }) /
184  double(hitChiSquareVec.size()));
185 
186  minimumDeltaPeakSig = std::max(minimumDeltaPeakSig, aveValue - rms);
187 
188  // std::cout << "===>> Calculating PCA, ave chiSquare: " << aveValue << ", rms: " << rms << ", cut: " << minimumDeltaPeakSig << std::endl;
189 
190  for (const auto& hit : hitPairVector) {
191  if (skeletonOnly && !((hit->getStatusBits() & reco::ClusterHit3D::SKELETONHIT) ==
193  continue;
194 
195  // Weight the hit by the peak time difference significance
196  double weight = std::max(minimumDeltaPeakSig,
197  double(hit->getHitChiSquare())); // hit->getDeltaPeakTime());
198  // ///hit->getSigmaPeakTime());
199 
200  meanPos(0) += hit->getPosition()[0] * weight;
201  meanPos(1) += hit->getPosition()[1] * weight;
202  meanPos(2) += hit->getPosition()[2] * weight;
203  numPairsInt++;
204 
205  meanWeightSum += weight;
206  }
207 
208  meanPos /= meanWeightSum;
209 
210  // Define elements of our covariance matrix
211  double xi2(0.);
212  double xiyi(0.);
213  double xizi(0.0);
214  double yi2(0.0);
215  double yizi(0.0);
216  double zi2(0.);
217  double weightSum(0.);
218 
219  // Back through the hits to build the matrix
220  for (const auto& hit : hitPairVector) {
221  if (skeletonOnly && !((hit->getStatusBits() & reco::ClusterHit3D::SKELETONHIT) ==
223  continue;
224 
225  double weight = 1. / std::max(minimumDeltaPeakSig,
226  double(hit->getHitChiSquare())); // hit->getDeltaPeakTime());
227  // ///hit->getSigmaPeakTime());
228  double x = (hit->getPosition()[0] - meanPos(0)) * weight;
229  double y = (hit->getPosition()[1] - meanPos(1)) * weight;
230  double z = (hit->getPosition()[2] - meanPos(2)) * weight;
231 
232  weightSum += weight * weight;
233 
234  xi2 += x * x;
235  xiyi += x * y;
236  xizi += x * z;
237  yi2 += y * y;
238  yizi += y * z;
239  zi2 += z * z;
240  }
241 
242  // Using Eigen package
243  Eigen::Matrix3d sig;
244 
245  sig << xi2, xiyi, xizi, xiyi, yi2, yizi, xizi, yizi, zi2;
246 
247  sig *= 1. / weightSum;
248 
249  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigenMat(sig);
250 
251  if (eigenMat.info() == Eigen::ComputationInfo::Success) {
252  // Now copy output
253  // The returned eigen values and vectors will be returned in an xyz system where x is the smallest spread,
254  // y is the next smallest and z is the largest. Adopt that convention going forward
255  reco::PrincipalComponents::EigenValues recobEigenVals = eigenMat.eigenvalues().cast<float>();
257  eigenMat.eigenvectors().transpose().cast<float>();
258 
259  // Check for a special case (which may have gone away with switch back to doubles for computation?)
260  if (std::isnan(recobEigenVals[0])) {
261  std::cout << "==> Third eigenvalue returns a nan" << std::endl;
262 
263  recobEigenVals[0] = 0.;
264 
265  // Assume the third axis is also kaput?
266  recobEigenVecs.row(0) = recobEigenVecs.row(1).cross(recobEigenVecs.row(2));
267  }
268 
269  // Store away
271  true, numPairsInt, recobEigenVals, recobEigenVecs, meanPos.cast<float>());
272  }
273  else {
274  mf::LogDebug("Cluster3D") << "PCA decompose failure, numPairs = " << numPairsInt << std::endl;
276  }
277 
278  return;
279  }
280 
282  const reco::HitPairListPtr& hitPairVector,
284  bool updateAvePos) const
285  {
286  // Once an axis has been found our goal is to refine it by using only the 2D hits
287  // We'll get 3D information for each of these by using the axis as a reference and use
288  // the point of closest approach as the 3D position
289 
290  // Define elements of our covariance matrix
291  float xi2(0.);
292  float xiyi(0.);
293  float xizi(0.);
294  float yi2(0.);
295  float yizi(0.);
296  float zi2(0.);
297 
298  float aveHitDoca(0.);
299  Eigen::Vector3f avePosUpdate(0., 0., 0.);
300  int nHits(0);
301 
302  // Recover existing line parameters for current cluster
303  const reco::PrincipalComponents& inputPca = pca;
304  Eigen::Vector3f avePosition(inputPca.getAvePosition());
305  Eigen::Vector3f axisDirVec(inputPca.getEigenVectors().row(0));
306 
307  // We float loop here so we can use this method for both the first time through
308  // and a second time through where we re-calculate the mean position
309  // So, we need to keep track of the poca which we do with a float vector
310  std::vector<Eigen::Vector3f> hitPosVec;
311 
312  // Outer loop over 3D hits
313  for (const auto& hit3D : hitPairVector) {
314  // Inner loop over 2D hits
315  for (const auto& hit : hit3D->getHits()) {
316  // Step one is to set up to determine the point of closest approach of this 2D hit to
317  // the cluster's current axis.
318  // Get this wire's geometry object
319  const geo::WireID& hitID = hit->WireID();
320  const geo::WireGeo& wire_geom = m_geometry->WireIDToWireGeo(hitID);
321 
322  // From this, get the parameters of the line for the wire
323  auto const wirePosArr = wire_geom.GetCenter();
324 
325  Eigen::Vector3f wireCenter(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
326  Eigen::Vector3f wireDirVec(
327  wire_geom.Direction().X(), wire_geom.Direction().Y(), wire_geom.Direction().Z());
328 
329  // Correct the wire position in x to set to correspond to the drift time
330  float hitPeak(hit->getHit()->PeakTime());
331 
332  Eigen::Vector3f wirePos(
333  detProp.ConvertTicksToX(hitPeak, hitID.Plane, hitID.TPC, hitID.Cryostat),
334  wireCenter[1],
335  wireCenter[2]);
336 
337  // Compute the wire plane normal for this view
338  Eigen::Vector3f xAxis(1., 0., 0.);
339  Eigen::Vector3f planeNormal =
340  xAxis.cross(wireDirVec); // This gives a normal vector in +z for a Y wire
341 
342  float docaInPlane(wirePos[0] - avePosition[0]);
343  float arcLenToPlane(0.);
344  float cosAxisToPlaneNormal = axisDirVec.dot(planeNormal);
345 
346  Eigen::Vector3f axisPlaneIntersection = wirePos;
347  Eigen::Vector3f hitPosTVec = wirePos;
348 
349  if (fabs(cosAxisToPlaneNormal) > 0.) {
350  Eigen::Vector3f deltaPos = wirePos - avePosition;
351 
352  arcLenToPlane = deltaPos.dot(planeNormal) / cosAxisToPlaneNormal;
353  axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
354  docaInPlane = wirePos[0] - axisPlaneIntersection[0];
355 
356  Eigen::Vector3f axisToInter = axisPlaneIntersection - wirePos;
357  float arcLenToDoca = axisToInter.dot(wireDirVec);
358 
359  hitPosTVec += arcLenToDoca * wireDirVec;
360  }
361 
362  // Get a vector from the wire position to our cluster's current average position
363  Eigen::Vector3f wVec = avePosition - wirePos;
364 
365  // Get the products we need to compute the arc lengths to the distance of closest approach
366  float a(axisDirVec.dot(axisDirVec));
367  float b(axisDirVec.dot(wireDirVec));
368  float c(wireDirVec.dot(wireDirVec));
369  float d(axisDirVec.dot(wVec));
370  float e(wireDirVec.dot(wVec));
371 
372  float den(a * c - b * b);
373  float arcLen1(0.);
374  float arcLen2(0.);
375 
376  // Parallel lines is a special case
377  if (fabs(den) > m_parallel) {
378  arcLen1 = (b * e - c * d) / den;
379  arcLen2 = (a * e - b * d) / den;
380  }
381  else {
382  mf::LogDebug("Cluster3D") << "** Parallel lines, need a solution here" << std::endl;
383  break;
384  }
385 
386  // Now get the hit position we'll use for the pca analysis
387  //float hitPos[] = {wirePos[0] + arcLen2 * wireDirVec[0],
388  // wirePos[1] + arcLen2 * wireDirVec[1],
389  // wirePos[2] + arcLen2 * wireDirVec[2]};
390  //float axisPos[] = {avePosition[0] + arcLen1 * axisDirVec[0],
391  // avePosition[1] + arcLen1 * axisDirVec[1],
392  // avePosition[2] + arcLen1 * axisDirVec[2]};
393  Eigen::Vector3f hitPos = wirePos + arcLen2 * wireDirVec;
394  Eigen::Vector3f axisPos = avePosition + arcLen1 * axisDirVec;
395  float deltaX = hitPos(0) - axisPos(0);
396  float deltaY = hitPos(1) - axisPos(1);
397  float deltaZ = hitPos(2) - axisPos(2);
398  float doca2 = deltaX * deltaX + deltaY * deltaY + deltaZ * deltaZ;
399  float doca = sqrt(doca2);
400 
401  docaInPlane = doca;
402 
403  aveHitDoca += fabs(docaInPlane);
404 
405  //Eigen::Vector3f deltaPos = hitPos - hitPosTVec;
406  //float deltaDoca = doca - docaInPlane;
407 
408  //if (fabs(deltaPos[0]) > 1. || fabs(deltaPos[1]) > 1. || fabs(deltaPos[2]) > 1. || fabs(deltaDoca) > 2.)
409  //{
410  // std::cout << "**************************************************************************************" << std::endl;
411  // std::cout << "Diff in doca: " << deltaPos[0] << "," << deltaPos[1] << "," << deltaPos[2] << ", deltaDoca: " << deltaDoca << std::endl;
412  // std::cout << "-- HitPosTVec: " << hitPosTVec[0] << "," << hitPosTVec[1] << "," << hitPosTVec[2] << std::endl;
413  // std::cout << "-- hitPos: " << hitPos[0] << "," << hitPos[1] << "," << hitPos[2] << std::endl;
414  // std::cout << "-- WirePos: " << wirePos[0] << "," << wirePos[1] << "," << wirePos[2] << std::endl;
415  //}
416 
417  // Set the hit's doca and arclen
418  hit->setDocaToAxis(fabs(docaInPlane));
419  hit->setArcLenToPoca(arcLenToPlane);
420 
421  // If this point is considered an outlier then we skip
422  // the accumulation for average position and covariance
423  if (hit->getStatusBits() & 0x80) { continue; }
424 
425  //hitPosTVec = hitPos;
426 
427  avePosUpdate += hitPosTVec;
428 
429  hitPosVec.push_back(hitPosTVec);
430 
431  nHits++;
432  }
433  }
434 
435  // Get updated average position
436  avePosUpdate /= float(nHits);
437 
438  // Get the average hit doca
439  aveHitDoca /= float(nHits);
440 
441  if (updateAvePos) { avePosition = avePosUpdate; }
442 
443  // Now loop through the hits and build out the covariance matrix
444  for (auto& hitPos : hitPosVec) {
445  // And increment the values in the covariance matrix
446  float x = hitPos[0] - avePosition[0];
447  float y = hitPos[1] - avePosition[1];
448  float z = hitPos[2] - avePosition[2];
449 
450  xi2 += x * x;
451  xiyi += x * y;
452  xizi += x * z;
453  yi2 += y * y;
454  yizi += y * z;
455  zi2 += z * z;
456  }
457 
458  // Accumulation done, now do the actual work
459 
460  // Using Eigen package
461  Eigen::Matrix3f sig;
462 
463  sig << xi2, xiyi, xizi, xiyi, yi2, yizi, xizi, yizi, zi2;
464 
465  sig *= 1. / (nHits - 1);
466 
467  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> eigenMat(sig);
468 
469  if (eigenMat.info() == Eigen::ComputationInfo::Success) {
470  // Now copy output
471  // The returned eigen values and vectors will be returned in an xyz system where x is the smallest spread,
472  // y is the next smallest and z is the largest. Adopt that convention going forward
473  reco::PrincipalComponents::EigenValues recobEigenVals = eigenMat.eigenvalues().cast<float>();
475  eigenMat.eigenvectors().transpose().cast<float>();
476 
477  // Store away
479  true, nHits, recobEigenVals, recobEigenVecs, avePosition, aveHitDoca);
480  }
481  else {
482  mf::LogDebug("Cluster3D") << "PCA decompose failure, numPairs = " << nHits << std::endl;
484  }
485 
486  return;
487  }
488 
490  const reco::PrincipalComponents& pca) const
491  {
492  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
493  // any outliers. Basically, any hit outside a scaled range of the average doca from the
494  // first pass is marked by setting the bit in the status word.
495 
496  // We'll need the current PCA axis to determine doca and arclen
497  Eigen::Vector3f avePosition(
498  pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
499  Eigen::Vector3f axisDirVec(pca.getEigenVectors().row(2));
500 
501  // We want to keep track of the average
502  float aveDoca3D(0.);
503 
504  // Outer loop over views
505  for (const auto* clusterHit3D : hitPairVector) {
506  // Always reset the existing status bit
507  clusterHit3D->clearStatusBits(0x80);
508 
509  // Now we need to calculate the doca and poca...
510  // Start by getting this hits position
511  Eigen::Vector3f clusPos(clusterHit3D->getPosition()[0],
512  clusterHit3D->getPosition()[1],
513  clusterHit3D->getPosition()[2]);
514 
515  // Form a TVector from this to the cluster average position
516  Eigen::Vector3f clusToHitVec = clusPos - avePosition;
517 
518  // With this we can get the arclength to the doca point
519  float arclenToPoca = clusToHitVec.dot(axisDirVec);
520 
521  // Get the coordinates along the axis for this point
522  Eigen::Vector3f docaPos = avePosition + arclenToPoca * axisDirVec;
523 
524  // Now get doca and poca
525  Eigen::Vector3f docaPosToClusPos = clusPos - docaPos;
526  float docaToAxis = docaPosToClusPos.norm();
527 
528  aveDoca3D += docaToAxis;
529 
530  // Ok, set the values in the hit
531  clusterHit3D->setDocaToAxis(docaToAxis);
532  clusterHit3D->setArclenToPoca(arclenToPoca);
533  }
534 
535  // Compute the average and store
536  aveDoca3D /= float(hitPairVector.size());
537 
538  pca.setAveHitDoca(aveDoca3D);
539 
540  return;
541  }
542 
544  const reco::PrincipalComponents& pca) const
545  {
546  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
547  // any outliers. Basically, any hit outside a scaled range of the average doca from the
548  // first pass is marked by setting the bit in the status word.
549 
550  // We'll need the current PCA axis to determine doca and arclen
551  Eigen::Vector3f avePosition(
552  pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
553  Eigen::Vector3f axisDirVec(pca.getEigenVectors().row(2));
554 
555  // Recover the principle eigen value for range constraints
556  float maxArcLen = 4. * sqrt(pca.getEigenValues()[0]);
557 
558  // We want to keep track of the average
559  float aveHitDoca(0.);
560 
561  // Outer loop over views
562  for (const auto* hit : hit2DListPtr) {
563  // Step one is to set up to determine the point of closest approach of this 2D hit to
564  // the cluster's current axis. We do that by finding the point of intersection of the
565  // cluster's axis with a plane defined by the wire the hit is associated with.
566  // Get this wire's geometry object
567  const geo::WireID& hitID = hit->WireID();
568  const geo::WireGeo& wire_geom = m_geometry->WireIDToWireGeo(hitID);
569 
570  // From this, get the parameters of the line for the wire
571  auto const wirePosArr = wire_geom.GetCenter();
572 
573  Eigen::Vector3f wireCenter(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
574  Eigen::Vector3f wireDirVec(
575  wire_geom.Direction().X(), wire_geom.Direction().Y(), wire_geom.Direction().Z());
576 
577  // Correct the wire position in x to set to correspond to the drift time
578  Eigen::Vector3f wirePos(hit->getXPosition(), wireCenter[1], wireCenter[2]);
579 
580  // Compute the wire plane normal for this view
581  Eigen::Vector3f xAxis(1., 0., 0.);
582  Eigen::Vector3f planeNormal =
583  xAxis.cross(wireDirVec); // This gives a normal vector in +z for a Y wire
584 
585  float arcLenToPlane(0.);
586  float docaInPlane(wirePos[0] - avePosition[0]);
587  float cosAxisToPlaneNormal = axisDirVec.dot(planeNormal);
588 
589  Eigen::Vector3f axisPlaneIntersection = wirePos;
590 
591  // If current cluster axis is not parallel to wire plane then find intersection point
592  if (fabs(cosAxisToPlaneNormal) > 0.) {
593  Eigen::Vector3f deltaPos = wirePos - avePosition;
594 
595  arcLenToPlane =
596  std::min(float(deltaPos.dot(planeNormal) / cosAxisToPlaneNormal), maxArcLen);
597  axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
598 
599  Eigen::Vector3f axisToInter = axisPlaneIntersection - wirePos;
600  float arcLenToDoca = axisToInter.dot(wireDirVec);
601 
602  // If the arc length along the wire to the poca is outside the TPC then reset
603  if (fabs(arcLenToDoca) > wire_geom.HalfL()) arcLenToDoca = wire_geom.HalfL();
604 
605  // If we were successful in getting to the wire plane then the doca is simply the
606  // difference in x coordinates... but we hvae to worry about the special cases so
607  // we calculate a 3D doca based on arclengths above...
608  Eigen::Vector3f docaVec = axisPlaneIntersection - (wirePos + arcLenToDoca * wireDirVec);
609  docaInPlane = docaVec.norm();
610  }
611 
612  aveHitDoca += fabs(docaInPlane);
613 
614  // Set the hit's doca and arclen
615  hit->setDocaToAxis(fabs(docaInPlane));
616  hit->setArcLenToPoca(arcLenToPlane);
617  }
618 
619  // Compute the average and store
620  aveHitDoca /= float(hit2DListPtr.size());
621 
622  pca.setAveHitDoca(aveHitDoca);
623 
624  return;
625  }
626 
628  const reco::HitPairListPtr& hitPairVector,
629  float maxDocaAllowed) const
630  {
631  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
632  // any outliers. Basically, any hit outside a scaled range of the average doca from the
633  // first pass is marked by setting the bit in the status word.
634 
635  // First get the average doca scaled by some appropriate factor
636  int numRejHits(0);
637 
638  // Outer loop over views
639  for (const auto& hit3D : hitPairVector) {
640  // Inner loop over hits in this view
641  for (const auto& hit : hit3D->getHits()) {
642  // Always reset the existing status bit
643  hit->clearStatusBits(0x80);
644 
645  if (hit->getDocaToAxis() > maxDocaAllowed) {
646  hit->setStatusBit(0x80);
647  numRejHits++;
648  }
649  }
650  }
651 
652  return numRejHits;
653  }
654 
656  const reco::HitPairListPtr& hitPairVector,
657  const reco::PrincipalComponents& pca,
658  float maxDocaAllowed) const
659  {
660  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
661  // any outliers. Basically, any hit outside a scaled range of the average doca from the
662  // first pass is marked by setting the bit in the status word.
663 
664  // First get the average doca scaled by some appropriate factor
665  int numRejHits(0);
666 
667  // We'll need the current PCA axis to determine doca and arclen
668  Eigen::Vector3f avePosition(
669  pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
670  Eigen::Vector3f axisDirVec(pca.getEigenVectors().row(2));
671 
672  // Outer loop over views
673  for (const auto* clusterHit3D : hitPairVector) {
674  // Always reset the existing status bit
675  clusterHit3D->clearStatusBits(0x80);
676 
677  // Now we need to calculate the doca and poca...
678  // Start by getting this hits position
679  Eigen::Vector3f clusPos(clusterHit3D->getPosition()[0],
680  clusterHit3D->getPosition()[1],
681  clusterHit3D->getPosition()[2]);
682 
683  // Form a TVector from this to the cluster average position
684  Eigen::Vector3f clusToHitVec = clusPos - avePosition;
685 
686  // With this we can get the arclength to the doca point
687  float arclenToPoca = clusToHitVec.dot(axisDirVec);
688 
689  // Get the coordinates along the axis for this point
690  Eigen::Vector3f docaPos = avePosition + arclenToPoca * axisDirVec;
691 
692  // Now get doca and poca
693  Eigen::Vector3f docaPosToClusPos = clusPos - docaPos;
694  float docaToAxis = docaPosToClusPos.norm();
695 
696  // Ok, set the values in the hit
697  clusterHit3D->setDocaToAxis(docaToAxis);
698  clusterHit3D->setArclenToPoca(arclenToPoca);
699 
700  // Check to see if this is a keeper
701  if (clusterHit3D->getDocaToAxis() > maxDocaAllowed) {
702  clusterHit3D->setStatusBit(0x80);
703  numRejHits++;
704  }
705  }
706 
707  return numRejHits;
708  }
709 
710 } // namespace lar_cluster3d
Float_t x
Definition: compare.C:6
float m_parallel
means lines are parallel
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:114
void setAveHitDoca(double doca) const
Definition: Cluster3D.h:248
std::list< const reco::ClusterHit2D * > Hit2DListPtr
export some data structure definitions
Definition: Cluster3D.h:325
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
int PCAAnalysis_reject3DOutliers(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca, float aveHitDoca) const
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
bool getSvdOK() const
Definition: Cluster3D.h:240
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.
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
PrincipalComponentsAlg(fhicl::ParameterSet const &pset)
Constructor.
T * get() const
Definition: ServiceHandle.h:69
Float_t den
Definition: plot.C:35
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
int PCAAnalysis_reject2DOutliers(const reco::HitPairListPtr &hitPairVector, float aveHitDoca) const
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
float getDocaToAxis() const
Definition: Cluster3D.h:166
float getAveHitDoca() const
Definition: Cluster3D.h:245
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:242
void PCAAnalysis_2D(const detinfo::DetectorPropertiesData &detProp, const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, bool updateAvePos=false) const
Vector_t Direction() const
Definition: WireGeo.h:289
Eigen::Vector3f EigenValues
Definition: Cluster3D.h:219
T get(std::string const &key) const
Definition: ParameterSet.h:314
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
Float_t d
Definition: plot.C:235
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
Definition of data types for geometry description.
Detector simulation of raw signals on wires.
This header file defines the interface to a principal components analysis designed to be used within ...
Encapsulate the geometry of a wire.
double ConvertTicksToX(double ticks, int p, int t, int c) const
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:172
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
double weight
Definition: plottest35.C:25
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
Eigen::Matrix3f EigenVectors
Definition: Cluster3D.h:220
float getArclenToPoca() const
Definition: Cluster3D.h:167
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
Hit is a "skeleton" hit.
Definition: Cluster3D.h:97
Float_t e
Definition: plot.C:35
void PCAAnalysis_calc2DDocas(const reco::Hit2DListPtr &hit2DVector, const reco::PrincipalComponents &pca) const
art framework interface to geometry description
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
void PCAAnalysis(const detinfo::DetectorPropertiesData &detProp, const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, float doca3DScl=3.) const
Run the Principal Components Analysis.