AnalogueClusteringToolImpl.icc

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef ANALOGUECLUSTERINGTOOLIMPL_ICC
#define ANALOGUECLUSTERINGTOOLIMPL_ICC
 
#include <limits>
 
#include <PixelReadoutGeometry/PixelModuleDesign.h>
#include "InDetReadoutGeometry/SiDetectorElement.h"
 
namespace ActsTrk::detail {
 
template <typename calib_data_t, typename traj_t>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::AnalogueClusteringToolImpl(
    const std::string& type,
    const std::string& name,
    const IInterface* parent)
    : base_class(type, name, parent)
{}
 
template <typename calib_data_t, typename traj_t>
StatusCode AnalogueClusteringToolImpl<calib_data_t, traj_t>::initialize()
{
    ATH_MSG_DEBUG("Initializing " << AthAlgTool::name() << " ...");
    ATH_CHECK(m_pixelDetEleCollKey.initialize());
    ATH_CHECK(m_clusterErrorKey.initialize()); 
    ATH_CHECK(m_lorentzAngleTool.retrieve());
    ATH_CHECK(AthAlgTool::detStore()->retrieve(m_pixelid, "PixelID"));
    ATH_MSG_DEBUG(AthAlgTool::name() << " successfully initialized");
 
    ATH_MSG_DEBUG( m_thickness );
    ATH_MSG_DEBUG( m_postCalibration );
    ATH_MSG_DEBUG( m_correctCovariance );
    ATH_MSG_DEBUG( m_calibratedCovarianceLowerBound );
 
    return StatusCode::SUCCESS;
}
 
template <typename calib_data_t, typename traj_t>
const InDetDD::SiDetectorElement&
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getDetectorElement(xAOD::DetectorIDHashType id) const
{
    const InDetDD::SiDetectorElement* element=nullptr;
    SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
    m_pixelDetEleCollKey,
    Gaudi::Hive::currentContext());
 
    const InDetDD::SiDetectorElementCollection* detElements(*pixelDetEleHandle);
    
    if (!pixelDetEleHandle.isValid() or detElements == nullptr) {
    ATH_MSG_ERROR(m_pixelDetEleCollKey.fullKey() << " is not available.");
    }
    else {
         element = detElements->getDetectorElement(id);
    }
    if (element == nullptr) {
    ATH_MSG_ERROR("No element corresponding to hash " << id << " for " << m_pixelDetEleCollKey.fullKey());
    throw std::runtime_error("SiDetectorElement is NULL");
    }
 
    return *element;
}
 
template <typename calib_data_t, typename traj_t>
std::pair<float, float>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::anglesOfIncidence(const InDetDD::SiDetectorElement& element,
                                                                    const Acts::Vector3& direction) const
{
 
    float projPhi = direction.dot(element.phiAxis());
    float projEta = direction.dot(element.etaAxis());
    float projNorm = direction.dot(element.normal());
 
    float anglePhi = std::atan2(projPhi, projNorm);
    float angleEta = std::atan2(projEta, projNorm);
 
    // Map the angles of inward-going tracks onto [-PI/2, PI/2]
    if (anglePhi > M_PI *0.5) {
    anglePhi -= M_PI;
    }
    if (anglePhi < -M_PI *0.5) {
    anglePhi += M_PI;
    }
 
    // settle the sign/pi periodicity issues
    float thetaloc;
    if (angleEta > -0.5 * M_PI && angleEta < M_PI / 2.) {
    thetaloc = M_PI_2 - angleEta;
    } else if (angleEta > M_PI_2 && angleEta < M_PI) {
    thetaloc = 1.5 * M_PI - angleEta;
    } else { // 3rd quadrant
    thetaloc = -M_PI_2 - angleEta;
    }
    angleEta = -1 * log(tan(thetaloc * 0.5));
 
    
    // Subtract the Lorentz angle effect
    float angleShift = m_lorentzAngleTool->getTanLorentzAngle(element.identifyHash(), Gaudi::Hive::currentContext());
    anglePhi = std::atan(std::tan(anglePhi) - element.design().readoutSide() * angleShift);
 
    return std::make_pair(anglePhi, angleEta);
}
 
 
template <typename calib_data_t, typename traj_t>
std::pair<float, float>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCentroid(
    const std::vector<Identifier>& rdos,
    const InDetDD::SiDetectorElement& element) const
{
 
    int rowmin = std::numeric_limits<int>::max();
    int rowmax = std::numeric_limits<int>::min();
    int colmin = std::numeric_limits<int>::max();
    int colmax = std::numeric_limits<int>::min();
    for (const Identifier& rid : rdos) {
        int row = m_pixelid->phi_index(rid);
        rowmin = std::min(row, rowmin);
        rowmax = std::max(row, rowmax);
 
        int col = m_pixelid->eta_index(rid);
        colmin = std::min(col, colmin);
        colmax = std::max(col, colmax);
    }
 
    const InDetDD::PixelModuleDesign& design =
    dynamic_cast<const InDetDD::PixelModuleDesign&>(element.design());
 
    InDetDD::SiLocalPosition pos1 =
    design.positionFromColumnRow(colmin, rowmin);
 
    InDetDD::SiLocalPosition pos2 =
    design.positionFromColumnRow(colmax, rowmin);
 
    InDetDD::SiLocalPosition pos3 =
    design.positionFromColumnRow(colmin, rowmax);
 
    InDetDD::SiLocalPosition pos4 =
    design.positionFromColumnRow(colmax, rowmax);
 
    InDetDD::SiLocalPosition centroid = 0.25 * (pos1 + pos2 + pos3 + pos4);
 
    double shift = m_lorentzAngleTool->getLorentzShift(element.identifyHash(), Gaudi::Hive::currentContext());
 
    return std::make_pair(centroid.xPhi() + shift, centroid.xEta());
}
 
template <typename calib_data_t, typename traj_t>
const typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::error_data_t*
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getErrorData() const
{
    SG::ReadCondHandle<calib_data_t> handle(
    m_clusterErrorKey,
    Gaudi::Hive::currentContext());
 
    if (!handle.isValid()) {
    ATH_MSG_ERROR(m_clusterErrorKey << " is not available.");
    return nullptr;
    }
 
    const error_data_t* data = handle->getClusterErrorData();
    if (data == nullptr) {
    ATH_MSG_ERROR("No cluster error data corresponding to " << m_clusterErrorKey);
    return nullptr;
    }
 
    return data;
}
 
template <typename calib_data_t, typename traj_t>
std::pair<std::optional<float>, std::optional<float>>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedPosition(
    const std::vector<Identifier>& rdos,
    const error_data_t& errorData,
    const InDetDD::SiDetectorElement& element,
    const std::pair<float, float>& angles,
    const xAOD::PixelCluster& cluster) const
{
    // TODO validate these angles
    auto& [anglePhi, angleEta] = angles;
 
    std::optional<float> posX;
    std::optional<float> posY;
    
    float omX = cluster.omegaX();
    float omY = cluster.omegaY();
    int nrows = cluster.channelsInPhi();
    int ncols = cluster.channelsInEta();
 
    if (omX > -0.5 && omY > -0.5 && (nrows > 1 || ncols > 1)) {
 
    std::pair<float, float> centroid = getCentroid(rdos, element);
 
        IdentifierHash idHash = element.identifyHash();
        
    std::pair<double, double> delta =
        errorData.getDelta(
        idHash,
        nrows,
        anglePhi,
        ncols,
        angleEta);
 
    if (nrows > 1)
        posX = centroid.first + delta.first * (omX - 0.5);
 
    if (ncols > 1)
        posY = centroid.second + delta.second * (omY - 0.5);
    }
 
    return std::make_pair(posX, posY);
}
 
template <typename calib_data_t, typename traj_t>
std::pair<std::optional<float>, std::optional<float>>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedError(
    const error_data_t& errorData,
    const InDetDD::SiDetectorElement& element,
    const std::pair<float, float>& angles,
    const xAOD::PixelCluster& cluster) const
{
  std::optional<float> errX {std::nullopt};
  std::optional<float> errY {std::nullopt};
  if (not m_correctCovariance) {
    return std::make_pair(errX, errY);
  }
  
    int nrows = cluster.channelsInPhi();
    int ncols = cluster.channelsInEta();
 
    auto& [anglePhi, angleEta] = angles;
 
    // Special case for very shallow tracks
    // Error estimated from geometrical projection of
    // the track path in silicon onto the module surface
    if (std::abs(anglePhi) > 1) {
    errX = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(anglePhi)) / std::sqrt(12);
    errY = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(angleEta));
    if (cluster.widthInEta() > errY) {
        errY = cluster.widthInEta() / std::sqrt(12);
    } else {
        errY = *errY / std::sqrt(12);
    }
    } else if (nrows > 1 or ncols > 1) {
        IdentifierHash idHash = element.identifyHash();
    auto [vX, vY] = errorData.getDeltaError(idHash);
    if (vX > 0)
      errX = vX;
    if (vY > 0)
      errY = vY;
    }
 
    return std::make_pair(errX, errY);
}
 
template <typename calib_data_t, typename traj_t>
std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
      typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
    const Acts::GeometryContext& gctx,
    const Acts::CalibrationContext& cctx,
    const xAOD::PixelCluster& cluster,
    const TrackStateProxy& state) const
{
    const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
    Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
    state.parameters()[Acts::eBoundPhi],
    state.parameters()[Acts::eBoundTheta]);
    return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
}
 
template <typename calib_data_t, typename traj_t>
std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
      typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
    const Acts::GeometryContext& gctx,
    const Acts::CalibrationContext& cctx,
    const xAOD::PixelCluster& cluster,
    const Acts::BoundTrackParameters& bound_parameters) const
{
    const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
    Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
    bound_parameters.parameters()[Acts::eBoundPhi],
    bound_parameters.parameters()[Acts::eBoundTheta]);
    return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
}
 
template <typename calib_data_t, typename traj_t>
std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
      typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
    const Acts::GeometryContext& /*gctx*/,
    const Acts::CalibrationContext& /*cctx*/,
    const xAOD::PixelCluster& cluster,
    const InDetDD::SiDetectorElement& detElement,
    const std::pair<float, float>& angles) const
{
    Pos pos = cluster.template localPosition<2>();
    Cov cov = cluster.template localCovariance<2>();
 
    assert(!cluster.rdoList().empty());
 
    const error_data_t *errorData = getErrorData();
    if (errorData == nullptr) {
    throw std::runtime_error("PixelClusterErrorData is NULL");
    }
 
 
    auto [posX, posY] = getCorrectedPosition(cluster.rdoList(), *errorData, detElement, angles, cluster);
    if (posX.has_value())
    pos[Acts::eBoundLoc0] = *posX;
    if (posY.has_value())
    pos[Acts::eBoundLoc1] = *posY;
 
    auto [errX, errY] = getCorrectedError(*errorData, detElement, angles, cluster);
    if (errX.has_value()) {
      double newErr = *errX;
      cov(0, 0) = std::max(newErr * newErr, cov(0, 0) * m_calibratedCovarianceLowerBound);
    }
    if (errY.has_value()) {
      double newErr = *errY;
      cov(1, 1) = std::max(newErr * newErr, cov(1, 1) * m_calibratedCovarianceLowerBound);
    }
    
    return std::make_pair(pos, cov);
}
 
template <typename calib_data_t, typename traj_t>
void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connect(OnTrackCalibrator<traj_t>& calibrator) const
{
   using CalibFuncPtr_t =
      std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
                typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
      (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
                                                             const Acts::CalibrationContext&,
                                                             const xAOD::PixelCluster&,
                                                             const TrackStateProxy&) const;
 
   calibrator.pixelCalibrator. template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
}
 
template <typename calib_data_t, typename traj_t>
void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connectPixelCalibrator(IOnBoundStateCalibratorTool::PixelCalibrator& pixelCalibrator) const
{
   using CalibFuncPtr_t =
      std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
                typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
      (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
                                                             const Acts::CalibrationContext&,
                                                             const xAOD::PixelCluster&,
                                                             const Acts::BoundTrackParameters&) const;
 
   pixelCalibrator.template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
}
 
template <typename calib_data_t, typename traj_t>
bool AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrateAfterMeasurementSelection() const
{ return m_postCalibration.value(); }
 
} // namespace ActsTrk
 
#endif