74{
75
76 SG::WriteCondHandle<InDet::TRT_TrackSegmentsToolCondData_xk> writeHandle{
m_writeKey, ctx};
79 return StatusCode::SUCCESS;
80 }
81
82 EventIDRange rangeTrt;
83
84 SG::ReadCondHandle<InDetDD::TRT_DetElementContainer> trtDetEleHandle(
m_trtDetEleContKey, ctx);
85 if (not trtDetEleHandle.isValid()) {
87 return StatusCode::FAILURE;
88 }
89
90 const InDetDD::TRT_Numerology* trtNum = trtDetEleHandle->getTRTNumerology();
91 if (trtNum==nullptr){
93 }
94
95 if (not trtDetEleHandle.range(rangeTrt)) {
96 ATH_MSG_FATAL(
"Failed to retrieve validity range for " << trtDetEleHandle.key());
97 return StatusCode::FAILURE;
98 }
99
100 auto writeCdo = std::make_unique<InDet::TRT_TrackSegmentsToolCondData_xk>();
101 if(writeCdo->m_ndzdr ) delete [] writeCdo->m_ndzdr ;
102 if(writeCdo->m_slope ) delete [] writeCdo->m_slope ;
103 if(writeCdo->m_islope ) delete [] writeCdo->m_islope ;
104
105 float RZ [4][200];
108
109
110
114
115 for(int ring = 0; ring!=Rings; ++ring) {
116
118
119 writeCdo->m_flayers[1][ring] = writeCdo->m_nlayers[1]; writeCdo->m_flayers[2][ring] = writeCdo->m_nlayers[2];
120 writeCdo->m_nlayers[1] += NSlayers ; writeCdo->m_nlayers[2] += NSlayers ;
121
122 for(int nsl=0; nsl!=NSlayers; ++nsl) {
123
124
125 for(
int f=0;
f!=NPhi; ++
f) {
126
127 const InDetDD::TRT_BaseElement* base1 = trtDetEleHandle->getBarrelDetElement(0,ring,f,nsl);
128 const InDetDD::TRT_BaseElement* base2 = trtDetEleHandle->getBarrelDetElement(1,ring,f,nsl);
129 if(!base1 || !base2) continue;
130 const InDetDD::TRT_BarrelElement* bael1 =
131 dynamic_cast<const InDetDD::TRT_BarrelElement*>(base1);
132 const Trk::RectangleBounds* rb1 =
133 dynamic_cast<const Trk::RectangleBounds*
>(&base1->
bounds());
134 const Trk::RectangleBounds* rb2 =
135 dynamic_cast<const Trk::RectangleBounds*
>(&base2->
bounds());
136 if(!bael1 || !rb1 || !rb2) continue;
137 float rmean = 0;
138 if(f==0) {
139
142 RZ [1][
n] = std::sqrt(C1.x()*C1.x()+C1.y()*C1.y());
143
146
147 RZ [2][
n] = std::sqrt(C2.x()*C2.x()+C2.y()*C2.y());
148
152 }
154 for(
int s=0;
s!=
ns; ++
s) {
155
156 const Identifier id1 =
m_trtid->straw_id(-1,f,ring,nsl,s );
160 if(f==0) rmean+=std::sqrt(sc1->x()*sc1->x()+sc1->y()*sc1->y());
161
162 if(!sc1 || !st1 || !tr1 ) {
ATH_MSG_ERROR(
"problem with TRT geometry");}
163 ++writeCdo->m_nstraws[1];
164 const Identifier
id2 =
m_trtid->straw_id(+1,f,ring,nsl,s );
168 if(!sc2 || !st2 || !tr2) {
ATH_MSG_ERROR(
"problem with TRT geometry");}
169 ++writeCdo->m_nstraws[2];
170 }
171 if(f==0) { RZ[1][
n-1] = RZ[2][
n-1] = rmean/
float(ns);}
172 }
173 }
174 }
175
176
177
178 int Wheels = trtNum->
getNEndcapWheels();
if(!Wheels)
return StatusCode::SUCCESS;
181 for(
int wh = 0;
wh!=Wheels; ++
wh) {
182
184 writeCdo->m_flayers[0][
wh] = writeCdo->m_nlayers[0]; writeCdo->m_flayers[3][
wh] = writeCdo->m_nlayers[3];
185 writeCdo->m_nlayers[0] += NSlayers ; writeCdo->m_nlayers[3] += NSlayers ;
186
187 for(int nsl = 0; nsl!=NSlayers; ++nsl) {
188
189 for(
int f=0;
f!=NPhi; ++
f) {
190
191 const InDetDD::TRT_BaseElement* base1 = trtDetEleHandle->getEndcapDetElement(0,wh,nsl,f);
192 const InDetDD::TRT_BaseElement* base2 = trtDetEleHandle->getEndcapDetElement(1,wh,nsl,f);
193 if(!base1 || !base2) continue;
194
195 const InDetDD::TRT_EndcapElement* enel1 =
196 dynamic_cast<const InDetDD::TRT_EndcapElement*>(base1);
197 const Trk::DiscBounds*
db1 =
198 dynamic_cast<const Trk::DiscBounds*
>(&base1->
bounds());
199 const Trk::DiscBounds*
db2 =
200 dynamic_cast<const Trk::DiscBounds*
>(&base2->
bounds());
201 if(!enel1 || !db1 || !db2) continue;
202 if(f==0) {
203
209
213
215 }
216
218
219 for(
int s=0;
s!=
ns; ++
s) {
220
221 const Identifier id1 =
m_trtid->straw_id(-2,f,wh,nsl,s );
225
226 if(!sc1 || !st1 || !tr1) {
ATH_MSG_ERROR(
"problem with TRT geometry");}
227 ++writeCdo->m_nstraws[0];
228
229 const Identifier
id2 =
m_trtid->straw_id(+2,f,wh,nsl,s );
233
234 if(!sc2 || !st2 || !tr2) {
ATH_MSG_ERROR(
"problem with TRT geometry");}
235 ++writeCdo->m_nstraws[3];
236 }
237 }
238 }
239 }
240
241 double zmax = RZ[3][writeCdo->m_nlayers[3]-1]+10.;
242 double rmax = RZ[2][writeCdo->m_nlayers[2]-1]+10.;
243 const Trk::CylinderBounds
CB(rmax,zmax);
244
245 float rapidity[26]=
246 {-2.05,-1.95,-1.84,-1.72,-1.62,-1.53,-1.43,-1.33,-1.21,-1.00,-.94, -.85,-.32,
247 .32, .85, .94, 1.00, 1.21, 1.33, 1.43, 1.53, 1.62, 1.72,1.84, 1.95,2.05};
248
249 std::deque<Amg::Vector3D>
G [26];
251 Trk::PerigeeSurface
ps(psv);
252
253 for(
int r=0;
r!=26; ++
r) {
254 writeCdo->m_dzdr[
r] = 1./std::tan(2.*std::atan(std::exp(-rapidity[
r])));
255 double pinv =-1./(
m_pTmin*std::sqrt(1.+writeCdo->m_dzdr[
r]*writeCdo->m_dzdr[
r]));
256 auto trackPar =
ps.createUniqueTrackParameters(
257 0.,
258 0.,
259 0.,
260 std::atan2(1.,
double(writeCdo->m_dzdr[
r])),
261 pinv,
262 std::nullopt);
265 }
266
268 for(
int b=0;
b!=4; ++
b) {
269 for(
unsigned int i=0;
i!=writeCdo->m_nlayers[
b]; ++
i) {
270 writeCdo->m_begin[
b][
i] =
n;
271 for(
float r : writeCdo->m_dzdr) {
272 if(
r >= Tmin[b][i] &&
r <= Tmax[b][i] ) {
273 writeCdo->m_end[
b][
i] =
n++;
274 }
275 }
276 }
277 }
278
279 writeCdo->m_ndzdr =
new unsigned int[
n];
280 writeCdo->m_islope =
new int[
n];
281 writeCdo->m_slope =
new float [
n];
282 writeCdo->m_cirsize = writeCdo->m_nstraws[0]+writeCdo->m_nstraws[1];
283
285 for (
int b = 0;
b != 4; ++
b) {
286
287 for (
unsigned int i = 0;
i != writeCdo->m_nlayers[
b]; ++
i) {
288
289 for (
int r = 0;
r != 26; ++
r) {
290 if (writeCdo->m_dzdr[
r] >= Tmin[b][i] &&
291 writeCdo->m_dzdr[
r] <= Tmax[b][i]) {
292 writeCdo->m_ndzdr[
n] =
r;
293 std::deque<Amg::Vector3D>::iterator gp0, gp1, gp =
G[
r].begin(),gpe =
G[
r].end();
294 if (b == 0 || b == 3) {
295
296 gp0 = gp;
297 for (++gp; gp != gpe; ++gp) {
298 if (b == 3 && (*gp).z() >= RZ[b][i])
299 break;
300 if (b == 0 && (*gp).z() <= RZ[b][i])
301 break;
302 gp1 = gp0;
303 gp0 = gp;
304 }
305 } else {
306 gp0 = gp;
307 for (++gp; gp != gpe; ++gp) {
308 if (std::sqrt((*gp).x() * (*gp).x() + (*gp).y() * (*gp).y()) >
309 RZ[b][i])
310 break;
311 gp1 = gp0;
312 gp0 = gp;
313 }
314 }
316 if (gp != gpe) {
323 double as = 1. / std::sqrt(ax * ax + ay * ay + az * az);
326 az *= as;
327 } else {
334 double as = 1. / std::sqrt(ax * ax + ay * ay + az * az);
337 az *= as;
338 }
340
341 if (b == 0 || b == 3) {
342 S = (RZ[
b][
i] -
z) / az;
343 } else {
344 double A = (
ax *
x +
ay *
y) * 2.;
345 double D = (RZ[
b][
i] -
x -
y) * (RZ[b][i] +
x +
y) + 2. *
x *
y;
348 double Sq =
A *
A + 2. * D *
B;
349 Sq > 0. ? Sq = std::sqrt(Sq) : Sq = 0.;
350 double S1 = -(
A + Sq) / B;
351 double S2 = -(
A - Sq) / B;
356 }
357 writeCdo->m_slope[
n] = std::atan2(
y + S * ay,
x + S * ax) *
m_A;
358 writeCdo->m_islope[
n] =
int(writeCdo->m_slope[n] *
m_Psi128);
360 }
361 }
362 }
363 }
364
365 if (writeHandle.
record(rangeTrt, std::move(writeCdo)).isFailure()) {
367 << " with EventRange " << rangeTrt
368 << " into Conditions Store");
369 return StatusCode::FAILURE;
370 }
372 << rangeTrt << " into Conditions Store");
373
374 return StatusCode::SUCCESS;
375
376}
struct TBPatternUnitContext S2
struct TBPatternUnitContext S1
unsigned int nStraws() const
Number of straws in the element.
virtual const Trk::SurfaceBounds & bounds() const override final
Straw layer bounds.
virtual const Trk::Surface & surface() const override final
Element Surface: access to the Surface (straw layer)
virtual const Amg::Vector3D & center() const override final
Element Surface: center of a straw layer.
const Amg::Transform3D & strawTransform(unsigned int straw) const
Straw transform - fast access in array, in Tracking frame: Amg.
const Amg::Vector3D & strawCenter(int straw) const
Straw Surface: Local -> global transform of the straw via integer.
unsigned int getNEndcapWheels() const
unsigned int getNEndcapPhi() const
unsigned int getNBarrelPhi() const
unsigned int getNBarrelRings() const
unsigned int getNBarrelLayers(unsigned int iMod) const
unsigned int getNEndcapLayers(unsigned int iWheel) const
Trk::MagneticFieldProperties m_fieldprop
ToolHandle< Trk::IPropagator > m_propTool
SG::WriteCondHandleKey< TRT_TrackSegmentsToolCondData_xk > m_writeKey
SG::ReadCondHandleKey< InDetDD::TRT_DetElementContainer > m_trtDetEleContKey
const std::string & key() const
StatusCode record(const EventIDRange &range, T *t)
record handle, with explicit range DEPRECATED
const DataObjID & fullKey() const
double halflengthY() const
for consitant naming
Eigen::Affine3d Transform3D
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Eigen::Matrix< double, 3, 1 > Vector3D