TMBTrack.cpp

Go to the documentation of this file.

#include <cmath>
#include "tmb_tree/TMBTrack.hpp"
#include "tmb_tree/TMBIsoTrack.hpp"
#include "tmb_tree/TMBVertex.hpp"
#include "kinem_util/AnglesUtil.hpp"

using std::abs;
using std::log;
using std::sqrt;
using std::sinh;
using std::sin;
using std::cos;
using std::atan2;

// Local functions.

namespace {

  // Constants.

  Double_t CLIGHT = 2.99792458e10; // Speed of light in cm/sec.
  Double_t EQ = 1.0e-13 * CLIGHT;  // Electron charge in units GeV/(Tesla-cm).
  Double_t MPI = 0.13957;          // Charged pion mass.

}

ClassImp(TMBTrack);


//_____________________________________________________________________________

// Default constructor.

TMBTrack::TMBTrack() :
  _q(0),
  _chi2_ndof(-1.),
  _smtdedx(0.),
  _smtededx(0.),
  _bz(0.),
  _rdca(0.),
  _z(0.),
  _xdca(0.),
  _ydca(0.),
  _x_ps(0.),
  _y_ps(0.),
  _z_ps(0.),
  _px_ps(0.),
  _py_ps(0.),
  _pz_ps(0.),
  _trpars_valid(false),
  _trpos_valid(false)
{
  for(int i=0; i<3; ++i)
    _hitmask[i] = 0;
  for(int i=0; i<15; ++i)
    _trerrs[i] = 0.;
}

// Initializing constructor.

TMBTrack::TMBTrack(const Double_t* trpars, const Double_t* trerrs,
                   Double_t bz,
                   const UInt_t* hitmask, Double_t chi2_ndof, 
                   Double_t smtdedx, Double_t smt_ededx,
                   const Double_t* ps_position, const Double_t* ps_momentum,
                   const Double_t* surfpars,
                   const TRef& isotref) :
  TMBLorentzVector(trpars[4] != 0. ? 1./abs(trpars[4]) : 0., // pT
                   asinh(trpars[3]),                         // eta
                   trpars[2],                                // phi
                   MPI,                                      // mpi+
                   kPtEtaPhiM),
  _q(trpars[4] > 0. ? 1 : (trpars[4] < 0. ? -1 : 0)),
  _chi2_ndof(chi2_ndof),
  _smtdedx(smtdedx),
  _smtededx(smt_ededx),
  _isotref(isotref),
  _bz(bz),
  _rdca(trpars[0]),
  _z(trpars[1]),
  _xdca(surfpars[0]),
  _ydca(surfpars[1]),
  _x_ps(ps_position[0]),
  _y_ps(ps_position[1]),
  _z_ps(ps_position[2]),
  _px_ps(ps_momentum[0]),
  _py_ps(ps_momentum[1]),
  _pz_ps(ps_momentum[2]),
  _trpars_valid(true),
  _trpos_valid(false)
{
  for(int i=0; i<3; ++i)
    _hitmask[i] = hitmask[i];
  for(int i=0; i<5; ++i)
    _trpars[i] = trpars[i];
  for(int i=0; i<15; ++i)
    _trerrs[i] = trerrs[i];
}

// Destructor.

TMBTrack::~TMBTrack()
{}

// Clear (resets transient attributes).

void TMBTrack::Clear(Option_t* opt)
{
  _trpars_valid = false;
  _trpos_valid = false;
  TMBLorentzVector::Clear(opt);
}

Int_t TMBTrack::nhit() const
{
  Int_t Nhit = 0;
  for (Int_t mask=0; mask<3; ++mask){
    for (Int_t bit=0; bit<32; ++bit) {
      if (_hitmask[mask] & (1 << bit) ) ++Nhit;
    }
  }
  return Nhit;
}

Int_t TMBTrack::ncft() const
{
  Int_t Ncft = 0;
  for (Int_t bit=16; bit<32; ++bit) {
    if (_hitmask[0] & (1 << bit) ) ++Ncft;
  }
  return Ncft;
}

const TMBIsoTrack* TMBTrack::getIsoTrk() const 
{
  return _isotref.IsValid() ? dynamic_cast<const TMBIsoTrack*>(_isotref.GetObject()) : 0;
}

// detector eta of track for external CFT layer
// code from B. Tuchming (tuchming@cea.fr) 
Double_t TMBTrack::det_etaCFT() const
{
  double theta = kinem::theta(Eta());
  double pi = kinem::PI;
  
  // flag whether track is forward or backward
  int forward = 1;
  if(theta > pi / 2.) forward = -1;
  double z_pos = 126. - (forward*_z);
  double tantheta = fabs(tan(theta));
  double height = z_pos * tantheta;

  //flag is track is central or forward
  bool central = (height > 51.430);
  double theta_det=0.;

  if(central) {
    theta_det = fabs( atan2 (51.430, (51.430/tantheta + forward*_z)));
    if(forward==-1) theta_det = pi - theta_det;
  }
  else {
    theta_det = fabs( atan2 (((126. - (forward*_z)) * tantheta), 126.));
    if(forward==-1) theta_det = pi - theta_det;
  }

  return kinem::eta(theta_det);
}

void TMBTrack::fill_trpars() const
{
  if(_last_p != *this) {
    _last_p = *this;
    _trpars_valid = false;
    _trpos_valid = false;
  }
  if(_trpars_valid)
    return;
  _trpars_valid = true;
  _trpars[0] = _rdca;
  _trpars[1] = _z;
  _trpars[2] = Phi();
  Double_t c = CosTheta();
  _trpars[3] = c / sqrt(1. - c*c);
  _trpars[4] = _q / Pt();
}

void TMBTrack::fill_trpos() const
{
  if(_last_p != *this) {
    _last_p = *this;
    _trpars_valid = false;
    _trpos_valid = false;
  }
  if(_trpos_valid)
    return;
  _trpos_valid = true;
  _trpos[0] = _xdca + _rdca * Py() / Pt();   // x = x0 + r*sin(phid)
  _trpos[1] = _ydca - _rdca * Px() / Pt();   // y = y0 - r*cos(phid)
  _trpos[2] = _z;
}

// Analysis methods.

// Propagate methods.  These implementations use the exact formulas for
// updating the track parameters, but the error matrix is not modified,
// which is correct in leading order in the short-distance approximation.

void TMBTrack::propagate(const TMBVertex* vert) 
{
  propagate(vert->vx(), vert->vy());
}

void TMBTrack::propagate(Double_t xdca1, Double_t ydca1)
{
  // Get original trf track parameters.  Tan(lambda) and q/pT are constants.

  Double_t r0 = rdca();
  Double_t z0 = z();
  Double_t phid0 = phid();
  Double_t tanlm = tlm();
  Double_t qpt0 = qpt();

  // Transverse curvature (1/radius).

  Double_t wt = EQ * qpt0 * _bz;
  assert(wt != 0.);

  // Difference between final and initial surface parameters 
  // (final minus initial).

  Double_t dx = xdca1 - xdca();
  Double_t dy = ydca1 - ydca();

  // Sine and cosine of the initial phi direction.

  Double_t sphid0 = sin(phid0);              // sin(phid0)
  Double_t cphid0 = cos(phid0);              // cos(phid0)

  // Sine and cosine of the bend angle.

  Double_t dxi = dx * cphid0 + dy * sphid0;
  Double_t drho = dx * sphid0 - dy * cphid0;
  Double_t sdphi = wt*dxi;                   // sin(dphi)
  Double_t cdphi = 1. - wt * (r0 - drho);    // cos(dphi)
  Double_t norm = sqrt(sdphi*sdphi + cdphi*cdphi);
  sdphi /= norm;
  cdphi /= norm;
  Double_t dphi = atan2(sdphi, cdphi);

  // Sine and cosine of the final phi direction.

  // Double_t sdphid1 = sphid0 * cdphi + cphid0 * sdphi;   // sin(phid0 + dphi)
  // Double_t cdphid1 = cphid0 * cdphi - sphid0 * sdphi;   // cos(phid0 + dphi)

  // Transverse path distance to destination dca.

  Double_t st = dphi / wt;

  // Calculate final trf track parameters.

  Double_t r1 = (r0 - drho - (sdphi*sdphi / (wt * (1. + cdphi)))) / cdphi;

  Double_t sinlm = tanlm / sqrt(1. + tanlm*tanlm);   // sin(lambda)
  Double_t z1 = z0 + sinlm * st;
  Double_t phid1 = phid0 + dphi;

  // Update class attributes.

  static_cast<TMBLorentzVector>(*this) = 
    TMBLorentzVector(Pt(), Eta(), phid1, M(), kPtEtaPhiM);
  _rdca = r1;
  _z = z1;
  _xdca = xdca1;
  _ydca = ydca1;
  assert(_trpars_valid);
  _trpars[0] = r1;
  _trpars[1] = z1;
  _trpars[2] = phid1;
}

// Impact parameter methods.

void TMBTrack::impact(const TMBVertex* vert, Double_t* ip, 
                      Double_t* iperr) const 
{
    impact(vert->vx(), vert->vy(), vert->vz(), ip, iperr);
}

void TMBTrack::impact(Double_t x, Double_t y, Double_t z, Double_t* ip, 
                      Double_t* iperr) const
{
  // Make a non-const copy of track, then call non-const method.

  TMBTrack track(*this);
  track.impact(x, y, z, ip, iperr);
}

void TMBTrack::impact(const TMBVertex* vert, Double_t* ip, Double_t* iperr)
{
    impact(vert->vx(), vert->vy(), vert->vz(), ip, iperr);
}

void TMBTrack::impact(Double_t x, Double_t y, Double_t z, Double_t* ip, 
                      Double_t* iperr)
{
  // Propagate to vertex dca surface.

  propagate(x, y);

  // Store result.

  ip[0] = _rdca;
  ip[1] = _z - z;
  if(iperr) {
    iperr[0] = _trerrs[0];  // (0,0)
    iperr[1] = _trerrs[1];  // (1,0)
    iperr[2] = _trerrs[2];  // (1,1)
  }
}

void TMBTrack::constrain(const TMBVertex* vert)
{
  // First constrain this track without vertex error.

  constrain(vert->vx(), vert->vy());
}

void TMBTrack::constrain(Double_t x, Double_t y)
{
  // First propagate to the beam dca surface.

  propagate(x, y);
  assert(false);
}

---