smash::GrandCanThermalizer Class Reference

#include <grandcan_thermalizer.h>

The GrandCanThermalizer class implements the following functionality:

1. Create a lattice and find the local rest frame energy density in each cell from the particles.
2. Remove particles from the cells, where the energy density is high enough. Save the energy, momentum and quantum numbers of the removed particles.
3. Sample new particles instead of the removed ones according to the grand-canonical thermal distribution, but with an additional constraint: the energy, momentum and quantum numbers should be the same as those of the removed particles.

The step 3. is a challenging task, so several algorithms are implemented that try to fulfil the requirements. The algorithms are a trade-off between mathematical rigour and computational speed. All of them are shown to reproduce the mean values of multiplicities correctly. However, this is not the case for multiplicity fluctuations. For details see Oliinychenko:2016vkg.

Definition at line 205 of file grandcan_thermalizer.h.

Collaboration diagram for smash::GrandCanThermalizer:

Public Member Functions
	GrandCanThermalizer (const std::array< double, 3 > lat_sizes, const std::array< int, 3 > n_cells, const std::array< double, 3 > origin, bool periodicity, double e_critical, double t_start, double delta_t, ThermalizationAlgorithm algo)
	Default constructor for the GranCanThermalizer to allocate the lattice. More...
	GrandCanThermalizer (Configuration &conf, const std::array< double, 3 > lat_sizes, const std::array< double, 3 > origin, bool periodicity)
bool	is_time_to_thermalize (const Clock &clock) const
	Check that the clock is close to n * period of thermalization, since the thermalization only happens at these times. More...
void	update_thermalizer_lattice (const Particles &particles, const DensityParameters &par, bool ignore_cells_under_threshold=true)
	Compute all the thermodynamical quantities on the lattice from particles. More...
ThreeVector	uniform_in_cell () const
void	renormalize_momenta (ParticleList &plist, const FourVector required_total_momentum)
	Changes energy and momenta of the particles in plist to match the required_total_momentum. More...
void	sample_multinomial (HadronClass particle_class, int N)
	The sample_multinomial function samples integer numbers n_i distributed according to the multinomial distribution with sum N: \( p(n_1, n_2, \dots) = \prod a_i^{n_i} \times \frac{N!}{n_1!n_2! \dots} \) if \( \sum n_i = N \) and \( p = 0 \) otherwise. More...
void	sample_in_random_cell_BF_algo (ParticleList &plist, const double time, size_t type_index)
	The total number of particles of species type_index is defined by mult_int_ array that is returned by. More...
void	thermalize_BF_algo (QuantumNumbers &conserved_initial, double time, int ntest)
	Samples particles according to the BF algorithm by making use of the. More...
template<typename F >
void	compute_N_in_cells_mode_algo (F &&condition)
	Computes average number of particles in each cell for the mode algorithm. More...
template<typename F >
ParticleData	sample_in_random_cell_mode_algo (const double time, F &&condition)
	Samples one particle and the species, cell, momentum and coordinate are chosen from the corresponding distributions. More...
void	thermalize_mode_algo (QuantumNumbers &conserved_initial, double time)
	Samples particles to the according to the mode algorithm. More...
void	thermalize (const Particles &particles, double time, int ntest)
	Main thermalize function, that chooses the algorithm to follow (BF or mode sampling). More...
void	print_statistics (const Clock &clock) const
	Generates standard output with information about the thermodynamic properties of the lattice, the thermalized region and the volume to be thermalized above the critical energy density. More...
RectangularLattice< ThermLatticeNode > &	lattice () const
	Getter function for the lattice. More...
double	e_crit () const
	Get the critical energy density. More...
ParticleList	particles_to_remove () const
	List of particles to be removed from the simulation. More...
ParticleList	particles_to_insert () const
	List of newly created particles to be inserted in the simulation. More...

Private Member Functions
ParticleTypePtrList	list_eos_particles () const
	Extracts the particles in the hadron gas equation of state from the complete list of particle types in SMASH. More...
HadronClass	get_class (size_t typelist_index) const
	Defines the class of hadrons by quantum numbers. More...
double	mult_class (const HadronClass cl) const

Private Attributes
std::vector< double >	N_in_cells_
	Number of particles to be sampled in one cell. More...
std::vector< size_t >	cells_to_sample_
	Cells above critical energy density. More...
HadronGasEos	eos_ = HadronGasEos(true)
	Hadron gas equation of state. More...
std::unique_ptr< RectangularLattice< ThermLatticeNode > >	lat_
	The lattice on which the thermodynamic quantities are calculated. More...
ParticleList	to_remove_
	Particles to be removed after this thermalization step. More...
ParticleList	sampled_list_
	Newly generated particles by thermalizer. More...
const ParticleTypePtrList	eos_typelist_
	List of particle types from which equation of state is computed. More...
const size_t	N_sorts_
	Number of different species to be sampled. More...
std::vector< double >	mult_sort_
	Real number multiplicity for each particle type. More...
std::vector< int >	mult_int_
	Integer multiplicity for each particle type. More...
std::array< double, 7 >	mult_classes_
	The different hadron species according to the enum defined in. More...
double	N_total_in_cells_
	Total number of particles over all cells in thermalization region. More...
double	cell_volume_
	Volume of a single cell, necessary to convert thermal densities to actual particle numbers. More...
const double	e_crit_
	Critical energy density above which cells are thermalized. More...
const double	t_start_
	Starting time of the simulation. More...
const double	period_
	Defines periodicity of the lattice in fm. More...
const ThermalizationAlgorithm	algorithm_
	Algorithm to choose for sampling of particles obeying conservation laws. More...

Constructor & Destructor Documentation

GrandCanThermalizer() [1/2]

smash::GrandCanThermalizer::GrandCanThermalizer	(	const std::array< double, 3 >	lat_sizes,
		const std::array< int, 3 >	n_cells,
		const std::array< double, 3 >	origin,
		bool	periodicity,
		double	e_critical,
		double	t_start,
		double	delta_t,
		ThermalizationAlgorithm	algo
	)

Default constructor for the GranCanThermalizer to allocate the lattice.

Parameters

[in]	lat_sizes	Size of lattice in x,y and z-direction in fm.
[in]	n_cells	Number of cells in x, y and z-direction.
[in]	origin	Coordinates of the left, down, near corner of the lattice in fm.
[in]	periodicity	Boolean to decide, if the lattice is periodically extended to infinity or not
[in]	e_critical	Critical energy density above which the cells are thermalized
[in]	t_start	Starting time of the simulation
[in]	delta_t	Timestep of the simulation
[in]	algo	Choice of algorithm for the canonical sampling

Definition at line 98 of file grandcan_thermalizer.cc.

    : eos_typelist_(list_eos_particles()),
      N_sorts_(eos_typelist_.size()),
      e_crit_(e_critical),
      t_start_(t_start),
      period_(delta_t),
      algorithm_(algo) {
  const LatticeUpdate upd = LatticeUpdate::EveryFixedInterval;
  lat_ = make_unique<RectangularLattice<ThermLatticeNode>>(
      lat_sizes, n_cells, origin, periodicity, upd);
  const std::array<double, 3> abc = lat_->cell_sizes();
  cell_volume_ = abc[0] * abc[1] * abc[2];
  cells_to_sample_.resize(50000);
  mult_sort_.resize(N_sorts_);
  mult_int_.resize(N_sorts_);
}

GrandCanThermalizer() [2/2]

smash::GrandCanThermalizer::GrandCanThermalizer ( Configuration &  conf, const std::array< double, 3 >  lat_sizes, const std::array< double, 3 >  origin, bool  periodicity  )

inline

See also

GrandCanThermalizer Exactly the same but taking values from config

Definition at line 227 of file grandcan_thermalizer.h.

      : GrandCanThermalizer(
            lat_sizes, conf.take({"Cell_Number"}), origin, periodicity,
            conf.take({"Critical_Edens"}), conf.take({"Start_Time"}),
            conf.take({"Timestep"}),
            conf.take({"Algorithm"}, ThermalizationAlgorithm::BiasedBF)) {}

Member Function Documentation

is_time_to_thermalize()

bool smash::GrandCanThermalizer::is_time_to_thermalize ( const Clock &  clock ) const

inline

Check that the clock is close to n * period of thermalization, since the thermalization only happens at these times.

Parameters

[in]

clock

Current system time

Definition at line 240 of file grandcan_thermalizer.h.

                                                       {
    const double t = clock.current_time();
    const int n = static_cast<int>(std::floor((t - t_start_) / period_));
    return (t > t_start_ &&
            t < t_start_ + n * period_ + clock.timestep_duration());
  }

Here is the call graph for this function:

update_thermalizer_lattice()

void smash::GrandCanThermalizer::update_thermalizer_lattice	(	const Particles &	particles,
		const DensityParameters &	par,
		bool	ignore_cells_under_threshold = true
	)

Compute all the thermodynamical quantities on the lattice from particles.

Parameters

[in]

particles

Current list of particles

See also

Particles

Parameters

[in]

par

Parameters necessary for density determination

See also

DensityParameters

Parameters

[in]

ignore_cells_under_threshold

Boolean that is true by default

Definition at line 120 of file grandcan_thermalizer.cc.

                                      {
  const DensityType dens_type = DensityType::Hadron;
  const LatticeUpdate update = LatticeUpdate::EveryFixedInterval;
  update_lattice(lat_.get(), update, dens_type, dens_par, particles);
  for (auto &node : *lat_) {
    /* If energy density is definitely below e_crit -
       no need to find T, mu, etc. So if e = T00 - T0i*vi <=
       T00 + sum abs(T0i) < e_crit, no efforts are necessary. */
    if (!ignore_cells_under_treshold ||
        node.Tmu0().x0() + std::abs(node.Tmu0().x1()) +
                std::abs(node.Tmu0().x2()) + std::abs(node.Tmu0().x3()) >=
            e_crit_) {
      node.compute_rest_frame_quantities(eos_);
    } else {
      node = ThermLatticeNode();
    }
  }
}

Here is the call graph for this function:

uniform_in_cell()

ThreeVector smash::GrandCanThermalizer::uniform_in_cell

(

)

const

Returns

3 vector uniformly sampled from the rectangular cell.

Definition at line 141 of file grandcan_thermalizer.cc.

                                                       {
  return ThreeVector(random::uniform(-0.5 * lat_->cell_sizes()[0],
                                     +0.5 * lat_->cell_sizes()[0]),
                     random::uniform(-0.5 * lat_->cell_sizes()[1],
                                     +0.5 * lat_->cell_sizes()[1]),
                     random::uniform(-0.5 * lat_->cell_sizes()[2],
                                     +0.5 * lat_->cell_sizes()[2]));
}

Here is the call graph for this function:

Here is the caller graph for this function:

renormalize_momenta()

void smash::GrandCanThermalizer::renormalize_momenta	(	ParticleList &	plist,
		const FourVector	required_total_momentum
	)

Changes energy and momenta of the particles in plist to match the required_total_momentum.

The procedure is described in Oliinychenko:2016vkg.

Parameters

[in]

plist

List of particles

See also

ParticleList

Parameters

[in]

required_total_momentum

The necessary total momentum of the cell

Definition at line 150 of file grandcan_thermalizer.cc.

                                                                   {
  const auto &log = logger<LogArea::GrandcanThermalizer>();

  // Centralize momenta
  QuantumNumbers conserved = QuantumNumbers(plist);
  log.info("Required 4-momentum: ", required_total_momentum);
  log.info("Sampled 4-momentum: ", conserved.momentum());
  const ThreeVector mom_to_add =
      (required_total_momentum.threevec() - conserved.momentum().threevec()) /
      plist.size();
  log.info("Adjusting momenta by ", mom_to_add);
  for (auto &particle : plist) {
    particle.set_4momentum(particle.type().mass(),
                           particle.momentum().threevec() + mom_to_add);
  }

  // Boost every particle to the common center of mass frame
  conserved = QuantumNumbers(plist);
  const ThreeVector beta_CM_generated = conserved.momentum().velocity();
  const ThreeVector beta_CM_required = required_total_momentum.velocity();

  double E = 0.0;
  double E_expected = required_total_momentum.abs();
  for (auto &particle : plist) {
    particle.boost_momentum(beta_CM_generated);
    E += particle.momentum().x0();
  }
  // Renorm. momenta by factor (1+a) to get the right energy, binary search
  const double tolerance = really_small;
  double a, a_min, a_max, er;
  const int max_iter = 50;
  int iter = 0;
  if (E_expected >= E) {
    a_min = 0.0;
    a_max = 1.0;
  } else {
    a_min = -1.0;
    a_max = 0.0;
  }
  do {
    a = 0.5 * (a_min + a_max);
    E = 0.0;
    for (const auto &particle : plist) {
      const double p2 = particle.momentum().threevec().sqr();
      const double E2 = particle.momentum().x0() * particle.momentum().x0();
      E += std::sqrt(E2 + a * (a + 2.0) * p2);
    }
    er = E - E_expected;
    if (er >= 0.0) {
      a_max = a;
    } else {
      a_min = a;
    }
    log.debug("Iteration ", iter, ": a = ", a, ", Δ = ", er);
    iter++;
  } while (std::abs(er) > tolerance && iter < max_iter);

  log.info("Renormalizing momenta by factor 1+a, a = ", a);
  for (auto &particle : plist) {
    particle.set_4momentum(particle.type().mass(),
                           (1 + a) * particle.momentum().threevec());
    particle.boost_momentum(-beta_CM_required);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample_multinomial()

void smash::GrandCanThermalizer::sample_multinomial	(	HadronClass	particle_class,
		int	N
	)

The sample_multinomial function samples integer numbers n_i distributed according to the multinomial distribution with sum N: \( p(n_1, n_2, \dots) = \prod a_i^{n_i} \times \frac{N!}{n_1!n_2! \dots} \) if \( \sum n_i = N \) and \( p = 0 \) otherwise.

Parameters

[in]

particle_class

A certain group of hadron species

See also

HadronClass

Parameters

[out]

N

Number of particles to be sampled

Todo:

(oliiny) what to do with this output?

Definition at line 216 of file grandcan_thermalizer.cc.

                                                              {
  /* The array mult_sort_ contains real numbers \f$ a_i \f$. The numbers \f$
   * n_i \f$ are saved in the mult_int_ array. Only particles of class
   * particle_class are sampled, where particle_class is defined by the
   * get_class function. */
  double sum = mult_class(particle_class);
  for (size_t i_type = 0; (i_type < N_sorts_) && (N_to_sample > 0); i_type++) {
    if (get_class(i_type) != particle_class) {
      continue;
    }
    const double p = mult_sort_[i_type] / sum;
    mult_int_[i_type] = random::binomial(N_to_sample, p);
    /*std::cout << eos_typelist_[i_type]->name() <<
             ": mult_sort = " << mult_sort_[i_type] <<
             ", sum = " << sum <<
             ", p = " << p <<
             ", N to sample = " << N_to_sample <<
             ", mult_int_ = " << mult_int_[i_type] << std::endl;*/
    sum -= mult_sort_[i_type];
    N_to_sample -= mult_int_[i_type];
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample_in_random_cell_BF_algo()

void smash::GrandCanThermalizer::sample_in_random_cell_BF_algo	(	ParticleList &	plist,
		const double	time,
		size_t	type_index
	)

The total number of particles of species type_index is defined by mult_int_ array that is returned by.

See also

sample_multinomial. This function samples mult_int_[type_index] particles. It chooses randomly the cell to sample and picks up momentum and coordinate from the corresponding distributions.

Parameters

[out]

plist

See also

ParticleList of newly produced particles

Parameters

[in]	time	Current time in the simulation to become zero component of sampled particles
[in]	type_index	Species that should be sampled

Definition at line 241 of file grandcan_thermalizer.cc.

                                                                           {
  N_in_cells_.clear();
  N_total_in_cells_ = 0.0;
  for (auto cell_index : cells_to_sample_) {
    const ThermLatticeNode cell = (*lat_)[cell_index];
    const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
    const double N_this_cell =
        cell_volume_ * gamma *
        HadronGasEos::partial_density(*eos_typelist_[type_index], cell.T(),
                                      cell.mub(), cell.mus());
    N_in_cells_.push_back(N_this_cell);
    N_total_in_cells_ += N_this_cell;
  }

  for (int i = 0; i < mult_int_[type_index]; i++) {
    // Choose random cell, probability = N_in_cell/N_total
    double r = random::uniform(0.0, N_total_in_cells_);
    double partial_sum = 0.0;
    int index_only_thermalized = -1;
    while (partial_sum < r) {
      index_only_thermalized++;
      partial_sum += N_in_cells_[index_only_thermalized];
    }
    const int cell_index = cells_to_sample_[index_only_thermalized];
    const ThermLatticeNode cell = (*lat_)[cell_index];
    const ThreeVector cell_center = lat_->cell_center(cell_index);

    ParticleData particle(*eos_typelist_[type_index]);
    // Note: it's pole mass for resonances!
    const double m = eos_typelist_[type_index]->mass();
    // Position
    particle.set_4position(FourVector(time, cell_center + uniform_in_cell()));
    // Momentum
    double momentum_radial = sample_momenta_from_thermal(cell.T(), m);
    Angles phitheta;
    phitheta.distribute_isotropically();
    particle.set_4momentum(m, phitheta.threevec() * momentum_radial);
    particle.boost_momentum(-cell.v());
    particle.set_formation_time(time);

    plist.push_back(particle);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

thermalize_BF_algo()

void smash::GrandCanThermalizer::thermalize_BF_algo	(	QuantumNumbers &	conserved_initial,
		double	time,
		int	ntest
	)

Samples particles according to the BF algorithm by making use of the.

See also

sample_in_random_cell_BF_algo. Quantum numbers of the sampled particles are required to be equal to the original particles in this region.

Parameters

[in]	conserved_initial	The quantum numbers of the total ensemble of of particles in the region to be thermalized
[in]	time	Current time of the simulation
[in]	ntest	Number of testparticles

Returns

Particle list with newly sampled particles according to Becattini-Feroni algorithm

Definition at line 287 of file grandcan_thermalizer.cc.

                                                                     {
  const auto &log = logger<LogArea::GrandcanThermalizer>();

  std::fill(mult_sort_.begin(), mult_sort_.end(), 0.0);
  for (auto cell_index : cells_to_sample_) {
    const ThermLatticeNode cell = (*lat_)[cell_index];
    const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
    for (size_t i = 0; i < N_sorts_; i++) {
      // N_i = n u^mu dsigma_mu = (isochronous hypersurface) n * V * gamma
      mult_sort_[i] += cell_volume_ * gamma * ntest *
                       HadronGasEos::partial_density(
                           *eos_typelist_[i], cell.T(), cell.mub(), cell.mus());
    }
  }

  std::fill(mult_classes_.begin(), mult_classes_.end(), 0.0);
  for (size_t i = 0; i < N_sorts_; i++) {
    mult_classes_[static_cast<size_t>(get_class(i))] += mult_sort_[i];
  }

  random::BesselSampler bessel_sampler_B(mult_class(HadronClass::Baryon),
                                         mult_class(HadronClass::Antibaryon),
                                         conserved_initial.baryon_number());

  while (true) {
    sampled_list_.clear();
    std::fill(mult_int_.begin(), mult_int_.end(), 0);
    const auto Nbar_antibar = bessel_sampler_B.sample();

    sample_multinomial(HadronClass::Baryon, Nbar_antibar.first);
    sample_multinomial(HadronClass::Antibaryon, Nbar_antibar.second);

    // Count strangeness of the sampled particles
    int S_sampled = 0;
    for (size_t i = 0; i < N_sorts_; i++) {
      S_sampled += eos_typelist_[i]->strangeness() * mult_int_[i];
    }

    std::pair<int, int> NS_antiS;
    if (algorithm_ == ThermalizationAlgorithm::BiasedBF) {
      random::BesselSampler bessel_sampler_S(
          mult_class(HadronClass::PositiveSMeson),
          mult_class(HadronClass::NegativeSMeson),
          conserved_initial.strangeness() - S_sampled);
      NS_antiS = bessel_sampler_S.sample();
    } else if (algorithm_ == ThermalizationAlgorithm::UnbiasedBF) {
      NS_antiS = std::make_pair(
          random::poisson(mult_class(HadronClass::PositiveSMeson)),
          random::poisson(mult_class(HadronClass::NegativeSMeson)));
      if (NS_antiS.first - NS_antiS.second !=
          conserved_initial.strangeness() - S_sampled) {
        continue;
      }
    }

    sample_multinomial(HadronClass::PositiveSMeson, NS_antiS.first);
    sample_multinomial(HadronClass::NegativeSMeson, NS_antiS.second);
    // Count charge of the sampled particles
    int ch_sampled = 0;
    for (size_t i = 0; i < N_sorts_; i++) {
      ch_sampled += eos_typelist_[i]->charge() * mult_int_[i];
    }

    std::pair<int, int> NC_antiC;
    if (algorithm_ == ThermalizationAlgorithm::BiasedBF) {
      random::BesselSampler bessel_sampler_C(
          mult_class(HadronClass::PositiveQZeroSMeson),
          mult_class(HadronClass::NegativeQZeroSMeson),
          conserved_initial.charge() - ch_sampled);
      NC_antiC = bessel_sampler_C.sample();
    } else if (algorithm_ == ThermalizationAlgorithm::UnbiasedBF) {
      NC_antiC = std::make_pair(
          random::poisson(mult_class(HadronClass::PositiveQZeroSMeson)),
          random::poisson(mult_class(HadronClass::NegativeQZeroSMeson)));
      if (NC_antiC.first - NC_antiC.second !=
          conserved_initial.charge() - ch_sampled) {
        continue;
      }
    }

    sample_multinomial(HadronClass::PositiveQZeroSMeson, NC_antiC.first);
    sample_multinomial(HadronClass::NegativeQZeroSMeson, NC_antiC.second);
    sample_multinomial(
        HadronClass::ZeroQZeroSMeson,
        random::poisson(mult_class(HadronClass::ZeroQZeroSMeson)));

    for (size_t itype = 0; itype < N_sorts_; itype++) {
      sample_in_random_cell_BF_algo(sampled_list_, time, itype);
    }
    double e_tot;
    const double e_init = conserved_initial.momentum().x0();
    e_tot = 0.0;
    for (auto &particle : sampled_list_) {
      e_tot += particle.momentum().x0();
    }
    if (std::abs(e_tot - e_init) > 0.01 * e_init) {
      log.debug("Rejecting: energy ", e_tot, " too far from e_init = ", e_init);
      continue;
    }
    break;
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

compute_N_in_cells_mode_algo()

template<typename F >

void smash::GrandCanThermalizer::compute_N_in_cells_mode_algo ( F &&  condition )

inline

Computes average number of particles in each cell for the mode algorithm.

Parameters

[in]

condition

Specifies the current mode (1 to 7)

Definition at line 314 of file grandcan_thermalizer.h.

                                                   {
    N_in_cells_.clear();
    N_total_in_cells_ = 0.0;
    for (auto cell_index : cells_to_sample_) {
      const ThermLatticeNode cell = (*lat_)[cell_index];
      const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
      double N_tot = 0.0;
      for (ParticleTypePtr i : eos_typelist_) {
        if (condition(i->strangeness(), i->baryon_number(), i->charge())) {
          // N_i = n u^mu dsigma_mu = (isochronous hypersurface) n * V * gamma
          N_tot += cell_volume_ * gamma *
                   HadronGasEos::partial_density(*i, cell.T(), cell.mub(),
                                                 cell.mus());
        }
      }
      N_in_cells_.push_back(N_tot);
      N_total_in_cells_ += N_tot;
    }
  }

Here is the call graph for this function:

Here is the caller graph for this function:

sample_in_random_cell_mode_algo()

template<typename F >

ParticleData smash::GrandCanThermalizer::sample_in_random_cell_mode_algo ( const double  time, F &&  condition  )

inline

Samples one particle and the species, cell, momentum and coordinate are chosen from the corresponding distributions.

The condition function limits the choice of possible species.

Condition is a function of the signature of quantum number S, B and Q. bool condition(int strangeness, int baryon_number, int charge);

Parameters

[in]	time	Current time in simulation
[in]	condition	Specifies the actual mode (1 to 7)

Definition at line 345 of file grandcan_thermalizer.h.

                                                              {
    // Choose random cell, probability = N_in_cell/N_total
    double r = random::uniform(0.0, N_total_in_cells_);
    double partial_sum = 0.0;
    int index_only_thermalized = -1;
    while (partial_sum < r) {
      index_only_thermalized++;
      partial_sum += N_in_cells_[index_only_thermalized];
    }
    const int cell_index = cells_to_sample_[index_only_thermalized];
    const ThermLatticeNode cell = (*lat_)[cell_index];
    const ThreeVector cell_center = lat_->cell_center(cell_index);
    const double gamma = 1.0 / std::sqrt(1.0 - cell.v().sqr());
    const double N_in_cell = N_in_cells_[index_only_thermalized];
    // Which sort to sample - probability N_i/N_tot
    r = random::uniform(0.0, N_in_cell);
    double N_sum = 0.0;
    ParticleTypePtr type_to_sample;
    for (ParticleTypePtr i : eos_typelist_) {
      if (!condition(i->strangeness(), i->baryon_number(), i->charge())) {
        continue;
      }
      N_sum +=
          cell_volume_ * gamma *
          HadronGasEos::partial_density(*i, cell.T(), cell.mub(), cell.mus());
      if (N_sum >= r) {
        type_to_sample = i;
        break;
      }
    }
    ParticleData particle(*type_to_sample);
    // Note: it's pole mass for resonances!
    const double m = type_to_sample->mass();
    // Position
    particle.set_4position(FourVector(time, cell_center + uniform_in_cell()));
    // Momentum
    double momentum_radial = sample_momenta_from_thermal(cell.T(), m);
    Angles phitheta;
    phitheta.distribute_isotropically();
    particle.set_4momentum(m, phitheta.threevec() * momentum_radial);
    particle.boost_momentum(-cell.v());
    particle.set_formation_time(time);
    return particle;
  }

Here is the call graph for this function:

Here is the caller graph for this function:

thermalize_mode_algo()

void smash::GrandCanThermalizer::thermalize_mode_algo	(	QuantumNumbers &	conserved_initial,
		double	time
	)

Samples particles to the according to the mode algorithm.

Quantum numbers of the sampled particles are required to be as in conserved_initial.

Parameters

[in]	conserved_initial	Quantum numbers of the original particles in the region to be thermalized
[in]	time	Current time of the simulation

Definition at line 391 of file grandcan_thermalizer.cc.

                                                    {
  double energy = 0.0;
  int S_plus = 0, S_minus = 0, B_plus = 0, B_minus = 0, E_plus = 0, E_minus = 0;
  // Mode 1: sample until energy is conserved, take only strangeness < 0
  auto condition1 = [](int, int, int) { return true; };
  compute_N_in_cells_mode_algo(condition1);
  while (conserved_initial.momentum().x0() > energy ||
         S_plus < conserved_initial.strangeness()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition1);
    energy += p.momentum().x0();
    if (p.pdgcode().strangeness() > 0) {
      sampled_list_.push_back(p);
      S_plus += p.pdgcode().strangeness();
    }
  }

  // Mode 2: sample until strangeness is conserved
  auto condition2 = [](int S, int, int) { return (S < 0); };
  compute_N_in_cells_mode_algo(condition2);
  while (S_plus + S_minus > conserved_initial.strangeness()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition2);
    const int s_part = p.pdgcode().strangeness();
    // Do not allow particles with S = -2 or -3 spoil the total sum
    if (S_plus + S_minus + s_part >= conserved_initial.strangeness()) {
      sampled_list_.push_back(p);
      S_minus += s_part;
    }
  }

  // Mode 3: sample non-strange baryons
  auto condition3 = [](int S, int, int) { return (S == 0); };
  QuantumNumbers conserved_remaining =
      conserved_initial - QuantumNumbers(sampled_list_);
  energy = 0.0;
  compute_N_in_cells_mode_algo(condition3);
  while (conserved_remaining.momentum().x0() > energy ||
         B_plus < conserved_remaining.baryon_number()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition3);
    energy += p.momentum().x0();
    if (p.pdgcode().baryon_number() > 0) {
      sampled_list_.push_back(p);
      B_plus += p.pdgcode().baryon_number();
    }
  }

  // Mode 4: sample non-strange anti-baryons
  auto condition4 = [](int S, int B, int) { return (S == 0) && (B < 0); };
  compute_N_in_cells_mode_algo(condition4);
  while (B_plus + B_minus > conserved_remaining.baryon_number()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition4);
    const int bar = p.pdgcode().baryon_number();
    if (B_plus + B_minus + bar >= conserved_remaining.baryon_number()) {
      sampled_list_.push_back(p);
      B_minus += bar;
    }
  }

  // Mode 5: sample non_strange mesons, but take only with charge > 0
  auto condition5 = [](int S, int B, int) { return (S == 0) && (B == 0); };
  conserved_remaining = conserved_initial - QuantumNumbers(sampled_list_);
  energy = 0.0;
  compute_N_in_cells_mode_algo(condition5);
  while (conserved_remaining.momentum().x0() > energy ||
         E_plus < conserved_remaining.charge()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition5);
    energy += p.momentum().x0();
    if (p.pdgcode().charge() > 0) {
      sampled_list_.push_back(p);
      E_plus += p.pdgcode().charge();
    }
  }

  // Mode 6: sample non_strange mesons to conserve charge
  auto condition6 = [](int S, int B, int C) {
    return (S == 0) && (B == 0) && (C < 0);
  };
  compute_N_in_cells_mode_algo(condition6);
  while (E_plus + E_minus > conserved_remaining.charge()) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition6);
    const int charge = p.pdgcode().charge();
    if (E_plus + E_minus + charge >= conserved_remaining.charge()) {
      sampled_list_.push_back(p);
      E_minus += charge;
    }
  }

  // Mode 7: sample neutral non-strange mesons to conserve energy
  auto condition7 = [](int S, int B, int C) {
    return (S == 0) && (B == 0) && (C == 0);
  };
  conserved_remaining = conserved_initial - QuantumNumbers(sampled_list_);
  energy = 0.0;
  compute_N_in_cells_mode_algo(condition7);
  while (conserved_remaining.momentum().x0() > energy) {
    ParticleData p = sample_in_random_cell_mode_algo(time, condition7);
    sampled_list_.push_back(p);
    energy += p.momentum().x0();
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

thermalize()

void smash::GrandCanThermalizer::thermalize	(	const Particles &	particles,
		double	time,
		int	ntest
	)

Main thermalize function, that chooses the algorithm to follow (BF or mode sampling).

Parameters

[out]	particles	List of sampled particles in thermalized region
[in]	time	Current time of the simulation
[in]	ntest	number of testparticles

Definition at line 492 of file grandcan_thermalizer.cc.

                                                {
  const auto &log = logger<LogArea::GrandcanThermalizer>();
  log.info("Starting forced thermalization, time ", time, " fm/c");
  to_remove_.clear();
  sampled_list_.clear();
  /* Remove particles from the cells with e > e_crit_,
   * sum up their conserved quantities */
  QuantumNumbers conserved_initial = QuantumNumbers();
  ThermLatticeNode node;
  for (auto &particle : particles) {
    const bool is_on_lattice =
        lat_->value_at(particle.position().threevec(), node);
    if (is_on_lattice && node.e() > e_crit_) {
      to_remove_.push_back(particle);
    }
  }
  /* Do not thermalize too small number of particles: for the number
   * of particles < 30 the algorithm tends to hang or crash too often. */
  if (to_remove_.size() > 30) {
    for (auto &particle : to_remove_) {
      conserved_initial.add_values(particle);
    }
  } else {
    to_remove_.clear();
    conserved_initial = QuantumNumbers();
  }
  log.info("Removed ", to_remove_.size(), " particles.");

  // Exit if there is nothing to thermalize
  if (conserved_initial == QuantumNumbers()) {
    return;
  }
  // Save the indices of cells inside the volume with e > e_crit_
  cells_to_sample_.clear();
  const size_t lattice_total_cells = lat_->size();
  for (size_t i = 0; i < lattice_total_cells; i++) {
    if ((*lat_)[i].e() > e_crit_) {
      cells_to_sample_.push_back(i);
    }
  }
  log.info("Number of cells in the thermalization region = ",
           cells_to_sample_.size(), ", its total volume [fm^3]: ",
           cells_to_sample_.size() * cell_volume_, ", in % of lattice: ",
           100.0 * cells_to_sample_.size() / lattice_total_cells);

  switch (algorithm_) {
    case ThermalizationAlgorithm::BiasedBF:
    case ThermalizationAlgorithm::UnbiasedBF:
      thermalize_BF_algo(conserved_initial, time, ntest);
      break;
    case ThermalizationAlgorithm::ModeSampling:
      thermalize_mode_algo(conserved_initial, time);
      break;
    default:
      throw std::invalid_argument(
          "This thermalization algorithm is"
          " not yet implemented");
  }
  log.info("Sampled ", sampled_list_.size(), " particles.");

  // Adjust momenta
  renormalize_momenta(sampled_list_, conserved_initial.momentum());
}

Here is the call graph for this function:

print_statistics()

void smash::GrandCanThermalizer::print_statistics

(

const Clock &

clock

)

const

Generates standard output with information about the thermodynamic properties of the lattice, the thermalized region and the volume to be thermalized above the critical energy density.

Parameters

[in]

clock

Current time of the simulation

Definition at line 557 of file grandcan_thermalizer.cc.

                                                                   {
  struct to_average {
    double T;
    double mub;
    double mus;
    double nb;
    double ns;
  };
  struct to_average on_lattice = {0.0, 0.0, 0.0, 0.0, 0.0};
  struct to_average in_therm_reg = {0.0, 0.0, 0.0, 0.0, 0.0};
  double e_sum_on_lattice = 0.0, e_sum_in_therm_reg = 0.0;
  int node_counter = 0;
  for (const auto &node : *lat_) {
    const double e = node.e();
    on_lattice.T += node.T() * e;
    on_lattice.mub += node.mub() * e;
    on_lattice.mus += node.mus() * e;
    on_lattice.nb += node.nb() * e;
    on_lattice.ns += node.ns() * e;
    e_sum_on_lattice += e;
    if (e >= e_crit_) {
      in_therm_reg.T += node.T() * e;
      in_therm_reg.mub += node.mub() * e;
      in_therm_reg.mus += node.mus() * e;
      in_therm_reg.nb += node.nb() * e;
      in_therm_reg.ns += node.ns() * e;
      e_sum_in_therm_reg += e;
      node_counter++;
    }
  }
  if (e_sum_on_lattice > really_small) {
    on_lattice.T /= e_sum_on_lattice;
    on_lattice.mub /= e_sum_on_lattice;
    on_lattice.mus /= e_sum_on_lattice;
    on_lattice.nb /= e_sum_on_lattice;
    on_lattice.ns /= e_sum_on_lattice;
  }
  if (e_sum_in_therm_reg > really_small) {
    in_therm_reg.T /= e_sum_in_therm_reg;
    in_therm_reg.mub /= e_sum_in_therm_reg;
    in_therm_reg.mus /= e_sum_in_therm_reg;
    in_therm_reg.nb /= e_sum_in_therm_reg;
    in_therm_reg.ns /= e_sum_in_therm_reg;
  }

  std::cout << "Current time [fm/c]: " << clock.current_time() << std::endl;
  std::cout << "Averages on the lattice - T[GeV], mub[GeV], mus[GeV], "
            << "nb[fm^-3], ns[fm^-3]: " << on_lattice.T << " " << on_lattice.mub
            << " " << on_lattice.mus << " " << on_lattice.nb << " "
            << on_lattice.ns << std::endl;
  std::cout << "Averages in therm. region - T[GeV], mub[GeV], mus[GeV], "
            << "nb[fm^-3], ns[fm^-3]: " << in_therm_reg.T << " "
            << in_therm_reg.mub << " " << in_therm_reg.mus << " "
            << in_therm_reg.nb << " " << in_therm_reg.ns << std::endl;
  std::cout << "Volume with e > e_crit [fm^3]: " << cell_volume_ * node_counter
            << std::endl;
}

Here is the call graph for this function:

lattice()

RectangularLattice<ThermLatticeNode>& smash::GrandCanThermalizer::lattice ( ) const

inline

Getter function for the lattice.

Definition at line 417 of file grandcan_thermalizer.h.

{ return *lat_; }

Here is the caller graph for this function:

e_crit()

double smash::GrandCanThermalizer::e_crit ( ) const

inline

Get the critical energy density.

Definition at line 419 of file grandcan_thermalizer.h.

{ return e_crit_; }

particles_to_remove()

ParticleList smash::GrandCanThermalizer::particles_to_remove ( ) const

inline

List of particles to be removed from the simulation.

Definition at line 421 of file grandcan_thermalizer.h.

{ return to_remove_; }

particles_to_insert()

ParticleList smash::GrandCanThermalizer::particles_to_insert ( ) const

inline

List of newly created particles to be inserted in the simulation.

Definition at line 423 of file grandcan_thermalizer.h.

{ return sampled_list_; }

list_eos_particles()

ParticleTypePtrList smash::GrandCanThermalizer::list_eos_particles ( ) const

inlineprivate

Extracts the particles in the hadron gas equation of state from the complete list of particle types in SMASH.

Definition at line 430 of file grandcan_thermalizer.h.

                                                 {
    ParticleTypePtrList res;
    for (const ParticleType& ptype : ParticleType::list_all()) {
      if (HadronGasEos::is_eos_particle(ptype)) {
        res.push_back(&ptype);
      }
    }
    return res;
  }

Here is the call graph for this function:

get_class()

HadronClass smash::GrandCanThermalizer::get_class ( size_t  typelist_index ) const

inlineprivate

Defines the class of hadrons by quantum numbers.

Parameters

[in]

typelist_index

Index for a certain quantum number

Definition at line 443 of file grandcan_thermalizer.h.

                                                     {
    const int B = eos_typelist_[typelist_index]->baryon_number();
    const int S = eos_typelist_[typelist_index]->strangeness();
    const int ch = eos_typelist_[typelist_index]->charge();
    // clang-format off
    return (B > 0) ? HadronClass::Baryon :
           (B < 0) ? HadronClass::Antibaryon :
           (S > 0) ? HadronClass::PositiveSMeson :
           (S < 0) ? HadronClass::NegativeSMeson :
           (ch > 0) ? HadronClass::PositiveQZeroSMeson :
           (ch < 0) ? HadronClass::NegativeQZeroSMeson :
                      HadronClass::ZeroQZeroSMeson;
    // clang-format on
  }

Here is the caller graph for this function:

mult_class()

double smash::GrandCanThermalizer::mult_class ( const HadronClass  cl ) const

inlineprivate

Parameters

[out]

cl

Multiplicity of the hadron class

Definition at line 458 of file grandcan_thermalizer.h.

                                                {
    return mult_classes_[static_cast<size_t>(cl)];
  }

Here is the caller graph for this function:

Member Data Documentation

N_in_cells_

std::vector<double> smash::GrandCanThermalizer::N_in_cells_

private

Number of particles to be sampled in one cell.

Definition at line 462 of file grandcan_thermalizer.h.

cells_to_sample_

std::vector<size_t> smash::GrandCanThermalizer::cells_to_sample_

private

Cells above critical energy density.

Definition at line 464 of file grandcan_thermalizer.h.

eos_

HadronGasEos smash::GrandCanThermalizer::eos_ = HadronGasEos(true)

private

Hadron gas equation of state.

Definition at line 466 of file grandcan_thermalizer.h.

lat_

std::unique_ptr<RectangularLattice<ThermLatticeNode> > smash::GrandCanThermalizer::lat_

private

The lattice on which the thermodynamic quantities are calculated.

Definition at line 468 of file grandcan_thermalizer.h.

to_remove_

ParticleList smash::GrandCanThermalizer::to_remove_

private

Particles to be removed after this thermalization step.

Definition at line 470 of file grandcan_thermalizer.h.

sampled_list_

ParticleList smash::GrandCanThermalizer::sampled_list_

private

Newly generated particles by thermalizer.

Definition at line 472 of file grandcan_thermalizer.h.

eos_typelist_

const ParticleTypePtrList smash::GrandCanThermalizer::eos_typelist_

private

List of particle types from which equation of state is computed.

Most particles are included, but not all of them. For example, photons and leptons are not included. Heavy hadrons, that can originate from Pythia, but do not interact in SMASH are not included. The latter are excluded to avoid violations of charm and bottomness conservation, when HadronGasEoS is used for forced thermalization.

Definition at line 481 of file grandcan_thermalizer.h.

N_sorts_

const size_t smash::GrandCanThermalizer::N_sorts_

private

Number of different species to be sampled.

Definition at line 483 of file grandcan_thermalizer.h.

mult_sort_

std::vector<double> smash::GrandCanThermalizer::mult_sort_

private

Real number multiplicity for each particle type.

Definition at line 485 of file grandcan_thermalizer.h.

mult_int_

std::vector<int> smash::GrandCanThermalizer::mult_int_

private

Integer multiplicity for each particle type.

Definition at line 487 of file grandcan_thermalizer.h.

mult_classes_

std::array<double, 7> smash::GrandCanThermalizer::mult_classes_

private

The different hadron species according to the enum defined in.

See also

HadronClass

Definition at line 492 of file grandcan_thermalizer.h.

N_total_in_cells_

double smash::GrandCanThermalizer::N_total_in_cells_

private

Total number of particles over all cells in thermalization region.

Definition at line 494 of file grandcan_thermalizer.h.

cell_volume_

double smash::GrandCanThermalizer::cell_volume_

private

Volume of a single cell, necessary to convert thermal densities to actual particle numbers.

Definition at line 499 of file grandcan_thermalizer.h.

e_crit_

const double smash::GrandCanThermalizer::e_crit_

private

Critical energy density above which cells are thermalized.

Definition at line 501 of file grandcan_thermalizer.h.

t_start_

const double smash::GrandCanThermalizer::t_start_

private

Starting time of the simulation.

Definition at line 503 of file grandcan_thermalizer.h.

period_

const double smash::GrandCanThermalizer::period_

private

Defines periodicity of the lattice in fm.

Definition at line 505 of file grandcan_thermalizer.h.

algorithm_

const ThermalizationAlgorithm smash::GrandCanThermalizer::algorithm_

private

Algorithm to choose for sampling of particles obeying conservation laws.

Definition at line 507 of file grandcan_thermalizer.h.

---

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

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