src/Utilities/ParametricParticleGenerator.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/Utilities/ParametricParticleGenerator.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/ParticleData.hpp"
0013 #include "Acts/Utilities/AngleHelpers.hpp"
0014
0015 #include <limits>
0016 #include <memory>
0017 #include <utility>
0018
0019 #include <HepMC3/Attribute.h>
0020 #include <HepMC3/FourVector.h>
0021 #include <HepMC3/GenEvent.h>
0022 #include <HepMC3/GenParticle.h>
0023 #include <HepMC3/GenVertex.h>
0024
0025 using namespace Acts::UnitLiterals;
0026
0027 namespace ActsExamples {
0028
0029 ParametricParticleGenerator::ParametricParticleGenerator(const Config& cfg)
0030     : m_cfg(cfg),
0031       m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg).value_or(0))) {
0032   m_pdgChoices = {
0033       m_cfg.pdg,
0034       static_cast<Acts::PdgParticle>(-m_cfg.pdg),
0035   };
0036
0037   // choose between particle/anti-particle if requested
0038   // the upper limit of the distribution is inclusive
0039   m_particleTypeChoice = UniformIndex(0u, m_cfg.randomizeCharge ? 1u : 0u);
0040   m_phiDist = UniformReal(m_cfg.phiMin, m_cfg.phiMax);
0041
0042   if (m_cfg.etaUniform) {
0043     double etaMin = Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMin);
0044     double etaMax = Acts::AngleHelpers::etaFromTheta(m_cfg.thetaMax);
0045
0046     UniformReal etaDist(etaMin, etaMax);
0047
0048     m_sinCosThetaDist =
0049         [=](RandomEngine& rng) mutable -> std::pair<double, double> {
0050       const double eta = etaDist(rng);
0051       const double theta = Acts::AngleHelpers::thetaFromEta(eta);
0052       return {std::sin(theta), std::cos(theta)};
0053     };
0054   } else {
0055     // since we want to draw the direction uniform on the unit sphere, we must
0056     // draw from cos(theta) instead of theta. see e.g.
0057     // https://mathworld.wolfram.com/SpherePointPicking.html
0058     double cosThetaMin = std::cos(m_cfg.thetaMin);
0059     // ensure upper bound is included. see e.g.
0060     // https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution
0061     double cosThetaMax = std::nextafter(std::cos(m_cfg.thetaMax),
0062                                         std::numeric_limits<double>::max());
0063
0064     UniformReal cosThetaDist(cosThetaMin, cosThetaMax);
0065
0066     m_sinCosThetaDist =
0067         [=](RandomEngine& rng) mutable -> std::pair<double, double> {
0068       const double cosTheta = cosThetaDist(rng);
0069       return {std::sqrt(1 - cosTheta * cosTheta), cosTheta};
0070     };
0071   }
0072
0073   if (m_cfg.pLogUniform) {
0074     // distributes p or pt uniformly in log space
0075     UniformReal dist(std::log(m_cfg.pMin), std::log(m_cfg.pMax));
0076
0077     m_somePDist = [=](RandomEngine& rng) mutable {
0078       return std::exp(dist(rng));
0079     };
0080   } else {
0081     // distributes p or pt uniformly
0082     UniformReal dist(m_cfg.pMin, m_cfg.pMax);
0083
0084     m_somePDist = [=](RandomEngine& rng) mutable { return dist(rng); };
0085   }
0086 }
0087
0088 std::shared_ptr<HepMC3::GenEvent> ParametricParticleGenerator::operator()(
0089     RandomEngine& rng) {
0090   auto event = std::make_shared<HepMC3::GenEvent>();
0091
0092   auto primaryVertex = std::make_shared<HepMC3::GenVertex>();
0093   primaryVertex->set_position(HepMC3::FourVector(0., 0., 0., 0.));
0094   event->add_vertex(primaryVertex);
0095
0096   primaryVertex->add_attribute("acts",
0097                                std::make_shared<HepMC3::BoolAttribute>(true));
0098
0099   // Produce pseudo beam particles
0100   auto beamParticle = std::make_shared<HepMC3::GenParticle>();
0101   beamParticle->set_momentum(HepMC3::FourVector(0., 0., 0., 0.));
0102   beamParticle->set_generated_mass(0.);
0103   beamParticle->set_pid(Acts::PdgParticle::eInvalid);
0104   beamParticle->set_status(4);
0105   primaryVertex->add_particle_in(beamParticle);
0106
0107   // counter will be reused as barcode particle number which must be non-zero.
0108   for (std::size_t ip = 1; ip <= m_cfg.numParticles; ++ip) {
0109     // draw parameters
0110     const unsigned int type = m_particleTypeChoice(rng);
0111     const Acts::PdgParticle pdg = m_pdgChoices[type];
0112     const double phi = m_phiDist(rng);
0113     const double someP = m_somePDist(rng);
0114
0115     const auto [sinTheta, cosTheta] = m_sinCosThetaDist(rng);
0116     // we already have sin/cos theta. they can be used directly
0117     const Acts::Vector3 dir = {sinTheta * std::cos(phi),
0118                                sinTheta * std::sin(phi), cosTheta};
0119
0120     const double p = someP * (m_cfg.pTransverse ? 1. / sinTheta : 1.);
0121
0122     // construct the particle;
0123     Acts::Vector3 momentum = p * dir;
0124     auto particle = std::make_shared<HepMC3::GenParticle>();
0125     HepMC3::FourVector hepMcMomentum(momentum.x() / 1_GeV, momentum.y() / 1_GeV,
0126                                      momentum.z() / 1_GeV,
0127                                      std::hypot(p, m_mass) / 1_GeV);
0128     particle->set_momentum(hepMcMomentum);
0129     particle->set_generated_mass(m_mass);
0130     particle->set_pid(pdg);
0131     particle->set_status(1);
0132
0133     event->add_particle(particle);
0134
0135     particle->add_attribute("pg_seq",
0136                             std::make_shared<HepMC3::UIntAttribute>(ip));
0137     particle->add_attribute("acts",
0138                             std::make_shared<HepMC3::BoolAttribute>(true));
0139
0140     primaryVertex->add_particle_out(particle);
0141   }
0142
0143   return event;
0144 }
0145
0146 }  // namespace ActsExamples