CoreStripSpacePointFormationTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "CoreStripSpacePointFormationTool.h"
#include "Acts/Utilities/SpacePointUtility.hpp"
#include "InDetIdentifier/SCT_ID.h"
#include "ReadoutGeometryBase/SiCellId.h"
#include "InDetCondTools/ISiLorentzAngleTool.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "SCT_ReadoutGeometry/StripStereoAnnulusDesign.h"
#include "xAODInDetMeasurement/StripClusterAuxContainer.h"
#include "xAODInDetMeasurement/ContainerAccessor.h"
#include "Acts/SpacePointFormation/SpacePointBuilderConfig.hpp"
#include "ActsCalibrators/xAODUncalibMeasSurfAcc.h"
#include "ActsCalibrators/xAODUncalibMeasCalibrator.h"
#include "ActsGeometry/ATLASSourceLink.h"
 
#include "StoreGate/WriteHandle.h"
namespace ActsTrk
{
 
 
  StatusCode CoreStripSpacePointFormationTool::initialize(){
 
    ATH_CHECK(detStore()->retrieve(m_stripId, "SCT_ID"));
    ATH_CHECK(m_lorentzAngleTool.retrieve());
    ATH_CHECK(m_trackingGeometryTool.retrieve());
 
    if(m_useSCTLayerDep_OverlapCuts)
      ATH_MSG_INFO("Use SCT SP overlap cuts based on layer number parity");
    
    return StatusCode::SUCCESS;
  }
 
  StatusCode CoreStripSpacePointFormationTool::produceSpacePoints(const EventContext &ctx,
                  const xAOD::StripClusterContainer &clusterContainer,
                  const InDet::SiElementPropertiesTable &properties,
                  const InDetDD::SiDetectorElementCollection &elements,
                  const Amg::Vector3D &beamSpotVertex,
                  std::vector<StripSP>& spacePoints,
                  std::vector<StripSP>& overlapSpacePoints,
                  bool processOverlaps,
                  const std::vector<IdentifierHash>& hashesToProcess,
                  const ContainerAccessor<xAOD::StripCluster, IdentifierHash, 1>& stripAccessor) const
  {


 
    auto spBuilderConfig = std::make_shared<Acts::SpacePointBuilderConfig>();
 
    detail::xAODUncalibMeasSurfAcc surfaceAccessor{m_trackingGeometryTool.get()};
    spBuilderConfig->slSurfaceAccessor.connect<&detail::xAODUncalibMeasSurfAcc::operator()>(&surfaceAccessor);
    spBuilderConfig->trackingGeometry = m_trackingGeometryTool->trackingGeometry();
 
    
    auto spConstructor = [this, &clusterContainer, &elements](const Acts::Vector3 &pos,
                                                              std::optional<double> /*t*/,
                                                              const Acts::Vector2 &cov,
                                                              std::optional<double> /*varT*/,
                                                              const boost::container::static_vector<Acts::SourceLink, 2> &slinks)
      -> StripSP {
      std::array<std::size_t,2> measIndices{ std::numeric_limits<std::size_t>::max(), std::numeric_limits<std::size_t>::max()};
      std::array<StripInformationHelper, 2> stripInfos; 
      size_t idx = 0;
      for (const auto& slink : slinks){
          const xAOD::UncalibratedMeasurement *hit =detail::xAODUncalibMeasCalibrator::unpack(slink);
          // Check if the cluster is in the cluster container
          const auto it = std::ranges::find(clusterContainer, dynamic_cast<const xAOD::StripCluster*>(hit));
          if (it != clusterContainer.end()){
            const auto cluster_index = std::distance(clusterContainer.begin(), it);
            const auto id = hit->identifierHash();
            const auto &element = elements.getDetectorElement(id);
            size_t stripIndex = 0;
            auto ends = this->getStripEnds(hit, element, stripIndex);
            auto vertex = Amg::Vector3D::Zero();
            StripInformationHelper stripInfo(id,ends.first, ends.second, vertex, ActsTrk::localXFromSourceLink(hit), cluster_index, stripIndex);
            assert( idx < measIndices.size() && idx < stripInfos.size());
            measIndices[idx] = cluster_index;
            stripInfos[idx] = std::move(stripInfo);
            ++idx;
          }
      }
      const auto& [firstInfo, secondInfo] = stripInfos;
      const auto topHalfStripLength = 0.5*firstInfo.stripDirection().norm();
      Eigen::Matrix<double, 3, 1> topStripDirection = -firstInfo.stripDirection()/(2.*topHalfStripLength);
      Eigen::Matrix<double, 3, 1> topStripCenter = 0.5*firstInfo.trajDirection();
 
      const auto bottomHalfStripLength = 0.5*secondInfo.stripDirection().norm();
      Eigen::Matrix<double, 3, 1> bottomStripDirection = -secondInfo.stripDirection()/(2.*bottomHalfStripLength);
      Eigen::Matrix<double, 3, 1> stripCenterDistance = firstInfo.stripCenter()  - secondInfo.stripCenter();
      
      StripSP sp;
      sp.idHashes = {firstInfo.idHash(), secondInfo.idHash()};
      sp.globPos = pos.cast<float>();
      sp.cov_r = cov(0,0);
      sp.cov_z = cov(1,0);
      sp.measurementIndexes = measIndices;
      sp.topHalfStripLength = topHalfStripLength;
      sp.bottomHalfStripLength = bottomHalfStripLength;
      sp.topStripDirection = topStripDirection.cast<float>();
      sp.bottomStripDirection = bottomStripDirection.cast<float>();
      sp.stripCenterDistance = stripCenterDistance.cast<float>();
      sp.topStripCenter = topStripCenter.cast<float>();
 
      return sp;
    };
 
    auto spBuilder = std::make_shared<Acts::SpacePointBuilder<StripSP>>(*spBuilderConfig, spConstructor);
 
    const auto hashesProc = (hashesToProcess.size() > 0 ? hashesToProcess : stripAccessor.allIdentifiers());
    for (auto &idHash : hashesProc) {
      const InDetDD::SiDetectorElement *thisElement = elements.getDetectorElement(idHash);
      if (thisElement->isStereo()) {
        continue;
      }
      // Retrieve the neighbours of the detector element
      const std::vector<IdentifierHash>& others = *properties.neighbours(idHash);
  
      if (others.empty()) {
        continue;
      }
      // This flag is use to trigger if the search should be performed.
      // In case there are no clusters on the neighbours of the selected
      // detector element, the flag stays false.
      bool search = false;
      size_t neighbour = 0;
      while (not search and neighbour < others.size()){
        search = stripAccessor.isIdentifierPresent( others.at(neighbour) );
        neighbour++;
      }
      if (not search) {
        continue;
      }
      // prepare clusters, indices and modules for space point formation
      std::array<std::vector<std::pair<const xAOD::StripCluster *, size_t>>, static_cast<size_t>(nNeighbours)> neighbourClusters{};
      std::array<std::vector<std::pair<ATLASUncalibSourceLink, size_t>>, static_cast<size_t>(nNeighbours)> neighbourSourceLinks{};
      std::array<const InDetDD::SiDetectorElement *, static_cast<size_t>(nNeighbours)> neighbourElements{};
 
      auto groupStart = clusterContainer.begin();
      // Get the detector element and range for the idHash
      neighbourElements[0] = thisElement;
      for (auto &this_range : stripAccessor.rangesForIdentifierDirect(idHash)){
        for (auto start = this_range.first; start != this_range.second; ++start){
          size_t position = std::distance(groupStart, start);
          neighbourClusters[0].push_back(std::make_pair(*start, position));
                if ((*start)->identifierHash() != thisElement->identifyHash()) {
                  throw std::logic_error("Identifier mismatch.");
                }
          auto slink = makeATLASUncalibSourceLink(&clusterContainer, (*start)->index());
          neighbourSourceLinks[0].emplace_back(std::make_pair(slink, position));
        }
  }
 
  Identifier thisId = thisElement->identify();
 
  // define overlap extends before building space points
  std::array<double, 14> overlapExtents{};
  //   Default case: you test the opposite element and the overlapping in phi (total 3 elements)
  int Nmax = 4;
 
  // In the barrel, test the eta overlaps as well (total 5 elements)
  if (m_stripId->is_barrel(thisId))
    Nmax = 6;
 
  // You can remove all the overlaps if requested.
  // Here you test only the opposite element
  if (not processOverlaps)
    Nmax = 2;
 
  float hwidth(properties.halfWidth(idHash));
  int n = 0;
 
  // The order of the elements in others is such that you first get the opposite element,
  // the overlapping in phi and then the overlapping in eta
  // For this reason you need to re-order the indices, since the SiSpacePointMakerTool will process
  // first the eta overlaps and then the phi ones
  const std::array<size_t, nNeighbours> neigbourIndices{ThisOne, Opposite, EtaMinus, EtaPlus, PhiMinus, PhiPlus};
 
  for (const auto &otherHash : others){
    if (++n == Nmax) break;
 
    if (not stripAccessor.isIdentifierPresent(otherHash))
      continue;
 
    const InDetDD::SiDetectorElement *otherElement = elements.getDetectorElement(otherHash);
 
    neighbourElements[neigbourIndices[n]] = otherElement;
    for (auto &this_range : stripAccessor.rangesForIdentifierDirect(otherHash)){
      for (auto start = this_range.first; start != this_range.second; ++start){
        size_t position = std::distance(groupStart, start);
        neighbourClusters[neigbourIndices[n]].push_back(std::make_pair(*start, position));
              if ((*start)->identifierHash() != otherElement->identifyHash()) {
                 throw std::logic_error("Identifier mismatch.");
              }
              auto slink = makeATLASUncalibSourceLink(&clusterContainer, (*start)->index());
        neighbourSourceLinks[neigbourIndices[n]].emplace_back(std::make_pair(slink, position));
      }
    }
 
    switch (n){
    case Opposite:
      {
        overlapExtents[0] = -m_overlapLimitOpposite;
        overlapExtents[1] = m_overlapLimitOpposite;
        break;
      }
    case PhiMinus:
      {
        overlapExtents[6] = -hwidth;
        overlapExtents[7] = -hwidth + m_overlapLimitPhi;
        overlapExtents[8] = hwidth - m_overlapLimitPhi;
        overlapExtents[9] = hwidth;
        break;
      }
    case PhiPlus:
      {
        overlapExtents[10] = hwidth - m_overlapLimitPhi;
        overlapExtents[11] = hwidth;
        overlapExtents[12] = -hwidth;
        overlapExtents[13] = -hwidth + m_overlapLimitPhi;
        break;
      }
    case EtaMinus:
      {
        overlapExtents[ 2] = m_overlapLimitEtaMin;
        overlapExtents[ 3] = m_overlapLimitEtaMax;
        if (m_useSCTLayerDep_OverlapCuts && (m_stripId->layer_disk(thisId) & 1) != 0) {
    overlapExtents[2] = -m_overlapLimitEtaMax;
    overlapExtents[3] = -m_overlapLimitEtaMin;
        }
        break;
      }
    default:
      {
        overlapExtents[ 4] = m_overlapLimitEtaMin;
        overlapExtents[ 5] = m_overlapLimitEtaMax;
        if (m_useSCTLayerDep_OverlapCuts && (m_stripId->layer_disk(thisId) & 1) == 0) {
    overlapExtents[4] = -m_overlapLimitEtaMax;
    overlapExtents[5] = -m_overlapLimitEtaMin;
        }
        break;
      }
    }
  }
 
  ATH_CHECK( fillSpacePoints(ctx, spBuilder, neighbourElements, neighbourSourceLinks, overlapExtents, beamSpotVertex,
                                   spacePoints, overlapSpacePoints) );
      }
 
    return StatusCode::SUCCESS;
  }
 
  StatusCode CoreStripSpacePointFormationTool::fillSpacePoints(const EventContext &ctx,
                     std::shared_ptr<Acts::SpacePointBuilder<StripSP>> spBuilder,
                     const std::array<const InDetDD::SiDetectorElement *,nNeighbours>& elements,
                     const std::array<std::vector<std::pair<ATLASUncalibSourceLink, size_t>>,nNeighbours>& sourceLinks,
                     const std::array<double, 14>& overlapExtents,
                     const Amg::Vector3D &beamSpotVertex,
                     std::vector<StripSP>& spacePoints,
                     std::vector<StripSP>& overlapSpacePoints ) const
  {
    // This function is called once all the needed quantities are collected.
    // It is used to build space points checking the compatibility of clusters on pairs of detector elements.
    // Detector elements and cluster collections are elements and clusters, respectively.
    // [0] is the trigger element
    // [1] is the opposite element
    // [2]-[3] are the elements tested for eta overlaps
    // [4]-[5] are the elements tested for phi overlaps
    //
    // To build space points:
    // - For the opposite element and the ones tested for eta overlaps, you have to check
    //   if clusters are compatible with the local position of the trigger cluster
    //   requiring that the distance between the two clusters in phi is withing a specified range.
    //   - overlapExtents[0], overlapExtents[1] are filled for the opposite element
    //   - overlapExtents[2], overlapExtents[3], overlapExtents[4], overlapExtents[5] are filled for the eta overlapping elements
    // - For the elements tested for phi overlaps, you have to check
    //   if clusters are compatible with the local position of the trigger cluster.
    //   This needs that the trigger cluster is at the edge of the trigger module
    //   and that the other cluster is on the compatible edge of its module
    //   - overlapExtents[6], overlapExtents[7], overlapExtents[10], overlapExtents[11]
    //     overlapExtents[8], overlapExtents[9], overlapExtents[12], overlapExtents[13] are filled for the phi overlapping elements
 
 
 
    Acts::Vector3 vertex(beamSpotVertex.x(), beamSpotVertex.y(), beamSpotVertex.z());
    constexpr int otherSideIndex{1};
    constexpr int maxEtaIndex{3};
    std::array<int,nNeighbours - 1> elementIndex{};
    int nElements = 0;
 
    // For the nNeighbours sides, fill elementIndex with the indices of the existing elements.
    // Same the number of elements in nElements to loop on the later on
    for (int n = 1; n != nNeighbours; ++n) {
      if (elements[n]){
  elementIndex[nElements++] = n;
      }
    }
    // return if all detector elements are nullptr
    if(!nElements) return StatusCode::SUCCESS;
 
    const InDetDD::SiDetectorElement *triggerElement = elements[0];
 
    bool isEndcap = triggerElement->isEndcap();
    std::vector<StripInformationHelper> stripInfos;
    stripInfos.reserve(sourceLinks[0].size());
 
    std::vector<ATLASUncalibSourceLink> triggerSlinks;
    triggerSlinks.reserve(sourceLinks[0].size());
 
    // loop on all clusters on the trigger detector element and save the related information
    for (auto &sourceLink_index : sourceLinks[0]){
      triggerSlinks.emplace_back(sourceLink_index.first);
    }
 
    double limit = 1. + m_stripLengthTolerance;
    double slimit = 0.;
 
    if (not m_allClusters){
      // Start processing the opposite side and the eta overlapping elements
      int n = 0;
      for (; n < nElements; ++n){
  int currentIndex = elementIndex[n];
  if (currentIndex > maxEtaIndex)
    break;
 
  // get the detector element and the IdentifierHash
  const InDetDD::SiDetectorElement *currentElement = elements[currentIndex];
 
  // retrieve the range
  double min = overlapExtents[currentIndex * 2 - 2];
  double max = overlapExtents[currentIndex * 2 - 1];
 
  size_t minStrip, maxStrip = 0;
 
  if (m_stripGapParameter != 0.){
    updateRange(*triggerElement, *currentElement, slimit, min, max);
    correctPolarRange(triggerElement, min, max, minStrip, maxStrip);
  }
 
  StripInformationHelper currentStripInfo;
  for (auto &sourceLink_index : sourceLinks[currentIndex]){
          float source_local_x = ActsTrk::localXFromSourceLink( sourceLink_index.first );
    const auto currentSlink = sourceLink_index.first;
    for (auto triggerSlink : triggerSlinks){
 
      double diff = source_local_x - ActsTrk::localXFromSourceLink( triggerSlink );
      // In negative endcap, local z is opposite of positive endcap
      // need to invert the difference for proper comparison
      if( m_stripId->barrel_ec(currentElement->identify())<0 ) diff = -diff;
 
      if (diff < min || diff > max)
        continue;
      if (currentIndex == otherSideIndex){
        ATH_CHECK( makeSpacePoint(ctx, spacePoints, spBuilder, triggerSlink, currentSlink, triggerElement,
               currentElement, limit, slimit, vertex));
 
      } else {
        ATH_CHECK(makeSpacePoint(ctx, overlapSpacePoints, spBuilder, triggerSlink, currentSlink, triggerElement,
              currentElement, limit, slimit, vertex));
      }
    }
  }
      }
      // process the phi overlapping elements
      // if possible n starts from 4
      for (; n < nElements; ++n){
  int currentIndex = elementIndex[n];
  const InDetDD::SiDetectorElement *currentElement = elements[currentIndex];
  double min = overlapExtents[4 * currentIndex - 10];
  double max = overlapExtents[4 * currentIndex - 9];
 
  size_t minStrip, maxStrip = 0;
 
  if (m_stripGapParameter != 0.){
    updateRange(*triggerElement, *currentElement, slimit, min, max);
    correctPolarRange(triggerElement, min, max, minStrip, maxStrip);
  }
 
  std::vector<ATLASUncalibSourceLink> triggerPhiSlinks;
  triggerSlinks.reserve(triggerSlinks.size());
  for (auto triggerSlink : triggerSlinks){
    auto centralValue = ActsTrk::localXFromSourceLink( triggerSlink );
    auto minValue = min;
    auto maxValue = max;
    if (isEndcap){
      size_t stripIndex = 0;
      getStripEnds(triggerSlink, triggerElement, stripIndex);
      centralValue = stripIndex;
      minValue = minStrip;
      maxValue = maxStrip;
    }
    if (minValue <= centralValue and centralValue <= maxValue){
      triggerPhiSlinks.emplace_back(triggerSlink);
    }
  }
  if (triggerPhiSlinks.empty())
    continue;
  min = overlapExtents[4 * currentIndex - 8];
  max = overlapExtents[4 * currentIndex - 7];
  if (m_stripGapParameter != 0.){
    updateRange(*triggerElement, *currentElement, slimit, min, max);
    correctPolarRange(currentElement, min, max, minStrip, maxStrip);
  }
 
  for (auto &sourceLink_index : sourceLinks[currentIndex]){
    const auto currentSlink = sourceLink_index.first;
 
    size_t currentStripIndex = 0;
    getStripEnds(currentSlink, currentElement, currentStripIndex);
    auto centralValue = ActsTrk::localXFromSourceLink( sourceLink_index.first );
    auto minValue = min;
    auto maxValue = max;
    if (isEndcap) {
      centralValue = currentStripIndex;
      minValue = minStrip;
      maxValue = maxStrip;
    }
    if (centralValue < minValue or centralValue > maxValue)
      continue;
    for (auto &triggerSlink : triggerPhiSlinks) {
      ATH_CHECK(makeSpacePoint(ctx, overlapSpacePoints, spBuilder, triggerSlink, currentSlink, triggerElement,
            currentElement, limit, slimit, vertex));
    }
  }
      }
      return StatusCode::SUCCESS;
 
    } // not m_allClusters
 
    for (int n = 0; n != nElements; ++n){
 
      int currentIndex = elementIndex[n];
      const InDetDD::SiDetectorElement *currentElement = elements[currentIndex];
 
      if (m_stripGapParameter != 0.){
  computeOffset(*triggerElement, *currentElement, slimit);
      }
 
      for (auto &sourceLink_index : sourceLinks[currentIndex]){
  size_t currentStripIndex = 0;
  getStripEnds(sourceLink_index.first, triggerElement, currentStripIndex);
  const auto currentSlink = sourceLink_index.first;
    
  for (auto triggerSlink : triggerSlinks){
    if (currentIndex == otherSideIndex){
      ATH_CHECK(makeSpacePoint(ctx, spacePoints, spBuilder, triggerSlink, currentSlink, triggerElement,
            currentElement, limit, slimit, vertex));
    }else{
      ATH_CHECK(makeSpacePoint(ctx, overlapSpacePoints, spBuilder, triggerSlink, currentSlink, triggerElement,
            currentElement, limit, slimit, vertex));
    }
  }
      }
    }
    return StatusCode::SUCCESS;
  }
 
 
  StatusCode CoreStripSpacePointFormationTool::makeSpacePoint(const EventContext &ctx,
                    std::vector<StripSP>& collection,
                    std::shared_ptr<Acts::SpacePointBuilder<StripSP>> spBuilder,
                    const ATLASUncalibSourceLink& currentSlink,
                    const ATLASUncalibSourceLink& anotherSlink,
                    const InDetDD::SiDetectorElement *currentElement,
                    const InDetDD::SiDetectorElement *anotherElement,
                    const double limit,
                    const double slimit,
                    const Acts::Vector3& vertex) const
  {
 
    auto tgContext = m_trackingGeometryTool->getGeometryContext(ctx).context();
 
    size_t stripIndex = 0;
    auto ends1 = getStripEnds(currentSlink, currentElement, stripIndex);
    auto ends2 = getStripEnds(anotherSlink, anotherElement, stripIndex);
    std::pair<Acts::Vector3, Acts::Vector3> ends1_acts;
    std::pair<Acts::Vector3, Acts::Vector3> ends2_acts;
    ends1_acts.first = ends1.first;
    ends1_acts.second = ends1.second;
    ends2_acts.first = ends2.first;
    ends2_acts.second = ends2.second;
    auto paramCovAccessor = [&](const Acts::SourceLink &slink) {
      const auto &atlasSLink = slink.get<ATLASUncalibSourceLink>();
      const xAOD::UncalibratedMeasurement *measurement = &getUncalibratedMeasurement(atlasSLink);
      Acts::BoundVector loc = Acts::BoundVector::Zero();
      Acts::BoundSquareMatrix cov = Acts::BoundMatrix::Zero();
      switch (measurement->type()) {
      case (xAOD::UncalibMeasType::StripClusterType):
         loc[Acts::eBoundLoc0] = measurement->localPosition<1>()[Trk::locX];
         cov.topLeftCorner<1, 1>() =
            measurement->localCovariance<1>().cast<double>();
         break;
      case (xAOD::UncalibMeasType::PixelClusterType):
         loc[Acts::eBoundLoc0] = measurement->localPosition<2>()[Trk::locX];
         loc[Acts::eBoundLoc1] = measurement->localPosition<2>()[Trk::locY];
         cov.topLeftCorner<2, 2>() =
            measurement->localCovariance<2>().cast<double>();
         break;
      default:
         throw std::domain_error(
                                 "Can only handle measurement type pixel or strip");
      }
      return std::make_pair(loc, cov);
    };
    std::vector<Acts::SourceLink> slinks;
 
    slinks.emplace_back(Acts::SourceLink{currentSlink});
    slinks.emplace_back(Acts::SourceLink{anotherSlink});
 
    Acts::SpacePointBuilderOptions spOpt{std::make_pair(ends1_acts, ends2_acts), paramCovAccessor};
    spOpt.vertex = vertex;
    spOpt.stripLengthTolerance = limit - 1;
    spOpt.stripLengthGapTolerance = slimit;
 
 
    spBuilder->buildSpacePoint(tgContext, slinks, spOpt,
                               std::back_inserter(collection));
    return StatusCode::SUCCESS;
  }
 
double CoreStripSpacePointFormationTool::computeOffset(const InDetDD::SiDetectorElement& element1,
                   const InDetDD::SiDetectorElement& element2,
                   double& stripLengthGapTolerance) const
{
    // Get transformation matrices and center positions of detector elements
    const Amg::Transform3D& t1 = element1.transform();
    const Amg::Transform3D& t2 = element2.transform();
    const Amg::Vector3D& c = element1.center();
 
    // Check if first element is an annulus
    bool isAnnulus = (element1.design().shape() == InDetDD::Annulus);
 
    // Compute x12 and radius of detector element
    double x12 = t1.linear().col(0).dot(t2.linear().col(0));
    double r = isAnnulus ? c.perp() : std::sqrt(t1(0, 3) * t1(0, 3) + t1(1, 3) * t1(1, 3));
 
    // Compute distance between detector elements in the direction of strips
    Amg::Vector3D dPos = t1.translation() - t2.translation();
    double s = dPos.dot(t1.linear().col(2));
 
    // Compute offset distance
    double dm = (m_stripGapParameter * r) * std::abs(s * x12);
    double d = isAnnulus ? dm / 0.04 : dm / std::sqrt((1. - x12) * (1. + x12));
 
    // Adjust offset distance for z-component of transformation matrix
    const double zComponentTolerance = 0.7;
    if (std::abs(t1(2, 2)) > zComponentTolerance)
        d *= (r / std::abs(t1(2, 3)));
 
    stripLengthGapTolerance = d;
 
    return dm;
}
 
  void CoreStripSpacePointFormationTool::updateRange(const InDetDD::SiDetectorElement &element1,
                                                     const InDetDD::SiDetectorElement &element2,
                                                     double &stripLengthGapTolerance,
                                                     double &min, double &max) const
  {
    double dm = computeOffset(element1, element2, stripLengthGapTolerance);
    min -= dm;
    max += dm;
  }
 
  void CoreStripSpacePointFormationTool::correctPolarRange(const InDetDD::SiDetectorElement *element,
                                                           double &min,
                                                           double &max,
                                                           size_t &minStrip,
                                                           size_t &maxStrip) const
  {
    if (element->isBarrel())
      return;
 
    // design for endcap modules
    const InDetDD::StripStereoAnnulusDesign *design = dynamic_cast<const InDetDD::StripStereoAnnulusDesign *>(&element->design());
    if (!design){
      ATH_MSG_FATAL("Invalid strip annulus design for module with identifier/identifierHash " << element->identify() << "/" << element->identifyHash());
      return;
    }
 
    // converting min and max from cartesian reference frame to polar frame
    auto firstPosition = (design->localPositionOfCell(design->strip1Dim(0, 0)) +
                          design->localPositionOfCell(design->strip1Dim(design->diodesInRow(0) - 1, 0))) *
                         0.5;
 
    double radius = firstPosition.xEta();
 
    InDetDD::SiCellId minCellId = element->cellIdOfPosition(InDetDD::SiLocalPosition(radius, min, 0.));
    InDetDD::SiCellId maxCellId = element->cellIdOfPosition(InDetDD::SiLocalPosition(radius, max, 0.));
 
    if (not minCellId.isValid()) 
      minCellId = InDetDD::SiCellId(0);
 
    if (not maxCellId.isValid())
      maxCellId = InDetDD::SiCellId(design->diodesInRow(0) - 1);
 
    minStrip = minCellId.strip();
    maxStrip = maxCellId.strip();
 
    // re-evaluate min and max in polar coordinate
    min = std::atan2(min, radius);
    max = std::atan2(max, radius);
  }
 
  std::pair<Amg::Vector3D, Amg::Vector3D>
  CoreStripSpacePointFormationTool::getStripEnds(const ATLASUncalibSourceLink &sourceLink,
             const InDetDD::SiDetectorElement *element,
             size_t &stripIndex) const
  {
    const xAOD::UncalibratedMeasurement &measurement = getUncalibratedMeasurement(sourceLink);
    auto cluster = dynamic_cast<const xAOD::StripCluster *>(&measurement);
    if(!cluster){
      ATH_MSG_FATAL("Could not cast UncalibratedMeasurement as StripCluster");
      return {};
    }
    float source_local_x = ActsTrk::localXFromSourceLink( sourceLink );
    if (element->isEndcap())  {
      // design for endcap modules
      const InDetDD::StripStereoAnnulusDesign *design = dynamic_cast<const InDetDD::StripStereoAnnulusDesign *>(&element->design());
      if (!design){
        ATH_MSG_FATAL("Invalid strip annulus design for module with identifier/identifierHash " << element->identify() << "/" << element->identifyHash());
        return {};
      }
 
      // calculate phi pitch for evaluating the strip index
      double phiPitchPhi = design->phiWidth()/design->diodesInRow(0);
      stripIndex = -std::floor(source_local_x / phiPitchPhi) + design->diodesInRow(0) *0.5 - 0.5;
 
      std::pair<Amg::Vector3D, Amg::Vector3D > ends = {
        element->globalPosition(design->stripPosAtR(stripIndex, 0, design->minR())),
        element->globalPosition(design->stripPosAtR(stripIndex, 0, design->maxR()))
      };
      return ends;
    }
 
    InDetDD::SiLocalPosition localPosition(0., source_local_x, 0.);
    std::pair<Amg::Vector3D, Amg::Vector3D> ends(element->endsOfStrip(localPosition));
    return ends;
  }
 
}