File indexing completed on 2024-04-06 12:31:44
0001 #include "TrackPropagation/RungeKutta/interface/RKPropagatorInZ.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003
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
0014
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
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
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
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
0053
0054
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); }