LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
ConnectionPathwayFeatureTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 #include "Pandora/StatusCodes.h"
11 
16 
19 
20 using namespace pandora;
21 
22 namespace lar_content
23 {
24 
25 InitialRegionFeatureTool::InitialRegionFeatureTool() :
26  m_defaultFloat(-10.f),
27  m_nHitsToConsider(10),
28  m_maxInitialGapSizeLimit(4.f),
29  m_minLargestGapSizeLimit(2.f)
30 {
31 }
32 
33 //------------------------------------------------------------------------------------------------------------------------------------------
34 
35 void InitialRegionFeatureTool::Run(LArMvaHelper::MvaFeatureVector &featureVector, const Algorithm *const pAlgorithm,
36  const ParticleFlowObject *const /*pShowerPfo*/, const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch,
37  const CartesianPointVector & /*showerStarts3D*/)
38 {
39  float initialGapSizeU(m_defaultFloat), initialGapSizeV(m_defaultFloat), initialGapSizeW(m_defaultFloat);
40  float largestGapSizeU(m_defaultFloat), largestGapSizeV(m_defaultFloat), largestGapSizeW(m_defaultFloat);
41  this->GetViewInitialRegionVariables(pAlgorithm, nuVertex3D, protoShowerMatch, TPC_VIEW_U, initialGapSizeU, largestGapSizeU);
42  this->GetViewInitialRegionVariables(pAlgorithm, nuVertex3D, protoShowerMatch, TPC_VIEW_V, initialGapSizeV, largestGapSizeV);
43  this->GetViewInitialRegionVariables(pAlgorithm, nuVertex3D, protoShowerMatch, TPC_VIEW_W, initialGapSizeW, largestGapSizeW);
44 
45  float minInitialGapSize(m_defaultFloat), middleInitialGapSize(m_defaultFloat), maxInitialGapSize(m_defaultFloat);
46  float minLargestGapSize(m_defaultFloat), middleLargestGapSize(m_defaultFloat), maxLargestGapSize(m_defaultFloat);
47  LArConnectionPathwayHelper::GetMinMiddleMax(initialGapSizeU, initialGapSizeV, initialGapSizeW, minInitialGapSize, middleInitialGapSize, maxInitialGapSize);
48  LArConnectionPathwayHelper::GetMinMiddleMax(largestGapSizeU, largestGapSizeV, largestGapSizeW, minLargestGapSize, middleLargestGapSize, maxLargestGapSize);
49 
50  featureVector.push_back(std::min(maxInitialGapSize, m_maxInitialGapSizeLimit));
51  featureVector.push_back(std::min(minLargestGapSize, m_minLargestGapSizeLimit));
52 }
53 
54 //------------------------------------------------------------------------------------------------------------------------------------------
55 
56 void InitialRegionFeatureTool::Run(LArMvaHelper::MvaFeatureMap &featureMap, StringVector &featureOrder, const std::string &featureToolName,
57  const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo, const CartesianVector &nuVertex3D,
58  const ProtoShowerMatch &protoShowerMatch, const CartesianPointVector &showerStarts3D)
59 {
60  LArMvaHelper::MvaFeatureVector toolFeatureVec;
61  this->Run(toolFeatureVec, pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, showerStarts3D);
62 
63  if (featureMap.find(featureToolName + "_initialGapSize") != featureMap.end())
64  {
65  std::cout << "Already wrote initialGapSize feature into map! Not writing again." << std::endl;
66  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
67  }
68 
69  featureOrder.push_back(featureToolName + "_initialGapSize");
70  featureMap[featureToolName + "_initialGapSize"] = toolFeatureVec[0];
71 
72  if (featureMap.find(featureToolName + "_largestGapSize") != featureMap.end())
73  {
74  std::cout << "Already wrote largestGapSize feature into map! Not writing again." << std::endl;
75  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
76  }
77 
78  featureOrder.push_back(featureToolName + "_largestGapSize");
79  featureMap[featureToolName + "_largestGapSize"] = toolFeatureVec[1];
80 }
81 
82 //------------------------------------------------------------------------------------------------------------------------------------------
83 
84 void InitialRegionFeatureTool::GetViewInitialRegionVariables(const Algorithm *const pAlgorithm, const CartesianVector &nuVertex3D,
85  const ProtoShowerMatch &protoShowerMatch, const HitType hitType, float &initialGapSize, float &largestGapSize) const
86 {
87  const ProtoShower &protoShower(hitType == TPC_VIEW_U
88  ? protoShowerMatch.GetProtoShowerU()
89  : (hitType == TPC_VIEW_V ? protoShowerMatch.GetProtoShowerV() : protoShowerMatch.GetProtoShowerW()));
90  const CartesianVector nuVertex2D(LArGeometryHelper::ProjectPosition(pAlgorithm->GetPandora(), nuVertex3D, hitType));
91  const CartesianVector &startDirection(protoShower.GetConnectionPathway().GetStartDirection());
92 
93  // Initial gap size
94  CaloHitVector spineHitVector(protoShower.GetSpineHitList().begin(), protoShower.GetSpineHitList().end());
95 
96  std::sort(spineHitVector.begin(), spineHitVector.end(), LArConnectionPathwayHelper::SortByDistanceToPoint(nuVertex2D));
97  initialGapSize = (nuVertex2D - spineHitVector.front()->GetPositionVector()).GetMagnitude();
98 
99  // Largest Gap Size
100  largestGapSize = m_defaultFloat;
101 
102  std::vector<float> longitudinalProjections;
103 
104  for (const CaloHit *const pCaloHit : protoShower.GetSpineHitList())
105  longitudinalProjections.push_back(startDirection.GetDotProduct(pCaloHit->GetPositionVector() - nuVertex2D));
106 
107  std::sort(longitudinalProjections.begin(), longitudinalProjections.end());
108 
109  const unsigned int nIterations(std::min(longitudinalProjections.size(), static_cast<long unsigned int>(m_nHitsToConsider)) - 1);
110 
111  for (unsigned int i = 0; i < nIterations; ++i)
112  largestGapSize = std::max(std::fabs(longitudinalProjections[i] - longitudinalProjections[i + 1]), largestGapSize);
113 }
114 
115 //------------------------------------------------------------------------------------------------------------------------------------------
116 
117 StatusCode InitialRegionFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
118 {
119  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "nHitsToConsider", m_nHitsToConsider));
120 
121  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DefaultFloat", m_defaultFloat));
122 
123  PANDORA_RETURN_RESULT_IF_AND_IF(
124  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxInitialGapSizeLimit", m_maxInitialGapSizeLimit));
125 
126  PANDORA_RETURN_RESULT_IF_AND_IF(
127  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinLargestGapSizeLimit", m_minLargestGapSizeLimit));
128 
129  return STATUS_CODE_SUCCESS;
130 }
131 
132 //------------------------------------------------------------------------------------------------------------------------------------------
133 //------------------------------------------------------------------------------------------------------------------------------------------
134 
136  m_defaultFloat(-10.f),
137  m_spineFitWindow(10),
138  m_nLayersHalfWindow(5),
139  m_pathwayLengthLimit(30.f),
140  m_pathwayScatteringAngle2DLimit(10.f)
141 {
142 }
143 
144 //------------------------------------------------------------------------------------------------------------------------------------------
145 
146 void ConnectionRegionFeatureTool::Run(LArMvaHelper::MvaFeatureVector &featureVector, const Algorithm *const pAlgorithm,
147  const ParticleFlowObject *const /*pShowerPfo*/, const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch,
148  const CartesianPointVector &showerStarts3D)
149 {
150  const float pathwayLength = (nuVertex3D - showerStarts3D.front()).GetMagnitude();
151  const float pathwayScatteringAngle2D = this->Get2DKink(pAlgorithm, protoShowerMatch, showerStarts3D.back());
152 
153  featureVector.push_back(std::min(pathwayLength, m_pathwayLengthLimit));
154  featureVector.push_back(std::min(pathwayScatteringAngle2D, m_pathwayScatteringAngle2DLimit));
155 }
156 
157 //------------------------------------------------------------------------------------------------------------------------------------------
158 
159 void ConnectionRegionFeatureTool::Run(LArMvaHelper::MvaFeatureMap &featureMap, StringVector &featureOrder,
160  const std::string &featureToolName, const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo,
161  const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const CartesianPointVector &showerStarts3D)
162 {
163  LArMvaHelper::MvaFeatureVector toolFeatureVec;
164  this->Run(toolFeatureVec, pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, showerStarts3D);
165 
166  if (featureMap.find(featureToolName + "_pathwayLength") != featureMap.end())
167  {
168  std::cout << "Already wrote pathwayLength feature into map! Not writing again." << std::endl;
169  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
170  }
171 
172  featureOrder.push_back(featureToolName + "_pathwayLength");
173  featureMap[featureToolName + "_pathwayLength"] = toolFeatureVec[0];
174 
175  if (featureMap.find(featureToolName + "_pathwayScatteringAngle2D") != featureMap.end())
176  {
177  std::cout << "Already wrote pathwayScatteringAngle2D feature into map! Not writing again." << std::endl;
178  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
179  }
180 
181  featureOrder.push_back(featureToolName + "_pathwayScatteringAngle2D");
182  featureMap[featureToolName + "_pathwayScatteringAngle2D"] = toolFeatureVec[1];
183 }
184 
185 //------------------------------------------------------------------------------------------------------------------------------------------
186 
188  const Algorithm *const pAlgorithm, const ProtoShowerMatch &protoShowerMatch, const CartesianVector &showerStart3D) const
189 {
190  try
191  {
192  CartesianPointVector spinePositionsU, spinePositionsV, spinePositionsW;
193 
194  for (HitType hitType : {TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W})
195  {
196  const ProtoShower &protoShower(hitType == TPC_VIEW_U
197  ? protoShowerMatch.GetProtoShowerU()
198  : (hitType == TPC_VIEW_V ? protoShowerMatch.GetProtoShowerV() : protoShowerMatch.GetProtoShowerW()));
199  CartesianPointVector &spinePositions(hitType == TPC_VIEW_U ? spinePositionsU : (hitType == TPC_VIEW_V ? spinePositionsV : spinePositionsW));
200 
201  for (const CaloHit *const pCaloHit : protoShower.GetSpineHitList())
202  spinePositions.push_back(pCaloHit->GetPositionVector());
203  }
204 
205  const TwoDSlidingFitResult spineFitU(&spinePositionsU, m_spineFitWindow, LArGeometryHelper::GetWirePitch(pAlgorithm->GetPandora(), TPC_VIEW_U));
206  const TwoDSlidingFitResult spineFitV(&spinePositionsV, m_spineFitWindow, LArGeometryHelper::GetWirePitch(pAlgorithm->GetPandora(), TPC_VIEW_V));
207  const TwoDSlidingFitResult spineFitW(&spinePositionsW, m_spineFitWindow, LArGeometryHelper::GetWirePitch(pAlgorithm->GetPandora(), TPC_VIEW_W));
208 
209  const float scatterAngleU(this->GetLargest2DKinkFromView(pAlgorithm, spineFitU, TPC_VIEW_U, showerStart3D));
210  const float scatterAngleV(this->GetLargest2DKinkFromView(pAlgorithm, spineFitV, TPC_VIEW_V, showerStart3D));
211  const float scatterAngleW(this->GetLargest2DKinkFromView(pAlgorithm, spineFitW, TPC_VIEW_W, showerStart3D));
212 
213  float minScatterAngle(m_defaultFloat), middleScatterAngle(m_defaultFloat), maxScatterAngle(m_defaultFloat);
214  LArConnectionPathwayHelper::GetMinMiddleMax(scatterAngleU, scatterAngleV, scatterAngleW, minScatterAngle, middleScatterAngle, maxScatterAngle);
215 
216  return middleScatterAngle;
217  }
218  catch (...)
219  {
220  }
221 
222  return m_defaultFloat;
223 }
224 
225 //------------------------------------------------------------------------------------------------------------------------------------------
226 
227 float ConnectionRegionFeatureTool::GetLargest2DKinkFromView(const Algorithm *const pAlgorithm, const TwoDSlidingFitResult &spineFit,
228  const HitType hitType, const CartesianVector &showerStart3D) const
229 {
230  const LayerFitResultMap &layerFitResultMap(spineFit.GetLayerFitResultMap());
231  const int minLayer(layerFitResultMap.begin()->first), maxLayer(layerFitResultMap.rbegin()->first);
232  const int nLayersSpanned(1 + maxLayer - minLayer);
233 
234  if (nLayersSpanned <= 2 * m_nLayersHalfWindow)
235  return m_defaultFloat;
236 
237  const CartesianVector projectedShowerStart(LArGeometryHelper::ProjectPosition(pAlgorithm->GetPandora(), showerStart3D, hitType));
238 
239  float showerStartL(0.f), showerStartT(0.f);
240  spineFit.GetLocalPosition(projectedShowerStart, showerStartL, showerStartT);
241 
242  float maxCentralLayer(spineFit.GetLayer(showerStartL) - m_nLayersHalfWindow);
243  float minCentralLayer(layerFitResultMap.begin()->first + m_nLayersHalfWindow + 1);
244 
245  float highestOpeningAngle(m_defaultFloat);
246  CartesianVector kinkPosition(0.f, 0.f, 0.f);
247  CartesianVector highestKinkPosition(0.f, 0.f, 0.f);
248 
249  for (LayerFitResultMap::const_iterator iter = layerFitResultMap.begin(), iterEnd = layerFitResultMap.end(); iter != iterEnd; ++iter)
250  {
251  const int layer(iter->first);
252 
253  if (layer < minCentralLayer)
254  continue;
255 
256  if (layer > maxCentralLayer)
257  continue;
258 
259  bool found(false);
260  float openingAngle2D(std::numeric_limits<float>::max());
261 
262  float thisHighestOpeningAngle(m_defaultFloat);
263  CartesianVector thisHighestKinkPosition(0.f, 0.f, 0.f);
264 
265  for (int i = 0; i < m_nLayersHalfWindow; ++i)
266  {
267  const int testLayer(layer + i);
268  const float rL(spineFit.GetL(testLayer));
269  const float rL1(spineFit.GetL(testLayer - m_nLayersHalfWindow));
270  const float rL2(spineFit.GetL(testLayer + m_nLayersHalfWindow));
271 
272  CartesianVector firstPosition(0.f, 0.f, 0.f), centralPosition(0.f, 0.f, 0.f), secondPosition(0.f, 0.f, 0.f);
273 
274  if ((STATUS_CODE_SUCCESS != spineFit.GetGlobalFitPosition(rL1, firstPosition)) ||
275  (STATUS_CODE_SUCCESS != spineFit.GetGlobalFitPosition(rL, centralPosition)) ||
276  (STATUS_CODE_SUCCESS != spineFit.GetGlobalFitPosition(rL2, secondPosition)))
277  {
278  continue;
279  }
280 
281  const CartesianVector firstDirection((centralPosition - firstPosition).GetUnitVector());
282  const CartesianVector secondDirection((secondPosition - centralPosition).GetUnitVector());
283 
284  const float testOpeningAngle2D(secondDirection.GetOpeningAngle(firstDirection) * 180.f / M_PI);
285 
286  if (testOpeningAngle2D < openingAngle2D)
287  {
288  openingAngle2D = testOpeningAngle2D;
289  found = true;
290  }
291 
292  if (testOpeningAngle2D > thisHighestOpeningAngle)
293  {
294  thisHighestOpeningAngle = openingAngle2D;
295  thisHighestKinkPosition = centralPosition;
296  }
297  }
298 
299  if (found)
300  {
301  if (openingAngle2D > highestOpeningAngle)
302  {
303  highestOpeningAngle = openingAngle2D;
304  highestKinkPosition = thisHighestKinkPosition;
305  }
306  }
307  }
308 
309  return highestOpeningAngle;
310 }
311 
312 //------------------------------------------------------------------------------------------------------------------------------------------
313 
314 StatusCode ConnectionRegionFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
315 {
316  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DefaultFloat", m_defaultFloat));
317 
318  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SpineFitWindow", m_spineFitWindow));
319 
320  PANDORA_RETURN_RESULT_IF_AND_IF(
321  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NLayersHalfWindow", m_nLayersHalfWindow));
322 
323  PANDORA_RETURN_RESULT_IF_AND_IF(
324  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PathwayLengthLimit", m_pathwayLengthLimit));
325 
326  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
327  XmlHelper::ReadValue(xmlHandle, "PathwayScatteringAngle2DLimit", m_pathwayScatteringAngle2DLimit));
328 
329  return STATUS_CODE_SUCCESS;
330 }
331 
332 //------------------------------------------------------------------------------------------------------------------------------------------
333 //------------------------------------------------------------------------------------------------------------------------------------------
334 
336  m_defaultFloat(-10.f),
337  m_defaultRatio(-0.5f),
338  m_spineFitWindow(20),
339  m_showerRadius(14.f),
340  m_showerFitWindow(1000),
341  m_edgeStep(2.f),
342  m_moliereFraction(0.9f),
343  m_maxNHitsLimit(2000.f),
344  m_maxFoundHitRatioLimit(1.5f),
345  m_maxScatterAngleLimit(40.f),
346  m_maxOpeningAngleLimit(20.f),
347  m_maxNuVertexEnergyWeightedMeanRadialDistanceLimit(20.f),
348  m_minShowerStartMoliereRadiusLimit(10.f)
349 {
350 }
351 
352 //------------------------------------------------------------------------------------------------------------------------------------------
353 
354 void ShowerRegionFeatureTool::Run(LArMvaHelper::MvaFeatureVector &featureVector, const Algorithm *const pAlgorithm,
355  const ParticleFlowObject *const pShowerPfo, const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch,
356  const CartesianPointVector &showerStarts3D)
357 {
358  float nHitsU(m_defaultFloat), foundHitRatioU(m_defaultRatio), scatterAngleU(m_defaultFloat), openingAngleU(m_defaultFloat),
359  nuVertexEnergyAsymmetryU(m_defaultRatio), nuVertexEnergyWeightedMeanRadialDistanceU(m_defaultFloat),
360  showerStartEnergyAsymmetryU(m_defaultRatio), showerStartMoliereRadiusU(m_defaultFloat);
361 
362  this->GetViewShowerRegionVariables(pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, TPC_VIEW_U, showerStarts3D.at(0), nHitsU,
363  foundHitRatioU, scatterAngleU, openingAngleU, nuVertexEnergyAsymmetryU, nuVertexEnergyWeightedMeanRadialDistanceU,
364  showerStartEnergyAsymmetryU, showerStartMoliereRadiusU);
365 
366  float nHitsV(m_defaultFloat), foundHitRatioV(m_defaultRatio), scatterAngleV(m_defaultFloat), openingAngleV(m_defaultFloat),
367  nuVertexEnergyAsymmetryV(m_defaultRatio), nuVertexEnergyWeightedMeanRadialDistanceV(m_defaultFloat),
368  showerStartEnergyAsymmetryV(m_defaultRatio), showerStartMoliereRadiusV(m_defaultFloat);
369 
370  this->GetViewShowerRegionVariables(pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, TPC_VIEW_V, showerStarts3D.at(0), nHitsV,
371  foundHitRatioV, scatterAngleV, openingAngleV, nuVertexEnergyAsymmetryV, nuVertexEnergyWeightedMeanRadialDistanceV,
372  showerStartEnergyAsymmetryV, showerStartMoliereRadiusV);
373 
374  float nHitsW(m_defaultFloat), foundHitRatioW(m_defaultRatio), scatterAngleW(m_defaultFloat), openingAngleW(m_defaultFloat),
375  nuVertexEnergyAsymmetryW(m_defaultRatio), nuVertexEnergyWeightedMeanRadialDistanceW(m_defaultFloat),
376  showerStartEnergyAsymmetryW(m_defaultRatio), showerStartMoliereRadiusW(m_defaultFloat);
377 
378  this->GetViewShowerRegionVariables(pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, TPC_VIEW_W, showerStarts3D.at(0), nHitsW,
379  foundHitRatioW, scatterAngleW, openingAngleW, nuVertexEnergyAsymmetryW, nuVertexEnergyWeightedMeanRadialDistanceW,
380  showerStartEnergyAsymmetryW, showerStartMoliereRadiusW);
381 
382  float minNHits(m_defaultFloat), minFoundHitRatio(m_defaultRatio), minScatterAngle(m_defaultFloat), minOpeningAngle(m_defaultFloat),
383  minNuVertexEnergyAsymmetry(m_defaultRatio), minNuVertexEnergyWeightedMeanRadialDistance(m_defaultFloat),
384  minShowerStartEnergyAsymmetry(m_defaultRatio), minShowerStartMoliereRadius(m_defaultFloat);
385 
386  float middleNHits(m_defaultFloat), middleFoundHitRatio(m_defaultRatio), middleScatterAngle(m_defaultFloat), middleOpeningAngle(m_defaultFloat),
387  middleNuVertexEnergyAsymmetry(m_defaultRatio), middleNuVertexEnergyWeightedMeanRadialDistance(m_defaultFloat),
388  middleShowerStartEnergyAsymmetry(m_defaultRatio), middleShowerStartMoliereRadius(m_defaultFloat);
389 
390  float maxNHits(m_defaultFloat), maxFoundHitRatio(m_defaultRatio), maxScatterAngle(m_defaultFloat), maxOpeningAngle(m_defaultFloat),
391  maxNuVertexEnergyAsymmetry(m_defaultRatio), maxNuVertexEnergyWeightedMeanRadialDistance(m_defaultFloat),
392  maxShowerStartEnergyAsymmetry(m_defaultRatio), maxShowerStartMoliereRadius(m_defaultFloat);
393 
394  LArConnectionPathwayHelper::GetMinMiddleMax(nHitsU, nHitsV, nHitsW, minNHits, middleNHits, maxNHits);
395  LArConnectionPathwayHelper::GetMinMiddleMax(foundHitRatioU, foundHitRatioV, foundHitRatioW, minFoundHitRatio, middleFoundHitRatio, maxFoundHitRatio);
396  LArConnectionPathwayHelper::GetMinMiddleMax(scatterAngleU, scatterAngleV, scatterAngleW, minScatterAngle, middleScatterAngle, maxScatterAngle);
397  LArConnectionPathwayHelper::GetMinMiddleMax(openingAngleU, openingAngleV, openingAngleW, minOpeningAngle, middleOpeningAngle, maxOpeningAngle);
398  LArConnectionPathwayHelper::GetMinMiddleMax(nuVertexEnergyAsymmetryU, nuVertexEnergyAsymmetryV, nuVertexEnergyAsymmetryW,
399  minNuVertexEnergyAsymmetry, middleNuVertexEnergyAsymmetry, maxNuVertexEnergyAsymmetry);
400  LArConnectionPathwayHelper::GetMinMiddleMax(nuVertexEnergyWeightedMeanRadialDistanceU, nuVertexEnergyWeightedMeanRadialDistanceV,
401  nuVertexEnergyWeightedMeanRadialDistanceW, minNuVertexEnergyWeightedMeanRadialDistance,
402  middleNuVertexEnergyWeightedMeanRadialDistance, maxNuVertexEnergyWeightedMeanRadialDistance);
403  LArConnectionPathwayHelper::GetMinMiddleMax(showerStartEnergyAsymmetryU, showerStartEnergyAsymmetryV, showerStartEnergyAsymmetryW,
404  minShowerStartEnergyAsymmetry, middleShowerStartEnergyAsymmetry, maxShowerStartEnergyAsymmetry);
405  LArConnectionPathwayHelper::GetMinMiddleMax(showerStartMoliereRadiusU, showerStartMoliereRadiusV, showerStartMoliereRadiusW,
406  minShowerStartMoliereRadius, middleShowerStartMoliereRadius, maxShowerStartMoliereRadius);
407 
408  featureVector.push_back(std::min(maxNHits, m_maxNHitsLimit));
409  featureVector.push_back(std::min(maxFoundHitRatio, m_maxFoundHitRatioLimit));
410  featureVector.push_back(std::min(maxScatterAngle, m_maxScatterAngleLimit));
411  featureVector.push_back(std::min(maxOpeningAngle, m_maxOpeningAngleLimit));
412  featureVector.push_back(maxNuVertexEnergyAsymmetry);
413  featureVector.push_back(std::min(maxNuVertexEnergyWeightedMeanRadialDistance, m_maxNuVertexEnergyWeightedMeanRadialDistanceLimit));
414  featureVector.push_back(maxShowerStartEnergyAsymmetry);
415  featureVector.push_back(std::min(minShowerStartMoliereRadius, m_minShowerStartMoliereRadiusLimit));
416 }
417 
418 //------------------------------------------------------------------------------------------------------------------------------------------
419 
420 void ShowerRegionFeatureTool::Run(LArMvaHelper::MvaFeatureMap &featureMap, StringVector &featureOrder, const std::string &featureToolName,
421  const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo, const CartesianVector &nuVertex3D,
422  const ProtoShowerMatch &protoShowerMatch, const CartesianPointVector &showerStarts3D)
423 {
424  LArMvaHelper::MvaFeatureVector toolFeatureVec;
425  this->Run(toolFeatureVec, pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, showerStarts3D);
426 
427  if (featureMap.find(featureToolName + "_nShowerHits") != featureMap.end())
428  {
429  std::cout << "Already wrote nShowerHits feature into map! Not writing again." << std::endl;
430  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
431  }
432 
433  featureOrder.push_back(featureToolName + "_nShowerHits");
434  featureMap[featureToolName + "_nShowerHits"] = toolFeatureVec[0];
435 
436  if (featureMap.find(featureToolName + "_foundHitRatio") != featureMap.end())
437  {
438  std::cout << "Already wrote foundHitRatio feature into map! Not writing again." << std::endl;
439  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
440  }
441 
442  featureOrder.push_back(featureToolName + "_foundHitRatio");
443  featureMap[featureToolName + "_foundHitRatio"] = toolFeatureVec[1];
444 
445  if (featureMap.find(featureToolName + "_scatterAngle") != featureMap.end())
446  {
447  std::cout << "Already wrote scatterAngle feature into map! Not writing again." << std::endl;
448  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
449  }
450 
451  featureOrder.push_back(featureToolName + "_scatterAngle");
452  featureMap[featureToolName + "_scatterAngle"] = toolFeatureVec[2];
453 
454  if (featureMap.find(featureToolName + "_openingAngle") != featureMap.end())
455  {
456  std::cout << "Already wrote openingAngle feature into map! Not writing again." << std::endl;
457  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
458  }
459 
460  featureOrder.push_back(featureToolName + "_openingAngle");
461  featureMap[featureToolName + "_openingAngle"] = toolFeatureVec[3];
462 
463  if (featureMap.find(featureToolName + "_nuVertexEnergyAsymmetry") != featureMap.end())
464  {
465  std::cout << "Already wrote nuVertexEnergyAsymmetry feature into map! Not writing again." << std::endl;
466  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
467  }
468 
469  featureOrder.push_back(featureToolName + "_nuVertexEnergyAsymmetry");
470  featureMap[featureToolName + "_nuVertexEnergyAsymmetry"] = toolFeatureVec[4];
471 
472  if (featureMap.find(featureToolName + "_nuVertexEnergyWeightedMeanRadialDistance") != featureMap.end())
473  {
474  std::cout << "Already wrote nuVertexEnergyWeightedMeanRadialDistance feature into map! Not writing again." << std::endl;
475  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
476  }
477 
478  featureOrder.push_back(featureToolName + "_nuVertexEnergyWeightedMeanRadialDistance");
479  featureMap[featureToolName + "_nuVertexEnergyWeightedMeanRadialDistance"] = toolFeatureVec[5];
480 
481  if (featureMap.find(featureToolName + "_showerStartEnergyAsymmetry") != featureMap.end())
482  {
483  std::cout << "Already wrote showerStartEnergyAsymmetry feature into map! Not writing again." << std::endl;
484  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
485  }
486 
487  featureOrder.push_back(featureToolName + "_showerStartEnergyAsymmetry");
488  featureMap[featureToolName + "_showerStartEnergyAsymmetry"] = toolFeatureVec[6];
489 
490  if (featureMap.find(featureToolName + "_showerStartMoliereRadius") != featureMap.end())
491  {
492  std::cout << "Already wrote showerStartMoliereRadius feature into map! Not writing again." << std::endl;
493  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
494  }
495 
496  featureOrder.push_back(featureToolName + "_showerStartMoliereRadius");
497  featureMap[featureToolName + "_showerStartMoliereRadius"] = toolFeatureVec[7];
498 }
499 
500 //------------------------------------------------------------------------------------------------------------------------------------------
501 
502 void ShowerRegionFeatureTool::GetViewShowerRegionVariables(const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo,
503  const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const HitType hitType, const CartesianVector &showerStart3D,
504  float &nHits, float &foundHitRatio, float &scatterAngle, float &openingAngle, float &nuVertexEnergyAsymmetry,
505  float &nuVertexEnergyWeightedMeanRadialDistance, float &showerStartEnergyAsymmetry, float &showerStartMoliereRadius)
506 {
507  CaloHitList viewHitList;
508  LArPfoHelper::GetCaloHits(pShowerPfo, hitType, viewHitList);
509 
510  const CartesianVector nuVertex2D(LArGeometryHelper::ProjectPosition(pAlgorithm->GetPandora(), nuVertex3D, hitType));
511  const CartesianVector showerStart2D(LArGeometryHelper::ProjectPosition(pAlgorithm->GetPandora(), showerStart3D, hitType));
512  const bool isDownstream(showerStart2D.GetZ() > nuVertex2D.GetZ());
513 
514  try
515  {
516  // Fit the spine to get the initial shower direction
517  const CaloHitList &spineHitList(hitType == TPC_VIEW_U ? protoShowerMatch.GetProtoShowerU().GetSpineHitList()
518  : (hitType == TPC_VIEW_V ? protoShowerMatch.GetProtoShowerV().GetSpineHitList()
519  : protoShowerMatch.GetProtoShowerW().GetSpineHitList()));
520 
521  CartesianPointVector spinePositions;
522  for (const CaloHit *const pCaloHit : spineHitList)
523  spinePositions.push_back(pCaloHit->GetPositionVector());
524 
525  const TwoDSlidingFitResult spineFitResult(&spinePositions, m_spineFitWindow, LArGeometryHelper::GetWirePitch(pAlgorithm->GetPandora(), hitType));
526 
527  // Now can build the shower
528  CaloHitList postShowerHitList;
529  CartesianPointVector postShowerPositions;
530 
531  this->BuildViewShower(pShowerPfo, spineFitResult, hitType, showerStart2D, nuVertex2D, postShowerHitList, postShowerPositions);
532 
533  this->GetShowerHitVariables(spineHitList, postShowerHitList, pShowerPfo, hitType, nHits, foundHitRatio);
534 
535  // Fit the shower
536  TwoDSlidingFitResult showerFitResult(
537  &postShowerPositions, m_showerFitWindow, LArGeometryHelper::GetWirePitch(pAlgorithm->GetPandora(), hitType));
538 
539  const bool isShowerDownstream((showerStart2D - showerFitResult.GetGlobalMinLayerPosition()).GetMagnitude() <
540  (showerStart2D - showerFitResult.GetGlobalMaxLayerPosition()).GetMagnitude());
541 
542  // Collect more hits
543  for (const CaloHit *const pCaloHit : viewHitList)
544  {
545  if (std::find(postShowerHitList.begin(), postShowerHitList.end(), pCaloHit) != postShowerHitList.end())
546  continue;
547 
548  const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
549  const CartesianVector displacement(hitPosition - showerStart2D);
550  const float l(showerFitResult.GetGlobalMinLayerDirection().GetDotProduct(displacement));
551  const float t(showerFitResult.GetGlobalMinLayerDirection().GetCrossProduct(displacement).GetMagnitude());
552 
553  if (((isDownstream && (l > 0.f)) || (!isDownstream && (l < 0.f))) && (t < m_showerRadius))
554  {
555  postShowerHitList.push_back(pCaloHit);
556  postShowerPositions.push_back(pCaloHit->GetPositionVector());
557  }
558  }
559 
560  this->GetShowerHitVariables(spineHitList, postShowerHitList, pShowerPfo, hitType, nHits, foundHitRatio);
561 
562  this->CalculateViewScatterAngle(nuVertex2D, spineFitResult, showerStart2D, showerFitResult, scatterAngle);
563 
564  this->CalculateViewOpeningAngle(showerFitResult, postShowerHitList, showerStart2D, openingAngle);
565 
567  spineFitResult, postShowerHitList, isDownstream, nuVertex2D, nuVertexEnergyAsymmetry, nuVertexEnergyWeightedMeanRadialDistance);
568 
570  showerFitResult, postShowerHitList, isShowerDownstream, showerStartEnergyAsymmetry, showerStartMoliereRadius);
571  }
572  catch (...)
573  {
574  }
575 }
576 
577 //------------------------------------------------------------------------------------------------------------------------------------------
578 
579 void ShowerRegionFeatureTool::BuildViewShower(const ParticleFlowObject *const pShowerPfo, const TwoDSlidingFitResult &spineFit, const HitType hitType,
580  const CartesianVector &showerStart2D, const CartesianVector &nuVertex2D, CaloHitList &postShowerHitList, CartesianPointVector &postShowerPositions)
581 {
582  // Find initial shower direction from spine fit
583  float lShowerStart(0.f), tShowerStart(0.f);
584  spineFit.GetLocalPosition(showerStart2D, lShowerStart, tShowerStart);
585 
586  const LayerFitResultMap &layerFitResultMap(spineFit.GetLayerFitResultMap());
587  const int showerStartLayer(spineFit.GetLayer(lShowerStart));
588  int closestLayer(std::numeric_limits<int>::max());
589 
590  for (const auto &entry : layerFitResultMap)
591  {
592  if (std::fabs(entry.first - showerStartLayer) < std::fabs(entry.first - closestLayer))
593  closestLayer = entry.first;
594  }
595 
596  const float gradient(layerFitResultMap.at(closestLayer).GetGradient());
597  CartesianVector showerDirection(0.f, 0.f, 0.f);
598 
599  spineFit.GetGlobalDirection(gradient, showerDirection);
600 
601  // Now find the shower
602  const bool isDownstream(showerStart2D.GetZ() > nuVertex2D.GetZ());
603  CaloHitList caloHitList;
604 
605  LArPfoHelper::GetCaloHits(pShowerPfo, hitType, caloHitList);
606 
607  for (const CaloHit *const pCaloHit : caloHitList)
608  {
609  const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
610  const CartesianVector displacement(hitPosition - showerStart2D);
611  const float l(showerDirection.GetDotProduct(displacement));
612  const float t(showerDirection.GetCrossProduct(displacement).GetMagnitude());
613 
614  if (((isDownstream && (l > 0.f)) || (!isDownstream && (l < 0.f))) && (t < m_showerRadius))
615  {
616  postShowerHitList.push_back(pCaloHit);
617  postShowerPositions.push_back(pCaloHit->GetPositionVector());
618  }
619  }
620 }
621 
622 //------------------------------------------------------------------------------------------------------------------------------------------
623 
624 void ShowerRegionFeatureTool::GetShowerHitVariables(const CaloHitList &spineHitList, const CaloHitList &postShowerHitList,
625  const ParticleFlowObject *const pShowerPfo, const HitType hitType, float &nHits, float &foundHitRatio)
626 {
627  int foundHits(spineHitList.size());
628 
629  for (const CaloHit *const pCaloHit : postShowerHitList)
630  {
631  if (std::find(spineHitList.begin(), spineHitList.end(), pCaloHit) == spineHitList.end())
632  ++foundHits;
633  }
634 
635  CaloHitList viewHitList;
636  LArPfoHelper::GetCaloHits(pShowerPfo, hitType, viewHitList);
637 
638  foundHitRatio = static_cast<float>(foundHits) / static_cast<float>(viewHitList.size());
639  nHits = postShowerHitList.size();
640 }
641 
642 //------------------------------------------------------------------------------------------------------------------------------------------
643 
644 void ShowerRegionFeatureTool::CalculateViewScatterAngle(const CartesianVector &nuVertex2D, const TwoDSlidingFitResult &spineFitResult,
645  const CartesianVector &showerStart2D, const TwoDSlidingFitResult &showerFitResult, float &scatterAngle)
646 {
647  const bool isDownstream(showerStart2D.GetZ() > nuVertex2D.GetZ());
648  const CartesianVector streamCorrectedConnectionPathwayDirection(
649  isDownstream ? spineFitResult.GetGlobalMinLayerDirection() : spineFitResult.GetGlobalMaxLayerDirection() * (-1.f));
650 
651  const bool isShowerDownstream((showerStart2D - showerFitResult.GetGlobalMinLayerPosition()).GetMagnitude() <
652  (showerStart2D - showerFitResult.GetGlobalMaxLayerPosition()).GetMagnitude());
653  const CartesianVector streamCorrectedShowerDirection(
654  isShowerDownstream ? showerFitResult.GetGlobalMinLayerDirection() : showerFitResult.GetGlobalMaxLayerDirection() * (-1.f));
655 
656  scatterAngle = streamCorrectedConnectionPathwayDirection.GetOpeningAngle(streamCorrectedShowerDirection) * 180.f / M_PI;
657 }
658 
659 //------------------------------------------------------------------------------------------------------------------------------------------
660 
661 void ShowerRegionFeatureTool::CalculateViewOpeningAngle(const TwoDSlidingFitResult &showerFitResult, const CaloHitList &postShowerHitList,
662  const CartesianVector &showerStart2D, float &openingAngle)
663 {
664  try
665  {
666  const CartesianVector directionAxis = showerFitResult.GetGlobalMinLayerDirection();
667  const CartesianVector orthoAxis = directionAxis.GetCrossProduct(CartesianVector(0.f, 1.f, 0.f));
668 
669  std::map<int, float> positiveEdges, negativeEdges;
670 
671  for (const CaloHit *const pCaloHit : postShowerHitList)
672  {
673  const CartesianVector position(pCaloHit->GetPositionVector() - showerStart2D);
674  const float thisT(directionAxis.GetCrossProduct(position).GetMagnitude());
675  const float thisL(directionAxis.GetDotProduct(position));
676  const float orthoL(orthoAxis.GetDotProduct(position));
677 
678  std::map<int, float> &edgeMap(orthoL > 0.f ? positiveEdges : negativeEdges);
679 
680  const int lIndex(std::floor(thisL / m_edgeStep));
681 
682  edgeMap[lIndex] = (edgeMap.find(lIndex) == edgeMap.end() ? thisT : std::max(edgeMap[lIndex], thisT));
683  }
684 
685  CartesianPointVector positiveEdgePositions, negativeEdgePositions;
686 
687  for (auto &entry : positiveEdges)
688  positiveEdgePositions.push_back(CartesianVector(entry.second, 0.f, entry.first));
689 
690  for (auto &entry : negativeEdges)
691  negativeEdgePositions.push_back(CartesianVector(entry.second, 0.f, entry.first));
692 
693  const TwoDSlidingFitResult positiveEdgeFit(&positiveEdgePositions, m_showerFitWindow, showerFitResult.GetLayerPitch());
694  const TwoDSlidingFitResult negativeEdgeFit(&negativeEdgePositions, m_showerFitWindow, showerFitResult.GetLayerPitch());
695 
696  const CartesianVector positiveMinLayer(positiveEdgeFit.GetGlobalMinLayerPosition());
697  const CartesianVector positiveMaxLayer(positiveEdgeFit.GetGlobalMaxLayerPosition());
698  const CartesianVector negativeMinLayer(negativeEdgeFit.GetGlobalMinLayerPosition());
699  const CartesianVector negativeMaxLayer(negativeEdgeFit.GetGlobalMaxLayerPosition());
700 
701  const CartesianVector globalPositiveMinLayer(
702  showerStart2D + (directionAxis * positiveMinLayer.GetZ()) + (orthoAxis * positiveMinLayer.GetX()));
703  const CartesianVector globalPositiveMaxLayer(
704  showerStart2D + (directionAxis * positiveMaxLayer.GetZ()) + (orthoAxis * positiveMaxLayer.GetX()));
705  const CartesianVector globalNegativeMinLayer(
706  showerStart2D + (directionAxis * negativeMinLayer.GetZ()) - (orthoAxis * negativeMinLayer.GetX()));
707  const CartesianVector globalNegativeMaxLayer(
708  showerStart2D + (directionAxis * negativeMaxLayer.GetZ()) - (orthoAxis * negativeMaxLayer.GetX()));
709 
710  const CartesianVector positiveEdgeVector((globalPositiveMaxLayer - globalPositiveMinLayer).GetUnitVector());
711  const CartesianVector negativeEdgeVector((globalNegativeMaxLayer - globalNegativeMinLayer).GetUnitVector());
712 
713  const float positiveOpeningAngle = directionAxis.GetOpeningAngle(positiveEdgeVector) * 180.f / M_PI;
714  const float negativeOpeningAngle = directionAxis.GetOpeningAngle(negativeEdgeVector) * 180.f / M_PI;
715 
716  openingAngle = std::max(positiveOpeningAngle, negativeOpeningAngle);
717  }
718  catch (...)
719  {
720  }
721 }
722 
723 //------------------------------------------------------------------------------------------------------------------------------------------
724 
725 void ShowerRegionFeatureTool::CalculateViewNuVertexConsistencyVariables(const TwoDSlidingFitResult &spineFitResult, const CaloHitList &postShowerHitList,
726  const bool isDownstream, const CartesianVector &nuVertex2D, float &nuVertexEnergyAsymmetry, float &nuVertexEnergyWeightedMeanRadialDistance)
727 {
728  const CartesianVector &directionAxis(isDownstream ? spineFitResult.GetGlobalMinLayerDirection() : spineFitResult.GetGlobalMaxLayerDirection());
729  const CartesianVector orthoAxis = directionAxis.GetCrossProduct(CartesianVector(0.f, 1.f, 0.f));
730 
731  float totalEnergy(0.f);
732 
733  // nuVertexEnergyAsymmetry
734  nuVertexEnergyAsymmetry = 0.f;
735 
736  for (const CaloHit *const pCaloHit : postShowerHitList)
737  {
738  const float hitEnergy(std::fabs(pCaloHit->GetElectromagneticEnergy()));
739 
740  totalEnergy += hitEnergy;
741 
742  const CartesianVector position(pCaloHit->GetPositionVector() - nuVertex2D);
743  const float thisL(orthoAxis.GetDotProduct(position));
744 
745  nuVertexEnergyAsymmetry += (thisL < 0.f) ? (-1.f * hitEnergy) : hitEnergy;
746  }
747 
748  nuVertexEnergyAsymmetry = (totalEnergy < std::numeric_limits<float>::epsilon()) ? m_defaultRatio : (nuVertexEnergyAsymmetry / totalEnergy);
749  nuVertexEnergyAsymmetry = std::fabs(nuVertexEnergyAsymmetry);
750 
751  // nuVertexEnergyWeightedMeanRadialDistance
752  nuVertexEnergyWeightedMeanRadialDistance = 0.f;
753 
754  for (const CaloHit *const pCaloHit : postShowerHitList)
755  {
756  const CartesianVector position(pCaloHit->GetPositionVector() - nuVertex2D);
757  const float hitEnergy(std::fabs(pCaloHit->GetElectromagneticEnergy()));
758 
759  nuVertexEnergyWeightedMeanRadialDistance += (directionAxis.GetCrossProduct(position).GetMagnitude() * hitEnergy);
760  }
761 
762  nuVertexEnergyWeightedMeanRadialDistance =
763  (totalEnergy < std::numeric_limits<float>::epsilon()) ? m_defaultFloat : nuVertexEnergyWeightedMeanRadialDistance / totalEnergy;
764 }
765 
766 //------------------------------------------------------------------------------------------------------------------------------------------
767 
769  const CaloHitList &postShowerHitList, const bool isShowerDownstream, float &showerStartEnergyAsymmetry, float &showerStartMoliereRadius)
770 {
771  const CartesianVector fitShowerStart(isShowerDownstream ? showerFitResult.GetGlobalMinLayerPosition() : showerFitResult.GetGlobalMaxLayerPosition());
772  const CartesianVector directionAxis(
773  isShowerDownstream ? showerFitResult.GetGlobalMinLayerDirection() : showerFitResult.GetGlobalMaxLayerDirection());
774  const CartesianVector orthoAxis = directionAxis.GetCrossProduct(CartesianVector(0.f, 1.f, 0.f));
775 
776  float totalEnergy(0.f);
777 
778  // showerStartEnergyAsymmetry
779  showerStartEnergyAsymmetry = 0.f;
780 
781  for (const CaloHit *const pCaloHit : postShowerHitList)
782  {
783  const float hitEnergy(std::fabs(pCaloHit->GetElectromagneticEnergy()));
784 
785  totalEnergy += hitEnergy;
786 
787  const CartesianVector position(pCaloHit->GetPositionVector() - fitShowerStart);
788  const float thisL(orthoAxis.GetDotProduct(position));
789 
790  showerStartEnergyAsymmetry += (thisL < 0.f) ? (-1.f * hitEnergy) : hitEnergy;
791  }
792 
793  showerStartEnergyAsymmetry = (totalEnergy < std::numeric_limits<float>::epsilon()) ? m_defaultRatio : (showerStartEnergyAsymmetry / totalEnergy);
794  showerStartEnergyAsymmetry = std::fabs(showerStartEnergyAsymmetry);
795 
796  // showerStartMoliereRadius
797  showerStartMoliereRadius = m_defaultFloat;
798 
799  CaloHitVector showerStartPostShowerHitVector(postShowerHitList.begin(), postShowerHitList.end());
800 
801  std::sort(showerStartPostShowerHitVector.begin(), showerStartPostShowerHitVector.end(),
802  [&fitShowerStart, &directionAxis](const CaloHit *const pCaloHitA, const CaloHit *const pCaloHitB) -> bool
803  {
804  const CartesianVector positionA(pCaloHitA->GetPositionVector() - fitShowerStart);
805  const CartesianVector positionB(pCaloHitB->GetPositionVector() - fitShowerStart);
806 
807  const float tA(directionAxis.GetCrossProduct(positionA).GetMagnitude());
808  const float tB(directionAxis.GetCrossProduct(positionB).GetMagnitude());
809 
810  return tA < tB;
811  });
812 
813  float showerStartRunningEnergySum(0.f);
814 
815  for (const CaloHit *const pCaloHit : showerStartPostShowerHitVector)
816  {
817  const float hitEnergy(std::fabs(pCaloHit->GetElectromagneticEnergy()));
818  showerStartRunningEnergySum += hitEnergy;
819 
820  if ((showerStartRunningEnergySum / totalEnergy) > m_moliereFraction)
821  {
822  const CartesianVector position(pCaloHit->GetPositionVector() - fitShowerStart);
823  showerStartMoliereRadius = directionAxis.GetCrossProduct(position).GetMagnitude();
824  break;
825  }
826  }
827 }
828 
829 //------------------------------------------------------------------------------------------------------------------------------------------
830 
831 StatusCode ShowerRegionFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
832 {
833  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DefaultFloat", m_defaultFloat));
834 
835  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DefaultRatio", m_defaultRatio));
836 
837  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SpineFitWindow", m_spineFitWindow));
838 
839  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ShowerRadius", m_showerRadius));
840 
841  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ShowerFitWindow", m_showerFitWindow));
842 
843  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "EdgeStep", m_edgeStep));
844 
845  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MoliereFraction", m_moliereFraction));
846 
847  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxNHitsLimit", m_maxNHitsLimit));
848 
849  PANDORA_RETURN_RESULT_IF_AND_IF(
850  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxFoundHitRatioLimit", m_maxFoundHitRatioLimit));
851 
852  PANDORA_RETURN_RESULT_IF_AND_IF(
853  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxScatterAngleLimit", m_maxScatterAngleLimit));
854 
855  PANDORA_RETURN_RESULT_IF_AND_IF(
856  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxOpeningAngleLimit", m_maxOpeningAngleLimit));
857 
858  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
859  XmlHelper::ReadValue(xmlHandle, "MaxNuVertexEnergyWeightedMeanRadialDistanceLimit", m_maxNuVertexEnergyWeightedMeanRadialDistanceLimit));
860 
861  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
862  XmlHelper::ReadValue(xmlHandle, "MinShowerStartMoliereRadiusLimit", m_minShowerStartMoliereRadiusLimit));
863 
864  return STATUS_CODE_SUCCESS;
865 }
866 
867 //------------------------------------------------------------------------------------------------------------------------------------------
868 //------------------------------------------------------------------------------------------------------------------------------------------
869 
871  m_defaultFloat(-10.f),
872  m_caloHitListNameU("CaloHitListU"),
873  m_caloHitListNameV("CaloHitListV"),
874  m_caloHitListNameW("CaloHitListW"),
875  m_maxTransverseDistance(0.75f),
876  m_maxSampleHits(3),
877  m_maxHitSeparation(1.f),
878  m_maxTrackFraction(0.8f)
879 {
880 }
881 
882 //------------------------------------------------------------------------------------------------------------------------------------------
883 
884 void AmbiguousRegionFeatureTool::Run(LArMvaHelper::MvaFeatureVector &featureVector, const Algorithm *const pAlgorithm,
885  const ParticleFlowObject *const /*pShowerPfo*/, const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch,
886  const CartesianPointVector & /*showerStarts3D*/)
887 {
888  float nAmbiguousViews(0.f);
889  this->CalculateNAmbiguousViews(protoShowerMatch, nAmbiguousViews);
890 
891  float maxUnaccountedEnergy(m_defaultFloat);
892  float unaccountedHitEnergyU(m_defaultFloat), unaccountedHitEnergyV(m_defaultFloat), unaccountedHitEnergyW(m_defaultFloat);
893 
894  if (this->GetViewAmbiguousHitVariables(pAlgorithm, protoShowerMatch, TPC_VIEW_U, nuVertex3D, unaccountedHitEnergyU))
895  maxUnaccountedEnergy = std::max(maxUnaccountedEnergy, unaccountedHitEnergyU);
896 
897  if (this->GetViewAmbiguousHitVariables(pAlgorithm, protoShowerMatch, TPC_VIEW_V, nuVertex3D, unaccountedHitEnergyV))
898  maxUnaccountedEnergy = std::max(maxUnaccountedEnergy, unaccountedHitEnergyV);
899 
900  if (this->GetViewAmbiguousHitVariables(pAlgorithm, protoShowerMatch, TPC_VIEW_W, nuVertex3D, unaccountedHitEnergyW))
901  maxUnaccountedEnergy = std::max(maxUnaccountedEnergy, unaccountedHitEnergyW);
902 
903  featureVector.push_back(nAmbiguousViews);
904  featureVector.push_back(maxUnaccountedEnergy);
905 }
906 
907 //------------------------------------------------------------------------------------------------------------------------------------------
908 
909 void AmbiguousRegionFeatureTool::Run(LArMvaHelper::MvaFeatureMap &featureMap, StringVector &featureOrder,
910  const std::string &featureToolName, const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo,
911  const CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const CartesianPointVector &showerStarts3D)
912 {
913  LArMvaHelper::MvaFeatureVector toolFeatureVec;
914  this->Run(toolFeatureVec, pAlgorithm, pShowerPfo, nuVertex3D, protoShowerMatch, showerStarts3D);
915 
916  if (featureMap.find(featureToolName + "_nAmbiguousViews") != featureMap.end())
917  {
918  std::cout << "Already wrote nAmbiguousViews feature into map! Not writing again." << std::endl;
919  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
920  }
921 
922  featureOrder.push_back(featureToolName + "_nAmbiguousViews");
923  featureMap[featureToolName + "_nAmbiguousViews"] = toolFeatureVec[0];
924 
925  if (featureMap.find(featureToolName + "_maxUnaccountedEnergy") != featureMap.end())
926  {
927  std::cout << "Already wrote maxUnaccountedEnergy feature into map! Not writing again." << std::endl;
928  throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
929  }
930 
931  featureOrder.push_back(featureToolName + "_maxUnaccountedEnergy");
932  featureMap[featureToolName + "_maxUnaccountedEnergy"] = toolFeatureVec[1];
933 }
934 
935 //------------------------------------------------------------------------------------------------------------------------------------------
936 
937 void AmbiguousRegionFeatureTool::CalculateNAmbiguousViews(const ProtoShowerMatch &protoShowerMatch, float &nAmbiguousViews)
938 {
939  nAmbiguousViews = 0.f;
940 
941  const int nAmbiguousHitsU(protoShowerMatch.GetProtoShowerU().GetAmbiguousHitList().size());
942  nAmbiguousViews += (nAmbiguousHitsU == 0) ? 0.f : 1.f;
943 
944  const int nAmbiguousHitsV(protoShowerMatch.GetProtoShowerV().GetAmbiguousHitList().size());
945  nAmbiguousViews += (nAmbiguousHitsV == 0) ? 0.f : 1.f;
946 
947  const int nAmbiguousHitsW(protoShowerMatch.GetProtoShowerW().GetAmbiguousHitList().size());
948  nAmbiguousViews += (nAmbiguousHitsW == 0) ? 0.f : 1.f;
949 }
950 
951 //------------------------------------------------------------------------------------------------------------------------------------------
952 
953 bool AmbiguousRegionFeatureTool::GetViewAmbiguousHitVariables(const Algorithm *const pAlgorithm, const ProtoShowerMatch &protoShowerMatch,
954  const HitType hitType, const CartesianVector &nuVertex3D, float &unaccountedHitEnergy)
955 {
956  std::map<int, CaloHitList> ambiguousHitSpines;
957  CaloHitList hitsToExcludeInEnergyCalcs; // to avoid double counting
958  const CartesianVector nuVertex2D(LArGeometryHelper::ProjectPosition(pAlgorithm->GetPandora(), nuVertex3D, hitType));
959  const ProtoShower &protoShower(hitType == TPC_VIEW_U
960  ? protoShowerMatch.GetProtoShowerU()
961  : (hitType == TPC_VIEW_V ? protoShowerMatch.GetProtoShowerV() : protoShowerMatch.GetProtoShowerW()));
962 
963  this->BuildAmbiguousSpines(pAlgorithm, hitType, protoShower, nuVertex2D, ambiguousHitSpines, hitsToExcludeInEnergyCalcs);
964 
965  if (ambiguousHitSpines.empty())
966  return false;
967 
968  const CartesianVector &startDirection(protoShower.GetConnectionPathway().GetStartDirection());
969  float startL(-std::numeric_limits<float>::max());
970  float ambiguousHitEnergyMean(0.f);
971 
972  // Get average energy of the ambiguous hits
973  for (const CaloHit *const pAmbiguousCaloHit : protoShower.GetAmbiguousHitList())
974  {
975  const float thisT(startDirection.GetCrossProduct(pAmbiguousCaloHit->GetPositionVector() - nuVertex2D).GetMagnitude());
976  const float thisL(startDirection.GetDotProduct(pAmbiguousCaloHit->GetPositionVector() - nuVertex2D));
977 
978  if ((thisL > startL) && (thisT < m_maxTransverseDistance))
979  startL = thisL;
980 
981  ambiguousHitEnergyMean += pAmbiguousCaloHit->GetElectromagneticEnergy() * 1000.f;
982  }
983 
984  ambiguousHitEnergyMean /= protoShower.GetAmbiguousHitList().size();
985 
986  // Get mean energy of other pathways, avoiding the float counting hits
987  float otherEnergyMeanSum(0.f);
988 
989  for (const auto &entry : ambiguousHitSpines)
990  {
991  int nOtherEnergyHits(0);
992  float otherEnergyMean(0.f);
993 
994  for (const CaloHit *const pOtherCaloHit : entry.second)
995  {
996  if (std::find(hitsToExcludeInEnergyCalcs.begin(), hitsToExcludeInEnergyCalcs.end(), pOtherCaloHit) != hitsToExcludeInEnergyCalcs.end())
997  continue;
998 
999  otherEnergyMean += pOtherCaloHit->GetElectromagneticEnergy() * 1000.f;
1000  ++nOtherEnergyHits;
1001  }
1002 
1003  if (nOtherEnergyHits == 0)
1004  continue;
1005 
1006  otherEnergyMean /= static_cast<float>(nOtherEnergyHits);
1007  otherEnergyMeanSum += otherEnergyMean;
1008  }
1009 
1010  // Get the spine mean energy, only consider a limited number of non-ambiguous hits
1011  float spineEnergyMean(0.f);
1012  unsigned int nSpineEnergyHits(0);
1013 
1014  for (const CaloHit *const pSpineCaloHit : protoShower.GetSpineHitList())
1015  {
1016  if (std::find(protoShower.GetAmbiguousHitList().begin(), protoShower.GetAmbiguousHitList().end(), pSpineCaloHit) !=
1017  protoShower.GetAmbiguousHitList().end())
1018  continue;
1019 
1020  const float thisL(startDirection.GetDotProduct(pSpineCaloHit->GetPositionVector() - nuVertex2D));
1021 
1022  if ((thisL > startL) && (nSpineEnergyHits < m_maxSampleHits))
1023  {
1024  spineEnergyMean += pSpineCaloHit->GetElectromagneticEnergy() * 1000.f;
1025  ++nSpineEnergyHits;
1026  }
1027  }
1028 
1029  if (nSpineEnergyHits == 0)
1030  return false;
1031 
1032  spineEnergyMean /= static_cast<float>(nSpineEnergyHits);
1033  unaccountedHitEnergy = ambiguousHitEnergyMean - otherEnergyMeanSum - spineEnergyMean;
1034 
1035  return true;
1036 }
1037 
1038 //------------------------------------------------------------------------------------------------------------------------------------------
1039 
1040 void AmbiguousRegionFeatureTool::BuildAmbiguousSpines(const Algorithm *const pAlgorithm, const HitType hitType, const ProtoShower &protoShower,
1041  const CartesianVector &nuVertex2D, std::map<int, CaloHitList> &ambiguousHitSpines, CaloHitList &hitsToExcludeInEnergyCalcs)
1042 {
1043  const CaloHitList *pCaloHitList;
1044 
1045  if (this->GetHitListOfType(pAlgorithm, hitType, pCaloHitList) != STATUS_CODE_SUCCESS)
1046  return;
1047 
1048  std::map<int, CaloHitList> ambiguousHitSpinesTemp;
1049 
1050  for (const CaloHit *const pCaloHit : *pCaloHitList)
1051  {
1052  if (std::find(protoShower.GetAmbiguousHitList().begin(), protoShower.GetAmbiguousHitList().end(), pCaloHit) !=
1053  protoShower.GetAmbiguousHitList().end())
1054  continue;
1055 
1056  int count(0);
1057 
1058  // A hit can be in more than one spine
1059  for (unsigned int i = 0; i < protoShower.GetAmbiguousDirectionVector().size(); ++i)
1060  {
1061  const CartesianVector &significantDirection(protoShower.GetAmbiguousDirectionVector()[i]);
1062  const CartesianVector displacement(pCaloHit->GetPositionVector() - nuVertex2D);
1063  const float thisT(significantDirection.GetCrossProduct(displacement).GetMagnitude());
1064  const float thisL(significantDirection.GetDotProduct(displacement));
1065 
1066  if ((thisL > 0.f) && (thisT < m_maxTransverseDistance))
1067  {
1068  ++count;
1069  ambiguousHitSpinesTemp[i].push_back(pCaloHit);
1070  }
1071 
1072  if (count == 2)
1073  hitsToExcludeInEnergyCalcs.push_back(pCaloHit);
1074  }
1075  }
1076 
1077  if (ambiguousHitSpinesTemp.empty())
1078  return;
1079 
1080  // Make sure the pathways are continuous
1081  for (const auto &entry : ambiguousHitSpinesTemp)
1082  {
1083  CaloHitList continuousSpine(this->FindAmbiguousContinuousSpine(entry.second, protoShower.GetAmbiguousHitList(), nuVertex2D));
1084 
1085  if (continuousSpine.size() > 0)
1086  ambiguousHitSpines[entry.first] = continuousSpine;
1087  }
1088 }
1089 
1090 //------------------------------------------------------------------------------------------------------------------------------------------
1091 
1092 StatusCode AmbiguousRegionFeatureTool::GetHitListOfType(const Algorithm *const pAlgorithm, const HitType hitType, const CaloHitList *&pCaloHitList) const
1093 {
1094  const std::string typeHitListName(hitType == TPC_VIEW_U ? m_caloHitListNameU
1095  : hitType == TPC_VIEW_V ? m_caloHitListNameV
1096  : m_caloHitListNameW);
1097 
1098  PANDORA_THROW_RESULT_IF_AND_IF(
1099  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*pAlgorithm, typeHitListName, pCaloHitList));
1100 
1101  if (!pCaloHitList || pCaloHitList->empty())
1102  return STATUS_CODE_NOT_INITIALIZED;
1103 
1104  return STATUS_CODE_SUCCESS;
1105 }
1106 
1107 //------------------------------------------------------------------------------------------------------------------------------------------
1108 
1110  const CaloHitList &caloHitList, const CaloHitList &ambiguousHitList, const CartesianVector &nuVertex2D)
1111 {
1112  CaloHitList continuousHitList;
1113 
1114  CaloHitVector caloHitVector(caloHitList.begin(), caloHitList.end());
1115  std::sort(caloHitVector.begin(), caloHitVector.end(), LArConnectionPathwayHelper::SortByDistanceToPoint(nuVertex2D));
1116 
1117  for (unsigned int i = 0; i < caloHitVector.size(); ++i)
1118  {
1119  CaloHitList connectedHitList;
1120  connectedHitList.push_back(caloHitVector[i]);
1121 
1122  if (LArClusterHelper::GetClosestDistance(connectedHitList.front()->GetPositionVector(), ambiguousHitList) > m_maxHitSeparation)
1123  continue;
1124 
1125  bool found(true);
1126 
1127  while (found)
1128  {
1129  found = false;
1130 
1131  for (unsigned int j = (i + 1); j < caloHitVector.size(); ++j)
1132  {
1133  const CaloHit *const pCaloHit(caloHitVector[j]);
1134 
1135  if (std::find(connectedHitList.begin(), connectedHitList.end(), pCaloHit) != connectedHitList.end())
1136  continue;
1137 
1138  if (LArClusterHelper::GetClosestDistance(pCaloHit->GetPositionVector(), connectedHitList) < m_maxHitSeparation)
1139  {
1140  // to avoid ends of tracks
1141  if (static_cast<float>(connectedHitList.size()) < static_cast<float>(caloHitVector.size() * m_maxTrackFraction))
1142  {
1143  found = true;
1144  connectedHitList.push_back(pCaloHit);
1145  }
1146 
1147  break;
1148  }
1149  }
1150  }
1151 
1152  if (connectedHitList.size() >= 2)
1153  {
1154  continuousHitList.insert(continuousHitList.begin(), connectedHitList.begin(), connectedHitList.end());
1155  break;
1156  }
1157  }
1158 
1159  return continuousHitList;
1160 }
1161 
1162 //------------------------------------------------------------------------------------------------------------------------------------------
1163 
1164 StatusCode AmbiguousRegionFeatureTool::ReadSettings(const TiXmlHandle xmlHandle)
1165 {
1166  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DefaultFloat", m_defaultFloat));
1167 
1168  PANDORA_RETURN_RESULT_IF_AND_IF(
1169  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CaloHitListNameU", m_caloHitListNameU));
1170 
1171  PANDORA_RETURN_RESULT_IF_AND_IF(
1172  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CaloHitListNameV", m_caloHitListNameV));
1173 
1174  PANDORA_RETURN_RESULT_IF_AND_IF(
1175  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CaloHitListNameW", m_caloHitListNameW));
1176 
1177  PANDORA_RETURN_RESULT_IF_AND_IF(
1178  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxTransverseDistance", m_maxTransverseDistance));
1179 
1180  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxSampleHits", m_maxSampleHits));
1181 
1182  PANDORA_RETURN_RESULT_IF_AND_IF(
1183  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxHitSeparation", m_maxHitSeparation));
1184 
1185  PANDORA_RETURN_RESULT_IF_AND_IF(
1186  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxTrackFraction", m_maxTrackFraction));
1187 
1188  return STATUS_CODE_SUCCESS;
1189 }
1190 
1191 } // namespace lar_content
float m_maxTransverseDistance
The max. proximity of a hits, included in a trajectory energy calcs.
Header file for the connection pathway feature tools.
const ProtoShower & GetProtoShowerU() const
Get the U view ProtoShower.
void CalculateNAmbiguousViews(const ProtoShowerMatch &protoShowerMatch, float &nAmbiguousViews)
Count the number of views with ambiguous hits.
float m_edgeStep
The binning of the shower boundaries.
std::string m_caloHitListNameW
The event W view hit list.
const ProtoShower & GetProtoShowerV() const
Get the V view ProtoShower.
void Run(LArMvaHelper::MvaFeatureVector &featureVector, const pandora::Algorithm *const pAlgorithm, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::CartesianPointVector &showerStarts3D)
float m_maxFoundHitRatioLimit
maxFoundHitRatio max. limit
void CalculateViewShowerStartConsistencyVariables(const TwoDSlidingFitResult &showerFitResult, const pandora::CaloHitList &postShowerHitList, const bool isShowerDownstream, float &showerStartEnergyAsymmetry, float &showerStartMoliereRadius)
Evaluate the shower start consistency variables.
Header file for the pfo helper class.
Header file for the ProtoShower class.
float m_maxInitialGapSizeLimit
maxInitialGapSizeLimit max. limit
const pandora::CartesianPointVector & GetAmbiguousDirectionVector() const
Get the ambiguous direction vector.
Header file for the connection pathway helper class.
bool GetViewAmbiguousHitVariables(const pandora::Algorithm *const pAlgorithm, const ProtoShowerMatch &protoShowerMatch, const pandora::HitType hitType, const pandora::CartesianVector &nuVertex3D, float &unaccountedHitEnergy)
Calculate the ambiguous region variables for the input view.
MvaTypes::MvaFeatureVector MvaFeatureVector
Definition: LArMvaHelper.h:75
const pandora::CaloHitList & GetAmbiguousHitList() const
Get the ambiguous hit list.
void CalculateViewNuVertexConsistencyVariables(const TwoDSlidingFitResult &spineFitResult, const pandora::CaloHitList &postShowerHitList, const bool isDownstream, const pandora::CartesianVector &nuVertex2D, float &nuVertexEnergyAsymmetry, float &nuVertexEnergyWeightedMeanRadialDistance)
Evaluate the neutrino vertex consistency variables.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_maxNuVertexEnergyWeightedMeanRadialDistanceLimit
maxNuVertexEnergyWeightedMeanRadialDistance max. limit
void Run(LArMvaHelper::MvaFeatureVector &featureVector, const pandora::Algorithm *const pAlgorithm, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::CartesianPointVector &showerStarts3D)
void GetViewShowerRegionVariables(const pandora::Algorithm *const pAlgorithm, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::HitType hitType, const pandora::CartesianVector &showerStart3D, float &nHits, float &foundHitRatio, float &scatterAngle, float &openingAngle, float &nuVertexEnergyAsymmetry, float &nuVertexEnergyWeightedMeanRadialDistance, float &showerStartEnergyAsymmetry, float &showerStartMoliereRadius)
Calculate the shower region variables for the input view.
float m_maxOpeningAngleLimit
maxOpeningAngle max. limit
const ConnectionPathway & GetConnectionPathway() const
Get the connection pathway.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
unsigned int m_spineFitWindow
The spine fit window.
float m_maxScatterAngleLimit
maxScatterAngle max. limit
intermediate_table::const_iterator const_iterator
unsigned int m_showerFitWindow
The shower fit window.
MvaTypes::MvaFeatureMap MvaFeatureMap
Definition: LArMvaHelper.h:78
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
const pandora::CaloHitList & GetSpineHitList() const
Get the spine hit list.
unsigned int m_spineFitWindow
The spine fit window.
ProtoShowerMatch class.
TFile f
Definition: plotHisto.C:6
std::string m_caloHitListNameU
The event U view hit list.
Header file for the geometry helper class.
float m_defaultRatio
Default float value for ratios.
float m_minLargestGapSizeLimit
minLargestGapSizeLimit max. limit
pandora::StatusCode GetHitListOfType(const pandora::Algorithm *const pAlgorithm, const pandora::HitType hitType, const pandora::CaloHitList *&pCaloHitList) const
Obtain the event hit list of a given view.
float m_showerRadius
The max. separation distance between a shower region hit and the shower core.
void GetShowerHitVariables(const pandora::CaloHitList &spineHitList, const pandora::CaloHitList &postShowerHitList, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::HitType hitType, float &nHits, float &foundHitRatio)
Evaluate the variables associated with the shower region hit multiplicity.
Header file for the cluster helper class.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
unsigned int m_nHitsToConsider
The number of hits which defines the initial region.
float m_pathwayScatteringAngle2DLimit
pathwayScatteringAngle2DLimit max. limit
void GetViewInitialRegionVariables(const pandora::Algorithm *const pAlgorithm, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::HitType hitType, float &initialGapSize, float &largestGapSize) const
Calculate the initial region variables for the input view.
float m_moliereFraction
The energy fraction which corresponds to minShowerStartMoliereRadius.
std::map< int, LayerFitResult > LayerFitResultMap
void Run(LArMvaHelper::MvaFeatureVector &featureVector, const pandora::Algorithm *const pAlgorithm, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::CartesianPointVector &showerStarts3D)
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float GetLayerPitch() const
Get the layer pitch, units cm.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
unsigned int m_maxSampleHits
The max. number of hits considered in the spine energy calcs.
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.
void BuildAmbiguousSpines(const pandora::Algorithm *const pAlgorithm, const pandora::HitType hitType, const ProtoShower &protoShower, const pandora::CartesianVector &nuVertex2D, std::map< int, pandora::CaloHitList > &ambiguousHitSpines, pandora::CaloHitList &hitsToExcludeInEnergyCalcs)
Determine the spine hits of the particles with which the ambiguous hits are shared.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
ProtoShower class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_pathwayLengthLimit
pathwayLengthLimit max. limit
HitType
Definition: HitType.h:12
pandora::CaloHitList FindAmbiguousContinuousSpine(const pandora::CaloHitList &caloHitList, const pandora::CaloHitList &ambiguousHitList, const pandora::CartesianVector &nuVertex2D)
Determine a continuous pathway of an ambigous particle&#39;s spine hits.
float m_maxHitSeparation
The max. separation of connected hits.
static void GetMinMiddleMax(const float value1, const float value2, const float value3, float &minValue, float &middleValue, float &maxValue)
Determine the lowest, median and highest value from an input of three numbers.
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.
std::string m_caloHitListNameV
The event V view hit list.
void CalculateViewOpeningAngle(const TwoDSlidingFitResult &showerFitResult, const pandora::CaloHitList &postShowerHitList, const pandora::CartesianVector &showerStart2D, float &openingAngle)
Calculate the opening angle of the shower region.
float m_maxTrackFraction
The fraction of found hits which are considered in the energy calcs.
void BuildViewShower(const pandora::ParticleFlowObject *const pShowerPfo, const TwoDSlidingFitResult &spineFit, const pandora::HitType hitType, const pandora::CartesianVector &showerStart2D, const pandora::CartesianVector &nuVertex2D, pandora::CaloHitList &postShowerHitList, pandora::CartesianPointVector &postShowerPositions)
Collect the shower region hits in a given view.
const ProtoShower & GetProtoShowerW() const
Get the W view ProtoShower.
static void GetCaloHits(const pandora::PfoList &pfoList, const pandora::HitType &hitType, pandora::CaloHitList &caloHitList)
Get a list of calo hits of a particular hit type from a list of pfos.
void Run(LArMvaHelper::MvaFeatureVector &featureVector, const pandora::Algorithm *const pAlgorithm, const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &nuVertex3D, const ProtoShowerMatch &protoShowerMatch, const pandora::CartesianPointVector &showerStarts3D)
int m_nLayersHalfWindow
The half window of each segment.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
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.
float Get2DKink(const pandora::Algorithm *const pAlgorithm, const ProtoShowerMatch &protoShowerMatch, const pandora::CartesianVector &showerStart3D) const
Obtain a cautious estimate of the largest 2D deflection of the connection pathway.
float GetLargest2DKinkFromView(const pandora::Algorithm *const pAlgorithm, const TwoDSlidingFitResult &spineFit, const pandora::HitType hitType, const pandora::CartesianVector &showerStart3D) const
Obtain a cautious estimate of the largest 2D deflection of a connection pathway in a given view...
void CalculateViewScatterAngle(const pandora::CartesianVector &nuVertex2D, const TwoDSlidingFitResult &spineFitResult, const pandora::CartesianVector &showerStart2D, const TwoDSlidingFitResult &showerFitResult, float &scatterAngle)
Calculate the connection pathway-shower region scatter angle.
const pandora::CartesianVector & GetStartDirection() const
Get the start direction of the connection pathway.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.
float m_minShowerStartMoliereRadiusLimit
minShowerStartMoliereRadius max. limit