RungeKutta/src/RKPropagatorInR.cc

0001 #include "TrackPropagation/RungeKutta/interface/RKPropagatorInR.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 // #include "CommonReco/RKPropagators/interface/RK4PreciseSolver.h"
0004 #include "RKCylindricalDistance.h"
0005 #include "CylindricalLorentzForce.h"
0006 #include "RKLocalFieldProvider.h"
0007 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0008 #include "RKAdaptiveSolver.h"
0009 #include "RKOne4OrderStep.h"
0010 #include "RKOneCashKarpStep.h"
0011 #include "CylindricalState.h"
0012 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
0013
0014 TrajectoryStateOnSurface RKPropagatorInR::myPropagate(const FreeTrajectoryState& ts, const Cylinder& cyl) const {
0015   //typedef RK4PreciseSolver<double,5>           Solver;
0016   typedef RKAdaptiveSolver<double, RKOne4OrderStep, 5> Solver;
0017   //typedef RKAdaptiveSolver<Scalar,RKOneCashKarpStep, 5>   Solver;
0018   typedef double Scalar;
0019   typedef Solver::Vector RKVector;
0020
0021   GlobalPoint pos(ts.position());
0022   GlobalVector mom(ts.momentum());
0023
0024   LocalPoint startpos = cyl.toLocal(pos);
0025   LocalVector startmom = cyl.toLocal(mom);
0026
0027   CylindricalState startState(startpos, startmom, ts.charge());
0028   const RKVector& start = startState.parameters();
0029
0030   RKLocalFieldProvider localField(*theVolume, cyl);
0031   CylindricalLorentzForce<double, 5> deriv(localField);
0032   RKCylindricalDistance<double, 5> dist;
0033   double eps = 1.e-5;
0034   Solver solver;
0035   try {
0036     Scalar step = cyl.radius() - startState.rho();
0037     RKVector rkresult = solver(startState.rho(), start, step, deriv, dist, eps);
0038     CylindricalState endState(cyl.radius(), rkresult, startState.prSign());
0039     return TrajectoryStateOnSurface(GlobalTrajectoryParameters(cyl.toGlobal(endState.position()),
0040                                                                cyl.toGlobal(endState.momentum()),
0041                                                                TrackCharge(endState.charge()),
0042                                                                theVolume),
0043                                     cyl);
0044   } catch (CylindricalLorentzForceException& e) {
0045     // the propagation failed due to momentum almost parallel to the plane.
0046     // This does not mean the propagation is impossible, but it should be done
0047     // in a different parametrization (e.g. s)
0048     return TrajectoryStateOnSurface();
0049   }
0050 }
0051
0052 TrajectoryStateOnSurface RKPropagatorInR::myPropagate(const FreeTrajectoryState&, const Plane&) const {
0053   return TrajectoryStateOnSurface();
0054 }
0055
0056 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInR::propagateWithPath(const FreeTrajectoryState&,
0057                                                                                const Plane&) const {
0058   return std::pair<TrajectoryStateOnSurface, double>();
0059 }
0060
0061 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInR::propagateWithPath(const FreeTrajectoryState&,
0062                                                                                const Cylinder&) const {
0063   return std::pair<TrajectoryStateOnSurface, double>();
0064 }
0065
0066 Propagator* RKPropagatorInR::clone() const { return new RKPropagatorInR(*this); }