trento/src/nucleon.cxx

0001 // TRENTO: Reduced Thickness Event-by-event Nuclear Topology
0002 // Copyright 2015 Jonah E. Bernhard, J. Scott Moreland
0003 // TRENTO3D: Three-dimensional extension of TRENTO by Weiyao Ke
0004 // MIT License
0005
0006 #include "nucleon.h"
0007
0008 #include <cmath>
0009 #include <limits>
0010 #include <random>
0011 #include <stdexcept>
0012
0013 #include <boost/math/constants/constants.hpp>
0014 #include <boost/math/special_functions/expint.hpp>
0015 #include <boost/math/tools/roots.hpp>
0016 #include <boost/program_options/variables_map.hpp>
0017
0018 #include "fwd_decl.h"
0019
0020 namespace trento {
0021
0022 namespace {
0023
0024 // Create ctor parameters for unit mean std::gamma_distribution.
0025 //   mean = alpha*beta == 1  ->  beta = 1/alpha
0026 // Used below in NucleonProfile ctor initializer list.
0027
0028 template <typename RealType> using param_type =
0029   typename std::gamma_distribution<RealType>::param_type;
0030
0031 template <typename RealType>
0032 param_type<RealType> gamma_param_unit_mean(RealType alpha = 1.) {
0033   return param_type<RealType>{alpha, 1./alpha};
0034 }
0035
0036 // These constants define distances in terms of the width of the nucleon profile
0037 // Gaussian thickness function.
0038
0039 // Truncation radius of the thickness function.
0040 constexpr double trunc_radius_widths = 5.;
0041
0042 // Maximum impact parameter for participation.
0043 constexpr double max_impact_widths = 6.;
0044
0045 // Trivial helper function.
0046 template <typename T>
0047 constexpr T sqr(T value) {
0048   return value * value;
0049 }
0050
0051 // Inelastic nucleon-nucleon cross section as function of beam energy sqrt(s)
0052 // Fit coefficients explained in the docs.
0053 double cross_sec_from_energy(double sqrts) {
0054   auto a = 3.1253;
0055   auto b = 0.1280;
0056   auto c = 2.0412;
0057   auto d = 1.8231;
0058   return a + b * pow(std::log(sqrts) - c, d);
0059 }
0060
0061 // Determine the cross section parameter for sampling participants.
0062 // See section "Fitting the cross section" in the online docs.
0063 double compute_cross_sec_param(const VarMap& var_map) {
0064   // Read parameters from the configuration.
0065
0066   // Use manual inelastic nucleon-nucleon cross section if specified.
0067   // Otherwise default to extrapolated cross section.
0068   auto sigma_nn = var_map["cross-section"].as<double>();
0069   if (sigma_nn < 0) {
0070     sigma_nn = cross_sec_from_energy(var_map["beam-energy"].as<double>());
0071   }
0072   auto width = var_map["nucleon-width"].as<double>();
0073
0074   // Initialize arguments for boost root finding function.
0075
0076   // Bracket min and max.
0077   auto a = -10.;
0078   auto b = 20.;
0079
0080   // Tolerance function.
0081   // Require 3/4 of double precision.
0082   math::tools::eps_tolerance<double> tol{
0083     (std::numeric_limits<double>::digits * 3) / 4};
0084
0085   // Maximum iterations.
0086   // This is overkill -- in testing only 10-20 iterations were required
0087   // (but no harm in overestimating).
0088   boost::uintmax_t max_iter = 1000;
0089
0090   // The right-hand side of the equation.
0091   auto rhs = sigma_nn / (4 * math::double_constants::pi * sqr(width));
0092
0093   // This quantity appears a couple times in the equation.
0094   auto c = sqr(max_impact_widths) / 4;
0095
0096   try {
0097     auto result = math::tools::toms748_solve(
0098       [&rhs, &c](double x) {
0099         using std::exp;
0100         using math::expint;
0101         return c - expint(-exp(x)) + expint(-exp(x-c)) - rhs;
0102       },
0103       a, b, tol, max_iter);
0104
0105     return .5*(result.first + result.second);
0106   }
0107   catch (const std::domain_error&) {
0108     // Root finding fails for very small nucleon widths, w^2/sigma_nn < ~0.01.
0109     throw std::domain_error{
0110       "unable to fit cross section -- nucleon width too small?"};
0111   }
0112 }
0113
0114 }  // unnamed namespace
0115
0116 NucleonProfile::NucleonProfile(const VarMap& var_map)
0117     : width_sqr_(sqr(var_map["nucleon-width"].as<double>())),
0118       trunc_radius_sqr_(sqr(trunc_radius_widths)*width_sqr_),
0119       max_impact_sqr_(sqr(max_impact_widths)*width_sqr_),
0120       neg_one_div_two_width_sqr_(-.5/width_sqr_),
0121       neg_one_div_four_width_sqr_(-.25/width_sqr_),
0122       one_div_four_pi_(0.5*math::double_constants::one_div_two_pi),
0123       cross_sec_param_(compute_cross_sec_param(var_map)),
0124       fast_exp_(-.5*sqr(trunc_radius_widths), 0., 1000),
0125       fluct_dist_(gamma_param_unit_mean(var_map["fluctuation"].as<double>())),
0126       prefactor_(math::double_constants::one_div_two_pi/width_sqr_),
0127       with_ncoll_(var_map["ncoll"].as<bool>())
0128 {}
0129
0130 }  // namespace trento