LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LongitudinalTrackHitsBaseTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
15 
16 using namespace pandora;
17 
18 namespace lar_content
19 {
20 
21 LongitudinalTrackHitsBaseTool::LongitudinalTrackHitsBaseTool() :
22  m_vtxDisplacementCutSquared(5.f * 5.f),
23  m_minTrackLengthSquared(7.5f * 7.5f)
24 {
25 }
26 
27 //------------------------------------------------------------------------------------------------------------------------------------------
28 
30  const CaloHitVector &inputTwoDHits, const MatchedSlidingFitMap &inputSlidingFitMap, ProtoHitVector &protoHitVector) const
31 {
32  MatchedSlidingFitMap matchedSlidingFitMap;
33  CartesianVector vtx3D(0.f, 0.f, 0.f), end3D(0.f, 0.f, 0.f);
34  this->GetVertexAndEndPositions(inputSlidingFitMap, matchedSlidingFitMap, vtx3D, end3D);
35 
36  for (const CaloHit *const pCaloHit2D : inputTwoDHits)
37  {
38  try
39  {
40  ProtoHit protoHit(pCaloHit2D);
41  this->GetLongitudinalTrackHit3D(matchedSlidingFitMap, vtx3D, end3D, protoHit);
42 
43  if (protoHit.IsPositionSet() && (protoHit.GetChi2() < m_chiSquaredCut))
44  protoHitVector.push_back(protoHit);
45  }
46  catch (StatusCodeException &)
47  {
48  }
49  }
50 }
51 
52 //------------------------------------------------------------------------------------------------------------------------------------------
53 
55  MatchedSlidingFitMap &outputSlidingFitMap, CartesianVector &bestVtx3D, CartesianVector &bestEnd3D) const
56 {
57  // TODO Tidy up: The code below is quite repetitive...
58  MatchedSlidingFitMap::const_iterator iterU = inputSlidingFitMap.find(TPC_VIEW_U);
59  const bool foundU(inputSlidingFitMap.end() != iterU);
60 
61  MatchedSlidingFitMap::const_iterator iterV = inputSlidingFitMap.find(TPC_VIEW_V);
62  const bool foundV(inputSlidingFitMap.end() != iterV);
63 
64  MatchedSlidingFitMap::const_iterator iterW = inputSlidingFitMap.find(TPC_VIEW_W);
65  const bool foundW(inputSlidingFitMap.end() != iterW);
66 
67  bool useU(false), useV(false), useW(false);
68  float bestChi2(std::numeric_limits<float>::max());
69 
70  for (unsigned int iPermutation = 0; iPermutation < 4; ++iPermutation)
71  {
72  const bool isForwardU((1 == iPermutation) ? false : true);
73  const bool isForwardV((2 == iPermutation) ? false : true);
74  const bool isForwardW((3 == iPermutation) ? false : true);
75 
76  CartesianVector vtxU(0.f, 0.f, 0.f), endU(0.f, 0.f, 0.f);
77  CartesianVector vtxV(0.f, 0.f, 0.f), endV(0.f, 0.f, 0.f);
78  CartesianVector vtxW(0.f, 0.f, 0.f), endW(0.f, 0.f, 0.f);
79 
80  if (foundU)
81  {
82  const TwoDSlidingFitResult &slidingFitResultU = iterU->second;
83  vtxU = (isForwardU ? slidingFitResultU.GetGlobalMinLayerPosition() : slidingFitResultU.GetGlobalMaxLayerPosition());
84  endU = (isForwardU ? slidingFitResultU.GetGlobalMaxLayerPosition() : slidingFitResultU.GetGlobalMinLayerPosition());
85  }
86 
87  if (foundV)
88  {
89  const TwoDSlidingFitResult &slidingFitResultV = iterV->second;
90  vtxV = (isForwardV ? slidingFitResultV.GetGlobalMinLayerPosition() : slidingFitResultV.GetGlobalMaxLayerPosition());
91  endV = (isForwardV ? slidingFitResultV.GetGlobalMaxLayerPosition() : slidingFitResultV.GetGlobalMinLayerPosition());
92  }
93 
94  if (foundW)
95  {
96  const TwoDSlidingFitResult &slidingFitResultW = iterW->second;
97  vtxW = (isForwardW ? slidingFitResultW.GetGlobalMinLayerPosition() : slidingFitResultW.GetGlobalMaxLayerPosition());
98  endW = (isForwardW ? slidingFitResultW.GetGlobalMaxLayerPosition() : slidingFitResultW.GetGlobalMinLayerPosition());
99  }
100 
101  CartesianVector vtx3D(0.f, 0.f, 0.f), end3D(0.f, 0.f, 0.f);
102  float vtxChi2(std::numeric_limits<float>::max()), endChi2(std::numeric_limits<float>::max());
103 
104  if (foundU && foundV)
105  {
106  this->UpdateBestPosition(TPC_VIEW_U, TPC_VIEW_V, vtxU, vtxV, vtx3D, vtxChi2);
107  this->UpdateBestPosition(TPC_VIEW_U, TPC_VIEW_V, endU, endV, end3D, endChi2);
108  }
109 
110  if (foundV && foundW)
111  {
112  this->UpdateBestPosition(TPC_VIEW_V, TPC_VIEW_W, vtxV, vtxW, vtx3D, vtxChi2);
113  this->UpdateBestPosition(TPC_VIEW_V, TPC_VIEW_W, endV, endW, end3D, endChi2);
114  }
115 
116  if (foundW && foundU)
117  {
118  this->UpdateBestPosition(TPC_VIEW_W, TPC_VIEW_U, vtxW, vtxU, vtx3D, vtxChi2);
119  this->UpdateBestPosition(TPC_VIEW_W, TPC_VIEW_U, endW, endU, end3D, endChi2);
120  }
121 
122  bool matchedU(false), matchedV(false), matchedW(false);
123  unsigned int matchedViews(0);
124 
125  if (foundU)
126  {
127  const CartesianVector projVtxU(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtx3D, TPC_VIEW_U));
128  const CartesianVector projEndU(LArGeometryHelper::ProjectPosition(this->GetPandora(), end3D, TPC_VIEW_U));
129 
130  if ((endU - vtxU).GetMagnitudeSquared() > m_minTrackLengthSquared &&
131  (projVtxU - vtxU).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projVtxU - endU).GetMagnitudeSquared()) &&
132  (projEndU - endU).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projEndU - vtxU).GetMagnitudeSquared()))
133  {
134  matchedU = true;
135  ++matchedViews;
136  }
137  }
138 
139  if (foundV)
140  {
141  const CartesianVector projVtxV(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtx3D, TPC_VIEW_V));
142  const CartesianVector projEndV(LArGeometryHelper::ProjectPosition(this->GetPandora(), end3D, TPC_VIEW_V));
143 
144  if ((endV - vtxV).GetMagnitudeSquared() > m_minTrackLengthSquared &&
145  (projVtxV - vtxV).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projVtxV - endV).GetMagnitudeSquared()) &&
146  (projEndV - endV).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projEndV - vtxV).GetMagnitudeSquared()))
147  {
148  matchedV = true;
149  ++matchedViews;
150  }
151  }
152 
153  if (foundW)
154  {
155  const CartesianVector projVtxW(LArGeometryHelper::ProjectPosition(this->GetPandora(), vtx3D, TPC_VIEW_W));
156  const CartesianVector projEndW(LArGeometryHelper::ProjectPosition(this->GetPandora(), end3D, TPC_VIEW_W));
157 
158  if ((endW - vtxW).GetMagnitudeSquared() > m_minTrackLengthSquared &&
159  (projVtxW - vtxW).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projVtxW - endW).GetMagnitudeSquared()) &&
160  (projEndW - endW).GetMagnitudeSquared() < std::min(m_vtxDisplacementCutSquared, (projEndW - vtxW).GetMagnitudeSquared()))
161  {
162  matchedW = true;
163  ++matchedViews;
164  }
165  }
166 
167  if (matchedViews < 2)
168  continue;
169 
170  if (vtxChi2 + endChi2 < bestChi2)
171  {
172  useU = matchedU;
173  useV = matchedV;
174  useW = matchedW;
175 
176  bestVtx3D = vtx3D;
177  bestEnd3D = end3D;
178  bestChi2 = vtxChi2 + endChi2;
179  }
180  }
181 
182  if (useU)
183  outputSlidingFitMap.insert(MatchedSlidingFitMap::value_type(iterU->first, iterU->second));
184 
185  if (useV)
186  outputSlidingFitMap.insert(MatchedSlidingFitMap::value_type(iterV->first, iterV->second));
187 
188  if (useW)
189  outputSlidingFitMap.insert(MatchedSlidingFitMap::value_type(iterW->first, iterW->second));
190 
191  if (outputSlidingFitMap.empty())
192  throw StatusCodeException(STATUS_CODE_NOT_FOUND);
193 }
194 
195 //------------------------------------------------------------------------------------------------------------------------------------------
196 
197 void LongitudinalTrackHitsBaseTool::UpdateBestPosition(const HitType hitType1, const HitType hitType2, const CartesianVector &vtx1,
198  const CartesianVector &vtx2, CartesianVector &bestVtx, float &bestChi2) const
199 {
200  CartesianVector mergedVtx(0.f, 0.f, 0.f);
201  float mergedChi2(std::numeric_limits<float>::max());
202 
203  LArGeometryHelper::MergeTwoPositions3D(this->GetPandora(), hitType1, hitType2, vtx1, vtx2, mergedVtx, mergedChi2);
204 
205  if (mergedChi2 < bestChi2)
206  {
207  bestVtx = mergedVtx;
208  bestChi2 = mergedChi2;
209  }
210 }
211 
212 //------------------------------------------------------------------------------------------------------------------------------------------
213 
214 StatusCode LongitudinalTrackHitsBaseTool::ReadSettings(const TiXmlHandle xmlHandle)
215 {
216  float vtxDisplacementCut = std::sqrt(m_vtxDisplacementCutSquared);
217  PANDORA_RETURN_RESULT_IF_AND_IF(
218  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "VertexDisplacementCut", vtxDisplacementCut));
219  m_vtxDisplacementCutSquared = vtxDisplacementCut * vtxDisplacementCut;
220 
221  float minTrackLength = std::sqrt(m_minTrackLengthSquared);
222  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinTrackLength", minTrackLength));
223  m_minTrackLengthSquared = minTrackLength * minTrackLength;
224 
225  return TrackHitsBaseTool::ReadSettings(xmlHandle);
226 }
227 
228 } // namespace lar_content
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Proto hits are temporary constructs to be used during iterative 3D hit procedure. ...
std::map< pandora::HitType, TwoDSlidingFitResult > MatchedSlidingFitMap
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.
virtual void GetLongitudinalTrackHit3D(const MatchedSlidingFitMap &matchedSlidingFitMap, const pandora::CartesianVector &vtx3D, const pandora::CartesianVector &end3D, ProtoHit &protoHit) const =0
Get the three dimensional position using a provided two dimensional calo hit and sliding linear fits ...
virtual void GetTrackHits3D(const pandora::CaloHitVector &inputTwoDHits, const MatchedSlidingFitMap &matchedSlidingFitMap, ProtoHitVector &protoHitVector) const
Calculate 3D hits from an input list of 2D hits.
ThreeDHitCreationAlgorithm::ProtoHitVector ProtoHitVector
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
bool IsPositionSet() const
Whether the proto hit position is set.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
double m_chiSquaredCut
The chi squared cut (accept only values below the cut value)
void UpdateBestPosition(const pandora::HitType hitType1, const pandora::HitType hitType2, const pandora::CartesianVector &vtx1, const pandora::CartesianVector &vtx2, pandora::CartesianVector &bestVtx, float &bestChi2) const
Combine two 2D coordinates to give a 3D coordinate.
TFile f
Definition: plotHisto.C:6
Header file for the three dimensional hit creation algorithm class.
Header file for the geometry helper class.
double GetChi2() const
Get the chi squared value.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
Header file for the longitudinal track hits base tool.
void GetVertexAndEndPositions(const MatchedSlidingFitMap &inputSlidingFitMap, MatchedSlidingFitMap &outputSlidingFitMap, pandora::CartesianVector &outputVtx3D, pandora::CartesianVector &outputEnd3D) const
Get reconstructed vertex and end positions for this 3D track.
HitType
Definition: HitType.h:12
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.