diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.cpp
index 184913d185d6a6aaabc579fbe514b417932b3158..79d6550dccdf1347d0e3901b60b52160e752f9c4 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.cpp
@@ -28,11 +28,12 @@ CaloFutureCorrectionBase::CaloFutureCorrectionBase( const std::string& type, con
m_cmLoc = LHCb::CaloFutureAlgUtils::CaloFutureIdLocation( "ClusterMatch" );
// allocate memory for m_params. params[area][type]
- CaloFutureCorrection::ParamVector vect;
- CaloFutureCorrection::Parameters typePar_pair = std::make_pair( CaloFutureCorrection::Empty, vect );
- std::vector<CaloFutureCorrection::Parameters> vect_pair( CaloFutureCorrection::nT, typePar_pair );
+ CaloFutureCorrection::ParamVector vect;
+ CaloFutureCorrection::PrecalculatedParams precalculated;
+ CaloFutureCorrection::CachedParameters typePar_tuple = {CaloFutureCorrection::Empty, vect, precalculated};
+ std::vector<CaloFutureCorrection::CachedParameters> vect_tuple( CaloFutureCorrection::nT, typePar_tuple );
m_params.reserve( area_size ); // area_size=3 (outer, middle, inner)
- for ( auto area = 0; area < area_size; ++area ) m_params.push_back( vect_pair );
+ for ( auto area = 0; area < area_size; ++area ) m_params.push_back( vect_tuple );
}
StatusCode CaloFutureCorrectionBase::initialize() {
@@ -84,10 +85,10 @@ StatusCode CaloFutureCorrectionBase::finalize() {
for ( auto area = 0; area < 3; ++area ) {
for ( auto itype = 0; itype < static_cast<int>( CaloFutureCorrection::nT ); ++itype ) {
- int func = m_params[area][itype].first;
+ int func = m_params[area][itype].function;
if ( func == CaloFutureCorrection::Empty ) continue;
const auto& type = CaloFutureCorrection::typeName[itype];
- const auto& vec = m_params[area][itype].second;
+ const auto& vec = m_params[area][itype].values;
if ( !vec.empty() ) {
@@ -178,55 +179,70 @@ void CaloFutureCorrectionBase::Params( const std::string& paramName, const CaloF
continue;
}
- m_params[area][type].first = static_cast<CaloFutureCorrection::Function>( func );
+ m_params[area][type].function = static_cast<CaloFutureCorrection::Function>( func );
- m_params[area][type].second.clear();
- m_params[area][type].second.reserve( dim );
+ m_params[area][type].values.clear();
+ m_params[area][type].values.reserve( dim );
for ( int i = 0; i < dim; ++i ) {
- m_params[area][type].second.push_back( params[pos] );
+ m_params[area][type].values.push_back( params[pos] );
pos += step;
}
+ // precalculated exponent, sinh and cosh values: see 'polynomial functions' and 'Sshape function' in
+ // CaloFutureCorrectionBase::getCorrection and CaloFutureCorrectionBase::getCorrectionDerivative
+ if ( !m_params[area][type].values.empty() ) {
+ const auto& b = m_params[area][type].values[0];
+ constexpr double delta = 0.5;
+ m_params[area][type].precalculated.myexp = myexp( b );
+ auto sinh_val = ( b != 0 ) ? mysinh( delta / b ) : INFINITY;
+ m_params[area][type].precalculated.mysinh = sinh_val;
+ m_params[area][type].precalculated.mycosh = sqrt( 1. + sinh_val * sinh_val );
+ }
}
}
-const CaloFutureCorrection::Parameters& CaloFutureCorrectionBase::getParams( const CaloFutureCorrection::Type type,
- const LHCb::CaloCellID id ) const {
+const CaloFutureCorrection::CachedParameters&
+CaloFutureCorrectionBase::getParams( const CaloFutureCorrection::Type type, const LHCb::CaloCellID id ) const {
auto area = id.area();
return m_params[area][type];
}
-double CaloFutureCorrectionBase::getCorrection( const CaloFutureCorrection::Type type, const LHCb::CaloCellID id,
- double var, double def ) const {
+std::optional<double> CaloFutureCorrectionBase::getCorrectionImpl( const CaloFutureCorrection::Type type,
+ const LHCb::CaloCellID id, const double var,
+ double* derivative ) const {
+ bool doDerivative = ( derivative != nullptr );
+ double cor = 0.;
const auto& pars = getParams( type, id );
if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) {
const auto& name = CaloFutureCorrection::typeName[type];
debug() << "Correction type " << name << " to be applied on cluster (seed = " << id << ") is a '"
- << CaloFutureCorrection::funcName[pars.first] << "' function with params = " << pars.second << endmsg;
+ << CaloFutureCorrection::funcName[pars.function] << "' function with params = " << pars.values << endmsg;
}
// compute correction
- if ( pars.first == CaloFutureCorrection::Unknown || pars.second.empty() ) return def;
+ if ( pars.function == CaloFutureCorrection::Unknown || pars.values.empty() ) return std::nullopt;
// list accessor - not correction :
- if ( pars.first == CaloFutureCorrection::ParamList || pars.first == CaloFutureCorrection::GlobalParamList ) {
+ if ( pars.function == CaloFutureCorrection::ParamList || pars.function == CaloFutureCorrection::GlobalParamList ) {
warning() << " Param accessor is a fake function - no correction to be applied - return default value" << endmsg;
- return def;
+ return std::nullopt;
}
- double cor = def;
// polynomial correction
- const auto& temp = pars.second;
+ const auto& temp = pars.values;
+ // and precalculated exponents of some of these values
+ auto& temp_precalc = pars.precalculated;
// polynomial functions
- if ( pars.first == CaloFutureCorrection::Polynomial || pars.first == CaloFutureCorrection::InversPolynomial ||
- pars.first == CaloFutureCorrection::ExpPolynomial || pars.first == CaloFutureCorrection::ReciprocalPolynomial ) {
+ if ( pars.function == CaloFutureCorrection::Polynomial || pars.function == CaloFutureCorrection::InversPolynomial ||
+ pars.function == CaloFutureCorrection::ExpPolynomial ||
+ pars.function == CaloFutureCorrection::ReciprocalPolynomial ) {
double v = 1.;
cor = 0.;
for ( auto i = temp.begin(); i != temp.end(); ++i ) {
cor += ( *i ) * v;
- if ( pars.first == CaloFutureCorrection::ReciprocalPolynomial )
+ if ( pars.function == CaloFutureCorrection::ReciprocalPolynomial )
v = ( var == 0 ) ? 0. : v / var;
else
v *= var;
@@ -235,36 +251,79 @@ double CaloFutureCorrectionBase::getCorrection( const CaloFutureCorrection::Type
if ( UNLIKELY( msgLevel( MSG::VERBOSE ) ) ) verbose() << "cor = " << cor << endmsg;
#endif
}
- if ( pars.first == CaloFutureCorrection::InversPolynomial ) cor = ( cor == 0 ) ? def : 1. / cor;
- if ( pars.first == CaloFutureCorrection::ExpPolynomial ) cor = ( cor == 0 ) ? def : myexp( cor );
+ if ( pars.function == CaloFutureCorrection::InversPolynomial )
+ if ( cor != 0 ) cor = 1. / cor;
+ if ( pars.function == CaloFutureCorrection::ExpPolynomial ) {
+ if ( cor != 0 ) {
+ if ( var == 0.0 )
+ // if var == 0.0 then cor == temp[0]
+ cor = temp_precalc.myexp;
+ else
+ cor = myexp( cor );
+ }
+ }
+
+ if ( doDerivative ) {
+ if ( pars.function == CaloFutureCorrection::Polynomial ||
+ pars.function == CaloFutureCorrection::InversPolynomial ||
+ pars.function == CaloFutureCorrection::ExpPolynomial ) {
+ *derivative = 0.;
+ double v = 1.;
+ int cnt = 0;
+ auto i = temp.begin();
+ for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
+ *derivative += ( *i ) * cnt * v;
+ v *= var;
+ }
+ if ( pars.function == CaloFutureCorrection::InversPolynomial ) *derivative = -( *derivative ) * cor * cor;
+ if ( pars.function == CaloFutureCorrection::ExpPolynomial ) *derivative = ( *derivative ) * cor;
+ } else if ( pars.function == CaloFutureCorrection::ReciprocalPolynomial ) {
+ *derivative = 0.;
+ if ( var != 0 ) {
+ auto v = 1. / ( var * var );
+ int cnt = 0;
+ auto i = temp.begin();
+ for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
+ *derivative -= ( *i ) * cnt * v;
+ v /= var;
+ }
+ }
+ }
+ }
}
// sigmoid function
- else if ( pars.first == CaloFutureCorrection::Sigmoid ) {
+ else if ( pars.function == CaloFutureCorrection::Sigmoid ) {
if ( temp.size() == 4 ) {
- const auto& a = temp[0];
- const auto& b = temp[1];
- const auto& c = temp[2];
- const auto& d = temp[3];
- cor = a + b * mytanh( c * ( var + d ) );
+ const auto& a = temp[0];
+ const auto& b = temp[1];
+ const auto& c = temp[2];
+ const auto& d = temp[3];
+ auto mytanh_val = mytanh( c * ( var + d ) );
+ cor = a + b * mytanh_val;
+ if ( doDerivative ) *derivative = b * c * ( 1. - mytanh_val * mytanh_val );
} else {
Warning( "The power sigmoid function must have 4 parameters" ).ignore();
}
}
// Sshape function
- else if ( pars.first == CaloFutureCorrection::Sshape || pars.first == CaloFutureCorrection::SshapeMod ) {
+ else if ( pars.function == CaloFutureCorrection::Sshape || pars.function == CaloFutureCorrection::SshapeMod ) {
if ( temp.size() == 1 ) {
const auto& b = temp[0];
constexpr double delta = 0.5;
if ( b > 0 ) {
double arg = var / delta;
- if ( pars.first == CaloFutureCorrection::SshapeMod ) {
- arg *= mysinh( delta / b );
- } else if ( pars.first == CaloFutureCorrection::Sshape ) {
- arg *= mycosh( delta / b );
+ double csh = 1.;
+ if ( pars.function == CaloFutureCorrection::SshapeMod ) {
+ arg *= temp_precalc.mysinh;
+ csh = temp_precalc.mysinh;
+ } else if ( pars.function == CaloFutureCorrection::Sshape ) {
+ arg *= temp_precalc.mycosh;
+ csh = temp_precalc.mycosh;
}
cor = b * mylog( arg + std::sqrt( arg * arg + 1.0 ) );
+ if ( doDerivative ) *derivative = b / delta * csh / std::sqrt( arg * arg + 1. );
}
} else {
Warning( "The Sshape function must have 1 parameter" ).ignore();
@@ -272,16 +331,30 @@ double CaloFutureCorrectionBase::getCorrection( const CaloFutureCorrection::Type
}
// Shower profile function
- else if ( pars.first == CaloFutureCorrection::ShowerProfile ) {
+ else if ( pars.function == CaloFutureCorrection::ShowerProfile ) {
if ( temp.size() == 10 ) {
if ( var > 0.5 ) {
- cor = temp[0] * myexp( -temp[1] * var );
- cor += temp[2] * myexp( -temp[3] * var );
- cor += temp[4] * myexp( -temp[5] * var );
+ auto tmp1 = temp[0] * myexp( -temp[1] * var );
+ auto tmp2 = temp[2] * myexp( -temp[3] * var );
+ auto tmp3 = temp[4] * myexp( -temp[5] * var );
+ cor = tmp1;
+ cor += tmp2;
+ cor += tmp3;
+ if ( doDerivative ) {
+ *derivative = -temp[1] * tmp1;
+ *derivative += -temp[3] * tmp2;
+ *derivative += -temp[5] * tmp3;
+ }
} else {
- cor = 2.;
- cor -= temp[6] * myexp( -temp[7] * var );
- cor -= temp[8] * myexp( -temp[9] * var );
+ cor = 2.;
+ auto tmp1 = temp[6] * myexp( -temp[7] * var );
+ auto tmp2 = temp[8] * myexp( -temp[9] * var );
+ cor -= tmp1;
+ cor -= tmp2;
+ if ( doDerivative ) {
+ *derivative += temp[7] * tmp1;
+ *derivative += temp[9] * tmp2;
+ }
}
} else {
Warning( "The ShowerProfile function must have 10 parameters" ).ignore();
@@ -289,10 +362,13 @@ double CaloFutureCorrectionBase::getCorrection( const CaloFutureCorrection::Type
}
// Sinusoidal function
- else if ( pars.first == CaloFutureCorrection::Sinusoidal ) {
+ else if ( pars.function == CaloFutureCorrection::Sinusoidal ) {
if ( temp.size() == 1 ) {
- const double& A = temp[0];
- cor = A * mysin( 2 * M_PI * var );
+ const double& A = temp[0];
+ constexpr double twopi = 2. * M_PI;
+ const auto sin_val = mysin( twopi * var );
+ cor = A * sin_val;
+ if ( doDerivative ) *derivative = A * twopi * sqrt( 1. - sin_val * sin_val );
} else {
Warning( "The Sinusoidal function must have 1 parameter" ).ignore();
}
@@ -310,10 +386,10 @@ void CaloFutureCorrectionBase::checkParams() {
}
for ( int itype = 0; itype < static_cast<int>( CaloFutureCorrection::nT ); ++itype ) {
- int func = m_params[area][itype].first;
+ int func = m_params[area][itype].function;
bool ok = true;
if ( func != CaloFutureCorrection::Empty ) {
- if ( m_params[area][itype].second.size() == 0 ) {
+ if ( m_params[area][itype].values.size() == 0 ) {
warning() << "Empty parameter vector of type " << CaloFutureCorrection::typeName[itype] << " of function "
<< CaloFutureCorrection::funcName[func] << endmsg;
ok = false;
@@ -326,163 +402,12 @@ void CaloFutureCorrectionBase::checkParams() {
}
if ( UNLIKELY( msgLevel( MSG::DEBUG ) && ok ) )
debug() << " o Will apply correction '" << CaloFutureCorrection::funcName[itype] << "' as a '"
- << CaloFutureCorrection::funcName[func] << "' function of " << m_params[area][itype].second.size()
+ << CaloFutureCorrection::funcName[func] << "' function of " << m_params[area][itype].values.size()
<< " parameters" << endmsg;
if ( !ok ) {
- m_params[area][itype].second.clear();
- m_params[area][itype].first = CaloFutureCorrection::Empty;
+ m_params[area][itype].values.clear();
+ m_params[area][itype].function = CaloFutureCorrection::Empty;
}
}
}
}
-
-double CaloFutureCorrectionBase::getCorrectionDerivative( const CaloFutureCorrection::Type type,
- const LHCb::CaloCellID id, double var, double def ) const {
- const auto& pars = getParams( type, id );
-
- if ( msgLevel( MSG::DEBUG ) ) {
- const auto name = CaloFutureCorrection::typeName[type];
- debug() << "Derivative for Correction type " << name << " to be calculated for cluster (seed = " << id << ") is a '"
- << CaloFutureCorrection::funcName[pars.first] << "' function with params = " << pars.second << endmsg;
- }
-
- // compute correction
- double cor = def;
- if ( pars.first == CaloFutureCorrection::Unknown || pars.second.empty() ) return cor;
-
- // polynomial correction
- const auto& temp = pars.second;
-
- // polynomial functions
- double ds = 0.;
- if ( pars.first == CaloFutureCorrection::Polynomial ) {
- ds = 0.;
- double v = 1.;
- int cnt = 0;
- auto i = temp.begin();
- for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
- ds += ( *i ) * cnt * v;
- v *= var;
- }
- }
-
- else if ( pars.first == CaloFutureCorrection::InversPolynomial ) {
- double v = 1.;
- cor = 0.;
- for ( auto i = temp.begin(); i != temp.end(); ++i ) {
- cor += ( *i ) * v;
- v *= var;
- }
- cor = ( cor == 0 ) ? def : 1. / cor;
-
- v = 1.;
- ds = 0.;
- int cnt = 0;
- auto i = temp.begin();
- for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
- ds += ( *i ) * cnt * v;
- v *= var;
- }
- ds *= -cor * cor;
- }
-
- else if ( pars.first == CaloFutureCorrection::ExpPolynomial ) {
- double v = 1.;
- cor = 0.;
- for ( auto i = temp.begin(); i != temp.end(); ++i ) {
- cor += ( *i ) * v;
- v *= var;
- }
- cor = ( cor == 0 ) ? def : myexp( cor );
-
- ds = 0.;
- v = 1.;
- int cnt = 0;
- auto i = temp.begin();
- for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
- ds += ( *i ) * cnt * v;
- v *= var;
- }
-
- ds *= cor;
- }
-
- else if ( pars.first == CaloFutureCorrection::ReciprocalPolynomial ) {
- ds = 0.;
- if ( var != 0 ) {
- auto v = 1. / ( var * var );
- int cnt = 0;
- auto i = temp.begin();
- for ( ++i, cnt++; i != temp.end(); ++i, cnt++ ) {
- ds -= ( *i ) * cnt * v;
- v /= var;
- }
- }
- }
-
- // sigmoid function
- else if ( pars.first == CaloFutureCorrection::Sigmoid ) {
- ds = 0.;
- if ( temp.size() == 4 ) {
- // double a = temp[0];
- const auto& b = temp[1];
- const auto& c = temp[2];
- const auto& d = temp[3];
- ds = b * c * ( 1. - std::pow( mytanh( c * ( var + d ) ), 2 ) );
- } else {
- Warning( "The power sigmoid function must have 4 parameters" ).ignore();
- }
- }
-
- // Sshape function
- else if ( pars.first == CaloFutureCorrection::Sshape || pars.first == CaloFutureCorrection::SshapeMod ) {
- ds = 0.;
- if ( temp.size() == 1 ) {
- const auto& b = temp[0];
- double delta = 0.5;
- if ( b > 0 ) {
- double csh = 1.;
- if ( pars.first == CaloFutureCorrection::SshapeMod ) {
- csh = mysinh( delta / b );
- } else if ( pars.first == CaloFutureCorrection::Sshape ) {
- csh = mycosh( delta / b );
- }
- const auto arg = var / delta * csh;
- ds = b / delta * csh / std::sqrt( arg * arg + 1. );
- }
- } else {
- Warning( "The Sshape function must have 1 parameter" ).ignore();
- }
- }
-
- // Shower profile function
- else if ( pars.first == CaloFutureCorrection::ShowerProfile ) {
- ds = 0.;
- if ( temp.size() == 10 ) {
- if ( var > 0.5 ) {
- ds = -temp[0] * temp[1] * myexp( -temp[1] * var );
- ds += -temp[2] * temp[3] * myexp( -temp[3] * var );
- ds += -temp[4] * temp[5] * myexp( -temp[5] * var );
- } else {
- ds = temp[6] * temp[7] * myexp( -temp[7] * var );
- ds += temp[8] * temp[9] * myexp( -temp[9] * var );
- }
- } else {
- Warning( "The ShowerProfile function must have 10 parameters" ).ignore();
- }
- }
-
- // Sinusoidal function
- else if ( pars.first == CaloFutureCorrection::Sinusoidal ) {
- ds = 0.;
- if ( temp.size() == 1 ) {
- const auto& A = temp[0];
- constexpr double twopi = 2. * M_PI;
- ds = A * twopi * mycos( twopi * var );
- } else {
- Warning( "The Sinusoidal function must have 1 parameter" ).ignore();
- }
- }
-
- return ds;
-}
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.h b/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.h
index d1a8e36f4ca0ca62bd36d488aff45edf48a3be32..7ac485a32a45c62bec5d75b49aff0957d6089cf5 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.h
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureCorrectionBase.h
@@ -72,8 +72,21 @@ namespace CaloFutureCorrection {
Empty, // for empty preallocated places in m_params
Unknown // type of correction from DB that is unspecified here. MUST be the last item.
};
+
using ParamVector = std::vector<double>;
- using Parameters = std::pair<CaloFutureCorrection::Function, ParamVector>;
+
+ struct PrecalculatedParams {
+ double myexp;
+ double mysinh;
+ double mycosh;
+ };
+
+ struct CachedParameters {
+ CaloFutureCorrection::Function function;
+ ParamVector values;
+ PrecalculatedParams precalculated;
+ };
+
enum Type {
// NB: numbering is not continuous due to removed PRS and SPD/Converted_photons -related parameters for Run 1-2
// CondDB compatibility
@@ -156,8 +169,6 @@ namespace CaloFutureCorrection {
inline const std::string funcName[nF] = {"InversPolynomial", "Polynomial", "ExpPolynomial", "ReciprocalPolynomial",
"Sigmoid", "Sshape", "ShowerProfile", "SshapeMod",
"Sinusoidal", "ParamList", "GlobalParamList", "Unknown"};
- using NewVector = std::vector<std::vector<std::pair<CaloFutureCorrection::Function, ParamVector>>>;
-
} // namespace CaloFutureCorrection
class CaloFutureCorrectionBase : public GaudiTool {
@@ -197,27 +208,37 @@ public:
}
// accessors
- const CaloFutureCorrection::Parameters& getParams( const CaloFutureCorrection::Type type,
- const LHCb::CaloCellID id = LHCb::CaloCellID() ) const;
- inline CaloFutureCorrection::ParamVector getParamVector( const CaloFutureCorrection::Type type,
- const LHCb::CaloCellID id = LHCb::CaloCellID() ) const {
- return getParams( type, id ).second;
+ const CaloFutureCorrection::CachedParameters& getParams( const CaloFutureCorrection::Type type,
+ const LHCb::CaloCellID id = LHCb::CaloCellID() ) const;
+ inline CaloFutureCorrection::ParamVector getParamVector( const CaloFutureCorrection::Type type,
+ const LHCb::CaloCellID id = LHCb::CaloCellID() ) const {
+ return getParams( type, id ).values;
}
double getParameter( CaloFutureCorrection::Type type, unsigned int i, const LHCb::CaloCellID id = LHCb::CaloCellID(),
double def = 0. ) const {
const auto& params = getParams( type, id );
- const auto& data = params.second;
+ const auto& data = params.values;
return ( i < data.size() ) ? data[i] : def;
}
- /// calculate analytic derivative for a given function type in cell id at a given point var with default value def
- double getCorrectionDerivative( const CaloFutureCorrection::Type type, const LHCb::CaloCellID id, double var = 0.,
- double def = 0. ) const;
- /// propagate cov.m. cov0 according to Jacobian jac: cov1 = (jac * cov * jac^T), see comments in
- /// CaloFutureECorrection.cpp and CaloFutureSCorrection.cpp
- // void recalculate_cov(const TMatrixD &jac, const TMatrixDSym &cov0, TMatrixDSym &cov1) const;
- double getCorrection( const CaloFutureCorrection::Type type, const LHCb::CaloCellID id, double var = 0.,
- double def = 1. ) const;
+ std::optional<double> getCorrection( const CaloFutureCorrection::Type type, LHCb::CaloCellID id,
+ double var = 0. ) const {
+ return getCorrectionImpl( type, id, var, nullptr );
+ }
+
+ struct CorrectionResult {
+ double value;
+ double derivative;
+ CorrectionResult( double _value, double _derivative ) : value( _value ), derivative( _derivative ) {}
+ };
+
+ std::optional<CorrectionResult> getCorrectionAndDerivative( const CaloFutureCorrection::Type type,
+ const LHCb::CaloCellID id, const double var = 0. ) const {
+ double d{0};
+ auto r = getCorrectionImpl( type, id, var, &d );
+ if ( !r ) return std::nullopt;
+ return std::optional<CorrectionResult>( {*r, d} );
+ }
protected:
Gaudi::Property<std::string> m_conditionName{this, "ConditionName", "none"};
@@ -333,10 +354,12 @@ private:
return *( a->second );
}
+ std::optional<double> getCorrectionImpl( const CaloFutureCorrection::Type type, const LHCb::CaloCellID id,
+ const double var, double* derivative ) const;
+
private:
- typedef std::vector<std::vector<std::pair<CaloFutureCorrection::Function, CaloFutureCorrection::ParamVector>>>
- NewVector;
- NewVector m_params;
+ typedef std::vector<std::vector<CaloFutureCorrection::CachedParameters>> NewVector;
+ NewVector m_params;
const int area_size = 3; // size of m_params (outer, middle, inner)
Gaudi::Property<std::map<std::string, std::vector<double>>> m_optParams{this, "Parameters"};
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureECorrection.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureECorrection.cpp
index df6e3c3748bac62aedaf7b34ec5218660d4b58d0..9bb7f1aaa16fa03c855009509eac1bc01abd37f3 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureECorrection.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureECorrection.cpp
@@ -300,29 +300,37 @@ CaloFutureECorrection::calcECorrection( const struct ECorrInputParams& _params )
// are useful for debugging in case of changes in the correction function code
//
//// aG = const(X,Y,E), no need to calculate derivatives
- double aG = getCorrection( CaloFutureCorrection::alphaG, cellID ); // global Ecal factor
+ double aG = getCorrection( CaloFutureCorrection::alphaG, cellID ).value_or( 1. ); // global Ecal factor
//// aE = alphaE(eEcal)
- double aE = getCorrection( CaloFutureCorrection::alphaE, cellID, eEcal ); // longitudinal leakage
- double aB = getCorrection( CaloFutureCorrection::alphaB, cellID, bDist ); // lateral leakage
+ const auto aECorDer = getCorrectionAndDerivative( CaloFutureCorrection::alphaE, cellID, eEcal )
+ .value_or( CorrectionResult{1., 0.} ); // longitudinal leakage
+ const auto aE = aECorDer.value;
+ const auto aBCorDer = getCorrectionAndDerivative( CaloFutureCorrection::alphaB, cellID, bDist )
+ .value_or( CorrectionResult{1., 0.} ); // lateral leakage
+ const auto aB = aBCorDer.value;
//// aX = alphaX(Asx1)
- double aX = getCorrection( CaloFutureCorrection::alphaX, cellID, Asx ); // module frame dead material X-direction
+ const auto aXCorDer = getCorrectionAndDerivative( CaloFutureCorrection::alphaX, cellID, Asx )
+ .value_or( CorrectionResult{1., 0.} ); // module frame dead material X-direction
+ const auto aX = aXCorDer.value;
//// aY = alphaY(Asy1)
- double aY = getCorrection( CaloFutureCorrection::alphaY, cellID, Asy ); // module frame dead material Y-direction
+ const auto aYCorDer = getCorrectionAndDerivative( CaloFutureCorrection::alphaY, cellID, Asy )
+ .value_or( CorrectionResult{1., 0.} ); // module frame dead material Y-direction
+ const auto aY = aYCorDer.value;
if ( m_correctCovariance ) {
- DaE = getCorrectionDerivative( CaloFutureCorrection::alphaE, cellID, eEcal );
- DaB = getCorrectionDerivative( CaloFutureCorrection::alphaB, cellID, bDist );
- DaX = getCorrectionDerivative( CaloFutureCorrection::alphaX, cellID, Asx );
- DaY = getCorrectionDerivative( CaloFutureCorrection::alphaY, cellID, Asy );
+ DaE = aECorDer.derivative;
+ DaB = aBCorDer.derivative;
+ DaX = aXCorDer.derivative;
+ DaY = aYCorDer.derivative;
}
// angular correction
// assume dtheta to be independent of X,Y,E, although it may still be implicitly a bit dependent on X,Y
- double gT = getCorrection( CaloFutureCorrection::globalT, cellID, dtheta ); // incidence angle (delta)
- double dT = getCorrection( CaloFutureCorrection::offsetT, cellID, dtheta, 0. ); // incidence angle (delta)
+ double gT = getCorrection( CaloFutureCorrection::globalT, cellID, dtheta ).value_or( 1. ); // incidence angle (delta)
+ double dT = getCorrection( CaloFutureCorrection::offsetT, cellID, dtheta ).value_or( 0. ); // incidence angle (delta)
// Energy offset
double sinT = m_det->cellSine( cellID );
- double offset = getCorrection( CaloFutureCorrection::offset, cellID, sinT, 0. );
+ double offset = getCorrection( CaloFutureCorrection::offset, cellID, sinT ).value_or( 0. );
// Apply Ecal leakage corrections
double alpha = aG * aE * aB * aX * aY;
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.cpp
index d2d0043d20e59d4b1dae52b91e0c2af562fa8452..b9c25e0a7b51bde43b25187ed60891b11ecc9eed 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.cpp
@@ -103,13 +103,13 @@ StatusCode CaloFutureLCorrection::process( LHCb::span<LHCb::CaloHypo* const> hyp
// Corrected angle
- double gamma0 = getCorrection( CaloFutureCorrection::gamma0, cellID );
- double delta0 = getCorrection( CaloFutureCorrection::delta0, cellID );
+ double gamma0 = getCorrection( CaloFutureCorrection::gamma0, cellID ).value_or( 1. );
+ double delta0 = getCorrection( CaloFutureCorrection::delta0, cellID ).value_or( 1. );
// NB: gammaP(ePrs = 0) != 0, deltaP(ePrs = 0) != 0 and depend on cellID.area()
// get gammaP and deltaP parameters (depending on cellID.area() for each cluster
- double gammaP = getCorrection( CaloFutureCorrection::gammaP, cellID, 0. );
- double deltaP = getCorrection( CaloFutureCorrection::deltaP, cellID, 0. );
+ double gammaP = getCorrection( CaloFutureCorrection::gammaP, cellID, 0. ).value_or( 1. );
+ double deltaP = getCorrection( CaloFutureCorrection::deltaP, cellID, 0. ).value_or( 1. );
double g = gamma0 - gammaP;
double d = delta0 + deltaP;
@@ -118,7 +118,7 @@ StatusCode CaloFutureLCorrection::process( LHCb::span<LHCb::CaloHypo* const> hyp
cos_theta = temp / std::sqrt( tan_theta * tan_theta + temp * temp );
const auto dz_fps = cos_theta * tg_fps;
- updateCounters( area, dz_fps, g, d );
+ if ( m_counterFlag.value() ) updateCounters( area, dz_fps, g, d );
// Recompute Z position and fill CaloFuturePosition
const auto zCor = z0 + dz_fps;
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.h b/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.h
index d6ef1c76671cca696ffddbaa0f9dcae510df2621..e7803568ba8e570ade9434525a348fce686d0553 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.h
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureLCorrection.h
@@ -60,6 +60,8 @@ private:
mutable std::vector<SCounter> m_counterPerCaloFutureAreaGamma;
mutable std::vector<SCounter> m_counterPerCaloFutureAreaDelta;
+ Gaudi::Property<bool> m_counterFlag{this, "EnableCounters", true, "Enable collection of counter-based statistics."};
+
bool isHypoNotValid( const LHCb::CaloHypo* hypo ) const;
bool isEnergyNegative( const LHCb::CaloHypo* hypo ) const;
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
index 08607200e5a40148a76ea30d8b4a5ab122fa970e..79a00247bbe50d9ca7033e0cd1e67e566d74b77d 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.cpp
@@ -65,7 +65,7 @@ StatusCode CaloFutureSCorrection::process( LHCb::span<LHCb::CaloHypo* const> hyp
if ( m_hypos.end() == h ) return Error( "Invalid hypothesis!", StatusCode::SUCCESS );
if ( hypo->e() < 0. ) {
- ++m_counterSkipNegativeEnergyCorrection;
+ if ( m_counterFlag.value() ) ++m_counterSkipNegativeEnergyCorrection;
continue;
}
@@ -131,8 +131,10 @@ StatusCode CaloFutureSCorrection::process( LHCb::span<LHCb::CaloHypo* const> hyp
if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) ) { printDebugInfo( hypo, params, xBar, yBar, xCor, yCor ); }
updatePosition( xCor, yCor, hypo );
- m_counterDeltaX += xCor - xBar;
- m_counterDeltaY += yCor - yBar;
+ if ( m_counterFlag.value() ) {
+ m_counterDeltaX += xCor - xBar;
+ m_counterDeltaY += yCor - yBar;
+ }
if ( m_correctCovariance ) { updateCovariance( dXhy_dXcl, dYhy_dYcl, hypo ); }
}
@@ -210,57 +212,61 @@ CaloFutureSCorrection::calculateSCorrections( const struct SCorrInputParams& par
double Asx = -( xBar - seedPos.x() ) / CellSize;
double Asy = -( yBar - seedPos.y() ) / CellSize;
- // cache intermediate values
- const double Asx0 = Asx;
- const double Asy0 = Asy;
-
// Sshape correction :
- Asx = getCorrection( CaloFutureCorrection::shapeX, cellID, Asx, Asx ); // Asx1
- Asy = getCorrection( CaloFutureCorrection::shapeY, cellID, Asy, Asy ); // Asy1
-
- // // cache intermediate values: necessary for debugging by calculating numeric derivatives dn_shapeX, dn_shapeY below
- // const double Asx1 = Asx;
- // const double Asy1 = Asy;
+ auto AsxCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::shapeX, cellID, Asx ).value_or( CorrectionResult{Asx, 1.} );
+ Asx = AsxCorDer.value; // Asx1
+ const auto DshapeX = AsxCorDer.derivative;
+ auto AsyCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::shapeY, cellID, Asy ).value_or( CorrectionResult{Asy, 1.} );
+ Asy = AsyCorDer.value; // Asy1
+ const auto DshapeY = AsyCorDer.derivative;
// Angular correction (if any) [ NEW - inserted between Sshape and residual correction ]
const double xs = seedPos.x() - Asx * CellSize; // xscor
const double ys = seedPos.y() - Asy * CellSize; // yscor
const double thx = myatan2( xs, z );
const double thy = myatan2( ys, z );
- const double daX = getCorrection( CaloFutureCorrection::angularX, cellID, thx, 0. );
- const double daY = getCorrection( CaloFutureCorrection::angularY, cellID, thy, 0. );
+ const auto [daX, DangularX] =
+ getCorrectionAndDerivative( CaloFutureCorrection::angularX, cellID, thx ).value_or( CorrectionResult{0., 0.} );
+ const auto [daY, DangularY] =
+ getCorrectionAndDerivative( CaloFutureCorrection::angularY, cellID, thy ).value_or( CorrectionResult{0., 0.} );
Asx -= daX;
Asy -= daY;
- // cache intermediate values
- const double Asx2 = Asx;
- const double Asy2 = Asy;
-
// residual correction (if any):
- bool residualX_flag = false;
- double dcX = getCorrection( CaloFutureCorrection::residual, cellID, Asx, 0. );
+ auto dcXCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::residual, cellID, Asx ).value_or( CorrectionResult{0., 0.} );
+ auto dcX = dcXCorDer.value;
if ( dcX == 0. ) {
- dcX = getCorrection( CaloFutureCorrection::residualX, cellID, Asx, 0. ); // check X-specific correction
- residualX_flag = true;
+ // check X-specific correction
+ dcXCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::residualX, cellID, Asx ).value_or( CorrectionResult{0., 0.} );
+ dcX = dcXCorDer.value;
}
- bool residualY_flag = false;
- double dcY = getCorrection( CaloFutureCorrection::residual, cellID, Asy, 0. );
+ const auto DresidualX = dcXCorDer.derivative;
+ auto dcYCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::residual, cellID, Asy ).value_or( CorrectionResult{0., 0.} );
+ auto dcY = dcYCorDer.value;
if ( dcY == 0. ) {
- dcY = getCorrection( CaloFutureCorrection::residualY, cellID, Asy, 0. ); // check Y-specific correction
- residualY_flag = true;
+ // check Y-specific correction
+ dcYCorDer =
+ getCorrectionAndDerivative( CaloFutureCorrection::residualY, cellID, Asy ).value_or( CorrectionResult{0., 0.} );
+ dcY = dcYCorDer.value;
}
+ const auto DresidualY = dcYCorDer.derivative;
Asx -= dcX;
Asy -= dcY;
- // cache intermediate values
- const double Asx3 = Asx;
- const double Asy3 = Asy;
-
// left/right - up/down asymmetries correction (if any) :
- double ddcX = ( xBar < 0 ) ? getCorrection( CaloFutureCorrection::asymM, cellID, Asx, 0. )
- : getCorrection( CaloFutureCorrection::asymP, cellID, Asx, 0. );
- double ddcY = ( yBar < 0 ) ? getCorrection( CaloFutureCorrection::asymM, cellID, Asy, 0. )
- : getCorrection( CaloFutureCorrection::asymP, cellID, Asy, 0. );
+ const auto [ddcX, DasymX] =
+ getCorrectionAndDerivative( ( xBar < 0 ) ? CaloFutureCorrection::asymM : CaloFutureCorrection::asymP, cellID,
+ Asx )
+ .value_or( CorrectionResult{0., 0.} );
+ const auto [ddcY, DasymY] =
+ getCorrectionAndDerivative( ( yBar < 0 ) ? CaloFutureCorrection::asymM : CaloFutureCorrection::asymP, cellID,
+ Asy )
+ .value_or( CorrectionResult{0., 0.} );
Asx += ddcX; // Asx4
Asy += ddcY; // Asy4
@@ -297,29 +303,12 @@ CaloFutureSCorrection::calculateSCorrections( const struct SCorrInputParams& par
*/
if ( m_correctCovariance ) {
-
// calculation of the analytic derivatives:
// NB: printouts comparing analytic calculations with numeric derivatives which are commented-out below
// are useful for debugging in case of changes in the correction function code
if ( UNLIKELY( msgLevel( MSG::DEBUG ) ) )
debug() << "---------- analytic derivatives of individual S-correction functions ---------------" << endmsg;
- double DshapeX = getCorrectionDerivative( CaloFutureCorrection::shapeX, cellID, Asx0, 1. );
- double DshapeY = getCorrectionDerivative( CaloFutureCorrection::shapeY, cellID, Asy0, 1. );
-
- double DangularX = getCorrectionDerivative( CaloFutureCorrection::angularX, cellID, thx, 0. );
- double DangularY = getCorrectionDerivative( CaloFutureCorrection::angularY, cellID, thy, 0. );
-
- double DresidualX = getCorrectionDerivative(
- ( residualX_flag ? CaloFutureCorrection::residualX : CaloFutureCorrection::residual ), cellID, Asx2, 0. );
- double DresidualY = getCorrectionDerivative(
- ( residualY_flag ? CaloFutureCorrection::residualY : CaloFutureCorrection::residual ), cellID, Asy2, 0. );
-
- double DasymX = ( xBar < 0 ) ? getCorrectionDerivative( CaloFutureCorrection::asymM, cellID, Asx3, 0. )
- : getCorrectionDerivative( CaloFutureCorrection::asymP, cellID, Asx3, 0. );
- double DasymY = ( yBar < 0 ) ? getCorrectionDerivative( CaloFutureCorrection::asymM, cellID, Asy3, 0. )
- : getCorrectionDerivative( CaloFutureCorrection::asymP, cellID, Asy3, 0. );
-
double tx = xs / z;
double ty = ys / z;
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
index 66c6ff179666a9fb8bd72a673b51d4cb87718600..2cb08e03cfa0efe81d9dad7cd84e8755e6655132 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureSCorrection.h
@@ -72,6 +72,8 @@ private:
mutable IncCounter m_counterSkipNegativeEnergyCorrection{this, "Skip negative energy correction"};
mutable SCounter m_counterDeltaX{this, "Delta(X)"};
mutable SCounter m_counterDeltaY{this, "Delta(Y)"};
+
+ Gaudi::Property<bool> m_counterFlag{this, "EnableCounters", true, "Enable collection of counter-based statistics."};
};
// ============================================================================
#endif // CALOFUTURERECO_CALOFUTURESCORRECTION_H
diff --git a/CaloFuture/CaloFutureReco/src/CaloFutureShowerOverlapTool.cpp b/CaloFuture/CaloFutureReco/src/CaloFutureShowerOverlapTool.cpp
index 815e255da21cfab229f0023429c715cd54632a74..08bdee56f76bb2e0ca01db380657d11d54a7c98f 100644
--- a/CaloFuture/CaloFutureReco/src/CaloFutureShowerOverlapTool.cpp
+++ b/CaloFuture/CaloFutureReco/src/CaloFutureShowerOverlapTool.cpp
@@ -81,7 +81,7 @@ namespace LHCb::Calo::Tools {
double showerFraction( const CaloFutureCorrectionBase& shape, const double d3d, const unsigned int area ) {
CaloCellID cellID( CaloCellCode::CaloIndex::EcalCalo, area, 0, 0 ); // fake cell
- return std::clamp( shape.getCorrection( CaloFutureCorrection::profile, cellID, d3d, 0. ), 0., 1. );
+ return std::clamp( shape.getCorrection( CaloFutureCorrection::profile, cellID, d3d ).value_or( 0. ), 0., 1. );
}
double fraction( const CaloFutureCorrectionBase& shape, const DeCalorimeter& det, const CaloCluster& cluster,