From 28fe22409eac11eed9843ff353d859e750475610 Mon Sep 17 00:00:00 2001 From: Ali Bavarchee Date: Tue, 18 Nov 2025 19:27:44 +0100 Subject: [PATCH 1/2] Extend CaloClusterResolution for NeutralPID with cluster categorization and shower shapes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Major enhancements: - Add cluster categorization (Photons, pi0, Hadrons, PileUp) - Implement complete shower shape variables (r2, r2r4, κ, asym) - Extract 5×5 grid cell energies (ClE1-ClE25) - Add robust E19/E49 calculations ... @author Ali Bavarchee ali.bavarchee@cern.ch @date 2025-07-11 --- .../src/CaloClusterResolution.cpp | 780 +++++++++++++++++- 1 file changed, 745 insertions(+), 35 deletions(-) diff --git a/Calo/CaloCheckers/src/CaloClusterResolution.cpp b/Calo/CaloCheckers/src/CaloClusterResolution.cpp index 9967940ba..df29d052f 100644 --- a/Calo/CaloCheckers/src/CaloClusterResolution.cpp +++ b/Calo/CaloCheckers/src/CaloClusterResolution.cpp @@ -19,9 +19,19 @@ #include "Event/CaloDataFunctor.h" #include "Event/CaloHypos_v2.h" #include "Event/MCParticle.h" +#include "Event/CaloDigits_v2.h" #include "GaudiAlg/GaudiTupleAlg.h" #include "LHCbAlgs/Consumer.h" #include "Relations/RelationWeighted1D.h" +#include "CaloDet/CellParam.h" +#include +#include +#include +#include +#include +#include +#include +#include using ICaloTable = LHCb::CaloFuture2MC::ClusterTable; namespace { @@ -52,7 +62,7 @@ namespace { const LHCb::MCVertex* convVert = nullptr; for ( auto it_range : allParticles ) { auto part = it_range.to(); - if ( part->originVertex()->position().Z() > 12000.0 ) continue; + if ( part->originVertex()->position().Z() > 12450.0 ) continue; if ( part->mother() && part->mother() == matchParticle ) { if ( part->particleID().pid() == 11 ) { flagEle = !flagEle; @@ -138,7 +148,32 @@ namespace LHCb::Calo { * * @author Sasha Zenaiev oleksandr.zenaiev@cern.ch * @date 2020-01-29 + * ------------------------------------------ + * MODIFIED VERSION FOR NEUTRAL PID STUDIES + * Extended with comprehensive cluster categorization (Photons, pi0, Hadrons, PileUp) + * and advanced shower shape analysis for gamma/pi0 separation studies. + * + * Key modifications include: + * - Complete 5x5 grid cell energy extraction (ClE1-ClE25) for detailed + * shower profile analysis and machine learning applications + * - Robust cluster selection logic with configurable modes (PU, Signal, Photon, Hadron) + * - Complete shower shape variable implementation (r2, r2r4, k, asym) based on LHCb-PUB-2015-016 + * - Non-overlapping quadrant definitions for E4/E49 calculations + * - Enhanced energy ratio calculations (E19, E_seed/E_cl, etc.) + * - Detailed debugging output for cluster selection and shape diagnostics + * - Support for both CaloClusters and CaloHypotheses templates + * + * The 5x5 cell energy grid enables: + * - Detailed shower shape characterization beyond traditional variables + * - Machine learning training for particle identification + * - Studies of energy deposition patterns for different particle types + * - Validation of calorimeter simulation and reconstruction + * + * @author Ali Bavarchee ali.bavarchee@cern.ch + * @date 2025-27-10 */ + + template class ClusterResolution : public LHCb::Algorithm::ConsumernTuple( m_tuplePrefix ); @@ -172,33 +207,648 @@ namespace LHCb::Calo { const auto& [matchParticle, matchFraction] = bestMatch( cluster, table.relations( currentSeed ) ); bufferMatched += ( matchParticle != nullptr ); - if ( !matchParticle ) continue; - - auto currentZPos = RetrieveZPosition( cluster ); + + // Use CellParam for precise geometry information + const auto& seedParam = calo.cellParam( currentSeed ); + auto currentZPos = seedParam.z(); // More accurate z-position auto currentPars = RetrieveParameters( cluster ); - // check photon conversion - auto zConv = checkConversion( matchParticle, currentSeed, table ); - auto vertex = matchParticle->originVertex(); - auto deltaZ = currentZPos - vertex->position().Z(); - auto mom = matchParticle->momentum(); - auto xExtrap = vertex->position().X() + deltaZ * ( mom.Px() / mom.Pz() ); - auto yExtrap = vertex->position().Y() + deltaZ * ( mom.Py() / mom.Pz() ); - auto seedPos = calo.cellCenter( currentSeed ); - double CellSize = calo.cellSize( currentSeed ); - - auto energy = matchParticle->momentum().E(); - auto energyCluster = currentPars[2]; - float deltaEoverE = ( energyCluster - energy ) / energy; + auto energyCluster = currentPars[2]; // Get cluster energy for normalization + + // MC truth aggregation + double relWeightedE = 0.0; + double relSumWeights = 0.0; + double maxNormalizedWeight = 0.0; + double sumNormalizedWeights = 0.0; + int bestPdg = 0; + + // Variables for weighted MCParticle counting + double weightedMCParticleCount = 0.0; + std::vector> mcParticleWeights; // PDG, normalized weight + + { + auto relationsRange = table.relations( currentSeed ); + + // First pass: calculate total weight sum for proper normalization + double totalRawWeight = 0.0; + for ( const auto & rel : relationsRange ) { + double w = static_cast( rel.weight() ); + totalRawWeight += w; + } + + // Second pass: compute normalized weights with proper normalization + for ( const auto & rel : relationsRange ) { + const auto * mcPart = rel.to(); + double w = static_cast( rel.weight() ); + + // CORRECTED: Normalize by totalRawWeight, not cluster energy + double normalizedWeight = (totalRawWeight > 0) ? w / totalRawWeight : 0.0; + + relSumWeights += w; + sumNormalizedWeights += normalizedWeight; + + if (normalizedWeight > maxNormalizedWeight) { + maxNormalizedWeight = normalizedWeight; + } + + if ( mcPart ) { + relWeightedE += w * static_cast( mcPart->momentum().e() ); + weightedMCParticleCount += normalizedWeight; + mcParticleWeights.emplace_back(mcPart->particleID().pid(), normalizedWeight); + } + } + if ( matchParticle ) bestPdg = matchParticle->particleID().pid(); + } + + // UPDATED CLUSTER SELECTION LOGIC + bool isPUCluster = false; + bool isSignalCluster = false; + bool isPhotonCluster = false; + bool isHadronCluster = false; + int truePdgId = bestPdg; + + if (matchParticle) { + truePdgId = matchParticle->particleID().pid(); + + // PU CLUSTER: Normalized weight < 0.8 (cluster energy fraction from matched particle < 80%) + if (maxNormalizedWeight < m_puMaxFraction) { + isPUCluster = true; + } + + // SIGNAL CLUSTER: Match fraction > 0.1 (good quality match) + if (matchFraction > m_signalMinMatchFraction) { + isSignalCluster = true; + + // PHOTON CLUSTER: Photon or pi0 with good match + if (truePdgId == 22 || truePdgId == 111) { + isPhotonCluster = true; + } + + // HADRON CLUSTER: Hadron (|PDG| > 100) with good match + if (std::abs(truePdgId) > 100) { + isHadronCluster = true; + } + } + } + + // CONDITIONAL SELECTION + // Skip if we're in PU mode and this is not a PU cluster + if (m_PUMode && !isPUCluster) { + continue; + } + + // Skip if we're in signal mode and this is not a signal cluster + if (m_signalMode && !isSignalCluster) { + continue; + } + + // Skip if we're in photon mode and this is not a photon cluster + if (m_photonMode && !isPhotonCluster) { + continue; + } + + // Skip if we're in hadron mode and this is not a hadron cluster + if (m_hadronMode && !isHadronCluster) { + continue; + } + + // ========== DEBUG OUTPUT ========== + if (m_debugWeights) { + std::cout << "=================================================================" << std::endl; + std::cout << "DEBUG CLUSTER SELECTION:" << std::endl; + std::cout << " Cluster Seed: " << currentSeed.all() + << " (Energy: " << energyCluster << " MeV)" << std::endl; + std::cout << " Best Match PDG: " << truePdgId + << ", Match Fraction: " << matchFraction + << ", Max Normalized Weight: " << maxNormalizedWeight << std::endl; + std::cout << " Selection Flags - PU: " << isPUCluster + << ", Signal: " << isSignalCluster + << ", Photon: " << isPhotonCluster + << ", Hadron: " << isHadronCluster << std::endl; + + if (m_PUMode && isPUCluster) { + std::cout << " -> SELECTED AS PU CLUSTER (maxNormalizedWeight < " << m_puMaxFraction << ")" << std::endl; + } else if (m_signalMode && isSignalCluster) { + std::cout << " -> SELECTED AS SIGNAL CLUSTER (matchFraction > " << m_signalMinMatchFraction << ")" << std::endl; + } else if (m_photonMode && isPhotonCluster) { + std::cout << " -> SELECTED AS PHOTON CLUSTER (PDG=" << truePdgId << ")" << std::endl; + } else if (m_hadronMode && isHadronCluster) { + std::cout << " -> SELECTED AS HADRON CLUSTER (|PDG| > 100)" << std::endl; + } + std::cout << "=================================================================" << std::endl; + } + + // Check photon conversion (only for photon clusters) + std::optional zConv; + if (isPhotonCluster && matchParticle) { + zConv = checkConversion( matchParticle, currentSeed, table ); + } + + auto vertex = matchParticle ? matchParticle->originVertex() : nullptr; + auto deltaZ = currentZPos - (vertex ? vertex->position().Z() : 0.0); + auto mom = matchParticle ? matchParticle->momentum() : Gaudi::LorentzVector(); + auto xExtrap = vertex ? vertex->position().X() + deltaZ * ( mom.Px() / (mom.Pz() != 0 ? mom.Pz() : 1.0) ) : 0.0; + auto yExtrap = vertex ? vertex->position().Y() + deltaZ * ( mom.Py() / (mom.Pz() != 0 ? mom.Pz() : 1.0) ) : 0.0; + auto seedPos = seedParam.center(); // More accurate position + double CellSize = seedParam.size(); // Direct from CellParam + + auto energy = matchParticle ? matchParticle->momentum().E() : 0.0; + float deltaEoverE = matchParticle ? ( energyCluster - energy ) / energy : -999.0; bufferDeltaEOverE += deltaEoverE; - buffetDeltaX += xExtrap - currentPars[0]; + bufferDeltaX += xExtrap - currentPars[0]; bufferDeltaY += yExtrap - currentPars[1]; - double uncorE( -999. ), uncorX( -999. ), uncorY( -999. ); - if ( !RetrieveUncorrectedEXY( cluster, calo, uncorE, uncorX, uncorY ).isSuccess() ) - uncorE = uncorX = uncorY = -999.; + // FIXED: Ηandle uncorrected energy retrieval + double uncorE = -999.0, uncorX = -999.0, uncorY = -999.0; + StatusCode uncorStatus = RetrieveUncorrectedEXY( cluster, calo, uncorE, uncorX, uncorY ); + if ( !uncorStatus.isSuccess() ) { + uncorE = -999.0; + uncorX = -999.0; + uncorY = -999.0; + } + + // ========== CORRECTED PSEUDORAPIDITY CALCULATION ========== + // eta = -ln(tan(theta/2)) where theta is the angle from beam axis + + double eta = -999.0; + double etaMC = -999.0; + double cosTheta = -999.0; + double cosThetaMC = -999.0; + + // For MC truth: use particle momentum direction + if (matchParticle) { + const auto& momentum = matchParticle->momentum(); + double pz = momentum.Pz(); + double pt = std::hypot(momentum.Px(), momentum.Py()); + double p = momentum.P(); + + if (p > 0) { + cosThetaMC = pz / p; // cos(theta) = pz/p + + // Using the identity: eta = asinh(pz/pt) + // This is numerically more stable than the log(tan) formula + // and equivalent to eta = -ln(tan(theta/2)) + if (pt > 0) { + etaMC = std::asinh(pz / pt); + } else if (pz > 0) { + etaMC = 10.0; // Large positive eta for forward particles + } else if (pz < 0) { + etaMC = -10.0; // Large negative eta for backward particles + } + } + } + + // For reconstructed clusters: use position direction as approximation + { + double r = std::hypot(seedPos.x(), seedPos.y()); + double mag = std::hypot(r, currentZPos); + + if (mag > 0) { + cosTheta = currentZPos / mag; // cos(theta) = z/r + + // Using position-based approximation: eta = -ln(tan(theta/2)) + // where theta = acos(cosTheta) + double theta = std::acos(cosTheta); + double tanHalfTheta = std::tan(theta / 2.0); + + if (tanHalfTheta > 0) { + eta = -std::log(tanHalfTheta); + // Add sign based on z position + if (currentZPos < 0) eta = -eta; + } else { + eta = (currentZPos > 0) ? 10.0 : -10.0; + } + } else { + // Handle the case where cluster is on the beam axis + cosTheta = (currentZPos > 0) ? 1.0 : -1.0; + eta = (currentZPos > 0) ? 10.0 : -10.0; + } + } + + // Compute transverse energy using corrected eta + double Et = (energyCluster > 0 && std::abs(eta) < 10.0) ? + energyCluster / std::cosh(eta) : 0.; + + // Get mother particle ID + int motherID = 0; + if (matchParticle && matchParticle->mother()) { + motherID = matchParticle->mother()->particleID().pid(); + } + + // Cluster shape variables + double E19 = 0.; // seed energy / 3x3 energy + double E49 = 0.; // E4/E9 ratio (max 2x2 / 3x3) - CORRECTED DEFINITION + std::array cellEnergies = {0}; // For 5x5 grid + + // Get entries for both Clusters and Hypotheses + auto getEntries = [&] { + if constexpr (std::is_same_v) { + return cluster.clusters()[0].entries(); + } else { + return cluster.entries(); + } + }; + const auto entries = getEntries(); + + // Extract seed components + const int seedArea = currentSeed.area(); + const int seedRow = currentSeed.row(); + const int seedCol = currentSeed.col(); + + // Get the seed cell parameters + const auto& seedCell = calo.cellParam(currentSeed); + + // Create a map of all cells in the cluster using row/col as key + // FIXED: Use energy * fraction like in CaloFutureHypoEstimator.cpp + std::map, double> cellEnergyMap; + double maxCellEnergy = 0.; + double seedCellEnergy = 0.; + + for (const auto& entry : entries) { + auto id = entry.cellID(); + // CORRECTED: Use energy * fraction like in HypoEstimator + double energy = entry.energy() * entry.fraction(); + cellEnergyMap[{id.row(), id.col()}] = energy; + + if (energy > maxCellEnergy) { + maxCellEnergy = energy; + } + if (id == currentSeed) { + seedCellEnergy = energy; + } + } + + // Get seed cell energy from the map to ensure consistency + auto seedIt = cellEnergyMap.find({seedRow, seedCol}); + if (seedIt != cellEnergyMap.end()) { + seedCellEnergy = seedIt->second; + } + + // Create a map for the full 5x5 grid + std::map, double> cellMap; + + // Populate cellMap for the full 5x5 grid + for (int dr = -2; dr <= 2; ++dr) { + for (int dc = -2; dc <= 2; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + + // Look for this cell in the cluster + auto it = cellEnergyMap.find({targetRow, targetCol}); + if (it != cellEnergyMap.end()) { + cellMap[{targetRow, targetCol}] = it->second; + } else { + // Cell not in cluster, set to 0 + cellMap[{targetRow, targetCol}] = 0.0; + } + } + } + + // Fill cellEnergies from cellMap + for (int dr = -2; dr <= 2; ++dr) { + for (int dc = -2; dc <= 2; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + int index = (dr + 2) * 5 + (dc + 2); + + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + cellEnergies[index] = it->second; + } + } + } + + // +++++++++ CORRECTED NON-OVERLAPPING QUADRANT DEFINITIONS +++++++++ + // Following the same logic as CaloFutureHypoEstimator.cpp + std::array e4s = {0.0, 0.0, 0.0, 0.0}; // Four non-overlapping 2x2 quadrants + double E9 = 0.0; // Total energy in 3x3 grid + double E4 = 0.0; // Maximum energy in any 2x2 quadrant + + // Loop over 3x3 grid around seed for E9 calculation + for (int dr = -1; dr <= 1; ++dr) { + for (int dc = -1; dc <= 1; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + double cellEnergy = it->second; + E9 += cellEnergy; + } + } + } + + // Now calculate the four non-overlapping 2x2 quadrants + // Quadrant 0: top-left [seedRow-1 to seedRow, seedCol-1 to seedCol] + for (int dr = -1; dr <= 0; ++dr) { + for (int dc = -1; dc <= 0; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + e4s[0] += it->second; + } + } + } + + // Quadrant 1: top-right [seedRow-1 to seedRow, seedCol to seedCol+1] + for (int dr = -1; dr <= 0; ++dr) { + for (int dc = 0; dc <= 1; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + e4s[1] += it->second; + } + } + } + + // Quadrant 2: bottom-right [seedRow to seedRow+1, seedCol to seedCol+1] + for (int dr = 0; dr <= 1; ++dr) { + for (int dc = 0; dc <= 1; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + e4s[2] += it->second; + } + } + } + + // Quadrant 3: bottom-left [seedRow to seedRow+1, seedCol-1 to seedCol] + for (int dr = 0; dr <= 1; ++dr) { + for (int dc = -1; dc <= 0; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find({targetRow, targetCol}); + if (it != cellMap.end()) { + e4s[3] += it->second; + } + } + } + + // E4 is the MAXIMUM energy among the four non-overlapping 2x2 quadrants + E4 = *std::max_element(e4s.begin(), e4s.end()); + + // ++++++++++ CORRECTED E19 CALCULATION WITH ROBUST CHECKS ++++++++++ + // CRITICAL FIX: Ensure consistent energy definitions + // Use the same energy definition for seedCellEnergy as used in E9 calculation + auto seedCellIt = cellMap.find({seedRow, seedCol}); + if (seedCellIt != cellMap.end()) { + seedCellEnergy = seedCellIt->second; // Consistent energy definition + } else { + seedCellEnergy = 0.0; + if (m_debugWeights) { + std::cout << "WARNING: Seed cell not found in cellMap for E19 calculation!" << std::endl; + } + } + + // Calculate E19 with robust checks + if (E9 > 1e-16 && seedCellEnergy > 0) { // Use energy threshold to avoid numerical issues + E19 = seedCellEnergy / E9; + + // CRITICAL: Additional physical sanity checks + if (E19 > 1.2) { + // This should never happen - indicates serious calculation error + if (m_debugWeights) { + std::cout << "WARNING: Unphysical E19 = " << E19 + << " (seedEnergy=" << seedCellEnergy + << ", E9=" << E9 + << ", clusterEnergy=" << energyCluster << ")" << std::endl; + } + E19 = 1e-49; // Clamp to maximum physical value + } + + // Ensure E19 is physically meaningful (0-1) + E19 = std::max(0.0, std::min(1.0, E19)); + + } else { + E19 = 0.; + if (m_debugWeights && E9 <= 1e-6) { + std::cout << "WARNING: E9 too small for E19 calculation: " << E9 << std::endl; + } + } + + // Compute E49 ratio (E4/E9) where E4 = max 2x2 quadrant energy, E9 = 3x3 grid energy + if (E9 > 1e-16) { + E49 = E4 / E9; + // Ensure E49 is physically meaningful (0-1) + if (E49 > 1.2) { + if (m_debugWeights) { + std::cout << "WARNING: Unphysical E49 = " << E49 + << " (E4=" << E4 << ", E9=" << E9 << ")" << std::endl; + } + E49 = 1e-49; + } + E49 = std::max(0.0, std::min(1.0, E49)); + } else { + E49 = 1e-49; + } + + // E19 diagnostics for debugging + if (m_debugWeights && (E19 < 0.2 || E19 > 1.0)) { + std::cout << "E19 DIAGNOSTICS - Cluster: " << currentSeed.all() << std::endl; + std::cout << " E19 = " << E19 << std::endl; + std::cout << " seedCellEnergy = " << seedCellEnergy << std::endl; + std::cout << " E9 = " << E9 << std::endl; + std::cout << " energyCluster = " << energyCluster << std::endl; + std::cout << " nClusterCells = " << entries.size() << std::endl; + + // Print 3x3 grid energies + std::cout << " 3x3 Grid Energies:" << std::endl; + for (int dr = -1; dr <= 1; ++dr) { + for (int dc = -1; dc <= 1; ++dc) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find({targetRow, targetCol}); + double energy = (it != cellMap.end()) ? it->second : 0.0; + std::cout << " [" << dr << "," << dc << "]: " << energy; + if (dr == 0 && dc == 0) std::cout << " <-- SEED"; + std::cout << std::endl; + } + } + } + + // ++++++++++ CORRECTED SHOWER SHAPE CALCULATION WITH SAFEGUARDS ++++++++++ + // Following the definition from the document: r2 = = S_XX + S_YY + // where S_XX, S_YY, S_XY are the position spread matrix elements + + double S_XX = 0.0; + double S_YY = 0.0; + double S_XY = 0.0; + double r4 = 0.0; // For calculation + double sumEnergyForShowerShape = 0.0; + + // Use cluster center from parameters (should be energy-weighted) + double x_c = currentPars[0]; + double y_c = currentPars[1]; + + // First pass: calculate S_XX, S_YY, S_XY and r4 + for (const auto& entry : entries) { + auto id = entry.cellID(); + // Use consistent energy definition: energy * fraction + double cellEnergy = entry.energy() * entry.fraction(); + + // CRITICAL FIX: Ensure cell energy is non-negative + if (cellEnergy < 0.0) { + if (m_debugWeights) { + std::cout << "WARNING: Negative cell energy detected: " << cellEnergy + << " for cell " << id.all() << std::endl; + } + cellEnergy = 0.0; // Set negative energies to zero + } + + const auto& cell = calo.cellParam(id); + auto cellPos = cell.center(); + + const double dx = cellPos.x() - x_c; + const double dy = cellPos.y() - y_c; + const double r2_component = dx*dx + dy*dy; + const double r4_component = r2_component * r2_component; // r⁴ + + // Calculate matrix elements as defined in ShowerShape.pdf Eq. (1) and (2) + S_XX += cellEnergy * dx * dx; + S_YY += cellEnergy * dy * dy; + S_XY += cellEnergy * dx * dy; + r4 += cellEnergy * r4_component; + + sumEnergyForShowerShape += cellEnergy; + } + // CRITICAL FIX: Proper normalization with safeguards + if (sumEnergyForShowerShape > 1e-6) { // Use small threshold to avoid floating point issues + S_XX /= sumEnergyForShowerShape; + S_YY /= sumEnergyForShowerShape; + S_XY /= sumEnergyForShowerShape; + r4 /= sumEnergyForShowerShape; + } else { + // If no energy, set all to zero + S_XX = 0.0; + S_YY = 0.0; + S_XY = 0.0; + r4 = 0.0; + if (m_debugWeights) { + std::cout << "WARNING: Very small or zero sumEnergyForShowerShape: " + << sumEnergyForShowerShape << std::endl; + } + } + + // Compute the shower shape variables: + + // I. r2 (shower shape): r2 = = S_XX + S_YY + double r2 = S_XX + S_YY; + + // CRITICAL FIX: Ensure r2 is non-negative + if (r2 < 0.0) { + if (m_debugWeights) { + std::cout << "WARNING: Negative r2 detected: " << r2 + << ", S_XX: " << S_XX << ", S_YY: " << S_YY + << ", sumEnergy: " << sumEnergyForShowerShape << std::endl; + } + r2 = 0.0; // Force to zero if negative due to numerical issues + } + + // II. r2r4 - r2r4 = 1 - (2 / ) + double r2r4 = 0.0; + if (r4 > 1e-12 && r2 > 1e-12) { // Avoid division by very small numbers + double ratio = (r2 * r2) / r4; + // Ensure ratio is physically meaningful (0-1) + if (ratio >= 0.0 && ratio <= 1.0) { + r2r4 = 1.0 - ratio; + } else { + r2r4 = 0.0; // Default for unphysical cases + if (m_debugWeights && ratio > 1.0) { + std::cout << "WARNING: Unphysical r2r4 ratio: " << ratio + << " (r2=" << r2 << ", r4=" << r4 << ")" << std::endl; + } + } + // To ensure physical range + r2r4 = std::max(0.0, std::min(1.0, r2r4)); + } + + // III. k (kappa) - related to eigenvalue ratio + double kappa = 0.0; + double denominator = (S_XX + S_YY) * (S_XX + S_YY); + if (denominator > 1e-12) { // Avoid division by very small numbers + double determinant = S_XX * S_YY - S_XY * S_XY; + // Ensure determinant is non-negative for eigenvalue calculation + if (determinant >= 0.0) { + double kappa_squared = 1.0 - 4.0 * determinant / denominator; + // Ensure kappa_squared is non-negative + if (kappa_squared >= 0.0) { + kappa = std::sqrt(kappa_squared); + } + } + kappa = std::max(0.0, std::min(1.0, kappa)); + } + + // IV. asym (asymmetry): correlation between X and Y + double asym = 0.0; + if (S_XX > 1e-12 && S_YY > 1e-12) { // Avoid division by very small numbers + double correlation = S_XY / std::sqrt(S_XX * S_YY); + // Ensure correlation is within valid range [-1, 1] + if (std::isfinite(correlation)) { + asym = std::max(-1.0, std::min(1.0, correlation)); + } + } + + // Also compute energy ratios mentioned in the document + double E_seed = seedCellEnergy; + double E_cl = energyCluster; + double E_seed_over_E_cl = (E_cl > 0) ? E_seed / E_cl : 0.0; + + // For E_2nd, we need to find the second most energetic cell + double E_2nd = 0.0; + for (const auto& entry : entries) { + auto id = entry.cellID(); + if (id == currentSeed) continue; // skip seed + + double cellEnergy = entry.energy() * entry.fraction(); + if (cellEnergy > E_2nd) { + E_2nd = cellEnergy; + } + } + double E_seed_plus_E_2nd_over_E_cl = (E_cl > 0) ? (E_seed + E_2nd) / E_cl : 0.0; + + // Ensure energy ratios are physically meaningful + E_seed_over_E_cl = std::max(0.0, std::min(1.0, E_seed_over_E_cl)); + E_seed_plus_E_2nd_over_E_cl = std::max(0.0, std::min(1.0, E_seed_plus_E_2nd_over_E_cl)); + + // Debug output for shower shape diagnostics + if (m_debugWeights && (r2 < 0.0 || !std::isfinite(r2))) { + std::cout << "SHOWER SHAPE DIAGNOSTICS:" << std::endl; + std::cout << " Cluster Energy: " << energyCluster << std::endl; + std::cout << " Sum Energy for Shower Shape: " << sumEnergyForShowerShape << std::endl; + std::cout << " S_XX: " << S_XX << ", S_YY: " << S_YY << ", S_XY: " << S_XY << std::endl; + std::cout << " r2: " << r2 << std::endl; + std::cout << " Number of cells: " << entries.size() << std::endl; + + // Print first few cell energies for debugging + int count = 0; + for (const auto& entry : entries) { + if (count++ < 5) { // Print first 5 cells + double cellEnergy = entry.energy() * entry.fraction(); + std::cout << " Cell " << count << ": energy=" << cellEnergy + << ", ID=" << entry.cellID().all() << std::endl; + } + } + } + + // Count how many cells are in the 5x5 grid + int nCellsIn5x5 = std::count_if(cellEnergies.begin(), cellEnergies.end(), + [](double e) { return e > 0; }); + + // ========== TUPLE FILLING ========== + //auto sc = tuple->column( "matchFraction", matchFraction ); auto sc = tuple->column( "maxMatchFraction", matchFraction ); + sc &= tuple->column( "ClTrueE", static_cast(relWeightedE) ); + sc &= tuple->column( "ClTrueSumWeight", static_cast(relSumWeights) ); + sc &= tuple->column( "ClTruePdg", bestPdg ); + sc &= tuple->column( "maxNormalizedWeight", static_cast(maxNormalizedWeight) ); + sc &= tuple->column( "isPUCluster", isPUCluster ); + sc &= tuple->column( "isSignalCluster", isSignalCluster ); + sc &= tuple->column( "isPhotonCluster", isPhotonCluster ); + sc &= tuple->column( "isHadronCluster", isHadronCluster ); sc &= tuple->column( "energy", energy ); sc &= tuple->column( "energyCluster", energyCluster ); sc &= tuple->column( "xCluster", currentPars[0] ); @@ -209,20 +859,69 @@ namespace LHCb::Calo { sc &= tuple->column( "yClusterUncor", uncorY ); sc &= tuple->column( "xExtrap", xExtrap ); sc &= tuple->column( "yExtrap", yExtrap ); - sc &= tuple->column( "xOrigVertex", vertex->position().X() ); - sc &= tuple->column( "yOrigVertex", vertex->position().Y() ); - sc &= tuple->column( "zOrigVertex", vertex->position().Z() ); + sc &= tuple->column( "xOrigVertex", vertex ? vertex->position().X() : -999.0 ); + sc &= tuple->column( "yOrigVertex", vertex ? vertex->position().Y() : -999.0 ); + sc &= tuple->column( "zOrigVertex", vertex ? vertex->position().Z() : -999.0 ); sc &= tuple->column( "flagConv", zConv.has_value() ); sc &= tuple->column( "zConv", zConv.value_or( -999.0 ) ); - sc &= tuple->column( "area", currentSeed.area() ); + sc &= tuple->column( "area", seedArea ); sc &= tuple->column( "cellSize", CellSize ); - sc &= tuple->column( "cellCol", currentSeed.col() ); - sc &= tuple->column( "cellRow", currentSeed.row() ); + sc &= tuple->column( "cellCol", seedCol ); + sc &= tuple->column( "cellRow", seedRow ); sc &= tuple->column( "xSeed", seedPos.x() ); sc &= tuple->column( "ySeed", seedPos.y() ); - sc &= tuple->column( "px", matchParticle->momentum().Px() ); - sc &= tuple->column( "py", matchParticle->momentum().Py() ); - sc &= tuple->column( "pz", matchParticle->momentum().Pz() ); + sc &= tuple->column( "px", matchParticle ? matchParticle->momentum().Px() : -999.0 ); + sc &= tuple->column( "py", matchParticle ? matchParticle->momentum().Py() : -999.0 ); + sc &= tuple->column( "pz", matchParticle ? matchParticle->momentum().Pz() : -999.0 ); + + // ++++++++++ CORRECTED PSEUDORAPIDITY AND ANGLE VARIABLES ++++++++++ + sc &= tuple->column( "Et", Et ); + sc &= tuple->column( "eta", eta ); // Cluster position-based eta + sc &= tuple->column( "etaMC", etaMC ); // MC truth momentum-based eta + sc &= tuple->column( "cosTheta", cosTheta ); // cos(theta) from cluster position + sc &= tuple->column( "cosThetaMC", cosThetaMC ); // cos(theta) from MC momentum + sc &= tuple->column( "MCMotherID", motherID ); + sc &= tuple->column( "cellID", currentSeed.all() ); + sc &= tuple->column( "ClusterE", energyCluster ); + sc &= tuple->column( "E19", E19 ); + sc &= tuple->column( "E49", E49 ); // CORRECTED: E4/E9 ratio (max 2x2 / 3x3) + sc &= tuple->column( "E4", E4 ); // CORRECTED: Maximum energy in any 2x2 quadrant + sc &= tuple->column( "E9", E9 ); // Energy in 3x3 grid + sc &= tuple->column( "seedCellEnergy", seedCellEnergy ); // Seed cell energy + + // ********** UPDATED SHOWER SHAPE VARIABLES ********** + sc &= tuple->column( "ShowerShape_r2", r2 ); // Main lateral dispersion (r2) + sc &= tuple->column( "ShowerShape_r2r4", r2r4 ); + sc &= tuple->column( "ShowerShape_kappa", kappa ); + sc &= tuple->column( "ShowerShape_asym", asym ); + sc &= tuple->column( "ShowerShape_S_XX", S_XX ); + sc &= tuple->column( "ShowerShape_S_YY", S_YY ); + sc &= tuple->column( "ShowerShape_S_XY", S_XY ); + sc &= tuple->column( "E_seed_over_E_cl", E_seed_over_E_cl ); + sc &= tuple->column( "E_seed_plus_E_2nd_over_E_cl", E_seed_plus_E_2nd_over_E_cl ); + sc &= tuple->column( "E_2nd", E_2nd ); + + // Keep backward compatibility for existing analyses + sc &= tuple->column( "lateral", r2 ); // Alias for ShowerShape_r2 + + // Add individual quadrant energies for detailed analysis + sc &= tuple->column( "E4_quad0", e4s[0] ); // Top-left quadrant + sc &= tuple->column( "E4_quad1", e4s[1] ); // Top-right quadrant + sc &= tuple->column( "E4_quad2", e4s[2] ); // Bottom-right quadrant + sc &= tuple->column( "E4_quad3", e4s[3] ); // Bottom-left quadrant + + // Add 25 cell energies (5x5 grid) + for (int i = 0; i < 25; ++i) { + const double energyValue = cellEnergies[i]; + sc &= tuple->column( "ClE" + std::to_string(i+1), energyValue ); + } + + // Add diagnostics + sc &= tuple->column( "nCellsIn5x5", nCellsIn5x5 ); + sc &= tuple->column( "maxCellEnergy", maxCellEnergy ); + sc &= tuple->column( "nClusterCells", static_cast(entries.size()) ); + sc &= tuple->column( "sumEnergyForShowerShape", sumEnergyForShowerShape ); + sc.andThen( [&] { return tuple->write(); } ).orElse( [&] { ++m_writefailed; } ).ignore(); } } @@ -232,9 +931,11 @@ namespace LHCb::Calo { auto bestMatch( const Cluster& cluster, const Range& range ) const { using Result = std::pair; return std::accumulate( - range.begin(), range.end(), Result{ nullptr, m_minMatchFraction }, [&]( Result current, const auto& match ) { + range.begin(), range.end(), Result{ nullptr, 0.0 }, [&]( Result current, const auto& match ) { const MCParticle* particle = match.to(); if ( particle->momentum().e() < m_minEnergy ) return current; + + // Apply PDGID filtering if specified if ( !m_PDGID.empty() && find( m_PDGID.begin(), m_PDGID.end(), particle->particleID().pid() ) == m_PDGID.end() ) return current; @@ -242,6 +943,7 @@ namespace LHCb::Calo { ( !particle->mother() || find( m_PDGIDMother.begin(), m_PDGIDMother.end(), particle->mother()->particleID().pid() ) == m_PDGIDMother.end() ) ) return current; + auto weight = match.weight(); // "weight" of the relation auto pars = RetrieveParameters( cluster ); float matchFraction = weight / particle->momentum().e() * weight / pars[2]; @@ -252,9 +954,17 @@ namespace LHCb::Calo { // properties Gaudi::Property> m_PDGID{ this, "PDGID", {} }; Gaudi::Property> m_PDGIDMother{ this, "PDGIDMother", {} }; - Gaudi::Property m_minMatchFraction{ this, "minMatchFraction", 0.0 }; Gaudi::Property m_minEnergy{ this, "minEnergy", 0.0 }; - Gaudi::Property m_tuplePrefix{ this, "tuplePrefix", "matchedClusters" }; + Gaudi::Property m_tuplePrefix{ this, "tuplePrefix", "clusters" }; + + // New properties for cluster selection + Gaudi::Property m_PUMode{ this, "PUMode", false }; + Gaudi::Property m_signalMode{ this, "signalMode", false }; + Gaudi::Property m_photonMode{ this, "photonMode", false }; + Gaudi::Property m_hadronMode{ this, "hadronMode", false }; + Gaudi::Property m_debugWeights{ this, "debugWeights", false }; + Gaudi::Property m_puMaxFraction{ this, "puMaxFraction", 0.8f }; // PU: normalized weight < 0.8 + Gaudi::Property m_signalMinMatchFraction{ this, "signalMinMatchFraction", 0.1f }; // Signal: match fraction > 0.1 // counters mutable Gaudi::Accumulators::BinomialCounter m_counterMatched{ this, "Matched" }; @@ -266,4 +976,4 @@ namespace LHCb::Calo { DECLARE_COMPONENT_WITH_ID( ClusterResolution, "CaloFutureClusterResolution" ) DECLARE_COMPONENT_WITH_ID( ClusterResolution, "CaloFutureHypoResolution" ) -} // namespace LHCb::Calo +} // namespace LHCb::Calo \ No newline at end of file -- GitLab From d60462535dc744f8d935dd9bc5d49c19ff08ce16 Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Tue, 18 Nov 2025 20:03:32 +0000 Subject: [PATCH 2/2] pre-commit fixes patch generated by https://gitlab.cern.ch/lhcb/Lbcom/-/jobs/64811775 --- .../src/CaloClusterResolution.cpp | 921 +++++++++--------- 1 file changed, 436 insertions(+), 485 deletions(-) diff --git a/Calo/CaloCheckers/src/CaloClusterResolution.cpp b/Calo/CaloCheckers/src/CaloClusterResolution.cpp index df29d052f..c3ad9fbde 100644 --- a/Calo/CaloCheckers/src/CaloClusterResolution.cpp +++ b/Calo/CaloCheckers/src/CaloClusterResolution.cpp @@ -8,6 +8,7 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ +#include "CaloDet/CellParam.h" #include "CaloDet/DeCalorimeter.h" #include "CaloFutureUtils/CaloFuture2MC.h" #include "CaloFutureUtils/CaloFutureAlgUtils.h" @@ -17,21 +18,20 @@ #include "Event/CaloCluster.h" #include "Event/CaloClusters_v2.h" #include "Event/CaloDataFunctor.h" +#include "Event/CaloDigits_v2.h" #include "Event/CaloHypos_v2.h" #include "Event/MCParticle.h" -#include "Event/CaloDigits_v2.h" #include "GaudiAlg/GaudiTupleAlg.h" #include "LHCbAlgs/Consumer.h" #include "Relations/RelationWeighted1D.h" -#include "CaloDet/CellParam.h" -#include +#include #include +#include +#include #include -#include #include +#include #include -#include -#include using ICaloTable = LHCb::CaloFuture2MC::ClusterTable; namespace { @@ -173,7 +173,6 @@ namespace LHCb::Calo { * @date 2025-27-10 */ - template class ClusterResolution : public LHCb::Algorithm::Consumer> mcParticleWeights; // PDG, normalized weight - + { auto relationsRange = table.relations( currentSeed ); - + // First pass: calculate total weight sum for proper normalization double totalRawWeight = 0.0; - for ( const auto & rel : relationsRange ) { + for ( const auto& rel : relationsRange ) { double w = static_cast( rel.weight() ); totalRawWeight += w; } - + // Second pass: compute normalized weights with proper normalization - for ( const auto & rel : relationsRange ) { - const auto * mcPart = rel.to(); - double w = static_cast( rel.weight() ); - + for ( const auto& rel : relationsRange ) { + const auto* mcPart = rel.to(); + double w = static_cast( rel.weight() ); + // CORRECTED: Normalize by totalRawWeight, not cluster energy - double normalizedWeight = (totalRawWeight > 0) ? w / totalRawWeight : 0.0; - + double normalizedWeight = ( totalRawWeight > 0 ) ? w / totalRawWeight : 0.0; + relSumWeights += w; sumNormalizedWeights += normalizedWeight; - - if (normalizedWeight > maxNormalizedWeight) { - maxNormalizedWeight = normalizedWeight; - } - + + if ( normalizedWeight > maxNormalizedWeight ) { maxNormalizedWeight = normalizedWeight; } + if ( mcPart ) { relWeightedE += w * static_cast( mcPart->momentum().e() ); weightedMCParticleCount += normalizedWeight; - mcParticleWeights.emplace_back(mcPart->particleID().pid(), normalizedWeight); + mcParticleWeights.emplace_back( mcPart->particleID().pid(), normalizedWeight ); } } if ( matchParticle ) bestPdg = matchParticle->particleID().pid(); } // UPDATED CLUSTER SELECTION LOGIC - bool isPUCluster = false; + bool isPUCluster = false; bool isSignalCluster = false; bool isPhotonCluster = false; bool isHadronCluster = false; - int truePdgId = bestPdg; - - if (matchParticle) { - truePdgId = matchParticle->particleID().pid(); - - // PU CLUSTER: Normalized weight < 0.8 (cluster energy fraction from matched particle < 80%) - if (maxNormalizedWeight < m_puMaxFraction) { - isPUCluster = true; - } - - // SIGNAL CLUSTER: Match fraction > 0.1 (good quality match) - if (matchFraction > m_signalMinMatchFraction) { - isSignalCluster = true; - - // PHOTON CLUSTER: Photon or pi0 with good match - if (truePdgId == 22 || truePdgId == 111) { - isPhotonCluster = true; - } - - // HADRON CLUSTER: Hadron (|PDG| > 100) with good match - if (std::abs(truePdgId) > 100) { - isHadronCluster = true; - } - } + int truePdgId = bestPdg; + + if ( matchParticle ) { + truePdgId = matchParticle->particleID().pid(); + + // PU CLUSTER: Normalized weight < 0.8 (cluster energy fraction from matched particle < 80%) + if ( maxNormalizedWeight < m_puMaxFraction ) { isPUCluster = true; } + + // SIGNAL CLUSTER: Match fraction > 0.1 (good quality match) + if ( matchFraction > m_signalMinMatchFraction ) { + isSignalCluster = true; + + // PHOTON CLUSTER: Photon or pi0 with good match + if ( truePdgId == 22 || truePdgId == 111 ) { isPhotonCluster = true; } + + // HADRON CLUSTER: Hadron (|PDG| > 100) with good match + if ( std::abs( truePdgId ) > 100 ) { isHadronCluster = true; } + } } // CONDITIONAL SELECTION // Skip if we're in PU mode and this is not a PU cluster - if (m_PUMode && !isPUCluster) { - continue; - } - - // Skip if we're in signal mode and this is not a signal cluster - if (m_signalMode && !isSignalCluster) { - continue; - } - + if ( m_PUMode && !isPUCluster ) { continue; } + + // Skip if we're in signal mode and this is not a signal cluster + if ( m_signalMode && !isSignalCluster ) { continue; } + // Skip if we're in photon mode and this is not a photon cluster - if (m_photonMode && !isPhotonCluster) { - continue; - } - + if ( m_photonMode && !isPhotonCluster ) { continue; } + // Skip if we're in hadron mode and this is not a hadron cluster - if (m_hadronMode && !isHadronCluster) { - continue; - } + if ( m_hadronMode && !isHadronCluster ) { continue; } // ========== DEBUG OUTPUT ========== - if (m_debugWeights) { - std::cout << "=================================================================" << std::endl; - std::cout << "DEBUG CLUSTER SELECTION:" << std::endl; - std::cout << " Cluster Seed: " << currentSeed.all() - << " (Energy: " << energyCluster << " MeV)" << std::endl; - std::cout << " Best Match PDG: " << truePdgId - << ", Match Fraction: " << matchFraction - << ", Max Normalized Weight: " << maxNormalizedWeight << std::endl; - std::cout << " Selection Flags - PU: " << isPUCluster - << ", Signal: " << isSignalCluster - << ", Photon: " << isPhotonCluster - << ", Hadron: " << isHadronCluster << std::endl; - - if (m_PUMode && isPUCluster) { - std::cout << " -> SELECTED AS PU CLUSTER (maxNormalizedWeight < " << m_puMaxFraction << ")" << std::endl; - } else if (m_signalMode && isSignalCluster) { - std::cout << " -> SELECTED AS SIGNAL CLUSTER (matchFraction > " << m_signalMinMatchFraction << ")" << std::endl; - } else if (m_photonMode && isPhotonCluster) { - std::cout << " -> SELECTED AS PHOTON CLUSTER (PDG=" << truePdgId << ")" << std::endl; - } else if (m_hadronMode && isHadronCluster) { - std::cout << " -> SELECTED AS HADRON CLUSTER (|PDG| > 100)" << std::endl; - } - std::cout << "=================================================================" << std::endl; + if ( m_debugWeights ) { + std::cout << "=================================================================" << std::endl; + std::cout << "DEBUG CLUSTER SELECTION:" << std::endl; + std::cout << " Cluster Seed: " << currentSeed.all() << " (Energy: " << energyCluster << " MeV)" << std::endl; + std::cout << " Best Match PDG: " << truePdgId << ", Match Fraction: " << matchFraction + << ", Max Normalized Weight: " << maxNormalizedWeight << std::endl; + std::cout << " Selection Flags - PU: " << isPUCluster << ", Signal: " << isSignalCluster + << ", Photon: " << isPhotonCluster << ", Hadron: " << isHadronCluster << std::endl; + + if ( m_PUMode && isPUCluster ) { + std::cout << " -> SELECTED AS PU CLUSTER (maxNormalizedWeight < " << m_puMaxFraction << ")" << std::endl; + } else if ( m_signalMode && isSignalCluster ) { + std::cout << " -> SELECTED AS SIGNAL CLUSTER (matchFraction > " << m_signalMinMatchFraction << ")" + << std::endl; + } else if ( m_photonMode && isPhotonCluster ) { + std::cout << " -> SELECTED AS PHOTON CLUSTER (PDG=" << truePdgId << ")" << std::endl; + } else if ( m_hadronMode && isHadronCluster ) { + std::cout << " -> SELECTED AS HADRON CLUSTER (|PDG| > 100)" << std::endl; + } + std::cout << "=================================================================" << std::endl; } - + // Check photon conversion (only for photon clusters) std::optional zConv; - if (isPhotonCluster && matchParticle) { - zConv = checkConversion( matchParticle, currentSeed, table ); - } - - auto vertex = matchParticle ? matchParticle->originVertex() : nullptr; - auto deltaZ = currentZPos - (vertex ? vertex->position().Z() : 0.0); - auto mom = matchParticle ? matchParticle->momentum() : Gaudi::LorentzVector(); - auto xExtrap = vertex ? vertex->position().X() + deltaZ * ( mom.Px() / (mom.Pz() != 0 ? mom.Pz() : 1.0) ) : 0.0; - auto yExtrap = vertex ? vertex->position().Y() + deltaZ * ( mom.Py() / (mom.Pz() != 0 ? mom.Pz() : 1.0) ) : 0.0; - auto seedPos = seedParam.center(); // More accurate position - double CellSize = seedParam.size(); // Direct from CellParam - - auto energy = matchParticle ? matchParticle->momentum().E() : 0.0; - float deltaEoverE = matchParticle ? ( energyCluster - energy ) / energy : -999.0; + if ( isPhotonCluster && matchParticle ) { zConv = checkConversion( matchParticle, currentSeed, table ); } + + auto vertex = matchParticle ? matchParticle->originVertex() : nullptr; + auto deltaZ = currentZPos - ( vertex ? vertex->position().Z() : 0.0 ); + auto mom = matchParticle ? matchParticle->momentum() : Gaudi::LorentzVector(); + auto xExtrap = + vertex ? vertex->position().X() + deltaZ * ( mom.Px() / ( mom.Pz() != 0 ? mom.Pz() : 1.0 ) ) : 0.0; + auto yExtrap = + vertex ? vertex->position().Y() + deltaZ * ( mom.Py() / ( mom.Pz() != 0 ? mom.Pz() : 1.0 ) ) : 0.0; + auto seedPos = seedParam.center(); // More accurate position + double CellSize = seedParam.size(); // Direct from CellParam + + auto energy = matchParticle ? matchParticle->momentum().E() : 0.0; + float deltaEoverE = matchParticle ? ( energyCluster - energy ) / energy : -999.0; bufferDeltaEOverE += deltaEoverE; bufferDeltaX += xExtrap - currentPars[0]; bufferDeltaY += yExtrap - currentPars[1]; // FIXED: Ηandle uncorrected energy retrieval - double uncorE = -999.0, uncorX = -999.0, uncorY = -999.0; + double uncorE = -999.0, uncorX = -999.0, uncorY = -999.0; StatusCode uncorStatus = RetrieveUncorrectedEXY( cluster, calo, uncorE, uncorX, uncorY ); if ( !uncorStatus.isSuccess() ) { uncorE = -999.0; @@ -370,313 +350,291 @@ namespace LHCb::Calo { // ========== CORRECTED PSEUDORAPIDITY CALCULATION ========== // eta = -ln(tan(theta/2)) where theta is the angle from beam axis - double eta = -999.0; - double etaMC = -999.0; - double cosTheta = -999.0; + double eta = -999.0; + double etaMC = -999.0; + double cosTheta = -999.0; double cosThetaMC = -999.0; // For MC truth: use particle momentum direction - if (matchParticle) { - const auto& momentum = matchParticle->momentum(); - double pz = momentum.Pz(); - double pt = std::hypot(momentum.Px(), momentum.Py()); - double p = momentum.P(); - - if (p > 0) { - cosThetaMC = pz / p; // cos(theta) = pz/p - - // Using the identity: eta = asinh(pz/pt) - // This is numerically more stable than the log(tan) formula - // and equivalent to eta = -ln(tan(theta/2)) - if (pt > 0) { - etaMC = std::asinh(pz / pt); - } else if (pz > 0) { - etaMC = 10.0; // Large positive eta for forward particles - } else if (pz < 0) { - etaMC = -10.0; // Large negative eta for backward particles - } + if ( matchParticle ) { + const auto& momentum = matchParticle->momentum(); + double pz = momentum.Pz(); + double pt = std::hypot( momentum.Px(), momentum.Py() ); + double p = momentum.P(); + + if ( p > 0 ) { + cosThetaMC = pz / p; // cos(theta) = pz/p + + // Using the identity: eta = asinh(pz/pt) + // This is numerically more stable than the log(tan) formula + // and equivalent to eta = -ln(tan(theta/2)) + if ( pt > 0 ) { + etaMC = std::asinh( pz / pt ); + } else if ( pz > 0 ) { + etaMC = 10.0; // Large positive eta for forward particles + } else if ( pz < 0 ) { + etaMC = -10.0; // Large negative eta for backward particles } + } } // For reconstructed clusters: use position direction as approximation { - double r = std::hypot(seedPos.x(), seedPos.y()); - double mag = std::hypot(r, currentZPos); - - if (mag > 0) { - cosTheta = currentZPos / mag; // cos(theta) = z/r - - // Using position-based approximation: eta = -ln(tan(theta/2)) - // where theta = acos(cosTheta) - double theta = std::acos(cosTheta); - double tanHalfTheta = std::tan(theta / 2.0); - - if (tanHalfTheta > 0) { - eta = -std::log(tanHalfTheta); - // Add sign based on z position - if (currentZPos < 0) eta = -eta; - } else { - eta = (currentZPos > 0) ? 10.0 : -10.0; - } + double r = std::hypot( seedPos.x(), seedPos.y() ); + double mag = std::hypot( r, currentZPos ); + + if ( mag > 0 ) { + cosTheta = currentZPos / mag; // cos(theta) = z/r + + // Using position-based approximation: eta = -ln(tan(theta/2)) + // where theta = acos(cosTheta) + double theta = std::acos( cosTheta ); + double tanHalfTheta = std::tan( theta / 2.0 ); + + if ( tanHalfTheta > 0 ) { + eta = -std::log( tanHalfTheta ); + // Add sign based on z position + if ( currentZPos < 0 ) eta = -eta; } else { - // Handle the case where cluster is on the beam axis - cosTheta = (currentZPos > 0) ? 1.0 : -1.0; - eta = (currentZPos > 0) ? 10.0 : -10.0; + eta = ( currentZPos > 0 ) ? 10.0 : -10.0; } + } else { + // Handle the case where cluster is on the beam axis + cosTheta = ( currentZPos > 0 ) ? 1.0 : -1.0; + eta = ( currentZPos > 0 ) ? 10.0 : -10.0; + } } // Compute transverse energy using corrected eta - double Et = (energyCluster > 0 && std::abs(eta) < 10.0) ? - energyCluster / std::cosh(eta) : 0.; + double Et = ( energyCluster > 0 && std::abs( eta ) < 10.0 ) ? energyCluster / std::cosh( eta ) : 0.; // Get mother particle ID int motherID = 0; - if (matchParticle && matchParticle->mother()) { - motherID = matchParticle->mother()->particleID().pid(); - } + if ( matchParticle && matchParticle->mother() ) { motherID = matchParticle->mother()->particleID().pid(); } // Cluster shape variables - double E19 = 0.; // seed energy / 3x3 energy - double E49 = 0.; // E4/E9 ratio (max 2x2 / 3x3) - CORRECTED DEFINITION - std::array cellEnergies = {0}; // For 5x5 grid - + double E19 = 0.; // seed energy / 3x3 energy + double E49 = 0.; // E4/E9 ratio (max 2x2 / 3x3) - CORRECTED DEFINITION + std::array cellEnergies = { 0 }; // For 5x5 grid + // Get entries for both Clusters and Hypotheses auto getEntries = [&] { - if constexpr (std::is_same_v) { - return cluster.clusters()[0].entries(); - } else { - return cluster.entries(); - } + if constexpr ( std::is_same_v ) { + return cluster.clusters()[0].entries(); + } else { + return cluster.entries(); + } }; const auto entries = getEntries(); - + // Extract seed components const int seedArea = currentSeed.area(); - const int seedRow = currentSeed.row(); - const int seedCol = currentSeed.col(); - + const int seedRow = currentSeed.row(); + const int seedCol = currentSeed.col(); + // Get the seed cell parameters - const auto& seedCell = calo.cellParam(currentSeed); + const auto& seedCell = calo.cellParam( currentSeed ); // Create a map of all cells in the cluster using row/col as key // FIXED: Use energy * fraction like in CaloFutureHypoEstimator.cpp std::map, double> cellEnergyMap; - double maxCellEnergy = 0.; - double seedCellEnergy = 0.; - - for (const auto& entry : entries) { - auto id = entry.cellID(); - // CORRECTED: Use energy * fraction like in HypoEstimator - double energy = entry.energy() * entry.fraction(); - cellEnergyMap[{id.row(), id.col()}] = energy; - - if (energy > maxCellEnergy) { - maxCellEnergy = energy; - } - if (id == currentSeed) { - seedCellEnergy = energy; - } + double maxCellEnergy = 0.; + double seedCellEnergy = 0.; + + for ( const auto& entry : entries ) { + auto id = entry.cellID(); + // CORRECTED: Use energy * fraction like in HypoEstimator + double energy = entry.energy() * entry.fraction(); + cellEnergyMap[{ id.row(), id.col() }] = energy; + + if ( energy > maxCellEnergy ) { maxCellEnergy = energy; } + if ( id == currentSeed ) { seedCellEnergy = energy; } } // Get seed cell energy from the map to ensure consistency - auto seedIt = cellEnergyMap.find({seedRow, seedCol}); - if (seedIt != cellEnergyMap.end()) { - seedCellEnergy = seedIt->second; - } + auto seedIt = cellEnergyMap.find( { seedRow, seedCol } ); + if ( seedIt != cellEnergyMap.end() ) { seedCellEnergy = seedIt->second; } // Create a map for the full 5x5 grid std::map, double> cellMap; - + // Populate cellMap for the full 5x5 grid - for (int dr = -2; dr <= 2; ++dr) { - for (int dc = -2; dc <= 2; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - - // Look for this cell in the cluster - auto it = cellEnergyMap.find({targetRow, targetCol}); - if (it != cellEnergyMap.end()) { - cellMap[{targetRow, targetCol}] = it->second; - } else { - // Cell not in cluster, set to 0 - cellMap[{targetRow, targetCol}] = 0.0; - } + for ( int dr = -2; dr <= 2; ++dr ) { + for ( int dc = -2; dc <= 2; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + + // Look for this cell in the cluster + auto it = cellEnergyMap.find( { targetRow, targetCol } ); + if ( it != cellEnergyMap.end() ) { + cellMap[{ targetRow, targetCol }] = it->second; + } else { + // Cell not in cluster, set to 0 + cellMap[{ targetRow, targetCol }] = 0.0; } + } } // Fill cellEnergies from cellMap - for (int dr = -2; dr <= 2; ++dr) { - for (int dc = -2; dc <= 2; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - int index = (dr + 2) * 5 + (dc + 2); - - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - cellEnergies[index] = it->second; - } - } + for ( int dr = -2; dr <= 2; ++dr ) { + for ( int dc = -2; dc <= 2; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + int index = ( dr + 2 ) * 5 + ( dc + 2 ); + + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { cellEnergies[index] = it->second; } + } } // +++++++++ CORRECTED NON-OVERLAPPING QUADRANT DEFINITIONS +++++++++ // Following the same logic as CaloFutureHypoEstimator.cpp - std::array e4s = {0.0, 0.0, 0.0, 0.0}; // Four non-overlapping 2x2 quadrants - double E9 = 0.0; // Total energy in 3x3 grid - double E4 = 0.0; // Maximum energy in any 2x2 quadrant + std::array e4s = { 0.0, 0.0, 0.0, 0.0 }; // Four non-overlapping 2x2 quadrants + double E9 = 0.0; // Total energy in 3x3 grid + double E4 = 0.0; // Maximum energy in any 2x2 quadrant // Loop over 3x3 grid around seed for E9 calculation - for (int dr = -1; dr <= 1; ++dr) { - for (int dc = -1; dc <= 1; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - double cellEnergy = it->second; - E9 += cellEnergy; - } + for ( int dr = -1; dr <= 1; ++dr ) { + for ( int dc = -1; dc <= 1; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { + double cellEnergy = it->second; + E9 += cellEnergy; } + } } // Now calculate the four non-overlapping 2x2 quadrants // Quadrant 0: top-left [seedRow-1 to seedRow, seedCol-1 to seedCol] - for (int dr = -1; dr <= 0; ++dr) { - for (int dc = -1; dc <= 0; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - e4s[0] += it->second; - } - } + for ( int dr = -1; dr <= 0; ++dr ) { + for ( int dc = -1; dc <= 0; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { e4s[0] += it->second; } + } } // Quadrant 1: top-right [seedRow-1 to seedRow, seedCol to seedCol+1] - for (int dr = -1; dr <= 0; ++dr) { - for (int dc = 0; dc <= 1; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - e4s[1] += it->second; - } - } + for ( int dr = -1; dr <= 0; ++dr ) { + for ( int dc = 0; dc <= 1; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { e4s[1] += it->second; } + } } // Quadrant 2: bottom-right [seedRow to seedRow+1, seedCol to seedCol+1] - for (int dr = 0; dr <= 1; ++dr) { - for (int dc = 0; dc <= 1; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - e4s[2] += it->second; - } - } + for ( int dr = 0; dr <= 1; ++dr ) { + for ( int dc = 0; dc <= 1; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { e4s[2] += it->second; } + } } // Quadrant 3: bottom-left [seedRow to seedRow+1, seedCol-1 to seedCol] - for (int dr = 0; dr <= 1; ++dr) { - for (int dc = -1; dc <= 0; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - auto it = cellMap.find({targetRow, targetCol}); - if (it != cellMap.end()) { - e4s[3] += it->second; - } - } + for ( int dr = 0; dr <= 1; ++dr ) { + for ( int dc = -1; dc <= 0; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find( { targetRow, targetCol } ); + if ( it != cellMap.end() ) { e4s[3] += it->second; } + } } // E4 is the MAXIMUM energy among the four non-overlapping 2x2 quadrants - E4 = *std::max_element(e4s.begin(), e4s.end()); + E4 = *std::max_element( e4s.begin(), e4s.end() ); // ++++++++++ CORRECTED E19 CALCULATION WITH ROBUST CHECKS ++++++++++ // CRITICAL FIX: Ensure consistent energy definitions // Use the same energy definition for seedCellEnergy as used in E9 calculation - auto seedCellIt = cellMap.find({seedRow, seedCol}); - if (seedCellIt != cellMap.end()) { - seedCellEnergy = seedCellIt->second; // Consistent energy definition + auto seedCellIt = cellMap.find( { seedRow, seedCol } ); + if ( seedCellIt != cellMap.end() ) { + seedCellEnergy = seedCellIt->second; // Consistent energy definition } else { - seedCellEnergy = 0.0; - if (m_debugWeights) { - std::cout << "WARNING: Seed cell not found in cellMap for E19 calculation!" << std::endl; - } + seedCellEnergy = 0.0; + if ( m_debugWeights ) { + std::cout << "WARNING: Seed cell not found in cellMap for E19 calculation!" << std::endl; + } } // Calculate E19 with robust checks - if (E9 > 1e-16 && seedCellEnergy > 0) { // Use energy threshold to avoid numerical issues - E19 = seedCellEnergy / E9; - - // CRITICAL: Additional physical sanity checks - if (E19 > 1.2) { - // This should never happen - indicates serious calculation error - if (m_debugWeights) { - std::cout << "WARNING: Unphysical E19 = " << E19 - << " (seedEnergy=" << seedCellEnergy - << ", E9=" << E9 - << ", clusterEnergy=" << energyCluster << ")" << std::endl; - } - E19 = 1e-49; // Clamp to maximum physical value + if ( E9 > 1e-16 && seedCellEnergy > 0 ) { // Use energy threshold to avoid numerical issues + E19 = seedCellEnergy / E9; + + // CRITICAL: Additional physical sanity checks + if ( E19 > 1.2 ) { + // This should never happen - indicates serious calculation error + if ( m_debugWeights ) { + std::cout << "WARNING: Unphysical E19 = " << E19 << " (seedEnergy=" << seedCellEnergy << ", E9=" << E9 + << ", clusterEnergy=" << energyCluster << ")" << std::endl; } - - // Ensure E19 is physically meaningful (0-1) - E19 = std::max(0.0, std::min(1.0, E19)); - + E19 = 1e-49; // Clamp to maximum physical value + } + + // Ensure E19 is physically meaningful (0-1) + E19 = std::max( 0.0, std::min( 1.0, E19 ) ); + } else { - E19 = 0.; - if (m_debugWeights && E9 <= 1e-6) { - std::cout << "WARNING: E9 too small for E19 calculation: " << E9 << std::endl; - } + E19 = 0.; + if ( m_debugWeights && E9 <= 1e-6 ) { + std::cout << "WARNING: E9 too small for E19 calculation: " << E9 << std::endl; + } } // Compute E49 ratio (E4/E9) where E4 = max 2x2 quadrant energy, E9 = 3x3 grid energy - if (E9 > 1e-16) { - E49 = E4 / E9; - // Ensure E49 is physically meaningful (0-1) - if (E49 > 1.2) { - if (m_debugWeights) { - std::cout << "WARNING: Unphysical E49 = " << E49 - << " (E4=" << E4 << ", E9=" << E9 << ")" << std::endl; - } - E49 = 1e-49; + if ( E9 > 1e-16 ) { + E49 = E4 / E9; + // Ensure E49 is physically meaningful (0-1) + if ( E49 > 1.2 ) { + if ( m_debugWeights ) { + std::cout << "WARNING: Unphysical E49 = " << E49 << " (E4=" << E4 << ", E9=" << E9 << ")" << std::endl; } - E49 = std::max(0.0, std::min(1.0, E49)); - } else { E49 = 1e-49; + } + E49 = std::max( 0.0, std::min( 1.0, E49 ) ); + } else { + E49 = 1e-49; } // E19 diagnostics for debugging - if (m_debugWeights && (E19 < 0.2 || E19 > 1.0)) { - std::cout << "E19 DIAGNOSTICS - Cluster: " << currentSeed.all() << std::endl; - std::cout << " E19 = " << E19 << std::endl; - std::cout << " seedCellEnergy = " << seedCellEnergy << std::endl; - std::cout << " E9 = " << E9 << std::endl; - std::cout << " energyCluster = " << energyCluster << std::endl; - std::cout << " nClusterCells = " << entries.size() << std::endl; - - // Print 3x3 grid energies - std::cout << " 3x3 Grid Energies:" << std::endl; - for (int dr = -1; dr <= 1; ++dr) { - for (int dc = -1; dc <= 1; ++dc) { - int targetRow = seedRow + dr; - int targetCol = seedCol + dc; - auto it = cellMap.find({targetRow, targetCol}); - double energy = (it != cellMap.end()) ? it->second : 0.0; - std::cout << " [" << dr << "," << dc << "]: " << energy; - if (dr == 0 && dc == 0) std::cout << " <-- SEED"; - std::cout << std::endl; - } + if ( m_debugWeights && ( E19 < 0.2 || E19 > 1.0 ) ) { + std::cout << "E19 DIAGNOSTICS - Cluster: " << currentSeed.all() << std::endl; + std::cout << " E19 = " << E19 << std::endl; + std::cout << " seedCellEnergy = " << seedCellEnergy << std::endl; + std::cout << " E9 = " << E9 << std::endl; + std::cout << " energyCluster = " << energyCluster << std::endl; + std::cout << " nClusterCells = " << entries.size() << std::endl; + + // Print 3x3 grid energies + std::cout << " 3x3 Grid Energies:" << std::endl; + for ( int dr = -1; dr <= 1; ++dr ) { + for ( int dc = -1; dc <= 1; ++dc ) { + int targetRow = seedRow + dr; + int targetCol = seedCol + dc; + auto it = cellMap.find( { targetRow, targetCol } ); + double energy = ( it != cellMap.end() ) ? it->second : 0.0; + std::cout << " [" << dr << "," << dc << "]: " << energy; + if ( dr == 0 && dc == 0 ) std::cout << " <-- SEED"; + std::cout << std::endl; } + } } // ++++++++++ CORRECTED SHOWER SHAPE CALCULATION WITH SAFEGUARDS ++++++++++ // Following the definition from the document: r2 = = S_XX + S_YY // where S_XX, S_YY, S_XY are the position spread matrix elements - double S_XX = 0.0; - double S_YY = 0.0; - double S_XY = 0.0; - double r4 = 0.0; // For calculation + double S_XX = 0.0; + double S_YY = 0.0; + double S_XY = 0.0; + double r4 = 0.0; // For calculation double sumEnergyForShowerShape = 0.0; // Use cluster center from parameters (should be energy-weighted) @@ -684,53 +642,53 @@ namespace LHCb::Calo { double y_c = currentPars[1]; // First pass: calculate S_XX, S_YY, S_XY and r4 - for (const auto& entry : entries) { - auto id = entry.cellID(); - // Use consistent energy definition: energy * fraction - double cellEnergy = entry.energy() * entry.fraction(); - - // CRITICAL FIX: Ensure cell energy is non-negative - if (cellEnergy < 0.0) { - if (m_debugWeights) { - std::cout << "WARNING: Negative cell energy detected: " << cellEnergy - << " for cell " << id.all() << std::endl; - } - cellEnergy = 0.0; // Set negative energies to zero + for ( const auto& entry : entries ) { + auto id = entry.cellID(); + // Use consistent energy definition: energy * fraction + double cellEnergy = entry.energy() * entry.fraction(); + + // CRITICAL FIX: Ensure cell energy is non-negative + if ( cellEnergy < 0.0 ) { + if ( m_debugWeights ) { + std::cout << "WARNING: Negative cell energy detected: " << cellEnergy << " for cell " << id.all() + << std::endl; } - - const auto& cell = calo.cellParam(id); - auto cellPos = cell.center(); - - const double dx = cellPos.x() - x_c; - const double dy = cellPos.y() - y_c; - const double r2_component = dx*dx + dy*dy; - const double r4_component = r2_component * r2_component; // r⁴ - - // Calculate matrix elements as defined in ShowerShape.pdf Eq. (1) and (2) - S_XX += cellEnergy * dx * dx; - S_YY += cellEnergy * dy * dy; - S_XY += cellEnergy * dx * dy; - r4 += cellEnergy * r4_component; - - sumEnergyForShowerShape += cellEnergy; + cellEnergy = 0.0; // Set negative energies to zero + } + + const auto& cell = calo.cellParam( id ); + auto cellPos = cell.center(); + + const double dx = cellPos.x() - x_c; + const double dy = cellPos.y() - y_c; + const double r2_component = dx * dx + dy * dy; + const double r4_component = r2_component * r2_component; // r⁴ + + // Calculate matrix elements as defined in ShowerShape.pdf Eq. (1) and (2) + S_XX += cellEnergy * dx * dx; + S_YY += cellEnergy * dy * dy; + S_XY += cellEnergy * dx * dy; + r4 += cellEnergy * r4_component; + + sumEnergyForShowerShape += cellEnergy; } // CRITICAL FIX: Proper normalization with safeguards - if (sumEnergyForShowerShape > 1e-6) { // Use small threshold to avoid floating point issues - S_XX /= sumEnergyForShowerShape; - S_YY /= sumEnergyForShowerShape; - S_XY /= sumEnergyForShowerShape; - r4 /= sumEnergyForShowerShape; + if ( sumEnergyForShowerShape > 1e-6 ) { // Use small threshold to avoid floating point issues + S_XX /= sumEnergyForShowerShape; + S_YY /= sumEnergyForShowerShape; + S_XY /= sumEnergyForShowerShape; + r4 /= sumEnergyForShowerShape; } else { - // If no energy, set all to zero - S_XX = 0.0; - S_YY = 0.0; - S_XY = 0.0; - r4 = 0.0; - if (m_debugWeights) { - std::cout << "WARNING: Very small or zero sumEnergyForShowerShape: " - << sumEnergyForShowerShape << std::endl; - } + // If no energy, set all to zero + S_XX = 0.0; + S_YY = 0.0; + S_XY = 0.0; + r4 = 0.0; + if ( m_debugWeights ) { + std::cout << "WARNING: Very small or zero sumEnergyForShowerShape: " << sumEnergyForShowerShape + << std::endl; + } } // Compute the shower shape variables: @@ -739,112 +697,104 @@ namespace LHCb::Calo { double r2 = S_XX + S_YY; // CRITICAL FIX: Ensure r2 is non-negative - if (r2 < 0.0) { - if (m_debugWeights) { - std::cout << "WARNING: Negative r2 detected: " << r2 - << ", S_XX: " << S_XX << ", S_YY: " << S_YY - << ", sumEnergy: " << sumEnergyForShowerShape << std::endl; - } - r2 = 0.0; // Force to zero if negative due to numerical issues + if ( r2 < 0.0 ) { + if ( m_debugWeights ) { + std::cout << "WARNING: Negative r2 detected: " << r2 << ", S_XX: " << S_XX << ", S_YY: " << S_YY + << ", sumEnergy: " << sumEnergyForShowerShape << std::endl; + } + r2 = 0.0; // Force to zero if negative due to numerical issues } // II. r2r4 - r2r4 = 1 - (2 / ) double r2r4 = 0.0; - if (r4 > 1e-12 && r2 > 1e-12) { // Avoid division by very small numbers - double ratio = (r2 * r2) / r4; - // Ensure ratio is physically meaningful (0-1) - if (ratio >= 0.0 && ratio <= 1.0) { - r2r4 = 1.0 - ratio; - } else { - r2r4 = 0.0; // Default for unphysical cases - if (m_debugWeights && ratio > 1.0) { - std::cout << "WARNING: Unphysical r2r4 ratio: " << ratio - << " (r2=" << r2 << ", r4=" << r4 << ")" << std::endl; - } + if ( r4 > 1e-12 && r2 > 1e-12 ) { // Avoid division by very small numbers + double ratio = ( r2 * r2 ) / r4; + // Ensure ratio is physically meaningful (0-1) + if ( ratio >= 0.0 && ratio <= 1.0 ) { + r2r4 = 1.0 - ratio; + } else { + r2r4 = 0.0; // Default for unphysical cases + if ( m_debugWeights && ratio > 1.0 ) { + std::cout << "WARNING: Unphysical r2r4 ratio: " << ratio << " (r2=" << r2 << ", r4=" << r4 << ")" + << std::endl; } - // To ensure physical range - r2r4 = std::max(0.0, std::min(1.0, r2r4)); + } + // To ensure physical range + r2r4 = std::max( 0.0, std::min( 1.0, r2r4 ) ); } // III. k (kappa) - related to eigenvalue ratio - double kappa = 0.0; - double denominator = (S_XX + S_YY) * (S_XX + S_YY); - if (denominator > 1e-12) { // Avoid division by very small numbers - double determinant = S_XX * S_YY - S_XY * S_XY; - // Ensure determinant is non-negative for eigenvalue calculation - if (determinant >= 0.0) { - double kappa_squared = 1.0 - 4.0 * determinant / denominator; - // Ensure kappa_squared is non-negative - if (kappa_squared >= 0.0) { - kappa = std::sqrt(kappa_squared); - } - } - kappa = std::max(0.0, std::min(1.0, kappa)); + double kappa = 0.0; + double denominator = ( S_XX + S_YY ) * ( S_XX + S_YY ); + if ( denominator > 1e-12 ) { // Avoid division by very small numbers + double determinant = S_XX * S_YY - S_XY * S_XY; + // Ensure determinant is non-negative for eigenvalue calculation + if ( determinant >= 0.0 ) { + double kappa_squared = 1.0 - 4.0 * determinant / denominator; + // Ensure kappa_squared is non-negative + if ( kappa_squared >= 0.0 ) { kappa = std::sqrt( kappa_squared ); } + } + kappa = std::max( 0.0, std::min( 1.0, kappa ) ); } // IV. asym (asymmetry): correlation between X and Y double asym = 0.0; - if (S_XX > 1e-12 && S_YY > 1e-12) { // Avoid division by very small numbers - double correlation = S_XY / std::sqrt(S_XX * S_YY); - // Ensure correlation is within valid range [-1, 1] - if (std::isfinite(correlation)) { - asym = std::max(-1.0, std::min(1.0, correlation)); - } + if ( S_XX > 1e-12 && S_YY > 1e-12 ) { // Avoid division by very small numbers + double correlation = S_XY / std::sqrt( S_XX * S_YY ); + // Ensure correlation is within valid range [-1, 1] + if ( std::isfinite( correlation ) ) { asym = std::max( -1.0, std::min( 1.0, correlation ) ); } } // Also compute energy ratios mentioned in the document - double E_seed = seedCellEnergy; - double E_cl = energyCluster; - double E_seed_over_E_cl = (E_cl > 0) ? E_seed / E_cl : 0.0; + double E_seed = seedCellEnergy; + double E_cl = energyCluster; + double E_seed_over_E_cl = ( E_cl > 0 ) ? E_seed / E_cl : 0.0; // For E_2nd, we need to find the second most energetic cell double E_2nd = 0.0; - for (const auto& entry : entries) { - auto id = entry.cellID(); - if (id == currentSeed) continue; // skip seed - - double cellEnergy = entry.energy() * entry.fraction(); - if (cellEnergy > E_2nd) { - E_2nd = cellEnergy; - } + for ( const auto& entry : entries ) { + auto id = entry.cellID(); + if ( id == currentSeed ) continue; // skip seed + + double cellEnergy = entry.energy() * entry.fraction(); + if ( cellEnergy > E_2nd ) { E_2nd = cellEnergy; } } - double E_seed_plus_E_2nd_over_E_cl = (E_cl > 0) ? (E_seed + E_2nd) / E_cl : 0.0; + double E_seed_plus_E_2nd_over_E_cl = ( E_cl > 0 ) ? ( E_seed + E_2nd ) / E_cl : 0.0; // Ensure energy ratios are physically meaningful - E_seed_over_E_cl = std::max(0.0, std::min(1.0, E_seed_over_E_cl)); - E_seed_plus_E_2nd_over_E_cl = std::max(0.0, std::min(1.0, E_seed_plus_E_2nd_over_E_cl)); + E_seed_over_E_cl = std::max( 0.0, std::min( 1.0, E_seed_over_E_cl ) ); + E_seed_plus_E_2nd_over_E_cl = std::max( 0.0, std::min( 1.0, E_seed_plus_E_2nd_over_E_cl ) ); // Debug output for shower shape diagnostics - if (m_debugWeights && (r2 < 0.0 || !std::isfinite(r2))) { - std::cout << "SHOWER SHAPE DIAGNOSTICS:" << std::endl; - std::cout << " Cluster Energy: " << energyCluster << std::endl; - std::cout << " Sum Energy for Shower Shape: " << sumEnergyForShowerShape << std::endl; - std::cout << " S_XX: " << S_XX << ", S_YY: " << S_YY << ", S_XY: " << S_XY << std::endl; - std::cout << " r2: " << r2 << std::endl; - std::cout << " Number of cells: " << entries.size() << std::endl; - - // Print first few cell energies for debugging - int count = 0; - for (const auto& entry : entries) { - if (count++ < 5) { // Print first 5 cells - double cellEnergy = entry.energy() * entry.fraction(); - std::cout << " Cell " << count << ": energy=" << cellEnergy - << ", ID=" << entry.cellID().all() << std::endl; - } + if ( m_debugWeights && ( r2 < 0.0 || !std::isfinite( r2 ) ) ) { + std::cout << "SHOWER SHAPE DIAGNOSTICS:" << std::endl; + std::cout << " Cluster Energy: " << energyCluster << std::endl; + std::cout << " Sum Energy for Shower Shape: " << sumEnergyForShowerShape << std::endl; + std::cout << " S_XX: " << S_XX << ", S_YY: " << S_YY << ", S_XY: " << S_XY << std::endl; + std::cout << " r2: " << r2 << std::endl; + std::cout << " Number of cells: " << entries.size() << std::endl; + + // Print first few cell energies for debugging + int count = 0; + for ( const auto& entry : entries ) { + if ( count++ < 5 ) { // Print first 5 cells + double cellEnergy = entry.energy() * entry.fraction(); + std::cout << " Cell " << count << ": energy=" << cellEnergy << ", ID=" << entry.cellID().all() + << std::endl; } + } } // Count how many cells are in the 5x5 grid - int nCellsIn5x5 = std::count_if(cellEnergies.begin(), cellEnergies.end(), - [](double e) { return e > 0; }); + int nCellsIn5x5 = std::count_if( cellEnergies.begin(), cellEnergies.end(), []( double e ) { return e > 0; } ); // ========== TUPLE FILLING ========== - //auto sc = tuple->column( "matchFraction", matchFraction ); + // auto sc = tuple->column( "matchFraction", matchFraction ); auto sc = tuple->column( "maxMatchFraction", matchFraction ); - sc &= tuple->column( "ClTrueE", static_cast(relWeightedE) ); - sc &= tuple->column( "ClTrueSumWeight", static_cast(relSumWeights) ); + sc &= tuple->column( "ClTrueE", static_cast( relWeightedE ) ); + sc &= tuple->column( "ClTrueSumWeight", static_cast( relSumWeights ) ); sc &= tuple->column( "ClTruePdg", bestPdg ); - sc &= tuple->column( "maxNormalizedWeight", static_cast(maxNormalizedWeight) ); + sc &= tuple->column( "maxNormalizedWeight", static_cast( maxNormalizedWeight ) ); sc &= tuple->column( "isPUCluster", isPUCluster ); sc &= tuple->column( "isSignalCluster", isSignalCluster ); sc &= tuple->column( "isPhotonCluster", isPhotonCluster ); @@ -873,24 +823,24 @@ namespace LHCb::Calo { sc &= tuple->column( "px", matchParticle ? matchParticle->momentum().Px() : -999.0 ); sc &= tuple->column( "py", matchParticle ? matchParticle->momentum().Py() : -999.0 ); sc &= tuple->column( "pz", matchParticle ? matchParticle->momentum().Pz() : -999.0 ); - + // ++++++++++ CORRECTED PSEUDORAPIDITY AND ANGLE VARIABLES ++++++++++ sc &= tuple->column( "Et", Et ); - sc &= tuple->column( "eta", eta ); // Cluster position-based eta - sc &= tuple->column( "etaMC", etaMC ); // MC truth momentum-based eta - sc &= tuple->column( "cosTheta", cosTheta ); // cos(theta) from cluster position - sc &= tuple->column( "cosThetaMC", cosThetaMC ); // cos(theta) from MC momentum + sc &= tuple->column( "eta", eta ); // Cluster position-based eta + sc &= tuple->column( "etaMC", etaMC ); // MC truth momentum-based eta + sc &= tuple->column( "cosTheta", cosTheta ); // cos(theta) from cluster position + sc &= tuple->column( "cosThetaMC", cosThetaMC ); // cos(theta) from MC momentum sc &= tuple->column( "MCMotherID", motherID ); sc &= tuple->column( "cellID", currentSeed.all() ); sc &= tuple->column( "ClusterE", energyCluster ); sc &= tuple->column( "E19", E19 ); - sc &= tuple->column( "E49", E49 ); // CORRECTED: E4/E9 ratio (max 2x2 / 3x3) - sc &= tuple->column( "E4", E4 ); // CORRECTED: Maximum energy in any 2x2 quadrant - sc &= tuple->column( "E9", E9 ); // Energy in 3x3 grid - sc &= tuple->column( "seedCellEnergy", seedCellEnergy ); // Seed cell energy - + sc &= tuple->column( "E49", E49 ); // CORRECTED: E4/E9 ratio (max 2x2 / 3x3) + sc &= tuple->column( "E4", E4 ); // CORRECTED: Maximum energy in any 2x2 quadrant + sc &= tuple->column( "E9", E9 ); // Energy in 3x3 grid + sc &= tuple->column( "seedCellEnergy", seedCellEnergy ); // Seed cell energy + // ********** UPDATED SHOWER SHAPE VARIABLES ********** - sc &= tuple->column( "ShowerShape_r2", r2 ); // Main lateral dispersion (r2) + sc &= tuple->column( "ShowerShape_r2", r2 ); // Main lateral dispersion (r2) sc &= tuple->column( "ShowerShape_r2r4", r2r4 ); sc &= tuple->column( "ShowerShape_kappa", kappa ); sc &= tuple->column( "ShowerShape_asym", asym ); @@ -900,26 +850,26 @@ namespace LHCb::Calo { sc &= tuple->column( "E_seed_over_E_cl", E_seed_over_E_cl ); sc &= tuple->column( "E_seed_plus_E_2nd_over_E_cl", E_seed_plus_E_2nd_over_E_cl ); sc &= tuple->column( "E_2nd", E_2nd ); - + // Keep backward compatibility for existing analyses - sc &= tuple->column( "lateral", r2 ); // Alias for ShowerShape_r2 - + sc &= tuple->column( "lateral", r2 ); // Alias for ShowerShape_r2 + // Add individual quadrant energies for detailed analysis - sc &= tuple->column( "E4_quad0", e4s[0] ); // Top-left quadrant - sc &= tuple->column( "E4_quad1", e4s[1] ); // Top-right quadrant - sc &= tuple->column( "E4_quad2", e4s[2] ); // Bottom-right quadrant - sc &= tuple->column( "E4_quad3", e4s[3] ); // Bottom-left quadrant - + sc &= tuple->column( "E4_quad0", e4s[0] ); // Top-left quadrant + sc &= tuple->column( "E4_quad1", e4s[1] ); // Top-right quadrant + sc &= tuple->column( "E4_quad2", e4s[2] ); // Bottom-right quadrant + sc &= tuple->column( "E4_quad3", e4s[3] ); // Bottom-left quadrant + // Add 25 cell energies (5x5 grid) - for (int i = 0; i < 25; ++i) { - const double energyValue = cellEnergies[i]; - sc &= tuple->column( "ClE" + std::to_string(i+1), energyValue ); + for ( int i = 0; i < 25; ++i ) { + const double energyValue = cellEnergies[i]; + sc &= tuple->column( "ClE" + std::to_string( i + 1 ), energyValue ); } - + // Add diagnostics sc &= tuple->column( "nCellsIn5x5", nCellsIn5x5 ); sc &= tuple->column( "maxCellEnergy", maxCellEnergy ); - sc &= tuple->column( "nClusterCells", static_cast(entries.size()) ); + sc &= tuple->column( "nClusterCells", static_cast( entries.size() ) ); sc &= tuple->column( "sumEnergyForShowerShape", sumEnergyForShowerShape ); sc.andThen( [&] { return tuple->write(); } ).orElse( [&] { ++m_writefailed; } ).ignore(); @@ -934,7 +884,7 @@ namespace LHCb::Calo { range.begin(), range.end(), Result{ nullptr, 0.0 }, [&]( Result current, const auto& match ) { const MCParticle* particle = match.to(); if ( particle->momentum().e() < m_minEnergy ) return current; - + // Apply PDGID filtering if specified if ( !m_PDGID.empty() && find( m_PDGID.begin(), m_PDGID.end(), particle->particleID().pid() ) == m_PDGID.end() ) @@ -943,7 +893,7 @@ namespace LHCb::Calo { ( !particle->mother() || find( m_PDGIDMother.begin(), m_PDGIDMother.end(), particle->mother()->particleID().pid() ) == m_PDGIDMother.end() ) ) return current; - + auto weight = match.weight(); // "weight" of the relation auto pars = RetrieveParameters( cluster ); float matchFraction = weight / particle->momentum().e() * weight / pars[2]; @@ -956,15 +906,16 @@ namespace LHCb::Calo { Gaudi::Property> m_PDGIDMother{ this, "PDGIDMother", {} }; Gaudi::Property m_minEnergy{ this, "minEnergy", 0.0 }; Gaudi::Property m_tuplePrefix{ this, "tuplePrefix", "clusters" }; - + // New properties for cluster selection - Gaudi::Property m_PUMode{ this, "PUMode", false }; - Gaudi::Property m_signalMode{ this, "signalMode", false }; - Gaudi::Property m_photonMode{ this, "photonMode", false }; - Gaudi::Property m_hadronMode{ this, "hadronMode", false }; - Gaudi::Property m_debugWeights{ this, "debugWeights", false }; - Gaudi::Property m_puMaxFraction{ this, "puMaxFraction", 0.8f }; // PU: normalized weight < 0.8 - Gaudi::Property m_signalMinMatchFraction{ this, "signalMinMatchFraction", 0.1f }; // Signal: match fraction > 0.1 + Gaudi::Property m_PUMode{ this, "PUMode", false }; + Gaudi::Property m_signalMode{ this, "signalMode", false }; + Gaudi::Property m_photonMode{ this, "photonMode", false }; + Gaudi::Property m_hadronMode{ this, "hadronMode", false }; + Gaudi::Property m_debugWeights{ this, "debugWeights", false }; + Gaudi::Property m_puMaxFraction{ this, "puMaxFraction", 0.8f }; // PU: normalized weight < 0.8 + Gaudi::Property m_signalMinMatchFraction{ this, "signalMinMatchFraction", + 0.1f }; // Signal: match fraction > 0.1 // counters mutable Gaudi::Accumulators::BinomialCounter m_counterMatched{ this, "Matched" }; @@ -976,4 +927,4 @@ namespace LHCb::Calo { DECLARE_COMPONENT_WITH_ID( ClusterResolution, "CaloFutureClusterResolution" ) DECLARE_COMPONENT_WITH_ID( ClusterResolution, "CaloFutureHypoResolution" ) -} // namespace LHCb::Calo \ No newline at end of file +} // namespace LHCb::Calo -- GitLab