TransverseImpactPointExtrapolator.cc

Go to the documentation of this file.

#include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "DataFormats/GeometrySurface/interface/Surface.h"
#include "boost/intrusive_ptr.hpp"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "DataFormats/GeometryVector/interface/Point2DBase.h"
#include "DataFormats/GeometryVector/interface/Vector2DBase.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"


TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator () :
  thePropagator(nullptr) {}


TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator (const MagneticField* field) :
  thePropagator(new AnalyticalPropagator(field, anyDirection)) {}

TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator (const Propagator& u) :
  thePropagator(u.clone())
{
  thePropagator->setPropagationDirection(anyDirection);
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::extrapolate (const FreeTrajectoryState& fts,
                        const GlobalPoint& vtx) const
{
  return doExtrapolation(fts, vtx, *thePropagator);
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::extrapolate (const TrajectoryStateOnSurface tsos,
                        const GlobalPoint& vtx) const
{
  if ( !tsos.isValid() )  return tsos;
  return doExtrapolation(tsos, vtx, *thePropagator);
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::extrapolate (const FreeTrajectoryState& fts,
                        const GlobalPoint& vtx,
                        const Propagator& p) const
{
  //
  // set propagator for bidirectional propagation
  //
   std::unique_ptr<Propagator> p_cloned = SetPropagationDirection(p, anyDirection);
   return doExtrapolation(fts, vtx, *(p_cloned.get()));
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::extrapolate (const TrajectoryStateOnSurface tsos,
                        const GlobalPoint& vtx,
                        const Propagator& p) const
{
  if ( !tsos.isValid() )  return tsos;
  //
  // set propagator for bidirectional propagation
  //
  std::unique_ptr<Propagator> p_cloned = SetPropagationDirection(p, anyDirection);
  return doExtrapolation(tsos, vtx, *(p_cloned.get()));
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::doExtrapolation (const TrajectoryStateOnSurface tsos,
                            const GlobalPoint& vtx,
                            const Propagator& p) const
{
  //
  // Compute tip surface
  //
  if (fabs(tsos.freeState()->transverseCurvature())<1.E-6){
    LogDebug("TransverseImpactPointExtrapolator")<< "negligeable curvature: using a trick to extrapolate:\n"<<tsos;

    //0T field probably
    //x is perpendicular to the momentum
    GlobalVector xLocal = GlobalVector(-tsos.globalMomentum().y(),tsos.globalMomentum().x(),0).unit();
    //y along global Z
    GlobalVector yLocal(0.,0.,1.);
    //z accordingly
    GlobalVector zLocal(xLocal.cross(yLocal));

    const Surface::PositionType& origin(vtx);
    Surface::RotationType rotation(xLocal,yLocal,zLocal);
    ReferenceCountingPointer<Plane> surface =  PlaneBuilder().plane(origin,rotation);

    return p.propagate(*tsos.freeState(),*surface);
  }else{
  ReferenceCountingPointer<Plane> surface =
    tipSurface(tsos.globalPosition(),tsos.globalMomentum(),
           1./tsos.transverseCurvature(),vtx);
  //
  // propagate
  //
  return p.propagate(tsos,*surface);
  }
}

TrajectoryStateOnSurface
TransverseImpactPointExtrapolator::doExtrapolation (const FreeTrajectoryState& fts,
                            const GlobalPoint& vtx,
                            const Propagator& p) const
{
  //
  // Compute tip surface
  //
  if (fabs(fts.transverseCurvature())<1.E-6){
    LogDebug("TransverseImpactPointExtrapolator")<< "negligeable curvature: using a trick to extrapolate:\n"<<fts;

    //0T field probably
    //x is perpendicular to the momentum
    GlobalVector xLocal = GlobalVector(-fts.momentum().y(),fts.momentum().x(),0).unit();
    //y along global Z
    GlobalVector yLocal(0.,0.,1.);
    //z accordingly
    GlobalVector zLocal(xLocal.cross(yLocal));

    const Surface::PositionType& origin(vtx);
    Surface::RotationType rotation(xLocal,yLocal,zLocal);
    ReferenceCountingPointer<Plane> surface =  PlaneBuilder().plane(origin,rotation);

    return p.propagate(fts,*surface);
  }else{
  ReferenceCountingPointer<Plane> surface =
    tipSurface(fts.position(),fts.momentum(),
           1./fts.transverseCurvature(),vtx);
  //
  // propagate
  //
  return p.propagate(fts,*surface);
  }
}

ReferenceCountingPointer<Plane>
TransverseImpactPointExtrapolator::tipSurface (const GlobalPoint& position,
                           const GlobalVector& momentum,
                           const double& signedTransverseRadius,
                           const GlobalPoint& vertex) const
{

  LogDebug("TransverseImpactPointExtrapolator")<< position<<"\n"
                            <<momentum<<"\n"
                            <<"signedTransverseRadius : "<<signedTransverseRadius<<"\n"
                            <<vertex;

  typedef Point2DBase<double,GlobalTag> PositionType2D;
  typedef Vector2DBase<double,GlobalTag> DirectionType2D;

  PositionType2D x0(position.x(),position.y());
  DirectionType2D t0(-momentum.y(),momentum.x());
  t0 = t0/t0.mag();

  PositionType2D xc(x0+signedTransverseRadius*t0);

  DirectionType2D vtxDirection(xc.x()-vertex.x(),xc.y()-vertex.y());
  double vtxDistance = vtxDirection.mag();

  const Surface::PositionType& origin(vertex);
  GlobalVector xLocal(vtxDirection.x()/vtxDistance,
              vtxDirection.y()/vtxDistance,
              0.);
  if ( vtxDistance<fabs(signedTransverseRadius) ) {
    LogDebug("TransverseImpactPointExtrapolator")<<"Inverting the x axis.";
    xLocal = -xLocal;
  }
  GlobalVector yLocal(0.,0.,1.);
  GlobalVector zLocal(xLocal.cross(yLocal));
  if ( zLocal.dot(momentum)<0. ) {
    LogDebug("TransverseImpactPointExtrapolator")<<"Inverting the y,z frame.";
    yLocal = -yLocal;
    zLocal = -zLocal;
  }
  Surface::RotationType rotation(xLocal,yLocal,zLocal);

  LogDebug("TransverseImpactPointExtrapolator")<<"plane center: "<<origin<<"\n"
                           <<"plane rotation axis:\n"
                           <<xLocal<<"\n"
                           <<yLocal<<"\n"
                           <<zLocal<<"\n"
                           <<"x0: "<<x0<<"\n"
                           <<"t0: "<<t0<<"\n"
                           <<"xc: "<<xc<<"\n"
                           <<"vtxDirection: "<<vtxDirection;

  return PlaneBuilder().plane(origin,rotation);
}