LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArTwoDSlidingFitResult.cc
Go to the documentation of this file.
1 
13 
15 
16 #include <algorithm>
17 #include <cmath>
18 #include <limits>
19 
20 using namespace pandora;
21 
22 namespace lar_content
23 {
24 
25 template <>
26 TwoDSlidingFitResult::TwoDSlidingFitResult(const Cluster *const pCluster, const unsigned int layerFitHalfWindow, const float layerPitch,
27  const float axisDeviationLimitForHitDivision) :
28  m_pCluster(pCluster),
29  m_layerFitHalfWindow(layerFitHalfWindow),
30  m_layerPitch(layerPitch),
31  m_axisIntercept(0.f, 0.f, 0.f),
32  m_axisDirection(0.f, 0.f, 0.f),
33  m_orthoDirection(0.f, 0.f, 0.f)
34 {
35  CartesianPointVector pointVector;
36  LArClusterHelper::GetCoordinateVector(pCluster, pointVector);
37 
38  this->CalculateAxes(pointVector, layerPitch);
39 
40  const CartesianVector xAxis(1.f, 0.f, 0.f);
41  const float cosOpeningAngle(xAxis.GetCosOpeningAngle(m_axisDirection));
42 
43  if (std::fabs(cosOpeningAngle) < axisDeviationLimitForHitDivision)
44  {
45  this->FillLayerFitContributionMap(pointVector);
46  }
47  else
48  {
49  // TODO Refactor hit splitting and ensure all parameters configurable
50  LArHitWidthHelper::ConstituentHitVector constituentHitVector(LArHitWidthHelper::GetConstituentHits(pCluster, 0.5f, 1.f, true));
51  const CartesianPointVector constituentHitPointVector(LArHitWidthHelper::GetConstituentHitPositionVector(constituentHitVector));
52  this->FillLayerFitContributionMap(constituentHitPointVector);
53  }
54 
56  this->FindSlidingFitSegments();
57 }
58 
59 template <>
60 TwoDSlidingFitResult::TwoDSlidingFitResult(const CartesianPointVector *const pPointVector, const unsigned int layerFitHalfWindow,
61  const float layerPitch, const float /*axisDeviationLimitForHitDivision*/) :
62  m_pCluster(nullptr),
63  m_layerFitHalfWindow(layerFitHalfWindow),
64  m_layerPitch(layerPitch),
65  m_axisIntercept(0.f, 0.f, 0.f),
66  m_axisDirection(0.f, 0.f, 0.f),
67  m_orthoDirection(0.f, 0.f, 0.f)
68 {
69  this->CalculateAxes(*pPointVector, layerPitch);
70  this->FillLayerFitContributionMap(*pPointVector);
72  this->FindSlidingFitSegments();
73 }
74 
75 //------------------------------------------------------------------------------------------------------------------------------------------
76 
77 template <>
78 TwoDSlidingFitResult::TwoDSlidingFitResult(const Cluster *const pCluster, const unsigned int layerFitHalfWindow, const float layerPitch,
79  const CartesianVector &axisIntercept, const CartesianVector &axisDirection, const CartesianVector &orthoDirection,
80  const float axisDeviationLimitForHitDivision) :
81  m_pCluster(pCluster),
82  m_layerFitHalfWindow(layerFitHalfWindow),
83  m_layerPitch(layerPitch),
84  m_axisIntercept(axisIntercept),
85  m_axisDirection(axisDirection),
86  m_orthoDirection(orthoDirection)
87 {
88  const CartesianVector xAxis(1.f, 0.f, 0.f);
89  const float cosOpeningAngle(xAxis.GetCosOpeningAngle(m_axisDirection));
90 
91  if (std::fabs(cosOpeningAngle) < axisDeviationLimitForHitDivision)
92  {
93  CartesianPointVector pointVector;
94  LArClusterHelper::GetCoordinateVector(pCluster, pointVector);
95  this->FillLayerFitContributionMap(pointVector);
96  }
97  else
98  {
99  // TODO Refactor hit splitting and ensure all parameters configurable
100  LArHitWidthHelper::ConstituentHitVector constituentHitVector(LArHitWidthHelper::GetConstituentHits(pCluster, 0.5f, 1.f, true));
101  const CartesianPointVector constituentHitPointVector(LArHitWidthHelper::GetConstituentHitPositionVector(constituentHitVector));
102  this->FillLayerFitContributionMap(constituentHitPointVector);
103  }
104 
105  this->PerformSlidingLinearFit();
106  this->FindSlidingFitSegments();
107 }
108 
109 template <>
110 TwoDSlidingFitResult::TwoDSlidingFitResult(const CartesianPointVector *const pPointVector, const unsigned int layerFitHalfWindow,
111  const float layerPitch, const CartesianVector &axisIntercept, const CartesianVector &axisDirection,
112  const CartesianVector &orthoDirection, const float /*axisDeviationLimitForHitDivision*/) :
113  m_pCluster(nullptr),
114  m_layerFitHalfWindow(layerFitHalfWindow),
115  m_layerPitch(layerPitch),
116  m_axisIntercept(axisIntercept),
117  m_axisDirection(axisDirection),
118  m_orthoDirection(orthoDirection)
119 {
120  this->FillLayerFitContributionMap(*pPointVector);
121  this->PerformSlidingLinearFit();
122  this->FindSlidingFitSegments();
123 }
124 
125 //------------------------------------------------------------------------------------------------------------------------------------------
126 
127 TwoDSlidingFitResult::TwoDSlidingFitResult(const unsigned int layerFitHalfWindow, const float layerPitch, const CartesianVector &axisIntercept,
128  const CartesianVector &axisDirection, const CartesianVector &orthoDirection, const LayerFitContributionMap &layerFitContributionMap) :
129  m_pCluster(nullptr),
130  m_layerFitHalfWindow(layerFitHalfWindow),
131  m_layerPitch(layerPitch),
132  m_axisIntercept(axisIntercept),
133  m_axisDirection(axisDirection),
134  m_orthoDirection(orthoDirection),
135  m_layerFitContributionMap(layerFitContributionMap)
136 {
137  this->PerformSlidingLinearFit();
138  this->FindSlidingFitSegments();
139 }
140 
141 //------------------------------------------------------------------------------------------------------------------------------------------
142 
143 const pandora::Cluster *TwoDSlidingFitResult::GetCluster() const
144 {
145  if (!m_pCluster)
146  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
147 
148  return m_pCluster;
149 }
150 
151 //------------------------------------------------------------------------------------------------------------------------------------------
152 
154 {
155  return (static_cast<float>(m_layerFitHalfWindow)) * m_layerPitch;
156 }
157 
158 //------------------------------------------------------------------------------------------------------------------------------------------
159 
161 {
162  if (m_layerFitResultMap.empty())
163  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
164 
165  return m_layerFitResultMap.begin()->first;
166 }
167 
168 //------------------------------------------------------------------------------------------------------------------------------------------
169 
171 {
172  if (m_layerFitResultMap.empty())
173  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
174 
175  return m_layerFitResultMap.rbegin()->first;
176 }
177 
178 //------------------------------------------------------------------------------------------------------------------------------------------
179 
180 int TwoDSlidingFitResult::GetLayer(const float rL) const
181 {
182  if (m_layerPitch < std::numeric_limits<float>::epsilon())
183  throw StatusCodeException(STATUS_CODE_FAILURE);
184 
185  return static_cast<int>(std::floor(rL / m_layerPitch));
186 }
187 
188 //------------------------------------------------------------------------------------------------------------------------------------------
189 
190 float TwoDSlidingFitResult::GetL(const int layer) const
191 {
192  return (static_cast<float>(layer) + 0.5f) * m_layerPitch;
193 }
194 
195 //------------------------------------------------------------------------------------------------------------------------------------------
196 
197 void TwoDSlidingFitResult::GetLocalPosition(const CartesianVector &position, float &rL, float &rT) const
198 {
199  const CartesianVector displacement(position - m_axisIntercept);
200 
201  rL = displacement.GetDotProduct(m_axisDirection);
202  rT = displacement.GetDotProduct(m_orthoDirection);
203 }
204 
205 //------------------------------------------------------------------------------------------------------------------------------------------
206 
207 void TwoDSlidingFitResult::GetLocalDirection(const CartesianVector &direction, float &dTdL) const
208 {
209  float pL(0.f), pT(0.f);
210  this->GetLocalPosition((direction + m_axisIntercept), pL, pT);
211 
212  if (std::fabs(pL) < std::numeric_limits<float>::epsilon())
213  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
214 
215  dTdL = pT / pL;
216 }
217 
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
220 void TwoDSlidingFitResult::GetGlobalPosition(const float rL, const float rT, CartesianVector &position) const
221 {
222  position = m_axisIntercept + (m_axisDirection * rL) + (m_orthoDirection * rT);
223 }
224 
225 //------------------------------------------------------------------------------------------------------------------------------------------
226 
227 void TwoDSlidingFitResult::GetGlobalDirection(const float dTdL, CartesianVector &direction) const
228 {
229  const float pL(1.f / std::sqrt(1.f + dTdL * dTdL));
230  const float pT(dTdL / std::sqrt(1.f + dTdL * dTdL));
231 
232  CartesianVector globalCoordinates(0.f, 0.f, 0.f);
233  this->GetGlobalPosition(pL, pT, globalCoordinates);
234  direction = (globalCoordinates - m_axisIntercept).GetUnitVector();
235 }
236 
237 //------------------------------------------------------------------------------------------------------------------------------------------
238 
240 {
241  if (m_layerFitResultMap.empty())
242  throw StatusCodeException(STATUS_CODE_FAILURE);
243 
245  CartesianVector position(0.f, 0.f, 0.f);
246  this->GetGlobalPosition(iter->second.GetL(), iter->second.GetFitT(), position);
247  return position;
248 }
249 
250 //------------------------------------------------------------------------------------------------------------------------------------------
251 
253 {
254  if (m_layerFitResultMap.empty())
255  throw StatusCodeException(STATUS_CODE_FAILURE);
256 
257  LayerFitResultMap::const_reverse_iterator iter = m_layerFitResultMap.rbegin();
258  CartesianVector position(0.f, 0.f, 0.f);
259  this->GetGlobalPosition(iter->second.GetL(), iter->second.GetFitT(), position);
260  return position;
261 }
262 
263 //------------------------------------------------------------------------------------------------------------------------------------------
264 
266 {
267  if (m_layerFitResultMap.empty())
268  throw StatusCodeException(STATUS_CODE_FAILURE);
269 
271  CartesianVector direction(0.f, 0.f, 0.f);
272  this->GetGlobalDirection(iter->second.GetGradient(), direction);
273  return direction;
274 }
275 
276 //------------------------------------------------------------------------------------------------------------------------------------------
277 
279 {
280  if (m_layerFitResultMap.empty())
281  throw StatusCodeException(STATUS_CODE_FAILURE);
282 
283  LayerFitResultMap::const_reverse_iterator iter = m_layerFitResultMap.rbegin();
284  CartesianVector direction(0.f, 0.f, 0.f);
285  this->GetGlobalDirection(iter->second.GetGradient(), direction);
286  return direction;
287 }
288 
289 //------------------------------------------------------------------------------------------------------------------------------------------
290 
292 {
293  if (m_layerFitResultMap.empty())
294  throw StatusCodeException(STATUS_CODE_FAILURE);
295 
297  return iter->second.GetRms();
298 }
299 
300 //------------------------------------------------------------------------------------------------------------------------------------------
301 
303 {
304  if (m_layerFitResultMap.empty())
305  throw StatusCodeException(STATUS_CODE_FAILURE);
306 
307  LayerFitResultMap::const_reverse_iterator iter = m_layerFitResultMap.rbegin();
308  return iter->second.GetRms();
309 }
310 
311 //------------------------------------------------------------------------------------------------------------------------------------------
312 
313 float TwoDSlidingFitResult::GetFitRms(const float rL) const
314 {
315  LayerInterpolation layerInterpolation;
316  const StatusCode statusCode(this->LongitudinalInterpolation(rL, layerInterpolation));
317 
318  if (STATUS_CODE_SUCCESS != statusCode)
319  throw StatusCodeException(statusCode);
320 
321  return this->GetFitRms(layerInterpolation);
322 }
323 
324 //------------------------------------------------------------------------------------------------------------------------------------------
325 
327 {
328  const float halfWindowLength(this->GetLayerFitHalfWindowLength());
329 
330  CartesianVector firstDirection(0.f, 0.f, 0.f);
331  CartesianVector secondDirection(0.f, 0.f, 0.f);
332 
333  const StatusCode firstStatusCode(this->GetGlobalFitDirection(rL - halfWindowLength, firstDirection));
334  const StatusCode secondStatusCode(this->GetGlobalFitDirection(rL + halfWindowLength, secondDirection));
335 
336  if (STATUS_CODE_SUCCESS != firstStatusCode)
337  throw StatusCodeException(firstStatusCode);
338 
339  if (STATUS_CODE_SUCCESS != secondStatusCode)
340  throw StatusCodeException(secondStatusCode);
341 
342  return firstDirection.GetDotProduct(secondDirection);
343 }
344 
345 //------------------------------------------------------------------------------------------------------------------------------------------
346 
347 StatusCode TwoDSlidingFitResult::GetGlobalFitPosition(const float rL, CartesianVector &position) const
348 {
349  LayerInterpolation layerInterpolation;
350  const StatusCode statusCode(this->LongitudinalInterpolation(rL, layerInterpolation));
351 
352  if (STATUS_CODE_SUCCESS != statusCode)
353  return statusCode;
354 
355  position = this->GetGlobalFitPosition(layerInterpolation);
356  return STATUS_CODE_SUCCESS;
357 }
358 
359 //------------------------------------------------------------------------------------------------------------------------------------------
360 
361 StatusCode TwoDSlidingFitResult::GetGlobalFitDirection(const float rL, CartesianVector &direction) const
362 {
363  LayerInterpolation layerInterpolation;
364  const StatusCode statusCode(this->LongitudinalInterpolation(rL, layerInterpolation));
365 
366  if (STATUS_CODE_SUCCESS != statusCode)
367  return statusCode;
368 
369  direction = this->GetGlobalFitDirection(layerInterpolation);
370  return STATUS_CODE_SUCCESS;
371 }
372 
373 //------------------------------------------------------------------------------------------------------------------------------------------
374 
375 StatusCode TwoDSlidingFitResult::GetGlobalFitPositionAtX(const float x, CartesianVector &position) const
376 {
377  LayerInterpolationList layerInterpolationList;
378  const StatusCode statusCode(this->TransverseInterpolation(x, layerInterpolationList));
379 
380  if (STATUS_CODE_SUCCESS != statusCode)
381  return statusCode;
382 
383  if (layerInterpolationList.size() != 1)
384  return STATUS_CODE_NOT_FOUND;
385 
386  position = this->GetGlobalFitPosition(layerInterpolationList.at(0));
387  return STATUS_CODE_SUCCESS;
388 }
389 
390 //------------------------------------------------------------------------------------------------------------------------------------------
391 
392 StatusCode TwoDSlidingFitResult::GetGlobalFitDirectionAtX(const float x, CartesianVector &direction) const
393 {
394  LayerInterpolationList layerInterpolationList;
395  const StatusCode statusCode(this->TransverseInterpolation(x, layerInterpolationList));
396 
397  if (STATUS_CODE_SUCCESS != statusCode)
398  return statusCode;
399 
400  if (layerInterpolationList.size() != 1)
401  return STATUS_CODE_NOT_FOUND;
402 
403  direction = this->GetGlobalFitDirection(layerInterpolationList.at(0));
404  return STATUS_CODE_SUCCESS;
405 }
406 
407 //------------------------------------------------------------------------------------------------------------------------------------------
408 
409 StatusCode TwoDSlidingFitResult::GetGlobalFitProjection(const CartesianVector &inputPosition, CartesianVector &projectedPosition) const
410 {
411  float rL(0.f), rT(0.f);
412  this->GetLocalPosition(inputPosition, rL, rT);
413  return this->GetGlobalFitPosition(rL, projectedPosition);
414 }
415 
416 //------------------------------------------------------------------------------------------------------------------------------------------
417 
418 StatusCode TwoDSlidingFitResult::GetGlobalFitPositionListAtX(const float x, CartesianPointVector &positionList) const
419 {
420  LayerInterpolationList layerInterpolationList;
421  const StatusCode statusCode(this->TransverseInterpolation(x, layerInterpolationList));
422 
423  if (STATUS_CODE_SUCCESS != statusCode)
424  return statusCode;
425 
426  for (LayerInterpolationList::const_iterator iter = layerInterpolationList.begin(), iterEnd = layerInterpolationList.end(); iter != iterEnd; ++iter)
427  {
428  positionList.push_back(this->GetGlobalFitPosition(*iter));
429  }
430 
431  return STATUS_CODE_SUCCESS;
432 }
433 
434 //------------------------------------------------------------------------------------------------------------------------------------------
435 
436 StatusCode TwoDSlidingFitResult::GetTransverseProjection(const float x, const FitSegment &fitSegment, CartesianVector &position) const
437 {
438  LayerInterpolation layerInterpolation;
439  const StatusCode statusCode(this->TransverseInterpolation(x, fitSegment, layerInterpolation));
440 
441  if (STATUS_CODE_SUCCESS != statusCode)
442  return statusCode;
443 
444  position = this->GetGlobalFitPosition(layerInterpolation);
445  return STATUS_CODE_SUCCESS;
446 }
447 
448 //------------------------------------------------------------------------------------------------------------------------------------------
449 
451  const float x, const FitSegment &fitSegment, CartesianVector &position, CartesianVector &direction) const
452 {
453  LayerInterpolation layerInterpolation;
454  const StatusCode statusCode(this->TransverseInterpolation(x, fitSegment, layerInterpolation));
455 
456  if (STATUS_CODE_SUCCESS != statusCode)
457  return statusCode;
458 
459  position = this->GetGlobalFitPosition(layerInterpolation);
460  direction = this->GetGlobalFitDirection(layerInterpolation);
461  return STATUS_CODE_SUCCESS;
462 }
463 
464 //------------------------------------------------------------------------------------------------------------------------------------------
465 
466 StatusCode TwoDSlidingFitResult::GetExtrapolatedPosition(const float rL, CartesianVector &position) const
467 {
468  const StatusCode statusCode(this->GetGlobalFitPosition(rL, position));
469 
470  if (STATUS_CODE_NOT_FOUND != statusCode)
471  return statusCode;
472 
473  const int thisLayer(this->GetLayer(rL));
474  const int minLayer(this->GetMinLayer());
475  const int maxLayer(this->GetMaxLayer());
476 
477  if (thisLayer <= minLayer)
478  {
479  position = (this->GetGlobalMinLayerPosition() + this->GetGlobalMinLayerDirection() * (rL - this->GetL(minLayer)));
480  }
481  else if (thisLayer >= maxLayer)
482  {
483  position = (this->GetGlobalMaxLayerPosition() + this->GetGlobalMaxLayerDirection() * (rL - this->GetL(maxLayer)));
484  }
485  else
486  {
487  return STATUS_CODE_FAILURE;
488  }
489 
490  return STATUS_CODE_SUCCESS;
491 }
492 
493 //------------------------------------------------------------------------------------------------------------------------------------------
494 
495 StatusCode TwoDSlidingFitResult::GetExtrapolatedDirection(const float rL, CartesianVector &direction) const
496 {
497  const StatusCode statusCode(this->GetGlobalFitDirection(rL, direction));
498 
499  if (STATUS_CODE_NOT_FOUND != statusCode)
500  return statusCode;
501 
502  const int thisLayer(this->GetLayer(rL));
503  const int minLayer(this->GetMinLayer());
504  const int maxLayer(this->GetMaxLayer());
505 
506  if (thisLayer <= minLayer)
507  {
508  direction = this->GetGlobalMinLayerDirection();
509  }
510  else if (thisLayer >= maxLayer)
511  {
512  direction = this->GetGlobalMaxLayerDirection();
513  }
514  else
515  {
516  return STATUS_CODE_FAILURE;
517  }
518 
519  return STATUS_CODE_SUCCESS;
520 }
521 
522 //------------------------------------------------------------------------------------------------------------------------------------------
523 
524 StatusCode TwoDSlidingFitResult::GetExtrapolatedPositionAtX(const float x, CartesianVector &position) const
525 {
526  const StatusCode statusCode(this->GetGlobalFitPositionAtX(x, position));
527 
528  if (STATUS_CODE_NOT_FOUND != statusCode)
529  return statusCode;
530 
531  const float minLayerX(this->GetGlobalMinLayerPosition().GetX());
532  const float maxLayerX(this->GetGlobalMaxLayerPosition().GetX());
533 
534  const int minLayer(this->GetMinLayer());
535  const int maxLayer(this->GetMaxLayer());
536 
537  const int innerLayer((minLayerX < maxLayerX) ? minLayer : maxLayer);
538  const int outerLayer((minLayerX < maxLayerX) ? maxLayer : minLayer);
539 
540  const CartesianVector innerVertex((innerLayer == minLayer) ? this->GetGlobalMinLayerPosition() : this->GetGlobalMaxLayerPosition());
541  const CartesianVector innerDirection((innerLayer == minLayer) ? this->GetGlobalMinLayerDirection() * -1.f : this->GetGlobalMaxLayerDirection());
542 
543  const CartesianVector outerVertex((outerLayer == minLayer) ? this->GetGlobalMinLayerPosition() : this->GetGlobalMaxLayerPosition());
544  const CartesianVector outerDirection((outerLayer == minLayer) ? this->GetGlobalMinLayerDirection() * -1.f : this->GetGlobalMaxLayerDirection());
545 
546  if (innerDirection.GetX() > -std::numeric_limits<float>::epsilon() || outerDirection.GetX() < +std::numeric_limits<float>::epsilon() ||
547  outerVertex.GetX() - innerVertex.GetX() < +std::numeric_limits<float>::epsilon())
548  {
549  return STATUS_CODE_NOT_FOUND;
550  }
551  else if (x >= outerVertex.GetX())
552  {
553  position = outerVertex + outerDirection * ((x - outerVertex.GetX()) / outerDirection.GetX());
554  }
555  else if (x <= innerVertex.GetX())
556  {
557  position = innerVertex + innerDirection * ((x - innerVertex.GetX()) / innerDirection.GetX());
558  }
559  else
560  {
561  return STATUS_CODE_NOT_FOUND;
562  }
563 
564  // TODO How to assign an uncertainty on the extrapolated position?
565  return STATUS_CODE_SUCCESS;
566 }
567 
568 //------------------------------------------------------------------------------------------------------------------------------------------
569 
571 {
572  int layer(this->GetLayer(rL));
573 
574  for (FitSegmentList::const_iterator iter = m_fitSegmentList.begin(), iterEnd = m_fitSegmentList.end(); iter != iterEnd; ++iter)
575  {
576  const FitSegment &fitSegment = *iter;
577 
578  if (layer >= fitSegment.GetStartLayer() && layer <= fitSegment.GetEndLayer())
579  return fitSegment;
580  }
581 
582  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
583 }
584 
585 // Private member functions start here
586 //------------------------------------------------------------------------------------------------------------------------------------------
587 
588 void TwoDSlidingFitResult::CalculateAxes(const CartesianPointVector &coordinateVector, const float layerPitch)
589 {
590  if (coordinateVector.size() < 2)
591  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
592 
593  CartesianVector centroid(0.f, 0.f, 0.f);
594  LArPcaHelper::EigenVectors eigenVecs;
595  LArPcaHelper::EigenValues eigenValues(0.f, 0.f, 0.f);
596  LArPcaHelper::RunPca(coordinateVector, centroid, eigenValues, eigenVecs);
597 
598  float minProjection(std::numeric_limits<float>::max());
599  CartesianVector fitDirection(eigenVecs.at(0));
600 
601  // ATTN TwoDSlidingFitResult convention is to point in direction of increasing z
602  if (fitDirection.GetZ() < 0.f)
603  fitDirection *= -1.f;
604 
605  for (const CartesianVector &coordinate : coordinateVector)
606  minProjection = std::min(minProjection, fitDirection.GetDotProduct(coordinate - centroid));
607 
608  // Define layers based on centroid rather than individual extremal hits
609  const float fitProjection(layerPitch * std::floor(minProjection / layerPitch));
610 
611  m_axisIntercept = centroid + (fitDirection * fitProjection);
612  m_axisDirection = fitDirection;
613 
614  // Use y-axis to generate an orthogonal axis (assuming that cluster occupies X-Z plane)
615  CartesianVector yAxis(0.f, 1.f, 0.f);
616  m_orthoDirection = yAxis.GetCrossProduct(m_axisDirection).GetUnitVector();
617 }
618 
619 //------------------------------------------------------------------------------------------------------------------------------------------
620 
621 void TwoDSlidingFitResult::FillLayerFitContributionMap(const CartesianPointVector &coordinateVector)
622 {
623  if (m_layerPitch < std::numeric_limits<float>::epsilon())
624  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
625 
626  if ((m_axisDirection.GetMagnitudeSquared() < std::numeric_limits<float>::epsilon()) ||
627  (m_orthoDirection.GetMagnitudeSquared() < std::numeric_limits<float>::epsilon()))
628  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
629 
630  if (!m_layerFitContributionMap.empty())
631  throw StatusCodeException(STATUS_CODE_FAILURE);
632 
633  for (CartesianPointVector::const_iterator iter = coordinateVector.begin(), iterEnd = coordinateVector.end(); iter != iterEnd; ++iter)
634  {
635  float rL(0.f), rT(0.f);
636  this->GetLocalPosition(*iter, rL, rT);
637  m_layerFitContributionMap[this->GetLayer(rL)].AddPoint(rL, rT);
638  }
639 }
640 
641 //------------------------------------------------------------------------------------------------------------------------------------------
642 
644 {
645  if (!m_layerFitResultMap.empty())
646  throw StatusCodeException(STATUS_CODE_FAILURE);
647 
648  if ((m_layerPitch < std::numeric_limits<float>::epsilon()) || (m_layerFitContributionMap.empty()))
649  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
650 
651  unsigned int slidingNPoints(0);
652  double slidingSumT(0.), slidingSumL(0.), slidingSumTT(0.), slidingSumLT(0.), slidingSumLL(0.);
653 
654  const LayerFitContributionMap &layerFitContributionMap(this->GetLayerFitContributionMap());
655  const int innerLayer(layerFitContributionMap.begin()->first);
656  const int layerFitHalfWindow(static_cast<int>(this->GetLayerFitHalfWindow()));
657 
658  for (int iLayer = innerLayer; iLayer < innerLayer + layerFitHalfWindow; ++iLayer)
659  {
660  LayerFitContributionMap::const_iterator lyrIter = layerFitContributionMap.find(iLayer);
661 
662  if (layerFitContributionMap.end() != lyrIter)
663  {
664  slidingSumT += lyrIter->second.GetSumT();
665  slidingSumL += lyrIter->second.GetSumL();
666  slidingSumTT += lyrIter->second.GetSumTT();
667  slidingSumLT += lyrIter->second.GetSumLT();
668  slidingSumLL += lyrIter->second.GetSumLL();
669  slidingNPoints += lyrIter->second.GetNPoints();
670  }
671  }
672 
673  const int outerLayer(layerFitContributionMap.rbegin()->first);
674 
675  for (int iLayer = innerLayer; iLayer <= outerLayer; ++iLayer)
676  {
677  const int fwdLayer(iLayer + layerFitHalfWindow);
678  LayerFitContributionMap::const_iterator fwdIter = layerFitContributionMap.find(fwdLayer);
679 
680  if (layerFitContributionMap.end() != fwdIter)
681  {
682  slidingSumT += fwdIter->second.GetSumT();
683  slidingSumL += fwdIter->second.GetSumL();
684  slidingSumTT += fwdIter->second.GetSumTT();
685  slidingSumLT += fwdIter->second.GetSumLT();
686  slidingSumLL += fwdIter->second.GetSumLL();
687  slidingNPoints += fwdIter->second.GetNPoints();
688  }
689 
690  const int bwdLayer(iLayer - layerFitHalfWindow - 1);
691  LayerFitContributionMap::const_iterator bwdIter = layerFitContributionMap.find(bwdLayer);
692 
693  if (layerFitContributionMap.end() != bwdIter)
694  {
695  slidingSumT -= bwdIter->second.GetSumT();
696  slidingSumL -= bwdIter->second.GetSumL();
697  slidingSumTT -= bwdIter->second.GetSumTT();
698  slidingSumLT -= bwdIter->second.GetSumLT();
699  slidingSumLL -= bwdIter->second.GetSumLL();
700  slidingNPoints -= bwdIter->second.GetNPoints();
701  }
702 
703  // require three points for meaningful results
704  if (slidingNPoints <= 2)
705  continue;
706 
707  // only fill the result map if there is an entry in the contribution map
708  if (layerFitContributionMap.end() == layerFitContributionMap.find(iLayer))
709  continue;
710 
711  const double denominator(slidingSumLL - slidingSumL * slidingSumL / static_cast<double>(slidingNPoints));
712 
713  if (std::fabs(denominator) < std::numeric_limits<float>::epsilon())
714  continue;
715 
716  const double gradient((slidingSumLT - slidingSumL * slidingSumT / static_cast<double>(slidingNPoints)) / denominator);
717  const double intercept(
718  (slidingSumLL * slidingSumT / static_cast<double>(slidingNPoints) - slidingSumL * slidingSumLT / static_cast<double>(slidingNPoints)) / denominator);
719  double variance((slidingSumTT - 2. * intercept * slidingSumT - 2. * gradient * slidingSumLT +
720  intercept * intercept * static_cast<double>(slidingNPoints) + 2. * gradient * intercept * slidingSumL +
721  gradient * gradient * slidingSumLL) /
722  (1. + gradient * gradient));
723 
724  if (variance < -std::numeric_limits<float>::epsilon())
725  continue;
726 
727  if (variance < std::numeric_limits<float>::epsilon())
728  variance = 0.;
729 
730  const double rms(std::sqrt(variance / static_cast<double>(slidingNPoints)));
731  const double l(this->GetL(iLayer));
732  const double fitT(intercept + gradient * l);
733 
734  const LayerFitResult layerFitResult(l, fitT, gradient, rms);
735  (void)m_layerFitResultMap.insert(LayerFitResultMap::value_type(iLayer, layerFitResult));
736  }
737 
738  if (m_layerFitResultMap.empty())
739  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
740 }
741 
742 //------------------------------------------------------------------------------------------------------------------------------------------
743 
745 {
746  if (!m_fitSegmentList.empty())
747  throw StatusCodeException(STATUS_CODE_FAILURE);
748 
749  unsigned int nSustainedSteps(0);
750  float sustainedDirectionStartX(0.f), sustainedDirectionEndX(0.f);
751 
752  CartesianVector previousPosition(0.f, 0.f, 0.f);
753  TransverseDirection previousDirection(UNKNOWN), sustainedDirection(UNKNOWN);
754 
755  LayerFitResultMap::const_iterator sustainedDirectionStartIter, sustainedDirectionEndIter;
756  const LayerFitResultMap &layerFitResultMap(this->GetLayerFitResultMap());
757 
758  for (LayerFitResultMap::const_iterator iter = layerFitResultMap.begin(), iterEnd = layerFitResultMap.end(); iter != iterEnd; ++iter)
759  {
760  CartesianVector position(0.f, 0.f, 0.f);
761  this->GetGlobalPosition(iter->second.GetL(), iter->second.GetFitT(), position);
762 
763  // TODO currentDirection could also be UNCHANGED_IN_X
764  const TransverseDirection currentDirection(((position - previousPosition).GetX() > 0.f) ? POSITIVE_IN_X : NEGATIVE_IN_X);
765 
766  if (previousDirection == currentDirection)
767  {
768  ++nSustainedSteps;
769 
770  if (nSustainedSteps > 2)
771  {
772  sustainedDirection = currentDirection;
773  sustainedDirectionEndIter = iter;
774  sustainedDirectionEndX = position.GetX();
775  }
776  }
777  else
778  {
779  if ((POSITIVE_IN_X == sustainedDirection) || (NEGATIVE_IN_X == sustainedDirection))
780  m_fitSegmentList.push_back(FitSegment(
781  sustainedDirectionStartIter->first, sustainedDirectionEndIter->first, sustainedDirectionStartX, sustainedDirectionEndX));
782 
783  nSustainedSteps = 0;
784  sustainedDirection = UNKNOWN;
785  sustainedDirectionStartIter = iter;
786  sustainedDirectionStartX = position.GetX();
787  }
788 
789  previousPosition = position;
790  previousDirection = currentDirection;
791  }
792 
793  if ((POSITIVE_IN_X == sustainedDirection) || (NEGATIVE_IN_X == sustainedDirection))
794  m_fitSegmentList.push_back(
795  FitSegment(sustainedDirectionStartIter->first, sustainedDirectionEndIter->first, sustainedDirectionStartX, sustainedDirectionEndX));
796 }
797 
798 //------------------------------------------------------------------------------------------------------------------------------------------
799 
800 void TwoDSlidingFitResult::GetMinAndMaxCoordinate(const bool isX, float &min, float &max) const
801 {
802  if (m_layerFitResultMap.empty())
803  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
804 
805  min = std::numeric_limits<float>::max();
806  max = -std::numeric_limits<float>::max();
807 
808  for (LayerFitResultMap::const_iterator iter = m_layerFitResultMap.begin(), iterEnd = m_layerFitResultMap.end(); iter != iterEnd; ++iter)
809  {
810  CartesianVector globalPosition(0.f, 0.f, 0.f);
811  this->GetGlobalPosition(iter->second.GetL(), iter->second.GetFitT(), globalPosition);
812  const float coordinate(isX ? globalPosition.GetX() : globalPosition.GetZ());
813  min = std::min(min, coordinate);
814  max = std::max(max, coordinate);
815  }
816 }
817 
818 //------------------------------------------------------------------------------------------------------------------------------------------
819 
820 CartesianVector TwoDSlidingFitResult::GetGlobalFitPosition(const LayerInterpolation &layerInterpolation) const
821 {
822  const LayerFitResultMap::const_iterator firstLayerIter(layerInterpolation.GetStartLayerIter());
823  const LayerFitResultMap::const_iterator secondLayerIter(layerInterpolation.GetEndLayerIter());
824 
825  const float firstWeight(layerInterpolation.GetStartLayerWeight());
826  const float secondWeight(layerInterpolation.GetEndLayerWeight());
827 
828  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
829  throw StatusCodeException(STATUS_CODE_FAILURE);
830 
831  CartesianVector firstLayerPosition(0.f, 0.f, 0.f);
832  this->GetGlobalPosition(firstLayerIter->second.GetL(), firstLayerIter->second.GetFitT(), firstLayerPosition);
833 
834  if (firstLayerIter == secondLayerIter)
835  return firstLayerPosition;
836 
837  CartesianVector secondLayerPosition(0.f, 0.f, 0.f);
838  this->GetGlobalPosition(secondLayerIter->second.GetL(), secondLayerIter->second.GetFitT(), secondLayerPosition);
839 
840  if (firstWeight + secondWeight < std::numeric_limits<float>::epsilon())
841  throw StatusCodeException(STATUS_CODE_FAILURE);
842 
843  return ((firstLayerPosition * firstWeight + secondLayerPosition * secondWeight) * (1.f / (firstWeight + secondWeight)));
844 }
845 
846 //------------------------------------------------------------------------------------------------------------------------------------------
847 
848 CartesianVector TwoDSlidingFitResult::GetGlobalFitDirection(const LayerInterpolation &layerInterpolation) const
849 {
850  const LayerFitResultMap::const_iterator firstLayerIter(layerInterpolation.GetStartLayerIter());
851  const LayerFitResultMap::const_iterator secondLayerIter(layerInterpolation.GetEndLayerIter());
852 
853  const float firstWeight(layerInterpolation.GetStartLayerWeight());
854  const float secondWeight(layerInterpolation.GetEndLayerWeight());
855 
856  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
857  throw StatusCodeException(STATUS_CODE_FAILURE);
858 
859  CartesianVector firstLayerDirection(0.f, 0.f, 0.f);
860  this->GetGlobalDirection(firstLayerIter->second.GetGradient(), firstLayerDirection);
861 
862  if (firstLayerIter == secondLayerIter)
863  return firstLayerDirection;
864 
865  CartesianVector secondLayerDirection(0.f, 0.f, 0.f);
866  this->GetGlobalDirection(secondLayerIter->second.GetGradient(), secondLayerDirection);
867 
868  if (firstWeight + secondWeight < std::numeric_limits<float>::epsilon())
869  throw StatusCodeException(STATUS_CODE_FAILURE);
870 
871  return ((firstLayerDirection * firstWeight + secondLayerDirection * secondWeight).GetUnitVector());
872 }
873 
874 //------------------------------------------------------------------------------------------------------------------------------------------
875 
876 float TwoDSlidingFitResult::GetFitRms(const LayerInterpolation &layerInterpolation) const
877 {
878  const LayerFitResultMap::const_iterator firstLayerIter(layerInterpolation.GetStartLayerIter());
879  const LayerFitResultMap::const_iterator secondLayerIter(layerInterpolation.GetEndLayerIter());
880 
881  const float firstWeight(layerInterpolation.GetStartLayerWeight());
882  const float secondWeight(layerInterpolation.GetEndLayerWeight());
883 
884  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
885  throw StatusCodeException(STATUS_CODE_FAILURE);
886 
887  const float firstLayerRms(firstLayerIter->second.GetRms());
888 
889  if (firstLayerIter == secondLayerIter)
890  return firstLayerRms;
891 
892  const float secondLayerRms(secondLayerIter->second.GetRms());
893 
894  if (firstWeight + secondWeight < std::numeric_limits<float>::epsilon())
895  throw StatusCodeException(STATUS_CODE_FAILURE);
896 
897  return ((firstLayerRms * firstWeight + secondLayerRms * secondWeight) / (firstWeight + secondWeight));
898 }
899 
900 //------------------------------------------------------------------------------------------------------------------------------------------
901 
902 StatusCode TwoDSlidingFitResult::LongitudinalInterpolation(const float rL, LayerInterpolation &layerInterpolation) const
903 {
904  double firstWeight(0.), secondWeight(0.);
905  LayerFitResultMap::const_iterator firstLayerIter, secondLayerIter;
906 
907  const StatusCode statusCode(this->GetLongitudinalSurroundingLayers(rL, firstLayerIter, secondLayerIter));
908 
909  if (STATUS_CODE_SUCCESS != statusCode)
910  return statusCode;
911 
912  this->GetLongitudinalInterpolationWeights(rL, firstLayerIter, secondLayerIter, firstWeight, secondWeight);
913  layerInterpolation = LayerInterpolation(firstLayerIter, secondLayerIter, firstWeight, secondWeight);
914 
915  return STATUS_CODE_SUCCESS;
916 }
917 
918 //------------------------------------------------------------------------------------------------------------------------------------------
919 
920 StatusCode TwoDSlidingFitResult::TransverseInterpolation(const float x, const FitSegment &fitSegment, LayerInterpolation &layerInterpolation) const
921 {
922  double firstWeight(0.), secondWeight(0.);
923  LayerFitResultMap::const_iterator firstLayerIter, secondLayerIter;
924 
925  const StatusCode statusCode(
926  this->GetTransverseSurroundingLayers(x, fitSegment.GetStartLayer(), fitSegment.GetEndLayer(), firstLayerIter, secondLayerIter));
927 
928  if (STATUS_CODE_SUCCESS != statusCode)
929  return statusCode;
930 
931  this->GetTransverseInterpolationWeights(x, firstLayerIter, secondLayerIter, firstWeight, secondWeight);
932  layerInterpolation = LayerInterpolation(firstLayerIter, secondLayerIter, firstWeight, secondWeight);
933 
934  return STATUS_CODE_SUCCESS;
935 }
936 
937 //------------------------------------------------------------------------------------------------------------------------------------------
938 
939 StatusCode TwoDSlidingFitResult::TransverseInterpolation(const float x, LayerInterpolationList &layerInterpolationList) const
940 {
941  for (FitSegmentList::const_iterator iter = m_fitSegmentList.begin(), iterEnd = m_fitSegmentList.end(); iter != iterEnd; ++iter)
942  {
943  LayerInterpolation layerInterpolation;
944  const StatusCode statusCode(this->TransverseInterpolation(x, *iter, layerInterpolation));
945 
946  if (STATUS_CODE_SUCCESS == statusCode)
947  {
948  layerInterpolationList.push_back(layerInterpolation);
949  }
950  else if (STATUS_CODE_NOT_FOUND != statusCode)
951  {
952  return statusCode;
953  }
954  }
955 
956  return STATUS_CODE_SUCCESS;
957 }
958 
959 //------------------------------------------------------------------------------------------------------------------------------------------
960 
962  const float rL, LayerFitResultMap::const_iterator &firstLayerIter, LayerFitResultMap::const_iterator &secondLayerIter) const
963 {
964  if (m_layerFitResultMap.empty())
965  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
966 
967  // Get minimum, maximum and input layers
968  const int minLayer(m_layerFitResultMap.begin()->first), maxLayer(m_layerFitResultMap.rbegin()->first);
969  const int thisLayer(this->GetLayer(rL));
970 
971  // Allow special case of single-layer sliding fit result
972  if (minLayer == thisLayer && thisLayer == maxLayer)
973  {
974  firstLayerIter = m_layerFitResultMap.find(minLayer);
975  secondLayerIter = m_layerFitResultMap.find(maxLayer);
976  return STATUS_CODE_SUCCESS;
977  }
978 
979  // For multi-layer sliding fit result, set start layer and bail out if out of range
980  const int startLayer((thisLayer >= maxLayer) ? thisLayer - 1 : thisLayer);
981 
982  if ((startLayer < minLayer) || (startLayer >= maxLayer))
983  return STATUS_CODE_NOT_FOUND;
984 
985  // First layer iterator
986  firstLayerIter = m_layerFitResultMap.end();
987 
988  for (int iLayer = startLayer; iLayer >= minLayer; --iLayer)
989  {
990  firstLayerIter = m_layerFitResultMap.find(iLayer);
991 
992  if (m_layerFitResultMap.end() != firstLayerIter)
993  break;
994  }
995 
996  if (m_layerFitResultMap.end() == firstLayerIter)
997  return STATUS_CODE_NOT_FOUND;
998 
999  // Second layer iterator
1000  secondLayerIter = m_layerFitResultMap.end();
1001 
1002  for (int iLayer = startLayer + 1; iLayer <= maxLayer; ++iLayer)
1003  {
1004  secondLayerIter = m_layerFitResultMap.find(iLayer);
1005 
1006  if (m_layerFitResultMap.end() != secondLayerIter)
1007  break;
1008  }
1009 
1010  if (m_layerFitResultMap.end() == secondLayerIter)
1011  return STATUS_CODE_NOT_FOUND;
1012 
1013  return STATUS_CODE_SUCCESS;
1014 }
1015 
1016 //------------------------------------------------------------------------------------------------------------------------------------------
1017 
1018 StatusCode TwoDSlidingFitResult::GetTransverseSurroundingLayers(const float x, const int minLayer, const int maxLayer,
1019  LayerFitResultMap::const_iterator &firstLayerIter, LayerFitResultMap::const_iterator &secondLayerIter) const
1020 {
1021  if (m_layerFitResultMap.empty())
1022  throw StatusCodeException(STATUS_CODE_NOT_INITIALIZED);
1023 
1024  LayerFitResultMap::const_iterator minLayerIter = m_layerFitResultMap.find(minLayer);
1025  if (m_layerFitResultMap.end() == minLayerIter)
1026  throw StatusCodeException(STATUS_CODE_FAILURE);
1027 
1028  LayerFitResultMap::const_iterator maxLayerIter = m_layerFitResultMap.find(maxLayer);
1029  if (m_layerFitResultMap.end() == maxLayerIter)
1030  throw StatusCodeException(STATUS_CODE_FAILURE);
1031 
1032  CartesianVector minPosition(0.f, 0.f, 0.f), maxPosition(0.f, 0.f, 0.f);
1033  this->GetGlobalPosition(minLayerIter->second.GetL(), minLayerIter->second.GetFitT(), minPosition);
1034  this->GetGlobalPosition(maxLayerIter->second.GetL(), maxLayerIter->second.GetFitT(), maxPosition);
1035 
1036  if ((std::fabs(maxPosition.GetX() - minPosition.GetX()) < std::numeric_limits<float>::epsilon()))
1037  return STATUS_CODE_NOT_FOUND;
1038 
1039  // Find start layer
1040  const float minL(minLayerIter->second.GetL());
1041  const float maxL(maxLayerIter->second.GetL());
1042  const float startL(minL + (maxL - minL) * (x - minPosition.GetX()) / (maxPosition.GetX() - minPosition.GetX()));
1043  const int startLayer(std::max(minLayer, std::min(maxLayer, this->GetLayer(startL))));
1044 
1045  // Find nearest layer iterator to start layer
1047  CartesianVector startLayerPosition(0.f, 0.f, 0.f);
1048 
1049  for (int iLayer = startLayer; iLayer <= maxLayer; ++iLayer)
1050  {
1051  startLayerIter = m_layerFitResultMap.find(iLayer);
1052 
1053  if (m_layerFitResultMap.end() != startLayerIter)
1054  break;
1055  }
1056 
1057  if (m_layerFitResultMap.end() == startLayerIter)
1058  return STATUS_CODE_NOT_FOUND;
1059 
1060  this->GetGlobalPosition(startLayerIter->second.GetL(), startLayerIter->second.GetFitT(), startLayerPosition);
1061 
1062  const bool startIsAhead((startLayerPosition.GetX() - x) > std::numeric_limits<float>::epsilon());
1063  const bool increasesWithLayers(maxPosition.GetX() > minPosition.GetX());
1064  const int increment = ((startIsAhead == increasesWithLayers) ? -1 : +1);
1065 
1066  // Find surrounding layer iterators
1067  // (Second layer iterator comes immediately after the fit has crossed the target X coordinate
1068  // and first layer iterator comes immediately before the second layer iterator).
1069  firstLayerIter = m_layerFitResultMap.end();
1070  secondLayerIter = m_layerFitResultMap.end();
1071 
1072  CartesianVector firstLayerPosition(0.f, 0.f, 0.f);
1073  CartesianVector secondLayerPosition(0.f, 0.f, 0.f);
1074 
1075  for (int iLayer = startLayerIter->first; (iLayer >= minLayer) && (iLayer <= maxLayer); iLayer += increment)
1076  {
1078  if (m_layerFitResultMap.end() == tempIter)
1079  continue;
1080 
1081  firstLayerIter = secondLayerIter;
1082  firstLayerPosition = secondLayerPosition;
1083  secondLayerIter = tempIter;
1084 
1085  this->GetGlobalPosition(secondLayerIter->second.GetL(), secondLayerIter->second.GetFitT(), secondLayerPosition);
1086  const bool isAhead(secondLayerPosition.GetX() > x);
1087 
1088  if (startIsAhead != isAhead)
1089  break;
1090 
1091  firstLayerIter = m_layerFitResultMap.end();
1092  }
1093 
1094  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
1095  return STATUS_CODE_NOT_FOUND;
1096 
1097  return STATUS_CODE_SUCCESS;
1098 }
1099 
1100 //------------------------------------------------------------------------------------------------------------------------------------------
1101 
1103  const LayerFitResultMap::const_iterator &secondLayerIter, double &firstWeight, double &secondWeight) const
1104 {
1105  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
1106  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
1107 
1108  const double deltaL(rL - firstLayerIter->second.GetL());
1109  const double deltaLLayers(secondLayerIter->second.GetL() - firstLayerIter->second.GetL());
1110 
1111  if (std::fabs(deltaLLayers) > std::numeric_limits<float>::epsilon())
1112  {
1113  firstWeight = 1. - deltaL / deltaLLayers;
1114  secondWeight = deltaL / deltaLLayers;
1115  }
1116  else
1117  {
1118  firstWeight = 0.5;
1119  secondWeight = 0.5;
1120  }
1121 }
1122 
1123 //------------------------------------------------------------------------------------------------------------------------------------------
1124 
1126  const LayerFitResultMap::const_iterator &secondLayerIter, double &firstWeight, double &secondWeight) const
1127 {
1128  if (m_layerFitResultMap.end() == firstLayerIter || m_layerFitResultMap.end() == secondLayerIter)
1129  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
1130 
1131  CartesianVector firstLayerPosition(0.f, 0.f, 0.f);
1132  CartesianVector secondLayerPosition(0.f, 0.f, 0.f);
1133 
1134  this->GetGlobalPosition(firstLayerIter->second.GetL(), firstLayerIter->second.GetFitT(), firstLayerPosition);
1135  this->GetGlobalPosition(secondLayerIter->second.GetL(), secondLayerIter->second.GetFitT(), secondLayerPosition);
1136 
1137  const double deltaP(x - firstLayerPosition.GetX());
1138  const double deltaPLayers(secondLayerPosition.GetX() - firstLayerPosition.GetX());
1139 
1140  if (std::fabs(deltaPLayers) > std::numeric_limits<float>::epsilon())
1141  {
1142  firstWeight = 1. - deltaP / deltaPLayers;
1143  secondWeight = deltaP / deltaPLayers;
1144  }
1145  else
1146  {
1147  firstWeight = 0.5;
1148  secondWeight = 0.5;
1149  }
1150 }
1151 
1152 } // namespace lar_content
Float_t x
Definition: compare.C:6
void FillLayerFitContributionMap(const pandora::CartesianPointVector &coordinateVector)
Fill the layer fit contribution map.
unsigned int GetLayerFitHalfWindow() const
Get the layer fit half window.
void GetMinAndMaxCoordinate(const bool isX, float &min, float &max) const
Get the minimum and maximum x or z coordinates associated with the sliding fit.
TwoDSlidingFitResult(const T *const pT, const unsigned int layerFitHalfWindow, const float layerPitch, const float axisDeviationLimitForHitDivision=0.95f)
Constructor using internal definition of primary axis.
pandora::CartesianVector m_axisIntercept
The axis intercept position.
int GetStartLayer() const
Get start layer.
void GetTransverseInterpolationWeights(const float x, const LayerFitResultMap::const_iterator &firstLayerIter, const LayerFitResultMap::const_iterator &secondLayerIter, double &firstWeight, double &secondWeight) const
Get interpolation weights for layers surrounding a specified transverse position. ...
pandora::StatusCode GetTransverseSurroundingLayers(const float x, const int minLayer, const int maxLayer, LayerFitResultMap::const_iterator &firstLayerIter, LayerFitResultMap::const_iterator &secondLayerIter) const
Get iterators for layers surrounding a specified transverse position.
Header file for the lar hit width helper class.
const pandora::Cluster * m_pCluster
The address of the cluster.
double GetEndLayerWeight() const
Get the end layer weight.
TransverseDirection
TransverseDirection enum.
void FindSlidingFitSegments()
Find sliding fit segments; sections with tramsverse direction.
pandora::StatusCode GetExtrapolatedPositionAtX(const float x, pandora::CartesianVector &position) const
Get extrapolated position (beyond span) for a given input x coordinate.
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.
pandora::CartesianVector m_axisDirection
The axis direction vector.
void GetLongitudinalInterpolationWeights(const float rL, const LayerFitResultMap::const_iterator &firstLayerIter, const LayerFitResultMap::const_iterator &secondLayerIter, double &firstWeight, double &secondWeight) const
Get interpolation weights for layers surrounding a specified longitudinal position.
static ConstituentHitVector GetConstituentHits(const pandora::Cluster *const pCluster, const float maxConstituentHitWidth, const float hitWidthScalingFactor, const bool isUniform)
Break up the cluster hits into constituent hits.
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
Header file for the principal curve analysis helper class.
intermediate_table::const_iterator const_iterator
pandora::StatusCode LongitudinalInterpolation(const float rL, LayerInterpolation &layerInterpolation) const
Get the pair of layers surrounding a specified longitudinal position.
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
pandora::StatusCode GetLongitudinalSurroundingLayers(const float rL, LayerFitResultMap::const_iterator &firstLayerIter, LayerFitResultMap::const_iterator &secondLayerIter) const
Get iterators for layers surrounding the specified longitudinal position.
std::map< int, LayerFitContribution > LayerFitContributionMap
float GetCosScatteringAngle(const float rL) const
Get scattering angle for a given longitudinal coordinate.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
unsigned int m_layerFitHalfWindow
The layer fit half window.
int GetMaxLayer() const
Get the maximum occupied layer in the sliding fit.
LayerFitResultMap::const_iterator GetEndLayerIter() const
Get the end layer iterator.
int GetMinLayer() const
Get the minimum occupied layer in the sliding fit.
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.
double GetStartLayerWeight() const
Get the start layer weight.
std::vector< ConstituentHit > ConstituentHitVector
Header file for the lar two dimensional sliding fit result class.
pandora::StatusCode TransverseInterpolation(const float x, const FitSegment &fitSegment, LayerInterpolation &layerInterpolation) const
Get the surrounding pair of layers for a specified transverse position and fit segment.
pandora::StatusCode GetExtrapolatedPosition(const float rL, pandora::CartesianVector &position) const
Get extrapolated position (beyond span) for a given input coordinate.
std::map< int, LayerFitResult > LayerFitResultMap
static pandora::CartesianPointVector GetConstituentHitPositionVector(const ConstituentHitVector &constituentHitVector)
Obtain a vector of the contituent hit central positions.
void CalculateAxes(const pandora::CartesianPointVector &coordinateVector, const float layerPitch)
Calculate the longitudinal and transverse axes.
const FitSegment & GetFitSegment(const float rL) const
Get fit segment for a given longitudinal coordinate.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
pandora::StatusCode GetExtrapolatedDirection(const float rL, pandora::CartesianVector &direction) const
Get extrapolated direction (beyond span) for a given input coordinate.
std::vector< LayerInterpolation > LayerInterpolationList
void GetGlobalDirection(const float dTdL, pandora::CartesianVector &direction) const
Get global direction coordinates for given sliding linear fit gradient.
const LayerFitResultMap & GetLayerFitResultMap() const
Get the layer fit result map.
float GetMinLayerRms() const
Get rms at minimum layer.
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.
pandora::StatusCode GetGlobalFitPositionAtX(const float x, pandora::CartesianVector &position) const
Get global fit position for a given input x coordinate.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
std::vector< pandora::CartesianVector > EigenVectors
Definition: LArPcaHelper.h:25
pandora::StatusCode GetGlobalFitProjection(const pandora::CartesianVector &inputPosition, pandora::CartesianVector &projectedPosition) const
Get projected position on global fit for a given position vector.
void PerformSlidingLinearFit()
Perform the sliding linear fit.
void GetGlobalPosition(const float rL, const float rT, pandora::CartesianVector &position) const
Get global coordinates for given sliding linear fit coordinates.
FitSegmentList m_fitSegmentList
The fit segment list.
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
float GetL(const int layer) const
Get longitudinal coordinate for a given sliding linear fit layer number.
LayerFitContributionMap m_layerFitContributionMap
The layer fit contribution map.
pandora::StatusCode GetGlobalFitDirectionAtX(const float x, pandora::CartesianVector &direction) const
Get global fit direction for a given input x coordinate.
pandora::StatusCode GetTransverseProjection(const float x, const FitSegment &fitSegment, pandora::CartesianVector &position) const
Get projected position for a given input x coordinate and fit segment.
static void GetCoordinateVector(const pandora::Cluster *const pCluster, pandora::CartesianPointVector &coordinateVector)
Get vector of hit coordinates from an input cluster.
LayerFitResultMap::const_iterator GetStartLayerIter() const
Get the start layer iterator.
float m_layerPitch
The layer pitch, units cm.
pandora::StatusCode GetGlobalFitPositionListAtX(const float x, pandora::CartesianPointVector &positionList) const
Get a list of projected positions for a given input x coordinate.
const LayerFitContributionMap & GetLayerFitContributionMap() const
Get the layer fit contribution map.
LayerFitResultMap m_layerFitResultMap
The layer fit result map.
int GetEndLayer() const
Get end layer.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
pandora::CartesianVector m_orthoDirection
The orthogonal direction vector.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
int GetLayer(const float rL) const
Get layer number for given sliding linear fit longitudinal coordinate.