From 62a73b082c975ee09c8581b4b7723dfd6da38c26 Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Tue, 30 Jun 2026 18:32:52 +0200 Subject: [PATCH 1/9] Addition of some WireMod diagnostic tools, addition of some flags (usage of SignalROIF, rounding of floats, fixed scaling factors to 1 at al. ) --- icaruscode/WireMod/CMakeLists.txt | 66 +- icaruscode/WireMod/WireModHitDiagnostic.fcl | 46 + .../WireMod/WireModHitDiagnostic_module.cc | 372 ++++++ .../WireMod/WireModScaleDebug_module.cc | 178 +++ icaruscode/WireMod/WireModifierXXW.fcl | 7 +- .../WireMod/WireModifierXXW_closure.fcl | 159 +++ icaruscode/WireMod/WireModifierXXW_module.cc | 1024 +++++++++++++++-- .../WireMod/stage1_icarus_wiremod_pops.fcl | 66 ++ .../stage1_icarus_wiremod_pops_hitdiag.fcl | 72 ++ .../WireMod/stage1_wiremod_2d_noscaling.fcl | 76 ++ .../WireMod/stage1_wiremod_2d_scaled.fcl | 76 ++ .../WireMod/stage1_wiremod_PT_noscaling.fcl | 91 ++ .../WireMod/stage1_wiremod_PT_scaled.fcl | 91 ++ .../WireMod/stage1_wiremod_diagnostic.fcl | 82 ++ .../stage1_wiremod_float_diagnostic.fcl | 73 ++ icaruscode/WireMod/wireMod_stage1.fcl | 47 + icaruscode/WireMod/wiremod_scale_debug.fcl | 26 + 17 files changed, 2422 insertions(+), 130 deletions(-) create mode 100644 icaruscode/WireMod/WireModHitDiagnostic.fcl create mode 100644 icaruscode/WireMod/WireModHitDiagnostic_module.cc create mode 100644 icaruscode/WireMod/WireModScaleDebug_module.cc create mode 100644 icaruscode/WireMod/WireModifierXXW_closure.fcl create mode 100644 icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl create mode 100644 icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_diagnostic.fcl create mode 100644 icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl create mode 100644 icaruscode/WireMod/wireMod_stage1.fcl create mode 100644 icaruscode/WireMod/wiremod_scale_debug.fcl diff --git a/icaruscode/WireMod/CMakeLists.txt b/icaruscode/WireMod/CMakeLists.txt index f81dbfe53..1c605d42f 100644 --- a/icaruscode/WireMod/CMakeLists.txt +++ b/icaruscode/WireMod/CMakeLists.txt @@ -1,4 +1,4 @@ -set(EXCLUDED_PLUGINS WireModMakeHists_plugin.cc ROICounter_module.cc HitComp_plugin.cc HitTTreeMaker_plugin.cc WireModifier_module.cc WireModifierXXW_module.cc) +set(EXCLUDED_PLUGINS WireModMakeHists_plugin.cc ROICounter_module.cc HitComp_plugin.cc HitTTreeMaker_plugin.cc WireModifier_module.cc WireModifierXXW_module.cc WireModHitDiagnostic_module.cc WireModScaleDebug_module.cc) art_make( LIB_LIBRARIES @@ -241,6 +241,70 @@ simple_plugin(WireModifierXXW "module" BASENAME_ONLY ) +simple_plugin(WireModHitDiagnostic "module" + ${ART_FRAMEWORK_CORE} + ${MF_MESSAGELOGGER} + ${Boost_SYSTEM_LIBRARY} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + art_root_io::TFileService_service + lardataobj::AnalysisBase + lardataobj::RawData + lardataobj::RecoBase + lardata::RecoObjects + larreco::RecoAlg + ${LARRECO_LIB} + ${LARDATA_LIB} + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${MF_MESSAGELOGGER} + ${FHICLCPP} + ${CLHEP} + ${ROOT_GEOM} + ${ROOT_XMLIO} + ${ROOT_GDML} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + BASENAME_ONLY +) + +simple_plugin(WireModScaleDebug "module" + ${ART_FRAMEWORK_CORE} + ${MF_MESSAGELOGGER} + ${Boost_SYSTEM_LIBRARY} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + art_root_io::TFileService_service + lardataobj::AnalysisBase + lardataobj::RawData + lardataobj::RecoBase + lardata::RecoObjects + larreco::RecoAlg + ${LARRECO_LIB} + ${LARDATA_LIB} + ${ART_FRAMEWORK_CORE} + ${ART_FRAMEWORK_PRINCIPAL} + ${ART_FRAMEWORK_SERVICES_REGISTRY} + ${ART_FRAMEWORK_SERVICES_OPTIONAL} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_RANDOMNUMBERGENERATOR_SERVICE} + ${ART_FRAMEWORK_SERVICES_OPTIONAL_TFILESERVICE_SERVICE} + ${MF_MESSAGELOGGER} + ${FHICLCPP} + ${CLHEP} + ${ROOT_GEOM} + ${ROOT_XMLIO} + ${ROOT_GDML} + ${ROOT_BASIC_LIB_LIST} + ${ROOT_HIST} + larcore::Geometry_Geometry_service + larcorealg::Geometry + BASENAME_ONLY +) + install_headers() install_fhicl() install_source() diff --git a/icaruscode/WireMod/WireModHitDiagnostic.fcl b/icaruscode/WireMod/WireModHitDiagnostic.fcl new file mode 100644 index 000000000..7e9e5eb19 --- /dev/null +++ b/icaruscode/WireMod/WireModHitDiagnostic.fcl @@ -0,0 +1,46 @@ +BEGIN_PROLOG + +WireModHitDiag: +{ + module_type: "WireModHitDiagnostic" + + BeforeHitLabels: [ + "dummyGaushit2dTPCEE", + "dummyGaushit2dTPCEW", + "dummyGaushit2dTPCWE", + "dummyGaushit2dTPCWW" + ] + + AfterHitLabels: [ + "gaushit2dTPCEE", + "gaushit2dTPCEW", + "gaushit2dTPCWE", + "gaushit2dTPCWW" + ] + + ScaleInfoLabels: [ + "WireModTPCEE", + "WireModTPCEW", + "WireModTPCWE", + "WireModTPCWW" + ] + + MatchWindowTicks: 10.0 + OnlyMC: true +} + +WireModROITree: { + SaveROITree: true + SaveTickHistograms: true + SavePerHitData: true + SaveByChannel: true +} + +WireModPerHitData: { + SaveROITree: false + SaveTickHistograms: false + SavePerHitData: true + SaveByChannel: false +} + +END_PROLOG diff --git a/icaruscode/WireMod/WireModHitDiagnostic_module.cc b/icaruscode/WireMod/WireModHitDiagnostic_module.cc new file mode 100644 index 000000000..1e1b5b45b --- /dev/null +++ b/icaruscode/WireMod/WireModHitDiagnostic_module.cc @@ -0,0 +1,372 @@ +// std includes +#include +#include +#include +#include + +// ROOT includes +#include "TTree.h" + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileService.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +// larsoft includes +#include "lardataobj/RecoBase/Hit.h" +#include "larcore/Geometry/WireReadout.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" +#include "larcore/CoreUtils/ServiceUtil.h" + +namespace wiremod { + +class WireModHitDiagnostic : public art::EDAnalyzer +{ +public: + explicit WireModHitDiagnostic(fhicl::ParameterSet const& pset); + void analyze(art::Event const& evt) override; + +private: + std::vector fBeforeHitLabels; // dummyGaushit2dTPC* + std::vector fAfterHitLabels; // gaushit2dTPC* + std::vector fScaleInfoLabels; // WireModTPC* (produce scale products) + float fMatchWindowTicks; + bool fOnlyMC; // if true, skip overlay hits (is_mc==0) when filling HitMatchTree + + const geo::WireReadoutGeom* fWireReadout = + &(art::ServiceHandle()->Get()); + + // HitMatchTree: one entry per before-hit + TTree* fHitMatchTree; + Int_t tm_run, tm_subrun, tm_event; + Int_t tm_channel, tm_cryo, tm_tpc, tm_plane, tm_wire; + Float_t tm_before_peak_time, tm_before_peak_amp, tm_before_integral, tm_before_rms; + Int_t tm_before_start_tick, tm_before_end_tick; + Float_t tm_before_gof; + Int_t tm_before_dof, tm_before_n_multi; + // MC / overlay flag from WireModifierXXW BackTracker + Int_t tm_is_mc; // 1 = MC hit, 0 = overlay data hit + // scale info from WireModifierXXW (-999 if overlay / no EDep match) + Float_t tm_scale_q; + Float_t tm_scale_sigma; + Float_t tm_truth_x; // truth X coordinate used for graph lookup (cm) + Float_t tm_theta_xw; // ThetaXW angle used for graph lookup (degrees) + Float_t tm_truth_dirx; // truth dxdr + Float_t tm_truth_diry; // truth dydr + Float_t tm_truth_dirz; // truth dzdr + Float_t tm_pitch; // wirePitch / cosG [cm] + Float_t tm_dqdx_before; // integral_before / pitch [ADC/cm], no lifetime correction + // match result + Int_t tm_is_matched; + Int_t tm_n_after_matches; + Float_t tm_after_peak_time, tm_after_peak_amp, tm_after_integral, tm_after_rms; + Int_t tm_after_start_tick, tm_after_end_tick; + Float_t tm_after_gof; + Int_t tm_after_dof; + Float_t tm_ratio_integral; + Float_t tm_ratio_rms; + + // NewHitTree: one entry per after-hit without a before-hit match + TTree* fNewHitTree; + Int_t tn_run, tn_subrun, tn_event; + Int_t tn_channel, tn_cryo, tn_tpc, tn_plane, tn_wire; + Float_t tn_peak_time, tn_peak_amp, tn_integral, tn_rms; + Int_t tn_start_tick, tn_end_tick; + Float_t tn_gof; + Int_t tn_dof; +}; + +//---------------------------------------------------------------------------- +WireModHitDiagnostic::WireModHitDiagnostic(fhicl::ParameterSet const& pset) + : EDAnalyzer{pset} + , fMatchWindowTicks{pset.get("MatchWindowTicks", 10.f)} + , fOnlyMC{pset.get("OnlyMC", false)} +{ + for (auto const& s : pset.get>("BeforeHitLabels")) + fBeforeHitLabels.emplace_back(s); + for (auto const& s : pset.get>("AfterHitLabels")) + fAfterHitLabels.emplace_back(s); + for (auto const& s : pset.get>("ScaleInfoLabels", + std::vector{})) + fScaleInfoLabels.emplace_back(s); + + art::ServiceHandle tfs; + + fHitMatchTree = tfs->make("HitMatchTree", + "WireMod before-after hit comparison"); + fHitMatchTree->Branch("run", &tm_run, "run/I"); + fHitMatchTree->Branch("subrun", &tm_subrun, "subrun/I"); + fHitMatchTree->Branch("event", &tm_event, "event/I"); + fHitMatchTree->Branch("channel", &tm_channel, "channel/I"); + fHitMatchTree->Branch("cryo", &tm_cryo, "cryo/I"); + fHitMatchTree->Branch("tpc", &tm_tpc, "tpc/I"); + fHitMatchTree->Branch("plane", &tm_plane, "plane/I"); + fHitMatchTree->Branch("wire", &tm_wire, "wire/I"); + fHitMatchTree->Branch("before_peak_time", &tm_before_peak_time, "before_peak_time/F"); + fHitMatchTree->Branch("before_peak_amp", &tm_before_peak_amp, "before_peak_amp/F"); + fHitMatchTree->Branch("before_integral", &tm_before_integral, "before_integral/F"); + fHitMatchTree->Branch("before_rms", &tm_before_rms, "before_rms/F"); + fHitMatchTree->Branch("before_start_tick", &tm_before_start_tick, "before_start_tick/I"); + fHitMatchTree->Branch("before_end_tick", &tm_before_end_tick, "before_end_tick/I"); + fHitMatchTree->Branch("before_gof", &tm_before_gof, "before_gof/F"); + fHitMatchTree->Branch("before_dof", &tm_before_dof, "before_dof/I"); + fHitMatchTree->Branch("before_n_multi", &tm_before_n_multi, "before_n_multi/I"); + fHitMatchTree->Branch("is_mc", &tm_is_mc, "is_mc/I"); + fHitMatchTree->Branch("scale_q", &tm_scale_q, "scale_q/F"); + fHitMatchTree->Branch("scale_sigma", &tm_scale_sigma, "scale_sigma/F"); + fHitMatchTree->Branch("truth_x", &tm_truth_x, "truth_x/F"); + fHitMatchTree->Branch("theta_xw", &tm_theta_xw, "theta_xw/F"); + fHitMatchTree->Branch("truth_dirx", &tm_truth_dirx, "truth_dirx/F"); + fHitMatchTree->Branch("truth_diry", &tm_truth_diry, "truth_diry/F"); + fHitMatchTree->Branch("truth_dirz", &tm_truth_dirz, "truth_dirz/F"); + fHitMatchTree->Branch("pitch", &tm_pitch, "pitch/F"); + fHitMatchTree->Branch("dqdx_before", &tm_dqdx_before, "dqdx_before/F"); + fHitMatchTree->Branch("is_matched", &tm_is_matched, "is_matched/I"); + fHitMatchTree->Branch("n_after_matches", &tm_n_after_matches, "n_after_matches/I"); + fHitMatchTree->Branch("after_peak_time", &tm_after_peak_time, "after_peak_time/F"); + fHitMatchTree->Branch("after_peak_amp", &tm_after_peak_amp, "after_peak_amp/F"); + fHitMatchTree->Branch("after_integral", &tm_after_integral, "after_integral/F"); + fHitMatchTree->Branch("after_rms", &tm_after_rms, "after_rms/F"); + fHitMatchTree->Branch("after_start_tick", &tm_after_start_tick, "after_start_tick/I"); + fHitMatchTree->Branch("after_end_tick", &tm_after_end_tick, "after_end_tick/I"); + fHitMatchTree->Branch("after_gof", &tm_after_gof, "after_gof/F"); + fHitMatchTree->Branch("after_dof", &tm_after_dof, "after_dof/I"); + fHitMatchTree->Branch("ratio_integral", &tm_ratio_integral, "ratio_integral/F"); + fHitMatchTree->Branch("ratio_rms", &tm_ratio_rms, "ratio_rms/F"); + + fNewHitTree = tfs->make("NewHitTree", + "After-mod hits with no before-mod counterpart"); + fNewHitTree->Branch("run", &tn_run, "run/I"); + fNewHitTree->Branch("subrun", &tn_subrun, "subrun/I"); + fNewHitTree->Branch("event", &tn_event, "event/I"); + fNewHitTree->Branch("channel", &tn_channel, "channel/I"); + fNewHitTree->Branch("cryo", &tn_cryo, "cryo/I"); + fNewHitTree->Branch("tpc", &tn_tpc, "tpc/I"); + fNewHitTree->Branch("plane", &tn_plane, "plane/I"); + fNewHitTree->Branch("wire", &tn_wire, "wire/I"); + fNewHitTree->Branch("peak_time", &tn_peak_time, "peak_time/F"); + fNewHitTree->Branch("peak_amp", &tn_peak_amp, "peak_amp/F"); + fNewHitTree->Branch("integral", &tn_integral, "integral/F"); + fNewHitTree->Branch("rms", &tn_rms, "rms/F"); + fNewHitTree->Branch("start_tick", &tn_start_tick, "start_tick/I"); + fNewHitTree->Branch("end_tick", &tn_end_tick, "end_tick/I"); + fNewHitTree->Branch("gof", &tn_gof, "gof/F"); + fNewHitTree->Branch("dof", &tn_dof, "dof/I"); +} + +// global TPC index: 0=EE (cryo0,tpcset0) 1=EW (cryo0,tpcset1) +// 2=WE (cryo1,tpcset0) 3=WW (cryo1,tpcset1) +int globalTPC(raw::ChannelID_t ch, geo::WireReadoutGeom const* wr) +{ + auto rop = wr->ChannelToROP(ch); + return static_cast(rop.Cryostat) * 2 + static_cast(rop.TPCset); +} + +//---------------------------------------------------------------------------- +void WireModHitDiagnostic::analyze(art::Event const& evt) +{ + tm_run = tn_run = static_cast(evt.run()); + tm_subrun = tn_subrun = static_cast(evt.subRun()); + tm_event = tn_event = static_cast(evt.event()); + + // Gather all after-hits once (searched by channel for every before-hit) + std::vector afterHits; + for (auto const& tag : fAfterHitLabels) + { + auto h = evt.getHandle>(tag); + if (!h.isValid()) + { + mf::LogWarning("WireModHitDiagnostic") + << "After-hit collection not found: " << tag; + continue; + } + for (auto const& hit : *h) + afterHits.push_back(&hit); + } + + // Build channel -> [after-hit indices] index + std::unordered_multimap afterByChannel; + afterByChannel.reserve(afterHits.size() * 2); + for (size_t i = 0; i < afterHits.size(); ++i) + afterByChannel.emplace(afterHits[i]->Channel(), i); + + std::vector afterMatched(afterHits.size(), false); + + // Process before-hit collections one at a time to keep the hit index + // aligned with the scale-info parallel vectors from WireModifierXXW. + for (size_t j = 0; j < fBeforeHitLabels.size(); ++j) + { + auto bHandle = evt.getHandle>(fBeforeHitLabels[j]); + if (!bHandle.isValid()) + { + mf::LogWarning("WireModHitDiagnostic") + << "Before-hit collection not found: " << fBeforeHitLabels[j]; + continue; + } + auto const& localHits = *bHandle; + + // Try to load scale-info and isMC parallel vectors for this TPC + std::vector const* pScaleQ = nullptr; + std::vector const* pScaleSigma = nullptr; + std::vector const* pTruthX = nullptr; + std::vector const* pThetaXW = nullptr; + std::vector const* pIsMC = nullptr; + std::vector const* pDirX = nullptr; + std::vector const* pDirY = nullptr; + std::vector const* pDirZ = nullptr; + std::vector const* pPitch = nullptr; + std::vector const* pDQdX = nullptr; + + if (j < fScaleInfoLabels.size()) + { + auto const& siLabel = fScaleInfoLabels[j]; + auto sqH = evt.getHandle>(art::InputTag(siLabel.label(), "scaleQ")); + auto ssH = evt.getHandle>(art::InputTag(siLabel.label(), "scaleSigma")); + auto txH = evt.getHandle>(art::InputTag(siLabel.label(), "truthX")); + auto thH = evt.getHandle>(art::InputTag(siLabel.label(), "thetaXW")); + auto imH = evt.getHandle> (art::InputTag(siLabel.label(), "isMC")); + auto dxH = evt.getHandle>(art::InputTag(siLabel.label(), "dirX")); + auto dyH = evt.getHandle>(art::InputTag(siLabel.label(), "dirY")); + auto dzH = evt.getHandle>(art::InputTag(siLabel.label(), "dirZ")); + auto ptH = evt.getHandle>(art::InputTag(siLabel.label(), "pitch")); + auto dqH = evt.getHandle>(art::InputTag(siLabel.label(), "dQdX")); + if (sqH.isValid()) pScaleQ = &(*sqH); + if (ssH.isValid()) pScaleSigma = &(*ssH); + if (txH.isValid()) pTruthX = &(*txH); + if (thH.isValid()) pThetaXW = &(*thH); + if (imH.isValid()) pIsMC = &(*imH); + if (dxH.isValid()) pDirX = &(*dxH); + if (dyH.isValid()) pDirY = &(*dyH); + if (dzH.isValid()) pDirZ = &(*dzH); + if (ptH.isValid()) pPitch = &(*ptH); + if (dqH.isValid()) pDQdX = &(*dqH); + } + + for (size_t i = 0; i < localHits.size(); ++i) + { + auto const& bh = localHits[i]; + auto const& wid = bh.WireID(); + + tm_channel = static_cast(bh.Channel()); + tm_cryo = static_cast(wid.Cryostat); + tm_tpc = globalTPC(bh.Channel(), fWireReadout); + tm_plane = static_cast(wid.Plane); + tm_wire = static_cast(wid.Wire); + tm_before_peak_time = bh.PeakTime(); + tm_before_peak_amp = bh.PeakAmplitude(); + tm_before_integral = bh.Integral(); + tm_before_rms = bh.RMS(); + tm_before_start_tick = bh.StartTick(); + tm_before_end_tick = bh.EndTick(); + tm_before_gof = bh.GoodnessOfFit(); + tm_before_dof = bh.DegreesOfFreedom(); + tm_before_n_multi = bh.Multiplicity(); + + // MC flag and scale info (parallel vectors indexed by hit position in this collection) + tm_is_mc = (pIsMC && i < pIsMC->size()) ? (*pIsMC)[i] : -1; + tm_scale_q = (pScaleQ && i < pScaleQ->size()) ? (*pScaleQ)[i] : -999.f; + tm_scale_sigma = (pScaleSigma && i < pScaleSigma->size()) ? (*pScaleSigma)[i] : -999.f; + tm_truth_x = (pTruthX && i < pTruthX->size()) ? (*pTruthX)[i] : -999.f; + tm_theta_xw = (pThetaXW && i < pThetaXW->size()) ? (*pThetaXW)[i] : -999.f; + tm_truth_dirx = (pDirX && i < pDirX->size()) ? (*pDirX)[i] : -999.f; + tm_truth_diry = (pDirY && i < pDirY->size()) ? (*pDirY)[i] : -999.f; + tm_truth_dirz = (pDirZ && i < pDirZ->size()) ? (*pDirZ)[i] : -999.f; + tm_pitch = (pPitch && i < pPitch->size()) ? (*pPitch)[i] : -999.f; + tm_dqdx_before = (pDQdX && i < pDQdX->size()) ? (*pDQdX)[i] : -999.f; + + // if OnlyMC is set, skip overlay hits + if (fOnlyMC && tm_is_mc != 1) continue; + + // Search for matching after-hit on the same channel + float winLo = static_cast(bh.StartTick()) - fMatchWindowTicks; + float winHi = static_cast(bh.EndTick()) + fMatchWindowTicks; + + size_t bestIdx = std::numeric_limits::max(); + float bestDist = std::numeric_limits::max(); + int nMatch = 0; + + auto range = afterByChannel.equal_range(bh.Channel()); + for (auto it = range.first; it != range.second; ++it) + { + float aPeak = afterHits[it->second]->PeakTime(); + if (aPeak < winLo || aPeak > winHi) continue; + ++nMatch; + float dist = std::abs(aPeak - bh.PeakTime()); + if (dist < bestDist) { bestDist = dist; bestIdx = it->second; } + } + + tm_n_after_matches = nMatch; + tm_is_matched = (bestIdx != std::numeric_limits::max()) ? 1 : 0; + + if (tm_is_matched) + { + for (auto it = range.first; it != range.second; ++it) + { + float aPeak = afterHits[it->second]->PeakTime(); + if (aPeak >= winLo && aPeak <= winHi) + afterMatched[it->second] = true; + } + auto const* ah = afterHits[bestIdx]; + tm_after_peak_time = ah->PeakTime(); + tm_after_peak_amp = ah->PeakAmplitude(); + tm_after_integral = ah->Integral(); + tm_after_rms = ah->RMS(); + tm_after_start_tick = ah->StartTick(); + tm_after_end_tick = ah->EndTick(); + tm_after_gof = ah->GoodnessOfFit(); + tm_after_dof = ah->DegreesOfFreedom(); + tm_ratio_integral = (bh.Integral() > 0.f) + ? ah->Integral() / bh.Integral() : -999.f; + tm_ratio_rms = (bh.RMS() > 0.f) + ? ah->RMS() / bh.RMS() : -999.f; + } + else + { + tm_after_peak_time = -999.f; + tm_after_peak_amp = -999.f; + tm_after_integral = -999.f; + tm_after_rms = -999.f; + tm_after_start_tick = -999; + tm_after_end_tick = -999; + tm_after_gof = -999.f; + tm_after_dof = -999; + tm_ratio_integral = -999.f; + tm_ratio_rms = -999.f; + } + + fHitMatchTree->Fill(); + } + } + + // NewHitTree: after-hits not claimed by any before-hit + for (size_t i = 0; i < afterHits.size(); ++i) + { + if (afterMatched[i]) continue; + + auto const* ah = afterHits[i]; + auto const& wid = ah->WireID(); + tn_channel = static_cast(ah->Channel()); + tn_cryo = static_cast(wid.Cryostat); + tn_tpc = globalTPC(ah->Channel(), fWireReadout); + tn_plane = static_cast(wid.Plane); + tn_wire = static_cast(wid.Wire); + tn_peak_time = ah->PeakTime(); + tn_peak_amp = ah->PeakAmplitude(); + tn_integral = ah->Integral(); + tn_rms = ah->RMS(); + tn_start_tick = ah->StartTick(); + tn_end_tick = ah->EndTick(); + tn_gof = ah->GoodnessOfFit(); + tn_dof = ah->DegreesOfFreedom(); + fNewHitTree->Fill(); + } +} + +DEFINE_ART_MODULE(WireModHitDiagnostic) + +} // namespace wiremod diff --git a/icaruscode/WireMod/WireModScaleDebug_module.cc b/icaruscode/WireMod/WireModScaleDebug_module.cc new file mode 100644 index 000000000..6f4b4f5c4 --- /dev/null +++ b/icaruscode/WireMod/WireModScaleDebug_module.cc @@ -0,0 +1,178 @@ +// std includes +#include +#include +#include +#include + +// ROOT includes +#include "TFile.h" +#include "TGraph2D.h" +#include "TMath.h" + +// art includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +// larsoft includes +#include "larcore/CoreUtils/ServiceUtil.h" +#include "larcore/Geometry/Geometry.h" +#include "larcore/Geometry/WireReadout.h" +#include "larcorealg/Geometry/WireReadoutGeom.h" + +namespace wiremod { + +class WireModScaleDebug : public art::EDAnalyzer { +public: + explicit WireModScaleDebug(fhicl::ParameterSet const& pset); + void beginJob() override; + void analyze(art::Event const&) override {} + +private: + std::string fRatioFileName_XXW; + std::vector fThetaEdges; + + static void shapeGraphAngle(TGraph2D& graph) + { + double* xPtr = graph.GetX(); + double* yPtr = graph.GetY(); + double* zPtr = graph.GetZ(); + for (int idx = 0; idx < graph.GetN(); ++idx) { + double x = *(xPtr + idx); + double y = *(yPtr + idx); + double z = *(zPtr + idx); + graph.SetPoint(idx, x, TMath::DegToRad() * y, z); + } + } + + static void shapeGraphPosNegate(TGraph2D& graph, bool negate) + { + if (!negate) return; + double* xPtr = graph.GetX(); + double* yPtr = graph.GetY(); + double* zPtr = graph.GetZ(); + for (int idx = 0; idx < graph.GetN(); ++idx) { + double x = *(xPtr + idx); + double y = *(yPtr + idx); + double z = *(zPtr + idx); + graph.SetPoint(idx, -x, y, z); + } + } +}; + +WireModScaleDebug::WireModScaleDebug(fhicl::ParameterSet const& pset) + : EDAnalyzer{pset} + , fRatioFileName_XXW{pset.get("RatioFileName_XXW", "Run4_Coarse_XW.root")} + , fThetaEdges{pset.get>("ThetaEdges", + {0, 6, 12, 18, 24, 30, 38, 46, 56, 90})} +{ +} + +void WireModScaleDebug::beginJob() +{ + const geo::GeometryCore* geometry = lar::providerFrom(); + const geo::WireReadoutGeom* wireReadout = &(art::ServiceHandle()->Get()); + + TFile* ratioFile = new TFile(("./" + fRatioFileName_XXW).c_str(), "READ"); + if (!ratioFile || ratioFile->IsZombie()) { + std::cout << "ERROR: could not open " << fRatioFileName_XXW + << " in current directory. Run lar from the directory containing it.\n"; + return; + } + + struct TPCSetInfo { unsigned cryo; unsigned tpcset; unsigned globalTPC; std::string name; }; + std::vector tpcsets = { + {0, 0, 0, "EE"}, {0, 1, 1, "EW"}, {1, 0, 2, "WE"}, {1, 1, 3, "WW"} + }; + + std::cout + << "TPC,plane,cathode_X,test_X,theta_xw_deg,theta_bin," + << "interp_raw_q,found_scale_q,interp_raw_sigma,found_scale_sigma," + << "nearest_node_theta_deg,expected_scale_q,expected_scale_sigma\n"; + + for (auto const& tpcinfo : tpcsets) { + readout::TPCsetID tpcsetid(tpcinfo.cryo, tpcinfo.tpcset); + geo::TPCID tpcid = wireReadout->TPCsetToTPCs(tpcsetid).front(); + const geo::TPCGeo* tpcGeom = geometry->TPCPtr(tpcid); + double cathodeX = tpcGeom->GetCathodeCenter().X(); + bool negate = (cathodeX <= 0); + + for (int plane = 0; plane < 3; ++plane) { + std::string nameQ = "TPC" + std::to_string(tpcinfo.globalTPC) + "_plane" + + std::to_string(plane) + "_ratio_integral"; + std::string nameS = "TPC" + std::to_string(tpcinfo.globalTPC) + "_plane" + + std::to_string(plane) + "_ratio_width"; + + TGraph2D* gQ_orig = dynamic_cast(ratioFile->Get(nameQ.c_str())); + TGraph2D* gS_orig = dynamic_cast(ratioFile->Get(nameS.c_str())); + if (!gQ_orig || !gS_orig) { + std::cout << "MISSING graph " << nameQ << " or " << nameS << "\n"; + continue; + } + TGraph2D* gQ = static_cast(gQ_orig->Clone()); + TGraph2D* gS = static_cast(gS_orig->Clone()); + + // mirrors WireModifierXXW_module.cc reconfigure() (InRadians:false, XAbs:true) + shapeGraphAngle(*gQ); + shapeGraphAngle(*gS); + shapeGraphPosNegate(*gQ, negate); + shapeGraphPosNegate(*gS, negate); + + double xmin = gQ->GetXmin(), xmax = gQ->GetXmax(); + double xmag = 0.5 * (std::abs(xmin) + std::abs(xmax)); + double testX = negate ? -xmag : xmag; + + double* xPtrQ = gQ->GetX(); + double* yPtrQ = gQ->GetY(); + double* zPtrQ = gQ->GetZ(); + double* zPtrS = gS->GetZ(); + int nPts = gQ->GetN(); + + for (double theta_deg = 0.0; theta_deg <= 90.0; theta_deg += 2.0) { + double theta_rad = theta_deg * TMath::DegToRad(); + + double raw_q = gQ->Interpolate(testX, theta_rad); + double raw_s = gS->Interpolate(testX, theta_rad); + // production threshold check (WireModUtility.cc): > 0.001 else default 1.0 + double found_q = (raw_q > 0.001) ? raw_q : 1.0; + double found_s = (raw_s > 0.001) ? raw_s : 1.0; + + // nearest real calibration grid node (2D nearest neighbor) + double bestD = 1e18; + int bestIdx = -1; + for (int i = 0; i < nPts; ++i) { + double dx = xPtrQ[i] - testX; + double dy = yPtrQ[i] - theta_rad; + double d = dx * dx * 1e-4 + dy * dy; // theta-dominated distance (X fixed at midpoint) + if (d < bestD) { bestD = d; bestIdx = i; } + } + double nearest_theta_deg = yPtrQ[bestIdx] * TMath::RadToDeg(); + double expected_q = zPtrQ[bestIdx]; + double expected_s = zPtrS[bestIdx]; + + int thetaBin = -1; + for (size_t b = 0; b + 1 < fThetaEdges.size(); ++b) + if (theta_deg >= fThetaEdges[b] && theta_deg < fThetaEdges[b + 1]) { thetaBin = (int)b; break; } + + std::cout << std::fixed + << tpcinfo.name << "," << plane << "," + << std::setprecision(1) << cathodeX << "," + << std::setprecision(1) << testX << "," + << std::setprecision(1) << theta_deg << "," << thetaBin << "," + << std::setprecision(4) << raw_q << "," << found_q << "," + << raw_s << "," << found_s << "," + << std::setprecision(1) << nearest_theta_deg << "," + << std::setprecision(4) << expected_q << "," << expected_s + << "\n"; + } + } + } + + ratioFile->Close(); +} + +DEFINE_ART_MODULE(WireModScaleDebug) + +} // namespace wiremod diff --git a/icaruscode/WireMod/WireModifierXXW.fcl b/icaruscode/WireMod/WireModifierXXW.fcl index 1598479ce..56abd6722 100644 --- a/icaruscode/WireMod/WireModifierXXW.fcl +++ b/icaruscode/WireMod/WireModifierXXW.fcl @@ -34,8 +34,11 @@ WireModXXW: LocalRatios: true EDepLabel: "filtersed" EDepShiftedLabel: "shiftedpriorSCE" - SaveByChannel: true - SaveByWire: false + SaveByChannel: false + SaveByWire: false + SaveROITree: false + SaveTickHistograms: false + SavePerHitData: false InRadians: false XAbs: true } diff --git a/icaruscode/WireMod/WireModifierXXW_closure.fcl b/icaruscode/WireMod/WireModifierXXW_closure.fcl new file mode 100644 index 000000000..fba253031 --- /dev/null +++ b/icaruscode/WireMod/WireModifierXXW_closure.fcl @@ -0,0 +1,159 @@ +#include "stage1_icarus_defs.fcl" +#include "wirechannelroiconverters_sbn.fcl" + +BEGIN_PROLOG + +dummyGaushit2dTPCEE: { @table::icarus_stage1_producers.gaushit2dTPCEE } +dummyGaushit2dTPCEW: { @table::icarus_stage1_producers.gaushit2dTPCEW } +dummyGaushit2dTPCWE: { @table::icarus_stage1_producers.gaushit2dTPCWE } +dummyGaushit2dTPCWW: { @table::icarus_stage1_producers.gaushit2dTPCWW } +dummyGaushit2dSeq: [dummyGaushit2dTPCEE, dummyGaushit2dTPCEW, dummyGaushit2dTPCWE, dummyGaushit2dTPCWW] +dummyGaushit2dTab: { dummyGaushit2dTPCEE: @local::dummyGaushit2dTPCEE + dummyGaushit2dTPCEW: @local::dummyGaushit2dTPCEW + dummyGaushit2dTPCWE: @local::dummyGaushit2dTPCWE + dummyGaushit2dTPCWW: @local::dummyGaushit2dTPCWW } + +WMROI2Wire: +{ + module_type: "ChannelROIToWire" + WireModuleLabelVec: [ "wire2channelroi2d:PHYSCRATEDATATPCEE", + "wire2channelroi2d:PHYSCRATEDATATPCEW", + "wire2channelroi2d:PHYSCRATEDATATPCWE", + "wire2channelroi2d:PHYSCRATEDATATPCWW" ] + OutInstanceLabelVec: [ "PHYSCRATEDATATPCEE", + "PHYSCRATEDATATPCEW", + "PHYSCRATEDATATPCWE", + "PHYSCRATEDATATPCWW" ] + DiagnosticOutput: false + UseFloatPrecision: true # avoid std::round() to 1 ADC steps +} + +WireMod: +{ + module_type: "WireModifierXXW" + RatioFileName_XXW: "Run4_Coarse_XW.root" + RatioFileName_YZ: "allRatio_Coarse_XW.root" + LocalRatios: true + EDepLabel: "filtersed" + EDepShiftedLabel: "shiftedpriorSCE" + SaveByChannel: false + SaveByWire: false + SaveROITree: false + SaveTickHistograms: false + SavePerHitData: false + InRadians: false + XAbs: true + UseChannelROIMode: false # true = read ChannelROI directly, bypassing Wire rounding + AdditiveModification: false # true = after = before + (q_mod - q_orig) + + SetNullScaleIntegral: false # true = force r_Q=1 + SetNullScaleWidth: false # true = force r_sigma=1 +} + +WireModTPCEE: @local::WireMod +WireModTPCEE.Cryo: 0 +WireModTPCEE.TPCset: 0 +WireModTPCEE.XXWScaleIntegral: [ "TPC0_plane0_ratio_integral", + "TPC0_plane1_ratio_integral", + "TPC0_plane2_ratio_integral"] +WireModTPCEE.XXWScaleWidth: [ "TPC0_plane0_ratio_width", + "TPC0_plane1_ratio_width", + "TPC0_plane2_ratio_width"] +WireModTPCEE.YZScaleIntegral: [ "TPC0_plane0_ratio_integral", + "TPC0_plane1_ratio_integral", + "TPC0_plane2_ratio_integral"] +WireModTPCEE.YZScaleWidth: [ "TPC0_plane0_ratio_width", + "TPC0_plane1_ratio_width", + "TPC0_plane2_ratio_width"] +WireModTPCEE.WireLabel: "WMROI2Wire:PHYSCRATEDATATPCEE" +WireModTPCEE.ChannelLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +WireModTPCEE.HitLabel: "dummyGaushit2dTPCEE" + +WireModTPCEW: @local::WireMod +WireModTPCEW.Cryo: 0 +WireModTPCEW.TPCset: 1 +WireModTPCEW.XXWScaleIntegral: [ "TPC1_plane0_ratio_integral", + "TPC1_plane1_ratio_integral", + "TPC1_plane2_ratio_integral"] +WireModTPCEW.XXWScaleWidth: [ "TPC1_plane0_ratio_width", + "TPC1_plane1_ratio_width", + "TPC1_plane2_ratio_width"] +WireModTPCEW.YZScaleIntegral: [ "TPC1_plane0_ratio_integral", + "TPC1_plane1_ratio_integral", + "TPC1_plane2_ratio_integral"] +WireModTPCEW.YZScaleWidth: [ "TPC1_plane0_ratio_width", + "TPC1_plane1_ratio_width", + "TPC1_plane2_ratio_width"] +WireModTPCEW.WireLabel: "WMROI2Wire:PHYSCRATEDATATPCEW" +WireModTPCEW.ChannelLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +WireModTPCEW.HitLabel: "dummyGaushit2dTPCEW" + +WireModTPCWE: @local::WireMod +WireModTPCWE.Cryo: 1 +WireModTPCWE.TPCset: 0 +WireModTPCWE.XXWScaleIntegral: [ "TPC2_plane0_ratio_integral", + "TPC2_plane1_ratio_integral", + "TPC2_plane2_ratio_integral"] +WireModTPCWE.XXWScaleWidth: [ "TPC2_plane0_ratio_width", + "TPC2_plane1_ratio_width", + "TPC2_plane2_ratio_width"] +WireModTPCWE.YZScaleIntegral: [ "TPC2_plane0_ratio_integral", + "TPC2_plane1_ratio_integral", + "TPC2_plane2_ratio_integral"] +WireModTPCWE.YZScaleWidth: [ "TPC2_plane0_ratio_width", + "TPC2_plane1_ratio_width", + "TPC2_plane2_ratio_width"] +WireModTPCWE.WireLabel: "WMROI2Wire:PHYSCRATEDATATPCWE" +WireModTPCWE.ChannelLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +WireModTPCWE.HitLabel: "dummyGaushit2dTPCWE" + +WireModTPCWW: @local::WireMod +WireModTPCWW.Cryo: 1 +WireModTPCWW.TPCset: 1 +WireModTPCWW.XXWScaleIntegral: [ "TPC3_plane0_ratio_integral", + "TPC3_plane1_ratio_integral", + "TPC3_plane2_ratio_integral"] +WireModTPCWW.XXWScaleWidth: [ "TPC3_plane0_ratio_width", + "TPC3_plane1_ratio_width", + "TPC3_plane2_ratio_width"] +WireModTPCWW.YZScaleIntegral: [ "TPC3_plane0_ratio_integral", + "TPC3_plane1_ratio_integral", + "TPC3_plane2_ratio_integral"] +WireModTPCWW.YZScaleWidth: [ "TPC3_plane0_ratio_width", + "TPC3_plane1_ratio_width", + "TPC3_plane2_ratio_width"] +WireModTPCWW.WireLabel: "WMROI2Wire:PHYSCRATEDATATPCWW" +WireModTPCWW.ChannelLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" +WireModTPCWW.HitLabel: "dummyGaushit2dTPCWW" + +WireModSeq: [ WireModTPCEE, WireModTPCEW, WireModTPCWE, WireModTPCWW ] +WireModTab: { WireModTPCEE: @local::WireModTPCEE + WireModTPCEW: @local::WireModTPCEW + WireModTPCWE: @local::WireModTPCWE + WireModTPCWW: @local::WireModTPCWW } + +WMWireToROI: +{ + module_type: "WireToChannelROI" + WireModuleLabelVec: [ "WireModTPCEE", + "WireModTPCEW", + "WireModTPCWE", + "WireModTPCWW" ] + OutInstanceLabelVec: [ "PHYSCRATEDATATPCEE", + "PHYSCRATEDATATPCEW", + "PHYSCRATEDATATPCWE", + "PHYSCRATEDATATPCWW" ] + DiagnosticOutput: false +} + +WireModSequence: [ @sequence::dummyGaushit2dSeq, WMROI2Wire, @sequence::WireModSeq, WMWireToROI ] +WireModTable: +{ + @table::dummyGaushit2dTab + WMROI2Wire: @local::WMROI2Wire + @table::WireModTab + WMWireToROI: @local::WMWireToROI +} + +END_PROLOG + diff --git a/icaruscode/WireMod/WireModifierXXW_module.cc b/icaruscode/WireMod/WireModifierXXW_module.cc index 0205badf9..0175d0412 100644 --- a/icaruscode/WireMod/WireModifierXXW_module.cc +++ b/icaruscode/WireMod/WireModifierXXW_module.cc @@ -1,4 +1,5 @@ // std inlcudes +#include #include #include @@ -7,6 +8,7 @@ #include "TGraph.h" #include "TH1D.h" #include "TH2F.h" +#include "TTree.h" // art includes #include "art/Framework/Core/EDProducer.h" @@ -63,21 +65,70 @@ namespace wiremod art::InputTag fEDepShiftedLabel; // which are the EDeps? unsigned int fCryo; // which Cryo are we in? unsigned int fTPCset; // which TPC are we in? - bool fLocalRatios; // is the ratio file local? - bool fSaveHistsByChannel; // save modified signals by channel? - bool fSaveHistsByWire; // save modified signals by wire? - bool fInRads; // is the TGraph2D angle axis in radians? - bool fXAbs; // is the TGraph2D x an absolute value? + bool fLocalRatios; // is the ratio file local? + bool fSaveHistsByChannel; // save modified signals by channel? + bool fSaveHistsByWire; // save modified signals by wire? + bool fSaveROITree; // fill ROITree TTree (before/after waveforms)? + bool fSaveTickHistograms; // fill EDep/Hit tick-vs-channel TH2Fs? + bool fSavePerHitData; // write per-hit scaleQ/scaleSigma/truthX/thetaXW/isMC art products? + bool fInRads; // is the TGraph2D angle axis in radians? + bool fXAbs; // is the TGraph2D x an absolute value? + bool fAdditiveModification; // additive (true) vs multiplicative (false) ROI modification + bool fSetNullScaleIntegral; // if true, force r_Q=1 for all sub-ROIs (disable integral scaling) + bool fSetNullScaleWidth; // if true, force r_sigma=1 for all sub-ROIs (disable width scaling) + bool fUseChannelROIMode; // true = read ChannelROI directly (float precision), false = Wire (nominal) + art::InputTag fChannelLabel; // ChannelROI input label, used only when fUseChannelROIMode is true double fOffset; // ad hoc offset - TH2F* fEdepTickInd1; // What tick are the Edeps projected to (Induction 1)? - TH2F* fEdepTickInd2; // What tick are the Edeps projected to (Induction 2)? - TH2F* fEdepTickColl; // What tick are the Edeps projected to (Collection)? - TH2F* fHitTickInd1; // What tick are the hits at (Induction 1)? - TH2F* fHitTickInd2; // What tick are the hits at (Induction 2)? - TH2F* fHitTickColl; // What tick are the hits at (Collection)? - TH2F* fEdepXvsChan; - TH2F* fHitXvsChan; - TH2F* fHitBTvsDetProp; // BackTracker vs DetProp check hist + TH2F* fEdepTickInd1 = nullptr; + TH2F* fEdepTickInd2 = nullptr; + TH2F* fEdepTickColl = nullptr; + TH2F* fHitTickInd1 = nullptr; + TH2F* fHitTickInd2 = nullptr; + TH2F* fHitTickColl = nullptr; + TH2F* fEdepXvsChan = nullptr; + TH2F* fHitXvsChan = nullptr; + TH2F* fHitBTvsDetProp = nullptr; + + // diagnostic tree for ROI modification studies + TTree* fROITree = nullptr; + Int_t tb_run; + Int_t tb_subrun; + Int_t tb_event; + Int_t tb_channel; + Int_t tb_wire; + Int_t tb_tpc; + Int_t tb_plane; + Int_t tb_roi_idx; + Float_t tb_roi_begin; + Float_t tb_roi_before_total_q; + Float_t tb_roi_before_center; + Float_t tb_roi_before_sigma; + Float_t tb_roi_after_total_q; + Float_t tb_roi_after_center; + Float_t tb_roi_after_sigma; + std::vector tb_roi_before; + std::vector tb_roi_after; + std::vector tb_hit_before_center; + std::vector tb_hit_before_integral; + std::vector tb_hit_before_sigma; + std::vector tb_hit_after_center; + std::vector tb_hit_after_integral; + std::vector tb_hit_after_sigma; + std::vector tb_scale_Q; + std::vector tb_scale_sigma; + std::vector tb_theta_xw; + std::vector tb_truth_x; // truth X coordinate used for XXW graph lookup + std::vector tb_track_id; + // 0=no EDep matched (overlay data hit) + // 1=energy veto (truth_E<0.3 && Q>80) + // 2=scaling applied (scale!=1) + // 3=EDep matched but map evaluates to 1 + std::vector tb_scale_reason; + char tb_xxw_graph_name[128]; // name of the sigma TGraph2D used for this ROI's plane + std::vector fNameVec_sigma_XXW; // kept from reconfigure for tree output + Int_t tb_roi_n_hits; // total hits matched to this ROI + Int_t tb_roi_n_pulse_train_hits; // hits with GoodnessOfFit==-1 && DoF==1 + Int_t tb_roi_is_pulse_train; // 1 if any hit is a pulse train }; // end WireModifierXXW class @@ -152,13 +203,26 @@ namespace wiremod fTPCset = pset.get("TPCset"); // what, if anything, are we putting in the histogram files - fSaveHistsByChannel = pset.get("SaveByChannel", false); - fSaveHistsByWire = pset.get("SaveByWire" , false); + fSaveHistsByChannel = pset.get("SaveByChannel", false); + fSaveHistsByWire = pset.get("SaveByWire", false); + fSaveROITree = pset.get("SaveROITree", false); + fSaveTickHistograms = pset.get("SaveTickHistograms", false); + fSavePerHitData = pset.get("SavePerHitData", false); // do we need to reshape the TGraph2D? fInRads = pset.get("InRadians", true); fXAbs = pset.get("XAbs", false); + // additive vs multiplicative ROI modification + fAdditiveModification = pset.get("AdditiveModification", false); + fSetNullScaleIntegral = pset.get("SetNullScaleIntegral", false); + fSetNullScaleWidth = pset.get("SetNullScaleWidth", false); + + // ChannelROI direct mode (float precision, bypasses Wire intermediate) + fUseChannelROIMode = pset.get("UseChannelROIMode", false); + if (fUseChannelROIMode) + fChannelLabel = pset.get("ChannelLabel"); + // try to read in the graphs/splines from a file // if that file does not exist then fake them fRatioFileName_XXW = pset.get("RatioFileName_XXW", "NOFILE"); @@ -240,6 +304,7 @@ namespace wiremod } } std::vector nameVec_sigma_XXW = pset.get>("XXWScaleWidth"); + fNameVec_sigma_XXW = nameVec_sigma_XXW; for (auto const& name : nameVec_sigma_XXW) { mf::LogDebug("WireModifierXXW") @@ -353,47 +418,104 @@ namespace wiremod } } - // we make these things - produces>(); + // we make these things (mode-dependent) + if (!fUseChannelROIMode) + produces>(); + else + produces>(); + // per-hit scale info: one entry per hit in fHitLabel, sentinel -999 if no EDep match + produces>("scaleQ"); + produces>("scaleSigma"); + produces>("truthX"); + produces>("thetaXW"); + // per-hit MC flag: 1 = MC hit (BackTracker matched), 0 = overlay data hit + produces>("isMC"); + // per-hit truth direction and derived geometry (filled only when SavePerHitData: true) + produces>("dirX"); // truth dxdr + produces>("dirY"); // truth dydr + produces>("dirZ"); // truth dzdr + produces>("pitch"); // wirePitch / cosG [cm] + produces>("dQdX"); // integral_before / pitch [ADC/cm], no lifetime correction - // tick check art::ServiceHandle tfs; - fEdepTickInd1 = tfs->make("EdepTickInd1", "1st Induction;Projected Tick;Projected Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fEdepTickInd2 = tfs->make("EdepTickInd2", "2nd Induction;Projected Tick;Projected Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fEdepTickColl = tfs->make("EdepTickColl", "Collection;Projected Tick;Projected Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fHitTickInd1 = tfs->make("HitTickInd1", "1st Induction;Tick;Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fHitTickInd2 = tfs->make("HitTickInd2", "2nd Induction;Tick;Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fHitTickColl = tfs->make("HitTickColl", "Collection;Tick;Channel", - 4096, -128, 4096, - 4603, 0, 55296); - fEdepXvsChan = tfs->make("EdepXvsChan", ";X (cm);Channel", - 800, -400, 400, - 4603, 0, 55296); - fHitXvsChan = tfs->make("HitXvsChan", ";X (cm);Channel", - 800, -400, 400, - 4603, 0, 55296); - fHitBTvsDetProp = tfs->make("HitBTvsDetProp", - ";X (cm) From BackTracker;X (cm) From DetectorProperties", - 800, -400, 400, - 800, -400, 400); - fEdepTickInd1->SetMarkerColor(kRed); - fEdepTickInd2->SetMarkerColor(kRed); - fEdepTickColl->SetMarkerColor(kRed); - fEdepXvsChan->SetMarkerColor(kRed); - fHitTickInd1->SetMarkerColor(kBlack); - fHitTickInd2->SetMarkerColor(kBlack); - fHitTickColl->SetMarkerColor(kBlack); - fHitXvsChan->SetMarkerColor(kBlack); + + if (fSaveTickHistograms) + { + fEdepTickInd1 = tfs->make("EdepTickInd1", "1st Induction;Projected Tick;Projected Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fEdepTickInd2 = tfs->make("EdepTickInd2", "2nd Induction;Projected Tick;Projected Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fEdepTickColl = tfs->make("EdepTickColl", "Collection;Projected Tick;Projected Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fHitTickInd1 = tfs->make("HitTickInd1", "1st Induction;Tick;Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fHitTickInd2 = tfs->make("HitTickInd2", "2nd Induction;Tick;Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fHitTickColl = tfs->make("HitTickColl", "Collection;Tick;Channel", + 4096, -128, 4096, + 4603, 0, 55296); + fEdepXvsChan = tfs->make("EdepXvsChan", ";X (cm);Channel", + 800, -400, 400, + 4603, 0, 55296); + fHitXvsChan = tfs->make("HitXvsChan", ";X (cm);Channel", + 800, -400, 400, + 4603, 0, 55296); + fHitBTvsDetProp = tfs->make("HitBTvsDetProp", + ";X (cm) From BackTracker;X (cm) From DetectorProperties", + 800, -400, 400, + 800, -400, 400); + fEdepTickInd1->SetMarkerColor(kRed); + fEdepTickInd2->SetMarkerColor(kRed); + fEdepTickColl->SetMarkerColor(kRed); + fEdepXvsChan->SetMarkerColor(kRed); + fHitTickInd1->SetMarkerColor(kBlack); + fHitTickInd2->SetMarkerColor(kBlack); + fHitTickColl->SetMarkerColor(kBlack); + fHitXvsChan->SetMarkerColor(kBlack); + } + + fROITree = tfs->make("ROITree", "ROI modification diagnostic"); + if (fSaveROITree) + { + fROITree->Branch("run", &tb_run, "run/I"); + fROITree->Branch("subrun", &tb_subrun, "subrun/I"); + fROITree->Branch("event", &tb_event, "event/I"); + fROITree->Branch("channel", &tb_channel, "channel/I"); + fROITree->Branch("wire", &tb_wire, "wire/I"); + fROITree->Branch("tpc", &tb_tpc, "tpc/I"); + fROITree->Branch("plane", &tb_plane, "plane/I"); + fROITree->Branch("roi_idx", &tb_roi_idx, "roi_idx/I"); + fROITree->Branch("roi_begin", &tb_roi_begin, "roi_begin/F"); + fROITree->Branch("roi_before_total_q", &tb_roi_before_total_q, "roi_before_total_q/F"); + fROITree->Branch("roi_before_center", &tb_roi_before_center, "roi_before_center/F"); + fROITree->Branch("roi_before_sigma", &tb_roi_before_sigma, "roi_before_sigma/F"); + fROITree->Branch("roi_after_total_q", &tb_roi_after_total_q, "roi_after_total_q/F"); + fROITree->Branch("roi_after_center", &tb_roi_after_center, "roi_after_center/F"); + fROITree->Branch("roi_after_sigma", &tb_roi_after_sigma, "roi_after_sigma/F"); + fROITree->Branch("roi_before", &tb_roi_before); + fROITree->Branch("roi_after", &tb_roi_after); + fROITree->Branch("hit_before_center", &tb_hit_before_center); + fROITree->Branch("hit_before_integral", &tb_hit_before_integral); + fROITree->Branch("hit_before_sigma", &tb_hit_before_sigma); + fROITree->Branch("hit_after_center", &tb_hit_after_center); + fROITree->Branch("hit_after_integral", &tb_hit_after_integral); + fROITree->Branch("hit_after_sigma", &tb_hit_after_sigma); + fROITree->Branch("scale_Q", &tb_scale_Q); + fROITree->Branch("scale_sigma", &tb_scale_sigma); + fROITree->Branch("theta_xw", &tb_theta_xw); + fROITree->Branch("truth_x", &tb_truth_x); + fROITree->Branch("track_id", &tb_track_id); + fROITree->Branch("scale_reason", &tb_scale_reason); + fROITree->Branch("xxw_graph_name", tb_xxw_graph_name, "xxw_graph_name/C"); + fROITree->Branch("roi_n_hits", &tb_roi_n_hits, "roi_n_hits/I"); + fROITree->Branch("roi_n_pulse_train_hits", &tb_roi_n_pulse_train_hits, "roi_n_pulse_train_hits/I"); + fROITree->Branch("roi_is_pulse_train", &tb_roi_is_pulse_train, "roi_is_pulse_train/I"); + } } //------------------------------------------------ @@ -406,12 +528,13 @@ namespace wiremod const detinfo::DetectorPropertiesData detProp = art::ServiceHandle()->DataFor(evt); - // get the things to do the things on - art::Handle< std::vector > wireHandle; - evt.getByLabel(fWireLabel, wireHandle); - auto const& wireVec(*wireHandle); - - art::FindManyP digit_assn(wireHandle, evt, fWireLabel); + // handles, read conditionally based on mode + art::Handle< std::vector > wireHandle; + art::Handle< std::vector > chanROIHandle; + if (!fUseChannelROIMode) + evt.getByLabel(fWireLabel, wireHandle); + else + evt.getByLabel(fChannelLabel, chanROIHandle); art::Handle< std::vector > edepHandle; evt.getByLabel(fEDepLabel, edepHandle); @@ -425,15 +548,27 @@ namespace wiremod evt.getByLabel(fHitLabel, hitHandle); auto const& hitVec(*hitHandle); - // put the new stuff somewhere - std::unique_ptr> new_wires(new std::vector()); + // per-hit scale products, only sized when SavePerHitData is on + const size_t perHitN = fSavePerHitData ? hitVec.size() : 0; + auto vec_scale_q = std::make_unique>(perHitN, -999.f); + auto vec_scale_sigma = std::make_unique>(perHitN, -999.f); + auto vec_truth_x = std::make_unique>(perHitN, -999.f); + auto vec_theta_xw = std::make_unique>(perHitN, -999.f); + auto vec_is_mc = std::make_unique> (perHitN, 0); + auto vec_dirx = std::make_unique>(perHitN, -999.f); + auto vec_diry = std::make_unique>(perHitN, -999.f); + auto vec_dirz = std::make_unique>(perHitN, -999.f); + auto vec_pitch = std::make_unique>(perHitN, -999.f); + auto vec_dqdx = std::make_unique>(perHitN, -999.f); - // Get the tick offset by comparting the hit ticks with the tick gotten by - // projecting the backtracked hit X to a tick. Yes this is silly thanks for asking + // Get the tick offset by comparing the hit ticks with the tick gotten by + // projecting the backtracked hit X to a tick. Yes this is silly thanks for asking. + // Also fills vec_is_mc: 1 if BackTracker finds a match (MC hit), 0 otherwise (overlay). double BT_Offset = 0; unsigned int BT_Hits = 0; - for (auto const& hit : hitVec) + for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) { + auto const& hit = hitVec[i_h]; double hitX = std::numeric_limits::max(); try { @@ -441,16 +576,16 @@ namespace wiremod } catch (...) { - // Overlay hit, ignore + // Overlay hit continue; } - double hitTick = hit.PeakTime(); - double projTick = - detProp.ConvertXToTicks(hitX, hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat); + if (fSavePerHitData) (*vec_is_mc)[i_h] = 1; + double hitTick = hit.PeakTime(); + double projTick = detProp.ConvertXToTicks(hitX, hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat); BT_Offset += (hitTick - projTick); ++BT_Hits; } - BT_Offset /= BT_Hits; + if (BT_Hits > 0) BT_Offset /= BT_Hits; sys::WireModUtility wmUtil(fGeometry, fWireReadout, detProp, false, // Channel Scale @@ -472,19 +607,24 @@ namespace wiremod wmUtil.graph2Ds_Sigma_XXW = fGraph_sigma_XXW; } + wmUtil.additiveModification = fAdditiveModification; + // add some debugging here mf::LogVerbatim("WireModifierXXW") << "DUMP CONFIG:" << '\n' << "---------------------------------------------------" << '\n' - << " applyChannelScale: " << wmUtil.applyChannelScale << '\n' - << " applyXScale: " << wmUtil.applyXScale << '\n' - << " applyYZScale: " << wmUtil.applyYZScale << '\n' - << " applyXXWScale: " << wmUtil.applyXXWScale << '\n' - << " applyXZAngleScale: " << wmUtil.applyXZAngleScale << '\n' - << " applyYZAngleScale: " << wmUtil.applyYZAngleScale << '\n' - << " applydEdXScale: " << wmUtil.applydEdXScale << '\n' - << " readoutWindowTicks: " << wmUtil.readoutWindowTicks << '\n' - << " tickOffset: " << wmUtil.tickOffset << '\n' + << " applyChannelScale: " << wmUtil.applyChannelScale << '\n' + << " applyXScale: " << wmUtil.applyXScale << '\n' + << " applyYZScale: " << wmUtil.applyYZScale << '\n' + << " applyXXWScale: " << wmUtil.applyXXWScale << '\n' + << " applyXZAngleScale: " << wmUtil.applyXZAngleScale << '\n' + << " applyYZAngleScale: " << wmUtil.applyYZAngleScale << '\n' + << " applydEdXScale: " << wmUtil.applydEdXScale << '\n' + << " additiveModification: " << wmUtil.additiveModification << '\n' + << " setNullScaleIntegral: " << fSetNullScaleIntegral << '\n' + << " setNullScaleWidth: " << fSetNullScaleWidth << '\n' + << " readoutWindowTicks: " << wmUtil.readoutWindowTicks << '\n' + << " tickOffset: " << wmUtil.tickOffset << '\n' << "---------------------------------------------------"; // tick check @@ -492,61 +632,71 @@ namespace wiremod << "edepVec.size() = " << edepVec.size() << '\n' << "edepShiftedVec.size() = " << edepShiftedVec.size() << '\n' << "hitVec.size() = " << hitVec.size(); - for (auto const& edep : edepVec) + if (fSaveTickHistograms) { - if (fWireReadout->FindTPCsetAtPosition(edep.MidPoint()) != readout::TPCsetID(fCryo, fTPCset)) - continue; - - geo::TPCGeo const* curTPCGeomPtr = fGeometry->PositionToTPCptr(edep.MidPoint()); - for (auto const& plane : fWireReadout->Iterate(curTPCGeomPtr->ID())) + for (auto const& edep : edepVec) { - TH2F* targetHist = (plane.View() == geo::kY) ? fEdepTickInd1 - : (plane.View() == geo::kV) ? fEdepTickInd2 - : (plane.View() == geo::kU) ? fEdepTickColl - : nullptr; - if (targetHist == nullptr) continue; + if (fWireReadout->FindTPCsetAtPosition(edep.MidPoint()) != readout::TPCsetID(fCryo, fTPCset)) + continue; - geo::WireID wireID = plane.NearestWireID(edep.MidPoint()); - float projTick = detProp.ConvertXToTicks(edep.X(), wireID.Plane, wireID.TPC, wireID.Cryostat) + wmUtil.tickOffset; - float projChan = fWireReadout->PlaneWireToChannel(wireID); - targetHist->Fill(projTick, projChan); - fEdepXvsChan->Fill(edep.X(), projChan); + geo::TPCGeo const* curTPCGeomPtr = fGeometry->PositionToTPCptr(edep.MidPoint()); + for (auto const& plane : fWireReadout->Iterate(curTPCGeomPtr->ID())) + { + TH2F* targetHist = (plane.View() == geo::kY) ? fEdepTickInd1 + : (plane.View() == geo::kV) ? fEdepTickInd2 + : (plane.View() == geo::kU) ? fEdepTickColl + : nullptr; + if (targetHist == nullptr) continue; + + geo::WireID wireID; + try { wireID = plane.NearestWireID(edep.MidPoint()); } + catch (...) { continue; } + float projTick = detProp.ConvertXToTicks(edep.X(), wireID.Plane, wireID.TPC, wireID.Cryostat) + wmUtil.tickOffset; + float projChan = fWireReadout->PlaneWireToChannel(wireID); + targetHist->Fill(projTick, projChan); + fEdepXvsChan->Fill(edep.X(), projChan); + } } - } - for (auto const& hit : hitVec) - { - TH2F* targetHist = (hit.View() == geo::kY) ? fHitTickInd1 - : (hit.View() == geo::kV) ? fHitTickInd2 - : (hit.View() == geo::kU) ? fHitTickColl - : nullptr; - if (targetHist == nullptr) continue; - - float tick = hit.PeakTime(); - float chan = hit.Channel(); - targetHist->Fill(tick, chan); - float hitX = detProp.ConvertTicksToX(tick, - hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat); - try + for (auto const& hit : hitVec) { - float hitXBT = fBT->HitToXYZ(fDetClocksData, hit).at(0); - fHitBTvsDetProp->Fill(hitXBT, hitX); - hitX = hitXBT; + TH2F* targetHist = (hit.View() == geo::kY) ? fHitTickInd1 + : (hit.View() == geo::kV) ? fHitTickInd2 + : (hit.View() == geo::kU) ? fHitTickColl + : nullptr; + if (targetHist == nullptr) continue; + + float tick = hit.PeakTime(); + float chan = hit.Channel(); + targetHist->Fill(tick, chan); + float hitX = detProp.ConvertTicksToX(tick, + hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat); + try + { + float hitXBT = fBT->HitToXYZ(fDetClocksData, hit).at(0); + fHitBTvsDetProp->Fill(hitXBT, hitX); + hitX = hitXBT; + } + catch (...) { } + fHitXvsChan->Fill(hitX, chan); } - catch (...) { } - fHitXvsChan->Fill(hitX, chan); } - // do the things - double offset_ADC = 0; // don't use an offset atm - mf::LogVerbatim("WireModifierXXW") - << "Get Edep Map"; + double offset_ADC = 0; + + if (!fUseChannelROIMode) { + // Wire (nominal) mode + auto const& wireVec(*wireHandle); + std::unique_ptr> new_wires(new std::vector()); + + mf::LogVerbatim("WireModifierXXW") << "Get Edep Map"; wmUtil.FillROIMatchedEdepMap(edepShiftedVec, wireVec, offset_ADC); - mf::LogVerbatim("WireModifierXXW") - << "Got Edep Map." << '\n' - << "Get Hit Map"; + mf::LogVerbatim("WireModifierXXW") << "Got Edep Map." << '\n' << "Get Hit Map"; wmUtil.FillROIMatchedHitMap(hitVec, wireVec); - mf::LogVerbatim("WireModifierXXW") - << "Got Hit Map."; + mf::LogVerbatim("WireModifierXXW") << "Got Hit Map."; + + tb_run = static_cast(evt.run()); + tb_subrun = static_cast(evt.subRun()); + tb_event = static_cast(evt.event()); // loop-de-loop for(size_t i_w = 0; i_w < wireVec.size(); ++i_w) @@ -591,6 +741,7 @@ namespace wiremod << " Checking ROI " << i_r; auto const& range = wire.SignalROI().get_ranges()[i_r]; + sys::WireModUtility::ROI_Key_t roi_key(wire.Channel(), i_r); std::vector modified_data(range.data()); @@ -669,20 +820,26 @@ namespace wiremod std::map SubROIMatchedScalesMap; + std::map SubROIMatchedTruthMap; + std::map SubROIScaleReasonMap; for (auto const& subroi_prop : subROIPropVec) { sys::WireModUtility::ScaleValues_t scale_vals; auto key = subroi_prop.key; auto key_it = SubROIMatchedEdepMap.find(key); + int scale_reason = 0; if (key_it != SubROIMatchedEdepMap.end() && key_it->second.size() > 0) { auto truth_vals = wmUtil.CalcPropertiesFromEdeps(key_it->second, offset_ADC); + SubROIMatchedTruthMap[key] = truth_vals; if (truth_vals.total_energy < 0.3 && subroi_prop.total_q > 80) { scale_vals.r_Q = 1.; scale_vals.r_sigma = 1.; + scale_reason = 1; // energy veto } else { scale_vals = wmUtil.GetScaleValues(truth_vals, roi_properties); @@ -692,18 +849,203 @@ namespace wiremod << "Scaling! Q scale: " << scale_vals.r_Q << " sigma scale: " << scale_vals.r_sigma; isModified = true; + scale_reason = 2; // scaling applied + } else { + scale_reason = 3; // EDep matched but map evaluates to 1 } } } else { scale_vals.r_Q = 1.; scale_vals.r_sigma = 1.; + scale_reason = 0; // no EDep matched (likely overlay data hit) } - SubROIMatchedScalesMap[key] = scale_vals; + SubROIMatchedScalesMap[key] = scale_vals; + SubROIScaleReasonMap[key] = scale_reason; } - wmUtil.ModifyROI(modified_data, roi_properties, subROIPropVec, SubROIMatchedScalesMap); + // null-scale overrides: force r_Q and/or r_sigma to 1 for all sub-ROIs. + // applied in-place so the diagnostic tree records what was actually used. + if (fSetNullScaleIntegral || fSetNullScaleWidth) { + for (auto& kv : SubROIMatchedScalesMap) { + if (fSetNullScaleIntegral) kv.second.r_Q = 1.0; + if (fSetNullScaleWidth) kv.second.r_sigma = 1.0; + } + } + + wmUtil.ModifyROI(modified_data, roi_properties, subROIPropVec, SubROIMatchedScalesMap, 1.5); + + // fill diagnostic tree + if (fSaveROITree) + { + // wire/geometry info + std::vector wireIDs = fWireReadout->ChannelToWire(wire.Channel()); + geo::WireID firstWireID = wireIDs.empty() ? geo::WireID{} : wireIDs.front(); + tb_channel = static_cast(wire.Channel()); + tb_wire = static_cast(firstWireID.Wire); + tb_tpc = static_cast(firstWireID.TPC); + tb_plane = static_cast(my_plane); + tb_roi_idx = static_cast(i_r); + + // ROI before + tb_roi_begin = roi_properties.begin; + tb_roi_before_total_q = roi_properties.total_q; + tb_roi_before_center = roi_properties.center; + tb_roi_before_sigma = roi_properties.sigma; + tb_roi_before.assign(range.data().begin(), range.data().end()); + + // ROI after: compute charge-weighted properties from modified waveform + float mod_total_q = 0.f, mod_center = 0.f, mod_sigma = 0.f; + for (size_t i_t = 0; i_t < modified_data.size(); ++i_t) + { + mod_total_q += modified_data[i_t]; + mod_center += modified_data[i_t] * (i_t + roi_properties.begin); + } + if (mod_total_q > 0.f) mod_center /= mod_total_q; + for (size_t i_t = 0; i_t < modified_data.size(); ++i_t) + { + float t = static_cast(i_t) + roi_properties.begin; + mod_sigma += modified_data[i_t] * (t - mod_center) * (t - mod_center); + } + if (mod_total_q > 0.f) mod_sigma = std::sqrt(mod_sigma / mod_total_q); + tb_roi_after_total_q = mod_total_q; + tb_roi_after_center = mod_center; + tb_roi_after_sigma = mod_sigma; + tb_roi_after = modified_data; + + // pulse train detection (GoodnessOfFit==-1 && DoF==1 is gaushit's sentinel) + { + int n_pt = 0; + for (auto const& hptr : matchedHitPtrVec) + if (hptr->GoodnessOfFit() == -1.0f && hptr->DegreesOfFreedom() == 1) + ++n_pt; + tb_roi_n_hits = static_cast(matchedHitPtrVec.size()); + tb_roi_n_pulse_train_hits = static_cast(n_pt); + tb_roi_is_pulse_train = (n_pt > 0) ? 1 : 0; + } + + // per-subROI hit candidate info and scale factors + tb_hit_before_center.clear(); + tb_hit_before_integral.clear(); + tb_hit_before_sigma.clear(); + tb_hit_after_center.clear(); + tb_hit_after_integral.clear(); + tb_hit_after_sigma.clear(); + tb_scale_Q.clear(); + tb_scale_sigma.clear(); + tb_theta_xw.clear(); + tb_truth_x.clear(); + tb_track_id.clear(); + tb_scale_reason.clear(); + { + std::string gname = (tb_plane < (int)fNameVec_sigma_XXW.size()) + ? fNameVec_sigma_XXW[tb_plane] : "unknown"; + std::strncpy(tb_xxw_graph_name, gname.c_str(), 127); + tb_xxw_graph_name[127] = '\0'; + } + + for (auto const& subroi_prop : subROIPropVec) + { + auto key = subroi_prop.key; + auto scale_vals = SubROIMatchedScalesMap.at(key); + + tb_hit_before_center.push_back(subroi_prop.center); + tb_hit_before_integral.push_back(subroi_prop.total_q); + tb_hit_before_sigma.push_back(subroi_prop.sigma); + tb_hit_after_center.push_back(subroi_prop.center); + tb_hit_after_integral.push_back(subroi_prop.total_q * static_cast(scale_vals.r_Q)); + tb_hit_after_sigma.push_back(subroi_prop.sigma * static_cast(scale_vals.r_sigma)); + tb_scale_Q.push_back(static_cast(scale_vals.r_Q)); + tb_scale_sigma.push_back(static_cast(scale_vals.r_sigma)); + + // dominant track ID: pick the TrackID with highest summed energy + int dominant_tid = -999; + auto edep_it = SubROIMatchedEdepMap.find(key); + if (edep_it != SubROIMatchedEdepMap.end() && !edep_it->second.empty()) + { + std::map tidEnergy; + for (auto const& ep : edep_it->second) + tidEnergy[ep->TrackID()] += ep->E(); + double maxE = 0.; + for (auto const& te : tidEnergy) + if (te.second > maxE) { maxE = te.second; dominant_tid = te.first; } + } + tb_track_id.push_back(dominant_tid); + tb_scale_reason.push_back(SubROIScaleReasonMap.at(key)); + + // angle in degrees (ThetaXW returns radians) + truth X for graph lookup + float theta_deg = -999.f; + float truth_x_val = -9999.f; + auto truth_it = SubROIMatchedTruthMap.find(key); + if (truth_it != SubROIMatchedTruthMap.end()) + { + auto const& tv = truth_it->second; + truth_x_val = static_cast(tv.x); + geo::TPCGeo const* tpcPtr = fGeometry->PositionToTPCptr( + geo::Point_t{(double)tv.x, (double)tv.y, (double)tv.z}); + if (tpcPtr != nullptr) + { + try + { + auto const& plane_obj = fWireReadout->Plane(tpcPtr->ID(), wire.View()); + theta_deg = static_cast(wmUtil.ThetaXW( + tv.dxdr, tv.dydr, tv.dzdr, plane_obj.ThetaZ()) * TMath::RadToDeg()); + } + catch (...) {} + } + } + tb_theta_xw.push_back(theta_deg); + tb_truth_x.push_back(truth_x_val); + } + + fROITree->Fill(); + } + + // Fill per-hit scale product vectors indexed by position in hitVec + if (fSavePerHitData && it_hit_map != wmUtil.ROIMatchedHitMap.end()) + { + for (auto const& subroi_prop : subROIPropVec) + { + size_t i_subroi = subroi_prop.key.second; + if (i_subroi >= it_hit_map->second.size()) continue; + size_t hitVec_idx = it_hit_map->second[i_subroi]; + + auto const& scale_vals = SubROIMatchedScalesMap.at(subroi_prop.key); + (*vec_scale_q )[hitVec_idx] = static_cast(scale_vals.r_Q); + (*vec_scale_sigma)[hitVec_idx] = static_cast(scale_vals.r_sigma); + + auto truth_it = SubROIMatchedTruthMap.find(subroi_prop.key); + if (truth_it != SubROIMatchedTruthMap.end()) + { + auto const& tv = truth_it->second; + (*vec_truth_x)[hitVec_idx] = static_cast(tv.x); + geo::TPCGeo const* tpcPtr = fGeometry->PositionToTPCptr( + geo::Point_t{tv.x, tv.y, tv.z}); + if (tpcPtr) + { + try { + auto const& plane_obj = fWireReadout->Plane( + hitVec[hitVec_idx].WireID().planeID()); + double thetaZ = plane_obj.ThetaZ(); + (*vec_theta_xw)[hitVec_idx] = static_cast( + wmUtil.ThetaXW(tv.dxdr, tv.dydr, tv.dzdr, thetaZ) * TMath::RadToDeg()); + (*vec_dirx)[hitVec_idx] = static_cast(tv.dxdr); + (*vec_diry)[hitVec_idx] = static_cast(tv.dydr); + (*vec_dirz)[hitVec_idx] = static_cast(tv.dzdr); + double cosG = std::abs(tv.dydr * std::sin(thetaZ) + tv.dzdr * std::cos(thetaZ)); + float wp = static_cast(plane_obj.WirePitch()); + float pval = (cosG > 1e-6) ? wp / static_cast(cosG) : -999.f; + (*vec_pitch)[hitVec_idx] = pval; + (*vec_dqdx)[hitVec_idx] = (pval > 0.f) + ? hitVec[hitVec_idx].Integral() / pval : -999.f; + } + catch (...) {} + } + } + } + } + new_rois.add_range(roi_properties.begin, modified_data); } @@ -773,6 +1115,434 @@ namespace wiremod } // end loop over wires evt.put(std::move(new_wires)); + + } else { + // ChannelROI (float precision) mode + auto const& chanROIVec(*chanROIHandle); + std::unique_ptr> new_chanROIs(new std::vector()); + + mf::LogVerbatim("WireModifierXXW") << "Get Edep Map"; + wmUtil.FillROIMatchedEdepMap(edepShiftedVec, chanROIVec, offset_ADC); + mf::LogVerbatim("WireModifierXXW") << "Got Edep Map." << '\n' << "Get Hit Map"; + wmUtil.FillROIMatchedHitMap(hitVec, chanROIVec); + mf::LogVerbatim("WireModifierXXW") << "Got Hit Map."; + + tb_run = static_cast(evt.run()); + tb_subrun = static_cast(evt.subRun()); + tb_event = static_cast(evt.event()); + + for(size_t i_c = 0; i_c < chanROIVec.size(); ++i_c) + { + mf::LogDebug("WireModifierXXW") << "Checking chanROI " << i_c; + + auto const& chanROI = chanROIVec.at(i_c); + recob::ChannelROI::RegionsOfInterest_t new_rois; + bool isModified = false; + + if (chanROI.NSignal() == 0) + continue; + + new_rois.resize(chanROI.NSignal()); + + geo::View_t chanView = fWireReadout->View(chanROI.Channel()); + unsigned int my_plane = geo::kUnknown; + if (chanView == fWireReadout->Plane(geo::PlaneID(0, 0, 0)).View()) my_plane = 0; + else if (chanView == fWireReadout->Plane(geo::PlaneID(0, 0, 1)).View()) my_plane = 1; + else if (chanView == fWireReadout->Plane(geo::PlaneID(0, 0, 2)).View()) my_plane = 2; + + if (my_plane == geo::kUnknown) + mf::LogDebug("WireModifierXXW") << "ChannelROI on unsupported plane. Skip."; + + // SignalROIF() returns by value, so a const& to its sub-objects would dangle. + // SignalROI() is a stable const ref instead. + auto const& roiShort = chanROI.SignalROI(); + const float adcScale = 1.0f / static_cast(chanROI.ADCScaleFactor()); + + for(size_t i_r = 0; i_r < roiShort.get_ranges().size(); ++i_r) + { + mf::LogDebug("WireModifierXXW") << " Checking ROI " << i_r; + + auto const& range_short = roiShort.get_ranges()[i_r]; + sys::WireModUtility::ROI_Key_t roi_key(chanROI.Channel(), i_r); + + // convert short -> float (same as SignalROIF() but without the temporary) + std::vector modified_data; + modified_data.reserve(range_short.data().size()); + for (auto s : range_short.data()) + modified_data.push_back(static_cast(s) * adcScale); + + auto toShort = [&chanROI](float v) { + return static_cast(std::round(v * chanROI.ADCScaleFactor())); + }; + + auto it_map = wmUtil.ROIMatchedEdepMap.find(roi_key); + if(it_map==wmUtil.ROIMatchedEdepMap.end()){ + std::vector pass_short(modified_data.size()); + std::transform(modified_data.begin(), modified_data.end(), pass_short.begin(), toShort); + new_rois.add_range(range_short.begin_index(), pass_short); + mf::LogDebug("WireModifierXXW") << " Could not find matching Edep. Skip"; + continue; + } + std::vector matchedEdepIdxVec = it_map->second; + if(matchedEdepIdxVec.size() == 0) + { + std::vector pass_short(modified_data.size()); + std::transform(modified_data.begin(), modified_data.end(), pass_short.begin(), toShort); + new_rois.add_range(range_short.begin_index(), pass_short); + mf::LogDebug("WireModifierXXW") << " No indices for Edep. Skip"; + continue; + } + std::vector matchedEdepPtrVec; + std::vector matchedEdepShiftedPtrVec; + for(auto i_e : matchedEdepIdxVec) + { + matchedEdepPtrVec.push_back(&edepVec[i_e]); + matchedEdepShiftedPtrVec.push_back(&edepShiftedVec[i_e]); + } + mf::LogDebug("WireModifierXXW") << " Found " << matchedEdepPtrVec.size() << " Edeps"; + + std::vector matchedHitPtrVec; + auto it_hit_map = wmUtil.ROIMatchedHitMap.find(roi_key); + if( it_hit_map != wmUtil.ROIMatchedHitMap.end() ) { + for( auto i_h : it_hit_map->second ) + matchedHitPtrVec.push_back(&hitVec[i_h]); + } + mf::LogDebug("WireModifierXXW") << " Found " << matchedHitPtrVec.size() << " matching hits"; + + auto roi_properties = wmUtil.CalcROIProperties(chanROI, i_r); + mf::LogDebug("WireModifierXXW") + << " ROI Properties:" << '\n' + << " key: (" << roi_properties.key.first + << ", " << roi_properties.key.second << ")" << '\n' + << " view: " << roi_properties.view << '\n' + << " begin: " << roi_properties.begin << '\n' + << " end: " << roi_properties.end << '\n' + << " total_q: " << roi_properties.total_q << '\n' + << " center: " << roi_properties.center << '\n' + << " sigma: " << roi_properties.sigma; + + auto subROIPropVec = wmUtil.CalcSubROIProperties(roi_properties, matchedHitPtrVec); + mf::LogDebug("WireModifierXXW") << " have " << subROIPropVec.size() << " SubROI"; + + auto SubROIMatchedEdepShiftedMap = + wmUtil.MatchEdepsToSubROIs(subROIPropVec, matchedEdepShiftedPtrVec, offset_ADC); + std::map> + SubROIMatchedEdepMap; + for (auto const& key_edepPtrVec_pair : SubROIMatchedEdepShiftedMap) + { + auto key = key_edepPtrVec_pair.first; + for (auto const& shifted_edep_ptr : key_edepPtrVec_pair.second) + { + for (size_t i_e = 0; i_e < matchedEdepShiftedPtrVec.size(); ++i_e) + { + if (shifted_edep_ptr == matchedEdepShiftedPtrVec[i_e]) + { + SubROIMatchedEdepMap[key].push_back(matchedEdepPtrVec[i_e]); + break; + } + } + } + } + + std::map SubROIMatchedScalesMap; + std::map SubROIMatchedTruthMap; + std::map SubROIScaleReasonMap; + for (auto const& subroi_prop : subROIPropVec) + { + sys::WireModUtility::ScaleValues_t scale_vals; + auto key = subroi_prop.key; + auto key_it = SubROIMatchedEdepMap.find(key); + int scale_reason = 0; + + if (key_it != SubROIMatchedEdepMap.end() && key_it->second.size() > 0) + { + auto truth_vals = wmUtil.CalcPropertiesFromEdeps(key_it->second, offset_ADC); + SubROIMatchedTruthMap[key] = truth_vals; + + if (truth_vals.total_energy < 0.3 && subroi_prop.total_q > 80) + { + scale_vals.r_Q = 1.; + scale_vals.r_sigma = 1.; + scale_reason = 1; + } else + { + scale_vals = wmUtil.GetScaleValues(truth_vals, roi_properties); + if (scale_vals.r_Q != 1. || scale_vals.r_sigma != 1.) + { + mf::LogDebug("WireModifierXXW") + << "Scaling! Q scale: " << scale_vals.r_Q + << " sigma scale: " << scale_vals.r_sigma; + isModified = true; + scale_reason = 2; + } else { + scale_reason = 3; + } + } + } else + { + scale_vals.r_Q = 1.; + scale_vals.r_sigma = 1.; + scale_reason = 0; + } + + SubROIMatchedScalesMap[key] = scale_vals; + SubROIScaleReasonMap[key] = scale_reason; + } + + if (fSetNullScaleIntegral || fSetNullScaleWidth) { + for (auto& kv : SubROIMatchedScalesMap) { + if (fSetNullScaleIntegral) kv.second.r_Q = 1.0; + if (fSetNullScaleWidth) kv.second.r_sigma = 1.0; + } + } + + // snapshot of float data before in-place modification (used by ROITree "before" branch) + std::vector modified_data_before = fSaveROITree ? modified_data : std::vector{}; + + wmUtil.ModifyROI(modified_data, roi_properties, subROIPropVec, SubROIMatchedScalesMap, 1.5); + + if (fSaveROITree) + { + std::vector wireIDs = fWireReadout->ChannelToWire(chanROI.Channel()); + geo::WireID firstWireID = wireIDs.empty() ? geo::WireID{} : wireIDs.front(); + tb_channel = static_cast(chanROI.Channel()); + tb_wire = static_cast(firstWireID.Wire); + tb_tpc = static_cast(firstWireID.TPC); + tb_plane = static_cast(my_plane); + tb_roi_idx = static_cast(i_r); + + tb_roi_begin = roi_properties.begin; + tb_roi_before_total_q = roi_properties.total_q; + tb_roi_before_center = roi_properties.center; + tb_roi_before_sigma = roi_properties.sigma; + // modified_data holds the float-converted original (before ModifyROI runs) + tb_roi_before = modified_data_before; + + float mod_total_q = 0.f, mod_center = 0.f, mod_sigma = 0.f; + for (size_t i_t = 0; i_t < modified_data.size(); ++i_t) + { + mod_total_q += modified_data[i_t]; + mod_center += modified_data[i_t] * (i_t + roi_properties.begin); + } + if (mod_total_q > 0.f) mod_center /= mod_total_q; + for (size_t i_t = 0; i_t < modified_data.size(); ++i_t) + { + float t = static_cast(i_t) + roi_properties.begin; + mod_sigma += modified_data[i_t] * (t - mod_center) * (t - mod_center); + } + if (mod_total_q > 0.f) mod_sigma = std::sqrt(mod_sigma / mod_total_q); + tb_roi_after_total_q = mod_total_q; + tb_roi_after_center = mod_center; + tb_roi_after_sigma = mod_sigma; + tb_roi_after = modified_data; + + { + int n_pt = 0; + for (auto const& hptr : matchedHitPtrVec) + if (hptr->GoodnessOfFit() == -1.0f && hptr->DegreesOfFreedom() == 1) + ++n_pt; + tb_roi_n_hits = static_cast(matchedHitPtrVec.size()); + tb_roi_n_pulse_train_hits = static_cast(n_pt); + tb_roi_is_pulse_train = (n_pt > 0) ? 1 : 0; + } + + tb_hit_before_center.clear(); + tb_hit_before_integral.clear(); + tb_hit_before_sigma.clear(); + tb_hit_after_center.clear(); + tb_hit_after_integral.clear(); + tb_hit_after_sigma.clear(); + tb_scale_Q.clear(); + tb_scale_sigma.clear(); + tb_theta_xw.clear(); + tb_truth_x.clear(); + tb_track_id.clear(); + tb_scale_reason.clear(); + { + std::string gname = (tb_plane < (int)fNameVec_sigma_XXW.size()) + ? fNameVec_sigma_XXW[tb_plane] : "unknown"; + std::strncpy(tb_xxw_graph_name, gname.c_str(), 127); + tb_xxw_graph_name[127] = '\0'; + } + + for (auto const& subroi_prop : subROIPropVec) + { + auto key = subroi_prop.key; + auto scale_vals = SubROIMatchedScalesMap.at(key); + + tb_hit_before_center.push_back(subroi_prop.center); + tb_hit_before_integral.push_back(subroi_prop.total_q); + tb_hit_before_sigma.push_back(subroi_prop.sigma); + tb_hit_after_center.push_back(subroi_prop.center); + tb_hit_after_integral.push_back(subroi_prop.total_q * static_cast(scale_vals.r_Q)); + tb_hit_after_sigma.push_back(subroi_prop.sigma * static_cast(scale_vals.r_sigma)); + tb_scale_Q.push_back(static_cast(scale_vals.r_Q)); + tb_scale_sigma.push_back(static_cast(scale_vals.r_sigma)); + + int dominant_tid = -999; + auto edep_it = SubROIMatchedEdepMap.find(key); + if (edep_it != SubROIMatchedEdepMap.end() && !edep_it->second.empty()) + { + std::map tidEnergy; + for (auto const& ep : edep_it->second) + tidEnergy[ep->TrackID()] += ep->E(); + double maxE = 0.; + for (auto const& te : tidEnergy) + if (te.second > maxE) { maxE = te.second; dominant_tid = te.first; } + } + tb_track_id.push_back(dominant_tid); + tb_scale_reason.push_back(SubROIScaleReasonMap.at(key)); + + float theta_deg = -999.f; + float truth_x_val = -9999.f; + auto truth_it = SubROIMatchedTruthMap.find(key); + if (truth_it != SubROIMatchedTruthMap.end()) + { + auto const& tv = truth_it->second; + truth_x_val = static_cast(tv.x); + geo::TPCGeo const* tpcPtr = fGeometry->PositionToTPCptr( + geo::Point_t{(double)tv.x, (double)tv.y, (double)tv.z}); + if (tpcPtr != nullptr) + { + try + { + auto const& plane_obj = fWireReadout->Plane(tpcPtr->ID(), roi_properties.view); + theta_deg = static_cast(wmUtil.ThetaXW( + tv.dxdr, tv.dydr, tv.dzdr, plane_obj.ThetaZ()) * TMath::RadToDeg()); + } + catch (...) {} + } + } + tb_theta_xw.push_back(theta_deg); + tb_truth_x.push_back(truth_x_val); + } + + fROITree->Fill(); + } + + if (fSavePerHitData && it_hit_map != wmUtil.ROIMatchedHitMap.end()) + { + for (auto const& subroi_prop : subROIPropVec) + { + size_t i_subroi = subroi_prop.key.second; + if (i_subroi >= it_hit_map->second.size()) continue; + size_t hitVec_idx = it_hit_map->second[i_subroi]; + + auto const& scale_vals = SubROIMatchedScalesMap.at(subroi_prop.key); + (*vec_scale_q )[hitVec_idx] = static_cast(scale_vals.r_Q); + (*vec_scale_sigma)[hitVec_idx] = static_cast(scale_vals.r_sigma); + + auto truth_it = SubROIMatchedTruthMap.find(subroi_prop.key); + if (truth_it != SubROIMatchedTruthMap.end()) + { + auto const& tv = truth_it->second; + (*vec_truth_x)[hitVec_idx] = static_cast(tv.x); + geo::TPCGeo const* tpcPtr = fGeometry->PositionToTPCptr( + geo::Point_t{tv.x, tv.y, tv.z}); + if (tpcPtr) + { + try { + auto const& plane_obj = fWireReadout->Plane( + hitVec[hitVec_idx].WireID().planeID()); + double thetaZ = plane_obj.ThetaZ(); + (*vec_theta_xw)[hitVec_idx] = static_cast( + wmUtil.ThetaXW(tv.dxdr, tv.dydr, tv.dzdr, thetaZ) * TMath::RadToDeg()); + (*vec_dirx)[hitVec_idx] = static_cast(tv.dxdr); + (*vec_diry)[hitVec_idx] = static_cast(tv.dydr); + (*vec_dirz)[hitVec_idx] = static_cast(tv.dzdr); + double cosG = std::abs(tv.dydr * std::sin(thetaZ) + tv.dzdr * std::cos(thetaZ)); + float wp = static_cast(plane_obj.WirePitch()); + float pval = (cosG > 1e-6) ? wp / static_cast(cosG) : -999.f; + (*vec_pitch)[hitVec_idx] = pval; + (*vec_dqdx)[hitVec_idx] = (pval > 0.f) + ? hitVec[hitVec_idx].Integral() / pval : -999.f; + } + catch (...) {} + } + } + } + } + + std::vector mod_short(modified_data.size()); + std::transform(modified_data.begin(), modified_data.end(), mod_short.begin(), toShort); + new_rois.add_range(roi_properties.begin, mod_short); + } // end ROI loop + + new_chanROIs->emplace_back(new_rois, chanROI.Channel(), chanROI.ADCScaleFactor()); + + if (fSaveHistsByChannel && isModified) + { + readout::ROPID ropID = fWireReadout->ChannelToROP(chanROI.Channel()); + std::string titleStr = "Cryo-" + std::to_string(ropID.Cryostat) + + "_TPCset-" + std::to_string(ropID.TPCset) + + "_ReadOutPlane-" + std::to_string(ropID.ROP) + + "_Channel-" + std::to_string(chanROI.Channel()); + TH1F* oldChannelHist = new TH1F(("Old_" + titleStr).c_str(), ";Sample;Arbitrary Units", + wmUtil.readoutWindowTicks, 0, wmUtil.readoutWindowTicks); + TH1F* newChannelHist = new TH1F(("New_" + titleStr).c_str(), ";Sample;Arbitrary Units", + wmUtil.readoutWindowTicks, 0, wmUtil.readoutWindowTicks); + for (size_t tick = 0; tick < wmUtil.readoutWindowTicks; ++tick) + { + float oldSample = (tick < chanROI.Signal().size()) ? chanROI.Signal().at(tick) : 0; + float newSample = (tick < new_chanROIs->back().Signal().size()) ? new_chanROIs->back().Signal().at(tick) : 0; + oldChannelHist->SetBinContent(tick + 1, oldSample); + newChannelHist->SetBinContent(tick + 1, newSample); + } + TH1F* oldChannelHist_toSave = tfs->make(*oldChannelHist); + TH1F* newChannelHist_toSave = tfs->make(*newChannelHist); + mf::LogDebug("WireModifierXXW") + << "Saved histograms " << oldChannelHist_toSave->GetName() << '\n' + << " and " << newChannelHist_toSave->GetName(); + } + + if (fSaveHistsByWire && isModified) + { + std::vector wireIDs = fWireReadout->ChannelToWire(chanROI.Channel()); + mf::LogDebug("WireModifierXXW") + << "Channel " << chanROI.Channel() << " has " << wireIDs.size() << " wire(s)"; + for (auto const& wireID : wireIDs) + { + std::string titleStr = "Cryo-" + std::to_string(wireID.Cryostat) + + "_TPC-" + std::to_string(wireID.TPC) + + "_Plane-" + std::to_string(wireID.Plane) + + "_Wire-" + std::to_string(wireID.Wire); + TH1F* oldWireHist = new TH1F(("Old_" + titleStr).c_str(), ";Sample;Arbitrary Units", + wmUtil.readoutWindowTicks, 0, wmUtil.readoutWindowTicks); + TH1F* newWireHist = new TH1F(("New_" + titleStr).c_str(), ";Sample;Arbitrary Units", + wmUtil.readoutWindowTicks, 0, wmUtil.readoutWindowTicks); + for (size_t tick = 0; tick < wmUtil.readoutWindowTicks; ++tick) + { + float oldSample = (tick < chanROI.Signal().size()) ? chanROI.Signal().at(tick) : 0; + float newSample = (tick < new_chanROIs->back().Signal().size()) ? new_chanROIs->back().Signal().at(tick) : 0; + oldWireHist->SetBinContent(tick + 1, oldSample); + newWireHist->SetBinContent(tick + 1, newSample); + } + TH1F* oldWireHist_toSave = tfs->make(*oldWireHist); + TH1F* newWireHist_toSave = tfs->make(*newWireHist); + mf::LogDebug("WireModifierXXW") + << "Saved histograms " << oldWireHist_toSave->GetName() << '\n' + << " and " << newWireHist_toSave->GetName(); + } + } + + } // end loop over chanROIs + + evt.put(std::move(new_chanROIs)); + + } // end else (ChannelROI mode) + + evt.put(std::move(vec_scale_q), "scaleQ"); + evt.put(std::move(vec_scale_sigma), "scaleSigma"); + evt.put(std::move(vec_truth_x), "truthX"); + evt.put(std::move(vec_theta_xw), "thetaXW"); + evt.put(std::move(vec_is_mc), "isMC"); + evt.put(std::move(vec_dirx), "dirX"); + evt.put(std::move(vec_diry), "dirY"); + evt.put(std::move(vec_dirz), "dirZ"); + evt.put(std::move(vec_pitch), "pitch"); + evt.put(std::move(vec_dqdx), "dQdX"); } DEFINE_ART_MODULE(WireModifierXXW) } // end namespace diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl new file mode 100644 index 000000000..a43af0a65 --- /dev/null +++ b/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl @@ -0,0 +1,66 @@ +#include "WireModifierXXW_closure.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +# add the WireMod sequence to the standard stage1 sequence +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +# produce the WireMod products +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +# tell gaushit2d to use the ROIs from WireMod +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# tell gaushitPT2d to use the ROIs from WireMod +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +# allRatio_Coarse_XW.root has X/theta axes, not Y/Z: disable YZ correction +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +# tell the caloskimmer to use the ROIs from WireMod +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# drop products that we don't need +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl new file mode 100644 index 000000000..6ccb100b4 --- /dev/null +++ b/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl @@ -0,0 +1,72 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: true +physics.producers.WireModTPCEE.SetNullScaleWidth: true +physics.producers.WireModTPCEW.SetNullScaleIntegral: true +physics.producers.WireModTPCEW.SetNullScaleWidth: true +physics.producers.WireModTPCWE.SetNullScaleIntegral: true +physics.producers.WireModTPCWE.SetNullScaleWidth: true +physics.producers.WireModTPCWW.SetNullScaleIntegral: true +physics.producers.WireModTPCWW.SetNullScaleWidth: true + +# allRatio_Coarse_XW.root has (X, theta) axes, not (Y, Z): disable YZ correction +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl b/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl new file mode 100644 index 000000000..8c35060fb --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl @@ -0,0 +1,76 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: true +physics.producers.WireModTPCEE.SetNullScaleWidth: true +physics.producers.WireModTPCEW.SetNullScaleIntegral: true +physics.producers.WireModTPCEW.SetNullScaleWidth: true +physics.producers.WireModTPCWE.SetNullScaleIntegral: true +physics.producers.WireModTPCWE.SetNullScaleWidth: true +physics.producers.WireModTPCWW.SetNullScaleIntegral: true +physics.producers.WireModTPCWW.SetNullScaleWidth: true + +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# force cluster3D (Pandora/TrackCaloSkimmer) onto gaushit2dTPCXX, same hit-finder +# WireMod/WireModHitDiag use (default points it at gaushitPT2dTPCXX instead) +physics.producers.cluster3DCryoE.Hit3DBuilderAlg.HitFinderTagVec: [ "gaushit2dTPCEW", "gaushit2dTPCEE" ] +physics.producers.cluster3DCryoW.Hit3DBuilderAlg.HitFinderTagVec: [ "gaushit2dTPCWW", "gaushit2dTPCWE" ] + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl b/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl new file mode 100644 index 000000000..6722f54b7 --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl @@ -0,0 +1,76 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# force cluster3D (Pandora/TrackCaloSkimmer) onto gaushit2dTPCXX, same hit-finder +# WireMod/WireModHitDiag use (default points it at gaushitPT2dTPCXX instead) +physics.producers.cluster3DCryoE.Hit3DBuilderAlg.HitFinderTagVec: [ "gaushit2dTPCEW", "gaushit2dTPCEE" ] +physics.producers.cluster3DCryoW.Hit3DBuilderAlg.HitFinderTagVec: [ "gaushit2dTPCWW", "gaushit2dTPCWE" ] + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl b/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl new file mode 100644 index 000000000..5a4836801 --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl @@ -0,0 +1,91 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + dummyGaushitPT2dTPCEE, dummyGaushitPT2dTPCEW, dummyGaushitPT2dTPCWE, dummyGaushitPT2dTPCWW, + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +# "before" hits at the same LongMaxHits=10 setting as gaushitPT2dTPCXX, reading +# the undisturbed wire (not redirected to WMWireToROI) +physics.producers.dummyGaushitPT2dTPCEE: { @table::icarus_stage1_producers.gaushitPT2dTPCEE } +physics.producers.dummyGaushitPT2dTPCEW: { @table::icarus_stage1_producers.gaushitPT2dTPCEW } +physics.producers.dummyGaushitPT2dTPCWE: { @table::icarus_stage1_producers.gaushitPT2dTPCWE } +physics.producers.dummyGaushitPT2dTPCWW: { @table::icarus_stage1_producers.gaushitPT2dTPCWW } + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: true +physics.producers.WireModTPCEE.SetNullScaleWidth: true +physics.producers.WireModTPCEW.SetNullScaleIntegral: true +physics.producers.WireModTPCEW.SetNullScaleWidth: true +physics.producers.WireModTPCWE.SetNullScaleIntegral: true +physics.producers.WireModTPCWE.SetNullScaleWidth: true +physics.producers.WireModTPCWW.SetNullScaleIntegral: true +physics.producers.WireModTPCWW.SetNullScaleWidth: true + +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.producers.WireModTPCEE.HitLabel: "dummyGaushitPT2dTPCEE" +physics.producers.WireModTPCEW.HitLabel: "dummyGaushitPT2dTPCEW" +physics.producers.WireModTPCWE.HitLabel: "dummyGaushitPT2dTPCWE" +physics.producers.WireModTPCWW.HitLabel: "dummyGaushitPT2dTPCWW" + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analyzers.WireModHitDiag.BeforeHitLabels: [ "dummyGaushitPT2dTPCEE", "dummyGaushitPT2dTPCEW", "dummyGaushitPT2dTPCWE", "dummyGaushitPT2dTPCWW" ] +physics.analyzers.WireModHitDiag.AfterHitLabels: [ "gaushitPT2dTPCEE", "gaushitPT2dTPCEW", "gaushitPT2dTPCWE", "gaushitPT2dTPCWW" ] +physics.analyzers.WireModHitDiag.ScaleInfoLabels: [ "WireModTPCEE", "WireModTPCEW", "WireModTPCWE", "WireModTPCWW" ] + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# Pandora/TrackCaloSkimmer already default to gaushitPT2dTPCXX, no override needed + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_dummyGaushitPT2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl b/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl new file mode 100644 index 000000000..824b5a4bd --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl @@ -0,0 +1,91 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + dummyGaushitPT2dTPCEE, dummyGaushitPT2dTPCEW, dummyGaushitPT2dTPCWE, dummyGaushitPT2dTPCWW, + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +# "before" hits at the same LongMaxHits=10 setting as gaushitPT2dTPCXX, reading +# the undisturbed wire (not redirected to WMWireToROI) +physics.producers.dummyGaushitPT2dTPCEE: { @table::icarus_stage1_producers.gaushitPT2dTPCEE } +physics.producers.dummyGaushitPT2dTPCEW: { @table::icarus_stage1_producers.gaushitPT2dTPCEW } +physics.producers.dummyGaushitPT2dTPCWE: { @table::icarus_stage1_producers.gaushitPT2dTPCWE } +physics.producers.dummyGaushitPT2dTPCWW: { @table::icarus_stage1_producers.gaushitPT2dTPCWW } + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.producers.WireModTPCEE.HitLabel: "dummyGaushitPT2dTPCEE" +physics.producers.WireModTPCEW.HitLabel: "dummyGaushitPT2dTPCEW" +physics.producers.WireModTPCWE.HitLabel: "dummyGaushitPT2dTPCWE" +physics.producers.WireModTPCWW.HitLabel: "dummyGaushitPT2dTPCWW" + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analyzers.WireModHitDiag.BeforeHitLabels: [ "dummyGaushitPT2dTPCEE", "dummyGaushitPT2dTPCEW", "dummyGaushitPT2dTPCWE", "dummyGaushitPT2dTPCWW" ] +physics.analyzers.WireModHitDiag.AfterHitLabels: [ "gaushitPT2dTPCEE", "gaushitPT2dTPCEW", "gaushitPT2dTPCWE", "gaushitPT2dTPCWW" ] +physics.analyzers.WireModHitDiag.ScaleInfoLabels: [ "WireModTPCEE", "WireModTPCEW", "WireModTPCWE", "WireModTPCWW" ] + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# Pandora/TrackCaloSkimmer already default to gaushitPT2dTPCXX, no override needed + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_dummyGaushitPT2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" diff --git a/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl b/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl new file mode 100644 index 000000000..875aeaa53 --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl @@ -0,0 +1,82 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +# produce the WireMod products +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +# tell gaushit2d to use the ROIs from WireMod +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# tell gaushitPT2d to use the ROIs from WireMod +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# tell the caloskimmer to use the ROIs from WireMod +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] +physics.analyzers.WireModHitDiag: @local::WireModHitDiag +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +# drop products that we don't need +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] +# use additive ROI modification (after = before + [q_mod - q_orig]) +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +# enable full diagnostic output (ROITree, tick histograms, per-hit scale vectors) +physics.producers.WireModTPCEE.SaveROITree: true +physics.producers.WireModTPCEE.SaveTickHistograms: true +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEE.SaveByChannel: true +physics.producers.WireModTPCEW.SaveROITree: true +physics.producers.WireModTPCEW.SaveTickHistograms: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCEW.SaveByChannel: true +physics.producers.WireModTPCWE.SaveROITree: true +physics.producers.WireModTPCWE.SaveTickHistograms: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWE.SaveByChannel: true +physics.producers.WireModTPCWW.SaveROITree: true +physics.producers.WireModTPCWW.SaveTickHistograms: true +physics.producers.WireModTPCWW.SavePerHitData: true +physics.producers.WireModTPCWW.SaveByChannel: true + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "DEBUG" diff --git a/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl b/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl new file mode 100644 index 000000000..5da6d2a1b --- /dev/null +++ b/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl @@ -0,0 +1,73 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +# float (ChannelROI-direct) mode: no WMWireToROI, WireModTPCxx produces ChannelROI directly +physics.reco: +[ + @sequence::dummyGaushit2dSeq, + WMROI2Wire, + @sequence::WireModSeq, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +physics.producers.WireModTPCEE.UseChannelROIMode: true +physics.producers.WireModTPCEW.UseChannelROIMode: true +physics.producers.WireModTPCWE.UseChannelROIMode: true +physics.producers.WireModTPCWW.UseChannelROIMode: true + +physics.producers.WireModTPCEE: { @table::physics.producers.WireModTPCEE @table::WireModROITree } +physics.producers.WireModTPCEW: { @table::physics.producers.WireModTPCEW @table::WireModROITree } +physics.producers.WireModTPCWE: { @table::physics.producers.WireModTPCWE @table::WireModROITree } +physics.producers.WireModTPCWW: { @table::physics.producers.WireModTPCWW @table::WireModROITree } + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +# WireModTPCxx produces ChannelROI with empty instance label in float mode +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WireModTPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WireModTPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WireModTPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WireModTPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WireModTPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WireModTPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WireModTPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WireModTPCWW" + +# WMWireToROI doesn't run in float mode, redirect to original ChannelROI +physics.analyzers.caloskimE.RawDigitproducers: [ "wire2channelroi2d:PHYSCRATEDATATPCEW", + "wire2channelroi2d:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "wire2channelroi2d:PHYSCRATEDATATPCWW", + "wire2channelroi2d:PHYSCRATEDATATPCWE" ] + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" +services.message.destinations.STDCOUT.threshold: "DEBUG" diff --git a/icaruscode/WireMod/wireMod_stage1.fcl b/icaruscode/WireMod/wireMod_stage1.fcl new file mode 100644 index 000000000..9cac89982 --- /dev/null +++ b/icaruscode/WireMod/wireMod_stage1.fcl @@ -0,0 +1,47 @@ +#include "WireModifierXXW_closure.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +# add the WireMod sequence to the standard stage1 sequence +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +# produce the WireMod products +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +# tell gaushit2d to use the ROIs from WireMod +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# tell gaushitPT2d to use the ROIs from WireMod +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# tell the caloskimmer to use the ROIs from WireMod +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +# drop products that we don't need +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModTPC*_*_*" +] +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "DEBUG" + diff --git a/icaruscode/WireMod/wiremod_scale_debug.fcl b/icaruscode/WireMod/wiremod_scale_debug.fcl new file mode 100644 index 000000000..cb0bded6a --- /dev/null +++ b/icaruscode/WireMod/wiremod_scale_debug.fcl @@ -0,0 +1,26 @@ +# EmptyEvent job: runs WireModScaleDebug once in beginJob(), no input file needed. +# Run from the directory containing Run4_Coarse_XW.root (LocalRatios convention). + +#include "stage1_icarus_wiremod_pops_hitdiag.fcl" + +source: { + module_type: EmptyEvent + maxEvents: 1 +} + +physics: { + analyzers: { + scaledebug: { + module_type: "WireModScaleDebug" + RatioFileName_XXW: "Run4_Coarse_XW.root" + ThetaEdges: [0, 6, 12, 18, 24, 30, 38, 46, 56, 90] + } + } + ana: [ scaledebug ] + trigger_paths: [] + end_paths: [ ana ] +} + +outputs: {} + +services.message.destinations.STDCOUT.threshold: "ERROR" From 1560a43371566f4265f31a80e178a1aebca2dbae Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Tue, 30 Jun 2026 17:02:42 -0500 Subject: [PATCH 2/9] Build fixes. --- icaruscode/WireMod/WireModifierXXW_module.cc | 49 +++++++------------- 1 file changed, 16 insertions(+), 33 deletions(-) diff --git a/icaruscode/WireMod/WireModifierXXW_module.cc b/icaruscode/WireMod/WireModifierXXW_module.cc index 542062e6a..7c57a53fa 100644 --- a/icaruscode/WireMod/WireModifierXXW_module.cc +++ b/icaruscode/WireMod/WireModifierXXW_module.cc @@ -679,28 +679,23 @@ namespace wiremod catch(...) { continue; // ignore depositions outside TPC } - for (auto const& plane : fWireReadout->Iterate(curTPCGeomID)) - { - if (fWireReadout->FindTPCsetAtPosition(edep.MidPoint()) != readout::TPCsetID(fCryo, fTPCset)) - continue; - geo::TPCGeo const* curTPCGeomPtr = fGeometry->PositionToTPCptr(edep.MidPoint()); - for (auto const& plane : fWireReadout->Iterate(curTPCGeomPtr->ID())) - { - TH2F* targetHist = (plane.View() == geo::kY) ? fEdepTickInd1 - : (plane.View() == geo::kV) ? fEdepTickInd2 - : (plane.View() == geo::kU) ? fEdepTickColl - : nullptr; - if (targetHist == nullptr) continue; - - geo::WireID wireID; - try { wireID = plane.NearestWireID(edep.MidPoint()); } - catch (...) { continue; } - float projTick = detProp.ConvertXToTicks(edep.X(), wireID.Plane, wireID.TPC, wireID.Cryostat) + wmUtil.tickOffset; - float projChan = fWireReadout->PlaneWireToChannel(wireID); - targetHist->Fill(projTick, projChan); - fEdepXvsChan->Fill(edep.X(), projChan); - } + geo::TPCGeo const* curTPCGeomPtr = fGeometry->PositionToTPCptr(edep.MidPoint()); + for (auto const& plane : fWireReadout->Iterate(curTPCGeomPtr->ID())) + { + TH2F* targetHist = (plane.View() == geo::kY) ? fEdepTickInd1 + : (plane.View() == geo::kV) ? fEdepTickInd2 + : (plane.View() == geo::kU) ? fEdepTickColl + : nullptr; + if (targetHist == nullptr) continue; + + geo::WireID wireID; + try { wireID = plane.NearestWireID(edep.MidPoint()); } + catch (...) { continue; } + float projTick = detProp.ConvertXToTicks(edep.X(), wireID.Plane, wireID.TPC, wireID.Cryostat) + wmUtil.tickOffset; + float projChan = fWireReadout->PlaneWireToChannel(wireID); + targetHist->Fill(projTick, projChan); + fEdepXvsChan->Fill(edep.X(), projChan); } } for (auto const& hit : hitVec) @@ -724,18 +719,6 @@ namespace wiremod } catch (...) { } fHitXvsChan->Fill(hitX, chan); - - geo::WireID wireID; - try { - wireID = plane.NearestWireID(edep.MidPoint()); - } - catch (...) { - continue; // don't fill non-active depositions - } - float projTick = detProp.ConvertXToTicks(edep.X(), wireID.Plane, wireID.TPC, wireID.Cryostat) + wmUtil.tickOffset; - float projChan = fWireReadout->PlaneWireToChannel(wireID); - targetHist->Fill(projTick, projChan); - fEdepXvsChan->Fill(edep.X(), projChan); } } From 4ae2154c5dfdbdec2619819cc8678dbbed83b71b Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Jul 2026 08:59:12 -0500 Subject: [PATCH 3/9] Change truth matching to use SimChannel's instead of SimEnergyDeposit's --- icaruscode/WireMod/WireModifierXXW_module.cc | 211 +++++++++++-------- 1 file changed, 118 insertions(+), 93 deletions(-) diff --git a/icaruscode/WireMod/WireModifierXXW_module.cc b/icaruscode/WireMod/WireModifierXXW_module.cc index 7c57a53fa..a83e1937b 100644 --- a/icaruscode/WireMod/WireModifierXXW_module.cc +++ b/icaruscode/WireMod/WireModifierXXW_module.cc @@ -28,6 +28,10 @@ #include "lardata/DetectorInfoServices/DetectorPropertiesService.h" #include "lardata/DetectorInfoServices/DetectorClocksService.h" #include "larsim/MCCheater/BackTrackerService.h" +#include "larevt/SpaceCharge/SpaceCharge.h" +#include "larevt/SpaceChargeServices/SpaceChargeService.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "nusimdata/SimulationBase/MCParticle.h" #include "lardata/Utilities/AssociationUtil.h" #include "lardataobj/RawData/RawDigit.h" #include "lardataobj/RecoBase/Wire.h" @@ -63,8 +67,10 @@ namespace wiremod std::vector fGraph_sigma_YZ; art::InputTag fWireLabel; // which wires are we pulling in? art::InputTag fHitLabel; // which hits are we pulling in? - art::InputTag fEDepLabel; // which are the EDeps? - art::InputTag fEDepShiftedLabel; // which are the EDeps? + art::InputTag fEDepLabel; // which are the EDeps? (legacy; now only the tick diagnostic) + art::InputTag fEDepShiftedLabel; // which are the EDeps? (legacy; now only the tick diagnostic) + art::InputTag fSimChannelLabel; // SimChannels for IDE-based scale truth (default "merge") + art::InputTag fG4Label; // MCParticles for trajectory direction (default "largeant") unsigned int fCryo; // which Cryo are we in? unsigned int fTPCset; // which TPC are we in? bool fLocalRatios; // is the ratio file local? @@ -201,6 +207,8 @@ namespace wiremod fHitLabel = pset.get("HitLabel"); fEDepLabel = pset.get("EDepLabel"); fEDepShiftedLabel = pset.get("EDepShiftedLabel"); + fSimChannelLabel = pset.get("SimChannelLabel", art::InputTag("merge")); + fG4Label = pset.get("G4Label", art::InputTag("largeant")); fCryo = pset.get("Cryo"); fTPCset = pset.get("TPCset"); @@ -460,6 +468,11 @@ namespace wiremod produces>("thetaXW"); // per-hit MC flag: 1 = MC hit (BackTracker matched), 0 = overlay data hit produces>("isMC"); + // per-hit SimChannel/BackTracker truth, computed identically to TrackCaloSkimmer::MakeHit + // (ChannelToTrackIDEs over the hit [StartTick,EndTick] window): total true energy [MeV] + // and number of ionization electrons summed over all TrackIDEs under the hit. + produces>("truthE"); + produces>("truthNelec"); // per-hit truth direction and derived geometry (filled only when SavePerHitData: true) produces>("dirX"); // truth dxdr produces>("dirY"); // truth dydr @@ -574,6 +587,15 @@ namespace wiremod evt.getByLabel(fEDepShiftedLabel, edepShiftedHandle); auto const& edepShiftedVec(*edepShiftedHandle); + // SimChannel + MCParticle inputs: the scale-value truth (position/direction/dE-dQ per + // sub-ROI) is derived from these, replacing the SimEnergyDeposit path. The SimChannel + // ionization gives the charge-weighted true position; the MCParticle trajectory gives + // the local track direction used by the ThetaXW/XXW scale lookup. + auto const& simchVec = *evt.getValidHandle>(fSimChannelLabel); + auto const& mcpVec = *evt.getValidHandle>(fG4Label); + std::map particleMap; + for (auto const& p : mcpVec) particleMap[p.TrackId()] = &p; + art::Handle< std::vector > hitHandle; evt.getByLabel(fHitLabel, hitHandle); auto const& hitVec(*hitHandle); @@ -590,15 +612,55 @@ namespace wiremod auto vec_dirz = std::make_unique>(perHitN, -999.f); auto vec_pitch = std::make_unique>(perHitN, -999.f); auto vec_dqdx = std::make_unique>(perHitN, -999.f); - - // Get the tick offset by comparing the hit ticks with the tick gotten by - // projecting the backtracked hit X to a tick. Yes this is silly thanks for asking. - // Also fills vec_is_mc: 1 if BackTracker finds a match (MC hit), 0 otherwise (overlay). + // SimChannel/BackTracker truth per hit (TrackCaloSkimmer-identical), -1 when no truth. + auto vec_truth_e = std::make_unique>(perHitN, -1.f); + auto vec_truth_nelec = std::make_unique>(perHitN, -1.f); + + // Per-hit truth matching, computed IDENTICALLY to TrackCaloSkimmer::MakeHit: + // ChannelToTrackIDEs(channel, StartTick, EndTick) summed over all TrackIDEs. + // isMC is now defined as "this window has any SimChannel ionization" (i.e. the hit is + // truth-matched), matching the TrackCaloSkimmer notion of a matched hit. A hit with no + // IDEs in its window is an overlay/data hit (isMC = 0, truthE/truthNelec = -1). + // + // The BackTracker HitToXYZ projection below is retained only to derive the legacy + // SimEnergyDeposit tick offset (BT_Offset) used by the edep scale path; it is removed + // once the scale path is migrated to SimChannel (see CalcPropertiesFromIDEs). double BT_Offset = 0; unsigned int BT_Hits = 0; for (size_t i_h = 0; i_h < hitVec.size(); ++i_h) { auto const& hit = hitVec[i_h]; + + // --- TrackCaloSkimmer-identical per-hit SimChannel truth --- + std::vector ides; + try + { + ides = fBT->ChannelToTrackIDEs(fDetClocksData, hit.Channel(), hit.StartTick(), hit.EndTick()); + } + catch (...) + { + ides.clear(); + } + if (fSavePerHitData) + { + if (!ides.empty()) + { + (*vec_is_mc)[i_h] = 1; + (*vec_truth_e)[i_h] = 0.f; + (*vec_truth_nelec)[i_h] = 0.f; + for (auto const& ide : ides) + { + (*vec_truth_e)[i_h] += ide.energy; + (*vec_truth_nelec)[i_h] += ide.numElectrons; + } + } + else + { + (*vec_is_mc)[i_h] = 0; // overlay/data hit (no SimChannel ionization in window) + } + } + + // --- legacy SimEnergyDeposit tick offset (BackTracker HitToXYZ) --- double hitX = std::numeric_limits::max(); try { @@ -606,16 +668,16 @@ namespace wiremod } catch (...) { - // Overlay hit + // Overlay hit -- no contribution to the tick offset continue; } - if (fSavePerHitData) (*vec_is_mc)[i_h] = 1; double hitTick = hit.PeakTime(); double projTick = detProp.ConvertXToTicks(hitX, hit.WireID().Plane, hit.WireID().TPC, hit.WireID().Cryostat); BT_Offset += (hitTick - projTick); ++BT_Hits; } if (BT_Hits > 0) BT_Offset /= BT_Hits; + (void)BT_Offset; // legacy edep tick offset; unused by the SimChannel/IDE scale path sys::WireModUtility wmUtil(fGeometry, fWireReadout, detProp, false, // Channel Scale @@ -628,7 +690,10 @@ namespace wiremod false, // X-dQ/dX false, // XZAngledQdX (fRatioFileName_XXW != "NOFILE"), // X-ThXW - BT_Offset); // Tick Offset + 0.0); // Tick Offset: IDEs carry their native readout tick (no offset) + // Space-charge provider for mapping IDE (at-the-wire) positions to the true trajectory + // frame inside CalcPropertiesFromIDEs. Null/disabled -> identity. + wmUtil.fSCE = lar::providerFrom(); if (fRatioFileName_YZ != "NOFILE") { wmUtil.graph2Ds_Charge_YZ = fGraph_charge_YZ; @@ -729,9 +794,9 @@ namespace wiremod auto const& wireVec(*wireHandle); std::unique_ptr> new_wires(new std::vector()); - mf::LogVerbatim("WireModifierXXW") << "Get Edep Map"; - wmUtil.FillROIMatchedEdepMap(edepShiftedVec, wireVec, offset_ADC); - mf::LogVerbatim("WireModifierXXW") << "Got Edep Map." << '\n' << "Get Hit Map"; + mf::LogVerbatim("WireModifierXXW") << "Get IDE Map"; + wmUtil.FillROIMatchedIDEMap(simchVec, wireVec, fDetClocksData, offset_ADC); + mf::LogVerbatim("WireModifierXXW") << "Got IDE Map." << '\n' << "Get Hit Map"; wmUtil.FillROIMatchedHitMap(hitVec, wireVec); mf::LogVerbatim("WireModifierXXW") << "Got Hit Map."; @@ -787,30 +852,26 @@ namespace wiremod std::vector modified_data(range.data()); - auto it_map = wmUtil.ROIMatchedEdepMap.find(roi_key); - if(it_map==wmUtil.ROIMatchedEdepMap.end()){ + auto it_map = wmUtil.ROIMatchedIDEMap.find(roi_key); + if(it_map==wmUtil.ROIMatchedIDEMap.end()){ new_rois.add_range(range.begin_index(), modified_data); mf::LogDebug("WireModifierXXW") - << " Could not find matching Edep. Skip"; + << " Could not find matching IDE. Skip"; continue; } - std::vector matchedEdepIdxVec = it_map->second; - if(matchedEdepIdxVec.size() == 0) + std::vector matchedIDEIdxVec = it_map->second; + if(matchedIDEIdxVec.size() == 0) { new_rois.add_range(range.begin_index(), modified_data); mf::LogDebug("WireModifierXXW") - << " No indices for Edep. Skip"; + << " No indices for IDE. Skip"; continue; } - std::vector matchedEdepPtrVec; - std::vector matchedEdepShiftedPtrVec; - for(auto i_e : matchedEdepIdxVec) - { - matchedEdepPtrVec.push_back(&edepVec[i_e]); - matchedEdepShiftedPtrVec.push_back(&edepShiftedVec[i_e]); - } + std::vector matchedIDEPtrVec; + for(auto i_e : matchedIDEIdxVec) + matchedIDEPtrVec.push_back(&wmUtil.fIDEVec[i_e]); mf::LogDebug("WireModifierXXW") - << " Found " << matchedEdepPtrVec.size() << " Edeps"; + << " Found " << matchedIDEPtrVec.size() << " IDEs"; std::vector matchedHitPtrVec; auto it_hit_map = wmUtil.ROIMatchedHitMap.find(roi_key); @@ -839,25 +900,8 @@ namespace wiremod mf::LogDebug("WireModifierXXW") << " have " << subROIPropVec.size() << " SubROI"; - auto SubROIMatchedEdepShiftedMap = - wmUtil.MatchEdepsToSubROIs(subROIPropVec, matchedEdepShiftedPtrVec, offset_ADC); - std::map> - SubROIMatchedEdepMap; - for (auto const& key_edepPtrVec_pair : SubROIMatchedEdepShiftedMap) - { - auto key = key_edepPtrVec_pair.first; - for (auto const& shifted_edep_ptr : key_edepPtrVec_pair.second) - { - for (size_t i_e = 0; i_e < matchedEdepShiftedPtrVec.size(); ++i_e) - { - if (shifted_edep_ptr == matchedEdepShiftedPtrVec[i_e]) - { - SubROIMatchedEdepMap[key].push_back(matchedEdepPtrVec[i_e]); - break; - } - } - } - } + auto SubROIMatchedIDEMap = + wmUtil.MatchIDEsToSubROIs(subROIPropVec, matchedIDEPtrVec); std::map SubROIMatchedScalesMap; @@ -868,12 +912,12 @@ namespace wiremod { sys::WireModUtility::ScaleValues_t scale_vals; auto key = subroi_prop.key; - auto key_it = SubROIMatchedEdepMap.find(key); + auto key_it = SubROIMatchedIDEMap.find(key); int scale_reason = 0; - if (key_it != SubROIMatchedEdepMap.end() && key_it->second.size() > 0) + if (key_it != SubROIMatchedIDEMap.end() && key_it->second.size() > 0) { - auto truth_vals = wmUtil.CalcPropertiesFromEdeps(key_it->second, offset_ADC); + auto truth_vals = wmUtil.CalcPropertiesFromIDEs(key_it->second, particleMap); SubROIMatchedTruthMap[key] = truth_vals; if (truth_vals.total_energy < 0.3 && subroi_prop.total_q > 80) @@ -1002,12 +1046,12 @@ namespace wiremod // dominant track ID: pick the TrackID with highest summed energy int dominant_tid = -999; - auto edep_it = SubROIMatchedEdepMap.find(key); - if (edep_it != SubROIMatchedEdepMap.end() && !edep_it->second.empty()) + auto edep_it = SubROIMatchedIDEMap.find(key); + if (edep_it != SubROIMatchedIDEMap.end() && !edep_it->second.empty()) { std::map tidEnergy; for (auto const& ep : edep_it->second) - tidEnergy[ep->TrackID()] += ep->E(); + tidEnergy[std::abs(ep->ide->trackID)] += ep->ide->energy; double maxE = 0.; for (auto const& te : tidEnergy) if (te.second > maxE) { maxE = te.second; dominant_tid = te.first; } @@ -1162,9 +1206,9 @@ namespace wiremod auto const& chanROIVec(*chanROIHandle); std::unique_ptr> new_chanROIs(new std::vector()); - mf::LogVerbatim("WireModifierXXW") << "Get Edep Map"; - wmUtil.FillROIMatchedEdepMap(edepShiftedVec, chanROIVec, offset_ADC); - mf::LogVerbatim("WireModifierXXW") << "Got Edep Map." << '\n' << "Get Hit Map"; + mf::LogVerbatim("WireModifierXXW") << "Get IDE Map"; + wmUtil.FillROIMatchedIDEMap(simchVec, chanROIVec, fDetClocksData, offset_ADC); + mf::LogVerbatim("WireModifierXXW") << "Got IDE Map." << '\n' << "Get Hit Map"; wmUtil.FillROIMatchedHitMap(hitVec, chanROIVec); mf::LogVerbatim("WireModifierXXW") << "Got Hit Map."; @@ -1216,31 +1260,27 @@ namespace wiremod return static_cast(std::round(v * chanROI.ADCScaleFactor())); }; - auto it_map = wmUtil.ROIMatchedEdepMap.find(roi_key); - if(it_map==wmUtil.ROIMatchedEdepMap.end()){ + auto it_map = wmUtil.ROIMatchedIDEMap.find(roi_key); + if(it_map==wmUtil.ROIMatchedIDEMap.end()){ std::vector pass_short(modified_data.size()); std::transform(modified_data.begin(), modified_data.end(), pass_short.begin(), toShort); new_rois.add_range(range_short.begin_index(), pass_short); - mf::LogDebug("WireModifierXXW") << " Could not find matching Edep. Skip"; + mf::LogDebug("WireModifierXXW") << " Could not find matching IDE. Skip"; continue; } - std::vector matchedEdepIdxVec = it_map->second; - if(matchedEdepIdxVec.size() == 0) + std::vector matchedIDEIdxVec = it_map->second; + if(matchedIDEIdxVec.size() == 0) { std::vector pass_short(modified_data.size()); std::transform(modified_data.begin(), modified_data.end(), pass_short.begin(), toShort); new_rois.add_range(range_short.begin_index(), pass_short); - mf::LogDebug("WireModifierXXW") << " No indices for Edep. Skip"; + mf::LogDebug("WireModifierXXW") << " No indices for IDE. Skip"; continue; } - std::vector matchedEdepPtrVec; - std::vector matchedEdepShiftedPtrVec; - for(auto i_e : matchedEdepIdxVec) - { - matchedEdepPtrVec.push_back(&edepVec[i_e]); - matchedEdepShiftedPtrVec.push_back(&edepShiftedVec[i_e]); - } - mf::LogDebug("WireModifierXXW") << " Found " << matchedEdepPtrVec.size() << " Edeps"; + std::vector matchedIDEPtrVec; + for(auto i_e : matchedIDEIdxVec) + matchedIDEPtrVec.push_back(&wmUtil.fIDEVec[i_e]); + mf::LogDebug("WireModifierXXW") << " Found " << matchedIDEPtrVec.size() << " IDEs"; std::vector matchedHitPtrVec; auto it_hit_map = wmUtil.ROIMatchedHitMap.find(roi_key); @@ -1265,25 +1305,8 @@ namespace wiremod auto subROIPropVec = wmUtil.CalcSubROIProperties(roi_properties, matchedHitPtrVec); mf::LogDebug("WireModifierXXW") << " have " << subROIPropVec.size() << " SubROI"; - auto SubROIMatchedEdepShiftedMap = - wmUtil.MatchEdepsToSubROIs(subROIPropVec, matchedEdepShiftedPtrVec, offset_ADC); - std::map> - SubROIMatchedEdepMap; - for (auto const& key_edepPtrVec_pair : SubROIMatchedEdepShiftedMap) - { - auto key = key_edepPtrVec_pair.first; - for (auto const& shifted_edep_ptr : key_edepPtrVec_pair.second) - { - for (size_t i_e = 0; i_e < matchedEdepShiftedPtrVec.size(); ++i_e) - { - if (shifted_edep_ptr == matchedEdepShiftedPtrVec[i_e]) - { - SubROIMatchedEdepMap[key].push_back(matchedEdepPtrVec[i_e]); - break; - } - } - } - } + auto SubROIMatchedIDEMap = + wmUtil.MatchIDEsToSubROIs(subROIPropVec, matchedIDEPtrVec); std::map SubROIMatchedScalesMap; @@ -1294,12 +1317,12 @@ namespace wiremod { sys::WireModUtility::ScaleValues_t scale_vals; auto key = subroi_prop.key; - auto key_it = SubROIMatchedEdepMap.find(key); + auto key_it = SubROIMatchedIDEMap.find(key); int scale_reason = 0; - if (key_it != SubROIMatchedEdepMap.end() && key_it->second.size() > 0) + if (key_it != SubROIMatchedIDEMap.end() && key_it->second.size() > 0) { - auto truth_vals = wmUtil.CalcPropertiesFromEdeps(key_it->second, offset_ADC); + auto truth_vals = wmUtil.CalcPropertiesFromIDEs(key_it->second, particleMap); SubROIMatchedTruthMap[key] = truth_vals; if (truth_vals.total_energy < 0.3 && subroi_prop.total_q > 80) @@ -1423,12 +1446,12 @@ namespace wiremod tb_scale_sigma.push_back(static_cast(scale_vals.r_sigma)); int dominant_tid = -999; - auto edep_it = SubROIMatchedEdepMap.find(key); - if (edep_it != SubROIMatchedEdepMap.end() && !edep_it->second.empty()) + auto edep_it = SubROIMatchedIDEMap.find(key); + if (edep_it != SubROIMatchedIDEMap.end() && !edep_it->second.empty()) { std::map tidEnergy; for (auto const& ep : edep_it->second) - tidEnergy[ep->TrackID()] += ep->E(); + tidEnergy[std::abs(ep->ide->trackID)] += ep->ide->energy; double maxE = 0.; for (auto const& te : tidEnergy) if (te.second > maxE) { maxE = te.second; dominant_tid = te.first; } @@ -1579,6 +1602,8 @@ namespace wiremod evt.put(std::move(vec_truth_x), "truthX"); evt.put(std::move(vec_theta_xw), "thetaXW"); evt.put(std::move(vec_is_mc), "isMC"); + evt.put(std::move(vec_truth_e), "truthE"); + evt.put(std::move(vec_truth_nelec), "truthNelec"); evt.put(std::move(vec_dirx), "dirX"); evt.put(std::move(vec_diry), "dirY"); evt.put(std::move(vec_dirz), "dirZ"); From 1904b78cdde972776e993902dadadeaaa0430831 Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Jul 2026 08:59:54 -0500 Subject: [PATCH 4/9] Add thxw resolution check to Diagnostic module. --- .../WireMod/WireModHitDiagnostic_module.cc | 123 ++++++++++++++++++ 1 file changed, 123 insertions(+) diff --git a/icaruscode/WireMod/WireModHitDiagnostic_module.cc b/icaruscode/WireMod/WireModHitDiagnostic_module.cc index 1e1b5b45b..34418b620 100644 --- a/icaruscode/WireMod/WireModHitDiagnostic_module.cc +++ b/icaruscode/WireMod/WireModHitDiagnostic_module.cc @@ -6,25 +6,45 @@ // ROOT includes #include "TTree.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TMath.h" // art includes #include "art/Framework/Core/EDAnalyzer.h" #include "art/Framework/Core/ModuleMacros.h" #include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" #include "art/Framework/Services/Registry/ServiceHandle.h" #include "art_root_io/TFileService.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" #include "canvas/Utilities/InputTag.h" #include "fhiclcpp/ParameterSet.h" #include "messagefacility/MessageLogger/MessageLogger.h" // larsoft includes #include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/TrackHitMeta.h" #include "larcore/Geometry/WireReadout.h" #include "larcorealg/Geometry/WireReadoutGeom.h" #include "larcore/CoreUtils/ServiceUtil.h" +#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" // util::pi namespace wiremod { +// WireMod ThetaXW (angle of the track to the drift/X direction in the wire-projected +// frame), replicated from sys::WireModUtility::ThetaXW so the reco-track angle is computed +// with exactly the same convention as the truth angle stored by WireModifierXXW. +inline double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) +{ + const double s = std::sin(planeAngle - util::pi()); + const double c = std::cos(planeAngle - util::pi()); + const double cosG = std::abs(dydr * s + dzdr * c); + return (cosG > 1e-9) ? std::abs(std::atan(dxdr / cosG)) : 0.5 * util::pi(); +} + class WireModHitDiagnostic : public art::EDAnalyzer { public: @@ -35,12 +55,20 @@ class WireModHitDiagnostic : public art::EDAnalyzer std::vector fBeforeHitLabels; // dummyGaushit2dTPC* std::vector fAfterHitLabels; // gaushit2dTPC* std::vector fScaleInfoLabels; // WireModTPC* (produce scale products) + std::vector fTrackLabels; // reco tracks for reco ThetaXW (pandoraTrackGaus*) float fMatchWindowTicks; bool fOnlyMC; // if true, skip overlay hits (is_mc==0) when filling HitMatchTree const geo::WireReadoutGeom* fWireReadout = &(art::ServiceHandle()->Get()); + // True (WireMod/MCParticle) vs reco (Pandora track) ThetaXW resolution histograms [deg]. + // Filled for MC hits that are on a reconstructed track: the true angle is the WireModifierXXW + // per-hit thetaXW; the reco angle is ThetaXW of the track direction at that hit. + TH2F* fThxwTrueVsReco = nullptr; // all planes + TH2F* fThxwTrueVsReco_p[3] = {nullptr, nullptr, nullptr}; + TH1F* fThxwResid = nullptr; // reco - true [deg] + // HitMatchTree: one entry per before-hit TTree* fHitMatchTree; Int_t tm_run, tm_subrun, tm_event; @@ -61,6 +89,10 @@ class WireModHitDiagnostic : public art::EDAnalyzer Float_t tm_truth_dirz; // truth dzdr Float_t tm_pitch; // wirePitch / cosG [cm] Float_t tm_dqdx_before; // integral_before / pitch [ADC/cm], no lifetime correction + // SimChannel/BackTracker truth from WireModifierXXW (TrackCaloSkimmer-identical recipe), + // -1 if the before-hit window has no SimChannel ionization (overlay/data hit) + Float_t tm_truth_e; // total true energy under the before-hit [MeV] + Float_t tm_truth_nelec; // total ionization electrons under the before-hit // match result Int_t tm_is_matched; Int_t tm_n_after_matches; @@ -94,9 +126,25 @@ WireModHitDiagnostic::WireModHitDiagnostic(fhicl::ParameterSet const& pset) for (auto const& s : pset.get>("ScaleInfoLabels", std::vector{})) fScaleInfoLabels.emplace_back(s); + for (auto const& s : pset.get>("TrackLabels", + std::vector{})) + fTrackLabels.emplace_back(s); art::ServiceHandle tfs; + // True-vs-reco ThetaXW resolution histograms (degrees) + const int nb = 90; + const double lo = 0.0, hi = 90.0; + fThxwTrueVsReco = tfs->make("ThxwTrueVsReco", + "True vs reco #theta_{XW} (all planes);true #theta_{XW} [deg];reco #theta_{XW} [deg]", + nb, lo, hi, nb, lo, hi); + for (int p = 0; p < 3; ++p) + fThxwTrueVsReco_p[p] = tfs->make(Form("ThxwTrueVsReco_p%d", p), + Form("True vs reco #theta_{XW} (plane %d);true #theta_{XW} [deg];reco #theta_{XW} [deg]", p), + nb, lo, hi, nb, lo, hi); + fThxwResid = tfs->make("ThxwResid", + "reco - true #theta_{XW};#Delta#theta_{XW} [deg];hits", 200, -45.0, 45.0); + fHitMatchTree = tfs->make("HitMatchTree", "WireMod before-after hit comparison"); fHitMatchTree->Branch("run", &tm_run, "run/I"); @@ -126,6 +174,8 @@ WireModHitDiagnostic::WireModHitDiagnostic(fhicl::ParameterSet const& pset) fHitMatchTree->Branch("truth_dirz", &tm_truth_dirz, "truth_dirz/F"); fHitMatchTree->Branch("pitch", &tm_pitch, "pitch/F"); fHitMatchTree->Branch("dqdx_before", &tm_dqdx_before, "dqdx_before/F"); + fHitMatchTree->Branch("truth_e", &tm_truth_e, "truth_e/F"); + fHitMatchTree->Branch("truth_nelec", &tm_truth_nelec, "truth_nelec/F"); fHitMatchTree->Branch("is_matched", &tm_is_matched, "is_matched/I"); fHitMatchTree->Branch("n_after_matches", &tm_n_after_matches, "n_after_matches/I"); fHitMatchTree->Branch("after_peak_time", &tm_after_peak_time, "after_peak_time/F"); @@ -189,6 +239,9 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) afterHits.push_back(&hit); } + // per-channel truth ThetaXW (from WireMod, on before-hits): channel -> (peaktime, plane, true_thxw_deg) + std::unordered_multimap> truthThxwByChannel; + // Build channel -> [after-hit indices] index std::unordered_multimap afterByChannel; afterByChannel.reserve(afterHits.size() * 2); @@ -221,6 +274,8 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) std::vector const* pDirZ = nullptr; std::vector const* pPitch = nullptr; std::vector const* pDQdX = nullptr; + std::vector const* pTruthE = nullptr; + std::vector const* pTruthNelec = nullptr; if (j < fScaleInfoLabels.size()) { @@ -235,6 +290,8 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) auto dzH = evt.getHandle>(art::InputTag(siLabel.label(), "dirZ")); auto ptH = evt.getHandle>(art::InputTag(siLabel.label(), "pitch")); auto dqH = evt.getHandle>(art::InputTag(siLabel.label(), "dQdX")); + auto teH = evt.getHandle>(art::InputTag(siLabel.label(), "truthE")); + auto tnH = evt.getHandle>(art::InputTag(siLabel.label(), "truthNelec")); if (sqH.isValid()) pScaleQ = &(*sqH); if (ssH.isValid()) pScaleSigma = &(*ssH); if (txH.isValid()) pTruthX = &(*txH); @@ -245,6 +302,8 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) if (dzH.isValid()) pDirZ = &(*dzH); if (ptH.isValid()) pPitch = &(*ptH); if (dqH.isValid()) pDQdX = &(*dqH); + if (teH.isValid()) pTruthE = &(*teH); + if (tnH.isValid()) pTruthNelec = &(*tnH); } for (size_t i = 0; i < localHits.size(); ++i) @@ -278,6 +337,14 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) tm_truth_dirz = (pDirZ && i < pDirZ->size()) ? (*pDirZ)[i] : -999.f; tm_pitch = (pPitch && i < pPitch->size()) ? (*pPitch)[i] : -999.f; tm_dqdx_before = (pDQdX && i < pDQdX->size()) ? (*pDQdX)[i] : -999.f; + tm_truth_e = (pTruthE && i < pTruthE->size()) ? (*pTruthE)[i] : -1.f; + tm_truth_nelec = (pTruthNelec && i < pTruthNelec->size()) ? (*pTruthNelec)[i] : -1.f; + + // Record the truth ThetaXW for this (MC) hit so the reco-track loop below can pair it + // with the reconstructed track angle on the same channel/tick. + if (tm_theta_xw > -900.f) + truthThxwByChannel.emplace(bh.Channel(), + std::make_tuple(tm_before_peak_time, tm_plane, tm_theta_xw)); // if OnlyMC is set, skip overlay hits if (fOnlyMC && tm_is_mc != 1) continue; @@ -365,6 +432,62 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) tn_dof = ah->DegreesOfFreedom(); fNewHitTree->Fill(); } + + // --------------------------------------------------------------------------- + // True-vs-reco ThetaXW: for every hit on a reconstructed track, compute the reco + // ThetaXW from the track direction at that hit, pair it with the WireMod truth + // ThetaXW on the same channel/tick, and fill the resolution histograms. + // --------------------------------------------------------------------------- + for (auto const& tag : fTrackLabels) + { + auto trackHandle = evt.getHandle>(tag); + if (!trackHandle.isValid()) + { + mf::LogWarning("WireModHitDiagnostic") << "Track collection not found: " << tag; + continue; + } + art::FindManyP hitsFromTracks(trackHandle, evt, tag); + if (!hitsFromTracks.isValid()) continue; + + for (size_t iTrk = 0; iTrk < trackHandle->size(); ++iTrk) + { + recob::Track const& trk = (*trackHandle)[iTrk]; + auto const& hits = hitsFromTracks.at(iTrk); + auto const& metas = hitsFromTracks.data(iTrk); + for (size_t iH = 0; iH < hits.size(); ++iH) + { + art::Ptr const& hit = hits[iH]; + size_t pt = metas[iH]->Index(); + if (pt == std::numeric_limits::max() || !trk.HasValidPoint(pt)) + continue; + + auto const& wid = hit->WireID(); + int plane = static_cast(wid.Plane); + double thetaZ = fWireReadout->Plane(wid.planeID()).ThetaZ(); + + auto dir = trk.DirectionAtPoint(pt); + if (std::abs(dir.X()) + std::abs(dir.Y()) + std::abs(dir.Z()) < 1e-6) continue; + double recoThxw = ThetaXW(dir.X(), dir.Y(), dir.Z(), thetaZ) * TMath::RadToDeg(); + + // find the truth ThetaXW on this channel within the tick window and same plane + float bestTrue = -1.f, bestDt = fMatchWindowTicks + 1.f; + auto range = truthThxwByChannel.equal_range(hit->Channel()); + for (auto it = range.first; it != range.second; ++it) + { + float tpeak; int tplane; float tthxw; + std::tie(tpeak, tplane, tthxw) = it->second; + if (tplane != plane) continue; + float dt = std::abs(tpeak - static_cast(hit->PeakTime())); + if (dt < bestDt) { bestDt = dt; bestTrue = tthxw; } + } + if (bestTrue < 0.f || bestDt > fMatchWindowTicks) continue; + + fThxwTrueVsReco->Fill(bestTrue, recoThxw); + if (plane >= 0 && plane < 3) fThxwTrueVsReco_p[plane]->Fill(bestTrue, recoThxw); + fThxwResid->Fill(recoThxw - bestTrue); + } + } + } } DEFINE_ART_MODULE(WireModHitDiagnostic) From 6a7c821c1cc14e3502f5de4de79941d1ce478add Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Wed, 1 Jul 2026 09:16:19 -0500 Subject: [PATCH 5/9] Fix dummy hit finder config to match the nominal hit finder. --- .../WireMod/stage1_icarus_wiremod_pops.fcl | 12 ++++++++++++ .../stage1_icarus_wiremod_pops_hitdiag.fcl | 12 ++++++++++++ .../stage1_run2_icarus_overlay_wiremodcorr.fcl | 12 ++++++++++++ ...e1_run2_larcv_icarus_overlay_wiremodcorr.fcl | 12 ++++++++++++ .../WireMod/stage1_wiremod_2d_noscaling.fcl | 12 ++++++++++++ icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl | 12 ++++++++++++ .../WireMod/stage1_wiremod_PT_noscaling.fcl | 17 +++++++++++------ icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl | 17 +++++++++++------ .../WireMod/stage1_wiremod_diagnostic.fcl | 12 ++++++++++++ .../WireMod/stage1_wiremod_float_diagnostic.fcl | 12 ++++++++++++ icaruscode/WireMod/wireMod_stage1.fcl | 12 ++++++++++++ 11 files changed, 130 insertions(+), 12 deletions(-) diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl index a43af0a65..29c758c50 100644 --- a/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl +++ b/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl @@ -27,6 +27,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + physics.producers.WireModTPCEE.AdditiveModification: true physics.producers.WireModTPCEW.AdditiveModification: true physics.producers.WireModTPCWE.AdditiveModification: true diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl index 6ccb100b4..dcc106191 100644 --- a/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl +++ b/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl @@ -24,6 +24,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + physics.producers.WireModTPCEE.AdditiveModification: true physics.producers.WireModTPCEW.AdditiveModification: true physics.producers.WireModTPCWE.AdditiveModification: true diff --git a/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl b/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl index fcb587e12..39f128396 100644 --- a/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl +++ b/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl @@ -27,6 +27,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] diff --git a/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl b/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl index 5ec51d08d..de6512f20 100644 --- a/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl +++ b/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl @@ -27,6 +27,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] diff --git a/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl b/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl index 8c35060fb..a01872f35 100644 --- a/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl +++ b/icaruscode/WireMod/stage1_wiremod_2d_noscaling.fcl @@ -24,6 +24,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + physics.producers.WireModTPCEE.AdditiveModification: true physics.producers.WireModTPCEW.AdditiveModification: true physics.producers.WireModTPCWE.AdditiveModification: true diff --git a/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl b/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl index 6722f54b7..ecd51ae46 100644 --- a/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl +++ b/icaruscode/WireMod/stage1_wiremod_2d_scaled.fcl @@ -24,6 +24,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + physics.producers.WireModTPCEE.AdditiveModification: true physics.producers.WireModTPCEW.AdditiveModification: true physics.producers.WireModTPCWE.AdditiveModification: true diff --git a/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl b/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl index 5a4836801..e6542fe45 100644 --- a/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl +++ b/icaruscode/WireMod/stage1_wiremod_PT_noscaling.fcl @@ -15,12 +15,17 @@ physics.producers: @table::physics.producers } -# "before" hits at the same LongMaxHits=10 setting as gaushitPT2dTPCXX, reading -# the undisturbed wire (not redirected to WMWireToROI) -physics.producers.dummyGaushitPT2dTPCEE: { @table::icarus_stage1_producers.gaushitPT2dTPCEE } -physics.producers.dummyGaushitPT2dTPCEW: { @table::icarus_stage1_producers.gaushitPT2dTPCEW } -physics.producers.dummyGaushitPT2dTPCWE: { @table::icarus_stage1_producers.gaushitPT2dTPCWE } -physics.producers.dummyGaushitPT2dTPCWW: { @table::icarus_stage1_producers.gaushitPT2dTPCWW } +# Keep the "before" hit finder (dummyGaushitPT2dTPC*) configuration identical to the "after" +# finder (gaushitPT2dTPC*): clone the fully-configured gaushitPT2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushitPT2dTPCEE: @local::physics.producers.gaushitPT2dTPCEE +physics.producers.dummyGaushitPT2dTPCEW: @local::physics.producers.gaushitPT2dTPCEW +physics.producers.dummyGaushitPT2dTPCWE: @local::physics.producers.gaushitPT2dTPCWE +physics.producers.dummyGaushitPT2dTPCWW: @local::physics.producers.gaushitPT2dTPCWW +physics.producers.dummyGaushitPT2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushitPT2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushitPT2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushitPT2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" diff --git a/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl b/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl index 824b5a4bd..997b62fab 100644 --- a/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl +++ b/icaruscode/WireMod/stage1_wiremod_PT_scaled.fcl @@ -15,12 +15,17 @@ physics.producers: @table::physics.producers } -# "before" hits at the same LongMaxHits=10 setting as gaushitPT2dTPCXX, reading -# the undisturbed wire (not redirected to WMWireToROI) -physics.producers.dummyGaushitPT2dTPCEE: { @table::icarus_stage1_producers.gaushitPT2dTPCEE } -physics.producers.dummyGaushitPT2dTPCEW: { @table::icarus_stage1_producers.gaushitPT2dTPCEW } -physics.producers.dummyGaushitPT2dTPCWE: { @table::icarus_stage1_producers.gaushitPT2dTPCWE } -physics.producers.dummyGaushitPT2dTPCWW: { @table::icarus_stage1_producers.gaushitPT2dTPCWW } +# Keep the "before" hit finder (dummyGaushitPT2dTPC*) configuration identical to the "after" +# finder (gaushitPT2dTPC*): clone the fully-configured gaushitPT2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushitPT2dTPCEE: @local::physics.producers.gaushitPT2dTPCEE +physics.producers.dummyGaushitPT2dTPCEW: @local::physics.producers.gaushitPT2dTPCEW +physics.producers.dummyGaushitPT2dTPCWE: @local::physics.producers.gaushitPT2dTPCWE +physics.producers.dummyGaushitPT2dTPCWW: @local::physics.producers.gaushitPT2dTPCWW +physics.producers.dummyGaushitPT2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushitPT2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushitPT2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushitPT2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" diff --git a/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl b/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl index 875aeaa53..ee1be4347 100644 --- a/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl +++ b/icaruscode/WireMod/stage1_wiremod_diagnostic.fcl @@ -27,6 +27,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] diff --git a/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl b/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl index 5da6d2a1b..8c3170b4b 100644 --- a/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl +++ b/icaruscode/WireMod/stage1_wiremod_float_diagnostic.fcl @@ -52,6 +52,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WireModTPCEW" physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WireModTPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WireModTPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + # WMWireToROI doesn't run in float mode, redirect to original ChannelROI physics.analyzers.caloskimE.RawDigitproducers: [ "wire2channelroi2d:PHYSCRATEDATATPCEW", "wire2channelroi2d:PHYSCRATEDATATPCEE" ] diff --git a/icaruscode/WireMod/wireMod_stage1.fcl b/icaruscode/WireMod/wireMod_stage1.fcl index 9cac89982..794bd9f24 100644 --- a/icaruscode/WireMod/wireMod_stage1.fcl +++ b/icaruscode/WireMod/wireMod_stage1.fcl @@ -27,6 +27,18 @@ physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDAT physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] From f636b16d029f120c1b817cbdb3c247694a02850c Mon Sep 17 00:00:00 2001 From: Gray Putnam Date: Thu, 2 Jul 2026 14:08:55 -0500 Subject: [PATCH 6/9] Fix thXW calculation. --- icaruscode/WireMod/WireModHitDiagnostic_module.cc | 4 ++-- icaruscode/WireMod/WireModifierXXW_module.cc | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/icaruscode/WireMod/WireModHitDiagnostic_module.cc b/icaruscode/WireMod/WireModHitDiagnostic_module.cc index 34418b620..e6a8e27b9 100644 --- a/icaruscode/WireMod/WireModHitDiagnostic_module.cc +++ b/icaruscode/WireMod/WireModHitDiagnostic_module.cc @@ -39,8 +39,8 @@ namespace wiremod { // with exactly the same convention as the truth angle stored by WireModifierXXW. inline double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) { - const double s = std::sin(planeAngle - util::pi()); - const double c = std::cos(planeAngle - util::pi()); + const double s = std::sin(planeAngle - 0.5 * util::pi()); + const double c = std::cos(planeAngle - 0.5 * util::pi()); const double cosG = std::abs(dydr * s + dzdr * c); return (cosG > 1e-9) ? std::abs(std::atan(dxdr / cosG)) : 0.5 * util::pi(); } diff --git a/icaruscode/WireMod/WireModifierXXW_module.cc b/icaruscode/WireMod/WireModifierXXW_module.cc index a83e1937b..fd5a8d4c2 100644 --- a/icaruscode/WireMod/WireModifierXXW_module.cc +++ b/icaruscode/WireMod/WireModifierXXW_module.cc @@ -1118,7 +1118,7 @@ namespace wiremod (*vec_dirx)[hitVec_idx] = static_cast(tv.dxdr); (*vec_diry)[hitVec_idx] = static_cast(tv.dydr); (*vec_dirz)[hitVec_idx] = static_cast(tv.dzdr); - double cosG = std::abs(tv.dydr * std::sin(thetaZ) + tv.dzdr * std::cos(thetaZ)); + double cosG = std::abs(tv.dydr * std::cos(thetaZ) - tv.dzdr * std::sin(thetaZ)); // pitch: project onto wire-normal (thetaZ - pi/2) float wp = static_cast(plane_obj.WirePitch()); float pval = (cosG > 1e-6) ? wp / static_cast(cosG) : -999.f; (*vec_pitch)[hitVec_idx] = pval; @@ -1516,7 +1516,7 @@ namespace wiremod (*vec_dirx)[hitVec_idx] = static_cast(tv.dxdr); (*vec_diry)[hitVec_idx] = static_cast(tv.dydr); (*vec_dirz)[hitVec_idx] = static_cast(tv.dzdr); - double cosG = std::abs(tv.dydr * std::sin(thetaZ) + tv.dzdr * std::cos(thetaZ)); + double cosG = std::abs(tv.dydr * std::cos(thetaZ) - tv.dzdr * std::sin(thetaZ)); // pitch: project onto wire-normal (thetaZ - pi/2) float wp = static_cast(plane_obj.WirePitch()); float pval = (cosG > 1e-6) ? wp / static_cast(cosG) : -999.f; (*vec_pitch)[hitVec_idx] = pval; From e1863b824129a95bc2eb3166d07acf0a22062b99 Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Wed, 8 Jul 2026 23:41:34 +0200 Subject: [PATCH 7/9] WireModHitDiagnostic: added some additional fields. Moved the TGraph2D interpolation to an option (UseGraph2DInterpolation), default is true. Renamed the fcls with _pops_ to _develop_ (bad practice on my local branch that got picked up) Added a fcl for null-wiremod modification. Added a fcl for the full diagnostic available in ntuples, including the modified ROIs. --- icaruscode/WireMod/WireModHitDiagnostic.fcl | 7 + .../WireMod/WireModHitDiagnostic_module.cc | 498 ++++++++++-------- .../WireMod/WireModifierXXW_closure.fcl | 1 + icaruscode/WireMod/WireModifierXXW_module.cc | 20 +- ....fcl => stage1_icarus_wiremod_develop.fcl} | 6 + ...stage1_icarus_wiremod_develop_hitdiag.fcl} | 21 +- ...e1_icarus_wiremod_develop_hitdiag_null.fcl | 13 + ...ge1_icarus_wiremod_develop_hitdiag_roi.fcl | 107 ++++ 8 files changed, 446 insertions(+), 227 deletions(-) rename icaruscode/WireMod/{stage1_icarus_wiremod_pops.fcl => stage1_icarus_wiremod_develop.fcl} (93%) rename icaruscode/WireMod/{stage1_icarus_wiremod_pops_hitdiag.fcl => stage1_icarus_wiremod_develop_hitdiag.fcl} (82%) create mode 100644 icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_null.fcl create mode 100644 icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_roi.fcl diff --git a/icaruscode/WireMod/WireModHitDiagnostic.fcl b/icaruscode/WireMod/WireModHitDiagnostic.fcl index 7e9e5eb19..604114f4b 100644 --- a/icaruscode/WireMod/WireModHitDiagnostic.fcl +++ b/icaruscode/WireMod/WireModHitDiagnostic.fcl @@ -25,6 +25,13 @@ WireModHitDiag: "WireModTPCWW" ] + TrackLabels: [ + "pandoraTrackGausCryoE", + "pandoraTrackGausCryoW" + ] + + DummyTrackLabels: [] # optional: Pandora on before-hits for PRE-WireMod reco theta + MatchWindowTicks: 10.0 OnlyMC: true } diff --git a/icaruscode/WireMod/WireModHitDiagnostic_module.cc b/icaruscode/WireMod/WireModHitDiagnostic_module.cc index e6a8e27b9..825559c2d 100644 --- a/icaruscode/WireMod/WireModHitDiagnostic_module.cc +++ b/icaruscode/WireMod/WireModHitDiagnostic_module.cc @@ -30,19 +30,75 @@ #include "larcore/Geometry/WireReadout.h" #include "larcorealg/Geometry/WireReadoutGeom.h" #include "larcore/CoreUtils/ServiceUtil.h" -#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" // util::pi namespace wiremod { -// WireMod ThetaXW (angle of the track to the drift/X direction in the wire-projected -// frame), replicated from sys::WireModUtility::ThetaXW so the reco-track angle is computed -// with exactly the same convention as the truth angle stored by WireModifierXXW. +// WireMod ThetaXW replicated from WireModUtility::ThetaXW so reco-track angles +// are computed with exactly the same convention as truth angles from WireModifierXXW. inline double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) { - const double s = std::sin(planeAngle - 0.5 * util::pi()); - const double c = std::cos(planeAngle - 0.5 * util::pi()); + const double s = std::sin(planeAngle - 0.5 * M_PI); + const double c = std::cos(planeAngle - 0.5 * M_PI); const double cosG = std::abs(dydr * s + dzdr * c); - return (cosG > 1e-9) ? std::abs(std::atan(dxdr / cosG)) : 0.5 * util::pi(); + return (cosG > 1e-9) ? std::abs(std::atan(dxdr / cosG)) : 0.5 * M_PI; +} + +// Map type: channel → (peakTime, plane, reco_thxw_deg) +// Used to carry pre-computed reco angles keyed by the hit they belong to. +using RecoThxwMap = std::unordered_multimap< + raw::ChannelID_t, + std::tuple>; + +static RecoThxwMap buildRecoThxwMap(art::Event const& evt, + std::vector const& tags, + geo::WireReadoutGeom const* wr) +{ + RecoThxwMap m; + for (auto const& tag : tags) { + auto h = evt.getHandle>(tag); + if (!h.isValid()) { + mf::LogWarning("WireModHitDiagnostic") + << "Track collection not found: " << tag; continue; + } + art::FindManyP hfm(h, evt, tag); + if (!hfm.isValid()) continue; + for (size_t iT = 0; iT < h->size(); ++iT) { + auto const& trk = (*h)[iT]; + auto const& hits = hfm.at(iT); + auto const& metas = hfm.data(iT); + for (size_t iH = 0; iH < hits.size(); ++iH) { + size_t pt = metas[iH]->Index(); + if (pt == std::numeric_limits::max() || !trk.HasValidPoint(pt)) + continue; + auto dir = trk.DirectionAtPoint(pt); + if (std::abs(dir.X()) + std::abs(dir.Y()) + std::abs(dir.Z()) < 1e-6) continue; + auto const& wid = hits[iH]->WireID(); + double thetaZ = wr->Plane(wid.planeID()).ThetaZ(); + float thxw_deg = static_cast( + ThetaXW(dir.X(), dir.Y(), dir.Z(), thetaZ) * TMath::RadToDeg()); + m.emplace(hits[iH]->Channel(), + std::make_tuple(hits[iH]->PeakTime(), static_cast(wid.Plane), thxw_deg)); + } + } + } + return m; +} + +// Returns -999 if no match within ±matchWindow ticks. +static float lookupRecoThxw(RecoThxwMap const& m, + raw::ChannelID_t ch, float peakTime, int plane, + float matchWindow) +{ + auto range = m.equal_range(ch); + float best = -999.f, bestDt = matchWindow + 1.f; + for (auto it = range.first; it != range.second; ++it) { + float tp; int pl; float thxw; + std::tie(tp, pl, thxw) = it->second; + if (pl != plane) continue; + float dt = std::abs(tp - peakTime); + if (dt < bestDt) { bestDt = dt; best = thxw; } + } + return best; } class WireModHitDiagnostic : public art::EDAnalyzer @@ -52,22 +108,24 @@ class WireModHitDiagnostic : public art::EDAnalyzer void analyze(art::Event const& evt) override; private: - std::vector fBeforeHitLabels; // dummyGaushit2dTPC* - std::vector fAfterHitLabels; // gaushit2dTPC* - std::vector fScaleInfoLabels; // WireModTPC* (produce scale products) - std::vector fTrackLabels; // reco tracks for reco ThetaXW (pandoraTrackGaus*) + std::vector fBeforeHitLabels; // dummyGaushit2dTPC* + std::vector fAfterHitLabels; // gaushit2dTPC* + std::vector fScaleInfoLabels; // WireModTPC* (scale products) + std::vector fTrackLabels; // pandoraTrackGaus* on after-hits + std::vector fDummyTrackLabels; // pandoraTrackGaus* on before-hits (optional) float fMatchWindowTicks; - bool fOnlyMC; // if true, skip overlay hits (is_mc==0) when filling HitMatchTree + bool fOnlyMC; const geo::WireReadoutGeom* fWireReadout = &(art::ServiceHandle()->Get()); - // True (WireMod/MCParticle) vs reco (Pandora track) ThetaXW resolution histograms [deg]. - // Filled for MC hits that are on a reconstructed track: the true angle is the WireModifierXXW - // per-hit thetaXW; the reco angle is ThetaXW of the track direction at that hit. - TH2F* fThxwTrueVsReco = nullptr; // all planes - TH2F* fThxwTrueVsReco_p[3] = {nullptr, nullptr, nullptr}; - TH1F* fThxwResid = nullptr; // reco - true [deg] + // ThetaXW 2D: truth vs reco (after-hit tracks) and truth vs dummy-reco (before-hit tracks) + TH2F* fThxwTrueVsReco = nullptr; + TH2F* fThxwTrueVsReco_p[3] = {nullptr, nullptr, nullptr}; + TH1F* fThxwResid = nullptr; + TH2F* fThxwTrueVsDummyReco = nullptr; + TH2F* fThxwTrueVsDummyReco_p[3] = {nullptr, nullptr, nullptr}; + TH1F* fThxwDummyResid = nullptr; // HitMatchTree: one entry per before-hit TTree* fHitMatchTree; @@ -77,31 +135,22 @@ class WireModHitDiagnostic : public art::EDAnalyzer Int_t tm_before_start_tick, tm_before_end_tick; Float_t tm_before_gof; Int_t tm_before_dof, tm_before_n_multi; - // MC / overlay flag from WireModifierXXW BackTracker - Int_t tm_is_mc; // 1 = MC hit, 0 = overlay data hit - // scale info from WireModifierXXW (-999 if overlay / no EDep match) - Float_t tm_scale_q; - Float_t tm_scale_sigma; - Float_t tm_truth_x; // truth X coordinate used for graph lookup (cm) - Float_t tm_theta_xw; // ThetaXW angle used for graph lookup (degrees) - Float_t tm_truth_dirx; // truth dxdr - Float_t tm_truth_diry; // truth dydr - Float_t tm_truth_dirz; // truth dzdr - Float_t tm_pitch; // wirePitch / cosG [cm] - Float_t tm_dqdx_before; // integral_before / pitch [ADC/cm], no lifetime correction - // SimChannel/BackTracker truth from WireModifierXXW (TrackCaloSkimmer-identical recipe), - // -1 if the before-hit window has no SimChannel ionization (overlay/data hit) - Float_t tm_truth_e; // total true energy under the before-hit [MeV] - Float_t tm_truth_nelec; // total ionization electrons under the before-hit - // match result - Int_t tm_is_matched; - Int_t tm_n_after_matches; + Int_t tm_is_mc; + Float_t tm_scale_q, tm_scale_sigma; + Float_t tm_truth_x; + Float_t tm_theta_xw; // truth ThetaXW from WireMod [deg] + Float_t tm_reco_theta_xw; // reco ThetaXW from after-hit track [deg] + Float_t tm_dummy_reco_theta_xw; // reco ThetaXW from before-hit track [deg] + Float_t tm_truth_dirx, tm_truth_diry, tm_truth_dirz; + Float_t tm_pitch; + Float_t tm_dqdx_before; + Float_t tm_truth_e, tm_truth_nelec; + Int_t tm_is_matched, tm_n_after_matches, tm_n_before_matches; Float_t tm_after_peak_time, tm_after_peak_amp, tm_after_integral, tm_after_rms; Int_t tm_after_start_tick, tm_after_end_tick; Float_t tm_after_gof; - Int_t tm_after_dof; - Float_t tm_ratio_integral; - Float_t tm_ratio_rms; + Int_t tm_after_dof, tm_after_n_multi; + Float_t tm_ratio_integral, tm_ratio_rms; // NewHitTree: one entry per after-hit without a before-hit match TTree* fNewHitTree; @@ -127,67 +176,92 @@ WireModHitDiagnostic::WireModHitDiagnostic(fhicl::ParameterSet const& pset) std::vector{})) fScaleInfoLabels.emplace_back(s); for (auto const& s : pset.get>("TrackLabels", - std::vector{})) + std::vector{})) fTrackLabels.emplace_back(s); + for (auto const& s : pset.get>("DummyTrackLabels", + std::vector{})) + fDummyTrackLabels.emplace_back(s); art::ServiceHandle tfs; - // True-vs-reco ThetaXW resolution histograms (degrees) - const int nb = 90; - const double lo = 0.0, hi = 90.0; + const int nb = 90; const double lo = 0., hi = 90.; + + // POST-WireMod reco vs truth fThxwTrueVsReco = tfs->make("ThxwTrueVsReco", - "True vs reco #theta_{XW} (all planes);true #theta_{XW} [deg];reco #theta_{XW} [deg]", + "Truth vs POST-WireMod reco #theta_{XW} (all planes);" + "truth #theta_{XW} [deg];reco #theta_{XW} [deg]", nb, lo, hi, nb, lo, hi); for (int p = 0; p < 3; ++p) fThxwTrueVsReco_p[p] = tfs->make(Form("ThxwTrueVsReco_p%d", p), - Form("True vs reco #theta_{XW} (plane %d);true #theta_{XW} [deg];reco #theta_{XW} [deg]", p), + Form("Truth vs POST-WireMod reco #theta_{XW} (plane %d);" + "truth #theta_{XW} [deg];reco #theta_{XW} [deg]", p), nb, lo, hi, nb, lo, hi); fThxwResid = tfs->make("ThxwResid", - "reco - true #theta_{XW};#Delta#theta_{XW} [deg];hits", 200, -45.0, 45.0); + "POST-WireMod reco - truth #theta_{XW};" + "#Delta#theta_{XW} [deg];hits", 200, -45., 45.); + + // PRE-WireMod reco vs truth (only filled when DummyTrackLabels are provided) + fThxwTrueVsDummyReco = tfs->make("ThxwTrueVsDummyReco", + "Truth vs PRE-WireMod reco #theta_{XW} (all planes);" + "truth #theta_{XW} [deg];reco #theta_{XW} [deg]", + nb, lo, hi, nb, lo, hi); + for (int p = 0; p < 3; ++p) + fThxwTrueVsDummyReco_p[p] = tfs->make(Form("ThxwTrueVsDummyReco_p%d", p), + Form("Truth vs PRE-WireMod reco #theta_{XW} (plane %d);" + "truth #theta_{XW} [deg];reco #theta_{XW} [deg]", p), + nb, lo, hi, nb, lo, hi); + fThxwDummyResid = tfs->make("ThxwDummyResid", + "PRE-WireMod reco - truth #theta_{XW};" + "#Delta#theta_{XW} [deg];hits", 200, -45., 45.); fHitMatchTree = tfs->make("HitMatchTree", "WireMod before-after hit comparison"); - fHitMatchTree->Branch("run", &tm_run, "run/I"); - fHitMatchTree->Branch("subrun", &tm_subrun, "subrun/I"); - fHitMatchTree->Branch("event", &tm_event, "event/I"); - fHitMatchTree->Branch("channel", &tm_channel, "channel/I"); - fHitMatchTree->Branch("cryo", &tm_cryo, "cryo/I"); - fHitMatchTree->Branch("tpc", &tm_tpc, "tpc/I"); - fHitMatchTree->Branch("plane", &tm_plane, "plane/I"); - fHitMatchTree->Branch("wire", &tm_wire, "wire/I"); - fHitMatchTree->Branch("before_peak_time", &tm_before_peak_time, "before_peak_time/F"); - fHitMatchTree->Branch("before_peak_amp", &tm_before_peak_amp, "before_peak_amp/F"); - fHitMatchTree->Branch("before_integral", &tm_before_integral, "before_integral/F"); - fHitMatchTree->Branch("before_rms", &tm_before_rms, "before_rms/F"); - fHitMatchTree->Branch("before_start_tick", &tm_before_start_tick, "before_start_tick/I"); - fHitMatchTree->Branch("before_end_tick", &tm_before_end_tick, "before_end_tick/I"); - fHitMatchTree->Branch("before_gof", &tm_before_gof, "before_gof/F"); - fHitMatchTree->Branch("before_dof", &tm_before_dof, "before_dof/I"); - fHitMatchTree->Branch("before_n_multi", &tm_before_n_multi, "before_n_multi/I"); - fHitMatchTree->Branch("is_mc", &tm_is_mc, "is_mc/I"); - fHitMatchTree->Branch("scale_q", &tm_scale_q, "scale_q/F"); - fHitMatchTree->Branch("scale_sigma", &tm_scale_sigma, "scale_sigma/F"); - fHitMatchTree->Branch("truth_x", &tm_truth_x, "truth_x/F"); - fHitMatchTree->Branch("theta_xw", &tm_theta_xw, "theta_xw/F"); - fHitMatchTree->Branch("truth_dirx", &tm_truth_dirx, "truth_dirx/F"); - fHitMatchTree->Branch("truth_diry", &tm_truth_diry, "truth_diry/F"); - fHitMatchTree->Branch("truth_dirz", &tm_truth_dirz, "truth_dirz/F"); - fHitMatchTree->Branch("pitch", &tm_pitch, "pitch/F"); - fHitMatchTree->Branch("dqdx_before", &tm_dqdx_before, "dqdx_before/F"); - fHitMatchTree->Branch("truth_e", &tm_truth_e, "truth_e/F"); - fHitMatchTree->Branch("truth_nelec", &tm_truth_nelec, "truth_nelec/F"); - fHitMatchTree->Branch("is_matched", &tm_is_matched, "is_matched/I"); - fHitMatchTree->Branch("n_after_matches", &tm_n_after_matches, "n_after_matches/I"); - fHitMatchTree->Branch("after_peak_time", &tm_after_peak_time, "after_peak_time/F"); - fHitMatchTree->Branch("after_peak_amp", &tm_after_peak_amp, "after_peak_amp/F"); - fHitMatchTree->Branch("after_integral", &tm_after_integral, "after_integral/F"); - fHitMatchTree->Branch("after_rms", &tm_after_rms, "after_rms/F"); - fHitMatchTree->Branch("after_start_tick", &tm_after_start_tick, "after_start_tick/I"); - fHitMatchTree->Branch("after_end_tick", &tm_after_end_tick, "after_end_tick/I"); - fHitMatchTree->Branch("after_gof", &tm_after_gof, "after_gof/F"); - fHitMatchTree->Branch("after_dof", &tm_after_dof, "after_dof/I"); - fHitMatchTree->Branch("ratio_integral", &tm_ratio_integral, "ratio_integral/F"); - fHitMatchTree->Branch("ratio_rms", &tm_ratio_rms, "ratio_rms/F"); + fHitMatchTree->Branch("run", &tm_run, "run/I"); + fHitMatchTree->Branch("subrun", &tm_subrun, "subrun/I"); + fHitMatchTree->Branch("event", &tm_event, "event/I"); + fHitMatchTree->Branch("channel", &tm_channel, "channel/I"); + fHitMatchTree->Branch("cryo", &tm_cryo, "cryo/I"); + fHitMatchTree->Branch("tpc", &tm_tpc, "tpc/I"); + fHitMatchTree->Branch("plane", &tm_plane, "plane/I"); + fHitMatchTree->Branch("wire", &tm_wire, "wire/I"); + fHitMatchTree->Branch("before_peak_time", &tm_before_peak_time, "before_peak_time/F"); + fHitMatchTree->Branch("before_peak_amp", &tm_before_peak_amp, "before_peak_amp/F"); + fHitMatchTree->Branch("before_integral", &tm_before_integral, "before_integral/F"); + fHitMatchTree->Branch("before_rms", &tm_before_rms, "before_rms/F"); + fHitMatchTree->Branch("before_start_tick", &tm_before_start_tick, "before_start_tick/I"); + fHitMatchTree->Branch("before_end_tick", &tm_before_end_tick, "before_end_tick/I"); + fHitMatchTree->Branch("before_gof", &tm_before_gof, "before_gof/F"); + fHitMatchTree->Branch("before_dof", &tm_before_dof, "before_dof/I"); + fHitMatchTree->Branch("before_n_multi", &tm_before_n_multi, "before_n_multi/I"); + fHitMatchTree->Branch("is_mc", &tm_is_mc, "is_mc/I"); + fHitMatchTree->Branch("scale_q", &tm_scale_q, "scale_q/F"); + fHitMatchTree->Branch("scale_sigma", &tm_scale_sigma, "scale_sigma/F"); + fHitMatchTree->Branch("truth_x", &tm_truth_x, "truth_x/F"); + // ThetaXW angles: all three perspectives + fHitMatchTree->Branch("theta_xw", &tm_theta_xw, "theta_xw/F"); + fHitMatchTree->Branch("reco_theta_xw", &tm_reco_theta_xw, "reco_theta_xw/F"); + fHitMatchTree->Branch("dummy_reco_theta_xw", &tm_dummy_reco_theta_xw, "dummy_reco_theta_xw/F"); + fHitMatchTree->Branch("truth_dirx", &tm_truth_dirx, "truth_dirx/F"); + fHitMatchTree->Branch("truth_diry", &tm_truth_diry, "truth_diry/F"); + fHitMatchTree->Branch("truth_dirz", &tm_truth_dirz, "truth_dirz/F"); + fHitMatchTree->Branch("pitch", &tm_pitch, "pitch/F"); + fHitMatchTree->Branch("dqdx_before", &tm_dqdx_before, "dqdx_before/F"); + fHitMatchTree->Branch("truth_e", &tm_truth_e, "truth_e/F"); + fHitMatchTree->Branch("truth_nelec", &tm_truth_nelec, "truth_nelec/F"); + fHitMatchTree->Branch("is_matched", &tm_is_matched, "is_matched/I"); + fHitMatchTree->Branch("n_after_matches", &tm_n_after_matches, "n_after_matches/I"); + fHitMatchTree->Branch("n_before_matches", &tm_n_before_matches, "n_before_matches/I"); + fHitMatchTree->Branch("after_peak_time", &tm_after_peak_time, "after_peak_time/F"); + fHitMatchTree->Branch("after_peak_amp", &tm_after_peak_amp, "after_peak_amp/F"); + fHitMatchTree->Branch("after_integral", &tm_after_integral, "after_integral/F"); + fHitMatchTree->Branch("after_rms", &tm_after_rms, "after_rms/F"); + fHitMatchTree->Branch("after_start_tick", &tm_after_start_tick, "after_start_tick/I"); + fHitMatchTree->Branch("after_end_tick", &tm_after_end_tick, "after_end_tick/I"); + fHitMatchTree->Branch("after_gof", &tm_after_gof, "after_gof/F"); + fHitMatchTree->Branch("after_dof", &tm_after_dof, "after_dof/I"); + fHitMatchTree->Branch("after_n_multi", &tm_after_n_multi, "after_n_multi/I"); + fHitMatchTree->Branch("ratio_integral", &tm_ratio_integral, "ratio_integral/F"); + fHitMatchTree->Branch("ratio_rms", &tm_ratio_rms, "ratio_rms/F"); fNewHitTree = tfs->make("NewHitTree", "After-mod hits with no before-mod counterpart"); @@ -224,46 +298,37 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) tm_subrun = tn_subrun = static_cast(evt.subRun()); tm_event = tn_event = static_cast(evt.event()); - // Gather all after-hits once (searched by channel for every before-hit) + RecoThxwMap recoThxwByChannel = buildRecoThxwMap(evt, fTrackLabels, fWireReadout); + RecoThxwMap dummyRecoThxwByChannel = buildRecoThxwMap(evt, fDummyTrackLabels, fWireReadout); + std::vector afterHits; - for (auto const& tag : fAfterHitLabels) - { + for (auto const& tag : fAfterHitLabels) { auto h = evt.getHandle>(tag); - if (!h.isValid()) - { + if (!h.isValid()) { mf::LogWarning("WireModHitDiagnostic") - << "After-hit collection not found: " << tag; - continue; + << "After-hit collection not found: " << tag; continue; } - for (auto const& hit : *h) - afterHits.push_back(&hit); + for (auto const& hit : *h) afterHits.push_back(&hit); } - // per-channel truth ThetaXW (from WireMod, on before-hits): channel -> (peaktime, plane, true_thxw_deg) - std::unordered_multimap> truthThxwByChannel; - - // Build channel -> [after-hit indices] index std::unordered_multimap afterByChannel; afterByChannel.reserve(afterHits.size() * 2); for (size_t i = 0; i < afterHits.size(); ++i) afterByChannel.emplace(afterHits[i]->Channel(), i); - std::vector afterMatched(afterHits.size(), false); + std::vector afterNBeforeMatches(afterHits.size(), 0); + RecoThxwMap truthThxwByChannel; - // Process before-hit collections one at a time to keep the hit index - // aligned with the scale-info parallel vectors from WireModifierXXW. for (size_t j = 0; j < fBeforeHitLabels.size(); ++j) { auto bHandle = evt.getHandle>(fBeforeHitLabels[j]); - if (!bHandle.isValid()) - { + if (!bHandle.isValid()) { mf::LogWarning("WireModHitDiagnostic") - << "Before-hit collection not found: " << fBeforeHitLabels[j]; - continue; + << "Before-hit collection not found: " << fBeforeHitLabels[j]; continue; } auto const& localHits = *bHandle; - // Try to load scale-info and isMC parallel vectors for this TPC + // scale-info parallel vectors for this TPC (indexed 1:1 with localHits) std::vector const* pScaleQ = nullptr; std::vector const* pScaleSigma = nullptr; std::vector const* pTruthX = nullptr; @@ -277,8 +342,7 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) std::vector const* pTruthE = nullptr; std::vector const* pTruthNelec = nullptr; - if (j < fScaleInfoLabels.size()) - { + if (j < fScaleInfoLabels.size()) { auto const& siLabel = fScaleInfoLabels[j]; auto sqH = evt.getHandle>(art::InputTag(siLabel.label(), "scaleQ")); auto ssH = evt.getHandle>(art::InputTag(siLabel.label(), "scaleSigma")); @@ -306,6 +370,24 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) if (tnH.isValid()) pTruthNelec = &(*tnH); } + // pre-pass: count before-hits per after-hit so n_before_matches is ready for tree fill. + // Window on PeakTime (not StartTick/EndTick) to avoid multi-Gaussian group mismatches. + for (size_t i = 0; i < localHits.size(); ++i) { + auto const& bh = localHits[i]; + if (fOnlyMC) { + int isMC = (pIsMC && i < pIsMC->size()) ? (*pIsMC)[i] : -1; + if (isMC != 1) continue; + } + float winLo = bh.PeakTime() - fMatchWindowTicks; + float winHi = bh.PeakTime() + fMatchWindowTicks; + auto range = afterByChannel.equal_range(bh.Channel()); + for (auto it = range.first; it != range.second; ++it) { + if (afterHits[it->second]->PeakTime() >= winLo && + afterHits[it->second]->PeakTime() <= winHi) + ++afterNBeforeMatches[it->second]; + } + } + for (size_t i = 0; i < localHits.size(); ++i) { auto const& bh = localHits[i]; @@ -326,40 +408,38 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) tm_before_dof = bh.DegreesOfFreedom(); tm_before_n_multi = bh.Multiplicity(); - // MC flag and scale info (parallel vectors indexed by hit position in this collection) - tm_is_mc = (pIsMC && i < pIsMC->size()) ? (*pIsMC)[i] : -1; - tm_scale_q = (pScaleQ && i < pScaleQ->size()) ? (*pScaleQ)[i] : -999.f; - tm_scale_sigma = (pScaleSigma && i < pScaleSigma->size()) ? (*pScaleSigma)[i] : -999.f; - tm_truth_x = (pTruthX && i < pTruthX->size()) ? (*pTruthX)[i] : -999.f; - tm_theta_xw = (pThetaXW && i < pThetaXW->size()) ? (*pThetaXW)[i] : -999.f; - tm_truth_dirx = (pDirX && i < pDirX->size()) ? (*pDirX)[i] : -999.f; - tm_truth_diry = (pDirY && i < pDirY->size()) ? (*pDirY)[i] : -999.f; - tm_truth_dirz = (pDirZ && i < pDirZ->size()) ? (*pDirZ)[i] : -999.f; - tm_pitch = (pPitch && i < pPitch->size()) ? (*pPitch)[i] : -999.f; - tm_dqdx_before = (pDQdX && i < pDQdX->size()) ? (*pDQdX)[i] : -999.f; - tm_truth_e = (pTruthE && i < pTruthE->size()) ? (*pTruthE)[i] : -1.f; - tm_truth_nelec = (pTruthNelec && i < pTruthNelec->size()) ? (*pTruthNelec)[i] : -1.f; - - // Record the truth ThetaXW for this (MC) hit so the reco-track loop below can pair it - // with the reconstructed track angle on the same channel/tick. + tm_is_mc = (pIsMC && i < pIsMC->size()) ? (*pIsMC)[i] : -1; + tm_scale_q = (pScaleQ && i < pScaleQ->size()) ? (*pScaleQ)[i] : -999.f; + tm_scale_sigma = (pScaleSigma && i < pScaleSigma->size()) ? (*pScaleSigma)[i] : -999.f; + tm_truth_x = (pTruthX && i < pTruthX->size()) ? (*pTruthX)[i] : -999.f; + tm_theta_xw = (pThetaXW && i < pThetaXW->size()) ? (*pThetaXW)[i] : -999.f; + tm_truth_dirx = (pDirX && i < pDirX->size()) ? (*pDirX)[i] : -999.f; + tm_truth_diry = (pDirY && i < pDirY->size()) ? (*pDirY)[i] : -999.f; + tm_truth_dirz = (pDirZ && i < pDirZ->size()) ? (*pDirZ)[i] : -999.f; + tm_pitch = (pPitch && i < pPitch->size()) ? (*pPitch)[i] : -999.f; + tm_dqdx_before = (pDQdX && i < pDQdX->size()) ? (*pDQdX)[i] : -999.f; + tm_truth_e = (pTruthE && i < pTruthE->size()) ? (*pTruthE)[i] : -1.f; + tm_truth_nelec = (pTruthNelec && i < pTruthNelec->size()) ? (*pTruthNelec)[i] : -1.f; + + tm_dummy_reco_theta_xw = lookupRecoThxw(dummyRecoThxwByChannel, + bh.Channel(), bh.PeakTime(), tm_plane, + fMatchWindowTicks); + if (tm_theta_xw > -900.f) truthThxwByChannel.emplace(bh.Channel(), - std::make_tuple(tm_before_peak_time, tm_plane, tm_theta_xw)); + std::make_tuple(bh.PeakTime(), tm_plane, tm_theta_xw)); - // if OnlyMC is set, skip overlay hits if (fOnlyMC && tm_is_mc != 1) continue; - // Search for matching after-hit on the same channel - float winLo = static_cast(bh.StartTick()) - fMatchWindowTicks; - float winHi = static_cast(bh.EndTick()) + fMatchWindowTicks; + float winLo = bh.PeakTime() - fMatchWindowTicks; + float winHi = bh.PeakTime() + fMatchWindowTicks; size_t bestIdx = std::numeric_limits::max(); float bestDist = std::numeric_limits::max(); int nMatch = 0; auto range = afterByChannel.equal_range(bh.Channel()); - for (auto it = range.first; it != range.second; ++it) - { + for (auto it = range.first; it != range.second; ++it) { float aPeak = afterHits[it->second]->PeakTime(); if (aPeak < winLo || aPeak > winHi) continue; ++nMatch; @@ -367,18 +447,13 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) if (dist < bestDist) { bestDist = dist; bestIdx = it->second; } } - tm_n_after_matches = nMatch; - tm_is_matched = (bestIdx != std::numeric_limits::max()) ? 1 : 0; - - if (tm_is_matched) - { - for (auto it = range.first; it != range.second; ++it) - { - float aPeak = afterHits[it->second]->PeakTime(); - if (aPeak >= winLo && aPeak <= winHi) - afterMatched[it->second] = true; - } - auto const* ah = afterHits[bestIdx]; + tm_n_after_matches = nMatch; + tm_is_matched = (bestIdx != std::numeric_limits::max()) ? 1 : 0; + tm_n_before_matches = (bestIdx != std::numeric_limits::max()) + ? afterNBeforeMatches[bestIdx] : 0; + + if (tm_is_matched) { + auto const* ah = afterHits[bestIdx]; tm_after_peak_time = ah->PeakTime(); tm_after_peak_amp = ah->PeakAmplitude(); tm_after_integral = ah->Integral(); @@ -387,34 +462,31 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) tm_after_end_tick = ah->EndTick(); tm_after_gof = ah->GoodnessOfFit(); tm_after_dof = ah->DegreesOfFreedom(); + tm_after_n_multi = ah->Multiplicity(); tm_ratio_integral = (bh.Integral() > 0.f) ? ah->Integral() / bh.Integral() : -999.f; tm_ratio_rms = (bh.RMS() > 0.f) ? ah->RMS() / bh.RMS() : -999.f; - } - else - { - tm_after_peak_time = -999.f; - tm_after_peak_amp = -999.f; - tm_after_integral = -999.f; - tm_after_rms = -999.f; - tm_after_start_tick = -999; - tm_after_end_tick = -999; - tm_after_gof = -999.f; - tm_after_dof = -999; - tm_ratio_integral = -999.f; - tm_ratio_rms = -999.f; + + tm_reco_theta_xw = lookupRecoThxw(recoThxwByChannel, + ah->Channel(), ah->PeakTime(), tm_plane, + fMatchWindowTicks); + } else { + tm_after_peak_time = -999.f; tm_after_peak_amp = -999.f; + tm_after_integral = -999.f; tm_after_rms = -999.f; + tm_after_start_tick = -999; tm_after_end_tick = -999; + tm_after_gof = -999.f; tm_after_dof = -999; + tm_after_n_multi = -999; + tm_ratio_integral = -999.f; tm_ratio_rms = -999.f; + tm_reco_theta_xw = -999.f; } fHitMatchTree->Fill(); } } - // NewHitTree: after-hits not claimed by any before-hit - for (size_t i = 0; i < afterHits.size(); ++i) - { - if (afterMatched[i]) continue; - + for (size_t i = 0; i < afterHits.size(); ++i) { + if (afterNBeforeMatches[i] > 0) continue; auto const* ah = afterHits[i]; auto const& wid = ah->WireID(); tn_channel = static_cast(ah->Channel()); @@ -433,58 +505,62 @@ void WireModHitDiagnostic::analyze(art::Event const& evt) fNewHitTree->Fill(); } - // --------------------------------------------------------------------------- - // True-vs-reco ThetaXW: for every hit on a reconstructed track, compute the reco - // ThetaXW from the track direction at that hit, pair it with the WireMod truth - // ThetaXW on the same channel/tick, and fill the resolution histograms. - // --------------------------------------------------------------------------- - for (auto const& tag : fTrackLabels) - { - auto trackHandle = evt.getHandle>(tag); - if (!trackHandle.isValid()) - { - mf::LogWarning("WireModHitDiagnostic") << "Track collection not found: " << tag; - continue; + for (auto const& tag : fTrackLabels) { + auto h = evt.getHandle>(tag); + if (!h.isValid()) continue; + art::FindManyP hfm(h, evt, tag); + if (!hfm.isValid()) continue; + for (size_t iT = 0; iT < h->size(); ++iT) { + auto const& trk = (*h)[iT]; + auto const& hits = hfm.at(iT); + auto const& metas = hfm.data(iT); + for (size_t iH = 0; iH < hits.size(); ++iH) { + size_t pt = metas[iH]->Index(); + if (pt == std::numeric_limits::max() || !trk.HasValidPoint(pt)) continue; + auto dir = trk.DirectionAtPoint(pt); + if (std::abs(dir.X()) + std::abs(dir.Y()) + std::abs(dir.Z()) < 1e-6) continue; + auto const& wid = hits[iH]->WireID(); + int plane = static_cast(wid.Plane); + double thetaZ = fWireReadout->Plane(wid.planeID()).ThetaZ(); + float recoThxw = static_cast( + ThetaXW(dir.X(), dir.Y(), dir.Z(), thetaZ) * TMath::RadToDeg()); + float trueThxw = lookupRecoThxw(truthThxwByChannel, + hits[iH]->Channel(), hits[iH]->PeakTime(), plane, + fMatchWindowTicks); + if (trueThxw < 0.f) continue; + fThxwTrueVsReco->Fill(trueThxw, recoThxw); + if (plane >= 0 && plane < 3) fThxwTrueVsReco_p[plane]->Fill(trueThxw, recoThxw); + fThxwResid->Fill(recoThxw - trueThxw); + } } - art::FindManyP hitsFromTracks(trackHandle, evt, tag); - if (!hitsFromTracks.isValid()) continue; + } - for (size_t iTrk = 0; iTrk < trackHandle->size(); ++iTrk) - { - recob::Track const& trk = (*trackHandle)[iTrk]; - auto const& hits = hitsFromTracks.at(iTrk); - auto const& metas = hitsFromTracks.data(iTrk); - for (size_t iH = 0; iH < hits.size(); ++iH) - { - art::Ptr const& hit = hits[iH]; + for (auto const& tag : fDummyTrackLabels) { + auto h = evt.getHandle>(tag); + if (!h.isValid()) continue; + art::FindManyP hfm(h, evt, tag); + if (!hfm.isValid()) continue; + for (size_t iT = 0; iT < h->size(); ++iT) { + auto const& trk = (*h)[iT]; + auto const& hits = hfm.at(iT); + auto const& metas = hfm.data(iT); + for (size_t iH = 0; iH < hits.size(); ++iH) { size_t pt = metas[iH]->Index(); - if (pt == std::numeric_limits::max() || !trk.HasValidPoint(pt)) - continue; - - auto const& wid = hit->WireID(); - int plane = static_cast(wid.Plane); - double thetaZ = fWireReadout->Plane(wid.planeID()).ThetaZ(); - + if (pt == std::numeric_limits::max() || !trk.HasValidPoint(pt)) continue; auto dir = trk.DirectionAtPoint(pt); if (std::abs(dir.X()) + std::abs(dir.Y()) + std::abs(dir.Z()) < 1e-6) continue; - double recoThxw = ThetaXW(dir.X(), dir.Y(), dir.Z(), thetaZ) * TMath::RadToDeg(); - - // find the truth ThetaXW on this channel within the tick window and same plane - float bestTrue = -1.f, bestDt = fMatchWindowTicks + 1.f; - auto range = truthThxwByChannel.equal_range(hit->Channel()); - for (auto it = range.first; it != range.second; ++it) - { - float tpeak; int tplane; float tthxw; - std::tie(tpeak, tplane, tthxw) = it->second; - if (tplane != plane) continue; - float dt = std::abs(tpeak - static_cast(hit->PeakTime())); - if (dt < bestDt) { bestDt = dt; bestTrue = tthxw; } - } - if (bestTrue < 0.f || bestDt > fMatchWindowTicks) continue; - - fThxwTrueVsReco->Fill(bestTrue, recoThxw); - if (plane >= 0 && plane < 3) fThxwTrueVsReco_p[plane]->Fill(bestTrue, recoThxw); - fThxwResid->Fill(recoThxw - bestTrue); + auto const& wid = hits[iH]->WireID(); + int plane = static_cast(wid.Plane); + double thetaZ = fWireReadout->Plane(wid.planeID()).ThetaZ(); + float dummyRecoThxw = static_cast( + ThetaXW(dir.X(), dir.Y(), dir.Z(), thetaZ) * TMath::RadToDeg()); + float trueThxw = lookupRecoThxw(truthThxwByChannel, + hits[iH]->Channel(), hits[iH]->PeakTime(), plane, + fMatchWindowTicks); + if (trueThxw < 0.f) continue; + fThxwTrueVsDummyReco->Fill(trueThxw, dummyRecoThxw); + if (plane >= 0 && plane < 3) fThxwTrueVsDummyReco_p[plane]->Fill(trueThxw, dummyRecoThxw); + fThxwDummyResid->Fill(dummyRecoThxw - trueThxw); } } } diff --git a/icaruscode/WireMod/WireModifierXXW_closure.fcl b/icaruscode/WireMod/WireModifierXXW_closure.fcl index fba253031..ed1dff124 100644 --- a/icaruscode/WireMod/WireModifierXXW_closure.fcl +++ b/icaruscode/WireMod/WireModifierXXW_closure.fcl @@ -48,6 +48,7 @@ WireMod: SetNullScaleIntegral: false # true = force r_Q=1 SetNullScaleWidth: false # true = force r_sigma=1 + UseGraph2DInterpolation: true # false = nearest bin center instead of Delaunay interpolation } WireModTPCEE: @local::WireMod diff --git a/icaruscode/WireMod/WireModifierXXW_module.cc b/icaruscode/WireMod/WireModifierXXW_module.cc index fd5a8d4c2..200241a06 100644 --- a/icaruscode/WireMod/WireModifierXXW_module.cc +++ b/icaruscode/WireMod/WireModifierXXW_module.cc @@ -81,9 +81,10 @@ namespace wiremod bool fSavePerHitData; // write per-hit scaleQ/scaleSigma/truthX/thetaXW/isMC art products? bool fInRads; // is the TGraph2D angle axis in radians? bool fXAbs; // is the TGraph2D x an absolute value? - bool fAdditiveModification; // additive (true) vs multiplicative (false) ROI modification - bool fSetNullScaleIntegral; // if true, force r_Q=1 for all sub-ROIs (disable integral scaling) - bool fSetNullScaleWidth; // if true, force r_sigma=1 for all sub-ROIs (disable width scaling) + bool fAdditiveModification; // additive (true) vs multiplicative (false) ROI modification + bool fSetNullScaleIntegral; // if true, force r_Q=1 for all sub-ROIs (disable integral scaling) + bool fSetNullScaleWidth; // if true, force r_sigma=1 for all sub-ROIs (disable width scaling) + bool fUseGraph2DInterpolation; // true = Delaunay interp (default), false = nearest bin center bool fUseChannelROIMode; // true = read ChannelROI directly (float precision), false = Wire (nominal) art::InputTag fChannelLabel; // ChannelROI input label, used only when fUseChannelROIMode is true double fOffset; // ad hoc offset @@ -224,9 +225,10 @@ namespace wiremod fXAbs = pset.get("XAbs", false); // additive vs multiplicative ROI modification - fAdditiveModification = pset.get("AdditiveModification", false); - fSetNullScaleIntegral = pset.get("SetNullScaleIntegral", false); - fSetNullScaleWidth = pset.get("SetNullScaleWidth", false); + fAdditiveModification = pset.get("AdditiveModification", false); + fSetNullScaleIntegral = pset.get("SetNullScaleIntegral", false); + fSetNullScaleWidth = pset.get("SetNullScaleWidth", false); + fUseGraph2DInterpolation = pset.get("UseGraph2DInterpolation", true); // ChannelROI direct mode (float precision, bypasses Wire intermediate) fUseChannelROIMode = pset.get("UseChannelROIMode", false); @@ -705,7 +707,8 @@ namespace wiremod wmUtil.graph2Ds_Sigma_XXW = fGraph_sigma_XXW; } - wmUtil.additiveModification = fAdditiveModification; + wmUtil.additiveModification = fAdditiveModification; + wmUtil.useGraph2DInterpolation = fUseGraph2DInterpolation; // add some debugging here mf::LogVerbatim("WireModifierXXW") @@ -718,7 +721,8 @@ namespace wiremod << " applyXZAngleScale: " << wmUtil.applyXZAngleScale << '\n' << " applyYZAngleScale: " << wmUtil.applyYZAngleScale << '\n' << " applydEdXScale: " << wmUtil.applydEdXScale << '\n' - << " additiveModification: " << wmUtil.additiveModification << '\n' + << " additiveModification: " << wmUtil.additiveModification << '\n' + << " useGraph2DInterpolation: " << wmUtil.useGraph2DInterpolation << '\n' << " setNullScaleIntegral: " << fSetNullScaleIntegral << '\n' << " setNullScaleWidth: " << fSetNullScaleWidth << '\n' << " readoutWindowTicks: " << wmUtil.readoutWindowTicks << '\n' diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_develop.fcl similarity index 93% rename from icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl rename to icaruscode/WireMod/stage1_icarus_wiremod_develop.fcl index 29c758c50..7e43eb6d1 100644 --- a/icaruscode/WireMod/stage1_icarus_wiremod_pops.fcl +++ b/icaruscode/WireMod/stage1_icarus_wiremod_develop.fcl @@ -53,6 +53,12 @@ physics.producers.WireModTPCWE.SetNullScaleWidth: false physics.producers.WireModTPCWW.SetNullScaleIntegral: false physics.producers.WireModTPCWW.SetNullScaleWidth: false +physics.producers.WireModTPCEE.UseGraph2DInterpolation: false +physics.producers.WireModTPCEW.UseGraph2DInterpolation: false +physics.producers.WireModTPCWE.UseGraph2DInterpolation: false +physics.producers.WireModTPCWW.UseGraph2DInterpolation: false + + # allRatio_Coarse_XW.root has X/theta axes, not Y/Z: disable YZ correction physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag.fcl similarity index 82% rename from icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl rename to icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag.fcl index dcc106191..2e67a701f 100644 --- a/icaruscode/WireMod/stage1_icarus_wiremod_pops_hitdiag.fcl +++ b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag.fcl @@ -41,14 +41,19 @@ physics.producers.WireModTPCEW.AdditiveModification: true physics.producers.WireModTPCWE.AdditiveModification: true physics.producers.WireModTPCWW.AdditiveModification: true -physics.producers.WireModTPCEE.SetNullScaleIntegral: true -physics.producers.WireModTPCEE.SetNullScaleWidth: true -physics.producers.WireModTPCEW.SetNullScaleIntegral: true -physics.producers.WireModTPCEW.SetNullScaleWidth: true -physics.producers.WireModTPCWE.SetNullScaleIntegral: true -physics.producers.WireModTPCWE.SetNullScaleWidth: true -physics.producers.WireModTPCWW.SetNullScaleIntegral: true -physics.producers.WireModTPCWW.SetNullScaleWidth: true +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.UseGraph2DInterpolation: false +physics.producers.WireModTPCEW.UseGraph2DInterpolation: false +physics.producers.WireModTPCWE.UseGraph2DInterpolation: false +physics.producers.WireModTPCWW.UseGraph2DInterpolation: false # allRatio_Coarse_XW.root has (X, theta) axes, not (Y, Z): disable YZ correction physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_null.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_null.fcl new file mode 100644 index 000000000..3260279ab --- /dev/null +++ b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_null.fcl @@ -0,0 +1,13 @@ +# Null-scale baseline: identical to stage1_icarus_wiremod_develop_hitdiag.fcl +# but with r_Q = r_sigma = 1 forced everywhere (no WireMod correction applied). + +#include "stage1_icarus_wiremod_develop_hitdiag.fcl" + +physics.producers.WireModTPCEE.SetNullScaleIntegral: true +physics.producers.WireModTPCEE.SetNullScaleWidth: true +physics.producers.WireModTPCEW.SetNullScaleIntegral: true +physics.producers.WireModTPCEW.SetNullScaleWidth: true +physics.producers.WireModTPCWE.SetNullScaleIntegral: true +physics.producers.WireModTPCWE.SetNullScaleWidth: true +physics.producers.WireModTPCWW.SetNullScaleIntegral: true +physics.producers.WireModTPCWW.SetNullScaleWidth: true diff --git a/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_roi.fcl b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_roi.fcl new file mode 100644 index 000000000..333fa6762 --- /dev/null +++ b/icaruscode/WireMod/stage1_icarus_wiremod_develop_hitdiag_roi.fcl @@ -0,0 +1,107 @@ +#include "WireModifierXXW_closure.fcl" +#include "WireModHitDiagnostic.fcl" +#include "stage1_run2_icarus_overlay.fcl" + +physics.reco: +[ + @sequence::WireModSequence, + @sequence::physics.reco +] + +physics.producers: +{ + @table::WireModTable + @table::physics.producers +} + +physics.producers.gaushit2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushit2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushit2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushit2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +physics.producers.gaushitPT2dTPCEE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEE" +physics.producers.gaushitPT2dTPCEW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCEW" +physics.producers.gaushitPT2dTPCWE.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWE" +physics.producers.gaushitPT2dTPCWW.CalDataModuleLabel: "WMWireToROI:PHYSCRATEDATATPCWW" + +# Keep the "before" hit finder (dummyGaushit2dTPC*) configuration identical to the "after" +# finder (gaushit2dTPC*): clone the fully-configured gaushit2d and reset ONLY the input ROIs +# to the un-modified wire2channelroi2d collection. +physics.producers.dummyGaushit2dTPCEE: @local::physics.producers.gaushit2dTPCEE +physics.producers.dummyGaushit2dTPCEW: @local::physics.producers.gaushit2dTPCEW +physics.producers.dummyGaushit2dTPCWE: @local::physics.producers.gaushit2dTPCWE +physics.producers.dummyGaushit2dTPCWW: @local::physics.producers.gaushit2dTPCWW +physics.producers.dummyGaushit2dTPCEE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEE" +physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCEW" +physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" +physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" + +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.UseGraph2DInterpolation: false +physics.producers.WireModTPCEW.UseGraph2DInterpolation: false +physics.producers.WireModTPCWE.UseGraph2DInterpolation: false +physics.producers.WireModTPCWW.UseGraph2DInterpolation: false + +# allRatio_Coarse_XW.root has (X, theta) axes, not (Y, Z): disable YZ correction +physics.producers.WireModTPCEE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCEW.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWE.RatioFileName_YZ: "NOFILE" +physics.producers.WireModTPCWW.RatioFileName_YZ: "NOFILE" + +# SavePerHitData: fills the per-hit WireModTPCxx/ROITree used by the notebook. +# SaveROITree: saves the full ROI waveform (roi_before) alongside the hit parameters. +# SaveTickHistograms / SaveByChannel: kept false to limit output file size. +physics.producers.WireModTPCEE.SavePerHitData: true +physics.producers.WireModTPCEW.SavePerHitData: true +physics.producers.WireModTPCWE.SavePerHitData: true +physics.producers.WireModTPCWW.SavePerHitData: true + +physics.producers.WireModTPCEE.SaveROITree: true +physics.producers.WireModTPCEW.SaveROITree: true +physics.producers.WireModTPCWE.SaveROITree: true +physics.producers.WireModTPCWW.SaveROITree: true + +physics.producers.WireModTPCEE.SaveTickHistograms: false +physics.producers.WireModTPCEW.SaveTickHistograms: false +physics.producers.WireModTPCWE.SaveTickHistograms: false +physics.producers.WireModTPCWW.SaveTickHistograms: false + +physics.producers.WireModTPCEE.SaveByChannel: false +physics.producers.WireModTPCEW.SaveByChannel: false +physics.producers.WireModTPCWE.SaveByChannel: false +physics.producers.WireModTPCWW.SaveByChannel: false + +physics.analyzers.WireModHitDiag: @local::WireModHitDiag + +physics.analysisPath: [ WireModHitDiag ] +physics.end_paths: [ outana, stream1, analysisPath ] + +physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", + "WMWireToROI:PHYSCRATEDATATPCEE" ] +physics.analyzers.caloskimW.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCWW", + "WMWireToROI:PHYSCRATEDATATPCWE" ] + +outputs.out1.outputCommands: +[ + @sequence::outputs.out1.outputCommands, + "drop *_dummyGaushit2dTPC*_*_*", + "drop *_WMROI2Wire_*_*", + "drop *_WireModXXWTPC*_*_*" +] + +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "merge" + +services.message.destinations.STDCOUT.threshold: "WARNING" From 13291e8393a74b1cb7a72628354eae97ef8192e5 Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Thu, 9 Jul 2026 00:16:39 +0200 Subject: [PATCH 8/9] Now the fixes to achieve closure are propagated to the production fcls --- ...stage1_run2_icarus_overlay_wiremodcorr.fcl | 19 +++++++++++++++++++ ..._run2_larcv_icarus_overlay_wiremodcorr.fcl | 19 +++++++++++++++++++ 2 files changed, 38 insertions(+) diff --git a/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl b/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl index 39f128396..3695afc0b 100644 --- a/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl +++ b/icaruscode/WireMod/stage1_run2_icarus_overlay_wiremodcorr.fcl @@ -39,6 +39,25 @@ physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHY physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.UseGraph2DInterpolation: false +physics.producers.WireModTPCEW.UseGraph2DInterpolation: false +physics.producers.WireModTPCWE.UseGraph2DInterpolation: false +physics.producers.WireModTPCWW.UseGraph2DInterpolation: false + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] diff --git a/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl b/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl index de6512f20..a773fdce4 100644 --- a/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl +++ b/icaruscode/WireMod/stage1_run2_larcv_icarus_overlay_wiremodcorr.fcl @@ -39,6 +39,25 @@ physics.producers.dummyGaushit2dTPCEW.CalDataModuleLabel: "wire2channelroi2d:PHY physics.producers.dummyGaushit2dTPCWE.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWE" physics.producers.dummyGaushit2dTPCWW.CalDataModuleLabel: "wire2channelroi2d:PHYSCRATEDATATPCWW" +physics.producers.WireModTPCEE.AdditiveModification: true +physics.producers.WireModTPCEW.AdditiveModification: true +physics.producers.WireModTPCWE.AdditiveModification: true +physics.producers.WireModTPCWW.AdditiveModification: true + +physics.producers.WireModTPCEE.SetNullScaleIntegral: false +physics.producers.WireModTPCEE.SetNullScaleWidth: false +physics.producers.WireModTPCEW.SetNullScaleIntegral: false +physics.producers.WireModTPCEW.SetNullScaleWidth: false +physics.producers.WireModTPCWE.SetNullScaleIntegral: false +physics.producers.WireModTPCWE.SetNullScaleWidth: false +physics.producers.WireModTPCWW.SetNullScaleIntegral: false +physics.producers.WireModTPCWW.SetNullScaleWidth: false + +physics.producers.WireModTPCEE.UseGraph2DInterpolation: false +physics.producers.WireModTPCEW.UseGraph2DInterpolation: false +physics.producers.WireModTPCWE.UseGraph2DInterpolation: false +physics.producers.WireModTPCWW.UseGraph2DInterpolation: false + # tell the caloskimmer to use the ROIs from WireMod physics.analyzers.caloskimE.RawDigitproducers: [ "WMWireToROI:PHYSCRATEDATATPCEW", "WMWireToROI:PHYSCRATEDATATPCEE" ] From 092ae11465bfe46adb8af8dfd4064149be8f1cac Mon Sep 17 00:00:00 2001 From: "poppi@baltig.infn.it" Date: Thu, 9 Jul 2026 21:31:27 +0200 Subject: [PATCH 9/9] this fixes corrects an issue with loss of precision in conversion fron channelROI to recob::Wire. There is a different approach which should be used in the future and can be tested with UseChannelROIMode: true, which Harry introduced to skip these conversions. UseChannelROIMode is ideal, but it was tested with a small sample, need higher stats check to verify that closure is achieved also in this way. --- icaruscode/WireMod/WireModifierXXW.fcl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/icaruscode/WireMod/WireModifierXXW.fcl b/icaruscode/WireMod/WireModifierXXW.fcl index bf53718cc..3f656a6bb 100644 --- a/icaruscode/WireMod/WireModifierXXW.fcl +++ b/icaruscode/WireMod/WireModifierXXW.fcl @@ -24,7 +24,8 @@ WMROI2Wire: "PHYSCRATEDATATPCEW", "PHYSCRATEDATATPCWE", "PHYSCRATEDATATPCWW" ] - DiagnosticOutput: false + DiagnosticOutput: false + UseFloatPrecision: true } WireMod: