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;
77 for (StringVector::const_iterator iter = inputClusterListNames.begin(), iterEnd = inputClusterListNames.end(); iter != iterEnd; ++iter)
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;
89 for (ClusterList::const_iterator cIter = pClusterList->begin(), cIterEnd = pClusterList->end(); cIter != cIterEnd; ++cIter)
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 for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
118 if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
125 if (pointingCluster.
GetLengthSquared() < std::numeric_limits<float>::epsilon())
128 if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
129 throw StatusCodeException(STATUS_CODE_FAILURE);
131 catch (StatusCodeException &statusCodeException)
133 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
134 throw statusCodeException;
149 for (ClusterVector::const_iterator cIter = inputClusters.begin(), cIterEnd = inputClusters.end(); cIter != cIterEnd; ++cIter)
151 const Cluster *
const pCluster = *cIter;
154 if (TPC_3D == hitType)
157 const CartesianVector vertexPosition((TPC_VIEW_U == hitType) ? vertexU : (TPC_VIEW_V == hitType) ? vertexV : vertexW);
159 TwoDSlidingFitResultMap::const_iterator sIter = slidingFitResultMap.find(pCluster);
160 if (slidingFitResultMap.end() == sIter)
166 for (
unsigned int iVtx = 0; iVtx < 2; ++iVtx)
170 float rL(0.f), rT(0.f);
175 outputClusters.push_back(pCluster);
189 ClusterVector availableClusters, clustersU, clustersV, clustersW;
196 const Cluster *pCluster1(NULL);
197 const Cluster *pCluster2(NULL);
198 const Cluster *pCluster3(NULL);
200 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, clustersW, pCluster1, pCluster2, pCluster3, chi2);
202 if (NULL == pCluster1 || NULL == pCluster2 || NULL == pCluster3)
209 const Cluster *
const pClusterU(
210 (TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : (TPC_VIEW_U == hitType3) ? pCluster3 : NULL);
211 const Cluster *
const pClusterV(
212 (TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : (TPC_VIEW_V == hitType3) ? pCluster3 : NULL);
213 const Cluster *
const pClusterW(
214 (TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : (TPC_VIEW_W == hitType3) ? pCluster3 : NULL);
216 particleList.push_back(
Particle(pClusterU, pClusterV, pClusterW));
218 vetoList.insert(pCluster1);
219 vetoList.insert(pCluster2);
220 vetoList.insert(pCluster3);
231 ClusterVector availableClusters, clustersU, clustersV, clustersW;
238 const Cluster *pCluster1(NULL);
239 const Cluster *pCluster2(NULL);
241 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersU, clustersV, pCluster1, pCluster2, chi2);
242 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersV, clustersW, pCluster1, pCluster2, chi2);
243 this->
GetBestChi2(pVertex, slidingFitResultMap, clustersW, clustersU, pCluster1, pCluster2, chi2);
245 if (NULL == pCluster1 || NULL == pCluster2)
251 const Cluster *
const pClusterU((TPC_VIEW_U == hitType1) ? pCluster1 : (TPC_VIEW_U == hitType2) ? pCluster2 : NULL);
252 const Cluster *
const pClusterV((TPC_VIEW_V == hitType1) ? pCluster1 : (TPC_VIEW_V == hitType2) ? pCluster2 : NULL);
253 const Cluster *
const pClusterW((TPC_VIEW_W == hitType1) ? pCluster1 : (TPC_VIEW_W == hitType2) ? pCluster2 : NULL);
255 particleList.push_back(
Particle(pClusterU, pClusterV, pClusterW));
257 vetoList.insert(pCluster1);
258 vetoList.insert(pCluster2);
266 const Cluster *&pBestCluster2,
const Cluster *&pBestCluster3,
float &bestChi2)
const
268 if (clusters1.empty() || clusters2.empty() || clusters3.empty())
272 for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
274 const Cluster *
const pCluster1 = *cIter1;
276 TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
277 if (slidingFitResultMap.end() == sIter1)
284 for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
286 const Cluster *
const pCluster2 = *cIter2;
288 TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
289 if (slidingFitResultMap.end() == sIter2)
296 for (ClusterVector::const_iterator cIter3 = clusters3.begin(), cIterEnd3 = clusters3.end(); cIter3 != cIterEnd3; ++cIter3)
298 const Cluster *
const pCluster3 = *cIter3;
300 TwoDSlidingFitResultMap::const_iterator sIter3 = slidingFitResultMap.find(pCluster3);
301 if (slidingFitResultMap.end() == sIter3)
308 const float thisChi2(this->
GetChi2(pVertex, pointingCluster1, pointingCluster2, pointingCluster3));
310 if (thisChi2 < bestChi2)
313 pBestCluster1 = pCluster1;
314 pBestCluster2 = pCluster2;
315 pBestCluster3 = pCluster3;
325 const ClusterVector &clusters1,
const ClusterVector &clusters2,
const Cluster *&pBestCluster1,
const Cluster *&pBestCluster2,
float &bestChi2)
const
327 if (clusters1.empty() || clusters2.empty())
331 for (ClusterVector::const_iterator cIter1 = clusters1.begin(), cIterEnd1 = clusters1.end(); cIter1 != cIterEnd1; ++cIter1)
333 const Cluster *
const pCluster1 = *cIter1;
335 TwoDSlidingFitResultMap::const_iterator sIter1 = slidingFitResultMap.find(pCluster1);
336 if (slidingFitResultMap.end() == sIter1)
339 const TwoDSlidingFitResult &slidingFitResult1 = sIter1->second;
340 const LArPointingCluster pointingCluster1(slidingFitResult1);
343 for (ClusterVector::const_iterator cIter2 = clusters2.begin(), cIterEnd2 = clusters2.end(); cIter2 != cIterEnd2; ++cIter2)
345 const Cluster *
const pCluster2 = *cIter2;
347 TwoDSlidingFitResultMap::const_iterator sIter2 = slidingFitResultMap.find(pCluster2);
348 if (slidingFitResultMap.end() == sIter2)
351 const TwoDSlidingFitResult &slidingFitResult2 = sIter2->second;
352 const LArPointingCluster pointingCluster2(slidingFitResult2);
355 const float thisChi2(this->
GetChi2(pVertex, pointingCluster1, pointingCluster2));
357 if (thisChi2 < bestChi2)
360 pBestCluster1 = pCluster1;
361 pBestCluster2 = pCluster2;
370 const Vertex *
const pVertex,
const LArPointingCluster &pointingCluster1,
const LArPointingCluster &pointingCluster2)
const
375 if (hitType1 == hitType2)
376 throw StatusCodeException(STATUS_CODE_FAILURE);
381 const LArPointingCluster::Vertex &pointingVertex1(this->
GetOuterVertex(vertex1, pointingCluster1));
382 const LArPointingCluster::Vertex &pointingVertex2(this->
GetOuterVertex(vertex2, pointingCluster2));
385 CartesianVector mergedPosition(0.f, 0.f, 0.f);
387 this->GetPandora(), hitType1, hitType2, pointingVertex1.GetPosition(), pointingVertex2.GetPosition(), mergedPosition, chi2);
395 const LArPointingCluster &pointingCluster2,
const LArPointingCluster &pointingCluster3)
const
401 if ((hitType1 == hitType2) || (hitType2 == hitType3) || (hitType3 == hitType1))
402 throw StatusCodeException(STATUS_CODE_FAILURE);
408 const LArPointingCluster::Vertex &pointingVertex1(this->
GetOuterVertex(vertex1, pointingCluster1));
409 const LArPointingCluster::Vertex &pointingVertex2(this->
GetOuterVertex(vertex2, pointingCluster2));
410 const LArPointingCluster::Vertex &pointingVertex3(this->
GetOuterVertex(vertex3, pointingCluster3));
413 CartesianVector mergedPosition(0.f, 0.f, 0.f);
415 pointingVertex2.GetPosition(), pointingVertex3.GetPosition(), mergedPosition, chi2);
424 for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
426 if (0 == vetoList.count(*iter))
427 outputVector.push_back(*iter);
435 for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
438 outputVector.push_back(*iter);
449 if (innerDistance < outerDistance)
471 if (particleList.empty())
474 const PfoList *pPfoList = NULL;
475 std::string pfoListName;
476 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*
this, pPfoList, pfoListName));
478 for (ParticleList::const_iterator iter = particleList.begin(), iterEnd = particleList.end(); iter != iterEnd; ++iter)
482 ClusterList clusterList;
483 const Cluster *
const pClusterU = particle.
m_pClusterU;
484 const Cluster *
const pClusterV = particle.
m_pClusterV;
485 const Cluster *
const pClusterW = particle.
m_pClusterW;
487 const bool isAvailableU((NULL != pClusterU) ? pClusterU->IsAvailable() :
true);
488 const bool isAvailableV((NULL != pClusterV) ? pClusterV->IsAvailable() :
true);
489 const bool isAvailableW((NULL != pClusterW) ? pClusterW->IsAvailable() :
true);
491 if (!(isAvailableU && isAvailableV && isAvailableW))
492 throw StatusCodeException(STATUS_CODE_FAILURE);
495 clusterList.push_back(pClusterU);
497 clusterList.push_back(pClusterV);
499 clusterList.push_back(pClusterW);
503 pfoParameters.m_particleId = MU_MINUS;
504 pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
505 pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
506 pfoParameters.m_energy = 0.f;
507 pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
508 pfoParameters.m_clusterList = clusterList;
510 const ParticleFlowObject *pPfo(NULL);
511 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*
this, pfoParameters, pPfo));
514 if (!pPfoList->empty())
515 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*
this,
m_outputPfoListName));
521 m_pClusterU(pClusterU),
522 m_pClusterV(pClusterV),
523 m_pClusterW(pClusterW)
526 throw StatusCodeException(STATUS_CODE_FAILURE);
532 if (!(TPC_VIEW_U == hitTypeU && TPC_VIEW_V == hitTypeV && TPC_VIEW_W == hitTypeW))
533 throw StatusCodeException(STATUS_CODE_FAILURE);
540 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
"InputClusterListNames",
m_inputClusterListNames));
542 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
"OutputPfoListName",
m_outputPfoListName));
544 PANDORA_RETURN_RESULT_IF_AND_IF(
545 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitHalfWindow",
m_slidingFitHalfWindow));
547 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
550 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
553 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"TwoViewChi2Cut",
m_twoViewChi2Cut));
555 PANDORA_RETURN_RESULT_IF_AND_IF(
556 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ThreeViewChi2Cut",
m_threeViewChi2Cut));
558 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.
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.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
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.
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.
fhicl::Table< sbnd::crt::CRTDetSimParams > Parameters
void SelectAvailableClusters(const pandora::ClusterSet &vetoList, const pandora::ClusterVector &inputVector, pandora::ClusterVector &outputVector) const
Select cluster which haven't been vetoed.
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
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.
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.
BEGIN_PROLOG could also be cout
const pandora::Cluster * m_pClusterU
Address of cluster in U view.