LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRaySplittingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 CosmicRaySplittingAlgorithm::CosmicRaySplittingAlgorithm() :
23  m_clusterMinLength(10.f),
24  m_halfWindowLayers(30),
25  m_samplingPitch(1.f),
26  m_maxCosSplittingAngle(0.9925f),
27  m_minCosMergingAngle(0.94f),
28  m_maxTransverseDisplacement(5.f),
29  m_maxLongitudinalDisplacement(30.f),
30  m_maxLongitudinalDisplacementSquared(m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)
31 {
32 }
33 
34 //------------------------------------------------------------------------------------------------------------------------------------------
35 
37 {
38  const ClusterList *pClusterList = NULL;
39  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pClusterList));
40 
41  // Get ordered list of clean clusters
42  ClusterVector clusterVector;
43  this->GetListOfCleanClusters(pClusterList, clusterVector);
44 
45  // Calculate sliding fit results for clean clusters
46  TwoDSlidingFitResultMap slidingFitResultMap;
47  this->BuildSlidingFitResultMap(clusterVector, slidingFitResultMap);
48 
49  // Loop over clusters, identify and perform splits
50  ClusterSet splitClusters;
51 
52  for (ClusterVector::const_iterator bIter = clusterVector.begin(), bIterEnd1 = clusterVector.end(); bIter != bIterEnd1; ++bIter)
53  {
54  if (splitClusters.count(*bIter) > 0)
55  continue;
56 
57  TwoDSlidingFitResultMap::const_iterator bFitIter = slidingFitResultMap.find(*bIter);
58 
59  if (slidingFitResultMap.end() == bFitIter)
60  continue;
61 
62  const TwoDSlidingFitResult &branchSlidingFitResult(bFitIter->second);
63 
64  // Find best split position for candidate branch cluster
65  CartesianVector splitPosition(0.f, 0.f, 0.f);
66  CartesianVector splitDirection1(0.f, 0.f, 0.f);
67  CartesianVector splitDirection2(0.f, 0.f, 0.f);
68 
69  if (STATUS_CODE_SUCCESS != this->FindBestSplitPosition(branchSlidingFitResult, splitPosition, splitDirection1, splitDirection2))
70  continue;
71 
72  // Find candidate replacement clusters to merge into branch cluster at the split position
73  TwoDSlidingFitResultMap::const_iterator bestReplacementIter1(slidingFitResultMap.end());
74  TwoDSlidingFitResultMap::const_iterator bestReplacementIter2(slidingFitResultMap.end());
75 
76  float bestLengthSquared1(m_maxLongitudinalDisplacementSquared);
77  float bestLengthSquared2(m_maxLongitudinalDisplacementSquared);
78 
79  for (ClusterVector::const_iterator rIter = clusterVector.begin(), rIterEnd = clusterVector.end(); rIter != rIterEnd; ++rIter)
80  {
81  if (splitClusters.count(*rIter) > 0)
82  continue;
83 
84  TwoDSlidingFitResultMap::const_iterator rFitIter = slidingFitResultMap.find(*rIter);
85 
86  if (slidingFitResultMap.end() == rFitIter)
87  continue;
88 
89  const TwoDSlidingFitResult &replacementSlidingFitResult(rFitIter->second);
90 
91  if (branchSlidingFitResult.GetCluster() == replacementSlidingFitResult.GetCluster())
92  continue;
93 
94  float lengthSquared1(std::numeric_limits<float>::max());
95  float lengthSquared2(std::numeric_limits<float>::max());
96 
97  if (STATUS_CODE_SUCCESS !=
98  this->ConfirmSplitPosition(branchSlidingFitResult, replacementSlidingFitResult, splitPosition, splitDirection1,
99  splitDirection2, lengthSquared1, lengthSquared2))
100  continue;
101 
102  if (lengthSquared1 < bestLengthSquared1)
103  {
104  bestLengthSquared1 = lengthSquared1;
105  bestReplacementIter1 = rFitIter;
106  }
107 
108  if (lengthSquared2 < bestLengthSquared2)
109  {
110  bestLengthSquared2 = lengthSquared2;
111  bestReplacementIter2 = rFitIter;
112  }
113  }
114 
115  const Cluster *pReplacementCluster1 = NULL;
116  const Cluster *pReplacementCluster2 = NULL;
117 
118  if (slidingFitResultMap.end() != bestReplacementIter1)
119  pReplacementCluster1 = bestReplacementIter1->first;
120 
121  if (slidingFitResultMap.end() != bestReplacementIter2)
122  pReplacementCluster2 = bestReplacementIter2->first;
123 
124  if (NULL == pReplacementCluster1 && NULL == pReplacementCluster2)
125  continue;
126 
127  const Cluster *const pBranchCluster = branchSlidingFitResult.GetCluster();
128 
129  // Crossed tracks
130  if (pReplacementCluster1 && pReplacementCluster2)
131  {
132  if (!(this->IdentifyCrossedTracks(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition)))
133  continue;
134 
135  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
136  this->PerformDoubleSplit(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition, splitDirection1, splitDirection2));
137  }
138  // Single scatter
139  else if (pReplacementCluster1)
140  {
141  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
142  this->PerformSingleSplit(pBranchCluster, pReplacementCluster1, splitPosition, splitDirection1, splitDirection2));
143  }
144  else if (pReplacementCluster2)
145  {
146  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
147  this->PerformSingleSplit(pBranchCluster, pReplacementCluster2, splitPosition, splitDirection2, splitDirection1));
148  }
149 
150  // Choose not to re-use clusters (for now)
151  if (pReplacementCluster1)
152  splitClusters.insert(pReplacementCluster1);
153 
154  if (pReplacementCluster2)
155  splitClusters.insert(pReplacementCluster2);
156 
157  splitClusters.insert(pBranchCluster);
158  }
159 
160  return STATUS_CODE_SUCCESS;
161 }
162 
163 //------------------------------------------------------------------------------------------------------------------------------------------
164 
165 void CosmicRaySplittingAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
166 {
167  for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
168  {
169  const Cluster *const pCluster = *iter;
170 
172  continue;
173 
174  clusterVector.push_back(pCluster);
175  }
176 
177  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
178 }
179 
180 //------------------------------------------------------------------------------------------------------------------------------------------
181 
183 {
184  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
185  {
186  if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
187  {
188  try
189  {
190  const float slidingFitPitch(LArGeometryHelper::GetWirePitch(this->GetPandora(), LArClusterHelper::GetClusterHitType(*iter)));
191  const TwoDSlidingFitResult slidingFitResult(*iter, m_halfWindowLayers, slidingFitPitch);
192 
193  if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
194  throw StatusCodeException(STATUS_CODE_FAILURE);
195  }
196  catch (StatusCodeException &statusCodeException)
197  {
198  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
199  throw statusCodeException;
200  }
201  }
202  }
203 }
204 
205 //------------------------------------------------------------------------------------------------------------------------------------------
206 
208  CartesianVector &splitPosition, CartesianVector &splitDirection1, CartesianVector &splitDirection2) const
209 {
210  // Find position of greatest scatter for this cluster
211  float splitCosTheta(m_maxCosSplittingAngle);
212  bool foundSplit(false);
213 
214  const CartesianVector &minPosition(branchSlidingFitResult.GetGlobalMinLayerPosition());
215  const CartesianVector &maxPosition(branchSlidingFitResult.GetGlobalMaxLayerPosition());
216  const float halfWindowLength(branchSlidingFitResult.GetLayerFitHalfWindowLength());
217 
218  float minL(0.f), maxL(0.f), minT(0.f), maxT(0.f);
219  branchSlidingFitResult.GetLocalPosition(minPosition, minL, minT);
220  branchSlidingFitResult.GetLocalPosition(maxPosition, maxL, maxT);
221 
222  const unsigned int nSamplingPoints = static_cast<unsigned int>((maxL - minL) / m_samplingPitch);
223 
224  for (unsigned int n = 0; n < nSamplingPoints; ++n)
225  {
226  const float alpha((0.5f + static_cast<float>(n)) / static_cast<float>(nSamplingPoints));
227  const float rL(minL + (maxL - minL) * alpha);
228 
229  CartesianVector centralPosition(0.f, 0.f, 0.f), forwardDirection(0.f, 0.f, 0.f), backwardDirection(0.f, 0.f, 0.f);
230 
231  if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL, centralPosition)) ||
232  (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL + halfWindowLength, forwardDirection)) ||
233  (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL - halfWindowLength, backwardDirection)))
234  {
235  continue;
236  }
237 
238  const float cosTheta(forwardDirection.GetDotProduct(backwardDirection));
239 
240  if (cosTheta < splitCosTheta)
241  {
242  CartesianVector forwardPosition(0.f, 0.f, 0.f);
243  CartesianVector backwardPosition(0.f, 0.f, 0.f);
244 
245  if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL + halfWindowLength, forwardPosition)) ||
246  (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL - halfWindowLength, backwardPosition)))
247  {
248  continue;
249  }
250 
251  CartesianVector forwardVectorOutwards(forwardPosition - centralPosition);
252  CartesianVector backwardVectorOutwards(backwardPosition - centralPosition);
253 
254  splitPosition = centralPosition;
255  splitDirection1 = (forwardDirection.GetDotProduct(forwardVectorOutwards) > 0.f) ? forwardDirection : forwardDirection * -1.f;
256  splitDirection2 = (backwardDirection.GetDotProduct(backwardVectorOutwards) > 0.f) ? backwardDirection : backwardDirection * -1.f;
257  splitCosTheta = cosTheta;
258  foundSplit = true;
259  }
260  }
261 
262  if (!foundSplit)
263  return STATUS_CODE_NOT_FOUND;
264 
265  return STATUS_CODE_SUCCESS;
266 }
267 
268 //------------------------------------------------------------------------------------------------------------------------------------------
269 
271  const TwoDSlidingFitResult &replacementSlidingFitResult, const CartesianVector &splitPosition, const CartesianVector &splitDirection1,
272  const CartesianVector &splitDirection2, float &bestLengthSquared1, float &bestLengthSquared2) const
273 {
274  // Check if the replacement cluster points to the split position on the branch cluster
275  bestLengthSquared1 = std::numeric_limits<float>::max();
276  bestLengthSquared2 = std::numeric_limits<float>::max();
277 
278  bool foundMatch(false);
279 
280  const CartesianVector minPosition(replacementSlidingFitResult.GetGlobalMinLayerPosition());
281  const CartesianVector maxPosition(replacementSlidingFitResult.GetGlobalMaxLayerPosition());
282  const CartesianVector minDirection(replacementSlidingFitResult.GetGlobalMinLayerDirection());
283  const CartesianVector maxDirection(replacementSlidingFitResult.GetGlobalMaxLayerDirection());
284 
285  const CartesianVector branchVertex1(branchSlidingFitResult.GetGlobalMinLayerPosition());
286  const CartesianVector branchVertex2(branchSlidingFitResult.GetGlobalMaxLayerPosition());
287  const CartesianVector branchDirection1(splitDirection1 * -1.f);
288  const CartesianVector branchDirection2(splitDirection2 * -1.f);
289 
290  const float cosSplittingAngle(-splitDirection1.GetDotProduct(splitDirection2));
291  const float branchLengthSquared((branchVertex2 - branchVertex1).GetMagnitudeSquared());
292  const float replacementLengthSquared((maxPosition - minPosition).GetMagnitudeSquared());
293 
294  // Loop over each end of the replacement cluster
295  for (unsigned int iFwd = 0; iFwd < 2; ++iFwd)
296  {
297  const CartesianVector pAxis((0 == iFwd) ? (maxPosition - minPosition) : (minPosition - maxPosition));
298  const CartesianVector vtxPosition((0 == iFwd) ? minPosition : maxPosition);
299  const CartesianVector endPosition((0 == iFwd) ? maxPosition : minPosition);
300  const CartesianVector vtxDirection((0 == iFwd) ? (pAxis.GetDotProduct(minDirection) > 0 ? minDirection : minDirection * -1.f)
301  : (pAxis.GetDotProduct(maxDirection) > 0 ? maxDirection : maxDirection * -1.f));
302 
303  // Choose the correct end of the replacement cluster and require proximity to the branch cluster
304  const float vtxDisplacement(LArClusterHelper::GetClosestDistance(vtxPosition, branchSlidingFitResult.GetCluster()));
305  const float endDisplacement(LArClusterHelper::GetClosestDistance(endPosition, branchSlidingFitResult.GetCluster()));
306 
307  const float lengthSquared((vtxPosition - splitPosition).GetMagnitudeSquared());
308  const float lengthSquared1((vtxPosition - branchVertex1).GetMagnitudeSquared());
309  const float lengthSquared2((vtxPosition - branchVertex2).GetMagnitudeSquared());
310 
311  if (vtxDisplacement > endDisplacement)
312  continue;
313 
314  if (std::min(lengthSquared, std::min(lengthSquared1, lengthSquared2)) > std::min(branchLengthSquared, replacementLengthSquared))
315  continue;
316 
317  // Require good pointing information between replacement cluster and candidate split position
318  float impactL(0.f), impactT(0.f);
319  LArPointingClusterHelper::GetImpactParameters(vtxPosition, vtxDirection, splitPosition, impactL, impactT);
320 
321  if (impactT > m_maxTransverseDisplacement || impactL > m_maxLongitudinalDisplacement || impactL < -1.f ||
322  impactT > std::max(1.5f, 0.577f * (1.f + impactL)))
323  continue;
324 
325  // Check the segment of the branch cluster above the split position
326  if (vtxDirection.GetDotProduct(branchDirection1) > std::max(m_minCosMergingAngle, cosSplittingAngle))
327  {
328  if (lengthSquared < bestLengthSquared1)
329  {
330  bestLengthSquared1 = lengthSquared;
331  foundMatch = true;
332  }
333  }
334 
335  // Check the segment of the branch cluster below the split position
336  if (vtxDirection.GetDotProduct(branchDirection2) > std::max(m_minCosMergingAngle, cosSplittingAngle))
337  {
338  if (lengthSquared < bestLengthSquared2)
339  {
340  bestLengthSquared2 = lengthSquared;
341  foundMatch = true;
342  }
343  }
344  }
345 
346  if (!foundMatch)
347  return STATUS_CODE_NOT_FOUND;
348 
349  return STATUS_CODE_SUCCESS;
350 }
351 
352 //------------------------------------------------------------------------------------------------------------------------------------------
353 
354 StatusCode CosmicRaySplittingAlgorithm::PerformSingleSplit(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster,
355  const CartesianVector &splitPosition, const CartesianVector &forwardDirection, const CartesianVector &backwardDirection) const
356 {
357  CaloHitList caloHitListToMove;
358  this->GetCaloHitListToMove(pBranchCluster, pReplacementCluster, splitPosition, forwardDirection, backwardDirection, caloHitListToMove);
359 
360  CaloHitList caloHitListToKeep;
361  this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove, caloHitListToKeep);
362 
363  if (caloHitListToKeep.empty())
364  return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster, pBranchCluster);
365 
366  return this->SplitCluster(pBranchCluster, pReplacementCluster, caloHitListToMove);
367 }
368 
369 //------------------------------------------------------------------------------------------------------------------------------------------
370 
371 StatusCode CosmicRaySplittingAlgorithm::PerformDoubleSplit(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster1,
372  const Cluster *const pReplacementCluster2, const CartesianVector &splitPosition, const CartesianVector &splitDirection1,
373  const CartesianVector &splitDirection2) const
374 {
375  CaloHitList caloHitListToMove1, caloHitListToMove2;
376  this->GetCaloHitListsToMove(pBranchCluster, splitPosition, splitDirection1, splitDirection2, caloHitListToMove1, caloHitListToMove2);
377 
378  CaloHitList caloHitListToKeep1;
379  this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove1, caloHitListToKeep1);
380 
381  if (caloHitListToKeep1.empty())
382  return STATUS_CODE_FAILURE;
383 
384  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->SplitCluster(pBranchCluster, pReplacementCluster1, caloHitListToMove1));
385 
386  CaloHitList caloHitListToKeep2;
387  this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove2, caloHitListToKeep2);
388 
389  if (caloHitListToKeep2.empty())
390  return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster2, pBranchCluster);
391 
392  return this->SplitCluster(pBranchCluster, pReplacementCluster2, caloHitListToMove2);
393 }
394 
395 //------------------------------------------------------------------------------------------------------------------------------------------
396 
397 void CosmicRaySplittingAlgorithm::GetCaloHitListToMove(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster,
398  const CartesianVector &splitPosition, const CartesianVector &forwardDirection, const CartesianVector &backwardDirection,
399  CaloHitList &caloHitListToMove) const
400 {
401  const CartesianVector forwardPosition(LArClusterHelper::GetClosestPosition(splitPosition, pBranchCluster));
402  const CartesianVector vtxPosition(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster));
403  const CartesianVector vtxDirection((forwardPosition - vtxPosition).GetUnitVector());
404 
405  CaloHitList caloHitListToSort, caloHitListToCheck;
406  pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
407 
408  for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)
409  {
410  const CaloHit *const pCaloHit = *iter;
411 
412  if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - forwardPosition) > 0.f)
413  {
414  caloHitListToMove.push_back(pCaloHit);
415  }
416  else if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - vtxPosition) > -1.25f)
417  {
418  caloHitListToCheck.push_back(pCaloHit);
419  }
420  }
421 
422  float closestLengthSquared(std::numeric_limits<float>::max());
423 
424  for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)
425  {
426  const CaloHit *const pCaloHit = *iter;
427 
428  if (vtxDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude() >
429  backwardDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude())
430  continue;
431 
432  const float lengthSquared((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared());
433 
434  if (lengthSquared < closestLengthSquared)
435  closestLengthSquared = lengthSquared;
436  }
437 
438  for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)
439  {
440  const CaloHit *const pCaloHit = *iter;
441 
442  if ((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared() >= closestLengthSquared)
443  caloHitListToMove.push_back(pCaloHit);
444  }
445 }
446 
447 //------------------------------------------------------------------------------------------------------------------------------------------
448 
449 void CosmicRaySplittingAlgorithm::GetCaloHitListsToMove(const Cluster *const pBranchCluster, const CartesianVector &splitPosition,
450  const CartesianVector &splitDirection1, const CartesianVector &splitDirection2, CaloHitList &caloHitListToMove1, CaloHitList &caloHitListToMove2) const
451 {
452  CaloHitList caloHitListToSort;
453  pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
454 
455  for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)
456  {
457  const CaloHit *const pCaloHit = *iter;
458 
459  const float displacement1(splitDirection1.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
460  const float displacement2(splitDirection2.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
461 
462  if (displacement1 > displacement2)
463  {
464  caloHitListToMove1.push_back(pCaloHit);
465  }
466  else
467  {
468  caloHitListToMove2.push_back(pCaloHit);
469  }
470  }
471 }
472 //------------------------------------------------------------------------------------------------------------------------------------------
473 
474 bool CosmicRaySplittingAlgorithm::IdentifyCrossedTracks(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster1,
475  const Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const
476 {
477  CartesianVector branchVertex1(0.f, 0.f, 0.f), branchVertex2(0.f, 0.f, 0.f);
478  LArClusterHelper::GetExtremalCoordinates(pBranchCluster, branchVertex1, branchVertex2);
479 
480  const CartesianVector replacementVertex1(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster1));
481  const CartesianVector replacementVertex2(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster2));
482 
483  if ((replacementVertex2 - replacementVertex1).GetMagnitudeSquared() > (branchVertex2 - branchVertex1).GetMagnitudeSquared())
484  return false;
485 
486  return true;
487 }
488 
489 //------------------------------------------------------------------------------------------------------------------------------------------
490 
492  const Cluster *const pBranchCluster, const CaloHitList &caloHitListToMove, CaloHitList &caloHitListToKeep) const
493 {
494  if (caloHitListToMove.empty())
495  return STATUS_CODE_FAILURE;
496 
497  CaloHitList caloHitList;
498  pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
499 
500  for (CaloHitList::const_iterator iter = caloHitList.begin(), iterEnd = caloHitList.end(); iter != iterEnd; ++iter)
501  {
502  const CaloHit *const pCaloHit = *iter;
503 
504  if (caloHitListToMove.end() == std::find(caloHitListToMove.begin(), caloHitListToMove.end(), pCaloHit))
505  caloHitListToKeep.push_back(pCaloHit);
506  }
507 
508  return STATUS_CODE_SUCCESS;
509 }
510 
511 //------------------------------------------------------------------------------------------------------------------------------------------
512 
514  const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster, const CaloHitList &caloHitListToMove) const
515 {
516  if (caloHitListToMove.empty())
517  return STATUS_CODE_FAILURE;
518 
519  for (CaloHitList::const_iterator iter = caloHitListToMove.begin(), iterEnd = caloHitListToMove.end(); iter != iterEnd; ++iter)
520  {
521  const CaloHit *const pCaloHit = *iter;
522  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::RemoveFromCluster(*this, pBranchCluster, pCaloHit));
523  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*this, pReplacementCluster, pCaloHit));
524  }
525 
526  return STATUS_CODE_SUCCESS;
527 }
528 
529 //------------------------------------------------------------------------------------------------------------------------------------------
530 
531 StatusCode CosmicRaySplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
532 {
533  PANDORA_RETURN_RESULT_IF_AND_IF(
534  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinLength", m_clusterMinLength));
535 
536  PANDORA_RETURN_RESULT_IF_AND_IF(
537  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitHalfWindow", m_halfWindowLayers));
538 
539  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SamplingPitch", m_samplingPitch));
540 
541  if (m_samplingPitch < std::numeric_limits<float>::epsilon())
542  return STATUS_CODE_INVALID_PARAMETER;
543 
544  PANDORA_RETURN_RESULT_IF_AND_IF(
545  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCosSplittingAngle", m_maxCosSplittingAngle));
546 
547  PANDORA_RETURN_RESULT_IF_AND_IF(
548  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosMergingAngle", m_minCosMergingAngle));
549 
550  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
551  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
552 
553  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
554  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
556 
557  return STATUS_CODE_SUCCESS;
558 }
559 
560 } // namespace lar_content
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position, then pulse-height.
void GetListOfCleanClusters(const pandora::ClusterList *const pClusterList, pandora::ClusterVector &clusterVector) const
Populate cluster vector with subset of cluster list, containing clusters judged to be clean...
void GetCaloHitListsToMove(const pandora::Cluster *const pBranchCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, pandora::CaloHitList &caloHitList1, pandora::CaloHitList &caloHitList2) const
Get lists of calo hits to move in order to split a branch cluster into two segments for case of two r...
bool IdentifyCrossedTracks(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const
Identify crossed tracks formed from the branch cluster and its replacement cluster.
Header file for the cosmic ray splitting algorithm class.
float m_maxLongitudinalDisplacement
maximum longitudinal displacement of associated clusters
void GetCaloHitListToMove(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection, pandora::CaloHitList &caloHitList) const
Get list of calo hits to move in order to split a branch cluster into two segments for case of one re...
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
intermediate_table::const_iterator const_iterator
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
float m_clusterMinLength
minimum length of clusters for this algorithm
TFile f
Definition: plotHisto.C:6
pandora::StatusCode SplitCluster(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CaloHitList &caloHitListToMove) const
Split the branch cluster and add hits to the replacement cluster.
Header file for the geometry helper class.
pandora::StatusCode PerformSingleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection) const
Split a branch cluster for case of one replacement cluster.
float m_minCosMergingAngle
largest relative angle allowed when merging clusters
if(nlines<=0)
Header file for the cluster helper class.
pandora::StatusCode GetCaloHitListToKeep(const pandora::Cluster *const pBranchCluster, const pandora::CaloHitList &caloHitListToMove, pandora::CaloHitList &caloHitListToKeep) const
Split the branch cluster and add hits to the replacement cluster.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
float m_maxCosSplittingAngle
smallest scatter angle allowed when splitting cluster
float m_maxLongitudinalDisplacementSquared
longitudinal displacement squared
pandora::StatusCode PerformDoubleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2) const
Split a branch cluster for case of two replacement clusters.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
static void GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z) ...
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
float m_maxTransverseDisplacement
maximum transverse displacement of associated clusters
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::StatusCode ConfirmSplitPosition(const TwoDSlidingFitResult &branchSlidingFitResult, const TwoDSlidingFitResult &replacementSlidingFitResult, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, float &lengthSquared1, float &lengthSquared2) const
Find a second replacement cluster that aligns with the scatter of the first branch cluster...
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
Char_t n[5]
std::vector< art::Ptr< recob::Cluster > > ClusterVector
pandora::StatusCode FindBestSplitPosition(const TwoDSlidingFitResult &slidingFitResult, pandora::CartesianVector &splitPosition, pandora::CartesianVector &splitDirection1, pandora::CartesianVector &splitDirection2) const
Find the position of greatest scatter along a sliding linear fit.
unsigned int m_halfWindowLayers
number of layers to use for half-window of sliding fit
static pandora::CartesianVector GetClosestPosition(const pandora::CartesianVector &position, const pandora::ClusterList &clusterList)
Get closest position in a list of clusters to a specified input position vector.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
float m_samplingPitch
sampling pitch for walking along sliding linear fit
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.