[PWGHF/HFL] Calculate mass of Zee

PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx

@@ -12,6 +12,7 @@
 /// \file taskElectronWeakBoson.cxx
 /// \brief task for WeakBoson (W/Z) based on electron in mid-rapidity
 /// \author S. Sakai & S. Ito (Univ. of Tsukuba)
+#include <vector>
 
 #include "Framework/runDataProcessing.h"
 #include "Framework/AnalysisTask.h"
@@ -30,6 +31,7 @@
 #include "Common/DataModel/PIDResponse.h"
 
 #include "PWGJE/DataModel/EMCALClusters.h"
+#include "PWGHF/Core/HfHelper.h"
 
 using namespace o2;
 using namespace o2::framework;
@@ -54,6 +56,7 @@ struct HfTaskElectronWeakBoson {
   Configurable<float> nclItsMin{"nclItsMin", 2.0f, "its # of cluster cut"};
   Configurable<float> nclTpcMin{"nclTpcMin", 100.0f, "tpc # if cluster cut"};
   Configurable<float> nclcrossTpcMin{"nclcrossTpcMin", 100.0f, "tpc # of crossedRows cut"};
+  Configurable<float> nsigTpcMinLose{"nsigTpcMinLose", -3.0, "tpc Nsig lose lower cut"};
   Configurable<float> nsigTpcMin{"nsigTpcMin", -1.0, "tpc Nsig lower cut"};
   Configurable<float> nsigTpcMax{"nsigTpcMax", 3.0, "tpc Nsig upper cut"};
 
@@ -72,6 +75,17 @@ struct HfTaskElectronWeakBoson {
   Configurable<float> energyIsolationMax{"energyIsolationMax", 0.1, "isolation cut on energy"};
   Configurable<int> trackIsolationMax{"trackIsolationMax", 3, "Maximum number of tracks in isolation cone"};
 
+  struct HfElectronCandidate {
+    float pt, eta, phi, energy;
+    int charge;
+    HfElectronCandidate(float ptr, float e, float ph, float en, int ch)
+      : pt(ptr), eta(e), phi(ph), energy(en), charge(ch) {}
+
+    int sign() const { return charge; }
+  };
+  std::vector<HfElectronCandidate> selectedElectronsIso;
+  std::vector<HfElectronCandidate> selectedElectronsAss;
+
   using SelectedClusters = o2::aod::EMCALClusters;
   // PbPb
   using TrackEle = o2::soa::Join<o2::aod::Tracks, o2::aod::FullTracks, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCFullEl>;
@@ -118,6 +132,8 @@ struct HfTaskElectronWeakBoson {
     const AxisSpec axisEMCtime{200, -100.0, 100, "EMC time"};
     const AxisSpec axisIsoEnergy{100, 0, 1, "Isolation energy(GeV/C)"};
     const AxisSpec axisIsoTrack{20, -0.5, 19.5, "Isolation Track"};
+    const AxisSpec axisInvMassZ{200, 0, 200, "M_{ee} (GeV/c^{2})"};
+    const AxisSpec axisInvMassDy{200, 0, 2, "M_{ee} (GeV/c^{2})"};
 
     // create registrygrams
     registry.add("hZvtx", "Z vertex", kTH1F, {axisZvtx});
@@ -144,6 +160,10 @@ struct HfTaskElectronWeakBoson {
     registry.add("hEMCtime", "EMC timing", kTH1F, {axisEMCtime});
     registry.add("hIsolationEnergy", "Isolation Energy", kTH2F, {{axisE}, {axisIsoEnergy}});
     registry.add("hIsolationTrack", "Isolation Track", kTH2F, {{axisE}, {axisIsoTrack}});
+    registry.add("hInvMassZeeLs", "invariant mass for Z LS pair", kTH2F, {{axisPt}, {axisInvMassZ}});
+    registry.add("hInvMassZeeUls", "invariant mass for Z ULS pair", kTH2F, {{axisPt}, {axisInvMassZ}});
+    registry.add("hInvMassDyLs", "invariant mass for DY LS pair", kTH2F, {{axisPt}, {axisInvMassDy}});
+    registry.add("hInvMassDyUls", "invariant mass for DY ULS pair", kTH2F, {{axisPt}, {axisInvMassDy}});
   }
   bool isIsolatedCluster(const o2::aod::EMCALCluster& cluster,
                          const SelectedClusters& clusters)
@@ -247,6 +267,17 @@ struct HfTaskElectronWeakBoson {
       registry.fill(HIST("hPt"), track.pt());
       registry.fill(HIST("hTPCNsigma"), track.p(), track.tpcNSigmaEl());
 
+      float energyTrk = 0.0;
+
+      if (track.tpcNSigmaEl() > nsigTpcMinLose && track.tpcNSigmaEl() < nsigTpcMax) {
+        selectedElectronsAss.emplace_back(
+          track.pt(),
+          track.eta(),
+          track.phi(),
+          energyTrk,
+          track.sign());
+      }
+
       // track - match
 
       //  continue;
@@ -279,6 +310,7 @@ struct HfTaskElectronWeakBoson {
           // LOG(info) << "tr phi0 = " << match.track_as<TrackEle>().phi();
           // LOG(info) << "tr phi1 = " << track.phi();
           // LOG(info) << "emc phi = " << phiEmc;
+
           if (nMatch == 0) {
             double dEta = match.track_as<TrackEle>().trackEtaEmcal() - etaEmc;
             double dPhi = match.track_as<TrackEle>().trackPhiEmcal() - phiEmc;
@@ -303,6 +335,7 @@ struct HfTaskElectronWeakBoson {
             const auto& cluster = match.emcalcluster_as<SelectedClusters>();
 
             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);
@@ -317,6 +350,13 @@ struct HfTaskElectronWeakBoson {
 
               if (isIsolated) {
                 registry.fill(HIST("hEopIsolation"), match.track_as<TrackEle>().pt(), eop);
+
+                selectedElectronsIso.emplace_back(
+                  match.track_as<TrackEle>().pt(),
+                  match.track_as<TrackEle>().eta(),
+                  match.track_as<TrackEle>().phi(),
+                  energyEmc,
+                  match.track_as<TrackEle>().sign());
               }
 
               if (isIsolatedTr) {
@@ -335,6 +375,38 @@ struct HfTaskElectronWeakBoson {
       }
 
     } // end of track loop
+
+    // calculate inv. mass
+    if (selectedElectronsIso.size() > 1) {
+      for (size_t i = 0; i < selectedElectronsIso.size(); i++) {
+        const auto& e1 = selectedElectronsIso[i];
+        for (size_t j = 0; j < selectedElectronsAss.size(); j++) {
+          const auto& e2 = selectedElectronsAss[j];
+
+          float ptIso = e1.pt;
+          float ptAss = e2.pt;
+          if (ptIso == ptAss)
+            continue;
+          auto arr1 = RecoDecayPtEtaPhi::pVector(e1.pt, e1.eta, e1.phi);
+          auto arr2 = RecoDecayPtEtaPhi::pVector(e2.pt, e2.eta, e2.phi);
+          double mass = RecoDecay::m(std::array{arr1, arr2}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron});
+          if (e1.sign() * e2.sign() > 0) {
+            registry.fill(HIST("hInvMassDyLs"), ptIso, mass);
+          } else {
+            registry.fill(HIST("hInvMassDyUls"), ptIso, mass);
+          }
+
+          if (ptAss < 20.0 && ptIso < 20.0)
+            continue;
+
+          if (e1.sign() * e2.sign() > 0) {
+            registry.fill(HIST("hInvMassZeeLs"), ptIso, mass);
+          } else {
+            registry.fill(HIST("hInvMassZeeUls"), ptIso, mass);
+          }
+        }
+      }
+    } // end of inv. mass calculation
   }
 };