AnalyticalImpactPointExtrapolator.cc

Go to the documentation of this file.

#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
#include "TrackingTools/GeomPropagators/interface/IterativeHelixExtrapolatorToLine.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
#include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"

// #include "CommonDet/DetUtilities/interface/DetailedDetTimer.h"

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

AnalyticalImpactPointExtrapolator::AnalyticalImpactPointExtrapolator (const Propagator& propagator,
                                      const MagneticField* field) :
  thePropagator(propagator.clone()),
  theField(field)
{
  thePropagator->setPropagationDirection(anyDirection);
}

TrajectoryStateOnSurface 
AnalyticalImpactPointExtrapolator::extrapolate (const FreeTrajectoryState& fts, 
                        const GlobalPoint& vtx) const
{
//   static TimingReport::Item& timer = detailedDetTimer("AnalyticalImpactPointExtrapolator");
//   TimeMe t(timer,false);

  return extrapolateSingleState(fts, vtx);
}

TrajectoryStateOnSurface 
AnalyticalImpactPointExtrapolator::extrapolate (const TrajectoryStateOnSurface tsos, 
                        const GlobalPoint& vtx) const
{
  if ( tsos.isValid() )  return extrapolateFullState(tsos,vtx);
  else  return tsos;
}

TrajectoryStateOnSurface 
AnalyticalImpactPointExtrapolator::extrapolateFullState (const TrajectoryStateOnSurface tsos, 
                             const GlobalPoint& vertex) const
{
  //
  // first determine IP plane using propagation with (single) FTS
  // could be optimised (will propagate errors even if duplicated below)
  //
  TrajectoryStateOnSurface singleState = 
    extrapolateSingleState(*tsos.freeTrajectoryState(),vertex);
  if ( !singleState.isValid() || tsos.singleState() )  return singleState;
  //
  // propagate multiTsos to plane found above
  //
  return thePropagator->propagate(tsos,singleState.surface());
}

TrajectoryStateOnSurface 
AnalyticalImpactPointExtrapolator::extrapolateSingleState (const FreeTrajectoryState& fts, 
                               const GlobalPoint& vertex) const
{
  //
  // initialisation of position, momentum and transverse curvature
  //
  GlobalPoint x(fts.position());
  GlobalVector p(fts.momentum());
  double rho = fts.transverseCurvature();
  //
  // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um 
  // difference in transversal position at 10m.
  //
  double s(0);
  if( fabs(rho)<1.e-10 ) {
    GlobalVector dx(p.dot(vertex-x)/p.mag2()*p);
    x += dx;
    float sign = p.dot(dx);
    s = sign>0 ? dx.mag() : -dx.mag();
  }
  //
  // Helix case 
  //
  else {
    HelixLineExtrapolation::PositionType helixPos(x);
    HelixLineExtrapolation::DirectionType helixDir(p);
    IterativeHelixExtrapolatorToLine extrapolator(helixPos,helixDir,rho,anyDirection);
    if ( !propagateWithHelix(extrapolator,vertex,x,p,s) )  return TrajectoryStateOnSurface();
  }
  //
  // Define target surface: origin on line, x_local from line 
  //   to helix at closest approach, z_local along the helix
  //   and y_local to complete right-handed system
  //
  GlobalVector zLocal(p.unit());
  GlobalVector yLocal(zLocal.cross(x-vertex).unit());
  GlobalVector xLocal(yLocal.cross(zLocal));
  Surface::RotationType rot(xLocal,yLocal,zLocal);
  PlaneBuilder::ReturnType surface = PlaneBuilder().plane(vertex,rot);
  //
  // Compute propagated state
  //
  GlobalTrajectoryParameters gtp(x,p,fts.charge(), theField);
  if (fts.hasError()) {
    //
    // compute jacobian
    //
    AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
    CurvilinearTrajectoryError cte( ROOT::Math::Similarity( analyticalJacobian.jacobian(), 
                                                            fts.curvilinearError().matrix()) );
    return TrajectoryStateOnSurface(gtp,cte,*surface);
  }
  else {
    //
    // return state without errors
    //
    return TrajectoryStateOnSurface(gtp,*surface);
  }
}

bool
AnalyticalImpactPointExtrapolator::propagateWithHelix (const IterativeHelixExtrapolatorToLine& extrapolator,
                               const GlobalPoint& vertex,
                               GlobalPoint& x, GlobalVector& p, double& s) const {
  //
  // save absolute value of momentum
  //
  double pmag(p.mag());
  //
  // get path length to solution
  //
  std::pair<bool,double> propResult = extrapolator.pathLength(vertex);
  if ( !propResult.first )  return false;
  s = propResult.second;
  // 
  // get point and (normalised) direction from path length
  //
  HelixLineExtrapolation::PositionType xGen = extrapolator.position(s);
  HelixLineExtrapolation::DirectionType pGen = extrapolator.direction(s);
  //
  // Fix normalisation and convert back to GlobalPoint / GlobalVector
  //
  x = GlobalPoint(xGen);
  pGen *= pmag/pGen.mag();
  p = GlobalVector(pGen);
  //
  return true;
}