LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArThreeDSlidingFitResult.cc
Go to the documentation of this file.
1 
9 #include "Helpers/ClusterFitHelper.h"
10 
11 #include "Objects/Cluster.h"
12 
15 
17 
18 #include <algorithm>
19 #include <cmath>
20 #include <limits>
21 
22 using namespace pandora;
23 
24 namespace lar_content
25 {
26 
27 template <typename T>
28 ThreeDSlidingFitResult::ThreeDSlidingFitResult(const T *const pT, const unsigned int layerWindow, const float layerPitch) :
29  m_primaryAxis(ThreeDSlidingFitResult::GetPrimaryAxis(pT, layerPitch)),
30  m_axisIntercept(m_primaryAxis.GetPosition()),
31  m_axisDirection(m_primaryAxis.GetMomentum()),
32  m_firstOrthoDirection(ThreeDSlidingFitResult::GetSeedDirection(m_axisDirection).GetCrossProduct(m_axisDirection).GetUnitVector()),
33  m_secondOrthoDirection(m_axisDirection.GetCrossProduct(m_firstOrthoDirection).GetUnitVector()),
34  m_firstFitResult(TwoDSlidingFitResult(pT, layerWindow, layerPitch, m_axisIntercept, m_axisDirection, m_firstOrthoDirection)),
35  m_secondFitResult(TwoDSlidingFitResult(pT, layerWindow, layerPitch, m_axisIntercept, m_axisDirection, m_secondOrthoDirection)),
36  m_minLayer(std::max(m_firstFitResult.GetMinLayer(), m_secondFitResult.GetMinLayer())),
37  m_maxLayer(std::min(m_firstFitResult.GetMaxLayer(), m_secondFitResult.GetMaxLayer())),
38  m_minLayerPosition(0.f, 0.f, 0.f),
39  m_maxLayerPosition(0.f, 0.f, 0.f),
40  m_minLayerDirection(0.f, 0.f, 0.f),
41  m_maxLayerDirection(0.f, 0.f, 0.f)
42 {
43  if (m_minLayer > m_maxLayer)
44  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
45 
46  const float minL(m_firstFitResult.GetL(m_minLayer));
47  const float maxL(m_firstFitResult.GetL(m_maxLayer));
48 
49  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetGlobalFitPosition(minL, m_minLayerPosition));
50  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetGlobalFitPosition(maxL, m_maxLayerPosition));
51  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetGlobalFitDirection(minL, m_minLayerDirection));
52  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetGlobalFitDirection(maxL, m_maxLayerDirection));
53 }
54 
55 //------------------------------------------------------------------------------------------------------------------------------------------
56 
57 const pandora::Cluster *ThreeDSlidingFitResult::GetCluster() const
58 {
60 }
61 
62 //------------------------------------------------------------------------------------------------------------------------------------------
63 
65 {
66  const int firstLayer(m_firstFitResult.GetMinLayer());
67  const int secondLayer(m_secondFitResult.GetMinLayer());
68 
69  return std::min(firstLayer, secondLayer);
70 }
71 
72 //------------------------------------------------------------------------------------------------------------------------------------------
73 
75 {
76  const int firstLayer(m_firstFitResult.GetMaxLayer());
77  const int secondLayer(m_secondFitResult.GetMaxLayer());
78 
79  return std::max(firstLayer, secondLayer);
80 }
81 
82 //------------------------------------------------------------------------------------------------------------------------------------------
83 
85 {
86  const float firstRms(m_firstFitResult.GetMinLayerRms());
87  const float secondRms(m_secondFitResult.GetMinLayerRms());
88 
89  return std::sqrt(firstRms * firstRms + secondRms * secondRms);
90 }
91 
92 //------------------------------------------------------------------------------------------------------------------------------------------
93 
95 {
96  const float firstRms(m_firstFitResult.GetMaxLayerRms());
97  const float secondRms(m_secondFitResult.GetMaxLayerRms());
98 
99  return std::sqrt(firstRms * firstRms + secondRms * secondRms);
100 }
101 
102 //------------------------------------------------------------------------------------------------------------------------------------------
103 
104 float ThreeDSlidingFitResult::GetFitRms(const float rL) const
105 {
106  const float firstRms(m_firstFitResult.GetFitRms(rL));
107  const float secondRms(m_secondFitResult.GetFitRms(rL));
108 
109  return std::sqrt(firstRms * firstRms + secondRms * secondRms);
110 }
111 
112 //------------------------------------------------------------------------------------------------------------------------------------------
113 
114 float ThreeDSlidingFitResult::GetLongitudinalDisplacement(const CartesianVector &position) const
115 {
116  return m_axisDirection.GetDotProduct(position - m_axisIntercept);
117 }
118 
119 //------------------------------------------------------------------------------------------------------------------------------------------
120 
121 StatusCode ThreeDSlidingFitResult::GetGlobalFitPosition(const float rL, CartesianVector &position) const
122 {
123  // Check that input coordinates are between first and last layers
124  const int layer1(m_firstFitResult.GetLayer(rL));
125  const int layer2(m_secondFitResult.GetLayer(rL));
126 
127  if (std::min(layer1, layer2) < m_minLayer || std::max(layer1, layer2) > m_maxLayer)
128  return STATUS_CODE_INVALID_PARAMETER;
129 
130  // Get local positions from each sliding fit (TODO: Make this more efficient)
131  CartesianVector firstPosition(0.f, 0.f, 0.f), secondPosition(0.f, 0.f, 0.f);
132  const StatusCode statusCode1(m_firstFitResult.GetGlobalFitPosition(rL, firstPosition));
133 
134  if (STATUS_CODE_SUCCESS != statusCode1)
135  return statusCode1;
136 
137  const StatusCode statusCode2(m_secondFitResult.GetGlobalFitPosition(rL, secondPosition));
138 
139  if (STATUS_CODE_SUCCESS != statusCode2)
140  return statusCode2;
141 
142  float rL1(0.f), rT1(0.f), rL2(0.f), rT2(0.f);
143  m_firstFitResult.GetLocalPosition(firstPosition, rL1, rT1);
144  m_secondFitResult.GetLocalPosition(secondPosition, rL2, rT2);
145 
146  // Combine local positions to give an overall global direction
147  this->GetGlobalPosition(rL, rT1, rT2, position);
148 
149  return STATUS_CODE_SUCCESS;
150 }
151 
152 //------------------------------------------------------------------------------------------------------------------------------------------
153 
154 StatusCode ThreeDSlidingFitResult::GetGlobalFitDirection(const float rL, CartesianVector &direction) const
155 {
156  // Check that input coordinates are between first and last layers
157  const int layer1(m_firstFitResult.GetLayer(rL));
158  const int layer2(m_secondFitResult.GetLayer(rL));
159 
160  if (std::min(layer1, layer2) < m_minLayer || std::max(layer1, layer2) > m_maxLayer)
161  return STATUS_CODE_INVALID_PARAMETER;
162 
163  // Get local directions from each sliding fit (TODO: Make this more efficient)
164  CartesianVector firstDirection(0.f, 0.f, 0.f), secondDirection(0.f, 0.f, 0.f);
165  const StatusCode statusCode1(m_firstFitResult.GetGlobalFitDirection(rL, firstDirection));
166 
167  if (STATUS_CODE_SUCCESS != statusCode1)
168  return statusCode1;
169 
170  const StatusCode statusCode2(m_secondFitResult.GetGlobalFitDirection(rL, secondDirection));
171 
172  if (STATUS_CODE_SUCCESS != statusCode2)
173  return statusCode2;
174 
175  float dTdL1(0.f), dTdL2(0.f);
176  m_firstFitResult.GetLocalDirection(firstDirection, dTdL1);
177  m_secondFitResult.GetLocalDirection(secondDirection, dTdL2);
178 
179  // Combine local directions to give an overall global direction
180  this->GetGlobalDirection(dTdL1, dTdL2, direction);
181 
182  return STATUS_CODE_SUCCESS;
183 }
184 
185 //------------------------------------------------------------------------------------------------------------------------------------------
186 
187 void ThreeDSlidingFitResult::GetGlobalPosition(const float rL, const float rT1, const float rT2, CartesianVector &position) const
188 {
190 }
191 
192 //------------------------------------------------------------------------------------------------------------------------------------------
193 
194 void ThreeDSlidingFitResult::GetGlobalDirection(const float dTdL1, const float dTdL2, CartesianVector &direction) const
195 {
196  const float pL(1.f / std::sqrt(1.f + dTdL1 * dTdL1 + dTdL2 * dTdL2));
197  const float pT1(dTdL1 / std::sqrt(1.f + dTdL1 * dTdL1 + dTdL2 * dTdL2));
198  const float pT2(dTdL2 / std::sqrt(1.f + dTdL1 * dTdL1 + dTdL2 * dTdL2));
199 
200  CartesianVector globalCoordinates(0.f, 0.f, 0.f);
201  this->GetGlobalPosition(pL, pT1, pT2, globalCoordinates);
202  direction = (globalCoordinates - m_axisIntercept).GetUnitVector();
203 }
204 
205 //------------------------------------------------------------------------------------------------------------------------------------------
206 
207 TrackState ThreeDSlidingFitResult::GetPrimaryAxis(const Cluster *const pCluster, const float layerPitch)
208 {
209  CartesianPointVector pointVector;
210  LArClusterHelper::GetCoordinateVector(pCluster, pointVector);
211  return ThreeDSlidingFitResult::GetPrimaryAxis(&pointVector, layerPitch);
212 }
213 
214 //------------------------------------------------------------------------------------------------------------------------------------------
215 
216 TrackState ThreeDSlidingFitResult::GetPrimaryAxis(const CartesianPointVector *const pPointVector, const float layerPitch)
217 {
218  if (pPointVector->size() < 2)
219  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
220 
221  CartesianVector centroid(0.f, 0.f, 0.f);
222  LArPcaHelper::EigenVectors eigenVecs;
223  LArPcaHelper::EigenValues eigenValues(0.f, 0.f, 0.f);
224  LArPcaHelper::RunPca(*pPointVector, centroid, eigenValues, eigenVecs);
225 
226  float minProjection(std::numeric_limits<float>::max());
227  CartesianVector fitDirection(eigenVecs.at(0));
228 
229  // By convention, the primary axis has a positive z-component.
230  // However, downstream algorithms should be independent of this convention.
231  if (fitDirection.GetZ() < 0.f)
232  fitDirection *= -1.f;
233 
234  for (const CartesianVector &coordinate : *pPointVector)
235  minProjection = std::min(minProjection, fitDirection.GetDotProduct(coordinate - centroid));
236 
237  // Define layers based on centroid rather than individual extremal hits
238  const float fitProjection(layerPitch * std::floor(minProjection / layerPitch));
239 
240  return TrackState(centroid + (fitDirection * fitProjection), fitDirection);
241 }
242 
243 //------------------------------------------------------------------------------------------------------------------------------------------
244 
245 CartesianVector ThreeDSlidingFitResult::GetSeedDirection(const CartesianVector &axisDirection)
246 {
247  const float px(std::fabs(axisDirection.GetX()));
248  const float py(std::fabs(axisDirection.GetY()));
249  const float pz(std::fabs(axisDirection.GetZ()));
250 
251  if (px < std::min(py, pz) + std::numeric_limits<float>::epsilon())
252  {
253  return CartesianVector(1.f, 0.f, 0.f);
254  }
255 
256  if (py < std::min(pz, px) + std::numeric_limits<float>::epsilon())
257  {
258  return CartesianVector(0.f, 1.f, 0.f);
259  }
260 
261  if (pz < std::min(px, py) + std::numeric_limits<float>::epsilon())
262  {
263  return CartesianVector(0.f, 0.f, 1.f);
264  }
265 
266  throw StatusCodeException(STATUS_CODE_FAILURE);
267 }
268 
269 //------------------------------------------------------------------------------------------------------------------------------------------
270 //------------------------------------------------------------------------------------------------------------------------------------------
271 
272 template ThreeDSlidingFitResult::ThreeDSlidingFitResult(const pandora::Cluster *const, const unsigned int, const float);
273 template ThreeDSlidingFitResult::ThreeDSlidingFitResult(const pandora::CartesianPointVector *const, const unsigned int, const float);
274 
275 } // namespace lar_content
const int m_maxLayer
The maximum combined layer.
static pandora::TrackState GetPrimaryAxis(const pandora::Cluster *const pCluster, const float slidingFitLayerPitch)
Calculate the position and direction of the primary axis.
pandora::CartesianVector EigenValues
Definition: LArPcaHelper.h:24
void GetLocalDirection(const pandora::CartesianVector &direction, float &dTdL) const
Get local sliding fit gradient for a given global direction.
float GetMaxLayerRms() const
Get rms at maximum layer.
int GetMaxLayer() const
Get the maximum occupied layer in the sliding fit.
static pandora::CartesianVector GetSeedDirection(const pandora::CartesianVector &axisDirection)
Generate a seed vector to be used in calculating the orthogonal axes.
float GetMaxLayerRms() const
Get rms at maximum layer.
STL namespace.
float GetMinLayerRms() const
Get rms at minimum layer.
Header file for the principal curve analysis helper class.
float GetLongitudinalDisplacement(const pandora::CartesianVector &position) const
Get longitudinal projection onto primary axis.
pandora::CartesianVector m_maxLayerPosition
The global position at the maximum combined layer.
void GetGlobalDirection(const float dTdL1, const float dTdL2, pandora::CartesianVector &direction) const
Get global direction coordinates for a given pair of sliding linear fit gradients.
TFile f
Definition: plotHisto.C:6
pandora::CartesianVector m_minLayerPosition
The global position at the minimum combined layer.
const int m_minLayer
The minimum combined layer.
int GetMaxLayer() const
Get the maximum occupied layer in the sliding fit.
const pandora::CartesianVector m_axisIntercept
The axis intercept position.
int GetMinLayer() const
Get the minimum occupied layer in the sliding fit.
const TwoDSlidingFitResult m_secondFitResult
The second sliding fit result.
Header file for the cluster helper class.
float GetFitRms(const float rL) const
Get fit rms for a given longitudinal coordinate.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
const TwoDSlidingFitResult m_firstFitResult
The first sliding fit result.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
const pandora::CartesianVector m_secondOrthoDirection
The orthogonal direction vector.
int GetMinLayer() const
Get the minimum occupied layer in the sliding fit.
const pandora::CartesianVector m_firstOrthoDirection
The orthogonal direction vector.
float GetMinLayerRms() const
Get rms at minimum layer.
Header file for the lar three dimensional sliding fit result class.
static void RunPca(const T &t, pandora::CartesianVector &centroid, EigenValues &outputEigenValues, EigenVectors &outputEigenVectors)
Run principal component analysis using input calo hits (TPC_VIEW_U,V,W or TPC_3D; all treated as 3D p...
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
void GetGlobalPosition(const float rL, const float rT1, const float rT2, pandora::CartesianVector &position) const
Get global coordinates for a given pair of sliding linear fit coordinates.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
float GetFitRms(const float rL) const
Get fit rms for a given longitudinal coordinate.
std::vector< pandora::CartesianVector > EigenVectors
Definition: LArPcaHelper.h:25
float GetL(const int layer) const
Get longitudinal coordinate for a given sliding linear fit layer number.
pandora::CartesianVector m_minLayerDirection
The global direction at the minimum combined layer.
static void GetCoordinateVector(const pandora::Cluster *const pCluster, pandora::CartesianPointVector &coordinateVector)
Get vector of hit coordinates from an input cluster.
pandora::CartesianVector m_maxLayerDirection
The global direction at the maximum combined layer.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
const pandora::CartesianVector m_axisDirection
The axis direction vector.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
ThreeDSlidingFitResult(const T *const pT, const unsigned int slidingFitWindow, const float slidingFitLayerPitch)
Constructor.
int GetLayer(const float rL) const
Get layer number for given sliding linear fit longitudinal coordinate.