lar_cluster3d::Cluster3D Class Reference

Definition of the Cluster3D class. More...

Inheritance diagram for lar_cluster3d::Cluster3D:

Classes
class	ArtOutputHandler

Public Member Functions
	Cluster3D (fhicl::ParameterSet const &pset)

Private Types
using	IdxToPCAMap = std::map< size_t, const reco::PrincipalComponents * >
	Special routine to handle creating and saving space points & edges PCA points. More...

Private Member Functions
void	beginJob () override
void	produce (art::Event &evt) override
void	PrepareEvent (const art::Event &evt)
	Event Preparation. More...
void	InitializeMonitoring ()
	Initialize the internal monitoring. More...
void	findTrackSeeds (art::Event &evt, reco::ClusterParameters &cluster, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, std::vector< recob::Seed > &seedVec, art::Assns< recob::Seed, recob::Hit > &seedHitAssns) const
	An interface to the seed finding algorithm. More...
void	splitClustersWithMST (reco::ClusterParameters &clusterParameters, reco::ClusterParametersList &clusterParametersList) const
	Attempt to split clusters by using a minimum spanning tree. More...
void	splitClustersWithHough (reco::ClusterParameters &clusterParameters, reco::ClusterParametersList &clusterParametersList) const
	Attempt to split clusters using the output of the Hough Filter. More...
void	MakeAndSaveSpacePoints (ArtOutputHandler &output, std::vector< recob::SpacePoint > &spacePointVec, art::Assns< recob::Hit, recob::SpacePoint > &spHitAssns, reco::HitPairListPtr &clusHitPairVector, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, const std::string &path="") const
	Special routine to handle creating and saving space points. More...
void	MakeAndSaveKinkPoints (ArtOutputHandler &output, reco::ConvexHullKinkTupleList &clusHitPairVector) const
	Special routine to handle creating and saving space points. More...
void	MakeAndSaveVertexPoints (ArtOutputHandler &, dcel2d::VertexList &, dcel2d::HalfEdgeList &) const
	Special routine to handle creating and saving space points & edges associated to voronoi diagrams. More...
void	MakeAndSavePCAPoints (ArtOutputHandler &, const reco::PrincipalComponents &, IdxToPCAMap &) const
size_t	FindAndStoreDaughters (util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, IdxToPCAMap &idxToPCAMap, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, Hit3DToSPPtrMap &best3DToSPPtrMap) const
	This will produce art output for daughters noting that it needs to be done recursively. More...
void	ProduceArtClusters (util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::HitPairList &hitPairList, reco::ClusterParametersList &clusterParametersList, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap) const
	Top level output routine, allows checking cluster status. More...
size_t	ConvertToArtOutput (util::GeometryUtilities const &gser, ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, IHit3DBuilder::RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, Hit3DToSPPtrMap &best3DToSPPtrMap) const
	Produces the art output from all the work done in this producer module. More...
bool	aParallelHitsCluster (const reco::PrincipalComponents &pca) const
	There are several places we will want to know if a candidate cluster is a "parallel hits" type of cluster so encapsulate that here. More...
size_t	countUltimateDaughters (reco::ClusterParameters &clusterParameters) const
	Count number of end of line daughters. More...

Private Attributes
bool	m_onlyMakSpacePoints
	If true we don't do the full cluster 3D processing. More...
bool	m_enableMonitoring
	Turn on monitoring of this algorithm. More...
float	m_parallelHitsCosAng
	Cut for PCA 3rd axis angle to X axis. More...
float	m_parallelHitsTransWid
	Cut on transverse width of cluster (PCA 2nd eigenvalue) More...
TTree *	m_pRecoTree
int	m_run
int	m_event
int	m_hits
	Keeps track of the number of hits seen. More...
int	m_hits3D
	Keeps track of the number of 3D hits made. More...
float	m_totalTime
	Keeps track of total execution time. More...
float	m_artHitsTime
	Keeps track of time to recover hits. More...
float	m_makeHitsTime
	Keeps track of time to build 3D hits. More...
float	m_buildNeighborhoodTime
	Keeps track of time to build epsilon neighborhood. More...
float	m_dbscanTime
	Keeps track of time to run DBScan. More...
float	m_clusterMergeTime
	Keeps track of the time to merge clusters. More...
float	m_pathFindingTime
	Keeps track of the path finding time. More...
float	m_finishTime
	Keeps track of time to run output module. More...
std::string	m_pathInstance
	Special instance for path points. More...
std::string	m_vertexInstance
	Special instance name for vertex points. More...
std::string	m_extremeInstance
	Instance name for the extreme points. More...
std::unique_ptr < lar_cluster3d::IHit3DBuilder >	m_hit3DBuilderAlg
	Builds the 3D hits to operate on. More...
std::unique_ptr < lar_cluster3d::IClusterAlg >	m_clusterAlg
	Algorithm to do 3D space point clustering. More...
std::unique_ptr < lar_cluster3d::IClusterModAlg >	m_clusterMergeAlg
	Algorithm to do cluster merging. More...
std::unique_ptr < lar_cluster3d::IClusterModAlg >	m_clusterPathAlg
	Algorithm to do cluster path finding. More...
std::unique_ptr < lar_cluster3d::IClusterParametersBuilder >	m_clusterBuilder
	Common cluster builder tool. More...
PrincipalComponentsAlg	m_pcaAlg
	Principal Components algorithm. More...
SkeletonAlg	m_skeletonAlg
	Skeleton point finder. More...
HoughSeedFinderAlg	m_seedFinderAlg
	Seed finder. More...
PCASeedFinderAlg	m_pcaSeedFinderAlg
	Use PCA axis to find seeds. More...
ParallelHitsSeedFinderAlg	m_parallelHitsAlg
	Deal with parallel hits clusters. More...

Detailed Description

Definition of the Cluster3D class.

Definition at line 109 of file Cluster3D_module.cc.

Member Typedef Documentation

using lar_cluster3d::Cluster3D::IdxToPCAMap = std::map<size_t, const reco::PrincipalComponents*>

private

Special routine to handle creating and saving space points & edges PCA points.

Parameters

output	the object containting the art output
clusterParamsList	List of clusters to get PCA's from

Definition at line 409 of file Cluster3D_module.cc.

Constructor & Destructor Documentation

lar_cluster3d::Cluster3D::Cluster3D ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 537 of file Cluster3D_module.cc.

    : EDProducer{pset}
    , m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
    , m_skeletonAlg(pset.get<fhicl::ParameterSet>("SkeletonAlg"))
    , m_seedFinderAlg(pset.get<fhicl::ParameterSet>("SeedFinderAlg"))
    , m_pcaSeedFinderAlg(pset.get<fhicl::ParameterSet>("PCASeedFinderAlg"))
    , m_parallelHitsAlg(pset.get<fhicl::ParameterSet>("ParallelHitsAlg"))
  {
    m_onlyMakSpacePoints = pset.get<bool>("MakeSpacePointsOnly", false);
    m_enableMonitoring = pset.get<bool>("EnableMonitoring", false);
    m_parallelHitsCosAng = pset.get<float>("ParallelHitsCosAng", 0.999);
    m_parallelHitsTransWid = pset.get<float>("ParallelHitsTransWid", 25.0);
    m_pathInstance = pset.get<std::string>("PathPointsName", "Path");
    m_vertexInstance = pset.get<std::string>("VertexPointsName", "Vertex");
    m_extremeInstance = pset.get<std::string>("ExtremePointsName", "Extreme");

    m_hit3DBuilderAlg = art::make_tool<lar_cluster3d::IHit3DBuilder>(
      pset.get<fhicl::ParameterSet>("Hit3DBuilderAlg"));
    m_clusterAlg =
      art::make_tool<lar_cluster3d::IClusterAlg>(pset.get<fhicl::ParameterSet>("ClusterAlg"));
    m_clusterMergeAlg = art::make_tool<lar_cluster3d::IClusterModAlg>(
      pset.get<fhicl::ParameterSet>("ClusterMergeAlg"));
    m_clusterPathAlg = art::make_tool<lar_cluster3d::IClusterModAlg>(
      pset.get<fhicl::ParameterSet>("ClusterPathAlg"));
    m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
      pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));

    // Handle special case of Space Point building outputting a new hit collection
    m_hit3DBuilderAlg->produces(producesCollector());

    produces<std::vector<recob::PCAxis>>();
    produces<std::vector<recob::PFParticle>>();
    produces<std::vector<recob::Cluster>>();
    produces<std::vector<recob::Seed>>();
    produces<std::vector<recob::Edge>>();

    produces<art::Assns<recob::PFParticle, recob::PCAxis>>();
    produces<art::Assns<recob::PFParticle, recob::Cluster>>();
    produces<art::Assns<recob::PFParticle, recob::Seed>>();
    produces<art::Assns<recob::Edge, recob::PFParticle>>();
    produces<art::Assns<recob::Seed, recob::Hit>>();
    produces<art::Assns<recob::Cluster, recob::Hit>>();

    produces<std::vector<recob::SpacePoint>>();
    produces<art::Assns<recob::SpacePoint, recob::PFParticle>>();
    produces<art::Assns<recob::Hit, recob::SpacePoint>>();
    produces<art::Assns<recob::SpacePoint, recob::Edge>>();

    produces<std::vector<recob::SpacePoint>>(m_pathInstance);
    produces<std::vector<recob::Edge>>(m_pathInstance);
    produces<art::Assns<recob::SpacePoint, recob::PFParticle>>(m_pathInstance);
    produces<art::Assns<recob::Edge, recob::PFParticle>>(m_pathInstance);
    produces<art::Assns<recob::Hit, recob::SpacePoint>>(m_pathInstance);
    produces<art::Assns<recob::SpacePoint, recob::Edge>>(m_pathInstance);

    produces<std::vector<recob::Edge>>(m_vertexInstance);
    produces<std::vector<recob::SpacePoint>>(m_vertexInstance);

    produces<std::vector<recob::SpacePoint>>(m_extremeInstance);
  }

Member Function Documentation

bool lar_cluster3d::Cluster3D::aParallelHitsCluster ( const reco::PrincipalComponents &  pca ) const

inlineprivate

There are several places we will want to know if a candidate cluster is a "parallel hits" type of cluster so encapsulate that here.

Parameters

pca	The Principal Components Analysis parameters for the cluster

Definition at line 469 of file Cluster3D_module.cc.

    {
      return fabs(pca.getEigenVectors().row(2)(0)) > m_parallelHitsCosAng &&
             3. * sqrt(pca.getEigenValues()(1)) > m_parallelHitsTransWid;
    }

void lar_cluster3d::Cluster3D::beginJob ( )

overrideprivate

beginJob will be tasked with initializing monitoring, in necessary, but also to init the geometry and detector services (and this probably needs to go in a "beginEvent" method?)

Definition at line 601 of file Cluster3D_module.cc.

  {
    /**
     *  @brief beginJob will be tasked with initializing monitoring, in necessary, but also to init the
     *         geometry and detector services (and this probably needs to go in a "beginEvent" method?)
     */
    if (m_enableMonitoring) this->InitializeMonitoring();
  }

size_t lar_cluster3d::Cluster3D::ConvertToArtOutput ( util::GeometryUtilities const &  gser, ArtOutputHandler &  output, reco::ClusterParameters &  clusterParameters, size_t  pfParticleParent, IHit3DBuilder::RecobHitToPtrMap &  hitToPtrMap, Hit3DToSPPtrMap &  hit3DToSPPtrMap, Hit3DToSPPtrMap &  best3DToSPPtrMap  ) const

private

Produces the art output from all the work done in this producer module.

Parameters

output	the object containting the art output
clusterParameters	Cluster info to output (in internal format)
pfParticleParent	The parent ID reference for the output PFParticle
hitToPtrMap	This maps our Cluster2D hits back to art Ptr's to reco Hits

Definition at line 1446 of file Cluster3D_module.cc.

  {
    // prepare the algorithm to compute the cluster characteristics;
    // we use the "standard" one here, except that we override selected items
    // (so, thanks to metaprogramming, we finally have wrappers of wrappers);
    // configuration would happen here, but we are using the default
    // configuration for that algorithm
    using OverriddenClusterParamsAlg_t =
      cluster::OverriddenClusterParamsAlg<cluster::StandardClusterParamsAlg>;

    cluster::ClusterParamsImportWrapper<OverriddenClusterParamsAlg_t> ClusterParamAlgo;

    // It should be straightforward at this point to transfer information from our vector of clusters
    // to the larsoft objects... of course we still have some work to do first, in particular to
    // find the candidate seeds and their seed hits

    // We keep track of 2 PCA axes, the first is the "full" PCA run over all the 3D hits in the
    // candidate cluster. The second will be that derived from just using the "skeleton" hits.
    // Make a copy of the full PCA to keep that, then get a reference for the skeleton PCA
    reco::PrincipalComponents& fullPCA = clusterParameters.getFullPCA();
    reco::PrincipalComponents& skeletonPCA = clusterParameters.getSkeletonPCA();

    size_t clusterStart = output.artClusterVector->size();

    // Start loop over views to build out the hit lists and the 2D cluster objects
    for (const auto& clusParametersPair : clusterParameters.getClusterParams()) {
      const reco::RecobClusterParameters& clusParams = clusParametersPair.second;

      // Protect against a missing view
      if (clusParams.m_view == geo::kUnknown) continue;

      // Get a vector of hit pointers to seed the cluster algorithm
      std::vector<const recob::Hit*> recobHitVec(clusParams.m_hitVector.size());

      std::transform(clusParams.m_hitVector.begin(),
                     clusParams.m_hitVector.end(),
                     recobHitVec.begin(),
                     [](const auto& hit2D) { return hit2D->getHit(); });

      // Get the tdc/wire slope... from the unit vector...
      double startWire(clusParams.m_startWire);
      double endWire(clusParams.m_endWire);
      double startTime(clusParams.m_startTime);
      double endTime(clusParams.m_endTime);

      // feed the algorithm with all the cluster hits
      ClusterParamAlgo.ImportHits(gser, recobHitVec);

      // create the recob::Cluster directly in the vector
      cluster::ClusterCreator artCluster(gser,
                                         ClusterParamAlgo,                // algo
                                         startWire,                       // start_wire
                                         0.,                              // sigma_start_wire
                                         startTime,                       // start_tick
                                         clusParams.m_sigmaStartTime,     // sigma_start_tick
                                         endWire,                         // end_wire
                                         0.,                              // sigma_end_wire,
                                         endTime,                         // end_tick
                                         clusParams.m_sigmaEndTime,       // sigma_end_tick
                                         output.artClusterVector->size(), // ID
                                         clusParams.m_view,               // view
                                         clusParams.m_plane,              // plane
                                         recob::Cluster::Sentry           // sentry
      );

      output.artClusterVector->emplace_back(artCluster.move());

      // We love looping. In this case, our list of hits is comprised of "ClusterHits" and we need to get a RecobHitVector instead...
      RecobHitVector recobHits;

      for (const auto& hit2D : clusParams.m_hitVector) {
        if (hit2D == nullptr || hit2D->getHit() == nullptr) continue;

        art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit2D->getHit()];
        recobHits.push_back(hitPtr);
      }

      output.makeClusterHitAssns(recobHits);
    }

    // Last, let's try to get seeds for tracking..
    // Keep track of how many we have so far
    size_t numSeedsStart = output.artSeedVector->size();

    // Empty daughter vector for now
    std::vector<size_t> nullVector;
    size_t pfParticleIdx(output.artPFParticleVector->size());

    recob::PFParticle pfParticle(13, pfParticleIdx, pfParticleParent, nullVector);
    output.artPFParticleVector->push_back(pfParticle);

    // Assume that if we are a daughter particle then there is a set of "best" 3D points and the convex hull will have been calculated from them
    // If there is no best list then assume the convex hull is from all of the points...
    std::string instance("");
    reco::HitPairListPtr* hit3DListPtr = &clusterParameters.getHitPairListPtr();
    Hit3DToSPPtrMap* hit3DToPtrMap = &hit3DToSPPtrMap;

    auto spVector = output.artSpacePointVector.get();
    auto edgeVector = output.artEdgeVector.get();
    auto pfPartEdgeAssns = output.artPFPartEdgeAssociations.get();
    auto spEdgeAssns = output.artEdgeSPAssociations.get();
    auto spHitAssns = output.artSPHitAssociations.get();
    auto pfPartSPAssns = output.artPFPartSPAssociations.get();

    // Make associations to all space points for this cluster
    int spacePointStart = spVector->size();

    MakeAndSaveSpacePoints(
      output, *spVector, *spHitAssns, *hit3DListPtr, hitToPtrMap, *hit3DToPtrMap);

    // And make sure to have associations to the PFParticle
    output.makePFPartSpacePointAssns(*spVector, *pfPartSPAssns, spacePointStart);

    // If they exist, make space points for the Path points
    if (!clusterParameters.getBestHitPairListPtr().empty()) {
      instance = m_pathInstance;
      hit3DListPtr = &clusterParameters.getBestHitPairListPtr();
      hit3DToPtrMap = &best3DToSPPtrMap;

      spVector = output.artPathPointVector.get();
      edgeVector = output.artPathEdgeVector.get();
      pfPartEdgeAssns = output.artPFPartPathEdgeAssociations.get();
      spEdgeAssns = output.artEdgePPAssociations.get();
      spHitAssns = output.artPPHitAssociations.get();
      pfPartSPAssns = output.artPFPartPPAssociations.get();

      spacePointStart = spVector->size();

      MakeAndSaveSpacePoints(
        output, *spVector, *spHitAssns, *hit3DListPtr, hitToPtrMap, *hit3DToPtrMap, instance);

      // And make sure to have associations to the PFParticle
      output.makePFPartSpacePointAssns(*spVector, *pfPartSPAssns, spacePointStart, instance);
    }

    // Now handle the edges according to whether associated with regular or "best" space points
    if (!clusterParameters.getBestEdgeList().empty()) {
      size_t edgeStart = edgeVector->size();

      for (const auto& edge : clusterParameters.getBestEdgeList()) {
        RecobSpacePointVector spacePointVec;

        try {
          spacePointVec.push_back(hit3DToPtrMap->at(std::get<0>(edge)));
          spacePointVec.push_back(hit3DToPtrMap->at(std::get<1>(edge)));
        }
        catch (...) {
          std::cout << "Caught exception in looking up space point ptr... " << std::get<0>(edge)
                    << ", " << std::get<1>(edge) << std::endl;
          continue;
        }

        edgeVector->emplace_back(
          std::get<2>(edge), spacePointVec[0].key(), spacePointVec[1].key(), edgeVector->size());

        output.makeEdgeSpacePointAssns(*edgeVector, spacePointVec, *spEdgeAssns, instance);
      }

      output.makePFPartEdgeAssns(*edgeVector, *pfPartEdgeAssns, edgeStart, instance);
    }

    // Look at making the PCAxis and associations - for both the skeleton (the first) and the full
    // First need some float to double conversion containers
    recob::PCAxis::EigenVectors eigenVecs;
    double eigenVals[] = {0., 0., 0.};
    double avePosition[] = {0., 0., 0.};

    eigenVecs.resize(3);

    for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
      avePosition[outerIdx] = skeletonPCA.getAvePosition()[outerIdx];
      eigenVals[outerIdx] = skeletonPCA.getEigenValues()[outerIdx];

      eigenVecs[outerIdx].resize(3);

      for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
        eigenVecs[outerIdx][innerIdx] = skeletonPCA.getEigenVectors().row(outerIdx)(innerIdx);
    }

    recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
                             skeletonPCA.getNumHitsUsed(),
                             eigenVals,
                             eigenVecs,
                             avePosition,
                             skeletonPCA.getAveHitDoca(),
                             output.artPCAxisVector->size());

    output.artPCAxisVector->push_back(skelPcAxis);

    for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
      avePosition[outerIdx] = fullPCA.getAvePosition()[outerIdx];
      eigenVals[outerIdx] = fullPCA.getEigenValues()[outerIdx];

      for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
        eigenVecs[outerIdx][innerIdx] = fullPCA.getEigenVectors().row(outerIdx)(innerIdx);
    }

    recob::PCAxis fullPcAxis(fullPCA.getSvdOK(),
                             fullPCA.getNumHitsUsed(),
                             eigenVals,   //fullPCA.getEigenValues(),
                             eigenVecs,   //fullPCA.getEigenVectors(),
                             avePosition, //fullPCA.getAvePosition(),
                             fullPCA.getAveHitDoca(),
                             output.artPCAxisVector->size());

    output.artPCAxisVector->push_back(fullPcAxis);

    // Create associations to the PFParticle
    output.makePFPartPCAAssns();
    output.makePFPartSeedAssns(numSeedsStart);
    output.makePFPartClusterAssns(clusterStart);

    return pfParticleIdx;
  }

size_t lar_cluster3d::Cluster3D::countUltimateDaughters ( reco::ClusterParameters &  clusterParameters ) const

private

Count number of end of line daughters.

Parameters

clusterParams

input cluster parameters to look at

Definition at line 1393 of file Cluster3D_module.cc.

  {
    size_t localCount(0);

    if (!clusterParameters.daughterList().empty()) {
      for (auto& clusterParams : clusterParameters.daughterList())
        localCount += countUltimateDaughters(clusterParams);
    }
    else
      localCount++;

    return localCount;
  }

size_t lar_cluster3d::Cluster3D::FindAndStoreDaughters ( util::GeometryUtilities const &  gser, ArtOutputHandler &  output, reco::ClusterParameters &  clusterParameters, size_t  pfParticleParent, IdxToPCAMap &  idxToPCAMap, IHit3DBuilder::RecobHitToPtrMap &  hitToPtrMap, Hit3DToSPPtrMap &  hit3DToSPPtrMap, Hit3DToSPPtrMap &  best3DToSPPtrMap  ) const

private

This will produce art output for daughters noting that it needs to be done recursively.

Parameters

output	the object containting the art output
clusterParameters	Cluster info to output (in internal format)
pfParticleParent	The parent ID reference for the output PFParticle
daughterList	List of PFParticle indices for stored daughters
hitToPtrMap	This maps our Cluster2D hits back to art Ptr's to reco Hits

Definition at line 1408 of file Cluster3D_module.cc.

  {
    // This is a recursive routine, we keep calling ourself as long as the daughter list is non empty
    if (!clusterParameters.daughterList().empty()) {
      for (auto& clusterParams : clusterParameters.daughterList())
        FindAndStoreDaughters(gser,
                              output,
                              clusterParams,
                              pfParticleParent,
                              idxToPCAMap,
                              hitToPtrMap,
                              hit3DToSPPtrMap,
                              best3DToSPPtrMap);
    }
    // Otherwise we want to store the information
    else {
      size_t daughterIdx = ConvertToArtOutput(gser,
                                              output,
                                              clusterParameters,
                                              pfParticleParent,
                                              hitToPtrMap,
                                              hit3DToSPPtrMap,
                                              best3DToSPPtrMap);

      idxToPCAMap[daughterIdx] = &clusterParameters.getFullPCA();
    }

    return idxToPCAMap.size();
  }

void lar_cluster3d::Cluster3D::findTrackSeeds ( art::Event &  evt, reco::ClusterParameters &  cluster, IHit3DBuilder::RecobHitToPtrMap &  hitToPtrMap, std::vector< recob::Seed > &  seedVec, art::Assns< recob::Seed, recob::Hit > &  seedHitAssns  ) const

private

An interface to the seed finding algorithm.

Parameters

evt	the ART event
cluster	structure of information representing a single cluster
hitToPtrMap	This maps our Cluster2D hits back to art Ptr's to reco Hits
seedVec	the output vector of candidate seeds
seedHitAssns	the associations between the seeds and the 2D hits making them

This method provides an interface to various algorithms for finding candiate recob::Seed objects and, as well, their candidate related seed hits

Definition at line 743 of file Cluster3D_module.cc.

  {
    /**
     *  @brief This method provides an interface to various algorithms for finding candiate
     *         recob::Seed objects and, as well, their candidate related seed hits
     */

    // Make sure we are using the right pca
    reco::PrincipalComponents& fullPCA = cluster.getFullPCA();
    reco::PrincipalComponents& skeletonPCA = cluster.getSkeletonPCA();
    reco::HitPairListPtr& hitPairListPtr = cluster.getHitPairListPtr();
    reco::HitPairListPtr skeletonListPtr;

    // We want to work with the "skeleton" hits so first step is to call the algorithm to
    // recover only these hits from the entire input collection
    m_skeletonAlg.GetSkeletonHits(hitPairListPtr, skeletonListPtr);

    // Skeleton hits are nice but we can do better if we then make a pass through to "average"
    // the skeleton hits position in the Y-Z plane
    m_skeletonAlg.AverageSkeletonPositions(skeletonListPtr);

    SeedHitPairListPairVec seedHitPairVec;

    // Some combination of the elements below will be used to determine which seed finding algorithm
    // to pursue below
    float eigenVal0 = 3. * sqrt(skeletonPCA.getEigenValues()[0]);
    float eigenVal1 = 3. * sqrt(skeletonPCA.getEigenValues()[1]);
    float eigenVal2 = 3. * sqrt(skeletonPCA.getEigenValues()[2]);
    float transRMS = std::hypot(eigenVal0, eigenVal1);

    bool foundGoodSeed(false);

    // Choose a method for finding the seeds based on the PCA that was run...
    // Currently we have an ad hoc if-else block which I hope will be improved soon!
    if (aParallelHitsCluster(fullPCA)) {
      // In this case we have a track moving relatively parallel to the wire plane with lots of
      // ambiguous 3D hits. Your best bet here is to use the "parallel hits" algorithm to get the
      // best axis and seeds
      // This algorithm does not fail (foundGoodSeed will always return true)
      foundGoodSeed = m_parallelHitsAlg.findTrackSeeds(hitPairListPtr, skeletonPCA, seedHitPairVec);
    }
    else if (eigenVal2 > 40. && transRMS < 5.) {
      // If the input cluster is relatively "straight" then chances are it is a single straight track,
      // probably a CR muon, and we can simply use the PCA to determine the seed
      // This algorithm will check both "ends" of the input hits and if the angles become inconsistent
      // then it will "fail"
      foundGoodSeed =
        m_pcaSeedFinderAlg.findTrackSeeds(skeletonListPtr, skeletonPCA, seedHitPairVec);
    }

    // In the event the above two methods failed then we hit it with the real seed finder
    if (!foundGoodSeed) {
      // If here then we have a complicated 3D cluster and we'll use the hough transform algorithm to
      // return a list of candidate seeds and seed hits
      m_seedFinderAlg.findTrackSeeds(skeletonListPtr, skeletonPCA, seedHitPairVec);
    }

    // Go through the returned lists and build out the art friendly seeds and hits
    for (const auto& seedHitPair : seedHitPairVec) {
      seedVec.push_back(seedHitPair.first);

      // We use a set here because our 3D hits can share 2D hits
      // The set will make sure we get unique combinations of 2D hits
      std::set<art::Ptr<recob::Hit>> seedHitSet;
      for (const auto& hit3D : seedHitPair.second) {
        for (const auto& hit2D : hit3D->getHits()) {
          if (!hit2D) continue;

          const recob::Hit* recobHit = hit2D->getHit();
          seedHitSet.insert(hitToPtrMap[recobHit]);
        }
      }

      RecobHitVector seedHitVec;
      for (const auto& hit2D : seedHitSet)
        seedHitVec.push_back(hit2D);

      util::CreateAssn(evt, seedVec, seedHitVec, seedHitAssns);
    }
  }

void lar_cluster3d::Cluster3D::InitializeMonitoring ( )

private

Initialize the internal monitoring.

Definition at line 702 of file Cluster3D_module.cc.

  {
    art::ServiceHandle<art::TFileService> tfs;
    m_pRecoTree = tfs->make<TTree>("monitoring", "LAr Reco");
    m_pRecoTree->Branch("run", &m_run, "run/I");
    m_pRecoTree->Branch("event", &m_event, "event/I");
    m_pRecoTree->Branch("hits", &m_hits, "hits/I");
    m_pRecoTree->Branch("hits3D", &m_hits3D, "hits3D/I");
    m_pRecoTree->Branch("totalTime", &m_totalTime, "time/F");
    m_pRecoTree->Branch("artHitsTime", &m_artHitsTime, "time/F");
    m_pRecoTree->Branch("makeHitsTime", &m_makeHitsTime, "time/F");
    m_pRecoTree->Branch("buildneigborhoodTime", &m_buildNeighborhoodTime, "time/F");
    m_pRecoTree->Branch("dbscanTime", &m_dbscanTime, "time/F");
    m_pRecoTree->Branch("clusterMergeTime", &m_clusterMergeTime, "time/F");
    m_pRecoTree->Branch("pathfindingtime", &m_pathFindingTime, "time/F");
    m_pRecoTree->Branch("finishTime", &m_finishTime, "time/F");

    m_clusterPathAlg->initializeHistograms(*tfs.get());
  }

void lar_cluster3d::Cluster3D::MakeAndSaveKinkPoints ( ArtOutputHandler &  output, reco::ConvexHullKinkTupleList &  clusHitPairVector  ) const

private

Special routine to handle creating and saving space points.

Parameters

output	the object containting the art output
clusHitPairVector	List of 3D hits to output as "extreme" space points

Definition at line 1726 of file Cluster3D_module.cc.

  {
    // Right now error matrix is uniform...
    double spError[] = {1., 0., 1., 0., 0., 1.};

    // Copy these hits to the vector to be stored with the event
    for (auto& kinkTuple : kinkTupleVec) {
      const reco::ClusterHit3D* hit = std::get<2>(std::get<0>(kinkTuple));

      double chisq = hit->getHitChiSquare(); // secret handshake...

      // Create and store the space point
      double spacePointPos[] = {
        hit->getPosition()[0], hit->getPosition()[1], hit->getPosition()[2]};

      spError[1] = hit->getTotalCharge();
      spError[3] = hit->getChargeAsymmetry();

      output.artExtremePointVector->push_back(
        recob::SpacePoint(spacePointPos, spError, chisq, output.artExtremePointVector->size()));
    }

    return;
  }

void lar_cluster3d::Cluster3D::MakeAndSavePCAPoints ( ArtOutputHandler &  output, const reco::PrincipalComponents &  fullPCA, IdxToPCAMap &  idxToPCAMap  ) const

private

Definition at line 1819 of file Cluster3D_module.cc.

  {
    // We actually do two things here:
    // 1) Create space points from the centroids of the PCA for each cluster
    // 2) Create the edges that link the space points together

    // The first task is to put the list of PCA's into some semblance of order... they may be
    // preordered by likely they are piecewise ordered so fix that here

    // We'll need the current PCA axis to determine doca and arclen
    Eigen::Vector3f avePosition(
      fullPCA.getAvePosition()[0], fullPCA.getAvePosition()[1], fullPCA.getAvePosition()[2]);
    Eigen::Vector3f axisDirVec(fullPCA.getEigenVectors().row(2));

    using DocaToPCAPair = std::pair<float, const reco::PrincipalComponents*>;
    using DocaToPCAVec = std::vector<DocaToPCAPair>;

    DocaToPCAVec docaToPCAVec;

    // Outer loop over views
    for (const auto& idxToPCA : idxToPCAMap) {
      const reco::PrincipalComponents* pca = idxToPCA.second;

      // Now we need to calculate the doca and poca...
      // Start by getting this hits position
      Eigen::Vector3f pcaPos(
        pca->getAvePosition()[0], pca->getAvePosition()[1], pca->getAvePosition()[2]);

      // Form a TVector from this to the cluster average position
      Eigen::Vector3f pcaToAveVec = pcaPos - avePosition;

      // With this we can get the arclength to the doca point
      float arclenToPoca = pcaToAveVec.dot(axisDirVec);

      docaToPCAVec.emplace_back(arclenToPoca, pca);
    }

    std::sort(docaToPCAVec.begin(), docaToPCAVec.end(), [](const auto& left, const auto& right) {
      return left.first < right.first;
    });

    // Set up the space point creation
    // Right now error matrix is uniform...
    double spError[] = {1., 0., 1., 0., 0., 1.};
    double chisq = 1.;

    const reco::PrincipalComponents* lastPCA(nullptr);

    // Set up to loop through the clusters
    for (const auto& docaToPCAPair : docaToPCAVec) {
      // Recover the PCA for this cluster
      const reco::PrincipalComponents* curPCA = docaToPCAPair.second;

      if (lastPCA) {
        double lastPointPos[] = {
          lastPCA->getAvePosition()[0], lastPCA->getAvePosition()[1], lastPCA->getAvePosition()[2]};
        size_t lastPointBin = output.artVertexPointVector->size();
        double curPointPos[] = {
          curPCA->getAvePosition()[0], curPCA->getAvePosition()[1], curPCA->getAvePosition()[2]};
        size_t curPointBin = lastPointBin + 1;

        output.artVertexPointVector->emplace_back(lastPointPos, spError, chisq, lastPointBin);
        output.artVertexPointVector->emplace_back(curPointPos, spError, chisq, curPointBin);

        Eigen::Vector3f distVec(curPointPos[0] - lastPointPos[0],
                                curPointPos[1] - lastPointPos[1],
                                curPointPos[2] - lastPointPos[2]);

        output.artVertexEdgeVector->emplace_back(
          distVec.norm(), lastPointBin, curPointBin, output.artEdgeVector->size());
      }

      lastPCA = curPCA;
    }
  }

void lar_cluster3d::Cluster3D::MakeAndSaveSpacePoints ( ArtOutputHandler &  output, std::vector< recob::SpacePoint > &  spacePointVec, art::Assns< recob::Hit, recob::SpacePoint > &  spHitAssns, reco::HitPairListPtr &  clusHitPairVector, IHit3DBuilder::RecobHitToPtrMap &  hitToPtrMap, Hit3DToSPPtrMap &  hit3DToSPPtrMap, const std::string &  path = ""  ) const

private

Special routine to handle creating and saving space points.

Parameters

output	the object containting the art output
clusterParameters	Cluster info to output (in internal format)
pfParticleParent	The parent ID reference for the output PFParticle
hitToPtrMap	This maps our Cluster2D hits back to art Ptr's to reco Hits

Definition at line 1668 of file Cluster3D_module.cc.

  {
    // Reserve space...
    spacePointVec.reserve(spacePointVec.size() + clusHitPairVector.size());

    // Right now error matrix is uniform...
    double spError[] = {1., 0., 1., 0., 0., 1.};

    // Copy these hits to the vector to be stored with the event
    for (auto& hitPair : clusHitPairVector) {
      // Skip those space points that have already been created
      if (hit3DToSPPtrMap.find(hitPair) != hit3DToSPPtrMap.end()) continue;

      // Don't make space point if this hit was "rejected"
      if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;

      double chisq = hitPair->getHitChiSquare(); // secret handshake...

      // Mark this hit pair as in use
      hitPair->setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);

      // Create and store the space point
      size_t spacePointID = spacePointVec.size();
      double spacePointPos[] = {
        hitPair->getPosition()[0], hitPair->getPosition()[1], hitPair->getPosition()[2]};

      spError[1] = hitPair->getTotalCharge();
      spError[3] = hitPair->getChargeAsymmetry();

      spacePointVec.emplace_back(spacePointPos, spError, chisq, spacePointID);

      // Update mapping
      hit3DToSPPtrMap[hitPair] = output.makeSpacePointPtr(spacePointID, instance);

      // space point hits associations
      RecobHitVector recobHits;

      for (const auto& hit : hitPair->getHits()) {
        if (!hit) continue;

        art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit->getHit()];
        recobHits.push_back(hitPtr);
      }

      if (!recobHits.empty())
        output.makeSpacePointHitAssns(spacePointVec, recobHits, spHitAssns, instance);
    }

    return;
  }

void lar_cluster3d::Cluster3D::MakeAndSaveVertexPoints ( ArtOutputHandler &  output, dcel2d::VertexList &  vertexList, dcel2d::HalfEdgeList &  halfEdgeList  ) const

private

Special routine to handle creating and saving space points & edges associated to voronoi diagrams.

Parameters

output	the object containting the art output
vertexList	list of vertices in the diagram
HalfEdgeList	list of half edges in the diagram

Definition at line 1753 of file Cluster3D_module.cc.

  {
    // We actually do two things here:
    // 1) Create space points to represent the vertex locations of the voronoi diagram
    // 2) Create the edges that link the space points together

    // Set up the space point creation
    // Right now error matrix is uniform...
    double spError[] = {1., 0., 1., 0., 0., 1.};
    double chisq = 1.;

    // Keep track of the vertex to space point association
    std::map<const dcel2d::Vertex*, size_t> vertexToSpacePointMap;

    // Copy these hits to the vector to be stored with the event
    for (auto& vertex : vertexList) {
      const dcel2d::Coords& coords = vertex.getCoords();

      // Create and store the space point
      double spacePointPos[] = {coords[0], coords[1], coords[2]};

      vertexToSpacePointMap[&vertex] = output.artVertexPointVector->size();

      output.artVertexPointVector->emplace_back(
        spacePointPos, spError, chisq, output.artVertexPointVector->size());
    }

    // Try to avoid double counting
    std::set<const dcel2d::HalfEdge*> halfEdgeSet;

    // Build the edges now
    for (const auto& halfEdge : halfEdgeList) {
      // Recover twin
      const dcel2d::HalfEdge* twin = halfEdge.getTwinHalfEdge();

      // It can happen that we have no twin... and also check that we've not been here before
      if (twin && halfEdgeSet.find(twin) == halfEdgeSet.end()) {
        // Recover the vertices
        const dcel2d::Vertex* fromVertex = twin->getTargetVertex();
        const dcel2d::Vertex* toVertex = halfEdge.getTargetVertex();

        // It can happen for the open edges that there is no target vertex
        if (!toVertex || !fromVertex) continue;

        if (vertexToSpacePointMap.find(fromVertex) == vertexToSpacePointMap.end() ||
            vertexToSpacePointMap.find(toVertex) == vertexToSpacePointMap.end())
          continue;

        // Need the distance between vertices
        Eigen::Vector3f distVec = toVertex->getCoords() - fromVertex->getCoords();

        output.artVertexEdgeVector->emplace_back(distVec.norm(),
                                                 vertexToSpacePointMap[fromVertex],
                                                 vertexToSpacePointMap[toVertex],
                                                 output.artEdgeVector->size());

        halfEdgeSet.insert(&halfEdge);
      }
    }

    return;
  }

void lar_cluster3d::Cluster3D::PrepareEvent ( const art::Event &  evt )

private

Event Preparation.

Parameters

evt	the ART event

Definition at line 725 of file Cluster3D_module.cc.

  {
    m_run = evt.run();
    m_event = evt.id().event();
    m_hits = 0;
    m_hits3D = 0;
    m_totalTime = 0.f;
    m_artHitsTime = 0.f;
    m_makeHitsTime = 0.f;
    m_buildNeighborhoodTime = 0.f;
    m_dbscanTime = 0.f;
    m_pathFindingTime = 0.f;
    m_finishTime = 0.f;
  }

void lar_cluster3d::Cluster3D::produce ( art::Event &  evt )

overrideprivate

Definition at line 613 of file Cluster3D_module.cc.

  {
    mf::LogInfo("Cluster3D") << " *** Cluster3D::produce(...)  [Run=" << evt.run()
                             << ", Event=" << evt.id().event() << "] Starting Now! *** "
                             << std::endl;

    // Set up for monitoring the timing... at some point this should be removed in favor of
    // external profilers
    cet::cpu_timer theClockTotal;
    cet::cpu_timer theClockFinish;

    if (m_enableMonitoring) theClockTotal.start();

    // This really only does anything if we are monitoring since it clears our tree variables
    this->PrepareEvent(evt);

    // Get instances of the primary data structures needed
    reco::ClusterParametersList clusterParametersList;
    IHit3DBuilder::RecobHitToPtrMap clusterHitToArtPtrMap;
    std::unique_ptr<reco::HitPairList> hitPairList(
      new reco::HitPairList); // Potentially lots of hits, use heap instead of stack

    // Call the algorithm that builds 3D hits and stores the hit collection
    m_hit3DBuilderAlg->Hit3DBuilder(evt, *hitPairList, clusterHitToArtPtrMap);

    // Only do the rest if we are not in the mode of only building space points (requested by ML folks)
    if (!m_onlyMakSpacePoints) {
      // Call the main workhorse algorithm for building the local version of candidate 3D clusters
      m_clusterAlg->Cluster3DHits(*hitPairList, clusterParametersList);

      // Try merging clusters
      m_clusterMergeAlg->ModifyClusters(clusterParametersList);

      // Run the path finding
      m_clusterPathAlg->ModifyClusters(clusterParametersList);
    }

    if (m_enableMonitoring) theClockFinish.start();

    // Get the art ouput object
    ArtOutputHandler output(evt, m_pathInstance, m_vertexInstance, m_extremeInstance);

    // Call the module that does the end processing (of which there is quite a bit of work!)
    // This goes here to insure that something is always written to the data store
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const detProp =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
    util::GeometryUtilities const gser{*lar::providerFrom<geo::Geometry>(), clockData, detProp};
    ProduceArtClusters(gser, output, *hitPairList, clusterParametersList, clusterHitToArtPtrMap);

    // Output to art
    output.outputObjects();

    // If monitoring then deal with the fallout
    if (m_enableMonitoring) {
      theClockFinish.stop();
      theClockTotal.stop();

      m_run = evt.run();
      m_event = evt.id().event();
      m_totalTime = theClockTotal.accumulated_real_time();
      m_artHitsTime = m_hit3DBuilderAlg->getTimeToExecute(IHit3DBuilder::COLLECTARTHITS);
      m_makeHitsTime = m_hit3DBuilderAlg->getTimeToExecute(IHit3DBuilder::BUILDTHREEDHITS);
      m_buildNeighborhoodTime = m_clusterAlg->getTimeToExecute(IClusterAlg::BUILDHITTOHITMAP);
      m_dbscanTime = m_clusterAlg->getTimeToExecute(IClusterAlg::RUNDBSCAN) +
                     m_clusterAlg->getTimeToExecute(IClusterAlg::BUILDCLUSTERINFO);
      m_clusterMergeTime = m_clusterMergeAlg->getTimeToExecute();
      m_pathFindingTime = m_clusterPathAlg->getTimeToExecute();
      m_finishTime = theClockFinish.accumulated_real_time();
      m_hits = static_cast<int>(clusterHitToArtPtrMap.size());
      m_hits3D = static_cast<int>(hitPairList->size());
      m_pRecoTree->Fill();

      mf::LogDebug("Cluster3D") << "*** Cluster3D total time: " << m_totalTime
                                << ", art: " << m_artHitsTime << ", make: " << m_makeHitsTime
                                << ", build: " << m_buildNeighborhoodTime
                                << ", clustering: " << m_dbscanTime
                                << ", merge: " << m_clusterMergeTime
                                << ", path: " << m_pathFindingTime << ", finish: " << m_finishTime
                                << std::endl;
    }

    // Will we ever get here? ;-)
    return;
  }

void lar_cluster3d::Cluster3D::ProduceArtClusters ( util::GeometryUtilities const &  gser, ArtOutputHandler &  output, reco::HitPairList &  hitPairList, reco::ClusterParametersList &  clusterParametersList, IHit3DBuilder::RecobHitToPtrMap &  hitToPtrMap  ) const

private

Top level output routine, allows checking cluster status.

Parameters

hitPairList	List of all 3D Hits in internal Cluster3D format
clusterParametersList	Data structure containing the cluster information to output
hitToPtrMap	This maps our Cluster2D hits back to art Ptr's to reco Hits

The workhorse to take the candidate 3D clusters and produce all of the necessary art output

Definition at line 1175 of file Cluster3D_module.cc.

  {
    /**
     *  @brief The workhorse to take the candidate 3D clusters and produce all of the necessary art output
     */

    mf::LogDebug("Cluster3D") << " *** Cluster3D::ProduceArtClusters() *** " << std::endl;

    // Make sure there is something to do here!
    if (!clusterParametersList.empty()) {
      // This is the loop over candidate 3D clusters
      // Note that it might be that the list of candidate clusters is modified by splitting
      // So we use the following construct to make sure we get all of them
      for (auto& clusterParameters : clusterParametersList) {
        // It should be straightforward at this point to transfer information from our vector of clusters
        // to the larsoft objects... of course we still have some work to do first, in particular to
        // find the candidate seeds and their seed hits

        // The chances of getting here and this condition not being true are probably zero... but check anyway
        if (!clusterParameters.getFullPCA().getSvdOK()) {
          mf::LogDebug("Cluster3D")
            << "--> no feature extraction done on this cluster!!" << std::endl;
          continue;
        }

        // Keep track of hit 3D to SP for when we do edges
        Hit3DToSPPtrMap hit3DToSPPtrMap;
        Hit3DToSPPtrMap best3DToSPPtrMap;

        // Do a special output of voronoi vertices here...
        dcel2d::VertexList& vertexList = clusterParameters.getVertexList();
        dcel2d::HalfEdgeList& halfEdgeList = clusterParameters.getHalfEdgeList();

        MakeAndSaveVertexPoints(output, vertexList, halfEdgeList);

        // Special case handling... if no daughters then call standard conversion routine to make sure space points
        // created, etc.
        if (clusterParameters.daughterList().empty()) {
          ConvertToArtOutput(gser,
                             output,
                             clusterParameters,
                             recob::PFParticle::kPFParticlePrimary,
                             hitToPtrMap,
                             hit3DToSPPtrMap,
                             best3DToSPPtrMap);

          // Get the extreme points
          MakeAndSaveKinkPoints(output,
                                clusterParameters.getConvexHull().getConvexHullKinkPoints());
        }
        // Otherwise, the cluster has daughters so we handle specially
        else {
          // Set up to keep track of parent/daughters
          IdxToPCAMap idxToPCAMap;
          size_t numTotalDaughters = countUltimateDaughters(clusterParameters);
          size_t pfParticleIdx(output.artPFParticleVector->size() + numTotalDaughters);

          FindAndStoreDaughters(gser,
                                output,
                                clusterParameters,
                                pfParticleIdx,
                                idxToPCAMap,
                                hitToPtrMap,
                                hit3DToSPPtrMap,
                                best3DToSPPtrMap);

          // Need to make a daughter vec from our map
          std::vector<size_t> daughterVec;

          for (auto& idxToPCA : idxToPCAMap)
            daughterVec.emplace_back(idxToPCA.first);

          // Now create/handle the parent PFParticle
          recob::PFParticle pfParticle(
            13, pfParticleIdx, recob::PFParticle::kPFParticlePrimary, daughterVec);
          output.artPFParticleVector->push_back(pfParticle);

          recob::PCAxis::EigenVectors eigenVecs;
          double eigenVals[] = {0., 0., 0.};
          double avePosition[] = {0., 0., 0.};

          eigenVecs.resize(3);

          reco::PrincipalComponents& skeletonPCA = clusterParameters.getSkeletonPCA();

          for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
            avePosition[outerIdx] = skeletonPCA.getAvePosition()[outerIdx];
            eigenVals[outerIdx] = skeletonPCA.getEigenValues()[outerIdx];

            eigenVecs[outerIdx].resize(3);

            // Be careful here... eigen stores in column major order buy default
            for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
              eigenVecs[outerIdx][innerIdx] = skeletonPCA.getEigenVectors().row(outerIdx)[innerIdx];
          }

          recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
                                   skeletonPCA.getNumHitsUsed(),
                                   eigenVals,
                                   eigenVecs,
                                   avePosition,
                                   skeletonPCA.getAveHitDoca(),
                                   output.artPCAxisVector->size());

          output.artPCAxisVector->push_back(skelPcAxis);

          reco::PrincipalComponents& fullPCA = clusterParameters.getFullPCA();

          for (size_t outerIdx = 0; outerIdx < 3; outerIdx++) {
            avePosition[outerIdx] = fullPCA.getAvePosition()[outerIdx];
            eigenVals[outerIdx] = fullPCA.getEigenValues()[outerIdx];

            for (size_t innerIdx = 0; innerIdx < 3; innerIdx++)
              eigenVecs[outerIdx][innerIdx] = fullPCA.getEigenVectors().row(outerIdx)(innerIdx);
          }

          recob::PCAxis fullPcAxis(fullPCA.getSvdOK(),
                                   fullPCA.getNumHitsUsed(),
                                   eigenVals,
                                   eigenVecs,
                                   avePosition,
                                   fullPCA.getAveHitDoca(),
                                   output.artPCAxisVector->size());

          output.artPCAxisVector->push_back(fullPcAxis);

          // Create associations to the PFParticle
          output.makePFPartPCAAssns();

          // Make associations to all space points for this cluster
          MakeAndSaveSpacePoints(output,
                                 *output.artSpacePointVector.get(),
                                 *output.artSPHitAssociations.get(),
                                 clusterParameters.getHitPairListPtr(),
                                 hitToPtrMap,
                                 hit3DToSPPtrMap);

          // Get the extreme points
          MakeAndSaveKinkPoints(output,
                                clusterParameters.getConvexHull().getConvexHullKinkPoints());

          // Build the edges now
          size_t edgeStart(output.artEdgeVector->size());

          for (const auto& edge : clusterParameters.getConvexHull().getConvexHullEdgeList()) {
            RecobSpacePointVector spacePointVec;

            try {
              spacePointVec.push_back(hit3DToSPPtrMap.at(std::get<0>(edge)));
              spacePointVec.push_back(hit3DToSPPtrMap.at(std::get<1>(edge)));
            }
            catch (...) {
              std::cout << "Caught exception in looking up space point ptr... " << std::get<0>(edge)
                        << ", " << std::get<1>(edge) << std::endl;
              continue;
            }

            output.artEdgeVector->emplace_back(std::get<2>(edge),
                                               spacePointVec[0].key(),
                                               spacePointVec[1].key(),
                                               output.artEdgeVector->size());

            output.makeEdgeSpacePointAssns(
              *output.artEdgeVector.get(), spacePointVec, *output.artEdgeSPAssociations.get());
          }

          output.makePFPartEdgeAssns(
            *output.artEdgeVector.get(), *output.artPFPartEdgeAssociations.get(), edgeStart);
        }
      }
    }

    // Right now error matrix is uniform...
    int nFreePoints(0);

    // Run through the HitPairVector and add any unused hit pairs to the list
    for (auto& hitPair : hitPairVector) {
      if (hitPair.bitsAreSet(reco::ClusterHit3D::MADESPACEPOINT)) continue;

      double spacePointPos[] = {
        hitPair.getPosition()[0], hitPair.getPosition()[1], hitPair.getPosition()[2]};
      double spacePointErr[] = {1., 0., 0., 1., 0., 1.};
      double chisq = hitPair.getHitChiSquare();

      RecobHitVector recobHits;

      for (const auto hit : hitPair.getHits()) {
        if (!hit) {
          chisq = -1000.;
          continue;
        }

        art::Ptr<recob::Hit> hitPtr = hitToPtrMap[hit->getHit()];
        recobHits.push_back(hitPtr);
      }

      nFreePoints++;

      spacePointErr[1] = hitPair.getTotalCharge();
      spacePointErr[3] = hitPair.getChargeAsymmetry();

      output.artSpacePointVector->push_back(
        recob::SpacePoint(spacePointPos, spacePointErr, chisq, output.artSpacePointVector->size()));

      if (!recobHits.empty())
        output.makeSpacePointHitAssns(
          *output.artSpacePointVector.get(), recobHits, *output.artSPHitAssociations.get());
    }

    std::cout << "++++>>>> total num hits: " << hitPairVector.size()
              << ", num free: " << nFreePoints << std::endl;
  }

void lar_cluster3d::Cluster3D::splitClustersWithHough ( reco::ClusterParameters &  clusterParameters, reco::ClusterParametersList &  clusterParametersList  ) const

private

Attempt to split clusters using the output of the Hough Filter.

Parameters

clusterParameters	The given cluster parameters object to try to split
clusterParametersList	The list of clusters

Definition at line 1031 of file Cluster3D_module.cc.

  {
    // @brief A method for splitted "crossed tracks" clusters into separate clusters
    //
    // If this routine is called then we believe we have a cluster which needs splitting.
    // The specific topology we are looking for is two long straight tracks which cross at some
    // point in close proximity so their hits were joined into a single 3D cluster. The method
    // to split this topology is to let the hough transform algorithm find the two leading candidates
    // and then to see if we use those to build two clusters instead of one.

    // Recover the hits we'll work on.
    // Note that we use on the skeleton hits so will need to recover them
    reco::HitPairListPtr& hitPairListPtr = clusterParameters.getHitPairListPtr();
    reco::HitPairListPtr skeletonListPtr;

    // We want to work with the "skeleton" hits so first step is to call the algorithm to
    // recover only these hits from the entire input collection
    m_skeletonAlg.GetSkeletonHits(hitPairListPtr, skeletonListPtr);

    // Skeleton hits are nice but we can do better if we then make a pass through to "average"
    // the skeleton hits position in the Y-Z plane
    m_skeletonAlg.AverageSkeletonPositions(skeletonListPtr);

    // Define the container for our lists of hits
    reco::HitPairListPtrList hitPairListPtrList;

    // Now feed this to the Hough Transform to find candidate straight lines
    m_seedFinderAlg.findTrackHits(
      skeletonListPtr, clusterParameters.getSkeletonPCA(), hitPairListPtrList);

    // We need at least two lists or else there is nothing to do
    if (hitPairListPtrList.size() < 2) return;

    // The game plan will be the following:
    // 1) Take the first list of hits and run the PCA on this to get an axis
    //    - Then calculate the 3d doca for ALL hits in the cluster to this axis
    //    - Move all hits within "3 sigam" of the axis to a new list
    // 2) run the PCA on the second list of hits to get that axis
    //    - Then calculate the 3d doca for all hits in our first list
    //    - Copy hits in the first list which are within 3 sigma of the new axis
    //      back into our original cluster - these are shared hits
    reco::HitPairListPtrList::iterator hitPairListIter = hitPairListPtrList.begin();
    reco::HitPairListPtr& firstHitList = *hitPairListIter++;
    reco::PrincipalComponents firstHitListPCA;

    m_pcaAlg.PCAAnalysis_3D(firstHitList, firstHitListPCA);

    // Make sure we have a successful calculation.
    if (firstHitListPCA.getSvdOK()) {
      // The fill routines below will expect to see unused 2D hits so we need to clear the
      // status bits... and I am not sure of a better way...
      for (const auto& hit3D : hitPairListPtr) {
        for (const auto& hit2D : hit3D->getHits())
          if (hit2D) hit2D->clearStatusBits(0x1);
      }

      // Calculate the 3D doca's for the hits which were used to make this PCA
      m_pcaAlg.PCAAnalysis_calc3DDocas(firstHitList, firstHitListPCA);

      // Divine from the ether some maximum allowed range for transfering hits
      float allowedHitRange = 6. * firstHitListPCA.getAveHitDoca();

      // Now go through and calculate the 3D doca's for ALL the hits in the original cluster
      m_pcaAlg.PCAAnalysis_calc3DDocas(hitPairListPtr, firstHitListPCA);

      // Let's make a new cluster to hold the hits
      clusterParametersList.push_back(reco::ClusterParameters());

      // Can we get a reference to what we just created?
      reco::ClusterParameters& newClusterParams = clusterParametersList.back();

      reco::HitPairListPtr& newClusterHitList = newClusterParams.getHitPairListPtr();

      newClusterHitList.resize(hitPairListPtr.size());

      // Do the actual copy of the hits we want
      reco::HitPairListPtr::iterator newListEnd = std::copy_if(hitPairListPtr.begin(),
                                                               hitPairListPtr.end(),
                                                               newClusterHitList.begin(),
                                                               CopyIfInRange(allowedHitRange));

      // Shrink to fit
      newClusterHitList.resize(std::distance(newClusterHitList.begin(), newListEnd));

      // And now remove these hits from the original cluster
      hitPairListPtr.remove_if(CopyIfInRange(allowedHitRange));

      // Get an empty hit to cluster map...
      reco::Hit2DToClusterMap hit2DToClusterMap;

      // Now "fill" the cluster parameters but turn off the hit rejection
      m_clusterBuilder->FillClusterParams(newClusterParams, hit2DToClusterMap, 0., 1.);

      // Set the skeleton pca to the value calculated above on input
      clusterParameters.getSkeletonPCA() = firstHitListPCA;

      // We are done with splitting out one track. Because the two tracks cross in
      // close proximity, this is the one case where we might consider sharing 3D hits
      // So let's make a little detour here to try to copy some of those hits back into
      // the main hit list
      reco::HitPairListPtr& secondHitList = *hitPairListIter;
      reco::PrincipalComponents secondHitListPCA;

      m_pcaAlg.PCAAnalysis_3D(secondHitList, secondHitListPCA);

      // Make sure we have a successful calculation.
      if (secondHitListPCA.getSvdOK()) {
        // Calculate the 3D doca's for the hits which were used to make this PCA
        m_pcaAlg.PCAAnalysis_calc3DDocas(secondHitList, secondHitListPCA);

        // Since this is the "other" cluster, we'll be a bit more generous in adding back hits
        float newAllowedHitRange = 6. * secondHitListPCA.getAveHitDoca();

        // Go through and calculate the 3D doca's for the hits in our new candidate cluster
        m_pcaAlg.PCAAnalysis_calc3DDocas(newClusterHitList, secondHitListPCA);

        // Create a temporary list to fill with the hits we might want to save
        reco::HitPairListPtr tempHitList(newClusterHitList.size());

        // Do the actual copy of the hits we want...
        reco::HitPairListPtr::iterator tempListEnd =
          std::copy_if(newClusterHitList.begin(),
                       newClusterHitList.end(),
                       tempHitList.begin(),
                       CopyIfInRange(newAllowedHitRange));

        hitPairListPtr.insert(hitPairListPtr.end(), tempHitList.begin(), tempListEnd);
      }

      // Of course, now we need to modify the original cluster parameters
      reco::ClusterParameters originalParams(hitPairListPtr);

      // Now "fill" the cluster parameters but turn off the hit rejection
      m_clusterBuilder->FillClusterParams(originalParams, hit2DToClusterMap, 0., 1.);

      // Overwrite original cluster parameters with our new values
      clusterParameters.getClusterParams() = originalParams.getClusterParams();
      clusterParameters.getFullPCA() = originalParams.getFullPCA();
      clusterParameters.getSkeletonPCA() = secondHitListPCA;
    }
  }

void lar_cluster3d::Cluster3D::splitClustersWithMST ( reco::ClusterParameters &  clusterParameters, reco::ClusterParametersList &  clusterParametersList  ) const

private

Attempt to split clusters by using a minimum spanning tree.

Parameters

clusterParameters	The given cluster parameters object to try to split
clusterParametersList	The list of clusters

Definition at line 838 of file Cluster3D_module.cc.

  {
    // This is being left in place for future development. Essentially, it was an attempt to implement
    // a Minimum Spanning Tree as a way to split a particular cluster topology, one where two straight
    // tracks cross closely enought to appear as one cluster. As of Feb 2, 2015 I think the idea is still
    // worth merit so am leaving this module in place for now.
    //
    // If this routine is called then we believe we have a cluster which needs splitting.
    // The way we will do this is to use a Minimum Spanning Tree algorithm to associate all
    // hits together by their distance apart. In theory, we should be able to split the cluster
    // by finding the largest distance and splitting at that point.
    //
    // Typedef some data structures that we will use.
    // Start with the adjacency map
    typedef std::pair<float, const reco::ClusterHit3D*> DistanceHit3DPair;
    typedef std::list<DistanceHit3DPair> DistanceHit3DPairList;
    typedef std::map<const reco::ClusterHit3D*, DistanceHit3DPairList> Hit3DToDistanceMap;

    // Now typedef the lists we'll keep
    typedef std::list<const reco::ClusterHit3D*> Hit3DList;
    typedef std::pair<Hit3DList::iterator, Hit3DList::iterator> Hit3DEdgePair;
    typedef std::pair<float, Hit3DEdgePair> DistanceEdgePair;
    typedef std::list<DistanceEdgePair> DistanceEdgePairList;

    struct DistanceEdgePairOrder {
      bool
      operator()(const DistanceEdgePair& left, const DistanceEdgePair& right) const
      {
        return left.first > right.first;
      }
    };

    // Recover the hits we'll work on.
    // Note that we use on the skeleton hits so will need to recover them
    reco::HitPairListPtr& hitPairListPtr = clusterParameters.getHitPairListPtr();
    reco::HitPairListPtr skeletonListPtr;

    // We want to work with the "skeleton" hits so first step is to call the algorithm to
    // recover only these hits from the entire input collection
    m_skeletonAlg.GetSkeletonHits(hitPairListPtr, skeletonListPtr);

    // Skeleton hits are nice but we can do better if we then make a pass through to "average"
    // the skeleton hits position in the Y-Z plane
    m_skeletonAlg.AverageSkeletonPositions(skeletonListPtr);

    // First task is to define and build the adjacency map
    Hit3DToDistanceMap hit3DToDistanceMap;

    for (reco::HitPairListPtr::const_iterator hit3DOuterItr = skeletonListPtr.begin();
         hit3DOuterItr != skeletonListPtr.end();) {
      const reco::ClusterHit3D* hit3DOuter = *hit3DOuterItr++;
      DistanceHit3DPairList& outerHitList = hit3DToDistanceMap[hit3DOuter];
      TVector3 outerPos(
        hit3DOuter->getPosition()[0], hit3DOuter->getPosition()[1], hit3DOuter->getPosition()[2]);

      for (reco::HitPairListPtr::const_iterator hit3DInnerItr = hit3DOuterItr;
           hit3DInnerItr != skeletonListPtr.end();
           hit3DInnerItr++) {
        const reco::ClusterHit3D* hit3DInner = *hit3DInnerItr;
        TVector3 innerPos(
          hit3DInner->getPosition()[0], hit3DInner->getPosition()[1], hit3DInner->getPosition()[2]);
        TVector3 deltaPos = innerPos - outerPos;
        float hitDistance(float(deltaPos.Mag()));

        if (hitDistance > 20.) continue;

        hit3DToDistanceMap[hit3DInner].emplace_back(hitDistance, hit3DOuter);
        outerHitList.emplace_back(hitDistance, hit3DInner);
      }

      // Make sure our membership bit is clear
      hit3DOuter->clearStatusBits(reco::ClusterHit3D::SELECTEDBYMST);
    }

    // Make pass through again to order each of the lists
    for (auto& mapPair : hit3DToDistanceMap) {
      mapPair.second.sort(Hit3DDistanceOrder());
    }

    // Get the containers for the MST to operate on/with
    Hit3DList hit3DList;
    DistanceEdgePairList distanceEdgePairList;

    // Initialize with first element
    hit3DList.emplace_back(skeletonListPtr.front());
    distanceEdgePairList.emplace_back(
      DistanceEdgePair(0., Hit3DEdgePair(hit3DList.begin(), hit3DList.begin())));

    skeletonListPtr.front()->setStatusBit(reco::ClusterHit3D::SELECTEDBYMST);

    float largestDistance(0.);
    float averageDistance(0.);

    // Now run the MST
    // Basically, we loop until the MST list is the same size as the input list
    while (hit3DList.size() < skeletonListPtr.size()) {
      Hit3DList::iterator bestHit3DIter = hit3DList.begin();
      float bestDist = 10000000.;

      // Loop through all hits currently in the list and look for closest hit not in the list
      for (Hit3DList::iterator hit3DIter = hit3DList.begin(); hit3DIter != hit3DList.end();
           hit3DIter++) {
        const reco::ClusterHit3D* hit3D = *hit3DIter;

        // For the given 3D hit, find the closest to it that is not already in the list
        DistanceHit3DPairList& nearestList = hit3DToDistanceMap[hit3D];

        while (!nearestList.empty()) {
          const reco::ClusterHit3D* hit3DToCheck = nearestList.front().second;

          if (!hit3DToCheck->bitsAreSet(reco::ClusterHit3D::SELECTEDBYMST)) {
            if (nearestList.front().first < bestDist) {
              bestHit3DIter = hit3DIter;
              bestDist = nearestList.front().first;
            }
            break;
          }

          nearestList.pop_front();
        }
      }

      if (bestDist > largestDistance) largestDistance = bestDist;

      averageDistance += bestDist;

      // Now we add the best hit not in the list to our list, keep track of the distance
      // to the object it was closest to
      const reco::ClusterHit3D* bestHit3D = *bestHit3DIter; // "best" hit already in the list
      const reco::ClusterHit3D* nextHit3D =
        hit3DToDistanceMap[bestHit3D].front().second; // "next" hit we are adding to the list

      Hit3DList::iterator nextHit3DIter = hit3DList.insert(hit3DList.end(), nextHit3D);

      distanceEdgePairList.emplace_back(bestDist, Hit3DEdgePair(bestHit3DIter, nextHit3DIter));

      nextHit3D->setStatusBit(reco::ClusterHit3D::SELECTEDBYMST);
    }

    averageDistance /= float(hit3DList.size());

    float thirdDist = 2. * sqrt(clusterParameters.getSkeletonPCA().getEigenValues()[2]);

    // Ok, find the largest distance in the iterator map
    distanceEdgePairList.sort(DistanceEdgePairOrder());

    DistanceEdgePairList::iterator largestDistIter = distanceEdgePairList.begin();

    for (DistanceEdgePairList::iterator edgeIter = distanceEdgePairList.begin();
         edgeIter != distanceEdgePairList.end();
         edgeIter++) {
      if (edgeIter->first < thirdDist) break;

      largestDistIter = edgeIter;
    }

    reco::HitPairListPtr::iterator breakIter = largestDistIter->second.second;
    reco::HitPairListPtr bestList;

    bestList.resize(std::distance(hit3DList.begin(), breakIter));

    std::copy(hit3DList.begin(), breakIter, bestList.begin());

    // Remove from the grand hit list and see what happens...
    // The pieces below are incomplete and were really for testing only.
    hitPairListPtr.sort();
    bestList.sort();

    reco::HitPairListPtr::iterator newListEnd = std::set_difference(hitPairListPtr.begin(),
                                                                    hitPairListPtr.end(),
                                                                    bestList.begin(),
                                                                    bestList.end(),
                                                                    hitPairListPtr.begin());

    hitPairListPtr.erase(newListEnd, hitPairListPtr.end());
  }

Member Data Documentation

float lar_cluster3d::Cluster3D::m_artHitsTime

private

Keeps track of time to recover hits.

Definition at line 499 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_buildNeighborhoodTime

private

Keeps track of time to build epsilon neighborhood.

Definition at line 501 of file Cluster3D_module.cc.

std::unique_ptr<lar_cluster3d::IClusterAlg> lar_cluster3d::Cluster3D::m_clusterAlg

private

Algorithm to do 3D space point clustering.

Definition at line 514 of file Cluster3D_module.cc.

std::unique_ptr<lar_cluster3d::IClusterParametersBuilder> lar_cluster3d::Cluster3D::m_clusterBuilder

private

Common cluster builder tool.

Definition at line 520 of file Cluster3D_module.cc.

std::unique_ptr<lar_cluster3d::IClusterModAlg> lar_cluster3d::Cluster3D::m_clusterMergeAlg

private

Algorithm to do cluster merging.

Definition at line 516 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_clusterMergeTime

private

Keeps track of the time to merge clusters.

Definition at line 503 of file Cluster3D_module.cc.

std::unique_ptr<lar_cluster3d::IClusterModAlg> lar_cluster3d::Cluster3D::m_clusterPathAlg

private

Algorithm to do cluster path finding.

Definition at line 518 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_dbscanTime

private

Keeps track of time to run DBScan.

Definition at line 502 of file Cluster3D_module.cc.

bool lar_cluster3d::Cluster3D::m_enableMonitoring

private

Turn on monitoring of this algorithm.

Definition at line 486 of file Cluster3D_module.cc.

int lar_cluster3d::Cluster3D::m_event

private

Definition at line 495 of file Cluster3D_module.cc.

std::string lar_cluster3d::Cluster3D::m_extremeInstance

private

Instance name for the extreme points.

Definition at line 508 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_finishTime

private

Keeps track of time to run output module.

Definition at line 505 of file Cluster3D_module.cc.

std::unique_ptr<lar_cluster3d::IHit3DBuilder> lar_cluster3d::Cluster3D::m_hit3DBuilderAlg

private

Builds the 3D hits to operate on.

Definition at line 512 of file Cluster3D_module.cc.

int lar_cluster3d::Cluster3D::m_hits

private

Keeps track of the number of hits seen.

Definition at line 496 of file Cluster3D_module.cc.

int lar_cluster3d::Cluster3D::m_hits3D

private

Keeps track of the number of 3D hits made.

Definition at line 497 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_makeHitsTime

private

Keeps track of time to build 3D hits.

Definition at line 500 of file Cluster3D_module.cc.

bool lar_cluster3d::Cluster3D::m_onlyMakSpacePoints

private

If true we don't do the full cluster 3D processing.

Algorithm parameters

Definition at line 485 of file Cluster3D_module.cc.

ParallelHitsSeedFinderAlg lar_cluster3d::Cluster3D::m_parallelHitsAlg

private

Deal with parallel hits clusters.

Definition at line 525 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_parallelHitsCosAng

private

Cut for PCA 3rd axis angle to X axis.

Definition at line 487 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_parallelHitsTransWid

private

Cut on transverse width of cluster (PCA 2nd eigenvalue)

Definition at line 488 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_pathFindingTime

private

Keeps track of the path finding time.

Definition at line 504 of file Cluster3D_module.cc.

std::string lar_cluster3d::Cluster3D::m_pathInstance

private

Special instance for path points.

Definition at line 506 of file Cluster3D_module.cc.

PrincipalComponentsAlg lar_cluster3d::Cluster3D::m_pcaAlg

private

Principal Components algorithm.

Definition at line 521 of file Cluster3D_module.cc.

PCASeedFinderAlg lar_cluster3d::Cluster3D::m_pcaSeedFinderAlg

private

Use PCA axis to find seeds.

Definition at line 524 of file Cluster3D_module.cc.

TTree* lar_cluster3d::Cluster3D::m_pRecoTree

private

Tree variables for output

Definition at line 493 of file Cluster3D_module.cc.

int lar_cluster3d::Cluster3D::m_run

private

Definition at line 494 of file Cluster3D_module.cc.

HoughSeedFinderAlg lar_cluster3d::Cluster3D::m_seedFinderAlg

private

Seed finder.

Definition at line 523 of file Cluster3D_module.cc.

SkeletonAlg lar_cluster3d::Cluster3D::m_skeletonAlg

private

Skeleton point finder.

Definition at line 522 of file Cluster3D_module.cc.

float lar_cluster3d::Cluster3D::m_totalTime

private

Keeps track of total execution time.

Definition at line 498 of file Cluster3D_module.cc.

std::string lar_cluster3d::Cluster3D::m_vertexInstance

private

Special instance name for vertex points.

Definition at line 507 of file Cluster3D_module.cc.

---

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

• Cluster3D_module.cc