Version: SMASH-3.1
smash::Experiment< Modus > Class Template Reference

#include <experiment.h>

template<typename Modus>
class smash::Experiment< Modus >

The main class, where the simulation of an experiment is executed.

The Experiment class owns all data (maybe indirectly) relevant for the execution of the experiment simulation. The experiment can be conducted in different running modi. Since the abstraction of these differences should not incur any overhead, the design is built around the Policy pattern.

The Policy pattern was defined by Andrei Alexandrescu in his book "Modern C++ Design: Generic Programming and Design Patterns Applied". Addison-Wesley:

A policy defines a class interface or a class template interface. The interface consists of one or all of the following: inner type definitions, member functions, and member variables.

The policy pattern can also be understood as a compile-time variant of the strategy pattern.

The Modus template parameter defines the "policy" of the Experiment class. It determines several aspects of the experiment execution at compile time. The original strategy pattern would select these differences at run time, thus incurring an overhead. This overhead becomes severe in cases where calls to strategy/policy functions are done very frequently. Using the policy pattern, the compiler can fully optimize: It creates a new instance of all functions in Experiment for all different Modus types.

Definition at line 188 of file experiment.h.

Inheritance diagram for smash::Experiment< Modus >:
smash::ExperimentBase

Public Member Functions

void run () override
 Runs the experiment. More...
 
 Experiment (Configuration &config, const std::filesystem::path &output_path)
 Create a new Experiment. More...
 
void initialize_new_event ()
 This is called in the beginning of each event. More...
 
void run_time_evolution (const double t_end, ParticleList &&add_plist={}, ParticleList &&remove_plist={})
 Runs the time evolution of an event with fixed-size time steps or without timesteps, from action to actions. More...
 
void do_final_decays ()
 Performs the final decays of an event. More...
 
void final_output ()
 Output at the end of an event. More...
 
Particlesfirst_ensemble ()
 Provides external access to SMASH particles. More...
 
std::vector< Particles > * all_ensembles ()
 Getter for all ensembles. More...
 
Modus * modus ()
 Provides external access to SMASH calculation modus. More...
 
void increase_event_number ()
 Increases the event number by one. More...
 
- Public Member Functions inherited from smash::ExperimentBase
 ExperimentBase ()=default
 
virtual ~ExperimentBase ()=default
 The virtual destructor avoids undefined behavior when destroying derived objects. More...
 

Private Member Functions

bool perform_action (Action &action, int i_ensemble, bool include_pauli_blocking=true)
 Perform the given action. More...
 
void create_output (const std::string &format, const std::string &content, const std::filesystem::path &output_path, const OutputParameters &par)
 Create a list of output files. More...
 
void propagate_and_shine (double to_time, Particles &particles)
 Propagate all particles until time to_time without any interactions and shine dileptons. More...
 
void run_time_evolution_timestepless (Actions &actions, int i_ensemble, const double end_time_propagation)
 Performs all the propagations and actions during a certain time interval neglecting the influence of the potentials. More...
 
void intermediate_output ()
 Intermediate output during an event. More...
 
void update_potentials ()
 Recompute potentials on lattices if necessary. More...
 
double compute_min_cell_length (double dt) const
 Calculate the minimal size for the grid cells such that the ScatterActionsFinder will find all collisions within the maximal transverse distance (which is determined by the maximal cross section). More...
 
double next_output_time () const
 Shortcut for next output time. More...
 
void count_nonempty_ensembles ()
 Counts the number of ensembles in wich interactions took place at the end of an event. More...
 
bool is_finished ()
 Checks wether the desired number events have been calculated. More...
 

Private Attributes

ExperimentParameters parameters_
 Struct of several member variables. More...
 
DensityParameters density_param_
 Structure to precalculate and hold parameters for density computations. More...
 
Modus modus_
 Instance of the Modus template parameter. More...
 
std::vector< Particlesensembles_
 Complete particle list, all ensembles in one vector. More...
 
std::unique_ptr< Potentialspotentials_
 An instance of potentials class, that stores parameters of potentials, calculates them and their gradients. More...
 
std::unique_ptr< PauliBlockerpauli_blocker_
 An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and computes phase-space density. More...
 
OutputsList outputs_
 A list of output formaters. More...
 
OutputPtr dilepton_output_
 The Dilepton output. More...
 
OutputPtr photon_output_
 The Photon output. More...
 
std::vector< bool > projectile_target_interact_
 Whether the projectile and the target collided. More...
 
std::vector< FourVectorbeam_momentum_ = {}
 The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen. More...
 
std::vector< std::unique_ptr< ActionFinderInterface > > action_finders_
 The Action finder objects. More...
 
std::unique_ptr< DecayActionsFinderDileptondilepton_finder_
 The Dilepton Action Finder. More...
 
std::unique_ptr< ActionFinderInterfacephoton_finder_
 The (Scatter) Actions Finder for Direct Photons. More...
 
int n_fractional_photons_
 Number of fractional photons produced per single reaction. More...
 
std::unique_ptr< DensityLatticej_QBS_lat_
 4-current for j_QBS lattice output More...
 
std::unique_ptr< DensityLatticejmu_B_lat_
 Baryon density on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_I3_lat_
 Isospin projection density on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_el_lat_
 Electric charge density on the lattice. More...
 
std::unique_ptr< FieldsLatticefields_lat_
 Mean-field A^mu on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_custom_lat_
 Custom density on the lattices. More...
 
DensityType dens_type_lattice_printout_ = DensityType::None
 Type of density for lattice printout. More...
 
std::unique_ptr< RectangularLattice< FourVector > > UB_lat_ = nullptr
 Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-velocity. More...
 
std::unique_ptr< RectangularLattice< FourVector > > UI3_lat_ = nullptr
 Lattices for symmetry potentials (evaluated in the local rest frame) times the isospin flow 4-velocity. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FB_lat_
 Lattices for the electric and magnetic components of the Skyrme or VDF force. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FI3_lat_
 Lattices for the electric and magnetic component of the symmetry force. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > EM_lat_
 Lattices for electric and magnetic field in fm^-2. More...
 
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > > Tmn_
 Lattices of energy-momentum tensors for printout. More...
 
std::unique_ptr< RectangularLattice< FourVector > > old_jmu_auxiliary_
 Auxiliary lattice for values of jmu at a time step t0. More...
 
std::unique_ptr< RectangularLattice< FourVector > > new_jmu_auxiliary_
 Auxiliary lattice for values of jmu at a time step t0 + dt. More...
 
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > four_gradient_auxiliary_
 Auxiliary lattice for calculating the four-gradient of jmu. More...
 
std::unique_ptr< RectangularLattice< FourVector > > old_fields_auxiliary_
 Auxiliary lattice for values of Amu at a time step t0. More...
 
std::unique_ptr< RectangularLattice< FourVector > > new_fields_auxiliary_
 Auxiliary lattice for values of Amu at a time step t0 + dt. More...
 
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > fields_four_gradient_auxiliary_
 Auxiliary lattice for calculating the four-gradient of Amu. More...
 
bool printout_rho_eckart_ = false
 Whether to print the Eckart rest frame density. More...
 
bool printout_tmn_ = false
 Whether to print the energy-momentum tensor. More...
 
bool printout_tmn_landau_ = false
 Whether to print the energy-momentum tensor in Landau frame. More...
 
bool printout_v_landau_ = false
 Whether to print the 4-velocity in Landau frame. More...
 
bool printout_j_QBS_ = false
 Whether to print the Q, B, S 4-currents. More...
 
bool printout_lattice_td_ = false
 Whether to print the thermodynamics quantities evaluated on the lattices. More...
 
bool printout_full_lattice_any_td_ = false
 Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in any format. More...
 
std::unique_ptr< GrandCanThermalizerthermalizer_
 Instance of class used for forced thermalization. More...
 
StringProcessprocess_string_ptr_
 Pointer to the string process class object, which is used to set the random seed for PYTHIA objects in each event. More...
 
const int nevents_ = 0
 Number of events. More...
 
int minimum_nonempty_ensembles_ = 0
 The number of ensembles, in which interactions take place, to be calculated. More...
 
EventCounting event_counting_ = EventCounting::Invalid
 The way in which the number of calculated events is specified. More...
 
int event_ = 0
 Current event. More...
 
int nonempty_ensembles_ = 0
 Number of ensembles containing an interaction. More...
 
int max_events_ = 0
 Maximum number of events to be calculated in order obtain the desired number of non-empty events using the MinimumNonemptyEnsembles option. More...
 
const double end_time_
 simulation time at which the evolution is stopped. More...
 
const double delta_time_startup_
 The clock's timestep size at start up. More...
 
const bool force_decays_
 This indicates whether we force all resonances to decay in the last timestep. More...
 
const bool use_grid_
 This indicates whether to use the grid. More...
 
const ExpansionProperties metric_
 This struct contains information on the metric to be used. More...
 
const bool dileptons_switch_
 This indicates whether dileptons are switched on. More...
 
const bool photons_switch_
 This indicates whether photons are switched on. More...
 
const bool bremsstrahlung_switch_
 This indicates whether bremsstrahlung is switched on. More...
 
const bool IC_output_switch_
 This indicates whether the IC output is enabled. More...
 
const TimeStepMode time_step_mode_
 This indicates whether to use time steps. More...
 
double max_transverse_distance_sqr_ = std::numeric_limits<double>::max()
 Maximal distance at which particles can interact in case of the geometric criterion, squared. More...
 
QuantumNumbers conserved_initial_
 The conserved quantities of the system. More...
 
double initial_mean_field_energy_
 The initial total mean field energy in the system. More...
 
SystemTimePoint time_start_ = SystemClock::now()
 system starting time of the simulation More...
 
DensityType dens_type_ = DensityType::None
 Type of density to be written to collision headers. More...
 
uint64_t interactions_total_ = 0
 Total number of interactions for current timestep. More...
 
uint64_t previous_interactions_total_ = 0
 Total number of interactions for previous timestep. More...
 
uint64_t wall_actions_total_ = 0
 Total number of wall-crossings for current timestep. More...
 
uint64_t previous_wall_actions_total_ = 0
 Total number of wall-crossings for previous timestep. More...
 
uint64_t total_pauli_blocked_ = 0
 Total number of Pauli-blockings for current timestep. More...
 
uint64_t total_hypersurface_crossing_actions_ = 0
 Total number of particles removed from the evolution in hypersurface crossing actions. More...
 
uint64_t discarded_interactions_total_ = 0
 Total number of discarded interactions, because they were invalidated before they could be performed. More...
 
double total_energy_removed_ = 0.0
 Total energy removed from the system in hypersurface crossing actions. More...
 
double total_energy_violated_by_Pythia_ = 0.0
 Total energy violation introduced by Pythia. More...
 
bool kinematic_cuts_for_IC_output_ = false
 This indicates whether kinematic cuts are enabled for the IC output. More...
 
int64_t seed_ = -1
 random seed for the next event. More...
 

Friends

class ExperimentBase
 
std::ostream & operator<< (std::ostream &out, const Experiment &e)
 Writes the initial state for the Experiment to the output stream. More...
 

Additional Inherited Members

- Static Public Member Functions inherited from smash::ExperimentBase
static std::unique_ptr< ExperimentBasecreate (Configuration &config, const std::filesystem::path &output_path)
 Factory method that creates and initializes a new Experiment<Modus>. More...
 

Constructor & Destructor Documentation

◆ Experiment()

template<typename Modus >
smash::Experiment< Modus >::Experiment ( Configuration config,
const std::filesystem::path &  output_path 
)
explicit

Create a new Experiment.

This constructor is only called from the ExperimentBase::create factory method.

Parameters
[in,out]configThe Configuration object contains all initial setup of the experiment. It is forwarded to the constructors of member variables as needed. Note that the object is passed by non-const reference. This is only necessary for bookkeeping: Values are not only read, but actually taken out of the object. Thus, all values that remain were not used.
[in]output_pathThe directory where the output files are written.

Member Function Documentation

◆ run()

template<typename Modus >
void smash::Experiment< Modus >::run
overridevirtual

Runs the experiment.

The constructor does the setup of the experiment. The run function executes the complete experiment.

Implements smash::ExperimentBase.

Definition at line 3073 of file experiment.h.

3073  {
3074  const auto &mainlog = logg[LMain];
3075  for (event_ = 0; !is_finished(); event_++) {
3076  mainlog.info() << "Event " << event_;
3077 
3078  // Sample initial particles, start clock, some printout and book-keeping
3080 
3082 
3083  if (force_decays_) {
3084  do_final_decays();
3085  }
3086 
3087  // Output at event end
3088  final_output();
3089  }
3090 }
const bool force_decays_
This indicates whether we force all resonances to decay in the last timestep.
Definition: experiment.h:595
void run_time_evolution(const double t_end, ParticleList &&add_plist={}, ParticleList &&remove_plist={})
Runs the time evolution of an event with fixed-size time steps or without timesteps,...
Definition: experiment.h:2282
void initialize_new_event()
This is called in the beginning of each event.
Definition: experiment.h:1911
bool is_finished()
Checks wether the desired number events have been calculated.
Definition: experiment.h:3042
int event_
Current event.
Definition: experiment.h:570
void do_final_decays()
Performs the final decays of an event.
Definition: experiment.h:2830
void final_output()
Output at the end of an event.
Definition: experiment.h:2879
const double end_time_
simulation time at which the evolution is stopped.
Definition: experiment.h:582
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
static constexpr int LMain
Definition: experiment.h:89

◆ initialize_new_event()

template<typename Modus >
void smash::Experiment< Modus >::initialize_new_event

This is called in the beginning of each event.

It initializes particles according to selected modus, resets the clock and saves the initial conserved quantities for subsequent sanity checks.

Definition at line 1911 of file experiment.h.

1911  {
1913  logg[LExperiment].info() << "random number seed: " << seed_;
1914  /* Set seed for the next event. It has to be positive, so it can be entered
1915  * in the config.
1916  *
1917  * We have to be careful about the minimal integer, whose absolute value
1918  * cannot be represented. */
1919  int64_t r = random::advance();
1920  while (r == INT64_MIN) {
1921  r = random::advance();
1922  }
1923  seed_ = std::abs(r);
1924  /* Set the random seed used in PYTHIA hadronization
1925  * to be same with the SMASH one.
1926  * In this way we ensure that the results are reproducible
1927  * for every event if one knows SMASH random seed. */
1928  if (process_string_ptr_ != NULL) {
1930  }
1931 
1932  for (Particles &particles : ensembles_) {
1933  particles.reset();
1934  }
1935 
1936  // Sample particles according to the initial conditions
1937  double start_time = -1.0;
1938 
1939  // Sample impact parameter only once per all ensembles
1940  // It should be the same for all ensembles
1941  if (modus_.is_collider()) {
1942  modus_.sample_impact();
1943  logg[LExperiment].info("Impact parameter = ", modus_.impact_parameter(),
1944  " fm");
1945  }
1946  for (Particles &particles : ensembles_) {
1947  start_time = modus_.initial_conditions(&particles, parameters_);
1948  }
1949  /* For box modus make sure that particles are in the box. In principle, after
1950  * a correct initialization they should be, so this is just playing it safe.
1951  */
1952  for (Particles &particles : ensembles_) {
1953  modus_.impose_boundary_conditions(&particles, outputs_);
1954  }
1955  // Reset the simulation clock
1956  double timestep = delta_time_startup_;
1957 
1958  switch (time_step_mode_) {
1959  case TimeStepMode::Fixed:
1960  break;
1961  case TimeStepMode::None:
1962  timestep = end_time_ - start_time;
1963  // Take care of the box modus + timestepless propagation
1964  const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
1965  if (max_dt > 0. && max_dt < timestep) {
1966  timestep = max_dt;
1967  }
1968  break;
1969  }
1970  std::unique_ptr<UniformClock> clock_for_this_event;
1971  if (modus_.is_list() && (timestep < 0.0)) {
1972  throw std::runtime_error(
1973  "Timestep for the given event is negative. \n"
1974  "This might happen if the formation times of the input particles are "
1975  "larger than the specified end time of the simulation.");
1976  }
1977  clock_for_this_event =
1978  std::make_unique<UniformClock>(start_time, timestep, end_time_);
1979  parameters_.labclock = std::move(clock_for_this_event);
1980 
1981  // Reset the output clock
1982  parameters_.outputclock->reset(start_time, true);
1983  // remove time before starting time in case of custom output times.
1984  parameters_.outputclock->remove_times_in_past(start_time);
1985 
1986  logg[LExperiment].debug(
1987  "Lab clock: t_start = ", parameters_.labclock->current_time(),
1988  ", dt = ", parameters_.labclock->timestep_duration());
1989 
1990  /* Save the initial conserved quantum numbers and total momentum in
1991  * the system for conservation checks */
1992  conserved_initial_ = QuantumNumbers(ensembles_);
1993  wall_actions_total_ = 0;
1995  interactions_total_ = 0;
2001  total_energy_removed_ = 0.0;
2003  // Print output headers
2004  logg[LExperiment].info() << hline;
2005  logg[LExperiment].info() << "Time[fm] Ekin[GeV] E_MF[GeV] ETotal[GeV] "
2006  << "ETot/N[GeV] D(ETot/N)[GeV] Scatt&Decays "
2007  << "Particles Comp.Time";
2008  logg[LExperiment].info() << hline;
2009  double E_mean_field = 0.0;
2010  if (potentials_) {
2011  // update_potentials();
2012  // if (parameters.outputclock->current_time() == 0.0 )
2013  // using the lattice is necessary
2014  if ((jmu_B_lat_ != nullptr)) {
2019  parameters_.labclock->timestep_duration(), true);
2020  // Because there was no lattice at t=-Delta_t, the time derivatives
2021  // drho_dt and dj^mu/dt at t=0 are huge, while they shouldn't be; we
2022  // overwrite the time derivative to zero by hand.
2023  for (auto &node : *jmu_B_lat_) {
2024  node.overwrite_drho_dt_to_zero();
2025  node.overwrite_djmu_dt_to_zero();
2026  }
2028  EM_lat_.get(), parameters_);
2029  }
2030  }
2031  initial_mean_field_energy_ = E_mean_field;
2034  parameters_.labclock->current_time(), E_mean_field,
2036 
2037  // Output at event start
2038  for (const auto &output : outputs_) {
2039  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2040  auto event_info = fill_event_info(
2041  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2043  output->at_eventstart(ensembles_[i_ens],
2044  // Pretend each ensemble is an independent event
2045  event_ * parameters_.n_ensembles + i_ens,
2046  event_info);
2047  }
2048  // For thermodynamic output
2049  output->at_eventstart(ensembles_, event_);
2050  // For thermodynamic lattice output
2052  if (printout_rho_eckart_) {
2053  switch (dens_type_lattice_printout_) {
2054  case DensityType::Baryon:
2055  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2057  break;
2059  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2061  break;
2062  case DensityType::None:
2063  break;
2064  default:
2065  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2067  *jmu_custom_lat_);
2068  }
2069  }
2070  if (printout_tmn_) {
2071  output->at_eventstart(event_, ThermodynamicQuantity::Tmn,
2073  }
2074  if (printout_tmn_landau_) {
2075  output->at_eventstart(event_, ThermodynamicQuantity::TmnLandau,
2077  }
2078  if (printout_v_landau_) {
2079  output->at_eventstart(event_, ThermodynamicQuantity::LandauVelocity,
2081  }
2082  if (printout_j_QBS_) {
2083  output->at_eventstart(event_, ThermodynamicQuantity::j_QBS,
2085  }
2086  }
2087  }
2088 
2089  /* In the ColliderModus, if Fermi motion is frozen, assign the beam momenta
2090  * to the nucleons in both the projectile and the target. Every ensemble
2091  * gets the same beam momenta, so no need to create beam_momenta_ vector
2092  * for every ensemble.
2093  */
2094  if (modus_.is_collider() && modus_.fermi_motion() == FermiMotion::Frozen) {
2095  for (ParticleData &particle : ensembles_[0]) {
2096  const double m = particle.effective_mass();
2097  double v_beam = 0.0;
2098  if (particle.belongs_to() == BelongsTo::Projectile) {
2099  v_beam = modus_.velocity_projectile();
2100  } else if (particle.belongs_to() == BelongsTo::Target) {
2101  v_beam = modus_.velocity_target();
2102  }
2103  const double gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
2104  beam_momentum_.emplace_back(
2105  FourVector(gamma * m, 0.0, 0.0, gamma * v_beam * m));
2106  } // loop over particles
2107  }
2108 }
double initial_mean_field_energy_
The initial total mean field energy in the system.
Definition: experiment.h:637
bool printout_tmn_
Whether to print the energy-momentum tensor.
Definition: experiment.h:508
QuantumNumbers conserved_initial_
The conserved quantities of the system.
Definition: experiment.h:631
DensityParameters density_param_
Structure to precalculate and hold parameters for density computations.
Definition: experiment.h:371
double total_energy_removed_
Total energy removed from the system in hypersurface crossing actions.
Definition: experiment.h:690
DensityType dens_type_lattice_printout_
Type of density for lattice printout.
Definition: experiment.h:456
double max_transverse_distance_sqr_
Maximal distance at which particles can interact in case of the geometric criterion,...
Definition: experiment.h:622
bool printout_j_QBS_
Whether to print the Q, B, S 4-currents.
Definition: experiment.h:517
const TimeStepMode time_step_mode_
This indicates whether to use time steps.
Definition: experiment.h:616
std::unique_ptr< DensityLattice > jmu_custom_lat_
Custom density on the lattices.
Definition: experiment.h:453
bool printout_full_lattice_any_td_
Whether to print the thermodynamics quantities evaluated on the lattices, point by point,...
Definition: experiment.h:524
std::unique_ptr< DensityLattice > jmu_B_lat_
Baryon density on the lattice.
Definition: experiment.h:435
std::vector< FourVector > beam_momentum_
The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen.
Definition: experiment.h:417
std::unique_ptr< RectangularLattice< FourVector > > new_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0 + dt.
Definition: experiment.h:491
const double delta_time_startup_
The clock's timestep size at start up.
Definition: experiment.h:589
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > > Tmn_
Lattices of energy-momentum tensors for printout.
Definition: experiment.h:486
std::vector< Particles > ensembles_
Complete particle list, all ensembles in one vector.
Definition: experiment.h:380
SystemTimePoint time_start_
system starting time of the simulation
Definition: experiment.h:640
uint64_t previous_wall_actions_total_
Total number of wall-crossings for previous timestep.
Definition: experiment.h:667
OutputsList outputs_
A list of output formaters.
Definition: experiment.h:398
bool printout_rho_eckart_
Whether to print the Eckart rest frame density.
Definition: experiment.h:505
bool printout_v_landau_
Whether to print the 4-velocity in Landau frame.
Definition: experiment.h:514
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > EM_lat_
Lattices for electric and magnetic field in fm^-2.
Definition: experiment.h:483
Modus modus_
Instance of the Modus template parameter.
Definition: experiment.h:377
std::unique_ptr< RectangularLattice< FourVector > > old_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0.
Definition: experiment.h:489
std::unique_ptr< DensityLattice > j_QBS_lat_
4-current for j_QBS lattice output
Definition: experiment.h:432
bool printout_tmn_landau_
Whether to print the energy-momentum tensor in Landau frame.
Definition: experiment.h:511
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of jmu.
Definition: experiment.h:494
StringProcess * process_string_ptr_
Pointer to the string process class object, which is used to set the random seed for PYTHIA objects i...
Definition: experiment.h:533
ExperimentParameters parameters_
Struct of several member variables.
Definition: experiment.h:368
std::unique_ptr< Potentials > potentials_
An instance of potentials class, that stores parameters of potentials, calculates them and their grad...
Definition: experiment.h:386
uint64_t wall_actions_total_
Total number of wall-crossings for current timestep.
Definition: experiment.h:661
uint64_t total_hypersurface_crossing_actions_
Total number of particles removed from the evolution in hypersurface crossing actions.
Definition: experiment.h:679
uint64_t interactions_total_
Total number of interactions for current timestep.
Definition: experiment.h:649
std::unique_ptr< DensityLattice > jmu_I3_lat_
Isospin projection density on the lattice.
Definition: experiment.h:438
uint64_t total_pauli_blocked_
Total number of Pauli-blockings for current timestep.
Definition: experiment.h:673
std::vector< bool > projectile_target_interact_
Whether the projectile and the target collided.
Definition: experiment.h:410
int64_t seed_
random seed for the next event.
Definition: experiment.h:701
uint64_t previous_interactions_total_
Total number of interactions for previous timestep.
Definition: experiment.h:655
uint64_t discarded_interactions_total_
Total number of discarded interactions, because they were invalidated before they could be performed.
Definition: experiment.h:685
bool kinematic_cuts_for_IC_output_
This indicates whether kinematic cuts are enabled for the IC output.
Definition: experiment.h:698
double total_energy_violated_by_Pythia_
Total energy violation introduced by Pythia.
Definition: experiment.h:695
void init_pythia_hadron_rndm()
Set PYTHIA random seeds to be desired values.
@ Frozen
Use fermi motion without potentials.
@ Fixed
Use fixed time step.
@ None
Don't use time steps; propagate from action to action.
@ EckartDensity
Density in the Eckart frame.
@ Tmn
Energy-momentum tensor in lab frame.
@ LandauVelocity
Velocity of the Landau rest frame.
@ j_QBS
Electric (Q), baryonic (B) and strange (S) currents.
@ TmnLandau
Energy-momentum tensor in Landau rest frame.
Engine::result_type advance()
Advance the engine's state and return the generated value.
Definition: random.h:78
void set_seed(T &&seed)
Sets the seed of the random number engine.
Definition: random.h:71
EventInfo fill_event_info(const std::vector< Particles > &ensembles, double E_mean_field, double modus_impact_parameter, const ExperimentParameters &parameters, bool projectile_target_interact, bool kinematic_cut_for_SMASH_IC)
Generate the EventInfo object which is passed to outputs_.
Definition: experiment.cc:617
std::string format_measurements(const std::vector< Particles > &ensembles, uint64_t scatterings_this_interval, const QuantumNumbers &conserved_initial, SystemTimePoint time_start, double time, double E_mean_field, double E_mean_field_initial)
Generate a string which will be printed to the screen when SMASH is running.
Definition: experiment.cc:349
void update_lattice(RectangularLattice< T > *lat, const LatticeUpdate update, const DensityType dens_type, const DensityParameters &par, const std::vector< Particles > &ensembles, const bool compute_gradient)
Updates the contents on the lattice.
Definition: density.h:542
double calculate_mean_field_energy(const Potentials &potentials, RectangularLattice< smash::DensityOnLattice > &jmu_B_lat, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *em_lattice, const ExperimentParameters &parameters)
Calculate the total mean field energy of the system; this will be printed to the screen when SMASH is...
Definition: experiment.cc:397
static constexpr int LExperiment
const std::string hline(113, '-')
String representing a horizontal line.
int n_ensembles
Number of parallel ensembles.
std::unique_ptr< Clock > outputclock
Output clock to keep track of the next output time.
std::unique_ptr< Clock > labclock
System clock (for simulation time keeping in the computational frame)

◆ run_time_evolution()

template<typename Modus >
void smash::Experiment< Modus >::run_time_evolution ( const double  t_end,
ParticleList &&  add_plist = {},
ParticleList &&  remove_plist = {} 
)

Runs the time evolution of an event with fixed-size time steps or without timesteps, from action to actions.

Within one timestep (fixed) evolution from action to action is invoked.

Parameters
[in]t_endTime until run_time_evolution is run, in SMASH this is the configured end_time, but it might differ if SMASH is used as an external library
[in]add_plistA by-default empty particle list which is added to the current particle content of the system
[in]remove_plistA by-default empty particle list which is removed from the current particle content of the system
Note
This function is meant to take over ownership of the to-be-added/removed particle lists and that's why these are passed by rvalue reference.

Definition at line 2282 of file experiment.h.

2284  {
2285  if (!add_plist.empty() || !remove_plist.empty()) {
2286  if (ensembles_.size() > 1) {
2287  throw std::runtime_error(
2288  "Adding or removing particles from SMASH is only possible when one "
2289  "ensemble is used.");
2290  }
2291  const double action_time = parameters_.labclock->current_time();
2292  /* Use two if statements. The first one is to check if the particles are
2293  * valid. Since this might remove all particles, a second if statement is
2294  * needed to avoid executing the action in that case.*/
2295  if (!add_plist.empty()) {
2297  }
2298  if (!add_plist.empty()) {
2299  // Create and perform action to add particle(s)
2300  auto action_add_particles = std::make_unique<FreeforallAction>(
2301  ParticleList{}, add_plist, action_time);
2302  perform_action(*action_add_particles, 0);
2303  }
2304  // Also here 2 if statements are needed as above.
2305  if (!remove_plist.empty()) {
2306  validate_and_adjust_particle_list(remove_plist);
2307  }
2308  if (!remove_plist.empty()) {
2309  ParticleList found_particles_to_remove;
2310  for (const auto &particle_to_remove : remove_plist) {
2311  const auto iterator_to_particle_to_be_removed_in_ensemble =
2312  std::find_if(
2313  ensembles_[0].begin(), ensembles_[0].end(),
2314  [&particle_to_remove, &action_time](const ParticleData &p) {
2316  particle_to_remove, p, action_time);
2317  });
2318  if (iterator_to_particle_to_be_removed_in_ensemble !=
2319  ensembles_[0].end())
2320  found_particles_to_remove.push_back(
2321  *iterator_to_particle_to_be_removed_in_ensemble);
2322  }
2323  // Sort the particles found to be removed according to their id and look
2324  // for duplicates (sorting is needed to call std::adjacent_find).
2325  std::sort(found_particles_to_remove.begin(),
2326  found_particles_to_remove.end(),
2327  [](const ParticleData &p1, const ParticleData &p2) {
2328  return p1.id() < p2.id();
2329  });
2330  const auto iterator_to_first_duplicate = std::adjacent_find(
2331  found_particles_to_remove.begin(), found_particles_to_remove.end(),
2332  [](const ParticleData &p1, const ParticleData &p2) {
2333  return p1.id() == p2.id();
2334  });
2335  if (iterator_to_first_duplicate != found_particles_to_remove.end()) {
2336  logg[LExperiment].error() << "The same particle has been asked to be "
2337  "removed multiple times:\n"
2338  << *iterator_to_first_duplicate;
2339  throw std::logic_error("Particle cannot be removed twice!");
2340  }
2341  if (auto delta = remove_plist.size() - found_particles_to_remove.size();
2342  delta > 0) {
2343  logg[LExperiment].warn(
2344  "When trying to remove particle(s) at the beginning ",
2345  "of the system evolution,\n", delta,
2346  " particle(s) could not be found and will be ignored.");
2347  }
2348  if (!found_particles_to_remove.empty()) {
2349  [[maybe_unused]] const auto number_particles_before_removal =
2350  ensembles_[0].size();
2351  // Create and perform action to remove particles
2352  auto action_remove_particles = std::make_unique<FreeforallAction>(
2353  found_particles_to_remove, ParticleList{}, action_time);
2354  perform_action(*action_remove_particles, 0);
2355 
2356  assert(number_particles_before_removal -
2357  found_particles_to_remove.size() ==
2358  ensembles_[0].size());
2359  }
2360  }
2361  }
2362 
2363  if (t_end > end_time_) {
2364  logg[LExperiment].fatal()
2365  << "Evolution asked to be run until " << t_end << " > " << end_time_
2366  << " and this cannot be done (because of how the clock works).";
2367  throw std::logic_error(
2368  "Experiment cannot evolve the system beyond End_Time.");
2369  }
2370  while (*(parameters_.labclock) < t_end) {
2371  const double dt = parameters_.labclock->timestep_duration();
2372  logg[LExperiment].debug("Timestepless propagation for next ", dt, " fm.");
2373 
2374  // Perform forced thermalization if required
2375  if (thermalizer_ &&
2376  thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
2377  const bool ignore_cells_under_treshold = true;
2378  // Thermodynamics in thermalizer is computed from all ensembles,
2379  // but thermalization actions act on each ensemble independently
2380  thermalizer_->update_thermalizer_lattice(ensembles_, density_param_,
2381  ignore_cells_under_treshold);
2382  const double current_t = parameters_.labclock->current_time();
2383  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2384  thermalizer_->thermalize(ensembles_[i_ens], current_t,
2386  ThermalizationAction th_act(*thermalizer_, current_t);
2387  if (th_act.any_particles_thermalized()) {
2388  perform_action(th_act, i_ens);
2389  }
2390  }
2391  }
2392 
2393  std::vector<Actions> actions(parameters_.n_ensembles);
2394  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2395  actions[i_ens].clear();
2396  if (ensembles_[i_ens].size() > 0 && action_finders_.size() > 0) {
2397  /* (1.a) Create grid. */
2398  const double min_cell_length = compute_min_cell_length(dt);
2399  logg[LExperiment].debug("Creating grid with minimal cell length ",
2400  min_cell_length);
2401  /* For the hyper-surface-crossing actions also unformed particles are
2402  * searched and therefore needed on the grid. */
2403  const bool include_unformed_particles = IC_output_switch_;
2404  const auto &grid =
2405  use_grid_ ? modus_.create_grid(ensembles_[i_ens], min_cell_length,
2406  dt, parameters_.coll_crit,
2407  include_unformed_particles)
2408  : modus_.create_grid(ensembles_[i_ens], min_cell_length,
2409  dt, parameters_.coll_crit,
2410  include_unformed_particles,
2412 
2413  const double gcell_vol = grid.cell_volume();
2414  /* (1.b) Iterate over cells and find actions. */
2415  grid.iterate_cells(
2416  [&](const ParticleList &search_list) {
2417  for (const auto &finder : action_finders_) {
2418  actions[i_ens].insert(finder->find_actions_in_cell(
2419  search_list, dt, gcell_vol, beam_momentum_));
2420  }
2421  },
2422  [&](const ParticleList &search_list,
2423  const ParticleList &neighbors_list) {
2424  for (const auto &finder : action_finders_) {
2425  actions[i_ens].insert(finder->find_actions_with_neighbors(
2426  search_list, neighbors_list, dt, beam_momentum_));
2427  }
2428  });
2429  }
2430  }
2431 
2432  /* \todo (optimizations) Adapt timestep size here */
2433 
2434  /* (2) Propagate from action to action until next output or timestep end */
2435  const double end_timestep_time = parameters_.labclock->next_time();
2436  while (next_output_time() < end_timestep_time) {
2437  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2438  run_time_evolution_timestepless(actions[i_ens], i_ens,
2439  next_output_time());
2440  }
2441  ++(*parameters_.outputclock);
2442 
2444  }
2445  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2446  run_time_evolution_timestepless(actions[i_ens], i_ens, end_timestep_time);
2447  }
2448 
2449  /* (3) Update potentials (if computed on the lattice) and
2450  * compute new momenta according to equations of motion */
2451  if (potentials_) {
2453  update_momenta(ensembles_, parameters_.labclock->timestep_duration(),
2454  *potentials_, FB_lat_.get(), FI3_lat_.get(), EM_lat_.get(),
2455  jmu_B_lat_.get());
2456  }
2457 
2458  /* (4) Expand universe if non-minkowskian metric; updates
2459  * positions and momenta according to the selected expansion */
2461  for (Particles &particles : ensembles_) {
2462  expand_space_time(&particles, parameters_, metric_);
2463  }
2464  }
2465 
2466  ++(*parameters_.labclock);
2467 
2468  /* (5) Check conservation laws.
2469  *
2470  * Check conservation of conserved quantities if potentials and string
2471  * fragmentation are off. If potentials are on then momentum is conserved
2472  * only in average. If string fragmentation is on, then energy and
2473  * momentum are only very roughly conserved in high-energy collisions. */
2476  std::string err_msg = conserved_initial_.report_deviations(ensembles_);
2477  if (!err_msg.empty()) {
2478  logg[LExperiment].error() << err_msg;
2479  throw std::runtime_error("Violation of conserved quantities!");
2480  }
2481  }
2482  }
2483 
2484  if (pauli_blocker_) {
2485  logg[LExperiment].info(
2486  "Interactions: Pauli-blocked/performed = ", total_pauli_blocked_, "/",
2488  }
2489 }
const ExpansionProperties metric_
This struct contains information on the metric to be used.
Definition: experiment.h:601
std::vector< std::unique_ptr< ActionFinderInterface > > action_finders_
The Action finder objects.
Definition: experiment.h:420
double next_output_time() const
Shortcut for next output time.
Definition: experiment.h:346
std::unique_ptr< GrandCanThermalizer > thermalizer_
Instance of class used for forced thermalization.
Definition: experiment.h:527
void run_time_evolution_timestepless(Actions &actions, int i_ensemble, const double end_time_propagation)
Performs all the propagations and actions during a certain time interval neglecting the influence of ...
Definition: experiment.h:2518
void intermediate_output()
Intermediate output during an event.
Definition: experiment.h:2582
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FI3_lat_
Lattices for the electric and magnetic component of the symmetry force.
Definition: experiment.h:479
double compute_min_cell_length(double dt) const
Calculate the minimal size for the grid cells such that the ScatterActionsFinder will find all collis...
Definition: experiment.h:338
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FB_lat_
Lattices for the electric and magnetic components of the Skyrme or VDF force.
Definition: experiment.h:475
std::unique_ptr< PauliBlocker > pauli_blocker_
An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and c...
Definition: experiment.h:392
bool perform_action(Action &action, int i_ensemble, bool include_pauli_blocking=true)
Perform the given action.
const bool use_grid_
This indicates whether to use the grid.
Definition: experiment.h:598
void update_potentials()
Recompute potentials on lattices if necessary.
Definition: experiment.h:2730
const bool IC_output_switch_
This indicates whether the IC output is enabled.
Definition: experiment.h:613
std::string report_deviations(const std::vector< Particles > &ensembles) const
Checks if the current particle list has still the same values and reports about differences.
constexpr int p
Proton.
void update_momenta(std::vector< Particles > &particles, double dt, const Potentials &pot, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *FB_lat, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *FI3_lat, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *EM_lat, DensityLattice *jB_lat)
Updates the momenta of all particles at the current time step according to the equations of motion:
Definition: propagation.cc:111
void expand_space_time(Particles *particles, const ExperimentParameters &parameters, const ExpansionProperties &metric)
Modifies positions and momentum of all particles to account for space-time deformation.
Definition: propagation.cc:86
bool are_particles_identical_at_given_time(const ParticleData &p1, const ParticleData &p2, double time)
Utility function to compare two ParticleData instances with respect to their PDG code,...
void validate_and_adjust_particle_list(ParticleList &particle_list)
Validate a particle list adjusting each particle to be a valid SMASH particle.
Definition: experiment.cc:639
@ Largest
Make cells as large as possible.
ExpansionMode mode_
Type of metric used.
Definition: propagation.h:28
const CollisionCriterion coll_crit
Employed collision criterion.
bool strings_switch
This indicates whether string fragmentation is switched on.
int testparticles
Number of test-particles.

◆ do_final_decays()

template<typename Modus >
void smash::Experiment< Modus >::do_final_decays

Performs the final decays of an event.

Exceptions
runtime_errorif found actions cannot be performed

Definition at line 2830 of file experiment.h.

2830  {
2831  /* At end of time evolution: Force all resonances to decay. In order to handle
2832  * decay chains, we need to loop until no further actions occur. */
2833  bool actions_performed, decays_found;
2834  uint64_t interactions_old;
2835  do {
2836  decays_found = false;
2837  interactions_old = interactions_total_;
2838  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2839  Actions actions;
2840 
2841  // Dileptons: shining of remaining resonances
2842  if (dilepton_finder_ != nullptr) {
2843  for (const auto &output : outputs_) {
2844  dilepton_finder_->shine_final(ensembles_[i_ens], output.get(), true);
2845  }
2846  }
2847  // Find actions.
2848  for (const auto &finder : action_finders_) {
2849  auto found_actions = finder->find_final_actions(ensembles_[i_ens]);
2850  if (!found_actions.empty()) {
2851  actions.insert(std::move(found_actions));
2852  decays_found = true;
2853  }
2854  }
2855  // Perform actions.
2856  while (!actions.is_empty()) {
2857  perform_action(*actions.pop(), i_ens, false);
2858  }
2859  }
2860  actions_performed = interactions_total_ > interactions_old;
2861  // Throw an error if actions were found but not performed
2862  if (decays_found && !actions_performed) {
2863  throw std::runtime_error("Final decays were found but not performed.");
2864  }
2865  // loop until no more decays occur
2866  } while (actions_performed);
2867 
2868  // Dileptons: shining of stable particles at the end
2869  if (dilepton_finder_ != nullptr) {
2870  for (const auto &output : outputs_) {
2871  for (Particles &particles : ensembles_) {
2872  dilepton_finder_->shine_final(particles, output.get(), false);
2873  }
2874  }
2875  }
2876 }
std::unique_ptr< DecayActionsFinderDilepton > dilepton_finder_
The Dilepton Action Finder.
Definition: experiment.h:423

◆ final_output()

template<typename Modus >
void smash::Experiment< Modus >::final_output

Output at the end of an event.

Definition at line 2879 of file experiment.h.

2879  {
2880  /* make sure the experiment actually ran (note: we should compare this
2881  * to the start time, but we don't know that. Therefore, we check that
2882  * the time is positive, which should heuristically be the same). */
2883  double E_mean_field = 0.0;
2884  if (likely(parameters_.labclock > 0)) {
2885  const uint64_t wall_actions_this_interval =
2887  const uint64_t interactions_this_interval = interactions_total_ -
2889  wall_actions_this_interval;
2890  if (potentials_) {
2891  // using the lattice is necessary
2892  if ((jmu_B_lat_ != nullptr)) {
2894  EM_lat_.get(), parameters_);
2895  }
2896  }
2897  if (std::abs(parameters_.labclock->current_time() - end_time_) >
2898  really_small) {
2899  logg[LExperiment].warn()
2900  << "SMASH not propagated until configured end time. Current time = "
2901  << parameters_.labclock->current_time()
2902  << "fm. End time = " << end_time_ << "fm.";
2903  } else {
2905  ensembles_, interactions_this_interval, conserved_initial_,
2907  }
2908  int total_particles = 0;
2909  for (const Particles &particles : ensembles_) {
2910  total_particles += particles.size();
2911  }
2912  if (IC_output_switch_ && (total_particles == 0)) {
2913  const double initial_system_energy_plus_Pythia_violations =
2915  const double fraction_of_total_system_energy_removed =
2916  initial_system_energy_plus_Pythia_violations / total_energy_removed_;
2917  // Verify there is no more energy in the system if all particles were
2918  // removed when crossing the hypersurface
2919  if (std::fabs(fraction_of_total_system_energy_removed - 1.) >
2920  really_small) {
2921  throw std::runtime_error(
2922  "There is remaining energy in the system although all particles "
2923  "were removed.\n"
2924  "E_remain = " +
2925  std::to_string((initial_system_energy_plus_Pythia_violations -
2927  " [GeV]");
2928  } else {
2929  logg[LExperiment].info() << hline;
2930  logg[LExperiment].info()
2931  << "Time real: " << SystemClock::now() - time_start_;
2932  logg[LExperiment].info()
2933  << "Interactions before reaching hypersurface: "
2936  logg[LExperiment].info()
2937  << "Total number of particles removed on hypersurface: "
2939  }
2940  } else {
2941  const double precent_discarded =
2943  ? static_cast<double>(discarded_interactions_total_) * 100.0 /
2945  : 0.0;
2946  std::stringstream msg_discarded;
2947  msg_discarded
2948  << "Discarded interaction number: " << discarded_interactions_total_
2949  << " (" << precent_discarded
2950  << "% of the total interaction number including wall crossings)";
2951 
2952  logg[LExperiment].info() << hline;
2953  logg[LExperiment].info()
2954  << "Time real: " << SystemClock::now() - time_start_;
2955  logg[LExperiment].debug() << msg_discarded.str();
2956 
2958  precent_discarded > 1.0) {
2959  // The choosen threshold of 1% is a heuristical value
2960  logg[LExperiment].warn()
2961  << msg_discarded.str()
2962  << "\nThe number of discarded interactions is large, which means "
2963  "the assumption for the stochastic criterion of\n1 interaction "
2964  "per particle per timestep is probably violated. Consider "
2965  "reducing the timestep size.";
2966  }
2967 
2968  logg[LExperiment].info() << "Final interaction number: "
2970  }
2971 
2972  // Check if there are unformed particles
2973  int unformed_particles_count = 0;
2974  for (const Particles &particles : ensembles_) {
2975  for (const ParticleData &particle : particles) {
2976  if (particle.formation_time() > end_time_) {
2977  unformed_particles_count++;
2978  }
2979  }
2980  }
2981  if (unformed_particles_count > 0) {
2982  logg[LExperiment].warn(
2983  "End time might be too small. ", unformed_particles_count,
2984  " unformed particles were found at the end of the evolution.");
2985  }
2986  }
2987 
2988  // Keep track of how many ensembles had interactions
2990 
2991  for (const auto &output : outputs_) {
2992  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2993  auto event_info = fill_event_info(
2994  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2996  output->at_eventend(ensembles_[i_ens],
2997  // Pretend each ensemble is an independent event
2998  event_ * parameters_.n_ensembles + i_ens, event_info);
2999  }
3000  // For thermodynamic output
3001  output->at_eventend(ensembles_, event_);
3002 
3003  // For thermodynamic lattice output
3004  if (printout_rho_eckart_) {
3006  output->at_eventend(event_, ThermodynamicQuantity::EckartDensity,
3008  output->at_eventend(ThermodynamicQuantity::EckartDensity);
3009  }
3010  }
3011  if (printout_tmn_) {
3012  output->at_eventend(event_, ThermodynamicQuantity::Tmn,
3014  output->at_eventend(ThermodynamicQuantity::Tmn);
3015  }
3016  if (printout_tmn_landau_) {
3017  output->at_eventend(event_, ThermodynamicQuantity::TmnLandau,
3019  output->at_eventend(ThermodynamicQuantity::TmnLandau);
3020  }
3021  if (printout_v_landau_) {
3022  output->at_eventend(event_, ThermodynamicQuantity::LandauVelocity,
3024  output->at_eventend(ThermodynamicQuantity::LandauVelocity);
3025  }
3026  if (printout_j_QBS_) {
3027  output->at_eventend(ThermodynamicQuantity::j_QBS);
3028  }
3029  }
3030 }
void count_nonempty_ensembles()
Counts the number of ensembles in wich interactions took place at the end of an event.
Definition: experiment.h:3033
double x0() const
Definition: fourvector.h:313
FourVector momentum() const
@ Stochastic
Stochastic Criteiron.
#define likely(x)
Tell the branch predictor that this expression is likely true.
Definition: macros.h:14
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37

◆ first_ensemble()

template<typename Modus >
Particles* smash::Experiment< Modus >::first_ensemble ( )
inline

Provides external access to SMASH particles.

This is helpful if SMASH is used as a 3rd-party library.

Definition at line 257 of file experiment.h.

257 { return &ensembles_[0]; }

◆ all_ensembles()

template<typename Modus >
std::vector<Particles>* smash::Experiment< Modus >::all_ensembles ( )
inline

Getter for all ensembles.

Definition at line 259 of file experiment.h.

259 { return &ensembles_; }

◆ modus()

template<typename Modus >
Modus* smash::Experiment< Modus >::modus ( )
inline

Provides external access to SMASH calculation modus.

This is helpful if SMASH is used as a 3rd-party library.

Definition at line 265 of file experiment.h.

265 { return &modus_; }

◆ increase_event_number()

template<typename Modus >
void smash::Experiment< Modus >::increase_event_number

Increases the event number by one.

This function is helpful if SMASH is used as a 3rd-party library.

Definition at line 3068 of file experiment.h.

3068  {
3069  event_++;
3070 }

◆ perform_action()

template<typename Modus >
bool smash::Experiment< Modus >::perform_action ( Action action,
int  i_ensemble,
bool  include_pauli_blocking = true 
)
private

Perform the given action.

Parameters
[in]actionThe action to perform
[in]i_ensembleindex of ensemble in which action is performed
[in]include_pauli_blockingwheter to take Pauli blocking into account. Skipping Pauli blocking is useful for example for final decays.
Returns
False if the action is rejected either due to invalidity or Pauli-blocking, or true if it's accepted and performed.

◆ create_output()

template<typename Modus >
void smash::Experiment< Modus >::create_output ( const std::string &  format,
const std::string &  content,
const std::filesystem::path &  output_path,
const OutputParameters par 
)
private

Create a list of output files.

Parameters
[in]formatFormat of the output file (e.g. Root, Oscar, Vtk)
[in]contentContent of the output (e.g. particles, collisions)
[in]output_pathPath of the output file
[in]parOutput options.(e.g. Extended)

Definition at line 720 of file experiment.h.

723  {
724  logg[LExperiment].info() << "Adding output " << content << " of format "
725  << format << std::endl;
726 
727  if (format == "VTK" && content == "Particles") {
728  outputs_.emplace_back(
729  std::make_unique<VtkOutput>(output_path, content, out_par));
730  } else if (format == "Root") {
731 #ifdef SMASH_USE_ROOT
732  if (content == "Initial_Conditions") {
733  outputs_.emplace_back(
734  std::make_unique<RootOutput>(output_path, "SMASH_IC", out_par));
735  } else {
736  outputs_.emplace_back(
737  std::make_unique<RootOutput>(output_path, content, out_par));
738  }
739 #else
740  logg[LExperiment].error(
741  "Root output requested, but Root support not compiled in");
742 #endif
743  } else if (format == "Binary") {
744  if (content == "Collisions" || content == "Dileptons" ||
745  content == "Photons") {
746  outputs_.emplace_back(std::make_unique<BinaryOutputCollisions>(
747  output_path, content, out_par));
748  } else if (content == "Particles") {
749  outputs_.emplace_back(std::make_unique<BinaryOutputParticles>(
750  output_path, content, out_par));
751  } else if (content == "Initial_Conditions") {
752  outputs_.emplace_back(std::make_unique<BinaryOutputInitialConditions>(
753  output_path, content, out_par));
754  }
755  } else if (format == "Oscar1999" || format == "Oscar2013") {
756  outputs_.emplace_back(
757  create_oscar_output(format, content, output_path, out_par));
758  } else if (content == "Thermodynamics" && format == "ASCII") {
759  outputs_.emplace_back(
760  std::make_unique<ThermodynamicOutput>(output_path, content, out_par));
761  } else if (content == "Thermodynamics" &&
762  (format == "Lattice_ASCII" || format == "Lattice_Binary")) {
764  outputs_.emplace_back(std::make_unique<ThermodynamicLatticeOutput>(
765  output_path, content, out_par, format == "Lattice_ASCII",
766  format == "Lattice_Binary"));
767  } else if (content == "Thermodynamics" && format == "VTK") {
768  printout_lattice_td_ = true;
769  outputs_.emplace_back(
770  std::make_unique<VtkOutput>(output_path, content, out_par));
771  } else if (content == "Initial_Conditions" && format == "ASCII") {
772  outputs_.emplace_back(
773  std::make_unique<ICOutput>(output_path, "SMASH_IC", out_par));
774  } else if ((format == "HepMC") || (format == "HepMC_asciiv3") ||
775  (format == "HepMC_treeroot")) {
776 #ifdef SMASH_USE_HEPMC
777  if (content == "Particles") {
778  if ((format == "HepMC") || (format == "HepMC_asciiv3")) {
779  outputs_.emplace_back(std::make_unique<HepMcOutput>(
780  output_path, "SMASH_HepMC_particles", false, "asciiv3"));
781  } else if (format == "HepMC_treeroot") {
782 #ifdef SMASH_USE_HEPMC_ROOTIO
783  outputs_.emplace_back(std::make_unique<HepMcOutput>(
784  output_path, "SMASH_HepMC_particles", false, "root"));
785 #else
786  logg[LExperiment].error(
787  "Requested HepMC_treeroot output not available, "
788  "ROOT or HepMC3 ROOTIO missing or not found by cmake.");
789 #endif
790  }
791  } else if (content == "Collisions") {
792  if ((format == "HepMC") || (format == "HepMC_asciiv3")) {
793  outputs_.emplace_back(std::make_unique<HepMcOutput>(
794  output_path, "SMASH_HepMC_collisions", true, "asciiv3"));
795  } else if (format == "HepMC_treeroot") {
796 #ifdef SMASH_USE_HEPMC_ROOTIO
797  outputs_.emplace_back(std::make_unique<HepMcOutput>(
798  output_path, "SMASH_HepMC_collisions", true, "root"));
799 #else
800  logg[LExperiment].error(
801  "Requested HepMC_treeroot output not available, "
802  "ROOT or HepMC3 ROOTIO missing or not found by cmake.");
803 #endif
804  }
805  } else {
806  logg[LExperiment].error(
807  "HepMC only available for Particles and "
808  "Collisions content. Requested for " +
809  content + ".");
810  }
811 #else
812  logg[LExperiment].error(
813  "HepMC output requested, but HepMC support not compiled in");
814 #endif
815  } else if (content == "Coulomb" && format == "VTK") {
816  outputs_.emplace_back(
817  std::make_unique<VtkOutput>(output_path, "Fields", out_par));
818  } else if (content == "Rivet") {
819 #ifdef SMASH_USE_RIVET
820  // flag to ensure that the Rivet format has not been already assigned
821  static bool rivet_format_already_selected = false;
822  // if the next check is true, then we are trying to assign the format twice
823  if (rivet_format_already_selected) {
824  logg[LExperiment].warn(
825  "Rivet output format can only be one, either YODA or YODA-full. "
826  "Only your first valid choice will be used.");
827  return;
828  }
829  if (format == "YODA") {
830  outputs_.emplace_back(std::make_unique<RivetOutput>(
831  output_path, "SMASH_Rivet", false, out_par.rivet_parameters));
832  rivet_format_already_selected = true;
833  } else if (format == "YODA-full") {
834  outputs_.emplace_back(std::make_unique<RivetOutput>(
835  output_path, "SMASH_Rivet_full", true, out_par.rivet_parameters));
836  rivet_format_already_selected = true;
837  } else {
838  logg[LExperiment].error("Rivet format " + format +
839  "not one of YODA or YODA-full");
840  }
841 #else
842  logg[LExperiment].error(
843  "Rivet output requested, but Rivet support not compiled in");
844 #endif
845  } else {
846  logg[LExperiment].error()
847  << "Unknown combination of format (" << format << ") and content ("
848  << content << "). Fix the config.";
849  }
850 }
bool printout_lattice_td_
Whether to print the thermodynamics quantities evaluated on the lattices.
Definition: experiment.h:520
FormattingHelper< T > format(const T &value, const char *unit, int width=-1, int precision=-1)
Acts as a stream modifier for std::ostream to output an object with an optional suffix string and wit...
Definition: logging.h:217
std::unique_ptr< OutputInterface > create_oscar_output(const std::string &format, const std::string &content, const std::filesystem::path &path, const OutputParameters &out_par)
Definition: oscaroutput.cc:811

◆ propagate_and_shine()

template<typename Modus >
void smash::Experiment< Modus >::propagate_and_shine ( double  to_time,
Particles particles 
)
private

Propagate all particles until time to_time without any interactions and shine dileptons.

Parameters
[in]to_timeTime at the end of propagation [fm]
[in,out]particlesParticles to be propagated

Definition at line 2492 of file experiment.h.

2493  {
2494  const double dt =
2495  propagate_straight_line(&particles, to_time, beam_momentum_);
2496  if (dilepton_finder_ != nullptr) {
2497  for (const auto &output : outputs_) {
2498  dilepton_finder_->shine(particles, output.get(), dt);
2499  }
2500  }
2501 }
double propagate_straight_line(Particles *particles, double to_time, const std::vector< FourVector > &beam_momentum)
Propagates the positions of all particles on a straight line to a given moment.
Definition: propagation.cc:44

◆ run_time_evolution_timestepless()

template<typename Modus >
void smash::Experiment< Modus >::run_time_evolution_timestepless ( Actions actions,
int  i_ensemble,
const double  end_time_propagation 
)
private

Performs all the propagations and actions during a certain time interval neglecting the influence of the potentials.

This function is called in either the time stepless cases or the cases with time steps.

Parameters
[in,out]actionsActions occur during a certain time interval. They provide the ending times of the propagations and are updated during the time interval.
[in]i_ensembleindex of ensemble to be evolved
[in]end_time_propagationtime until propagation should be performed

Definition at line 2518 of file experiment.h.

2519  {
2520  Particles &particles = ensembles_[i_ensemble];
2521  logg[LExperiment].debug(
2522  "Timestepless propagation: ", "Actions size = ", actions.size(),
2523  ", end time = ", end_time_propagation);
2524 
2525  // iterate over all actions
2526  while (!actions.is_empty()) {
2527  if (actions.earliest_time() > end_time_propagation) {
2528  break;
2529  }
2530  // get next action
2531  ActionPtr act = actions.pop();
2532  if (!act->is_valid(particles)) {
2534  logg[LExperiment].debug(~einhard::DRed(), "✘ ", act,
2535  " (discarded: invalid)");
2536  continue;
2537  }
2538  logg[LExperiment].debug(~einhard::Green(), "✔ ", act,
2539  ", action time = ", act->time_of_execution());
2540 
2541  /* (1) Propagate to the next action. */
2542  propagate_and_shine(act->time_of_execution(), particles);
2543 
2544  /* (2) Perform action.
2545  *
2546  * Update the positions of the incoming particles, because the information
2547  * in the action object will be outdated as the particles have been
2548  * propagated since the construction of the action. */
2549  act->update_incoming(particles);
2550  const bool performed = perform_action(*act, i_ensemble);
2551 
2552  /* No need to update actions for outgoing particles
2553  * if the action is not performed. */
2554  if (!performed) {
2555  continue;
2556  }
2557 
2558  /* (3) Update actions for newly-produced particles. */
2559 
2560  const double end_time_timestep = parameters_.labclock->next_time();
2561  // New actions are always search until the end of the current timestep
2562  const double time_left = end_time_timestep - act->time_of_execution();
2563  const ParticleList &outgoing_particles = act->outgoing_particles();
2564  // Grid cell volume set to zero, since there is no grid
2565  const double gcell_vol = 0.0;
2566  for (const auto &finder : action_finders_) {
2567  // Outgoing particles can still decay, cross walls...
2568  actions.insert(finder->find_actions_in_cell(outgoing_particles, time_left,
2569  gcell_vol, beam_momentum_));
2570  // ... and collide with other particles.
2571  actions.insert(finder->find_actions_with_surrounding_particles(
2572  outgoing_particles, particles, time_left, beam_momentum_));
2573  }
2574 
2576  }
2577 
2578  propagate_and_shine(end_time_propagation, particles);
2579 }
A stream modifier that allows to colorize the log output.
Definition: einhard.hpp:143
void propagate_and_shine(double to_time, Particles &particles)
Propagate all particles until time to_time without any interactions and shine dileptons.
Definition: experiment.h:2492
void check_interactions_total(uint64_t interactions_total)
Make sure interactions_total can be represented as a 32-bit integer.
Definition: experiment.h:2510

◆ intermediate_output()

template<typename Modus >
void smash::Experiment< Modus >::intermediate_output
private

Intermediate output during an event.

Auxiliary variable to communicate the time in the computational frame at the functions printing the thermodynamics lattice output

Definition at line 2582 of file experiment.h.

2582  {
2583  const uint64_t wall_actions_this_interval =
2586  const uint64_t interactions_this_interval = interactions_total_ -
2588  wall_actions_this_interval;
2590  double E_mean_field = 0.0;
2593  double computational_frame_time = 0.0;
2594  if (potentials_) {
2595  // using the lattice is necessary
2596  if ((jmu_B_lat_ != nullptr)) {
2598  EM_lat_.get(), parameters_);
2599  /*
2600  * Mean field calculated in a box should remain approximately constant if
2601  * the system is in equilibrium, and so deviations from its original value
2602  * may signal a phase transition or other dynamical process. This
2603  * comparison only makes sense in the Box Modus, hence the condition.
2604  */
2605  if (modus_.is_box()) {
2606  double tmp = (E_mean_field - initial_mean_field_energy_) /
2607  (E_mean_field + initial_mean_field_energy_);
2608  /*
2609  * This is displayed when the system evolves away from its initial
2610  * configuration (which is when the total mean field energy in the box
2611  * deviates from its initial value).
2612  */
2613  if (std::abs(tmp) > 0.01) {
2614  logg[LExperiment].info()
2615  << "\n\n\n\t The mean field at t = "
2616  << parameters_.outputclock->current_time()
2617  << " [fm] differs from the mean field at t = 0:"
2618  << "\n\t\t initial_mean_field_energy_ = "
2619  << initial_mean_field_energy_ << " [GeV]"
2620  << "\n\t\t abs[(E_MF - E_MF(t=0))/(E_MF + E_MF(t=0))] = "
2621  << std::abs(tmp)
2622  << "\n\t\t E_MF/E_MF(t=0) = "
2623  << E_mean_field / initial_mean_field_energy_ << "\n\n";
2624  }
2625  }
2626  }
2627  }
2628 
2630  ensembles_, interactions_this_interval, conserved_initial_, time_start_,
2631  parameters_.outputclock->current_time(), E_mean_field,
2633  const LatticeUpdate lat_upd = LatticeUpdate::AtOutput;
2634 
2635  // save evolution data
2636  if (!(modus_.is_box() && parameters_.outputclock->current_time() <
2637  modus_.equilibration_time())) {
2638  for (const auto &output : outputs_) {
2639  if (output->is_dilepton_output() || output->is_photon_output() ||
2640  output->is_IC_output()) {
2641  continue;
2642  }
2643  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2644  auto event_info = fill_event_info(
2645  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2647 
2648  output->at_intermediate_time(ensembles_[i_ens], parameters_.outputclock,
2649  density_param_, event_info);
2650  computational_frame_time = event_info.current_time;
2651  }
2652  // For thermodynamic output
2653  output->at_intermediate_time(ensembles_, parameters_.outputclock,
2654  density_param_);
2655 
2656  // Thermodynamic output on the lattice versus time
2657  if (printout_rho_eckart_) {
2658  switch (dens_type_lattice_printout_) {
2659  case DensityType::Baryon:
2661  density_param_, ensembles_, false);
2662  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2664  output->thermodynamics_lattice_output(*jmu_B_lat_,
2665  computational_frame_time);
2666  break;
2668  update_lattice(jmu_I3_lat_.get(), lat_upd,
2670  ensembles_, false);
2671  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2673  *jmu_I3_lat_);
2674  output->thermodynamics_lattice_output(*jmu_I3_lat_,
2675  computational_frame_time);
2676  break;
2677  case DensityType::None:
2678  break;
2679  default:
2680  update_lattice(jmu_custom_lat_.get(), lat_upd,
2682  ensembles_, false);
2683  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2685  *jmu_custom_lat_);
2686  output->thermodynamics_lattice_output(*jmu_custom_lat_,
2687  computational_frame_time);
2688  }
2689  }
2692  density_param_, ensembles_, false);
2693  if (printout_tmn_) {
2694  output->thermodynamics_output(ThermodynamicQuantity::Tmn,
2696  output->thermodynamics_lattice_output(
2697  ThermodynamicQuantity::Tmn, *Tmn_, computational_frame_time);
2698  }
2699  if (printout_tmn_landau_) {
2700  output->thermodynamics_output(ThermodynamicQuantity::TmnLandau,
2702  output->thermodynamics_lattice_output(
2704  computational_frame_time);
2705  }
2706  if (printout_v_landau_) {
2707  output->thermodynamics_output(ThermodynamicQuantity::LandauVelocity,
2709  output->thermodynamics_lattice_output(
2711  computational_frame_time);
2712  }
2713  }
2714  if (EM_lat_) {
2715  output->fields_output("Efield", "Bfield", *EM_lat_);
2716  }
2717  if (printout_j_QBS_) {
2718  output->thermodynamics_lattice_output(
2719  *j_QBS_lat_, computational_frame_time, ensembles_, density_param_);
2720  }
2721 
2722  if (thermalizer_) {
2723  output->thermodynamics_output(*thermalizer_);
2724  }
2725  }
2726  }
2727 }
LatticeUpdate
Enumerator option for lattice updates.
Definition: lattice.h:36

◆ update_potentials()

template<typename Modus >
void smash::Experiment< Modus >::update_potentials
private

Recompute potentials on lattices if necessary.

Definition at line 2730 of file experiment.h.

2730  {
2731  if (potentials_) {
2732  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
2737  parameters_.labclock->timestep_duration(), true);
2738  }
2739  if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
2740  jmu_B_lat_ != nullptr) {
2745  parameters_.labclock->timestep_duration(), true);
2746  const size_t UBlattice_size = UB_lat_->size();
2747  for (size_t i = 0; i < UBlattice_size; i++) {
2748  auto jB = (*jmu_B_lat_)[i];
2749  const FourVector flow_four_velocity_B =
2750  std::abs(jB.rho()) > very_small_double ? jB.jmu_net() / jB.rho()
2751  : FourVector();
2752  double baryon_density = jB.rho();
2753  ThreeVector baryon_grad_j0 = jB.grad_j0();
2754  ThreeVector baryon_dvecj_dt = jB.dvecj_dt();
2755  ThreeVector baryon_curl_vecj = jB.curl_vecj();
2756  if (potentials_->use_skyrme()) {
2757  (*UB_lat_)[i] =
2758  flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
2759  (*FB_lat_)[i] =
2760  potentials_->skyrme_force(baryon_density, baryon_grad_j0,
2761  baryon_dvecj_dt, baryon_curl_vecj);
2762  }
2763  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
2764  auto jI3 = (*jmu_I3_lat_)[i];
2765  const FourVector flow_four_velocity_I3 =
2766  std::abs(jI3.rho()) > very_small_double
2767  ? jI3.jmu_net() / jI3.rho()
2768  : FourVector();
2769  (*UI3_lat_)[i] = flow_four_velocity_I3 *
2770  potentials_->symmetry_pot(jI3.rho(), baryon_density);
2771  (*FI3_lat_)[i] = potentials_->symmetry_force(
2772  jI3.rho(), jI3.grad_j0(), jI3.dvecj_dt(), jI3.curl_vecj(),
2773  baryon_density, baryon_grad_j0, baryon_dvecj_dt,
2774  baryon_curl_vecj);
2775  }
2776  }
2777  }
2778  if (potentials_->use_coulomb()) {
2781  for (size_t i = 0; i < EM_lat_->size(); i++) {
2782  ThreeVector electric_field = {0., 0., 0.};
2783  ThreeVector position = jmu_el_lat_->cell_center(i);
2784  jmu_el_lat_->integrate_volume(electric_field,
2786  potentials_->coulomb_r_cut(), position);
2787  ThreeVector magnetic_field = {0., 0., 0.};
2788  jmu_el_lat_->integrate_volume(magnetic_field,
2790  potentials_->coulomb_r_cut(), position);
2791  (*EM_lat_)[i] = std::make_pair(electric_field, magnetic_field);
2792  }
2793  } // if ((potentials_->use_skyrme() || ...
2794  if (potentials_->use_vdf() && jmu_B_lat_ != nullptr) {
2799  parameters_.labclock->timestep_duration(), true);
2802  fields_lat_.get(), old_fields_auxiliary_.get(),
2805  parameters_.labclock->timestep_duration());
2806  }
2807  const size_t UBlattice_size = UB_lat_->size();
2808  for (size_t i = 0; i < UBlattice_size; i++) {
2809  auto jB = (*jmu_B_lat_)[i];
2810  (*UB_lat_)[i] = potentials_->vdf_pot(jB.rho(), jB.jmu_net());
2813  (*FB_lat_)[i] = potentials_->vdf_force(
2814  jB.rho(), jB.drho_dxnu().x0(), jB.drho_dxnu().threevec(),
2815  jB.grad_rho_cross_vecj(), jB.jmu_net().x0(), jB.grad_j0(),
2816  jB.jmu_net().threevec(), jB.dvecj_dt(), jB.curl_vecj());
2817  break;
2819  auto Amu = (*fields_lat_)[i];
2820  (*FB_lat_)[i] = potentials_->vdf_force(
2821  Amu.grad_A0(), Amu.dvecA_dt(), Amu.curl_vecA());
2822  break;
2823  }
2824  } // for (size_t i = 0; i < UBlattice_size; i++)
2825  } // if potentials_->use_vdf()
2826  }
2827 }
std::unique_ptr< RectangularLattice< FourVector > > new_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0 + dt.
Definition: experiment.h:499
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > fields_four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of Amu.
Definition: experiment.h:502
std::unique_ptr< FieldsLattice > fields_lat_
Mean-field A^mu on the lattice.
Definition: experiment.h:444
std::unique_ptr< RectangularLattice< FourVector > > UB_lat_
Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-vel...
Definition: experiment.h:462
std::unique_ptr< RectangularLattice< FourVector > > old_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0.
Definition: experiment.h:497
std::unique_ptr< DensityLattice > jmu_el_lat_
Electric charge density on the lattice.
Definition: experiment.h:441
static ThreeVector B_field_integrand(ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
Integrand for calculating the magnetic field using the Biot-Savart formula.
Definition: potentials.h:338
static ThreeVector E_field_integrand(ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
Integrand for calculating the electric field.
Definition: potentials.h:320
constexpr double very_small_double
A very small double, used to avoid division by zero.
Definition: constants.h:40
void update_fields_lattice(RectangularLattice< FieldsOnLattice > *fields_lat, RectangularLattice< FourVector > *old_fields, RectangularLattice< FourVector > *new_fields, RectangularLattice< std::array< FourVector, 4 >> *fields_four_grad_lattice, DensityLattice *jmu_B_lat, const LatticeUpdate fields_lat_update, const Potentials &potentials, const double time_step)
Updates the contents on the lattice of FieldsOnLattice type.
Definition: fields.cc:14
FieldDerivativesMode field_derivatives_mode
mode of calculating field derivatives

◆ compute_min_cell_length()

template<typename Modus >
double smash::Experiment< Modus >::compute_min_cell_length ( double  dt) const
inlineprivate

Calculate the minimal size for the grid cells such that the ScatterActionsFinder will find all collisions within the maximal transverse distance (which is determined by the maximal cross section).

Parameters
[in]dtThe current time step size [fm]
Returns
The minimal required size of cells

Definition at line 338 of file experiment.h.

338  {
341  }
342  return std::sqrt(4 * dt * dt + max_transverse_distance_sqr_);
343  }
double fixed_min_cell_length
Fixed minimal grid cell length (in fm).

◆ next_output_time()

template<typename Modus >
double smash::Experiment< Modus >::next_output_time ( ) const
inlineprivate

Shortcut for next output time.

Definition at line 346 of file experiment.h.

346  {
347  return parameters_.outputclock->next_time();
348  }

◆ count_nonempty_ensembles()

template<typename Modus >
void smash::Experiment< Modus >::count_nonempty_ensembles
private

Counts the number of ensembles in wich interactions took place at the end of an event.

Definition at line 3033 of file experiment.h.

3033  {
3034  for (bool has_interaction : projectile_target_interact_) {
3035  if (has_interaction) {
3037  }
3038  }
3039 }
int nonempty_ensembles_
Number of ensembles containing an interaction.
Definition: experiment.h:573

◆ is_finished()

template<typename Modus >
bool smash::Experiment< Modus >::is_finished
private

Checks wether the desired number events have been calculated.

Returns
wether the experiment is is_finished

Definition at line 3042 of file experiment.h.

3042  {
3044  return event_ >= nevents_;
3045  }
3048  return true;
3049  }
3050  if (event_ >= max_events_) {
3051  logg[LExperiment].warn()
3052  << "Maximum number of events (" << max_events_
3053  << ") exceeded. Stopping calculation. "
3054  << "The fraction of empty ensembles is "
3055  << (1.0 - static_cast<double>(nonempty_ensembles_) /
3057  << ". If this fraction is expected, try increasing the "
3058  "Maximum_Ensembles_Run.";
3059  return true;
3060  }
3061  return false;
3062  }
3063  throw std::runtime_error("Event counting option is invalid");
3064  return false;
3065 }
int minimum_nonempty_ensembles_
The number of ensembles, in which interactions take place, to be calculated.
Definition: experiment.h:559
const int nevents_
Number of events.
Definition: experiment.h:549
int max_events_
Maximum number of events to be calculated in order obtain the desired number of non-empty events usin...
Definition: experiment.h:579
EventCounting event_counting_
The way in which the number of calculated events is specified.
Definition: experiment.h:567
@ FixedNumber
The desired number of events is simulated disregarding of whether an interaction took place.
@ MinimumNonEmpty
Events are simulated until there are at least a given number of ensembles in which an interaction too...

Friends And Related Function Documentation

◆ ExperimentBase

template<typename Modus >
friend class ExperimentBase
friend

Definition at line 189 of file experiment.h.

Member Data Documentation

◆ parameters_

template<typename Modus >
ExperimentParameters smash::Experiment< Modus >::parameters_
private

Struct of several member variables.

These variables are combined into a struct for efficient input to functions outside of this class.

Definition at line 368 of file experiment.h.

◆ density_param_

template<typename Modus >
DensityParameters smash::Experiment< Modus >::density_param_
private

Structure to precalculate and hold parameters for density computations.

Definition at line 371 of file experiment.h.

◆ modus_

template<typename Modus >
Modus smash::Experiment< Modus >::modus_
private

Instance of the Modus template parameter.

May store modus-specific data and contains modus-specific function implementations.

Definition at line 377 of file experiment.h.

◆ ensembles_

template<typename Modus >
std::vector<Particles> smash::Experiment< Modus >::ensembles_
private

Complete particle list, all ensembles in one vector.

Definition at line 380 of file experiment.h.

◆ potentials_

template<typename Modus >
std::unique_ptr<Potentials> smash::Experiment< Modus >::potentials_
private

An instance of potentials class, that stores parameters of potentials, calculates them and their gradients.

Definition at line 386 of file experiment.h.

◆ pauli_blocker_

template<typename Modus >
std::unique_ptr<PauliBlocker> smash::Experiment< Modus >::pauli_blocker_
private

An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and computes phase-space density.

Definition at line 392 of file experiment.h.

◆ outputs_

template<typename Modus >
OutputsList smash::Experiment< Modus >::outputs_
private

A list of output formaters.

They will be called to write the state of the particles to file.

Definition at line 398 of file experiment.h.

◆ dilepton_output_

template<typename Modus >
OutputPtr smash::Experiment< Modus >::dilepton_output_
private

The Dilepton output.

Definition at line 401 of file experiment.h.

◆ photon_output_

template<typename Modus >
OutputPtr smash::Experiment< Modus >::photon_output_
private

The Photon output.

Definition at line 404 of file experiment.h.

◆ projectile_target_interact_

template<typename Modus >
std::vector<bool> smash::Experiment< Modus >::projectile_target_interact_
private

Whether the projectile and the target collided.

One value for each ensemble.

Definition at line 410 of file experiment.h.

◆ beam_momentum_

template<typename Modus >
std::vector<FourVector> smash::Experiment< Modus >::beam_momentum_ = {}
private

The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen.

It's only valid in the ColliderModus, so is set as an empty vector by default.

Definition at line 417 of file experiment.h.

◆ action_finders_

template<typename Modus >
std::vector<std::unique_ptr<ActionFinderInterface> > smash::Experiment< Modus >::action_finders_
private

The Action finder objects.

Definition at line 420 of file experiment.h.

◆ dilepton_finder_

template<typename Modus >
std::unique_ptr<DecayActionsFinderDilepton> smash::Experiment< Modus >::dilepton_finder_
private

The Dilepton Action Finder.

Definition at line 423 of file experiment.h.

◆ photon_finder_

template<typename Modus >
std::unique_ptr<ActionFinderInterface> smash::Experiment< Modus >::photon_finder_
private

The (Scatter) Actions Finder for Direct Photons.

Definition at line 426 of file experiment.h.

◆ n_fractional_photons_

template<typename Modus >
int smash::Experiment< Modus >::n_fractional_photons_
private

Number of fractional photons produced per single reaction.

Definition at line 429 of file experiment.h.

◆ j_QBS_lat_

template<typename Modus >
std::unique_ptr<DensityLattice> smash::Experiment< Modus >::j_QBS_lat_
private

4-current for j_QBS lattice output

Definition at line 432 of file experiment.h.

◆ jmu_B_lat_

template<typename Modus >
std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_B_lat_
private

Baryon density on the lattice.

Definition at line 435 of file experiment.h.

◆ jmu_I3_lat_

template<typename Modus >
std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_I3_lat_
private

Isospin projection density on the lattice.

Definition at line 438 of file experiment.h.

◆ jmu_el_lat_

template<typename Modus >
std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_el_lat_
private

Electric charge density on the lattice.

Definition at line 441 of file experiment.h.

◆ fields_lat_

template<typename Modus >
std::unique_ptr<FieldsLattice> smash::Experiment< Modus >::fields_lat_
private

Mean-field A^mu on the lattice.

Definition at line 444 of file experiment.h.

◆ jmu_custom_lat_

template<typename Modus >
std::unique_ptr<DensityLattice> smash::Experiment< Modus >::jmu_custom_lat_
private

Custom density on the lattices.

In the config user asks for some kind of density for printout. Baryon and isospin projection density are anyway needed for potentials. If user asks for some other density type for printout, it will be handled using jmu_custom variable.

Definition at line 453 of file experiment.h.

◆ dens_type_lattice_printout_

template<typename Modus >
DensityType smash::Experiment< Modus >::dens_type_lattice_printout_ = DensityType::None
private

Type of density for lattice printout.

Definition at line 456 of file experiment.h.

◆ UB_lat_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::UB_lat_ = nullptr
private

Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-velocity.

Definition at line 462 of file experiment.h.

◆ UI3_lat_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::UI3_lat_ = nullptr
private

Lattices for symmetry potentials (evaluated in the local rest frame) times the isospin flow 4-velocity.

Definition at line 468 of file experiment.h.

◆ FB_lat_

template<typename Modus >
std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::FB_lat_
private

Lattices for the electric and magnetic components of the Skyrme or VDF force.

Definition at line 475 of file experiment.h.

◆ FI3_lat_

template<typename Modus >
std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::FI3_lat_
private

Lattices for the electric and magnetic component of the symmetry force.

Definition at line 479 of file experiment.h.

◆ EM_lat_

template<typename Modus >
std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector> > > smash::Experiment< Modus >::EM_lat_
private

Lattices for electric and magnetic field in fm^-2.

Definition at line 483 of file experiment.h.

◆ Tmn_

template<typename Modus >
std::unique_ptr<RectangularLattice<EnergyMomentumTensor> > smash::Experiment< Modus >::Tmn_
private

Lattices of energy-momentum tensors for printout.

Definition at line 486 of file experiment.h.

◆ old_jmu_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::old_jmu_auxiliary_
private

Auxiliary lattice for values of jmu at a time step t0.

Definition at line 489 of file experiment.h.

◆ new_jmu_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::new_jmu_auxiliary_
private

Auxiliary lattice for values of jmu at a time step t0 + dt.

Definition at line 491 of file experiment.h.

◆ four_gradient_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<std::array<FourVector, 4> > > smash::Experiment< Modus >::four_gradient_auxiliary_
private

Auxiliary lattice for calculating the four-gradient of jmu.

Definition at line 494 of file experiment.h.

◆ old_fields_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::old_fields_auxiliary_
private

Auxiliary lattice for values of Amu at a time step t0.

Definition at line 497 of file experiment.h.

◆ new_fields_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<FourVector> > smash::Experiment< Modus >::new_fields_auxiliary_
private

Auxiliary lattice for values of Amu at a time step t0 + dt.

Definition at line 499 of file experiment.h.

◆ fields_four_gradient_auxiliary_

template<typename Modus >
std::unique_ptr<RectangularLattice<std::array<FourVector, 4> > > smash::Experiment< Modus >::fields_four_gradient_auxiliary_
private

Auxiliary lattice for calculating the four-gradient of Amu.

Definition at line 502 of file experiment.h.

◆ printout_rho_eckart_

template<typename Modus >
bool smash::Experiment< Modus >::printout_rho_eckart_ = false
private

Whether to print the Eckart rest frame density.

Definition at line 505 of file experiment.h.

◆ printout_tmn_

template<typename Modus >
bool smash::Experiment< Modus >::printout_tmn_ = false
private

Whether to print the energy-momentum tensor.

Definition at line 508 of file experiment.h.

◆ printout_tmn_landau_

template<typename Modus >
bool smash::Experiment< Modus >::printout_tmn_landau_ = false
private

Whether to print the energy-momentum tensor in Landau frame.

Definition at line 511 of file experiment.h.

◆ printout_v_landau_

template<typename Modus >
bool smash::Experiment< Modus >::printout_v_landau_ = false
private

Whether to print the 4-velocity in Landau frame.

Definition at line 514 of file experiment.h.

◆ printout_j_QBS_

template<typename Modus >
bool smash::Experiment< Modus >::printout_j_QBS_ = false
private

Whether to print the Q, B, S 4-currents.

Definition at line 517 of file experiment.h.

◆ printout_lattice_td_

template<typename Modus >
bool smash::Experiment< Modus >::printout_lattice_td_ = false
private

Whether to print the thermodynamics quantities evaluated on the lattices.

Definition at line 520 of file experiment.h.

◆ printout_full_lattice_any_td_

template<typename Modus >
bool smash::Experiment< Modus >::printout_full_lattice_any_td_ = false
private

Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in any format.

Definition at line 524 of file experiment.h.

◆ thermalizer_

template<typename Modus >
std::unique_ptr<GrandCanThermalizer> smash::Experiment< Modus >::thermalizer_
private

Instance of class used for forced thermalization.

Definition at line 527 of file experiment.h.

◆ process_string_ptr_

template<typename Modus >
StringProcess* smash::Experiment< Modus >::process_string_ptr_
private

Pointer to the string process class object, which is used to set the random seed for PYTHIA objects in each event.

Definition at line 533 of file experiment.h.

◆ nevents_

template<typename Modus >
const int smash::Experiment< Modus >::nevents_ = 0
private

Number of events.

Event is a single simulation of a physical phenomenon: elementary particle or nucleus-nucleus collision. Result of a single SMASH event is random (by construction) as well as result of one collision in nature. To compare simulation with experiment one has to take ensemble averages, i.e. perform simulation and real experiment many times and compare average results.

nevents_ is number of times single phenomenon (particle or nucleus-nucleus collision) will be simulated.

Definition at line 549 of file experiment.h.

◆ minimum_nonempty_ensembles_

template<typename Modus >
int smash::Experiment< Modus >::minimum_nonempty_ensembles_ = 0
private

The number of ensembles, in which interactions take place, to be calculated.

Can be specified as an inout instead of the number of events. In this case events will be calculated until this number of ensembles is reached.

Definition at line 559 of file experiment.h.

◆ event_counting_

template<typename Modus >
EventCounting smash::Experiment< Modus >::event_counting_ = EventCounting::Invalid
private

The way in which the number of calculated events is specified.

Can be either a fixed number of simulated events or a minimum number of events that contain interactions.

Definition at line 567 of file experiment.h.

◆ event_

template<typename Modus >
int smash::Experiment< Modus >::event_ = 0
private

Current event.

Definition at line 570 of file experiment.h.

◆ nonempty_ensembles_

template<typename Modus >
int smash::Experiment< Modus >::nonempty_ensembles_ = 0
private

Number of ensembles containing an interaction.

Definition at line 573 of file experiment.h.

◆ max_events_

template<typename Modus >
int smash::Experiment< Modus >::max_events_ = 0
private

Maximum number of events to be calculated in order obtain the desired number of non-empty events using the MinimumNonemptyEnsembles option.

Definition at line 579 of file experiment.h.

◆ end_time_

template<typename Modus >
const double smash::Experiment< Modus >::end_time_
private

simulation time at which the evolution is stopped.

Definition at line 582 of file experiment.h.

◆ delta_time_startup_

template<typename Modus >
const double smash::Experiment< Modus >::delta_time_startup_
private

The clock's timestep size at start up.

Stored here so that the next event will remember this.

Definition at line 589 of file experiment.h.

◆ force_decays_

template<typename Modus >
const bool smash::Experiment< Modus >::force_decays_
private

This indicates whether we force all resonances to decay in the last timestep.

Definition at line 595 of file experiment.h.

◆ use_grid_

template<typename Modus >
const bool smash::Experiment< Modus >::use_grid_
private

This indicates whether to use the grid.

Definition at line 598 of file experiment.h.

◆ metric_

template<typename Modus >
const ExpansionProperties smash::Experiment< Modus >::metric_
private

This struct contains information on the metric to be used.

Definition at line 601 of file experiment.h.

◆ dileptons_switch_

template<typename Modus >
const bool smash::Experiment< Modus >::dileptons_switch_
private

This indicates whether dileptons are switched on.

Definition at line 604 of file experiment.h.

◆ photons_switch_

template<typename Modus >
const bool smash::Experiment< Modus >::photons_switch_
private

This indicates whether photons are switched on.

Definition at line 607 of file experiment.h.

◆ bremsstrahlung_switch_

template<typename Modus >
const bool smash::Experiment< Modus >::bremsstrahlung_switch_
private

This indicates whether bremsstrahlung is switched on.

Definition at line 610 of file experiment.h.

◆ IC_output_switch_

template<typename Modus >
const bool smash::Experiment< Modus >::IC_output_switch_
private

This indicates whether the IC output is enabled.

Definition at line 613 of file experiment.h.

◆ time_step_mode_

template<typename Modus >
const TimeStepMode smash::Experiment< Modus >::time_step_mode_
private

This indicates whether to use time steps.

Definition at line 616 of file experiment.h.

◆ max_transverse_distance_sqr_

template<typename Modus >
double smash::Experiment< Modus >::max_transverse_distance_sqr_ = std::numeric_limits<double>::max()
private

Maximal distance at which particles can interact in case of the geometric criterion, squared.

Definition at line 622 of file experiment.h.

◆ conserved_initial_

template<typename Modus >
QuantumNumbers smash::Experiment< Modus >::conserved_initial_
private

The conserved quantities of the system.

This struct carries the sums of the single particle's various quantities as measured at the beginning of the evolution and can be used to regularly check if they are still good.

Definition at line 631 of file experiment.h.

◆ initial_mean_field_energy_

template<typename Modus >
double smash::Experiment< Modus >::initial_mean_field_energy_
private

The initial total mean field energy in the system.

Note: will only be calculated if lattice is on.

Definition at line 637 of file experiment.h.

◆ time_start_

template<typename Modus >
SystemTimePoint smash::Experiment< Modus >::time_start_ = SystemClock::now()
private

system starting time of the simulation

Definition at line 640 of file experiment.h.

◆ dens_type_

template<typename Modus >
DensityType smash::Experiment< Modus >::dens_type_ = DensityType::None
private

Type of density to be written to collision headers.

Definition at line 643 of file experiment.h.

◆ interactions_total_

template<typename Modus >
uint64_t smash::Experiment< Modus >::interactions_total_ = 0
private

Total number of interactions for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 649 of file experiment.h.

◆ previous_interactions_total_

template<typename Modus >
uint64_t smash::Experiment< Modus >::previous_interactions_total_ = 0
private

Total number of interactions for previous timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 655 of file experiment.h.

◆ wall_actions_total_

template<typename Modus >
uint64_t smash::Experiment< Modus >::wall_actions_total_ = 0
private

Total number of wall-crossings for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 661 of file experiment.h.

◆ previous_wall_actions_total_

template<typename Modus >
uint64_t smash::Experiment< Modus >::previous_wall_actions_total_ = 0
private

Total number of wall-crossings for previous timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 667 of file experiment.h.

◆ total_pauli_blocked_

template<typename Modus >
uint64_t smash::Experiment< Modus >::total_pauli_blocked_ = 0
private

Total number of Pauli-blockings for current timestep.

For timestepless mode the whole run time is considered as one timestep.

Definition at line 673 of file experiment.h.

◆ total_hypersurface_crossing_actions_

template<typename Modus >
uint64_t smash::Experiment< Modus >::total_hypersurface_crossing_actions_ = 0
private

Total number of particles removed from the evolution in hypersurface crossing actions.

Definition at line 679 of file experiment.h.

◆ discarded_interactions_total_

template<typename Modus >
uint64_t smash::Experiment< Modus >::discarded_interactions_total_ = 0
private

Total number of discarded interactions, because they were invalidated before they could be performed.

Definition at line 685 of file experiment.h.

◆ total_energy_removed_

template<typename Modus >
double smash::Experiment< Modus >::total_energy_removed_ = 0.0
private

Total energy removed from the system in hypersurface crossing actions.

Definition at line 690 of file experiment.h.

◆ total_energy_violated_by_Pythia_

template<typename Modus >
double smash::Experiment< Modus >::total_energy_violated_by_Pythia_ = 0.0
private

Total energy violation introduced by Pythia.

Definition at line 695 of file experiment.h.

◆ kinematic_cuts_for_IC_output_

template<typename Modus >
bool smash::Experiment< Modus >::kinematic_cuts_for_IC_output_ = false
private

This indicates whether kinematic cuts are enabled for the IC output.

Definition at line 698 of file experiment.h.

◆ seed_

template<typename Modus >
int64_t smash::Experiment< Modus >::seed_ = -1
private

random seed for the next event.

Definition at line 701 of file experiment.h.


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