#include <AnalyticalPropagator.h>

Inheritance diagram for AnalyticalPropagator:

Public Member Functions
	AnalyticalPropagator (const MagneticField *field, PropagationDirection dir=alongMomentum, float maxDPhi=1.6, bool isOld=true)

virtual AnalyticalPropagator *	clone () const

virtual bool	setMaxDirectionChange (float phiMax)

void	setMaxRelativeChangeInBz (const float maxDBz)

	~AnalyticalPropagator ()

Public Member Functions inherited from Propagator
template<typename STA , typename SUR >
TrajectoryStateOnSurface	propagate (STA const &state, SUR const &surface) const

virtual FreeTrajectoryState	propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const final

virtual FreeTrajectoryState	propagate (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const final

virtual FreeTrajectoryState	propagate (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const final

virtual std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const FreeTrajectoryState &, const Surface &) const final

virtual std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const TrajectoryStateOnSurface &tsos, const Surface &sur) const final

virtual std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const TrajectoryStateOnSurface &tsos, const Plane &sur) const

virtual std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const TrajectoryStateOnSurface &tsos, const Cylinder &sur) const

virtual std::pair < FreeTrajectoryState, double >	propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest) const

virtual std::pair < FreeTrajectoryState, double >	propagateWithPath (const FreeTrajectoryState &ftsStart, const GlobalPoint &pDest1, const GlobalPoint &pDest2) const
	Propagate to PCA to a line (given by 2 points) given a starting point. More...

virtual std::pair < FreeTrajectoryState, double >	propagateWithPath (const FreeTrajectoryState &ftsStart, const reco::BeamSpot &beamSpot) const
	Propagate to PCA to a line (given by beamSpot position and slope) given a starting point. More...

virtual PropagationDirection	propagationDirection () const final

	Propagator (PropagationDirection dir=alongMomentum)

virtual void	setPropagationDirection (PropagationDirection dir)

virtual	~Propagator ()

Private Types
typedef std::pair < TrajectoryStateOnSurface, double >	TsosWP

Private Member Functions
virtual const MagneticField *	magneticField () const

std::pair < TrajectoryStateOnSurface, double >	propagatedStateWithPath (const FreeTrajectoryState &fts, const Surface &surface, const GlobalTrajectoryParameters &gtp, const double &s) const
	propagation of errors (if needed) and generation of a new TSOS More...

bool	propagateParametersOnCylinder (const FreeTrajectoryState &fts, const Cylinder &cylinder, GlobalPoint &x, GlobalVector &p, double &s) const
	parameter propagation to cylinder (returns position, momentum and path length) More...

bool	propagateParametersOnPlane (const FreeTrajectoryState &fts, const Plane &plane, GlobalPoint &x, GlobalVector &p, double &s) const
	parameter propagation to plane (returns position, momentum and path length) More...

bool	propagateWithHelixCrossing (HelixPlaneCrossing &, const Plane &, const float, GlobalPoint &, GlobalVector &, double &s) const
	helix parameter propagation to a plane using HelixPlaneCrossing More...

bool	propagateWithLineCrossing (const GlobalPoint &, const GlobalVector &, const Plane &, GlobalPoint &, double &) const
	straight line parameter propagation to a plane More...

bool	propagateWithLineCrossing (const GlobalPoint &, const GlobalVector &, const Cylinder &, GlobalPoint &, double &) const
	straight line parameter propagation to a cylinder More...

std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const FreeTrajectoryState &fts, const Plane &plane) const override
	propagation to plane with path length More...

std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const FreeTrajectoryState &fts, const Cylinder &cylinder) const override
	propagation to cylinder with path length More...

Private Attributes
bool	isOldPropagationType

const MagneticField *	theField

float	theMaxDBzRatio

float	theMaxDPhi2

Detailed Description

(Mostly) analytical helix propagation to cylindrical or planar surfaces. Based on GtfGeometricalPropagator with successive replacement of components (currently: propagation to arbitrary plane).

Definition at line 23 of file AnalyticalPropagator.h.

Member Typedef Documentation

typedef std::pair<TrajectoryStateOnSurface,double> AnalyticalPropagator::TsosWP

private

Definition at line 119 of file AnalyticalPropagator.h.

Constructor & Destructor Documentation

AnalyticalPropagator::AnalyticalPropagator	(	const MagneticField *	field,
		PropagationDirection	dir = `alongMomentum`,
		float	maxDPhi = `1.6`,
		bool	isOld = `true`
	)

inline

Definition at line 27 of file AnalyticalPropagator.h.

Referenced by clone().

                                                  :
     Propagator(dir),
     theMaxDPhi2(maxDPhi*maxDPhi),
     theMaxDBzRatio(0.5),
     theField(field),
     isOldPropagationType(isOld) {}

AnalyticalPropagator::~AnalyticalPropagator ( )

inline

Definition at line 36 of file AnalyticalPropagator.h.

36 {}

Member Function Documentation

virtual AnalyticalPropagator* AnalyticalPropagator::clone ( void ) const

inlinevirtual

Implements Propagator.

Definition at line 70 of file AnalyticalPropagator.h.

References AnalyticalPropagator().

   {
     return new AnalyticalPropagator(*this);
   }

virtual const MagneticField* AnalyticalPropagator::magneticField ( ) const

inlineprivatevirtual

Implements Propagator.

Definition at line 116 of file AnalyticalPropagator.h.

References theField.

116 {return theField;}

AnalyticalPropagator::theField

const MagneticField * theField

Definition: AnalyticalPropagator.h:122

std::pair< TrajectoryStateOnSurface, double > AnalyticalPropagator::propagatedStateWithPath	(	const FreeTrajectoryState &	fts,
		const Surface &	surface,
		const GlobalTrajectoryParameters &	gtp,
		const double &	s
	)		const

private

propagation of errors (if needed) and generation of a new TSOS

Definition at line 109 of file AnalyticalPropagator.cc.

References SurfaceSideDefinition::afterSurface, alongMomentum, SurfaceSideDefinition::beforeSurface, FreeTrajectoryState::curvilinearError(), FreeTrajectoryState::hasError(), AnalyticalCurvilinearJacobian::jacobian(), CurvilinearTrajectoryError::matrix(), GlobalTrajectoryParameters::momentum(), FreeTrajectoryState::parameters(), GlobalTrajectoryParameters::position(), and alignCSCRings::s.

 {
   //
   // for forward propagation: state is before surface,
   // for backward propagation: state is after surface
   //
   SurfaceSide side = PropagationDirectionFromPath()(s,propagationDirection())==alongMomentum 
     ? beforeSurface : afterSurface;
   // 
   //
   // error propagation (if needed) and conversion to a TrajectoryStateOnSurface
   //
   if (fts.hasError()) {
     //
     // compute jacobian
     //
     AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
     const AlgebraicMatrix55 &jacobian = analyticalJacobian.jacobian();
     // CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()));
     return TsosWP(TrajectoryStateOnSurface(gtp,
                        ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()),
                        surface,side),s);
   }
   else {
     //
     // return state without errors
     //
     return TsosWP(TrajectoryStateOnSurface(gtp,surface,side),s);
   }
 }

bool AnalyticalPropagator::propagateParametersOnCylinder	(	const FreeTrajectoryState &	fts,
		const Cylinder &	cylinder,
		GlobalPoint &	x,
		GlobalVector &	p,
		double &	s
	)		const

private

parameter propagation to cylinder (returns position, momentum and path length)

Definition at line 143 of file AnalyticalPropagator.cc.

References funct::abs(), constexpr, alignCSCRings::e, f, PV3DBase< T, PVType, FrameType >::mag(), FreeTrajectoryState::momentum(), AlCaHLTBitMon_ParallelJobs::p, GloballyPositioned< T >::position(), FreeTrajectoryState::position(), alignCSCRings::s, FreeTrajectoryState::transverseCurvature(), Vector3DBase< T, FrameTag >::unit(), unlikely, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

 {
 
   GlobalPoint const & sp = cylinder.position();
   if unlikely(sp.x()!=0. || sp.y()!=0.) {
     throw PropagationException("Cannot propagate to an arbitrary cylinder");
   }
   // preset output
   x = fts.position();
   p = fts.momentum();
   s = 0;
   // (transverse) curvature
   auto rho = fts.transverseCurvature();
   //
   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
   // difference in transversal position at 10m.
   //
   if unlikely( std::abs(rho)<1.e-10f )
     return propagateWithLineCrossing(fts.position(),p,cylinder,x,s);
   //
   // Helix case
   //
   // check for possible intersection
   constexpr float tolerance = 1.e-4; // 1 micron distance
   auto rdiff = x.perp() - cylinder.radius();
   if ( std::abs(rdiff) < tolerance )  return true;
   //
   // Instantiate HelixBarrelCylinderCrossing and get solutions
   //
   HelixBarrelCylinderCrossing cylinderCrossing(fts.position(),fts.momentum(),rho,
                            propagationDirection(),cylinder);
   if unlikely( !cylinderCrossing.hasSolution() )  return false;
   // path length
   s = cylinderCrossing.pathLength();
   // point
   x = cylinderCrossing.position();
   // direction (renormalised)
   p = cylinderCrossing.direction().unit()*fts.momentum().mag();
   return true;
 }

bool AnalyticalPropagator::propagateParametersOnPlane	(	const FreeTrajectoryState &	fts,
		const Plane &	plane,
		GlobalPoint &	x,
		GlobalVector &	p,
		double &	s
	)		const

private

parameter propagation to plane (returns position, momentum and path length)

Definition at line 187 of file AnalyticalPropagator.cc.

References funct::abs(), alignCSCRings::e, f, likely, LogDebug, PV3DBase< T, PVType, FrameType >::mag(), Basic3DVector< T >::mag(), FreeTrajectoryState::momentum(), oppositeToMomentum, AlCaHLTBitMon_ParallelJobs::p, perp(), Geom::Phi< T >::phi(), PV3DBase< T, PVType, FrameType >::phi(), FreeTrajectoryState::position(), alignCSCRings::s, indexGen::s2, FreeTrajectoryState::transverseCurvature(), and unlikely.

 {
   // initialisation of position, momentum and path length
   x = fts.position();
   p = fts.momentum();
   s = 0;
   // (transverse) curvature
   auto rho = fts.transverseCurvature();
   //
   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
   // difference in transversal position at 10m.
   //
   if unlikely( std::abs(rho)<1.e-10f )
     return propagateWithLineCrossing(fts.position(),p,plane,x,s);
   //
   // Helix case 
   //
 
   //
   // Frame-independant point and vector are created explicitely to 
   // avoid confusing gcc (refuses to compile with temporary objects
   // in the constructor).
   //
   HelixPlaneCrossing::PositionType helixPos(x);
   HelixPlaneCrossing::DirectionType helixDir(p);
   if likely(isOldPropagationType) {
       OptimalHelixPlaneCrossing planeCrossing(plane,helixPos,helixDir,rho,propagationDirection());
       return propagateWithHelixCrossing(*planeCrossing,plane,fts.momentum().mag(),x,p,s);
     }
 
 
 
 
   //--- Alternative implementation to be used for the propagation of the parameters  of looping
   //    particles that cross twice the (infinite) surface of the plane. It is not trivial to determine
   //    which of the two intersections has to be returned.
   
   //---- FIXME: WHAT FOLLOWS HAS TO BE REWRITTEN IN A CLEANER (AND CPU-OPTIMIZED) WAY ---------
   LogDebug("AnalyticalPropagator") << "In AnaliticalProp, calling HAPC " << "\n"
                    << "plane is centered in xyz: " 
                    << plane.position().x() << " , "
                    << plane.position().y() << " , " 
                    << plane.position().z() << "\n";
   
   
   GlobalPoint gp1 = fts.position();
   GlobalVector gm1 = fts.momentum();
   double s1 = 0;
   double rho1 = fts.transverseCurvature();
   HelixPlaneCrossing::PositionType helixPos1(gp1);
   HelixPlaneCrossing::DirectionType helixDir1(gm1);
   LogDebug("AnalyticalPropagator") << "gp1 before calling planeCrossing1: " << gp1 << "\n";
   OptimalHelixPlaneCrossing planeCrossing1(plane,helixPos1,helixDir1,rho1,propagationDirection());
   
   HelixPlaneCrossing::PositionType xGen;
   HelixPlaneCrossing::DirectionType pGen;
       
   double tolerance(0.0050);
   if(propagationDirection()==oppositeToMomentum)
     tolerance *=-1;
   
   bool check1 = propagateWithHelixCrossing(*planeCrossing1,plane,fts.momentum().mag(),gp1,gm1,s1);
   double dphi1 = fabs(fts.momentum().phi()-gm1.phi());
   LogDebug("AnalyticalPropagator") << "check1, s1, dphi, gp1: " 
                    << check1 << " , "
                    << s1 << " , " 
                    << dphi1 << " , "
                    << gp1 << "\n";
   
   //move forward a bit to avoid that the propagator doesn't propagate because the state is already on surface.
   //we want to go to the other point of intersection between the helix and the plane
   xGen = (*planeCrossing1).position(s1+tolerance);
   pGen = (*planeCrossing1).direction(s1+tolerance);
   
   /*
     if(!check1 || s1>170 ){
     //PropagationDirection newDir = (propagationDirection() == alongMomentum) ? oppositeToMomentum : alongMomentum;
     PropagationDirection newDir = anyDirection;
     HelixArbitraryPlaneCrossing  planeCrossing1B(helixPos1,helixDir1,rho1,newDir);
     check1 = propagateWithHelixCrossing(planeCrossing1B,plane,fts.momentum().mag(),gp1,gm1,s1);
     LogDebug("AnalyticalPropagator") << "after second attempt, check1, s1,gp1: "
     << check1 << " , "
     << s1 << " , " << gp1 << "\n";
     
     xGen = planeCrossing1B.position(s1+tolerance);
     pGen = planeCrossing1B.direction(s1+tolerance);
     }
       */
   
   if(!check1){
     LogDebug("AnalyticalPropagator") << "failed also second attempt. No idea what to do, then bailout" << "\n";
   }
   
   
   pGen *= gm1.mag()/pGen.mag();
   GlobalPoint gp2(xGen);
   GlobalVector gm2(pGen);
   double s2 = 0;
   double rho2 = rho1;
   HelixPlaneCrossing::PositionType helixPos2(gp2);
   HelixPlaneCrossing::DirectionType helixDir2(gm2);
   OptimalHelixPlaneCrossing planeCrossing2(plane,helixPos2,helixDir2,rho2,propagationDirection());
   
   bool check2 = propagateWithHelixCrossing(*planeCrossing2,plane,gm2.mag(),gp2,gm2,s2);
   
   if(!check2){
     x = gp1;
     p = gm1;
     s = s1;
     return check1;
   }
   
   if(!check1){
     edm::LogError("AnalyticalPropagator") << "LOGIC ERROR: I should not have entered here!" << "\n";
     return false;
   }
   
   
   LogDebug("AnalyticalPropagator")  << "check2, s2, gp2: " 
                     << check2 << " , "
                     << s2 << " , " << gp2 << "\n";
   
   
   double dist1 = (plane.position()-gp1).perp();
   double dist2 = (plane.position()-gp2).perp();
   
   
   LogDebug("AnalyticalPropagator") << "propDir, dist1, dist2: " 
                    << propagationDirection() << " , "
                    << dist1 << " , " 
                    << dist2 << "\n";
   
   //If there are two solutions, the one which is the closest to the module's center is chosen
   if(dist1<2*dist2){
     x = gp1;
     p = gm1;
     s = s1;
     return check1;
   }else if(dist2<2*dist1){
     x = gp2;
     p = gm2;
     s = s1+s2+tolerance;
     return check2;
   }else{
     if(fabs(s1)<fabs(s2)){
       x = gp1;
       p = gm1;
       s = s1;
       return check1;
     }else{
       x = gp2;
       p = gm2;
       s = s1+s2+tolerance;
       return check2;
     }   
   }
   
   //-------- END of ugly piece of code  ---------------
   
 }

bool AnalyticalPropagator::propagateWithHelixCrossing	(	HelixPlaneCrossing &	planeCrossing,
		const Plane &	plane,
		const float	pmag,
		GlobalPoint &	x,
		GlobalVector &	p,
		double &	s
	)		const

private

helix parameter propagation to a plane using HelixPlaneCrossing

Definition at line 403 of file AnalyticalPropagator.cc.

References HelixPlaneCrossing::direction(), PV3DBase< T, PVType, FrameType >::mag(), HelixPlaneCrossing::pathLength(), HelixPlaneCrossing::position(), and unlikely.

                                            {
   // get solution
   std::pair<bool,double> propResult = planeCrossing.pathLength(plane);
   if unlikely( !propResult.first )  return false;
 
   s = propResult.second;
   x = GlobalPoint(planeCrossing.position(s));
   // direction (reconverted to GlobalVector, renormalised)
   GlobalVector pGen =  GlobalVector(planeCrossing.direction(s));
   pGen *= pmag/pGen.mag();
   p = pGen;
   //
   return true;
 }

bool AnalyticalPropagator::propagateWithLineCrossing	(	const GlobalPoint &	,
		const GlobalVector &	,
		const Plane &	,
		GlobalPoint &	,
		double &
	)		const

private

straight line parameter propagation to a plane

bool AnalyticalPropagator::propagateWithLineCrossing	(	const GlobalPoint &	,
		const GlobalVector &	,
		const Cylinder &	,
		GlobalPoint &	,
		double &
	)		const

private

straight line parameter propagation to a cylinder

std::pair<TrajectoryStateOnSurface,double> AnalyticalPropagator::propagateWithPath	(	const FreeTrajectoryState &	fts,
		const Plane &	plane
	)		const

overrideprivatevirtual

propagation to plane with path length

Implements Propagator.

Referenced by DualBzeroTrajectoryFactory::propagateExternal(), and DualTrajectoryFactory::propagateExternal().

std::pair<TrajectoryStateOnSurface,double> AnalyticalPropagator::propagateWithPath	(	const FreeTrajectoryState &	fts,
		const Cylinder &	cylinder
	)		const

overrideprivatevirtual

propagation to cylinder with path length

Implements Propagator.

virtual bool AnalyticalPropagator::setMaxDirectionChange ( float phiMax )

inlinevirtual

limitation of change in transverse direction (to avoid loops).

Reimplemented from Propagator.

Definition at line 65 of file AnalyticalPropagator.h.

References theMaxDPhi2.

                                                     { 
     theMaxDPhi2 = phiMax*phiMax;
     return true;
   }

void AnalyticalPropagator::setMaxRelativeChangeInBz ( const float maxDBz )

inline

Set the maximum relative change in Bz (Bz_at_end-Bz_at_start)/Bz_at_start for a single propagation. The default is no limit. NB: this propagator assumes constant, non-zero magnetic field parallel to the z-axis!

Definition at line 79 of file AnalyticalPropagator.h.

References theMaxDBzRatio.

                                                      {
     theMaxDBzRatio = maxDBz;
   }

Member Data Documentation

bool AnalyticalPropagator::isOldPropagationType

private

Definition at line 123 of file AnalyticalPropagator.h.

const MagneticField* AnalyticalPropagator::theField

private

Definition at line 122 of file AnalyticalPropagator.h.

Referenced by magneticField().

float AnalyticalPropagator::theMaxDBzRatio

private

Definition at line 121 of file AnalyticalPropagator.h.

Referenced by setMaxRelativeChangeInBz().

float AnalyticalPropagator::theMaxDPhi2

private

Definition at line 120 of file AnalyticalPropagator.h.

Referenced by setMaxDirectionChange().

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation