RungeKutta/src/RKPropagatorInZ.cc

0001 #include "TrackPropagation/RungeKutta/interface/RKPropagatorInZ.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 // #include "CommonReco/RKPropagators/interface/RK4PreciseSolver.h"
0004 #include "RKCurvilinearDistance.h"
0005 #include "CurvilinearLorentzForce.h"
0006 #include "RKLocalFieldProvider.h"
0007 #include "DataFormats/GeometrySurface/interface/Plane.h"
0008 #include "RKAdaptiveSolver.h"
0009 #include "RKOne4OrderStep.h"
0010 #include "RKOneCashKarpStep.h"
0011
0012 TrajectoryStateOnSurface RKPropagatorInZ::myPropagate(const FreeTrajectoryState& ts, const Plane& plane) const {
0013   //typedef RK4PreciseSolver<double,5>           Solver;
0014   //typedef RKAdaptiveSolver<double,RKOne4OrderStep, 5>   Solver;
0015   typedef RKAdaptiveSolver<double, RKOneCashKarpStep, 5> Solver;
0016   typedef Solver::Vector RKVector;
0017
0018   GlobalPoint pos(ts.position());
0019   GlobalVector mom(ts.momentum());
0020
0021   // cout << "RKPropagatorInZ: starting from FTS " << ts << endl;
0022
0023   LocalPoint startpos = plane.toLocal(pos);
0024   LocalVector startmom = plane.toLocal(mom);
0025   double pzSign = startmom.z() > 0 ? 1.0 : -1.0;
0026
0027   // cout << "In local plane coordinates: " << startpos << ", momentum " << startmom << endl;
0028
0029   RKVector start;
0030   start(0) = startpos.x();
0031   start(1) = startpos.y();
0032   start(2) = startmom.x() / startmom.z();
0033   start(3) = startmom.y() / startmom.z();
0034   start(4) = pzSign * ts.charge() / startmom.mag();
0035
0036   // cout << "RKPropagatorInZ: Solving with par " <<  startpos.z() << " and state " << start << endl;
0037
0038   RKLocalFieldProvider localField(*theVolume, plane);
0039
0040   CurvilinearLorentzForce<double, 5> deriv(localField);
0041   RKCurvilinearDistance<double, 5> dist;
0042   double eps = 1.e-5;
0043   Solver solver;
0044   try {
0045     RKVector rkresult = solver(startpos.z(), start, -startpos.z(), deriv, dist, eps);
0046
0047     return TrajectoryStateOnSurface(
0048         LocalTrajectoryParameters(rkresult(4), rkresult(2), rkresult(3), rkresult(0), rkresult(1), pzSign),
0049         plane,
0050         theVolume);
0051   } catch (CurvilinearLorentzForceException& e) {
0052     // the propagation failed due to momentum almost parallel to the plane.
0053     // This does not mean the propagation is impossible, but it should be done
0054     // in a different parametrization (e.g. s)
0055     return TrajectoryStateOnSurface();
0056   }
0057 }
0058
0059 TrajectoryStateOnSurface RKPropagatorInZ::myPropagate(const FreeTrajectoryState&, const Cylinder&) const {
0060   return TrajectoryStateOnSurface();
0061 }
0062
0063 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInZ::propagateWithPath(const FreeTrajectoryState&,
0064                                                                                const Plane&) const {
0065   return std::pair<TrajectoryStateOnSurface, double>();
0066 }
0067
0068 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInZ::propagateWithPath(const FreeTrajectoryState&,
0069                                                                                const Cylinder&) const {
0070   return std::pair<TrajectoryStateOnSurface, double>();
0071 }
0072
0073 Propagator* RKPropagatorInZ::clone() const { return new RKPropagatorInZ(*this); }