HelixBarrelCylinderCrossing.cc

Go to the documentation of this file.

#include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
#include "TrackingTools/GeomPropagators/src/RealQuadEquation.h"
#include "TrackingTools/GeomPropagators/interface/StraightLineCylinderCrossing.h"

#include "DataFormats/GeometrySurface/interface/Cylinder.h"

#include <iostream>
#include <cmath>

#include <tuple>

template <typename T>
inline T sqr(const T& t) {
  return t * t;
}

HelixBarrelCylinderCrossing::HelixBarrelCylinderCrossing(const GlobalPoint& startingPos,
                                                         const GlobalVector& startingDir,
                                                         double rho,
                                                         PropagationDirection propDir,
                                                         const Cylinder& cyl,
                                                         Solution sol) {
  // assumes the cylinder is centered at 0,0
  double R = cyl.radius();

  // protect for zero curvature case
  const double sraightLineCutoff = 1.e-7;
  if (fabs(rho) * R < sraightLineCutoff && fabs(rho) * startingPos.perp() < sraightLineCutoff) {
    // switch to straight line case
    StraightLineCylinderCrossing slc(cyl.toLocal(startingPos), cyl.toLocal(startingDir), propDir);
    std::pair<bool, double> pl = slc.pathLength(cyl);
    if (pl.first) {
      theSolExists = true;
      theS = pl.second;
      thePos = cyl.toGlobal(slc.position(theS));
      theDir = startingDir;
    } else
      theSolExists = false;
    return;  // all needed data members have been set
  }

  double R2cyl = R * R;
  double pt = startingDir.perp();
  Point center(startingPos.x() - startingDir.y() / (pt * rho), startingPos.y() + startingDir.x() / (pt * rho));
  double p2 = startingPos.perp2();
  bool solveForX;
  double B, C, E, F;
  if (fabs(center.x()) > fabs(center.y())) {
    solveForX = false;
    E = (R2cyl - p2) / (2. * center.x());
    F = center.y() / center.x();
    B = 2. * (startingPos.y() - F * startingPos.x() - E * F);
    C = 2. * E * startingPos.x() + E * E + p2 - R2cyl;
  } else {
    solveForX = true;
    E = (R2cyl - p2) / (2. * center.y());
    F = center.x() / center.y();
    B = 2. * (startingPos.x() - F * startingPos.y() - E * F);
    C = 2. * E * startingPos.y() + E * E + p2 - R2cyl;
  }

  RealQuadEquation eq(1 + F * F, B, C);
  if (!eq.hasSolution) {
    theSolExists = false;
    return;
  }

  Vector d1, d2;
  ;
  if (solveForX) {
    d1 = Point(eq.first, E - F * eq.first);
    d2 = Point(eq.second, E - F * eq.second);
  } else {
    d1 = Point(E - F * eq.first, eq.first);
    d2 = Point(E - F * eq.second, eq.second);
  }

  Vector theD;
  int theActualDir;

  std::tie(theD, theActualDir) = chooseSolution(d1, d2, startingPos, startingDir, propDir);
  if (!theSolExists)
    return;

  float ipabs = 1.f / startingDir.mag();
  float sinTheta = float(pt) * ipabs;
  float cosTheta = startingDir.z() * ipabs;

  // -------

  auto dMag = theD.mag();
  float tmp = 0.5f * float(dMag * rho);
  if (std::abs(tmp) > 1.f)
    tmp = std::copysign(1.f, tmp);
  theS = theActualDir * 2.f * std::asin(tmp) / (float(rho) * sinTheta);
  thePos = GlobalPoint(startingPos.x() + theD.x(), startingPos.y() + theD.y(), startingPos.z() + theS * cosTheta);

  if (sol == onlyPos)
    return;

  if (theS < 0)
    tmp = -tmp;
  auto sinPhi = 2.f * tmp * sqrt(1.f - tmp * tmp);
  auto cosPhi = 1.f - 2.f * tmp * tmp;
  theDir = DirectionType(startingDir.x() * cosPhi - startingDir.y() * sinPhi,
                         startingDir.x() * sinPhi + startingDir.y() * cosPhi,
                         startingDir.z());

  if (sol != bothSol)
    return;

  //-----  BM
  //double momProj1 = startingDir.x()*d1.x() + startingDir.y()*d1.y();
  //double momProj2 = startingDir.x()*d2.x() + startingDir.y()*d2.y();

  int theActualDir1 = propDir == alongMomentum ? 1 : -1;
  int theActualDir2 = propDir == alongMomentum ? 1 : -1;

  auto dMag1 = d1.mag();
  auto tmp1 = 0.5f * dMag1 * float(rho);
  if (std::abs(tmp1) > 1.f)
    tmp1 = std::copysign(1.f, tmp1);
  auto theS1 = theActualDir1 * 2.f * std::asin(tmp1) / (rho * sinTheta);
  thePos1 = GlobalPoint(startingPos.x() + d1.x(), startingPos.y() + d1.y(), startingPos.z() + theS1 * cosTheta);

  auto dMag2 = d2.mag();
  auto tmp2 = 0.5f * dMag2 * float(rho);
  if (std::abs(tmp2) > 1.f)
    tmp2 = std::copysign(1.f, tmp2);
  auto theS2 = theActualDir2 * 2.f * std::asin(tmp2) / (float(rho) * sinTheta);
  thePos2 = GlobalPoint(startingPos.x() + d2.x(), startingPos.y() + d2.y(), startingPos.z() + theS2 * cosTheta);
}

std::pair<HelixBarrelCylinderCrossing::Vector, int> HelixBarrelCylinderCrossing::chooseSolution(
    const Vector& d1,
    const Vector& d2,
    const PositionType& startingPos,
    const DirectionType& startingDir,
    PropagationDirection propDir) {
  Vector theD;
  int theActualDir;

  auto momProj1 = startingDir.x() * d1.x() + startingDir.y() * d1.y();
  auto momProj2 = startingDir.x() * d2.x() + startingDir.y() * d2.y();

  if (propDir == anyDirection) {
    theSolExists = true;
    if (d1.mag2() < d2.mag2()) {
      theD = d1;
      theActualDir = (momProj1 > 0) ? 1 : -1;
    } else {
      theD = d2;
      theActualDir = (momProj2 > 0) ? 1 : -1;
    }
  } else {
    int propSign = propDir == alongMomentum ? 1 : -1;
    if (momProj1 * momProj2 < 0) {
      // if different signs return the positive one
      theSolExists = true;
      theD = (momProj1 * propSign > 0) ? d1 : d2;
      theActualDir = propSign;
    } else if (momProj1 * propSign > 0) {
      // if both positive, return the shortest
      theSolExists = true;
      theD = (d1.mag2() < d2.mag2()) ? d1 : d2;
      theActualDir = propSign;
    } else
      theSolExists = false;
  }

  return std::pair<Vector, int>(theD, theActualDir);
}