Definition of the Cluster3D class. More...

Inheritance diagram for lar_cluster3d::Cluster3D:

Classes
class	ArtOutputHandler

Public Types
using	ModuleType = EDProducer

template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

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

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

void	fillProductDescriptions ()

void	registerProducts (ProductDescriptions &productsToRegister)

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

std::unique_ptr< Worker >	makeWorker (WorkerParams const &wp)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Protected Member Functions
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

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	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 110 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 390 of file Cluster3D_module.cc.

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 17 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::detail::Producer::Table = Modifier::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 26 of file Producer.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 517 of file Cluster3D_module.cc.

References m_clusterAlg, m_clusterBuilder, m_clusterMergeAlg, m_clusterPathAlg, m_enableMonitoring, m_extremeInstance, m_hit3DBuilderAlg, m_onlyMakSpacePoints, m_parallelHitsAlg, m_parallelHitsCosAng, m_parallelHitsTransWid, m_pathInstance, m_pcaAlg, m_pcaSeedFinderAlg, m_seedFinderAlg, m_skeletonAlg, m_vertexInstance, and art::ProductRegistryHelper::producesCollector().

     : 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 449 of file Cluster3D_module.cc.

References countUltimateDaughters(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), m_parallelHitsCosAng, and m_parallelHitsTransWid.

Referenced by findTrackSeeds().

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

void lar_cluster3d::Cluster3D::beginJob ( )

overrideprivatevirtual

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?)

Reimplemented from art::EDProducer.

Definition at line 580 of file Cluster3D_module.cc.

References InitializeMonitoring(), and m_enableMonitoring.

   {
     if (m_enableMonitoring) this->InitializeMonitoring();
   }

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::consumes ( InputTag const & tag )

protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

   {
     return collector_.consumes<T, BT>(tag);
   }

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

   {
     collector_.consumesMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::consumesView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const & tag )

inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

   {
     return collector_.consumesView<T, BT>(tag);
   }

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 1238 of file Cluster3D_module.cc.

Referenced by FindAndStoreDaughters(), and ProduceArtClusters().

   {
     // 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 1187 of file Cluster3D_module.cc.

References reco::ClusterParameters::daughterList().

Referenced by aParallelHitsCluster(), and ProduceArtClusters().

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

void art::detail::Producer::doBeginJob ( SharedResources const & resources )

inherited

Definition at line 22 of file Producer.cc.

References art::detail::Producer::beginJobWithFrame(), and art::detail::Producer::setupQueues().

   {
     setupQueues(resources);
     ProcessingFrame const frame{ScheduleID{}};
     beginJobWithFrame(frame);
   }

bool art::detail::Producer::doBeginRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 65 of file Producer.cc.

References art::detail::Producer::beginRunWithFrame(), art::RangeSet::forRun(), art::RunPrincipal::makeRun(), r, art::RunPrincipal::runID(), and art::ModuleContext::scheduleID().

   {
     auto r = rp.makeRun(mc, RangeSet::forRun(rp.runID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doBeginSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 85 of file Producer.cc.

References art::detail::Producer::beginSubRunWithFrame(), art::RangeSet::forSubRun(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::SubRunPrincipal::subRunID().

   {
     auto sr = srp.makeSubRun(mc, RangeSet::forSubRun(srp.subRunID()));
     ProcessingFrame const frame{mc.scheduleID()};
     beginSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

void art::detail::Producer::doEndJob ( )

inherited

Definition at line 30 of file Producer.cc.

References art::detail::Producer::endJobWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     endJobWithFrame(frame);
   }

bool art::detail::Producer::doEndRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 75 of file Producer.cc.

References art::detail::Producer::endRunWithFrame(), art::RunPrincipal::makeRun(), r, art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto r = rp.makeRun(mc, rp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endRunWithFrame(r, frame);
     r.commitProducts();
     return true;
   }

bool art::detail::Producer::doEndSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 95 of file Producer.cc.

References art::detail::Producer::endSubRunWithFrame(), art::SubRunPrincipal::makeSubRun(), art::ModuleContext::scheduleID(), and art::Principal::seenRanges().

   {
     auto sr = srp.makeSubRun(mc, srp.seenRanges());
     ProcessingFrame const frame{mc.scheduleID()};
     endSubRunWithFrame(sr, frame);
     sr.commitProducts();
     return true;
   }

bool art::detail::Producer::doEvent	(	EventPrincipal &	ep,
		ModuleContext const &	mc,
		std::atomic< std::size_t > &	counts_run,
		std::atomic< std::size_t > &	counts_passed,
		std::atomic< std::size_t > &	counts_failed
	)

inherited

Definition at line 105 of file Producer.cc.

References art::detail::Producer::checkPutProducts_, e, art::EventPrincipal::makeEvent(), art::detail::Producer::produceWithFrame(), and art::ModuleContext::scheduleID().

   {
     auto e = ep.makeEvent(mc);
     ++counts_run;
     ProcessingFrame const frame{mc.scheduleID()};
     produceWithFrame(e, frame);
     e.commitProducts(checkPutProducts_, &expectedProducts<InEvent>());
     ++counts_passed;
     return true;
   }

void art::detail::Producer::doRespondToCloseInputFile ( FileBlock const & fb )

inherited

Definition at line 44 of file Producer.cc.

References art::detail::Producer::respondToCloseInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToCloseOutputFiles ( FileBlock const & fb )

inherited

Definition at line 58 of file Producer.cc.

References art::detail::Producer::respondToCloseOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseOutputFilesWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenInputFile ( FileBlock const & fb )

inherited

Definition at line 37 of file Producer.cc.

References art::detail::Producer::respondToOpenInputFileWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenInputFileWithFrame(fb, frame);
   }

void art::detail::Producer::doRespondToOpenOutputFiles ( FileBlock const & fb )

inherited

Definition at line 51 of file Producer.cc.

References art::detail::Producer::respondToOpenOutputFilesWithFrame().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenOutputFilesWithFrame(fb, frame);
   }

void art::Modifier::fillProductDescriptions ( )

inherited

Definition at line 10 of file Modifier.cc.

References art::ProductRegistryHelper::fillDescriptions(), and art::ModuleBase::moduleDescription().

   {
     ProductRegistryHelper::fillDescriptions(moduleDescription());
   }

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 1201 of file Cluster3D_module.cc.

References ConvertToArtOutput(), reco::ClusterParameters::daughterList(), and reco::ClusterParameters::getFullPCA().

Referenced by ProduceArtClusters().

   {
     // 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 721 of file Cluster3D_module.cc.

References aParallelHitsCluster(), lar_cluster3d::SkeletonAlg::AverageSkeletonPositions(), util::CreateAssn(), lar_cluster3d::ParallelHitsSeedFinderAlg::findTrackSeeds(), lar_cluster3d::PCASeedFinderAlg::findTrackSeeds(), lar_cluster3d::HoughSeedFinderAlg::findTrackSeeds(), reco::PrincipalComponents::getEigenValues(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHitPairListPtr(), lar_cluster3d::SkeletonAlg::GetSkeletonHits(), reco::ClusterParameters::getSkeletonPCA(), m_parallelHitsAlg, m_pcaSeedFinderAlg, m_seedFinderAlg, m_skeletonAlg, and art::PtrVector< T >::push_back().

   {
     // 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);
     }
   }

std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const

inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

   {
     return collector_.getConsumables();
   }

void lar_cluster3d::Cluster3D::InitializeMonitoring ( )

private

Initialize the internal monitoring.

Definition at line 682 of file Cluster3D_module.cc.

References art::ServiceHandle< T, SCOPE >::get(), m_artHitsTime, m_buildNeighborhoodTime, m_clusterMergeTime, m_clusterPathAlg, m_dbscanTime, m_event, m_finishTime, m_hits, m_hits3D, m_makeHitsTime, m_pathFindingTime, m_pRecoTree, m_run, and m_totalTime.

Referenced by beginJob().

   {
     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 1516 of file Cluster3D_module.cc.

References lar_cluster3d::Cluster3D::ArtOutputHandler::artExtremePointVector, reco::ClusterHit3D::getChargeAsymmetry(), reco::ClusterHit3D::getHitChiSquare(), reco::ClusterHit3D::getPosition(), and reco::ClusterHit3D::getTotalCharge().

Referenced by ProduceArtClusters().

   {
     // 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 1607 of file Cluster3D_module.cc.

References lar_cluster3d::Cluster3D::ArtOutputHandler::artEdgeVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artVertexEdgeVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artVertexPointVector, reco::PrincipalComponents::getAvePosition(), reco::PrincipalComponents::getEigenVectors(), art::left(), and art::right().

   {
     // 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 1459 of file Cluster3D_module.cc.

References art::PtrVector< T >::empty(), reco::ClusterHit3D::MADESPACEPOINT, lar_cluster3d::Cluster3D::ArtOutputHandler::makeSpacePointHitAssns(), lar_cluster3d::Cluster3D::ArtOutputHandler::makeSpacePointPtr(), art::PtrVector< T >::push_back(), and reco::ClusterHit3D::REJECTEDHIT.

Referenced by ConvertToArtOutput(), and ProduceArtClusters().

   {
     // 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 1542 of file Cluster3D_module.cc.

References lar_cluster3d::Cluster3D::ArtOutputHandler::artEdgeVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artVertexEdgeVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artVertexPointVector, dcel2d::Vertex::getCoords(), dcel2d::HalfEdge::getTargetVertex(), and dcel2d::HalfEdge::getTwinHalfEdge().

Referenced by ProduceArtClusters().

   {
     // 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;
   }

std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const & wp )

inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

   {
     return doMakeWorker(wp);
   }

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::mayConsume ( InputTag const & tag )

protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

   {
     return collector_.mayConsume<T, BT>(tag);
   }

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

   {
     collector_.mayConsumeMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const & tag )

inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

   {
     return collector_.mayConsumeView<T, BT>(tag);
   }

ModuleDescription const & art::ModuleBase::moduleDescription ( ) const

inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

   {
     if (md_.has_value()) {
       return *md_;
     }
 
     throw Exception{errors::LogicError,
                     "There was an error while calling moduleDescription().\n"}
       << "The moduleDescription() base-class member function cannot be called\n"
          "during module construction.  To determine which module is "
          "responsible\n"
          "for calling it, find the '<module type>:<module "
          "label>@Construction'\n"
          "tag in the message prefix above.  Please contact artists@fnal.gov\n"
          "for guidance.\n";
   }

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

private

Event Preparation.

Parameters

evt	the ART event

Definition at line 704 of file Cluster3D_module.cc.

References art::EventID::event(), art::Event::id(), m_artHitsTime, m_buildNeighborhoodTime, m_dbscanTime, m_event, m_finishTime, m_hits, m_hits3D, m_makeHitsTime, m_pathFindingTime, m_run, m_totalTime, and art::Event::run().

Referenced by produce().

   {
     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 )

overrideprivatevirtual

Implements art::EDProducer.

Definition at line 591 of file Cluster3D_module.cc.

References lar_cluster3d::IClusterAlg::BUILDCLUSTERINFO, lar_cluster3d::IClusterAlg::BUILDHITTOHITMAP, lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS, lar_cluster3d::IHit3DBuilder::COLLECTARTHITS, art::EventID::event(), art::Event::id(), m_artHitsTime, m_buildNeighborhoodTime, m_clusterAlg, m_clusterMergeAlg, m_clusterMergeTime, m_clusterPathAlg, m_dbscanTime, m_enableMonitoring, m_event, m_extremeInstance, m_finishTime, m_hit3DBuilderAlg, m_hits, m_hits3D, m_makeHitsTime, m_onlyMakSpacePoints, m_pathFindingTime, m_pathInstance, m_pRecoTree, m_run, m_totalTime, m_vertexInstance, lar_cluster3d::Cluster3D::ArtOutputHandler::outputObjects(), PrepareEvent(), ProduceArtClusters(), art::Event::run(), and lar_cluster3d::IClusterAlg::RUNDBSCAN.

   {
     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>(),
                                        art::ServiceHandle<geo::WireReadout const>()->Get(),
                                        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 970 of file Cluster3D_module.cc.

Referenced by produce().

   {
     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 art::Modifier::registerProducts ( ProductDescriptions & productsToRegister )

inherited

Definition at line 16 of file Modifier.cc.

References art::ModuleBase::moduleDescription(), and art::ProductRegistryHelper::registerProducts().

   {
     ProductRegistryHelper::registerProducts(productsToRegister,
                                             moduleDescription());
   }

void art::ModuleBase::setModuleDescription ( ModuleDescription const & md )

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

void art::ModuleBase::sortConsumables ( std::string const & current_process_name )

inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

   {
     // Now that we know we have seen all the consumes declarations,
     // sort the results for fast lookup later.
     collector_.sortConsumables(current_process_name);
   }

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 827 of file Cluster3D_module.cc.

References lar_cluster3d::SkeletonAlg::AverageSkeletonPositions(), lar_cluster3d::HoughSeedFinderAlg::findTrackHits(), reco::PrincipalComponents::getAveHitDoca(), reco::ClusterParameters::getClusterParams(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHitPairListPtr(), lar_cluster3d::SkeletonAlg::GetSkeletonHits(), reco::ClusterParameters::getSkeletonPCA(), reco::PrincipalComponents::getSvdOK(), m_clusterBuilder, m_pcaAlg, m_seedFinderAlg, m_skeletonAlg, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas(), and x1.

   {
     // @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;
     }
   }

Member Data Documentation

float lar_cluster3d::Cluster3D::m_artHitsTime

private

Keeps track of time to recover hits.

Definition at line 479 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

float lar_cluster3d::Cluster3D::m_buildNeighborhoodTime

private

Keeps track of time to build epsilon neighborhood.

Definition at line 481 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

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

private

Algorithm to do 3D space point clustering.

Definition at line 494 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

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

private

Common cluster builder tool.

Definition at line 500 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and splitClustersWithHough().

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

private

Algorithm to do cluster merging.

Definition at line 496 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

float lar_cluster3d::Cluster3D::m_clusterMergeTime

private

Keeps track of the time to merge clusters.

Definition at line 483 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), and produce().

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

private

Algorithm to do cluster path finding.

Definition at line 498 of file Cluster3D_module.cc.

Referenced by Cluster3D(), InitializeMonitoring(), and produce().

float lar_cluster3d::Cluster3D::m_dbscanTime

private

Keeps track of time to run DBScan.

Definition at line 482 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

bool lar_cluster3d::Cluster3D::m_enableMonitoring

private

Turn on monitoring of this algorithm.

Definition at line 466 of file Cluster3D_module.cc.

Referenced by beginJob(), Cluster3D(), and produce().

int lar_cluster3d::Cluster3D::m_event

private

Definition at line 475 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

std::string lar_cluster3d::Cluster3D::m_extremeInstance

private

Instance name for the extreme points.

Definition at line 488 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

float lar_cluster3d::Cluster3D::m_finishTime

private

Keeps track of time to run output module.

Definition at line 485 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

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

private

Builds the 3D hits to operate on.

Definition at line 492 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

int lar_cluster3d::Cluster3D::m_hits

private

Keeps track of the number of hits seen.

Definition at line 476 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

int lar_cluster3d::Cluster3D::m_hits3D

private

Keeps track of the number of 3D hits made.

Definition at line 477 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

float lar_cluster3d::Cluster3D::m_makeHitsTime

private

Keeps track of time to build 3D hits.

Definition at line 480 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

bool lar_cluster3d::Cluster3D::m_onlyMakSpacePoints

private

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

Algorithm parameters

Definition at line 465 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

ParallelHitsSeedFinderAlg lar_cluster3d::Cluster3D::m_parallelHitsAlg

private

Deal with parallel hits clusters.

Definition at line 505 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and findTrackSeeds().

float lar_cluster3d::Cluster3D::m_parallelHitsCosAng

private

Cut for PCA 3rd axis angle to X axis.

Definition at line 467 of file Cluster3D_module.cc.

Referenced by aParallelHitsCluster(), and Cluster3D().

float lar_cluster3d::Cluster3D::m_parallelHitsTransWid

private

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

Definition at line 468 of file Cluster3D_module.cc.

Referenced by aParallelHitsCluster(), and Cluster3D().

float lar_cluster3d::Cluster3D::m_pathFindingTime

private

Keeps track of the path finding time.

Definition at line 484 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

std::string lar_cluster3d::Cluster3D::m_pathInstance

private

Special instance for path points.

Definition at line 486 of file Cluster3D_module.cc.

Referenced by Cluster3D(), ConvertToArtOutput(), and produce().

PrincipalComponentsAlg lar_cluster3d::Cluster3D::m_pcaAlg

private

Principal Components algorithm.

Definition at line 501 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and splitClustersWithHough().

PCASeedFinderAlg lar_cluster3d::Cluster3D::m_pcaSeedFinderAlg

private

Use PCA axis to find seeds.

Definition at line 504 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and findTrackSeeds().

TTree* lar_cluster3d::Cluster3D::m_pRecoTree

private

Tree variables for output

Definition at line 473 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), and produce().

int lar_cluster3d::Cluster3D::m_run

private

Definition at line 474 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

HoughSeedFinderAlg lar_cluster3d::Cluster3D::m_seedFinderAlg

private

Seed finder.

Definition at line 503 of file Cluster3D_module.cc.

Referenced by Cluster3D(), findTrackSeeds(), and splitClustersWithHough().

SkeletonAlg lar_cluster3d::Cluster3D::m_skeletonAlg

private

Skeleton point finder.

Definition at line 502 of file Cluster3D_module.cc.

Referenced by Cluster3D(), findTrackSeeds(), and splitClustersWithHough().

float lar_cluster3d::Cluster3D::m_totalTime

private

Keeps track of total execution time.

Definition at line 478 of file Cluster3D_module.cc.

Referenced by InitializeMonitoring(), PrepareEvent(), and produce().

std::string lar_cluster3d::Cluster3D::m_vertexInstance

private

Special instance name for vertex points.

Definition at line 487 of file Cluster3D_module.cc.

Referenced by Cluster3D(), and produce().

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

larreco/v10_01_07/source/larreco/ClusterFinder/Cluster3D_module.cc

Classes

Public Types

Public Member Functions

Protected Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation