TrackTrajectory.cxx

Go to the documentation of this file.

//==============================================
// TrackTrajectory class implementation
//==============================================

#include "AtlfastEvent/TrackTrajectory.h"
#include <iostream>
#include <cmath>
#include "GaudiKernel/MsgStream.h"
#include <iomanip>

namespace Atlfast {
  using std::abs;  
  using std::setw;
  using std::setprecision;
  using std::setiosflags;
  using std::ios;

std::ostream& operator << ( std::ostream& s, const TrackTrajectory& traj ) 
{
  TrackParameters param = traj.parameters();

    s << "TrackTrajectory"<<std::endl; 

    s<<setiosflags(ios::fixed);
    s<<setprecision(7);
    s<<setw(20);

  s   <<setw(20)<< param.pT()
      << "        pT"<<std::endl;
  s    <<setw(20)<< param.eta()
      << "        eta"<<std::endl;

  s << "----------------------------" << std::endl;
 
  s    <<setw(20)<< param.impactParameter()
      << "        d0"<<std::endl; 
  s    <<setw(20)<< param.zPerigee()
      << "        z0" << std::endl;  
  s   <<setw(20)<< param.phi().val()
      << "       phi" <<std::endl;
  s    <<setw(20)<<param.cotTheta()
      << " cot Theta" << std::endl;
  s    <<setw(20)<< param.invPtCharge()
      << "  1/(q.PT)" <<std::endl;
  s << "----------------------------" << std::endl;
  s    <<setw(20)<< traj.radius()
      << "    radius" <<std::endl;
  s    <<setw(20)<< traj.curvature()
      << " curvature" <<std::endl<<std::endl;

  return s;
}

MsgStream& operator << ( MsgStream& s, const TrackTrajectory& traj ) 
{
  TrackParameters param = traj.parameters();
  
  s << "TrackTrajectory"<<endreq; 

    s<<setiosflags(ios::fixed);
    s<<setprecision(7);
    s<<setw(20);

  s   <<setw(20)<< param.pT()
      << "        pT"<<endreq;
  s    <<setw(20)<< param.eta()
      << "        eta"<<endreq;

  s << "----------------------------" << endreq;
 
  s    <<setw(20)<< param.impactParameter()
      << "        d0"<<endreq; 
  s    <<setw(20)<< param.zPerigee()
      << "        z0" << endreq;  
  s   <<setw(20)<< param.phi().val()
      << "       phi" <<endreq;
  s    <<setw(20)<<param.cotTheta()
      << " cot Theta" << endreq;
  s    <<setw(20)<< param.invPtCharge()
      << "  1/(q.PT)" <<endreq;
  s << "----------------------------" << endreq;
  s    <<setw(20)<< traj.radius()
      << "    radius" <<endreq;
  s    <<setw(20)<< traj.curvature()
      << " curvature" <<endreq<<endreq;

  return s;
}

  // Explicit Constructor 
    TrackTrajectory::TrackTrajectory(TrackParameters p,
                                     ::HepPoint3D startPoint,
                                     double curvature
                                     ): 
    m_parameters(p), m_startPoint(startPoint), m_curvature(curvature)
  {  }
  
  // Copy Constructor
  TrackTrajectory::TrackTrajectory( const TrackTrajectory& other ) 
  {
    (*this) = other;
  }
  // Assignment 
  TrackTrajectory& TrackTrajectory::operator=( const TrackTrajectory& other ) {
    m_parameters      = other.m_parameters ;
    m_startPoint      = other.m_startPoint;
    m_curvature       = other.m_curvature;
    return *this ;
  } 
  
  // Destructor
  TrackTrajectory::~TrackTrajectory() {}   

  // Constructor from a momentum, vertex and charge
  TrackTrajectory::TrackTrajectory(  
                                   const Hep3Vector& momentum, 
                                   const Hep3Vector& vertex,
                                   double charge,
                                   double bField 
                                   )  
  {
    double pt = momentum.perp();
    double q  = charge;
    q /= abs(q); //  q = |1|
    
    // Calculate helix radius of curvature and centre

    // Working in GeV and cm 
    //double radius = pt/(2*0.00149898*bField);

    // Working in MeV and mm 
    double radius = pt/(2*0.149898*bField);
    double x0  = vertex.x() + q*radius*sin(momentum.phi());
    double y0  = vertex.y() - q*radius*cos(momentum.phi());
    Hep3Vector hCentre(x0,y0);
    
    // calculate parameters     
    // Tidier version of Fortran code zvhepa.F - Tarta
    
    double impactParameter = -q*(hCentre.perp() - radius);
    double theta = momentum.theta();
    double cotTheta = 1/tan(theta);
    double eta = momentum.pseudoRapidity();
    double zPerigee;
    double phi;
    
    //  Calculate Phi
    
    if (abs(vertex.x()) < 1e-4 && abs(vertex.y()) < 1e-4)
      {phi = momentum.phi();} else 
        {if ( (vertex.phi() - hCentre.phi()) > 0 )
          { phi =  hCentre.phi() + M_PI/2; } else   
            { phi = hCentre.phi() - M_PI/2;} }
    
    double x1 =  sin(phi)*q*hCentre.perp()/(radius*radius);
    double y1 = -cos(phi)*q*hCentre.perp()/(radius*radius);
    if (x1*x0 < 0 || y1*y0 < 0) {phi += M_PI;}
    if (phi < 0) {phi += 2.*M_PI;} else 
      { if(phi > 2.*M_PI) {phi -= 2.*M_PI;} }
    
    // z of perigee
    
    double dot = vertex.x()*momentum.x() + vertex.y()*momentum.y();
    double dotsign = dot/abs(dot) ;
    double dV = (vertex.perp()*vertex.perp()) - 
      (impactParameter*impactParameter);
    double denom = hCentre.perp()/radius;
    
    // NAN protection
    if (dV < 0) {dV = 0;}    
    if (abs(dotsign)!=1) dotsign=0;
    
    zPerigee = vertex.z() - (2*q*radius)*dotsign*cotTheta*
      asin((sqrt(dV/(denom)))*q/(2*radius) );
    
    // set parameters
    m_parameters = TrackParameters
      (eta, phi, pt, impactParameter, zPerigee, cotTheta, q/pt );
    m_startPoint = HepPoint3D(vertex);
    m_curvature = radius;  
  }

  // simple query services
  TrackParameters TrackTrajectory::parameters() const {
    return m_parameters;
  }
  HepPoint3D TrackTrajectory::startPoint() const
  { return m_startPoint ; }
  double TrackTrajectory::radius() const
  { return m_startPoint.perp() ; }
  
  double TrackTrajectory::curvature () const
  { return m_curvature ; }
  
  int    TrackTrajectory::signOfCharge() const              
  { return int( abs(m_parameters.invPtCharge() )/m_parameters.invPtCharge() ) ;  }
  
}  // End namespace

---