LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
BranchSplittingAlgorithm.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 BranchSplittingAlgorithm::BranchSplittingAlgorithm() :
21  m_maxTransverseDisplacement(1.5f),
22  m_maxLongitudinalDisplacement(10.f),
23  m_minLongitudinalExtension(3.f),
24  m_minCosRelativeAngle(0.966f),
25  m_projectionAngularAllowance(20.f)
26 {
27 }
28 
29 //------------------------------------------------------------------------------------------------------------------------------------------
30 
31 void BranchSplittingAlgorithm::FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFit, const TwoDSlidingFitResult &principalSlidingFit,
32  CartesianVector &principalStartPosition, CartesianVector &branchSplitPosition, CartesianVector &branchSplitDirection) const
33 {
34  // Conventions:
35  // (1) Delta ray is split from the branch cluster
36  // (2) Delta ray occurs where the vertex of the principal cluster meets the vertex of the branch cluster
37  // Method loops over the inner and outer positions of the principal and branch clusters, trying all
38  // possible assignments of vertex and end position until a split is found
39  for (unsigned int principalForward = 0; principalForward < 2; ++principalForward)
40  {
41  const CartesianVector principalVertexPosition(
42  1 == principalForward ? principalSlidingFit.GetGlobalMinLayerPosition() : principalSlidingFit.GetGlobalMaxLayerPosition());
43  const CartesianVector principalEndPosition(
44  1 != principalForward ? principalSlidingFit.GetGlobalMinLayerPosition() : principalSlidingFit.GetGlobalMaxLayerPosition());
45  const CartesianVector principalVertexDirection(1 == principalForward ? principalSlidingFit.GetGlobalMinLayerDirection()
46  : principalSlidingFit.GetGlobalMaxLayerDirection() * -1.f);
47 
48  CartesianVector projectedBranchPosition(0.f, 0.f, 0.f);
49  bool projectedPositionFound(false), projectedPositionFail(false);
50 
51  for (unsigned int branchForward = 0; branchForward < 2; ++branchForward)
52  {
53  const CartesianVector branchVertexPosition(
54  1 == branchForward ? branchSlidingFit.GetGlobalMinLayerPosition() : branchSlidingFit.GetGlobalMaxLayerPosition());
55  const CartesianVector branchEndPosition(
56  1 != branchForward ? branchSlidingFit.GetGlobalMinLayerPosition() : branchSlidingFit.GetGlobalMaxLayerPosition());
57  const CartesianVector branchEndDirection(
58  1 != branchForward ? branchSlidingFit.GetGlobalMinLayerDirection() : branchSlidingFit.GetGlobalMaxLayerDirection() * -1.f);
59 
60  if (principalVertexDirection.GetDotProduct(branchEndDirection) < 0.5f)
61  continue;
62 
63  if ((principalEndPosition - branchEndPosition).GetMagnitudeSquared() < (principalVertexPosition - branchVertexPosition).GetMagnitudeSquared())
64  continue;
65 
66  // Project the principal vertex onto the branch cluster
67  try
68  {
69  if (!projectedPositionFound && !projectedPositionFail)
70  {
71  projectedBranchPosition = LArPointingClusterHelper::GetProjectedPosition(
72  principalVertexPosition, principalVertexDirection, branchSlidingFit.GetCluster(), m_projectionAngularAllowance);
73  projectedPositionFound = true;
74  }
75  }
76  catch (StatusCodeException &)
77  {
78  projectedPositionFail = true;
79  }
80 
81  if (!projectedPositionFound || projectedPositionFail)
82  continue;
83 
84  const float projectedDistanceSquared((projectedBranchPosition - principalVertexPosition).GetMagnitudeSquared());
85 
86  if (projectedDistanceSquared > m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)
87  continue;
88 
89  const float commonDistanceSquared((projectedBranchPosition - branchEndPosition).GetMagnitudeSquared());
90 
91  if (projectedDistanceSquared > commonDistanceSquared)
92  continue;
93 
94  const float replacementDistanceSquared((projectedBranchPosition - principalEndPosition).GetMagnitudeSquared());
95 
96  if (replacementDistanceSquared < m_minLongitudinalExtension * m_minLongitudinalExtension)
97  continue;
98 
99  const float branchDistanceSquared((projectedBranchPosition - branchVertexPosition).GetMagnitudeSquared());
100 
101  if (branchDistanceSquared > 4.f * replacementDistanceSquared)
102  continue;
103 
104  // Require that principal vertex and branch projection have good (and improved) pointing
105  bool foundSplit(false);
106 
107  const float halfWindowLength(branchSlidingFit.GetLayerFitHalfWindowLength());
108  const float deltaL(1 == branchForward ? +halfWindowLength : -halfWindowLength);
109 
110  float localL(0.f), localT(0.f);
111  CartesianVector forwardDirection(0.f, 0.f, 0.f);
112  branchSlidingFit.GetLocalPosition(projectedBranchPosition, localL, localT);
113 
114  if (STATUS_CODE_SUCCESS != branchSlidingFit.GetGlobalFitDirection(localL + deltaL, forwardDirection))
115  continue;
116 
117  CartesianVector projectedBranchDirection(1 == branchForward ? forwardDirection : forwardDirection * -1.f);
118  const float cosTheta(-projectedBranchDirection.GetDotProduct(principalVertexDirection));
119 
120  try
121  {
122  const float currentCosTheta(branchSlidingFit.GetCosScatteringAngle(localL));
123 
124  if (cosTheta < currentCosTheta)
125  continue;
126  }
127  catch (StatusCodeException &)
128  {
129  }
130 
131  float rT1(0.f), rL1(0.f), rT2(0.f), rL2(0.f);
132  LArPointingClusterHelper::GetImpactParameters(projectedBranchPosition, projectedBranchDirection, principalVertexPosition, rL1, rT1);
133  LArPointingClusterHelper::GetImpactParameters(principalVertexPosition, principalVertexDirection, projectedBranchPosition, rL2, rT2);
134 
136  {
137  foundSplit = true;
138  principalStartPosition = principalVertexPosition;
139  branchSplitPosition = projectedBranchPosition;
140  branchSplitDirection = projectedBranchDirection * -1.f;
141  }
142 
143  if (foundSplit)
144  return;
145  }
146  }
147 
148  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
149 }
150 
151 //------------------------------------------------------------------------------------------------------------------------------------------
152 
153 StatusCode BranchSplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
154 {
155  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
156  XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
157 
158  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
159  XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
160 
161  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
162  XmlHelper::ReadValue(xmlHandle, "MinLongitudinalExtension", m_minLongitudinalExtension));
163 
164  PANDORA_RETURN_RESULT_IF_AND_IF(
165  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosRelativeAngle", m_minCosRelativeAngle));
166 
167  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
168  XmlHelper::ReadValue(xmlHandle, "ProjectionAngularAllowance", m_projectionAngularAllowance));
169 
171 }
172 
173 } // namespace lar_content
static pandora::CartesianVector GetProjectedPosition(const pandora::CartesianVector &initialPosition, const pandora::CartesianVector &initialDirection, const pandora::Cluster *const pCluster, const float projectionAngularAllowance)
Get projected position on a cluster from a specified position and direction.
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
float GetCosScatteringAngle(const float rL) const
Get scattering angle for a given longitudinal coordinate.
TFile f
Definition: plotHisto.C:6
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
Header file for the branch splitting algorithm class.
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.
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.
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.