smash::Nucleus Class Reference

#include <nucleus.h>

A nucleus is a collection of particles that are initialized, before the beginning of the simulation and all have the same velocity.

Definition at line 27 of file nucleus.h.

Inheritance diagram for smash::Nucleus:

Classes
struct	TestparticleConfusion

Public Member Functions
	Nucleus ()=default
	default constructor More...
	Nucleus (Configuration &config, int nTest)
	Constructor for Nucleus, that needs the configuration parameters from the inputfile and the number of testparticles. More...
	Nucleus (const std::map< PdgCode, int > &particle_list, int nTest)
	Constructor which directly initializes the Nucleus with particles and respective counts. More...
virtual	~Nucleus ()=default
double	mass () const
virtual ThreeVector	distribute_nucleon ()
	The distribution of return values from this function is according to a spherically symmetric Woods-Saxon distribution suitable for this nucleus. More...
double	woods_saxon (double x)
	Woods-Saxon distribution. More...
virtual void	arrange_nucleons ()
	Sets the positions of the nucleons inside a nucleus. More...
virtual void	set_parameters_automatic ()
	Sets the deformation parameters of the Woods-Saxon distribution according to the current mass number. More...
virtual void	set_parameters_from_config (Configuration &config)
	Sets the parameters of the Woods-Saxon according to manually added values in the configuration file. More...
virtual void	generate_fermi_momenta ()
	Generates momenta according to Fermi motion for the nucleons. More...
void	boost (double beta_scalar)
	Boosts the nuclei into the computational frame, such that the nucleons have the appropriate momentum and the nuclei are lorentz-contracted. More...
void	fill_from_list (const std::map< PdgCode, int > &particle_list, int testparticles)
	Adds particles from a map PDG code => Number_of_particles_with_that_PDG_code to the nucleus. More...
void	shift (double z_offset, double x_offset, double simulation_time)
	Shifts the nucleus to correct impact parameter and z displacement. More...
virtual void	rotate ()
	Rotates the nucleus. More...
void	copy_particles (Particles *particles)
	Copies the particles from this nucleus into the particle list. More...
size_t	size () const
	Number of numerical (=test-)particles in the nucleus: More...
size_t	number_of_particles () const
	Number of physical particles in the nucleus: More...
size_t	number_of_protons () const
	Number of physical protons in the nucleus: More...
FourVector	center () const
	Calculate geometrical center of the nucleus. More...
void	set_label (BelongsTo label)
	Sets target / projectile labels on nucleons. More...
void	align_center ()
	Shifts the nucleus so that its center is at (0,0,0) More...
virtual double	nucleon_density (double r, double, double) const
	Return the Woods-Saxon probability density for the given position. More...
virtual double	nucleon_density_unnormalized (double r, double, double) const
	Return the unnormalized Woods-Saxon distribution for the given position without deformation. More...
virtual double	calculate_saturation_density () const
virtual void	set_saturation_density (double density)
	Sets the saturation density of the nucleus. More...
std::vector< ParticleData >::iterator	begin ()
	For iterators over the particle list: More...
std::vector< ParticleData >::iterator	end ()
	For iterators over the particle list: More...
std::vector< ParticleData >::const_iterator	cbegin () const
	For const iterators over the particle list: More...
std::vector< ParticleData >::const_iterator	cend () const
	For const iterators over the particle list: More...
void	set_diffusiveness (double diffuse)
	Sets the diffusiveness of the nucleus. More...
double	get_diffusiveness () const
double	get_saturation_density () const
double	default_nuclear_radius ()
	Default nuclear radius calculated as: More...
void	set_nuclear_radius (double rad)
	Sets the nuclear radius. More...
double	get_nuclear_radius () const

Protected Member Functions
void	random_euler_angles ()
	Randomly generate Euler angles. More...

Protected Attributes
std::vector< ParticleData >	particles_
	Particles associated with this nucleus. More...
double	saturation_density_ = nuclear_density
	Saturation density of this nucleus. More...
double	euler_phi_
	Euler angel phi. More...
double	euler_theta_
	Euler angel theta. More...
double	euler_psi_
	Euler angel psi. More...

Private Attributes
double	diffusiveness_
	Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0, we obtain a hard sphere). More...
double	nuclear_radius_
	Nuclear radius of this nucleus. More...
double	proton_radius_ = 1.2
	Single proton radius in fm. More...
size_t	testparticles_ = 1
	Number of testparticles per physical particle. More...

Friends
std::ostream &	operator<< (std::ostream &, const Nucleus &)
	Writes the state of the Nucleus object to the output stream. More...

Constructor & Destructor Documentation

Nucleus() [1/3]

smash::Nucleus::Nucleus ( )

default

default constructor

Nucleus() [2/3]

smash::Nucleus::Nucleus	(	Configuration &	config,
		int	nTest
	)

Constructor for Nucleus, that needs the configuration parameters from the inputfile and the number of testparticles.

Parameters

[in]	config	contains the parameters from the inputfile on the numbers of particles with a certain PDG code
[in]	nTest	number of testparticles

Definition at line 31 of file nucleus.cc.

                                                 {
  // Fill nuclei with particles.
  std::map<PdgCode, int> part = config.take({"Particles"});
  fill_from_list(part, nTest);
  // Look for user-defined values or take the default parameters.
  if (config.has_value({"Diffusiveness"}) && config.has_value({"Radius"}) &&
      config.has_value({"Saturation_Density"})) {
    set_parameters_from_config(config);
  } else if (!config.has_value({"Diffusiveness"}) &&
             !config.has_value({"Radius"}) &&
             !config.has_value({"Saturation_Density"})) {
    set_parameters_automatic();
    set_saturation_density(calculate_saturation_density());
  } else {
    throw std::invalid_argument(
        "Diffusiveness, Radius and Saturation_Density "
        "required to manually configure the Woods-Saxon"
        " distribution. Only one/two were provided. \n"
        "Providing none of the above mentioned "
        "parameters automatically configures the "
        "distribution based on the atomic number.");
  }
}

Nucleus() [3/3]

smash::Nucleus::Nucleus	(	const std::map< PdgCode, int > &	particle_list,
		int	nTest
	)

Constructor which directly initializes the Nucleus with particles and respective counts.

Only used for testing.

Parameters

[in]	particle_list	std::map, which maps PdgCode and count of this particle.
[in]	nTest	Number of test particles.

Definition at line 25 of file nucleus.cc.

                                                                     {
  fill_from_list(particle_list, nTest);
  set_parameters_automatic();
  set_saturation_density(calculate_saturation_density());
}

~Nucleus()

virtual smash::Nucleus::~Nucleus ( )

virtualdefault

Member Function Documentation

mass()

double smash::Nucleus::mass

(

)

const

Returns

Mass of the nucleus [GeV]. It needs to be double to allow for calculations at LHC energies.

Definition at line 55 of file nucleus.cc.

                           {
  double total_mass = 0.;
  for (auto i = cbegin(); i != cend(); i++) {
    total_mass += i->momentum().abs();
  }
  return total_mass / (testparticles_ + 0.0);
}

distribute_nucleon()

ThreeVector smash::Nucleus::distribute_nucleon ( )

virtual

The distribution of return values from this function is according to a spherically symmetric Woods-Saxon distribution suitable for this nucleus.

\(\frac{dN}{dr} = \frac{r^2}{\exp\left(\frac{r-r_0}{d}\right) + 1}\) where \(d\) is the diffusiveness_ parameter and \(r_0\) is nuclear_radius_.

Returns

Woods-Saxon distributed position.

Woods-Saxon-distribution

The distribution

Nucleons in nuclei are distributed according to a Woods-Saxon-distribution (see Woods:1954zz [62])

\[\frac{dN}{d^3r} = \frac{\rho_0}{\exp\left(\frac{r-r_0}{d}\right) +1},\]

where \(d\) is the diffusiveness of the nucleus. For \(d=0\), the nucleus is a hard sphere. \(\rho_0\) and \(r_0\) are, in this limit, the nuclear ground state density and nuclear radius, respectively. For small \(d\), this is still approximately true.

This distribution is obviously spherically symmetric, hence we can rewrite \(d^3r = 4\pi r^2 dr\) and obtain

\[\frac{dN}{4\pi\rho_0dr} = \frac{r^2}{\exp\left(\frac{r-r_0}{d}\right) + 1}.\]

Let us rewrite that in units of \(d\) (that's the diffusiveness) and drop any constraints on normalization (since in the end we only care about relative probabilities: we create as many nucleons as we need). Now, \(p(B)\) is the un-normalized probability to obtain a point at \(r = Bd\) (with \(R = r_0/d\)):

\[p(B) = \frac{B^2}{\exp(B-R) + 1}.\]

Splitting it up in two regimes

We shift the distribution so that \(B-R\) is 0 at \(t = 0\): \(t = B-R\):

\[p^{(1)}(t)= \frac{(t+R)^2}{\exp(t)+1}\]

and observe

\[\frac{1}{\exp(x)+1} = \frac{e^{-x}}{e^{-x}e^{x}+e^{-x}} = \frac{e^{-x}}{e^{-x}+1}.\]

The distribution function can now be split into two cases. For negative t (first case), \(-|t| = t\), and for positive t (second case), \(-|t| = -t\):

\[p^{(1)}(t) = \frac{1}{e^{-|t|}+1} \cdot (t+R)^2 \cdot \begin{cases} 1 & -R \le t < 0 \\ e^{-t} & t \ge 0 \end{cases}.\]

Apart from the first term, all that remains here can easily and exactly be generated from unrejected uniform random numbers (see below). The first term itself - \((1+e^{-|t|})^{-1}\) - is a number between 1/2 and 1.

If we now have a variable \(t\) distributed according to the remainder, \(p^{(2)}(t)\), and reject \(t\) with a probability \(p^{(rej)}(t) = 1 - p^{(survive)}(t) = 1 - (1+e^{-|t|})^{-1}\), the resulting distribution is \(p^{(combined)}(t) = p^{(2)}(t) \cdot p^{(survive)}(t)\). Hence, we need to generate \(p^{(2)}(t)\), which we can normalize to

\[\tilde{p}^{(2)}(t) = \frac{1}{1+3/R+6/R^2+6/R^3} \cdot \begin{cases} \frac{3}{R^3} (t+R)^2 & -R \le t < 0 \\ e^{-t} \left( \frac{3}{R}+\frac{6}{R^2}t+\frac{6}{R^3}\frac{1}{2}t^2 \right) & t \ge 0 \end{cases}.\]

(the tilde \(\tilde{p}\) means that this is normalized).

Four parts inside the rejection

Let \(c_1 = 1+3/R+6/R^2+6/R^3\). The above means:

\[\mbox{Choose: } \begin{cases} \tilde p^{({\rm I})} = \frac{3}{R^3}(t+R)^2 \Theta(-t) \Theta(t+R) \\ \tilde p^{({\rm II})}= e^{-t}\Theta(t) \\ \tilde p^{({\rm III})}=e^{-t}\Theta(t) t \\ \tilde p^{({\rm IV})} =e^{-t}\Theta(t) \frac{1}{2} t^2 \end{cases} \mbox{ with a probability of }\begin{cases} \frac{1}{c_1} \cdot 1 \\ \frac{1}{c_1} \cdot \frac{3}{R} \\ \frac{1}{c_1} \cdot \frac{6}{R^2} \\ \frac{1}{c_1} \cdot \frac{6}{R^3} \end{cases}.\]

Let us see how those are generated. \(\chi_i\) are uniformly distributed numbers between 0 and 1.

\[p(\chi_i) = \Theta(\chi_i)\Theta(1-\chi_i)\]

For simple distributions (only one \(\chi\) involved), we invert \(t(\chi)\), derive it w.r.t. \(t\) and normalize.

Case I: \f$p^{({\rm I})}\f$

Simply from one random number:

\[t = R\left( \sqrt[ 3 ]{\chi} - 1 \right)\]

\[\tilde p^{({\rm I})} = \frac{3}{R^3}(t+R)^2 \mbox{ for } -R \le t \le 0\]

Case II: \f$p^{({\rm II})}\f$

Again, from one only:

\[t = -\log(\chi)\]

\[p(t) = \frac{d\chi}{dt}\]

\[p^{({\rm II})} = e^{-t} \mbox{ for } t > 0\]

Case III: \f$p^{({\rm III})}\f$

Here, we need two variables:

\[t = -\log{\chi_1} -\log{\chi_2}\]

\(p^{({\rm III})}\) is now the folding of \(p^{({\rm II})}\) with itself[1]:

\[p^{({\rm III})} = \int_{-\infty}^{\infty} d\tau e^{-\tau} e^{-(t-\tau)} \Theta(\tau) \Theta(t-\tau) = t e^{-t} \mbox{ for } t > 0\]

Case IV: \f$p^{({\rm IV})}\f$

Three variables needed:

\[t = -\log{\chi_1} -\log{\chi_2} -\log{\chi_3}\]

\(p^{({\rm IV})}\) is now the folding of \(p^{({\rm II})}\) with \(p^{({\rm III})}\):

\[p^{({\rm IV})} = \int_{ - \infty}^{\infty} d\tau e^{- \tau} \left( t - \tau \right) e^{ - (t - \tau)} \Theta(\tau) \Theta(t - \tau) = \frac{1}{2} t^2 e^{ -t} \mbox{ for } t > 0\]

[1]: This is [the probability to find a \(\tau\)] times [the probability to find the value \(\tau_2 = t-\tau\) that added to \(\tau\) yields \(t\)], integrated over all possible combinations that have that property.

From the beginning

So, the algorithm needs to do all this from the end:

• Decide which branch \(\tilde p^{({\rm I - IV})}\) to go into

• Generate \(t\) from the distribution in the respective branches
• Reject that number with a probability \(1-(1+\exp(-|t|))^{-1}\) (the efficiency of this should be \(\gg \frac{1}{2}\))

• Shift and rescale \(t\) to \(r = d\cdot t + r_0\)

Reimplemented in smash::DeformedNucleus, and smash::CustomNucleus.

Definition at line 215 of file nucleus.cc.

                                        {
  // Get the solid angle of the nucleon.
  Angles dir;
  dir.distribute_isotropically();
  // diffusiveness_ zero or negative? Use hard sphere.
  if (almost_equal(diffusiveness_, 0.)) {
    return dir.threevec() * nuclear_radius_ * std::cbrt(random::canonical());
  }
  if (almost_equal(nuclear_radius_, 0.)) {
    return smash::ThreeVector();
  }
  double radius_scaled = nuclear_radius_ / diffusiveness_;
  double prob_range1 = 1.0;
  double prob_range2 = 3. / radius_scaled;
  double prob_range3 = 2. * prob_range2 / radius_scaled;
  double prob_range4 = 1. * prob_range3 / radius_scaled;
  double ranges234 = prob_range2 + prob_range3 + prob_range4;
  double t;
  do {
    double which_range = random::uniform(-prob_range1, ranges234);
    if (which_range < 0.0) {
      t = radius_scaled * (std::cbrt(random::canonical()) - 1.);
    } else {
      t = -std::log(random::canonical());
      if (which_range >= prob_range2) {
        t -= std::log(random::canonical());
        if (which_range >= prob_range2 + prob_range3) {
          t -= std::log(random::canonical());
        }
      }
    }
  } while (random::canonical() > 1. / (1. + std::exp(-std::abs(t))));
  double position_scaled = t + radius_scaled;
  double position = position_scaled * diffusiveness_;
  return dir.threevec() * position;
}

woods_saxon()

double smash::Nucleus::woods_saxon

(

double

x

)

Woods-Saxon distribution.

Parameters

[in]

x

the position at which to evaluate the function

Returns

un-normalized Woods-saxon probability

Definition at line 261 of file nucleus.cc.

                                    {
  return r * r / (std::exp((r - nuclear_radius_) / diffusiveness_) + 1);
}

arrange_nucleons()

void smash::Nucleus::arrange_nucleons ( )

virtual

Sets the positions of the nucleons inside a nucleus.

Reimplemented in smash::CustomNucleus.

Definition at line 265 of file nucleus.cc.

                               {
  for (auto i = begin(); i != end(); i++) {
    // Initialize momentum
    i->set_4momentum(i->pole_mass(), 0.0, 0.0, 0.0);
    /* Sampling the Woods-Saxon, get the radial
     * position and solid angle for the nucleon. */
    ThreeVector pos = distribute_nucleon();
 
    // Set the position of the nucleon.
    i->set_4position(FourVector(0.0, pos));
  }
 
  // Recenter and rotate
  align_center();
  rotate();
}

set_parameters_automatic()

void smash::Nucleus::set_parameters_automatic ( )

virtual

Sets the deformation parameters of the Woods-Saxon distribution according to the current mass number.

The values are taken from DeVries:1987atn [19] and Loizides:2014vua [35]. They are in agreement with MC-Glauber models such as GLISSANDO (see Rybczynski:2013yba [46]) and TGlauber MC (see Loizides:2017ack [36]).

Definition at line 282 of file nucleus.cc.

                                       {
  int A = Nucleus::number_of_particles();
  int Z = Nucleus::number_of_protons();
  if (A == 1) {  // single particle
    /* In case of testparticles, an infinite reaction loop will be
     * avoided by a small finite spread according to a single particles
     * 'nucleus'. The proper solution will be to introduce parallel
     * ensembles. */
    set_nuclear_radius(
        testparticles_ == 1 ? 0. : 1. - std::exp(-(testparticles_ - 1.) * 0.1));
    set_diffusiveness(testparticles_ == 1 ? -1. : 0.02);
  } else if ((A == 238) && (Z == 92)) {  // Uranium
    // Default values.
    set_diffusiveness(0.556);
    set_nuclear_radius(6.86);
  } else if ((A == 208) && (Z == 82)) {  // Lead
    // Default values.
    set_diffusiveness(0.54);
    set_nuclear_radius(6.67);
  } else if ((A == 197) && (Z == 79)) {  // Gold
    // Default values from \iref{Schopper:2004qco}
    set_diffusiveness(0.523);
    set_nuclear_radius(6.55);
  } else if ((A == 129) && (Z == 54)) {  // Xenon
    // Default values.
    set_diffusiveness(0.59);
    set_nuclear_radius(5.36);
  } else if ((A == 63) && (Z == 29)) {  // Copper
    // Default values.
    set_diffusiveness(0.5977);
    set_nuclear_radius(4.20641);
  } else if (A == 96) {
    if (Z == 40) {  // Zirconium
      // Default values.
      set_diffusiveness(0.46);
      set_nuclear_radius(5.02);
    } else if (Z == 44) {  // Ruthenium
      // Default values.
      set_diffusiveness(0.46);
      set_nuclear_radius(5.085);
    } else {
      // radius and diffusiveness taken from \iref{Rybczynski:2013yba}
      set_diffusiveness(0.54);
      set_nuclear_radius(1.12 * std::pow(A, 1.0 / 3.0) -
                         0.86 * std::pow(A, -1.0 / 3.0));
    }
  } else {
    // saturation density already has reasonable default
    set_nuclear_radius(default_nuclear_radius());
    if (A <= 16) {
      set_diffusiveness(0.545);
    } else {
      // diffusiveness taken from \iref{Rybczynski:2013yba}
      set_diffusiveness(0.54);
    }
  }
}

set_parameters_from_config()

void smash::Nucleus::set_parameters_from_config ( Configuration &  config )

virtual

Sets the parameters of the Woods-Saxon according to manually added values in the configuration file.

Parameters

config

The configuration for this nucleus (projectile or target).

Definition at line 340 of file nucleus.cc.

                                                              {
  set_diffusiveness(static_cast<double>(config.take({"Diffusiveness"})));
  set_nuclear_radius(static_cast<double>(config.take({"Radius"})));
  // Saturation density (normalization for accept/reject sampling)
  set_saturation_density(
      static_cast<double>(config.take({"Saturation_Density"})));
}

generate_fermi_momenta()

void smash::Nucleus::generate_fermi_momenta ( )

virtual

Generates momenta according to Fermi motion for the nucleons.

For neutrons and protons Fermi momenta are calculated as \( p_{F} = (3 \pi^2 \rho)^{1/3}\), where \( rho \) is neutron density for neutrons and proton density for protons. The actual momenta \(p_x\), \(p_y\), \(p_z\) are uniformly distributed in the sphere with radius \(p_F\).

Reimplemented in smash::CustomNucleus.

Definition at line 348 of file nucleus.cc.

                                     {
  const int N_n = std::count_if(begin(), end(), [](const ParticleData i) {
    return i.pdgcode() == pdg::n;
  });
  const int N_p = std::count_if(begin(), end(), [](const ParticleData i) {
    return i.pdgcode() == pdg::p;
  });
  const FourVector nucleus_center = center();
  const int A = N_n + N_p;
  constexpr double pi2_3 = 3.0 * M_PI * M_PI;
  logg[LNucleus].debug() << N_n << " neutrons, " << N_p << " protons.";
 
  ThreeVector ptot = ThreeVector(0.0, 0.0, 0.0);
  for (auto i = begin(); i != end(); i++) {
    // Only protons and neutrons get Fermi momenta
    if (i->pdgcode() != pdg::p && i->pdgcode() != pdg::n) {
      if (i->is_baryon()) {
        logg[LNucleus].warn() << "No rule to calculate Fermi momentum "
                              << "for particle " << i->pdgcode();
      }
      continue;
    }
    const double r = (i->position() - nucleus_center).abs3();
    const double theta = (i->position().threevec().get_theta());
    const double phi = (i->position().threevec().get_phi());
    double rho = nucleon_density(r, std::cos(theta), phi);
 
    if (i->pdgcode() == pdg::p) {
      rho = rho * N_p / A;
    }
    if (i->pdgcode() == pdg::n) {
      rho = rho * N_n / A;
    }
    const double p =
        hbarc * std::pow(pi2_3 * rho * random::uniform(0.0, 1.0), 1.0 / 3.0);
    Angles phitheta;
    phitheta.distribute_isotropically();
    const ThreeVector ith_3momentum = phitheta.threevec() * p;
    ptot += ith_3momentum;
    i->set_3momentum(ith_3momentum);
    logg[LNucleus].debug() << "Particle: " << *i << ", pF[GeV]: "
                           << hbarc * std::pow(pi2_3 * rho, 1.0 / 3.0)
                           << " r[fm]: " << r
                           << " Nuclear radius[fm]: " << nuclear_radius_;
  }
  if (A == 0) {
    // No Fermi momenta should be assigned
    assert(ptot.x1() == 0.0 && ptot.x2() == 0.0 && ptot.x3() == 0.0);
  } else {
    /* Ensure zero total momentum of nucleus - redistribute ptot equally
     * among protons and neutrons */
    const ThreeVector centralizer = ptot / A;
    for (auto i = begin(); i != end(); i++) {
      if (i->pdgcode() == pdg::p || i->pdgcode() == pdg::n) {
        i->set_4momentum(i->pole_mass(),
                         i->momentum().threevec() - centralizer);
      }
    }
  }
}

boost()

void smash::Nucleus::boost

(

double

beta_scalar

)

Boosts the nuclei into the computational frame, such that the nucleons have the appropriate momentum and the nuclei are lorentz-contracted.

Note that the usual boost cannot be applied for nuclei, since the particles would end up with different times and the binding energy needs to be taken into account.

Parameters

[in]

beta_scalar

velocity in z-direction used for boost.

Definition at line 409 of file nucleus.cc.

                                      {
  double beta_squared = beta_scalar * beta_scalar;
  double one_over_gamma = std::sqrt(1.0 - beta_squared);
  double gamma = 1.0 / one_over_gamma;
  /* We are talking about a /passive/ lorentz transformation here, as
   * far as I can see, so we need to boost in the direction opposite to
   * where we want to go
   *     ( The vector we transform - p - stays unchanged, but we go into
   *       a system that moves with -beta. Now in this frame, it seems
   *       like p has been accelerated with +beta.
   *     ) */
  for (auto i = begin(); i != end(); i++) {
    /* a real Lorentz Transformation would leave the particles at
     * different times here, which we would then have to propagate back
     * to equal times. Since we know the result, we can simply multiply
     * the z-value with 1/gamma. */
    FourVector this_position = i->position();
    this_position.set_x3(this_position.x3() * one_over_gamma);
    i->set_4position(this_position);
    /* The simple Lorentz transformation of momenta does not take into account
     * that nucleus has binding energy. Here we apply the method used
     * in the JAM code \iref{Nara:1999dz}: p' = p_beam + gamma*p_F.
     * This formula is derived under assumption that all nucleons have
     * the same binding energy. */
    FourVector mom_i = i->momentum();
    i->set_4momentum(i->pole_mass(), mom_i.x1(), mom_i.x2(),
                     gamma * (beta_scalar * mom_i.x0() + mom_i.x3()));
  }
}

fill_from_list()

void smash::Nucleus::fill_from_list	(	const std::map< PdgCode, int > &	particle_list,
		int	testparticles
	)

Adds particles from a map PDG code => Number_of_particles_with_that_PDG_code to the nucleus.

E.g., the map [2212: 6, 2112: 7] initializes C-13 (6 protons and 7 neutrons). The particles are only created, no position or momenta are yet assigned. It is also possible to use any other PDG code, in addition to nucleons.

Parameters

[out]	particle_list	The particle slots that are created.
[in]	testparticles	Number of test particles to use.

Definition at line 439 of file nucleus.cc.

                                                {
  testparticles_ = testparticles;
  for (auto n = particle_list.cbegin(); n != particle_list.cend(); ++n) {
    const ParticleType &current_type = ParticleType::find(n->first);
    double current_mass = current_type.mass();
    for (unsigned int i = 0; i < n->second * testparticles_; i++) {
      // append particle to list and set its PDG code.
      particles_.emplace_back(current_type);
      particles_.back().set_4momentum(current_mass, 0.0, 0.0, 0.0);
    }
  }
}

shift()

void smash::Nucleus::shift	(	double	z_offset,
		double	x_offset,
		double	simulation_time
	)

Shifts the nucleus to correct impact parameter and z displacement.

Parameters

[in]	z_offset	is the shift in z-direction
[in]	x_offset	is the shift in x-direction
[in]	simulation_time	set the time and formation_time of each particle to this value.

Definition at line 453 of file nucleus.cc.

                                                                            {
  // Move the nucleus in z and x directions, and set the time.
  for (auto i = begin(); i != end(); i++) {
    FourVector this_position = i->position();
    this_position.set_x3(this_position.x3() + z_offset);
    this_position.set_x1(this_position.x1() + x_offset);
    this_position.set_x0(simulation_time);
    i->set_4position(this_position);
    i->set_formation_time(simulation_time);
  }
}

rotate()

virtual void smash::Nucleus::rotate ( )

inlinevirtual

Rotates the nucleus.

(Due to spherical symmetry of nondeformed nuclei, there is nothing to do.)

Reimplemented in smash::DeformedNucleus.

Definition at line 147 of file nucleus.h.

{}

copy_particles()

void smash::Nucleus::copy_particles

(

Particles *

particles

)

Copies the particles from this nucleus into the particle list.

Parameters

[out]

particles

Particle list with all constituents of a nucleus

Definition at line 465 of file nucleus.cc.

                                                          {
  for (auto p = begin(); p != end(); p++) {
    external_particles->insert(*p);
  }
}

size()

size_t smash::Nucleus::size ( ) const

inline

Number of numerical (=test-)particles in the nucleus:

Definition at line 157 of file nucleus.h.

{ return particles_.size(); }

number_of_particles()

size_t smash::Nucleus::number_of_particles ( ) const

inline

Number of physical particles in the nucleus:

Exceptions

TestparticleConfusion

if the number of the nucleons is not a multiple of testparticles_.

Definition at line 165 of file nucleus.h.

                                            {
    size_t nop = particles_.size() / testparticles_;
    /* If size() is not a multiple of testparticles_, this will throw an
     * error. */
    if (nop * testparticles_ != particles_.size()) {
      throw TestparticleConfusion(
          "Number of test particles and test particles"
          "per particle are incompatible.");
    }
    return nop;
  }

number_of_protons()

size_t smash::Nucleus::number_of_protons ( ) const

inline

Number of physical protons in the nucleus:

Returns

number of protons

Exceptions

Testparticleconfusion

if the number of the protons is not a multiple of testparticles_.

Definition at line 184 of file nucleus.h.

                                          {
    size_t proton_counter = 0;
    /* If n_protons is not a multiple of testparticles_, this will throw an
     * error. */
    for (auto &particle : particles_) {
      if (particle.type().pdgcode() == pdg::p) {
        proton_counter++;
      }
    }
 
    size_t n_protons = proton_counter / testparticles_;
 
    if (n_protons * testparticles_ != proton_counter) {
      throw TestparticleConfusion(
          "Number of test protons and test particles"
          "per proton are incompatible.");
    }
 
    return n_protons;
  }

center()

FourVector smash::Nucleus::center

(

)

const

Calculate geometrical center of the nucleus.

Returns

\(\mathbf{r}_s = \frac{1}{N} \sum_{i=1}^N \mathbf{r}_i\) (for a nucleus with N particles that are at the positions \(\mathbf{r}_i\)).

Definition at line 471 of file nucleus.cc.

                                 {
  FourVector centerpoint(0.0, 0.0, 0.0, 0.0);
  for (auto p = cbegin(); p != cend(); p++) {
    centerpoint += p->position();
  }
  centerpoint /= size();
  return centerpoint;
}

set_label()

void smash::Nucleus::set_label ( BelongsTo  label )

inline

Sets target / projectile labels on nucleons.

Definition at line 213 of file nucleus.h.

                                  {
    for (ParticleData &data : particles_) {
      data.set_belongs_to(label);
    }
  }

align_center()

void smash::Nucleus::align_center ( )

inline

Shifts the nucleus so that its center is at (0,0,0)

See also

center()

Definition at line 223 of file nucleus.h.

                      {
    FourVector centerpoint = center();
    for (auto p = particles_.begin(); p != particles_.end(); ++p) {
      p->set_4position(p->position() - centerpoint);
    }
  }

nucleon_density()

double smash::Nucleus::nucleon_density ( double  r, double  , double    ) const

virtual

Return the Woods-Saxon probability density for the given position.

This corresponds to the nuclear density at the very same position.

Parameters

[in]

r

The radius at which to sample

Returns

The Woods-Saxon density

Reimplemented in smash::DeformedNucleus.

Definition at line 487 of file nucleus.cc.

                                                              {
  return get_saturation_density() /
         (std::exp((r - nuclear_radius_) / diffusiveness_) + 1.);
}

nucleon_density_unnormalized()

double smash::Nucleus::nucleon_density_unnormalized ( double  r, double  , double    ) const

virtual

Return the unnormalized Woods-Saxon distribution for the given position without deformation.

Parameters

[in]

r

The radius

Returns

The unnormalized Woods-Saxon distribution

Reimplemented in smash::DeformedNucleus.

Definition at line 492 of file nucleus.cc.

                                                                           {
  return 1.0 / (std::exp((r - nuclear_radius_) / diffusiveness_) + 1.);
}

calculate_saturation_density()

double smash::Nucleus::calculate_saturation_density ( ) const

virtual

Returns

the normalized ground state density for the corresponding Woods-Saxon parameter. This is done by integrating the Woods-Saxon distribution and setting the normalization such that the integral of the Woods-Saxon distribution yields the number of particles in the nucleus \(\int\rho(r)d^3r = N_{particles}\).

Reimplemented in smash::DeformedNucleus.

Definition at line 496 of file nucleus.cc.

                                                   {
  Integrator2d integrate;
  // Transform integral from (0, oo) to (0, 1) via r = (1 - t) / t.
  // To prevent overflow, the integration is only performed to t = 0.01 which
  // corresponds to r = 99fm. Additionally the precision settings in the
  // Integrator2d scheme are equally important. However both these point affect
  // the result only after the seventh digit which should not be relevant here.
  const auto result = integrate(0.01, 1, -1, 1, [&](double t, double cosx) {
    const double r = (1 - t) / t;
    return twopi * std::pow(r, 2.0) *
           nucleon_density_unnormalized(r, cosx, 0.0) / std::pow(t, 2.0);
  });
  const auto rho0 = number_of_particles() / result.value();
  return rho0;
}

set_saturation_density()

virtual void smash::Nucleus::set_saturation_density ( double  density )

inlinevirtual

Sets the saturation density of the nucleus.

See also

saturation_density_

Definition at line 261 of file nucleus.h.

                                                      {
    saturation_density_ = density;
  }

random_euler_angles()

void smash::Nucleus::random_euler_angles ( )

protected

Randomly generate Euler angles.

Necessary for rotation of deformed and custom nuclei, whenever a new nucleus of this kind is initialized.

Definition at line 480 of file nucleus.cc.

                                  {
  // Sample euler_theta_ such that cos(theta) is uniform
  euler_phi_ = twopi * random::uniform(0., 1.);
  euler_theta_ = std::acos(2 * random::uniform(0., 1.) - 1);
  euler_psi_ = twopi * random::uniform(0., 1.);
}

begin()

std::vector<ParticleData>::iterator smash::Nucleus::begin ( )

inline

For iterators over the particle list:

Definition at line 309 of file nucleus.h.

                                                 {
    return particles_.begin();
  }

end()

std::vector<ParticleData>::iterator smash::Nucleus::end ( )

inline

For iterators over the particle list:

Definition at line 313 of file nucleus.h.

{ return particles_.end(); }

cbegin()

std::vector<ParticleData>::const_iterator smash::Nucleus::cbegin ( ) const

inline

For const iterators over the particle list:

Definition at line 315 of file nucleus.h.

                                                              {
    return particles_.cbegin();
  }

cend()

std::vector<ParticleData>::const_iterator smash::Nucleus::cend ( ) const

inline

For const iterators over the particle list:

Definition at line 319 of file nucleus.h.

                                                            {
    return particles_.cend();
  }

set_diffusiveness()

void smash::Nucleus::set_diffusiveness ( double  diffuse )

inline

Sets the diffusiveness of the nucleus.

See also

diffusiveness_

Definition at line 326 of file nucleus.h.

{ diffusiveness_ = diffuse; }

get_diffusiveness()

double smash::Nucleus::get_diffusiveness ( ) const

inline

Returns

the diffusiveness of the nucleus

See also

diffusiveness_

Definition at line 331 of file nucleus.h.

{ return diffusiveness_; }

get_saturation_density()

double smash::Nucleus::get_saturation_density ( ) const

inline

Returns

the saturation density of the nucleus

See also

saturation_density_

Definition at line 336 of file nucleus.h.

{ return saturation_density_; }

default_nuclear_radius()

double smash::Nucleus::default_nuclear_radius ( )

inline

Default nuclear radius calculated as:

• \( r = r_\mathrm{proton} \ A^{1/3} \qquad \qquad \qquad \ \) for A <= 16
• \( r = 1.12 \ A^{1/3} - 0.86 \ A^{-1/3} \qquad \) for A > 16

Returns

default radius for the nucleus in fm

Definition at line 344 of file nucleus.h.

                                         {
    int A = number_of_particles();
 
    if (A <= 16) {
      // radius: rough guess for all nuclei not listed explicitly with A <= 16
      return (proton_radius_ * std::cbrt(A));
    } else {
      // radius taken from \iref{Rybczynski:2013yba}
      return (1.12 * std::pow(A, 1.0 / 3.0) - 0.86 * std::pow(A, -1.0 / 3.0));
    }
  }

set_nuclear_radius()

void smash::Nucleus::set_nuclear_radius ( double  rad )

inline

Sets the nuclear radius.

See also

nuclear_radius

Definition at line 359 of file nucleus.h.

{ nuclear_radius_ = rad; }

get_nuclear_radius()

double smash::Nucleus::get_nuclear_radius ( ) const

inline

Returns

the nuclear radius

See also

nuclear_radius

Definition at line 364 of file nucleus.h.

{ return nuclear_radius_; }

Member Data Documentation

diffusiveness_

double smash::Nucleus::diffusiveness_

private

Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0, we obtain a hard sphere).

Definition at line 275 of file nucleus.h.

nuclear_radius_

double smash::Nucleus::nuclear_radius_

private

Nuclear radius of this nucleus.

Definition at line 277 of file nucleus.h.

proton_radius_

double smash::Nucleus::proton_radius_ = 1.2

private

Single proton radius in fm.

See also

default_nuclear_radius

Definition at line 282 of file nucleus.h.

testparticles_

size_t smash::Nucleus::testparticles_ = 1

private

Number of testparticles per physical particle.

Definition at line 284 of file nucleus.h.

particles_

std::vector<ParticleData> smash::Nucleus::particles_

protected

Particles associated with this nucleus.

Definition at line 288 of file nucleus.h.

saturation_density_

double smash::Nucleus::saturation_density_ = nuclear_density

protected

Saturation density of this nucleus.

Definition at line 292 of file nucleus.h.

euler_phi_

double smash::Nucleus::euler_phi_

protected

Euler angel phi.

Definition at line 301 of file nucleus.h.

euler_theta_

double smash::Nucleus::euler_theta_

protected

Euler angel theta.

Definition at line 303 of file nucleus.h.

euler_psi_

double smash::Nucleus::euler_psi_

protected

Euler angel psi.

Definition at line 305 of file nucleus.h.

---

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

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