RungeKutta/src/RKPropagatorInS.cc

0001 #include "TrackPropagation/RungeKutta/interface/RKPropagatorInS.h"
0002 #include "RKCartesianDistance.h"
0003 #include "CartesianLorentzForce.h"
0004 #include "RKLocalFieldProvider.h"
0005 #include "DataFormats/GeometrySurface/interface/Plane.h"
0006 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0007 #include "RKAdaptiveSolver.h"
0008 #include "RKOne4OrderStep.h"
0009 #include "RKOneCashKarpStep.h"
0010 #include "PathToPlane2Order.h"
0011 #include "CartesianStateAdaptor.h"
0012 #include "TrackingTools/GeomPropagators/interface/StraightLineCylinderCrossing.h"
0013 #include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.h"
0014 #include "MagneticField/VolumeGeometry/interface/MagVolume.h"
0015 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionFromPath.h"
0016 #include "FrameChanger.h"
0017 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0018 #include "AnalyticalErrorPropagation.h"
0019 #include "GlobalParametersWithPath.h"
0020 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0021
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/Utilities/interface/Likely.h"
0024
0025 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInS::propagateWithPath(const FreeTrajectoryState& fts,
0026                                                                                const Plane& plane) const {
0027   GlobalParametersWithPath gp = propagateParametersOnPlane(fts, plane);
0028   if UNLIKELY (!gp)
0029     return TsosWP(TrajectoryStateOnSurface(), 0.);
0030
0031   SurfaceSideDefinition::SurfaceSide side =
0032       PropagationDirectionFromPath()(gp.s(), propagationDirection()) == alongMomentum
0033           ? SurfaceSideDefinition::beforeSurface
0034           : SurfaceSideDefinition::afterSurface;
0035   return analyticalErrorPropagation(fts, plane, side, gp.parameters(), gp.s());
0036 }
0037
0038 std::pair<TrajectoryStateOnSurface, double> RKPropagatorInS::propagateWithPath(const FreeTrajectoryState& fts,
0039                                                                                const Cylinder& cyl) const {
0040   GlobalParametersWithPath gp = propagateParametersOnCylinder(fts, cyl);
0041   if UNLIKELY (!gp)
0042     return TsosWP(TrajectoryStateOnSurface(), 0.);
0043
0044   SurfaceSideDefinition::SurfaceSide side =
0045       PropagationDirectionFromPath()(gp.s(), propagationDirection()) == alongMomentum
0046           ? SurfaceSideDefinition::beforeSurface
0047           : SurfaceSideDefinition::afterSurface;
0048   return analyticalErrorPropagation(fts, cyl, side, gp.parameters(), gp.s());
0049 }
0050
0051 GlobalParametersWithPath RKPropagatorInS::propagateParametersOnPlane(const FreeTrajectoryState& ts,
0052                                                                      const Plane& plane) const {
0053   GlobalPoint gpos(ts.position());
0054   GlobalVector gmom(ts.momentum());
0055   double startZ = plane.localZ(gpos);
0056   // (transverse) curvature
0057   double rho = ts.transverseCurvature();
0058   //
0059   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0060   // difference in transversal position at 10m.
0061   //
0062   if UNLIKELY (fabs(rho) < 1.e-10) {
0063     //
0064     // Instantiate auxiliary object for finding intersection.
0065     // Frame-independant point and vector are created explicitely to
0066     // avoid confusing gcc (refuses to compile with temporary objects
0067     // in the constructor).
0068     //
0069     LogDebug("RKPropagatorInS") << " startZ = " << startZ;
0070
0071     if UNLIKELY (fabs(startZ) < 1e-5) {
0072       LogDebug("RKPropagatorInS") << "Propagation is not performed: state is already on final surface.";
0073       GlobalTrajectoryParameters res(gpos, gmom, ts.charge(), theVolume);
0074       return GlobalParametersWithPath(res, 0.0);
0075     }
0076
0077     StraightLinePlaneCrossing::PositionType pos(gpos);
0078     StraightLinePlaneCrossing::DirectionType dir(gmom);
0079     StraightLinePlaneCrossing planeCrossing(pos, dir, propagationDirection());
0080     //
0081     // get solution
0082     //
0083     std::pair<bool, double> propResult = planeCrossing.pathLength(plane);
0084     if LIKELY (propResult.first && theVolume != nullptr) {
0085       double s = propResult.second;
0086       // point (reconverted to GlobalPoint)
0087       GlobalPoint x(planeCrossing.position(s));
0088       GlobalTrajectoryParameters res(x, gmom, ts.charge(), theVolume);
0089       return GlobalParametersWithPath(res, s);
0090     }
0091     //do someting
0092     LogDebug("RKPropagatorInS") << "Straight line propgation to plane failed !!";
0093     return GlobalParametersWithPath();
0094   }
0095
0096 #ifdef EDM_ML_DEBUG
0097   if (theVolume != 0) {
0098     LogDebug("RKPropagatorInS") << "RKPropagatorInS: starting prop to plane in volume with pos "
0099                                 << theVolume->position() << " Z axis " << theVolume->toGlobal(LocalVector(0, 0, 1));
0100
0101     LogDebug("RKPropagatorInS") << "The starting position is " << ts.position() << " (global) "
0102                                 << theVolume->toLocal(ts.position()) << " (local) ";
0103
0104     FrameChanger changer;
0105     auto localPlane = changer.transformPlane(plane, *theVolume);
0106     LogDebug("RKPropagatorInS") << "The plane position is " << plane.position() << " (global) " << localPlane.position()
0107                                 << " (local) ";
0108
0109     LogDebug("RKPropagatorInS") << "The initial distance to plane is " << plane.localZ(ts.position());
0110
0111     StraightLinePlaneCrossing cross(ts.position().basicVector(), ts.momentum().basicVector());
0112     std::pair<bool, double> res3 = cross.pathLength(plane);
0113     LogDebug("RKPropagatorInS") << "straight line distance " << res3.first << " " << res3.second;
0114   }
0115 #endif
0116
0117   typedef RKAdaptiveSolver<double, RKOneCashKarpStep, 6> Solver;
0118   typedef Solver::Vector RKVector;
0119
0120   RKLocalFieldProvider field(fieldProvider());
0121   PathToPlane2Order pathLength(field, &field.frame());
0122   CartesianLorentzForce deriv(field, ts.charge());
0123
0124   RKCartesianDistance dist;
0125   double eps = theTolerance;
0126   Solver solver;
0127   double stot = 0;
0128   PropagationDirection currentDirection = propagationDirection();
0129
0130   // in magVolume frame
0131   RKVector start(CartesianStateAdaptor::rkstate(rkPosition(gpos), rkMomentum(gmom)));
0132   int safeGuard = 0;
0133   while (safeGuard++ < 100) {
0134     CartesianStateAdaptor startState(start);
0135
0136     std::pair<bool, double> path =
0137         pathLength(plane, startState.position(), startState.momentum(), (double)ts.charge(), currentDirection);
0138     if UNLIKELY (!path.first) {
0139       LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to plane failed!"
0140                                   << "...distance to plane " << plane.localZ(globalPosition(startState.position()))
0141                                   << "...Local starting position in volume " << startState.position()
0142                                   << "...Magnetic field " << field.inTesla(startState.position());
0143
0144       return GlobalParametersWithPath();
0145     }
0146
0147     LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to plane is " << path.second;
0148
0149     double sstep = path.second;
0150     if UNLIKELY (std::abs(sstep) < eps) {
0151       LogDebug("RKPropagatorInS") << "On-surface accuracy not reached, but pathLength calculation says we are there! "
0152                                   << "path " << path.second << " distance to plane is " << startZ;
0153       GlobalTrajectoryParameters res(gtpFromVolumeLocal(startState, ts.charge()));
0154       return GlobalParametersWithPath(res, stot);
0155     }
0156
0157     LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep << " current distance to plane is "
0158                                 << startZ;
0159
0160     RKVector rkresult = solver(0, start, sstep, deriv, dist, eps);
0161     stot += sstep;
0162     CartesianStateAdaptor cur(rkresult);
0163     double remainingZ = plane.localZ(globalPosition(cur.position()));
0164
0165     if (fabs(remainingZ) < eps) {
0166       LogDebug("RKPropagatorInS") << "On-surface accuracy reached! " << remainingZ;
0167       GlobalTrajectoryParameters res(gtpFromVolumeLocal(cur, ts.charge()));
0168       return GlobalParametersWithPath(res, stot);
0169     }
0170
0171     start = rkresult;
0172
0173     if (remainingZ * startZ > 0) {
0174       LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again " << remainingZ;
0175     } else {
0176       LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction " << remainingZ;
0177       currentDirection = invertDirection(currentDirection);
0178     }
0179     startZ = remainingZ;
0180   }
0181
0182   edm::LogError("FailedPropagation") << " too many iterations trying to reach plane ";
0183   return GlobalParametersWithPath();
0184 }
0185
0186 GlobalParametersWithPath RKPropagatorInS::propagateParametersOnCylinder(const FreeTrajectoryState& ts,
0187                                                                         const Cylinder& cyl) const {
0188   typedef RKAdaptiveSolver<double, RKOneCashKarpStep, 6> Solver;
0189   typedef Solver::Vector RKVector;
0190
0191   const GlobalPoint& sp = cyl.position();
0192   if UNLIKELY (sp.x() != 0. || sp.y() != 0.) {
0193     throw PropagationException("Cannot propagate to an arbitrary cylinder");
0194   }
0195
0196   GlobalPoint gpos(ts.position());
0197   GlobalVector gmom(ts.momentum());
0198   LocalPoint pos(cyl.toLocal(gpos));
0199   LocalVector mom(cyl.toLocal(gmom));
0200   double startR = cyl.radius() - pos.perp();
0201
0202   // LogDebug("RKPropagatorInS")  << "RKPropagatorInS: starting from FTS " << ts ;
0203
0204   // (transverse) curvature
0205   double rho = ts.transverseCurvature();
0206   //
0207   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0208   // difference in transversal position at 10m.
0209   //
0210   if UNLIKELY (fabs(rho) < 1.e-10) {
0211     //
0212     // Instantiate auxiliary object for finding intersection.
0213     // Frame-independant point and vector are created explicitely to
0214     // avoid confusing gcc (refuses to compile with temporary objects
0215     // in the constructor).
0216     //
0217
0218     StraightLineBarrelCylinderCrossing cylCrossing(gpos, gmom, propagationDirection());
0219
0220     //
0221     // get solution
0222     //
0223     std::pair<bool, double> propResult = cylCrossing.pathLength(cyl);
0224     if LIKELY (propResult.first && theVolume != nullptr) {
0225       double s = propResult.second;
0226       // point (reconverted to GlobalPoint)
0227       GlobalPoint x(cylCrossing.position(s));
0228       GlobalTrajectoryParameters res(x, gmom, ts.charge(), theVolume);
0229       LogDebug("RKPropagatorInS") << "Straight line propagation to cylinder succeeded !!";
0230       return GlobalParametersWithPath(res, s);
0231     }
0232
0233     //do someting
0234     edm::LogError("RKPropagatorInS") << "Straight line propagation to cylinder failed !!";
0235     return GlobalParametersWithPath();
0236   }
0237
0238   RKLocalFieldProvider field(fieldProvider(cyl));
0239   // StraightLineCylinderCrossing pathLength( pos, mom, propagationDirection());
0240   CartesianLorentzForce deriv(field, ts.charge());
0241
0242   RKCartesianDistance dist;
0243   double eps = theTolerance;
0244   Solver solver;
0245   double stot = 0;
0246   PropagationDirection currentDirection = propagationDirection();
0247
0248   RKVector start(CartesianStateAdaptor::rkstate(pos.basicVector(), mom.basicVector()));
0249   int safeGuard = 0;
0250   while (safeGuard++ < 100) {
0251     CartesianStateAdaptor startState(start);
0252     StraightLineCylinderCrossing pathLength(
0253         LocalPoint(startState.position()), LocalVector(startState.momentum()), currentDirection, eps);
0254
0255     std::pair<bool, double> path = pathLength.pathLength(cyl);
0256     if UNLIKELY (!path.first) {
0257       LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path length calculation to cylinder failed!"
0258                                   << "Radius " << cyl.radius() << " pos.perp() "
0259                                   << LocalPoint(startState.position()).perp();
0260       return GlobalParametersWithPath();
0261     }
0262
0263     LogDebug("RKPropagatorInS") << "RKPropagatorInS: Path lenght to cylinder is " << path.second << " from point (R,z) "
0264                                 << startState.position().perp() << ", " << startState.position().z() << " to R "
0265                                 << cyl.radius();
0266
0267     double sstep = path.second;
0268     if UNLIKELY (std::abs(sstep) < eps) {
0269       LogDebug("RKPropagatorInS") << "accuracy not reached, but pathLength calculation says we are there! "
0270                                   << path.second;
0271
0272       GlobalTrajectoryParameters res(gtpFromLocal(startState.position(), startState.momentum(), ts.charge(), cyl));
0273       return GlobalParametersWithPath(res, stot);
0274     }
0275
0276     LogDebug("RKPropagatorInS") << "RKPropagatorInS: Solving for " << sstep << " current distance to cylinder is "
0277                                 << startR;
0278
0279     RKVector rkresult = solver(0, start, sstep, deriv, dist, eps);
0280     stot += sstep;
0281     CartesianStateAdaptor cur(rkresult);
0282     double remainingR = cyl.radius() - cur.position().perp();
0283
0284     if (fabs(remainingR) < eps) {
0285       LogDebug("RKPropagatorInS") << "Accuracy reached! " << remainingR;
0286       GlobalTrajectoryParameters res(gtpFromLocal(cur.position(), cur.momentum(), ts.charge(), cyl));
0287       return GlobalParametersWithPath(res, stot);
0288     }
0289
0290     start = rkresult;
0291     if (remainingR * startR > 0) {
0292       LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in same direction again " << remainingR;
0293     } else {
0294       LogDebug("RKPropagatorInS") << "Accuracy not reached yet, trying in opposite direction " << remainingR;
0295       currentDirection = invertDirection(currentDirection);
0296     }
0297     startR = remainingR;
0298   }
0299
0300   edm::LogError("FailedPropagation") << " too many iterations trying to reach cylinder ";
0301   return GlobalParametersWithPath();
0302 }
0303
0304 Propagator* RKPropagatorInS::clone() const { return new RKPropagatorInS(*this); }
0305
0306 GlobalTrajectoryParameters RKPropagatorInS::gtpFromLocal(const Basic3DVector<float>& lpos,
0307                                                          const Basic3DVector<float>& lmom,
0308                                                          TrackCharge ch,
0309                                                          const Surface& surf) const {
0310   return GlobalTrajectoryParameters(surf.toGlobal(LocalPoint(lpos)), surf.toGlobal(LocalVector(lmom)), ch, theVolume);
0311 }
0312
0313 RKLocalFieldProvider RKPropagatorInS::fieldProvider() const { return RKLocalFieldProvider(*theVolume); }
0314
0315 RKLocalFieldProvider RKPropagatorInS::fieldProvider(const Cylinder& cyl) const {
0316   return RKLocalFieldProvider(*theVolume, cyl);
0317 }
0318
0319 PropagationDirection RKPropagatorInS::invertDirection(PropagationDirection dir) const {
0320   if (dir == anyDirection)
0321     return dir;
0322   return (dir == alongMomentum ? oppositeToMomentum : alongMomentum);
0323 }
0324
0325 Basic3DVector<double> RKPropagatorInS::rkPosition(const GlobalPoint& pos) const {
0326   if (theVolume != nullptr)
0327     return theVolume->toLocal(pos).basicVector();
0328   else
0329     return pos.basicVector();
0330 }
0331
0332 Basic3DVector<double> RKPropagatorInS::rkMomentum(const GlobalVector& mom) const {
0333   if (theVolume != nullptr)
0334     return theVolume->toLocal(mom).basicVector();
0335   else
0336     return mom.basicVector();
0337 }
0338
0339 GlobalPoint RKPropagatorInS::globalPosition(const Basic3DVector<float>& pos) const {
0340   if (theVolume != nullptr)
0341     return theVolume->toGlobal(LocalPoint(pos));
0342   else
0343     return GlobalPoint(pos);
0344 }
0345
0346 GlobalVector RKPropagatorInS::globalMomentum(const Basic3DVector<float>& mom) const
0347
0348 {
0349   if (theVolume != nullptr)
0350     return theVolume->toGlobal(LocalVector(mom));
0351   else
0352     return GlobalVector(mom);
0353 }
0354
0355 GlobalTrajectoryParameters RKPropagatorInS::gtpFromVolumeLocal(const CartesianStateAdaptor& state,
0356                                                                TrackCharge charge) const {
0357   return GlobalTrajectoryParameters(
0358       globalPosition(state.position()), globalMomentum(state.momentum()), charge, theVolume);
0359 }