Skip to content

Commit

Permalink
Fix segmentation position calculation for Allegro ECal.
Browse files Browse the repository at this point in the history 
Fix the calculation of cell position for
FCCSWGridModuleThetaMerged_k4geo.
Also add a layer() method to that class to extract the layer from a cell ID.
  • Loading branch information
sss committed Jan 22, 2025
1 parent d7fec08 commit bd6c23a
Show file tree
Hide file tree
Showing 2 changed files with 142 additions and 8 deletions.
39 changes: 36 additions & 3 deletions detectorSegmentations/include/detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

// FCCSW
#include "detectorSegmentations/GridTheta_k4geo.h"
#include <atomic>

/** FCCSWGridModuleThetaMerged_k4geo Detector/DetSegmentation/DetSegmentation/FCCSWGridModuleThetaMerged_k4geo.h FCCSWGridModuleThetaMerged_k4geo.h
*
Expand All @@ -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();
Expand All @@ -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.
Expand Down Expand Up @@ -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;
Expand All @@ -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;
};
}
}
Expand Down
111 changes: 106 additions & 5 deletions detectorSegmentations/src/FCCSWGridModuleThetaMerged_k4geo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include <iostream>
#include "DD4hep/Detector.h"
#include "DD4hep/VolumeManager.h"

namespace dd4hep {
namespace DDSegmentation {
Expand Down Expand Up @@ -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 {
Expand All @@ -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
Expand All @@ -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);
Expand Down Expand Up @@ -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
Expand All @@ -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);
Expand All @@ -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;
}

}
}

0 comments on commit bd6c23a

Please sign in to comment.