120 {
121
122 ATH_MSG_DEBUG(
"Transforming the TRT segment into a track...");
123
124
125
126
129 return nullptr;
130 }
132
133
134
135
138 const Trk::TrackParameters* segPar =
139 tS.associatedSurface()
140 .createUniqueTrackParameters(par[Trk::loc1],
141 par[Trk::loc2],
144 par[Trk::qOverP],
145 std::move(ep))
146 .release();
147
150 } else {
151 ATH_MSG_DEBUG(
"Could not get initial TRT segment parameters! ");
152
153 return nullptr;
154 }
155
156
157 auto ntsos = std::make_unique<Trk::TrackStates>();
158
159
160
161
163
164 const Trk::TrackStateOnSurface* par_tsos = nullptr;
165
166
168 Trk::PerigeeSurface perigeeSurface(perigeePosition);
169
170 std::unique_ptr<Trk::TrackParameters>
tmp =
172 std::unique_ptr<Trk::Perigee> perParm = nullptr;
173
176 }
177 if (perParm) {
179 } else {
180 ATH_MSG_DEBUG(
"Failed to build perigee parameters.Discard...");
181 ntsos->clear();
182 delete segPar;
183 segPar = nullptr;
184 return nullptr;
185 }
186
187
188 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
189 typePattern;
191 par_tsos = new Trk::TrackStateOnSurface(
192 nullptr, std::move(perParm), nullptr, typePattern);
193
194 ntsos->push_back(par_tsos);
195 }
196
197
198
199
200
201
202 int npseudo = 0;
203 double pseudotheta = 0;
204 const Trk::MeasurementBase* pseudo = nullptr;
205
206 const Trk::Surface *firstsurf = nullptr, *firstecsurf = nullptr,
207 *lastsurf = nullptr;
208 const Trk::MeasurementBase* firstmeas = nullptr;
209
210 int nbarrel = 0, nendcap = 0;
211
212 std::vector<std::pair<double, double>> points;
213 points.reserve(40);
214 double oldphi = 0;
215
216 const Trk::Surface *surf_zmin=nullptr;
217 const Trk::Surface *surf_zmax=nullptr;
220
221
223
224
225 const Trk::TrackStateOnSurface* seg_tsos = nullptr;
226
227
228 if (dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(
230
231 if (pseudo) {
232
233 pseudotheta =
236 }
237
239
240 npseudo++;
241
243
244 seg_tsos =
246 }
247
248 } else {
249
250
251
252
253 seg_tsos =
255
256
257
258
259
260 if (!firstmeas)
262 if (!firstsurf)
264
266
267
269
270 nendcap++;
271 if (!firstecsurf)
273
274 double tmpphi =
276 if (!points.empty() &&
277 std::abs(tmpphi - oldphi) >
M_PI) {
278 if (tmpphi < 0)
280 else
282 }
283
284 oldphi = tmpphi;
285
289 }
293 }
294
295 points.emplace_back(
297
298 } else
299 nbarrel++;
300
301
302
303
304 }
305
306
307 if (seg_tsos)
308 ntsos->push_back(seg_tsos);
309 }
310
311
314
315
317 return new Trk::Track(info, std::move(ntsos), std::move(fq));
318 } else {
319
320
321
322
323
324
325 if (npseudo == 1)
327
328
329 double myqoverp = 0, myphi = 0, myd0 = 0, myz0 = 0, mytheta = pseudotheta;
330
331 if (nendcap < 4) {
332
333
334
335
336
337 myqoverp =
par[4] * std::sin(pseudotheta) / std::sin(par[3]);
338
339
341 Trk::PerigeeSurface perigeeSurface(perigeePosition);
342
343 std::unique_ptr<const Trk::TrackParameters>
tmp =
345 std::unique_ptr<const Trk::Perigee> tempper = nullptr;
348 }
349 if (!tempper) {
351 delete segPar;
352 segPar = nullptr;
353 return nullptr;
354 }
355
356
357 myd0 = tempper->parameters()[
Trk::d0];
358 myphi = tempper->parameters()[
Trk::phi0];
359
360 } else {
361
362
363
364
365
366 double sx = 0,
sy = 0, sxx = 0, sxy = 0,
d = 0;
367 float zmin = points.empty() ? 0 : points.begin()->first,
zmax = 0;
368
369
370 for (auto & point : points) {
373 sxy += point.first * point.second;
374 sxx += point.first * point.first;
375 if (std::abs(point.first) > std::abs(zmax)) {
377 }
378 if (std::abs(point.first) < std::abs(zmin)) {
380 }
381 }
382
383 if (std::abs(pseudotheta) < 1.e-6) {
385
390 double diff_z (meas_zmax.z() - meas_zmin.z());
391 if (std::abs(diff_z)>1e-6) {
392 double diff_r( meas_zmax.perp() - meas_zmin.perp());
393 pseudotheta = std::atan2(diff_r, diff_z);
394 ATH_MSG_DEBUG(
"Initial pseudo theta is zero. Compute pseudo from inner- and outermost endcap measurements. delta R: " <<
395 meas_zmax.perp() << " - " << meas_zmin.perp() << " = " << diff_r
396 << " , diff Z: " << meas_zmax.z() << " - " << meas_zmin.z()
397 << " = " << diff_z
398 << " " << pseudotheta);
399 }
400 }
401
402 if (std::abs(pseudotheta) < 1.e-6) {
403 constexpr std::array<float,2> boundary_r {644., 1004.};
404
405 bool is_across_sides=std::signbit(zmin)!=std::signbit(zmax);
406 double r_path=is_across_sides ? 2 * boundary_r[1] : boundary_r[1]-boundary_r[0];
409 if (surf_zmin && surf_zmax) {
410 bool is_reverse=std::abs(surf_zmin->
center().z())>std::abs(surf_zmax->
center().z());
411
412
413
414
415
416
417 std::array<const Trk::Surface *,2> boundary_surface{
418 is_reverse ? surf_zmax : surf_zmin,
419 is_reverse ? surf_zmin : surf_zmax};
420
421 std::array<Amg::Vector2D,2> translation;
422 std::array<Amg::Vector2D,2> height;
423 for (unsigned int boundary_i=boundary_surface.size(); boundary_i-->0;) {
426 const Trk::CylinderBounds &
427 cylinder_bounds = dynamic_cast<const Trk::CylinderBounds &>(boundary_surface[boundary_i]->bounds());
428 height[boundary_i]= project2D( boundary_surface[boundary_i]->
transform().
rotation()
430 translation[boundary_i]=project2D( boundary_surface[boundary_i]->
transform().translation() );
431 }
432 }
433 r1 = (translation[1] + height[1] - ( translation[0] + height[0])).
perp();
434 r2 = (translation[1] + height[1] - ( translation[0] - height[0])).
perp();
435 r_path=std::max(r1,r2);
436
437 }
438 pseudotheta = std::atan2(r_path, zmax - zmin);
439 ATH_MSG_DEBUG(
"Initial pseudo theta is zero. Pseudo theta from inner- and outermost surfaces. Deleta r " << r2 <<
" - " << r1 <<
" -> " << r_path
440 << " , delta Z " << zmax << " - " << zmin << " -> " << (zmax-zmin)
441 << " " << pseudotheta);
442 }
443 }
444
445
446 d = (points.size() * sxx - sx * sx);
447 double dphidz = ((points.size() * sxy - sy * sx) / d);
448 myqoverp = (std::abs(pseudotheta) > 1
e-6)
449 ? -dphidz / (0.6 * std::tan(pseudotheta))
450 : 1000.;
451
452
453 double halfz = 200.;
454 const Trk::CylinderBounds* cylb =
455 dynamic_cast<const Trk::CylinderBounds*
>(&firstsurf->
bounds());
456 if (!cylb)
458 else
460 const Trk::CylinderBounds* cylb2 =
461 dynamic_cast<const Trk::CylinderBounds*>(&lastsurf->bounds());
462 double halfz2 = 200.;
463 if (!cylb2)
465 else
468 Amg::Vector3D strawdir2 = -lastsurf->transform().rotation().col(2);
471
472
473 if (std::abs(lastsurf->center().z()) < 2650 * mm) {
474 pos2 = lastsurf->center() + halfz2 * strawdir2;
475 if (nbarrel == 0) {
476 double dr = std::abs(std::tan(pseudotheta) * (lastsurf->center().z() -
477 firstsurf->
center().z()));
478 pos1 = firstsurf->
center() + (halfz -
dr) * strawdir1;
479 } else {
480 double dz = std::abs((pos2.perp() - firstsurf->
center().perp()) /
481 std::tan(pseudotheta));
482 if (pos2.z() > 0)
483 dz = -dz;
484 double z1 = pos2.z() + dz;
487 }
488 } else {
489 double dr = std::abs(std::tan(pseudotheta) *
490 (lastsurf->center().z() - firstsurf->
center().z()));
491 pos2 = lastsurf->center() + (
dr - halfz2) * strawdir2;
492 pos1 = firstsurf->
center() - halfz * strawdir1;
493 }
494
495
496
497 if (nbarrel == 0 &&
498 std::abs(std::tan(pseudotheta) * (firstsurf->
center().z() -
499 lastsurf->center().z())) < 250 * mm &&
500 std::abs(firstsurf->
center().z()) > 1000 * mm) {
501
502
503 const Trk::MeasurementBase* firstmeas =
504 (**ntsos->begin()).measurementOnTrack();
506 newcov.setZero();
508 newcov(1, 1) = (myqoverp != 0) ? .0001 * myqoverp * myqoverp : 1.e-8;
509 Trk::LocalParameters newpar(std::make_pair(0,
Trk::locZ),
511 auto newpseudo =
512 std::make_unique<Trk::PseudoMeasurementOnTrack>(
513 std::move(newpar),
514 std::move(newcov),
516
517 ntsos->erase(ntsos->begin());
518 ntsos->insert(ntsos->begin(),
519 new Trk::TrackStateOnSurface(std::move(newpseudo), nullptr));
520 }
521
523
524 MagField::AtlasFieldCache fieldCache;
525
526
527 SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{
529 };
530 const AtlasFieldCacheCondObj* fieldCondObj{ *readHandle };
531 if (fieldCondObj == nullptr) {
533 "segToTrack: Failed to retrieve AtlasFieldCacheCondObj with key "
535 return nullptr;
536 }
538
539
541 field1.data());
542
544
545 double phideflection =
546 -.3 * (pos2 - pos1).
perp() * field1.z() * myqoverp / std::sin(pseudotheta);
547 double precisephi = (nbarrel == 0)
548 ? (pos2 - pos1).phi() - .5 * phideflection
549 : (pos2 - pos1).
phi() + .5 * phideflection;
550 double radius = (myqoverp != 0. && field1.z() != 0.)
551 ? -std::sin(pseudotheta) / (.3 * field1.z() * myqoverp)
552 : 1.e6;
553 double precisetheta =
554 (myqoverp != 0.)
555 ? std::atan2(std::abs(radius * phideflection), pos2.z() - pos1.z())
556 : pseudotheta;
557 if (precisetheta < 0)
558 precisetheta +=
M_PI;
559 if (precisephi < -
M_PI)
560 precisephi += 2 *
M_PI;
561 if (precisephi >
M_PI)
562 precisephi -= 2 *
M_PI;
563
564
565
566 const Trk::StraightLineSurface* surfforpar = nullptr;
567 if (nbarrel == 0)
568 surfforpar = dynamic_cast<const Trk::StraightLineSurface*>(firstsurf);
569 else
570 surfforpar = dynamic_cast<const Trk::StraightLineSurface*>(lastsurf);
571 if (!surfforpar)
572 ATH_MSG_ERROR(
"Cast of surface failed, should never happen");
573
575 precisephi,
576 precisetheta,
577 myqoverp,
578 *surfforpar);
579 Trk::PerigeeSurface persurf;
581 m_extrapolator->extrapolateDirectly(ctx, ataline, persurf).release();
582
583
584 if (extrappar) {
585 if (nendcap >= 4) {
586 myphi = extrappar->parameters()[
Trk::phi0];
587 myd0 = extrappar->parameters()[
Trk::d0];
588 }
589
590
591 double z0 = extrappar->parameters()[
Trk::z0];
592 if (nbarrel == 0)
593 mytheta = std::atan(std::tan(extrappar->parameters()[
Trk::theta]) *
594 std::abs((z0 - pos1.z()) / pos1.z()));
595 else
596 mytheta = std::atan(std::tan(extrappar->parameters()[
Trk::theta]) *
597 std::abs((z0 - pos2.z()) / pos2.z()));
598
599 if (mytheta < 0)
601
602 delete extrappar;
603 extrappar = nullptr;
604 }
605 }
608 while (myphi < -
M_PI)
610
611 double P[5] = { myd0, myz0, myphi, mytheta, myqoverp };
612
613
614
615 auto per =
616 std::make_unique<Trk::Perigee>(P[0], P[1], P[2], P[3], P[4], Trk::PerigeeSurface());
617 std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
618 typePattern;
620 Trk::TrackStateOnSurface* seg_tsos = new Trk::TrackStateOnSurface(
621 nullptr, std::move(per), nullptr, typePattern);
622 ntsos->insert(ntsos->begin(), seg_tsos);
623
624
625
626
627
628
629 Trk::Track newTrack (info, std::move(ntsos), std::move(fq));
630 Trk::Track* fitTrack =
632
633
634 if (segPar) {
635 delete segPar;
636 segPar = nullptr;
637 }
638
639 if (!fitTrack) {
641 return nullptr;
642 }
643
644
645
646
647
651 do {
652
653 if ((*parit)->covariance() &&
655 firstmeaspar = *parit;
656 ++parit;
657 } while (firstmeaspar == nullptr &&
659
660
661
662
663
665
666 if (!perTrack || !perTrack->covariance()) {
667 ATH_MSG_ERROR(
"Cast of perigee fails, should never happen !");
668 return nullptr;
669 } else {
670 ATH_MSG_VERBOSE (
"Perigee after refit with fudges to make it converge : " << (*perTrack) );
671
672 if(firstmeaspar && firstmeaspar->
position().perp()<2000*mm && std::abs(firstmeaspar->
position().z())<3000*mm){
673
674
676
677 double scaleZ = std::sqrt(tS.localCovariance()(1,1))/std::sqrt( (fcovmat)(1,1));
678 double scaleTheta = std::sqrt(tS.localCovariance()(3,3))/std::sqrt( (fcovmat)(3,3));
679
680 fcovmat(1,0)=scaleZ*((fcovmat)(1,0));
681 fcovmat(0,1) = (fcovmat)(1,0);
682 fcovmat(1,1)=tS.localCovariance()(1,1);
683 fcovmat(2,1)=scaleZ*((fcovmat)(2,1));
684 fcovmat(1,2) = (fcovmat)(2,1);
685 fcovmat(3,1)=scaleZ*scaleTheta*((fcovmat)(3,1));
686 fcovmat(1,3) = (fcovmat)(3,1);
687 fcovmat(4,1)=scaleZ*((fcovmat)(4,1));
688 fcovmat(1,4) = (fcovmat)(4,1);
689 fcovmat(3,0)=scaleTheta*((fcovmat)(3,0));
690 fcovmat(0,3) = (fcovmat)(3,0);
691 fcovmat(3,2)=scaleTheta*((fcovmat)(3,2));
692 fcovmat(2,3) = (fcovmat)(3,2);
693 fcovmat(3,3)=tS.localCovariance()(3,3);
694 fcovmat(4,3)=scaleTheta*((fcovmat)(4,3));
695 fcovmat(3,4) = (fcovmat)(4,3);
696
697
698 const
AmgVector(5)& par = firstmeaspar->parameters();
699 const Trk::TrackParameters* updatedPars =
700 firstmeaspar->associatedSurface().createUniqueTrackParameters(
701 par[Trk::loc1],
702 par[Trk::loc2],
705 par[Trk::qOverP],
706 std::move(fcovmat)).release();
707
708
709 const Trk::TrackParameters* newperpar =
711 *updatedPars,
712 perTrack->associatedSurface(),
713 Trk::anyDirection,
714 false,
715 Trk::nonInteracting).release();
716 delete updatedPars; updatedPars = nullptr;
717
718 if (!newperpar || !newperpar->covariance()) {
719 ATH_MSG_WARNING (
"Can not hack perigee parameters, extrapolation failed");
720 delete newperpar; newperpar = nullptr;
721 } else {
722
723
724
726
727 errmat = *newperpar->covariance();
728 delete newperpar; newperpar = nullptr;
729
731 if( CM(1,1)==0.||CM(3,3)==0. ) {
732 ATH_MSG_DEBUG (
"Hacked perigee covariance is CRAP, reject track");
733 delete fitTrack; return nullptr;
734 } else {
735 ATH_MSG_VERBOSE (
"Perigee after fit with scaled covariance matrix : " << *perTrack);
736 }
737 }
738 }
739 }
740
741 return fitTrack;
742 }
743 }
Scalar perp() const
perp method - perpendicular length
Scalar phi() const
phi method
Scalar theta() const
theta method
#define ATH_MSG_VERBOSE(x)
#define AmgSymMatrix(dim)
char data[hepevt_bytes_allocation_ATLAS]
#define ATLAS_THREAD_SAFE
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
DataModel_detail::const_iterator< DataVector > const_iterator
Standard const_iterator.
void getField(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field value at given position xyz[3] is in mm, bxyz[3] is in kT if deriv[9] is given,...
double halflengthZ() const
This method returns the halflengthZ.
std::unique_ptr< FitQuality > uniqueClone() const
NVI uniqueClone.
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
virtual const Surface & associatedSurface() const =0
Interface method to get the associated Surface.
virtual const Amg::Vector3D & globalPosition() const =0
Interface method to get the global Position.
const Amg::MatrixX & localCovariance() const
Interface method to get the localError.
std::unique_ptr< MeasurementBase > uniqueClone() const
NVI Clone giving up unique pointer.
const Amg::Vector3D & position() const
Access method for the position.
const FitQuality * fitQuality() const
return the FitQuality object, returns NULL if no FitQuality is defined
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Identifier associatedDetectorElementIdentifier() const
return Identifier of the associated Detector Element
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
const Amg::Vector3D & center() const
Returns the center position of the Surface.
@ TRTStandalone
TRT Standalone.
const Surface & associatedSurface() const override final
returns the surface for the local to global transformation
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
const DataVector< const TrackParameters > * trackParameters() const
Return a pointer to a vector of TrackParameters.
const Perigee * perigeeParameters() const
return Perigee.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Eigen::Matrix< double, 3, 1 > Vector3D
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
ParametersT< TrackParametersDim, Charged, StraightLineSurface > AtaStraightLine
ParametersBase< TrackParametersDim, Charged > TrackParameters
1.14.0