/Users/jmonk/Physics/ForIA/src/Track.cxx

Go to the documentation of this file.

#include "ForIA/Track.hh"

#include <gsl/gsl_cdf.h>

#include <cmath>
#include <stdexcept>

namespace ForIA{
 
  Track::Track(): IFourMomentum(), m_gotChi2_dof(false), m_gotProbability(false), 
  m_author(BAD_AUTHOR),
  m_doE(true), m_isMC(false), 
  m_haveThetaOrigins(false), m_haveTV(false),
  m_haveMatchedVertices(false){
    
  }
  double Track::eta() const{
    return m_eta;
  }
  
  double Track::theta(Origin org) const{
    
    switch(org){
      case ORIGIN:
        return m_theta;
        break;
      case BEAM_SPOT:
        if(!m_haveThetaOrigins){
          throw std::runtime_error("Track: Do not have information on theta relative to beam spot!");
        }
        return m_theta_wrt_bs;
        break;
      case PRIMARY_VERTEX:
        
        if(m_haveThetaOrigins){
          return m_theta_wrt_pvtx;
        }else{
          return m_theta;
        }
        
        break;
        
      case TRACK_VERTEX:
        if(!m_haveTV){
          throw std::runtime_error("Track: Do not have information on associated vertex!");
        }
        return m_theta_wrt_tv;
        break;
    }
    
    return m_theta;
  }
  
  double Track::theta()const{
    return theta(PRIMARY_VERTEX);
  }
  
  double Track::phi(Origin org) const{
    
    switch(org){
      case ORIGIN:
        return m_phi;
        break;
      case BEAM_SPOT:
        return m_phi_wrt_bs;
        break;
      case PRIMARY_VERTEX:
        return m_phi_wrt_pvtx;
        break;
        
      case TRACK_VERTEX:
        if(!m_haveTV){
          throw std::runtime_error("Track: Do not have information on associated vertex!");
        }
        return m_phi_wrt_tv;
        break;

    }
    return m_phi;
  }
  
  double Track::phi()const{
    return phi(PRIMARY_VERTEX);
  }
  
  double Track::ET() const{
    return m_PT;
  }
  
  double Track::E() const{
    if(!m_doE) return m_e;
    m_e =  ET() / fabs(cos(theta()));
    m_doE = false;
    return m_e;
  }
  
  double Track::PT() const{
    return m_PT;
  }
  
  double Track::px(Origin org) const{
    return PT() * sin(phi(org));
  }
  
  double Track::py(Origin org)const{
    return PT() * cos(phi(org));
  }
  
  double Track::pz(Origin org) const{
    double tanp = tan(theta(org));
    if(fabs(tanp) < 1.e-12){
      std::cout<<"Warning: treating track as massless"<<std::endl;
      return E();
    }
    
    return PT() / tanp;
  }
  
  double Track::px()const{
    return px(PRIMARY_VERTEX);
  }
  
  double Track::py()const{
    return py(PRIMARY_VERTEX);
  }
  
  double Track::pz()const{
    return pz(PRIMARY_VERTEX);
  }
  double Track::d0(Origin org) const{
    
    switch(org){
      case ORIGIN:
        return m_d0;
        break;
      case BEAM_SPOT:
        return m_d0_wrt_bs;
        break;
      case PRIMARY_VERTEX:
        return m_d0_wrt_pvtx;
        break;
      case TRACK_VERTEX:
        if(!m_haveTV){
          throw std::runtime_error("Track: Do not have information on associated vertex!");
        }
        return m_d0_wrt_tv;
        break;
    }
    
    return m_d0;
  }
  
  double Track::z0(Origin org) const{
    
    switch(org){
      case ORIGIN:
        return m_z0;
        break;
      case BEAM_SPOT:
        return m_z0_wrt_bs;
        break;
      case PRIMARY_VERTEX:
        return m_z0_wrt_pvtx;
        break;
      case TRACK_VERTEX:
        if(!m_haveTV){
          throw std::runtime_error("Track: Do not have information on associated vertex!");
        }
        return m_z0_wrt_tv;
        break;
    }
    
    return m_z0;
  }
  
  double Track::chi2()const{
    
    // fudge factor of 1.3 to account between track hit error estimate differences
    // between MC and data
    if(isMC()) return m_chi2 * 1.3;
    
    return m_chi2;
  }
  
  int Track::DOF()const{
    return m_dof;
  }
  
  double Track::chi2_DOF()const{
    if(m_gotChi2_dof) return m_chi2_dof;
    m_chi2_dof = chi2() / ((double)m_dof);
    m_gotChi2_dof = true;
    return m_chi2_dof;
  }
  
  double Track::fitProbability()const{
    if(m_gotProbability) return m_probability;
    
    // see:
    //  http://www.gnu.org/software/gsl/manual/html_node/The-Chi_002dsquared-Distribution.html
    //  http://mathworld.wolfram.com/Chi-SquaredDistribution.html
    
    // this function is the integral of the chi^2 probability density above the given value
    // i.e. it is the probability thst the chi^2 is equal to or larger than the given value
    m_probability = gsl_cdf_chisq_Q(chi2(), (double)m_dof);
    
    m_gotProbability = true;
    return m_probability;
  }
  
  int Track::nPix() const{
    return m_nPix;
  }
  
  int Track::nSCT() const{
    return m_nSCT;
  }
  
  int Track::nTRT()const{
    return m_nTRT;
  }
  
  bool Track::expectBLayer()const{
    return m_expectBLayer;
  }
  
  int Track::nBLayer()const{
    return m_nBLayer;
  }
  
  Track::TrackAuthor Track::author()const{
    return m_author;
  }
  
  bool Track::isMC()const{
    return m_isMC;
  }
  
  int Track::id()const{
    return m_id;
  }
  
  VertexVector Track::associatedVertices()const{
    
    if(m_haveMatchedVertices) return VertexVector(m_vertices.begin(), m_vertices.end());
    m_event->matchVerticesToTracks();
    
    return VertexVector(m_vertices.begin(), m_vertices.end());
  }
  
  std::ostream &operator << (std::ostream &out, const Track &trk){
    out<<"Track: { PT: "<<trk.m_PT<<", eta: "<<trk.m_eta<<", phi: "<<trk.m_phi;
    out<<", d0: "<<trk.m_d0<<", z0: "<<trk.m_z0;
    out<<", #Pix: "<<trk.m_nPix<<", #SCT: "<<trk.m_nSCT<<"}";
    return out;
  }
}