Upgrade to WireTracker_info

CMakeLists.txt

@@ -27,6 +27,7 @@ find_package(ROOT COMPONENTS RIO Tree MathCore)
 find_package(podio 1.2 REQUIRED)
 find_package(EDM4HEP REQUIRED)
 find_package(k4FWCore 1.3 REQUIRED)
+find_package(k4geo REQUIRED)
 find_package(Gaudi REQUIRED)
 find_package(Delphes REQUIRED)
 find_package(GSL REQUIRED)

DCHdigi/CMakeLists.txt

@@ -60,6 +60,7 @@ gaudi_add_module(${PackageName}
   LINK
   k4FWCore::k4FWCore
   k4FWCore::k4Interface
+  k4geo::detectorCommon
   Gaudi::GaudiKernel
   EDM4HEP::edm4hep
   extensionDict

DCHdigi/include/DCHdigi_v01.h

@@ -55,13 +55,13 @@
 #include "DD4hep/Detector.h" // for dd4hep::VolumeManager
 #include "DDSegmentation/BitFieldCoder.h"
 
+// k4geo
+#include "detectorCommon/WireTracker_info.h"
+
 // STL
 #include <random>
 #include <string>
 
-// data extension for detector DCH_v2
-#include "DDRec/DCH_info.h"
-
 // ROOT headers
 #include "TFile.h"
 #include "TH1D.h"
@@ -71,6 +71,8 @@
 // Class developed by Walaa for the CLS
 #include "AlgData.h"
 
+using Vector3D = dd4hep::rec::WireTracker_info::Vector3D;
+
 /// constant to convert from mm (EDM4hep) to DD4hep (cm)
 
 struct DCHdigi_v01 final
@@ -142,10 +144,11 @@ struct DCHdigi_v01 final
 
   int CalculateNphiFromCellID(dd4hep::DDSegmentation::CellID id) const { return m_decoder->get(id, "nphi"); }
 
-  TVector3 Convert_EDM4hepVector_to_TVector3(const edm4hep::Vector3d& v, double scale) const {
+  /// Convert EDM4hep Vector3d to Vector3D as defined in WireTracker_info
+  Vector3D Convert_EDM4hepVector_to_Vector3D(const edm4hep::Vector3d& v, double scale) const {
     return {v[0] * scale, v[1] * scale, v[2] * scale};
   };
-  edm4hep::Vector3d Convert_TVector3_to_EDM4hepVector(const TVector3& v, double scale) const {
+  edm4hep::Vector3d Convert_Vector3D_to_EDM4hepVector(const  Vector3D& v, double scale) const {
     return {v.x() * scale, v.y() * scale, v.z() * scale};
   };
 

DCHdigi/include/DCHdigi_v02.h

@@ -64,6 +64,9 @@
 // k4FWCore
 #include "k4FWCore/Transformer.h"
 
+// k4geo
+#include "detectorCommon/WireTracker_info.h"
+
 // k4Interface
 #include <k4Interface/IGeoSvc.h>
 #include <k4Interface/IUniqueIDGenSvc.h>
@@ -77,15 +80,14 @@
 // DD4hep
 #include "DDSegmentation/BitFieldCoder.h"
 
-// DDRec
-#include "DDRec/DCH_info.h"
-
 // delphes
 #include "TrackCovariance/TrkUtil.h"
 
 // ROOT
 #include "TRandom3.h"
 
+using Vector3D = dd4hep::rec::WireTracker_info::Vector3D;
+
 class DCHdigi_v02 final
     : public k4FWCore::MultiTransformer<
           std::tuple<edm4hep::SenseWireHitCollection, edm4hep::TrackerHitSimTrackerHitLinkCollection>(
@@ -166,10 +168,10 @@ class DCHdigi_v02 final
       "Together with ReadoutWindowStartTime_ns, defines the readout window. Any DigiHits with arrival time after "
       "ReadoutWindowStartTime_ns + ReadoutWindowDuration_ns are discarded."};
 
-  /// Convert EDM4hep Vector3d to TVector3
-  TVector3 toTVector3(const edm4hep::Vector3d& v) const { return {v[0], v[1], v[2]}; };
-  /// Convert TVector3 to EDM4hep Vector3d
-  edm4hep::Vector3d toEDM4hepVector(const TVector3& v) const { return {v.x(), v.y(), v.z()}; };
+  /// Convert EDM4hep Vector3d to Vector3D as defined in WireTracker_info
+  Vector3D toVector3D(const edm4hep::Vector3d& v) const { return {v[0], v[1], v[2]}; };
+  /// Convert Vector3D as defined in WireTracker_info to EDM4hep Vector3d
+  edm4hep::Vector3d toEDM4hepVector(const Vector3D& v) const { return {v.x(), v.y(), v.z()}; };
 
   // /// Function to calculate the drift time from the distance to the wire
   double get_drift_time_ns(double distance_to_wire_mm) const;

DCHdigi/src/DCHdigi_v01.cpp

@@ -48,7 +48,8 @@ StatusCode DCHdigi_v01::initialize() {
   ///////////////////////////////////////////////////////////////////////////////////
   ///////////////////////////  retrieve data extension     //////////////////////////
   ///////////////////////////////////////////////////////////////////////////////////
-  this->dch_data = DCH_DE.extension<dd4hep::rec::DCH_info>();
+  auto wt_info =  DCH_DE.extension<dd4hep::rec::WireTracker_info_struct>();
+  this->dch_data = dynamic_cast<dd4hep::rec::DCH_info*>(wt_info);
   if (not dch_data->IsValid())
     ThrowException("No valid data extension was found for detector <<" + DCH_name + ">>.");
 
@@ -131,19 +132,20 @@ DCHdigi_v01::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits,
   // loop over hit collection
   for (const auto& input_sim_hit : input_sim_hits) {
     dd4hep::DDSegmentation::CellID cellid = input_sim_hit.getCellID();
+    int superlayer = this->CalculateLayerFromCellID(cellid);
     int ilayer = this->CalculateLayerFromCellID(cellid);
     int nphi = this->CalculateNphiFromCellID(cellid);
-    auto hit_position = Convert_EDM4hepVector_to_TVector3(input_sim_hit.getPosition(), MM_TO_CM);
+    auto hit_position = Convert_EDM4hepVector_to_Vector3D(input_sim_hit.getPosition(), MM_TO_CM);
 
     // -------------------------------------------------------------------------
     //      calculate hit position projection into the wire
-    TVector3 hit_to_wire_vector = this->dch_data->Calculate_hitpos_to_wire_vector(ilayer, nphi, hit_position);
-    TVector3 hit_projection_on_the_wire = hit_position + hit_to_wire_vector;
+    auto hit_to_wire_vector = this->dch_data->Calculate_hitpos_to_wire_vector(superlayer,ilayer, /*isector=*/0, nphi, hit_position);
+    auto hit_projection_on_the_wire = hit_position + hit_to_wire_vector;
     if (m_create_debug_histos.value()) {
-      double distance_hit_wire = hit_to_wire_vector.Mag();
+      double distance_hit_wire = hit_to_wire_vector.R();
       hDpw->Fill(distance_hit_wire);
     }
-    TVector3 wire_direction_ez = this->dch_data->Calculate_wire_vector_ez(ilayer, nphi);
+    auto wire_direction_ez = this->dch_data->Calculate_wire_vector_ez(superlayer, ilayer, /*isector=*/0, nphi);
 
     // -------------------------------------------------------------------------
     //       smear the position
@@ -156,15 +158,15 @@ DCHdigi_v01::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits,
     hit_projection_on_the_wire += smearing_z * (wire_direction_ez.Unit());
     if (m_create_debug_histos.value()) {
       // the distance from the hit projection and the wire should be zero
-      TVector3 dummy_vector = this->dch_data->Calculate_hitpos_to_wire_vector(ilayer, nphi, hit_projection_on_the_wire);
-      hDww->Fill(dummy_vector.Mag());
+      auto dummy_vector = this->dch_data->Calculate_hitpos_to_wire_vector(superlayer, ilayer, /*isector=*/0, nphi, hit_projection_on_the_wire);
+      hDww->Fill(dummy_vector.R());
     }
 
     //       smear position perpendicular to the wire
     double smearing_xy = gauss_xy_cm(rng_engine);
     if (m_create_debug_histos.value())
       hSxy->Fill(smearing_xy);
-    float distanceToWire_real = hit_to_wire_vector.Mag();
+    float distanceToWire_real = hit_to_wire_vector.R();
 
     // protect against negative values
     float distanceToWire_smeared = std::max(0.0, distanceToWire_real + smearing_xy);
@@ -173,16 +175,16 @@ DCHdigi_v01::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits,
     std::int32_t quality = 0;
     float eDepError = 0;
     // length units back to mm
-    auto positionSW = Convert_TVector3_to_EDM4hepVector(hit_projection_on_the_wire, 1. / MM_TO_CM);
-    // auto  directionSW    = Convert_TVector3_to_EDM4hepVector(wire_direction_ez, 1. / MM_TO_CM);
+    auto positionSW = Convert_Vector3D_to_EDM4hepVector(hit_projection_on_the_wire, 1. / MM_TO_CM);
+    // auto  directionSW    = Convert_XYZVector_to_EDM4hepVector(wire_direction_ez, 1. / MM_TO_CM);
     float distanceToWire = distanceToWire_smeared / MM_TO_CM;
 
     // The direction of the sense wires can be calculated as:
     //   RotationZ(WireAzimuthalAngle) * RotationX(stereoangle)
     // One point of the wire is for example the following:
     //   RotationZ(WireAzimuthalAngle) * Position(cell_rave_z0, 0 , 0)
     // variables aredefined below
-    auto WireAzimuthalAngle = this->dch_data->Get_cell_phi_angle(ilayer, nphi);
+    auto WireAzimuthalAngle = this->dch_data->Get_cell_phi_angle(superlayer, ilayer, /*isector=*/0, nphi);
     float WireStereoAngle = 0;
     {
       auto l = this->dch_data->database.at(ilayer);
@@ -191,7 +193,7 @@ DCHdigi_v01::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits,
       // when building the twisted tube, the twist angle is defined as:
       //     cell_twistangle    = l.StereoSign() * DCH_i->twist_angle
       // which forces the stereoangle of the wire to have the oposite sign
-      WireStereoAngle = (-1.) * l.StereoSign() * dch_data->stereoangle_z0(cell_rave_z0);
+      WireStereoAngle = (-1.) * dch_data->StereoSign(l) * dch_data->stereoangle_z0(cell_rave_z0);
     }
 
     extension::MutableSenseWireHit oDCHdigihit;

DCHdigi/src/DCHdigi_v02.cpp

@@ -10,6 +10,7 @@
 #include "DD4hep/Detector.h"
 
 // STL
+#include <detectorCommon/WireTracker_info.h>
 #include <ranges>
 
 namespace {
@@ -58,7 +59,8 @@ StatusCode DCHdigi_v02::initialize() {
   // Retrieve the detector element
   dd4hep::DetElement dch_detelem = m_geoSvc->getDetector()->detectors().at(dch_name);
   // Retrieve the DCH_info data extension for the drift chamber
-  m_dch_info = dch_detelem.extension<dd4hep::rec::DCH_info>();
+  auto wt_info =  dch_detelem.extension<dd4hep::rec::WireTracker_info_struct>();
+  m_dch_info = dynamic_cast<dd4hep::rec::DCH_info*>(wt_info);
   if (not m_dch_info->IsValid()) {
     error() << "No valid data extension was found for detector <<" << dch_name << ">>." << endmsg;
     return StatusCode::FAILURE;
@@ -118,8 +120,8 @@ DCHdigi_v02::operator()(const edm4hep::SimTrackerHitCollection& input,
   for (const auto& [cellID, simhits] : cell_map) {
 
     // Some geometry values needed for the calculations below
-    int layer = m_dch_info->CalculateILayerFromCellIDFields(m_decoder->get(cellID, "layer"),
-                                                            m_decoder->get(cellID, "superlayer"));
+    int superlayer = m_decoder->get(cellID, "superlayer"); 
+    int layer = m_dch_info->CalculateILayerFromCellIDFields(m_decoder->get(cellID, "layer"), superlayer);
     int nphi = m_decoder->get(cellID, "nphi");
 
     /* THE FOLLOWING CALCULATION OF WIRE ANGLES HAS BEEN COPIED AS IS FROM DCHdigi_v01! */
@@ -128,16 +130,14 @@ DCHdigi_v02::operator()(const edm4hep::SimTrackerHitCollection& input,
     // One point of the wire is for example the following:
     //   RotationZ(WireAzimuthalAngle) * Position(cell_rave_z0, 0 , 0)
     // variables are defined below
-    auto WireAzimuthalAngle = this->m_dch_info->Get_cell_phi_angle(layer, nphi);
+    auto WireAzimuthalAngle = this->m_dch_info->Get_cell_phi_angle(superlayer, layer, /*sector=*/0, nphi);
     float WireStereoAngle = 0;
     {
       auto l = this->m_dch_info->database.at(layer);
-      // radial middle point of the cell at Z=0
-      auto cell_rave_z0 = 0.5 * (l.radius_fdw_z0 + l.radius_fuw_z0);
       // when building the twisted tube, the twist angle is defined as:
-      //     cell_twistangle    = l.StereoSign() * DCH_i->twist_angle
+      //     cell_twistangle    = DCH_i->StereoSign(l) * DCH_i->twist_angle
       // which forces the stereoangle of the wire to have the oposite sign
-      WireStereoAngle = (-1.) * l.StereoSign() * m_dch_info->stereoangle_z0(cell_rave_z0);
+      WireStereoAngle = (-1.) * l.stereo_sw_z0;
     }
     /* END OF COPYING WIRE ANGLE CALCULATION FROM DCHdigi_v01 */
 
@@ -153,14 +153,14 @@ DCHdigi_v02::operator()(const edm4hep::SimTrackerHitCollection& input,
       //////////////////////////
 
       // Get hit position to calculate distance to wire
-      // Need to convert to TVector3 to use the DCH_info methods
+      // Need to convert to Vector3D to use the WireInfo_info methods
       // Use dd4hep:mm as scale to convert into the dd4hep default units (_ddu)
-      auto simhit_position_ddu = this->toTVector3(simhit.getPosition()) * dd4hep::mm;
+      auto simhit_position_ddu = this->toVector3D(simhit.getPosition()) * dd4hep::mm;
 
-      auto hit_to_wire_vector_ddu = m_dch_info->Calculate_hitpos_to_wire_vector(layer, nphi, simhit_position_ddu);
+      auto hit_to_wire_vector_ddu = m_dch_info->Calculate_hitpos_to_wire_vector(superlayer, layer, /*isector=*/0, nphi, simhit_position_ddu);
       auto hit_projection_on_the_wire_ddu = simhit_position_ddu + hit_to_wire_vector_ddu;
       double distance_to_wire_mm =
-          hit_to_wire_vector_ddu.Mag() / dd4hep::mm; // Explicitly cast to mm, no matter what the default unit is
+          hit_to_wire_vector_ddu.R() / dd4hep::mm; // Explicitly cast to mm, no matter what the default unit is
 
       ////////////////////////////////////////////
       // POSITION SMEARING AND TIME COMPUTATION //
@@ -172,7 +172,7 @@ DCHdigi_v02::operator()(const edm4hep::SimTrackerHitCollection& input,
 
       // z smearing
       double smearing_z_ddu = random_engine.Gaus(0.0, m_z_resolution_mm.value() * dd4hep::mm);
-      TVector3 wire_direction_ez_ddu = (m_dch_info->Calculate_wire_vector_ez(layer, nphi)).Unit();
+      auto wire_direction_ez_ddu = (m_dch_info->Calculate_wire_vector_ez(superlayer, layer,/*sector=*/0, nphi)).Unit();
       hit_projection_on_the_wire_ddu +=
           smearing_z_ddu *
           wire_direction_ez_ddu; // Need to multiply smearing_z_ddu with dd4hep::mm to cast into default units

Tracking/CMakeLists.txt

@@ -37,6 +37,7 @@ gaudi_add_module(${PackageName}
   EDM4HEP::edm4hep
   k4FWCore::k4FWCore
   k4FWCore::k4Interface
+  k4geo::detectorCommon
   extensionDict
   onnxruntime::onnxruntime
   torch

Tracking/components/GenfitTrackFitter.cpp

@@ -154,9 +154,9 @@ struct GenfitTrackFitter final
     // If the detector doesn't have a drift chamber, this part will be skipped
     try {
 
-      std::string DCH_name(m_DCH_name.value());
-      dd4hep::DetElement DCH_DE = m_geoSvc->getDetector()->detectors().at(DCH_name);
-      m_dch_info = DCH_DE.extension<dd4hep::rec::DCH_info>();
+      std::string DCH_name(m_WireTracker_name.value());
+      dd4hep::DetElement WireTracker_DE = m_geoSvc->getDetector()->detectors().at(DCH_name);
+      m_wire_info = WireTracker_DE.extension<dd4hep::rec::WireTracker_info_struct>();
       dd4hep::SensitiveDetector dch_sd = m_geoSvc->getDetector()->sensitiveDetector(DCH_name);
       dd4hep::Readout dch_readout = dch_sd.readout();
       m_dc_decoder = dch_readout.idSpec().decoder();
@@ -454,12 +454,12 @@ struct GenfitTrackFitter final
   GenfitInterface::GenfitMaterialInterface* m_geoMaterial;
 
   dd4hep::rec::SurfaceManager* m_surfMan;
-  dd4hep::rec::DCH_info* m_dch_info;
+  dd4hep::rec::WireTracker_info_struct* m_wire_info;
   dd4hep::DDSegmentation::BitFieldCoder* m_dc_decoder;
 
-  Gaudi::Property<std::string> m_DCH_name{
-      this, "DCHName", "DCH_v2",
-      "DCHName in the detector description (used to retrieve DCH geometry and material information)"};
+  Gaudi::Property<std::string> m_WireTracker_name{
+      this, "WireTrackerName", "DCH_v2",
+      "WireTrackerName in the detector description (used to retrieve Wire Tracker geometry and material information)"};
 
   // ====================== ECAL Geometry Parameters ======================
   double m_eCalBarrelInnerR;
@@ -616,7 +616,7 @@ struct GenfitTrackFitter final
                     edm4hep::TrackCollection& FittedTracks, edm4hep::TrackCollection& FittedTracksWithFilteredHits,
                     edm4hep::TrackerHitPlaneCollection& FittedHits, bool runCalorimeterExtrapolation) const {
 
-    GenfitInterface::GenfitTrack track_interface(track, m_skipTrackOrdering, m_dch_info, m_dc_decoder,
+    GenfitInterface::GenfitTrack track_interface(track, m_skipTrackOrdering, m_wire_info, m_dc_decoder,
                                                  m_genfitField.get());
 
     TVector3 Init_position(m_init_position.value()[0] * dd4hep::mm, m_init_position.value()[1] * dd4hep::mm,
@@ -773,7 +773,7 @@ struct GenfitTrackFitter final
 
     for (int pdgCode : m_particleHypothesis) {
 
-      GenfitInterface::GenfitTrack track_interface(track, m_skipTrackOrdering, m_dch_info, m_dc_decoder,
+      GenfitInterface::GenfitTrack track_interface(track, m_skipTrackOrdering, m_wire_info, m_dc_decoder,
                                                    m_genfitField.get());
 
       track_interface.InitializeTrack(m_RadialThresholdPromptTrack.value(), m_useFirstHitAsReference, LimitHits,

Tracking/include/genfit_interfaces/GenfitTrack.h

@@ -47,6 +47,8 @@
 #include "DDRec/Vector3D.h"
 #include "DDSegmentation/BitFieldCoder.h"
 
+#include "detectorCommon/WireTracker_info.h"
+
 #include "edm4hep/MutableTrack.h"
 #include "edm4hep/TrackCollection.h"
 
@@ -69,7 +71,7 @@
  *    - execute the GENFIT fitting procedure
  *
  *  The class maintains access to both the EDM4hep and GENFIT representations of the track, and stores
- *  relevant geometry information such as the drift chamber (DCH) description and segmentation decoder.
+ *  relevant geometry information such as the wire tracker info description and segmentation decoder.
  *
  *
  *  Author: Andrea De Vita
@@ -82,7 +84,7 @@ namespace GenfitInterface {
 class GenfitTrack {
 public:
   GenfitTrack(const edm4hep::Track& track, const bool skipTrackOrdering = false,
-              const dd4hep::rec::DCH_info* dch_info = nullptr,
+              const dd4hep::rec::WireTracker_info_struct* wire_info = nullptr,
               const dd4hep::DDSegmentation::BitFieldCoder* decoder = nullptr,
               const GenfitInterface::GenfitField* fieldMap = nullptr);
 
@@ -174,7 +176,7 @@ class GenfitTrack {
 
   TVector3 m_VP_referencePoint{0., 0., 0.};
 
-  const dd4hep::rec::DCH_info* m_dch_info;
+  const dd4hep::rec::WireTracker_info_struct* m_wire_info;
   const dd4hep::DDSegmentation::BitFieldCoder* m_dc_decoder;
   const GenfitInterface::GenfitField* m_fieldMap;
 };

Tracking/include/genfit_interfaces/GenfitWireMeasurement.h

@@ -30,9 +30,10 @@
 #include <TrackPoint.h>
 #include <WirePointMeasurement.h>
 
-#include "DDRec/DCH_info.h"
 #include "DDSegmentation/BitFieldCoder.h"
 
+#include "detectorCommon/WireTracker_info.h"
+
 #include "edm4hep/SenseWireHitCollection.h"
 
 /** @class WireMeasurement
@@ -56,9 +57,9 @@ namespace GenfitInterface {
 
 class WireMeasurement {
 public:
-  WireMeasurement(const edm4hep::SenseWireHit& hit, const dd4hep::rec::DCH_info* dch_info,
-                  const dd4hep::DDSegmentation::BitFieldCoder* decoder, const int det_idx, const int hit_idx,
-                  const int debug_lvl);
+  WireMeasurement(const edm4hep::SenseWireHit& hit, const dd4hep::rec::WireTracker_info_struct* wire_info,
+                  const dd4hep::DDSegmentation::BitFieldCoder* decoder, const bool has_sectors, const int det_idx,
+                  const int hit_idx, const int debug_lvl);
 
   genfit::WirePointMeasurement* getGenFit() const { return m_genfitHit; };