TrackParticlexAODHelpers.h

Go to the documentation of this file.

// Dear emacs, this is -*- c++ -*-
 
/*
  Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
*/
 
// $Id: $
#ifndef XAOD_TRACKPARTICLEXAODHELPERS_H
#define XAOD_TRACKPARTICLEXAODHELPERS_H
 
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/Vertex.h"
#include <cmath>
 
 
namespace xAOD {
 
  namespace TrackingHelpers {
    inline
    double sqr(double a) { return a*a; }
 
    double d0significance(const xAOD::TrackParticle *tp);
 
    inline
    double d0significanceUnsafe(const xAOD::TrackParticle *tp) {
      double d0 = tp->d0();
      // elements in definingParametersCovMatrixDiagVec should be : sigma_d0^2, sigma_z0^2
      double sigma_d0 = std::sqrt( tp->definingParametersCovMatrixDiagVec()[0] );
      return d0/sigma_d0;
    }
 
    inline 
    double d0UncertaintyBeamSpot2(double track_phi0,  double beam_sigma_x,double beam_sigma_y, double beam_sigma_xy) {
      double sin_phi = sin(track_phi0);
      double cos_phi = cos(track_phi0);
      double d0_uncert2= sin_phi * ( sin_phi * sqr(beam_sigma_x)
                                    -cos_phi * beam_sigma_xy)
                        +cos_phi * ( cos_phi * sqr(beam_sigma_y)
                                    -sin_phi * beam_sigma_xy);
      return d0_uncert2;
    }
 
 
    double d0significance(const xAOD::TrackParticle *tp, double beam_sigma_x, double beam_sigma_y, double beam_sigma_xy);
 
    // In case beam_sigma_xy is larger than the square of beam_sigma_x and the square of beam_sigma_y the result is undefined.
    inline
    double d0significanceUnsafe(const xAOD::TrackParticle *tp, double beam_sigma_x,double beam_sigma_y, double beam_sigma_xy) {
      double d0 = tp->d0();
      double sigma_d0 = std::sqrt( tp->definingParametersCovMatrixDiagVec()[0] + d0UncertaintyBeamSpot2(tp->phi(),beam_sigma_x, beam_sigma_y, beam_sigma_xy) );
      return d0/sigma_d0;
    }
 
    double z0significance(const xAOD::TrackParticle *tp, const xAOD::Vertex *vx=NULL);
 
    double z0sinthetasignificance(const xAOD::TrackParticle *tp, const xAOD::Vertex *vx=NULL);
 
    inline
    double z0significanceUnsafe(const xAOD::TrackParticle *tp) {
      double z0 = tp->z0() + tp->vz();
      double sigma_z0 = std::sqrt( tp->definingParametersCovMatrixDiagVec()[1] );
      return z0/sigma_z0;
    }
 
    inline
    double z0significanceUnsafe(const xAOD::TrackParticle *tp, const xAOD::Vertex *vx) {
      // use z0 relative to the given primary vertex.
      double z0 = tp->z0() + tp->vz() - vx->z();
 
      double sigma_z0 = std::sqrt( tp->definingParametersCovMatrixDiagVec()[1] );
      return z0/sigma_z0;
    }
 
 
    bool hasValidCov(const xAOD::TrackParticle *tp);
 
    inline
    bool hasValidCovD0(const xAOD::TrackParticle *tp ) {
          if (hasValidCov(tp)) {
         if  (tp->definingParametersCovMatrixDiagVec()[0]>0.) {
                return true;
             }
          }
          return false;
    }
 
    inline
    bool hasValidCovZ0(const xAOD::TrackParticle *tp) {
          if (hasValidCov(tp)) {
             if  ( tp->definingParametersCovMatrixDiagVec()[1] > 0. ) {
                return true;
             }
          }
          return false;
    }
 
    inline
    bool hasValidCovD0andZ0( const xAOD::TrackParticle *tp) {
          if (hasValidCov(tp)) {
             if  ( ( tp->definingParametersCovMatrixDiagVec()[0] > 0. )
                   && ( tp->definingParametersCovMatrixDiagVec()[1] > 0. ) ) {
                return true;
             }
          }
          return false;
    }
 
 
    inline
    bool checkPVReference(const xAOD::TrackParticle *tp, const xAOD::Vertex *vx, const double max_pv_dxy_sqr=0.5*0.5) {
        if (hasValidCovD0(tp) && vx) {
         return std::abs( sqr(vx->x())+ sqr(vx->y())) <  max_pv_dxy_sqr * tp->definingParametersCovMatrixDiagVec()[0];
        }
        return false;
    }
 
    inline
    bool checkBeamSpotSigma(double beam_sigma_x,double beam_sigma_y, double beam_sigma_xy) {
      return sqr(beam_sigma_x)+sqr(beam_sigma_y)>=2*std::abs(beam_sigma_xy);
    }
 
    double pTErr2(const xAOD::TrackParticle *tp);
 
    inline
    double pTErr(const xAOD::TrackParticle *tp) {
      return sqrt( pTErr2(tp) );
    }
 
 
    inline
    double pTErr2Unsafe(const xAOD::TrackParticle *tp) {
 
      //            /   d                       \2    /   d                        \2          d       d
      // pt_err^2 =|  ------   (pt) *sigma_q/p  |  + | --------   pt * sigma_theta |    +    ----- pt ------- pt * sigma_theta_qp
      //            \  d q/p                    /     \ d theta                    /          d qp    d theta
 
      //                d             / d pt           2   d pt                     \   / d pt                  \2
      //          =   ------   (pt) * | ----- sigma_q/p  + ------- * sigma_theta_qp | + | ------- * sigma_theta |
      //               d q/p          \ d q/p              d theta                  /   \ d theta               /
 
      //   d                pt
      // ------ pt  =  - -------------
      //  d q/p          fabs(qOverP)
 
      //   d                pt
      //  ------- pt =   -----------
      //  d theta         tan(theta)
 
      double pt = tp->pt();
      double diff_qp =  - pt / std::abs(tp->qOverP());
 
      double diff_theta = pt / tan( tp->theta() );
 
      // since 1/tan (thata) defined for all floats 0..M_PI, don't have to protect :
      //
      // double diff_theta = ( std::abs(tp->theta()-M_PI_2) > std::numeric_limits<float>::epsilon() * 2
      //                      ? pt / tan( tp->theta() )
      //                      : 0 )
 
      auto cov = tp->definingParametersCovMatrix();
      // elements in definingParametersCovMatrix should be : 0: sigma_d0^2,
      //                                                     1: sigma_d0_z0,  sigma_z0^2,
      //                                                     3: sigma_d0_phi, sigma_z0_phi, sigma_phi^2
      //                                                     6: sigma_d0_th,  sigma_z0_th,  sigma_phi_th, sigma_th^2
      //                                                    10: sigma_d0_qp,  sigma_z0_qp,  sigma_phi_qp, sigma_th_qp, sigma_qp^2
 
      double pt_err2 = diff_qp * (diff_qp * cov(4, 4) + diff_theta * cov(3, 4) ) + sqr(diff_theta) * cov(3, 3);
      return pt_err2;
    }
 
    inline
    double pTErrUnsafe(const xAOD::TrackParticle *tp) {
      return sqrt(pTErr2Unsafe(tp));
    }
 
    inline
    bool hasValidCovQoverP(const xAOD::TrackParticle *tp ) {
          if (hasValidCov(tp)) {
             if  (std::abs(tp->qOverP())>0.) {
                return true;
             }
          }
          return false;
    }
 
  }// TrackParticleHelpers
 
} // namespace xAOD
  
#endif // XAOD_TRACKPARTICLEXAODHELPERS_H