LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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((pVertexList && (pVertexList->size() == 1) && (VERTEX_3D == (*(pVertexList->begin()))->GetVertexType())) ? *(pVertexList->begin()) : NULL);
39 
40  if (!pSelectedVertex)
41  {
42  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
43  std::cout << "VertexBasedPfoRecoveryAlgorithm: unable to find vertex in current list " << std::endl;
44 
45  return STATUS_CODE_SUCCESS;
46  }
47 
48  // Get the available clusters from each view
49  ClusterVector availableClusters;
50  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNames, availableClusters));
51 
52  // Build a set of sliding fit results
53  TwoDSlidingFitResultMap slidingFitResultMap;
54  this->BuildSlidingFitResultMap(availableClusters, slidingFitResultMap);
55 
56  // Select seed clusters (adjacent to vertex)
57  ClusterVector selectedClusters;
58  this->SelectVertexClusters(pSelectedVertex, slidingFitResultMap, availableClusters, selectedClusters);
59 
60  // Match the cluster end points
61  ClusterSet vetoList;
62  ParticleList particleList;
63  this->MatchThreeViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
64  this->MatchTwoViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
65 
66  // Build new particles
67  this->BuildParticles(particleList);
68 
69  return STATUS_CODE_SUCCESS;
70 }
71 
72 //------------------------------------------------------------------------------------------------------------------------------------------
73 
74 StatusCode VertexBasedPfoRecoveryAlgorithm::GetAvailableClusters(const StringVector inputClusterListNames, ClusterVector &clusterVector) const
75 {
76  for (StringVector::const_iterator iter = inputClusterListNames.begin(), iterEnd = inputClusterListNames.end(); iter != iterEnd; ++iter)
77  {
78  const ClusterList *pClusterList(NULL);
79  PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, *iter, pClusterList));
80 
81  if (NULL == pClusterList)
82  {
83  if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
84  std::cout << "VertexBasedPfoRecoveryAlgorithm: could not find cluster list " << *iter << std::endl;
85  continue;
86  }
87 
88  for (ClusterList::const_iterator cIter = pClusterList->begin(), cIterEnd = pClusterList->end(); cIter != cIterEnd; ++cIter)
89  {
90  const Cluster *const pCluster = *cIter;
91 
92  if (!pCluster->IsAvailable())
93  continue;
94 
95  if (pCluster->GetNCaloHits() <=1)
96  continue;
97 
98  if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
99  continue;
100 
101  clusterVector.push_back(pCluster);
102  }
103  }
104 
105  std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
106 
107  return STATUS_CODE_SUCCESS;
108 }
109 //------------------------------------------------------------------------------------------------------------------------------------------
110 
112  TwoDSlidingFitResultMap &slidingFitResultMap) const
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((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 :
210  (TPC_VIEW_U == hitType3) ? pCluster3 : NULL);
211  const Cluster *const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 :
212  (TPC_VIEW_V == hitType3) ? pCluster3 : NULL);
213  const Cluster *const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 :
214  (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,
266  const Cluster *&pBestCluster1, 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 
369 float VertexBasedPfoRecoveryAlgorithm::GetChi2(const Vertex *const pVertex, const LArPointingCluster &pointingCluster1,
370  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);
386  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), hitType1, hitType2,
387  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,
415  pointingVertex1.GetPosition(), pointingVertex2.GetPosition(), pointingVertex3.GetPosition(), mergedPosition, chi2);
416 
417  return chi2;
418 }
419 
420 //------------------------------------------------------------------------------------------------------------------------------------------
421 
422 void VertexBasedPfoRecoveryAlgorithm::SelectAvailableClusters(const ClusterSet &vetoList, const ClusterVector &inputVector,
423  ClusterVector &outputVector) const
424 {
425  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
426  {
427  if (0 == vetoList.count(*iter))
428  outputVector.push_back(*iter);
429  }
430 }
431 
432 //------------------------------------------------------------------------------------------------------------------------------------------
433 
434 void VertexBasedPfoRecoveryAlgorithm::SelectClusters(const HitType hitType, const ClusterVector &inputVector, ClusterVector &outputVector) const
435 {
436  for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
437  {
438  if (hitType == LArClusterHelper::GetClusterHitType(*iter))
439  outputVector.push_back(*iter);
440  }
441 }
442 
443 //------------------------------------------------------------------------------------------------------------------------------------------
444 
446  const LArPointingCluster &cluster) const
447 {
448  const float innerDistance((vertex - cluster.GetInnerVertex().GetPosition()).GetMagnitudeSquared());
449  const float outerDistance((vertex - cluster.GetOuterVertex().GetPosition()).GetMagnitudeSquared());
450 
451  if (innerDistance < outerDistance)
452  return cluster.GetInnerVertex();
453  else
454  return cluster.GetOuterVertex();
455 }
456 
457 //------------------------------------------------------------------------------------------------------------------------------------------
458 
460  const LArPointingCluster &cluster) const
461 {
462  const LArPointingCluster::Vertex &innerVertex = this->GetInnerVertex(vertex, cluster);
463 
464  if (innerVertex.IsInnerVertex())
465  return cluster.GetOuterVertex();
466  else
467  return cluster.GetInnerVertex();
468 }
469 
470 //------------------------------------------------------------------------------------------------------------------------------------------
471 
473 {
474  if (particleList.empty())
475  return;
476 
477  const PfoList *pPfoList = NULL; std::string pfoListName;
478  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));
479 
480  for (ParticleList::const_iterator iter = particleList.begin(), iterEnd = particleList.end(); iter != iterEnd; ++iter)
481  {
482  const Particle &particle = *iter;
483 
484  ClusterList clusterList;
485  const Cluster *const pClusterU = particle.m_pClusterU;
486  const Cluster *const pClusterV = particle.m_pClusterV;
487  const Cluster *const pClusterW = particle.m_pClusterW;
488 
489  const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() : true);
490  const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() : true);
491  const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() : true);
492 
493  if(!(isAvailableU && isAvailableV && isAvailableW))
494  throw StatusCodeException(STATUS_CODE_FAILURE);
495 
496  if (pClusterU) clusterList.push_back(pClusterU);
497  if (pClusterV) clusterList.push_back(pClusterV);
498  if (pClusterW) clusterList.push_back(pClusterW);
499 
500  // TODO Correct these placeholder parameters
501  PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
502  pfoParameters.m_particleId = MU_MINUS; // TRACK
503  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
504  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
505  pfoParameters.m_energy = 0.f;
506  pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
507  pfoParameters.m_clusterList = clusterList;
508 
509  const ParticleFlowObject *pPfo(NULL);
510  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
511  }
512 
513  if (!pPfoList->empty())
514  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
515 }
516 
517 //------------------------------------------------------------------------------------------------------------------------------------------
518 
519 VertexBasedPfoRecoveryAlgorithm::Particle::Particle(const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW) :
520  m_pClusterU(pClusterU),
521  m_pClusterV(pClusterV),
522  m_pClusterW(pClusterW)
523 {
524  if (NULL == m_pClusterU && NULL == m_pClusterV && NULL == m_pClusterW)
525  throw StatusCodeException(STATUS_CODE_FAILURE);
526 
527  const HitType hitTypeU(NULL == m_pClusterU ? TPC_VIEW_U : LArClusterHelper::GetClusterHitType(m_pClusterU));
528  const HitType hitTypeV(NULL == m_pClusterV ? TPC_VIEW_V : LArClusterHelper::GetClusterHitType(m_pClusterV));
529  const HitType hitTypeW(NULL == m_pClusterW ? TPC_VIEW_W : LArClusterHelper::GetClusterHitType(m_pClusterW));
530 
531  if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
532  throw StatusCodeException(STATUS_CODE_FAILURE);
533 }
534 
535 //------------------------------------------------------------------------------------------------------------------------------------------
536 
537 StatusCode VertexBasedPfoRecoveryAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
538 {
539  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
540  "InputClusterListNames", m_inputClusterListNames));
541 
542  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
543  "OutputPfoListName", m_outputPfoListName));
544 
545  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
546  "SlidingFitHalfWindow", m_slidingFitHalfWindow));
547 
548  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
549  "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
550 
551  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
552  "MaxTransverseDisplacement", m_maxTransverseDisplacement));
553 
554  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
555  "TwoViewChi2Cut", m_twoViewChi2Cut));
556 
557  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
558  "ThreeViewChi2Cut", m_threeViewChi2Cut));
559 
560  return STATUS_CODE_SUCCESS;
561 }
562 
563 } // 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.
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.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
pandora::StringVector m_inputClusterListNames
The list of input cluster list names.
intermediate_table::const_iterator const_iterator
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::vector< art::Ptr< recob::Cluster > > ClusterVector
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.
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.
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:178
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