Version: SMASH-1.5
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:
[legend]
Collaboration diagram for smash::Nucleus:
[legend]

Classes

struct  TestparticleConfusion
 

Public Member Functions

 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)
 
virtual ~Nucleus ()=default
 
double mass () const
 
virtual ThreeVector distribute_nucleon () const
 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...
 
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...
 
FourVector center () const
 Calculate geometrical center of the nucleus. More...
 
void align_center ()
 Shifts the nucleus so that its center is at (0,0,0) 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
 
void set_saturation_density (double density)
 Sets the saturation density of the nucleus. More...
 
double get_saturation_density () const
 
double default_nuclear_radius ()
 
void set_nuclear_radius (double rad)
 Sets the nuclear radius. More...
 
double get_nuclear_radius () const
 

Private Attributes

double diffusiveness_
 Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0, we obtain a hard sphere). More...
 
double saturation_density_ = nuclear_density
 Saturation density of this nucleus. 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...
 
std::vector< ParticleDataparticles_
 Particles associated with this nucleus. 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/2]

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]configcontains the parameters from the inputfile on the numbers of particles with a certain PDG code
[in]nTestnumber of testparticles

Definition at line 32 of file nucleus.cc.

32  {
33  // Fill nuclei with particles.
34  std::map<PdgCode, int> part = config.take({"Particles"});
35  fill_from_list(part, nTest);
36  // Look for user-defined values or take the default parameters.
37  if (config.has_value({"Diffusiveness"}) && config.has_value({"Radius"}) &&
38  config.has_value({"Saturation_Density"})) {
40  } else if (!config.has_value({"Diffusiveness"}) &&
41  !config.has_value({"Radius"}) &&
42  !config.has_value({"Saturation_Density"})) {
44  } else {
45  throw std::invalid_argument(
46  "Diffussiveness, Radius and Saturation_Density "
47  "required to manually configure the Woods-Saxon"
48  " distribution. Only one/two were provided. \n"
49  "Providing none of the above mentioned "
50  "parameters automatically configures the "
51  "distribution based on the atomic number.");
52  }
53 }
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...
Definition: nucleus.cc:433
virtual void set_parameters_automatic()
Sets the deformation parameters of the Woods-Saxon distribution according to the current mass number...
Definition: nucleus.cc:282
virtual void set_parameters_from_config(Configuration &config)
Sets the parameters of the Woods-Saxon according to manually added values in the configuration file...
Definition: nucleus.cc:334
Here is the call graph for this function:

◆ Nucleus() [2/2]

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

Definition at line 27 of file nucleus.cc.

27  {
28  fill_from_list(particle_list, nTest);
30 }
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...
Definition: nucleus.cc:433
virtual void set_parameters_automatic()
Sets the deformation parameters of the Woods-Saxon distribution according to the current mass number...
Definition: nucleus.cc:282
Here is the call graph for this function:

◆ ~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.

55  {
56  double total_mass = 0.;
57  for (auto i = cbegin(); i != cend(); i++) {
58  total_mass += i->momentum().abs();
59  }
60  return total_mass / (testparticles_ + 0.0);
61 }
std::vector< ParticleData >::const_iterator cend() const
For const iterators over the particle list:
Definition: nucleus.h:219
std::vector< ParticleData >::const_iterator cbegin() const
For const iterators over the particle list:
Definition: nucleus.h:215
size_t testparticles_
Number of testparticles per physical particle.
Definition: nucleus.h:203
Here is the call graph for this function:

◆ distribute_nucleon()

ThreeVector smash::Nucleus::distribute_nucleon ( ) const
virtual

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

Woods-Saxon-distribution

\(\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.

The distribution

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

\[\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: \(p^{({\rm I})}\)

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: \(p^{({\rm II})}\)

Again, from one only:

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

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

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

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

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: \(p^{({\rm IV})}\)

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.

Definition at line 215 of file nucleus.cc.

215  {
216  // Get the solid angle of the nucleon.
217  Angles dir;
218  dir.distribute_isotropically();
219  // diffusiveness_ zero or negative? Use hard sphere.
220  if (almost_equal(diffusiveness_, 0.)) {
221  return dir.threevec() * nuclear_radius_ * std::cbrt(random::canonical());
222  }
223  if (almost_equal(nuclear_radius_, 0.)) {
224  return smash::ThreeVector();
225  }
226  double radius_scaled = nuclear_radius_ / diffusiveness_;
227  double prob_range1 = 1.0;
228  double prob_range2 = 3. / radius_scaled;
229  double prob_range3 = 2. * prob_range2 / radius_scaled;
230  double prob_range4 = 1. * prob_range3 / radius_scaled;
231  double ranges234 = prob_range2 + prob_range3 + prob_range4;
232  double t;
234  do {
235  double which_range = random::uniform(-prob_range1, ranges234);
236  if (which_range < 0.0) {
237  t = radius_scaled * (std::cbrt(random::canonical()) - 1.);
238  } else {
239  t = -std::log(random::canonical());
240  if (which_range >= prob_range2) {
241  t -= std::log(random::canonical());
242  if (which_range >= prob_range2 + prob_range3) {
243  t -= std::log(random::canonical());
244  }
245  }
246  }
254  } while (random::canonical() > 1. / (1. + std::exp(-std::abs(t))));
256  double position_scaled = t + radius_scaled;
257  double position = position_scaled * diffusiveness_;
258  return dir.threevec() * position;
259 }
The ThreeVector class represents a physical three-vector with the components .
Definition: threevector.h:30
bool almost_equal(const N x, const N y)
Checks two numbers for relative approximate equality.
Definition: numerics.h:42
T canonical()
Definition: random.h:110
T uniform(T min, T max)
Definition: random.h:85
double nuclear_radius_
Nuclear radius of this nucleus.
Definition: nucleus.h:196
double diffusiveness_
Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0...
Definition: nucleus.h:192
Here is the call graph for this function:
Here is the caller graph for this function:

◆ woods_saxon()

double smash::Nucleus::woods_saxon ( double  x)

Woods-Saxon distribution.

Parameters
[in]xthe position at which to evaluate the function
Returns
un-normalized Woods-saxon probability

Definition at line 261 of file nucleus.cc.

261  {
262  return r * r / (std::exp((r - nuclear_radius_) / diffusiveness_) + 1);
263 }
double nuclear_radius_
Nuclear radius of this nucleus.
Definition: nucleus.h:196
double diffusiveness_
Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0...
Definition: nucleus.h:192

◆ arrange_nucleons()

void smash::Nucleus::arrange_nucleons ( )

Sets the positions of the nucleons inside a nucleus.

Definition at line 265 of file nucleus.cc.

265  {
266  for (auto i = begin(); i != end(); i++) {
267  // Initialize momentum
268  i->set_4momentum(i->pole_mass(), 0.0, 0.0, 0.0);
269  /* Sampling the Woods-Saxon, get the radial
270  * position and solid angle for the nucleon. */
271  ThreeVector pos = distribute_nucleon();
272 
273  // Set the position of the nucleon.
274  i->set_4position(FourVector(0.0, pos));
275  }
276 
277  // Recenter and rotate
278  align_center();
279  rotate();
280 }
virtual void rotate()
Rotates the nucleus.
Definition: nucleus.h:134
virtual ThreeVector distribute_nucleon() const
The distribution of return values from this function is according to a spherically symmetric Woods-Sa...
Definition: nucleus.cc:215
std::vector< ParticleData >::iterator begin()
For iterators over the particle list:
Definition: nucleus.h:209
void align_center()
Shifts the nucleus so that its center is at (0,0,0)
Definition: nucleus.h:175
std::vector< ParticleData >::iterator end()
For iterators over the particle list:
Definition: nucleus.h:213
Here is the call graph for this function:

◆ 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 DeJager:1987qc. They are in agreement with MC-Glauber models such as GLISSANDO (see Rybczynski:2013yba) and TGlauber MC (see Loizides:2017ack)

Definition at line 282 of file nucleus.cc.

282  {
284  switch (A) {
285  case 1: // single particle
286  /* In case of testparticles, an infinite reaction loop will be
287  * avoided by a small finite spread according to a single particles
288  * 'nucleus'. The proper solution will be to introduce parallel
289  * ensembles. */
291  ? 0.
292  : 1. - std::exp(-(testparticles_ - 1.) * 0.1));
293  set_diffusiveness(testparticles_ == 1 ? -1. : 0.02);
294  break;
295  case 238: // Uranium
296  // Default values.
297  set_diffusiveness(0.556);
298  set_nuclear_radius(6.86);
299  set_saturation_density(0.166);
300  break;
301  case 208: // Lead
302  // Default values.
303  set_diffusiveness(0.54);
304  set_nuclear_radius(6.67);
305  set_saturation_density(0.161);
306  break;
307  case 197: // Gold
308  // Default values.
309  set_diffusiveness(0.535);
310  set_nuclear_radius(6.38);
311  set_saturation_density(0.1695);
312  break;
313  case 63: // Copper
314  // Default values.
315  set_diffusiveness(0.5977);
316  set_nuclear_radius(4.20641);
317  set_saturation_density(0.1686);
318  break;
319  default:
320  if (A <= 16) {
321  // radius: rough guess for all nuclei not listed explicitly with A <= 16
322  // saturation density already has reasonable default
323  set_nuclear_radius(1.2 * std::cbrt(A));
324  set_diffusiveness(0.545);
325  } else {
326  // radius and diffusiveness taken from \iref{Rybczynski:2013yba}
327  set_diffusiveness(0.54);
328  set_nuclear_radius(1.12 * pow(std::cbrt(A), 1.0 / 3.0) -
329  0.86 * pow(std::cbrt(A), -1.0 / 3.0));
330  }
331  }
332 }
void set_diffusiveness(double diffuse)
Sets the diffusiveness of the nucleus.
Definition: nucleus.h:226
void set_nuclear_radius(double rad)
Sets the nuclear radius.
Definition: nucleus.h:256
void set_saturation_density(double density)
Sets the saturation density of the nucleus.
Definition: nucleus.h:236
size_t testparticles_
Number of testparticles per physical particle.
Definition: nucleus.h:203
size_t number_of_particles() const
Number of physical particles in the nucleus:
Definition: nucleus.h:152
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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
configThe configuration for this nucleus (projectile or target).

Definition at line 334 of file nucleus.cc.

334  {
335  set_diffusiveness(static_cast<double>(config.take({"Diffusiveness"})));
336  set_nuclear_radius(static_cast<double>(config.take({"Radius"})));
337  // Saturation density (normalization for accept/reject sampling)
339  static_cast<double>(config.take({"Saturation_Density"})));
340 }
void set_diffusiveness(double diffuse)
Sets the diffusiveness of the nucleus.
Definition: nucleus.h:226
void set_nuclear_radius(double rad)
Sets the nuclear radius.
Definition: nucleus.h:256
void set_saturation_density(double density)
Sets the saturation density of the nucleus.
Definition: nucleus.h:236
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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::DeformedNucleus.

Definition at line 342 of file nucleus.cc.

342  {
343  const int N_n = std::count_if(begin(), end(), [](const ParticleData i) {
344  return i.pdgcode() == pdg::n;
345  });
346  const int N_p = std::count_if(begin(), end(), [](const ParticleData i) {
347  return i.pdgcode() == pdg::p;
348  });
349  const FourVector nucleus_center = center();
350  const int A = N_n + N_p;
351  constexpr double pi2_3 = 3.0 * M_PI * M_PI;
352  const auto &log = logger<LogArea::Nucleus>();
353 
354  log.debug() << N_n << " neutrons, " << N_p << " protons.";
355 
356  ThreeVector ptot = ThreeVector(0.0, 0.0, 0.0);
357  for (auto i = begin(); i != end(); i++) {
358  // Only protons and neutrons get Fermi momenta
359  if (i->pdgcode() != pdg::p && i->pdgcode() != pdg::n) {
360  if (i->is_baryon()) {
361  log.warn() << "No rule to calculate Fermi momentum "
362  << "for particle " << i->pdgcode();
363  }
364  continue;
365  }
366  const double r = (i->position() - nucleus_center).abs3();
367  double rho = nuclear_density /
368  (std::exp((r - nuclear_radius_) / diffusiveness_) + 1.);
369  if (i->pdgcode() == pdg::p) {
370  rho = rho * N_p / A;
371  }
372  if (i->pdgcode() == pdg::n) {
373  rho = rho * N_n / A;
374  }
375  const double p =
376  hbarc * std::pow(pi2_3 * rho * random::uniform(0.0, 1.0), 1.0 / 3.0);
377  Angles phitheta;
378  phitheta.distribute_isotropically();
379  const ThreeVector ith_3momentum = phitheta.threevec() * p;
380  ptot += ith_3momentum;
381  i->set_3momentum(ith_3momentum);
382  log.debug() << "Particle: " << *i
383  << ", pF[GeV]: " << hbarc * std::pow(pi2_3 * rho, 1.0 / 3.0)
384  << " r[fm]: " << r
385  << " Nuclear radius[fm]: " << nuclear_radius_;
386  }
387  if (A == 0) {
388  // No Fermi momenta should be assigned
389  assert(ptot.x1() == 0.0 && ptot.x2() == 0.0 && ptot.x3() == 0.0);
390  } else {
391  /* Ensure zero total momentum of nucleus - redistribute ptot equally
392  * among protons and neutrons */
393  const ThreeVector centralizer = ptot / A;
394  for (auto i = begin(); i != end(); i++) {
395  if (i->pdgcode() == pdg::p || i->pdgcode() == pdg::n) {
396  i->set_4momentum(i->pole_mass(),
397  i->momentum().threevec() - centralizer);
398  }
399  }
400  }
401 }
FourVector center() const
Calculate geometrical center of the nucleus.
Definition: nucleus.cc:465
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:25
T uniform(T min, T max)
Definition: random.h:85
double nuclear_radius_
Nuclear radius of this nucleus.
Definition: nucleus.h:196
double diffusiveness_
Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0...
Definition: nucleus.h:192
constexpr int p
Proton.
constexpr double nuclear_density
Ground state density of symmetric nuclear matter [fm^-3].
Definition: constants.h:42
std::vector< ParticleData >::iterator begin()
For iterators over the particle list:
Definition: nucleus.h:209
constexpr int n
Neutron.
std::vector< ParticleData >::iterator end()
For iterators over the particle list:
Definition: nucleus.h:213
Here is the call graph for this function:

◆ 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_scalarvelocity in z-direction used for boost.

Definition at line 403 of file nucleus.cc.

403  {
404  double beta_squared = beta_scalar * beta_scalar;
405  double one_over_gamma = std::sqrt(1.0 - beta_squared);
406  double gamma = 1.0 / one_over_gamma;
407  /* We are talking about a /passive/ lorentz transformation here, as
408  * far as I can see, so we need to boost in the direction opposite to
409  * where we want to go
410  * ( The vector we transform - p - stays unchanged, but we go into
411  * a system that moves with -beta. Now in this frame, it seems
412  * like p has been accelerated with +beta.
413  * ) */
414  for (auto i = begin(); i != end(); i++) {
415  /* a real Lorentz Transformation would leave the particles at
416  * different times here, which we would then have to propagate back
417  * to equal times. Since we know the result, we can simply multiply
418  * the z-value with 1/gamma. */
419  FourVector this_position = i->position();
420  this_position.set_x3(this_position.x3() * one_over_gamma);
421  i->set_4position(this_position);
422  /* The simple Lorentz transformation of momenta does not take into account
423  * that nucleus has binding energy. Here we apply the method used
424  * in the JAM code \iref{Nara:1999dz}: p' = p_beam + gamma*p_F.
425  * This formula is derived under assumption that all nucleons have
426  * the same binding energy. */
427  FourVector mom_i = i->momentum();
428  i->set_4momentum(i->pole_mass(), mom_i.x1(), mom_i.x2(),
429  gamma * (beta_scalar * mom_i.x0() + mom_i.x3()));
430  }
431 }
std::vector< ParticleData >::iterator begin()
For iterators over the particle list:
Definition: nucleus.h:209
std::vector< ParticleData >::iterator end()
For iterators over the particle list:
Definition: nucleus.h:213
Here is the call graph for this function:

◆ 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_listThe particle slots that are created.
[in]testparticlesNumber of test particles to use.

Definition at line 433 of file nucleus.cc.

434  {
435  testparticles_ = testparticles;
436  for (auto n = particle_list.cbegin(); n != particle_list.cend(); ++n) {
437  const ParticleType &current_type = ParticleType::find(n->first);
438  double current_mass = current_type.mass();
439  for (unsigned int i = 0; i < n->second * testparticles_; i++) {
440  // append particle to list and set its PDG code.
441  particles_.emplace_back(current_type);
442  particles_.back().set_4momentum(current_mass, 0.0, 0.0, 0.0);
443  }
444  }
445 }
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
constexpr int n
Neutron.
size_t testparticles_
Number of testparticles per physical particle.
Definition: nucleus.h:203
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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_offsetis the shift in z-direction
[in]x_offsetis the shift in x-direction
[in]simulation_timeset the time and formation_time of each particle to this value.

Definition at line 447 of file nucleus.cc.

447  {
448  // Move the nucleus in z and x directions, and set the time.
449  for (auto i = begin(); i != end(); i++) {
450  FourVector this_position = i->position();
451  this_position.set_x3(this_position.x3() + z_offset);
452  this_position.set_x1(this_position.x1() + x_offset);
453  this_position.set_x0(simulation_time);
454  i->set_4position(this_position);
455  i->set_formation_time(simulation_time);
456  }
457 }
std::vector< ParticleData >::iterator begin()
For iterators over the particle list:
Definition: nucleus.h:209
std::vector< ParticleData >::iterator end()
For iterators over the particle list:
Definition: nucleus.h:213
Here is the call graph for this function:

◆ 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 134 of file nucleus.h.

134 {}
Here is the caller graph for this function:

◆ copy_particles()

void smash::Nucleus::copy_particles ( Particles particles)

Copies the particles from this nucleus into the particle list.

Parameters
[out]particlesParticle list with all constituents of a nucleus

Definition at line 459 of file nucleus.cc.

459  {
460  for (auto p = begin(); p != end(); p++) {
461  external_particles->insert(*p);
462  }
463 }
constexpr int p
Proton.
std::vector< ParticleData >::iterator begin()
For iterators over the particle list:
Definition: nucleus.h:209
std::vector< ParticleData >::iterator end()
For iterators over the particle list:
Definition: nucleus.h:213
Here is the call graph for this function:

◆ size()

size_t smash::Nucleus::size ( ) const
inline

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

Definition at line 144 of file nucleus.h.

144 { return particles_.size(); }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the caller graph for this function:

◆ number_of_particles()

size_t smash::Nucleus::number_of_particles ( ) const
inline

Number of physical particles in the nucleus:

Exceptions
TestparticleConfusionif the number of the nucleons is not a multiple of testparticles_.

Definition at line 152 of file nucleus.h.

152  {
153  size_t nop = particles_.size() / testparticles_;
154  /* if size() is not a multiple of testparticles_, this will throw an
155  * error. */
156  if (nop * testparticles_ != particles_.size()) {
157  throw TestparticleConfusion(
158  "Number of test particles and test particles"
159  "per particle are incompatible");
160  }
161  return nop;
162  }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
size_t testparticles_
Number of testparticles per physical particle.
Definition: nucleus.h:203
Here is the caller graph for this function:

◆ center()

FourVector smash::Nucleus::center ( ) const

Calculate geometrical center of the nucleus.

Returns
\(\vec r_s = \frac{1}{N} \sum_{i=1}^N \vec r_i\) (for a nucleus with N particles that are at the positions \(\vec r_i\)).

Definition at line 465 of file nucleus.cc.

465  {
466  FourVector centerpoint(0.0, 0.0, 0.0, 0.0);
467  for (auto p = cbegin(); p != cend(); p++) {
468  centerpoint += p->position();
469  }
470  centerpoint /= size();
471  return centerpoint;
472 }
size_t size() const
Number of numerical (=test-)particles in the nucleus:
Definition: nucleus.h:144
std::vector< ParticleData >::const_iterator cend() const
For const iterators over the particle list:
Definition: nucleus.h:219
std::vector< ParticleData >::const_iterator cbegin() const
For const iterators over the particle list:
Definition: nucleus.h:215
constexpr int p
Proton.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 175 of file nucleus.h.

175  {
176  FourVector centerpoint = center();
177  for (auto p = particles_.begin(); p != particles_.end(); ++p) {
178  p->set_4position(p->position() - centerpoint);
179  }
180  }
FourVector center() const
Calculate geometrical center of the nucleus.
Definition: nucleus.cc:465
constexpr int p
Proton.
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the call graph for this function:
Here is the caller graph for this function:

◆ begin()

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

For iterators over the particle list:

Definition at line 209 of file nucleus.h.

209  {
210  return particles_.begin();
211  }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the caller graph for this function:

◆ end()

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

For iterators over the particle list:

Definition at line 213 of file nucleus.h.

213 { return particles_.end(); }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the caller graph for this function:

◆ cbegin()

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

For const iterators over the particle list:

Definition at line 215 of file nucleus.h.

215  {
216  return particles_.cbegin();
217  }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the caller graph for this function:

◆ cend()

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

For const iterators over the particle list:

Definition at line 219 of file nucleus.h.

219  {
220  return particles_.cend();
221  }
std::vector< ParticleData > particles_
Particles associated with this nucleus.
Definition: nucleus.h:205
Here is the caller graph for this function:

◆ set_diffusiveness()

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

Sets the diffusiveness of the nucleus.

See also
diffusiveness_

Definition at line 226 of file nucleus.h.

226 { diffusiveness_ = diffuse; }
double diffusiveness_
Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0...
Definition: nucleus.h:192
Here is the caller graph for this function:

◆ get_diffusiveness()

double smash::Nucleus::get_diffusiveness ( ) const
inline
Returns
the diffusiveness of the nucleus
See also
diffusiveness_

Definition at line 231 of file nucleus.h.

231 { return diffusiveness_; }
double diffusiveness_
Diffusiveness of Woods-Saxon distribution of this nucleus in fm (for diffusiveness_ == 0...
Definition: nucleus.h:192
Here is the caller graph for this function:

◆ set_saturation_density()

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

Sets the saturation density of the nucleus.

See also
saturation_density_

Definition at line 236 of file nucleus.h.

236  {
237  saturation_density_ = density;
238  }
double saturation_density_
Saturation density of this nucleus.
Definition: nucleus.h:194
Here is the caller graph for this function:

◆ 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 243 of file nucleus.h.

243 { return saturation_density_; }
double saturation_density_
Saturation density of this nucleus.
Definition: nucleus.h:194
Here is the caller graph for this function:

◆ default_nuclear_radius()

double smash::Nucleus::default_nuclear_radius ( )
inline
Returns
a default radius for the nucleus Nuclear radius is calculated with the proton radius times the third root of the number of nucleons.

Definition at line 249 of file nucleus.h.

249  {
250  return proton_radius_ * std::pow(number_of_particles(), 1. / 3.);
251  }
double proton_radius_
Single proton radius in fm.
Definition: nucleus.h:201
size_t number_of_particles() const
Number of physical particles in the nucleus:
Definition: nucleus.h:152
Here is the call graph for this function:

◆ set_nuclear_radius()

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

Sets the nuclear radius.

See also
nuclear_radius

Definition at line 256 of file nucleus.h.

256 { nuclear_radius_ = rad; }
double nuclear_radius_
Nuclear radius of this nucleus.
Definition: nucleus.h:196
Here is the caller graph for this function:

◆ get_nuclear_radius()

double smash::Nucleus::get_nuclear_radius ( ) const
inline
Returns
the nuclear radius
See also
nuclear_radius

Definition at line 261 of file nucleus.h.

261 { return nuclear_radius_; }
double nuclear_radius_
Nuclear radius of this nucleus.
Definition: nucleus.h:196
Here is the caller graph for this function:

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 192 of file nucleus.h.

◆ saturation_density_

double smash::Nucleus::saturation_density_ = nuclear_density
private

Saturation density of this nucleus.

Definition at line 194 of file nucleus.h.

◆ nuclear_radius_

double smash::Nucleus::nuclear_radius_
private

Nuclear radius of this nucleus.

Definition at line 196 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 201 of file nucleus.h.

◆ testparticles_

size_t smash::Nucleus::testparticles_ = 1
private

Number of testparticles per physical particle.

Definition at line 203 of file nucleus.h.

◆ particles_

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

Particles associated with this nucleus.

Definition at line 205 of file nucleus.h.


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