diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/TrkExRungeKuttaPropagator/RungeKuttaPropagator.h b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/TrkExRungeKuttaPropagator/RungeKuttaPropagator.h
new file mode 100755
index 0000000000000000000000000000000000000000..67e42103f4ea1151b27ec703db128623f560f8ad
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/TrkExRungeKuttaPropagator/RungeKuttaPropagator.h
@@ -0,0 +1,445 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/////////////////////////////////////////////////////////////////////////////////
+// Header file for class RungeKuttaPropagator, (c) ATLAS Detector software
+/////////////////////////////////////////////////////////////////////////////////
+
+#ifndef RungeKuttaPropagator_H
+#define RungeKuttaPropagator_H
+
+#include <list>
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ServiceHandle.h"
+#include "MagFieldInterfaces/IMagFieldSvc.h"
+#include "TrkExInterfaces/IPropagator.h"
+#include "TrkExInterfaces/IPatternParametersPropagator.h"
+#include "TrkEventPrimitives/ParticleHypothesis.h"
+#include "TrkEventPrimitives/PropDirection.h"
+#include "TrkParameters/TrackParameters.h"
+#include "TrkNeutralParameters/NeutralParameters.h"
+
+namespace Trk {
+
+ class Surface ;
+ class MagneticFieldProperties ;
+ class CylinderBounds ;
+ class CylinderSurface ;
+ class PatternTrackParameters ;
+
+ /**
+ @class RungeKuttaPropagator
+
+ Trk::RungeKuttaPropagator is algorithm for track parameters propagation through
+ magnetic field with or without jacobian of transformation. This algorithm
+ contains three steps.
+
+ 1.The first step of the algorithm is track parameters transformation from
+ local presentation for given start surface to global Runge Kutta coordinates.
+
+ 2.The second step is propagation through magnetic field with or without jacobian.
+
+ 3.Third step is transformation from global Runge Kutta presentation to local
+ presentation of given output surface.
+
+
+ AtaPlane AtaStraightLine AtaDisc AtaCylinder Perigee
+ | | | | |
+ | | | | |
+ V V V V V
+ -----------------------------------------------------------------
+ | Local->Global transformation
+ V
+ Global position (Runge Kutta presentation)
+ |
+ |
+ Propagation to next surface with or without jacobian
+ using Nystroem algorithm
+ (See Handbook Net. Bur. of Standards, procedure 25.5.20)
+ |
+ V Global->Local transformation
+ ----------------------------------------------------------------
+ | | | | |
+ | | | | |
+ V V V V V
+ PlaneSurface StraightLineSurface DiscSurface CylinderSurface PerigeeSurface
+
+ For propagation using Runge Kutta method we use global coordinate, direction,
+ inverse momentum and Jacobian of transformation. All this parameters we save
+ in array P[42].
+ /dL0 /dL1 /dPhi /dThe /dCM
+ X ->P[0] dX / P[ 7] P[14] P[21] P[28] P[35]
+ Y ->P[1] dY / P[ 8] P[15] P[22] P[29] P[36]
+ Z ->P[2] dZ / P[ 9] P[16] P[23] P[30] P[37]
+ Ax ->P[3] dAx/ P[10] P[17] P[24] P[31] P[38]
+ Ay ->P[4] dAy/ P[11] P[18] P[25] P[32] P[39]
+ Az ->P[5] dAz/ P[12] P[19] P[26] P[33] P[40]
+ CM ->P[6] dCM/ P[13] P[20] P[27] P[34] P[41]
+
+ where
+ in case local presentation
+
+ L0 - first local coordinate (surface dependent)
+ L1 - second local coordinate (surface dependent)
+ Phi - Azimuthal angle
+ The - Polar angle
+ CM - charge/momentum
+
+ in case global presentation
+
+ X - global x-coordinate = surface dependent
+ Y - global y-coordinate = surface dependent
+ Z - global z-coordinate = sutface dependent
+ Ax - direction cosine to x-axis = Sin(The)*Cos(Phi)
+ Ay - direction cosine to y-axis = Sin(The)*Sin(Phi)
+ Az - direction cosine to z-axis = Cos(The)
+ CM - charge/momentum = local CM
+
+ Comment:
+ if pointer to const * = 0 algorithm will propagate track
+ parameters and jacobian of transformation according straight line model
+
+
+ @author Igor.Gavrilenko@cern.ch
+ */
+
+ class RungeKuttaPropagator : public AthAlgTool, virtual public IPropagator,
+ virtual public IPatternParametersPropagator
+ {
+ /////////////////////////////////////////////////////////////////////////////////
+ // Public methods:
+ /////////////////////////////////////////////////////////////////////////////////
+
+ public:
+
+ RungeKuttaPropagator(const std::string&,const std::string&,const IInterface*);
+
+ virtual ~RungeKuttaPropagator ();
+
+ /** AlgTool initailize method.*/
+ StatusCode initialize();
+
+ /** AlgTool finalize method */
+ StatusCode finalize();
+
+ /** Main propagation mehtod NeutralParameters */
+
+ const NeutralParameters* propagate
+ (const NeutralParameters &,
+ const Surface &,
+ PropDirection ,
+ BoundaryCheck ,
+ bool ) const;
+
+ /** Main propagation mehtod without transport jacobian production*/
+
+ const TrackParameters* propagate
+ (const TrackParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ const MagneticFieldProperties &,
+ ParticleHypothesis ,
+ bool ,
+ const TrackingVolume* ) const;
+
+ /** Main propagation mehtod with transport jacobian production*/
+
+ const TrackParameters* propagate
+ (const TrackParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ const MagneticFieldProperties &,
+ TransportJacobian *&,
+ ParticleHypothesis ,
+ bool ,
+ const TrackingVolume* ) const;
+
+ /** The propagation method finds the closest surface */
+
+ const TrackParameters* propagate
+ (const TrackParameters &,
+ std::vector<DestSurf> &,
+ PropDirection ,
+ const MagneticFieldProperties &,
+ ParticleHypothesis ,
+ std::vector<unsigned int> &,
+ double &,
+ bool ,
+ bool ,
+ const TrackingVolume* ) const;
+
+ /** Main propagation mehtod for parameters only without transport jacobian productio*/
+
+ const TrackParameters* propagateParameters
+ (const TrackParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ const MagneticFieldProperties &,
+ ParticleHypothesis ,
+ bool ,
+ const TrackingVolume* ) const;
+
+
+ /** Main propagation mehtod for parameters only with transport jacobian productio*/
+
+ const TrackParameters* propagateParameters
+ (const TrackParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ const MagneticFieldProperties &,
+ TransportJacobian *&,
+ ParticleHypothesis ,
+ bool ,
+ const TrackingVolume* ) const;
+
+ /** Global position together with direction of the trajectory on the surface */
+
+ const IntersectionSolution* intersect
+ (const TrackParameters &,
+ const Surface &,
+ const MagneticFieldProperties &,
+ ParticleHypothesis particle=pion,
+ const TrackingVolume* tvol=0 ) const;
+
+ /** GlobalPositions list interface:*/
+
+ void globalPositions
+ (std::list<Amg::Vector3D> &,
+ const TrackParameters &,
+ const MagneticFieldProperties &,
+ const CylinderBounds& ,
+ double ,
+ ParticleHypothesis particle=pion,
+ const TrackingVolume* tvol=0 ) const;
+
+ /** Test quality Jacobian calculation */
+
+ void JacobianTest
+ (const TrackParameters &,
+ const Surface &,
+ const MagneticFieldProperties&) const;
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Public methods for Trk::PatternTrackParameters (from IPattern'Propagator)
+ /////////////////////////////////////////////////////////////////////////////////
+
+ /** Main propagation mehtod */
+
+ bool propagate
+ (PatternTrackParameters &,
+ const Surface &,
+ PatternTrackParameters &,
+ PropDirection ,
+ const MagneticFieldProperties &,
+ ParticleHypothesis particle=pion) const;
+
+ /** Main propagation mehtod with step to surface calculation*/
+
+ bool propagate
+ (PatternTrackParameters &,
+ const Surface &,
+ PatternTrackParameters &,
+ PropDirection ,
+ const MagneticFieldProperties &,
+ double &,
+ ParticleHypothesis particle=pion) const;
+
+ /** Main propagation mehtod for parameters only */
+
+ bool propagateParameters
+ (PatternTrackParameters &,
+ const Surface &,
+ PatternTrackParameters &,
+ PropDirection ,
+ const MagneticFieldProperties &,
+ ParticleHypothesis particle=pion) const;
+
+ /** Main propagation mehtod for parameters only with step to surface calculation*/
+
+ bool propagateParameters
+ (PatternTrackParameters &,
+ const Surface &,
+ PatternTrackParameters &,
+ PropDirection ,
+ const MagneticFieldProperties &,
+ double &,
+ ParticleHypothesis particle=pion) const;
+
+ /** GlobalPositions list interface:*/
+
+ void globalPositions
+ (std::list<Amg::Vector3D> &,
+ const PatternTrackParameters &,
+ const MagneticFieldProperties &,
+ const CylinderBounds &,
+ double ,
+ ParticleHypothesis particle=pion) const;
+
+ /** GlobalPostions and steps for set surfaces */
+
+ void globalPositions
+ (const PatternTrackParameters &,
+ std::list<const Surface*> &,
+ std::list< std::pair<Amg::Vector3D,double> > &,
+ const MagneticFieldProperties &,
+ ParticleHypothesis particle=pion ) const;
+
+ private:
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Private methods:
+ /////////////////////////////////////////////////////////////////////////////////
+
+ /** Internal RungeKutta propagation method for charge track parameters*/
+
+ const TrackParameters* propagateRungeKutta
+ (bool ,
+ const TrackParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ const MagneticFieldProperties&,
+ double *,
+ bool returnCurv) const;
+
+ /** Internal RungeKutta propagation method for neutral track parameters*/
+
+ const NeutralParameters* propagateStraightLine
+ (bool ,
+ const NeutralParameters &,
+ const Surface &,
+ const PropDirection ,
+ BoundaryCheck ,
+ double *,
+ bool returnCurv) const;
+
+ /** Internal RungeKutta propagation method */
+
+ bool propagateRungeKutta
+ (bool ,
+ PatternTrackParameters &,
+ const Surface &,
+ PatternTrackParameters &,
+ PropDirection ,
+ const MagneticFieldProperties&,
+ double &) const;
+
+ /** Internal RungeKutta propagation method for propation with jacobian*/
+
+ bool propagateWithJacobian
+ (bool ,
+ int ,
+ double *,
+ double *,
+ double &) const;
+
+ /** Propagation methods runge kutta step*/
+
+ double rungeKuttaStep
+ (bool ,
+ double ,
+ double *,
+ bool &) const;
+
+ /** Propagation methods runge kutta step*/
+
+ double rungeKuttaStepWithGradient
+ (double ,
+ double *,
+ bool &) const;
+
+ /** Propagation methods straight line step*/
+
+ double straightLineStep
+ (bool ,
+ double ,
+ double*) const;
+
+ /** Test new propagation to cylinder boundary */
+
+ bool newCrossPoint
+ (const CylinderSurface&,
+ const double *,
+ const double *) const;
+
+ /** Step estimator with directions correction */
+
+ double stepEstimatorWithCurvature(int,double*,const double*,bool&) const;
+
+ /** Step reduction */
+ double stepReduction(const double*) const;
+
+ /** Build new track parameters without propagation */
+
+ const TrackParameters* buildTrackParametersWithoutPropagation
+ (const TrackParameters &,double*) const;
+
+ const NeutralParameters* buildTrackParametersWithoutPropagation
+ (const NeutralParameters&,double*) const;
+
+ void globalOneSidePositions
+ (std::list<Amg::Vector3D> &,
+ const double *,
+ const MagneticFieldProperties &,
+ const CylinderBounds &,
+ double ,
+ ParticleHypothesis particle=pion) const;
+
+ void globalTwoSidePositions
+ (std::list<Amg::Vector3D> &,
+ const double *,
+ const MagneticFieldProperties &,
+ const CylinderBounds &,
+ double ,
+ ParticleHypothesis particle=pion) const;
+
+ const Trk::TrackParameters* crossPoint
+ (const TrackParameters &,
+ std::vector<DestSurf> &,
+ std::vector<unsigned int>&,
+ double *,
+ std::pair<double,int> &) const;
+
+ void getField (double*,double* ) const;
+ void getFieldGradient(double*,double*,double*) const;
+
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Private data members:
+ /////////////////////////////////////////////////////////////////////////////////
+
+ double m_dlt ; // accuracy parameter
+ double m_helixStep ; // max step whith helix model
+ double m_straightStep ; // max step whith srtaight line model
+ bool m_usegradient ; // use magnetif field gradient
+ mutable double m_direction ;
+ mutable bool m_mcondition ;
+ mutable bool m_solenoid ;
+ mutable bool m_needgradient ;
+ ServiceHandle<MagField::IMagFieldSvc> m_fieldServiceHandle;
+ MagField::IMagFieldSvc* m_fieldService ;
+ };
+
+ /////////////////////////////////////////////////////////////////////////////////
+ // Inline methods for magnetic field information
+ /////////////////////////////////////////////////////////////////////////////////
+
+ inline void RungeKuttaPropagator::getField(double* R,double* H) const
+ {
+ if(m_solenoid) return m_fieldService->getFieldZR(R,H);
+ return m_fieldService->getField (R,H);
+ }
+
+ inline void RungeKuttaPropagator::getFieldGradient(double* R,double* H,double* dH) const
+ {
+ if(m_solenoid) return m_fieldService->getFieldZR(R,H,dH);
+ return m_fieldService->getField (R,H,dH);
+ }
+}
+
+#endif // RungeKuttaPropagator_H
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/cmt/requirements b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/cmt/requirements
new file mode 100755
index 0000000000000000000000000000000000000000..7719e0da1877d4c7dd7960c85926c3de2591873c
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/cmt/requirements
@@ -0,0 +1,26 @@
+package TrkExRungeKuttaPropagator
+
+author Igor Gavrilenko <Igor.Gavrilenko@cern.ch>
+manager Igor Gavrilenko <Igor.Gavrilenko@cern.ch>
+
+private
+
+use TrkPatternParameters TrkPatternParameters-* Tracking/TrkEvent
+use TrkExUtils TrkExUtils-* Tracking/TrkExtrapolation
+use TrkSurfaces TrkSurfaces-* Tracking/TrkDetDescr
+use TrkGeometry TrkGeometry-* Tracking/TrkDetDescr
+
+public
+
+use GaudiInterface GaudiInterface-* External
+use AtlasPolicy AtlasPolicy-*
+use MagFieldInterfaces * MagneticField
+use AthenaBaseComps AthenaBaseComps-* Control
+use TrkExInterfaces TrkExInterfaces-* Tracking/TrkExtrapolation
+use TrkEventPrimitives TrkEventPrimitives-* Tracking/TrkEvent
+use TrkParameters TrkParameters-* Tracking/TrkEvent
+use TrkNeutralParameters TrkNeutralParameters-* Tracking/TrkEvent
+
+library TrkExRungeKuttaPropagator *.cxx components/*.cxx
+apply_pattern component_library
+
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/doc/mainpage.h b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/doc/mainpage.h
new file mode 100755
index 0000000000000000000000000000000000000000..e5bcf45bbfe79fd0763eeb1a427ff8b9a9ee7a5c
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/doc/mainpage.h
@@ -0,0 +1,29 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+@mainpage The TrkExRungeKutta package
+
+Propagation of TrackParameters and their associated covariance matrices in a Runge-Kutta integrated method
+
+@section TrkExRkIntroduction Introduction
+
+The RungeKuttaPropagator implements the IPropagator interface with a Runge-Kutta
+integrated method to be used in a non-homogenous magnetic field.
+In case of absence of a magnetic field the RungeKuttaPropagator would propagate the track-parameters
+according to a straight line, in case of an homogenous magnetic field, a helical propagation is performed.
+
+Common transformations with the STEP_Propagator can be found in the TrkExUtils package.
+
+The tool implements two abstract interfaces:
+- under Trk::IPropagator it provides propagation of EDM Trk::TrackParameters
+- under Trk::IPatternParametersPropagator it provides propagation of parameters for internal use, the Trk::PatternTrackParameters
+
+@section TrkExRkComments Comments
+
+Please let me know of any errors, or if anything is unclear.
+@author Igor.Gavrilenko@cern.ch
+
+*/
+
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..9a01e5d43576577b10743648eab5f52e460e4b05
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/RungeKuttaPropagator.cxx
@@ -0,0 +1,1729 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/////////////////////////////////////////////////////////////////////////////////
+// RungeKuttaPropagator.cxx, (c) ATLAS Detector software
+/////////////////////////////////////////////////////////////////////////////////
+
+#include "TrkExUtils/RungeKuttaUtils.h"
+#include "TrkExRungeKuttaPropagator/RungeKuttaPropagator.h"
+#include "TrkSurfaces/ConeSurface.h"
+#include "TrkSurfaces/DiscSurface.h"
+#include "TrkSurfaces/PlaneSurface.h"
+#include "TrkSurfaces/PerigeeSurface.h"
+#include "TrkSurfaces/CylinderSurface.h"
+#include "TrkSurfaces/StraightLineSurface.h"
+#include "TrkGeometry/MagneticFieldProperties.h"
+#include "TrkExUtils/IntersectionSolution.h"
+#include "TrkExUtils/TransportJacobian.h"
+#include "TrkPatternParameters/PatternTrackParameters.h"
+
+/////////////////////////////////////////////////////////////////////////////////
+// Constructor
+/////////////////////////////////////////////////////////////////////////////////
+
+Trk::RungeKuttaPropagator::RungeKuttaPropagator
+(const std::string& p,const std::string& n,const IInterface* t) :
+ AthAlgTool(p,n,t),
+ m_fieldServiceHandle("AtlasFieldSvc",n)
+{
+ m_dlt = .000200;
+ m_helixStep = 1. ;
+ m_straightStep = .01 ;
+ m_usegradient = false ;
+ m_fieldService = 0 ;
+ m_solenoid = true ;
+
+ declareInterface<Trk::IPropagator>(this);
+ declareInterface<Trk::IPatternParametersPropagator>(this);
+ declareProperty("AccuracyParameter" ,m_dlt );
+ declareProperty("MaxHelixStep" ,m_helixStep );
+ declareProperty("MaxStraightLineStep", m_straightStep);
+ declareProperty("IncludeBgradients" , m_usegradient );
+ declareProperty("MagFieldSvc" , m_fieldServiceHandle);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// initialize
+/////////////////////////////////////////////////////////////////////////////////
+
+StatusCode Trk::RungeKuttaPropagator::initialize()
+{
+ if( !m_fieldServiceHandle.retrieve() ){
+ ATH_MSG_FATAL("Failed to retrieve " << m_fieldServiceHandle );
+ return StatusCode::FAILURE;
+ }
+ ATH_MSG_DEBUG("Retrieved " << m_fieldServiceHandle );
+ m_fieldService = &*m_fieldServiceHandle;
+
+ msg(MSG::INFO) << name() <<" initialize() successful" << endreq;
+ return StatusCode::SUCCESS;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// finalize
+/////////////////////////////////////////////////////////////////////////////////
+
+StatusCode Trk::RungeKuttaPropagator::finalize()
+{
+ msg(MSG::INFO) << name() <<" finalize() successful" << endreq;
+ return StatusCode::SUCCESS;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Destructor
+/////////////////////////////////////////////////////////////////////////////////
+
+Trk::RungeKuttaPropagator::~RungeKuttaPropagator(){}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for NeutralParameters propagation
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::NeutralParameters* Trk::RungeKuttaPropagator::propagate
+(const Trk::NeutralParameters & Tp,
+ const Trk::Surface & Su,
+ Trk::PropDirection D ,
+ Trk::BoundaryCheck B ,
+ bool returnCurv) const
+{
+ double J[25];
+ return propagateStraightLine(true,Tp,Su,D,B,J,returnCurv);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for track parameters and covariance matrix propagation
+// without transport Jacobian production
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagate
+(const Trk::TrackParameters & Tp,
+ const Trk::Surface & Su,
+ Trk::PropDirection D,
+ Trk::BoundaryCheck B,
+ const MagneticFieldProperties& M,
+ ParticleHypothesis ,
+ bool returnCurv,
+ const TrackingVolume* ) const
+{
+ double J[25];
+ return propagateRungeKutta(true,Tp,Su,D,B,M,J,returnCurv);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for track parameters and covariance matrix propagation
+// with transport Jacobian production
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagate
+(const Trk::TrackParameters & Tp ,
+ const Trk::Surface& Su ,
+ Trk::PropDirection D ,
+ Trk::BoundaryCheck B ,
+ const MagneticFieldProperties& M ,
+ TransportJacobian *& Jac,
+ ParticleHypothesis ,
+ bool returnCurv,
+ const TrackingVolume* ) const
+{
+ double J[25];
+
+ const Trk::TrackParameters* Tpn = propagateRungeKutta(true,Tp,Su,D,B,M,J,returnCurv);
+
+ if(Tpn) {
+ J[24]=J[20]; J[23]=0.; J[22]=0.; J[21]=0.; J[20]=0.;
+ Jac = new Trk::TransportJacobian(J);
+ }
+ else Jac = 0;
+ return Tpn;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function to finds the closest surface
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagate
+(const TrackParameters & Tp ,
+ std::vector<DestSurf> & DS ,
+ PropDirection D ,
+ const MagneticFieldProperties& M ,
+ ParticleHypothesis ,
+ std::vector<unsigned int> & Sol ,
+ double & Path,
+ bool usePathLim,
+ bool ,
+ const TrackingVolume* ) const
+{
+ Sol.erase(Sol.begin(),Sol.end()); Path = 0.; if(DS.empty()) return 0;
+ m_direction = D;
+
+ // Test is it measured track parameters
+ //
+ bool useJac; Tp.covariance() ? useJac = true : useJac = false;
+
+ // Magnetic field information preparation
+ //
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ (useJac && m_usegradient) ? m_needgradient = true : m_needgradient = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ // Transform to global presentation
+ //
+ Trk::RungeKuttaUtils utils;
+
+ double Po[45],Pn[45];
+ if(!utils.transformLocalToGlobal(useJac,Tp,Po)) return 0; Po[42]=Po[43]=Po[44]=0.;
+
+ // Straight line track propagation for small step
+ //
+ if(D!=0) {
+ double S= m_straightStep; if(D < 0) S = -S; S = straightLineStep(useJac,S,Po);
+ }
+
+ double Wmax = 50000. ; // Max pass
+ double W = 0. ; // Current pass
+ double Smax = 100. ; // Max step
+ if(D < 0) Smax = -Smax;
+ if(usePathLim) Wmax = fabs(Path);
+
+ std::multimap<double,int> DN; double Scut[3];
+ int Nveto = utils.fillDistancesMap(DS,DN,Po,W,&Tp.associatedSurface(),Scut);
+
+ // Test conditions tor start propagation and chocse direction if D == 0
+ //
+ if(DN.empty()) return 0;
+
+ if(D == 0 && fabs(Scut[0]) < fabs(Scut[1])) Smax = -Smax;
+
+ if(Smax < 0. && Scut[0] > Smax) Smax = Scut[0];
+ if(Smax > 0. && Scut[1] < Smax) Smax = Scut[1];
+ if(Wmax > 3.*Scut[2] ) Wmax = 3.*Scut[2];
+
+ double Sl = Smax ;
+ double St = Smax ;
+ bool InS = false;
+ const TrackParameters* To = 0 ;
+
+ for(int i=0; i!=45; ++i) Pn[i]=Po[i];
+
+ //----------------------------------Niels van Eldik patch
+ double last_St = 0. ;
+ bool last_InS = !InS ;
+ bool reverted_P = false;
+ //----------------------------------
+ while (fabs(W) < Wmax) {
+
+ std::pair<double,int> SN;
+ double S;
+
+ if(m_mcondition) {
+
+ //----------------------------------Niels van Eldik patch
+ if (reverted_P && St == last_St && InS == last_InS /*&& condition_fulfiled*/) {
+ // inputs are not changed will get same result.
+ break;
+ }
+ last_St = St;
+ last_InS = InS;
+ //----------------------------------
+
+ if(!m_needgradient) S = rungeKuttaStep (useJac,St,Pn,InS);
+ else S = rungeKuttaStepWithGradient( St,Pn,InS);
+ }
+ else {
+
+ //----------------------------------Niels van Eldik patch
+ if (reverted_P && St == last_St /*&& !condition_fulfiled*/) {
+ // inputs are not changed will get same result.
+ break;
+ }
+ last_St = St;
+ last_InS = InS;
+ //----------------------------------
+
+ S = straightLineStep(useJac,St,Pn);
+ }
+ //----------------------------------Niels van Eldik patch
+ reverted_P=false;
+ //----------------------------------
+
+ bool next; SN=utils.stepEstimator(DS,DN,Po,Pn,W,m_straightStep,Nveto,next);
+
+ if(next) {for(int i=0; i!=45; ++i) Po[i]=Pn[i]; W+=S; Nveto=-1; }
+ else {for(int i=0; i!=45; ++i) Pn[i]=Po[i]; reverted_P=true;}
+
+ if (fabs(S)+1. < fabs(St)) Sl=S;
+ InS ? St = 2.*S : St = S;
+
+ if(SN.second >= 0) {
+
+ double Sa = fabs(SN.first);
+
+ if(Sa > m_straightStep) {
+ if(fabs(St) > Sa) St = SN.first;
+ }
+ else {
+ Path = W+SN.first;
+ if((To = crossPoint(Tp,DS,Sol,Pn,SN))) return To;
+ Nveto = SN.second; St = Sl;
+ }
+ }
+ }
+ return 0;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for track parameters propagation without covariance matrix
+// without transport Jacobian production
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagateParameters
+(const Trk::TrackParameters & Tp,
+ const Trk::Surface & Su,
+ Trk::PropDirection D,
+ Trk::BoundaryCheck B,
+ const MagneticFieldProperties& M,
+ ParticleHypothesis ,
+ bool returnCurv,
+ const TrackingVolume* ) const
+{
+ double J[25];
+ return propagateRungeKutta(false,Tp,Su,D,B,M,J,returnCurv);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for track parameters propagation without covariance matrix
+// with transport Jacobian production
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagateParameters
+(const Trk::TrackParameters & Tp ,
+ const Trk::Surface & Su ,
+ Trk::PropDirection D ,
+ Trk::BoundaryCheck B ,
+ const MagneticFieldProperties& M ,
+ TransportJacobian *& Jac,
+ ParticleHypothesis ,
+ bool returnCurv,
+ const TrackingVolume* ) const
+{
+ double J[25];
+ const Trk::TrackParameters* Tpn = propagateRungeKutta (true,Tp,Su,D,B,M,J,returnCurv);
+
+ if(Tpn) {
+ J[24]=J[20]; J[23]=0.; J[22]=0.; J[21]=0.; J[20]=0.;
+ Jac = new Trk::TransportJacobian(J);
+ }
+ else Jac = 0;
+ return Tpn;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for neutral track parameters propagation with or without jacobian
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::NeutralParameters* Trk::RungeKuttaPropagator::propagateStraightLine
+(bool useJac,
+ const Trk::NeutralParameters & Tp ,
+ const Trk::Surface & Su ,
+ Trk::PropDirection D ,
+ Trk::BoundaryCheck B ,
+ double * Jac ,
+ bool returnCurv) const
+{
+ const Trk::Surface* su = &Su;
+ if(su == &Tp.associatedSurface()) return buildTrackParametersWithoutPropagation(Tp,Jac);
+
+ m_direction = D ;
+ m_mcondition = false;
+
+
+ Trk::RungeKuttaUtils utils;
+
+ double P[64],Step = 0; if(!utils.transformLocalToGlobal(useJac,Tp,P)) return 0;
+
+ const Amg::Transform3D& T = Su.transform();
+ int ty = Su.type();
+
+ if (ty == Trk::Surface::Plane ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Line ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Disc ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cylinder ) {
+
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(su);
+
+ double r0[3] = {P[0],P[1],P[2]};
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),m_direction,0.};
+
+ if(!propagateWithJacobian(useJac,2,s,P,Step)) return 0;
+
+ // For cylinder we do test for next cross point
+ //
+ if(cyl->bounds().halfPhiSector() < 3.1 && newCrossPoint(*cyl,r0,P)) {
+ s[8] = 0.; if(!propagateWithJacobian(useJac,2,s,P,Step)) return 0;
+ }
+ }
+ else if (ty == Trk::Surface::Perigee ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cone ) {
+
+ double k = static_cast<const Trk::ConeSurface*>(su)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,m_direction,0.};
+ if(!propagateWithJacobian(useJac,3,s,P,Step)) return 0;
+ }
+ else return 0;
+
+ if(m_direction!=0. && (m_direction*Step)<0.) return 0;
+
+ // Common transformation for all surfaces (angles and momentum)
+ //
+ if(useJac) {
+ double p=1./P[6]; P[35]*=p; P[36]*=p; P[37]*=p; P[38]*=p; P[39]*=p; P[40]*=p;
+ }
+
+ bool uJ = useJac; if(returnCurv) uJ = false;
+
+ double p[5]; utils.transformGlobalToLocal(su,uJ,P,p,Jac);
+
+ if(B) {Amg::Vector2D L(p[0],p[1]); if(!Su.insideBounds(L,0.)) return 0;}
+
+
+ // Transformation to curvilinear presentation
+ //
+ if(returnCurv) utils.transformGlobalToCurvilinear(useJac,P,p,Jac);
+
+
+ if(!useJac || !Tp.covariance()) {
+
+ if(!returnCurv) {
+ return Su.createNeutralParameters(p[0],p[1],p[2],p[3],p[4],0);
+ }
+ else {
+ Amg::Vector3D gp(P[0],P[1],P[2]);
+ return new Trk::NeutralCurvilinearParameters(gp,p[2],p[3],p[4]);
+ }
+ }
+
+ AmgSymMatrix(5)* e = utils.newCovarianceMatrix(Jac,*Tp.covariance());
+ AmgSymMatrix(5)& cv = *e;
+
+ if(cv(0,0)<=0. || cv(1,1)<=0. || cv(2,2)<=0. || cv(3,3)<=0. || cv(4,4)<=0.) {
+ delete e; return 0;
+ }
+
+ if(!returnCurv) {
+ return Su.createNeutralParameters(p[0],p[1],p[2],p[3],p[4],e);
+ }
+ else {
+ Amg::Vector3D gp(P[0],P[1],P[2]);
+ return new Trk::NeutralCurvilinearParameters(gp,p[2],p[3],p[4],e);
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for charged track parameters propagation with or without jacobian
+/////////////////////////////////////////////////////////////////////////////////
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::propagateRungeKutta
+(bool useJac,
+ const Trk::TrackParameters & Tp ,
+ const Trk::Surface & Su ,
+ Trk::PropDirection D ,
+ Trk::BoundaryCheck B ,
+ const MagneticFieldProperties& M ,
+ double * Jac ,
+ bool returnCurv) const
+{
+ const Trk::Surface* su = &Su;
+
+ if(!&Tp || !su) return 0;
+
+ m_direction = D ;
+
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ (useJac && m_usegradient)? m_needgradient = true : m_needgradient = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ if(su == &Tp.associatedSurface()) return buildTrackParametersWithoutPropagation(Tp,Jac);
+
+ Trk::RungeKuttaUtils utils;
+
+ double P[64],Step = 0.; if(!utils.transformLocalToGlobal(useJac,Tp,P)) return 0;
+
+ const Amg::Transform3D& T = Su.transform();
+ int ty = Su.type();
+
+ if (ty == Trk::Surface::Plane ) {
+
+ double s[4], d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return 0;
+
+ }
+ else if (ty == Trk::Surface::Line ) {
+
+ double s[6] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Disc ) {
+
+ double s[4], d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cylinder ) {
+
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(su);
+
+ double r0[3] = {P[0],P[1],P[2]};
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),m_direction,0.};
+ if(!propagateWithJacobian(useJac,2,s,P,Step)) return 0;
+
+ // For cylinder we do test for next cross point
+ //
+ if(cyl->bounds().halfPhiSector() < 3.1 && newCrossPoint(*cyl,r0,P)) {
+ s[8] = 0.; if(!propagateWithJacobian(useJac,2,s,P,Step)) return 0;
+ }
+ }
+ else if (ty == Trk::Surface::Perigee ) {
+
+ double s[6] = {T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cone ) {
+
+ double k = static_cast<const Trk::ConeSurface*>(su)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,m_direction,0.};
+ if(!propagateWithJacobian(useJac,3,s,P,Step)) return 0;
+ }
+ else return 0;
+
+ if(m_direction && (m_direction*Step)<0.) {return 0;}
+
+ // Common transformation for all surfaces (angles and momentum)
+ //
+ if(useJac) {
+ double p=1./P[6]; P[35]*=p; P[36]*=p; P[37]*=p; P[38]*=p; P[39]*=p; P[40]*=p;
+ }
+
+ bool uJ = useJac; if(returnCurv) uJ = false;
+ double p[5]; utils.transformGlobalToLocal(su,uJ,P,p,Jac);
+
+ if(B) {Amg::Vector2D L(p[0],p[1]); if(!Su.insideBounds(L,0.)) return 0;}
+
+ // Transformation to curvilinear presentation
+ //
+ if(returnCurv) utils.transformGlobalToCurvilinear(useJac,P,p,Jac);
+
+ if(!useJac || !Tp.covariance()) {
+
+ if(!returnCurv) {
+ return Su.createTrackParameters(p[0],p[1],p[2],p[3],p[4],0);
+ }
+ else {
+ Amg::Vector3D gp(P[0],P[1],P[2]);
+ return new Trk::CurvilinearParameters(gp,p[2],p[3],p[4]);
+ }
+ }
+
+ AmgSymMatrix(5)* e = utils.newCovarianceMatrix(Jac,*Tp.covariance());
+ AmgSymMatrix(5)& cv = *e;
+
+ if(cv(0,0)<=0. || cv(1,1)<=0. || cv(2,2)<=0. || cv(3,3)<=0. || cv(4,4)<=0.) {
+ delete e; return 0;
+ }
+
+ if(!returnCurv) {
+ return Su.createTrackParameters(p[0],p[1],p[2],p[3],p[4],e);
+ }
+ else {
+ Amg::Vector3D gp(P[0],P[1],P[2]);
+ return new Trk::CurvilinearParameters(gp,p[2],p[3],p[4],e);
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Global positions calculation inside CylinderBounds
+// where mS - max step allowed if mS > 0 propagate along momentum
+// mS < 0 propogate opposite momentum
+/////////////////////////////////////////////////////////////////////////////////
+
+void Trk::RungeKuttaPropagator::globalPositions
+(std::list<Amg::Vector3D> & GP,
+ const TrackParameters & Tp,
+ const MagneticFieldProperties & M ,
+ const CylinderBounds & CB,
+ double mS,
+ ParticleHypothesis ,
+ const TrackingVolume* ) const
+{
+ Trk::RungeKuttaUtils utils;
+ double P[45]; if(!utils.transformLocalToGlobal(false,Tp,P)) return;
+
+ m_direction = fabs(mS);
+ if(mS > 0.) globalOneSidePositions(GP,P,M,CB, mS);
+ else globalTwoSidePositions(GP,P,M,CB,-mS);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Global position together with direction of the trajectory on the surface
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::IntersectionSolution* Trk::RungeKuttaPropagator::intersect
+( const Trk::TrackParameters & Tp,
+ const Trk::Surface & Su,
+ const MagneticFieldProperties& M ,
+ ParticleHypothesis ,
+ const TrackingVolume* ) const
+{
+ bool nJ = false;
+ const Trk::Surface* su = &Su;
+ m_direction = 0. ;
+
+ m_needgradient = false;
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ Trk::RungeKuttaUtils utils;
+
+ double P[64]; if(!utils.transformLocalToGlobal(false,Tp,P)) return 0;
+ double Step = 0.;
+
+ const Amg::Transform3D& T = Su.transform();
+ int ty = Su.type();
+
+ if (ty == Trk::Surface::Plane ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(nJ,1,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Line ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ if(!propagateWithJacobian(nJ,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Disc ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(nJ,1,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cylinder ) {
+
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(su);
+
+ double r0[3] = {P[0],P[1],P[2]};
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),m_direction,0.};
+
+ if(!propagateWithJacobian(nJ,2,s,P,Step)) return 0;
+
+ // For cylinder we do test for next cross point
+ //
+ if(cyl->bounds().halfPhiSector() < 3.1 && newCrossPoint(*cyl,r0,P)) {
+ s[8] = 0.; if(!propagateWithJacobian(nJ,2,s,P,Step)) return 0;
+ }
+ }
+ else if (ty == Trk::Surface::Perigee ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ if(!propagateWithJacobian(nJ,0,s,P,Step)) return 0;
+ }
+ else if (ty == Trk::Surface::Cone ) {
+
+ double k = static_cast<const Trk::ConeSurface*>(su)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,m_direction,0.};
+ if(!propagateWithJacobian(nJ,3,s,P,Step)) return 0;
+ }
+ else return 0;
+
+ Amg::Vector3D Glo(P[0],P[1],P[2]);
+ Amg::Vector3D Dir(P[3],P[4],P[5]);
+ Trk::IntersectionSolution* Int = new Trk::IntersectionSolution();
+ Int->push_back(new Trk::TrackSurfaceIntersection(Glo,Dir,Step));
+ return Int;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Runge Kutta main program for propagation with or without Jacobian
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::propagateWithJacobian
+(bool Jac ,
+ int kind ,
+ double * Su,
+ double * P ,
+ double & W ) const
+{
+ const double Smax = 1000. ; // max. step allowed
+ double Wmax = 100000. ; // Max way allowed
+ double Wwrong = 500. ; // Max way with wrong direction
+ double* R = &P[ 0] ; // Start coordinates
+ double* A = &P[ 3] ; // Start directions
+ double* SA = &P[42] ; SA[0]=SA[1]=SA[2]=0.;
+
+ if(m_mcondition && fabs(P[6]) > .1) return false;
+
+ // Step estimation until surface
+ //
+ Trk::RungeKuttaUtils utils;
+ bool Q; double S, Step=utils.stepEstimator(kind,Su,P,Q); if(!Q) return false;
+
+ bool dir = true;
+ if(m_mcondition && m_direction && m_direction*Step < 0.) {
+ Step = -Step; dir = false;
+ }
+
+ Step>Smax ? S=Smax : Step<-Smax ? S=-Smax : S=Step;
+ double So = fabs(S); int iS = 0;
+
+ bool InS = false;
+
+ // Rkuta extrapolation
+ //
+ int niter = 0;
+ while(fabs(Step) > m_straightStep) {
+
+ if(++niter > 10000) return false;
+
+ if(m_mcondition) {
+
+ if(!m_needgradient) W+=(S=rungeKuttaStep (Jac,S,P,InS));
+ else W+=(S=rungeKuttaStepWithGradient( S,P,InS));
+ }
+ else {
+ W+=(S=straightLineStep(Jac,S,P));
+ }
+
+ Step = stepEstimatorWithCurvature(kind,Su,P,Q); if(!Q) return false;
+
+ if(!dir) {
+ if(m_direction && m_direction*Step < 0.) Step = -Step;
+ else dir = true;
+ }
+
+ if(S*Step<0.) {S = -S; ++iS;}
+
+ double aS = fabs(S );
+ double aStep = fabs(Step);
+ if ( aS > aStep ) S = Step;
+ else if(!iS && InS && aS*2. < aStep) S*=2. ;
+ if(!dir && fabs(W) > Wwrong ) return false;
+ if(iS > 10 || (iS>3 && fabs(S)>=So) || fabs(W)>Wmax ) {if(!kind) break; return false;}
+ So=fabs(S);
+ }
+
+ // Output track parameteres
+ //
+ W+=Step;
+
+ if(fabs(Step) < .001) return true;
+
+ A [0]+=(SA[0]*Step);
+ A [1]+=(SA[1]*Step);
+ A [2]+=(SA[2]*Step);
+ double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
+ R[0]+=Step*(A[0]-.5*Step*SA[0]); A[0]*=CBA;
+ R[1]+=Step*(A[1]-.5*Step*SA[1]); A[1]*=CBA;
+ R[2]+=Step*(A[2]-.5*Step*SA[2]); A[2]*=CBA;
+ return true;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Runge Kutta trajectory model (units->mm,MeV,kGauss)
+// Uses Nystroem algorithm (See Handbook Net. Bur. ofStandards, procedure 25.5.20)
+/////////////////////////////////////////////////////////////////////////////////
+
+double Trk::RungeKuttaPropagator::rungeKuttaStep
+ (bool Jac,
+ double S ,
+ double * P ,
+ bool & InS) const
+{
+ const double C33 = 1./3. ;
+ double* R = &P[ 0]; // Coordinates
+ double* A = &P[ 3]; // Directions
+ double* sA = &P[42];
+ double Pi = 149.89626*P[6]; // Invert mometum/2.
+ double dltm = m_dlt*.03 ;
+
+ double f0[3],f[3]; getField(R,f0);
+
+ bool Helix = false; if(fabs(S) < m_helixStep) Helix = true;
+ while(S != 0.) {
+
+ double S3=C33*S, S4=.25*S, PS2=Pi*S;
+
+ // First point
+ //
+ double H0[3] = {f0[0]*PS2, f0[1]*PS2, f0[2]*PS2};
+ double A0 = A[1]*H0[2]-A[2]*H0[1] ;
+ double B0 = A[2]*H0[0]-A[0]*H0[2] ;
+ double C0 = A[0]*H0[1]-A[1]*H0[0] ;
+ double A2 = A[0]+A0 ;
+ double B2 = A[1]+B0 ;
+ double C2 = A[2]+C0 ;
+ double A1 = A2+A[0] ;
+ double B1 = B2+A[1] ;
+ double C1 = C2+A[2] ;
+
+ // Second point
+ //
+ if(!Helix) {
+ double gP[3]={R[0]+A1*S4, R[1]+B1*S4, R[2]+C1*S4};
+ getField(gP,f);
+ }
+ else {f[0]=f0[0]; f[1]=f0[1]; f[2]=f0[2];}
+
+ double H1[3] = {f[0]*PS2,f[1]*PS2,f[2]*PS2};
+ double A3 = B2*H1[2]-C2*H1[1]+A[0] ;
+ double B3 = C2*H1[0]-A2*H1[2]+A[1] ;
+ double C3 = A2*H1[1]-B2*H1[0]+A[2] ;
+ double A4 = B3*H1[2]-C3*H1[1]+A[0] ;
+ double B4 = C3*H1[0]-A3*H1[2]+A[1] ;
+ double C4 = A3*H1[1]-B3*H1[0]+A[2] ;
+ double A5 = A4-A[0]+A4 ;
+ double B5 = B4-A[1]+B4 ;
+ double C5 = C4-A[2]+C4 ;
+
+ // Last point
+ //
+ if(!Helix) {
+ double gP[3]={R[0]+S*A4, R[1]+S*B4, R[2]+S*C4};
+ getField(gP,f);
+ }
+ else {f[0]=f0[0]; f[1]=f0[1]; f[2]=f0[2];}
+
+ double H2[3] = {f[0]*PS2,f[1]*PS2,f[2]*PS2};
+ double A6 = B5*H2[2]-C5*H2[1] ;
+ double B6 = C5*H2[0]-A5*H2[2] ;
+ double C6 = A5*H2[1]-B5*H2[0] ;
+
+
+ // Test approximation quality on give step and possible step reduction
+ //
+ //double dE[4] = {(A1+A6)-(A3+A4),(B1+B6)-(B3+B4),(C1+C6)-(C3+C4),S};
+ //double cS = stepReduction(dE);
+ //if (cSn < 1.) {S*=cS; continue;}
+ //cSn == 1. InS = false : InS = true;
+ //
+
+ // Test approximation quality on give step and possible step reduction
+ //
+ double EST = fabs((A1+A6)-(A3+A4))+fabs((B1+B6)-(B3+B4))+fabs((C1+C6)-(C3+C4));
+ if(EST>m_dlt) {S*=.5; dltm = 0.; continue;} EST<dltm ? InS = true : InS = false;
+
+ // Parameters calculation
+ //
+ double A00 = A[0], A11=A[1], A22=A[2];
+ R[0]+=(A2+A3+A4)*S3; A[0] = ((A0+2.*A3)+(A5+A6))*C33;
+ R[1]+=(B2+B3+B4)*S3; A[1] = ((B0+2.*B3)+(B5+B6))*C33;
+ R[2]+=(C2+C3+C4)*S3; A[2] = ((C0+2.*C3)+(C5+C6))*C33;
+ double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
+ A[0]*=CBA; A[1]*=CBA; A[2]*=CBA;
+
+ double Sl = 2./S;
+ sA[0] = A6*Sl;
+ sA[1] = B6*Sl;
+ sA[2] = C6*Sl;
+
+ if(!Jac) return S;
+
+ // Jacobian calculation
+ //
+ for(int i=21; i<35; i+=7) {
+
+ double* dR = &P[i] ;
+ double* dA = &P[i+3] ;
+ double dA0 = H0[ 2]*dA[1]-H0[ 1]*dA[2] ;
+ double dB0 = H0[ 0]*dA[2]-H0[ 2]*dA[0] ;
+ double dC0 = H0[ 1]*dA[0]-H0[ 0]*dA[1] ;
+ double dA2 = dA0+dA[0] ;
+ double dB2 = dB0+dA[1] ;
+ double dC2 = dC0+dA[2] ;
+ double dA3 = dA[0]+dB2*H1[2]-dC2*H1[1] ;
+ double dB3 = dA[1]+dC2*H1[0]-dA2*H1[2] ;
+ double dC3 = dA[2]+dA2*H1[1]-dB2*H1[0] ;
+ double dA4 = dA[0]+dB3*H1[2]-dC3*H1[1] ;
+ double dB4 = dA[1]+dC3*H1[0]-dA3*H1[2] ;
+ double dC4 = dA[2]+dA3*H1[1]-dB3*H1[0] ;
+ double dA5 = dA4+dA4-dA[0] ;
+ double dB5 = dB4+dB4-dA[1] ;
+ double dC5 = dC4+dC4-dA[2] ;
+ double dA6 = dB5*H2[2]-dC5*H2[1] ;
+ double dB6 = dC5*H2[0]-dA5*H2[2] ;
+ double dC6 = dA5*H2[1]-dB5*H2[0] ;
+ dR[0]+=(dA2+dA3+dA4)*S3; dA[0]=((dA0+2.*dA3)+(dA5+dA6))*C33 ;
+ dR[1]+=(dB2+dB3+dB4)*S3; dA[1]=((dB0+2.*dB3)+(dB5+dB6))*C33 ;
+ dR[2]+=(dC2+dC3+dC4)*S3; dA[2]=((dC0+2.*dC3)+(dC5+dC6))*C33 ;
+ }
+
+ double* dR = &P[35] ;
+ double* dA = &P[38] ;
+ double dA0 = H0[ 2]*dA[1]-H0[ 1]*dA[2]+A0 ;
+ double dB0 = H0[ 0]*dA[2]-H0[ 2]*dA[0]+B0 ;
+ double dC0 = H0[ 1]*dA[0]-H0[ 0]*dA[1]+C0 ;
+ double dA2 = dA0+dA[0] ;
+ double dB2 = dB0+dA[1] ;
+ double dC2 = dC0+dA[2] ;
+ double dA3 = (dA[0]+dB2*H1[2]-dC2*H1[1])+(A3-A00) ;
+ double dB3 = (dA[1]+dC2*H1[0]-dA2*H1[2])+(B3-A11) ;
+ double dC3 = (dA[2]+dA2*H1[1]-dB2*H1[0])+(C3-A22) ;
+ double dA4 = (dA[0]+dB3*H1[2]-dC3*H1[1])+(A4-A00) ;
+ double dB4 = (dA[1]+dC3*H1[0]-dA3*H1[2])+(B4-A11) ;
+ double dC4 = (dA[2]+dA3*H1[1]-dB3*H1[0])+(C4-A22) ;
+ double dA5 = dA4+dA4-dA[0] ;
+ double dB5 = dB4+dB4-dA[1] ;
+ double dC5 = dC4+dC4-dA[2] ;
+ double dA6 = dB5*H2[2]-dC5*H2[1]+A6 ;
+ double dB6 = dC5*H2[0]-dA5*H2[2]+B6 ;
+ double dC6 = dA5*H2[1]-dB5*H2[0]+C6 ;
+ dR[0]+=(dA2+dA3+dA4)*S3; dA[0]=((dA0+2.*dA3)+(dA5+dA6))*C33 ;
+ dR[1]+=(dB2+dB3+dB4)*S3; dA[1]=((dB0+2.*dB3)+(dB5+dB6))*C33 ;
+ dR[2]+=(dC2+dC3+dC4)*S3; dA[2]=((dC0+2.*dC3)+(dC5+dC6))*C33 ;
+ return S;
+ }
+ return S;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Runge Kutta trajectory model (units->mm,MeV,kGauss)
+// Uses Nystroem algorithm (See Handbook Net. Bur. ofStandards, procedure 25.5.20)
+// Where magnetic field information iS
+// f[ 0],f[ 1],f[ 2] - Hx , Hy and Hz of the magnetic field
+// f[ 3],f[ 4],f[ 5] - dHx/dx, dHx/dy and dHx/dz
+// f[ 6],f[ 7],f[ 8] - dHy/dx, dHy/dy and dHy/dz
+// f[ 9],f[10],f[11] - dHz/dx, dHz/dy and dHz/dz
+//
+/////////////////////////////////////////////////////////////////////////////////
+
+double Trk::RungeKuttaPropagator::rungeKuttaStepWithGradient
+(double S ,
+ double * P ,
+ bool & InS) const
+{
+ const double C33 = 1./3. ;
+ double* R = &P[ 0]; // Coordinates
+ double* A = &P[ 3]; // Directions
+ double* sA = &P[42];
+ double Pi = 149.89626*P[6]; // Invert mometum/2.
+ double dltm = m_dlt*.03 ;
+
+ double f0[3],f1[3],f2[3],g0[9],g1[9],g2[9],H0[12],H1[12],H2[12];
+ getFieldGradient(R,f0,g0);
+
+ while(S != 0.) {
+
+
+ double S3=C33*S, S4=.25*S, PS2=Pi*S;
+
+ // First point
+ //
+ H0[0] = f0[0]*PS2; H0[1] = f0[1]*PS2; H0[2] = f0[2]*PS2;
+ double A0 = A[1]*H0[2]-A[2]*H0[1] ;
+ double B0 = A[2]*H0[0]-A[0]*H0[2] ;
+ double C0 = A[0]*H0[1]-A[1]*H0[0] ;
+ double A2 = A[0]+A0 ;
+ double B2 = A[1]+B0 ;
+ double C2 = A[2]+C0 ;
+ double A1 = A2+A[0] ;
+ double B1 = B2+A[1] ;
+ double C1 = C2+A[2] ;
+
+ // Second point
+ //
+ double gP1[3]={R[0]+A1*S4, R[1]+B1*S4, R[2]+C1*S4};
+ getFieldGradient(gP1,f1,g1);
+
+ H1[0] = f1[0]*PS2; H1[1] = f1[1]*PS2; H1[2] = f1[2]*PS2;
+ double A3 = B2*H1[2]-C2*H1[1]+A[0] ;
+ double B3 = C2*H1[0]-A2*H1[2]+A[1] ;
+ double C3 = A2*H1[1]-B2*H1[0]+A[2] ;
+ double A4 = B3*H1[2]-C3*H1[1]+A[0] ;
+ double B4 = C3*H1[0]-A3*H1[2]+A[1] ;
+ double C4 = A3*H1[1]-B3*H1[0]+A[2] ;
+ double A5 = A4-A[0]+A4 ;
+ double B5 = B4-A[1]+B4 ;
+ double C5 = C4-A[2]+C4 ;
+
+ // Last point
+ //
+ double gP2[3]={R[0]+S*A4, R[1]+S*B4, R[2]+S*C4};
+ getFieldGradient(gP2,f2,g2);
+
+ H2[0] = f2[0]*PS2; H2[1] = f2[1]*PS2; H2[2] = f2[2]*PS2;
+ double A6 = B5*H2[2]-C5*H2[1] ;
+ double B6 = C5*H2[0]-A5*H2[2] ;
+ double C6 = A5*H2[1]-B5*H2[0] ;
+
+
+ // Test approximation quality on give step and possible step reduction
+ /*
+ double dE[4] = {(A1+A6)-(A3+A4),(B1+B6)-(B3+B4),(C1+C6)-(C3+C4),S};
+ double cS = stepReduction(dE);
+ if (cSn < 1.) {S*=cS; continue;}
+ cSn == 1. InS = false : InS = true;
+ */
+
+ // Test approximation quality on give step and possible step reduction
+ //
+ double EST = fabs((A1+A6)-(A3+A4))+fabs((B1+B6)-(B3+B4))+fabs((C1+C6)-(C3+C4));
+ if(EST>m_dlt) {S*=.5; dltm = 0.; continue;} EST<dltm ? InS = true : InS = false;
+
+ // Parameters calculation
+ //
+ double A00 = A[0], A11=A[1], A22=A[2];
+ R[0]+=(A2+A3+A4)*S3; A[0] = ((A0+2.*A3)+(A5+A6))*C33;
+ R[1]+=(B2+B3+B4)*S3; A[1] = ((B0+2.*B3)+(B5+B6))*C33;
+ R[2]+=(C2+C3+C4)*S3; A[2] = ((C0+2.*C3)+(C5+C6))*C33;
+ double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
+ A[0]*=CBA; A[1]*=CBA; A[2]*=CBA;
+
+ double Sl = 2./S;
+ sA[0] = A6*Sl;
+ sA[1] = B6*Sl;
+ sA[2] = C6*Sl;
+
+ // Jacobian calculation
+ //
+ for(int i=3; i!=12; ++i) {H0[i]=g0[i-3]*PS2; H1[i]=g1[i-3]*PS2; H2[i]=g2[i-3]*PS2;}
+
+ for(int i=7; i<35; i+=7) {
+
+ double* dR = &P[i] ;
+ double* dA = &P[i+3] ;
+ double dH0 = H0[ 3]*dR[0]+H0[ 4]*dR[1]+H0[ 5]*dR[2] ; // dHx/dp
+ double dH1 = H0[ 6]*dR[0]+H0[ 7]*dR[1]+H0[ 8]*dR[2] ; // dHy/dp
+ double dH2 = H0[ 9]*dR[0]+H0[10]*dR[1]+H0[11]*dR[2] ; // dHz/dp
+ double dA0 =(H0[ 2]*dA[1]-H0[ 1]*dA[2])+(A[1]*dH2-A[2]*dH1) ; // dA0/dp
+ double dB0 =(H0[ 0]*dA[2]-H0[ 2]*dA[0])+(A[2]*dH0-A[0]*dH2) ; // dB0/dp
+ double dC0 =(H0[ 1]*dA[0]-H0[ 0]*dA[1])+(A[0]*dH1-A[1]*dH0) ; // dC0/dp
+ double dA2 = dA0+dA[0], dX = dR[0]+(dA2+dA[0])*S4 ; // dX /dp
+ double dB2 = dB0+dA[1], dY = dR[1]+(dB2+dA[1])*S4 ; // dY /dp
+ double dC2 = dC0+dA[2], dZ = dR[2]+(dC2+dA[2])*S4 ; // dZ /dp
+ dH0 = H1[ 3]*dX +H1[ 4]*dY +H1[ 5]*dZ ; // dHx/dp
+ dH1 = H1[ 6]*dX +H1[ 7]*dY +H1[ 8]*dZ ; // dHy/dp
+ dH2 = H1[ 9]*dX +H1[10]*dY +H1[11]*dZ ; // dHz/dp
+ double dA3 =(dA[0]+dB2*H1[2]-dC2*H1[1])+(B2*dH2-C2*dH1) ; // dA3/dp
+ double dB3 =(dA[1]+dC2*H1[0]-dA2*H1[2])+(C2*dH0-A2*dH2) ; // dB3/dp
+ double dC3 =(dA[2]+dA2*H1[1]-dB2*H1[0])+(A2*dH1-B2*dH0) ; // dC3/dp
+ double dA4 =(dA[0]+dB3*H1[2]-dC3*H1[1])+(B3*dH2-C3*dH1) ; // dA4/dp
+ double dB4 =(dA[1]+dC3*H1[0]-dA3*H1[2])+(C3*dH0-A3*dH2) ; // dB4/dp
+ double dC4 =(dA[2]+dA3*H1[1]-dB3*H1[0])+(A3*dH1-B3*dH0) ; // dC4/dp
+ double dA5 = dA4+dA4-dA[0]; dX = dR[0]+dA4*S ; // dX /dp
+ double dB5 = dB4+dB4-dA[1]; dY = dR[1]+dB4*S ; // dY /dp
+ double dC5 = dC4+dC4-dA[2]; dZ = dR[2]+dC4*S ; // dZ /dp
+ dH0 = H2[ 3]*dX +H2[ 4]*dY +H2[ 5]*dZ ; // dHx/dp
+ dH1 = H2[ 6]*dX +H2[ 7]*dY +H2[ 8]*dZ ; // dHy/dp
+ dH2 = H2[ 9]*dX +H2[10]*dY +H2[11]*dZ ; // dHz/dp
+ double dA6 =(dB5*H2[2]-dC5*H2[1])+(B5*dH2-C5*dH1) ; // dA6/dp
+ double dB6 =(dC5*H2[0]-dA5*H2[2])+(C5*dH0-A5*dH2) ; // dB6/dp
+ double dC6 =(dA5*H2[1]-dB5*H2[0])+(A5*dH1-B5*dH0) ; // dC6/dp
+ dR[0]+=(dA2+dA3+dA4)*S3; dA[0]=((dA0+2.*dA3)+(dA5+dA6))*C33 ;
+ dR[1]+=(dB2+dB3+dB4)*S3; dA[1]=((dB0+2.*dB3)+(dB5+dB6))*C33 ;
+ dR[2]+=(dC2+dC3+dC4)*S3; dA[2]=((dC0+2.*dC3)+(dC5+dC6))*C33 ;
+ }
+
+ double* dR = &P[35] ;
+ double* dA = &P[38] ;
+
+ double dH0 = H0[ 3]*dR[0]+H0[ 4]*dR[1]+H0[ 5]*dR[2] ; // dHx/dp
+ double dH1 = H0[ 6]*dR[0]+H0[ 7]*dR[1]+H0[ 8]*dR[2] ; // dHy/dp
+ double dH2 = H0[ 9]*dR[0]+H0[10]*dR[1]+H0[11]*dR[2] ; // dHz/dp
+ double dA0 =(H0[ 2]*dA[1]-H0[ 1]*dA[2])+(A[1]*dH2-A[2]*dH1+A0) ; // dA0/dp
+ double dB0 =(H0[ 0]*dA[2]-H0[ 2]*dA[0])+(A[2]*dH0-A[0]*dH2+B0) ; // dB0/dp
+ double dC0 =(H0[ 1]*dA[0]-H0[ 0]*dA[1])+(A[0]*dH1-A[1]*dH0+C0) ; // dC0/dp
+ double dA2 = dA0+dA[0], dX = dR[0]+(dA2+dA[0])*S4 ; // dX /dp
+ double dB2 = dB0+dA[1], dY = dR[1]+(dB2+dA[1])*S4 ; // dY /dp
+ double dC2 = dC0+dA[2], dZ = dR[2]+(dC2+dA[2])*S4 ; // dZ /dp
+ dH0 = H1[ 3]*dX +H1[ 4]*dY +H1[ 5]*dZ ; // dHx/dp
+ dH1 = H1[ 6]*dX +H1[ 7]*dY +H1[ 8]*dZ ; // dHy/dp
+ dH2 = H1[ 9]*dX +H1[10]*dY +H1[11]*dZ ; // dHz/dp
+ double dA3 =(dA[0]+dB2*H1[2]-dC2*H1[1])+((B2*dH2-C2*dH1)+(A3-A00)) ; // dA3/dp
+ double dB3 =(dA[1]+dC2*H1[0]-dA2*H1[2])+((C2*dH0-A2*dH2)+(B3-A11)) ; // dB3/dp
+ double dC3 =(dA[2]+dA2*H1[1]-dB2*H1[0])+((A2*dH1-B2*dH0)+(C3-A22)) ; // dC3/dp
+ double dA4 =(dA[0]+dB3*H1[2]-dC3*H1[1])+((B3*dH2-C3*dH1)+(A4-A00)) ; // dA4/dp
+ double dB4 =(dA[1]+dC3*H1[0]-dA3*H1[2])+((C3*dH0-A3*dH2)+(B4-A11)) ; // dB4/dp
+ double dC4 =(dA[2]+dA3*H1[1]-dB3*H1[0])+((A3*dH1-B3*dH0)+(C4-A22)) ; // dC4/dp
+ double dA5 = dA4+dA4-dA[0]; dX = dR[0]+dA4*S ; // dX /dp
+ double dB5 = dB4+dB4-dA[1]; dY = dR[1]+dB4*S ; // dY /dp
+ double dC5 = dC4+dC4-dA[2]; dZ = dR[2]+dC4*S ; // dZ /dp
+ dH0 = H2[ 3]*dX +H2[ 4]*dY +H2[ 5]*dZ ; // dHx/dp
+ dH1 = H2[ 6]*dX +H2[ 7]*dY +H2[ 8]*dZ ; // dHy/dp
+ dH2 = H2[ 9]*dX +H2[10]*dY +H2[11]*dZ ; // dHz/dp
+ double dA6 =(dB5*H2[2]-dC5*H2[1])+(B5*dH2-C5*dH1+A6) ; // dA6/dp
+ double dB6 =(dC5*H2[0]-dA5*H2[2])+(C5*dH0-A5*dH2+B6) ; // dB6/dp
+ double dC6 =(dA5*H2[1]-dB5*H2[0])+(A5*dH1-B5*dH0+C6) ; // dC6/dp
+ dR[0]+=(dA2+dA3+dA4)*S3; dA[0]=((dA0+2.*dA3)+(dA5+dA6))*C33 ;
+ dR[1]+=(dB2+dB3+dB4)*S3; dA[1]=((dB0+2.*dB3)+(dB5+dB6))*C33 ;
+ dR[2]+=(dC2+dC3+dC4)*S3; dA[2]=((dC0+2.*dC3)+(dC5+dC6))*C33 ;
+ return S;
+ }
+ return S;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Straight line trajectory model
+/////////////////////////////////////////////////////////////////////////////////
+
+double Trk::RungeKuttaPropagator::straightLineStep
+(bool Jac,
+ double S ,
+ double* P ) const
+{
+ double* R = &P[ 0]; // Start coordinates
+ double* A = &P[ 3]; // Start directions
+ double* sA = &P[42];
+
+ // Track parameters in last point
+ //
+ R[0]+=(A[0]*S); R[1]+=(A[1]*S); R[2]+=(A[2]*S); if(!Jac) return S;
+
+ // Derivatives of track parameters in last point
+ //
+ for(int i=7; i<42; i+=7) {
+
+ double* dR = &P[i];
+ double* dA = &P[i+3];
+ dR[0]+=(dA[0]*S); dR[1]+=(dA[1]*S); dR[2]+=(dA[2]*S);
+ }
+ sA[0]=sA[1]=sA[2]=0.; return S;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for simple track parameters and covariance matrix propagation
+// Ta->Su = Tb
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::propagate
+(Trk::PatternTrackParameters & Ta,
+ const Trk::Surface & Su,
+ Trk::PatternTrackParameters & Tb,
+ Trk::PropDirection D ,
+ const MagneticFieldProperties& M ,
+ ParticleHypothesis ) const
+{
+ double S;
+ return propagateRungeKutta(true,Ta,Su,Tb,D,M,S);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for simple track parameters and covariance matrix propagation
+// Ta->Su = Tb with step to surface calculation
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::propagate
+(Trk::PatternTrackParameters & Ta,
+ const Trk::Surface & Su,
+ Trk::PatternTrackParameters & Tb,
+ Trk::PropDirection D ,
+ const MagneticFieldProperties& M ,
+ double & S ,
+ ParticleHypothesis ) const
+{
+ return propagateRungeKutta(true,Ta,Su,Tb,D,M,S);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for simple track parameters propagation without covariance matrix
+// Ta->Su = Tb
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::propagateParameters
+(Trk::PatternTrackParameters & Ta,
+ const Trk::Surface & Su,
+ Trk::PatternTrackParameters & Tb,
+ Trk::PropDirection D ,
+ const MagneticFieldProperties& M ,
+ ParticleHypothesis ) const
+{
+ double S;
+ return propagateRungeKutta(false,Ta,Su,Tb,D,M,S);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for simple track parameters propagation without covariance matrix
+// Ta->Su = Tb with step calculation
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::propagateParameters
+(Trk::PatternTrackParameters & Ta,
+ const Trk::Surface & Su,
+ Trk::PatternTrackParameters & Tb,
+ Trk::PropDirection D ,
+ const MagneticFieldProperties& M ,
+ double & S ,
+ ParticleHypothesis ) const
+{
+ return propagateRungeKutta(false,Ta,Su,Tb,D,M,S);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Global positions calculation inside CylinderBounds
+// where mS - max step allowed if mS > 0 propagate along momentum
+// mS < 0 propogate opposite momentum
+/////////////////////////////////////////////////////////////////////////////////
+
+void Trk::RungeKuttaPropagator::globalPositions
+(std::list<Amg::Vector3D> & GP,
+ const Trk::PatternTrackParameters& Tp,
+ const MagneticFieldProperties & M,
+ const CylinderBounds & CB,
+ double mS,
+ ParticleHypothesis ) const
+{
+ Trk::RungeKuttaUtils utils;
+ double P[45]; if(!utils.transformLocalToGlobal(false,Tp,P)) return;
+
+ m_direction = fabs(mS);
+ if(mS > 0.) globalOneSidePositions(GP,P,M,CB, mS);
+ else globalTwoSidePositions(GP,P,M,CB,-mS);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// GlobalPostions and steps for set surfaces
+/////////////////////////////////////////////////////////////////////////////////
+
+void Trk::RungeKuttaPropagator::globalPositions
+(const PatternTrackParameters & Tp,
+ std::list<const Trk::Surface*> & SU,
+ std::list< std::pair<Amg::Vector3D,double> > & GP,
+ const Trk::MagneticFieldProperties & M ,
+ ParticleHypothesis ) const
+{
+ m_direction = 0. ;
+ m_mcondition = false;
+
+ m_needgradient = false;
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ Trk::RungeKuttaUtils utils;
+
+ double Step = 0.,P[64]; if(!utils.transformLocalToGlobal(false,Tp,P)) return;
+
+
+ std::list<const Trk::Surface*>::iterator su = SU.begin(), sue = SU.end();
+
+ // Loop trough all input surfaces
+ //
+ for(; su!=sue; ++su) {
+
+ const Amg::Transform3D& T = (*su)->transform();
+ int ty = (*su)->type();
+
+ if( ty == Trk::Surface::Plane ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(false,1,s,P,Step)) return;
+ }
+ else if(ty == Trk::Surface::Line ) {
+
+ double s[6]={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ if(!propagateWithJacobian(false,0,s,P,Step)) return;
+ }
+ else if (ty == Trk::Surface::Disc ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(false,1,s,P,Step)) return;
+ }
+ else if (ty == Trk::Surface::Cylinder ) {
+
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(*su);
+ double r0[3] = {P[0],P[1],P[2]};
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),m_direction,0.};
+
+ if(!propagateWithJacobian(false,2,s,P,Step)) return;
+
+ // For cylinder we do test for next cross point
+ //
+ if(cyl->bounds().halfPhiSector() < 3.1 && newCrossPoint(*cyl,r0,P)) {
+ s[8] = 0.; if(!propagateWithJacobian(false,2,s,P,Step)) return;
+ }
+ }
+ else if (ty == Trk::Surface::Perigee ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ if(!propagateWithJacobian(false,0,s,P,Step)) return ;
+ }
+ else if (ty == Trk::Surface::Cone ) {
+
+ double k = static_cast<const Trk::ConeSurface*>(*su)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,m_direction,0.};
+ if(!propagateWithJacobian(false,3,s,P,Step)) return;
+ }
+ else return;
+
+ Amg::Vector3D gp(P[0],P[1],P[2]); GP.push_back(std::make_pair(gp,Step));
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Main function for simple track propagation with or without jacobian
+// Ta->Su = Tb for pattern track parameters
+/////////////////////////////////////////////////////////////////////////////////
+bool Trk::RungeKuttaPropagator::propagateRungeKutta
+(bool useJac,
+ Trk::PatternTrackParameters & Ta ,
+ const Trk::Surface & Su ,
+ Trk::PatternTrackParameters & Tb ,
+ Trk::PropDirection D ,
+ const MagneticFieldProperties& M ,
+ double & Step ) const
+{
+ const Trk::Surface* su = &Su;
+ if(su == Ta.associatedSurface()) {Tb = Ta; return true;}
+
+ m_direction = D ;
+
+ if(useJac && !Ta.iscovariance()) useJac = false;
+
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ (useJac && m_usegradient) ? m_needgradient = true : m_needgradient = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ Trk::RungeKuttaUtils utils;
+
+ double P[45]; if(!utils.transformLocalToGlobal(useJac,Ta,P)) return false; Step = 0.;
+
+ const Amg::Transform3D& T = Su.transform();
+ int ty = Su.type();
+
+ if (ty == Trk::Surface::Plane ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return false;
+ }
+ else if (ty == Trk::Surface::Line ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2)};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return false;
+ }
+ else if (ty == Trk::Surface::Disc ) {
+
+ double s[4],d = T(0,3)*T(0,2)+T(1,3)*T(1,2)+T(2,3)*T(2,2);
+
+ if(d>=0.) {s[0]= T(0,2); s[1]= T(1,2); s[2]= T(2,2); s[3]= d;}
+ else {s[0]=-T(0,2); s[1]=-T(1,2); s[2]=-T(2,2); s[3]=-d;}
+ if(!propagateWithJacobian(useJac,1,s,P,Step)) return false;
+ }
+ else if (ty == Trk::Surface::Cylinder ) {
+
+ const Trk::CylinderSurface* cyl = static_cast<const Trk::CylinderSurface*>(su);
+
+ double r0[3] = {P[0],P[1],P[2]};
+ double s [9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),cyl->bounds().r(),m_direction,0.};
+
+ if(!propagateWithJacobian(useJac,2,s,P,Step)) return false;
+
+ // For cylinder we do test for next cross point
+ //
+ if(cyl->bounds().halfPhiSector() < 3.1 && newCrossPoint(*cyl,r0,P)) {
+ s[8] = 0.; if(!propagateWithJacobian(useJac,2,s,P,Step)) return false;
+ }
+ }
+ else if (ty == Trk::Surface::Perigee ) {
+
+ double s[6] ={T(0,3),T(1,3),T(2,3),0.,0.,1.};
+ if(!propagateWithJacobian(useJac,0,s,P,Step)) return false;
+ }
+ else if (ty == Trk::Surface::Cone ) {
+
+ double k = static_cast<const Trk::ConeSurface*>(su)->bounds().tanAlpha(); k = k*k+1.;
+ double s[9] = {T(0,3),T(1,3),T(2,3),T(0,2),T(1,2),T(2,2),k,m_direction,0.};
+ if(!propagateWithJacobian(useJac,3,s,P,Step)) return false;
+ }
+ else return false;
+
+ if(m_direction && (m_direction*Step)<0.) return false;
+
+ // Common transformation for all surfaces (angles and momentum)
+ //
+ if(useJac) {
+ double p=1./P[6]; P[35]*=p; P[36]*=p; P[37]*=p; P[38]*=p; P[39]*=p; P[40]*=p;
+ }
+
+ double p[5],Jac[21]; utils.transformGlobalToLocal(su,useJac,P,p,Jac);
+
+ // New simple track parameters production
+ //
+ Tb.setParameters(&Su,p);
+ if(useJac) {
+ Tb.newCovarianceMatrix(Ta,Jac);
+ const double* cv = Tb.cov();
+ if( cv[0]<=0. || cv[2]<=0. || cv[5]<=0. || cv[9]<=0. || cv[14]<=0.) return false;
+ }
+ return true;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Test new cross point
+/////////////////////////////////////////////////////////////////////////////////
+
+bool Trk::RungeKuttaPropagator::newCrossPoint
+(const Trk::CylinderSurface& Su,
+ const double * Ro,
+ const double * P ) const
+{
+ const double pi = 3.1415927, pi2=2.*pi;
+ const Amg::Transform3D& T = Su.transform();
+ double Ax[3] = {T(0,0),T(1,0),T(2,0)};
+ double Ay[3] = {T(0,1),T(1,1),T(2,1)};
+
+ double R = Su.bounds().r();
+ double x = Ro[0]-T(0,3);
+ double y = Ro[1]-T(1,3);
+ double z = Ro[2]-T(2,3);
+
+ double RC = x*Ax[0]+y*Ax[1]+z*Ax[2];
+ double RS = x*Ay[0]+y*Ay[1]+z*Ay[2];
+
+ if( (RC*RC+RS*RS) <= (R*R) ) return false;
+
+ x = P[0]-T(0,3);
+ y = P[1]-T(1,3);
+ z = P[2]-T(2,3);
+ RC = x*Ax[0]+y*Ax[1]+z*Ax[2];
+ RS = x*Ay[0]+y*Ay[1]+z*Ay[2];
+ double dF = fabs(atan2(RS,RC)-Su.bounds().averagePhi());
+ if(dF > pi) dF = pi2-pi;
+ if(dF <= Su.bounds().halfPhiSector()) return false;
+ return true;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Build new track parameters without propagation
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::buildTrackParametersWithoutPropagation
+(const Trk::TrackParameters& Tp,double* Jac) const
+{
+ Jac[0]=Jac[6]=Jac[12]=Jac[18]=Jac[20]=1.;
+ Jac[1]=Jac[2]=Jac[3]=Jac[4]=Jac[5]=Jac[7]=Jac[8]=Jac[9]=Jac[10]=Jac[11]=Jac[13]=Jac[14]=Jac[15]=Jac[16]=Jac[17]=Jac[19]=0.;
+ return Tp.clone();
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Build new neutral track parameters without propagation
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::NeutralParameters* Trk::RungeKuttaPropagator::buildTrackParametersWithoutPropagation
+(const Trk::NeutralParameters& Tp,double* Jac) const
+{
+ Jac[0]=Jac[6]=Jac[12]=Jac[18]=Jac[20]=1.;
+ Jac[1]=Jac[2]=Jac[3]=Jac[4]=Jac[5]=Jac[7]=Jac[8]=Jac[9]=Jac[10]=Jac[11]=Jac[13]=Jac[14]=Jac[15]=Jac[16]=Jac[17]=Jac[19]=0.;
+ return Tp.clone();
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Step estimator take into accout curvature
+/////////////////////////////////////////////////////////////////////////////////
+
+double Trk::RungeKuttaPropagator::stepEstimatorWithCurvature
+(int kind,double* Su,const double* P,bool& Q) const
+{
+ // Straight step estimation
+ //
+ Trk::RungeKuttaUtils utils;
+ double Step = utils.stepEstimator(kind,Su,P,Q); if(!Q) return 0.;
+ double AStep = fabs(Step);
+ if( kind || AStep < m_straightStep || !m_mcondition ) return Step;
+
+ const double* SA = &P[42]; // Start direction
+ double S = .5*Step;
+
+ double Ax = P[3]+S*SA[0];
+ double Ay = P[4]+S*SA[1];
+ double Az = P[5]+S*SA[2];
+ double As = 1./sqrt(Ax*Ax+Ay*Ay+Az*Az);
+ double PN[6] = {P[0],P[1],P[2],Ax*As,Ay*As,Az*As};
+ double StepN = utils.stepEstimator(kind,Su,PN,Q); if(!Q) {Q = true; return Step;}
+ if(fabs(StepN) < AStep) return StepN; return Step;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Global positions calculation inside CylinderBounds (one side)
+// where mS - max step allowed if mS > 0 propagate along momentum
+// mS < 0 propogate opposite momentum
+/////////////////////////////////////////////////////////////////////////////////
+
+void Trk::RungeKuttaPropagator::globalOneSidePositions
+(std::list<Amg::Vector3D> & GP,
+ const double * P,
+ const MagneticFieldProperties & M ,
+ const CylinderBounds & CB,
+ double mS,
+ ParticleHypothesis ) const
+{
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ double Pm[45]; for(int i=0; i!=7; ++i) Pm[i]=P[i];
+
+ double W = 0. ; // way
+ double R2max = CB.r()*CB.r() ; // max. radius**2 of region
+ double Zmax = CB.halflengthZ() ; // max. Z of region
+ double R2 = P[0]*P[0]+P[1]*P[1]; // Start radius**2
+ double Dir = P[0]*P[3]+P[1]*P[4]; // Direction
+ double S = mS ; // max step allowed
+ double R2m = R2 ;
+
+ if(m_mcondition && fabs(P[6]) > .1) return;
+
+ // Test position of the track
+ //
+ if(fabs(P[2]) > Zmax || R2 > R2max) return;
+
+ Amg::Vector3D g0(P[0],P[1],P[2]); GP.push_back(g0);
+
+ bool per = false; if(fabs(Dir)<.00001) per = true;
+ bool InS = false;
+
+ int niter = 0;
+ int sm = 1;
+ for(int s = 0; s!=2; ++s) {
+
+ if(s) {if(per) break; S = -S;}
+ double p[45]; for(int i=0; i!=7; ++i) p[i]=P[i]; p[42]=p[43]=p[44]=0.;
+
+ while(W<100000. && ++niter < 1000) {
+
+ if(m_mcondition) {
+ W+=(S=rungeKuttaStep (0,S,p,InS));
+ }
+ else {
+ W+=(S=straightLineStep(0, S,p));
+ }
+ if(InS && fabs(2.*S)<mS) S*=2.;
+
+ Amg::Vector3D g(p[0],p[1],p[2]);
+ if(!s) GP.push_back(g); else GP.push_front(g);
+
+ // Test position of the track
+ //
+ R2 = p[0]*p[0]+p[1]*p[1];
+
+ if(R2 < R2m) {
+ Pm[0]=p[0]; Pm[1]=p[1]; Pm[2]=p[2];
+ Pm[3]=p[3]; Pm[4]=p[4]; Pm[5]=p[5]; R2m = R2; sm = s;
+ }
+ if(fabs(p[2]) > Zmax || R2 > R2max) break;
+ if(!s && P[3]*p[3]+P[4]*p[4] < 0. ) break;
+
+ // Test perigee
+ //
+ if((p[0]*p[3]+p[1]*p[4])*Dir < 0.) {
+ if(s) break; per = true;
+ }
+ }
+ }
+
+ if(R2m < 400.) return;
+
+ // Propagate to perigee
+ //
+ Pm[42]=Pm[43]=Pm[44]=0.;
+ per = false;
+
+ for(int s = 0; s!=3; ++s) {
+
+ double A = (1.-Pm[5])*(1.+Pm[5]); if(A==0.) break;
+ S = -(Pm[0]*Pm[3]+Pm[1]*Pm[4])/A;
+ if(fabs(S) < 1. || ++niter > 1000) break;
+
+ if(m_mcondition) {
+ W+=(S=rungeKuttaStep (0,S,Pm,InS));
+ }
+ else {
+ W+=(S=straightLineStep(0,S,Pm));
+ }
+ per = true;
+ }
+
+ if(per) {
+ if(sm) {Amg::Vector3D gf(Pm[0],Pm[1],Pm[2]); GP.front() = gf;}
+ else {Amg::Vector3D gf(Pm[0],Pm[1],Pm[2]); GP.back () = gf;}
+ }
+ else {
+ double x = GP.front().x() , y = GP.front().y();
+ if( (x*x+y*y) > (Pm[0]*Pm[0]+Pm[1]*Pm[1]) ) {
+ if(sm) GP.pop_front();
+ else GP.pop_back ();
+ }
+ }
+ return;
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Global positions calculation inside CylinderBounds (one side)
+// where mS - max step allowed if mS > 0 propagate along momentum
+// mS < 0 propogate opposite momentum
+/////////////////////////////////////////////////////////////////////////////////
+
+void Trk::RungeKuttaPropagator::globalTwoSidePositions
+(std::list<Amg::Vector3D> & GP,
+ const double * P,
+ const MagneticFieldProperties & M ,
+ const CylinderBounds & CB,
+ double mS,
+ ParticleHypothesis ) const
+{
+ M.magneticFieldMode()==2 ? m_solenoid = true : m_solenoid = false;
+ M.magneticFieldMode()!=0 ? m_mcondition = true : m_mcondition = false;
+
+ double W = 0. ; // way
+ double R2max = CB.r()*CB.r() ; // max. radius**2 of region
+ double Zmax = CB.halflengthZ() ; // max. Z of region
+ double R2 = P[0]*P[0]+P[1]*P[1]; // Start radius**2
+ double S = mS ; // max step allowed
+
+ if(m_mcondition && fabs(P[6]) > .1) return;
+
+ // Test position of the track
+ //
+ if(fabs(P[2]) > Zmax || R2 > R2max) return;
+
+ Amg::Vector3D g0(P[0],P[1],P[2]); GP.push_back(g0);
+
+ bool InS = false;
+
+ int niter = 0;
+ for(int s = 0; s!=2; ++s) {
+
+ if(s) {S = -S;}
+ double p[45]; for(int i=0; i!=7; ++i) p[i]=P[i]; p[42]=p[43]=p[44]=0.;
+
+ while(W<100000. && ++niter < 1000) {
+
+ if(m_mcondition) {
+ W+=(S=rungeKuttaStep (0,S,p,InS));
+ }
+ else {
+ W+=(S=straightLineStep(0, S,p));
+ }
+ if(InS && fabs(2.*S)<mS) S*=2.;
+
+ Amg::Vector3D g(p[0],p[1],p[2]);
+ if(!s) GP.push_back(g); else GP.push_front(g);
+
+ // Test position of the track
+ //
+ R2 = p[0]*p[0]+p[1]*p[1];
+ if(fabs(p[2]) > Zmax || R2 > R2max) break;
+ }
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Track parameters in cross point preparation
+/////////////////////////////////////////////////////////////////////////////////
+
+const Trk::TrackParameters* Trk::RungeKuttaPropagator::crossPoint
+(const TrackParameters & Tp,
+ std::vector<DestSurf> & SU,
+ std::vector<unsigned int>& So,
+ double * P ,
+ std::pair<double,int> & SN) const
+{
+ double* R = &P[ 0] ; // Start coordinates
+ double* A = &P[ 3] ; // Start directions
+ double* SA = &P[42] ; // d(directions)/dStep
+ double Step = SN.first ;
+ int N = SN.second;
+
+ double As[3],Rs[3];
+
+ As[0] = A[0]+SA[0]*Step;
+ As[1] = A[1]+SA[1]*Step;
+ As[2] = A[2]+SA[2]*Step;
+
+ double CBA = 1./sqrt(As[0]*As[0]+As[1]*As[1]+As[2]*As[2]);
+
+ Rs[0] = R[0]+Step*(As[0]-.5*Step*SA[0]); As[0]*=CBA;
+ Rs[1] = R[1]+Step*(As[1]-.5*Step*SA[1]); As[1]*=CBA;
+ Rs[2] = R[2]+Step*(As[2]-.5*Step*SA[2]); As[2]*=CBA;
+
+ Amg::Vector3D pos(Rs[0],Rs[1],Rs[2]);
+ Amg::Vector3D dir(As[0],As[1],As[2]);
+
+ Trk::DistanceSolution ds = SU[N].first->straightLineDistanceEstimate(pos,dir,SU[N].second);
+ if(ds.currentDistance(false) > .010) return 0;
+
+ P[0] = Rs[0]; A[0] = As[0];
+ P[1] = Rs[1]; A[1] = As[1];
+ P[2] = Rs[2]; A[2] = As[2];
+
+ So.push_back(N);
+
+ // Transformation track parameters
+ //
+ Trk::RungeKuttaUtils utils;
+ bool useJac; Tp.covariance() ? useJac = true : useJac = false;
+
+ if(useJac) {
+ double d=1./P[6]; P[35]*=d; P[36]*=d; P[37]*=d; P[38]*=d; P[39]*=d; P[40]*=d;
+ }
+ double p[5],Jac[25];
+ utils.transformGlobalToLocal(SU[N].first,useJac,P,p,Jac);
+
+ if(!useJac) return SU[N].first->createTrackParameters(p[0],p[1],p[2],p[3],p[4],0);
+
+ AmgSymMatrix(5)* e = utils.newCovarianceMatrix(Jac,*Tp.covariance());
+ AmgSymMatrix(5)& cv = *e;
+
+ if(cv(0,0)<=0. || cv(1,1)<=0. || cv(2,2)<=0. || cv(3,3)<=0. || cv(4,4)<=0.) {
+ delete e; return 0;
+ }
+ return SU[N].first->createTrackParameters(p[0],p[1],p[2],p[3],p[4],e);
+}
+
+/////////////////////////////////////////////////////////////////////////////////
+// Step reduction from STEP propagator
+// Input information ERRx,ERRy,ERRz,current step
+/////////////////////////////////////////////////////////////////////////////////
+
+double Trk::RungeKuttaPropagator::stepReduction(const double* E) const
+{
+ double dlt2 = m_dlt*m_dlt;
+ double dR2 = (E[0]*E[0]+E[1]*E[1]+E[2]*E[2])*E[3]*E[3]*4.; if(dR2 < 1.e-40) dR2 = 1.e-40;
+ double cS = std::pow(dlt2/dR2,0.125);
+ if(cS < .25) cS = .25;
+ if((dR2 > 16.*dlt2 && cS < 1.) || cS >=3.) return cS;
+ return 1.;
+}
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_entries.cxx b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_entries.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..3163d9b61f48ef6cbd3efe000ce29e36490c33c9
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_entries.cxx
@@ -0,0 +1,13 @@
+#include "GaudiKernel/DeclareFactoryEntries.h"
+#include "TrkExRungeKuttaPropagator/RungeKuttaPropagator.h"
+
+using namespace Trk;
+
+DECLARE_TOOL_FACTORY( RungeKuttaPropagator )
+
+/** factory entries need to have the name of the package */
+DECLARE_FACTORY_ENTRIES( TrkExRungeKuttaPropagator )
+{
+ DECLARE_TOOL( RungeKuttaPropagator )
+}
+
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_load.cxx b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_load.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..88ab7fe92bfd6fe9a2057ec7eca8a7a7b4a44808
--- /dev/null
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaPropagator/src/components/TrkExRungeKuttaPropagator_load.cxx
@@ -0,0 +1,3 @@
+#include "GaudiKernel/LoadFactoryEntries.h"
+
+LOAD_FACTORY_ENTRIES( TrkExRungeKuttaPropagator )