ATLAS Offline Software
SiDetElementsLayer_xk.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Implementation file for class SiDetElemenstLauerLxk
8 // (c) ATLAS Detector software
10 //
12 // Version 1.0 21/04/2004 I.Gavrilenko
14 
15 #include <cmath>
16 
17 #include "SiDetElementsRoadTool_xk/SiDetElementsComparison.h"
18 #include "SiDetElementsRoadTool_xk/SiDetElementsLayer_xk.h"
19 #include "CxxUtils/trapping_fp.h"
20 
22 // Get barrel detector elements
23 // Input parameters: startPoint[0] - X searchDirection[0] - Ax
24 // startPoint[1] - Y searchDirection[1] - Ay
25 // startPoint[2] - Z searchDirection[2] - Az
26 // startPoint[3] - R
27 // startPoint[4] - width
28 // startPoint[5] - step
30 
31 void InDet::SiDetElementsLayer_xk::getBarrelDetElements
32 (const std::array<float,6> & startingPoint,
33  const std::array<float,3> & searchDirection,
34  std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
35  std::vector<bool> &used) const
36 {
37  // Tell clang to optimize assuming that FP exceptions can trap.
38  // Otherwise, it can vectorize the division, which can lead to
39  // spurious division-by-zero traps from unused vector lanes.
40  CXXUTILS_TRAPPING_FP;
41 
45  float minusB =2*( searchDirection[0]*startingPoint[0]
46  +searchDirection[1]*startingPoint[1]);
47  float C = (m_r-startingPoint[0]-startingPoint[1])
48  *(m_r+startingPoint[0]+startingPoint[1])
49  +2.*startingPoint[0]*startingPoint[1];
50  float twoA = 2.*( searchDirection[0]*searchDirection[0]
51  +searchDirection[1]*searchDirection[1]);
52  if(twoA == 0.) return;
53  float sq = minusB*minusB+2.*C*twoA;
54  if (sq > 0){
55  sq=std::sqrt(sq);
56  }
57  else {
58  sq=0.;
59  }
61  float s1 =-(minusB+sq)/twoA;
62  float s2 =-(minusB-sq)/twoA;
63  float s;
67  if((s1*s2) > 0.) {
68  s = (std::abs(s1) < std::abs(s2) ? s1 : s2);
69  }
71  else{
72  s = (s1 > 0. ? s1 : s2);
73  }
75  float zc = startingPoint[2]+searchDirection[2]*s;
77  float At = std::sqrt(1.-searchDirection[2]*searchDirection[2]);
78 
84  if(At != 0. && std::abs(zc-m_z) > (m_dz+(m_dr*std::abs(searchDirection[2])+startingPoint[4])/At)) return;
86  float phiCrossing = std::atan2(startingPoint[1]+searchDirection[1]*s,startingPoint[0]+searchDirection[0]*s);
88  float reducedRoadWidth = startingPoint[4]/m_r;
89  getDetElements(startingPoint,searchDirection,phiCrossing,reducedRoadWidth,lDE,used);
90 
91 }
92 
93 
94 void InDet::SiDetElementsLayer_xk::getITkBarrelDetElements
95 (const std::array<float,6> & startingPoint,
96  const std::array<float,3> & searchDirection,
97  std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
98  std::vector<bool> &used) const
99 {
100  constexpr float pi = M_PI;
101  constexpr float pi2 = 2.*pi;
102 
106  float minusB =2*( searchDirection[0]*startingPoint[0]
107  +searchDirection[1]*startingPoint[1]);
108  float C = (m_r-startingPoint[0]-startingPoint[1])
109  *(m_r+startingPoint[0]+startingPoint[1])
110  +2.*startingPoint[0]*startingPoint[1];
111  float twoA = 2.*( searchDirection[0]*searchDirection[0]
112  +searchDirection[1]*searchDirection[1]);
113  if(twoA == 0.) return;
114  float sq = minusB*minusB+2.*C*twoA;
115  if (sq > 0){
116  sq=std::sqrt(sq);
117  }
118  else {
119  sq=0.;
120  }
122  float s1 =-(minusB+sq)/twoA;
123  float s2 =-(minusB-sq)/twoA;
124  float s;
128  if((s1*s2) > 0.) {
129  s = (std::abs(s1) < std::abs(s2) ? s1 : s2);
130  }
132  else{
133  s = (s1 > 0. ? s1 : s2);
134  }
136  float zc = startingPoint[2]+searchDirection[2]*s;
138  float At = std::sqrt(1.-searchDirection[2]*searchDirection[2]);
139  float Sz = (m_dr*std::abs(searchDirection[2])+startingPoint[4])/At;
140 
146  if(At != 0. && std::abs(zc-m_z) > (m_dz+Sz)) return;
148  float reducedRoadWidth = startingPoint[4]/m_r + m_dfe;
150  float phiCrossing = std::atan2(startingPoint[1]+searchDirection[1]*s,startingPoint[0]+searchDirection[0]*s) - reducedRoadWidth;
151  if(phiCrossing<-pi) phiCrossing += pi2;
152 
153  int ie = m_elements.size();
154  int i1 = ie-1;
155 
156  if (phiCrossing <= m_elements[ 0].phi()) i1 = 0 ;
157  else if(phiCrossing > m_elements[i1].phi()) i1 = ie;
158  else {
159  int i0=0;
160  while((i1-i0)>1) {
161  int i = (i0+i1)/2;
162  if (m_elements[i].phi() > phiCrossing) i1=i;
163  else i0=i;
164  }
165  }
166 
167  phiCrossing += 2.*reducedRoadWidth;
168 
169  for(int i = i1; i!=ie; ++i) {
170 
171  if( m_elements[i].phi() > phiCrossing ) return;
172 
173  if(std::abs(zc-m_elements[i].z()) < (m_elements[i].dz()+Sz)
174  && m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
175  lDE.emplace_back(&m_elements[i],s,0); //distance not used for ITk?
176  used[i]=true;
177  }
178  }
179 
180  phiCrossing -= pi2;
181 
182  for(int i=0; i!=i1; ++i) {
183 
184  if( m_elements[i].phi() > phiCrossing) return;
185 
186  if(std::abs(zc-m_elements[i].z()) < (m_elements[i].dz()+Sz)
187  && m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
188  lDE.emplace_back(&m_elements[i],s,0);
189  used[i]=true;
190  }
191  }
192 
193 }
194 
195 
197 // Get endcap detector elements
198 // Input parameters: startingPoint[0] - X searchDirection[0] - Ax
199 // startingPoint[1] - Y searchDirection[1] - Ay
200 // startingPoint[2] - Z searchDirection[2] - Az
201 // startingPoint[3] - R
202 // startingPoint[4] - width
203 // startingPoint[5] - step
205 
206 void InDet::SiDetElementsLayer_xk::getEndcapDetElements
207 (const std::array<float,6> & startingPoint,
208  const std::array<float,3> & searchDirection,
209  std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
210  std::vector<bool> &used) const
211 {
214  float s =(m_z-startingPoint[2])/searchDirection[2];
216  float xc = startingPoint[0]+searchDirection[0]*s;
217  float yc = startingPoint[1]+searchDirection[1]*s;
218  float rc = std::sqrt(xc*xc+yc*yc);
220  float A23 = searchDirection[2]*startingPoint[3];
225  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;
226  float phiCrossing = std::atan2(yc,xc);
227  float reducedRoadWidth = startingPoint[4]/rc;
228  getDetElements(startingPoint,searchDirection,phiCrossing,reducedRoadWidth,lDE,used);
229 
230 }
231 
232 
233 void InDet::SiDetElementsLayer_xk::getITkEndcapDetElements
234 (const std::array<float,6> & startingPoint,
235  const std::array<float,3> & searchDirection,
236  std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
237  std::vector<bool> &used) const
238 {
239  constexpr float pi = M_PI;
240  constexpr float pi2 = 2.*pi;
241 
244  float s =(m_z-startingPoint[2])/searchDirection[2];
246  float xc = startingPoint[0]+searchDirection[0]*s;
247  float yc = startingPoint[1]+searchDirection[1]*s;
248  float rc = std::sqrt(xc*xc+yc*yc);
250  float A23 = searchDirection[2]*startingPoint[3];
255  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;
256  float reducedRoadWidth = startingPoint[4]/rc + m_dfe;
257  float phiCrossing = std::atan2(yc,xc)-reducedRoadWidth;
258  if(phiCrossing<-pi) phiCrossing += pi2;
259 
260  int ie = m_elements.size();
261  int i1 = ie-1;
262 
263  if (phiCrossing <= m_elements[ 0].phi()) i1 = 0 ;
264  else if(phiCrossing > m_elements[i1].phi()) i1 = ie;
265  else {
266  int i0=0;
267  while((i1-i0)>1) {
268  int i = (i0+i1)/2;
269  if (m_elements[i].phi() > phiCrossing) i1=i;
270  else i0=i;
271  }
272  }
273 
274  phiCrossing += 2.*reducedRoadWidth;
275 
276  for(int i = i1; i!=ie; ++i) {
277 
278  if( m_elements[i].phi() > phiCrossing ) return;
279 
280  if(m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
281  lDE.emplace_back(&m_elements[i],s,0); //distance not used for ITk?
282  used[i]=true;
283  }
284  }
285 
286  phiCrossing -= pi2;
287 
288  for(int i=0; i!=i1; ++i) {
289 
290  if( m_elements[i].phi() > phiCrossing) return;
291 
292  if(m_elements[i].intersectITk(&(startingPoint[0]),&(searchDirection[0]),s)) {
293  lDE.emplace_back(&m_elements[i],s,0);
294  used[i]=true;
295  }
296  }
297 
298 }
299 
300 
302 // Get detector elements
303 // Input parameters: startingPoint[0] - X searchDirection[0] - Ax
304 // startingPoint[1] - Y searchDirection[1] - Ay
305 // startingPoint[2] - Z searchDirection[2] - Az
306 // startingPoint[3] - R
307 // startingPoint[4] - width
308 // startingPoint[5] - step
310 
311 void InDet::SiDetElementsLayer_xk::getDetElements
312 (const std::array<float,6> & startingPoint,
313  const std::array<float,3> & searchDirection,
314  float phiCrossing,
315  float reducedRoadWidth,
316  std::vector<InDet::SiDetElementLink_xk::ElementWay> &lDE,
317  std::vector<bool> &used) const
318 {
319  constexpr float pi = M_PI;
320  constexpr float pi2 = 2.*pi;
321  int im = int(m_elements.size())-1;
322  if(im<0) return;
323  int i0 = 0;
324  int i1 = im;
325 
327  if (phiCrossing> m_elements[i0].phi() && phiCrossing< m_elements[i1].phi()) {
328  while((i1-i0)>1) {
329  int i = (i0+i1)/2;
330  if (m_elements[i].phi() > phiCrossing){
331  i1=i;
332  }
333  else{
334  i0=i;
335  }
336  }
337  i0 = i1;
338  }
339  //
340  std::array<float,3> intersectionOutcome{};
341  int i = i0;
343  while(1) {
344  assert( static_cast<unsigned int>(i)<m_elements.size() );
346  if(!used[i]) {
348  float dPhi =std::abs(m_elements[i].phi()-phiCrossing);
349  if(dPhi>pi) dPhi=std::abs(dPhi-pi2);
350 
353  if((dPhi-reducedRoadWidth)>m_dfe) break;
354 
359  m_elements[i].intersect(&(startingPoint[0]),&(searchDirection[0]),&(intersectionOutcome[0]));
360 
362  if ( intersectionOutcome[0]<=startingPoint[4]
363  && (intersectionOutcome[1]<=startingPoint[4])
364  ) {
366  lDE.emplace_back(&m_elements[i],startingPoint[5]+intersectionOutcome[2],std::max(intersectionOutcome[0],intersectionOutcome[1]));
367  used[i]=true;
368  }
369  }
370  ++i;
371  if(i>im) i=0;
372  if(i==i0) return;
373  }
375  i1 = i; i = i0;
376  while(1) {
378  --i;
380  if(i<0) i=im;
382  if(i==i1) return;
383  assert( static_cast<unsigned int>(i)<m_elements.size() );
384  if(!used[i]) {
385  float dPhi =std::abs(m_elements[i].phi()-phiCrossing);
386  if(dPhi>pi) dPhi=std::abs(dPhi-pi2);
387  if((dPhi-reducedRoadWidth)>m_dfe) return;
388  m_elements[i].intersect(&(startingPoint[0]),&(searchDirection[0]),&(intersectionOutcome[0]));
389 
390  if((intersectionOutcome[0]-startingPoint[4])<=0 && (intersectionOutcome[1]-startingPoint[4])<=0.) {
391  lDE.emplace_back(&m_elements[i],startingPoint[5]+intersectionOutcome[2],std::max(intersectionOutcome[0],intersectionOutcome[1]));
392  used[i]=true;
393  }
394  }
395  }
396 }
397 
399 // Sort detector elements in azimuthal angle order
401 
402 void InDet::SiDetElementsLayer_xk::sortDetectorElements()
403 {
404  std::sort(m_elements.begin(),m_elements.end(),InDet::compDetElements_Azimuthal());
405 }
406 
used
InDet::SiDetElementsLayer_xk::getDetElements
void 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
internal helper which resolves the phi-multiplicity of elements within a layer.
Definition: SiDetElementsLayer_xk.cxx:312
CXXUTILS_TRAPPING_FP
#define CXXUTILS_TRAPPING_FP
Definition: trapping_fp.h:24
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
max
#define max(a, b)
Definition: cfImp.cxx:41
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
InDet::SiDetElementsLayer_xk::getBarrelDetElements
void 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
Get barrel detector elements Input parameters: startPoint[0] - X searchDirection[0] - Ax startPoint[1...
Definition: SiDetElementsLayer_xk.cxx:32
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
InDet::compDetElements_Azimuthal
Definition: SiDetElementsComparison.h:107
InDet::SiDetElementsLayer_xk::m_z
float m_z
Definition: SiDetElementsLayer_xk.h:113
SiDetElementsLayer_xk.h
InDet::SiDetElementsLayer_xk::getITkEndcapDetElements
void 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
Definition: SiDetElementsLayer_xk.cxx:234
python.PhysicalConstants.pi2
float pi2
Definition: PhysicalConstants.py:52
M_PI
#define M_PI
Definition: ActiveFraction.h:11
python.atlas_oh.im
im
Definition: atlas_oh.py:167
InDet::SiDetElementsLayer_xk::dz
const float & dz() const
Definition: SiDetElementsLayer_xk.h:55
pi
#define pi
Definition: TileMuonFitter.cxx:65
PlotCalibFromCool.ie
ie
Definition: PlotCalibFromCool.py:420
InDet::SiDetElementsLayer_xk::getITkBarrelDetElements
void 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
Definition: SiDetElementsLayer_xk.cxx:95
lumiFormat.i
int i
Definition: lumiFormat.py:92
sq
#define sq(x)
Definition: CurvedSegmentFinder.cxx:6
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:530
InDet::SiDetElementsLayer_xk::getEndcapDetElements
void 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
Definition: SiDetElementsLayer_xk.cxx:207
InDet::SiDetElementsLayer_xk::sortDetectorElements
void sortDetectorElements()
Definition: SiDetElementsLayer_xk.cxx:402
InDet::SiDetElementsLayer_xk::m_dz
float m_dz
Definition: SiDetElementsLayer_xk.h:114
SiDetElementsComparison.h
trapping_fp.h
Tell the compiler to optimize assuming that FP may trap.
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
InDet::SiDetElementsLayer_xk::m_dfe
float m_dfe
Definition: SiDetElementsLayer_xk.h:117
InDet::SiDetElementsLayer_xk::z
const float & z() const
Definition: SiDetElementsLayer_xk.h:54
InDet::SiDetElementsLayer_xk::m_elements
std::vector< SiDetElementLink_xk > m_elements
Definition: SiDetElementsLayer_xk.h:118
InDet::SiDetElementsLayer_xk::m_dr
float m_dr
Definition: SiDetElementsLayer_xk.h:116
InDet::SiDetElementsLayer_xk::m_r
float m_r
Definition: SiDetElementsLayer_xk.h:115
Generated on Fri Jul 5 2024 21:26:54 for ATLAS Offline Software by   doxygen 1.8.18