#include <HelixBarrelCylinderCrossing.h>

Public Types
typedef GlobalVector	DirectionType
typedef GlobalPoint	PositionType
Public Member Functions
DirectionType	direction () const
bool	hasSolution () const
	HelixBarrelCylinderCrossing (const GlobalPoint &startingPos, const GlobalVector &startingDir, double rho, PropagationDirection propDir, const Cylinder &cyl)
double	pathLength () const
PositionType	position () const
PositionType	position1 () const
	Method to access separately each solution of the helix-cylinder crossing equations.
PositionType	position2 () const
	Method to access separately each solution of the helix-cylinder crossing equations.
Private Types
typedef Basic2DVector< TmpType >	Point
typedef double	TmpType
typedef Basic2DVector< TmpType >	Vector
Private Member Functions
void	chooseSolution (const Point &p1, const Point &p2, const PositionType &startingPos, const DirectionType &startingDir, PropagationDirection propDir) dso_internal
Private Attributes
int	theActualDir
Vector	theD
DirectionType	theDir
PositionType	thePos
PositionType	thePos1
PositionType	thePos2
double	theS
bool	theSolExists

Detailed Description

Calculates the crossing of a helix with a barrel cylinder.

Definition at line 16 of file HelixBarrelCylinderCrossing.h.

Member Typedef Documentation

typedef GlobalVector HelixBarrelCylinderCrossing::DirectionType

Definition at line 25 of file HelixBarrelCylinderCrossing.h.

typedef Basic2DVector<TmpType> HelixBarrelCylinderCrossing::Point [private]

Definition at line 19 of file HelixBarrelCylinderCrossing.h.

typedef GlobalPoint HelixBarrelCylinderCrossing::PositionType

Definition at line 24 of file HelixBarrelCylinderCrossing.h.

typedef double HelixBarrelCylinderCrossing::TmpType [private]

Definition at line 18 of file HelixBarrelCylinderCrossing.h.

typedef Basic2DVector<TmpType> HelixBarrelCylinderCrossing::Vector [private]

Definition at line 20 of file HelixBarrelCylinderCrossing.h.

Constructor & Destructor Documentation

HelixBarrelCylinderCrossing::HelixBarrelCylinderCrossing	(	const GlobalPoint &	startingPos,
		const GlobalVector &	startingDir,
		double	rho,
		PropagationDirection	propDir,
		const Cylinder &	cyl
	)

Definition at line 14 of file HelixBarrelCylinderCrossing.cc.

References abs, alongMomentum, funct::C, chooseSolution(), Vispa::Plugins::EdmBrowser::EdmDataAccessor::eq(), F(), RealQuadEquation::first, RealQuadEquation::hasSolution, Basic2DVector< T >::mag(), PV3DBase< T, PVType, FrameType >::mag(), p2, StraightLineCylinderCrossing::pathLength(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::perp2(), dttmaxenums::R, Cylinder::radius(), rho, RealQuadEquation::second, mathSSE::sqrt(), theActualDir, theD, theDir, thePos, thePos1, thePos2, theS, theSolExists, tmp, Surface::toGlobal(), GloballyPositioned< T >::toLocal(), PV3DBase< T, PVType, FrameType >::x(), Basic2DVector< T >::x(), PV3DBase< T, PVType, FrameType >::y(), Basic2DVector< T >::y(), and PV3DBase< T, PVType, FrameType >::z().

{
  // 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);
  }
  
  chooseSolution(d1, d2, startingPos, startingDir, propDir);
  if (!theSolExists) return;

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


  //-----  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;


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


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

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

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

Member Function Documentation

void HelixBarrelCylinderCrossing::chooseSolution	(	const Point &	p1,
		const Point &	p2,
		const PositionType &	startingPos,
		const DirectionType &	startingDir,
		PropagationDirection	propDir
	)		`[private]`

Definition at line 129 of file HelixBarrelCylinderCrossing.cc.

References alongMomentum, anyDirection, Basic2DVector< T >::mag2(), theActualDir, theD, theSolExists, PV3DBase< T, PVType, FrameType >::x(), Basic2DVector< T >::x(), PV3DBase< T, PVType, FrameType >::y(), and Basic2DVector< T >::y().

Referenced by HelixBarrelCylinderCrossing().

{
  double momProj1 = startingDir.x()*d1.x() + startingDir.y()*d1.y();
  double 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;
  }
}

DirectionType HelixBarrelCylinderCrossing::direction ( ) const [inline]

Returns the direction along the helix that corresponds to path length "s" from the starting point. As for position, the direction of the crossing with a cylinder (if it exists!) is given by direction( pathLength( cylinder)).

Definition at line 60 of file HelixBarrelCylinderCrossing.h.

References theDir.

{ return theDir;}

bool HelixBarrelCylinderCrossing::hasSolution ( ) const [inline]

Definition at line 32 of file HelixBarrelCylinderCrossing.h.

References theSolExists.

{ return theSolExists;}

double HelixBarrelCylinderCrossing::pathLength ( ) const [inline]

Propagation status (true if valid) and (signed) path length along the helix from the starting point to the cylinder. The starting point and the cylinder are given in the constructor.

Definition at line 38 of file HelixBarrelCylinderCrossing.h.

References theS.

{ return theS;}

PositionType HelixBarrelCylinderCrossing::position ( void ) const [inline]

Returns the position along the helix that corresponds to path length "s" from the starting point. If s is obtained from the pathLength method the position is the destination point, i.e. the position of the crossing with a cylinder (if it exists!) is given by position( pathLength( cylinder)).

Definition at line 46 of file HelixBarrelCylinderCrossing.h.

References thePos.