ActsFatras/Utilities/LandauDistribution.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include <limits>
0012 #include <random>
0013
0014 namespace ActsFatras {
0015
0016 /// Draw random numbers from a Landau distribution.
0017 ///
0018 /// Implements the same interface as the standard library distributions.
0019 class LandauDistribution {
0020  public:
0021   /// Parameter struct that contains all distribution parameters.
0022   struct param_type {
0023     /// Parameters must link back to the host distribution.
0024     using distribution_type = LandauDistribution;
0025
0026     /// Location parameter.
0027     ///
0028     /// @warning This is neither the mean nor the most probable value.
0029     double location = 0.0;
0030     /// Scale parameter.
0031     double scale = 1.0;
0032
0033     /// Construct from parameters.
0034     param_type(double location_, double scale_)
0035         : location(location_), scale(scale_) {}
0036     // Explicitlely defaulted construction and assignment
0037     param_type() = default;
0038     param_type(const param_type &) = default;
0039     param_type(param_type &&) = default;
0040     param_type &operator=(const param_type &) = default;
0041     param_type &operator=(param_type &&) = default;
0042
0043     /// Parameters should be EqualityComparable
0044     friend bool operator==(const param_type &lhs, const param_type &rhs) {
0045       return (lhs.location == rhs.location) && (lhs.scale == rhs.scale);
0046     }
0047     friend bool operator!=(const param_type &lhs, const param_type &rhs) {
0048       return !(lhs == rhs);
0049     }
0050   };
0051   /// The type of the generated values.
0052   using result_type = double;
0053
0054   /// Construct directly from the distribution parameters.
0055   LandauDistribution(double location, double scale) : m_cfg(location, scale) {}
0056   /// Construct from a parameter object.
0057   LandauDistribution(const param_type &cfg) : m_cfg(cfg) {}
0058   // Explicitlely defaulted construction and assignment
0059   LandauDistribution() = default;
0060   LandauDistribution(const LandauDistribution &) = default;
0061   LandauDistribution(LandauDistribution &&) = default;
0062   LandauDistribution &operator=(const LandauDistribution &) = default;
0063   LandauDistribution &operator=(LandauDistribution &&) = default;
0064
0065   /// Reset any possible internal state. Noop, since there is no internal state.
0066   void reset() {}
0067   /// Return the currently configured distribution parameters.
0068   param_type param() const { return m_cfg; }
0069   /// Set the distribution parameters.
0070   void param(const param_type &cfg) { m_cfg = cfg; }
0071
0072   /// The minimum value the distribution generates.
0073   result_type min() const { return -std::numeric_limits<double>::infinity(); }
0074   /// The maximum value the distribution generates.
0075   result_type max() const { return std::numeric_limits<double>::infinity(); }
0076
0077   /// Generate a random number from the configured Landau distribution.
0078   template <typename Generator>
0079   result_type operator()(Generator &generator) {
0080     return (*this)(generator, m_cfg);
0081   }
0082   /// Generate a random number from the given Landau distribution.
0083   template <typename Generator>
0084   result_type operator()(Generator &generator, const param_type &params) {
0085     const auto z = std::uniform_real_distribution<double>()(generator);
0086     return params.location + params.scale * quantile(z);
0087   }
0088
0089   /// Provide standard comparison operators
0090   friend bool operator==(const LandauDistribution &lhs,
0091                          const LandauDistribution &rhs) {
0092     return lhs.m_cfg == rhs.m_cfg;
0093   }
0094   friend bool operator!=(const LandauDistribution &lhs,
0095                          const LandauDistribution &rhs) {
0096     return !(lhs == rhs);
0097   }
0098
0099  private:
0100   param_type m_cfg;
0101
0102   static double quantile(double z);
0103 };
0104
0105 }  // namespace ActsFatras