StitchingCosmicRayMergingTool class. More...

#include <StitchingCosmicRayMergingTool.h>

Inheritance diagram for lar_content::StitchingCosmicRayMergingTool:

Classes
class	PfoAssociation
	PfoAssociation class. More...

Public Member Functions
	StitchingCosmicRayMergingTool ()
	Default constructor. More...

void	Run (const MasterAlgorithm const pAlgorithm, const pandora::PfoList const pMultiPfoList, PfoToLArTPCMap &pfoToLArTPCMap, PfoToFloatMap &stitchedPfosToX0Map)
	Run the algorithm tool. More...

Private Types
typedef std::unordered_map < const pandora::ParticleFlowObject *, LArPointingCluster >	ThreeDPointingClusterMap

typedef std::unordered_map < const pandora::LArTPC *, pandora::PfoList >	LArTPCToPfoMap

typedef std::unordered_map < const pandora::ParticleFlowObject *, PfoAssociation >	PfoAssociationMap

typedef std::unordered_map < const pandora::ParticleFlowObject *, PfoAssociationMap >	PfoAssociationMatrix

typedef std::unordered_map < const pandora::ParticleFlowObject *, pandora::PfoList >	PfoMergeMap

typedef std::unordered_map < const pandora::ParticleFlowObject *, LArPointingCluster::Vertex >	PfoToPointingVertexMap

typedef std::unordered_map < const pandora::ParticleFlowObject *, PfoToPointingVertexMap >	PfoToPointingVertexMatrix

Private Member Functions
pandora::StatusCode	ReadSettings (const pandora::TiXmlHandle xmlHandle)

void	SelectPrimaryPfos (const pandora::PfoList *pInputPfoList, const PfoToLArTPCMap &pfoToLArTPCMap, pandora::PfoList &outputPfoList) const
	Select primary Pfos from the input list of Pfos. More...

void	BuildPointingClusterMaps (const pandora::PfoList &inputPfoList, const PfoToLArTPCMap &pfoToLArTPCMap, ThreeDPointingClusterMap &pointingClusterMap) const
	Build a 3D pointing cluster for each Pfo. More...

void	BuildTPCMaps (const pandora::PfoList &inputPfoList, const PfoToLArTPCMap &pfoToLArTPCMap, LArTPCToPfoMap &larTPCToPfoMap) const
	Build a list of Pfos for each tpc. More...

void	CreatePfoMatches (const LArTPCToPfoMap &larTPCToPfoMap, const ThreeDPointingClusterMap &pointingClusterMap, PfoAssociationMatrix &pfoAssociationMatrix) const
	Create associations between Pfos using 3D pointing clusters. More...

void	CreatePfoMatches (const pandora::LArTPC &larTPC1, const pandora::LArTPC &larTPC2, const pandora::ParticleFlowObject const pPfo1, const pandora::ParticleFlowObject const pPfo2, const ThreeDPointingClusterMap &pointingClusterMap, PfoAssociationMatrix &pfoAssociationMatrix) const
	Create associations between Pfos using 3D pointing clusters. More...

void	SelectPfoMatches (const PfoAssociationMatrix &pfoAssociationMatrix, PfoMergeMap &pfoSelectedMatches) const
	Select the best associations between Pfos; create a mapping between associated Pfos, handling any ambiguities. More...

void	SelectPfoMerges (const PfoMergeMap &pfoMatches, PfoMergeMap &pfoMerges) const
	Create an initial map of Pfo merges to be made. More...

void	CollectAssociatedPfos (const pandora::ParticleFlowObject const pSeedPfo, const pandora::ParticleFlowObject const pCurrentPfo, const PfoMergeMap &pfoMerges, const pandora::PfoSet &vetoSet, pandora::PfoList &associatedList) const
	Collect up associations between Pfos. More...

void	OrderPfoMerges (const PfoToLArTPCMap &pfoToLArTPCMap, const ThreeDPointingClusterMap &pointingClusterMap, const PfoMergeMap &inputPfoMerges, PfoMergeMap &outputPfoMerges) const
	Identify the vertex Pfo and then re-order the map of merges so that the vertex Pfo will be enlarged. More...

void	StitchPfos (const MasterAlgorithm *const pAlgorithm, const ThreeDPointingClusterMap &pointingClusterMap, const PfoMergeMap &pfoMerges, PfoToLArTPCMap &pfoToLArTPCMap, PfoToFloatMap &stitchedPfosToX0Map) const
	Apply X0 corrections, and then stitch together Pfos. More...

void	ShiftPfo (const MasterAlgorithm const pAlgorithm, const pandora::ParticleFlowObject const pPfoToShift, const pandora::ParticleFlowObject *const pMatchedPfo, const float x0, const PfoToLArTPCMap &pfoToLArTPCMap, const PfoToPointingVertexMatrix &pfoToPointingVertexMatrix) const
	Shift a pfo given its pfo stitching pair. More...

bool	CalculateX0 (const PfoToLArTPCMap &pfoToLArTPCMap, const ThreeDPointingClusterMap &pointingClusterMap, const pandora::PfoVector &pfoVector, float &x0, PfoToPointingVertexMatrix &pfoToPointingVertexMatrix) const
	Calculate x0 shift for a group of associated Pfos. More...

Private Attributes
bool	m_useXcoordinate

bool	m_alwaysApplyT0Calculation

int	m_halfWindowLayers

float	m_minLengthSquared

float	m_minCosRelativeAngle

float	m_relaxMinLongitudinalDisplacement
	The minimum value of the longitudinal impact parameter for association if both verticies fall in the detector gap. More...

float	m_maxLongitudinalDisplacementX

float	m_maxTransverseDisplacement

float	m_relaxCosRelativeAngle

float	m_relaxTransverseDisplacement

unsigned int	m_minNCaloHits3D

float	m_maxX0FractionalDeviation
	The maximum allowed fractional difference of an X0 contribution for matches to be stitched. More...

float	m_boundaryToleranceWidth
	The distance from the APA/CPA boundary inside which the deviation consideration is ignored. More...

Detailed Description

StitchingCosmicRayMergingTool class.

Definition at line 23 of file StitchingCosmicRayMergingTool.h.

Member Typedef Documentation

typedef std::unordered_map<const pandora::LArTPC *, pandora::PfoList> lar_content::StitchingCosmicRayMergingTool::LArTPCToPfoMap

private

Definition at line 110 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, PfoAssociation> lar_content::StitchingCosmicRayMergingTool::PfoAssociationMap

private

Definition at line 121 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, PfoAssociationMap> lar_content::StitchingCosmicRayMergingTool::PfoAssociationMatrix

private

Definition at line 122 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, pandora::PfoList> lar_content::StitchingCosmicRayMergingTool::PfoMergeMap

private

Definition at line 148 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, LArPointingCluster::Vertex> lar_content::StitchingCosmicRayMergingTool::PfoToPointingVertexMap

private

Definition at line 201 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, PfoToPointingVertexMap> lar_content::StitchingCosmicRayMergingTool::PfoToPointingVertexMatrix

private

Definition at line 202 of file StitchingCosmicRayMergingTool.h.

typedef std::unordered_map<const pandora::ParticleFlowObject *, LArPointingCluster> lar_content::StitchingCosmicRayMergingTool::ThreeDPointingClusterMap

private

Definition at line 98 of file StitchingCosmicRayMergingTool.h.

Constructor & Destructor Documentation

lar_content::StitchingCosmicRayMergingTool::StitchingCosmicRayMergingTool ( )

Default constructor.

Definition at line 25 of file StitchingCosmicRayMergingTool.cc.

                                                              :
     m_useXcoordinate(false),
     m_alwaysApplyT0Calculation(true),
     m_halfWindowLayers(30),
     m_minLengthSquared(50.f),
     m_minCosRelativeAngle(0.966),
     m_relaxMinLongitudinalDisplacement(-5.f),
     m_maxLongitudinalDisplacementX(15.f),
     m_maxTransverseDisplacement(5.f),
     m_relaxCosRelativeAngle(0.906),
     m_relaxTransverseDisplacement(2.5f),
     m_minNCaloHits3D(0),
     m_maxX0FractionalDeviation(0.3f),
     m_boundaryToleranceWidth(10.f)
 {
 }

Member Function Documentation

void lar_content::StitchingCosmicRayMergingTool::BuildPointingClusterMaps	(	const pandora::PfoList &	inputPfoList,
		const PfoToLArTPCMap &	pfoToLArTPCMap,
		ThreeDPointingClusterMap &	pointingClusterMap
	)		const

private

Build a 3D pointing cluster for each Pfo.

Parameters

inputPfoList	the input list of Pfos
pfoToLArTPCMap	the input mapping between Pfos and tpc
pointingClusterMap	the mapping between Pfos and their corresponding 3D pointing clusters

Definition at line 98 of file StitchingCosmicRayMergingTool.cc.

 {
     for (const ParticleFlowObject *const pPfo : inputPfoList)
     {
         try
         {
             PfoToLArTPCMap::const_iterator tpcIter(pfoToLArTPCMap.find(pPfo));
 
             if (pfoToLArTPCMap.end() == tpcIter)
                 throw StatusCodeException(STATUS_CODE_NOT_FOUND);
 
             const float slidingFitPitch(tpcIter->second->GetWirePitchW());
 
             ClusterList clusterList;
             LArPfoHelper::GetThreeDClusterList(pPfo, clusterList);
 
             if (1 != clusterList.size())
                 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
 
             const ThreeDSlidingFitResult slidingFitResult(clusterList.front(), m_halfWindowLayers, slidingFitPitch);
             (void)pointingClusterMap.insert(ThreeDPointingClusterMap::value_type(pPfo, LArPointingCluster(slidingFitResult)));
         }
         catch (const StatusCodeException &)
         {
         }
     }
 }

void lar_content::StitchingCosmicRayMergingTool::BuildTPCMaps	(	const pandora::PfoList &	inputPfoList,
		const PfoToLArTPCMap &	pfoToLArTPCMap,
		LArTPCToPfoMap &	larTPCToPfoMap
	)		const

private

Build a list of Pfos for each tpc.

Parameters

inputPfoList	the input list of Pfos
pfoToLArTPCMap	the input mapping between Pfos and tpc
larTPCToPfoMap	the output mapping between tpc and Pfos

Definition at line 129 of file StitchingCosmicRayMergingTool.cc.

 {
     for (const ParticleFlowObject *const pPfo : inputPfoList)
     {
         PfoToLArTPCMap::const_iterator iter(pfoToLArTPCMap.find(pPfo));
 
         if (pfoToLArTPCMap.end() != iter)
             larTPCToPfoMap[iter->second].push_back(pPfo);
     }
 }

bool lar_content::StitchingCosmicRayMergingTool::CalculateX0	(	const PfoToLArTPCMap &	pfoToLArTPCMap,
		const ThreeDPointingClusterMap &	pointingClusterMap,
		const pandora::PfoVector &	pfoVector,
		float &	x0,
		PfoToPointingVertexMatrix &	pfoToPointingVertexMatrix
	)		const

private

Calculate x0 shift for a group of associated Pfos.

Parameters

pfoToLArTPCMap	the mapping between pfos and tpc
pointingClusterMap	the mapping between Pfos and their corresponding 3D pointing clusters
pfoVector	the vector of parent Pfos to stitch together
x0	the output x0 value
pfoToPointingVertexMatrix	map of pfo to a map of matched pfo and the corresponding pointing vertex used in stitching

Returns: bool true if all x0 contributions were consistent, false otherwise

Definition at line 662 of file StitchingCosmicRayMergingTool.cc.

 {
     float sumX(0.f), sumN(0.f);
 
     for (PfoVector::const_iterator iter1 = pfoVector.begin(), iterEnd = pfoVector.end(); iter1 != iterEnd; ++iter1)
     {
         const ParticleFlowObject *const pPfo1(*iter1);
         PfoToLArTPCMap::const_iterator tpcIter1 = pfoToLArTPCMap.find(pPfo1);
         ThreeDPointingClusterMap::const_iterator pointingIter1 = pointingClusterMap.find(pPfo1);
 
         if (pfoToLArTPCMap.end() == tpcIter1 || pointingClusterMap.end() == pointingIter1)
             throw StatusCodeException(STATUS_CODE_FAILURE);
 
         const LArTPC *const pLArTPC1(tpcIter1->second);
         const LArPointingCluster &pointingCluster1(pointingIter1->second);
 
         for (PfoVector::const_iterator iter2 = std::next(iter1); iter2 != iterEnd; ++iter2)
         {
             const ParticleFlowObject *const pPfo2(*iter2);
             PfoToLArTPCMap::const_iterator tpcIter2 = pfoToLArTPCMap.find(pPfo2);
             ThreeDPointingClusterMap::const_iterator pointingIter2 = pointingClusterMap.find(pPfo2);
 
             if (pfoToLArTPCMap.end() == tpcIter2 || pointingClusterMap.end() == pointingIter2)
                 throw StatusCodeException(STATUS_CODE_FAILURE);
 
             const LArTPC *const pLArTPC2(tpcIter2->second);
             const LArPointingCluster &pointingCluster2(pointingIter2->second);
 
             if (!LArStitchingHelper::CanTPCsBeStitched(*pLArTPC1, *pLArTPC2))
                 continue;
 
             // Calculate X0 for the closest pair of vertices
             LArPointingCluster::Vertex pointingVertex1, pointingVertex2;
             try
             {
                 LArStitchingHelper::GetClosestVertices(*pLArTPC1, *pLArTPC2, pointingCluster1, pointingCluster2, pointingVertex1, pointingVertex2);
 
                 // Record pfo1 stitching vertex for pfo1<->pfo2 match, used to determine shifting direction in later step
                 const PfoToPointingVertexMatrix::iterator pfoToPointingVertexMatrixIter1(pfoToPointingVertexMatrix.find(pPfo1));
                 if (pfoToPointingVertexMatrixIter1 == pfoToPointingVertexMatrix.end())
                 {
                     // No matches present in map, add pfo1<->pfo2 match
                     const PfoToPointingVertexMap pfoToPointingVertexMap({{pPfo2, pointingVertex1}});
                     (void)pfoToPointingVertexMatrix.insert(PfoToPointingVertexMatrix::value_type(pPfo1, pfoToPointingVertexMap));
                 }
                 else
                 {
                     // ATTN: another match for a different TPC boundary may be present, add pfo1<->pfo2 match
                     PfoToPointingVertexMap &pfoToPointingVertexMap(pfoToPointingVertexMatrixIter1->second);
                     const PfoToPointingVertexMap::iterator pfoToPointingVertexMapIter(pfoToPointingVertexMap.find(pPfo2));
                     if (pfoToPointingVertexMapIter == pfoToPointingVertexMap.end())
                     {
                         (void)pfoToPointingVertexMap.insert(PfoToPointingVertexMap::value_type(pPfo2, pointingVertex1));
                     }
                     else
                     {
                         if ((pfoToPointingVertexMapIter->second.GetPosition() - pointingVertex1.GetPosition()).GetMagnitude() >
                             std::numeric_limits<float>::epsilon())
                             throw StatusCodeException(STATUS_CODE_FAILURE);
                     }
                 }
 
                 // Record pfo2 stitching vertex for pfo1<->pfo2 match, used to determine shifting direction in later step
                 const PfoToPointingVertexMatrix::iterator pfoToPointingVertexMatrixIter2(pfoToPointingVertexMatrix.find(pPfo2));
                 if (pfoToPointingVertexMatrixIter2 == pfoToPointingVertexMatrix.end())
                 {
                     // No matches present in map, add pfo1<->pfo2 match
                     const PfoToPointingVertexMap pfoToPointingVertexMap({{pPfo1, pointingVertex2}});
                     (void)pfoToPointingVertexMatrix.insert(PfoToPointingVertexMatrix::value_type(pPfo2, pfoToPointingVertexMap));
                 }
                 else
                 {
                     // ATTN: another match for a different TPC boundary may be present, add pfo1<->pfo2 match
                     PfoToPointingVertexMap &pfoToPointingVertexMap(pfoToPointingVertexMatrixIter2->second);
                     const PfoToPointingVertexMap::iterator pfoToPointingVertexMapIter(pfoToPointingVertexMap.find(pPfo1));
                     if (pfoToPointingVertexMapIter == pfoToPointingVertexMap.end())
                     {
                         (void)pfoToPointingVertexMap.insert(PfoToPointingVertexMap::value_type(pPfo1, pointingVertex2));
                     }
                     else
                     {
                         if ((pfoToPointingVertexMapIter->second.GetPosition() - pointingVertex2.GetPosition()).GetMagnitude() >
                             std::numeric_limits<float>::epsilon())
                             throw StatusCodeException(STATUS_CODE_FAILURE);
                     }
                 }
 
                 const float tpcBoundaryCenterX(LArStitchingHelper::GetTPCBoundaryCenterX(*pLArTPC1, *pLArTPC2));
                 const bool isCPAStitch(pLArTPC1->GetCenterX() < tpcBoundaryCenterX ? !pLArTPC1->IsDriftInPositiveX() : !pLArTPC2->IsDriftInPositiveX());
                 float thisX0(LArStitchingHelper::CalculateX0(*pLArTPC1, *pLArTPC2, pointingVertex1, pointingVertex2));
 
                 thisX0 *= isCPAStitch ? -1.f : 1.f;
 
                 // If multiple boundaries identified, check if stitching contribution is consistent
                 if ((sumN > std::numeric_limits<float>::epsilon()) && (sumX > std::numeric_limits<float>::epsilon()))
                 {
                     const float fractionalDiff(std::fabs((sumX - (thisX0 * sumN)) / sumX));
 
                     if ((fractionalDiff > m_maxX0FractionalDeviation) && (std::fabs(sumX / sumN) > m_boundaryToleranceWidth))
                         return false;
                 }
 
                 sumX += thisX0;
                 sumN += 1.f;
             }
             catch (const pandora::StatusCodeException &statusCodeException)
             {
                 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
                     std::cout << "StitchingCosmicRayMergingTool: Attempting to stitch a pfo with multiple vertices for the same match" << std::endl;
             }
         }
     }
 
     if ((sumN < std::numeric_limits<float>::epsilon()) || (std::fabs(sumX) < std::numeric_limits<float>::epsilon()))
         throw StatusCodeException(STATUS_CODE_NOT_FOUND);
 
     x0 = (sumX / sumN);
 
     return true;
 }

void lar_content::StitchingCosmicRayMergingTool::CollectAssociatedPfos	(	const pandora::ParticleFlowObject *const	pSeedPfo,
		const pandora::ParticleFlowObject *const	pCurrentPfo,
		const PfoMergeMap &	pfoMerges,
		const pandora::PfoSet &	vetoSet,
		pandora::PfoList &	associatedList
	)		const

private

Collect up associations between Pfos.

Parameters

pSeedPfo	the seed Pfo (will be enlarged after merging)
pCurrentPfo	the target Pfo (will be deleted after merging)
pfoMerges	the list of associations between Pfos
vetoSet	the set of associated Pfos that have already been considered
associatedList	the output list of associated Pfos

Definition at line 432 of file StitchingCosmicRayMergingTool.cc.

 {
     if (vetoSet.count(pCurrentPfo))
         return;
 
     PfoMergeMap::const_iterator iter1 = pfoMergeMap.find(pCurrentPfo);
 
     if (pfoMergeMap.end() == iter1)
         return;
 
     for (PfoList::const_iterator iter2 = iter1->second.begin(), iterEnd2 = iter1->second.end(); iter2 != iterEnd2; ++iter2)
     {
         const ParticleFlowObject *const pAssociatedPfo = *iter2;
 
         if (pAssociatedPfo == pSeedPfo)
             continue;
 
         if (associatedList.end() != std::find(associatedList.begin(), associatedList.end(), pAssociatedPfo))
             continue;
 
         associatedList.push_back(pAssociatedPfo);
 
         this->CollectAssociatedPfos(pSeedPfo, pAssociatedPfo, pfoMergeMap, vetoSet, associatedList);
     }
 }

void lar_content::StitchingCosmicRayMergingTool::CreatePfoMatches	(	const LArTPCToPfoMap &	larTPCToPfoMap,
		const ThreeDPointingClusterMap &	pointingClusterMap,
		PfoAssociationMatrix &	pfoAssociationMatrix
	)		const

private

Create associations between Pfos using 3D pointing clusters.

Parameters

larTPCToPfoMap	the input mapping between tpc and Pfos
pointingClusterMap	the input mapping between Pfos and their corresponding 3D pointing clusters
pfoAssociationMatrix	the output matrix of associations between Pfos

Definition at line 142 of file StitchingCosmicRayMergingTool.cc.

 {
     LArTPCVector larTPCVector;
     for (const auto &mapEntry : larTPCToPfoMap)
         larTPCVector.push_back(mapEntry.first);
     std::sort(larTPCVector.begin(), larTPCVector.end(), LArStitchingHelper::SortTPCs);
 
     for (LArTPCVector::const_iterator tpcIter1 = larTPCVector.begin(), tpcIterEnd = larTPCVector.end(); tpcIter1 != tpcIterEnd; ++tpcIter1)
     {
         const LArTPC *const pLArTPC1(*tpcIter1);
         const PfoList &pfoList1(larTPCToPfoMap.at(pLArTPC1));
 
         for (LArTPCVector::const_iterator tpcIter2 = tpcIter1; tpcIter2 != tpcIterEnd; ++tpcIter2)
         {
             const LArTPC *const pLArTPC2(*tpcIter2);
             const PfoList &pfoList2(larTPCToPfoMap.at(pLArTPC2));
 
             if (!LArStitchingHelper::CanTPCsBeStitched(*pLArTPC1, *pLArTPC2))
                 continue;
 
             for (const ParticleFlowObject *const pPfo1 : pfoList1)
             {
                 for (const ParticleFlowObject *const pPfo2 : pfoList2)
                     this->CreatePfoMatches(*pLArTPC1, *pLArTPC2, pPfo1, pPfo2, pointingClusterMap, pfoAssociationMatrix);
             }
         }
     }
 }

void lar_content::StitchingCosmicRayMergingTool::CreatePfoMatches	(	const pandora::LArTPC &	larTPC1,
		const pandora::LArTPC &	larTPC2,
		const pandora::ParticleFlowObject *const	pPfo1,
		const pandora::ParticleFlowObject *const	pPfo2,
		const ThreeDPointingClusterMap &	pointingClusterMap,
		PfoAssociationMatrix &	pfoAssociationMatrix
	)		const

private

Create associations between Pfos using 3D pointing clusters.

Parameters

larTPC1	the tpc description for the first Pfo
larTPC2	the tpc description for the second Pfo
pPfo1	the first Pfo
pPfo2	the second Pfo
pointingClusterMap	the input mapping between Pfos and their corresponding 3D pointing clusters
pfoAssociationMatrix	the output matrix of associations between Pfos

void lar_content::StitchingCosmicRayMergingTool::OrderPfoMerges	(	const PfoToLArTPCMap &	pfoToLArTPCMap,
		const ThreeDPointingClusterMap &	pointingClusterMap,
		const PfoMergeMap &	inputPfoMerges,
		PfoMergeMap &	outputPfoMerges
	)		const

private

Identify the vertex Pfo and then re-order the map of merges so that the vertex Pfo will be enlarged.

Parameters

pfoToLArTPCMap	the mapping between pfos and tpc
pointingClusterMap	the mapping between Pfos and their corresponding 3D pointing clusters
inputPfoMerges	the input map of Pfo merges
outputPfoMerges	the re-ordered map of Pfo merges

Definition at line 461 of file StitchingCosmicRayMergingTool.cc.

 {
     PfoVector inputPfoVector;
     for (const auto &mapEntry : inputPfoMerges)
         inputPfoVector.push_back(mapEntry.first);
     std::sort(inputPfoVector.begin(), inputPfoVector.end(), LArPfoHelper::SortByNHits);
 
     for (const ParticleFlowObject *const pInputPfo : inputPfoVector)
     {
         const PfoList &pfoList(inputPfoMerges.at(pInputPfo));
 
         float bestLength(0.f);
         const ParticleFlowObject *pVertexPfo = nullptr;
 
         for (PfoList::const_iterator iter1 = pfoList.begin(), iterEnd = pfoList.end(); iter1 != iterEnd; ++iter1)
         {
             const ParticleFlowObject *const pPfo1(*iter1);
             PfoToLArTPCMap::const_iterator tpcIter1 = pfoToLArTPCMap.find(pPfo1);
             ThreeDPointingClusterMap::const_iterator pointingIter1 = pointingClusterMap.find(pPfo1);
 
             if (pfoToLArTPCMap.end() == tpcIter1 || pointingClusterMap.end() == pointingIter1)
                 throw StatusCodeException(STATUS_CODE_FAILURE);
 
             const LArTPC *const pLArTPC1(tpcIter1->second);
             const LArPointingCluster &pointingCluster1(pointingIter1->second);
 
             for (PfoList::const_iterator iter2 = iter1; iter2 != iterEnd; ++iter2)
             {
                 const ParticleFlowObject *const pPfo2(*iter2);
                 PfoToLArTPCMap::const_iterator tpcIter2 = pfoToLArTPCMap.find(pPfo2);
                 ThreeDPointingClusterMap::const_iterator pointingIter2 = pointingClusterMap.find(pPfo2);
 
                 if (pfoToLArTPCMap.end() == tpcIter2 || pointingClusterMap.end() == pointingIter2)
                     throw StatusCodeException(STATUS_CODE_FAILURE);
 
                 const LArTPC *const pLArTPC2(tpcIter2->second);
                 const LArPointingCluster &pointingCluster2(pointingIter2->second);
 
                 if (pLArTPC1 == pLArTPC2)
                     continue;
 
                 const float thisLength(LArStitchingHelper::GetTPCDisplacement(*pLArTPC1, *pLArTPC2));
 
                 if (thisLength < bestLength)
                     continue;
 
                 bestLength = thisLength;
 
                 try
                 {
                     pVertexPfo = nullptr;
 
                     LArPointingCluster::Vertex nearVertex1, nearVertex2;
                     LArStitchingHelper::GetClosestVertices(*pLArTPC1, *pLArTPC2, pointingCluster1, pointingCluster2, nearVertex1, nearVertex2);
 
                     const LArPointingCluster::Vertex &farVertex1(
                         nearVertex1.IsInnerVertex() ? pointingCluster1.GetOuterVertex() : pointingCluster1.GetInnerVertex());
                     const LArPointingCluster::Vertex &farVertex2(
                         nearVertex2.IsInnerVertex() ? pointingCluster2.GetOuterVertex() : pointingCluster2.GetInnerVertex());
                     const float deltaY(farVertex1.GetPosition().GetY() - farVertex2.GetPosition().GetY());
 
                     if (std::fabs(deltaY) < std::numeric_limits<float>::epsilon())
                         throw StatusCodeException(STATUS_CODE_NOT_FOUND);
 
                     pVertexPfo = ((deltaY > 0.f) ? pPfo1 : pPfo2);
                 }
                 catch (const pandora::StatusCodeException &)
                 {
                 }
             }
         }
 
         if (pVertexPfo)
             outputPfoMerges[pVertexPfo].insert(outputPfoMerges[pVertexPfo].begin(), pfoList.begin(), pfoList.end());
     }
 }

StatusCode lar_content::StitchingCosmicRayMergingTool::ReadSettings ( const pandora::TiXmlHandle xmlHandle )

private

Definition at line 818 of file StitchingCosmicRayMergingTool.cc.

 {
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ThreeDStitchingMode", m_useXcoordinate));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "AlwaysApplyT0Calculation", m_alwaysApplyT0Calculation));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HalfWindowLayers", m_halfWindowLayers));
 
     float minLength(std::sqrt(m_minLengthSquared));
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinPfoLength", minLength));
     m_minLengthSquared = minLength * minLength;
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosRelativeAngle", m_minCosRelativeAngle));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "RelaxMinLongitudinalDisplacement", m_relaxMinLongitudinalDisplacement));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacementX", m_maxLongitudinalDisplacementX));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "LooserMinCosRelativeAngle", m_relaxCosRelativeAngle));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "LooserMaxTransverseDisplacement", m_relaxTransverseDisplacement));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinNCaloHits3D", m_minNCaloHits3D));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
         XmlHelper::ReadValue(xmlHandle, "MaxX0FractionalDeviation", m_maxX0FractionalDeviation));
 
     PANDORA_RETURN_RESULT_IF_AND_IF(
         STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "BoundaryToleranceWidth", m_boundaryToleranceWidth));
 
     return STATUS_CODE_SUCCESS;
 }

void lar_content::StitchingCosmicRayMergingTool::Run	(	const MasterAlgorithm *const	pAlgorithm,
		const pandora::PfoList *const	pMultiPfoList,
		PfoToLArTPCMap &	pfoToLArTPCMap,
		PfoToFloatMap &	stitchedPfosToX0Map
	)

virtual

Run the algorithm tool.

Parameters

pAlgorithm	address of the calling algorithm
pMultiPfoList	the list of pfos in multiple lar tpcs
pfoToLArTPCMap	the pfo to lar tpc map
stitchedPfosToX0Map	a map of cosmic-ray pfos that have been stitched between lar tpcs to the X0 shift

Implements lar_content::StitchingBaseTool.

Definition at line 42 of file StitchingCosmicRayMergingTool.cc.

 {
     if (PandoraContentApi::GetSettings(*pAlgorithm)->ShouldDisplayAlgorithmInfo())
         std::cout << "----> Running Algorithm Tool: " << this->GetInstanceName() << ", " << this->GetType() << std::endl;
 
     if (this->GetPandora().GetGeometry()->GetLArTPCMap().size() < 2)
         return;
 
     if (pfoToLArTPCMap.empty())
         throw StatusCodeException(STATUS_CODE_NOT_FOUND);
 
     PfoList primaryPfos;
     this->SelectPrimaryPfos(pMultiPfoList, pfoToLArTPCMap, primaryPfos);
 
     ThreeDPointingClusterMap pointingClusterMap;
     this->BuildPointingClusterMaps(primaryPfos, pfoToLArTPCMap, pointingClusterMap);
 
     LArTPCToPfoMap larTPCToPfoMap;
     this->BuildTPCMaps(primaryPfos, pfoToLArTPCMap, larTPCToPfoMap);
 
     PfoAssociationMatrix pfoAssociationMatrix;
     this->CreatePfoMatches(larTPCToPfoMap, pointingClusterMap, pfoAssociationMatrix);
 
     PfoMergeMap pfoSelectedMatches;
     this->SelectPfoMatches(pfoAssociationMatrix, pfoSelectedMatches);
 
     PfoMergeMap pfoSelectedMerges;
     this->SelectPfoMerges(pfoSelectedMatches, pfoSelectedMerges);
 
     PfoMergeMap pfoOrderedMerges;
     this->OrderPfoMerges(pfoToLArTPCMap, pointingClusterMap, pfoSelectedMerges, pfoOrderedMerges);
 
     this->StitchPfos(pAlgorithm, pointingClusterMap, pfoOrderedMerges, pfoToLArTPCMap, stitchedPfosToX0Map);
 }

void lar_content::StitchingCosmicRayMergingTool::SelectPfoMatches	(	const PfoAssociationMatrix &	pfoAssociationMatrix,
		PfoMergeMap &	pfoSelectedMatches
	)		const

private

Select the best associations between Pfos; create a mapping between associated Pfos, handling any ambiguities.

Parameters

pfoAssociationMatrix	the input list of all associations between Pfos
pfoSelectedMatches	the output list of good associations between pfos (candidates for merging)

Definition at line 291 of file StitchingCosmicRayMergingTool.cc.

 {
     // First step: loop over association matrix and find best associations A -> X and B -> Y
     // =====================================================================================
     PfoAssociationMatrix bestAssociationMatrix;
 
     PfoVector pfoVector1;
     for (const auto &mapEntry : pfoAssociationMatrix)
         pfoVector1.push_back(mapEntry.first);
     std::sort(pfoVector1.begin(), pfoVector1.end(), LArPfoHelper::SortByNHits);
 
     for (const ParticleFlowObject *const pPfo1 : pfoVector1)
     {
         const PfoAssociationMap &pfoAssociationMap(pfoAssociationMatrix.at(pPfo1));
 
         const ParticleFlowObject *pBestPfoInner(nullptr);
         PfoAssociation bestAssociationInner(PfoAssociation::UNDEFINED, PfoAssociation::UNDEFINED, 0.f);
 
         const ParticleFlowObject *pBestPfoOuter(nullptr);
         PfoAssociation bestAssociationOuter(PfoAssociation::UNDEFINED, PfoAssociation::UNDEFINED, 0.f);
 
         PfoVector pfoVector2;
         for (const auto &mapEntry : pfoAssociationMap)
             pfoVector2.push_back(mapEntry.first);
         std::sort(pfoVector2.begin(), pfoVector2.end(), LArPfoHelper::SortByNHits);
 
         for (const ParticleFlowObject *const pPfo2 : pfoVector2)
         {
             const PfoAssociation &pfoAssociation(pfoAssociationMap.at(pPfo2));
 
             // Inner associations
             if (pfoAssociation.GetParent() == PfoAssociation::INNER)
             {
                 if (pfoAssociation.GetFigureOfMerit() > bestAssociationInner.GetFigureOfMerit())
                 {
                     bestAssociationInner = pfoAssociation;
                     pBestPfoInner = pPfo2;
                 }
             }
 
             // Outer associations
             if (pfoAssociation.GetParent() == PfoAssociation::OUTER)
             {
                 if (pfoAssociation.GetFigureOfMerit() > bestAssociationOuter.GetFigureOfMerit())
                 {
                     bestAssociationOuter = pfoAssociation;
                     pBestPfoOuter = pPfo2;
                 }
             }
         }
 
         if (pBestPfoInner)
             (void)bestAssociationMatrix[pPfo1].insert(PfoAssociationMap::value_type(pBestPfoInner, bestAssociationInner));
 
         if (pBestPfoOuter)
             (void)bestAssociationMatrix[pPfo1].insert(PfoAssociationMap::value_type(pBestPfoOuter, bestAssociationOuter));
     }
 
     // Second step: make the merge if A -> X and B -> Y is in fact A -> B and B -> A
     // =============================================================================
     PfoVector pfoVector3;
     for (const auto &mapEntry : bestAssociationMatrix)
         pfoVector3.push_back(mapEntry.first);
     std::sort(pfoVector3.begin(), pfoVector3.end(), LArPfoHelper::SortByNHits);
 
     for (const ParticleFlowObject *const pParentPfo : pfoVector3)
     {
         const PfoAssociationMap &parentAssociationMap(bestAssociationMatrix.at(pParentPfo));
 
         PfoVector pfoVector4;
         for (const auto &mapEntry : parentAssociationMap)
             pfoVector4.push_back(mapEntry.first);
         std::sort(pfoVector4.begin(), pfoVector4.end(), LArPfoHelper::SortByNHits);
 
         for (const ParticleFlowObject *const pDaughterPfo : pfoVector4)
         {
             const PfoAssociation &parentToDaughterAssociation(parentAssociationMap.at(pDaughterPfo));
             PfoAssociationMatrix::const_iterator iter5 = bestAssociationMatrix.find(pDaughterPfo);
 
             if (bestAssociationMatrix.end() == iter5)
                 continue;
 
             const PfoAssociationMap &daughterAssociationMap(iter5->second);
 
             PfoAssociationMap::const_iterator iter6 = daughterAssociationMap.find(pParentPfo);
             if (daughterAssociationMap.end() == iter6)
                 continue;
 
             const PfoAssociation &daughterToParentAssociation(iter6->second);
 
             if (parentToDaughterAssociation.GetParent() == daughterToParentAssociation.GetDaughter() &&
                 parentToDaughterAssociation.GetDaughter() == daughterToParentAssociation.GetParent())
             {
                 pfoMatches[pParentPfo].push_back(pDaughterPfo);
             }
         }
     }
 }

void lar_content::StitchingCosmicRayMergingTool::SelectPfoMerges	(	const PfoMergeMap &	pfoMatches,
		PfoMergeMap &	pfoMerges
	)		const

private

Create an initial map of Pfo merges to be made.

Parameters

pfoMatches	the input list of associated Pfos
pfoMerges	the output mapping between Pfos to enlarge and Pfos to delete

Definition at line 392 of file StitchingCosmicRayMergingTool.cc.

 {
     PfoSet vetoSet;
 
     PfoVector inputPfoVector;
     for (const auto &mapEntry : pfoMatches)
         inputPfoVector.push_back(mapEntry.first);
     std::sort(inputPfoVector.begin(), inputPfoVector.end(), LArPfoHelper::SortByNHits);
 
     for (const ParticleFlowObject *const pInputPfo : inputPfoVector)
     {
         const PfoList &pfoList(pfoMatches.at(pInputPfo));
 
         for (const ParticleFlowObject *const pSeedPfo : pfoList)
         {
             if (vetoSet.count(pSeedPfo))
                 continue;
 
             PfoList mergeList;
             this->CollectAssociatedPfos(pSeedPfo, pSeedPfo, pfoMatches, vetoSet, mergeList);
 
             vetoSet.insert(pSeedPfo);
             PfoList &selectedPfoList(pfoMerges[pSeedPfo]);
             selectedPfoList.push_back(pSeedPfo);
 
             for (const ParticleFlowObject *const pAssociatedPfo : mergeList)
             {
                 // Check if particle has already been counted
                 if (vetoSet.count(pAssociatedPfo) || (selectedPfoList.end() != std::find(selectedPfoList.begin(), selectedPfoList.end(), pAssociatedPfo)))
                     throw StatusCodeException(STATUS_CODE_FAILURE);
 
                 vetoSet.insert(pAssociatedPfo);
                 selectedPfoList.push_back(pAssociatedPfo);
             }
         }
     }
 }

void lar_content::StitchingCosmicRayMergingTool::SelectPrimaryPfos	(	const pandora::PfoList *	pInputPfoList,
		const PfoToLArTPCMap &	pfoToLArTPCMap,
		pandora::PfoList &	outputPfoList
	)		const

private

Select primary Pfos from the input list of Pfos.

Parameters

pInputPfoList	the input list of Pfos
pfoToLArTPCMap	the input mapping between Pfos and tpc
outputPfoList	the output list of Pfos

Definition at line 80 of file StitchingCosmicRayMergingTool.cc.

 {
     for (const ParticleFlowObject *const pPfo : *pInputPfoList)
     {
         if (!LArPfoHelper::IsFinalState(pPfo) || !LArPfoHelper::IsTrack(pPfo))
             continue;
 
         if (!pfoToLArTPCMap.count(pPfo))
             continue;
 
         outputPfoList.push_back(pPfo);
     }
 
     outputPfoList.sort(LArPfoHelper::SortByNHits);
 }

void lar_content::StitchingCosmicRayMergingTool::ShiftPfo	(	const MasterAlgorithm *const	pAlgorithm,
		const pandora::ParticleFlowObject *const	pPfoToShift,
		const pandora::ParticleFlowObject *const	pMatchedPfo,
		const float	x0,
		const PfoToLArTPCMap &	pfoToLArTPCMap,
		const PfoToPointingVertexMatrix &	pfoToPointingVertexMatrix
	)		const

private

Shift a pfo given its pfo stitching pair.

Parameters

pPfoToShift	the pfo of the stitching pair to shift
pMatchedPfo	the pfo of the stitching pair to remain stationary
x0	the distance by which pPfoToShift is to be shifted (direction of shift is determined in method)
pfoToLArTPCMap	the pfo to lar tpc map
pfoToPointingVertexMatrix	the map [pfo -> map [matched pfo -> pfo stitching vertex]]

Definition at line 618 of file StitchingCosmicRayMergingTool.cc.

 {
     // get stitching vertex for the pfo to be shifted
     const PfoToPointingVertexMatrix::const_iterator pfoToPointingVertexMatrixIter(pfoToPointingVertexMatrix.find(pPfoToShift));
     const LArPointingCluster::Vertex stitchingVertex(pfoToPointingVertexMatrixIter->second.at(pMatchedPfo));
 
     const LArTPC *const pShiftLArTPC(pfoToLArTPCMap.at(pPfoToShift));
     const LArTPC *const pMatchedLArTPC(pfoToLArTPCMap.at(pMatchedPfo));
 
     // determine shift sign from the relative position of stitching vertex and the relevant TPC boundary position
     const float tpcBoundaryCenterX(LArStitchingHelper::GetTPCBoundaryCenterX(*pShiftLArTPC, *pMatchedLArTPC));
     float tpcBoundaryX(0.f);
 
     if (pShiftLArTPC->GetCenterX() < tpcBoundaryCenterX)
     {
         tpcBoundaryX = pShiftLArTPC->GetCenterX() + (pShiftLArTPC->GetWidthX() / 2.f);
     }
     else
     {
         tpcBoundaryX = pShiftLArTPC->GetCenterX() - (pShiftLArTPC->GetWidthX() / 2.f);
     }
 
     const float positionShiftSign = stitchingVertex.GetPosition().GetX() < tpcBoundaryX ? 1.f : -1.f;
 
     // ATTN: No CPA/APA sign needed since x0 calculation corresponds to an APA
     object_creation::ParticleFlowObject::Metadata metadata;
     metadata.m_propertiesToAdd["X0"] = x0;
 
     // ATTN: Set the X0 shift for all particles in hierarchy
     PfoList downstreamPfoList;
     LArPfoHelper::GetAllDownstreamPfos(pPfoToShift, downstreamPfoList);
 
     for (const ParticleFlowObject *const pHierarchyPfo : downstreamPfoList)
         PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*pAlgorithm, pHierarchyPfo, metadata));
 
     const float signedX0(std::fabs(x0) * positionShiftSign);
 
     pAlgorithm->ShiftPfoHierarchy(pPfoToShift, pfoToLArTPCMap, signedX0);
 }

void lar_content::StitchingCosmicRayMergingTool::StitchPfos	(	const MasterAlgorithm *const	pAlgorithm,
		const ThreeDPointingClusterMap &	pointingClusterMap,
		const PfoMergeMap &	pfoMerges,
		PfoToLArTPCMap &	pfoToLArTPCMap,
		PfoToFloatMap &	stitchedPfosToX0Map
	)		const

private

Apply X0 corrections, and then stitch together Pfos.

Parameters

pAlgorithm	the address of the parent stitching algorithm
pointingClusterMap	the mapping between Pfos and their corresponding 3D pointing clusters
pfoMerges	the input map of Pfo merges
pfoToLArTPCMap	the pfo to lar tpc map
stitchedPfosToX0Map	a map of cosmic-ray pfos that have been stitched between lar tpcs to the X0 shift

Definition at line 541 of file StitchingCosmicRayMergingTool.cc.

 {
     PfoVector pfoVectorToEnlarge;
     for (const auto &mapEntry : pfoMerges)
         pfoVectorToEnlarge.push_back(mapEntry.first);
     std::sort(pfoVectorToEnlarge.begin(), pfoVectorToEnlarge.end(), LArPfoHelper::SortByNHits);
 
     for (const ParticleFlowObject *const pPfoToEnlarge : pfoVectorToEnlarge)
     {
         const PfoList &pfoList(pfoMerges.at(pPfoToEnlarge));
         const PfoVector pfoVector(pfoList.begin(), pfoList.end());
         PfoToPointingVertexMatrix pfoToPointingVertexMatrix;
         float x0(0.f);
 
         if (!m_useXcoordinate || m_alwaysApplyT0Calculation)
         {
             try
             {
                 // If stitching contributions are inconsistent, abort
                 if (!this->CalculateX0(pfoToLArTPCMap, pointingClusterMap, pfoVector, x0, pfoToPointingVertexMatrix))
                     continue;
             }
             catch (const pandora::StatusCodeException &)
             {
                 continue;
             }
         }
 
         // ATTN: shift the pfos one at a time
         PfoSet shiftedPfos;
         for (PfoVector::const_iterator iterI = pfoVector.begin(); iterI != pfoVector.end(); ++iterI)
         {
             const ParticleFlowObject *const pPfoI(*iterI);
             const LArTPC *const pLArTPCI(pfoToLArTPCMap.at(pPfoI));
 
             for (PfoVector::const_iterator iterJ = std::next(iterI); iterJ != pfoVector.end(); ++iterJ)
             {
                 const ParticleFlowObject *const pPfoJ(*iterJ);
                 const LArTPC *const pLArTPCJ(pfoToLArTPCMap.at(pPfoJ));
 
                 if (!LArStitchingHelper::CanTPCsBeStitched(*pLArTPCI, *pLArTPCJ))
                     continue;
 
                 if (std::find(shiftedPfos.begin(), shiftedPfos.end(), pPfoI) == shiftedPfos.end())
                 {
                     if (!m_useXcoordinate || m_alwaysApplyT0Calculation)
                         this->ShiftPfo(pAlgorithm, pPfoI, pPfoJ, x0, pfoToLArTPCMap, pfoToPointingVertexMatrix);
 
                     shiftedPfos.insert(pPfoI);
                 }
 
                 if (std::find(shiftedPfos.begin(), shiftedPfos.end(), pPfoJ) == shiftedPfos.end())
                 {
                     if (!m_useXcoordinate || m_alwaysApplyT0Calculation)
                         this->ShiftPfo(pAlgorithm, pPfoJ, pPfoI, x0, pfoToLArTPCMap, pfoToPointingVertexMatrix);
 
                     shiftedPfos.insert(pPfoJ);
                 }
             }
         }
 
         // now merge all pfos
         for (const ParticleFlowObject *const pPfoToDelete : shiftedPfos)
         {
             if (pPfoToDelete == pPfoToEnlarge)
                 continue;
 
             pAlgorithm->StitchPfos(pPfoToEnlarge, pPfoToDelete, pfoToLArTPCMap);
         }
 
         stitchedPfosToX0Map.insert(PfoToFloatMap::value_type(pPfoToEnlarge, x0));
     }
 }