StraightLineCylinderCrossing.cc

Go to the documentation of this file.

#include "TrackingTools/GeomPropagators/interface/StraightLineCylinderCrossing.h"

#include "DataFormats/GeometrySurface/interface/Cylinder.h"
#include "TrackingTools/GeomPropagators/src/RealQuadEquation.h"

#include <cmath>
#include <iostream>
using namespace std;

StraightLineCylinderCrossing::
StraightLineCylinderCrossing( const LocalPoint& startingPos, const LocalVector& startingDir,
                  const PropagationDirection propDir, double tolerance) :
  theX0(startingPos),
  theP0(startingDir.unit()),
  thePropDir(propDir), 
  theTolerance(tolerance) {}

std::pair<bool,double>
StraightLineCylinderCrossing::pathLength (const Cylinder& cylinder) const
{
  //
  // radius of cylinder and transversal position relative to axis
  //
  double R(cylinder.radius());
  PositionType2D xt2d(theX0.x(),theX0.y());
  //
  // transverse direction
  // 
  DirectionType2D pt2d(theP0.x(),theP0.y());
  //
  // solution of quadratic equation for s - assume |theP0|=1
  //
  RealQuadEquation eq(pt2d.mag2(),2.*xt2d.dot(pt2d),xt2d.mag2()-R*R);
  if ( !eq.hasSolution ) {
    /*
    double A=   pt2d.mag2(); 
    double B=   2.*xt2d.dot(pt2d);
    double C=   xt2d.mag2()-R*R;
    cout << "A= " << pt2d.mag2() 
     << " B= " << 2.*xt2d.dot(pt2d)
     << " C= " << xt2d.mag2()-R*R
     << " D= " << B*B - 4*A*C
     << endl;
    */
    return std::pair<bool,double>(false,0.);
  }
  //
  // choice of solution and verification of direction
  //
  return chooseSolution(eq.first,eq.second);
}

std::pair<bool,double>
StraightLineCylinderCrossing::chooseSolution (const double s1, 
                          const double s2) const
{
  //
  // follows the logic implemented in HelixBarrelCylinderCrossing
  //
  if ( thePropDir==anyDirection ) {
    return std::pair<bool,double>(true,(fabs(s1)<fabs(s2)?s1:s2));
  }
  else {
    int propSign = thePropDir==alongMomentum ? 1 : -1;
    if ( s1*s2 < 0) {
      // if different signs return the positive one
      return std::pair<bool,double>(true,((s1*propSign>0)?s1:s2));
    }
    else if ( s1*propSign>0 ) {
      // if both positive, return the shortest
      return std::pair<bool,double>(true,(fabs(s1)<fabs(s2)?s1:s2));
    }
    else {
      // if both negative, check if the smaller (abs value) is smaller then tolerance
      double shorter = std::min( fabs(s1), fabs(s2));
      if (shorter < theTolerance) return std::pair<bool,double>(true,0);
      else                        return std::pair<bool,double>(false,0.);
    }
  }
}