LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrainedVertexSelectionAlgorithm.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
17 
24 
26 
28 
29 #include <random>
30 
31 using namespace pandora;
32 
33 namespace lar_content
34 {
35 
36 TrainedVertexSelectionAlgorithm::TrainedVertexSelectionAlgorithm() :
38  m_trainingSetMode(false),
39  m_allowClassifyDuringTraining(false),
40  m_mcVertexXCorrection(0.f),
41  m_minClusterCaloHits(12),
42  m_slidingFitWindow(100),
43  m_minShowerSpineLength(15.f),
44  m_beamDeweightingConstant(0.4),
45  m_localAsymmetryConstant(3.f),
46  m_globalAsymmetryConstant(1.f),
47  m_showerAsymmetryConstant(1.f),
48  m_energyKickConstant(0.06),
49  m_showerClusteringDistance(3.f),
50  m_minShowerClusterHits(1),
51  m_useShowerClusteringApproximation(false),
52  m_regionRadius(10.f),
53  m_rPhiFineTuningRadius(2.f),
54  m_maxTrueVertexRadius(1.f),
55  m_useRPhiFeatureForRegion(false),
56  m_dropFailedRPhiFastScoreCandidates(true),
57  m_testBeamMode(false),
58  m_legacyEventShapes(true),
59  m_legacyVariables(true)
60 {
61 }
62 
63 //------------------------------------------------------------------------------------------------------------------------------------------
64 
65 void TrainedVertexSelectionAlgorithm::CalculateShowerClusterList(const ClusterList &inputClusterList, ShowerClusterList &showerClusterList) const
66 {
67  ClusterEndPointsMap clusterEndPointsMap;
68  ClusterList showerLikeClusters;
69  this->GetShowerLikeClusterEndPoints(inputClusterList, showerLikeClusters, clusterEndPointsMap);
70 
71  const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
72  ClusterList availableShowerLikeClusters(showerLikeClusters.begin(), showerLikeClusters.end());
73 
74  HitKDTree2D kdTree;
75  HitToClusterMap hitToClusterMap;
76 
78  this->PopulateKdTree(availableShowerLikeClusters, kdTree, hitToClusterMap);
79 
80  while (!availableShowerLikeClusters.empty())
81  {
82  ClusterList showerCluster;
83  showerCluster.push_back(availableShowerLikeClusters.back());
84  availableShowerLikeClusters.pop_back();
85 
86  bool addedCluster(true);
87  while (addedCluster && !availableShowerLikeClusters.empty())
88  {
89  addedCluster = false;
90  for (const Cluster *const pCluster : showerCluster)
91  {
93  {
94  addedCluster = this->AddClusterToShower(kdTree, hitToClusterMap, availableShowerLikeClusters, pCluster, showerCluster);
95  }
96  else
97  {
98  addedCluster = this->AddClusterToShower(clusterEndPointsMap, availableShowerLikeClusters, pCluster, showerCluster);
99  }
100 
101  if (addedCluster)
102  break;
103  }
104  }
105 
106  unsigned int totHits(0);
107  for (const Cluster *const pCluster : showerCluster)
108  totHits += pCluster->GetNCaloHits();
109 
110  if (totHits < m_minClusterCaloHits)
111  continue;
112 
113  showerClusterList.emplace_back(showerCluster, slidingFitPitch, m_slidingFitWindow);
114  }
115 }
116 
117 //------------------------------------------------------------------------------------------------------------------------------------------
118 
120  const ClusterList &clusterList, ClusterList &showerLikeClusters, ClusterEndPointsMap &clusterEndPointsMap) const
121 {
122  for (const Cluster *const pCluster : clusterList)
123  {
124  if (pCluster->GetNCaloHits() < m_minShowerClusterHits)
125  continue;
126 
127  if (this->IsClusterShowerLike(pCluster))
128  showerLikeClusters.push_back(pCluster);
129 
130  CaloHitList clusterCaloHitList;
131  pCluster->GetOrderedCaloHitList().FillCaloHitList(clusterCaloHitList);
132 
133  CaloHitVector clusterCaloHitVector(clusterCaloHitList.begin(), clusterCaloHitList.end());
134  std::sort(clusterCaloHitVector.begin(), clusterCaloHitVector.end(), LArClusterHelper::SortHitsByPosition);
135 
136  if (clusterCaloHitVector.empty())
137  continue;
138 
139  ClusterEndPoints clusterEndPoints(clusterCaloHitVector.front()->GetPositionVector(), clusterCaloHitVector.back()->GetPositionVector());
140  clusterEndPointsMap.emplace(pCluster, clusterEndPoints);
141  }
142 }
143 
144 //------------------------------------------------------------------------------------------------------------------------------------------
145 
146 void TrainedVertexSelectionAlgorithm::PopulateKdTree(const ClusterList &clusterList, HitKDTree2D &kdTree, HitToClusterMap &hitToClusterMap) const
147 {
148  CaloHitList allCaloHits;
149 
150  for (const Cluster *const pCluster : clusterList)
151  {
152  CaloHitList daughterHits;
153  pCluster->GetOrderedCaloHitList().FillCaloHitList(daughterHits);
154  allCaloHits.insert(allCaloHits.end(), daughterHits.begin(), daughterHits.end());
155 
156  for (const CaloHit *const pCaloHit : daughterHits)
157  (void)hitToClusterMap.insert(HitToClusterMap::value_type(pCaloHit, pCluster));
158  }
159 
160  HitKDNode2DList hitKDNode2DList;
161  KDTreeBox hitsBoundingRegion2D(fill_and_bound_2d_kd_tree(allCaloHits, hitKDNode2DList));
162  kdTree.build(hitKDNode2DList, hitsBoundingRegion2D);
163 }
164 
165 //------------------------------------------------------------------------------------------------------------------------------------------
166 
168  ClusterList &availableShowerLikeClusters, const Cluster *const pCluster, ClusterList &showerCluster) const
169 {
170  const auto existingEndPointsIter(clusterEndPointsMap.find(pCluster));
171  if (existingEndPointsIter == clusterEndPointsMap.end())
172  return false;
173 
174  const ClusterEndPoints &existingClusterEndPoints(existingEndPointsIter->second);
175 
176  for (auto iter = availableShowerLikeClusters.begin(); iter != availableShowerLikeClusters.end(); ++iter)
177  {
178  const Cluster *const pAvailableShowerLikeCluster(*iter);
179  const auto &newEndPointsIter(clusterEndPointsMap.find(pAvailableShowerLikeCluster));
180 
181  if (newEndPointsIter == clusterEndPointsMap.end())
182  continue;
183 
184  const ClusterEndPoints &newClusterEndPoints(newEndPointsIter->second);
185  const float startStartDistance((newClusterEndPoints.first - existingClusterEndPoints.first).GetMagnitude());
186  const float startEndDistance((newClusterEndPoints.first - existingClusterEndPoints.second).GetMagnitude());
187  const float endStartDistance((newClusterEndPoints.second - existingClusterEndPoints.first).GetMagnitude());
188  const float endEndDistance((newClusterEndPoints.second - existingClusterEndPoints.second).GetMagnitude());
189 
190  const float smallestDistance(std::min(std::min(startStartDistance, startEndDistance), std::min(endStartDistance, endEndDistance)));
191 
192  if (smallestDistance < m_showerClusteringDistance)
193  {
194  showerCluster.push_back(pAvailableShowerLikeCluster);
195  availableShowerLikeClusters.erase(iter); // Now must return, after invalidating current iterator
196  return true;
197  }
198  }
199 
200  return false;
201 }
202 
203 //------------------------------------------------------------------------------------------------------------------------------------------
204 
206  ClusterList &availableShowerLikeClusters, const Cluster *const pCluster, ClusterList &showerCluster) const
207 {
208  ClusterSet nearbyClusters;
209  CaloHitList daughterHits;
210  pCluster->GetOrderedCaloHitList().FillCaloHitList(daughterHits);
211 
212  for (const CaloHit *const pCaloHit : daughterHits)
213  {
215 
216  HitKDNode2DList found;
217  kdTree.search(searchRegionHits, found);
218 
219  for (const auto &hit : found)
220  (void)nearbyClusters.insert(hitToClusterMap.at(hit.data));
221  }
222 
223  for (auto iter = availableShowerLikeClusters.begin(); iter != availableShowerLikeClusters.end(); ++iter)
224  {
225  const Cluster *const pAvailableShowerLikeCluster(*iter);
226 
227  if ((pAvailableShowerLikeCluster != pCluster) && nearbyClusters.count(pAvailableShowerLikeCluster))
228  {
229  showerCluster.push_back(pAvailableShowerLikeCluster);
230  availableShowerLikeClusters.erase(iter); // Now must return, after invalidating current iterator
231  return true;
232  }
233  }
234 
235  return false;
236 }
237 
238 //------------------------------------------------------------------------------------------------------------------------------------------
239 
241  const ClusterList &clusterListU, const ClusterList &clusterListV, const ClusterList &clusterListW, const VertexVector &vertexVector) const
242 {
243  float eventEnergy(0.f);
244  unsigned int nShoweryHits(0), nHits(0);
245 
246  this->IncrementShoweryParameters(clusterListU, nShoweryHits, nHits, eventEnergy);
247  this->IncrementShoweryParameters(clusterListV, nShoweryHits, nHits, eventEnergy);
248  this->IncrementShoweryParameters(clusterListW, nShoweryHits, nHits, eventEnergy);
249 
250  const unsigned int nClusters(clusterListU.size() + clusterListV.size() + clusterListW.size());
251  const float eventShoweryness((nHits > 0) ? static_cast<float>(nShoweryHits) / static_cast<float>(nHits) : 0.f);
252 
253  ClusterList allClusters(clusterListU);
254  allClusters.insert(allClusters.end(), clusterListV.begin(), clusterListV.end());
255  allClusters.insert(allClusters.end(), clusterListW.begin(), clusterListW.end());
256 
257  const ClusterListMap clusterListMap{{TPC_VIEW_U, clusterListU}, {TPC_VIEW_V, clusterListV}, {TPC_VIEW_W, clusterListW}};
258 
259  float eventArea(0.f), longitudinality(0.f);
261  this->GetLegacyEventShapeFeatures(allClusters, eventArea, longitudinality);
262  else
263  this->GetEventShapeFeatures(clusterListMap, eventArea, longitudinality);
264 
265  return EventFeatureInfo(eventShoweryness, eventEnergy, eventArea, longitudinality, nHits, nClusters, vertexVector.size());
266 }
267 
268 //------------------------------------------------------------------------------------------------------------------------------------------
269 
271  const ClusterList &clusterList, unsigned int &nShoweryHits, unsigned int &nHits, float &eventEnergy) const
272 {
273  for (const Cluster *const pCluster : clusterList)
274  {
275  if (this->IsClusterShowerLike(pCluster))
276  nShoweryHits += pCluster->GetNCaloHits();
277 
278  eventEnergy += pCluster->GetElectromagneticEnergy();
279  nHits += pCluster->GetNCaloHits();
280  }
281 }
282 
283 //------------------------------------------------------------------------------------------------------------------------------------------
284 
285 inline bool TrainedVertexSelectionAlgorithm::IsClusterShowerLike(const Cluster *const pCluster) const
286 {
287  return (pCluster->GetParticleId() == E_MINUS && LArClusterHelper::GetLength(pCluster) < m_minShowerSpineLength);
288 }
289 
290 //------------------------------------------------------------------------------------------------------------------------------------------
291 
292 void TrainedVertexSelectionAlgorithm::GetLegacyEventShapeFeatures(const ClusterList &clusterList, float &eventVolume, float &longitudinality) const
293 {
294  InputFloat xMin, yMin, zMin, xMax, yMax, zMax;
295 
296  for (const Cluster *const pCluster : clusterList)
297  {
298  CartesianVector minPosition(0.f, 0.f, 0.f), maxPosition(0.f, 0.f, 0.f);
299  LArClusterHelper::GetClusterBoundingBox(pCluster, minPosition, maxPosition);
300 
301  this->UpdateSpanCoordinate(minPosition.GetX(), maxPosition.GetX(), xMin, xMax);
302  this->UpdateSpanCoordinate(minPosition.GetY(), maxPosition.GetY(), yMin, yMax);
303  this->UpdateSpanCoordinate(minPosition.GetZ(), maxPosition.GetZ(), zMin, zMax);
304  }
305 
306  const float xSpan(this->GetCoordinateSpan(xMax, xMin));
307  const float ySpan(this->GetCoordinateSpan(yMax, zMin));
308  const float zSpan(this->GetCoordinateSpan(yMax, zMin));
309 
310  // Calculate the volume and longitudinality of the event (ySpan often 0 - to be investigated).
311  if ((xSpan > std::numeric_limits<float>::epsilon()) && (ySpan > std::numeric_limits<float>::epsilon()))
312  {
313  eventVolume = xSpan * ySpan * zSpan;
314  longitudinality = zSpan / (xSpan + ySpan);
315  }
316 
317  else if (ySpan > std::numeric_limits<float>::epsilon())
318  {
319  eventVolume = ySpan * ySpan * zSpan;
320  longitudinality = zSpan / (ySpan + ySpan);
321  }
322 
323  else if (xSpan > std::numeric_limits<float>::epsilon())
324  {
325  eventVolume = xSpan * xSpan * zSpan;
326  longitudinality = zSpan / (xSpan + xSpan);
327  }
328 }
329 
330 //------------------------------------------------------------------------------------------------------------------------------------------
331 
332 void TrainedVertexSelectionAlgorithm::GetEventShapeFeatures(const ClusterListMap &clusterListMap, float &eventArea, float &longitudinality) const
333 {
334  float xSpanU(0.f), zSpanU(0.f), xSpanV(0.f), zSpanV(0.f), xSpanW(0.f), zSpanW(0.f);
335 
336  this->Get2DSpan(clusterListMap.at(TPC_VIEW_U), xSpanU, zSpanU);
337  this->Get2DSpan(clusterListMap.at(TPC_VIEW_V), xSpanV, zSpanV);
338  this->Get2DSpan(clusterListMap.at(TPC_VIEW_W), xSpanW, zSpanW);
339 
340  const float xSpan = (xSpanU + xSpanV + xSpanW) / 3.f;
341  const float zSpan = (zSpanU + zSpanV + zSpanW) / 3.f;
342 
343  if ((xSpan > std::numeric_limits<float>::epsilon()) && (zSpan > std::numeric_limits<float>::epsilon()))
344  {
345  eventArea = xSpan * zSpan;
346  longitudinality = zSpan / (xSpan + zSpan);
347  }
348 }
349 
350 //------------------------------------------------------------------------------------------------------------------------------------------
351 
352 void TrainedVertexSelectionAlgorithm::Get2DSpan(const ClusterList &clusterList, float &xSpan, float &zSpan) const
353 {
354  FloatVector xPositions, zPositions;
355 
356  for (const Cluster *const pCluster : clusterList)
357  {
358  const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
359 
360  for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(), iterEnd = orderedCaloHitList.end(); iter != iterEnd; ++iter)
361  {
362  for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
363  {
364  xPositions.push_back((*hitIter)->GetPositionVector().GetX());
365  zPositions.push_back((*hitIter)->GetPositionVector().GetZ());
366  }
367  }
368  }
369 
370  std::sort(xPositions.begin(), xPositions.end());
371  std::sort(zPositions.begin(), zPositions.end());
372 
373  if (xPositions.empty())
374  {
375  xSpan = 0;
376  zSpan = 0;
377  }
378  else
379  {
380  const int low = std::round(0.05 * xPositions.size());
381  const int high = std::round(0.95 * xPositions.size()) - 1;
382 
383  xSpan = xPositions[high] - xPositions[low];
384  zSpan = zPositions[high] - zPositions[low];
385  }
386 }
387 
388 //------------------------------------------------------------------------------------------------------------------------------------------
389 
391  const float minPositionCoord, const float maxPositionCoord, InputFloat &minCoord, InputFloat &maxCoord) const
392 {
393  if (!minCoord.IsInitialized() || minPositionCoord < minCoord.Get())
394  minCoord = minPositionCoord;
395 
396  if (!maxCoord.IsInitialized() || maxPositionCoord > maxCoord.Get())
397  maxCoord = maxPositionCoord;
398 }
399 
400 //------------------------------------------------------------------------------------------------------------------------------------------
401 
402 inline float TrainedVertexSelectionAlgorithm::GetCoordinateSpan(const InputFloat &minCoord, const InputFloat &maxCoord) const
403 {
404  if (minCoord.IsInitialized() && maxCoord.IsInitialized())
405  return std::fabs(maxCoord.Get() - minCoord.Get());
406 
407  return 0.f;
408 }
409 
410 //------------------------------------------------------------------------------------------------------------------------------------------
411 
413 {
414  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_eventShoweryness));
415  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_eventEnergy));
416  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_eventArea));
417  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_longitudinality));
418  if (this->IsBeamModeOn())
419  {
420  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_nHits));
421  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_nClusters));
422  featureVector.push_back(static_cast<double>(eventFeatureInfo.m_nCandidates));
423  }
424 }
425 
426 //------------------------------------------------------------------------------------------------------------------------------------------
427 
429  const SlidingFitDataListMap &slidingFitDataListMap, const ShowerClusterListMap &showerClusterListMap, const KDTreeMap &kdTreeMap,
430  const Vertex *const pVertex, VertexFeatureInfoMap &vertexFeatureInfoMap) const
431 {
432  float bestFastScore(-std::numeric_limits<float>::max()); // not actually used - artefact of toolizing RPhi score and still using performance trick
433 
434  // ATTN - If beam mode is false GetBeamDeweightingScore will fail, so have a default value that we'll ignore when poplating the feature vector
435  double tempBeamDeweight{0.f};
436  if (this->IsBeamModeOn())
437  tempBeamDeweight = this->GetBeamDeweightingScore(beamConstants, pVertex);
438 
439  const double beamDeweighting(tempBeamDeweight);
440 
441  const double energyKick(LArMvaHelper::CalculateFeaturesOfType<EnergyKickFeatureTool>(m_featureToolVector, this, pVertex,
442  slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore)
443  .at(0)
444  .Get());
445 
446  const double localAsymmetry(LArMvaHelper::CalculateFeaturesOfType<LocalAsymmetryFeatureTool>(m_featureToolVector, this, pVertex,
447  slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore)
448  .at(0)
449  .Get());
450 
451  const double globalAsymmetry(LArMvaHelper::CalculateFeaturesOfType<GlobalAsymmetryFeatureTool>(m_featureToolVector, this, pVertex,
452  slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore)
453  .at(0)
454  .Get());
455 
456  const double showerAsymmetry(LArMvaHelper::CalculateFeaturesOfType<ShowerAsymmetryFeatureTool>(m_featureToolVector, this, pVertex,
457  slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore)
458  .at(0)
459  .Get());
460 
461  //const double rPhiFeature(LArMvaHelper::CalculateFeaturesOfType<RPhiFeatureTool>(m_featureToolVector, this, pVertex,
462  // slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore).at(0).Get());
463 
464  double dEdxAsymmetry(0.f), vertexEnergy(0.f);
465 
466  if (!m_legacyVariables)
467  {
468  dEdxAsymmetry = LArMvaHelper::CalculateFeaturesOfType<EnergyDepositionAsymmetryFeatureTool>(m_featureToolVector, this, pVertex,
469  slidingFitDataListMap, clusterListMap, kdTreeMap, showerClusterListMap, beamDeweighting, bestFastScore)
470  .at(0)
471  .Get();
472 
473  vertexEnergy = this->GetVertexEnergy(pVertex, kdTreeMap);
474  }
475 
476  VertexFeatureInfo vertexFeatureInfo(beamDeweighting, 0.f, energyKick, localAsymmetry, globalAsymmetry, showerAsymmetry, dEdxAsymmetry, vertexEnergy);
477  vertexFeatureInfoMap.emplace(pVertex, vertexFeatureInfo);
478 }
479 
480 //------------------------------------------------------------------------------------------------------------------------------------------
481 
483  VertexFeatureInfoMap &vertexFeatureInfoMap, const Vertex *const pVertex, VertexScoreList &initialScoreList) const
484 {
485  VertexFeatureInfo vertexFeatureInfo = vertexFeatureInfoMap.at(pVertex);
486 
487  const float beamDeweightingScore(vertexFeatureInfo.m_beamDeweighting / m_beamDeweightingConstant);
488  const float energyKickScore(-vertexFeatureInfo.m_energyKick / m_energyKickConstant);
489  const float localAsymmetryScore(vertexFeatureInfo.m_localAsymmetry / m_localAsymmetryConstant);
490  const float globalAsymmetryScore(vertexFeatureInfo.m_globalAsymmetry / m_globalAsymmetryConstant);
491  const float showerAsymmetryScore(vertexFeatureInfo.m_showerAsymmetry / m_showerAsymmetryConstant);
492 
493  initialScoreList.emplace_back(pVertex, beamDeweightingScore + energyKickScore + localAsymmetryScore + globalAsymmetryScore + showerAsymmetryScore);
494 }
495 
496 //------------------------------------------------------------------------------------------------------------------------------------------
497 
499 {
500  std::sort(initialScoreList.begin(), initialScoreList.end());
501 
502  for (const VertexScore &vertexScore : initialScoreList)
503  {
504  const Vertex *const pVertex(vertexScore.GetVertex());
505  bool farEnoughAway(true);
506 
507  for (const Vertex *const pRegionVertex : bestRegionVertices)
508  {
509  if (pRegionVertex == pVertex)
510  {
511  farEnoughAway = false;
512  break;
513  }
514 
515  const float distance = (pRegionVertex->GetPosition() - pVertex->GetPosition()).GetMagnitude();
516 
517  if (distance <= m_regionRadius)
518  {
519  farEnoughAway = false;
520  break;
521  }
522  }
523 
524  if (farEnoughAway)
525  bestRegionVertices.push_back(pVertex);
526  }
527 }
528 
529 //------------------------------------------------------------------------------------------------------------------------------------------
530 
531 void TrainedVertexSelectionAlgorithm::ProduceTrainingSets(const VertexVector &vertexVector, const VertexVector &bestRegionVertices,
532  VertexFeatureInfoMap &vertexFeatureInfoMap, const LArMvaHelper::MvaFeatureVector &eventFeatureList, const KDTreeMap &kdTreeMap) const
533 {
534  if (vertexVector.empty())
535  return;
536 
537  // Create a distribution for random coin flips.
538  std::random_device device;
539  std::mt19937 generator(device());
540  std::bernoulli_distribution coinFlip(0.5);
541 
542  const std::string interactionType(this->GetInteractionType());
543 
544  // Produce training examples for the vertices representing regions.
545  const Vertex *const pBestRegionVertex(this->ProduceTrainingExamples(bestRegionVertices, vertexFeatureInfoMap, coinFlip, generator,
546  interactionType, m_trainingOutputFileRegion, eventFeatureList, kdTreeMap, m_regionRadius, m_useRPhiFeatureForRegion));
547 
548  // Get all the vertices in the best region.
549  VertexVector regionalVertices{pBestRegionVertex};
550  for (const Vertex *const pVertex : vertexVector)
551  {
552  if (pVertex == pBestRegionVertex)
553  continue;
554 
555  if ((pBestRegionVertex->GetPosition() - pVertex->GetPosition()).GetMagnitude() < m_regionRadius)
556  regionalVertices.push_back(pVertex);
557  }
558 
559  this->CalculateRPhiScores(regionalVertices, vertexFeatureInfoMap, kdTreeMap);
560 
561  // Produce training examples for the final vertices within the best region.
562  if (!regionalVertices.empty())
563  {
564  this->ProduceTrainingExamples(regionalVertices, vertexFeatureInfoMap, coinFlip, generator, interactionType,
565  m_trainingOutputFileVertex, eventFeatureList, kdTreeMap, m_maxTrueVertexRadius, true);
566  }
567 }
568 
569 //------------------------------------------------------------------------------------------------------------------------------------------
570 
572  VertexVector &vertexVector, VertexFeatureInfoMap &vertexFeatureInfoMap, const KDTreeMap &kdTreeMap) const
573 {
574  float bestFastScore(-std::numeric_limits<float>::max());
575 
576  for (auto iter = vertexVector.begin(); iter != vertexVector.end(); /* no increment */)
577  {
578  VertexFeatureInfo &vertexFeatureInfo = vertexFeatureInfoMap.at(*iter);
579  vertexFeatureInfo.m_rPhiFeature = static_cast<float>(LArMvaHelper::CalculateFeaturesOfType<RPhiFeatureTool>(m_featureToolVector, this,
580  *iter, SlidingFitDataListMap(), ClusterListMap(), kdTreeMap, ShowerClusterListMap(), vertexFeatureInfo.m_beamDeweighting, bestFastScore)
581  .at(0)
582  .Get());
583 
584  if (m_dropFailedRPhiFastScoreCandidates && (vertexFeatureInfo.m_rPhiFeature <= std::numeric_limits<float>::epsilon()))
585  iter = vertexVector.erase(iter);
586 
587  else
588  ++iter;
589  }
590 }
591 
592 //------------------------------------------------------------------------------------------------------------------------------------------
593 
595 {
596  // Extract input collections
597  const MCParticleList *pMCParticleList(nullptr);
598  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));
599 
600  const CaloHitList *pCaloHitList(nullptr);
601  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_caloHitListName, pCaloHitList));
602 
603  LArMCParticleHelper::MCContributionMap targetMCParticlesToGoodHitsMap;
604 
605  if (m_testBeamMode)
606  {
608  LArMCParticleHelper::IsTriggeredBeamParticle, targetMCParticlesToGoodHitsMap);
609  }
610  else
611  {
612  // ATTN Assumes single neutrino is parent of all neutrino-induced mc particles
614  LArMCParticleHelper::IsBeamNeutrinoFinalState, targetMCParticlesToGoodHitsMap);
615  }
616 
617  MCParticleList mcPrimaryList;
618  for (const auto &mapEntry : targetMCParticlesToGoodHitsMap)
619  mcPrimaryList.push_back(mapEntry.first);
620  mcPrimaryList.sort(LArMCParticleHelper::SortByMomentum);
621 
623  return LArInteractionTypeHelper::ToString(interactionType);
624 }
625 
626 //------------------------------------------------------------------------------------------------------------------------------------------
627 
629  const VertexFeatureInfoMap &vertexFeatureInfoMap, std::bernoulli_distribution &coinFlip, std::mt19937 &generator,
630  const std::string &interactionType, const std::string &trainingOutputFile, const LArMvaHelper::MvaFeatureVector &eventFeatureList,
631  const KDTreeMap &kdTreeMap, const float maxRadius, const bool useRPhi) const
632 {
633  const Vertex *pBestVertex(nullptr);
634  float bestVertexDr(std::numeric_limits<float>::max());
635 
636  LArMvaHelper::MvaFeatureVector bestVertexFeatureList;
637  this->GetBestVertex(vertexVector, pBestVertex, bestVertexDr);
638 
639  VertexFeatureInfo bestVertexFeatureInfo(vertexFeatureInfoMap.at(pBestVertex));
640  this->AddVertexFeaturesToVector(bestVertexFeatureInfo, bestVertexFeatureList, useRPhi);
641 
642  for (const Vertex *const pVertex : vertexVector)
643  {
644  if (pVertex == pBestVertex)
645  continue;
646 
647  LArMvaHelper::MvaFeatureVector featureList;
648  VertexFeatureInfo vertexFeatureInfo(vertexFeatureInfoMap.at(pVertex));
649  this->AddVertexFeaturesToVector(vertexFeatureInfo, featureList, useRPhi);
650 
651  if (!m_legacyVariables)
652  {
653  LArMvaHelper::MvaFeatureVector sharedFeatureList;
654  float separation(0.f), axisHits(0.f);
655  this->GetSharedFeatures(pVertex, pBestVertex, kdTreeMap, separation, axisHits);
656  VertexSharedFeatureInfo sharedFeatureInfo(separation, axisHits);
657  this->AddSharedFeaturesToVector(sharedFeatureInfo, sharedFeatureList);
658 
659  if (pBestVertex && (bestVertexDr < maxRadius))
660  {
661  if (coinFlip(generator))
662  LArMvaHelper::ProduceTrainingExample(trainingOutputFile + "_" + interactionType + ".txt", true,
663  LArMvaHelper::ConcatenateFeatureLists(eventFeatureList, bestVertexFeatureList, featureList, sharedFeatureList));
664  else
665  LArMvaHelper::ProduceTrainingExample(trainingOutputFile + "_" + interactionType + ".txt", false,
666  LArMvaHelper::ConcatenateFeatureLists(eventFeatureList, featureList, bestVertexFeatureList, sharedFeatureList));
667  }
668  }
669  else
670  {
671  if (pBestVertex && (bestVertexDr < maxRadius))
672  {
673  if (coinFlip(generator))
674  LArMvaHelper::ProduceTrainingExample(trainingOutputFile + "_" + interactionType + ".txt", true,
675  LArMvaHelper::ConcatenateFeatureLists(eventFeatureList, bestVertexFeatureList, featureList));
676  else
677  LArMvaHelper::ProduceTrainingExample(trainingOutputFile + "_" + interactionType + ".txt", false,
678  LArMvaHelper::ConcatenateFeatureLists(eventFeatureList, featureList, bestVertexFeatureList));
679  }
680  }
681  }
682 
683  return pBestVertex;
684 }
685 
686 //------------------------------------------------------------------------------------------------------------------------------------------
687 
689  const Vertex *const pVertex1, const Vertex *const pVertex2, const KDTreeMap &kdTreeMap, float &separation, float &axisHits) const
690 {
691  separation = (pVertex1->GetPosition() - pVertex2->GetPosition()).GetMagnitude();
692 
693  this->IncrementSharedAxisValues(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex1->GetPosition(), TPC_VIEW_U),
694  LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex2->GetPosition(), TPC_VIEW_U), kdTreeMap.at(TPC_VIEW_U), axisHits);
695 
696  this->IncrementSharedAxisValues(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex1->GetPosition(), TPC_VIEW_V),
697  LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex2->GetPosition(), TPC_VIEW_V), kdTreeMap.at(TPC_VIEW_V), axisHits);
698 
699  this->IncrementSharedAxisValues(LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex1->GetPosition(), TPC_VIEW_W),
700  LArGeometryHelper::ProjectPosition(this->GetPandora(), pVertex2->GetPosition(), TPC_VIEW_W), kdTreeMap.at(TPC_VIEW_W), axisHits);
701 
702  axisHits = separation > std::numeric_limits<float>::epsilon() ? axisHits / separation : 0.f;
703 }
704 
705 //------------------------------------------------------------------------------------------------------------------------------------------
706 
708  const CartesianVector pos1, const CartesianVector pos2, HitKDTree2D &kdTree, float &axisHits) const
709 {
710  if (pos1 == pos2)
711  return;
712 
713  // Define the axis and perpendicular directions
714  const CartesianVector unitAxis = (pos1 - pos2).GetUnitVector();
715  const CartesianVector unitPerp(unitAxis.GetZ(), 0, -unitAxis.GetX());
716 
717  // Define the corners of the search box
718  const CartesianVector point1 = pos1 + unitPerp;
719  const CartesianVector point2 = pos1 - unitPerp;
720  const CartesianVector point3 = pos2 + unitPerp;
721  const CartesianVector point4 = pos2 - unitPerp;
722 
723  // Find the total coordinate span these points cover
724  const float xMin{std::min({point1.GetX(), point2.GetX(), point3.GetX(), point4.GetX()})};
725  const float xMax{std::max({point1.GetX(), point2.GetX(), point3.GetX(), point4.GetX()})};
726  const float zMin{std::min({point1.GetZ(), point2.GetZ(), point3.GetZ(), point4.GetZ()})};
727  const float zMax{std::max({point1.GetZ(), point2.GetZ(), point3.GetZ(), point4.GetZ()})};
728 
729  // Use a kd search to find the hits in the 'wide' area
730  KDTreeBox searchBox(xMin, zMin, xMax, zMax);
731  HitKDNode2DList found;
732  kdTree.search(searchBox, found);
733 
734  // Use IsHitInBox method to check the 'wide' area hits for those in the search box
735  for (auto f : found)
736  {
737  const CartesianVector &hitPos = f.data->GetPositionVector();
738  bool inBox = this->IsHitInBox(hitPos, point1, point2, point3, point4);
739 
740  if (inBox)
741  ++axisHits;
742  }
743 }
744 
745 //------------------------------------------------------------------------------------------------------------------------------------------
746 
747 bool TrainedVertexSelectionAlgorithm::IsHitInBox(const CartesianVector &hitPos, const CartesianVector &point1,
748  const CartesianVector &point2, const CartesianVector &point3, const CartesianVector &point4) const
749 {
750  bool b1 = std::signbit(((point2 - point1).GetCrossProduct(point2 - hitPos)).GetY());
751  bool b2 = std::signbit(((point4 - point3).GetCrossProduct(point4 - hitPos)).GetY());
752 
753  if (!(b1 xor b2))
754  return false;
755 
756  bool b3 = std::signbit(((point3 - point1).GetCrossProduct(point3 - hitPos)).GetY());
757  bool b4 = std::signbit(((point4 - point2).GetCrossProduct(point4 - hitPos)).GetY());
758 
759  return (b3 xor b4);
760 }
761 
762 //------------------------------------------------------------------------------------------------------------------------------------------
763 
764 void TrainedVertexSelectionAlgorithm::GetBestVertex(const VertexVector &vertexVector, const Vertex *&pBestVertex, float &bestVertexDr) const
765 {
766  // Extract input collections
767  const MCParticleList *pMCParticleList(nullptr);
768  PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));
769 
770  // Obtain vector: true targets
771  MCParticleVector mcTargetVector;
772 
773  if (m_testBeamMode)
774  {
775  LArMCParticleHelper::GetTrueTestBeamParticles(pMCParticleList, mcTargetVector);
776  }
777  else
778  {
779  LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, mcTargetVector);
780  }
781 
782  for (const Vertex *const pVertex : vertexVector)
783  {
784  float mcVertexDr(std::numeric_limits<float>::max());
785  for (const MCParticle *const pMCTarget : mcTargetVector)
786  {
787  const CartesianVector &mcTargetPosition(
788  (m_testBeamMode && std::fabs(pMCTarget->GetParticleId()) == 11) ? pMCTarget->GetVertex() : pMCTarget->GetEndpoint());
789  const CartesianVector mcTargetCorrectedPosition(
790  mcTargetPosition.GetX() + m_mcVertexXCorrection, mcTargetPosition.GetY(), mcTargetPosition.GetZ());
791  const float dr = (mcTargetCorrectedPosition - pVertex->GetPosition()).GetMagnitude();
792 
793  if (dr < mcVertexDr)
794  mcVertexDr = dr;
795  }
796 
797  if (mcVertexDr < bestVertexDr)
798  {
799  bestVertexDr = mcVertexDr;
800  pBestVertex = pVertex;
801  }
802  }
803 }
804 
805 //------------------------------------------------------------------------------------------------------------------------------------------
806 
808  const VertexFeatureInfo &vertexFeatureInfo, LArMvaHelper::MvaFeatureVector &featureVector, const bool useRPhi) const
809 {
810  if (this->IsBeamModeOn())
811  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_beamDeweighting));
812  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_energyKick));
813  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_globalAsymmetry));
814  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_localAsymmetry));
815  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_showerAsymmetry));
816 
817  if (!m_legacyVariables)
818  {
819  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_dEdxAsymmetry));
820  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_vertexEnergy));
821  }
822 
823  if (useRPhi)
824  featureVector.push_back(static_cast<double>(vertexFeatureInfo.m_rPhiFeature));
825 }
826 
827 //------------------------------------------------------------------------------------------------------------------------------------------
828 
830  const VertexSharedFeatureInfo &vertexSharedFeatureInfo, LArMvaHelper::MvaFeatureVector &featureVector) const
831 {
832  featureVector.push_back(static_cast<double>(vertexSharedFeatureInfo.m_separation));
833  featureVector.push_back(static_cast<double>(vertexSharedFeatureInfo.m_axisHits));
834 }
835 
836 //------------------------------------------------------------------------------------------------------------------------------------------
837 
839  const Vertex *const pFavouriteVertex, const VertexVector &vertexVector, VertexScoreList &finalVertexScoreList) const
840 {
841  if (pFavouriteVertex)
842  {
843  const CartesianVector vertexPosition(pFavouriteVertex->GetPosition());
844 
845  for (const Vertex *const pVertex : vertexVector)
846  {
847  if ((pVertex->GetPosition() - vertexPosition).GetMagnitude() < m_rPhiFineTuningRadius)
848  {
849  const float rPhiScore(vertexFeatureInfoMap.at(pVertex).m_rPhiFeature);
850  finalVertexScoreList.emplace_back(pVertex, rPhiScore);
851  }
852  }
853  }
854 }
855 
856 //------------------------------------------------------------------------------------------------------------------------------------------
857 //------------------------------------------------------------------------------------------------------------------------------------------
858 
859 StatusCode TrainedVertexSelectionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
860 {
861  AlgorithmToolVector algorithmToolVector;
862  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ProcessAlgorithmToolList(*this, xmlHandle, "FeatureTools", algorithmToolVector));
863 
864  for (AlgorithmTool *const pAlgorithmTool : algorithmToolVector)
865  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, LArMvaHelper::AddFeatureToolToVector(pAlgorithmTool, m_featureToolVector));
866 
867  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "TrainingSetMode", m_trainingSetMode));
868 
869  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
870  XmlHelper::ReadValue(xmlHandle, "AllowClassifyDuringTraining", m_allowClassifyDuringTraining));
871 
872  PANDORA_RETURN_RESULT_IF_AND_IF(
873  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MCVertexXCorrection", m_mcVertexXCorrection));
874 
875  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
876  XmlHelper::ReadValue(xmlHandle, "TrainingOutputFileRegion", m_trainingOutputFileRegion));
877 
878  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
879  XmlHelper::ReadValue(xmlHandle, "TrainingOutputFileVertex", m_trainingOutputFileVertex));
880 
882  {
883  std::cout << "TrainedVertexSelectionAlgorithm: TrainingOutputFileRegion and TrainingOutputFileVertex are required for training set "
884  << "mode" << std::endl;
885  return STATUS_CODE_INVALID_PARAMETER;
886  }
887 
888  PANDORA_RETURN_RESULT_IF_AND_IF(
889  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MCParticleListName", m_mcParticleListName));
890 
891  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CaloHitListName", m_caloHitListName));
892 
893  if (m_trainingSetMode && (m_mcParticleListName.empty() || m_caloHitListName.empty()))
894  {
895  std::cout << "TrainedVertexSelectionAlgorithm: MCParticleListName and CaloHitListName are required for training set mode" << std::endl;
896  return STATUS_CODE_INVALID_PARAMETER;
897  }
898 
899  PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));
900 
901  PANDORA_RETURN_RESULT_IF_AND_IF(
902  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterCaloHits", m_minClusterCaloHits));
903 
904  PANDORA_RETURN_RESULT_IF_AND_IF(
905  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitWindow", m_slidingFitWindow));
906 
907  PANDORA_RETURN_RESULT_IF_AND_IF(
908  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinShowerSpineLength", m_minShowerSpineLength));
909 
910  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
911  XmlHelper::ReadValue(xmlHandle, "BeamDeweightingConstant", m_beamDeweightingConstant));
912 
913  PANDORA_RETURN_RESULT_IF_AND_IF(
914  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LocalAsymmetryConstant", m_localAsymmetryConstant));
915 
916  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
917  XmlHelper::ReadValue(xmlHandle, "GlobalAsymmetryConstant", m_globalAsymmetryConstant));
918 
919  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
920  XmlHelper::ReadValue(xmlHandle, "ShowerAsymmetryConstant", m_showerAsymmetryConstant));
921 
922  PANDORA_RETURN_RESULT_IF_AND_IF(
923  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "EnergyKickConstant", m_energyKickConstant));
924 
925  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
926  XmlHelper::ReadValue(xmlHandle, "ShowerClusteringDistance", m_showerClusteringDistance));
927 
928  PANDORA_RETURN_RESULT_IF_AND_IF(
929  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinShowerClusterHits", m_minShowerClusterHits));
930 
931  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
932  XmlHelper::ReadValue(xmlHandle, "UseShowerClusteringApproximation", m_useShowerClusteringApproximation));
933 
934  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "RegionRadius", m_regionRadius));
935 
936  PANDORA_RETURN_RESULT_IF_AND_IF(
937  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "RPhiFineTuningRadius", m_rPhiFineTuningRadius));
938 
939  PANDORA_RETURN_RESULT_IF_AND_IF(
940  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxTrueVertexRadius", m_maxTrueVertexRadius));
941 
942  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
943  XmlHelper::ReadValue(xmlHandle, "UseRPhiFeatureForRegion", m_useRPhiFeatureForRegion));
944 
945  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
946  XmlHelper::ReadValue(xmlHandle, "DropFailedRPhiFastScoreCandidates", m_dropFailedRPhiFastScoreCandidates));
947 
948  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "TestBeamMode", m_testBeamMode));
949 
950  PANDORA_RETURN_RESULT_IF_AND_IF(
951  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LegacyEventShapes", m_legacyEventShapes));
952 
953  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LegacyVariables", m_legacyVariables));
954 
956  std::cout << "TrainedVertexSelectionAlgorithm: WARNING -- Producing training sample using incorrect legacy event shapes, consider turning LegacyEventShapes off"
957  << std::endl;
958 
960 }
961 
962 } // namespace lar_content
void AddEventFeaturesToVector(const EventFeatureInfo &eventFeatureInfo, LArMvaHelper::MvaFeatureVector &featureVector) const
Add the event features to a vector in the correct order.
void GetLegacyEventShapeFeatures(const pandora::ClusterList &clusterList, float &eventVolume, float &longitudinality) const
Get the event shape features.
std::string m_trainingOutputFileVertex
The training output file for the vertex mva.
unsigned int m_nCandidates
The total number of vertex candidates.
bool m_useShowerClusteringApproximation
Whether to use the shower clustering distance approximation.
bool m_useRPhiFeatureForRegion
Whether to use the r/phi feature for the region vertex.
Header file for the kd tree linker algo template class.
void GetBestRegionVertices(VertexScoreList &initialScoreList, pandora::VertexVector &bestRegionVertices) const
Get the list of top-N separated vertices.
std::unordered_map< const pandora::CaloHit *, const pandora::Cluster * > HitToClusterMap
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
MvaTypes::MvaFeatureVector MvaFeatureVector
Definition: LArMvaHelper.h:75
unsigned int m_minClusterCaloHits
The min number of hits parameter in the energy score.
Header file for the interaction type helper class.
void PopulateInitialScoreList(VertexFeatureInfoMap &vertexFeatureInfoMap, const pandora::Vertex *const pVertex, VertexScoreList &initialScoreList) const
Populate the initial vertex score list for a given vertex.
Header file for the energy deposition asymmetry feature tool class.
static pandora::StatusCode ProduceTrainingExample(const std::string &trainingOutputFile, const bool result, TCONTAINER &&featureContainer)
Produce a training example with the given features and result.
Definition: LArMvaHelper.h:285
VertexFeatureTool::FeatureToolVector m_featureToolVector
The feature tool vector.
Box structure used to define 2D field. It&#39;s used in KDTree building step to divide the detector space...
void PopulateFinalVertexScoreList(const VertexFeatureInfoMap &vertexFeatureInfoMap, const pandora::Vertex *const pFavouriteVertex, const pandora::VertexVector &vertexVector, VertexScoreList &finalVertexScoreList) const
Populate the final vertex score list using the r/phi score to find the best vertex in the vicinity...
bool IsBeamModeOn() const
Whether algorithm is running in beam mode, assuming neutrinos travel in positive z-direction.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
bool m_legacyEventShapes
Whether to use the old event shapes calculation.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
intermediate_table::const_iterator const_iterator
void AddVertexFeaturesToVector(const VertexFeatureInfo &vertexFeatureInfo, LArMvaHelper::MvaFeatureVector &featureVector, const bool useRPhi) const
Add the vertex features to a vector in the correct order.
void IncrementSharedAxisValues(const pandora::CartesianVector pos1, const pandora::CartesianVector pos2, HitKDTree2D &kdTree, float &axisHits) const
Increments the axis hits information for one view.
const pandora::Vertex * ProduceTrainingExamples(const pandora::VertexVector &vertexVector, const VertexFeatureInfoMap &vertexFeatureInfoMap, std::bernoulli_distribution &coinFlip, std::mt19937 &generator, const std::string &interactionType, const std::string &trainingOutputFile, const LArMvaHelper::MvaFeatureVector &eventFeatureList, const KDTreeMap &kdTreeMap, const float maxRadius, const bool useRPhi) const
Produce a set of training examples for a binary classifier.
bool IsClusterShowerLike(const pandora::Cluster *const pCluster) const
Find whether a cluster is shower-like.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static InteractionType GetInteractionType(const pandora::MCParticleList &mcPrimaryList)
Get the interaction type of an event.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
void CalculateRPhiScores(pandora::VertexVector &vertexVector, VertexFeatureInfoMap &vertexFeatureInfoMap, const KDTreeMap &kdTreeMap) const
Calculate the r/phi scores for the vertices in a vector, possibly erasing those that fail the fast sc...
Header file for the trained vertex selection algorithm class.
void Get2DSpan(const pandora::ClusterList &clusterList, float &xSpan, float &zSpan) const
Get the coordinate span in one view.
float m_localAsymmetryConstant
The local asymmetry constant for the initial region score list.
bool IsHitInBox(const pandora::CartesianVector &hitPos, const pandora::CartesianVector &point1, const pandora::CartesianVector &point2, const pandora::CartesianVector &point3, const pandora::CartesianVector &point4) const
Determines whether a hit lies within the box defined by four other positions.
unsigned int m_minShowerClusterHits
The minimum number of shower cluster hits.
TFile f
Definition: plotHisto.C:6
Header file for the geometry helper class.
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
static void GetClusterBoundingBox(const pandora::Cluster *const pCluster, pandora::CartesianVector &minimumCoordinate, pandora::CartesianVector &maximumCoordinate)
Get minimum and maximum X, Y and Z positions of the calo hits in a cluster.
void ProduceTrainingSets(const pandora::VertexVector &vertexVector, const pandora::VertexVector &bestRegionVertices, VertexFeatureInfoMap &vertexFeatureInfoMap, const LArMvaHelper::MvaFeatureVector &eventFeatureList, const KDTreeMap &kdTreeMap) const
Produce the region and vertex training sets.
static void SelectReconstructableMCParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const PrimaryParameters &parameters, std::function< bool(const pandora::MCParticle *const)> fCriteria, MCContributionMap &selectedMCParticlesToHitsMap)
Select target, reconstructable mc particles that match given criteria.
std::map< pandora::HitType, const ShowerClusterList > ShowerClusterListMap
Map of shower cluster lists for passing to tools.
std::string m_mcParticleListName
The MC particle list for creating training examples.
void PopulateVertexFeatureInfoMap(const BeamConstants &beamConstants, const ClusterListMap &clusterListMap, const SlidingFitDataListMap &slidingFitDataListMap, const ShowerClusterListMap &showerClusterListMap, const KDTreeMap &kdTreeMap, const pandora::Vertex *const pVertex, VertexFeatureInfoMap &vertexFeatureInfoMap) const
Populate the vertex feature info map for a given vertex.
float m_globalAsymmetryConstant
The global asymmetry constant for the initial region score list.
Header file for the lar monte carlo particle helper helper class.
std::pair< pandora::CartesianVector, pandora::CartesianVector > ClusterEndPoints
static bool SortHitsByPosition(const pandora::CaloHit *const pLhs, const pandora::CaloHit *const pRhs)
Sort calo hits by their position (use Z, followed by X, followed by Y)
Header file for the cluster helper class.
static std::string ToString(const InteractionType interactionType)
Get a string representation of an interaction type.
EventFeatureInfo CalculateEventFeatures(const pandora::ClusterList &clusterListU, const pandora::ClusterList &clusterListV, const pandora::ClusterList &clusterListW, const pandora::VertexVector &vertexVector) const
Calculate the event parameters.
void build(std::vector< KDTreeNodeInfoT< DATA, DIM >> &eltList, const KDTreeBoxT< DIM > &region)
Build the KD tree from the "eltList" in the space define by "region".
float m_maxTrueVertexRadius
The maximum distance at which a vertex candidate can be considered the &#39;true&#39; vertex.
std::map< pandora::HitType, const pandora::ClusterList & > ClusterListMap
Map array of cluster lists for passing to tools.
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
static float GetLength(const pandora::Cluster *const pCluster)
Get length of cluster.
static pandora::StatusCode AddFeatureToolToVector(pandora::AlgorithmTool *const pFeatureTool, MvaFeatureToolVector< Ts... > &featureToolVector)
Add a feature tool to a vector of feature tools.
Definition: LArMvaHelper.h:451
float m_showerAsymmetryConstant
The shower asymmetry constant for the initial region score list.
static bool IsTriggeredBeamParticle(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary triggered beam MCParticle.
Header file for the local asymmetry feature tool class.
Header file for the file helper class.
void PopulateKdTree(const pandora::ClusterList &clusterList, HitKDTree2D &kdTree, HitToClusterMap &hitToClusterMap) const
Populate kd tree with information about hits in a provided list of clusters.
Detector simulation of raw signals on wires.
float m_beamDeweightingConstant
The beam deweighting constant for the initial region score list.
bool m_allowClassifyDuringTraining
Whether classification is allowed during training.
void CalculateShowerClusterList(const pandora::ClusterList &inputClusterList, ShowerClusterList &showerClusterList) const
Calculate the shower cluster map for a cluster list.
std::string m_trainingOutputFileRegion
The training output file for the region mva.
bool m_dropFailedRPhiFastScoreCandidates
Whether to drop candidates that fail the r/phi fast score test.
void AddSharedFeaturesToVector(const VertexSharedFeatureInfo &vertexSharedFeatureInfo, LArMvaHelper::MvaFeatureVector &featureVector) const
Add the shared features to a vector in the correct order.
float m_energyKickConstant
The energy kick constant for the initial region score list.
void GetSharedFeatures(const pandora::Vertex *const pVertex1, const pandora::Vertex *const pVertex2, const KDTreeMap &kdTreeMap, float &separation, float &axisHits) const
Calculates the shared features of a pair of vertex candidates.
static void GetTrueNeutrinos(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueNeutrinos)
Get neutrino MC particles from an input MC particle list.
std::map< pandora::HitType, const SlidingFitDataList > SlidingFitDataListMap
Map of sliding fit data lists for passing to tools.
float GetVertexEnergy(const pandora::Vertex *const pVertex, const KDTreeMap &kdTreeMap) const
Calculate the energy of a vertex candidate by summing values from all three planes.
void UpdateSpanCoordinate(const float minPositionCoord, const float maxPositionCoord, pandora::InputFloat &minCoord, pandora::InputFloat &maxCoord) const
Update the min/max coordinate spans.
float m_mcVertexXCorrection
The correction to the x-coordinate of the MC vertex position.
void IncrementShoweryParameters(const pandora::ClusterList &clusterList, unsigned int &nShoweryHits, unsigned int &nHits, float &eventEnergy) const
Increment the showery hit parameters for a cluster list.
float GetCoordinateSpan(const pandora::InputFloat &minCoord, const pandora::InputFloat &maxCoord) const
Get the coordinate span.
bool m_legacyVariables
Whether to only use the old variables.
unsigned int m_slidingFitWindow
The layer window for the sliding linear fits.
std::map< const pandora::Cluster *const, ClusterEndPoints > ClusterEndPointsMap
void GetEventShapeFeatures(const ClusterListMap &clusterListMap, float &eventArea, float &longitudinality) const
Get the event shape features.
std::string m_caloHitListName
The 2D CaloHit list name.
std::vector< art::Ptr< recob::Vertex > > VertexVector
KDTreeBox fill_and_bound_2d_kd_tree(const MANAGED_CONTAINER< const T * > &points, std::vector< KDTreeNodeInfoT< const T *, 2 >> &nodes)
fill_and_bound_2d_kd_tree
Header file for the shower asymmetry feature tool class.
bool AddClusterToShower(const ClusterEndPointsMap &clusterEndPointsMap, pandora::ClusterList &availableShowerLikeClusters, const pandora::Cluster *const pCluster, pandora::ClusterList &showerCluster) const
Try to add an available cluster to a given shower cluster, using shower clustering approximation...
float m_minShowerSpineLength
The minimum length at which all are considered to be tracks.
static MvaFeatureVector ConcatenateFeatureLists()
Recursively concatenate vectors of features (terminating method)
Definition: LArMvaHelper.h:541
void GetShowerLikeClusterEndPoints(const pandora::ClusterList &clusterList, pandora::ClusterList &showerLikeClusters, ClusterEndPointsMap &clusterEndPointsMap) const
Add the endpoints of any shower-like clusters to the map.
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
Definition: ModuleGraph.h:25
std::map< pandora::HitType, const std::reference_wrapper< HitKDTree2D > > KDTreeMap
Map array of hit kd trees for passing to tools.
float GetBeamDeweightingScore(const BeamConstants &beamConstants, const pandora::Vertex *const pVertex) const
Get the beam deweighting score for a vertex.
static void GetTrueTestBeamParticles(const pandora::MCParticleList *const pMCParticleList, pandora::MCParticleVector &trueTestBeamParticles)
Get triggered test beam MC particles from an input MC particle list.
Header file for the energy kick feature tool class.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_showerClusteringDistance
The shower clustering distance.
float m_axisHits
The hit density along the axis between the two vertices.
static bool IsBeamNeutrinoFinalState(const pandora::MCParticle *const pMCParticle)
Returns true if passed a primary neutrino final state MCParticle.
KDTreeBox build_2d_kd_search_region(const pandora::CaloHit *const point, const float x_span, const float z_span)
build_2d_kd_search_region
pandora::StringVector m_inputClusterListNames
The list of cluster list names.
void search(const KDTreeBoxT< DIM > &searchBox, std::vector< KDTreeNodeInfoT< DATA, DIM >> &resRecHitList)
Search in the KDTree for all points that would be contained in the given searchbox The founded points...
std::map< const pandora::Vertex *const, VertexFeatureInfo > VertexFeatureInfoMap
float m_rPhiFineTuningRadius
The maximum distance the r/phi tune can move a vertex.
void GetBestVertex(const pandora::VertexVector &vertexVector, const pandora::Vertex *&pBestVertex, float &bestVertexDr) const
Use the MC information to get the best vertex from a list.
Header file for the r/phi feature tool class.
Header file for the global asymmetry feature tool class.
std::string GetInteractionType() const
Get the interaction type string.