From 0624cd6bb42030ff4c13751dab0b5f3693ec633f Mon Sep 17 00:00:00 2001 From: Goetz Gaycken <Goetz.Gaycken@cern.ch> Date: Tue, 14 Jul 2015 15:45:19 +0200 Subject: [PATCH] Create tag xAODTracking-00-13-08-01 (xAODTracking-00-13-08-01) * Added missing inline keywords for TrackParticle helper methods. 2015-05-28 Goetz Gaycken * added helper to compute the d0 significance which takes the width of the beamspot into account. --- .../Root/TrackParticlexAODHelpers.cxx | 32 +++-- ...Tracking_TrackParticlexAODHelpers_test.cxx | 124 ++++++++++++++++++ .../xAODTracking/TrackParticlexAODHelpers.h | 82 ++++++++++-- 3 files changed, 219 insertions(+), 19 deletions(-) diff --git a/Event/xAOD/xAODTracking/Root/TrackParticlexAODHelpers.cxx b/Event/xAOD/xAODTracking/Root/TrackParticlexAODHelpers.cxx index c34ff8df059..4d9375aea40 100644 --- a/Event/xAOD/xAODTracking/Root/TrackParticlexAODHelpers.cxx +++ b/Event/xAOD/xAODTracking/Root/TrackParticlexAODHelpers.cxx @@ -17,7 +17,7 @@ namespace xAOD { static SG::AuxElement::Accessor< std::vector<float> > acc( "definingParametersCovMatrix" ); if( !acc.isAvailable( *tp ) ) { throw std::runtime_error("TrackParticle without covariance matrix for the defining parameters."); - } + } } } @@ -28,16 +28,28 @@ namespace xAOD { return true; } - double TrackingHelpers::d0significance(const xAOD::TrackParticle *tp) { - checkTPAndDefiningParamCov(tp); - double d0 = tp->d0(); - // elements in definingParametersCovMatrixVec should be : sigma_d0^2, sigma_d0_z0, sigma_z0^2 - double sigma_d0 = tp->definingParametersCovMatrixVec().at(0); - if (sigma_d0<=0.) { - throw std::runtime_error("TrackParticle with zero or negative d0 uncertainty."); + namespace TrackingHelpers { + inline double _d0significance(const xAOD::TrackParticle *tp, double d0_uncert_beam_spot_2) { + checkTPAndDefiningParamCov(tp); + double d0 = tp->d0(); + // elements in definingParametersCovMatrixVec should be : sigma_d0^2, sigma_d0_z0, sigma_z0^2 + double sigma_d0 = tp->definingParametersCovMatrixVec().at(0); + if (sigma_d0<=0.) { + throw std::runtime_error("TrackParticle with zero or negative d0 uncertainty."); + } + return d0/sqrt(sigma_d0+d0_uncert_beam_spot_2); } - return d0/sqrt(sigma_d0); + } + double TrackingHelpers::d0significance(const xAOD::TrackParticle *tp) { + return _d0significance(tp,0.); + } + + double TrackingHelpers::d0significance(const xAOD::TrackParticle *tp, double beam_sigma_x, double beam_sigma_y, double beam_sigma_xy) { + if (!checkBeamSpotSigma(beam_sigma_x, beam_sigma_y, beam_sigma_xy)) { + throw std::runtime_error("Beamspot covariance matrix is invalid."); + } + return _d0significance(tp,d0UncertaintyBeamSpot2(tp->phi(),beam_sigma_x, beam_sigma_y, beam_sigma_xy)); } double TrackingHelpers::z0significance(const xAOD::TrackParticle *tp, const xAOD::Vertex *vx) { @@ -63,7 +75,7 @@ namespace xAOD { if (sigma_z0<=0.) { throw std::runtime_error("TrackParticle with zero or negative z0 uncertainty."); } - return z0/sqrt(sigma_z0); + return z0/sqrt(sigma_z0); } double TrackingHelpers::pTErr2(const xAOD::TrackParticle *tp) { diff --git a/Event/xAOD/xAODTracking/test/xAODTracking_TrackParticlexAODHelpers_test.cxx b/Event/xAOD/xAODTracking/test/xAODTracking_TrackParticlexAODHelpers_test.cxx index b401c93c112..b296d1a9c27 100644 --- a/Event/xAOD/xAODTracking/test/xAODTracking_TrackParticlexAODHelpers_test.cxx +++ b/Event/xAOD/xAODTracking/test/xAODTracking_TrackParticlexAODHelpers_test.cxx @@ -112,6 +112,128 @@ void testPtErr() { } +namespace { + class Cov_t + { + public: + Cov_t(double cov_x, double cov_y, double cov_xy) { + m_cov[0]=cov_x; + m_cov[1]=cov_y; + m_cov[2]=cov_xy; + } + + double x() const { return m_cov[0]; } + double y() const { return m_cov[1]; } + double xy() const { return m_cov[2]; } + private: + double m_cov[3]; + }; + + template <class T> + bool checkEqual(T ref, T value ) { + if (std::abs(ref) * std::numeric_limits<T>::epsilon() < 4*std::numeric_limits<T>::min() ) { + return std::abs( value - ref ) < 4*std::numeric_limits<T>::epsilon(); + } + else { + return std::abs( value - ref ) < 4*std::numeric_limits<T>::epsilon() * ref; + } + } + + double toRad(double deg) { + return deg * M_PI/180; + } + + double sqr( double a) { + return a*a; + } +} + +void test_d0significance() { + xAOD::TrackParticleAuxContainer aux; + xAOD::TrackParticleContainer tpc; + tpc.setStore( &aux ); + + std::vector<Cov_t> beamspot_sigma_list; + std::vector<float> phi_list; + std::vector<bool> valid; + + phi_list.push_back(0.) ; + phi_list.push_back(toRad(60)) ; + phi_list.push_back(toRad(90)) ; + phi_list.push_back(toRad(210)) ; + phi_list.push_back(M_PI+acos(1/sqrt(2))) ; + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, 0.015*0.01 * 0.) ); + valid.push_back(true); + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, 0.015*0.01 ) ); + valid.push_back(true); + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, 0.015*0.01 *.6) ); + valid.push_back(true); + + beamspot_sigma_list.push_back( Cov_t(0.01, 0.015, -0.015*0.01 *.6) ); + valid.push_back(true); + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, -0.015*0.01) ); + valid.push_back(true); + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, -0.015*0.01*(1.+1e-6)) ); + valid.push_back(false); + + beamspot_sigma_list.push_back( Cov_t(0.015, 0.01, 0.015*0.01*(1.+1e-6)) ); + valid.push_back(false); + + float ref_d0=1; + float ref_z0=10; + float ref_d0_uncert=.05; + float ref_z0_uncert=.5; + float ref_qoverp = 1/30e3; + float ref_theta=60/180*M_PI; + std::vector<bool>::const_iterator valid_iter = valid.begin(); + for (const Cov_t &beamspot_sigma: beamspot_sigma_list) { + for( const float &phi: phi_list) { + assert( valid_iter != valid.end()); + xAOD::TrackParticle* p = new xAOD::TrackParticle(); + tpc.push_back( p ); + p->setDefiningParameters(ref_d0, ref_z0, phi, ref_theta, ref_qoverp); + xAOD::ParametersCovMatrix_t cov; + cov(0,0)=xAOD::TrackingHelpers::sqr(ref_d0_uncert); + cov(1,1)=xAOD::TrackingHelpers::sqr(ref_z0_uncert); + p->setDefiningParametersCovMatrix(cov); + + double vx_cov=beamspot_sigma.x(); + double vy_cov=beamspot_sigma.y(); + double vxy_cov=beamspot_sigma.xy(); + double expected_d0_uncert_2=sqr(ref_d0_uncert) + sin(phi)*(sin(phi)*sqr(vx_cov)-cos(phi)*vxy_cov)+cos(phi)*(cos(phi)*sqr(vy_cov)-sin(phi)*vxy_cov); + double _phi =p->phi(); + assert( checkEqual<float>(_phi, phi)); + assert( checkEqual<float>(ref_d0, p->d0())); + double _d0_uncert = p->definingParametersCovMatrixVec()[0]; + assert( checkEqual<float>(sqr(ref_d0_uncert), _d0_uncert)); + expected_d0_uncert_2=p->definingParametersCovMatrixVec()[0] + sin(_phi)*(sin(_phi)*sqr(vx_cov)-cos(_phi)*vxy_cov)+cos(_phi)*(cos(_phi)*sqr(vy_cov)-sin(_phi)*vxy_cov); + + assert( expected_d0_uncert_2 >= 0.); + double expected_d0_significance = ref_d0 / sqrt( expected_d0_uncert_2); + expected_d0_significance = p->d0() / sqrt( expected_d0_uncert_2); + + try { + if (*valid_iter) { + double d0_significance = xAOD::TrackingHelpers::d0significanceUnsafe(p, beamspot_sigma.x(), beamspot_sigma.y(), beamspot_sigma.xy()); + assert( checkEqual<float>( d0_significance, expected_d0_significance ) ); + } + double d0_significance = xAOD::TrackingHelpers::d0significance(p, beamspot_sigma.x(), beamspot_sigma.y(), beamspot_sigma.xy()); + assert( checkEqual<float>( d0_significance, expected_d0_significance ) ); + } + catch( std::exception &) { + assert( !( *valid_iter ) ); + } + } + assert( valid_iter != valid.end()); + ++valid_iter; + } +} + int main() { std::cout << "Run:" << __FILE__ << std::endl; // Create the main containers to test: @@ -269,6 +391,8 @@ int main() { testPtErr(); + test_d0significance(); + std::cout << "PASSED:" << __FILE__ << std::endl; return 0; } diff --git a/Event/xAOD/xAODTracking/xAODTracking/TrackParticlexAODHelpers.h b/Event/xAOD/xAODTracking/xAODTracking/TrackParticlexAODHelpers.h index 53da5cb3348..33055a0ca36 100644 --- a/Event/xAOD/xAODTracking/xAODTracking/TrackParticlexAODHelpers.h +++ b/Event/xAOD/xAODTracking/xAODTracking/TrackParticlexAODHelpers.h @@ -16,7 +16,10 @@ namespace xAOD { namespace TrackingHelpers { - + ///@brief convenience method to calculate the square of a value. + inline + double sqr(double a) { return a*a; } + ///@brief Get the impact parameter significance of a track particle in the r-phi plane. ///@return the impact parameter significance of the track particle. ///@throw this method may throw an exception in case the uncertainty is zero or the covariance matrix does not exist. @@ -28,6 +31,7 @@ namespace xAOD { ///In case the covariance matrix does not exist or the uncertainty is zero, this method will return an undefined value. ///In case the track particle does not exist the behaviour is undefined (segmentation fault or random return value). ///To validated the inputs use @ref hasValidCovD0, @ref hasValidCovD0andZ0, or @ref checkPVReference. + 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 @@ -35,6 +39,51 @@ namespace xAOD { return d0/sigma_d0; } + ///@brief calculate the squared d0 uncertainty component due to the size of the beam spot. + ///@param track_phi0 the phi angle of the track particle at the perigee wrt. the average beamspot position. + ///@param beam_sigma_x the width of the beamspot in the x-direction e.g. @ref IBeamCondSvc:: beamSigma(0). + ///@param beam_sigma_y the width of the beamspot in the y-direction e.g. @ref IBeamCondSvc:: beamSigma(1). + ///@param beam_sigma_xy the beamspot xy correlation e.g. @ref IBeamCondSvc:: beamSigmaXY. + ///@return the squared d0 uncertainty component due to the size of the beam spot. + 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; + } + + + ///@brief Get the impact parameter significance of a track particle in the r-phi plane where the d0 uncertainty takes the finite beamspot width into account. + ///@param beam_sigma_x the width of the beamspot in the x-direction e.g. @ref IBeamCondSvc:: beamSigma(0). + ///@param beam_sigma_y the width of the beamspot in the y-direction e.g. @ref IBeamCondSvc:: beamSigma(1). + ///@param beam_sigma_xy the beamspot xy correlation e.g. @ref IBeamCondSvc:: beamSigmaXY. + ///@return the impact parameter significance of the track particle where the d0 uncertainty takes the finite beamspot width into account. + ///@throw this method may throw an exception in case the uncertainty is zero or the covariance matrix does not exist or the beamspot covariance matrix is invalid. + ///The impact parameter and uncertainty are those stored in the track particle. Will perform input (tp, cov) validity checks. + double d0significance(const xAOD::TrackParticle *tp, double beam_sigma_x, double beam_sigma_y, double beam_sigma_xy); + + ///@brief Unsafe version of @ref d0significance with beam spot uncertainty. + ///@param tp pointer to a track particle. + ///@param beam_sigma_x the width of the beamspot in the x-direction e.g. @ref IBeamCondSvc:: beamSigma(0). + ///@param beam_sigma_y the width of the beamspot in the y-direction e.g. @ref IBeamCondSvc:: beamSigma(1). + ///@param beam_sigma_xy the beamspot xy correlation e.g. @ref IBeamCondSvc:: beamSigmaXY. + ///@return the impact parameter significance of the track particle where the d0 uncertainty takes the finite beamspot width into account, or an undefined value + ///In case the covariance matrix does not exist or the uncertainty is zero, this method will return an undefined value. + ///In case the track particle does not exist the behaviour is undefined (segmentation fault or random return value). + // 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. + ///To validated the inputs use @ref hasValidCovD0, @ref hasValidCovD0andZ0, or @ref checkPVReference or @ref checkBeamSpotSigma. + 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; + } + ///@brief Get the impact parameter significance of a track particle in the z direction. ///@param tp a pointer to a track particle. ///@param vx a pointer to a primary vertex with respect to which z0 is expressed or NULL. @@ -49,6 +98,7 @@ namespace xAOD { ///@return the impact parameter significance of the track particle or an undefined value ///In case the covariance matrix does not exist or the uncertainty is zero or smaller, this method will return an undefined value. ///In case the track particle does not exist the behaviour is undefined (segmentation fault or random return value). + 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 @@ -58,17 +108,18 @@ namespace xAOD { ///@brief Unsafe version of @ref z0significance, which uses z0 relative to a given primary vertex. ///@param tp a valid pointer to a track particle. - ///@param vx a valid pointer to a primary vertex with respect to which z0 is expressed. + ///@param vx a valid pointer to a primary vertex with respect to which z0 is expressed. ///@return the impact parameter significance of the track particle or an undefined value ///In case the covariance matrix does not exist or the uncertainty is zero or smaller, this method will return an undefined value. ///In case the track particle or vertex does not exist the behaviour is undefined (segmentation fault or random return value). ///To validated the inputs use @ref hasValidCovZ0 or @ref hasValidCovD0andZ0, and @ref checkPVReference. ///If the given vertex results from a fit which includes the track particle tp, this z0 impact parameter significance ///will be biased. + 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. + // assume that z0 was expressed relative to the associated vertex. if (tp->vertex()) { z0 += tp->vertex()->z(); } @@ -84,18 +135,20 @@ namespace xAOD { ///@brief Check whether the given track particle is valid and has a valid d0 uncertainty. ///@return true if the track particle is valid, has a covariance matrix, and a valid d0 uncertainty. + inline bool hasValidCovD0(const xAOD::TrackParticle *tp ) { if (hasValidCov(tp)) { if ( !tp->definingParametersCovMatrixVec().empty() && tp->definingParametersCovMatrixVec()[0]>0.) { return true; - } + } } return false; } ///@brief Check whether the given track particle is valid and has a valid z0 uncertainty. ///@return true if the track particle is valid, has a covariance matrix, and a valid z0 uncertainty. + inline bool hasValidCovZ0(const xAOD::TrackParticle *tp) { if (hasValidCov(tp)) { if ( ! ( tp->definingParametersCovMatrixVec().size() > 2 ) @@ -108,6 +161,7 @@ namespace xAOD { ///@brief Check whether the given track particle is valid and has a valid d0 and z0 uncertainty. ///@return true if the track particle is valid, has a covariance matrix, and a valid d0 and z0 uncertainty. + inline bool hasValidCovD0andZ0( const xAOD::TrackParticle *tp) { if (hasValidCov(tp)) { if ( ! ( tp->definingParametersCovMatrixVec().size() > 2 ) @@ -119,8 +173,6 @@ namespace xAOD { return false; } - ///@brief convenience method to calculate the square of a value. - double sqr(double a) { return a*a; } ///@brief test whether the given primary vertex has a significant displacement in r-phi wrt. d0 uncertainty. ///@param tp pointer to a track particle. @@ -129,6 +181,7 @@ namespace xAOD { ///@return true if the given primary vertex has no significant displacement in r-phi wrt. d0 uncertainty. ///The method will also return false if the track particle or vertex are invalid or if the track particle ///does not have a valid d0 uncertainty. + 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]; @@ -136,15 +189,25 @@ namespace xAOD { return false; } + ///@brief check that the beamspot covariance matrix is valid + ///@param beam_sigma_x the width of the beamspot in the x-direction e.g. @ref IBeamCondSvc:: beamSigma(0). + ///@param beam_sigma_y the width of the beamspot in the y-direction e.g. @ref IBeamCondSvc:: beamSigma(1). + ///@param beam_sigma_xy the beamspot xy correlation e.g. @ref IBeamCondSvc:: beamSigmaXY. + ///@return true if the beamspot covariance matrix is valid. + 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); + } ///@brief compute the uncertainty of pt squared. ///@param tp a pointer to a track particle. ///@throw will throw an exception if the track particle is not valid, no covariance matrix of the defining parameters is set or the covariance matrix has wrong dimension. double pTErr2(const xAOD::TrackParticle *tp); - + ///@brief compute the uncertainty of pt. ///@param tp a pointer to a track particle. ///@throw will throw an exception if the track particle is not valid, no covariance matrix of the defining parameters is set or the covariance matrix has wrong dimension. + inline double pTErr(const xAOD::TrackParticle *tp) { return sqrt( pTErr2(tp) ); } @@ -153,6 +216,7 @@ namespace xAOD { ///@brief compute the uncertainty of pt squared. ///@param tp a valid pointer to a track particle for which the defining parameters covariance matrix is set and valid ///undefined behaviour if tp is invalid or no valid covariance matrix is set for the defining parameters. + inline double pTErr2Unsafe(const xAOD::TrackParticle *tp) { // / d \2 / d \2 d d @@ -197,12 +261,14 @@ namespace xAOD { ///@brief compute the uncertainty of pt. ///@param tp a valid pointer to a track particle for which the defining parameters covariance matrix is set and valid ///undefined behaviour if tp is invalid or no valid covariance matrix is set for the defining parameters. + inline double pTErrUnsafe(const xAOD::TrackParticle *tp) { return sqrt(pTErr2Unsafe(tp)); } /// return true if the covariance matrix of the defining parameters is set, has enough elements and the q/p is valid /// @TODO also check theta ? + inline bool hasValidCovQoverP(const xAOD::TrackParticle *tp ) { if (hasValidCov(tp)) { if ( tp->definingParametersCovMatrixVec().size()>=15 && std::abs(tp->qOverP())>0.) { @@ -212,8 +278,6 @@ namespace xAOD { return false; } - - }// TrackParticleHelpers } // namespace xAOD -- GitLab