9 #include "Helpers/MCParticleHelper.h" 15 #include "Pandora/PdgTable.h" 17 #include "Objects/CaloHit.h" 18 #include "Objects/ParticleFlowObject.h" 27 m_maxBremsstrahlungSeparation(2.5
f)
40 const LArMCParticle *
const pLArMCParticle(dynamic_cast<const LArMCParticle *>(pLeadingParticle));
43 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
76 throw StatusCodeException(STATUS_CODE_FAILURE);
79 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
81 const MCParticle *pParentMCParticle = pMCParticle;
85 const MCParticleList &parentList(pParentMCParticle->GetParentList());
87 if (parentList.size() != 1)
88 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
90 pParentMCParticle = parentList.front();
93 return pParentMCParticle;
100 for (
const MCParticle *
const pMCParticle : *pMCParticleList)
108 if (pMCParticle == pParentMCParticle)
110 mcToLeadingMap[pMCParticle] = pMCParticle;
115 mcToLeadingMap[pMCParticle] = pLeadingMCParticle;
130 CaloHitList selectedCaloHitList;
155 CaloHitList &selectedCaloHitList,
const bool selectInputHits,
const float minHitSharingFraction,
const CaloHitList &recoMuonHitList,
158 if (!selectInputHits)
160 selectedCaloHitList.insert(selectedCaloHitList.end(), pCaloHitList->begin(), pCaloHitList->end());
164 for (
const CaloHit *
const pCaloHit : *pCaloHitList)
168 const MCParticle *
const pHitParticle(MCParticleHelper::GetMainMCParticle(pCaloHit));
170 if (mcToTargetMCMap.find(pHitParticle) == mcToTargetMCMap.end())
176 if (std::find(recoMuonHitList.begin(), recoMuonHitList.end(), pCaloHit) != recoMuonHitList.end())
181 for (
const auto &mapEntry : pCaloHit->GetMCParticleWeightMap())
182 mcParticleContributionVector.push_back(mapEntry.first);
184 std::sort(mcParticleContributionVector.begin(), mcParticleContributionVector.end(), PointerLessThan<MCParticle>());
186 MCParticleWeightMap targetWeightMap;
187 for (
const MCParticle *
const pMCParticle : mcParticleContributionVector)
189 const float weight(pCaloHit->GetMCParticleWeightMap().at(pMCParticle));
192 if (mcToTargetMCMap.end() != mcIter)
193 targetWeightMap[mcIter->second] +=
weight;
197 for (
const auto &mapEntry : targetWeightMap)
198 mcTargetContributionVector.push_back(mapEntry.first);
199 std::sort(mcTargetContributionVector.begin(), mcTargetContributionVector.end(), PointerLessThan<MCParticle>());
201 float bestTargetWeight(0.
f), targetWeightSum(0.
f);
203 for (
const MCParticle *
const pTargetMCParticle : mcTargetContributionVector)
205 const float targetWeight(targetWeightMap.at(pTargetMCParticle));
206 targetWeightSum += targetWeight;
208 if (targetWeight > bestTargetWeight)
210 bestTargetWeight = targetWeight;
214 if ((targetWeightSum < std::numeric_limits<float>::epsilon()) || ((bestTargetWeight / targetWeightSum) < minHitSharingFraction))
221 selectedCaloHitList.push_back(pCaloHit);
223 catch (
const StatusCodeException &)
234 const MCParticle *
const pHitMCParticle(MCParticleHelper::GetMainMCParticle(pCaloHit));
236 MCParticleList ancestorMCParticleList;
239 bool isPostBremsstrahlung(
false);
240 const MCParticle *leadingMCParticle(
nullptr);
242 for (
const MCParticle *
const pAncestorMCParticle : ancestorMCParticleList)
247 leadingMCParticle = pAncestorMCParticle;
250 if (pAncestorMCParticle->GetParticleId() == PHOTON)
251 isPostBremsstrahlung =
true;
254 if (isPostBremsstrahlung && leadingMCParticle)
256 leadingMCParticleToPostBremsstrahlungHitList[leadingMCParticle].push_back(pCaloHit);
269 for (
auto &entry : leadingMCParticleToPostBremsstrahlungHitList)
270 leadingMCParticleVector.push_back(entry.first);
272 for (
const MCParticle *
const pLeadingMCParticle : leadingMCParticleVector)
275 if (leadingMCToTrueHitListMap.find(pLeadingMCParticle) == leadingMCToTrueHitListMap.end())
279 maxBremsstrahlungSeparation, leadingMCToTrueHitListMap, TPC_VIEW_U);
281 maxBremsstrahlungSeparation, leadingMCToTrueHitListMap, TPC_VIEW_V);
283 maxBremsstrahlungSeparation, leadingMCToTrueHitListMap, TPC_VIEW_W);
293 CaloHitList leadingViewHitList;
294 for (
const CaloHit *
const pCaloHit : leadingMCToTrueHitListMap.at(pLeadingMCParticle))
296 if (pCaloHit->GetHitType() == tpcView)
297 leadingViewHitList.push_back(pCaloHit);
300 if (leadingViewHitList.empty())
303 CaloHitList postBremsstrahlungViewHitList;
304 for (
const CaloHit *
const pCaloHit : leadingMCParticleToPostBremsstrahlungHitList.at(pLeadingMCParticle))
306 if (pCaloHit->GetHitType() == tpcView)
307 postBremsstrahlungViewHitList.push_back(pCaloHit);
310 if (postBremsstrahlungViewHitList.empty())
313 bool hitsAdded(
false);
319 for (
const CaloHit *
const pPostBremsstrahlungHit : postBremsstrahlungViewHitList)
321 if (std::find(leadingViewHitList.begin(), leadingViewHitList.end(), pPostBremsstrahlungHit) != leadingViewHitList.end())
326 if (separationDistance < maxBremsstrahlungSeparation)
328 leadingViewHitList.push_back(pPostBremsstrahlungHit);
335 CaloHitList &leadingHitList(leadingMCToTrueHitListMap.at(pLeadingMCParticle));
336 for (
const CaloHit *
const pCaloHit : leadingViewHitList)
338 if (std::find(leadingHitList.begin(), leadingHitList.end(), pCaloHit) == leadingHitList.end())
339 leadingHitList.push_back(pCaloHit);
347 for (
const MCParticle *
const pMCParticle : *pMCParticleList)
354 if (pMCParticle == pParentMCParticle)
356 selectedParticles.push_back(pMCParticle);
361 selectedParticles.push_back(pMCParticle);
371 CaloHitList &parentTrackHits, CaloHitList &otherTrackHits, CaloHitList &otherShowerHits)
375 for (
const CaloHit *
const pCaloHit : matchedPfoHitList)
377 const MCParticle *
const pHitParticle(MCParticleHelper::GetMainMCParticle(pCaloHit));
381 (pHitParticle == pParentCosmicRay) ? parentTrackHits.push_back(pCaloHit) : otherTrackHits.push_back(pCaloHit);
387 if (pHitLeadingParticle != pLeadingParticle)
388 otherShowerHits.push_back(pCaloHit);
396 const CaloHitList &cosmicRayPfoHitList,
const CaloHitList &leadingMCHitList, CaloHitList &leadingHitsInParentCosmicRay)
398 for (
const CaloHit *
const pCaloHit : cosmicRayPfoHitList)
400 if (std::find(leadingMCHitList.begin(), leadingMCHitList.end(), pCaloHit) != leadingMCHitList.end())
401 leadingHitsInParentCosmicRay.push_back(pCaloHit);
409 float closestDistance(std::numeric_limits<float>::max());
411 CaloHitList caloHitList;
412 pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
414 for (
const CaloHit *
const pCaloHit : caloHitList)
418 if (distance < closestDistance)
419 closestDistance = distance;
422 return closestDistance;
429 float shortestDistanceSquared(std::numeric_limits<float>::max());
430 const CartesianVector referencePoint(pCaloHit->GetPositionVector());
432 for (
const CartesianVector &testPosition : cartesianPointVector)
434 const float separationSquared((testPosition - referencePoint).GetMagnitudeSquared());
436 if (separationSquared < shortestDistanceSquared)
437 shortestDistanceSquared = separationSquared;
440 return std::sqrt(shortestDistanceSquared);
446 const Cluster *
const pCluster,
const float maxDistanceToCluster,
const float maxDistanceToReferencePoint, CartesianVector &closestPosition)
449 float shortestDistanceSquared(std::numeric_limits<float>::max());
451 for (
const CartesianVector &testPosition : cartesianPointVector)
456 const float separationSquared((testPosition - referencePoint).GetMagnitude());
458 if (separationSquared > maxDistanceToReferencePoint)
461 if (separationSquared < shortestDistanceSquared)
463 shortestDistanceSquared = separationSquared;
464 closestPosition = testPosition;
469 return found ? STATUS_CODE_SUCCESS : STATUS_CODE_NOT_FOUND;
475 CartesianVector &outputPosition1, CartesianVector &outputPosition2)
477 bool distanceFound(
false);
478 float minDistanceSquared(std::numeric_limits<float>::max());
480 CartesianVector closestPosition1(0.
f, 0.
f, 0.
f);
481 CartesianVector closestPosition2(0.
f, 0.
f, 0.
f);
483 CaloHitList caloHitList2;
484 pCluster2->GetOrderedCaloHitList().FillCaloHitList(caloHitList2);
486 for (
const CartesianVector &positionVector1 : cartesianPointVector1)
488 for (
const CaloHit *
const pCaloHit : caloHitList2)
490 const CartesianVector &positionVector2(pCaloHit->GetPositionVector());
492 const float distanceSquared((positionVector1 - positionVector2).GetMagnitudeSquared());
494 if (distanceSquared < minDistanceSquared)
496 minDistanceSquared = distanceSquared;
497 closestPosition1 = positionVector1;
498 closestPosition2 = positionVector2;
499 distanceFound =
true;
505 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
507 outputPosition1 = closestPosition1;
508 outputPosition2 = closestPosition2;
static MCProcess GetLeadingProcess(const pandora::MCParticle *const pMCParticle)
Return the MCProcess of the leading particle (tier 1) in the delta ray/michel hierarchy.
static float GetClosestDistance(const pandora::Cluster *const pCluster, const pandora::CartesianPointVector &cartesianPointVector)
Get closest distance between a specified cluster and list of positions.
static void AddInPostBremsstrahlungHits(const LeadingMCParticleToPostBremsstrahlungHitList &leadingMCParticleToPostBremsstrahlungHitList, const float maxBremsstrahlungSeparation, LArMCParticleHelper::MCContributionMap &leadingMCToTrueHitListMap)
Identify the reconstructable post-bremsstrahlung radiation hits.
static void GetMuonPfoContaminationContribution(const pandora::CaloHitList &cosmicRayPfoHitList, const pandora::CaloHitList &leadingMCHitList, pandora::CaloHitList &leadingHitsInParentCosmicRay)
Determine the leading MCParticle hits within a cosmic ray pfo hit list.
static void SelectReconstructableLeadingParticles(const pandora::MCParticleList *pMCParticleList, const pandora::CaloHitList *pCaloHitList, const ValidationParameters ¶meters, const pandora::CaloHitList &recoMuonHitList, LArMCParticleHelper::MCContributionMap &selectedMCParticlesToHitsMap)
Select target, reconstructable mc particles in the cosmic ray hierarchy.
std::unordered_map< const pandora::MCParticle *, pandora::CaloHitList > MCContributionMap
bool m_selectInputHits
whether to select input hits
static bool RejectBremsstrahlungHits(const pandora::CaloHit *const pCaloHit, LeadingMCParticleToPostBremsstrahlungHitList &leadingMCParticleToPostBremsstrahlungHitList)
Identify and record the hits that are post-bremsstralung radiation in a cosmic ray hierarchy...
static const pandora::MCParticle * GetLeadingParticle(const pandora::MCParticle *const pMCParticle)
Return leading particle in the delta ray/michel hierarchy i.e. tier below cosmic ray.
MCProcess GetProcess() const
Get the process.
std::unordered_map< const pandora::CaloHit *, const pandora::MCParticle * > CaloHitToMCMap
static int GetHierarchyTier(const pandora::MCParticle *const pMCParticle)
Determine the position in the hierarchy for the MCParticle.
static void SelectLeadingMCParticles(const pandora::MCParticleList *pMCParticleList, pandora::MCParticleVector &selectedParticles)
Select all tier 0 and tier 1 MCParticles in cosmic ray hierarchies from an input list.
static void SelectParticlesByHitCount(const pandora::MCParticleVector &candidateTargets, const MCContributionMap &mcToTrueHitListMap, const MCRelationMap &mcToTargetMCMap, const PrimaryParameters ¶meters, MCContributionMap &selectedMCParticlesToHitsMap)
Filter an input vector of MCParticles to ensure they have sufficient good hits to be reconstructable...
static bool SortByMomentum(const pandora::MCParticle *const pLhs, const pandora::MCParticle *const pRhs)
Sort mc particles by their momentum.
static bool IsCosmicRay(const pandora::MCParticle *const pMCParticle)
Return true if passed a primary cosmic ray MCParticle.
float m_maxBremsstrahlungSeparation
The maximum separation of a reconstructable post-bremsstrahlung hit from the pre-radiation hits...
std::map< const pandora::MCParticle *, pandora::CaloHitList > LeadingMCParticleToPostBremsstrahlungHitList
static bool IsMichel(const pandora::MCParticle *const pMCParticle)
Return true if input MCParticle is a child of a cosmic ray and has 'decay' process tag...
Header file for the lar monte carlo particle helper helper class.
Header file for the cluster helper class.
Header file for the muon leading helper class.
static void GetMCParticleToCaloHitMatches(const pandora::CaloHitList *const pCaloHitList, const MCRelationMap &mcToTargetMCMap, CaloHitToMCMap &hitToMCMap, MCContributionMap &mcToTrueHitListMap)
Match calo hits to their parent particles.
static void GetClosestPositions(const pandora::CartesianPointVector &cartesianPointVector1, const pandora::Cluster *const pCluster2, pandora::CartesianVector &outputPosition1, pandora::CartesianVector &outputPosition2)
Get the closest positions between a list of positions and a cluster.
static void GetPfoMatchContamination(const pandora::MCParticle *const pLeadingParticle, const pandora::CaloHitList &matchedPfoHitList, pandora::CaloHitList &parentTrackHits, pandora::CaloHitList &otherTrackHits, pandora::CaloHitList &otherShowerHits)
Separate a leading pfo hit list according to the true owner of the hit e.g. other shower...
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
static const pandora::MCParticle * GetParentMCParticle(const pandora::MCParticle *const pMCParticle)
Get the parent mc particle.
static void SelectCaloHits(const pandora::CaloHitList *const pCaloHitList, const LArMCParticleHelper::MCRelationMap &mcToTargetMCMap, pandora::CaloHitList &selectedCaloHitList, const bool selectInputHits, const float minHitSharingFraction, const pandora::CaloHitList &recoMuonHitList, LeadingMCParticleToPostBremsstrahlungHitList &leadingMCParticleToPostBremsstrahlungHitList)
Select a subset of calo hits representing those that represent "reconstructable" regions of the event...
float m_minHitSharingFraction
the minimum Hit sharing fraction
static void GetMCToLeadingMap(const pandora::MCParticleList *const pMCParticleList, LArMCParticleHelper::MCRelationMap &mcToLeadingMap)
Construct the hierarchy folding map (cosmic rays folded to themselves, delta ray/michel hierarchy fol...
static pandora::StatusCode GetClosestPosition(const pandora::CartesianVector &referencePoint, const pandora::CartesianPointVector &cartesianPointVector, const pandora::Cluster *const pCluster, const float maxDistanceToCluster, const float maxDistanceToReferencePoint, pandora::CartesianVector &closestPosition)
Get the closest position from an input list of projected positions that lies close to both a referenc...
static void GetAllAncestorMCParticles(const pandora::MCParticle *const pMCParticle, pandora::MCParticleList &ancestorMCParticleList)
Get all ancestor mc particles.
LArMCParticleHelper class.
ValidationParameters()
Constructor.
static bool IsMuonLeading(const pandora::MCParticle *const pMCParticle)
Return true if input MCParticle is in tier 1 of the cosmic ray hierarchy.
std::unordered_map< const pandora::MCParticle *, const pandora::MCParticle * > MCRelationMap
static bool IsDeltaRay(const pandora::MCParticle *const pMCParticle)
Return true if input MCParticle is a child of a cosmic ray and has 'delta ray' process tag...
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.
ValidationParameters class.