smash::detail Namespace Reference

specialize a type for all of the STL containers. More...

Classes
struct	KeyTraits
	Class template to store Key traits outside the Key class, allowing for reuse both in the Key class itself and in helper implementation details. More...
struct	is_stl_container
	Implementation of the type trait to infer if a type is an STL container. More...
struct	is_stl_container< std::vector< Args... > >
	Trait specialization for std::vector. More...
struct	is_stl_container< std::set< Args... > >
	Trait specialization for std::set. More...
struct	is_stl_container< std::map< Args... > >
	Trait specialization for std::map. More...

Functions
void	sample_manybody_phasespace_impl (double sqrts, const ParticleTypePtrList &types, std::vector< FourVector > &sampled_momenta)
	Implementation of the full n-body phase-space sampling (masses, momenta, angles) in the center-of-mass frame for the final state particles, using the M-method from CERN-68-15, paragraph 9.6. More...
void	sample_manybody_phasespace_MCMC (const ParticleTypePtrList &types, std::vector< FourVector > &sampled_momenta)
	Metropolis–Hastings sampling of the many body phase space, starting from an initial guess in sampled_momenta that is assumed appropriate. More...
double	interpolate_trilinear (double ax, double ay, double az, double f1, double f2, double f3, double f4, double f5, double f6, double f7, double f8)
	Perform a trilinear 1st order interpolation. More...
template<typename T >
const KeyTraits< T >::validator_type &	get_default_validator () noexcept
	Function template to get a default trivial validator. More...
template<typename T >
constexpr auto	type_name ()
	Get type of variable as string in a human-readable way. More...
template<typename N , typename = std::enable_if_t<std::is_floating_point_v<N>>>
bool	almost_equal_knuthish (const N x, const N y, const N epsilon, const N threshold=N{0.0}) noexcept
	Compare whether two floating-point numbers are approximately equal à la Knuth up to a given tolerance. More...
template<typename Converter , class Range , std::enable_if_t< std::is_same_v< Range, Particles >||std::is_same_v< Range, ParticleList >, bool > = true>
void	write_in_chunk_impl (const Range &particles, const OutputFormatter< Converter > &formatter, std::function< void(const typename Converter::type &)> write, std::size_t max_buffer_bytes=1 '000 '000 '000)
	Writes particle data in multiple chunks if the total buffer size exceeds a predefined maximum. More...

Detailed Description

specialize a type for all of the STL containers.

Function Documentation

sample_manybody_phasespace_impl()

void smash::detail::sample_manybody_phasespace_impl	(	double	sqrts,
		const ParticleTypePtrList &	types,
		std::vector< FourVector > &	sampled_momenta
	)

Implementation of the full n-body phase-space sampling (masses, momenta, angles) in the center-of-mass frame for the final state particles, using the M-method from CERN-68-15, paragraph 9.6.

The algorithm proceeds in two stages:

1. Generate invariant masses \(M_{12}, M_{123}, M_{1234}, \ldots\) from the measure

   \[ dM_{12}\, dM_{123}\, dM_{1234}\, \cdots, \]

   while respecting non-trivial kinematic limits.

   Introduce shifted variables

   \[ T_{12} = M_{12} - (m_1 + m_2),\qquad T_{123} = M_{123} - (m_1 + m_2 + m_3),\ \ldots \]

   and sample uniformly under the ordering constraint

   \[ 0 \le T_{12} \le T_{123} \le T_{1234} \le \cdots \le \sqrt{s} - \sum_i m_i. \]

   A practical trick is to draw values uniformly in \([0,\,\sqrt{s} - \sum_i m_i]\) and sort them.

2. Accept or reject each invariant-mass configuration with weight proportional to

   \[ R_2(\sqrt{s}, M_{n-1}, m_n) \times R_2(M_{n-1}, M_{n-2}, m_{n-1}) \times \cdots \times R_2(M_2, m_1, m_2) \times \prod_i M_i. \]

The maximum weight is estimated heuristically; following an idea by Scott Pratt, it is expected near

\[ T_{12} = T_{123} = T_{1234} = \cdots = \frac{\sqrt{s} - \sum_i m_i}{n - 1}. \]

Definition at line 410 of file action.cc.

                                                                             {
  const size_t n = types.size();
  assert(n > 1);
  sampled_momenta.resize(n);
  std::vector<double> masses{}, masses_sum(n), Minv(n);
  // Maximum estimate is rough and can be wrong. We multiply it by additional
  // factor to be on the safer side, and increase it if needed.
  double safety_factor = 1.1 + (n - 2) * 0.2;
  int rejection_counter = 0;
  constexpr int rejection_limit = 200;
  double available_energy =
      sqrts - std::accumulate(types.begin(), types.end(), 0.0,
                              [](double sum, const ParticleTypePtr &type) {
                                return sum + type->min_mass_spectral();
                              });
  double acceptance = 1;
  do {
    double random_01;
    // This loop increases the maximum if the safety_factor is too small
    safety_factor *= std::sqrt(acceptance);
    do {
      // Mass sampling from spectral functions
      double weight_sqr_max = safety_factor * safety_factor;
      masses.clear();
      masses_sum.clear();
      for (const auto &type : types) {
        if (type->is_stable()) {
          masses.push_back(type->mass());
        } else {
          masses.push_back(
              type->sample_breit_wigner_spectral_function(available_energy));
          const double max_ratio = std::max(
              type->max_ratio_spectral_full_to_breit_wigner(),
              type->ratio_spectral_full_to_breit_wigner(available_energy));
          weight_sqr_max *= max_ratio * max_ratio;
        }
      }
      // Arrange a convenient vector of m1, m1 + m2, m1 + m2 + m3, ...
      std::partial_sum(masses.begin(), masses.end(),
                       std::back_inserter(masses_sum));
      const double masses_sum_all = masses_sum[n - 1];
      const double Ekin_share = (sqrts - masses_sum_all) / (n - 1);
      for (size_t i = 1; i < n; i++) {
        // This maximum estimate idea is due Scott Pratt: maximum should be
        // roughly at equal kinetic energies
        weight_sqr_max *=
            pCM_sqr(i * Ekin_share + masses_sum[i],
                    (i - 1) * Ekin_share + masses_sum[i - 1], masses[i]);
      }
      // Generate invariant masses of 1, 12, 123, 1234, etc.
      // Minv = {m1, M12, M123, ..., M123n-1, sqrts}
      Minv[0] = 0.0;
      Minv[n - 1] = sqrts - masses_sum_all;
      for (size_t i = 1; i < n - 1; i++) {
        Minv[i] = random::uniform(0.0, sqrts - masses_sum_all);
      }
      std::sort(Minv.begin(), Minv.end());
      for (size_t i = 0; i < n; i++) {
        Minv[i] += masses_sum[i];
      }
 
      double weight_sqr = 1;
      for (size_t i = 0; i < n; i++) {
        const double ratio =
            types[i]->is_stable()
                ? 1
                : types[i]->ratio_spectral_full_to_breit_wigner(masses[i]);
        weight_sqr *= ratio * ratio;
      }
      for (size_t i = 1; i < n; i++) {
        weight_sqr *= pCM_sqr(Minv[i], Minv[i - 1], masses[i]);
      }
      acceptance = weight_sqr / weight_sqr_max;
      rejection_counter++;
      random_01 = random::canonical();
    } while (acceptance < random_01 * random_01 &&
             rejection_counter < rejection_limit);
    if (acceptance > 1) {
      logg[LAction].debug()
          << "sample_manybody_phasespace_impl: alarm, weight > 1, w^2 = "
          << acceptance << ". Increasing safety factor.";
    }
  } while (acceptance > 1 && rejection_counter < rejection_limit);
 
  // Boost particles to the right frame
  std::vector<ThreeVector> beta(n);
  for (size_t i = n - 1; i > 0; i--) {
    const double pcm = pCM(Minv[i], Minv[i - 1], masses[i]);
    Angles phitheta;
    phitheta.distribute_isotropically();
    const ThreeVector isotropic_unitvector = phitheta.threevec();
    sampled_momenta[i] =
        FourVector(std::sqrt(masses[i] * masses[i] + pcm * pcm),
                   pcm * isotropic_unitvector);
    if (i >= 2) {
      beta[i - 2] = pcm * isotropic_unitvector /
                    std::sqrt(pcm * pcm + Minv[i - 1] * Minv[i - 1]);
    }
    if (i == 1) {
      sampled_momenta[0] =
          FourVector(std::sqrt(masses[0] * masses[0] + pcm * pcm),
                     -pcm * isotropic_unitvector);
    }
  }
 
  for (size_t i = 0; i < n - 2; i++) {
    // After each boost except the last one the sum of 3-momenta should be 0
    FourVector ptot = FourVector(0.0, 0.0, 0.0, 0.0);
    for (size_t j = 0; j <= i + 1; j++) {
      ptot += sampled_momenta[j];
    }
    logg[LAction].debug() << "Total momentum of 0.." << i + 1 << " = "
                          << ptot.threevec() << " and should be (0, 0, 0).";
 
    // Boost the first i+1 particles to the next CM frame
    for (size_t j = 0; j <= i + 1; j++) {
      sampled_momenta[j] = sampled_momenta[j].lorentz_boost(beta[i]);
    }
  }
 
  FourVector ptot_all = FourVector(0.0, 0.0, 0.0, 0.0);
  for (size_t j = 0; j < n; j++) {
    ptot_all += sampled_momenta[j];
  }
  logg[LAction].debug() << "Total 4-momentum = " << ptot_all << ", should be ("
                        << sqrts << ", 0, 0, 0)";
  if (rejection_counter >= rejection_limit) {
    logg[LAction].warn()
        << "Failed to sample kinematically correct 4-momenta of "
        << std::accumulate(types.begin(), types.end(), std::string{},
                           [](const std::string &a, const ParticleTypePtr &b) {
                             return a + b->name();
                           })
        << " with energy " << sqrts
        << " GeV.\n Using MCMC fallback, which conserves energy and momentum "
           "but may deviate slightly from the correct spectral distribution.";
    sample_manybody_phasespace_MCMC(types, sampled_momenta);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample_manybody_phasespace_MCMC()

void smash::detail::sample_manybody_phasespace_MCMC	(	const ParticleTypePtrList &	types,
		std::vector< FourVector > &	sampled_momenta
	)

Metropolis–Hastings sampling of the many body phase space, starting from an initial guess in sampled_momenta that is assumed appropriate.

The algorithm works by repeatedly picking a random pair of particles, resampling their masses from the spectral functions, and adjusting their momenta accordingly in the CM frame of the pair, which conserves energy and momentum. This is done for a fixed number of iterations heuristically chosen (200), but no systematic analysis was done.

The function is for now used as a fallback for the rejection algorithm in sample_manybody_phasespace_impl and takes its initial guess from there.

Definition at line 552 of file action.cc.

                                                                             {
  constexpr int monte_carlo_iterations = 200;
  // sampled_momenta already contains an initial guess
  const int n = types.size();
  int idx1, idx2;
 
  assert(sampled_momenta.size() == static_cast<size_t>(n));
  // It is enough to check that the energy is a positive number
  const bool all_zero =
      std::all_of(sampled_momenta.begin(), sampled_momenta.end(),
                  [](const FourVector &p) { return p[0] < really_small; });
  if (all_zero) {
    throw std::runtime_error(
        "All initial momenta for the MCMC algorithm manybody phase space "
        "sampling are zero, which should not happen.");
  }
 
  for (int i = 0; i < monte_carlo_iterations; i++) {
    // Pick random pair of particles and resample their masses
    do {
      idx1 = random::uniform_int(0, n - 1);
      do {
        idx2 = random::uniform_int(0, n - 1);
      } while (idx2 == idx1);
    } while (types[idx1]->is_stable() && types[idx2]->is_stable());
 
    // Energy of the pair and CM frame velocity
    const FourVector p_sum = sampled_momenta[idx1] + sampled_momenta[idx2];
    const ThreeVector beta = p_sum.velocity();
    const double sqrts_12 = p_sum.abs();
    double m1, m2;
    do {
      m1 = types[idx1]->sample_full_spectral_function(sqrts_12);
      m2 = types[idx2]->sample_full_spectral_function(sqrts_12);
    } while (sqrts_12 < m1 + m2);
    const double pcm = pCM(sqrts_12, m1, m2);
    Angles phitheta;
    phitheta.distribute_isotropically();
    /*
     * Metropolis-Hastings acceptance. As the Breit Wigner is biased towards
     * lower masses compared to the full spectral function, we generally want to
     * accept more often if the ratio increases, and less if it decreases.
     */
    double acc = 1.0;
    if (!types[idx1]->is_stable()) {
      acc *= types[idx1]->ratio_spectral_full_to_breit_wigner(m1) /
             types[idx1]->ratio_spectral_full_to_breit_wigner(
                 sampled_momenta[idx1].abs());
    }
    if (!types[idx2]->is_stable()) {
      acc *= types[idx2]->ratio_spectral_full_to_breit_wigner(m2) /
             types[idx2]->ratio_spectral_full_to_breit_wigner(
                 sampled_momenta[idx2].abs());
    }
 
    if (random::canonical() < acc) {
      // Accept and impose momentum conservation with back-to-back particles
      sampled_momenta[idx1] =
          FourVector(std::sqrt(m1 * m1 + pcm * pcm), pcm * phitheta.threevec())
              .lorentz_boost(-beta);
      sampled_momenta[idx2] =
          FourVector(std::sqrt(m2 * m2 + pcm * pcm), -pcm * phitheta.threevec())
              .lorentz_boost(-beta);
    }
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

interpolate_trilinear()

double smash::detail::interpolate_trilinear	(	double	ax,
		double	ay,
		double	az,
		double	f1,
		double	f2,
		double	f3,
		double	f4,
		double	f5,
		double	f6,
		double	f7,
		double	f8
	)

Perform a trilinear 1st order interpolation.

Assume, we seek the value of a function \( f \) at position \((x, y, z)\). We know the position \((x, y, z)\) lies within a 3D cube, for which the values of the function f are known at each corner \((f_1, ..., f_8)\). We can now interpolate those values trilinearly to obtain an estimate of \( f \) at position \((x, y, z)\).

For this interpolation, linear functions are used in each direction \(x\), \(y\), and \(z\) respectively with \( a_x \), \( a_y \), and \( a_z \) as the slope parameters, e.g., \( f_1 + a_y \cdot (f_3 - f_1) \) for an approximation between the corners \( f_1 \) and \( f_3 \). For the \(y\)-direction, the linear functions are based on the cube's corners. The \(x\)-direction then uses the four obtained values from the interpolation in \(y\)-direction and finally the \(z\)-direction interpolation is based on the two values obtained by the combined interpolations in \(x\)- and \(y\)-direction. Since the position \((x, y, z)\) of the wanted value of function \( f \) is within the cube, the allowed values for \( a_x \), \( a_y \), and \( a_z \) are between 0 and 1.

Positional placement of the cube:

• \(x\)-direction: lower left front to lower right front corner \( (f_1 \) to \( f_2) \)
• \(y\)-direction: lower left front to upper left front corner \( (f_1 \) to \( f_3) \)
• \(z\)-direction: lower left front to lower left back corner \( (f_1 \) to \( f_5) \)

Note

\( a_x \), \( a_y \), and \( a_z \) have to be chosen in a way that the linear interpolations reflect the position \( (x, y, z) \), i.e. \( (x, y, z) = (x_1 + a_x \cdot x_2, y_1 + a_y \cdot y_3, z_1 + a_z \cdot z_5) \) with \( (x_i, y_i, z_i) \) representing the coordinates of the corner at \( f_i \).

Parameters

[in]	ax	Fraction of the step in x-direction and used as slope parameter
[in]	ay	Fraction of the step in y-direction and used as slope parameter
[in]	az	Fraction of the step in z-direction and used as slope parameter
[in]	f1	Value at the lower left front corner of the cube
[in]	f2	Value at the lower right front corner of the cube
[in]	f3	Value at the upper left front corner of the cube
[in]	f4	Value at the upper right front corner of the cube
[in]	f5	Value at the lower left back corner of the cube
[in]	f6	Value at the lower right back corner of the cube
[in]	f7	Value at the upper left back corner of the cube
[in]	f8	Value at the upper right back corner of the cube

Returns

Interpolated value

Definition at line 658 of file hadgas_eos.cc.

                                                              {
  assert(ax >= 0 && ax <= 1);
  assert(ay >= 0 && ay <= 1);
  assert(az >= 0 && az <= 1);
  double res = az * (ax * (ay * f8 + (1.0 - ay) * f6) +
                     (1.0 - ax) * (ay * f7 + (1.0 - ay) * f5)) +
               (1 - az) * (ax * (ay * f4 + (1.0 - ay) * f2) +
                           (1.0 - ax) * (ay * f3 + (1.0 - ay) * f1));
  return res;
}

Here is the caller graph for this function:

get_default_validator()

template<typename T >

const KeyTraits<T>::validator_type& smash::detail::get_default_validator ( )

noexcept

Function template to get a default trivial validator.

Returns

A const reference to a functor that always returns true .

Attention

It might look unnecessary to have a function returning the functor and you might think that the functor as a constant global variable template would be enough. However, this would be in general wrong because this functor is used in the Key constructors which are used by the InputKeys class, that is a collection of static Keys. Hence, since initialization order of static/global objects in C++ is undefined, we need to do something else. We use therefore the "construct on first use idiom", making the functor a static object in a function scope. For more information, refer for example to ISO C++ FAQ.

Definition at line 103 of file key.h.

                                                                            {
  static const typename KeyTraits<T>::validator_type always_true =
      [](const T&) noexcept { return true; };
  return always_true;
}

type_name()

template<typename T >

constexpr auto smash::detail::type_name ( )

constexpr

Get type of variable as string in a human-readable way.

Template Parameters

T	The type to be returned.

Returns

A std::string containing the name of the type.

Definition at line 27 of file numeric_cast.h.

                           {
  std::string_view name, prefix, suffix;
#ifdef __clang__
  name = __PRETTY_FUNCTION__;
  prefix = "auto smash::detail::type_name() [T = ";
  suffix = "]";
#elif defined(__GNUC__)
  name = __PRETTY_FUNCTION__;
  prefix = "constexpr auto smash::detail::type_name() [with T = ";
  suffix = "]";
#elif defined(_MSC_VER)
  name = __FUNCSIG__;
  prefix = "auto __cdecl smash::detail::type_name<";
  suffix = ">(void)";
#else
  name = "UNKNOWN";
  prefix = "";
  suffix = "";
#endif
  name.remove_prefix(prefix.size());
  name.remove_suffix(suffix.size());
  return std::string{name};
}

almost_equal_knuthish()

template<typename N , typename = std::enable_if_t<std::is_floating_point_v<N>>>

bool smash::detail::almost_equal_knuthish ( const N  x, const N  y, const N  epsilon, const N  threshold = N{0.0}  )

noexcept

Compare whether two floating-point numbers are approximately equal à la Knuth up to a given tolerance.

On top of Knuth's tolerance predicate some corner cases are treated and the caller can specify a threshold as last parameter to make the test consider numbers equal if the absolute value of their difference is below of it.

Parameters

[in]	x	First of the two numbers.
[in]	y	Second of the two numbers.
[in]	epsilon	The relative tolerance for the test.
[in]	threshold	Threshold for the number comparison. By default this is zero, implying no threshold is considered.

Returns

false if either x or y is not a finite number, provided that the type supports non-numeric representations (i.e. is infinite or NAN);

true if x == y;

true if x == 0 and if \( |x| \le \varepsilon\);

true if y == 0 and if \( |y| \le \varepsilon\);

true if \( |x - y| \le M_\mathrm{threshold} \);

true if \( |x - y| \le \varepsilon \cdot \max(|x|, |y|) \) (Knuth's tolerance predicate);

false otherwise.

Definition at line 59 of file numerics.h.

                                                {0.0}) noexcept {
  assert(epsilon > 0);
  assert(threshold >= 0);
  if constexpr (std::numeric_limits<N>::is_iec559) {
    if (!std::isfinite(x) || !std::isfinite(y)) {
      return false;
    }
  }
  if (x == y)
    return true;
  else if (x == 0)
    return std::abs(y) <= epsilon;
  else if (y == 0)
    return std::abs(x) <= epsilon;
  else
    return std::abs(x - y) <= threshold ||
           std::abs(x - y) <= epsilon * std::max(std::abs(x), std::abs(y));
}

write_in_chunk_impl()

template<typename Converter , class Range , std::enable_if_t< std::is_same_v< Range, Particles >||std::is_same_v< Range, ParticleList >, bool > = true>

void smash::detail::write_in_chunk_impl	(	const Range &	particles,
		const OutputFormatter< Converter > &	formatter,
		std::function< void(const typename Converter::type &)>	write,
		std::size_t	max_buffer_bytes = 1'000'000'000
	)

Writes particle data in multiple chunks if the total buffer size exceeds a predefined maximum.

This method avoids creating a single excessively large binary buffer when writing many particles at once. Instead, it splits the write into several smaller chunks. This can prevent excessive memory usage and improve stability on systems or filesystems that may have trouble with very large write calls.

The maximum buffer size is currently set to 1 GB (10^9 bytes), but this can be adapted in the future if needed. If the total data size of the particle block is below this threshold, the method simply delegates to the write function which should perform a single write call.

Otherwise, the data is accumulated particle by particle until the buffer reaches the threshold. The buffer is then given to the write function which should flush to disk.

Note

This utility does not strictly belong to this file. At the moment, the objects using it do not have a clean hierarchy that would allow both OscarOutput and BinaryOutput to inherit a shared implementation, hence it lives here temporarily.

Todo:

Once the hierarchy is cleaned up, move this into a common base class that OscarOutput and BinaryOutput inherit from.

Template Parameters

Converter	Converter used by OutputFormatter to produce the buffer type (must define Converter::type).
Range	Container type — enforced to be either Particles or ParticleList.

Parameters

[in]	particles	Container of particles whose particle_line representation is to be written.
[in]	formatter	Formatter responsible for converting particles into the corresponding Converter::type buffer representation.
[in]	write	Callable that receives each filled buffer chunk and performs the actual write to the underlying output (e.g. file).
[in]	max_buffer_bytes	Maximum buffer size in bytes before the accumulated data is flushed via write (default: 1'000'000'000).

Exceptions

std::runtime_error

If the estimated size of a single particle line exceeds half of max_buffer_bytes. In that case, only one particle would fit per chunk, which defeats the purpose of chunked writing, and the caller must increase max_buffer_bytes accordingly.

Definition at line 695 of file outputformatter.h.

                                                {
  if (particles.size() == 0)
    return;
 
  const std::size_t bytes_per_particle =
      formatter.compute_single_size(particles.front());
 
  if (2.0 * bytes_per_particle > max_buffer_bytes) {
    throw std::runtime_error(
        "write_in_chunk_impl: the estimated size of a single particle line "
        "exceeds half of the configured max_buffer_bytes.\n"
        "This effectively means only one particle would fit per chunk, "
        "which defeats the purpose of chunked writing.\n"
        "Increase max_buffer_bytes to at least twice the particle line size "
        "to use this function correctly.");
  }
 
  if (particles.size() * bytes_per_particle <= max_buffer_bytes) {
    write(formatter.particles_data_chunk(particles));
    return;
  }
 
  using Buffer = typename Converter::type;
  Buffer buffer;
  buffer.reserve(max_buffer_bytes);
  std::size_t current_size = 0;
 
  for (const auto& particle : particles) {
    Buffer line = formatter.single_particle_data(particle);
    const std::size_t line_size = line.size();
    if (current_size + line_size > max_buffer_bytes) {
      write(buffer);
      buffer.clear();
      current_size = 0;
    }
    buffer.insert(buffer.end(), std::make_move_iterator(line.begin()),
                  std::make_move_iterator(line.end()));
    current_size += line_size;
  }
 
  if (!buffer.empty()) {
    write(buffer);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function: