diff --git a/CommonMC/src/MakeSurfaceSteps_module.cc b/CommonMC/src/MakeSurfaceSteps_module.cc index b82d0095b6..f8e07bd676 100644 --- a/CommonMC/src/MakeSurfaceSteps_module.cc +++ b/CommonMC/src/MakeSurfaceSteps_module.cc @@ -60,7 +60,17 @@ 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"); + + 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) { GeomHandle det; @@ -72,6 +82,7 @@ namespace mu2e { for(auto const& vdspmc : vdspmccol) { // only some VDs are kept auto isid = vdmap_.find(vdspmc.virtualDetectorId()); + if(debug_ > 0) std::cout<<" VID "<emplace_back(isid->second,vdspmc,det); // no aggregation of VD hits } auto nvdsteps = ssc->size(); diff --git a/DataProducts/inc/SurfaceId.hh b/DataProducts/inc/SurfaceId.hh index d70c661bb8..000ff49c28 100644 --- a/DataProducts/inc/SurfaceId.hh +++ b/DataProducts/inc/SurfaceId.hh @@ -20,7 +20,10 @@ 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 + 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_FrontPanel = 320 // calorimeter front panel passive materials (foam, carbon) }; static std::string const& typeName(); @@ -38,7 +41,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 4d454e0922..512d3748b7 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -34,7 +34,17 @@ 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::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"), + 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::make_pair(SurfaceIdEnum::EMC_FrontPanel, "EMC_FrontPanel") }; std::map const& SurfaceIdDetail::names(){ return nam; 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..b6cfdf5fe1 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,62 @@ 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) + 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 + 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; + + // 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 +121,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 +138,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 +163,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/CMakeLists.txt b/KinKalGeom/CMakeLists.txt index c600163619..985f3f48eb 100644 --- a/KinKalGeom/CMakeLists.txt +++ b/KinKalGeom/CMakeLists.txt @@ -1,6 +1,7 @@ cet_make_library( SOURCE src/CRV.cc + src/Calo.cc src/KinKalGeom.cc LIBRARIES PUBLIC KinKal::Geometry diff --git a/KinKalGeom/inc/Calo.hh b/KinKalGeom/inc/Calo.hh new file mode 100644 index 0000000000..571e75bbbf --- /dev/null +++ b/KinKalGeom/inc/Calo.hh @@ -0,0 +1,56 @@ +// +// 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_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 KKGeom { + class Calo { + public: + using CylPtr = std::shared_ptr; + using DiskPtr = std::shared_ptr; + // 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_;} + 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_;} + + // 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 + }; + } +} + +#endif diff --git a/KinKalGeom/inc/KinKalGeom.hh b/KinKalGeom/inc/KinKalGeom.hh index 05770e6730..b6dc492cf1 100644 --- a/KinKalGeom/inc/KinKalGeom.hh +++ b/KinKalGeom/inc/KinKalGeom.hh @@ -10,6 +10,7 @@ #include "Offline/KinKalGeom/inc/DetectorSolenoid.hh" #include "Offline/KinKalGeom/inc/CRV.hh" #include "Offline/KinKalGeom/inc/TestCRV.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" #include "KinKal/Geometry/Surface.hh" #include "Offline/Mu2eInterfaces/inc/Detector.hh" @@ -33,17 +34,18 @@ namespace mu2e { void surfaces(SurfaceId const& sid, std::vector& surfs) const; // find all surfaces that match a vector of Ids void surfaces(std::vector const& ids, std::vector& surfs) const; - // hierarchical accessors - auto const& DS() const {return ds_; } - auto const& ST() const {return st_; } - auto const& tracker() const {return tracker_; } -// auto const& CRV() const {return crv_; } - auto const& TCRV() const {return tcrv_; } + // hierarchical accessors: return the underlying objects, not the unique_ptrs + 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 std::unique_ptr tracker_; std::unique_ptr ds_; std::unique_ptr st_; + std::unique_ptr calo_; //KKGeom::CRV crv_; std::unique_ptr tcrv_; // the map used to find surfaces by Id diff --git a/KinKalGeom/src/Calo.cc b/KinKalGeom/src/Calo.cc new file mode 100644 index 0000000000..6c7152162a --- /dev/null +++ b/KinKalGeom/src/Calo.cc @@ -0,0 +1,22 @@ +#include "Offline/KinKalGeom/inc/Calo.hh" +namespace mu2e { + namespace KKGeom { + using KinKal::VEC3; + using KinKal::Cylinder; + using KinKal::Disk; + + 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,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)}, + 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/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 c827389d0f..7488e8a4ac 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 @@ -391,6 +395,9 @@ Mu2eKinKal : { Upstream : true BackToTracker : false ToOPA : true + ToCaloD0 : true + ToCaloD1 : true + ToCaloMaterial : true } LHSEEDXTRAP : { @@ -401,6 +408,9 @@ Mu2eKinKal : { Upstream : false BackToTracker : false ToOPA : false + ToCaloD0 : true + ToCaloD1 : true + ToCaloMaterial : true } CENTRALHELIX : { @@ -435,9 +445,18 @@ Mu2eKinKal : { } LINEEXTRAPOLATION : { - MaxDt : 200.0 # (ns) - MinV : 1e-5 - ToCRV : true + Debug : 2 + MaxDt : 200.0 + BCorrTolerance : 1e-5 + ToTrackerEnds : true + Upstream : true + BackToTracker : false + ToTrackerEnds : false + ToOPA : false + ToCaloD0 : true + ToCaloD1 : true + ToCaloMaterial : true + IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } @@ -461,6 +480,7 @@ Mu2eKinKal : { } ExtrapolationSettings : @local::Mu2eKinKal.LHSEEDXTRAP UsePDGCharge: false + #MakeValidationPlots : false HelixMask: { MinHelixMom : 0 } @@ -486,6 +506,7 @@ Mu2eKinKal : { } ExtrapolationSettings : @local::Mu2eKinKal.LHDRIFTXTRAP UsePDGCharge: false + #MakeValidationPlots : false HelixMask: { MinHelixMom : 0 } @@ -543,6 +564,7 @@ Mu2eKinKal : { @table::Mu2eKinKal.CENTRALHELIX @table::Mu2eKinKal.KKPrecursors } + MakeValidationPlots : false } CHDriftFit : { @@ -559,6 +581,7 @@ Mu2eKinKal : { @table::Mu2eKinKal.CENTRALHELIX @table::Mu2eKinKal.KKPrecursors } + MakeValidationPlots : false } } diff --git a/Mu2eKinKal/inc/ExtrapolateCalo.hh b/Mu2eKinKal/inc/ExtrapolateCalo.hh new file mode 100644 index 0000000000..1906c73dd7 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCalo.hh @@ -0,0 +1,161 @@ +// 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 KKGeom::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]; + 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; + } + 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/ExtrapolateCaloMaterial.hh b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh new file mode 100644 index 0000000000..d84e583c0c --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCaloMaterial.hh @@ -0,0 +1,173 @@ +// predicate to extrapolate to passive calorimeter front panel materials (foam, carbon) +// Sophie Middleton (2026) +#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 +#include "TH2F.h" +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()), disk_(0), + 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; } + + // 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 + + private: + double maxDt_; // maximum extrapolation time + 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 { + // 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(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; + } + return false; // stop, we found material crossing + } else { + if(debug_ == -300) { + std::cout << " No intersection found with annulus in Z-range" << std::endl; + } + return false; // stop - panel in range but no intersection means we're past it + } + } +} +#endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index ab5a4205c7..8c544a90a9 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -12,12 +12,16 @@ #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/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" +#include "TProfile.h" namespace mu2e { using KinKal::VEC3; using KinKal::TimeDir; @@ -37,16 +41,45 @@ 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; 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_; double btol_, intertol_, maxdt_; KKMaterial const& kkmat_; + AnnPtr tsdaptr_; + DiskPtr trkfrontptr_, trkmidptr_, trkbackptr_; + FruPtr opaptr_; bool backToTracker_, toOPA_, toTrackerEnds_, upstream_; + 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 (disk 0) + mutable ExtrapolateCaloMaterial extrapCaloMatD1_; // extrapolation through passive calo materials (disk 1) + mutable bool geom_initialized_ = false; + 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) + 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) : @@ -55,29 +88,97 @@ namespace mu2e { intertol_(extrapconfig.interTol()), maxdt_(extrapconfig.MaxDt()), kkmat_(kkmat), + tsdaptr_(nullptr), + trkfrontptr_(nullptr), + trkmidptr_(nullptr), + trkbackptr_(nullptr), + opaptr_(nullptr), backToTracker_(extrapconfig.BackToTracker()), toOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), - upstream_(extrapconfig.Upstream()) - {} + upstream_(extrapconfig.Upstream()), + toCaloD0_(extrapconfig.ToCaloD0()), + toCaloD1_(extrapconfig.ToCaloD1()), + toCaloMaterial_(extrapconfig.ToCaloMaterial()), + geom_initialized_(false) + { + // geometry initialization is deferred to first use via initGeometry() + } + + inline void KKExtrap::initGeometry() const { + if(geom_initialized_) return; // already initialized + GeomHandle kkg_h; + auto const& kkg = *kkg_h; + + const_cast(tsdaptr_) = kkg.DS().upstreamAbsorberPtr(); + const_cast(trkfrontptr_) = kkg.tracker().frontPtr(); + const_cast(trkmidptr_) = kkg.tracker().middlePtr(); + const_cast(trkbackptr_) = kkg.tracker().backPtr(); + const_cast(opaptr_) = kkg.DS().outerProtonAbsorberPtr(); + + // calo surfaces + const_cast(calod0frontptr_) = kkg.calo().EMC_Disk_0_FrontPtr(); + const_cast(calod0backptr_) = kkg.calo().EMC_Disk_0_BackPtr(); + const_cast(calod1frontptr_) = kkg.calo().EMC_Disk_1_FrontPtr(); + const_cast(calod1backptr_) = kkg.calo().EMC_Disk_1_BackPtr(); + + // initialize predicates + const_cast(TSDA_) = ExtrapolateToZ(maxdt_,btol_,kkg.DS().upstreamAbsorber().center().Z(),debug_); + const_cast(trackerFront_) = ExtrapolateToZ(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); + 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 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_); + 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; + const_cast(geom_initialized_) = true; + } template void KKExtrap::extrapolate(KKTrack& ktrk) const { - GeomHandle kkg_h; + initGeometry(); // initialize geometry on first use + GeomHandle kkg_h; auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); // 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(toCaloD1_)toCaloD1(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) + 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 @@ -99,14 +200,29 @@ 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 { - GeomHandle kkg_h; + initGeometry(); + GeomHandle kkg_h; auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); - ExtrapolateToZ trackerBack(maxdt_,btol_,kkg.tracker()->back().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); + ExtrapolateToZ trackerBack(maxdt_,btol_,kkg.tracker().back().center().Z(),debug_); // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't // reflect inside the tracker volume auto const& ftraj = ktrk.fitTraj(); @@ -121,18 +237,17 @@ namespace mu2e { static const SurfaceId tt_back("TT_Back"); // start with the middle - auto midinter = KinKal::intersect(ftraj,kkg.tracker()->middle(),ftraj.range(),intertol_); + 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); + auto frontinter = KinKal::intersect(fhel,kkg.tracker().front(),fhel.range(),intertol_,fronttdir); if(!frontinter.good()){ // start from the middle TimeRange frange = ftraj.range(); if(midinter.good())frange = fronttdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - frontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),frange,intertol_,fronttdir); + frontinter = KinKal::intersect(ftraj,kkg.tracker().front(),frange,intertol_,fronttdir); } if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); } @@ -140,19 +255,92 @@ namespace mu2e { // start from the middle TimeRange brange = ftraj.range(); if(midinter.good())brange = backtdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - auto backinter = KinKal::intersect(ftraj,kkg.tracker()->back(),brange,intertol_,backtdir); + auto backinter = KinKal::intersect(ftraj,kkg.tracker().back(),brange,intertol_,backtdir); if(backinter.good())ktrk.addIntersection(tt_back,backinter); } } + 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; + if(debug_ > 0){ + std::cout<<"toCaloD0 DEBUG:"< 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_); + // record the standard tracker intersections + static const SurfaceId d1_front("EMC_Disk_1_Front"); + static const SurfaceId d1_back("EMC_Disk_1_Back"); + static const SurfaceId d1_inner("EMC_Disk_1_Inner"); + static const SurfaceId d1_outer("EMC_Disk_1_Outer"); + + if(tocalofront){ + TimeRange frange = ftraj.range(); + auto frontinter = KinKal::intersect(ftraj,*calod1frontptr_,frange,intertol_,fronttdir); + if(frontinter.good()) ktrk.addIntersection(d1_front,frontinter); + } + if(tocaloback){ + TimeRange brange = ftraj.range(); + auto backinter = KinKal::intersect(ftraj,*calod1backptr_,brange,intertol_,backtdir); + if(backinter.good())ktrk.addIntersection(d1_back,backinter); + } + // get intersections with inner and outer cylinders + auto innerinter = KinKal::intersect(ftraj,kkg.calo().EMC_Disk_1_Inner(),ftraj.range(),intertol_); + if(innerinter.good()) ktrk.addIntersection(d1_inner,innerinter); + auto outerinter = KinKal::intersect(ftraj,kkg.calo().EMC_Disk_1_Outer(),ftraj.range(),intertol_); + if(outerinter.good()) ktrk.addIntersection(d1_outer,outerinter); + } + template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { - GeomHandle kkg_h; + initGeometry(); + GeomHandle kkg_h; auto const& kkg = *kkg_h; using KKIPAXING = KKShellXing; using KKIPAXINGPTR = std::shared_ptr; // extraplate the fit through the IPA. This will add material effects for each intersection. It will continue till the // track exits the IPA - ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg.DS()->innerProtonAbsorberPtr(),debug_); + ExtrapolateIPA extrapIPA(maxdt_,btol_,intertol_,kkg.DS().innerProtonAbsorberPtr(),debug_); if(extrapIPA.debug() > 0)std::cout << "extrapolating to IPA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId IPASID("IPA"); @@ -163,7 +351,7 @@ namespace mu2e { if(extrapIPA.intersection().good()){ // we have a good intersection. Use this to create a Shell material Xing auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - auto const& IPA = kkg.DS()->innerProtonAbsorberPtr(); + auto const& IPA = kkg.DS().innerProtonAbsorberPtr(); KKIPAXINGPTR ipaxingptr = std::make_shared(IPA,IPASID,*kkmat_.IPAMaterial(),extrapIPA.intersection(),reftrajptr,ipathick_,extrapIPA.interTolerance()); if(extrapIPA.debug() > 0){ double dmom, paramomvar, perpmomvar; @@ -190,12 +378,13 @@ namespace mu2e { template bool KKExtrap::extrapolateST(KKTrack& ktrk,TimeDir tdir) const { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; - GeomHandle kkg_h; + initGeometry(); + 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 // track exits the ST in Z - ExtrapolateST extrapST(maxdt_,btol_,intertol_,*kkg.ST(),debug_); + ExtrapolateST extrapST(maxdt_,btol_,intertol_,kkg.ST(),debug_); auto const& ftraj = ktrk.fitTraj(); double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); auto startdir = ftraj.direction(starttime); @@ -230,16 +419,17 @@ namespace mu2e { } template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { - GeomHandle kkg_h; + initGeometry(); + GeomHandle kkg_h; auto const& kkg = *kkg_h; - ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker()->front().center().Z(),debug_); + ExtrapolateToZ trackerFront(maxdt_,btol_,kkg.tracker().front().center().Z(),debug_); if(trackerFront.debug() > 0)std::cout << "extrapolating to Tracker " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TrackerSID("TT_Front"); ktrk.extrapolate(tdir,trackerFront); // the last piece appended should cover the necessary range auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto trkfrontinter = KinKal::intersect(ftraj,kkg.tracker()->front(),ktraj.range(),intertol_,tdir); + auto trkfrontinter = KinKal::intersect(ftraj,kkg.tracker().front(),ktraj.range(),intertol_,tdir); if(trkfrontinter.onsurface_){ // dont worry about bounds here ktrk.addIntersection(TrackerSID,trkfrontinter); return true; @@ -247,10 +437,153 @@ 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) 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 BOTH calo front panel disks" << 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); + + // 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 === Disk " << disk_id << " Crossing ===" << std::endl; + std::cout << " Found intersection with calo front panel Annulus" << std::endl; + 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 = 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; + } + + // Create the material xing using combined front panel thickness + KKCALOMATXINGPTR matxingptr = std::make_shared( + annulusptr, CaloMatSID, *kkmat_.CaloFrontFoamMaterial(), + 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){ + 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 + 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 + 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) { + 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 { - GeomHandle kkg_h; + initGeometry(); + GeomHandle kkg_h; auto const& kkg = *kkg_h; - ExtrapolateToZ TSDA(maxdt_,btol_,kkg.DS()->upstreamAbsorber().center().Z(),debug_); + ExtrapolateToZ TSDA(maxdt_,btol_,kkg.DS().upstreamAbsorber().center().Z(),debug_); if(TSDA.debug() > 0)std::cout << "extrapolating to TSDA " << std::endl; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId TSDASID("TSDA"); @@ -261,22 +594,38 @@ namespace mu2e { bool retval = epos.Z() < TSDA.zVal(); if(retval){ auto const& ktraj = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto tsdainter = KinKal::intersect(ftraj,kkg.DS()->upstreamAbsorber(),ktraj.range(),intertol_,tdir); + auto tsdainter = KinKal::intersect(ftraj,kkg.DS().upstreamAbsorber(),ktraj.range(),intertol_,tdir); if(tsdainter.onsurface_)ktrk.addIntersection(TSDASID,tsdainter); } return retval; } template void KKExtrap::toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const { - GeomHandle kkg_h; + initGeometry(); + GeomHandle kkg_h; auto const& kkg = *kkg_h; auto const& ftraj = ktrk.fitTraj(); static const SurfaceId OPASID("OPA"); TimeRange trange = tdir == TimeDir::forwards ? TimeRange(tstart,ftraj.range().end()) : TimeRange(ftraj.range().begin(),tstart); - auto opainter = KinKal::intersect(ftraj,kkg.DS()->outerProtonAbsorber(),trange,intertol_,tdir); + auto opainter = KinKal::intersect(ftraj,kkg.DS().outerProtonAbsorber(),trange,intertol_,tdir); if(opainter.good()){ 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/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 3a881c64f6..ac7c0cb015 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -114,6 +114,9 @@ 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 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/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index d45316dd45..23b19e49d9 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -7,12 +7,14 @@ #include "fhiclcpp/types/Atom.h" #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/Table.h" +#include "fhiclcpp/types/OptionalTable.h" #include "fhiclcpp/types/Tuple.h" #include "fhiclcpp/types/OptionalAtom.h" #include "art/Framework/Core/EDProducer.h" #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" @@ -60,8 +62,11 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" #include "Offline/Mu2eKinKal/inc/KKConstantBField.hh" #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // root #include "TH1F.h" +#include "TH2F.h" +#include "TProfile.h" #include "TTree.h" #include "Math/AxisAngle.h" // C++ @@ -130,6 +135,8 @@ namespace mu2e { fhicl::Table extSettings { Name("ExtensionSettings") }; 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 { @@ -137,11 +144,13 @@ 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: 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 @@ -171,9 +180,22 @@ 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 + 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}, fitflag_(TrkFitFlag::KKCentralHelix), @@ -195,7 +217,9 @@ namespace mu2e { useFitCharge_(settings().modSettings().useFitCharge()), minCenterRho_(settings().modSettings().minCenterRho()), sampleinrange_(settings().modSettings().sampleInRange()), - sampleinbounds_(settings().modSettings().sampleInBounds()) + sampleinbounds_(settings().modSettings().sampleInBounds()), + extrapolate_(false), + makeValidationPlots_(settings().makeValidationPlots()) { // collection handling for(const auto& cseedtag : settings().modSettings().seedCollections()) { cseedCols_.emplace_back(consumes(cseedtag)); } @@ -223,8 +247,67 @@ namespace mu2e { for(auto const& sidname : settings().modSettings().sampleSurfaces()) { 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 + KinKalGeom smap; + smap.surfaces(ssids_,surfaces_to_sample_); + // configure extrapolation + if(settings().Extrapolation()){ + extrapolate_ = true; + // create KKExtrap + extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); + } + } + + 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(); @@ -240,7 +323,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 ) { @@ -340,6 +423,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_); @@ -365,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()); } @@ -400,7 +511,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); @@ -424,6 +535,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) diff --git a/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc new file mode 100644 index 0000000000..98c225298a --- /dev/null +++ b/Mu2eKinKal/src/ExtrapolateCaloMaterial.cc @@ -0,0 +1,108 @@ +// 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; + 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()), + 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(); + + // 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; + } + } + + // 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/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/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index c973b93ba9..2dd573ba10 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -53,6 +53,7 @@ #include "Offline/Mu2eKinKal/inc/KKFitUtilities.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateTCRV.hh" +#include "Offline/Mu2eKinKal/inc/KKExtrap.hh" // root #include "TH1F.h" #include "TTree.h" @@ -120,21 +121,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\"")}; }; @@ -172,12 +165,13 @@ 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_; 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 Config config_; // initial fit configuration object Config exconfig_; // extension configuration object }; @@ -224,17 +218,18 @@ namespace mu2e { }else{ throw cet::exception("RECO")<<"mu2e::KinematicLineFit: Parameter constraint configuration error"<< endl; } + // 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 } // configure extrapolation if(settings().Extrapolation()){ extrapolate_ = true; - toCRV_ = settings().Extrapolation()->ToCRV(); + // create KKExtrap for calorimeter and upstream extrapolations + extrap_ = std::make_unique(*settings().Extrapolation(),kkmat_); // 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(); } @@ -258,7 +253,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_,sample_); } void KinematicLineFit::produce(art::Event& event ) { @@ -384,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); @@ -413,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"); @@ -440,6 +435,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) diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 5e62035125..ed31614c96 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" @@ -133,12 +136,14 @@ 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 { 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 +190,17 @@ 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; + 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}, @@ -201,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; @@ -244,6 +261,55 @@ namespace mu2e { } } + void LoopHelixFit::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 LoopHelixFit::beginRun(art::Run& run) { // setup things that rely on data related to beginRun auto const& ptable = GlobalConstantsHandle(); @@ -419,6 +485,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 + 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); 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 diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index 49aea5119f..d4ab01dfaa 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/TrkReco/fcl/prolog.fcl @@ -8,6 +8,7 @@ BEGIN_PROLOG # # may use a more careful optimization, but this is it for now #------------------------------------------------------------------------------ + TrkReco: { McUtils : { tool_type : "TrkRecoMcUtils" comboHitCollTag : "makeSH"