LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArPointingClusterHelper.cc
Go to the documentation of this file.
1 
11 
12 using namespace pandora;
13 
14 namespace lar_content
15 {
16 
17 float LArPointingClusterHelper::GetLengthSquared(const LArPointingCluster &pointingCluster)
18 {
19  const LArPointingCluster::Vertex &innerVertex(pointingCluster.GetInnerVertex());
20  const LArPointingCluster::Vertex &outerVertex(pointingCluster.GetOuterVertex());
21  return (innerVertex.GetPosition() - outerVertex.GetPosition()).GetMagnitudeSquared();
22 }
23 
24 //------------------------------------------------------------------------------------------------------------------------------------------
25 
26 float LArPointingClusterHelper::GetLength(const LArPointingCluster &pointingCluster)
27 {
28  return std::sqrt(LArPointingClusterHelper::GetLengthSquared(pointingCluster));
29 }
30 
31 //------------------------------------------------------------------------------------------------------------------------------------------
32 
33 bool LArPointingClusterHelper::IsNode(const CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex,
34  const float minLongitudinalDistance, const float maxTransverseDistance)
35 {
36  float rL(0.f), rT(0.f);
37  LArPointingClusterHelper::GetImpactParameters(daughterVertex.GetPosition(), daughterVertex.GetDirection(), parentVertex, rL, rT);
38 
39  if (std::fabs(rL) > std::fabs(minLongitudinalDistance) || rT > maxTransverseDistance)
40  return false;
41 
42  return true;
43 }
44 
45 //------------------------------------------------------------------------------------------------------------------------------------------
46 
47 bool LArPointingClusterHelper::IsEmission(const CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex,
48  const float minLongitudinalDistance, const float maxLongitudinalDistance, const float maxTransverseDistance, const float angularAllowance)
49 {
50  float rL(0.f), rT(0.f);
51  LArPointingClusterHelper::GetImpactParameters(daughterVertex.GetPosition(), daughterVertex.GetDirection(), parentVertex, rL, rT);
52 
53  if (std::fabs(rL) > std::fabs(minLongitudinalDistance) && (rL < 0 || rL > maxLongitudinalDistance))
54  return false;
55 
56  const float tanSqTheta(std::pow(std::tan(M_PI * angularAllowance / 180.f), 2.0));
57 
58  if (rT * rT > maxTransverseDistance * maxTransverseDistance + rL * rL * tanSqTheta)
59  return false;
60 
61  return true;
62 }
63 
64 //------------------------------------------------------------------------------------------------------------------------------------------
65 
66 CartesianVector LArPointingClusterHelper::GetProjectedPosition(const CartesianVector &vertexPosition,
67  const CartesianVector &vertexDirection, const pandora::Cluster *const pCluster, const float projectionAngularAllowance)
68 {
69  const CaloHit *pClosestCaloHit(nullptr);
70  float closestDistanceSquared(std::numeric_limits<float>::max());
71  const float minCosTheta(std::cos(M_PI * projectionAngularAllowance / 180.f));
72 
73  for (const OrderedCaloHitList::value_type &layerEntry : pCluster->GetOrderedCaloHitList())
74  {
75  for (const CaloHit *const pCaloHit : *layerEntry.second)
76  {
77  const CartesianVector hitProjection(pCaloHit->GetPositionVector() - vertexPosition);
78  const float distanceSquared(hitProjection.GetMagnitudeSquared());
79 
80  if (distanceSquared > std::numeric_limits<float>::epsilon())
81  {
82  // TODO Try to give more weight to on-axis projections
83  if (distanceSquared < closestDistanceSquared)
84  {
85  if (-hitProjection.GetUnitVector().GetDotProduct(vertexDirection) > minCosTheta)
86  {
87  pClosestCaloHit = pCaloHit;
88  closestDistanceSquared = distanceSquared;
89  }
90  }
91  }
92  else
93  {
94  return pCaloHit->GetPositionVector();
95  }
96  }
97  }
98 
99  if (pClosestCaloHit)
100  return pClosestCaloHit->GetPositionVector();
101 
102  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
103 }
104 
105 //------------------------------------------------------------------------------------------------------------------------------------------
106 
107 void LArPointingClusterHelper::GetClosestVertices(const bool useX, const bool useY, const bool useZ, const LArPointingCluster &pointingClusterI,
108  const LArPointingCluster &pointingClusterJ, LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
109 {
110  if (pointingClusterI.GetCluster() == pointingClusterJ.GetCluster())
111  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
112 
113  if (!useX && !useY && !useZ)
114  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
115 
116  for (unsigned int useInnerI = 0; useInnerI < 2; ++useInnerI)
117  {
118  const LArPointingCluster::Vertex &vtxI(useInnerI == 1 ? pointingClusterI.GetInnerVertex() : pointingClusterI.GetOuterVertex());
119  const LArPointingCluster::Vertex &endI(useInnerI == 0 ? pointingClusterI.GetInnerVertex() : pointingClusterI.GetOuterVertex());
120 
121  for (unsigned int useInnerJ = 0; useInnerJ < 2; ++useInnerJ)
122  {
123  const LArPointingCluster::Vertex &vtxJ(useInnerJ == 1 ? pointingClusterJ.GetInnerVertex() : pointingClusterJ.GetOuterVertex());
124  const LArPointingCluster::Vertex &endJ(useInnerJ == 0 ? pointingClusterJ.GetInnerVertex() : pointingClusterJ.GetOuterVertex());
125 
126  const float vtxI_vtxJ_dx(useX ? (vtxI.GetPosition().GetX() - vtxJ.GetPosition().GetX()) : 0.f);
127  const float vtxI_vtxJ_dy(useY ? (vtxI.GetPosition().GetY() - vtxJ.GetPosition().GetY()) : 0.f);
128  const float vtxI_vtxJ_dz(useZ ? (vtxI.GetPosition().GetZ() - vtxJ.GetPosition().GetZ()) : 0.f);
129  const float vtxI_vtxJ(vtxI_vtxJ_dx * vtxI_vtxJ_dx + vtxI_vtxJ_dy * vtxI_vtxJ_dy + vtxI_vtxJ_dz * vtxI_vtxJ_dz);
130 
131  const float vtxI_endJ_dx(useX ? (vtxI.GetPosition().GetX() - endJ.GetPosition().GetX()) : 0.f);
132  const float vtxI_endJ_dy(useY ? (vtxI.GetPosition().GetY() - endJ.GetPosition().GetY()) : 0.f);
133  const float vtxI_endJ_dz(useZ ? (vtxI.GetPosition().GetZ() - endJ.GetPosition().GetZ()) : 0.f);
134  const float vtxI_endJ(vtxI_endJ_dx * vtxI_endJ_dx + vtxI_endJ_dy * vtxI_endJ_dy + vtxI_endJ_dz * vtxI_endJ_dz);
135 
136  const float endI_vtxJ_dx(useX ? (endI.GetPosition().GetX() - vtxJ.GetPosition().GetX()) : 0.f);
137  const float endI_vtxJ_dy(useY ? (endI.GetPosition().GetY() - vtxJ.GetPosition().GetY()) : 0.f);
138  const float endI_vtxJ_dz(useZ ? (endI.GetPosition().GetZ() - vtxJ.GetPosition().GetZ()) : 0.f);
139  const float endI_vtxJ(endI_vtxJ_dx * endI_vtxJ_dx + endI_vtxJ_dy * endI_vtxJ_dy + endI_vtxJ_dz * endI_vtxJ_dz);
140 
141  const float endI_endJ_dx(useX ? (endI.GetPosition().GetX() - endJ.GetPosition().GetX()) : 0.f);
142  const float endI_endJ_dy(useY ? (endI.GetPosition().GetY() - endJ.GetPosition().GetY()) : 0.f);
143  const float endI_endJ_dz(useZ ? (endI.GetPosition().GetZ() - endJ.GetPosition().GetZ()) : 0.f);
144  const float endI_endJ(endI_endJ_dx * endI_endJ_dx + endI_endJ_dy * endI_endJ_dy + endI_endJ_dz * endI_endJ_dz);
145 
146  if ((vtxI_vtxJ < std::min(vtxI_endJ, std::min(endI_vtxJ, endI_endJ))) && (endI_endJ > std::max(vtxI_endJ, std::max(endI_vtxJ, vtxI_vtxJ))))
147  {
148  closestVertexI = vtxI;
149  closestVertexJ = vtxJ;
150  return;
151  }
152  }
153  }
154 
155  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
156 }
157 
158 //------------------------------------------------------------------------------------------------------------------------------------------
159 
160 void LArPointingClusterHelper::GetClosestVertices(const LArPointingCluster &pointingClusterI, const LArPointingCluster &pointingClusterJ,
161  LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
162 {
163  return LArPointingClusterHelper::GetClosestVertices(true, true, true, pointingClusterI, pointingClusterJ, closestVertexI, closestVertexJ);
164 }
165 
166 //------------------------------------------------------------------------------------------------------------------------------------------
167 
168 void LArPointingClusterHelper::GetClosestVerticesInX(const LArPointingCluster &pointingClusterI, const LArPointingCluster &pointingClusterJ,
169  LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
170 {
171  return LArPointingClusterHelper::GetClosestVertices(true, false, false, pointingClusterI, pointingClusterJ, closestVertexI, closestVertexJ);
172 }
173 
174 //------------------------------------------------------------------------------------------------------------------------------------------
175 
176 void LArPointingClusterHelper::GetClosestVerticesInYZ(const LArPointingCluster &pointingClusterI,
177  const LArPointingCluster &pointingClusterJ, LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
178 {
179  return LArPointingClusterHelper::GetClosestVertices(false, true, true, pointingClusterI, pointingClusterJ, closestVertexI, closestVertexJ);
180 }
181 
182 //------------------------------------------------------------------------------------------------------------------------------------------
183 
184 void LArPointingClusterHelper::GetImpactParametersInYZ(
185  const LArPointingCluster::Vertex &initialVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
186 {
187  if (std::fabs(initialVertex.GetDirection().GetX()) > 1.f - std::numeric_limits<float>::epsilon())
188  throw pandora::StatusCodeException(pandora::STATUS_CODE_NOT_FOUND);
189 
190  const pandora::CartesianVector initialPosition(0.f, initialVertex.GetPosition().GetY(), initialVertex.GetPosition().GetZ());
191  const pandora::CartesianVector initialDirection(0.f, initialVertex.GetDirection().GetY(), initialVertex.GetDirection().GetZ());
192  const pandora::CartesianVector targetPosition(0.f, targetVertex.GetPosition().GetY(), targetVertex.GetPosition().GetZ());
193 
194  return LArPointingClusterHelper::GetImpactParameters(initialPosition, initialDirection.GetUnitVector(), targetPosition, longitudinal, transverse);
195 }
196 
197 //------------------------------------------------------------------------------------------------------------------------------------------
198 
199 void LArPointingClusterHelper::GetImpactParameters(
200  const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
201 {
202  return LArPointingClusterHelper::GetImpactParameters(
203  pointingVertex.GetPosition(), pointingVertex.GetDirection(), targetVertex.GetPosition(), longitudinal, transverse);
204 }
205 
206 //------------------------------------------------------------------------------------------------------------------------------------------
207 
208 void LArPointingClusterHelper::GetImpactParameters(
209  const LArPointingCluster::Vertex &pointingVertex, const CartesianVector &targetPosition, float &longitudinal, float &transverse)
210 {
211  return LArPointingClusterHelper::GetImpactParameters(
212  pointingVertex.GetPosition(), pointingVertex.GetDirection(), targetPosition, longitudinal, transverse);
213 }
214 
215 //------------------------------------------------------------------------------------------------------------------------------------------
216 
217 void LArPointingClusterHelper::GetImpactParameters(const CartesianVector &initialPosition, const CartesianVector &initialDirection,
218  const CartesianVector &targetPosition, float &longitudinal, float &transverse)
219 {
220  // sign convention for longitudinal distance:
221  // -positive value means initial position is downstream of target position
222  transverse = initialDirection.GetCrossProduct(targetPosition - initialPosition).GetMagnitude();
223  longitudinal = -initialDirection.GetDotProduct(targetPosition - initialPosition);
224 }
225 
226 //------------------------------------------------------------------------------------------------------------------------------------------
227 
228 void LArPointingClusterHelper::GetAverageDirection(
229  const LArPointingCluster::Vertex &firstVertex, const LArPointingCluster::Vertex &secondVertex, CartesianVector &averageDirection)
230 {
231  const Cluster *const pFirstCluster(firstVertex.GetCluster());
232  const Cluster *const pSecondCluster(secondVertex.GetCluster());
233 
234  if (pFirstCluster == pSecondCluster)
235  throw pandora::StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
236 
237  const float energy1(pFirstCluster->GetHadronicEnergy());
238  const float energy2(pSecondCluster->GetHadronicEnergy());
239 
240  if ((energy1 < std::numeric_limits<float>::epsilon()) || (energy2 < std::numeric_limits<float>::epsilon()))
241  throw pandora::StatusCodeException(STATUS_CODE_NOT_ALLOWED);
242 
243  averageDirection = (firstVertex.GetDirection() * energy1 + secondVertex.GetDirection() * energy2).GetUnitVector();
244 }
245 
246 //------------------------------------------------------------------------------------------------------------------------------------------
247 
249  const LArPointingCluster::Vertex &secondVertex, CartesianVector &intersectPosition, float &firstDisplacement, float &secondDisplacement)
250 {
251  const Cluster *const pFirstCluster(firstVertex.GetCluster());
252  const Cluster *const pSecondCluster(secondVertex.GetCluster());
253 
254  if (pFirstCluster == pSecondCluster)
255  throw pandora::StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
256 
257  LArPointingClusterHelper::GetIntersection(firstVertex.GetPosition(), firstVertex.GetDirection(), secondVertex.GetPosition(),
258  secondVertex.GetDirection(), intersectPosition, firstDisplacement, secondDisplacement);
259 }
260 
261 //------------------------------------------------------------------------------------------------------------------------------------------
262 
263 void LArPointingClusterHelper::GetIntersection(const CartesianVector &a1, const CartesianVector &a2, const CartesianVector &b1,
264  const CartesianVector &b2, CartesianVector &intersectPosition, float &firstDisplacement, float &secondDisplacement)
265 {
266  // note: input lines are r = a1 + P * a2 and r = b1 + Q * b2
267  const float cosTheta = a2.GetDotProduct(b2);
268 
269  // lines must be non-parallel
270  if (1.f - std::fabs(cosTheta) < std::numeric_limits<float>::epsilon())
271  throw pandora::StatusCodeException(STATUS_CODE_NOT_ALLOWED);
272 
273  // calculate the intersection (by minimising the distance between the lines)
274  const float P = ((a2 - b2 * cosTheta).GetDotProduct(b1 - a1)) / (1.f - cosTheta * cosTheta);
275  const float Q = ((a2 * cosTheta - b2).GetDotProduct(b1 - a1)) / (1.f - cosTheta * cosTheta);
276 
277  // position of intersection (or point of closest approach)
278  intersectPosition = (a1 + a2 * P + b1 + b2 * Q) * 0.5f;
279 
280  // displacements of intersection from input vertices
281  firstDisplacement = P;
282  secondDisplacement = Q;
283 }
284 
285 //------------------------------------------------------------------------------------------------------------------------------------------
286 
287 void LArPointingClusterHelper::GetIntersection(const LArPointingCluster::Vertex &vertexCluster, const Cluster *const pTargetCluster,
288  CartesianVector &intersectPosition, float &displacementL, float &displacementT)
289 {
290  displacementT = +std::numeric_limits<float>::max();
291  displacementL = -std::numeric_limits<float>::max();
292 
293  float rL(0.f), rT(0.f);
294  float figureOfMerit(std::numeric_limits<float>::max());
295  float foundIntersection(false);
296 
297  const OrderedCaloHitList &orderedCaloHitList(pTargetCluster->GetOrderedCaloHitList());
298 
299  for (OrderedCaloHitList::const_iterator iter1 = orderedCaloHitList.begin(), iterEnd1 = orderedCaloHitList.end(); iter1 != iterEnd1; ++iter1)
300  {
301  for (CaloHitList::const_iterator iter2 = iter1->second->begin(), iterEnd2 = iter1->second->end(); iter2 != iterEnd2; ++iter2)
302  {
303  const CartesianVector &hitPosition = (*iter2)->GetPositionVector();
304 
305  LArPointingClusterHelper::GetImpactParameters(vertexCluster.GetPosition(), vertexCluster.GetDirection(), hitPosition, rL, rT);
306 
307  if (rT < figureOfMerit)
308  {
309  figureOfMerit = rT;
310 
311  displacementL = rL;
312  displacementT = rT;
313  intersectPosition = hitPosition;
314  foundIntersection = true;
315  }
316  }
317  }
318 
319  if (!foundIntersection)
320  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
321 }
322 
323 //------------------------------------------------------------------------------------------------------------------------------------------
324 
325 LArPointingCluster::Vertex LArPointingClusterHelper::GetBestVertexEstimate(const LArPointingClusterVertexList &vertexList,
326  const LArPointingClusterList &pointingClusterList, const float minLongitudinalDistance, const float maxLongitudinalDistance,
327  const float maxTransverseDistance, const float angularAllowance)
328 {
329  float bestAssociatedEnergy(0.f);
330  LArPointingClusterVertexList::const_iterator bestVertexIter(vertexList.end());
331 
332  for (LArPointingClusterVertexList::const_iterator iter = vertexList.begin(), iterEnd = vertexList.end(); iter != iterEnd; ++iter)
333  {
334  const LArPointingCluster::Vertex &vertex(*iter);
335 
336  LArPointingClusterVertexList associatedVertices;
337  LArPointingClusterHelper::CollectAssociatedClusters(vertex, pointingClusterList, minLongitudinalDistance, maxLongitudinalDistance,
338  maxTransverseDistance, angularAllowance, associatedVertices);
339 
340  const float associatedEnergy(LArPointingClusterHelper::GetAssociatedEnergy(vertex, associatedVertices));
341 
342  if (associatedEnergy > bestAssociatedEnergy)
343  {
344  bestVertexIter = iter;
345  bestAssociatedEnergy = associatedEnergy;
346  }
347  }
348 
349  if (vertexList.end() == bestVertexIter)
350  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
351 
352  return (*bestVertexIter);
353 }
354 
355 //------------------------------------------------------------------------------------------------------------------------------------------
356 
357 void LArPointingClusterHelper::CollectAssociatedClusters(const LArPointingCluster::Vertex &vertex, const LArPointingClusterList &inputList,
358  const float minLongitudinalDistance, const float maxLongitudinalDistance, const float maxTransverseDistance,
359  const float angularAllowance, LArPointingClusterVertexList &outputList)
360 {
361  for (LArPointingClusterList::const_iterator iter = inputList.begin(), iterEnd = inputList.end(); iter != iterEnd; ++iter)
362  {
363  const LArPointingCluster &pointingCluster = *iter;
364  const LArPointingCluster::Vertex &innerVertex = pointingCluster.GetInnerVertex();
365  const LArPointingCluster::Vertex &outerVertex = pointingCluster.GetOuterVertex();
366 
367  const float innerDistanceSquared = (innerVertex.GetPosition() - vertex.GetPosition()).GetMagnitudeSquared();
368  const float outerDistanceSquared = (outerVertex.GetPosition() - vertex.GetPosition()).GetMagnitudeSquared();
369 
370  const LArPointingCluster::Vertex &chosenVertex((innerDistanceSquared < outerDistanceSquared) ? innerVertex : outerVertex);
371 
372  if (LArPointingClusterHelper::IsNode(vertex.GetPosition(), chosenVertex, minLongitudinalDistance, maxTransverseDistance) ||
373  LArPointingClusterHelper::IsEmission(vertex.GetPosition(), chosenVertex, minLongitudinalDistance, maxLongitudinalDistance,
374  maxTransverseDistance, angularAllowance))
375  {
376  outputList.push_back(chosenVertex);
377  }
378  }
379 }
380 
381 //------------------------------------------------------------------------------------------------------------------------------------------
382 
383 float LArPointingClusterHelper::GetAssociatedEnergy(const LArPointingCluster::Vertex &vertex, const LArPointingClusterVertexList &associatedVertices)
384 {
385  float associatedEnergy(0.f);
386 
387  for (LArPointingClusterVertexList::const_iterator iter = associatedVertices.begin(), iterEnd = associatedVertices.end(); iter != iterEnd; ++iter)
388  {
389  const LArPointingCluster::Vertex &clusterVertex(*iter);
390  const Cluster *const pCluster(clusterVertex.GetCluster());
391 
392  const float clusterEnergy(LArClusterHelper::GetEnergyFromLength(pCluster));
393  const float clusterLength(LArClusterHelper::GetLength(pCluster));
394  const float deltaLength(clusterVertex.GetDirection().GetDotProduct(vertex.GetPosition() - clusterVertex.GetPosition()));
395 
396  if (deltaLength < std::numeric_limits<float>::epsilon())
397  {
398  associatedEnergy += clusterEnergy;
399  }
400  else if (deltaLength < clusterLength)
401  {
402  associatedEnergy += clusterEnergy * (1.f - (deltaLength / clusterLength));
403  }
404  }
405 
406  return associatedEnergy;
407 }
408 
409 } // namespace lar_content
std::vector< LArPointingCluster > LArPointingClusterList
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
std::pair< float, float > GetIntersection(double Ax, double Ay, double Bx, double By, double Cx, double Cy, double Dx, double Dy)
Definition: Polygon2D.cxx:36
std::vector< LArPointingCluster::Vertex > LArPointingClusterVertexList
TFile f
Definition: plotHisto.C:6
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
#define a2
Header file for the cluster helper class.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
const pandora::CartesianVector & GetDirection() const
Get the vertex direction.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
#define a1
vertex reconstruction