GeomPropagators/src/AnalyticalPropagator.cc

0001 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0004 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0005 #include "DataFormats/GeometrySurface/interface/Plane.h"
0006
0007 #include "MagneticField/Engine/interface/MagneticField.h"
0008
0009 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0010 #include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
0011 #include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.h"
0012 #include "TrackingTools/GeomPropagators/interface/OptimalHelixPlaneCrossing.h"
0013 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
0014 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0015 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionFromPath.h"
0016 #include "TrackingTools/TrajectoryState/interface/SurfaceSideDefinition.h"
0017 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/Utilities/interface/Likely.h"
0021
0022 #include <cmath>
0023
0024 using namespace SurfaceSideDefinition;
0025
0026 std::pair<TrajectoryStateOnSurface, double> AnalyticalPropagator::propagateWithPath(const FreeTrajectoryState& fts,
0027                                                                                     const Plane& plane) const {
0028   // check curvature
0029   float rho = fts.transverseCurvature();
0030
0031   // propagate parameters
0032   GlobalPoint x;
0033   GlobalVector p;
0034   double s;
0035
0036   // check if already on plane
0037   if LIKELY (plane.localZclamped(fts.position()) != 0) {
0038     // propagate
0039     bool parametersOK = this->propagateParametersOnPlane(fts, plane, x, p, s);
0040     // check status and deltaPhi limit
0041     float dphi2 = float(s) * rho;
0042     dphi2 = dphi2 * dphi2 * fts.momentum().perp2();
0043     if UNLIKELY (!parametersOK || dphi2 > theMaxDPhi2 * fts.momentum().mag2())
0044       return TsosWP(TrajectoryStateOnSurface(), 0.);
0045   } else {
0046     LogDebug("AnalyticalPropagator") << "not going anywhere. Already on surface.\n"
0047                                      << "plane.localZ(fts.position()): " << plane.localZ(fts.position()) << "\n"
0048                                      << "plane.position().mag(): " << plane.position().mag() << "\n"
0049                                      << "plane.posPrec: " << plane.posPrec();
0050     x = fts.position();
0051     p = fts.momentum();
0052     s = 0;
0053   }
0054   //
0055   // Compute propagated state and check change in curvature
0056   //
0057   GlobalTrajectoryParameters gtp(x, p, fts.charge(), theField);
0058   if UNLIKELY (std::abs(gtp.transverseCurvature() - rho) > theMaxDBzRatio * std::abs(rho))
0059     return TsosWP(TrajectoryStateOnSurface(), 0.);
0060   //
0061   // construct TrajectoryStateOnSurface
0062   //
0063   return propagatedStateWithPath(fts, plane, gtp, s);
0064 }
0065
0066 std::pair<TrajectoryStateOnSurface, double> AnalyticalPropagator::propagateWithPath(const FreeTrajectoryState& fts,
0067                                                                                     const Cylinder& cylinder) const {
0068   // check curvature
0069   auto rho = fts.transverseCurvature();
0070
0071   // propagate parameters
0072   GlobalPoint x;
0073   GlobalVector p;
0074   double s = 0;
0075
0076   bool parametersOK = this->propagateParametersOnCylinder(fts, cylinder, x, p, s);
0077   // check status and deltaPhi limit
0078   float dphi2 = s * rho;
0079   dphi2 = dphi2 * dphi2 * fts.momentum().perp2();
0080   if UNLIKELY (!parametersOK || dphi2 > theMaxDPhi2 * fts.momentum().mag2())
0081     return TsosWP(TrajectoryStateOnSurface(), 0.);
0082   //
0083   // Compute propagated state and check change in curvature
0084   //
0085   GlobalTrajectoryParameters gtp(x, p, fts.charge(), theField);
0086   if UNLIKELY (std::abs(gtp.transverseCurvature() - rho) > theMaxDBzRatio * std::abs(rho))
0087     return TsosWP(TrajectoryStateOnSurface(), 0.);
0088   //
0089   // create result TSOS on TangentPlane (local parameters & errors are better defined)
0090   //
0091
0092   //try {
0093   ConstReferenceCountingPointer<TangentPlane> plane(
0094       cylinder.tangentPlane(x));  // need to be here until tsos is created!
0095   return propagatedStateWithPath(fts, *plane, gtp, s);
0096   /*
0097   } catch(...) {
0098     std::cout << "wrong tangent to cylinder " << x
0099               << " pos, rad " << cylinder.position() << " " << cylinder.radius()
0100               << std::endl;
0101     return TsosWP(TrajectoryStateOnSurface(),0.);
0102   }
0103   */
0104 }
0105
0106 std::pair<TrajectoryStateOnSurface, double> AnalyticalPropagator::propagatedStateWithPath(
0107     const FreeTrajectoryState& fts,
0108     const Surface& surface,
0109     const GlobalTrajectoryParameters& gtp,
0110     const double& s) const {
0111   //
0112   // for forward propagation: state is before surface,
0113   // for backward propagation: state is after surface
0114   //
0115   SurfaceSide side =
0116       PropagationDirectionFromPath()(s, propagationDirection()) == alongMomentum ? beforeSurface : afterSurface;
0117   //
0118   //
0119   // error propagation (if needed) and conversion to a TrajectoryStateOnSurface
0120   //
0121   if (fts.hasError()) {
0122     //
0123     // compute jacobian
0124     //
0125     AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
0126     const AlgebraicMatrix55& jacobian = analyticalJacobian.jacobian();
0127     // CurvilinearTrajectoryError cte(ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()));
0128     return TsosWP(
0129         TrajectoryStateOnSurface(gtp, ROOT::Math::Similarity(jacobian, fts.curvilinearError().matrix()), surface, side),
0130         s);
0131   } else {
0132     //
0133     // return state without errors
0134     //
0135     return TsosWP(TrajectoryStateOnSurface(gtp, surface, side), s);
0136   }
0137 }
0138
0139 bool AnalyticalPropagator::propagateParametersOnCylinder(
0140     const FreeTrajectoryState& fts, const Cylinder& cylinder, GlobalPoint& x, GlobalVector& p, double& s) const {
0141   GlobalPoint const& sp = cylinder.position();
0142   if UNLIKELY (sp.x() != 0. || sp.y() != 0.) {
0143     throw PropagationException("Cannot propagate to an arbitrary cylinder");
0144   }
0145   // preset output
0146   x = fts.position();
0147   p = fts.momentum();
0148   s = 0;
0149   // (transverse) curvature
0150   auto rho = fts.transverseCurvature();
0151   //
0152   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0153   // difference in transversal position at 10m.
0154   //
0155   if UNLIKELY (std::abs(rho) < 1.e-10f)
0156     return propagateWithLineCrossing(fts.position(), p, cylinder, x, s);
0157   //
0158   // Helix case
0159   //
0160   // check for possible intersection
0161   constexpr float tolerance = 1.e-4;  // 1 micron distance
0162   auto rdiff = x.perp() - cylinder.radius();
0163   if (std::abs(rdiff) < tolerance)
0164     return true;
0165   //
0166   // Instantiate HelixBarrelCylinderCrossing and get solutions
0167   //
0168   HelixBarrelCylinderCrossing cylinderCrossing(fts.position(), fts.momentum(), rho, propagationDirection(), cylinder);
0169   if UNLIKELY (!cylinderCrossing.hasSolution())
0170     return false;
0171   // path length
0172   s = cylinderCrossing.pathLength();
0173   // point
0174   x = cylinderCrossing.position();
0175   // direction (renormalised)
0176   p = cylinderCrossing.direction().unit() * fts.momentum().mag();
0177   return true;
0178 }
0179
0180 bool AnalyticalPropagator::propagateParametersOnPlane(
0181     const FreeTrajectoryState& fts, const Plane& plane, GlobalPoint& x, GlobalVector& p, double& s) const {
0182   // initialisation of position, momentum and path length
0183   x = fts.position();
0184   p = fts.momentum();
0185   s = 0;
0186   // (transverse) curvature
0187   auto rho = fts.transverseCurvature();
0188   //
0189   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0190   // difference in transversal position at 10m.
0191   //
0192   if UNLIKELY (std::abs(rho) < 1.e-10f)
0193     return propagateWithLineCrossing(fts.position(), p, plane, x, s);
0194   //
0195   // Helix case
0196   //
0197
0198   //
0199   // Frame-independant point and vector are created explicitely to
0200   // avoid confusing gcc (refuses to compile with temporary objects
0201   // in the constructor).
0202   //
0203   HelixPlaneCrossing::PositionType helixPos(x);
0204   HelixPlaneCrossing::DirectionType helixDir(p);
0205   if LIKELY (isOldPropagationType) {
0206     OptimalHelixPlaneCrossing planeCrossing(plane, helixPos, helixDir, rho, propagationDirection());
0207     return propagateWithHelixCrossing(*planeCrossing, plane, fts.momentum().mag(), x, p, s);
0208   }
0209
0210   //--- Alternative implementation to be used for the propagation of the parameters  of looping
0211   //    particles that cross twice the (infinite) surface of the plane. It is not trivial to determine
0212   //    which of the two intersections has to be returned.
0213
0214   //---- FIXME: WHAT FOLLOWS HAS TO BE REWRITTEN IN A CLEANER (AND CPU-OPTIMIZED) WAY ---------
0215   LogDebug("AnalyticalPropagator") << "In AnaliticalProp, calling HAPC "
0216                                    << "\n"
0217                                    << "plane is centered in xyz: " << plane.position().x() << " , "
0218                                    << plane.position().y() << " , " << plane.position().z() << "\n";
0219
0220   GlobalPoint gp1 = fts.position();
0221   GlobalVector gm1 = fts.momentum();
0222   double s1 = 0;
0223   double rho1 = fts.transverseCurvature();
0224   HelixPlaneCrossing::PositionType helixPos1(gp1);
0225   HelixPlaneCrossing::DirectionType helixDir1(gm1);
0226   LogDebug("AnalyticalPropagator") << "gp1 before calling planeCrossing1: " << gp1 << "\n";
0227   OptimalHelixPlaneCrossing planeCrossing1(plane, helixPos1, helixDir1, rho1, propagationDirection());
0228
0229   HelixPlaneCrossing::PositionType xGen;
0230   HelixPlaneCrossing::DirectionType pGen;
0231
0232   double tolerance(0.0050);
0233   if (propagationDirection() == oppositeToMomentum)
0234     tolerance *= -1;
0235
0236   bool check1 = propagateWithHelixCrossing(*planeCrossing1, plane, fts.momentum().mag(), gp1, gm1, s1);
0237   double dphi1 = fabs(fts.momentum().phi() - gm1.phi());
0238   LogDebug("AnalyticalPropagator") << "check1, s1, dphi, gp1: " << check1 << " , " << s1 << " , " << dphi1 << " , "
0239                                    << gp1 << "\n";
0240
0241   //move forward a bit to avoid that the propagator doesn't propagate because the state is already on surface.
0242   //we want to go to the other point of intersection between the helix and the plane
0243   xGen = (*planeCrossing1).position(s1 + tolerance);
0244   pGen = (*planeCrossing1).direction(s1 + tolerance);
0245
0246   /*
0247     if(!check1 || s1>170 ){
0248     //PropagationDirection newDir = (propagationDirection() == alongMomentum) ? oppositeToMomentum : alongMomentum;
0249     PropagationDirection newDir = anyDirection;
0250     HelixArbitraryPlaneCrossing  planeCrossing1B(helixPos1,helixDir1,rho1,newDir);
0251     check1 = propagateWithHelixCrossing(planeCrossing1B,plane,fts.momentum().mag(),gp1,gm1,s1);
0252     LogDebug("AnalyticalPropagator") << "after second attempt, check1, s1,gp1: "
0253     << check1 << " , "
0254     << s1 << " , " << gp1 << "\n";
0255
0256     xGen = planeCrossing1B.position(s1+tolerance);
0257     pGen = planeCrossing1B.direction(s1+tolerance);
0258     }
0259       */
0260
0261   if (!check1) {
0262     LogDebug("AnalyticalPropagator") << "failed also second attempt. No idea what to do, then bailout"
0263                                      << "\n";
0264   }
0265
0266   pGen *= gm1.mag() / pGen.mag();
0267   GlobalPoint gp2(xGen);
0268   GlobalVector gm2(pGen);
0269   double s2 = 0;
0270   double rho2 = rho1;
0271   HelixPlaneCrossing::PositionType helixPos2(gp2);
0272   HelixPlaneCrossing::DirectionType helixDir2(gm2);
0273   OptimalHelixPlaneCrossing planeCrossing2(plane, helixPos2, helixDir2, rho2, propagationDirection());
0274
0275   bool check2 = propagateWithHelixCrossing(*planeCrossing2, plane, gm2.mag(), gp2, gm2, s2);
0276
0277   if (!check2) {
0278     x = gp1;
0279     p = gm1;
0280     s = s1;
0281     return check1;
0282   }
0283
0284   if (!check1) {
0285     edm::LogError("AnalyticalPropagator") << "LOGIC ERROR: I should not have entered here!"
0286                                           << "\n";
0287     return false;
0288   }
0289
0290   LogDebug("AnalyticalPropagator") << "check2, s2, gp2: " << check2 << " , " << s2 << " , " << gp2 << "\n";
0291
0292   double dist1 = (plane.position() - gp1).perp();
0293   double dist2 = (plane.position() - gp2).perp();
0294
0295   LogDebug("AnalyticalPropagator") << "propDir, dist1, dist2: " << propagationDirection() << " , " << dist1 << " , "
0296                                    << dist2 << "\n";
0297
0298   //If there are two solutions, the one which is the closest to the module's center is chosen
0299   if (dist1 < 2 * dist2) {
0300     x = gp1;
0301     p = gm1;
0302     s = s1;
0303     return check1;
0304   } else if (dist2 < 2 * dist1) {
0305     x = gp2;
0306     p = gm2;
0307     s = s1 + s2 + tolerance;
0308     return check2;
0309   } else {
0310     if (fabs(s1) < fabs(s2)) {
0311       x = gp1;
0312       p = gm1;
0313       s = s1;
0314       return check1;
0315     } else {
0316       x = gp2;
0317       p = gm2;
0318       s = s1 + s2 + tolerance;
0319       return check2;
0320     }
0321   }
0322
0323   //-------- END of ugly piece of code  ---------------
0324 }
0325
0326 bool AnalyticalPropagator::propagateWithLineCrossing(
0327     const GlobalPoint& x0, const GlobalVector& p0, const Plane& plane, GlobalPoint& x, double& s) const {
0328   //
0329   // Instantiate auxiliary object for finding intersection.
0330   // Frame-independant point and vector are created explicitely to
0331   // avoid confusing gcc (refuses to compile with temporary objects
0332   // in the constructor).
0333   //
0334   StraightLinePlaneCrossing::PositionType pos(x0);
0335   StraightLinePlaneCrossing::DirectionType dir(p0);
0336   StraightLinePlaneCrossing planeCrossing(pos, dir, propagationDirection());
0337   //
0338   // get solution
0339   //
0340   std::pair<bool, double> propResult = planeCrossing.pathLength(plane);
0341   if (!propResult.first)
0342     return false;
0343   s = propResult.second;
0344   // point (reconverted to GlobalPoint)
0345   x = GlobalPoint(planeCrossing.position(s));
0346   //
0347   return true;
0348 }
0349
0350 bool AnalyticalPropagator::propagateWithLineCrossing(
0351     const GlobalPoint& x0, const GlobalVector& p0, const Cylinder& cylinder, GlobalPoint& x, double& s) const {
0352   //
0353   // Instantiate auxiliary object for finding intersection.
0354   // Frame-independant point and vector are created explicitely to
0355   // avoid confusing gcc (refuses to compile with temporary objects
0356   // in the constructor).
0357   //
0358   StraightLineBarrelCylinderCrossing cylCrossing(x0, p0, propagationDirection());
0359   //
0360   // get solution
0361   //
0362   std::pair<bool, double> propResult = cylCrossing.pathLength(cylinder);
0363   if (!propResult.first)
0364     return false;
0365   s = propResult.second;
0366   // point (reconverted to GlobalPoint)
0367   x = cylCrossing.position(s);
0368   //
0369   return true;
0370 }
0371
0372 bool AnalyticalPropagator::propagateWithHelixCrossing(HelixPlaneCrossing& planeCrossing,
0373                                                       const Plane& plane,
0374                                                       const float pmag,
0375                                                       GlobalPoint& x,
0376                                                       GlobalVector& p,
0377                                                       double& s) const {
0378   // get solution
0379   std::pair<bool, double> propResult = planeCrossing.pathLength(plane);
0380   if UNLIKELY (!propResult.first)
0381     return false;
0382
0383   s = propResult.second;
0384   x = GlobalPoint(planeCrossing.position(s));
0385   // direction (reconverted to GlobalVector, renormalised)
0386   GlobalVector pGen = GlobalVector(planeCrossing.direction(s));
0387   pGen *= pmag / pGen.mag();
0388   p = pGen;
0389   //
0390   return true;
0391 }