LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
ShowerStartFinderTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 #include "Pandora/AlgorithmTool.h"
11 
14 
18 
20 
21 using namespace pandora;
22 
23 namespace lar_content
24 {
25 
26 ShowerStartFinderTool::ShowerStartFinderTool() :
27  m_spineSlidingFitWindow(10),
28  m_longitudinalCoordinateBinSize(1.f),
29  m_nInitialEnergyBins(5),
30  m_minSigmaDeviation(1.f),
31  m_maxEdgeGap(3.f),
32  m_longitudinalShowerFraction(0.85f),
33  m_minShowerOpeningAngle(2.f),
34  m_molliereRadius(4.5f),
35  m_showerSlidingFitWindow(1000),
36  m_maxLayerSeparation(5)
37 {
38 }
39 
40 //------------------------------------------------------------------------------------------------------------------------------------------
41 
42 StatusCode ShowerStartFinderTool::Run(const ParticleFlowObject *const pShowerPfo, const CartesianVector &peakDirection,
43  const HitType hitType, const CaloHitList &showerSpineHitList, CartesianVector &showerStartPosition, CartesianVector &showerStartDirection)
44 {
45  try
46  {
47  // Fit the shower spine
48  CartesianPointVector hitPositions;
49  for (const CaloHit *const pCaloHit : showerSpineHitList)
50  hitPositions.push_back(pCaloHit->GetPositionVector());
51 
52  const TwoDSlidingFitResult spineTwoDSlidingFit(
53  &hitPositions, m_spineSlidingFitWindow, LArGeometryHelper::GetWirePitch(this->GetPandora(), hitType));
54 
55  // First obtain the longitudinal position of spine hits
56  LongitudinalPositionMap longitudinalPositionMap;
57  this->ObtainLongitudinalDecomposition(spineTwoDSlidingFit, showerSpineHitList, longitudinalPositionMap);
58 
59  // Obtain spine energy profile
60  EnergySpectrumMap energySpectrumMap;
61  this->GetEnergyDistribution(showerSpineHitList, longitudinalPositionMap, energySpectrumMap);
62 
63  // Now find the shower start position/direction
64  const bool isEndDownstream(peakDirection.GetZ() > 0.f);
65 
66  this->FindShowerStartAndDirection(pShowerPfo, hitType, spineTwoDSlidingFit, energySpectrumMap, showerSpineHitList, isEndDownstream,
67  showerStartPosition, showerStartDirection);
68  }
69  catch (...)
70  {
71  return STATUS_CODE_FAILURE;
72  }
73 
74  return STATUS_CODE_SUCCESS;
75 }
76 
77 //------------------------------------------------------------------------------------------------------------------------------------------
78 
80  const CaloHitList &showerSpineHitList, LongitudinalPositionMap &longitudinalPositionMap) const
81 {
82  // Find hits in each layer
83  LayerToHitMap layerToHitMap;
84 
85  for (const CaloHit *const pCaloHit : showerSpineHitList)
86  {
87  float hitL(0.f), hitT(0.f);
88  const CartesianVector hitPosition(pCaloHit->GetPositionVector());
89 
90  spineTwoDSlidingFit.GetLocalPosition(hitPosition, hitL, hitT);
91  layerToHitMap[spineTwoDSlidingFit.GetLayer(hitL)].push_back(pCaloHit);
92  }
93 
94  // Walk through layers and store the longitudinal distance of each hit
95  // IMPORTANT: layer positions correspond to the middle of the layer
96  // IMPORTANT: the longitudinal distance is measured relative to the first layer position (where l = 0)
97  // IMPORTANT: at any point, the running distance is that at the lowest neighbour central position
98  float runningDistance(0);
99  const LayerFitResultMap &layerFitResultMap(spineTwoDSlidingFit.GetLayerFitResultMap());
100 
101  for (auto iter = layerToHitMap.begin(); iter != layerToHitMap.end(); ++iter)
102  {
103  const int layer(iter->first);
104  const float layerL(layerFitResultMap.at(layer).GetL());
105 
106  // Get global positions of layer and its neighbours
107  const int higherLayer(std::next(iter) == layerToHitMap.end() ? layer : std::next(iter)->first);
108  const int middleLayer(layer);
109  const int lowerLayer(iter == layerToHitMap.begin() ? layer : std::prev(iter)->first);
110  CartesianVector lowerLayerPosition(0.f, 0.f, 0.f), middleLayerPosition(0.f, 0.f, 0.f), higherLayerPosition(0.f, 0.f, 0.f);
111 
112  spineTwoDSlidingFit.GetGlobalPosition(layerFitResultMap.at(lowerLayer).GetL(), layerFitResultMap.at(lowerLayer).GetFitT(), lowerLayerPosition);
113  spineTwoDSlidingFit.GetGlobalPosition(layerFitResultMap.at(middleLayer).GetL(), layerFitResultMap.at(middleLayer).GetFitT(), middleLayerPosition);
114  spineTwoDSlidingFit.GetGlobalPosition(layerFitResultMap.at(higherLayer).GetL(), layerFitResultMap.at(higherLayer).GetFitT(), higherLayerPosition);
115 
116  const float layerLength = std::next(iter) == layerToHitMap.end() ? 0.f
117  : iter == layerToHitMap.begin() ? 0.f
118  : (middleLayerPosition - lowerLayerPosition).GetMagnitude();
119 
120  for (const CaloHit *const pCaloHit : iter->second)
121  {
122  const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
123 
124  float hitL(0.f), hitT(0.f);
125  spineTwoDSlidingFit.GetLocalPosition(hitPosition, hitL, hitT);
126 
127  // Is the hit below or above the middle of the layer?
128  // Therefore, determine which layer positions the hit is between
129  CartesianVector localLowerLayerPosition(0.f, 0.f, 0.f), localHigherLayerPosition(0.f, 0.f, 0.f);
130 
131  if (hitL < layerL)
132  {
133  localLowerLayerPosition = lowerLayerPosition;
134  localHigherLayerPosition = iter == layerToHitMap.begin() ? higherLayerPosition : middleLayerPosition;
135  }
136  else
137  {
138  localLowerLayerPosition = std::next(iter) == layerToHitMap.end() ? lowerLayerPosition : middleLayerPosition;
139  localHigherLayerPosition = higherLayerPosition;
140  }
141 
142  // Determine the axis to project onto
143  // Calculate the distance of axis between its two nearest neighbours
144  const CartesianVector displacement((higherLayerPosition - lowerLayerPosition).GetUnitVector());
145  float longitudinalDisplacement = (displacement.GetDotProduct(hitPosition - lowerLayerPosition) + runningDistance);
146 
147  // If we've passed the central position, we need to add on the layer length
148  longitudinalDisplacement += (hitL > layerL ? layerLength : 0.f);
149 
150  longitudinalPositionMap[pCaloHit] = longitudinalDisplacement;
151  }
152 
153  runningDistance += layerLength;
154  }
155 }
156 
157 //------------------------------------------------------------------------------------------------------------------------------------------
158 
159 void ShowerStartFinderTool::GetEnergyDistribution(const CaloHitList &showerSpineHitList,
160  const LongitudinalPositionMap &longitudinalPositionMap, EnergySpectrumMap &energySpectrumMap) const
161 {
162  float e0(0.f);
163 
164  for (const CaloHit *pCaloHit : showerSpineHitList)
165  e0 += pCaloHit->GetElectromagneticEnergy();
166 
167  for (const CaloHit *pCaloHit : showerSpineHitList)
168  {
169  const float fractionalEnergy(pCaloHit->GetElectromagneticEnergy() / e0);
170  const float projection(longitudinalPositionMap.at(pCaloHit));
171  const int longitudinalIndex = std::floor(projection / m_longitudinalCoordinateBinSize);
172 
173  if (energySpectrumMap.find(longitudinalIndex) == energySpectrumMap.end())
174  energySpectrumMap[longitudinalIndex] = fractionalEnergy;
175  else
176  energySpectrumMap[longitudinalIndex] += fractionalEnergy;
177  }
178 }
179 
180 //------------------------------------------------------------------------------------------------------------------------------------------
181 
182 void ShowerStartFinderTool::FindShowerStartAndDirection(const ParticleFlowObject *const pShowerPfo, const HitType hitType,
183  const TwoDSlidingFitResult &spineTwoDSlidingFit, const EnergySpectrumMap &energySpectrumMap, const CaloHitList &showerSpineHitList,
184  const bool isEndDownstream, CartesianVector &showerStartPosition, CartesianVector &showerStartDirection) const
185 {
186  // Search for shower start at significant energy deviations
187  int longitudinalStartBin(0);
188 
189  if (isEndDownstream)
190  {
191  longitudinalStartBin = this->FindShowerStartLongitudinalCoordinate(pShowerPfo, hitType, spineTwoDSlidingFit, energySpectrumMap,
192  showerSpineHitList, isEndDownstream, energySpectrumMap.begin(), std::prev(energySpectrumMap.end()));
193  }
194  else
195  {
196  longitudinalStartBin = this->FindShowerStartLongitudinalCoordinate(pShowerPfo, hitType, spineTwoDSlidingFit, energySpectrumMap,
197  showerSpineHitList, isEndDownstream, energySpectrumMap.rbegin(), std::prev(energySpectrumMap.rend()));
198  }
199 
200  const float longitudinalStartCoordinate(longitudinalStartBin * m_longitudinalCoordinateBinSize);
201 
202  this->ConvertLongitudinalProjectionToGlobal(spineTwoDSlidingFit, longitudinalStartCoordinate, showerStartPosition, showerStartDirection);
203 
204  showerStartDirection = isEndDownstream ? showerStartDirection : showerStartDirection * (-1.f);
205 }
206 
207 //------------------------------------------------------------------------------------------------------------------------------------------
208 
209 template <typename T>
210 int ShowerStartFinderTool::FindShowerStartLongitudinalCoordinate(const ParticleFlowObject *const pShowerPfo, const HitType hitType,
211  const TwoDSlidingFitResult &spineTwoDSlidingFit, const EnergySpectrumMap &energySpectrumMap, const CaloHitList &showerSpineHitList,
212  const bool isEndDownstream, const T startIter, const T endIter) const
213 {
214  // Characterise initial energy
215  float meanEnergy(0.f), energySigma(0.f);
216 
217  if (energySpectrumMap.size() > m_nInitialEnergyBins)
218  this->CharacteriseInitialEnergy(energySpectrumMap, isEndDownstream, meanEnergy, energySigma);
219 
220  // Now loop through energySpectrumMap
221  T iter(startIter);
222 
223  std::advance(iter, (energySpectrumMap.size() > m_nInitialEnergyBins) ? m_nInitialEnergyBins : (energySpectrumMap.size() - 1));
224 
225  for (; iter != endIter; iter++)
226  {
227  const float longitudinalCoordinate(iter->first * m_longitudinalCoordinateBinSize);
228  const float energyDeviation((iter->second - meanEnergy) / energySigma);
229 
230  try
231  {
232  // Use energy and local topology to assess whether we are at the shower start
233  if ((energyDeviation > m_minSigmaDeviation) &&
234  this->IsShowerTopology(pShowerPfo, hitType, spineTwoDSlidingFit, longitudinalCoordinate, showerSpineHitList, isEndDownstream))
235  {
236  break;
237  }
238  }
239  catch (...)
240  {
241  continue;
242  }
243  }
244 
245  return iter->first;
246 }
247 
248 //------------------------------------------------------------------------------------------------------------------------------------------
249 
251  const EnergySpectrumMap &energySpectrumMap, const bool isEndDownstream, float &meanEnergy, float &energySigma) const
252 {
253  auto energySpectrumIter(isEndDownstream ? energySpectrumMap.begin() : std::prev(energySpectrumMap.end()));
254 
255  for (unsigned int i = 0; i < m_nInitialEnergyBins; ++i)
256  {
257  meanEnergy += energySpectrumIter->second;
258  isEndDownstream ? ++energySpectrumIter : --energySpectrumIter;
259  }
260 
261  meanEnergy /= static_cast<float>(m_nInitialEnergyBins);
262 
263  energySpectrumIter = isEndDownstream ? energySpectrumMap.begin() : std::prev(energySpectrumMap.end());
264 
265  for (unsigned int i = 0; i < m_nInitialEnergyBins; ++i)
266  {
267  energySigma += std::pow(energySpectrumIter->second - meanEnergy, 2);
268  isEndDownstream ? ++energySpectrumIter : --energySpectrumIter;
269  }
270 
271  energySigma = std::sqrt(energySigma / m_nInitialEnergyBins);
272 }
273 
274 //------------------------------------------------------------------------------------------------------------------------------------------
275 
276 bool ShowerStartFinderTool::IsShowerTopology(const ParticleFlowObject *const pShowerPfo, const HitType hitType, const TwoDSlidingFitResult &spineTwoDSlidingFit,
277  const float longitudinalDistance, const CaloHitList &showerSpineHitList, const bool isEndDownstream) const
278 {
279  CartesianVector showerStartPosition(0.f, 0.f, 0.f), showerStartDirection(0.f, 0.f, 0.f);
280  this->ConvertLongitudinalProjectionToGlobal(spineTwoDSlidingFit, longitudinalDistance, showerStartPosition, showerStartDirection);
281 
282  // Build shower
283  CartesianPointVector showerRegionPositionVector;
284  if (this->BuildShowerRegion(pShowerPfo, hitType, showerSpineHitList, showerStartPosition, showerStartDirection, isEndDownstream,
285  showerRegionPositionVector) != STATUS_CODE_SUCCESS)
286  {
287  return false;
288  }
289 
290  // Characterise the shower
291  bool isBetween(false);
292  CartesianVector positiveEdgeStart(0.f, 0.f, 0.f), positiveEdgeEnd(0.f, 0.f, 0.f), positiveEdgeDirection(0.f, 0.f, 0.f);
293  CartesianVector negativeEdgeStart(0.f, 0.f, 0.f), negativeEdgeEnd(0.f, 0.f, 0.f), negativeEdgeDirection(0.f, 0.f, 0.f);
294 
295  if (this->CharacteriseShowerTopology(showerRegionPositionVector, showerStartPosition, hitType, isEndDownstream, showerStartDirection,
296  positiveEdgeStart, positiveEdgeEnd, negativeEdgeStart, negativeEdgeEnd, isBetween) != STATUS_CODE_SUCCESS)
297  {
298  return false;
299  }
300 
301  if (!isBetween)
302  return false;
303 
304  positiveEdgeStart = showerStartPosition;
305  negativeEdgeStart = showerStartPosition;
306  positiveEdgeDirection = positiveEdgeEnd - positiveEdgeStart;
307  negativeEdgeDirection = negativeEdgeEnd - negativeEdgeStart;
308 
309  const float showerOpeningAngle(positiveEdgeDirection.GetOpeningAngle(negativeEdgeDirection) * 180.f / 3.14);
310 
311  if (showerOpeningAngle < m_minShowerOpeningAngle)
312  return false;
313 
314  return true;
315 }
316 
317 //------------------------------------------------------------------------------------------------------------------------------------------
318 
320  const float longitudinalDistance, CartesianVector &globalPosition, CartesianVector &globalDirection) const
321 {
322  const LayerFitResultMap &layerFitResultMap(spineTwoDSlidingFit.GetLayerFitResultMap());
323 
324  float runningDistance(0.f);
325  int showerStartLayer(0);
326  CartesianVector previousLayerPosition(0.f, 0.f, 0.f);
327  spineTwoDSlidingFit.GetGlobalPosition(layerFitResultMap.begin()->second.GetL(), layerFitResultMap.begin()->second.GetFitT(), previousLayerPosition);
328 
329  for (auto iter = std::next(layerFitResultMap.begin()); iter != layerFitResultMap.end(); ++iter)
330  {
331  const int layer(iter->first);
332  showerStartLayer = layer;
333 
334  CartesianVector layerPosition(0.f, 0.f, 0.f);
335  spineTwoDSlidingFit.GetGlobalPosition(layerFitResultMap.at(layer).GetL(), layerFitResultMap.at(layer).GetFitT(), layerPosition);
336 
337  runningDistance += (layerPosition - previousLayerPosition).GetMagnitude();
338 
339  if (runningDistance > longitudinalDistance)
340  break;
341 
342  previousLayerPosition = layerPosition;
343  }
344 
345  const float lCoordinate(layerFitResultMap.at(showerStartLayer).GetL()), tCoordinate(layerFitResultMap.at(showerStartLayer).GetFitT());
346  const float localGradient(layerFitResultMap.at(showerStartLayer).GetGradient());
347 
348  spineTwoDSlidingFit.GetGlobalPosition(lCoordinate, tCoordinate, globalPosition);
349  spineTwoDSlidingFit.GetGlobalDirection(localGradient, globalDirection);
350 }
351 
352 //------------------------------------------------------------------------------------------------------------------------------------------
353 
354 StatusCode ShowerStartFinderTool::BuildShowerRegion(const ParticleFlowObject *const pShowerPfo, const HitType hitType,
355  const CaloHitList &showerSpineHitList, const CartesianVector &showerStartPosition, const CartesianVector &showerStartDirection,
356  const bool isEndDownstream, CartesianPointVector &showerRegionPositionVector) const
357 {
358  CaloHitList viewShowerHitList;
359  LArPfoHelper::GetCaloHits(pShowerPfo, hitType, viewShowerHitList);
360 
361  // Find shower region hits
362  CaloHitList showerRegionHitList;
363 
364  for (const CaloHit *const pCaloHit : viewShowerHitList)
365  {
366  if (std::find(showerSpineHitList.begin(), showerSpineHitList.end(), pCaloHit) != showerSpineHitList.end())
367  continue;
368 
369  const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
370  float l(showerStartDirection.GetDotProduct(hitPosition - showerStartPosition));
371 
372  l *= isEndDownstream ? 1.f : -1.f;
373 
374  if (l < 0.f)
375  continue;
376 
377  if (showerStartDirection.GetCrossProduct(hitPosition - showerStartPosition).GetMagnitude() > m_molliereRadius)
378  continue;
379 
380  showerRegionPositionVector.push_back(pCaloHit->GetPositionVector());
381  showerRegionHitList.push_back(pCaloHit);
382  }
383 
384  // Searching for a continuous shower, so truncate at first gap
385  CartesianPointVector coordinateListP, coordinateListN;
386  try
387  {
388  // Perform shower sliding fit
389  const TwoDSlidingShowerFitResult twoDShowerSlidingFit(
390  &showerRegionPositionVector, m_showerSlidingFitWindow, LArGeometryHelper::GetWirePitch(this->GetPandora(), hitType));
391  const LayerFitResultMap &layerFitResultMapS(twoDShowerSlidingFit.GetShowerFitResult().GetLayerFitResultMap());
392  const LayerFitResultMap &layerFitResultMapP(twoDShowerSlidingFit.GetPositiveEdgeFitResult().GetLayerFitResultMap());
393  const LayerFitResultMap &layerFitResultMapN(twoDShowerSlidingFit.GetNegativeEdgeFitResult().GetLayerFitResultMap());
394 
395  float showerStartL(0.f), showerStartT(0.f);
396  twoDShowerSlidingFit.GetShowerFitResult().GetLocalPosition(showerStartPosition, showerStartL, showerStartT);
397 
398  const int startLayer(twoDShowerSlidingFit.GetShowerFitResult().GetLayer(showerStartL));
399 
400  // Check for gaps in the shower
401  for (LayerFitResultMap::const_iterator iterS = layerFitResultMapS.begin(); iterS != layerFitResultMapS.end(); ++iterS)
402  {
403  const int layer(iterS->first);
404 
405  if (isEndDownstream && (layer < startLayer))
406  continue;
407 
408  if (!isEndDownstream && (layer > startLayer))
409  continue;
410 
411  LayerFitResultMap::const_iterator iterP = layerFitResultMapP.find(layer);
412  LayerFitResultMap::const_iterator iterN = layerFitResultMapN.find(layer);
413 
414  if (layerFitResultMapP.end() != iterP)
415  {
416  CartesianVector positiveEdgePosition(0.f, 0.f, 0.f);
417  twoDShowerSlidingFit.GetShowerFitResult().GetGlobalPosition(iterP->second.GetL(), iterP->second.GetFitT(), positiveEdgePosition);
418  coordinateListP.push_back(positiveEdgePosition);
419  }
420 
421  if (layerFitResultMapN.end() != iterN)
422  {
423  CartesianVector negativeEdgePosition(0.f, 0.f, 0.f);
424  twoDShowerSlidingFit.GetShowerFitResult().GetGlobalPosition(iterN->second.GetL(), iterN->second.GetFitT(), negativeEdgePosition);
425  coordinateListN.push_back(negativeEdgePosition);
426  }
427  }
428 
429  if ((coordinateListP.size() == 0) || (coordinateListN.size() == 0))
430  return STATUS_CODE_FAILURE;
431 
432  std::sort(coordinateListP.begin(), coordinateListP.end(), LArConnectionPathwayHelper::SortByDistanceToPoint(showerStartPosition));
433  std::sort(coordinateListN.begin(), coordinateListN.end(), LArConnectionPathwayHelper::SortByDistanceToPoint(showerStartPosition));
434 
435  float pMaximumL(std::numeric_limits<float>::max());
436  CartesianVector previousCoordinateP(coordinateListP.front());
437 
438  for (auto iterP = std::next(coordinateListP.begin()); iterP != coordinateListP.end(); ++iterP)
439  {
440  const CartesianVector &coordinateP(*iterP);
441  const float separationSquared((coordinateP - previousCoordinateP).GetMagnitudeSquared());
442 
443  if (separationSquared > (m_maxEdgeGap * m_maxEdgeGap))
444  break;
445 
446  float thisT(0.f);
447  twoDShowerSlidingFit.GetShowerFitResult().GetLocalPosition(coordinateP, pMaximumL, thisT);
448 
449  previousCoordinateP = coordinateP;
450  }
451 
452  float nMaximumL(std::numeric_limits<float>::max());
453  CartesianVector previousCoordinateN(coordinateListN.front());
454 
455  for (auto iterN = std::next(coordinateListN.begin()); iterN != coordinateListN.end(); ++iterN)
456  {
457  const CartesianVector &coordinateN(*iterN);
458  const float separationSquared((coordinateN - previousCoordinateN).GetMagnitudeSquared());
459 
460  if (separationSquared > (m_maxEdgeGap * m_maxEdgeGap))
461  break;
462 
463  float thisT(0.f);
464  twoDShowerSlidingFit.GetShowerFitResult().GetLocalPosition(coordinateN, nMaximumL, thisT);
465 
466  previousCoordinateN = coordinateN;
467  }
468 
469  // Now truncate the halo hit position vector
470  showerRegionPositionVector.clear();
471 
472  for (const CaloHit *const pCaloHit : showerRegionHitList)
473  {
474  float thisL(0.f), thisT(0.f);
475  const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
476 
477  twoDShowerSlidingFit.GetShowerFitResult().GetLocalPosition(hitPosition, thisL, thisT);
478 
479  // because the end of showers can get really narrow and thus fail the opening angle check
480  if (isEndDownstream && (thisL > (m_longitudinalShowerFraction * std::max(pMaximumL, nMaximumL))))
481  continue;
482 
483  if (!isEndDownstream && (thisL < (m_longitudinalShowerFraction * std::min(pMaximumL, nMaximumL))))
484  continue;
485 
486  showerRegionPositionVector.push_back(pCaloHit->GetPositionVector());
487  }
488  }
489  catch (const StatusCodeException &)
490  {
491  return STATUS_CODE_FAILURE;
492  }
493 
494  return STATUS_CODE_SUCCESS;
495 }
496 
497 //------------------------------------------------------------------------------------------------------------------------------------------
498 
499 StatusCode ShowerStartFinderTool::CharacteriseShowerTopology(const CartesianPointVector &showerRegionPositionVector, const CartesianVector &showerStartPosition,
500  const HitType hitType, const bool isEndDownstream, const CartesianVector &showerStartDirection, CartesianVector &positiveEdgeStart,
501  CartesianVector &positiveEdgeEnd, CartesianVector &negativeEdgeStart, CartesianVector &negativeEdgeEnd, bool &isBetween) const
502 {
503  try
504  {
505  // Perform shower sliding fit
506  const TwoDSlidingShowerFitResult twoDShowerSlidingFit(
507  &showerRegionPositionVector, m_showerSlidingFitWindow, LArGeometryHelper::GetWirePitch(this->GetPandora(), hitType));
508  const LayerFitResultMap &layerFitResultMapS(twoDShowerSlidingFit.GetShowerFitResult().GetLayerFitResultMap());
509  const LayerFitResultMap &layerFitResultMapP(twoDShowerSlidingFit.GetPositiveEdgeFitResult().GetLayerFitResultMap());
510  const LayerFitResultMap &layerFitResultMapN(twoDShowerSlidingFit.GetNegativeEdgeFitResult().GetLayerFitResultMap());
511 
512  float showerStartL(0.f), showerStartT(0.f);
513  twoDShowerSlidingFit.GetShowerFitResult().GetLocalPosition(showerStartPosition, showerStartL, showerStartT);
514 
515  const int startLayer(twoDShowerSlidingFit.GetShowerFitResult().GetLayer(showerStartL));
516  const int endLayer(twoDShowerSlidingFit.GetShowerFitResult().GetMaxLayer());
517 
518  // Does the shower look to originate from the shower start position
519  int layerCount(0);
520  bool isFirstBetween(false), isLastBetween(false);
521  CartesianPointVector coordinateListP, coordinateListN;
522 
523  for (LayerFitResultMap::const_iterator iterS = layerFitResultMapS.begin(); iterS != layerFitResultMapS.end(); ++iterS)
524  {
525  int layer(iterS->first);
526 
527  if (isEndDownstream && (layer < startLayer))
528  continue;
529 
530  if (!isEndDownstream && (layer > startLayer))
531  continue;
532 
533  LayerFitResultMap::const_iterator iterP = layerFitResultMapP.find(layer);
534  LayerFitResultMap::const_iterator iterN = layerFitResultMapN.find(layer);
535 
536  CartesianVector positiveEdgePosition(0.f, 0.f, 0.f), negativeEdgePosition(0.f, 0.f, 0.f);
537  if (layerFitResultMapP.end() != iterP)
538  {
539  twoDShowerSlidingFit.GetShowerFitResult().GetGlobalPosition(iterP->second.GetL(), iterP->second.GetFitT(), positiveEdgePosition);
540  coordinateListP.push_back(positiveEdgePosition);
541  }
542 
543  if (layerFitResultMapN.end() != iterN)
544  {
545  twoDShowerSlidingFit.GetShowerFitResult().GetGlobalPosition(iterN->second.GetL(), iterN->second.GetFitT(), negativeEdgePosition);
546  coordinateListN.push_back(negativeEdgePosition);
547  }
548 
549  if ((layerFitResultMapP.end() != iterP) && (layerFitResultMapN.end() != iterN))
550  {
551  const CartesianVector positiveDisplacement((positiveEdgePosition - showerStartPosition).GetUnitVector());
552  const bool isPositiveClockwise(this->IsClockwiseRotation(showerStartDirection, positiveDisplacement));
553 
554  const CartesianVector negativeDisplacement((negativeEdgePosition - showerStartPosition).GetUnitVector());
555  const bool isNegativeClockwise(this->IsClockwiseRotation(showerStartDirection, negativeDisplacement));
556 
557  ++layerCount;
558 
559  if (layerCount == 1)
560  isFirstBetween = (isPositiveClockwise != isNegativeClockwise);
561 
562  isLastBetween = (isPositiveClockwise != isNegativeClockwise);
563  }
564  }
565 
566  isBetween = (isFirstBetween || isLastBetween);
567 
568  // Find the extremal points of the shower
569  if (this->GetBoundaryExtremalPoints(twoDShowerSlidingFit, layerFitResultMapP, showerStartPosition, startLayer, endLayer,
570  positiveEdgeStart, positiveEdgeEnd) != STATUS_CODE_SUCCESS)
571  {
572  return STATUS_CODE_FAILURE;
573  }
574 
575  if (this->GetBoundaryExtremalPoints(twoDShowerSlidingFit, layerFitResultMapN, showerStartPosition, startLayer, endLayer,
576  negativeEdgeStart, negativeEdgeEnd) != STATUS_CODE_SUCCESS)
577  {
578  return STATUS_CODE_FAILURE;
579  }
580  }
581  catch (const StatusCodeException &)
582  {
583  return STATUS_CODE_FAILURE;
584  }
585 
586  return STATUS_CODE_SUCCESS;
587 }
588 
589 //------------------------------------------------------------------------------------------------------------------------------------------
590 
591 bool ShowerStartFinderTool::IsClockwiseRotation(const CartesianVector &showerStartDirection, const CartesianVector &displacementVector) const
592 {
593  // Clockwise with respect to +ve Z
594  const float openingAngleFromCore(showerStartDirection.GetOpeningAngle(displacementVector));
595 
596  const CartesianVector clockwiseRotation(
597  displacementVector.GetZ() * std::sin(openingAngleFromCore) + displacementVector.GetX() * std::cos(openingAngleFromCore), 0.f,
598  displacementVector.GetZ() * std::cos(openingAngleFromCore) - displacementVector.GetX() * std::sin(openingAngleFromCore));
599 
600  const CartesianVector anticlockwiseRotation(
601  displacementVector.GetZ() * std::sin(-1.f * openingAngleFromCore) + displacementVector.GetX() * std::cos(-1.f * openingAngleFromCore), 0.f,
602  displacementVector.GetZ() * std::cos(-1.f * openingAngleFromCore) - displacementVector.GetX() * std::sin(-1.f * openingAngleFromCore));
603 
604  const float clockwiseT((clockwiseRotation - showerStartDirection).GetMagnitude());
605  const float anticlockwiseT((anticlockwiseRotation - showerStartDirection).GetMagnitude());
606  const bool isClockwise(clockwiseT < anticlockwiseT);
607 
608  return isClockwise;
609 }
610 
611 //------------------------------------------------------------------------------------------------------------------------------------------
612 
614  const LayerFitResultMap &layerFitResultMap, const CartesianVector &showerStartPosition, const int showerStartLayer,
615  const int showerEndLayer, CartesianVector &boundaryEdgeStart, CartesianVector &boundaryEdgeEnd) const
616 {
617  int boundaryStartLayer(std::numeric_limits<int>::max());
618  int boundaryEndLayer(std::numeric_limits<int>::max());
619 
620  for (auto &entry : layerFitResultMap)
621  {
622  const int bestStartSeparation(std::abs(showerStartLayer - boundaryStartLayer));
623  const int thisStartSeparation(std::abs(showerStartLayer - entry.first));
624 
625  if (((bestStartSeparation - thisStartSeparation) > 0) && (thisStartSeparation < bestStartSeparation))
626  boundaryStartLayer = entry.first;
627 
628  const int bestEndSeparation(std::abs(showerEndLayer - boundaryEndLayer));
629  const int thisEndSeparation(std::abs(showerEndLayer - entry.first));
630 
631  if (((bestEndSeparation - thisEndSeparation) > 0) && (thisEndSeparation < bestEndSeparation))
632  boundaryEndLayer = entry.first;
633  }
634 
635  if (std::abs(showerStartLayer - boundaryStartLayer) > m_maxLayerSeparation)
636  return STATUS_CODE_FAILURE;
637 
638  const float showerStartBoundaryLocalL(layerFitResultMap.at(boundaryStartLayer).GetL()),
639  showerStartBoundaryLocalT(layerFitResultMap.at(boundaryStartLayer).GetFitT());
640  const float showerEndBoundaryLocalL(layerFitResultMap.at(boundaryEndLayer).GetL()),
641  showerEndBoundaryLocalT(layerFitResultMap.at(boundaryEndLayer).GetFitT());
642 
643  showerTwoDSlidingFit.GetShowerFitResult().GetGlobalPosition(showerStartBoundaryLocalL, showerStartBoundaryLocalT, boundaryEdgeStart);
644  showerTwoDSlidingFit.GetShowerFitResult().GetGlobalPosition(showerEndBoundaryLocalL, showerEndBoundaryLocalT, boundaryEdgeEnd);
645 
646  // Swap start and end if needs be
647  if ((showerStartPosition - boundaryEdgeEnd).GetMagnitudeSquared() < (showerStartPosition - boundaryEdgeStart).GetMagnitude())
648  {
649  const CartesianVector startTemp(boundaryEdgeStart), endTemp(boundaryEdgeEnd);
650 
651  boundaryEdgeStart = endTemp;
652  boundaryEdgeEnd = startTemp;
653  }
654 
655  return STATUS_CODE_SUCCESS;
656 }
657 
658 //------------------------------------------------------------------------------------------------------------------------------------------
659 
660 StatusCode ShowerStartFinderTool::ReadSettings(const TiXmlHandle xmlHandle)
661 {
662  PANDORA_RETURN_RESULT_IF_AND_IF(
663  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SpineSlidingFitWindow", m_spineSlidingFitWindow));
664 
665  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
666  XmlHelper::ReadValue(xmlHandle, "LongitudinalCoordinateBinSize", m_longitudinalCoordinateBinSize));
667 
668  PANDORA_RETURN_RESULT_IF_AND_IF(
669  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "nInitialEnergyBins", m_nInitialEnergyBins));
670 
671  PANDORA_RETURN_RESULT_IF_AND_IF(
672  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinSigmaDeviation", m_minSigmaDeviation));
673 
674  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxEdgeGap", m_maxEdgeGap));
675 
676  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
677  XmlHelper::ReadValue(xmlHandle, "LongitudinalShowerFraction", m_longitudinalShowerFraction));
678 
679  PANDORA_RETURN_RESULT_IF_AND_IF(
680  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinShowerOpeningAngle", m_minShowerOpeningAngle));
681 
682  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MolliereRadius", m_molliereRadius));
683 
684  PANDORA_RETURN_RESULT_IF_AND_IF(
685  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ShowerSlidingFitWindow", m_showerSlidingFitWindow));
686 
687  PANDORA_RETURN_RESULT_IF_AND_IF(
688  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxLayerSeparation", m_maxLayerSeparation));
689 
690  return STATUS_CODE_SUCCESS;
691 }
692 
693 } // namespace lar_content
std::map< int, pandora::CaloHitList > LayerToHitMap
float m_longitudinalShowerFraction
The shower region fraction considered.
Header file for the lar two dimensional sliding shower fit result class.
void CharacteriseInitialEnergy(const EnergySpectrumMap &energySpectrumMap, const bool isEndDownstream, float &meanEnergy, float &energySigma) const
Find the mean and standard deviation of the energy depositions in the initial region.
Header file for the pfo helper class.
Header file for the connection pathway helper class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_minShowerOpeningAngle
The min. opening angle of a sensible shower.
void ObtainLongitudinalDecomposition(const TwoDSlidingFitResult &spineTwoDSlidingFit, const pandora::CaloHitList &showerSpineHitList, LongitudinalPositionMap &longitudinalPositionMap) const
Create the [shower spine hit -> shower spine fit longitudinal projection] map.
float m_longitudinalCoordinateBinSize
The longitudinal coordinate bin size.
Header file for the shower start finder tool class.
int m_maxLayerSeparation
The max. allowed separation between the shower start and boundary start layers.
const TwoDSlidingFitResult & GetShowerFitResult() const
Get the sliding fit result for the full shower cluster.
constexpr auto abs(T v)
Returns the absolute value of the argument.
unsigned int m_showerSlidingFitWindow
The sliding window used to fit the shower region.
pandora::StatusCode BuildShowerRegion(const pandora::ParticleFlowObject *const pShowerPfo, const pandora::HitType hitType, const pandora::CaloHitList &showerSpineHitList, const pandora::CartesianVector &showerStartPosition, const pandora::CartesianVector &showerStartDirection, const bool isEndDownstream, pandora::CartesianPointVector &showerRegionPositionVector) const
Build the downstream &#39;shower region&#39; at a given longitudinal distance along the spine.
intermediate_table::const_iterator const_iterator
TFile f
Definition: plotHisto.C:6
unsigned int m_nInitialEnergyBins
The number of longitudinal bins that define the initial region.
Header file for the geometry helper class.
std::map< const pandora::CaloHit *, float > LongitudinalPositionMap
void GetEnergyDistribution(const pandora::CaloHitList &showerSpineHitList, const LongitudinalPositionMap &longitudinalPositionMap, EnergySpectrumMap &energySpectrumMap) const
Create the longituidnal energy distribution.
int GetMaxLayer() const
Get the maximum occupied layer in the sliding fit.
float m_maxEdgeGap
The max. allowed layer gap in a shower boundary.
Header file for the lar two dimensional sliding fit result class.
float m_molliereRadius
The max. distance from the shower core of a collected shower region hit.
std::map< int, LayerFitResult > LayerFitResultMap
pandora::StatusCode GetBoundaryExtremalPoints(const TwoDSlidingShowerFitResult &showerTwoDSlidingFit, const LayerFitResultMap &layerFitResultMap, const pandora::CartesianVector &showerStartPosition, const int showerStartLayer, const int showerEndLayer, pandora::CartesianVector &boundaryEdgeStart, pandora::CartesianVector &boundaryEdgeEnd) const
Determine the start and end positions of a shower boundary.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
const TwoDSlidingFitResult & GetPositiveEdgeFitResult() const
Get the sliding fit result for the positive shower edge.
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.
bool IsShowerTopology(const pandora::ParticleFlowObject *const pShowerPfo, const pandora::HitType hitType, const TwoDSlidingFitResult &spineTwoDSlidingFit, const float longitudinalDistance, const pandora::CaloHitList &showerSpineHitList, const bool isEndDownstream) const
Whether a sensible shower cascade looks to originate at a given position.
pandora::StatusCode CharacteriseShowerTopology(const pandora::CartesianPointVector &showerRegionPositionVector, const pandora::CartesianVector &showerStartPosition, const pandora::HitType hitType, const bool isEndDownstream, const pandora::CartesianVector &showerStartDirection, pandora::CartesianVector &positiveEdgeStart, pandora::CartesianVector &positiveEdgeEnd, pandora::CartesianVector &negativeEdgeStart, pandora::CartesianVector &negativeEdgeEnd, bool &isBetween) const
Parameterise the topological structure of the shower region.
bool IsClockwiseRotation(const pandora::CartesianVector &showerStartDirection, const pandora::CartesianVector &displacementVector) const
Determine whether a point lies on the RHS or LHS (wrt +ve Z) of the shower core.
const TwoDSlidingFitResult & GetNegativeEdgeFitResult() const
Get the sliding fit result for the negative shower edge.
HitType
Definition: HitType.h:12
void GetGlobalPosition(const float rL, const float rT, pandora::CartesianVector &position) const
Get global coordinates for given sliding linear fit coordinates.
float m_minSigmaDeviation
The min. average energy deviation of a candidate shower start.
void ConvertLongitudinalProjectionToGlobal(const TwoDSlidingFitResult &spineTwoDSlidingFit, const float longitudinalDistance, pandora::CartesianVector &globalPosition, pandora::CartesianVector &globalDirection) const
Determine the (X,Y,Z) position and direction at a given longitudinal distance along the spine...
pandora::StatusCode Run(const pandora::ParticleFlowObject *const pShowerPfo, const pandora::CartesianVector &peakDirection, const pandora::HitType hitType, const pandora::CaloHitList &showerSpineHitList, pandora::CartesianVector &showerStartPosition, pandora::CartesianVector &showerStartDirection)
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 GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
int GetLayer(const float rL) const
Get layer number for given sliding linear fit longitudinal coordinate.
unsigned int m_spineSlidingFitWindow
The sliding window used to fit the shower spine.
int FindShowerStartLongitudinalCoordinate(const pandora::ParticleFlowObject *const pShowerPfo, const pandora::HitType hitType, const TwoDSlidingFitResult &spineTwoDSlidingFit, const EnergySpectrumMap &energySpectrumMap, const pandora::CaloHitList &showerSpineHitList, const bool isEndDownstream, const T startIter, const T endIter) const
Find the longitudinal bin which corresponds to the start position of the shower cascade.
void FindShowerStartAndDirection(const pandora::ParticleFlowObject *const pShowerPfo, const pandora::HitType hitType, const TwoDSlidingFitResult &spineTwoDSlidingFit, const EnergySpectrumMap &energySpectrumMap, const pandora::CaloHitList &showerSpineHitList, const bool isEndDownstream, pandora::CartesianVector &showerStartPosition, pandora::CartesianVector &showerStartDirection) const
Find the position at which the shower cascade looks to originate, and its initial direction...