9 #include "Pandora/AlgorithmHeaders.h" 22 VertexBasedPfoRecoveryAlgorithm::VertexBasedPfoRecoveryAlgorithm() :
23 m_slidingFitHalfWindow(10),
24 m_maxLongitudinalDisplacement(5.
f),
25 m_maxTransverseDisplacement(2.
f),
26 m_twoViewChi2Cut(5.
f),
27 m_threeViewChi2Cut(5.
f)
36 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetCurrentList(*
this, pVertexList));
38 const Vertex *
const pSelectedVertex(
39 (pVertexList && (pVertexList->size() == 1) && (VERTEX_3D == (*(pVertexList->begin()))->GetVertexType())) ? *(pVertexList->begin()) : NULL);
43 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
44 std::cout <<
"VertexBasedPfoRecoveryAlgorithm: unable to find vertex in current list " << std::endl;
46 return STATUS_CODE_SUCCESS;
59 this->
SelectVertexClusters(pSelectedVertex, slidingFitResultMap, availableClusters, selectedClusters);
64 this->
MatchThreeViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
65 this->
MatchTwoViews(pSelectedVertex, slidingFitResultMap, selectedClusters, vetoList, particleList);
70 return STATUS_CODE_SUCCESS;
79 const ClusterList *pClusterList(NULL);
80 PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*
this, *iter, pClusterList));
82 if (NULL == pClusterList)
84 if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
85 std::cout <<
"VertexBasedPfoRecoveryAlgorithm: could not find cluster list " << *iter << std::endl;
91 const Cluster *
const pCluster = *cIter;
93 if (!pCluster->IsAvailable())
96 if (pCluster->GetNCaloHits() <= 1)
102 clusterVector.push_back(pCluster);
108 return STATUS_CODE_SUCCESS;
116 if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
124 if (pointingCluster.
GetLengthSquared() < std::numeric_limits<float>::epsilon())
127 if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
128 throw StatusCodeException(STATUS_CODE_FAILURE);
130 catch (StatusCodeException &statusCodeException)
132 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
133 throw statusCodeException;
150 const Cluster *
const pCluster = *cIter;
153 if (TPC_3D == hitType)
156 const CartesianVector vertexPosition((TPC_VIEW_U == hitType) ? vertexU : (TPC_VIEW_V == hitType) ? vertexV : vertexW);
159 if (slidingFitResultMap.end() == sIter)
165 for (
unsigned int iVtx = 0; iVtx < 2; ++iVtx)
169 float rL(0.
f), rT(0.
f);
174 outputClusters.push_back(pCluster);
188 ClusterVector availableClusters, clustersU, clustersV, clustersW;
195 const Cluster *pCluster1(NULL);
196 const Cluster *pCluster2(NULL);
197 const Cluster *pCluster3(NULL);
199 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, clustersW, pCluster1, pCluster2, pCluster3, chi2);
201 if (NULL == pCluster1 || NULL == pCluster2 || NULL == pCluster3)
208 const Cluster *
const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1
209 : (TPC_VIEW_U == hitType2) ? pCluster2
210 : (TPC_VIEW_U == hitType3) ? pCluster3
212 const Cluster *
const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1
213 : (TPC_VIEW_V == hitType2) ? pCluster2
214 : (TPC_VIEW_V == hitType3) ? pCluster3
216 const Cluster *
const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1
217 : (TPC_VIEW_W == hitType2) ? pCluster2
218 : (TPC_VIEW_W == hitType3) ? pCluster3
221 particleList.push_back(
Particle(pClusterU, pClusterV, pClusterW));
223 vetoList.insert(pCluster1);
224 vetoList.insert(pCluster2);
225 vetoList.insert(pCluster3);
236 ClusterVector availableClusters, clustersU, clustersV, clustersW;
243 const Cluster *pCluster1(NULL);
244 const Cluster *pCluster2(NULL);
246 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, pCluster1, pCluster2, chi2);
247 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersV, clustersW, pCluster1, pCluster2, chi2);
248 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersW, clustersU, pCluster1, pCluster2, chi2);
250 if (NULL == pCluster1 || NULL == pCluster2)
256 const Cluster *
const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
257 const Cluster *
const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
258 const Cluster *
const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
260 particleList.push_back(
Particle(pClusterU, pClusterV, pClusterW));
262 vetoList.insert(pCluster1);
263 vetoList.insert(pCluster2);
271 const Cluster *&pBestCluster2,
const Cluster *&pBestCluster3,
float &bestChi2)
const 273 if (clusters1.empty() || clusters2.empty() || clusters3.empty())
279 const Cluster *
const pCluster1 = *cIter1;
282 if (slidingFitResultMap.end() == sIter1)
291 const Cluster *
const pCluster2 = *cIter2;
294 if (slidingFitResultMap.end() == sIter2)
303 const Cluster *
const pCluster3 = *cIter3;
306 if (slidingFitResultMap.end() == sIter3)
313 const float thisChi2(this->
GetChi2(pVertex, pointingCluster1, pointingCluster2, pointingCluster3));
315 if (thisChi2 < bestChi2)
318 pBestCluster1 = pCluster1;
319 pBestCluster2 = pCluster2;
320 pBestCluster3 = pCluster3;
330 const ClusterVector &clusters1,
const ClusterVector &clusters2,
const Cluster *&pBestCluster1,
const Cluster *&pBestCluster2,
float &bestChi2)
const 332 if (clusters1.empty() || clusters2.empty())
338 const Cluster *
const pCluster1 = *cIter1;
341 if (slidingFitResultMap.end() == sIter1)
350 const Cluster *
const pCluster2 = *cIter2;
353 if (slidingFitResultMap.end() == sIter2)
360 const float thisChi2(this->
GetChi2(pVertex, pointingCluster1, pointingCluster2));
362 if (thisChi2 < bestChi2)
365 pBestCluster1 = pCluster1;
366 pBestCluster2 = pCluster2;
380 if (hitType1 == hitType2)
381 throw StatusCodeException(STATUS_CODE_FAILURE);
390 CartesianVector mergedPosition(0.
f, 0.
f, 0.
f);
392 this->GetPandora(), hitType1, hitType2, pointingVertex1.
GetPosition(), pointingVertex2.
GetPosition(), mergedPosition, chi2);
406 if ((hitType1 == hitType2) || (hitType2 == hitType3) || (hitType3 == hitType1))
407 throw StatusCodeException(STATUS_CODE_FAILURE);
418 CartesianVector mergedPosition(0.
f, 0.
f, 0.
f);
431 if (0 == vetoList.count(*iter))
432 outputVector.push_back(*iter);
443 outputVector.push_back(*iter);
454 if (innerDistance < outerDistance)
476 if (particleList.empty())
479 const PfoList *pPfoList = NULL;
480 std::string pfoListName;
481 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pPfoList, pfoListName));
487 ClusterList clusterList;
488 const Cluster *
const pClusterU = particle.
m_pClusterU;
489 const Cluster *
const pClusterV = particle.
m_pClusterV;
490 const Cluster *
const pClusterW = particle.
m_pClusterW;
492 const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() :
true);
493 const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() :
true);
494 const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() :
true);
496 if (!(isAvailableU && isAvailableV && isAvailableW))
497 throw StatusCodeException(STATUS_CODE_FAILURE);
500 clusterList.push_back(pClusterU);
502 clusterList.push_back(pClusterV);
504 clusterList.push_back(pClusterW);
507 PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
508 pfoParameters.m_particleId = MU_MINUS;
509 pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
510 pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
511 pfoParameters.m_energy = 0.f;
512 pfoParameters.m_momentum = CartesianVector(0.
f, 0.
f, 0.
f);
513 pfoParameters.m_clusterList = clusterList;
515 const ParticleFlowObject *pPfo(NULL);
516 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*
this, pfoParameters, pPfo));
519 if (!pPfoList->empty())
520 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*
this,
m_outputPfoListName));
526 m_pClusterU(pClusterU),
527 m_pClusterV(pClusterV),
528 m_pClusterW(pClusterW)
531 throw StatusCodeException(STATUS_CODE_FAILURE);
537 if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
538 throw StatusCodeException(STATUS_CODE_FAILURE);
545 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
"InputClusterListNames",
m_inputClusterListNames));
547 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputPfoListName",
m_outputPfoListName));
549 PANDORA_RETURN_RESULT_IF_AND_IF(
550 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitHalfWindow",
m_slidingFitHalfWindow));
552 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
555 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
558 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"TwoViewChi2Cut",
m_twoViewChi2Cut));
560 PANDORA_RETURN_RESULT_IF_AND_IF(
561 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ThreeViewChi2Cut",
m_threeViewChi2Cut));
563 return STATUS_CODE_SUCCESS;
std::string m_outputPfoListName
The name of the output pfo list.
pandora::StatusCode Run()
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 LArPointingCluster::Vertex & GetInnerVertex(const pandora::CartesianVector &vertex, const LArPointingCluster &cluster) const
Find nearest end of pointing cluster to a specified position vector.
void GetBestChi2(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &clusters1, const pandora::ClusterVector &clusters2, const pandora::ClusterVector &clusters3, const pandora::Cluster *&pBestCluster1, const pandora::Cluster *&pBestCluster2, const pandora::Cluster *&pBestCluster3, float &chi2) const
Get best-matched triplet of clusters from a set of input cluster vectors.
static void MergeTwoPositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from two views to give unified 3D position.
std::vector< Particle > ParticleList
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
float m_maxLongitudinalDisplacement
Particle(const pandora::Cluster *const pClusterU, const pandora::Cluster *const pClusterV, const pandora::Cluster *const pClusterW)
Constructor.
void SelectVertexClusters(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &inputClusters, pandora::ClusterVector &outputClusters) const
Select clusters in proximity to reconstructed vertex.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void MatchTwoViews(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &selectedClusters, pandora::ClusterSet &vetoList, ParticleList &particleList) const
Match clusters from two views.
static pandora::CartesianVector ProjectPosition(const pandora::Pandora &pandora, const pandora::CartesianVector &position3D, const pandora::HitType view)
Project 3D position into a given 2D view.
LArPointingCluster class.
Cluster finding and building.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
void MatchThreeViews(const pandora::Vertex *const pVertex, const TwoDSlidingFitResultMap &slidingFitResultMap, const pandora::ClusterVector &selectedClusters, pandora::ClusterSet &vetoList, ParticleList &particleList) const
Match clusters from three views.
const pandora::Cluster * GetCluster() const
Get the address of the cluster.
Header file for the geometry helper class.
pandora::StringVector m_inputClusterListNames
The list of input cluster list names.
Header file for the cluster helper class.
const Vertex & GetOuterVertex() const
Get the outer vertex.
const Vertex & GetInnerVertex() const
Get the inner vertex.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
unsigned int m_slidingFitHalfWindow
static void MergeThreePositions3D(const pandora::Pandora &pandora, const pandora::HitType view1, const pandora::HitType view2, const pandora::HitType view3, const pandora::CartesianVector &position1, const pandora::CartesianVector &position2, const pandora::CartesianVector &position3, pandora::CartesianVector &position3D, float &chiSquared)
Merge 2D positions from three views to give unified 3D position.
float GetChi2(const pandora::Vertex *const pVertex, const LArPointingCluster &pointingCluster1, const LArPointingCluster &pointingCluster2) const
Merge two pointing clusters and return chi-squared metric giving consistency of matching.
float GetLengthSquared() const
Get length squared of pointing cluster.
pandora::StatusCode GetAvailableClusters(const pandora::StringVector inputClusterListName, pandora::ClusterVector &clusterVector) const
Get a vector of available clusters.
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
const LArPointingCluster::Vertex & GetOuterVertex(const pandora::CartesianVector &vertex, const LArPointingCluster &cluster) const
Find furthest end of pointing cluster from a specified position vector.
void SelectAvailableClusters(const pandora::ClusterSet &vetoList, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select cluster which haven't been vetoed.
const pandora::Cluster * m_pClusterV
Address of cluster in V view.
float m_maxTransverseDisplacement
void SelectClusters(const pandora::HitType hitType, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select clusters of a specified hit type.
boost::graph_traits< ModuleGraph >::vertex_descriptor Vertex
bool IsInnerVertex() const
Is this the inner vertex.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
Header file for the vertex-based particle recovery algorithm.
const pandora::Cluster * m_pClusterW
Address of cluster in W view.
std::list< Vertex > VertexList
const pandora::CartesianVector & GetPosition() const
Get the vertex position.
void BuildParticles(const ParticleList &particleList)
Build particle flow objects from matched clusters.
TwoDSlidingFitResult class.
const pandora::Cluster * m_pClusterU
Address of cluster in U view.