RKPropagatorInZ.cc

Go to the documentation of this file.

#include "TrackPropagation/RungeKutta/interface/RKPropagatorInZ.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
// #include "CommonReco/RKPropagators/interface/RK4PreciseSolver.h"
#include "RKCurvilinearDistance.h"
#include "CurvilinearLorentzForce.h"
#include "RKLocalFieldProvider.h"
#include "DataFormats/GeometrySurface/interface/Plane.h"
#include "RKAdaptiveSolver.h"
#include "RKOne4OrderStep.h"
#include "RKOneCashKarpStep.h"

TrajectoryStateOnSurface RKPropagatorInZ::myPropagate(const FreeTrajectoryState& ts, const Plane& plane) const {
  //typedef RK4PreciseSolver<double,5>           Solver;
  //typedef RKAdaptiveSolver<double,RKOne4OrderStep, 5>   Solver;
  typedef RKAdaptiveSolver<double, RKOneCashKarpStep, 5> Solver;
  typedef Solver::Vector RKVector;

  GlobalPoint pos(ts.position());
  GlobalVector mom(ts.momentum());

  // cout << "RKPropagatorInZ: starting from FTS " << ts << endl;

  LocalPoint startpos = plane.toLocal(pos);
  LocalVector startmom = plane.toLocal(mom);
  double pzSign = startmom.z() > 0 ? 1.0 : -1.0;

  // cout << "In local plane coordinates: " << startpos << ", momentum " << startmom << endl;

  RKVector start;
  start(0) = startpos.x();
  start(1) = startpos.y();
  start(2) = startmom.x() / startmom.z();
  start(3) = startmom.y() / startmom.z();
  start(4) = pzSign * ts.charge() / startmom.mag();

  // cout << "RKPropagatorInZ: Solving with par " <<  startpos.z() << " and state " << start << endl;

  RKLocalFieldProvider localField(*theVolume, plane);

  CurvilinearLorentzForce<double, 5> deriv(localField);
  RKCurvilinearDistance<double, 5> dist;
  double eps = 1.e-5;
  Solver solver;
  try {
    RKVector rkresult = solver(startpos.z(), start, -startpos.z(), deriv, dist, eps);

    return TrajectoryStateOnSurface(
        LocalTrajectoryParameters(rkresult(4), rkresult(2), rkresult(3), rkresult(0), rkresult(1), pzSign),
        plane,
        theVolume);
  } catch (CurvilinearLorentzForceException& e) {
    // the propagation failed due to momentum almost parallel to the plane.
    // This does not mean the propagation is impossible, but it should be done
    // in a different parametrization (e.g. s)
    return TrajectoryStateOnSurface();
  }
}

TrajectoryStateOnSurface RKPropagatorInZ::myPropagate(const FreeTrajectoryState&, const Cylinder&) const {
  return TrajectoryStateOnSurface();
}

std::pair<TrajectoryStateOnSurface, double> RKPropagatorInZ::propagateWithPath(const FreeTrajectoryState&,
                                                                               const Plane&) const {
  return std::pair<TrajectoryStateOnSurface, double>();
}

std::pair<TrajectoryStateOnSurface, double> RKPropagatorInZ::propagateWithPath(const FreeTrajectoryState&,
                                                                               const Cylinder&) const {
  return std::pair<TrajectoryStateOnSurface, double>();
}

Propagator* RKPropagatorInZ::clone() const { return new RKPropagatorInZ(*this); }