LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ThreeViewLongitudinalTracksAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
13 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 ThreeViewLongitudinalTracksAlgorithm::ThreeViewLongitudinalTracksAlgorithm() :
22  m_nMaxTensorToolRepeats(1000),
23  m_vertexChi2Cut(10.f),
24  m_reducedChi2Cut(5.f),
25  m_samplingPitch(1.f)
26 {
27 }
28 
29 //------------------------------------------------------------------------------------------------------------------------------------------
30 
32  const Cluster *const pClusterU, const Cluster *const pClusterV, const Cluster *const pClusterW)
33 {
34  LongitudinalOverlapResult overlapResult;
35  this->CalculateOverlapResult(pClusterU, pClusterV, pClusterW, overlapResult);
36 
37  if (overlapResult.IsInitialized())
38  this->GetMatchingControl().GetOverlapTensor().SetOverlapResult(pClusterU, pClusterV, pClusterW, overlapResult);
39 }
40 
41 //------------------------------------------------------------------------------------------------------------------------------------------
42 
43 void ThreeViewLongitudinalTracksAlgorithm::CalculateOverlapResult(const Cluster *const pClusterU, const Cluster *const pClusterV,
44  const Cluster *const pClusterW, LongitudinalOverlapResult &longitudinalOverlapResult)
45 {
46  const TwoDSlidingFitResult &slidingFitResultU(this->GetCachedSlidingFitResult(pClusterU));
47  const TwoDSlidingFitResult &slidingFitResultV(this->GetCachedSlidingFitResult(pClusterV));
48  const TwoDSlidingFitResult &slidingFitResultW(this->GetCachedSlidingFitResult(pClusterW));
49 
50  // Loop over possible permutations of cluster direction
51  TrackOverlapResult bestOverlapResult;
52 
53  for (unsigned int iPermutation = 0; iPermutation < 4; ++iPermutation)
54  {
55  const bool isForwardU((1 == iPermutation) ? false : true);
56  const bool isForwardV((2 == iPermutation) ? false : true);
57  const bool isForwardW((3 == iPermutation) ? false : true);
58 
59  // Get 2D start and end positions from each sliding fit for this permutation
60  const CartesianVector vtxU((isForwardU) ? slidingFitResultU.GetGlobalMinLayerPosition() : slidingFitResultU.GetGlobalMaxLayerPosition());
61  const CartesianVector endU((!isForwardU) ? slidingFitResultU.GetGlobalMinLayerPosition() : slidingFitResultU.GetGlobalMaxLayerPosition());
62 
63  const CartesianVector vtxV((isForwardV) ? slidingFitResultV.GetGlobalMinLayerPosition() : slidingFitResultV.GetGlobalMaxLayerPosition());
64  const CartesianVector endV((!isForwardV) ? slidingFitResultV.GetGlobalMinLayerPosition() : slidingFitResultV.GetGlobalMaxLayerPosition());
65 
66  const CartesianVector vtxW((isForwardW) ? slidingFitResultW.GetGlobalMinLayerPosition() : slidingFitResultW.GetGlobalMaxLayerPosition());
67  const CartesianVector endW((!isForwardW) ? slidingFitResultW.GetGlobalMinLayerPosition() : slidingFitResultW.GetGlobalMaxLayerPosition());
68 
69  // Merge start and end positions (three views)
70  const float halfLengthSquaredU(0.25 * (vtxU - endU).GetMagnitudeSquared());
71  const float halfLengthSquaredV(0.25 * (vtxV - endV).GetMagnitudeSquared());
72  const float halfLengthSquaredW(0.25 * (vtxW - endW).GetMagnitudeSquared());
73 
74  CartesianVector vtxMergedU(0.f, 0.f, 0.f), vtxMergedV(0.f, 0.f, 0.f), vtxMergedW(0.f, 0.f, 0.f);
75  CartesianVector endMergedU(0.f, 0.f, 0.f), endMergedV(0.f, 0.f, 0.f), endMergedW(0.f, 0.f, 0.f);
76 
77  float vtxChi2(std::numeric_limits<float>::max());
78  float endChi2(std::numeric_limits<float>::max());
79 
81  this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W, vtxU, vtxV, vtxW, vtxMergedU, vtxMergedV, vtxMergedW, vtxChi2);
82 
83  if (vtxChi2 > m_vertexChi2Cut)
84  continue;
85 
86  if (((vtxMergedU - vtxU).GetMagnitudeSquared() > std::min(halfLengthSquaredU, (vtxMergedU - endU).GetMagnitudeSquared())) ||
87  ((vtxMergedV - vtxV).GetMagnitudeSquared() > std::min(halfLengthSquaredV, (vtxMergedV - endV).GetMagnitudeSquared())) ||
88  ((vtxMergedW - vtxW).GetMagnitudeSquared() > std::min(halfLengthSquaredW, (vtxMergedW - endW).GetMagnitudeSquared())))
89  continue;
90 
92  this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W, endU, endV, endW, endMergedU, endMergedV, endMergedW, endChi2);
93 
94  if (endChi2 > m_vertexChi2Cut)
95  continue;
96 
97  if (((endMergedU - endU).GetMagnitudeSquared() > std::min(halfLengthSquaredU, (endMergedU - vtxU).GetMagnitudeSquared())) ||
98  ((endMergedV - endV).GetMagnitudeSquared() > std::min(halfLengthSquaredV, (endMergedV - vtxV).GetMagnitudeSquared())) ||
99  ((endMergedW - endW).GetMagnitudeSquared() > std::min(halfLengthSquaredW, (endMergedW - vtxW).GetMagnitudeSquared())))
100  continue;
101 
102  // Merge start and end positions (two views)
103  float chi2(0.f);
104  CartesianVector position3D(0.f, 0.f, 0.f);
105  CartesianPointVector vtxList3D, endList3D;
106 
107  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, vtxU, vtxV, position3D, chi2);
108  vtxList3D.push_back(position3D);
109 
110  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_V, TPC_VIEW_W, vtxV, vtxW, position3D, chi2);
111  vtxList3D.push_back(position3D);
112 
113  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_W, TPC_VIEW_U, vtxW, vtxU, position3D, chi2);
114  vtxList3D.push_back(position3D);
115 
116  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, endU, endV, position3D, chi2);
117  endList3D.push_back(position3D);
118 
119  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_V, TPC_VIEW_W, endV, endW, position3D, chi2);
120  endList3D.push_back(position3D);
121 
122  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), TPC_VIEW_W, TPC_VIEW_U, endW, endU, position3D, chi2);
123  endList3D.push_back(position3D);
124 
125  // Find best matched 3D trajactory
126  for (CartesianPointVector::const_iterator iterI = vtxList3D.begin(), iterEndI = vtxList3D.end(); iterI != iterEndI; ++iterI)
127  {
128  const CartesianVector &vtxMerged3D(*iterI);
129 
130  for (CartesianPointVector::const_iterator iterJ = endList3D.begin(), iterEndJ = endList3D.end(); iterJ != iterEndJ; ++iterJ)
131  {
132  const CartesianVector &endMerged3D(*iterJ);
133 
134  TrackOverlapResult overlapResult;
135  this->CalculateOverlapResult(slidingFitResultU, slidingFitResultV, slidingFitResultW, vtxMerged3D, endMerged3D, overlapResult);
136 
137  if (overlapResult.IsInitialized() && (overlapResult.GetNMatchedSamplingPoints() > 0) && (overlapResult > bestOverlapResult))
138  {
139  bestOverlapResult = overlapResult;
140  longitudinalOverlapResult = LongitudinalOverlapResult(overlapResult, vtxChi2, endChi2);
141  }
142  }
143  }
144  }
145 }
146 
147 //------------------------------------------------------------------------------------------------------------------------------------------
148 
150  const TwoDSlidingFitResult &slidingFitResultV, const TwoDSlidingFitResult &slidingFitResultW, const CartesianVector &vtxMerged3D,
151  const CartesianVector &endMerged3D, TrackOverlapResult &overlapResult) const
152 {
153  // Calculate start and end positions of linear trajectory
154  const CartesianVector vtxMergedU(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtxMerged3D, TPC_VIEW_U));
155  const CartesianVector vtxMergedV(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtxMerged3D, TPC_VIEW_V));
156  const CartesianVector vtxMergedW(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtxMerged3D, TPC_VIEW_W));
157 
158  const CartesianVector endMergedU(LArGeometryHelper::ProjectPosition(this->GetPandora(), endMerged3D, TPC_VIEW_U));
159  const CartesianVector endMergedV(LArGeometryHelper::ProjectPosition(this->GetPandora(), endMerged3D, TPC_VIEW_V));
160  const CartesianVector endMergedW(LArGeometryHelper::ProjectPosition(this->GetPandora(), endMerged3D, TPC_VIEW_W));
161 
162  const unsigned int nSamplingPoints = static_cast<unsigned int>((endMerged3D - vtxMerged3D).GetMagnitude() / m_samplingPitch);
163 
164  if (0 == nSamplingPoints)
165  return;
166 
167  // Loop over sampling points and calculate track overlap result
168  float deltaChi2(0.f), totalChi2(0.f);
169  unsigned int nMatchedSamplingPoints(0);
170 
171  for (unsigned int n = 0; n < nSamplingPoints; ++n)
172  {
173  const float alpha((0.5f + static_cast<float>(n)) / static_cast<float>(nSamplingPoints));
174  const CartesianVector linearU(vtxMergedU + (endMergedU - vtxMergedU) * alpha);
175  const CartesianVector linearV(vtxMergedV + (endMergedV - vtxMergedV) * alpha);
176  const CartesianVector linearW(vtxMergedW + (endMergedW - vtxMergedW) * alpha);
177 
178  CartesianVector posU(0.f, 0.f, 0.f), posV(0.f, 0.f, 0.f), posW(0.f, 0.f, 0.f);
179  if ((STATUS_CODE_SUCCESS != slidingFitResultU.GetGlobalFitProjection(linearU, posU)) ||
180  (STATUS_CODE_SUCCESS != slidingFitResultV.GetGlobalFitProjection(linearV, posV)) ||
181  (STATUS_CODE_SUCCESS != slidingFitResultW.GetGlobalFitProjection(linearW, posW)))
182  {
183  continue;
184  }
185 
186  CartesianVector mergedU(0.f, 0.f, 0.f), mergedV(0.f, 0.f, 0.f), mergedW(0.f, 0.f, 0.f);
187  LArGeometryHelper::MergeThreePositions(this->GetPandora(), posU, posV, posW, mergedU, mergedV, mergedW, deltaChi2);
188 
189  if (deltaChi2 < m_reducedChi2Cut)
190  ++nMatchedSamplingPoints;
191 
192  totalChi2 += deltaChi2;
193  }
194 
195  if (nMatchedSamplingPoints > 0)
196  overlapResult = TrackOverlapResult(nMatchedSamplingPoints, nSamplingPoints, totalChi2);
197 }
198 
199 //------------------------------------------------------------------------------------------------------------------------------------------
200 
202 {
203  unsigned int repeatCounter(0);
204 
205  for (TensorToolVector::const_iterator iter = m_algorithmToolVector.begin(), iterEnd = m_algorithmToolVector.end(); iter != iterEnd;)
206  {
207  if ((*iter)->Run(this, this->GetMatchingControl().GetOverlapTensor()))
208  {
209  iter = m_algorithmToolVector.begin();
210 
211  if (++repeatCounter > m_nMaxTensorToolRepeats)
212  break;
213  }
214  else
215  {
216  ++iter;
217  }
218  }
219 }
220 
221 //------------------------------------------------------------------------------------------------------------------------------------------
222 
223 StatusCode ThreeViewLongitudinalTracksAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
224 {
225  AlgorithmToolVector algorithmToolVector;
226  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*this, xmlHandle, "TrackTools", algorithmToolVector));
227 
228  for (AlgorithmToolVector::const_iterator iter = algorithmToolVector.begin(), iterEnd = algorithmToolVector.end(); iter != iterEnd; ++iter)
229  {
230  LongitudinalTensorTool *const pLongitudinalTensorTool(dynamic_cast<LongitudinalTensorTool *>(*iter));
231 
232  if (!pLongitudinalTensorTool)
233  return STATUS_CODE_INVALID_PARAMETER;
234 
235  m_algorithmToolVector.push_back(pLongitudinalTensorTool);
236  }
237 
238  PANDORA_RETURN_RESULT_IF_AND_IF(
239  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NMaxTensorToolRepeats", m_nMaxTensorToolRepeats));
240 
241  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "VertexChi2Cut", m_vertexChi2Cut));
242 
243  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ReducedChi2Cut", m_reducedChi2Cut));
244 
245  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SamplingPitch", m_samplingPitch));
246 
247  if (m_samplingPitch < std::numeric_limits<float>::epsilon())
248  return STATUS_CODE_INVALID_PARAMETER;
249 
250  return BaseAlgorithm::ReadSettings(xmlHandle);
251 }
252 
253 } // namespace lar_content
LongitudinalOverlapResult class.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static void MergeTwoPositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from two views to give unified 3D position.
TensorType & GetOverlapTensor()
Get the overlap tensor.
unsigned int GetNMatchedSamplingPoints() const
Get the number of matched sampling points.
void ExamineOverlapContainer()
Examine contents of overlap container, collect together best-matching 2D particles and modify cluster...
bool IsInitialized() const
Whether the track overlap result has been initialized.
unsigned int m_nMaxTensorToolRepeats
The maximum number of repeat loops over tensor tools.
const TwoDSlidingFitResult & GetCachedSlidingFitResult(const pandora::Cluster *const pCluster) const
Get a sliding fit result from the algorithm cache.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
intermediate_table::const_iterator const_iterator
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
void SetOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW, const OverlapResult &overlapResult)
Set overlap result.
Header file for the cluster helper class.
void CalculateOverlapResult(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW)
Calculate cluster overlap result and store in container.
float m_samplingPitch
Pitch used to generate sampling points along tracks.
TensorToolVector m_algorithmToolVector
The algorithm tool vector.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
static void MergeThreePositions(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 &outputU, pandora::CartesianVector &outputV, pandora::CartesianVector &outputW, float &chiSquared)
Merge 2D positions from three views to give unified 2D positions for each view.
TrackOverlapResult class.
float m_vertexChi2Cut
The maximum allowed chi2 for associating end points from three views.
pandora::StatusCode GetGlobalFitProjection(const pandora::CartesianVector &inputPosition, pandora::CartesianVector &projectedPosition) const
Get projected position on global fit for a given position vector.
Header file for the three view longitudinal tracks algorithm class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Char_t n[5]
float m_reducedChi2Cut
The maximum allowed chi2 for associating hit positions from three views.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.