From 1ae30cf90f5da577ae4367998c1393bf6983afde Mon Sep 17 00:00:00 2001 From: Tong Li <118776996+dasphy@users.noreply.github.com> Date: Wed, 24 Jul 2024 17:32:50 +0200 Subject: [PATCH] Photon ID shape parameter: correct very small negative values caused by computational precision in width of theta/module (#90) * 1st * 1st * add cluster shape variables * add cluster shape variables * add flags for shape var in EMB only * fix issues & add more variables * update * fix of width calculation * fix Ratio_E vs phi * fix fill empty cells & m_do_pi0_photon_shapeVar flag * up * up * apply crash fix from Giovanni * use edm4hep::labels::ShapeParameterNames * Update RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp Co-authored-by: Giovanni Marchiori <39376142+giovannimarchiori@users.noreply.github.com> * small fix * width theta logE weights * width theta logE weights 2 * adjust var definition * correct very small negative values caused by computational precision * fix * fix * fix small negative values Co-authored-by: Juraj Smiesko <34742917+kjvbrt@users.noreply.github.com> * fix small negative values --------- Co-authored-by: Tong Li Co-authored-by: Tong Li Co-authored-by: Giovanni Marchiori <39376142+giovannimarchiori@users.noreply.github.com> Co-authored-by: Juraj Smiesko <34742917+kjvbrt@users.noreply.github.com> --- .../src/components/AugmentClustersFCCee.cpp | 105 +++++++++++++----- 1 file changed, 80 insertions(+), 25 deletions(-) diff --git a/RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp b/RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp index c30b0f0f..8cf0900c 100644 --- a/RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp +++ b/RecFCCeeCalorimeter/src/components/AugmentClustersFCCee.cpp @@ -633,16 +633,33 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev // do pi0/photon shape var only for EMB if (m_do_photon_shapeVar && systemID == systemID_EMB) { - double w_theta; if (m_do_widthTheta_logE_weights) { - w_theta = (sumWeightLayer[layer+startPositionToFill] != 0.) ? sqrt(theta2_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill], 2)) : 0. ; + double w_theta2 = theta2_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumWeightLayer[layer+startPositionToFill], 2); + // Very small but negative values caused by computational precision are allowed, + // otherwise crash. + if (w_theta2 < -1e-9) { + error() << "w_theta2 in theta width calculation is negative: " << w_theta2 << endmsg; + return StatusCode::FAILURE; + } + width_theta[layer+startPositionToFill] = (sumWeightLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_theta2)) : 0. ; } else { - w_theta = (sumEnLayer[layer+startPositionToFill] > 0.) ? sqrt(theta2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2)) : 0. ; + double w_theta2 = theta2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(theta_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2); + // Very small but negative values caused by computational precision are allowed, + // otherwise crash. + if (w_theta2 < -1e-9) { + error() << "w_theta2 in theta width calculation is negative: " << w_theta2 << endmsg; + return StatusCode::FAILURE; + } + width_theta[layer+startPositionToFill] = (sumEnLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_theta2)) : 0. ; } - width_theta[layer+startPositionToFill] = w_theta; - - double w_module = (sumEnLayer[layer+startPositionToFill] > 0.) ? sqrt(module2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(module_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2)) : 0. ; - width_module[layer+startPositionToFill] = w_module; + double w_module2 = module2_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill] - std::pow(module_E_layer[layer+startPositionToFill] / sumEnLayer[layer+startPositionToFill], 2); + // Very small but negative values caused by computational precision are allowed, + // otherwise crash. + if (w_module2 < -1e-9) { + error() << "w_module2 in module width calculation is negative: " << w_module2 << endmsg; + return StatusCode::FAILURE; + } + width_module[layer+startPositionToFill] = (sumEnLayer[layer+startPositionToFill] > 0.) ? std::sqrt(std::fabs(w_module2)) : 0. ; double Ratio_E = (E_cell_Max[layer+startPositionToFill] - E_cell_secMax[layer+startPositionToFill]) / (E_cell_Max[layer+startPositionToFill] + E_cell_secMax[layer+startPositionToFill]); @@ -766,15 +783,34 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev double theta2_E_9Bin = theta2_E_7Bin + theta_m4 * theta_m4 * weightLog_m4 + theta_p4 * theta_p4 * weightLog_p4; double theta_E_9Bin = theta_E_7Bin + theta_m4 * weightLog_m4 + theta_p4 * weightLog_p4; - double _w_theta_3Bin = sqrt(theta2_E_3Bin / sum_weightLog_3Bin - std::pow(theta_E_3Bin / sum_weightLog_3Bin, 2)); - double _w_theta_5Bin = sqrt(theta2_E_5Bin / sum_weightLog_5Bin - std::pow(theta_E_5Bin / sum_weightLog_5Bin, 2)); - double _w_theta_7Bin = sqrt(theta2_E_7Bin / sum_weightLog_7Bin - std::pow(theta_E_7Bin / sum_weightLog_7Bin, 2)); - double _w_theta_9Bin = sqrt(theta2_E_9Bin / sum_weightLog_9Bin - std::pow(theta_E_9Bin / sum_weightLog_9Bin, 2)); - - width_theta_3Bin[layer+startPositionToFill] = (sum_weightLog_3Bin > 0.) ? _w_theta_3Bin : 0. ; - width_theta_5Bin[layer+startPositionToFill] = (sum_weightLog_5Bin > 0.) ? _w_theta_5Bin : 0. ; - width_theta_7Bin[layer+startPositionToFill] = (sum_weightLog_7Bin > 0.) ? _w_theta_7Bin : 0. ; - width_theta_9Bin[layer+startPositionToFill] = (sum_weightLog_9Bin > 0.) ? _w_theta_9Bin : 0. ; + double _w_theta_3Bin2 = theta2_E_3Bin / sum_weightLog_3Bin - std::pow(theta_E_3Bin / sum_weightLog_3Bin, 2); + double _w_theta_5Bin2 = theta2_E_5Bin / sum_weightLog_5Bin - std::pow(theta_E_5Bin / sum_weightLog_5Bin, 2); + double _w_theta_7Bin2 = theta2_E_7Bin / sum_weightLog_7Bin - std::pow(theta_E_7Bin / sum_weightLog_7Bin, 2); + double _w_theta_9Bin2 = theta2_E_9Bin / sum_weightLog_9Bin - std::pow(theta_E_9Bin / sum_weightLog_9Bin, 2); + + // Very small but negative values caused by computational precision are allowed, + // otherwise crash. + if (_w_theta_3Bin2 < -1e-9) { + error() << "_w_theta_3Bin2 in theta width calculation is negative: " << _w_theta_3Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_5Bin2 < -1e-9) { + error() << "_w_theta_5Bin2 in theta width calculation is negative: " << _w_theta_5Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_7Bin2 < -1e-9) { + error() << "_w_theta_7Bin2 in theta width calculation is negative: " << _w_theta_7Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_9Bin2 < -1e-9) { + error() << "_w_theta_9Bin2 in theta width calculation is negative: " << _w_theta_9Bin2 << endmsg; + return StatusCode::FAILURE; + } + + width_theta_3Bin[layer+startPositionToFill] = (sum_weightLog_3Bin > 0.) ? std::sqrt(std::fabs(_w_theta_3Bin2)) : 0. ; + width_theta_5Bin[layer+startPositionToFill] = (sum_weightLog_5Bin > 0.) ? std::sqrt(std::fabs(_w_theta_5Bin2)) : 0. ; + width_theta_7Bin[layer+startPositionToFill] = (sum_weightLog_7Bin > 0.) ? std::sqrt(std::fabs(_w_theta_7Bin2)) : 0. ; + width_theta_9Bin[layer+startPositionToFill] = (sum_weightLog_9Bin > 0.) ? std::sqrt(std::fabs(_w_theta_9Bin2)) : 0. ; } else { double theta2_E_3Bin = theta_m1 * theta_m1 * E_m1 + E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max[layer+startPositionToFill] + theta_p1 * theta_p1 * E_p1; double theta_E_3Bin = theta_m1 * E_m1 + E_cell_Max_theta[layer+startPositionToFill] * E_cell_Max[layer+startPositionToFill] + theta_p1 * E_p1; @@ -785,15 +821,34 @@ StatusCode AugmentClustersFCCee::execute([[maybe_unused]] const EventContext &ev double theta2_E_9Bin = theta2_E_7Bin + theta_m4 * theta_m4 * E_m4 + theta_p4 * theta_p4 * E_p4; double theta_E_9Bin = theta_E_7Bin + theta_m4 * E_m4 + theta_p4 * E_p4; - double _w_theta_3Bin = sqrt(theta2_E_3Bin / sum_E_3Bin - std::pow(theta_E_3Bin / sum_E_3Bin, 2)); - double _w_theta_5Bin = sqrt(theta2_E_5Bin / sum_E_5Bin - std::pow(theta_E_5Bin / sum_E_5Bin, 2)); - double _w_theta_7Bin = sqrt(theta2_E_7Bin / sum_E_7Bin - std::pow(theta_E_7Bin / sum_E_7Bin, 2)); - double _w_theta_9Bin = sqrt(theta2_E_9Bin / sum_E_9Bin - std::pow(theta_E_9Bin / sum_E_9Bin, 2)); - - width_theta_3Bin[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? _w_theta_3Bin : 0. ; - width_theta_5Bin[layer+startPositionToFill] = (sum_E_5Bin > 0.) ? _w_theta_5Bin : 0. ; - width_theta_7Bin[layer+startPositionToFill] = (sum_E_7Bin > 0.) ? _w_theta_7Bin : 0. ; - width_theta_9Bin[layer+startPositionToFill] = (sum_E_9Bin > 0.) ? _w_theta_9Bin : 0. ; + double _w_theta_3Bin2 = theta2_E_3Bin / sum_E_3Bin - std::pow(theta_E_3Bin / sum_E_3Bin, 2); + double _w_theta_5Bin2 = theta2_E_5Bin / sum_E_5Bin - std::pow(theta_E_5Bin / sum_E_5Bin, 2); + double _w_theta_7Bin2 = theta2_E_7Bin / sum_E_7Bin - std::pow(theta_E_7Bin / sum_E_7Bin, 2); + double _w_theta_9Bin2 = theta2_E_9Bin / sum_E_9Bin - std::pow(theta_E_9Bin / sum_E_9Bin, 2); + + // Very small but negative values caused by computational precision are allowed, + // otherwise crash. + if (_w_theta_3Bin2 < -1e-9) { + error() << "_w_theta_3Bin2 in theta width calculation is negative: " << _w_theta_3Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_5Bin2 < -1e-9) { + error() << "_w_theta_5Bin2 in theta width calculation is negative: " << _w_theta_5Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_7Bin2 < -1e-9) { + error() << "_w_theta_7Bin2 in theta width calculation is negative: " << _w_theta_7Bin2 << endmsg; + return StatusCode::FAILURE; + } + if (_w_theta_9Bin2 < -1e-9) { + error() << "_w_theta_9Bin2 in theta width calculation is negative: " << _w_theta_9Bin2 << endmsg; + return StatusCode::FAILURE; + } + + width_theta_3Bin[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? std::sqrt(std::fabs(_w_theta_3Bin2)) : 0. ; + width_theta_5Bin[layer+startPositionToFill] = (sum_E_5Bin > 0.) ? std::sqrt(std::fabs(_w_theta_5Bin2)) : 0. ; + width_theta_7Bin[layer+startPositionToFill] = (sum_E_7Bin > 0.) ? std::sqrt(std::fabs(_w_theta_7Bin2)) : 0. ; + width_theta_9Bin[layer+startPositionToFill] = (sum_E_9Bin > 0.) ? std::sqrt(std::fabs(_w_theta_9Bin2)) : 0. ; } E_fr_side_pm2[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? (sum_E_5Bin / sum_E_3Bin - 1.) : 0. ; E_fr_side_pm3[layer+startPositionToFill] = (sum_E_3Bin > 0.) ? (sum_E_7Bin / sum_E_3Bin - 1.) : 0. ;