diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/TrkExRungeKuttaIntersector/IntersectorWrapper.h b/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/TrkExRungeKuttaIntersector/IntersectorWrapper.h
index c8b2115419ae40dd028b41de9a0ae533f8f162e7..c91108a7c3cf6f51776e8743fcf6f9ce5c922ab9 100755
--- a/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/TrkExRungeKuttaIntersector/IntersectorWrapper.h
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/TrkExRungeKuttaIntersector/IntersectorWrapper.h
@@ -1,214 +1,224 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+ */
+
+//////////////////////////////////////////////////////////////////////
+// wrapper to IIntersector tool to provide IPropagator functionality
+// default configuration wraps the RungeKutta intersector
+// (c) ATLAS Detector software
+//////////////////////////////////////////////////////////////////////
+
+#ifndef TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
+#define TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "GeoPrimitives/GeoPrimitives.h"
+#include "TrkExInterfaces/IPropagator.h"
+
+namespace Trk
+{
+class IIntersector;
+class TrackSurfaceIntersection;
+
+class IntersectorWrapper: public AthAlgTool,
+ virtual public IPropagator
+ {
-//////////////////////////////////////////////////////////////////////
-// wrapper to IIntersector tool to provide IPropagator functionality
-// default configuration wraps the RungeKutta intersector
-// (c) ATLAS Detector software
-//////////////////////////////////////////////////////////////////////
-
-#ifndef TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
-#define TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
-
-#include "AthenaBaseComps/AthAlgTool.h"
-#include "GaudiKernel/ToolHandle.h"
-#include "GeoPrimitives/GeoPrimitives.h"
-#include "TrkExInterfaces/IPropagator.h"
-
-namespace Trk
-{
-class IIntersector;
-class TrackSurfaceIntersection;
-
-class IntersectorWrapper: public AthAlgTool,
- virtual public IPropagator
-{
-
-public:
+ public:
using IPropagator::propagate;
using IPropagator::propagateT;
- IntersectorWrapper (const std::string& type,
- const std::string& name,
- const IInterface* parent);
- ~IntersectorWrapper (void); // destructor
-
- StatusCode initialize();
- StatusCode finalize();
-
-
-
- /** N 0) <b>Neutral parameters method </b>
- - returns a ParametersBase object as well, 0 if the extrapolation did not succeed
- */
- /// implemented
- const NeutralParameters* propagate(const NeutralParameters&,
- const Surface&,
- PropDirection,
- BoundaryCheck,
- bool) const;
-
- /** [TrackParameters] --------------------------------------------------------- */
-
- /** Propagation interface:
-
- The propagation method called by the TrkExtrapolator. The extrapolator
- is responsible for the underlying logic of which surface to go to.
- */
- /// implemented
- const TrackParameters* propagate( const TrackParameters& parm,
- const Surface& sf,
- PropDirection dir,
- BoundaryCheck bcheck,
- const MagneticFieldProperties& mprop,
- ParticleHypothesis particle,
- bool returnCurv,
- const TrackingVolume*) const;
-
-
- /** Propagation interface:
-
- The propagation method called by the TrkExtrapolator. The propagator
- 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 { return 0; }
-
- /** Propagation interface:
-
- The propagation method called by the TrkExtrapolator. The propagator
- finds the closest surface. Timing included.
- */
- const TrackParameters* propagateT( const TrackParameters&,
- std::vector<DestSurf>&,
- PropDirection,
- const MagneticFieldProperties&,
- ParticleHypothesis,
- std::vector<unsigned int>&,
- PathLimit&, TimeLimit&,
- bool,
- const TrackingVolume*,
- std::vector<Trk::HitInfo>*&) const { return 0; }
-
-
- /** Propagation interface:
-
- The propagation method including the return of the TransportJacobian matrix.
-
- */
- /// implemented
- const TrackParameters* propagate( const TrackParameters&,
- const Surface&,
- PropDirection,
- BoundaryCheck,
- const MagneticFieldProperties&,
- TransportJacobian*&,
- double&,
- ParticleHypothesis,
- bool,
- const TrackingVolume*) const;
-
-
- /** Propagation interface without Covariance matrix propagation
- the pathlength has to be returned for eventual following propagateCovariance
- */
- /// implemented
- const TrackParameters* propagateParameters( const TrackParameters& parm,
- const Surface& sf,
- PropDirection dir,
- BoundaryCheck bcheck,
- const MagneticFieldProperties& mprop,
- ParticleHypothesis particle=pion,
- bool returnCurv = false,
- const TrackingVolume* tVol=0) const;
-
- /// implemented
- const TrackParameters* propagateParameters( const TrackParameters& parm,
- const Surface& sf,
- PropDirection dir,
- BoundaryCheck bcheck,
- const MagneticFieldProperties& mprop,
- TransportJacobian*&,
- ParticleHypothesis particle=pion,
- bool returnCurv = false,
- const TrackingVolume* tVol=0) const;
-
-
- /** Intersection interface:
-
- The intersection interface might be used by the material service as well to estimate
- the surfaces (sensitive and nonesensitive) while propagation
- */
- const IntersectionSolution* intersect( const TrackParameters& parm,
- const Surface& sf,
- const MagneticFieldProperties& mprop,
- ParticleHypothesis particle=pion,
- const TrackingVolume* tVol=0) const;
-
- /** GlobalPositions list interface:
- This is used mostly in pattern recognition in the road finder, the propagation direction is intrinsically given
- by the sign of the stepSize.
-
- To avoid memory fragmentation in multiple use of pattern recognition processes and respecting the possible iterative
- filling of the positions list, the list of GlobalPositions is given by reference through the signature and a void
- method has been chosen.
- */
- void globalPositions(std::list<Amg::Vector3D>& positionslist,
- const TrackParameters& parm,
- const MagneticFieldProperties& mprop,
- const CylinderBounds& cylbo,
- double stepSize,
- ParticleHypothesis particle=pion,
- const TrackingVolume* tVol=0) const;
-
- //placeholder for compatibility with new interface
- const TrackSurfaceIntersection* intersectSurface(const Surface&,
- const TrackSurfaceIntersection*,
- const double,
- const MagneticFieldProperties&,
- ParticleHypothesis) const {return 0;}
-
- /** Validation Action:
- Can be implemented optionally, outside access to internal validation steps */
- void validationAction() const {}
-
-
-
-private:
-
- // private methods
- void
- createParameters (const Surface& surface,
- BoundaryCheck boundsCheck,
- bool curvilinear) const;
- void
- findIntersection (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir = Trk::anyDirection) const;
-
- // helpers, managers, tools
- ToolHandle<IIntersector> m_intersector;
- ToolHandle<IPropagator> m_linePropagator;
-
- // current parameter cache
- mutable double m_charge;
- mutable const TrackSurfaceIntersection* m_intersection;
- mutable Amg::Vector3D m_momentum;
- mutable const TrackParameters* m_parameters;
- mutable Amg::Vector3D m_position;
- mutable double m_qOverP;
-};
-
-} // end of namespace
-
-
-#endif // TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
+ IntersectorWrapper (const std::string& type,
+ const std::string& name,
+ const IInterface* parent);
+ ~IntersectorWrapper (void); // destructor
+
+ StatusCode initialize();
+ StatusCode finalize();
+
+
+
+ /** N 0) <b>Neutral parameters method </b>
+ - returns a ParametersBase object as well, 0 if the extrapolation did not succeed
+ */
+ /// implemented
+ virtual const NeutralParameters* propagate(const NeutralParameters&,
+ const Surface&,
+ PropDirection,
+ BoundaryCheck,
+ bool) const override;
+
+ /** [TrackParameters] --------------------------------------------------------- */
+
+ /** Propagation interface:
+
+ The propagation method called by the TrkExtrapolator. The extrapolator
+ is responsible for the underlying logic of which surface to go to.
+ */
+ /// implemented
+ virtual const TrackParameters* propagate( const TrackParameters& parm,
+ const Surface& sf,
+ PropDirection dir,
+ BoundaryCheck bcheck,
+ const MagneticFieldProperties& mprop,
+ ParticleHypothesis particle,
+ bool returnCurv,
+ const TrackingVolume*) const override;
+
+
+ /** Propagation interface:
+
+ The propagation method called by the TrkExtrapolator. The propagator
+ finds the closest surface.
+ */
+ virtual const TrackParameters* propagate( const TrackParameters&,
+ std::vector<DestSurf>&,
+ PropDirection,
+ const MagneticFieldProperties&,
+ ParticleHypothesis,
+ std::vector<unsigned int>&,
+ double&,
+ bool,
+ bool,
+ const TrackingVolume*) const override{ return 0; }
+
+ /** Propagation interface:
+
+ The propagation method called by the TrkExtrapolator. The propagator
+ finds the closest surface. Timing included.
+ */
+ virtual const TrackParameters* propagateT( const TrackParameters&,
+ std::vector<DestSurf>&,
+ PropDirection,
+ const MagneticFieldProperties&,
+ ParticleHypothesis,
+ std::vector<unsigned int>&,
+ PathLimit&, TimeLimit&,
+ bool,
+ const TrackingVolume*,
+ std::vector<Trk::HitInfo>*&) const override{ return 0; }
+
+
+ /** Propagation interface:
+
+ The propagation method including the return of the TransportJacobian matrix.
+
+*/
+ /// implemented
+ virtual const TrackParameters* propagate( const TrackParameters&,
+ const Surface&,
+ PropDirection,
+ BoundaryCheck,
+ const MagneticFieldProperties&,
+ TransportJacobian*&,
+ double&,
+ ParticleHypothesis,
+ bool,
+ const TrackingVolume*) const override;
+
+
+ /** Propagation interface without Covariance matrix propagation
+ the pathlength has to be returned for eventual following propagateCovariance
+ */
+ /// implemented
+ virtual const TrackParameters* propagateParameters( const TrackParameters& parm,
+ const Surface& sf,
+ PropDirection dir,
+ BoundaryCheck bcheck,
+ const MagneticFieldProperties& mprop,
+ ParticleHypothesis particle=pion,
+ bool returnCurv = false,
+ const TrackingVolume* tVol=0) const override;
+
+ /// implemented
+ virtual const TrackParameters* propagateParameters( const TrackParameters& parm,
+ const Surface& sf,
+ PropDirection dir,
+ BoundaryCheck bcheck,
+ const MagneticFieldProperties& mprop,
+ TransportJacobian*&,
+ ParticleHypothesis particle=pion,
+ bool returnCurv = false,
+ const TrackingVolume* tVol=0) const override;
+
+
+ /** Intersection interface:
+
+ The intersection interface might be used by the material service as well to estimate
+ the surfaces (sensitive and nonesensitive) while propagation
+ */
+ virtual const IntersectionSolution* intersect( const TrackParameters& parm,
+ const Surface& sf,
+ const MagneticFieldProperties& mprop,
+ ParticleHypothesis particle=pion,
+ const TrackingVolume* tVol=0) const override;
+
+ /** GlobalPositions list interface:
+ This is used mostly in pattern recognition in the road finder, the propagation direction is intrinsically given
+ by the sign of the stepSize.
+
+ To avoid memory fragmentation in multiple use of pattern recognition processes and respecting the possible iterative
+ filling of the positions list, the list of GlobalPositions is given by reference through the signature and a void
+ method has been chosen.
+ */
+ virtual void globalPositions(std::list<Amg::Vector3D>& positionslist,
+ const TrackParameters& parm,
+ const MagneticFieldProperties& mprop,
+ const CylinderBounds& cylbo,
+ double stepSize,
+ ParticleHypothesis particle=pion,
+ const TrackingVolume* tVol=0) const override;
+
+ //placeholder for compatibility with new interface
+ virtual const TrackSurfaceIntersection* intersectSurface(const Surface&,
+ const TrackSurfaceIntersection*,
+ const double,
+ const MagneticFieldProperties&,
+ ParticleHypothesis) const override{return 0;}
+
+ /** Validation Action:
+ Can be implemented optionally, outside access to internal validation steps */
+ virtual void validationAction() const override{}
+
+
+
+ private:
+ struct Cache{
+ double m_charge;
+ double m_qOverP;
+ std::unique_ptr<const TrackSurfaceIntersection> m_intersection;
+ const TrackParameters* m_parameters;
+ Amg::Vector3D m_position;
+ Amg::Vector3D m_momentum;
+
+ Cache():
+ m_charge{},
+ m_qOverP{},
+ m_intersection{nullptr},
+ m_parameters{nullptr},
+ m_position{},
+ m_momentum{}{
+ }
+ };
+ // private methods
+ void
+ createParameters (Cache& cache,
+ const Surface& surface,
+ BoundaryCheck boundsCheck,
+ bool curvilinear) const;
+ void
+ findIntersection (Cache& cache,
+ const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir = Trk::anyDirection) const;
+
+ // helpers, managers, tools
+ ToolHandle<IIntersector> m_intersector;
+ ToolHandle<IPropagator> m_linePropagator;
+ };
+
+} // end of namespace
+
+
+#endif // TRKEXRUNGEKUTTAINTERSECTOR_INTERSECTORWRAPPER_H
diff --git a/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/src/IntersectorWrapper.cxx b/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/src/IntersectorWrapper.cxx
index db3c53eddd4ad35c25fee262182b33de42655feb..d1d40d98d0a6d710edb510d8e8c78cc175919ff0 100755
--- a/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/src/IntersectorWrapper.cxx
+++ b/Tracking/TrkExtrapolation/TrkExRungeKuttaIntersector/src/IntersectorWrapper.cxx
@@ -1,258 +1,260 @@
/*
- Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
-*/
+ Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+ */
+
+//////////////////////////////////////////////////////////////////////
+// wrapper to IIntersector tool to provide IPropagator functionality
+// default configuration wraps the RungeKutta intersector
+// (c) ATLAS Detector software
+//////////////////////////////////////////////////////////////////////
+
+#include <iomanip>
+#include <vector>
+#include "GaudiKernel/SystemOfUnits.h"
+#include "TrkExInterfaces/IIntersector.h"
+#include "TrkExRungeKuttaIntersector/IntersectorWrapper.h"
+#include "TrkExUtils/IntersectionSolution.h"
+#include "TrkExUtils/TrackSurfaceIntersection.h"
+#include "TrkGeometry/MagneticFieldProperties.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"
+#include "TrkParameters/TrackParameters.h"
+
+namespace Trk
+{
+
+//<<<<<< CLASS STRUCTURE INITIALIZATION >>>>>>
+
+IntersectorWrapper::IntersectorWrapper (const std::string& type,
+ const std::string& name,
+ const IInterface* parent)
+ :AthAlgTool (type, name, parent),
+ m_intersector ("Trk::RungeKuttaIntersector/RungeKuttaIntersector"),
+ m_linePropagator ("")
+ {
+ declareInterface<Trk::IPropagator>(this);
+ declareProperty("Intersector", m_intersector);
+ declareProperty("LinePropagator", m_linePropagator);
+ }
+
+IntersectorWrapper::~IntersectorWrapper (void)
+{}
+
+//<<<<<< PUBLIC MEMBER FUNCTION DEFINITIONS >>>>>>
+
+StatusCode
+IntersectorWrapper::initialize(){
+ // get the Tools
+
+ ATH_CHECK(m_intersector.retrieve());
+
+ ATH_MSG_VERBOSE( "Retrieved tool " << m_intersector );
+
+ if (!m_linePropagator.empty()) {
+ ATH_CHECK( m_linePropagator.retrieve());
+ }
+ return StatusCode::SUCCESS;
+}
+
+
+StatusCode IntersectorWrapper::finalize()
+{
+ return StatusCode::SUCCESS;
+}
+
+const NeutralParameters*
+IntersectorWrapper::propagate (const NeutralParameters& parameters,
+ const Surface& surface,
+ PropDirection dir,
+ BoundaryCheck boundsCheck,
+ bool curvilinear) const
+{
+ return m_linePropagator->propagate(parameters,surface,dir,boundsCheck,curvilinear);
+}
+
+const TrackParameters*
+IntersectorWrapper::propagate (const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir,
+ BoundaryCheck boundsCheck,
+ const MagneticFieldProperties& /*magProperties*/,
+ ParticleHypothesis /*particle*/,
+ bool curvilinear,
+ const TrackingVolume*) const
+{
+ Cache cache{};
+ findIntersection(cache,parameters,surface, dir);
+ createParameters(cache,surface,boundsCheck,curvilinear);
+ return cache.m_parameters;
+}
+
+const TrackParameters*
+IntersectorWrapper::propagate (const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir,
+ BoundaryCheck boundsCheck,
+ const MagneticFieldProperties& /*magProperties*/,
+ TransportJacobian*& /*transportJac*/,
+ double&,
+ ParticleHypothesis /*particle*/,
+ bool curvilinear,
+ const TrackingVolume*) const
+{
+ Cache cache{};
+ findIntersection(cache,parameters,surface,dir);
+ createParameters(cache,surface,boundsCheck,curvilinear);
+ return cache.m_parameters;
+}
+
+const TrackParameters*
+IntersectorWrapper::propagateParameters (const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir,
+ BoundaryCheck boundsCheck,
+ const MagneticFieldProperties& /*magProperties*/,
+ ParticleHypothesis /*particle*/,
+ bool curvilinear,
+ const TrackingVolume*) const
+{
+
+ Cache cache{};
+ findIntersection(cache,parameters,surface,dir);
+ createParameters(cache,surface,boundsCheck,curvilinear);
+ return cache.m_parameters;
+}
+
+const TrackParameters*
+IntersectorWrapper::propagateParameters (const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir,
+ BoundaryCheck boundsCheck,
+ const MagneticFieldProperties& /*magProperties*/,
+ TransportJacobian*& /*transportJac*/,
+ ParticleHypothesis /*particle*/,
+ bool curvilinear,
+ const TrackingVolume*) const
+{
+ Cache cache{};
+ findIntersection(cache,parameters,surface,dir);
+ createParameters(cache,surface,boundsCheck,curvilinear);
+ return cache.m_parameters;
+}
+
+const IntersectionSolution*
+IntersectorWrapper::intersect (const TrackParameters& parameters,
+ const Surface& surface,
+ const MagneticFieldProperties& /*magProperties*/,
+ ParticleHypothesis /*particle*/,
+ const TrackingVolume*) const
+{
+ Cache cache{};
+ findIntersection(cache,parameters,surface);
+ IntersectionSolution* solution = new IntersectionSolution;
+ if (cache.m_intersection.get()) {
+ solution->push_back(cache.m_intersection.release());
+ }
+ return solution;
+}
+
+void
+IntersectorWrapper::globalPositions (std::list<Amg::Vector3D>&,
+ const TrackParameters&,
+ const MagneticFieldProperties&,
+ const CylinderBounds&,
+ double,
+ ParticleHypothesis,
+ const TrackingVolume* ) const
+{
+ ATH_MSG_WARNING( "globalPositions method from the IPropagator interface"
+ << " has NOT been implememnted" );
+}
+
+//<<<<<< PRIVATE MEMBER FUNCTION DEFINITIONS >>>>>>
+
+void
+IntersectorWrapper::createParameters (Cache& cache,
+ const Surface& surface,
+ BoundaryCheck /*boundsCheck*/,
+ bool curvilinear) const
+{
+ cache.m_parameters = nullptr;
+ if (! cache.m_intersection) return;
+
+ // curvilinear special (simple) case
+ if (curvilinear)
+ {
+ cache.m_parameters=new CurvilinearParameters(cache.m_intersection->position(),
+ cache.m_intersection->direction().phi(),
+ cache.m_intersection->direction().theta(),
+ cache.m_qOverP);
+ return;
+ }
+
+ cache.m_parameters=surface.createTrackParameters(cache.m_intersection->position(),
+ cache.m_intersection->direction(),
+ cache.m_charge,nullptr);
+ // unrecognized Surface
+ if( !cache.m_parameters ) ATH_MSG_WARNING( " Failed to create parameters " );
+}
+
+void
+IntersectorWrapper::findIntersection (Cache& cache,
+ const TrackParameters& parameters,
+ const Surface& surface,
+ PropDirection dir) const
+{
+ cache.m_charge = parameters.charge();
+ cache.m_momentum = parameters.momentum();
+ cache.m_parameters=¶meters;
+ cache.m_position = Amg::Vector3D(parameters.position());
+ cache.m_qOverP = 1./cache.m_momentum.mag();
+ cache.m_intersection.reset(new TrackSurfaceIntersection(cache.m_position,cache.m_momentum*cache.m_qOverP,0.));
+ cache.m_qOverP *= cache.m_charge;
+
+ const TrackSurfaceIntersection* oldIntersection = cache.m_intersection.release();
+ cache.m_intersection.reset( m_intersector->intersectSurface(surface,
+ oldIntersection,
+ cache.m_qOverP));
+ delete oldIntersection;
+ if (! cache.m_intersection.get())
+ {
+ ATH_MSG_DEBUG( " no intersection found " );
+ return;
+ }
+
+ // check for correct propagation direction
+ if ((dir == Trk::alongMomentum && cache.m_intersection->pathlength() < 0.)
+ || (dir == Trk::oppositeMomentum && cache.m_intersection->pathlength() > 0.))
+ {
+ if (msgLvl(MSG::DEBUG))
+ {
+ if (dir == Trk::alongMomentum)
+ {
+ ATH_MSG_DEBUG( " requested alongMomentum, but pathlength "
+ << cache.m_intersection->pathlength() );
+ }
+ else if (dir == Trk::oppositeMomentum)
+ {
+ ATH_MSG_DEBUG( " requested oppositeMomentum, but pathlength "
+ << cache.m_intersection->pathlength() );
+ }
+ }
+ cache.m_intersection=nullptr;
+ return;
+ }
+
+ ATH_MSG_DEBUG( std::setiosflags(std::ios::fixed)
+ << " intersection at r,phi,z "
+ << std::setw(10) << std::setprecision(1) << cache.m_intersection->position().perp()
+ << std::setw(9) << std::setprecision(4) << cache.m_intersection->position().phi()
+ << std::setw(10) << std::setprecision(1) << cache.m_intersection->position().z()
+ << " momentum "
+ << std::setw(9) << std::setprecision(3) << 1./fabs(cache.m_qOverP*Gaudi::Units::GeV) );
+}
+
+} // end of namespace
-//////////////////////////////////////////////////////////////////////
-// wrapper to IIntersector tool to provide IPropagator functionality
-// default configuration wraps the RungeKutta intersector
-// (c) ATLAS Detector software
-//////////////////////////////////////////////////////////////////////
-
-#include <iomanip>
-#include <vector>
-#include "GaudiKernel/SystemOfUnits.h"
-#include "TrkExInterfaces/IIntersector.h"
-#include "TrkExRungeKuttaIntersector/IntersectorWrapper.h"
-#include "TrkExUtils/IntersectionSolution.h"
-#include "TrkExUtils/TrackSurfaceIntersection.h"
-#include "TrkGeometry/MagneticFieldProperties.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"
-#include "TrkParameters/TrackParameters.h"
-
-namespace Trk
-{
-
-//<<<<<< CLASS STRUCTURE INITIALIZATION >>>>>>
-
-IntersectorWrapper::IntersectorWrapper (const std::string& type,
- const std::string& name,
- const IInterface* parent)
- : AthAlgTool (type, name, parent),
- m_intersector ("Trk::RungeKuttaIntersector/RungeKuttaIntersector"),
- m_linePropagator (""),
- m_charge{},
- m_intersection (nullptr),
- m_momentum{},
- m_parameters (nullptr),
- m_position{},
- m_qOverP{}
-{
- declareInterface<Trk::IPropagator>(this);
- declareProperty("Intersector", m_intersector);
- declareProperty("LinePropagator", m_linePropagator);
-}
-
-IntersectorWrapper::~IntersectorWrapper (void)
-{}
-
-//<<<<<< PUBLIC MEMBER FUNCTION DEFINITIONS >>>>>>
-
-StatusCode
-IntersectorWrapper::initialize(){
- // get the Tools
- ATH_CHECK(m_intersector.retrieve());
- ATH_MSG_VERBOSE( "Retrieved tool " << m_intersector );
- if (!m_linePropagator.empty()) ATH_CHECK( m_linePropagator.retrieve());
- return StatusCode::SUCCESS;
-}
-
-StatusCode
-IntersectorWrapper::finalize()
-{
- delete m_intersection;
-
- return StatusCode::SUCCESS;
-}
-
-const NeutralParameters*
-IntersectorWrapper::propagate (const NeutralParameters& parameters,
- const Surface& surface,
- PropDirection dir,
- BoundaryCheck boundsCheck,
- bool curvilinear) const
-{
- return m_linePropagator->propagate(parameters,surface,dir,boundsCheck,curvilinear);
-}
-
-const TrackParameters*
-IntersectorWrapper::propagate (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir,
- BoundaryCheck boundsCheck,
- const MagneticFieldProperties& /*magProperties*/,
- ParticleHypothesis /*particle*/,
- bool curvilinear,
- const TrackingVolume*) const
-{
- findIntersection(parameters,surface, dir);
- createParameters(surface,boundsCheck,curvilinear);
- return m_parameters;
-}
-
-const TrackParameters*
-IntersectorWrapper::propagate (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir,
- BoundaryCheck boundsCheck,
- const MagneticFieldProperties& /*magProperties*/,
- TransportJacobian*& /*transportJac*/,
- double&,
- ParticleHypothesis /*particle*/,
- bool curvilinear,
- const TrackingVolume*) const
-{
- findIntersection(parameters,surface,dir);
- createParameters(surface,boundsCheck,curvilinear);
- return m_parameters;
-}
-
-const TrackParameters*
-IntersectorWrapper::propagateParameters (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir,
- BoundaryCheck boundsCheck,
- const MagneticFieldProperties& /*magProperties*/,
- ParticleHypothesis /*particle*/,
- bool curvilinear,
- const TrackingVolume*) const
-{
-
- findIntersection(parameters,surface,dir);
- createParameters(surface,boundsCheck,curvilinear);
- return m_parameters;
-}
-
-const TrackParameters*
-IntersectorWrapper::propagateParameters (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir,
- BoundaryCheck boundsCheck,
- const MagneticFieldProperties& /*magProperties*/,
- TransportJacobian*& /*transportJac*/,
- ParticleHypothesis /*particle*/,
- bool curvilinear,
- const TrackingVolume*) const
-{
- findIntersection(parameters,surface,dir);
- createParameters(surface,boundsCheck,curvilinear);
- return m_parameters;
-}
-
-const IntersectionSolution*
-IntersectorWrapper::intersect (const TrackParameters& parameters,
- const Surface& surface,
- const MagneticFieldProperties& /*magProperties*/,
- ParticleHypothesis /*particle*/,
- const TrackingVolume*) const
-{
- findIntersection(parameters,surface);
- IntersectionSolution* solution = new IntersectionSolution;
- if (m_intersection) solution->push_back(m_intersection);
- return solution;
-}
-
-void
-IntersectorWrapper::globalPositions (std::list<Amg::Vector3D>&,
- const TrackParameters&,
- const MagneticFieldProperties&,
- const CylinderBounds&,
- double,
- ParticleHypothesis,
- const TrackingVolume* ) const
-{
- ATH_MSG_WARNING( "globalPositions method from the IPropagator interface"
- << " has NOT been implememnted" );
-}
-
-//<<<<<< PRIVATE MEMBER FUNCTION DEFINITIONS >>>>>>
-
-void
-IntersectorWrapper::createParameters (const Surface& surface,
- BoundaryCheck /*boundsCheck*/,
- bool curvilinear) const
-{
- m_parameters = 0;
- if (! m_intersection) return;
-
- // curvilinear special (simple) case
- if (curvilinear)
- {
- m_parameters = new CurvilinearParameters(m_intersection->position(),
- m_intersection->direction().phi(),
- m_intersection->direction().theta(),
- m_qOverP);
- return;
- }
-
- m_parameters = surface.createTrackParameters(m_intersection->position(),m_intersection->direction(),m_charge,nullptr);
- // unrecognized Surface
- if( !m_parameters ) ATH_MSG_WARNING( " Failed to create parameters " );
-}
-
-void
-IntersectorWrapper::findIntersection (const TrackParameters& parameters,
- const Surface& surface,
- PropDirection dir) const
-{
- //if (! m_intersection || ¶meters != m_parameters)
- {
- delete m_intersection;
- m_charge = parameters.charge();
- m_momentum = parameters.momentum();
- m_parameters = ¶meters;
- m_position = parameters.position();
- m_position = Amg::Vector3D(parameters.position());
- m_qOverP = 1./m_momentum.mag();
- m_intersection = new TrackSurfaceIntersection(m_position,m_momentum*m_qOverP,0.);
- m_qOverP *= m_charge;
- }
-
- const TrackSurfaceIntersection* oldIntersection = m_intersection;
- m_intersection = m_intersector->intersectSurface(surface,
- oldIntersection,
- m_qOverP);
- delete oldIntersection;
- if (! m_intersection)
- {
- ATH_MSG_DEBUG( " no intersection found " );
- return;
- }
-
- // check for correct propagation direction
- if ((dir == Trk::alongMomentum && m_intersection->pathlength() < 0.)
- || (dir == Trk::oppositeMomentum && m_intersection->pathlength() > 0.))
- {
- if (msgLvl(MSG::DEBUG))
- {
- if (dir == Trk::alongMomentum)
- {
- ATH_MSG_DEBUG( " requested alongMomentum, but pathlength "
- << m_intersection->pathlength() );
- }
- else if (dir == Trk::oppositeMomentum)
- {
- ATH_MSG_DEBUG( " requested oppositeMomentum, but pathlength "
- << m_intersection->pathlength() );
- }
- }
- delete m_intersection;
- m_intersection = 0;
- return;
- }
-
- ATH_MSG_DEBUG( std::setiosflags(std::ios::fixed)
- << " intersection at r,phi,z "
- << std::setw(10) << std::setprecision(1) << m_intersection->position().perp()
- << std::setw(9) << std::setprecision(4) << m_intersection->position().phi()
- << std::setw(10) << std::setprecision(1) << m_intersection->position().z()
- << " momentum "
- << std::setw(9) << std::setprecision(3) << 1./fabs(m_qOverP*Gaudi::Units::GeV) );
-}
-
-} // end of namespace
-