GeomPropagators/src/HelixBarrelCylinderCrossing.cc

0001 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
0002 #include "TrackingTools/GeomPropagators/src/RealQuadEquation.h"
0003 #include "TrackingTools/GeomPropagators/interface/StraightLineCylinderCrossing.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0006
0007 #include <iostream>
0008 #include <cmath>
0009
0010 #include <tuple>
0011
0012 template <typename T>
0013 inline T sqr(const T& t) {
0014   return t * t;
0015 }
0016
0017 HelixBarrelCylinderCrossing::HelixBarrelCylinderCrossing(const GlobalPoint& startingPos,
0018                                                          const GlobalVector& startingDir,
0019                                                          double rho,
0020                                                          PropagationDirection propDir,
0021                                                          const Cylinder& cyl,
0022                                                          Solution sol) {
0023   // assumes the cylinder is centered at 0,0
0024   double R = cyl.radius();
0025
0026   // protect for zero curvature case
0027   const double sraightLineCutoff = 1.e-7;
0028   if (fabs(rho) * R < sraightLineCutoff && fabs(rho) * startingPos.perp() < sraightLineCutoff) {
0029     // switch to straight line case
0030     StraightLineCylinderCrossing slc(cyl.toLocal(startingPos), cyl.toLocal(startingDir), propDir);
0031     std::pair<bool, double> pl = slc.pathLength(cyl);
0032     if (pl.first) {
0033       theSolExists = true;
0034       theS = pl.second;
0035       thePos = cyl.toGlobal(slc.position(theS));
0036       theDir = startingDir;
0037     } else
0038       theSolExists = false;
0039     return;  // all needed data members have been set
0040   }
0041
0042   double R2cyl = R * R;
0043   double pt = startingDir.perp();
0044   Point center(startingPos.x() - startingDir.y() / (pt * rho), startingPos.y() + startingDir.x() / (pt * rho));
0045   double p2 = startingPos.perp2();
0046   bool solveForX;
0047   double B, C, E, F;
0048   if (fabs(center.x()) > fabs(center.y())) {
0049     solveForX = false;
0050     E = (R2cyl - p2) / (2. * center.x());
0051     F = center.y() / center.x();
0052     B = 2. * (startingPos.y() - F * startingPos.x() - E * F);
0053     C = 2. * E * startingPos.x() + E * E + p2 - R2cyl;
0054   } else {
0055     solveForX = true;
0056     E = (R2cyl - p2) / (2. * center.y());
0057     F = center.x() / center.y();
0058     B = 2. * (startingPos.x() - F * startingPos.y() - E * F);
0059     C = 2. * E * startingPos.y() + E * E + p2 - R2cyl;
0060   }
0061
0062   RealQuadEquation eq(1 + F * F, B, C);
0063   if (!eq.hasSolution) {
0064     theSolExists = false;
0065     return;
0066   }
0067
0068   Vector d1, d2;
0069   ;
0070   if (solveForX) {
0071     d1 = Point(eq.first, E - F * eq.first);
0072     d2 = Point(eq.second, E - F * eq.second);
0073   } else {
0074     d1 = Point(E - F * eq.first, eq.first);
0075     d2 = Point(E - F * eq.second, eq.second);
0076   }
0077
0078   Vector theD;
0079   int theActualDir;
0080
0081   std::tie(theD, theActualDir) = chooseSolution(d1, d2, startingPos, startingDir, propDir);
0082   if (!theSolExists)
0083     return;
0084
0085   float ipabs = 1.f / startingDir.mag();
0086   float sinTheta = float(pt) * ipabs;
0087   float cosTheta = startingDir.z() * ipabs;
0088
0089   // -------
0090
0091   auto dMag = theD.mag();
0092   float tmp = 0.5f * float(dMag * rho);
0093   if (std::abs(tmp) > 1.f)
0094     tmp = std::copysign(1.f, tmp);
0095   theS = theActualDir * 2.f * std::asin(tmp) / (float(rho) * sinTheta);
0096   thePos = GlobalPoint(startingPos.x() + theD.x(), startingPos.y() + theD.y(), startingPos.z() + theS * cosTheta);
0097
0098   if (sol == onlyPos)
0099     return;
0100
0101   if (theS < 0)
0102     tmp = -tmp;
0103   auto sinPhi = 2.f * tmp * sqrt(1.f - tmp * tmp);
0104   auto cosPhi = 1.f - 2.f * tmp * tmp;
0105   theDir = DirectionType(startingDir.x() * cosPhi - startingDir.y() * sinPhi,
0106                          startingDir.x() * sinPhi + startingDir.y() * cosPhi,
0107                          startingDir.z());
0108
0109   if (sol != bothSol)
0110     return;
0111
0112   //-----  BM
0113   //double momProj1 = startingDir.x()*d1.x() + startingDir.y()*d1.y();
0114   //double momProj2 = startingDir.x()*d2.x() + startingDir.y()*d2.y();
0115
0116   int theActualDir1 = propDir == alongMomentum ? 1 : -1;
0117   int theActualDir2 = propDir == alongMomentum ? 1 : -1;
0118
0119   auto dMag1 = d1.mag();
0120   auto tmp1 = 0.5f * dMag1 * float(rho);
0121   if (std::abs(tmp1) > 1.f)
0122     tmp1 = std::copysign(1.f, tmp1);
0123   auto theS1 = theActualDir1 * 2.f * std::asin(tmp1) / (rho * sinTheta);
0124   thePos1 = GlobalPoint(startingPos.x() + d1.x(), startingPos.y() + d1.y(), startingPos.z() + theS1 * cosTheta);
0125
0126   auto dMag2 = d2.mag();
0127   auto tmp2 = 0.5f * dMag2 * float(rho);
0128   if (std::abs(tmp2) > 1.f)
0129     tmp2 = std::copysign(1.f, tmp2);
0130   auto theS2 = theActualDir2 * 2.f * std::asin(tmp2) / (float(rho) * sinTheta);
0131   thePos2 = GlobalPoint(startingPos.x() + d2.x(), startingPos.y() + d2.y(), startingPos.z() + theS2 * cosTheta);
0132 }
0133
0134 std::pair<HelixBarrelCylinderCrossing::Vector, int> HelixBarrelCylinderCrossing::chooseSolution(
0135     const Vector& d1,
0136     const Vector& d2,
0137     const PositionType& startingPos,
0138     const DirectionType& startingDir,
0139     PropagationDirection propDir) {
0140   Vector theD;
0141   int theActualDir;
0142
0143   auto momProj1 = startingDir.x() * d1.x() + startingDir.y() * d1.y();
0144   auto momProj2 = startingDir.x() * d2.x() + startingDir.y() * d2.y();
0145
0146   if (propDir == anyDirection) {
0147     theSolExists = true;
0148     if (d1.mag2() < d2.mag2()) {
0149       theD = d1;
0150       theActualDir = (momProj1 > 0) ? 1 : -1;
0151     } else {
0152       theD = d2;
0153       theActualDir = (momProj2 > 0) ? 1 : -1;
0154     }
0155   } else {
0156     int propSign = propDir == alongMomentum ? 1 : -1;
0157     if (momProj1 * momProj2 < 0) {
0158       // if different signs return the positive one
0159       theSolExists = true;
0160       theD = (momProj1 * propSign > 0) ? d1 : d2;
0161       theActualDir = propSign;
0162     } else if (momProj1 * propSign > 0) {
0163       // if both positive, return the shortest
0164       theSolExists = true;
0165       theD = (d1.mag2() < d2.mag2()) ? d1 : d2;
0166       theActualDir = propSign;
0167     } else
0168       theSolExists = false;
0169   }
0170
0171   return std::pair<Vector, int>(theD, theActualDir);
0172 }