LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
OvershootTracksTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
14 
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 OvershootTracksTool::OvershootTracksTool() :
24  m_splitMode(true),
25  m_maxVertexXSeparation(2.f),
26  m_cosThetaCutForKinkSearch(0.94f)
27 {
28 }
29 
30 //------------------------------------------------------------------------------------------------------------------------------------------
31 
33  ThreeViewTransverseTracksAlgorithm *const pAlgorithm, const IteratorList &iteratorList, ModificationList &modificationList) const
34 {
35  for (IteratorList::const_iterator iIter1 = iteratorList.begin(), iIter1End = iteratorList.end(); iIter1 != iIter1End; ++iIter1)
36  {
37  for (IteratorList::const_iterator iIter2 = iIter1; iIter2 != iIter1End; ++iIter2)
38  {
39  if (iIter1 == iIter2)
40  continue;
41 
42  try
43  {
44  const unsigned int nMatchedSamplingPoints1((*iIter1)->GetOverlapResult().GetNMatchedSamplingPoints());
45  const unsigned int nMatchedSamplingPoints2((*iIter2)->GetOverlapResult().GetNMatchedSamplingPoints());
46  IteratorList::const_iterator iIterA((nMatchedSamplingPoints1 >= nMatchedSamplingPoints2) ? iIter1 : iIter2);
47  IteratorList::const_iterator iIterB((nMatchedSamplingPoints1 >= nMatchedSamplingPoints2) ? iIter2 : iIter1);
48 
49  Particle particle(*(*iIterA), *(*iIterB));
50  const LArPointingCluster pointingClusterA1(pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA1));
51  const LArPointingCluster pointingClusterB1(pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB1));
52  const LArPointingCluster pointingClusterA2(pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA2));
53  const LArPointingCluster pointingClusterB2(pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB2));
54 
55  LArPointingCluster::Vertex vertexA1, vertexB1, vertexA2, vertexB2;
56  LArPointingClusterHelper::GetClosestVerticesInX(pointingClusterA1, pointingClusterB1, vertexA1, vertexB1);
57  LArPointingClusterHelper::GetClosestVerticesInX(pointingClusterA2, pointingClusterB2, vertexA2, vertexB2);
58 
59  if (!this->PassesVertexCuts(vertexA1, vertexB1) || !this->PassesVertexCuts(vertexA2, vertexB2))
60  continue;
61 
62  this->SetSplitPosition(vertexA1, vertexA2, vertexB1, vertexB2, particle);
63 
64  const bool isA1LowestInX(this->IsALowestInX(pointingClusterA1, pointingClusterB1));
65  const bool isA2LowestInX(this->IsALowestInX(pointingClusterA2, pointingClusterB2));
66  const bool isThreeDKink(m_majorityRulesMode ? true : this->IsThreeDKink(pAlgorithm, particle, isA1LowestInX, isA2LowestInX));
67 
68  if (isThreeDKink != m_splitMode)
69  continue;
70 
71  // Construct the modification object
72  Modification modification;
73 
74  if (m_splitMode)
75  {
76  modification.m_splitPositionMap[particle.m_pCommonCluster].push_back(particle.m_splitPosition);
77  }
78  else
79  {
80  const bool vertex1AIsLowX(vertexA1.GetPosition().GetX() < vertexB1.GetPosition().GetX());
81  const Cluster *const pLowXCluster1(vertex1AIsLowX ? particle.m_pClusterA1 : particle.m_pClusterB1);
82  const Cluster *const pHighXCluster1(vertex1AIsLowX ? particle.m_pClusterB1 : particle.m_pClusterA1);
83  modification.m_clusterMergeMap[pLowXCluster1].push_back(pHighXCluster1);
84 
85  const bool vertex2AIsLowX(vertexA2.GetPosition().GetX() < vertexB2.GetPosition().GetX());
86  const Cluster *const pLowXCluster2(vertex2AIsLowX ? particle.m_pClusterA2 : particle.m_pClusterB2);
87  const Cluster *const pHighXCluster2(vertex2AIsLowX ? particle.m_pClusterB2 : particle.m_pClusterA2);
88  modification.m_clusterMergeMap[pLowXCluster2].push_back(pHighXCluster2);
89  }
90 
91  modification.m_affectedClusters.push_back(particle.m_pCommonCluster);
92  modification.m_affectedClusters.push_back(particle.m_pClusterA1);
93  modification.m_affectedClusters.push_back(particle.m_pClusterA2);
94  modification.m_affectedClusters.push_back(particle.m_pClusterB1);
95  modification.m_affectedClusters.push_back(particle.m_pClusterB2);
96 
97  modificationList.push_back(modification);
98  }
99  catch (StatusCodeException &statusCodeException)
100  {
101  if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
102  throw statusCodeException;
103 
104  continue;
105  }
106  }
107  }
108 }
109 
110 //------------------------------------------------------------------------------------------------------------------------------------------
111 
113 {
114  float longitudinalAB(-std::numeric_limits<float>::max()), transverseAB(std::numeric_limits<float>::max());
115  LArPointingClusterHelper::GetImpactParameters(vertexA, vertexB, longitudinalAB, transverseAB);
116 
117  float longitudinalBA(-std::numeric_limits<float>::max()), transverseBA(std::numeric_limits<float>::max());
118  LArPointingClusterHelper::GetImpactParameters(vertexB, vertexA, longitudinalBA, transverseBA);
119 
120  if (std::min(longitudinalAB, longitudinalBA) < m_minLongitudinalImpactParameter)
121  return false;
122 
123  return true;
124 }
125 
126 //------------------------------------------------------------------------------------------------------------------------------------------
127 
129  const LArPointingCluster::Vertex &vertexB1, const LArPointingCluster::Vertex &vertexB2, Particle &particle) const
130 {
131  bool splitAtElementA(false), splitAtElementB(false);
132 
133  if (std::fabs(vertexA1.GetPosition().GetX() - vertexA2.GetPosition().GetX()) < m_maxVertexXSeparation)
134  {
135  splitAtElementA = true;
136  }
137  else if (std::fabs(vertexB1.GetPosition().GetX() - vertexB2.GetPosition().GetX()) < m_maxVertexXSeparation)
138  {
139  splitAtElementB = true;
140  }
141 
142  if (!splitAtElementA && !splitAtElementB)
143  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
144 
145  particle.m_splitPosition1 = splitAtElementA ? vertexA1.GetPosition() : vertexB1.GetPosition();
146  particle.m_splitPosition2 = splitAtElementA ? vertexA2.GetPosition() : vertexB2.GetPosition();
147 
148  CartesianVector splitPosition(0.f, 0.f, 0.f);
149  float chiSquared(std::numeric_limits<float>::max());
151  LArClusterHelper::GetClusterHitType(particle.m_pClusterA2), particle.m_splitPosition1, particle.m_splitPosition2, splitPosition, chiSquared);
152 
153  particle.m_splitPosition = splitPosition;
154 }
155 
156 //------------------------------------------------------------------------------------------------------------------------------------------
157 
159  ThreeViewTransverseTracksAlgorithm *const pAlgorithm, const Particle &particle, const bool isA1LowestInX, const bool isA2LowestInX) const
160 {
161  try
162  {
163  const TwoDSlidingFitResult &lowXFitResult1(isA1LowestInX ? pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA1)
164  : pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB1));
165  const TwoDSlidingFitResult &highXFitResult1(isA1LowestInX ? pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB1)
166  : pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA1));
167  const TwoDSlidingFitResult &lowXFitResult2(isA2LowestInX ? pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA2)
168  : pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB2));
169  const TwoDSlidingFitResult &highXFitResult2(isA2LowestInX ? pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterB2)
170  : pAlgorithm->GetCachedSlidingFitResult(particle.m_pClusterA2));
171  const TwoDSlidingFitResult &fitResultCommon3(pAlgorithm->GetCachedSlidingFitResult(particle.m_pCommonCluster));
172 
173  const float minusX(this->GetXSamplingPoint(particle.m_splitPosition, false, fitResultCommon3, lowXFitResult1, lowXFitResult2));
174  const float plusX(this->GetXSamplingPoint(particle.m_splitPosition, true, fitResultCommon3, highXFitResult1, highXFitResult2));
175 
176  CartesianVector minus1(0.f, 0.f, 0.f), split1(particle.m_splitPosition1), plus1(0.f, 0.f, 0.f);
177  CartesianVector minus2(0.f, 0.f, 0.f), split2(particle.m_splitPosition2), plus2(0.f, 0.f, 0.f);
178  CartesianVector minus3(0.f, 0.f, 0.f), split3(particle.m_splitPosition), plus3(0.f, 0.f, 0.f);
179 
180  if ((STATUS_CODE_SUCCESS != lowXFitResult1.GetGlobalFitPositionAtX(minusX, minus1)) ||
181  (STATUS_CODE_SUCCESS != highXFitResult1.GetGlobalFitPositionAtX(plusX, plus1)) ||
182  (STATUS_CODE_SUCCESS != lowXFitResult2.GetGlobalFitPositionAtX(minusX, minus2)) ||
183  (STATUS_CODE_SUCCESS != highXFitResult2.GetGlobalFitPositionAtX(plusX, plus2)) ||
184  (STATUS_CODE_SUCCESS != fitResultCommon3.GetGlobalFitPositionAtX(minusX, minus3)) ||
185  (STATUS_CODE_SUCCESS != fitResultCommon3.GetGlobalFitPositionAtX(plusX, plus3)))
186  {
187  return true; // majority rules, by default
188  }
189 
190  // Extract results
191  const HitType hitType1(LArClusterHelper::GetClusterHitType(particle.m_pClusterA1));
192  const HitType hitType2(LArClusterHelper::GetClusterHitType(particle.m_pClusterA2));
194 
195  CartesianVector minus(0.f, 0.f, 0.f), split(0.f, 0.f, 0.f), plus(0.f, 0.f, 0.f);
196  float chi2Minus(std::numeric_limits<float>::max()), chi2Split(std::numeric_limits<float>::max()),
197  chi2Plus(std::numeric_limits<float>::max());
198  LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), hitType1, hitType2, hitType3, minus1, minus2, minus3, minus, chi2Minus);
199  LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), hitType1, hitType2, hitType3, split1, split2, split3, split, chi2Split);
200  LArGeometryHelper::MergeThreePositions3D(this->GetPandora(), hitType1, hitType2, hitType3, plus1, plus2, plus3, plus, chi2Plus);
201 
202  // Apply final cuts
203  const CartesianVector minusToSplit((split - minus).GetUnitVector());
204  const CartesianVector splitToPlus((plus - split).GetUnitVector());
205  const float dotProduct(minusToSplit.GetDotProduct(splitToPlus));
206 
207  if (dotProduct > m_cosThetaCutForKinkSearch)
208  return false;
209  }
210  catch (StatusCodeException &)
211  {
212  }
213 
214  return true;
215 }
216 
217 //------------------------------------------------------------------------------------------------------------------------------------------
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
220 OvershootTracksTool::Particle::Particle(const TensorType::Element &elementA, const TensorType::Element &elementB) :
221  m_splitPosition(0.f, 0.f, 0.f),
222  m_splitPosition1(0.f, 0.f, 0.f),
223  m_splitPosition2(0.f, 0.f, 0.f)
224 {
225  const HitType commonView((elementA.GetClusterU() == elementB.GetClusterU()) ? TPC_VIEW_U
226  : (elementA.GetClusterV() == elementB.GetClusterV()) ? TPC_VIEW_V
227  : (elementA.GetClusterW() == elementB.GetClusterW()) ? TPC_VIEW_W
228  : HIT_CUSTOM);
229 
230  if (HIT_CUSTOM == commonView)
231  throw StatusCodeException(STATUS_CODE_FAILURE);
232 
233  m_pCommonCluster = (TPC_VIEW_U == commonView) ? elementA.GetClusterU()
234  : (TPC_VIEW_V == commonView) ? elementA.GetClusterV()
235  : elementA.GetClusterW();
236  m_pClusterA1 = (TPC_VIEW_U == commonView) ? elementA.GetClusterV() : elementA.GetClusterU();
237  m_pClusterA2 = (TPC_VIEW_U == commonView) ? elementA.GetClusterW()
238  : (TPC_VIEW_V == commonView) ? elementA.GetClusterW()
239  : elementA.GetClusterV();
240  m_pClusterB1 = (TPC_VIEW_U == commonView) ? elementB.GetClusterV() : elementB.GetClusterU();
241  m_pClusterB2 = (TPC_VIEW_U == commonView) ? elementB.GetClusterW()
242  : (TPC_VIEW_V == commonView) ? elementB.GetClusterW()
243  : elementB.GetClusterV();
244 
246  throw StatusCodeException(STATUS_CODE_FAILURE);
247 }
248 
249 //------------------------------------------------------------------------------------------------------------------------------------------
250 //------------------------------------------------------------------------------------------------------------------------------------------
251 
252 StatusCode OvershootTracksTool::ReadSettings(const TiXmlHandle xmlHandle)
253 {
254  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SplitMode", m_splitMode));
255 
256  PANDORA_RETURN_RESULT_IF_AND_IF(
257  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxVertexXSeparation", m_maxVertexXSeparation));
258 
259  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
260  XmlHelper::ReadValue(xmlHandle, "CosThetaCutForKinkSearch", m_cosThetaCutForKinkSearch));
261 
262  return ThreeDKinkBaseTool::ReadSettings(xmlHandle);
263 }
264 
265 } // namespace lar_content
pandora::ClusterList m_affectedClusters
The list of affected clusters.
std::vector< Modification > ModificationList
float GetXSamplingPoint(const pandora::CartesianVector &splitPosition1, const bool isForwardInX, const TwoDSlidingFitResult &fitResult1, const TwoDSlidingFitResult &fitResult2, const TwoDSlidingFitResult &fitResult3) const
Get a sampling point in x that is common to sliding linear fit objects in three views.
pandora::CartesianVector m_splitPosition2
The candidate split position in view 2.
bool IsThreeDKink(ThreeViewTransverseTracksAlgorithm *const pAlgorithm, const Particle &particle, const bool isA1LowestInX, const bool isA2LowestInX) const
Whether the provided particle is consistent with being a kink, when examined in three dimensions at t...
std::vector< TensorType::ElementList::const_iterator > IteratorList
void GetIteratorListModifications(ThreeViewTransverseTracksAlgorithm *const pAlgorithm, const IteratorList &iteratorList, ModificationList &modificationList) const
Get modification objects for a specific elements of the tensor, identifying required splits and merge...
bool m_splitMode
Whether to run in cluster splitting mode, as opposed to cluster merging mode.
const TwoDSlidingFitResult & GetCachedSlidingFitResult(const pandora::Cluster *const pCluster) const
Get a sliding fit result from the algorithm cache.
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
const pandora::Cluster * m_pCommonCluster
Address of the common cluster.
void SetSplitPosition(const LArPointingCluster::Vertex &vertexA1, const LArPointingCluster::Vertex &vertexA2, const LArPointingCluster::Vertex &vertexB1, const LArPointingCluster::Vertex &vertexB2, Particle &particle) const
Set split position for a provided particle.
intermediate_table::const_iterator const_iterator
LArPointingCluster class.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
ClusterMergeMap m_clusterMergeMap
The cluster merge map.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
TFile f
Definition: plotHisto.C:6
const pandora::Cluster * m_pClusterA1
Address of cluster in element A, view 1.
Header file for the geometry helper class.
float m_maxVertexXSeparation
The max separation between accompanying clusters vertex x positions to make split.
Header file for the cluster helper class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
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).
static void MergeThreePositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::HitType view3, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, const pandora::CartesianVector &position3, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from three views to give unified 3D position.
const pandora::Cluster * m_pClusterB2
Address of cluster in element B, view 2.
static bool IsALowestInX(const LArPointingCluster &pointingClusterA, const LArPointingCluster &pointingClusterB)
Whether pointing cluster labelled A extends to lowest x positions (as opposed to that labelled B) ...
Header file for the overshoot tracks tool class.
float m_cosThetaCutForKinkSearch
The cos theta cut used for the kink search in three dimensions.
pandora::CartesianVector m_splitPosition
The candidate split position for the common cluster.
Particle(const TensorType::Element &elementA, const TensorType::Element &elementB)
Constructor.
HitType
Definition: HitType.h:12
bool PassesVertexCuts(const LArPointingCluster::Vertex &vertexA, const LArPointingCluster::Vertex &vertexB) const
Whether a pair of vertices pass longitudinal projection cuts.
pandora::CartesianVector m_splitPosition1
The candidate split position in view 1.
ThreeDKinkBaseTool class.
static void GetClosestVerticesInX(const LArPointingCluster &pointingClusterI, const LArPointingCluster &pointingClusterJ, LArPointingCluster::Vertex &closestVertexI, LArPointingCluster::Vertex &closestVertexJ)
Given a pair of pointing clusters, find the pair of vertices with smallest x-separation.
bool m_majorityRulesMode
Whether to run in majority rules mode (always split overshoots, always merge undershoots) ...
float m_minLongitudinalImpactParameter
The min longitudinal impact parameter for connecting accompanying clusters.
const pandora::Cluster * m_pClusterA2
Address of cluster in element A, view 2.
const pandora::Cluster * m_pClusterB1
Address of cluster in element B, view 1.
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
SplitPositionMap m_splitPositionMap
The split position map.