LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
CosmicRayTaggingTool.cc
Go to the documentation of this file.
1 
9 #include "Pandora/AlgorithmHeaders.h"
10 
12 
15 
17 
18 using namespace pandora;
19 
20 namespace lar_content
21 {
22 
23 CosmicRayTaggingTool::CosmicRayTaggingTool() :
24  m_cutMode("nominal"),
25  m_angularUncertainty(5.f),
26  m_positionalUncertainty(3.f),
27  m_maxAssociationDist(3.f * 18.f),
28  m_minimumHits(15),
29  m_inTimeMargin(5.f),
30  m_inTimeMaxX0(1.f),
31  m_marginY(20.f),
32  m_marginZ(10.f),
33  m_maxNeutrinoCosTheta(0.2f),
34  m_minCosmicCosTheta(0.6f),
35  m_maxCosmicCurvature(0.04f),
36  m_face_Xa(0.f),
37  m_face_Xc(0.f),
38  m_face_Yb(0.f),
39  m_face_Yt(0.f),
40  m_face_Zu(0.f),
41  m_face_Zd(0.f)
42 {
43 }
44 
45 //------------------------------------------------------------------------------------------------------------------------------------------
46 
48 {
49  if ("nominal" == m_cutMode)
50  {
51  // Nominal values set in constructor
52  }
53  else if ("cautious" == m_cutMode)
54  {
55  m_minCosmicCosTheta = std::numeric_limits<double>::max();
56  m_maxCosmicCurvature = -std::numeric_limits<double>::max();
57  }
58  else if ("aggressive" == m_cutMode)
59  {
60  m_minCosmicCosTheta = 0.6f;
61  m_maxCosmicCurvature = 0.1f;
62  }
63  else
64  {
65  std::cout << "CosmicRayTaggingTool - invalid cut mode " << m_cutMode << std::endl;
66  return STATUS_CODE_INVALID_PARAMETER;
67  }
68 
69  return STATUS_CODE_SUCCESS;
70 }
71 
72 //------------------------------------------------------------------------------------------------------------------------------------------
73 
74 void CosmicRayTaggingTool::FindAmbiguousPfos(const PfoList &parentCosmicRayPfos, PfoList &ambiguousPfos, const MasterAlgorithm *const /*pAlgorithm*/)
75 {
76  if (this->GetPandora().GetSettings()->ShouldDisplayAlgorithmInfo())
77  std::cout << "----> Running Algorithm Tool: " << this->GetInstanceName() << ", " << this->GetType() << std::endl;
78 
79  // TODO First time only, TODO Refactor with master algorithm
80  const LArTPCMap &larTPCMap(this->GetPandora().GetGeometry()->GetLArTPCMap());
81  const LArTPC *const pFirstLArTPC(larTPCMap.begin()->second);
82 
83  float parentMinX(pFirstLArTPC->GetCenterX() - 0.5f * pFirstLArTPC->GetWidthX());
84  float parentMaxX(pFirstLArTPC->GetCenterX() + 0.5f * pFirstLArTPC->GetWidthX());
85  float parentMinY(pFirstLArTPC->GetCenterY() - 0.5f * pFirstLArTPC->GetWidthY());
86  float parentMaxY(pFirstLArTPC->GetCenterY() + 0.5f * pFirstLArTPC->GetWidthY());
87  float parentMinZ(pFirstLArTPC->GetCenterZ() - 0.5f * pFirstLArTPC->GetWidthZ());
88  float parentMaxZ(pFirstLArTPC->GetCenterZ() + 0.5f * pFirstLArTPC->GetWidthZ());
89 
90  for (const LArTPCMap::value_type &mapEntry : larTPCMap)
91  {
92  const LArTPC *const pLArTPC(mapEntry.second);
93  parentMinX = std::min(parentMinX, pLArTPC->GetCenterX() - 0.5f * pLArTPC->GetWidthX());
94  parentMaxX = std::max(parentMaxX, pLArTPC->GetCenterX() + 0.5f * pLArTPC->GetWidthX());
95  parentMinY = std::min(parentMinY, pLArTPC->GetCenterY() - 0.5f * pLArTPC->GetWidthY());
96  parentMaxY = std::max(parentMaxY, pLArTPC->GetCenterY() + 0.5f * pLArTPC->GetWidthY());
97  parentMinZ = std::min(parentMinZ, pLArTPC->GetCenterZ() - 0.5f * pLArTPC->GetWidthZ());
98  parentMaxZ = std::max(parentMaxZ, pLArTPC->GetCenterZ() + 0.5f * pLArTPC->GetWidthZ());
99  }
100 
101  m_face_Xa = parentMinX;
102  m_face_Xc = parentMaxX;
103  m_face_Yb = parentMinY;
104  m_face_Yt = parentMaxY;
105  m_face_Zu = parentMinZ;
106  m_face_Zd = parentMaxZ;
107 
108  PfoToPfoListMap pfoAssociationMap;
109  this->GetPfoAssociations(parentCosmicRayPfos, pfoAssociationMap);
110 
111  PfoToSliceIdMap pfoToSliceIdMap;
112  this->SliceEvent(parentCosmicRayPfos, pfoAssociationMap, pfoToSliceIdMap);
113 
114  CRCandidateList candidates;
115  this->GetCRCandidates(parentCosmicRayPfos, pfoToSliceIdMap, candidates);
116 
117  PfoToBoolMap pfoToInTimeMap;
118  this->CheckIfInTime(candidates, pfoToInTimeMap);
119 
120  PfoToBoolMap pfoToIsContainedMap;
121  this->CheckIfContained(candidates, pfoToIsContainedMap);
122 
123  PfoToBoolMap pfoToIsTopToBottomMap;
124  this->CheckIfTopToBottom(candidates, pfoToIsTopToBottomMap);
125 
126  UIntSet neutrinoSliceSet;
127  this->GetNeutrinoSlices(candidates, pfoToInTimeMap, pfoToIsContainedMap, neutrinoSliceSet);
128 
129  PfoToBoolMap pfoToIsLikelyCRMuonMap;
130  this->TagCRMuons(candidates, pfoToInTimeMap, pfoToIsTopToBottomMap, neutrinoSliceSet, pfoToIsLikelyCRMuonMap);
131 
132  for (const ParticleFlowObject *const pPfo : parentCosmicRayPfos)
133  {
134  if (!pfoToIsLikelyCRMuonMap.at(pPfo))
135  ambiguousPfos.push_back(pPfo);
136  }
137 }
138 
139 //------------------------------------------------------------------------------------------------------------------------------------------
140 
141 bool CosmicRayTaggingTool::GetValid3DCluster(const ParticleFlowObject *const pPfo, const Cluster *&pCluster3D) const
142 {
143  pCluster3D = nullptr;
144  ClusterList clusters3D;
145  LArPfoHelper::GetThreeDClusterList(pPfo, clusters3D);
146 
147  // ATTN Only uses first cluster with hits of type TPC_3D
148  if (clusters3D.empty() || (clusters3D.front()->GetNCaloHits() < m_minimumHits))
149  return false;
150 
151  pCluster3D = clusters3D.front();
152  return true;
153 }
154 
155 //------------------------------------------------------------------------------------------------------------------------------------------
156 
157 void CosmicRayTaggingTool::GetPfoAssociations(const PfoList &parentCosmicRayPfos, PfoToPfoListMap &pfoAssociationMap) const
158 {
159  // ATTN If wire w pitches vary between TPCs, exception will be raised in initialisation of lar pseudolayer plugin
160  const LArTPC *const pFirstLArTPC(this->GetPandora().GetGeometry()->GetLArTPCMap().begin()->second);
161  const float layerPitch(pFirstLArTPC->GetWirePitchW());
162 
163  PfoToSlidingFitsMap pfoToSlidingFitsMap;
164 
165  for (const ParticleFlowObject *const pPfo : parentCosmicRayPfos)
166  {
167  const pandora::Cluster *pCluster(nullptr);
168  if (!this->GetValid3DCluster(pPfo, pCluster) || !pCluster)
169  continue;
170 
171  (void)pfoToSlidingFitsMap.insert(PfoToSlidingFitsMap::value_type(pPfo,
172  std::make_pair(ThreeDSlidingFitResult(pCluster, 5, layerPitch), ThreeDSlidingFitResult(pCluster, 100, layerPitch)))); // TODO Configurable
173  }
174 
175  for (const ParticleFlowObject *const pPfo1 : parentCosmicRayPfos)
176  {
177  PfoToSlidingFitsMap::const_iterator iter1(pfoToSlidingFitsMap.find(pPfo1));
178  if (pfoToSlidingFitsMap.end() == iter1)
179  continue;
180 
181  const ThreeDSlidingFitResult &fitPos1(iter1->second.first), &fitDir1(iter1->second.second);
182 
183  for (const ParticleFlowObject *const pPfo2 : parentCosmicRayPfos)
184  {
185  if (pPfo1 == pPfo2)
186  continue;
187 
188  PfoToSlidingFitsMap::const_iterator iter2(pfoToSlidingFitsMap.find(pPfo2));
189  if (pfoToSlidingFitsMap.end() == iter2)
190  continue;
191 
192  const ThreeDSlidingFitResult &fitPos2(iter2->second.first), &fitDir2(iter2->second.second);
193 
194  // TODO Use existing LArPointingClusters and IsEmission/IsNode logic, for consistency
195  if (!(this->CheckAssociation(fitPos1.GetGlobalMinLayerPosition(), fitDir1.GetGlobalMinLayerDirection() * -1.f,
196  fitPos2.GetGlobalMinLayerPosition(), fitDir2.GetGlobalMinLayerDirection() * -1.f) ||
197  this->CheckAssociation(fitPos1.GetGlobalMinLayerPosition(), fitDir1.GetGlobalMinLayerDirection() * -1.f,
198  fitPos2.GetGlobalMaxLayerPosition(), fitDir2.GetGlobalMaxLayerDirection()) ||
199  this->CheckAssociation(fitPos1.GetGlobalMaxLayerPosition(), fitDir1.GetGlobalMaxLayerDirection(),
200  fitPos2.GetGlobalMinLayerPosition(), fitDir2.GetGlobalMinLayerDirection() * -1.f) ||
201  this->CheckAssociation(fitPos1.GetGlobalMaxLayerPosition(), fitDir1.GetGlobalMaxLayerDirection(),
202  fitPos2.GetGlobalMaxLayerPosition(), fitDir2.GetGlobalMaxLayerDirection())))
203  {
204  continue;
205  }
206 
207  PfoList &pfoList1(pfoAssociationMap[pPfo1]), &pfoList2(pfoAssociationMap[pPfo2]);
208 
209  if (pfoList1.end() == std::find(pfoList1.begin(), pfoList1.end(), pPfo2))
210  pfoList1.push_back(pPfo2);
211 
212  if (pfoList2.end() == std::find(pfoList2.begin(), pfoList2.end(), pPfo1))
213  pfoList2.push_back(pPfo1);
214  }
215  }
216 }
217 
218 //------------------------------------------------------------------------------------------------------------------------------------------
219 
221  const CartesianVector &endPoint1, const CartesianVector &endDir1, const CartesianVector &endPoint2, const CartesianVector &endDir2) const
222 {
223  // TODO This function needs significant tidying and checks (variable names must be self describing, check deltaTheta, etc.)
224  const CartesianVector n(endDir1.GetUnitVector());
225  const CartesianVector m(endDir2.GetUnitVector());
226  const CartesianVector a(endPoint2 - endPoint1);
227  const float b(n.GetDotProduct(m));
228 
229  // Parallel lines never meet
230  if (std::fabs(b - 1.f) < std::numeric_limits<float>::epsilon())
231  return false;
232 
233  // Distance from endPoint1 along endDir1 to the point of closest approach
234  const float lambda((n - m * b).GetDotProduct(a) / (1.f - b * b));
235 
236  // Distance from endPoint2 along endDir2 to the point of closest approach
237  const float mu((n * b - m).GetDotProduct(a) / (1.f - b * b));
238 
239  // Calculate the maximum "vertex uncertainty"
240  const float deltaTheta(m_angularUncertainty * M_PI / 180.f);
241  const float maxVertexUncertainty(m_maxAssociationDist * std::sin(deltaTheta) + m_positionalUncertainty);
242 
243  // Ensure that the distances to the point of closest approch are within the limits
244  if ((lambda < -maxVertexUncertainty) || (mu < -maxVertexUncertainty) || (lambda > m_maxAssociationDist + maxVertexUncertainty) ||
245  (mu > m_maxAssociationDist + maxVertexUncertainty))
246  {
247  return false;
248  }
249 
250  // Ensure point of closest approach is within detector
251  const CartesianVector impactPosition((endPoint1 + n * lambda + endPoint2 + m * mu) * 0.5f);
252 
253  if ((impactPosition.GetX() < m_face_Xa - maxVertexUncertainty) || (impactPosition.GetX() > m_face_Xc + maxVertexUncertainty) ||
254  (impactPosition.GetY() < m_face_Yb - maxVertexUncertainty) || (impactPosition.GetY() > m_face_Yt + maxVertexUncertainty) ||
255  (impactPosition.GetZ() < m_face_Zu - maxVertexUncertainty) || (impactPosition.GetZ() > m_face_Zd + maxVertexUncertainty))
256  {
257  return false;
258  }
259 
260  // Compare distance of closest approach and max uncertainty in impact parameter
261  const float maxImpactDist(std::sin(deltaTheta) * (std::fabs(mu) + std::fabs(lambda)) + m_positionalUncertainty);
262  const CartesianVector d(a - n * lambda + m * mu);
263 
264  return (d.GetMagnitude() < maxImpactDist);
265 }
266 
267 //------------------------------------------------------------------------------------------------------------------------------------------
268 
269 void CosmicRayTaggingTool::SliceEvent(const PfoList &parentCosmicRayPfos, const PfoToPfoListMap &pfoAssociationMap, PfoToSliceIdMap &pfoToSliceIdMap) const
270 {
271  SliceList sliceList;
272 
273  for (const ParticleFlowObject *const pPfo : parentCosmicRayPfos)
274  {
275  bool isAlreadyInSlice(false);
276 
277  for (const PfoList &slice : sliceList)
278  {
279  if (std::find(slice.begin(), slice.end(), pPfo) != slice.end())
280  {
281  isAlreadyInSlice = true;
282  break;
283  }
284  }
285 
286  if (!isAlreadyInSlice)
287  {
288  sliceList.push_back(PfoList());
289  this->FillSlice(pPfo, pfoAssociationMap, sliceList.back());
290  }
291  }
292 
293  unsigned int sliceId(0);
294  for (const PfoList &slice : sliceList)
295  {
296  for (const ParticleFlowObject *const pPfo : slice)
297  {
298  if (!pfoToSliceIdMap.insert(PfoToSliceIdMap::value_type(pPfo, sliceId)).second)
299  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
300  }
301 
302  ++sliceId;
303  }
304 }
305 
306 //------------------------------------------------------------------------------------------------------------------------------------------
307 
308 void CosmicRayTaggingTool::FillSlice(const ParticleFlowObject *const pPfo, const PfoToPfoListMap &pfoAssociationMap, PfoList &slice) const
309 {
310  if (std::find(slice.begin(), slice.end(), pPfo) != slice.end())
311  return;
312 
313  slice.push_back(pPfo);
314 
315  PfoToPfoListMap::const_iterator iter(pfoAssociationMap.find(pPfo));
316 
317  if (pfoAssociationMap.end() != iter)
318  {
319  for (const ParticleFlowObject *const pAssociatedPfo : iter->second)
320  this->FillSlice(pAssociatedPfo, pfoAssociationMap, slice);
321  }
322 }
323 
324 //------------------------------------------------------------------------------------------------------------------------------------------
325 
326 void CosmicRayTaggingTool::GetCRCandidates(const PfoList &parentCosmicRayPfos, const PfoToSliceIdMap &pfoToSliceIdMap, CRCandidateList &candidates) const
327 {
328  for (const ParticleFlowObject *const pPfo : parentCosmicRayPfos)
329  {
330  if (!LArPfoHelper::IsFinalState(pPfo))
331  throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
332 
333  candidates.push_back(CRCandidate(this->GetPandora(), pPfo, pfoToSliceIdMap.at(pPfo)));
334  }
335 }
336 
337 //------------------------------------------------------------------------------------------------------------------------------------------
338 
339 void CosmicRayTaggingTool::CheckIfInTime(const CRCandidateList &candidates, PfoToBoolMap &pfoToInTimeMap) const
340 {
341  const LArTPCMap &larTPCMap(this->GetPandora().GetGeometry()->GetLArTPCMap());
342 
343  for (const CRCandidate &candidate : candidates)
344  {
345  // Cosmic-ray muons extending outside of (single) physical volume if given t0 is that of the beam particle
346  float minX(std::numeric_limits<float>::max()), maxX(-std::numeric_limits<float>::max());
347 
348  if (candidate.m_canFit)
349  {
350  minX = ((candidate.m_endPoint1.GetX() < candidate.m_endPoint2.GetX()) ? candidate.m_endPoint1.GetX() : candidate.m_endPoint2.GetX());
351  maxX = ((candidate.m_endPoint1.GetX() > candidate.m_endPoint2.GetX()) ? candidate.m_endPoint1.GetX() : candidate.m_endPoint2.GetX());
352  }
353  else
354  {
355  // Handle any particles with small numbers of 3D hits, for which no 3D sliding fit information is available
356  for (const Cluster *const pCluster : candidate.m_pPfo->GetClusterList())
357  {
358  float clusterMinX(std::numeric_limits<float>::max()), clusterMaxX(-std::numeric_limits<float>::max());
359  pCluster->GetClusterSpanX(clusterMinX, clusterMaxX);
360  minX = std::min(clusterMinX, minX);
361  maxX = std::max(clusterMaxX, maxX);
362  }
363  }
364 
365  bool isInTime((minX > m_face_Xa - m_inTimeMargin) && (maxX < m_face_Xc + m_inTimeMargin));
366 
367  // Cosmic-ray muons that have been shifted and stitched across mid plane between volumes
368  if (isInTime)
369  {
370  try
371  {
372  if (std::fabs(LArPfoHelper::GetVertex(candidate.m_pPfo)->GetX0()) > m_inTimeMaxX0)
373  isInTime = false;
374  }
375  catch (const StatusCodeException &)
376  {
377  }
378  }
379 
380  // Cosmic-ray muons extending outside of (any individual) physical volume if given t0 is that of the beam particle
381  if (isInTime)
382  {
383  CaloHitList caloHitList;
384  LArPfoHelper::GetCaloHits(candidate.m_pPfo, TPC_VIEW_U, caloHitList);
385  LArPfoHelper::GetCaloHits(candidate.m_pPfo, TPC_VIEW_V, caloHitList);
386  LArPfoHelper::GetCaloHits(candidate.m_pPfo, TPC_VIEW_W, caloHitList);
387 
388  bool isFirstHit(true);
389  bool isInSingleVolume(true);
390  unsigned int volumeId(std::numeric_limits<unsigned int>::max());
391 
392  for (const CaloHit *const pCaloHit : caloHitList)
393  {
394  const LArCaloHit *const pLArCaloHit(dynamic_cast<const LArCaloHit *>(pCaloHit));
395 
396  if (!pLArCaloHit)
397  continue;
398 
399  if (isFirstHit)
400  {
401  isFirstHit = false;
402  volumeId = pLArCaloHit->GetLArTPCVolumeId();
403  }
404  else if (volumeId != pLArCaloHit->GetLArTPCVolumeId())
405  {
406  isInSingleVolume = false;
407  break;
408  }
409  }
410 
411  LArTPCMap::const_iterator tpcIter(larTPCMap.find(volumeId));
412 
413  if (isInSingleVolume && (larTPCMap.end() != tpcIter))
414  {
415  const float thisFaceXLow(tpcIter->second->GetCenterX() - 0.5f * tpcIter->second->GetWidthX());
416  const float thisFaceXHigh(tpcIter->second->GetCenterX() + 0.5f * tpcIter->second->GetWidthX());
417 
418  if (!((minX > thisFaceXLow - m_inTimeMargin) && (maxX < thisFaceXHigh + m_inTimeMargin)))
419  isInTime = false;
420  }
421  }
422 
423  if (!pfoToInTimeMap.insert(PfoToBoolMap::value_type(candidate.m_pPfo, isInTime)).second)
424  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
425  }
426 }
427 
428 //------------------------------------------------------------------------------------------------------------------------------------------
429 
430 void CosmicRayTaggingTool::CheckIfContained(const CRCandidateList &candidates, PfoToBoolMap &pfoToIsContainedMap) const
431 {
432  for (const CRCandidate &candidate : candidates)
433  {
434  const float upperY(
435  (candidate.m_endPoint1.GetY() > candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetY() : candidate.m_endPoint2.GetY());
436  const float lowerY(
437  (candidate.m_endPoint1.GetY() < candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetY() : candidate.m_endPoint2.GetY());
438 
439  const float zAtUpperY(
440  (candidate.m_endPoint1.GetY() > candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetZ() : candidate.m_endPoint2.GetZ());
441  const float zAtLowerY(
442  (candidate.m_endPoint1.GetY() < candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetZ() : candidate.m_endPoint2.GetZ());
443 
444  const bool isContained((upperY < m_face_Yt - m_marginY) && (upperY > m_face_Yb + m_marginY) && (lowerY < m_face_Yt - m_marginY) &&
445  (lowerY > m_face_Yb + m_marginY) && (zAtUpperY < m_face_Zd - m_marginZ) && (zAtUpperY > m_face_Zu + m_marginZ) &&
446  (zAtLowerY < m_face_Zd - m_marginZ) && (zAtLowerY > m_face_Zu + m_marginZ));
447 
448  if (!pfoToIsContainedMap.insert(PfoToBoolMap::value_type(candidate.m_pPfo, isContained)).second)
449  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
450  }
451 }
452 
453 //------------------------------------------------------------------------------------------------------------------------------------------
454 
455 void CosmicRayTaggingTool::CheckIfTopToBottom(const CRCandidateList &candidates, PfoToBoolMap &pfoToIsTopToBottomMap) const
456 {
457  for (const CRCandidate &candidate : candidates)
458  {
459  const float upperY(
460  (candidate.m_endPoint1.GetY() > candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetY() : candidate.m_endPoint2.GetY());
461  const float lowerY(
462  (candidate.m_endPoint1.GetY() < candidate.m_endPoint2.GetY()) ? candidate.m_endPoint1.GetY() : candidate.m_endPoint2.GetY());
463 
464  const bool isTopToBottom((upperY > m_face_Yt - m_marginY) && (lowerY < m_face_Yb + m_marginY));
465 
466  if (!pfoToIsTopToBottomMap.insert(PfoToBoolMap::value_type(candidate.m_pPfo, isTopToBottom)).second)
467  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
468  }
469 }
470 
471 //------------------------------------------------------------------------------------------------------------------------------------------
472 
473 void CosmicRayTaggingTool::GetNeutrinoSlices(const CRCandidateList &candidates, const PfoToBoolMap &pfoToInTimeMap,
474  const PfoToBoolMap &pfoToIsContainedMap, UIntSet &neutrinoSliceSet) const
475 {
476  IntBoolMap sliceIdToIsInTimeMap;
477 
478  for (const CRCandidate &candidate : candidates)
479  {
480  if (sliceIdToIsInTimeMap.find(candidate.m_sliceId) == sliceIdToIsInTimeMap.end())
481  sliceIdToIsInTimeMap.insert(std::make_pair(candidate.m_sliceId, true));
482 
483  if (!pfoToInTimeMap.at(candidate.m_pPfo))
484  sliceIdToIsInTimeMap.at(candidate.m_sliceId) = false;
485  }
486 
487  for (const CRCandidate &candidate : candidates)
488  {
489  if (neutrinoSliceSet.count(candidate.m_sliceId))
490  continue;
491 
492  const bool likelyNeutrino(candidate.m_canFit && sliceIdToIsInTimeMap.at(candidate.m_sliceId) &&
493  (candidate.m_theta < m_maxNeutrinoCosTheta || pfoToIsContainedMap.at(candidate.m_pPfo)));
494 
495  if (likelyNeutrino)
496  (void)neutrinoSliceSet.insert(candidate.m_sliceId);
497  }
498 }
499 
500 //------------------------------------------------------------------------------------------------------------------------------------------
501 
502 void CosmicRayTaggingTool::TagCRMuons(const CRCandidateList &candidates, const PfoToBoolMap &pfoToInTimeMap,
503  const PfoToBoolMap &pfoToIsTopToBottomMap, const UIntSet &neutrinoSliceSet, PfoToBoolMap &pfoToIsLikelyCRMuonMap) const
504 {
505  for (const CRCandidate &candidate : candidates)
506  {
507  const bool likelyCRMuon(!neutrinoSliceSet.count(candidate.m_sliceId) &&
508  (!pfoToInTimeMap.at(candidate.m_pPfo) ||
509  (candidate.m_canFit &&
510  (pfoToIsTopToBottomMap.at(candidate.m_pPfo) ||
511  ((candidate.m_theta > m_minCosmicCosTheta) && (candidate.m_curvature < m_maxCosmicCurvature))))));
512 
513  if (!pfoToIsLikelyCRMuonMap.insert(PfoToBoolMap::value_type(candidate.m_pPfo, likelyCRMuon)).second)
514  throw StatusCodeException(STATUS_CODE_ALREADY_PRESENT);
515  }
516 }
517 
518 //------------------------------------------------------------------------------------------------------------------------------------------
519 //------------------------------------------------------------------------------------------------------------------------------------------
520 
521 CosmicRayTaggingTool::CRCandidate::CRCandidate(const Pandora &pandora, const ParticleFlowObject *const pPfo, const unsigned int sliceId) :
522  m_pPfo(pPfo),
523  m_sliceId(sliceId),
524  m_canFit(false),
525  m_endPoint1(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()),
526  m_endPoint2(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(), std::numeric_limits<float>::max()),
527  m_length(std::numeric_limits<float>::max()),
528  m_curvature(std::numeric_limits<float>::max()),
529  m_theta(std::numeric_limits<float>::max())
530 {
531  ClusterList clusters3D;
532  LArPfoHelper::GetThreeDClusterList(pPfo, clusters3D);
533 
534  if (!clusters3D.empty() && (clusters3D.front()->GetNCaloHits() > 15)) // TODO Configurable
535  {
536  m_canFit = true;
537  const LArTPC *const pFirstLArTPC(pandora.GetGeometry()->GetLArTPCMap().begin()->second);
538  const ThreeDSlidingFitResult slidingFitResult(clusters3D.front(), 5, pFirstLArTPC->GetWirePitchW()); // TODO Configurable
539  this->CalculateFitVariables(slidingFitResult);
540  }
541 }
542 
543 //------------------------------------------------------------------------------------------------------------------------------------------
544 
546 {
547  m_endPoint1 = slidingFitResult.GetGlobalMinLayerPosition();
548  m_endPoint2 = slidingFitResult.GetGlobalMaxLayerPosition();
549  m_length = (m_endPoint2 - m_endPoint1).GetMagnitude();
550 
551  if (std::fabs(m_length) > std::numeric_limits<float>::epsilon())
552  m_theta = std::fabs(m_endPoint2.GetY() - m_endPoint1.GetY()) / m_length;
553 
554  const float layerPitch(slidingFitResult.GetFirstFitResult().GetLayerPitch());
555 
556  CartesianPointVector directionList;
557  for (int i = slidingFitResult.GetMinLayer(); i < slidingFitResult.GetMaxLayer(); ++i)
558  {
559  CartesianVector direction(0.f, 0.f, 0.f);
560  if (STATUS_CODE_SUCCESS == slidingFitResult.GetGlobalFitDirection(static_cast<float>(i) * layerPitch, direction))
561  directionList.push_back(direction);
562  }
563 
564  CartesianVector meanDirection(0.f, 0.f, 0.f);
565  for (const CartesianVector &direction : directionList)
566  meanDirection += direction;
567 
568  if (!directionList.empty() > 0)
569  meanDirection *= 1.f / static_cast<float>(directionList.size());
570 
571  m_curvature = 0.f;
572  for (const CartesianVector &direction : directionList)
573  m_curvature += (direction - meanDirection).GetMagnitude();
574 
575  if (!directionList.empty() > 0)
576  m_curvature *= 1.f / static_cast<float>(directionList.size());
577 }
578 
579 //------------------------------------------------------------------------------------------------------------------------------------------
580 
581 StatusCode CosmicRayTaggingTool::ReadSettings(const TiXmlHandle xmlHandle)
582 {
583  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CutMode", m_cutMode));
584  std::transform(m_cutMode.begin(), m_cutMode.end(), m_cutMode.begin(), ::tolower);
585 
586  PANDORA_RETURN_RESULT_IF_AND_IF(
587  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "AngularUncertainty", m_angularUncertainty));
588 
589  PANDORA_RETURN_RESULT_IF_AND_IF(
590  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PositionalUncertainty", m_positionalUncertainty));
591 
592  PANDORA_RETURN_RESULT_IF_AND_IF(
593  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxAssociationDist", m_maxAssociationDist));
594 
595  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HitThreshold", m_minimumHits));
596 
597  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "InTimeMargin", m_inTimeMargin));
598 
599  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "InTimeMaxX0", m_inTimeMaxX0));
600 
601  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MarginY", m_marginY));
602 
603  PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MarginZ", m_marginZ));
604 
605  PANDORA_RETURN_RESULT_IF_AND_IF(
606  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxNeutrinoCosTheta", m_maxNeutrinoCosTheta));
607 
608  PANDORA_RETURN_RESULT_IF_AND_IF(
609  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosmicCosTheta", m_minCosmicCosTheta));
610 
611  PANDORA_RETURN_RESULT_IF_AND_IF(
612  STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCosmicCurvature", m_maxCosmicCurvature));
613 
614  return STATUS_CODE_SUCCESS;
615 }
616 
617 } // namespace lar_content
std::unordered_map< const pandora::ParticleFlowObject *, pandora::PfoList > PfoToPfoListMap
std::unordered_map< const pandora::ParticleFlowObject *, bool > PfoToBoolMap
void FindAmbiguousPfos(const pandora::PfoList &parentCosmicRayPfos, pandora::PfoList &ambiguousPfos, const MasterAlgorithm *const pAlgorithm)
Find the list of ambiguous pfos (could represent cosmic-ray muons or neutrinos)
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void TagCRMuons(const CRCandidateList &candidates, const PfoToBoolMap &pfoToInTimeMap, const PfoToBoolMap &pfoToIsTopToBottomMap, const UIntSet &neutrinoSliceSet, PfoToBoolMap &pfoToIsLikelyCRMuonMap) const
Tag Pfos which are likely to be a CR muon.
Header file for the pfo helper class.
std::unordered_map< const pandora::ParticleFlowObject *, unsigned int > PfoToSliceIdMap
float m_maxNeutrinoCosTheta
The maximum cos(theta) that a Pfo can have to be classified as a likely neutrino. ...
const pandora::CartesianVector & GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
Header file for the lar calo hit class.
void SliceEvent(const pandora::PfoList &parentCosmicRayPfos, const PfoToPfoListMap &pfoAssociationMap, PfoToSliceIdMap &pfoToSliceIdMap) const
Break the event up into slices of associated Pfos.
bool CheckAssociation(const pandora::CartesianVector &endPoint1, const pandora::CartesianVector &endDir1, const pandora::CartesianVector &endPoint2, const pandora::CartesianVector &endDir2) const
Check whethe two Pfo endpoints are associated by distance of closest approach.
float m_positionalUncertainty
The uncertainty in cm for the position of Pfo endpoint in 3D.
void CheckIfTopToBottom(const CRCandidateList &candidates, PfoToBoolMap &pfoToIsTopToBottomMap) const
Check if each candidate is "top to bottom".
int GetMaxLayer() const
Get the maximum occupied layer in the sliding fit.
static const pandora::Vertex * GetVertex(const pandora::ParticleFlowObject *const pPfo)
Get the pfo vertex.
double m_theta
Direction made with vertical.
bool m_canFit
If there are a sufficient number of 3D hits to perform a fitting.
Class to encapsulate the logic required determine if a Pfo should or shouldn&#39;t be tagged as a cosmic ...
STL namespace.
intermediate_table::const_iterator const_iterator
void CheckIfContained(const CRCandidateList &candidates, PfoToBoolMap &pfoToIsContainedMap) const
Check if each candidate is "contained" (contained = no associations to Y or Z detector faces...
unsigned int GetLArTPCVolumeId() const
Get the lar tpc volume id.
Definition: LArCaloHit.h:167
float m_angularUncertainty
The uncertainty in degrees for the angle of a Pfo.
LAr calo hit class.
Definition: LArCaloHit.h:39
CRCandidate(const pandora::Pandora &pandora, const pandora::ParticleFlowObject *const pPfo, const unsigned int sliceId)
Constructor.
void FillSlice(const pandora::ParticleFlowObject *const pPfo, const PfoToPfoListMap &pfoAssociationMap, pandora::PfoList &slice) const
Fill a slice iteratively using Pfo associations.
std::string m_cutMode
Choose a set of cuts using a keyword - "cautious" = remove as few neutrinos as possible "nominal" = o...
std::list< CRCandidate > CRCandidateList
TFile f
Definition: plotHisto.C:6
double m_length
Straight line length of the linear fit.
float m_inTimeMaxX0
The maximum pfo x0 (determined from shifted vertex) to allow pfo to still be considered in time...
Header file for the cluster helper class.
void CheckIfInTime(const CRCandidateList &candidates, PfoToBoolMap &pfoToInTimeMap) const
Check if each candidate is "in time".
pandora::CartesianVector m_endPoint2
Second fitted end point in 3D.
const pandora::CartesianVector & GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
float m_inTimeMargin
The maximum distance outside of the physical detector volume that a Pfo may be to still be considered...
double m_curvature
Measure of the curvature of the track.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
Float_t d
Definition: plot.C:235
void CalculateFitVariables(const ThreeDSlidingFitResult &slidingFitResult)
Calculate all variables which require a fit.
float GetLayerPitch() const
Get the layer pitch, units cm.
void GetPfoAssociations(const pandora::PfoList &parentCosmicRayPfos, PfoToPfoListMap &pfoAssociationMap) const
Get mapping between Pfos that are associated with it other by pointing.
int GetMinLayer() const
Get the minimum occupied layer in the sliding fit.
static void GetThreeDClusterList(const pandora::ParticleFlowObject *const pPfo, pandora::ClusterList &clusterList)
Get the list of 3D clusters from an input pfo.
const pandora::CartesianVector & GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
float m_minCosmicCosTheta
The minimum cos(theta) that a Pfo can have to be classified as a likely CR muon.
bool GetValid3DCluster(const pandora::ParticleFlowObject *const pPfo, const pandora::Cluster *&pCluster3D) const
Get the 3D calo hit cluster associated with a given Pfo, and check if it has sufficient hits...
std::vector< pandora::PfoList > SliceList
std::unordered_map< int, bool > IntBoolMap
float m_marginZ
The minimum distance from a detector Z-face for a Pfo to be associated.
float m_marginY
The minimum distance from a detector Y-face for a Pfo to be associated.
float m_maxCosmicCurvature
The maximum curvature that a Pfo can have to be classified as a likely CR muon.
MasterAlgorithm class.
static bool IsFinalState(const pandora::ParticleFlowObject *const pPfo)
Whether a pfo is a primary parent particle.
const TwoDSlidingFitResult & GetFirstFitResult() const
Get the first sliding fit result for this cluster.
Header file for the cosmic-ray tagging tool class.
float m_maxAssociationDist
The maximum distance from endpoint to point of closest approach, typically a multiple of LAr radiatio...
Char_t n[5]
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
unsigned int m_minimumHits
The minimum number of hits for a Pfo to be considered.
std::unordered_map< const pandora::ParticleFlowObject *, SlidingFitPair > PfoToSlidingFitsMap
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:82
void GetCRCandidates(const pandora::PfoList &parentCosmicRayPfos, const PfoToSliceIdMap &pfoToSliceIdMap, CRCandidateList &candidates) const
Make a list of CRCandidates.
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 GetNeutrinoSlices(const CRCandidateList &candidates, const PfoToBoolMap &pfoToInTimeMap, const PfoToBoolMap &pfoToIsContainedMap, UIntSet &neutrinoSliceSet) const
Get the slice indices which contain a likely neutrino Pfo.
const pandora::CartesianVector & GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
pandora::CartesianVector m_endPoint1
First fitted end point in 3D.