Version: SMASH-3.4
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 193 of file experiment.h.

Inheritance diagram for smash::Experiment< Modus >:
[legend]
Collaboration diagram for smash::Experiment< Modus >:
[legend]

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_interactions ()
 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...
 
bool printout_coulomb_vtk_ = false
 Whether to write the electric and magnetic fields to VTK files. 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...
 
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 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 dileptons_bremsstrahlung_switch_
 This indicates whether dilepton production via bremsstrahlung is switched on. More...
 
const bool photons_switch_
 This indicates whether photons are switched on. More...
 
const bool photons_bremsstrahlung_switch_
 This indicates whether bremsstrahlung is switched on. More...
 
const bool IC_switch_
 This indicates whether the experiment will be used as initial condition for hydrodynamics. More...
 
const bool IC_dynamic_
 This indicates if the IC is dynamic. 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 3313 of file experiment.h.

3313  {
3314  const auto &mainlog = logg[LMain];
3315  for (event_ = 0; !is_finished(); event_++) {
3316  mainlog.info() << "Event " << event_;
3317 
3318  // Sample initial particles, start clock, some printout and book-keeping
3320 
3322 
3324 
3325  // Output at event end
3326  final_output();
3327  }
3328 }
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:2514
void initialize_new_event()
This is called in the beginning of each event.
Definition: experiment.h:2096
bool is_finished()
Checks wether the desired number events have been calculated.
Definition: experiment.h:3282
int event_
Current event.
Definition: experiment.h:578
void do_final_interactions()
Performs the final decays of an event.
Definition: experiment.h:3081
void final_output()
Output at the end of an event.
Definition: experiment.h:3129
const double end_time_
simulation time at which the evolution is stopped.
Definition: experiment.h:590
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > & logg
An array that stores all pre-configured Logger objects.
Definition: logging.h:245
static constexpr int LMain
Definition: experiment.h:94

◆ 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 2096 of file experiment.h.

2096  {
2098  logg[LExperiment].info() << "random number seed: " << seed_;
2099  /* Set seed for the next event. It has to be positive, so it can be entered
2100  * in the config.
2101  *
2102  * We have to be careful about the minimal integer, whose absolute value
2103  * cannot be represented. */
2104  int64_t r = random::advance();
2105  while (r == INT64_MIN) {
2106  r = random::advance();
2107  }
2108  seed_ = std::abs(r);
2109  /* Set the random seed used in PYTHIA hadronization
2110  * to be same with the SMASH one.
2111  * In this way we ensure that the results are reproducible
2112  * for every event if one knows SMASH random seed. */
2113  if (process_string_ptr_ != NULL) {
2115  }
2116 
2117  for (Particles &particles : ensembles_) {
2118  particles.reset();
2119  }
2120 
2121  // Sample particles according to the initial conditions
2122  double start_time = -1.0;
2123 
2124  // Sample impact parameter only once per all ensembles
2125  // It should be the same for all ensembles
2126  if (modus_.is_collider()) {
2127  modus_.sample_impact();
2128  logg[LExperiment].info("Impact parameter = ", modus_.impact_parameter(),
2129  " fm");
2130  }
2131  for (Particles &particles : ensembles_) {
2132  start_time = modus_.initial_conditions(&particles, parameters_);
2133  }
2134  /* For box modus make sure that particles are in the box. In principle, after
2135  * a correct initialization they should be, so this is just playing it safe.
2136  */
2137  for (Particles &particles : ensembles_) {
2138  modus_.impose_boundary_conditions(&particles, outputs_);
2139  }
2140  // Reset the simulation clock
2141  double timestep = delta_time_startup_;
2142 
2143  switch (time_step_mode_) {
2144  case TimeStepMode::Fixed:
2145  break;
2146  case TimeStepMode::None:
2147  timestep = end_time_ - start_time;
2148  // Take care of the box modus + timestepless propagation
2149  const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
2150  if (max_dt > 0. && max_dt < timestep) {
2151  timestep = max_dt;
2152  }
2153  break;
2154  }
2155  std::unique_ptr<UniformClock> clock_for_this_event;
2156  if (modus_.is_list() && (timestep < 0.0)) {
2157  throw std::runtime_error(
2158  "Timestep for the given event is negative. \n"
2159  "This might happen if the formation times of the input particles are "
2160  "larger than the specified end time of the simulation.");
2161  }
2162  clock_for_this_event =
2163  std::make_unique<UniformClock>(start_time, timestep, end_time_);
2164  parameters_.labclock = std::move(clock_for_this_event);
2165 
2166  // Reset the output clock
2167  parameters_.outputclock->reset(start_time, true);
2168  // remove time before starting time in case of custom output times.
2169  parameters_.outputclock->remove_times_in_past(start_time);
2170 
2171  logg[LExperiment].debug(
2172  "Lab clock: t_start = ", parameters_.labclock->current_time(),
2173  ", dt = ", parameters_.labclock->timestep_duration());
2174 
2175  /* Save the initial conserved quantum numbers and total momentum in
2176  * the system for conservation checks */
2177  conserved_initial_ = QuantumNumbers(ensembles_);
2178  wall_actions_total_ = 0;
2180  interactions_total_ = 0;
2186  total_energy_removed_ = 0.0;
2188  // Print output headers
2189  logg[LExperiment].info() << hline;
2190  logg[LExperiment].info() << "Time[fm] Ekin[GeV] E_MF[GeV] ETotal[GeV] "
2191  << "ETot/N[GeV] D(ETot/N)[GeV] Scatt&Decays "
2192  << "Particles Comp.Time";
2193  logg[LExperiment].info() << hline;
2194  double E_mean_field = 0.0;
2195  if (potentials_) {
2196  // update_potentials();
2197  // if (parameters.outputclock->current_time() == 0.0 )
2198  // using the lattice is necessary
2199  if ((jmu_B_lat_ != nullptr)) {
2204  parameters_.labclock->timestep_duration(), true);
2205  // Because there was no lattice at t=-Delta_t, the time derivatives
2206  // drho_dt and dj^mu/dt at t=0 are huge, while they shouldn't be; we
2207  // overwrite the time derivative to zero by hand.
2208  for (auto &node : *jmu_B_lat_) {
2209  node.overwrite_drho_dt_to_zero();
2210  node.overwrite_djmu_dt_to_zero();
2211  }
2213  EM_lat_.get(), parameters_);
2214  }
2215  }
2216  initial_mean_field_energy_ = E_mean_field;
2219  parameters_.labclock->current_time(), E_mean_field,
2221 
2222  // Output at event start
2223  for (const auto &output : outputs_) {
2224  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2225  auto event_info = fill_event_info(
2226  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2228  output->at_eventstart(ensembles_[i_ens], {event_, i_ens}, event_info);
2229  }
2230  // For thermodynamic output
2231  output->at_eventstart(ensembles_, event_);
2232  // For thermodynamic lattice output
2234  if (printout_rho_eckart_) {
2235  switch (dens_type_lattice_printout_) {
2236  case DensityType::Baryon:
2239  break;
2243  break;
2244  case DensityType::None:
2245  break;
2246  default:
2249  *jmu_custom_lat_);
2250  }
2251  }
2252  if (printout_tmn_) {
2253  output->at_eventstart(event_, ThermodynamicQuantity::Tmn,
2255  }
2256  if (printout_tmn_landau_) {
2259  }
2260  if (printout_v_landau_) {
2263  }
2264  if (printout_j_QBS_) {
2265  output->at_eventstart(event_, ThermodynamicQuantity::j_QBS,
2267  }
2268  }
2269  }
2270 
2271  /* In the ColliderModus, if Fermi motion is frozen, assign the beam momenta
2272  * to the nucleons in both the projectile and the target. Every ensemble
2273  * gets the same beam momenta, so no need to create beam_momenta_ vector
2274  * for every ensemble.
2275  */
2276  if (modus_.is_collider() && modus_.fermi_motion() == FermiMotion::Frozen) {
2277  for (ParticleData &particle : ensembles_[0]) {
2278  const double m = particle.effective_mass();
2279  double v_beam = 0.0;
2280  if (particle.belongs_to() == BelongsTo::Projectile) {
2281  v_beam = modus_.velocity_projectile();
2282  } else if (particle.belongs_to() == BelongsTo::Target) {
2283  v_beam = modus_.velocity_target();
2284  }
2285  const double gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
2286  beam_momentum_.emplace_back(
2287  FourVector(gamma * m, 0.0, 0.0, gamma * v_beam * m));
2288  } // loop over particles
2289  }
2290 }
double initial_mean_field_energy_
The initial total mean field energy in the system.
Definition: experiment.h:651
bool printout_tmn_
Whether to print the energy-momentum tensor.
Definition: experiment.h:513
QuantumNumbers conserved_initial_
The conserved quantities of the system.
Definition: experiment.h:645
DensityParameters density_param_
Structure to precalculate and hold parameters for density computations.
Definition: experiment.h:376
double total_energy_removed_
Total energy removed from the system in hypersurface crossing actions.
Definition: experiment.h:704
DensityType dens_type_lattice_printout_
Type of density for lattice printout.
Definition: experiment.h:461
double max_transverse_distance_sqr_
Maximal distance at which particles can interact in case of the geometric criterion,...
Definition: experiment.h:636
bool printout_j_QBS_
Whether to print the Q, B, S 4-currents.
Definition: experiment.h:522
const TimeStepMode time_step_mode_
This indicates whether to use time steps.
Definition: experiment.h:630
std::unique_ptr< DensityLattice > jmu_custom_lat_
Custom density on the lattices.
Definition: experiment.h:458
bool printout_full_lattice_any_td_
Whether to print the thermodynamics quantities evaluated on the lattices, point by point,...
Definition: experiment.h:529
std::unique_ptr< DensityLattice > jmu_B_lat_
Baryon density on the lattice.
Definition: experiment.h:440
std::vector< FourVector > beam_momentum_
The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen.
Definition: experiment.h:422
std::unique_ptr< RectangularLattice< FourVector > > new_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0 + dt.
Definition: experiment.h:496
const double delta_time_startup_
The clock's timestep size at start up.
Definition: experiment.h:597
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > > Tmn_
Lattices of energy-momentum tensors for printout.
Definition: experiment.h:491
std::vector< Particles > ensembles_
Complete particle list, all ensembles in one vector.
Definition: experiment.h:385
SystemTimePoint time_start_
system starting time of the simulation
Definition: experiment.h:654
uint64_t previous_wall_actions_total_
Total number of wall-crossings for previous timestep.
Definition: experiment.h:681
OutputsList outputs_
A list of output formaters.
Definition: experiment.h:403
bool printout_rho_eckart_
Whether to print the Eckart rest frame density.
Definition: experiment.h:510
bool printout_v_landau_
Whether to print the 4-velocity in Landau frame.
Definition: experiment.h:519
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > EM_lat_
Lattices for electric and magnetic field in fm^-2.
Definition: experiment.h:488
Modus modus_
Instance of the Modus template parameter.
Definition: experiment.h:382
std::unique_ptr< RectangularLattice< FourVector > > old_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0.
Definition: experiment.h:494
std::unique_ptr< DensityLattice > j_QBS_lat_
4-current for j_QBS lattice output
Definition: experiment.h:437
bool printout_tmn_landau_
Whether to print the energy-momentum tensor in Landau frame.
Definition: experiment.h:516
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of jmu.
Definition: experiment.h:499
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:541
ExperimentParameters parameters_
Struct of several member variables.
Definition: experiment.h:373
std::unique_ptr< Potentials > potentials_
An instance of potentials class, that stores parameters of potentials, calculates them and their grad...
Definition: experiment.h:391
uint64_t wall_actions_total_
Total number of wall-crossings for current timestep.
Definition: experiment.h:675
uint64_t total_hypersurface_crossing_actions_
Total number of particles removed from the evolution in hypersurface crossing actions.
Definition: experiment.h:693
uint64_t interactions_total_
Total number of interactions for current timestep.
Definition: experiment.h:663
std::unique_ptr< DensityLattice > jmu_I3_lat_
Isospin projection density on the lattice.
Definition: experiment.h:443
uint64_t total_pauli_blocked_
Total number of Pauli-blockings for current timestep.
Definition: experiment.h:687
std::vector< bool > projectile_target_interact_
Whether the projectile and the target collided.
Definition: experiment.h:415
int64_t seed_
random seed for the next event.
Definition: experiment.h:715
uint64_t previous_interactions_total_
Total number of interactions for previous timestep.
Definition: experiment.h:669
uint64_t discarded_interactions_total_
Total number of discarded interactions, because they were invalidated before they could be performed.
Definition: experiment.h:699
bool kinematic_cuts_for_IC_output_
This indicates whether kinematic cuts are enabled for the IC output.
Definition: experiment.h:712
double total_energy_violated_by_Pythia_
Total energy violation introduced by Pythia.
Definition: experiment.h:709
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.
constexpr Section output
Section for the output information.
Definition: input_keys.h:205
Engine::result_type advance()
Advance the engine's state and return the generated value.
Definition: random.h:81
void set_seed(T &&seed)
Sets the seed of the random number engine.
Definition: random.h:74
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:586
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:318
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:366
static constexpr int LExperiment
void update_lattice(RectangularLattice< DensityOnLattice > *lat, RectangularLattice< FourVector > *old_jmu, RectangularLattice< FourVector > *new_jmu, RectangularLattice< std::array< FourVector, 4 >> *four_grad_lattice, const LatticeUpdate update, const DensityType dens_type, const DensityParameters &par, const std::vector< Particles > &ensembles, const double time_step, const bool compute_gradient)
Updates the contents on the lattice of DensityOnLattice type.
Definition: density.cc:191
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)
Here is the call graph for this function:

◆ 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 2514 of file experiment.h.

2516  {
2517  if (!add_plist.empty() || !remove_plist.empty()) {
2518  if (ensembles_.size() > 1) {
2519  throw std::runtime_error(
2520  "Adding or removing particles from SMASH is only possible when one "
2521  "ensemble is used.");
2522  }
2523  const double action_time = parameters_.labclock->current_time();
2524  /* Use two if statements. The first one is to check if the particles are
2525  * valid. Since this might remove all particles, a second if statement is
2526  * needed to avoid executing the action in that case.*/
2527  if (!add_plist.empty()) {
2529  }
2530  if (!add_plist.empty()) {
2531  // Create and perform action to add particle(s)
2532  auto action_add_particles = std::make_unique<FreeforallAction>(
2533  ParticleList{}, add_plist, action_time);
2534  perform_action(*action_add_particles, 0);
2535  }
2536  // Also here 2 if statements are needed as above.
2537  if (!remove_plist.empty()) {
2538  validate_and_adjust_particle_list(remove_plist);
2539  }
2540  if (!remove_plist.empty()) {
2541  ParticleList found_particles_to_remove;
2542  for (const auto &particle_to_remove : remove_plist) {
2543  const auto iterator_to_particle_to_be_removed_in_ensemble =
2544  std::find_if(
2545  ensembles_[0].begin(), ensembles_[0].end(),
2546  [&particle_to_remove, &action_time](const ParticleData &p) {
2548  particle_to_remove, p, action_time);
2549  });
2550  if (iterator_to_particle_to_be_removed_in_ensemble !=
2551  ensembles_[0].end())
2552  found_particles_to_remove.push_back(
2553  *iterator_to_particle_to_be_removed_in_ensemble);
2554  }
2555  // Sort the particles found to be removed according to their id and look
2556  // for duplicates (sorting is needed to call std::adjacent_find).
2557  std::sort(found_particles_to_remove.begin(),
2558  found_particles_to_remove.end(),
2559  [](const ParticleData &p1, const ParticleData &p2) {
2560  return p1.id() < p2.id();
2561  });
2562  const auto iterator_to_first_duplicate = std::adjacent_find(
2563  found_particles_to_remove.begin(), found_particles_to_remove.end(),
2564  [](const ParticleData &p1, const ParticleData &p2) {
2565  return p1.id() == p2.id();
2566  });
2567  if (iterator_to_first_duplicate != found_particles_to_remove.end()) {
2568  logg[LExperiment].error() << "The same particle has been asked to be "
2569  "removed multiple times:\n"
2570  << *iterator_to_first_duplicate;
2571  throw std::logic_error("Particle cannot be removed twice!");
2572  }
2573  if (auto delta = remove_plist.size() - found_particles_to_remove.size();
2574  delta > 0) {
2575  logg[LExperiment].warn(
2576  "When trying to remove particle(s) at the beginning ",
2577  "of the system evolution,\n", delta,
2578  " particle(s) could not be found and will be ignored.");
2579  }
2580  if (!found_particles_to_remove.empty()) {
2581  [[maybe_unused]] const auto number_particles_before_removal =
2582  ensembles_[0].size();
2583  // Create and perform action to remove particles
2584  auto action_remove_particles = std::make_unique<FreeforallAction>(
2585  found_particles_to_remove, ParticleList{}, action_time);
2586  perform_action(*action_remove_particles, 0);
2587 
2588  assert(number_particles_before_removal -
2589  found_particles_to_remove.size() ==
2590  ensembles_[0].size());
2591  }
2592  }
2593  }
2594 
2595  if (t_end > end_time_) {
2596  logg[LExperiment].fatal()
2597  << "Evolution asked to be run until " << t_end << " > " << end_time_
2598  << " and this cannot be done (because of how the clock works).";
2599  throw std::logic_error(
2600  "Experiment cannot evolve the system beyond End_Time.");
2601  }
2602  while (*(parameters_.labclock) < t_end) {
2603  const double dt = parameters_.labclock->timestep_duration();
2604  logg[LExperiment].debug("Timestepless propagation for next ", dt, " fm.");
2605 
2606  // Perform forced thermalization if required
2607  if (thermalizer_ &&
2608  thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
2609  const bool ignore_cells_under_treshold = true;
2610  // Thermodynamics in thermalizer is computed from all ensembles,
2611  // but thermalization actions act on each ensemble independently
2612  thermalizer_->update_thermalizer_lattice(ensembles_, density_param_,
2613  ignore_cells_under_treshold);
2614  const double current_t = parameters_.labclock->current_time();
2615  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2616  thermalizer_->thermalize(ensembles_[i_ens], current_t,
2618  ThermalizationAction th_act(*thermalizer_, current_t);
2619  if (th_act.any_particles_thermalized()) {
2620  perform_action(th_act, i_ens);
2621  }
2622  }
2623  }
2624 
2625  if (IC_dynamic_) {
2626  modus_.build_fluidization_lattice(parameters_.labclock->current_time(),
2628  }
2629 
2630  std::vector<Actions> actions(parameters_.n_ensembles);
2631  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2632  actions[i_ens].clear();
2633  if (ensembles_[i_ens].size() > 0 && action_finders_.size() > 0) {
2634  /* (1.a) Create grid. */
2635  const double min_cell_length = compute_min_cell_length(dt);
2636  logg[LExperiment].debug("Creating grid with minimal cell length ",
2637  min_cell_length);
2638  /* For the hyper-surface-crossing actions also unformed particles are
2639  * searched and therefore needed on the grid. */
2640  const bool include_unformed_particles = IC_switch_;
2641  const auto &grid =
2642  use_grid_ ? modus_.create_grid(ensembles_[i_ens], min_cell_length,
2643  dt, parameters_.coll_crit,
2644  include_unformed_particles)
2645  : modus_.create_grid(ensembles_[i_ens], min_cell_length,
2646  dt, parameters_.coll_crit,
2647  include_unformed_particles,
2649 
2650  const double gcell_vol = grid.cell_volume();
2651  /* (1.b) Iterate over cells and find actions. */
2652  grid.iterate_cells(
2653  [&](const ParticleList &search_list) {
2654  for (const auto &finder : action_finders_) {
2655  actions[i_ens].insert(finder->find_actions_in_cell(
2656  search_list, dt, gcell_vol, beam_momentum_));
2657  }
2658  },
2659  [&](const ParticleList &search_list,
2660  const ParticleList &neighbors_list) {
2661  for (const auto &finder : action_finders_) {
2662  actions[i_ens].insert(finder->find_actions_with_neighbors(
2663  search_list, neighbors_list, dt, beam_momentum_));
2664  }
2665  });
2666  }
2667  }
2668 
2669  /* \todo (optimizations) Adapt timestep size here */
2670 
2671  /* (2) Propagate from action to action until next output or timestep end */
2672  const double end_timestep_time = parameters_.labclock->next_time();
2673  while (next_output_time() < end_timestep_time) {
2674  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2675  run_time_evolution_timestepless(actions[i_ens], i_ens,
2676  next_output_time());
2677  }
2678  ++(*parameters_.outputclock);
2679 
2681  }
2682  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2683  run_time_evolution_timestepless(actions[i_ens], i_ens, end_timestep_time);
2684  }
2685 
2686  /* (3) Update potentials (if computed on the lattice) and
2687  * compute new momenta according to equations of motion */
2688  if (potentials_) {
2690  update_momenta(ensembles_, parameters_.labclock->timestep_duration(),
2691  *potentials_, FB_lat_.get(), FI3_lat_.get(), EM_lat_.get(),
2692  jmu_B_lat_.get());
2693  }
2694 
2695  /* (4) Expand universe if non-minkowskian metric; updates
2696  * positions and momenta according to the selected expansion */
2698  for (Particles &particles : ensembles_) {
2699  expand_space_time(&particles, parameters_, metric_);
2700  }
2701  }
2702 
2703  ++(*parameters_.labclock);
2704 
2705  /* (5) Check conservation laws.
2706  *
2707  * Check conservation of conserved quantities if potentials and string
2708  * fragmentation are off. If potentials are on then momentum is conserved
2709  * only in average. If string fragmentation is on, then energy and
2710  * momentum are only very roughly conserved in high-energy collisions. */
2713  std::string err_msg = conserved_initial_.report_deviations(ensembles_);
2714  if (!err_msg.empty()) {
2715  logg[LExperiment].error() << err_msg;
2716  throw std::runtime_error("Violation of conserved quantities!");
2717  }
2718  }
2719  }
2720 
2721  /* Increment once more the output clock in order to have it prepared for the
2722  * final_output() call. Once close to the end time, the while-loop above to
2723  * produce intermediate output is not entered as the next_output_time() is
2724  * never strictly smaller than end_timestep_time (they are usually equal).
2725  * Since in the final_output() function the current time of the output clock
2726  * is used to produce the output, this has to be incremented before producing
2727  * the final output and it makes sense to do it here.
2728  */
2729  ++(*parameters_.outputclock);
2730 
2731  if (pauli_blocker_) {
2732  logg[LExperiment].info(
2733  "Interactions: Pauli-blocked/performed = ", total_pauli_blocked_, "/",
2735  }
2736 }
const ExpansionProperties metric_
This struct contains information on the metric to be used.
Definition: experiment.h:603
std::vector< std::unique_ptr< ActionFinderInterface > > action_finders_
The Action finder objects.
Definition: experiment.h:425
double next_output_time() const
Shortcut for next output time.
Definition: experiment.h:351
const bool IC_dynamic_
This indicates if the IC is dynamic.
Definition: experiment.h:627
std::unique_ptr< GrandCanThermalizer > thermalizer_
Instance of class used for forced thermalization.
Definition: experiment.h:535
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:2765
void intermediate_output()
Intermediate output during an event.
Definition: experiment.h:2829
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FI3_lat_
Lattices for the electric and magnetic component of the symmetry force.
Definition: experiment.h:484
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:343
const bool IC_switch_
This indicates whether the experiment will be used as initial condition for hydrodynamics.
Definition: experiment.h:624
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:480
std::unique_ptr< PauliBlocker > pauli_blocker_
An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and c...
Definition: experiment.h:397
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:600
void update_potentials()
Recompute potentials on lattices if necessary.
Definition: experiment.h:2980
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:131
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:106
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:608
@ 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.
Here is the call graph for this function:

◆ do_final_interactions()

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

Performs the final decays of an event.

Exceptions
runtime_errorif found actions cannot be performed

Definition at line 3081 of file experiment.h.

3081  {
3082  /* At end of time evolution: Force all resonances to decay. In order to handle
3083  * decay chains, we need to loop until no further actions occur. */
3084  bool actions_performed, actions_found;
3085  uint64_t interactions_old;
3086  do {
3087  actions_found = false;
3088  interactions_old = interactions_total_;
3089  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
3090  Actions actions;
3091  // Dileptons: shining of remaining resonances
3092  if (dilepton_finder_ != nullptr) {
3093  for (const auto &output : outputs_) {
3094  dilepton_finder_->shine_final(ensembles_[i_ens], output.get(), true);
3095  }
3096  }
3097  // Find actions.
3098  for (const auto &finder : action_finders_) {
3099  auto found_actions = finder->find_final_actions(ensembles_[i_ens]);
3100  if (!found_actions.empty()) {
3101  actions.insert(std::move(found_actions));
3102  actions_found = true;
3103  }
3104  }
3105  // Perform actions.
3106  while (!actions.is_empty()) {
3107  perform_action(*actions.pop(), i_ens, false);
3108  }
3109  }
3110  actions_performed = interactions_total_ > interactions_old;
3111  // Throw an error if actions were found but not performed
3112  if (actions_found && !actions_performed) {
3113  throw std::runtime_error("Final actions were found but not performed.");
3114  }
3115  // loop until no more decays occur
3116  } while (actions_performed);
3117 
3118  // Dileptons: shining of stable particles at the end
3119  if (dilepton_finder_ != nullptr) {
3120  for (const auto &output : outputs_) {
3121  for (Particles &particles : ensembles_) {
3122  dilepton_finder_->shine_final(particles, output.get(), false);
3123  }
3124  }
3125  }
3126 }
std::unique_ptr< DecayActionsFinderDilepton > dilepton_finder_
The Dilepton Action Finder.
Definition: experiment.h:428
Here is the call graph for this function:

◆ final_output()

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

Output at the end of an event.

Definition at line 3129 of file experiment.h.

3129  {
3130  /* make sure the experiment actually ran (note: we should compare this
3131  * to the start time, but we don't know that. Therefore, we check that
3132  * the time is positive, which should heuristically be the same). */
3133  double E_mean_field = 0.0;
3134  if (likely(parameters_.labclock > 0)) {
3135  const uint64_t wall_actions_this_interval =
3137  const uint64_t interactions_this_interval = interactions_total_ -
3139  wall_actions_this_interval;
3140  if (potentials_) {
3141  // using the lattice is necessary
3142  if ((jmu_B_lat_ != nullptr)) {
3144  EM_lat_.get(), parameters_);
3145  }
3146  }
3147  if (std::abs(parameters_.labclock->current_time() - end_time_) >
3148  really_small) {
3149  logg[LExperiment].warn()
3150  << "SMASH not propagated until configured end time. Current time = "
3151  << parameters_.labclock->current_time()
3152  << "fm. End time = " << end_time_ << "fm.";
3153  } else {
3155  ensembles_, interactions_this_interval, conserved_initial_,
3157  }
3158  int total_particles = 0;
3159  for (const Particles &particles : ensembles_) {
3160  total_particles += particles.size();
3161  }
3162  if (IC_switch_ && (total_particles == 0)) {
3163  const double initial_system_energy_plus_Pythia_violations =
3165  const double fraction_of_total_system_energy_removed =
3166  initial_system_energy_plus_Pythia_violations / total_energy_removed_;
3167  // Verify there is no more energy in the system if all particles were
3168  // removed when crossing the hypersurface
3169  if (std::fabs(fraction_of_total_system_energy_removed - 1.) >
3170  really_small) {
3171  throw std::runtime_error(
3172  "There is remaining energy in the system although all particles "
3173  "were removed.\n"
3174  "E_remain = " +
3175  std::to_string((initial_system_energy_plus_Pythia_violations -
3177  " [GeV]");
3178  } else {
3179  logg[LExperiment].info() << hline;
3180  logg[LExperiment].info()
3181  << "Time real: " << SystemClock::now() - time_start_;
3182  logg[LExperiment].info()
3183  << "Interactions before reaching hypersurface: "
3186  logg[LExperiment].info()
3187  << "Total number of particles removed on hypersurface: "
3189  }
3190  } else {
3191  const double precent_discarded =
3193  ? static_cast<double>(discarded_interactions_total_) * 100.0 /
3195  : 0.0;
3196  std::stringstream msg_discarded;
3197  msg_discarded
3198  << "Discarded interaction number: " << discarded_interactions_total_
3199  << " (" << precent_discarded
3200  << "% of the total interaction number including wall crossings)";
3201 
3202  logg[LExperiment].info() << hline;
3203  logg[LExperiment].info()
3204  << "Time real: " << SystemClock::now() - time_start_;
3205  logg[LExperiment].debug() << msg_discarded.str();
3206 
3208  precent_discarded > 1.0) {
3209  // The chosen threshold of 1% is a heuristical value
3210  logg[LExperiment].warn()
3211  << msg_discarded.str()
3212  << "\nThe number of discarded interactions is large, which means "
3213  "the assumption for the stochastic criterion of\n1 interaction "
3214  "per particle per timestep is probably violated. Consider "
3215  "reducing the timestep size.";
3216  }
3217 
3218  logg[LExperiment].info() << "Final interaction number: "
3220  }
3221 
3222  // Check if there are unformed particles
3223  int unformed_particles_count = 0;
3224  for (const Particles &particles : ensembles_) {
3225  for (const ParticleData &particle : particles) {
3226  if (particle.formation_time() > end_time_) {
3227  unformed_particles_count++;
3228  }
3229  }
3230  }
3231  if (unformed_particles_count > 0) {
3232  logg[LExperiment].warn(
3233  "End time might be too small. ", unformed_particles_count,
3234  " unformed particles were found at the end of the evolution.");
3235  }
3236  }
3237 
3238  // Keep track of how many ensembles had interactions
3240 
3241  for (const auto &output : outputs_) {
3242  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
3243  auto event_info = fill_event_info(
3244  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
3246  output->at_eventend(ensembles_[i_ens], {event_, i_ens}, event_info);
3247  }
3248  // For thermodynamic output
3249  output->at_eventend(ensembles_, event_);
3250 
3251  // For thermodynamic lattice output
3252  if (printout_rho_eckart_) {
3255  }
3256  }
3257  if (printout_tmn_) {
3258  output->at_eventend(ThermodynamicQuantity::Tmn);
3259  }
3260  if (printout_tmn_landau_) {
3262  }
3263  if (printout_v_landau_) {
3265  }
3266  if (printout_j_QBS_) {
3267  output->at_eventend(ThermodynamicQuantity::j_QBS);
3268  }
3269  }
3270 }
void count_nonempty_ensembles()
Counts the number of ensembles in wich interactions took place at the end of an event.
Definition: experiment.h:3273
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
std::string to_string(ThermodynamicQuantity quantity)
Convert a ThermodynamicQuantity enum value to its corresponding string.
Definition: stringify.cc:26
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:41
Here is the call graph for this function:

◆ 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 262 of file experiment.h.

262 { return &ensembles_[0]; }

◆ all_ensembles()

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

Getter for all ensembles.

Definition at line 264 of file experiment.h.

264 { 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 270 of file experiment.h.

270 { 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 3308 of file experiment.h.

3308  {
3309  event_++;
3310 }

◆ 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 734 of file experiment.h.

737  {
738  // Disable output which do not properly work with multiple ensembles
739  if (ensembles_.size() > 1) {
740  auto abort_because_of = [](const std::string &s) {
741  throw std::invalid_argument(
742  s + " output is not available with multiple parallel ensembles.");
743  };
744  if (content == "Initial_Conditions") {
745  abort_because_of("Initial_Conditions");
746  }
747  if ((format == "HepMC") || (format == "HepMC_asciiv3") ||
748  (format == "HepMC_treeroot")) {
749  abort_because_of("HepMC");
750  }
751  if (content == "Rivet") {
752  abort_because_of("Rivet");
753  }
754  if (content == "Collisions") {
755  logg[LExperiment].warn(
756  "Information coming from different ensembles in 'Collisions' output "
757  "is not distinguishable.\nSuch an output with multiple parallel "
758  "ensembles should only be used if later in the data analysis\nit is "
759  "not necessary to trace back which data belongs to which ensemble.");
760  }
761  }
762 
763  if (format == "VTK" && content == "Particles") {
764  outputs_.emplace_back(
765  std::make_unique<VtkOutput>(output_path, content, out_par));
766  } else if (format == "Root") {
767 #ifdef SMASH_USE_ROOT
768  if (content == "Initial_Conditions") {
769  outputs_.emplace_back(
770  std::make_unique<RootOutput>(output_path, "SMASH_IC", out_par));
771  } else {
772  outputs_.emplace_back(
773  std::make_unique<RootOutput>(output_path, content, out_par));
774  }
775 #else
776  logg[LExperiment].error(
777  "Root output requested, but Root support not compiled in");
778 #endif
779  } else if ((format == "Binary" || format == "Oscar2013_bin") &&
780  (content == "Collisions" || content == "Particles" ||
781  content == "Dileptons" || content == "Photons" ||
782  content == "Initial_Conditions")) {
783  outputs_.emplace_back(
784  create_binary_output(format, content, output_path, out_par));
785  } else if (format == "Oscar1999" || format == "Oscar2013") {
786  outputs_.emplace_back(
787  create_oscar_output(format, content, output_path, out_par));
788  } else if (format == "ASCII" &&
789  (content == "Particles" || content == "Collisions" ||
790  content == "Dileptons" || content == "Photons" ||
791  content == "Initial_Conditions")) {
792  outputs_.emplace_back(
793  create_oscar_output(format, content, output_path, out_par));
794  } else if (content == "Thermodynamics" && format == "ASCII") {
795  outputs_.emplace_back(
796  std::make_unique<ThermodynamicOutput>(output_path, content, out_par));
797  } else if (content == "Thermodynamics" &&
798  (format == "Lattice_ASCII" || format == "Lattice_Binary")) {
800  outputs_.emplace_back(std::make_unique<ThermodynamicLatticeOutput>(
801  output_path, content, out_par, format == "Lattice_ASCII",
802  format == "Lattice_Binary"));
803  } else if (content == "Thermodynamics" && format == "VTK") {
804  printout_lattice_td_ = true;
805  outputs_.emplace_back(
806  std::make_unique<VtkOutput>(output_path, content, out_par));
807  } else if (content == "Initial_Conditions" && format == "For_vHLLE") {
808  if (IC_dynamic_) {
809  throw std::invalid_argument(
810  "Dynamic initial conditions are only available in Oscar2013 and "
811  "Binary formats.");
812  }
813  outputs_.emplace_back(
814  std::make_unique<ICOutput>(output_path, "SMASH_IC_For_vHLLE", out_par));
815  } else if ((format == "HepMC") || (format == "HepMC_asciiv3") ||
816  (format == "HepMC_treeroot")) {
817 #ifdef SMASH_USE_HEPMC
818  if (content == "Particles") {
819  if ((format == "HepMC") || (format == "HepMC_asciiv3")) {
820  outputs_.emplace_back(std::make_unique<HepMcOutput>(
821  output_path, "SMASH_HepMC_particles", false, "asciiv3"));
822  } else if (format == "HepMC_treeroot") {
823 #ifdef SMASH_USE_HEPMC_ROOTIO
824  outputs_.emplace_back(std::make_unique<HepMcOutput>(
825  output_path, "SMASH_HepMC_particles", false, "root"));
826 #else
827  logg[LExperiment].error(
828  "Requested HepMC_treeroot output not available, "
829  "ROOT or HepMC3_ROOTIO missing or not found by cmake.");
830 #endif
831  }
832  } else if (content == "Collisions") {
833  if ((format == "HepMC") || (format == "HepMC_asciiv3")) {
834  outputs_.emplace_back(std::make_unique<HepMcOutput>(
835  output_path, "SMASH_HepMC_collisions", true, "asciiv3"));
836  } else if (format == "HepMC_treeroot") {
837 #ifdef SMASH_USE_HEPMC_ROOTIO
838  outputs_.emplace_back(std::make_unique<HepMcOutput>(
839  output_path, "SMASH_HepMC_collisions", true, "root"));
840 #else
841  logg[LExperiment].error(
842  "Requested HepMC_treeroot output not available, "
843  "ROOT or HepMC3_ROOTIO missing or not found by cmake.");
844 #endif
845  }
846  } else {
847  logg[LExperiment].error(
848  "HepMC only available for Particles and "
849  "Collisions content. Requested for " +
850  content + ".");
851  }
852 #else
853  logg[LExperiment].error(
854  "HepMC output requested, but HepMC support not compiled in");
855 #endif
856  } else if (content == "Coulomb" && format == "VTK") {
857  printout_coulomb_vtk_ = true;
858  outputs_.emplace_back(
859  std::make_unique<VtkOutput>(output_path, "Fields", out_par));
860  } else if (content == "Rivet") {
861 #ifdef SMASH_USE_RIVET
862  // flag to ensure that the Rivet format has not been already assigned
863  static bool rivet_format_already_selected = false;
864  // if the next check is true, then we are trying to assign the format twice
865  if (rivet_format_already_selected) {
866  logg[LExperiment].warn(
867  "Rivet output format can only be one, either YODA or YODA-full. "
868  "Only your first valid choice will be used.");
869  return;
870  }
871  if (format == "YODA") {
872  outputs_.emplace_back(std::make_unique<RivetOutput>(
873  output_path, "SMASH_Rivet", false, out_par.rivet_parameters));
874  rivet_format_already_selected = true;
875  } else if (format == "YODA-full") {
876  outputs_.emplace_back(std::make_unique<RivetOutput>(
877  output_path, "SMASH_Rivet_full", true, out_par.rivet_parameters));
878  rivet_format_already_selected = true;
879  } else {
880  logg[LExperiment].error("Rivet format " + format +
881  "not one of YODA or YODA-full");
882  }
883 #else
884  logg[LExperiment].error(
885  "Rivet output requested, but Rivet support not compiled in");
886 #endif
887  } else {
888  logg[LExperiment].error()
889  << "Unknown combination of format (" << format << ") and content ("
890  << content << "). Fix the config.";
891  }
892 
893  logg[LExperiment].info() << "Added output " << content << " of format "
894  << format << "\n";
895 }
bool printout_coulomb_vtk_
Whether to write the electric and magnetic fields to VTK files.
Definition: experiment.h:532
bool printout_lattice_td_
Whether to print the thermodynamics quantities evaluated on the lattices.
Definition: experiment.h:525
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:926
std::unique_ptr< OutputInterface > create_binary_output(const std::string &format, const std::string &content, const std::filesystem::path &path, const OutputParameters &out_par)
Create a binary output object.
Here is the call graph for this function:

◆ 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 2739 of file experiment.h.

2740  {
2741  const double dt =
2742  propagate_straight_line(&particles, to_time, beam_momentum_);
2743  if (dilepton_finder_ != nullptr) {
2744  for (const auto &output : outputs_) {
2745  dilepton_finder_->shine(particles, output.get(), dt);
2746  }
2747  }
2748 }
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
Here is the call graph for this function:

◆ 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 2765 of file experiment.h.

2766  {
2767  Particles &particles = ensembles_[i_ensemble];
2768  logg[LExperiment].debug(
2769  "Timestepless propagation: ", "Actions size = ", actions.size(),
2770  ", end time = ", end_time_propagation);
2771 
2772  // iterate over all actions
2773  while (!actions.is_empty()) {
2774  if (actions.earliest_time() > end_time_propagation) {
2775  break;
2776  }
2777  // get next action
2778  ActionPtr act = actions.pop();
2779  if (!act->is_valid(particles)) {
2781  logg[LExperiment].debug(~einhard::DRed(), "✘ ", act,
2782  " (discarded: invalid)");
2783  continue;
2784  }
2785  logg[LExperiment].debug(~einhard::Green(), "✔ ", act,
2786  ", action time = ", act->time_of_execution());
2787 
2788  /* (1) Propagate to the next action. */
2789  propagate_and_shine(act->time_of_execution(), particles);
2790 
2791  /* (2) Perform action.
2792  *
2793  * Update the positions of the incoming particles, because the information
2794  * in the action object will be outdated as the particles have been
2795  * propagated since the construction of the action. */
2796  act->update_incoming(particles);
2797  const bool performed = perform_action(*act, i_ensemble);
2798 
2799  /* No need to update actions for outgoing particles
2800  * if the action is not performed. */
2801  if (!performed) {
2802  continue;
2803  }
2804 
2805  /* (3) Update actions for newly-produced particles. */
2806 
2807  const double end_time_timestep = parameters_.labclock->next_time();
2808  // New actions are always search until the end of the current timestep
2809  const double time_left = end_time_timestep - act->time_of_execution();
2810  const ParticleList &outgoing_particles = act->outgoing_particles();
2811  // Grid cell volume set to zero, since there is no grid
2812  const double gcell_vol = 0.0;
2813  for (const auto &finder : action_finders_) {
2814  // Outgoing particles can still decay, cross walls...
2815  actions.insert(finder->find_actions_in_cell(outgoing_particles, time_left,
2816  gcell_vol, beam_momentum_));
2817  // ... and collide with other particles.
2818  actions.insert(finder->find_actions_with_surrounding_particles(
2819  outgoing_particles, particles, time_left, beam_momentum_));
2820  }
2821 
2823  }
2824 
2825  propagate_and_shine(end_time_propagation, particles);
2826 }
A stream modifier that allows to colorize the log output.
Definition: einhard.hpp:147
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:2739
void check_interactions_total(uint64_t interactions_total)
Make sure interactions_total can be represented as a 32-bit integer.
Definition: experiment.h:2757
Here is the call graph for this function:

◆ 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 2829 of file experiment.h.

2829  {
2830  const uint64_t wall_actions_this_interval =
2833  const uint64_t interactions_this_interval = interactions_total_ -
2835  wall_actions_this_interval;
2837  double E_mean_field = 0.0;
2838  /// Auxiliary variable to communicate the time in the computational frame
2839  /// at the functions printing the thermodynamics lattice output
2840  double computational_frame_time = 0.0;
2841  if (potentials_) {
2842  // using the lattice is necessary
2843  if ((jmu_B_lat_ != nullptr)) {
2845  EM_lat_.get(), parameters_);
2846  /*
2847  * Mean field calculated in a box should remain approximately constant if
2848  * the system is in equilibrium, and so deviations from its original value
2849  * may signal a phase transition or other dynamical process. This
2850  * comparison only makes sense in the Box Modus, hence the condition.
2851  */
2852  if (modus_.is_box()) {
2853  double tmp = (E_mean_field - initial_mean_field_energy_) /
2854  (E_mean_field + initial_mean_field_energy_);
2855  /*
2856  * This is displayed when the system evolves away from its initial
2857  * configuration (which is when the total mean field energy in the box
2858  * deviates from its initial value).
2859  */
2860  if (std::abs(tmp) > 0.01) {
2861  logg[LExperiment].info()
2862  << "\n\n\n\t The mean field at t = "
2863  << parameters_.outputclock->current_time()
2864  << " [fm] differs from the mean field at t = 0:"
2865  << "\n\t\t initial_mean_field_energy_ = "
2866  << initial_mean_field_energy_ << " [GeV]"
2867  << "\n\t\t abs[(E_MF - E_MF(t=0))/(E_MF + E_MF(t=0))] = "
2868  << std::abs(tmp)
2869  << "\n\t\t E_MF/E_MF(t=0) = "
2870  << E_mean_field / initial_mean_field_energy_ << "\n\n";
2871  }
2872  }
2873  }
2874  }
2875 
2877  ensembles_, interactions_this_interval, conserved_initial_, time_start_,
2878  parameters_.outputclock->current_time(), E_mean_field,
2880  const LatticeUpdate lat_upd = LatticeUpdate::AtOutput;
2881 
2882  // save evolution data
2883  if (!(modus_.is_box() && parameters_.outputclock->current_time() <
2884  modus_.equilibration_time())) {
2885  for (const auto &output : outputs_) {
2886  if (output->is_dilepton_output() || output->is_photon_output() ||
2887  output->is_IC_output()) {
2888  continue;
2889  }
2890  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2891  auto event_info = fill_event_info(
2892  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2894 
2895  output->at_intermediate_time(ensembles_[i_ens], parameters_.outputclock,
2896  density_param_, {event_, i_ens},
2897  event_info);
2898  computational_frame_time = event_info.current_time;
2899  }
2900  // For thermodynamic output
2901  output->at_intermediate_time(ensembles_, parameters_.outputclock,
2902  density_param_);
2903 
2904  // Thermodynamic output on the lattice versus time
2905  if (printout_rho_eckart_) {
2906  switch (dens_type_lattice_printout_) {
2907  case DensityType::Baryon:
2910  ensembles_, false);
2911  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2913  output->thermodynamics_lattice_output(*jmu_B_lat_,
2914  computational_frame_time);
2915  break;
2918  jmu_I3_lat_.get(), lat_upd, DensityType::BaryonicIsospin,
2919  density_param_, ensembles_, false);
2920  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2922  *jmu_I3_lat_);
2923  output->thermodynamics_lattice_output(*jmu_I3_lat_,
2924  computational_frame_time);
2925  break;
2926  case DensityType::None:
2927  break;
2928  default:
2931  density_param_, ensembles_, false);
2932  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2934  *jmu_custom_lat_);
2935  output->thermodynamics_lattice_output(*jmu_custom_lat_,
2936  computational_frame_time);
2937  }
2938  }
2942  ensembles_, false);
2943  if (printout_tmn_) {
2944  output->thermodynamics_output(ThermodynamicQuantity::Tmn,
2946  output->thermodynamics_lattice_output(
2947  ThermodynamicQuantity::Tmn, *Tmn_, computational_frame_time);
2948  }
2949  if (printout_tmn_landau_) {
2950  output->thermodynamics_output(ThermodynamicQuantity::TmnLandau,
2952  output->thermodynamics_lattice_output(
2954  computational_frame_time);
2955  }
2956  if (printout_v_landau_) {
2957  output->thermodynamics_output(ThermodynamicQuantity::LandauVelocity,
2959  output->thermodynamics_lattice_output(
2961  computational_frame_time);
2962  }
2963  }
2964  if (EM_lat_) {
2965  output->fields_output("Efield", "Bfield", *EM_lat_);
2966  }
2967  if (printout_j_QBS_) {
2968  output->thermodynamics_lattice_output(
2969  *j_QBS_lat_, computational_frame_time, ensembles_, density_param_);
2970  }
2971 
2972  if (thermalizer_) {
2973  output->thermodynamics_output(*thermalizer_);
2974  }
2975  }
2976  }
2977 }
void update_lattice_accumulating_ensembles(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 when ensembles are used.
Definition: density.h:654
LatticeUpdate
Enumerator option for lattice updates.
Definition: lattice.h:38
Here is the call graph for this function:

◆ update_potentials()

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

Recompute potentials on lattices if necessary.

Definition at line 2980 of file experiment.h.

2980  {
2981  if (potentials_) {
2982  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
2987  parameters_.labclock->timestep_duration(), true);
2988  }
2989  if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
2990  jmu_B_lat_ != nullptr) {
2995  parameters_.labclock->timestep_duration(), true);
2996  const size_t UBlattice_size = UB_lat_->size();
2997  for (size_t i = 0; i < UBlattice_size; i++) {
2998  auto jB = (*jmu_B_lat_)[i];
2999  const FourVector flow_four_velocity_B =
3000  std::abs(jB.rho()) > very_small_double ? jB.jmu_net() / jB.rho()
3001  : FourVector();
3002  double baryon_density = jB.rho();
3003  ThreeVector baryon_grad_j0 = jB.grad_j0();
3004  ThreeVector baryon_dvecj_dt = jB.dvecj_dt();
3005  ThreeVector baryon_curl_vecj = jB.curl_vecj();
3006  if (potentials_->use_skyrme()) {
3007  (*UB_lat_)[i] =
3008  flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
3009  (*FB_lat_)[i] =
3010  potentials_->skyrme_force(baryon_density, baryon_grad_j0,
3011  baryon_dvecj_dt, baryon_curl_vecj);
3012  }
3013  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
3014  auto jI3 = (*jmu_I3_lat_)[i];
3015  const FourVector flow_four_velocity_I3 =
3016  std::abs(jI3.rho()) > very_small_double
3017  ? jI3.jmu_net() / jI3.rho()
3018  : FourVector();
3019  (*UI3_lat_)[i] = flow_four_velocity_I3 *
3020  potentials_->symmetry_pot(jI3.rho(), baryon_density);
3021  (*FI3_lat_)[i] = potentials_->symmetry_force(
3022  jI3.rho(), jI3.grad_j0(), jI3.dvecj_dt(), jI3.curl_vecj(),
3023  baryon_density, baryon_grad_j0, baryon_dvecj_dt,
3024  baryon_curl_vecj);
3025  }
3026  }
3027  }
3028  if (potentials_->use_coulomb()) {
3031  density_param_, ensembles_, true);
3032  for (size_t i = 0; i < EM_lat_->size(); i++) {
3033  ThreeVector electric_field = {0., 0., 0.};
3034  ThreeVector position = jmu_el_lat_->cell_center(i);
3035  jmu_el_lat_->integrate_volume(electric_field,
3037  potentials_->coulomb_r_cut(), position);
3038  ThreeVector magnetic_field = {0., 0., 0.};
3039  jmu_el_lat_->integrate_volume(magnetic_field,
3041  potentials_->coulomb_r_cut(), position);
3042  (*EM_lat_)[i] = std::make_pair(electric_field, magnetic_field);
3043  }
3044  } // if ((potentials_->use_skyrme() || ...
3045  if (potentials_->use_vdf() && jmu_B_lat_ != nullptr) {
3050  parameters_.labclock->timestep_duration(), true);
3053  fields_lat_.get(), old_fields_auxiliary_.get(),
3056  parameters_.labclock->timestep_duration());
3057  }
3058  const size_t UBlattice_size = UB_lat_->size();
3059  for (size_t i = 0; i < UBlattice_size; i++) {
3060  auto jB = (*jmu_B_lat_)[i];
3061  (*UB_lat_)[i] = potentials_->vdf_pot(jB.rho(), jB.jmu_net());
3064  (*FB_lat_)[i] = potentials_->vdf_force(
3065  jB.rho(), jB.drho_dxnu().x0(), jB.drho_dxnu().threevec(),
3066  jB.grad_rho_cross_vecj(), jB.jmu_net().x0(), jB.grad_j0(),
3067  jB.jmu_net().threevec(), jB.dvecj_dt(), jB.curl_vecj());
3068  break;
3070  auto Amu = (*fields_lat_)[i];
3071  (*FB_lat_)[i] = potentials_->vdf_force(
3072  Amu.grad_A0(), Amu.dvecA_dt(), Amu.curl_vecA());
3073  break;
3074  }
3075  } // for (size_t i = 0; i < UBlattice_size; i++)
3076  } // if potentials_->use_vdf()
3077  }
3078 }
std::unique_ptr< RectangularLattice< FourVector > > new_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0 + dt.
Definition: experiment.h:504
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > fields_four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of Amu.
Definition: experiment.h:507
std::unique_ptr< FieldsLattice > fields_lat_
Mean-field A^mu on the lattice.
Definition: experiment.h:449
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:467
std::unique_ptr< RectangularLattice< FourVector > > old_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0.
Definition: experiment.h:502
std::unique_ptr< DensityLattice > jmu_el_lat_
Electric charge density on the lattice.
Definition: experiment.h:446
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:383
static ThreeVector E_field_integrand(ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
Integrand for calculating the electric field.
Definition: potentials.h:365
constexpr double very_small_double
A very small double, used to avoid division by zero.
Definition: constants.h:44
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
Here is the call graph for this function:

◆ 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 343 of file experiment.h.

343  {
346  }
347  return std::sqrt(4 * dt * dt + max_transverse_distance_sqr_);
348  }
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 351 of file experiment.h.

351  {
352  return parameters_.outputclock->next_time();
353  }

◆ 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 3273 of file experiment.h.

3273  {
3274  for (bool has_interaction : projectile_target_interact_) {
3275  if (has_interaction) {
3277  }
3278  }
3279 }
int nonempty_ensembles_
Number of ensembles containing an interaction.
Definition: experiment.h:581

◆ 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 3282 of file experiment.h.

3282  {
3284  return event_ >= nevents_;
3285  }
3288  return true;
3289  }
3290  if (event_ >= max_events_) {
3291  logg[LExperiment].warn()
3292  << "Maximum number of events (" << max_events_
3293  << ") exceeded. Stopping calculation. "
3294  << "The fraction of empty ensembles is "
3295  << (1.0 - static_cast<double>(nonempty_ensembles_) /
3297  << ". If this fraction is expected, try increasing the "
3298  "Maximum_Ensembles_Run.";
3299  return true;
3300  }
3301  return false;
3302  }
3303  throw std::runtime_error("Event counting option is invalid");
3304  return false;
3305 }
int minimum_nonempty_ensembles_
The number of ensembles, in which interactions take place, to be calculated.
Definition: experiment.h:567
int max_events_
Maximum number of events to be calculated in order obtain the desired number of non-empty events usin...
Definition: experiment.h:587
EventCounting event_counting_
The way in which the number of calculated events is specified.
Definition: experiment.h:575
int nevents_
Number of events.
Definition: experiment.h:557
@ 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 194 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 373 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 376 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 382 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 385 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 391 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 397 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 403 of file experiment.h.

◆ dilepton_output_

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

The Dilepton output.

Definition at line 406 of file experiment.h.

◆ photon_output_

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

The Photon output.

Definition at line 409 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 415 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 422 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 425 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 428 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 431 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 434 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 437 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 440 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 443 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 446 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 449 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 458 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 461 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 467 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 473 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 480 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 484 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 488 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 491 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 494 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 496 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 499 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 502 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 504 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 507 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 510 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 513 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 516 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 519 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 522 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 525 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 529 of file experiment.h.

◆ printout_coulomb_vtk_

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

Whether to write the electric and magnetic fields to VTK files.

Definition at line 532 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 535 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 541 of file experiment.h.

◆ nevents_

template<typename Modus >
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 557 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 567 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 575 of file experiment.h.

◆ event_

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

Current event.

Definition at line 578 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 581 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 587 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 590 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 597 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 600 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 603 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 606 of file experiment.h.

◆ dileptons_bremsstrahlung_switch_

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

This indicates whether dilepton production via bremsstrahlung is switched on.

Definition at line 612 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 615 of file experiment.h.

◆ photons_bremsstrahlung_switch_

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

This indicates whether bremsstrahlung is switched on.

Definition at line 618 of file experiment.h.

◆ IC_switch_

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

This indicates whether the experiment will be used as initial condition for hydrodynamics.

Currently only the Collider modus can achieve this.

Definition at line 624 of file experiment.h.

◆ IC_dynamic_

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

This indicates if the IC is dynamic.

Definition at line 627 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 630 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 636 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 645 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 651 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 654 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 657 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 663 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 669 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 675 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 681 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 687 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 693 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 699 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 704 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 709 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 712 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 715 of file experiment.h.


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