PWGHF:updated to calculate mass of Zee

Author: sashingo

PWGHF/Core/HfHelper.h

@@ -17,10 +17,8 @@
 #ifndef PWGHF_CORE_HFHELPER_H_
 #define PWGHF_CORE_HFHELPER_H_
 
-#include <vector>
-
-#include "Math/GenVector/Boost.h"
-#include "Math/Vector4D.h"
+#include <Math/GenVector/Boost.h>
+#include <Math/Vector4D.h>
 #include <TPDGCode.h>
 
 #include "CommonConstants/PhysicsConstants.h"
@@ -41,6 +39,7 @@ class HfHelper
 
   // 2-prong
 
+
   // D0(bar) → π± K∓
 
   template <typename T>
@@ -177,6 +176,16 @@ class HfHelper
     return candidate.cosThetaStar(std::array{o2::constants::physics::MassD0, o2::constants::physics::MassPiPlus}, o2::constants::physics::MassBPlus, 1);
   }
 
+  // Z → e+ e- decay
+
+  template <typename T>
+  auto invMassZtoEE(const T& e1, const T& e2)
+  {
+    auto arr1 = std::array{e1.px(), e1.py(), e1.pz()};
+    auto arr2 = std::array{e2.px(), e2.py(), e2.pz()};
+    return RecoDecay::m(std::array{arr1, arr2}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron});
+  }
+
   // 3-prong
 
   // D± → π± K∓ π±

PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx

@@ -30,6 +30,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 +55,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,10 +74,25 @@ 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 ElectronCandidate {
+    float pt, eta, phi, energy;
+    ElectronCandidate(float p, float e, float ph, float en)
+      : pt(p), eta(e), phi(ph), energy(en) {}
+
+    // 運動量成分を計算する関数を追加
+    float px() const { return pt * std::cos(phi); }
+    float py() const { return pt * std::sin(phi); }
+    float pz() const { return pt * std::sinh(eta); }
+  };
+  std::vector<ElectronCandidate> selectedElectrons_Iso;
+  std::vector<ElectronCandidate> selectedElectrons_Ass;
+
   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>;
 
+  HfHelper hfHelper;
+
   // pp
   // using TrackEle = o2::soa::Filtered<o2::soa::Join<o2::aod::Tracks, o2::aod::FullTracks, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCEl, o2::aod::pidTOFEl>>;
 
@@ -118,6 +135,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 +163,8 @@ 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("hInvMassZee", "invariant mass for Z", kTH2F, {{axisPt}, {axisInvMassZ}});
+    registry.add("hInvMassDy", "invariant mass for DY", kTH2F, {{axisPt}, {axisInvMassDy}});
   }
   bool isIsolatedCluster(const o2::aod::EMCALCluster& cluster,
                          const SelectedClusters& clusters)
@@ -247,6 +268,16 @@ 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) {
+        selectedElectrons_Ass.emplace_back(
+          track.pt(),
+          track.eta(),
+          track.phi(),
+          energyTrk);
+      }
+
       // 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,12 @@ struct HfTaskElectronWeakBoson {
 
               if (isIsolated) {
                 registry.fill(HIST("hEopIsolation"), match.track_as<TrackEle>().pt(), eop);
+
+                selectedElectrons_Iso.emplace_back(
+                  match.track_as<TrackEle>().pt(),
+                  match.track_as<TrackEle>().eta(),
+                  match.track_as<TrackEle>().phi(),
+                  energyEmc);
               }
 
               if (isIsolatedTr) {
@@ -335,6 +374,26 @@ struct HfTaskElectronWeakBoson {
       }
 
     } // end of track loop
+
+    // calculate inv. mass
+    if (selectedElectrons_Iso.size() > 1) {
+      for (size_t i = 0; i < selectedElectrons_Iso.size() - 1; ++i) {
+        const auto& e1 = selectedElectrons_Iso[i];
+        for (size_t j = 0; j < selectedElectrons_Ass.size() - 1; ++j) {
+          const auto& e2 = selectedElectrons_Ass[j];
+
+          if (e1.px() == e2.px())
+            continue;
+          auto mass = hfHelper.invMassZtoEE(e1, e2);
+          double ptIso = sqrt(pow(e1.px(), 2) + pow(e1.py(), 2));
+          double ptAss = sqrt(pow(e2.px(), 2) + pow(e2.py(), 2));
+          registry.fill(HIST("hInvMassDy"), ptIso, mass);
+          if (ptAss < 20.0 && ptIso < 20.0)
+            continue;
+          registry.fill(HIST("hInvMassZee"), ptIso, mass);
+        }
+      }
+    } // end of inv. mass calculation
   }
 };