Reopening CaloSurfaceStep extraction with material effects (draft)

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<DetectorSystem> 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 "<<vdspmc.virtualDetectorId().name()<<std::endl;
       if(isid != vdmap_.end())ssc->emplace_back(isid->second,vdspmc,det); // no aggregation of VD hits
     }
     auto nvdsteps = ssc->size();

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_; }

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<SurfaceIdEnum::enum_type,std::string> const& SurfaceIdDetail::names(){
     return nam;

GeometryService/inc/KinKalGeomMaker.hh

@@ -13,6 +13,7 @@ namespace mu2e {
       void makeTracker();
       void makeDS();
       void makeTarget();
+      void makeCalo();
       void makeTCRV();
       std::unique_ptr<KinKalGeom> kkg_;
   };

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<mu2e::Calorimeter>());
+    auto const& tracker = *(GeomHandle<mu2e::Tracker>());
+
+    // 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<KKGeom::Calo>(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<Surface>(calo.EMC_Disk_0_OuterPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Inner),std::static_pointer_cast<Surface>(calo.EMC_Disk_0_InnerPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Inner),std::static_pointer_cast<Surface>(calo.EMC_Disk_1_InnerPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Outer),std::static_pointer_cast<Surface>(calo.EMC_Disk_1_OuterPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Front),std::static_pointer_cast<Surface>(calo.EMC_Disk_0_FrontPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Front),std::static_pointer_cast<Surface>(calo.EMC_Disk_1_FrontPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_0_Back),std::static_pointer_cast<Surface>(calo.EMC_Disk_0_BackPtr())));
+    kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::EMC_Disk_1_Back),std::static_pointer_cast<Surface>(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<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,zMidLocal),orvd,halfLen);
     auto inner = std::make_shared<Cylinder>(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<Disk>(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<Disk>(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<Disk>(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<KinKal::Disk>(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<KinKal::Disk>(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<KinKal::Disk>(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<Surface>(front)));
     kkg_->map_.emplace(std::make_pair(SurfaceId(SurfaceIdEnum::TT_Mid),std::static_pointer_cast<Surface>(mid)));
@@ -82,11 +138,11 @@ namespace mu2e {
     // currently use hard-coded geometry
     auto inner= std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),950,5450);
     auto outer= std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-1482),1328,5450); // bounding surfaces
-    auto front= std::make_shared<Disk>(outer->frontDisk());
-    auto back= std::make_shared<Disk>(outer->backDisk());
+    auto front= std::make_shared<KinKal::Disk>(outer->frontDisk());
+    auto back= std::make_shared<KinKal::Disk>(outer->backDisk());
     auto ipa= std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-2770),300.0,500.0);
-    auto ipafront= std::make_shared<Disk>(ipa->frontDisk());
-    auto ipaback= std::make_shared<Disk>(ipa->backDisk());
+    auto ipafront= std::make_shared<KinKal::Disk>(ipa->frontDisk());
+    auto ipaback= std::make_shared<KinKal::Disk>(ipa->backDisk());
     auto opa= std::make_shared<Frustrum>(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<Annulus>(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<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-4300),75,400.0);
     auto inner= std::make_shared<Cylinder>(VEC3(0.0,0.0,1.0),VEC3(0.0,0.0,-4300),21.5,400.0);
-    auto front= std::make_shared<Disk>(outer->frontDisk());
-    auto back= std::make_shared<Disk>(outer->backDisk());
+    auto front= std::make_shared<KinKal::Disk>(outer->frontDisk());
+    auto back= std::make_shared<KinKal::Disk>(outer->backDisk());
     double startz = -4700;
     double endz = -3900;
     double dz = (endz-startz)/36.0;

KinKalGeom/CMakeLists.txt

@@ -1,6 +1,7 @@
 cet_make_library(
     SOURCE
       src/CRV.cc
+      src/Calo.cc
       src/KinKalGeom.cc
     LIBRARIES PUBLIC
       KinKal::Geometry

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 <memory>
+namespace mu2e {
+  namespace KKGeom {
+    class Calo {
+      public:
+        using CylPtr = std::shared_ptr<KinKal::Cylinder>;
+        using DiskPtr = std::shared_ptr<KinKal::Disk>;
+        // 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

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<SurfacePair>& surfs) const;
       // find all surfaces that match a vector of Ids
       void surfaces(std::vector<SurfaceId> const& ids, std::vector<SurfacePair>& 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<KKGeom::Tracker> tracker_;
       std::unique_ptr<KKGeom::DetectorSolenoid> ds_;
       std::unique_ptr<KKGeom::StoppingTarget> st_;
+      std::unique_ptr<KKGeom::Calo> calo_;
       //KKGeom::CRV crv_;
       std::unique_ptr<KKGeom::TestCRV> tcrv_;
       // the map used to find surfaces by Id

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<Cylinder>(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<Cylinder>(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<Cylinder>(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<Cylinder>(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<Disk>(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<Disk>(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<Disk>(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<Disk>(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)
+    {}
+  }
+}

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