diff --git a/MC/config/PWGLF/ini/GeneratorLF_antid_and_highpt.ini b/MC/config/PWGLF/ini/GeneratorLF_antid_and_highpt.ini new file mode 100644 index 000000000..b53e7d34b --- /dev/null +++ b/MC/config/PWGLF/ini/GeneratorLF_antid_and_highpt.ini @@ -0,0 +1,6 @@ +[GeneratorExternal] +fileName = ${O2DPG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_antid_and_highpt.C +funcName = generateAntidAndHighPt(0.3,5.0) + +[GeneratorPythia8] +config=${O2_ROOT}/share/Generators/egconfig/pythia8_inel.cfg diff --git a/MC/config/PWGLF/ini/tests/GeneratorLF_antid_and_highpt.C b/MC/config/PWGLF/ini/tests/GeneratorLF_antid_and_highpt.C new file mode 100644 index 000000000..ab733edb5 --- /dev/null +++ b/MC/config/PWGLF/ini/tests/GeneratorLF_antid_and_highpt.C @@ -0,0 +1,25 @@ +int External() { + std::string path{"o2sim_Kine.root"}; + + 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{}; + tree->SetBranchAddress("MCTrack", &tracks); + + auto nEvents = tree->GetEntries(); + auto nSelected = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 1000010020"); + if (nSelected == 0) { + std::cerr << "No event of interest\n"; + return 1; + } + return 0; +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_antid_and_highpt.C b/MC/config/PWGLF/pythia8/generator_pythia8_antid_and_highpt.C new file mode 100644 index 000000000..327b5aa84 --- /dev/null +++ b/MC/config/PWGLF/pythia8/generator_pythia8_antid_and_highpt.C @@ -0,0 +1,167 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "FairGenerator.h" +#include "FairPrimaryGenerator.h" +#include "Generators/GeneratorPythia8.h" +#include "Pythia8/Pythia.h" +#include "TDatabasePDG.h" +#include "TMath.h" +#include "TParticlePDG.h" +#include "TRandom3.h" +#include "TSystem.h" +#include "fairlogger/Logger.h" +#include +#include +#include +#include +using namespace Pythia8; +#endif + +/// Pythia8 event generator for pp collisions +/// Selection of events with leading particle (pt>ptThreshold) and containing at +/// least one (anti)deuteron produced by simple coalescence + +class GeneratorPythia8AntidAndHighPt : public o2::eventgen::GeneratorPythia8 { +public: + /// Constructor + GeneratorPythia8AntidAndHighPt(double p0 = 0.3, double pt_leading = 5.0) + : o2::eventgen::GeneratorPythia8() { + mP0 = p0; + mPt_leading = pt_leading; + } + /// Destructor + ~GeneratorPythia8AntidAndHighPt() = default; + + bool Init() override { + addSubGenerator(0, "Pythia8 with (anti)deuterons and high pt particle"); + return o2::eventgen::GeneratorPythia8::Init(); + } + +protected: + bool generateEvent() override { + fmt::printf(">> Generating event %d\n", mGeneratedEvents); + + bool genOk = false; + int localCounter{0}; + while (!genOk) { + if (GeneratorPythia8::generateEvent()) { + genOk = selectEvent(mPythia.event); + } + localCounter++; + } + fmt::printf( + ">> Generation of event of interest successful after %i iterations\n", + localCounter); + std::cout << std::endl << std::endl; + notifySubGenerator(0); + + mGeneratedEvents++; + + return true; + } + + bool selectEvent(Pythia8::Event &event) { + + bool has_particle_of_interest = false; + + // Deuteron Mass [GeV] + double md = 1.87561294257; + + // Protons and Neutrons + vector proton_ID; + vector neutron_ID; + vector proton_status; + vector neutron_status; + + // ptMax + double pt_max{0}; + + for (auto iPart{0}; iPart < event.size(); ++iPart) { + + // Only final-state particles + if (event[iPart].status() <= 0) { + continue; + } + + //(Anti)Proton selection + if (abs(event[iPart].id()) == 2212) { + proton_ID.push_back(iPart); + proton_status.push_back(0); + } + + //(Anti)Neutron selection + if (abs(event[iPart].id()) == 2112) { + neutron_ID.push_back(iPart); + neutron_status.push_back(0); + } + + if (std::abs(event[iPart].eta()) < 0.8 && (!event[iPart].isNeutral()) && + event[iPart].pT() > pt_max) + pt_max = event[iPart].pT(); + } + + // Skip Events with no leading particle + if (pt_max < mPt_leading) + return false; + + if (proton_ID.size() < 1 || neutron_ID.size() < 1) + return false; + + for (uint32_t ip = 0; ip < proton_ID.size(); ++ip) { + + // Skip used protons + if (proton_status[ip] < 0) { + continue; + } + for (uint32_t in = 0; in < neutron_ID.size(); ++in) { + + // Skip used neutrons + if (neutron_status[in] < 0) { + continue; + } + + auto sign_p = + event[proton_ID[ip]].id() / abs(event[proton_ID[ip]].id()); + auto sign_n = + event[neutron_ID[in]].id() / abs(event[neutron_ID[in]].id()); + + auto p1 = event[proton_ID[ip]].p(); + auto p2 = event[neutron_ID[in]].p(); + auto p = p1 + p2; + p1.bstback(p); + p2.bstback(p); + + // Coalescence + if (p1.pAbs() <= mP0 && p2.pAbs() <= mP0 && sign_p == sign_n) { + p.e(std::hypot(p.pAbs(), md)); + event.append(sign_p * 1000010020, 121, 0, 0, 0, 0, 0, 0, p.px(), + p.py(), p.pz(), p.e(), md); + event[proton_ID[ip]].statusNeg(); + event[proton_ID[ip]].daughter1(event.size() - 1); + event[neutron_ID[in]].statusNeg(); + event[neutron_ID[in]].daughter1(event.size() - 1); + proton_status[ip] = -1; + neutron_status[in] = -1; + has_particle_of_interest = true; + } + } + } + + return has_particle_of_interest; + } + +private: + double mP0 = 0.3; + double mPt_leading = 5.0; + uint64_t mGeneratedEvents = 0; +}; + +///___________________________________________________________ +FairGenerator *generateAntidAndHighPt(double p0 = 0.3, + double pt_leading = 5.0) { + + auto myGenerator = new GeneratorPythia8AntidAndHighPt(p0, pt_leading); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGenerator->readString("Random:setSeed on"); + myGenerator->readString("Random:seed " + std::to_string(seed)); + return myGenerator; +} diff --git a/MC/run/PWGLF/run_GeneratorLF_antid_and_highpt.sh b/MC/run/PWGLF/run_GeneratorLF_antid_and_highpt.sh new file mode 100644 index 000000000..ce8c90be3 --- /dev/null +++ b/MC/run/PWGLF/run_GeneratorLF_antid_and_highpt.sh @@ -0,0 +1,28 @@ +#!/bin/bash + + +# make sure O2DPG + O2 is loaded +[ ! "${O2DPG_ROOT}" ] && echo "Error: This needs O2DPG loaded" && exit 1 +[ ! "${O2_ROOT}" ] && echo "Error: This needs O2 loaded" && exit 1 + +# ----------- LOAD UTILITY FUNCTIONS -------------------------- +. ${O2_ROOT}/share/scripts/jobutils.sh + +# ----------- START ACTUAL JOB ----------------------------- + +NWORKERS=${NWORKERS:-8} +MODULES="--skipModules ZDC" +SIMENGINE=${SIMENGINE:-TGeant4} +NSIGEVENTS=${NSIGEVENTS:-1} +NTIMEFRAMES=${NTIMEFRAMES:-1} +INTRATE=${INTRATE:-500000} +SYSTEM=${SYSTEM:-pp} +ENERGY=${ENERGY:-13600} +[[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED="" + +# create workflow +${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py -eCM ${ENERGY} -col ${SYSTEM} -gen external -j ${NWORKERS} -ns ${NSIGEVENTS} -tf ${NTIMEFRAMES} -interactionRate ${INTRATE} -confKey "Diamond.width[0]=0.005;Diamond.width[1]=0.005;Diamond.width[2]=6." -e ${SIMENGINE} ${SEED} -mod "--skipModules ZDC" \ + -ini $O2DPG_ROOT/MC/config/PWGLF/ini/GeneratorLF_antid_and_highpt.ini + +# run workflow +${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py -f workflow.json -tt aod --cpu-limit 32