[PWGHF] added a function to isolation e (#8786)

Author: sashingo

PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx

@@ -39,9 +39,9 @@ struct HfTaskElectronWeakBoson {
 
   // configurable parameters
   Configurable<int> nBinsPt{"nBinsPt", 100, "N bins in pt registry"};
-  Configurable<float> BinPtmax{"BinPtmax", 100.0, "maximum pt registry"};
+  Configurable<float> binPtmax{"binPtmax", 100.0, "maximum pt registry"};
   Configurable<int> nBinsE{"nBinsE", 100, "N bins in E registry"};
-  Configurable<float> BinEmax{"BinEmax", 100.0, "maximum E registry"};
+  Configurable<float> binEmax{"binEmax", 100.0, "maximum E registry"};
 
   Configurable<float> vtxZ{"vtxZ", 10.f, ""};
 
@@ -64,7 +64,10 @@ struct HfTaskElectronWeakBoson {
   Configurable<float> timeEmcMax{"timeEmcMax", +20., "Maximum EMCcluster timing"};
   Configurable<float> m02Min{"m02Min", 0.1, "Minimum M02"};
   Configurable<float> m02Max{"m02Max", 0.9, "Maximum M02"};
-  Configurable<float> rMatchMax{"rMatchMax", 0.1, "cluster - track matching cut"};
+  Configurable<float> rMatchMax{"rMatchMax", 0.05, "cluster - track matching cut"};
+
+  Configurable<float> rIsolation{"rIsolation", 0.3, "cone radius for isolation cut"};
+  Configurable<float> 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<aod::Collisions>::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<TrackEle>().trackPhiEmcal() - phiEmc;
             dPhi = RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI);
 
-            registry.fill(HIST("hMatchPhi"), phiEmc, match.track_as<TrackEle>().phi());
-            registry.fill(HIST("hMatchEta"), etaEmc, match.track_as<TrackEle>().eta());
+            registry.fill(HIST("hMatchPhi"), phiEmc, match.track_as<TrackEle>().trackPhiEmcal());
+            registry.fill(HIST("hMatchEta"), etaEmc, match.track_as<TrackEle>().trackEtaEmcal());
 
             double r = RecoDecay::sqrtSumOfSquares(dPhi, dEta);
             if (r < rMin) {
@@ -223,16 +262,23 @@ 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<SelectedClusters>();
+            bool isIsolated = isIsolatedCluster(cluster, emcClusters);
+
             double eop = energyEmc / match.track_as<TrackEle>().p();
             // LOG(info) << "E/p" << eop;
             registry.fill(HIST("hEopNsigTPC"), match.track_as<TrackEle>().tpcNSigmaEl(), eop);
             registry.fill(HIST("hM02"), match.track_as<TrackEle>().tpcNSigmaEl(), m02Emc);
             registry.fill(HIST("hM20"), match.track_as<TrackEle>().tpcNSigmaEl(), m20Emc);
             if (match.track_as<TrackEle>().tpcNSigmaEl() > nsigTpcMin && match.track_as<TrackEle>().tpcNSigmaEl() < nsigTpcMax) {
               registry.fill(HIST("hEop"), match.track_as<TrackEle>().pt(), eop);
+
+              if (isIsolated) {
+                registry.fill(HIST("hEopIsolation"), match.track_as<TrackEle>().pt(), eop);
+              }
             }
           }