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 : (TPC_VIEW_V == hitType) ? inputClusterListV : inputClusterListW);
89 clusterList.push_back(pCluster);
100 ClusterList vertexClusterList;
116 const Cluster *
const pCluster = *iter;
118 if (!pCluster->IsAvailable())
127 float xMin(0.
f), xMax(0.
f);
133 selectedClusterList.push_back(pCluster);
141 CartesianPointVector vertexList;
147 if (!(*iter)->IsAvailable())
154 catch (StatusCodeException &) {}
161 const Cluster *
const pCluster = *iter;
163 if (!pCluster->IsAvailable())
173 selectedClusterList.push_back(pCluster);
178 catch (StatusCodeException &) {}
201 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
203 if ((0 == pCluster1->GetNCaloHits()) || (0 == pCluster2->GetNCaloHits()))
204 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
206 float xMin1(0.
f), xMax1(0.
f), xMin2(0.
f), xMax2(0.
f);
210 const float xSpan1(xMax1 - xMin1), xSpan2(xMax2 - xMin2);
212 if ((xSpan1 < std::numeric_limits<float>::epsilon()) || (xSpan2 < std::numeric_limits<float>::epsilon()))
213 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
217 float xOverlapFraction1(xOverlap / xSpan1), xOverlapFraction2(xOverlap / xSpan2);
231 const Cluster *
const pCluster2,
const float xMin2,
const float xMax2,
float &xOverlapFraction1,
float &xOverlapFraction2)
const 233 if (PandoraContentApi::GetGeometry(*this)->GetDetectorGapList().empty())
236 const float xMin(
std::min(xMin1, xMin2));
237 const float xMax(
std::max(xMax1, xMax2));
238 float xMinEff1(xMin1), xMaxEff1(xMax1), xMinEff2(xMin2), xMaxEff2(xMax2);
243 const float effectiveXSpan1(xMaxEff1 - xMinEff1), effectiveXSpan2(xMaxEff2 - xMinEff2);
244 const float effectiveXOverlapSpan(
std::min(xMaxEff1, xMaxEff2) -
std::max(xMinEff1, xMinEff2));
246 if ((effectiveXSpan1 > std::numeric_limits<float>::epsilon()) && (effectiveXSpan2 > std::numeric_limits<float>::epsilon()))
248 xOverlapFraction1 =
std::min(1.
f, (effectiveXOverlapSpan / effectiveXSpan1));
249 xOverlapFraction2 =
std::min(1.
f, (effectiveXOverlapSpan / effectiveXSpan2));
264 const int nSamplingPointsLeft(1 + static_cast<int>((xMinEff - xMin) /
m_sampleStepSize));
265 const int nSamplingPointsRight(1 + static_cast<int>((xMax - xMaxEff) /
m_sampleStepSize));
266 float dxMin(0.
f), dxMax(0.
f);
268 for (
int iSample = 1; iSample <= nSamplingPointsLeft; ++iSample)
275 dxMin = xMinEff - xSample;
278 for (
int iSample = 1; iSample <= nSamplingPointsRight; ++iSample)
285 dxMax = xSample - xMaxEff;
288 xMinEff = xMinEff - dxMin;
289 xMaxEff = xMaxEff + dxMax;
291 catch (StatusCodeException &) {}
301 for (
const Cluster *
const pKeyCluster : sortedKeyClusters)
303 ClusterList clusterListU, clusterListV, clusterListW;
305 overlapTensor.
GetConnectedElements(pKeyCluster,
true, clusterListU, clusterListV, clusterListW);
306 const unsigned int nU(clusterListU.size()), nV(clusterListV.size()), nW(clusterListW.size());
308 if ((0 == nU * nV) && (0 == nV * nW) && (0 == nW * nU))
311 ClusterList clusterListAll;
312 clusterListAll.insert(clusterListAll.end(), clusterListU.begin(), clusterListU.end());
313 clusterListAll.insert(clusterListAll.end(), clusterListV.begin(), clusterListV.end());
314 clusterListAll.insert(clusterListAll.end(), clusterListW.begin(), clusterListW.end());
316 if ((1 == nU * nV * nW) && this->
CheckConsistency(*(clusterListU.begin()), *(clusterListV.begin()), *(clusterListW.begin())))
320 else if ((0 == nU * nV * nW) && ((1 == nU && 1 == nV) || (1 == nV && 1 == nW) || (1 == nW && 1 == nU)))
337 float xMinU(0.
f), xMinV(0.
f), xMinW(0.
f), xMaxU(0.
f), xMaxV(0.
f), xMaxW(0.
f);
344 const float xOverlap(xMax - xMin);
346 if (xOverlap < std::numeric_limits<float>::epsilon())
363 const float pseudoChi2(((u - vw2u) * (u - vw2u) + (v - wu2v) * (v - wu2v) + (w - uv2w) * (w - uv2w)) / 3.
f);
375 if (clusterList.size() < 2)
376 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
378 const PfoList *pPfoList = NULL; std::string pfoListName;
379 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pPfoList, pfoListName));
382 PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
383 pfoParameters.m_particleId = MU_MINUS;
384 pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
385 pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
386 pfoParameters.m_energy = 0.f;
387 pfoParameters.m_momentum = CartesianVector(0.
f, 0.
f, 0.
f);
388 pfoParameters.m_clusterList = clusterList;
390 const ParticleFlowObject *pPfo(NULL);
391 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*
this, pfoParameters, pPfo));
393 if (!pPfoList->empty())
395 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*
this,
m_outputPfoListName));
396 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*
this,
m_outputPfoListName));
408 const Cluster *
const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
409 const Cluster *
const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
410 const Cluster *
const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
412 if (pClusterU && pClusterV && !pClusterW)
414 m_clusterNavigationMapUV[pClusterU].push_back(pClusterV);
416 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterU))
417 m_keyClusters.push_back(pClusterU);
419 else if (!pClusterU && pClusterV && pClusterW)
421 m_clusterNavigationMapVW[pClusterV].push_back(pClusterW);
423 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterV))
424 m_keyClusters.push_back(pClusterV);
426 else if (pClusterU && !pClusterV && pClusterW)
428 m_clusterNavigationMapWU[pClusterW].push_back(pClusterU);
430 if (m_keyClusters.end() == std::find(m_keyClusters.begin(), m_keyClusters.end(), pClusterW))
431 m_keyClusters.push_back(pClusterW);
435 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
442 ClusterList &clusterListU, ClusterList &clusterListV, ClusterList &clusterListW)
const 444 if (ignoreUnavailable && !pCluster->IsAvailable())
449 if (!((TPC_VIEW_U == hitType) || (TPC_VIEW_V == hitType) || (TPC_VIEW_W == hitType)))
450 throw StatusCodeException(STATUS_CODE_FAILURE);
452 ClusterList &clusterList((TPC_VIEW_U == hitType) ? clusterListU : (TPC_VIEW_V == hitType) ? clusterListV : clusterListW);
453 const ClusterNavigationMap &navigationMap((TPC_VIEW_U == hitType) ? m_clusterNavigationMapUV : (TPC_VIEW_V == hitType) ? m_clusterNavigationMapVW : m_clusterNavigationMapWU);
455 if (clusterList.end() != std::find(clusterList.begin(), clusterList.end(), pCluster))
458 clusterList.push_back(pCluster);
462 if (navigationMap.end() == iter)
466 this->GetConnectedElements(*cIter, ignoreUnavailable, clusterListU, clusterListV, clusterListW);
474 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
"InputClusterListNames",
m_inputClusterListNames));
475 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputPfoListName",
m_outputPfoListName));
477 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
480 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
484 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
485 "MinClusterLength", minClusterLength));
488 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
491 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
494 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
497 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
500 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
503 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
506 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
511 std::cout <<
"ParticleRecoveryAlgorithm: Invalid value for SampleStepSize " <<
m_sampleStepSize << std::endl;
512 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
515 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
518 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
521 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.
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.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
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).
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 ...
std::vector< art::Ptr< recob::Cluster > > ClusterVector
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.
std::unordered_map< const pandora::Cluster *, pandora::ClusterList > ClusterNavigationMap
static void GetClusterSpanX(const pandora::Cluster *const pCluster, float &xmin, float &xmax)
Get minimum and maximum X positions of the calo hits in a cluster.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
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.