#include <ImpactPoint3dEstimator.h>

Inheritance diagram for Trk::ImpactPoint3dEstimator:

Collaboration diagram for Trk::ImpactPoint3dEstimator:

Public Member Functions
virtual StatusCode	initialize () override

virtual StatusCode	finalize () override

	ImpactPoint3dEstimator (const std::string &t, const std::string &n, const IInterface *p)
	Default constructor due to Athena interface. More...

virtual	~ImpactPoint3dEstimator ()

virtual std::unique_ptr< PlaneSurface >	Estimate3dIP (const Trk::TrackParameters trackPerigee, const Amg::Vector3D theVertex, double &distance) const override

virtual std::unique_ptr< PlaneSurface >	Estimate3dIP (const Trk::NeutralParameters neutralPerigee, const Amg::Vector3D theVertex, double &distance) const override

virtual bool	addIP3dAtaPlane (VxTrackAtVertex &, const Amg::Vector3D &vertex) const override
	Actual estimate method, changing the state of Trk::VxTrackAtVertex. More...

virtual const Trk::AtaPlane *	IP3dAtaPlane (VxTrackAtVertex &vtxTrack, const Amg::Vector3D &vertex) const override
	This method creates the ImpactPoint3dAtaPlane as the parameters of the track at the point of closest approach in 3d to the given vertex. More...

virtual const Trk::NeutralAtaPlane *	IP3dNeutralAtaPlane (const NeutralParameters *initNeutPerigee, const Amg::Vector3D &vertex) const override

Private Member Functions
template<typename T >
std::unique_ptr< PlaneSurface >	Estimate3dIPNoCurvature (const T , const Amg::Vector3D theVertex, double &distance) const

Private Attributes
ToolHandle< Trk::IExtrapolator >	m_extrapolator

SG::ReadCondHandleKey< AtlasFieldCacheCondObj >	m_fieldCacheCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}

int	m_maxiterations

double	m_precision

Detailed Description

This object calculates the point of minimum distance to the vertex in 3d. Consider that this point is different from the usually used point of closest approach on the transverse plane.

This algorithm makes use of a simple iterative Newton process

Author: N. Giacinto Piacquadio (for the Freiburg Group)

Changes:

David Shope david.nosp@m..ric.nosp@m.hard..nosp@m.shop.nosp@m.e@cer.nosp@m.n.ch (2016-03-18)

EDM Migration to xAOD - move Trk::Vertex to Amg::Vector3D

Definition at line 39 of file ImpactPoint3dEstimator.h.

Constructor & Destructor Documentation

◆ ImpactPoint3dEstimator()

Trk::ImpactPoint3dEstimator::ImpactPoint3dEstimator	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p
	)

Default constructor due to Athena interface.

Definition at line 28 of file ImpactPoint3dEstimator.cxx.

   {
     declareProperty("Extrapolator",m_extrapolator);
     declareProperty("MaxIterations",m_maxiterations);
     declareProperty("Precision",m_precision);
   }
  
   ImpactPoint3dEstimator::~ImpactPoint3dEstimator() = default;

◆ ~ImpactPoint3dEstimator()

Trk::ImpactPoint3dEstimator::~ImpactPoint3dEstimator ( )

virtualdefault

Member Function Documentation

◆ addIP3dAtaPlane()

bool Trk::ImpactPoint3dEstimator::addIP3dAtaPlane	(	VxTrackAtVertex &	vtxTrack,
		const Amg::Vector3D &	vertex
	)		const

overridevirtual

Actual estimate method, changing the state of Trk::VxTrackAtVertex.

Definition at line 287 of file ImpactPoint3dEstimator.cxx.

                                    {
       const AtaPlane* myPlane=IP3dAtaPlane(vtxTrack,vertex);
       if (myPlane)
         {
           vtxTrack.setImpactPoint3dAtaPlane(myPlane);
           return true;
         }
     } else { //for neutrals
       const NeutralAtaPlane* myPlane=IP3dNeutralAtaPlane(vtxTrack.initialNeutralPerigee(),vertex);
       if (myPlane)      {
         ATH_MSG_VERBOSE ("Adding plane: " << myPlane->associatedSurface() );
         vtxTrack.setImpactPoint3dNeutralAtaPlane(myPlane);
         return true;
       }
     }
     return false;
   }
  
  

◆ Estimate3dIP() [1/2]

std::unique_ptr< PlaneSurface > Trk::ImpactPoint3dEstimator::Estimate3dIP	(	const Trk::NeutralParameters *	neutralPerigee,
		const Amg::Vector3D *	theVertex,
		double &	distance
	)		const

overridevirtual

Definition at line 108 of file ImpactPoint3dEstimator.cxx.

   {
     ATH_MSG_DEBUG("Neutral particle --  propagate like a straight line");
     return Estimate3dIPNoCurvature(neutralPerigee, theVertex, distance);
   }
  
   std::unique_ptr<PlaneSurface>

◆ Estimate3dIP() [2/2]

std::unique_ptr< PlaneSurface > Trk::ImpactPoint3dEstimator::Estimate3dIP	(	const Trk::TrackParameters *	trackPerigee,
		const Amg::Vector3D *	theVertex,
		double &	distance
	)		const

overridevirtual

Definition at line 117 of file ImpactPoint3dEstimator.cxx.

   {
     SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, Gaudi::Hive::currentContext()};
     const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
  
     MagField::AtlasFieldCache fieldCache;
     fieldCondObj->getInitializedCache (fieldCache);
  
     double magnFieldVect[3];
     fieldCache.getField(trackPerigee->associatedSurface().center().data(),magnFieldVect);
     if(magnFieldVect[2] == 0 ){
       ATH_MSG_DEBUG("Magnetic field in the Z direction is 0 --  propagate like a straight line");
       return Estimate3dIPNoCurvature(trackPerigee, theVertex, distance);
     }
     
     const Trk::Perigee* thePerigee=dynamic_cast<const Trk::Perigee*>(trackPerigee);
     if (thePerigee==nullptr)
     {
       ATH_MSG_DEBUG(
         " ImpactPoint3dEstimator didn't get a Perigee* as ParametersBase*: "
         "cast not possible. Need to EXTRAPOLATE...");
       Trk::PerigeeSurface perigeeSurface(*theVertex);
       std::unique_ptr<const Trk::TrackParameters> tmp =
         m_extrapolator->extrapolateDirectly(
           Gaudi::Hive::currentContext(), *trackPerigee, perigeeSurface);
       if (tmp && tmp->associatedSurface().type() == Trk::SurfaceType::Perigee) {
         thePerigee = static_cast<const Trk::Perigee*>(tmp.release());
       }
       if (thePerigee == nullptr){
         return nullptr;
       }
     }
  
     ATH_MSG_VERBOSE( " Now running ImpactPoint3dEstimator::Estimate3dIP" );
     double phi0=thePerigee->parameters()[Trk::phi0];
     double dCosPhi0=-std::sin(phi0);
     double dSinPhi0=std::cos(phi0);
     double theta=thePerigee->parameters()[Trk::theta];
     double cotTheta=1./std::tan(thePerigee->parameters()[Trk::theta]);
     double d0=thePerigee->parameters()[Trk::d0];
  
     //I need the radius (magnetic field...)
     double Bz=magnFieldVect[2]*299.792;
     double Rt=std::sin(theta)/(Bz*thePerigee->parameters()[Trk::qOverP]);
  
     double x0=thePerigee->associatedSurface().center().x()+(d0-Rt)*dCosPhi0;
     double y0=thePerigee->associatedSurface().center().y()+(d0-Rt)*dSinPhi0;
     double z0=thePerigee->associatedSurface().center().z()+thePerigee->parameters()[Trk::z0]+Rt*phi0*cotTheta;
  
     if (thePerigee!=trackPerigee) {
       delete thePerigee;
       thePerigee=nullptr;
     }
  
     double xc=theVertex->x();
     double yc=theVertex->y();
     double zc=theVertex->z();
  
     double phiActual=phi0;
     double dCosPhiActual=-std::sin(phiActual);
     double dSinPhiActual=std::cos(phiActual);
  
     double secderivative=0.;
     double derivative=0.;
  
     int ncycle=0;
     bool isok=false;
  
     double deltaphi=0.;
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
     std::cout << std::setprecision(25) << "actual distance  before cycle is: " << std::hypot(x0-xc+Rt*dCosPhiActual, 
                                                                                   y0-yc+Rt*dSinPhiActual, 
                                                                                   z0-zc-Rt*cotTheta*phiActual) << std::endl;
 #endif
  
     do {
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
       ATH_MSG_VERBOSE( "Cycle number: " << ncycle << " old phi: " << phiActual  );
 #endif
      
  
       derivative=(x0-xc)*(-Rt*dSinPhiActual)+(y0-yc)*Rt*dCosPhiActual+(z0-zc-Rt*phiActual*cotTheta)*(-Rt*cotTheta);
       secderivative=Rt*(-(x0-xc)*dCosPhiActual-(y0-yc)*dSinPhiActual+Rt*cotTheta*cotTheta);
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
       ATH_MSG_VERBOSE( "derivative is: " << derivative << " sec derivative is: " << secderivative  );
 #endif
  
       deltaphi=-derivative/secderivative;
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
       std::cout << std::setprecision(25) << "deltaphi: " << deltaphi << std::endl;
 #endif
  
       phiActual+=deltaphi;
       dCosPhiActual=-std::sin(phiActual);
       dSinPhiActual=std::cos(phiActual);
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
       ATH_MSG_VERBOSE( "derivative is: " << derivative << " sec derivative is: " << secderivative  );
       std::cout << std::setprecision(25) << std::hypot(x0-xc+Rt*dCosPhiActual, y0-yc+Rt*dSinPhiActual,
                                                       z0-zc-Rt*cotTheta*phiActual) << std::endl;
       ATH_MSG_VERBOSE( "actual distance is: " << std::hypot(x0-xc+Rt*dCosPhiActual,
                                                             y0-yc+Rt*dSinPhiActual,
                                                             z0-zc-Rt*cotTheta*phiActual));
 #endif
  
       if (secderivative<0) throw error::ImpactPoint3dEstimatorProblem("Second derivative is negative");
  
       if (ncycle>m_maxiterations) throw error::ImpactPoint3dEstimatorProblem("Too many loops: could not find minimum distance to vertex");
  
       ncycle+=1;
       if (ncycle>m_maxiterations||std::abs(deltaphi)<m_precision) {
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
         ATH_MSG_VERBOSE( "found minimum at: " << phiActual  );
 #endif
         isok=true;
       }
  
     } while (!isok);
  
     //now you have to construct the plane with PlaneSurface
     //first vector at 3d impact point
     Amg::Vector3D MomentumDir(std::cos(phiActual)*std::sin(theta),std::sin(phiActual)*std::sin(theta),std::cos(theta));
     Amg::Vector3D DeltaR(x0-xc+Rt*dCosPhiActual,y0-yc+Rt*dSinPhiActual,z0-zc-Rt*cotTheta*phiActual);
     distance=DeltaR.mag();
     if (distance==0.){
       ATH_MSG_WARNING("DeltaR is zero in ImpactPoint3dEstimator::Estimate3dIP, returning nullptr");
       return nullptr;
     }
     DeltaR=DeltaR.unit();
  
  
     //correct DeltaR from small deviations from orthogonality to DeltaR
  
     Amg::Vector3D DeltaRcorrected=DeltaR-(DeltaR.dot(MomentumDir))*MomentumDir;
  
     if ((DeltaR-DeltaRcorrected).mag()>1e-4)
     {
       ATH_MSG_WARNING( " DeltaR and MomentumDir are not orthogonal "  );
       ATH_MSG_DEBUG( std::setprecision(10) << " DeltaR-DeltaRcorrected: "  << (DeltaR-DeltaRcorrected).mag()  );
     }
  
     Amg::Vector3D YDir=MomentumDir.cross(DeltaRcorrected);
  
     //store the impact 3d point
     Amg::Vector3D vertex(x0+Rt*dCosPhiActual,y0+Rt*dSinPhiActual,z0-Rt*cotTheta*phiActual);
  
     ATH_MSG_VERBOSE( "final minimal distance is: " << std::hypot(x0-xc+Rt*dCosPhiActual,
                                                                  y0-yc+Rt*dSinPhiActual,
                                                                  z0-zc-Rt*cotTheta*phiActual));
  
     ATH_MSG_DEBUG( "POCA in 3D is: " << vertex  );
  
  
     //store the plane...
     ATH_MSG_VERBOSE( "plane to which to extrapolate X " << DeltaRcorrected << " Y " << YDir << " Z " << MomentumDir  );
  
     Amg::Transform3D thePlane(DeltaRcorrected, YDir, MomentumDir, *theVertex);
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
     std::cout << "the translation is, directly from Transform3d: " << thePlane.getTranslation() << endmsg;
 #endif
     return std::make_unique<PlaneSurface>(thePlane);
   }//end of estimate 3dIP method
  
   bool

◆ Estimate3dIPNoCurvature()

template<typename T >

std::unique_ptr< PlaneSurface > Trk::ImpactPoint3dEstimator::Estimate3dIPNoCurvature	(	const T *	thePerigee,
		const Amg::Vector3D *	theVertex,
		double &	distance
	)		const

private

Definition at line 63 of file ImpactPoint3dEstimator.cxx.

   {
       const Amg::Vector3D  momentumUnit = thePerigee->momentum().unit();
       double pathLength  =  ( *theVertex  - thePerigee->position() ).dot( momentumUnit )
                                                                        / (  momentumUnit.dot( momentumUnit )) ;
        //first vector at 3d impact point
  
  
       Amg::Vector3D POCA   =  thePerigee->position()  + pathLength * momentumUnit;// Position of closest approach
       Amg::Vector3D DeltaR =  *theVertex  - POCA;
       distance=DeltaR.mag();
       DeltaR=DeltaR.unit();
  
  
       //correct DeltaR from small deviations from orthogonality to DeltaR -- DeltaR.dot(momentumUnit) should equal 0 if the above is correct
  
       Amg::Vector3D DeltaRcorrected=DeltaR-(DeltaR.dot(momentumUnit))*momentumUnit;
  
       if ((DeltaR-DeltaRcorrected).mag()>1e-4)
       {
         ATH_MSG_WARNING( " DeltaR and MomentumDir are not orthogonal " );
         ATH_MSG_DEBUG( std::setprecision(10) << " DeltaR-DeltaRcorrected: "  << (DeltaR-DeltaRcorrected).mag() );
       }
  
       Amg::Vector3D YDir=momentumUnit.cross(DeltaRcorrected);
  
       ATH_MSG_VERBOSE( "final minimal distance is: " << distance);
       ATH_MSG_DEBUG( "POCA in 3D is: " << POCA );
  
       //store the plane...
       ATH_MSG_VERBOSE( "plane to which to extrapolate X " << DeltaRcorrected << " Y " << YDir << " Z " << momentumUnit);
  
       Amg::Transform3D thePlane(DeltaRcorrected, YDir, momentumUnit, *theVertex);
  
 #ifdef IMPACTPOINT3DESTIMATOR_DEBUG
       std::cout << "the translation is, directly from Transform3d: " << thePlane->getTranslation() << endmsg;
 #endif
  
       return std::make_unique<PlaneSurface>(thePlane);
  
   }
  
   std::unique_ptr<PlaneSurface>

◆ finalize()

StatusCode Trk::ImpactPoint3dEstimator::finalize ( )

overridevirtual

Definition at line 54 of file ImpactPoint3dEstimator.cxx.

◆ initialize()

StatusCode Trk::ImpactPoint3dEstimator::initialize ( )

overridevirtual

Definition at line 41 of file ImpactPoint3dEstimator.cxx.

                                   {
       ATH_CHECK( m_extrapolator.retrieve() );
     } else {
       m_extrapolator.disable();
     }
     ATH_CHECK( m_fieldCacheCondObjInputKey.initialize() );
  
     ATH_MSG_DEBUG( "Initialize successful"  );
     return StatusCode::SUCCESS;
   }
  
   StatusCode ImpactPoint3dEstimator::finalize()

◆ IP3dAtaPlane()

const Trk::AtaPlane * Trk::ImpactPoint3dEstimator::IP3dAtaPlane	(	VxTrackAtVertex &	vtxTrack,
		const Amg::Vector3D &	vertex
	)		const

overridevirtual

This method creates the ImpactPoint3dAtaPlane as the parameters of the track at the point of closest approach in 3d to the given vertex.

The parameters and errors are defined on the plane intersecting the track at point of closest approach, with track ortogonal to the plane and center of the plane defined as the given vertex.

Definition at line 309 of file ImpactPoint3dEstimator.cxx.

     {
       double distance = 0;
       theSurfaceAtIP = Estimate3dIP(vtxTrack.initialPerigee(),&vertex,distance);
     }
     catch (error::ImpactPoint3dEstimatorProblem err)
     {
       ATH_MSG_WARNING( " ImpactPoint3dEstimator failed to find minimum distance between track and vertex seed: " << err.p  );
       return nullptr;
     }
     if(!theSurfaceAtIP){
       ATH_MSG_WARNING( " ImpactPoint3dEstimator failed to find minimum distance and returned 0 " );
       return nullptr;
     }
 #ifdef ImpactPoint3dAtaPlaneFactory_DEBUG
     ATH_MSG_VERBOSE( "Original perigee was: " << *(vtxTrack.initialPerigee())  );
     ATH_MSG_VERBOSE( "The resulting surface is: " << *theSurfaceAtIP  );
 #endif
    const auto* pTrackPar = m_extrapolator->extrapolate(
      Gaudi::Hive::currentContext(),
      *(vtxTrack.initialPerigee()),
      *theSurfaceAtIP).release();
    if (const Trk::AtaPlane* res = dynamic_cast<const Trk::AtaPlane *>(pTrackPar); res){
      return res;
    }
    ATH_MSG_WARNING("TrackParameters ptr returned from extrapolate could not be cast to Trk::AtaPlane* in IP3dAtaPlane(..)");
    return nullptr;
   }
  
  

◆ IP3dNeutralAtaPlane()

const Trk::NeutralAtaPlane * Trk::ImpactPoint3dEstimator::IP3dNeutralAtaPlane	(	const NeutralParameters *	initNeutPerigee,
		const Amg::Vector3D &	vertex
	)		const

overridevirtual

Definition at line 345 of file ImpactPoint3dEstimator.cxx.

     {
         double distance = 0;
         theSurfaceAtIP = Estimate3dIP(initNeutPerigee,&vertex,distance);
     }
     catch (error::ImpactPoint3dEstimatorProblem err)
     {
       ATH_MSG_WARNING( " ImpactPoint3dEstimator failed to find minimum distance between track and vertex seed: " << err.p  );
       return nullptr;
     }
     if(!theSurfaceAtIP){
       ATH_MSG_WARNING( " ImpactPoint3dEstimator failed to find minimum distance and returned 0 " );
       return nullptr;
     } 
 #ifdef ImpactPoint3dAtaPlaneFactory_DEBUG
     ATH_MSG_VERBOSE( "Original neutral perigee was: " << *initNeutPerigee  );
     ATH_MSG_VERBOSE( "The resulting surface is: " << *theSurfaceAtIP  );
 #endif
  
     const Trk::NeutralAtaPlane* res = nullptr;
     std::unique_ptr<const Trk::NeutralParameters> tmp =  m_extrapolator->extrapolate(*initNeutPerigee,*theSurfaceAtIP);
     if(dynamic_cast<const Trk::NeutralAtaPlane*> (tmp.get())){
       res = static_cast<const Trk::NeutralAtaPlane*> (tmp.release());
     }
    return res;
   }
  
  

Member Data Documentation

◆ m_extrapolator

ToolHandle< Trk::IExtrapolator > Trk::ImpactPoint3dEstimator::m_extrapolator

private

Definition at line 88 of file ImpactPoint3dEstimator.h.

◆ m_fieldCacheCondObjInputKey

SG::ReadCondHandleKey<AtlasFieldCacheCondObj> Trk::ImpactPoint3dEstimator::m_fieldCacheCondObjInputKey {this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"}

private

Definition at line 89 of file ImpactPoint3dEstimator.h.

◆ m_maxiterations

int Trk::ImpactPoint3dEstimator::m_maxiterations

private

Definition at line 92 of file ImpactPoint3dEstimator.h.

◆ m_precision

double Trk::ImpactPoint3dEstimator::m_precision

private

Definition at line 93 of file ImpactPoint3dEstimator.h.

The documentation for this class was generated from the following files:

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ ImpactPoint3dEstimator()

◆ ~ImpactPoint3dEstimator()

Member Function Documentation

◆ addIP3dAtaPlane()

◆ Estimate3dIP() [1/2]

◆ Estimate3dIP() [2/2]

◆ Estimate3dIPNoCurvature()

◆ finalize()

◆ initialize()

◆ IP3dAtaPlane()

◆ IP3dNeutralAtaPlane()

Member Data Documentation

◆ m_extrapolator

◆ m_fieldCacheCondObjInputKey

◆ m_maxiterations

◆ m_precision