SimplifiedGeometryPropagator/src/HelixTrajectory.cc

0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/HelixTrajectory.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/StraightTrajectory.h"
0003 #include "FastSimulation/SimplifiedGeometryPropagator/interface/BarrelSimplifiedGeometry.h"
0004 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ForwardSimplifiedGeometry.h"
0005 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Particle.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0008 #include <cmath>
0009
0010 // helix phi definition
0011 // ranges from 0 to 2PI
0012 // 0 corresponds to the positive x direction
0013 // phi increases counterclockwise
0014
0015 fastsim::HelixTrajectory::HelixTrajectory(const fastsim::Particle& particle, double magneticFieldZ)
0016     : Trajectory(particle)
0017       // exact: r = gamma*beta*m_0*c / (q*e*B) = p_T / (q * e * B)
0018       // momentum in units of GeV/c: r = p_T * 10^9 / (c * q * B)
0019       // in cmssw units: r = p_T / (c * 10^-4 * q * B)
0020       ,
0021       radius_(
0022           std::abs(momentum_.Pt() / (fastsim::Constants::speedOfLight * 1e-4 * particle.charge() * magneticFieldZ))),
0023       phi_(std::atan(momentum_.Py() / momentum_.Px()) +
0024            (momentum_.Px() * particle.charge() < 0 ? 3. * M_PI / 2. : M_PI / 2.))
0025       // maybe consider (for -pi/2<x<pi/2)
0026       // cos(atan(x)) = 1 / sqrt(x^2+1)
0027       // -> cos(atan(x) + pi/2)  = - x / sqrt(x^2+1)
0028       // -> cos(atan(x) +3*pi/2) = + x / sqrt(x^2+1)
0029       // sin(atan(x)) = x / sqrt(x^2+1)
0030       // -> sin(atan(x) + pi/2)  = + 1 / sqrt(x^2+1)
0031       // -> sin(atan(x) +3*pi/2) = - 1 / sqrt(x^2+1)
0032       ,
0033       centerX_(position_.X() -
0034                radius_ * (momentum_.Py() / momentum_.Px()) /
0035                    std::sqrt((momentum_.Py() / momentum_.Px()) * (momentum_.Py() / momentum_.Px()) + 1) *
0036                    (momentum_.Px() * particle.charge() < 0 ? 1. : -1.)),
0037       centerY_(position_.Y() -
0038                radius_ * 1 / std::sqrt((momentum_.Py() / momentum_.Px()) * (momentum_.Py() / momentum_.Px()) + 1) *
0039                    (momentum_.Px() * particle.charge() < 0 ? -1. : 1.))
0040       //, centerX_(position_.X() - radius_*std::cos(phi_))
0041       //, centerY_(position_.Y() - radius_*std::sin(phi_))
0042       ,
0043       centerR_(std::sqrt(centerX_ * centerX_ + centerY_ * centerY_)),
0044       minR_(std::abs(centerR_ - radius_)),
0045       maxR_(centerR_ + radius_)
0046       // omega = q * e * B / (gamma * m) = q * e *B / (E / c^2) = q * e * B * c^2 / E
0047       // omega: negative for negative q -> seems to be what we want.
0048       // energy in units of GeV: omega = q * B * c^2 / (E * 10^9)
0049       // in cmssw units: omega[1/ns] = q * B * c^2 * 10^-4 / E
0050       ,
0051       phiSpeed_(-particle.charge() * magneticFieldZ * fastsim::Constants::speedOfLight *
0052                 fastsim::Constants::speedOfLight * 1e-4 / momentum_.E()) {
0053   ;
0054 }
0055
0056 bool fastsim::HelixTrajectory::crosses(const BarrelSimplifiedGeometry& layer) const {
0057   return (minR_ < layer.getRadius() && maxR_ > layer.getRadius());
0058 }
0059
0060 double fastsim::HelixTrajectory::nextCrossingTimeC(const BarrelSimplifiedGeometry& layer, bool onLayer) const {
0061   if (!crosses(layer))
0062     return -1;
0063
0064   // solve the following equation for sin(phi)
0065   // (x^2 + y^2 = R_L^2)     (1)      the layer
0066   // x = x_c + R_H*cos(phi)  (2)      the helix in the xy plane
0067   // y = y_c + R_H*sin(phi)  (3)      the helix in the xy plane
0068   // with
0069   // R_L: the radius of the layer
0070   // x_c,y_c the center of the helix in xy plane
0071   // R_H, the radius of the helix
0072   // phi, the phase of the helix
0073   //
0074   // substitute (2) and (3) in (1)
0075   // =>
0076   //   x_c^2 + 2*x_c*R_H*cos(phi) + R_H^2*cos^2(phi)
0077   // + y_c^2 + 2*y_c*R_H*sin(phi) + R_H^2*sin^2(phi)
0078   // = R_L^2
0079   // =>
0080   // (x_c^2 + y_c^2 + R_H^2 - R_L^2) + (2*y_c*R_H)*sin(phi) = -(2*x_c*R_H)*cos(phi)
0081   //
0082   // rewrite
0083   //               E                 +       F    *sin(phi) =      G     *cos(phi)
0084   // =>
0085   // E^2 + 2*E*F*sin(phi) + F^2*sin^2(phi) = G^2*(1-sin^2(phi))
0086   // rearrange
0087   // sin^2(phi)*(F^2 + G^2) + sin(phi)*(2*E*F) + (E^2 - G^2) = 0
0088   //
0089   // rewrite
0090   // sin^2(phi)*     a      + sin(phi)*   b    +      c      = 0
0091   // => sin(phi) = (-b +/- sqrt(b^2 - 4*ac)) / (2*a)
0092   // with
0093   // a = F^2 + G^2
0094   // b = 2*E*F
0095   // c = E^2 - G^2
0096
0097   double E = centerX_ * centerX_ + centerY_ * centerY_ + radius_ * radius_ - layer.getRadius() * layer.getRadius();
0098   double F = 2 * centerY_ * radius_;
0099   double G = 2 * centerX_ * radius_;
0100
0101   double a = F * F + G * G;
0102   double b = 2 * E * F;
0103   double c = E * E - G * G;
0104
0105   double delta = b * b - 4 * a * c;
0106
0107   // case of no solution
0108   if (delta < 0) {
0109     // Should not be reached: Full Propagation does always have a solution "if(crosses(layer)) == -1"
0110     // Even if particle is outside all layers -> can turn around in magnetic field
0111     return -1;
0112   }
0113
0114   // Uses a numerically more stable procedure:
0115   // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
0116   double sqrtDelta = sqrt(delta);
0117   double phi1 = 0, phi2 = 0;
0118   if (b < 0) {
0119     phi1 = std::asin((2. * c) / (-b + sqrtDelta));
0120     phi2 = std::asin((-b + sqrtDelta) / (2. * a));
0121   } else {
0122     phi1 = std::asin((-b - sqrtDelta) / (2. * a));
0123     phi2 = std::asin((2. * c) / (-b - sqrtDelta));
0124   }
0125
0126   // asin is ambiguous, make sure to have the right solution
0127   if (std::abs(layer.getRadius() - getRadParticle(phi1)) > 1.0e-2) {
0128     phi1 = -phi1 + M_PI;
0129   }
0130   if (std::abs(layer.getRadius() - getRadParticle(phi2)) > 1.0e-2) {
0131     phi2 = -phi2 + M_PI;
0132   }
0133
0134   // another ambiguity
0135   if (phi1 < 0) {
0136     phi1 += 2. * M_PI;
0137   }
0138   if (phi2 < 0) {
0139     phi2 += 2. * M_PI;
0140   }
0141
0142   // find the corresponding times when the intersection occurs
0143   // make sure they are positive
0144   double t1 = (phi1 - phi_) / phiSpeed_;
0145   while (t1 < 0) {
0146     t1 += 2 * M_PI / std::abs(phiSpeed_);
0147   }
0148   double t2 = (phi2 - phi_) / phiSpeed_;
0149   while (t2 < 0) {
0150     t2 += 2 * M_PI / std::abs(phiSpeed_);
0151   }
0152
0153   // Check if propagation successful (numerical reasons): both solutions (phi1, phi2) have to be on the layer (same radius)
0154   // Can happen due to numerical instabilities of geometrical function (if momentum is almost parallel to x/y axis)
0155   // Get crossingTimeC from StraightTrajectory as good approximation
0156   if (std::abs(layer.getRadius() - getRadParticle(phi1)) > 1.0e-2 ||
0157       std::abs(layer.getRadius() - getRadParticle(phi2)) > 1.0e-2) {
0158     StraightTrajectory traj(*this);
0159     return traj.nextCrossingTimeC(layer, onLayer);
0160   }
0161
0162   // if the particle is already on the layer, we need to make sure the 2nd solution is picked.
0163   // happens if particle turns around in the magnetic field instead of hitting the next layer
0164   if (onLayer) {
0165     bool particleMovesInwards = momentum_.X() * position_.X() + momentum_.Y() * position_.Y() < 0;
0166
0167     double posX1 = centerX_ + radius_ * std::cos(phi1);
0168     double posY1 = centerY_ + radius_ * std::sin(phi1);
0169     double momX1 = momentum_.X() * std::cos(phi1 - phi_) - momentum_.Y() * std::sin(phi1 - phi_);
0170     double momY1 = momentum_.X() * std::sin(phi1 - phi_) + momentum_.Y() * std::cos(phi1 - phi_);
0171     bool particleMovesInwards1 = momX1 * posX1 + momY1 * posY1 < 0;
0172
0173     double posX2 = centerX_ + radius_ * std::cos(phi2);
0174     double posY2 = centerY_ + radius_ * std::sin(phi2);
0175     double momX2 = momentum_.X() * std::cos(phi2 - phi_) - momentum_.Y() * std::sin(phi2 - phi_);
0176     double momY2 = momentum_.X() * std::sin(phi2 - phi_) + momentum_.Y() * std::cos(phi2 - phi_);
0177     bool particleMovesInwards2 = momX2 * posX2 + momY2 * posY2 < 0;
0178
0179     if (particleMovesInwards1 != particleMovesInwards) {
0180       return t1 * fastsim::Constants::speedOfLight;
0181     } else if (particleMovesInwards2 != particleMovesInwards) {
0182       return t2 * fastsim::Constants::speedOfLight;
0183     }
0184     // try to catch numerical issues again..
0185     else {
0186       return -1;
0187     }
0188   }
0189
0190   return std::min(t1, t2) * fastsim::Constants::speedOfLight;
0191 }
0192
0193 void fastsim::HelixTrajectory::move(double deltaTimeC) {
0194   double deltaT = deltaTimeC / fastsim::Constants::speedOfLight;
0195   double deltaPhi = phiSpeed_ * deltaT;
0196   position_.SetXYZT(centerX_ + radius_ * std::cos(phi_ + deltaPhi),
0197                     centerY_ + radius_ * std::sin(phi_ + deltaPhi),
0198                     position_.Z() + momentum_.Z() / momentum_.E() * deltaTimeC,
0199                     position_.T() + deltaT);
0200   // Rotation defined by
0201   // x' = x cos θ - y sin θ
0202   // y' = x sin θ + y cos θ
0203   momentum_.SetXYZT(momentum_.X() * std::cos(deltaPhi) - momentum_.Y() * std::sin(deltaPhi),
0204                     momentum_.X() * std::sin(deltaPhi) + momentum_.Y() * std::cos(deltaPhi),
0205                     momentum_.Z(),
0206                     momentum_.E());
0207 }
0208
0209 double fastsim::HelixTrajectory::getRadParticle(double phi) const {
0210   return sqrt((centerX_ + radius_ * std::cos(phi)) * (centerX_ + radius_ * std::cos(phi)) +
0211               (centerY_ + radius_ * std::sin(phi)) * (centerY_ + radius_ * std::sin(phi)));
0212 }