LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArPfoHelper.cc
Go to the documentation of this file.
1 
9 #include "Pandora/PdgTable.h"
10 
15 
17 
18 #include <algorithm>
19 #include <cmath>
20 #include <limits>
21 
22 using namespace pandora;
23 
24 namespace lar_content
25 {
26 
27 void LArPfoHelper::GetCoordinateVector(const ParticleFlowObject *const pPfo, const HitType &hitType, CartesianPointVector &coordinateVector)
28 {
29  ClusterList clusterList;
30  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
31 
32  for (const Cluster *const pCluster : clusterList)
33  LArClusterHelper::GetCoordinateVector(pCluster, coordinateVector);
34 }
35 
36 //------------------------------------------------------------------------------------------------------------------------------------------
37 
38 void LArPfoHelper::GetCaloHits(const PfoList &pfoList, const HitType &hitType, CaloHitList &caloHitList)
39 {
40  for (const ParticleFlowObject *const pPfo : pfoList)
41  LArPfoHelper::GetCaloHits(pPfo, hitType, caloHitList);
42 }
43 
44 //------------------------------------------------------------------------------------------------------------------------------------------
45 
46 void LArPfoHelper::GetCaloHits(const ParticleFlowObject *const pPfo, const HitType &hitType, CaloHitList &caloHitList)
47 {
48  ClusterList clusterList;
49  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
50 
51  for (const Cluster *const pCluster : clusterList)
52  pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
53 }
54 
55 //------------------------------------------------------------------------------------------------------------------------------------------
56 
57 void LArPfoHelper::GetIsolatedCaloHits(const PfoList &pfoList, const HitType &hitType, CaloHitList &caloHitList)
58 {
59  for (const ParticleFlowObject *const pPfo : pfoList)
60  LArPfoHelper::GetIsolatedCaloHits(pPfo, hitType, caloHitList);
61 }
62 
63 //------------------------------------------------------------------------------------------------------------------------------------------
64 
65 void LArPfoHelper::GetIsolatedCaloHits(const ParticleFlowObject *const pPfo, const HitType &hitType, CaloHitList &caloHitList)
66 {
67  ClusterList clusterList;
68  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
69 
70  for (const Cluster *const pCluster : clusterList)
71  caloHitList.insert(caloHitList.end(), pCluster->GetIsolatedCaloHitList().begin(), pCluster->GetIsolatedCaloHitList().end());
72 }
73 
74 //------------------------------------------------------------------------------------------------------------------------------------------
75 
76 void LArPfoHelper::GetAllCaloHits(const ParticleFlowObject *const pPfo, CaloHitList &caloHitList)
77 {
78  LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_U, caloHitList);
79  LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_V, caloHitList);
80  LArPfoHelper::GetCaloHits(pPfo, TPC_VIEW_W, caloHitList);
81  LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_U, caloHitList);
82  LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_V, caloHitList);
83  LArPfoHelper::GetIsolatedCaloHits(pPfo, TPC_VIEW_W, caloHitList);
84 }
85 
86 //------------------------------------------------------------------------------------------------------------------------------------------
87 
88 void LArPfoHelper::GetClusters(const PfoList &pfoList, const HitType &hitType, ClusterList &clusterList)
89 {
90  for (const ParticleFlowObject *const pPfo : pfoList)
91  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
92 }
93 
94 //------------------------------------------------------------------------------------------------------------------------------------------
95 
96 void LArPfoHelper::GetClusters(const ParticleFlowObject *const pPfo, const HitType &hitType, ClusterList &clusterList)
97 {
98  for (const Cluster *const pCluster : pPfo->GetClusterList())
99  {
100  if (hitType != LArClusterHelper::GetClusterHitType(pCluster))
101  continue;
102 
103  clusterList.push_back(pCluster);
104  }
105 }
106 
107 //------------------------------------------------------------------------------------------------------------------------------------------
108 
109 unsigned int LArPfoHelper::GetNumberOfTwoDHits(const ParticleFlowObject *const pPfo)
110 {
111  unsigned int totalHits(0);
112 
113  ClusterList clusterList;
114  LArPfoHelper::GetTwoDClusterList(pPfo, clusterList);
115  for (const Cluster *const pCluster : clusterList)
116  totalHits += pCluster->GetNCaloHits();
117 
118  return totalHits;
119 }
120 
121 //------------------------------------------------------------------------------------------------------------------------------------------
122 
123 void LArPfoHelper::GetTwoDClusterList(const ParticleFlowObject *const pPfo, ClusterList &clusterList)
124 {
125  for (const Cluster *const pCluster : pPfo->GetClusterList())
126  {
127  if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
128  continue;
129 
130  clusterList.push_back(pCluster);
131  }
132 }
133 
134 //------------------------------------------------------------------------------------------------------------------------------------------
135 
136 void LArPfoHelper::GetThreeDClusterList(const ParticleFlowObject *const pPfo, ClusterList &clusterList)
137 {
138  for (const Cluster *const pCluster : pPfo->GetClusterList())
139  {
140  if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
141  continue;
142 
143  clusterList.push_back(pCluster);
144  }
145 }
146 
147 //------------------------------------------------------------------------------------------------------------------------------------------
148 
149 void LArPfoHelper::GetAllConnectedPfos(const PfoList &inputPfoList, PfoList &outputPfoList)
150 {
151  for (const ParticleFlowObject *const pPfo : inputPfoList)
152  LArPfoHelper::GetAllConnectedPfos(pPfo, outputPfoList);
153 }
154 
155 //------------------------------------------------------------------------------------------------------------------------------------------
156 
157 void LArPfoHelper::GetAllConnectedPfos(const ParticleFlowObject *const pPfo, PfoList &outputPfoList)
158 {
159  if (outputPfoList.end() != std::find(outputPfoList.begin(), outputPfoList.end(), pPfo))
160  return;
161 
162  outputPfoList.push_back(pPfo);
163  LArPfoHelper::GetAllConnectedPfos(pPfo->GetParentPfoList(), outputPfoList);
164  LArPfoHelper::GetAllConnectedPfos(pPfo->GetDaughterPfoList(), outputPfoList);
165 }
166 
167 //------------------------------------------------------------------------------------------------------------------------------------------
168 
169 void LArPfoHelper::GetAllDownstreamPfos(const PfoList &inputPfoList, PfoList &outputPfoList)
170 {
171  for (const ParticleFlowObject *const pPfo : inputPfoList)
172  LArPfoHelper::GetAllDownstreamPfos(pPfo, outputPfoList);
173 }
174 
175 //------------------------------------------------------------------------------------------------------------------------------------------
176 
177 void LArPfoHelper::GetAllDownstreamPfos(const ParticleFlowObject *const pPfo, PfoList &outputPfoList)
178 {
179  if (outputPfoList.end() != std::find(outputPfoList.begin(), outputPfoList.end(), pPfo))
180  return;
181 
182  outputPfoList.push_back(pPfo);
183  LArPfoHelper::GetAllDownstreamPfos(pPfo->GetDaughterPfoList(), outputPfoList);
184 }
185 
186 //------------------------------------------------------------------------------------------------------------------------------------------
187 
188 void LArPfoHelper::GetAllDownstreamPfos(
189  const pandora::ParticleFlowObject *const pPfo, pandora::PfoList &outputTrackPfoList, pandora::PfoList &outputLeadingShowerPfoList)
190 {
191  if (LArPfoHelper::IsTrack(pPfo))
192  {
193  outputTrackPfoList.emplace_back(pPfo);
194  for (const ParticleFlowObject *pChild : pPfo->GetDaughterPfoList())
195  {
196  if (std::find(outputTrackPfoList.begin(), outputTrackPfoList.end(), pChild) == outputTrackPfoList.end())
197  {
198  const int pdg{std::abs(pChild->GetParticleId())};
199  if (pdg == E_MINUS)
200  {
201  outputLeadingShowerPfoList.emplace_back(pChild);
202  }
203  else
204  {
205  outputTrackPfoList.emplace_back(pChild);
206  LArPfoHelper::GetAllDownstreamPfos(pChild, outputTrackPfoList, outputLeadingShowerPfoList);
207  }
208  }
209  }
210  }
211  else
212  {
213  outputLeadingShowerPfoList.emplace_back(pPfo);
214  }
215 }
216 
217 //------------------------------------------------------------------------------------------------------------------------------------------
218 
219 int LArPfoHelper::GetHierarchyTier(const ParticleFlowObject *const pPfo)
220 {
221  const ParticleFlowObject *pParentPfo = pPfo;
222  int tier(0);
223 
224  while (pParentPfo->GetParentPfoList().empty() == false)
225  {
226  if (1 != pParentPfo->GetParentPfoList().size())
227  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
228 
229  pParentPfo = *(pParentPfo->GetParentPfoList().begin());
230  ++tier;
231  }
232 
233  return tier;
234 }
235 
236 //------------------------------------------------------------------------------------------------------------------------------------------
237 
238 float LArPfoHelper::GetTwoDLengthSquared(const ParticleFlowObject *const pPfo)
239 {
240  if (!LArPfoHelper::IsTwoD(pPfo))
241  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
242 
243  float lengthSquared(0.f);
244 
245  for (const Cluster *const pCluster : pPfo->GetClusterList())
246  {
247  if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
248  continue;
249 
250  lengthSquared += LArClusterHelper::GetLengthSquared(pCluster);
251  }
252 
253  return lengthSquared;
254 }
255 
256 //------------------------------------------------------------------------------------------------------------------------------------------
257 
258 float LArPfoHelper::GetThreeDLengthSquared(const ParticleFlowObject *const pPfo)
259 {
260  if (!LArPfoHelper::IsThreeD(pPfo))
261  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
262 
263  float lengthSquared(0.f);
264 
265  for (const Cluster *const pCluster : pPfo->GetClusterList())
266  {
267  if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
268  continue;
269 
270  lengthSquared += LArClusterHelper::GetLengthSquared(pCluster);
271  }
272 
273  return lengthSquared;
274 }
275 
276 //------------------------------------------------------------------------------------------------------------------------------------------
277 
278 float LArPfoHelper::GetClosestDistance(const ParticleFlowObject *const pPfo, const Cluster *const pCluster)
279 {
280  const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
281 
282  ClusterList clusterList;
283  LArPfoHelper::GetClusters(pPfo, hitType, clusterList);
284 
285  if (clusterList.empty())
286  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
287 
288  float bestDistance(std::numeric_limits<float>::max());
289 
290  for (const Cluster *const pPfoCluster : clusterList)
291  {
292  const float thisDistance(LArClusterHelper::GetClosestDistance(pCluster, pPfoCluster));
293 
294  if (thisDistance < bestDistance)
295  bestDistance = thisDistance;
296  }
297 
298  return bestDistance;
299 }
300 
301 //------------------------------------------------------------------------------------------------------------------------------------------
302 
303 float LArPfoHelper::GetThreeDSeparation(const ParticleFlowObject *const pPfo1, const ParticleFlowObject *const pPfo2)
304 {
305  ClusterList clusterList1, clusterList2;
306 
307  LArPfoHelper::GetClusters(pPfo1, TPC_3D, clusterList1);
308  LArPfoHelper::GetClusters(pPfo2, TPC_3D, clusterList2);
309 
310  if (clusterList1.empty() || clusterList2.empty())
311  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
312 
313  return LArClusterHelper::GetClosestDistance(clusterList1, clusterList2);
314 }
315 
316 //------------------------------------------------------------------------------------------------------------------------------------------
317 
318 bool LArPfoHelper::IsTwoD(const ParticleFlowObject *const pPfo)
319 {
320  for (const Cluster *const pCluster : pPfo->GetClusterList())
321  {
322  if (TPC_3D != LArClusterHelper::GetClusterHitType(pCluster))
323  return true;
324  }
325 
326  return false;
327 }
328 
329 //------------------------------------------------------------------------------------------------------------------------------------------
330 
331 bool LArPfoHelper::IsThreeD(const ParticleFlowObject *const pPfo)
332 {
333  for (const Cluster *const pCluster : pPfo->GetClusterList())
334  {
335  if (TPC_3D == LArClusterHelper::GetClusterHitType(pCluster))
336  return true;
337  }
338 
339  return false;
340 }
341 
342 //------------------------------------------------------------------------------------------------------------------------------------------
343 
344 bool LArPfoHelper::IsTrack(const ParticleFlowObject *const pPfo)
345 {
346  const int pdg(pPfo->GetParticleId());
347 
348  // muon, pion, proton, kaon
349  return ((MU_MINUS == std::abs(pdg)) || (PI_PLUS == std::abs(pdg)) || (PROTON == std::abs(pdg)) || (K_PLUS == std::abs(pdg)));
350 }
351 
352 //------------------------------------------------------------------------------------------------------------------------------------------
353 
354 bool LArPfoHelper::IsShower(const ParticleFlowObject *const pPfo)
355 {
356  const int pdg(pPfo->GetParticleId());
357 
358  // electron, photon
359  return ((E_MINUS == std::abs(pdg)) || (PHOTON == std::abs(pdg)));
360 }
361 
362 //------------------------------------------------------------------------------------------------------------------------------------------
363 
364 int LArPfoHelper::GetPrimaryNeutrino(const ParticleFlowObject *const pPfo)
365 {
366  try
367  {
368  const ParticleFlowObject *const pParentPfo = LArPfoHelper::GetParentNeutrino(pPfo);
369  return pParentPfo->GetParticleId();
370  }
371  catch (const StatusCodeException &)
372  {
373  return 0;
374  }
375 }
376 
377 //------------------------------------------------------------------------------------------------------------------------------------------
378 
379 bool LArPfoHelper::IsFinalState(const ParticleFlowObject *const pPfo)
380 {
381  if (pPfo->GetParentPfoList().empty() && !LArPfoHelper::IsNeutrino(pPfo))
382  return true;
383 
384  if (LArPfoHelper::IsNeutrinoFinalState(pPfo))
385  return true;
386 
387  return false;
388 }
389 
390 //------------------------------------------------------------------------------------------------------------------------------------------
391 
392 bool LArPfoHelper::IsNeutrinoFinalState(const ParticleFlowObject *const pPfo)
393 {
394  return ((pPfo->GetParentPfoList().size() == 1) && (LArPfoHelper::IsNeutrino(*(pPfo->GetParentPfoList().begin()))));
395 }
396 
397 //------------------------------------------------------------------------------------------------------------------------------------------
398 
399 bool LArPfoHelper::IsNeutrino(const ParticleFlowObject *const pPfo)
400 {
401  const int absoluteParticleId(std::abs(pPfo->GetParticleId()));
402 
403  if ((NU_E == absoluteParticleId) || (NU_MU == absoluteParticleId) || (NU_TAU == absoluteParticleId))
404  return true;
405 
406  return false;
407 }
408 
409 //------------------------------------------------------------------------------------------------------------------------------------------
410 
411 bool LArPfoHelper::IsTestBeamFinalState(const ParticleFlowObject *const pPfo)
412 {
413  return LArPfoHelper::IsTestBeam(LArPfoHelper::GetParentPfo(pPfo));
414 }
415 
416 //------------------------------------------------------------------------------------------------------------------------------------------
417 
418 bool LArPfoHelper::IsTestBeam(const ParticleFlowObject *const pPfo)
419 {
420  const PropertiesMap &properties(pPfo->GetPropertiesMap());
421  const PropertiesMap::const_iterator iter(properties.find("IsTestBeam"));
422 
423  if (iter != properties.end())
424  return ((iter->second > 0.f) ? true : false);
425 
426  return false;
427 }
428 
429 //------------------------------------------------------------------------------------------------------------------------------------------
430 
431 void LArPfoHelper::GetRecoNeutrinos(const PfoList *const pPfoList, PfoList &recoNeutrinos)
432 {
433  if (!pPfoList)
434  return;
435 
436  for (const ParticleFlowObject *const pPfo : *pPfoList)
437  {
438  if (LArPfoHelper::IsNeutrino(pPfo))
439  recoNeutrinos.push_back(pPfo);
440  }
441 }
442 
443 //------------------------------------------------------------------------------------------------------------------------------------------
444 
445 const ParticleFlowObject *LArPfoHelper::GetParentPfo(const ParticleFlowObject *const pPfo)
446 {
447  const ParticleFlowObject *pParentPfo = pPfo;
448 
449  while (pParentPfo->GetParentPfoList().empty() == false)
450  {
451  pParentPfo = *(pParentPfo->GetParentPfoList().begin());
452  }
453 
454  return pParentPfo;
455 }
456 
457 //------------------------------------------------------------------------------------------------------------------------------------------
458 
459 const ParticleFlowObject *LArPfoHelper::GetParentNeutrino(const ParticleFlowObject *const pPfo)
460 {
461  const ParticleFlowObject *const pParentPfo = LArPfoHelper::GetParentPfo(pPfo);
462 
463  if (!LArPfoHelper::IsNeutrino(pParentPfo))
464  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
465 
466  return pParentPfo;
467 }
468 
469 //------------------------------------------------------------------------------------------------------------------------------------------
470 
471 const Vertex *LArPfoHelper::GetVertex(const ParticleFlowObject *const pPfo)
472 {
473  if (pPfo->GetVertexList().empty())
474  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
475 
476  const Vertex *pVertex(nullptr);
477 
478  // ATTN : Test beam parent pfos contain an interaction and start vertex
479  if (LArPfoHelper::IsTestBeam(pPfo) && pPfo->GetParentPfoList().empty())
480  {
481  pVertex = LArPfoHelper::GetVertexWithLabel(pPfo->GetVertexList(), VERTEX_START);
482  }
483  else
484  {
485  if (pPfo->GetVertexList().size() != 1)
486  throw StatusCodeException(STATUS_CODE_FAILURE);
487 
488  pVertex = *(pPfo->GetVertexList().begin());
489  }
490 
491  if (VERTEX_3D != pVertex->GetVertexType())
492  throw StatusCodeException(STATUS_CODE_FAILURE);
493 
494  return pVertex;
495 }
496 
497 //------------------------------------------------------------------------------------------------------------------------------------------
498 
499 const Vertex *LArPfoHelper::GetTestBeamInteractionVertex(const ParticleFlowObject *const pPfo)
500 {
501  if (pPfo->GetVertexList().empty() || !pPfo->GetParentPfoList().empty() || !LArPfoHelper::IsTestBeam(pPfo))
502  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
503 
504  const Vertex *pInteractionVertex(LArPfoHelper::GetVertexWithLabel(pPfo->GetVertexList(), VERTEX_INTERACTION));
505 
506  if (VERTEX_3D != pInteractionVertex->GetVertexType())
507  throw StatusCodeException(STATUS_CODE_FAILURE);
508 
509  return pInteractionVertex;
510 }
511 
512 //------------------------------------------------------------------------------------------------------------------------------------------
513 
514 const Vertex *LArPfoHelper::GetVertexWithLabel(const VertexList &vertexList, const VertexLabel vertexLabel)
515 {
516  const Vertex *pTargetVertex(nullptr);
517 
518  for (const Vertex *pCandidateVertex : vertexList)
519  {
520  if (pCandidateVertex->GetVertexLabel() == vertexLabel)
521  {
522  if (!pTargetVertex)
523  {
524  pTargetVertex = pCandidateVertex;
525  }
526  else
527  {
528  throw StatusCodeException(STATUS_CODE_FAILURE);
529  }
530  }
531  }
532 
533  if (!pTargetVertex)
534  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
535 
536  return pTargetVertex;
537 }
538 
539 //------------------------------------------------------------------------------------------------------------------------------------------
540 
541 void LArPfoHelper::GetSlidingFitTrajectory(const CartesianPointVector &pointVector, const CartesianVector &vertexPosition,
542  const unsigned int layerWindow, const float layerPitch, LArTrackStateVector &trackStateVector, IntVector *const pIndexVector)
543 {
544  LArPfoHelper::SlidingFitTrajectoryImpl(&pointVector, vertexPosition, layerWindow, layerPitch, trackStateVector, pIndexVector);
545 }
546 
547 //------------------------------------------------------------------------------------------------------------------------------------------
548 
549 void LArPfoHelper::GetSlidingFitTrajectory(const ParticleFlowObject *const pPfo, const Vertex *const pVertex,
550  const unsigned int layerWindow, const float layerPitch, LArTrackStateVector &trackStateVector)
551 {
552  CaloHitList caloHitList;
553  LArPfoHelper::GetCaloHits(pPfo, TPC_3D, caloHitList);
554  LArPfoHelper::SlidingFitTrajectoryImpl(&caloHitList, pVertex->GetPosition(), layerWindow, layerPitch, trackStateVector);
555 }
556 
557 //------------------------------------------------------------------------------------------------------------------------------------------
558 
559 LArShowerPCA LArPfoHelper::GetPrincipalComponents(const CartesianPointVector &pointVector, const CartesianVector &vertexPosition)
560 {
561  // Run the PCA analysis
562  CartesianVector centroid(0.f, 0.f, 0.f);
563  LArPcaHelper::EigenVectors eigenVecs;
564  LArPcaHelper::EigenValues eigenValues(0.f, 0.f, 0.f);
565  LArPcaHelper::RunPca(pointVector, centroid, eigenValues, eigenVecs);
566 
567  // Require that principal eigenvalue should always be positive
568  if (eigenValues.GetX() < std::numeric_limits<float>::epsilon())
569  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
570 
571  // By convention, principal axis should always point away from vertex
572  const float testProjection(eigenVecs.at(0).GetDotProduct(vertexPosition - centroid));
573  const float directionScaleFactor((testProjection > std::numeric_limits<float>::epsilon()) ? -1.f : 1.f);
574 
575  const CartesianVector primaryAxis(eigenVecs.at(0) * directionScaleFactor);
576  const CartesianVector secondaryAxis(eigenVecs.at(1) * directionScaleFactor);
577  const CartesianVector tertiaryAxis(eigenVecs.at(2) * directionScaleFactor);
578 
579  return LArShowerPCA(centroid, primaryAxis, secondaryAxis, tertiaryAxis, eigenValues);
580 }
581 
582 //------------------------------------------------------------------------------------------------------------------------------------------
583 
584 LArShowerPCA LArPfoHelper::GetPrincipalComponents(const ParticleFlowObject *const pPfo, const Vertex *const pVertex)
585 {
586  CartesianPointVector pointVector;
587  LArPfoHelper::GetCoordinateVector(pPfo, TPC_3D, pointVector);
588  return LArPfoHelper::GetPrincipalComponents(pointVector, pVertex->GetPosition());
589 }
590 
591 //------------------------------------------------------------------------------------------------------------------------------------------
592 
593 bool LArPfoHelper::SortByHitProjection(const LArTrackTrajectoryPoint &lhs, const LArTrackTrajectoryPoint &rhs)
594 {
595  if (lhs.first != rhs.first)
596  return (lhs.first < rhs.first);
597 
598  // ATTN Removed to support use with CartesianVector only (no CaloHit) input
599  // if (lhs.second.GetCaloHit() && rhs.second.GetCaloHit())
600  // return (lhs.second.GetCaloHit()->GetInputEnergy() > rhs.second.GetCaloHit()->GetInputEnergy());
601 
602  const float dx(lhs.second.GetPosition().GetX() - rhs.second.GetPosition().GetX());
603  const float dy(lhs.second.GetPosition().GetY() - rhs.second.GetPosition().GetY());
604  const float dz(lhs.second.GetPosition().GetZ() - rhs.second.GetPosition().GetZ());
605  return (dx + dy + dz > 0.f);
606 }
607 
608 //------------------------------------------------------------------------------------------------------------------------------------------
609 
610 bool LArPfoHelper::SortByNHits(const ParticleFlowObject *const pLhs, const ParticleFlowObject *const pRhs)
611 {
612  unsigned int nTwoDHitsLhs(0), nThreeDHitsLhs(0);
613  float energyLhs(0.f);
614  for (ClusterList::const_iterator iter = pLhs->GetClusterList().begin(), iterEnd = pLhs->GetClusterList().end(); iter != iterEnd; ++iter)
615  {
616  const Cluster *const pClusterLhs = *iter;
617 
618  if (TPC_3D != LArClusterHelper::GetClusterHitType(pClusterLhs))
619  nTwoDHitsLhs += pClusterLhs->GetNCaloHits();
620  else
621  nThreeDHitsLhs += pClusterLhs->GetNCaloHits();
622 
623  energyLhs += pClusterLhs->GetHadronicEnergy();
624  }
625 
626  unsigned int nTwoDHitsRhs(0), nThreeDHitsRhs(0);
627  float energyRhs(0.f);
628  for (ClusterList::const_iterator iter = pRhs->GetClusterList().begin(), iterEnd = pRhs->GetClusterList().end(); iter != iterEnd; ++iter)
629  {
630  const Cluster *const pClusterRhs = *iter;
631 
632  if (TPC_3D != LArClusterHelper::GetClusterHitType(pClusterRhs))
633  nTwoDHitsRhs += pClusterRhs->GetNCaloHits();
634  else
635  nThreeDHitsRhs += pClusterRhs->GetNCaloHits();
636 
637  energyRhs += pClusterRhs->GetHadronicEnergy();
638  }
639 
640  if (nTwoDHitsLhs != nTwoDHitsRhs)
641  return (nTwoDHitsLhs > nTwoDHitsRhs);
642 
643  if (nThreeDHitsLhs != nThreeDHitsRhs)
644  return (nThreeDHitsLhs > nThreeDHitsRhs);
645 
646  // ATTN Need an efficient (balance with well-motivated) tie-breaker here. Pfo length, for instance, is extremely slow.
647  return (energyLhs > energyRhs);
648 }
649 
650 //------------------------------------------------------------------------------------------------------------------------------------------
651 
652 void LArPfoHelper::GetBreadthFirstHierarchyRepresentation(const pandora::ParticleFlowObject *const pPfo, pandora::PfoList &pfoList)
653 {
654  const ParticleFlowObject *pRoot{pPfo};
655  PfoList parents{pRoot->GetParentPfoList()};
656  while (!parents.empty())
657  {
658  if (parents.size() > 1)
659  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
660  pRoot = parents.front();
661  parents = pRoot->GetParentPfoList();
662  }
663  PfoList queue;
664  pfoList.emplace_back(pRoot);
665  queue.emplace_back(pRoot);
666 
667  while (!queue.empty())
668  {
669  const PfoList &daughters{queue.front()->GetDaughterPfoList()};
670  queue.pop_front();
671  for (const ParticleFlowObject *pDaughter : daughters)
672  {
673  pfoList.emplace_back(pDaughter);
674  queue.emplace_back(pDaughter);
675  }
676  }
677 }
678 
679 //------------------------------------------------------------------------------------------------------------------------------------------
680 
681 template <typename T>
682 void LArPfoHelper::SlidingFitTrajectoryImpl(const T *const pT, const CartesianVector &vertexPosition, const unsigned int layerWindow,
683  const float layerPitch, LArTrackStateVector &trackStateVector, IntVector *const pIndexVector)
684 {
685  CartesianPointVector pointVector;
686 
687  for (const auto &nextPoint : *pT)
688  pointVector.push_back(LArObjectHelper::TypeAdaptor::GetPosition(nextPoint));
689 
690  if (pointVector.empty())
691  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
692 
693  std::sort(pointVector.begin(), pointVector.end(), LArClusterHelper::SortCoordinatesByPosition);
694 
695  LArTrackTrajectory trackTrajectory;
696  IntVector indicesWithoutSpacePoints;
697  if (pIndexVector)
698  pIndexVector->clear();
699 
700  try
701  {
702  const ThreeDSlidingFitResult slidingFitResult(&pointVector, layerWindow, layerPitch);
703  const CartesianVector minPosition(slidingFitResult.GetGlobalMinLayerPosition());
704  const CartesianVector maxPosition(slidingFitResult.GetGlobalMaxLayerPosition());
705 
706  if ((maxPosition - minPosition).GetMagnitudeSquared() < std::numeric_limits<float>::epsilon())
707  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
708 
709  const CartesianVector seedPosition((maxPosition + minPosition) * 0.5f);
710  const CartesianVector seedDirection((maxPosition - minPosition).GetUnitVector());
711 
712  const float scaleFactor((seedDirection.GetDotProduct(seedPosition - vertexPosition) > 0.f) ? +1.f : -1.f);
713 
714  int index(-1);
715  for (const auto &nextPoint : *pT)
716  {
717  ++index;
718 
719  try
720  {
721  const float rL(slidingFitResult.GetLongitudinalDisplacement(LArObjectHelper::TypeAdaptor::GetPosition(nextPoint)));
722 
723  CartesianVector position(0.f, 0.f, 0.f);
724  const StatusCode positionStatusCode(slidingFitResult.GetGlobalFitPosition(rL, position));
725 
726  if (positionStatusCode != STATUS_CODE_SUCCESS)
727  throw StatusCodeException(positionStatusCode);
728 
729  CartesianVector direction(0.f, 0.f, 0.f);
730  const StatusCode directionStatusCode(slidingFitResult.GetGlobalFitDirection(rL, direction));
731 
732  if (directionStatusCode != STATUS_CODE_SUCCESS)
733  throw StatusCodeException(directionStatusCode);
734 
735  const float projection(seedDirection.GetDotProduct(position - seedPosition));
736 
737  trackTrajectory.push_back(LArTrackTrajectoryPoint(projection * scaleFactor,
738  LArTrackState(position, direction * scaleFactor, LArObjectHelper::TypeAdaptor::GetCaloHit(nextPoint)), index));
739  }
740  catch (const StatusCodeException &statusCodeException1)
741  {
742  indicesWithoutSpacePoints.push_back(index);
743 
744  if (STATUS_CODE_FAILURE == statusCodeException1.GetStatusCode())
745  throw statusCodeException1;
746  }
747  }
748  }
749  catch (const StatusCodeException &statusCodeException2)
750  {
751  if (STATUS_CODE_FAILURE == statusCodeException2.GetStatusCode())
752  throw statusCodeException2;
753  }
754 
755  // Sort trajectory points by distance along track
756  std::sort(trackTrajectory.begin(), trackTrajectory.end(), LArPfoHelper::SortByHitProjection);
757 
758  for (const LArTrackTrajectoryPoint &larTrackTrajectoryPoint : trackTrajectory)
759  {
760  trackStateVector.push_back(larTrackTrajectoryPoint.second);
761  if (pIndexVector)
762  pIndexVector->push_back(larTrackTrajectoryPoint.GetIndex());
763  }
764 
765  // Store indices of spacepoints with no associated trajectory point at the end of pIndexVector
766  if (pIndexVector)
767  {
768  for (const int index : indicesWithoutSpacePoints)
769  pIndexVector->push_back(index);
770  }
771 }
772 
773 //------------------------------------------------------------------------------------------------------------------------------------------
774 
775 template void LArPfoHelper::SlidingFitTrajectoryImpl(
776  const CartesianPointVector *const, const CartesianVector &, const unsigned int, const float, LArTrackStateVector &, IntVector *const);
777 template void LArPfoHelper::SlidingFitTrajectoryImpl(
778  const CaloHitList *const, const CartesianVector &, const unsigned int, const float, LArTrackStateVector &, IntVector *const);
779 
780 } // namespace lar_content
Header file for the pfo helper class.
pandora::CartesianVector EigenValues
Definition: LArPcaHelper.h:24
LArTrackState class.
Definition: LArPfoObjects.h:29
constexpr auto abs(T v)
Returns the absolute value of the argument.
Header file for the principal curve analysis helper class.
intermediate_table::const_iterator const_iterator
std::vector< int > IntVector
float GetLongitudinalDisplacement(const pandora::CartesianVector &position) const
Get longitudinal projection onto primary axis.
std::vector< LArTrackTrajectoryPoint > LArTrackTrajectory
TFile f
Definition: plotHisto.C:6
Header file for the object helper class.
Header file for the cluster helper class.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
const pandora::CartesianVector & GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
Header file for the lar three dimensional sliding fit result class.
std::vector< pandora::CartesianVector > EigenVectors
Definition: LArPcaHelper.h:25
HitType
Definition: HitType.h:12
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
std::vector< LArTrackState > LArTrackStateVector
Definition: LArPfoObjects.h:67
LArTrackTrajectoryPoint class.
Definition: LArPfoObjects.h:74
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
LArShowerPCA class.
std::list< Vertex > VertexList
Definition: DCEL.h:169
const pandora::CartesianVector & GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.