Definition of the Cluster3D class. More...

Inheritance diagram for lar_cluster3d::Cluster3D:

Classes
class	ArtOutputHandler

Public Types
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

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

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

virtual	~Cluster3D ()
	Destructor. More...

void	beginJob ()
	declare the standard art functions that we'll implement in this producer module More...

void	endJob ()

void	produce (art::Event &evt)

void	reconfigure (fhicl::ParameterSet const &pset)

template<typename PROD , BranchType B = InEvent>
ProductID	getProductID (std::string const &instanceName={}) const

template<typename PROD , BranchType B>
ProductID	getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const

bool	modifiesEvent () const

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

template<typename T , art::BranchType BT>
art::ProductToken< T >	consumes (InputTag const &it)

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

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

template<typename T , art::BranchType BT>
art::ViewToken< T >	consumesView (InputTag const &it)

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

template<typename T , art::BranchType BT>
art::ProductToken< T >	mayConsume (InputTag const &it)

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

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

template<typename T , art::BranchType BT>
art::ViewToken< T >	mayConsumeView (InputTag const &it)

base_engine_t &	createEngine (seed_t seed)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make)

base_engine_t &	createEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)

seed_t	get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)

Static Public Member Functions
static cet::exempt_ptr< Consumer >	non_module_context ()

Protected Member Functions
CurrentProcessingContext const *	currentContext () const

void	validateConsumedProduct (BranchType const bt, ProductInfo const &pi)

void	prepareForJob (fhicl::ParameterSet const &pset)

void	showMissingConsumes () const

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	PrepareEvent (const art::Event &evt)
	Event Preparation. More...

void	InitializeMonitoring ()
	Initialize the internal monitoring. More...

void	findTrackSeeds (art::Event &evt, reco::ClusterParameters &cluster, 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...

size_t	ConvertToArtOutput (ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap) const
	Produces the art output from all the work done in this producer module. More...

void	MakeAndSaveSpacePoints (ArtOutputHandler &output, reco::HitPairListPtr &clusHitPairVector, RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap, int spacePointStart) 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 (ArtOutputHandler &output, reco::ClusterParameters &clusterParameters, size_t pfParticleParent, IdxToPCAMap &idxToPCAMap, RecobHitToPtrMap &hitToPtrMap, Hit3DToSPPtrMap &hit3DToSPPtrMap) const
	This will produce art output for daughters noting that it needs to be done recursively. More...

void	ProduceArtClusters (ArtOutputHandler &output, reco::HitPairList &hitPairList, reco::ClusterParametersList &clusterParametersList, RecobHitToPtrMap &hitToPtrMap) const
	Top level output routine, allows checking cluster status. 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_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...

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_pathFindingTime
	Keeps track of the path finding time. More...

float	m_finishTime
	Keeps track of time to run output module. More...

std::string	m_spacePointInstance
	Special instance name for vertex points. More...

geo::Geometry *	m_geometry
	pointer to the Geometry service More...

const detinfo::DetectorProperties *	m_detector
	Pointer to the detector properties. 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...

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

using art::EDProducer::ModuleType = EDProducer

inherited

Definition at line 34 of file EDProducer.h.

template<typename UserConfig , typename KeysToIgnore = void>

using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>

inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>

inherited

Definition at line 35 of file EDProducer.h.

Constructor & Destructor Documentation

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

Constructor.

Parameters

pset	- reference to the parameters used by this module and its algorithms

Definition at line 433 of file Cluster3D_module.cc.

References m_spacePointInstance, and reconfigure().

                                                   :
     m_clusterBuilder(pset.get<fhicl::ParameterSet>("ClusterParamsBuilder")),
     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"))
 {
     this->reconfigure(pset);
     
     m_spacePointInstance = "Voronoi";
 
     produces< std::vector<recob::PCAxis>>();
     produces< std::vector<recob::PFParticle>>();
     produces< std::vector<recob::Cluster>>();
     produces< std::vector<recob::SpacePoint>>();
     produces< std::vector<recob::SpacePoint>>(m_spacePointInstance);
     produces< std::vector<recob::Seed>>();
     produces< std::vector<recob::Edge>>();
     produces< std::vector<recob::Edge>>(m_spacePointInstance);
     produces< art::Assns<recob::PFParticle, recob::PCAxis>>();
     produces< art::Assns<recob::PFParticle, recob::Cluster>>();
     produces< art::Assns<recob::PFParticle, recob::SpacePoint>>();
     produces< art::Assns<recob::PFParticle, recob::Seed>>();
     produces< art::Assns<recob::PFParticle, recob::Edge>>();
     produces< art::Assns<recob::Seed,       recob::Hit>>();
     produces< art::Assns<recob::Cluster,    recob::Hit>>();
     produces< art::Assns<recob::SpacePoint, recob::Hit>>();
     produces< art::Assns<recob::Edge,       recob::SpacePoint>>();
 }

lar_cluster3d::Cluster3D::~Cluster3D ( )

virtual

Destructor.

Definition at line 466 of file Cluster3D_module.cc.

467 {

468 }

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 370 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()[2][0]) > m_parallelHitsCosAng && 3. * sqrt(pca.getEigenValues()[1]) > m_parallelHitsTransWid;
     }

void lar_cluster3d::Cluster3D::beginJob ( )

virtual

declare the standard art functions that we'll implement in this producer module

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

References InitializeMonitoring(), m_detector, m_enableMonitoring, and m_geometry.

 {
     if (m_enableMonitoring)
         this->InitializeMonitoring();
     
     art::ServiceHandle<geo::Geometry> geometry;
     
     m_geometry = &*geometry;
     m_detector = lar::providerFrom<detinfo::DetectorPropertiesService>();
 }

template<typename T , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::consumes ( InputTag const & it )

inherited

Definition at line 147 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::consumesMany ( )

inherited

Definition at line 162 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::consumesView ( InputTag const & it )

inherited

Definition at line 172 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

size_t lar_cluster3d::Cluster3D::ConvertToArtOutput	(	ArtOutputHandler &	output,
		reco::ClusterParameters &	clusterParameters,
		size_t	pfParticleParent,
		RecobHitToPtrMap &	hitToPtrMap,
		Hit3DToSPPtrMap &	hit3DToSPPtrMap
	)		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 1289 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();
     
     // As tracks become more parallel to the wire plane the number of "ambiguous" 3D hits can increase
     // rapidly. Now that we have more information we can go back through these hits and do a better job
     // selecting "the right ones". Here we call the "medial skeleton" algorithm which uses a modification
     // of a standard medial skeleton procedure to get the 3D hits we want
     // But note that even this is hopeless in the worst case and, in fact, it can be a time waster
     // So bypass when you recognize that condition
     /*
      if (!aParallelHitsCluster(fullPCA))
      {
      int nSkeletonPoints = m_skeletonAlg.FindMedialSkeleton(clusterParameters.getHitPairListPtr());
      
      // If enough skeleton points then rerun pca with only those
      if (nSkeletonPoints > 10)
      {
      // Now rerun the principal components axis on just those points
      m_pcaAlg.PCAAnalysis_3D(clusterParameters.getHitPairListPtr(), skeletonPCA, true);
      
      // If there was a failure (can that happen?) then restore the full PCA
      if (!skeletonPCA.getSvdOK()) skeletonPCA = fullPCA;
      }
      
      // Here we can try to handle a specific case. It can happen that two tracks (think CR muons here) pass so
      // close together at some point to get merged into one cluster. Now that we have skeletonized the hits and
      // have run the PCA on the skeleton points we can try to divide these two tracks. The signature will be that
      // their are a large number of total hits, that the PCA will have a large spread in two dimensions. The
      // spread in the third dimension will be an indicator of the actual separation between the two tracks
      // which we might try to exploit in the actual algorithm.
      // hardwire for now to see what is going on...
      if (skeletonPCA.getNumHitsUsed() > 1000 && skeletonPCA.getEigenValues()[1] > 100. && fabs(skeletonPCA.getEigenVectors()[2][0]) < m_parallelHitsCosAng)
      {
      mf::LogDebug("Cluster3D") << "--> Detected crossed axes!! Total # hits: " << fullPCA.getNumHitsUsed() <<
      "\n    Skeleton PCA # hits: " << skeletonPCA.getNumHitsUsed() << ", eigenValues: " <<
      skeletonPCA.getEigenValues()[0] << ", " <<skeletonPCA.getEigenValues()[1] << ", " <<skeletonPCA.getEigenValues()[2] << std::endl;
      
      splitClustersWithHough(clusterParameters, clusterParametersList);
      }
      }
      */
     size_t clusterStart = output.artClusterVector->size();
     
     // Start loop over views to build out the hit lists and the 2D cluster objects
     for(reco::PlaneToClusterParamsMap::const_iterator planeItr = clusterParameters.getClusterParams().begin(); planeItr != clusterParameters.getClusterParams().end(); planeItr++)
     {
         const reco::RecobClusterParameters& clusParams = planeItr->second;
         
         // Protect against a missing view
         if (clusParams.m_view == geo::kUnknown) continue;
         
         // 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(reco::ClusterHit2DVec::const_iterator hitItr = clusParams.m_hitVector.begin(); hitItr != clusParams.m_hitVector.end(); hitItr++)
         {
             art::Ptr<recob::Hit> hitPtr = hitToPtrMap[&(*hitItr)->getHit()];
             recobHits.push_back(hitPtr);
         }
         
         // And sorting! Sorting is good for the mind, soul and body
         // ooopsss... don't do this else event display will look funky
         //                std::sort(recobHits.begin(), recobHits.end());
         
         // 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);
         
         // plane ID is not a part of clusParams... get the one from the first hit
         geo::PlaneID plane; // invalid by default
         if (!recobHits.empty())
             plane = recobHits.front()->WireID().planeID();
         
         // feed the algorithm with all the cluster hits
         ClusterParamAlgo.ImportHits(recobHits);
         
         // create the recob::Cluster directly in the vector
         cluster::ClusterCreator artCluster(
                                            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
                                            plane,                                // plane
                                            recob::Cluster::Sentry                // sentry
                                            );
         
         output.artClusterVector->emplace_back(artCluster.move());
         
         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();
     
     // Call the magical algorith to do the dirty work
     //            findTrackSeeds(evt, clusterParameters, hitToPtrMap, *artSeedVector, *artSeedHitAssociations);
 
     // Deal with converting the Hit Pairs to art
     // Recover the hit pairs and start looping! Love to loop!
     reco::HitPairListPtr& clusHitPairVector = clusterParameters.getHitPairListPtr();
     //            reco::HitPairListPtr& clusHitPairVector = clusterParameters.getBestHitPairListPtr();
     
     // Keep track of current start for space points
     int spacePointStart(output.artSpacePointVector->size());
     
     // Copy these hits to the vector to be stored with the event
     for (auto& hitPair : clusHitPairVector)
     {
         // Don't make space point if this hit was "rejected"
         if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
         
         double chisq = hitPair->getHitChiSquare();    // secret handshake...
         
 //        if      ( hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) && !hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -1.;  // pure skeleton point
 //        else if (!hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) &&  hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -2.;  // pure edge point
 //        else if ( hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) &&  hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -3.;  // skeleton and edge point
 //
 //        if      (hitPair->bitsAreSet(reco::ClusterHit3D::SEEDHIT)                                                          ) chisq = -4.;  // Seed point
 //
 //        if ((hitPair->getStatusBits() & 0x7) != 0x7) chisq = -10.;
         
         // Mark this hit pair as in use
         hitPair->setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);
 
         // Create and store the space point
         size_t spacePointID    = output.artSpacePointVector->size();
         double spacePointPos[] = {hitPair->getPosition()[0],hitPair->getPosition()[1],hitPair->getPosition()[2]};
         double spacePointErr[] = {m_detector->GetXTicksCoefficient() * hitPair->getSigmaPeakTime(), 0., 0., 0.15, 0., 0.15};
         output.artSpacePointVector->push_back(recob::SpacePoint(spacePointPos, spacePointErr, chisq, output.artSpacePointVector->size()));
         
         // Update mapping
         hit3DToSPPtrMap[hitPair] = spacePointID;
         
         // 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(recobHits);
     }
 
     // Build the edges now
     size_t edgeStart(output.artEdgeVector->size());
     
     for(const auto& edge : clusterParameters.getBestEdgeList())
         output.artEdgeVector->emplace_back(std::get<2>(edge), hit3DToSPPtrMap[std::get<0>(edge)], hit3DToSPPtrMap[std::get<1>(edge)], output.artEdgeVector->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);
     
     // 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()[outerIdx][innerIdx];
     }
     
     
     recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
                              skeletonPCA.getNumHitsUsed(),
                              eigenVals,                      //skeletonPCA.getEigenValues(),
                              eigenVecs,                      //skeletonPCA.getEigenVectors(),
                              avePosition,                    //skeletonPCA.getAvePosition(),
                              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()[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);
     output.makePFPartSpacePointAssns(spacePointStart);
     output.makePFPartEdgeAssns(edgeStart);
 
     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 1251 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;
 }

EngineCreator::base_engine_t & EngineCreator::createEngine ( seed_t seed )

inherited

Definition at line 26 of file EngineCreator.cc.

References art::EngineCreator::rng().

Referenced by evgen::CosmicsGen::CosmicsGen(), rndm::NuRandomService::createEngine(), cluster::fuzzyCluster::fuzzyCluster(), cluster::HoughLineFinder::HoughLineFinder(), art::MixFilter< T >::initEngine_(), larg4::LArG4::LArG4(), evgen::LightSource::LightSource(), evgen::NeutronOsc::NeutronOsc(), evgen::NucleonDecay::NucleonDecay(), opdet::OpMCDigi::OpMCDigi(), opdet::OptDetDigitizer::OptDetDigitizer(), phot::PhotonLibraryPropagation::PhotonLibraryPropagation(), detsim::SimDriftElectrons::SimDriftElectrons(), evgen::SingleGen::SingleGen(), evgen::SNNueAr40CCGen::SNNueAr40CCGen(), ToyOneShowerGen::ToyOneShowerGen(), and trkf::Track3DKalman::Track3DKalman().

 {
   return rng()->createEngine(placeholder_schedule_id(), seed);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make
	)

inherited

Definition at line 32 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make);
 }

EngineCreator::base_engine_t & EngineCreator::createEngine	(	seed_t	seed,
		std::string const &	kind_of_engine_to_make,
		label_t const &	engine_label
	)

inherited

Definition at line 40 of file EngineCreator.cc.

References art::EngineCreator::rng().

 {
   return rng()->createEngine(
     placeholder_schedule_id(), seed, kind_of_engine_to_make, engine_label);
 }

CurrentProcessingContext const * art::EDProducer::currentContext ( ) const

protectedinherited

Definition at line 120 of file EDProducer.cc.

References art::EDProducer::current_context_.

   {
     return current_context_.get();
   }

void lar_cluster3d::Cluster3D::endJob ( )

virtual

Reimplemented from art::EDProducer.

Definition at line 509 of file Cluster3D_module.cc.

510 {

511 }

size_t lar_cluster3d::Cluster3D::FindAndStoreDaughters	(	ArtOutputHandler &	output,
		reco::ClusterParameters &	clusterParameters,
		size_t	pfParticleParent,
		IdxToPCAMap &	idxToPCAMap,
		RecobHitToPtrMap &	hitToPtrMap,
		Hit3DToSPPtrMap &	hit3DToSPPtrMap
	)		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 1265 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(output, clusterParams, pfParticleParent, idxToPCAMap, hitToPtrMap, hit3DToSPPtrMap);
     }
     // Otherwise we want to store the information
     else
     {
         size_t daughterIdx = ConvertToArtOutput(output, clusterParameters, pfParticleParent, hitToPtrMap, hit3DToSPPtrMap);
         
         idxToPCAMap[daughterIdx] = &clusterParameters.getFullPCA();
     }
         
     return idxToPCAMap.size();
 }

void lar_cluster3d::Cluster3D::findTrackSeeds	(	art::Event &	evt,
		reco::ClusterParameters &	cluster,
		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 631 of file Cluster3D_module.cc.

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

 {
     // 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  = sqrt(std::pow(eigenVal1,2) + std::pow(eigenVal2,2));
     
     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 (eigenVal0 > 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(*this, evt, seedVec, seedHitVec, seedHitAssns);
     }
     
     return;
 }

EngineCreator::seed_t EngineCreator::get_seed_value	(	fhicl::ParameterSet const &	pset,
		char const	key[] = `"seed"`,
		seed_t const	implicit_seed = `-1`
	)

inherited

Definition at line 49 of file EngineCreator.cc.

References fhicl::ParameterSet::get().

Referenced by art::MixFilter< T >::initEngine_().

 {
   auto const& explicit_seeds = pset.get<std::vector<int>>(key, {});
   return explicit_seeds.empty() ? implicit_seed : explicit_seeds.front();
 }

template<typename PROD , BranchType B>

ProductID art::EDProducer::getProductID ( std::string const & instanceName = {} ) const

inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

   {
     return ProducerBase::getProductID<PROD, B>(moduleDescription_,
                                                instanceName);
   }

template<typename PROD , BranchType B>

ProductID art::ProducerBase::getProductID	(	ModuleDescription const &	moduleDescription,
		std::string const &	instanceName
	)		const

inherited

Definition at line 56 of file ProducerBase.h.

References B, and art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

   {
     auto const& pd =
       get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
     return pd.productID();
   }

void lar_cluster3d::Cluster3D::InitializeMonitoring ( )

private

Initialize the internal monitoring.

Definition at line 597 of file Cluster3D_module.cc.

References m_artHitsTime, m_buildNeighborhoodTime, m_dbscanTime, m_event, m_finishTime, m_hits, m_makeHitsTime, m_pathFindingTime, m_pRecoTree, m_run, m_totalTime, and art::TFileDirectory::make().

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("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("pathfindingtime",      &m_pathFindingTime,       "time/F");
     m_pRecoTree->Branch("finishTime",           &m_finishTime,            "time/F");
 }

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

private

Definition at line 1658 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()[0][0], fullPCA.getEigenVectors()[0][1], fullPCA.getEigenVectors()[0][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(DocaToPCAPair(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(NULL);
     
     // 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;
     }
     
     return;
 }

void lar_cluster3d::Cluster3D::MakeAndSaveSpacePoints	(	ArtOutputHandler &	output,
		reco::HitPairListPtr &	clusHitPairVector,
		RecobHitToPtrMap &	hitToPtrMap,
		Hit3DToSPPtrMap &	hit3DToSPPtrMap,
		int	spacePointStart
	)		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 1538 of file Cluster3D_module.cc.

References lar_cluster3d::Cluster3D::ArtOutputHandler::artSpacePointVector, reco::ClusterHit3D::MADESPACEPOINT, lar_cluster3d::Cluster3D::ArtOutputHandler::makePFPartSpacePointAssns(), lar_cluster3d::Cluster3D::ArtOutputHandler::makeSpacePointHitAssns(), and reco::ClusterHit3D::REJECTEDHIT.

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& 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...
         
 //        if      ( hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) && !hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -1.;  // pure skeleton point
 //        else if (!hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) &&  hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -2.;  // pure edge point
 //        else if ( hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) &&  hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) chisq = -3.;  // skeleton and edge point
         
 //        if      (hitPair->bitsAreSet(reco::ClusterHit3D::SEEDHIT)                                                          ) chisq = -4.;  // Seed point
         
 //        if ((hitPair->getStatusBits() & 0x7) != 0x7) chisq = -10.;
         
         // Mark this hit pair as in use
         hitPair->setStatusBit(reco::ClusterHit3D::MADESPACEPOINT);
         
         // Create and store the space point
         size_t spacePointID    = output.artSpacePointVector->size();
         double spacePointPos[] = {hitPair->getPosition()[0],hitPair->getPosition()[1],hitPair->getPosition()[2]};
         output.artSpacePointVector->push_back(recob::SpacePoint(spacePointPos, spError, chisq, output.artSpacePointVector->size()));
         
         // Update mapping
         hit3DToSPPtrMap[hitPair] = spacePointID;
         
         // 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(recobHits);
     }
     
     output.makePFPartSpacePointAssns(spacePointStart);
 
     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 1595 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.at(fromVertex), vertexToSpacePointMap.at(toVertex), output.artEdgeVector->size());
             
             halfEdgeSet.insert(&halfEdge);
         }
     }
 
     return;
 }

template<typename T , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ProductToken<T> art::Consumer::mayConsume ( InputTag const & it )

inherited

Definition at line 190 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ProductToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::Product,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ProductToken<T>{it};
 }

template<typename T , art::BranchType BT>

void art::Consumer::mayConsumeMany ( )

inherited

Definition at line 205 of file Consumer.h.

 {
   if (!moduleContext_)
     return;
 
   consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
 }

template<typename Element , BranchType = InEvent>

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

inherited

template<typename T , art::BranchType BT>

art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const & it )

inherited

Definition at line 215 of file Consumer.h.

References art::InputTag::instance(), art::InputTag::label(), and art::InputTag::process().

 {
   if (!moduleContext_)
     return ViewToken<T>::invalid();
 
   consumables_[BT].emplace_back(ConsumableType::ViewElement,
                                 TypeID{typeid(T)},
                                 it.label(),
                                 it.instance(),
                                 it.process());
   return ViewToken<T>{it};
 }

bool art::ProducerBase::modifiesEvent ( ) const

inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID().

     {
       return true;
     }

cet::exempt_ptr< art::Consumer > art::Consumer::non_module_context ( )

staticinherited

Definition at line 76 of file Consumer.cc.

 {
   static Consumer invalid{InvalidTag{}};
   return &invalid;
 }

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

private

Event Preparation.

Parameters

evt	the ART event

Definition at line 615 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_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_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 art::Consumer::prepareForJob ( fhicl::ParameterSet const & pset )

protectedinherited

Definition at line 89 of file Consumer.cc.

References fhicl::ParameterSet::get_if_present().

Referenced by art::EDProducer::doBeginJob(), art::EDFilter::doBeginJob(), and art::EDAnalyzer::doBeginJob().

 {
   if (!moduleContext_)
     return;
 
   pset.get_if_present("errorOnMissingConsumes", requireConsumes_);
   for (auto& consumablesPerBranch : consumables_) {
     cet::sort_all(consumablesPerBranch);
   }
 }

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

virtual

Producer method for reovering the 2D hits and driving the 3D reconstruciton

Implements art::EDProducer.

Definition at line 515 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_clusterPathAlg, m_dbscanTime, m_enableMonitoring, m_event, m_finishTime, m_hit3DBuilderAlg, m_hits, m_makeHitsTime, m_pathFindingTime, m_pRecoTree, m_run, m_spacePointInstance, m_totalTime, lar_cluster3d::Cluster3D::ArtOutputHandler::outputObjects(), lar_cluster3d::IClusterAlg::PATHFINDING, 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
     m_hit3DBuilderAlg->Hit3DBuilder(evt, *hitPairList, clusterHitToArtPtrMap);
     
     std::cout << "++> Produced: " << hitPairList->size() << " hits" << std::endl;
         
     // Call the main workhorse algorithm for building the local version of candidate 3D clusters
     m_clusterAlg->Cluster3DHits(*hitPairList, clusterParametersList);
     
     std::cout << "++> Produced: " << clusterParametersList.size() << " clusters" << std::endl;
         
     // 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(*this, evt, m_spacePointInstance);
     
     std::cout << "++> Outputting clusters" << std::endl;
 
     // 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
     ProduceArtClusters(output, *hitPairList, clusterParametersList, clusterHitToArtPtrMap);
     
     // Output to art
     output.outputObjects();
     
     if (m_enableMonitoring) theClockFinish.stop();
     
     // If monitoring then deal with the fallout
     if (m_enableMonitoring)
     {
         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_pathFindingTime       = m_clusterAlg->getTimeToExecute(IClusterAlg::PATHFINDING);
         m_finishTime            = theClockFinish.accumulated_real_time();
         m_hits                  = static_cast<int>(clusterHitToArtPtrMap.size());
         m_pRecoTree->Fill();
         
         mf::LogDebug("Cluster3D") << "*** Cluster3D total time: " << m_totalTime << ", art: " << m_artHitsTime << ", make: " << m_makeHitsTime
         << ", build: " << m_buildNeighborhoodTime << ", clustering: " << m_dbscanTime << ", path: " << m_pathFindingTime << ", finish: " << m_finishTime << std::endl;
     }
     
     // Will we ever get here? ;-)
     return;
 }

void lar_cluster3d::Cluster3D::ProduceArtClusters	(	ArtOutputHandler &	output,
		reco::HitPairList &	hitPairList,
		reco::ClusterParametersList &	clusterParametersList,
		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 1065 of file Cluster3D_module.cc.

References lar_cluster3d::Cluster3D::ArtOutputHandler::artEdgeVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artPCAxisVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artPFParticleVector, lar_cluster3d::Cluster3D::ArtOutputHandler::artSpacePointVector, ConvertToArtOutput(), countUltimateDaughters(), FindAndStoreDaughters(), reco::PrincipalComponents::getAveHitDoca(), reco::PrincipalComponents::getAvePosition(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), reco::PrincipalComponents::getNumHitsUsed(), reco::PrincipalComponents::getSvdOK(), recob::PFParticle::kPFParticlePrimary, reco::ClusterHit3D::MADESPACEPOINT, MakeAndSaveSpacePoints(), MakeAndSaveVertexPoints(), lar_cluster3d::Cluster3D::ArtOutputHandler::makePFPartEdgeAssns(), lar_cluster3d::Cluster3D::ArtOutputHandler::makePFPartPCAAssns(), and lar_cluster3d::Cluster3D::ArtOutputHandler::makeSpacePointHitAssns().

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;
             
             // Keep track of current start for space points
             int spacePointStart(output.artSpacePointVector->size());
             
             // Do a special output of voronoi vertices here...
             dcel2d::VertexList&   vertexList   = clusterParameters.getVertexList();
             dcel2d::HalfEdgeList& halfEdgeList = clusterParameters.getHalfEdgeList();
             
             std::cout << "Preparing to save the vertex point list, size: " << vertexList.size() << ", half edges: " << halfEdgeList.size() << std::endl;
             
             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(output, clusterParameters, recob::PFParticle::kPFParticlePrimary, hitToPtrMap, hit3DToSPPtrMap);
             }
             // 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(output, clusterParameters, pfParticleIdx, idxToPCAMap, hitToPtrMap, hit3DToSPPtrMap);
                 
                 // Now make the piecewise curve
 //                MakeAndSavePCAPoints(output, clusterParameters.getFullPCA(), idxToPCAMap);
 
                 // 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);
                     
                     for(size_t innerIdx = 0; innerIdx < 3; innerIdx++) eigenVecs[outerIdx][innerIdx] = skeletonPCA.getEigenVectors()[outerIdx][innerIdx];
                 }
                 
                 
                 recob::PCAxis skelPcAxis(skeletonPCA.getSvdOK(),
                                          skeletonPCA.getNumHitsUsed(),
                                          eigenVals,                      //skeletonPCA.getEigenValues(),
                                          eigenVecs,                      //skeletonPCA.getEigenVectors(),
                                          avePosition,                    //skeletonPCA.getAvePosition(),
                                          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()[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();
                 
                 // Make associations to all space points for this cluster
                 MakeAndSaveSpacePoints(output, clusterParameters.getHitPairListPtr(), hitToPtrMap, hit3DToSPPtrMap, spacePointStart);
 
                 // Build the edges now
                 size_t edgeStart(output.artEdgeVector->size());
                 
                 for(const auto& edge : clusterParameters.getBestEdgeList())
                 {
                     Hit3DToSPPtrMap::iterator hit0Itr = hit3DToSPPtrMap.find(std::get<0>(edge));
                     Hit3DToSPPtrMap::iterator hit1Itr = hit3DToSPPtrMap.find(std::get<1>(edge));
                     
                     bool hit0Found = hit0Itr != hit3DToSPPtrMap.end();
                     bool hit1Found = hit1Itr != hit3DToSPPtrMap.end();
                     
                     if (!hit0Found || !hit1Found) std::cout << "<<<<< Did not find matching space point " << hit0Found << ", " << hit1Found << " >>>>>>" << std::endl;
 
                     output.artEdgeVector->push_back(recob::Edge(std::get<2>(edge), hit3DToSPPtrMap[std::get<0>(edge)], hit3DToSPPtrMap[std::get<1>(edge)], output.artEdgeVector->size()));
                 }
                 
                 output.makePFPartEdgeAssns(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(-100.);
         
         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++;
         
         output.artSpacePointVector->push_back(recob::SpacePoint(spacePointPos, spacePointErr, chisq, output.artSpacePointVector->size()));
         
         if (!recobHits.empty()) output.makeSpacePointHitAssns(recobHits);
     }
     
     std::cout << "++++>>>> total num hits: " << hitPairVector.size() << ", num free: " << nFreePoints << std::endl;
 
     return;
 }

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

Definition at line 472 of file Cluster3D_module.cc.

References fhicl::ParameterSet::get(), m_clusterAlg, m_clusterMergeAlg, m_clusterPathAlg, m_enableMonitoring, m_hit3DBuilderAlg, m_parallelHitsAlg, m_parallelHitsCosAng, m_parallelHitsTransWid, m_pcaAlg, m_pcaSeedFinderAlg, m_seedFinderAlg, m_skeletonAlg, lar_cluster3d::HoughSeedFinderAlg::reconfigure(), lar_cluster3d::SkeletonAlg::reconfigure(), lar_cluster3d::PCASeedFinderAlg::reconfigure(), lar_cluster3d::ParallelHitsSeedFinderAlg::reconfigure(), and lar_cluster3d::PrincipalComponentsAlg::reconfigure().

Referenced by Cluster3D().

 {
     m_enableMonitoring     = pset.get<bool>       ("EnableMonitoring",         false);
     m_parallelHitsCosAng   = pset.get<float>      ("ParallelHitsCosAng",       0.999);
     m_parallelHitsTransWid = pset.get<float>      ("ParallelHitsTransWid",      25.0);
     
     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_pcaAlg.reconfigure(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"));
     m_skeletonAlg.reconfigure(pset.get<fhicl::ParameterSet>("SkeletonAlg"));
     m_seedFinderAlg.reconfigure(pset.get<fhicl::ParameterSet>("SeedFinderAlg"));
     m_pcaSeedFinderAlg.reconfigure(pset.get<fhicl::ParameterSet>("PCASeedFinderAlg"));
     m_parallelHitsAlg.reconfigure(pset.get<fhicl::ParameterSet>("ParallelHitsAlg"));
 }

void art::Consumer::showMissingConsumes ( ) const

protectedinherited

Definition at line 125 of file Consumer.cc.

Referenced by art::EDProducer::doEndJob(), art::EDFilter::doEndJob(), art::EDAnalyzer::doEndJob(), and art::RootOutput::endJob().

 {
   if (!moduleContext_)
     return;
 
   // If none of the branches have missing consumes statements, exit early.
   if (std::all_of(cbegin(missingConsumes_),
                   cend(missingConsumes_),
                   [](auto const& perBranch) { return perBranch.empty(); }))
     return;
 
   constexpr cet::HorizontalRule rule{60};
   mf::LogPrint log{"MTdiagnostics"};
   log << '\n'
       << rule('=') << '\n'
       << "The following consumes (or mayConsume) statements are missing from\n"
       << module_context(moduleDescription_) << '\n'
       << rule('-') << '\n';
 
   cet::for_all_with_index(
     missingConsumes_, [&log](std::size_t const i, auto const& perBranch) {
       for (auto const& pi : perBranch) {
         log << "  "
             << assemble_consumes_statement(static_cast<BranchType>(i), pi)
             << '\n';
       }
     });
   log << rule('=');
 }

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

References lar_cluster3d::SkeletonAlg::AverageSkeletonPositions(), lar_cluster3d::ClusterParamsBuilder::FillClusterParams(), 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;
     }
 
     return;
 }

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

References lar_cluster3d::SkeletonAlg::AverageSkeletonPositions(), reco::ClusterHit3D::bitsAreSet(), reco::ClusterHit3D::clearStatusBits(), reco::PrincipalComponents::getEigenValues(), reco::ClusterParameters::getHitPairListPtr(), reco::ClusterHit3D::getPosition(), lar_cluster3d::SkeletonAlg::GetSkeletonHits(), reco::ClusterParameters::getSkeletonPCA(), art::left(), m_skeletonAlg, art::right(), reco::ClusterHit3D::SELECTEDBYMST, and reco::ClusterHit3D::setStatusBit().

 {
     // 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(DistanceHit3DPair(hitDistance,hit3DOuter));
             outerHitList.emplace_back(DistanceHit3DPair(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;
                 }
                 else 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(DistanceEdgePair(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());
     
     return;
 }

void art::Consumer::validateConsumedProduct	(	BranchType const	bt,
		ProductInfo const &	pi
	)

protectedinherited

Definition at line 101 of file Consumer.cc.

References art::errors::ProductRegistrationFailure.

 {
   // Early exits if consumes tracking has been disabled or if the
   // consumed product is an allowed consumable.
   if (!moduleContext_)
     return;
 
   if (cet::binary_search_all(consumables_[bt], pi))
     return;
 
   if (requireConsumes_) {
     throw Exception(errors::ProductRegistrationFailure,
                     "Consumer: an error occurred during validation of a "
                     "retrieved product\n\n")
       << "The following consumes (or mayConsume) statement is missing from\n"
       << module_context(moduleDescription_) << ":\n\n"
       << "  " << assemble_consumes_statement(bt, pi) << "\n\n";
   }
 
   missingConsumes_[bt].insert(pi);
 }