GFlash/src/GflashTrajectory.cc

0001 // Most methods of this class are excerpted from CDF Helix and Trajectory class
0002
0003 #include "SimGeneral/GFlash/interface/GflashTrajectory.h"
0004
0005 GflashTrajectory::GflashTrajectory()
0006     : _cotTheta(0.0),
0007       _curvature(0.0),
0008       _z0(0.0),
0009       _d0(0.0),
0010       _phi0(0.0),
0011       _isStale(true),
0012       _sinPhi0(2),
0013       _cosPhi0(2),
0014       _sinTheta(2),
0015       _cosTheta(2),
0016       _s(-999.999),
0017       _aa(2),
0018       _ss(2),
0019       _cc(2) {
0020   // detault constructor
0021 }
0022
0023 // GflashTrajectory::GflashTrajectory(const HepGeom::Vector3D<double>  &
0024 // MomentumGev, const HepGeom::Point3D<double>   & PositionCm,       double q,
0025 // double
0026 // BFieldTesla)
0027 void GflashTrajectory::initializeTrajectory(const HepGeom::Vector3D<double> &MomentumGev,
0028                                             const HepGeom::Point3D<double> &PositionCm,
0029                                             double q,
0030                                             double BFieldTesla) {
0031   double CotTheta = 0.0;
0032   double W = 0;
0033   double Z0 = 0;
0034   double D0 = 0;
0035   double Phi0 = 0;
0036
0037   if (BFieldTesla != 0.0 && q != 0.0) {
0038     double CurvatureConstant = 0.0029979;
0039     double Helicity = -1.0 * fabs(BFieldTesla) * fabs(q) / (BFieldTesla * q);
0040     double Radius = fabs(MomentumGev.perp() / (CurvatureConstant * BFieldTesla * q));
0041
0042     if (Radius == 0.0)
0043       W = HUGE_VAL;
0044     else
0045       W = Helicity / Radius;
0046     double phi1 = MomentumGev.phi();
0047     double x = PositionCm.x(), y = PositionCm.y(), z = PositionCm.z();
0048     double sinPhi1 = sin(phi1), cosPhi1 = cos(phi1);
0049     double gamma = atan((x * cosPhi1 + y * sinPhi1) / (x * sinPhi1 - y * cosPhi1 - 1 / W));
0050     Phi0 = phi1 + gamma;
0051     if (Phi0 > M_PI)
0052       Phi0 = Phi0 - 2.0 * M_PI;
0053     if (Phi0 < -M_PI)
0054       Phi0 = Phi0 + 2.0 * M_PI;
0055     D0 = ((1 / W + y * cosPhi1 - x * sinPhi1) / cos(gamma) - 1 / W);
0056     CotTheta = MomentumGev.z() / MomentumGev.perp();
0057     Z0 = z + gamma * CotTheta / W;
0058   } else {
0059     CLHEP::Hep3Vector direction = MomentumGev.unit();
0060     CLHEP::Hep3Vector projectedDirection = CLHEP::Hep3Vector(direction.x(), direction.y(), 0.0).unit();
0061     double s = projectedDirection.dot(PositionCm);
0062     double sprime = s / sin(direction.theta());
0063     Z0 = (PositionCm - sprime * direction).z();
0064     Phi0 = MomentumGev.phi();
0065     CotTheta = MomentumGev.z() / MomentumGev.perp();
0066     W = 0.0;
0067     D0 = (PositionCm.y() * cos(Phi0) - PositionCm.x() * sin(Phi0));
0068   }
0069
0070   _cotTheta = CotTheta;
0071   _curvature = W / 2;
0072   _z0 = Z0;
0073   _d0 = D0;
0074   _phi0 = Phi0;
0075
0076   _isStale = true;
0077   _s = -999.999;
0078   _aa = -999.999;
0079   _ss = -999.999;
0080   _cc = -999.999;
0081   _sinPhi0 = 1.0;
0082   _cosPhi0 = 1.0;
0083   _sinTheta = 1.0;
0084   _cosTheta = 1.0;
0085 }
0086
0087 GflashTrajectory::~GflashTrajectory() {}
0088
0089 void GflashTrajectory::setCotTheta(double cotTheta) {
0090   _cotTheta = cotTheta;
0091   _isStale = true;
0092 }
0093
0094 void GflashTrajectory::setCurvature(double curvature) {
0095   _curvature = curvature;
0096   _isStale = true;
0097 }
0098
0099 void GflashTrajectory::setZ0(double z0) {
0100   _z0 = z0;
0101   _isStale = true;
0102 }
0103
0104 void GflashTrajectory::setD0(double d0) {
0105   _d0 = d0;
0106   _isStale = true;
0107 }
0108
0109 void GflashTrajectory::setPhi0(double phi0) {
0110   _phi0 = phi0;
0111   _isStale = true;
0112 }
0113
0114 double GflashTrajectory::getSinPhi0() const {
0115   _refreshCache();
0116   return _sinPhi0;
0117 }
0118 double GflashTrajectory::getCosPhi0() const {
0119   _refreshCache();
0120   return _cosPhi0;
0121 }
0122 double GflashTrajectory::getSinTheta() const {
0123   _refreshCache();
0124   return _sinTheta;
0125 }
0126 double GflashTrajectory::getCosTheta() const {
0127   _refreshCache();
0128   return _cosTheta;
0129 }
0130
0131 HepGeom::Point3D<double> GflashTrajectory::getPosition(double s) const {
0132   _cacheSinesAndCosines(s);
0133   if (s == 0.0 || _curvature == 0.0) {
0134     return HepGeom::Point3D<double>(
0135         -_d0 * _sinPhi0 + s * _cosPhi0 * _sinTheta, _d0 * _cosPhi0 + s * _sinPhi0 * _sinTheta, _z0 + s * _cosTheta);
0136   } else {
0137     return HepGeom::Point3D<double>(
0138         (_cosPhi0 * _ss - _sinPhi0 * (2.0 * _curvature * _d0 + 1.0 - _cc)) / (2.0 * _curvature),
0139         (_sinPhi0 * _ss + _cosPhi0 * (2.0 * _curvature * _d0 + 1.0 - _cc)) / (2.0 * _curvature),
0140         _s * _cosTheta + _z0);
0141   }
0142 }
0143
0144 HepGeom::Vector3D<double> GflashTrajectory::getDirection(double s) const {
0145   _cacheSinesAndCosines(s);
0146   if (s == 0.0) {
0147     return HepGeom::Vector3D<double>(_cosPhi0 * _sinTheta, _sinPhi0 * _sinTheta, _cosTheta);
0148   }
0149   double xtan = _sinTheta * (_cosPhi0 * _cc - _sinPhi0 * _ss);
0150   double ytan = _sinTheta * (_cosPhi0 * _ss + _sinPhi0 * _cc);
0151   double ztan = _cosTheta;
0152   return HepGeom::Vector3D<double>(xtan, ytan, ztan);
0153 }
0154
0155 void GflashTrajectory::getGflashTrajectoryPoint(GflashTrajectoryPoint &point, double s) const {
0156   _cacheSinesAndCosines(s);
0157
0158   double cP0sT = _cosPhi0 * _sinTheta, sP0sT = _sinPhi0 * _sinTheta;
0159   if (s && _curvature) {
0160     point.getPosition().set((_cosPhi0 * _ss - _sinPhi0 * (2.0 * _curvature * _d0 + 1.0 - _cc)) / (2.0 * _curvature),
0161                             (_sinPhi0 * _ss + _cosPhi0 * (2.0 * _curvature * _d0 + 1.0 - _cc)) / (2.0 * _curvature),
0162                             s * _cosTheta + _z0);
0163
0164     point.getMomentum().set(cP0sT * _cc - sP0sT * _ss, cP0sT * _ss + sP0sT * _cc, _cosTheta);
0165     point.setPathLength(s);
0166   } else {
0167     point.getPosition().set(-_d0 * _sinPhi0 + s * cP0sT, _d0 * _cosPhi0 + s * sP0sT, _z0 + s * _cosTheta);
0168
0169     point.getMomentum().set(cP0sT, sP0sT, _cosTheta);
0170
0171     point.setPathLength(s);
0172   }
0173 }
0174
0175 double GflashTrajectory::getPathLengthAtRhoEquals(double rho) const {
0176   return (getSinTheta() ? (getL2DAtR(rho) / getSinTheta()) : 0.0);
0177 }
0178
0179 double GflashTrajectory::getPathLengthAtZ(double z) const {
0180   return (getCosTheta() ? (z - getZ0()) / getCosTheta() : 0.0);
0181 }
0182
0183 double GflashTrajectory::getZAtR(double rho) const { return _z0 + getCotTheta() * getL2DAtR(rho); }
0184
0185 double GflashTrajectory::getL2DAtR(double rho) const {
0186   double L2D;
0187
0188   double c = getCurvature();
0189   double d = getD0();
0190
0191   if (c != 0.0) {
0192     double rad = (rho * rho - d * d) / (1.0 + 2.0 * c * d);
0193     double rprime;
0194     if (rad < 0.0) {
0195       rprime = 0.0;
0196     } else {
0197       rprime = sqrt(rad);
0198     }
0199     if (c * rprime > 1.0 || c * rprime < -1.0) {
0200       L2D = c * rprime > 0. ? M_PI / c : -M_PI / c;
0201     } else
0202       L2D = asin(c * rprime) / c;
0203   } else {
0204     double rad = rho * rho - d * d;
0205     double rprime;
0206     if (rad < 0.0)
0207       rprime = 0.0;
0208     else
0209       rprime = sqrt(rad);
0210
0211     L2D = rprime;
0212   }
0213   return L2D;
0214 }
0215
0216 // Update _sinTheta,_cosTheta,_sinPhi0, and _cosPhi0
0217 void GflashTrajectory::_refreshCache() const {
0218   if (_isStale) {
0219     _isStale = false;
0220     double theta;
0221     if (_cotTheta == 0.0) {
0222       theta = M_PI / 2.0;
0223     } else {
0224       theta = atan(1.0 / _cotTheta);
0225       if (theta < 0.0)
0226         theta += M_PI;
0227     }
0228     if (theta == 0.0) {
0229       _sinTheta = 0.0;
0230       _cosTheta = 1.0;
0231     } else {
0232       _cosTheta = cos(theta);
0233       _sinTheta = sqrt(1 - _cosTheta * _cosTheta);
0234     }
0235     if (_phi0 == 0.0) {
0236       _sinPhi0 = 0.0;
0237       _cosPhi0 = 1.0;
0238     } else {
0239       _cosPhi0 = cos(_phi0);
0240       _sinPhi0 = sin(_phi0);
0241       //      _sinPhi0 = sqrt(1.0-_cosPhi0*_cosPhi0);
0242       //      if (_phi0>M_PI) _sinPhi0 = -_sinPhi0;
0243     }
0244   }
0245 }
0246 // Update _s, _aa, _ss, and _cc if the arclength has changed.
0247 void GflashTrajectory::_cacheSinesAndCosines(double s) const {
0248   _refreshCache();
0249   if (_s != s) {
0250     _s = s;
0251     _aa = 2.0 * _s * _curvature * _sinTheta;
0252     if (_aa == 0.0) {
0253       _ss = 0.0;
0254       _cc = 1.0;
0255     } else {
0256       _ss = sin(_aa);
0257       _cc = cos(_aa);
0258     }
0259   }
0260 }