ActsExamples/Generators/Pythia8ProcessGenerator.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Generators/Pythia8ProcessGenerator.hpp"
0010
0011 #include "Acts/Utilities/MathHelpers.hpp"
0012 #include "ActsExamples/EventData/SimVertex.hpp"
0013 #include "ActsFatras/EventData/Barcode.hpp"
0014
0015 #include <algorithm>
0016 #include <iterator>
0017 #include <ostream>
0018 #include <random>
0019 #include <utility>
0020
0021 #include <Pythia8/Pythia.h>
0022
0023 namespace ActsExamples {
0024
0025 struct Pythia8RandomEngineWrapper : public Pythia8::RndmEngine {
0026   RandomEngine* rng{nullptr};
0027
0028   struct {
0029     std::size_t numUniformRandomNumbers = 0;
0030     double first = std::numeric_limits<double>::quiet_NaN();
0031     double last = std::numeric_limits<double>::quiet_NaN();
0032   } statistics;
0033
0034   Pythia8RandomEngineWrapper() = default;
0035
0036   double flat() override {
0037     if (rng == nullptr) {
0038       throw std::runtime_error(
0039           "Pythia8RandomEngineWrapper: no random engine set");
0040     }
0041
0042     double value = std::uniform_real_distribution<double>(0.0, 1.0)(*rng);
0043     if (statistics.numUniformRandomNumbers == 0) {
0044       statistics.first = value;
0045     }
0046     statistics.last = value;
0047     statistics.numUniformRandomNumbers++;
0048     return value;
0049   }
0050
0051   void setRandomEngine(RandomEngine& rng_) { rng = &rng_; }
0052   void clearRandomEngine() { rng = nullptr; }
0053 };
0054
0055 Pythia8Generator::Pythia8Generator(const Config& cfg, Acts::Logging::Level lvl)
0056     : m_cfg(cfg),
0057       m_logger(Acts::getDefaultLogger("Pythia8Generator", lvl)),
0058       m_pythia8(std::make_unique<Pythia8::Pythia>("", false)) {
0059   // disable all output by default but allow re-enable via config
0060   m_pythia8->settings.flag("Print:quiet", true);
0061   for (const auto& setting : m_cfg.settings) {
0062     ACTS_VERBOSE("use Pythia8 setting '" << setting << "'");
0063     m_pythia8->readString(setting.c_str());
0064   }
0065   m_pythia8->settings.mode("Beams:idA", m_cfg.pdgBeam0);
0066   m_pythia8->settings.mode("Beams:idB", m_cfg.pdgBeam1);
0067   m_pythia8->settings.mode("Beams:frameType", 1);
0068   m_pythia8->settings.parm("Beams:eCM",
0069                            m_cfg.cmsEnergy / Acts::UnitConstants::GeV);
0070
0071   m_pythia8RndmEngine = std::make_shared<Pythia8RandomEngineWrapper>();
0072
0073 #if PYTHIA_VERSION_INTEGER >= 8310
0074   m_pythia8->setRndmEnginePtr(m_pythia8RndmEngine);
0075 #else
0076   m_pythia8->setRndmEnginePtr(m_pythia8RndmEngine.get());
0077 #endif
0078
0079   RandomEngine rng{m_cfg.initializationSeed};
0080   m_pythia8RndmEngine->setRandomEngine(rng);
0081   m_pythia8->init();
0082   m_pythia8RndmEngine->clearRandomEngine();
0083 }
0084
0085 // needed to allow unique_ptr of forward-declared Pythia class
0086 Pythia8Generator::~Pythia8Generator() {
0087   ACTS_INFO("Pythia8Generator produced "
0088             << m_pythia8RndmEngine->statistics.numUniformRandomNumbers
0089             << " uniform random numbers");
0090   ACTS_INFO(
0091       "                 first = " << m_pythia8RndmEngine->statistics.first);
0092   ACTS_INFO(
0093       "                  last = " << m_pythia8RndmEngine->statistics.last);
0094 }
0095
0096 std::pair<SimVertexContainer, SimParticleContainer>
0097 Pythia8Generator::operator()(RandomEngine& rng) {
0098   using namespace Acts::UnitLiterals;
0099
0100   SimVertexContainer::sequence_type vertices;
0101   SimParticleContainer::sequence_type particles;
0102
0103   // pythia8 is not thread safe and generation needs to be protected
0104   std::lock_guard<std::mutex> lock(m_pythia8Mutex);
0105   // use per-thread random engine also in pythia
0106
0107   m_pythia8RndmEngine->setRandomEngine(rng);
0108
0109   {
0110     Acts::FpeMonitor mon{0};  // disable all FPEs while we're in Pythia8
0111     m_pythia8->next();
0112   }
0113
0114   if (m_cfg.printShortEventListing) {
0115     m_pythia8->process.list();
0116   }
0117   if (m_cfg.printLongEventListing) {
0118     m_pythia8->event.list();
0119   }
0120
0121   // create the primary vertex
0122   vertices.emplace_back(SimVertexBarcode{0}, Acts::Vector4(0., 0., 0., 0.));
0123
0124   // convert generated final state particles into internal format
0125   for (int ip = 0; ip < m_pythia8->event.size(); ++ip) {
0126     const auto& genParticle = m_pythia8->event[ip];
0127
0128     // ignore beam particles
0129     if (genParticle.statusHepMC() == 4) {
0130       continue;
0131     }
0132     // only interested in final, visible particles
0133     if (!genParticle.isFinal()) {
0134       continue;
0135     }
0136     if (!genParticle.isVisible()) {
0137       continue;
0138     }
0139
0140     // production vertex. Pythia8 time uses units mm/c, and we use c=1
0141     Acts::Vector4 pos4(genParticle.xProd() * 1_mm, genParticle.yProd() * 1_mm,
0142                        genParticle.zProd() * 1_mm, genParticle.tProd() * 1_mm);
0143
0144     // define the particle identifier including possible secondary vertices
0145
0146     SimBarcode particleId(0u);
0147     // ensure particle identifier component is non-zero
0148     particleId.setParticle(1u + particles.size());
0149     // only secondaries have a defined vertex position
0150     if (m_cfg.labelSecondaries && genParticle.hasVertex()) {
0151       // either add to existing secondary vertex if exists or create new one
0152
0153       // check if an existing vertex is close enough
0154       auto it =
0155           std::ranges::find_if(vertices, [&pos4, this](const SimVertex& v) {
0156             return (pos4.head<3>() - v.position()).norm() <
0157                    m_cfg.spatialVertexThreshold;
0158           });
0159
0160       if (it != vertices.end()) {
0161         particleId.setVertexSecondary(std::distance(vertices.begin(), it));
0162         it->outgoing.insert(particleId);
0163       } else {
0164         // no matching secondary vertex exists -> create new one
0165         particleId.setVertexSecondary(vertices.size());
0166         auto& vertex = vertices.emplace_back(
0167             static_cast<SimVertexBarcode>(particleId.vertexId()), pos4);
0168         vertex.outgoing.insert(particleId);
0169         ACTS_VERBOSE("created new secondary vertex " << pos4.transpose());
0170       }
0171     } else {
0172       auto& primaryVertex = vertices.front();
0173       primaryVertex.outgoing.insert(particleId);
0174     }
0175
0176     // construct internal particle
0177     const auto pdg = static_cast<Acts::PdgParticle>(genParticle.id());
0178     const auto charge = genParticle.charge() * 1_e;
0179     const auto mass = genParticle.m0() * 1_GeV;
0180     SimParticleState particle(particleId, pdg, charge, mass);
0181     particle.setPosition4(pos4);
0182     // normalization/ units are not import for the direction
0183     particle.setDirection(genParticle.px(), genParticle.py(), genParticle.pz());
0184     particle.setAbsoluteMomentum(
0185         Acts::fastHypot(genParticle.px(), genParticle.py(), genParticle.pz()) *
0186         1_GeV);
0187
0188     particles.push_back(SimParticle(particle, particle));
0189   }
0190
0191   std::pair<SimVertexContainer, SimParticleContainer> out;
0192   out.first.insert(vertices.begin(), vertices.end());
0193   out.second.insert(particles.begin(), particles.end());
0194
0195   m_pythia8RndmEngine->clearRandomEngine();
0196
0197   return out;
0198 }
0199
0200 }  // namespace ActsExamples