ActsExamples/Generators/Pythia8ProcessGenerator.cpp

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