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 "ActsPlugins/FpeMonitoring/FpeMonitor.hpp"
0012
0013 #include <ostream>
0014 #include <random>
0015
0016 #include <HepMC3/WriterAscii.h>
0017 #include <Pythia8/Pythia.h>
0018 #include <Pythia8Plugins/HepMC3.h>
0019
0020 namespace ActsExamples {
0021
0022 struct Pythia8RandomEngineWrapper : public Pythia8::RndmEngine {
0023   RandomEngine* rng{nullptr};
0024
0025   struct {
0026     std::size_t numUniformRandomNumbers = 0;
0027     double first = std::numeric_limits<double>::quiet_NaN();
0028     double last = std::numeric_limits<double>::quiet_NaN();
0029   } statistics;
0030
0031   Pythia8RandomEngineWrapper() = default;
0032
0033   double flat() override {
0034     if (rng == nullptr) {
0035       throw std::runtime_error(
0036           "Pythia8RandomEngineWrapper: no random engine set");
0037     }
0038
0039     double value = std::uniform_real_distribution<double>(0.0, 1.0)(*rng);
0040     if (statistics.numUniformRandomNumbers == 0) {
0041       statistics.first = value;
0042     }
0043     statistics.last = value;
0044     statistics.numUniformRandomNumbers++;
0045     return value;
0046   }
0047
0048   void setRandomEngine(RandomEngine& rng_) { rng = &rng_; }
0049   void clearRandomEngine() { rng = nullptr; }
0050 };
0051
0052 struct Pythia8GeneratorImpl {
0053   std::unique_ptr<HepMC3::Writer> m_hepMC3Writer;
0054   std::unique_ptr<HepMC3::Pythia8ToHepMC3> m_hepMC3Converter;
0055   std::shared_ptr<Pythia8RandomEngineWrapper> m_pythia8RndmEngine;
0056 };
0057
0058 Pythia8Generator::Pythia8Generator(const Config& cfg, Acts::Logging::Level lvl)
0059     : m_cfg(cfg),
0060       m_logger(Acts::getDefaultLogger("Pythia8Generator", lvl)),
0061       m_pythia8(std::make_unique<Pythia8::Pythia>("", false)) {
0062   ACTS_INFO("Pythia8Generator: init");
0063   // disable all output by default but allow re-enable via config
0064   m_pythia8->settings.flag("Print:quiet", true);
0065   for (const auto& setting : m_cfg.settings) {
0066     ACTS_VERBOSE("use Pythia8 setting '" << setting << "'");
0067     m_pythia8->readString(setting.c_str());
0068   }
0069   m_pythia8->settings.mode("Beams:idA", m_cfg.pdgBeam0);
0070   m_pythia8->settings.mode("Beams:idB", m_cfg.pdgBeam1);
0071   m_pythia8->settings.mode("Beams:frameType", 1);
0072   m_pythia8->settings.parm("Beams:eCM",
0073                            m_cfg.cmsEnergy / Acts::UnitConstants::GeV);
0074
0075   m_impl = std::make_unique<Pythia8GeneratorImpl>();
0076
0077   m_impl->m_pythia8RndmEngine = std::make_shared<Pythia8RandomEngineWrapper>();
0078
0079 #if PYTHIA_VERSION_INTEGER >= 8310
0080   m_pythia8->setRndmEnginePtr(m_impl->m_pythia8RndmEngine);
0081 #else
0082   m_pythia8->setRndmEnginePtr(m_impl->m_pythia8RndmEngine.get());
0083 #endif
0084
0085   RandomEngine rng{m_cfg.initializationSeed};
0086   m_impl->m_pythia8RndmEngine->setRandomEngine(rng);
0087   m_pythia8->init();
0088   m_impl->m_pythia8RndmEngine->clearRandomEngine();
0089
0090   m_impl->m_hepMC3Converter = std::make_unique<HepMC3::Pythia8ToHepMC3>();
0091   if (m_cfg.writeHepMC3.has_value()) {
0092     ACTS_DEBUG("Initializing HepMC3 output to: " << m_cfg.writeHepMC3.value());
0093     m_impl->m_hepMC3Writer =
0094         std::make_unique<HepMC3::WriterAscii>(m_cfg.writeHepMC3.value());
0095   }
0096 }
0097
0098 // needed to allow unique_ptr of forward-declared Pythia class
0099 Pythia8Generator::~Pythia8Generator() {
0100   if (m_impl->m_hepMC3Writer) {
0101     m_impl->m_hepMC3Writer->close();
0102   }
0103
0104   ACTS_DEBUG("Pythia8Generator produced "
0105              << m_impl->m_pythia8RndmEngine->statistics.numUniformRandomNumbers
0106              << " uniform random numbers");
0107   ACTS_DEBUG("                 first = "
0108              << m_impl->m_pythia8RndmEngine->statistics.first);
0109   ACTS_DEBUG("                  last = "
0110              << m_impl->m_pythia8RndmEngine->statistics.last);
0111 }
0112
0113 std::shared_ptr<HepMC3::GenEvent> Pythia8Generator::operator()(
0114     RandomEngine& rng) {
0115   using namespace Acts::UnitLiterals;
0116
0117   // pythia8 is not thread safe and generation needs to be protected
0118   std::lock_guard<std::mutex> lock(m_pythia8Mutex);
0119   // use per-thread random engine also in pythia
0120
0121   m_impl->m_pythia8RndmEngine->setRandomEngine(rng);
0122
0123   {
0124     ActsPlugins::FpeMonitor mon{0};  // disable all FPEs while we're in Pythia8
0125     m_pythia8->next();
0126   }
0127
0128   if (m_cfg.printShortEventListing) {
0129     m_pythia8->process.list();
0130   }
0131   if (m_cfg.printLongEventListing) {
0132     m_pythia8->event.list();
0133   }
0134
0135   auto genEvent = std::make_shared<HepMC3::GenEvent>();
0136   genEvent->set_units(HepMC3::Units::GEV, HepMC3::Units::MM);
0137
0138   assert(m_impl->m_hepMC3Converter != nullptr);
0139   m_impl->m_hepMC3Converter->fill_next_event(*m_pythia8, genEvent.get(),
0140                                              genEvent->event_number());
0141
0142   if (m_impl->m_hepMC3Converter && m_impl->m_hepMC3Writer) {
0143     m_impl->m_hepMC3Writer->write_event(*genEvent);
0144   }
0145
0146   m_impl->m_pythia8RndmEngine->clearRandomEngine();
0147
0148   return genEvent;
0149 }
0150
0151 }  // namespace ActsExamples