Skip to content

feat: Harmonize acts and traccc parameter smearing #1171

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

feat: Harmonize acts and traccc parameter smearing #1171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
142 changes: 126 additions & 16 deletions core/include/traccc/utils/seed_generator.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
// Library include(s).
#include "traccc/edm/track_parameters.hpp"
#include "traccc/utils/particle.hpp"
#include "traccc/utils/trigonometric_helpers.hpp"

// detray include(s).
#include <detray/geometry/barcode.hpp>
Expand All @@ -28,14 +29,53 @@ struct seed_generator {
using algebra_type = typename detector_t::algebra_type;
using ctx_t = typename detector_t::geometry_context;

/// Configure the seed generator
struct config {
// Smearing parameters
/// Constant term of the loc0 resolution.
double sigma_loc0 = 20 * traccc::unit<scalar>::um;
/// Pt-dependent loc0 resolution of the form sigma_loc0 =
/// A*exp(-1.*abs(B)*pt).
double sigma_loc0_pT_a = 30 * traccc::unit<scalar>::um;
double sigma_loc0_pT_b = 0.3 / traccc::unit<scalar>::GeV;
/// Constant term of the loc1 resolution.
double sigma_loc1 = 20 * traccc::unit<scalar>::um;
/// Pt-dependent loc1 resolution of the form sigma_loc1 =
/// A*exp(-1.*abs(B)*pt).
double sigma_loc1_pT_a = 30 * traccc::unit<scalar>::um;
double sigma_loc1_pT_b = 0.3 / traccc::unit<scalar>::GeV;
/// Time resolution.
double sigma_time = 1 * traccc::unit<scalar>::ns;
/// Phi angular resolution.
double sigma_phi = 1 * traccc::unit<scalar>::degree;
/// Theta angular resolution.
double sigma_theta = 1 * traccc::unit<scalar>::degree;
/// Relative transverse momentum resolution.
double sigma_pT_rel = 0.05;

/// Optional. Initial sigmas for the track parameters which overwrites
/// the smearing params if set.
std::optional<std::array<double, e_bound_size>> initial_sigmas;
/// Initial sigma(q/pt) for the track parameters.
/// @note The resulting q/p sigma is added to the one in `initialSigmas`
double initialsigma_qopt =
0.1 * traccc::unit<scalar>::e / traccc::unit<scalar>::GeV;
/// Initial sigma(pt)/pt for the track parameters.
/// @note The resulting q/p sigma is added to the one in `initialSigmas`
double initialsigma_pT_rel = 0.1;

/// Inflation factors for the variances
std::array<scalar, e_bound_size> cov_inflation{1., 1., 1., 1., 1., 1.};
};

/// Constructor with detector
///
/// @param det input detector
/// @param stddevs standard deviations for parameter smearing
seed_generator(const detector_t& det,
const std::array<scalar, e_bound_size>& stddevs,
const std::size_t sd = 0, ctx_t ctx = {})
: m_detector(det), m_stddevs(stddevs), m_ctx(ctx) {
explicit seed_generator(const detector_t& det, const config& cfg = {},
const std::size_t sd = std::mt19937::default_seed,
ctx_t ctx = {})
: m_detector(det), m_ctx(ctx), m_cfg(cfg) {
m_generator.seed(static_cast<std::mt19937::result_type>(sd));
}

Expand Down Expand Up @@ -73,32 +113,102 @@ struct seed_generator {
sf);
}

for (std::size_t i = 0; i < e_bound_size; i++) {
// Call the smearing operator
(*this)(bound_param, ptc_type);

if (m_stddevs[i] != scalar{0}) {
bound_param[i] = std::normal_distribution<scalar>(
bound_param[i], m_stddevs[i])(m_generator);
}
return bound_param;
}

/// Seed generator operation
///
/// @param vertex vertex of particle
/// @param stddevs standard deviations for track parameter smearing
void operator()(bound_track_parameters<algebra_type>& bound_param,
const traccc::pdg_particle<scalar>& ptc_type) {

const scalar q{ptc_type.charge()};
const auto pos = bound_param.bound_local();
const auto time = bound_param.time();
const auto phi = bound_param.phi();
const auto theta = bound_param.theta();
const auto pt = bound_param.pT(q);
const auto qop = bound_param.qop();

// Compute momentum-dependent resolutions
const double sigma_loc_0 =
m_cfg.sigma_loc0 +
m_cfg.sigma_loc0_pT_a *
math::exp(-1.0 * math::fabs(m_cfg.sigma_loc0_pT_b) * pt);
const double sigma_loc_1 =
m_cfg.sigma_loc1 +
m_cfg.sigma_loc1_pT_a *
math::exp(-1.0 * math::fabs(m_cfg.sigma_loc1_pT_b) * pt);
// Shortcuts for other resolutions
const double sigma_qop = math::sqrt(
math::pow(m_cfg.sigma_pT_rel * qop, 2) +
math::pow(m_cfg.sigma_theta * (qop * math::tan(theta)), 2));

// Smear the position/time
// Note that we smear d0 and z0 in the perigee frame
const scalar smeared_loc0 = std::normal_distribution<scalar>(
pos[0], static_cast<scalar>(sigma_loc_0))(m_generator);
const scalar smeared_loc1 = std::normal_distribution<scalar>(
pos[1], static_cast<scalar>(sigma_loc_1))(m_generator);
bound_param.set_bound_local({smeared_loc0, smeared_loc1});

// Time
bound_param.set_time(std::normal_distribution<scalar>(
time, static_cast<scalar>(m_cfg.sigma_time))(m_generator));

// Smear direction angles phi,theta ensuring correct bounds
const scalar smeared_phi = std::normal_distribution<scalar>(
phi, static_cast<scalar>(m_cfg.sigma_phi))(m_generator);
const scalar smeared_theta = std::normal_distribution<scalar>(
theta, static_cast<scalar>(m_cfg.sigma_theta))(m_generator);
const auto [new_phi, new_theta] =
detail::wrap_phi_theta(smeared_phi, smeared_theta);
bound_param.set_phi(new_phi);
bound_param.set_theta(new_theta);

// Compute smeared q/p
bound_param.set_qop(std::normal_distribution<scalar>(
qop, static_cast<scalar>(sigma_qop))(m_generator));

// Build the track covariance matrix using the smearing sigmas
std::array<scalar, e_bound_size> sigmas{
static_cast<scalar>(sigma_loc_0),
static_cast<scalar>(sigma_loc_1),
static_cast<scalar>(m_cfg.sigma_phi),
static_cast<scalar>(m_cfg.sigma_theta),
static_cast<scalar>(sigma_qop),
static_cast<scalar>(m_cfg.sigma_time)};

for (std::size_t i = e_bound_loc0; i < e_bound_size; ++i) {
double sigma = sigmas[i];
double variance = sigma * sigma;

// Inflate the initial covariance
variance *= m_cfg.cov_inflation[i];

getter::element(bound_param.covariance(), i, i) =
m_stddevs[i] * m_stddevs[i];
static_cast<scalar>(variance);
}

assert(!bound_param.is_invalid());

return bound_param;
}

private:
// Random generator
/// Random generator
std::random_device m_rd{};
std::mt19937 m_generator{m_rd()};

// Detector object
/// Detector object
const detector_t& m_detector;
/// Standard deviations for parameter smearing
std::array<scalar, e_bound_size> m_stddevs;
/// Geometry context
ctx_t m_ctx;

/// Seed generator configuration
config m_cfg{};
};

} // namespace traccc
13 changes: 2 additions & 11 deletions examples/run/cpu/misaligned_truth_fitting_example.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,15 +110,6 @@ int main(int argc, char* argv[]) {
* Do the reconstruction
*****************************/

/// Standard deviations for seed track parameters
static constexpr std::array<scalar, e_bound_size> stddevs = {
0.03f * traccc::unit<scalar>::mm,
0.03f * traccc::unit<scalar>::mm,
0.017f,
0.017f,
0.001f / traccc::unit<scalar>::GeV,
1.f * traccc::unit<scalar>::ns};

// Fitting algorithm objects
// Alg0
traccc::fitting_config fit_cfg0(fitting_opts);
Expand All @@ -135,9 +126,9 @@ int main(int argc, char* argv[]) {

// Seed generators
traccc::seed_generator<host_detector_type> sg0(
host_det, stddevs, 0, fit_cfg0.propagation.context);
host_det, {}, 0, fit_cfg0.propagation.context);
traccc::seed_generator<host_detector_type> sg1(
host_det, stddevs, 0, fit_cfg1.propagation.context);
host_det, {}, 0, fit_cfg1.propagation.context);

// Iterate over events
for (auto event = input_opts.skip;
Expand Down
11 changes: 1 addition & 10 deletions examples/run/cpu/truth_finding_example.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -90,15 +90,6 @@ int seq_run(const traccc::opts::track_finding& finding_opts,
* Do the reconstruction
*****************************/

// Standard deviations for seed track parameters
static constexpr std::array<traccc::scalar, traccc::e_bound_size> stddevs =
{1e-4f * traccc::unit<traccc::scalar>::mm,
1e-4f * traccc::unit<traccc::scalar>::mm,
1e-3f,
1e-3f,
1e-4f / traccc::unit<traccc::scalar>::GeV,
1e-4f * traccc::unit<traccc::scalar>::ns};

// Propagation configuration
detray::propagation::config propagation_config(propagation_opts);

Expand Down Expand Up @@ -137,7 +128,7 @@ int seq_run(const traccc::opts::track_finding& finding_opts,
[&]<typename detector_traits_t>(
const typename detector_traits_t::host& det) {
traccc::seed_generator<typename detector_traits_t::host> sg(
det, stddevs);
det);
evt_data.generate_truth_candidates(
truth_track_candidates, truth_measurements, sg, host_mr,
truth_finding_opts.m_pT_min);
Expand Down
11 changes: 1 addition & 10 deletions examples/run/cpu/truth_fitting_example.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,15 +91,6 @@ int main(int argc, char* argv[]) {
* Do the reconstruction
*****************************/

/// Standard deviations for seed track parameters
static constexpr std::array<scalar, e_bound_size> stddevs = {
0.03f * traccc::unit<scalar>::mm,
0.03f * traccc::unit<scalar>::mm,
0.017f,
0.017f,
0.001f / traccc::unit<scalar>::GeV,
1.f * traccc::unit<scalar>::ns};

// Fitting algorithm object
traccc::fitting_config fit_cfg(fitting_opts);
fit_cfg.propagation = propagation_opts;
Expand Down Expand Up @@ -127,7 +118,7 @@ int main(int argc, char* argv[]) {
const typename detector_traits_t::host& det) {
// Seed generator
traccc::seed_generator<typename detector_traits_t::host> sg(
det, stddevs);
det);
evt_data.generate_truth_candidates(
truth_track_candidates, truth_measurements, sg, host_mr);
});
Expand Down
11 changes: 1 addition & 10 deletions examples/run/cuda/truth_finding_example_cuda.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,15 +122,6 @@ int seq_run(const traccc::opts::track_finding& finding_opts,
traccc::buffer_from_host_detector(polymorphic_detector, device_mr,
async_copy);

// Standard deviations for seed track parameters
static constexpr std::array<traccc::scalar, traccc::e_bound_size> stddevs =
{1e-4f * traccc::unit<traccc::scalar>::mm,
1e-4f * traccc::unit<traccc::scalar>::mm,
1e-3f,
1e-3f,
1e-4f / traccc::unit<traccc::scalar>::GeV,
1e-4f * traccc::unit<traccc::scalar>::ns};

// Propagation configuration
detray::propagation::config propagation_config(propagation_opts);

Expand Down Expand Up @@ -175,7 +166,7 @@ int seq_run(const traccc::opts::track_finding& finding_opts,
const typename detector_traits_t::host& det) {
// Seed generator
traccc::seed_generator<typename detector_traits_t::host> sg(
det, stddevs);
det);
evt_data.generate_truth_candidates(
truth_track_candidates, truth_measurements, sg, host_mr,
truth_finding_opts.m_pT_min);
Expand Down
11 changes: 1 addition & 10 deletions examples/run/cuda/truth_fitting_example_cuda.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -120,15 +120,6 @@ int main(int argc, char* argv[]) {
traccc::buffer_from_host_detector(polymorphic_detector, device_mr,
async_copy);

/// Standard deviations for seed track parameters
static constexpr std::array<scalar, e_bound_size> stddevs = {
0.03f * traccc::unit<scalar>::mm,
0.03f * traccc::unit<scalar>::mm,
0.017f,
0.017f,
0.001f / traccc::unit<scalar>::GeV,
1.f * traccc::unit<scalar>::ns};

// Fitting algorithm object
traccc::fitting_config fit_cfg(fitting_opts);
fit_cfg.propagation = propagation_opts;
Expand Down Expand Up @@ -160,7 +151,7 @@ int main(int argc, char* argv[]) {
const typename detector_traits_t::host& det) {
// Seed generator
traccc::seed_generator<typename detector_traits_t::host> sg(
det, stddevs);
det);
evt_data.generate_truth_candidates(
truth_track_candidates, truth_measurements, sg, host_mr);
});
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_ckf_combinatorics_telescope.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -114,8 +114,7 @@ TEST_P(CpuCkfCombinatoricsTelescopeTests, Run) {
*****************************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Finding algorithm configuration
traccc::finding_config cfg_no_limit;
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_ckf_sparse_tracks_telescope.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,8 +126,7 @@ TEST_P(CkfSparseTrackTelescopeTests, Run) {
*****************************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Finding algorithm configuration
typename traccc::finding_config cfg;
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_kalman_fitter_hole_count.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,8 +121,7 @@ TEST_P(KalmanFittingHoleCountTests, Run) {
***************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Fitting algorithm object
traccc::fitting_config fit_cfg;
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_kalman_fitter_momentum_resolution.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -142,8 +142,7 @@ TEST_P(KalmanFittingMomentumResolutionTests, Run) {
***************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Fitting algorithm object
traccc::fitting_config fit_cfg;
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_kalman_fitter_telescope.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,8 +122,7 @@ TEST_P(KalmanFittingTelescopeTests, Run) {
***************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Fitting algorithm object
traccc::fitting_config fit_cfg;
Expand Down
3 changes: 1 addition & 2 deletions tests/cpu/test_kalman_fitter_wire_chamber.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -124,8 +124,7 @@ TEST_P(KalmanFittingWireChamberTests, Run) {
***************/

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Fitting algorithm object
traccc::fitting_config fit_cfg;
Expand Down
3 changes: 1 addition & 2 deletions tests/cuda/test_ckf_combinatorics_telescope.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,8 +131,7 @@ TEST_P(CudaCkfCombinatoricsTelescopeTests, Run) {
vecmem::cuda::async_copy copy{stream.cudaStream()};

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Finding algorithm configuration
typename traccc::cuda::combinatorial_kalman_filter_algorithm::config_type
Expand Down
3 changes: 1 addition & 2 deletions tests/cuda/test_ckf_toy_detector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,8 +133,7 @@ TEST_P(CkfToyDetectorTests, Run) {
vecmem::cuda::async_copy copy{stream.cudaStream()};

// Seed generator
seed_generator<host_detector_type> sg(detector.as<detector_traits>(),
stddevs);
seed_generator<host_detector_type> sg(detector.as<detector_traits>());

// Finding algorithm configuration
typename traccc::cuda::combinatorial_kalman_filter_algorithm::config_type
Expand Down
2 changes: 1 addition & 1 deletion tests/cuda/test_kalman_fitter_telescope.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,7 +141,7 @@ TEST_P(KalmanFittingTelescopeTests, Run) {

// Seed generator
seed_generator<host_detector_type> sg(
polymorphic_detector.as<detector_traits>(), stddevs);
polymorphic_detector.as<detector_traits>());

// Fitting algorithm object
traccc::cuda::kalman_fitting_algorithm::config_type fit_cfg;
Expand Down
Loading
Loading