From a3dcedde3ad07f51c8d1fbfabdef4220597419c8 Mon Sep 17 00:00:00 2001 From: nzardosh Date: Fri, 13 Feb 2026 11:19:43 +0100 Subject: [PATCH 1/2] separating UPC tables --- PWGJE/DataModel/JetReducedData.h | 12 +++++---- PWGJE/TableProducer/derivedDataProducer.cxx | 10 +++++--- PWGJE/TableProducer/derivedDataWriter.cxx | 28 ++++++++++++++------- 3 files changed, 32 insertions(+), 18 deletions(-) diff --git a/PWGJE/DataModel/JetReducedData.h b/PWGJE/DataModel/JetReducedData.h index 69d8bc18d5e..42e994551ef 100644 --- a/PWGJE/DataModel/JetReducedData.h +++ b/PWGJE/DataModel/JetReducedData.h @@ -139,11 +139,6 @@ DECLARE_SOA_TABLE_STAGED(JCollisions, "JCOLLISION", jcollision::CentFT0C, jcollision::CentFT0M, jcollision::CentFT0CVariant1, - jcollision::AmplitudesFV0, - jcollision::AmplitudesFT0A, - jcollision::AmplitudesFT0C, - jcollision::AmplitudesFDDA, - jcollision::AmplitudesFDDC, jcollision::HadronicRate, jcollision::TrackOccupancyInTimeRange, jcollision::Alias, @@ -154,6 +149,13 @@ DECLARE_SOA_TABLE_STAGED(JCollisions, "JCOLLISION", using JCollision = JCollisions::iterator; using StoredJCollision = StoredJCollisions::iterator; +DECLARE_SOA_TABLE_STAGED(JCollisionUPCs, "JCOLLISIONUPC", + jcollision::AmplitudesFV0, + jcollision::AmplitudesFT0A, + jcollision::AmplitudesFT0C, + jcollision::AmplitudesFDDA, + jcollision::AmplitudesFDDC); + DECLARE_SOA_TABLE_STAGED(JCollisionMcInfos, "JCOLLISIONMCINFO", jcollision::Weight, jcollision::GetSubGeneratorId); diff --git a/PWGJE/TableProducer/derivedDataProducer.cxx b/PWGJE/TableProducer/derivedDataProducer.cxx index c3b78cdac73..fde62e9ccf5 100644 --- a/PWGJE/TableProducer/derivedDataProducer.cxx +++ b/PWGJE/TableProducer/derivedDataProducer.cxx @@ -78,6 +78,7 @@ struct JetDerivedDataProducerTask { Produces jBCsTable; Produces jBCParentIndexTable; Produces jCollisionsTable; + Produces jCollisionUPCsTable; Produces jCollisionMcInfosTable; Produces jCollisionsParentIndexTable; Produces jCollisionsBunchCrossingIndexTable; @@ -323,9 +324,10 @@ struct JetDerivedDataProducerTask { amplitudesFDDC.clear(); } } + products.jCollisionUPCsTable(amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC); } - products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), collision.centFV0A(), -1.0, collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), collision.centFT0CVariant1(), amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC, hadronicRate, collision.trackOccupancyInTimeRange(), collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision, upcGapResult), collision.rct_raw(), triggerBit); // note change multFT0C to multFT0M when problems with multFT0A are fixed + products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), collision.centFV0A(), -1.0, collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), collision.centFT0CVariant1(), hadronicRate, collision.trackOccupancyInTimeRange(), collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision, upcGapResult), collision.rct_raw(), triggerBit); // note change multFT0C to multFT0M when problems with multFT0A are fixed products.jCollisionsParentIndexTable(collision.globalIndex()); products.jCollisionsBunchCrossingIndexTable(collision.bcId()); } @@ -339,7 +341,7 @@ struct JetDerivedDataProducerTask { triggerDecider.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), jetderiveddatautilities::JTriggerMasks); triggerBit = jetderiveddatautilities::setTriggerSelectionBit(triggerDecider.getTriggerOfInterestResults(bc.globalBC())); } - products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC, -1.0, -1, collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision), collision.rct_raw(), triggerBit); + products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1, collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision), collision.rct_raw(), triggerBit); products.jCollisionsParentIndexTable(collision.globalIndex()); products.jCollisionsBunchCrossingIndexTable(collision.bcId()); } @@ -347,7 +349,7 @@ struct JetDerivedDataProducerTask { void processCollisionsRun2(soa::Join::iterator const& collision) { - products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, collision.centRun2V0A(), collision.centRun2V0M(), -1.0, -1.0, -1.0, -1.0, amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC, 1.0, -1, collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision), collision.rct_raw(), 0); // note change multFT0C to multFT0M when problems with multFT0A are fixed + products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, collision.centRun2V0A(), collision.centRun2V0M(), -1.0, -1.0, -1.0, -1.0, 1.0, -1, collision.alias_raw(), jetderiveddatautilities::setEventSelectionBit(collision), collision.rct_raw(), 0); // note change multFT0C to multFT0M when problems with multFT0A are fixed products.jCollisionsParentIndexTable(collision.globalIndex()); products.jCollisionsBunchCrossingIndexTable(collision.bcId()); } @@ -355,7 +357,7 @@ struct JetDerivedDataProducerTask { void processCollisionsALICE3(aod::Collision const& collision) { - products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC, -1.0, -1, -1.0, 0, 0, 0); + products.jCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1, -1.0, 0, 0, 0); products.jCollisionsParentIndexTable(collision.globalIndex()); products.jCollisionsBunchCrossingIndexTable(-1); } diff --git a/PWGJE/TableProducer/derivedDataWriter.cxx b/PWGJE/TableProducer/derivedDataWriter.cxx index 0da868140ef..d6208242541 100644 --- a/PWGJE/TableProducer/derivedDataWriter.cxx +++ b/PWGJE/TableProducer/derivedDataWriter.cxx @@ -63,6 +63,7 @@ struct JetDerivedDataWriter { Produces storedJBCsTable; Produces storedJBCParentIndexTable; Produces storedJCollisionsTable; + Produces storedJCollisionUPCsTable; Produces storedJCollisionMcInfosTable; Produces storedJCollisionsParentIndexTable; Produces storedJCollisionsBunchCrossingIndexTable; @@ -463,12 +464,27 @@ struct JetDerivedDataWriter { collisionMapping.clear(); collisionMapping.resize(collisions.size(), -1); + for (auto const& collision : collisions) { + if (collision.isCollisionSelected()) { + products.storedJCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), collision.centFV0A(), collision.centFV0M(), collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), collision.centFT0CVariant1(), collision.hadronicRate(), collision.trackOccupancyInTimeRange(), collision.alias_raw(), collision.eventSel(), collision.rct_raw(), collision.triggerSel()); + collisionMapping[collision.globalIndex()] = products.storedJCollisionsTable.lastIndex(); + products.storedJCollisionMcInfosTable(collision.weight(), collision.getSubGeneratorId()); + products.storedJCollisionsParentIndexTable(collision.collisionId()); + if (doprocessBCs) { + products.storedJCollisionsBunchCrossingIndexTable(bcMapping[collision.bcId()]); + } + } + } + } + PROCESS_SWITCH(JetDerivedDataWriter, processColllisons, "write out output tables for collisions", true); + + void processUPCCollisionInfos(soa::Join const& collisions) + { std::vector amplitudesFV0; std::vector amplitudesFT0A; std::vector amplitudesFT0C; std::vector amplitudesFDDA; std::vector amplitudesFDDC; - for (auto const& collision : collisions) { if (collision.isCollisionSelected()) { amplitudesFV0.clear(); @@ -486,17 +502,11 @@ struct JetDerivedDataWriter { std::copy(amplitudesFT0CSpan.begin(), amplitudesFT0CSpan.end(), std::back_inserter(amplitudesFT0C)); std::copy(amplitudesFDDASpan.begin(), amplitudesFDDASpan.end(), std::back_inserter(amplitudesFDDA)); std::copy(amplitudesFDDCSpan.begin(), amplitudesFDDCSpan.end(), std::back_inserter(amplitudesFDDC)); - products.storedJCollisionsTable(collision.posX(), collision.posY(), collision.posZ(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), collision.centFV0A(), collision.centFV0M(), collision.centFT0A(), collision.centFT0C(), collision.centFT0M(), collision.centFT0CVariant1(), amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC, collision.hadronicRate(), collision.trackOccupancyInTimeRange(), collision.alias_raw(), collision.eventSel(), collision.rct_raw(), collision.triggerSel()); - collisionMapping[collision.globalIndex()] = products.storedJCollisionsTable.lastIndex(); - products.storedJCollisionMcInfosTable(collision.weight(), collision.getSubGeneratorId()); - products.storedJCollisionsParentIndexTable(collision.collisionId()); - if (doprocessBCs) { - products.storedJCollisionsBunchCrossingIndexTable(bcMapping[collision.bcId()]); - } + products.storedJCollisionUPCsTable(amplitudesFV0, amplitudesFT0A, amplitudesFT0C, amplitudesFDDA, amplitudesFDDC); } } } - PROCESS_SWITCH(JetDerivedDataWriter, processColllisons, "write out output tables for collisions", true); + PROCESS_SWITCH(JetDerivedDataWriter, processUPCCollisionInfos, "write out table for upc collision info", false); void processTracks(soa::Join const& collisions, soa::Join const& tracks) { From e4228e0576597a6d41eba1743cb2165eb459b6c3 Mon Sep 17 00:00:00 2001 From: nzardosh Date: Fri, 13 Feb 2026 12:07:11 +0100 Subject: [PATCH 2/2] adding possibility to select multiple cluster definitions for hadronic correction --- .../emcalClusterHadronicCorrectionTask.cxx | 260 +++++++++--------- 1 file changed, 137 insertions(+), 123 deletions(-) diff --git a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx index 588030e132c..fd2b8f84a98 100644 --- a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx @@ -45,7 +45,7 @@ struct EmcalClusterHadronicCorrectionTask { PresliceUnsorted perTrackMatchedTrack = aod::jemctrack::trackId; // define configurables here - Configurable clusterDefinitionS{"clusterDefinition", "kV3Default", "cluster definition to be selected, e.g. V3Default"}; + Configurable clusterDefinitions{"clusterDefinitions", "kV3Default", "cluster definitions to be selected, e.g. V3Default"}; Configurable minTime{"minTime", -25., "Minimum cluster time for time cut"}; Configurable maxTime{"maxTime", 20., "Maximum cluster time for time cut"}; @@ -75,6 +75,8 @@ struct EmcalClusterHadronicCorrectionTask { Configurable useFraction1{"useFraction1", false, "Fractional momentum subtraction for clusterE1 and clusterEAll1"}; Configurable useFraction2{"useFraction2", false, "Fractional momentum subtraction for clusterE2 and clusterEAll2"}; + std::vector clusterDefinitionsVec; + void init(o2::framework::InitContext&) { // Event histograms @@ -97,15 +99,20 @@ struct EmcalClusterHadronicCorrectionTask { // Matched-Track histograms registry.add("h_matchedtracks", "Total matched tracks; track status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); - } - aod::EMCALClusterDefinition clusterDefinition = aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionS.value); - Filter clusterDefinitionSelection = (aod::jcluster::definition == static_cast(clusterDefinition)); + std::string clusterDefinitionsString = clusterDefinitions.value; + size_t start = 0; + size_t end; + while ((end = clusterDefinitionsString.find(',', start)) != std::string::npos) { + clusterDefinitionsVec.push_back(static_cast(aod::emcalcluster::getClusterDefinitionFromString(clusterDefinitionsString.substr(start, end - start)))); + start = end + 1; + } + } // The matching of clusters and tracks is already centralised in the EMCAL framework. // One only needs to apply a filter on matched clusters // Here looping over all collisions matched to EMCAL clusters - void processMatchedCollisions(aod::JetCollision const&, soa::Filtered> const& clusters, aod::JEMCTracks const& emcTracks, aod::JetTracks const&) + void processMatchedCollisions(aod::JetCollision const&, soa::Join const& clusters, aod::JEMCTracks const& emcTracks, aod::JetTracks const&) { registry.fill(HIST("h_allcollisions"), 1); @@ -114,35 +121,132 @@ struct EmcalClusterHadronicCorrectionTask { return; } - // Looping over all clusters matched to the collision - for (const auto& cluster : clusters) { - registry.fill(HIST("h_matchedclusters"), 1); - - double clusterE1; - double clusterE2; - double clusterEAll1; - double clusterEAll2; - clusterE1 = clusterE2 = clusterEAll1 = clusterEAll2 = cluster.energy(); + // Looping over all cluster definitions + for (const auto& clusterDefinition : clusterDefinitionsVec) { + // Looping over all clusters matched to the collision + for (const auto& cluster : clusters) { + if (cluster.definition() != clusterDefinition) { + continue; // Skip clusters that do not match the current cluster definition + } - registry.fill(HIST("h_ClsE"), cluster.energy()); - registry.fill(HIST("h_ClsM02"), cluster.m02()); - registry.fill(HIST("h_ClsTime"), cluster.time()); + registry.fill(HIST("h_matchedclusters"), 1); + + double clusterE1; + double clusterE2; + double clusterEAll1; + double clusterEAll2; + clusterE1 = clusterE2 = clusterEAll1 = clusterEAll2 = cluster.energy(); + + registry.fill(HIST("h_ClsE"), cluster.energy()); + registry.fill(HIST("h_ClsM02"), cluster.m02()); + registry.fill(HIST("h_ClsTime"), cluster.time()); + + int nMatches = 0; // counter for closest matched track + double closestTrkP = 0.0; // closest track momentum + double totalTrkP = 0.0; // total track momentum + + // pT-dependent track-matching instead of PID based track-matching to be adapted from Run 2 - suggested by Markus Fasel + + TF1 funcPtDepEta("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])"); + funcPtDepEta.SetParameters(eta0, eta1, eta2); + TF1 funcPtDepPhi("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])"); + funcPtDepEta.SetParameters(phi0, phi1, phi2); + + // No matched tracks (trackless case) + if (cluster.matchedTracks().size() == 0) { + // Use original cluster energy values, no subtraction needed. + registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), 0); + registry.fill(HIST("h_Ecluster1"), clusterE1); + registry.fill(HIST("h_Ecluster2"), clusterE2); + registry.fill(HIST("h_EclusterAll1"), clusterEAll1); + registry.fill(HIST("h_EclusterAll2"), clusterEAll2); + registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), clusterE1); + registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2); + registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1); + registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2); + clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2); + continue; + } - int nMatches = 0; // counter for closest matched track - double closestTrkP = 0.0; // closest track momentum - double totalTrkP = 0.0; // total track momentum + // Looping over all matched tracks for the cluster + // Total number of matched tracks = 20 (hard-coded) + for (const auto& matchedTrack : cluster.matchedTracks_as()) { + if (matchedTrack.pt() < minTrackPt) { + continue; + } + double mom = std::abs(matchedTrack.p()); + registry.fill(HIST("h_matchedtracks"), 1); + + // CASE 1: skip tracks with a very low pT + constexpr double kMinMom = 1e-6; + if (mom < kMinMom) { + continue; + } // end CASE 1 + + // CASE 2: + // a) If one matched track -> 100% energy subtraction + // b) If more than one matched track -> 100% energy subtraction + // c) If you want to do systematic studies -> perform the above two checks a) and b), and then subtract 70% energy instead of 100% + + // Perform dEta/dPhi matching + auto emcTrack = (emcTracks.sliceBy(perTrackMatchedTrack, matchedTrack.globalIndex())).iteratorAt(0); + double dEta = emcTrack.etaDiff(); + double dPhi = emcTrack.phiDiff(); + + // Apply the eta and phi matching thresholds + // dEta and dPhi cut : ensures that the matched track is within the desired eta/phi window + + // Do pT-dependent track matching + if (doMomDepMatching) { + auto trackEtaMax = funcPtDepEta.Eval(mom); + auto trackPhiHigh = +funcPtDepPhi.Eval(mom); + auto trackPhiLow = -funcPtDepPhi.Eval(mom); + + if ((dPhi < trackPhiHigh && dPhi > trackPhiLow) && std::fabs(dEta) < trackEtaMax) { + if (nMatches == 0) { + closestTrkP = mom; + } + totalTrkP += mom; + nMatches++; + } + } else { + // Do fixed dEta/dPhi matching (non-pT dependent) + if (std::fabs(dEta) >= minDEta || std::fabs(dPhi) >= minDPhi) { + continue; // Skip this track if outside the fixed cut region + } - // pT-dependent track-matching instead of PID based track-matching to be adapted from Run 2 - suggested by Markus Fasel + // If track passes fixed dEta/dPhi cuts, process it + if (nMatches == 0) { + closestTrkP = mom; // Closest track match + } + totalTrkP += mom; // Accumulate momentum + nMatches++; // Count this match + } - TF1 funcPtDepEta("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])"); - funcPtDepEta.SetParameters(eta0, eta1, eta2); - TF1 funcPtDepPhi("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])"); - funcPtDepEta.SetParameters(phi0, phi1, phi2); + } // End of track loop + registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), nMatches); + + if (nMatches >= 1) { + if (useM02SubtractionScheme1) { + // Do M02-based correction if enabled + clusterE1 = subtractM02ClusterEnergy(cluster.m02(), clusterE1, nMatches, closestTrkP, hadCorr1, useFraction1); + clusterEAll1 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll1, nMatches, totalTrkP, hadCorralltrks1, useFraction1); + } else { + // Default energy subtraction (100% and 70%) + clusterE1 = subtractClusterEnergy(clusterE1, closestTrkP, hadCorr1, nMatches, useFraction1); + clusterEAll1 = subtractClusterEnergy(clusterEAll1, totalTrkP, hadCorralltrks1, nMatches, useFraction1); + } - // No matched tracks (trackless case) - if (cluster.matchedTracks().size() == 0) { - // Use original cluster energy values, no subtraction needed. - registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), 0); + if (useM02SubtractionScheme2) { + // Do M02-based correction if enabled + clusterE2 = subtractM02ClusterEnergy(cluster.m02(), clusterE2, nMatches, closestTrkP, hadCorr2, useFraction2); + clusterEAll2 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll2, nMatches, totalTrkP, hadCorralltrks2, useFraction2); + } else { + // Default energy subtraction (100% and 70%) + clusterE2 = subtractClusterEnergy(clusterE2, closestTrkP, hadCorr2, nMatches, useFraction2); + clusterEAll2 = subtractClusterEnergy(clusterEAll2, totalTrkP, hadCorralltrks2, nMatches, useFraction2); + } + } registry.fill(HIST("h_Ecluster1"), clusterE1); registry.fill(HIST("h_Ecluster2"), clusterE2); registry.fill(HIST("h_EclusterAll1"), clusterEAll1); @@ -151,102 +255,12 @@ struct EmcalClusterHadronicCorrectionTask { registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2); registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1); registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2); - clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2); - continue; - } - - // Looping over all matched tracks for the cluster - // Total number of matched tracks = 20 (hard-coded) - for (const auto& matchedTrack : cluster.matchedTracks_as()) { - if (matchedTrack.pt() < minTrackPt) { - continue; - } - double mom = std::abs(matchedTrack.p()); - registry.fill(HIST("h_matchedtracks"), 1); - - // CASE 1: skip tracks with a very low pT - constexpr double kMinMom = 1e-6; - if (mom < kMinMom) { - continue; - } // end CASE 1 - - // CASE 2: - // a) If one matched track -> 100% energy subtraction - // b) If more than one matched track -> 100% energy subtraction - // c) If you want to do systematic studies -> perform the above two checks a) and b), and then subtract 70% energy instead of 100% - - // Perform dEta/dPhi matching - auto emcTrack = (emcTracks.sliceBy(perTrackMatchedTrack, matchedTrack.globalIndex())).iteratorAt(0); - double dEta = emcTrack.etaDiff(); - double dPhi = emcTrack.phiDiff(); - - // Apply the eta and phi matching thresholds - // dEta and dPhi cut : ensures that the matched track is within the desired eta/phi window - - // Do pT-dependent track matching - if (doMomDepMatching) { - auto trackEtaMax = funcPtDepEta.Eval(mom); - auto trackPhiHigh = +funcPtDepPhi.Eval(mom); - auto trackPhiLow = -funcPtDepPhi.Eval(mom); - - if ((dPhi < trackPhiHigh && dPhi > trackPhiLow) && std::fabs(dEta) < trackEtaMax) { - if (nMatches == 0) { - closestTrkP = mom; - } - totalTrkP += mom; - nMatches++; - } - } else { - // Do fixed dEta/dPhi matching (non-pT dependent) - if (std::fabs(dEta) >= minDEta || std::fabs(dPhi) >= minDPhi) { - continue; // Skip this track if outside the fixed cut region - } - - // If track passes fixed dEta/dPhi cuts, process it - if (nMatches == 0) { - closestTrkP = mom; // Closest track match - } - totalTrkP += mom; // Accumulate momentum - nMatches++; // Count this match - } - } // End of track loop - registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), nMatches); - - if (nMatches >= 1) { - if (useM02SubtractionScheme1) { - // Do M02-based correction if enabled - clusterE1 = subtractM02ClusterEnergy(cluster.m02(), clusterE1, nMatches, closestTrkP, hadCorr1, useFraction1); - clusterEAll1 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll1, nMatches, totalTrkP, hadCorralltrks1, useFraction1); - } else { - // Default energy subtraction (100% and 70%) - clusterE1 = subtractClusterEnergy(clusterE1, closestTrkP, hadCorr1, nMatches, useFraction1); - clusterEAll1 = subtractClusterEnergy(clusterEAll1, totalTrkP, hadCorralltrks1, nMatches, useFraction1); - } + // Fill the table with all four corrected energies + clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2); - if (useM02SubtractionScheme2) { - // Do M02-based correction if enabled - clusterE2 = subtractM02ClusterEnergy(cluster.m02(), clusterE2, nMatches, closestTrkP, hadCorr2, useFraction2); - clusterEAll2 = subtractM02ClusterEnergy(cluster.m02(), clusterEAll2, nMatches, totalTrkP, hadCorralltrks2, useFraction2); - } else { - // Default energy subtraction (100% and 70%) - clusterE2 = subtractClusterEnergy(clusterE2, closestTrkP, hadCorr2, nMatches, useFraction2); - clusterEAll2 = subtractClusterEnergy(clusterEAll2, totalTrkP, hadCorralltrks2, nMatches, useFraction2); - } - } - registry.fill(HIST("h_Ecluster1"), clusterE1); - registry.fill(HIST("h_Ecluster2"), clusterE2); - registry.fill(HIST("h_EclusterAll1"), clusterEAll1); - registry.fill(HIST("h_EclusterAll2"), clusterEAll2); - registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), clusterE1); - registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), clusterE2); - registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), clusterEAll1); - registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), clusterEAll2); - - // Fill the table with all four corrected energies - clusterEnergyCorrectedTable(clusterE1, clusterE2, clusterEAll1, clusterEAll2); - - } // End of cluster loop + } // End of cluster loop + } // End of cluster definition loop } // process function ends PROCESS_SWITCH(EmcalClusterHadronicCorrectionTask, processMatchedCollisions, "hadronic correction", true);