LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
VertexBasedPfoRecoveryAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 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  for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
115  {
116  if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
117  {
118  try
119  {
120  const float slidingFitPitch(LArGeometryHelper::GetWirePitch(this->GetPandora(), LArClusterHelper::GetClusterHitType(*iter)));
121  const TwoDSlidingFitResult slidingFitResult(*iter, m_slidingFitHalfWindow, slidingFitPitch);
122  const LArPointingCluster pointingCluster(slidingFitResult);
123 
124  if (pointingCluster.GetLengthSquared() < std::numeric_limits<float>::epsilon())
125  continue;
126 
127  if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
128  throw StatusCodeException(STATUS_CODE_FAILURE);
129  }
130  catch (StatusCodeException &statusCodeException)
131  {
132  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
133  throw statusCodeException;
134  }
135  }
136  }
137 }
138 
139 //------------------------------------------------------------------------------------------------------------------------------------------
140 
141 void VertexBasedPfoRecoveryAlgorithm::SelectVertexClusters(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
142  const ClusterVector &inputClusters, ClusterVector &outputClusters) const
143 {
144  const CartesianVector vertexU(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_U));
145  const CartesianVector vertexV(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_V));
146  const CartesianVector vertexW(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), TPC_VIEW_W));
147 
148  for (ClusterVector::const_iterator cIter = inputClusters.begin(), cIterEnd = inputClusters.end(); cIter != cIterEnd; ++cIter)
149  {
150  const Cluster *const pCluster = *cIter;
151  const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
152 
153  if (TPC_3D == hitType)
154  continue;
155 
156  const CartesianVector vertexPosition((TPC_VIEW_U == hitType) ? vertexU : (TPC_VIEW_V == hitType) ? vertexV : vertexW);
157 
158  TwoDSlidingFitResultMap::const_iterator sIter = slidingFitResultMap.find(pCluster);
159  if (slidingFitResultMap.end() == sIter)
160  continue;
161 
162  const TwoDSlidingFitResult &slidingFitResult = sIter->second;
163  const LArPointingCluster pointingCluster(slidingFitResult);
164 
165  for (unsigned int iVtx = 0; iVtx < 2; ++iVtx)
166  {
167  const LArPointingCluster::Vertex &pointingVertex((0 == iVtx) ? pointingCluster.GetInnerVertex() : pointingCluster.GetOuterVertex());
168 
169  float rL(0.f), rT(0.f);
170  LArPointingClusterHelper::GetImpactParameters(pointingVertex, vertexPosition, rL, rT);
171 
172  if (rL > -1.f && rL < m_maxLongitudinalDisplacement && rT < m_maxTransverseDisplacement)
173  {
174  outputClusters.push_back(pCluster);
175  break;
176  }
177  }
178  }
179 }
180 
181 //------------------------------------------------------------------------------------------------------------------------------------------
182 
183 void VertexBasedPfoRecoveryAlgorithm::MatchThreeViews(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
184  const ClusterVector &inputClusters, ClusterSet &vetoList, ParticleList &particleList) const
185 {
186  while (true)
187  {
188  ClusterVector availableClusters, clustersU, clustersV, clustersW;
189  this->SelectAvailableClusters(vetoList, inputClusters, availableClusters);
190  this->SelectClusters(TPC_VIEW_U, availableClusters, clustersU);
191  this->SelectClusters(TPC_VIEW_V, availableClusters, clustersV);
192  this->SelectClusters(TPC_VIEW_W, availableClusters, clustersW);
193 
194  float chi2(m_threeViewChi2Cut);
195  const Cluster *pCluster1(NULL);
196  const Cluster *pCluster2(NULL);
197  const Cluster *pCluster3(NULL);
198 
199  this->GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, clustersW, pCluster1, pCluster2, pCluster3, chi2);
200 
201  if (NULL == pCluster1 || NULL == pCluster2 || NULL == pCluster3)
202  return;
203 
204  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
205  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
206  const HitType hitType3(LArClusterHelper::GetClusterHitType(pCluster3));
207 
208  const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1
209  : (TPC_VIEW_U == hitType2) ? pCluster2
210  : (TPC_VIEW_U == hitType3) ? pCluster3
211  : NULL);
212  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1
213  : (TPC_VIEW_V == hitType2) ? pCluster2
214  : (TPC_VIEW_V == hitType3) ? pCluster3
215  : NULL);
216  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1
217  : (TPC_VIEW_W == hitType2) ? pCluster2
218  : (TPC_VIEW_W == hitType3) ? pCluster3
219  : NULL);
220 
221  particleList.push_back(Particle(pClusterU, pClusterV, pClusterW));
222 
223  vetoList.insert(pCluster1);
224  vetoList.insert(pCluster2);
225  vetoList.insert(pCluster3);
226  }
227 }
228 
229 //------------------------------------------------------------------------------------------------------------------------------------------
230 
231 void VertexBasedPfoRecoveryAlgorithm::MatchTwoViews(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
232  const ClusterVector &inputClusters, ClusterSet &vetoList, ParticleList &particleList) const
233 {
234  while (true)
235  {
236  ClusterVector availableClusters, clustersU, clustersV, clustersW;
237  this->SelectAvailableClusters(vetoList, inputClusters, availableClusters);
238  this->SelectClusters(TPC_VIEW_U, availableClusters, clustersU);
239  this->SelectClusters(TPC_VIEW_V, availableClusters, clustersV);
240  this->SelectClusters(TPC_VIEW_W, availableClusters, clustersW);
241 
242  float chi2(m_twoViewChi2Cut);
243  const Cluster *pCluster1(NULL);
244  const Cluster *pCluster2(NULL);
245 
246  this->GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, pCluster1, pCluster2, chi2);
247  this->GetBestChi2(pVertex, slidingFitResultMap, clustersV, clustersW, pCluster1, pCluster2, chi2);
248  this->GetBestChi2(pVertex, slidingFitResultMap, clustersW, clustersU, pCluster1, pCluster2, chi2);
249 
250  if (NULL == pCluster1 || NULL == pCluster2)
251  return;
252 
253  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
254  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
255 
256  const Cluster *const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
257  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
258  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
259 
260  particleList.push_back(Particle(pClusterU, pClusterV, pClusterW));
261 
262  vetoList.insert(pCluster1);
263  vetoList.insert(pCluster2);
264  }
265 }
266 
267 //------------------------------------------------------------------------------------------------------------------------------------------
268 
269 void VertexBasedPfoRecoveryAlgorithm::GetBestChi2(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
270  const ClusterVector &clusters1, const ClusterVector &clusters2, const ClusterVector &clusters3, const Cluster *&pBestCluster1,
271  const Cluster *&pBestCluster2, const Cluster *&pBestCluster3, float &bestChi2) const
272 {
273  if (clusters1.empty() || clusters2.empty() || clusters3.empty())
274  return;
275 
276  // First loop
277  for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
278  {
279  const Cluster *const pCluster1 = *cIter1;
280 
281  TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
282  if (slidingFitResultMap.end() == sIter1)
283  continue;
284 
285  const TwoDSlidingFitResult &slidingFitResult1 = sIter1->second;
286  const LArPointingCluster pointingCluster1(slidingFitResult1);
287 
288  // Second loop
289  for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
290  {
291  const Cluster *const pCluster2 = *cIter2;
292 
293  TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
294  if (slidingFitResultMap.end() == sIter2)
295  continue;
296 
297  const TwoDSlidingFitResult &slidingFitResult2 = sIter2->second;
298  const LArPointingCluster pointingCluster2(slidingFitResult2);
299 
300  // Third loop
301  for (ClusterVector::const_iterator cIter3 = clusters3.begin(), cIterEnd3 = clusters3.end(); cIter3 != cIterEnd3; ++cIter3)
302  {
303  const Cluster *const pCluster3 = *cIter3;
304 
305  TwoDSlidingFitResultMap::const_iterator sIter3 = slidingFitResultMap.find(pCluster3);
306  if (slidingFitResultMap.end() == sIter3)
307  continue;
308 
309  const TwoDSlidingFitResult &slidingFitResult3 = sIter3->second;
310  const LArPointingCluster pointingCluster3(slidingFitResult3);
311 
312  // Calculate chi-squared
313  const float thisChi2(this->GetChi2(pVertex, pointingCluster1, pointingCluster2, pointingCluster3));
314 
315  if (thisChi2 < bestChi2)
316  {
317  bestChi2 = thisChi2;
318  pBestCluster1 = pCluster1;
319  pBestCluster2 = pCluster2;
320  pBestCluster3 = pCluster3;
321  }
322  }
323  }
324  }
325 }
326 
327 //------------------------------------------------------------------------------------------------------------------------------------------
328 
329 void VertexBasedPfoRecoveryAlgorithm::GetBestChi2(const Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap,
330  const ClusterVector &clusters1, const ClusterVector &clusters2, const Cluster *&pBestCluster1, const Cluster *&pBestCluster2, float &bestChi2) const
331 {
332  if (clusters1.empty() || clusters2.empty())
333  return;
334 
335  // First loop
336  for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
337  {
338  const Cluster *const pCluster1 = *cIter1;
339 
340  TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
341  if (slidingFitResultMap.end() == sIter1)
342  continue;
343 
344  const TwoDSlidingFitResult &slidingFitResult1 = sIter1->second;
345  const LArPointingCluster pointingCluster1(slidingFitResult1);
346 
347  // Second loop
348  for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
349  {
350  const Cluster *const pCluster2 = *cIter2;
351 
352  TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
353  if (slidingFitResultMap.end() == sIter2)
354  continue;
355 
356  const TwoDSlidingFitResult &slidingFitResult2 = sIter2->second;
357  const LArPointingCluster pointingCluster2(slidingFitResult2);
358 
359  // Calculate chi-squared
360  const float thisChi2(this->GetChi2(pVertex, pointingCluster1, pointingCluster2));
361 
362  if (thisChi2 < bestChi2)
363  {
364  bestChi2 = thisChi2;
365  pBestCluster1 = pCluster1;
366  pBestCluster2 = pCluster2;
367  }
368  }
369  }
370 }
371 
372 //------------------------------------------------------------------------------------------------------------------------------------------
373 
375  const Vertex *const pVertex, const LArPointingCluster &pointingCluster1, const LArPointingCluster &pointingCluster2) const
376 {
377  const HitType hitType1(LArClusterHelper::GetClusterHitType(pointingCluster1.GetCluster()));
378  const HitType hitType2(LArClusterHelper::GetClusterHitType(pointingCluster2.GetCluster()));
379 
380  if (hitType1 == hitType2)
381  throw StatusCodeException(STATUS_CODE_FAILURE);
382 
383  const CartesianVector vertex1(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType1));
384  const CartesianVector vertex2(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType2));
385 
386  const LArPointingCluster::Vertex &pointingVertex1(this->GetOuterVertex(vertex1, pointingCluster1));
387  const LArPointingCluster::Vertex &pointingVertex2(this->GetOuterVertex(vertex2, pointingCluster2));
388 
389  float chi2(0.f);
390  CartesianVector mergedPosition(0.f, 0.f, 0.f);
392  this->GetPandora(), hitType1, hitType2, pointingVertex1.GetPosition(), pointingVertex2.GetPosition(), mergedPosition, chi2);
393 
394  return chi2;
395 }
396 
397 //------------------------------------------------------------------------------------------------------------------------------------------
398 
399 float VertexBasedPfoRecoveryAlgorithm::GetChi2(const Vertex *const pVertex, const LArPointingCluster &pointingCluster1,
400  const LArPointingCluster &pointingCluster2, const LArPointingCluster &pointingCluster3) const
401 {
402  const HitType hitType1(LArClusterHelper::GetClusterHitType(pointingCluster1.GetCluster()));
403  const HitType hitType2(LArClusterHelper::GetClusterHitType(pointingCluster2.GetCluster()));
404  const HitType hitType3(LArClusterHelper::GetClusterHitType(pointingCluster3.GetCluster()));
405 
406  if ((hitType1 == hitType2) || (hitType2 == hitType3) || (hitType3 == hitType1))
407  throw StatusCodeException(STATUS_CODE_FAILURE);
408 
409  const CartesianVector vertex1(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType1));
410  const CartesianVector vertex2(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType2));
411  const CartesianVector vertex3(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex->GetPosition(), hitType3));
412 
413  const LArPointingCluster::Vertex &pointingVertex1(this->GetOuterVertex(vertex1, pointingCluster1));
414  const LArPointingCluster::Vertex &pointingVertex2(this->GetOuterVertex(vertex2, pointingCluster2));
415  const LArPointingCluster::Vertex &pointingVertex3(this->GetOuterVertex(vertex3, pointingCluster3));
416 
417  float chi2(0.f);
418  CartesianVector mergedPosition(0.f, 0.f, 0.f);
419  LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), hitType1, hitType2, hitType3, pointingVertex1.GetPosition(),
420  pointingVertex2.GetPosition(), pointingVertex3.GetPosition(), mergedPosition, chi2);
421 
422  return chi2;
423 }
424 
425 //------------------------------------------------------------------------------------------------------------------------------------------
426 
427 void VertexBasedPfoRecoveryAlgorithm::SelectAvailableClusters(const ClusterSet &vetoList, const ClusterVector &inputVector, ClusterVector &outputVector) const
428 {
429  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
430  {
431  if (0 == vetoList.count(*iter))
432  outputVector.push_back(*iter);
433  }
434 }
435 
436 //------------------------------------------------------------------------------------------------------------------------------------------
437 
438 void VertexBasedPfoRecoveryAlgorithm::SelectClusters(const HitType hitType, const ClusterVector &inputVector, ClusterVector &outputVector) const
439 {
440  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
441  {
442  if (hitType == LArClusterHelper::GetClusterHitType(*iter))
443  outputVector.push_back(*iter);
444  }
445 }
446 
447 //------------------------------------------------------------------------------------------------------------------------------------------
448 
450 {
451  const float innerDistance((vertex - cluster.GetInnerVertex().GetPosition()).GetMagnitudeSquared());
452  const float outerDistance((vertex - cluster.GetOuterVertex().GetPosition()).GetMagnitudeSquared());
453 
454  if (innerDistance < outerDistance)
455  return cluster.GetInnerVertex();
456  else
457  return cluster.GetOuterVertex();
458 }
459 
460 //------------------------------------------------------------------------------------------------------------------------------------------
461 
463 {
464  const LArPointingCluster::Vertex &innerVertex = this->GetInnerVertex(vertex, cluster);
465 
466  if (innerVertex.IsInnerVertex())
467  return cluster.GetOuterVertex();
468  else
469  return cluster.GetInnerVertex();
470 }
471 
472 //------------------------------------------------------------------------------------------------------------------------------------------
473 
475 {
476  if (particleList.empty())
477  return;
478 
479  const PfoList *pPfoList = NULL;
480  std::string pfoListName;
481  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
482 
483  for (ParticleList::const_iterator iter = particleList.begin(), iterEnd = particleList.end(); iter != iterEnd; ++iter)
484  {
485  const Particle &particle = *iter;
486 
487  ClusterList clusterList;
488  const Cluster *const pClusterU = particle.m_pClusterU;
489  const Cluster *const pClusterV = particle.m_pClusterV;
490  const Cluster *const pClusterW = particle.m_pClusterW;
491 
492  const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() : true);
493  const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() : true);
494  const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() : true);
495 
496  if (!(isAvailableU && isAvailableV && isAvailableW))
497  throw StatusCodeException(STATUS_CODE_FAILURE);
498 
499  if (pClusterU)
500  clusterList.push_back(pClusterU);
501  if (pClusterV)
502  clusterList.push_back(pClusterV);
503  if (pClusterW)
504  clusterList.push_back(pClusterW);
505 
506  // TODO Correct these placeholder parameters
507  PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
508  pfoParameters.m_particleId = MU_MINUS; // TRACK
509  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
510  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
511  pfoParameters.m_energy = 0.f;
512  pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
513  pfoParameters.m_clusterList = clusterList;
514 
515  const ParticleFlowObject *pPfo(NULL);
516  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
517  }
518 
519  if (!pPfoList->empty())
520  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
521 }
522 
523 //------------------------------------------------------------------------------------------------------------------------------------------
524 
525 VertexBasedPfoRecoveryAlgorithm::Particle::Particle(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW) :
526  m_pClusterU(pClusterU),
527  m_pClusterV(pClusterV),
528  m_pClusterW(pClusterW)
529 {
530  if (NULL == m_pClusterU && NULL == m_pClusterV && NULL == m_pClusterW)
531  throw StatusCodeException(STATUS_CODE_FAILURE);
532 
533  const HitType hitTypeU(NULL == m_pClusterU ? TPC_VIEW_U : LArClusterHelper::GetClusterHitType(m_pClusterU));
534  const HitType hitTypeV(NULL == m_pClusterV ? TPC_VIEW_V : LArClusterHelper::GetClusterHitType(m_pClusterV));
535  const HitType hitTypeW(NULL == m_pClusterW ? TPC_VIEW_W : LArClusterHelper::GetClusterHitType(m_pClusterW));
536 
537  if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
538  throw StatusCodeException(STATUS_CODE_FAILURE);
539 }
540 
541 //------------------------------------------------------------------------------------------------------------------------------------------
542 
543 StatusCode VertexBasedPfoRecoveryAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
544 {
545  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));
546 
547  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));
548 
549  PANDORA_RETURN_RESULT_IF_AND_IF(
550  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitHalfWindow", m_slidingFitHalfWindow));
551 
552  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
553  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
554 
555  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
556  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
557 
558  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "TwoViewChi2Cut", m_twoViewChi2Cut));
559 
560  PANDORA_RETURN_RESULT_IF_AND_IF(
561  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ThreeViewChi2Cut", m_threeViewChi2Cut));
562 
563  return STATUS_CODE_SUCCESS;
564 }
565 
566 } // namespace lar_content
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.
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.
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
Cluster finding and building.
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.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
TFile f
Definition: plotHisto.C:6
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 float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
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.
void SelectAvailableClusters(const pandora::ClusterSet &vetoList, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select cluster which haven&#39;t been vetoed.
const pandora::Cluster * m_pClusterV
Address of cluster in V view.
HitType
Definition: HitType.h:12
void SelectClusters(const pandora::HitType hitType, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select clusters of a specified hit type.
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
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:169
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
void BuildParticles(const ParticleList &particleList)
Build particle flow objects from matched clusters.
const pandora::Cluster * m_pClusterU
Address of cluster in U view.
vertex reconstruction