From d8fa84026fa7763642fe0cd635245540d33f525c Mon Sep 17 00:00:00 2001 From: Archita-Dash Date: Fri, 20 Dec 2024 15:09:33 +0100 Subject: [PATCH 1/4] implementing hadronic correction --- .../emcalClusterHadronicCorrectionTask.h | 54 +++ PWGJE/TableProducer/CMakeLists.txt | 5 + .../emcalClusterHadronicCorrectionTask.cxx | 307 ++++++++++++++++++ 3 files changed, 366 insertions(+) create mode 100644 PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h create mode 100644 PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx diff --git a/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h b/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h new file mode 100644 index 00000000000..f6459d70138 --- /dev/null +++ b/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h @@ -0,0 +1,54 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// **Hadronic Correction in the EMCAL framework: to avoid the double counting of the charged particles' contribution in jets** +/// \author Archita Rani Dash + +#ifndef PWGJE_DATAMODEL_EMCALCLUSTERHADRONICCORRECTION_H_ +#define PWGJE_DATAMODEL_EMCALCLUSTERHADRONICCORRECTION_H_ + +#include +#include "Framework/AnalysisDataModel.h" +#include "EMCALClusterDefinition.h" + +namespace o2::aod +{ + +namespace emcalhadroniccorrection +{ +// four different columns for the 'Cluster Energies after Hadronic Correction' to make it more flexible +// 2 for the closest match and 2 for all matched tracks + +// cluster values +// DECLARE_SOA_COLUMN(HadCorrEnergy, hadcorrEnergy, float); //! cluster energy (GeV) after hadronic correction + +// hadronic corrected energy values +DECLARE_SOA_COLUMN(HadCorrOneTrack1, hadCorrOneTrack1, float); //! with hadronic correction fraction (100%) for one matched track +DECLARE_SOA_COLUMN(HadCorrOneTrack2, hadCorrOneTrack2, float); //! with hadronic correction fraction (70%) for one matched track - systematic studies +DECLARE_SOA_COLUMN(HadCorrAllTracks1, hadCorrAllTracks1, float); //! with hadronic correction fraction (100%) for all matched tracks +DECLARE_SOA_COLUMN(HadCorrAllTracks2, hadCorrAllTracks2, float); //! with hadronic correction fraction (70%) for all matched tracks - for systematic studies + + +} // namespace emcalhadroniccorrection + +//Table Definitions - define what needs to be written into the tables produced by this tableproducer task +DECLARE_SOA_TABLE(EmcalHCs, "AOD", "EMCALHCS", //! + o2::soa::Index<>, //! + emcalhadroniccorrection::HadCorrOneTrack1, // corrected cluster energy for 1 matched track (f = 100%) + emcalhadroniccorrection::HadCorrOneTrack2, // corrected cluster energy for 1 matched track (f = 70%) + emcalhadroniccorrection::HadCorrAllTracks1, // corrected cluster energy for all matched tracks (f = 100%) + emcalhadroniccorrection::HadCorrAllTracks2 // corrected cluster energy for all matched tracks (f = 70%) + ) + +using EmcalHC = EmcalHCs::iterator; + +} //namespace o2::aod +#endif // PWGJE_DATAMODEL_EMCALCLUSTERHADRONICCORRECTION_H_ diff --git a/PWGJE/TableProducer/CMakeLists.txt b/PWGJE/TableProducer/CMakeLists.txt index 4e8141566f3..3cb21f46f8e 100644 --- a/PWGJE/TableProducer/CMakeLists.txt +++ b/PWGJE/TableProducer/CMakeLists.txt @@ -110,3 +110,8 @@ o2physics_add_dpl_workflow(emcal-matchedtracks-writer SOURCES emcalMatchedTracksTask.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(emcal-cluster-hadronic-correction-task + SOURCES emcalClusterHadronicCorrectionTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2::EMCALBase O2::EMCALReconstruction + COMPONENT_NAME Analysis) diff --git a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx new file mode 100644 index 00000000000..1467369e7ef --- /dev/null +++ b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx @@ -0,0 +1,307 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// **Hadronic Correction in the EMCAL framework: to avoid the double counting of the charged particles' contribution in jets** +/// \author Archita Rani Dash + +#include +#include +#include +#include +#include +#include +#include + +#include "CCDB/BasicCCDBManager.h" +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoA.h" +#include "Framework/HistogramRegistry.h" + +#include "DetectorsBase/GeometryManager.h" + +#include "PWGJE/DataModel/EMCALClusters.h" +#include "PWGJE/DataModel/EMCALMatchedTracks.h" +#include "PWGJE/DataModel/EMCALMatchedCollisions.h" +#include "PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h" + +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/PIDResponse.h" +#include "DataFormatsEMCAL/Cell.h" +#include "DataFormatsEMCAL/Constants.h" +#include "DataFormatsEMCAL/AnalysisCluster.h" +#include "EMCALBase/Geometry.h" +#include "EMCALBase/ClusterFactory.h" +#include "EMCALBase/NonlinearityHandler.h" +#include "EMCALReconstruction/Clusterizer.h" +#include "PWGJE/Core/JetUtilities.h" +#include "TVector2.h" + +#include "CommonDataFormat/InteractionRecord.h" + + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using collisionEvSelIt = o2::soa::Join::iterator; +using selectedClusters = o2::soa::Filtered; +using myTracks = o2::soa::Filtered>; + + +struct EmcalClusterHadronicCorrectionTask { + Produces hadroniccorrectedclusters; + + HistogramRegistry registry; + //Configurables for Histogram Binning + Preslice perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId; // looking at clusterID column in the EMCALMatchedTracks for every cluster + + //define configurables here + Configurable mClusterDefinition{"clusterDefinition", 10, "cluster definition to be selected, e.g. 10=kV3Default"}; + Configurable minTime{"minTime", -25., "Minimum cluster time for time cut"}; + Configurable maxTime{"maxTime", 20., "Maximum cluster time for time cut"}; + Configurable minM02{"minM02", 0.1, "Minimum M02 for M02 cut"}; + Configurable maxM02{"maxM02", 0.9, "Maximum M02 for M02 cut"}; + Configurable minTrackPt{"minTrackPt", 0.15, "Minimum pT for tracks"}; + Configurable fHadCorr1{"HadCorr1", 1., "hadronic correction fraction for complete cluster energy subtraction for one matched track" }; //100% - default + Configurable fHadCorr2{"HadCorr2", 0.7, "hadronic correction fraction for systematic studies for one matched track"}; //70% + Configurable fHadCorralltrks1{"HadCorralltrks1", 1., "hadronic correction fraction for complete cluster energy subtraction for all matched tracks" }; //100% - all tracks + Configurable fHadCorralltrks2{"HadCorralltrks2", 0.7, "hadronic correction fraction for systematic studies for all matched tracks"}; //70% + Configurable minDEta{"minDEta", 0.01, "Minimum dEta between track and cluster"}; + Configurable minDPhi{"minDPhi", 0.01, "Minimum dPhi between track and cluster"}; + Configurable fConstantSubtractionValue{"ConstantSubtractionValue", 0.236, "Value to be used for constant MIP subtraction (only applicable if using constant subtraction in M02 scheme)"}; + + //pT-dependent track-matching configurables + Configurable Eta0{"eta0", 0.04, "Param 0 in eta for pt-dependent matching"}; + Configurable Eta1{"eta1", 0.010, "Param 1 in eta for pt-dependent matching"}; + Configurable Eta2{"eta2", 2.5, "Param 2 in eta for pt-dependent matching"}; + Configurable Phi0{"phi0", 0.09, "Param 0 in phi for pt-dependent matching"}; + Configurable Phi1{"phi1", 0.015, "Param 1 in phi for pt-dependent matching"}; + Configurable Phi2{"phi2", 2.0, "Param 2 in phi for pt-dependent matching"}; + Configurable PhiMatch{"phiMatch", 0.050, "phi match value in pp"}; + Configurable EtaMatch{"etaMatch", 0.025, "eta match value in pp"}; + + Configurable doHadCorrSyst{"doHadCorrSyst", false, "Do hadronic correction for systematic studies"}; + Configurable doMomDepMatching{"doMomDepMatching", true, "Do momentum dependent track matching"}; // to be always set to true in Run 3 + Configurable UseM02SubtractionScheme1{"UseM02SubtractionScheme1", false, "Flag to enable hadronic correction scheme using cluster M02 value for Ecluster1 and EclusterAll1"}; + Configurable UseM02SubtractionScheme2{"UseM02SubtractionScheme2", false, "Flag to enable hadronic correction scheme using cluster M02 value for Ecluster2 and EclusterAll2"}; + Configurable UseFraction1{"UseFraction1", false, "Fractional momentum subtraction for Ecluster1 and EclusterAll1"}; + Configurable UseFraction2{"UseFraction2", false, "Fractional momentum subtraction for Ecluster2 and EclusterAll2"}; + + void init(o2::framework::InitContext&) { + + + //Event histograms + registry.add("h_allcollisions", "Total events; event status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); + + //Matched-Cluster histograms + registry.add("h_matchedclusters", "Total matched clusters; cluster status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); + registry.add("h_ClsE", "; Cls E w/o correction (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_Ecluster1", "; Ecluster1 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_Ecluster2", "; Ecluster2 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_EclusterAll1", "; EclusterAll1 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_EclusterAll2", "; EclusterAll2 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); + registry.add("h_ClsTime", "Cluster time distribution of uncorrected cluster E; #it{t}_{cls} (ns); entries", {HistType::kTH1F, {{500, -250., 250.}}}); + registry.add("h_ClsM02", "Cluster M02 distribution of uncorrected cluster E; #it{M}_{02}; entries", {HistType::kTH1F, {{400, 0., 5.}}}); + registry.add("h2_ClsEvsNmatches", "Original cluster energy vs Nmatches; Cls E w/o correction (GeV); Nmatches", {HistType::kTH2F,{{350, 0., 350.},{100, -0.5, 21.}}}); + registry.add("h2_ClsEvsEcluster1", "; Cls E w/o correction (GeV); Ecluster1 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); + registry.add("h2_ClsEvsEcluster2", "; Cls E w/o correction (GeV); Ecluster2 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); + registry.add("h2_ClsEvsEclusterAll1", "; Cls E w/o correction (GeV); EclusterAll1 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); + registry.add("h2_ClsEvsEclusterAll2", "; Cls E w/o correction (GeV); EclusterAll2 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); + + //Matched-Track histograms + registry.add("h_matchedtracks", "Total matched tracks; track status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); + } + + Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == mClusterDefinition) && (o2::aod::emcalcluster::time >= minTime) && (o2::aod::emcalcluster::time <= maxTime) && (o2::aod::emcalcluster::m02 > minM02) && (o2::aod::emcalcluster::m02 < maxM02); + Filter trackSelection = (o2::aod::track::pt >= minTrackPt); + //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(collisionEvSelIt const&, selectedClusters const& clusters, o2::aod::EMCALMatchedTracks const& matchedtracks, myTracks const&) + { + registry.fill(HIST("h_allcollisions"), 1); + + // skip events with no clusters + if (clusters.size() == 0) { + return; + } + + //Looping over all clusters matched to the collision + for (const auto& cluster : clusters) { + registry.fill(HIST("h_matchedclusters"), 1); + + double Ecluster1; double Ecluster2; double EclusterAll1; double EclusterAll2; + Ecluster1 = Ecluster2 = EclusterAll1 = EclusterAll2 = cluster.energy(); + + registry.fill(HIST("h_ClsE"), cluster.energy()); + registry.fill(HIST("h_ClsM02"), cluster.m02()); + registry.fill(HIST("h_ClsTime"), cluster.time()); + //selecting ALL MATCHED TRACKS after slicing all entries in perClusterMatchedTracks by the cluster globalIndex + auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex()); + + 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 (tracksofcluster.size() == 0) { + // Use original cluster energy values, no subtraction needed. + registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), 0); + registry.fill(HIST("h_Ecluster1"), Ecluster1); + registry.fill(HIST("h_Ecluster2"), Ecluster2); + registry.fill(HIST("h_EclusterAll1"), EclusterAll1); + registry.fill(HIST("h_EclusterAll2"), EclusterAll2); + registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), Ecluster1); + registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), Ecluster2); + registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), EclusterAll1); + registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), EclusterAll2); + hadroniccorrectedclusters(Ecluster1, Ecluster2, EclusterAll1, EclusterAll2); + continue; + } + + //Looping over all matched tracks for the cluster + //Total number of matched tracks = 20 (hard-coded) + for (const auto& match : tracksofcluster) { + + double mom = abs(match.track_as().p()); + registry.fill(HIST("h_matchedtracks"), 1); + + //CASE 1: skip tracks with a very low pT + if (mom < 1e-6) { + 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 + double dEta = match.track_as().eta() - cluster.eta(); + double dPhi = TVector2::Phi_mpi_pi(match.track_as().phi() - cluster.phi()); + + // 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) && fabs(dEta) < trackEtaMax) { + if (Nmatches == 0) { + closestTrkP = mom; + } + totalTrkP += mom; + Nmatches++; + } + } else { + // Do fixed dEta/dPhi matching (non-pT dependent) + if (fabs(dEta) >= minDEta || 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 + Ecluster1 = subtractM02ClusterEnergy(cluster.m02(), Ecluster1, Nmatches, closestTrkP, fHadCorr1, UseFraction1); + EclusterAll1 = subtractM02ClusterEnergy(cluster.m02(), EclusterAll1, Nmatches, totalTrkP, fHadCorralltrks1, UseFraction1); + } else { + // Default energy subtraction (100% and 70%) + Ecluster1 = subtractClusterEnergy(Ecluster1, closestTrkP, fHadCorr1, Nmatches, UseFraction1); + EclusterAll1 = subtractClusterEnergy(EclusterAll1, totalTrkP, fHadCorralltrks1, Nmatches, UseFraction1); + } + + if (UseM02SubtractionScheme2) { + // Do M02-based correction if enabled + Ecluster2 = subtractM02ClusterEnergy(cluster.m02(), Ecluster2, Nmatches, closestTrkP, fHadCorr2, UseFraction2); + EclusterAll2 = subtractM02ClusterEnergy(cluster.m02(), EclusterAll2, Nmatches, totalTrkP, fHadCorralltrks2, UseFraction2); + } else { + // Default energy subtraction (100% and 70%) + Ecluster2 = subtractClusterEnergy(Ecluster2, closestTrkP, fHadCorr2, Nmatches, UseFraction2); + EclusterAll2 = subtractClusterEnergy(EclusterAll2, totalTrkP, fHadCorralltrks2, Nmatches, UseFraction2); + } + } + registry.fill(HIST("h_Ecluster1"), Ecluster1); + registry.fill(HIST("h_Ecluster2"), Ecluster2); + registry.fill(HIST("h_EclusterAll1"), EclusterAll1); + registry.fill(HIST("h_EclusterAll2"), EclusterAll2); + registry.fill(HIST("h2_ClsEvsEcluster1"), cluster.energy(), Ecluster1); + registry.fill(HIST("h2_ClsEvsEcluster2"), cluster.energy(), Ecluster2); + registry.fill(HIST("h2_ClsEvsEclusterAll1"), cluster.energy(), EclusterAll1); + registry.fill(HIST("h2_ClsEvsEclusterAll2"), cluster.energy(), EclusterAll2); + + // Fill the table with all four corrected energies + hadroniccorrectedclusters(Ecluster1, Ecluster2, EclusterAll1, EclusterAll2); + + } // End of cluster loop + } // process function ends + + // Helper function to prevent negative energy subtraction + double subtractClusterEnergy(double Ecluster, double mom, double correctionFactor, int Nmatches, bool UseFraction) { + double Ecorr = Ecluster; + + // if (UseConstantSubtractionValue) { + if(!UseFraction) { + Ecorr = Ecluster- fConstantSubtractionValue* Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) + } else { + Ecorr = Ecluster - correctionFactor * mom; // Fractional momentum subtraction + } + return (Ecorr < 0) ? 0 : Ecorr; + } + + // Helper function for M02-based energy subtraction (based on cluster-shape) + double subtractM02ClusterEnergy(double m02, double Ecluster, int Nmatches, double totalTrkP, double correctionFactor, bool UseFraction) { + double Ecorr = Ecluster; + + // For M02 in the single photon region, the signal is primarily: Single photons, single electrons, single MIPs + if (m02 > 0.1 && m02 < 0.4) { //circular clusters(electron/photon) + Ecorr = Ecluster; // Single electron, single MIP + } else if (m02 > 0.4) { + // Large M02 region (M02 > 0.4), more complex overlaps and hadronic showers. + // The signal is primarily: Single hadronic shower, photon-photon overlap, photon-MIP overlap, MIP-MIP overlap, + // MIP-hadronic shower overlap, hadronic shower - hadronic shower overlap) + + if (!UseFraction) { + Ecorr = Ecluster- fConstantSubtractionValue* Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) + } else { + Ecorr = Ecluster - correctionFactor * totalTrkP; // Fractional momentum subtraction + } + } + return (Ecorr < 0) ? 0 : Ecorr; // Prevent negative energy + } + + PROCESS_SWITCH(EmcalClusterHadronicCorrectionTask, processMatchedCollisions, "Process matched clusters from collision", true); +}; //end of struct + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"emcal-cluster-hadronic-correction-task"})}; } From 3f68b35d21e319cab8dca1ce799ecb148b7873a1 Mon Sep 17 00:00:00 2001 From: Archita-Dash Date: Fri, 20 Dec 2024 18:34:57 +0100 Subject: [PATCH 2/4] removing filtered clusters and instead using all the EMCAL clusters --- PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx index 1467369e7ef..2b6acd64f7f 100644 --- a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx @@ -54,7 +54,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using collisionEvSelIt = o2::soa::Join::iterator; -using selectedClusters = o2::soa::Filtered; +// using selectedClusters = o2::soa::Filtered; using myTracks = o2::soa::Filtered>; @@ -127,7 +127,7 @@ struct EmcalClusterHadronicCorrectionTask { //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(collisionEvSelIt const&, selectedClusters const& clusters, o2::aod::EMCALMatchedTracks const& matchedtracks, myTracks const&) + void processMatchedCollisions(collisionEvSelIt const&, o2::aod::EMCALClusters const& clusters, o2::aod::EMCALMatchedTracks const& matchedtracks, myTracks const&) { registry.fill(HIST("h_allcollisions"), 1); From f3cdfa816d61b88b3b13cce741ed89cfb0da8beb Mon Sep 17 00:00:00 2001 From: Archita-Dash Date: Fri, 20 Dec 2024 18:36:10 +0100 Subject: [PATCH 3/4] removing filtered clusters and instead using all the EMCAL clusters --- PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx index 2b6acd64f7f..43542d28919 100644 --- a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx @@ -54,7 +54,6 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using collisionEvSelIt = o2::soa::Join::iterator; -// using selectedClusters = o2::soa::Filtered; using myTracks = o2::soa::Filtered>; From e8a2c6b3cb60de934c91e94f3f540b5957e36d84 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 20 Dec 2024 17:38:52 +0000 Subject: [PATCH 4/4] Please consider the following formatting changes --- .../emcalClusterHadronicCorrectionTask.h | 29 +++-- .../emcalClusterHadronicCorrectionTask.cxx | 112 +++++++++--------- 2 files changed, 71 insertions(+), 70 deletions(-) diff --git a/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h b/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h index f6459d70138..42a6219b159 100644 --- a/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h +++ b/PWGJE/DataModel/emcalClusterHadronicCorrectionTask.h @@ -31,24 +31,23 @@ namespace emcalhadroniccorrection // DECLARE_SOA_COLUMN(HadCorrEnergy, hadcorrEnergy, float); //! cluster energy (GeV) after hadronic correction // hadronic corrected energy values -DECLARE_SOA_COLUMN(HadCorrOneTrack1, hadCorrOneTrack1, float); //! with hadronic correction fraction (100%) for one matched track -DECLARE_SOA_COLUMN(HadCorrOneTrack2, hadCorrOneTrack2, float); //! with hadronic correction fraction (70%) for one matched track - systematic studies -DECLARE_SOA_COLUMN(HadCorrAllTracks1, hadCorrAllTracks1, float); //! with hadronic correction fraction (100%) for all matched tracks -DECLARE_SOA_COLUMN(HadCorrAllTracks2, hadCorrAllTracks2, float); //! with hadronic correction fraction (70%) for all matched tracks - for systematic studies - +DECLARE_SOA_COLUMN(HadCorrOneTrack1, hadCorrOneTrack1, float); //! with hadronic correction fraction (100%) for one matched track +DECLARE_SOA_COLUMN(HadCorrOneTrack2, hadCorrOneTrack2, float); //! with hadronic correction fraction (70%) for one matched track - systematic studies +DECLARE_SOA_COLUMN(HadCorrAllTracks1, hadCorrAllTracks1, float); //! with hadronic correction fraction (100%) for all matched tracks +DECLARE_SOA_COLUMN(HadCorrAllTracks2, hadCorrAllTracks2, float); //! with hadronic correction fraction (70%) for all matched tracks - for systematic studies } // namespace emcalhadroniccorrection -//Table Definitions - define what needs to be written into the tables produced by this tableproducer task -DECLARE_SOA_TABLE(EmcalHCs, "AOD", "EMCALHCS", //! - o2::soa::Index<>, //! - emcalhadroniccorrection::HadCorrOneTrack1, // corrected cluster energy for 1 matched track (f = 100%) - emcalhadroniccorrection::HadCorrOneTrack2, // corrected cluster energy for 1 matched track (f = 70%) - emcalhadroniccorrection::HadCorrAllTracks1, // corrected cluster energy for all matched tracks (f = 100%) - emcalhadroniccorrection::HadCorrAllTracks2 // corrected cluster energy for all matched tracks (f = 70%) - ) +// Table Definitions - define what needs to be written into the tables produced by this tableproducer task +DECLARE_SOA_TABLE(EmcalHCs, "AOD", "EMCALHCS", //! + o2::soa::Index<>, //! + emcalhadroniccorrection::HadCorrOneTrack1, // corrected cluster energy for 1 matched track (f = 100%) + emcalhadroniccorrection::HadCorrOneTrack2, // corrected cluster energy for 1 matched track (f = 70%) + emcalhadroniccorrection::HadCorrAllTracks1, // corrected cluster energy for all matched tracks (f = 100%) + emcalhadroniccorrection::HadCorrAllTracks2 // corrected cluster energy for all matched tracks (f = 70%) +) using EmcalHC = EmcalHCs::iterator; -} //namespace o2::aod -#endif // PWGJE_DATAMODEL_EMCALCLUSTERHADRONICCORRECTION_H_ +} // namespace o2::aod +#endif // PWGJE_DATAMODEL_EMCALCLUSTERHADRONICCORRECTION_H_ diff --git a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx index 43542d28919..a3bf8324e69 100644 --- a/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx +++ b/PWGJE/TableProducer/emcalClusterHadronicCorrectionTask.cxx @@ -49,37 +49,35 @@ #include "CommonDataFormat/InteractionRecord.h" - using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; using collisionEvSelIt = o2::soa::Join::iterator; -using myTracks = o2::soa::Filtered>; - +using myTracks = o2::soa::Filtered>; struct EmcalClusterHadronicCorrectionTask { Produces hadroniccorrectedclusters; HistogramRegistry registry; - //Configurables for Histogram Binning - Preslice perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId; // looking at clusterID column in the EMCALMatchedTracks for every cluster + // Configurables for Histogram Binning + Preslice perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId; // looking at clusterID column in the EMCALMatchedTracks for every cluster - //define configurables here + // define configurables here Configurable mClusterDefinition{"clusterDefinition", 10, "cluster definition to be selected, e.g. 10=kV3Default"}; Configurable minTime{"minTime", -25., "Minimum cluster time for time cut"}; Configurable maxTime{"maxTime", 20., "Maximum cluster time for time cut"}; Configurable minM02{"minM02", 0.1, "Minimum M02 for M02 cut"}; Configurable maxM02{"maxM02", 0.9, "Maximum M02 for M02 cut"}; Configurable minTrackPt{"minTrackPt", 0.15, "Minimum pT for tracks"}; - Configurable fHadCorr1{"HadCorr1", 1., "hadronic correction fraction for complete cluster energy subtraction for one matched track" }; //100% - default - Configurable fHadCorr2{"HadCorr2", 0.7, "hadronic correction fraction for systematic studies for one matched track"}; //70% - Configurable fHadCorralltrks1{"HadCorralltrks1", 1., "hadronic correction fraction for complete cluster energy subtraction for all matched tracks" }; //100% - all tracks - Configurable fHadCorralltrks2{"HadCorralltrks2", 0.7, "hadronic correction fraction for systematic studies for all matched tracks"}; //70% + Configurable fHadCorr1{"HadCorr1", 1., "hadronic correction fraction for complete cluster energy subtraction for one matched track"}; // 100% - default + Configurable fHadCorr2{"HadCorr2", 0.7, "hadronic correction fraction for systematic studies for one matched track"}; // 70% + Configurable fHadCorralltrks1{"HadCorralltrks1", 1., "hadronic correction fraction for complete cluster energy subtraction for all matched tracks"}; // 100% - all tracks + Configurable fHadCorralltrks2{"HadCorralltrks2", 0.7, "hadronic correction fraction for systematic studies for all matched tracks"}; // 70% Configurable minDEta{"minDEta", 0.01, "Minimum dEta between track and cluster"}; Configurable minDPhi{"minDPhi", 0.01, "Minimum dPhi between track and cluster"}; Configurable fConstantSubtractionValue{"ConstantSubtractionValue", 0.236, "Value to be used for constant MIP subtraction (only applicable if using constant subtraction in M02 scheme)"}; - //pT-dependent track-matching configurables + // pT-dependent track-matching configurables Configurable Eta0{"eta0", 0.04, "Param 0 in eta for pt-dependent matching"}; Configurable Eta1{"eta1", 0.010, "Param 1 in eta for pt-dependent matching"}; Configurable Eta2{"eta2", 2.5, "Param 2 in eta for pt-dependent matching"}; @@ -96,13 +94,13 @@ struct EmcalClusterHadronicCorrectionTask { Configurable UseFraction1{"UseFraction1", false, "Fractional momentum subtraction for Ecluster1 and EclusterAll1"}; Configurable UseFraction2{"UseFraction2", false, "Fractional momentum subtraction for Ecluster2 and EclusterAll2"}; - void init(o2::framework::InitContext&) { - + void init(o2::framework::InitContext&) + { - //Event histograms + // Event histograms registry.add("h_allcollisions", "Total events; event status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); - //Matched-Cluster histograms + // Matched-Cluster histograms registry.add("h_matchedclusters", "Total matched clusters; cluster status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); registry.add("h_ClsE", "; Cls E w/o correction (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); registry.add("h_Ecluster1", "; Ecluster1 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); @@ -111,21 +109,21 @@ struct EmcalClusterHadronicCorrectionTask { registry.add("h_EclusterAll2", "; EclusterAll2 (GeV); entries", {HistType::kTH1F, {{350, 0., 350.}}}); registry.add("h_ClsTime", "Cluster time distribution of uncorrected cluster E; #it{t}_{cls} (ns); entries", {HistType::kTH1F, {{500, -250., 250.}}}); registry.add("h_ClsM02", "Cluster M02 distribution of uncorrected cluster E; #it{M}_{02}; entries", {HistType::kTH1F, {{400, 0., 5.}}}); - registry.add("h2_ClsEvsNmatches", "Original cluster energy vs Nmatches; Cls E w/o correction (GeV); Nmatches", {HistType::kTH2F,{{350, 0., 350.},{100, -0.5, 21.}}}); - registry.add("h2_ClsEvsEcluster1", "; Cls E w/o correction (GeV); Ecluster1 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); - registry.add("h2_ClsEvsEcluster2", "; Cls E w/o correction (GeV); Ecluster2 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); - registry.add("h2_ClsEvsEclusterAll1", "; Cls E w/o correction (GeV); EclusterAll1 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); - registry.add("h2_ClsEvsEclusterAll2", "; Cls E w/o correction (GeV); EclusterAll2 (GeV)", {HistType::kTH2F,{{350, 0., 350.},{350, 0., 350.}}}); + registry.add("h2_ClsEvsNmatches", "Original cluster energy vs Nmatches; Cls E w/o correction (GeV); Nmatches", {HistType::kTH2F, {{350, 0., 350.}, {100, -0.5, 21.}}}); + registry.add("h2_ClsEvsEcluster1", "; Cls E w/o correction (GeV); Ecluster1 (GeV)", {HistType::kTH2F, {{350, 0., 350.}, {350, 0., 350.}}}); + registry.add("h2_ClsEvsEcluster2", "; Cls E w/o correction (GeV); Ecluster2 (GeV)", {HistType::kTH2F, {{350, 0., 350.}, {350, 0., 350.}}}); + registry.add("h2_ClsEvsEclusterAll1", "; Cls E w/o correction (GeV); EclusterAll1 (GeV)", {HistType::kTH2F, {{350, 0., 350.}, {350, 0., 350.}}}); + registry.add("h2_ClsEvsEclusterAll2", "; Cls E w/o correction (GeV); EclusterAll2 (GeV)", {HistType::kTH2F, {{350, 0., 350.}, {350, 0., 350.}}}); - //Matched-Track histograms + // Matched-Track histograms registry.add("h_matchedtracks", "Total matched tracks; track status;entries", {HistType::kTH1F, {{1, 0.5, 1.5}}}); } Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == mClusterDefinition) && (o2::aod::emcalcluster::time >= minTime) && (o2::aod::emcalcluster::time <= maxTime) && (o2::aod::emcalcluster::m02 > minM02) && (o2::aod::emcalcluster::m02 < maxM02); Filter trackSelection = (o2::aod::track::pt >= minTrackPt); - //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 + // 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(collisionEvSelIt const&, o2::aod::EMCALClusters const& clusters, o2::aod::EMCALMatchedTracks const& matchedtracks, myTracks const&) { registry.fill(HIST("h_allcollisions"), 1); @@ -135,22 +133,25 @@ struct EmcalClusterHadronicCorrectionTask { return; } - //Looping over all clusters matched to the collision + // Looping over all clusters matched to the collision for (const auto& cluster : clusters) { registry.fill(HIST("h_matchedclusters"), 1); - double Ecluster1; double Ecluster2; double EclusterAll1; double EclusterAll2; + double Ecluster1; + double Ecluster2; + double EclusterAll1; + double EclusterAll2; Ecluster1 = Ecluster2 = EclusterAll1 = EclusterAll2 = cluster.energy(); registry.fill(HIST("h_ClsE"), cluster.energy()); registry.fill(HIST("h_ClsM02"), cluster.m02()); registry.fill(HIST("h_ClsTime"), cluster.time()); - //selecting ALL MATCHED TRACKS after slicing all entries in perClusterMatchedTracks by the cluster globalIndex + // selecting ALL MATCHED TRACKS after slicing all entries in perClusterMatchedTracks by the cluster globalIndex auto tracksofcluster = matchedtracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex()); int Nmatches = 0; // counter for closest matched track double closestTrkP = 0.0; // closest track momentum - double totalTrkP = 0.0; // total 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 @@ -175,22 +176,22 @@ struct EmcalClusterHadronicCorrectionTask { continue; } - //Looping over all matched tracks for the cluster - //Total number of matched tracks = 20 (hard-coded) + // Looping over all matched tracks for the cluster + // Total number of matched tracks = 20 (hard-coded) for (const auto& match : tracksofcluster) { double mom = abs(match.track_as().p()); registry.fill(HIST("h_matchedtracks"), 1); - //CASE 1: skip tracks with a very low pT + // CASE 1: skip tracks with a very low pT if (mom < 1e-6) { 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% + // 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 double dEta = match.track_as().eta() - cluster.eta(); @@ -199,12 +200,11 @@ struct EmcalClusterHadronicCorrectionTask { // 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); + // Do pT-dependent track matching + if (doMomDepMatching) { + auto trackEtaMax = funcPtDepEta.Eval(mom); auto trackPhiHigh = +funcPtDepPhi.Eval(mom); - auto trackPhiLow = -funcPtDepPhi.Eval(mom); + auto trackPhiLow = -funcPtDepPhi.Eval(mom); if ((dPhi < trackPhiHigh && dPhi > trackPhiLow) && fabs(dEta) < trackEtaMax) { if (Nmatches == 0) { @@ -221,13 +221,13 @@ struct EmcalClusterHadronicCorrectionTask { // If track passes fixed dEta/dPhi cuts, process it if (Nmatches == 0) { - closestTrkP = mom; // Closest track match + closestTrkP = mom; // Closest track match } - totalTrkP += mom; // Accumulate momentum - Nmatches++; // Count this match + totalTrkP += mom; // Accumulate momentum + Nmatches++; // Count this match } - } // End of track loop + } // End of track loop registry.fill(HIST("h2_ClsEvsNmatches"), cluster.energy(), Nmatches); if (Nmatches >= 1) { @@ -267,40 +267,42 @@ struct EmcalClusterHadronicCorrectionTask { } // process function ends // Helper function to prevent negative energy subtraction - double subtractClusterEnergy(double Ecluster, double mom, double correctionFactor, int Nmatches, bool UseFraction) { + double subtractClusterEnergy(double Ecluster, double mom, double correctionFactor, int Nmatches, bool UseFraction) + { double Ecorr = Ecluster; // if (UseConstantSubtractionValue) { - if(!UseFraction) { - Ecorr = Ecluster- fConstantSubtractionValue* Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) + if (!UseFraction) { + Ecorr = Ecluster - fConstantSubtractionValue * Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) } else { - Ecorr = Ecluster - correctionFactor * mom; // Fractional momentum subtraction + Ecorr = Ecluster - correctionFactor * mom; // Fractional momentum subtraction } return (Ecorr < 0) ? 0 : Ecorr; } // Helper function for M02-based energy subtraction (based on cluster-shape) - double subtractM02ClusterEnergy(double m02, double Ecluster, int Nmatches, double totalTrkP, double correctionFactor, bool UseFraction) { + double subtractM02ClusterEnergy(double m02, double Ecluster, int Nmatches, double totalTrkP, double correctionFactor, bool UseFraction) + { double Ecorr = Ecluster; // For M02 in the single photon region, the signal is primarily: Single photons, single electrons, single MIPs - if (m02 > 0.1 && m02 < 0.4) { //circular clusters(electron/photon) - Ecorr = Ecluster; // Single electron, single MIP + if (m02 > 0.1 && m02 < 0.4) { // circular clusters(electron/photon) + Ecorr = Ecluster; // Single electron, single MIP } else if (m02 > 0.4) { // Large M02 region (M02 > 0.4), more complex overlaps and hadronic showers. // The signal is primarily: Single hadronic shower, photon-photon overlap, photon-MIP overlap, MIP-MIP overlap, // MIP-hadronic shower overlap, hadronic shower - hadronic shower overlap) if (!UseFraction) { - Ecorr = Ecluster- fConstantSubtractionValue* Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) + Ecorr = Ecluster - fConstantSubtractionValue * Nmatches; // Use constant value for MIP-subtraction (regardless of the cluster-shape) } else { - Ecorr = Ecluster - correctionFactor * totalTrkP; // Fractional momentum subtraction + Ecorr = Ecluster - correctionFactor * totalTrkP; // Fractional momentum subtraction } } - return (Ecorr < 0) ? 0 : Ecorr; // Prevent negative energy + return (Ecorr < 0) ? 0 : Ecorr; // Prevent negative energy } PROCESS_SWITCH(EmcalClusterHadronicCorrectionTask, processMatchedCollisions, "Process matched clusters from collision", true); -}; //end of struct +}; // end of struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"emcal-cluster-hadronic-correction-task"})}; }