smash::QuantumSampling Class Reference

#include <quantumsampling.h>

This class:

• Calculates chemical potentials given density of particle species
• Calculates maxima of a Juttner distribution for these chemical potentials
• Samples Juttner distribution. This is the main intent of this class, while previous points are auxiliary calculations for it.

Definition at line 30 of file quantumsampling.h.

Collaboration diagram for smash::QuantumSampling:

Classes
struct	ParametersForMaximumMomentumRootFinder
	Struct object that holds the parameters relevant to finding the momentum for which the maximum of the distribution occurs. More...

Public Member Functions
	QuantumSampling (const std::map< PdgCode, int > &initial_multiplicities, double volume, double temperature)
	Constructor of a QuantumSampling object. More...
double	sample (const PdgCode pdg)
	Sampling radial momenta of given particle species from Boltzmann, Bose, or Fermi distribution. More...

Static Public Member Functions
static double	p_max_root_equation (double p, double mass, double temperature, double effective_chemical_potential, double statistics)
	Root equation for finding the radial momentum at which the Juttner distribution function has its maximum. More...
static int	p_max_root_equation_for_GSL (const gsl_vector *roots_array, void *parameters, gsl_vector *function)
	Root equation for finding the radial momentum at which the Juttner distribution function has its maximum, suited for the GSL root finding procedure. More...
static void	print_state_p_max (unsigned int iter, gsl_multiroot_fsolver *solver)
	A GSL utility which allows for printing out the status of the solver during the root finding procedure. More...
static int	find_p_at_maximum_of_the_distribution (double mass, double temperature, double effective_chemical_potential, double statistics, double p_max_initial_guess, double solution_precision, double *p_max)
	A GSL root solver for finding the radial momentum value at which the maximum of the given Juttner distribution function occurs. More...
static double	maximum_of_the_distribution (double mass, double temperature, double effective_chemical_potential, double statistics, double solution_precision)
	A convenience wrapper for finding the maximum value of the Juttner distribution, returning the value of the distribution for the momentum at which the maximum occurs, identified by find_p_at_maximum_of_the_distribution(). More...

Private Attributes
std::map< PdgCode, double >	effective_chemical_potentials_
	Tabulated effective chemical potentials for every particle species. More...
std::map< PdgCode, double >	distribution_function_maximums_
	Tabulated distribution function maxima for every particle species. More...
const double	volume_
	Volume [fm^3] in which particles sre sampled. More...
const double	temperature_
	Temperature [GeV]. More...

Constructor & Destructor Documentation

QuantumSampling()

smash::QuantumSampling::QuantumSampling	(	const std::map< PdgCode, int > &	initial_multiplicities,
		double	volume,
		double	temperature
	)

Constructor of a QuantumSampling object.

Parameters

[in]	initial_multiplicities	a map of pdg codes of samples particle species and corresponding multiplicities
[in]	volume	volume V in which the particles are sampled [fm^3], needed to calculate the density of the species
[in]	temperature	temperature T of the system [GeV]

Definition at line 210 of file quantumsampling.cc.

    : volume_(volume), temperature_(temperature) {
  /*
   * This is the precision which we expect from the solution; note that
   * solution precision also goes into the precision of calculating the
   * integrals involved etc. Recommended precision is at least 1e-7.
   */
  constexpr double solution_precision = 1e-8;
 
  for (const auto &pdg_and_mult : initial_multiplicities) {
    const PdgCode pdg = pdg_and_mult.first;
    const int number_of_particles = pdg_and_mult.second;
    const double V_in_GeV = volume_ / (hbarc * hbarc * hbarc);
    const double number_density = number_of_particles / V_in_GeV;
    const double spin_degeneracy = pdg.spin_degeneracy();
    // '+' for fermions, '-' for bosons
    const double quantum_statistics = (pdg.spin() % 2 == 0) ? -1.0 : 1.0;
    const ParticleType &ptype = ParticleType::find(pdg);
    const double particle_mass = ptype.mass();
    double chemical_potential = 0.0;
    ChemicalPotentialSolver mu_solver;
    // Calling the wrapper for the GSL chemical potential finder
    chemical_potential = mu_solver.effective_chemical_potential(
        spin_degeneracy, particle_mass, number_density, temperature_,
        quantum_statistics, solution_precision);
    effective_chemical_potentials_[pdg] = chemical_potential;
    const double distribution_function_maximum = maximum_of_the_distribution(
        particle_mass, temperature_, chemical_potential, quantum_statistics,
        solution_precision);
    distribution_function_maximums_[pdg] = distribution_function_maximum;
  }
}

Here is the call graph for this function:

Member Function Documentation

p_max_root_equation()

double smash::QuantumSampling::p_max_root_equation ( double  p, double  mass, double  temperature, double  effective_chemical_potential, double  statistics  )

static

Root equation for finding the radial momentum at which the Juttner distribution function has its maximum.

Parameters

[in]	p	radial momentum, i.e., length of the momentum vector [GeV]
[in]	mass	(pole) mass m of the particle species [GeV]
[in]	temperature	temperature T of the system [GeV]
[in]	effective_chemical_potential	effective chemical potential mu of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)

Returns

the extremum equation for the maximum of the Juttner distribution

Definition at line 33 of file quantumsampling.cc.

                                                               {
  const double Ekin = std::sqrt(p * p + mass * mass);
  const double term1 =
      2 * (1 + statistics * std::exp(-(Ekin - effective_chemical_potential) /
                                     temperature));
  const double term2 = (p * p) / (temperature * Ekin);
 
  return term1 - term2;
}

Here is the caller graph for this function:

p_max_root_equation_for_GSL()

int smash::QuantumSampling::p_max_root_equation_for_GSL ( const gsl_vector *  roots_array, void *  parameters, gsl_vector *  function  )

static

Root equation for finding the radial momentum at which the Juttner distribution function has its maximum, suited for the GSL root finding procedure.

Parameters

[in]	roots_array	an array holding the current best estimate of the roots of the solved equation
[in]	parameters	refers to the parameters as provided in the GSL root solving procedure
[in]	function	refers to the root equation(s) as provided in the GSL root solving procedure (where it's called "function")

Returns

the root equation suited for GSL root finding procedure

Definition at line 46 of file quantumsampling.cc.

                                                                       {
  struct ParametersForMaximumMomentumRootFinder *par =
      static_cast<struct ParametersForMaximumMomentumRootFinder *>(parameters);
 
  const double mass = (par->mass);
  const double temperature = (par->temperature);
  const double effective_chemical_potential =
      (par->effective_chemical_potential);
  const double statistics = (par->statistics);
 
  const double p_radial = gsl_vector_get(roots_array, 0);
 
  gsl_vector_set(function, 0,
                 p_max_root_equation(p_radial, mass, temperature,
                                     effective_chemical_potential, statistics));
 
  return GSL_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

print_state_p_max()

void smash::QuantumSampling::print_state_p_max ( unsigned int  iter, gsl_multiroot_fsolver *  solver  )

static

A GSL utility which allows for printing out the status of the solver during the root finding procedure.

Parameters

[in]	iter	variable keeping track of how many steps in the root solving procedure have been taken
[in]	solver	GSL solver object, which has acces to the current best estimate of the roots and the corresponding function values

Returns

message about the current state of the solver

Definition at line 67 of file quantumsampling.cc.

                                                                       {
  printf(
      "\n***\nfind_p_at_maximum_of_the_distribution(): iter = %3u \t"
      "x = % .3f \t"
      "f(x) = % .3e \n",
      iter, gsl_vector_get(solver->x, 0), gsl_vector_get(solver->f, 0));
}

find_p_at_maximum_of_the_distribution()

int smash::QuantumSampling::find_p_at_maximum_of_the_distribution ( double  mass, double  temperature, double  effective_chemical_potential, double  statistics, double  p_max_initial_guess, double  solution_precision, double *  p_max  )

static

A GSL root solver for finding the radial momentum value at which the maximum of the given Juttner distribution function occurs.

For the value of the distribution at the maximum, one shoud use the function maximum_of_the_distribution().

Parameters

[in]	mass	(pole) mass m of the particle species [GeV]
[in]	temperature	temperature T of the system [GeV]
[in]	effective_chemical_potential	effective chemical potential mu of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	p_max_initial_guess	the initial guess for the value of the solution [GeV]
[in]	solution_precision	precision with which the solution is found
[out]	p_max	the solution (momentum for which the distribution takes on the maximum value) stored in an array object [GeV]

Definition at line 76 of file quantumsampling.cc.

                   {
  const gsl_multiroot_fsolver_type *Solver_name;
  gsl_multiroot_fsolver *Root_finder;
 
  int status;
  size_t iter = 0;
  size_t initial_guess_update = 0;
 
  const size_t problem_dimension = 1;
 
  struct ParametersForMaximumMomentumRootFinder parameters = {
      mass, temperature, effective_chemical_potential, statistics};
 
  gsl_multiroot_function MaximumOfDistribution = {
      &p_max_root_equation_for_GSL, problem_dimension, &parameters};
 
  double roots_array_initial[1] = {p_max_initial_guess};
 
  gsl_vector *roots_array = gsl_vector_alloc(problem_dimension);
  gsl_vector_set(roots_array, 0, roots_array_initial[0]);
 
  Solver_name = gsl_multiroot_fsolver_hybrids;
  Root_finder = gsl_multiroot_fsolver_alloc(Solver_name, problem_dimension);
  gsl_multiroot_fsolver_set(Root_finder, &MaximumOfDistribution, roots_array);
 
  // print_state_p_max (iter, Root_finder);
 
  do {
    iter++;
 
    /*
     * gsl_multiroot_fsolver_iterate returns either 0 for a correct behavior,
     * or an error code (a positive integer) when the solver is stuck
     */
    status = gsl_multiroot_fsolver_iterate(Root_finder);
 
    // print_state_p_max (iter, Root_finder);
 
    /*
     * Check whether the solver is stuck
     */
    if (status) {
      if (initial_guess_update < 100) {
        /*
         * In case the solution is not found for the (somewhat small) default
         * initial guess of p_max_initial_guess = 0.05 [GeV], the value of the
         * p_max_initial_guess is increased and the root solving procedure is
         * restarted. This can take place up to a 100 times, with the largest
         * possible initial guess of p_max_initial_guess = 5.05 [GeV].
         */
        p_max_initial_guess += 0.05;
        initial_guess_update++;
        roots_array_initial[0] = p_max_initial_guess;
        gsl_vector_set(roots_array, 0, roots_array_initial[0]);
        gsl_multiroot_fsolver_set(Root_finder, &MaximumOfDistribution,
                                  roots_array);
        iter = 0;
      } else {
        std::cout << "\n\nGSL error message:\n"
                  << gsl_strerror(status) << std::endl;
        logg[LogArea::Distributions::id].warn(
            "\n\nThe GSL solver"
            "\nfind_p_at_maximum_of_the_distribution\nis stuck!"
            "\n\nInput parameters:"
            "\n                  mass [GeV] = ",
            mass, "\n           temperature [GeV] = ", temperature,
            "\neffective_chemical_potential = ", effective_chemical_potential,
            "\n                  statistics = ", statistics,
            "\n          solution_precision = ", solution_precision,
            "\n\n"
            "Initialization cannot sample the momenta without "
            "calculating the distribution maximum."
            "\nTry adjusting the initial guess (which is "
            "looped over in the GSL procedure) or the "
            "solution precision."
            "\nUncomment print_state_p_max to check solver progress.\n\n\n");
        throw std::runtime_error(
            "QuantumSampling::find_p_at_maximum_of_the_distribution returned "
            "no result.\n\n");
        continue;
      }
    }
 
    status = gsl_multiroot_test_residual(Root_finder->f, solution_precision);
 
    if (status == GSL_SUCCESS) {
      p_max[0] = gsl_vector_get(Root_finder->x, 0);
    }
  } while (status == GSL_CONTINUE && iter < 100000);
 
  gsl_multiroot_fsolver_free(Root_finder);
  gsl_vector_free(roots_array);
 
  return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function:

maximum_of_the_distribution()

double smash::QuantumSampling::maximum_of_the_distribution ( double  mass, double  temperature, double  effective_chemical_potential, double  statistics, double  solution_precision  )

static

A convenience wrapper for finding the maximum value of the Juttner distribution, returning the value of the distribution for the momentum at which the maximum occurs, identified by find_p_at_maximum_of_the_distribution().

Parameters

[in]	mass	(pole) mass m of the particle species [GeV]
[in]	temperature	temperature T of the system [GeV]
[in]	effective_chemical_potential	effective chemical potential mu of the system [GeV]
[in]	statistics	quantum statistics of the particles species (+1 for Fermi, -1 for Bose, 0 for Boltzmann)
[in]	solution_precision	precision with which the solution is found

Definition at line 175 of file quantumsampling.cc.

                                                  {
  /*
   * Momentum at which the distribution function has its maximum.
   */
  double p_max[1];
  p_max[0] = 0.0;
 
  /*
   * Initial guess for the value of p_max, in GeV. This value is
   * looped over within the GSL solver, so that many initial guesses
   * are probed if the solution is not found.
   */
  double initial_guess_p_max = 0.050;
 
  /*
   * Calling the GSL distribution maximum finder
   */
  find_p_at_maximum_of_the_distribution(
      mass, temperature, effective_chemical_potential, statistics,
      initial_guess_p_max, solution_precision, p_max);
 
  double distribution_function_maximum =
      p_max[0] * p_max[0] *
      juttner_distribution_func(p_max[0], mass, temperature,
                                effective_chemical_potential, statistics);
 
  return distribution_function_maximum;
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample()

double smash::QuantumSampling::sample

(

const PdgCode

pdg

)

Sampling radial momenta of given particle species from Boltzmann, Bose, or Fermi distribution.

This sampler uses the simplest rejection sampling.

Parameters

[in]

pdg

the pdg code of the sampled particle species return the sampled momentum [GeV]

Definition at line 256 of file quantumsampling.cc.

                                                {
  const ParticleType &ptype = ParticleType::find(pdg);
  const double mass = ptype.mass();
  const double mu = effective_chemical_potentials_.find(pdg)->second;
  const double distr_max = distribution_function_maximums_.find(pdg)->second;
  /*
   * The variable maximum_momentum denotes the "far right" boundary of the
   * sampled region; we assume that no particle has momentum larger than 10 GeV
   */
  constexpr double maximum_momentum = 10.0;  // in [GeV]
  const double statistics = (pdg.spin() % 2 == 0) ? -1.0 : 1.0;
  double sampled_momentum = 0.0, sampled_ratio = 0.0;
 
  do {
    sampled_momentum = random::uniform(0.0, maximum_momentum);
    double distribution_at_sampled_p =
        sampled_momentum * sampled_momentum *
        juttner_distribution_func(sampled_momentum, mass, temperature_, mu,
                                  statistics);
    sampled_ratio = distribution_at_sampled_p / distr_max;
  } while (random::canonical() > sampled_ratio);
 
  return sampled_momentum;
}

Here is the call graph for this function:

Member Data Documentation

effective_chemical_potentials_

std::map<PdgCode, double> smash::QuantumSampling::effective_chemical_potentials_

private

Tabulated effective chemical potentials for every particle species.

Definition at line 155 of file quantumsampling.h.

distribution_function_maximums_

std::map<PdgCode, double> smash::QuantumSampling::distribution_function_maximums_

private

Tabulated distribution function maxima for every particle species.

Definition at line 157 of file quantumsampling.h.

volume_

const double smash::QuantumSampling::volume_

private

Volume [fm^3] in which particles sre sampled.

Definition at line 159 of file quantumsampling.h.

temperature_

const double smash::QuantumSampling::temperature_

private

Temperature [GeV].

Definition at line 161 of file quantumsampling.h.

---

The documentation for this class was generated from the following files:

• src/include/smash/quantumsampling.h
• src/quantumsampling.cc