All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
VertexBasedPfoRecoveryAlgorithm.cc
Go to the documentation of this file.
1 /**
2  * @file larpandoracontent/LArThreeDReco/LArPfoRecovery/VertexBasedPfoRecoveryAlgorithm.cc
3  *
4  * @brief Implementation of the vertex-based particle recovery algorithm
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 VertexBasedPfoRecoveryAlgorithm::VertexBasedPfoRecoveryAlgorithm() :
23  m_slidingFitHalfWindow(10),
24  m_maxLongitudinalDisplacement(5.f),
25  m_maxTransverseDisplacement(2.f),
26  m_twoViewChi2Cut(5.f),
27  m_threeViewChi2Cut(5.f)
28 {
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
34 {
35  const VertexList *pVertexList = NULL;
36  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetCurrentList(*this, pVertexList));
37 
38  const Vertex *const pSelectedVertex(
39  (pVertexList && (pVertexList->size() == 1) && (VERTEX_3D == (*(pVertexList->begin()))->GetVertexType())) ? *(pVertexList->begin()) : NULL);
40 
41  if (!pSelectedVertex)
42  {
43  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
44  std::cout << "VertexBasedPfoRecoveryAlgorithm: unable to find vertex in current list " << std::endl;
45 
46  return STATUS_CODE_SUCCESS;
47  }
48 
49  // Get the available clusters from each view
50  ClusterVector availableClusters;
51  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNames, availableClusters));
52 
53  // Build a set of sliding fit results
54  TwoDSlidingFitResultMap slidingFitResultMap;
55  this->BuildSlidingFitResultMap(availableClusters, slidingFitResultMap);
56 
57  // Select seed clusters (adjacent to vertex)
58  ClusterVector selectedClusters;
59  this->SelectVertexClusters(pSelectedVertex, slidingFitResultMap, availableClusters, selectedClusters);
60 
61  // Match the cluster end points
62  ClusterSet vetoList;
63  ParticleList particleList;
64  this->MatchThreeViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
65  this->MatchTwoViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
66 
67  // Build new particles
68  this->BuildParticles(particleList);
69 
70  return STATUS_CODE_SUCCESS;
71 }
72 
73 //------------------------------------------------------------------------------------------------------------------------------------------
74 
75 StatusCode VertexBasedPfoRecoveryAlgorithm::GetAvailableClusters(const StringVector inputClusterListNames, ClusterVector &clusterVector) const
76 {
77  for (StringVector::const_iterator iter = inputClusterListNames.begin(), iterEnd = inputClusterListNames.end(); iter != iterEnd; ++iter)
78  {
79  const ClusterList *pClusterList(NULL);
80  PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, *iter, pClusterList));
81 
82  if (NULL == pClusterList)
83  {
84  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
85  std::cout << "VertexBasedPfoRecoveryAlgorithm: could not find cluster list " << *iter << std::endl;
86  continue;
87  }
88 
89  for (ClusterList::const_iterator cIter = pClusterList->begin(), cIterEnd = pClusterList->end(); cIter != cIterEnd; ++cIter)
90  {
91  const Cluster *const pCluster = *cIter;
92 
93  if (!pCluster->IsAvailable())
94  continue;
95 
96  if (pCluster->GetNCaloHits() <= 1)
97  continue;
98 
99  if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
100  continue;
101 
102  clusterVector.push_back(pCluster);
103  }
104  }
105 
106  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
107 
108  return STATUS_CODE_SUCCESS;
109 }
110 //------------------------------------------------------------------------------------------------------------------------------------------
111 
113 {
114  const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
115 
116  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
117  {
118  if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
119  {
120  try
121  {
122  const TwoDSlidingFitResult slidingFitResult(*iter, m_slidingFitHalfWindow, slidingFitPitch);
123  const LArPointingCluster pointingCluster(slidingFitResult);
124 
125  if (pointingCluster.GetLengthSquared() < std::numeric_limits<float>::epsilon())
126  continue;
127 
128  if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
129  throw StatusCodeException(STATUS_CODE_FAILURE);
130  }
131  catch (StatusCodeException &statusCodeException)
132  {
133  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
134  throw statusCodeException;
135  }
136  }
137  }
138 }
139 
140 //------------------------------------------------------------------------------------------------------------------------------------------
141 
142 void VertexBasedPfoRecoveryAlgorithm::SelectVertexClusters(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
143  const ClusterVector &inputClusters, ClusterVector &outputClusters) const
144 {
145  const CartesianVector vertexU(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_U));
146  const CartesianVector vertexV(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_V));
147  const CartesianVector vertexW(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_W));
148 
149  for (ClusterVector::const_iterator cIter = inputClusters.begin(), cIterEnd = inputClusters.end(); cIter != cIterEnd; ++cIter)
150  {
151  const Cluster *const pCluster = *cIter;
152  const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
153 
154  if (TPC_3D == hitType)
155  continue;
156 
157  const CartesianVector vertexPosition((TPC_VIEW_U == hitType) ? vertexU : (TPC_VIEW_V == hitType) ? vertexV : vertexW);
158 
159  TwoDSlidingFitResultMap::const_iterator sIter = slidingFitResultMap.find(pCluster);
160  if (slidingFitResultMap.end() == sIter)
161  continue;
162 
163  const TwoDSlidingFitResult &slidingFitResult = sIter->second;
164  const LArPointingCluster pointingCluster(slidingFitResult);
165 
166  for (unsigned int iVtx = 0; iVtx < 2; ++iVtx)
167  {
168  const LArPointingCluster::Vertex &pointingVertex((0 == iVtx) ? pointingCluster.GetInnerVertex() : pointingCluster.GetOuterVertex());
169 
170  float rL(0.f), rT(0.f);
171  LArPointingClusterHelper::GetImpactParameters(pointingVertex, vertexPosition, rL, rT);
172 
173  if (rL > -1.f && rL < m_maxLongitudinalDisplacement && rT < m_maxTransverseDisplacement)
174  {
175  outputClusters.push_back(pCluster);
176  break;
177  }
178  }
179  }
180 }
181 
182 //------------------------------------------------------------------------------------------------------------------------------------------
183 
184 void VertexBasedPfoRecoveryAlgorithm::MatchThreeViews(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
185  const ClusterVector &inputClusters, ClusterSet &vetoList, ParticleList &particleList) const
186 {
187  while (true)
188  {
189  ClusterVector availableClusters, clustersU, clustersV, clustersW;
190  this->SelectAvailableClusters(vetoList, inputClusters, availableClusters);
191  this->SelectClusters(TPC_VIEW_U, availableClusters, clustersU);
192  this->SelectClusters(TPC_VIEW_V, availableClusters, clustersV);
193  this->SelectClusters(TPC_VIEW_W, availableClusters, clustersW);
194 
195  float chi2(m_threeViewChi2Cut);
196  const Cluster *pCluster1(NULL);
197  const Cluster *pCluster2(NULL);
198  const Cluster *pCluster3(NULL);
199 
200  this->GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, clustersW, pCluster1, pCluster2, pCluster3, chi2);
201 
202  if (NULL == pCluster1 || NULL == pCluster2 || NULL == pCluster3)
203  return;
204 
205  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
206  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
207  const HitType hitType3(LArClusterHelper::GetClusterHitType(pCluster3));
208 
209  const Cluster *const pClusterU(
210  (TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : (TPC_VIEW_U == hitType3) ? pCluster3 : NULL);
211  const Cluster *const pClusterV(
212  (TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : (TPC_VIEW_V == hitType3) ? pCluster3 : NULL);
213  const Cluster *const pClusterW(
214  (TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : (TPC_VIEW_W == hitType3) ? pCluster3 : NULL);
215 
216  particleList.push_back(Particle(pClusterU, pClusterV, pClusterW));
217 
218  vetoList.insert(pCluster1);
219  vetoList.insert(pCluster2);
220  vetoList.insert(pCluster3);
221  }
222 }
223 
224 //------------------------------------------------------------------------------------------------------------------------------------------
225 
226 void VertexBasedPfoRecoveryAlgorithm::MatchTwoViews(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
227  const ClusterVector &inputClusters, ClusterSet &vetoList, ParticleList &particleList) const
228 {
229  while (true)
230  {
231  ClusterVector availableClusters, clustersU, clustersV, clustersW;
232  this->SelectAvailableClusters(vetoList, inputClusters, availableClusters);
233  this->SelectClusters(TPC_VIEW_U, availableClusters, clustersU);
234  this->SelectClusters(TPC_VIEW_V, availableClusters, clustersV);
235  this->SelectClusters(TPC_VIEW_W, availableClusters, clustersW);
236 
237  float chi2(m_twoViewChi2Cut);
238  const Cluster *pCluster1(NULL);
239  const Cluster *pCluster2(NULL);
240 
241  this->GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, pCluster1, pCluster2, chi2);
242  this->GetBestChi2(pVertex, slidingFitResultMap, clustersV, clustersW, pCluster1, pCluster2, chi2);
243  this->GetBestChi2(pVertex, slidingFitResultMap, clustersW, clustersU, pCluster1, pCluster2, chi2);
244 
245  if (NULL == pCluster1 || NULL == pCluster2)
246  return;
247 
248  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
249  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
250 
251  const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
252  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
253  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
254 
255  particleList.push_back(Particle(pClusterU, pClusterV, pClusterW));
256 
257  vetoList.insert(pCluster1);
258  vetoList.insert(pCluster2);
259  }
260 }
261 
262 //------------------------------------------------------------------------------------------------------------------------------------------
263 
264 void VertexBasedPfoRecoveryAlgorithm::GetBestChi2(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
265  const ClusterVector &clusters1, const ClusterVector &clusters2, const ClusterVector &clusters3, const Cluster *&pBestCluster1,
266  const Cluster *&pBestCluster2, const Cluster *&pBestCluster3, float &bestChi2) const
267 {
268  if (clusters1.empty() || clusters2.empty() || clusters3.empty())
269  return;
270 
271  // First loop
272  for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
273  {
274  const Cluster *const pCluster1 = *cIter1;
275 
276  TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
277  if (slidingFitResultMap.end() == sIter1)
278  continue;
279 
280  const TwoDSlidingFitResult &slidingFitResult1 = sIter1->second;
281  const LArPointingCluster pointingCluster1(slidingFitResult1);
282 
283  // Second loop
284  for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
285  {
286  const Cluster *const pCluster2 = *cIter2;
287 
288  TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
289  if (slidingFitResultMap.end() == sIter2)
290  continue;
291 
292  const TwoDSlidingFitResult &slidingFitResult2 = sIter2->second;
293  const LArPointingCluster pointingCluster2(slidingFitResult2);
294 
295  // Third loop
296  for (ClusterVector::const_iterator cIter3 = clusters3.begin(), cIterEnd3 = clusters3.end(); cIter3 != cIterEnd3; ++cIter3)
297  {
298  const Cluster *const pCluster3 = *cIter3;
299 
300  TwoDSlidingFitResultMap::const_iterator sIter3 = slidingFitResultMap.find(pCluster3);
301  if (slidingFitResultMap.end() == sIter3)
302  continue;
303 
304  const TwoDSlidingFitResult &slidingFitResult3 = sIter3->second;
305  const LArPointingCluster pointingCluster3(slidingFitResult3);
306 
307  // Calculate chi-squared
308  const float thisChi2(this->GetChi2(pVertex, pointingCluster1, pointingCluster2, pointingCluster3));
309 
310  if (thisChi2 < bestChi2)
311  {
312  bestChi2 = thisChi2;
313  pBestCluster1 = pCluster1;
314  pBestCluster2 = pCluster2;
315  pBestCluster3 = pCluster3;
316  }
317  }
318  }
319  }
320 }
321 
322 //------------------------------------------------------------------------------------------------------------------------------------------
323 
324 void VertexBasedPfoRecoveryAlgorithm::GetBestChi2(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
325  const ClusterVector &clusters1, const ClusterVector &clusters2, const Cluster *&pBestCluster1, const Cluster *&pBestCluster2, float &bestChi2) const
326 {
327  if (clusters1.empty() || clusters2.empty())
328  return;
329 
330  // First loop
331  for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
332  {
333  const Cluster *const pCluster1 = *cIter1;
334 
335  TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
336  if (slidingFitResultMap.end() == sIter1)
337  continue;
338 
339  const TwoDSlidingFitResult &slidingFitResult1 = sIter1->second;
340  const LArPointingCluster pointingCluster1(slidingFitResult1);
341 
342  // Second loop
343  for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
344  {
345  const Cluster *const pCluster2 = *cIter2;
346 
347  TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
348  if (slidingFitResultMap.end() == sIter2)
349  continue;
350 
351  const TwoDSlidingFitResult &slidingFitResult2 = sIter2->second;
352  const LArPointingCluster pointingCluster2(slidingFitResult2);
353 
354  // Calculate chi-squared
355  const float thisChi2(this->GetChi2(pVertex, pointingCluster1, pointingCluster2));
356 
357  if (thisChi2 < bestChi2)
358  {
359  bestChi2 = thisChi2;
360  pBestCluster1 = pCluster1;
361  pBestCluster2 = pCluster2;
362  }
363  }
364  }
365 }
366 
367 //------------------------------------------------------------------------------------------------------------------------------------------
368 
370  const Vertex *const pVertex, const LArPointingCluster &pointingCluster1, const LArPointingCluster &pointingCluster2) const
371 {
372  const HitType hitType1(LArClusterHelper::GetClusterHitType(pointingCluster1.GetCluster()));
373  const HitType hitType2(LArClusterHelper::GetClusterHitType(pointingCluster2.GetCluster()));
374 
375  if (hitType1 == hitType2)
376  throw StatusCodeException(STATUS_CODE_FAILURE);
377 
378  const CartesianVector vertex1(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType1));
379  const CartesianVector vertex2(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType2));
380 
381  const LArPointingCluster::Vertex &pointingVertex1(this->GetOuterVertex(vertex1, pointingCluster1));
382  const LArPointingCluster::Vertex &pointingVertex2(this->GetOuterVertex(vertex2, pointingCluster2));
383 
384  float chi2(0.f);
385  CartesianVector mergedPosition(0.f, 0.f, 0.f);
387  this->GetPandora(), hitType1, hitType2, pointingVertex1.GetPosition(), pointingVertex2.GetPosition(), mergedPosition, chi2);
388 
389  return chi2;
390 }
391 
392 //------------------------------------------------------------------------------------------------------------------------------------------
393 
394 float VertexBasedPfoRecoveryAlgorithm::GetChi2(const Vertex *const pVertex, const LArPointingCluster &pointingCluster1,
395  const LArPointingCluster &pointingCluster2, const LArPointingCluster &pointingCluster3) const
396 {
397  const HitType hitType1(LArClusterHelper::GetClusterHitType(pointingCluster1.GetCluster()));
398  const HitType hitType2(LArClusterHelper::GetClusterHitType(pointingCluster2.GetCluster()));
399  const HitType hitType3(LArClusterHelper::GetClusterHitType(pointingCluster3.GetCluster()));
400 
401  if ((hitType1 == hitType2) || (hitType2 == hitType3) || (hitType3 == hitType1))
402  throw StatusCodeException(STATUS_CODE_FAILURE);
403 
404  const CartesianVector vertex1(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType1));
405  const CartesianVector vertex2(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType2));
406  const CartesianVector vertex3(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType3));
407 
408  const LArPointingCluster::Vertex &pointingVertex1(this->GetOuterVertex(vertex1, pointingCluster1));
409  const LArPointingCluster::Vertex &pointingVertex2(this->GetOuterVertex(vertex2, pointingCluster2));
410  const LArPointingCluster::Vertex &pointingVertex3(this->GetOuterVertex(vertex3, pointingCluster3));
411 
412  float chi2(0.f);
413  CartesianVector mergedPosition(0.f, 0.f, 0.f);
414  LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), hitType1, hitType2, hitType3, pointingVertex1.GetPosition(),
415  pointingVertex2.GetPosition(), pointingVertex3.GetPosition(), mergedPosition, chi2);
416 
417  return chi2;
418 }
419 
420 //------------------------------------------------------------------------------------------------------------------------------------------
421 
422 void VertexBasedPfoRecoveryAlgorithm::SelectAvailableClusters(const ClusterSet &vetoList, const ClusterVector &inputVector, ClusterVector &outputVector) const
423 {
424  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
425  {
426  if (0 == vetoList.count(*iter))
427  outputVector.push_back(*iter);
428  }
429 }
430 
431 //------------------------------------------------------------------------------------------------------------------------------------------
432 
433 void VertexBasedPfoRecoveryAlgorithm::SelectClusters(const HitType hitType, const ClusterVector &inputVector, ClusterVector &outputVector) const
434 {
435  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
436  {
437  if (hitType == LArClusterHelper::GetClusterHitType(*iter))
438  outputVector.push_back(*iter);
439  }
440 }
441 
442 //------------------------------------------------------------------------------------------------------------------------------------------
443 
445 {
446  const float innerDistance((vertex - cluster.GetInnerVertex().GetPosition()).GetMagnitudeSquared());
447  const float outerDistance((vertex - cluster.GetOuterVertex().GetPosition()).GetMagnitudeSquared());
448 
449  if (innerDistance < outerDistance)
450  return cluster.GetInnerVertex();
451  else
452  return cluster.GetOuterVertex();
453 }
454 
455 //------------------------------------------------------------------------------------------------------------------------------------------
456 
458 {
459  const LArPointingCluster::Vertex &innerVertex = this->GetInnerVertex(vertex, cluster);
460 
461  if (innerVertex.IsInnerVertex())
462  return cluster.GetOuterVertex();
463  else
464  return cluster.GetInnerVertex();
465 }
466 
467 //------------------------------------------------------------------------------------------------------------------------------------------
468 
470 {
471  if (particleList.empty())
472  return;
473 
474  const PfoList *pPfoList = NULL;
475  std::string pfoListName;
476  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
477 
478  for (ParticleList::const_iterator iter = particleList.begin(), iterEnd = particleList.end(); iter != iterEnd; ++iter)
479  {
480  const Particle &particle = *iter;
481 
482  ClusterList clusterList;
483  const Cluster *const pClusterU = particle.m_pClusterU;
484  const Cluster *const pClusterV = particle.m_pClusterV;
485  const Cluster *const pClusterW = particle.m_pClusterW;
486 
487  const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() : true);
488  const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() : true);
489  const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() : true);
490 
491  if (!(isAvailableU && isAvailableV && isAvailableW))
492  throw StatusCodeException(STATUS_CODE_FAILURE);
493 
494  if (pClusterU)
495  clusterList.push_back(pClusterU);
496  if (pClusterV)
497  clusterList.push_back(pClusterV);
498  if (pClusterW)
499  clusterList.push_back(pClusterW);
500 
501  // TODO Correct these placeholder parameters
503  pfoParameters.m_particleId = MU_MINUS; // TRACK
504  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
505  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
506  pfoParameters.m_energy = 0.f;
507  pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
508  pfoParameters.m_clusterList = clusterList;
509 
510  const ParticleFlowObject *pPfo(NULL);
511  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
512  }
513 
514  if (!pPfoList->empty())
515  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
516 }
517 
518 //------------------------------------------------------------------------------------------------------------------------------------------
519 
520 VertexBasedPfoRecoveryAlgorithm::Particle::Particle(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW) :
521  m_pClusterU(pClusterU),
522  m_pClusterV(pClusterV),
523  m_pClusterW(pClusterW)
524 {
525  if (NULL == m_pClusterU && NULL == m_pClusterV && NULL == m_pClusterW)
526  throw StatusCodeException(STATUS_CODE_FAILURE);
527 
528  const HitType hitTypeU(NULL == m_pClusterU ? TPC_VIEW_U : LArClusterHelper::GetClusterHitType(m_pClusterU));
529  const HitType hitTypeV(NULL == m_pClusterV ? TPC_VIEW_V : LArClusterHelper::GetClusterHitType(m_pClusterV));
530  const HitType hitTypeW(NULL == m_pClusterW ? TPC_VIEW_W : LArClusterHelper::GetClusterHitType(m_pClusterW));
531 
532  if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
533  throw StatusCodeException(STATUS_CODE_FAILURE);
534 }
535 
536 //------------------------------------------------------------------------------------------------------------------------------------------
537 
538 StatusCode VertexBasedPfoRecoveryAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
539 {
540  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));
541 
542  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));
543 
544  PANDORA_RETURN_RESULT_IF_AND_IF(
545  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitHalfWindow", m_slidingFitHalfWindow));
546 
547  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
548  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
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, !=, XmlHelper::ReadValue(xmlHandle, "TwoViewChi2Cut", m_twoViewChi2Cut));
554 
555  PANDORA_RETURN_RESULT_IF_AND_IF(
556  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ThreeViewChi2Cut", m_threeViewChi2Cut));
557 
558  return STATUS_CODE_SUCCESS;
559 }
560 
561 } // namespace lar_content
process_name vertex
Definition: cheaterreco.fcl:51
std::string m_outputPfoListName
The name of the output pfo list.
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.
const LArPointingCluster::Vertex & GetInnerVertex(const pandora::CartesianVector &vertex, const LArPointingCluster &cluster) const
Find nearest end of pointing cluster to a specified position vector.
void GetBestChi2(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &clusters1, const pandora::ClusterVector &clusters2, const pandora::ClusterVector &clusters3, const pandora::Cluster *&pBestCluster1, const pandora::Cluster *&pBestCluster2, const pandora::Cluster *&pBestCluster3, float &chi2) const
Get best-matched triplet of clusters from a set of input cluster vectors.
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.
process_name cluster
Definition: cheaterreco.fcl:51
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
Particle(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW)
Constructor.
void SelectVertexClusters(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &inputClusters, pandora::ClusterVector &outputClusters) const
Select clusters in proximity to reconstructed vertex.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void MatchTwoViews(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &selectedClusters, pandora::ClusterSet &vetoList, ParticleList &particleList) const
Match clusters from two views.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
LArPointingCluster class.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
void MatchThreeViews(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &selectedClusters, pandora::ClusterSet &vetoList, ParticleList &particleList) const
Match clusters from three views.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
Header file for the geometry helper class.
pandora::StringVector m_inputClusterListNames
The list of input cluster list names.
Header file for the cluster helper class.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
static void MergeThreePositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::HitType view3, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, const pandora::CartesianVector &position3, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from three views to give unified 3D position.
float GetChi2(const pandora::Vertex *const pVertex, const LArPointingCluster &pointingCluster1, const LArPointingCluster &pointingCluster2) const
Merge two pointing clusters and return chi-squared metric giving consistency of matching.
float GetLengthSquared() const
Get length squared of pointing cluster.
pandora::StatusCode GetAvailableClusters(const pandora::StringVector inputClusterListName, pandora::ClusterVector &clusterVector) const
Get a vector of available clusters.
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
const LArPointingCluster::Vertex & GetOuterVertex(const pandora::CartesianVector &vertex, const LArPointingCluster &cluster) const
Find furthest end of pointing cluster from a specified position vector.
fhicl::Table< sbnd::crt::CRTDetSimParams > Parameters
void SelectAvailableClusters(const pandora::ClusterSet &vetoList, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select cluster which haven&#39;t been vetoed.
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
Definition: reco_sbnd.fcl:182
const pandora::Cluster * m_pClusterV
Address of cluster in V view.
void SelectClusters(const pandora::HitType hitType, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select clusters of a specified hit type.
bool IsInnerVertex() const
Is this the inner vertex.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
Header file for the vertex-based particle recovery algorithm.
const pandora::Cluster * m_pClusterW
Address of cluster in W view.
std::list< Vertex > VertexList
Definition: DCEL.h:182
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
void BuildParticles(const ParticleList &particleList)
Build particle flow objects from matched clusters.
BEGIN_PROLOG could also be cout
const pandora::Cluster * m_pClusterU
Address of cluster in U view.