From cee3459cdfa17e6a61923b95ad797217ded6f7f4 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 10:34:17 -0600 Subject: [PATCH 01/31] initial development --- DataProducts/inc/SurfaceId.hh | 1 + DataProducts/src/SurfaceId.cc | 10 +- KinKalGeom/inc/SurfaceMap.hh | 1 + TrkReco/fcl/prolog.fcl | 8 ++ .../src/InterdetectorTimeBootstrap_module.cc | 129 ++++++++++++++++++ 5 files changed, 148 insertions(+), 1 deletion(-) create mode 100644 TrkReco/src/InterdetectorTimeBootstrap_module.cc diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index d70c661bb8..b9fbde90f0 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -21,6 +21,7 @@ namespace mu2e { OPA=95, TSDA, // Absorbers in the DS ST_Front=100,ST_Back, ST_Inner, ST_Outer, ST_Foils, ST_Wires, // stopping target bounding surfaces and components TCRV=200 // CRV test planes + CalD0_Front=300, CalD0_Back, CalD1_Front, CalD1_Back, CalD0_Inner, CalD0_Outer, CalD1_Inner, CalD1_Outer //calo VD }; static std::string const& typeName(); diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 4d454e0922..a4dc3363ea 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -34,7 +34,15 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Inner, "ST_Inner"), std::make_pair(SurfaceIdEnum::ST_Outer, "ST_Outer"), std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), - std::make_pair(SurfaceIdEnum::TCRV, "TCRV") + std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), + std::make_pair(SurfaceIdEnum::CalD0_Front, "CalD0_Front"), + std::make_pair(SurfaceIdEnum::CalD0_Back, "CalD0_Back"), + std::make_pair(SurfaceIdEnum::CalD1_Front, "CalD1_Front"), + std::make_pair(SurfaceIdEnum::CalD1_Back, "CalD1_Back"), + std::make_pair(SurfaceIdEnum::CalD0_Inner, "CalD0_Inner"), + std::make_pair(SurfaceIdEnum::CalD0_Outer "CalD0_Outer"), + std::make_pair(SurfaceIdEnum::CalD1_Inner, "CalD1_Inner"), + std::make_pair(SurfaceIdEnum::CalD1_Outer, "CalD1_Outer") }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/KinKalGeom/inc/SurfaceMap.hh b/KinKalGeom/inc/SurfaceMap.hh index ea9d455f9e..c533775d01 100644 --- a/KinKalGeom/inc/SurfaceMap.hh +++ b/KinKalGeom/inc/SurfaceMap.hh @@ -31,6 +31,7 @@ namespace mu2e { auto const& DS() const {return ds_; } auto const& ST() const {return st_; } auto const& tracker() const {return tracker_; } + auto const& calo() const {return calo_; } auto const& TCRV() const {return tcrv_; } private: // local copy of detector objects; these hold the actual (typed) surface objects diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index 49aea5119f..aef1c10024 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/TrkReco/fcl/prolog.fcl @@ -8,6 +8,14 @@ BEGIN_PROLOG # # may use a more careful optimization, but this is it for now #------------------------------------------------------------------------------ +InterdetectorTimeBootstrap : { + kalseedModuleLabel : "KKLine" + caloClusterModuleLabel : "CaloClusterMaker" + crvClusterModuleLabel : "SelectReco" + maxDeltaR : 1000. + maxDeltaT : 1000. +} + TrkReco: { McUtils : { tool_type : "TrkRecoMcUtils" comboHitCollTag : "makeSH" diff --git a/TrkReco/src/InterdetectorTimeBootstrap_module.cc b/TrkReco/src/InterdetectorTimeBootstrap_module.cc new file mode 100644 index 0000000000..0651ad9b1c --- /dev/null +++ b/TrkReco/src/InterdetectorTimeBootstrap_module.cc @@ -0,0 +1,129 @@ +// An art analyzer module for the initial calibration of interdetector timing. + +#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 "Offline/RecoDataProducts/inc/KalSeed.hh" +#include "Offline/RecoDataProducts/inc/CaloCluster.hh" +#include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" +#include "Offline/GeometryService/inc/GeometryService.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" +#include "Offline/CosmicRayShieldGeom/inc/CosmicRayShield.hh" +#include "Offline/TrackerGeom/inc/Tracker.hh" +#include "CLHEP/Vector/ThreeVector.h" +#include "TTree.h" +#include "art_root_io/TFileService.h" +#include +#include + +using namespace std; + +namespace mu2e { + + class InterdetectorTimeBootstrap : public art::EDAnalyzer { + + public: + + struct Config { + using Name=fhicl::Name; + using Comment=fhicl::Comment; + fhicl::Atom diagLevel{Name("diagLevel"),Comment("diag level"),0}; + fhicl::Atom kalTag{Name("kalseedModuleLabel"), Comment("kal seed info")}; + fhicl::Atom caloTag{Name("caloClusterModuleLabel"), Comment("calo cluster info")}; + fhicl::Atom crvTag{Name("crvClusterModuleLabel"), Comment("crv info")}; + fhicl::Atom maxDeltaR{Name("maxDeltaR"),Comment("maxDeltaR"),1000.}; + fhicl::Atom maxDeltaT{Name("maxDeltaT"),Comment("maxDeltaT"),1000.}; + }; + typedef art::EDAnalyzer::Table Parameters; + + explicit InterdetectorTimeBootstrap(const Parameters& conf); + virtual ~InterdetectorTimeBootstrap() {} + + + virtual void beginJob(); + virtual void endJob(); + virtual void analyze(const art::Event& e) override; + + private: + Config _conf; + int _diagLevel; + art::InputTag _kalseedToken; + art::InputTag _caloClusterToken; + art::InputTag _crvClusterToken; + double _maxDeltaR; + double _maxDeltaT; + TTree* _timingTree; + double _deltaT_trkcal; + double _deltaT_trkcrv; +}; + + + InterdetectorTimeBootstrap::InterdetectorTimeBootstrap(const Parameters& conf): + art::EDAnalyzer(conf), + _diagLevel(conf().diagLevel()), + _kalseedToken(conf().kalTag()), + _caloClusterToken(conf().caloTag()), + _crvClusterToken(conf().crvTag()), + _maxDeltaR(conf().maxDeltaR()), + _maxDeltaT(conf().maxDeltaT()) + {} + + //InterdetectorTimeBootstrap::~InterdetectorTimeBootstrap() {} + + void InterdetectorTimeBootstrap::beginJob() { + art::ServiceHandle tfs; + _timingTree = tfs->make("timingTree", "Timing Differences"); + _timingTree->Branch("deltaT_trkcal", &_deltaT_trkcal, "deltaT_trkcal/F"); + _timingTree->Branch("deltaT_trkcrv", &_deltaT_trkcrv, "deltaT_trkcrv/F"); + } + + void InterdetectorTimeBootstrap::analyze(const art::Event& event) { + // Access event data products for all three main detectors: + art::Handle tracksHandle; + event.getByLabel(_kalseedToken, tracksHandle); + const KalSeedCollection& tracks = *tracksHandle; + + art::Handle clustersHandle; + event.getByLabel(_caloClusterToken, clustersHandle); + const CaloClusterCollection& clusters = *clustersHandle; + + art::Handle crvHandle; + event.getByLabel(_crvClusterToken, crvHandle); + const CrvCoincidenceClusterCollection& CRVclusters = *crvHandle; + + // Get geometry handles + //GeomHandle calo; + //GeomHandle tracker; + //GeomHandle crv; + + for (const auto& track : tracks) { + // Use the track' parameters for extrapolation + CLHEP::Hep3Vector track_pos_at_calo; // Position of track extrapolated to the calo face + + // Use specific Mu2e geometry functions to find the intersection + bool intersects = true; // Use geometry to check intersection + + if (intersects) { + for (const auto& cluster : clusters) { + const CLHEP::Hep3Vector& cluster_pos = cluster.cog3Vector(); + // Simple geometric match: check proximity + if ((track_pos_at_calo - cluster_pos).mag() < _maxDeltaR) { + // Calculate time difference + // Time of flight correction would be applied here in a real scenario + double time_calo = cluster.time(); + // The track time needs to be adjusted for the time of flight from the track origin to the calo + double time_track_extrapolated = track.time() + (track_pos_at_calo.mag() / 299.792458); // TOF in ns + _deltaT_trkcal = time_calo - time_track_extrapolated; + + // Fill the NTuple for later analysis (e.g., histogramming the mean) + _timingTree->Fill(); + } + } + } + } + } +} // namespace mu2e + +DEFINE_ART_MODULE(mu2e::InterdetectorTimeBootstrap) From 9162348925f3163eaadae823c40143a889d473cb Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 10:53:36 -0600 Subject: [PATCH 02/31] added new Kinkal vols --- KinKalGeom/inc/Calo.hh | 50 ++++++++++++++++++++++++++++++++++++ KinKalGeom/inc/SurfaceMap.hh | 1 + KinKalGeom/src/Calo.cc | 22 ++++++++++++++++ 3 files changed, 73 insertions(+) create mode 100644 KinKalGeom/inc/Calo.hh create mode 100644 KinKalGeom/src/Calo.cc diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh new file mode 100644 index 0000000000..ad1b8e31db --- /dev/null +++ b/KinKalGeom/inc/Calo.hh @@ -0,0 +1,50 @@ +// +// Define the nominal calorimeter boundary and reference surfaces, used to extrapolate and sample KinKal track fits, and to build +// the passive materials in the fit +// original author: Sophie Middleton (2025) +// +#ifndef KinKalGeom_Calorimeter_hh +#define KinKalGeom_Calorimeter_hh +#include "KinKal/Geometry/Cylinder.hh" +#include "KinKal/Geometry/Disk.hh" +#include "KinKal/Geometry/Annulus.hh" +#include +namespace mu2e { + namespace KinKalGeom { + class Calorimeter { + public: + using CylPtr = std::shared_ptr; + using DiskPtr = std::shared_ptr; + // default constructor with nominal geometry + Calorimeter(); + // accessors - d0 - disk 0, d1 - disk 1 + // return by reference + auto const& d0_outer() const { return *d0_outer_; } + auto const& d0_inner() const { return *d0_inner_; } + auto const& d0_front() const { return *d0_front_; } + auto const& d0_back() const { return *d0_back_; } + + auto const& d1_outer() const { return *d1_outer_; } + auto const& d1_inner() const { return *d1_inner_; } + auto const& d1_front() const { return *d1_front_; } + auto const& d1_back() const { return *d1_back_; } + + // return by ptr + auto const& d0outerPtr() const { return d0_outer_; } + auto const& d0innerPtr() const { return d0_inner_; } + auto const& d0frontPtr() const { return d0_front_; } + auto const& d0backPtr() const { return d0_back_; } + + auto const& d1outerPtr() const { return d1_outer_; } + auto const& d1innerPtr() const { return d1_inner_; } + auto const& d1frontPtr() const { return d1_front_; } + auto const& d1backPtr() const { return d1_back_; } + + private: + CylPtr d0_outer_, d0_inner_, d1_outer_, d1_inner_; // active volume boundary + DiskPtr d0_front_, d0_back_, d1_front_, d1_back_; // disk front and back + }; + } +} + +#endif diff --git a/KinKalGeom/inc/SurfaceMap.hh b/KinKalGeom/inc/SurfaceMap.hh index c533775d01..2cf935a14d 100644 --- a/KinKalGeom/inc/SurfaceMap.hh +++ b/KinKalGeom/inc/SurfaceMap.hh @@ -6,6 +6,7 @@ #define KinKalGeom_SurfaceMap_hh #include "Offline/KinKalGeom/inc/Tracker.hh" #include "Offline/KinKalGeom/inc/StoppingTarget.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" #include "Offline/KinKalGeom/inc/DetectorSolenoid.hh" #include "Offline/KinKalGeom/inc/TestCRV.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc new file mode 100644 index 0000000000..5e60587912 --- /dev/null +++ b/KinKalGeom/src/Calo.cc @@ -0,0 +1,22 @@ +#include "Offline/KinKalGeom/inc/Calo.hh" +namespace mu2e { + namespace KinKalGeom { + using KinKal::VEC3; + using KinKal::Cylinder; + using KinKal::Disk; + // currently use hard-coded geometry. Note: these are only comparable to the MC virtual detector positions + // at the ~10 um level, as the G4 virtual detector steppoints are random by roughly that amount (step tolerance) + //FIXME need to update values for the Calorimeter + Calo::Calo() : + d0_outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, + d0_inner_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),376.9,1635.11)}, + d0_front_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + d0_back_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1639.10),950.)} + + d1_outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, + d1_inner_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),376.9,1635.11)}, + d1_front_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + d1_back_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1639.10),950.)} + {} + } +} From 3b9699e1aa9c420551cde4a610d6e99b047978b7 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 13:04:45 -0600 Subject: [PATCH 03/31] added many surfaces, need positions --- DataProducts/inc/SurfaceId.hh | 9 +++-- DataProducts/src/SurfaceId.cc | 27 +++++++++----- KinKalGeom/inc/Calo.hh | 67 +++++++++++++++++++++++------------ KinKalGeom/inc/SurfaceMap.hh | 1 + KinKalGeom/src/Calo.cc | 13 ++----- KinKalGeom/src/SurfaceMap.cc | 6 +++- 6 files changed, 78 insertions(+), 45 deletions(-) diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index b9fbde90f0..2467589aee 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -21,7 +21,12 @@ namespace mu2e { OPA=95, TSDA, // Absorbers in the DS ST_Front=100,ST_Back, ST_Inner, ST_Outer, ST_Foils, ST_Wires, // stopping target bounding surfaces and components TCRV=200 // CRV test planes - CalD0_Front=300, CalD0_Back, CalD1_Front, CalD1_Back, CalD0_Inner, CalD0_Outer, CalD1_Inner, CalD1_Outer //calo VD + // front/back = z positions + // edge = top/bottom positions + // surf = left/right positions + EMC_Disk_0_SurfIn = 300, EMC_Disk_0_SurfOut, EMC_Disk_1_SurfIn, EMC_Disk_1_SurfOut, + EMC_Disk_0_EdgeIn, EMC_Disk_0_EdgeOut,EMC_Disk_1_EdgeIn, EMC_Disk_1_EdgeOut, EMC_0_FrontIn, EMC_0_FrontOut, EMC_1_FrontIn, EMC_1_FrontOut, EMC_2_FrontIn, EMC_2_FrontOut, EMC_3_FrontIn, EMC_3_FrontOut + }; static std::string const& typeName(); @@ -39,7 +44,7 @@ namespace mu2e { // forward some accessors auto const& id() const { return sid_; } int index() const { return index_; } - auto const& name() const { return sid_.name(); } + auto const& name() const { return sid_.name();} bool indexMatch(SurfaceId const& other) const { return index_ == other.index_ || index_ < 0 || other.index_ < 0; } bool indexCompare(SurfaceId const& other) const { return index_<0 || other.index_ < 0 ? false : index_ < other.index_; } diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index a4dc3363ea..75ca3bae38 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -35,14 +35,25 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Outer, "ST_Outer"), std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), - std::make_pair(SurfaceIdEnum::CalD0_Front, "CalD0_Front"), - std::make_pair(SurfaceIdEnum::CalD0_Back, "CalD0_Back"), - std::make_pair(SurfaceIdEnum::CalD1_Front, "CalD1_Front"), - std::make_pair(SurfaceIdEnum::CalD1_Back, "CalD1_Back"), - std::make_pair(SurfaceIdEnum::CalD0_Inner, "CalD0_Inner"), - std::make_pair(SurfaceIdEnum::CalD0_Outer "CalD0_Outer"), - std::make_pair(SurfaceIdEnum::CalD1_Inner, "CalD1_Inner"), - std::make_pair(SurfaceIdEnum::CalD1_Outer, "CalD1_Outer") + + std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfIn, "EMC_Disk_0_SurfIn"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfOut, "EMC_Disk_0_SurfOut"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfIn, "EMC_Disk_1_SurfIn"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfOut, "EMC_Disk_1_SurfOut"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeIn, "EMC_Disk_0_EdgeIn"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeOut "EMC_Disk_0_EdgeOut"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeIn, "EMC_Disk_1_EdgeIn"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeOut, "EMC_Disk_1_EdgeOut"), + + std::make_pair(SurfaceIdEnum::EMC_0_FrontIn, "EMC_0_FrontIn"), + std::make_pair(SurfaceIdEnum::EMC_0_FrontOut, "EMC_0_FrontOut"), + std::make_pair(SurfaceIdEnum::EMC_1_FrontIn, "EMC_1_FrontIn"), + std::make_pair(SurfaceIdEnum::EMC_1_FrontOut, "EMC_1_FrontOut"), + std::make_pair(SurfaceIdEnum::EMC_0_FrontIn, "EMC_2_FrontIn"), + std::make_pair(SurfaceIdEnum::EMC_0_FrontOut, "EMC_2_FrontOut"), + std::make_pair(SurfaceIdEnum::EMC_1_FrontIn, "EMC_3_FrontIn"), + std::make_pair(SurfaceIdEnum::EMC_1_FrontOut, "EMC_3_FrontOut"), + }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh index ad1b8e31db..8dd2643fc6 100644 --- a/KinKalGeom/inc/Calo.hh +++ b/KinKalGeom/inc/Calo.hh @@ -17,32 +17,53 @@ namespace mu2e { using DiskPtr = std::shared_ptr; // default constructor with nominal geometry Calorimeter(); - // accessors - d0 - disk 0, d1 - disk 1 // return by reference - auto const& d0_outer() const { return *d0_outer_; } - auto const& d0_inner() const { return *d0_inner_; } - auto const& d0_front() const { return *d0_front_; } - auto const& d0_back() const { return *d0_back_; } - - auto const& d1_outer() const { return *d1_outer_; } - auto const& d1_inner() const { return *d1_inner_; } - auto const& d1_front() const { return *d1_front_; } - auto const& d1_back() const { return *d1_back_; } - - // return by ptr - auto const& d0outerPtr() const { return d0_outer_; } - auto const& d0innerPtr() const { return d0_inner_; } - auto const& d0frontPtr() const { return d0_front_; } - auto const& d0backPtr() const { return d0_back_; } - - auto const& d1outerPtr() const { return d1_outer_; } - auto const& d1innerPtr() const { return d1_inner_; } - auto const& d1frontPtr() const { return d1_front_; } - auto const& d1backPtr() const { return d1_back_; } + auto const& EMC_Disk_0_SurfIn() const { return *EMC_Disk_0_SurfIn_; } + auto const& EMC_Disk_0_SurfOut() const { return *EMC_Disk_0_SurfOut_; } + auto const& EMC_Disk_1_SurfIn() const { return *EMC_Disk_1_SurfIn_; } + auto const& EMC_Disk_1_SurfOut() const { return *EMC_Disk_1_SurfOut_; } + + auto const& EMC_Disk_0_EdgeIn() const { return *EMC_Disk_0_EdgeIn_; } + auto const& EMC_Disk_0_EdgeOut() const { return *EMC_Disk_0_EdgeOut_; } + auto const& EMC_Disk_1_EdgeIn() const { return *EMC_Disk_1_EdgeIn_; } + auto const& EMC_Disk_1_EdgeOut() const { return *EMC_Disk_1_EdgeOut_; } + + auto const& EMC_0_FrontIn() const { return *EMC_0_FrontIn_; } + auto const& EMC_0_FrontOut() const { return *EMC_0_FrontOut_; } + auto const& EMC_1_FrontIn() const { return *EMC_1_FrontIn_; } + auto const& EMC_1_FrontOut() const { return *EMC_1_FrontOut_; } + + auto const& EMC_2_FrontIn() const { return *EMC_2_FrontIn_; } + auto const& EMC_2_FrontOut() const { return *EMC_2_FrontOut_; } + auto const& EMC_3_FrontIn() const { return *EMC_3_FrontIn_; } + auto const& EMC_3_FrontOut() const { return *EMC_3_FrontOut_; } + +// return by ptr + auto const& EMC_Disk_0_SurfInPtr() const { return EMC_Disk_0_SurfIn_;} + auto const& EMC_Disk_0_SurfOutPtr() const { return EMC_Disk_0_SurfOut_;} + auto const& EMC_Disk_1_SurfInPtr() const { return EMC_Disk_1_SurfIn_;} + auto const& EMC_Disk_1_SurfOutPtr() const { return EMC_Disk_1_SurfOut_;} + + auto const& EMC_Disk_0_EdgeInPtr() const { return EMC_Disk_0_EdgeIn_;} + auto const& EMC_Disk_0_EdgeOutPtr() const { return EMC_Disk_0_EdgeOut_;} + auto const& EMC_Disk_1_EdgeInPtr() const { return EMC_Disk_1_EdgeIn_;} + auto const& EMC_Disk_1_EdgeOutPtr() const { return EMC_Disk_1_EdgeOut_;} + + auto const& EMC_0_FrontInPtr() const { return EMC_0_FrontIn_;} + auto const& EMC_0_FrontOutPtr() const { return EMC_0_FrontOut_;} + auto const& EMC_1_FrontInPtr() const { return EMC_1_FrontIn_;} + auto const& EMC_1_FrontOutPtr() const { return EMC_1_FrontOut_;} + + auto const& EMC_2_FrontInPtr() const { return EMC_2_FrontIn_;} + auto const& EMC_2_FrontOutPtr() const { return EMC_2_FrontOut_;} + auto const& EMC_3_FrontInPtr() const { return EMC_3_FrontIn_;} + auto const& EMC_3_FrontOutPtr() const { return EMC_3_FrontOut_;} private: - CylPtr d0_outer_, d0_inner_, d1_outer_, d1_inner_; // active volume boundary - DiskPtr d0_front_, d0_back_, d1_front_, d1_back_; // disk front and back + CylPtr EMC_Disk_0_SurfIn_, EMC_Disk_0_SurfOut_, EMC_Disk_1_SurfIn_, EMC_Disk_1_SurfOut_, EMC_Disk_0_EdgeIn_, EMC_Disk_0_EdgeOut_,EMC_Disk_1_EdgeIn_, EMC_Disk_1_EdgeOut_; // active volume boundary in XY amd YZ + + DiskPtr EMC_0_FrontIn_, EMC_0_FrontOut_, EMC_1_FrontIn_, EMC_1_FrontOut_, EMC_2_FrontIn_, EMC_2_FrontOut_, EMC_3_FrontIn_, EMC_3_FrontOut_; + }; } } diff --git a/KinKalGeom/inc/SurfaceMap.hh b/KinKalGeom/inc/SurfaceMap.hh index 2cf935a14d..e9c4826b66 100644 --- a/KinKalGeom/inc/SurfaceMap.hh +++ b/KinKalGeom/inc/SurfaceMap.hh @@ -40,6 +40,7 @@ namespace mu2e { KinKalGeom::StoppingTarget st_; KinKalGeom::DetectorSolenoid ds_; KinKalGeom::TestCRV tcrv_; + KinKalGeom::Calo calo_; // the actual map std::multimap map_; }; diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index 5e60587912..782f9315af 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -4,19 +4,10 @@ namespace mu2e { using KinKal::VEC3; using KinKal::Cylinder; using KinKal::Disk; - // currently use hard-coded geometry. Note: these are only comparable to the MC virtual detector positions - // at the ~10 um level, as the G4 virtual detector steppoints are random by roughly that amount (step tolerance) - //FIXME need to update values for the Calorimeter + Calo::Calo() : - d0_outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, - d0_inner_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),376.9,1635.11)}, - d0_front_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - d0_back_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1639.10),950.)} + EMC_Disk_0_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, - d1_outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, - d1_inner_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),376.9,1635.11)}, - d1_front_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - d1_back_{ std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1639.10),950.)} {} } } diff --git a/KinKalGeom/src/SurfaceMap.cc b/KinKalGeom/src/SurfaceMap.cc index de9e45d930..235cc27954 100644 --- a/KinKalGeom/src/SurfaceMap.cc +++ b/KinKalGeom/src/SurfaceMap.cc @@ -27,10 +27,14 @@ namespace mu2e { for(size_t ifoil=0;ifoil < st_.foils().size();++ifoil){ map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::ST_Foils,ifoil),std::static_pointer_cast(st_.foilPtr(ifoil)))); } - // test CRV; Planes are numbered by their vertical (y) position +// test CRV; Planes are numbered by their vertical (y) position map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TCRV,0),std::static_pointer_cast(tcrv_.t1Ptr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TCRV,1),std::static_pointer_cast(tcrv_.ex1Ptr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TCRV,2),std::static_pointer_cast(tcrv_.t2Ptr()))); + + //calo FIXME - need to remove the frontPtr + map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_SurfIn),std::static_pointer_cast(calo_.EMC_Disk_0_SurfInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_SurfOut),std::static_pointer_cast(calo_.EMC_Disk_0_SurfOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_SurfIn),std::static_pointer_cast(calo_.EMC_Disk_1_SurfInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_SurfOut),std::static_pointer_cast(calo_.EMC_Disk_1_SurfOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_EdgeIn),std::static_pointer_cast(calo_.EMC_Disk_0_EdgeInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_EdgeOut),std::static_pointer_cast(calo_.EMC_Disk_0_EdgeOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_EdgeIn),std::static_pointer_cast(calo_.EMC_Disk_1_EdgeInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_EdgeOut),std::static_pointer_cast(calo_.EMC_Disk_1_EdgeOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_0_FrontIn),std::static_pointer_cast(calo_.EMC_0_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_0_FrontOut),std::static_pointer_cast(calo_.EMC_0_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_1_FrontIn),std::static_pointer_cast(calo_.EMC_1_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_1_FrontOut),std::static_pointer_cast(calo_.EMC_1_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_2_FrontIn),std::static_pointer_cast(calo_.EMC_2_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_2_FrontOut),std::static_pointer_cast(calo_.EMC_2_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_3_FrontIn),std::static_pointer_cast(calo_.EMC_3_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_3_FrontOut),std::static_pointer_cast(calo_.EMC_3_FrontOutPtr()))); + } void SurfaceMap::surfaces(SurfaceIdCollection const& ids,SurfacePairCollection& surfs) const { surfs.clear(); From 8219beb4ef7940d753cb088399d7dd5f22fd3b56 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 14:59:55 -0600 Subject: [PATCH 04/31] updated coords --- KinKalGeom/src/Calo.cc | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index 782f9315af..8ff2e1ae7e 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -5,8 +5,41 @@ namespace mu2e { using KinKal::Cylinder; using KinKal::Disk; + // D0 = 11842 -10175.0 = 1667 + // D1 = 13220 -10175.0 = 3045 + /* + double calorimeter.caloDiskRadiusIn = 335; + double calorimeter.caloDiskRadiusOut = 719; + half length = 192.295 + + D0 front = 1667 + D0 back = + D1 front = + D1 back = + + Cylinder(VEC3 const& axis, VEC3 const& center, double radius, double halflen ); + Disk(VEC3 const& norm,VEC3 const& udir, VEC3 const& center, double radius) + */ Calo::Calo() : - EMC_Disk_0_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,4.0),850.11,1635.11)}, + EMC_Disk_0_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, + EMC_Disk_0_SurfOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, + EMC_Disk_1_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, + EMC_Disk_1_SurfOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, + + EMC_Disk_0_EdgeIn_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,1667),719,192.295)}, + EMC_Disk_0_EdgeOut_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,1667),719,192.295)}, + EMC_Disk_1_EdgeIn_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,3045),719,192.295)}, + EMC_Disk_1_EdgeOut_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,3045),719,192.295)}, + + EMC_0_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_0_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_1_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_1_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + + EMC_2_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_2_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_3_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_3_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, {} } From 20c7dc82c65d17b347f38d245be80a0d106e5d0b Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 15:40:29 -0600 Subject: [PATCH 05/31] updated surfaces --- DataProducts/inc/SurfaceId.hh | 4 +-- DataProducts/src/SurfaceId.cc | 26 +++++--------- KinKalGeom/inc/Calo.hh | 65 +++++++++++------------------------ KinKalGeom/src/Calo.cc | 30 +++++----------- KinKalGeom/src/SurfaceMap.cc | 12 +++++-- 5 files changed, 48 insertions(+), 89 deletions(-) diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index 2467589aee..87d7f6b955 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -24,8 +24,8 @@ namespace mu2e { // front/back = z positions // edge = top/bottom positions // surf = left/right positions - EMC_Disk_0_SurfIn = 300, EMC_Disk_0_SurfOut, EMC_Disk_1_SurfIn, EMC_Disk_1_SurfOut, - EMC_Disk_0_EdgeIn, EMC_Disk_0_EdgeOut,EMC_Disk_1_EdgeIn, EMC_Disk_1_EdgeOut, EMC_0_FrontIn, EMC_0_FrontOut, EMC_1_FrontIn, EMC_1_FrontOut, EMC_2_FrontIn, EMC_2_FrontOut, EMC_3_FrontIn, EMC_3_FrontOut + EMC_Disk_0_Outer = 300, EMC_Disk_0_Inner, EMC_Disk_1_Inner, EMC_Disk_1_Outer, + EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back }; diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 75ca3bae38..96442b88b5 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -36,24 +36,14 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfIn, "EMC_Disk_0_SurfIn"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfOut, "EMC_Disk_0_SurfOut"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfIn, "EMC_Disk_1_SurfIn"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfOut, "EMC_Disk_1_SurfOut"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeIn, "EMC_Disk_0_EdgeIn"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeOut "EMC_Disk_0_EdgeOut"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeIn, "EMC_Disk_1_EdgeIn"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeOut, "EMC_Disk_1_EdgeOut"), - - std::make_pair(SurfaceIdEnum::EMC_0_FrontIn, "EMC_0_FrontIn"), - std::make_pair(SurfaceIdEnum::EMC_0_FrontOut, "EMC_0_FrontOut"), - std::make_pair(SurfaceIdEnum::EMC_1_FrontIn, "EMC_1_FrontIn"), - std::make_pair(SurfaceIdEnum::EMC_1_FrontOut, "EMC_1_FrontOut"), - std::make_pair(SurfaceIdEnum::EMC_0_FrontIn, "EMC_2_FrontIn"), - std::make_pair(SurfaceIdEnum::EMC_0_FrontOut, "EMC_2_FrontOut"), - std::make_pair(SurfaceIdEnum::EMC_1_FrontIn, "EMC_3_FrontIn"), - std::make_pair(SurfaceIdEnum::EMC_1_FrontOut, "EMC_3_FrontOut"), - + std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfIn, "EMC_Disk_0_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfOut, "EMC_Disk_0_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfIn, "EMC_Disk_1_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfOut, "EMC_Disk_1_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeIn, "EMC_Disk_0_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeOut "EMC_Disk_1_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeIn, "EMC_Disk_0_Back"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeOut, "EMC_Disk_1_Back") }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh index 8dd2643fc6..3504dbd1c0 100644 --- a/KinKalGeom/inc/Calo.hh +++ b/KinKalGeom/inc/Calo.hh @@ -18,52 +18,29 @@ namespace mu2e { // default constructor with nominal geometry Calorimeter(); // return by reference - auto const& EMC_Disk_0_SurfIn() const { return *EMC_Disk_0_SurfIn_; } - auto const& EMC_Disk_0_SurfOut() const { return *EMC_Disk_0_SurfOut_; } - auto const& EMC_Disk_1_SurfIn() const { return *EMC_Disk_1_SurfIn_; } - auto const& EMC_Disk_1_SurfOut() const { return *EMC_Disk_1_SurfOut_; } - - auto const& EMC_Disk_0_EdgeIn() const { return *EMC_Disk_0_EdgeIn_; } - auto const& EMC_Disk_0_EdgeOut() const { return *EMC_Disk_0_EdgeOut_; } - auto const& EMC_Disk_1_EdgeIn() const { return *EMC_Disk_1_EdgeIn_; } - auto const& EMC_Disk_1_EdgeOut() const { return *EMC_Disk_1_EdgeOut_; } - - auto const& EMC_0_FrontIn() const { return *EMC_0_FrontIn_; } - auto const& EMC_0_FrontOut() const { return *EMC_0_FrontOut_; } - auto const& EMC_1_FrontIn() const { return *EMC_1_FrontIn_; } - auto const& EMC_1_FrontOut() const { return *EMC_1_FrontOut_; } - - auto const& EMC_2_FrontIn() const { return *EMC_2_FrontIn_; } - auto const& EMC_2_FrontOut() const { return *EMC_2_FrontOut_; } - auto const& EMC_3_FrontIn() const { return *EMC_3_FrontIn_; } - auto const& EMC_3_FrontOut() const { return *EMC_3_FrontOut_; } - -// return by ptr - auto const& EMC_Disk_0_SurfInPtr() const { return EMC_Disk_0_SurfIn_;} - auto const& EMC_Disk_0_SurfOutPtr() const { return EMC_Disk_0_SurfOut_;} - auto const& EMC_Disk_1_SurfInPtr() const { return EMC_Disk_1_SurfIn_;} - auto const& EMC_Disk_1_SurfOutPtr() const { return EMC_Disk_1_SurfOut_;} - - auto const& EMC_Disk_0_EdgeInPtr() const { return EMC_Disk_0_EdgeIn_;} - auto const& EMC_Disk_0_EdgeOutPtr() const { return EMC_Disk_0_EdgeOut_;} - auto const& EMC_Disk_1_EdgeInPtr() const { return EMC_Disk_1_EdgeIn_;} - auto const& EMC_Disk_1_EdgeOutPtr() const { return EMC_Disk_1_EdgeOut_;} - - auto const& EMC_0_FrontInPtr() const { return EMC_0_FrontIn_;} - auto const& EMC_0_FrontOutPtr() const { return EMC_0_FrontOut_;} - auto const& EMC_1_FrontInPtr() const { return EMC_1_FrontIn_;} - auto const& EMC_1_FrontOutPtr() const { return EMC_1_FrontOut_;} - - auto const& EMC_2_FrontInPtr() const { return EMC_2_FrontIn_;} - auto const& EMC_2_FrontOutPtr() const { return EMC_2_FrontOut_;} - auto const& EMC_3_FrontInPtr() const { return EMC_3_FrontIn_;} - auto const& EMC_3_FrontOutPtr() const { return EMC_3_FrontOut_;} + auto const& EMC_Disk_0_Outer() const { return *EMC_Disk_0_Outer_;} + auto const& EMC_Disk_0_Inner() const { return *EMC_Disk_0_Inner_;} + auto const& EMC_Disk_1_Inner() const { return *EMC_Disk_1_Inner_;} + auto const& EMC_Disk_1_Outer() const { return *EMC_Disk_1_Outer_;} + + auto const& EMC_Disk_0_Front() const { return *EMC_Disk_0_Front_;} + auto const& EMC_Disk_1_Front() const { return *EMC_Disk_1_Front_;} + auto const& EMC_Disk_0_Back() const { return *EMC_Disk_0_Back_;} + auto const& EMC_Disk_1_Back() const { return *EMC_Disk_1_Back_;} + + // return by ptr + auto const& EMC_Disk_0_OuterPtr() const { return EMC_Disk_0_Outer_;} + auto const& EMC_Disk_0_InnerPtr() const { return EMC_Disk_0_Inner_;} + auto const& EMC_Disk_1_InnerPtr() const { return EMC_Disk_1_Inner_;} + auto const& EMC_Disk_1_OuterPtr() const { return EMC_Disk_1_Outer_;} + auto const& EMC_Disk_0_FrontPtr() const { return EMC_Disk_0_Front_;} + auto const& EMC_Disk_1_FrontPtr() const { return EMC_Disk_1_Front_;} + auto const& EMC_Disk_0_BackPtr() const { return EMC_Disk_0_Back_;} + auto const& EMC_Disk_1_BackPtr() const { return EMC_Disk_1_Back_;} private: - CylPtr EMC_Disk_0_SurfIn_, EMC_Disk_0_SurfOut_, EMC_Disk_1_SurfIn_, EMC_Disk_1_SurfOut_, EMC_Disk_0_EdgeIn_, EMC_Disk_0_EdgeOut_,EMC_Disk_1_EdgeIn_, EMC_Disk_1_EdgeOut_; // active volume boundary in XY amd YZ - - DiskPtr EMC_0_FrontIn_, EMC_0_FrontOut_, EMC_1_FrontIn_, EMC_1_FrontOut_, EMC_2_FrontIn_, EMC_2_FrontOut_, EMC_3_FrontIn_, EMC_3_FrontOut_; - + CylPtr EMC_Disk_0_Outer, EMC_Disk_0_Inner, EMC_Disk_1_Inner, EMC_Disk_1_Outer; + DiskPtr EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back; }; } } diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index 8ff2e1ae7e..d649450b54 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -12,34 +12,20 @@ namespace mu2e { double calorimeter.caloDiskRadiusOut = 719; half length = 192.295 - D0 front = 1667 - D0 back = - D1 front = - D1 back = Cylinder(VEC3 const& axis, VEC3 const& center, double radius, double halflen ); Disk(VEC3 const& norm,VEC3 const& udir, VEC3 const& center, double radius) */ Calo::Calo() : - EMC_Disk_0_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, - EMC_Disk_0_SurfOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, - EMC_Disk_1_SurfIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, - EMC_Disk_1_SurfOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, + EMC_Disk_0_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),335,192.295)}, + EMC_Disk_0_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, + EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),335,192.295)}, + EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, - EMC_Disk_0_EdgeIn_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,1667),719,192.295)}, - EMC_Disk_0_EdgeOut_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,1667),719,192.295)}, - EMC_Disk_1_EdgeIn_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,3045),719,192.295)}, - EMC_Disk_1_EdgeOut_ { std::make_shared(VEC3(0.0,1.0,0.0),VEC3(0.0,0.0,3045),719,192.295)}, - - EMC_0_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_0_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_1_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_1_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - - EMC_2_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_2_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_3_FrontIn_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, - EMC_3_FrontOut_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-1631.12),950.)}, + EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667),719.)}, + EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+2*192.295),719.)}, + EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045),719.)}, + EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045 + +2*192.295 ),719.)}, {} } diff --git a/KinKalGeom/src/SurfaceMap.cc b/KinKalGeom/src/SurfaceMap.cc index 235cc27954..3a093d0de4 100644 --- a/KinKalGeom/src/SurfaceMap.cc +++ b/KinKalGeom/src/SurfaceMap.cc @@ -32,9 +32,15 @@ namespace mu2e { map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TCRV,1),std::static_pointer_cast(tcrv_.ex1Ptr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TCRV,2),std::static_pointer_cast(tcrv_.t2Ptr()))); - //calo FIXME - need to remove the frontPtr - map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_SurfIn),std::static_pointer_cast(calo_.EMC_Disk_0_SurfInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_SurfOut),std::static_pointer_cast(calo_.EMC_Disk_0_SurfOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_SurfIn),std::static_pointer_cast(calo_.EMC_Disk_1_SurfInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_SurfOut),std::static_pointer_cast(calo_.EMC_Disk_1_SurfOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_EdgeIn),std::static_pointer_cast(calo_.EMC_Disk_0_EdgeInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_EdgeOut),std::static_pointer_cast(calo_.EMC_Disk_0_EdgeOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_EdgeIn),std::static_pointer_cast(calo_.EMC_Disk_1_EdgeInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_EdgeOut),std::static_pointer_cast(calo_.EMC_Disk_1_EdgeOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_0_FrontIn),std::static_pointer_cast(calo_.EMC_0_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_0_FrontOut),std::static_pointer_cast(calo_.EMC_0_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_1_FrontIn),std::static_pointer_cast(calo_.EMC_1_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_1_FrontOut),std::static_pointer_cast(calo_.EMC_1_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_2_FrontIn),std::static_pointer_cast(calo_.EMC_2_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_2_FrontOut),std::static_pointer_cast(calo_.EMC_2_FrontOutPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_3_FrontIn),std::static_pointer_cast(calo_.EMC_3_FrontInPtr()))); map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_3_FrontOut),std::static_pointer_cast(calo_.EMC_3_FrontOutPtr()))); - + //calo + map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Outer),std::static_pointer_cast(calo_.EMC_Disk_0_OuterPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Inner),std::static_pointer_cast(calo_.EMC_Disk_0_InnerPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Inner),std::static_pointer_cast(calo_.EMC_Disk_1_InnerPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Outer),std::static_pointer_cast(calo_.EMC_Disk_1_OuterPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Front),std::static_pointer_cast(calo_.EMC_Disk_0_FrontPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Front),std::static_pointer_cast(calo_.EMC_Disk_1_FrontPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Back),std::static_pointer_cast(calo_.EMC_Disk_0_BackPtr()))); +map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Back),std::static_pointer_cast(calo_.EMC_Disk_1_BackPtr()))); } void SurfaceMap::surfaces(SurfaceIdCollection const& ids,SurfacePairCollection& surfs) const { surfs.clear(); From 0a5a9857ba27adc2f816b4301107e57f1a17f593 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 25 Nov 2025 16:39:34 -0600 Subject: [PATCH 06/31] builds --- DataProducts/inc/SurfaceId.hh | 5 +---- DataProducts/src/SurfaceId.cc | 16 ++++++++-------- KinKalGeom/CMakeLists.txt | 1 + KinKalGeom/inc/Calo.hh | 12 ++++++------ KinKalGeom/src/Calo.cc | 3 +-- TrkReco/src/InterdetectorTimeBootstrap_module.cc | 4 ++-- 6 files changed, 19 insertions(+), 22 deletions(-) diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index 87d7f6b955..e870cd3664 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -20,10 +20,7 @@ namespace mu2e { IPA=90, IPA_Front, IPA_Back, OPA=95, TSDA, // Absorbers in the DS ST_Front=100,ST_Back, ST_Inner, ST_Outer, ST_Foils, ST_Wires, // stopping target bounding surfaces and components - TCRV=200 // CRV test planes - // front/back = z positions - // edge = top/bottom positions - // surf = left/right positions + TCRV=200, // CRV test planes EMC_Disk_0_Outer = 300, EMC_Disk_0_Inner, EMC_Disk_1_Inner, EMC_Disk_1_Outer, EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 96442b88b5..68c0bc058a 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -36,14 +36,14 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfIn, "EMC_Disk_0_Outer"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_SurfOut, "EMC_Disk_0_Inner"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfIn, "EMC_Disk_1_Inner"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_SurfOut, "EMC_Disk_1_Outer"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeIn, "EMC_Disk_0_Front"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_EdgeOut "EMC_Disk_1_Front"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeIn, "EMC_Disk_0_Back"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_EdgeOut, "EMC_Disk_1_Back") + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Outer, "EMC_Disk_0_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Inner, "EMC_Disk_0_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Inner, "EMC_Disk_1_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Outer, "EMC_Disk_1_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Front, "EMC_Disk_0_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Front, "EMC_Disk_1_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Back, "EMC_Disk_0_Back"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Back, "EMC_Disk_1_Back") }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/KinKalGeom/CMakeLists.txt b/KinKalGeom/CMakeLists.txt index eb8b46431c..a9e6b172ab 100644 --- a/KinKalGeom/CMakeLists.txt +++ b/KinKalGeom/CMakeLists.txt @@ -6,6 +6,7 @@ cet_make_library( src/SurfaceMap.cc src/TestCRV.cc src/Tracker.cc + src/Calo.cc LIBRARIES PUBLIC KinKal::Geometry Offline::GeneralUtilities diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh index 3504dbd1c0..a064d5113a 100644 --- a/KinKalGeom/inc/Calo.hh +++ b/KinKalGeom/inc/Calo.hh @@ -3,20 +3,20 @@ // the passive materials in the fit // original author: Sophie Middleton (2025) // -#ifndef KinKalGeom_Calorimeter_hh -#define KinKalGeom_Calorimeter_hh +#ifndef KinKalGeom_Calo_hh +#define KinKalGeom_Calo_hh #include "KinKal/Geometry/Cylinder.hh" #include "KinKal/Geometry/Disk.hh" #include "KinKal/Geometry/Annulus.hh" #include namespace mu2e { namespace KinKalGeom { - class Calorimeter { + class Calo { public: using CylPtr = std::shared_ptr; using DiskPtr = std::shared_ptr; // default constructor with nominal geometry - Calorimeter(); + Calo(); // return by reference auto const& EMC_Disk_0_Outer() const { return *EMC_Disk_0_Outer_;} auto const& EMC_Disk_0_Inner() const { return *EMC_Disk_0_Inner_;} @@ -39,8 +39,8 @@ namespace mu2e { auto const& EMC_Disk_1_BackPtr() const { return EMC_Disk_1_Back_;} private: - CylPtr EMC_Disk_0_Outer, EMC_Disk_0_Inner, EMC_Disk_1_Inner, EMC_Disk_1_Outer; - DiskPtr EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back; + CylPtr EMC_Disk_0_Inner_, EMC_Disk_0_Outer_ , EMC_Disk_1_Inner_, EMC_Disk_1_Outer_; + DiskPtr EMC_Disk_0_Front_, EMC_Disk_0_Back_, EMC_Disk_1_Front_, EMC_Disk_1_Back_; }; } } diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index d649450b54..f3c663e5f8 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -25,8 +25,7 @@ namespace mu2e { EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667),719.)}, EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+2*192.295),719.)}, EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045),719.)}, - EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045 + +2*192.295 ),719.)}, - + EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045 + +2*192.295 ),719.)} {} } } diff --git a/TrkReco/src/InterdetectorTimeBootstrap_module.cc b/TrkReco/src/InterdetectorTimeBootstrap_module.cc index 0651ad9b1c..6f0e40c3ae 100644 --- a/TrkReco/src/InterdetectorTimeBootstrap_module.cc +++ b/TrkReco/src/InterdetectorTimeBootstrap_module.cc @@ -1,4 +1,4 @@ -// An art analyzer module for the initial calibration of interdetector timing. +/* An art analyzer module for the initial calibration of interdetector timing. #include "art/Framework/Core/EDAnalyzer.h" #include "art/Framework/Core/ModuleMacros.h" @@ -126,4 +126,4 @@ namespace mu2e { } } // namespace mu2e -DEFINE_ART_MODULE(mu2e::InterdetectorTimeBootstrap) +DEFINE_ART_MODULE(mu2e::InterdetectorTimeBootstrap)*/ From 10d9916bf1ae98241ddc8a139d09d42c71336648 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Wed, 26 Nov 2025 13:02:01 -0600 Subject: [PATCH 07/31] compiles --- Mu2eKinKal/fcl/prolog.fcl | 2 + Mu2eKinKal/inc/ExtrapolateCalo.hh | 160 ++++++++++++++++++ Mu2eKinKal/inc/KKExtrap.hh | 38 ++++- Mu2eKinKal/inc/KKFitSettings.hh | 1 + .../src/InterdetectorTimeBootstrap_module.cc | 1 - mu2e.fcl | 5 + 6 files changed, 204 insertions(+), 3 deletions(-) create mode 100644 Mu2eKinKal/inc/ExtrapolateCalo.hh create mode 100644 mu2e.fcl diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index ec8343d1fa..2235f9c8af 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -387,6 +387,7 @@ Mu2eKinKal : { Upstream : true BackToTracker : false ToOPA : true + ToCaloD0 : true } LHSEEDXTRAP : { @@ -397,6 +398,7 @@ Mu2eKinKal : { Upstream : false BackToTracker : false ToOPA : false + ToCaloD0 : true } CENTRALHELIX : { diff --git a/Mu2eKinKal/inc/ExtrapolateCalo.hh b/Mu2eKinKal/inc/ExtrapolateCalo.hh new file mode 100644 index 0000000000..68f745c32e --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCalo.hh @@ -0,0 +1,160 @@ +// predicate to extrapolate to calo +// Sophie Middleton (2025) +#ifndef Mu2eKinKal_ExtrapolateCalo_hh +#define Mu2eKinKal_ExtrapolateCalo_hh +#include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" +#include "cetlib_except/exception.h" +#include +#include +#include +namespace mu2e { + using KinKal::TimeDir; + using KinKal::TimeRange; + using KinKalGeom::Calo; + using KinKal::Annulus; + using KinKal::Intersection; + using AnnPtr = std::shared_ptr; + using CaloDisk = std::vector; + using CylPtr = std::shared_ptr; + class ExtrapolateCalo { + public: + ExtrapolateCalo() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10), + d0zmin_(std::numeric_limits::max()), + d0zmax_(-std::numeric_limits::max()), + d1zmin_(std::numeric_limits::max()), + d1zmax_(-std::numeric_limits::max()), + rmin_(std::numeric_limits::max()), + rmax_(-std::numeric_limits::max()), + debug_(0){} + + ExtrapolateCalo(double maxdt, double dptol, double intertol, Calo const& calo, int debug=0) : + maxDt_(maxdt), dptol_(dptol), intertol_(intertol), + d0zmin_( (calo.EMC_Disk_0_Outer().center() - calo.EMC_Disk_0_Outer().axis()*calo.EMC_Disk_0_Outer().halfLength()).Z()), + d0zmax_( (calo.EMC_Disk_0_Outer().center() + calo.EMC_Disk_0_Outer().axis()*calo.EMC_Disk_0_Outer().halfLength()).Z()), + d1zmin_( (calo.EMC_Disk_1_Outer().center() - calo.EMC_Disk_1_Outer().axis()*calo.EMC_Disk_1_Outer().halfLength()).Z()), + d1zmax_( (calo.EMC_Disk_1_Outer().center() + calo.EMC_Disk_1_Outer().axis()*calo.EMC_Disk_1_Outer().halfLength()).Z()), + rmin_( calo.EMC_Disk_0_Inner().radius()), rmax_(calo.EMC_Disk_0_Outer().radius()), + disks_(2),debug_(debug) {} + // interface for extrapolation + double maxDt() const { return maxDt_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + double d0zmin() const { return d0zmin_; } + double d0zmax() const { return d0zmax_; } + double d1zmin() const { return d1zmin_; } + double d1zmax() const { return d1zmax_; } + double rmin() const { return rmin_; } + double rmax() const { return rmax_; } + auto const& disks () const { return disks_; } + auto const& intersection() const { return inter_; } + + int debug() const { return debug_; } + // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + // reset between tracks + void reset() const { inter_ = Intersection(); sid_ = SurfaceId(); ann_ = AnnPtr();} + // find the nearest disk to a z positionin a given z direction + size_t nearestDisk(double zpos, double zdir) const; + private: + double maxDt_; // maximum extrapolation time + double dptol_; // fractional momentum tolerance + double intertol_; // intersection tolerance (mm) + mutable Intersection inter_; // cache of most recent intersection + mutable SurfaceId sid_; // cache of most recent intersection disk SID + mutable AnnPtr ann_; // cache of most recent intersection disk surface + // cache of front and back Z positions + double d0zmin_, d0zmax_; // z range of disk0 volume + double d1zmin_, d1zmax_; // z range of disk1 volume + double rmin_, rmax_; // inner and outer radii of the anuli + CaloDisk disks_; // disks + int debug_; // debug level + }; + + template bool ExtrapolateCalo::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { + // we are answering the question: did the segment last added to this extrapolated track hit a calo disk or not? If so we are done + // extrapolating (for now) and we want to find all the intersections in that piece. If not, and if we're still inside or heading towards the + // disks, keep going. + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + // add a small buffer to the test range to prevent re-intersection with the same piece + static const double epsilon(1e-7); // small difference to avoid re-intersecting + if(ktraj.range().range() <= epsilon) return true; // keep going if the step is very small + auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; + auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + auto vel = ktraj.velocity(stime); // physical velocity + if(tdir == TimeDir::backwards) vel *= -1.0; + auto spos = ktraj.position3(stime); + auto epos = ktraj.position3(etime); + if(debug_ > 2)std::cout << "calo extrap tdir " << tdir << " start z " << spos.Z() << " end z " << epos.Z() << " zvel " << vel.Z() << " rho " << spos.Rho() << std::endl; + // stop if the particle is heading away from the calo + if( (vel.Z() > 0 && spos.Z() > d1zmax_ ) || (vel.Z() < 0 && spos.Z() < d0zmin_)){ + reset(); // clear any cache + if(debug_ > 1)std::cout << "Heading away from calo: done" << std::endl; + return false; + } + // if the particle is going in the right direction but haven't yet reached the calo in Z just keep going + if( (vel.Z() > 0 && epos.Z() < d0zmin_) || (vel.Z() < 0 && epos.Z() > d1zmax_) ){ + reset(); + if(debug_ > 2)std::cout << "Heading towards calo, z " << spos.Z()<< std::endl; + return true; + } + // if we get to here we are in the correct Z range. Test disks. + int idisk = nearestDisk(spos.Z(),vel.Z()); + if(idisk >= (int)disks_.size())return true; + if(debug_ > 2)std::cout << "Looping on disks " << std::endl; + int ddisk = vel.Z() > 0.0 ? 1 : -1; // iteration direction + auto trange = tdir == TimeDir::forwards ? TimeRange(stime,ktraj.range().end()) : TimeRange(ktraj.range().begin(),stime); + // loop over disks in the z range of this piece + while(idisk >= 0 && idisk < (int)disks_.size() && (disks_[idisk]->center().Z() - epos.Z())*ddisk < 0.0){ + auto diskptr = disks_[idisk]; + if(debug_ > 2)std::cout << "disk " << idisk << " z " << diskptr->center().Z() << std::endl; + auto newinter = KinKal::intersect(ktraj,*diskptr,trange,intertol_,tdir); + if(debug_ > 2)std::cout << "calo disk inter " << newinter << std::endl; + if(newinter.good()){ + // update the cache + inter_ = newinter; + ann_ = disks_[idisk]; + //sid_ = SurfaceId(SurfaceIdEnum::calo_Foils,idisk); //FIXME + if(debug_ > 0)std::cout << "Good calo disk " << newinter << " sid " << sid_ << std::endl; + return false; + } + idisk += ddisk; // otherwise continue loopin on disks + } + // no more intersections: keep extending in Z till we clear the calo + reset(); + if(debug_ > 1)std::cout << "Extrapolating to calo edge, z " << spos.Z() << std::endl; + if(vel.Z() > 0.0) + return spos.Z() < d1zmax_; + else + return spos.Z() > d0zmin_; + } + + size_t ExtrapolateCalo::nearestDisk(double zpos, double zvel) const { + size_t retval = disks_.size(); + if(zvel > 0.0){ // going forwards in z + for(auto idisk= disks_.begin(); idisk != disks_.end(); idisk++){ + auto const& diskptr = *idisk; + if(diskptr->center().Z() > zpos){ + retval = std::distance(disks_.begin(),idisk); + break; + } + } + } else { + for(auto idisk= disks_.rbegin(); idisk != disks_.rend(); idisk++){ + auto const& diskptr = *idisk; + if(diskptr->center().Z() < zpos){ + auto jdisk = idisk.base()-1; // points to the equivalent forwards object + retval = std::distance(disks_.begin(),jdisk); + break; + } + } + } + return retval; + } + +} +#endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index c475b04e1f..624c928511 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -12,6 +12,7 @@ #include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateCalo.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/Mu2eKinKal/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" @@ -34,6 +35,7 @@ namespace mu2e { template bool extrapolateIPA(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; + template void toCaloD0(KKTrack& ktrk) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; @@ -45,10 +47,11 @@ namespace mu2e { AnnPtr tsdaptr_; DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; FruPtr opaptr_; - bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; - ExtrapolateToZ TSDA_, trackerFront_, trackerBack_; // extrapolation predicate based on Z values + bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCaloD0_; + ExtrapolateToZ TSDA_, trackerFront_, trackerBack_, calod0Front_, calod0Back_; // extrapolation predicate based on Z values ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA ExtrapolateST extrapST_; // extrapolation to intersections with the ST + //ExtrapolateCalo ; // extrapolation to intersections with the Calo double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO }; @@ -68,9 +71,12 @@ namespace mu2e { toOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), upstream_(extrapconfig.Upstream()), + toCaloD0_(extrapconfig.ToCaloD0()), TSDA_(maxdt_,btol_,smap_.DS().upstreamAbsorber().center().Z(),debug_), trackerFront_(maxdt_,btol_,smap_.tracker().front().center().Z(),debug_), trackerBack_(maxdt_,btol_,smap_.tracker().back().center().Z(),debug_), + calod0Front_(maxdt_,btol_,smap_.calo().EMC_Disk_0_Front().center().Z(),debug_), + calod0Back_(maxdt_,btol_,smap_.calo().EMC_Disk_0_Back().center().Z(),debug_), extrapIPA_(maxdt_,btol_,intertol_,smap_.DS().innerProtonAbsorberPtr(),debug_), extrapST_(maxdt_,btol_,intertol_,smap_.ST(),debug_) { @@ -82,12 +88,14 @@ namespace mu2e { // define the time direction according to the fit direction inside the tracker auto const& ftraj = ktrk.fitTraj(); if(toTrackerEnds_)toTrackerEnds(ktrk); + if(toCaloD0_)toCaloD0(ktrk); if(upstream_){ auto dir0 = ftraj.direction(ftraj.t0()); TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); // extrapolate through the IPA in this direction. bool exitsIPA = extrapolateIPA(ktrk,tdir); + if(exitsIPA){ // if it exits out the back, extrapolate through the target bool exitsST = extrapolateST(ktrk,tdir); if(exitsST) { // if it exits out the back, extrapolate to the TSDA (DS rear absorber) @@ -153,6 +161,32 @@ namespace mu2e { } } + template void KKExtrap::toCaloD0(KKTrack& ktrk) const { + auto const& ftraj = ktrk.fitTraj(); + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; + TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; + auto tofront = ktrk.extrapolate(fronttdir,calod0Back_); + auto toback = ktrk.extrapolate(backtdir,calod0Front_); + // record the standard tracker intersections + static const SurfaceId d0_front("EMC_Disk_0_Front"); + static const SurfaceId d0_back("EMC_Disk_0_Back"); + + if(tofront){ + // check the front piece first; that is usually correct + // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection + auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); + auto frontinter = KinKal::intersect(fhel,*trkfrontptr_,fhel.range(),intertol_,fronttdir); + if(frontinter.good()) ktrk.addIntersection(d0_front,frontinter); + } + if(toback){ + // start from the middle + TimeRange brange = ftraj.range(); + auto backinter = KinKal::intersect(ftraj,*trkbackptr_,brange,intertol_,backtdir); + if(backinter.good())ktrk.addIntersection(d0_back,backinter); + } + } + template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 9a03aee67c..a13022aca8 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -114,6 +114,7 @@ namespace mu2e { fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit (ns)") }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; + fhicl::Atom ToCaloD0 { Name("ToCaloD0"), Comment("Extrapolate tracks to the disk0 ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; }; diff --git a/TrkReco/src/InterdetectorTimeBootstrap_module.cc b/TrkReco/src/InterdetectorTimeBootstrap_module.cc index 6f0e40c3ae..02413fae0d 100644 --- a/TrkReco/src/InterdetectorTimeBootstrap_module.cc +++ b/TrkReco/src/InterdetectorTimeBootstrap_module.cc @@ -111,7 +111,6 @@ namespace mu2e { // Simple geometric match: check proximity if ((track_pos_at_calo - cluster_pos).mag() < _maxDeltaR) { // Calculate time difference - // Time of flight correction would be applied here in a real scenario double time_calo = cluster.time(); // The track time needs to be adjusted for the time of flight from the track origin to the calo double time_track_extrapolated = track.time() + (track_pos_at_calo.mag() / 299.792458); // TOF in ns diff --git a/mu2e.fcl b/mu2e.fcl new file mode 100644 index 0000000000..2030ff9d26 --- /dev/null +++ b/mu2e.fcl @@ -0,0 +1,5 @@ +#include "Production/JobConfig/recoMC/OnSpill.fcl" +services.DbService.purpose: MDC2025_best +services.DbService.version: v1_4 +outputs.LoopHelixOutput.fileName: "mcs.mu2e.test.v1.sequencer.art" +//services.GeometryService.inputFile: "Offline/Mu2eG4/geom/geom_run1_a.txt" From 55253949707d37d3583f0be07eab7a0e93ed2dd0 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Wed, 26 Nov 2025 13:11:34 -0600 Subject: [PATCH 08/31] prints out, not tested --- Mu2eKinKal/inc/KKExtrap.hh | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 624c928511..ae93ecd6be 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -171,19 +171,21 @@ namespace mu2e { // record the standard tracker intersections static const SurfaceId d0_front("EMC_Disk_0_Front"); static const SurfaceId d0_back("EMC_Disk_0_Back"); - + std::cout<<"extrapolating to calo d0"< Date: Mon, 1 Dec 2025 10:00:07 -0600 Subject: [PATCH 09/31] added surface steps from front and back of calo --- CommonMC/src/MakeSurfaceSteps_module.cc | 7 +++++ Mu2eKinKal/fcl/prolog.fcl | 1 + Mu2eKinKal/inc/KKExtrap.hh | 41 ++++++++++++++++++++++--- Mu2eKinKal/inc/KKFitSettings.hh | 1 + 4 files changed, 46 insertions(+), 4 deletions(-) diff --git a/CommonMC/src/MakeSurfaceSteps_module.cc b/CommonMC/src/MakeSurfaceSteps_module.cc index b82d0095b6..ed614802b0 100644 --- a/CommonMC/src/MakeSurfaceSteps_module.cc +++ b/CommonMC/src/MakeSurfaceSteps_module.cc @@ -60,9 +60,15 @@ namespace mu2e { vdmap_[VirtualDetectorId(VirtualDetectorId::TT_Back)] = SurfaceId("TT_Back"); vdmap_[VirtualDetectorId(VirtualDetectorId::TT_OutSurf)] = SurfaceId("TT_Outer"); vdmap_[VirtualDetectorId(VirtualDetectorId::TT_InSurf)] = SurfaceId("TT_Inner"); + + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_0_SurfIn)] = SurfaceId("EMC_Disk_0_Front"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_SurfIn)] = SurfaceId("EMC_Disk_1_Front"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_0_SurfOut)] = SurfaceId("EMC_Disk_0_Back"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_SurfOut)] = SurfaceId("EMC_Disk_1_Back"); } void MakeSurfaceSteps::produce(art::Event& event) { + std::cout<<"MakeSurfaceSteps ---> produce"< det; // create output std::unique_ptr ssc(new SurfaceStepCollection); @@ -72,6 +78,7 @@ namespace mu2e { for(auto const& vdspmc : vdspmccol) { // only some VDs are kept auto isid = vdmap_.find(vdspmc.virtualDetectorId()); + std::cout<<" VID "<emplace_back(isid->second,vdspmc,det); // no aggregation of VD hits } auto nvdsteps = ssc->size(); diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 2235f9c8af..0e396b5f63 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -388,6 +388,7 @@ Mu2eKinKal : { BackToTracker : false ToOPA : true ToCaloD0 : true + ToCaloD1 : true } LHSEEDXTRAP : { diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index ae93ecd6be..ba27939a98 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -36,6 +36,7 @@ namespace mu2e { template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; template void toCaloD0(KKTrack& ktrk) const; + template void toCaloD1(KKTrack& ktrk) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; @@ -47,8 +48,8 @@ namespace mu2e { AnnPtr tsdaptr_; DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; FruPtr opaptr_; - bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCaloD0_; - ExtrapolateToZ TSDA_, trackerFront_, trackerBack_, calod0Front_, calod0Back_; // extrapolation predicate based on Z values + bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCaloD0_, toCaloD1_; + ExtrapolateToZ TSDA_, trackerFront_, trackerBack_, calod0Front_, calod0Back_, calod1Front_, calod1Back_; // extrapolation predicate based on Z values ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA ExtrapolateST extrapST_; // extrapolation to intersections with the ST //ExtrapolateCalo ; // extrapolation to intersections with the Calo @@ -72,11 +73,14 @@ namespace mu2e { toTrackerEnds_(extrapconfig.ToTrackerEnds()), upstream_(extrapconfig.Upstream()), toCaloD0_(extrapconfig.ToCaloD0()), + toCaloD1_(extrapconfig.ToCaloD1()), TSDA_(maxdt_,btol_,smap_.DS().upstreamAbsorber().center().Z(),debug_), trackerFront_(maxdt_,btol_,smap_.tracker().front().center().Z(),debug_), trackerBack_(maxdt_,btol_,smap_.tracker().back().center().Z(),debug_), calod0Front_(maxdt_,btol_,smap_.calo().EMC_Disk_0_Front().center().Z(),debug_), calod0Back_(maxdt_,btol_,smap_.calo().EMC_Disk_0_Back().center().Z(),debug_), + calod1Front_(maxdt_,btol_,smap_.calo().EMC_Disk_1_Front().center().Z(),debug_), + calod1Back_(maxdt_,btol_,smap_.calo().EMC_Disk_1_Back().center().Z(),debug_), extrapIPA_(maxdt_,btol_,intertol_,smap_.DS().innerProtonAbsorberPtr(),debug_), extrapST_(maxdt_,btol_,intertol_,smap_.ST(),debug_) { @@ -89,6 +93,7 @@ namespace mu2e { auto const& ftraj = ktrk.fitTraj(); if(toTrackerEnds_)toTrackerEnds(ktrk); if(toCaloD0_)toCaloD0(ktrk); + if(toCaloD1_)toCaloD1(ktrk); if(upstream_){ auto dir0 = ftraj.direction(ftraj.t0()); TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; @@ -166,8 +171,8 @@ namespace mu2e { auto dir0 = ftraj.direction(ftraj.t0()); TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - auto tofront = ktrk.extrapolate(fronttdir,calod0Back_); - auto toback = ktrk.extrapolate(backtdir,calod0Front_); + auto tofront = ktrk.extrapolate(fronttdir,calod0Front_); + auto toback = ktrk.extrapolate(backtdir,calod0Back_); // record the standard tracker intersections static const SurfaceId d0_front("EMC_Disk_0_Front"); static const SurfaceId d0_back("EMC_Disk_0_Back"); @@ -189,6 +194,34 @@ namespace mu2e { } } + template void KKExtrap::toCaloD1(KKTrack& ktrk) const { + auto const& ftraj = ktrk.fitTraj(); + auto dir0 = ftraj.direction(ftraj.t0()); + TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; + TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; + auto tofront = ktrk.extrapolate(fronttdir,calod1Front_); + auto toback = ktrk.extrapolate(backtdir,calod1Back_); + // record the standard tracker intersections + static const SurfaceId d1_front("EMC_Disk_1_Front"); + static const SurfaceId d1_back("EMC_Disk_1_Back"); + std::cout<<"extrapolating to calo d1"< bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index a13022aca8..16de69d4e1 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -115,6 +115,7 @@ namespace mu2e { fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom ToCaloD0 { Name("ToCaloD0"), Comment("Extrapolate tracks to the disk0 ends") }; + fhicl::Atom ToCaloD1 { Name("ToCaloD1"), Comment("Extrapolate tracks to the disk1 ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; }; From 7c5bb422e670963312826a44ceb2e77f0cbe84d0 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Mon, 1 Dec 2025 11:16:47 -0600 Subject: [PATCH 10/31] added --- CommonMC/src/MakeSurfaceSteps_module.cc | 7 ++++++- Mu2eKinKal/inc/KKExtrap.hh | 20 ++++++++++++-------- 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/CommonMC/src/MakeSurfaceSteps_module.cc b/CommonMC/src/MakeSurfaceSteps_module.cc index ed614802b0..8203758697 100644 --- a/CommonMC/src/MakeSurfaceSteps_module.cc +++ b/CommonMC/src/MakeSurfaceSteps_module.cc @@ -65,7 +65,12 @@ namespace mu2e { vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_SurfIn)] = SurfaceId("EMC_Disk_1_Front"); vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_0_SurfOut)] = SurfaceId("EMC_Disk_0_Back"); vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_SurfOut)] = SurfaceId("EMC_Disk_1_Back"); - } + + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_0_EdgeIn)] = SurfaceId("EMC_Disk_0_Inner"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_EdgeIn)] = SurfaceId("EMC_Disk_1_Inner"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_0_EdgeOut)] = SurfaceId("EMC_Disk_0_Outer"); + vdmap_[VirtualDetectorId(VirtualDetectorId::EMC_Disk_1_EdgeOut)] = SurfaceId("EMC_Disk_1_Outer"); +} void MakeSurfaceSteps::produce(art::Event& event) { std::cout<<"MakeSurfaceSteps ---> produce"< 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - auto tofront = ktrk.extrapolate(fronttdir,calod0Front_); - auto toback = ktrk.extrapolate(backtdir,calod0Back_); + std::cout<<"calod0Front_ "< 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - auto tofront = ktrk.extrapolate(fronttdir,calod1Front_); - auto toback = ktrk.extrapolate(backtdir,calod1Back_); + auto tocalofront = ktrk.extrapolate(fronttdir,calod1Front_); + auto tocaloback = ktrk.extrapolate(backtdir,calod1Back_); + std::cout<<"calod1Front_ "< Date: Mon, 1 Dec 2025 12:42:18 -0600 Subject: [PATCH 11/31] works --- KinKalGeom/src/Calo.cc | 2 +- Mu2eKinKal/inc/KKExtrap.hh | 15 +++++++++++---- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index f3c663e5f8..8282b1148b 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -25,7 +25,7 @@ namespace mu2e { EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667),719.)}, EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+2*192.295),719.)}, EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045),719.)}, - EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045 + +2*192.295 ),719.)} + EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045+2*192.295 ),719.)} {} } } diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 23bd832e95..4b20260589 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -47,6 +47,8 @@ namespace mu2e { KKMaterial const& kkmat_; AnnPtr tsdaptr_; DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; + DiskPtr calod0frontptr_, calod0backptr_; + DiskPtr calod1frontptr_, calod1backptr_; FruPtr opaptr_; bool backToTracker_, toOPA_, toTrackerEnds_, upstream_, toCaloD0_, toCaloD1_; ExtrapolateToZ TSDA_, trackerFront_, trackerBack_, calod0Front_, calod0Back_, calod1Front_, calod1Back_; // extrapolation predicate based on Z values @@ -67,6 +69,10 @@ namespace mu2e { trkfrontptr_(smap_.tracker().frontPtr()), trkmidptr_(smap_.tracker().middlePtr()), trkbackptr_(smap_.tracker().backPtr()), + calod0frontptr_(smap_.calo().EMC_Disk_0_FrontPtr()), + calod0backptr_(smap_.calo().EMC_Disk_0_BackPtr()), + calod1frontptr_(smap_.calo().EMC_Disk_1_FrontPtr()), + calod1backptr_(smap_.calo().EMC_Disk_1_BackPtr()), opaptr_(smap_.DS().outerProtonAbsorberPtr()), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), @@ -81,6 +87,7 @@ namespace mu2e { calod0Back_(maxdt_,btol_,smap_.calo().EMC_Disk_0_Back().center().Z(),debug_), calod1Front_(maxdt_,btol_,smap_.calo().EMC_Disk_1_Front().center().Z(),debug_), calod1Back_(maxdt_,btol_,smap_.calo().EMC_Disk_1_Back().center().Z(),debug_), + extrapIPA_(maxdt_,btol_,intertol_,smap_.DS().innerProtonAbsorberPtr(),debug_), extrapST_(maxdt_,btol_,intertol_,smap_.ST(),debug_) { @@ -183,14 +190,14 @@ namespace mu2e { // check the front piece first; that is usually correct // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto frontinter = KinKal::intersect(fhel,*trkfrontptr_,fhel.range(),intertol_,fronttdir); + auto frontinter = KinKal::intersect(fhel,*calod0frontptr_,fhel.range(),intertol_,fronttdir); if(frontinter.good()) ktrk.addIntersection(d0_front,frontinter); std::cout<<"to front "< Date: Mon, 1 Dec 2025 14:00:13 -0600 Subject: [PATCH 12/31] autoplay --- Mu2eKinKal/fcl/prolog.fcl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 0e396b5f63..9fbbd4c06a 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -400,6 +400,7 @@ Mu2eKinKal : { BackToTracker : false ToOPA : false ToCaloD0 : true + ToCaloD1 : true } CENTRALHELIX : { @@ -437,6 +438,8 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MinV : 1e-5 ToCRV : true + ToCaloD0 : true + ToCaloD1 : true } } From 23ab8d1abd4a46ef86f09aeb13f3d5fed42580e4 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 20 Jan 2026 09:16:59 -0600 Subject: [PATCH 13/31] updated Calo --- KinKalGeom/src/Calo.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index 8282b1148b..693627720f 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -22,10 +22,10 @@ namespace mu2e { EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),335,192.295)}, EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, - EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667),719.)}, - EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+2*192.295),719.)}, - EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045),719.)}, - EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045+2*192.295 ),719.)} + EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667-192.295),719.)}, + EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+192.295),719.)}, + EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045-192.295),719.)}, + EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045+192.295 ),719.)} {} } } From 2610df839b9af5b9f1741874b1b9d284764fa681 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 12 Mar 2026 13:12:05 -0500 Subject: [PATCH 14/31] removed --- .../src/InterdetectorTimeBootstrap_module.cc | 128 ------------------ mu2e.fcl | 5 - 2 files changed, 133 deletions(-) delete mode 100644 TrkReco/src/InterdetectorTimeBootstrap_module.cc delete mode 100644 mu2e.fcl diff --git a/TrkReco/src/InterdetectorTimeBootstrap_module.cc b/TrkReco/src/InterdetectorTimeBootstrap_module.cc deleted file mode 100644 index 02413fae0d..0000000000 --- a/TrkReco/src/InterdetectorTimeBootstrap_module.cc +++ /dev/null @@ -1,128 +0,0 @@ -/* An art analyzer module for the initial calibration of interdetector timing. - -#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 "Offline/RecoDataProducts/inc/KalSeed.hh" -#include "Offline/RecoDataProducts/inc/CaloCluster.hh" -#include "Offline/RecoDataProducts/inc/CrvCoincidenceCluster.hh" -#include "Offline/GeometryService/inc/GeometryService.hh" -#include "Offline/GeometryService/inc/GeomHandle.hh" -#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" -#include "Offline/CosmicRayShieldGeom/inc/CosmicRayShield.hh" -#include "Offline/TrackerGeom/inc/Tracker.hh" -#include "CLHEP/Vector/ThreeVector.h" -#include "TTree.h" -#include "art_root_io/TFileService.h" -#include -#include - -using namespace std; - -namespace mu2e { - - class InterdetectorTimeBootstrap : public art::EDAnalyzer { - - public: - - struct Config { - using Name=fhicl::Name; - using Comment=fhicl::Comment; - fhicl::Atom diagLevel{Name("diagLevel"),Comment("diag level"),0}; - fhicl::Atom kalTag{Name("kalseedModuleLabel"), Comment("kal seed info")}; - fhicl::Atom caloTag{Name("caloClusterModuleLabel"), Comment("calo cluster info")}; - fhicl::Atom crvTag{Name("crvClusterModuleLabel"), Comment("crv info")}; - fhicl::Atom maxDeltaR{Name("maxDeltaR"),Comment("maxDeltaR"),1000.}; - fhicl::Atom maxDeltaT{Name("maxDeltaT"),Comment("maxDeltaT"),1000.}; - }; - typedef art::EDAnalyzer::Table Parameters; - - explicit InterdetectorTimeBootstrap(const Parameters& conf); - virtual ~InterdetectorTimeBootstrap() {} - - - virtual void beginJob(); - virtual void endJob(); - virtual void analyze(const art::Event& e) override; - - private: - Config _conf; - int _diagLevel; - art::InputTag _kalseedToken; - art::InputTag _caloClusterToken; - art::InputTag _crvClusterToken; - double _maxDeltaR; - double _maxDeltaT; - TTree* _timingTree; - double _deltaT_trkcal; - double _deltaT_trkcrv; -}; - - - InterdetectorTimeBootstrap::InterdetectorTimeBootstrap(const Parameters& conf): - art::EDAnalyzer(conf), - _diagLevel(conf().diagLevel()), - _kalseedToken(conf().kalTag()), - _caloClusterToken(conf().caloTag()), - _crvClusterToken(conf().crvTag()), - _maxDeltaR(conf().maxDeltaR()), - _maxDeltaT(conf().maxDeltaT()) - {} - - //InterdetectorTimeBootstrap::~InterdetectorTimeBootstrap() {} - - void InterdetectorTimeBootstrap::beginJob() { - art::ServiceHandle tfs; - _timingTree = tfs->make("timingTree", "Timing Differences"); - _timingTree->Branch("deltaT_trkcal", &_deltaT_trkcal, "deltaT_trkcal/F"); - _timingTree->Branch("deltaT_trkcrv", &_deltaT_trkcrv, "deltaT_trkcrv/F"); - } - - void InterdetectorTimeBootstrap::analyze(const art::Event& event) { - // Access event data products for all three main detectors: - art::Handle tracksHandle; - event.getByLabel(_kalseedToken, tracksHandle); - const KalSeedCollection& tracks = *tracksHandle; - - art::Handle clustersHandle; - event.getByLabel(_caloClusterToken, clustersHandle); - const CaloClusterCollection& clusters = *clustersHandle; - - art::Handle crvHandle; - event.getByLabel(_crvClusterToken, crvHandle); - const CrvCoincidenceClusterCollection& CRVclusters = *crvHandle; - - // Get geometry handles - //GeomHandle calo; - //GeomHandle tracker; - //GeomHandle crv; - - for (const auto& track : tracks) { - // Use the track' parameters for extrapolation - CLHEP::Hep3Vector track_pos_at_calo; // Position of track extrapolated to the calo face - - // Use specific Mu2e geometry functions to find the intersection - bool intersects = true; // Use geometry to check intersection - - if (intersects) { - for (const auto& cluster : clusters) { - const CLHEP::Hep3Vector& cluster_pos = cluster.cog3Vector(); - // Simple geometric match: check proximity - if ((track_pos_at_calo - cluster_pos).mag() < _maxDeltaR) { - // Calculate time difference - double time_calo = cluster.time(); - // The track time needs to be adjusted for the time of flight from the track origin to the calo - double time_track_extrapolated = track.time() + (track_pos_at_calo.mag() / 299.792458); // TOF in ns - _deltaT_trkcal = time_calo - time_track_extrapolated; - - // Fill the NTuple for later analysis (e.g., histogramming the mean) - _timingTree->Fill(); - } - } - } - } - } -} // namespace mu2e - -DEFINE_ART_MODULE(mu2e::InterdetectorTimeBootstrap)*/ diff --git a/mu2e.fcl b/mu2e.fcl deleted file mode 100644 index 2030ff9d26..0000000000 --- a/mu2e.fcl +++ /dev/null @@ -1,5 +0,0 @@ -#include "Production/JobConfig/recoMC/OnSpill.fcl" -services.DbService.purpose: MDC2025_best -services.DbService.version: v1_4 -outputs.LoopHelixOutput.fileName: "mcs.mu2e.test.v1.sequencer.art" -//services.GeometryService.inputFile: "Offline/Mu2eG4/geom/geom_run1_a.txt" From 7b108aee5c9d80a44f9b6b5f522b868892a65f1f Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 9 Apr 2026 13:32:47 -0500 Subject: [PATCH 15/31] surface id --- DataProducts/src/SurfaceId.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 68c0bc058a..484a74af4f 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -35,7 +35,6 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Outer, "ST_Outer"), std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), - std::make_pair(SurfaceIdEnum::EMC_Disk_0_Outer, "EMC_Disk_0_Outer"), std::make_pair(SurfaceIdEnum::EMC_Disk_0_Inner, "EMC_Disk_0_Inner"), std::make_pair(SurfaceIdEnum::EMC_Disk_1_Inner, "EMC_Disk_1_Inner"), From 67899a542648c17a9c1a7bd64db27fe1356f629e Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 23 Apr 2026 11:24:52 -0500 Subject: [PATCH 16/31] works with geom service --- KinKalGeom/inc/Calo.hh | 12 ++++++++++-- KinKalGeom/src/Calo.cc | 31 +++++++++++-------------------- Mu2eKinKal/inc/KKExtrap.hh | 10 +++++----- 3 files changed, 26 insertions(+), 27 deletions(-) diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh index 91198a478d..571e75bbbf 100644 --- a/KinKalGeom/inc/Calo.hh +++ b/KinKalGeom/inc/Calo.hh @@ -15,8 +15,9 @@ namespace mu2e { public: using CylPtr = std::shared_ptr; using DiskPtr = std::shared_ptr; - // default constructor with nominal geometry - Calo(); + // constructor with geometry parameters from Calorimeter service + Calo(double z0, double z1, double r0_inner, double r0_outer, double r1_inner, double r1_outer, + double z0_front, double z0_back, double z1_front, double z1_back); // return by reference auto const& EMC_Disk_0_Outer() const { return *EMC_Disk_0_Outer_;} auto const& EMC_Disk_0_Inner() const { return *EMC_Disk_0_Inner_;} @@ -38,9 +39,16 @@ namespace mu2e { auto const& EMC_Disk_0_BackPtr() const { return EMC_Disk_0_Back_;} auto const& EMC_Disk_1_BackPtr() const { return EMC_Disk_1_Back_;} + // accessors for local Z positions (relative to tracker center) + double EMC_Disk_0_Front_Z() const { return z0_front_; } + double EMC_Disk_0_Back_Z() const { return z0_back_; } + double EMC_Disk_1_Front_Z() const { return z1_front_; } + double EMC_Disk_1_Back_Z() const { return z1_back_; } + private: CylPtr EMC_Disk_0_Inner_, EMC_Disk_0_Outer_ , EMC_Disk_1_Inner_, EMC_Disk_1_Outer_; DiskPtr EMC_Disk_0_Front_, EMC_Disk_0_Back_, EMC_Disk_1_Front_, EMC_Disk_1_Back_; + double z0_front_, z0_back_, z1_front_, z1_back_; // local Z positions }; } } diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index cd4e9de84b..ccd0871edf 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -5,27 +5,18 @@ namespace mu2e { using KinKal::Cylinder; using KinKal::Disk; - // D0 = 11842 -10175.0 = 1667 - // D1 = 13220 -10175.0 = 3045 - /* - double calorimeter.caloDiskRadiusIn = 335; - double calorimeter.caloDiskRadiusOut = 719; - half length = 192.295 + Calo::Calo(double z0, double z1, double r0_inner, double r0_outer, double r1_inner, double r1_outer, + double z0_front, double z0_back, double z1_front, double z1_back) : + EMC_Disk_0_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z0),r0_inner, 0.5*(z0_back-z0_front))}, + EMC_Disk_0_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z0),r0_outer, 0.5*(z0_back-z0_front))}, + EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z1),r1_inner, 0.5*(z1_back-z1_front))}, + EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z1),r1_outer, 0.5*(z1_back-z1_front))}, - - Cylinder(VEC3 const& axis, VEC3 const& center, double radius, double halflen ); - Disk(VEC3 const& norm,VEC3 const& udir, VEC3 const& center, double radius) - */ - Calo::Calo() : - EMC_Disk_0_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),335,192.295)}, - EMC_Disk_0_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,1667),719,192.295)}, - EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),335,192.295)}, - EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,3045),719,192.295)}, - - EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667-192.295),719.)}, - EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,1667+192.295),719.)}, - EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045-192.295),719.)}, - EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,3045+192.295 ),719.)} + EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z0_front),r0_outer)}, + EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z0_back),r0_outer)}, + EMC_Disk_1_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z1_front),r1_outer)}, + EMC_Disk_1_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z1_back),r1_outer)}, + z0_front_(z0_front), z0_back_(z0_back), z1_front_(z1_front), z1_back_(z1_back) {} } } diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 097ad697b5..c9a7633eab 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -110,11 +110,11 @@ namespace mu2e { const_cast(extrapIPA_) = ExtrapolateIPA(maxdt_,btol_,intertol_,kkg.DS().innerProtonAbsorberPtr(),debug_); const_cast(extrapST_) = ExtrapolateST(maxdt_,btol_,intertol_,kkg.ST(),debug_); - // calo predicates - const_cast(calod0Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Front().center().Z(),debug_); - const_cast(calod0Back_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Back().center().Z(),debug_); - const_cast(calod1Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_1_Front().center().Z(),debug_); - const_cast(calod1Back_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_1_Back().center().Z(),debug_); + // calo predicates - use local Z values stored in Calo class + const_cast(calod0Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Front_Z(),debug_); + const_cast(calod0Back_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Back_Z(),debug_); + const_cast(calod1Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_1_Front_Z(),debug_); + const_cast(calod1Back_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_1_Back_Z(),debug_); if(debug_ > 0)std::cout << "IPA limits z " << extrapIPA_.zmin() << " " << extrapIPA_.zmax() << std::endl; if(debug_ > 0)std::cout << "ST limits z " << extrapST_.zmin() << " " << extrapST_.zmax() << " r " << extrapST_.rmin() << " " << extrapST_.rmax() << std::endl; From 08e83c4724613c2e23f13d11ae07c3f3d25667dd Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 23 Apr 2026 11:30:23 -0500 Subject: [PATCH 17/31] fixed the calo front --- Mu2eKinKal/inc/KKExtrap.hh | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index c9a7633eab..0f32496705 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -229,10 +229,9 @@ namespace mu2e { static const SurfaceId d0_back("EMC_Disk_0_Back"); std::cout<<"extrapolating to calo d0"< Date: Thu, 23 Apr 2026 13:55:20 -0500 Subject: [PATCH 18/31] seems to work --- Mu2eKinKal/fcl/prolog.fcl | 11 +++++-- Mu2eKinKal/inc/KKExtrap.hh | 36 ++++++++++++--------- Mu2eKinKal/src/KinematicLineFit_module.cc | 38 +++++++++-------------- 3 files changed, 43 insertions(+), 42 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 50d35aeb7f..6e87770a5d 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -439,11 +439,16 @@ Mu2eKinKal : { } LINEEXTRAPOLATION : { - MaxDt : 200.0 # (ns) - MinV : 1e-5 - ToCRV : true + MaxDt : 200.0 + BCorrTolerance : 1e-5 + ToTrackerEnds : true + Upstream : true + BackToTracker : false + ToTrackerEnds : false + ToOPA : false ToCaloD0 : true ToCaloD1 : true + IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 0f32496705..bcbcbb362f 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -212,21 +212,21 @@ namespace mu2e { template void KKExtrap::toCaloD0(KKTrack& ktrk) const { initGeometry(); + GeomHandle kkg_h; + auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - { - GeomHandle kkg_h; - auto const& kkg = *kkg_h; - std::cout<<"calod0Front_ "< void KKExtrap::toCaloD1(KKTrack& ktrk) const { initGeometry(); + GeomHandle kkg_h; + auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; auto tocalofront = ktrk.extrapolate(fronttdir,calod1Front_); auto tocaloback = ktrk.extrapolate(backtdir,calod1Back_); - { - GeomHandle kkg_h; - auto const& kkg = *kkg_h; - std::cout<<"calod1Front_ "< bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index a6aa9f5a9f..ece632de6e 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -52,6 +52,7 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateTCRV.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // root #include "TH1F.h" #include "TTree.h" @@ -119,21 +120,13 @@ namespace mu2e { fhicl::Atom sampleTBuff { Name("SampleTimeBuffer"), Comment("Time buffer for sample intersections (nsec)") }; }; - // Extrapolation configuration - struct KKExtrapConfig { - fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit") }; - fhicl::Atom MinV { Name("MinV"), Comment("Minimum Y vel to extrapolate a fit") }; - fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV") }; - fhicl::Atom Debug { Name("Debug"), Comment("Debug level"), 0 }; - }; - struct GlobalConfig { fhicl::Table modSettings { Name("ModuleSettings") }; fhicl::Table mu2eSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; fhicl::Table matSettings { Name("MaterialSettings") }; - fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; + fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; }; @@ -171,10 +164,12 @@ namespace mu2e { double intertol_; // surface intersection tolerance (mm) double sampletbuff_; // simple time buffer; replace this with extrapolation TODO bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface + SurfaceIdCollection ssids_; // surface IDs to sample KinKalGeom::SurfacePairCollection sample_; // surfaces to sample the fit bool extrapolate_, toCRV_; ExtrapolateTCRV TCRV_; // extrapolation predicate based on Z values double tcrvthick_ = 0.1056; // st foil thickness: should come from geometry service TODO + std::unique_ptr extrap_; // calorimeter and other extrapolations Config config_; // initial fit configuration object Config exconfig_; // extension configuration object }; @@ -221,26 +216,15 @@ namespace mu2e { }else{ throw cet::exception("RECO")<<"mu2e::KinematicLineFit: Parameter constraint configuration error"<< endl; } - SurfaceIdCollection ssids; + // store surface IDs to be resolved when geometry is available for(auto const& sidname : settings().modSettings().sampleSurfaces()) { - ssids.push_back(SurfaceId(sidname,-1)); // match all elements + ssids_.push_back(SurfaceId(sidname,-1)); // match all elements } - // translate the sample and extend surface names to actual surfaces using the KinKalGeom. This should come from the - // geometry service eventually, TODO - GeomHandle kkg_h; - auto const& kkg = *kkg_h; - kkg.surfaces(ssids,sample_); // configure extrapolation if(settings().Extrapolation()){ extrapolate_ = true; - toCRV_ = settings().Extrapolation()->ToCRV(); - // global configs - double maxdt = settings().Extrapolation()->MaxDt(); - double btol = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension - double minv = settings().Extrapolation()->MinV(); - int debug = settings().Extrapolation()->Debug(); - // extrapolate to the front of the tracker - TCRV_ = ExtrapolateTCRV(maxdt,btol,intertol_,minv,kkg.TCRV(),debug); + // create KKExtrap for calorimeter and upstream extrapolations + extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); } if(print_ > 0) std::cout << config_; @@ -259,6 +243,10 @@ namespace mu2e { GeomHandle bfmgr; GeomHandle det; kkbf_ = std::make_unique(*bfmgr,*det); + // now that geometry is available, resolve sample surface IDs to actual surfaces + GeomHandle kkg_h; + auto const& kkg = *kkg_h; + kkg.surfaces(ssids_,sample_); } void KinematicLineFit::produce(art::Event& event ) { @@ -435,6 +423,8 @@ namespace mu2e { ktrk.addTCRVXing(crvxingptr,tdir); } } while(hadintersection); + // apply additional extrapolations (calorimeter, tracker ends, IPA, ST, etc) via KKExtrap if configured + if(extrap_) extrap_->extrapolate(ktrk); } } DEFINE_ART_MODULE(mu2e::KinematicLineFit) From 7564a99c38c002c90fbcacce0de49c71d6bd027e Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Fri, 24 Apr 2026 11:20:04 -0500 Subject: [PATCH 19/31] calo intersections in correct location - calo mother --- GeometryService/inc/KinKalGeomMaker.hh | 1 + GeometryService/src/KinKalGeomMaker.cc | 95 +++++++++++++++++++++++--- KinKalGeom/src/Calo.cc | 8 +-- Mu2eKinKal/inc/KKExtrap.hh | 7 +- 4 files changed, 95 insertions(+), 16 deletions(-) diff --git a/GeometryService/inc/KinKalGeomMaker.hh b/GeometryService/inc/KinKalGeomMaker.hh index d0adbe8dbc..77209b50ee 100644 --- a/GeometryService/inc/KinKalGeomMaker.hh +++ b/GeometryService/inc/KinKalGeomMaker.hh @@ -13,6 +13,7 @@ namespace mu2e { void makeTracker(); void makeDS(); void makeTarget(); + void makeCalo(); void makeTCRV(); std::unique_ptr kkg_; }; diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index 5c955ef72e..aeec645dc8 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -10,14 +10,17 @@ #include "Offline/KinKalGeom/inc/DetectorSolenoid.hh" #include "Offline/KinKalGeom/inc/StoppingTarget.hh" #include "Offline/KinKalGeom/inc/TestCRV.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" #include "Offline/BeamlineGeom/inc/Beamline.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/DetectorSolenoidGeom/inc/DetectorSolenoid.hh" +#include "Offline/CalorimeterGeom/inc/Calorimeter.hh" namespace mu2e { using KinKal::VEC3; using KinKal::Cylinder; - using KinKal::Disk; + // NOTE: Do NOT use `using KinKal::Disk` because CalorimeterGeom also defines mu2e::Disk + // Use explicit KinKal::Disk instead using KinKal::Surface; using KinKal::Frustrum; using KinKal::Annulus; @@ -36,9 +39,81 @@ namespace mu2e { makeDS(); makeTarget(); makeTCRV(); + makeCalo(); return kkg_; } + void KinKalGeomMaker::makeCalo() { + // construct calorimeter geometry from GeometryService + auto const& calo_det = *(GeomHandle()); + auto const& tracker = *(GeomHandle()); + + // Extract geometry from the first two disks (D0 and D1) + if (calo_det.nDisks() < 2) { + throw cet::exception("KinKalGeomMaker") + << "Calorimeter has fewer than 2 disks; cannot construct KinKal geometry\n"; + } + + auto const& disk0 = calo_det.disk(0); + auto const& disk1 = calo_det.disk(1); + auto const& geom0 = disk0.geomInfo(); + auto const& geom1 = disk1.geomInfo(); + + // Extract Z positions and radii from geometry (in global coordinates) + double z0_global = geom0.origin().z(); + double z1_global = geom1.origin().z(); + double r0_inner = geom0.innerEnvelopeR(); + double r0_outer = geom0.outerEnvelopeR(); + double r1_inner = geom1.innerEnvelopeR(); + double r1_outer = geom1.outerEnvelopeR(); + + // Use true disk geometry: compute symmetric front/back from disk center + // geomInfo().size().z() gives the full disk extent (front plate + crystals + back plate) + // This matches how VirtualDetectors are defined in VirtualDetectorMaker + double diskHalfLen_0 = 0.5 * geom0.size().z(); + double diskHalfLen_1 = 0.5 * geom1.size().z(); + + double z0_front_global = z0_global - diskHalfLen_0; + double z0_back_global = z0_global + diskHalfLen_0; + double z1_front_global = z1_global - diskHalfLen_1; + double z1_back_global = z1_global + diskHalfLen_1; + + // Convert from global to local coordinates (relative to tracker center) + double tracker_z0 = tracker.g4Tracker()->z0(); + double z0 = z0_global - tracker_z0; + double z1 = z1_global - tracker_z0; + double z0_front = z0_front_global - tracker_z0; + double z0_back = z0_back_global - tracker_z0; + double z1_front = z1_front_global - tracker_z0; + double z1_back = z1_back_global - tracker_z0; + + std::cout << "KinKalGeomMaker::makeCalo() DEBUG:" << std::endl; + std::cout << " Tracker z0 offset: " << tracker_z0 << std::endl; + std::cout << " Disk 0 (global): z=" << z0_global << " z_front=" << z0_front_global << " z_back=" << z0_back_global << std::endl; + std::cout << " Disk 0 (local): z=" << z0 << " z_front=" << z0_front << " z_back=" << z0_back << std::endl; + std::cout << " Disk 0 thickness check: front_to_back = " << (z0_back_global - z0_front_global) << " mm (should be positive)" << std::endl; + std::cout << " Disk 0 center check: center z = " << z0_global << ", midpoint of front/back = " << (0.5*(z0_front_global + z0_back_global)) << std::endl; + std::cout << " Disk 1 (global): z=" << z1_global << " z_front=" << z1_front_global << " z_back=" << z1_back_global << std::endl; + std::cout << " Disk 1 (local): z=" << z1 << " z_front=" << z1_front << " z_back=" << z1_back << std::endl; + std::cout << " Disk 1 thickness check: front_to_back = " << (z1_back_global - z1_front_global) << " mm (should be positive)" << std::endl; + std::cout << " Disk 1 center check: center z = " << z1_global << ", midpoint of front/back = " << (0.5*(z1_front_global + z1_back_global)) << std::endl; + std::cout << " Radii: r0_inner=" << r0_inner << " r0_outer=" << r0_outer << " r1_inner=" << r1_inner << " r1_outer=" << r1_outer << std::endl; + + // Construct Calo with geometry in local coordinates + kkg_->calo_ = std::make_unique(z0, z1, r0_inner, r0_outer, r1_inner, r1_outer, z0_front, z0_back, z1_front, z1_back); + auto const& calo = *kkg_->calo_; + + // Add all calo surfaces to the map + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Outer),std::static_pointer_cast(calo.EMC_Disk_0_OuterPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Inner),std::static_pointer_cast(calo.EMC_Disk_0_InnerPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Inner),std::static_pointer_cast(calo.EMC_Disk_1_InnerPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Outer),std::static_pointer_cast(calo.EMC_Disk_1_OuterPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Front),std::static_pointer_cast(calo.EMC_Disk_0_FrontPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Front),std::static_pointer_cast(calo.EMC_Disk_1_FrontPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Back),std::static_pointer_cast(calo.EMC_Disk_0_BackPtr()))); + kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Back),std::static_pointer_cast(calo.EMC_Disk_1_BackPtr()))); + } + void KinKalGeomMaker::makeTracker() { // surfaces need to match with virtual detectors. The following is extracted from VirtualDetectorMaker and needs to be updated if that changes. // Note that these are placed at the center of the VDs, which have half-thickness of 0.01mm. Since the VD hits are recorded where the SimParticle @@ -65,9 +140,9 @@ namespace mu2e { auto outer = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,zMidLocal),orvd,halfLen); auto inner = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,zMidLocal),irvd,halfLen); // expand the disk radii to the DS - auto front = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zFrontLocal),irds); - auto mid = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zMidLocal),irds); - auto back = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zBackLocal),irds); + auto front = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zFrontLocal),irds); + auto mid = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zMidLocal),irds); + auto back = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,zBackLocal),irds); // add all these to the map kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TT_Front),std::static_pointer_cast(front))); kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TT_Mid),std::static_pointer_cast(mid))); @@ -82,11 +157,11 @@ namespace mu2e { // currently use hard-coded geometry auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),950,5450); auto outer= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),1328,5450); // bounding surfaces - auto front= std::make_shared(outer->frontDisk()); - auto back= std::make_shared(outer->backDisk()); + auto front= std::make_shared(outer->frontDisk()); + auto back= std::make_shared(outer->backDisk()); auto ipa= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-2770),300.0,500.0); - auto ipafront= std::make_shared(ipa->frontDisk()); - auto ipaback= std::make_shared(ipa->backDisk()); + auto ipafront= std::make_shared(ipa->frontDisk()); + auto ipaback= std::make_shared(ipa->backDisk()); auto opa= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-3766),454.0,728.4,2125.0); // inner surface auto tsda= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,-5967),235.0,525.0); // back surface @@ -107,8 +182,8 @@ namespace mu2e { // currently use hard-coded geometry auto outer = std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-4300),75,400.0); auto inner= std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-4300),21.5,400.0); - auto front= std::make_shared(outer->frontDisk()); - auto back= std::make_shared(outer->backDisk()); + auto front= std::make_shared(outer->frontDisk()); + auto back= std::make_shared(outer->backDisk()); double startz = -4700; double endz = -3900; double dz = (endz-startz)/36.0; diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc index ccd0871edf..6c7152162a 100644 --- a/KinKalGeom/src/Calo.cc +++ b/KinKalGeom/src/Calo.cc @@ -7,10 +7,10 @@ namespace mu2e { Calo::Calo(double z0, double z1, double r0_inner, double r0_outer, double r1_inner, double r1_outer, double z0_front, double z0_back, double z1_front, double z1_back) : - EMC_Disk_0_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z0),r0_inner, 0.5*(z0_back-z0_front))}, - EMC_Disk_0_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z0),r0_outer, 0.5*(z0_back-z0_front))}, - EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z1),r1_inner, 0.5*(z1_back-z1_front))}, - EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,z1),r1_outer, 0.5*(z1_back-z1_front))}, + EMC_Disk_0_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,0.5*(z0_front+z0_back)),r0_inner, 0.5*(z0_back-z0_front))}, + EMC_Disk_0_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,0.5*(z0_front+z0_back)),r0_outer, 0.5*(z0_back-z0_front))}, + EMC_Disk_1_Inner_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,0.5*(z1_front+z1_back)),r1_inner, 0.5*(z1_back-z1_front))}, + EMC_Disk_1_Outer_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,0.5*(z1_front+z1_back)),r1_outer, 0.5*(z1_back-z1_front))}, EMC_Disk_0_Front_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z0_front),r0_outer)}, EMC_Disk_0_Back_ { std::make_shared(VEC3(0.0,0.0,1.0),VEC3(1.0,0.0,0.0),VEC3(0.0,0.0,z0_back),r0_outer)}, diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index bcbcbb362f..0c8b371913 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -218,8 +218,11 @@ namespace mu2e { auto dir0 = ftraj.direction(ftraj.t0()); TimeDir fronttdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - std::cout<<"calod0Front_ "< 0){ + std::cout<<"toCaloD0 DEBUG:"< fitSettings { Name("FitSettings") }; fhicl::Table extSettings { Name("ExtensionSettings") }; fhicl::Table matSettings { Name("MaterialSettings") }; -// helix module specific config + fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; }; class CentralHelixFit : public art::EDProducer { @@ -142,6 +144,7 @@ namespace mu2e { protected: bool goodFit(KKTRK const& ktrk) const; void sampleFit(KKTRK& kktrk) const; + void extrapolate(KKTRK& ktrk) const; TrkFitFlag fitflag_; // parameter-specific functions that need to be overridden in subclasses // data payload @@ -172,6 +175,8 @@ namespace mu2e { bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface KinKalGeom::SurfacePairCollection sample_; // surfaces to sample the fit std::array paramconstraints_; + bool extrapolate_; + std::unique_ptr extrap_; // calorimeter and other extrapolations }; CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings}, @@ -194,7 +199,8 @@ namespace mu2e { useFitCharge_(settings().modSettings().useFitCharge()), minCenterRho_(settings().modSettings().minCenterRho()), sampleinrange_(settings().modSettings().sampleInRange()), - sampleinbounds_(settings().modSettings().sampleInBounds()) + sampleinbounds_(settings().modSettings().sampleInBounds()), + extrapolate_(false) { // collection handling for(const auto& cseedtag : settings().modSettings().seedCollections()) { cseedCols_.emplace_back(consumes(cseedtag)); } @@ -226,6 +232,12 @@ namespace mu2e { // geometry service eventually, TODO KinKalGeom smap; smap.surfaces(ssids,sample_); + // configure extrapolation + if(settings().Extrapolation()){ + extrapolate_ = true; + // create KKExtrap for calorimeter and upstream extrapolations + extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); + } } void CentralHelixFit::beginRun(art::Run& run) { @@ -339,6 +351,8 @@ namespace mu2e { } if(print_>1)kktrk->printFit(std::cout,print_); + // extrapolate as required + if(goodfit && extrapolate_) extrapolate(*kktrk); if(goodfit || saveall_){ TrkFitFlag fitflag;//(hptr->status()); fitflag.merge(fitflag_); @@ -423,6 +437,11 @@ namespace mu2e { } } + void CentralHelixFit::extrapolate(KKTRK& ktrk) const { + // apply extrapolations (calorimeter, tracker ends, IPA, ST, etc) via KKExtrap if configured + if(extrap_) extrap_->extrapolate(ktrk); + } + } // mu2e using mu2e::CentralHelixFit; DEFINE_ART_MODULE(CentralHelixFit) From 28174e8dbf22fed99ec0267811275a0d28ddcd2f Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Fri, 24 Apr 2026 13:55:03 -0500 Subject: [PATCH 21/31] added code --- GeometryService/src/KinKalGeomMaker.cc | 7 ------- Mu2eKinKal/inc/KKExtrap.hh | 3 --- Mu2eKinKal/src/CentralHelixFit_module.cc | 4 ++-- 3 files changed, 2 insertions(+), 12 deletions(-) diff --git a/GeometryService/src/KinKalGeomMaker.cc b/GeometryService/src/KinKalGeomMaker.cc index 63cdfd3e56..b6cfdf5fe1 100644 --- a/GeometryService/src/KinKalGeomMaker.cc +++ b/GeometryService/src/KinKalGeomMaker.cc @@ -49,11 +49,6 @@ namespace mu2e { auto const& tracker = *(GeomHandle()); // Extract geometry from the first two disks (D0 and D1) - if (calo_det.nDisks() < 2) { - throw cet::exception("KinKalGeomMaker") - << "Calorimeter has fewer than 2 disks; cannot construct KinKal geometry\n"; - } - auto const& disk0 = calo_det.disk(0); auto const& disk1 = calo_det.disk(1); auto const& geom0 = disk0.geomInfo(); @@ -68,8 +63,6 @@ namespace mu2e { double r1_outer = geom1.outerEnvelopeR(); // Use true disk geometry: compute symmetric front/back from disk center - // geomInfo().size().z() gives the full disk extent (front plate + crystals + back plate) - // This matches how VirtualDetectors are defined in VirtualDetectorMaker double diskHalfLen_0 = 0.5 * geom0.size().z(); double diskHalfLen_1 = 0.5 * geom1.size().z(); diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index dfe405b456..b1b9f5fd8b 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -189,7 +189,6 @@ namespace mu2e { auto midinter = KinKal::intersect(ftraj,kkg.tracker().middle(),ftraj.range(),intertol_); if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); if(tofront){ - // check the front piece first; that is usually correct // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); auto frontinter = KinKal::intersect(fhel,kkg.tracker().front(),fhel.range(),intertol_,fronttdir); @@ -266,13 +265,11 @@ namespace mu2e { static const SurfaceId d1_outer("EMC_Disk_1_Outer"); if(tocalofront){ - // use full trajectory range for intersection calculation TimeRange frange = ftraj.range(); auto frontinter = KinKal::intersect(ftraj,*calod1frontptr_,frange,intertol_,fronttdir); if(frontinter.good()) ktrk.addIntersection(d1_front,frontinter); } if(tocaloback){ - // start from the middle TimeRange brange = ftraj.range(); auto backinter = KinKal::intersect(ftraj,*calod1backptr_,brange,intertol_,backtdir); if(backinter.good())ktrk.addIntersection(d1_back,backinter); diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index 07ae94e7c5..053b24f27a 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -176,7 +176,7 @@ namespace mu2e { KinKalGeom::SurfacePairCollection sample_; // surfaces to sample the fit std::array paramconstraints_; bool extrapolate_; - std::unique_ptr extrap_; // calorimeter and other extrapolations + std::unique_ptr extrap_; // extrapolations }; CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings}, @@ -235,7 +235,7 @@ namespace mu2e { // configure extrapolation if(settings().Extrapolation()){ extrapolate_ = true; - // create KKExtrap for calorimeter and upstream extrapolations + // create KKExtrap extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); } } From fedf7fe610bc82fdb048a9d5b9e2d075b34ce04a Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 28 Apr 2026 08:36:05 -0500 Subject: [PATCH 22/31] merged --- Mu2eKinKal/src/KinematicLineFit_module.cc | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index e428a64c3b..342ca85a8d 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -229,6 +229,12 @@ namespace mu2e { extrapolate_ = true; // create KKExtrap for calorimeter and upstream extrapolations extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); + toCRV_ = settings().Extrapolation()->ToCRV(); + // global configs + maxdt_ = settings().Extrapolation()->MaxDt(); + btol_ = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension + minv_ = settings().Extrapolation()->MinV(); + extrapdebug_ = settings().Extrapolation()->Debug(); } @@ -248,6 +254,14 @@ namespace mu2e { GeomHandle bfmgr; GeomHandle det; kkbf_ = std::make_unique(*bfmgr,*det); + // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService + GeomHandle kkg_h; + auto const& kkg = *kkg_h; + kkg.surfaces(ssids_,surfacess_to_sample_); + // now that geometry is available, resolve sample surface IDs to actual surfaces + GeomHandle kkg_h; + auto const& kkg = *kkg_h; + kkg.surfaces(ssids_,sample_); } void KinematicLineFit::produce(art::Event& event ) { From 5aba2839bfa54c0bf456cf2db8ed09ab8d769780 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 28 Apr 2026 08:39:21 -0500 Subject: [PATCH 23/31] merged --- CommonMC/src/MakeSurfaceSteps_module.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/CommonMC/src/MakeSurfaceSteps_module.cc b/CommonMC/src/MakeSurfaceSteps_module.cc index 8203758697..f8e07bd676 100644 --- a/CommonMC/src/MakeSurfaceSteps_module.cc +++ b/CommonMC/src/MakeSurfaceSteps_module.cc @@ -73,7 +73,6 @@ namespace mu2e { } void MakeSurfaceSteps::produce(art::Event& event) { - std::cout<<"MakeSurfaceSteps ---> produce"< det; // create output std::unique_ptr ssc(new SurfaceStepCollection); @@ -83,7 +82,7 @@ namespace mu2e { for(auto const& vdspmc : vdspmccol) { // only some VDs are kept auto isid = vdmap_.find(vdspmc.virtualDetectorId()); - std::cout<<" VID "< 0) std::cout<<" VID "<emplace_back(isid->second,vdspmc,det); // no aggregation of VD hits } auto nvdsteps = ssc->size(); From e1d169e12a5d3d7deceda9c62729a5c3057df556 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 28 Apr 2026 09:06:57 -0500 Subject: [PATCH 24/31] builds with merge --- Mu2eKinKal/src/CentralHelixFit_module.cc | 8 ++++---- Mu2eKinKal/src/KinematicLineFit_module.cc | 14 +++----------- 2 files changed, 7 insertions(+), 15 deletions(-) diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index 28c5d3fe40..0a90870f6d 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -175,7 +175,7 @@ namespace mu2e { double minCenterRho_; // min center distance to z axis bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface SurfaceIdCollection ssids_; - KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit + KinKalGeom::SurfacePairCollection surfaces_to_sample_; // surfaces to sample the fit std::array paramconstraints_; bool extrapolate_; std::unique_ptr extrap_; // extrapolations @@ -233,7 +233,7 @@ namespace mu2e { // translate the sample and extend surface names to actual surfaces using the KinKalGeom. This should come from the // geometry service eventually, TODO KinKalGeom smap; - smap.surfaces(ssids,sample_); + smap.surfaces(ssids_,surfaces_to_sample_); // configure extrapolation if(settings().Extrapolation()){ extrapolate_ = true; @@ -257,7 +257,7 @@ namespace mu2e { // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService GeomHandle kkg_h; auto const& kkg = *kkg_h; - kkg.surfaces(ssids_,surfacess_to_sample_); + kkg.surfaces(ssids_,surfaces_to_sample_); } void CentralHelixFit::produce(art::Event& event ) { @@ -419,7 +419,7 @@ namespace mu2e { void CentralHelixFit::sampleFit(KKTRK& kktrk) const { auto const& ftraj = kktrk.fitTraj(); double tbeg = ftraj.range().begin(); - for(auto const& surf : surfacess_to_sample_){ + for(auto const& surf : surfaces_to_sample_){ // search for intersections with each surface from the begining double tstart = tbeg - sampletbuff_; bool goodinter(true); diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 342ca85a8d..2dd573ba10 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -169,8 +169,6 @@ namespace mu2e { KinKalGeom::SurfacePairCollection sample_; // surfaces to sample the fit bool extrapolate_, toCRV_; double maxdt_ = 0.0, btol_ = 0.0, minv_ = 0.0; - SurfaceIdCollection ssids_; - KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit int extrapdebug_ = 0; double tcrvthick_ = 150.0; // CRV sector thickness: should come from geometry service TODO std::unique_ptr extrap_; // calorimeter and other extrapolations @@ -222,18 +220,16 @@ namespace mu2e { } // store surface IDs to be resolved when geometry is available for(auto const& sidname : settings().modSettings().sampleSurfaces()) { - ssids__.push_back(SurfaceId(sidname,-1)); // match all elements + ssids_.push_back(SurfaceId(sidname,-1)); // match all elements } // configure extrapolation if(settings().Extrapolation()){ extrapolate_ = true; // create KKExtrap for calorimeter and upstream extrapolations extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); - toCRV_ = settings().Extrapolation()->ToCRV(); // global configs maxdt_ = settings().Extrapolation()->MaxDt(); btol_ = settings().extSettings().btol(); // use the same BField cor. tolerance as in fit extension - minv_ = settings().Extrapolation()->MinV(); extrapdebug_ = settings().Extrapolation()->Debug(); } @@ -257,10 +253,6 @@ namespace mu2e { // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService GeomHandle kkg_h; auto const& kkg = *kkg_h; - kkg.surfaces(ssids_,surfacess_to_sample_); - // now that geometry is available, resolve sample surface IDs to actual surfaces - GeomHandle kkg_h; - auto const& kkg = *kkg_h; kkg.surfaces(ssids_,sample_); } @@ -387,7 +379,7 @@ namespace mu2e { kktrk.extendTraj(extrange); double tbeg = ftraj.range().begin(); - for(auto const& surf : surfacess_to_sample_){ + for(auto const& surf : sample_){ // search for intersections with each surface from the begining double tstart = tbeg; bool goodinter(true); @@ -416,7 +408,7 @@ namespace mu2e { GeomHandle kkg_h; auto const& kkg = *kkg_h; // extrapolate to the extracted CRV. This function should be migrated to KKExtrap TODO - auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,*kkg.TCRV(),extrapdebug_); + auto TCRV = ExtrapolateTCRV(maxdt_,btol_,intertol_,minv_,kkg.TCRV(),extrapdebug_); auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TCRVSID("TCRV"); From 9b47730de4700d59dfafba41fe08071a26762687 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 30 Apr 2026 18:59:09 -0500 Subject: [PATCH 25/31] remvoed old module --- TrkReco/fcl/prolog.fcl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index aef1c10024..d4ab01dfaa 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/TrkReco/fcl/prolog.fcl @@ -8,13 +8,6 @@ BEGIN_PROLOG # # may use a more careful optimization, but this is it for now #------------------------------------------------------------------------------ -InterdetectorTimeBootstrap : { - kalseedModuleLabel : "KKLine" - caloClusterModuleLabel : "CaloClusterMaker" - crvClusterModuleLabel : "SelectReco" - maxDeltaR : 1000. - maxDeltaT : 1000. -} TrkReco: { McUtils : { tool_type : "TrkRecoMcUtils" From 4c3710296281a43221a22fcc3eed53a1f872b559 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 5 May 2026 10:02:25 -0500 Subject: [PATCH 26/31] material effects --- DataProducts/inc/SurfaceId.hh | 4 +- DataProducts/src/SurfaceId.cc | 4 +- Mu2eKinKal/CMakeLists.txt | 1 + Mu2eKinKal/fcl/prolog.fcl | 7 ++ Mu2eKinKal/inc/ExtrapolateCalo.hh | 3 +- Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh | 136 +++++++++++++++++++++ Mu2eKinKal/inc/KKExtrap.hh | 140 ++++++++++++++++++++-- Mu2eKinKal/inc/KKFitSettings.hh | 1 + Mu2eKinKal/inc/KKMaterial.hh | 7 ++ Mu2eKinKal/inc/KKTrack.hh | 16 +++ Mu2eKinKal/src/ExtrapolateCaloMaterial.cc | 56 +++++++++ Mu2eKinKal/src/KKMaterial.cc | 5 +- TrackerConditions/data/MaterialsList.data | 3 + 13 files changed, 366 insertions(+), 17 deletions(-) create mode 100644 Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh create mode 100644 Mu2eKinKal/src/ExtrapolateCaloMaterial.cc diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index e870cd3664..000ff49c28 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -22,8 +22,8 @@ namespace mu2e { ST_Front=100,ST_Back, ST_Inner, ST_Outer, ST_Foils, ST_Wires, // stopping target bounding surfaces and components TCRV=200, // CRV test planes EMC_Disk_0_Outer = 300, EMC_Disk_0_Inner, EMC_Disk_1_Inner, EMC_Disk_1_Outer, - EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back - + EMC_Disk_0_Front, EMC_Disk_1_Front, EMC_Disk_0_Back, EMC_Disk_1_Back, + EMC_FrontPanel = 320 // calorimeter front panel passive materials (foam, carbon) }; static std::string const& typeName(); diff --git a/DataProducts/src/SurfaceId.cc b/DataProducts/src/SurfaceId.cc index 484a74af4f..512d3748b7 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -34,6 +34,7 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Inner, "ST_Inner"), std::make_pair(SurfaceIdEnum::ST_Outer, "ST_Outer"), std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), + std::make_pair(SurfaceIdEnum::ST_Wires, "ST_Wires"), std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), std::make_pair(SurfaceIdEnum::EMC_Disk_0_Outer, "EMC_Disk_0_Outer"), std::make_pair(SurfaceIdEnum::EMC_Disk_0_Inner, "EMC_Disk_0_Inner"), @@ -42,7 +43,8 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::EMC_Disk_0_Front, "EMC_Disk_0_Front"), std::make_pair(SurfaceIdEnum::EMC_Disk_1_Front, "EMC_Disk_1_Front"), std::make_pair(SurfaceIdEnum::EMC_Disk_0_Back, "EMC_Disk_0_Back"), - std::make_pair(SurfaceIdEnum::EMC_Disk_1_Back, "EMC_Disk_1_Back") + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Back, "EMC_Disk_1_Back"), + std::make_pair(SurfaceIdEnum::EMC_FrontPanel, "EMC_FrontPanel") }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/Mu2eKinKal/CMakeLists.txt b/Mu2eKinKal/CMakeLists.txt index 51fbad3f3e..b402f3002f 100644 --- a/Mu2eKinKal/CMakeLists.txt +++ b/Mu2eKinKal/CMakeLists.txt @@ -4,6 +4,7 @@ cet_make_library( src/CADSHU.cc src/Chi2SHU.cc src/DriftANNSHU.cc + src/ExtrapolateCaloMaterial.cc src/KKBField.cc src/KKConstantBField.cc src/KKFitSettings.cc diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 6e87770a5d..13be7b979e 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -12,6 +12,9 @@ Mu2eKinKal : { strawWireMaterialName : "straw-wire" IPAMaterialName : "HDPE" STMaterialName : "Target" + CaloFrontFoamMaterialName : "AluminumHoneycomb" + CaloFrontCarbonMaterialName : "CarbonFiber" + CaloBackPlateMaterialName : "Polyetheretherketone" IonizationEnergyLossMode : 1 # Moyal mean SolidScatteringFraction : 0.999999 GasScatteringFraction : 0.9999999 @@ -384,6 +387,7 @@ Mu2eKinKal : { } LHDRIFTXTRAP : { + Debug : -300 IntersectionTolerance : 0.1 MaxDt : 200.0 # (ns) BCorrTolerance : 1e-4 # momemntum fraction @@ -393,6 +397,7 @@ Mu2eKinKal : { ToOPA : true ToCaloD0 : true ToCaloD1 : true + ToCaloMaterial : true } LHSEEDXTRAP : { @@ -439,6 +444,7 @@ Mu2eKinKal : { } LINEEXTRAPOLATION : { + Debug : 2 MaxDt : 200.0 BCorrTolerance : 1e-5 ToTrackerEnds : true @@ -448,6 +454,7 @@ Mu2eKinKal : { ToOPA : false ToCaloD0 : true ToCaloD1 : true + ToCaloMaterial : true IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } diff --git a/Mu2eKinKal/inc/ExtrapolateCalo.hh b/Mu2eKinKal/inc/ExtrapolateCalo.hh index d9d3aa7bd2..1906c73dd7 100644 --- a/Mu2eKinKal/inc/ExtrapolateCalo.hh +++ b/Mu2eKinKal/inc/ExtrapolateCalo.hh @@ -118,7 +118,8 @@ namespace mu2e { // update the cache inter_ = newinter; ann_ = disks_[idisk]; - //sid_ = SurfaceId(SurfaceIdEnum::calo_Foils,idisk); //FIXME + SurfaceIdEnum::enum_type diskEnum = (idisk == 0) ? SurfaceIdEnum::EMC_Disk_0_Outer : SurfaceIdEnum::EMC_Disk_1_Outer; + sid_ = SurfaceId(diskEnum, idisk); if(debug_ > 0)std::cout << "Good calo disk " << newinter << " sid " << sid_ << std::endl; return false; } diff --git a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh new file mode 100644 index 0000000000..e3a575f662 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh @@ -0,0 +1,136 @@ +// predicate to extrapolate to passive calorimeter front panel materials (foam, carbon) +// Sophie Middleton (2025) +#ifndef Mu2eKinKal_ExtrapolateCaloMaterial_hh +#define Mu2eKinKal_ExtrapolateCaloMaterial_hh +#include "KinKal/Trajectory/ParticleTrajectory.hh" +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" +#include "cetlib_except/exception.h" +#include +namespace mu2e { + using KinKal::TimeDir; + using KinKal::TimeRange; + using KinKal::Intersection; + using KKGeom::Calo; + using KinKal::Annulus; + using AnnPtr = std::shared_ptr; + + class ExtrapolateCaloMaterial { + public: + // Surface types for calo passive materials + enum SurfaceType { FrontPanelFoam = 0, FrontPanelCarbon = 1, Unknown = 2 }; + + ExtrapolateCaloMaterial() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10), + fp_z_(std::numeric_limits::max()), + inter_(), surftype_(Unknown), intersection_found_(false), fpann_(nullptr), debug_(0) {} + + // Constructor that initializes from calorimeter geometry + // depth: which disk (0 or 1) + ExtrapolateCaloMaterial(double maxdt, double dptol, double intertol, Calo const& calo, + int depth, int debug=0); + + // interface for extrapolation + double maxDt() const { return maxDt_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + auto const& frontPanelAnnulus() const { return fpann_; } + auto const& intersection() const { return inter_; } + auto surfaceType() const { return surftype_; } + int debug() const { return debug_; } + + // extrapolation predicate: returns true if track should continue through material region + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + // reset between tracks + void reset() const { inter_ = Intersection(); surftype_ = Unknown; intersection_found_ = false; } + + // thickness of passive materials (should come from geometry TODO) + static constexpr double foamThickness_ = 21.75; // mm, front panel foam + static constexpr double carbonThickness_ = 3.0; // mm, front panel carbon + + private: + double maxDt_; // maximum extrapolation time + double dptol_; // fractional momentum tolerance + double intertol_; // intersection tolerance (mm) + double fp_z_; // z position of front panel + mutable Intersection inter_; // cache of most recent intersection + mutable SurfaceType surftype_; // type of surface that was intersected + mutable bool intersection_found_; // flag to prevent finding same intersection twice + AnnPtr fpann_; // annulus surface for front panel + int debug_; // debug level + }; + + template bool ExtrapolateCaloMaterial::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { + // Calorimeter has only ONE surface (unlike multi-foil predicates like ST). + // Once we've found the intersection, we must stop - don't search again! + if(debug_ == -300) { + std::cout << " [needsExtrapolation] CHECK: intersection_found_ = " << intersection_found_ << std::endl; + } + if(intersection_found_) { + if(debug_ == -300) std::cout << " [needsExtrapolation] Already found intersection, returning false" << std::endl; + return false; // Already found the one intersection - stop + } + + // Check if we need to continue extrapolating to find calorimeter front panel intersection + // Use FULL trajectory range (particles already past calorimeter position after tracker fit) + auto const& trange = fittraj.range(); + auto stime = trange.begin(); + auto etime = trange.end(); + + auto spos = fittraj.position3(stime); + auto epos = fittraj.position3(etime); + + // Get velocity at start + auto vel_unit = fittraj.direction(stime); + double zvel = vel_unit.Z(); + double zvel_signed = zvel * (tdir == TimeDir::forwards ? 1.0 : -1.0); + + if(debug_ == -300) { + std::cout << "\n[needsExtrapolation] Checking front panel crossing (full trajectory):" << std::endl; + std::cout << " Trajectory Z range: [" << spos.Z() << ", " << epos.Z() << "] mm" << std::endl; + std::cout << " Z velocity: " << zvel << " (signed: " << zvel_signed << ")" << std::endl; + std::cout << " Front panel Z: " << fp_z_ << " mm" << std::endl; + std::cout << " Time range: [" << trange.begin() << ", " << trange.end() << "] ns" << std::endl; + } + + // Quick Z range check - is the front panel even in our trajectory volume? + double zmin = std::min(spos.Z(), epos.Z()); + double zmax = std::max(spos.Z(), epos.Z()); + + if(debug_ == -300) { + std::cout << " [needsExtrapolation] Z-range check: zmin=" << zmin << " zmax=" << zmax << " fp_z_=" << fp_z_ << std::endl; + std::cout << " [needsExtrapolation] Panel in range? " << (fp_z_ >= zmin && fp_z_ <= zmax) << std::endl; + } + + if(fp_z_ < zmin || fp_z_ > zmax) { + // Front panel not in trajectory range yet + if(debug_ == -300) { + std::cout << " Front panel Z " << fp_z_ << " outside trajectory Z range [" << zmin << ", " << zmax << "]" << std::endl; + } + return true; // keep going + } + + // Front panel is within trajectory Z range - test for intersection with annulus + if(debug_ == -300) { + std::cout << " Testing intersection with annulus (full trajectory):" << std::endl; + } + Intersection newinter = KinKal::intersect(fittraj, *fpann_, trange, intertol_); + + if(newinter.good()){ + inter_ = newinter; + intersection_found_ = true; // Mark that we've found THE intersection + if(debug_ == -300) { + std::cout << " Good intersection with front panel FOUND" << std::endl; + } + return false; // stop, we found material crossing + } else { + if(debug_ == -300) { + std::cout << " No intersection found" << std::endl; + } + return true; // keep going + } + } +} +#endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index b1b9f5fd8b..dec35f51d0 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -13,10 +13,12 @@ #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateCalo.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/Mu2eKinKal/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" +#include "cetlib_except/exception.h" #include "Offline/GeometryService/inc/GeometryService.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" namespace mu2e { @@ -38,6 +40,7 @@ namespace mu2e { template bool extrapolateIPA(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; + template bool extrapolateCaloMaterial(KKTrack& ktrk,TimeDir tdir) const; template void toCaloD0(KKTrack& ktrk) const; template void toCaloD1(KKTrack& ktrk) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; @@ -51,17 +54,21 @@ namespace mu2e { DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; FruPtr opaptr_; bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; - bool toCaloD0_, toCaloD1_; + bool toCaloD0_, toCaloD1_, toCaloMaterial_; // calo surfaces and predicates DiskPtr calod0frontptr_, calod0backptr_, calod1frontptr_, calod1backptr_; mutable ExtrapolateToZ calod0Front_, calod0Back_, calod1Front_, calod1Back_; mutable ExtrapolateToZ TSDA_, trackerFront_, trackerBack_; // extrapolation predicate based on Z values mutable ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA mutable ExtrapolateST extrapST_; // extrapolation to intersections with the ST + mutable ExtrapolateCaloMaterial extrapCaloMat_; // extrapolation through passive calo materials mutable bool geom_initialized_ = false; void initGeometry() const; // lazy initialization method double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO + // calorimeter front panel passive material thicknesses (should come from geometry service TODO) + double calofmthick_ = 21.75; // calo front panel foam thickness (mm) + double calocarb_thick_ = 3.0; // calo front panel carbon thickness (mm) }; KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig,KKMaterial const& kkmat) : @@ -81,6 +88,7 @@ namespace mu2e { upstream_(extrapconfig.Upstream()), toCaloD0_(extrapconfig.ToCaloD0()), toCaloD1_(extrapconfig.ToCaloD1()), + toCaloMaterial_(extrapconfig.ToCaloMaterial()), geom_initialized_(false) { // geometry initialization is deferred to first use via initGeometry() @@ -88,7 +96,7 @@ namespace mu2e { inline void KKExtrap::initGeometry() const { if(geom_initialized_) return; // already initialized - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; const_cast(tsdaptr_) = kkg.DS().upstreamAbsorberPtr(); @@ -109,6 +117,8 @@ namespace mu2e { const_cast(trackerBack_) = ExtrapolateToZ(maxdt_,btol_,kkg.tracker().back().center().Z(),debug_); const_cast(extrapIPA_) = ExtrapolateIPA(maxdt_,btol_,intertol_,kkg.DS().innerProtonAbsorberPtr(),debug_); const_cast(extrapST_) = ExtrapolateST(maxdt_,btol_,intertol_,kkg.ST(),debug_); + // Initialize calo material predicates (depth 0 = disk 0 front panel) + const_cast(extrapCaloMat_) = ExtrapolateCaloMaterial(maxdt_,btol_,intertol_,kkg.calo(),0,debug_); // calo predicates - use local Z values stored in Calo class const_cast(calod0Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Front_Z(),debug_); @@ -123,7 +133,7 @@ namespace mu2e { template void KKExtrap::extrapolate(KKTrack& ktrk) const { initGeometry(); // lazy initialization: initialize geometry on first use - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); @@ -141,7 +151,13 @@ namespace mu2e { if(exitsIPA){ // if it exits out the back, extrapolate through the target bool exitsST = extrapolateST(ktrk,tdir); - if(exitsST) { // if it exits out the back, extrapolate to the TSDA (DS rear absorber) + if(exitsST) { // if it exits out the back, optionally extrapolate through calorimeter materials + if(toCaloMaterial_) { + if(!extrapolateCaloMaterial(ktrk,tdir)) { + if(debug_ > 0) std::cout << "Calo material extrapolation did not complete (particle exited backwards)" << std::endl; + } + } + // then extrapolate to the TSDA (DS rear absorber) bool hitTSDA = extrapolateTSDA(ktrk,tdir); // if we hit the TSDA we are done. Otherwise if we reflected, go back through the ST if(!hitTSDA){ // reflection upstream of the target: go back through the target @@ -167,7 +183,7 @@ namespace mu2e { template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); @@ -211,7 +227,7 @@ namespace mu2e { template void KKExtrap::toCaloD0(KKTrack& ktrk) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); @@ -250,7 +266,7 @@ namespace mu2e { template void KKExtrap::toCaloD1(KKTrack& ktrk) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); @@ -283,7 +299,7 @@ namespace mu2e { template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; @@ -328,7 +344,7 @@ namespace mu2e { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; // extraplate the fit through the ST. This will add material effects for each foil intersection. It will continue till the @@ -369,7 +385,7 @@ namespace mu2e { template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); if(trackerFront.debug() > 0)std::cout << "extrapolating to Tracker " << std::endl; @@ -386,9 +402,109 @@ namespace mu2e { return false; } + template bool KKExtrap::extrapolateCaloMaterial(KKTrack& ktrk,TimeDir tdir) const { + using KKCALOMATXING = KKShellXing; + using KKCALOMATXINGPTR = std::shared_ptr; + initGeometry(); + + // Extrapolate through passive calorimeter front panel materials (foam, carbon) + // This adds material effects (energy loss, scattering) for material crossing + // Pattern mirrors IPA and ST extrapolation + + if(debug_ == -300) { + std::cout << "\n### KKExtrap::extrapolateCaloMaterial - START ###" << std::endl; + std::cout << " Direction: " << (tdir == TimeDir::forwards ? "FORWARD" : "BACKWARD") << std::endl; + std::cout << " About to extrapolate through calo front panel passive materials" << std::endl; + std::cout << " Material: AluminumHoneycomb (21.75mm) + CarbonFiber (3.0mm) = 24.75mm total" << std::endl; + } + + auto const& ftraj = ktrk.fitTraj(); + + // Front panel combined thickness: 21.75 mm (foam) + 3.0 mm (carbon) = 24.75 mm + double fp_combined_thick = calofmthick_ + calocarb_thick_; + + // Static surface ID for calo material crossings + static const SurfaceId CaloMatSID(SurfaceIdDetail::EMC_FrontPanel); + + // Check for intersection with front panel (thin surface, only one crossing expected) + int xing_count = 0; + ktrk.extrapolate(tdir,extrapCaloMat_); + if(extrapCaloMat_.intersection().good()){ + xing_count++; + if(debug_ == -300) { + std::cout << "\n === Crossing #" << xing_count << " ===" << std::endl; + std::cout << " Found intersection with calo front panel Annulus" << std::endl; + auto const& inter = extrapCaloMat_.intersection(); + std::cout << " Intersection: " << inter << std::endl; + } + + // We have a good intersection. Create a shell material Xing for the annulus surface + auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + auto const& annulusptr = extrapCaloMat_.frontPanelAnnulus(); + + if(debug_ == -300) { + std::cout << " Before material: momentum = " << reftrajptr->momentum() << " MeV/c" << std::endl; + } + + // Create the material xing using combined front panel thickness + KKCALOMATXINGPTR matxingptr = std::make_shared( + annulusptr, CaloMatSID, *kkmat_.CaloFrontFoamMaterial(), + extrapCaloMat_.intersection(), reftrajptr, fp_combined_thick, extrapCaloMat_.interTolerance()); + + if(debug_ == -300){ + double dmom, paramomvar, perpmomvar; + matxingptr->materialEffects(dmom,paramomvar,perpmomvar); + std::cout << " Material Effects:" << std::endl; + std::cout << " dE/dx momentum loss: " << dmom << " MeV" << std::endl; + std::cout << " Parallel scattering sigma: " << std::sqrt(paramomvar) << " MeV" << std::endl; + std::cout << " Perpendicular scattering sigma: " << std::sqrt(perpmomvar) << " MeV" << std::endl; + } + + // Try to add the xing. Note: This may fail if the surface is in the trajectory past + // In that case, we just track the intersection for analysis + try { + ktrk.addCaloMatXing(matxingptr,tdir); + if(debug_ == -300){ + auto const& newtrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + std::cout << " After material: momentum = " << newtrajptr->momentum() << " MeV/c" << std::endl; + } + // Clear the cached intersection after successfully adding it + extrapCaloMat_.reset(); + } catch (cet::exception const& e) { + if(debug_ > 0) { + std::cout << "Warning: Could not add calo material xing: " << e.what() << std::endl; + std::cout << "Intersection will be tracked but not applied to fit" << std::endl; + } + // Still clear the cached intersection to prevent re-finding it + extrapCaloMat_.reset(); + } + } + + if(debug_ == -300) { + std::cout << "\n === Loop Complete ===" << std::endl; + std::cout << " Total crossings found: " << xing_count << std::endl; + } + + // check if the particle exited in the same physical direction or not (reflection) + double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto startdir = ftraj.direction(starttime); + double endtime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); + auto enddir = ftraj.direction(endtime); + + bool no_reflection = (enddir.Z() * startdir.Z() > 0.0); + + if(debug_ == -300) { + std::cout << " Reflection check: startdir.Z()=" << startdir.Z() << " enddir.Z()=" << enddir.Z() << std::endl; + std::cout << " No reflection (particle exited in same Z direction): " << no_reflection << std::endl; + std::cout << "### KKExtrap::extrapolateCaloMaterial - END ###\n" << std::endl; + } + + return no_reflection; + } + template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; ExtrapolateToZ TSDA(maxdt_,btol_,kkg.DS().upstreamAbsorber().center().Z(),debug_); if(TSDA.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; @@ -409,7 +525,7 @@ namespace mu2e { template void KKExtrap::toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const { initGeometry(); - GeomHandle kkg_h; + GeomHandle kkg_h; auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId OPASID("OPA"); diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 36710ed780..ac7c0cb015 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -116,6 +116,7 @@ namespace mu2e { fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom ToCaloD0 { Name("ToCaloD0"), Comment("Extrapolate tracks to the disk0 ends") }; fhicl::Atom ToCaloD1 { Name("ToCaloD1"), Comment("Extrapolate tracks to the disk1 ends") }; + fhicl::Atom ToCaloMaterial { Name("ToCaloMaterial"), Comment("Extrapolate through calorimeter front panel passive materials"), false }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; }; diff --git a/Mu2eKinKal/inc/KKMaterial.hh b/Mu2eKinKal/inc/KKMaterial.hh index 127d59642d..89588c7d0e 100644 --- a/Mu2eKinKal/inc/KKMaterial.hh +++ b/Mu2eKinKal/inc/KKMaterial.hh @@ -30,6 +30,9 @@ namespace mu2e { fhicl::Atom strawWireMaterialName{ Name("strawWireMaterialName"), Comment("strawWireMaterialName") }; fhicl::Atom IPAMaterialName{ Name("IPAMaterialName"), Comment("IPA MaterialName") }; fhicl::Atom STMaterialName{ Name("STMaterialName"), Comment("Stopping Target MaterialName") }; + fhicl::Atom CaloFrontFoamMaterialName{ Name("CaloFrontFoamMaterialName"), Comment("Calo front panel foam MaterialName") }; + fhicl::Atom CaloFrontCarbonMaterialName{ Name("CaloFrontCarbonMaterialName"), Comment("Calo front panel carbon MaterialName") }; + fhicl::Atom CaloBackPlateMaterialName{ Name("CaloBackPlateMaterialName"), Comment("Calo back plate MaterialName") }; fhicl::Atom elossMode { Name("IonizationEnergyLossMode"), Comment( "Ionization energy loss mode") }; fhicl::Atom solidScatter{ Name("SolidScatteringFraction"), Comment("DahlLynch Scattering model cutoff Fraction for solids") }; fhicl::Atom gasScatter{ Name("GasScatteringFraction"), Comment("DahlLynch Scattering model cutoff Fraction for gases") }; @@ -40,9 +43,13 @@ namespace mu2e { KKStrawMaterial const& strawMaterial() const; auto IPAMaterial() const { return matdbinfo_->findDetMaterial(ipamatname_); } auto STMaterial() const { return matdbinfo_->findDetMaterial(stmatname_); } + auto CaloFrontFoamMaterial() const { return matdbinfo_->findDetMaterial(calofoammatname_); } + auto CaloFrontCarbonMaterial() const { return matdbinfo_->findDetMaterial(calocarbonmatname_); } + auto CaloBackPlateMaterial() const { return matdbinfo_->findDetMaterial(calobackplatematname_); } private: KKFileFinder filefinder_; // used to find material info std::string wallmatname_, gasmatname_, wirematname_,ipamatname_, stmatname_; + std::string calofoammatname_, calocarbonmatname_, calobackplatematname_; mutable std::unique_ptr matdbinfo_; // material database mutable std::unique_ptr smat_; // straw material }; diff --git a/Mu2eKinKal/inc/KKTrack.hh b/Mu2eKinKal/inc/KKTrack.hh index c6dc11848c..8ee2f3f67c 100644 --- a/Mu2eKinKal/inc/KKTrack.hh +++ b/Mu2eKinKal/inc/KKTrack.hh @@ -48,6 +48,9 @@ namespace mu2e { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; using KKSTXINGCOL = std::vector; + using KKCALOMAT_XING = KKShellXing; + using KKCALOMAT_XINGPTR = std::shared_ptr; + using KKCALOMAT_XINGCOL = std::vector; using KKCRVXING = KKShellXing; using KKCRVXINGPTR = std::shared_ptr; using KKCRVXINGCOL = std::vector; @@ -125,6 +128,8 @@ namespace mu2e { void addIPAXing(KKIPAXINGPTR const& ipaxing,TimeDir const& tdir); // add ST Xing void addSTXing(KKSTXINGPTR const& stxing,TimeDir const& tdir); + // add calo material Xing + void addCaloMatXing(KKCALOMAT_XINGPTR const& calomatxing,TimeDir const& tdir); // add TCRV Xing void addTCRVXing(KKCRVXINGPTR const& crvxing,TimeDir const& tdir); // add intersections @@ -137,6 +142,7 @@ namespace mu2e { KKSTRAWXINGCOL const& strawXings() const { return strawxings_; } KKIPAXINGCOL const& IPAXings() const { return ipaxings_; } KKSTXINGCOL const& STXings() const { return stxings_; } + KKCALOMAT_XINGCOL const& CaloMatXings() const { return calomatxings_; } KKCRVXINGCOL const& CRVXings() const { return crvxings_; } KKINTERCOL const& intersections() const { return inters_; } KKCALOHITCOL const& caloHits() const { return calohits_; } @@ -155,6 +161,7 @@ namespace mu2e { KKCALOHITCOL calohits_; // calo hits used in this fit KKIPAXINGCOL ipaxings_; // ipa material crossings used in extrapolation KKSTXINGCOL stxings_; // stopping target material crossings used in extrapolation + KKCALOMAT_XINGCOL calomatxings_; // calorimeter passive material crossings used in extrapolation KKCRVXINGCOL crvxings_; // crv crossings using in extrapolation KKINTERCOL inters_; // other recorded intersections KKSTRAWHITCLUSTERCOL strawhitclusters_; // straw hit clusters used in this fit @@ -339,6 +346,15 @@ namespace mu2e { stxings_.push_back(stxingptr); } + template void KKTrack::addCaloMatXing(KKCALOMAT_XINGPTR const& calomatxingptr,TimeDir const& tdir) { + // convert to a generic Xing + std::shared_ptr> exptr = std::static_pointer_cast>(calomatxingptr); + // extrapolate the fit through this xing + if(!this->extrapolate(exptr,tdir))throw cet::exception("RECO")<<"mu2e::KKTrack: Calo Material extrapolation failure"<< std::endl; + // store the xing + calomatxings_.push_back(calomatxingptr); + } + template void KKTrack::addTCRVXing(KKCRVXINGPTR const& crvxingptr,TimeDir const& tdir) { // convert to a generic Xing std::shared_ptr> exptr = std::static_pointer_cast>(crvxingptr); diff --git a/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc new file mode 100644 index 0000000000..0a188c123a --- /dev/null +++ b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc @@ -0,0 +1,56 @@ +#include "Offline/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" +#include "KinKal/Geometry/Annulus.hh" +#include +#include + +namespace mu2e { + using KinKal::VEC3; + using KKGeom::Calo; + + ExtrapolateCaloMaterial::ExtrapolateCaloMaterial(double maxdt, double dptol, double intertol, + Calo const& calo, int depth, int debug) : + maxDt_(maxdt), dptol_(dptol), intertol_(intertol), fp_z_(std::numeric_limits::max()), + inter_(), surftype_(Unknown), intersection_found_(false), fpann_(nullptr), debug_(debug) + { + // Get disk 0 or 1 based on depth + auto const& diskref = (depth == 0) ? calo.EMC_Disk_0_Front() : calo.EMC_Disk_1_Front(); + + // Z position of front panel + // Front panel dimensions from geometry: + // FPInnerRadius = 336, FPOuterRadius = 680 + // Total FP thickness: foam (21.75) + carbon (3.0) = 24.75 mm + double fp_inner_r = 336.0; + double fp_outer_r = 680.0; + double disk_z = diskref.center().Z(); + + // Front panel surfaces: positioned ~50mm before disk (FOA position) + // This should be read from geometry service TODO + fp_z_ = disk_z - 50.0; // approximate position + + // Create annulus surface for front panel (combined foam + carbon layers) + // Annulus constructor: Annulus(norm, udir, center, innerrad, outerrad) + VEC3 norm(0, 0, 1); // normal to disk (Z-direction) + VEC3 udir(1, 0, 0); // u-direction along disk (radial) + VEC3 center_fp(0, 0, fp_z_); + fpann_ = std::make_shared(norm, udir, center_fp, + fp_inner_r, fp_outer_r); + + if(debug_ == -300) { + std::cout << "\n=== ExtrapolateCaloMaterial Initialization ==="<< std::endl; + std::cout << " Disk index: " << depth << std::endl; + std::cout << " Disk center Z: " << disk_z << " mm" << std::endl; + std::cout << " Front panel Z: " << fp_z_ << " mm (offset: " << (disk_z - fp_z_) << " mm upstream)" << std::endl; + std::cout << " Annulus geometry:" << std::endl; + std::cout << " Inner radius: " << fp_inner_r << " mm" << std::endl; + std::cout << " Outer radius: " << fp_outer_r << " mm" << std::endl; + std::cout << " Normal: (" << norm.X() << ", " << norm.Y() << ", " << norm.Z() << ")" << std::endl; + std::cout << " U-direction: (" << udir.X() << ", " << udir.Y() << ", " << udir.Z() << ")" << std::endl; + std::cout << " Material thickness: 21.75 mm foam + 3.0 mm carbon = 24.75 mm total" << std::endl; + std::cout << " Max extrapolation time: " << maxDt_ << " (unbounded if < 0)" << std::endl; + std::cout << " Intersection tolerance: " << intertol_ << " mm" << std::endl; + std::cout << "=== Initialization Complete ===\n" << std::endl; + } + } + +} diff --git a/Mu2eKinKal/src/KKMaterial.cc b/Mu2eKinKal/src/KKMaterial.cc index 126193fc9b..a69d90e289 100644 --- a/Mu2eKinKal/src/KKMaterial.cc +++ b/Mu2eKinKal/src/KKMaterial.cc @@ -14,7 +14,10 @@ namespace mu2e { gasmatname_(matconfig.strawGasMaterialName()), wirematname_(matconfig.strawWireMaterialName()), ipamatname_(matconfig.IPAMaterialName()), - stmatname_(matconfig.STMaterialName()) { + stmatname_(matconfig.STMaterialName()), + calofoammatname_(matconfig.CaloFrontFoamMaterialName()), + calocarbonmatname_(matconfig.CaloFrontCarbonMaterialName()), + calobackplatematname_(matconfig.CaloBackPlateMaterialName()) { MatEnv::DetMaterialConfig dmconf; dmconf.elossmode_ = (DetMaterial::energylossmode)matconfig.elossMode(); dmconf.scatterfrac_solid_ = matconfig.solidScatter(); diff --git a/TrackerConditions/data/MaterialsList.data b/TrackerConditions/data/MaterialsList.data index 9d2b64e7e6..bad5f37027 100644 --- a/TrackerConditions/data/MaterialsList.data +++ b/TrackerConditions/data/MaterialsList.data @@ -23,6 +23,9 @@ Mylar 1.4 0. 0. -3 10 Carbon 0 4 Oxygen 0 8 Hydrogen 0 straw-wire 19.300000 0. 0. +1 100.0e-2 Tungsten 0 -10 -20 -30 20.0 1.0 solid # 80:20 (by volume) Argon:CO2 straw-gas 0.0016980 0. 0. +2 0.78304 Argon 0 0.21696 CO2 1 -10 -20 -30 20.0 1.0 gas 0.01 +# Calorimeter materials +AluminumHoneycomb 0.03 0. 0. +1 100.0e-2 Aluminum 0 -10 -20 -30 20.0 1.0 solid +CarbonFiber 1.60 0. 0. +1 100.0e-2 Carbon 0 -10 -20 -30 20.0 1.0 solid straw-wall 1.4325 0. 0. +3 96.95e-2 Mylar 1 1.80e-2 Gold 0 1.25e-2 Aluminum 0 -10 -20 -30 20.0 1.0 solid 1.5 #straw-wall 1.48 0. 0. +3 97.0e-2 Mylar 1 2.0e-2 Copper 0 1.0e-2 Aluminum 0 -10 -20 -30 20.0 1.0 solid cathode 5.2 0. 0. +2 86.4e-2 Copper 0 13.6e-2 straw-wall 1 -10 -20 -30 20.0 1.0 solid From e8a2b455fc50d619d4b6b90e7011949aaa8a0a71 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 5 May 2026 10:33:23 -0500 Subject: [PATCH 27/31] works --- Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh index e3a575f662..143ae266e9 100644 --- a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh +++ b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh @@ -127,9 +127,9 @@ namespace mu2e { return false; // stop, we found material crossing } else { if(debug_ == -300) { - std::cout << " No intersection found" << std::endl; + std::cout << " No intersection found with annulus in Z-range" << std::endl; } - return true; // keep going + return false; // stop - panel in range but no intersection means we're past it } } } From 52ab0e6fef76d21b96139bfe79761ab38ba37500 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Tue, 5 May 2026 10:51:12 -0500 Subject: [PATCH 28/31] works --- Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh index 143ae266e9..299c579002 100644 --- a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh +++ b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh @@ -1,5 +1,5 @@ // predicate to extrapolate to passive calorimeter front panel materials (foam, carbon) -// Sophie Middleton (2025) +// Sophie Middleton (2026) #ifndef Mu2eKinKal_ExtrapolateCaloMaterial_hh #define Mu2eKinKal_ExtrapolateCaloMaterial_hh #include "KinKal/Trajectory/ParticleTrajectory.hh" From 3df021bab933afe6ed60fd7d20bdf259ad879ed1 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Wed, 6 May 2026 15:30:15 -0500 Subject: [PATCH 29/31] removed --- Mu2eKinKal/fcl/prolog.fcl | 5 +- Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh | 39 ++++++- Mu2eKinKal/inc/KKExtrap.hh | 126 +++++++++++++++++++--- Mu2eKinKal/src/ExtrapolateCaloMaterial.cc | 54 +++++++++- Mu2eKinKal/src/LoopHelixFit_module.cc | 85 +++++++++++++++ 5 files changed, 288 insertions(+), 21 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 13be7b979e..2b3e63b35f 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -397,7 +397,7 @@ Mu2eKinKal : { ToOPA : true ToCaloD0 : true ToCaloD1 : true - ToCaloMaterial : true + ToCaloMaterial : false } LHSEEDXTRAP : { @@ -410,6 +410,7 @@ Mu2eKinKal : { ToOPA : false ToCaloD0 : true ToCaloD1 : true + ToCaloMaterial : false } CENTRALHELIX : { @@ -454,7 +455,7 @@ Mu2eKinKal : { ToOPA : false ToCaloD0 : true ToCaloD1 : true - ToCaloMaterial : true + ToCaloMaterial : false IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } diff --git a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh index 299c579002..d84e583c0c 100644 --- a/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh +++ b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh @@ -10,6 +10,7 @@ #include "Offline/KinKalGeom/inc/Calo.hh" #include "cetlib_except/exception.h" #include +#include "TH2F.h" namespace mu2e { using KinKal::TimeDir; using KinKal::TimeRange; @@ -24,7 +25,7 @@ namespace mu2e { enum SurfaceType { FrontPanelFoam = 0, FrontPanelCarbon = 1, Unknown = 2 }; ExtrapolateCaloMaterial() : maxDt_(-1.0), dptol_(1e10), intertol_(1e10), - fp_z_(std::numeric_limits::max()), + fp_z_(std::numeric_limits::max()), disk_(0), inter_(), surftype_(Unknown), intersection_found_(false), fpann_(nullptr), debug_(0) {} // Constructor that initializes from calorimeter geometry @@ -46,6 +47,16 @@ namespace mu2e { // reset between tracks void reset() const { inter_ = Intersection(); surftype_ = Unknown; intersection_found_ = false; } + // Validation plot methods - histograms created externally via TFileService + void setValidationHistograms(TH2F* h_eff, TH2F* h_pos, TH1F* h_deg = nullptr, TH1F* h_scat100 = nullptr, TH1F* h_scat80 = nullptr) { + h_intersection_efficiency_ = h_eff; + h_frontpanel_hits_ = h_pos; + h_momentum_degradation_ = h_deg; + h_scatter_100mev_ = h_scat100; + h_scatter_80mev_ = h_scat80; + } + void fillValidationPlots(double momentum, int disk, double pos_x, double pos_y, double dmom = 0.0, double scatter_angle = 0.0) const; + // thickness of passive materials (should come from geometry TODO) static constexpr double foamThickness_ = 21.75; // mm, front panel foam static constexpr double carbonThickness_ = 3.0; // mm, front panel carbon @@ -55,11 +66,19 @@ namespace mu2e { double dptol_; // fractional momentum tolerance double intertol_; // intersection tolerance (mm) double fp_z_; // z position of front panel + int disk_; // disk ID (0 or 1) mutable Intersection inter_; // cache of most recent intersection mutable SurfaceType surftype_; // type of surface that was intersected mutable bool intersection_found_; // flag to prevent finding same intersection twice AnnPtr fpann_; // annulus surface for front panel int debug_; // debug level + + // Validation histograms (managed externally by TFileService) + mutable TH2F* h_intersection_efficiency_ = nullptr; // Plot 1: efficiency vs momentum and disk + mutable TH2F* h_frontpanel_hits_ = nullptr; // Plot 2: hit distribution on front panel (phi vs r) + mutable TH1F* h_momentum_degradation_ = nullptr; // Plot 3: momentum degradation distribution + mutable TH1F* h_scatter_100mev_ = nullptr; // Plot 4: scattering angle at 100 MeV + mutable TH1F* h_scatter_80mev_ = nullptr; // Plot 4: scattering angle at 80 MeV }; template bool ExtrapolateCaloMaterial::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { @@ -118,9 +137,27 @@ namespace mu2e { } Intersection newinter = KinKal::intersect(fittraj, *fpann_, trange, intertol_); + if(debug_ > 0 || debug_ == -300) { + std::cout << "[ExtrapolateCaloMaterial::needsExtrapolation] Intersection test result:" << std::endl; + std::cout << " good() = " << newinter.good() << std::endl; + if(newinter.good()) { + std::cout << " *** INTERSECTION FOUND ***" << std::endl; + std::cout << " Position: (" << newinter.pos_.X() << ", " \ + << newinter.pos_.Y() << ", " << newinter.pos_.Z() << ")" << std::endl; + } else { + std::cout << " No intersection detected" << std::endl; + } + } + if(newinter.good()){ inter_ = newinter; intersection_found_ = true; // Mark that we've found THE intersection + + // Fill validation plots + double momentum = fittraj.momentum(newinter.time_); + auto pos = newinter.pos_; + fillValidationPlots(momentum, disk_, pos.X(), pos.Y()); + if(debug_ == -300) { std::cout << " Good intersection with front panel FOUND" << std::endl; } diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index dec35f51d0..1704356b17 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -21,6 +21,7 @@ #include "cetlib_except/exception.h" #include "Offline/GeometryService/inc/GeometryService.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "TProfile.h" namespace mu2e { using KinKal::VEC3; using KinKal::TimeDir; @@ -45,6 +46,8 @@ namespace mu2e { template void toCaloD1(KKTrack& ktrk) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; + // validation histogram support + void setValidationHistograms(class TH2F* h_eff, class TH2F* h_pos, class TH1F* h_deg = nullptr, class TH1F* h_scat100 = nullptr, class TH1F* h_scat80 = nullptr, class TProfile* h_res = nullptr) const; private: int debug_; @@ -61,7 +64,8 @@ namespace mu2e { mutable ExtrapolateToZ TSDA_, trackerFront_, trackerBack_; // extrapolation predicate based on Z values mutable ExtrapolateIPA extrapIPA_; // extrapolation to intersections with the IPA mutable ExtrapolateST extrapST_; // extrapolation to intersections with the ST - mutable ExtrapolateCaloMaterial extrapCaloMat_; // extrapolation through passive calo materials + mutable ExtrapolateCaloMaterial extrapCaloMat_; // extrapolation through passive calo materials (disk 0) + mutable ExtrapolateCaloMaterial extrapCaloMatD1_; // extrapolation through passive calo materials (disk 1) mutable bool geom_initialized_ = false; void initGeometry() const; // lazy initialization method double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO @@ -69,6 +73,13 @@ namespace mu2e { // calorimeter front panel passive material thicknesses (should come from geometry service TODO) double calofmthick_ = 21.75; // calo front panel foam thickness (mm) double calocarb_thick_ = 3.0; // calo front panel carbon thickness (mm) + // Validation histogram pointers - stored to re-apply after geometry initialization + mutable class TH2F* h_intersection_efficiency_ptr_ = nullptr; + mutable class TH2F* h_frontpanel_hits_ptr_ = nullptr; + mutable class TH1F* h_momentum_degradation_ptr_ = nullptr; + mutable class TH1F* h_scatter_100mev_ptr_ = nullptr; + mutable class TH1F* h_scatter_80mev_ptr_ = nullptr; + mutable class TProfile* h_resolution_vs_momentum_ptr_ = nullptr; }; KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig,KKMaterial const& kkmat) : @@ -117,8 +128,18 @@ namespace mu2e { const_cast(trackerBack_) = ExtrapolateToZ(maxdt_,btol_,kkg.tracker().back().center().Z(),debug_); const_cast(extrapIPA_) = ExtrapolateIPA(maxdt_,btol_,intertol_,kkg.DS().innerProtonAbsorberPtr(),debug_); const_cast(extrapST_) = ExtrapolateST(maxdt_,btol_,intertol_,kkg.ST(),debug_); - // Initialize calo material predicates (depth 0 = disk 0 front panel) + // Initialize calo material predicates for both disks const_cast(extrapCaloMat_) = ExtrapolateCaloMaterial(maxdt_,btol_,intertol_,kkg.calo(),0,debug_); + const_cast(extrapCaloMatD1_) = ExtrapolateCaloMaterial(maxdt_,btol_,intertol_,kkg.calo(),1,debug_); + + // Re-apply validation histograms if they were previously set + if(h_intersection_efficiency_ptr_ || h_frontpanel_hits_ptr_ || h_momentum_degradation_ptr_ || h_scatter_100mev_ptr_ || h_scatter_80mev_ptr_ || h_resolution_vs_momentum_ptr_) { + const_cast(extrapCaloMat_).setValidationHistograms(h_intersection_efficiency_ptr_, h_frontpanel_hits_ptr_, h_momentum_degradation_ptr_, h_scatter_100mev_ptr_, h_scatter_80mev_ptr_); + const_cast(extrapCaloMatD1_).setValidationHistograms(h_intersection_efficiency_ptr_, h_frontpanel_hits_ptr_, h_momentum_degradation_ptr_, h_scatter_100mev_ptr_, h_scatter_80mev_ptr_); + if(debug_ > 0) { + std::cout << "Re-applied validation histograms to both calo material extrapolators after geometry initialization" << std::endl; + } + } // calo predicates - use local Z values stored in Calo class const_cast(calod0Front_) = ExtrapolateToZ(maxdt_,btol_,kkg.calo().EMC_Disk_0_Front_Z(),debug_); @@ -179,6 +200,20 @@ namespace mu2e { // optionally test for intersection with the OPA if(toOPA_)toOPA(ktrk,starttime,tdir); } + + // Fill momentum resolution histogram after all extrapolations + if(h_resolution_vs_momentum_ptr_) { + auto const& ftraj = ktrk.fitTraj(); + double t_mid = ftraj.range().mid(); + double p_fit = ftraj.momentum(t_mid); + if(p_fit > 0.0) { + // Get momentum variance from the state estimate + auto state = ftraj.stateEstimate(t_mid); + double p_err = std::sqrt(state.momentumVariance()); // radial momentum uncertainty + double rel_err = p_err / p_fit; + h_resolution_vs_momentum_ptr_->Fill(p_fit, rel_err); + } + } } template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { @@ -407,14 +442,14 @@ namespace mu2e { using KKCALOMATXINGPTR = std::shared_ptr; initGeometry(); - // Extrapolate through passive calorimeter front panel materials (foam, carbon) + // Extrapolate through passive calorimeter front panel materials (foam, carbon) for BOTH disks // This adds material effects (energy loss, scattering) for material crossing // Pattern mirrors IPA and ST extrapolation if(debug_ == -300) { std::cout << "\n### KKExtrap::extrapolateCaloMaterial - START ###" << std::endl; std::cout << " Direction: " << (tdir == TimeDir::forwards ? "FORWARD" : "BACKWARD") << std::endl; - std::cout << " About to extrapolate through calo front panel passive materials" << std::endl; + std::cout << " About to extrapolate through BOTH calo front panel disks" << std::endl; std::cout << " Material: AluminumHoneycomb (21.75mm) + CarbonFiber (3.0mm) = 24.75mm total" << std::endl; } @@ -426,21 +461,39 @@ namespace mu2e { // Static surface ID for calo material crossings static const SurfaceId CaloMatSID(SurfaceIdDetail::EMC_FrontPanel); - // Check for intersection with front panel (thin surface, only one crossing expected) - int xing_count = 0; - ktrk.extrapolate(tdir,extrapCaloMat_); - if(extrapCaloMat_.intersection().good()){ - xing_count++; + // Lambda function to process a single disk + auto processDiskMaterial = [&](ExtrapolateCaloMaterial& extrapCaloMatDisk, int disk_id) { + // Check for intersection with this disk's front panel + if(!ktrk.extrapolate(tdir, extrapCaloMatDisk)) { + if(debug_ == -300) { + std::cout << " [Disk " << disk_id << "] Did not reach front panel" << std::endl; + } + return false; // Particle didn't reach this disk + } + + if(!extrapCaloMatDisk.intersection().good()) { + if(debug_ == -300) { + std::cout << " [Disk " << disk_id << "] No intersection with annulus" << std::endl; + } + return false; // No intersection + } + if(debug_ == -300) { - std::cout << "\n === Crossing #" << xing_count << " ===" << std::endl; + std::cout << "\n === Disk " << disk_id << " Crossing ===" << std::endl; std::cout << " Found intersection with calo front panel Annulus" << std::endl; - auto const& inter = extrapCaloMat_.intersection(); + auto const& inter = extrapCaloMatDisk.intersection(); std::cout << " Intersection: " << inter << std::endl; } // We have a good intersection. Create a shell material Xing for the annulus surface auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - auto const& annulusptr = extrapCaloMat_.frontPanelAnnulus(); + auto const& annulusptr = extrapCaloMatDisk.frontPanelAnnulus(); + + // Fill validation histograms + auto const& inter = extrapCaloMatDisk.intersection(); + double momentum = reftrajptr->momentum(); + double xpos = inter.pos_.X(); + double ypos = inter.pos_.Y(); if(debug_ == -300) { std::cout << " Before material: momentum = " << reftrajptr->momentum() << " MeV/c" << std::endl; @@ -449,11 +502,24 @@ namespace mu2e { // Create the material xing using combined front panel thickness KKCALOMATXINGPTR matxingptr = std::make_shared( annulusptr, CaloMatSID, *kkmat_.CaloFrontFoamMaterial(), - extrapCaloMat_.intersection(), reftrajptr, fp_combined_thick, extrapCaloMat_.interTolerance()); + extrapCaloMatDisk.intersection(), reftrajptr, fp_combined_thick, extrapCaloMatDisk.interTolerance()); + + // Extract momentum degradation from material effects + double dmom = 0.0, paramomvar = 0.0, perpmomvar = 0.0; + matxingptr->materialEffects(dmom, paramomvar, perpmomvar); + double dmom_abs = std::abs(dmom); + + // Calculate scattering angle from perpendicular momentum uncertainty + double scatter_angle = 0.0; + if (perpmomvar > 0.0 && momentum > 0.0) { + double p_perp_err = std::sqrt(perpmomvar); + scatter_angle = p_perp_err / momentum; // scattering angle in radians + } + + // Fill validation histograms with momentum degradation and scattering angle + extrapCaloMatDisk.fillValidationPlots(momentum, disk_id, xpos, ypos, dmom_abs, scatter_angle); if(debug_ == -300){ - double dmom, paramomvar, perpmomvar; - matxingptr->materialEffects(dmom,paramomvar,perpmomvar); std::cout << " Material Effects:" << std::endl; std::cout << " dE/dx momentum loss: " << dmom << " MeV" << std::endl; std::cout << " Parallel scattering sigma: " << std::sqrt(paramomvar) << " MeV" << std::endl; @@ -469,15 +535,26 @@ namespace mu2e { std::cout << " After material: momentum = " << newtrajptr->momentum() << " MeV/c" << std::endl; } // Clear the cached intersection after successfully adding it - extrapCaloMat_.reset(); + extrapCaloMatDisk.reset(); + return true; } catch (cet::exception const& e) { if(debug_ > 0) { std::cout << "Warning: Could not add calo material xing: " << e.what() << std::endl; std::cout << "Intersection will be tracked but not applied to fit" << std::endl; } // Still clear the cached intersection to prevent re-finding it - extrapCaloMat_.reset(); + extrapCaloMatDisk.reset(); + return false; } + }; + + // Process both disk 0 and disk 1 + int xing_count = 0; + if(processDiskMaterial(const_cast(extrapCaloMat_), 0)) { + xing_count++; + } + if(processDiskMaterial(const_cast(extrapCaloMatD1_), 1)) { + xing_count++; } if(debug_ == -300) { @@ -535,5 +612,20 @@ namespace mu2e { ktrk.addIntersection(OPASID,opainter); } } + + inline void KKExtrap::setValidationHistograms(class TH2F* h_eff, class TH2F* h_pos, class TH1F* h_deg, class TH1F* h_scat100, class TH1F* h_scat80, class TProfile* h_res) const { + // Store pointers for later re-application after geometry initialization + const_cast(this)->h_intersection_efficiency_ptr_ = h_eff; + const_cast(this)->h_frontpanel_hits_ptr_ = h_pos; + const_cast(this)->h_momentum_degradation_ptr_ = h_deg; + const_cast(this)->h_scatter_100mev_ptr_ = h_scat100; + const_cast(this)->h_scatter_80mev_ptr_ = h_scat80; + const_cast(this)->h_resolution_vs_momentum_ptr_ = h_res; + // If geometry is already initialized, apply to extrapCaloMat_ immediately + if(geom_initialized_) { + const_cast(extrapCaloMat_).setValidationHistograms(h_eff, h_pos, h_deg, h_scat100, h_scat80); + const_cast(extrapCaloMatD1_).setValidationHistograms(h_eff, h_pos, h_deg, h_scat100, h_scat80); + } + } } #endif diff --git a/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc index 0a188c123a..98c225298a 100644 --- a/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc +++ b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc @@ -1,8 +1,11 @@ +// ExtrapolateCaloMaterial.cc +// Sophie Middleton (Caltech), 2026 #include "Offline/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh" #include "Offline/KinKalGeom/inc/Calo.hh" #include "KinKal/Geometry/Annulus.hh" #include #include +#include namespace mu2e { using KinKal::VEC3; @@ -11,7 +14,7 @@ namespace mu2e { ExtrapolateCaloMaterial::ExtrapolateCaloMaterial(double maxdt, double dptol, double intertol, Calo const& calo, int depth, int debug) : maxDt_(maxdt), dptol_(dptol), intertol_(intertol), fp_z_(std::numeric_limits::max()), - inter_(), surftype_(Unknown), intersection_found_(false), fpann_(nullptr), debug_(debug) + disk_(depth), inter_(), surftype_(Unknown), intersection_found_(false), fpann_(nullptr), debug_(debug) { // Get disk 0 or 1 based on depth auto const& diskref = (depth == 0) ? calo.EMC_Disk_0_Front() : calo.EMC_Disk_1_Front(); @@ -53,4 +56,53 @@ namespace mu2e { } } + // Fill validation plots with intersection data + void ExtrapolateCaloMaterial::fillValidationPlots(double momentum, int disk, double pos_x, double pos_y, double dmom, double scatter_angle) const { + // Histograms are created externally via TFileService + // Only fill if they have been provided + if (!h_intersection_efficiency_ && !h_frontpanel_hits_ && !h_momentum_degradation_ && !h_scatter_100mev_ && !h_scatter_80mev_) { + if(debug_ > 0) { + std::cout << "[fillValidationPlots] WARNING: No histograms available!" << std::endl; + } + return; + } + + // Plot 1: Intersection Efficiency Map + if (h_intersection_efficiency_) { + h_intersection_efficiency_->Fill(momentum, disk); + } + + // Plot 2: Front Panel Hit Positions + // Convert (x, y) to cylindrical coordinates (r, phi) + if (h_frontpanel_hits_) { + double r = std::sqrt(pos_x * pos_x + pos_y * pos_y); + double phi = std::atan2(pos_y, pos_x); + if (phi < 0) phi += 2.0 * M_PI; + h_frontpanel_hits_->Fill(phi, r); + } + + // Plot 3: Momentum Degradation Distribution + if (h_momentum_degradation_ && dmom > 0.0) { + h_momentum_degradation_->Fill(dmom); + } + + // Plot 4: Multiple Scattering Angle Distribution + if (scatter_angle > 0.0) { + double scatter_mrad = scatter_angle * 1000.0; // convert to milliradians + if (h_scatter_100mev_ && momentum > 95.0 && momentum < 105.0) { + h_scatter_100mev_->Fill(scatter_mrad); + } + if (h_scatter_80mev_ && momentum > 75.0 && momentum < 85.0) { + h_scatter_80mev_->Fill(scatter_mrad); + } + } + + if(debug_ > 0) { + std::cout << "[fillValidationPlots] FILLED: momentum=" << momentum + << " disk=" << disk + << " dmom=" << dmom + << " scatter_angle=" << scatter_angle << std::endl; + } + } + } diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 5e62035125..c64365176b 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -14,6 +14,9 @@ #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Run.h" #include "art/Framework/Principal/Handle.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileService.h" +#include "TH2F.h" // conditions #include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" #include "Offline/ProditionsService/inc/ProditionsHandle.hh" @@ -139,6 +142,7 @@ namespace mu2e { public: using Parameters = art::EDProducer::Table; explicit LoopHelixFit(const Parameters& settings); + void beginJob() override; void beginRun(art::Run& run) override; void produce(art::Event& event) override; void endJob() override; @@ -185,6 +189,16 @@ namespace mu2e { int nAmbiguous_ = 0; int nDownstream_ = 0; int nUpstream_ = 0; + // validation histograms for calorimeter material intersection + TH2F* h_intersection_efficiency_ = nullptr; + TH2F* h_frontpanel_hits_ = nullptr; + TH1F* h_momentum_degradation_ = nullptr; + TH1F* h_scatter_100mev_ = nullptr; + TH1F* h_scatter_80mev_ = nullptr; + TProfile* h_resolution_vs_momentum_ = nullptr; + TH2F* h_pid_e_vs_p_ = nullptr; + TH1F* h_pid_ep_ratio_ = nullptr; + TH2F* h_dedx_vs_momentum_ = nullptr; }; LoopHelixFit::LoopHelixFit(const Parameters& settings) : art::EDProducer{settings}, @@ -244,6 +258,53 @@ namespace mu2e { } } + void LoopHelixFit::beginJob() { + // Create validation histograms for calorimeter material effects + art::ServiceHandle tfs; + if(extrap_) { + h_intersection_efficiency_ = tfs->make("h_intersection_efficiency", + "Calorimeter front panel intersection efficiency;Momentum (MeV/c);Disk ID", + 100, 0.0, 150.0, 2, -0.5, 1.5); + + h_frontpanel_hits_ = tfs->make("h_frontpanel_hits", + "Track hits on calorimeter front panel;Phi (rad);Radius (mm)", + 32, 0.0, 2.0*M_PI, 60, 300.0, 600.0); + + h_momentum_degradation_ = tfs->make("h_momentum_degradation", + "Momentum degradation due to calorimeter material;dP (MeV/c);Count", + 100, 0.0, 2.0); + + h_scatter_100mev_ = tfs->make("h_scatter_100mev", + "Multiple scattering angle (100 MeV);Scattering angle (mrad);Count", + 100, 10.0, 50.0); + + h_scatter_80mev_ = tfs->make("h_scatter_80mev", + "Multiple scattering angle (80 MeV);Scattering angle (mrad);Count", + 100, 10.0, 50.0); + + h_resolution_vs_momentum_ = tfs->make("h_resolution_vs_momentum", + "Momentum resolution (relative error) vs. fitted momentum;Fitted Momentum (MeV/c);Relative Error (dP/P)", + 20, 50.0, 150.0); + + // Pass histograms to extrapolation module + extrap_->setValidationHistograms(h_intersection_efficiency_, h_frontpanel_hits_, h_momentum_degradation_, h_scatter_100mev_, h_scatter_80mev_, h_resolution_vs_momentum_); + } + + // Create PID discriminant histograms + h_pid_e_vs_p_ = tfs->make("h_pid_e_vs_p", + "Particle ID: Energy vs. Momentum;E/p ratio;Momentum (MeV/c)", + 50, 0.0, 2.0, 60, 50.0, 150.0); + + h_pid_ep_ratio_ = tfs->make("h_pid_ep_ratio", + "Particle ID: E/p ratio distribution;E/p ratio;Count", + 100, 0.0, 2.0); + + // Create energy loss vs momentum histogram (Plot 9) + h_dedx_vs_momentum_ = tfs->make("h_dedx_vs_momentum", + "Energy Loss vs. Momentum;Momentum (MeV/c);-dE/dx (MeV/mm)", + 60, 50.0, 150.0, 100, 0.0, 5.0); + } + void LoopHelixFit::beginRun(art::Run& run) { // setup things that rely on data related to beginRun auto const& ptable = GlobalConstantsHandle(); @@ -419,6 +480,30 @@ namespace mu2e { // convert to seed output format auto kkseed = kkfit_.createSeed(*ktrk,fitflag,*calo_h,*nominalTracker_h); if(print_>0) print_track_info(kkseed, *ktrk); + + // Fill PID histogram if calorimeter intersection exists + if((h_pid_e_vs_p_ || h_pid_ep_ratio_) && hseed.caloCluster().isNonnull()) { + auto const& seg = kkseed.segments()[0]; + double p_fit = seg.loopHelix().momentum(seg.tmin()); + double caloE = hseed.caloCluster()->energyDep(); + if(p_fit > 0.0) { + double ep_ratio = caloE / p_fit; + if(h_pid_e_vs_p_) h_pid_e_vs_p_->Fill(ep_ratio, p_fit); + if(h_pid_ep_ratio_) h_pid_ep_ratio_->Fill(ep_ratio); + } + } + + // Fill energy loss vs momentum histogram (Plot 9) + if(h_dedx_vs_momentum_ && hseed.caloCluster().isNonnull() && kkseed.segments().size() > 0) { + auto const& seg = kkseed.segments()[0]; + double p_fit = seg.loopHelix().momentum(seg.tmin()); + double E_loss = hseed.caloCluster()->energyDep(); // energy deposited + double thickness_mm = 24.75; // foam (21.75 mm) + carbon (3.0 mm) + if(p_fit > 0.0 && thickness_mm > 0.0 && E_loss > 0.0) { + double dedx = E_loss / thickness_mm; // MeV/mm + h_dedx_vs_momentum_->Fill(p_fit, dedx); + } + } kkseedcol->push_back(kkseed); // fill assns with the helix seed auto kseedptr = art::Ptr(KalSeedCollectionPID,kkseedcol->size()-1,KalSeedCollectionGetter); From bb79d78fcc3417980b07c4f58c28d6848e3606ce Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 14 May 2026 10:11:07 -0500 Subject: [PATCH 30/31] updates --- Mu2eKinKal/fcl/prolog.fcl | 10 ++- Mu2eKinKal/src/CentralHelixFit_module.cc | 94 +++++++++++++++++++++++- Mu2eKinKal/src/LoopHelixFit_module.cc | 33 +++++---- 3 files changed, 119 insertions(+), 18 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 2b3e63b35f..7488e8a4ac 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -397,7 +397,7 @@ Mu2eKinKal : { ToOPA : true ToCaloD0 : true ToCaloD1 : true - ToCaloMaterial : false + ToCaloMaterial : true } LHSEEDXTRAP : { @@ -410,7 +410,7 @@ Mu2eKinKal : { ToOPA : false ToCaloD0 : true ToCaloD1 : true - ToCaloMaterial : false + ToCaloMaterial : true } CENTRALHELIX : { @@ -455,7 +455,7 @@ Mu2eKinKal : { ToOPA : false ToCaloD0 : true ToCaloD1 : true - ToCaloMaterial : false + ToCaloMaterial : true IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } @@ -480,6 +480,7 @@ Mu2eKinKal : { } ExtrapolationSettings : @local::Mu2eKinKal.LHSEEDXTRAP UsePDGCharge: false + #MakeValidationPlots : false HelixMask: { MinHelixMom : 0 } @@ -505,6 +506,7 @@ Mu2eKinKal : { } ExtrapolationSettings : @local::Mu2eKinKal.LHDRIFTXTRAP UsePDGCharge: false + #MakeValidationPlots : false HelixMask: { MinHelixMom : 0 } @@ -562,6 +564,7 @@ Mu2eKinKal : { @table::Mu2eKinKal.CENTRALHELIX @table::Mu2eKinKal.KKPrecursors } + MakeValidationPlots : false } CHDriftFit : { @@ -578,6 +581,7 @@ Mu2eKinKal : { @table::Mu2eKinKal.CENTRALHELIX @table::Mu2eKinKal.KKPrecursors } + MakeValidationPlots : false } } diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index 0a90870f6d..23b19e49d9 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -14,6 +14,7 @@ #include "art/Framework/Principal/Event.h" #include "art/Framework/Principal/Run.h" #include "art/Framework/Principal/Handle.h" +#include "art_root_io/TFileService.h" // conditions #include "Offline/GlobalConstantsService/inc/GlobalConstantsHandle.hh" #include "Offline/ProditionsService/inc/ProditionsHandle.hh" @@ -64,6 +65,8 @@ #include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // root #include "TH1F.h" +#include "TH2F.h" +#include "TProfile.h" #include "TTree.h" #include "Math/AxisAngle.h" // C++ @@ -133,6 +136,7 @@ namespace mu2e { fhicl::Table matSettings { Name("MaterialSettings") }; // helix module specific config fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; + fhicl::Atom makeValidationPlots { Name("MakeValidationPlots"), Comment("Enable creation of validation plots for calorimeter material effects") }; }; class CentralHelixFit : public art::EDProducer { @@ -140,6 +144,7 @@ namespace mu2e { using Parameters = art::EDProducer::Table; explicit CentralHelixFit(const Parameters& settings); virtual ~CentralHelixFit() {} + void beginJob() override; void beginRun(art::Run& run) override; void produce(art::Event& event) override; protected: @@ -179,6 +184,17 @@ namespace mu2e { std::array paramconstraints_; bool extrapolate_; std::unique_ptr extrap_; // extrapolations + bool makeValidationPlots_ = false; + // validation histograms for calorimeter material intersection + TH2F* h_intersection_efficiency_ = nullptr; + TH2F* h_frontpanel_hits_ = nullptr; + TH1F* h_momentum_degradation_ = nullptr; + TH1F* h_scatter_100mev_ = nullptr; + TH1F* h_scatter_80mev_ = nullptr; + TProfile* h_resolution_vs_momentum_ = nullptr; + TH2F* h_pid_e_vs_p_ = nullptr; + TH1F* h_pid_ep_ratio_ = nullptr; + TH2F* h_dedx_vs_momentum_ = nullptr; }; CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings}, @@ -202,7 +218,8 @@ namespace mu2e { minCenterRho_(settings().modSettings().minCenterRho()), sampleinrange_(settings().modSettings().sampleInRange()), sampleinbounds_(settings().modSettings().sampleInBounds()), - extrapolate_(false) + extrapolate_(false), + makeValidationPlots_(settings().makeValidationPlots()) { // collection handling for(const auto& cseedtag : settings().modSettings().seedCollections()) { cseedCols_.emplace_back(consumes(cseedtag)); } @@ -242,6 +259,55 @@ namespace mu2e { } } + void CentralHelixFit::beginJob() { + // Create validation histograms for calorimeter material effects + art::ServiceHandle tfs; + if(makeValidationPlots_ && extrap_) { + h_intersection_efficiency_ = tfs->make("h_intersection_efficiency", + "Calorimeter front panel intersection efficiency;Momentum (MeV/c);Disk ID", + 100, 0.0, 150.0, 2, -0.5, 1.5); + + h_frontpanel_hits_ = tfs->make("h_frontpanel_hits", + "Track hits on calorimeter front panel;Phi (rad);Radius (mm)", + 32, 0.0, 2.0*M_PI, 60, 300.0, 600.0); + + h_momentum_degradation_ = tfs->make("h_momentum_degradation", + "Momentum degradation due to calorimeter material;dP (MeV/c);Count", + 100, 0.0, 2.0); + + h_scatter_100mev_ = tfs->make("h_scatter_100mev", + "Multiple scattering angle (100 MeV);Scattering angle (mrad);Count", + 100, 10.0, 50.0); + + h_scatter_80mev_ = tfs->make("h_scatter_80mev", + "Multiple scattering angle (80 MeV);Scattering angle (mrad);Count", + 100, 10.0, 50.0); + + h_resolution_vs_momentum_ = tfs->make("h_resolution_vs_momentum", + "Momentum resolution (relative error) vs. fitted momentum;Fitted Momentum (MeV/c);Relative Error (dP/P)", + 20, 50.0, 150.0); + + // Pass histograms to extrapolation module + extrap_->setValidationHistograms(h_intersection_efficiency_, h_frontpanel_hits_, h_momentum_degradation_, h_scatter_100mev_, h_scatter_80mev_, h_resolution_vs_momentum_); + } + + if(makeValidationPlots_) { + // Create PID discriminant histograms + h_pid_e_vs_p_ = tfs->make("h_pid_e_vs_p", + "Particle ID: Energy vs. Momentum;E/p ratio;Momentum (MeV/c)", + 50, 0.0, 2.0, 60, 50.0, 150.0); + + h_pid_ep_ratio_ = tfs->make("h_pid_ep_ratio", + "Particle ID: E/p ratio distribution;E/p ratio;Count", + 100, 0.0, 2.0); + + // Create energy loss vs momentum histogram + h_dedx_vs_momentum_ = tfs->make("h_dedx_vs_momentum", + "Energy Loss vs. Momentum;Momentum (MeV/c);-dE/dx (MeV/mm)", + 60, 50.0, 150.0, 100, 0.0, 5.0); + } + } + void CentralHelixFit::beginRun(art::Run& run) { // setup things that rely on data related to beginRun auto const& ptable = GlobalConstantsHandle(); @@ -384,7 +450,33 @@ namespace mu2e { if(!degen){ sampleFit(*kktrk); auto kkseed = kkfit_.createSeed(*kktrk,fitflag,*calo_h,*nominalTracker_h); + // Fill PID histograms if enabled + if(makeValidationPlots_ && cc_H.isValid() && cc_H->size() > 0) { + auto const& cluster = cc_H->at(0); + if(makeValidationPlots_ && kkseed.segments().size() > 0) { + auto const& seg = kkseed.segments()[0]; + double p_fit = seg.loopHelix().momentum(seg.tmin()); + double E_cluster = cluster.energyDep(); + if(p_fit > 0.0 && E_cluster > 0.0) { + double ep_ratio = E_cluster / p_fit; + if(h_pid_e_vs_p_) h_pid_e_vs_p_->Fill(ep_ratio, p_fit); + if(h_pid_ep_ratio_) h_pid_ep_ratio_->Fill(ep_ratio); + } + } + } kkseedcol->push_back(kkseed); + // Fill energy loss vs momentum histogram + if(makeValidationPlots_ && cc_H.isValid() && cc_H->size() > 0 && kkseed.segments().size() > 0) { + auto const& cluster = cc_H->at(0); + auto const& seg = kkseed.segments()[0]; + double p_fit = seg.loopHelix().momentum(seg.tmin()); + double E_loss = cluster.energyDep(); // energy deposited + double thickness_mm = 24.75; // foam (21.75 mm) + carbon (3.0 mm) + if(p_fit > 0.0 && thickness_mm > 0.0 && E_loss > 0.0) { + double dedx = E_loss / thickness_mm; // MeV/mm + if(h_dedx_vs_momentum_) h_dedx_vs_momentum_->Fill(p_fit, dedx); + } + } // save (unpersistable) KKTrk in the event kktrkcol->push_back(kktrk.release()); } diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index c64365176b..ed31614c96 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -136,6 +136,7 @@ namespace mu2e { fhicl::OptionalAtom fitDirection { Name("FitDirection"), Comment("Particle direction to fit, either \"upstream\" or \"downstream\"")}; fhicl::Atom pdgCharge { Name("UsePDGCharge"), Comment("Use particle charge from fitParticle")}; fhicl::OptionalTable HelixMask { Name("HelixMask"), Comment("Selections applied to input helices")}; + //fhicl::Atom makeValidationPlots { Name("MakeValidationPlots"), Comment("Enable creation of validation plots for calorimeter material effects") }; }; class LoopHelixFit : public art::EDProducer { @@ -189,6 +190,7 @@ namespace mu2e { int nAmbiguous_ = 0; int nDownstream_ = 0; int nUpstream_ = 0; + bool makeValidationPlots_ = false; // validation histograms for calorimeter material intersection TH2F* h_intersection_efficiency_ = nullptr; TH2F* h_frontpanel_hits_ = nullptr; @@ -215,7 +217,8 @@ namespace mu2e { kkmat_(settings().matSettings()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), - fixedfield_(false) + fixedfield_(false)//, + //makeValidationPlots_(settings().makeValidationPlots()) { std::string fdir; if(settings().fitDirection(fdir))fdir_ = fdir; @@ -261,7 +264,7 @@ namespace mu2e { void LoopHelixFit::beginJob() { // Create validation histograms for calorimeter material effects art::ServiceHandle tfs; - if(extrap_) { + if(makeValidationPlots_ && extrap_) { h_intersection_efficiency_ = tfs->make("h_intersection_efficiency", "Calorimeter front panel intersection efficiency;Momentum (MeV/c);Disk ID", 100, 0.0, 150.0, 2, -0.5, 1.5); @@ -290,19 +293,21 @@ namespace mu2e { extrap_->setValidationHistograms(h_intersection_efficiency_, h_frontpanel_hits_, h_momentum_degradation_, h_scatter_100mev_, h_scatter_80mev_, h_resolution_vs_momentum_); } - // Create PID discriminant histograms - h_pid_e_vs_p_ = tfs->make("h_pid_e_vs_p", - "Particle ID: Energy vs. Momentum;E/p ratio;Momentum (MeV/c)", - 50, 0.0, 2.0, 60, 50.0, 150.0); + if(makeValidationPlots_) { + // Create PID discriminant histograms + h_pid_e_vs_p_ = tfs->make("h_pid_e_vs_p", + "Particle ID: Energy vs. Momentum;E/p ratio;Momentum (MeV/c)", + 50, 0.0, 2.0, 60, 50.0, 150.0); - h_pid_ep_ratio_ = tfs->make("h_pid_ep_ratio", - "Particle ID: E/p ratio distribution;E/p ratio;Count", - 100, 0.0, 2.0); + h_pid_ep_ratio_ = tfs->make("h_pid_ep_ratio", + "Particle ID: E/p ratio distribution;E/p ratio;Count", + 100, 0.0, 2.0); - // Create energy loss vs momentum histogram (Plot 9) - h_dedx_vs_momentum_ = tfs->make("h_dedx_vs_momentum", - "Energy Loss vs. Momentum;Momentum (MeV/c);-dE/dx (MeV/mm)", - 60, 50.0, 150.0, 100, 0.0, 5.0); + // Create energy loss vs momentum histogram + h_dedx_vs_momentum_ = tfs->make("h_dedx_vs_momentum", + "Energy Loss vs. Momentum;Momentum (MeV/c);-dE/dx (MeV/mm)", + 60, 50.0, 150.0, 100, 0.0, 5.0); + } } void LoopHelixFit::beginRun(art::Run& run) { @@ -493,7 +498,7 @@ namespace mu2e { } } - // Fill energy loss vs momentum histogram (Plot 9) + // Fill energy loss vs momentum histogram if(h_dedx_vs_momentum_ && hseed.caloCluster().isNonnull() && kkseed.segments().size() > 0) { auto const& seg = kkseed.segments()[0]; double p_fit = seg.loopHelix().momentum(seg.tmin()); From cfb40c2f8de889b5871b7eebe0f60c3683e508a0 Mon Sep 17 00:00:00 2001 From: sophieMu2e Date: Thu, 14 May 2026 10:34:40 -0500 Subject: [PATCH 31/31] updates --- Mu2eKinKal/inc/KKExtrap.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 1704356b17..8c544a90a9 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -67,7 +67,7 @@ namespace mu2e { mutable ExtrapolateCaloMaterial extrapCaloMat_; // extrapolation through passive calo materials (disk 0) mutable ExtrapolateCaloMaterial extrapCaloMatD1_; // extrapolation through passive calo materials (disk 1) mutable bool geom_initialized_ = false; - void initGeometry() const; // lazy initialization method + void initGeometry() const; // new initialization method double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO // calorimeter front panel passive material thicknesses (should come from geometry service TODO) @@ -153,7 +153,7 @@ namespace mu2e { } template void KKExtrap::extrapolate(KKTrack& ktrk) const { - initGeometry(); // lazy initialization: initialize geometry on first use + initGeometry(); // initialize geometry on first use GeomHandle kkg_h; auto const& kkg = *kkg_h; ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_);