9 #include "Pandora/AlgorithmHeaders.h" 21 TwoDSlidingFitSplittingAndSplicingAlgorithm::TwoDSlidingFitSplittingAndSplicingAlgorithm() :
22 m_shortHalfWindowLayers(10),
23 m_longHalfWindowLayers(20),
24 m_minClusterLength(7.5
f),
25 m_vetoDisplacement(1.5
f),
26 m_runCosmicMode(false)
34 const ClusterList *pClusterList = NULL;
35 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*
this, pClusterList));
39 unsigned int nIterations(0);
41 while (++nIterations < 100)
63 this->
BuildClusterExtensionList(clusterVector, branchSlidingFitResultMap, replacementSlidingFitResultMap, intermediateList);
68 if (STATUS_CODE_SUCCESS != this->
RunSplitAndExtension(splitList, branchSlidingFitResultMap, replacementSlidingFitResultMap))
72 return STATUS_CODE_SUCCESS;
81 const Cluster *
const pCluster = *iter;
86 clusterVector.push_back(pCluster);
99 if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
106 if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
107 throw StatusCodeException(STATUS_CODE_FAILURE);
109 catch (StatusCodeException &statusCodeException)
111 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
112 throw statusCodeException;
127 const Cluster *
const pClusterI = *iterI;
131 const Cluster *
const pClusterJ = *iterJ;
133 if (pClusterI == pClusterJ)
143 if (branchSlidingFitResultMap.end() == iterBranchI || branchSlidingFitResultMap.end() == iterBranchJ ||
144 replacementSlidingFitResultMap.end() == iterReplacementI || replacementSlidingFitResultMap.end() == iterReplacementJ)
157 float branchChisqI(0.
f);
158 CartesianVector branchSplitPositionI(0.
f, 0.
f, 0.
f);
159 CartesianVector branchSplitDirectionI(0.
f, 0.
f, 0.
f);
160 CartesianVector replacementStartPositionJ(0.
f, 0.
f, 0.
f);
164 this->
FindBestSplitPosition(branchSlidingFitI, replacementSlidingFitJ, replacementStartPositionJ, branchSplitPositionI, branchSplitDirectionI);
165 branchChisqI = this->
CalculateBranchChi2(pClusterI, branchSplitPositionI, branchSplitDirectionI);
167 catch (StatusCodeException &)
172 float branchChisqJ(0.
f);
173 CartesianVector branchSplitPositionJ(0.
f, 0.
f, 0.
f);
174 CartesianVector branchSplitDirectionJ(0.
f, 0.
f, 0.
f);
175 CartesianVector replacementStartPositionI(0.
f, 0.
f, 0.
f);
179 this->
FindBestSplitPosition(branchSlidingFitJ, replacementSlidingFitI, replacementStartPositionI, branchSplitPositionJ, branchSplitDirectionJ);
180 branchChisqJ = this->
CalculateBranchChi2(pClusterJ, branchSplitPositionJ, branchSplitDirectionJ);
182 catch (StatusCodeException &)
187 if (branchChisqI > 0.
f && branchChisqJ > 0.
f)
189 const CartesianVector relativeDirection((branchSplitPositionJ - branchSplitPositionI).GetUnitVector());
191 if (branchSplitDirectionI.GetDotProduct(relativeDirection) > 0.f && branchSplitDirectionJ.GetDotProduct(relativeDirection) < 0.f)
195 const float newBranchChisqI(this->
CalculateBranchChi2(pClusterI, branchSplitPositionI, relativeDirection));
196 const float newBranchChisqJ(this->
CalculateBranchChi2(pClusterJ, branchSplitPositionJ, relativeDirection * -1.
f));
197 branchChisqI = newBranchChisqI;
198 branchChisqJ = newBranchChisqJ;
200 catch (StatusCodeException &)
207 if (branchChisqI > branchChisqJ)
209 clusterExtensionList.push_back(
210 ClusterExtension(pClusterI, pClusterJ, replacementStartPositionJ, branchSplitPositionI, branchSplitDirectionI));
213 else if (branchChisqJ > branchChisqI)
215 clusterExtensionList.push_back(
216 ClusterExtension(pClusterJ, pClusterI, replacementStartPositionI, branchSplitPositionJ, branchSplitDirectionJ));
227 ClusterList branchList;
228 for (
const auto &mapEntry : branchMap)
229 branchList.push_back(mapEntry.first);
232 ClusterList replacementList;
233 for (
const auto &mapEntry : replacementMap)
234 replacementList.push_back(mapEntry.first);
239 const CartesianVector &branchVertex = thisSplit.GetBranchVertex();
240 const CartesianVector &replacementVertex = thisSplit.GetReplacementVertex();
242 const float distanceSquared((branchVertex - replacementVertex).GetMagnitudeSquared());
245 bool branchVeto(
false), replacementVeto(
false);
248 for (
const Cluster *
const pBranchCluster : branchList)
252 if (slidingFit.GetCluster() == thisSplit.GetReplacementCluster() || slidingFit.GetCluster() == thisSplit.GetBranchCluster())
255 const float minDistanceSquared((replacementVertex - slidingFit.GetGlobalMinLayerPosition()).GetMagnitudeSquared());
256 const float maxDistanceSquared((replacementVertex - slidingFit.GetGlobalMaxLayerPosition()).GetMagnitudeSquared());
258 if (std::min(minDistanceSquared, maxDistanceSquared) < std::max(distanceSquared, vetoDistanceSquared))
266 for (
const Cluster *
const pReplacementCluster : replacementList)
270 if (slidingFit.GetCluster() == thisSplit.GetReplacementCluster() || slidingFit.GetCluster() == thisSplit.GetBranchCluster())
273 const float minDistanceSquared((branchVertex - slidingFit.GetGlobalMinLayerPosition()).GetMagnitudeSquared());
274 const float maxDistanceSquared((branchVertex - slidingFit.GetGlobalMaxLayerPosition()).GetMagnitudeSquared());
276 if (std::min(minDistanceSquared, maxDistanceSquared) < std::max(distanceSquared, vetoDistanceSquared))
278 replacementVeto =
true;
283 if (branchVeto || replacementVeto)
286 outputList.push_back(thisSplit);
293 const Cluster *
const pCluster,
const CartesianVector &splitPosition,
const CartesianVector &splitDirection)
const 295 CaloHitList principalCaloHitList, branchCaloHitList;
297 this->
SplitBranchCluster(pCluster, splitPosition, splitDirection, principalCaloHitList, branchCaloHitList);
299 float totalChi2(0.
f);
300 float totalHits(0.
f);
304 const CaloHit *
const pCaloHit = *iter;
306 const CartesianVector hitPosition(pCaloHit->GetPositionVector());
307 const CartesianVector projectedPosition(splitPosition + splitDirection * splitDirection.GetDotProduct(hitPosition - splitPosition));
309 totalChi2 += (hitPosition - projectedPosition).GetMagnitudeSquared();
314 return std::sqrt(totalChi2 / totalHits);
316 throw StatusCodeException(STATUS_CODE_NOT_ALLOWED);
322 const CartesianVector &splitDirection, CaloHitList &principalCaloHitList, CaloHitList &branchCaloHitList)
const 325 CaloHitList caloHitsToDistribute;
326 pCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitsToDistribute);
328 for (
CaloHitList::const_iterator iter = caloHitsToDistribute.begin(), iterEnd = caloHitsToDistribute.end(); iter != iterEnd; ++iter)
330 const CaloHit *
const pCaloHit = *iter;
332 if (splitDirection.GetDotProduct((pCaloHit->GetPositionVector() - splitPosition)) > 0.f)
334 branchCaloHitList.push_back(pCaloHit);
338 principalCaloHitList.push_back(pCaloHit);
342 if (branchCaloHitList.empty())
343 throw StatusCodeException(STATUS_CODE_NOT_ALLOWED);
351 bool foundSplit(
false);
359 const CartesianVector &branchSplitPosition = thisSplit.
GetBranchVertex();
368 if (branchResultMap.end() == iterBranch1 || branchResultMap.end() == iterBranch2 ||
369 replacementResultMap.end() == iterReplacement1 || replacementResultMap.end() == iterReplacement2)
372 PANDORA_RETURN_RESULT_IF(
373 STATUS_CODE_SUCCESS, !=, this->
ReplaceBranch(pBranchCluster, pReplacementCluster, branchSplitPosition, branchSplitDirection));
374 branchResultMap.erase(iterBranch1);
375 branchResultMap.erase(iterBranch2);
377 replacementResultMap.erase(iterReplacement1);
378 replacementResultMap.erase(iterReplacement2);
384 return STATUS_CODE_SUCCESS;
386 return STATUS_CODE_NOT_FOUND;
392 const Cluster *
const pReplacementCluster,
const CartesianVector &branchSplitPosition,
const CartesianVector &branchSplitDirection)
const 394 ClusterList clusterList;
395 clusterList.push_back(pBranchCluster);
396 clusterList.push_back(pReplacementCluster);
398 std::string clusterListToSaveName, clusterListToDeleteName;
399 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
400 PandoraContentApi::InitializeFragmentation(*
this, clusterList, clusterListToDeleteName, clusterListToSaveName));
403 PandoraContentApi::Cluster::Parameters principalParameters;
404 pReplacementCluster->GetOrderedCaloHitList().FillCaloHitList(principalParameters.m_caloHitList);
407 PandoraContentApi::Cluster::Parameters residualParameters;
409 pBranchCluster, branchSplitPosition, branchSplitDirection, principalParameters.m_caloHitList, residualParameters.m_caloHitList);
411 const Cluster *pPrincipalCluster(NULL), *pResidualCluster(NULL);
412 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
this, principalParameters, pPrincipalCluster));
413 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
this, residualParameters, pResidualCluster));
414 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::EndFragmentation(*
this, clusterListToSaveName, clusterListToDeleteName));
416 return STATUS_CODE_SUCCESS;
423 PANDORA_RETURN_RESULT_IF_AND_IF(
424 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ShortHalfWindow",
m_shortHalfWindowLayers));
426 PANDORA_RETURN_RESULT_IF_AND_IF(
427 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"LongHalfWindow",
m_longHalfWindowLayers));
429 PANDORA_RETURN_RESULT_IF_AND_IF(
430 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterLength",
m_minClusterLength));
432 PANDORA_RETURN_RESULT_IF_AND_IF(
433 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"VetoDisplacement",
m_vetoDisplacement));
435 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"CosmicMode",
m_runCosmicMode));
437 return STATUS_CODE_SUCCESS;
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position, then pulse-height.
const pandora::Cluster * GetReplacementCluster() const
return the address of the replacement Cluster
virtual void FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFit, const TwoDSlidingFitResult &replacementSlidingFit, pandora::CartesianVector &replacementStartPosition, pandora::CartesianVector &branchSplitPosition, pandora::CartesianVector &branchSplitDirection) const =0
Output the best split positions in branch and replacement clusters.
virtual pandora::StatusCode Run()
std::vector< ClusterExtension > ClusterExtensionList
unsigned int m_longHalfWindowLayers
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, const unsigned int halfWindowLayers, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
void BuildClusterExtensionList(const pandora::ClusterVector &clusterVector, const TwoDSlidingFitResultMap &branchResultMap, const TwoDSlidingFitResultMap &replacementResultMap, ClusterExtensionList &clusterExtensionList) const
Build a list of candidate splits.
const pandora::Cluster * GetBranchCluster() const
return the address of the branch Cluster
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
const pandora::CartesianVector & GetBranchVertex() const
return the split position of the branch cluster
const pandora::CartesianVector & GetBranchDirection() const
return the split direction of the branch cluster
Header file for the geometry helper class.
void GetListOfCleanClusters(const pandora::ClusterList *const pClusterList, pandora::ClusterVector &clusterVector) const
Populate cluster vector with subset of cluster list, containing clusters judged to be clean...
pandora::StatusCode RunSplitAndExtension(const ClusterExtensionList &splitList, TwoDSlidingFitResultMap &branchResultMap, TwoDSlidingFitResultMap &replacementResultMap) const
Run the machinary that performs the cluster splitting and extending.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::StatusCode ReplaceBranch(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &branchSplitPosition, const pandora::CartesianVector &branchSplitDirection) const
Remove a branch from a cluster and replace it with a second cluster.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
Header file for the two dimensional sliding fit splitting and splicing algorithm class.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
void SplitBranchCluster(const pandora::Cluster *const pCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection, pandora::CaloHitList &principalCaloHitList, pandora::CaloHitList &branchCaloHitList) const
Separate cluster into the branch hits to be split from the primary cluster.
float CalculateBranchChi2(const pandora::Cluster *const pCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection) const
Calculate RMS deviation of branch hits relative to the split direction.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
unsigned int m_shortHalfWindowLayers
void PruneClusterExtensionList(const ClusterExtensionList &inputList, const TwoDSlidingFitResultMap &branchResultMap, const TwoDSlidingFitResultMap &replacementResultMap, ClusterExtensionList &outputList) const
Finalize the list of candidate splits.
TwoDSlidingFitResult class.