LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRayTrackMatchingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 CosmicRayTrackMatchingAlgorithm::CosmicRayTrackMatchingAlgorithm() :
22  m_clusterMinLength(10.f),
23  m_vtxXOverlap(3.f),
24  m_minXOverlap(3.f),
25  m_minXOverlapFraction(0.8f),
26  m_maxDisplacement(10.f)
27 {
28 }
29 
30 //------------------------------------------------------------------------------------------------------------------------------------------
31 
33 {
34  ClusterVector clusterVector;
35 
36  // Select long clusters
37  for (const Cluster *const pCluster : inputVector)
38  {
40  continue;
41 
42  clusterVector.push_back(pCluster);
43  }
44 
45  // Remove long delta rays
46  for (const Cluster *const pCluster : clusterVector)
47  {
48  const float lengthSquared(LArClusterHelper::GetLengthSquared(pCluster));
49  CartesianVector innerVertex(0.f, 0.f, 0.f), outerVertex(0.f, 0.f, 0.f);
50  LArClusterHelper::GetExtremalCoordinates(pCluster, innerVertex, outerVertex);
51 
52  bool isDeltaRay(false);
53 
54  for (const Cluster *const pClusterCheck : clusterVector)
55  {
56  if (pCluster == pClusterCheck)
57  continue;
58 
59  if ((LArClusterHelper::GetLengthSquared(pClusterCheck) > 10.f * lengthSquared) &&
60  (LArClusterHelper::GetClosestDistance(innerVertex, pClusterCheck) < m_vtxXOverlap ||
61  LArClusterHelper::GetClosestDistance(outerVertex, pClusterCheck) < m_vtxXOverlap))
62  {
63  isDeltaRay = true;
64  break;
65  }
66  }
67 
68  if (!isDeltaRay)
69  outputVector.push_back(pCluster);
70  }
71 }
72 
73 //------------------------------------------------------------------------------------------------------------------------------------------
74 
75 bool CosmicRayTrackMatchingAlgorithm::MatchClusters(const Cluster *const pCluster1, const Cluster *const pCluster2) const
76 {
77  // Use start and end points
78  CartesianVector innerVertex1(0.f, 0.f, 0.f), outerVertex1(0.f, 0.f, 0.f);
79  CartesianVector innerVertex2(0.f, 0.f, 0.f), outerVertex2(0.f, 0.f, 0.f);
80 
81  LArClusterHelper::GetExtremalCoordinates(pCluster1, innerVertex1, outerVertex1);
82  LArClusterHelper::GetExtremalCoordinates(pCluster2, innerVertex2, outerVertex2);
83 
84  const float dxA(std::fabs(innerVertex2.GetX() - innerVertex1.GetX())); // inner1 <-> inner2
85  const float dxB(std::fabs(outerVertex2.GetX() - outerVertex1.GetX())); // outer1 <-> outer2
86 
87  const float dxC(std::fabs(outerVertex2.GetX() - innerVertex1.GetX())); // inner1 <-> outer2
88  const float dxD(std::fabs(innerVertex2.GetX() - outerVertex1.GetX())); // inner2 <-> outer1
89 
90  const float xVertex(std::min(std::max(dxA, dxB), std::max(dxC, dxD)));
91 
92  if (xVertex < m_vtxXOverlap)
93  return true;
94 
95  // use X overlap
96  float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
97  pCluster1->GetClusterSpanX(xMin1, xMax1);
98  pCluster2->GetClusterSpanX(xMin2, xMax2);
99 
100  const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
101  const float xSpan(std::max(xMax1, xMax2) - std::min(xMin1, xMin2));
102 
103  if (xSpan < std::numeric_limits<float>::epsilon())
104  return false;
105 
106  if (xOverlap > m_minXOverlap && xOverlap / xSpan > m_minXOverlapFraction)
107  return true;
108 
109  return false;
110 }
111 
112 //------------------------------------------------------------------------------------------------------------------------------------------
113 
115  const Cluster *const pCluster1, const Cluster *const pCluster2, const Cluster *const pCluster3) const
116 {
117  // Check that three clusters have a consistent 3D position
118  const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
119  const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
120  const HitType hitType3(LArClusterHelper::GetClusterHitType(pCluster3));
121 
122  if (hitType1 == hitType2 || hitType2 == hitType3 || hitType3 == hitType1)
123  throw StatusCodeException(STATUS_CODE_FAILURE);
124 
125  CartesianVector innerVertex1(0.f, 0.f, 0.f), outerVertex1(0.f, 0.f, 0.f);
126  CartesianVector innerVertex2(0.f, 0.f, 0.f), outerVertex2(0.f, 0.f, 0.f);
127  CartesianVector innerVertex3(0.f, 0.f, 0.f), outerVertex3(0.f, 0.f, 0.f);
128 
129  LArClusterHelper::GetExtremalCoordinates(pCluster1, innerVertex1, outerVertex1);
130  LArClusterHelper::GetExtremalCoordinates(pCluster2, innerVertex2, outerVertex2);
131  LArClusterHelper::GetExtremalCoordinates(pCluster3, innerVertex3, outerVertex3);
132 
133  for (unsigned int n = 0; n < 4; ++n)
134  {
135  CartesianVector vtx1(1 == n ? outerVertex1 : innerVertex1);
136  CartesianVector end1(1 == n ? innerVertex1 : outerVertex1);
137 
138  CartesianVector vtx2(2 == n ? outerVertex2 : innerVertex2);
139  CartesianVector end2(2 == n ? innerVertex2 : outerVertex2);
140 
141  CartesianVector vtx3(3 == n ? outerVertex3 : innerVertex3);
142  CartesianVector end3(3 == n ? innerVertex3 : outerVertex3);
143 
144  if (std::fabs(vtx1.GetX() - vtx2.GetX()) < std::max(m_vtxXOverlap, std::fabs(vtx1.GetX() - end2.GetX())) &&
145  std::fabs(end1.GetX() - end2.GetX()) < std::max(m_vtxXOverlap, std::fabs(end1.GetX() - vtx2.GetX())) &&
146  std::fabs(vtx2.GetX() - vtx3.GetX()) < std::max(m_vtxXOverlap, std::fabs(vtx2.GetX() - end3.GetX())) &&
147  std::fabs(end2.GetX() - end3.GetX()) < std::max(m_vtxXOverlap, std::fabs(end2.GetX() - vtx3.GetX())) &&
148  std::fabs(vtx3.GetX() - vtx1.GetX()) < std::max(m_vtxXOverlap, std::fabs(vtx3.GetX() - end1.GetX())) &&
149  std::fabs(end3.GetX() - end1.GetX()) < std::max(m_vtxXOverlap, std::fabs(end3.GetX() - vtx1.GetX())))
150  {
151  float chi2(0.f);
152  CartesianVector projVtx1(0.f, 0.f, 0.f), projEnd1(0.f, 0.f, 0.f);
153  CartesianVector projVtx2(0.f, 0.f, 0.f), projEnd2(0.f, 0.f, 0.f);
154  CartesianVector projVtx3(0.f, 0.f, 0.f), projEnd3(0.f, 0.f, 0.f);
155 
156  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType1, hitType2, vtx1, vtx2, projVtx3, chi2);
157  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType1, hitType2, end1, end2, projEnd3, chi2);
158  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType2, hitType3, vtx2, vtx3, projVtx1, chi2);
159  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType2, hitType3, end2, end3, projEnd1, chi2);
160  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType3, hitType1, vtx3, vtx1, projVtx2, chi2);
161  LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType3, hitType1, end3, end1, projEnd2, chi2);
162 
163  const bool matchedVtx1(LArClusterHelper::GetClosestDistance(projVtx1, pCluster1) < m_maxDisplacement);
164  const bool matchedVtx2(LArClusterHelper::GetClosestDistance(projVtx2, pCluster2) < m_maxDisplacement);
165  const bool matchedVtx3(LArClusterHelper::GetClosestDistance(projVtx3, pCluster3) < m_maxDisplacement);
166 
167  const bool matchedEnd1(LArClusterHelper::GetClosestDistance(projEnd1, pCluster1) < m_maxDisplacement);
168  const bool matchedEnd2(LArClusterHelper::GetClosestDistance(projEnd2, pCluster2) < m_maxDisplacement);
169  const bool matchedEnd3(LArClusterHelper::GetClosestDistance(projEnd3, pCluster3) < m_maxDisplacement);
170 
171  const bool matchedCluster1(matchedVtx1 || matchedEnd1);
172  const bool matchedCluster2(matchedVtx2 || matchedEnd2);
173  const bool matchedCluster3(matchedVtx3 || matchedEnd3);
174  const bool matchedVtx(matchedVtx1 || matchedVtx2 || matchedVtx3);
175  const bool matchedEnd(matchedEnd1 || matchedEnd2 || matchedEnd3);
176 
177  if (matchedCluster1 && matchedCluster2 && matchedCluster3 && matchedVtx && matchedEnd)
178  return true;
179  }
180  }
181 
182  return false;
183 }
184 
185 //------------------------------------------------------------------------------------------------------------------------------------------
186 
187 void CosmicRayTrackMatchingAlgorithm::SetPfoParameters(const Particle &particle, PandoraContentApi::ParticleFlowObject::Parameters &pfoParameters) const
188 {
189  ClusterList clusterList;
190  if (particle.m_pClusterU)
191  clusterList.push_back(particle.m_pClusterU);
192  if (particle.m_pClusterV)
193  clusterList.push_back(particle.m_pClusterV);
194  if (particle.m_pClusterW)
195  clusterList.push_back(particle.m_pClusterW);
196 
197  // TODO Correct these placeholder parameters
198  pfoParameters.m_particleId = MU_MINUS; // TRACK
199  pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
200  pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
201  pfoParameters.m_energy = 0.f;
202  pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
203  pfoParameters.m_clusterList = clusterList;
204 }
205 
206 //------------------------------------------------------------------------------------------------------------------------------------------
207 
208 StatusCode CosmicRayTrackMatchingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
209 {
210  PANDORA_RETURN_RESULT_IF_AND_IF(
211  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ClusterMinLength", m_clusterMinLength));
212 
213  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "VtxXOverlap", m_vtxXOverlap));
214 
215  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlap", m_minXOverlap));
216 
217  PANDORA_RETURN_RESULT_IF_AND_IF(
218  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlapFraction", m_minXOverlapFraction));
219 
220  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDisplacement", m_maxDisplacement));
221 
223 }
224 
225 } // namespace lar_content
void SetPfoParameters(const Particle &protoParticle, PandoraContentApi::ParticleFlowObject::Parameters &pfoParameters) const
Calculate Pfo properties from proto particle.
void SelectCleanClusters(const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select a set of clusters judged to be clean.
float m_vtxXOverlap
requirement on X overlap of start/end positions
float m_clusterMinLength
minimum length of clusters for this algorithm
const pandora::Cluster * m_pClusterV
Address of cluster in V view.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
float m_maxDisplacement
requirement on 3D consistency checks
Header file for the cluster helper class.
const pandora::Cluster * m_pClusterW
Address of cluster in W view.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
const pandora::Cluster * m_pClusterU
Address of cluster in U view.
static float MergeTwoPositions(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const float position1, const float position2)
Merge two views (U,V) to give a third view (Z).
float m_minXOverlap
requirement on minimum X overlap for associated clusters
bool MatchClusters(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2) const
Match a pair of clusters from two views.
bool CheckMatchedClusters3D(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, const pandora::Cluster *const pCluster3) const
Check that three clusters have a consistent 3D position.
static void GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z) ...
float m_minXOverlapFraction
requirement on minimum X overlap fraction for associated clusters
HitType
Definition: HitType.h:12
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
Char_t n[5]
std::vector< art::Ptr< recob::Cluster > > ClusterVector
Header file for the cosmic ray track matching algorithm class.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.