diff --git a/Calo/CaloCheckers/src/CaloClusterResolution.cpp b/Calo/CaloCheckers/src/CaloClusterResolution.cpp index 9967940bacc364d95ceeeee6a21354d1e00f31a5..c3ad9fbdef5e4e41253bad0565be2a56bf6adfc0 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,11 +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 "GaudiAlg/GaudiTupleAlg.h" #include "LHCbAlgs/Consumer.h" #include "Relations/RelationWeighted1D.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,31 @@ 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 +206,599 @@ namespace LHCb::Calo { const auto& [matchParticle, matchFraction] = bestMatch( cluster, table.relations( currentSeed ) ); bufferMatched += ( matchParticle != nullptr ); - if ( !matchParticle ) continue; - - auto currentZPos = RetrieveZPosition( cluster ); - 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; + + // Use CellParam for precise geometry information + const auto& seedParam = calo.cellParam( currentSeed ); + auto currentZPos = seedParam.z(); // More accurate z-position + auto currentPars = RetrieveParameters( cluster ); + 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 +809,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 +881,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 +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]; @@ -252,9 +904,18 @@ 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" };