Trk::ImpactPoint3dEstimator Class Reference

This object calculates the point of minimum distance to the vertex in 3d. More...

#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.
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.
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.
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 26 of file ImpactPoint3dEstimator.cxx.

                                                                                                             :
    base_class(t,n,p),
    m_extrapolator(""),
    m_maxiterations(20),
    m_precision(1e-10)//DeltaPhi
  {
    declareProperty("Extrapolator",m_extrapolator);
    declareProperty("MaxIterations",m_maxiterations);
    declareProperty("Precision",m_precision);
  }

~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 285 of file ImpactPoint3dEstimator.cxx.

  {
    if (vtxTrack.initialPerigee()) {
      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 106 of file ImpactPoint3dEstimator.cxx.

  {
    ATH_MSG_DEBUG("Neutral particle --  propagate like a straight line");
    return Estimate3dIPNoCurvature(neutralPerigee, theVertex, distance);
  }

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 115 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

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 61 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);
 
  }

finalize()

StatusCode Trk::ImpactPoint3dEstimator::finalize ( )

overridevirtual

Definition at line 52 of file ImpactPoint3dEstimator.cxx.

  {
    ATH_MSG_DEBUG( "Finalize successful"  );
    return StatusCode::SUCCESS;
  }

initialize()

StatusCode Trk::ImpactPoint3dEstimator::initialize ( )

overridevirtual

Definition at line 39 of file ImpactPoint3dEstimator.cxx.

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

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 307 of file ImpactPoint3dEstimator.cxx.

  {
    if (!vtxTrack.initialPerigee() && vtxTrack.initialNeutralPerigee())
      ATH_MSG_WARNING( "Calling ImpactPoint3dEstimator::IP3dAtaPlane cannot return NeutralAtaPlane"  );
    std::unique_ptr<PlaneSurface> theSurfaceAtIP;
    try
    {
      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_DEBUG("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 343 of file ImpactPoint3dEstimator.cxx.

  {
    std::unique_ptr<PlaneSurface> theSurfaceAtIP;
 
    try
    {
        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.

{this, "AtlasFieldCacheCondObj", "fieldCondObj", "Name of the Magnetic Field conditions object key"};

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:

• ImpactPoint3dEstimator.h
• ImpactPoint3dEstimator.cxx