In This Package:

CircleTraj.cpp

Go to the documentation of this file.

// $Id: CircleTraj.cpp,v 1.12 2007/10/16 11:50:59 wouter Exp $
// Include files

// local
#include "Kernel/CircleTraj.h"

#include "Math/GenVector/AxisAngle.h"
#include <math.h>

using namespace LHCb;
using namespace ROOT::Math;
using namespace Gaudi;

std::auto_ptr<Trajectory> CircleTraj::clone() const
{
  return std::auto_ptr<Trajectory>(new CircleTraj(*this));
}

CircleTraj::CircleTraj( const Point& origin,
                        const Vector& dir1,
                        const Vector& dir2,
                        double radius)
  : Trajectory(0.,radius*std::asin((dir1.unit()).Cross(dir2.unit()).r())),
    m_origin(origin),
    m_normal(dir1.Cross(dir2).unit()),
    m_dirStart(dir1.unit()),
    m_radius(radius)
{
};

CircleTraj::CircleTraj( const Point& origin,
                        const Vector& normal,
                        const Vector& origin2point,
                        const Range& range)
  : Trajectory(range),
    m_origin(origin),
    m_normal(normal.unit()),
    m_dirStart(origin2point-origin2point.Dot(m_normal)*m_normal),
    m_radius(m_dirStart.r())
{
};

Trajectory::Point CircleTraj::position( double s ) const
{
  return m_origin+m_radius*AxisAngle(m_normal,s/m_radius)(m_dirStart);
};

Trajectory::Vector CircleTraj::direction( double s ) const
{
  return m_normal.Cross(AxisAngle(m_normal,s/m_radius)(m_dirStart));
};

Trajectory::Vector CircleTraj::curvature( double s ) const 
{
  return (-1.0/m_radius)*AxisAngle(m_normal,s/m_radius)(m_dirStart);
};

void CircleTraj::expansion( double s,
                            Point& p,
                            Vector& dp,
                            Vector& ddp ) const
{
  Vector r(  AxisAngle(m_normal,s/m_radius)(m_dirStart) );
  ddp =  (-1.0/m_radius)*r;
  dp  = m_normal.Cross(r);
  p   = m_origin+m_radius*r;
};

double CircleTraj::muEstimate( const Point& point ) const
{
  // get vector from origin, to point after projecting it 
  // into the plane of the circle. (i.e. this vector is normal
  // to m_normal)
  Vector r( (point - m_normal.Dot(point-m_origin)*m_normal)-m_origin );

  // Determine delta phi angle between arclength=0 angle and angle of r
  double dphi = r.phi() - m_dirStart.phi();
  
  // Check whether angle outside of [-pi/2,+pi/2]
  if( m_dirStart.Dot( r ) < 0) {
    if( dphi > M_PI ) dphi -= 2*M_PI;
    else              dphi += 2*M_PI;
  }

  return m_radius * dphi;
};

double CircleTraj::distTo1stError( double /*arclen*/, double tolerance, int /*direction*/) const 
{
  // require 2nd order term to be less than tolerance...
  return std::sqrt(2*tolerance*m_radius);
};

double CircleTraj::distTo2ndError( double /*arclen*/, double tolerance , int /*direction*/ ) const
{
  // require 3rd order term to be less than tolerance
  // cbrt is in the C99 standard -- hope it is available on all platforms...
  // return cbrt(6*tolerance*m_radius*m_radius);
  // cbrt is NOT available on windows at this time...
  return pow(6*tolerance*m_radius*m_radius,double(1)/3);
};

---