ATLAS Offline Software
SiDetElementsRoadMaker_xk.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Implementation file for class InDet::SiDetElementsRoadMaker_xk
8 // (c) ATLAS Detector software
11 // Version 1.0 21/04/2004 I.Gavrilenko
13 
14 #include "SiDetElementsRoadTool_xk/SiDetElementsRoadMaker_xk.h"
15 #include "SiSPSeededTrackFinderData/SiDetElementRoadMakerData_xk.h"
16 
17 #include "SiDetElementsRoadUtils_xk.h"
18 
19 #include "GaudiKernel/ContextSpecificPtr.h"
20 #include "PixelReadoutGeometry/PixelDetectorManager.h"
21 #include "SCT_ReadoutGeometry/SCT_DetectorManager.h"
22 #include "SiDetElementsRoadTool_xk/SiDetElementsComparison.h"
23 #include "StoreGate/ReadCondHandle.h"
24 #include "TrkPrepRawData/PrepRawData.h"
25 #include <cmath>
26 
27 #include <ostream>
28 #include <iomanip>
29 
31 // Constructor
33 
34 InDet::SiDetElementsRoadMaker_xk::SiDetElementsRoadMaker_xk
35 (const std::string& t, const std::string& n, const IInterface* p)
36  : base_class(t, n, p)
37 {
38 }
39 
41 // Initialisation
43 
44 StatusCode InDet::SiDetElementsRoadMaker_xk::initialize()
45 {
46  //Class optimization checks
47  static_assert(std::is_trivially_copyable<SiDetElementLink_xk::ElementWay>::value);
48  static_assert(std::is_trivially_destructible<SiDetElementLink_xk::ElementWay>::value);
49  static_assert(std::is_trivially_copyable<SiDetElementLink_xk>::value);
50  static_assert(std::is_trivially_destructible<SiDetElementLink_xk>::value);
51  static_assert(std::is_nothrow_move_constructible<SiDetElementsLayer_xk>::value);
52 
53 
54  if (!m_usePIX && !m_useSCT) {
55  ATH_MSG_FATAL("Please don't call this tool if usePixel and useSCT are false");
56  return StatusCode::FAILURE;
57  }
58 
59  if (m_fieldmode == "NoField") m_fieldModeEnum = Trk::NoField;
60  else if (m_fieldmode == "MapSolenoid") m_fieldModeEnum = Trk::FastField;
61  else m_fieldModeEnum = Trk::FullField;
62 
63  // Get propagator tool
64  //
65  ATH_CHECK(m_proptool.retrieve());
66 
67  // Get output print level
68  //
69  m_outputlevel = msg().level()-MSG::DEBUG;
70 
71  computeBounds();
72 
73  ATH_CHECK(m_layerVecKey.initialize());
74  ATH_CHECK(m_fieldCondObjInputKey.initialize());
75 
76  return StatusCode::SUCCESS;
77 }
78 
80 // Finalize
82 
83 StatusCode InDet::SiDetElementsRoadMaker_xk::finalize()
84 {
85  return StatusCode::SUCCESS;
86 }
87 
88 
90 // Dumps relevant information into the MsgStream
92 
93 MsgStream& InDet::SiDetElementsRoadMaker_xk::dump(MsgStream& out) const
94 {
95  out<<"\n";
96  return dumpConditions(out);
97 }
98 
100 // Dumps conditions information into the MsgStream
102 
103 MsgStream& InDet::SiDetElementsRoadMaker_xk::dumpConditions(MsgStream& out) const
104 {
105  int n = 62-m_proptool.type().size();
106  std::string s1;
107  for (int i=0; i<n; ++i) s1.append(" ");
108  s1.append("|");
109 
110  std::string fieldmode[9] = {"NoField" , "ConstantField", "SolenoidalField",
111  "ToroidalField" , "Grid3DField" , "RealisticField" ,
112  "UndefinedField", "AthenaField" , "?????" };
113 
114  Trk::MagneticFieldMode fieldModeEnum(m_fieldModeEnum);
115 
116  SG::ReadCondHandle<AtlasFieldCacheCondObj> fieldHandle(m_fieldCondObjInputKey);
117  const AtlasFieldCacheCondObj* fieldCondObj(*fieldHandle);
118  if (fieldCondObj) {
119  MagField::AtlasFieldCache fieldCache;
120  fieldCondObj->getInitializedCache(fieldCache);
121  if (!fieldCache.solenoidOn()) fieldModeEnum = Trk::NoField;
122  }
123 
124  Trk::MagneticFieldProperties fieldprop(fieldModeEnum);
125 
126  int mode = fieldprop.magneticFieldMode();
127  if (mode<0 || mode>8) mode = 8;
128 
129  n = 62-fieldmode[mode].size();
130  std::string s3;
131  for (int i=0; i<n; ++i) s3.append(" ");
132  s3.append("|");
133 
134  n = 62-m_sct.size();
135  std::string s5;
136  for (int i=0; i<n; ++i) s5.append(" ");
137  s5.append("|");
138 
139  n = 62-m_pix.size();
140  std::string s6;
141  for (int i=0; i<n; ++i) s6.append(" ");
142  s6.append("|");
143 
144  const EventContext& ctx = Gaudi::Hive::currentContext();
145 
146  const SiDetElementsLayerVectors_xk &layer = *getLayers(ctx);
147 
148  int maps = 0;
149  if (!layer[0].empty()) ++maps;
150  if (!layer[1].empty()) ++maps;
151  if (!layer[2].empty()) ++maps;
152  out<<"|----------------------------------------------------------------------"
153  <<"-------------------|"
154  <<"\n";
155  if (m_useSCT) {
156  out<<"| SCT detector manager | "<<m_sct <<s5<<"\n";
157  }
158  if (m_usePIX) {
159  out<<"| Pixel detector manager | "<<m_pix <<s6<<"\n";
160  }
161  out<<"| Tool for propagation | "<<m_proptool.type() <<s1<<"\n";
162  out<<"| Magnetic field mode | "<<fieldmode[mode] <<s3<<"\n";
163  out<<"| Width of the road (mm) | "
164  <<std::setw(12)<<std::setprecision(5)<<m_width
165  <<" |"<<"\n";
166  out<<"|----------------------------------------------------------------------"
167  <<"-------------------|"
168  <<"\n";
169 
170  if (!maps || m_outputlevel==0) return out;
171 
172  if (!layer[1].empty()) {
173  int nl = layer[1].size();
174  int nc = 0;
175  for (const auto & i : layer[1]) nc+=i.nElements();
176  out<<"|----------------------------------------------------------------|"
177  <<"\n";
178  out<<"| Barrel map containt "
179  <<std::setw(3)<<nl<<" layers and"
180  <<std::setw(5)<<nc<<" elements |"
181  <<"\n";
182  out<<"|------|-----------|------------|------------|------------|------|"
183  <<"\n";
184  out<<"| n | R | Z min | Z max | max dF | nEl |"
185  <<"\n";
186  out<<"|------|-----------|------------|------------|------------|------|"
187  <<"\n";
188  for (unsigned int i=0; i!=layer[1].size(); ++i) {
189  double zmin = layer[1].at(i).z()-layer[1].at(i).dz();
190  double zmax = layer[1].at(i).z()+layer[1].at(i).dz();
191  out<<"| "
192  <<std::setw(4)<<i<<" |"
193  <<std::setw(10)<<std::setprecision(4)<< layer[1].at(i).r ()<<" | "
194  <<std::setw(10)<<std::setprecision(4)<< zmin<<" | "
195  <<std::setw(10)<<std::setprecision(4)<< zmax<<" | "
196  <<std::setw(10)<<std::setprecision(4)<< layer[1].at(i).dfe()<<" | "
197  <<std::setw(4)<<layer[1].at(i).nElements()<<" | "
198  <<"\n";
199  }
200  out<<"|------|-----------|------------|------------|------------|------|"
201  <<"\n";
202 
203  }
204  if (!layer[0].empty()) {
205  int nl = layer[0].size();
206  int nc = 0;
207  for (const auto & i : layer[0]) nc+=i.nElements();
208  out<<"|----------------------------------------------------------------|"
209  <<"\n";
210  out<<"| L.Endcap map containt"
211  <<std::setw(3)<<nl<<" layers and"
212  <<std::setw(5)<<nc<<" elements |"
213  <<"\n";
214 
215  out<<"|------|-----------|------------|------------|------------|------|"
216  <<"\n";
217  out<<"| n | Z | R min | R max | max dF | nEl |"
218  <<"\n";
219  out<<"|------|-----------|------------|------------|------------|------|"
220  <<"\n";
221  for (unsigned int i=0; i!=layer[0].size(); ++i) {
222  double rmin = layer[0].at(i).r()-layer[0].at(i).dr();
223  double rmax = layer[0].at(i).r()+layer[0].at(i).dr();
224  out<<"| "
225  <<std::setw(4)<<i<<" |"
226  <<std::setw(10)<<std::setprecision(4)<< layer[0].at(i).z()<<" | "
227  <<std::setw(10)<<std::setprecision(4)<< rmin<<" | "
228  <<std::setw(10)<<std::setprecision(4)<< rmax<<" | "
229  <<std::setw(10)<<std::setprecision(4)<<layer[0].at(i).dfe()<<" | "
230  <<std::setw(4)<<layer[0].at(i).nElements()<<" | "
231  <<"\n";
232  }
233  out<<"|------|-----------|------------|------------|------------|------|"
234  <<"\n";
235  }
236  if (!layer[2].empty()) {
237  int nl = layer[2].size();
238  int nc = 0;
239  for (const auto & i : layer[2]) nc+=i.nElements();
240  out<<"|----------------------------------------------------------------|"
241  <<"\n";
242  out<<"| R.Endcap map containt"
243  <<std::setw(3)<<nl<<" layers and"
244  <<std::setw(5)<<nc<<" elements |"
245  <<"\n";
246  out<<"|------|-----------|------------|------------|------------|------|"
247  <<"\n";
248  out<<"| n | Z | R min | R max | max dF | nEl |"
249  <<"\n";
250  out<<"|------|-----------|------------|------------|------------|------|"
251  <<"\n";
252  for (unsigned int i=0; i!=layer[2].size(); ++i) {
253  double rmin = layer[2].at(i).r()-layer[0].at(i).dr();
254  double rmax = layer[2].at(i).r()+layer[0].at(i).dr();
255  out<<"| "
256  <<std::setw(4)<<i<<" |"
257  <<std::setw(10)<<std::setprecision(4)<< layer[2].at(i).z()<<" | "
258  <<std::setw(10)<<std::setprecision(4)<< rmin<<" | "
259  <<std::setw(10)<<std::setprecision(4)<< rmax<<" | "
260  <<std::setw(10)<<std::setprecision(4)<<layer[2].at(i).dfe()<<" | "
261  <<std::setw(4)<<layer[2].at(i).nElements()<<" | "
262  <<"\n";
263  }
264  out<<"|------|-----------|------------|------------|------------|------|"
265  <<"\n";
266  }
267  out<<"\n";
268  return out;
269 }
270 
272 // Dumps relevant information into the ostream
274 
275 std::ostream& InDet::SiDetElementsRoadMaker_xk::dump(std::ostream& out) const
276 {
277  return out;
278 }
279 
281 // Overload of << operator MsgStream
283 
284 MsgStream& InDet::operator <<
285 (MsgStream& sl, const InDet::SiDetElementsRoadMaker_xk& se)
286 {
287  return se.dump(sl);
288 }
289 
291 // Overload of << operator std::ostream
293 
294 std::ostream& InDet::operator <<
295 (std::ostream& sl, const InDet::SiDetElementsRoadMaker_xk& se)
296 {
297  return se.dump(sl);
298 }
299 
300 
302 // Main methods for road builder using input list global positions
304 
305 void InDet::SiDetElementsRoadMaker_xk::detElementsRoad
306 (std::deque<Amg::Vector3D>& globalPositions,
307  std::vector<const InDetDD::SiDetectorElement*>& Road,
308  bool testDirection,
309  SiDetElementRoadMakerData_xk & roadMakerData,
310  const EventContext& ctx) const
311 {
312  if (!m_usePIX && !m_useSCT) return;
313 
325  const SiDetElementsLayerVectors_xk &layer = *getLayers(ctx);
326 
328  std::deque<Amg::Vector3D>::iterator currentPosition=globalPositions.begin(), endPositions=globalPositions.end();
329 
332  std::array<float,6> par_startingPoint{static_cast<float>((*currentPosition).x()), // x of first position
333  static_cast<float>((*currentPosition).y()), // y of first position
334  static_cast<float>((*currentPosition).z()), // Z of first position
335  static_cast<float>(sqrt((*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y())), // r of first position
336  m_width, // road width
337  0.};
338 
342  int n0 = 0;
343  for (; n0!=static_cast<int>(layer[0].size()); ++n0) {
344  if (par_startingPoint[2] > layer[0][n0].z()) break;
345  }
346 
350  int n1 = 0;
351  for (; n1!=static_cast<int>(layer[1].size()); ++n1) {
352  if (par_startingPoint[3] < layer[1][n1].r()) break;
353  }
356  int n2 = 0;
357  for (; n2!=static_cast<int>(layer[2].size()); ++n2) {
358  if (par_startingPoint[2] < layer[2][n2].z()) break;
359  }
360 
361 
362 
366  if (!roadMakerData.isInitialized){
367  bookUsageTracker(roadMakerData,layer);
368  }
369  else{
372  roadMakerData.resetUsageTracker();
373  }
374 
375  std::vector<InDet::SiDetElementLink_xk::ElementWay> lDE;
376  lDE.reserve(8); //reasonable minimum guess
378  ++currentPosition;
379  while (currentPosition!=endPositions) {
381  std::array<float,4> par_targetPoint{static_cast<float>((*currentPosition).x()),
382  static_cast<float>((*currentPosition).y()),
383  static_cast<float>((*currentPosition).z()),
384  static_cast<float>(sqrt((*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y()))
385  };
386 
388  float dx = par_targetPoint[0]-par_startingPoint[0];
389  float dy = par_targetPoint[1]-par_startingPoint[1];
390  float dz = par_targetPoint[2]-par_startingPoint[2];
391  float dist3D = std::sqrt(dx*dx+dy*dy+dz*dz);
392  if (dist3D <=0.) {
393  ++currentPosition;
394  continue;
395  }
396  float inverseDistance = 1./dist3D;
397  std::array<float,3> searchDirection{dx*inverseDistance, dy*inverseDistance, dz*inverseDistance};
399 
403  float unitSepTransverseComp = searchDirection[0]*searchDirection[0]+searchDirection[1]*searchDirection[1];
404  float dr = 0. ;
405  if (unitSepTransverseComp!=0.) {
408  float sm = -( searchDirection[0]*par_startingPoint[0] +
409  searchDirection[1]*par_startingPoint[1])
410  /unitSepTransverseComp;
411 
417  if (sm > 1. && sm < dist3D) {
418  par_targetPoint[0] = par_startingPoint[0]+searchDirection[0]*sm;
419  par_targetPoint[1] = par_startingPoint[1]+searchDirection[1]*sm;
420  par_targetPoint[2] = par_startingPoint[2]+searchDirection[2]*sm;
421  par_targetPoint[3] = std::sqrt(par_targetPoint[0]*par_targetPoint[0]+par_targetPoint[1]*par_targetPoint[1]);
424  dr = 20.;
425  } else {
426  ++currentPosition;
427  }
428  } else {
429  ++currentPosition;
430  }
431 
433 
435 
437  if (par_targetPoint[3]>par_startingPoint[3]) {
439  for (; n1<static_cast<int>(layer[1].size()); ++n1) {
441  if (par_targetPoint[3] < layer[1][n1].r()) break;
442  assert( roadMakerData.elementUsageTracker[1].size() > static_cast<unsigned int>(n1) );
444  if(m_ITkGeometry) layer[1][n1].getITkBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
445  else layer[1][n1].getBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
446  }
448  } else {
449  for (--n1; n1>=0; --n1) {
451  if (par_targetPoint[3] > layer[1][n1].r()+dr) break;
452  assert( roadMakerData.elementUsageTracker[1].size() > static_cast<unsigned int>(n1) );
454  if(m_ITkGeometry) layer[1][n1].getITkBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
455  else layer[1][n1].getBarrelDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[1][n1]);
456 
457  }
458  ++n1;
459  }
460 
463  if (par_targetPoint[2]>par_startingPoint[2]) {
464  for (; n2<static_cast<int>(layer[2].size()); ++n2) {
465  if (par_targetPoint[2] < layer[2][n2].z()) break;
466  assert( roadMakerData.elementUsageTracker[2].size() > static_cast<unsigned int>(n2) );
468  if(m_ITkGeometry) layer[2][n2].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[2][n2]);
469  else layer[2][n2].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[2][n2]);
470  }
471  } else {
472  for (--n2; n2>=0; --n2) {
473  if (par_targetPoint[2] > layer[2][n2].z()) break;
474  assert( roadMakerData.elementUsageTracker[2].size() > static_cast<unsigned int>(n2) );
476  if(m_ITkGeometry) layer[2][n2].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[2][n2]);
477  else layer[2][n2].getEndcapDetElements(par_startingPoint, searchDirection, lDE, roadMakerData.elementUsageTracker[2][n2]);
478  }
479  ++n2;
480  }
481 
484  if (par_targetPoint[2]<par_startingPoint[2]) {
485  for (; n0<static_cast<int>(layer[0].size()); ++n0) {
486  if (par_targetPoint[2] > layer[0][n0].z()) break;
487  assert( roadMakerData.elementUsageTracker[0].size() > static_cast<unsigned int>(n0) );
489  if(m_ITkGeometry) layer[0][n0].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
490  else layer[0][n0].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
491  }
492  } else {
493  for (--n0; n0>=0; --n0) {
494  if (par_targetPoint[2] < layer[0][n0].z()) break;
495  assert( roadMakerData.elementUsageTracker[0].size() > static_cast<unsigned int>(n0) );
497  if(m_ITkGeometry) layer[0][n0].getITkEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
498  else layer[0][n0].getEndcapDetElements(par_startingPoint, searchDirection, lDE,roadMakerData.elementUsageTracker[0][n0]);
499  }
500  ++n0;
501  }
503  par_startingPoint[0] = par_targetPoint[0];
504  par_startingPoint[1] = par_targetPoint[1];
505  par_startingPoint[2] = par_targetPoint[2];
506  par_startingPoint[3] = par_targetPoint[3];
508  par_startingPoint[5]+= dist3D;
509  }
510  auto vec2 = lDE;
511  std::sort(lDE.begin(),lDE.end(),InDet::compDetElementWays());
512  // Fill pointers to detector elements
513  Road.reserve(lDE.size());
514  for (auto & d : lDE){
515  if (testDirection && d.way() < 0) {continue;}
516  Road.push_back(d.link()->detElement());
517  }
518 }
519 
521 void InDet::SiDetElementsRoadMaker_xk::bookUsageTracker(InDet::SiDetElementRoadMakerData_xk & data, const SiDetElementsLayerVectors_xk &layers) {
522 
525  for ( unsigned int side_i=0; side_i<3; ++side_i) {
526  data.elementUsageTracker[side_i].resize( layers[side_i].size() );
527  for (unsigned int layer_i=0; layer_i < layers[side_i].size(); ++layer_i) {
529  data.elementUsageTracker[side_i][layer_i].resize( layers[side_i][layer_i].nElements() );
530  }
531  }
532  data.isInitialized=true;
533 }
534 
536 // Main methods for road builder using track parameters and direction
538 
539 void InDet::SiDetElementsRoadMaker_xk::detElementsRoad
540 (const EventContext& ctx,
541  MagField::AtlasFieldCache& fieldCache,
542  const Trk::TrackParameters& Tp,
543  Trk::PropDirection direction,
544  std::vector<const InDetDD::SiDetectorElement*>& Road,
545  SiDetElementRoadMakerData_xk & roadMakerData) const
546 {
547  if (!m_usePIX && !m_useSCT) return;
549  double qp = fabs(500.*Tp.parameters()[4]);
551  if (qp < 1.e-10) qp = 1.e-10;
553  double S = m_step/qp;
555  if (S > 1000. ) S = 1000. ;
556 
557  bool testDirection = true;
558  if (direction<0) {
559  testDirection = false;
560  S=-S;
561  }
562 
563  Trk::MagneticFieldMode fieldModeEnum(m_fieldModeEnum);
564  if (!fieldCache.solenoidOn()) fieldModeEnum = Trk::NoField;
565  Trk::MagneticFieldProperties fieldprop(fieldModeEnum);
566 
567  // Note: could also give fieldCache directly to propagator if it would be more efficient - would
568  // need to add interface RDS 2020/03
569  std::deque<Amg::Vector3D> G;
570 
573  m_proptool->globalPositions(ctx, G, Tp, fieldprop,getBound(fieldCache, Tp), S, Trk::pion);
575  if (G.size()<2) return;
576 
579  if (direction > 0) {
580  std::deque<Amg::Vector3D>::iterator currentPosition=G.begin(), nextPosition, endPositions=G.end();
581  float r0 = (*currentPosition).x()*(*currentPosition).x()+(*currentPosition).y()*(*currentPosition).y();
582 
583  while (currentPosition!=endPositions) {
584  nextPosition = currentPosition;
585  if (++nextPosition == endPositions) break;
586 
587  float r = (*nextPosition).x()*(*nextPosition).x()+(*nextPosition).y()*(*nextPosition).y();
589  if (r < r0) {
590  r0 = r;
591  currentPosition = G.erase(currentPosition);
592  } else {
593  break;
594  }
595  }
596  }
598  detElementsRoad(G, Road,testDirection, roadMakerData,ctx);
599 }
600 
601 
603 // Map of detector elements production
605 
606 void InDet::SiDetElementsRoadMaker_xk::computeBounds()
607 {
608 
609  StatusCode sc;
610 
611  // Get Pixel Detector Manager
612  //
613  const InDetDD::PixelDetectorManager* pixmgr = nullptr;
614  if (m_usePIX) {
615  sc = detStore()->retrieve(pixmgr, m_pix);
616  if (sc.isFailure() || !pixmgr) {
617  ATH_MSG_FATAL("Could not get PixelDetectorManager !");
618  return;
619  }
620  }
621 
622  // Get SCT Detector Manager
623  //
624  const InDetDD::SCT_DetectorManager* sctmgr = nullptr;
625  if (m_useSCT) {
626  sc = detStore()->retrieve(sctmgr, m_sct);
627  if (sc.isFailure() || !sctmgr) {
628  ATH_MSG_FATAL("Could not get SCT_DetectorManager !");
629  return;
630  }
631  }
632 
633  const PixelID* IDp = nullptr;
634  const SCT_ID* IDs = nullptr;
635 
636  if (m_usePIX && detStore()->retrieve(IDp, "PixelID").isFailure()) {
637  ATH_MSG_FATAL("Could not get Pixel ID helper");
638  }
639 
640  if (m_useSCT && detStore()->retrieve(IDs, "SCT_ID").isFailure()) {
641  ATH_MSG_FATAL("Could not get SCT ID helper");
642  }
643 
644 
645  if (!IDs && !IDp) return;
646 
647  InDetDD::SiDetectorElementCollection::const_iterator s, se;
648  std::vector<InDetDD::SiDetectorElement*> pW[3];
649 
650  if (IDp) {
651  // Loop over each wafer of pixels
652  //
653  s = pixmgr->getDetectorElementBegin();
654  se = pixmgr->getDetectorElementEnd ();
655 
656  for (; s!=se; ++s) {
657  if ((*s)->isBarrel() ) pW[1].push_back((*s)); // Barrel
658  else if ((*s)->center().z() > 0.) pW[2].push_back((*s)); // Right endcap
659  else pW[0].push_back((*s)); // Left endcap
660  }
661  }
662 
663  if (IDs) {
664  // Loop over each wafer of sct
665  //
666  s = sctmgr->getDetectorElementBegin();
667  se = sctmgr->getDetectorElementEnd ();
668 
669  for (; s!=se; ++s) {
670  if ((*s)->isBarrel() ) pW[1].push_back((*s)); // Barrel
671  else if ((*s)->center().z() > 0.) pW[2].push_back((*s)); // Right endcap
672  else pW[0].push_back((*s)); // Left endcap
673  }
674  }
675 
676  int nel = pW[0].size()+pW[1].size()+pW[2].size();
677  if (!nel) return;
678 
679  std::sort(pW[1].begin(), pW[1].end(), InDet::compDetElements_RAZ());
680  std::sort(pW[0].begin(), pW[0].end(), InDet::compDetElements_ZRA());
681  std::sort(pW[2].begin(), pW[2].end(), InDet::compDetElements_ZRA());
682 
683  double mzmin [3]; // min Z coordinate
684  double mzmax [3]; // max Z coordinate
685  double mrmin [3]; // min radius
686  double mrmax [3]; // max radius
687  bool has[3] {false,false,false};
688 
689  for (int N=0; N!=3; ++N) {
690  double P[40];
691  int im = static_cast<int>(pW[N].size()-1);
692  int If = 0 ;
693  double z0 = 0. ;
694  double r0 = 0. ;
695  mrmin[N] = 100000.;
696  mrmax[N] =-100000.;
697  mzmin[N] = 100000.;
698  mzmax[N] =-100000.;
699 
700  for (int i = 0; i<= im; ++i) {
701  InDet::SiDetElementsRoadUtils_xk::detElementInformation(*(pW[N][i]), P);
702 
703  if (P[ 9] < mrmin[N]) mrmin[N] = P[ 9];
704  if (P[10] > mrmax[N]) mrmax[N] = P[10];
705  if (P[11] < mzmin[N]) mzmin[N] = P[11];
706  if (P[12] > mzmax[N]) mzmax[N] = P[12];
707 
708  double r = P[0];
709  double z = P[1];
710  bool newl = false;
711  if (N==1) {
712  if (fabs(r-r0) > 10.) {
713  newl=true;
714  r0=r;
715  }
716  } else {
717  if (fabs(z-z0) > 10.) {
718  newl=true;
719  r0=r;
720  z0=z;
721  }
722  }
723 
724  if (newl || i==im) {
725  int Il = i-1;
726  if (i==im) ++Il;
727 
728  if (If<=Il) {
729  has[N]=true;
730  }
731  If = i;
732  }
733  }
734  }
735 
736  // CylinderBounds production
737  //
738  double zmi = +100000;
739  double zma = -100000;
740  double rma = -100000;
741  for (int i=0; i!=3; ++i) {
742  if (has[i]) {
743  if (mzmin[i]<zmi) zmi=mzmin[i];
744  if (mzmax[i]>zma) zma=mzmax[i];
745  if (mrmax[i]>rma) rma=mrmax[i];
746  }
747  }
748 
749  double hz = fabs(zma);
750  if (hz<fabs(zmi)) hz = fabs(zmi);
751  const Trk::CylinderBounds CB(rma+20., hz+20.);
752  m_bounds = CB;
753 }
754 
755 
757 // Distance to detector element according stright line model
759 
760 float InDet::SiDetElementsRoadMaker_xk::stepToDetElement
761 (const InDetDD::SiDetectorElement*& de, Amg::Vector3D& r, Amg::Vector3D& a)
762 {
763  Amg::Vector3D R = de->center();
764  Amg::Vector3D A = de->normal();
765  double D = a.x()*A.x()+a.y()*A.y()+a.z()*A.z();
766  if (D==0.) return static_cast<float>(D);
767  return static_cast<float>((A.x()*(R.x()-r.x())+A.y()*(R.y()-r.y())+A.z()*(R.z()-r.z()))/D);
768 }
769 
771 // Cylinder bounds parameters estimation
773 
774 Trk::CylinderBounds InDet::SiDetElementsRoadMaker_xk::getBound
775 (MagField::AtlasFieldCache& fieldCache,
776  const Trk::TrackParameters& Tp) const
777 {
778  const double cor = 1.;
779 
780  double zfield = 0.;
781  if (m_fieldModeEnum!=Trk::NoField && fieldCache.solenoidOn()) {
782  const Amg::Vector3D& pos = Tp.position();
783  double f[3], p[3] = {pos[Amg::x], pos[Amg::y], pos[Amg::z]};
784 
785  fieldCache.getFieldZR(p, f);
786 
787  zfield = 299.7925*f[2];
788  }
789 
790  if (fabs(zfield) < .0000001) return m_bounds;
791 
792  const AmgVector(5)& Vp = Tp.parameters();
793 
794  double cur = zfield*Vp[4]/sin(Vp[3]);
795 
796  if (fabs(cur)*m_bounds.r() < cor) return m_bounds;
797 
798  double rad = 1./cur;
799  if (cor*fabs(rad) > m_bounds.r() ) return m_bounds;
800 
801  const Amg::Vector3D& Gp = Tp.position();
802 
803  double S, C;
804  sincos(Vp[2], &S, &C);
805 
806  double xc = Gp.x()+S*rad;
807  double yc = Gp.y()-C*rad;
808  double rm = (sqrt(xc*xc+yc*yc)+fabs(rad))*cor;
809  if (rm > m_bounds.r()) return m_bounds;
810  Trk::CylinderBounds CB(rm, m_bounds.halflengthZ());
811  return CB;
812 }
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
SiDetElementRoadMakerData_xk.h
beamspotman.r
def r
Definition: beamspotman.py:676
InDet::SiDetElementsRoadMaker_xk::m_proptool
PublicToolHandle< Trk::IPropagator > m_proptool
Definition: SiDetElementsRoadMaker_xk.h:135
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
D3PDMakerTestInstan::vec2
std::vector< D3PDTest::MyVec2 > vec2
Definition: D3PDMakerTestDict.h:14
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
SiDetElementsRoadUtils_xk.h
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
InDetDD::SCT_DetectorManager
Definition: SCT_DetectorManager.h:49
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
Trk::MagneticFieldProperties
Definition: MagneticFieldProperties.h:31
AtlasFieldCacheCondObj
Definition: AtlasFieldCacheCondObj.h:19
InDet::SiDetElementsRoadMaker_xk::initialize
virtual StatusCode initialize() override
Definition: SiDetElementsRoadMaker_xk.cxx:44
MagField::AtlasFieldCache::getFieldZR
void getFieldZR(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field valaue on the z-r plane at given position works only inside the solenoid.
Definition: AtlasFieldCache.cxx:86
beamspotman.cur
def cur
Definition: beamspotman.py:671
InDet::SiDetElementRoadMakerData_xk::elementUsageTracker
ElementUsageTracker elementUsageTracker
Definition: SiDetElementRoadMakerData_xk.h:56
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
PixelAthClusterMonAlgCfg.zmin
zmin
Definition: PixelAthClusterMonAlgCfg.py:176
hist_file_dump.d
d
Definition: hist_file_dump.py:137
InDetDD::SolidStateDetectorElementBase::center
virtual const Amg::Vector3D & center() const override final
Center in global coordinates.
Amg::y
@ y
Definition: GeoPrimitives.h:35
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
InDet::SiDetElementsRoadMaker_xk::finalize
virtual StatusCode finalize() override
Definition: SiDetElementsRoadMaker_xk.cxx:83
module_driven_slicing.layers
layers
Definition: module_driven_slicing.py:114
InDetDD::SCT_DetectorManager::getDetectorElementBegin
virtual SiDetectorElementCollection::const_iterator getDetectorElementBegin() const override
Definition: SCT_DetectorManager.cxx:93
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
python.atlas_oh.im
im
Definition: atlas_oh.py:167
InDet::SiDetElementsRoadMaker_xk::m_ITkGeometry
BooleanProperty m_ITkGeometry
Definition: SiDetElementsRoadMaker_xk.h:157
athena.value
value
Definition: athena.py:122
InDet::SiDetElementsRoadMaker_xk::m_useSCT
BooleanProperty m_useSCT
Definition: SiDetElementsRoadMaker_xk.h:151
InDet::SiDetElementsRoadMaker_xk::stepToDetElement
static float stepToDetElement(const InDetDD::SiDetectorElement *&, Amg::Vector3D &, Amg::Vector3D &)
Definition: SiDetElementsRoadMaker_xk.cxx:761
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
InDet::SiDetElementsRoadMaker_xk::m_fieldmode
StringProperty m_fieldmode
Definition: SiDetElementsRoadMaker_xk.h:156
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
InDet::SiDetElementsRoadMaker_xk
Definition: SiDetElementsRoadMaker_xk.h:61
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
InDet::SiDetElementsRoadMaker_xk::dump
MsgStream & dump(MsgStream &out) const override
Definition: SiDetElementsRoadMaker_xk.cxx:93
keylayer_zslicemap.se
se
Definition: keylayer_zslicemap.py:194
JetTiledMap::S
@ S
Definition: TiledEtaPhiMap.h:44
empty
bool empty(TH1 *h)
Definition: computils.cxx:294
InDet::SiDetElementsRoadMaker_xk::m_pix
StringProperty m_pix
Definition: SiDetElementsRoadMaker_xk.h:154
PrepRawData.h
ReadCondHandle.h
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
InDet::SiDetElementRoadMakerData_xk
Definition: SiDetElementRoadMakerData_xk.h:28
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:93
InDet::SiDetElementsRoadMaker_xk::bookUsageTracker
static void bookUsageTracker(InDet::SiDetElementRoadMakerData_xk &data, const SiDetElementsLayerVectors_xk &layers)
this method is used to initialize the detector element usage tracker member of the event data struct ...
Definition: SiDetElementsRoadMaker_xk.cxx:521
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition: MagneticFieldMode.h:20
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
InDet::SiDetElementsRoadMaker_xk::m_layerVecKey
SG::ReadCondHandleKey< SiDetElementsLayerVectors_xk > m_layerVecKey
Created by SiDetElementsRoadCondAlg_xk.
Definition: SiDetElementsRoadMaker_xk.h:141
Amg::z
@ z
Definition: GeoPrimitives.h:36
Trk::PropDirection
PropDirection
Definition: PropDirection.h:19
InDet::compDetElements_ZRA
Definition: SiDetElementsComparison.h:62
InDet::SiDetElementsRoadMaker_xk::m_sct
StringProperty m_sct
Definition: SiDetElementsRoadMaker_xk.h:155
AtlasFieldCacheCondObj::getInitializedCache
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d (solenoid) or 3-d (toroid) position.
Definition: AtlasFieldCacheCondObj.h:52
InDet::SiDetElementsRoadMaker_xk::m_bounds
Trk::CylinderBounds m_bounds
Definition: SiDetElementsRoadMaker_xk.h:162
PixelDetectorManager.h
InDet::SiDetElementsRoadMaker_xk::m_fieldCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
Definition: SiDetElementsRoadMaker_xk.h:144
lumiFormat.i
int i
Definition: lumiFormat.py:92
z
#define z
G
#define G(x, y, z)
Definition: MD5.cxx:113
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Amg::x
@ x
Definition: GeoPrimitives.h:34
Trk::CylinderBounds
Definition: CylinderBounds.h:46
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
Trk::pion
@ pion
Definition: ParticleHypothesis.h:29
PixelAthClusterMonAlgCfg.zmax
zmax
Definition: PixelAthClusterMonAlgCfg.py:176
InDetDD::PixelDetectorManager::getDetectorElementBegin
virtual SiDetectorElementCollection::const_iterator getDetectorElementBegin() const override
Definition: PixelDetectorManager.cxx:110
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Preparation.mode
mode
Definition: Preparation.py:95
InDet::SiDetElementsRoadMaker_xk::m_outputlevel
int m_outputlevel
Definition: SiDetElementsRoadMaker_xk.h:164
InDet::SiDetElementsRoadMaker_xk::SiDetElementsRoadMaker_xk
SiDetElementsRoadMaker_xk(const std::string &, const std::string &, const IInterface *)
Definition: SiDetElementsRoadMaker_xk.cxx:35
Trk::ParametersBase
Definition: ParametersBase.h:55
InDet::compDetElementWays
Definition: SiDetElementsComparison.h:90
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
InDet::SiDetElementsRoadUtils_xk::detElementInformation
void detElementInformation(const InDetDD::SiDetectorElement &Si, double *P)
Definition: SiDetElementsRoadUtils_xk.cxx:42
InDetDD::SolidStateDetectorElementBase::normal
virtual const Amg::Vector3D & normal() const override final
Get reconstruction local normal axes in global frame.
MagField::AtlasFieldCache::solenoidOn
bool solenoidOn() const
status of the magnets
Trk::FullField
@ FullField
Field is set to be realistic, but within a given Volume.
Definition: MagneticFieldMode.h:21
SiDetElementsComparison.h
InDet::SiDetElementsRoadMaker_xk::dumpConditions
MsgStream & dumpConditions(MsgStream &out) const
Definition: SiDetElementsRoadMaker_xk.cxx:103
Trk::NoField
@ NoField
Field is set to 0., 0., 0.,.
Definition: MagneticFieldMode.h:18
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
SiDetElementsRoadMaker_xk.h
InDetDD::PixelDetectorManager
Definition: PixelDetectorManager.h:47
ReadCellNoiseFromCoolCompare.s3
s3
Definition: ReadCellNoiseFromCoolCompare.py:380
InDetDD::SCT_DetectorManager::getDetectorElementEnd
virtual SiDetectorElementCollection::const_iterator getDetectorElementEnd() const override
Definition: SCT_DetectorManager.cxx:98
InDet::SiDetElementsLayerVectors_xk
std::vector< std::vector< InDet::SiDetElementsLayer_xk > > SiDetElementsLayerVectors_xk
Definition: SiDetElementsLayerVectors_xk.h:22
InDetDD::SiDetectorElement
Definition: SiDetectorElement.h:109
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::MagneticFieldProperties::magneticFieldMode
MagneticFieldMode magneticFieldMode() const
Returns the MagneticFieldMode as specified.
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
InDet::SiDetElementsRoadMaker_xk::computeBounds
void computeBounds()
Definition: SiDetElementsRoadMaker_xk.cxx:606
Rtt_histogram.n1
n1
Definition: Rtt_histogram.py:21
InDet::SiDetElementsRoadMaker_xk::detElementsRoad
virtual void detElementsRoad(std::deque< Amg::Vector3D > &globalPositions, std::vector< const InDetDD::SiDetectorElement * > &Road, bool testDirection, SiDetElementRoadMakerData_xk &roadMakerData, const EventContext &ctx) const override
This signature assumes you already have a list of positions along the trajectory.
Definition: SiDetElementsRoadMaker_xk.cxx:306
makeTRTBarrelCans.dy
tuple dy
Definition: makeTRTBarrelCans.py:21
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
SCT_ID
Definition: SCT_ID.h:68
InDet::SiDetElementsRoadMaker_xk::getLayers
const SiDetElementsLayerVectors_xk * getLayers(const EventContext &ctx) const
Definition: SiDetElementsRoadMaker_xk.h:181
DEBUG
#define DEBUG
Definition: page_access.h:11
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition: AtlasFieldCache.h:43
InDet::SiDetElementsRoadMaker_xk::m_width
FloatProperty m_width
Definition: SiDetElementsRoadMaker_xk.h:152
InDetDD::PixelDetectorManager::getDetectorElementEnd
virtual SiDetectorElementCollection::const_iterator getDetectorElementEnd() const override
Definition: PixelDetectorManager.cxx:115
fitman.hz
def hz
Definition: fitman.py:516
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
SCT_DetectorManager.h
InDet::SiDetElementsRoadMaker_xk::getBound
Trk::CylinderBounds getBound(MagField::AtlasFieldCache &fieldCache, const Trk::TrackParameters &) const
Definition: SiDetElementsRoadMaker_xk.cxx:775
Trk::MagneticFieldMode
MagneticFieldMode
Definition: MagneticFieldMode.h:17
InDet::compDetElements_RAZ
Definition: SiDetElementsComparison.h:37
PixelID
Definition: PixelID.h:67
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
plotBeamSpotMon.nc
int nc
Definition: plotBeamSpotMon.py:83
InDet::SiDetElementsRoadMaker_xk::m_usePIX
BooleanProperty m_usePIX
Definition: SiDetElementsRoadMaker_xk.h:150
InDet::SiDetElementsRoadMaker_xk::m_step
DoubleProperty m_step
Definition: SiDetElementsRoadMaker_xk.h:153
python.SystemOfUnits.rad
int rad
Definition: SystemOfUnits.py:111
Trk::CylinderBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
InDet::SiDetElementRoadMakerData_xk::resetUsageTracker
void resetUsageTracker()
method to reset the flags stored in the elementUsageTracker (below) when building a new search road.
Definition: SiDetElementRoadMakerData_xk.cxx:4
Trk::CylinderBounds::r
virtual double r() const override final
This method returns the radius.
InDet::SiDetElementRoadMakerData_xk::isInitialized
bool isInitialized
Flag to check if the event data was already initialized by the client tool.
Definition: SiDetElementRoadMakerData_xk.h:59
InDet::SiDetElementsRoadMaker_xk::m_fieldModeEnum
Trk::MagneticFieldMode m_fieldModeEnum
Definition: SiDetElementsRoadMaker_xk.h:163
Generated on Fri Jul 5 2024 21:26:55 for ATLAS Offline Software by   doxygen 1.8.18