diff --git a/MC/config/PWGLF/hybrid/hybrid_coalescence.json b/MC/config/PWGLF/hybrid/hybrid_coalescence.json new file mode 100644 index 000000000..32c3f54ae --- /dev/null +++ b/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] +} diff --git a/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini b/MC/config/PWGLF/ini/GeneratorAntideuteronCoalescenceJet.ini new file mode 100644 index 000000000..1a6394c65 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C b/MC/config/PWGLF/ini/tests/GeneratorAntideuteronCoalescenceJet.C new file mode 100644 index 000000000..f456070a6 --- /dev/null +++ b/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 +#include +#include +#include +#endif // !defined(__CLING__) || defined(__ROOTCLING__) + +int External() +{ + std::string path{"o2sim_Kine.root"}; + + std::vector 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* tracks = nullptr; + tree->SetBranchAddress("MCTrack", &tracks); + + std::vector nSignal(injectedPDGs.size(), 0); + std::vector nEventsWithSignal(injectedPDGs.size(), 0); + + const auto nEvents = tree->GetEntries(); + + for (Long64_t i = 0; i < nEvents; ++i) { + tree->GetEntry(i); + + std::vector 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; +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_antideuteron_injection.C b/MC/config/PWGLF/pythia8/generator_pythia8_antideuteron_injection.C new file mode 100644 index 000000000..c4f17ea52 --- /dev/null +++ b/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 +#include +#include + +#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; +}