SiDetElementsLayer_xk.cxx

Go to the documentation of this file.

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

//   Implementation file for class SiDetElemenstLauerLxk
// (c) ATLAS Detector software
// 
// Version 1.0 21/04/2004 I.Gavrilenko
 
#include <cmath>
 
#include "SiDetElementsRoadTool_xk/SiDetElementsComparison.h"
#include "SiDetElementsRoadTool_xk/SiDetElementsLayer_xk.h"
#include "CxxUtils/trapping_fp.h"

// Get barrel detector elements
// Input parameters: startPoint[0] - X     searchDirection[0] - Ax 
//                   startPoint[1] - Y     searchDirection[1] - Ay
//                   startPoint[2] - Z     searchDirection[2] - Az
//                   startPoint[3] - R
//                   startPoint[4] - width
//                   startPoint[5] - step
 
void InDet::SiDetElementsLayer_xk::getBarrelDetElements
(const std::array<float,6> & startingPoint,
 const std::array<float,3> & searchDirection,
 std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
 std::vector<bool>   &used) const
{
  // Tell clang to optimize assuming that FP exceptions can trap.
  // Otherwise, it can vectorize the division, which can lead to
  // spurious division-by-zero traps from unused vector lanes.
  CXXUTILS_TRAPPING_FP;

  float minusB  =2*( searchDirection[0]*startingPoint[0]
                    +searchDirection[1]*startingPoint[1]); 
  float C  = (m_r-startingPoint[0]-startingPoint[1])
            *(m_r+startingPoint[0]+startingPoint[1])
            +2.*startingPoint[0]*startingPoint[1];
  float twoA  = 2.*( searchDirection[0]*searchDirection[0]
                    +searchDirection[1]*searchDirection[1]); 
  if(twoA == 0.) return;
  float sq = minusB*minusB+2.*C*twoA;  
  if (sq > 0){ 
    sq=std::sqrt(sq);
  }
  else {
    sq=0.;
  }
  float s1 =-(minusB+sq)/twoA;
  float s2 =-(minusB-sq)/twoA; 
  float s;
  if((s1*s2) > 0.) {
    s = (std::abs(s1) < std::abs(s2) ? s1 : s2); 
  }
  else{     
    s = (s1  > 0. ? s1 : s2); 
  }  
  float zc   = startingPoint[2]+searchDirection[2]*s;
  float At   = std::sqrt(1.-searchDirection[2]*searchDirection[2]);
  
  if(At != 0. && std::abs(zc-m_z) > (m_dz+(m_dr*std::abs(searchDirection[2])+startingPoint[4])/At)) return;
  float phiCrossing   = std::atan2(startingPoint[1]+searchDirection[1]*s,startingPoint[0]+searchDirection[0]*s);
  float reducedRoadWidth   = startingPoint[4]/m_r;
  getDetElements(startingPoint,searchDirection,phiCrossing,reducedRoadWidth,lDE,used);
 
}
 
 
void InDet::SiDetElementsLayer_xk::getITkBarrelDetElements
(const std::array<float,6> & startingPoint,
 const std::array<float,3> & searchDirection,
 std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
 std::vector<bool>   &used) const
{
  constexpr float pi = M_PI;
  constexpr float pi2 = 2.*pi;

  float minusB  =2*( searchDirection[0]*startingPoint[0]
                    +searchDirection[1]*startingPoint[1]);
  float C  = (m_r-startingPoint[0]-startingPoint[1])
            *(m_r+startingPoint[0]+startingPoint[1])
            +2.*startingPoint[0]*startingPoint[1];
  float twoA  = 2.*( searchDirection[0]*searchDirection[0]
                    +searchDirection[1]*searchDirection[1]);
  if(twoA == 0.) return;
  float sq = minusB*minusB+2.*C*twoA;
  if (sq > 0){
    sq=std::sqrt(sq);
  }
  else {
    sq=0.;
  }
  float s1 =-(minusB+sq)/twoA;
  float s2 =-(minusB-sq)/twoA;
  float s;
  if((s1*s2) > 0.) {
    s = (std::abs(s1) < std::abs(s2) ? s1 : s2);
  }
  else{
    s = (s1  > 0. ? s1 : s2);
  }
  float zc   = startingPoint[2]+searchDirection[2]*s;
  float At   = std::sqrt(1.-searchDirection[2]*searchDirection[2]);
  float Sz = (m_dr*std::abs(searchDirection[2])+startingPoint[4])/At;

  if(At != 0. && std::abs(zc-m_z) > (m_dz+Sz)) return;
  float reducedRoadWidth = startingPoint[4]/m_r + m_dfe;
  float phiCrossing = std::atan2(startingPoint[1]+searchDirection[1]*s,startingPoint[0]+searchDirection[0]*s) - reducedRoadWidth;
  if(phiCrossing<-pi) phiCrossing += pi2;
 
  int   ie   = m_elements.size();
  int   i1   = ie-1;
 
  if     (phiCrossing <= m_elements[ 0].phi()) i1 = 0 ;
  else if(phiCrossing >  m_elements[i1].phi()) i1 = ie;
  else  {
    int i0=0;
    while((i1-i0)>1) {
      int i = (i0+i1)/2;
      if (m_elements[i].phi() > phiCrossing) i1=i;
      else i0=i;
    }
  }
 
  phiCrossing += 2.*reducedRoadWidth;
 
  for(int i = i1; i!=ie; ++i) {
 
    if( m_elements[i].phi() > phiCrossing ) return;
 
    if(std::abs(zc-m_elements[i].z()) < (m_elements[i].dz()+Sz)
       &&  m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
      lDE.emplace_back(&m_elements[i],s,0); //distance not used for ITk?
      used[i]=true;
    }
  }
 
  phiCrossing -= pi2;
 
  for(int i=0; i!=i1; ++i)  {
 
    if( m_elements[i].phi() > phiCrossing) return;
 
    if(std::abs(zc-m_elements[i].z()) < (m_elements[i].dz()+Sz)
       &&  m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
      lDE.emplace_back(&m_elements[i],s,0);
      used[i]=true;
    }
  }
 
}
 

// Get endcap detector elements
// Input parameters: startingPoint[0] - X     searchDirection[0] - Ax 
//                   startingPoint[1] - Y     searchDirection[1] - Ay
//                   startingPoint[2] - Z     searchDirection[2] - Az
//                   startingPoint[3] - R
//                   startingPoint[4] - width
//                   startingPoint[5] - step
 
void InDet::SiDetElementsLayer_xk::getEndcapDetElements
(const std::array<float,6> & startingPoint,
 const std::array<float,3> & searchDirection,
 std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
 std::vector<bool>   &used) const
{
  float s   =(m_z-startingPoint[2])/searchDirection[2];
  float xc  = startingPoint[0]+searchDirection[0]*s;
  float yc  = startingPoint[1]+searchDirection[1]*s;
  float rc  = std::sqrt(xc*xc+yc*yc);
  float A23 = searchDirection[2]*startingPoint[3];
  if(A23 != 0. && std::abs(rc-m_r) > m_dr+std::abs(2.*(startingPoint[0]*searchDirection[0]+startingPoint[1]*searchDirection[1])*m_dz/A23)+startingPoint[4]) return;
  float phiCrossing  = std::atan2(yc,xc);
  float reducedRoadWidth  = startingPoint[4]/rc;
  getDetElements(startingPoint,searchDirection,phiCrossing,reducedRoadWidth,lDE,used);
 
}
 
 
void InDet::SiDetElementsLayer_xk::getITkEndcapDetElements
(const std::array<float,6> & startingPoint,
 const std::array<float,3> & searchDirection,
 std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
 std::vector<bool>   &used) const
{
  constexpr float pi = M_PI;
  constexpr float pi2 = 2.*pi;

  float s   =(m_z-startingPoint[2])/searchDirection[2];
  float xc  = startingPoint[0]+searchDirection[0]*s;
  float yc  = startingPoint[1]+searchDirection[1]*s;
  float rc  = std::sqrt(xc*xc+yc*yc);
  float A23 = searchDirection[2]*startingPoint[3];
  if(A23 != 0. && std::abs(rc-m_r) > m_dr+std::abs(2.*(startingPoint[0]*searchDirection[0]+startingPoint[1]*searchDirection[1])*m_dz/A23)+startingPoint[4]) return;
  float reducedRoadWidth  = startingPoint[4]/rc + m_dfe;
  float phiCrossing  = std::atan2(yc,xc)-reducedRoadWidth;
  if(phiCrossing<-pi) phiCrossing += pi2;
 
  int   ie   = m_elements.size();
  int   i1   = ie-1;
 
  if     (phiCrossing <= m_elements[ 0].phi()) i1 = 0 ;
  else if(phiCrossing >  m_elements[i1].phi()) i1 = ie;
  else  {
    int i0=0;
    while((i1-i0)>1) {
      int i = (i0+i1)/2;
      if (m_elements[i].phi() > phiCrossing) i1=i;
      else i0=i;
    }
  }
 
  phiCrossing += 2.*reducedRoadWidth;
 
  for(int i = i1; i!=ie; ++i) {
 
    if( m_elements[i].phi() > phiCrossing ) return;
 
    if(m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
      lDE.emplace_back(&m_elements[i],s,0); //distance not used for ITk?
      used[i]=true;
    }
  }
 
  phiCrossing -= pi2;
 
  for(int i=0; i!=i1; ++i)  {
 
    if( m_elements[i].phi() > phiCrossing) return;
 
    if(m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
      lDE.emplace_back(&m_elements[i],s,0);
      used[i]=true;
    }
  }
 
}
 

// Get detector elements
// Input parameters: startingPoint[0] - X     searchDirection[0] - Ax 
//                   startingPoint[1] - Y     searchDirection[1] - Ay
//                   startingPoint[2] - Z     searchDirection[2] - Az
//                   startingPoint[3] - R
//                   startingPoint[4] - width
//                   startingPoint[5] - step
 
void InDet::SiDetElementsLayer_xk::getDetElements
(const std::array<float,6> & startingPoint,
 const std::array<float,3> & searchDirection,
 float phiCrossing,
 float reducedRoadWidth,
 std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
 std::vector<bool>   &used) const
{
  constexpr float pi = M_PI;
  constexpr float pi2 = 2.*pi; 
  int im  = int(m_elements.size())-1; 
  if(im<0) return;
  int i0  = 0;
  int i1  = im;

  if (phiCrossing> m_elements[i0].phi() && phiCrossing< m_elements[i1].phi()) {
    while((i1-i0)>1) {
      int i = (i0+i1)/2; 
      if (m_elements[i].phi() > phiCrossing){
        i1=i;
      }
      else{
        i0=i;
      }
    }
    i0 = i1;
  }
  //
  std::array<float,3> intersectionOutcome{};
  int i = i0;
  while(1) {
    assert( static_cast<unsigned int>(i)<m_elements.size() );
    if(!used[i]) {
      float dPhi =std::abs(m_elements[i].phi()-phiCrossing); 
      if(dPhi>pi) dPhi=std::abs(dPhi-pi2);    

      if((dPhi-reducedRoadWidth)>m_dfe) break;

      m_elements[i].intersect(&(startingPoint[0]),&(searchDirection[0]),&(intersectionOutcome[0]));

      if  (     intersectionOutcome[0]<=startingPoint[4] 
            && (intersectionOutcome[1]<=startingPoint[4])
          ) {
         lDE.emplace_back(&m_elements[i],startingPoint[5]+intersectionOutcome[2],std::max(intersectionOutcome[0],intersectionOutcome[1])); 
         used[i]=true;
      }
    }
    ++i; 
    if(i>im) i=0; 
    if(i==i0) return; 
  }
  i1 = i; i = i0;
  while(1) {
    --i; 
    if(i<0) i=im; 
    if(i==i1) return;
    assert( static_cast<unsigned int>(i)<m_elements.size() );
    if(!used[i]) {
      float dPhi =std::abs(m_elements[i].phi()-phiCrossing); 
      if(dPhi>pi) dPhi=std::abs(dPhi-pi2);
      if((dPhi-reducedRoadWidth)>m_dfe) return;
      m_elements[i].intersect(&(startingPoint[0]),&(searchDirection[0]),&(intersectionOutcome[0]));
 
      if((intersectionOutcome[0]-startingPoint[4])<=0 && (intersectionOutcome[1]-startingPoint[4])<=0.) {
         lDE.emplace_back(&m_elements[i],startingPoint[5]+intersectionOutcome[2],std::max(intersectionOutcome[0],intersectionOutcome[1])); 
         used[i]=true;
      }
    }
  }
}

// Sort detector elements in azimuthal angle order
 
void InDet::SiDetElementsLayer_xk::sortDetectorElements()
{
  std::sort(m_elements.begin(),m_elements.end(),InDet::compDetElements_Azimuthal());
}