PrincipalComponentsAlg.cxx
Go to the documentation of this file.
1 /**
2  * @file Cluster3D_module.cc
3  *
4  * @brief Producer module to create 3D clusters from input hits
5  *
6  */
7 
8 // Framework Includes
9 
12 
13 // LArSoft includes
19 #include "lardata/RecoObjects/Cluster3D.h"
23 
24 // ROOT includes
25 #include "TVector3.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 // implementation follows
38 
39 namespace lar_cluster3d {
40 
42 {
43  this->reconfigure(pset);
44 }
45 
46 //------------------------------------------------------------------------------------------------------------------------------------------
47 
49 {
50 }
51 
53 {
55 
56  m_parallel = pset.get<double>("ParallelLines", 0.00001);
57  m_geometry = &*geometry;
58 }
59 
60 void PrincipalComponentsAlg::getHit2DPocaToAxis(const TVector3& axisPos,
61  const TVector3& axisDir,
62  const reco::ClusterHit2D* hit2D,
63  TVector3& poca,
64  double& arcLenAxis,
65  double& arcLenWire,
66  double& doca)
67 {
68  // Step one is to set up to determine the point of closest approach of this 2D hit to
69  // the cluster's current axis.
70  // Get this wire's geometry object
71  const geo::WireID& hitID = hit2D->getHit().WireID();
72  const geo::WireGeo& wire_geom = m_geometry->WireIDToWireGeo(hitID);
73 
74  // From this, get the parameters of the line for the wire
75  double wirePos[3] = {0.,0.,0.};
76  TVector3 wireDir(wire_geom.Direction());
77 
78  wire_geom.GetCenter(wirePos);
79 
80  // Correct the wire position in x to set to correspond to the drift time
81  double hitPeak(hit2D->getHit().PeakTime());
82 
83  wirePos[0] = m_detector->ConvertTicksToX(hitPeak, hitID.Plane, hitID.TPC, hitID.Cryostat);
84 
85  // Get a vector from the wire position to our cluster's current average position
86  TVector3 wVec(axisPos.X()-wirePos[0], axisPos.Y()-wirePos[1], axisPos.Z()-wirePos[2]);
87 
88  // Get the products we need to compute the arc lengths to the distance of closest approach
89  double a(axisDir.Dot(axisDir));
90  double b(axisDir.Dot(wireDir));
91  double c(wireDir.Dot(wireDir));
92  double d(axisDir.Dot(wVec));
93  double e(wireDir.Dot(wVec));
94 
95  double den(a*c - b*b);
96 
97  // Parallel lines is a special case
98  if (fabs(den) > m_parallel)
99  {
100  arcLenAxis = (b*e - c*d) / den;
101  arcLenWire = (a*e - b*d) / den;
102  }
103  else
104  {
105  mf::LogDebug("Cluster3D") << "** Parallel lines, need a solution here" << std::endl;
106  arcLenAxis = 0.;
107  arcLenWire = 0.;
108  }
109 
110  // Now get the hit position we'll use for the pca analysis
111  poca = TVector3(wirePos[0] + arcLenWire * wireDir[0],
112  wirePos[1] + arcLenWire * wireDir[1],
113  wirePos[2] + arcLenWire * wireDir[2]);
114  TVector3 axisPocaPos(axisPos[0] + arcLenAxis * axisDir[0],
115  axisPos[1] + arcLenAxis * axisDir[1],
116  axisPos[2] + arcLenAxis * axisDir[2]);
117 
118  double deltaX(poca.X() - axisPocaPos.X());
119  double deltaY(poca.Y() - axisPocaPos.Y());
120  double deltaZ(poca.Z() - axisPocaPos.Z());
121  double doca2(deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ);
122 
123  doca = sqrt(doca2);
124 
125  return;
126 }
127 
128 
130 {
132  {
133  return left->getDocaToAxis() < right->getDocaToAxis();
134  }
135 
136 };
137 
139 {
141  {
142  return left->getArclenToPoca() < right->getArclenToPoca();
143  }
144 
145 };
146 
148 {
150  {
151  return fabs(left->getArclenToPoca()) < fabs(right->getArclenToPoca());
152  }
153 
154 };
155 
156 void PrincipalComponentsAlg::PCAAnalysis(const reco::HitPairListPtr& hitPairVector, reco::PrincipalComponents& pca, double doca3DScl) const
157 {
158  // This is the controlling outside function for running
159  // a Principal Components Analysis on the hits in our
160  // input cluster.
161  // There is a bootstrap process to be followed
162  // 1) Get the initial results from all 3D hits associated
163  // with the cluster
164  // 2) Refine the axis by using only the 2D hits on wires
165 
166  // Get the first axis from all 3D hits
167  PCAAnalysis_3D(hitPairVector, pca);
168 
169  // Make sure first pass was good
170  if (!pca.getSvdOK()) return;
171 
172  // First attempt to refine it using only 2D information
173  reco::PrincipalComponents pcaLoop = pca;
174 
175  PCAAnalysis_2D(hitPairVector, pcaLoop);
176 
177  // If valid result then go to next steps
178  if (pcaLoop.getSvdOK())
179  {
180  // Let's check the angle between the original and the updated axis
181  double cosAngle = pcaLoop.getEigenVectors()[0][0] * pca.getEigenVectors()[0][0]
182  + pcaLoop.getEigenVectors()[0][1] * pca.getEigenVectors()[0][1]
183  + pcaLoop.getEigenVectors()[0][2] * pca.getEigenVectors()[0][2];
184 
185  // Set the scale factor for the outlier rejection
186  double sclFctr(3.);
187 
188  // If we had a significant change then let's do some outlier rejection, etc.
189  if (cosAngle < 1.0) // pretty much everyone takes a turn
190  {
191  int maxIterations(3);
192  double maxRange = 3.*sqrt(pcaLoop.getEigenValues()[1]);
193  double aveDoca = pcaLoop.getAveHitDoca(); // was 0.2
194 
195  maxRange = sclFctr * 0.5*(maxRange+aveDoca); // was std::max(maxRange, aveDoca);
196 
197  int numRejHits = PCAAnalysis_reject2DOutliers(hitPairVector, pcaLoop, maxRange);
198  int totalRejects(numRejHits);
199  int maxRejects(0.4*hitPairVector.size());
200 
201  // Try looping to see if we improve things
202  while(maxIterations-- && numRejHits > 0 && totalRejects < maxRejects)
203  {
204  // Run the PCA
205  PCAAnalysis_2D(hitPairVector, pcaLoop, true);
206 
207  maxRange = sclFctr * 0.5*(3.*sqrt(pcaLoop.getEigenValues()[1])+pcaLoop.getAveHitDoca());
208 
209  numRejHits = PCAAnalysis_reject2DOutliers(hitPairVector, pcaLoop, maxRange);
210  }
211 
212  }
213 
214  // Ok at this stage copy the latest results back into the cluster
215  pca = pcaLoop;
216 
217  // Now we make one last pass through the 3D hits to reject outliers there
218  PCAAnalysis_reject3DOutliers(hitPairVector, pca, doca3DScl * pca.getAveHitDoca());
219  }
220  else pca = pcaLoop;
221 
222  return;
223 }
224 
225 void PrincipalComponentsAlg::PCAAnalysis_3D(const reco::HitPairListPtr& hitPairVector, reco::PrincipalComponents& pca, bool skeletonOnly) const
226 {
227  // We want to run a PCA on the input TkrVecPoints...
228  // The steps are:
229  // 1) do a mean normalization of the input vec points
230  // 2) compute the covariance matrix
231  // 3) run the SVD
232  // 4) extract the eigen vectors and values
233  // see what happens
234 
235  // Run through the HitPairList and get the mean position of all the hits
236  double meanPos[] = {0.,0.,0.};
237  int numPairsInt(0);
238 
239  for (const auto& hit : hitPairVector)
240  {
241  if (skeletonOnly && !((hit->getStatusBits() & 0x10000000) == 0x10000000)) continue;
242 
243  meanPos[0] += hit->getPosition()[0];
244  meanPos[1] += hit->getPosition()[1];
245  meanPos[2] += hit->getPosition()[2];
246  numPairsInt++;
247  }
248 
249  double numPairs = double(numPairsInt);
250 
251  meanPos[0] /= numPairs;
252  meanPos[1] /= numPairs;
253  meanPos[2] /= numPairs;
254 
255  // Define elements of our covariance matrix
256  double xi2 = 0.;
257  double xiyi = 0.;
258  double xizi = 0.;
259  double yi2 = 0.;
260  double yizi = 0.;
261  double zi2 = 0.;
262 
263  // Back through the hits to build the matrix
264  for (const auto& hit : hitPairVector)
265  {
266  if (skeletonOnly && !((hit->getStatusBits() & 0x10000000) == 0x10000000)) continue;
267 
268  double x = hit->getPosition()[0] - meanPos[0];
269  double y = hit->getPosition()[1] - meanPos[1];
270  double z = hit->getPosition()[2] - meanPos[2];
271 
272  xi2 += x * x;
273  xiyi += x * y;
274  xizi += x * z;
275  yi2 += y * y;
276  yizi += y * z;
277  zi2 += z * z;
278  }
279 
280  // Create the actual matrix
281  TMatrixD sigma(3, 3);
282 
283  sigma(0,0) = xi2;
284  sigma(0,1) = sigma(1,0) = xiyi;
285  sigma(0,2) = sigma(2,0) = xizi;
286  sigma(1,1) = yi2;
287  sigma(1,2) = sigma(2,1) = yizi;
288  sigma(2,2) = zi2;
289 
290  // Scale by number of pairs
291  sigma *= (1./(numPairs - 1.));
292 
293  // Set up the SVD
294  TDecompSVD rootSVD(sigma);
295 
296  // run the decomposition
297  bool svdOk(false);
298 
299  try
300  {
301  svdOk = rootSVD.Decompose();
302  }
303  catch(...)
304  {
305  svdOk = false;
306  mf::LogDebug("Cluster3D") << "PCA decompose failure, numPairs = " << numPairs << std::endl;
307  }
308 
309  if (svdOk)
310  {
311  // Extract results
312  TVectorD eigenVals = rootSVD.GetSig();
313  TMatrixD eigenVecs = rootSVD.GetU();
314 
315  // Get the eigen values
316  double recobEigenVals[] = {eigenVals[0], eigenVals[1], eigenVals[2]};
317 
318  // Grab the principle axes
320  std::vector<double> tempVec;
321 
322  // Get the first column vector
323  tempVec.push_back(eigenVecs(0, 0));
324  tempVec.push_back(eigenVecs(1, 0));
325  tempVec.push_back(eigenVecs(2, 0));
326  recobEigenVecs.push_back(tempVec);
327 
328  // Now the second
329  tempVec.clear();
330  tempVec.push_back(eigenVecs(0, 1));
331  tempVec.push_back(eigenVecs(1, 1));
332  tempVec.push_back(eigenVecs(2, 1));
333  recobEigenVecs.push_back(tempVec);
334 
335  // And the last
336  tempVec.clear();
337  tempVec.push_back(eigenVecs(0, 2));
338  tempVec.push_back(eigenVecs(1, 2));
339  tempVec.push_back(eigenVecs(2, 2));
340  recobEigenVecs.push_back(tempVec);
341 
342  // Store away
343  pca = reco::PrincipalComponents(svdOk, numPairsInt, recobEigenVals, recobEigenVecs, meanPos);
344  }
345  else
346  {
348  }
349 
350  return;
351 }
352 
353 void PrincipalComponentsAlg::PCAAnalysis_2D(const reco::HitPairListPtr& hitPairVector, reco::PrincipalComponents& pca, bool updateAvePos) const
354 {
355  // Once an axis has been found our goal is to refine it by using only the 2D hits
356  // We'll get 3D information for each of these by using the axis as a reference and use
357  // the point of closest approach as the 3D position
358 
359  // Define elements of our covariance matrix
360  double xi2(0.);
361  double xiyi(0.);
362  double xizi(0.);
363  double yi2(0.);
364  double yizi(0.);
365  double zi2(0.);
366 
367  double aveHitDoca(0.);
368  TVector3 avePosUpdate(0.,0.,0.);
369  int nHits(0);
370 
371  // Recover existing line parameters for current cluster
372  const reco::PrincipalComponents& inputPca = pca;
373  TVector3 avePosition(inputPca.getAvePosition()[0], inputPca.getAvePosition()[1], inputPca.getAvePosition()[2]);
374  TVector3 axisDirVec(inputPca.getEigenVectors()[0][0], inputPca.getEigenVectors()[0][1], inputPca.getEigenVectors()[0][2]);
375 
376  // We double loop here so we can use this method for both the first time through
377  // and a second time through where we re-calculate the mean position
378  // So, we need to keep track of the poca which we do with a double vector
379  std::vector<TVector3> hitPosVec;
380 
381  // Outer loop over 3D hits
382  for (const auto& hit3D : hitPairVector)
383  {
384  // Inner loop over 2D hits
385  for (const auto& hit : hit3D->getHits())
386  {
387  // Step one is to set up to determine the point of closest approach of this 2D hit to
388  // the cluster's current axis.
389  // Get this wire's geometry object
390  const geo::WireID& hitID = hit->getHit().WireID();
391  const geo::WireGeo& wire_geom = m_geometry->WireIDToWireGeo(hitID);
392 
393  // From this, get the parameters of the line for the wire
394  double wirePosArr[3] = {0.,0.,0.};
395  wire_geom.GetCenter(wirePosArr);
396 
397  TVector3 wireCenter(wirePosArr[0], wirePosArr[1], wirePosArr[2]);
398  TVector3 wireDirVec(wire_geom.Direction());
399 
400  // Correct the wire position in x to set to correspond to the drift time
401  double hitPeak(hit->getHit().PeakTime());
402 
403  TVector3 wirePos(m_detector->ConvertTicksToX(hitPeak, hitID.Plane, hitID.TPC, hitID.Cryostat), wireCenter[1], wireCenter[2]);
404 
405  // Compute the wire plane normal for this view
406  TVector3 xAxis(1.,0.,0.);
407  TVector3 planeNormal = xAxis.Cross(wireDirVec); // This gives a normal vector in +z for a Y wire
408 
409  double docaInPlane(wirePos[0] - avePosition[0]);
410  double arcLenToPlane(0.);
411  double cosAxisToPlaneNormal = axisDirVec.Dot(planeNormal);
412 
413  TVector3 axisPlaneIntersection = wirePos;
414  TVector3 hitPosTVec = wirePos;
415 
416  if (fabs(cosAxisToPlaneNormal) > 0.)
417  {
418  TVector3 deltaPos = wirePos - avePosition;
419 
420  arcLenToPlane = deltaPos.Dot(planeNormal) / cosAxisToPlaneNormal;
421  axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
422  docaInPlane = wirePos[0] - axisPlaneIntersection[0];
423 
424  TVector3 axisToInter = axisPlaneIntersection - wirePos;
425  double arcLenToDoca = axisToInter.Dot(wireDirVec);
426 
427  hitPosTVec += arcLenToDoca * wireDirVec;
428  }
429 
430  // Get a vector from the wire position to our cluster's current average position
431  TVector3 wVec = avePosition - wirePos;
432 
433  // Get the products we need to compute the arc lengths to the distance of closest approach
434  double a(axisDirVec.Dot(axisDirVec));
435  double b(axisDirVec.Dot(wireDirVec));
436  double c(wireDirVec.Dot(wireDirVec));
437  double d(axisDirVec.Dot(wVec));
438  double e(wireDirVec.Dot(wVec));
439 
440  double den(a*c - b*b);
441  double arcLen1(0.);
442  double arcLen2(0.);
443 
444  // Parallel lines is a special case
445  if (fabs(den) > m_parallel)
446  {
447  arcLen1 = (b*e - c*d) / den;
448  arcLen2 = (a*e - b*d) / den;
449  }
450  else
451  {
452  mf::LogDebug("Cluster3D") << "** Parallel lines, need a solution here" << std::endl;
453  break;
454  }
455 
456  // Now get the hit position we'll use for the pca analysis
457  //double hitPos[] = {wirePos[0] + arcLen2 * wireDirVec[0],
458  // wirePos[1] + arcLen2 * wireDirVec[1],
459  // wirePos[2] + arcLen2 * wireDirVec[2]};
460  //double axisPos[] = {avePosition[0] + arcLen1 * axisDirVec[0],
461  // avePosition[1] + arcLen1 * axisDirVec[1],
462  // avePosition[2] + arcLen1 * axisDirVec[2]};
463  TVector3 hitPos = wirePos + arcLen2 * wireDirVec;
464  TVector3 axisPos = avePosition + arcLen1 * axisDirVec;
465  double deltaX = hitPos[0] - axisPos[0];
466  double deltaY = hitPos[1] - axisPos[1];
467  double deltaZ = hitPos[2] - axisPos[2];
468  double doca2 = deltaX*deltaX + deltaY*deltaY + deltaZ*deltaZ;
469  double doca = sqrt(doca2);
470 
471  docaInPlane = doca;
472 
473  aveHitDoca += fabs(docaInPlane);
474 
475  //TVector3 deltaPos = hitPos - hitPosTVec;
476  //double deltaDoca = doca - docaInPlane;
477 
478  //if (fabs(deltaPos[0]) > 1. || fabs(deltaPos[1]) > 1. || fabs(deltaPos[2]) > 1. || fabs(deltaDoca) > 2.)
479  //{
480  // std::cout << "**************************************************************************************" << std::endl;
481  // std::cout << "Diff in doca: " << deltaPos[0] << "," << deltaPos[1] << "," << deltaPos[2] << ", deltaDoca: " << deltaDoca << std::endl;
482  // std::cout << "-- HitPosTVec: " << hitPosTVec[0] << "," << hitPosTVec[1] << "," << hitPosTVec[2] << std::endl;
483  // std::cout << "-- hitPos: " << hitPos[0] << "," << hitPos[1] << "," << hitPos[2] << std::endl;
484  // std::cout << "-- WirePos: " << wirePos[0] << "," << wirePos[1] << "," << wirePos[2] << std::endl;
485  //}
486 
487  // Set the hit's doca and arclen
488  hit->setDocaToAxis(fabs(docaInPlane));
489  hit->setArcLenToPoca(arcLenToPlane);
490 
491  // If this point is considered an outlier then we skip
492  // the accumulation for average position and covariance
493  if (hit->getStatusBits() & 0x80)
494  {
495  continue;
496  }
497 
498  //hitPosTVec = hitPos;
499 
500  avePosUpdate += hitPosTVec;
501 
502  hitPosVec.push_back(hitPosTVec);
503 
504  nHits++;
505  }
506  }
507 
508  // Get updated average position
509  avePosUpdate[0] /= double(nHits);
510  avePosUpdate[1] /= double(nHits);
511  avePosUpdate[2] /= double(nHits);
512 
513  // Get the average hit doca
514  aveHitDoca /= double(nHits);
515 
516  if (updateAvePos)
517  {
518  avePosition = avePosUpdate;
519  }
520 
521  // Now loop through the hits and build out the covariance matrix
522  for(auto& hitPos : hitPosVec)
523  {
524  // And increment the values in the covariance matrix
525  double x = hitPos[0] - avePosition[0];
526  double y = hitPos[1] - avePosition[1];
527  double z = hitPos[2] - avePosition[2];
528 
529  xi2 += x * x;
530  xiyi += x * y;
531  xizi += x * z;
532  yi2 += y * y;
533  yizi += y * z;
534  zi2 += z * z;
535  }
536 
537  // Accumulation done, now do the actual work
538 
539  // Create the actual matrix
540  TMatrixD sigma(3, 3);
541 
542  sigma(0,0) = xi2;
543  sigma(0,1) = sigma(1,0) = xiyi;
544  sigma(0,2) = sigma(2,0) = xizi;
545  sigma(1,1) = yi2;
546  sigma(1,2) = sigma(2,1) = yizi;
547  sigma(2,2) = zi2;
548 
549  // Scale by number of pairs
550  sigma *= (1./(nHits - 1.));
551 
552  // Set up the SVD
553  TDecompSVD rootSVD(sigma);
554 
555  // run the decomposition
556  bool svdOk(false);
557 
558  try
559  {
560  svdOk = rootSVD.Decompose();
561  }
562  catch(...)
563  {
564  svdOk = false;
565  mf::LogDebug("Cluster3D") << "PCA decompose failure, nhits = " << nHits << std::endl;
566  }
567 
568  if (svdOk)
569  {
570  // Extract results
571  TVectorD eigenVals = rootSVD.GetSig();
572  TMatrixD eigenVecs = rootSVD.GetU();
573 
574  // Get the eigen values
575  double recobEigenVals[] = {eigenVals[0], eigenVals[1], eigenVals[2]};
576 
577  // Grab the principle axes
579  std::vector<double> tempVec;
580 
581  // Get the first column vector
582  tempVec.push_back(eigenVecs(0, 0));
583  tempVec.push_back(eigenVecs(1, 0));
584  tempVec.push_back(eigenVecs(2, 0));
585  recobEigenVecs.push_back(tempVec);
586 
587  // Now the second
588  tempVec.clear();
589  tempVec.push_back(eigenVecs(0, 1));
590  tempVec.push_back(eigenVecs(1, 1));
591  tempVec.push_back(eigenVecs(2, 1));
592  recobEigenVecs.push_back(tempVec);
593 
594  // And the last
595  tempVec.clear();
596  tempVec.push_back(eigenVecs(0, 2));
597  tempVec.push_back(eigenVecs(1, 2));
598  tempVec.push_back(eigenVecs(2, 2));
599  recobEigenVecs.push_back(tempVec);
600 
601  // Save the average position
602  double avePosToSave[] = {avePosition[0],avePosition[1],avePosition[2]};
603 
604  // Store away
605  pca = reco::PrincipalComponents(svdOk, nHits, recobEigenVals, recobEigenVecs, avePosToSave, aveHitDoca);
606  }
607  else
608  {
610  }
611 
612  return;
613 }
614 
616  const reco::PrincipalComponents& pca) const
617 {
618  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
619  // any outliers. Basically, any hit outside a scaled range of the average doca from the
620  // first pass is marked by setting the bit in the status word.
621 
622  // We'll need the current PCA axis to determine doca and arclen
623  TVector3 avePosition(pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
624  TVector3 axisDirVec(pca.getEigenVectors()[0][0], pca.getEigenVectors()[0][1], pca.getEigenVectors()[0][2]);
625 
626  // We want to keep track of the average
627  double aveDoca3D(0.);
628 
629  // Outer loop over views
630  for (const auto* clusterHit3D : hitPairVector)
631  {
632  // Always reset the existing status bit
633  clusterHit3D->clearStatusBits(0x80);
634 
635  // Now we need to calculate the doca and poca...
636  // Start by getting this hits position
637  TVector3 clusPos(clusterHit3D->getPosition()[0],clusterHit3D->getPosition()[1],clusterHit3D->getPosition()[2]);
638 
639  // Form a TVector from this to the cluster average position
640  TVector3 clusToHitVec = clusPos - avePosition;
641 
642  // With this we can get the arclength to the doca point
643  double arclenToPoca = clusToHitVec.Dot(axisDirVec);
644 
645  // Get the coordinates along the axis for this point
646  TVector3 docaPos = avePosition + arclenToPoca * axisDirVec;
647 
648  // Now get doca and poca
649  TVector3 docaPosToClusPos = clusPos - docaPos;
650  double docaToAxis = docaPosToClusPos.Mag();
651 
652  aveDoca3D += docaToAxis;
653 
654  // Ok, set the values in the hit
655  clusterHit3D->setDocaToAxis(docaToAxis);
656  clusterHit3D->setArclenToPoca(arclenToPoca);
657  }
658 
659  // Compute the average and store
660  aveDoca3D /= double(hitPairVector.size());
661 
662  pca.setAveHitDoca(aveDoca3D);
663 
664  return;
665 }
666 
668  const reco::PrincipalComponents& pca ) const
669 {
670  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
671  // any outliers. Basically, any hit outside a scaled range of the average doca from the
672  // first pass is marked by setting the bit in the status word.
673 
674  // We'll need the current PCA axis to determine doca and arclen
675  TVector3 avePosition(pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
676  TVector3 axisDirVec(pca.getEigenVectors()[0][0], pca.getEigenVectors()[0][1], pca.getEigenVectors()[0][2]);
677 
678  // Recover the principle eigen value for range constraints
679  double maxArcLen = 4.*sqrt(pca.getEigenValues()[0]);
680 
681  // We want to keep track of the average
682  double aveHitDoca(0.);
683 
684  // Outer loop over views
685  for (const auto* hit : hit2DListPtr)
686  {
687  // Step one is to set up to determine the point of closest approach of this 2D hit to
688  // the cluster's current axis. We do that by finding the point of intersection of the
689  // cluster's axis with a plane defined by the wire the hit is associated with.
690  // Get this wire's geometry object
691  const geo::WireID& hitID = hit->getHit().WireID();
692  const geo::WireGeo& wire_geom = m_geometry->WireIDToWireGeo(hitID);
693 
694  // From this, get the parameters of the line for the wire
695  double wirePosArr[3] = {0.,0.,0.};
696  wire_geom.GetCenter(wirePosArr);
697 
698  TVector3 wireCenter(wirePosArr[0], wirePosArr[1], wirePosArr[2]);
699  TVector3 wireDirVec(wire_geom.Direction());
700 
701  // Correct the wire position in x to set to correspond to the drift time
702  TVector3 wirePos(hit->getXPosition(), wireCenter[1], wireCenter[2]);
703 
704  // Compute the wire plane normal for this view
705  TVector3 xAxis(1.,0.,0.);
706  TVector3 planeNormal = xAxis.Cross(wireDirVec); // This gives a normal vector in +z for a Y wire
707 
708  double arcLenToPlane(0.);
709  double docaInPlane(wirePos[0] - avePosition[0]);
710  double cosAxisToPlaneNormal = axisDirVec.Dot(planeNormal);
711 
712  TVector3 axisPlaneIntersection = wirePos;
713 
714  // If current cluster axis is not parallel to wire plane then find intersection point
715  if (fabs(cosAxisToPlaneNormal) > 0.)
716  {
717  TVector3 deltaPos = wirePos - avePosition;
718 
719  arcLenToPlane = std::min(deltaPos.Dot(planeNormal) / cosAxisToPlaneNormal, maxArcLen);
720  axisPlaneIntersection = avePosition + arcLenToPlane * axisDirVec;
721  TVector3 axisToInter = axisPlaneIntersection - wirePos;
722  double arcLenToDoca = axisToInter.Dot(wireDirVec);
723 
724  // If the arc length along the wire to the poca is outside the TPC then reset
725  if (fabs(arcLenToDoca) > wire_geom.HalfL()) arcLenToDoca = wire_geom.HalfL();
726 
727  // If we were successful in getting to the wire plane then the doca is simply the
728  // difference in x coordinates... but we hvae to worry about the special cases so
729  // we calculate a 3D doca based on arclengths above...
730  TVector3 docaVec = axisPlaneIntersection - (wirePos + arcLenToDoca * wireDirVec);
731  docaInPlane = docaVec.Mag();
732  }
733 
734  aveHitDoca += fabs(docaInPlane);
735 
736  // Set the hit's doca and arclen
737  hit->setDocaToAxis(fabs(docaInPlane));
738  hit->setArcLenToPoca(arcLenToPlane);
739  }
740 
741  // Compute the average and store
742  aveHitDoca /= double(hit2DListPtr.size());
743 
744  pca.setAveHitDoca(aveHitDoca);
745 
746  return;
747 }
748 
751  double maxDocaAllowed) const
752 {
753  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
754  // any outliers. Basically, any hit outside a scaled range of the average doca from the
755  // first pass is marked by setting the bit in the status word.
756 
757  // First get the average doca scaled by some appropriate factor
758  int numRejHits(0);
759 
760  // Outer loop over views
761  for (const auto& hit3D : hitPairVector)
762  {
763  // Inner loop over hits in this view
764  for (const auto& hit : hit3D->getHits())
765  {
766  // Always reset the existing status bit
767  hit->clearStatusBits(0x80);
768 
769  if (hit->getDocaToAxis() > maxDocaAllowed)
770  {
771  hit->setStatusBit(0x80);
772  numRejHits++;
773  }
774  }
775  }
776 
777  return numRejHits;
778 }
779 
781  const reco::PrincipalComponents& pca,
782  double maxDocaAllowed) const
783 {
784  // Our mission, should we choose to accept it, is to scan through the 2D hits and reject
785  // any outliers. Basically, any hit outside a scaled range of the average doca from the
786  // first pass is marked by setting the bit in the status word.
787 
788  // First get the average doca scaled by some appropriate factor
789  int numRejHits(0);
790 
791  // We'll need the current PCA axis to determine doca and arclen
792  TVector3 avePosition(pca.getAvePosition()[0], pca.getAvePosition()[1], pca.getAvePosition()[2]);
793  TVector3 axisDirVec(pca.getEigenVectors()[0][0], pca.getEigenVectors()[0][1], pca.getEigenVectors()[0][2]);
794 
795  // Outer loop over views
796  for (const auto* clusterHit3D : hitPairVector)
797  {
798  // Always reset the existing status bit
799  clusterHit3D->clearStatusBits(0x80);
800 
801  // Now we need to calculate the doca and poca...
802  // Start by getting this hits position
803  TVector3 clusPos(clusterHit3D->getPosition()[0],clusterHit3D->getPosition()[1],clusterHit3D->getPosition()[2]);
804 
805  // Form a TVector from this to the cluster average position
806  TVector3 clusToHitVec = clusPos - avePosition;
807 
808  // With this we can get the arclength to the doca point
809  double arclenToPoca = clusToHitVec.Dot(axisDirVec);
810 
811  // Get the coordinates along the axis for this point
812  TVector3 docaPos = avePosition + arclenToPoca * axisDirVec;
813 
814  // Now get doca and poca
815  TVector3 docaPosToClusPos = clusPos - docaPos;
816  double docaToAxis = docaPosToClusPos.Mag();
817 
818  // Ok, set the values in the hit
819  clusterHit3D->setDocaToAxis(docaToAxis);
820  clusterHit3D->setArclenToPoca(arclenToPoca);
821 
822  // Check to see if this is a keeper
823  if (clusterHit3D->getDocaToAxis() > maxDocaAllowed)
824  {
825  clusterHit3D->setStatusBit(0x80);
826  numRejHits++;
827  }
828  }
829 
830  return numRejHits;
831 }
832 
833 
834 
835 } // namespace lar_cluster3d
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
void setAveHitDoca(double doca) const
Definition: Cluster3D.h:252
void PCAAnalysis_2D(const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, bool updateAvePos=false) const
std::list< const reco::ClusterHit2D * > Hit2DListPtr
export some data structure definitions
Definition: Cluster3D.h:334
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
PrincipalComponentsAlg(fhicl::ParameterSet const &pset)
Constructor.
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
geo::WireID WireID() const
Definition: Hit.h:233
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
int PCAAnalysis_reject2DOutliers(const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, double aveHitDoca) const
float getDocaToAxis() const
Definition: Cluster3D.h:169
art framework interface to geometry description
const recob::Hit * getHit() const
Definition: Cluster3D.h:78
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
const double e
void PCAAnalysis(const reco::HitPairListPtr &hitPairVector, reco::PrincipalComponents &pca, double doca3DScl=3.) const
Run the Principal Components Analysis.
const double a
T get(std::string const &key) const
Definition: ParameterSet.h:271
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
double m_parallel
means lines are parallel
int PCAAnalysis_reject3DOutliers(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca, double aveHitDoca) const
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Vector Direction() const
Definition: WireGeo.h:587
Detector simulation of raw signals on wires.
Encapsulate the geometry of a wire.
double HalfL() const
Returns half the length of the wire [cm].
Definition: WireGeo.h:128
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
Declaration of signal hit object.
void getHit2DPocaToAxis(const TVector3 &axisPos, const TVector3 &axisDir, const reco::ClusterHit2D *hit2D, TVector3 &poca, double &arcLenAxis, double &arcLenWire, double &doca)
This is used to get the poca, doca and arclen along cluster axis to 2D hit.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
Encapsulate the construction of a single detector plane.
const float getAveHitDoca() const
Definition: Cluster3D.h:249
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Eigen::Matrix3f EigenVectors
Definition: Cluster3D.h:227
float getArclenToPoca() const
Definition: Cluster3D.h:170
static bool * b
Definition: config.cpp:1043
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
list x
Definition: train.py:276
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
void PCAAnalysis_calc2DDocas(const reco::Hit2DListPtr &hit2DVector, const reco::PrincipalComponents &pca) const
const detinfo::DetectorProperties * m_detector
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247
QTextStream & endl(QTextStream &s)
WireGeo const & WireIDToWireGeo(geo::WireID const &wireid) const