diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
index 67dd6513807147a43a6e6ff4b4840d6754c51750..8340f1a2016c8013dd129f1ca57de40150c757dd 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 94224bcc98782982bade1151291f54e1f128c6d3..1ef2b7b75f8379f3ce80c7fc9609e8f5ee4193a0 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 ed9b945fd7666d51c4696af3ae0592c981a261cb..e22da5001735f61591b92f6c498f47418c906d1c 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 21236565d5b8b139c0d055b0afb51cdfeed15012..3dc8e769667d14ca16f5b9fadfd3160c09d83289 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 a39a899db30b6fe10d306662872735892fd36af6..dea5a22d9103975fc2c2eb0d509ee58bd2c23d35 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 305f0173d80e76b6190b54109fe5c47f84100bd0..24348ab7dbe290620b4e3864a7c4c7705419c04a 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 8e972966b1d0fa8becfd881baad5c031308e84d3..b505fe9d4581da7428311dc74e378ae77d673592 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);