SiDetElementsRoadMaker_xk.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/

//   Implementation file for class InDet::SiDetElementsRoadMaker_xk
// (c) ATLAS Detector software
// Version 1.0 21/04/2004 I.Gavrilenko
 
#include "SiDetElementsRoadTool_xk/SiDetElementsRoadMaker_xk.h"
#include "SiSPSeededTrackFinderData/SiDetElementRoadMakerData_xk.h"
 
#include "SiDetElementsRoadUtils_xk.h"
 
#include "GaudiKernel/ContextSpecificPtr.h"
#include "PixelReadoutGeometry/PixelDetectorManager.h"
#include "SCT_ReadoutGeometry/SCT_DetectorManager.h"
#include "SiDetElementsRoadTool_xk/SiDetElementsComparison.h"
#include "StoreGate/ReadCondHandle.h"
#include "TrkPrepRawData/PrepRawData.h"
#include <cmath>
 
#include <ostream>
#include <iomanip>

// Constructor
 
InDet::SiDetElementsRoadMaker_xk::SiDetElementsRoadMaker_xk
(const std::string& t, const std::string& n, const IInterface* p)
  : base_class(t, n, p)
{
}

// Initialisation
 
StatusCode InDet::SiDetElementsRoadMaker_xk::initialize()
{
  //Class optimization checks
  static_assert(std::is_trivially_copyable<SiDetElementLink_xk::ElementWay>::value);
  static_assert(std::is_trivially_destructible<SiDetElementLink_xk::ElementWay>::value);
  static_assert(std::is_trivially_copyable<SiDetElementLink_xk>::value);
  static_assert(std::is_trivially_destructible<SiDetElementLink_xk>::value);
  static_assert(std::is_nothrow_move_constructible<SiDetElementsLayer_xk>::value);
 
 
  if (!m_usePIX && !m_useSCT) {
    ATH_MSG_FATAL("Please don't call this tool if usePixel and useSCT are false");
    return StatusCode::FAILURE;
  }
 
  if (m_fieldmode == "NoField") m_fieldModeEnum = Trk::NoField;
  else if (m_fieldmode == "MapSolenoid") m_fieldModeEnum = Trk::FastField;
  else m_fieldModeEnum = Trk::FullField;
 
  // Get propagator tool
  //
  ATH_CHECK(m_proptool.retrieve());
 
  // Get output print level
  //
  m_outputlevel = msg().level()-MSG::DEBUG;
 
  computeBounds();
 
  ATH_CHECK(m_layerVecKey.initialize());
  ATH_CHECK(m_fieldCondObjInputKey.initialize());
 
  return StatusCode::SUCCESS;
}

// Finalize
 
StatusCode InDet::SiDetElementsRoadMaker_xk::finalize()
{
  return StatusCode::SUCCESS;
}
 

// Dumps relevant information into the MsgStream
 
MsgStream& InDet::SiDetElementsRoadMaker_xk::dump(MsgStream& out) const
{
  out<<"\n";
  return dumpConditions(out);
}

// Dumps conditions information into the MsgStream
 
MsgStream& InDet::SiDetElementsRoadMaker_xk::dumpConditions(MsgStream& out) const
{
  int n = 62-m_proptool.type().size();
  std::string s1;
  for (int i=0; i<n; ++i) s1.append(" ");
  s1.append("|");
 
  std::string fieldmode[9] = {"NoField"       , "ConstantField", "SolenoidalField",
                              "ToroidalField" , "Grid3DField"  , "RealisticField" ,
                              "UndefinedField", "AthenaField"  , "?????"         };
 
  Trk::MagneticFieldMode fieldModeEnum(m_fieldModeEnum);
 
  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldHandle(m_fieldCondObjInputKey);
  const AtlasFieldCacheCondObj* fieldCondObj(*fieldHandle);
  if (fieldCondObj) {
    MagField::AtlasFieldCache fieldCache;
    fieldCondObj->getInitializedCache(fieldCache);
    if (!fieldCache.solenoidOn()) fieldModeEnum = Trk::NoField;
  }
 
  Trk::MagneticFieldProperties fieldprop(fieldModeEnum);
 
  int mode = fieldprop.magneticFieldMode();
  if (mode<0 || mode>8) mode = 8;
 
  n = 62-fieldmode[mode].size();
  std::string s3;
  for (int i=0; i<n; ++i) s3.append(" ");
  s3.append("|");
 
  n = 62-m_sct.size();
  std::string s5;
  for (int i=0; i<n; ++i) s5.append(" ");
  s5.append("|");
 
  n = 62-m_pix.size();
  std::string s6;
  for (int i=0; i<n; ++i) s6.append(" ");
  s6.append("|");
 
  const EventContext& ctx = Gaudi::Hive::currentContext();
 
  const SiDetElementsLayerVectors_xk &layer = *getLayers(ctx);
 
  int maps = 0;
  if (!layer[0].empty()) ++maps;
  if (!layer[1].empty()) ++maps;
  if (!layer[2].empty()) ++maps;
  auto prec = out.precision();
  out<<"|----------------------------------------------------------------------"
     <<"-------------------|"
     <<"\n";
  if (m_useSCT) {
    out<<"| SCT   detector manager  | "<<m_sct             <<s5<<"\n";
  }
  if (m_usePIX) {
    out<<"| Pixel detector manager  | "<<m_pix             <<s6<<"\n";
  }
  out<<"| Tool for propagation    | "<<m_proptool.type() <<s1<<"\n";
  out<<"| Magnetic field mode     | "<<fieldmode[mode]   <<s3<<"\n";
  out<<"| Width of the road (mm)  | "
     <<std::setw(12)<<std::setprecision(5)<<m_width
     <<"                                                  |"<<"\n";
  out<<"|----------------------------------------------------------------------"
     <<"-------------------|"
     <<"\n";
 
  if (!maps || m_outputlevel==0) return out;
 
  if (!layer[1].empty()) {
    int nl = layer[1].size();
    int nc = 0;
    for (const auto & i : layer[1]) nc+=i.nElements();
    out<<"|----------------------------------------------------------------|"
       <<"\n";
    out<<"| Barrel map contains "
       <<std::setw(3)<<nl<<" layers and"
       <<std::setw(5)<<nc<<" elements               |"
       <<"\n";
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    out<<"|   n  |     R     |   Z min    |   Z max    |  max dF    | nEl  |"
       <<"\n";
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    for (unsigned int i=0; i!=layer[1].size(); ++i) {
      double zmin = layer[1].at(i).z()-layer[1].at(i).dz();
      double zmax = layer[1].at(i).z()+layer[1].at(i).dz();
      out<<"| "
     <<std::setw(4)<<i<<" |"
     <<std::setw(10)<<std::setprecision(4)<<  layer[1].at(i).r ()<<" | "
     <<std::setw(10)<<std::setprecision(4)<<                zmin<<" | "
     <<std::setw(10)<<std::setprecision(4)<<                zmax<<" | "
     <<std::setw(10)<<std::setprecision(4)<< layer[1].at(i).dfe()<<" | "
     <<std::setw(4)<<layer[1].at(i).nElements()<<" | "
     <<"\n";
    }
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
 
  }
  if (!layer[0].empty()) {
    int nl = layer[0].size();
    int nc = 0;
    for (const auto & i : layer[0]) nc+=i.nElements();
    out<<"|----------------------------------------------------------------|"
       <<"\n";
    out<<"| L.Endcap map contains"
       <<std::setw(3)<<nl<<" layers and"
       <<std::setw(5)<<nc<<" elements              |"
       <<"\n";
 
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    out<<"|   n  |     Z     |   R min    |   R max    |  max dF    | nEl  |"
       <<"\n";
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    for (unsigned int i=0; i!=layer[0].size(); ++i) {
      double rmin = layer[0].at(i).r()-layer[0].at(i).dr();
      double rmax = layer[0].at(i).r()+layer[0].at(i).dr();
      out<<"| "
     <<std::setw(4)<<i<<" |"
     <<std::setw(10)<<std::setprecision(4)<<  layer[0].at(i).z()<<" | "
     <<std::setw(10)<<std::setprecision(4)<<               rmin<<" | "
     <<std::setw(10)<<std::setprecision(4)<<               rmax<<" | "
     <<std::setw(10)<<std::setprecision(4)<<layer[0].at(i).dfe()<<" | "
     <<std::setw(4)<<layer[0].at(i).nElements()<<" | "
     <<"\n";
    }
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
  }
  if (!layer[2].empty()) {
    int nl = layer[2].size();
    int nc = 0;
    for (const auto & i : layer[2]) nc+=i.nElements();
    out<<"|----------------------------------------------------------------|"
       <<"\n";
    out<<"| R.Endcap map contains"
       <<std::setw(3)<<nl<<" layers and"
       <<std::setw(5)<<nc<<" elements              |"
       <<"\n";
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    out<<"|   n  |     Z     |   R min    |   R max    |  max dF    | nEl  |"
       <<"\n";
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
    for (unsigned int i=0; i!=layer[2].size(); ++i) {
      double rmin = layer[2].at(i).r()-layer[2].at(i).dr();
      double rmax = layer[2].at(i).r()+layer[2].at(i).dr();
      out<<"| "
     <<std::setw(4)<<i<<" |"
     <<std::setw(10)<<std::setprecision(4)<<  layer[2].at(i).z()<<" | "
     <<std::setw(10)<<std::setprecision(4)<<               rmin<<" | "
     <<std::setw(10)<<std::setprecision(4)<<               rmax<<" | "
     <<std::setw(10)<<std::setprecision(4)<<layer[2].at(i).dfe()<<" | "
     <<std::setw(4)<<layer[2].at(i).nElements()<<" | "
     <<"\n";
    }
    out<<"|------|-----------|------------|------------|------------|------|"
       <<"\n";
  }
  out<<"\n";
  out.precision(prec);
  return out;
}

// Dumps relevant information into the ostream
 
std::ostream& InDet::SiDetElementsRoadMaker_xk::dump(std::ostream& out) const
{
  return out;
}

// Overload of << operator MsgStream
 
MsgStream& InDet::operator <<
(MsgStream& sl, const InDet::SiDetElementsRoadMaker_xk& se)
{
  return se.dump(sl);
}

// Overload of << operator std::ostream
 
std::ostream& InDet::operator <<
(std::ostream& sl, const InDet::SiDetElementsRoadMaker_xk& se)
{
  return se.dump(sl);
}
 

// Main methods for road builder using input list global positions
 
void InDet::SiDetElementsRoadMaker_xk::detElementsRoad
(std::deque<Amg::Vector3D>& globalPositions,
 std::vector<const InDetDD::SiDetectorElement*>& Road,
 bool testDirection,
 SiDetElementRoadMakerData_xk & roadMakerData,
 const EventContext& ctx) const
{
  if (!m_usePIX && !m_useSCT) return;

  const SiDetElementsLayerVectors_xk &layer = *getLayers(ctx);

  std::deque<Amg::Vector3D>::iterator currentPosition=globalPositions.begin(), endPositions=globalPositions.end();

  std::array<float,6> par_startingPoint{static_cast<float>((*currentPosition).x()),            // x of first position
                static_cast<float>((*currentPosition).y()),             // y of first position
                static_cast<float>((*currentPosition).z()),             // Z of first position
                static_cast<float>(sqrt((*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y())),    // r of first position
                m_width,    // road width
                0.};

  int n0 = 0;
  for (; n0!=static_cast<int>(layer[0].size()); ++n0) {
    if (par_startingPoint[2] > layer[0][n0].z()) break;
  }

  int n1 = 0;
  for (; n1!=static_cast<int>(layer[1].size()); ++n1) {
    if (par_startingPoint[3] < layer[1][n1].r()) break;
  }
  int n2 = 0;
  for (; n2!=static_cast<int>(layer[2].size()); ++n2) {
    if (par_startingPoint[2] < layer[2][n2].z()) break;
  }
 
 

  if (!roadMakerData.isInitialized){
    bookUsageTracker(roadMakerData,layer);
  }
  else{
    roadMakerData.resetUsageTracker();
  }
 
  std::vector<InDet::SiDetElementLink_xk::ElementWay> lDE;
  lDE.reserve(8); //reasonable minimum guess
  ++currentPosition;
  while (currentPosition!=endPositions) {
    std::array<float,4>  par_targetPoint{static_cast<float>((*currentPosition).x()),         
                    static_cast<float>((*currentPosition).y()),         
                    static_cast<float>((*currentPosition).z()),         
                        static_cast<float>(sqrt((*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y()))   
                    };

    float dx = par_targetPoint[0]-par_startingPoint[0];         
    float dy = par_targetPoint[1]-par_startingPoint[1];         
    float dz = par_targetPoint[2]-par_startingPoint[2];         
    float dist3D = std::sqrt(dx*dx+dy*dy+dz*dz);   
    if (dist3D <=0.) {                    
      ++currentPosition;                  
      continue;
    }
    float inverseDistance = 1./dist3D;                   
    std::array<float,3> searchDirection{dx*inverseDistance, dy*inverseDistance, dz*inverseDistance};

    float unitSepTransverseComp = searchDirection[0]*searchDirection[0]+searchDirection[1]*searchDirection[1];     
    float dr = 0.                 ;
    if (unitSepTransverseComp!=0.) {
      float sm = -( searchDirection[0]*par_startingPoint[0] +
                    searchDirection[1]*par_startingPoint[1])
                    /unitSepTransverseComp;

      if (sm > 1. && sm < dist3D) {   
        par_targetPoint[0] = par_startingPoint[0]+searchDirection[0]*sm;
        par_targetPoint[1] = par_startingPoint[1]+searchDirection[1]*sm;
        par_targetPoint[2] = par_startingPoint[2]+searchDirection[2]*sm;
        par_targetPoint[3] = std::sqrt(par_targetPoint[0]*par_targetPoint[0]+par_targetPoint[1]*par_targetPoint[1]);
          dr    = 20.;  
      } else {
        ++currentPosition;
      }
    } else {
      ++currentPosition;
    }



    if (par_targetPoint[3]>par_startingPoint[3]) {
      for (; n1<static_cast<int>(layer[1].size()); ++n1) {
          if (par_targetPoint[3] < layer[1][n1].r()) break;
        assert( roadMakerData.elementUsageTracker[1].size() > static_cast<unsigned int>(n1) );
        if(m_ITkGeometry) layer[1][n1].getITkBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
        else layer[1][n1].getBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
      }
    } else {
      for (--n1; n1>=0; --n1) {
          if (par_targetPoint[3] > layer[1][n1].r()+dr) break;
        assert( roadMakerData.elementUsageTracker[1].size() > static_cast<unsigned int>(n1) );
        if(m_ITkGeometry) layer[1][n1].getITkBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
        else layer[1][n1].getBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
 
      }
      ++n1;
    }

    if (par_targetPoint[2]>par_startingPoint[2]) {
      for (; n2<static_cast<int>(layer[2].size()); ++n2) {
          if (par_targetPoint[2] < layer[2][n2].z()) break;
        assert( roadMakerData.elementUsageTracker[2].size() > static_cast<unsigned int>(n2) );
        if(m_ITkGeometry) layer[2][n2].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[2][n2]);
        else layer[2][n2].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[2][n2]);
      }
    } else {
      for (--n2; n2>=0; --n2) {
          if (par_targetPoint[2] > layer[2][n2].z()) break;
        assert( roadMakerData.elementUsageTracker[2].size() > static_cast<unsigned int>(n2) );
        if(m_ITkGeometry) layer[2][n2].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[2][n2]);
        else layer[2][n2].getEndcapDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[2][n2]);
      }
      ++n2;
    }

    if (par_targetPoint[2]<par_startingPoint[2]) {
      for (; n0<static_cast<int>(layer[0].size()); ++n0) {
          if (par_targetPoint[2] > layer[0][n0].z()) break;
        assert( roadMakerData.elementUsageTracker[0].size() > static_cast<unsigned int>(n0) );
        if(m_ITkGeometry) layer[0][n0].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
        else layer[0][n0].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
      }
    } else {
      for (--n0; n0>=0; --n0) {
          if (par_targetPoint[2] < layer[0][n0].z()) break;
        assert( roadMakerData.elementUsageTracker[0].size() > static_cast<unsigned int>(n0) );
        if(m_ITkGeometry) layer[0][n0].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
        else layer[0][n0].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
      }
      ++n0;
    }
    par_startingPoint[0] = par_targetPoint[0];
    par_startingPoint[1] = par_targetPoint[1];
    par_startingPoint[2] = par_targetPoint[2];
    par_startingPoint[3] = par_targetPoint[3];
    par_startingPoint[5]+= dist3D;
  }
  auto vec2 = lDE;
  std::sort(lDE.begin(),lDE.end(),InDet::compDetElementWays());
  // Fill pointers to detector elements
  Road.reserve(lDE.size());
  for (auto & d : lDE){
    if (testDirection && d.way() < 0) {continue;}
    Road.push_back(d.link()->detElement());
  }
}

void InDet::SiDetElementsRoadMaker_xk::bookUsageTracker(InDet::SiDetElementRoadMakerData_xk & data, const SiDetElementsLayerVectors_xk &layers) {

    for ( unsigned int side_i=0; side_i<3; ++side_i) {
      data.elementUsageTracker[side_i].resize( layers[side_i].size() ); 
      for (unsigned int layer_i=0; layer_i < layers[side_i].size(); ++layer_i) {
          data.elementUsageTracker[side_i][layer_i].resize( layers[side_i][layer_i].nElements() );
      }
    }
    data.isInitialized=true;
}

// Main methods for road builder using track parameters and direction
 
void InDet::SiDetElementsRoadMaker_xk::detElementsRoad
(const EventContext& ctx,
 MagField::AtlasFieldCache& fieldCache,
 const Trk::TrackParameters& Tp,
 Trk::PropDirection direction,
 std::vector<const InDetDD::SiDetectorElement*>& Road,
 SiDetElementRoadMakerData_xk & roadMakerData) const
{
  if (!m_usePIX && !m_useSCT) return;
  double qp   = fabs(500.*Tp.parameters()[4]);
  if (qp < 1.e-10) qp = 1.e-10;
  double S    = m_step/qp;
  if (S  > 1000. ) S  = 1000. ;
 
  bool testDirection = true;
  if (direction<0) {
    testDirection = false;
    S=-S;
  }
 
  Trk::MagneticFieldMode fieldModeEnum(m_fieldModeEnum);
  if (!fieldCache.solenoidOn()) fieldModeEnum = Trk::NoField;
  Trk::MagneticFieldProperties fieldprop(fieldModeEnum);
 
  // Note: could also give fieldCache directly to propagator if it would be more efficient - would
  // need to add interface RDS 2020/03
  std::deque<Amg::Vector3D> G;

  m_proptool->globalPositions(ctx, G, Tp, fieldprop,getBound(fieldCache, Tp), S, Trk::pion);
  if (G.size()<2) return;

  if (direction > 0) {
    std::deque<Amg::Vector3D>::iterator currentPosition=G.begin(), nextPosition, endPositions=G.end();
    float r0 = (*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y();
 
    while (currentPosition!=endPositions) {
      nextPosition = currentPosition;
      if (++nextPosition == endPositions) break;
 
      float r = (*nextPosition).x()*(*nextPosition).x()+(*nextPosition).y()*(*nextPosition).y();
      if (r < r0) {
        r0 = r;
        currentPosition = G.erase(currentPosition);
      } else {
        break;
      }
    }
  }
  detElementsRoad(G, Road,testDirection, roadMakerData,ctx);
}
 

// Map of detector elements production
 
void InDet::SiDetElementsRoadMaker_xk::computeBounds()
{
 
  StatusCode sc;
 
  // Get  Pixel Detector Manager
  //
  const InDetDD::PixelDetectorManager* pixmgr = nullptr;
  if (m_usePIX) {
    sc = detStore()->retrieve(pixmgr, m_pix);
    if (sc.isFailure() || !pixmgr) {
      ATH_MSG_FATAL("Could not get PixelDetectorManager  !");
      return;
    }
  }
 
  // Get  SCT Detector Manager
  //
  const InDetDD::SCT_DetectorManager* sctmgr = nullptr;
  if (m_useSCT) {
    sc = detStore()->retrieve(sctmgr, m_sct);
    if (sc.isFailure() || !sctmgr) {
      ATH_MSG_FATAL("Could not get SCT_DetectorManager !");
      return;
    }
  }
 
  const PixelID* IDp = nullptr;
  const SCT_ID*  IDs = nullptr;
 
  if (m_usePIX && detStore()->retrieve(IDp, "PixelID").isFailure()) {
    ATH_MSG_FATAL("Could not get Pixel ID helper");
  }
 
  if (m_useSCT && detStore()->retrieve(IDs, "SCT_ID").isFailure()) {
    ATH_MSG_FATAL("Could not get SCT ID helper");
  }
 
 
  if (!IDs && !IDp) return;
 
  InDetDD::SiDetectorElementCollection::const_iterator s, se;
  std::vector<InDetDD::SiDetectorElement const*> pW[3];
 
  if (IDp) {
    // Loop over each wafer of pixels
    //
    s  =  pixmgr->getDetectorElementBegin();
    se =  pixmgr->getDetectorElementEnd  ();
 
    for (; s!=se; ++s) {
      if      ((*s)->isBarrel()       ) pW[1].push_back((*s)); // Barrel
      else if ((*s)->center().z() > 0.) pW[2].push_back((*s)); // Right endcap
      else                              pW[0].push_back((*s)); // Left  endcap
    }
  }
 
  if (IDs) {
    // Loop over each wafer of sct
    //
    s  = sctmgr->getDetectorElementBegin();
    se = sctmgr->getDetectorElementEnd  ();
 
    for (; s!=se; ++s) {
      if      ((*s)->isBarrel()       ) pW[1].push_back((*s)); // Barrel
      else if ((*s)->center().z() > 0.) pW[2].push_back((*s)); // Right endcap
      else                              pW[0].push_back((*s)); // Left  endcap
    }
  }
 
  int nel = pW[0].size()+pW[1].size()+pW[2].size();
  if (!nel) return;
 
  std::sort(pW[1].begin(), pW[1].end(), InDet::compDetElements_RAZ());
  std::sort(pW[0].begin(), pW[0].end(), InDet::compDetElements_ZRA());
  std::sort(pW[2].begin(), pW[2].end(), InDet::compDetElements_ZRA());
 
  double   mzmin [3];  // min Z coordinate
  double   mzmax [3];  // max Z coordinate
  double   mrmin [3];  // min radius
  double   mrmax [3];  // max radius
  bool     has[3] {false,false,false};
 
  for (int N=0; N!=3; ++N) {
    double P[40];
    int im    = static_cast<int>(pW[N].size()-1);
    int If    = 0      ;
    double z0 = 0.     ;
    double r0 = 0.     ;
    mrmin[N] = 100000.;
    mrmax[N] =-100000.;
    mzmin[N] = 100000.;
    mzmax[N] =-100000.;
 
    for (int i = 0; i<= im; ++i) {
      InDet::SiDetElementsRoadUtils_xk::detElementInformation(*(pW[N][i]), P);
 
      if (P[ 9] < mrmin[N]) mrmin[N] = P[ 9];
      if (P[10] > mrmax[N]) mrmax[N] = P[10];
      if (P[11] < mzmin[N]) mzmin[N] = P[11];
      if (P[12] > mzmax[N]) mzmax[N] = P[12];
 
      double r    = P[0];
      double z    = P[1];
      bool   newl = false;
      if (N==1) {
    if (fabs(r-r0) > 10.) {
          newl=true;
          r0=r;
        }
      } else {
    if (fabs(z-z0) > 10.) {
          newl=true;
          r0=r;
          z0=z;
        }
      }
 
      if (newl || i==im) {
    int Il = i-1;
        if (i==im) ++Il;
 
    if (If<=Il) {
          has[N]=true;
    }
    If = i;
      }
    }
  }
 
  // CylinderBounds production
  //
  double zmi = +100000;
  double zma = -100000;
  double rma = -100000;
  for (int i=0; i!=3; ++i) {
    if (has[i]) {
      if (mzmin[i]<zmi) zmi=mzmin[i];
      if (mzmax[i]>zma) zma=mzmax[i];
      if (mrmax[i]>rma) rma=mrmax[i];
    }
  }
 
  double hz = fabs(zma);
  if (hz<fabs(zmi)) hz = fabs(zmi);
  const Trk::CylinderBounds CB(rma+20., hz+20.);
  m_bounds  = CB;
}
 

// Distance to detector element according stright line model
 
float InDet::SiDetElementsRoadMaker_xk::stepToDetElement
(const InDetDD::SiDetectorElement*& de, Amg::Vector3D& r, Amg::Vector3D& a)
{
  Amg::Vector3D R = de->center();
  Amg::Vector3D A = de->normal();
  double D = a.x()*A.x()+a.y()*A.y()+a.z()*A.z();
  if (D==0.) return static_cast<float>(D);
  return static_cast<float>((A.x()*(R.x()-r.x())+A.y()*(R.y()-r.y())+A.z()*(R.z()-r.z()))/D);
}

// Cylinder bounds parameters estimation
 
Trk::CylinderBounds InDet::SiDetElementsRoadMaker_xk::getBound
(MagField::AtlasFieldCache& fieldCache,
 const Trk::TrackParameters& Tp) const
{
  const double cor = 1.;
 
  double zfield = 0.;
  if (m_fieldModeEnum!=Trk::NoField && fieldCache.solenoidOn()) {
    const Amg::Vector3D& pos = Tp.position();
    double f[3], p[3] = {pos[Amg::x], pos[Amg::y], pos[Amg::z]};
 
    fieldCache.getFieldZR(p, f);
 
    zfield = 299.7925*f[2];
  }
 
  if (fabs(zfield) < .0000001) return m_bounds;
 
  const AmgVector(5)& Vp = Tp.parameters();
 
  double cur = zfield*Vp[4]/sin(Vp[3]);
 
  if (fabs(cur)*m_bounds.r() < cor) return m_bounds;
 
  double rad = 1./cur;
  if (cor*fabs(rad) > m_bounds.r()  ) return m_bounds;
 
  const  Amg::Vector3D& Gp = Tp.position();
 
  double S, C;
  sincos(Vp[2], &S, &C);
 
  double xc = Gp.x()+S*rad;
  double yc = Gp.y()-C*rad;
  double rm = (sqrt(xc*xc+yc*yc)+fabs(rad))*cor;
  if (rm > m_bounds.r()) return m_bounds;
  Trk::CylinderBounds CB(rm, m_bounds.halflengthZ());
  return CB;
}