AnalyticalTrajectoryExtrapolatorToLine.cc

Go to the documentation of this file.

#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.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/PlaneBuilder.h"
#include "DataFormats/GeometrySurface/interface/Line.h"

#include <cmath>

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

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

TrajectoryStateOnSurface AnalyticalTrajectoryExtrapolatorToLine::extrapolate(const FreeTrajectoryState& fts,
                                                                             const Line& line) const {
  return extrapolateSingleState(fts, line);
}

TrajectoryStateOnSurface AnalyticalTrajectoryExtrapolatorToLine::extrapolate(const TrajectoryStateOnSurface tsos,
                                                                             const Line& line) const {
  if (tsos.isValid())
    return extrapolateFullState(tsos, line);
  else
    return tsos;
}

TrajectoryStateOnSurface AnalyticalTrajectoryExtrapolatorToLine::extrapolateFullState(
    const TrajectoryStateOnSurface tsos, const Line& line) 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(), line);
  if (!singleState.isValid() || tsos.singleState())
    return singleState;
  //
  // propagate multiTsos to plane found above
  //
  return thePropagator->propagate(tsos, singleState.surface());
}

TrajectoryStateOnSurface AnalyticalTrajectoryExtrapolatorToLine::extrapolateSingleState(const FreeTrajectoryState& fts,
                                                                                        const Line& line) const {
  //   static TimingReport::Item& timer = detailedDetTimer("AnalyticalTrajectoryExtrapolatorToLine");
  //   TimeMe t(timer,false);
  //
  // 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) {
    Line tangent(x, p);
    GlobalPoint xold(x);
    x = tangent.closerPointToLine(line);
    GlobalVector dx(x - xold);
    float sign = p.dot(x - xold);
    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, line, 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
  //
  GlobalPoint origin(line.closerPointToLine(Line(x, p)));
  GlobalVector zLocal(p.unit());
  GlobalVector yLocal(zLocal.cross(x - origin).unit());
  GlobalVector xLocal(yLocal.cross(zLocal));
  Surface::RotationType rot(xLocal, yLocal, zLocal);
  PlaneBuilder::ReturnType surface = PlaneBuilder().plane(origin, rot);
  //
  // Compute propagated state
  //
  GlobalTrajectoryParameters gtp(x, p, fts.charge(), thePropagator->magneticField());
  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 TrajectoryStateOnSurface(gtp, cte, *surface);
  } else {
    //
    // return state without errors
    //
    return TrajectoryStateOnSurface(gtp, *surface);
  }
}

bool AnalyticalTrajectoryExtrapolatorToLine::propagateWithHelix(const IterativeHelixExtrapolatorToLine& extrapolator,
                                                                const Line& line,
                                                                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(line);
  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;
}