change jet definition

Author: alcaliva

MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C

@@ -1,13 +1,3 @@
-#include "TFile.h"
-#include "TTree.h"
-#include "TMath.h"
-#include "fastjet/ClusterSequence.hh"
-#include <vector>
-#include <iostream>
-#include "DataFormatsMC/MCTrack.h"
-
-using namespace fastjet;
-
 int External() {
   std::string path{"o2sim_Kine.root"};
 
@@ -22,73 +12,16 @@ int External() {
     std::cerr << "Cannot find tree o2sim in file " << path << "\n";
     return 1;
   }
-
   std::vector<o2::MCTrack> *tracks{};
   tree->SetBranchAddress("MCTrack", &tracks);
 
-  // Jet parameters
-  const double ptJetThreshold = 10.0;
-  const double jetR = 0.4;
-  const int gapSize = 4; // 4 events auto-accepted, 5th needs strange-in-jet
-
-  // Xi and Omega PDG codes
-  const std::vector<int> pdgXiOmega = {3312, -3312, 3334, -3334};
-
-  Long64_t nEntries = tree->GetEntries();
-
-  for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
-    tree->GetEntry(iEntry);
-    if (!tracks || tracks->empty()) continue;
-
-    bool acceptEvent = false;
-
-    // Gap-trigger logic
-    if (iEntry % (gapSize + 1) < gapSize) {
-      // Accept event automatically
-      acceptEvent = true;
-      std::cout << "Gap-trigger accepted event " << iEntry << " (no Xi/Omega check)\n";
-    } else {
-      // Require Xi/Omega inside a jet
-      std::vector<PseudoJet> fjParticles;
-      std::vector<int> fjIndexMap;
-      for (size_t i = 0; i < tracks->size(); ++i) {
-        const auto &t = tracks->at(i);
-        if (t.GetPdgCode() == 0) continue;
-        if (std::abs(t.GetEta()) > 0.8) continue; // acceptance cut
-        fjParticles.emplace_back(t.GetPx(), t.GetPy(), t.GetPz(), t.GetE());
-        fjIndexMap.push_back(i);
-      }
-
-      if (!fjParticles.empty()) {
-        JetDefinition jetDef(antikt_algorithm, jetR);
-        ClusterSequence cs(fjParticles, jetDef);
-        std::vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets(ptJetThreshold));
-
-        for (auto &jet : jets) {
-          auto constituents = jet.constituents();
-          for (auto &c : constituents) {
-            int trackIndex = fjIndexMap[c.user_index()];
-            int pdg = tracks->at(trackIndex).GetPdgCode();
-            if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) {
-              acceptEvent = true;
-              std::cout << "Accepted event " << iEntry
-                        << ": Jet pt = " << jet.pt()
-                        << ", eta = " << jet.eta()
-                        << ", phi = " << jet.phi()
-                        << ", contains PDG " << pdg << "\n";
-              break;
-            }
-          }
-          if (acceptEvent) break;
-        }
-      }
-    }
-
-    if (acceptEvent) {
-      // Here you could break if you only want the first accepted event,
-      // or continue to process all events following the gap-trigger pattern
-    }
+  auto nXi = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 3312");
+  auto nOmega = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 3334");
+  auto nStr = nXi+nOmega;
+    
+  if (nStr == 0) {
+    std::cerr << "No event of interest\n";
+    return 1;
   }
-
   return 0;
 }

MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C

@@ -1,106 +1,126 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
 #include "FairGenerator.h"
 #include "FairPrimaryGenerator.h"
 #include "Generators/GeneratorPythia8.h"
-
 #include "Pythia8/Pythia.h"
-
-#include <fastjet/PseudoJet.hh>
-#include <fastjet/JetDefinition.hh>
-#include <fastjet/ClusterSequence.hh>
-
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include "TParticlePDG.h"
+#include "TRandom3.h"
+#include "TSystem.h"
+#include "TVector2.h"
+#include "fairlogger/Logger.h"
 #include <cmath>
-#include <vector>
+#include <fstream>
 #include <string>
-#include <algorithm>
-
+#include <vector>
 using namespace Pythia8;
-using namespace fastjet;
+#endif
 
-/// Pythia8 event generator for pp collisions
-/// Selection of events with Xi or Omega inside jets with pt > 10 GeV
-/// Jets built from physical primaries OR HF decay products
+/// Event generator using Pythia ropes for Ξ and Ω production in jets.
+/// Jets are defined as particles within a cone around the highest-pT particle
+/// (approximate jet axis) and include charged final-state particles that are
+/// either physical primaries or from heavy-flavor decays.
+/// Jets must be fully within the acceptance and have pT > 10 GeV.
+/// One event of interest is generated after 4 minimum-bias events.
 
-class GeneratorPythia8StrangeInJet : public o2::eventgen::GeneratorPythia8 {
+class GeneratorPythia8StrangenessInJet : public o2::eventgen::GeneratorPythia8 {
 public:
   /// Constructor
-  GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0,
-                              double jetR = 0.4,
-                              int gapSize = 4)
+  GeneratorPythia8StrangenessInJet(double ptJetMin = 10.0, double Rjet = 0.4, int gapSize = 4)
     : o2::eventgen::GeneratorPythia8(),
-      mPtJetThreshold(ptJetThreshold),
-      mJetR(jetR),
+      mptJetMin(ptJetMin),
+      mRjet(Rjet),
       mGapSize(gapSize)
   {
     fmt::printf(
-      ">> Pythia8 generator: Xi/Omega inside jets with ptJet > %.1f GeV, R = %.1f, gap = %d\n",
-      ptJetThreshold, jetR, gapSize);
+      ">> Pythia8 generator: Xi/Omega inside jets with pt_{jet} > %.1f GeV, R_{jet} = %.1f, gap = %d\n",
+      ptJetMin, Rjet, gapSize);
   }
-
-  ~GeneratorPythia8StrangeInJet() = default;
+  /// Destructor
+  ~GeneratorPythia8StrangenessInJet() = default;
 
   bool Init() override {
     addSubGenerator(0, "Pythia8 events with Xi/Omega inside jets");
     return o2::eventgen::GeneratorPythia8::Init();
   }
 
 protected:
-  /// Check if particle is physical primary OR from HF decay
-  bool isPhysicalPrimaryOrFromHF(const Pythia8::Particle& p,
-                                const Pythia8::Event& event)
+  /// Check if particle is physical primary or from HF decay
+  bool isPhysicalPrimaryOrFromHF(const Pythia8::Particle& p, const Pythia8::Event& event)
   {
+    // Select only final-state particles
     if (!p.isFinal()) {
       return false;
     }
 
-    const int absPdg = std::abs(p.id());
-
     // Particle species selection
-    const bool isAcceptedSpecies = (absPdg == 211 || absPdg == 321 || absPdg == 2212 || absPdg == 1000010020 || absPdg == 11 || absPdg == 13);
-
-    if (!isAcceptedSpecies) {
+    const int absPdg = std::abs(p.id());
+    if (absPdg!=211 && absPdg!=321 && absPdg!= 2212 && absPdg!=1000010020 && absPdg!=11 && absPdg!=13)
       return false;
-    }
 
     // Walk up ancestry
-    int motherIdx = p.mother1();
+    int motherId = p.mother1();
 
-    while (motherIdx > 0) {
-      const auto& mother = event[motherIdx];
+    while (motherId > 0) {
+
+      // Get mother
+      const auto& mother = event[motherId];
       const int absMotherPdg = std::abs(mother.id());
 
-      // Charm or beauty hadron → accept (HF decay)
+      // Check if particle is from HF decay
       if ((absMotherPdg / 100 == 4) || (absMotherPdg / 100 == 5) || (absMotherPdg / 1000 == 4) || (absMotherPdg / 1000 == 5)) {
         return true;
       }
 
-      // Weakly decaying hadron → reject (non-physical primary)
+      // Reject non-physical primary hadrons
       if (mother.isHadron() && mother.tau0() > 1.0) {
         return false;
       }
-      motherIdx = mother.mother1();
+      motherId = mother.mother1();
     }
     return true;
   }
 
+  // Compute delta phi
+  double getDeltaPhi(double a1, double a2)
+  {
+    double deltaPhi{0.0};
+    double phi1 = TVector2::Phi_0_2pi(a1);
+    double phi2 = TVector2::Phi_0_2pi(a2);
+    double diff = std::fabs(phi1 - phi2);
+
+    if (diff <= M_PI)
+      deltaPhi = diff;
+    if (diff > M_PI)
+      deltaPhi = 2.0 * M_PI - diff;
+
+    return deltaPhi;
+  }
+
   bool generateEvent() override {
-    fmt::printf(">> Generating event %lu\n", mGeneratedEvents);
+    fmt::printf(">> Generating event %d\n", mGeneratedEvents);
 
     bool genOk = false;
     int localCounter{0};
+    constexpr int kMaxTries{100000};
 
     // Accept mGapSize events unconditionally, then one triggered event
     if (mGeneratedEvents % (mGapSize + 1) < mGapSize) {
       genOk = GeneratorPythia8::generateEvent();
       fmt::printf(">> Gap-trigger accepted event (no strangeness check)\n");
     } else {
-      while (!genOk) {
+      while (!genOk && localCounter < kMaxTries) {
         if (GeneratorPythia8::generateEvent()) {
           genOk = selectEvent(mPythia.event);
         }
         localCounter++;
       }
-      fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n",
-                  localCounter);
+      if (!genOk) {
+        fmt::printf("Failed to generate triggered event after %d tries\n",kMaxTries);
+        return false;
+      }
+      fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n", localCounter);
     }
 
     notifySubGenerator(0);
@@ -109,56 +129,54 @@ protected:
   }
 
   bool selectEvent(Pythia8::Event &event) {
-    const std::vector<int> pdgXiOmega = {3312, -3312, 3334, -3334};
-    const double mpi = 0.1395704;
 
-    std::vector<fastjet::PseudoJet> fjParticles;
+    double etaJet{-999.0}, phiJet{-999.0}, ptMax{0.0};
+    std::vector<int> particleID;
+    bool containsXiOrOmega{false};
 
-    for (int i = 0; i < event.size(); ++i) {
+    for (int i = 0; i < event.size(); i++) {
       const auto& p = event[i];
+      if (std::abs(p.id()) == 3312 || std::abs(p.id()) == 3334)
+        containsXiOrOmega = true;
 
-      // --- Jet input selection ---
-      if (!p.isFinal()) continue;
-      if (!p.isCharged()) continue;
+      if (std::abs(p.eta()) > 0.8 || p.pT() < 0.1) continue;
       if (!isPhysicalPrimaryOrFromHF(p, event)) continue;
-      if (std::abs(p.eta()) > 0.8) continue;
+      particleID.emplace_back(i);
 
-      double pt = std::sqrt(p.px() * p.px() + p.py() * p.py());
-      if (pt < 0.1) continue;
-
-      double energy = std::sqrt(p.p() * p.p() + mpi * mpi);
-
-      fastjet::PseudoJet pj(p.px(), p.py(), p.pz(), energy);
-      pj.set_user_index(i);   // map back to Pythia index
-      fjParticles.push_back(pj);
+      if (p.pT() > ptMax) {
+        ptMax = p.pT();
+        etaJet = p.eta();
+        phiJet = p.phi();
+      }
     }
+    if (ptMax == 0.0)
+      return false;
+    if (std::abs(etaJet) + Rjet > 0.8)
+      return false;
+    if (!containsXiOrOmega)
+      return false;
 
-    if (fjParticles.empty()) return false;
-
-    fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, mJetR);
-    fastjet::ClusterSequence cs(fjParticles, jetDef);
-    auto jets = fastjet::sorted_by_pt(cs.inclusive_jets(mPtJetThreshold));
-
-    for (const auto& jet : jets) {
-      for (const auto& c : jet.constituents()) {
-        int idx = c.user_index();
-        int pdg = event[idx].id();
+    double ptJet{0.0};
+    for (int i = 0 ; i < particleID.size() ; i++) {
+      int id = particleID[i];
+      const auto& p = event[id];
 
-        if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) {
-          fmt::printf(
-            ">> Accepted jet: pt = %.2f, eta = %.2f, phi = %.2f, contains PDG %d\n",
-            jet.pt(), jet.eta(), jet.phi(), pdg);
-          return true;
-        }
+      double deltaEta = std::abs(p.eta() - etaJet);
+      double deltaPhi = getDeltaPhi(p.phi(),phiJet);
+      double deltaR = std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi);
+      if (deltaR < Rjet) {
+        ptJet += p.pT();
       }
     }
+    if (ptJet < ptJetMin)
+      return false;
 
-    return false;
+    return true;
   }
 
 private:
-  double   mPtJetThreshold{10.0};
-  double   mJetR{0.4};
+  double   mptJetMin{10.0};
+  double   mRjet{0.4};
   int      mGapSize{4};
   uint64_t mGeneratedEvents{0};
 };