SnippetHit3DBuilder_tool.cc

Go to the documentation of this file.

// Framework Includes
#include "art/Framework/Core/ProducesCollector.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/DetectorClocksService.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 HitStartEndPair = std::pair<raw::TDCtick_t, raw::TDCtick_t>;
  using SnippetHitMap = std::map<HitStartEndPair, HitVector>;
  using PlaneToSnippetHitMap = std::map<geo::PlaneID, SnippetHitMap>;
  using TPCToPlaneToSnippetHitMap = std::map<geo::TPCID, PlaneToSnippetHitMap>;
  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>;
  using SnippetHitMapItrPair = std::pair<SnippetHitMap::iterator, SnippetHitMap::iterator>;
  using PlaneSnippetHitMapItrPairVec = std::vector<SnippetHitMapItrPair>;

  class SnippetHit3DBuilder : public IHit3DBuilder {
  public:
    explicit SnippetHit3DBuilder(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(PlaneToSnippetHitMap& planeToHitVectorMap,
                           reco::HitPairList& hitPairList) const;

    size_t BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec& planeSnippetHitMapItrPairVec,
                                reco::HitPairList& hitPairList) const;

    using HitMatchTriplet =
      std::tuple<const reco::ClusterHit2D*, const reco::ClusterHit2D*, const reco::ClusterHit3D>;
    using HitMatchTripletVec = std::vector<HitMatchTriplet>;
    using HitMatchTripletVecMap = std::map<geo::WireID, HitMatchTripletVec>;

    int findGoodHitPairs(SnippetHitMap::iterator&,
                         SnippetHitMap::iterator&,
                         SnippetHitMap::iterator&,
                         HitMatchTripletVecMap&) const;

    void findGoodTriplets(HitMatchTripletVecMap&, HitMatchTripletVecMap&, 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;

    bool WireIDsIntersect(const geo::WireID&, const geo::WireID&, geo::WireIDIntersection&) const;

    float closestApproach(const Eigen::Vector3f&,
                          const Eigen::Vector3f&,
                          const Eigen::Vector3f&,
                          const Eigen::Vector3f&,
                          float&,
                          float&) 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_hitWidthSclFctr;
    float m_deltaPeakTimeSig;
    float m_rangeNumSig;
    float m_LongHitStretchFctr;
    float m_pulseHeightFrac;
    float m_PHLowSelection;
    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 PlaneToSnippetHitMap m_planeToSnippetHitMap;
    mutable PlaneToWireToHitSetMap m_planeToWireToHitSetMap;

    mutable ChannelStatusByPlaneVec m_channelStatus;
    mutable size_t m_numBadChannels;

    mutable bool m_weHaveAllBeenHereBefore = false;

    const geo::Geometry* m_geometry; //< pointer to the Geometry service
    const geo::WireReadoutGeom* m_wireReadoutGeom;
    const lariov::ChannelStatusProvider* m_channelFilter;
  };

  SnippetHit3DBuilder::SnippetHit3DBuilder(fhicl::ParameterSet const& pset)
    : 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_hitWidthSclFctr = pset.get<float>("HitWidthScaleFactor", 6.);
    m_rangeNumSig = pset.get<float>("RangeNumSigma", 3.);
    m_LongHitStretchFctr = pset.get<float>("LongHitsStretchFactor", 1.5);
    m_pulseHeightFrac = pset.get<float>("PulseHeightFraction", 0.5);
    m_PHLowSelection = pset.get<float>("PHLowSelection", 20.);
    m_deltaPeakTimeSig = pset.get<float>("DeltaPeakTimeSig", 1.7);
    m_zPosOffset = pset.get<float>("ZPosOffset", 0.0);
    m_invalidTPCVec = pset.get<std::vector<int>>("InvalidTPCVec", std::vector<int>());
    m_wirePitchScaleFactor = pset.get<float>("WirePitchScaleFactor", 1.9);
    m_maxHit3DChiSquare = pset.get<float>("MaxHitChiSquare", 6.0);
    m_outputHistograms = pset.get<bool>("OutputHistograms", false);

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

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

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

      clear();

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

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

  void SnippetHit3DBuilder::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 SnippetHit3DBuilder::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 SnippetHit3DBuilder::BuildChannelStatusVec() const
  {
    // This is called each event, clear out the previous version and start over
    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 SnippetHit3DBuilder::Hit3DBuilder(art::Event& evt,
                                         reco::HitPairList& hitPairList,
                                         RecobHitToPtrMap& clusterHitToArtPtrMap)
  {
    // Clear the internal data structures
    m_clusterHit2DMasterList.clear();
    m_planeToSnippetHitMap.clear();
    m_planeToWireToHitSetMap.clear();

    m_timeVector.resize(NUMTIMEVALUES, 0.);

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

    // Recover the 2D hits and then organize them into data structures which will be used in the
    // DBscan algorithm for building the 3D clusters
    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 SnippetHit3DBuilder::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_planeToSnippetHitMap, hitPairList);

    if (m_enableMonitoring) {
      theClockMakeHits.stop();

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

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

  //------------------------------------------------------------------------------------------------------------------------------------------
  struct SetStartTimeOrder {
    bool operator()(const SnippetHitMapItrPair& left, const SnippetHitMapItrPair& right) const
    {
      // Special case handling, there is nothing to compare for the left or right
      if (left.first == left.second) return false;
      if (right.first == right.second) return true;

      // de-referencing a bunch of pairs here...
      return left.first->first.first < right.first->first.first;
    }
  };

  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 SnippetHit3DBuilder::BuildHitPairMap(PlaneToSnippetHitMap& planeToSnippetHitMap,
                                              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++) {
        PlaneToSnippetHitMap::iterator mapItr0 =
          planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0));
        PlaneToSnippetHitMap::iterator mapItr1 =
          planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1));
        PlaneToSnippetHitMap::iterator mapItr2 =
          planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2));

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

        if (nPlanesWithHits < 2) continue;

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

        PlaneSnippetHitMapItrPairVec hitItrVec = {
          SnippetHitMapItrPair(snippetHitMap0.begin(), snippetHitMap0.end()),
          SnippetHitMapItrPair(snippetHitMap1.begin(), snippetHitMap1.end()),
          SnippetHitMapItrPair(snippetHitMap2.begin(), snippetHitMap2.end())};

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

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

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

    return hitPairList.size();
  }

  size_t SnippetHit3DBuilder::BuildHitPairMapByTPC(
    PlaneSnippetHitMapItrPairVec& snippetHitMapItrVec,
    reco::HitPairList& hitPairList) const
  {
    // Define functions to set start/end iterators in the loop below
    auto SetStartIterator =
      [](SnippetHitMap::iterator startItr, SnippetHitMap::iterator endItr, float startTime) {
        while (startItr != endItr) {
          if (startItr->first.second < startTime)
            startItr++;
          else
            break;
        }
        return startItr;
      };

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

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

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

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

      if (nPlanesWithHits < 2) break;

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

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

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

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

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

      snippetHitMapItrVec.front().first++;
    }

    return hitPairList.size();
  }

  int SnippetHit3DBuilder::findGoodHitPairs(SnippetHitMap::iterator& firstSnippetItr,
                                            SnippetHitMap::iterator& startItr,
                                            SnippetHitMap::iterator& endItr,
                                            HitMatchTripletVecMap& hitMatchMap) const
  {
    int numPairs(0);

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

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

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

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

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

          reco::ClusterHit3D pair;

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

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

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

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

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

          numPairs++;
        }

        secondHitItr++;
      }
    }

    return numPairs;
  }

  void SnippetHit3DBuilder::findGoodTriplets(HitMatchTripletVecMap& pair12Map,
                                             HitMatchTripletVecMap& 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[&std::get<2>(hit2Dhit3DPair)] = false;
      }

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

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

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

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

            for (const auto& hit2Dhit3DPair13 : pair13.second) {
              // Protect against double counting
              if (std::get<0>(hit2Dhit3DPair12) != std::get<0>(hit2Dhit3DPair13)) continue;

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

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

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

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

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

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

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

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

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

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

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

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

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

    return;
  }

  bool SnippetHit3DBuilder::makeHitPair(reco::ClusterHit3D& hitPair,
                                        const reco::ClusterHit2D* hit1,
                                        const reco::ClusterHit2D* hit2,
                                        float hitWidthSclFctr,
                                        size_t hitPairCntr) const
  {
    // Assume failure
    bool result(false);

    // Start by checking time consistency since this is fastest
    // Wires intersect so now we can check the timing
    float hit1Peak = hit1->getTimeTicks();
    float hit1Sigma = hit1->getHit()->RMS();

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

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

    // Basically, allow the range to extend to the nearest end of the snippet
    if (hit1NDF < 2)
      hit1Sigma *= m_LongHitStretchFctr; // This sets the range to the width of the pulse
    if (hit2NDF < 2) hit2Sigma *= m_LongHitStretchFctr;

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

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

      // delta peak time consistency check here
      if (deltaPeakTime <
          m_deltaPeakTimeSig * sigmaPeakTime) // 2 sigma consistency? (do this way to avoid divide)
      {
        // We assume in this routine that we are looking at hits in different views
        // The first mission is to check that the wires intersect
        const geo::WireID& hit1WireID = hit1->WireID();
        const geo::WireID& hit2WireID = hit2->WireID();

        geo::WireIDIntersection widIntersect;

        if (WireIDsIntersect(hit1WireID, hit2WireID, widIntersect)) {
          float oneOverWghts = hit1SigSq * hit2SigSq / (hit1SigSq + hit2SigSq);
          float avePeakTime = (hit1Peak / hit1SigSq + hit2Peak / hit2SigSq) * oneOverWghts;
          float totalCharge = hit1->getHit()->Integral() + hit2->getHit()->Integral();
          float hitChiSquare = std::pow((hit1Peak - avePeakTime), 2) / hit1SigSq +
                               std::pow((hit2Peak - avePeakTime), 2) / hit2SigSq;

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

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

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

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

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

          reco::ClusterHit2DVec hitVector;

          hitVector.resize(3, NULL);

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

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

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

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

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

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

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

          result = true;
        }
      }
    }

    // Send it back
    return result;
  }

  bool SnippetHit3DBuilder::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 *= m_LongHitStretchFctr;

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

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

          avePeakTime /= weightSum;
          xPosition /= weightSum;

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

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

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

            // Basically, allow the range to extend to the nearest end of the snippet
            if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr;

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

            hitChiSquare += hitSig * hitSig;

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

          if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare);

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

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

            // Basically, allow the range to extend to the nearest end of the snippet
            if (hit2D->getHit()->DegreesOfFreedom() < 2) range *= m_LongHitStretchFctr;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

            result = true;
          }
        }
      }
    }

    // return success/fail
    return result;
  }

  bool SnippetHit3DBuilder::WireIDsIntersect(const geo::WireID& wireID0,
                                             const geo::WireID& wireID1,
                                             geo::WireIDIntersection& widIntersection) const
  {
    bool success(false);

    // Do quick check that things are in the same logical TPC
    if (wireID0.Cryostat != wireID1.Cryostat || wireID0.TPC != wireID1.TPC ||
        wireID0.Plane == wireID1.Plane)
      return success;

    // Recover wire geometry information for each wire
    const geo::WireGeo& wireGeo0 = m_wireReadoutGeom->Wire(wireID0);
    const geo::WireGeo& wireGeo1 = m_wireReadoutGeom->Wire(wireID1);

    // Get wire position and direction for first wire
    auto wirePosArr = wireGeo0.GetCenter();

    Eigen::Vector3f wirePos0(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
    Eigen::Vector3f wireDir0(
      wireGeo0.Direction().X(), wireGeo0.Direction().Y(), wireGeo0.Direction().Z());

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

    Eigen::Vector3f wirePos1(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
    Eigen::Vector3f wireDir1(
      wireGeo1.Direction().X(), wireGeo1.Direction().Y(), wireGeo1.Direction().Z());

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

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

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

        success = true;
      }
    }

    return success;
  }

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

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

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

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

  float SnippetHit3DBuilder::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 SnippetHit3DBuilder::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* SnippetHit3DBuilder::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 SnippetHit3DBuilder::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 SnippetHit3DBuilder::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 SnippetHit3DBuilder::DistanceFromPointToHitWire(const Eigen::Vector3f& position,
                                                        const geo::WireID& wireIDIn) const
  {
    float distance = std::numeric_limits<float>::max();

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

      // Get wire position and direction for first wire
      auto const wirePosArr = wireGeo.GetCenter();

      Eigen::Vector3f wirePos(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z());
      Eigen::Vector3f wireDir(
        wireGeo.Direction().X(), wireGeo.Direction().Y(), wireGeo.Direction().Z());

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

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

      // Make sure arclen is in range
      if (abs(arcLen) < wireGeo.HalfL()) {
        Eigen::Vector3f docaVec = hitPosition - (wirePos + arcLen * wireDir);

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

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

    return distance;
  }

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

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

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

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

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

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

    cet::cpu_timer theClockMakeHits;

    if (m_enableMonitoring) theClockMakeHits.start();

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

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

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

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

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

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

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

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

      std::cout << debugMessage << std::endl;

      m_weHaveAllBeenHereBefore = true;
    }

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

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

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

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

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

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

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

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

        m_planeToSnippetHitMap[planeID][hitStartEndPair].emplace_back(
          &m_clusterHit2DMasterList.back());
        m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back());
      }
    }

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

      wireAssns.addSingle(chanWireItr->second, hitPtrPair.second);
    }
  }

  void SnippetHit3DBuilder::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()) { continue; }

      rawDigitAssns.addSingle(chanRawDigitItr->second, hitPtrPair.second);
    }
  }

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

  DEFINE_ART_CLASS_TOOL(SnippetHit3DBuilder)
} // namespace lar_cluster3d