lar_cluster3d::ClusterMergeAlg Class Reference

Inheritance diagram for lar_cluster3d::ClusterMergeAlg:

Public Member Functions
	ClusterMergeAlg (const fhicl::ParameterSet &)
	Constructor. More...
	~ClusterMergeAlg ()
	Destructor. More...
void	configure (fhicl::ParameterSet const &pset) override
void	initializeHistograms (art::TFileDirectory &) override
	Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder. More...
void	ModifyClusters (reco::ClusterParametersList &) const override
	Scan an input collection of clusters and modify those according to the specific implementing algorithm. More...
float	getTimeToExecute () const override
	If monitoring, recover the time to execute a particular function. More...
Public Member Functions inherited from lar_cluster3d::IClusterModAlg
virtual	~IClusterModAlg () noexcept=default
	Virtual Destructor. More...
virtual void	configure (const fhicl::ParameterSet &)=0
	Interface for configuring the particular algorithm tool. More...

Private Member Functions
bool	linearClusters (reco::ClusterParameters &, reco::ClusterParameters &) const
bool	mergeClusters (reco::ClusterParameters &, reco::ClusterParameters &) const
float	closestApproach (const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, Eigen::Vector3f &, Eigen::Vector3f &, Eigen::Vector3f &) const
const reco::ClusterHit3D *	findClosestHit3D (const Eigen::Vector3f &, const Eigen::Vector3f &, const reco::HitPairListPtr &) const
const reco::ClusterHit3D *	findFurthestHit3D (const Eigen::Vector3f &, const Eigen::Vector3f &, const reco::HitPairListPtr &) const

Private Attributes
bool	fEnableMonitoring
	Data members to follow. More...
float	fMinTransEigenVal
	Set a mininum allowed value for the transverse eigen values. More...
float	fMinEigenToProcess
	Don't look anymore at clusters below this size. More...
bool	fOutputHistograms
	Take the time to create and fill some histograms for diagnostics. More...
float	fTimeToProcess
	Keep track of how long it took to run this algorithm. More...
std::vector< TH1F * >	fFirstEigenValueHists
	First Cluster eigen value. More...
std::vector< TH1F * >	fNextEigenValueHists
	Next Cluster eigen value. More...
TH1F *	fNumMergedClusters
	How many clusters were merged? More...
TH1F *	f1stTo2ndPosLenHist
	Distance between cluster centers. More...
TH1F *	fRMaxFirstHist
	radius of "cylinder" first cluster More...
TH1F *	fCosMaxFirstHist
	Cos angle beteeen cylinder axis and edge. More...
TH1F *	fCosFirstAxisHist
	Cos angle between vector between centers and first axis. More...
TH1F *	f1stTo2ndProjEigenHist
	arc length along vector between centers from first center to cylinder More...
TH1F *	fRMaxNextHist
	radius of "cylinder" next cluster More...
TH1F *	fCosMaxNextHist
	Cos angle beteeen cylinder axis and edge. More...
TH1F *	fCosNextAxisHist
	Cos angle between vector between centers and next axis. More...
TH1F *	f2ndTo1stProjEigenHist
	arc length along vector between centers from next center to cylinder More...
TH1F *	fGapBetweenClusHist
	Gap between clusters. More...
TH1F *	fGapRatToLenHist
	Ratio of gap to distance between centers. More...
TH1F *	fProjEndPointLenHist
	Projection of vector between the endpoints. More...
TH1F *	fProjEndPointRatHist
	Ratio of projection to vector length. More...
TH1F *	fAxesDocaHist
	Closest distance between two primary axes. More...
TH1F *	f1stDocaArcLRatHist
	arc length to POCA for DOCA first cluster More...
TH1F *	f2ndDocaArcLRatHist
	arc length to POCA for DOCA second cluster More...
TH1F *	f1stTo2ndPosLenRatHist
	Ratio of distance between centers and first proj eigen. More...
TH1F *	fGapRatHist
	Ratio of gap and next proj eigen. More...
PrincipalComponentsAlg	fPCAAlg

Detailed Description

Definition at line 34 of file ClusterMergeAlg_tool.cc.

Constructor & Destructor Documentation

lar_cluster3d::ClusterMergeAlg::ClusterMergeAlg ( const fhicl::ParameterSet &  )

explicit

Constructor.

Parameters

pset

lar_cluster3d::ClusterMergeAlg::~ClusterMergeAlg

(

)

Destructor.

Definition at line 134 of file ClusterMergeAlg_tool.cc.

{
}

Member Function Documentation

float lar_cluster3d::ClusterMergeAlg::closestApproach ( const Eigen::Vector3f &  P0, const Eigen::Vector3f &  u0, const Eigen::Vector3f &  P1, const Eigen::Vector3f &  u1, Eigen::Vector3f &  poca0, Eigen::Vector3f &  poca1, Eigen::Vector3f &  firstNextUnit  ) const

private

Definition at line 580 of file ClusterMergeAlg_tool.cc.

{
    // Technique is to compute the arclength to each point of closest approach
    Eigen::Vector3f w0 = P0 - P1;
    float a(1.);
    float b(u0.dot(u1));
    float c(1.);
    float d(u0.dot(w0));
    float e(u1.dot(w0));
    float den(a * c - b * b);

    // Make sure lines are not colinear
    if (std::abs(den) > 10. * std::numeric_limits<float>::epsilon())
    {
        float arcLen0 = (b * e - c * d) / den;
        float arcLen1 = (a * e - b * d) / den;

        poca0 = P0 + arcLen0 * u0;
        poca1 = P1 + arcLen1 * u1;
    }
    // Otherwise, take the poca to be the distance to the line at the second point
    else
    {
        float arcLenToNextPoint = w0.dot(u0);

        poca0 = P0 + arcLenToNextPoint * u0;
        poca1 = P1;
    }

    firstNextUnit = poca1 - poca0;

    float docaDist = firstNextUnit.norm();

    firstNextUnit = firstNextUnit.normalized();

    return docaDist;
}

void lar_cluster3d::ClusterMergeAlg::configure ( fhicl::ParameterSet const &  pset )

override

Definition at line 140 of file ClusterMergeAlg_tool.cc.

{
    fEnableMonitoring     = pset.get<bool> ("EnableMonitoring",     true  );
    fMinTransEigenVal     = pset.get<float>("MinTransEigenVal",     0.09  );
    fMinEigenToProcess    = pset.get<float>("MinEigenToProcess",    2.0   );
    fOutputHistograms     = pset.get<bool> ("OutputHistograms",     false );

    fTimeToProcess = 0.;

    // If asked, define some histograms
    if (fOutputHistograms)
    {
        // Access ART's TFileService, which will handle creating and writing
        // histograms and n-tuples for us.
        art::ServiceHandle<art::TFileService> tfs;

        // Make a directory for these histograms
        art::TFileDirectory dir = tfs->mkdir("MergeClusters");

        fFirstEigenValueHists.resize(3,nullptr);
        fNextEigenValueHists.resize(3,nullptr);

        std::vector<float> maxValsVec = {20., 50., 250.};

        for(size_t idx : {0, 1, 2})
        {
            fFirstEigenValueHists[idx] = dir.make<TH1F>(Form("FEigen1st%1zu",idx),"Eigen Val", 200, 0., maxValsVec[idx]);
            fNextEigenValueHists[idx]  = dir.make<TH1F>(Form("FEigen2nd%1zu",idx),"Eigen Val", 200, 0., maxValsVec[idx]);
        }

        fNumMergedClusters      = dir.make<TH1F>("NumMergedClus",      "Number Merged",             200,    0., 1000.);

        f1stTo2ndPosLenHist     = dir.make<TH1F>("1stTo2ndPosLen",     "Distance between Clusters", 250,    0., 1000.);

        fRMaxFirstHist          = dir.make<TH1F>("rMaxFirst",          "Radius of First Cluster",   200,    0., 100.);
        fCosMaxFirstHist        = dir.make<TH1F>("CosMaxFirst",        "Cos Angle First Cyl/Axis",  200,    0., 1.  );
        fCosFirstAxisHist       = dir.make<TH1F>("CosFirstAxis",       "Cos Angle First Next/Axis", 200,    0., 1.  );
        f1stTo2ndProjEigenHist  = dir.make<TH1F>("1stTo2ndProjEigen",  "Projected Distance First",  200,    0., 200.);

        fRMaxNextHist           = dir.make<TH1F>("rMaxNext",           "Radius of Next Cluster",    200,    0., 100.);
        fCosMaxNextHist         = dir.make<TH1F>("CosMaxNext",         "Cos Angle Next Cyl/Axis",   200,    0., 1.  );
        fCosNextAxisHist        = dir.make<TH1F>("CosNextAxis",        "Cos Angle Next Next/Axis",  200,    0., 1.  );
        f2ndTo1stProjEigenHist  = dir.make<TH1F>("2ndTo1stProjEigen",  "Projected Distance Next",   200,    0., 200.);

        fGapBetweenClusHist     = dir.make<TH1F>("ClusterGap",         "Gap Between Clusters",      400, -200., 200.);
        fGapRatToLenHist        = dir.make<TH1F>("GapRatToLen",        "Ratio Gap to Distance",     100,   -8.,   2.);
        fProjEndPointLenHist    = dir.make<TH1F>("ProjEndPointLen",    "Projected End Point Len",   200, -100., 100.);
        fProjEndPointRatHist    = dir.make<TH1F>("ProjEndPointRat",    "Projected End Point Ratio", 100,    0.,   1.);

        fAxesDocaHist           = dir.make<TH1F>("AxesDocaHist",       "DOCA",                      200,    0.,  25.);
        f1stDocaArcLRatHist     = dir.make<TH1F>("ALenPOCA1Hist",      "Arc Len to POCA 1",         400,  -50.,  50.);
        f2ndDocaArcLRatHist     = dir.make<TH1F>("ALenPOCA2Hist",      "Arc Len to POCA 2",         400,  -50.,  50.);

        f1stTo2ndPosLenRatHist  = dir.make<TH1F>("1stTo2ndPosLenRat",  "Ratio clus dist to eigen",  200.,   0.,  20.);
        fGapRatHist             = dir.make<TH1F>("GapRat",             "Ratio Gap to Next Eigen",   400,  -20.,  20.);
    }

    return;
}

const reco::ClusterHit3D * lar_cluster3d::ClusterMergeAlg::findClosestHit3D ( const Eigen::Vector3f &  refPoint, const Eigen::Vector3f &  refVector, const reco::HitPairListPtr &  hitList  ) const

private

Definition at line 622 of file ClusterMergeAlg_tool.cc.

{
    const reco::ClusterHit3D* nearestHit3D(hitList.front());
    float                     closest(std::numeric_limits<float>::max());

    for(const auto& hit3D : hitList)
    {
        Eigen::Vector3f refToHitVec = hit3D->getPosition() - refPoint;
        float           arcLenToHit = refToHitVec.dot(refVector);

        if (arcLenToHit < closest)
        {
            nearestHit3D = hit3D;
            closest      = arcLenToHit;
        }
    }

    return nearestHit3D;
}

const reco::ClusterHit3D * lar_cluster3d::ClusterMergeAlg::findFurthestHit3D ( const Eigen::Vector3f &  refPoint, const Eigen::Vector3f &  refVector, const reco::HitPairListPtr &  hitList  ) const

private

Definition at line 642 of file ClusterMergeAlg_tool.cc.

{
    const reco::ClusterHit3D* nearestHit3D(hitList.front());
    float                     furthest(-std::numeric_limits<float>::max());

    for(const auto& hit3D : hitList)
    {
        Eigen::Vector3f refToHitVec = hit3D->getPosition() - refPoint;
        float           arcLenToHit = refToHitVec.dot(refVector);

        if (arcLenToHit > furthest)
        {
            nearestHit3D = hit3D;
            furthest     = arcLenToHit;
        }
    }

    return nearestHit3D;
}

float lar_cluster3d::ClusterMergeAlg::getTimeToExecute ( ) const

inlineoverridevirtual

If monitoring, recover the time to execute a particular function.

Implements lar_cluster3d::IClusterModAlg.

Definition at line 70 of file ClusterMergeAlg_tool.cc.

{return fTimeToProcess;}

void lar_cluster3d::ClusterMergeAlg::initializeHistograms ( art::TFileDirectory &  )

overridevirtual

Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder.

Parameters

TFileDirectory

- the folder to store the hists in

Implements lar_cluster3d::IClusterModAlg.

Definition at line 200 of file ClusterMergeAlg_tool.cc.

{
    return;
}

bool lar_cluster3d::ClusterMergeAlg::linearClusters ( reco::ClusterParameters &  firstCluster, reco::ClusterParameters &  nextCluster  ) const

private

Definition at line 324 of file ClusterMergeAlg_tool.cc.

{
    // Assume failure
    bool consistent(false);

    // The goal here is to compare the two input PCA's and determine if they are effectively colinear and
    // within reasonable range to consider merging them. Note that a key assumption is that the first input
    // PCA is from the "bigger" cluster, the second is "smaller" and may have a less reliable PCA.

    // Dereference the PCA's
    const reco::PrincipalComponents& firstPCA = firstCluster.getFullPCA();
    const reco::PrincipalComponents& nextPCA  = nextCluster.getFullPCA();

    // Recover the positions of the centers of the two clusters
    const Eigen::Vector3f& firstCenter = firstPCA.getAvePosition();
    const Eigen::Vector3f& nextCenter  = nextPCA.getAvePosition();

    // And form a vector between the two centers
    Eigen::Vector3f firstPosToNextPosVec  = nextCenter - firstCenter;
    Eigen::Vector3f firstPosToNextPosUnit = firstPosToNextPosVec.normalized();
    float           firstPosToNextPosLen  = firstPosToNextPosVec.norm();

    // Now get the first PCA's primary axis and since we'll use them get all of them at once...
    Eigen::Vector3f firstAxis0(firstPCA.getEigenVectors().row(0));
    Eigen::Vector3f firstAxis1(firstPCA.getEigenVectors().row(1));
    Eigen::Vector3f firstAxis2(firstPCA.getEigenVectors().row(2));

    // Will want the distance of closest approach of the next cluser's center to the primary, start by finding arc length
    float arcLenToNextDoca = firstPosToNextPosVec.dot(firstAxis2);

    // Adopt the convention that the cluster axis is in same direction as vector from first to next centers
    // And preserve the overall orientation of the PCA by flipping all if we flip the first one
    if (arcLenToNextDoca < 0.)
    {
        firstAxis0       = -firstAxis0;
        firstAxis1       = -firstAxis1;
        firstAxis2       = -firstAxis2;
        arcLenToNextDoca = -arcLenToNextDoca;
    }

    // Recover the eigen values of the first and second PCAs for selection cuts
    Eigen::Vector3f firstEigenVals(std::min(1.5 * sqrt(std::max(firstPCA.getEigenValues()[0],fMinTransEigenVal)),  20.),
                                   std::min(1.5 * sqrt(std::max(firstPCA.getEigenValues()[1],fMinTransEigenVal)),  50.),
                                   std::min(1.5 * sqrt(         firstPCA.getEigenValues()[2]),                    250.));

    // We treat the PCA as defining an elliptical tube and we use the projection of the vector between centers to
    // get the radius to the edge of the tube, from which we can determine the projected length inside the tube
    Eigen::Vector2f firstProj01Unit = Eigen::Vector2f(firstPosToNextPosUnit.dot(firstAxis0),firstPosToNextPosUnit.dot(firstAxis1)).normalized();
    float           firstEigen0Proj = firstEigenVals[0] * firstProj01Unit[0];
    float           firstEigen1Proj = firstEigenVals[1] * firstProj01Unit[1];
    float           rMaxFirst       = std::sqrt(firstEigen0Proj * firstEigen0Proj + firstEigen1Proj * firstEigen1Proj);
    float           cosMaxFirst     = firstEigenVals[2] / std::sqrt(firstEigenVals[2] * firstEigenVals[2] + rMaxFirst * rMaxFirst);
    float           cosFirstAxis    = firstAxis2.dot(firstPosToNextPosUnit);

    // Now calculate a measure of the length inside the tube along the vector between the cluster centers
    float firstEigenVal2       = firstEigenVals[2];
    float firstToNextProjEigen = firstEigenVal2;

    // There are two cases to consider, the first is that the vector exits out the side of the tube
    if (cosFirstAxis < cosMaxFirst)
    {
        // In this case we need the sign of the angle between the vector and the cluster primary axis
        float sinFirstAxis = std::sqrt(1. - cosFirstAxis * cosFirstAxis);

        // Here the length will be given by the radius, not the primary eigen value
        firstToNextProjEigen = rMaxFirst / sinFirstAxis;
    }
    else firstToNextProjEigen /= cosFirstAxis;

    // Get scale factor for selecting this pair
    float firstPosToNextPosScaleFactor = 8.;

    // A brief interlude to fill a few histograms
    if (fOutputHistograms)
    {
        f1stTo2ndPosLenHist->Fill(firstPosToNextPosLen, 1.);
        fRMaxFirstHist->Fill(rMaxFirst, 1.);
        fCosMaxFirstHist->Fill(cosMaxFirst, 1.);
        fCosFirstAxisHist->Fill(cosFirstAxis, 1.);
        f1stTo2ndProjEigenHist->Fill(firstToNextProjEigen, 1.);
    }

    // This makes first selection, it should eliminate most of the junk cases:
    if (firstPosToNextPosLen < firstPosToNextPosScaleFactor * firstToNextProjEigen)
    {
        // Recover the axes for the next PCA and make sure pointing convention is observed
        Eigen::Vector3f nextAxis0(nextPCA.getEigenVectors().row(0));
        Eigen::Vector3f nextAxis1(nextPCA.getEigenVectors().row(1));
        Eigen::Vector3f nextAxis2(nextPCA.getEigenVectors().row(2));

        // Recover the cos of the angle between the primary axes...
        float cosFirstNextAxis = firstAxis2.dot(nextAxis2);

        // And in this case we want to choose the sign of the next cluster's axes so that the
        // primary axes "point" in the same direction
        if (cosFirstNextAxis < 0.)
        {
            nextAxis0         = -nextAxis0;
            nextAxis1         = -nextAxis1;
            nextAxis2         = -nextAxis2;
            cosFirstNextAxis  = -cosFirstNextAxis;
        }

        // Get the eigen values again
        Eigen::Vector3f nextEigenVals(std::min(1.5 * sqrt(std::max(nextPCA.getEigenValues()[0],fMinTransEigenVal)),  20.),
                                      std::min(1.5 * sqrt(std::max(nextPCA.getEigenValues()[1],fMinTransEigenVal)),  50.),
                                      std::min(1.5 * sqrt(         nextPCA.getEigenValues()[2]),                    250.));

        // Repeat the calculation of the length of the vector through the cluster "tube"...
        Eigen::Vector2f nextProj01Unit = Eigen::Vector2f(firstPosToNextPosUnit.dot(nextAxis0),firstPosToNextPosUnit.dot(nextAxis1)).normalized();
        float           nextEigen0Proj = nextEigenVals[0] * nextProj01Unit[0];
        float           nextEigen1Proj = nextEigenVals[1] * nextProj01Unit[1];
        float           rMaxNext       = std::sqrt(nextEigen0Proj * nextEigen0Proj + nextEigen1Proj * nextEigen1Proj);
        float           cosMaxNext     = nextEigenVals[2] / std::sqrt(nextEigenVals[2] * nextEigenVals[2] + rMaxNext * rMaxNext);
        float           cosNextAxis    = std::abs(nextAxis2.dot(firstPosToNextPosUnit));

        // Now calculate a measure of the length inside the cylider along the vector between the cluster centers
        float nextToFirstProjEigen = nextEigenVals[2];

        // There are two cases to consider, the first is that the vector exits out the side of the cylinder
        if (cosNextAxis < cosMaxNext)
        {
            // In this case we need the sign of the angle between the vector and the cluster primary axis
            float sinNextAxis = std::sqrt(1. - cosNextAxis * cosNextAxis);

            // Here the length will be given by the radius, not the primary eigen value
            nextToFirstProjEigen = rMaxNext / sinNextAxis;
        }
        else nextToFirstProjEigen /= cosNextAxis;

        // Allow a generous gap but significantly derate as angle to next cluster increases
        //float nextToFirstScaleFactor = 6. * cosFirstNextAxis;
        float nextToFirstScaleFactor = 8. * cosFirstNextAxis;

        // Form the "gap" along the axis between centers from the projections of the eigen values
        // Note the gap can be negative
        float gapFirstToNext = firstPosToNextPosLen - firstToNextProjEigen - nextToFirstProjEigen;

        // Look at the presumed nearest endpoints of the two clusters
        Eigen::Vector3f firstEndPoint = firstCenter + firstEigenVals[2] * firstAxis2;
        Eigen::Vector3f nextTailPoint = nextCenter  - nextEigenVals[2]  * nextAxis2;
        Eigen::Vector3f endToTailVec  = nextTailPoint - firstEndPoint;

        // Get projection along vector between centers
        float endPointProjection = endToTailVec.dot(firstPosToNextPosUnit);
        float endToTailLen       = endToTailVec.norm();

        // Another brief interlude to fill some histograms
        if (fOutputHistograms)
        {
            for(size_t idx = 0; idx < 3; idx++) fNextEigenValueHists[idx]->Fill(nextEigenVals[idx],1.);

            fRMaxNextHist->Fill(rMaxNext, 1.);
            fCosMaxNextHist->Fill(cosMaxNext, 1.);
            fCosNextAxisHist->Fill(cosNextAxis, 1.);
            f2ndTo1stProjEigenHist->Fill(nextToFirstProjEigen, 1.);
            fGapBetweenClusHist->Fill(gapFirstToNext, 1.);
            fProjEndPointLenHist->Fill(endPointProjection, 1.);
            fProjEndPointRatHist->Fill(std::abs(endPointProjection)/endToTailLen, 1.);
        }

        // We mirror here the first selection above but now operate on the distance between centers less
        // the first cluster's projected eigen value
        if (gapFirstToNext < nextToFirstScaleFactor * nextToFirstProjEigen || (std::abs(endPointProjection) < nextToFirstProjEigen && endToTailLen < nextEigenVals[2]))
        {
            // Now check the distance of closest approach betweent the two vectors
            // Closest approach calculaiton results vectors
            Eigen::Vector3f firstPoca;
            Eigen::Vector3f nextPoca;
            Eigen::Vector3f firstToNextUnit;

            // Recover the doca of the two axes and their points of closest approach along each axis
            float lineDoca = closestApproach(firstCenter, firstAxis2, nextCenter, nextAxis2, firstPoca, nextPoca, firstToNextUnit);

            // Get the range through the first clusters "tube" for this doca vec
            Eigen::Vector3f firstPOCAProjUnit     = Eigen::Vector3f(firstToNextUnit.dot(firstAxis0),firstToNextUnit.dot(firstAxis1),firstToNextUnit.dot(firstAxis1)).normalized();
            float           firstPOCAVecProjEigen = std::sqrt(std::pow(firstEigenVals[0] * firstPOCAProjUnit[0],2)
                                                  +           std::pow(firstEigenVals[1] * firstPOCAProjUnit[1],2)
                                                  +           std::pow(firstEigenVals[2] * firstPOCAProjUnit[2],2));

            // Similarly for the next vector
            Eigen::Vector3f nextPOCAProjUnit     = Eigen::Vector3f(firstToNextUnit.dot(firstAxis0),firstToNextUnit.dot(firstAxis1),firstToNextUnit.dot(firstAxis1)).normalized();
            float           nextPOCAVecProjEigen = std::sqrt(std::pow(nextEigenVals[0] * nextPOCAProjUnit[0],2)
                                                 +           std::pow(nextEigenVals[1] * nextPOCAProjUnit[1],2)
                                                 +           std::pow(nextEigenVals[2] * nextPOCAProjUnit[2],2));

            if (fOutputHistograms)
            {
                // The below returned the signed arc lengths to their respective pocas
                float arcLenToFirstPoca = (firstPoca - firstCenter).dot(firstAxis2);
                float arcLenToNextPoca  = (nextPoca  - nextCenter ).dot(nextAxis2);

                fAxesDocaHist->Fill(lineDoca, 1.);
                f1stDocaArcLRatHist->Fill(arcLenToFirstPoca/firstEigenVals[2],1.);
                f2ndDocaArcLRatHist->Fill(arcLenToNextPoca/nextEigenVals[2],1.);
            }

            // Scaling factor to increase doca distance as distance grows
            float rMaxScaleFactor = 1.2;

            if (lineDoca < rMaxScaleFactor * (firstPOCAVecProjEigen + nextPOCAVecProjEigen))
            {
                consistent = true;

                if (fOutputHistograms)
                {
                    f1stTo2ndPosLenRatHist->Fill(firstPosToNextPosLen/firstToNextProjEigen, 1.);
                }
            }
        }
    }

    return consistent;
}

bool lar_cluster3d::ClusterMergeAlg::mergeClusters ( reco::ClusterParameters &  firstClusterParams, reco::ClusterParameters &  nextClusterParams  ) const

private

Definition at line 539 of file ClusterMergeAlg_tool.cc.

{
    bool merged(false);

    // Merge the next cluster into the first one
    // Get the hits
    reco::HitPairListPtr& hitPairListPtr = firstClusterParams.getHitPairListPtr();

    // We copy the hits from the old to new but note that we need to update the parameters for each 2D hit we add
    for(const auto* hit : nextClusterParams.getHitPairListPtr())
    {
        hitPairListPtr.push_back(hit);

        for(const auto* hit2D : hit->getHits())
            if (hit2D) firstClusterParams.UpdateParameters(hit2D);
    }

    // Recalculate the PCA
    fPCAAlg.PCAAnalysis_3D(firstClusterParams.getHitPairListPtr(), firstClusterParams.getFullPCA());

    // Must have a valid pca
    if (firstClusterParams.getFullPCA().getSvdOK())
    {
        // Finish out the merging here
        reco::Hit3DToEdgeMap& firstEdgeMap = firstClusterParams.getHit3DToEdgeMap();
        reco::Hit3DToEdgeMap& nextEdgeMap  = nextClusterParams.getHit3DToEdgeMap();

        for(const auto& pair : nextEdgeMap) firstEdgeMap[pair.first] = pair.second;

        // Set the skeleton PCA to make sure it has some value
        firstClusterParams.getSkeletonPCA() = firstClusterParams.getFullPCA();

        // Zap the cluster we merged into the new one...
        nextClusterParams = reco::ClusterParameters();

        merged = true;
    }

    return merged;
}

void lar_cluster3d::ClusterMergeAlg::ModifyClusters ( reco::ClusterParametersList &  clusterParametersList ) const

overridevirtual

Scan an input collection of clusters and modify those according to the specific implementing algorithm.

Parameters

clusterParametersList

A list of cluster objects (parameters from associated hits)

Top level interface for algorithm to consider pairs of clusters from the input list and determine if they are consistent with each other and, therefore, should be merged. This is done by looking at the PCA for each cluster and looking at the projection of the primary axis along the vector connecting their centers.

Implements lar_cluster3d::IClusterModAlg.

Definition at line 205 of file ClusterMergeAlg_tool.cc.

{
    /**
     *  @brief Top level interface for algorithm to consider pairs of clusters from the input
     *         list and determine if they are consistent with each other and, therefore, should
     *         be merged. This is done by looking at the PCA for each cluster and looking at the
     *         projection of the primary axis along the vector connecting their centers.
     */

    // Initial clustering is done, now trim the list and get output parameters
    cet::cpu_timer theClockBuildClusters;

    // Start clocks if requested
    if (fEnableMonitoring) theClockBuildClusters.start();

    // Resort by the largest first eigen value (so the "longest" cluster)
    clusterParametersList.sort([](auto& left, auto& right){return left.getFullPCA().getEigenValues()[2] > right.getFullPCA().getEigenValues()[2];});

    // The idea is to continually loop through all clusters until we get to the point where we are no longer doing any merging.
    // If two clusters are merged then we need to recycle through again because it may be that the new cluster can match when
    // it might not have in the past
    size_t                                lastClusterListCount = clusterParametersList.size() + 1;
    reco::ClusterParametersList::iterator lastFirstClusterItr  = clusterParametersList.begin();

    int numMergedClusters(0);
    int nOutsideLoops(0);

    while(clusterParametersList.size() != lastClusterListCount)
    {
        // Update the last count
        lastClusterListCount = clusterParametersList.size();

        // Keep track of the first cluster iterator each pass through
        reco::ClusterParametersList::iterator firstClusterItr = lastFirstClusterItr++;

        // Loop through the clusters
        while(firstClusterItr != clusterParametersList.end())
        {
            reco::ClusterParameters&              firstClusterParams = *firstClusterItr;
            reco::ClusterParametersList::iterator nextClusterItr     = firstClusterItr;

            // Take a brief interlude to do some histogramming
            if (fOutputHistograms)
            {
                const reco::PrincipalComponents& firstPCA = firstClusterParams.getFullPCA();

                Eigen::Vector3f firstEigenVals(std::min(1.5 * std::sqrt(std::max(firstPCA.getEigenValues()[0],fMinTransEigenVal)),  50.),
                                               std::min(1.5 * std::sqrt(std::max(firstPCA.getEigenValues()[1],fMinTransEigenVal)), 100.),
                                                        1.5 * std::sqrt(         firstPCA.getEigenValues()[2]));

                for(size_t idx = 0; idx < 3; idx++) fFirstEigenValueHists[idx]->Fill(firstEigenVals[idx],1.);
            }

            // Once you get down to the smallest clusters if they haven't already been absorbed there is no need to check them
            if (firstClusterParams.getFullPCA().getEigenValues()[2] < fMinEigenToProcess) break;

            while(++nextClusterItr != clusterParametersList.end())
            {
                reco::ClusterParameters& nextClusterParams = *nextClusterItr;

                // On any given loop through here it **might** be that the first cluster has been modified. So can't cache
                // the parameters, need the curret ones
                if (linearClusters(firstClusterParams,nextClusterParams))
                {
                    if (mergeClusters(firstClusterParams, nextClusterParams))
                    {
                        // Now remove the "next" cluster
                        nextClusterItr = clusterParametersList.erase(nextClusterItr);

                        // Our new cluster may be larger than those before it so we should try to move backwards to make sure to
                        // give now smaller clusters the chance to get merged as well
                        reco::ClusterParametersList::iterator biggestItr = firstClusterItr;

                        // Step backwards until we hit the beginning of the list
                        while(biggestItr-- != clusterParametersList.begin())
                        {
                            // The list has already been sorted largest to smallest so if we find a cluster that is "larger" then
                            // we have hit the limit
                            if ((*biggestItr).getFullPCA().getEigenValues()[2] > (*firstClusterItr).getFullPCA().getEigenValues()[2])
                            {
                                // Want to insert after the cluster with the larger primary axes.. but there is
                                // no point in doing anything if the new cluster is already in the right spot
                                if (++biggestItr != firstClusterItr) clusterParametersList.splice(biggestItr, clusterParametersList, firstClusterItr);

                                break;
                            }
                        }

                        // Restart loop?
                        nextClusterItr = firstClusterItr;

                        numMergedClusters++;
                    }
                }
            }

            firstClusterItr++;
        }

        nOutsideLoops++;
    }

    std::cout << "==> # merged: " << numMergedClusters << ", in " << nOutsideLoops << " outside loops " << std::endl;

    if (fOutputHistograms) fNumMergedClusters->Fill(numMergedClusters, 1.);


    if (fEnableMonitoring)
    {
        theClockBuildClusters.stop();

        fTimeToProcess = theClockBuildClusters.accumulated_real_time();
    }

    mf::LogDebug("Cluster3D") << ">>>>> Merge clusters done, found " << clusterParametersList.size() << " clusters" << std::endl;

    return;
}

Member Data Documentation

TH1F* lar_cluster3d::ClusterMergeAlg::f1stDocaArcLRatHist

private

arc length to POCA for DOCA first cluster

Definition at line 117 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::f1stTo2ndPosLenHist

private

Distance between cluster centers.

Definition at line 99 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::f1stTo2ndPosLenRatHist

private

Ratio of distance between centers and first proj eigen.

Definition at line 120 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::f1stTo2ndProjEigenHist

private

arc length along vector between centers from first center to cylinder

Definition at line 104 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::f2ndDocaArcLRatHist

private

arc length to POCA for DOCA second cluster

Definition at line 118 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::f2ndTo1stProjEigenHist

private

arc length along vector between centers from next center to cylinder

Definition at line 109 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fAxesDocaHist

private

Closest distance between two primary axes.

Definition at line 116 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fCosFirstAxisHist

private

Cos angle between vector between centers and first axis.

Definition at line 103 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fCosMaxFirstHist

private

Cos angle beteeen cylinder axis and edge.

Definition at line 102 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fCosMaxNextHist

private

Cos angle beteeen cylinder axis and edge.

Definition at line 107 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fCosNextAxisHist

private

Cos angle between vector between centers and next axis.

Definition at line 108 of file ClusterMergeAlg_tool.cc.

bool lar_cluster3d::ClusterMergeAlg::fEnableMonitoring

private

Data members to follow.

If true then turn on monitoring (e.g. timing)

Definition at line 87 of file ClusterMergeAlg_tool.cc.

std::vector<TH1F*> lar_cluster3d::ClusterMergeAlg::fFirstEigenValueHists

private

First Cluster eigen value.

Definition at line 95 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fGapBetweenClusHist

private

Gap between clusters.

Definition at line 111 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fGapRatHist

private

Ratio of gap and next proj eigen.

Definition at line 121 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fGapRatToLenHist

private

Ratio of gap to distance between centers.

Definition at line 112 of file ClusterMergeAlg_tool.cc.

float lar_cluster3d::ClusterMergeAlg::fMinEigenToProcess

private

Don't look anymore at clusters below this size.

Definition at line 89 of file ClusterMergeAlg_tool.cc.

float lar_cluster3d::ClusterMergeAlg::fMinTransEigenVal

private

Set a mininum allowed value for the transverse eigen values.

Definition at line 88 of file ClusterMergeAlg_tool.cc.

std::vector<TH1F*> lar_cluster3d::ClusterMergeAlg::fNextEigenValueHists

private

Next Cluster eigen value.

Definition at line 96 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fNumMergedClusters

private

How many clusters were merged?

Definition at line 97 of file ClusterMergeAlg_tool.cc.

bool lar_cluster3d::ClusterMergeAlg::fOutputHistograms

private

Take the time to create and fill some histograms for diagnostics.

Definition at line 90 of file ClusterMergeAlg_tool.cc.

PrincipalComponentsAlg lar_cluster3d::ClusterMergeAlg::fPCAAlg

private

Definition at line 123 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fProjEndPointLenHist

private

Projection of vector between the endpoints.

Definition at line 113 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fProjEndPointRatHist

private

Ratio of projection to vector length.

Definition at line 114 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fRMaxFirstHist

private

radius of "cylinder" first cluster

Definition at line 101 of file ClusterMergeAlg_tool.cc.

TH1F* lar_cluster3d::ClusterMergeAlg::fRMaxNextHist

private

radius of "cylinder" next cluster

Definition at line 106 of file ClusterMergeAlg_tool.cc.

float lar_cluster3d::ClusterMergeAlg::fTimeToProcess

mutableprivate

Keep track of how long it took to run this algorithm.

Definition at line 93 of file ClusterMergeAlg_tool.cc.

---

The documentation for this class was generated from the following file:

• larreco/larreco/RecoAlg/Cluster3DAlgs/ClusterMergeAlg_tool.cc