LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
DeltaRaySplittingAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
14 
15 using namespace pandora;
16 
17 namespace lar_content
18 {
19 
20 DeltaRaySplittingAlgorithm::DeltaRaySplittingAlgorithm() :
21  m_stepSize(1.f),
22  m_maxTransverseDisplacement(1.5f),
23  m_maxLongitudinalDisplacement(10.f),
24  m_minCosRelativeAngle(0.985f)
25 {
26 }
27 
28 //------------------------------------------------------------------------------------------------------------------------------------------
29 
30 void DeltaRaySplittingAlgorithm::FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFit, const TwoDSlidingFitResult &principalSlidingFit,
31  CartesianVector &principalStartPosition, CartesianVector &branchSplitPosition, CartesianVector &branchSplitDirection) const
32 {
33  // Conventions:
34  // (1) Delta ray is split from the branch cluster
35  // (2) Delta ray occurs where the vertex of the principal cluster meets the vertex of the branch cluster
36  // Method loops over the inner and outer positions of the principal and branch clusters, trying all
37  // possible assignments of vertex and end position until a split is found
38 
39  for (unsigned int principalForward = 0; principalForward < 2; ++principalForward)
40  {
41  const CartesianVector principalVertex(
42  1 == principalForward ? principalSlidingFit.GetGlobalMinLayerPosition() : principalSlidingFit.GetGlobalMaxLayerPosition());
43  const CartesianVector principalEnd(
44  1 == principalForward ? principalSlidingFit.GetGlobalMaxLayerPosition() : principalSlidingFit.GetGlobalMinLayerPosition());
45  const CartesianVector principalDirection(1 == principalForward ? principalSlidingFit.GetGlobalMinLayerDirection()
46  : principalSlidingFit.GetGlobalMaxLayerDirection() * -1.f);
47 
48  if (LArClusterHelper::GetClosestDistance(principalVertex, branchSlidingFit.GetCluster()) > m_maxLongitudinalDisplacement)
49  continue;
50 
51  for (unsigned int branchForward = 0; branchForward < 2; ++branchForward)
52  {
53  const CartesianVector branchVertex(
54  1 == branchForward ? branchSlidingFit.GetGlobalMinLayerPosition() : branchSlidingFit.GetGlobalMaxLayerPosition());
55  const CartesianVector branchEnd(
56  1 == branchForward ? branchSlidingFit.GetGlobalMaxLayerPosition() : branchSlidingFit.GetGlobalMinLayerPosition());
57  const CartesianVector branchDirection(
58  1 == branchForward ? branchSlidingFit.GetGlobalMinLayerDirection() : branchSlidingFit.GetGlobalMaxLayerDirection() * -1.f);
59 
60  // Require vertices to be closest two ends
61  const float vertex_to_vertex((principalVertex - branchVertex).GetMagnitudeSquared());
62  const float vertex_to_end((principalVertex - branchEnd).GetMagnitudeSquared());
63  const float end_to_vertex((principalEnd - branchVertex).GetMagnitudeSquared());
64  const float end_to_end((principalEnd - branchEnd).GetMagnitudeSquared());
65 
66  // (sign convention for vertexProjection: positive means that clusters overlap)
67  const float vertexProjection(+branchDirection.GetDotProduct(principalVertex - branchVertex));
68  const float cosRelativeAngle(-branchDirection.GetDotProduct(principalDirection));
69 
70  if (vertex_to_vertex > std::min(end_to_end, std::min(vertex_to_end, end_to_vertex)))
71  continue;
72 
73  if (end_to_end < std::max(vertex_to_vertex, std::max(vertex_to_end, end_to_vertex)))
74  continue;
75 
76  if (vertexProjection < 0.f && cosRelativeAngle > m_minCosRelativeAngle)
77  continue;
78 
79  if (cosRelativeAngle < 0.f)
80  continue;
81 
82  // Serach for a split by winding back the branch cluster sliding fit
83  bool foundSplit(false);
84 
85  const float halfWindowLength(branchSlidingFit.GetLayerFitHalfWindowLength());
86  const float deltaL(1 == branchForward ? +halfWindowLength : -halfWindowLength);
87 
88  float branchDistance(std::max(0.f, vertexProjection) + 0.5f * m_stepSize);
89 
90  while (!foundSplit)
91  {
92  branchDistance += m_stepSize;
93 
94  const CartesianVector linearProjection(branchVertex + branchDirection * branchDistance);
95 
96  if (principalDirection.GetDotProduct(linearProjection - principalVertex) < -m_maxLongitudinalDisplacement)
97  break;
98 
99  if ((linearProjection - branchVertex).GetMagnitudeSquared() > (linearProjection - branchEnd).GetMagnitudeSquared())
100  break;
101 
102  float localL(0.f), localT(0.f);
103  CartesianVector truncatedPosition(0.f, 0.f, 0.f);
104  CartesianVector forwardDirection(0.f, 0.f, 0.f);
105  branchSlidingFit.GetLocalPosition(linearProjection, localL, localT);
106 
107  if ((STATUS_CODE_SUCCESS != branchSlidingFit.GetGlobalFitPosition(localL, truncatedPosition)) ||
108  (STATUS_CODE_SUCCESS != branchSlidingFit.GetGlobalFitDirection(localL + deltaL, forwardDirection)))
109  {
110  continue;
111  }
112 
113  CartesianVector truncatedDirection(1 == branchForward ? forwardDirection : forwardDirection * -1.f);
114  const float cosTheta(-truncatedDirection.GetDotProduct(principalDirection));
115 
116  float rT1(0.f), rL1(0.f), rT2(0.f), rL2(0.f);
117  LArPointingClusterHelper::GetImpactParameters(truncatedPosition, truncatedDirection, principalVertex, rL1, rT1);
118  LArPointingClusterHelper::GetImpactParameters(principalVertex, principalDirection, truncatedPosition, rL2, rT2);
119 
120  if ((cosTheta > m_minCosRelativeAngle) && (rT1 < m_maxTransverseDisplacement) && (rT2 < m_maxTransverseDisplacement))
121  {
122  foundSplit = true;
123  principalStartPosition = principalVertex;
124  branchSplitPosition = truncatedPosition;
125  branchSplitDirection = truncatedDirection * -1.f;
126  }
127  }
128 
129  if (foundSplit)
130  return;
131  }
132  }
133 
134  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
135 }
136 
137 //------------------------------------------------------------------------------------------------------------------------------------------
138 
139 StatusCode DeltaRaySplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
140 {
141  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "StepSize", m_stepSize));
142 
143  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
144  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
145 
146  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
147  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
148 
149  PANDORA_RETURN_RESULT_IF_AND_IF(
150  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosRelativeAngle", m_minCosRelativeAngle));
151 
153 }
154 
155 } // namespace lar_content
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
void FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFit, const TwoDSlidingFitResult &replacementSlidingFit, pandora::CartesianVector &replacementStartPosition, pandora::CartesianVector &branchSplitPosition, pandora::CartesianVector &branchSplitDirection) const
Output the best split positions in branch and replacement clusters.
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
TFile f
Definition: plotHisto.C:6
Header file for the delta ray splitting algorithm class.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.