All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CosmicRayTrackRecoveryAlgorithm.cc
Go to the documentation of this file.
1 /**
2  * @file larpandoracontent/LArTwoDReco/LArCosmicRay/CosmicRayTrackRecoveryAlgorithm.cc
3  *
4  * @brief Implementation of the cosmic ray longitudinal track recovery algorithm class.
5  *
6  * $Log: $
7  */
8 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 CosmicRayTrackRecoveryAlgorithm::CosmicRayTrackRecoveryAlgorithm() :
23  m_clusterMinLength(10.f),
24  m_clusterMinSpanZ(2.f),
25  m_clusterMinOverlapX(6.f),
26  m_clusterMaxDeltaX(3.f),
27  m_clusterMinHits(3)
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
34 {
35  // Get the available clusters for each view
36  ClusterVector availableClustersU, availableClustersV, availableClustersW;
37  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameU, availableClustersU));
38  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameV, availableClustersV));
39  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameW, availableClustersW));
40 
41  // Select clean clusters in each view
42  ClusterVector cleanClustersU, cleanClustersV, cleanClustersW;
43  this->SelectCleanClusters(availableClustersU, cleanClustersU);
44  this->SelectCleanClusters(availableClustersV, cleanClustersV);
45  this->SelectCleanClusters(availableClustersW, cleanClustersW);
46 
47  // Calculate sliding fit results for clean clusters
48  TwoDSlidingFitResultMap slidingFitResultMap;
49  this->BuildSlidingFitResultMap(cleanClustersU, slidingFitResultMap);
50  this->BuildSlidingFitResultMap(cleanClustersV, slidingFitResultMap);
51  this->BuildSlidingFitResultMap(cleanClustersW, slidingFitResultMap);
52 
53  // Match clusters between pairs of views (using start/end information)
54  ClusterAssociationMap matchedClustersUV, matchedClustersVW, matchedClustersWU;
55  this->MatchViews(cleanClustersU, cleanClustersV, slidingFitResultMap, matchedClustersUV);
56  this->MatchViews(cleanClustersV, cleanClustersW, slidingFitResultMap, matchedClustersVW);
57  this->MatchViews(cleanClustersW, cleanClustersU, slidingFitResultMap, matchedClustersWU);
58 
59  // Create candidate particles using one, two and three primary views
60  ParticleList candidateParticles;
61 
62  this->MatchThreeViews(cleanClustersU, cleanClustersV, cleanClustersW, matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
63 
64  this->MatchTwoViews(cleanClustersU, cleanClustersV, cleanClustersW, matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
65 
66  this->MatchOneView(cleanClustersU, cleanClustersV, cleanClustersW, matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
67 
68  // Build particle flow objects from candidate particles
69  this->BuildParticles(candidateParticles);
70 
71  return STATUS_CODE_SUCCESS;
72 }
73 
74 //------------------------------------------------------------------------------------------------------------------------------------------
75 
76 StatusCode CosmicRayTrackRecoveryAlgorithm::GetAvailableClusters(const std::string &inputClusterListName, ClusterVector &clusterVector) const
77 {
78  const ClusterList *pClusterList = NULL;
79  PANDORA_RETURN_RESULT_IF_AND_IF(
80  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, inputClusterListName, pClusterList))
81 
82  if (!pClusterList || pClusterList->empty())
83  {
84  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
85  std::cout << "CosmicRayTrackRecoveryAlgorithm: unable to find cluster list " << inputClusterListName << std::endl;
86 
87  return STATUS_CODE_SUCCESS;
88  }
89 
90  for (const Cluster *const pCluster : *pClusterList)
91  {
92  if (PandoraContentApi::IsAvailable(*this, pCluster))
93  clusterVector.push_back(pCluster);
94  }
95 
96  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
97 
98  return STATUS_CODE_SUCCESS;
99 }
100 
101 //------------------------------------------------------------------------------------------------------------------------------------------
102 
104 {
105  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
106  {
107  const Cluster *const pCluster = *iter;
108 
109  // Remove clusters with insufficient hits for a sliding fit
110  if (pCluster->GetNCaloHits() < m_clusterMinHits)
111  continue;
112  // Remove clusters below a minimum length
114  continue;
115 
116  // Remove clusters nearly parallel to Z or X
117  CartesianVector minCoordinate(0.f, 0.f, 0.f);
118  CartesianVector maxCoordinate(0.f, 0.f, 0.f);
119  LArClusterHelper::GetClusterBoundingBox(pCluster, minCoordinate, maxCoordinate);
120 
121  const CartesianVector deltaCoordinate(maxCoordinate - minCoordinate);
122  if (std::fabs(deltaCoordinate.GetZ()) < m_clusterMinSpanZ || std::fabs(deltaCoordinate.GetX()) < m_clusterMinOverlapX)
123  continue;
124 
125  outputVector.push_back(pCluster);
126  }
127 }
128 
129 //------------------------------------------------------------------------------------------------------------------------------------------
130 
132 {
133  const unsigned int m_halfWindowLayers(25);
134  const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
135 
136  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
137  {
138  if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
139  {
140  try
141  {
142  const TwoDSlidingFitResult slidingFitResult(*iter, m_halfWindowLayers, slidingFitPitch);
143 
144  if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
145  throw StatusCodeException(STATUS_CODE_FAILURE);
146  }
147  catch (StatusCodeException &statusCodeException)
148  {
149  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
150  throw statusCodeException;
151  }
152  }
153  }
154 }
155 
156 //------------------------------------------------------------------------------------------------------------------------------------------
157 
158 void CosmicRayTrackRecoveryAlgorithm::MatchViews(const ClusterVector &clusterVector1, const ClusterVector &clusterVector2,
159  const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
160 {
161  for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
162  this->MatchClusters(*iter1, clusterVector2, slidingFitResultMap, clusterAssociationMap);
163 
164  for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
165  this->MatchClusters(*iter2, clusterVector1, slidingFitResultMap, clusterAssociationMap);
166 }
167 
168 //------------------------------------------------------------------------------------------------------------------------------------------
169 
170 void CosmicRayTrackRecoveryAlgorithm::MatchClusters(const Cluster *const pSeedCluster, const ClusterVector &targetClusters,
171  const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
172 {
173  // Match seed cluster to target clusters according to alignment in X position of start/end positions
174  // Two possible matches: (a) one-to-one associations where both the track start and end positions match up
175  // (b) one-to-two associations where the track is split into two clusters in one view
176  // Require overlap in X (according to clusterMinOverlapX) and alignment in X (according to clusterMaxDeltaX)
177 
178  TwoDSlidingFitResultMap::const_iterator fsIter = slidingFitResultMap.find(pSeedCluster);
179 
180  if (slidingFitResultMap.end() == fsIter)
181  throw StatusCodeException(STATUS_CODE_FAILURE);
182 
183  const TwoDSlidingFitResult &slidingFitResult1(fsIter->second);
184  const CartesianVector &innerVertex1(slidingFitResult1.GetGlobalMinLayerPosition());
185  const CartesianVector &outerVertex1(slidingFitResult1.GetGlobalMaxLayerPosition());
186  const float xSpan1(std::fabs(outerVertex1.GetX() - innerVertex1.GetX()));
187 
188  const Cluster *pBestClusterInner(NULL);
189  const Cluster *pBestClusterOuter(NULL);
190  const Cluster *pBestCluster(NULL);
191 
192  float bestDisplacementInner(m_clusterMaxDeltaX);
193  float bestDisplacementOuter(m_clusterMaxDeltaX);
194  float bestDisplacement(2.f * m_clusterMaxDeltaX);
195 
196  for (ClusterVector::const_iterator tIter = targetClusters.begin(), tIterEnd = targetClusters.end(); tIter != tIterEnd; ++tIter)
197  {
198  const Cluster *const pTargetCluster = *tIter;
199 
200  UIntSet daughterVolumeIntersection;
201  LArGeometryHelper::GetCommonDaughterVolumes(pSeedCluster, pTargetCluster, daughterVolumeIntersection);
202 
203  if (daughterVolumeIntersection.empty())
204  continue;
205 
207  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
208 
209  TwoDSlidingFitResultMap::const_iterator ftIter = slidingFitResultMap.find(*tIter);
210 
211  if (slidingFitResultMap.end() == ftIter)
212  throw StatusCodeException(STATUS_CODE_FAILURE);
213 
214  const TwoDSlidingFitResult &slidingFitResult2(ftIter->second);
215  const CartesianVector &innerVertex2(slidingFitResult2.GetGlobalMinLayerPosition());
216  const CartesianVector &outerVertex2(slidingFitResult2.GetGlobalMaxLayerPosition());
217  const float xSpan2(std::fabs(outerVertex2.GetX() - innerVertex2.GetX()));
218 
219  if (xSpan2 > 1.5f * xSpan1)
220  continue;
221 
222  const float xMin1(std::min(innerVertex1.GetX(), outerVertex1.GetX()));
223  const float xMax1(std::max(innerVertex1.GetX(), outerVertex1.GetX()));
224  const float xMin2(std::min(innerVertex2.GetX(), outerVertex2.GetX()));
225  const float xMax2(std::max(innerVertex2.GetX(), outerVertex2.GetX()));
226  const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
227 
228  if (xOverlap < m_clusterMinOverlapX)
229  continue;
230 
231  const float dxMin(std::fabs(xMin2 - xMin1));
232  const float dxMax(std::fabs(xMax2 - xMax1));
233 
234  if (dxMin < bestDisplacementInner)
235  {
236  pBestClusterInner = pTargetCluster;
237  bestDisplacementInner = dxMin;
238  }
239 
240  if (dxMax < bestDisplacementOuter)
241  {
242  pBestClusterOuter = pTargetCluster;
243  bestDisplacementOuter = dxMax;
244  }
245 
246  if (dxMin + dxMax < bestDisplacement)
247  {
248  pBestCluster = pTargetCluster;
249  bestDisplacement = dxMin + dxMax;
250  }
251  }
252 
253  if (pBestCluster)
254  {
255  ClusterList &seedList(clusterAssociationMap[pSeedCluster]);
256 
257  if (seedList.end() == std::find(seedList.begin(), seedList.end(), pBestCluster))
258  seedList.push_back(pBestCluster);
259 
260  ClusterList &bestList(clusterAssociationMap[pBestCluster]);
261 
262  if (bestList.end() == std::find(bestList.begin(), bestList.end(), pSeedCluster))
263  bestList.push_back(pSeedCluster);
264  }
265  else if (pBestClusterInner && pBestClusterOuter)
266  {
267  TwoDSlidingFitResultMap::const_iterator iterInner = slidingFitResultMap.find(pBestClusterInner);
268  TwoDSlidingFitResultMap::const_iterator iterOuter = slidingFitResultMap.find(pBestClusterOuter);
269 
270  if (slidingFitResultMap.end() == iterInner || slidingFitResultMap.end() == iterOuter)
271  throw StatusCodeException(STATUS_CODE_FAILURE);
272 
273  const LArPointingCluster pointingClusterInner(iterInner->second);
274  const LArPointingCluster pointingClusterOuter(iterOuter->second);
275 
276  LArPointingCluster::Vertex pointingVertexInner, pointingVertexOuter;
277 
278  try
279  {
280  LArPointingClusterHelper::GetClosestVertices(pointingClusterInner, pointingClusterOuter, pointingVertexInner, pointingVertexOuter);
281  }
282  catch (StatusCodeException &)
283  {
284  return;
285  }
286 
287  const LArPointingCluster::Vertex pointingEndInner(
288  pointingVertexInner.IsInnerVertex() ? pointingClusterInner.GetOuterVertex() : pointingClusterInner.GetInnerVertex());
289  const LArPointingCluster::Vertex pointingEndOuter(
290  pointingVertexOuter.IsInnerVertex() ? pointingClusterOuter.GetOuterVertex() : pointingClusterOuter.GetInnerVertex());
291 
292  const float rSpan((pointingEndInner.GetPosition() - pointingEndOuter.GetPosition()).GetMagnitude());
293 
294  if (LArPointingClusterHelper::IsEmission(pointingVertexInner.GetPosition(), pointingVertexOuter, -1.f, 0.75f * rSpan, 5.f, 10.f) &&
295  LArPointingClusterHelper::IsEmission(pointingVertexOuter.GetPosition(), pointingVertexInner, -1.f, 0.75f * rSpan, 5.f, 10.f))
296  {
297  ClusterList &bestInnerList(clusterAssociationMap[pBestClusterInner]);
298 
299  if (bestInnerList.end() == std::find(bestInnerList.begin(), bestInnerList.end(), pSeedCluster))
300  bestInnerList.push_back(pSeedCluster);
301 
302  ClusterList &bestOuterList(clusterAssociationMap[pBestClusterOuter]);
303 
304  if (bestOuterList.end() == std::find(bestOuterList.begin(), bestOuterList.end(), pSeedCluster))
305  bestOuterList.push_back(pSeedCluster);
306  }
307  }
308 }
309 
310 //------------------------------------------------------------------------------------------------------------------------------------------
311 
312 void CosmicRayTrackRecoveryAlgorithm::MatchThreeViews(const ClusterVector &clusterVectorU, const ClusterVector &clusterVectorV,
313  const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW,
314  const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
315 {
316  ClusterSet vetoList;
317  this->BuildVetoList(particleList, vetoList);
318 
319  ParticleList newParticleList;
320 
321  const ClusterVector &clusterVector1(clusterVectorU);
322  const ClusterVector &clusterVector2(clusterVectorV);
323  const ClusterVector &clusterVector3(clusterVectorW);
324 
325  const ClusterAssociationMap &matchedClusters12(matchedClustersUV);
326  const ClusterAssociationMap &matchedClusters23(matchedClustersVW);
327  const ClusterAssociationMap &matchedClusters31(matchedClustersWU);
328 
329  for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
330  {
331  const Cluster *const pCluster1 = *iter1;
332 
333  if (vetoList.count(pCluster1))
334  continue;
335 
336  const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
337  const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());
338 
339  const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
340  const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());
341 
342  for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEndV = clusterVector2.end(); iter2 != iterEndV; ++iter2)
343  {
344  const Cluster *const pCluster2 = *iter2;
345 
346  if (vetoList.count(pCluster2))
347  continue;
348 
349  const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
350  const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());
351 
352  const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
353  const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());
354 
355  for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
356  {
357  const Cluster *const pCluster3 = *iter3;
358 
359  if (vetoList.count(pCluster3))
360  continue;
361 
362  const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
363  const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());
364 
365  const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
366  const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());
367 
368  const bool match12(
369  (matchedClusters12_pCluster1.size() + matchedClusters12_pCluster2.size() > 0) &&
370  ((matchedClusters12_pCluster1.size() == 1 && std::find(matchedClusters12_pCluster1.begin(), matchedClusters12_pCluster1.end(),
371  pCluster2) != matchedClusters12_pCluster1.end()) ||
372  (matchedClusters12_pCluster1.size() == 0)) &&
373  ((matchedClusters12_pCluster2.size() == 1 && std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(),
374  pCluster1) != matchedClusters12_pCluster2.end()) ||
375  (matchedClusters12_pCluster2.size() == 0)));
376 
377  const bool match23(
378  (matchedClusters23_pCluster2.size() + matchedClusters23_pCluster3.size() > 0) &&
379  ((matchedClusters23_pCluster2.size() == 1 && std::find(matchedClusters23_pCluster2.begin(), matchedClusters23_pCluster2.end(),
380  pCluster3) != matchedClusters23_pCluster2.end()) ||
381  (matchedClusters23_pCluster2.size() == 0)) &&
382  ((matchedClusters23_pCluster3.size() == 1 && std::find(matchedClusters23_pCluster3.begin(), matchedClusters23_pCluster3.end(),
383  pCluster2) != matchedClusters23_pCluster3.end()) ||
384  (matchedClusters23_pCluster3.size() == 0)));
385 
386  const bool match31(
387  (matchedClusters31_pCluster3.size() + matchedClusters31_pCluster1.size() > 0) &&
388  ((matchedClusters31_pCluster3.size() == 1 && std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(),
389  pCluster1) != matchedClusters31_pCluster3.end()) ||
390  (matchedClusters31_pCluster3.size() == 0)) &&
391  ((matchedClusters31_pCluster1.size() == 1 && std::find(matchedClusters31_pCluster1.begin(), matchedClusters31_pCluster1.end(),
392  pCluster3) != matchedClusters31_pCluster1.end()) ||
393  (matchedClusters31_pCluster1.size() == 0)));
394 
395  if (match12 && match23 && match31)
396  {
397  Particle newParticle;
398  newParticle.m_clusterList.push_back(pCluster1);
399  newParticle.m_clusterList.push_back(pCluster2);
400  newParticle.m_clusterList.push_back(pCluster3);
401  newParticleList.push_back(newParticle);
402  }
403  }
404  }
405  }
406 
407  return this->RemoveAmbiguities(newParticleList, particleList);
408 }
409 
410 //------------------------------------------------------------------------------------------------------------------------------------------
411 
412 void CosmicRayTrackRecoveryAlgorithm::MatchTwoViews(const ClusterVector &clusterVectorU, const ClusterVector &clusterVectorV,
413  const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW,
414  const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
415 {
416  ClusterSet vetoList;
417  this->BuildVetoList(particleList, vetoList);
418 
419  ParticleList newParticleList;
420 
421  for (unsigned int iView = 0; iView < 3; ++iView)
422  {
423  const ClusterVector &clusterVector1((0 == iView) ? clusterVectorU : (1 == iView) ? clusterVectorV : clusterVectorW);
424  const ClusterVector &clusterVector2((0 == iView) ? clusterVectorV : (1 == iView) ? clusterVectorW : clusterVectorU);
425  const ClusterVector &clusterVector3((0 == iView) ? clusterVectorW : (1 == iView) ? clusterVectorU : clusterVectorV);
426 
427  const ClusterAssociationMap &matchedClusters12(((0 == iView) ? matchedClustersUV : (1 == iView) ? matchedClustersVW : matchedClustersWU));
428  const ClusterAssociationMap &matchedClusters23(((0 == iView) ? matchedClustersVW : (1 == iView) ? matchedClustersWU : matchedClustersUV));
429  const ClusterAssociationMap &matchedClusters31(((0 == iView) ? matchedClustersWU : (1 == iView) ? matchedClustersUV : matchedClustersVW));
430 
431  for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
432  {
433  const Cluster *const pCluster1 = *iter1;
434 
435  if (vetoList.count(pCluster1))
436  continue;
437 
438  const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
439  const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());
440 
441  const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
442  const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());
443 
444  for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
445  {
446  const Cluster *const pCluster2 = *iter2;
447 
448  if (vetoList.count(pCluster2))
449  continue;
450 
451  const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
452  const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());
453 
454  const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
455  const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());
456 
457  const bool match12(
458  (matchedClusters12_pCluster1.size() == 1 && std::find(matchedClusters12_pCluster1.begin(), matchedClusters12_pCluster1.end(),
459  pCluster2) != matchedClusters12_pCluster1.end()) &&
460  (matchedClusters12_pCluster2.size() == 1 && std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(),
461  pCluster1) != matchedClusters12_pCluster2.end()) &&
462  (matchedClusters23_pCluster2.size() == 0 && matchedClusters31_pCluster1.size() == 0));
463 
464  if (!match12)
465  continue;
466 
467  Particle newParticle;
468  newParticle.m_clusterList.push_back(pCluster1);
469  newParticle.m_clusterList.push_back(pCluster2);
470 
471  for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
472  {
473  const Cluster *const pCluster3 = *iter3;
474 
475  if (vetoList.count(pCluster3))
476  continue;
477 
478  const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
479  const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());
480 
481  const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
482  const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());
483 
484  const bool match3((matchedClusters31_pCluster3.size() + matchedClusters23_pCluster3.size() > 0) &&
485  ((matchedClusters31_pCluster3.size() == 1 &&
486  std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(), pCluster1) !=
487  matchedClusters31_pCluster3.end()) ||
488  (matchedClusters31_pCluster3.size() == 0)) &&
489  ((matchedClusters23_pCluster3.size() == 1 &&
490  std::find(matchedClusters23_pCluster3.begin(), matchedClusters23_pCluster3.end(), pCluster2) !=
491  matchedClusters23_pCluster3.end()) ||
492  (matchedClusters23_pCluster3.size() == 0)));
493 
494  if (match3)
495  newParticle.m_clusterList.push_back(pCluster3);
496  }
497 
498  newParticleList.push_back(newParticle);
499  }
500  }
501  }
502 
503  return this->RemoveAmbiguities(newParticleList, particleList);
504 }
505 
506 //------------------------------------------------------------------------------------------------------------------------------------------
507 
508 void CosmicRayTrackRecoveryAlgorithm::MatchOneView(const ClusterVector &clusterVectorU, const ClusterVector &clusterVectorV,
509  const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW,
510  const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
511 {
512  ClusterSet vetoList;
513  this->BuildVetoList(particleList, vetoList);
514 
515  ParticleList newParticleList;
516 
517  for (unsigned int iView = 0; iView < 3; ++iView)
518  {
519  const ClusterVector &clusterVector1((0 == iView) ? clusterVectorU : (1 == iView) ? clusterVectorV : clusterVectorW);
520  const ClusterVector &clusterVector2((0 == iView) ? clusterVectorV : (1 == iView) ? clusterVectorW : clusterVectorU);
521  const ClusterVector &clusterVector3((0 == iView) ? clusterVectorW : (1 == iView) ? clusterVectorU : clusterVectorV);
522 
523  const ClusterAssociationMap &matchedClusters12(((0 == iView) ? matchedClustersUV : (1 == iView) ? matchedClustersVW : matchedClustersWU));
524  const ClusterAssociationMap &matchedClusters23(((0 == iView) ? matchedClustersVW : (1 == iView) ? matchedClustersWU : matchedClustersUV));
525  const ClusterAssociationMap &matchedClusters31(((0 == iView) ? matchedClustersWU : (1 == iView) ? matchedClustersUV : matchedClustersVW));
526 
527  for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
528  {
529  const Cluster *const pCluster1 = *iter1;
530 
531  if (vetoList.count(pCluster1))
532  continue;
533 
534  const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
535  const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());
536 
537  const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
538  const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());
539 
540  if (matchedClusters12_pCluster1.size() + matchedClusters31_pCluster1.size() > 0)
541  continue;
542 
543  Particle newParticle;
544  newParticle.m_clusterList.push_back(pCluster1);
545 
546  for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
547  {
548  const Cluster *const pCluster2 = *iter2;
549 
550  if (vetoList.count(pCluster2))
551  continue;
552 
553  const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
554  const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());
555 
556  const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
557  const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());
558 
559  if (matchedClusters12_pCluster2.size() == 1 &&
560  std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(), pCluster1) !=
561  matchedClusters12_pCluster2.end() &&
562  matchedClusters23_pCluster2.size() == 0)
563  newParticle.m_clusterList.push_back(pCluster2);
564  }
565 
566  for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
567  {
568  const Cluster *const pCluster3 = *iter3;
569 
570  if (vetoList.count(pCluster3))
571  continue;
572 
573  const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
574  const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());
575 
576  const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
577  const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());
578 
579  if (matchedClusters31_pCluster3.size() == 1 &&
580  std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(), pCluster1) !=
581  matchedClusters31_pCluster3.end() &&
582  matchedClusters23_pCluster3.size() == 0)
583  newParticle.m_clusterList.push_back(pCluster3);
584  }
585 
586  if (newParticle.m_clusterList.size() > 1)
587  newParticleList.push_back(newParticle);
588  }
589  }
590 
591  return this->RemoveAmbiguities(newParticleList, particleList);
592 }
593 
594 //------------------------------------------------------------------------------------------------------------------------------------------
595 
596 void CosmicRayTrackRecoveryAlgorithm::BuildVetoList(const ParticleList &particleList, ClusterSet &vetoList) const
597 {
598  for (ParticleList::const_iterator pIter = particleList.begin(), pIterEnd = particleList.end(); pIter != pIterEnd; ++pIter)
599  {
600  const Particle &particle = *pIter;
601  vetoList.insert(particle.m_clusterList.begin(), particle.m_clusterList.end());
602  }
603 }
604 
605 //------------------------------------------------------------------------------------------------------------------------------------------
606 
607 void CosmicRayTrackRecoveryAlgorithm::RemoveAmbiguities(const ParticleList &inputParticleList, ParticleList &outputParticleList) const
608 {
609  for (ParticleList::const_iterator pIter1 = inputParticleList.begin(), pIterEnd1 = inputParticleList.end(); pIter1 != pIterEnd1; ++pIter1)
610  {
611  const Particle &particle1 = *pIter1;
612  const ClusterList &clusterList1 = particle1.m_clusterList;
613 
614  try
615  {
616  bool isUnique(true);
617 
618  for (ParticleList::const_iterator pIter2 = pIter1, pIterEnd2 = inputParticleList.end(); pIter2 != pIterEnd2; ++pIter2)
619  {
620  const Particle &particle2 = *pIter2;
621  const ClusterList &clusterList2 = particle2.m_clusterList;
622 
623  ClusterSet duplicateSet;
624 
625  for (ClusterList::const_iterator cIter1 = clusterList1.begin(), cIterEnd1 = clusterList1.end(); cIter1 != cIterEnd1; ++cIter1)
626  {
627  const Cluster *pCluster = *cIter1;
628 
629  if (clusterList2.end() != std::find(clusterList2.begin(), clusterList2.end(), pCluster))
630  duplicateSet.insert(pCluster);
631  }
632 
633  if (duplicateSet.size() > 0 && clusterList1.size() != clusterList2.size())
634  {
635  isUnique = false;
636  break;
637  }
638  }
639 
640  if (!isUnique)
641  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
642 
643  outputParticleList.push_back(particle1);
644  }
645  catch (StatusCodeException &)
646  {
647  std::cout << " Warning in CosmicRayTrackRecoveryAlgorithm: found duplicate particles in candidate list " << std::endl;
648  }
649  }
650 }
651 
652 //------------------------------------------------------------------------------------------------------------------------------------------
653 
655 {
656  if (particleList.empty())
657  return;
658 
659  const PfoList *pPfoList(nullptr);
660  std::string pfoListName;
661  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
662 
663  for (const Particle &particle : particleList)
664  {
665  if (!PandoraContentApi::IsAvailable(*this, &particle.m_clusterList))
666  continue;
667 
668  ClusterList clusterList;
669  this->MergeClusters(particle.m_clusterList, clusterList);
670 
671  if (clusterList.empty())
672  throw StatusCodeException(STATUS_CODE_FAILURE);
673 
675  pfoParameters.m_particleId = MU_MINUS; // TRACK
676  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
677  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
678  pfoParameters.m_energy = 0.f;
679  pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
680  pfoParameters.m_clusterList = clusterList;
681 
682  const ParticleFlowObject *pPfo(nullptr);
683  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
684  }
685 
686  if (!pPfoList->empty())
687  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
688 }
689 
690 //------------------------------------------------------------------------------------------------------------------------------------------
691 
692 void CosmicRayTrackRecoveryAlgorithm::MergeClusters(const ClusterList &inputClusterList, ClusterList &outputClusterList) const
693 {
694  ClusterList clusterListU, clusterListV, clusterListW;
695  LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_U, clusterListU);
696  LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_V, clusterListV);
697  LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_W, clusterListW);
698 
699  for (unsigned int iView = 0; iView < 3; ++iView)
700  {
701  const ClusterList clusterList((0 == iView) ? clusterListU : (1 == iView) ? clusterListV : clusterListW);
702  const std::string inputClusterListName((0 == iView) ? m_inputClusterListNameU : (1 == iView) ? m_inputClusterListNameV : m_inputClusterListNameW);
703 
704  if (clusterList.empty())
705  continue;
706 
707  const Cluster *const pSeedCluster(clusterList.front());
708 
709  if (!PandoraContentApi::IsAvailable(*this, pSeedCluster))
710  throw StatusCodeException(STATUS_CODE_FAILURE);
711 
712  for (const Cluster *const pAssociatedCluster : clusterList)
713  {
714  if (pAssociatedCluster == pSeedCluster)
715  continue;
716 
717  if (!PandoraContentApi::IsAvailable(*this, pAssociatedCluster))
718  continue;
719 
720  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
721  PandoraContentApi::MergeAndDeleteClusters(*this, pSeedCluster, pAssociatedCluster, inputClusterListName, inputClusterListName));
722  }
723 
724  outputClusterList.push_back(pSeedCluster);
725  }
726 }
727 
728 //------------------------------------------------------------------------------------------------------------------------------------------
729 
730 StatusCode CosmicRayTrackRecoveryAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
731 {
732  PANDORA_RETURN_RESULT_IF_AND_IF(
733  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinLength", m_clusterMinLength));
734 
735  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinSpanZ", m_clusterMinSpanZ));
736 
737  PANDORA_RETURN_RESULT_IF_AND_IF(
738  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinOverlapX", m_clusterMinOverlapX));
739 
740  PANDORA_RETURN_RESULT_IF_AND_IF(
741  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMaxDeltaX", m_clusterMaxDeltaX));
742 
743  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinHits", m_clusterMinHits));
744 
745  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameU", m_inputClusterListNameU));
746  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameV", m_inputClusterListNameV));
747  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameW", m_inputClusterListNameW));
748  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));
749 
750  return STATUS_CODE_SUCCESS;
751 }
752 
753 } // namespace lar_content
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
void RemoveAmbiguities(const ParticleList &inputParticleList, ParticleList &outputParticleList) const
Remove particles with duplicate clusters.
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 MatchTwoViews(const pandora::ClusterVector &clusterVectorU, const pandora::ClusterVector &clusterVectorV, const pandora::ClusterVector &clusterVectorW, const ClusterAssociationMap &clusterAssociationMapUV, const ClusterAssociationMap &clusterAssociationMapVW, const ClusterAssociationMap &clusterAssociationMapWU, ParticleList &particleList) const
Create candidate particles using two primary clusters and one pair of broken clusters.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterAssociationMap
void BuildParticles(const ParticleList &particleList)
Build particle flow objects.
static bool IsEmission(const pandora::CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex, const float minLongitudinalDistance, const float maxLongitudinalDistance, const float maxTransverseDistance, const float angularAllowance)
Whether pointing vertex is emitted from a given position.
void MatchOneView(const pandora::ClusterVector &clusterVectorU, const pandora::ClusterVector &clusterVectorV, const pandora::ClusterVector &clusterVectorW, const ClusterAssociationMap &clusterAssociationMapUV, const ClusterAssociationMap &clusterAssociationMapVW, const ClusterAssociationMap &clusterAssociationMapWU, ParticleList &particleList) const
Create candidate particles using one primary cluster and one pair of broken clusters.
void MatchClusters(const pandora::Cluster *const pSeedCluster, const pandora::ClusterVector &targetClusters, const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
Match a seed cluster with a list of target clusters and populate the cluster association map...
void BuildVetoList(const ParticleList &particleList, pandora::ClusterSet &vetoList) const
Build the list of clusters already used to create particles.
static void GetClustersByHitType(const pandora::ClusterList &inputClusters, const pandora::HitType hitType, pandora::ClusterList &clusterList)
Get the subset of clusters, from a provided list, that match the specified hit type.
LArPointingCluster class.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
static void GetClosestVertices(const bool useX, const bool useY, const bool useZ, const LArPointingCluster &pointingClusterI, const LArPointingCluster &pointingClusterJ, LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
Given a pair of pointing clusters, receive the closest or farthest pair of vertices.
Header file for the geometry helper class.
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.
void MatchThreeViews(const pandora::ClusterVector &clusterVectorU, const pandora::ClusterVector &clusterVectorV, const pandora::ClusterVector &clusterVectorW, const ClusterAssociationMap &clusterAssociationMapUV, const ClusterAssociationMap &clusterAssociationMapVW, const ClusterAssociationMap &clusterAssociationMapWU, ParticleList &particleList) const
Create candidate particles using three primary clusters.
pandora::StatusCode GetAvailableClusters(const std::string &inputClusterListName, pandora::ClusterVector &clusterVector) const
Get a vector of available clusters.
Header file for the cluster helper class.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
void SelectCleanClusters(const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select a set of clusters judged to be clean.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
fhicl::Table< sbnd::crt::CRTDetSimParams > Parameters
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
Definition: reco_sbnd.fcl:182
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
bool IsInnerVertex() const
Is this the inner vertex.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void MatchViews(const pandora::ClusterVector &clusterVector1, const pandora::ClusterVector &clusterVector2, const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
Match a pair of cluster vectors and populate the cluster association map.
void MergeClusters(const pandora::ClusterList &inputClusterList, pandora::ClusterList &outputClusterList) const
Merge broken clusters into a single cluster.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
BEGIN_PROLOG could also be cout
static void GetCommonDaughterVolumes(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, UIntSet &intersect)
Return the set of common daughter volumes between two 2D clusters.
Header file for the cosmic ray longitudinal track recovery algorithm class.