lar_cluster3d::StandardHit3DBuilder Class Reference

StandardHit3DBuilder class definiton. More...

Inheritance diagram for lar_cluster3d::StandardHit3DBuilder:

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
	StandardHit3DBuilder (fhicl::ParameterSet const &pset)
	Constructor. More...
void	produces (art::ProducesCollector &) override
	Each algorithm may have different objects it wants "produced" so use this to let the top level producer module "know" what it is outputting. More...
void	Hit3DBuilder (art::Event &, reco::HitPairList &, RecobHitToPtrMap &) override
	Given a set of recob hits, run DBscan to form 3D clusters. More...
float	getTimeToExecute (IHit3DBuilder::TimeValues index) const override
	If monitoring, recover the time to execute a particular function. More...

Private Types
using	PlaneHitVectorItrPairVec = std::vector< std::pair< HitVector::iterator, HitVector::iterator >>
	Given the ClusterHit2D objects, build the HitPairMap. More...
using	HitMatchPair = std::pair< const reco::ClusterHit2D *, reco::ClusterHit3D >
	This builds a list of candidate hit pairs from lists of hits on two planes. More...
using	HitMatchPairVec = std::vector< HitMatchPair >
using	HitMatchPairVecMap = std::map< geo::WireID, HitMatchPairVec >
using	ChannelStatusVec = std::vector< size_t >
	define data structure for keeping track of channel status More...
using	ChannelStatusByPlaneVec = std::vector< ChannelStatusVec >

Private Member Functions
void	CollectArtHits (const art::Event &evt) const
	Extract the ART hits and the ART hit-particle relationships. More...
void	BuildHit3D (reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...
void	CreateNewRecobHitCollection (art::Event &, reco::HitPairList &, std::vector< recob::Hit > &, RecobHitToPtrMap &)
	Create a new 2D hit collection from hits associated to 3D space points. More...
void	makeWireAssns (const art::Event &, art::Assns< recob::Wire, recob::Hit > &, RecobHitToPtrMap &) const
	Create recob::Wire to recob::Hit associations. More...
void	makeRawDigitAssns (const art::Event &, art::Assns< raw::RawDigit, recob::Hit > &, RecobHitToPtrMap &) const
	Create raw::RawDigit to recob::Hit associations. More...
size_t	BuildHitPairMap (PlaneToHitVectorMap &planeToHitVectorMap, reco::HitPairList &hitPairList) const
	Given the ClusterHit2D objects, build the HitPairMap. More...
size_t	BuildHitPairMapByTPC (PlaneHitVectorItrPairVec &planeHitVectorItrPairVec, reco::HitPairList &hitPairList) const
int	findGoodHitPairs (const reco::ClusterHit2D *, HitVector::iterator &, HitVector::iterator &, HitMatchPairVecMap &) const
void	findGoodTriplets (HitMatchPairVecMap &, HitMatchPairVecMap &, reco::HitPairList &) 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...
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_numSigmaPeakTime
float	m_hitWidthSclFctr
float	m_deltaPeakTimeSig
std::vector< int >	m_invalidTPCVec
float	m_wirePitchScaleFactor
	Scaling factor to determine max distance allowed between candidate pairs. More...
float	m_maxHit3DChiSquare
	Provide ability to select hits based on "chi square". More...
bool	m_outputHistograms
	Take the time to create and fill some histograms for diagnostics. More...
bool	m_enableMonitoring
float	m_wirePitch [3]
std::vector< float >	m_timeVector
float	m_zPosOffset
TTree *	m_tupleTree
	output analysis tree More...
std::vector< float >	m_deltaTimeVec
std::vector< float >	m_chiSquare3DVec
std::vector< float >	m_maxPullVec
std::vector< float >	m_overlapFractionVec
std::vector< float >	m_overlapRangeVec
std::vector< float >	m_maxDeltaPeakVec
std::vector< float >	m_maxSideVecVec
std::vector< float >	m_pairWireDistVec
std::vector< float >	m_smallChargeDiffVec
std::vector< int >	m_smallIndexVec
std::vector< float >	m_qualityMetricVec
std::vector< float >	m_spacePointChargeVec
std::vector< float >	m_hitAsymmetryVec
Hit2DList	m_clusterHit2DMasterList
PlaneToHitVectorMap	m_planeToHitVectorMap
PlaneToWireToHitSetMap	m_planeToWireToHitSetMap
ChannelStatusByPlaneVec	m_channelStatus
size_t	m_numBadChannels
bool	m_weHaveAllBeenHereBefore = false
const geo::Geometry *	m_geometry
const geo::WireReadoutGeom *	m_wireReadoutGeom
const lariov::ChannelStatusProvider *	m_channelFilter

Detailed Description

StandardHit3DBuilder class definiton.

Definition at line 76 of file StandardHit3DBuilder_tool.cc.

Member Typedef Documentation

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

private

Definition at line 237 of file StandardHit3DBuilder_tool.cc.

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

private

define data structure for keeping track of channel status

Definition at line 236 of file StandardHit3DBuilder_tool.cc.

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

private

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

Definition at line 160 of file StandardHit3DBuilder_tool.cc.

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

private

Definition at line 161 of file StandardHit3DBuilder_tool.cc.

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

private

Definition at line 162 of file StandardHit3DBuilder_tool.cc.

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

private

Given the ClusterHit2D objects, build the HitPairMap.

Definition at line 152 of file StandardHit3DBuilder_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 43 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 56 of file IHit3DBuilder.h.

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

Constructor & Destructor Documentation

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

explicit

Constructor.

Parameters

pset

Definition at line 295 of file StandardHit3DBuilder_tool.cc.

References clear(), Get, art::ServiceHandle< T, SCOPE >::get(), m_channelFilter, m_chiSquare3DVec, m_deltaPeakTimeSig, m_deltaTimeVec, m_enableMonitoring, m_geometry, m_hitAsymmetryVec, m_hitFinderTagVec, m_hitWidthSclFctr, m_invalidTPCVec, m_maxDeltaPeakVec, m_maxHit3DChiSquare, m_maxPullVec, m_maxSideVecVec, m_numSigmaPeakTime, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, m_smallChargeDiffVec, m_smallIndexVec, m_spacePointChargeVec, m_tupleTree, m_wirePitch, m_wirePitchScaleFactor, m_wireReadoutGeom, m_zPosOffset, and geo::WireReadoutGeom::Plane().

    : m_geometry(art::ServiceHandle<geo::Geometry const>{}.get())
    , m_channelFilter(&art::ServiceHandle<lariov::ChannelStatusService const>()->GetProvider())
  {
    m_hitFinderTagVec = pset.get<std::vector<art::InputTag>>("HitFinderTagVec",
                                                             std::vector<art::InputTag>{"gaushit"});
    m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
    m_numSigmaPeakTime = pset.get<float>("NumSigmaPeakTime", 3.);
    m_hitWidthSclFctr = pset.get<float>("HitWidthScaleFactor", 6.);
    m_deltaPeakTimeSig = pset.get<float>("DeltaPeakTimeSig", 1.7);
    m_zPosOffset = pset.get<float>("ZPosOffset", 0.0);
    m_invalidTPCVec = pset.get<std::vector<int>>("InvalidTPCVec", std::vector<int>{});
    m_wirePitchScaleFactor = pset.get<float>("WirePitchScaleFactor", 1.9);
    m_maxHit3DChiSquare = pset.get<float>("MaxHitChiSquare", 6.0);
    m_outputHistograms = pset.get<bool>("OutputHistograms", false);

    m_geometry = art::ServiceHandle<geo::Geometry const>{}.get();
    m_wireReadoutGeom = &art::ServiceHandle<geo::WireReadout 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_wireReadoutGeom->Plane({tpcid, 0}).WirePitch();
    m_wirePitch[1] = m_wireReadoutGeom->Plane({tpcid, 1}).WirePitch();
    m_wirePitch[2] = m_wireReadoutGeom->Plane({tpcid, 2}).WirePitch();

    if (m_outputHistograms) {
      // Access ART's TFileService, which will handle creating and writing
      // histograms and n-tuples for us.
      art::ServiceHandle<art::TFileService> tfs;

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

      clear();

      m_tupleTree->Branch("DeltaTime2D", "std::vector<float>", &m_deltaTimeVec);
      m_tupleTree->Branch("ChiSquare3D", "std::vector<float>", &m_chiSquare3DVec);
      m_tupleTree->Branch("MaxPullValue", "std::vector<float>", &m_maxPullVec);
      m_tupleTree->Branch("OverlapFraction", "std::vector<float>", &m_overlapFractionVec);
      m_tupleTree->Branch("OverlapRange", "std::vector<float>", &m_overlapRangeVec);
      m_tupleTree->Branch("MaxDeltaPeak", "std::vector<float>", &m_maxDeltaPeakVec);
      m_tupleTree->Branch("MaxSideVec", "std::vector<float>", &m_maxSideVecVec);
      m_tupleTree->Branch("PairWireDistVec", "std::vector<float>", &m_pairWireDistVec);
      m_tupleTree->Branch("SmallChargeDiff", "std::vector<float>", &m_smallChargeDiffVec);
      m_tupleTree->Branch("SmallChargeIdx", "std::vector<int>", &m_smallIndexVec);
      m_tupleTree->Branch("QualityMetric", "std::vector<float>", &m_qualityMetricVec);
      m_tupleTree->Branch("SPCharge", "std::vector<float>", &m_spacePointChargeVec);
      m_tupleTree->Branch("HitAsymmetry", "std::vector<float>", &m_hitAsymmetryVec);
    }
  }

Member Function Documentation

void lar_cluster3d::StandardHit3DBuilder::BuildChannelStatusVec ( ) const

private

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

Definition at line 373 of file StandardHit3DBuilder_tool.cc.

References geo::WireReadoutGeom::ChannelToWire(), art::left(), m_channelFilter, m_channelStatus, m_numBadChannels, m_wireReadoutGeom, geo::WireReadoutGeom::Nchannels(), geo::WireReadoutGeom::Nplanes(), geo::WireReadoutGeom::Nwires(), geo::PlaneID::Plane, art::right(), lar_cluster3d::SetPeakHitPairIteratorOrder(), and geo::WireID::Wire.

Referenced by BuildHit3D().

  {
    // This is called each event, clear out the previous version and start over
    if (!m_channelStatus.empty()) m_channelStatus.clear();

    m_numBadChannels = 0;
    m_channelStatus.resize(m_wireReadoutGeom->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_wireReadoutGeom->Nwires(geo::PlaneID{tpcid, idx}), 5);
    }

    // Loop through the channels and mark those that are "bad"
    for (size_t channel = 0; channel < m_wireReadoutGeom->Nchannels(); channel++) {
      if (m_channelFilter->IsGood(channel)) continue;

      std::vector<geo::WireID> wireIDVec = m_wireReadoutGeom->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::StandardHit3DBuilder::BuildHit3D ( reco::HitPairList &  hitPairList ) const

private

Given the ClusterHit2D objects, build the HitPairMap.

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

Definition at line 472 of file StandardHit3DBuilder_tool.cc.

References BuildChannelStatusVec(), BuildHitPairMap(), lar_cluster3d::IHit3DBuilder::BUILDTHREEDHITS, reco::ClusterHit2D::getHit(), reco::ClusterHit2D::getTimeTicks(), m_enableMonitoring, m_numRMS, m_planeToHitVectorMap, m_timeVector, and recob::Hit::RMS().

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_planeToHitVectorMap, hitPairList);

    if (m_enableMonitoring) {
      theClockMakeHits.stop();

      m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time();
    }

    mf::LogDebug("Cluster3D") << ">>>>> 3D hit building done, found " << numHitPairs << " 3D Hits"
                              << std::endl;
  }

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

private

Given the ClusterHit2D objects, build the HitPairMap.

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

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

Definition at line 549 of file StandardHit3DBuilder_tool.cc.

References BuildHitPairMapByTPC(), m_geometry, m_numSigmaPeakTime, 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++) {
        PlaneToHitVectorMap::iterator mapItr0 =
          planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0));
        PlaneToHitVectorMap::iterator mapItr1 =
          planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1));
        PlaneToHitVectorMap::iterator mapItr2 =
          planeToHitVectorMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2));

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

        if (nPlanesWithHits < 2) continue;

        HitVector& hitVector0 = mapItr0->second;
        HitVector& hitVector1 = mapItr1->second;
        HitVector& hitVector2 = mapItr2->second;

        // We are going to resort the hits into "start time" order...
        std::sort(hitVector0.begin(),
                  hitVector0.end(),
                  SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);
        std::sort(hitVector1.begin(),
                  hitVector1.end(),
                  SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);
        std::sort(hitVector2.begin(),
                  hitVector2.end(),
                  SetHitEarliestTimeOrder(m_numSigmaPeakTime)); //SetHitStartTimeOrder);

        PlaneHitVectorItrPairVec hitItrVec = {
          HitVectorItrPair(hitVector0.begin(), hitVector0.end()),
          HitVectorItrPair(hitVector1.begin(), hitVector1.end()),
          HitVectorItrPair(hitVector2.begin(), hitVector2.end())};

        totalNumHits += BuildHitPairMapByTPC(hitItrVec, hitPairList);
      }
    }

    // Return the hit pair list but sorted by z and y positions (faster traversal in next steps)
    hitPairList.sort(SetPairStartTimeOrder);

    // Where are we?
    mf::LogDebug("Cluster3D") << "Total number hits: " << totalNumHits << std::endl;
    mf::LogDebug("Cluster3D") << "Created a total of " << hitPairList.size()
                              << " hit pairs, counted: " << hitPairCntr << std::endl;
    mf::LogDebug("Cluster3D") << "-- Triplets: " << nTriplets
                              << ", dead channel pairs: " << nDeadChanHits << std::endl;

    return hitPairList.size();
  }

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

private

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

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

Definition at line 621 of file StandardHit3DBuilder_tool.cc.

References findGoodHitPairs(), findGoodTriplets(), reco::ClusterHit2D::getTimeTicks(), and m_numSigmaPeakTime.

Referenced by BuildHitPairMap().

  {
    // Define functions to set start/end iterators in the loop below
    auto SetStartIterator =
      [](HitVector::iterator startItr, HitVector::iterator endItr, float rms, float startTime) {
        while (startItr != endItr) {
          float numRMS(rms);
          if ((*startItr)->getTimeTicks() + numRMS * (*startItr)->getHit()->RMS() < startTime)
            startItr++;
          else
            break;
        }
        return startItr;
      };

    auto SetEndIterator =
      [](HitVector::iterator firstItr, HitVector::iterator endItr, float rms, float endTime) {
        while (firstItr != endItr) {
          float numRMS(rms);
          if ((*firstItr)->getTimeTicks() - numRMS * (*firstItr)->getHit()->RMS() < endTime)
            firstItr++;
          else
            break;
        }
        return firstItr;
      };

    //*********************************************************************************
    // Basically, we try to loop until done...
    while (1) {
      // Sort so that the earliest hit time will be the first element, etc.
      std::sort(hitItrVec.begin(), hitItrVec.end(), SetStartTimeOrder(m_numSigmaPeakTime));

      // This loop iteration's golden hit
      const reco::ClusterHit2D* goldenHit = *hitItrVec[0].first;

      // The range of history... (for this hit)
      float goldenTimeStart = goldenHit->getTimeTicks() -
                              m_numSigmaPeakTime * goldenHit->getHit()->RMS() -
                              std::numeric_limits<float>::epsilon();
      float goldenTimeEnd = goldenHit->getTimeTicks() +
                            m_numSigmaPeakTime * goldenHit->getHit()->RMS() +
                            std::numeric_limits<float>::epsilon();

      // Set iterators to insure we'll be in the overlap ranges
      HitVector::iterator hitItr1Start = SetStartIterator(
        hitItrVec[1].first, hitItrVec[1].second, m_numSigmaPeakTime, goldenTimeStart);
      HitVector::iterator hitItr1End =
        SetEndIterator(hitItr1Start, hitItrVec[1].second, m_numSigmaPeakTime, goldenTimeEnd);
      HitVector::iterator hitItr2Start = SetStartIterator(
        hitItrVec[2].first, hitItrVec[2].second, m_numSigmaPeakTime, goldenTimeStart);
      HitVector::iterator hitItr2End =
        SetEndIterator(hitItr2Start, hitItrVec[2].second, m_numSigmaPeakTime, goldenTimeEnd);

      // Since we'll use these many times in the internal loops, pre make the pairs for the second set of hits
      HitMatchPairVecMap pair12Map;
      HitMatchPairVecMap pair13Map;

      size_t n12Pairs = findGoodHitPairs(goldenHit, hitItr1Start, hitItr1End, pair12Map);
      size_t n13Pairs = findGoodHitPairs(goldenHit, hitItr2Start, hitItr2End, pair13Map);

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

      hitItrVec[0].first++;

      int nPlanesWithHits(0);

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

      if (nPlanesWithHits < 2) break;
    }

    return hitPairList.size();
  }

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

private

Perform charge integration between limits.

Definition at line 1236 of file StandardHit3DBuilder_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::StandardHit3DBuilder::clear ( )

private

clear the tuple vectors before processing next event

Definition at line 356 of file StandardHit3DBuilder_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 Hit3DBuilder(), and StandardHit3DBuilder().

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

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

private

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

Parameters

evt	- the ART event

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

Definition at line 1462 of file StandardHit3DBuilder_tool.cc.

References geo::WireReadoutGeom::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_planeToHitVectorMap, m_planeToWireToHitSetMap, m_timeVector, m_weHaveAllBeenHereBefore, m_wireReadoutGeom, geo::GeometryCore::Ncryostats(), geo::GeometryCore::NTPC(), geo::PlaneID::Plane, lar_cluster3d::SetHitTimeOrder(), 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;

    auto const clock_data =
      art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
    auto const det_prop =
      art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);

    // Loop through the list of input sources
    for (const auto& inputTag : m_hitFinderTagVec) {
      art::Handle<std::vector<recob::Hit>> recobHitHandle;
      evt.getByLabel(inputTag, recobHitHandle);

      if (!recobHitHandle.isValid() || recobHitHandle->size() == 0) continue;

      recobHitVec.reserve(recobHitVec.size() + recobHitHandle->size());

      for (const auto& hit : *recobHitHandle)
        recobHitVec.push_back(&hit);
    }

    // If the vector is empty there is nothing to do
    if (recobHitVec.empty()) return;

    cet::cpu_timer theClockMakeHits;

    if (m_enableMonitoring) theClockMakeHits.start();

    // We'll want to correct the hit times for the plane offsets
    // (note this is already taken care of when converting to position)
    std::map<geo::PlaneID, double> planeOffsetMap;

    // Try to output a formatted string
    std::string debugMessage("");

    // Initialize the plane to hit vector map
    for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) {
      for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) {
        m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = HitVector();
        m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] = HitVector();
        m_planeToHitVectorMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] = HitVector();

        // What we want here are the relative offsets between the planes
        // Note that plane 0 is assumed the "first" plane and is the reference
        planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)] = 0.;
        planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)] =
          det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1)) -
          det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));
        planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] =
          det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2)) -
          det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0));

        // Should we provide output?
        if (!m_weHaveAllBeenHereBefore) {
          std::ostringstream outputString;

          outputString << "***> plane 0 offset: "
                       << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 0)]
                       << ", plane 1: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 1)]
                       << ", plane    2: " << planeOffsetMap[geo::PlaneID(cryoIdx, tpcIdx, 2)]
                       << "\n";
          debugMessage += outputString.str();
          outputString << "     Det prop plane 0: "
                       << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0))
                       << ", plane 1: "
                       << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 1))
                       << ", plane 2: "
                       << det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 2))
                       << ", Trig: " << trigger_offset(clock_data) << "\n";
          debugMessage += outputString.str();
        }
      }
    }

    if (!m_weHaveAllBeenHereBefore) {
      mf::LogDebug("Cluster3D") << debugMessage << std::endl;

      m_weHaveAllBeenHereBefore = true;
    }

    // Cycle through the recob hits to build ClusterHit2D objects and insert
    // them into the map
    for (const auto& recobHit : recobHitVec) {
      // Reject hits with negative charge, these are misreconstructed
      if (recobHit->Integral() < 0.) continue;

      // For some detectors we can have multiple wire ID's associated to a given channel.
      // So we recover the list of these wire IDs
      const std::vector<geo::WireID>& wireIDs =
        m_wireReadoutGeom->ChannelToWire(recobHit->Channel());

      // And then loop over all possible to build out our maps
      for (const auto& wireID : wireIDs) {
        // Check if this is an invalid TPC
        // (for example, in protoDUNE there are logical TPC's which see no signal)
        if (std::find(m_invalidTPCVec.begin(), m_invalidTPCVec.end(), wireID.TPC) !=
            m_invalidTPCVec.end())
          continue;

        // Note that a plane ID will define cryostat, TPC and plane
        const geo::PlaneID& planeID = wireID.planeID();

        double hitPeakTime(recobHit->PeakTime() - planeOffsetMap[planeID]);
        double xPosition(det_prop.ConvertTicksToX(
          recobHit->PeakTime(), planeID.Plane, planeID.TPC, planeID.Cryostat));

        m_clusterHit2DMasterList.emplace_back(0, 0., 0., xPosition, hitPeakTime, wireID, recobHit);

        m_planeToHitVectorMap[planeID].push_back(&m_clusterHit2DMasterList.back());
        m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back());
      }
    }

    // Make a loop through to sort the recover hits in time order
    for (auto& hitVectorMap : m_planeToHitVectorMap)
      std::sort(hitVectorMap.second.begin(), hitVectorMap.second.end(), SetHitTimeOrder);

    if (m_enableMonitoring) {
      theClockMakeHits.stop();

      m_timeVector[COLLECTARTHITS] = theClockMakeHits.accumulated_real_time();
    }

    mf::LogDebug("Cluster3D") << ">>>>> Number of ART hits: " << m_clusterHit2DMasterList.size()
                              << std::endl;
  }

void lar_cluster3d::StandardHit3DBuilder::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 1598 of file StandardHit3DBuilder_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::StandardHit3DBuilder::DistanceFromPointToHitWire ( const Eigen::Vector3f &  position, const geo::WireID &  wireID  ) const

private

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

Definition at line 1408 of file StandardHit3DBuilder_tool.cc.

References geo::vect::dot(), reco::ClusterHit2D::getHit(), m_wireReadoutGeom, lar_cluster3d::Hit2DSetCompare::operator()(), recob::Hit::PeakTime(), lar_cluster3d::SetHitTimeOrder(), and geo::WireReadoutGeom::WireEndPoints().

Referenced by makeHitTriplet().

  {
    float distance;

    // Embed the call to the geometry's services nearest wire id method in a try-catch block
    try {
      // Get the wire endpoints
      Eigen::Vector3d wireStart;
      Eigen::Vector3d wireEnd;

      m_wireReadoutGeom->WireEndPoints(wireIDIn, &wireStart[0], &wireEnd[0]);

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

      // Want the wire direction
      Eigen::Vector3d wireDir = wireEnd - wireStart;

      wireDir.normalize();

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

      Eigen::Vector3d docaVec = hitPosition - (wireStart + arcLen * wireDir);

      distance = docaVec.norm();
    }
    catch (std::exception& exc) {
      // This can happen, almost always because the coordinates are **just** out of range
      mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - "
                                  << exc.what() << std::endl;

      // Assume extremum for wire number depending on z coordinate
      distance = 0.;
    }

    return distance;
  }

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

private

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

Definition at line 1325 of file StandardHit3DBuilder_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::StandardHit3DBuilder::findGoodHitPairs ( const reco::ClusterHit2D *  goldenHit, HitVector::iterator &  startItr, HitVector::iterator &  endItr, HitMatchPairVecMap &  hitMatchMap  ) const

private

Definition at line 711 of file StandardHit3DBuilder_tool.cc.

References m_hitWidthSclFctr, makeHitPair(), and reco::ClusterHit2D::WireID().

Referenced by BuildHitPairMapByTPC().

  {
    int numPairs(0);

    // Loop through the input secon hits and make pairs
    while (startItr != endItr) {
      const reco::ClusterHit2D* hit = *startItr++;
      reco::ClusterHit3D pair;

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

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

      hitMatchMap[wireID].emplace_back(hit, pair);

      numPairs++;
    }

    return numPairs;
  }

void lar_cluster3d::StandardHit3DBuilder::findGoodTriplets ( HitMatchPairVecMap &  pair12Map, HitMatchPairVecMap &  pair13Map, reco::HitPairList &  hitPairList  ) const

private

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

Definition at line 736 of file StandardHit3DBuilder_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[&hit2Dhit3DPair.second] = false;
      }

      // The outer loop is over all hit pairs made from the first two plane combinations
      for (const auto& pair12 : pair12Map) {
        if (pair12.second.empty()) continue;

        // This loop is over hit pairs that share the same first two plane wires but may have different
        // hit times on those wires
        for (const auto& hit2Dhit3DPair12 : pair12.second) {
          const reco::ClusterHit3D& pair1 = hit2Dhit3DPair12.second;

          // populate the map with initial value
          usedPairMap[&pair1] = false;

          // The simplest approach here is to loop over all possibilities and let the triplet builder weed out the weak candidates
          for (const auto& pair13 : pair13Map) {
            if (pair13.second.empty()) continue;

            for (const auto& hit2Dhit3DPair13 : pair13.second) {
              const reco::ClusterHit2D* hit2 = hit2Dhit3DPair13.first;
              const reco::ClusterHit3D& pair2 = hit2Dhit3DPair13.second;

              // If success try for the triplet
              reco::ClusterHit3D triplet;

              if (makeHitTriplet(triplet, pair1, hit2)) {
                triplet.setID(hitPairList.size());
                hitPairList.emplace_back(triplet);
                usedPairMap[&pair1] = true;
                usedPairMap[&pair2] = true;
              }
            }
          }
        }
      }

      // One more loop through the other pairs to check for sick channels
      if (m_numBadChannels > 0) {
        for (const auto& pairMapPair : usedPairMap) {
          if (pairMapPair.second) continue;

          const reco::ClusterHit3D* pair = pairMapPair.first;

          // Here we look to see if we failed to make a triplet because the partner wire was dead/noisy/sick
          if (makeDeadChannelPair(deadChanPair, *pair, 4, 0))
            tempDeadChanVec.emplace_back(deadChanPair);
        }

        // Handle the dead wire triplets
        if (!tempDeadChanVec.empty()) {
          // If we have many then see if we can trim down a bit by keeping those with time significance
          if (tempDeadChanVec.size() > 1) {
            // Sort by "significance" of agreement
            std::sort(tempDeadChanVec.begin(),
                      tempDeadChanVec.end(),
                      [](const auto& left, const auto& right) {
                        return left.getDeltaPeakTime() / left.getSigmaPeakTime() <
                               right.getDeltaPeakTime() / right.getSigmaPeakTime();
                      });

            // What is the range of "significance" from first to last?
            float firstSig = tempDeadChanVec.front().getDeltaPeakTime() /
                             tempDeadChanVec.front().getSigmaPeakTime();
            float lastSig =
              tempDeadChanVec.back().getDeltaPeakTime() / tempDeadChanVec.back().getSigmaPeakTime();
            float sigRange = lastSig - firstSig;

            if (lastSig > 0.5 * m_deltaPeakTimeSig && sigRange > 0.5) {
              // Declare a maximum of 1.5 * the average of the first and last pairs...
              float maxSignificance = std::max(0.75 * (firstSig + lastSig), 1.0);

              std::vector<reco::ClusterHit3D>::iterator firstBadElem = std::find_if(
                tempDeadChanVec.begin(),
                tempDeadChanVec.end(),
                [&maxSignificance](const auto& pair) {
                  return pair.getDeltaPeakTime() / pair.getSigmaPeakTime() > maxSignificance;
                });

              // But only keep the best 10?
              if (std::distance(tempDeadChanVec.begin(), firstBadElem) > 20)
                firstBadElem = tempDeadChanVec.begin() + 20;
              // Keep at least one hit...
              else if (firstBadElem == tempDeadChanVec.begin())
                firstBadElem++;

              tempDeadChanVec.resize(std::distance(tempDeadChanVec.begin(), firstBadElem));
            }
          }

          for (auto& pair : tempDeadChanVec) {
            pair.setID(hitPairList.size());
            hitPairList.emplace_back(pair);
          }
        }
      }
    }

    return;
  }

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

private

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

Definition at line 1354 of file StandardHit3DBuilder_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;
  }

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

inlineoverridevirtual

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

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 102 of file StandardHit3DBuilder_tool.cc.

References tca::evt.

    {
      return m_timeVector[index];
    }

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

overridevirtual

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

Parameters

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

Associations with wires.

Associations with raw digits.

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 418 of file StandardHit3DBuilder_tool.cc.

References BuildHit3D(), clear(), CollectArtHits(), CreateNewRecobHitCollection(), m_clusterHit2DMasterList, m_outputHistograms, m_planeToHitVectorMap, 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_planeToHitVectorMap.clear();
    m_planeToWireToHitSetMap.clear();

    m_timeVector.resize(NUMTIMEVALUES, 0.);

    // Get a hit refiner
    std::unique_ptr<std::vector<recob::Hit>> outputHitPtrVec(new std::vector<recob::Hit>);

    // Recover the 2D hits and then organize them into data structures which will be used in the
    // DBscan algorithm for building the 3D clusters
    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
      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::StandardHit3DBuilder::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 1252 of file StandardHit3DBuilder_tool.cc.

References geo::CryostatID::Cryostat, reco::ClusterHit3D::getHits(), reco::ClusterHit3D::getPosition(), reco::ClusterHit2D::getStatusBits(), m_channelStatus, m_wireReadoutGeom, 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(), and geo::WireReadoutGeom::WireIDsIntersect().

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
      if (auto widIntersect0 = m_wireReadoutGeom->WireIDsIntersect(wireID0, wireID)) {
        if (auto widIntersect1 = m_wireReadoutGeom->WireIDsIntersect(wireID1, wireID)) {
          Eigen::Vector3f newPosition(
            pair.getPosition()[0], pair.getPosition()[1], pair.getPosition()[2]);

          newPosition[1] = (newPosition[1] + widIntersect0->y + widIntersect1->y) / 3.;
          newPosition[2] =
            (newPosition[2] + widIntersect0->z + widIntersect1->z - 2. * m_zPosOffset) / 3.;

          pairOut = pair;
          pairOut.setWireID(wireID);
          pairOut.setPosition(newPosition);

          if (hit0->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET)
            hit0->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);
          if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET)
            hit1->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET);

          hit0->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);
          hit1->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET);

          result = true;
        }
      }
    }

    return result;
  }

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

private

Make a HitPair object by checking two hits.

Definition at line 853 of file StandardHit3DBuilder_tool.cc.

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

Referenced by findGoodHitPairs(), and makeHitTriplet().

  {
    // Assume failure
    bool result(false);

    // We assume in this routine that we are looking at hits in different views
    // The first mission is to check that the wires intersect
    const geo::WireID& hit1WireID = hit1->WireID();
    const geo::WireID& hit2WireID = hit2->WireID();

    if (auto widIntersect = m_wireReadoutGeom->WireIDsIntersect(hit1WireID, hit2WireID)) {
      // Wires intersect so now we can check the timing
      float hit1Peak = hit1->getTimeTicks();
      float hit1Sigma = hit1->getHit()->RMS();

      float hit2Peak = hit2->getTimeTicks();
      float hit2Sigma = hit2->getHit()->RMS();

      // "Long hits" are an issue... so we deal with these differently
      int hit1NDF = hit1->getHit()->DegreesOfFreedom();
      int hit2NDF = hit2->getHit()->DegreesOfFreedom();

      // Basically, allow the range to extend to the nearest end of the snippet
      if (hit1NDF < 2)
        hit1Sigma = std::min(hit1Peak - float(hit1->getHit()->StartTick()),
                             float(hit1->getHit()->EndTick()) - hit1Peak);
      if (hit2NDF < 2)
        hit2Sigma = std::min(hit2Peak - float(hit2->getHit()->StartTick()),
                             float(hit2->getHit()->EndTick()) - hit2Peak);

      // The "hit sigma" is the gaussian fit sigma of the hit, we need to expand a bit to allow hit overlap efficiency
      float hit1Width = hitWidthSclFctr * hit1Sigma;
      float hit2Width = hitWidthSclFctr * hit2Sigma;

      // Coarse check hit times are "in range"
      if (fabs(hit1Peak - hit2Peak) <= (hit1Width + hit2Width)) {
        // Check to see that hit peak times are consistent with each other
        float hit1SigSq = hit1Sigma * hit1Sigma;
        float hit2SigSq = hit2Sigma * hit2Sigma;
        float deltaPeakTime = std::fabs(hit1Peak - hit2Peak);
        float sigmaPeakTime = std::sqrt(hit1SigSq + hit2SigSq);

        // delta peak time consistency check here
        if (deltaPeakTime < m_deltaPeakTimeSig *
                              sigmaPeakTime) // 2 sigma consistency? (do this way to avoid divide)
        {
          float oneOverWghts = hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
          float avePeakTime = (hit1Peak / hit1SigSq + hit2Peak / hit2SigSq) * oneOverWghts;
          float totalCharge = hit1->getHit()->Integral() + hit2->getHit()->Integral();
          float hitChiSquare = std::pow((hit1Peak - avePeakTime), 2) / hit1SigSq +
                               std::pow((hit2Peak - avePeakTime), 2) / hit2SigSq;

          float xPositionHit1(hit1->getXPosition());
          float xPositionHit2(hit2->getXPosition());
          float xPosition = (xPositionHit1 / hit1SigSq + xPositionHit2 / hit2SigSq) * hit1SigSq *
                            hit2SigSq / (hit1SigSq + hit2SigSq);

          Eigen::Vector3f position(
            xPosition, float(widIntersect->y), float(widIntersect->z) - m_zPosOffset);

          // If to here then we need to sort out the hit pair code telling what views are used
          unsigned statusBits = 1 << hit1->WireID().Plane | 1 << hit2->WireID().Plane;

          // handle status bits for the 2D hits
          if (hit1->getStatusBits() & reco::ClusterHit2D::USEDINPAIR)
            hit1->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);
          if (hit2->getStatusBits() & reco::ClusterHit2D::USEDINPAIR)
            hit2->setStatusBit(reco::ClusterHit2D::SHAREDINPAIR);

          hit1->setStatusBit(reco::ClusterHit2D::USEDINPAIR);
          hit2->setStatusBit(reco::ClusterHit2D::USEDINPAIR);

          reco::ClusterHit2DVec hitVector;

          hitVector.resize(3, NULL);

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

          unsigned int cryostatIdx = hit1->WireID().Cryostat;
          unsigned int tpcIdx = hit1->WireID().TPC;

          // Initialize the wireIdVec
          std::vector<geo::WireID> wireIDVec = {geo::WireID(cryostatIdx, tpcIdx, 0, 0),
                                                geo::WireID(cryostatIdx, tpcIdx, 1, 0),
                                                geo::WireID(cryostatIdx, tpcIdx, 2, 0)};

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

          // For compiling at the moment
          std::vector<float> hitDelTSigVec = {0., 0., 0.};

          hitDelTSigVec[hit1->WireID().Plane] = deltaPeakTime / sigmaPeakTime;
          hitDelTSigVec[hit2->WireID().Plane] = deltaPeakTime / sigmaPeakTime;

          // Create the 3D cluster hit
          hitPair.initialize(hitPairCntr,
                             statusBits,
                             position,
                             totalCharge,
                             avePeakTime,
                             deltaPeakTime,
                             sigmaPeakTime,
                             hitChiSquare,
                             0.,
                             0.,
                             0.,
                             0.,
                             hitVector,
                             hitDelTSigVec,
                             wireIDVec);

          result = true;
        }
      }
    }

    // Send it back
    return result;
  }

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

private

Make a 3D HitPair object by checking two hits.

Definition at line 979 of file StandardHit3DBuilder_tool.cc.

References util::abs(), chargeIntegral(), geo::CryostatID::Cryostat, recob::Hit::DegreesOfFreedom(), DistanceFromPointToHitWire(), recob::Hit::EndTick(), 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_maxDeltaPeakVec, m_maxPullVec, m_maxSideVecVec, m_outputHistograms, m_overlapFractionVec, m_overlapRangeVec, m_pairWireDistVec, m_qualityMetricVec, 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, recob::Hit::StartTick(), 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 = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                                float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());

            float weight = 1. / (hitRMS * hitRMS);

            avePeakTime += peakTime * weight;
            xPosition += hit2D->getXPosition() * weight;
            weightSum += weight;
          }

          avePeakTime /= weightSum;
          xPosition /= weightSum;

          Eigen::Vector2f pair0hYZVec(pair0h.getPosition()[1], pair0h.getPosition()[2]);
          Eigen::Vector2f pair1hYZVec(pair1h.getPosition()[1], pair1h.getPosition()[2]);
          Eigen::Vector2f pairYZVec(pair.getPosition()[1], pair.getPosition()[2]);
          Eigen::Vector3f position(xPosition,
                                   float((pairYZVec[0] + pair0hYZVec[0] + pair1hYZVec[0]) / 3.),
                                   float((pairYZVec[1] + pair0hYZVec[1] + pair1hYZVec[1]) / 3.));

          // Armed with the average peak time, now get hitChiSquare and the sig vec
          float hitChiSquare(0.);
          float sigmaPeakTime(std::sqrt(1. / weightSum));
          std::vector<float> hitDelTSigVec;

          for (const auto& hit2D : hitVector) {
            float hitRMS = hit2D->getHit()->RMS();

            // Basically, allow the range to extend to the nearest end of the snippet
            if (hit2D->getHit()->DegreesOfFreedom() < 2)
              hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                                float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());

            float combRMS = std::sqrt(hitRMS * hitRMS - sigmaPeakTime * sigmaPeakTime);
            float peakTime = hit2D->getTimeTicks();
            float deltaTime = peakTime - avePeakTime;
            float hitSig = deltaTime / combRMS;

            hitChiSquare += hitSig * hitSig;

            hitDelTSigVec.emplace_back(std::fabs(hitSig));
          }

          if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare);

          int lowMinIndex(std::numeric_limits<int>::max());
          int lowMaxIndex(std::numeric_limits<int>::min());
          int hiMinIndex(std::numeric_limits<int>::max());
          int hiMaxIndex(std::numeric_limits<int>::min());

          // First task is to get the min/max values for the common overlap region
          for (const auto& hit2D : hitVector) {
            float range = 2. * hit2D->getHit()->RMS();

            // Basically, allow the range to extend to the nearest end of the snippet
            if (hit2D->getHit()->DegreesOfFreedom() < 2)
              range = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),
                               float(hit2D->getHit()->EndTick()) - hit2D->getTimeTicks());

            int hitStart = hit2D->getHit()->PeakTime() - range - 0.5;
            int hitStop = hit2D->getHit()->PeakTime() + range + 0.5;

            lowMinIndex = std::min(hitStart, lowMinIndex);
            lowMaxIndex = std::max(hitStart, lowMaxIndex);
            hiMinIndex = std::min(hitStop + 1, hiMinIndex);
            hiMaxIndex = std::max(hitStop + 1, hiMaxIndex);
          }

          // Keep only "good" hits...
          if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex) {
            // One more pass through hits to get charge
            std::vector<float> chargeVec;

            for (const auto& hit2D : hitVector)
              chargeVec.push_back(chargeIntegral(hit2D->getHit()->PeakTime(),
                                                 hit2D->getHit()->PeakAmplitude(),
                                                 hit2D->getHit()->RMS(),
                                                 lowMaxIndex,
                                                 hiMinIndex));

            float totalCharge =
              std::accumulate(chargeVec.begin(), chargeVec.end(), 0.) / float(chargeVec.size());
            float overlapRange = float(hiMinIndex - lowMaxIndex);
            float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex);

            // Set up to compute the charge asymmetry
            std::vector<float> smallestChargeDiffVec;
            std::vector<float> chargeAveVec;
            float smallestDiff(std::numeric_limits<float>::max());
            float maxDeltaPeak(0.);
            size_t chargeIndex(0);

            for (size_t idx = 0; idx < 3; idx++) {
              size_t leftIdx = (idx + 2) % 3;
              size_t rightIdx = (idx + 1) % 3;

              smallestChargeDiffVec.push_back(std::abs(chargeVec[leftIdx] - chargeVec[rightIdx]));
              chargeAveVec.push_back(float(0.5 * (chargeVec[leftIdx] + chargeVec[rightIdx])));

              if (smallestChargeDiffVec.back() < smallestDiff) {
                smallestDiff = smallestChargeDiffVec.back();
                chargeIndex = idx;
              }

              // Take opportunity to look at peak time diff
              float deltaPeakTime =
                hitVector[leftIdx]->getTimeTicks() - hitVector[rightIdx]->getTimeTicks();

              if (std::abs(deltaPeakTime) > maxDeltaPeak) maxDeltaPeak = std::abs(deltaPeakTime);

              if (m_outputHistograms) m_deltaTimeVec.push_back(deltaPeakTime);
            }

            float chargeAsymmetry = (chargeAveVec[chargeIndex] - chargeVec[chargeIndex]) /
                                    (chargeAveVec[chargeIndex] + chargeVec[chargeIndex]);

            // If this is true there has to be a negative charge that snuck in somehow
            if (chargeAsymmetry < -1. || chargeAsymmetry > 1.) {
              const geo::WireID& hitWireID = hitVector[chargeIndex]->WireID();

              std::cout << "============> Charge asymmetry out of range: " << chargeAsymmetry
                        << " <============" << std::endl;
              std::cout << "     hit C: " << hitWireID.Cryostat << ", TPC: " << hitWireID.TPC
                        << ", Plane: " << hitWireID.Plane << ", Wire: " << hitWireID.Wire
                        << std::endl;
              std::cout << "     charge: " << chargeVec[0] << ", " << chargeVec[1] << ", "
                        << chargeVec[2] << std::endl;
              std::cout << "     index: " << chargeIndex << ", smallest diff: " << smallestDiff
                        << std::endl;
              return result;
            }

            // Usurping "deltaPeakTime" to be the maximum pull
            float deltaPeakTime = *std::max_element(hitDelTSigVec.begin(), hitDelTSigVec.end());

            if (m_outputHistograms) {
              m_smallChargeDiffVec.push_back(smallestDiff);
              m_smallIndexVec.push_back(chargeIndex);
              m_maxPullVec.push_back(deltaPeakTime);
              m_qualityMetricVec.push_back(hitChiSquare);
              m_spacePointChargeVec.push_back(totalCharge);
              m_overlapFractionVec.push_back(overlapFraction);
              m_overlapRangeVec.push_back(overlapRange);
              m_maxDeltaPeakVec.push_back(maxDeltaPeak);
              m_hitAsymmetryVec.push_back(chargeAsymmetry);
            }

            // Try to weed out cases where overlap doesn't match peak separation
            if (maxDeltaPeak > overlapRange) return result;

            // Create the 3D cluster hit
            hitTriplet.initialize(0,
                                  statusBits,
                                  position,
                                  totalCharge,
                                  avePeakTime,
                                  deltaPeakTime,
                                  sigmaPeakTime,
                                  hitChiSquare,
                                  overlapFraction,
                                  chargeAsymmetry,
                                  0.,
                                  0.,
                                  hitVector,
                                  hitDelTSigVec,
                                  wireIDVec);

            result = true;
          }
        }
      }
    }

    // return success/fail
    return result;
  }

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

private

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

Definition at line 1708 of file StandardHit3DBuilder_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::StandardHit3DBuilder::makeWireAssns ( const art::Event &  evt, art::Assns< recob::Wire, recob::Hit > &  wireAssns, RecobHitToPtrMap &  recobHitPtrMap  ) const

private

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

Definition at line 1662 of file StandardHit3DBuilder_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::StandardHit3DBuilder::NearestWireID ( const Eigen::Vector3f &  position, const geo::WireID &  wireID  ) const

private

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

Definition at line 1380 of file StandardHit3DBuilder_tool.cc.

References geo::WireID::asPlaneID(), geo::TPCGeo::Length(), m_geometry, m_wireReadoutGeom, geo::WireReadoutGeom::Nwires(), geo::WireReadoutGeom::Plane(), geo::vect::toPoint(), geo::GeometryCore::TPC(), 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_wireReadoutGeom->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->TPC().Length())
        wireID.Wire = 0;
      else
        wireID.Wire = m_wireReadoutGeom->Nwires(wireIDIn.asPlaneID()) - 1;
    }

    return wireID;
  }

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

overridevirtual

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

Implements lar_cluster3d::IHit3DBuilder.

Definition at line 347 of file StandardHit3DBuilder_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>>();
  }

Member Data Documentation

const lariov::ChannelStatusProvider* lar_cluster3d::StandardHit3DBuilder::m_channelFilter

private

Definition at line 292 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and StandardHit3DBuilder().

ChannelStatusByPlaneVec lar_cluster3d::StandardHit3DBuilder::m_channelStatus

mutableprivate

Definition at line 285 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and makeDeadChannelPair().

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

mutableprivate

Definition at line 267 of file StandardHit3DBuilder_tool.cc.

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

Hit2DList lar_cluster3d::StandardHit3DBuilder::m_clusterHit2DMasterList

mutableprivate

Definition at line 281 of file StandardHit3DBuilder_tool.cc.

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

float lar_cluster3d::StandardHit3DBuilder::m_deltaPeakTimeSig

private

Definition at line 250 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 266 of file StandardHit3DBuilder_tool.cc.

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

bool lar_cluster3d::StandardHit3DBuilder::m_enableMonitoring

private

Definition at line 257 of file StandardHit3DBuilder_tool.cc.

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

const geo::Geometry* lar_cluster3d::StandardHit3DBuilder::m_geometry

private

Definition at line 290 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildHitPairMap(), CollectArtHits(), NearestWireID(), and StandardHit3DBuilder().

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

mutableprivate

Definition at line 278 of file StandardHit3DBuilder_tool.cc.

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

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

private

Data members to follow.

Definition at line 247 of file StandardHit3DBuilder_tool.cc.

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

float lar_cluster3d::StandardHit3DBuilder::m_hitWidthSclFctr

private

Definition at line 249 of file StandardHit3DBuilder_tool.cc.

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

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

private

Definition at line 251 of file StandardHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and StandardHit3DBuilder().

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

mutableprivate

Definition at line 271 of file StandardHit3DBuilder_tool.cc.

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

float lar_cluster3d::StandardHit3DBuilder::m_maxHit3DChiSquare

private

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

Definition at line 254 of file StandardHit3DBuilder_tool.cc.

Referenced by StandardHit3DBuilder().

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

mutableprivate

Definition at line 268 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 272 of file StandardHit3DBuilder_tool.cc.

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

size_t lar_cluster3d::StandardHit3DBuilder::m_numBadChannels

mutableprivate

Definition at line 286 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), and findGoodTriplets().

float lar_cluster3d::StandardHit3DBuilder::m_numSigmaPeakTime

private

Definition at line 248 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildHitPairMap(), BuildHitPairMapByTPC(), and StandardHit3DBuilder().

bool lar_cluster3d::StandardHit3DBuilder::m_outputHistograms

private

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

Definition at line 255 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 269 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 270 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 273 of file StandardHit3DBuilder_tool.cc.

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

PlaneToHitVectorMap lar_cluster3d::StandardHit3DBuilder::m_planeToHitVectorMap

mutableprivate

Definition at line 282 of file StandardHit3DBuilder_tool.cc.

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

PlaneToWireToHitSetMap lar_cluster3d::StandardHit3DBuilder::m_planeToWireToHitSetMap

mutableprivate

Definition at line 283 of file StandardHit3DBuilder_tool.cc.

Referenced by CollectArtHits(), and Hit3DBuilder().

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

mutableprivate

Definition at line 276 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 274 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 275 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 277 of file StandardHit3DBuilder_tool.cc.

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

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

mutableprivate

Definition at line 259 of file StandardHit3DBuilder_tool.cc.

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

TTree* lar_cluster3d::StandardHit3DBuilder::m_tupleTree

private

output analysis tree

Definition at line 264 of file StandardHit3DBuilder_tool.cc.

Referenced by Hit3DBuilder(), and StandardHit3DBuilder().

bool lar_cluster3d::StandardHit3DBuilder::m_weHaveAllBeenHereBefore = false

mutableprivate

Definition at line 288 of file StandardHit3DBuilder_tool.cc.

Referenced by CollectArtHits().

float lar_cluster3d::StandardHit3DBuilder::m_wirePitch[3]

private

Definition at line 258 of file StandardHit3DBuilder_tool.cc.

Referenced by makeHitTriplet(), and StandardHit3DBuilder().

float lar_cluster3d::StandardHit3DBuilder::m_wirePitchScaleFactor

private

Scaling factor to determine max distance allowed between candidate pairs.

Definition at line 253 of file StandardHit3DBuilder_tool.cc.

Referenced by makeHitTriplet(), and StandardHit3DBuilder().

const geo::WireReadoutGeom* lar_cluster3d::StandardHit3DBuilder::m_wireReadoutGeom

private

Definition at line 291 of file StandardHit3DBuilder_tool.cc.

Referenced by BuildChannelStatusVec(), CollectArtHits(), DistanceFromPointToHitWire(), makeDeadChannelPair(), makeHitPair(), NearestWireID(), and StandardHit3DBuilder().

float lar_cluster3d::StandardHit3DBuilder::m_zPosOffset

private

Definition at line 261 of file StandardHit3DBuilder_tool.cc.

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

---

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

• larreco/v10_01_07/source/larreco/RecoAlg/Cluster3DAlgs/StandardHit3DBuilder_tool.cc