File indexing completed on 2024-10-04 22:55:08
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
0057 double rho = ts.transverseCurvature();
0058
0059
0060
0061
0062 if UNLIKELY (fabs(rho) < 1.e-10) {
0063
0064
0065
0066
0067
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
0082
0083 std::pair<bool, double> propResult = planeCrossing.pathLength(plane);
0084 if LIKELY (propResult.first && theVolume != nullptr) {
0085 double s = propResult.second;
0086
0087 GlobalPoint x(planeCrossing.position(s));
0088 GlobalTrajectoryParameters res(x, gmom, ts.charge(), theVolume);
0089 return GlobalParametersWithPath(res, s);
0090 }
0091
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
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
0203
0204
0205 double rho = ts.transverseCurvature();
0206
0207
0208
0209
0210 if UNLIKELY (fabs(rho) < 1.e-10) {
0211
0212
0213
0214
0215
0216
0217
0218 StraightLineBarrelCylinderCrossing cylCrossing(gpos, gmom, propagationDirection());
0219
0220
0221
0222
0223 std::pair<bool, double> propResult = cylCrossing.pathLength(cyl);
0224 if LIKELY (propResult.first && theVolume != nullptr) {
0225 double s = propResult.second;
0226
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
0234 edm::LogError("RKPropagatorInS") << "Straight line propagation to cylinder failed !!";
0235 return GlobalParametersWithPath();
0236 }
0237
0238 RKLocalFieldProvider field(fieldProvider(cyl));
0239
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 {
0314 assert(theVolume);
0315 return RKLocalFieldProvider(*theVolume);
0316 }
0317
0318 RKLocalFieldProvider RKPropagatorInS::fieldProvider(const Cylinder& cyl) const {
0319 return RKLocalFieldProvider(*theVolume, cyl);
0320 }
0321
0322 PropagationDirection RKPropagatorInS::invertDirection(PropagationDirection dir) const {
0323 if (dir == anyDirection)
0324 return dir;
0325 return (dir == alongMomentum ? oppositeToMomentum : alongMomentum);
0326 }
0327
0328 Basic3DVector<double> RKPropagatorInS::rkPosition(const GlobalPoint& pos) const {
0329 if (theVolume != nullptr)
0330 return theVolume->toLocal(pos).basicVector();
0331 else
0332 return pos.basicVector();
0333 }
0334
0335 Basic3DVector<double> RKPropagatorInS::rkMomentum(const GlobalVector& mom) const {
0336 if (theVolume != nullptr)
0337 return theVolume->toLocal(mom).basicVector();
0338 else
0339 return mom.basicVector();
0340 }
0341
0342 GlobalPoint RKPropagatorInS::globalPosition(const Basic3DVector<float>& pos) const {
0343 if (theVolume != nullptr)
0344 return theVolume->toGlobal(LocalPoint(pos));
0345 else
0346 return GlobalPoint(pos);
0347 }
0348
0349 GlobalVector RKPropagatorInS::globalMomentum(const Basic3DVector<float>& mom) const
0350
0351 {
0352 if (theVolume != nullptr)
0353 return theVolume->toGlobal(LocalVector(mom));
0354 else
0355 return GlobalVector(mom);
0356 }
0357
0358 GlobalTrajectoryParameters RKPropagatorInS::gtpFromVolumeLocal(const CartesianStateAdaptor& state,
0359 TrackCharge charge) const {
0360 return GlobalTrajectoryParameters(
0361 globalPosition(state.position()), globalMomentum(state.momentum()), charge, theVolume);
0362 }