/cvmfs/atlas.cern.ch/repo/sw/ASG/AnalysisBase/2.4.31/xAODTracking/xAODTracking/TrackParticlexAODHelpers.h

// Dear emacs, this is -*- c++ -*-
// $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 definingParametersCovMatrixVec should be : sigma_d0^2, sigma_d0_z0, sigma_z0^2
      double sigma_d0 = std::sqrt( tp->definingParametersCovMatrixVec()[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();
      // elements in definingParametersCovMatrixVec should be : sigma_d0^2, sigma_d0_z0, sigma_z0^2
      double sigma_d0 = std::sqrt( tp->definingParametersCovMatrixVec()[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);

    inline
    double z0significanceUnsafe(const xAOD::TrackParticle *tp) {
      double z0 = tp->z0();
      // elements in definingParametersCovMatrixVec should be : sigma_z0^2, sigma_d0_z0, sigma_z0^2
      double sigma_z0 = std::sqrt( tp->definingParametersCovMatrixVec()[2] );
      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() - vx->z();
      // assume that z0 was expressed relative to the associated vertex.
      if (tp->vertex()) {
    z0 += tp->vertex()->z();
      }
      // elements in definingParametersCovMatrixVec should be : sigma_z0^2, sigma_d0_z0, sigma_z0^2
      double sigma_z0 = std::sqrt( tp->definingParametersCovMatrixVec()[2] );
      return z0/sigma_z0;
    }


    bool hasValidCov(const xAOD::TrackParticle *tp);

    inline
    bool hasValidCovD0(const xAOD::TrackParticle *tp ) {
          if (hasValidCov(tp)) {
         if  ( !tp->definingParametersCovMatrixVec().empty()
                   && tp->definingParametersCovMatrixVec()[0]>0.) {
                return true;
             }
          }
          return false;
    }

    inline
    bool hasValidCovZ0(const xAOD::TrackParticle *tp) {
          if (hasValidCov(tp)) {
             if  ( ! ( tp->definingParametersCovMatrixVec().size() > 2 )
                   && ( tp->definingParametersCovMatrixVec()[2] > 0. ) ) {
                return true;
             }
          }
          return false;
    }

    inline
    bool hasValidCovD0andZ0( const xAOD::TrackParticle *tp) {
          if (hasValidCov(tp)) {
             if  ( ! ( tp->definingParametersCovMatrixVec().size() > 2 ) 
                   && ( tp->definingParametersCovMatrixVec()[0] > 0. )
                   && ( tp->definingParametersCovMatrixVec()[2] > 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->definingParametersCovMatrixVec()[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 )

      
      const std::vector<float>&cov= tp->definingParametersCovMatrixVec();
      // elements in definingParametersCovMatrixVec 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[14] + diff_theta * cov[13] ) + sqr(diff_theta) * cov[9];
      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  ( tp->definingParametersCovMatrixVec().size()>=15 && std::abs(tp->qOverP())>0.) {
                return true;
             }
          }
          return false;
    }

  }// TrackParticleHelpers

} // namespace xAOD
  
#endif // XAOD_TRACKPARTICLEXAODHELPERS_H

---