LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArThreeDSlidingConeFitResult.cc
Go to the documentation of this file.
1 
9 #include "Objects/Cluster.h"
10 
12 
14 
15 #include <iterator>
16 
17 using namespace pandora;
18 
19 namespace lar_content
20 {
21 
22 float SimpleCone::GetMeanRT(const Cluster *const pCluster) const
23 {
24  CartesianPointVector hitPositionVector;
25  LArClusterHelper::GetCoordinateVector(pCluster, hitPositionVector);
26 
27  float rTSum(0.f);
28  const unsigned int nClusterHits(pCluster->GetNCaloHits());
29 
30  for (const CartesianVector &hitPosition : hitPositionVector)
31  {
32  const CartesianVector displacement(hitPosition - this->GetConeApex());
33  const float rT(displacement.GetCrossProduct(this->GetConeDirection()).GetMagnitude());
34  rTSum += rT;
35  }
36 
37  return ((nClusterHits > 0) ? rTSum / static_cast<float>(nClusterHits) : 0.f);
38 }
39 
40 //------------------------------------------------------------------------------------------------------------------------------------------
41 
42 float SimpleCone::GetBoundedHitFraction(const Cluster *const pCluster, const float coneLength, const float coneTanHalfAngle) const
43 {
44  CartesianPointVector hitPositionVector;
45  LArClusterHelper::GetCoordinateVector(pCluster, hitPositionVector);
46 
47  unsigned int nMatchedHits(0);
48  const unsigned int nClusterHits(pCluster->GetNCaloHits());
49 
50  for (const CartesianVector &hitPosition : hitPositionVector)
51  {
52  const CartesianVector displacement(hitPosition - this->GetConeApex());
53  const float rL(displacement.GetDotProduct(this->GetConeDirection()));
54 
55  if ((rL < 0.f) || (rL > coneLength))
56  continue;
57 
58  const float rT(displacement.GetCrossProduct(this->GetConeDirection()).GetMagnitude());
59 
60  if (rL * coneTanHalfAngle > rT)
61  ++nMatchedHits;
62  }
63 
64  return ((nClusterHits > 0) ? static_cast<float>(nMatchedHits) / static_cast<float>(nClusterHits) : 0.f);
65 }
66 
67 //------------------------------------------------------------------------------------------------------------------------------------------
68 //------------------------------------------------------------------------------------------------------------------------------------------
69 
70 template <typename T>
71 ThreeDSlidingConeFitResult::ThreeDSlidingConeFitResult(const T *const pT, const unsigned int slidingFitWindow, const float slidingFitLayerPitch) :
72  m_slidingFitResult(ThreeDSlidingFitResult(pT, slidingFitWindow, slidingFitLayerPitch))
73 {
74  const CartesianVector &minLayerPosition3D(m_slidingFitResult.GetGlobalMinLayerPosition());
75  const CartesianVector &maxLayerPosition3D(m_slidingFitResult.GetGlobalMaxLayerPosition());
76 
79 
80  const LayerFitContributionMap &contributionMap1(fitResult1.GetLayerFitContributionMap());
81  const LayerFitContributionMap &contributionMap2(fitResult2.GetLayerFitContributionMap());
82 
83  const int nSteps(static_cast<int>((maxLayerPosition3D - minLayerPosition3D).GetMagnitude() / slidingFitLayerPitch));
84 
85  for (int iStep = 0; iStep <= nSteps; ++iStep)
86  {
87  try
88  {
89  const float rL((static_cast<float>(iStep) + 0.5f) * slidingFitLayerPitch);
90  CartesianVector fitPosition3D(0.f, 0.f, 0.f), fitDirection3D(0.f, 0.f, 0.f);
91 
92  if ((STATUS_CODE_SUCCESS != m_slidingFitResult.GetGlobalFitPosition(rL, fitPosition3D)) ||
93  (STATUS_CODE_SUCCESS != m_slidingFitResult.GetGlobalFitDirection(rL, fitDirection3D)))
94  {
95  continue;
96  }
97 
98  // ATTN Do not include layers without contributions in track state map (contain only interpolations)
99  if (!contributionMap1.count(fitResult1.GetLayer(rL)) && !contributionMap2.count(fitResult2.GetLayer(rL)))
100  continue;
101 
102  (void)m_trackStateMap.insert(TrackStateMap::value_type(iStep, TrackState(fitPosition3D, fitDirection3D)));
103  }
104  catch (const StatusCodeException &)
105  {
106  /* Deliberately empty */
107  }
108  }
109 }
110 
111 //------------------------------------------------------------------------------------------------------------------------------------------
112 
114  const unsigned int nLayersForConeFit, const unsigned int nCones, const ConeSelection coneSelection, SimpleConeList &simpleConeList) const
115 {
116  const TrackStateMap &trackStateMap(this->GetTrackStateMap());
117  const unsigned int nLayers(trackStateMap.size());
118 
119  if (nLayers + 1 < nLayersForConeFit + nCones)
120  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
121 
122  // Calculate intervals and offsets such that cones are evenly spaced along sliding fit and are equidistant from the extremal layers
123  const unsigned int coneInterval((nCones > 1) ? (nLayers - nLayersForConeFit) / (nCones - 1) : 1);
124  const ThreeDSlidingFitResult &slidingFitResult(this->GetSlidingFitResult());
125  const CartesianVector coneDisplacement(slidingFitResult.GetGlobalMaxLayerPosition() - slidingFitResult.GetGlobalMinLayerPosition());
126  const bool isForward(coneDisplacement.GetZ() > std::numeric_limits<float>::epsilon());
127  const unsigned int coneOffset1((nLayers - nLayersForConeFit - (nCones - 1) * coneInterval) / 2);
128  const unsigned int coneOffset2((1 == nLayers % 2) && isForward ? 1 : 0);
129  const unsigned int coneOffset(coneOffset1 + coneOffset2);
130 
131  TrackStateLinkedList trackStateList;
132  CartesianVector directionSum(0.f, 0.f, 0.f);
133  const float clusterLength((trackStateMap.begin()->second.GetPosition() - trackStateMap.rbegin()->second.GetPosition()).GetMagnitude());
134 
135  unsigned int nConeSamplingSteps(0);
136 
137  for (TrackStateMap::const_iterator iter = trackStateMap.begin(), iterEnd = trackStateMap.end(); iter != iterEnd; ++iter)
138  {
139  if (nConeSamplingSteps >= nCones)
140  return;
141 
142  trackStateList.push_back(iter->second);
143  directionSum += iter->second.GetMomentum();
144 
145  const unsigned int beginDistance(static_cast<unsigned int>(std::distance(trackStateMap.begin(), iter)));
146 
147  if (beginDistance + 1 < nLayersForConeFit)
148  continue;
149 
150  const TrackState &maxLayerTrackState(trackStateList.back());
151  const TrackState &minLayerTrackState(trackStateList.front());
152 
153  if ((beginDistance + 1 >= nLayersForConeFit + coneOffset) && (beginDistance + 1 - nLayersForConeFit - coneOffset) % coneInterval == 0)
154  {
155  const CartesianVector &minLayerApex(minLayerTrackState.GetPosition());
156  const CartesianVector &maxLayerApex(maxLayerTrackState.GetPosition());
157 
158  const CartesianVector minLayerDirection(directionSum.GetUnitVector());
159  const CartesianVector maxLayerDirection(directionSum.GetUnitVector() * -1.f);
160 
161  // TODO Estimate cone length and angle here too, maybe by projecting positions onto direction and looking at rT distribution?
162  ++nConeSamplingSteps;
163  const float placeHolderTanHalfAngle(0.5f);
164 
165  if ((CONE_FORWARD_ONLY == coneSelection) || (CONE_BOTH_DIRECTIONS == coneSelection))
166  simpleConeList.push_back(SimpleCone(minLayerApex, minLayerDirection, clusterLength, placeHolderTanHalfAngle));
167 
168  if ((CONE_BACKWARD_ONLY == coneSelection) || (CONE_BOTH_DIRECTIONS == coneSelection))
169  simpleConeList.push_back(SimpleCone(maxLayerApex, maxLayerDirection, clusterLength, placeHolderTanHalfAngle));
170  }
171 
172  directionSum -= minLayerTrackState.GetMomentum();
173  trackStateList.pop_front();
174  }
175 }
176 
177 //------------------------------------------------------------------------------------------------------------------------------------------
178 //------------------------------------------------------------------------------------------------------------------------------------------
179 
180 template ThreeDSlidingConeFitResult::ThreeDSlidingConeFitResult(const pandora::Cluster *const, const unsigned int, const float);
181 template ThreeDSlidingConeFitResult::ThreeDSlidingConeFitResult(const pandora::CartesianPointVector *const, const unsigned int, const float);
182 
183 } // namespace lar_content
const ThreeDSlidingFitResult m_slidingFitResult
The sliding fit result for the full cluster.
ThreeDSlidingConeFitResult(const T *const pT, const unsigned int slidingFitWindow, const float slidingFitLayerPitch)
Constructor.
const TwoDSlidingFitResult & GetSecondFitResult() const
Get the second sliding fit result for this cluster.
std::vector< SimpleCone > SimpleConeList
intermediate_table::const_iterator const_iterator
const ThreeDSlidingFitResult & GetSlidingFitResult() const
Get the sliding fit result for the full cluster.
std::list< pandora::TrackState > TrackStateLinkedList
The track state linked list typedef.
Header file for the lar three dimensional sliding cone fit result class.
void GetSimpleConeList(const unsigned int nLayersForConeFit, const unsigned int nCones, const ConeSelection coneSelection, SimpleConeList &simpleConeList) const
Get the list of simple cones fitted to the three dimensional cluster.
std::map< int, LayerFitContribution > LayerFitContributionMap
TFile f
Definition: plotHisto.C:6
std::map< int, pandora::TrackState > TrackStateMap
Header file for the cluster helper class.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
const pandora::CartesianVector & GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
const TrackStateMap & GetTrackStateMap() const
Get the track state map, which caches results from the sliding fit result.
const TwoDSlidingFitResult & GetFirstFitResult() const
Get the first sliding fit result for this cluster.
ConeSelection
ConeSelection enum.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
TrackStateMap m_trackStateMap
The track state map.
const pandora::CartesianVector & GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.