diff --git a/CommonMC/src/MakeSurfaceSteps_module.cc b/CommonMC/src/MakeSurfaceSteps_module.cc index b82d0095b6..0d015c4370 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..e870cd3664 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 + }; 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..484a74af4f 100644 --- a/DataProducts/src/SurfaceId.cc +++ b/DataProducts/src/SurfaceId.cc @@ -34,7 +34,15 @@ namespace mu2e { std::make_pair(SurfaceIdEnum::ST_Inner, "ST_Inner"), std::make_pair(SurfaceIdEnum::ST_Outer, "ST_Outer"), std::make_pair(SurfaceIdEnum::ST_Foils, "ST_Foils"), - std::make_pair(SurfaceIdEnum::TCRV, "TCRV") + std::make_pair(SurfaceIdEnum::TCRV, "TCRV"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Outer, "EMC_Disk_0_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Inner, "EMC_Disk_0_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Inner, "EMC_Disk_1_Inner"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Outer, "EMC_Disk_1_Outer"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Front, "EMC_Disk_0_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Front, "EMC_Disk_1_Front"), + std::make_pair(SurfaceIdEnum::EMC_Disk_0_Back, "EMC_Disk_0_Back"), + std::make_pair(SurfaceIdEnum::EMC_Disk_1_Back, "EMC_Disk_1_Back") }; std::map const& SurfaceIdDetail::names(){ return nam; diff --git a/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/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index c827389d0f..d0cb140dab 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -391,6 +391,8 @@ Mu2eKinKal : { Upstream : true BackToTracker : false ToOPA : true + ToCaloD0 : true + ToCaloD1 : true } LHSEEDXTRAP : { @@ -401,6 +403,8 @@ Mu2eKinKal : { Upstream : false BackToTracker : false ToOPA : false + ToCaloD0 : true + ToCaloD1 : true } CENTRALHELIX : { @@ -435,9 +439,15 @@ Mu2eKinKal : { } LINEEXTRAPOLATION : { - MaxDt : 200.0 # (ns) - MinV : 1e-5 - ToCRV : true + MaxDt : 200.0 + BCorrTolerance : 1e-5 + Upstream : true + BackToTracker : false + ToTrackerEnds : true + ToOPA : false + ToCaloD0 : true + ToCaloD1 : true + IntersectionTolerance : 0.1 # tolerance for intersections (mm) } } diff --git a/Mu2eKinKal/inc/ExtrapolateCalo.hh b/Mu2eKinKal/inc/ExtrapolateCalo.hh new file mode 100644 index 0000000000..d9d3aa7bd2 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCalo.hh @@ -0,0 +1,160 @@ +// predicate to extrapolate to calo +// Sophie Middleton (2025) +#ifndef Mu2eKinKal_ExtrapolateCalo_hh +#define Mu2eKinKal_ExtrapolateCalo_hh +#include "Offline/Mu2eKinKal/inc/KKTrack.hh" +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/Calo.hh" +#include "cetlib_except/exception.h" +#include +#include +#include +namespace mu2e { + using KinKal::TimeDir; + using KinKal::TimeRange; + using 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]; + //sid_ = SurfaceId(SurfaceIdEnum::calo_Foils,idisk); //FIXME + if(debug_ > 0)std::cout << "Good calo disk " << newinter << " sid " << sid_ << std::endl; + return false; + } + idisk += ddisk; // otherwise continue loopin on disks + } + // no more intersections: keep extending in Z till we clear the calo + reset(); + if(debug_ > 1)std::cout << "Extrapolating to calo edge, z " << spos.Z() << std::endl; + if(vel.Z() > 0.0) + return spos.Z() < d1zmax_; + else + return spos.Z() > d0zmin_; + } + + size_t ExtrapolateCalo::nearestDisk(double zpos, double zvel) const { + size_t retval = disks_.size(); + if(zvel > 0.0){ // going forwards in z + for(auto idisk= disks_.begin(); idisk != disks_.end(); idisk++){ + auto const& diskptr = *idisk; + if(diskptr->center().Z() > zpos){ + retval = std::distance(disks_.begin(),idisk); + break; + } + } + } else { + for(auto idisk= disks_.rbegin(); idisk != disks_.rend(); idisk++){ + auto const& diskptr = *idisk; + if(diskptr->center().Z() < zpos){ + auto jdisk = idisk.base()-1; // points to the equivalent forwards object + retval = std::distance(disks_.begin(),jdisk); + break; + } + } + } + return retval; + } + +} +#endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index ab5a4205c7..b1b9f5fd8b 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -12,6 +12,7 @@ #include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateCalo.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/Mu2eKinKal/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" @@ -37,6 +38,8 @@ namespace mu2e { template bool extrapolateIPA(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; + template void toCaloD0(KKTrack& ktrk) const; + template void toCaloD1(KKTrack& ktrk) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; template void toOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; @@ -44,7 +47,19 @@ namespace mu2e { 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_; + // 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 bool geom_initialized_ = false; + void initGeometry() const; // lazy initialization method double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO }; @@ -55,26 +70,75 @@ 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()), + 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_); + + // 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 { + initGeometry(); // lazy initialization: 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) @@ -102,11 +166,12 @@ namespace mu2e { } template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { + 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 +186,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 +204,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 { + 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 +300,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 +327,13 @@ namespace mu2e { template bool KKExtrap::extrapolateST(KKTrack& ktrk,TimeDir tdir) const { using KKSTXING = KKShellXing; using KKSTXINGPTR = std::shared_ptr; + 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 +368,17 @@ namespace mu2e { } template bool KKExtrap::extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const { + 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; @@ -248,9 +387,10 @@ namespace mu2e { } template bool KKExtrap::extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const { + 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,19 +401,20 @@ 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 { + 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); } diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 3a881c64f6..36710ed780 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -114,6 +114,8 @@ 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 Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; }; diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index d45316dd45..0a90870f6d 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -7,6 +7,7 @@ #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" @@ -60,6 +61,7 @@ #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 "TTree.h" @@ -130,6 +132,7 @@ namespace mu2e { fhicl::Table extSettings { Name("ExtensionSettings") }; fhicl::Table matSettings { Name("MaterialSettings") }; // helix module specific config + fhicl::OptionalTable Extrapolation { Name("Extrapolation") }; }; class CentralHelixFit : public art::EDProducer { @@ -142,6 +145,7 @@ namespace mu2e { protected: bool goodFit(KKTRK const& ktrk) const; void sampleFit(KKTRK& kktrk) const; + void extrapolate(KKTRK& ktrk) const; TrkFitFlag fitflag_; // parameter-specific functions that need to be overridden in subclasses // data payload @@ -171,9 +175,11 @@ 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 + }; CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings}, fitflag_(TrkFitFlag::KKCentralHelix), @@ -195,7 +201,8 @@ namespace mu2e { useFitCharge_(settings().modSettings().useFitCharge()), minCenterRho_(settings().modSettings().minCenterRho()), sampleinrange_(settings().modSettings().sampleInRange()), - sampleinbounds_(settings().modSettings().sampleInBounds()) + sampleinbounds_(settings().modSettings().sampleInBounds()), + extrapolate_(false) { // collection handling for(const auto& cseedtag : settings().modSettings().seedCollections()) { cseedCols_.emplace_back(consumes(cseedtag)); } @@ -223,6 +230,16 @@ 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::beginRun(art::Run& run) { @@ -240,7 +257,7 @@ namespace mu2e { // translate the sample surface names to actual surfaces using the KinKalGeom. This must be done after construction as the KKGeom object now comes from GeometryService GeomHandle kkg_h; auto const& kkg = *kkg_h; - kkg.surfaces(ssids_,surfacess_to_sample_); + kkg.surfaces(ssids_,surfaces_to_sample_); } void CentralHelixFit::produce(art::Event& event ) { @@ -340,6 +357,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_); @@ -400,7 +419,7 @@ namespace mu2e { void CentralHelixFit::sampleFit(KKTRK& kktrk) const { auto const& ftraj = kktrk.fitTraj(); double tbeg = ftraj.range().begin(); - for(auto const& surf : surfacess_to_sample_){ + for(auto const& surf : surfaces_to_sample_){ // search for intersections with each surface from the begining double tstart = tbeg - sampletbuff_; bool goodinter(true); @@ -424,6 +443,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/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/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"