9 #include "Pandora/AlgorithmHeaders.h" 26 ParticleRecoveryAlgorithm::ParticleRecoveryAlgorithm() :
28 m_minClusterCaloHits(20),
29 m_minClusterLengthSquared(5.
f * 5.
f),
30 m_minClusterXSpan(0.25
f),
31 m_vertexClusterMode(false),
32 m_minVertexLongitudinalDistance(-2.5
f),
33 m_maxVertexTransverseDistance(1.5
f),
34 m_minXOverlapFraction(0.75
f),
35 m_minXOverlapFractionGaps(0.75
f),
36 m_sampleStepSize(0.5
f),
37 m_slidingFitHalfWindow(10),
46 ClusterList inputClusterListU, inputClusterListV, inputClusterListW;
47 this->
GetInputClusters(inputClusterListU, inputClusterListV, inputClusterListW);
49 ClusterList selectedClusterListU, selectedClusterListV, selectedClusterListW;
55 this->
FindOverlaps(selectedClusterListU, selectedClusterListV, overlapTensor);
56 this->
FindOverlaps(selectedClusterListV, selectedClusterListW, overlapTensor);
57 this->
FindOverlaps(selectedClusterListW, selectedClusterListU, overlapTensor);
60 return STATUS_CODE_SUCCESS;
69 const ClusterList *pClusterList(NULL);
70 PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*
this, *iter, pClusterList));
72 if (!pClusterList || pClusterList->empty())
74 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
75 std::cout <<
"ParticleRecoveryAlgorithm: unable to find cluster list " << *iter << std::endl;
82 const Cluster *
const pCluster(*cIter);
85 if ((TPC_VIEW_U != hitType) && (TPC_VIEW_V != hitType) && (TPC_VIEW_W != hitType))
86 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
88 ClusterList &clusterList((TPC_VIEW_U == hitType) ? inputClusterListU
89 : (TPC_VIEW_V == hitType) ? inputClusterListV
91 clusterList.push_back(pCluster);
102 ClusterList vertexClusterList;
118 const Cluster *
const pCluster = *iter;
120 if (!pCluster->IsAvailable())
129 float xMin(0.
f), xMax(0.
f);
130 pCluster->GetClusterSpanX(xMin, xMax);
135 selectedClusterList.push_back(pCluster);
143 CartesianPointVector vertexList;
149 if (!(*iter)->IsAvailable())
156 catch (StatusCodeException &)
165 const Cluster *
const pCluster = *iter;
167 if (!pCluster->IsAvailable())
177 selectedClusterList.push_back(pCluster);
182 catch (StatusCodeException &)
207 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
209 if ((0 == pCluster1->GetNCaloHits()) || (0 == pCluster2->GetNCaloHits()))
210 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
212 float xMin1(0.
f), xMax1(0.
f), xMin2(0.
f), xMax2(0.
f);
213 pCluster1->GetClusterSpanX(xMin1, xMax1);
214 pCluster2->GetClusterSpanX(xMin2, xMax2);
216 const float xSpan1(xMax1 - xMin1), xSpan2(xMax2 - xMin2);
218 if ((xSpan1 < std::numeric_limits<float>::epsilon()) || (xSpan2 < std::numeric_limits<float>::epsilon()))
219 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
221 const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
223 float xOverlapFraction1(xOverlap / xSpan1), xOverlapFraction2(xOverlap / xSpan2);
237 const Cluster *
const pCluster2,
const float xMin2,
const float xMax2,
float &xOverlapFraction1,
float &xOverlapFraction2)
const 239 if (PandoraContentApi::GetGeometry(*this)->GetDetectorGapList().empty())
242 const float xMin(std::min(xMin1, xMin2));
243 const float xMax(std::max(xMax1, xMax2));
244 float xMinEff1(xMin1), xMaxEff1(xMax1), xMinEff2(xMin2), xMaxEff2(xMax2);
249 const float effectiveXSpan1(xMaxEff1 - xMinEff1), effectiveXSpan2(xMaxEff2 - xMinEff2);
250 const float effectiveXOverlapSpan(std::min(xMaxEff1, xMaxEff2) - std::max(xMinEff1, xMinEff2));
252 if ((effectiveXSpan1 > std::numeric_limits<float>::epsilon()) && (effectiveXSpan2 > std::numeric_limits<float>::epsilon()))
254 xOverlapFraction1 = std::min(1.
f, (effectiveXOverlapSpan / effectiveXSpan1));
255 xOverlapFraction2 = std::min(1.
f, (effectiveXOverlapSpan / effectiveXSpan2));
262 const pandora::Cluster *
const pCluster,
const float xMin,
const float xMax,
float &xMinEff,
float &xMaxEff)
const 271 const int nSamplingPointsLeft(1 + static_cast<int>((xMinEff - xMin) /
m_sampleStepSize));
272 const int nSamplingPointsRight(1 + static_cast<int>((xMax - xMaxEff) /
m_sampleStepSize));
273 float dxMin(0.
f), dxMax(0.
f);
275 for (
int iSample = 1; iSample <= nSamplingPointsLeft; ++iSample)
277 const float xSample(std::max(xMin, xMinEff - static_cast<float>(iSample) *
m_sampleStepSize));
282 dxMin = xMinEff - xSample;
285 for (
int iSample = 1; iSample <= nSamplingPointsRight; ++iSample)
287 const float xSample(std::min(xMax, xMaxEff + static_cast<float>(iSample) *
m_sampleStepSize));
292 dxMax = xSample - xMaxEff;
295 xMinEff = xMinEff - dxMin;
296 xMaxEff = xMaxEff + dxMax;
298 catch (StatusCodeException &)
310 for (
const Cluster *
const pKeyCluster : sortedKeyClusters)
312 ClusterList clusterListU, clusterListV, clusterListW;
314 overlapTensor.
GetConnectedElements(pKeyCluster,
true, clusterListU, clusterListV, clusterListW);
315 const unsigned int nU(clusterListU.size()), nV(clusterListV.size()), nW(clusterListW.size());
317 if ((0 == nU * nV) && (0 == nV * nW) && (0 == nW * nU))
320 ClusterList clusterListAll;
321 clusterListAll.insert(clusterListAll.end(), clusterListU.begin(), clusterListU.end());
322 clusterListAll.insert(clusterListAll.end(), clusterListV.begin(), clusterListV.end());
323 clusterListAll.insert(clusterListAll.end(), clusterListW.begin(), clusterListW.end());
325 if ((1 == nU * nV * nW) && this->
CheckConsistency(*(clusterListU.begin()), *(clusterListV.begin()), *(clusterListW.begin())))
329 else if ((0 == nU * nV * nW) && ((1 == nU && 1 == nV) || (1 == nV && 1 == nW) || (1 == nW && 1 == nU)))
346 float xMinU(0.
f), xMinV(0.
f), xMinW(0.
f), xMaxU(0.
f), xMaxV(0.
f), xMaxW(0.
f);
347 pClusterU->GetClusterSpanX(xMinU, xMaxU);
348 pClusterV->GetClusterSpanX(xMinV, xMaxV);
349 pClusterW->GetClusterSpanX(xMinW, xMaxW);
351 const float xMin(std::max(xMinU, std::max(xMinV, xMinW)));
352 const float xMax(std::min(xMaxU, std::min(xMaxV, xMaxW)));
353 const float xOverlap(xMax - xMin);
355 if (xOverlap < std::numeric_limits<float>::epsilon())
359 float u(std::numeric_limits<float>::max()), v(std::numeric_limits<float>::max()),
w(std::numeric_limits<float>::max());
372 const float pseudoChi2(((u - vw2u) * (u - vw2u) + (v - wu2v) * (v - wu2v) + (w - uv2w) * (w - uv2w)) / 3.
f);
384 if (clusterList.size() < 2)
385 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
387 const PfoList *pPfoList = NULL;
388 std::string pfoListName;
389 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pPfoList, pfoListName));
392 PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
393 pfoParameters.m_particleId = MU_MINUS;
394 pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
395 pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
396 pfoParameters.m_energy = 0.f;
397 pfoParameters.m_momentum = CartesianVector(0.
f, 0.
f, 0.
f);
398 pfoParameters.m_clusterList = clusterList;
400 const ParticleFlowObject *pPfo(NULL);
401 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*
this, pfoParameters, pPfo));
403 if (!pPfoList->empty())
405 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*
this,
m_outputPfoListName));
406 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*
this,
m_outputPfoListName));
418 const Cluster *
const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
419 const Cluster *
const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
420 const Cluster *
const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
422 if (pClusterU && pClusterV && !pClusterW)
424 m_clusterNavigationMapUV[pClusterU].push_back(pClusterV);
426 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterU))
427 m_keyClusters.push_back(pClusterU);
429 else if (!pClusterU && pClusterV && pClusterW)
431 m_clusterNavigationMapVW[pClusterV].push_back(pClusterW);
433 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterV))
434 m_keyClusters.push_back(pClusterV);
436 else if (pClusterU && !pClusterV && pClusterW)
438 m_clusterNavigationMapWU[pClusterW].push_back(pClusterU);
440 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterW))
441 m_keyClusters.push_back(pClusterW);
445 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
452 ClusterList &clusterListU, ClusterList &clusterListV, ClusterList &clusterListW)
const 454 if (ignoreUnavailable && !pCluster->IsAvailable())
459 if (!((TPC_VIEW_U == hitType) || (TPC_VIEW_V == hitType) || (TPC_VIEW_W == hitType)))
460 throw StatusCodeException(STATUS_CODE_FAILURE);
462 ClusterList &clusterList((TPC_VIEW_U == hitType) ? clusterListU : (TPC_VIEW_V == hitType) ? clusterListV : clusterListW);
464 : (TPC_VIEW_V == hitType) ? m_clusterNavigationMapVW
465 : m_clusterNavigationMapWU);
467 if (clusterList.end() != std::find(clusterList.begin(), clusterList.end(), pCluster))
470 clusterList.push_back(pCluster);
474 if (navigationMap.end() == iter)
478 this->GetConnectedElements(*cIter, ignoreUnavailable, clusterListU, clusterListV, clusterListW);
486 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
"InputClusterListNames",
m_inputClusterListNames));
487 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputPfoListName",
m_outputPfoListName));
489 PANDORA_RETURN_RESULT_IF_AND_IF(
490 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterCaloHits",
m_minClusterCaloHits));
492 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"CheckGaps",
m_checkGaps));
495 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterLength", minClusterLength));
498 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterXSpan",
m_minClusterXSpan));
500 PANDORA_RETURN_RESULT_IF_AND_IF(
501 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"VertexClusterMode",
m_vertexClusterMode));
503 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
506 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
509 PANDORA_RETURN_RESULT_IF_AND_IF(
510 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinXOverlapFraction",
m_minXOverlapFraction));
512 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
515 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SampleStepSize",
m_sampleStepSize));
519 std::cout <<
"ParticleRecoveryAlgorithm: Invalid value for SampleStepSize " <<
m_sampleStepSize << std::endl;
520 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
523 PANDORA_RETURN_RESULT_IF_AND_IF(
524 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitHalfWindow",
m_slidingFitHalfWindow));
526 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"PseudoChi2Cut",
m_pseudoChi2Cut));
528 return STATUS_CODE_SUCCESS;
pandora::StringVector m_inputClusterListNames
The list of cluster list names.
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.
unsigned int m_minClusterCaloHits
The min number of hits in base cluster selection method.
void GetInputClusters(pandora::ClusterList &inputClusterListU, pandora::ClusterList &inputClusterListV, pandora::ClusterList &inputClusterListW) const
Get the input cluster lists for processing in this algorithm.
const pandora::ClusterList & GetKeyClusters() const
Get the list of key clusters.
Header file for the lar pointing cluster class.
float m_minClusterLengthSquared
The min length (squared) in base cluster selection method.
void ExamineTensor(const SimpleOverlapTensor &overlapTensor) const
Identify unambiguous cluster overlaps and resolve ambiguous overlaps, creating new track particles...
float m_maxVertexTransverseDistance
Vertex association check: max transverse distance cut.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterNavigationMap
bool CheckConsistency(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW) const
Whether a trio of clusters are consistent with representing projections of the same 3d trajectory...
LArPointingCluster class.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
static pandora::StatusCode GetAverageZ(const pandora::Cluster *const pCluster, const float xmin, const float xmax, float &averageZ)
Get average Z positions of the calo hits in a cluster in range xmin to xmax.
float m_minXOverlapFractionGaps
The min x overlap fraction when there are gaps involved.
void SelectInputClusters(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters for processing in this algorithm.
Header file for the geometry helper class.
void CalculateEffectiveSpan(const pandora::Cluster *const pCluster, const float xMin, const float xMax, float &xMinEff, float &xMaxEff) const
Calculate effective span for a given clsuter taking gaps into account.
void CreateTrackParticle(const pandora::ClusterList &clusterList) const
Create and save a track particle containing the provided clusters.
bool IsOverlap(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2) const
Whether two clusters overlap convincingly in x.
void AddAssociation(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2)
Add an association between two clusters to the simple overlap tensor.
Header file for the cluster helper class.
void FindOverlaps(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2, SimpleOverlapTensor &overlapTensor) const
Find cluster overlaps and record these in the overlap tensor.
const Vertex & GetOuterVertex() const
Get the outer vertex.
pandora::StatusCode Run()
const Vertex & GetInnerVertex() const
Get the inner vertex.
float m_minVertexLongitudinalDistance
Vertex association check: min longitudinal distance cut.
static float MergeTwoPositions(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const float position1, const float position2)
Merge two views (U,V) to give a third view (Z).
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
float m_minClusterXSpan
The min x span required in order to consider a cluster.
float m_pseudoChi2Cut
The selection cut on the matched chi2.
void CalculateEffectiveOverlapFractions(const pandora::Cluster *const pCluster1, const float xMin1, const float xMax1, const pandora::Cluster *const pCluster2, const float xMin2, const float xMax2, float &xOverlapFraction1, float &xOverlapFraction2) const
Calculate effective overlap fractions taking into account gaps.
static bool IsXSamplingPointInGap(const pandora::Pandora &pandora, const float xSample, const TwoDSlidingFitResult &slidingFitResult, const float gapTolerance=0.f)
Whether there is a gap in a cluster (described via its sliding fit result) at a specified x sampling ...
unsigned int m_slidingFitHalfWindow
The half window for the fit sliding result constructor.
float m_minXOverlapFraction
The min x overlap fraction required in order to id overlapping clusters.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
SimpleOverlapTensor class.
bool m_checkGaps
Whether to check for gaps in the calculation of the overlap.
bool m_vertexClusterMode
Whether to demand clusters are associated with vertices of existing particles.
Header file for the track recovery algorithm class.
std::string m_outputPfoListName
The output pfo list name.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void StandardClusterSelection(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters for processing in this algorithm.
float m_sampleStepSize
The sampling step size used in association checks, units cm.
void VertexClusterSelection(const pandora::ClusterList &inputClusterList, pandora::ClusterList &selectedClusterList) const
Select a subset of input clusters nodally associated with the vertices of existing particles...
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
TwoDSlidingFitResult class.
static bool IsNode(const pandora::CartesianVector &parentVertex, const LArPointingCluster::Vertex &daughterVertex, const float minLongitudinalDistance, const float maxTransverseDistance)
Whether pointing vertex is adjacent to a given position.
void GetConnectedElements(const pandora::Cluster *const pCluster, const bool ignoreUnavailable, pandora::ClusterList &clusterListU, pandora::ClusterList &clusterListV, pandora::ClusterList &clusterListW) const
Get elements connected to a specified cluster.