StandardHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::StandardHit3DBuilder:

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

void	produces (art::ProducesCollector &) override
	Each algorithm may have different objects it wants "produced" so use this to let the top level producer module "know" what it is outputting. More...

void	configure (const fhicl::ParameterSet &) override
	Interface for configuring the particular algorithm tool. More...

void	Hit3DBuilder (art::Event &, reco::HitPairList &, RecobHitToPtrMap &) override
	Given a set of recob hits, run DBscan to form 3D clusters. More...

float	getTimeToExecute (IHit3DBuilder::TimeValues index) const override
	If monitoring, recover the time to execute a particular function. More...

Public Member Functions inherited from lar_cluster3d::IHit3DBuilder
virtual	~IHit3DBuilder () noexcept=default
	Virtual Destructor. More...

Private Types
using	PlaneHitVectorItrPairVec = std::vector< std::pair< HitVector::iterator, HitVector::iterator >>
	Given the ClusterHit2D objects, build the HitPairMap. More...

using	HitMatchPair = std::pair< const reco::ClusterHit2D *, reco::ClusterHit3D >
	This builds a list of candidate hit pairs from lists of hits on two planes. More...

using	HitMatchPairVec = std::vector< HitMatchPair >

using	HitMatchPairVecMap = std::map< geo::WireID, HitMatchPairVec >

using	ChannelStatusVec = std::vector< size_t >
	define data structure for keeping track of channel status More...

using	ChannelStatusByPlaneVec = std::vector< ChannelStatusVec >

Private Member Functions
void	CollectArtHits (const art::Event &evt) const
	Extract the ART hits and the ART hit-particle relationships. More...

void	BuildHit3D (reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...

void	CreateNewRecobHitCollection (art::Event &, reco::HitPairList &, std::vector< recob::Hit > &, RecobHitToPtrMap &)
	Create a new 2D hit collection from hits associated to 3D space points. More...

void	makeWireAssns (const art::Event &, art::Assns< recob::Wire, recob::Hit > &, RecobHitToPtrMap &) const
	Create recob::Wire to recob::Hit associations. More...

void	makeRawDigitAssns (const art::Event &, art::Assns< raw::RawDigit, recob::Hit > &, RecobHitToPtrMap &) const
	Create raw::RawDigit to recob::Hit associations. More...

size_t	BuildHitPairMap (PlaneToHitVectorMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...

size_t	BuildHitPairMapByTPC (PlaneHitVectorItrPairVec &planeHitVectorItrPairVec, reco::HitPairList &hitPairList) const

int	findGoodHitPairs (const reco::ClusterHit2D *, HitVector::iterator &, HitVector::iterator &, HitMatchPairVecMap &) const

void	findGoodTriplets (HitMatchPairVecMap &, HitMatchPairVecMap &, reco::HitPairList &, bool=false) const
	This algorithm takes lists of hit pairs and finds good triplets. More...

bool	makeHitPair (reco::ClusterHit3D &pairOut, const reco::ClusterHit2D hit1, const reco::ClusterHit2D hit2, float hitWidthSclFctr=1., size_t hitPairCntr=0) const
	Make a HitPair object by checking two hits. More...

bool	makeHitTriplet (reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pairIn, const reco::ClusterHit2D *hit2) const
	Make a 3D HitPair object by checking two hits. More...

bool	makeDeadChannelPair (reco::ClusterHit3D &pairOut, const reco::ClusterHit3D &pair, size_t maxStatus=4, size_t minStatus=0, float minOverlap=0.2) const
	Make a 3D HitPair object from a valid pair and a dead channel in the missing plane. More...

const reco::ClusterHit2D *	FindBestMatchingHit (const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, float pairDeltaTimeLimits) const
	A utility routine for finding a 2D hit closest in time to the given pair. More...

int	FindNumberInRange (const Hit2DSet &hit2DSet, const reco::ClusterHit3D &pair, float range) const
	A utility routine for returning the number of 2D hits from the list in a given range. More...

geo::WireID	NearestWireID (const Eigen::Vector3f &position, const geo::WireID &wireID) const
	Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range. More...

float	DistanceFromPointToHitWire (const Eigen::Vector3f &position, const geo::WireID &wireID) const
	Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range. More...

void	BuildChannelStatusVec (PlaneToWireToHitSetMap &planeToWiretoHitSetMap) const
	Create the internal channel status vector (assume will eventually be event-by-event) More...

float	chargeIntegral (float, float, float, float, int, int) const
	Perform charge integration between limits. More...

void	clear ()
	clear the tuple vectors before processing next event More...

Private Attributes
std::vector< art::InputTag >	m_hitFinderTagVec
	Data members to follow. More...

float	m_numSigmaPeakTime

float	m_hitWidthSclFctr

float	m_deltaPeakTimeSig

std::vector< int >	m_invalidTPCVec

float	m_wirePitchScaleFactor
	Scaling factor to determine max distance allowed between candidate pairs. More...

float	m_maxHit3DChiSquare
	Provide ability to select hits based on "chi square". More...

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

bool	m_enableMonitoring

float	m_wirePitch [3]

std::vector< float >	m_timeVector

float	m_zPosOffset

TTree *	m_tupleTree
	output analysis tree More...

std::vector< float >	m_deltaTimeVec

std::vector< float >	m_chiSquare3DVec

std::vector< float >	m_maxPullVec

std::vector< float >	m_overlapFractionVec

std::vector< float >	m_overlapRangeVec

std::vector< float >	m_maxDeltaPeakVec

std::vector< float >	m_maxSideVecVec

std::vector< float >	m_pairWireDistVec

std::vector< float >	m_smallChargeDiffVec

std::vector< int >	m_smallIndexVec

std::vector< float >	m_qualityMetricVec

std::vector< float >	m_spacePointChargeVec

std::vector< float >	m_hitAsymmetryVec

Hit2DList	m_clusterHit2DMasterList

PlaneToHitVectorMap	m_planeToHitVectorMap

PlaneToWireToHitSetMap	m_planeToWireToHitSetMap

ChannelStatusByPlaneVec	m_channelStatus

size_t	m_numBadChannels

bool	m_weHaveAllBeenHereBefore = false

const geo::Geometry *	m_geometry

const lariov::ChannelStatusProvider *	m_channelFilter

Additional Inherited Members
Public Types inherited from lar_cluster3d::IHit3DBuilder
enum	TimeValues { COLLECTARTHITS = 0, BUILDTHREEDHITS = 1, BUILDNEWHITS = 2, NUMTIMEVALUES }
	enumerate the possible values for time checking if monitoring timing More...

using	RecobHitToPtrMap = std::unordered_map< const recob::Hit *, art::Ptr< recob::Hit >>
	Defines a structure mapping art representation to internal. More...

Detailed Description

StandardHit3DBuilder class definiton.

Definition at line 77 of file StandardHit3DBuilder_tool.cc.

Member Typedef Documentation

using lar_cluster3d::StandardHit3DBuilder::ChannelStatusByPlaneVec = std::vector<ChannelStatusVec>

private

Definition at line 245 of file StandardHit3DBuilder_tool.cc.

using lar_cluster3d::StandardHit3DBuilder::ChannelStatusVec = std::vector<size_t>

private

define data structure for keeping track of channel status

Definition at line 244 of file StandardHit3DBuilder_tool.cc.

using lar_cluster3d::StandardHit3DBuilder::HitMatchPair = std::pair<const reco::ClusterHit2D*, reco::ClusterHit3D>

private

This builds a list of candidate hit pairs from lists of hits on two planes.

Definition at line 164 of file StandardHit3DBuilder_tool.cc.

using lar_cluster3d::StandardHit3DBuilder::HitMatchPairVec = std::vector<HitMatchPair>

private

Definition at line 165 of file StandardHit3DBuilder_tool.cc.

using lar_cluster3d::StandardHit3DBuilder::HitMatchPairVecMap = std::map<geo::WireID, HitMatchPairVec>

private

Definition at line 166 of file StandardHit3DBuilder_tool.cc.

using lar_cluster3d::StandardHit3DBuilder::PlaneHitVectorItrPairVec = std::vector<std::pair<HitVector::iterator, HitVector::iterator>>

private

Given the ClusterHit2D objects, build the HitPairMap.

Definition at line 156 of file StandardHit3DBuilder_tool.cc.

Constructor & Destructor Documentation

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

explicit

Constructor.

Parameters

pset

Definition at line 302 of file StandardHit3DBuilder_tool.cc.

     : m_geometry(art::ServiceHandle<geo::Geometry const>{}.get())
     , m_channelFilter(&art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider())
   {
     this->configure(pset);
   }

Member Function Documentation

void lar_cluster3d::StandardHit3DBuilder::BuildChannelStatusVec ( PlaneToWireToHitSetMap & planeToWiretoHitSetMap ) const

private

Create the internal channel status vector (assume will eventually be event-by-event)

Definition at line 391 of file StandardHit3DBuilder_tool.cc.

   {
     // This is called each event, clear out the previous version and start over
     if (!m_channelStatus.empty()) m_channelStatus.clear();
 
     m_numBadChannels = 0;
     m_channelStatus.resize(m_geometry->Nplanes());
 
     // Loop through views/planes to set the wire length vectors
     for (size_t idx = 0; idx < m_channelStatus.size(); idx++) {
       m_channelStatus[idx] = ChannelStatusVec(m_geometry->Nwires(idx), 5);
     }
 
     // Loop through the channels and mark those that are "bad"
     for (size_t channel = 0; channel < m_geometry->Nchannels(); channel++) {
       if (!m_channelFilter->IsGood(channel)) {
         std::vector<geo::WireID> wireIDVec = m_geometry->ChannelToWire(channel);
         geo::WireID wireID = wireIDVec[0];
         lariov::ChannelStatusProvider::Status_t chanStat = m_channelFilter->Status(channel);
 
         m_channelStatus[wireID.Plane][wireID.Wire] = chanStat;
         m_numBadChannels++;
       }
     }
 
     // add quiet wires in U plane for microboone (this will done "correctly" in near term)
     //    PlaneToWireToHitSetMap::iterator plane0HitItr = planeToWireToHitSetMap.find(geo::PlaneID(0,0,0));
 
     //    if (plane0HitItr != planeToWireToHitSetMap.end())
     //    {
     ////        WireToHitSetMap& wireToHitSetMap = uPlaneHitItr->second;
 
     //        for(size_t idx = 2016; idx < 2096; idx++)  m_channelStatus[0][idx] = 3;
     //        for(size_t idx = 2192; idx < 2304; idx++)  m_channelStatus[0][idx] = 3;
     //        for(size_t idx = 2352; idx < 2384; idx++)  m_channelStatus[0][idx] = 3;
     //        //for(size_t idx = 2016; idx < 2096; idx++) if (wireToHitSetMap.find(idx) == wireToHitSetMap.end()) m_channelStatus[0][idx] = 3;
     //        //for(size_t idx = 2192; idx < 2304; idx++) if (wireToHitSetMap.find(idx) == wireToHitSetMap.end()) m_channelStatus[0][idx] = 3;
     //        //for(size_t idx = 2352; idx < 2384; idx++) if (wireToHitSetMap.find(idx) == wireToHitSetMap.end()) m_channelStatus[0][idx] = 3;
     ////      for(size_t idx = 2016; idx < 2384; idx++) m_channelStatus[0][idx] = 3;
     //    }
 
     return;
   }

void lar_cluster3d::StandardHit3DBuilder::BuildHit3D ( reco::HitPairList & hitPairList ) const

private

Given the ClusterHit2D objects, build the HitPairMap.

Driver for processing input 2D hits, transforming to 3D hits and building lists of associated 3D hits (candidate 3D clusters)

Definition at line 512 of file StandardHit3DBuilder_tool.cc.

   {
     /**
      *  @brief Driver for processing input 2D hits, transforming to 3D hits and building lists
      *         of associated 3D hits (candidate 3D clusters)
      */
     cet::cpu_timer theClockMakeHits;
 
     if (m_enableMonitoring) theClockMakeHits.start();
 
     // The first task is to take the lists of input 2D hits (a map of view to sorted lists of 2D hits)
     // and then to build a list of 3D hits to be used in downstream processing
     BuildChannelStatusVec(m_planeToWireToHitSetMap);
 
     size_t numHitPairs = BuildHitPairMap(m_planeToHitVectorMap, hitPairList);
 
     if (m_enableMonitoring) {
       theClockMakeHits.stop();
 
       m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
     }
 
     mf::LogDebug("Cluster3D") << ">>>>> 3D hit building done, found " << numHitPairs << " 3D Hits"
                               << std::endl;
 
     return;
   }

size_t lar_cluster3d::StandardHit3DBuilder::BuildHitPairMap	(	PlaneToHitVectorMap &	planeToHitVectorMap,
		reco::HitPairList &	hitPairList
	)		const

private

Given the ClusterHit2D objects, build the HitPairMap.

Given input 2D hits, build out the lists of possible 3D hits

   The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
   However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
   be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
   and then look for the match in the V plane. In the event we don't find the match in the V plane then we
   will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.

Definition at line 595 of file StandardHit3DBuilder_tool.cc.

   {
     /**
      *  @brief Given input 2D hits, build out the lists of possible 3D hits
      *
      *         The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
      *         However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
      *         be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
      *         and then look for the match in the V plane. In the event we don't find the match in the V plane then we
      *         will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
      */
     size_t totalNumHits(0);
     size_t hitPairCntr(0);
 
     size_t nTriplets(0);
     size_t nDeadChanHits(0);
 
     // Set up to loop over cryostats and tpcs...
     for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) {
       for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) {
         PlaneToHitVectorMap::iterator mapItr0 =
           planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0));
         PlaneToHitVectorMap::iterator mapItr1 =
           planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1));
         PlaneToHitVectorMap::iterator mapItr2 =
           planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2));
 
         size_t nPlanesWithHits =
           (mapItr0 != planeToHitVectorMap.end() && !mapItr0->second.empty() ? 1 : 0) +
           (mapItr1 != planeToHitVectorMap.end() && !mapItr1->second.empty() ? 1 : 0) +
           (mapItr2 != planeToHitVectorMap.end() && !mapItr2->second.empty() ? 1 : 0);
 
         if (nPlanesWithHits < 2) continue;
 
         HitVector& hitVector0 = mapItr0->second;
         HitVector& hitVector1 = mapItr1->second;
         HitVector& hitVector2 = mapItr2->second;
 
         // We are going to resort the hits into "start time" order...
         std::sort(hitVector0.begin(),
                   hitVector0.end(),
                   SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);
         std::sort(hitVector1.begin(),
                   hitVector1.end(),
                   SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);
         std::sort(hitVector2.begin(),
                   hitVector2.end(),
                   SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);
 
         PlaneHitVectorItrPairVec hitItrVec = {
           HitVectorItrPair(hitVector0.begin(), hitVector0.end()),
           HitVectorItrPair(hitVector1.begin(), hitVector1.end()),
           HitVectorItrPair(hitVector2.begin(), hitVector2.end())};
 
         totalNumHits += BuildHitPairMapByTPC(hitItrVec, hitPairList);
       }
     }
 
     // Return the hit pair list but sorted by z and y positions (faster traversal in next steps)
     hitPairList.sort(SetPairStartTimeOrder);
 
     // Where are we?
     mf::LogDebug("Cluster3D") << "Total number hits: " << totalNumHits << std::endl;
     mf::LogDebug("Cluster3D") << "Created a total of " << hitPairList.size()
                               << " hit pairs, counted: " << hitPairCntr << std::endl;
     mf::LogDebug("Cluster3D") << "-- Triplets: " << nTriplets
                               << ", dead channel pairs: " << nDeadChanHits << std::endl;
 
     return hitPairList.size();
   }

size_t lar_cluster3d::StandardHit3DBuilder::BuildHitPairMapByTPC	(	PlaneHitVectorItrPairVec &	planeHitVectorItrPairVec,
		reco::HitPairList &	hitPairList
	)		const

private

Given input 2D hits, build out the lists of possible 3D hits

   The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
   However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
   be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
   and then look for the match in the V plane. In the event we don't find the match in the V plane then we
   will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.

Definition at line 668 of file StandardHit3DBuilder_tool.cc.

   {
     /**
      *  @brief Given input 2D hits, build out the lists of possible 3D hits
      *
      *         The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view
      *         However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be
      *         be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes
      *         and then look for the match in the V plane. In the event we don't find the match in the V plane then we
      *         will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high.
      */
 
     // Define functions to set start/end iterators in the loop below
     auto SetStartIterator =
       [](HitVector::iterator startItr, HitVector::iterator endItr, float rms, float startTime) {
         while (startItr != endItr) {
           float numRMS(rms);
           if ((*startItr)->getTimeTicks() + numRMS * (*startItr)->getHit()->RMS() < startTime)
             startItr++;
           else
             break;
         }
         return startItr;
       };
 
     auto SetEndIterator =
       [](HitVector::iterator firstItr, HitVector::iterator endItr, float rms, float endTime) {
         while (firstItr != endItr) {
           float numRMS(rms);
           if ((*firstItr)->getTimeTicks() - numRMS * (*firstItr)->getHit()->RMS() < endTime)
             firstItr++;
           else
             break;
         }
         return firstItr;
       };
 
     size_t nTriplets(0);
     size_t nDeadChanHits(0);
 
     //*********************************************************************************
     // Basically, we try to loop until done...
     while (1) {
       // Sort so that the earliest hit time will be the first element, etc.
       std::sort(hitItrVec.begin(), hitItrVec.end(), SetStartTimeOrder(m_numSigmaPeakTime));
 
       // This loop iteration's golden hit
       const reco::ClusterHit2D* goldenHit = *hitItrVec[0].first;
 
       // The range of history... (for this hit)
       float goldenTimeStart = goldenHit->getTimeTicks() -
                               m_numSigmaPeakTime * goldenHit->getHit()->RMS() -
                               std::numeric_limits<float>::epsilon();
       float goldenTimeEnd = goldenHit->getTimeTicks() +
                             m_numSigmaPeakTime * goldenHit->getHit()->RMS() +
                             std::numeric_limits<float>::epsilon();
 
       // Set iterators to insure we'll be in the overlap ranges
       HitVector::iterator hitItr1Start = SetStartIterator(
         hitItrVec[1].first, hitItrVec[1].second, m_numSigmaPeakTime, goldenTimeStart);
       HitVector::iterator hitItr1End =
         SetEndIterator(hitItr1Start, hitItrVec[1].second, m_numSigmaPeakTime, goldenTimeEnd);
       HitVector::iterator hitItr2Start = SetStartIterator(
         hitItrVec[2].first, hitItrVec[2].second, m_numSigmaPeakTime, goldenTimeStart);
       HitVector::iterator hitItr2End =
         SetEndIterator(hitItr2Start, hitItrVec[2].second, m_numSigmaPeakTime, goldenTimeEnd);
 
       // Since we'll use these many times in the internal loops, pre make the pairs for the second set of hits
       size_t curHitListSize(hitPairList.size());
       HitMatchPairVecMap pair12Map;
       HitMatchPairVecMap pair13Map;
 
       size_t n12Pairs = findGoodHitPairs(goldenHit, hitItr1Start, hitItr1End, pair12Map);
       size_t n13Pairs = findGoodHitPairs(goldenHit, hitItr2Start, hitItr2End, pair13Map);
 
       nDeadChanHits += hitPairList.size() - curHitListSize;
       curHitListSize = hitPairList.size();
 
       if (n12Pairs > n13Pairs)
         findGoodTriplets(pair12Map, pair13Map, hitPairList);
       else
         findGoodTriplets(pair13Map, pair12Map, hitPairList);
 
       nTriplets += hitPairList.size() - curHitListSize;
 
       hitItrVec[0].first++;
 
       int nPlanesWithHits(0);
 
       for (auto& pair : hitItrVec)
         if (pair.first != pair.second) nPlanesWithHits++;
 
       if (nPlanesWithHits < 2) break;
     }
 
     return hitPairList.size();
   }

float lar_cluster3d::StandardHit3DBuilder::chargeIntegral	(	float	peakMean,
		float	peakAmp,
		float	peakSigma,
		float	areaNorm,
		int	low,
		int	hi
	)		const

private

Perform charge integration between limits.

Definition at line 1301 of file StandardHit3DBuilder_tool.cc.

   {
     float integral(0);
 
     for (int sigPos = low; sigPos < hi; sigPos++) {
       float arg = (float(sigPos) - peakMean + 0.5) / peakSigma;
       integral += peakAmp * std::exp(-0.5 * arg * arg);
     }
 
     return integral;
   }

void lar_cluster3d::StandardHit3DBuilder::clear ( )

private

clear the tuple vectors before processing next event

Definition at line 371 of file StandardHit3DBuilder_tool.cc.

   {
     m_deltaTimeVec.clear();
     m_chiSquare3DVec.clear();
     m_maxPullVec.clear();
     m_overlapFractionVec.clear();
     m_overlapRangeVec.clear();
     m_maxDeltaPeakVec.clear();
     m_maxSideVecVec.clear();
     m_pairWireDistVec.clear();
     m_smallChargeDiffVec.clear();
     m_smallIndexVec.clear();
     m_qualityMetricVec.clear();
     m_spacePointChargeVec.clear();
     m_hitAsymmetryVec.clear();
 
     return;
   }

void lar_cluster3d::StandardHit3DBuilder::CollectArtHits ( const art::Event & evt ) const

private

Extract the ART hits and the ART hit-particle relationships.

Parameters

evt	- the ART event

Recover the 2D hits from art and fill out the local data structures for the 3D clustering

Definition at line 1538 of file StandardHit3DBuilder_tool.cc.

   {
     /**
      *  @brief Recover the 2D hits from art and fill out the local data structures for the 3D clustering
      */
 
     // Start by getting a vector of valid, non empty hit collections to make sure we really have something to do here...
     // Here is a container for the hits...
     std::vector<const recob::Hit*> recobHitVec;
 
     auto const clock_data =
       art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
     auto const det_prop =
       art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
 
     // Loop through the list of input sources
     for (const auto& inputTag : m_hitFinderTagVec) {
       art::Handle<std::vector<recob::Hit>> recobHitHandle;
       evt.getByLabel(inputTag, recobHitHandle);
 
       if (!recobHitHandle.isValid() || recobHitHandle->size() == 0) continue;
 
       recobHitVec.reserve(recobHitVec.size() + recobHitHandle->size());
 
       for (const auto& hit : *recobHitHandle)
         recobHitVec.push_back(&hit);
     }
 
     // If the vector is empty there is nothing to do
     if (recobHitVec.empty()) return;
 
     cet::cpu_timer theClockMakeHits;
 
     if (m_enableMonitoring) theClockMakeHits.start();
 
     // We'll want to correct the hit times for the plane offsets
     // (note this is already taken care of when converting to position)
     std::map<geo::PlaneID, double> planeOffsetMap;
 
     // Try to output a formatted string
     std::string debugMessage("");
 
     // Initialize the plane to hit vector map
     for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) {
       for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) {
         m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = HitVector();
         m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] = HitVector();
         m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] = HitVector();
 
         // What we want here are the relative offsets between the planes
         // Note that plane 0 is assumed the "first" plane and is the reference
         planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = 0.;
         planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] =
           det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1)) -
           det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));
         planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] =
           det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2)) -
           det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));
 
         // Should we provide output?
         if (!m_weHaveAllBeenHereBefore) {
           std::ostringstream outputString;
 
           outputString << "***> plane 0 offset: "
                        << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)]
                        << ", plane 1: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)]
                        << ", plane    2: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)]
                        << "\n";
           debugMessage += outputString.str();
           outputString << "     Det prop plane 0: "
                        << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0))
                        << ", plane 1: "
                        << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1))
                        << ", plane 2: "
                        << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2))
                        << ", Trig: " << trigger_offset(clock_data) << "\n";
           debugMessage += outputString.str();
         }
       }
     }
 
     if (!m_weHaveAllBeenHereBefore) {
       mf::LogDebug("Cluster3D") << debugMessage << std::endl;
 
       m_weHaveAllBeenHereBefore = true;
     }
 
     // Cycle through the recob hits to build ClusterHit2D objects and insert
     // them into the map
     for (const auto& recobHit : recobHitVec) {
       // Reject hits with negative charge, these are misreconstructed
       if (recobHit->Integral() < 0.) continue;
 
       // For some detectors we can have multiple wire ID's associated to a given channel.
       // So we recover the list of these wire IDs
       const std::vector<geo::WireID>& wireIDs = m_geometry->ChannelToWire(recobHit->Channel());
 
       // And then loop over all possible to build out our maps
       for (const auto& wireID : wireIDs) {
         // Check if this is an invalid TPC
         // (for example, in protoDUNE there are logical TPC's which see no signal)
         if (std::find(m_invalidTPCVec.begin(), m_invalidTPCVec.end(), wireID.TPC) !=
             m_invalidTPCVec.end())
           continue;
 
         // Note that a plane ID will define cryostat, TPC and plane
         const geo::PlaneID& planeID = wireID.planeID();
 
         double hitPeakTime(recobHit->PeakTime() - planeOffsetMap[planeID]);
         double xPosition(det_prop.ConvertTicksToX(
           recobHit->PeakTime(), planeID.Plane, planeID.TPC, planeID.Cryostat));
 
         m_clusterHit2DMasterList.emplace_back(0, 0., 0., xPosition, hitPeakTime, wireID, recobHit);
 
         m_planeToHitVectorMap[planeID].push_back(&m_clusterHit2DMasterList.back());
         m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back());
       }
     }
 
     // Make a loop through to sort the recover hits in time order
     for (auto& hitVectorMap : m_planeToHitVectorMap)
       std::sort(hitVectorMap.second.begin(), hitVectorMap.second.end(), SetHitTimeOrder);
 
     if (m_enableMonitoring) {
       theClockMakeHits.stop();
 
       m_timeVector[COLLECTARTHITS] = theClockMakeHits.accumulated_real_time();
     }
 
     mf::LogDebug("Cluster3D") << ">>>>> Number of ART hits: " << m_clusterHit2DMasterList.size()
                               << std::endl;
   }

void lar_cluster3d::StandardHit3DBuilder::configure ( const fhicl::ParameterSet & )

overridevirtual

Interface for configuring the particular algorithm tool.

Parameters

ParameterSet The input set of parameters for configuration

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 322 of file StandardHit3DBuilder_tool.cc.

   {
     m_hitFinderTagVec = pset.get<std::vector<art::InputTag>>(
       "HitFinderTagVec", std::vector<art::InputTag>() = {"gaushit"});
     m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
     m_numSigmaPeakTime = pset.get<float>("NumSigmaPeakTime", 3.);
     m_hitWidthSclFctr = pset.get<float>("HitWidthScaleFactor", 6.);
     m_deltaPeakTimeSig = pset.get<float>("DeltaPeakTimeSig", 1.7);
     m_zPosOffset = pset.get<float>("ZPosOffset", 0.0);
     m_invalidTPCVec = pset.get<std::vector<int>>("InvalidTPCVec", std::vector<int>());
     m_wirePitchScaleFactor = pset.get<float>("WirePitchScaleFactor", 1.9);
     m_maxHit3DChiSquare = pset.get<float>("MaxHitChiSquare", 6.0);
     m_outputHistograms = pset.get<bool>("OutputHistograms", false);
 
     m_geometry = art::ServiceHandle<geo::Geometry const>{}.get();
 
     // Returns the wire pitch per plane assuming they will be the same for all TPCs
     m_wirePitch[0] = m_geometry->WirePitch(0);
     m_wirePitch[1] = m_geometry->WirePitch(1);
     m_wirePitch[2] = m_geometry->WirePitch(2);
 
     // Access ART's TFileService, which will handle creating and writing
     // histograms and n-tuples for us.
     art::ServiceHandle<art::TFileService> tfs;
 
     if (m_outputHistograms) {
       m_tupleTree = tfs->make<TTree>("Hit3DBuilderTree", "Tree by StandardHit3DBuilder");
 
       clear();
 
       m_tupleTree->Branch("DeltaTime2D", "std::vector<float>", &m_deltaTimeVec);
       m_tupleTree->Branch("ChiSquare3D", "std::vector<float>", &m_chiSquare3DVec);
       m_tupleTree->Branch("MaxPullValue", "std::vector<float>", &m_maxPullVec);
       m_tupleTree->Branch("OverlapFraction", "std::vector<float>", &m_overlapFractionVec);
       m_tupleTree->Branch("OverlapRange", "std::vector<float>", &m_overlapRangeVec);
       m_tupleTree->Branch("MaxDeltaPeak", "std::vector<float>", &m_maxDeltaPeakVec);
       m_tupleTree->Branch("MaxSideVec", "std::vector<float>", &m_maxSideVecVec);
       m_tupleTree->Branch("PairWireDistVec", "std::vector<float>", &m_pairWireDistVec);
       m_tupleTree->Branch("SmallChargeDiff", "std::vector<float>", &m_smallChargeDiffVec);
       m_tupleTree->Branch("SmallChargeIdx", "std::vector<int>", &m_smallIndexVec);
       m_tupleTree->Branch("QualityMetric", "std::vector<float>", &m_qualityMetricVec);
       m_tupleTree->Branch("SPCharge", "std::vector<float>", &m_spacePointChargeVec);
       m_tupleTree->Branch("HitAsymmetry", "std::vector<float>", &m_hitAsymmetryVec);
     }
 
     return;
   }

void lar_cluster3d::StandardHit3DBuilder::CreateNewRecobHitCollection	(	art::Event &	event,
		reco::HitPairList &	hitPairList,
		std::vector< recob::Hit > &	hitPtrVec,
		RecobHitToPtrMap &	recobHitToPtrMap
	)

private

Create a new 2D hit collection from hits associated to 3D space points.

Definition at line 1674 of file StandardHit3DBuilder_tool.cc.

   {
     // Set up the timing
     cet::cpu_timer theClockBuildNewHits;
 
     if (m_enableMonitoring) theClockBuildNewHits.start();
 
     // We want to build a unique list of hits from the input 3D points which, we know, reuse 2D points frequently
     // At the same time we need to create a new 2D hit with the "correct" WireID and replace the old 2D hit in the
     // ClusterHit2D object with this new one... while keeping track of the use of the old ones. My head is spinning...
     // Declare a set which will allow us to keep track of those CusterHit2D objects we have seen already
     std::set<const reco::ClusterHit2D*> visitedHit2DSet;
 
     // Use this handy art utility to make art::Ptr objects to the new recob::Hits for use in the output phase
     art::PtrMaker<recob::Hit> ptrMaker(event);
 
     // Reserve enough memory to replace every recob::Hit we have considered (this is upper limit)
     hitPtrVec.reserve(m_clusterHit2DMasterList.size());
 
     // Scheme is to loop through all 3D hits, then through each associated ClusterHit2D object
     for (reco::ClusterHit3D& hit3D : hitPairList) {
       reco::ClusterHit2DVec& hit2DVec = hit3D.getHits();
 
       // The loop is over the index so we can recover the correct WireID to associate to the new hit when made
       for (size_t idx = 0; idx < hit3D.getHits().size(); idx++) {
         const reco::ClusterHit2D* hit2D = hit2DVec[idx];
 
         // Have we seen this 2D hit already?
         if (visitedHit2DSet.find(hit2D) == visitedHit2DSet.end()) {
           visitedHit2DSet.insert(hit2D);
 
           // Create and save the new recob::Hit with the correct WireID
           hitPtrVec.emplace_back(
             recob::HitCreator(*hit2D->getHit(), hit3D.getWireIDs()[idx]).copy());
 
           // Recover a pointer to it...
           recob::Hit* newHit = &hitPtrVec.back();
 
           // Create a mapping from this hit to an art Ptr representing it
           recobHitToPtrMap[newHit] = ptrMaker(hitPtrVec.size() - 1);
 
           // And set the pointer to this hit in the ClusterHit2D object
           const_cast<reco::ClusterHit2D*>(hit2D)->setHit(newHit);
         }
       }
     }
 
     size_t numNewHits = hitPtrVec.size();
 
     if (m_enableMonitoring) {
       theClockBuildNewHits.stop();
 
       m_timeVector[BUILDNEWHITS] = theClockBuildNewHits.accumulated_real_time();
     }
 
     mf::LogDebug("Cluster3D") << ">>>>> New output recob::Hit size: " << numNewHits << " (vs "
                               << m_clusterHit2DMasterList.size() << " input)" << std::endl;
 
     return;
   }

float lar_cluster3d::StandardHit3DBuilder::DistanceFromPointToHitWire	(	const Eigen::Vector3f &	position,
		const geo::WireID &	wireID
	)		const

private

Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range.

Definition at line 1482 of file StandardHit3DBuilder_tool.cc.

   {
     float distance;
 
     // Embed the call to the geometry's services nearest wire id method in a try-catch block
     try {
       // Get the wire endpoints
       Eigen::Vector3d wireStart;
       Eigen::Vector3d wireEnd;
 
       m_geometry->WireEndPoints(wireIDIn, &wireStart[0], &wireEnd[0]);
 
       // Want the hit position to have same x value as wire coordinates
       Eigen::Vector3d hitPosition(wireStart[0], position[1], position[2]);
 
       // Want the wire direction
       Eigen::Vector3d wireDir = wireEnd - wireStart;
 
       wireDir.normalize();
 
       // Get arc length to doca
       double arcLen = (hitPosition - wireStart).dot(wireDir);
 
       Eigen::Vector3d docaVec = hitPosition - (wireStart + arcLen * wireDir);
 
       distance = docaVec.norm();
     }
     catch (std::exception& exc) {
       // This can happen, almost always because the coordinates are **just** out of range
       mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - "
                                   << exc.what() << std::endl;
 
       // Assume extremum for wire number depending on z coordinate
       distance = 0.;
     }
 
     return distance;
   }

const reco::ClusterHit2D * lar_cluster3d::StandardHit3DBuilder::FindBestMatchingHit	(	const Hit2DSet &	hit2DSet,
		const reco::ClusterHit3D &	pair,
		float	pairDeltaTimeLimits
	)		const

private

A utility routine for finding a 2D hit closest in time to the given pair.

Definition at line 1398 of file StandardHit3DBuilder_tool.cc.

   {
     static const float minCharge(0.);
 
     const reco::ClusterHit2D* bestVHit(0);
 
     float pairAvePeakTime(pair.getAvePeakTime());
 
     // Idea is to loop through the input set of hits and look for the best combination
     for (const auto& hit2D : hit2DSet) {
       if (hit2D->getHit()->Integral() < minCharge) continue;
 
       float hitVPeakTime(hit2D->getTimeTicks());
       float deltaPeakTime(pairAvePeakTime - hitVPeakTime);
 
       if (deltaPeakTime > pairDeltaTimeLimits) continue;
 
       if (deltaPeakTime < -pairDeltaTimeLimits) break;
 
       pairDeltaTimeLimits = fabs(deltaPeakTime);
       bestVHit = hit2D;
     }
 
     return bestVHit;
   }

int lar_cluster3d::StandardHit3DBuilder::findGoodHitPairs	(	const reco::ClusterHit2D *	goldenHit,
		HitVector::iterator &	startItr,
		HitVector::iterator &	endItr,
		HitMatchPairVecMap &	hitMatchMap
	)		const

private

Definition at line 768 of file StandardHit3DBuilder_tool.cc.

   {
     int numPairs(0);
 
     // Loop through the input secon hits and make pairs
     while (startItr != endItr) {
       const reco::ClusterHit2D* hit = *startItr++;
       reco::ClusterHit3D pair;
 
       // pair returned with a negative ave time is signal of failure
       if (!makeHitPair(pair, goldenHit, hit, m_hitWidthSclFctr)) continue;
 
       geo::WireID wireID = hit->WireID();
 
       hitMatchMap[wireID].emplace_back(hit, pair);
 
       numPairs++;
     }
 
     return numPairs;
   }

void lar_cluster3d::StandardHit3DBuilder::findGoodTriplets	(	HitMatchPairVecMap &	pair12Map,
		HitMatchPairVecMap &	pair13Map,
		reco::HitPairList &	hitPairList,
		bool	tagged = `false`
	)		const

private

This algorithm takes lists of hit pairs and finds good triplets.

Definition at line 794 of file StandardHit3DBuilder_tool.cc.

   {
     // Build triplets from the two lists of hit pairs
     if (!pair12Map.empty()) {
       // temporary container for dead channel hits
       std::vector<reco::ClusterHit3D> tempDeadChanVec;
       reco::ClusterHit3D deadChanPair;
 
       // Keep track of which third plane hits have been used
       std::map<const reco::ClusterHit3D*, bool> usedPairMap;
 
       // Initial population of this map with the pair13Map hits
       for (const auto& pair13 : pair13Map) {
         for (const auto& hit2Dhit3DPair : pair13.second)
           usedPairMap[&hit2Dhit3DPair.second] = false;
       }
 
       // The outer loop is over all hit pairs made from the first two plane combinations
       for (const auto& pair12 : pair12Map) {
         if (pair12.second.empty()) continue;
 
         // This loop is over hit pairs that share the same first two plane wires but may have different
         // hit times on those wires
         for (const auto& hit2Dhit3DPair12 : pair12.second) {
           const reco::ClusterHit3D& pair1 = hit2Dhit3DPair12.second;
 
           // populate the map with initial value
           usedPairMap[&pair1] = false;
 
           // The simplest approach here is to loop over all possibilities and let the triplet builder weed out the weak candidates
           for (const auto& pair13 : pair13Map) {
             if (pair13.second.empty()) continue;
 
             for (const auto& hit2Dhit3DPair13 : pair13.second) {
               const reco::ClusterHit2D* hit2 = hit2Dhit3DPair13.first;
               const reco::ClusterHit3D& pair2 = hit2Dhit3DPair13.second;
 
               // If success try for the triplet
               reco::ClusterHit3D triplet;
 
               if (makeHitTriplet(triplet, pair1, hit2)) {
                 triplet.setID(hitPairList.size());
                 hitPairList.emplace_back(triplet);
                 usedPairMap[&pair1] = true;
                 usedPairMap[&pair2] = true;
               }
             }
           }
         }
       }
 
       // One more loop through the other pairs to check for sick channels
       if (m_numBadChannels > 0) {
         for (const auto& pairMapPair : usedPairMap) {
           if (pairMapPair.second) continue;
 
           const reco::ClusterHit3D* pair = pairMapPair.first;
 
           // Here we look to see if we failed to make a triplet because the partner wire was dead/noisy/sick
           if (makeDeadChannelPair(deadChanPair, *pair, 4, 0, 0.))
             tempDeadChanVec.emplace_back(deadChanPair);
         }
 
         // Handle the dead wire triplets
         if (!tempDeadChanVec.empty()) {
           // If we have many then see if we can trim down a bit by keeping those with time significance
           if (tempDeadChanVec.size() > 1) {
             // Sort by "significance" of agreement
             std::sort(tempDeadChanVec.begin(),
                       tempDeadChanVec.end(),
                       [](const auto& left, const auto& right) {
                         return left.getDeltaPeakTime() / left.getSigmaPeakTime() <
                                right.getDeltaPeakTime() / right.getSigmaPeakTime();
                       });
 
             // What is the range of "significance" from first to last?
             float firstSig = tempDeadChanVec.front().getDeltaPeakTime() /
                              tempDeadChanVec.front().getSigmaPeakTime();
             float lastSig =
               tempDeadChanVec.back().getDeltaPeakTime() / tempDeadChanVec.back().getSigmaPeakTime();
             float sigRange = lastSig - firstSig;
 
             if (lastSig > 0.5 * m_deltaPeakTimeSig && sigRange > 0.5) {
               // Declare a maximum of 1.5 * the average of the first and last pairs...
               float maxSignificance = std::max(0.75 * (firstSig + lastSig), 1.0);
 
               std::vector<reco::ClusterHit3D>::iterator firstBadElem = std::find_if(
                 tempDeadChanVec.begin(),
                 tempDeadChanVec.end(),
                 [&maxSignificance](const auto& pair) {
                   return pair.getDeltaPeakTime() / pair.getSigmaPeakTime() > maxSignificance;
                 });
 
               // But only keep the best 10?
               if (std::distance(tempDeadChanVec.begin(), firstBadElem) > 20)
                 firstBadElem = tempDeadChanVec.begin() + 20;
               // Keep at least one hit...
               else if (firstBadElem == tempDeadChanVec.begin())
                 firstBadElem++;
 
               tempDeadChanVec.resize(std::distance(tempDeadChanVec.begin(), firstBadElem));
             }
           }
 
           for (auto& pair : tempDeadChanVec) {
             pair.setID(hitPairList.size());
             hitPairList.emplace_back(pair);
           }
         }
       }
     }
 
     return;
   }

int lar_cluster3d::StandardHit3DBuilder::FindNumberInRange	(	const Hit2DSet &	hit2DSet,
		const reco::ClusterHit3D &	pair,
		float	range
	)		const

private

A utility routine for returning the number of 2D hits from the list in a given range.

Definition at line 1427 of file StandardHit3DBuilder_tool.cc.

   {
     static const float minCharge(0.);
 
     int numberInRange(0);
     float pairAvePeakTime(pair.getAvePeakTime());
 
     // Idea is to loop through the input set of hits and look for the best combination
     for (const auto& hit2D : hit2DSet) {
       if (hit2D->getHit()->Integral() < minCharge) continue;
 
       float hitVPeakTime(hit2D->getTimeTicks());
       float deltaPeakTime(pairAvePeakTime - hitVPeakTime);
 
       if (deltaPeakTime > range) continue;
 
       if (deltaPeakTime < -range) break;
 
       numberInRange++;
     }
 
     return numberInRange;
   }

float lar_cluster3d::StandardHit3DBuilder::getTimeToExecute ( IHit3DBuilder::TimeValues index ) const

inlineoverridevirtual

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

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 106 of file StandardHit3DBuilder_tool.cc.

     {
       return m_timeVector[index];
     }

void lar_cluster3d::StandardHit3DBuilder::Hit3DBuilder	(	art::Event &	evt,
		reco::HitPairList &	hitPairList,
		RecobHitToPtrMap &	clusterHitToArtPtrMap
	)

overridevirtual

Given a set of recob hits, run DBscan to form 3D clusters.

Parameters

hitPairList	The input list of 3D hits to run clustering on
clusterParametersList	A list of cluster objects (parameters from associated hits)

Associations with wires.

Associations with raw digits.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 455 of file StandardHit3DBuilder_tool.cc.

   {
     // Clear the internal data structures
     m_clusterHit2DMasterList.clear();
     m_planeToHitVectorMap.clear();
     m_planeToWireToHitSetMap.clear();
 
     m_timeVector.resize(NUMTIMEVALUES, 0.);
 
     // Get a hit refiner
     std::unique_ptr<std::vector<recob::Hit>> outputHitPtrVec(new std::vector<recob::Hit>);
 
     // Recover the 2D hits and then organize them into data structures which will be used in the
     // DBscan algorithm for building the 3D clusters
     this->CollectArtHits(evt);
 
     // If there are no hits in our view/wire data structure then do not proceed with the full analysis
     if (!m_planeToWireToHitSetMap.empty()) {
       // Call the algorithm that builds 3D hits
       this->BuildHit3D(hitPairList);
 
       // If we built 3D points then attempt to output a new hit list as well
       if (!hitPairList.empty())
         CreateNewRecobHitCollection(evt, hitPairList, *outputHitPtrVec, clusterHitToArtPtrMap);
     }
 
     // Set up to make the associations (if desired)
     /// Associations with wires.
     std::unique_ptr<art::Assns<recob::Wire, recob::Hit>> wireAssns(
       new art::Assns<recob::Wire, recob::Hit>);
 
     makeWireAssns(evt, *wireAssns, clusterHitToArtPtrMap);
 
     /// Associations with raw digits.
     std::unique_ptr<art::Assns<raw::RawDigit, recob::Hit>> rawDigitAssns(
       new art::Assns<raw::RawDigit, recob::Hit>);
 
     makeRawDigitAssns(evt, *rawDigitAssns, clusterHitToArtPtrMap);
 
     // Move everything into the event
     evt.put(std::move(outputHitPtrVec));
     evt.put(std::move(wireAssns));
     evt.put(std::move(rawDigitAssns));
 
     // Handle tree output too
     if (m_outputHistograms) {
       m_tupleTree->Fill();
 
       clear();
     }
 
     return;
   }

bool lar_cluster3d::StandardHit3DBuilder::makeDeadChannelPair	(	reco::ClusterHit3D &	pairOut,
		const reco::ClusterHit3D &	pair,
		size_t	maxStatus = `4`,
		size_t	minStatus = `0`,
		float	minOverlap = `0.2`
	)		const

private

Make a 3D HitPair object from a valid pair and a dead channel in the missing plane.

Definition at line 1319 of file StandardHit3DBuilder_tool.cc.

   {
     // Assume failure (most common result)
     bool result(false);
 
     const reco::ClusterHit2D* hit0 = pair.getHits()[0];
     const reco::ClusterHit2D* hit1 = pair.getHits()[1];
 
     size_t missPlane(2);
 
     // u plane hit is missing
     if (!hit0) {
       hit0 = pair.getHits()[2];
       missPlane = 0;
     }
     // v plane hit is missing
     else if (!hit1) {
       hit1 = pair.getHits()[2];
       missPlane = 1;
     }
 
     // Which plane is missing?
     geo::WireID wireID0 = hit0->WireID();
     geo::WireID wireID1 = hit1->WireID();
 
     // Ok, recover the wireID expected in the third plane...
     geo::WireID wireIn(wireID0.Cryostat, wireID0.TPC, missPlane, 0);
     geo::WireID wireID = NearestWireID(pair.getPosition(), wireIn);
 
     // There can be a round off issue so check the next wire as well
     bool wireStatus = m_channelStatus[wireID.Plane][wireID.Wire] < maxChanStatus &&
                       m_channelStatus[wireID.Plane][wireID.Wire] >= minChanStatus;
     bool wireOneStatus = m_channelStatus[wireID.Plane][wireID.Wire + 1] < maxChanStatus &&
                          m_channelStatus[wireID.Plane][wireID.Wire + 1] >= minChanStatus;
 
     // Make sure they are of at least the minimum status
     if (wireStatus || wireOneStatus) {
       // Sort out which is the wire we're dealing with
       if (!wireStatus) wireID.Wire += 1;
 
       // Want to refine position since we "know" the missing wire
       geo::WireIDIntersection widIntersect0;
 
       if (m_geometry->WireIDsIntersect(wireID0, wireID, widIntersect0)) {
         geo::WireIDIntersection widIntersect1;
 
         if (m_geometry->WireIDsIntersect(wireID1, wireID, widIntersect1)) {
           Eigen::Vector3f newPosition(
             pair.getPosition()[0], pair.getPosition()[1], pair.getPosition()[2]);
 
           newPosition[1] = (newPosition[1] + widIntersect0.y + widIntersect1.y) / 3.;
           newPosition[2] =
             (newPosition[2] + widIntersect0.z + widIntersect1.z - 2. * m_zPosOffset) / 3.;
 
           pairOut = pair;
           pairOut.setWireID(wireID);
           pairOut.setPosition(newPosition);
 
           if (hit0->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET)
             hit0->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
           if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET)
             hit1->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
 
           hit0->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
           hit1->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
 
           result = true;
         }
       }
     }
 
     return result;
   }

bool lar_cluster3d::StandardHit3DBuilder::makeHitPair	(	reco::ClusterHit3D &	pairOut,
		const reco::ClusterHit2D *	hit1,
		const reco::ClusterHit2D *	hit2,
		float	hitWidthSclFctr = `1.`,
		size_t	hitPairCntr = `0`
	)		const

private

Make a HitPair object by checking two hits.

Definition at line 913 of file StandardHit3DBuilder_tool.cc.

   {
     // Assume failure
     bool result(false);
 
     // We assume in this routine that we are looking at hits in different views
     // The first mission is to check that the wires intersect
     const geo::WireID& hit1WireID = hit1->WireID();
     const geo::WireID& hit2WireID = hit2->WireID();
 
     geo::WireIDIntersection widIntersect;
 
     if (m_geometry->WireIDsIntersect(hit1WireID, hit2WireID, widIntersect)) {
       // Wires intersect so now we can check the timing
       float hit1Peak = hit1->getTimeTicks();
       float hit1Sigma = hit1->getHit()->RMS();
 
       float hit2Peak = hit2->getTimeTicks();
       float hit2Sigma = hit2->getHit()->RMS();
 
       // "Long hits" are an issue... so we deal with these differently
       int hit1NDF = hit1->getHit()->DegreesOfFreedom();
       int hit2NDF = hit2->getHit()->DegreesOfFreedom();
 
       // Basically, allow the range to extend to the nearest end of the snippet
       if (hit1NDF < 2)
         hit1Sigma = std::min(hit1Peak - float(hit1->getHit()->StartTick()),
                              float(hit1->getHit()->EndTick()) - hit1Peak);
       if (hit2NDF < 2)
         hit2Sigma = std::min(hit2Peak - float(hit2->getHit()->StartTick()),
                              float(hit2->getHit()->EndTick()) - hit2Peak);
 
       // The "hit sigma" is the gaussian fit sigma of the hit, we need to expand a bit to allow hit overlap efficiency
       float hit1Width = hitWidthSclFctr * hit1Sigma;
       float hit2Width = hitWidthSclFctr * hit2Sigma;
 
       // Coarse check hit times are "in range"
       if (fabs(hit1Peak - hit2Peak) <= (hit1Width + hit2Width)) {
         // Check to see that hit peak times are consistent with each other
         float hit1SigSq = hit1Sigma * hit1Sigma;
         float hit2SigSq = hit2Sigma * hit2Sigma;
         float deltaPeakTime = std::fabs(hit1Peak - hit2Peak);
         float sigmaPeakTime = std::sqrt(hit1SigSq + hit2SigSq);
 
         // delta peak time consistency check here
         if (deltaPeakTime < m_deltaPeakTimeSig *
                               sigmaPeakTime) // 2 sigma consistency? (do this way to avoid divide)
         {
           float oneOverWghts = hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
           float avePeakTime = (hit1Peak / hit1SigSq + hit2Peak / hit2SigSq) * oneOverWghts;
           float totalCharge = hit1->getHit()->Integral() + hit2->getHit()->Integral();
           float hitChiSquare = std::pow((hit1Peak - avePeakTime), 2) / hit1SigSq +
                                std::pow((hit2Peak - avePeakTime), 2) / hit2SigSq;
 
           float xPositionHit1(hit1->getXPosition());
           float xPositionHit2(hit2->getXPosition());
           float xPosition = (xPositionHit1 / hit1SigSq + xPositionHit2 / hit2SigSq) * hit1SigSq *
                             hit2SigSq / (hit1SigSq + hit2SigSq);
 
           Eigen::Vector3f position(
             xPosition, float(widIntersect.y), float(widIntersect.z) - m_zPosOffset);
 
           // If to here then we need to sort out the hit pair code telling what views are used
           unsigned statusBits = 1 << hit1->WireID().Plane | 1 << hit2->WireID().Plane;
 
           // handle status bits for the 2D hits
           if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINPAIR)
             hit1->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);
           if (hit2->getStatusBits() & reco::ClusterHit2D::USEDINPAIR)
             hit2->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);
 
           hit1->setStatusBit(reco::ClusterHit2D::USEDINPAIR);
           hit2->setStatusBit(reco::ClusterHit2D::USEDINPAIR);
 
           reco::ClusterHit2DVec hitVector;
 
           hitVector.resize(3, NULL);
 
           hitVector[hit1->WireID().Plane] = hit1;
           hitVector[hit2->WireID().Plane] = hit2;
 
           unsigned int cryostatIdx = hit1->WireID().Cryostat;
           unsigned int tpcIdx = hit1->WireID().TPC;
 
           // Initialize the wireIdVec
           std::vector<geo::WireID> wireIDVec = {geo::WireID(cryostatIdx, tpcIdx, 0, 0),
                                                 geo::WireID(cryostatIdx, tpcIdx, 1, 0),
                                                 geo::WireID(cryostatIdx, tpcIdx, 2, 0)};
 
           wireIDVec[hit1->WireID().Plane] = hit1->WireID();
           wireIDVec[hit2->WireID().Plane] = hit2->WireID();
 
           // For compiling at the moment
           std::vector<float> hitDelTSigVec = {0., 0., 0.};
 
           hitDelTSigVec[hit1->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
           hitDelTSigVec[hit2->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
 
           // Create the 3D cluster hit
           hitPair.initialize(hitPairCntr,
                              statusBits,
                              position,
                              totalCharge,
                              avePeakTime,
                              deltaPeakTime,
                              sigmaPeakTime,
                              hitChiSquare,
                              0.,
                              0.,
                              0.,
                              0.,
                              hitVector,
                              hitDelTSigVec,
                              wireIDVec);
 
           result = true;
         }
       }
     }
 
     // Send it back
     return result;
   }

bool lar_cluster3d::StandardHit3DBuilder::makeHitTriplet	(	reco::ClusterHit3D &	pairOut,
		const reco::ClusterHit3D &	pairIn,
		const reco::ClusterHit2D *	hit2
	)		const

private

Make a 3D HitPair object by checking two hits.

Definition at line 1042 of file StandardHit3DBuilder_tool.cc.

   {
     // Assume failure
     bool result(false);
 
     // We are going to force the wire pitch here, some time in the future we need to fix
     static const double wirePitch =
       0.5 * m_wirePitchScaleFactor * *std::max_element(m_wirePitch, m_wirePitch + 3);
 
     // Recover hit info
     float hitTimeTicks = hit->getTimeTicks();
     float hitSigma = hit->getHit()->RMS();
 
     // Special case check
     if (hitSigma > 2. * hit->getHit()->PeakAmplitude())
       hitSigma = 2. * hit->getHit()->PeakAmplitude();
 
     // Let's do a quick consistency check on the input hits to make sure we are in range...
     // Require the W hit to be "in range" with the UV Pair
     if (fabs(hitTimeTicks - pair.getAvePeakTime()) <
         m_hitWidthSclFctr * (pair.getSigmaPeakTime() + hitSigma)) {
       // Check the distance from the point to the wire the hit is on
       float hitWireDist = DistanceFromPointToHitWire(pair.getPosition(), hit->WireID());
 
       if (m_outputHistograms) m_maxSideVecVec.push_back(hitWireDist);
 
       // Reject hits that are not within range
       if (hitWireDist < wirePitch) {
         if (m_outputHistograms) m_pairWireDistVec.push_back(hitWireDist);
 
         // Use the existing code to see the U and W hits are willing to pair with the V hit
         reco::ClusterHit3D pair0h;
         reco::ClusterHit3D pair1h;
 
         // Recover all the hits involved
         const reco::ClusterHit2DVec& pairHitVec = pair.getHits();
         const reco::ClusterHit2D* hit0 = pairHitVec[0];
         const reco::ClusterHit2D* hit1 = pairHitVec[1];
 
         if (!hit0)
           hit0 = pairHitVec[2];
         else if (!hit1)
           hit1 = pairHitVec[2];
 
         // If good pairs made here then we can try to make a triplet
         if (makeHitPair(pair0h, hit0, hit, m_hitWidthSclFctr) &&
             makeHitPair(pair1h, hit1, hit, m_hitWidthSclFctr)) {
           // Get a copy of the input hit vector (note the order is by plane - by definition)
           reco::ClusterHit2DVec hitVector = pair.getHits();
 
           // include the new hit
           hitVector[hit->WireID().Plane] = hit;
 
           // Set up to get average peak time, hitChiSquare, etc.
           unsigned int statusBits(0x7);
           float avePeakTime(0.);
           float weightSum(0.);
           float xPosition(0.);
 
           // And get the wire IDs
           std::vector<geo::WireID> wireIDVec = {geo::WireID(), geo::WireID(), geo::WireID()};
 
           // First loop through the hits to get WireIDs and calculate the averages
           for (size_t planeIdx = 0; planeIdx < 3; planeIdx++) {
             const reco::ClusterHit2D* hit2D = hitVector[planeIdx];
 
             wireIDVec[planeIdx] = hit2D->WireID();
 
             if (hit2D->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET)
               hit2D->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
 
             hit2D->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
 
             float hitRMS = hit2D->getHit()->RMS();
             float peakTime = hit2D->getTimeTicks();
 
             // Basically, allow the range to extend to the nearest end of the snippet
             if (hit2D->getHit()->DegreesOfFreedom() < 2)
               hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                                 float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());
 
             float weight = 1. / (hitRMS * hitRMS);
 
             avePeakTime += peakTime * weight;
             xPosition += hit2D->getXPosition() * weight;
             weightSum += weight;
           }
 
           avePeakTime /= weightSum;
           xPosition /= weightSum;
 
           Eigen::Vector2f pair0hYZVec(pair0h.getPosition()[1], pair0h.getPosition()[2]);
           Eigen::Vector2f pair1hYZVec(pair1h.getPosition()[1], pair1h.getPosition()[2]);
           Eigen::Vector2f pairYZVec(pair.getPosition()[1], pair.getPosition()[2]);
           Eigen::Vector3f position(xPosition,
                                    float((pairYZVec[0] + pair0hYZVec[0] + pair1hYZVec[0]) / 3.),
                                    float((pairYZVec[1] + pair0hYZVec[1] + pair1hYZVec[1]) / 3.));
 
           // Armed with the average peak time, now get hitChiSquare and the sig vec
           float hitChiSquare(0.);
           float sigmaPeakTime(std::sqrt(1. / weightSum));
           std::vector<float> hitDelTSigVec;
 
           for (const auto& hit2D : hitVector) {
             float hitRMS = hit2D->getHit()->RMS();
 
             // Basically, allow the range to extend to the nearest end of the snippet
             if (hit2D->getHit()->DegreesOfFreedom() < 2)
               hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                                 float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());
 
             float combRMS = std::sqrt(hitRMS * hitRMS - sigmaPeakTime * sigmaPeakTime);
             float peakTime = hit2D->getTimeTicks();
             float deltaTime = peakTime - avePeakTime;
             float hitSig = deltaTime / combRMS;
 
             hitChiSquare += hitSig * hitSig;
 
             hitDelTSigVec.emplace_back(std::fabs(hitSig));
           }
 
           if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare);
 
           int lowMinIndex(std::numeric_limits<int>::max());
           int lowMaxIndex(std::numeric_limits<int>::min());
           int hiMinIndex(std::numeric_limits<int>::max());
           int hiMaxIndex(std::numeric_limits<int>::min());
 
           // First task is to get the min/max values for the common overlap region
           for (const auto& hit2D : hitVector) {
             float range = 2. * hit2D->getHit()->RMS();
 
             // Basically, allow the range to extend to the nearest end of the snippet
             if (hit2D->getHit()->DegreesOfFreedom() < 2)
               range = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                                float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());
 
             int hitStart = hit2D->getHit()->PeakTime() - range - 0.5;
             int hitStop = hit2D->getHit()->PeakTime() + range + 0.5;
 
             lowMinIndex = std::min(hitStart, lowMinIndex);
             lowMaxIndex = std::max(hitStart, lowMaxIndex);
             hiMinIndex = std::min(hitStop + 1, hiMinIndex);
             hiMaxIndex = std::max(hitStop + 1, hiMaxIndex);
           }
 
           // Keep only "good" hits...
           if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex) {
             // One more pass through hits to get charge
             std::vector<float> chargeVec;
 
             for (const auto& hit2D : hitVector)
               chargeVec.push_back(chargeIntegral(hit2D->getHit()->PeakTime(),
                                                  hit2D->getHit()->PeakAmplitude(),
                                                  hit2D->getHit()->RMS(),
                                                  1.,
                                                  lowMaxIndex,
                                                  hiMinIndex));
 
             float totalCharge =
               std::accumulate(chargeVec.begin(), chargeVec.end(), 0.) / float(chargeVec.size());
             float overlapRange = float(hiMinIndex - lowMaxIndex);
             float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex);
 
             // Set up to compute the charge asymmetry
             std::vector<float> smallestChargeDiffVec;
             std::vector<float> chargeAveVec;
             float smallestDiff(std::numeric_limits<float>::max());
             float maxDeltaPeak(0.);
             size_t chargeIndex(0);
 
             for (size_t idx = 0; idx < 3; idx++) {
               size_t leftIdx = (idx + 2) % 3;
               size_t rightIdx = (idx + 1) % 3;
 
               smallestChargeDiffVec.push_back(std::abs(chargeVec[leftIdx] - chargeVec[rightIdx]));
               chargeAveVec.push_back(float(0.5 * (chargeVec[leftIdx] + chargeVec[rightIdx])));
 
               if (smallestChargeDiffVec.back() < smallestDiff) {
                 smallestDiff = smallestChargeDiffVec.back();
                 chargeIndex = idx;
               }
 
               // Take opportunity to look at peak time diff
               float deltaPeakTime =
                 hitVector[leftIdx]->getTimeTicks() - hitVector[rightIdx]->getTimeTicks();
 
               if (std::abs(deltaPeakTime) > maxDeltaPeak) maxDeltaPeak = std::abs(deltaPeakTime);
 
               if (m_outputHistograms) m_deltaTimeVec.push_back(deltaPeakTime);
             }
 
             float chargeAsymmetry = (chargeAveVec[chargeIndex] - chargeVec[chargeIndex]) /
                                     (chargeAveVec[chargeIndex] + chargeVec[chargeIndex]);
 
             // If this is true there has to be a negative charge that snuck in somehow
             if (chargeAsymmetry < -1. || chargeAsymmetry > 1.) {
               const geo::WireID& hitWireID = hitVector[chargeIndex]->WireID();
 
               std::cout << "============> Charge asymmetry out of range: " << chargeAsymmetry
                         << " <============" << std::endl;
               std::cout << "     hit C: " << hitWireID.Cryostat << ", TPC: " << hitWireID.TPC
                         << ", Plane: " << hitWireID.Plane << ", Wire: " << hitWireID.Wire
                         << std::endl;
               std::cout << "     charge: " << chargeVec[0] << ", " << chargeVec[1] << ", "
                         << chargeVec[2] << std::endl;
               std::cout << "     index: " << chargeIndex << ", smallest diff: " << smallestDiff
                         << std::endl;
               return result;
             }
 
             // Usurping "deltaPeakTime" to be the maximum pull
             float deltaPeakTime = *std::max_element(hitDelTSigVec.begin(), hitDelTSigVec.end());
 
             if (m_outputHistograms) {
               m_smallChargeDiffVec.push_back(smallestDiff);
               m_smallIndexVec.push_back(chargeIndex);
               m_maxPullVec.push_back(deltaPeakTime);
               m_qualityMetricVec.push_back(hitChiSquare);
               m_spacePointChargeVec.push_back(totalCharge);
               m_overlapFractionVec.push_back(overlapFraction);
               m_overlapRangeVec.push_back(overlapRange);
               m_maxDeltaPeakVec.push_back(maxDeltaPeak);
               m_hitAsymmetryVec.push_back(chargeAsymmetry);
             }
 
             // Try to weed out cases where overlap doesn't match peak separation
             if (maxDeltaPeak > overlapRange) return result;
 
             // Create the 3D cluster hit
             hitTriplet.initialize(0,
                                   statusBits,
                                   position,
                                   totalCharge,
                                   avePeakTime,
                                   deltaPeakTime,
                                   sigmaPeakTime,
                                   hitChiSquare,
                                   overlapFraction,
                                   chargeAsymmetry,
                                   0.,
                                   0.,
                                   hitVector,
                                   hitDelTSigVec,
                                   wireIDVec);
 
             result = true;
           }
         }
       }
     }
 
     // return success/fail
     return result;
   }

void lar_cluster3d::StandardHit3DBuilder::makeRawDigitAssns	(	const art::Event &	evt,
		art::Assns< raw::RawDigit, recob::Hit > &	rawDigitAssns,
		RecobHitToPtrMap &	recobHitPtrMap
	)		const

private

Create raw::RawDigit to recob::Hit associations.

Definition at line 1786 of file StandardHit3DBuilder_tool.cc.

   {
     // Let's make sure the input associations container is empty
     rawDigitAssns = art::Assns<raw::RawDigit, recob::Hit>();
 
     // First task is to recover all of the previous wire <--> hit associations and map them by channel number
     // Create the temporary container
     std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>> channelToRawDigitMap;
 
     // Go through the list of input sources and fill out the map
     for (const auto& inputTag : m_hitFinderTagVec) {
       art::ValidHandle<std::vector<recob::Hit>> hitHandle =
         evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
 
       art::FindOneP<raw::RawDigit> hitToRawDigitAssns(hitHandle, evt, inputTag);
 
       if (hitToRawDigitAssns.isValid()) {
         for (size_t rawDigitIdx = 0; rawDigitIdx < hitToRawDigitAssns.size(); rawDigitIdx++) {
           art::Ptr<raw::RawDigit> rawDigit = hitToRawDigitAssns.at(rawDigitIdx);
 
           channelToRawDigitMap[rawDigit->Channel()] = rawDigit;
         }
       }
     }
 
     // Now fill the container
     for (const auto& hitPtrPair : recobHitPtrMap) {
       raw::ChannelID_t channel = hitPtrPair.first->Channel();
 
       std::unordered_map<raw::ChannelID_t, art::Ptr<raw::RawDigit>>::iterator chanRawDigitItr =
         channelToRawDigitMap.find(channel);
 
       if (!(chanRawDigitItr != channelToRawDigitMap.end())) {
         mf::LogDebug("Cluster3D") << "** Did not find channel to wire match! Skipping..."
                                   << std::endl;
         continue;
       }
 
       rawDigitAssns.addSingle(chanRawDigitItr->second, hitPtrPair.second);
     }
 
     return;
   }

void lar_cluster3d::StandardHit3DBuilder::makeWireAssns	(	const art::Event &	evt,
		art::Assns< recob::Wire, recob::Hit > &	wireAssns,
		RecobHitToPtrMap &	recobHitPtrMap
	)		const

private

Create recob::Wire to recob::Hit associations.

Definition at line 1739 of file StandardHit3DBuilder_tool.cc.

   {
     // Let's make sure the input associations container is empty
     wireAssns = art::Assns<recob::Wire, recob::Hit>();
 
     // First task is to recover all of the previous wire <--> hit associations and map them by channel number
     // Create the temporary container
     std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>> channelToWireMap;
 
     // Go through the list of input sources and fill out the map
     for (const auto& inputTag : m_hitFinderTagVec) {
       art::ValidHandle<std::vector<recob::Hit>> hitHandle =
         evt.getValidHandle<std::vector<recob::Hit>>(inputTag);
 
       art::FindOneP<recob::Wire> hitToWireAssns(hitHandle, evt, inputTag);
 
       if (hitToWireAssns.isValid()) {
         for (size_t wireIdx = 0; wireIdx < hitToWireAssns.size(); wireIdx++) {
           art::Ptr<recob::Wire> wire = hitToWireAssns.at(wireIdx);
 
           channelToWireMap[wire->Channel()] = wire;
         }
       }
     }
 
     // Now fill the container
     for (const auto& hitPtrPair : recobHitPtrMap) {
       raw::ChannelID_t channel = hitPtrPair.first->Channel();
 
       std::unordered_map<raw::ChannelID_t, art::Ptr<recob::Wire>>::iterator chanWireItr =
         channelToWireMap.find(channel);
 
       if (!(chanWireItr != channelToWireMap.end())) {
         mf::LogDebug("Cluster3D") << "** Did not find channel to wire match! Skipping..."
                                   << std::endl;
         continue;
       }
 
       wireAssns.addSingle(chanWireItr->second, hitPtrPair.second);
     }
 
     return;
   }

geo::WireID lar_cluster3d::StandardHit3DBuilder::NearestWireID	(	const Eigen::Vector3f &	position,
		const geo::WireID &	wireID
	)		const

private

Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range.

Definition at line 1454 of file StandardHit3DBuilder_tool.cc.

   {
     geo::WireID wireID = wireIDIn;
 
     // Embed the call to the geometry's services nearest wire id method in a try-catch block
     try {
       // Switch from NearestWireID to this method to avoid the roundoff error issues...
       double distanceToWire = m_geometry->Plane(wireIDIn).WireCoordinate(position.data());
 
       wireID.Wire = int(distanceToWire);
     }
     catch (std::exception& exc) {
       // This can happen, almost always because the coordinates are **just** out of range
       mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - "
                                   << exc.what() << std::endl;
 
       // Assume extremum for wire number depending on z coordinate
       if (position[2] < 0.5 * m_geometry->DetLength())
         wireID.Wire = 0;
       else
         wireID.Wire = m_geometry->Nwires(wireIDIn.Plane) - 1;
     }
 
     return wireID;
   }

void lar_cluster3d::StandardHit3DBuilder::produces ( art::ProducesCollector & collector )

overridevirtual

Each algorithm may have different objects it wants "produced" so use this to let the top level producer module "know" what it is outputting.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 312 of file StandardHit3DBuilder_tool.cc.

   {
     collector.produces<std::vector<recob::Hit>>();
     collector.produces<art::Assns<recob::Wire, recob::Hit>>();
     collector.produces<art::Assns<raw::RawDigit, recob::Hit>>();
   }