CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/TrackPropagation/RungeKutta/src/RKPropagatorInR.cc

Go to the documentation of this file.
00001 #include "TrackPropagation/RungeKutta/interface/RKPropagatorInR.h"
00002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00003 // #include "CommonReco/RKPropagators/interface/RK4PreciseSolver.h"
00004 #include "RKCylindricalDistance.h"
00005 #include "CylindricalLorentzForce.h"
00006 #include "RKLocalFieldProvider.h"
00007 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
00008 #include "RKAdaptiveSolver.h"
00009 #include "RKOne4OrderStep.h"
00010 #include "RKOneCashKarpStep.h"
00011 #include "CylindricalState.h"
00012 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
00013 
00014 TrajectoryStateOnSurface 
00015 RKPropagatorInR::propagate (const FreeTrajectoryState& ts, const Cylinder& cyl) const
00016 {
00017   //typedef RK4PreciseSolver<double,5>           Solver;
00018     typedef RKAdaptiveSolver<double,RKOne4OrderStep, 5>     Solver;
00019     //typedef RKAdaptiveSolver<Scalar,RKOneCashKarpStep, 5>   Solver;
00020     typedef double                                          Scalar;
00021     typedef Solver::Vector                                  RKVector;
00022 
00023     GlobalPoint pos( ts.position());
00024     GlobalVector mom( ts.momentum());
00025 
00026     LocalPoint startpos = cyl.toLocal(pos);
00027     LocalVector startmom = cyl.toLocal(mom);
00028 
00029     CylindricalState startState( startpos, startmom, ts.charge());
00030     RKVector start = startState.parameters();
00031 
00032     RKLocalFieldProvider localField( *theVolume, cyl);
00033     CylindricalLorentzForce<double,5> deriv(localField);
00034     RKCylindricalDistance<double,5> dist;
00035     double eps = 1.e-5;
00036     Solver solver;
00037     try {
00038         Scalar step = cyl.radius() - startState.rho();
00039         RKVector rkresult = solver( startState.rho(), start, step, deriv, dist, eps);
00040         CylindricalState endState( cyl.radius(), rkresult, startState.prSign());
00041         return TrajectoryStateOnSurface( GlobalTrajectoryParameters( cyl.toGlobal( endState.position()), 
00042                                                                      cyl.toGlobal( endState.momentum()),
00043                                                                      TrackCharge(endState.charge()), 
00044                                                                      theVolume),
00045                                          cyl);
00046     }
00047     catch (CylindricalLorentzForceException& e) {
00048         // the propagation failed due to momentum almost parallel to the plane.
00049         // This does not mean the propagation is impossible, but it should be done
00050         // in a different parametrization (e.g. s)  
00051         return TrajectoryStateOnSurface();
00052     }
00053 }
00054 
00055 TrajectoryStateOnSurface 
00056 RKPropagatorInR::propagate (const FreeTrajectoryState&, const Plane&) const
00057 {
00058     return TrajectoryStateOnSurface();
00059 }
00060 
00061 std::pair< TrajectoryStateOnSurface, double> 
00062 RKPropagatorInR::propagateWithPath (const FreeTrajectoryState&, const Plane&) const
00063 {
00064     return std::pair< TrajectoryStateOnSurface, double>();
00065 }
00066 
00067 std::pair< TrajectoryStateOnSurface, double> 
00068 RKPropagatorInR::propagateWithPath (const FreeTrajectoryState&, const Cylinder&) const
00069 {
00070     return std::pair< TrajectoryStateOnSurface, double>();
00071 }
00072 
00073 Propagator * RKPropagatorInR::clone() const
00074 {
00075     return new RKPropagatorInR(*this);
00076 }