diff --git a/MC/config/PWGLF/ini/GeneratorDoubleLambdaTriggered.ini b/MC/config/PWGLF/ini/GeneratorDoubleLambdaTriggered.ini new file mode 100644 index 000000000..f1f484e81 --- /dev/null +++ b/MC/config/PWGLF/ini/GeneratorDoubleLambdaTriggered.ini @@ -0,0 +1,9 @@ +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_doubleLambdas.C +funcName=generateDoubleLambda(4, 0.2, 10, 0.8) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_ropes_136tev.cfg + +[DecayerPythia8] +config[0]=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/decayer/base.cfg diff --git a/MC/config/PWGLF/ini/tests/GeneratorDoubleLambdaTriggered.C b/MC/config/PWGLF/ini/tests/GeneratorDoubleLambdaTriggered.C new file mode 100644 index 000000000..96fcd8bbd --- /dev/null +++ b/MC/config/PWGLF/ini/tests/GeneratorDoubleLambdaTriggered.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 nLambda = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 3122"); + + if (nLambda == 0) { + std::cerr << "No event of interest\n"; + return 1; + } + return 0; +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_doubleLambdas.C b/MC/config/PWGLF/pythia8/generator_pythia8_doubleLambdas.C new file mode 100644 index 000000000..b03a57ef0 --- /dev/null +++ b/MC/config/PWGLF/pythia8/generator_pythia8_doubleLambdas.C @@ -0,0 +1,162 @@ +#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 "TVector2.h" +#include "fairlogger/Logger.h" +#include +#include +#include +#include +using namespace Pythia8; +#endif + +/// Event generator using Pythia ropes +/// One event of interest is triggered after n minimum-bias events. +/// Event of interest : an event that contains at least two generated (anti-)Lambdas + +class GeneratorPythia8DoubleLambda : public o2::eventgen::GeneratorPythia8 +{ +public: + /// Constructor + GeneratorPythia8DoubleLambda(int gapSize = 4, double minPt = 0.2, double maxPt = 10, double maxEta = 0.8) + : o2::eventgen::GeneratorPythia8(), + mGapSize(gapSize), + mMinPt(minPt), + mMaxPt(maxPt), + mMaxEta(maxEta) + { + fmt::printf( + ">> Pythia8 generator: two (anti-)Lambdas, gap = %d, minPtLambda = %f, maxPtLambda = %f, |etaLambda| < %f\n", gapSize, minPt, maxPt, maxEta); + } + /// Destructor + ~GeneratorPythia8DoubleLambda() = default; + + bool Init() override + { + addSubGenerator(0, "Pythia8 events with two (anti-)Lambdas"); + return o2::eventgen::GeneratorPythia8::Init(); + } + +protected: + /// Check if particle is physical primary or from HF decay + bool isLambdaPhysicalPrimaryOrFromHF(const Pythia8::Particle &p, const Pythia8::Event &event) + { + // Select only final-state particles + if (!p.isFinal()) + { + return false; + } + + // Walk up ancestry + int motherId = p.mother1(); + + while (motherId > 0) + { + + // Get mother + const auto &mother = event[motherId]; + const int absMotherPdg = std::abs(mother.id()); + + // Check if particle is from HF decay + if ((absMotherPdg / 100 == 4) || (absMotherPdg / 100 == 5) || (absMotherPdg / 1000 == 4) || (absMotherPdg / 1000 == 5)) + { + return true; + } + + // Reject non-physical primary hadrons + if (mother.isHadron() && mother.tau0() > 1.0) + { + return false; + } + motherId = mother.mother1(); + } + return true; + } + + bool generateEvent() override + { + 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-event (no strangeness check)\n"); + } + else + { + while (!genOk && localCounter < kMaxTries) + { + if (GeneratorPythia8::generateEvent()) + { + genOk = selectEvent(mPythia.event); + } + localCounter++; + } + if (!genOk) + { + fmt::printf("Failed to generate triggered event after %d tries\n", kMaxTries); + return false; + } + fmt::printf(">> Triggered event: event accepted after %d iterations (double (anti-)Lambdas)\n", localCounter); + } + + notifySubGenerator(0); + mGeneratedEvents++; + return true; + } + + bool selectEvent(Pythia8::Event &event) + { + + std::vector particleID; + int nLambda{0}; + + for (int i = 0; i < event.size(); i++) + { + const auto &p = event[i]; + if (std::abs(p.eta()) > mMaxEta || p.pT() < mMinPt || p.pT() > mMaxPt) + continue; + if (!isLambdaPhysicalPrimaryOrFromHF(p, event)) + continue; + + if (std::abs(p.id()) == 3122) + nLambda++; + + particleID.emplace_back(i); + } + if (nLambda < 2) + return false; + + return true; + } + +private: + int mGapSize{4}; + double mMinPt{0.2}; + double mMaxPt{10.0}; + double mMaxEta{0.8}; + uint64_t mGeneratedEvents{0}; +}; + +///___________________________________________________________ +FairGenerator *generateDoubleLambda(int gap = 4, double minPt = 0.2, double maxPt = 10, double maxEta = 0.8) +{ + + auto myGenerator = new GeneratorPythia8DoubleLambda(gap, minPt, maxPt, maxEta); + 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_DoubleLambdaTriggered.sh b/MC/run/PWGLF/run_DoubleLambdaTriggered.sh new file mode 100755 index 000000000..061ae86e5 --- /dev/null +++ b/MC/run/PWGLF/run_DoubleLambdaTriggered.sh @@ -0,0 +1,39 @@ +#!/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 + +# ----------- CONFIGURE -------------------------- +export IGNORE_VALIDITYCHECK_OF_CCDB_LOCALCACHE=1 +#export ALICEO2_CCDB_LOCALCACHE=.ccdb + + +# ----------- START ACTUAL JOB ----------------------------- + +NWORKERS=${NWORKERS:-8} +SIMENGINE=${SIMENGINE:-TGeant4} +NSIGEVENTS=${NSIGEVENTS:-100} +NTIMEFRAMES=${NTIMEFRAMES:-1} +INTRATE=${INTRATE:-50000} +SYSTEM=${SYSTEM:-pp} +ENERGY=${ENERGY:-13600} +CFGINIFILE=${CFGINIFILE:-"${O2DPG_ROOT}/MC/config/PWGLF/ini/GeneratorDoubleLambdaTriggered.ini"} +SEED="-seed 1995" +#[[ ${SPLITID} != "" ]] && SEED="-seed ${SPLITID}" || SEED="" + +echo "NWORKERS = $NWORKERS" + +# create workflow +O2_SIM_WORKFLOW=${O2_SIM_WORKFLOW:-"${O2DPG_ROOT}/MC/bin/o2dpg_sim_workflow.py"} +$O2_SIM_WORKFLOW -eCM ${ENERGY} -col ${SYSTEM} -gen external \ + -j ${NWORKERS} \ + -ns ${NSIGEVENTS} -tf ${NTIMEFRAMES} -interactionRate ${INTRATE} \ + -confKey "Diamond.width[2]=6." \ + ${SEED} \ + -e ${SIMENGINE} \ + -ini $CFGINIFILE + +# run workflow +O2_SIM_WORKFLOW_RUNNER=${O2_SIM_WORKFLOW_RUNNER:-"${O2DPG_ROOT}/MC/bin/o2_dpg_workflow_runner.py"} +$O2_SIM_WORKFLOW_RUNNER -f workflow.json -tt aod --cpu-limit $NWORKERS