LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ThreeViewTrackFragmentsAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 ThreeViewTrackFragmentsAlgorithm::ThreeViewTrackFragmentsAlgorithm() :
22  m_nMaxTensorToolRepeats(1000),
23  m_minXOverlap(3.f),
24  m_minXOverlapFraction(0.8f),
25  m_maxPointDisplacementSquared(1.5f * 1.5f),
26  m_minMatchedSamplingPointFraction(0.5f),
27  m_minMatchedHits(5)
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
33 void ThreeViewTrackFragmentsAlgorithm::UpdateForNewCluster(const Cluster *const pNewCluster)
34 {
35  try
36  {
37  this->AddToSlidingFitCache(pNewCluster);
38  }
39  catch (StatusCodeException &statusCodeException)
40  {
41  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
42  throw statusCodeException;
43 
44  return;
45  }
46 
47  const HitType hitType(LArClusterHelper::GetClusterHitType(pNewCluster));
48 
49  if (!((TPC_VIEW_U == hitType) || (TPC_VIEW_V == hitType) || (TPC_VIEW_W == hitType)))
50  throw StatusCodeException(STATUS_CODE_FAILURE);
51 
52  // ATTN This is non-standard usage, supported here only (for legacy purposes)
53  MatchingType &matchingControl(this->GetMatchingControl());
54  ClusterList &clusterList((TPC_VIEW_U == hitType) ? matchingControl.m_clusterListU
55  : (TPC_VIEW_V == hitType) ? matchingControl.m_clusterListV
56  : matchingControl.m_clusterListW);
57 
58  if (clusterList.end() != std::find(clusterList.begin(), clusterList.end(), pNewCluster))
59  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
60 
61  clusterList.push_back(pNewCluster);
62 
63  ClusterList clusterList1(this->GetSelectedClusterList((TPC_VIEW_U == hitType) ? TPC_VIEW_V : TPC_VIEW_U));
64  ClusterList clusterList2(this->GetSelectedClusterList((TPC_VIEW_W == hitType) ? TPC_VIEW_V : TPC_VIEW_W));
65  clusterList1.sort(LArClusterHelper::SortByNHits);
66  clusterList2.sort(LArClusterHelper::SortByNHits);
67 
68  for (const Cluster *const pCluster1 : clusterList1)
69  {
70  if (TPC_VIEW_U == hitType)
71  {
72  this->CalculateOverlapResult(pNewCluster, pCluster1, nullptr);
73  }
74  else if (TPC_VIEW_V == hitType)
75  {
76  this->CalculateOverlapResult(pCluster1, pNewCluster, nullptr);
77  }
78  else
79  {
80  this->CalculateOverlapResult(pCluster1, nullptr, pNewCluster);
81  }
82  }
83 
84  for (const Cluster *const pCluster2 : clusterList2)
85  {
86  if (TPC_VIEW_U == hitType)
87  {
88  this->CalculateOverlapResult(pNewCluster, nullptr, pCluster2);
89  }
90  else if (TPC_VIEW_V == hitType)
91  {
92  this->CalculateOverlapResult(nullptr, pNewCluster, pCluster2);
93  }
94  else
95  {
96  this->CalculateOverlapResult(nullptr, pCluster2, pNewCluster);
97  }
98  }
99 }
100 
101 //------------------------------------------------------------------------------------------------------------------------------------------
102 
103 void ThreeViewTrackFragmentsAlgorithm::RebuildClusters(const ClusterList &rebuildList, ClusterList &newClusters) const
104 {
105  const ClusterList *pNewClusterList = nullptr;
106  std::string oldClusterListName, newClusterListName;
107 
108  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::InitializeReclustering(*this, TrackList(), rebuildList, oldClusterListName));
109  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
110  PandoraContentApi::RunClusteringAlgorithm(*this, m_reclusteringAlgorithmName, pNewClusterList, newClusterListName));
111 
112  newClusters.insert(newClusters.end(), pNewClusterList->begin(), pNewClusterList->end());
113  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::EndReclustering(*this, newClusterListName));
114 }
115 
116 //------------------------------------------------------------------------------------------------------------------------------------------
117 
119 {
120  ClusterList clusterListU(this->GetSelectedClusterList(TPC_VIEW_U));
121  ClusterList clusterListV(this->GetSelectedClusterList(TPC_VIEW_V));
122  ClusterList clusterListW(this->GetSelectedClusterList(TPC_VIEW_W));
123  clusterListU.sort(LArClusterHelper::SortByNHits);
124  clusterListV.sort(LArClusterHelper::SortByNHits);
125  clusterListW.sort(LArClusterHelper::SortByNHits);
126 
127  for (const Cluster *const pClusterU : clusterListU)
128  {
129  for (const Cluster *const pClusterV : clusterListV)
130  this->CalculateOverlapResult(pClusterU, pClusterV, nullptr);
131  }
132 
133  for (const Cluster *const pClusterU : clusterListU)
134  {
135  for (const Cluster *const pClusterW : clusterListW)
136  this->CalculateOverlapResult(pClusterU, nullptr, pClusterW);
137  }
138 
139  for (const Cluster *const pClusterV : clusterListV)
140  {
141  for (const Cluster *const pClusterW : clusterListW)
142  this->CalculateOverlapResult(nullptr, pClusterV, pClusterW);
143  }
144 }
145 
146 //------------------------------------------------------------------------------------------------------------------------------------------
147 
148 void ThreeViewTrackFragmentsAlgorithm::CalculateOverlapResult(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW)
149 {
150  const HitType missingHitType(((nullptr != pClusterU) && (nullptr != pClusterV) && (nullptr == pClusterW)) ? TPC_VIEW_W
151  : ((nullptr != pClusterU) && (nullptr == pClusterV) && (nullptr != pClusterW)) ? TPC_VIEW_V
152  : ((nullptr == pClusterU) && (nullptr != pClusterV) && (nullptr != pClusterW)) ? TPC_VIEW_U
153  : HIT_CUSTOM);
154 
155  if (HIT_CUSTOM == missingHitType)
156  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
157 
158  // Calculate new overlap result and replace old overlap result where necessary
159  FragmentOverlapResult oldOverlapResult, newOverlapResult;
160  const Cluster *pMatchedClusterU(nullptr), *pMatchedClusterV(nullptr), *pMatchedClusterW(nullptr);
161 
162  const TwoDSlidingFitResult &fitResult1(
163  (TPC_VIEW_U == missingHitType) ? this->GetCachedSlidingFitResult(pClusterV) : this->GetCachedSlidingFitResult(pClusterU));
164 
165  const TwoDSlidingFitResult &fitResult2((TPC_VIEW_U == missingHitType) ? this->GetCachedSlidingFitResult(pClusterW)
166  : (TPC_VIEW_V == missingHitType) ? this->GetCachedSlidingFitResult(pClusterW)
167  : this->GetCachedSlidingFitResult(pClusterV));
168 
169  const ClusterList &inputClusterList(this->GetInputClusterList(missingHitType));
170 
171  const Cluster *pBestMatchedCluster(nullptr);
172  const StatusCode statusCode(this->CalculateOverlapResult(fitResult1, fitResult2, inputClusterList, pBestMatchedCluster, newOverlapResult));
173 
174  if ((STATUS_CODE_SUCCESS != statusCode) && (STATUS_CODE_NOT_FOUND != statusCode))
175  throw StatusCodeException(statusCode);
176 
177  if (!newOverlapResult.IsInitialized())
178  return;
179 
180  MatchingType::TensorType &overlapTensor(this->GetMatchingControl().GetOverlapTensor());
181 
182  if (STATUS_CODE_SUCCESS == statusCode)
183  {
184  pMatchedClusterU = ((nullptr != pClusterU) ? pClusterU : pBestMatchedCluster);
185  pMatchedClusterV = ((nullptr != pClusterV) ? pClusterV : pBestMatchedCluster);
186  pMatchedClusterW = ((nullptr != pClusterW) ? pClusterW : pBestMatchedCluster);
187 
188  if (nullptr == pMatchedClusterU || nullptr == pMatchedClusterV || nullptr == pMatchedClusterW)
189  throw StatusCodeException(STATUS_CODE_FAILURE);
190 
191  try
192  {
193  oldOverlapResult = overlapTensor.GetOverlapResult(pMatchedClusterU, pMatchedClusterV, pMatchedClusterW);
194  }
195  catch (StatusCodeException &)
196  {
197  }
198  }
199 
200  if (!oldOverlapResult.IsInitialized())
201  {
202  overlapTensor.SetOverlapResult(pMatchedClusterU, pMatchedClusterV, pMatchedClusterW, newOverlapResult);
203  }
204  else if (newOverlapResult.GetFragmentCaloHitList().size() > oldOverlapResult.GetFragmentCaloHitList().size())
205  {
206  overlapTensor.ReplaceOverlapResult(pMatchedClusterU, pMatchedClusterV, pMatchedClusterW, newOverlapResult);
207  }
208  else if (newOverlapResult.GetFragmentCaloHitList().size() == oldOverlapResult.GetFragmentCaloHitList().size())
209  {
210  float newEnergySum(0.f), oldEnergySum(0.f);
211  for (const CaloHit *const pCaloHit : newOverlapResult.GetFragmentCaloHitList())
212  newEnergySum += pCaloHit->GetHadronicEnergy();
213  for (const CaloHit *const pCaloHit : oldOverlapResult.GetFragmentCaloHitList())
214  oldEnergySum += pCaloHit->GetHadronicEnergy();
215 
216  if (newEnergySum > oldEnergySum)
217  overlapTensor.ReplaceOverlapResult(pMatchedClusterU, pMatchedClusterV, pMatchedClusterW, newOverlapResult);
218  }
219 }
220 
221 //------------------------------------------------------------------------------------------------------------------------------------------
222 
224  const ClusterList &inputClusterList, const Cluster *&pBestMatchedCluster, FragmentOverlapResult &fragmentOverlapResult) const
225 {
226  const Cluster *const pCluster1(fitResult1.GetCluster());
227  const Cluster *const pCluster2(fitResult2.GetCluster());
228 
229  float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
230  pCluster1->GetClusterSpanX(xMin1, xMax1);
231  pCluster2->GetClusterSpanX(xMin2, xMax2);
232 
233  const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
234 
235  if (xOverlap < std::numeric_limits<float>::epsilon())
236  return STATUS_CODE_NOT_FOUND;
237 
238  CartesianPointVector projectedPositions;
239  const StatusCode statusCode1(this->GetProjectedPositions(fitResult1, fitResult2, projectedPositions));
240 
241  if (STATUS_CODE_SUCCESS != statusCode1)
242  return statusCode1;
243 
244  CaloHitList matchedHits;
245  ClusterList matchedClusters;
246  HitToClusterMap hitToClusterMap;
247  const StatusCode statusCode2(this->GetMatchedHits(inputClusterList, projectedPositions, hitToClusterMap, matchedHits));
248 
249  if (STATUS_CODE_SUCCESS != statusCode2)
250  return statusCode2;
251 
252  const StatusCode statusCode3(this->GetMatchedClusters(matchedHits, hitToClusterMap, matchedClusters, pBestMatchedCluster));
253 
254  if (STATUS_CODE_SUCCESS != statusCode3)
255  return statusCode3;
256 
257  if (!this->CheckMatchedClusters(projectedPositions, matchedClusters))
258  return STATUS_CODE_NOT_FOUND;
259 
260  FragmentOverlapResult overlapResult;
261  this->GetFragmentOverlapResult(projectedPositions, matchedHits, matchedClusters, overlapResult);
262 
263  if (!this->CheckOverlapResult(overlapResult))
264  return STATUS_CODE_NOT_FOUND;
265 
266  fragmentOverlapResult = overlapResult;
267  return STATUS_CODE_SUCCESS;
268 }
269 
270 //------------------------------------------------------------------------------------------------------------------------------------------
271 
273  const TwoDSlidingFitResult &fitResult1, const TwoDSlidingFitResult &fitResult2, CartesianPointVector &projectedPositions) const
274 {
275  const Cluster *const pCluster1(fitResult1.GetCluster());
276  const Cluster *const pCluster2(fitResult2.GetCluster());
277 
278  // Check hit types
279  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
280  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
281  const HitType hitType3((TPC_VIEW_U != hitType1 && TPC_VIEW_U != hitType2) ? TPC_VIEW_U
282  : (TPC_VIEW_V != hitType1 && TPC_VIEW_V != hitType2) ? TPC_VIEW_V
283  : (TPC_VIEW_W != hitType1 && TPC_VIEW_W != hitType2) ? TPC_VIEW_W
284  : HIT_CUSTOM);
285 
286  if (HIT_CUSTOM == hitType3)
287  return STATUS_CODE_INVALID_PARAMETER;
288 
289  // Check absolute and fractional overlap in x coordinate
290  float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
291  pCluster1->GetClusterSpanX(xMin1, xMax1);
292  pCluster2->GetClusterSpanX(xMin2, xMax2);
293 
294  const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
295  const float xSpan(std::max(xMax1, xMax2) - std::min(xMin1, xMin2));
296 
297  if ((xOverlap < m_minXOverlap) || (xSpan < std::numeric_limits<float>::epsilon()) || ((xOverlap / xSpan) < m_minXOverlapFraction))
298  return STATUS_CODE_NOT_FOUND;
299 
300  // Identify vertex and end positions (2D)
301  const CartesianVector minPosition1(fitResult1.GetGlobalMinLayerPosition());
302  const CartesianVector maxPosition1(fitResult1.GetGlobalMaxLayerPosition());
303  const CartesianVector minPosition2(fitResult2.GetGlobalMinLayerPosition());
304  const CartesianVector maxPosition2(fitResult2.GetGlobalMaxLayerPosition());
305 
306  const float dx_A(std::fabs(minPosition2.GetX() - minPosition1.GetX()));
307  const float dx_B(std::fabs(maxPosition2.GetX() - maxPosition1.GetX()));
308  const float dx_C(std::fabs(maxPosition2.GetX() - minPosition1.GetX()));
309  const float dx_D(std::fabs(minPosition2.GetX() - maxPosition1.GetX()));
310 
311  if (std::min(dx_C, dx_D) > std::max(dx_A, dx_B) && std::min(dx_A, dx_B) > std::max(dx_C, dx_D))
312  return STATUS_CODE_NOT_FOUND;
313 
314  const CartesianVector &vtxPosition1(minPosition1);
315  const CartesianVector &endPosition1(maxPosition1);
316  const CartesianVector &vtxPosition2((dx_A < dx_C) ? minPosition2 : maxPosition2);
317  const CartesianVector &endPosition2((dx_A < dx_C) ? maxPosition2 : minPosition2);
318 
319  // Calculate vertex and end positions (3D)
320  float vtxChi2(0.f);
321  CartesianVector vtxPosition3D(0.f, 0.f, 0.f);
322  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), hitType1, hitType2, vtxPosition1, vtxPosition2, vtxPosition3D, vtxChi2);
323 
324  float endChi2(0.f);
325  CartesianVector endPosition3D(0.f, 0.f, 0.f);
326  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), hitType1, hitType2, endPosition1, endPosition2, endPosition3D, endChi2);
327 
328  const CartesianVector vtxProjection3(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtxPosition3D, hitType3));
329  const CartesianVector endProjection3(LArGeometryHelper::ProjectPosition(this->GetPandora(), endPosition3D, hitType3));
330 
331  const float pitchU{LArGeometryHelper::GetWirePitch(this->GetPandora(), TPC_VIEW_U)};
332  const float pitchV{LArGeometryHelper::GetWirePitch(this->GetPandora(), TPC_VIEW_V)};
333  const float pitchW{LArGeometryHelper::GetWirePitch(this->GetPandora(), TPC_VIEW_W)};
334  const float pitchMax{std::max({pitchU, pitchV, pitchW})};
335  const float samplingPitch(0.5f * pitchMax);
336  const float nSamplingPoints((endProjection3 - vtxProjection3).GetMagnitude() / samplingPitch);
337 
338  if (nSamplingPoints < 1.f)
339  return STATUS_CODE_NOT_FOUND;
340 
341  // Loop over trajectory points
342  bool foundLastPosition(false);
343  CartesianVector lastPosition(0.f, 0.f, 0.f);
344 
345  for (float iSample = 0.5f; iSample < nSamplingPoints; iSample += 1.f)
346  {
347  const CartesianVector linearPosition3D(vtxPosition3D + (endPosition3D - vtxPosition3D) * (iSample / nSamplingPoints));
348  const CartesianVector linearPosition1(LArGeometryHelper::ProjectPosition(this->GetPandora(), linearPosition3D, hitType1));
349  const CartesianVector linearPosition2(LArGeometryHelper::ProjectPosition(this->GetPandora(), linearPosition3D, hitType2));
350 
351  float chi2(0.f);
352  CartesianVector fitPosition1(0.f, 0.f, 0.f), fitPosition2(0.f, 0.f, 0.f);
353 
354  if ((STATUS_CODE_SUCCESS != fitResult1.GetGlobalFitProjection(linearPosition1, fitPosition1)) ||
355  (STATUS_CODE_SUCCESS != fitResult2.GetGlobalFitProjection(linearPosition2, fitPosition2)))
356  {
357  continue;
358  }
359 
360  float rL1(0.f), rL2(0.f), rT1(0.f), rT2(0.f);
361  fitResult1.GetLocalPosition(fitPosition1, rL1, rT1);
362  fitResult2.GetLocalPosition(fitPosition2, rL2, rT2);
363 
364  const float x(0.5 * (fitPosition1.GetX() + fitPosition2.GetX()));
365  CartesianVector position1(0.f, 0.f, 0.f), position2(0.f, 0.f, 0.f), position3(0.f, 0.f, 0.f);
366 
367  try
368  {
369  const FitSegment &fitSegment1 = fitResult1.GetFitSegment(rL1);
370  const FitSegment &fitSegment2 = fitResult2.GetFitSegment(rL2);
371 
372  if ((STATUS_CODE_SUCCESS != fitResult1.GetTransverseProjection(x, fitSegment1, position1)) ||
373  (STATUS_CODE_SUCCESS != fitResult2.GetTransverseProjection(x, fitSegment2, position2)))
374  {
375  continue;
376  }
377  }
378  catch (StatusCodeException &statusCodeException)
379  {
380  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
381  throw statusCodeException;
382 
383  continue;
384  }
385 
386  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType1, hitType2, position1, position2, position3, chi2);
387 
388  // TODO For highly multi-valued x, projected positions can be unreliable. Need to make interpolation more robust for these cases.
389  if (foundLastPosition)
390  {
391  const float thisDisplacement((lastPosition - position3).GetMagnitude());
392  if (thisDisplacement > 2.f * samplingPitch)
393  {
394  const float nExtraPoints(thisDisplacement / samplingPitch);
395  for (float iExtra = 0.5f; iExtra < nExtraPoints; iExtra += 1.f)
396  {
397  const CartesianVector extraPosition(position3 + (lastPosition - position3) * (iExtra / nExtraPoints));
398  projectedPositions.push_back(extraPosition);
399  }
400  }
401  }
402 
403  projectedPositions.push_back(position3);
404 
405  lastPosition = position3;
406  foundLastPosition = true;
407  }
408 
409  // Bail out if list of projected positions is empty
410  if (projectedPositions.empty())
411  return STATUS_CODE_NOT_FOUND;
412 
413  return STATUS_CODE_SUCCESS;
414 }
415 
416 //------------------------------------------------------------------------------------------------------------------------------------------
417 
418 StatusCode ThreeViewTrackFragmentsAlgorithm::GetMatchedHits(const ClusterList &inputClusterList,
419  const CartesianPointVector &projectedPositions, HitToClusterMap &hitToClusterMap, CaloHitList &matchedHits) const
420 {
421  CaloHitVector availableCaloHits;
422 
423  for (ClusterList::const_iterator iter = inputClusterList.begin(), iterEnd = inputClusterList.end(); iter != iterEnd; ++iter)
424  {
425  const Cluster *const pCluster = *iter;
426 
427  if (!pCluster->IsAvailable())
428  continue;
429 
430  CaloHitList caloHitList;
431  pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
432  availableCaloHits.insert(availableCaloHits.end(), caloHitList.begin(), caloHitList.end());
433 
434  for (CaloHitList::const_iterator hIter = caloHitList.begin(), hIterEnd = caloHitList.end(); hIter != hIterEnd; ++hIter)
435  hitToClusterMap.insert(HitToClusterMap::value_type(*hIter, pCluster));
436  }
437 
438  std::sort(availableCaloHits.begin(), availableCaloHits.end(), LArClusterHelper::SortHitsByPosition);
439 
440  for (const CartesianVector &projectedPosition : projectedPositions)
441  {
442  const CaloHit *pClosestCaloHit(nullptr);
443  float closestDistanceSquared(std::numeric_limits<float>::max()), tieBreakerBestEnergy(0.f);
444 
445  for (const CaloHit *const pCaloHit : availableCaloHits)
446  {
447  const float distanceSquared((pCaloHit->GetPositionVector() - projectedPosition).GetMagnitudeSquared());
448 
449  if ((distanceSquared < closestDistanceSquared) ||
450  ((std::fabs(distanceSquared - closestDistanceSquared) < std::numeric_limits<float>::epsilon()) &&
451  (pCaloHit->GetHadronicEnergy() > tieBreakerBestEnergy)))
452  {
453  pClosestCaloHit = pCaloHit;
454  closestDistanceSquared = distanceSquared;
455  tieBreakerBestEnergy = pCaloHit->GetHadronicEnergy();
456  }
457  }
458 
459  if ((closestDistanceSquared < m_maxPointDisplacementSquared) && (nullptr != pClosestCaloHit) &&
460  (matchedHits.end() == std::find(matchedHits.begin(), matchedHits.end(), pClosestCaloHit)))
461  matchedHits.push_back(pClosestCaloHit);
462  }
463 
464  if (matchedHits.empty())
465  return STATUS_CODE_NOT_FOUND;
466 
467  return STATUS_CODE_SUCCESS;
468 }
469 
470 //------------------------------------------------------------------------------------------------------------------------------------------
471 
472 StatusCode ThreeViewTrackFragmentsAlgorithm::GetMatchedClusters(const CaloHitList &matchedHits, const HitToClusterMap &hitToClusterMap,
473  ClusterList &matchedClusters, const Cluster *&pBestMatchedCluster) const
474 {
475  ClusterToMatchedHitsMap clusterToMatchedHitsMap;
476 
477  for (CaloHitList::const_iterator iter = matchedHits.begin(), iterEnd = matchedHits.end(); iter != iterEnd; ++iter)
478  {
479  const CaloHit *const pCaloHit = *iter;
480  HitToClusterMap::const_iterator cIter = hitToClusterMap.find(pCaloHit);
481 
482  if (hitToClusterMap.end() == cIter)
483  throw StatusCodeException(STATUS_CODE_FAILURE);
484 
485  ++clusterToMatchedHitsMap[cIter->second];
486 
487  if (matchedClusters.end() == std::find(matchedClusters.begin(), matchedClusters.end(), cIter->second))
488  matchedClusters.push_back(cIter->second);
489  }
490 
491  if (matchedClusters.empty())
492  return STATUS_CODE_NOT_FOUND;
493 
494  pBestMatchedCluster = nullptr;
495  unsigned int bestClusterMatchedHits(0);
496  float tieBreakerBestEnergy(0.f);
497 
498  ClusterVector sortedClusters;
499  for (const auto &mapEntry : clusterToMatchedHitsMap)
500  sortedClusters.push_back(mapEntry.first);
501  std::sort(sortedClusters.begin(), sortedClusters.end(), LArClusterHelper::SortByNHits);
502 
503  for (const Cluster *const pCluster : sortedClusters)
504  {
505  const unsigned int nMatchedHits(clusterToMatchedHitsMap.at(pCluster));
506 
507  if ((nMatchedHits > bestClusterMatchedHits) || ((nMatchedHits == bestClusterMatchedHits) && (pCluster->GetHadronicEnergy() > tieBreakerBestEnergy)))
508  {
509  pBestMatchedCluster = pCluster;
510  bestClusterMatchedHits = nMatchedHits;
511  tieBreakerBestEnergy = pCluster->GetHadronicEnergy();
512  }
513  }
514 
515  if (nullptr == pBestMatchedCluster)
516  return STATUS_CODE_NOT_FOUND;
517 
518  return STATUS_CODE_SUCCESS;
519 }
520 
521 //------------------------------------------------------------------------------------------------------------------------------------------
522 
523 void ThreeViewTrackFragmentsAlgorithm::GetFragmentOverlapResult(const CartesianPointVector &projectedPositions,
524  const CaloHitList &matchedHits, const ClusterList &matchedClusters, FragmentOverlapResult &fragmentOverlapResult) const
525 {
526  float chi2Sum(0.f);
527  unsigned int nMatchedSamplingPoints(0);
528 
529  CaloHitVector sortedMatchedHits(matchedHits.begin(), matchedHits.end());
530  std::sort(sortedMatchedHits.begin(), sortedMatchedHits.end(), LArClusterHelper::SortHitsByPosition);
531 
532  for (const CartesianVector &projectedPosition : projectedPositions)
533  {
534  float closestDistanceSquared(std::numeric_limits<float>::max());
535 
536  for (const CaloHit *const pCaloHit : matchedHits)
537  {
538  const float distanceSquared((pCaloHit->GetPositionVector() - projectedPosition).GetMagnitudeSquared());
539 
540  if (distanceSquared < closestDistanceSquared)
541  closestDistanceSquared = distanceSquared;
542  }
543 
544  if (closestDistanceSquared < m_maxPointDisplacementSquared)
545  {
546  ++nMatchedSamplingPoints;
547  chi2Sum += closestDistanceSquared;
548  }
549  }
550 
551  fragmentOverlapResult = FragmentOverlapResult(nMatchedSamplingPoints, projectedPositions.size(), chi2Sum, matchedHits, matchedClusters);
552 }
553 
554 //------------------------------------------------------------------------------------------------------------------------------------------
555 
556 bool ThreeViewTrackFragmentsAlgorithm::CheckMatchedClusters(const CartesianPointVector &projectedPositions, const ClusterList &matchedClusters) const
557 {
558  if (projectedPositions.empty() || matchedClusters.empty())
559  return false;
560 
561  // Calculate X and Z span of projected positions
562  float minXproj(+std::numeric_limits<float>::max());
563  float maxXproj(-std::numeric_limits<float>::max());
564  float minZproj(+std::numeric_limits<float>::max());
565  float maxZproj(-std::numeric_limits<float>::max());
566 
567  for (const CartesianVector &projectedPosition : projectedPositions)
568  {
569  minXproj = std::min(minXproj, projectedPosition.GetX());
570  maxXproj = std::max(maxXproj, projectedPosition.GetX());
571  minZproj = std::min(minZproj, projectedPosition.GetZ());
572  maxZproj = std::max(maxZproj, projectedPosition.GetZ());
573  }
574 
575  const float dXproj(maxXproj - minXproj);
576  const float dZproj(maxZproj - minZproj);
577  const float projectedLengthSquared(dXproj * dXproj + dZproj * dZproj);
578 
579  // Calculate X and Z span of matched clusters
580  float minXcluster(+std::numeric_limits<float>::max());
581  float maxXcluster(-std::numeric_limits<float>::max());
582  float minZcluster(+std::numeric_limits<float>::max());
583  float maxZcluster(-std::numeric_limits<float>::max());
584 
585  for (const Cluster *const pCluster : matchedClusters)
586  {
587  CartesianVector minPosition(0.f, 0.f, 0.f);
588  CartesianVector maxPosition(0.f, 0.f, 0.f);
589 
590  LArClusterHelper::GetClusterBoundingBox(pCluster, minPosition, maxPosition);
591 
592  minXcluster = std::min(minXcluster, minPosition.GetX());
593  maxXcluster = std::max(maxXcluster, maxPosition.GetX());
594  minZcluster = std::min(minZcluster, minPosition.GetZ());
595  maxZcluster = std::max(maxZcluster, maxPosition.GetZ());
596  }
597 
598  const float dXcluster(maxXcluster - minXcluster);
599  const float dZcluster(maxZcluster - minZcluster);
600  const float clusterLengthSquared(dXcluster * dXcluster + dZcluster * dZcluster);
601 
602  // Require that the span of the matched clusters is no larger than twice the span of the projected positions
603  if (clusterLengthSquared > 4.f * projectedLengthSquared)
604  return false;
605 
606  return true;
607 }
608 
609 //------------------------------------------------------------------------------------------------------------------------------------------
610 
612 {
613  // ATTN This method is currently mirrored in ClearTrackFragments tool
615  return false;
616 
617  if (overlapResult.GetFragmentCaloHitList().size() < m_minMatchedHits)
618  return false;
619 
620  return true;
621 }
622 
623 //------------------------------------------------------------------------------------------------------------------------------------------
624 
626 {
627  unsigned int repeatCounter(0);
628 
629  for (TensorToolVector::const_iterator iter = m_algorithmToolVector.begin(), iterEnd = m_algorithmToolVector.end(); iter != iterEnd;)
630  {
631  if ((*iter)->Run(this, this->GetMatchingControl().GetOverlapTensor()))
632  {
633  iter = m_algorithmToolVector.begin();
634 
635  if (++repeatCounter > m_nMaxTensorToolRepeats)
636  break;
637  }
638  else
639  {
640  ++iter;
641  }
642  }
643 }
644 
645 //------------------------------------------------------------------------------------------------------------------------------------------
646 
647 StatusCode ThreeViewTrackFragmentsAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
648 {
649  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithm(*this, xmlHandle, "ClusterRebuilding", m_reclusteringAlgorithmName));
650 
651  AlgorithmToolVector algorithmToolVector;
652  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*this, xmlHandle, "TrackTools", algorithmToolVector));
653 
654  for (AlgorithmToolVector::const_iterator iter = algorithmToolVector.begin(), iterEnd = algorithmToolVector.end(); iter != iterEnd; ++iter)
655  {
656  FragmentTensorTool *const pFragmentTensorTool(dynamic_cast<FragmentTensorTool *>(*iter));
657 
658  if (!pFragmentTensorTool)
659  return STATUS_CODE_INVALID_PARAMETER;
660 
661  m_algorithmToolVector.push_back(pFragmentTensorTool);
662  }
663 
664  PANDORA_RETURN_RESULT_IF_AND_IF(
665  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NMaxTensorToolRepeats", m_nMaxTensorToolRepeats));
666 
667  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlap", m_minXOverlap));
668 
669  PANDORA_RETURN_RESULT_IF_AND_IF(
670  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlapFraction", m_minXOverlapFraction));
671 
672  float maxPointDisplacement = std::sqrt(m_maxPointDisplacementSquared);
673  PANDORA_RETURN_RESULT_IF_AND_IF(
674  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxPointDisplacement", maxPointDisplacement));
675  m_maxPointDisplacementSquared = maxPointDisplacement * maxPointDisplacement;
676 
677  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
678  XmlHelper::ReadValue(xmlHandle, "MinMatchedSamplingPointFraction", m_minMatchedSamplingPointFraction));
679 
680  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinMatchedHits", m_minMatchedHits));
681 
682  return BaseAlgorithm::ReadSettings(xmlHandle);
683 }
684 
685 } // namespace lar_content
Float_t x
Definition: compare.C:6
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::ClusterList m_clusterListW
The selected modified cluster list W.
pandora::ClusterList m_clusterListU
The selected modified cluster list U.
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.
FragmentOverlapResult class.
void AddToSlidingFitCache(const pandora::Cluster *const pCluster)
Add a new sliding fit result, for the specified cluster, to the algorithm cache.
void RebuildClusters(const pandora::ClusterList &rebuildList, pandora::ClusterList &newClusters) const
Rebuild clusters after fragmentation.
void GetFragmentOverlapResult(const pandora::CartesianPointVector &projectedPositions, const pandora::CaloHitList &matchedHits, const pandora::ClusterList &matchedClusters, FragmentOverlapResult &fragmentOverlapResult) const
Get the populated fragment overlap result.
static void MergeTwoPositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from two views to give unified 3D position.
pandora::ClusterList m_clusterListV
The selected modified cluster list V.
void PerformMainLoop()
Main loop over cluster combinations in order to populate the overlap container. Responsible for calli...
bool IsInitialized() const
Whether the track overlap result has been initialized.
pandora::StatusCode GetMatchedClusters(const pandora::CaloHitList &matchedHits, const HitToClusterMap &hitToClusterMap, pandora::ClusterList &matchedClusters, const pandora::Cluster *&pBestMatchedCluster) const
Get the list of the relevant clusters and the address of the single best matched cluster.
const TwoDSlidingFitResult & GetCachedSlidingFitResult(const pandora::Cluster *const pCluster) const
Get a sliding fit result from the algorithm cache.
pandora::StatusCode GetProjectedPositions(const TwoDSlidingFitResult &fitResult1, const TwoDSlidingFitResult &fitResult2, pandora::CartesianPointVector &projectedPositions) const
Get the list of projected positions, in the third view, corresponding to a pair of sliding fit result...
std::unordered_map< const pandora::Cluster *, unsigned int > ClusterToMatchedHitsMap
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
unsigned int m_minMatchedHits
minimum number of matched calo hits
intermediate_table::const_iterator const_iterator
void ReplaceOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW, const OverlapResult &overlapResult)
SetReplace an existing overlap result.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
void ExamineOverlapContainer()
Examine contents of overlap container, collect together best-matching 2D particles and modify cluster...
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void UpdateForNewCluster(const pandora::Cluster *const pNewCluster)
Update to reflect addition of a new cluster to the problem space.
static void GetClusterBoundingBox(const pandora::Cluster *const pCluster, pandora::CartesianVector &minimumCoordinate, pandora::CartesianVector &maximumCoordinate)
Get minimum and maximum X, Y and Z positions of the calo hits in a cluster.
const pandora::CaloHitList & GetFragmentCaloHitList() const
Get the list of fragment-associated hits.
void SetOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW, const OverlapResult &overlapResult)
Set overlap result.
bool CheckOverlapResult(const FragmentOverlapResult &overlapResult) const
Whether the matched clusters and hits pass the algorithm quality cuts.
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.
std::string m_reclusteringAlgorithmName
Name of daughter algorithm to use for cluster re-building.
float GetMatchedFraction() const
Get the fraction of sampling points resulting in a match.
const pandora::ClusterList & GetInputClusterList(const pandora::HitType hitType) const
static float MergeTwoPositions(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const float position1, const float position2)
Merge two views (U,V) to give a third view (Z).
const FitSegment & GetFitSegment(const float rL) const
Get fit segment for a given longitudinal coordinate.
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.
Header file for the three view fragments algorithm base class.
void CalculateOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW)
Calculate cluster overlap result and store in container.
bool CheckMatchedClusters(const pandora::CartesianPointVector &projectedPositions, const pandora::ClusterList &matchedClusters) const
Whether the matched clusters are consistent with the projected positions.
const pandora::ClusterList & GetSelectedClusterList(const pandora::HitType hitType) const
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
const OverlapResult & GetOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW) const
Get the overlap result for a specified trio of clusters.
HitType
Definition: HitType.h:12
pandora::StatusCode GetGlobalFitProjection(const pandora::CartesianVector &inputPosition, pandora::CartesianVector &projectedPosition) const
Get projected position on global fit for a given position vector.
TensorToolVector m_algorithmToolVector
The algorithm tool list.
pandora::StatusCode GetTransverseProjection(const float x, const FitSegment &fitSegment, pandora::CartesianVector &position) const
Get projected position for a given input x coordinate and fit segment.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
std::unordered_map< const pandora::CaloHit *, const pandora::Cluster * > HitToClusterMap
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
unsigned int m_nMaxTensorToolRepeats
The maximum number of repeat loops over tensor tools.
float m_minXOverlapFraction
requirement on minimum X overlap fraction for associated clusters
float m_maxPointDisplacementSquared
maximum allowed distance (squared) between projected points and associated hits
float m_minXOverlap
requirement on minimum X overlap for associated clusters
float m_minMatchedSamplingPointFraction
minimum fraction of matched sampling points
pandora::StatusCode GetMatchedHits(const pandora::ClusterList &inputClusterList, const pandora::CartesianPointVector &projectedPositions, HitToClusterMap &hitToClusterMap, pandora::CaloHitList &matchedCaloHits) const
Get the list of hits associated with the projected positions and a useful hit to cluster map...