lar_cluster3d::SnippetHit3DBuilder Class Reference

SnippetHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::SnippetHit3DBuilder:

Public Member Functions
	SnippetHit3DBuilder (fhicl::ParameterSet const &pset)
	Constructor. More...
virtual 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...
virtual void	configure (const fhicl::ParameterSet &) override
	Interface for configuring the particular algorithm tool. More...
virtual void	Hit3DBuilder (art::Event &, reco::HitPairList &, RecobHitToPtrMap &) override
	Given a set of recob hits, run DBscan to form 3D clusters. More...
virtual 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	HitMatchTriplet = std::tuple< const reco::ClusterHit2D *, const reco::ClusterHit2D *, const reco::ClusterHit3D >
	This builds a list of candidate hit pairs from lists of hits on two planes. More...
using	HitMatchTripletVec = std::vector< HitMatchTriplet >
using	HitMatchTripletVecMap = std::map< geo::WireID, HitMatchTripletVec >
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 (PlaneToSnippetHitMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...
size_t	BuildHitPairMapByTPC (PlaneSnippetHitMapItrPairVec &planeSnippetHitMapItrPairVec, reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...
int	findGoodHitPairs (SnippetHitMap::iterator &, SnippetHitMap::iterator &, SnippetHitMap::iterator &, HitMatchTripletVecMap &) const
void	findGoodTriplets (HitMatchTripletVecMap &, HitMatchTripletVecMap &, 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...
bool	WireIDsIntersect (const geo::WireID &, const geo::WireID &, geo::WireIDIntersection &) const
	function to detemine if two wires "intersect" (in the 2D sense) More...
float	closestApproach (const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, float &, float &) const
	function to compute the distance of closest approach and the arc length to the points of closest approach 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_hitWidthSclFctr
float	m_deltaPeakTimeSig
float	m_rangeNumSig
float	m_LongHitStretchFctr
float	m_pulseHeightFrac
float	m_PHLowSelection
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
PlaneToSnippetHitMap	m_planeToSnippetHitMap
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

SnippetHit3DBuilder class definiton.

Definition at line 81 of file SnippetHit3DBuilder_tool.cc.

Member Typedef Documentation

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

private

Definition at line 231 of file SnippetHit3DBuilder_tool.cc.

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

private

define data structure for keeping track of channel status

Definition at line 230 of file SnippetHit3DBuilder_tool.cc.

using lar_cluster3d::SnippetHit3DBuilder::HitMatchTriplet = std::tuple<const reco::ClusterHit2D*,const reco::ClusterHit2D*,const reco::ClusterHit3D>

private

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

Definition at line 154 of file SnippetHit3DBuilder_tool.cc.

using lar_cluster3d::SnippetHit3DBuilder::HitMatchTripletVec = std::vector<HitMatchTriplet>

private

Definition at line 155 of file SnippetHit3DBuilder_tool.cc.

using lar_cluster3d::SnippetHit3DBuilder::HitMatchTripletVecMap = std::map<geo::WireID,HitMatchTripletVec>

private

Definition at line 156 of file SnippetHit3DBuilder_tool.cc.

Constructor & Destructor Documentation

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

explicit

Constructor.

Parameters

pset

Definition at line 291 of file SnippetHit3DBuilder_tool.cc.

                                                                      :
    m_channelFilter(&art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider())

{
    this->configure(pset);
}

Member Function Documentation

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

private

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

Definition at line 379 of file SnippetHit3DBuilder_tool.cc.

{
    // This is called each event, clear out the previous version and start over
    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++;
        }
    }

    return;
}

void lar_cluster3d::SnippetHit3DBuilder::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 482 of file SnippetHit3DBuilder_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_planeToSnippetHitMap, 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::SnippetHit3DBuilder::BuildHitPairMap ( PlaneToSnippetHitMap &  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 538 of file SnippetHit3DBuilder_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++)
        {
            //************************************
            // Kludge
//            if (!(cryoIdx == 1 && tpcIdx == 0)) continue;

            PlaneToSnippetHitMap::iterator mapItr0 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,0));
            PlaneToSnippetHitMap::iterator mapItr1 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,1));
            PlaneToSnippetHitMap::iterator mapItr2 = planeToSnippetHitMap.find(geo::PlaneID(cryoIdx,tpcIdx,2));

            size_t nPlanesWithHits = (mapItr0 != planeToSnippetHitMap.end() && !mapItr0->second.empty() ? 1 : 0)
                                   + (mapItr1 != planeToSnippetHitMap.end() && !mapItr1->second.empty() ? 1 : 0)
                                   + (mapItr2 != planeToSnippetHitMap.end() && !mapItr2->second.empty() ? 1 : 0);

            if (nPlanesWithHits < 2) continue;

            SnippetHitMap& snippetHitMap0 = mapItr0->second;
            SnippetHitMap& snippetHitMap1 = mapItr1->second;
            SnippetHitMap& snippetHitMap2 = mapItr2->second;

            PlaneSnippetHitMapItrPairVec hitItrVec = {SnippetHitMapItrPair(snippetHitMap0.begin(),snippetHitMap0.end()),
                                                      SnippetHitMapItrPair(snippetHitMap1.begin(),snippetHitMap1.end()),
                                                      SnippetHitMapItrPair(snippetHitMap2.begin(),snippetHitMap2.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::SnippetHit3DBuilder::BuildHitPairMapByTPC ( PlaneSnippetHitMapItrPairVec &  planeSnippetHitMapItrPairVec, 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 597 of file SnippetHit3DBuilder_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 = [](SnippetHitMap::iterator startItr, SnippetHitMap::iterator endItr, float startTime)
    {
        while(startItr != endItr)
        {
            if (startItr->first.second < startTime) startItr++;
            else break;
        }
        return startItr;
    };

    auto SetEndIterator = [](SnippetHitMap::iterator lastItr, SnippetHitMap::iterator endItr, float endTime)
    {
        while(lastItr != endItr)
        {
            if (lastItr->first.first < endTime) lastItr++;
            else break;
        }
        return lastItr;
    };

    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(snippetHitMapItrVec.begin(),snippetHitMapItrVec.end(),SetStartTimeOrder());

        // Make sure there are still hits on at least
        int nPlanesWithHits(0);

        for(auto& pair : snippetHitMapItrVec)
            if (pair.first != pair.second) nPlanesWithHits++;

        if (nPlanesWithHits < 2) break;

        // End condition: no more hit snippets
//        if (snippetHitMapItrVec.front().first == snippetHitMapItrVec.front().second) break;

        // This loop iteration's snippet iterator
        SnippetHitMap::iterator firstSnippetItr = snippetHitMapItrVec.front().first;

        // Set iterators to insure we'll be in the overlap ranges
        SnippetHitMap::iterator snippetHitMapItr1Start = SetStartIterator(snippetHitMapItrVec[1].first, snippetHitMapItrVec[1].second, firstSnippetItr->first.first);
        SnippetHitMap::iterator snippetHitMapItr1End   = SetEndIterator( snippetHitMapItr1Start,        snippetHitMapItrVec[1].second, firstSnippetItr->first.second);
        SnippetHitMap::iterator snippetHitMapItr2Start = SetStartIterator(snippetHitMapItrVec[2].first, snippetHitMapItrVec[2].second, firstSnippetItr->first.first);
        SnippetHitMap::iterator snippetHitMapItr2End   = SetEndIterator( snippetHitMapItr2Start,        snippetHitMapItrVec[2].second, firstSnippetItr->first.second);

        // 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());
        HitMatchTripletVecMap pair12Map;
        HitMatchTripletVecMap pair13Map;

        size_t n12Pairs = findGoodHitPairs(firstSnippetItr, snippetHitMapItr1Start, snippetHitMapItr1End, pair12Map);
        size_t n13Pairs = findGoodHitPairs(firstSnippetItr, snippetHitMapItr2Start, snippetHitMapItr2End, pair13Map);

        nDeadChanHits  += hitPairList.size() - curHitListSize;
        curHitListSize  = hitPairList.size();

        if (n12Pairs > n13Pairs) findGoodTriplets(pair12Map, pair13Map, hitPairList);
        else                     findGoodTriplets(pair13Map, pair12Map, hitPairList);

        nTriplets += hitPairList.size() - curHitListSize;

        snippetHitMapItrVec.front().first++;
    }

    return hitPairList.size();
}

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

private

Perform charge integration between limits.

Definition at line 1302 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::clear ( )

private

clear the tuple vectors before processing next event

Definition at line 360 of file SnippetHit3DBuilder_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;
}

float lar_cluster3d::SnippetHit3DBuilder::closestApproach ( const Eigen::Vector3f &  P0, const Eigen::Vector3f &  u0, const Eigen::Vector3f &  P1, const Eigen::Vector3f &  u1, float &  arcLen0, float &  arcLen1  ) const

private

function to compute the distance of closest approach and the arc length to the points of closest approach

Definition at line 1277 of file SnippetHit3DBuilder_tool.cc.

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

    arcLen0 = (b * e - c * d) / den;
    arcLen1 = (a * e - b * d) / den;

    Eigen::Vector3f poca0 = P0 + arcLen0 * u0;
    Eigen::Vector3f poca1 = P1 + arcLen1 * u1;

    return (poca0 - poca1).norm();
}

void lar_cluster3d::SnippetHit3DBuilder::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 1535 of file SnippetHit3DBuilder_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;

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

    // Need the detector properties which needs the clocks
    auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const det_prop   = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);

    // 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_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,0)] = SnippetHitMap();
            m_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,1)] = SnippetHitMap();
            m_planeToSnippetHitMap[geo::PlaneID(cryoIdx,tpcIdx,2)] = SnippetHitMap();

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

        std::cout << 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());

        // Start/End ticks to identify the snippet
        HitStartEndPair hitStartEndPair(recobHit->StartTick(),recobHit->EndTick());

        // Can this really happen?
        if (hitStartEndPair.second <= hitStartEndPair.first)
        {
            std::cout << "Yes, found a hit with end time less than start time: " << hitStartEndPair.first << "/" << hitStartEndPair.second << ", mult: " << recobHit->Multiplicity() << std::endl;
            continue;
        }

        // And then loop over all possible to build out our maps
        //for(const auto& wireID : wireIDs)
        for(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_planeToSnippetHitMap[planeID][hitStartEndPair].emplace_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_planeToSnippetHitMap)
//        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::SnippetHit3DBuilder::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 309 of file SnippetHit3DBuilder_tool.cc.

{
    m_hitFinderTagVec      = pset.get<std::vector<art::InputTag>>("HitFinderTagVec",       std::vector<art::InputTag>()={"gaushit"});
    m_enableMonitoring     = pset.get<bool                      >("EnableMonitoring",      true);
    m_hitWidthSclFctr      = pset.get<float                     >("HitWidthScaleFactor",   6.  );
    m_rangeNumSig          = pset.get<float                     >("RangeNumSigma",         3.  );
    m_LongHitStretchFctr   = pset.get<float                     >("LongHitsStretchFactor", 1.5 );
    m_pulseHeightFrac      = pset.get<float                     >("PulseHeightFraction",   0.5 );
    m_PHLowSelection       = pset.get<float                     >("PHLowSelection",        20. );
    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.
    if (m_outputHistograms)
    {
        art::ServiceHandle<art::TFileService> tfs;

        m_tupleTree = tfs->make<TTree>("Hit3DBuilderTree", "Tree by SnippetHit3DBuilder");

        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::SnippetHit3DBuilder::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 1673 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::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 1474 of file SnippetHit3DBuilder_tool.cc.

{
    float distance = std::numeric_limits<float>::max();

    // Embed the call to the geometry's services nearest wire id method in a try-catch block
    try
    {
        // Recover wire geometry information for each wire
        const geo::WireGeo& wireGeo = m_geometry->WireIDToWireGeo(wireIDIn);

        // Get wire position and direction for first wire
        double wirePosArr[3] = {0.,0.,0.};
        wireGeo.GetCenter(wirePosArr);

        Eigen::Vector3f wirePos(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
        Eigen::Vector3f wireDir(wireGeo.Direction().X(),wireGeo.Direction().Y(),wireGeo.Direction().Z());

        //*********************************
        // Kludge
//        if (wireIDIn.Plane > 0) wireDir[2] = -wireDir[2];


        // Want the hit position to have same x value as wire coordinates
        Eigen::Vector3f hitPosition(wirePos[0],position[1],position[2]);

        // Get arc length to doca
        double arcLen = (hitPosition - wirePos).dot(wireDir);

        // Make sure arclen is in range
        if (abs(arcLen) < wireGeo.HalfL())
        {
            Eigen::Vector3f docaVec = hitPosition - (wirePos + 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::SnippetHit3DBuilder::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 1397 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::findGoodHitPairs ( SnippetHitMap::iterator &  firstSnippetItr, SnippetHitMap::iterator &  startItr, SnippetHitMap::iterator &  endItr, HitMatchTripletVecMap &  hitMatchMap  ) const

private

Definition at line 682 of file SnippetHit3DBuilder_tool.cc.

{
    int numPairs(0);

    HitVector::iterator firstMaxItr = std::max_element(firstSnippetItr->second.begin(),firstSnippetItr->second.end(),[](const auto& left, const auto& right){return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();});
    float               firstPHCut  = firstMaxItr != firstSnippetItr->second.end() ? m_pulseHeightFrac * (*firstMaxItr)->getHit()->PeakAmplitude() : 4096.;

    // Loop through the hits on the first snippet
    for(const auto& hit1 : firstSnippetItr->second)
    {
        // Let's focus on the largest hit in the chain
        if (hit1->getHit()->DegreesOfFreedom() > 1 && hit1->getHit()->PeakAmplitude() < firstPHCut && hit1->getHit()->PeakAmplitude() < m_PHLowSelection) continue;

        // Inside loop iterator
        SnippetHitMap::iterator secondHitItr = startItr;

        // Loop through the input secon hits and make pairs
        while(secondHitItr != endItr)
        {
            HitVector::iterator secondMaxItr = std::max_element(secondHitItr->second.begin(),secondHitItr->second.end(),[](const auto& left, const auto& right){return left->getHit()->PeakAmplitude() < right->getHit()->PeakAmplitude();});
            float               secondPHCut  = secondMaxItr != secondHitItr->second.end() ? m_pulseHeightFrac * (*secondMaxItr)->getHit()->PeakAmplitude() : 0.;

            for(const auto& hit2 : secondHitItr->second)
            {
                // Again, focus on the large hits
                if (hit2->getHit()->DegreesOfFreedom() > 1 && hit2->getHit()->PeakAmplitude() < secondPHCut && hit2->getHit()->PeakAmplitude() < m_PHLowSelection) continue;

                reco::ClusterHit3D  pair;

                // pair returned with a negative ave time is signal of failure
                if (!makeHitPair(pair, hit1, hit2, m_hitWidthSclFctr)) continue;

                std::vector<const recob::Hit*> recobHitVec = {nullptr,nullptr,nullptr};

                recobHitVec[hit1->WireID().Plane] = hit1->getHit();
                recobHitVec[hit2->WireID().Plane] = hit2->getHit();

                geo::WireID wireID = hit2->WireID();

                hitMatchMap[wireID].emplace_back(hit1,hit2,pair);

                numPairs++;
            }

            secondHitItr++;
        }
    }

    return numPairs;
}

void lar_cluster3d::SnippetHit3DBuilder::findGoodTriplets ( HitMatchTripletVecMap &  pair12Map, HitMatchTripletVecMap &  pair13Map, reco::HitPairList &  hitPairList, bool  tagged = false  ) const

private

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

Definition at line 736 of file SnippetHit3DBuilder_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[&std::get<2>(hit2Dhit3DPair)] = 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  = std::get<2>(hit2Dhit3DPair12);

                // 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)
                    {
                        // Protect against double counting
                        if (std::get<0>(hit2Dhit3DPair12) != std::get<0>(hit2Dhit3DPair13)) continue;

                        const reco::ClusterHit2D* hit2  = std::get<1>(hit2Dhit3DPair13);
                        const reco::ClusterHit3D& pair2 = std::get<2>(hit2Dhit3DPair13);

                        // 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::SnippetHit3DBuilder::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 1424 of file SnippetHit3DBuilder_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;
}

virtual float lar_cluster3d::SnippetHit3DBuilder::getTimeToExecute ( IHit3DBuilder::TimeValues  index ) const

inlineoverridevirtual

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

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 110 of file SnippetHit3DBuilder_tool.cc.

{return m_timeVector[index];}

void lar_cluster3d::SnippetHit3DBuilder::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 428 of file SnippetHit3DBuilder_tool.cc.

{
    // Clear the internal data structures
    m_clusterHit2DMasterList.clear();
    m_planeToSnippetHitMap.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::SnippetHit3DBuilder::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 1320 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::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 851 of file SnippetHit3DBuilder_tool.cc.

{
    // Assume failure
    bool result(false);

    // Start by checking time consistency since this is fastest
    // 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 *= m_LongHitStretchFctr; // This sets the range to the width of the pulse
    if (hit2NDF < 2) hit2Sigma *= m_LongHitStretchFctr;

    // 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)
        {
            // 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 (WireIDsIntersect(hit1WireID, hit2WireID, widIntersect))
            {
                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;
            }
        }
//        else std::cout << "-MakeHitPair, deltaPeakTime: " << deltaPeakTime << ", scl fctr: " << m_deltaPeakTimeSig << ", sigmaPeakTime: " << sigmaPeakTime << std::endl;
    }
//    else std::cout << "-MakeHitPair, delta peak: " << hit1Peak - hit2Peak << ", hit1Width: " << hit1Width << ", hit2Width: " << hit2Width << std::endl;

    // Send it back
    return result;
}

bool lar_cluster3d::SnippetHit3DBuilder::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 976 of file SnippetHit3DBuilder_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 *= m_LongHitStretchFctr;
                        //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 *= m_LongHitStretchFctr;
                        //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 = m_rangeNumSig * hit2D->getHit()->RMS();

                    // Basically, allow the range to extend to the nearest end of the snippet
                    if (hit2D->getHit()->DegreesOfFreedom() < 2) range *= m_LongHitStretchFctr;
                        //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 > 3. * 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;
                }
//                else std::cout << "-Rejecting triple with chiSquare: " << hitChiSquare << " and hiMinIndex: " << hiMinIndex << ", loMaxIndex: " << lowMaxIndex << std::endl;
            }
        }
    }
//    else std::cout << "-MakeTriplet hit cut, delta: " << hitTimeTicks - pair.getAvePeakTime() << ", min scale fctr: " <<m_hitWidthSclFctr << ", pair sig: " << pair.getSigmaPeakTime() << ", hitSigma: " << hitSigma << std::endl;

    // return success/fail
    return result;
}

void lar_cluster3d::SnippetHit3DBuilder::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 1785 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::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 SnippetHit3DBuilder_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::SnippetHit3DBuilder::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 1449 of file SnippetHit3DBuilder_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::SnippetHit3DBuilder::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 300 of file SnippetHit3DBuilder_tool.cc.

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

bool lar_cluster3d::SnippetHit3DBuilder::WireIDsIntersect ( const geo::WireID &  wireID0, const geo::WireID &  wireID1, geo::WireIDIntersection &  widIntersection  ) const

private

function to detemine if two wires "intersect" (in the 2D sense)

Definition at line 1224 of file SnippetHit3DBuilder_tool.cc.

{
    bool success(false);

    // Do quick check that things are in the same logical TPC
    if (wireID0.Cryostat != wireID1.Cryostat || wireID0.TPC != wireID1.TPC || wireID0.Plane == wireID1.Plane) return success;
        
    // Recover wire geometry information for each wire
    const geo::WireGeo& wireGeo0 = m_geometry->WireIDToWireGeo(wireID0);
    const geo::WireGeo& wireGeo1 = m_geometry->WireIDToWireGeo(wireID1);

    // Get wire position and direction for first wire
    double wirePosArr[3] = {0.,0.,0.};
    wireGeo0.GetCenter(wirePosArr);

    Eigen::Vector3f wirePos0(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
    Eigen::Vector3f wireDir0(wireGeo0.Direction().X(),wireGeo0.Direction().Y(),wireGeo0.Direction().Z());

    //*********************************
    // Kludge
//    if (wireID0.Plane > 0) wireDir0[2] = -wireDir0[2];

    // And now the second one
    wireGeo1.GetCenter(wirePosArr);

    Eigen::Vector3f wirePos1(wirePosArr[0],wirePosArr[1],wirePosArr[2]);
    Eigen::Vector3f wireDir1(wireGeo1.Direction().X(),wireGeo1.Direction().Y(),wireGeo1.Direction().Z());

    //**********************************
    // Kludge
//    if (wireID1.Plane > 0) wireDir1[2] = -wireDir1[2];

    // Get the distance of closest approach
    float arcLen0;
    float arcLen1;

    if (closestApproach(wirePos0, wireDir0, wirePos1, wireDir1, arcLen0, arcLen1))
    {
        // Now check that arc lengths are within range
        if (std::abs(arcLen0) < wireGeo0.HalfL() && std::abs(arcLen1) < wireGeo1.HalfL())
        {
            Eigen::Vector3f poca0 = wirePos0 + arcLen0 * wireDir0;

            widIntersection.y = poca0[1];
            widIntersection.z = poca0[2];

            success = true;
        }
    }

    return success;
}

Member Data Documentation

const lariov::ChannelStatusProvider* lar_cluster3d::SnippetHit3DBuilder::m_channelFilter

private

Definition at line 288 of file SnippetHit3DBuilder_tool.cc.

ChannelStatusByPlaneVec lar_cluster3d::SnippetHit3DBuilder::m_channelStatus

mutableprivate

Definition at line 282 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_chiSquare3DVec

mutableprivate

Definition at line 263 of file SnippetHit3DBuilder_tool.cc.

Hit2DList lar_cluster3d::SnippetHit3DBuilder::m_clusterHit2DMasterList

mutableprivate

Definition at line 277 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_deltaPeakTimeSig

private

Definition at line 243 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_deltaTimeVec

mutableprivate

Definition at line 262 of file SnippetHit3DBuilder_tool.cc.

bool lar_cluster3d::SnippetHit3DBuilder::m_enableMonitoring

private

Definition at line 253 of file SnippetHit3DBuilder_tool.cc.

const geo::Geometry* lar_cluster3d::SnippetHit3DBuilder::m_geometry

private

Definition at line 287 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_hitAsymmetryVec

mutableprivate

Definition at line 274 of file SnippetHit3DBuilder_tool.cc.

std::vector<art::InputTag> lar_cluster3d::SnippetHit3DBuilder::m_hitFinderTagVec

private

Data members to follow.

Definition at line 241 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_hitWidthSclFctr

private

Definition at line 242 of file SnippetHit3DBuilder_tool.cc.

std::vector<int> lar_cluster3d::SnippetHit3DBuilder::m_invalidTPCVec

private

Definition at line 248 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_LongHitStretchFctr

private

Definition at line 245 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxDeltaPeakVec

mutableprivate

Definition at line 267 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_maxHit3DChiSquare

private

Provide ability to select hits based on "chi square".

Definition at line 250 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxPullVec

mutableprivate

Definition at line 264 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_maxSideVecVec

mutableprivate

Definition at line 268 of file SnippetHit3DBuilder_tool.cc.

size_t lar_cluster3d::SnippetHit3DBuilder::m_numBadChannels

mutableprivate

Definition at line 283 of file SnippetHit3DBuilder_tool.cc.

bool lar_cluster3d::SnippetHit3DBuilder::m_outputHistograms

private

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

Definition at line 251 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_overlapFractionVec

mutableprivate

Definition at line 265 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_overlapRangeVec

mutableprivate

Definition at line 266 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_pairWireDistVec

mutableprivate

Definition at line 269 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_PHLowSelection

private

Definition at line 247 of file SnippetHit3DBuilder_tool.cc.

PlaneToSnippetHitMap lar_cluster3d::SnippetHit3DBuilder::m_planeToSnippetHitMap

mutableprivate

Definition at line 278 of file SnippetHit3DBuilder_tool.cc.

PlaneToWireToHitSetMap lar_cluster3d::SnippetHit3DBuilder::m_planeToWireToHitSetMap

mutableprivate

Definition at line 279 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_pulseHeightFrac

private

Definition at line 246 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_qualityMetricVec

mutableprivate

Definition at line 272 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_rangeNumSig

private

Definition at line 244 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_smallChargeDiffVec

mutableprivate

Definition at line 270 of file SnippetHit3DBuilder_tool.cc.

std::vector<int> lar_cluster3d::SnippetHit3DBuilder::m_smallIndexVec

mutableprivate

Definition at line 271 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_spacePointChargeVec

mutableprivate

Definition at line 273 of file SnippetHit3DBuilder_tool.cc.

std::vector<float> lar_cluster3d::SnippetHit3DBuilder::m_timeVector

mutableprivate

Definition at line 255 of file SnippetHit3DBuilder_tool.cc.

TTree* lar_cluster3d::SnippetHit3DBuilder::m_tupleTree

private

output analysis tree

Definition at line 260 of file SnippetHit3DBuilder_tool.cc.

bool lar_cluster3d::SnippetHit3DBuilder::m_weHaveAllBeenHereBefore = false

mutableprivate

Definition at line 285 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitch[3]

private

Definition at line 254 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitchScaleFactor

private

Scaling factor to determine max distance allowed between candidate pairs.

Definition at line 249 of file SnippetHit3DBuilder_tool.cc.

float lar_cluster3d::SnippetHit3DBuilder::m_zPosOffset

private

Definition at line 257 of file SnippetHit3DBuilder_tool.cc.

---

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

• SnippetHit3DBuilder_tool.cc