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