AnalogueClusteringToolImpl.icc

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
#ifndef ANALOGUECLUSTERINGTOOLIMPL_ICC
#define ANALOGUECLUSTERINGTOOLIMPL_ICC
 
#include <limits>
#include "ActsGeometry/ActsDetectorElement.h"
#include <PixelReadoutGeometry/PixelModuleDesign.h>
#include "AthenaBaseComps/AthMsgStreamMacros.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "InDetMeasurementUtilities/Helpers.h"
 
#include <stdexcept>
 
namespace ActsTrk::detail {
 
template <typename calib_data_t, typename traj_t>
StatusCode AnalogueClusteringToolImpl<calib_data_t, traj_t>::initialize()
{
    ATH_MSG_DEBUG("Initializing " << AthAlgTool::name() << " ...");
 
    ATH_CHECK( BASE::initialize() );
    ATH_CHECK(m_clusterErrorKey.initialize()); 
    ATH_MSG_DEBUG(AthAlgTool::name() << " successfully initialized");
    
    ATH_MSG_DEBUG( this->calibrateAfterMeasurementSelection() );
    ATH_MSG_DEBUG( m_correctCovariance );
    ATH_MSG_DEBUG( m_calibratedCovarianceLowerBound );
    ATH_MSG_DEBUG( m_errorStrategy );
 
    return StatusCode::SUCCESS;
}
 
template <typename calib_data_t, typename traj_t>
std::pair<float, float>
AnalogueClusteringCalibrator<calib_data_t, traj_t>::getCentroid(
    const EventContext& ctx,
    const xAOD::PixelCluster& cluster,
    const InDetDD::SiDetectorElement& element) const
{
  // Reproduce the base position of the legacy ITk::PixelClusterOnTrackTool,
  // i.e. the geometric 4-corner bounding-box centroid plus the Lorentz shift. This
  // is the reference the getDelta charge-interpolation constants are calibrated
  // against, so the analogue correction must be applied on top of it.
  SG::ConstAccessor<SG::JaggedVecElt<Identifier::value_type> >::element_type
     rdos = cluster.rdoList();
    
  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();
  // convert Identifier::value_type into Identifier
  for (Identifier::value_type rid_value : rdos) {
    Identifier rid(rid_value);
    int row = this->pixelID().phi_index(rid);
    rowmin = std::min(row, rowmin);
    rowmax = std::max(row, rowmax);
    
    int col = this->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 = this->getLorentzShift(element.identifyHash(), ctx);
  
  return std::make_pair(centroid.xPhi() + shift, centroid.xEta());
}
 
template <typename calib_data_t, typename traj_t>
const typename AnalogueClusteringCalibrator<calib_data_t,traj_t>::error_data_t*
AnalogueClusteringToolImpl<calib_data_t, traj_t>::getErrorData(const EventContext &ctx) const
{
    SG::ReadCondHandle<calib_data_t> handle(m_clusterErrorKey,ctx);
 
    if (!handle.isValid()) {
      ATH_MSG_ERROR(m_clusterErrorKey << " is not available.");
      return nullptr;
    }
 
    const typename AnalogueClusteringCalibrator<calib_data_t,traj_t>::error_data_t* data = handle->getClusterErrorData();
    if (data == nullptr) {
      ATH_MSG_ERROR("No cluster error data corresponding to " << m_clusterErrorKey);
      return nullptr;
    }
 
    return data;
}
 
namespace {
   // negative limit of eta when computing eta as eta(T x) {return -log(tan(x*0.5)); } with x=M_PI-epsilon
   // was eta<T>(static_cast<T>(M_PI)-std::numeric_limits<T>::epsilon());
   // with T eta(T theta) { return static_cast<float>(-1*log(tan(theta*0.5))); }, and T=float
   constexpr float etaMin() { return -16.3991603f; }
   // positive limit of eta when computing eta as eta(T x) {return -log(tan(x*0.5)); } with x=0+epsilon
   // was eta<T>(static_cast<T>(0)-std::numeric_limits<T>::epsilon()), with T=float
   constexpr float etaMax() { return  16.6355323f; }
 
   // compute eta from tan_theta
   double computeEta(double tan_theta) {
      // eta = - ln ( tan theta/2)
      // tan theta /2 = tan theta / (1 + sec theta)
      // 1+tan^2 theta = sec^2 theta
      // sec theta = sqrt ( 1+ tan^2 theta)
      double one_plus_sec_theta = 1.+std::copysign(std::sqrt ( 1. +  std::pow(tan_theta,2) ), tan_theta );
 
      // reproduce roughly the original range;
      // @TODO what would be the allowed range for getDelta ?
      static constexpr double min_eta = etaMin();
      static constexpr double max_eta = etaMax();
 
      return std::max(min_eta,
                      std::min( max_eta,
                               (std::abs(one_plus_sec_theta)>tan_theta*std::numeric_limits<double>::epsilon()
                                ? -std::log( tan_theta / ( one_plus_sec_theta ))
                                :min_eta ) ));
   }
}
 
 
template <typename calib_data_t, typename traj_t>
std::pair<std::optional<float>, std::optional<float>>
AnalogueClusteringCalibrator<calib_data_t, traj_t>::getCorrectedPosition(
    const EventContext& ctx,
    const xAOD::PixelCluster& cluster,
    const error_data_t& errorData,
    const InDetDD::SiDetectorElement& element,
    const std::pair<float, float>& tan_angles) const
{
    auto& [tan_anglePhi, tan_angleEta] = tan_angles; // Phi and Eta denote the direction of the axis
 
    std::optional<float> posX {std::nullopt};
    std::optional<float> posY {std::nullopt};
 
    int nrows = cluster.channelsInPhi();
    int ncols = cluster.channelsInEta();
 
    // If size of cluster (in any direction) is higher than 1, we may need to apply a correction
    if (nrows > 1 or ncols > 1) {
      const auto& [omX, omY] = TrackingUtilities::computeOmegas(cluster,
                                                                this->pixelID());
      // if the omegas are valid, we can apply the correction
      if (omX > -0.5 and omY > -0.5) {
        std::pair<float, float> centroid = getCentroid(ctx, cluster, element);
        IdentifierHash idHash = element.identifyHash();
 
        // getDelta expects the pseudo-rapidity built from the polar angle
        // measured with respect to the module *plane* (theta = pi/2 - angleFromNormal),
        // exactly as the legacy ITk::PixelClusterOnTrackTool does. tan_angleEta is
        // the tangent of the incidence angle with respect to the module *normal*, so
        // the polar-angle tangent is its cotangent: tan(pi/2 - b) = 1/tan(b).
        // (For all quadrants of boweta the legacy thetaloc satisfies tan(thetaloc) = cot(boweta).)
        const double local_eta = tan_angleEta != 0.f ? computeEta(1.0 / tan_angleEta) : 0.0;
        std::pair<double, double> delta =
          errorData.getDelta(idHash,
                             nrows,
                             atan(tan_anglePhi),
                             ncols,
                             local_eta);
        
        if (nrows > 1)
          posX = centroid.first + delta.first * (omX - 0.5);
        
        if (ncols > 1)
          posY = centroid.second + delta.second * (omY - 0.5);
      } // check on omega values
    } // check on cluster sizes
    
    return std::make_pair(posX, posY);
}
 
  template <typename calib_data_t, typename traj_t>
  std::pair<std::optional<float>, std::optional<float>>
  AnalogueClusteringCalibrator<calib_data_t, traj_t>::getCorrectedError(
    const typename AnalogueClusteringCalibrator<calib_data_t, traj_t>::error_data_t& errorData,
    [[maybe_unused]] const InDetDD::SiDetectorElement& element,
    const std::pair<float, float>& tan_angles,
    const xAOD::PixelCluster& cluster) const
  {
    std::optional<float> errX {std::nullopt};
    std::optional<float> errY {std::nullopt};
    if (not m_options.m_correctCovariance) {
      return std::make_pair(errX, errY);
    }
    
    typename BASE::Cov cov = cluster.template localCovariance<2>();
    int nrows = cluster.channelsInPhi();
    int ncols = cluster.channelsInEta();
 
    auto& [tan_anglePhi, tan_angleEta] = tan_angles;
 
    if (m_options.m_errorStrategy == 0) {
    
      // Special case for very shallow tracks
      // Error estimated from geometrical projection of
      // the track path in silicon onto the module surface
      static constexpr float tan_1 = 1.55740772f; // std::tan(1.f);
      if (std::abs(tan_anglePhi) > tan_1) {
        static const double one_over_sqrt12 = 1./std::sqrt(12);
        double thickness = element.design().thickness();
        errX = thickness * std::abs(tan_anglePhi) * one_over_sqrt12;
        errY = thickness * std::abs(tan_angleEta);
        if (cluster.widthInEta() > errY) {
          errY = cluster.widthInEta()  * one_over_sqrt12;
        } else {
          errY = *errY * one_over_sqrt12;
        }
      } else if (nrows > 1 or ncols > 1) {
        IdentifierHash idHash = cluster.identifierHash();
        auto [vX, vY] = errorData.getDeltaError(idHash);
        if (nrows > 1 and vX > 0)
          errX = vX;
        if (ncols > 1 and vY > 0)
          errY = vY;
      }
    } else if (m_options.m_errorStrategy == 1) { // this should check if we are using broad errors in the clustering tool
 
      //This assumes that we ran broad errors in the clustering tool
      errX = std::sqrt(cov(0,0)) / nrows;
      errY = std::sqrt(cov(1,1)) / ncols;
      
    } else  { // throw exception for the moment
      throw std::runtime_error("Only errorStrategy 0 or 1 is supported for the moment");
    }
 
    return std::make_pair(errX, errY);
}
 
template <typename calib_data_t, typename traj_t>
std::tuple<typename AnalogueClusteringCalibrator<calib_data_t, traj_t>::BASE::Pos,
           typename AnalogueClusteringCalibrator<calib_data_t, traj_t>::BASE::Cov,
           unsigned int>
AnalogueClusteringCalibrator<calib_data_t, traj_t>::calibrate(
    const EventContext& ctx,
    const Acts::GeometryContext& /*gctx*/,
    const Acts::CalibrationContext& /*cctx*/,
    const xAOD::PixelCluster& cluster,
    const InDetDD::SiDetectorElement& detElement,
    const std::pair<float, float>& tan_incidence_angles) const
{  
    typename BASE::Pos pos = cluster.template localPosition<2>();
    typename BASE::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(ctx, cluster, *errorData, detElement, tan_incidence_angles);
    if (posX.has_value())
      pos[Acts::eBoundLoc0] = *posX;
    if (posY.has_value())
      pos[Acts::eBoundLoc1] = *posY;
    
    auto [errX, errY] = getCorrectedError(*errorData, detElement, tan_incidence_angles, cluster);
    if (errX.has_value()) {
      double newErr = *errX;
      cov(0, 0) = std::max(newErr * newErr, cov(0, 0) * m_options.m_calibratedCovarianceLowerBound);
    }
    if (errY.has_value()) {
      double newErr = *errY;
      cov(1, 1) = std::max(newErr * newErr, cov(1, 1) * m_options.m_calibratedCovarianceLowerBound);
    }
    
    return std::make_tuple(pos, cov, 0u);
}
 
 
} // namespace ActsTrk
 
#endif