diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
new file mode 100755
index 0000000000000000000000000000000000000000..7de8fe2ae0330c200601a38a1a89c0ff3bda6b36
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidParametrization.h
@@ -0,0 +1,312 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+ Fast (approximate) methods for solenoidal field properties
+ ----------------------------------------------------------
+ Copyright (C) 2000 by ATLAS Collaboration
+ ***************************************************************************/
+
+#ifndef TRKEXSOLENOIDALINTERSECTOR_SOLENOIDPARAMETRIZATION_H
+#define TRKEXSOLENOIDALINTERSECTOR_SOLENOIDPARAMETRIZATION_H
+
+//<<<<<< INCLUDES >>>>>>
+
+#include <cassert>
+#include <cmath>
+// #include <iostream>
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+
+//<<<<<< CLASS DECLARATIONS >>>>>>
+
+// class SolenoidParametrization
+//
+// This is a singleton class so that any object can obtain access to the
+// approximate field without repeating the expensive parametrisation method
+//
+
+namespace Trk
+{
+
+class SolenoidParametrization
+{
+
+public:
+ // private constructors (as singleton)
+ ~SolenoidParametrization (void);
+
+ // forbidden copy constructor
+ // forbidden assignment operator
+
+ static SolenoidParametrization* instance (void); // access to singleton
+ 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 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?
+
+private:
+ 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);
+
+ 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;
+ static const double s_maximumImpactAtOrigin;
+ static const double s_maximumZatOrigin;
+ static const int s_numberParameters;
+ static const double s_rInner;
+ 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;
+
+ // copy, assignment: no semantics, no implementation
+ SolenoidParametrization (const SolenoidParametrization&);
+ SolenoidParametrization &operator= (const SolenoidParametrization&);
+};
+
+//<<<<<< INLINE PRIVATE MEMBER FUNCTIONS >>>>>>
+
+inline int
+SolenoidParametrization::fieldKey(void)
+{
+ double z = m_zAtAxis - s_binZeroZ;
+ int zBin = static_cast<int>(s_binInvSizeZ*z);
+ m_interpolateZ = z*s_binInvSizeZ - double(zBin);
+ if (zBin < 0)
+ {
+ m_interpolateZ = 0.;
+ zBin = 0;
+ }
+ else if (zBin > s_maxBinZ - 2)
+ {
+ m_interpolateZ = 0.;
+ zBin = s_maxBinZ - 1;
+ }
+ m_complementZ = 1. - m_interpolateZ;
+
+ int thetaBin = static_cast<int>(s_binInvSizeTheta*m_cotTheta);
+ m_interpolateTheta = m_cotTheta*s_binInvSizeTheta - double(thetaBin);
+ if (thetaBin > s_maxBinTheta - 3)
+ {
+ m_interpolateTheta = 0.;
+ thetaBin = s_maxBinTheta - 2;
+ }
+ m_complementTheta = 1. - m_interpolateTheta;
+ // std::cout << " m_zAtAxis " << m_zAtAxis
+ // << " m_cotTheta " << m_cotTheta
+ // << " zBin " << zBin
+ // << " interpolateZ " << m_interpolateZ
+ // << " thetaBin " << thetaBin
+ // << " interpolateTheta " << m_interpolateTheta << std::endl;
+
+ return 2*s_numberParameters*(s_maxBinTheta*zBin + thetaBin);
+}
+
+inline void
+SolenoidParametrization::integrate(double& firstIntegral,
+ double& secondIntegral,
+ double zBegin,
+ double zEnd) const
+{
+ double zDiff = zEnd - zBegin;
+ double zBeg2 = zBegin*zBegin;
+ double zBeg3 = zBeg2*zBegin;
+ double zEnd2 = zEnd*zEnd;
+ 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;
+ 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);
+}
+
+inline double
+SolenoidParametrization::interpolate(int key1,
+ int key2,
+ int key3,
+ int key4) 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);
+}
+
+inline void
+SolenoidParametrization::setTerms(int key1)
+{
+ int key2 = key1 + s_numberParameters;
+ int key3 = key2 + s_numberParameters;
+ int key4 = key3 + s_numberParameters;
+
+ assert (key1 >= 0);
+ assert (key4 < 14688);
+ assert (s_parameters[key1] != 0.);
+ assert (s_parameters[key3] != 0.);
+ if (m_cotTheta < 7.)
+ {
+ assert (s_parameters[key2] != 0.);
+ assert (s_parameters[key4] != 0.);
+ }
+
+
+ m_fieldAtOrigin = interpolate(key1++,key2++,key3++,key4++);
+ m_quadraticTerm = interpolate(key1++,key2++,key3++,key4++);
+ m_cubicTerm = interpolate(key1,key2,key3,key4);
+}
+
+//<<<<<< INLINE PUBLIC MEMBER FUNCTIONS >>>>>>
+
+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; }
+
+inline double
+SolenoidParametrization::fieldComponent (double z) const
+{
+ double z_local = m_signTheta*z - m_zAtAxis;
+ double z_squared = z_local*z_local;
+ double value = m_fieldAtOrigin +
+ m_quadraticTerm*z_squared +
+ m_cubicTerm*z_squared*z_local;
+ return value;
+}
+
+inline double
+SolenoidParametrization::fieldComponent (double r, double z, double cotTheta) const
+{
+ Amg::Vector3D field;
+ Amg::Vector3D position(r, 0., z);
+ m_magFieldSvc->getField(&position,&field);
+ return s_lightSpeed*(field.z() - field.perp()*cotTheta);
+}
+
+inline void
+SolenoidParametrization::fieldIntegrals(double& firstIntegral,
+ double& secondIntegral,
+ double zBegin,
+ double zEnd)
+{
+ zBegin = m_signTheta*zBegin;
+ zEnd = m_signTheta*zEnd;
+ if (zEnd < s_zInner || zBegin > s_zInner)
+ {
+ integrate(firstIntegral,secondIntegral,zBegin-m_zAtAxis,zEnd-m_zAtAxis);
+ }
+ 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);
+ }
+
+ // std::cout << " zBegin < 0. " << zBegin
+ // << " zEnd " << zEnd
+ // << " m_signTheta " << m_signTheta
+ // << " m_zAtAxis " << m_zAtAxis << std::endl;
+}
+
+inline bool
+SolenoidParametrization::hasParametrization (void) const
+{ return m_hasParametrization; }
+
+inline double
+SolenoidParametrization::maximumR (void) const
+{ return s_rInner; }
+
+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
+{ return (origin.perp() < s_maximumImpactAtOrigin
+ && std::abs(origin.z()) < s_maximumZatOrigin); }
+
+} // end of namespace
+
+#endif // TRKEXSOLENOIDALINTERSECTOR_SOLENOIDPARAMETRIZATION_H
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h
new file mode 100755
index 0000000000000000000000000000000000000000..3af54b96b50cd1be278f761834f160ea5505fbc1
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/TrkExSolenoidalIntersector/SolenoidalIntersector.h
@@ -0,0 +1,283 @@
+/*
+ Copyright (C) 2002-2017 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, 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 Intersection* intersectSurface(const Surface& surface,
+ const Intersection* trackIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : PerigeeSurface */
+ const Intersection* approachPerigeeSurface(const PerigeeSurface& surface,
+ const Intersection* trackIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : StraightLineSurface */
+ const Intersection* approachStraightLineSurface(const StraightLineSurface& surface,
+ const Intersection* trackIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : CylinderSurface */
+ const Intersection* intersectCylinderSurface (const CylinderSurface& surface,
+ const Intersection* trackIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : DiscSurface */
+ const Intersection* intersectDiscSurface (const DiscSurface& surface,
+ const Intersection* trackIntersection,
+ const double qOverP);
+
+ /**IIntersector interface method for specific Surface type : PlaneSurface */
+ const Intersection* intersectPlaneSurface(const PlaneSurface& surface,
+ const Intersection* trackIntersection,
+ 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 Intersection* 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
+ int m_countExtrapolations;
+ int 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 m_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.)
+ {
+ m_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.)
+ {
+ m_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 (m_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)
+{
+ SurfaceIntersection SLIntersect = surface.straightLineIntersection(m_position, m_direction, false, false);
+ if (! SLIntersect.valid) return 0;
+
+ const Intersection* intersection = new Intersection(SLIntersect.intersection,
+ 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 Intersection* trackIntersection, double qOverP)
+{
+ if (trackIntersection->serialNumber() != m_intersectionNumber || qOverP != m_qOverP)
+ {
+ // ATH_MSG_INFO(" initialize parameters. Diff: " << trackIntersection->serialNumber()-m_intersectionNumber
+ // << " at R,Z: " << trackIntersection->position().perp() << ", " << trackIntersection->position().z());
+ ++m_countExtrapolations;
+ m_position.x() = trackIntersection->position().x();
+ m_position.y() = trackIntersection->position().y();
+ m_position.z() = trackIntersection->position().z();
+ m_radius = m_position.perp();
+ m_direction.x() = trackIntersection->direction().x();
+ m_direction.y() = trackIntersection->direction().y();
+ m_direction.z() = trackIntersection->direction().z();
+ m_sinTheta = m_direction.perp();
+ m_oneOverSinTheta = 1./m_sinTheta;
+ m_cotTheta = m_direction.z() * m_oneOverSinTheta;
+ m_pathLength = trackIntersection->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
+
+
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/cmt/requirements b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/cmt/requirements
new file mode 100755
index 0000000000000000000000000000000000000000..0b1443eb39085804e25a6465954865a123ae8a8e
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/cmt/requirements
@@ -0,0 +1,25 @@
+package TrkExSolenoidalIntersector
+
+author Alan Poppleton <Alan.Poppleton@cern.ch>
+
+private
+use EventPrimitives EventPrimitives-* Event
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+
+public
+use AtlasPolicy AtlasPolicy-*
+use AthenaBaseComps AthenaBaseComps-* Control
+use GaudiInterface GaudiInterface-* External
+use GeoPrimitives GeoPrimitives-* DetectorDescription
+use MagFieldInterfaces MagFieldInterfaces-* MagneticField
+use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
+use TrkExUtils TrkExUtils-* Tracking/TrkExtrapolation
+
+library TrkExSolenoidalIntersector SolenoidalIntersector.cxx \
+ SolenoidParametrization.cxx \
+ components/*.cxx
+
+apply_pattern component_library
+
+private
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..7e0c2974b451ef9c1d5307da8d72cf27be4de52e
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidParametrization.cxx
@@ -0,0 +1,476 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/***************************************************************************
+ Fast (approximate) methods for solenoidal field properties
+ ----------------------------------------------------------
+ Copyright (C) 2000 by ATLAS Collaboration
+ ***************************************************************************/
+
+//<<<<<< INCLUDES >>>>>>
+
+#include <algorithm>
+#include <iomanip>
+#include <iostream>
+#include "EventPrimitives/EventPrimitives.h"
+#include "GaudiKernel/SystemOfUnits.h"
+#include "TrkExSolenoidalIntersector/SolenoidParametrization.h"
+
+namespace Trk
+{
+
+//<<<<<< PRIVATE VARIABLE DEFINITIONS >>>>>>
+
+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;
+const double SolenoidParametrization::s_maximumImpactAtOrigin= 30.*Gaudi::Units::mm;
+const double SolenoidParametrization::s_maximumZatOrigin = 250.*Gaudi::Units::mm;
+const int SolenoidParametrization::s_numberParameters = 6;
+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.};
+
+//<<<<<< CLASS STRUCTURE INITIALIZATION >>>>>>
+
+// note: both constructors are private
+SolenoidParametrization::SolenoidParametrization(void)
+ : m_hasParametrization (false)
+{
+ // get central field value in units required for fast tracking
+ s_centralField = fieldComponent(0.,0.,0.);
+
+ // note field = B*c
+ std::cout << " SolenoidParametrization: centralField "
+ << std::setw(6) << std::setprecision(3) << s_centralField/s_lightSpeed /Gaudi::Units::tesla
+ << "T" << std::endl;
+}
+
+SolenoidParametrization::SolenoidParametrization(MagField::IMagFieldSvc* magFieldSvc)
+ : m_hasParametrization (false),
+ m_magFieldSvc (magFieldSvc)
+{
+ // get central field value in units required for fast tracking (i.e. field = B*c)
+ s_centralField = fieldComponent(0.,0.,0.);
+ std::cout << " SolenoidParametrization: centralField "
+ << std::setw(6) << std::setprecision(3) << s_centralField/s_lightSpeed /Gaudi::Units::tesla
+ << "T";
+
+ // now parametrise field - if requested
+ if (magFieldSvc)
+ {
+ std::cout << " - please be patient while the solenoid is parametrised !! ";
+ parametrizeSolenoid();
+ }
+ std::cout << std::endl;
+}
+
+//<<<<<< PRIVATE MEMBER FUNCTION DEFINITIONS >>>>>>
+
+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 ?
+
+ for (int binZ = 0; binZ < s_maxBinZ; ++binZ)
+ {
+ m_cotTheta = smallOffset;
+ for (int binTheta = 0; binTheta < s_maxBinTheta - 1; ++binTheta)
+ {
+ double r = 0.;
+ double z = m_zAtAxis;
+ int n = 200;
+ double dr;
+ if (m_cotTheta < s_zOuter/s_rOuter)
+ {
+ dr = s_rOuter/double(n);
+ }
+ else
+ {
+ dr = s_zOuter/(m_cotTheta*double(n));
+ }
+
+ Amg::VectorX difference(n);
+ Amg::MatrixX derivative = Amg::MatrixX::Zero(n,s_numberParameters);
+ for (int k = 0; k < n; ++k)
+ {
+ r += dr;
+ z += dr*m_cotTheta;
+ double w = (n - k)*(n - k);
+ double zLocal = z - m_zAtAxis;
+ if (r > s_rInner || z > s_zInner)
+ {
+ // derivative(k,0) = 0.;
+ // derivative(k,1) = 0.;
+ // derivative(k,2) = 0.;
+ derivative(k,3) = w;
+ derivative(k,4) = w*zLocal*zLocal;
+ derivative(k,5) = w*zLocal*zLocal*zLocal;
+ }
+ else
+ {
+ derivative(k,0) = w;
+ derivative(k,1) = w*zLocal*zLocal;
+ derivative(k,2) = w*zLocal*zLocal*zLocal;
+ // derivative(k,3) = 0.;
+ // derivative(k,4) = 0.;
+ // derivative(k,5) = 0.;
+ }
+ difference(k) = w*(fieldComponent(r,z,m_cotTheta) - s_centralField);
+ }
+
+ // solve for parametrization coefficients
+ Amg::VectorX solution = derivative.colPivHouseholderQr().solve(difference);
+
+ 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);
+
+ // 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);
+ key -= s_numberParameters;
+ }
+
+ // duplicate next to previous z-bin for contiguous neighbour lookup
+ 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);
+ }
+
+ // some debug print
+// key -= 6;
+// if ( key < s_numberParameters*s_maxBinTheta
+// || key >= s_numberParameters*s_maxBinTheta*(s_maxBinZ - 2))
+// {
+// double z_max;
+// if (m_cotTheta < s_zInner/s_rInner)
+// {
+// z_max = s_rInner*m_cotTheta + m_zAtAxis;
+// }
+// else
+// {
+// z_max = s_zInner;
+// }
+// cout << std::setiosflags(std::ios::fixed) << key
+// << " cotTheta " << std::setw(7) << std::setprecision(2) << m_cotTheta
+// << " inner terms: z0 "<< std::setw(6) << std::setprecision(3)
+// << s_parameters[key]/s_centralField
+// << " z^2 "<< std::setw(6) << std::setprecision(3)
+// << s_parameters[key+1]*z_max*z_max/s_centralField
+// << " z^3 " << std::setw(6) << std::setprecision(3)
+// << s_parameters[key+2]*z_max*z_max*z_max/s_centralField
+// << " outer terms: z0 "<< std::setw(6) << std::setprecision(3)
+// << s_parameters[key+3]/s_centralField
+// << " z^2 "<< std::setw(6) << std::setprecision(3)
+// << s_parameters[key+4]*z_max*z_max/s_centralField
+// << " z^3 " << std::setw(6) << std::setprecision(3)
+// << s_parameters[key+5]*z_max*z_max*z_max/s_centralField
+// << std::resetiosflags(std::ios::fixed) << endl;
+// }
+ m_cotTheta += 1./s_binInvSizeTheta;
+ }
+ m_zAtAxis += 1./s_binInvSizeZ;
+ }
+
+ // duplicate end theta-bins for contiguous neighbour lookup
+ m_zAtAxis = s_binZeroZ; // + smallOffset ??
+ for (int binZ = 0; binZ < s_maxBinZ; ++binZ)
+ {
+ m_cotTheta = double(s_maxBinTheta)/s_binInvSizeTheta;
+ int key = fieldKey();
+ 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];
+ ++key;
+ }
+ m_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
+{
+ // integrate along lines of const eta from origin to r = 1m or |z| = 2.65m
+ // direction normalised st transverse component = 1 (equiv to a fixed pt)
+ std::cout << std::setiosflags(std::ios::fixed)
+ << " eta rEnd mean(Bz) max(dBz/dR) mean(Bt) max(dBt/dR) "
+ << "min(Bt) max(Bt) reverse-bend(z) integrals: Bt.dR Bl.dR"
+ << " asymm: x y z"
+ << " " << std::endl
+ << " m T T/m T T/m "
+ << " T T m T.m T.m"
+ << " T T T"
+ << std::endl;
+
+ double maxR = 1000.*Gaudi::Units::mm;
+ double maxZ = 2650.*Gaudi::Units::mm;
+ int numSteps = 1000;
+
+ // step through eta-range
+ double eta = 0.;
+ for (int i = 0; i != 31; ++i)
+ {
+ double phi = 0.;
+ double cotTheta = sinh(eta);
+ Amg::Vector3D direction(cos(phi),sin(phi),cotTheta);
+ double rEnd = maxR;
+ if (std::abs(cotTheta) > maxZ/maxR) rEnd = maxZ/direction.z();
+ double step = rEnd/static_cast<double>(numSteps); // radial step in mm
+ Amg::Vector3D position(0.,0.,0.);
+ double meanBL = 0.;
+ double meanBT = 0.;
+ double meanBZ = 0.;
+ double maxBT = 0.;
+ double minBT = 9999.;
+ double maxGradBT = 0.;
+ double maxGradBZ = 0.;
+ double prevBT = 0.;
+ double prevBZ = 0.;
+ double reverseZ = 0.;
+
+ // look up field along eta-line
+ for (int j = 0; j != numSteps; ++j)
+ {
+ position += 0.5*step*direction;
+ Amg::Vector3D field;
+ m_magFieldSvc->getField(&position,&field);
+ Amg::Vector3D vCrossB = direction.cross(field);
+ position += 0.5*step*direction;
+ double BZ = field.z();
+ double BT = vCrossB.x()*direction.y() - vCrossB.y()*direction.x();
+ double BL = std::sqrt(vCrossB.mag2() - BT*BT);
+ meanBL += BL;
+ meanBT += BT;
+ meanBZ += BZ;
+ if (BT > maxBT) maxBT = BT;
+ if (BT < minBT) minBT = BT;
+ if (j > 0)
+ {
+ if (BT*prevBT < 0.) reverseZ = position.z() - step*direction.z();
+ double grad = std::abs(BT - prevBT);
+ if (grad > maxGradBT) maxGradBT = grad;
+ grad = std::abs(BZ - prevBZ);
+ if (grad > maxGradBZ) maxGradBZ = grad;
+ }
+ prevBT = BT;
+ prevBZ = BZ;
+ }
+
+ // normalize
+ maxGradBT *= static_cast<double>(numSteps)/step;
+ maxGradBZ *= static_cast<double>(numSteps)/step;
+ meanBL /= static_cast<double>(numSteps);
+ meanBT /= static_cast<double>(numSteps);
+ meanBZ /= static_cast<double>(numSteps);
+ double integralBL = meanBL*rEnd;
+ double integralBT = meanBT*rEnd;
+
+ std::cout << std::setw(6) << std::setprecision(2) << eta
+ << std::setw(8) << std::setprecision(3) << rEnd /Gaudi::Units::meter
+ << std::setw(11) << std::setprecision(4) << meanBZ /Gaudi::Units::tesla
+ << std::setw(12) << std::setprecision(3) << maxGradBZ /Gaudi::Units::tesla
+ << std::setw(13) << std::setprecision(4) << meanBT /Gaudi::Units::tesla
+ << std::setw(12) << std::setprecision(3) << maxGradBT /Gaudi::Units::tesla
+ << std::setw(13) << std::setprecision(4) << minBT /Gaudi::Units::tesla
+ << std::setw(8) << std::setprecision(4) << maxBT /Gaudi::Units::tesla;
+ if (reverseZ > 0.)
+ {
+ std::cout << std::setw(17) << std::setprecision(2) << reverseZ /Gaudi::Units::meter
+ << std::setw(18) << std::setprecision(4) << integralBT /(Gaudi::Units::tesla*Gaudi::Units::meter)
+ << std::setw(8) << std::setprecision(4) << integralBL /(Gaudi::Units::tesla*Gaudi::Units::meter)
+ << " ";
+ }
+ else
+ {
+ std::cout << std::setw(35) << std::setprecision(4) << integralBT /(Gaudi::Units::tesla*Gaudi::Units::meter)
+ << std::setw(8) << std::setprecision(4) << integralBL /(Gaudi::Units::tesla*Gaudi::Units::meter)
+ << " ";
+ }
+
+ // check symmetry (reflect in each axis)
+ for (int k = 0; k != 3; ++k)
+ {
+ if (k == 0) direction = Amg::Vector3D(-cos(phi),sin(phi),cotTheta);
+ if (k == 1) direction = Amg::Vector3D(cos(phi),-sin(phi),cotTheta);
+ if (k == 2) direction = Amg::Vector3D(cos(phi),sin(phi),-cotTheta);
+ position = Amg::Vector3D(0.,0.,0.);
+ double asymm = 0.;
+ // look up field along eta-line
+ for (int j = 0; j != numSteps; ++j)
+ {
+ position += 0.5*step*direction;
+ Amg::Vector3D field;
+ m_magFieldSvc->getField(&position,&field);
+ Amg::Vector3D vCrossB = direction.cross(field);
+ position += 0.5*step*direction;
+ double BT = vCrossB.x()*direction.y() - vCrossB.y()*direction.x();
+ asymm += BT;
+ }
+ asymm = asymm/static_cast<double>(numSteps) - meanBT;
+ std::cout << std::setw(9) << std::setprecision(4) << asymm /Gaudi::Units::tesla;
+ }
+ std::cout << std::endl;
+ eta += 0.1;
+ }
+}
+
+void
+SolenoidParametrization::printParametersForEtaLine (double eta, double z_origin)
+{
+ m_signTheta = 1.;
+ m_zAtAxis = z_origin;
+ m_cotTheta = 1./std::tan(2.*std::atan(1./std::exp(eta)));
+ int key = fieldKey();
+ double z_max;
+ if (m_cotTheta < s_zInner/s_rInner)
+ {
+ z_max = s_rInner*m_cotTheta;
+ }
+ else
+ {
+ z_max = s_zInner;
+ }
+ std::cout << std::setiosflags(std::ios::fixed)
+ << "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
+ << " z^2 "<< std::setw(6) << std::setprecision(3)
+ << s_parameters[key+1]*z_max*z_max/s_centralField
+ << " z^3 " << std::setw(6) << std::setprecision(3)
+ << s_parameters[key+2]*z_max*z_max*z_max/s_centralField
+ << " outer terms: z0 "<< std::setw(6) << std::setprecision(3)
+ << s_parameters[key+3]/s_centralField
+ << " z^2 "<< std::setw(6) << std::setprecision(3)
+ << s_parameters[key+4]*z_max*z_max/s_centralField
+ << " z^3 " << std::setw(6) << std::setprecision(3)
+ << s_parameters[key+5]*z_max*z_max*z_max/s_centralField
+ << std::resetiosflags(std::ios::fixed) << std::endl;
+}
+
+void
+SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin)
+{
+ 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.;
+ int n = 200;
+ double dr;
+ if (cotTheta < s_zOuter/s_rOuter)
+ {
+ dr = s_rOuter/double(n);
+ }
+ else
+ {
+ dr = s_zOuter/(cotTheta*double(n));
+ }
+ double chiSquareIn = 0.;
+ double chiSquareOut = 0.;
+ double nIn = 0.;
+ double nOut = 0.;
+ double worstBCalc = 0.;
+ double worstBTrue = 0.;
+ double worstDiff = -1.;
+ double worstR = 0.;
+ double worstZ = 0.;
+ for (int k = 0; k < n; ++k)
+ {
+ double b = fieldComponent(r,z,cotTheta);
+ setParameters(r,z,cotTheta);
+ double diff = (fieldComponent(z) - b)/s_lightSpeed;
+
+ if (r > s_rInner || z > s_zInner)
+ {
+ chiSquareOut+= diff*diff;
+ nOut += 1.;
+ }
+ else
+ {
+ chiSquareIn += diff*diff;
+ nIn += 1.;
+ }
+
+ if (std::abs(diff) > worstDiff)
+ {
+ worstDiff = std::abs(diff);
+ worstBCalc = fieldComponent(z);
+ worstBTrue = b;
+ worstR = r;
+ worstZ = z;
+ }
+ r += dr;
+ z += dr*cotTheta;
+ }
+
+ std::cout << std::setiosflags(std::ios::fixed)
+ << "SolenoidParametrization: line with eta " << std::setw(6) << std::setprecision(2) << eta
+ << " from (r,z) 0.0, " << std::setw(6) << std::setprecision(1) << zOrigin
+ << " rms diff inner/outer " << std::setw(6) << std::setprecision(3)
+ << std::sqrt(chiSquareIn/nIn) /Gaudi::Units::tesla << " " << std::setw(6) << std::setprecision(3)
+ << std::sqrt(chiSquareOut/nOut) /Gaudi::Units::tesla << std::endl
+ << " worst residual at: (r,z) "
+ << std::setw(6) << std::setprecision(1) << worstR
+ << ", " << std::setw(6) << std::setprecision(1) << worstZ
+ << " with B true/calc " << std::setw(6) << std::setprecision(3)
+ << worstBTrue/s_lightSpeed /Gaudi::Units::tesla
+ << " " << std::setw(6) << std::setprecision(3) << worstBCalc/s_lightSpeed /Gaudi::Units::tesla
+ << std::resetiosflags(std::ios::fixed) << std::endl;
+}
+
+} // end of namespace
+
+
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..c4322033ec6e65637b5310fa70495fe94973c670
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/SolenoidalIntersector.cxx
@@ -0,0 +1,466 @@
+/*
+ Copyright (C) 2002-2017 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 Intersection* 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 Intersection* 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 Intersection* 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 Intersection* 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 Intersection* 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 Intersection* 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
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_entries.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_entries.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..48593d54165849334d2076dd9226c83e652a85f8
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_entries.cxx
@@ -0,0 +1,10 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "TrkExSolenoidalIntersector/SolenoidalIntersector.h"
+
+typedef Trk::SolenoidalIntersector TrkSolenoidalIntersector;
+DECLARE_TOOL_FACTORY( TrkSolenoidalIntersector )
+
+DECLARE_FACTORY_ENTRIES( TrkExSolenoidalIntersector )
+{
+ DECLARE_NAMESPACE_TOOL( Trk, SolenoidalIntersector )
+}
diff --git a/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_load.cxx b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_load.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..97fdfcbf98d2c4ae9084a77530a2f41831e15f13
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExSolenoidalIntersector/src/components/TrkExSolenoidalIntersector_load.cxx
@@ -0,0 +1,4 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( TrkExSolenoidalIntersector )
+