ActsFatras/Utilities/LandauDistribution.hpp

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 #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 location_ Location parameter value
0035     /// @param scale_ Scale parameter value
0036     param_type(double location_, double scale_)
0037         : location(location_), scale(scale_) {}
0038     // Explicitly defaulted construction and assignment
0039     param_type() = default;
0040     /// @brief Copy constructor
0041     param_type(const param_type &) = default;
0042     /// @brief Move constructor
0043     param_type(param_type &&) = default;
0044     /// @brief Copy assignment operator
0045     /// @return Reference to this parameter object
0046     param_type &operator=(const param_type &) = default;
0047     /// @brief Move assignment operator
0048     /// @return Reference to this parameter object
0049     param_type &operator=(param_type &&) = default;
0050
0051     /// Parameters should be EqualityComparable
0052     friend bool operator==(const param_type &lhs, const param_type &rhs) {
0053       return (lhs.location == rhs.location) && (lhs.scale == rhs.scale);
0054     }
0055   };
0056   /// The type of the generated values.
0057   using result_type = double;
0058
0059   /// Construct directly from the distribution parameters.
0060   /// @param location Location parameter of the distribution
0061   /// @param scale Scale parameter of the distribution
0062   LandauDistribution(double location, double scale) : m_cfg(location, scale) {}
0063   /// Construct from a parameter object.
0064   /// @param cfg Parameter configuration object
0065   explicit LandauDistribution(const param_type &cfg) : m_cfg(cfg) {}
0066   // Explicitlely defaulted construction and assignment
0067   LandauDistribution() = default;
0068   /// @brief Copy constructor
0069   LandauDistribution(const LandauDistribution &) = default;
0070   /// @brief Move constructor
0071   LandauDistribution(LandauDistribution &&) = default;
0072   /// @brief Copy assignment operator
0073   /// @return Reference to this distribution object
0074   LandauDistribution &operator=(const LandauDistribution &) = default;
0075   /// @brief Move assignment operator
0076   /// @return Reference to this distribution object
0077   LandauDistribution &operator=(LandauDistribution &&) = default;
0078
0079   /// Construct from energy loss and FWHM
0080   ///
0081   /// @param energyLoss the (unscaled) most probable energy loss value
0082   /// @param energyLossFwhm the full width at half maximum of the energy loss
0083   /// @return A LandauDistribution configured with the derived location and scale
0084   static LandauDistribution fromFwhm(double energyLoss, double energyLossFwhm) {
0085     double scale = energyLossFwhm / kScaleFactorFwhm;
0086     double location = energyLoss - scale * kScaleFactorMPV;
0087     return LandauDistribution(location, scale);
0088   }
0089
0090   /// Reset any possible internal state. Noop, since there is no internal state.
0091   void reset() {}
0092   /// Return the currently configured distribution parameters.
0093   /// @return Current parameter configuration
0094   param_type param() const { return m_cfg; }
0095   /// Set the distribution parameters.
0096   /// @param cfg New parameter configuration to use
0097   void param(const param_type &cfg) { m_cfg = cfg; }
0098
0099   /// The minimum value the distribution generates.
0100   /// @return Minimum possible value (negative infinity)
0101   result_type min() const { return -std::numeric_limits<double>::infinity(); }
0102   /// The maximum value the distribution generates.
0103   /// @return Maximum possible value (positive infinity)
0104   result_type max() const { return std::numeric_limits<double>::infinity(); }
0105
0106   /// Generate a random number from the configured Landau distribution.
0107   /// @param generator Random number generator to use
0108   /// @return Random value from the Landau distribution
0109   template <typename Generator>
0110   result_type operator()(Generator &generator) {
0111     return (*this)(generator, m_cfg);
0112   }
0113   /// Generate a random number from the given Landau distribution.
0114   /// @param generator Random number generator to use
0115   /// @param params Distribution parameters to use for this generation
0116   /// @return Random value from the Landau distribution with given parameters
0117   template <typename Generator>
0118   result_type operator()(Generator &generator, const param_type &params) {
0119     const auto z = std::uniform_real_distribution<double>()(generator);
0120     return params.location + params.scale * quantile(z);
0121   }
0122
0123   /// Provide standard comparison operators
0124   friend bool operator==(const LandauDistribution &lhs,
0125                          const LandauDistribution &rhs) {
0126     return lhs.m_cfg == rhs.m_cfg;
0127   }
0128
0129  private:
0130   param_type m_cfg;
0131
0132   static double quantile(double z);
0133
0134   /// Scaling for most probable value
0135   static constexpr double kScaleFactorMPV = -0.222783;
0136   /// Scaling for FWHM
0137   static constexpr double kScaleFactorFwhm = 4.018646;
0138 };
0139
0140 }  // namespace ActsFatras