Skip to content

Commit

Permalink
draft
Browse files Browse the repository at this point in the history
  • Loading branch information
andiwand committed Oct 19, 2023
1 parent 25331d0 commit 689493e
Show file tree
Hide file tree
Showing 7 changed files with 152 additions and 77 deletions.
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
// This file is part of the Acts project.
//
// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
// Copyright (C) 2017-2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
Expand All @@ -18,22 +18,23 @@
#include <random>
#include <utility>

namespace {
double etaToCosTheta(double eta) {
const double theta = 2 * std::atan(std::exp(-eta));
return std::cos(theta);
}
} // namespace

ActsExamples::ParametricParticleGenerator::ParametricParticleGenerator(
const Config& cfg)
: m_cfg(cfg),
m_charge(cfg.charge.value_or(Acts::findCharge(m_cfg.pdg).value_or(0))),
m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg).value_or(0))),
// since we want to draw the direction uniform on the unit sphere, we must
// draw from cos(theta) instead of theta. see e.g.
// https://mathworld.wolfram.com/SpherePointPicking.html
m_cosThetaMin(std::cos(m_cfg.thetaMin)),
// ensure upper bound is included. see e.g.
// https://en.cppreference.com/w/cpp/numeric/random/uniform_real_distribution
m_cosThetaMax(std::nextafter(std::cos(m_cfg.thetaMax),
std::numeric_limits<double>::max())),
// in case we force uniform eta generation
m_etaMin(-std::log(std::tan(0.5 * m_cfg.thetaMin))),
m_etaMax(-std::log(std::tan(0.5 * m_cfg.thetaMax))) {}
m_cosThetaMin(etaToCosTheta(m_cfg.etaMin.get())),
m_cosThetaMax(etaToCosTheta(m_cfg.etaMax.get())),
m_charge(cfg.charge.value_or(Acts::findCharge(m_cfg.pdg.get()).value())),
m_mass(cfg.mass.value_or(Acts::findMass(m_cfg.pdg.get()).value())) {}

ActsExamples::SimParticleContainer
ActsExamples::ParametricParticleGenerator::operator()(RandomEngine& rng) {
Expand All @@ -45,17 +46,17 @@ ActsExamples::ParametricParticleGenerator::operator()(RandomEngine& rng) {
UniformIndex particleTypeChoice(0u, m_cfg.randomizeCharge ? 1u : 0u);
// (anti-)particle choice is one random draw but defines two properties
const Acts::PdgParticle pdgChoices[] = {
m_cfg.pdg,
static_cast<Acts::PdgParticle>(-m_cfg.pdg),
m_cfg.pdg.get(),
static_cast<Acts::PdgParticle>(-m_cfg.pdg.get()),
};
const double qChoices[] = {
m_charge,
-m_charge,
};
UniformReal phiDist(m_cfg.phiMin, m_cfg.phiMax);
UniformReal cosThetaDist(m_cosThetaMin, m_cosThetaMax);
UniformReal etaDist(m_etaMin, m_etaMax);
UniformReal pDist(m_cfg.pMin, m_cfg.pMax);
UniformReal etaDist(m_cfg.etaMin.get(), m_cfg.etaMax.get());
UniformReal pDist(m_cfg.pMin.get(), m_cfg.pMax.get());

SimParticleContainer particles;
particles.reserve(m_cfg.numParticles);
Expand Down
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
// This file is part of the Acts project.
//
// Copyright (C) 2017-2020 CERN for the benefit of the Acts project
// Copyright (C) 2017-2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
Expand All @@ -12,6 +12,7 @@
#include "Acts/Definitions/Units.hpp"
#include "ActsExamples/EventData/SimParticle.hpp"
#include "ActsExamples/Framework/RandomNumbers.hpp"
#include "ActsExamples/Framework/Required.hpp"
#include "ActsExamples/Generators/EventGenerator.hpp"

#include <array>
Expand All @@ -34,27 +35,28 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator {
/// Low, high (exclusive) for the transverse direction angle.
double phiMin = -M_PI;
double phiMax = M_PI;
/// Low, high (inclusive) for the longitudinal direction angle.
///
/// This intentionally uses theta instead of eta so it can represent the
/// full direction space with finite values.

/// Low, high (exclusive) for the pseudo rapidity.
///
/// @note This is the standard generation, for detector performance
/// classification, where a flat distribution in eta can be useful,
/// this can be set by the etaUniform flag;
///
double thetaMin = std::numeric_limits<double>::min();
double thetaMax = M_PI - std::numeric_limits<double>::epsilon();
bool etaUniform = false;
Required<double> etaMin;
Required<double> etaMax;
Required<bool> etaUniform;

/// Low, high (exclusive) for absolute/transverse momentum.
double pMin = 1 * Acts::UnitConstants::GeV;
double pMax = 10 * Acts::UnitConstants::GeV;
/// Indicate if the momentum referse to transverse momentum
bool pTransverse = false;
Required<double> pMin;
Required<double> pMax;
/// Indicate if the momentum referse to transverse momentum.
Required<bool> pTransverse;

/// (Absolute) PDG particle number to identify the particle type.
Acts::PdgParticle pdg = Acts::PdgParticle::eMuon;
Required<Acts::PdgParticle> pdg;
/// Randomize the charge and flip the PDG particle number sign accordingly.
bool randomizeCharge = false;
Required<bool> randomizeCharge;

/// Number of particles.
size_t numParticles = 1;

Expand All @@ -71,13 +73,12 @@ class ParametricParticleGenerator : public EventGenerator::ParticlesGenerator {

private:
Config m_cfg;
// will be automatically set from PDG data tables
double m_charge;
double m_mass;
double m_cosThetaMin;
double m_cosThetaMax;
double m_etaMin;
double m_etaMax;

double m_cosThetaMin{};
double m_cosThetaMax{};

double m_charge{};
double m_mass{};
};

} // namespace ActsExamples
86 changes: 86 additions & 0 deletions Examples/Framework/include/ActsExamples/Framework/Required.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
// This file is part of the Acts project.
//
// Copyright (C) 2023 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include <iostream>
#include <optional>
#include <stdexcept>

namespace ActsExamples {

/// A simple std::optional like container which can only be overwritten by a
/// defined value.
template <typename T>
class Required {
public:
using value_type = T;

constexpr Required() = default;

constexpr Required(const T& data) : m_data(data) {}
constexpr Required(T&& data) : m_data(std::move(data)) {}

constexpr Required(const Required& other) {
std::cout << "copy construct" << std::endl;
check(other.m_data);
m_data = other.m_data;
}
constexpr Required(Required&& other) {
std::cout << "move construct" << std::endl;
check(other.m_data);
m_data = std::move(other.m_data);
}

Required& operator=(const Required& other) {
std::cout << "copy assignment" << std::endl;
check(other.m_data);
m_data = other.m_data;
return *this;
}
Required& operator=(Required&& other) {
std::cout << "move assignment" << std::endl;
check(other.m_data);
m_data = std::move(other.m_data);
return *this;
}

T& operator*() noexcept { return get(); }
const T& operator*() const noexcept { return get(); }
T* operator->() noexcept { return &get(); }
const T* operator->() const noexcept { return &get(); }

constexpr explicit operator bool() const noexcept { return has(); }

constexpr bool has() const noexcept { return m_data.has_value(); }

constexpr T& get() & { return m_data.value(); }
constexpr const T& get() const& { return m_data.value(); }
constexpr T get() && { return m_data.value(); }

template <class... Args>
constexpr T& emplace(Args&&... args) {
return m_data.emplace(std::forward<Args>(args)...);
}

template <class U, class... Args>
constexpr T& emplace(std::initializer_list<U> ilist, Args&&... args) {
return m_data.emplace(ilist, std::forward<Args>(args)...);
}

private:
std::optional<T> m_data;

static void check(const std::optional<T>& data) {
if (!data) {
throw std::runtime_error("data required");
}
}
};

} // namespace ActsExamples
10 changes: 10 additions & 0 deletions Examples/Python/include/Acts/Plugins/Python/Utilities.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,23 @@
#pragma once

#include "Acts/Utilities/TypeTraits.hpp"
#include "ActsExamples/Framework/Required.hpp"

#include <string>
#include <unordered_map>

#include <boost/preprocessor/seq/for_each.hpp>
#include <boost/preprocessor/variadic/to_seq.hpp>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>

namespace PYBIND11_NAMESPACE {
namespace detail {
template <typename T>
struct type_caster<ActsExamples::Required<T>>
: optional_caster<ActsExamples::Required<T>> {};
} // namespace detail
} // namespace PYBIND11_NAMESPACE

namespace Acts::Python {

Expand Down
53 changes: 15 additions & 38 deletions Examples/Python/src/Generators.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@
#include "Acts/Utilities/Logger.hpp"
#include "ActsExamples/EventData/SimParticle.hpp"
#include "ActsExamples/Framework/RandomNumbers.hpp"
#include "ActsExamples/Framework/Required.hpp"
#include "ActsExamples/Generators/EventGenerator.hpp"
#include "ActsExamples/Generators/MultiplicityGenerators.hpp"
#include "ActsExamples/Generators/ParametricParticleGenerator.hpp"
Expand All @@ -36,16 +37,6 @@ class IReader;

namespace py = pybind11;

namespace {
double thetaToEta(double theta) {
assert(theta != 0);
return -1 * std::log(std::tan(theta / 2.));
}
double etaToTheta(double eta) {
return 2 * std::atan(std::exp(-eta));
}
} // namespace

namespace Acts::Python {

void addGenerators(Context& ctx) {
Expand Down Expand Up @@ -135,8 +126,8 @@ void addGenerators(Context& ctx) {
.def(py::init<>())
.def_readwrite("phiMin", &Config::phiMin)
.def_readwrite("phiMax", &Config::phiMax)
.def_readwrite("thetaMin", &Config::thetaMin)
.def_readwrite("thetaMax", &Config::thetaMax)
.def_readwrite("etaMin", &Config::etaMin)
.def_readwrite("etaMax", &Config::etaMax)
.def_readwrite("etaUniform", &Config::etaUniform)
.def_readwrite("pMin", &Config::pMin)
.def_readwrite("pMax", &Config::pMax)
Expand All @@ -147,41 +138,26 @@ void addGenerators(Context& ctx) {
.def_readwrite("mass", &Config::mass)
.def_readwrite("charge", &Config::charge)
.def_property(
"p",
[](Config& cfg) {
return std::pair{cfg.pMin, cfg.pMax};
},
[](Config& cfg, std::pair<double, double> value) {
"p", [](Config& cfg) { return std::tie(cfg.pMin, cfg.pMax); },
[](Config& cfg, std::pair<ActsExamples::Required<double>,
ActsExamples::Required<double>>
value) {
cfg.pMin = value.first;
cfg.pMax = value.second;
})
.def_property(
"phi",
[](Config& cfg) {
return std::pair{cfg.phiMin, cfg.phiMax};
},
"phi", [](Config& cfg) { return std::tie(cfg.phiMin, cfg.phiMax); },
[](Config& cfg, std::pair<double, double> value) {
cfg.phiMin = value.first;
cfg.phiMax = value.second;
})
.def_property(
"theta",
[](Config& cfg) {
return std::pair{cfg.thetaMin, cfg.thetaMax};
},
[](Config& cfg, std::pair<double, double> value) {
cfg.thetaMin = value.first;
cfg.thetaMax = value.second;
})
.def_property(
"eta",
[](Config& cfg) {
return std::pair{thetaToEta(cfg.thetaMin),
thetaToEta(cfg.thetaMax)};
},
[](Config& cfg, std::pair<double, double> value) {
cfg.thetaMin = etaToTheta(value.first);
cfg.thetaMax = etaToTheta(value.second);
"eta", [](Config& cfg) { return std::tie(cfg.etaMin, cfg.etaMax); },
[](Config& cfg, std::pair<ActsExamples::Required<double>,
ActsExamples::Required<double>>
value) {
cfg.etaMin = value.first;
cfg.etaMax = value.second;
});
}

Expand Down Expand Up @@ -211,4 +187,5 @@ void addGenerators(Context& ctx) {
py::arg("mean"))
.def_readwrite("mean", &ActsExamples::PoissonMultiplicityGenerator::mean);
}

} // namespace Acts::Python
4 changes: 2 additions & 2 deletions Examples/Run/Common/src/ParticleGunOptions.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,8 +77,8 @@ ActsExamples::Options::readParticleGunOptions(const Variables& vars) {
getRange("gen-eta", 1.0, etaMin, etaMax);

pgCfg.etaUniform = vars["gen-eta-uniform"].template as<bool>();
pgCfg.thetaMin = 2 * std::atan(std::exp(-etaMin));
pgCfg.thetaMax = 2 * std::atan(std::exp(-etaMax));
pgCfg.etaMin = etaMin;
pgCfg.etaMax = etaMax;
getRange("gen-mom-gev", 1_GeV, pgCfg.pMin, pgCfg.pMax);
pgCfg.pTransverse = vars["gen-mom-transverse"].template as<bool>();
pgCfg.pdg =
Expand Down
2 changes: 1 addition & 1 deletion Examples/Scripts/Python/full_chain_odd.py
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@
addParticleGun(
s,
MomentumConfig(1.0 * u.GeV, 10.0 * u.GeV, transverse=True),
EtaConfig(-3.0, 3.0),
EtaConfig(-3.0, 3.0, uniform=True),
PhiConfig(0.0, 360.0 * u.degree),
ParticleConfig(4, acts.PdgParticle.eMuon, randomizeCharge=True),
vtxGen=acts.examples.GaussianVertexGenerator(
Expand Down

0 comments on commit 689493e

Please sign in to comment.