From 5ae0739bed301d6bf5b48702dde28ef91bcd84c4 Mon Sep 17 00:00:00 2001 From: Wojciech Krzemien <wojciech.krzemien@ncbj.gov.pl> Date: Tue, 15 Oct 2019 15:12:21 +0200 Subject: [PATCH] Add updateCovariance and updatePosition methods for hypo Also, remove some unecessary comments and code --- .../src/CaloFutureSCorrection.cpp | 210 +++++++++--------- .../src/CaloFutureSCorrection.h | 3 + 2 files changed, 103 insertions(+), 110 deletions(-) diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp index 8f8774f07..af223589a 100644 --- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp +++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp @@ -54,13 +54,11 @@ CaloFutureSCorrection::CaloFutureSCorrection( const std::string& type, const std StatusCode CaloFutureSCorrection::finalize() { m_hypos.clear(); - // finalize the base class return CaloFutureCorrectionBase::finalize(); } // ============================================================================ StatusCode CaloFutureSCorrection::initialize() { - // first initialize the base class StatusCode sc = CaloFutureCorrectionBase::initialize(); if ( sc.isFailure() ) { return Error( "Unable initialize the base class CaloFutureCorrectionBase!", sc ); } if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) debug() << "Condition name : " << m_conditionName << endmsg; @@ -78,11 +76,9 @@ StatusCode CaloFutureSCorrection::correct( } for ( auto* hypo : hypos ) { - // check the Hypo auto h = std::find( m_hypos.begin(), m_hypos.end(), hypo->hypothesis() ); if ( m_hypos.end() == h ) return Error( "Invalid hypothesis!", StatusCode::SUCCESS ); - // No correction for negative energy : if ( hypo->e() < 0. ) { ++m_counterSkipNegativeEnergyCorrection; continue; @@ -113,13 +109,6 @@ StatusCode CaloFutureSCorrection::correct( continue; } - // int area = cellID.area(); // currently not used - - // Somewhat inelegant way of data sharing between this const method process() and calcSCorrection(). Use of private - // fields of this CaloFutureSCorrection class for the same purpose would break the constness of the process() - // interface. - // Note, that calling private calcSCorrection(..., &results) is thread-safe b/c that's a ptr to a stack variable. - const LHCb::CaloPosition& position = MainCluster->position(); const double xBar = position.x(); const double yBar = position.y(); @@ -128,8 +117,6 @@ StatusCode CaloFutureSCorrection::correct( auto seedPos = m_det->cellCenter( cellID ); SCorrInputParams params{cellID, seedPos, xBar, yBar, z}; - // SCorrOutputParams results = {0, 0}; // output parameters: just d(Xhypo)/d(Xcluster) and d(Yhypo)/d(Ycluster) - /** here all information is available * * ( ) Ecal energy in 3x3 : ( not used ) @@ -141,30 +128,28 @@ StatusCode CaloFutureSCorrection::correct( * (7) Position of seed cell : seedPos */ - // passing OUTPUT parameter &results as ptr to calcSCorrection() is thread-safe b/c that's a local stack variable - // calcSCorrection( xBar, yBar, xCor, yCor, params, - //&results ); // non-nullptr results means calculate and update results - struct SCorrResults results = getSCorrections( params ); double xCor = results.xCor; double yCor = results.yCor; - ; // corrected hypo position - const double& dXhy_dXcl = results.dXhy_dXcl; - const double& dYhy_dYcl = results.dYhy_dYcl; + // corrected hypo position + double dXhy_dXcl = results.dXhy_dXcl; + double dYhy_dYcl = results.dYhy_dYcl; // protection against unphysical d(Xhypo)/d(Xcluster) == 0 or d(Yhypo)/d(Ycluster) == 0 if ( fabs( dXhy_dXcl ) < 1e-10 ) { warning() << "unphysical d(Xhypo)/d(Xcluster) = " << dXhy_dXcl << " reset to 1 as if Xhypo = Xcluster" << endmsg; - const_cast<double&>( dXhy_dXcl ) = 1.; + dXhy_dXcl = 1.; } if ( fabs( dYhy_dYcl ) < 1e-10 ) { warning() << "unphysical d(Yhypo)/d(Ycluster) = " << dYhy_dYcl << " reset to 1 as if Yhypo = Ycluster" << endmsg; - const_cast<double&>( dYhy_dYcl ) = 1.; + dYhy_dYcl = 1.; } // numeric partial derivatives w.r.t. X and Y, necessary to check after any change to the S-corrections if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) && m_correctCovariance ) { - debugDerivativesCalculation( xBar, yBar, xCor, yCor, dXhy_dXcl, dYhy_dYcl, params ); + const double& dXhy_dXcl_test = dXhy_dXcl; + const double& dYhy_dYcl_test = dYhy_dYcl; + debugDerivativesCalculation( xBar, yBar, xCor, yCor, dXhy_dXcl_test, dYhy_dYcl_test, params ); } const LHCb::CaloPosition* pos = hypo->position(); @@ -179,70 +164,59 @@ StatusCode CaloFutureSCorrection::correct( debug() << "xcel/ycel " << params.seedPos.x() << "/" << params.seedPos.y() << endmsg; } - // update position - LHCb::CaloPosition::Parameters& parameters = hypo->position()->parameters(); - parameters( LHCb::CaloPosition::Index::X ) = xCor; - parameters( LHCb::CaloPosition::Index::Y ) = yCor; + updatePosition(xCor, yCor, hypo); m_counterDeltaX += xCor - xBar; m_counterDeltaY += yCor - yBar; - // update cov.m.: error propagation due to the S-correction - if ( m_correctCovariance ) { - LHCb::CaloPosition::Covariance& covariance = hypo->position()->covariance(); - - if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "before s-cor cov.m. = \n" << covariance << endmsg; } - - // cov.m packing in double array[5] following ROOT::Math::SMatrix<double,3,3>::Array() - // for row/column indices (X:0, Y:1, E:2), see comments in CaloFutureECorrection::process() - double c1[6]; - - c1[0] = covariance( LHCb::CaloPosition::Index::X, - LHCb::CaloPosition::Index::X ); // arr[0] not relying on LHCb::CaloPosition::Index::X == 0 - c1[2] = covariance( LHCb::CaloPosition::Index::Y, - LHCb::CaloPosition::Index::Y ); // arr[2] not relying on LHCb::CaloPosition::Index::Y == 1 - c1[5] = covariance( LHCb::CaloPosition::Index::E, - LHCb::CaloPosition::Index::E ); // arr[5] not relying on LHCb::CaloPosition::Index::E == 2 - c1[1] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ); // arr[1] - c1[3] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ); // arr[3] - c1[4] = covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ); // arr[4] - - // cov1 = (J * cov0 * J^T) for the special case of diagonal Jacobian for (X,Y,E) -> (X1=X1(X), Y1=Y1(Y), E1=E) - c1[0] *= dXhy_dXcl * dXhy_dXcl; - c1[1] *= dXhy_dXcl * dYhy_dYcl; - c1[2] *= dYhy_dYcl * dYhy_dYcl; - c1[3] *= dXhy_dXcl; - c1[4] *= dYhy_dYcl; - // c1[5] remains unchanged (energy is not chaged by S-correction) - - // alternatively, a code fragment for a general-form Jacobian (cf. a similar comment in - // CaloFutureECorrection::process()) TMatrixD jac(3, 3); // just a diagonal Jacobian in case of (X,Y,E) -> (X1(X), - // Y1(Y), E) transformation jac(0,0) = dXhy_dXcl; jac(1,1) = dYhy_dYcl; jac(2,2) = 1.; if ( msgLevel( MSG::DEBUG) - // ){ debug() << "s-cor jacobian = " << endmsg; jac.Print(); } TMarixDSym cov0(3) = ... // to be - // initilized from hypo->position()->covariance() TMarixDSym cov1(3); // resulting extrapolated - // cov.m. recalculate_cov(jac, cov0, cov1); // calculate: cov1 = (J * cov0 * J^T) - - // finally update CaloHypo::position()->covariance() - covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = c1[0]; // cov1(0,0); - covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = c1[2]; // cov1(1,1); - covariance( LHCb::CaloPosition::Index::E, LHCb::CaloPosition::Index::E ) = c1[5]; // cov1(2,2); - covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ) = c1[1]; // cov1(0,1); - covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ) = c1[3]; // cov1(0,2); - covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ) = c1[4]; // cov1(1,2); - - // // DG: my little paranoia, should be always ok since Covariance is SMatrix<3,3,double> internally represented - // as double array[5] assert( covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::Y) - // == covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::X)); assert( - // covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::E) == - // covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::X)); assert( - // covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::E) == - // covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::Y)); - - if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "after s-cor cov.m. = \n" << covariance << endmsg; } + if ( m_correctCovariance) { + updateCovariance(dXhy_dXcl, dYhy_dYcl, hypo); } } return StatusCode::SUCCESS; } +void CaloFutureSCorrection::updatePosition(double xCor, double yCor, LHCb::CaloHypo* hypo) const { + LHCb::CaloPosition::Parameters& parameters = hypo->position()->parameters(); + parameters( LHCb::CaloPosition::Index::X ) = xCor; + parameters( LHCb::CaloPosition::Index::Y ) = yCor; +} +void CaloFutureSCorrection::updateCovariance(double dXhy_dXcl, double dYhy_dYcl, LHCb::CaloHypo* hypo) const { + LHCb::CaloPosition::Covariance& covariance = hypo->position()->covariance(); + + if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "before s-cor cov.m. = \n" << covariance << endmsg; } + + // cov.m packing in double array[5] following ROOT::Math::SMatrix<double,3,3>::Array() + // for row/column indices (X:0, Y:1, E:2), see comments in CaloFutureECorrection::process() + double c1[6]; + + c1[0] = covariance( LHCb::CaloPosition::Index::X, + LHCb::CaloPosition::Index::X ); // arr[0] not relying on LHCb::CaloPosition::Index::X == 0 + c1[2] = covariance( LHCb::CaloPosition::Index::Y, + LHCb::CaloPosition::Index::Y ); // arr[2] not relying on LHCb::CaloPosition::Index::Y == 1 + c1[5] = covariance( LHCb::CaloPosition::Index::E, + LHCb::CaloPosition::Index::E ); // arr[5] not relying on LHCb::CaloPosition::Index::E == 2 + c1[1] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ); // arr[1] + c1[3] = covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ); // arr[3] + c1[4] = covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ); // arr[4] + + // cov1 = (J * cov0 * J^T) for the special case of diagonal Jacobian for (X,Y,E) -> (X1=X1(X), Y1=Y1(Y), E1=E) + c1[0] *= dXhy_dXcl * dXhy_dXcl; + c1[1] *= dXhy_dXcl * dYhy_dYcl; + c1[2] *= dYhy_dYcl * dYhy_dYcl; + c1[3] *= dXhy_dXcl; + c1[4] *= dYhy_dYcl; + // c1[5] remains unchanged (energy is not changed by S-correction) + + covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::X ) = c1[0]; // cov1(0,0); + covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::Y ) = c1[2]; // cov1(1,1); + covariance( LHCb::CaloPosition::Index::E, LHCb::CaloPosition::Index::E ) = c1[5]; // cov1(2,2); + covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::Y ) = c1[1]; // cov1(0,1); + covariance( LHCb::CaloPosition::Index::X, LHCb::CaloPosition::Index::E ) = c1[3]; // cov1(0,2); + covariance( LHCb::CaloPosition::Index::Y, LHCb::CaloPosition::Index::E ) = c1[4]; // cov1(1,2); + + if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { debug() << "after s-cor cov.m. = \n" << covariance << endmsg; } +} + // ============================================================================ struct CaloFutureSCorrection::SCorrResults @@ -348,37 +322,6 @@ CaloFutureSCorrection::getSCorrections( const struct SCorrInputParams& params ) */ if ( m_correctCovariance ) { - // // ---- calculation of numeric derivatives of individual correction functions, important for debugging in case of - // code changes --- debug() << "---------- numeric derivatives of individual S-correction functions ---------------" - // << endmsg; double tmpd = ( fabs(Asx0) > 1.e-5 ) ? Asx0*2.e-2 : 2.e-7; double dn_shapeX = ( - // getCorrection(CaloFutureCorrection::shapeX, cellID, Asx0 + tmpd, Asx0 + tmpd) - Asx1 )/tmpd; tmpd = ( fabs(Asy0) - // > 1.e-5 ) ? Asy0*2.e-2 : 2.e-7; double dn_shapeY = ( getCorrection(CaloFutureCorrection::shapeY, cellID, Asy0 - // + tmpd, Asy0 + tmpd) - Asy1 )/tmpd; - // - // double dn_angularX = ( getCorrection(CaloFutureCorrection::angularX, cellID, thx*1.002, 0.) - daX )/thx/2e-3; - // double dn_angularY = ( getCorrection(CaloFutureCorrection::angularY, cellID, thy*1.002, 0.) - daY )/thy/2e-3; - // - // tmpd = ( fabs(Asx2) > 1.e-5 ) ? Asx2*2.e-3 : 2.e-8; - // double dn_residualX = ( getCorrection((residualX_flag ? CaloFutureCorrection::residualX : - // CaloFutureCorrection::residual), - // cellID, Asx2 + tmpd, 0.) - // - dcX )/tmpd; - // tmpd = ( fabs(Asy2) > 1.e-5 ) ? Asy2*2.e-3 : 2.e-8; - // double dn_residualY = ( getCorrection((residualY_flag ? CaloFutureCorrection::residualY : - // CaloFutureCorrection::residual), - // cellID, Asy2 + tmpd, 0.) - // - dcY )/tmpd; - // tmpd = ( fabs(Asx3) > 1.e-5 ) ? Asx3*2.e-3 : 2.e-8; - // double dn_asymX = (xBar < 0 ) ? - // ( getCorrection(CaloFutureCorrection::asymM , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd : - // ( getCorrection(CaloFutureCorrection::asymP , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd ; - // - // tmpd = ( fabs(Asy3) > 1.e-5 ) ? Asy3*2.e-3 : 2.e-8; - // double dn_asymY = (yBar < 0 ) ? - // ( getCorrection(CaloFutureCorrection::asymM , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd : - // ( getCorrection(CaloFutureCorrection::asymP , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd ; - // // - // ------------------------------------------------------------------------------------------------------------------------------- // calculation of the analytic derivatives: // NB: printouts comparing analytic calculations with numeric derivatives which are commented-out below @@ -610,3 +553,50 @@ void CaloFutureSCorrection::debugDerivativesCalculation( const double& xBar, con return; } + +// WK: This was part of the comment - maybe it should be moved as kind of unit tests? +// // ---- calculation of numeric derivatives of individual correction functions, important for debugging in case of +// code changes --- debug() << "---------- numeric derivatives of individual S-correction functions ---------------" +// << endmsg; double tmpd = ( fabs(Asx0) > 1.e-5 ) ? Asx0*2.e-2 : 2.e-7; double dn_shapeX = ( +// getCorrection(CaloFutureCorrection::shapeX, cellID, Asx0 + tmpd, Asx0 + tmpd) - Asx1 )/tmpd; tmpd = ( fabs(Asy0) +// > 1.e-5 ) ? Asy0*2.e-2 : 2.e-7; double dn_shapeY = ( getCorrection(CaloFutureCorrection::shapeY, cellID, Asy0 +// + tmpd, Asy0 + tmpd) - Asy1 )/tmpd; +// +// double dn_angularX = ( getCorrection(CaloFutureCorrection::angularX, cellID, thx*1.002, 0.) - daX )/thx/2e-3; +// double dn_angularY = ( getCorrection(CaloFutureCorrection::angularY, cellID, thy*1.002, 0.) - daY )/thy/2e-3; +// +// tmpd = ( fabs(Asx2) > 1.e-5 ) ? Asx2*2.e-3 : 2.e-8; +// double dn_residualX = ( getCorrection((residualX_flag ? CaloFutureCorrection::residualX : +// CaloFutureCorrection::residual), +// cellID, Asx2 + tmpd, 0.) +// - dcX )/tmpd; +// tmpd = ( fabs(Asy2) > 1.e-5 ) ? Asy2*2.e-3 : 2.e-8; +// double dn_residualY = ( getCorrection((residualY_flag ? CaloFutureCorrection::residualY : +// CaloFutureCorrection::residual), +// cellID, Asy2 + tmpd, 0.) +// - dcY )/tmpd; +// tmpd = ( fabs(Asx3) > 1.e-5 ) ? Asx3*2.e-3 : 2.e-8; +// double dn_asymX = (xBar < 0 ) ? +// ( getCorrection(CaloFutureCorrection::asymM , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd : +// ( getCorrection(CaloFutureCorrection::asymP , cellID , Asx2 + tmpd , 0.) - ddcX )/tmpd ; +// +// tmpd = ( fabs(Asy3) > 1.e-5 ) ? Asy3*2.e-3 : 2.e-8; +// double dn_asymY = (yBar < 0 ) ? +// ( getCorrection(CaloFutureCorrection::asymM , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd : +// ( getCorrection(CaloFutureCorrection::asymP , cellID , Asy2 + tmpd , 0.) - ddcY )/tmpd ; +// // +// ------------------------------------------------------------------------------------------------------------------------------- +// +// // DG: my little paranoia, should be always ok since Covariance is SMatrix<3,3,double> internally represented +// as double array[5] assert( covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::Y) +// == covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::X)); assert( +// covariance(LHCb::CaloFuturePosition::Index::X, LHCb::CaloFuturePosition::Index::E) == +// covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::X)); assert( +// covariance(LHCb::CaloFuturePosition::Index::Y, LHCb::CaloFuturePosition::Index::E) == +// covariance(LHCb::CaloFuturePosition::Index::E, LHCb::CaloFuturePosition::Index::Y)); +// alternatively, a code fragment for a general-form Jacobian (cf. a similar comment in +// CaloFutureECorrection::process()) TMatrixD jac(3, 3); // just a diagonal Jacobian in case of (X,Y,E) -> (X1(X), +// Y1(Y), E) transformation jac(0,0) = dXhy_dXcl; jac(1,1) = dYhy_dYcl; jac(2,2) = 1.; if ( msgLevel( MSG::DEBUG) +// ){ debug() << "s-cor jacobian = " << endmsg; jac.Print(); } TMarixDSym cov0(3) = ... // to be +// initilized from hypo->position()->covariance() TMarixDSym cov1(3); // resulting extrapolated +// cov.m. recalculate_cov(jac, cov0, cov1); // calculate: cov1 = (J * cov0 * J^T) diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h index 85ee7f3bc..2616ecc44 100644 --- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h +++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h @@ -77,6 +77,9 @@ private: struct SCorrResults getSCorrections( const struct SCorrInputParams& params ) const; + void updateCovariance(double dXhy_dXcl, double dYhy_dYcl, LHCb::CaloHypo* hypo) const; + void updatePosition(double xCor, double yCor, LHCb::CaloHypo* hypo) const; + private: using IncCounter = Gaudi::Accumulators::Counter<>; using SCounter = Gaudi::Accumulators::StatCounter<float>; -- GitLab