From ba7248012a80ee81de560e6ee303e469e0ec314a Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Thu, 16 Jan 2025 18:12:15 +0100 Subject: [PATCH] Implemented basic pythia8 generator with deep trigger forwarding --- .../external/generator/pythia8_deep.C | 51 +++++++++++++ MC/config/examples/ini/pythia8_impactb.ini | 16 ++++ .../examples/ini/tests/pythia8_impactb.C | 75 +++++++++++++++++++ .../examples/pythia8/generator/pythia8_hi.cfg | 21 ++++++ .../examples/trigger/trigger_impactb.macro | 45 +++++++++++ 5 files changed, 208 insertions(+) create mode 100644 MC/config/examples/external/generator/pythia8_deep.C create mode 100644 MC/config/examples/ini/pythia8_impactb.ini create mode 100644 MC/config/examples/ini/tests/pythia8_impactb.C create mode 100644 MC/config/examples/pythia8/generator/pythia8_hi.cfg create mode 100644 MC/config/examples/trigger/trigger_impactb.macro diff --git a/MC/config/examples/external/generator/pythia8_deep.C b/MC/config/examples/external/generator/pythia8_deep.C new file mode 100644 index 000000000..72a7ed95b --- /dev/null +++ b/MC/config/examples/external/generator/pythia8_deep.C @@ -0,0 +1,51 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "Pythia8/Pythia.h" +#include "FairGenerator.h" +#include "FairPrimaryGenerator.h" +#include "Generators/Trigger.h" +#include +#include "Generators/GeneratorPythia8.h" +#include +#include "CommonUtils/ConfigurationMacroHelper.h" + +using namespace Pythia8; +#endif + +// Basic implementation of deep-triggered Pythia8 as external generator + +class GeneratorPythia8Deep : public o2::eventgen::GeneratorPythia8 +{ + // Settings are fed via the configuration file specified in the .ini file + // Triggers need to be handled like this otherwise the simulation with hybrid will + // not be able to recognise the provided triggers. + +public : + + GeneratorPythia8Deep() : o2::eventgen::GeneratorPythia8() + { + mInterface = reinterpret_cast(&mPythia); + mInterfaceName = "pythia8"; + o2::eventgen::DeepTrigger deeptrigger = nullptr; + + // external trigger via ini file + auto ¶ms = o2::eventgen::TriggerExternalParam::Instance(); + LOG(info) << "Setting up external trigger for Pythia8 with following parameters"; + LOG(info) << params; + auto external_trigger_filename = params.fileName; + auto external_trigger_func = params.funcName; + deeptrigger = o2::conf::GetFromMacro(external_trigger_filename, external_trigger_func, "o2::eventgen::DeepTrigger", "deeptrigger"); + if (!deeptrigger) + { + LOG(fatal) << "Failed to retrieve \'external trigger\': problem with configuration "; + } else { + LOG(info) << "External trigger for Pythia8 is set"; + addDeepTrigger(deeptrigger); + setTriggerMode(o2::eventgen::Generator::kTriggerOR); + } + } +}; + +FairGenerator *generator_pythia8_deep() +{ + return new GeneratorPythia8Deep(); +} \ No newline at end of file diff --git a/MC/config/examples/ini/pythia8_impactb.ini b/MC/config/examples/ini/pythia8_impactb.ini new file mode 100644 index 000000000..a37740e5f --- /dev/null +++ b/MC/config/examples/ini/pythia8_impactb.ini @@ -0,0 +1,16 @@ +### The setup uses an external event generator trigger which is +### defined in the following file and it is retrieved and configured +### according to the specified function call + +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/examples/external/generator/pythia8_deep.C +funcName=generator_pythia8_deep() + +[TriggerExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/examples/trigger/trigger_impactb.macro +funcName=trigger_impactb_pythia8(0.,15.) + +### This part configures Pythia8 + +[GeneratorPythia8] +config = ${O2DPG_MC_CONFIG_ROOT}/MC/config/examples/pythia8/generator/pythia8_hi.cfg diff --git a/MC/config/examples/ini/tests/pythia8_impactb.C b/MC/config/examples/ini/tests/pythia8_impactb.C new file mode 100644 index 000000000..5dbe1eccf --- /dev/null +++ b/MC/config/examples/ini/tests/pythia8_impactb.C @@ -0,0 +1,75 @@ +{ + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + 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); + + o2::dataformats::MCEventHeader *mcheader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &mcheader); + + // Check if all events are filled + auto nEvents = tree->GetEntries(); + for (Long64_t i = 0; i < nEvents; ++i) + { + tree->GetEntry(i); + if (tracks->empty()) + { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + } + // Check if there are 100 events, as simulated in the o2dpg-test + if (nEvents != 100) + { + std::cerr << "Expected 100 events, got " << nEvents << "\n"; + return 1; + } + // check if each event has two lead ions with the correct energy (based on Pythia simulation) + // exits if the particle is not lead 208 + bool isvalid; + for (int i = 0; i < nEvents; i++) + { + auto check = tree->GetEntry(i); + int count = 0; + for (int idxMCTrack = 0; idxMCTrack < tracks->size(); ++idxMCTrack) + { + auto track = tracks->at(idxMCTrack); + double energy = track.GetEnergy(); + // Check if track energy is right for the lead ions (a tolerance of 100 MeV is considered, straight equality does not work due to floating point precision) + if (std::abs(energy - 547158) < 1e-1) // Lead ion energy is 547158 MeV + { + if (track.GetPdgCode() != 1000822080) + { + std::cerr << "Found 547158 GeV particle with pdgID " << track.GetPdgCode() << "\n"; + return 1; + } + count++; + } + } + if (count < 2) + { + std::cerr << "Event " << i << " has less than 2 lead ions at 547158 GeV\n"; + return 1; + } + // Check if event impact parameter is < 15 fm + double impactParameter = mcheader->getInfo(o2::dataformats::MCInfoKeys::impactParameter, isvalid); + if (impactParameter > 15) + { + std::cerr << "Event " << i << " has impact parameter " << impactParameter << " fm outside range\n"; + return 1; + } + } + return 0; +} \ No newline at end of file diff --git a/MC/config/examples/pythia8/generator/pythia8_hi.cfg b/MC/config/examples/pythia8/generator/pythia8_hi.cfg new file mode 100644 index 000000000..d064b2dc4 --- /dev/null +++ b/MC/config/examples/pythia8/generator/pythia8_hi.cfg @@ -0,0 +1,21 @@ +### random +Random:setSeed = on +Random:seed = 0 + +### beams +Beams:idA = 1000822080 +Beams:idB = 1000822080 +Beams:eCM = 5300.000000 + +### heavy-ion settings (valid for Pb-Pb 5520 only) +HeavyIon:SigFitNGen = 0 +HeavyIon:SigFitDefPar = 13.88,1.84,0.22,0.0,0.0,0.0,0.0,0.0 +HeavyIon:bWidth = 14.48 + +### decays +ParticleDecays:limitTau0 = on +ParticleDecays:tau0Max = 10. + +### phase space cuts +PhaseSpace:pTHatMin = 0.000000 +PhaseSpace:pTHatMax = -1.000000 \ No newline at end of file diff --git a/MC/config/examples/trigger/trigger_impactb.macro b/MC/config/examples/trigger/trigger_impactb.macro new file mode 100644 index 000000000..076af058c --- /dev/null +++ b/MC/config/examples/trigger/trigger_impactb.macro @@ -0,0 +1,45 @@ +#include "Generators/Trigger.h" +#include "TParticle.h" +#include + +// a very simple trigger example, examining generated particles +o2::eventgen::Trigger trigger() +{ + // + return [](const std::vector& particles) -> bool { + std::cout << "Running trigger on event with size " << particles.size() << "\n"; + if (particles.size() > 10000) { + return true; + } + return false; + }; +} + +#include "Pythia8/Pythia.h" +#include "Pythia8/HIInfo.h" +#include +// a deep trigger example, looking into the internal generator state +o2::eventgen::DeepTrigger + trigger_impactb_pythia8(double bmin = 5., double bmax = 10.) +{ + return [bmin, bmax](void* interface, std::string name) -> bool { + if (!name.compare("pythia8")) { + auto py8 = reinterpret_cast(interface); +#if PYTHIA_VERSION_INTEGER < 8300 + auto hiinfo = py8->info.hiinfo; +#else + auto hiinfo = py8->info.hiInfo; +#endif + if (!hiinfo) { + LOG(fatal) << "Cannot define impact parameter: is \'pythia8\' running in heavy-ion mode?"; + } + auto b = hiinfo->b(); + auto selected = (b > bmin && b < bmax); + LOG(info) << "Impact parameter = " << b << " fm: " << (selected ? "selected" : "rejected"); + return selected; + } else { + LOG(fatal) << "Cannot define impact parameter for generator interface \'" << name << "\'"; + } + return false; + }; +}