lar_cluster3d::SnippetHit3DBuilder Class Reference

SnippetHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::SnippetHit3DBuilder:

Public Types
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...

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...

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 &) 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, size_t minStatus) 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 () const
	Create the internal channel status vector (assume will eventually be event-by-event) More...
float	chargeIntegral (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

Detailed Description

SnippetHit3DBuilder class definiton.

Definition at line 80 of file SnippetHit3DBuilder_tool.cc.

Member Typedef Documentation

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

private

Definition at line 257 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 256 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 164 of file SnippetHit3DBuilder_tool.cc.

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

private

Definition at line 165 of file SnippetHit3DBuilder_tool.cc.

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

private

Definition at line 166 of file SnippetHit3DBuilder_tool.cc.

using lar_cluster3d::IHit3DBuilder::RecobHitToPtrMap = std::unordered_map<const recob::Hit*, art::Ptr<recob::Hit>>

inherited

Defines a structure mapping art representation to internal.

Definition at line 60 of file IHit3DBuilder.h.

Member Enumeration Documentation

enum lar_cluster3d::IHit3DBuilder::TimeValues

inherited

enumerate the possible values for time checking if monitoring timing

Enumerator
COLLECTARTHITS
BUILDTHREEDHITS
BUILDNEWHITS
NUMTIMEVALUES

Definition at line 73 of file IHit3DBuilder.h.

{ COLLECTARTHITS = 0, BUILDTHREEDHITS = 1, BUILDNEWHITS = 2, NUMTIMEVALUES };

Constructor & Destructor Documentation

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

explicit

Constructor.

Parameters

pset

Definition at line 317 of file SnippetHit3DBuilder_tool.cc.

References configure().

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

  {
    this->configure(pset);
  }

Member Function Documentation

void lar_cluster3d::SnippetHit3DBuilder::BuildChannelStatusVec ( ) const

private

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

Definition at line 402 of file SnippetHit3DBuilder_tool.cc.

References geo::GeometryCore::ChannelToWire(), m_channelFilter, m_channelStatus, m_geometry, m_numBadChannels, geo::GeometryCore::Nchannels(), geo::GeometryCore::Nplanes(), geo::GeometryCore::Nwires(), geo::PlaneID::Plane, and geo::WireID::Wire.

Referenced by BuildHit3D().

  {
    // 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
    constexpr geo::TPCID tpcid{0, 0};
    for (unsigned int idx = 0; idx < m_channelStatus.size(); idx++) {
      m_channelStatus[idx] = ChannelStatusVec(m_geometry->Nwires(geo::PlaneID{tpcid, 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)) continue;

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

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 500 of file SnippetHit3DBuilder_tool.cc.

References BuildChannelStatusVec(), BuildHitPairMap(), lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS, m_enableMonitoring, m_planeToSnippetHitMap, and m_timeVector.

Referenced by Hit3DBuilder().

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

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

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 548 of file SnippetHit3DBuilder_tool.cc.

References BuildHitPairMapByTPC(), m_geometry, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), and lar_cluster3d::SetPairStartTimeOrder().

Referenced by BuildHit3D().

  {
    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++) {
        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 609 of file SnippetHit3DBuilder_tool.cc.

References findGoodHitPairs(), and findGoodTriplets().

Referenced by BuildHitPairMap().

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

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

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

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

      snippetHitMapItrVec.front().first++;
    }

    return hitPairList.size();
  }

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

private

Perform charge integration between limits.

Definition at line 1351 of file SnippetHit3DBuilder_tool.cc.

Referenced by makeHitTriplet().

  {
    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 385 of file SnippetHit3DBuilder_tool.cc.

References m_chiSquare3DVec, m_deltaTimeVec, m_hitAsymmetryVec, m_maxDeltaPeakVec, m_maxPullVec, m_maxSideVecVec, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, m_smallChargeDiffVec, m_smallIndexVec, and m_spacePointChargeVec.

Referenced by configure(), and Hit3DBuilder().

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

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 1326 of file SnippetHit3DBuilder_tool.cc.

References d, den, and e.

Referenced by WireIDsIntersect().

  {
    // 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 1587 of file SnippetHit3DBuilder_tool.cc.

References geo::GeometryCore::ChannelToWire(), lar_cluster3d::IHit3DBuilder::COLLECTARTHITS, geo::CryostatID::Cryostat, art::ProductRetriever::getByLabel(), art::Handle< T >::isValid(), m_clusterHit2DMasterList, m_enableMonitoring, m_geometry, m_hitFinderTagVec, m_invalidTPCVec, m_planeToSnippetHitMap, m_planeToWireToHitSetMap, m_timeVector, m_weHaveAllBeenHereBefore, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), geo::PlaneID::Plane, geo::TPCID::TPC, and detinfo::trigger_offset().

Referenced by Hit3DBuilder().

  {
    // 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 335 of file SnippetHit3DBuilder_tool.cc.

References clear(), art::ServiceHandle< T, SCOPE >::get(), fhicl::ParameterSet::get(), m_chiSquare3DVec, m_deltaPeakTimeSig, m_deltaTimeVec, m_enableMonitoring, m_geometry, m_hitAsymmetryVec, m_hitFinderTagVec, m_hitWidthSclFctr, m_invalidTPCVec, m_LongHitStretchFctr, m_maxDeltaPeakVec, m_maxHit3DChiSquare, m_maxPullVec, m_maxSideVecVec, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_PHLowSelection, m_pulseHeightFrac, m_qualityMetricVec, m_rangeNumSig, m_smallChargeDiffVec, m_smallIndexVec, m_spacePointChargeVec, m_tupleTree, m_wirePitch, m_wirePitchScaleFactor, and m_zPosOffset.

Referenced by SnippetHit3DBuilder().

  {
    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
    constexpr geo::TPCID tpcid{0, 0};
    m_wirePitch[0] = m_geometry->WirePitch(geo::PlaneID{tpcid, 0});
    m_wirePitch[1] = m_geometry->WirePitch(geo::PlaneID{tpcid, 1});
    m_wirePitch[2] = m_geometry->WirePitch(geo::PlaneID{tpcid, 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);
    }
  }

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 1738 of file SnippetHit3DBuilder_tool.cc.

References lar_cluster3d::IHit3DBuilder::BUILDNEWHITS, reco::ClusterHit2D::getHit(), m_clusterHit2DMasterList, m_enableMonitoring, and m_timeVector.

Referenced by Hit3DBuilder().

  {
    // 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 1527 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), geo::WireGeo::Direction(), geo::vect::dot(), geo::WireGeo::GetCenter(), geo::WireGeo::HalfL(), m_geometry, and geo::GeometryCore::WireIDToWireGeo().

Referenced by makeHitTriplet().

  {
    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
      auto const wirePosArr = wireGeo.GetCenter();

      Eigen::Vector3f wirePos(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
      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 1444 of file SnippetHit3DBuilder_tool.cc.

References reco::ClusterHit3D::getAvePeakTime().

  {
    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 696 of file SnippetHit3DBuilder_tool.cc.

References art::left(), m_hitWidthSclFctr, m_PHLowSelection, m_pulseHeightFrac, makeHitPair(), art::right(), and geo::WireID::WireID().

Referenced by BuildHitPairMapByTPC().

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

private

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

Definition at line 766 of file SnippetHit3DBuilder_tool.cc.

References art::left(), m_deltaPeakTimeSig, m_numBadChannels, makeDeadChannelPair(), makeHitTriplet(), art::right(), and reco::ClusterHit3D::setID().

Referenced by BuildHitPairMapByTPC().

  {
    // 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))
            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 1473 of file SnippetHit3DBuilder_tool.cc.

References reco::ClusterHit3D::getAvePeakTime().

  {
    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 108 of file SnippetHit3DBuilder_tool.cc.

References tca::evt.

    {
      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 446 of file SnippetHit3DBuilder_tool.cc.

References BuildHit3D(), clear(), CollectArtHits(), CreateNewRecobHitCollection(), m_clusterHit2DMasterList, m_outputHistograms, m_planeToSnippetHitMap, m_planeToWireToHitSetMap, m_timeVector, m_tupleTree, makeRawDigitAssns(), makeWireAssns(), lar_cluster3d::IHit3DBuilder::NUMTIMEVALUES, and art::Event::put().

  {
    // 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)
    std::unique_ptr<art::Assns<recob::Wire, recob::Hit>> wireAssns(
      new art::Assns<recob::Wire, recob::Hit>);

    makeWireAssns(evt, *wireAssns, clusterHitToArtPtrMap);

    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();
    }
  }

bool lar_cluster3d::SnippetHit3DBuilder::makeDeadChannelPair ( reco::ClusterHit3D &  pairOut, const reco::ClusterHit3D &  pair, size_t  maxStatus, size_t  minStatus  ) const

private

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

Definition at line 1367 of file SnippetHit3DBuilder_tool.cc.

References geo::CryostatID::Cryostat, reco::ClusterHit3D::getHits(), reco::ClusterHit3D::getPosition(), reco::ClusterHit2D::getStatusBits(), m_channelStatus, m_geometry, m_zPosOffset, NearestWireID(), geo::PlaneID::Plane, reco::ClusterHit3D::setPosition(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit3D::setWireID(), reco::ClusterHit2D::SHAREDINTRIPLET, geo::TPCID::TPC, reco::ClusterHit2D::USEDINTRIPLET, geo::WireID::Wire, reco::ClusterHit2D::WireID(), geo::GeometryCore::WireIDsIntersect(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by findGoodTriplets().

  {
    // 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 886 of file SnippetHit3DBuilder_tool.cc.

References geo::CryostatID::Cryostat, recob::Hit::DegreesOfFreedom(), reco::ClusterHit2D::getHit(), reco::ClusterHit2D::getStatusBits(), reco::ClusterHit2D::getTimeTicks(), reco::ClusterHit2D::getXPosition(), reco::ClusterHit3D::initialize(), recob::Hit::Integral(), m_deltaPeakTimeSig, m_LongHitStretchFctr, m_zPosOffset, geo::PlaneID::Plane, recob::Hit::RMS(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit2D::SHAREDINPAIR, geo::TPCID::TPC, reco::ClusterHit2D::USEDINPAIR, reco::ClusterHit2D::WireID(), WireIDsIntersect(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by findGoodHitPairs(), and makeHitTriplet().

  {
    // 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 1014 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), chargeIntegral(), geo::CryostatID::Cryostat, recob::Hit::DegreesOfFreedom(), DistanceFromPointToHitWire(), reco::ClusterHit3D::getAvePeakTime(), reco::ClusterHit2D::getHit(), reco::ClusterHit3D::getHits(), reco::ClusterHit3D::getPosition(), reco::ClusterHit3D::getSigmaPeakTime(), reco::ClusterHit2D::getStatusBits(), reco::ClusterHit2D::getTimeTicks(), reco::ClusterHit2D::getXPosition(), reco::ClusterHit3D::initialize(), m_chiSquare3DVec, m_deltaTimeVec, m_hitAsymmetryVec, m_hitWidthSclFctr, m_LongHitStretchFctr, m_maxDeltaPeakVec, m_maxPullVec, m_maxSideVecVec, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, m_rangeNumSig, m_smallChargeDiffVec, m_smallIndexVec, m_spacePointChargeVec, m_wirePitch, m_wirePitchScaleFactor, makeHitPair(), recob::Hit::PeakAmplitude(), geo::PlaneID::Plane, recob::Hit::RMS(), reco::ClusterHit2D::setStatusBit(), reco::ClusterHit2D::SHAREDINTRIPLET, geo::TPCID::TPC, reco::ClusterHit2D::USEDINTRIPLET, weight, geo::WireID::Wire, reco::ClusterHit2D::WireID(), and geo::WireID::WireID().

Referenced by findGoodTriplets().

  {
    // 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(),
                                                 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 1847 of file SnippetHit3DBuilder_tool.cc.

References art::Assns< L, R, D >::addSingle(), raw::RawDigit::Channel(), DEFINE_ART_CLASS_TOOL, art::Assns< L, R, D >::end(), art::ProductRetriever::getValidHandle(), and m_hitFinderTagVec.

Referenced by Hit3DBuilder().

  {
    // 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 1802 of file SnippetHit3DBuilder_tool.cc.

References art::Assns< L, R, D >::addSingle(), recob::Wire::Channel(), art::Assns< L, R, D >::end(), art::ProductRetriever::getValidHandle(), and m_hitFinderTagVec.

Referenced by Hit3DBuilder().

  {
    // 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 1499 of file SnippetHit3DBuilder_tool.cc.

References geo::PlaneID::asPlaneID(), geo::GeometryCore::DetLength(), m_geometry, geo::GeometryCore::Nwires(), geo::GeometryCore::Plane(), geo::vect::toPoint(), geo::WireID::Wire, and geo::PlaneGeo::WireCoordinate().

Referenced by makeDeadChannelPair().

  {
    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(geo::vect::toPoint(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.asPlaneID()) - 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 326 of file SnippetHit3DBuilder_tool.cc.

References art::ProducesCollector::produces().

  {
    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 1270 of file SnippetHit3DBuilder_tool.cc.

References util::abs(), closestApproach(), geo::CryostatID::Cryostat, geo::WireGeo::Direction(), geo::WireGeo::GetCenter(), geo::WireGeo::HalfL(), m_geometry, geo::PlaneID::Plane, geo::TPCID::TPC, geo::GeometryCore::WireIDToWireGeo(), geo::WireIDIntersection::y, and geo::WireIDIntersection::z.

Referenced by makeHitPair().

  {
    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
    auto wirePosArr = wireGeo0.GetCenter();

    Eigen::Vector3f wirePos0(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
    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
    wirePosArr = wireGeo1.GetCenter();

    Eigen::Vector3f wirePos1(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
    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 314 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec().

ChannelStatusByPlaneVec lar_cluster3d::SnippetHit3DBuilder::m_channelStatus

mutableprivate

Definition at line 308 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and makeDeadChannelPair().

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

mutableprivate

Definition at line 290 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

Hit2DList lar_cluster3d::SnippetHit3DBuilder::m_clusterHit2DMasterList

mutableprivate

Definition at line 304 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), CreateNewRecobHitCollection(), and Hit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_deltaPeakTimeSig

private

Definition at line 269 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), findGoodTriplets(), and makeHitPair().

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

mutableprivate

Definition at line 289 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

bool lar_cluster3d::SnippetHit3DBuilder::m_enableMonitoring

private

Definition at line 280 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildHit3D(), CollectArtHits(), configure(), and CreateNewRecobHitCollection().

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

private

Definition at line 313 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), BuildHitPairMap(), CollectArtHits(), configure(), DistanceFromPointToHitWire(), makeDeadChannelPair(), NearestWireID(), and WireIDsIntersect().

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

mutableprivate

Definition at line 301 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

private

Data members to follow.

Definition at line 267 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), configure(), makeRawDigitAssns(), and makeWireAssns().

float lar_cluster3d::SnippetHit3DBuilder::m_hitWidthSclFctr

private

Definition at line 268 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), findGoodHitPairs(), and makeHitTriplet().

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

private

Definition at line 274 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and configure().

float lar_cluster3d::SnippetHit3DBuilder::m_LongHitStretchFctr

private

Definition at line 271 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), makeHitPair(), and makeHitTriplet().

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

mutableprivate

Definition at line 294 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

float lar_cluster3d::SnippetHit3DBuilder::m_maxHit3DChiSquare

private

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

Definition at line 277 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure().

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

mutableprivate

Definition at line 291 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 295 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

size_t lar_cluster3d::SnippetHit3DBuilder::m_numBadChannels

mutableprivate

Definition at line 309 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and findGoodTriplets().

bool lar_cluster3d::SnippetHit3DBuilder::m_outputHistograms

private

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

Definition at line 278 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), Hit3DBuilder(), and makeHitTriplet().

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

mutableprivate

Definition at line 292 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 293 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 296 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

float lar_cluster3d::SnippetHit3DBuilder::m_PHLowSelection

private

Definition at line 273 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and findGoodHitPairs().

PlaneToSnippetHitMap lar_cluster3d::SnippetHit3DBuilder::m_planeToSnippetHitMap

mutableprivate

Definition at line 305 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildHit3D(), CollectArtHits(), and Hit3DBuilder().

PlaneToWireToHitSetMap lar_cluster3d::SnippetHit3DBuilder::m_planeToWireToHitSetMap

mutableprivate

Definition at line 306 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and Hit3DBuilder().

float lar_cluster3d::SnippetHit3DBuilder::m_pulseHeightFrac

private

Definition at line 272 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and findGoodHitPairs().

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

mutableprivate

Definition at line 299 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

float lar_cluster3d::SnippetHit3DBuilder::m_rangeNumSig

private

Definition at line 270 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 297 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 298 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 300 of file SnippetHit3DBuilder_tool.cc.

Referenced by clear(), configure(), and makeHitTriplet().

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

mutableprivate

Definition at line 282 of file SnippetHit3DBuilder_tool.cc.

Referenced by BuildHit3D(), CollectArtHits(), CreateNewRecobHitCollection(), and Hit3DBuilder().

TTree* lar_cluster3d::SnippetHit3DBuilder::m_tupleTree

private

output analysis tree

Definition at line 287 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and Hit3DBuilder().

bool lar_cluster3d::SnippetHit3DBuilder::m_weHaveAllBeenHereBefore = false

mutableprivate

Definition at line 311 of file SnippetHit3DBuilder_tool.cc.

Referenced by CollectArtHits().

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitch[3]

private

Definition at line 281 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and makeHitTriplet().

float lar_cluster3d::SnippetHit3DBuilder::m_wirePitchScaleFactor

private

Scaling factor to determine max distance allowed between candidate pairs.

Definition at line 276 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), and makeHitTriplet().

float lar_cluster3d::SnippetHit3DBuilder::m_zPosOffset

private

Definition at line 284 of file SnippetHit3DBuilder_tool.cc.

Referenced by configure(), makeDeadChannelPair(), and makeHitPair().

---

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

• larreco/v09_25_00/source/larreco/RecoAlg/Cluster3DAlgs/SnippetHit3DBuilder_tool.cc