From ad1d73f901655e0555b09dbae17e7c7ec3b1f3bb Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Fri, 12 Dec 2025 07:41:51 -0600 Subject: [PATCH 01/10] Add Supera as a submodule --- .gitmodules | 4 ++++ dunereco/CMakeLists.txt | 11 +++++++++++ dunereco/Supera | 1 + ups/product_deps | 1 + 4 files changed, 17 insertions(+) create mode 100644 .gitmodules create mode 160000 dunereco/Supera diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..5e846003 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,4 @@ +[submodule "dunereco/Supera"] + path = dunereco/Supera + url = https://github.com/nvpm-lab/Supera.git + branch = develop diff --git a/dunereco/CMakeLists.txt b/dunereco/CMakeLists.txt index f11ff943..70f234de 100644 --- a/dunereco/CMakeLists.txt +++ b/dunereco/CMakeLists.txt @@ -16,6 +16,17 @@ add_subdirectory(TransformerCVN) add_subdirectory(VLNets) add_subdirectory(PointResTree) add_subdirectory(WireFilter) + +if ( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Supera/.git" ) + execute_process(COMMAND git submodule update --init --recursive + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) +endif() + +# We need to run Supera/setup.sh first, otherwise it won't have a CMakeLists.txt. +execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh pdvd + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Supera ) +add_subdirectory(Supera) + # these two were not built in dunetpc at the time of the split # add_subdirectory(SNSlicer) # add_subdirectory(SNUtils) diff --git a/dunereco/Supera b/dunereco/Supera new file mode 160000 index 00000000..bd14cd18 --- /dev/null +++ b/dunereco/Supera @@ -0,0 +1 @@ +Subproject commit bd14cd18aa98d56a275822ae07b18a79c6469591 diff --git a/ups/product_deps b/ups/product_deps index cc58949d..ba6f4816 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -181,6 +181,7 @@ wpdir product_dir wire-cell-cfg product version cetmodules v3_24_01 - only_for_build dunecore v10_14_01d00 +larcv2 v2_3_2 v2_3_2 end_product_list From 6474c05ae36f99df40b138fca74586738be3da62 Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Fri, 12 Dec 2025 16:55:15 -0600 Subject: [PATCH 02/10] fix product_deps --- ups/product_deps | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/ups/product_deps b/ups/product_deps index ba6f4816..5951b16c 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -181,7 +181,7 @@ wpdir product_dir wire-cell-cfg product version cetmodules v3_24_01 - only_for_build dunecore v10_14_01d00 -larcv2 v2_3_2 v2_3_2 +larcv2 v2_2_6 v2_2_6 end_product_list @@ -236,11 +236,11 @@ end_product_list # case it is optional. # #################################### -qualifier dunecore notes -c14:debug c14:debug -c14:prof c14:prof -e26:debug e26:debug -e26:prof e26:prof +qualifier dunecore larcv2 notes +c14:debug c14:debug c14:debug:p3915 +c14:prof c14:prof c14:p3915:prof +e26:debug e26:debug debug:e26:p3915 +e26:prof e26:prof e26:p3915:prof end_qualifier_list #################################### From bc676675d66015b70bcd3e77df3aca88e0461fcb Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Mon, 15 Dec 2025 05:53:32 -0600 Subject: [PATCH 03/10] change pdvd to dune --- dunereco/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dunereco/CMakeLists.txt b/dunereco/CMakeLists.txt index 70f234de..f2cdf323 100644 --- a/dunereco/CMakeLists.txt +++ b/dunereco/CMakeLists.txt @@ -23,7 +23,7 @@ if ( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Supera/.git" ) endif() # We need to run Supera/setup.sh first, otherwise it won't have a CMakeLists.txt. -execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh pdvd +execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh dune WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Supera ) add_subdirectory(Supera) From 08ea462907b5f5dd64d784a58d67d8a9b85b8b4b Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Mon, 15 Dec 2025 06:27:22 -0600 Subject: [PATCH 04/10] Update Supera submodule to nvpm-lab/Supera develop --- dunereco/Supera | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dunereco/Supera b/dunereco/Supera index bd14cd18..9e18035a 160000 --- a/dunereco/Supera +++ b/dunereco/Supera @@ -1 +1 @@ -Subproject commit bd14cd18aa98d56a275822ae07b18a79c6469591 +Subproject commit 9e18035a780a4d597b918bf459f8c6ca92611670 From 29d1415d9c2ee9f3a02d0f8cbe46b73933a9d50b Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Mon, 15 Dec 2025 06:47:46 -0600 Subject: [PATCH 05/10] Update Supera submodule --- dunereco/CMakeLists.txt | 2 +- dunereco/Supera | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/dunereco/CMakeLists.txt b/dunereco/CMakeLists.txt index f2cdf323..70f234de 100644 --- a/dunereco/CMakeLists.txt +++ b/dunereco/CMakeLists.txt @@ -23,7 +23,7 @@ if ( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Supera/.git" ) endif() # We need to run Supera/setup.sh first, otherwise it won't have a CMakeLists.txt. -execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh dune +execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh pdvd WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Supera ) add_subdirectory(Supera) diff --git a/dunereco/Supera b/dunereco/Supera index 9e18035a..e02249cb 160000 --- a/dunereco/Supera +++ b/dunereco/Supera @@ -1 +1 @@ -Subproject commit 9e18035a780a4d597b918bf459f8c6ca92611670 +Subproject commit e02249cbb97e4c4a97c758b9d165ff1c3025c8c5 From 367675577d7b52fd9855c495bfdc8028f5be9231 Mon Sep 17 00:00:00 2001 From: Nived Puthumana Meleppattu Date: Wed, 17 Dec 2025 05:57:20 -0600 Subject: [PATCH 06/10] update submodule --- .gitmodules | 4 ++-- dunereco/CMakeLists.txt | 2 +- dunereco/Supera | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/.gitmodules b/.gitmodules index 5e846003..3b56343f 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,4 +1,4 @@ [submodule "dunereco/Supera"] path = dunereco/Supera - url = https://github.com/nvpm-lab/Supera.git - branch = develop + url = https://github.com/DeepLearnPhysics/Supera.git + branch = dune diff --git a/dunereco/CMakeLists.txt b/dunereco/CMakeLists.txt index 70f234de..f2cdf323 100644 --- a/dunereco/CMakeLists.txt +++ b/dunereco/CMakeLists.txt @@ -23,7 +23,7 @@ if ( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Supera/.git" ) endif() # We need to run Supera/setup.sh first, otherwise it won't have a CMakeLists.txt. -execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh pdvd +execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh dune WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Supera ) add_subdirectory(Supera) diff --git a/dunereco/Supera b/dunereco/Supera index e02249cb..905f1e73 160000 --- a/dunereco/Supera +++ b/dunereco/Supera @@ -1 +1 @@ -Subproject commit e02249cbb97e4c4a97c758b9d165ff1c3025c8c5 +Subproject commit 905f1e7304badec465a66c6a39f2802dabf21e30 From b39c9c07197fb9504ef15cc1eeffaffc663e8d46 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Sun, 21 Dec 2025 13:00:01 -0800 Subject: [PATCH 07/10] Temporary addition of the Cluster3D tool to build space points. This tool has been vetted and is actively used by both SBN experiments. Ultimately it belongs in larreco but that PR appears to be stalled and so adding here temporarily to facilite MC production for upcoming SPINE workshop. --- dunereco/Cluster3D/CMakeLists.txt | 64 + .../Cluster3D/SnippetHit3DBuilderDUNE_tool.cc | 2097 +++++++++++++++++ .../Cluster3D/cluster3dalgorithms_dune.fcl | 28 + 3 files changed, 2189 insertions(+) create mode 100644 dunereco/Cluster3D/CMakeLists.txt create mode 100644 dunereco/Cluster3D/SnippetHit3DBuilderDUNE_tool.cc create mode 100644 dunereco/Cluster3D/cluster3dalgorithms_dune.fcl diff --git a/dunereco/Cluster3D/CMakeLists.txt b/dunereco/Cluster3D/CMakeLists.txt new file mode 100644 index 00000000..dc093414 --- /dev/null +++ b/dunereco/Cluster3D/CMakeLists.txt @@ -0,0 +1,64 @@ +# Following the prescription to build from sbncode +# where should the scripts/..xml file be installed? Perhaps in bin? + +art_make_library( + LIBRARIES + larevt::Filters + lardataalg::DetectorInfo + lardataobj::RecoBase + larcorealg::Geometry + lardata::ArtDataHelper + lardata::RecoObjects + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support + art::Utilities + ROOT::Core + ROOT::Gpad + ROOT::Hist + ROOT::Matrix + ROOT::Physics + ROOT::Tree + art_root_io::tfile_support + art_root_io::TFileService_service + art::Persistency_Provenance + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib + cetlib_except::cetlib_except) +set( TOOL_LIBRARIES + larevt::Filters + larreco::RecoAlg_Cluster3DAlgs + lardataalg::DetectorInfo + lardataobj::RecoBase + lardata::ArtDataHelper + art::Framework_Services_Registry + art_root_io::tfile_support + art_root_io::TFileService_service + ROOT::Tree) +set( MODULE_LIBRARIES + larreco::Calorimetry + larreco::RecoAlg_Cluster3DAlgs + larreco::RecoAlg_TCAlg + larreco::RecoAlg + larreco::ClusterFinder + larsim::MCCheater_ParticleInventoryService_service + lardataobj::AnalysisBase + lardataobj::RecoBase + nurandom::RandomUtils_NuRandomService_service + nusimdata::SimulationBase + art::Framework_Services_Registry + art_root_io::tfile_support + ROOT::Core + ROOT::Physics + art_root_io::TFileService_service + messagefacility::MF_MessageLogger + ) +cet_build_plugin(SnippetHit3DBuilderDUNE art::tool LIBRARIES ${TOOL_LIBRARIES} SOURCE SnippetHit3DBuilderDUNE_tool.cc) +set_property(SOURCE SnippetHit3DBuilderDUNE_tool.cc APPEND PROPERTY COMPILE_DEFINITIONS EIGEN_FFTW_DEFAULT) + +install_headers() +install_fhicl() +install_source() \ No newline at end of file diff --git a/dunereco/Cluster3D/SnippetHit3DBuilderDUNE_tool.cc b/dunereco/Cluster3D/SnippetHit3DBuilderDUNE_tool.cc new file mode 100644 index 00000000..d2880ae2 --- /dev/null +++ b/dunereco/Cluster3D/SnippetHit3DBuilderDUNE_tool.cc @@ -0,0 +1,2097 @@ +/** + * @file SnippetHit3DBuilderDUNE_tool.cc + * + * @brief This tool provides "standard" 3D hits built (by this tool) from 2D hits + * + */ + +// 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 + +// std includes +#include +#include +#include // std::accumulate +#include + +// Ack! +#include "TH1F.h" +#include "TTree.h" + +//------------------------------------------------------------------------------------------------------------------------------------------ +// implementation follows + +namespace lar_cluster3d { + + /** + * @brief What follows are several highly useful typedefs which we + * want to expose to the outside world + */ + + // 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; + using HitStartEndPair = std::pair; + using SnippetHitMap = std::map; + using PlaneToSnippetHitMap = std::map; + using TPCToPlaneToSnippetHitMap = std::map; + using Hit2DList = std::list; + using Hit2DSet = std::set; + using WireToHitSetMap = std::map; + using PlaneToWireToHitSetMap = std::map; + using TPCToPlaneToWireToHitSetMap = std::map; + using HitVectorMap = std::map; + using SnippetHitMapItrPair = std::pair; + using PlaneSnippetHitMapItrPairVec = std::vector; + + /** + * @brief SnippetHit3DBuilderDUNE class definiton + */ + class SnippetHit3DBuilderDUNE : public IHit3DBuilder { + public: + /** + * @brief Constructor + * + * @param pset + */ + explicit SnippetHit3DBuilderDUNE(fhicl::ParameterSet const& pset); + + /** + * @brief Each algorithm may have different objects it wants "produced" so use this to + * let the top level producer module "know" what it is outputting + */ + void produces(art::ProducesCollector&) override; + + void configure(const fhicl::ParameterSet&); + + /** + * @brief Given a set of recob hits, run DBscan to form 3D clusters + * + * @param hitPairList The input list of 3D hits to run clustering on + * @param clusterParametersList A list of cluster objects (parameters from associated hits) + */ + void Hit3DBuilder(art::Event&, reco::HitPairList&, RecobHitToPtrMap&) override; + + /** + * @brief If monitoring, recover the time to execute a particular function + */ + float getTimeToExecute(IHit3DBuilder::TimeValues index) const override + { + return m_timeVector[index]; + } + + private: + /** + * @brief Extract the ART hits and the ART hit-particle relationships + * + * @param evt - the ART event + */ + void CollectArtHits(const art::Event& evt) const; + + /** + * @brief Given the ClusterHit2D objects, build the HitPairMap + */ + void BuildHit3D(reco::HitPairList& hitPairList) const; + + /** + * @brief Create a new 2D hit collection from hits associated to 3D space points + */ + void CreateNewRecobHitCollection(art::Event&, + reco::HitPairList&, + std::vector&, + RecobHitToPtrMap&); + + /** + * @brief Create recob::Wire to recob::Hit associations + */ + void makeWireAssns(const art::Event&, + art::Assns&, + RecobHitToPtrMap&) const; + + /** + * @brief Create raw::RawDigit to recob::Hit associations + */ + void makeRawDigitAssns(const art::Event&, + art::Assns&, + RecobHitToPtrMap&) const; + + /** + * @brief Given the ClusterHit2D objects, build the HitPairMap + */ + size_t BuildHitPairMap(PlaneToSnippetHitMap& planeToHitVectorMap, + reco::HitPairList& hitPairList) const; + + /** + * @brief This defines a structure allowing us to see which hits have already been matched. Generally, if two + * hits have been matched then there is no point considering them anymore as it simply means duplicate + * space points will be created + */ + using MatchedHitMap = + std::unordered_map>; + + /** + * @brief Given the ClusterHit2D objects, build the HitPairMap + */ + size_t BuildHitPairMapByTPC(PlaneSnippetHitMapItrPairVec& planeSnippetHitMapItrPairVec, + reco::HitPairList& hitPairList) const; + + /** + * @brief This builds a list of candidate hit pairs from lists of hits on two planes + */ + using HitMatchTriplet = + std::tuple; + using HitMatchTripletVec = std::vector; + using HitMatchTripletVecMap = std::map; + + int findGoodHitPairs(SnippetHitMap::iterator&, + SnippetHitMap::iterator&, + SnippetHitMap::iterator&, + HitMatchTripletVecMap&) const; + + /** + * @brief This algorithm takes lists of hit pairs and finds good triplets + */ + void findGoodTriplets(HitMatchTripletVecMap&, + HitMatchTripletVecMap&, + MatchedHitMap&, + reco::HitPairList&, + bool = false) const; + + /** + * @brief This will look at storing pair "orphans" where the 2D hits are otherwise unused + */ + + int saveOrphanPairs(HitMatchTripletVecMap&, MatchedHitMap&, reco::HitPairList&) const; + + /** + * @brief This will look at storing pairs where the third channel is consistent with a "bad" channel + */ + + int saveBadChannelPairs(HitMatchTripletVecMap&, MatchedHitMap&, reco::HitPairList&) const; + + /** + * @brief Make a HitPair object by checking two hits + */ + bool makeHitPair(reco::ClusterHit3D& pairOut, + const reco::ClusterHit2D* hit1, + const reco::ClusterHit2D* hit2, + float hitWidthSclFctr = 1., + size_t hitPairCntr = 0) const; + + /** + * @brief Make a 3D HitPair object by checking two hits + */ + bool makeHitTriplet(reco::ClusterHit3D& pairOut, + MatchedHitMap& matchedHitMap, + const reco::ClusterHit3D& pairIn, + const reco::ClusterHit2D* hit2) const; + + /** + * @brief function to detemine if two wires "intersect" (in the 2D sense) + */ + + bool WireIDsIntersect(const geo::WireID&, const geo::WireID&, geo::WireIDIntersection&) const; + + /** + * @brief function to compute the distance of closest approach and the arc length to the points of closest approach + */ + float closestApproach(const Eigen::Vector3f&, + const Eigen::Vector3f&, + const Eigen::Vector3f&, + const Eigen::Vector3f&, + float&, + float&) const; + + /** + * @brief A utility routine for finding a 2D hit closest in time to the given pair + */ + const reco::ClusterHit2D* FindBestMatchingHit(const Hit2DSet& hit2DSet, + const reco::ClusterHit3D& pair, + float pairDeltaTimeLimits) const; + + /** + * @brief A utility routine for returning the number of 2D hits from the list in a given range + */ + int FindNumberInRange(const Hit2DSet& hit2DSet, + const reco::ClusterHit3D& pair, + float range) const; + + /** + * @brief Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range + */ + geo::WireID NearestWireID(const Eigen::Vector3f& position, const geo::WireID& wireID) const; + + /** + * @brief Jacket the calls to finding the nearest wire in order to intercept the exceptions if out of range + */ + float DistanceFromPointToHitWire(const Eigen::Vector3f& position, + const geo::WireID& wireID) const; + + /** + * @brief Perform charge integration between limits + */ + float chargeIntegral(float, float, float, float, int, int) const; + + /** + * @brief clear the tuple vectors before processing next event + */ + void clear(); + + /** + * @brief Data members to follow + */ + using TickCorrectionArray = std::vector>>; + + std::vector m_hitFinderTagVec; + float m_hitWidthSclFctr; + float m_deltaPeakTimeSig; + float m_rangeNumSig; + float m_LongHitStretchFctr; + float m_pulseHeightFrac; + float m_PHLowSelection; + float m_minPHFor2HitPoints; ///< Set a minimum pulse height for 2 hit space point candidates + std::vector m_invalidTPCVec; + float + m_wirePitchScaleFactor; ///< Scaling factor to determine max distance allowed between candidate pairs + float m_maxHit3DChiSquare; ///< Provide ability to select hits based on "chi square" + bool m_saveBadChannelPairs; ///< Should we save pairs consistent with bad channels? + bool m_saveMythicalPoints; ///< Should we save valid 2 hit space points? + float m_maxMythicalChiSquare; ///< Selection cut on mythical points + bool m_useT0Offsets; ///< If true then we will use the LArSoft interplane offsets + bool m_outputHistograms; ///< Take the time to create and fill some histograms for diagnostics + bool m_makeAssociations; ///< Do we make wire/rawdigit associations to space points? + + bool m_enableMonitoring; ///< + float m_wirePitch[3]; + mutable std::vector m_timeVector; ///< + + float m_zPosOffset; + + using PlaneToT0OffsetMap = std::map; + + PlaneToT0OffsetMap m_PlaneToT0OffsetMap; + + // Define some basic histograms + TTree* m_tupleTree; ///< output analysis tree + + mutable std::vector m_deltaPeakTimePlane0Vec; + mutable std::vector m_deltaPeakSigmaPlane0Vec; + mutable std::vector m_deltaPeakTimePlane1Vec; + mutable std::vector m_deltaPeakSigmaPlane1Vec; + mutable std::vector m_deltaPeakTimePlane2Vec; + mutable std::vector m_deltaPeakSigmaPlane2Vec; + + mutable std::vector m_deltaTimeVec; + mutable std::vector m_deltaTime0Vec; + mutable std::vector m_deltaSigma0Vec; + mutable std::vector m_deltaTime1Vec; + mutable std::vector m_deltaSigma1Vec; + mutable std::vector m_deltaTime2Vec; + mutable std::vector m_deltaSigma2Vec; + mutable std::vector m_chiSquare3DVec; + mutable std::vector m_maxPullVec; + mutable std::vector m_overlapFractionVec; + mutable std::vector m_overlapRangeVec; + mutable std::vector m_maxDeltaPeakVec; + mutable std::vector m_maxSideVecVec; + mutable std::vector m_pairWireDistVec; + mutable std::vector m_smallChargeDiffVec; + mutable std::vector m_smallIndexVec; + mutable std::vector m_qualityMetricVec; + mutable std::vector m_spacePointChargeVec; + mutable std::vector m_hitAsymmetryVec; + mutable std::vector m_2hit1stPHVec; + mutable std::vector m_2hit2ndPHVec; + mutable std::vector m_2hitDeltaPHVec; + mutable std::vector m_2hitSumPHVec; + + // Get instances of the primary data structures needed + mutable Hit2DList m_clusterHit2DMasterList; + mutable PlaneToSnippetHitMap m_planeToSnippetHitMap; + mutable PlaneToWireToHitSetMap m_planeToWireToHitSetMap; + + mutable bool m_weHaveAllBeenHereBefore = false; + + const geo::Geometry* m_geometry; //< pointer to the Geometry service + const geo::WireReadoutGeom* m_wireReadoutAlg; //< The new wire readout service + const lariov::ChannelStatusProvider* m_channelFilter; + }; + + SnippetHit3DBuilderDUNE::SnippetHit3DBuilderDUNE(fhicl::ParameterSet const& pset) + : m_channelFilter(&art::ServiceHandle()->GetProvider()) + + { + this->configure(pset); + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + + void SnippetHit3DBuilderDUNE::produces(art::ProducesCollector& collector) + { + collector.produces>(); + collector.produces>(); + collector.produces>(); + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + + void SnippetHit3DBuilderDUNE::configure(fhicl::ParameterSet const& pset) + { + m_hitFinderTagVec = pset.get>("HitFinderTagVec", {"gaushit"}); + m_enableMonitoring = pset.get("EnableMonitoring", true); + m_hitWidthSclFctr = pset.get("HitWidthScaleFactor", 6.); + m_rangeNumSig = pset.get("RangeNumSigma", 3.); + m_LongHitStretchFctr = pset.get("LongHitsStretchFactor", 1.5); + m_pulseHeightFrac = pset.get("PulseHeightFraction", 0.5); + m_PHLowSelection = pset.get("PHLowSelection", 20.); + m_minPHFor2HitPoints = pset.get("MinPHFor2HitPoints", 15.); + m_deltaPeakTimeSig = pset.get("DeltaPeakTimeSig", 1.7); + m_zPosOffset = pset.get("ZPosOffset", 0.0); + m_invalidTPCVec = pset.get>("InvalidTPCVec", std::vector()); + m_wirePitchScaleFactor = pset.get("WirePitchScaleFactor", 1.9); + m_maxHit3DChiSquare = pset.get("MaxHitChiSquare", 6.0); + m_saveBadChannelPairs = pset.get("SaveBadChannelPairs", true); + m_saveMythicalPoints = pset.get("SaveMythicalPoints", true); + m_maxMythicalChiSquare = pset.get("MaxMythicalChiSquare", 10.); + m_useT0Offsets = pset.get("UseT0Offsets", true); + m_outputHistograms = pset.get("OutputHistograms", false); + m_makeAssociations = pset.get("MakeAssociations", false); + + m_geometry = art::ServiceHandle{}.get(); + m_wireReadoutAlg = &art::ServiceHandle + { + } -> 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_wireReadoutAlg->Plane(geo::PlaneID{tpcid, 0}).WirePitch(); + m_wirePitch[1] = m_wireReadoutAlg->Plane(geo::PlaneID{tpcid, 1}).WirePitch(); + m_wirePitch[2] = m_wireReadoutAlg->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 tfs; + + m_tupleTree = tfs->make("Hit3DBuilderTree", "Tree by SnippetHit3DBuilderDUNE"); + + clear(); + + m_tupleTree->Branch("DeltaPeakTimePair0", "std::vector", &m_deltaPeakTimePlane0Vec); + m_tupleTree->Branch("DeltaPeakSigmaPair0", "std::vector", &m_deltaPeakSigmaPlane0Vec); + m_tupleTree->Branch("DeltaPeakTimePair1", "std::vector", &m_deltaPeakTimePlane1Vec); + m_tupleTree->Branch("DeltaPeakSigmaPair1", "std::vector", &m_deltaPeakSigmaPlane1Vec); + m_tupleTree->Branch("DeltaPeakTimePair2", "std::vector", &m_deltaPeakTimePlane2Vec); + m_tupleTree->Branch("DeltaPeakSigmaPair2", "std::vector", &m_deltaPeakSigmaPlane2Vec); + + m_tupleTree->Branch("DeltaTime2D", "std::vector", &m_deltaTimeVec); + m_tupleTree->Branch("DeltaTime2D0", "std::vector", &m_deltaTime0Vec); + m_tupleTree->Branch("DeltaSigma2D0", "std::vector", &m_deltaSigma0Vec); + m_tupleTree->Branch("DeltaTime2D1", "std::vector", &m_deltaTime1Vec); + m_tupleTree->Branch("DeltaSigma2D1", "std::vector", &m_deltaSigma1Vec); + m_tupleTree->Branch("DeltaTime2D2", "std::vector", &m_deltaTime2Vec); + m_tupleTree->Branch("DeltaSigma2D2", "std::vector", &m_deltaSigma2Vec); + m_tupleTree->Branch("ChiSquare3D", "std::vector", &m_chiSquare3DVec); + m_tupleTree->Branch("MaxPullValue", "std::vector", &m_maxPullVec); + m_tupleTree->Branch("OverlapFraction", "std::vector", &m_overlapFractionVec); + m_tupleTree->Branch("OverlapRange", "std::vector", &m_overlapRangeVec); + m_tupleTree->Branch("MaxDeltaPeak", "std::vector", &m_maxDeltaPeakVec); + m_tupleTree->Branch("MaxSideVec", "std::vector", &m_maxSideVecVec); + m_tupleTree->Branch("PairWireDistVec", "std::vector", &m_pairWireDistVec); + m_tupleTree->Branch("SmallChargeDiff", "std::vector", &m_smallChargeDiffVec); + m_tupleTree->Branch("SmallChargeIdx", "std::vector", &m_smallIndexVec); + m_tupleTree->Branch("QualityMetric", "std::vector", &m_qualityMetricVec); + m_tupleTree->Branch("SPCharge", "std::vector", &m_spacePointChargeVec); + m_tupleTree->Branch("HitAsymmetry", "std::vector", &m_hitAsymmetryVec); + + m_tupleTree->Branch("2hit1stPH", "std::vector", &m_2hit1stPHVec); + m_tupleTree->Branch("2hit2ndPH", "std::vector", &m_2hit2ndPHVec); + m_tupleTree->Branch("2hitDeltaPH", "std::vector", &m_2hitDeltaPHVec); + m_tupleTree->Branch("2hitSumPH", "std::vector", &m_2hitSumPHVec); + } + + return; + } + + void SnippetHit3DBuilderDUNE::clear() + { + m_deltaPeakTimePlane0Vec.clear(); + m_deltaPeakSigmaPlane0Vec.clear(); + m_deltaPeakTimePlane1Vec.clear(); + m_deltaPeakSigmaPlane1Vec.clear(); + m_deltaPeakTimePlane1Vec.clear(); + m_deltaPeakSigmaPlane1Vec.clear(); + + m_deltaTimeVec.clear(); + m_deltaTime0Vec.clear(); + m_deltaSigma0Vec.clear(); + m_deltaTime1Vec.clear(); + m_deltaSigma1Vec.clear(); + m_deltaTime2Vec.clear(); + m_deltaSigma2Vec.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(); + + m_2hit1stPHVec.clear(); + m_2hit2ndPHVec.clear(); + m_2hitDeltaPHVec.clear(); + m_2hitSumPHVec.clear(); + + return; + } + + 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 SnippetHit3DBuilderDUNE::Hit3DBuilder(art::Event& evt, + reco::HitPairList& hitPairList, + RecobHitToPtrMap& clusterHitToArtPtrMap) + { + // Clear the internal data structures + m_clusterHit2DMasterList.clear(); + m_planeToSnippetHitMap.clear(); + m_planeToWireToHitSetMap.clear(); + + mf::LogDebug("SnippetHit3D") << "SnippetHit3DBuilderDUNE entered" << std::endl; + + // Do the one time initialization of the tick offsets. + if (m_PlaneToT0OffsetMap.empty()) { + // Need the detector properties which needs the clocks + auto const clock_data = + art::ServiceHandle()->DataFor(evt); + auto const det_prop = + art::ServiceHandle()->DataFor(evt, clock_data); + + // 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++) { + for (size_t planeIdx = 0; planeIdx < m_wireReadoutAlg->Nplanes(); planeIdx++) { + geo::PlaneID planeID(cryoIdx, tpcIdx, planeIdx); + + if (m_useT0Offsets) + m_PlaneToT0OffsetMap[planeID] = + det_prop.GetXTicksOffset(planeID) - + det_prop.GetXTicksOffset(geo::PlaneID(cryoIdx, tpcIdx, 0)); + else + m_PlaneToT0OffsetMap[planeID] = 0.; + } + } + } + } + + m_timeVector.resize(NUMTIMEVALUES, 0.); + + // Get a hit refiner + std::unique_ptr> outputHitPtrVec(new std::vector); + + // Recover the 2D hits and then organize them into data structures which will be used in the + // DBscan algorithm for building the 3D clusters + this->CollectArtHits(evt); + + mf::LogDebug("SnippetHit3D") << "SnippetHit3DBuilderDUNE building space points"; + + // If there are no hits in our view/wire data structure then do not proceed with the full analysis + if (!m_planeToWireToHitSetMap.empty()) { + // Call the algorithm that builds 3D hits + this->BuildHit3D(hitPairList); + + // If we built 3D points then attempt to output a new hit list as well + if (!hitPairList.empty()) + CreateNewRecobHitCollection(evt, hitPairList, *outputHitPtrVec, clusterHitToArtPtrMap); + } + + // Set up to make the associations (if desired) + /// Associations with wires. + std::unique_ptr> wireAssns( + new art::Assns); + + if (m_makeAssociations) makeWireAssns(evt, *wireAssns, clusterHitToArtPtrMap); + + /// Associations with raw digits. + std::unique_ptr> rawDigitAssns( + new art::Assns); + + if (m_makeAssociations) 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(); + } + + mf::LogDebug("SnippetHit3D") << "SnippetHit3DBuilderDUNE exited"; + + return; + } + + void SnippetHit3DBuilderDUNE::BuildHit3D(reco::HitPairList& hitPairList) const + { + /** + * @brief Driver for processing input 2D hits, transforming to 3D hits and building lists + * of associated 3D hits (candidate 3D clusters) + */ + cet::cpu_timer theClockMakeHits; + + if (m_enableMonitoring) theClockMakeHits.start(); + + size_t numHitPairs = BuildHitPairMap(m_planeToSnippetHitMap, hitPairList); + + if (m_enableMonitoring) { + theClockMakeHits.stop(); + + m_timeVector[BUILDTHREEDHITS] = theClockMakeHits.accumulated_real_time(); + } + + mf::LogDebug("SnippetHit3D") << ">>>>> 3D hit building done, found " << numHitPairs + << " 3D Hits" << std::endl; + + return; + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + class SetStartTimeOrder { + public: + 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; + } + + private: + }; + + 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 SnippetHit3DBuilderDUNE::BuildHitPairMap(PlaneToSnippetHitMap& planeToSnippetHitMap, + reco::HitPairList& hitPairList) const + { + /** + * @brief Given input 2D hits, build out the lists of possible 3D hits + * + * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view + * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be + * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes + * and then look for the match in the V plane. In the event we don't find the match in the V plane then we + * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high. + */ + size_t totalNumHits(0); + size_t hitPairCntr(0); + + size_t nTriplets(0); + + // Set up to loop over cryostats and tpcs... + for (size_t cryoIdx = 0; cryoIdx < m_geometry->Ncryostats(); cryoIdx++) { + for (size_t tpcIdx = 0; tpcIdx < m_geometry->NTPC(); tpcIdx++) { + //************************************ + // Kludge + // if (!(cryoIdx == 1 && tpcIdx == 0)) continue; + + PlaneToSnippetHitMap::iterator mapItr0 = + planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0)); + PlaneToSnippetHitMap::iterator mapItr1 = + planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1)); + PlaneToSnippetHitMap::iterator mapItr2 = + planeToSnippetHitMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2)); + + size_t nPlanesWithHits = + (mapItr0 != planeToSnippetHitMap.end() && !mapItr0->second.empty() ? 1 : 0) + + (mapItr1 != planeToSnippetHitMap.end() && !mapItr1->second.empty() ? 1 : 0) + + (mapItr2 != planeToSnippetHitMap.end() && !mapItr2->second.empty() ? 1 : 0); + + if (nPlanesWithHits < 2) continue; + + SnippetHitMap& snippetHitMap0 = mapItr0->second; + SnippetHitMap& snippetHitMap1 = mapItr1->second; + SnippetHitMap& snippetHitMap2 = mapItr2->second; + + PlaneSnippetHitMapItrPairVec hitItrVec = { + SnippetHitMapItrPair(snippetHitMap0.begin(), snippetHitMap0.end()), + SnippetHitMapItrPair(snippetHitMap1.begin(), snippetHitMap1.end()), + SnippetHitMapItrPair(snippetHitMap2.begin(), snippetHitMap2.end())}; + + totalNumHits += BuildHitPairMapByTPC(hitItrVec, hitPairList); + } + } + + // Return the hit pair list but sorted by z and y positions (faster traversal in next steps) + hitPairList.sort(SetPairStartTimeOrder); + + // Where are we? + mf::LogDebug("SnippetHit3D") << "Total number hits: " << totalNumHits << std::endl; + mf::LogDebug("SnippetHit3D") << "Created a total of " << hitPairList.size() + << " hit pairs, counted: " << hitPairCntr << std::endl; + mf::LogDebug("SnippetHit3D") << "-- Triplets: " << nTriplets << std::endl; + + return hitPairList.size(); + } + + size_t SnippetHit3DBuilderDUNE::BuildHitPairMapByTPC( + PlaneSnippetHitMapItrPairVec& snippetHitMapItrVec, + reco::HitPairList& hitPairList) const + { + /** + * @brief Given input 2D hits, build out the lists of possible 3D hits + * + * The current strategy: ideally all 3D hits would be comprised of a triplet of 2D hits, one from each view + * However, we have concern that, in particular, the v-plane may have some inefficiency which we have to be + * be prepared to deal with. The idea, then, is to first make the association of hits in the U and W planes + * and then look for the match in the V plane. In the event we don't find the match in the V plane then we + * will evaluate the situation and in some instances keep the U-W pairs in order to keep efficiency high. + */ + + // Define functions to set start/end iterators in the loop below + auto SetStartIterator = + [](SnippetHitMap::iterator startItr, SnippetHitMap::iterator endItr, float startTime) { + while (startItr != endItr) { + if (startItr->first.second < startTime) + startItr++; + else + break; + } + return startItr; + }; + + auto SetEndIterator = + [](SnippetHitMap::iterator lastItr, SnippetHitMap::iterator endItr, float endTime) { + while (lastItr != endItr) { + if (lastItr->first.first < endTime) + lastItr++; + else + break; + } + return lastItr; + }; + + size_t nTriplets(0); + size_t nOrphanPairs(0); + + // Structure to keep track of hit associations + MatchedHitMap matchedHitMap; + + //********************************************************************************* + // Basically, we try to loop until done... + while (1) { + // Sort so that the earliest hit time will be the first element, etc. + std::sort(snippetHitMapItrVec.begin(), snippetHitMapItrVec.end(), SetStartTimeOrder()); + + // Make sure there are still hits on at least + int nPlanesWithHits(0); + + for (auto& pair : snippetHitMapItrVec) + if (pair.first != pair.second) nPlanesWithHits++; + + if (nPlanesWithHits < 2) break; + + // End condition: no more hit snippets + // if (snippetHitMapItrVec.front().first == snippetHitMapItrVec.front().second) break; + + // This loop iteration's snippet iterator + SnippetHitMap::iterator firstSnippetItr = snippetHitMapItrVec.front().first; + + // Set iterators to insure we'll be in the overlap ranges + SnippetHitMap::iterator snippetHitMapItr1Start = SetStartIterator( + snippetHitMapItrVec[1].first, snippetHitMapItrVec[1].second, firstSnippetItr->first.first); + SnippetHitMap::iterator snippetHitMapItr1End = SetEndIterator( + snippetHitMapItr1Start, snippetHitMapItrVec[1].second, firstSnippetItr->first.second); + SnippetHitMap::iterator snippetHitMapItr2Start = SetStartIterator( + snippetHitMapItrVec[2].first, snippetHitMapItrVec[2].second, firstSnippetItr->first.first); + SnippetHitMap::iterator snippetHitMapItr2End = SetEndIterator( + snippetHitMapItr2Start, snippetHitMapItrVec[2].second, firstSnippetItr->first.second); + + // Since we'll use these many times in the internal loops, pre make the pairs for the second set of hits + size_t curHitListSize(hitPairList.size()); + HitMatchTripletVecMap pair12Map; + HitMatchTripletVecMap pair13Map; + + size_t n12Pairs = + findGoodHitPairs(firstSnippetItr, snippetHitMapItr1Start, snippetHitMapItr1End, pair12Map); + size_t n13Pairs = + findGoodHitPairs(firstSnippetItr, snippetHitMapItr2Start, snippetHitMapItr2End, pair13Map); + + curHitListSize = hitPairList.size(); + + if (n12Pairs > n13Pairs) + findGoodTriplets(pair12Map, pair13Map, matchedHitMap, hitPairList); + else + findGoodTriplets(pair13Map, pair12Map, matchedHitMap, hitPairList); + + if (m_saveBadChannelPairs) { + nOrphanPairs += saveBadChannelPairs(pair12Map, matchedHitMap, hitPairList); + nOrphanPairs += saveBadChannelPairs(pair13Map, matchedHitMap, hitPairList); + } + + if (m_saveMythicalPoints) { + nOrphanPairs += saveOrphanPairs(pair12Map, matchedHitMap, hitPairList); + nOrphanPairs += saveOrphanPairs(pair13Map, matchedHitMap, hitPairList); + } + + nTriplets += hitPairList.size() - curHitListSize; + + snippetHitMapItrVec.front().first++; + } + + mf::LogDebug("SnippetHit3D") << "--> Created " << nTriplets << " triplets of which " + << nOrphanPairs << " are orphans" << std::endl; + + return hitPairList.size(); + } + + int SnippetHit3DBuilderDUNE::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 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 SnippetHit3DBuilderDUNE::findGoodTriplets(HitMatchTripletVecMap& pair12Map, + HitMatchTripletVecMap& pair13Map, + MatchedHitMap& matchedHitMap, + reco::HitPairList& hitPairList, + bool tagged) const + { + // Build triplets from the two lists of hit pairs + if (!pair12Map.empty() && !pair13Map.empty()) { + // temporary container for dead channel hits + std::vector tempDeadChanVec; + reco::ClusterHit3D deadChanPair; + + // Keep track of which third plane hits have been used + std::map 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, matchedHitMap, pair1, hit2)) { + triplet.setID(hitPairList.size()); + hitPairList.emplace_back(triplet); + usedPairMap[&pair1] = true; + usedPairMap[&pair2] = true; + } + } + } + } + } + } + + return; + } + + int SnippetHit3DBuilderDUNE::saveOrphanPairs(HitMatchTripletVecMap& pairMap, + MatchedHitMap& matchedHitMap, + reco::HitPairList& hitPairList) const + { + int curTripletCount = hitPairList.size(); + + // Build triplets from the two lists of hit pairs + if (!pairMap.empty()) { + // Initial population of this map with the pair13Map hits + for (const auto& pair : pairMap) { + if (pair.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& hit2Dhit3DPair : pair.second) { + const reco::ClusterHit3D& hit3D = std::get<2>(hit2Dhit3DPair); + + // No point considering a 3D hit that has been used to make a space point already + if (hit3D.getStatusBits() & reco::ClusterHit3D::MADESPACEPOINT) continue; + + const reco::ClusterHit2D* hit1 = std::get<0>(hit2Dhit3DPair); + const reco::ClusterHit2D* hit2 = std::get<1>(hit2Dhit3DPair); + + // Have these already been used in some fashion? + if (matchedHitMap[hit1].find(hit2) != matchedHitMap[hit1].end()) continue; + + if (m_outputHistograms) { + m_2hit1stPHVec.emplace_back(hit1->getHit()->PeakAmplitude()); + m_2hit2ndPHVec.emplace_back(hit2->getHit()->PeakAmplitude()); + m_2hitDeltaPHVec.emplace_back(hit2->getHit()->PeakAmplitude() - + hit1->getHit()->PeakAmplitude()); + m_2hitSumPHVec.emplace_back(hit2->getHit()->PeakAmplitude() + + hit1->getHit()->PeakAmplitude()); + } + + // If both hits already appear in a triplet then there is no gain here so reject + if (hit1->getHit()->PeakAmplitude() < m_minPHFor2HitPoints || + hit2->getHit()->PeakAmplitude() < m_minPHFor2HitPoints) + continue; + + // Require that one of the hits is on the collection plane + if (hit1->WireID().Plane == 2 || hit2->WireID().Plane == 2) { + // Allow cut on the quality of the space point + if (hit3D.getHitChiSquare() < m_maxMythicalChiSquare) { + // Add to the list + hitPairList.emplace_back(hit3D); + hitPairList.back().setID(hitPairList.size() - 1); + matchedHitMap[hit1].insert(hit2); + matchedHitMap[hit2].insert(hit1); + } + } + } + } + } + + return hitPairList.size() - curTripletCount; + } + + int SnippetHit3DBuilderDUNE::saveBadChannelPairs(HitMatchTripletVecMap& pairMap, + MatchedHitMap& matchedHitMap, + reco::HitPairList& hitPairList) const + { + int curTripletCount = hitPairList.size(); + + // Assuming we havebad channels and the pairmap is not empty... + if (!m_channelFilter->NoisyChannels().empty() && !pairMap.empty()) { + // Initial population of this map with the pairMap hits + for (const auto& pair : pairMap) { + if (pair.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& hit2Dhit3DPair : pair.second) { + const reco::ClusterHit3D& hit3D = std::get<2>(hit2Dhit3DPair); + + // No point considering a 3D hit that has been used to make a space point already + if (hit3D.getStatusBits() & reco::ClusterHit3D::MADESPACEPOINT) continue; + + const reco::ClusterHit2D* hit0 = std::get<0>(hit2Dhit3DPair); + const reco::ClusterHit2D* hit1 = std::get<1>(hit2Dhit3DPair); + + // Have both of these hits already been used in some fashion? + if (matchedHitMap.find(hit0) != matchedHitMap.end() && + matchedHitMap[hit0].find(hit1) != matchedHitMap[hit0].end()) + continue; + + // Which plane is missing? + geo::WireID wireID0 = hit0->WireID(); + geo::WireID wireID1 = hit1->WireID(); + + // Determine the missing plane + int missPlane = 3 - (wireID0.Plane + wireID1.Plane); + + // Ok, recover the wireID expected in the third plane... + geo::WireID wireIn(wireID0.Cryostat, wireID0.TPC, missPlane, 0); + geo::WireID wireID = NearestWireID(hit3D.getPosition(), wireIn); + + raw::ChannelID_t channel = m_wireReadoutAlg->PlaneWireToChannel(wireID); + + geo::WireID wireIDNext(wireID.Cryostat, wireID.TPC, wireID.Plane, wireID.Wire + 1); + raw::ChannelID_t chanNext = m_wireReadoutAlg->PlaneWireToChannel(wireIDNext); + + bool wireStatus = + m_channelFilter->IsPresent(channel) && !m_channelFilter->IsGood(channel); + bool wireOneStatus = + m_channelFilter->IsPresent(chanNext) && !m_channelFilter->IsGood(chanNext); + + // 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 = wireIDNext; + + // Want to refine position since we "know" the missing wire + if (auto widIntersect0 = m_wireReadoutAlg->WireIDsIntersect(wireID0, wireID)) { + if (auto widIntersect1 = m_wireReadoutAlg->WireIDsIntersect(wireID1, wireID)) { + Eigen::Vector3f newPosition( + hit3D.getPosition()[0], hit3D.getPosition()[1], hit3D.getPosition()[2]); + + newPosition[1] = (newPosition[1] + widIntersect0->y + widIntersect1->y) / 3.; + newPosition[2] = + (newPosition[2] + widIntersect0->z + widIntersect1->z - 2. * m_zPosOffset) / 3.; + + hit3D.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); + + // Add to the list + hitPairList.emplace_back(hit3D); + hitPairList.back().setID(hitPairList.size() - 1); + matchedHitMap[hit0].insert(hit1); + matchedHitMap[hit1].insert(hit0); + } + } + } + } + } + } + + return hitPairList.size() - curTripletCount; + } + + bool SnippetHit3DBuilderDUNE::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 = std::pow(hit1Sigma, 2); + float hit2SigSq = std::pow(hit2Sigma, 2); + float deltaPeakTime = std::fabs(hit1Peak - hit2Peak); + float sigmaPeakTime = std::sqrt(hit1SigSq + hit2SigSq); + + if (m_outputHistograms) { + // brute force... sigh... + int plane1 = hit1->WireID().Plane; + int plane2 = hit2->WireID().Plane; + int planeIdx = (plane1 + plane2) - 1; // should be 0 for 0-1, 1 for 0-2 and 2 for 1-2 + + float deltaTicks = hit2Peak - hit1Peak; + + if (plane1 > plane2) deltaTicks = -deltaTicks; + + if (planeIdx == 0) { + m_deltaPeakTimePlane0Vec.emplace_back(deltaTicks); + m_deltaPeakSigmaPlane0Vec.emplace_back(sigmaPeakTime); + } + else if (planeIdx == 1) { + m_deltaPeakTimePlane1Vec.emplace_back(deltaTicks); + m_deltaPeakSigmaPlane1Vec.emplace_back(sigmaPeakTime); + } + else { + m_deltaPeakTimePlane2Vec.emplace_back(deltaTicks); + m_deltaPeakSigmaPlane2Vec.emplace_back(sigmaPeakTime); + } + } + + // 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 averageCharge = 0.5 * (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 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 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, + averageCharge, + avePeakTime, + deltaPeakTime, + sigmaPeakTime, + hitChiSquare, + 0., + 0., + 0., + 0., + hitVector, + hitDelTSigVec, + wireIDVec); + + result = true; + } + } + } + + // Send it back + return result; + } + + bool SnippetHit3DBuilderDUNE::makeHitTriplet(reco::ClusterHit3D& hitTriplet, + MatchedHitMap& matchedHitMap, + 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); + + // Let's mark that this pair had a valid 3 wire combination + pair.setStatusBit(reco::ClusterHit3D::MADESPACEPOINT); + + // 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]; + + bool notMatched = matchedHitMap[hit].find(hit0) == matchedHitMap[hit].end() && + matchedHitMap[hit].find(hit1) == matchedHitMap[hit].end() && + matchedHitMap[hit0].find(hit1) == matchedHitMap[hit0].end(); + + // If good pairs made here then we can try to make a triplet + if (notMatched && 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 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(); + + float hitRMS = hit2D->getHit()->RMS(); + float peakTime = hit2D->getTimeTicks(); + + // Basically, allow the range to extend to the nearest end of the snippet + if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr; + //hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks()); + + float weight = 1. / (hitRMS * hitRMS); + + avePeakTime += peakTime * weight; + xPosition += hit2D->getXPosition() * weight; + weightSum += weight; + } + + avePeakTime /= weightSum; + xPosition /= weightSum; + + Eigen::Vector2f pair0hYZVec(pair0h.getPosition()[1], pair0h.getPosition()[2]); + Eigen::Vector2f pair1hYZVec(pair1h.getPosition()[1], pair1h.getPosition()[2]); + Eigen::Vector2f pairYZVec(pair.getPosition()[1], pair.getPosition()[2]); + Eigen::Vector3f position(xPosition, + float((pairYZVec[0] + pair0hYZVec[0] + pair1hYZVec[0]) / 3.), + float((pairYZVec[1] + pair0hYZVec[1] + pair1hYZVec[1]) / 3.)); + + // Armed with the average peak time, now get hitChiSquare and the sig vec + float hitChiSquare(0.); + float sigmaPeakTime(std::sqrt(1. / weightSum)); + std::vector hitDelTSigVec; + + for (const auto& hit2D : hitVector) { + float hitRMS = hit2D->getHit()->RMS(); + + // Basically, allow the range to extend to the nearest end of the snippet + if (hit2D->getHit()->DegreesOfFreedom() < 2) hitRMS *= m_LongHitStretchFctr; + //hitRMS = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks()); + + float combRMS = std::sqrt(hitRMS * hitRMS - sigmaPeakTime * sigmaPeakTime); + float peakTime = hit2D->getTimeTicks(); + float deltaTime = peakTime - avePeakTime; + float hitSig = deltaTime / combRMS; + + hitChiSquare += hitSig * hitSig; + + hitDelTSigVec.emplace_back(std::fabs(hitSig)); + } + + if (m_outputHistograms) m_chiSquare3DVec.push_back(hitChiSquare); + + int lowMinIndex(std::numeric_limits::max()); + int lowMaxIndex(std::numeric_limits::min()); + int hiMinIndex(std::numeric_limits::max()); + int hiMaxIndex(std::numeric_limits::min()); + + // First task is to get the min/max values for the common overlap region + for (const auto& hit2D : hitVector) { + float range = m_rangeNumSig * hit2D->getHit()->RMS(); + + // Basically, allow the range to extend to the nearest end of the snippet + if (hit2D->getHit()->DegreesOfFreedom() < 2) range *= m_LongHitStretchFctr; + //range = std::min(hit2D->getTimeTicks() - float(hit2D->getHit()->StartTick()),float(hit2D->getHit()->EndTick())-hit2D->getTimeTicks()); + + int hitStart = hit2D->getHit()->PeakTime() - range - 0.5; + int hitStop = hit2D->getHit()->PeakTime() + range + 0.5; + + lowMinIndex = std::min(hitStart, lowMinIndex); + lowMaxIndex = std::max(hitStart, lowMaxIndex); + hiMinIndex = std::min(hitStop + 1, hiMinIndex); + hiMaxIndex = std::max(hitStop + 1, hiMaxIndex); + } + + // Keep only "good" hits... + if (hitChiSquare < m_maxHit3DChiSquare && hiMinIndex > lowMaxIndex) { + // One more pass through hits to get charge + std::vector chargeVec; + + for (const auto& hit2D : hitVector) + chargeVec.push_back(chargeIntegral(hit2D->getHit()->PeakTime(), + hit2D->getHit()->PeakAmplitude(), + hit2D->getHit()->RMS(), + 1., + lowMaxIndex, + hiMinIndex)); + + float totalCharge = + std::accumulate(chargeVec.begin(), chargeVec.end(), 0.) / float(chargeVec.size()); + float overlapRange = float(hiMinIndex - lowMaxIndex); + float overlapFraction = overlapRange / float(hiMaxIndex - lowMinIndex); + + // Set up to compute the charge asymmetry + std::vector smallestChargeDiffVec; + std::vector chargeAveVec; + float smallestDiff(std::numeric_limits::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[rightIdx]->getTimeTicks() - hitVector[leftIdx]->getTimeTicks(); + + if (std::abs(deltaPeakTime) > maxDeltaPeak) maxDeltaPeak = std::abs(deltaPeakTime); + + if (m_outputHistograms) { + int deltaTimeIdx = + (hitVector[leftIdx]->WireID().Plane + hitVector[rightIdx]->WireID().Plane) - 1; + float combRMS = std::sqrt(std::pow(hitVector[leftIdx]->getHit()->RMS(), 2) + + std::pow(hitVector[rightIdx]->getHit()->RMS(), 2)); + + m_deltaTimeVec.push_back(deltaPeakTime); + + // Want to get the sign of the difference correct... + if (deltaTimeIdx == 0) // This is planes 1 and 0 + { + m_deltaTime0Vec.emplace_back( + float(hitVector[1]->getTimeTicks() - hitVector[0]->getTimeTicks())); + m_deltaSigma0Vec.emplace_back(combRMS); + } + else if (deltaTimeIdx == 1) // This is planes 0 and 2 + { + m_deltaTime1Vec.emplace_back( + float(hitVector[2]->getTimeTicks() - hitVector[0]->getTimeTicks())); + m_deltaSigma1Vec.emplace_back(combRMS); + } + else // must be planes 1 and 2 + { + m_deltaTime2Vec.emplace_back( + float(hitVector[2]->getTimeTicks() - hitVector[1]->getTimeTicks())); + m_deltaSigma2Vec.emplace_back(combRMS); + } + } + } + + 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); + + // Since we are keeping the triplet, mark the hits as used + for (size_t hit2DIdx = 0; hit2DIdx < hitVector.size(); hit2DIdx++) { + const reco::ClusterHit2D* hit2D = hitVector[hit2DIdx]; + + matchedHitMap[hit2D].insert(hitVector[(hit2DIdx + 1) % hitVector.size()]); + matchedHitMap[hit2D].insert(hitVector[(hit2DIdx + 2) % hitVector.size()]); + + if (hit2D->getStatusBits() & reco::ClusterHit2D::USEDINTRIPLET) + hit2D->setStatusBit(reco::ClusterHit2D::SHAREDINTRIPLET); + + hit2D->setStatusBit(reco::ClusterHit2D::USEDINTRIPLET); + } + + // Mark the input pair + pair.setStatusBit(reco::ClusterHit3D::MADESPACEPOINT); + + result = true; + } + // else mf::LogDebug("SnippetHit3D") << "-Rejecting triple with chiSquare: " << hitChiSquare << " and hiMinIndex: " << hiMinIndex << ", loMaxIndex: " << lowMaxIndex; + } + } + } + // else mf::LogDebug("SnippetHit3D") << "-MakeTriplet hit cut, delta: " << hitTimeTicks - pair.getAvePeakTime() << ", min scale fctr: " <Wire(wireID0); + const geo::WireGeo& wireGeo1 = m_wireReadoutAlg->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()); + + //********************************* + // Kludge + // if (wireID0.Plane > 0) wireDir0[2] = -wireDir0[2]; + + // And now the second one + wirePosArr = wireGeo1.GetCenter(); + + Eigen::Vector3f wirePos1(wirePosArr.X(), wirePosArr.Y(), wirePosArr.Z()); + Eigen::Vector3f wireDir1( + wireGeo1.Direction().X(), wireGeo1.Direction().Y(), wireGeo1.Direction().Z()); + + //********************************** + // Kludge + // if (wireID1.Plane > 0) wireDir1[2] = -wireDir1[2]; + + // Get the distance of closest approach + float arcLen0; + float arcLen1; + + if (closestApproach(wirePos0, wireDir0, wirePos1, wireDir1, arcLen0, arcLen1)) { + // Now check that arc lengths are within range + if (std::abs(arcLen0) < wireGeo0.HalfL() && std::abs(arcLen1) < wireGeo1.HalfL()) { + Eigen::Vector3f poca0 = wirePos0 + arcLen0 * wireDir0; + + widIntersection.y = poca0[1]; + widIntersection.z = poca0[2]; + + success = true; + } + } + + return success; + } + + float SnippetHit3DBuilderDUNE::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 SnippetHit3DBuilderDUNE::chargeIntegral(float peakMean, + float peakAmp, + float peakSigma, + float areaNorm, + 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; + } + + const reco::ClusterHit2D* SnippetHit3DBuilderDUNE::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 SnippetHit3DBuilderDUNE::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 SnippetHit3DBuilderDUNE::NearestWireID(const Eigen::Vector3f& position, + const geo::WireID& wireIDIn) const + { + int wire(0); + + // Embed the call to the geometry's services nearest wire id method in a try-catch block + try { + // Not sure the thinking above but is wrong... switch back to NearestWireID... + wire = (m_wireReadoutAlg->NearestWireID(geo::vect::toPoint(position.data()), wireIDIn)).Wire; + } + catch (std::exception& exc) { + // This can happen, almost always because the coordinates are **just** out of range + mf::LogDebug("Cluster3D") << "Exception caught finding nearest wire, position - " + << exc.what() << std::endl; + + // Assume extremum for wire number depending on z coordinate + if (position[2] < m_geometry->TPC({0, 0}).ActiveLength() * 0.5) + wire = 0; + else + wire = m_wireReadoutAlg->Nwires(wireIDIn.asPlaneID()) - 1; + } + + geo::WireID wireID(wireIDIn.Cryostat, wireIDIn.TPC, wireIDIn.Plane, wire); + + return wireID; + } + + float SnippetHit3DBuilderDUNE::DistanceFromPointToHitWire(const Eigen::Vector3f& position, + const geo::WireID& wireIDIn) const + { + float distance = std::numeric_limits::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_wireReadoutAlg->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()); + + //********************************* + // Kludge + // if (wireIDIn.Plane > 0) wireDir[2] = -wireDir[2]; + + // Want the hit position to have same x value as wire coordinates + Eigen::Vector3f hitPosition(wirePos[0], position[1], position[2]); + + // Get arc length to doca + double arcLen = (hitPosition - wirePos).dot(wireDir); + + // Make sure arclen is in range + if (abs(arcLen) < wireGeo.HalfL()) { + Eigen::Vector3f docaVec = hitPosition - (wirePos + arcLen * wireDir); + + distance = docaVec.norm(); + } + } + catch (std::exception& exc) { + // This can happen, almost always because the coordinates are **just** out of range + mf::LogWarning("Cluster3D") << "Exception caught finding nearest wire, position - " + << exc.what() << std::endl; + + // Assume extremum for wire number depending on z coordinate + distance = 0.; + } + + return distance; + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + 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 SnippetHit3DBuilderDUNE::CollectArtHits(const art::Event& evt) const + { + /** + * @brief Recover the 2D hits from art and fill out the local data structures for the 3D clustering + */ + + // Start by getting a vector of valid, non empty hit collections to make sure we really have something to do here... + // Here is a container for the hits... + std::vector recobHitVec; + + // Loop through the list of input sources + for (const auto& inputTag : m_hitFinderTagVec) { + art::Handle> 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(); + + // Need the detector properties which needs the clocks + auto const clock_data = + art::ServiceHandle()->DataFor(evt); + auto const det_prop = + art::ServiceHandle()->DataFor(evt, clock_data); + + // Try to output a formatted string + std::string debugMessage(""); + + // Keep track of x position limits + std::map planeIDToPositionMap; + + // 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(); + + // Should we provide output? + if (!m_weHaveAllBeenHereBefore) { + std::ostringstream outputString; + + outputString << "***> plane 0 offset: " + << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx, tpcIdx, 0))->second + << ", plane 1: " + << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx, tpcIdx, 1))->second + << ", plane 2: " + << m_PlaneToT0OffsetMap.find(geo::PlaneID(cryoIdx, tpcIdx, 2))->second + << "\n"; + 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() + "\n"; + } + + double xPosition(det_prop.ConvertTicksToX(0., 2, tpcIdx, cryoIdx)); + + planeIDToPositionMap[geo::PlaneID(cryoIdx, tpcIdx, 2)] = xPosition; + } + } + + if (!m_weHaveAllBeenHereBefore) { + for (const auto& planeToPositionPair : planeIDToPositionMap) { + std::ostringstream outputString; + + outputString << "***> Plane " << planeToPositionPair.first + << " has time=0 position: " << planeToPositionPair.second << "\n"; + + debugMessage += outputString.str(); + } + + mf::LogDebug("SnippetHit3D") << 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& wireIDs = + m_wireReadoutAlg->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) { + mf::LogDebug("SnippetHit3D") + << "Yes, found a hit with end time less than start time: " << hitStartEndPair.first << "/" + << hitStartEndPair.second << ", mult: " << recobHit->Multiplicity(); + continue; + } + + // And then loop over all possible to build out our maps + //for(const auto& wireID : wireIDs) + for (auto wireID : wireIDs) { + // Check if this is an invalid TPC + // (for example, in protoDUNE there are logical TPC's which see no signal) + if (std::find(m_invalidTPCVec.begin(), m_invalidTPCVec.end(), wireID.TPC) != + m_invalidTPCVec.end()) + continue; + + // Note that a plane ID will define cryostat, TPC and plane + const geo::PlaneID& planeID = wireID.planeID(); + + double hitPeakTime(recobHit->PeakTime() - m_PlaneToT0OffsetMap.find(planeID)->second); + double xPosition(det_prop.ConvertTicksToX( + recobHit->PeakTime(), planeID.Plane, planeID.TPC, planeID.Cryostat)); + + m_clusterHit2DMasterList.emplace_back(0, 0., 0., xPosition, hitPeakTime, wireID, recobHit); + + m_planeToSnippetHitMap[planeID][hitStartEndPair].emplace_back( + &m_clusterHit2DMasterList.back()); + m_planeToWireToHitSetMap[planeID][wireID.Wire].insert(&m_clusterHit2DMasterList.back()); + } + } + + // Make a loop through to sort the recover hits in time order + // for(auto& hitVectorMap : m_planeToSnippetHitMap) + // std::sort(hitVectorMap.second.begin(), hitVectorMap.second.end(), SetHitTimeOrder); + + if (m_enableMonitoring) { + theClockMakeHits.stop(); + + m_timeVector[COLLECTARTHITS] = theClockMakeHits.accumulated_real_time(); + } + + mf::LogDebug("SnippetHit3D") << ">>>>> Number of ART hits: " << m_clusterHit2DMasterList.size() + << std::endl; + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + + void SnippetHit3DBuilderDUNE::CreateNewRecobHitCollection(art::Event& event, + reco::HitPairList& hitPairList, + std::vector& 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 visitedHit2DSet; + + // Use this handy art utility to make art::Ptr objects to the new recob::Hits for use in the output phase + art::PtrMaker 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]; + + if (!hit2D) continue; + + // 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(hit2D)->setHit(newHit); + } + } + } + + size_t numNewHits = hitPtrVec.size(); + + if (m_enableMonitoring) { + theClockBuildNewHits.stop(); + + m_timeVector[BUILDNEWHITS] = theClockBuildNewHits.accumulated_real_time(); + } + + mf::LogDebug("SnippetHit3D") << ">>>>> New output recob::Hit size: " << numNewHits << " (vs " + << m_clusterHit2DMasterList.size() << " input)" << std::endl; + + return; + } + + void SnippetHit3DBuilderDUNE::makeWireAssns(const art::Event& evt, + art::Assns& wireAssns, + RecobHitToPtrMap& recobHitPtrMap) const + { + // Let's make sure the input associations container is empty + wireAssns = art::Assns(); + + // 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> channelToWireMap; + + // Go through the list of input sources and fill out the map + for (const auto& inputTag : m_hitFinderTagVec) { + art::ValidHandle> hitHandle = + evt.getValidHandle>(inputTag); + + art::FindOneP hitToWireAssns(hitHandle, evt, inputTag); + + if (hitToWireAssns.isValid()) { + for (size_t wireIdx = 0; wireIdx < hitToWireAssns.size(); wireIdx++) { + art::Ptr 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>::iterator chanWireItr = + channelToWireMap.find(channel); + + if (!(chanWireItr != channelToWireMap.end())) { + //mf::LogDebug("SnippetHit3D") << "** Did not find channel to wire match! Skipping..." << std::endl; + continue; + } + + wireAssns.addSingle(chanWireItr->second, hitPtrPair.second); + } + + return; + } + + void SnippetHit3DBuilderDUNE::makeRawDigitAssns(const art::Event& evt, + art::Assns& rawDigitAssns, + RecobHitToPtrMap& recobHitPtrMap) const + { + // Let's make sure the input associations container is empty + rawDigitAssns = art::Assns(); + + // 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> channelToRawDigitMap; + + // Go through the list of input sources and fill out the map + for (const auto& inputTag : m_hitFinderTagVec) { + art::ValidHandle> hitHandle = + evt.getValidHandle>(inputTag); + + art::FindOneP hitToRawDigitAssns(hitHandle, evt, inputTag); + + if (hitToRawDigitAssns.isValid()) { + for (size_t rawDigitIdx = 0; rawDigitIdx < hitToRawDigitAssns.size(); rawDigitIdx++) { + art::Ptr 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>::iterator chanRawDigitItr = + channelToRawDigitMap.find(channel); + + if (chanRawDigitItr == channelToRawDigitMap.end()) { + //mf::LogDebug("SnippetHit3D") << "** Did not find channel to wire match! Skipping..." << std::endl; + continue; + } + + rawDigitAssns.addSingle(chanRawDigitItr->second, hitPtrPair.second); + } + + return; + } + + //------------------------------------------------------------------------------------------------------------------------------------------ + + DEFINE_ART_CLASS_TOOL(SnippetHit3DBuilderDUNE) +} // namespace lar_cluster3d diff --git a/dunereco/Cluster3D/cluster3dalgorithms_dune.fcl b/dunereco/Cluster3D/cluster3dalgorithms_dune.fcl new file mode 100644 index 00000000..3147127a --- /dev/null +++ b/dunereco/Cluster3D/cluster3dalgorithms_dune.fcl @@ -0,0 +1,28 @@ +# This defines the local DUNE version of the 3d space point builder tool for Cluster3D + +BEGIN_PROLOG + +standard_snippethit3dbuilderDUNE: +{ + tool_type: SnippetHit3DBuilderDUNE + HitFinderTagVec: ["gaushit"] + EnableMonitoring: true # enable monitoring of functions + HitWidthScaleFactor: 3.0 # + RangeNumSigma: 3.0 # + LongHitsStretchFactor: 1.5 # Allows to stretch long hits widths if desired + PulseHeightFraction: 0.4 # If multiple hits are under this fraction of the peak hit amplitude consider rejecting + PHLowSelection: 8. # If matching the above then require hits larger than this + MinPHFor2HitPoints: 4. # Don't consider hits less than this for 2 hit space points + DeltaPeakTimeSig: 1.75 # "Significance" of agreement between 2 hit peak times + ZPosOffset: 0.00 + MaxHitChiSquare: 6.00 + WirePitchScaleFactor: 1.9 + SaveBadChannelPairs: true # Create and save space points with one channel "bad" + SaveMythicalPoints: true # Create and save 2 hit combinations + MaxMythicalChiSquare: 25. # Cut on 2 hit points + UseT0Offsets: false + OutputHistograms: false + MakeAssociations: false +} + +END_PROLOG From 13ffe35898853dfbd0446213a558aaff081e0ec9 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Sun, 21 Dec 2025 13:00:44 -0800 Subject: [PATCH 08/10] Switch to DepoFluxWriter for the output of SimChannel info. Needed for SPINE. --- dunereco/DUNEWireCell/wirecell_dune.fcl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dunereco/DUNEWireCell/wirecell_dune.fcl b/dunereco/DUNEWireCell/wirecell_dune.fcl index 860d507c..a3c44c6c 100644 --- a/dunereco/DUNEWireCell/wirecell_dune.fcl +++ b/dunereco/DUNEWireCell/wirecell_dune.fcl @@ -710,9 +710,9 @@ dunefd_horizdrift_1x2x6_sim_nfsp : { configs: ["pgrapher/experiment/dune10kt-1x2x6/wcls-sim-drift-simchannel-nf-sp.jsonnet"] # Contract note: these exact "type:name" must be used to identify # the configuration data structures for these components in the Jsonnet. - inputers: ["wclsSimDepoSource:"] + inputers: ["wclsSimDepoSetSource:"] outputers: [ - "wclsSimChannelSink:postdrift", + "wclsDepoFluxWriter:postdrift", "wclsFrameSaver:spsignals" ] # Make available parameters via Jsonnet's std.extVar() From 2335c1114b55b1d523f6d707473d49cdefc72422 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Sun, 21 Dec 2025 16:11:16 -0800 Subject: [PATCH 09/10] Add the Cluster3D directory --- dunereco/CMakeLists.txt | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/dunereco/CMakeLists.txt b/dunereco/CMakeLists.txt index f11ff943..745c2faa 100644 --- a/dunereco/CMakeLists.txt +++ b/dunereco/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(AnaUtils) +add_subdirectory(Cluster3D) add_subdirectory(ClusterFinderDUNE) add_subdirectory(CVN) add_subdirectory(DUNEPandora) @@ -16,6 +17,32 @@ add_subdirectory(TransformerCVN) add_subdirectory(VLNets) add_subdirectory(PointResTree) add_subdirectory(WireFilter) + +# Supera +# +# Why are we checking explicitly for git submodule initialization? +# ---------------------------------------------------------------- +# Because it is stored as a git submodule, it should be initialized first before +# any build is attempted. `mrb gitCheckout` should run `--recurse-submodules` to +# this effect. However, it is implemented in a hacky way: the option for submodules +# is only added to `mrb gitCheckout` command if we are cloning from SBNSoftware/sbncode +# and nothing else. This means cloning from SBNSoftware/sbncode#somePRnumber will +# NOT carry this option. This type of cloning can happen in CI builds for example. +# In other words, we need to fix the hack with another hack: if Supera does not +# look initialized, we take the matter in our hands and run the git submodule +# initialization command. +# +# (This explanation above might be relevant for any other git submodules.) +if ( NOT EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/Supera/.git" ) + execute_process(COMMAND git submodule update --init --recursive + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) +endif() + +# We need to run Supera/setup.sh first, otherwise it won't have a CMakeLists.txt. +execute_process(COMMAND sh ${CMAKE_CURRENT_SOURCE_DIR}/Supera/setup.sh dune + WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/Supera ) + +add_subdirectory(Supera) # these two were not built in dunetpc at the time of the split # add_subdirectory(SNSlicer) # add_subdirectory(SNUtils) From d40bf6d9eae2dfe36d9b9e094683b599a94f0593 Mon Sep 17 00:00:00 2001 From: Tracy Usher Date: Mon, 22 Dec 2025 10:50:29 -0800 Subject: [PATCH 10/10] Need to declare packages that are needed in build (I am not sure why I needed to add "duneprototypes" to top level CMakeLists.txt file but not "larcv2") --- CMakeLists.txt | 1 + ups/product_deps | 14 ++++++++------ 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d4b39f6..d6ec3a44 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,6 +45,7 @@ find_package( canvas REQUIRED EXPORT) find_package(ROOT REQUIRED) find_package( dunecore REQUIRED EXPORT ) find_package( dunepdlegacy REQUIRED EXPORT ) +find_package( duneprototypes REQUIRED EXPORT ) find_package( larcore REQUIRED EXPORT ) find_package( larcoreobj REQUIRED EXPORT ) find_package( larcorealg REQUIRED EXPORT ) diff --git a/ups/product_deps b/ups/product_deps index 070f0eb6..a6d169c0 100644 --- a/ups/product_deps +++ b/ups/product_deps @@ -179,8 +179,10 @@ wpdir product_dir wire-cell-cfg # Add the dependent product and version product version -cetmodules v3_24_01 - only_for_build +cetmodules v3_24_01 - only_for_build dunecore v10_15_00d00 +duneprototypes v10_15_00d00 +larcv2 v2_2_6 end_product_list @@ -235,11 +237,11 @@ end_product_list # case it is optional. # #################################### -qualifier dunecore notes -c14:debug c14:debug -c14:prof c14:prof -e26:debug e26:debug -e26:prof e26:prof +qualifier dunecore duneprototypes larcv2 notes +c14:debug c14:debug c14:debug c14:debug:p3915 +c14:prof c14:prof c14:prof c14:prof:p3915 +e26:debug e26:debug e26:debug e26:debug:p3915 +e26:prof e26:prof e26:prof e26:prof:p3915 end_qualifier_list ####################################