diff --git a/PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx b/PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx index e86fe36813f..55f586dcbe7 100644 --- a/PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx +++ b/PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx @@ -39,9 +39,9 @@ struct HfTaskElectronWeakBoson { // configurable parameters Configurable nBinsPt{"nBinsPt", 100, "N bins in pt registry"}; - Configurable BinPtmax{"BinPtmax", 100.0, "maximum pt registry"}; + Configurable binPtmax{"binPtmax", 100.0, "maximum pt registry"}; Configurable nBinsE{"nBinsE", 100, "N bins in E registry"}; - Configurable BinEmax{"BinEmax", 100.0, "maximum E registry"}; + Configurable binEmax{"binEmax", 100.0, "maximum E registry"}; Configurable vtxZ{"vtxZ", 10.f, ""}; @@ -64,7 +64,10 @@ struct HfTaskElectronWeakBoson { Configurable timeEmcMax{"timeEmcMax", +20., "Maximum EMCcluster timing"}; Configurable m02Min{"m02Min", 0.1, "Minimum M02"}; Configurable m02Max{"m02Max", 0.9, "Maximum M02"}; - Configurable rMatchMax{"rMatchMax", 0.1, "cluster - track matching cut"}; + Configurable rMatchMax{"rMatchMax", 0.05, "cluster - track matching cut"}; + + Configurable rIsolation{"rIsolation", 0.3, "cone radius for isolation cut"}; + Configurable energyIsolationMax{"energyIsolationMax", 0.1, "isolation cut on energy"}; using SelectedClusters = o2::aod::EMCALClusters; // PbPb @@ -79,7 +82,7 @@ struct HfTaskElectronWeakBoson { Filter etafilter = (aod::track::eta < etaTrUp) && (aod::track::eta > etaTrLow); Filter dcaxyfilter = (nabs(aod::track::dcaXY) < dcaxyMax); - Filter filter_globalTr = requireGlobalTrackInFilter(); + Filter filterGlobalTr = requireGlobalTrackInFilter(); Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == clusterDefinition) && (o2::aod::emcalcluster::time >= timeEmcMin) && (o2::aod::emcalcluster::time <= timeEmcMax) && (o2::aod::emcalcluster::m02 > m02Min) && (o2::aod::emcalcluster::m02 < m02Max); @@ -98,9 +101,9 @@ struct HfTaskElectronWeakBoson { const AxisSpec axisZvtx{400, -20, 20, "Zvtx"}; const AxisSpec axisCounter{1, 0, 1, "events"}; const AxisSpec axisEta{200, -1.0, 1.0, "#eta"}; - const AxisSpec axisPt{nBinsPt, 0, BinPtmax, "p_{T}"}; + const AxisSpec axisPt{nBinsPt, 0, binPtmax, "p_{T}"}; const AxisSpec axisNsigma{100, -5, 5, "N#sigma"}; - const AxisSpec axisE{nBinsE, 0, BinEmax, "Energy"}; + const AxisSpec axisE{nBinsE, 0, binEmax, "Energy"}; const AxisSpec axisM02{100, 0, 1, "M02"}; const AxisSpec axisdPhi{200, -1, 1, "dPhi"}; const AxisSpec axisdEta{200, -1, 1, "dEta"}; @@ -110,6 +113,7 @@ struct HfTaskElectronWeakBoson { const AxisSpec axisCluster{100, 0.0, 200.0, "counts"}; const AxisSpec axisITSNCls{20, 0.0, 20, "counts"}; const AxisSpec axisEMCtime{200, -100.0, 100, "EMC time"}; + const AxisSpec axisIsoEnergy{100, 0, 1, "Isolation energy(GeV/C)"}; // create registrygrams registry.add("hZvtx", "Z vertex", kTH1F, {axisZvtx}); @@ -130,12 +134,47 @@ struct HfTaskElectronWeakBoson { registry.add("hMatchPhi", "Match in Phi", kTH2F, {{axisPhi}, {axisPhi}}); registry.add("hMatchEta", "Match in Eta", kTH2F, {{axisEta}, {axisEta}}); registry.add("hEop", "energy momentum match", kTH2F, {{axisPt}, {axisEop}}); + registry.add("hEopIsolation", "energy momentum match after isolation", kTH2F, {{axisPt}, {axisEop}}); registry.add("hEopNsigTPC", "Eop vs. Nsigma", kTH2F, {{axisNsigma}, {axisEop}}); registry.add("hEMCtime", "EMC timing", kTH1F, {axisEMCtime}); + registry.add("hIsolationEnergy", "Isolation Energy", kTH2F, {{axisE}, {axisIsoEnergy}}); + } + bool isIsolatedCluster(const o2::aod::EMCALCluster& cluster, + const SelectedClusters& clusters) + { + float energySum = 0.0; + float isoEnergy = 10.0; + float etaAssCluster = cluster.eta(); + float phiAssCluster = cluster.phi(); + + for (const auto& associateCluster : clusters) { + // Calculate angular distances + double dEta = associateCluster.eta() - etaAssCluster; + double dPhi = associateCluster.phi() - phiAssCluster; + + // Normalize φ difference + dPhi = RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI); + + // Calculate ΔR + double deltaR = std::sqrt(dEta * dEta + dPhi * dPhi); + + // Sum energy within isolation cone + if (deltaR < rIsolation) { + energySum += associateCluster.energy(); + } + } + + if (energySum > 0) { + isoEnergy = energySum / cluster.energy() - 1.0; + } + + registry.fill(HIST("hIsolationEnergy"), cluster.energy(), isoEnergy); + + return (isoEnergy < energyIsolationMax); } void process(soa::Filtered::iterator const& collision, - SelectedClusters const&, + SelectedClusters const& emcClusters, TrackEle const& tracks, o2::aod::EMCALMatchedTracks const& matchedtracks) { @@ -210,8 +249,8 @@ struct HfTaskElectronWeakBoson { double dPhi = match.track_as().trackPhiEmcal() - phiEmc; dPhi = RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI); - registry.fill(HIST("hMatchPhi"), phiEmc, match.track_as().phi()); - registry.fill(HIST("hMatchEta"), etaEmc, match.track_as().eta()); + registry.fill(HIST("hMatchPhi"), phiEmc, match.track_as().trackPhiEmcal()); + registry.fill(HIST("hMatchEta"), etaEmc, match.track_as().trackEtaEmcal()); double r = RecoDecay::sqrtSumOfSquares(dPhi, dEta); if (r < rMin) { @@ -223,9 +262,12 @@ struct HfTaskElectronWeakBoson { registry.fill(HIST("hEMCtime"), timeEmc); registry.fill(HIST("hEnergy"), energyEmc); - if (r < rMatchMax) + if (r > rMatchMax) continue; + const auto& cluster = match.emcalcluster_as(); + bool isIsolated = isIsolatedCluster(cluster, emcClusters); + double eop = energyEmc / match.track_as().p(); // LOG(info) << "E/p" << eop; registry.fill(HIST("hEopNsigTPC"), match.track_as().tpcNSigmaEl(), eop); @@ -233,6 +275,10 @@ struct HfTaskElectronWeakBoson { registry.fill(HIST("hM20"), match.track_as().tpcNSigmaEl(), m20Emc); if (match.track_as().tpcNSigmaEl() > nsigTpcMin && match.track_as().tpcNSigmaEl() < nsigTpcMax) { registry.fill(HIST("hEop"), match.track_as().pt(), eop); + + if (isIsolated) { + registry.fill(HIST("hEopIsolation"), match.track_as().pt(), eop); + } } }