Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:31:56

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); }