Add hybrid genberator for antideuterons in jets

MC/config/PWGLF/hybrid/hybrid_coalescence.json

@@ -0,0 +1,24 @@
+
+{
+  "mode": "sequential",
+  "generators": [
+    {
+      "name": "pythia8",
+      "config": {
+        "config": "${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg",
+        "hooksFileName": "",
+        "hooksFuncName": "",
+        "includePartonEvent": false,
+        "particleFilter": "",
+        "verbose": 0
+      }
+    },
+    {
+      "name": "external",
+      "config": {
+        "iniFile": "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini"
+      }
+    }
+  ],
+  "fractions": [4, 1]
+}

MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini

@@ -0,0 +1,6 @@
+[GeneratorExternal]
+fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_antideuteron_injection.C
+funcName = generateAntideuterons(0.8)
+
+[GeneratorPythia8]
+config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_jet.cfg

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

@@ -0,0 +1,85 @@
+
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "SimulationDataFormat/MCTrack.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include <algorithm>
+#include <iostream>
+#include <string>
+#include <vector>
+#endif // !defined(__CLING__) || defined(__ROOTCLING__)
+
+int External()
+{
+  std::string path{"o2sim_Kine.root"};
+
+  std::vector<int> injectedPDGs = {
+    -1000010020 // anti-deuteron
+  };
+
+  TFile file(path.c_str(), "READ");
+  if (file.IsZombie()) {
+    std::cerr << "Cannot open ROOT file " << path << "\n";
+    return 1;
+  }
+
+  auto* tree = (TTree*)file.Get("o2sim");
+  if (!tree) {
+    std::cerr << "Cannot find tree o2sim in file " << path << "\n";
+    return 1;
+  }
+
+  std::vector<o2::MCTrack>* tracks = nullptr;
+  tree->SetBranchAddress("MCTrack", &tracks);
+
+  std::vector<int> nSignal(injectedPDGs.size(), 0);
+  std::vector<int> nEventsWithSignal(injectedPDGs.size(), 0);
+
+  const auto nEvents = tree->GetEntries();
+
+  for (Long64_t i = 0; i < nEvents; ++i) {
+    tree->GetEntry(i);
+
+    std::vector<bool> found(injectedPDGs.size(), false);
+
+    for (const auto& track : *tracks) {
+      const int pdg = track.GetPdgCode();
+
+      auto it = std::find(
+          injectedPDGs.begin(),
+          injectedPDGs.end(),
+          pdg);
+
+      if (it != injectedPDGs.end()) {
+        const auto index = std::distance(injectedPDGs.begin(), it);
+        ++nSignal[index];
+        found[index] = true;
+      }
+    }
+
+    for (size_t j = 0; j < injectedPDGs.size(); ++j) {
+      if (found[j]) {
+        ++nEventsWithSignal[j];
+      }
+    }
+  }
+
+  std::cout << "--------------------------------\n";
+  std::cout << "# Events: " << nEvents << "\n";
+
+  for (size_t i = 0; i < injectedPDGs.size(); ++i) {
+    std::cout << "PDG: " << injectedPDGs[i] << "\n";
+    std::cout << "# Generated particles: " << nSignal[i] << "\n";
+    std::cout << "# Events containing particle: "
+              << nEventsWithSignal[i] << "\n";
+
+    if (nSignal[i] == 0) {
+      std::cerr << "No generated particles with PDG "
+                << injectedPDGs[i] << "\n";
+      return 1;
+    }
+  }
+
+  return 0;
+}

MC/config/PWGLF/pythia8/generator_pythia8_antideuteron_injection.C

@@ -0,0 +1,114 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+
+#include "FairGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/Pythia.h"
+#include "TRandom.h"
+#include "fairlogger/Logger.h"
+
+#include <cmath>
+#include <cstdint>
+#include <string>
+
+#endif // !defined(__CLING__) || defined(__ROOTCLING__)
+
+/// PYTHIA8 event generator with antideuteron injection.
+/// Events are generated until a final-state antiproton within the specified etaMax is found.
+/// The antiproton is then replaced by an antideuteron with twice its momentum.
+
+class GeneratorPythia8AntideuteronInjection
+    : public o2::eventgen::GeneratorPythia8 {
+
+ public:
+
+  GeneratorPythia8AntideuteronInjection(double etaMax = 0.8)
+      : o2::eventgen::GeneratorPythia8(),
+        mEtaMax(etaMax) {}
+
+  ~GeneratorPythia8AntideuteronInjection() override = default;
+
+  bool Init() override {
+    addSubGenerator(0, "Pythia8 events with antideuterons");
+    return o2::eventgen::GeneratorPythia8::Init();
+  }
+
+ protected:
+
+
+  bool generateEvent() override {
+    constexpr int maxAttempts = 10000;
+    for (int attempt = 0; attempt < maxAttempts; ++attempt) {
+      if (!GeneratorPythia8::generateEvent()) {
+        return false;
+      }
+
+      if (EventHasAntideuteron(mPythia.event)) {
+        LOG(debug) << ">> An antideuteron was formed!";
+
+        notifySubGenerator(0);
+        ++mGeneratedEvents;
+
+        return true;
+      }
+    }
+
+    LOG(error) << "Failed to generate an antideuteron after " << maxAttempts << " attempts";
+    return false;
+  }
+
+  bool EventHasAntideuteron(Pythia8::Event& event) {
+
+    const double md = 1.87561294257; // Deuteron mass in GeV/c^2
+    const int pdgAntid = -1000010020; // Antideuteron PDG code
+
+    for (int iPart = 0; iPart < event.size(); ++iPart) {
+
+      if (!event[iPart].isFinal()) {
+        continue;
+      }
+
+      if (std::abs(event[iPart].eta()) > mEtaMax) {
+        continue;
+      }
+
+      if (event[iPart].id() != -2212) {
+        continue;
+      }
+
+      double px = 2.0 * event[iPart].px();
+      double py = 2.0 * event[iPart].py();
+      double pz = 2.0 * event[iPart].pz();
+      double energy = std::sqrt(md * md + px * px + py * py + pz * pz);
+
+      // Remove antiproton from final state
+      event[iPart].statusNeg();
+
+      // Add antideuteron
+      const int iAntid = event.size();
+      event.append(pdgAntid, 121, iPart, 0, 0, 0, 0, 0, px, py, pz, energy, md);
+      // Update daughter indices
+      event[iPart].daughter1(iAntid);
+      event[iPart].daughter2(iAntid);
+      return true;
+    }
+    return false;
+  }
+
+ private:
+
+  double mEtaMax = 0.8;
+  uint64_t mGeneratedEvents = 0;
+};
+
+/// Generator factory
+
+FairGenerator* generateAntideuterons(double etaMax = 0.8) {
+
+  auto* myGenerator = new GeneratorPythia8AntideuteronInjection(etaMax);
+  auto seed = gRandom->TRandom::GetSeed() % 900000000;
+
+  myGenerator->readString("Random:setSeed = on");
+  myGenerator->readString("Random:seed = " + std::to_string(seed));
+
+  return myGenerator;
+}