All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TwoViewCosmicRayRemovalTool.cc
Go to the documentation of this file.
1 /**
2  * @file larpandoracontent/LArThreeDReco/LArTransverseTrackMatching/TwoViewCosmicRayRemovalTool.cc
3  *
4  * @brief Implementation of the two view cosmic removal tool class.
5  *
6  * $Log: $
7  */
8 
10 #include "Pandora/AlgorithmHeaders.h"
11 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 TwoViewCosmicRayRemovalTool::TwoViewCosmicRayRemovalTool() :
23  m_minSeparation(2.f),
24  m_slidingFitWindow(10000),
25  m_distanceToLine(0.5f),
26  m_minContaminationLength(3.f),
27  m_maxDistanceToHit(1.f),
28  m_minRemnantClusterSize(3),
29  m_maxDistanceToTrack(2.f)
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
36 {
37  m_pParentAlgorithm = pAlgorithm;
38 
39  if (PandoraContentApi::GetSettings(*m_pParentAlgorithm)->ShouldDisplayAlgorithmInfo())
40  std::cout << "----> Running Algorithm Tool: " << this->GetInstanceName() << ", " << this->GetType() << std::endl;
41 
42  bool changesMade(false);
43 
44  ClusterVector sortedKeyClusters;
45  overlapMatrix.GetSortedKeyClusters(sortedKeyClusters);
46 
47  ClusterSet usedKeyClusters;
48  for (const Cluster *const pKeyCluster : sortedKeyClusters)
49  {
50  if (usedKeyClusters.count(pKeyCluster))
51  continue;
52 
53  MatrixType::ElementList elementList;
54  overlapMatrix.GetConnectedElements(pKeyCluster, true, elementList);
55 
56  for (const MatrixType::Element &element : elementList)
57  usedKeyClusters.insert(element.GetCluster1());
58 
59  const bool changesMadeInIteration = this->RemoveCosmicRayHits(elementList);
60 
61  changesMade = (changesMade ? changesMade : changesMadeInIteration);
62  }
63 
64  return changesMade;
65 }
66 
67 //------------------------------------------------------------------------------------------------------------------------------------------
68 
70 {
71  ClusterSet modifiedClusters, checkedClusters;
72 
73  for (const MatrixType::Element &element : elementList)
74  {
75  for (const HitType hitType : m_pParentAlgorithm->GetHitTypeVector())
76  {
77  const Cluster *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType));
78 
79  if (checkedClusters.count(pDeltaRayCluster))
80  continue;
81 
82  // ATTN: The underlying matrix will update during this loop, do not proceed if element has been modified
83  if ((modifiedClusters.count(element.GetCluster1())) || (modifiedClusters.count(element.GetCluster2())))
84  continue;
85 
86  if (!this->PassElementChecks(element, hitType))
87  continue;
88 
89  if (!this->IsContaminated(element, hitType))
90  continue;
91 
92  if (!this->IsBestElement(element, hitType, elementList))
93  continue;
94 
95  checkedClusters.insert(pDeltaRayCluster);
96 
97  // Attempt to pull delta ray hits out of cluster
98  CaloHitList deltaRayHits;
99  this->CreateSeed(element, hitType, deltaRayHits);
100 
101  if (deltaRayHits.empty())
102  continue;
103 
104  // ATTN: If seed cannot be grown, abort
105  CaloHitList deltaRayRemnantHits;
106  if (this->GrowSeed(element, hitType, deltaRayHits, deltaRayRemnantHits) != STATUS_CODE_SUCCESS)
107  continue;
108 
109  if (deltaRayHits.size() == pDeltaRayCluster->GetNCaloHits())
110  continue;
111 
112  modifiedClusters.insert(pDeltaRayCluster);
113 
114  this->SplitDeltaRayCluster(element, hitType, deltaRayHits, deltaRayRemnantHits);
115  }
116  }
117 
118  return !modifiedClusters.empty();
119 }
120 
121 //------------------------------------------------------------------------------------------------------------------------------------------
122 
123 bool TwoViewCosmicRayRemovalTool::PassElementChecks(const MatrixType::Element &element, const HitType hitType) const
124 {
125  // ATTN: Avoid endpoints, topological michel reconstruction is very ambiguous
126  if (this->IsMuonEndpoint(element, hitType))
127  return false;
128 
129  const Cluster *pMuonCluster(nullptr), *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType));
130 
131  if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
132  return false;
133 
134  const float separation(LArClusterHelper::GetClosestDistance(pDeltaRayCluster, pMuonCluster));
135 
136  return separation <= m_minSeparation;
137 }
138 
139 //------------------------------------------------------------------------------------------------------------------------------------------
140 
141 bool TwoViewCosmicRayRemovalTool::IsMuonEndpoint(const MatrixType::Element &element, const HitType hitType) const
142 {
143  for (const HitType &otherHitType : m_pParentAlgorithm->GetHitTypeVector())
144  {
145  if (otherHitType == hitType)
146  continue;
147 
148  const Cluster *pMuonCluster(nullptr), *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType));
149 
150  if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
151  return false;
152 
153  float xMinDR(+std::numeric_limits<float>::max()), xMaxDR(-std::numeric_limits<float>::max());
154  pDeltaRayCluster->GetClusterSpanX(xMinDR, xMaxDR);
155 
156  float xMinCR(-std::numeric_limits<float>::max()), xMaxCR(+std::numeric_limits<float>::max());
157  pMuonCluster->GetClusterSpanX(xMinCR, xMaxCR);
158 
159  if ((xMinDR < xMinCR) || (xMaxDR > xMaxCR))
160  return false;
161 
162  float zMinDR(+std::numeric_limits<float>::max()), zMaxDR(-std::numeric_limits<float>::max());
163  pDeltaRayCluster->GetClusterSpanZ(xMinDR, xMaxDR, zMinDR, zMaxDR);
164 
165  float zMinCR(-std::numeric_limits<float>::max()), zMaxCR(+std::numeric_limits<float>::max());
166  pMuonCluster->GetClusterSpanZ(xMinCR, xMaxCR, zMinCR, zMaxCR);
167 
168  if ((zMinDR < zMinCR) || (zMaxDR > zMaxCR))
169  return false;
170  }
171 
172  return true;
173 }
174 
175 //------------------------------------------------------------------------------------------------------------------------------------------
176 
177 bool TwoViewCosmicRayRemovalTool::IsContaminated(const MatrixType::Element &element, const HitType hitType) const
178 {
179  const Cluster *pMuonCluster(nullptr), *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType));
180 
181  if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
182  return false;
183 
184  const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
185  const TwoDSlidingFitResult slidingFitResult(pMuonCluster, m_slidingFitWindow, slidingFitPitch);
186 
187  CartesianVector muonDirection(0.f, 0.f, 0.f);
188  slidingFitResult.GetGlobalDirection(slidingFitResult.GetLayerFitResultMap().begin()->second.GetGradient(), muonDirection);
189 
190  CartesianVector deltaRayVertex(0.f, 0.f, 0.f), muonVertex(0.f, 0.f, 0.f);
191  LArClusterHelper::GetClosestPositions(pDeltaRayCluster, pMuonCluster, deltaRayVertex, muonVertex);
192 
193  // Find furthest point in DR cluster that lies on the projected muon track
194  float furthestSeparation(0.f);
195  CartesianVector extendedPoint(0.f, 0.f, 0.f);
196 
197  CaloHitList deltaRayHitList;
198  pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
199 
200  for (const CaloHit *const pCaloHit : deltaRayHitList)
201  {
202  const CartesianVector &position(pCaloHit->GetPositionVector());
203  const float separation((position - muonVertex).GetMagnitude());
204 
205  if (separation > furthestSeparation)
206  {
207  if (this->IsCloseToLine(position, muonVertex, muonVertex + muonDirection, m_distanceToLine))
208  {
209  furthestSeparation = separation;
210  extendedPoint = position;
211  }
212  }
213  }
214 
215  // Check if delta ray has significant track contamination
216  if (furthestSeparation < m_minContaminationLength)
217  return false;
218 
219  // ATTN: Avoid cases where opening angle is small - it is easy to make mistakes in these instances
220  CaloHitList muonHitList;
221  pMuonCluster->GetOrderedCaloHitList().FillCaloHitList(muonHitList);
222 
223  for (const CaloHit *const pCaloHit : muonHitList)
224  {
225  if (this->IsInLineSegment(deltaRayVertex, extendedPoint, pCaloHit->GetPositionVector()))
226  return false;
227  }
228 
229  return true;
230 }
231 
232 //------------------------------------------------------------------------------------------------------------------------------------------
233 
235  const CartesianVector &hitPosition, const CartesianVector &lineStart, const CartesianVector &lineEnd, const float distanceToLine) const
236 {
237  CartesianVector lineDirection(lineStart - lineEnd);
238  lineDirection = lineDirection.GetUnitVector();
239 
240  const float transverseDistanceFromLine(lineDirection.GetCrossProduct(hitPosition - lineStart).GetMagnitude());
241 
242  if (transverseDistanceFromLine > distanceToLine)
243  return false;
244 
245  return true;
246 }
247 
248 //------------------------------------------------------------------------------------------------------------------------------------------
249 
251  const CartesianVector &lowerBoundary, const CartesianVector &upperBoundary, const CartesianVector &point) const
252 {
253  const float segmentBoundaryGradient = (-1.f) * (upperBoundary.GetX() - lowerBoundary.GetX()) / (upperBoundary.GetZ() - lowerBoundary.GetZ());
254  const float xPointOnUpperLine((point.GetZ() - upperBoundary.GetZ()) / segmentBoundaryGradient + upperBoundary.GetX());
255  const float xPointOnLowerLine((point.GetZ() - lowerBoundary.GetZ()) / segmentBoundaryGradient + lowerBoundary.GetX());
256 
257  if (std::fabs(xPointOnUpperLine - point.GetX()) < std::numeric_limits<float>::epsilon())
258  return true;
259 
260  if (std::fabs(xPointOnLowerLine - point.GetX()) < std::numeric_limits<float>::epsilon())
261  return true;
262 
263  if ((point.GetX() > xPointOnUpperLine) && (point.GetX() > xPointOnLowerLine))
264  return false;
265 
266  if ((point.GetX() < xPointOnUpperLine) && (point.GetX() < xPointOnLowerLine))
267  return false;
268 
269  return true;
270 }
271 
272 //------------------------------------------------------------------------------------------------------------------------------------------
273 
274 bool TwoViewCosmicRayRemovalTool::IsBestElement(const MatrixType::Element &element, const HitType hitType, const MatrixType::ElementList &elementList) const
275 {
276  const float chiSquared(element.GetOverlapResult().GetReducedChiSquared());
277  const unsigned int hitSum(element.GetCluster1()->GetNCaloHits() + element.GetCluster2()->GetNCaloHits());
278 
279  for (const MatrixType::Element &testElement : elementList)
280  {
281  if (m_pParentAlgorithm->GetCluster(testElement, hitType) != m_pParentAlgorithm->GetCluster(element, hitType))
282  continue;
283 
284  if ((testElement.GetCluster1() == element.GetCluster1()) && (testElement.GetCluster2() == element.GetCluster2()))
285  continue;
286 
287  const float testChiSquared(testElement.GetOverlapResult().GetReducedChiSquared());
288  const unsigned int testHitSum(testElement.GetCluster1()->GetNCaloHits() + testElement.GetCluster2()->GetNCaloHits());
289 
290  if ((testHitSum < hitSum) || ((testHitSum == hitSum) && (testChiSquared > chiSquared)))
291  continue;
292 
293  if (this->PassElementChecks(testElement, hitType))
294  return false;
295  }
296 
297  return true;
298 }
299 
300 //------------------------------------------------------------------------------------------------------------------------------------------
301 
302 void TwoViewCosmicRayRemovalTool::CreateSeed(const MatrixType::Element &element, const HitType hitType, CaloHitList &collectedHits) const
303 {
304  // To avoid fluctuations, parameterise the muon track
305  CartesianVector positionOnMuon(0.f, 0.f, 0.f), muonDirection(0.f, 0.f, 0.f);
306  if (m_pParentAlgorithm->ParameteriseMuon(element.GetOverlapResult().GetCommonMuonPfoList().front(),
307  m_pParentAlgorithm->GetCluster(element, hitType), positionOnMuon, muonDirection) != STATUS_CODE_SUCCESS)
308  return;
309 
310  CartesianPointVector deltaRayProjectedPositions;
311  if (this->ProjectDeltaRayPositions(element, hitType, deltaRayProjectedPositions) != STATUS_CODE_SUCCESS)
312  return;
313 
314  CaloHitList deltaRayHitList;
315  m_pParentAlgorithm->GetCluster(element, hitType)->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
316 
317  CaloHitVector deltaRayHitVector(deltaRayHitList.begin(), deltaRayHitList.end());
318  std::sort(deltaRayHitVector.begin(), deltaRayHitVector.end(), LArClusterHelper::SortHitsByPosition);
319 
320  for (const CaloHit *const pCaloHit : deltaRayHitVector)
321  {
322  const CartesianVector &position(pCaloHit->GetPositionVector());
323 
324  for (const CartesianVector &projectedPosition : deltaRayProjectedPositions)
325  {
326  const float distanceToProjectionSquared((position - projectedPosition).GetMagnitudeSquared());
327 
328  if (distanceToProjectionSquared < m_maxDistanceToHit * m_maxDistanceToHit)
329  {
330  // To prevent gappy seeds
331  if ((!collectedHits.empty()) && (LArClusterHelper::GetClosestDistance(position, collectedHits) > m_maxDistanceToHit))
332  continue;
333 
334  const float distanceToMuonHitsSquared(muonDirection.GetCrossProduct(position - positionOnMuon).GetMagnitudeSquared());
335 
336  if (distanceToMuonHitsSquared < m_maxDistanceToHit * m_maxDistanceToHit)
337  continue;
338 
339  collectedHits.push_back(pCaloHit);
340 
341  break;
342  }
343  }
344  }
345 }
346 
347 //------------------------------------------------------------------------------------------------------------------------------------------
348 
350  const MatrixType::Element &element, const HitType hitType, CaloHitList &deltaRayHits, CaloHitList &remnantHits) const
351 {
352  const Cluster *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType)), *pMuonCluster(nullptr);
353 
354  if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
355  return STATUS_CODE_NOT_FOUND;
356 
357  // To avoid fluctuatuions, parameterise the muon track
358  CartesianVector positionOnMuon(0.f, 0.f, 0.f), muonDirection(0.f, 0.f, 0.f);
359  if (m_pParentAlgorithm->ParameteriseMuon(element.GetOverlapResult().GetCommonMuonPfoList().front(), pDeltaRayCluster, positionOnMuon,
360  muonDirection) != STATUS_CODE_SUCCESS)
361  return STATUS_CODE_NOT_FOUND;
362 
363  // Identify delta ray hits
364  this->CollectHitsFromDeltaRay(positionOnMuon, muonDirection, pDeltaRayCluster, m_maxDistanceToHit, true, deltaRayHits, deltaRayHits);
365 
366  // Identify remnant hits
367  this->CollectHitsFromDeltaRay(positionOnMuon, muonDirection, pDeltaRayCluster, m_maxDistanceToHit, false, deltaRayHits, remnantHits);
368 
369  return STATUS_CODE_SUCCESS;
370 }
371 
372 //------------------------------------------------------------------------------------------------------------------------------------------
373 
374 void TwoViewCosmicRayRemovalTool::CollectHitsFromDeltaRay(const CartesianVector &positionOnMuon, const CartesianVector &muonDirection,
375  const Cluster *const pDeltaRayCluster, const float &minDistanceFromMuon, const bool demandCloserToCollected,
376  const CaloHitList &protectedHits, CaloHitList &collectedHits) const
377 {
378  CaloHitList deltaRayHitList;
379  pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
380 
381  bool hitsAdded(false);
382 
383  do
384  {
385  hitsAdded = false;
386 
387  for (const CaloHit *const pCaloHit : deltaRayHitList)
388  {
389  if (std::find(protectedHits.begin(), protectedHits.end(), pCaloHit) != protectedHits.end())
390  continue;
391 
392  if (std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end())
393  continue;
394 
395  const float distanceToCollectedHits(collectedHits.empty()
396  ? std::numeric_limits<float>::max()
397  : LArClusterHelper::GetClosestDistance(pCaloHit->GetPositionVector(), collectedHits));
398  const float distanceToMuonHits(muonDirection.GetCrossProduct(pCaloHit->GetPositionVector() - positionOnMuon).GetMagnitude());
399 
400  if ((distanceToMuonHits < minDistanceFromMuon) || (demandCloserToCollected && (distanceToCollectedHits > distanceToMuonHits)))
401  continue;
402 
403  collectedHits.push_back(pCaloHit);
404  hitsAdded = true;
405  }
406  } while (hitsAdded);
407 }
408 
409 //------------------------------------------------------------------------------------------------------------------------------------------
410 
412  const MatrixType::Element &element, const HitType hitType, CaloHitList &collectedHits, CaloHitList &deltaRayRemnantHits) const
413 {
414  const Cluster *const pDeltaRayCluster(m_pParentAlgorithm->GetCluster(element, hitType)), *pMuonCluster(nullptr);
415 
416  if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
417  return;
418 
419  m_pParentAlgorithm->UpdateUponDeletion(pMuonCluster);
420  m_pParentAlgorithm->UpdateUponDeletion(pDeltaRayCluster);
421 
422  std::string originalListName, fragmentListName;
423  ClusterList originalClusterList(1, pDeltaRayCluster);
424 
425  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
426  PandoraContentApi::ReplaceCurrentList<Cluster>(*m_pParentAlgorithm, m_pParentAlgorithm->GetClusterListName(hitType)));
427  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
428  PandoraContentApi::InitializeFragmentation(*m_pParentAlgorithm, originalClusterList, originalListName, fragmentListName));
429 
430  CaloHitList deltaRayHitList;
431  pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
432 
433  const Cluster *pDeltaRay(nullptr), *pDeltaRayRemnant(nullptr);
434 
435  for (const CaloHit *const pCaloHit : deltaRayHitList)
436  {
437  const bool isDeltaRay(std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end());
438  const bool isDeltaRayRemnant(
439  !isDeltaRay && (std::find(deltaRayRemnantHits.begin(), deltaRayRemnantHits.end(), pCaloHit) != deltaRayRemnantHits.end()));
440  const Cluster *&pCluster(isDeltaRay ? pDeltaRay : isDeltaRayRemnant ? pDeltaRayRemnant : pMuonCluster);
441 
442  if (!pCluster)
443  {
445  parameters.m_caloHitList.push_back(pCaloHit);
446  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*m_pParentAlgorithm, parameters, pCluster));
447  }
448  else
449  {
450  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*m_pParentAlgorithm, pCluster, pCaloHit));
451  }
452  }
453 
454  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::EndFragmentation(*m_pParentAlgorithm, fragmentListName, originalListName));
455 
456  ClusterVector clusterVector;
457  PfoVector pfoVector;
458  if (pDeltaRayRemnant)
459  this->ReclusterRemnant(hitType, pMuonCluster, pDeltaRayRemnant, clusterVector, pfoVector);
460 
461  clusterVector.push_back(pMuonCluster);
462  pfoVector.push_back(element.GetOverlapResult().GetCommonMuonPfoList().front());
463  clusterVector.push_back(pDeltaRay);
464  pfoVector.push_back(nullptr);
465 
466  m_pParentAlgorithm->UpdateForNewClusters(clusterVector, pfoVector);
467 }
468 
469 //------------------------------------------------------------------------------------------------------------------------------------------
470 
471 void TwoViewCosmicRayRemovalTool::ReclusterRemnant(const HitType hitType, const Cluster *const pMuonCluster,
472  const Cluster *const pDeltaRayRemnant, ClusterVector &clusterVector, PfoVector &pfoVector) const
473 {
474  std::string caloHitListName(hitType == TPC_VIEW_U ? "CaloHitListU" : hitType == TPC_VIEW_V ? "CaloHitListV" : "CaloHitListW");
475 
476  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<CaloHit>(*m_pParentAlgorithm, caloHitListName));
477  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
478  PandoraContentApi::ReplaceCurrentList<Cluster>(*m_pParentAlgorithm, m_pParentAlgorithm->GetClusterListName(hitType)));
479  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Delete<Cluster>(*m_pParentAlgorithm, pDeltaRayRemnant));
480 
481  const ClusterList *pClusterList(nullptr);
482  std::string newClusterListName;
483  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
484  PandoraContentApi::RunClusteringAlgorithm(*m_pParentAlgorithm, m_pParentAlgorithm->GetClusteringAlgName(), pClusterList, newClusterListName));
485 
486  const ClusterList remnantClusters(pClusterList->begin(), pClusterList->end());
487 
488  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
489  PandoraContentApi::SaveList<Cluster>(*m_pParentAlgorithm, newClusterListName, m_pParentAlgorithm->GetClusterListName(hitType)));
490  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
491  PandoraContentApi::ReplaceCurrentList<Cluster>(*m_pParentAlgorithm, m_pParentAlgorithm->GetClusterListName(hitType)));
492 
493  for (const Cluster *const pRemnant : remnantClusters)
494  {
495  if (pRemnant->GetNCaloHits() < m_minRemnantClusterSize)
496  {
497  if (LArClusterHelper::GetClosestDistance(pRemnant, pMuonCluster) < m_maxDistanceToTrack)
498  {
499  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*m_pParentAlgorithm, pMuonCluster, pRemnant));
500  continue;
501  }
502  }
503 
504  clusterVector.push_back(pRemnant);
505  pfoVector.push_back(nullptr);
506  }
507 }
508 
509 //------------------------------------------------------------------------------------------------------------------------------------------
510 
512  const MatrixType::Element &element, const HitType hitType, CartesianPointVector &projectedPositions) const
513 {
514  HitTypeVector hitTypes({TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W});
515 
516  hitTypes.erase(std::find(hitTypes.begin(), hitTypes.end(), hitType));
517 
518  const Cluster *const pCluster1(m_pParentAlgorithm->GetCluster(element, hitTypes[0]));
519  const Cluster *const pCluster2(m_pParentAlgorithm->GetCluster(element, hitTypes[1]));
520 
521  if (!pCluster1 || !pCluster2)
522  return STATUS_CODE_NOT_FOUND;
523 
524  return m_pParentAlgorithm->GetProjectedPositions(pCluster1, pCluster2, projectedPositions);
525 }
526 
527 //------------------------------------------------------------------------------------------------------------------------------------------
528 
529 StatusCode TwoViewCosmicRayRemovalTool::ReadSettings(const TiXmlHandle xmlHandle)
530 {
531  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinSeparation", m_minSeparation));
532 
533  PANDORA_RETURN_RESULT_IF_AND_IF(
534  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitWindow", m_slidingFitWindow));
535 
536  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DistanceToLine", m_distanceToLine));
537 
538  PANDORA_RETURN_RESULT_IF_AND_IF(
539  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinContaminationLength", m_minContaminationLength));
540 
541  PANDORA_RETURN_RESULT_IF_AND_IF(
542  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDistanceToHit", m_maxDistanceToHit));
543 
544  PANDORA_RETURN_RESULT_IF_AND_IF(
545  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinRemnantClusterSize", m_minRemnantClusterSize));
546 
547  PANDORA_RETURN_RESULT_IF_AND_IF(
548  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDistanceToTrack", m_maxDistanceToTrack));
549 
550  return STATUS_CODE_SUCCESS;
551 }
552 
553 } // namespace lar_content
void SplitDeltaRayCluster(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits, pandora::CaloHitList &deltaRayRemnantHits) const
Fragment the delta ray cluster, refining the matched delta ray cluster perhaps creating significant r...
void CollectHitsFromDeltaRay(const pandora::CartesianVector &positionOnMuon, const pandora::CartesianVector &muonDirection, const pandora::Cluster *const pDeltaRayCluster, const float &minDistanceFromMuon, const bool demandCloserToCollected, const pandora::CaloHitList &protectedHits, pandora::CaloHitList &collectedHits) const
Collect hits from the delta ray cluster to either keep (demandCloserToCollected == true) or separate ...
unsigned int m_slidingFitWindow
The sliding fit window used in cosmic ray parameterisations.
Header file for the pfo helper class.
pandora::StatusCode GetProjectedPositions(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, pandora::CartesianPointVector &projectedPositions) const
Use two clusters from different views to calculate projected positions in the remaining third view...
bool RemoveCosmicRayHits(const MatrixType::ElementList &elementList) const
Remove hits from delta ray clusters that belong to the parent muon.
void CreateSeed(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits) const
Collect hits in the delta ray cluster that lie close to calculated projected positions forming a seed...
Header file for the cosmic ray removal tool class.
void UpdateForNewClusters(const pandora::ClusterVector &newClusterVector, const pandora::PfoVector &pfoVector)
Add a new cluster to algorithm ownership maps and, if it a delta ray cluster, to the underlying match...
const pandora::Cluster * GetCluster(const MatrixType::Element &element, const pandora::HitType hitType)
Get the address of the given hit type cluster.
static void GetClosestPositions(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, pandora::CartesianVector &position1, pandora::CartesianVector &position2)
Get pair of closest positions for a pair of clusters.
void GetSortedKeyClusters(pandora::ClusterVector &sortedKeyClusters) const
Get a sorted vector of key clusters (view 1 clusters with current implementation) ...
unsigned int m_minRemnantClusterSize
The minimum hit number of a remnant cluster for it to be considered significant.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
pandora::StatusCode ProjectDeltaRayPositions(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CartesianPointVector &projectedPositions) const
Use two views of a delta ray pfo to calculate projected positions in a given third view...
void UpdateUponDeletion(const pandora::Cluster *const pDeletedCluster)
Update to reflect cluster deletion.
HitTypeVector GetHitTypeVector()
Obtain the HitTypeVector of input views.
Header file for the geometry helper class.
bool PassElementChecks(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether element satifies simple checks.
float m_distanceToLine
The maximum perpendicular distance of a position to a line for it to be considered close...
bool IsInLineSegment(const pandora::CartesianVector &lowerBoundary, const pandora::CartesianVector &upperBoundary, const pandora::CartesianVector &point) const
Whether the projection of a given position lies on a defined line.
static bool SortHitsByPosition(const pandora::CaloHit *const pLhs, const pandora::CaloHit *const pRhs)
Sort calo hits by their position (use Z, followed by X, followed by Y)
Header file for the cluster helper class.
pandora::StatusCode GrowSeed(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits, pandora::CaloHitList &deltaRayRemantHits) const
Examine remaining hits in the delta ray cluster adding them to the delta ray seed or parent muon if a...
pandora::StatusCode ParameteriseMuon(const pandora::ParticleFlowObject *const pParentMuon, const pandora::Cluster *const pDeltaRayCluster, pandora::CartesianVector &positionOnMuon, pandora::CartesianVector &muonDirection) const
Parameterise the projection of a cosmic ray track in order to avoid poor/sparse projections.
Header file for the lar two dimensional sliding fit result class.
void ReclusterRemnant(const pandora::HitType hitType, const pandora::Cluster *const pMuonCluster, const pandora::Cluster *const pDeltaRayRemnant, pandora::ClusterVector &clusterVector, pandora::PfoVector &pfoVector) const
Reculster the remnant cluster, merging insignificant created clusters into the parent muon (if proxim...
bool IsBestElement(const MatrixType::Element &element, const pandora::HitType hitType, const MatrixType::ElementList &elementList) const
Determine whether the input element is the best to use to modify the contaminated cluster (best is de...
float m_maxDistanceToTrack
The minimum distance of an insignificant remnant cluster from the cosmic ray track for it to be merge...
TwoViewDeltaRayMatchingAlgorithm * m_pParentAlgorithm
Address of the parent matching algorithm.
float m_minContaminationLength
The minimum projected length of a delta ray cluster onto the muon track for it to be considered conta...
bool Run(TwoViewDeltaRayMatchingAlgorithm *const pAlgorithm, MatrixType &overlapMatrix)
Run the algorithm tool.
bool IsContaminated(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether the cluster under investigation has muon contamination.
void GetGlobalDirection(const float dTdL, pandora::CartesianVector &direction) const
Get global direction coordinates for given sliding linear fit gradient.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
const LayerFitResultMap & GetLayerFitResultMap() const
Get the layer fit result map.
bool IsCloseToLine(const pandora::CartesianVector &hitPosition, const pandora::CartesianVector &lineStart, const pandora::CartesianVector &lineEnd, const float distanceToLine) const
Whether a given position is close to a defined line.
fhicl::Table< sbnd::crt::CRTDetSimParams > Parameters
bool IsMuonEndpoint(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether the matched clusters suggest that the delta ray is at the endpoint of the cosmic ra...
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
Definition: reco_sbnd.fcl:182
float m_minSeparation
The minimum delta ray - parent muon cluster separation required to investigate delta ray cluster...
void GetConnectedElements(const pandora::Cluster *const pCluster, const bool ignoreUnavailable, ElementList &elementList) const
Get a list of elements connected to a specified cluster.
float m_maxDistanceToHit
The maximum distance of a hit from the cosmic ray track for it to be added into the CR...
const std::string & GetClusterListName(const pandora::HitType hitType) const
Get the cluster list name corresponding to a specified hit type.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
pandora::StatusCode GetMuonCluster(const pandora::PfoList &commonMuonPfoList, const pandora::HitType hitType, const pandora::Cluster *&pMuonCluster) const
Return the cluster of the common cosmic ray pfo in a given view (function demands there to be only on...
const std::string & GetClusteringAlgName() const
Get the name of the clustering algorithm to be used to recluster created delta ray remnants...
BEGIN_PROLOG could also be cout
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.