In This Package:

ParabolaTraj.cpp

Go to the documentation of this file.

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

// local
#include "Kernel/ParabolaTraj.h"
#include "GaudiKernel/SystemOfUnits.h"
using namespace LHCb;
using namespace ROOT::Math;

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

ParabolaTraj::ParabolaTraj( const Point& point,
                            const Vector& dir,
                            const Vector& curv,
                            const Range& range)
  : Trajectory(range),
    m_pos(point),
    m_dir(dir.unit()),
    m_curv(curv)
{
};

Trajectory::Point ParabolaTraj::position( double arclength ) const
{
  return m_pos + arclength * (m_dir + 0.5 * arclength * m_curv);
};

Trajectory::Vector ParabolaTraj::direction( double arclength ) const
{
  return m_dir + arclength * m_curv;
};

Trajectory::Vector ParabolaTraj::curvature( double /* arclength */) const 
{
  return m_curv;
};

void ParabolaTraj::expansion( double arclength,
                              Point& p,
                              Vector& dp,
                              Vector& ddp ) const
{
  ddp = m_curv;
  dp  = m_dir + arclength*m_curv;
  p   = m_pos + arclength* (m_dir + 0.5 * arclength * m_curv);
};

double ParabolaTraj::muEstimate( const Point& point ) const
{
  Vector r = point - m_pos;
  if (m_curv.R()<0.01*m_dir.R()) { // small curvature limit: neglect curvature
     return r.Dot(m_dir);
  }
  // get vector from m_pos to point projected into plane of parabola
  Vector normal = m_dir.Cross(m_curv).unit();
  r -= normal.Dot(r)*normal;
  // get normalized 'x' and 'y' coordinates of this vector by projecting onto the
  // axis. In terms of these, the parabola is parameterized as (arclen, arclen^2/2)
  // (i.e. arclen is actually the distance along the 'x' coordinate!)
  double x = r.Dot(m_dir);
  double y = r.Dot(m_curv);
  //  now we need to minimize the distance between (x,y) and (s,s*s/2)
  //  where s is the arclen (well, not quite, but that is what we really
  //  use in 'point(arclen)' ;-)
  //  This requires solving a 3rd order polynomial, so we assume m_curve<<1
  //  and solve a linear equation instead.
  return x/(1-y);
};

double ParabolaTraj::distTo1stError( double , double tolerance , int ) const 
{
  return std::sqrt(2*tolerance/m_curv.R());
};

double ParabolaTraj::distTo2ndError( double , double , int ) const
{
  return 10*Gaudi::Units::km;  
};

---