StandardHit3DBuilder_tool.cc

Go to the documentation of this file.

// Framework Includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Persistency/Common/PtrMaker.h"
#include "art/Utilities/ToolMacros.h"
#include "art_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/InputTag.h"
#include "cetlib/cpu_timer.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "lardata/ArtDataHelper/HitCreator.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
#include "larreco/RecoAlg/Cluster3DAlgs/IHit3DBuilder.h"

// Eigen
#include <Eigen/Core>

// std includes
#include <iostream>
#include <memory>
#include <numeric> // std::accumulate
#include <string>

// Ack!
#include "TH1F.h"
#include "TTree.h"

//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows

namespace lar_cluster3d {

  // forward declaration to define an ordering function for our hit set
  struct Hit2DSetCompare {
    bool operator()(const reco::ClusterHit2D*, const reco::ClusterHit2D*) const;
  };

  using HitVector = std::vector<const reco::ClusterHit2D*>;
  using PlaneToHitVectorMap = std::map<geo::PlaneID, HitVector>;
  using TPCToPlaneToHitVectorMap = std::map<geo::TPCID, PlaneToHitVectorMap>;
  using Hit2DList = std::list<reco::ClusterHit2D>;
  using Hit2DSet = std::set<const reco::ClusterHit2D*, Hit2DSetCompare>;
  using WireToHitSetMap = std::map<unsigned int, Hit2DSet>;
  using PlaneToWireToHitSetMap = std::map<geo::PlaneID, WireToHitSetMap>;
  using TPCToPlaneToWireToHitSetMap = std::map<geo::TPCID, PlaneToWireToHitSetMap>;
  using HitVectorMap = std::map<size_t, HitVector>;

  class StandardHit3DBuilder : public IHit3DBuilder {
  public:
    explicit StandardHit3DBuilder(fhicl::ParameterSet const& pset);

    void produces(art::ProducesCollector&) override;

    void Hit3DBuilder(art::Event&, reco::HitPairList&, RecobHitToPtrMap&) override;

    float getTimeToExecute(IHit3DBuilder::TimeValues index) const override
    {
      return m_timeVector[index];
    }

  private:
    void CollectArtHits(const art::Event& evt) const;

    void BuildHit3D(reco::HitPairList& hitPairList) const;

    void CreateNewRecobHitCollection(art::Event&,
                                     reco::HitPairList&,
                                     std::vector<recob::Hit>&,
                                     RecobHitToPtrMap&);

    void makeWireAssns(const art::Event&,
                       art::Assns<recob::Wire, recob::Hit>&,
                       RecobHitToPtrMap&) const;

    void makeRawDigitAssns(const art::Event&,
                           art::Assns<raw::RawDigit, recob::Hit>&,
                           RecobHitToPtrMap&) const;

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

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

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

    using HitMatchPair = std::pair<const reco::ClusterHit2D*, reco::ClusterHit3D>;
    using HitMatchPairVec = std::vector<HitMatchPair>;
    using HitMatchPairVecMap = std::map<geo::WireID, HitMatchPairVec>;

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

    void findGoodTriplets(HitMatchPairVecMap&, HitMatchPairVecMap&, reco::HitPairList&) const;

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

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

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

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

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

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

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

    void BuildChannelStatusVec() const;

    float chargeIntegral(float, float, float, int, int) const;

    using ChannelStatusVec = std::vector<size_t>;
    using ChannelStatusByPlaneVec = std::vector<ChannelStatusVec>;

    void clear();

    std::vector<art::InputTag> m_hitFinderTagVec;
    float m_numSigmaPeakTime;
    float m_hitWidthSclFctr;
    float m_deltaPeakTimeSig;
    std::vector<int> m_invalidTPCVec;
    float
      m_wirePitchScaleFactor; 
    float m_maxHit3DChiSquare; 
    bool m_outputHistograms;   

    bool m_enableMonitoring; 
    float m_wirePitch[3];
    mutable std::vector<float> m_timeVector; 

    float m_zPosOffset;

    // Define some basic histograms
    TTree* m_tupleTree; 

    mutable std::vector<float> m_deltaTimeVec;
    mutable std::vector<float> m_chiSquare3DVec;
    mutable std::vector<float> m_maxPullVec;
    mutable std::vector<float> m_overlapFractionVec;
    mutable std::vector<float> m_overlapRangeVec;
    mutable std::vector<float> m_maxDeltaPeakVec;
    mutable std::vector<float> m_maxSideVecVec;
    mutable std::vector<float> m_pairWireDistVec;
    mutable std::vector<float> m_smallChargeDiffVec;
    mutable std::vector<int> m_smallIndexVec;
    mutable std::vector<float> m_qualityMetricVec;
    mutable std::vector<float> m_spacePointChargeVec;
    mutable std::vector<float> m_hitAsymmetryVec;

    // Get instances of the primary data structures needed
    mutable Hit2DList m_clusterHit2DMasterList;
    mutable PlaneToHitVectorMap m_planeToHitVectorMap;
    mutable PlaneToWireToHitSetMap m_planeToWireToHitSetMap;

    mutable ChannelStatusByPlaneVec m_channelStatus;
    mutable size_t m_numBadChannels;

    mutable bool m_weHaveAllBeenHereBefore = false;

    const geo::Geometry* m_geometry;
    const geo::WireReadoutGeom* m_wireReadoutGeom;
    const lariov::ChannelStatusProvider* m_channelFilter;
  };

  StandardHit3DBuilder::StandardHit3DBuilder(fhicl::ParameterSet const& pset)
    : 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);
    }
  }

  //------------------------------------------------------------------------------------------------------------------------------------------

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

  //------------------------------------------------------------------------------------------------------------------------------------------

  void StandardHit3DBuilder::clear()
  {
    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 StandardHit3DBuilder::BuildChannelStatusVec() const
  {
    // 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;
    }
  }

  bool SetPeakHitPairIteratorOrder(const reco::HitPairList::iterator& left,
                                   const reco::HitPairList::iterator& right)
  {
    return (*left).getAvePeakTime() < (*right).getAvePeakTime();
  }

  struct HitPairClusterOrder {
    bool operator()(const reco::HitPairClusterMap::iterator& left,
                    const reco::HitPairClusterMap::iterator& right)
    {
      // Watch out for the case where two clusters can have the same number of hits!
      if (left->second.size() == right->second.size()) return left->first < right->first;

      return left->second.size() > right->second.size();
    }
  };

  void StandardHit3DBuilder::Hit3DBuilder(art::Event& evt,
                                          reco::HitPairList& hitPairList,
                                          RecobHitToPtrMap& clusterHitToArtPtrMap)
  {
    // 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();
    }
  }

  void StandardHit3DBuilder::BuildHit3D(reco::HitPairList& hitPairList) const
  {
    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;
  }

  //------------------------------------------------------------------------------------------------------------------------------------------
  namespace {

    class SetHitEarliestTimeOrder {
    public:
      SetHitEarliestTimeOrder() : m_numRMS(1.) {}
      SetHitEarliestTimeOrder(float numRMS) : m_numRMS(numRMS) {}

      bool operator()(const reco::ClusterHit2D* left, const reco::ClusterHit2D* right) const
      {
        return left->getTimeTicks() - m_numRMS * left->getHit()->RMS() <
               right->getTimeTicks() - m_numRMS * right->getHit()->RMS();
      }

    public:
      float m_numRMS;
    };

    using HitVectorItrPair = std::pair<HitVector::iterator, HitVector::iterator>;

    class SetStartTimeOrder {
    public:
      SetStartTimeOrder() : m_numRMS(1.) {}
      SetStartTimeOrder(float numRMS) : m_numRMS(numRMS) {}

      bool operator()(const HitVectorItrPair& left, const HitVectorItrPair& right) const
      {
        // Protect against possible issue?
        if (left.first != left.second && right.first != right.second) {
          // Sort by "modified start time" of pulse
          return (*left.first)->getTimeTicks() - m_numRMS * (*left.first)->getHit()->RMS() <
                 (*right.first)->getTimeTicks() - m_numRMS * (*right.first)->getHit()->RMS();
        }

        return left.first != left.second;
      }

    private:
      float m_numRMS;
    };

    bool SetPairStartTimeOrder(const reco::ClusterHit3D& left, const reco::ClusterHit3D& right)
    {
      // Sort by "modified start time" of pulse
      return left.getAvePeakTime() - left.getSigmaPeakTime() <
             right.getAvePeakTime() - right.getSigmaPeakTime();
    }
  }

  //------------------------------------------------------------------------------------------------------------------------------------------

  size_t StandardHit3DBuilder::BuildHitPairMap(PlaneToHitVectorMap& planeToHitVectorMap,
                                               reco::HitPairList& hitPairList) const
  {
    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 StandardHit3DBuilder::BuildHitPairMapByTPC(PlaneHitVectorItrPairVec& hitItrVec,
                                                    reco::HitPairList& hitPairList) const
  {
    // 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();
  }

  int StandardHit3DBuilder::findGoodHitPairs(const reco::ClusterHit2D* goldenHit,
                                             HitVector::iterator& startItr,
                                             HitVector::iterator& endItr,
                                             HitMatchPairVecMap& hitMatchMap) const
  {
    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 StandardHit3DBuilder::findGoodTriplets(HitMatchPairVecMap& pair12Map,
                                              HitMatchPairVecMap& pair13Map,
                                              reco::HitPairList& hitPairList) const
  {
    // 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;
  }

  bool StandardHit3DBuilder::makeHitPair(reco::ClusterHit3D& hitPair,
                                         const reco::ClusterHit2D* hit1,
                                         const reco::ClusterHit2D* hit2,
                                         float hitWidthSclFctr,
                                         size_t hitPairCntr) const
  {
    // 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 StandardHit3DBuilder::makeHitTriplet(reco::ClusterHit3D& hitTriplet,
                                            const reco::ClusterHit3D& pair,
                                            const reco::ClusterHit2D* hit) const
  {
    // 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;
  }

  float StandardHit3DBuilder::chargeIntegral(float peakMean,
                                             float peakAmp,
                                             float peakSigma,
                                             int low,
                                             int hi) const
  {
    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;
  }

  bool StandardHit3DBuilder::makeDeadChannelPair(reco::ClusterHit3D& pairOut,
                                                 const reco::ClusterHit3D& pair,
                                                 size_t maxChanStatus,
                                                 size_t minChanStatus) const
  {
    // 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;
  }

  const reco::ClusterHit2D* StandardHit3DBuilder::FindBestMatchingHit(
    const Hit2DSet& hit2DSet,
    const reco::ClusterHit3D& pair,
    float pairDeltaTimeLimits) const
  {
    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 StandardHit3DBuilder::FindNumberInRange(const Hit2DSet& hit2DSet,
                                              const reco::ClusterHit3D& pair,
                                              float range) const
  {
    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;
  }

  geo::WireID StandardHit3DBuilder::NearestWireID(const Eigen::Vector3f& position,
                                                  const geo::WireID& wireIDIn) const
  {
    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;
  }

  float StandardHit3DBuilder::DistanceFromPointToHitWire(const Eigen::Vector3f& position,
                                                         const geo::WireID& wireIDIn) const
  {
    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;
  }

  //------------------------------------------------------------------------------------------------------------------------------------------
  bool SetHitTimeOrder(const reco::ClusterHit2D* left, const reco::ClusterHit2D* right)
  {
    // Sort by "modified start time" of pulse
    return left->getHit()->PeakTime() < right->getHit()->PeakTime();
  }

  bool Hit2DSetCompare::operator()(const reco::ClusterHit2D* left,
                                   const reco::ClusterHit2D* right) const
  {
    return left->getHit()->PeakTime() < right->getHit()->PeakTime();
  }

  //------------------------------------------------------------------------------------------------------------------------------------------
  void StandardHit3DBuilder::CollectArtHits(const art::Event& evt) const
  {
    // 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 StandardHit3DBuilder::CreateNewRecobHitCollection(art::Event& event,
                                                         reco::HitPairList& hitPairList,
                                                         std::vector<recob::Hit>& hitPtrVec,
                                                         RecobHitToPtrMap& recobHitToPtrMap)
  {
    // 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;
  }

  void StandardHit3DBuilder::makeWireAssns(const art::Event& evt,
                                           art::Assns<recob::Wire, recob::Hit>& wireAssns,
                                           RecobHitToPtrMap& recobHitPtrMap) const
  {
    // 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;
  }

  void StandardHit3DBuilder::makeRawDigitAssns(const art::Event& evt,
                                               art::Assns<raw::RawDigit, recob::Hit>& rawDigitAssns,
                                               RecobHitToPtrMap& recobHitPtrMap) const
  {
    // 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;
  }

  //------------------------------------------------------------------------------------------------------------------------------------------

  DEFINE_ART_CLASS_TOOL(StandardHit3DBuilder)
} // namespace lar_cluster3d