9 #include "Pandora/AlgorithmHeaders.h" 21 TrackClusterCreationAlgorithm::TrackClusterCreationAlgorithm() :
22 m_mergeBackFilteredHits(true),
24 m_maxCaloHitSeparationSquared(1.3
f * 1.3
f),
25 m_minCaloHitSeparationSquared(0.4
f * 0.4
f),
26 m_closeSeparationSquared(0.9
f * 0.9
f),
35 const CaloHitList *pCaloHitList = NULL;
36 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*
this, pCaloHitList));
38 OrderedCaloHitList selectedCaloHitList, rejectedCaloHitList;
39 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->
FilterCaloHits(pCaloHitList, selectedCaloHitList, rejectedCaloHitList));
47 this->
IdentifyJoins(selectedCaloHitList, forwardHitAssociationMap, backwardHitAssociationMap, hitJoinMap);
48 this->
CreateClusters(selectedCaloHitList, hitJoinMap, hitToClusterMap);
51 this->
CreateClusters(rejectedCaloHitList, hitJoinMap, hitToClusterMap);
53 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->
AddFilteredCaloHits(selectedCaloHitList, rejectedCaloHitList, hitToClusterMap));
55 return STATUS_CODE_SUCCESS;
61 const CaloHitList *
const pCaloHitList, OrderedCaloHitList &selectedCaloHitList, OrderedCaloHitList &rejectedCaloHitList)
const 63 CaloHitList availableHitList;
65 for (
const CaloHit *
const pCaloHit : *pCaloHitList)
67 if (PandoraContentApi::IsAvailable(*
this, pCaloHit) && pCaloHit->GetMipEquivalentEnergy() >=
m_minMipFraction)
68 availableHitList.push_back(pCaloHit);
71 if (availableHitList.empty())
72 return STATUS_CODE_SUCCESS;
74 HitType view{availableHitList.front()->GetHitType()};
78 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, selectedCaloHitList.Add(availableHitList));
82 CaloHitVector caloHits(iter->second->begin(), iter->second->end());
85 for (
const CaloHit *
const pCaloHitI : caloHits)
87 bool useCaloHit(
true);
89 for (
const CaloHit *
const pCaloHitJ : caloHits)
91 if (pCaloHitI == pCaloHitJ)
94 if ((pCaloHitI->GetMipEquivalentEnergy() < pCaloHitJ->GetMipEquivalentEnergy()) &&
95 ((pCaloHitI->GetPositionVector() - pCaloHitJ->GetPositionVector()).GetMagnitudeSquared() < minSeparationSquaredAdjusted))
103 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, rejectedCaloHitList.Add(pCaloHitI));
107 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, selectedCaloHitList.Remove(rejectedCaloHitList));
108 return STATUS_CODE_SUCCESS;
114 const OrderedCaloHitList &selectedCaloHitList,
const OrderedCaloHitList &rejectedCaloHitList,
HitToClusterMap &hitToClusterMap)
const 118 CaloHitList *pCaloHitList = NULL;
119 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, selectedCaloHitList.GetCaloHitsInPseudoLayer(iter->first, pCaloHitList));
121 if (!pCaloHitList || pCaloHitList->empty())
123 HitType view{pCaloHitList->front()->GetHitType()};
127 CaloHitSet unavailableHits;
129 CaloHitVector inputAvailableHits(iter->second->begin(), iter->second->end());
132 CaloHitVector clusteredHits(pCaloHitList->begin(), pCaloHitList->end());
140 CaloHitVector newClusteredHits;
142 for (
const CaloHit *
const pCaloHitI : inputAvailableHits)
144 if (unavailableHits.count(pCaloHitI))
147 if (hitToClusterMap.end() != hitToClusterMap.find(pCaloHitI))
150 const CaloHit *pClosestHit = NULL;
151 float closestSeparationSquared(minSeparationSquaredAdjusted);
153 for (
const CaloHit *
const pCaloHitJ : clusteredHits)
155 if (pCaloHitI->GetMipEquivalentEnergy() > pCaloHitJ->GetMipEquivalentEnergy())
158 const float separationSquared((pCaloHitI->GetPositionVector() - pCaloHitJ->GetPositionVector()).GetMagnitudeSquared());
160 if (separationSquared < closestSeparationSquared)
162 closestSeparationSquared = separationSquared;
163 pClosestHit = pCaloHitJ;
172 if (hitToClusterMap.end() == mapIter)
173 throw StatusCodeException(STATUS_CODE_FAILURE);
175 const Cluster *
const pCluster = mapIter->second;
176 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*
this, pCluster, pCaloHitI));
177 (void)hitToClusterMap.insert(HitToClusterMap::value_type(pCaloHitI, pCluster));
179 newClusteredHits.push_back(pCaloHitI);
183 for (
const CaloHit *
const pCaloHit : newClusteredHits)
185 clusteredHits.push_back(pCaloHit);
186 unavailableHits.insert(pCaloHit);
191 return STATUS_CODE_SUCCESS;
201 unsigned int nLayersConsidered(0);
203 CaloHitVector caloHitsI(iterI->second->begin(), iterI->second->end());
207 (nLayersConsidered++ <=
m_maxGapLayers + 1) && (iterJ != iterJEnd); ++iterJ)
209 if (iterJ->first == iterI->first || iterJ->first > iterI->first +
m_maxGapLayers + 1)
212 CaloHitVector caloHitsJ(iterJ->second->begin(), iterJ->second->end());
215 for (
const CaloHit *
const pCaloHitI : caloHitsI)
217 for (
const CaloHit *
const pCaloHitJ : caloHitsJ)
231 CaloHitVector caloHits(iter->second->begin(), iter->second->end());
234 for (
const CaloHit *
const pCaloHit : caloHits)
237 const CaloHit *
const pForwardHit((forwardHitAssociationMap.end() == fwdIter) ? NULL : fwdIter->second.GetPrimaryTarget());
240 const CaloHit *
const pForwardHitCheck((backwardHitAssociationMap.end() == fwdCheckIter) ? NULL : fwdCheckIter->second.GetPrimaryTarget());
242 if ((NULL != pForwardHit) && (pForwardHitCheck != pCaloHit))
246 const CaloHit *
const pBackwardHit((backwardHitAssociationMap.end() == bwdIter) ? NULL : bwdIter->second.GetPrimaryTarget());
249 const CaloHit *
const pBackwardHitCheck((forwardHitAssociationMap.end() == bwdCheckIter) ? NULL : bwdCheckIter->second.GetPrimaryTarget());
251 if ((NULL != pBackwardHit) && (pBackwardHitCheck != pCaloHit))
264 CaloHitVector caloHits(iter->second->begin(), iter->second->end());
267 for (
const CaloHit *
const pCaloHit : caloHits)
269 const CaloHit *
const pForwardJoinHit = this->
GetJoinHit(pCaloHit, forwardHitAssociationMap, backwardHitAssociationMap);
270 const CaloHit *
const pBackwardJoinHit = this->
GetJoinHit(pForwardJoinHit, backwardHitAssociationMap, forwardHitAssociationMap);
272 if ((NULL == pForwardJoinHit) || (NULL == pBackwardJoinHit) || (pBackwardJoinHit != pCaloHit))
277 if (hitJoinMap.end() == joinIter)
278 hitJoinMap.insert(HitJoinMap::value_type(pCaloHit, pForwardJoinHit));
280 if ((hitJoinMap.end() != joinIter) && (joinIter->second != pForwardJoinHit))
281 throw StatusCodeException(STATUS_CODE_FAILURE);
293 CaloHitVector caloHits(iter->second->begin(), iter->second->end());
296 for (
const CaloHit *
const pCaloHit : caloHits)
298 const Cluster *pCluster = NULL;
302 if (hitToClusterMap.end() == mapIter)
304 PandoraContentApi::Cluster::Parameters parameters;
305 parameters.m_caloHitList.push_back(pCaloHit);
306 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
this, parameters, pCluster));
307 hitToClusterMap.insert(HitToClusterMap::value_type(pCaloHit, pCluster));
311 pCluster = mapIter->second;
316 if (hitJoinMap.end() == joinIter)
319 if (hitToClusterMap.end() != hitToClusterMap.find(joinIter->second))
320 throw StatusCodeException(STATUS_CODE_FAILURE);
322 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*
this, pCluster, joinIter->second));
323 hitToClusterMap.insert(HitToClusterMap::value_type(joinIter->second, pCluster));
333 HitType view{pCaloHitI->GetHitType()};
337 const float distanceSquared((pCaloHitJ->GetPositionVector() - pCaloHitI->GetPositionVector()).GetMagnitudeSquared());
339 if (distanceSquared > maxSeparationSquaredAdjusted)
344 if (forwardHitAssociationMap.end() == forwardIter)
346 forwardHitAssociationMap.insert(HitAssociationMap::value_type(pCaloHitI,
HitAssociation(pCaloHitJ, distanceSquared)));
348 else if (distanceSquared < forwardIter->
second.GetPrimaryDistanceSquared())
355 if (backwardHitAssociationMap.end() == backwardIter)
357 backwardHitAssociationMap.insert(HitAssociationMap::value_type(pCaloHitJ,
HitAssociation(pCaloHitI, distanceSquared)));
359 else if (distanceSquared < backwardIter->
second.GetPrimaryDistanceSquared())
361 backwardIter->second =
HitAssociation(pCaloHitI, distanceSquared);
370 HitType view{pCaloHitI->GetHitType()};
377 if ((forwardHitAssociationMap.end() == forwardIter) || (backwardHitAssociationMap.end() == backwardIter))
409 if (hitAssociationMapI.end() == iterI)
412 const CaloHit *
const pPrimaryTarget = iterI->second.GetPrimaryTarget();
413 const CaloHit *
const pSecondaryTarget = iterI->second.GetSecondaryTarget();
415 if (NULL == pSecondaryTarget)
416 return pPrimaryTarget;
418 unsigned int primaryNSteps(0), secondaryNSteps(0);
419 const CaloHit *
const pPrimaryTrace = this->
TraceHitAssociation(pPrimaryTarget, hitAssociationMapI, hitAssociationMapJ, primaryNSteps);
420 const CaloHit *
const pSecondaryTrace = this->
TraceHitAssociation(pSecondaryTarget, hitAssociationMapI, hitAssociationMapJ, secondaryNSteps);
422 if ((pPrimaryTrace == pSecondaryTrace) || (secondaryNSteps < 5))
423 return pPrimaryTarget;
434 const CaloHit *pThisHit = pCaloHit;
435 const CaloHit *pLastHit = pCaloHit;
443 if (hitAssociationMapI.end() == iterI)
446 pLastHit = iterI->second.GetPrimaryTarget();
449 if (hitAssociationMapJ.end() == iterJ)
452 if (iterJ->second.GetPrimaryTarget() != pThisHit)
463 PANDORA_RETURN_RESULT_IF_AND_IF(
464 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MergeBackFilteredHits",
m_mergeBackFilteredHits));
466 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxGapLayers",
m_maxGapLayers));
469 PANDORA_RETURN_RESULT_IF_AND_IF(
470 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxCaloHitSeparation", maxCaloHitSeparation));
474 PANDORA_RETURN_RESULT_IF_AND_IF(
475 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinCaloHitSeparation", minCaloHitSeparation));
479 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"CloseSeparation", closeSeparation));
482 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinMipFraction",
m_minMipFraction));
484 return STATUS_CODE_SUCCESS;
std::unordered_map< const pandora::CaloHit *, const pandora::Cluster * > HitToClusterMap
pandora::StatusCode FilterCaloHits(const pandora::CaloHitList *const pCaloHitList, pandora::OrderedCaloHitList &selectedCaloHitList, pandora::OrderedCaloHitList &rejectedCaloHitList) const
Filter out low pulse height hits in close proximity to high pulse height hits.
static float GetWirePitchRatio(const pandora::Pandora &pandora, const pandora::HitType view)
Return the ratio of the wire pitch of the specified view to the minimum wire pitch for the detector...
float m_maxCaloHitSeparationSquared
Square of maximum calo hit separation.
pandora::StatusCode AddFilteredCaloHits(const pandora::OrderedCaloHitList &selectedCaloHitList, const pandora::OrderedCaloHitList &rejectedCaloHitList, HitToClusterMap &hitToClusterMap) const
Merge previously filtered hits back into their associated clusters.
void IdentifyJoins(const pandora::OrderedCaloHitList &orderedCaloHitList, const HitAssociationMap &forwardHitAssociationMap, const HitAssociationMap &backwardHitAssociationMap, HitJoinMap &hitJoinMap) const
Identify final hit joins for use in cluster formation.
pandora::StatusCode Run()
Header file for the geometry helper class.
float m_minMipFraction
Minimum fraction of a MIP to consider a hit.
float GetSecondaryDistanceSquared() const
Get the secondary distance squared.
Header file for the cluster creation algorithm class.
std::unordered_map< const pandora::CaloHit *, HitAssociation > HitAssociationMap
void MakePrimaryAssociations(const pandora::OrderedCaloHitList &orderedCaloHitList, HitAssociationMap &forwardHitAssociationMap, HitAssociationMap &backwardHitAssociationMap) const
Control primary association formation.
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.
const pandora::CaloHit * GetPrimaryTarget() const
Get the primary target.
std::unordered_map< const pandora::CaloHit *, const pandora::CaloHit * > HitJoinMap
float GetPrimaryDistanceSquared() const
Get the primary distance squared.
void SetSecondaryTarget(const pandora::CaloHit *const pSecondaryTarget, const float secondaryDistanceSquared)
Set secondary target.
void CreateClusters(const pandora::OrderedCaloHitList &orderedCaloHitList, const HitJoinMap &hitJoinMap, HitToClusterMap &hitToClusterMap) const
Final cluster formation.
void CreatePrimaryAssociation(const pandora::CaloHit *const pCaloHitI, const pandora::CaloHit *const pCaloHitJ, HitAssociationMap &forwardHitAssociationMap, HitAssociationMap &backwardHitAssociationMap) const
Create primary association if appropriate, hitI<->hitJ.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
unsigned int m_maxGapLayers
Maximum number of layers for a gap.
const pandora::CaloHit * GetJoinHit(const pandora::CaloHit *const pCaloHit, const HitAssociationMap &hitAssociationMapI, const HitAssociationMap &hitAssociationMapJ) const
Get hit to join by tracing associations via map I, checking via map J.
bool m_mergeBackFilteredHits
Merge rejected hits into their associated clusters.
void MakeSecondaryAssociations(const pandora::OrderedCaloHitList &orderedCaloHitList, HitAssociationMap &forwardHitAssociationMap, HitAssociationMap &backwardHitAssociationMap) const
Control secondary association formation.
second_as<> second
Type of time stored in seconds, in double precision.
void CreateSecondaryAssociation(const pandora::CaloHit *const pCaloHitI, const pandora::CaloHit *const pCaloHitJ, HitAssociationMap &forwardHitAssociationMap, HitAssociationMap &backwardHitAssociationMap) const
Create secondary association if appropriate, hitI<->hitJ.
float m_minCaloHitSeparationSquared
Square of minimum calo hit separation.
float m_closeSeparationSquared
Length scale (squared) for close hit separation.
const pandora::CaloHit * TraceHitAssociation(const pandora::CaloHit *const pCaloHit, const HitAssociationMap &hitAssociationMapI, const HitAssociationMap &hitAssociationMapJ, unsigned int &nSteps) const
Get last hit obtained by tracing associations via map I, checking via map J.