Fix segmentation position calculation for Allegro ECal.

detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h

@@ -3,6 +3,7 @@
 
 // FCCSW
 #include "detectorSegmentations/GridTheta_k4geo.h"
+#include <atomic>
 
 /** FCCSWGridModuleThetaMerged_k4geo Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged_k4geo.h FCCSWGridModuleThetaMerged_k4geo.h
  *
@@ -21,7 +22,7 @@ class FCCSWGridModuleThetaMerged_k4geo : public GridTheta_k4geo {
   FCCSWGridModuleThetaMerged_k4geo(const BitFieldCoder* decoder);
 
   /// destructor
-  virtual ~FCCSWGridModuleThetaMerged_k4geo() = default;
+  virtual ~FCCSWGridModuleThetaMerged_k4geo();
 
   /// read n(modules) from detector metadata
   void GetNModulesFromGeom();
@@ -33,15 +34,15 @@ class FCCSWGridModuleThetaMerged_k4geo : public GridTheta_k4geo {
    *   @param[in] aCellId ID of a cell.
    *   return Position (relative to R, phi of Geant4 volume it belongs to, scaled for R=1).
    */
-  virtual Vector3D position(const CellID& aCellID) const;
+  virtual Vector3D position(const CellID& aCellID) const override;
   /**  Determine the cell ID based on the position.
    *   @param[in] aLocalPosition (not used).
    *   @param[in] aGlobalPosition
    *   @param[in] aVolumeId ID of the Geant4 volume
    *   return Cell ID.
    */
   virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition,
-                        const VolumeID& aVolumeID) const;
+                        const VolumeID& aVolumeID) const override;
   /**  Determine the azimuthal angle (relative to the G4 volume) based on the cell ID.
    *   @param[in] aCellId ID of a cell.
    *   return Phi.
@@ -94,6 +95,20 @@ class FCCSWGridModuleThetaMerged_k4geo : public GridTheta_k4geo {
    */
   inline const std::string& fieldNameModule() const { return m_moduleID; }
 
+  /// Extract the layer index fom a cell ID.
+  int layer(const CellID& aCellID) const;
+
+  /// Determine the volume ID from the full cell ID by removing all local fields
+  virtual VolumeID volumeID(const CellID& cellID) const override;
+
+  /// Return true if this segmentation can have cells that span multiple
+  /// volumes.  That is, points from multiple distinct volumes may
+  /// be assigned to the same cell.
+  virtual bool cellsSpanVolumes() const override
+  {
+    return true;
+  }
+
 protected:
   /// the field name used for layer
   std::string m_layerID;
@@ -110,6 +125,24 @@ class FCCSWGridModuleThetaMerged_k4geo : public GridTheta_k4geo {
   /// number of layers (from the geometry)
   int m_nLayers;
 
+private:
+  /// Tabulate the cylindrical radii of all layers, as well as the
+  /// local x and z components needed for the proper phi offset.
+  struct LayerInfo {
+    LayerInfo(double the_rho, double the_xloc, double the_zloc)
+      : rho(the_rho), xloc(the_xloc), zloc(the_zloc)
+    {}
+    double rho;
+    double xloc;
+    double zloc;
+  };
+  std::vector<LayerInfo> initLayerInfo(const CellID& cID) const;
+
+  // The vector of tabulated values, indexed by layer number.
+  // We can't build this in the constructor --- the volumes won't have
+  // been created yet.  Instead, build it lazily the first time it's needed.
+  // Since that's in a const method, make it thread-safe.
+  mutable std::atomic<const std::vector<LayerInfo>*> m_layerInfo = nullptr;
 };
 }
 }

detectorSegmentations/src/FCCSWGridModuleThetaMerged_k4geo.cpp

@@ -2,6 +2,7 @@
 
 #include <iostream>
 #include "DD4hep/Detector.h"
+#include "DD4hep/VolumeManager.h"
 
 namespace dd4hep {
 namespace DDSegmentation {
@@ -35,6 +36,12 @@ FCCSWGridModuleThetaMerged_k4geo::FCCSWGridModuleThetaMerged_k4geo(const BitFiel
   GetNLayersFromGeom();
 }
 
+FCCSWGridModuleThetaMerged_k4geo::~FCCSWGridModuleThetaMerged_k4geo()
+{
+  delete m_layerInfo;
+}
+
+
 void FCCSWGridModuleThetaMerged_k4geo::GetNModulesFromGeom() {
   dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance());
   try {
@@ -59,14 +66,96 @@ void FCCSWGridModuleThetaMerged_k4geo::GetNLayersFromGeom() {
   std::cout << "Number of layers read from detector metadata and used in readout class: " << m_nLayers << std::endl;
 }
 
+/// Tabulate the cylindrical radii of all layers, as well as the
+/// local x and z components needed for the proper phi offset.
+std::vector<FCCSWGridModuleThetaMerged_k4geo::LayerInfo>
+FCCSWGridModuleThetaMerged_k4geo::initLayerInfo(const CellID& cID) const
+{
+  dd4hep::Detector* dd4hepgeo = &(dd4hep::Detector::getInstance());
+  VolumeManager vman = VolumeManager::getVolumeManager(*dd4hepgeo);
+
+  std::vector<LayerInfo> out;
+  out.reserve(m_nLayers);
+  VolumeID vID = cID;
+  _decoder->set(vID, m_thetaID, 0);
+  for (int l = 0; l < m_nLayers; l++) {
+
+    // Look up a volume in layer l in the volume manager, and find its radius
+    // by transforming the origin in the local coordinate system to global
+    // coordinates.
+    _decoder->set(vID, m_layerID, l);
+    VolumeManagerContext* vc = vman.lookupContext(vID);
+    Position wpos = vc->localToWorld({0,0,0});
+    double rho = wpos.Rho();
+
+    // If different modules are ganged together, we want to put hits
+    // in the center of all the ganged modules.  This represents
+    // a rotation in phi (in global coordinates).  Convert this rotation
+    // to displacement in local coordinates.  This will be in the x-z plane,
+    // and will be constant for all modules in a layer, but will be different
+    // for different layers (even with identical ganging).
+    double xloc = 0;
+    double zloc = 0;
+    double phioff = phi(vID);
+    if (phioff > 0) {
+      // We need to apply a phi offset.  Calculate it by rotating
+      // the global position in phi and converting back to local
+      // coordinates.
+      //
+      // For the Allegro calorimeter, this can also be calculated as:
+      //
+      //   d = rho * 2 * sin(phioff/2)
+      //   beta = pi/2 - phioff/2 - asin(sin(InclinationAngle)*EMBarrel_rmin/rho)
+      //   xloc = - d * sin(beta)
+      //   yloc =   d * cos(beta)
+      //
+      // but we prefer to do the calculation via explicit rotations because
+      // it's easier to see that that was is correct, and it also avoids
+      // the explicit dependencies on the geometry parameters.
+      Position wpos2 = RotateZ(wpos, phioff);
+      Position lpos2 = vc->worldToLocal(wpos2);
+      xloc = lpos2.X();
+      zloc = lpos2.Z();
+    }
+
+    out.emplace_back(rho, xloc, zloc);
+  }
+  return out;
+}
+
+
 /// determine the local position based on the cell ID
 Vector3D FCCSWGridModuleThetaMerged_k4geo::position(const CellID& cID) const {
 
+  // Get the vector of layer info.  If it hasn't been made yet,
+  // calculate it now.
+  const std::vector<LayerInfo>* liv = m_layerInfo.load();
+  if (!liv) {
+    auto liv_new = new std::vector<LayerInfo>(initLayerInfo(cID));
+    if (m_layerInfo.compare_exchange_strong(liv, liv_new)) {
+      liv = liv_new;
+    }
+    else {
+      delete liv_new;
+    }
+  }
+
+  VolumeID vID = cID;
+  _decoder->set(vID, m_thetaID, 0);
+  int layer = this->layer(vID);
+
   // debug
   // std::cout << "cellID: " << cID << std::endl;
 
-  // cannot return for R=0 otherwise will lose phi info, return for R=1
-  return positionFromRThetaPhi(1.0, theta(cID), phi(cID));
+  // Calculate the position in local coordinates.
+  // The volume here has the cross-section of a cell in the x-z plane;
+  // it extends the length of the calorimeter along the y-axis.
+  // We set the y-coordinate based on the theta bin, and x and z
+  // based on the phi offset required for this layer.
+  const LayerInfo& li = liv->at(layer);
+  return Vector3D(li.xloc,
+                  -li.rho / tan(theta(cID)),
+                  li.zloc);
 }
 
 /// determine the cell ID based on the global position
@@ -75,7 +164,7 @@ CellID FCCSWGridModuleThetaMerged_k4geo::cellID(const Vector3D& /* localPosition
   CellID cID = vID;
 
   // retrieve layer (since merging depends on layer)
-  int layer = _decoder->get(vID, m_layerID);
+  int layer = this->layer(vID);
 
   // retrieve theta
   double lTheta = thetaFromXYZ(globalPosition);
@@ -112,7 +201,7 @@ CellID FCCSWGridModuleThetaMerged_k4geo::cellID(const Vector3D& /* localPosition
 double FCCSWGridModuleThetaMerged_k4geo::phi(const CellID& cID) const {
 
   // retrieve layer
-  int layer = _decoder->get(cID, m_layerID);
+  int layer = this->layer(cID);
 
   // calculate phi offset due to merging
   // assume that m_mergedModules[layer]>=1
@@ -131,7 +220,7 @@ double FCCSWGridModuleThetaMerged_k4geo::phi(const CellID& cID) const {
 double FCCSWGridModuleThetaMerged_k4geo::theta(const CellID& cID) const {
 
   // retrieve layer
-  int layer = _decoder->get(cID, m_layerID);
+  int layer = this->layer(cID);
 
   // retrieve theta bin from cellID and determine theta position
   CellID thetaValue = _decoder->get(cID, m_thetaID);
@@ -152,5 +241,17 @@ double FCCSWGridModuleThetaMerged_k4geo::theta(const CellID& cID) const {
   return _theta;
 }
 
+/// Extract the layer index fom a cell ID.
+int FCCSWGridModuleThetaMerged_k4geo::layer(const CellID& cID) const {
+  return _decoder->get(cID, m_layerID);
+}
+
+/// Determine the volume ID from the full cell ID by removing all local fields
+VolumeID FCCSWGridModuleThetaMerged_k4geo::volumeID(const CellID& cID) const {
+  VolumeID vID = cID;
+  _decoder->set(vID, m_thetaID, 0);
+  return vID;
+}
+
 }
 }