From 9656e128b280a413ecf33c40c8de438fbcc11b0f Mon Sep 17 00:00:00 2001
From: scott snyder <sss@karma>
Date: Thu, 21 Feb 2019 15:47:41 +0100
Subject: [PATCH] TrkExSolenoidalIntersector: Preparing to make interfaces
const.
The tool was caching event information in member variables.
Rework to remove some of this information and assocciate the rest
with the TrackSurfaceIntersection object.
Results are identical, as per RunTier0Tests.
Also changed CRLF line endings to LF.
---
.../SolenoidParametrization.h | 262 ++---
.../SolenoidalIntersector.h | 563 ++++-----
.../share/SolenoidalIntersector_test.ref | 5 +-
.../src/SolenoidParametrization.cxx | 215 ++--
.../src/SolenoidalIntersector.cxx | 1005 +++++++++--------
.../test/SolenoidParametrization_test.cxx | 14 +-
.../test/SolenoidalIntersector_test.cxx | 2 -
7 files changed, 1074 insertions(+), 992 deletions(-)
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
index 67dd6513807..8340f1a2016 100755
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
/***************************************************************************
@@ -33,44 +33,75 @@ class SolenoidParametrization
{
public:
- // private constructors (as singleton)
- ~SolenoidParametrization (void);
+ class BinParameters
+ {
+ public:
+ BinParameters (const double zAtAxis, const double cotTheta)
+ : m_cotTheta (cotTheta), m_zAtAxis (zAtAxis) {}
+ BinParameters (const double r, const double z, const double cotTheta);
+ double m_signTheta = 1;
+ double m_cotTheta = 0;
+ double m_zAtAxis = 0;
+ double m_interpolateZ = 0;
+ double m_complementZ = 0;
+ double m_interpolateTheta = 0;
+ double m_complementTheta = 0;
+ };
+ class Parameters
+ : public BinParameters
+ {
+ public:
+ Parameters (const SolenoidParametrization& spar,
+ const double r, const double z, const double cotTheta);
+ double m_fieldAtOrigin;
+ double m_quadraticTerm;
+ double m_cubicTerm;
+ };
+
+ SolenoidParametrization (MagField::IMagFieldSvc* magFieldSvc); // configuration from Intersector
+ ~SolenoidParametrization() = default;
// forbidden copy constructor
// forbidden assignment operator
- static SolenoidParametrization* instance (MagField::IMagFieldSvc* magFieldSvc); // initialize singleton
- static void clearInstance(void); // clear (to be used before reinitialization)
-
- double centralField (void) const;
- double fieldComponent (double) const; // parametrized - perp to track in rz-plane
+ double centralField() const;
+ double fieldComponent (double z,
+ const Parameters& parms) const; // parametrized - perp to track in rz-plane
double fieldComponent (double r, double z, double cotTheta) const; // from MagFieldSvc
- void fieldIntegrals (double&, double&, double, double);
- bool hasParametrization (void) const;
- double maximumR (void) const; // param valid to maximumR
- double maximumZ (void) const; // param valid to maximumZ
- void printFieldIntegrals (void) const;
- void printParametersForEtaLine (double, double);
- void printResidualForEtaLine (double, double);
- void setParameters (double, double, double); // set line in rz-plane
- bool validOrigin(Amg::Vector3D origin) const; // param valid for this origin?
+ void fieldIntegrals (double& firstIntegral,
+ double& secondIntegral,
+ double zBegin,
+ double zEnd,
+ Parameters& parms) const;
+ double maximumR() const; // param valid to maximumR
+ double maximumZ() const; // param valid to maximumZ
+ void printFieldIntegrals() const;
+ void printParametersForEtaLine(double eta, double z_origin) const;
+ void printResidualForEtaLine (double eta, double zOrigin) const;
+ bool validOrigin(const Amg::Vector3D& origin) const; // param valid for this origin?
+
+ // OK to use this parametrization for CURRENT?
+ bool currentMatches (double current) const;
private:
- // implementation does not work without magnetic field service
- // SolenoidParametrization (void); // singleton constructor
- SolenoidParametrization (MagField::IMagFieldSvc* magFieldSvc); // configuration from Intersector
- int fieldKey(void);
- void integrate(double&, double&, double, double) const;
- double interpolate(int, int, int, int) const;
- void parametrizeSolenoid (void);
- void setTerms(int);
+ friend class Parameters;
+
+ static int fieldKey(BinParameters& parms);
+ void integrate(double& firstIntegral,
+ double& secondIntegral,
+ double zBegin,
+ double zEnd,
+ const Parameters& parms) const;
+ double interpolate(int key1, int key2,
+ int key3, int key4,
+ const Parameters& parms) const;
+ void parametrizeSolenoid();
+ void setTerms(int, Parameters& parms) const;
- static SolenoidParametrization* s_instance;
static const double s_binInvSizeTheta;
static const double s_binInvSizeZ;
static const double s_binZeroTheta;
static const double s_binZeroZ;
- static double s_centralField;
static const double s_lightSpeed;
static const int s_maxBinTheta;
static const int s_maxBinZ;
@@ -81,60 +112,52 @@ private:
static const double s_rOuter;
static const double s_zInner;
static const double s_zOuter;
- static double s_parameters[14688];
- double m_complementTheta;
- double m_complementZ;
- double m_cotTheta;
- double m_cubicTerm;
- double m_fieldAtOrigin;
- bool m_hasParametrization;
- double m_interpolateTheta;
- double m_interpolateZ;
MagField::IMagFieldSvc* m_magFieldSvc;
- double m_quadraticTerm;
- double m_signTheta;
- double m_zAtAxis;
+ double m_currentMin;
+ double m_currentMax;
+ double m_centralField;
+ double m_parameters[14688];
// copy, assignment: no semantics, no implementation
- SolenoidParametrization (const SolenoidParametrization&);
- SolenoidParametrization &operator= (const SolenoidParametrization&);
+ SolenoidParametrization (const SolenoidParametrization&) = delete;
+ SolenoidParametrization &operator= (const SolenoidParametrization&) = delete;
};
//<<<<<< INLINE PRIVATE MEMBER FUNCTIONS >>>>>>
inline int
-SolenoidParametrization::fieldKey(void)
+SolenoidParametrization::fieldKey(BinParameters& parms)
{
- double z = m_zAtAxis - s_binZeroZ;
+ double z = parms.m_zAtAxis - s_binZeroZ;
int zBin = static_cast<int>(s_binInvSizeZ*z);
- m_interpolateZ = z*s_binInvSizeZ - double(zBin);
+ parms.m_interpolateZ = z*s_binInvSizeZ - double(zBin);
if (zBin < 0)
{
- m_interpolateZ = 0.;
+ parms.m_interpolateZ = 0.;
zBin = 0;
}
else if (zBin > s_maxBinZ - 2)
{
- m_interpolateZ = 0.;
+ parms.m_interpolateZ = 0.;
zBin = s_maxBinZ - 1;
}
- m_complementZ = 1. - m_interpolateZ;
+ parms.m_complementZ = 1. - parms.m_interpolateZ;
- int thetaBin = static_cast<int>(s_binInvSizeTheta*m_cotTheta);
- m_interpolateTheta = m_cotTheta*s_binInvSizeTheta - double(thetaBin);
+ int thetaBin = static_cast<int>(s_binInvSizeTheta*parms.m_cotTheta);
+ parms.m_interpolateTheta = parms.m_cotTheta*s_binInvSizeTheta - double(thetaBin);
if (thetaBin > s_maxBinTheta - 3)
{
- m_interpolateTheta = 0.;
+ parms.m_interpolateTheta = 0.;
thetaBin = s_maxBinTheta - 2;
}
- m_complementTheta = 1. - m_interpolateTheta;
- // std::cout << " m_zAtAxis " << m_zAtAxis
- // << " m_cotTheta " << m_cotTheta
+ parms.m_complementTheta = 1. - parms.m_interpolateTheta;
+ // std::cout << " zAtAxis " << zAtAxis
+ // << " cotTheta " << cotTheta
// << " zBin " << zBin
- // << " interpolateZ " << m_interpolateZ
+ // << " interpolateZ " << parms.m_interpolateZ
// << " thetaBin " << thetaBin
- // << " interpolateTheta " << m_interpolateTheta << std::endl;
+ // << " interpolateTheta " << parms.m_interpolateTheta << std::endl;
return 2*s_numberParameters*(s_maxBinTheta*zBin + thetaBin);
}
@@ -143,7 +166,8 @@ inline void
SolenoidParametrization::integrate(double& firstIntegral,
double& secondIntegral,
double zBegin,
- double zEnd) const
+ double zEnd,
+ const Parameters& parms) const
{
double zDiff = zEnd - zBegin;
double zBeg2 = zBegin*zBegin;
@@ -152,29 +176,30 @@ SolenoidParametrization::integrate(double& firstIntegral,
double zEnd3 = zEnd2*zEnd;
double zDiff4 = 0.25*(zEnd2 + zBeg2)*(zEnd2 - zBeg2);
- firstIntegral += m_fieldAtOrigin*zDiff +
- m_quadraticTerm*(zEnd3 - zBeg3)*0.333333333333 +
- m_cubicTerm*zDiff4;
+ firstIntegral += parms.m_fieldAtOrigin*zDiff +
+ parms.m_quadraticTerm*(zEnd3 - zBeg3)*0.333333333333 +
+ parms.m_cubicTerm*zDiff4;
double zDiffInv = 1./zDiff;
- secondIntegral += m_fieldAtOrigin*zDiff +
- m_quadraticTerm*(zDiffInv*zDiff4 - zBeg3)*0.666666666667 +
- m_cubicTerm*(0.1*zDiffInv*(zEnd2*zEnd3 - zBeg2*zBeg3) - 0.5*zBeg2*zBeg2);
+ secondIntegral += parms.m_fieldAtOrigin*zDiff +
+ parms.m_quadraticTerm*(zDiffInv*zDiff4 - zBeg3)*0.666666666667 +
+ parms.m_cubicTerm*(0.1*zDiffInv*(zEnd2*zEnd3 - zBeg2*zBeg3) - 0.5*zBeg2*zBeg2);
}
inline double
SolenoidParametrization::interpolate(int key1,
- int key2,
- int key3,
- int key4) const
+ int key2,
+ int key3,
+ int key4,
+ const Parameters& parms) const
{
- return ((s_parameters[key1]*m_complementZ +
- s_parameters[key2]*m_interpolateZ)*m_complementTheta +
- (s_parameters[key3]*m_complementZ +
- s_parameters[key4]*m_interpolateZ)*m_interpolateTheta);
+ return ((m_parameters[key1]*parms.m_complementZ +
+ m_parameters[key2]*parms.m_interpolateZ)*parms.m_complementTheta +
+ (m_parameters[key3]*parms.m_complementZ +
+ m_parameters[key4]*parms.m_interpolateZ)*parms.m_interpolateTheta);
}
inline void
-SolenoidParametrization::setTerms(int key1)
+SolenoidParametrization::setTerms(int key1, Parameters& parms) const
{
int key2 = key1 + s_numberParameters;
int key3 = key2 + s_numberParameters;
@@ -182,52 +207,34 @@ SolenoidParametrization::setTerms(int key1)
assert (key1 >= 0);
assert (key4 < 14688);
- assert (s_parameters[key1] != 0.);
- assert (s_parameters[key3] != 0.);
- if (m_cotTheta < 7.)
+ assert (m_parameters[key1] != 0.);
+ assert (m_parameters[key3] != 0.);
+ if (parms.m_cotTheta < 7.)
{
- assert (s_parameters[key2] != 0.);
- assert (s_parameters[key4] != 0.);
+ assert (m_parameters[key2] != 0.);
+ assert (m_parameters[key4] != 0.);
}
- m_fieldAtOrigin = interpolate(key1++,key2++,key3++,key4++);
- m_quadraticTerm = interpolate(key1++,key2++,key3++,key4++);
- m_cubicTerm = interpolate(key1,key2,key3,key4);
+ parms.m_fieldAtOrigin = interpolate(key1++,key2++,key3++,key4++,parms);
+ parms.m_quadraticTerm = interpolate(key1++,key2++,key3++,key4++,parms);
+ parms.m_cubicTerm = interpolate(key1,key2,key3,key4,parms);
}
//<<<<<< INLINE PUBLIC MEMBER FUNCTIONS >>>>>>
-// class does not work without magnetic field service
-//inline SolenoidParametrization*
-//SolenoidParametrization::instance()
-//{
-// if (s_instance == 0) s_instance = new SolenoidParametrization();
-// return s_instance;
-//}
-
-inline SolenoidParametrization*
-SolenoidParametrization::instance(MagField::IMagFieldSvc* magFieldSvc)
-{
- // this method is provided to configure the initialization.
- // It should not be called more than once without an intervening clearInstance().
- assert (s_instance == 0);
- s_instance = new SolenoidParametrization(magFieldSvc);
- return s_instance;
-}
-
inline double
SolenoidParametrization::centralField (void) const
-{ return s_centralField; }
+{ return m_centralField; }
inline double
-SolenoidParametrization::fieldComponent (double z) const
+SolenoidParametrization::fieldComponent (double z, const Parameters& parms) const
{
- double z_local = m_signTheta*z - m_zAtAxis;
+ double z_local = parms.m_signTheta*z - parms.m_zAtAxis;
double z_squared = z_local*z_local;
- double value = m_fieldAtOrigin +
- m_quadraticTerm*z_squared +
- m_cubicTerm*z_squared*z_local;
+ double value = parms.m_fieldAtOrigin +
+ parms.m_quadraticTerm*z_squared +
+ parms.m_cubicTerm*z_squared*z_local;
return value;
}
@@ -244,32 +251,29 @@ inline void
SolenoidParametrization::fieldIntegrals(double& firstIntegral,
double& secondIntegral,
double zBegin,
- double zEnd)
+ double zEnd,
+ Parameters& parms) const
{
- zBegin = m_signTheta*zBegin;
- zEnd = m_signTheta*zEnd;
+ zBegin = parms.m_signTheta*zBegin;
+ zEnd = parms.m_signTheta*zEnd;
if (zEnd < s_zInner || zBegin > s_zInner)
{
- integrate(firstIntegral,secondIntegral,zBegin-m_zAtAxis,zEnd-m_zAtAxis);
+ integrate(firstIntegral,secondIntegral,zBegin-parms.m_zAtAxis,zEnd-parms.m_zAtAxis, parms);
}
else
{
- integrate(firstIntegral,secondIntegral,zBegin-m_zAtAxis,s_zInner);
- int key = fieldKey() + s_numberParameters/2;
- setTerms(key);
- integrate(firstIntegral,secondIntegral,s_zInner,zEnd-m_zAtAxis);
+ integrate(firstIntegral,secondIntegral,zBegin-parms.m_zAtAxis,s_zInner, parms);
+ int key = fieldKey(parms) + s_numberParameters/2;
+ setTerms(key, parms);
+ integrate(firstIntegral,secondIntegral,s_zInner,zEnd-parms.m_zAtAxis,parms);
}
// std::cout << " zBegin < 0. " << zBegin
// << " zEnd " << zEnd
- // << " m_signTheta " << m_signTheta
- // << " m_zAtAxis " << m_zAtAxis << std::endl;
+ // << " m_signTheta " << parms.m_signTheta
+ // << " m_zAtAxis " << parms.m_zAtAxis << std::endl;
}
-inline bool
-SolenoidParametrization::hasParametrization (void) const
-{ return m_hasParametrization; }
-
inline double
SolenoidParametrization::maximumR (void) const
{ return s_rInner; }
@@ -278,32 +282,8 @@ inline double
SolenoidParametrization::maximumZ (void) const
{ return s_zInner; }
-inline void
-SolenoidParametrization::setParameters (double r, double z, double cotTheta)
-{
- if (cotTheta > 0.)
- {
- m_signTheta = 1.;
- m_cotTheta = cotTheta;
- m_zAtAxis = z - r*cotTheta;
- }
- else
- {
- m_signTheta = -1.;
- m_cotTheta = -cotTheta;
- m_zAtAxis = r*cotTheta - z;
- }
- int key = fieldKey();
- if (r > s_rInner || m_signTheta*z > s_zInner)
- {
- key += s_numberParameters/2;
- }
-// cout << "r " << r << " key " << key << endl;
- setTerms(key);
-}
-
inline bool
-SolenoidParametrization::validOrigin(Amg::Vector3D origin) const
+SolenoidParametrization::validOrigin(const Amg::Vector3D& origin) const
{ return (origin.perp() < s_maximumImpactAtOrigin
&& std::abs(origin.z()) < s_maximumZatOrigin); }
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h
index 94224bcc987..1ef2b7b75f8 100755
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h
@@ -2,282 +2,287 @@
Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
-///////////////////////////////////////////////////////////////////
-// SolenoidalIntersector.h, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#ifndef TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
-#define TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
-
-#include <cmath>
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/IIncidentListener.h"
-#include "GaudiKernel/ServiceHandle.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "GeoPrimitives/GeoPrimitives.h"
-#include "MagFieldInterfaces/IMagFieldSvc.h"
-#include "TrkExInterfaces/IIntersector.h"
-#include "TrkExUtils/TrackSurfaceIntersection.h"
-#include "TrkExSolenoidalIntersector/SolenoidParametrization.h"
-
-class IIncidentSvc;
-namespace Trk
-{
-
-class SolenoidalIntersector: public AthAlgTool,
- virtual public IIntersector, virtual public IIncidentListener
-{
-
-public:
- SolenoidalIntersector (const std::string& type,
- const std::string& name,
- const IInterface* parent);
- ~SolenoidalIntersector (void); // destructor
-
- StatusCode initialize();
- StatusCode finalize();
-
- /** handle for incident service */
- void handle(const Incident& inc) ;
-
- /**IIntersector interface method for general Surface type */
- const TrackSurfaceIntersection* intersectSurface(const Surface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method for specific Surface type : PerigeeSurface */
- const TrackSurfaceIntersection* approachPerigeeSurface(const PerigeeSurface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method for specific Surface type : StraightLineSurface */
- const TrackSurfaceIntersection* approachStraightLineSurface(const StraightLineSurface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method for specific Surface type : CylinderSurface */
- const TrackSurfaceIntersection* intersectCylinderSurface (const CylinderSurface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method for specific Surface type : DiscSurface */
- const TrackSurfaceIntersection* intersectDiscSurface (const DiscSurface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method for specific Surface type : PlaneSurface */
- const TrackSurfaceIntersection* intersectPlaneSurface(const PlaneSurface& surface,
- const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
- const double qOverP);
-
- /**IIntersector interface method to check validity of parametrization within extrapolation range */
- bool isValid (Amg::Vector3D startPosition,
- Amg::Vector3D endPosition) const;
-
- /** tabulate parametrization details */
- void validationAction() const;
-
-private:
- double circularArcLength(double, double, double, double, double,
- double, double&, double&);
- double linearArcLength(double);
- bool extrapolateToR(double endRadius);
- bool extrapolateToZ(double endZ);
- const TrackSurfaceIntersection* intersection(const Surface& surface);
- void setParameters(const TrackSurfaceIntersection* intersection,
- double qOverP);
-
- // services and tools:
- ServiceHandle<IIncidentSvc> m_incidentSvc; //!< IncidentSvc to catch begin of event
- ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc;
- ToolHandle<IIntersector> m_rungeKuttaIntersector;
-
- double m_cotTheta;
- float m_currentMax;
- float m_currentMin;
- double m_deltaPhiTolerance;
- Amg::Vector3D m_direction;
- unsigned m_intersectionNumber;
- double m_oneOverSinTheta;
- double m_pathLength;
- Amg::Vector3D m_position;
- double m_qOverP;
- double m_qOverPt;
- double m_radius;
- double m_sinTheta;
- SolenoidParametrization* m_solenoidParametrization;
- double m_surfaceTolerance;
- double m_validRadius;
- double m_validZ;
-
- // counters
- unsigned long long m_countExtrapolations;
- unsigned long long m_countRKSwitches;
-
-};
-
-
-// arc length to intersect of 2 circles: circular track and circle at (0,0) with radius endRadius
-inline double
-SolenoidalIntersector::circularArcLength(double endRadius,
- double radiusOfCurvature,
- double xCentre,
- double yCentre,
- double cosPhi,
- double sinPhi,
- double& cosPhiIntersect,
- double& sinPhiIntersect)
-{
- int trapped = 0;
- double radiusSquared = xCentre*xCentre + yCentre*yCentre;
- double term = 0.5*(radiusSquared +
- radiusOfCurvature*radiusOfCurvature -
- endRadius*endRadius)/(radiusSquared*radiusOfCurvature);
- if (std::abs(xCentre) < std::abs(yCentre))
- {
- double dx2 = yCentre*yCentre * (1./(radiusSquared*term*term) - 1.);
- if (dx2 < 0.)
- {
- trapped = 1;
- }
- else
- {
- if (yCentre*term > 0.)
- {
- sinPhiIntersect = term*(-xCentre+std::sqrt(dx2));
- }
- else
- {
- sinPhiIntersect = term*(-xCentre-std::sqrt(dx2));
- }
- cosPhiIntersect = (sinPhiIntersect*xCentre + radiusSquared*term)/yCentre;
- }
- }
- else
- {
- double dy2 = xCentre*xCentre * (1./(radiusSquared*term*term) - 1.);
- if (dy2 < 0.)
- {
- trapped = 1;
- }
- else
- {
- if (xCentre*term > 0.)
- {
- cosPhiIntersect = term*(yCentre+std::sqrt(dy2));
- }
- else
- {
- cosPhiIntersect = term*(yCentre-std::sqrt(dy2));
- }
- sinPhiIntersect = (cosPhiIntersect*yCentre - radiusSquared*term)/xCentre;
- }
- }
- if (trapped == 0)
- {
- double deltaPhi;
- double sinDeltaPhi = sinPhiIntersect*cosPhi - cosPhiIntersect*sinPhi;
- if (std::abs(sinDeltaPhi) > 0.1)
- {
- deltaPhi = asin(sinDeltaPhi);
- }
- else
- {
- deltaPhi = sinDeltaPhi*(1. + 0.166667*sinDeltaPhi*sinDeltaPhi);
- }
- return (radiusOfCurvature*deltaPhi);
- }
- else
- {
- cosPhiIntersect = cosPhi;
- sinPhiIntersect = sinPhi;
- return 0.;
- }
-}
-
-// arc length to intersect of a line to a circle of radius endRadius centred at (0,0)
-// +ve (-ve) endRadius selects the solution on the same (opposite) side of (0,0)
-inline double
-SolenoidalIntersector::linearArcLength(double endRadius)
-{
- double arcLength = (-m_direction.x()*m_position.x() - m_direction.y()*m_position.y()) *
- m_oneOverSinTheta;
- double radiusSquared = endRadius*endRadius - m_radius*m_radius + arcLength*arcLength;
- if (radiusSquared > 0.)
- {
- if (endRadius > 0.)
- {
- arcLength += std::sqrt(radiusSquared);
- }
- else
- {
- arcLength -= std::sqrt(radiusSquared);
- }
- }
- return arcLength;
-}
-
-inline const TrackSurfaceIntersection*
-SolenoidalIntersector::intersection(const Surface& surface)
-{
- Intersection SLIntersect = surface.straightLineIntersection(m_position, m_direction, false, false);
- if (! SLIntersect.valid) return 0;
-
- const TrackSurfaceIntersection* intersection = new TrackSurfaceIntersection(SLIntersect.position,
- m_direction,
- m_pathLength);
- // // validate
- // if (! intersection)
- // {
- // ATH_MSG_WARNING(" this should never fail");
- // return 0;
- // }
-
- // const Amg::Vector2D* localPos = surface.positionOnSurface(intersection->position(), false);
- // if (! localPos)
- // {
- // ATH_MSG_INFO(" localPos fails surface type " << surface.type()
- // << " at R,Z " << m_position.perp() << ", " << m_position.z()
- // << " surface R " << surface.globalReferencePoint().perp());
- // return 0;
- // }
-
-
- // ATH_MSG_INFO(" serial diff " << intersection->serialNumber() - m_intersectionNumber
- // << " at R,Z: " << m_radius << ", " << m_position.z());
-
- m_intersectionNumber = intersection->serialNumber();
- return intersection;
-}
-
-inline void
-SolenoidalIntersector::setParameters(const TrackSurfaceIntersection* trackTrackSurfaceIntersection, double qOverP)
-{
- if (trackTrackSurfaceIntersection->serialNumber() != m_intersectionNumber || qOverP != m_qOverP)
- {
- // ATH_MSG_INFO(" initialize parameters. Diff: " << trackTrackSurfaceIntersection->serialNumber()-m_intersectionNumber
- // << " at R,Z: " << trackTrackSurfaceIntersection->position().perp() << ", " << trackTrackSurfaceIntersection->position().z());
- ++m_countExtrapolations;
- m_position.x() = trackTrackSurfaceIntersection->position().x();
- m_position.y() = trackTrackSurfaceIntersection->position().y();
- m_position.z() = trackTrackSurfaceIntersection->position().z();
- m_radius = m_position.perp();
- m_direction.x() = trackTrackSurfaceIntersection->direction().x();
- m_direction.y() = trackTrackSurfaceIntersection->direction().y();
- m_direction.z() = trackTrackSurfaceIntersection->direction().z();
- m_sinTheta = m_direction.perp();
- m_oneOverSinTheta = 1./m_sinTheta;
- m_cotTheta = m_direction.z() * m_oneOverSinTheta;
- m_pathLength = trackTrackSurfaceIntersection->pathlength();
- m_qOverP = qOverP;
- m_qOverPt = qOverP * m_oneOverSinTheta;
- m_solenoidParametrization->setParameters(m_radius,m_position.z(),m_cotTheta);
- }
-}
-
-} // end of namespace
-
-
-#endif // TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
-
-
+///////////////////////////////////////////////////////////////////
+// SolenoidalIntersector.h, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#ifndef TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
+#define TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GaudiKernel/ContextSpecificPtr.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+#include "TrkExInterfaces/IIntersector.h"
+#include "TrkExUtils/TrackSurfaceIntersection.h"
+#include "TrkExSolenoidalIntersector/SolenoidParametrization.h"
+#include "CxxUtils/checker_macros.h"
+#include <mutex>
+#include <cmath>
+
+class IIncidentSvc;
+namespace Trk
+{
+
+class SolenoidalIntersector: public extends<AthAlgTool, IIntersector>
+{
+
+public:
+ /**
+ * @brief Constants of motion and other cached values.
+ *
+ * There is some data we want to persist across calls to intersectSurface()
+ * for a single trajectory. We do this by attaching this structure
+ * as a cache block to TrackSurfaceIntersection. Some of these are
+ * actually constants of motion for a helical trajectory. We also
+ * include here the solenoidal parametrization. This gets updated
+ * as we step through the field; including it here avoids a reevaluation
+ * at the start of the next call. We also include the position at the
+ * end of the stepping phase, m_lastPosition. This is because before
+ * after we finish stepping but before we return, we improve the
+ * estimate of the intersection by calculating the straight line
+ * intersection with the surface. But if we then calcuate a further
+ * intersection with the same trajectory, we want to start at the
+ * point at the end of the previous stepping, not the position
+ * we returned.
+ */
+ struct Constants
+ : public TrackSurfaceIntersection::IIntersectionCache
+ {
+ Constants (const SolenoidParametrization& solpar,
+ const TrackSurfaceIntersection& trackTrackSurfaceIntersection,
+ const double qOverP);
+ virtual std::unique_ptr<IIntersectionCache> clone() const override
+ { return std::make_unique<Constants> (*this); }
+
+ double m_sinTheta;
+ double m_oneOverSinTheta;
+ double m_cotTheta;
+ double m_qOverPt;
+ const SolenoidParametrization& m_solPar;
+ Amg::Vector3D m_lastPosition;
+ SolenoidParametrization::Parameters m_solParams;
+ };
+ SolenoidalIntersector (const std::string& type,
+ const std::string& name,
+ const IInterface* parent);
+ ~SolenoidalIntersector (void); // destructor
+
+ StatusCode initialize();
+ StatusCode finalize();
+
+ /**IIntersector interface method for general Surface type */
+ const TrackSurfaceIntersection* intersectSurface(const Surface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : PerigeeSurface */
+ const TrackSurfaceIntersection* approachPerigeeSurface(const PerigeeSurface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : StraightLineSurface */
+ const TrackSurfaceIntersection* approachStraightLineSurface(const StraightLineSurface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : CylinderSurface */
+ const TrackSurfaceIntersection* intersectCylinderSurface (const CylinderSurface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : DiscSurface */
+ const TrackSurfaceIntersection* intersectDiscSurface (const DiscSurface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : PlaneSurface */
+ const TrackSurfaceIntersection* intersectPlaneSurface(const PlaneSurface& surface,
+ const TrackSurfaceIntersection* trackTrackSurfaceIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method to check validity of parametrization within extrapolation range */
+ bool isValid (Amg::Vector3D startPosition,
+ Amg::Vector3D endPosition) const;
+
+ /** tabulate parametrization details */
+ void validationAction() const;
+
+private:
+ const SolenoidParametrization* getSolenoidParametrization() const;
+
+ double circularArcLength(double, double, double, double, double,
+ double, double&, double&) const;
+ double linearArcLength(const TrackSurfaceIntersection& isect,
+ const Constants& com,
+ const double radius2,
+ const double endRadius) const;
+ bool extrapolateToR(TrackSurfaceIntersection& isect,
+ double& radius2,
+ Constants& com,
+ const double endRadius) const;
+ bool extrapolateToZ(TrackSurfaceIntersection& isect,
+ Constants& com,
+ const double endZ) const;
+ const TrackSurfaceIntersection* intersection(std::unique_ptr<TrackSurfaceIntersection> isect,
+ Constants& com,
+ const Surface& surface) const;
+
+ std::unique_ptr<TrackSurfaceIntersection>
+ newIntersection (const TrackSurfaceIntersection& oldIsect,
+ const SolenoidParametrization& solpar,
+ const double qOverP,
+ Constants*& com) const;
+
+ ServiceHandle<MagField::IMagFieldSvc> m_magFieldSvc;
+ ToolHandle<IIntersector> m_rungeKuttaIntersector;
+
+ double m_deltaPhiTolerance;
+ double m_surfaceTolerance;
+
+ // counters
+ mutable std::atomic<unsigned long long> m_countExtrapolations;
+ mutable std::atomic<unsigned long long> m_countRKSwitches;
+
+ mutable std::mutex m_mutex;
+ mutable Gaudi::Hive::ContextSpecificPtr<const SolenoidParametrization> m_lastSolenoidParametrization ATLAS_THREAD_SAFE;
+ // List of active solenoid parametrizations. Second element of the pair
+ // is a use count.
+ typedef std::list<std::pair<SolenoidParametrization, int> > Parmlist_t;
+ mutable Parmlist_t m_solenoidParametrizations ATLAS_THREAD_SAFE;
+};
+
+
+// arc length to intersect of 2 circles: circular track and circle at (0,0) with radius endRadius
+inline double
+SolenoidalIntersector::circularArcLength(double endRadius,
+ double radiusOfCurvature,
+ double xCentre,
+ double yCentre,
+ double cosPhi,
+ double sinPhi,
+ double& cosPhiIntersect,
+ double& sinPhiIntersect) const
+{
+ int trapped = 0;
+ double radiusSquared = xCentre*xCentre + yCentre*yCentre;
+ double term = 0.5*(radiusSquared +
+ radiusOfCurvature*radiusOfCurvature -
+ endRadius*endRadius)/(radiusSquared*radiusOfCurvature);
+ if (std::abs(xCentre) < std::abs(yCentre))
+ {
+ double dx2 = yCentre*yCentre * (1./(radiusSquared*term*term) - 1.);
+ if (dx2 < 0.)
+ {
+ trapped = 1;
+ }
+ else
+ {
+ if (yCentre*term > 0.)
+ {
+ sinPhiIntersect = term*(-xCentre+std::sqrt(dx2));
+ }
+ else
+ {
+ sinPhiIntersect = term*(-xCentre-std::sqrt(dx2));
+ }
+ cosPhiIntersect = (sinPhiIntersect*xCentre + radiusSquared*term)/yCentre;
+ }
+ }
+ else
+ {
+ double dy2 = xCentre*xCentre * (1./(radiusSquared*term*term) - 1.);
+ if (dy2 < 0.)
+ {
+ trapped = 1;
+ }
+ else
+ {
+ if (xCentre*term > 0.)
+ {
+ cosPhiIntersect = term*(yCentre+std::sqrt(dy2));
+ }
+ else
+ {
+ cosPhiIntersect = term*(yCentre-std::sqrt(dy2));
+ }
+ sinPhiIntersect = (cosPhiIntersect*yCentre - radiusSquared*term)/xCentre;
+ }
+ }
+ if (trapped == 0)
+ {
+ double deltaPhi;
+ double sinDeltaPhi = sinPhiIntersect*cosPhi - cosPhiIntersect*sinPhi;
+ if (std::abs(sinDeltaPhi) > 0.1)
+ {
+ deltaPhi = asin(sinDeltaPhi);
+ }
+ else
+ {
+ deltaPhi = sinDeltaPhi*(1. + 0.166667*sinDeltaPhi*sinDeltaPhi);
+ }
+ return (radiusOfCurvature*deltaPhi);
+ }
+ else
+ {
+ cosPhiIntersect = cosPhi;
+ sinPhiIntersect = sinPhi;
+ return 0.;
+ }
+}
+
+// arc length to intersect of a line to a circle of radius endRadius centred at (0,0)
+// +ve (-ve) endRadius selects the solution on the same (opposite) side of (0,0)
+inline double
+SolenoidalIntersector::linearArcLength(const TrackSurfaceIntersection& isect,
+ const Constants& com,
+ const double radius2,
+ const double endRadius) const
+{
+ const Amg::Vector3D& pos = isect.position();
+ const Amg::Vector3D& dir = isect.direction();
+
+ double arcLength = (-dir.x()*pos.x() - dir.y()*pos.y()) *
+ com.m_oneOverSinTheta;
+ double radiusSquared = endRadius*endRadius - radius2 + arcLength*arcLength;
+ if (radiusSquared > 0.)
+ {
+ if (endRadius > 0.)
+ {
+ arcLength += std::sqrt(radiusSquared);
+ }
+ else
+ {
+ arcLength -= std::sqrt(radiusSquared);
+ }
+ }
+ return arcLength;
+}
+
+inline const TrackSurfaceIntersection*
+SolenoidalIntersector::intersection(std::unique_ptr<TrackSurfaceIntersection> isect,
+ Constants& com,
+ const Surface& surface) const
+{
+ // Improve the estimate of the intersection by calculating
+ // the straight-line intersection.
+ Intersection SLIntersect = surface.straightLineIntersection(isect->position(), isect->direction(), false, false);
+ if (SLIntersect.valid)
+ {
+ // But first save our current position, so that we can
+ // start from here on the next call.
+ com.m_lastPosition = isect->position();
+ isect->position() = SLIntersect.position;
+ return isect.release();
+ }
+
+ return nullptr;
+}
+
+} // end of namespace
+
+
+#endif // TRKEXSOLENOIDALINTERSECTOR_SOLENOIDALINTERSECTOR_H
+
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/share/SolenoidalIntersector_test.ref b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/share/SolenoidalIntersector_test.ref
index ed9b945fd76..e22da500173 100644
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/share/SolenoidalIntersector_test.ref
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/share/SolenoidalIntersector_test.ref
@@ -10,7 +10,7 @@ JobOptionsSvc INFO Job options successfully read in from ../share/Solenoi
ApplicationMgr SUCCESS
====================================================================================================================================
Welcome to ApplicationMgr (GaudiCoreSvc v27r1p99)
- running on karma on Tue Jan 29 12:11:05 2019
+ running on karma on Sun Feb 3 21:41:27 2019
====================================================================================================================================
ApplicationMgr INFO Application Manager Configured successfully
AtlasFieldSvc INFO initialize() ...
@@ -34,9 +34,8 @@ AtlasFieldSvc INFO Initialized the field map from /home/sss/nobackup/atla
AtlasFieldSvc INFO Currents imported and map initialized
AtlasFieldSvc INFO BeginRun incident handled
AtlasFieldSvc INFO incidents handled successfully
- SolenoidParametrization: centralField 2T - please be patient while the solenoid is parametrised !!
-ToolSvc.Trk::So... INFO SolenoidParametrization current: 7730 valid radius,Z : 570, 2150
test_plane
+ SolenoidParametrization: centralField 2T - please be patient while the solenoid is parametrised !!
test_line
test_cylinder
test_disc
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx
index 21236565d5b..3dc8e769667 100755
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2017, 2019 CERN for the benefit of the ATLAS collaboration
*/
/***************************************************************************
@@ -46,8 +46,6 @@ const double SolenoidParametrization::s_binInvSizeTheta = 1./0.1;
const double SolenoidParametrization::s_binInvSizeZ = 1./20.*Gaudi::Units::mm;
const double SolenoidParametrization::s_binZeroTheta = 0.;
const double SolenoidParametrization::s_binZeroZ = -160.*Gaudi::Units::mm;
-double SolenoidParametrization::s_centralField = 0.;
-SolenoidParametrization* SolenoidParametrization::s_instance = 0;
const double SolenoidParametrization::s_lightSpeed = -1.*299792458*Gaudi::Units::m/Gaudi::Units::s;
const int SolenoidParametrization::s_maxBinTheta = 72;
const int SolenoidParametrization::s_maxBinZ = 17;
@@ -58,38 +56,74 @@ const double SolenoidParametrization::s_rInner = 570.*Gaudi::Units::mm;
const double SolenoidParametrization::s_rOuter = 1050.*Gaudi::Units::mm;
const double SolenoidParametrization::s_zInner = 2150.0*Gaudi::Units::mm; // just after wheel #7
const double SolenoidParametrization::s_zOuter = 2800.0*Gaudi::Units::mm; // just after wheel #9
-double SolenoidParametrization::s_parameters[] = {14688*0.};
+
+
+SolenoidParametrization::BinParameters::BinParameters (const double r,
+ const double z,
+ const double cotTheta)
+{
+ if (cotTheta > 0) {
+ m_signTheta = 1;
+ m_cotTheta = cotTheta;
+ m_zAtAxis = z - r*cotTheta;
+ }
+ else {
+ m_signTheta = -1;
+ m_cotTheta = -cotTheta;
+ m_zAtAxis = r*cotTheta - z;
+ }
+}
+
+
+SolenoidParametrization::Parameters::Parameters (const SolenoidParametrization& spar,
+ const double r,
+ const double z,
+ const double cotTheta)
+ : BinParameters (r, z, cotTheta)
+{
+ int key = fieldKey(*this);
+ if (r > s_rInner || m_signTheta*z > s_zInner)
+ {
+ key += s_numberParameters/2;
+ }
+ spar.setTerms (key, *this);
+}
+
//<<<<<< CLASS STRUCTURE INITIALIZATION >>>>>>
// note: both constructors are private
// SolenoidParametrization::SolenoidParametrization(void)
-// : m_hasParametrization (false)
// {
// // field component not defined without a m_magFieldSvc
// // get central field value in units required for fast tracking
-// s_centralField = fieldComponent(0.,0.,0.);
+// m_centralField = fieldComponent(0.,0.,0.);
//
// // note field = B*c
// RestoreIOSFlags restore(std::cout);
// std::cout << " SolenoidParametrization: centralField "
-// << std::setw(6) << std::setprecision(3) << s_centralField/s_lightSpeed /Gaudi::Units::tesla
+// << std::setw(6) << std::setprecision(3) << m_centralField/s_lightSpeed /Gaudi::Units::tesla
// << "T" << std::endl;
//}
SolenoidParametrization::SolenoidParametrization(MagField::IMagFieldSvc* magFieldSvc)
- : m_hasParametrization (false),
- m_magFieldSvc (magFieldSvc)
+ : m_magFieldSvc (magFieldSvc),
+ m_parameters ()
{
- // get central field value in units required for fast tracking (i.e. field = B*c)
+ // allow 0.1% fluctuation in current.
+ double current = magFieldSvc->solenoidCurrent();
+ m_currentMax = 1.001*current;
+ m_currentMin = 0.999*current;
+
+ // get central field value in units required for fast tracking (i.e. field = B*c)
if (!m_magFieldSvc)
throw std::logic_error("fieldComponent not defined without magnetic field service.");
- s_centralField = fieldComponent(0.,0.,0.);
+ m_centralField = fieldComponent(0.,0.,0.);
RestoreIOSFlags restore(std::cout);
std::cout << " SolenoidParametrization: centralField "
- << std::setw(6) << std::setprecision(3) << s_centralField/s_lightSpeed /Gaudi::Units::tesla
+ << std::setw(6) << std::setprecision(3) << m_centralField/s_lightSpeed /Gaudi::Units::tesla
<< "T";
-
+
// now parametrise field - if requested
//if (magFieldSvc)
{
@@ -104,31 +138,28 @@ SolenoidParametrization::SolenoidParametrization(MagField::IMagFieldSvc* magFiel
void
SolenoidParametrization::parametrizeSolenoid(void)
{
- m_hasParametrization = true;
-
// set parametrisation granularity (up to cotTheta = 7.)
// get value of cubic term for approx: Bz = Bcentral*(1 - term * z^3)
// 'fit' to average over cotTheta lines
- m_signTheta = 1.;
double smallOffset = 0.0000000000001; // avoid FPE
- m_zAtAxis = s_binZeroZ; // + smallOffset ?
+ double zAtAxis = s_binZeroZ; // + smallOffset ?
for (int binZ = 0; binZ < s_maxBinZ; ++binZ)
{
- m_cotTheta = smallOffset;
+ double cotTheta = smallOffset;
for (int binTheta = 0; binTheta < s_maxBinTheta - 1; ++binTheta)
{
double r = 0.;
- double z = m_zAtAxis;
+ double z = zAtAxis;
int n = 200;
double dr;
- if (m_cotTheta < s_zOuter/s_rOuter)
+ if (cotTheta < s_zOuter/s_rOuter)
{
dr = s_rOuter/double(n);
}
else
{
- dr = s_zOuter/(m_cotTheta*double(n));
+ dr = s_zOuter/(cotTheta*double(n));
}
Amg::VectorX difference(n);
@@ -136,9 +167,9 @@ SolenoidParametrization::parametrizeSolenoid(void)
for (int k = 0; k < n; ++k)
{
r += dr;
- z += dr*m_cotTheta;
+ z += dr*cotTheta;
double w = (n - k)*(n - k);
- double zLocal = z - m_zAtAxis;
+ double zLocal = z - zAtAxis;
if (r > s_rInner || z > s_zInner)
{
// derivative(k,0) = 0.;
@@ -157,31 +188,32 @@ SolenoidParametrization::parametrizeSolenoid(void)
// derivative(k,4) = 0.;
// derivative(k,5) = 0.;
}
- difference(k) = w*(fieldComponent(r,z,m_cotTheta) - s_centralField);
+ difference(k) = w*(fieldComponent(r,z,cotTheta) - m_centralField);
}
// solve for parametrization coefficients
Amg::VectorX solution = derivative.colPivHouseholderQr().solve(difference);
+ BinParameters parms (zAtAxis, cotTheta);
- int key = fieldKey();
- assert (s_parameters[key] == 0.);
- s_parameters[key++] = s_centralField + solution(0);
- s_parameters[key++] = solution(1);
- s_parameters[key++] = solution(2);
- s_parameters[key++] = s_centralField + solution(3);
- s_parameters[key++] = solution(4);
- s_parameters[key++] = solution(5);
+ int key = fieldKey(parms);
+ assert (m_parameters[key] == 0.);
+ m_parameters[key++] = m_centralField + solution(0);
+ m_parameters[key++] = solution(1);
+ m_parameters[key++] = solution(2);
+ m_parameters[key++] = m_centralField + solution(3);
+ m_parameters[key++] = solution(4);
+ m_parameters[key++] = solution(5);
// duplicate last z-bin for contiguous neighbour lookup
if (binZ == s_maxBinZ - 1)
{
- assert (s_parameters[key] == 0.);
- s_parameters[key++] = s_centralField + solution(0);
- s_parameters[key++] = solution(1);
- s_parameters[key++] = solution(2);
- s_parameters[key++] = s_centralField + solution(3);
- s_parameters[key++] = solution(4);
- s_parameters[key++] = solution(5);
+ assert (m_parameters[key] == 0.);
+ m_parameters[key++] = m_centralField + solution(0);
+ m_parameters[key++] = solution(1);
+ m_parameters[key++] = solution(2);
+ m_parameters[key++] = m_centralField + solution(3);
+ m_parameters[key++] = solution(4);
+ m_parameters[key++] = solution(5);
key -= s_numberParameters;
}
@@ -189,13 +221,13 @@ SolenoidParametrization::parametrizeSolenoid(void)
if (binZ > 0)
{
key -= 2*s_numberParameters*s_maxBinTheta;
- assert (s_parameters[key] == 0.);
- s_parameters[key++] = s_centralField + solution(0);
- s_parameters[key++] = solution(1);
- s_parameters[key++] = solution(2);
- s_parameters[key++] = s_centralField + solution(3);
- s_parameters[key++] = solution(4);
- s_parameters[key++] = solution(5);
+ assert (m_parameters[key] == 0.);
+ m_parameters[key++] = m_centralField + solution(0);
+ m_parameters[key++] = solution(1);
+ m_parameters[key++] = solution(2);
+ m_parameters[key++] = m_centralField + solution(3);
+ m_parameters[key++] = solution(4);
+ m_parameters[key++] = solution(5);
}
// some debug print
@@ -204,66 +236,54 @@ SolenoidParametrization::parametrizeSolenoid(void)
// || key >= s_numberParameters*s_maxBinTheta*(s_maxBinZ - 2))
// {
// double z_max;
-// if (m_cotTheta < s_zInner/s_rInner)
+// if (cotTheta < s_zInner/s_rInner)
// {
-// z_max = s_rInner*m_cotTheta + m_zAtAxis;
+// z_max = s_rInner*cotTheta + zAtAxis;
// }
// else
// {
// z_max = s_zInner;
// }
// cout << std::setiosflags(std::ios::fixed) << key
-// << " cotTheta " << std::setw(7) << std::setprecision(2) << m_cotTheta
+// << " cotTheta " << std::setw(7) << std::setprecision(2) << cotTheta
// << " inner terms: z0 "<< std::setw(6) << std::setprecision(3)
-// << s_parameters[key]/s_centralField
+// << m_parameters[key]/m_centralField
// << " z^2 "<< std::setw(6) << std::setprecision(3)
-// << s_parameters[key+1]*z_max*z_max/s_centralField
+// << m_parameters[key+1]*z_max*z_max/m_centralField
// << " z^3 " << std::setw(6) << std::setprecision(3)
-// << s_parameters[key+2]*z_max*z_max*z_max/s_centralField
+// << m_parameters[key+2]*z_max*z_max*z_max/m_centralField
// << " outer terms: z0 "<< std::setw(6) << std::setprecision(3)
-// << s_parameters[key+3]/s_centralField
+// << m_parameters[key+3]/m_centralField
// << " z^2 "<< std::setw(6) << std::setprecision(3)
-// << s_parameters[key+4]*z_max*z_max/s_centralField
+// << m_parameters[key+4]*z_max*z_max/m_centralField
// << " z^3 " << std::setw(6) << std::setprecision(3)
-// << s_parameters[key+5]*z_max*z_max*z_max/s_centralField
+// << m_parameters[key+5]*z_max*z_max*z_max/m_centralField
// << std::resetiosflags(std::ios::fixed) << endl;
// }
- m_cotTheta += 1./s_binInvSizeTheta;
+ cotTheta += 1./s_binInvSizeTheta;
}
- m_zAtAxis += 1./s_binInvSizeZ;
+ zAtAxis += 1./s_binInvSizeZ;
}
// duplicate end theta-bins for contiguous neighbour lookup
- m_zAtAxis = s_binZeroZ; // + smallOffset ??
+ zAtAxis = s_binZeroZ; // + smallOffset ??
for (int binZ = 0; binZ < s_maxBinZ; ++binZ)
{
- m_cotTheta = double(s_maxBinTheta)/s_binInvSizeTheta;
- int key = fieldKey();
+ double cotTheta = double(s_maxBinTheta)/s_binInvSizeTheta;
+ BinParameters parms (zAtAxis, cotTheta);
+ int key = fieldKey(parms);
for (int k = 0; k < 2*s_numberParameters; ++k)
{
- assert (s_parameters[key+2*s_numberParameters] == 0.);
- s_parameters[key+2*s_numberParameters] = s_parameters[key];
+ assert (m_parameters[key+2*s_numberParameters] == 0.);
+ m_parameters[key+2*s_numberParameters] = m_parameters[key];
++key;
}
- m_zAtAxis += 1./s_binInvSizeZ;
+ zAtAxis += 1./s_binInvSizeZ;
}
}
//<<<<<< PUBLIC MEMBER FUNCTION DEFINITIONS >>>>>>
-SolenoidParametrization::~SolenoidParametrization()
-{}
-
-void
-SolenoidParametrization::clearInstance(void)
-{
- delete s_instance;
- s_instance = 0;
- std::fill (s_parameters,
- s_parameters + sizeof(s_parameters)/sizeof(s_parameters[0]),
- 0);
-}
-
void
SolenoidParametrization::printFieldIntegrals (void) const
{
@@ -392,16 +412,15 @@ SolenoidParametrization::printFieldIntegrals (void) const
}
void
-SolenoidParametrization::printParametersForEtaLine (double eta, double z_origin)
+SolenoidParametrization::printParametersForEtaLine (double eta, double z_origin) const
{
- m_signTheta = 1.;
- m_zAtAxis = z_origin;
- m_cotTheta = 1./std::tan(2.*std::atan(1./std::exp(eta)));
- int key = fieldKey();
+ double cotTheta = 1./std::tan(2.*std::atan(1./std::exp(eta)));
+ BinParameters parms (z_origin, cotTheta);
+ int key = fieldKey(parms);
double z_max;
- if (m_cotTheta < s_zInner/s_rInner)
+ if (cotTheta < s_zInner/s_rInner)
{
- z_max = s_rInner*m_cotTheta;
+ z_max = s_rInner*cotTheta;
}
else
{
@@ -411,25 +430,23 @@ SolenoidParametrization::printParametersForEtaLine (double eta, double z_origin)
<< "SolenoidParametrization: line with eta " << std::setw(6) << std::setprecision(2) << eta
<< " from (r,z) 0.0," << std::setw(6) << std::setprecision(1) << z_origin
<< " inner terms: z0 "<< std::setw(6) << std::setprecision(2)
- << s_parameters[key]/s_centralField
+ << m_parameters[key]/m_centralField
<< " z^2 "<< std::setw(6) << std::setprecision(3)
- << s_parameters[key+1]*z_max*z_max/s_centralField
+ << m_parameters[key+1]*z_max*z_max/m_centralField
<< " z^3 " << std::setw(6) << std::setprecision(3)
- << s_parameters[key+2]*z_max*z_max*z_max/s_centralField
+ << m_parameters[key+2]*z_max*z_max*z_max/m_centralField
<< " outer terms: z0 "<< std::setw(6) << std::setprecision(3)
- << s_parameters[key+3]/s_centralField
+ << m_parameters[key+3]/m_centralField
<< " z^2 "<< std::setw(6) << std::setprecision(3)
- << s_parameters[key+4]*z_max*z_max/s_centralField
+ << m_parameters[key+4]*z_max*z_max/m_centralField
<< " z^3 " << std::setw(6) << std::setprecision(3)
- << s_parameters[key+5]*z_max*z_max*z_max/s_centralField
+ << m_parameters[key+5]*z_max*z_max*z_max/m_centralField
<< std::resetiosflags(std::ios::fixed) << std::endl;
}
void
-SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin)
+SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin) const
{
- m_signTheta = 1.;
- m_zAtAxis = zOrigin;
double cotTheta = 1./std::tan(2.*std::atan(1./std::exp(std::abs(eta))));
double z = zOrigin;
double r = 0.;
@@ -455,8 +472,8 @@ SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin)
for (int k = 0; k < n; ++k)
{
double b = fieldComponent(r,z,cotTheta);
- setParameters(r,z,cotTheta);
- double diff = (fieldComponent(z) - b)/s_lightSpeed;
+ Parameters parms (*this, r, z, cotTheta);
+ double diff = (fieldComponent(z, parms) - b)/s_lightSpeed;
if (r > s_rInner || z > s_zInner)
{
@@ -472,7 +489,7 @@ SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin)
if (std::abs(diff) > worstDiff)
{
worstDiff = std::abs(diff);
- worstBCalc = fieldComponent(z);
+ worstBCalc = fieldComponent(z, parms);
worstBTrue = b;
worstR = r;
worstZ = z;
@@ -496,6 +513,14 @@ SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin)
<< std::resetiosflags(std::ios::fixed) << std::endl;
}
+
+bool SolenoidParametrization::currentMatches (double current) const
+{
+ // allow 0.1% fluctuation
+ return current >= m_currentMin && current < m_currentMax;
+}
+
+
} // end of namespace
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx
index a39a899db30..dea5a22d910 100755
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx
@@ -1,466 +1,545 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2017, 2019 CERN for the benefit of the ATLAS collaboration
*/
-///////////////////////////////////////////////////////////////////
-// SolenoidalIntersector.cxx, (c) ATLAS Detector software
-///////////////////////////////////////////////////////////////////
-
-#include "GaudiKernel/Incident.h"
-#include "GaudiKernel/SystemOfUnits.h"
-#include "TrkExSolenoidalIntersector/SolenoidalIntersector.h"
-#include "TrkParameters/TrackParameters.h"
-#include "TrkSurfaces/CylinderSurface.h"
-#include "TrkSurfaces/DiscSurface.h"
-#include "TrkSurfaces/PerigeeSurface.h"
-#include "TrkSurfaces/PlaneSurface.h"
-#include "TrkSurfaces/StraightLineSurface.h"
-#include "TrkSurfaces/Surface.h"
-
-namespace Trk
-{
-SolenoidalIntersector::SolenoidalIntersector (const std::string& type,
- const std::string& name,
- const IInterface* parent)
- : AthAlgTool (type, name, parent),
- m_incidentSvc ("IncidentSvc", name),
- m_magFieldSvc ("MagField::AtlasFieldSvc/AtlasFieldSvc", name),
- m_rungeKuttaIntersector ("Trk::RungeKuttaIntersector/RungeKuttaIntersector"),
- m_cotTheta (0.),
- m_currentMax (0.),
- m_currentMin (0.),
- m_deltaPhiTolerance (0.01), // upper limit for small angle approx
- m_intersectionNumber (0),
- m_oneOverSinTheta (1.),
- m_pathLength (0.),
- m_qOverP (0.),
- m_qOverPt (0.),
- m_radius (0.),
- m_sinTheta (1.),
- m_solenoidParametrization (0),
- m_surfaceTolerance (2.0*Gaudi::Units::micrometer),
- m_validRadius (0.),
- m_validZ (0.),
- m_countExtrapolations (0),
- m_countRKSwitches (0)
-{
- declareInterface<Trk::IIntersector>(this);
- declareProperty("MagFieldSvc", m_magFieldSvc );
- declareProperty("RungeKuttaIntersector", m_rungeKuttaIntersector);
- declareProperty("SurfaceTolerance", m_surfaceTolerance);
-}
-
-SolenoidalIntersector::~SolenoidalIntersector (void)
-{}
-
-StatusCode
-SolenoidalIntersector::initialize()
-{
- // print name and package version
- ATH_MSG_INFO( "SolenoidalIntersector::initialize() - package version " << PACKAGE_VERSION );
-
- if (m_incidentSvc.retrieve().isFailure())
- {
- ATH_MSG_FATAL( "Failed to retrieve service " << m_incidentSvc );
- return StatusCode::FAILURE;
- }
- // register to the incident service:
- // BeginEvent needed to ensure the solenoid field parametrization is still valid
- m_incidentSvc->addListener( this, IncidentType::BeginEvent);
-
- if (! m_magFieldSvc.empty())
- {
- if (m_magFieldSvc.retrieve().isFailure())
- {
- ATH_MSG_FATAL( "Failed to retrieve service " << m_magFieldSvc );
- return StatusCode::FAILURE;
- }
- else
- {
- ATH_MSG_INFO( "Retrieved service " << m_magFieldSvc );
- }
- }
-
- if (m_rungeKuttaIntersector.retrieve().isFailure())
- {
- ATH_MSG_FATAL( "Failed to retrieve tool " << m_rungeKuttaIntersector );
- return StatusCode::FAILURE;
- }
- else
- {
- ATH_MSG_INFO( "Retrieved tool " << m_rungeKuttaIntersector );
- }
-
- return StatusCode::SUCCESS;
-}
-
-StatusCode
-SolenoidalIntersector::finalize()
-{
- ATH_MSG_INFO( "finalized after " << m_countExtrapolations << " extrapolations with "
- << m_countRKSwitches << " switches to RK integration");
-
- return StatusCode::SUCCESS;
-}
-
-/** handle for incident service */
-void
-SolenoidalIntersector::handle(const Incident& inc)
-{
- // (re)parametrize the solenoidal field whenever necessary at BeginEvent
- if (inc.type() == IncidentType::BeginEvent )
- {
- // intersector invalid when field off
- if (! m_magFieldSvc->solenoidOn())
- {
- m_currentMax = 0.;
- m_validRadius = 0.;
- m_validZ = 0.;
- return;
- }
-
- if (m_magFieldSvc->solenoidCurrent() > m_currentMax
- || m_magFieldSvc->solenoidCurrent() < m_currentMin)
- {
- // allow 0.1% fluctuation
- m_currentMax = 1.001*m_magFieldSvc->solenoidCurrent();
- m_currentMin = 0.999*m_magFieldSvc->solenoidCurrent();
-
- // instantiate the field parametrisation
- SolenoidParametrization::clearInstance();
- m_solenoidParametrization = SolenoidParametrization::instance(&*m_magFieldSvc);
-
- // set limit of validity (cylinder bounds)
- if (! m_solenoidParametrization)
- {
- m_validRadius = 0.;
- m_validZ = 0.;
- ATH_MSG_WARNING(" mag field parametrization fails" );
- return;
- }
- m_validRadius = m_solenoidParametrization->maximumR();
- m_validZ = m_solenoidParametrization->maximumZ();
- ATH_MSG_INFO( " SolenoidParametrization current: " << m_magFieldSvc->solenoidCurrent()
- << " valid radius,Z : " << m_validRadius << ", " << m_validZ);
-
- if (msgLvl(MSG::DEBUG)) validationAction();
- }
- }
-}
-
-/**IIntersector interface method for general Surface type */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::intersectSurface(const Surface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{
- const PlaneSurface* plane = dynamic_cast<const PlaneSurface*>(&surface);
- if (plane) return intersectPlaneSurface(*plane,trackIntersection,qOverP);
-
- const StraightLineSurface* straightLine = dynamic_cast<const StraightLineSurface*>(&surface);
- if (straightLine) return m_rungeKuttaIntersector->approachStraightLineSurface
- (*straightLine, trackIntersection, qOverP);
-
- const CylinderSurface* cylinder = dynamic_cast<const CylinderSurface*>(&surface);
- if (cylinder) return intersectCylinderSurface(*cylinder,trackIntersection,qOverP);
-
- const DiscSurface* disc = dynamic_cast<const DiscSurface*>(&surface);
- if (disc) return intersectDiscSurface(*disc,trackIntersection,qOverP);
-
- const PerigeeSurface* perigee = dynamic_cast<const PerigeeSurface*>(&surface);
- if (perigee) return m_rungeKuttaIntersector->approachPerigeeSurface
- (*perigee, trackIntersection, qOverP);
-
- ATH_MSG_WARNING( " unrecognized Surface" );
- return 0;
-}
-
-/**IIntersector interface method for specific Surface type : PerigeeSurface */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::approachPerigeeSurface(const PerigeeSurface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{ return m_rungeKuttaIntersector->approachPerigeeSurface(surface, trackIntersection, qOverP); }
-
-/**IIntersector interface method for specific Surface type : StraightLineSurface */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::approachStraightLineSurface(const StraightLineSurface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{ return m_rungeKuttaIntersector->approachStraightLineSurface(surface, trackIntersection, qOverP); }
-
-/**IIntersector interface method for specific Surface type : CylinderSurface */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::intersectCylinderSurface(const CylinderSurface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{
- double endRadius = surface.globalReferencePoint().perp();
- if (endRadius > m_validRadius)
- return m_rungeKuttaIntersector->intersectCylinderSurface(surface, trackIntersection, qOverP);
-
- // set member data unless extrapolation from existing parameters
- setParameters(trackIntersection, qOverP);
-
- if (std::abs(endRadius -trackIntersection->position().perp()) > m_surfaceTolerance
- && ! extrapolateToR(endRadius)) return 0;
- return intersection(surface);
-}
-
-/**IIntersector interface method for specific Surface type : DiscSurface */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::intersectDiscSurface (const DiscSurface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{
- double endZ = surface.center().z();
- if (std::abs(endZ) > m_validZ)
- return m_rungeKuttaIntersector->intersectDiscSurface(surface, trackIntersection, qOverP);
-
- // set member data unless extrapolation from existing parameters
- setParameters(trackIntersection, qOverP);
- if (std::abs(endZ -trackIntersection->position().z()) > m_surfaceTolerance
- && ! extrapolateToZ(endZ))
- {
- return 0;
- }
-
- return intersection(surface);
-}
-
-/**IIntersector interface method for specific Surface type : PlaneSurface */
-const Trk::TrackSurfaceIntersection*
-SolenoidalIntersector::intersectPlaneSurface(const PlaneSurface& surface,
- const TrackSurfaceIntersection* trackIntersection,
- const double qOverP)
-{
- if (std::abs(surface.center().z()) > m_validZ
- || surface.center().perp() > m_validRadius)
- return m_rungeKuttaIntersector->intersectPlaneSurface(surface, trackIntersection, qOverP);
-
- // set member data unless extrapolation from existing parameters
- setParameters(trackIntersection, qOverP);
-
- // step until sufficiently near to plane surface
- // this gives wrong answer!
- // DistanceSolution solution = surface.straightLineDistanceEstimate(m_position,m_direction);
- // if (! solution.signedDistance()) return 0;
- // double distance = solution.currentDistance(true);
-
- int numberSteps = 0;
- double dot = surface.normal().dot(m_direction);
- double offset = surface.normal().dot(surface.center() - m_position);
-
- while (std::abs(offset) > m_surfaceTolerance*std::abs(dot))
- {
- // take care if grazing incidence - quit if looper
- if (std::abs(dot) < 0.0001) return 0;
- double distance = offset/dot;
-
- // extrapolate
- if (m_sinTheta < 0.9)
- {
- if (! extrapolateToZ(m_position.z()+distance*m_direction.z())) return 0;
- }
- else
- {
- if (! extrapolateToR((m_position+distance*m_direction).perp())) return 0;
- }
-
- // check we are getting closer to the plane, switch to RK in case of difficulty
- dot = surface.normal().dot(m_direction);
- offset = surface.normal().dot(surface.center() - m_position);
- if (std::abs(offset) > m_surfaceTolerance*std::abs(dot)
- && (++numberSteps > 5 || std::abs(offset) > 0.5*std::abs(distance*dot)))
- {
- ++m_countRKSwitches;
- ATH_MSG_DEBUG(" switch to RK after " << numberSteps << " steps at offset "
- << offset << " dot " << surface.normal().dot(m_direction));
-
- return m_rungeKuttaIntersector->intersectPlaneSurface(surface, trackIntersection, qOverP);
- }
- };
-
- return intersection(surface);
-}
-
-/**IIntersector interface method to check validity of parametrization within extrapolation range */
-bool
-SolenoidalIntersector::isValid (Amg::Vector3D startPosition, Amg::Vector3D endPosition) const
-{
- // check cylinder bounds for valid parametrization
- if (std::abs(endPosition.z()) < m_validZ
- && endPosition.perp() < m_validRadius
- && m_solenoidParametrization->validOrigin(startPosition))
- {
- // ATH_MSG_INFO(" choose solenoidal");
-
- return true;
- }
-
- // ATH_MSG_INFO(" choose rungeKutta");
-
- return false;
-}
-
-/** tabulate parametrization details */
-void
-SolenoidalIntersector::validationAction() const
-{
- // validate parametrization
- if (m_solenoidParametrization)
- {
- for (int ieta = 0; ieta != 27; ++ieta)
- {
- double eta = 0.05 + 0.1*static_cast<double>(ieta);
- m_solenoidParametrization->printParametersForEtaLine(+eta,0.);
- m_solenoidParametrization->printParametersForEtaLine(-eta,0.);
- }
- for (int ieta = 0; ieta != 27; ++ieta)
- {
- double eta = 0.05 + 0.1*static_cast<double>(ieta);
- m_solenoidParametrization->printResidualForEtaLine(+eta,0.);
- m_solenoidParametrization->printResidualForEtaLine(-eta,0.);
- }
- m_solenoidParametrization->printFieldIntegrals();
- }
-}
-
-// private methods
-
-
-bool
-SolenoidalIntersector::extrapolateToR(double endR)
-{
- // ATH_MSG_INFO(" extrapolateToR endR " << endR << " from " << m_position.z()
- // << " r " << m_radius << " cotTheta " << m_cotTheta);
-
- double fieldComponent = m_solenoidParametrization->fieldComponent(m_position.z());
- double curvature = fieldComponent*m_qOverPt;
- double arcLength = linearArcLength(endR);
- if (std::abs(arcLength*curvature) > m_deltaPhiTolerance)
- {
- double radiusOfCurvature = 1./curvature;
- double xCentre = m_position.x() -
- radiusOfCurvature*m_direction.y()*m_oneOverSinTheta;
- double yCentre = m_position.y() +
- radiusOfCurvature*m_direction.x()*m_oneOverSinTheta;
- double cosPhi;
- double sinPhi;
- arcLength = circularArcLength(endR,
- radiusOfCurvature,
- xCentre,
- yCentre,
- m_direction.x() * m_oneOverSinTheta,
- m_direction.y() * m_oneOverSinTheta,
- cosPhi,
- sinPhi);
- if (std::abs(arcLength*m_cotTheta) < m_surfaceTolerance)
- {
- m_position.x() = xCentre + radiusOfCurvature*sinPhi;
- m_position.y() = yCentre - radiusOfCurvature*cosPhi;
- m_position.z() += arcLength*m_cotTheta;
- m_radius = endR;
- m_direction.x() = cosPhi*m_sinTheta;
- m_direction.y() = sinPhi*m_sinTheta;
- m_pathLength += arcLength*m_oneOverSinTheta;
- arcLength = 0.;
- }
- }
-
- double deltaZ = arcLength*m_cotTheta;
- if (std::abs(deltaZ) < m_surfaceTolerance) // avoid FPE with constant curvature parabolic approx
- {
- double cosPhi = m_direction.x()*m_oneOverSinTheta;
- double sinPhi = m_direction.y()*m_oneOverSinTheta;
- if (std::abs(arcLength) > m_surfaceTolerance)
- {
- double sinDPhi = 0.5*arcLength*curvature;
- double cosDPhi = 1. - 0.5*sinDPhi*sinDPhi *
- (1.0+0.25*sinDPhi*sinDPhi);
- double temp = cosPhi;
- cosPhi = temp*cosDPhi - sinPhi*sinDPhi;
- sinPhi = temp*sinDPhi + sinPhi*cosDPhi;
- m_direction.x() = (cosPhi*cosDPhi - sinPhi*sinDPhi)*m_sinTheta;
- m_direction.y() = (sinPhi*cosDPhi + cosPhi*sinDPhi)*m_sinTheta;
- }
-
- m_position.x() += arcLength*cosPhi;
- m_position.y() += arcLength*sinPhi;
- m_position.z() += arcLength*m_cotTheta;
- m_radius = endR;
- m_pathLength += arcLength*m_oneOverSinTheta;
- }
- else
- {
- extrapolateToZ(m_position.z() + deltaZ);
- if (std::abs(endR - m_radius) > m_surfaceTolerance)
- {
- deltaZ = linearArcLength(endR) * m_cotTheta;
- extrapolateToZ(m_position.z() + deltaZ);
- }
- }
-
- return true;
-}
-
-bool
-SolenoidalIntersector::extrapolateToZ(double endZ)
-{
- double firstIntegral = 0.;
- double secondIntegral = 0.;
- m_solenoidParametrization->fieldIntegrals(firstIntegral,
- secondIntegral,
- m_position.z(),
- endZ);
- // ATH_MSG_INFO(" extrapolateToZ firstIntegral, secondIntegral " << 1.E6*firstIntegral
- // << ", " << 1.E6*secondIntegral);
- double DFiMax = 0.1;
- double cosPhi = m_direction.x()*m_oneOverSinTheta;
- double sinPhi = m_direction.y()*m_oneOverSinTheta;
- double cosBeta;
- double sinBeta;
- double deltaPhi2 = secondIntegral*m_qOverPt/std::abs(m_cotTheta);
- if (std::abs(deltaPhi2) < DFiMax)
- {
- double deltaPhi2Squared= deltaPhi2*deltaPhi2;
- sinBeta = 1. - 0.166667*deltaPhi2Squared;
- cosBeta = -0.5*deltaPhi2*(1.-0.083333*deltaPhi2Squared);
- }
- else if (2.*std::abs(deltaPhi2) < M_PI)
- {
- sinBeta = std::sin(deltaPhi2) / deltaPhi2;
- cosBeta = (std::cos(deltaPhi2) - 1.) / deltaPhi2;
- }
- else
- {
- return false;
- }
-
- double radialDistance = (endZ - m_position.z()) / m_cotTheta;
- m_position.x() += radialDistance*(sinPhi*cosBeta + cosPhi*sinBeta);
- m_position.y() += radialDistance*(sinPhi*sinBeta - cosPhi*cosBeta);
- m_position.z() = endZ;
- m_radius = m_position.perp();
- m_pathLength += radialDistance*m_oneOverSinTheta;
-
- double cosAlpha;
- double sinAlpha;
- double deltaPhi1 = firstIntegral*m_qOverPt / std::abs(m_cotTheta);
- if (std::abs(deltaPhi1) < DFiMax)
- {
- double deltaPhi1Squared= deltaPhi1*deltaPhi1;
- sinAlpha = deltaPhi1*(1. - 0.166667*deltaPhi1Squared);
- cosAlpha = 1. - 0.5*deltaPhi1Squared*(1.-0.083333*deltaPhi1Squared);
- }
- else
- {
- sinAlpha = std::sin(deltaPhi1);
- cosAlpha = std::cos(deltaPhi1);
- }
- m_direction.x() = (cosPhi*cosAlpha - sinPhi*sinAlpha)*m_sinTheta;
- m_direction.y() = (sinPhi*cosAlpha + cosPhi*sinAlpha)*m_sinTheta;;
- return true;
-}
-
-} // end of namespace
+///////////////////////////////////////////////////////////////////
+// SolenoidalIntersector.cxx, (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+
+#include "GaudiKernel/SystemOfUnits.h"
+#include "TrkExSolenoidalIntersector/SolenoidalIntersector.h"
+#include "TrkParameters/TrackParameters.h"
+#include "TrkSurfaces/CylinderSurface.h"
+#include "TrkSurfaces/DiscSurface.h"
+#include "TrkSurfaces/PerigeeSurface.h"
+#include "TrkSurfaces/PlaneSurface.h"
+#include "TrkSurfaces/StraightLineSurface.h"
+#include "TrkSurfaces/Surface.h"
+
+namespace Trk
+{
+
+SolenoidalIntersector::Constants::Constants (const SolenoidParametrization& solpar,
+ const TrackSurfaceIntersection& trackTrackSurfaceIntersection,
+ const double qOverP)
+ : m_sinTheta (trackTrackSurfaceIntersection.direction().perp()),
+ m_oneOverSinTheta (1./m_sinTheta),
+ m_cotTheta (trackTrackSurfaceIntersection.direction().z() * m_oneOverSinTheta),
+ m_qOverPt (qOverP * m_oneOverSinTheta),
+ m_solPar (solpar),
+ m_lastPosition (trackTrackSurfaceIntersection.position()),
+ m_solParams (solpar,
+ trackTrackSurfaceIntersection.position().perp(),
+ trackTrackSurfaceIntersection.position().z(),
+ m_cotTheta)
+{
+}
+
+
+SolenoidalIntersector::SolenoidalIntersector (const std::string& type,
+ const std::string& name,
+ const IInterface* parent)
+ : base_class (type, name, parent),
+ m_magFieldSvc ("MagField::AtlasFieldSvc/AtlasFieldSvc", name),
+ m_rungeKuttaIntersector ("Trk::RungeKuttaIntersector/RungeKuttaIntersector"),
+ m_deltaPhiTolerance (0.01), // upper limit for small angle approx
+ m_surfaceTolerance (2.0*Gaudi::Units::micrometer),
+ m_countExtrapolations (0),
+ m_countRKSwitches (0)
+{
+ declareInterface<Trk::IIntersector>(this);
+ declareProperty("MagFieldSvc", m_magFieldSvc );
+ declareProperty("RungeKuttaIntersector", m_rungeKuttaIntersector);
+ declareProperty("SurfaceTolerance", m_surfaceTolerance);
+}
+
+SolenoidalIntersector::~SolenoidalIntersector (void)
+{}
+
+StatusCode
+SolenoidalIntersector::initialize()
+{
+ // print name and package version
+ ATH_MSG_INFO( "SolenoidalIntersector::initialize() - package version " << PACKAGE_VERSION );
+
+ if (! m_magFieldSvc.empty())
+ {
+ if (m_magFieldSvc.retrieve().isFailure())
+ {
+ ATH_MSG_FATAL( "Failed to retrieve service " << m_magFieldSvc );
+ return StatusCode::FAILURE;
+ }
+ else
+ {
+ ATH_MSG_INFO( "Retrieved service " << m_magFieldSvc );
+ }
+ }
+
+ if (m_rungeKuttaIntersector.retrieve().isFailure())
+ {
+ ATH_MSG_FATAL( "Failed to retrieve tool " << m_rungeKuttaIntersector );
+ return StatusCode::FAILURE;
+ }
+ else
+ {
+ ATH_MSG_INFO( "Retrieved tool " << m_rungeKuttaIntersector );
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode
+SolenoidalIntersector::finalize()
+{
+ ATH_MSG_INFO( "finalized after " << m_countExtrapolations << " extrapolations with "
+ << m_countRKSwitches << " switches to RK integration");
+
+ return StatusCode::SUCCESS;
+}
+
+
+const SolenoidParametrization*
+SolenoidalIntersector::getSolenoidParametrization() const
+{
+ double current = m_magFieldSvc->solenoidCurrent();
+
+ // Check to see if the last one we used is ok.
+ const SolenoidParametrization* lastpar = m_lastSolenoidParametrization.get();
+ if (lastpar && lastpar->currentMatches (current)) {
+ return lastpar;
+ }
+
+ std::lock_guard<std::mutex> lock (m_mutex);
+
+ // Search for a new one, and release the reference count on the old one.
+ const SolenoidParametrization* thispar = nullptr;
+ Parmlist_t::iterator todel = m_solenoidParametrizations.end();
+ for (Parmlist_t::iterator it = m_solenoidParametrizations.begin();
+ it != m_solenoidParametrizations.end() && (lastpar || !thispar);
+ ++it)
+ {
+ if (&it->first == lastpar) {
+ lastpar = nullptr;
+ if (--it->second <= 0) {
+ todel = it;
+ }
+ }
+ else if (!thispar && it->first.currentMatches (current)) {
+ thispar = &it->first;
+ ++it->second;
+ }
+ }
+
+ if (todel != m_solenoidParametrizations.end()) {
+ m_solenoidParametrizations.erase (todel);
+ }
+
+ if (!thispar) {
+ // Didn't find one; make a new one.
+ m_solenoidParametrizations.emplace_back (&*m_magFieldSvc, 1);
+ thispar = &m_solenoidParametrizations.back().first;
+ }
+
+ // Remember which one we used last.
+ m_lastSolenoidParametrization.set (thispar);
+
+ return thispar;
+}
+
+
+/**IIntersector interface method for general Surface type */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::intersectSurface(const Surface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{
+ const PlaneSurface* plane = dynamic_cast<const PlaneSurface*>(&surface);
+ if (plane) return intersectPlaneSurface(*plane,trackIntersection,qOverP);
+
+ const StraightLineSurface* straightLine = dynamic_cast<const StraightLineSurface*>(&surface);
+ if (straightLine) return m_rungeKuttaIntersector->approachStraightLineSurface
+ (*straightLine, trackIntersection, qOverP);
+
+ const CylinderSurface* cylinder = dynamic_cast<const CylinderSurface*>(&surface);
+ if (cylinder) return intersectCylinderSurface(*cylinder,trackIntersection,qOverP);
+
+ const DiscSurface* disc = dynamic_cast<const DiscSurface*>(&surface);
+ if (disc) return intersectDiscSurface(*disc,trackIntersection,qOverP);
+
+ const PerigeeSurface* perigee = dynamic_cast<const PerigeeSurface*>(&surface);
+ if (perigee) return m_rungeKuttaIntersector->approachPerigeeSurface
+ (*perigee, trackIntersection, qOverP);
+
+ ATH_MSG_WARNING( " unrecognized Surface" );
+ return 0;
+}
+
+/**IIntersector interface method for specific Surface type : PerigeeSurface */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::approachPerigeeSurface(const PerigeeSurface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{ return m_rungeKuttaIntersector->approachPerigeeSurface(surface, trackIntersection, qOverP); }
+
+/**IIntersector interface method for specific Surface type : StraightLineSurface */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::approachStraightLineSurface(const StraightLineSurface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{ return m_rungeKuttaIntersector->approachStraightLineSurface(surface, trackIntersection, qOverP); }
+
+/**IIntersector interface method for specific Surface type : CylinderSurface */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::intersectCylinderSurface(const CylinderSurface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{
+ const SolenoidParametrization* solenoidParametrization =
+ getSolenoidParametrization();
+
+ double endRadius = surface.globalReferencePoint().perp();
+ if (!solenoidParametrization || endRadius > solenoidParametrization->maximumR())
+ return m_rungeKuttaIntersector->intersectCylinderSurface(surface, trackIntersection, qOverP);
+
+ Constants* com = nullptr;
+ std::unique_ptr<TrackSurfaceIntersection> isect =
+ newIntersection (*trackIntersection,
+ *solenoidParametrization,
+ qOverP,
+ com);
+ ++m_countExtrapolations;
+
+ double radius2 = isect->position().perp2();
+ if (std::abs(endRadius - sqrt(radius2)) > m_surfaceTolerance
+ && ! extrapolateToR(*isect, radius2, *com, endRadius)) return 0;
+ return intersection(std::move(isect), *com, surface);
+}
+
+/**IIntersector interface method for specific Surface type : DiscSurface */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::intersectDiscSurface (const DiscSurface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{
+ const SolenoidParametrization* solenoidParametrization =
+ getSolenoidParametrization();
+
+ double endZ = surface.center().z();
+ if (!solenoidParametrization || std::abs(endZ) > solenoidParametrization->maximumZ())
+ return m_rungeKuttaIntersector->intersectDiscSurface(surface, trackIntersection, qOverP);
+
+ Constants* com = nullptr;
+ std::unique_ptr<TrackSurfaceIntersection> isect =
+ newIntersection (*trackIntersection,
+ *solenoidParametrization,
+ qOverP,
+ com);
+
+ ++m_countExtrapolations;
+
+ if (std::abs(endZ -trackIntersection->position().z()) > m_surfaceTolerance
+ && ! extrapolateToZ(*isect, *com, endZ))
+ {
+ return 0;
+ }
+
+ return intersection(std::move(isect), *com, surface);
+}
+
+/**IIntersector interface method for specific Surface type : PlaneSurface */
+const Trk::TrackSurfaceIntersection*
+SolenoidalIntersector::intersectPlaneSurface(const PlaneSurface& surface,
+ const TrackSurfaceIntersection* trackIntersection,
+ const double qOverP)
+{
+ const SolenoidParametrization* solenoidParametrization =
+ getSolenoidParametrization();
+
+ if (!solenoidParametrization ||
+ std::abs(surface.center().z()) > solenoidParametrization->maximumZ()
+ || surface.center().perp() > solenoidParametrization->maximumR())
+ return m_rungeKuttaIntersector->intersectPlaneSurface(surface, trackIntersection, qOverP);
+
+ Constants* com = nullptr;
+ std::unique_ptr<TrackSurfaceIntersection> isect =
+ newIntersection (*trackIntersection,
+ *solenoidParametrization,
+ qOverP,
+ com);
+ Amg::Vector3D& pos = isect->position();
+ Amg::Vector3D& dir = isect->direction();
+ ++m_countExtrapolations;
+ double radius2 = pos.perp2();
+
+ // step until sufficiently near to plane surface
+ // this gives wrong answer!
+ // DistanceSolution solution = surface.straightLineDistanceEstimate(pos,dir);
+ // if (! solution.signedDistance()) return 0;
+ // double distance = solution.currentDistance(true);
+
+ int numberSteps = 0;
+ double dot = surface.normal().dot(dir);
+ double offset = surface.normal().dot(surface.center() - pos);
+
+ while (std::abs(offset) > m_surfaceTolerance*std::abs(dot))
+ {
+ // take care if grazing incidence - quit if looper
+ if (std::abs(dot) < 0.0001) return 0;
+ double distance = offset/dot;
+
+ // extrapolate
+ if (com->m_sinTheta < 0.9)
+ {
+ if (! extrapolateToZ(*isect, *com, pos.z()+distance*dir.z())) return 0;
+ radius2 = pos.perp2();
+ }
+ else
+ {
+ if (! extrapolateToR(*isect, radius2, *com, (pos+distance*dir).perp())) return 0;
+ }
+
+ // check we are getting closer to the plane, switch to RK in case of difficulty
+ dot = surface.normal().dot(dir);
+ offset = surface.normal().dot(surface.center() - pos);
+ if (std::abs(offset) > m_surfaceTolerance*std::abs(dot)
+ && (++numberSteps > 5 || std::abs(offset) > 0.5*std::abs(distance*dot)))
+ {
+ ++m_countRKSwitches;
+ ATH_MSG_DEBUG(" switch to RK after " << numberSteps << " steps at offset "
+ << offset << " dot " << surface.normal().dot(dir));
+
+ return m_rungeKuttaIntersector->intersectPlaneSurface(surface, trackIntersection, qOverP);
+ }
+ };
+
+ return intersection(std::move(isect), *com, surface);
+}
+
+/**IIntersector interface method to check validity of parametrization within extrapolation range */
+bool
+SolenoidalIntersector::isValid (Amg::Vector3D startPosition, Amg::Vector3D endPosition) const
+{
+ const SolenoidParametrization* solenoidParametrization =
+ getSolenoidParametrization();
+
+ // check cylinder bounds for valid parametrization
+ if (solenoidParametrization &&
+ std::abs(endPosition.z()) < solenoidParametrization->maximumZ()
+ && endPosition.perp() < solenoidParametrization->maximumR()
+ && getSolenoidParametrization()->validOrigin(startPosition))
+ {
+ // ATH_MSG_INFO(" choose solenoidal");
+
+ return true;
+ }
+
+ // ATH_MSG_INFO(" choose rungeKutta");
+
+ return false;
+}
+
+/** tabulate parametrization details */
+void
+SolenoidalIntersector::validationAction() const
+{
+ const SolenoidParametrization* solenoidParametrization =
+ getSolenoidParametrization();
+
+ // validate parametrization
+ if (solenoidParametrization)
+ {
+ for (int ieta = 0; ieta != 27; ++ieta)
+ {
+ double eta = 0.05 + 0.1*static_cast<double>(ieta);
+ solenoidParametrization->printParametersForEtaLine(+eta,0.);
+ solenoidParametrization->printParametersForEtaLine(-eta,0.);
+ }
+ for (int ieta = 0; ieta != 27; ++ieta)
+ {
+ double eta = 0.05 + 0.1*static_cast<double>(ieta);
+ solenoidParametrization->printResidualForEtaLine(+eta,0.);
+ solenoidParametrization->printResidualForEtaLine(-eta,0.);
+ }
+ solenoidParametrization->printFieldIntegrals();
+ }
+}
+
+// private methods
+
+
+bool
+SolenoidalIntersector::extrapolateToR(TrackSurfaceIntersection& isect,
+ double& radius2,
+ Constants& com,
+ const double endR) const
+{
+ Amg::Vector3D& pos = isect.position();
+ Amg::Vector3D& dir = isect.direction();
+
+ // ATH_MSG_INFO(" extrapolateToR endR " << endR << " from " << pos.z()
+ // << " r " << sqrt(radius2) << " cotTheta " << com.m_cotTheta);
+
+ double fieldComponent = com.m_solPar.fieldComponent(pos.z(), com.m_solParams);
+ double curvature = fieldComponent*com.m_qOverPt;
+ double arcLength = linearArcLength(isect, com, radius2, endR);
+ if (std::abs(arcLength*curvature) > m_deltaPhiTolerance)
+ {
+ double radiusOfCurvature = 1./curvature;
+ double xCentre = pos.x() -
+ radiusOfCurvature*dir.y()*com.m_oneOverSinTheta;
+ double yCentre = pos.y() +
+ radiusOfCurvature*dir.x()*com.m_oneOverSinTheta;
+ double cosPhi;
+ double sinPhi;
+ arcLength = circularArcLength(endR,
+ radiusOfCurvature,
+ xCentre,
+ yCentre,
+ dir.x() * com.m_oneOverSinTheta,
+ dir.y() * com.m_oneOverSinTheta,
+ cosPhi,
+ sinPhi);
+ if (std::abs(arcLength*com.m_cotTheta) < m_surfaceTolerance)
+ {
+ pos.x() = xCentre + radiusOfCurvature*sinPhi;
+ pos.y() = yCentre - radiusOfCurvature*cosPhi;
+ pos.z() += arcLength*com.m_cotTheta;
+ radius2 = endR*endR;
+ dir.x() = cosPhi*com.m_sinTheta;
+ dir.y() = sinPhi*com.m_sinTheta;
+ isect.pathlength() += arcLength*com.m_oneOverSinTheta;
+ arcLength = 0.;
+ }
+ }
+
+ double deltaZ = arcLength*com.m_cotTheta;
+ if (std::abs(deltaZ) < m_surfaceTolerance) // avoid FPE with constant curvature parabolic approx
+ {
+ double cosPhi = dir.x()*com.m_oneOverSinTheta;
+ double sinPhi = dir.y()*com.m_oneOverSinTheta;
+ if (std::abs(arcLength) > m_surfaceTolerance)
+ {
+ double sinDPhi = 0.5*arcLength*curvature;
+ double cosDPhi = 1. - 0.5*sinDPhi*sinDPhi *
+ (1.0+0.25*sinDPhi*sinDPhi);
+ double temp = cosPhi;
+ cosPhi = temp*cosDPhi - sinPhi*sinDPhi;
+ sinPhi = temp*sinDPhi + sinPhi*cosDPhi;
+ dir.x() = (cosPhi*cosDPhi - sinPhi*sinDPhi)*com.m_sinTheta;
+ dir.y() = (sinPhi*cosDPhi + cosPhi*sinDPhi)*com.m_sinTheta;
+ }
+
+ pos.x() += arcLength*cosPhi;
+ pos.y() += arcLength*sinPhi;
+ pos.z() += arcLength*com.m_cotTheta;
+ radius2 = endR*endR;
+ isect.pathlength() += arcLength*com.m_oneOverSinTheta;
+ }
+ else
+ {
+ extrapolateToZ(isect, com, pos.z() + deltaZ);
+ radius2 = pos.perp2();
+ if (std::abs(endR - sqrt(radius2)) > m_surfaceTolerance)
+ {
+ deltaZ = linearArcLength(isect, com, radius2, endR) * com.m_cotTheta;
+ extrapolateToZ(isect, com, pos.z() + deltaZ);
+ radius2 = pos.perp2();
+ }
+ }
+
+ return true;
+}
+
+bool
+SolenoidalIntersector::extrapolateToZ(TrackSurfaceIntersection& isect,
+ Constants& com,
+ const double endZ) const
+{
+ Amg::Vector3D& pos = isect.position();
+ Amg::Vector3D& dir = isect.direction();
+
+ double firstIntegral = 0.;
+ double secondIntegral = 0.;
+ com.m_solPar.fieldIntegrals(firstIntegral,
+ secondIntegral,
+ pos.z(),
+ endZ,
+ com.m_solParams);
+ // ATH_MSG_INFO(" extrapolateToZ firstIntegral, secondIntegral " << 1.E6*firstIntegral
+ // << ", " << 1.E6*secondIntegral);
+ double DFiMax = 0.1;
+ double cosPhi = dir.x()*com.m_oneOverSinTheta;
+ double sinPhi = dir.y()*com.m_oneOverSinTheta;
+ double cosBeta;
+ double sinBeta;
+ double deltaPhi2 = secondIntegral*com.m_qOverPt/std::abs(com.m_cotTheta);
+ if (std::abs(deltaPhi2) < DFiMax)
+ {
+ double deltaPhi2Squared= deltaPhi2*deltaPhi2;
+ sinBeta = 1. - 0.166667*deltaPhi2Squared;
+ cosBeta = -0.5*deltaPhi2*(1.-0.083333*deltaPhi2Squared);
+ }
+ else if (2.*std::abs(deltaPhi2) < M_PI)
+ {
+ sinBeta = std::sin(deltaPhi2) / deltaPhi2;
+ cosBeta = (std::cos(deltaPhi2) - 1.) / deltaPhi2;
+ }
+ else
+ {
+ return false;
+ }
+
+ double radialDistance = (endZ - pos.z()) / com.m_cotTheta;
+ pos.x() += radialDistance*(sinPhi*cosBeta + cosPhi*sinBeta);
+ pos.y() += radialDistance*(sinPhi*sinBeta - cosPhi*cosBeta);
+ pos.z() = endZ;
+ isect.pathlength() += radialDistance*com.m_oneOverSinTheta;
+
+ double cosAlpha;
+ double sinAlpha;
+ double deltaPhi1 = firstIntegral*com.m_qOverPt / std::abs(com.m_cotTheta);
+ if (std::abs(deltaPhi1) < DFiMax)
+ {
+ double deltaPhi1Squared= deltaPhi1*deltaPhi1;
+ sinAlpha = deltaPhi1*(1. - 0.166667*deltaPhi1Squared);
+ cosAlpha = 1. - 0.5*deltaPhi1Squared*(1.-0.083333*deltaPhi1Squared);
+ }
+ else
+ {
+ sinAlpha = std::sin(deltaPhi1);
+ cosAlpha = std::cos(deltaPhi1);
+ }
+ dir.x() = (cosPhi*cosAlpha - sinPhi*sinAlpha)*com.m_sinTheta;
+ dir.y() = (sinPhi*cosAlpha + cosPhi*sinAlpha)*com.m_sinTheta;
+ return true;
+}
+
+
+std::unique_ptr<TrackSurfaceIntersection>
+SolenoidalIntersector::newIntersection (const TrackSurfaceIntersection& isect,
+ const SolenoidParametrization& solpar,
+ const double qOverP,
+ Constants*& com) const
+{
+ const TrackSurfaceIntersection::IIntersectionCache* oldCache = isect.cache();
+ std::unique_ptr<Constants> cache;
+ const Amg::Vector3D* lastPosition = nullptr;
+ if (oldCache && typeid(*oldCache) == typeid(Constants)) {
+ cache = std::make_unique<Constants> (*static_cast<const Constants*> (oldCache));
+ lastPosition = &cache->m_lastPosition;
+ }
+ else {
+ cache = std::make_unique<Constants> (solpar, isect, qOverP);
+ }
+
+ com = cache.get();
+ auto newIsect = std::make_unique<TrackSurfaceIntersection> (isect,
+ std::move (cache));
+ if (lastPosition) {
+ newIsect->position() = *lastPosition;
+ }
+ return newIsect;
+}
+
+
+} // end of namespace
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidParametrization_test.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidParametrization_test.cxx
index 305f0173d80..24348ab7dbe 100644
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidParametrization_test.cxx
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidParametrization_test.cxx
@@ -32,9 +32,7 @@ void test1 (Trk::SolenoidParametrization& sol)
{
std::cout << "test1\n";
- assert (sol.hasParametrization());
- sol.setParameters (/*r=*/30*cm, /*z=*/123*cm, /*cotTheta=*/1);
- assert (sol.hasParametrization());
+ Trk::SolenoidParametrization::Parameters parms (sol, /*r=*/30*cm, /*z=*/123*cm, /*cotTheta=*/1);
assert( sol.validOrigin (Amg::Vector3D ({10, 10, 10})) );
assert( !sol.validOrigin (Amg::Vector3D ({50, 10, 10})) );
@@ -44,14 +42,14 @@ void test1 (Trk::SolenoidParametrization& sol)
assert( Athena_test::isEqual (sol.maximumZ(), 2150) );
assert( Athena_test::isEqual (sol.centralField(), -0.598746) );
- assert( Athena_test::isEqual (sol.fieldComponent (123*cm), -0.592178) );
+ assert( Athena_test::isEqual (sol.fieldComponent (123*cm, parms), -0.592178) );
assert( Athena_test::isEqual (sol.fieldComponent (30*cm, 123*cm, 1),
-0.551903) );
double firstIntegral = 0;
double secondIntegral = 0;
sol.fieldIntegrals (firstIntegral, secondIntegral,
- 100*cm, 150*cm);
+ 100*cm, 150*cm, parms);
assert( Athena_test::isEqual (firstIntegral, -295.526) );
assert( Athena_test::isEqual (secondIntegral, -296.768) );
}
@@ -66,10 +64,8 @@ int main()
Incident inc_br ("test", IncidentType::BeginRun);
dynamic_cast<IIncidentListener*>(&*field)->handle (inc_br);
- Trk::SolenoidParametrization* sol =
- Trk::SolenoidParametrization::instance (&*field);
-
- test1 (*sol);
+ Trk::SolenoidParametrization sol (&*field);
+ test1 (sol);
return 0;
}
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidalIntersector_test.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidalIntersector_test.cxx
index 8e972966b1d..b505fe9d458 100644
--- a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidalIntersector_test.cxx
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/test/SolenoidalIntersector_test.cxx
@@ -290,8 +290,6 @@ int main()
ServiceHandle<MagField::IMagFieldSvc> field ("MagField::AtlasFieldSvc/AtlasFieldSvc", "test");
Incident inc_br ("test", IncidentType::BeginRun);
dynamic_cast<IIncidentListener*>(&*field)->handle (inc_br);
- Incident inc_be ("test", IncidentType::BeginEvent);
- dynamic_cast<IIncidentListener*>(&*tool)->handle (inc_be);
test_plane (*tool);
test_line (*tool);
--
GitLab