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

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

Public Member Functions

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

Private Member Functions

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

Private Attributes

ExperimentParameters parameters_
 Struct of several member variables. More...
 
DensityParameters density_param_
 Structure to precalculate and hold parameters for density computations. More...
 
Modus modus_
 Instance of the Modus template parameter. More...
 
std::vector< Particlesensembles_
 Complete particle list, all ensembles in one vector. More...
 
std::unique_ptr< Potentialspotentials_
 An instance of potentials class, that stores parameters of potentials, calculates them and their gradients. More...
 
std::unique_ptr< PauliBlockerpauli_blocker_
 An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and computes phase-space density. More...
 
OutputsList outputs_
 A list of output formaters. More...
 
OutputPtr dilepton_output_
 The Dilepton output. More...
 
OutputPtr photon_output_
 The Photon output. More...
 
std::vector< bool > projectile_target_interact_
 Whether the projectile and the target collided. More...
 
std::vector< FourVectorbeam_momentum_ = {}
 The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen. More...
 
std::vector< std::unique_ptr< ActionFinderInterface > > action_finders_
 The Action finder objects. More...
 
std::unique_ptr< DecayActionsFinderDileptondilepton_finder_
 The Dilepton Action Finder. More...
 
std::unique_ptr< ActionFinderInterfacephoton_finder_
 The (Scatter) Actions Finder for Direct Photons. More...
 
int n_fractional_photons_
 Number of fractional photons produced per single reaction. More...
 
std::unique_ptr< DensityLatticej_QBS_lat_
 4-current for j_QBS lattice output More...
 
std::unique_ptr< DensityLatticejmu_B_lat_
 Baryon density on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_I3_lat_
 Isospin projection density on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_el_lat_
 Electric charge density on the lattice. More...
 
std::unique_ptr< FieldsLatticefields_lat_
 Mean-field A^mu on the lattice. More...
 
std::unique_ptr< DensityLatticejmu_custom_lat_
 Custom density on the lattices. More...
 
DensityType dens_type_lattice_printout_ = DensityType::None
 Type of density for lattice printout. More...
 
std::unique_ptr< RectangularLattice< FourVector > > UB_lat_ = nullptr
 Lattices for Skyrme or VDF potentials (evaluated in the local rest frame) times the baryon flow 4-velocity. More...
 
std::unique_ptr< RectangularLattice< FourVector > > UI3_lat_ = nullptr
 Lattices for symmetry potentials (evaluated in the local rest frame) times the isospin flow 4-velocity. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FB_lat_
 Lattices for the electric and magnetic components of the Skyrme or VDF force. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FI3_lat_
 Lattices for the electric and magnetic component of the symmetry force. More...
 
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > EM_lat_
 Lattices for electric and magnetic field in fm^-2. More...
 
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > > Tmn_
 Lattices of energy-momentum tensors for printout. More...
 
std::unique_ptr< RectangularLattice< FourVector > > old_jmu_auxiliary_
 Auxiliary lattice for values of jmu at a time step t0. More...
 
std::unique_ptr< RectangularLattice< FourVector > > new_jmu_auxiliary_
 Auxiliary lattice for values of jmu at a time step t0 + dt. More...
 
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > four_gradient_auxiliary_
 Auxiliary lattice for calculating the four-gradient of jmu. More...
 
std::unique_ptr< RectangularLattice< FourVector > > old_fields_auxiliary_
 Auxiliary lattice for values of Amu at a time step t0. More...
 
std::unique_ptr< RectangularLattice< FourVector > > new_fields_auxiliary_
 Auxiliary lattice for values of Amu at a time step t0 + dt. More...
 
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > fields_four_gradient_auxiliary_
 Auxiliary lattice for calculating the four-gradient of Amu. More...
 
bool printout_rho_eckart_ = false
 Whether to print the Eckart rest frame density. More...
 
bool printout_tmn_ = false
 Whether to print the energy-momentum tensor. More...
 
bool printout_tmn_landau_ = false
 Whether to print the energy-momentum tensor in Landau frame. More...
 
bool printout_v_landau_ = false
 Whether to print the 4-velocity in Landau frame. More...
 
bool printout_j_QBS_ = false
 Whether to print the Q, B, S 4-currents. More...
 
bool printout_lattice_td_ = false
 Whether to print the thermodynamics quantities evaluated on the lattices. More...
 
bool printout_full_lattice_any_td_ = false
 Whether to print the thermodynamics quantities evaluated on the lattices, point by point, in any format. More...
 
std::unique_ptr< GrandCanThermalizerthermalizer_
 Instance of class used for forced thermalization. More...
 
StringProcessprocess_string_ptr_
 Pointer to the string process class object, which is used to set the random seed for PYTHIA objects in each event. More...
 
int nevents_ = 0
 Number of events. More...
 
int minimum_nonempty_ensembles_ = 0
 The number of ensembles, in which interactions take place, to be calculated. More...
 
EventCounting event_counting_ = EventCounting::Invalid
 The way in which the number of calculated events is specified. More...
 
int event_ = 0
 Current event. More...
 
int nonempty_ensembles_ = 0
 Number of ensembles containing an interaction. More...
 
int max_events_ = 0
 Maximum number of events to be calculated in order obtain the desired number of non-empty events using the MinimumNonemptyEnsembles option. More...
 
const double end_time_
 simulation time at which the evolution is stopped. More...
 
const double delta_time_startup_
 The clock's timestep size at start up. More...
 
const bool force_decays_
 This indicates whether we force all resonances to decay in the last timestep. More...
 
const bool use_grid_
 This indicates whether to use the grid. More...
 
const ExpansionProperties metric_
 This struct contains information on the metric to be used. More...
 
const bool dileptons_switch_
 This indicates whether dileptons are switched on. More...
 
const bool photons_switch_
 This indicates whether photons are switched on. More...
 
const bool bremsstrahlung_switch_
 This indicates whether bremsstrahlung is switched on. More...
 
const bool IC_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 3239 of file experiment.h.

3239  {
3240  const auto &mainlog = logg[LMain];
3241  for (event_ = 0; !is_finished(); event_++) {
3242  mainlog.info() << "Event " << event_;
3243 
3244  // Sample initial particles, start clock, some printout and book-keeping
3246 
3248 
3249  if (force_decays_) {
3250  do_final_decays();
3251  }
3252 
3253  // Output at event end
3254  final_output();
3255  }
3256 }
const bool force_decays_
This indicates whether we force all resonances to decay in the last timestep.
Definition: experiment.h:598
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:2439
void initialize_new_event()
This is called in the beginning of each event.
Definition: experiment.h:2071
bool is_finished()
Checks wether the desired number events have been calculated.
Definition: experiment.h:3208
int event_
Current event.
Definition: experiment.h:573
void do_final_decays()
Performs the final decays of an event.
Definition: experiment.h:3006
void final_output()
Output at the end of an event.
Definition: experiment.h:3055
const double end_time_
simulation time at which the evolution is stopped.
Definition: experiment.h:585
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:40
static constexpr int LMain
Definition: experiment.h:92

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

2071  {
2073  logg[LExperiment].info() << "random number seed: " << seed_;
2074  /* Set seed for the next event. It has to be positive, so it can be entered
2075  * in the config.
2076  *
2077  * We have to be careful about the minimal integer, whose absolute value
2078  * cannot be represented. */
2079  int64_t r = random::advance();
2080  while (r == INT64_MIN) {
2081  r = random::advance();
2082  }
2083  seed_ = std::abs(r);
2084  /* Set the random seed used in PYTHIA hadronization
2085  * to be same with the SMASH one.
2086  * In this way we ensure that the results are reproducible
2087  * for every event if one knows SMASH random seed. */
2088  if (process_string_ptr_ != NULL) {
2090  }
2091 
2092  for (Particles &particles : ensembles_) {
2093  particles.reset();
2094  }
2095 
2096  // Sample particles according to the initial conditions
2097  double start_time = -1.0;
2098 
2099  // Sample impact parameter only once per all ensembles
2100  // It should be the same for all ensembles
2101  if (modus_.is_collider()) {
2102  modus_.sample_impact();
2103  logg[LExperiment].info("Impact parameter = ", modus_.impact_parameter(),
2104  " fm");
2105  }
2106  for (Particles &particles : ensembles_) {
2107  start_time = modus_.initial_conditions(&particles, parameters_);
2108  }
2109  /* For box modus make sure that particles are in the box. In principle, after
2110  * a correct initialization they should be, so this is just playing it safe.
2111  */
2112  for (Particles &particles : ensembles_) {
2113  modus_.impose_boundary_conditions(&particles, outputs_);
2114  }
2115  // Reset the simulation clock
2116  double timestep = delta_time_startup_;
2117 
2118  switch (time_step_mode_) {
2119  case TimeStepMode::Fixed:
2120  break;
2121  case TimeStepMode::None:
2122  timestep = end_time_ - start_time;
2123  // Take care of the box modus + timestepless propagation
2124  const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
2125  if (max_dt > 0. && max_dt < timestep) {
2126  timestep = max_dt;
2127  }
2128  break;
2129  }
2130  std::unique_ptr<UniformClock> clock_for_this_event;
2131  if (modus_.is_list() && (timestep < 0.0)) {
2132  throw std::runtime_error(
2133  "Timestep for the given event is negative. \n"
2134  "This might happen if the formation times of the input particles are "
2135  "larger than the specified end time of the simulation.");
2136  }
2137  clock_for_this_event =
2138  std::make_unique<UniformClock>(start_time, timestep, end_time_);
2139  parameters_.labclock = std::move(clock_for_this_event);
2140 
2141  // Reset the output clock
2142  parameters_.outputclock->reset(start_time, true);
2143  // remove time before starting time in case of custom output times.
2144  parameters_.outputclock->remove_times_in_past(start_time);
2145 
2146  logg[LExperiment].debug(
2147  "Lab clock: t_start = ", parameters_.labclock->current_time(),
2148  ", dt = ", parameters_.labclock->timestep_duration());
2149 
2150  /* Save the initial conserved quantum numbers and total momentum in
2151  * the system for conservation checks */
2152  conserved_initial_ = QuantumNumbers(ensembles_);
2153  wall_actions_total_ = 0;
2155  interactions_total_ = 0;
2161  total_energy_removed_ = 0.0;
2163  // Print output headers
2164  logg[LExperiment].info() << hline;
2165  logg[LExperiment].info() << "Time[fm] Ekin[GeV] E_MF[GeV] ETotal[GeV] "
2166  << "ETot/N[GeV] D(ETot/N)[GeV] Scatt&Decays "
2167  << "Particles Comp.Time";
2168  logg[LExperiment].info() << hline;
2169  double E_mean_field = 0.0;
2170  if (potentials_) {
2171  // update_potentials();
2172  // if (parameters.outputclock->current_time() == 0.0 )
2173  // using the lattice is necessary
2174  if ((jmu_B_lat_ != nullptr)) {
2179  parameters_.labclock->timestep_duration(), true);
2180  // Because there was no lattice at t=-Delta_t, the time derivatives
2181  // drho_dt and dj^mu/dt at t=0 are huge, while they shouldn't be; we
2182  // overwrite the time derivative to zero by hand.
2183  for (auto &node : *jmu_B_lat_) {
2184  node.overwrite_drho_dt_to_zero();
2185  node.overwrite_djmu_dt_to_zero();
2186  }
2188  EM_lat_.get(), parameters_);
2189  }
2190  }
2191  initial_mean_field_energy_ = E_mean_field;
2194  parameters_.labclock->current_time(), E_mean_field,
2196 
2197  // Output at event start
2198  for (const auto &output : outputs_) {
2199  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2200  auto event_info = fill_event_info(
2201  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2203  output->at_eventstart(ensembles_[i_ens], {event_, i_ens}, event_info);
2204  }
2205  // For thermodynamic output
2206  output->at_eventstart(ensembles_, event_);
2207  // For thermodynamic lattice output
2209  if (printout_rho_eckart_) {
2210  switch (dens_type_lattice_printout_) {
2211  case DensityType::Baryon:
2212  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2214  break;
2216  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2218  break;
2219  case DensityType::None:
2220  break;
2221  default:
2222  output->at_eventstart(event_, ThermodynamicQuantity::EckartDensity,
2224  *jmu_custom_lat_);
2225  }
2226  }
2227  if (printout_tmn_) {
2228  output->at_eventstart(event_, ThermodynamicQuantity::Tmn,
2230  }
2231  if (printout_tmn_landau_) {
2232  output->at_eventstart(event_, ThermodynamicQuantity::TmnLandau,
2234  }
2235  if (printout_v_landau_) {
2236  output->at_eventstart(event_, ThermodynamicQuantity::LandauVelocity,
2238  }
2239  if (printout_j_QBS_) {
2240  output->at_eventstart(event_, ThermodynamicQuantity::j_QBS,
2242  }
2243  }
2244  }
2245 
2246  /* In the ColliderModus, if Fermi motion is frozen, assign the beam momenta
2247  * to the nucleons in both the projectile and the target. Every ensemble
2248  * gets the same beam momenta, so no need to create beam_momenta_ vector
2249  * for every ensemble.
2250  */
2251  if (modus_.is_collider() && modus_.fermi_motion() == FermiMotion::Frozen) {
2252  for (ParticleData &particle : ensembles_[0]) {
2253  const double m = particle.effective_mass();
2254  double v_beam = 0.0;
2255  if (particle.belongs_to() == BelongsTo::Projectile) {
2256  v_beam = modus_.velocity_projectile();
2257  } else if (particle.belongs_to() == BelongsTo::Target) {
2258  v_beam = modus_.velocity_target();
2259  }
2260  const double gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
2261  beam_momentum_.emplace_back(
2262  FourVector(gamma * m, 0.0, 0.0, gamma * v_beam * m));
2263  } // loop over particles
2264  }
2265 }
double initial_mean_field_energy_
The initial total mean field energy in the system.
Definition: experiment.h:646
bool printout_tmn_
Whether to print the energy-momentum tensor.
Definition: experiment.h:511
QuantumNumbers conserved_initial_
The conserved quantities of the system.
Definition: experiment.h:640
DensityParameters density_param_
Structure to precalculate and hold parameters for density computations.
Definition: experiment.h:374
double total_energy_removed_
Total energy removed from the system in hypersurface crossing actions.
Definition: experiment.h:699
DensityType dens_type_lattice_printout_
Type of density for lattice printout.
Definition: experiment.h:459
double max_transverse_distance_sqr_
Maximal distance at which particles can interact in case of the geometric criterion,...
Definition: experiment.h:631
bool printout_j_QBS_
Whether to print the Q, B, S 4-currents.
Definition: experiment.h:520
const TimeStepMode time_step_mode_
This indicates whether to use time steps.
Definition: experiment.h:625
std::unique_ptr< DensityLattice > jmu_custom_lat_
Custom density on the lattices.
Definition: experiment.h:456
bool printout_full_lattice_any_td_
Whether to print the thermodynamics quantities evaluated on the lattices, point by point,...
Definition: experiment.h:527
std::unique_ptr< DensityLattice > jmu_B_lat_
Baryon density on the lattice.
Definition: experiment.h:438
std::vector< FourVector > beam_momentum_
The initial nucleons in the ColliderModus propagate with beam_momentum_, if Fermi motion is frozen.
Definition: experiment.h:420
std::unique_ptr< RectangularLattice< FourVector > > new_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0 + dt.
Definition: experiment.h:494
const double delta_time_startup_
The clock's timestep size at start up.
Definition: experiment.h:592
std::unique_ptr< RectangularLattice< EnergyMomentumTensor > > Tmn_
Lattices of energy-momentum tensors for printout.
Definition: experiment.h:489
std::vector< Particles > ensembles_
Complete particle list, all ensembles in one vector.
Definition: experiment.h:383
SystemTimePoint time_start_
system starting time of the simulation
Definition: experiment.h:649
uint64_t previous_wall_actions_total_
Total number of wall-crossings for previous timestep.
Definition: experiment.h:676
OutputsList outputs_
A list of output formaters.
Definition: experiment.h:401
bool printout_rho_eckart_
Whether to print the Eckart rest frame density.
Definition: experiment.h:508
bool printout_v_landau_
Whether to print the 4-velocity in Landau frame.
Definition: experiment.h:517
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > EM_lat_
Lattices for electric and magnetic field in fm^-2.
Definition: experiment.h:486
Modus modus_
Instance of the Modus template parameter.
Definition: experiment.h:380
std::unique_ptr< RectangularLattice< FourVector > > old_jmu_auxiliary_
Auxiliary lattice for values of jmu at a time step t0.
Definition: experiment.h:492
std::unique_ptr< DensityLattice > j_QBS_lat_
4-current for j_QBS lattice output
Definition: experiment.h:435
bool printout_tmn_landau_
Whether to print the energy-momentum tensor in Landau frame.
Definition: experiment.h:514
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of jmu.
Definition: experiment.h:497
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:536
ExperimentParameters parameters_
Struct of several member variables.
Definition: experiment.h:371
std::unique_ptr< Potentials > potentials_
An instance of potentials class, that stores parameters of potentials, calculates them and their grad...
Definition: experiment.h:389
uint64_t wall_actions_total_
Total number of wall-crossings for current timestep.
Definition: experiment.h:670
uint64_t total_hypersurface_crossing_actions_
Total number of particles removed from the evolution in hypersurface crossing actions.
Definition: experiment.h:688
uint64_t interactions_total_
Total number of interactions for current timestep.
Definition: experiment.h:658
std::unique_ptr< DensityLattice > jmu_I3_lat_
Isospin projection density on the lattice.
Definition: experiment.h:441
uint64_t total_pauli_blocked_
Total number of Pauli-blockings for current timestep.
Definition: experiment.h:682
std::vector< bool > projectile_target_interact_
Whether the projectile and the target collided.
Definition: experiment.h:413
int64_t seed_
random seed for the next event.
Definition: experiment.h:710
uint64_t previous_interactions_total_
Total number of interactions for previous timestep.
Definition: experiment.h:664
uint64_t discarded_interactions_total_
Total number of discarded interactions, because they were invalidated before they could be performed.
Definition: experiment.h:694
bool kinematic_cuts_for_IC_output_
This indicates whether kinematic cuts are enabled for the IC output.
Definition: experiment.h:707
double total_energy_violated_by_Pythia_
Total energy violation introduced by Pythia.
Definition: experiment.h:704
void init_pythia_hadron_rndm()
Set PYTHIA random seeds to be desired values.
@ Frozen
Use fermi motion without potentials.
@ Fixed
Use fixed time step.
@ None
Don't use time steps; propagate from action to action.
@ EckartDensity
Density in the Eckart frame.
@ Tmn
Energy-momentum tensor in lab frame.
@ LandauVelocity
Velocity of the Landau rest frame.
@ j_QBS
Electric (Q), baryonic (B) and strange (S) currents.
@ TmnLandau
Energy-momentum tensor in Landau rest frame.
Engine::result_type advance()
Advance the engine's state and return the generated value.
Definition: random.h:78
void set_seed(T &&seed)
Sets the seed of the random number engine.
Definition: random.h:71
EventInfo fill_event_info(const std::vector< Particles > &ensembles, double E_mean_field, double modus_impact_parameter, const ExperimentParameters &parameters, bool projectile_target_interact, bool kinematic_cut_for_SMASH_IC)
Generate the EventInfo object which is passed to outputs_.
Definition: experiment.cc:595
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:327
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:375
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:186
const std::string hline(113, '-')
String representing a horizontal line.
int n_ensembles
Number of parallel ensembles.
std::unique_ptr< Clock > outputclock
Output clock to keep track of the next output time.
std::unique_ptr< Clock > labclock
System clock (for simulation time keeping in the computational frame)

◆ run_time_evolution()

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

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

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

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

Definition at line 2439 of file experiment.h.

2441  {
2442  if (!add_plist.empty() || !remove_plist.empty()) {
2443  if (ensembles_.size() > 1) {
2444  throw std::runtime_error(
2445  "Adding or removing particles from SMASH is only possible when one "
2446  "ensemble is used.");
2447  }
2448  const double action_time = parameters_.labclock->current_time();
2449  /* Use two if statements. The first one is to check if the particles are
2450  * valid. Since this might remove all particles, a second if statement is
2451  * needed to avoid executing the action in that case.*/
2452  if (!add_plist.empty()) {
2454  }
2455  if (!add_plist.empty()) {
2456  // Create and perform action to add particle(s)
2457  auto action_add_particles = std::make_unique<FreeforallAction>(
2458  ParticleList{}, add_plist, action_time);
2459  perform_action(*action_add_particles, 0);
2460  }
2461  // Also here 2 if statements are needed as above.
2462  if (!remove_plist.empty()) {
2463  validate_and_adjust_particle_list(remove_plist);
2464  }
2465  if (!remove_plist.empty()) {
2466  ParticleList found_particles_to_remove;
2467  for (const auto &particle_to_remove : remove_plist) {
2468  const auto iterator_to_particle_to_be_removed_in_ensemble =
2469  std::find_if(
2470  ensembles_[0].begin(), ensembles_[0].end(),
2471  [&particle_to_remove, &action_time](const ParticleData &p) {
2473  particle_to_remove, p, action_time);
2474  });
2475  if (iterator_to_particle_to_be_removed_in_ensemble !=
2476  ensembles_[0].end())
2477  found_particles_to_remove.push_back(
2478  *iterator_to_particle_to_be_removed_in_ensemble);
2479  }
2480  // Sort the particles found to be removed according to their id and look
2481  // for duplicates (sorting is needed to call std::adjacent_find).
2482  std::sort(found_particles_to_remove.begin(),
2483  found_particles_to_remove.end(),
2484  [](const ParticleData &p1, const ParticleData &p2) {
2485  return p1.id() < p2.id();
2486  });
2487  const auto iterator_to_first_duplicate = std::adjacent_find(
2488  found_particles_to_remove.begin(), found_particles_to_remove.end(),
2489  [](const ParticleData &p1, const ParticleData &p2) {
2490  return p1.id() == p2.id();
2491  });
2492  if (iterator_to_first_duplicate != found_particles_to_remove.end()) {
2493  logg[LExperiment].error() << "The same particle has been asked to be "
2494  "removed multiple times:\n"
2495  << *iterator_to_first_duplicate;
2496  throw std::logic_error("Particle cannot be removed twice!");
2497  }
2498  if (auto delta = remove_plist.size() - found_particles_to_remove.size();
2499  delta > 0) {
2500  logg[LExperiment].warn(
2501  "When trying to remove particle(s) at the beginning ",
2502  "of the system evolution,\n", delta,
2503  " particle(s) could not be found and will be ignored.");
2504  }
2505  if (!found_particles_to_remove.empty()) {
2506  [[maybe_unused]] const auto number_particles_before_removal =
2507  ensembles_[0].size();
2508  // Create and perform action to remove particles
2509  auto action_remove_particles = std::make_unique<FreeforallAction>(
2510  found_particles_to_remove, ParticleList{}, action_time);
2511  perform_action(*action_remove_particles, 0);
2512 
2513  assert(number_particles_before_removal -
2514  found_particles_to_remove.size() ==
2515  ensembles_[0].size());
2516  }
2517  }
2518  }
2519 
2520  if (t_end > end_time_) {
2521  logg[LExperiment].fatal()
2522  << "Evolution asked to be run until " << t_end << " > " << end_time_
2523  << " and this cannot be done (because of how the clock works).";
2524  throw std::logic_error(
2525  "Experiment cannot evolve the system beyond End_Time.");
2526  }
2527  while (*(parameters_.labclock) < t_end) {
2528  const double dt = parameters_.labclock->timestep_duration();
2529  logg[LExperiment].debug("Timestepless propagation for next ", dt, " fm.");
2530 
2531  // Perform forced thermalization if required
2532  if (thermalizer_ &&
2533  thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
2534  const bool ignore_cells_under_treshold = true;
2535  // Thermodynamics in thermalizer is computed from all ensembles,
2536  // but thermalization actions act on each ensemble independently
2537  thermalizer_->update_thermalizer_lattice(ensembles_, density_param_,
2538  ignore_cells_under_treshold);
2539  const double current_t = parameters_.labclock->current_time();
2540  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2541  thermalizer_->thermalize(ensembles_[i_ens], current_t,
2543  ThermalizationAction th_act(*thermalizer_, current_t);
2544  if (th_act.any_particles_thermalized()) {
2545  perform_action(th_act, i_ens);
2546  }
2547  }
2548  }
2549 
2550  if (IC_dynamic_) {
2551  modus_.build_fluidization_lattice(parameters_.labclock->current_time(),
2553  }
2554 
2555  std::vector<Actions> actions(parameters_.n_ensembles);
2556  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2557  actions[i_ens].clear();
2558  if (ensembles_[i_ens].size() > 0 && action_finders_.size() > 0) {
2559  /* (1.a) Create grid. */
2560  const double min_cell_length = compute_min_cell_length(dt);
2561  logg[LExperiment].debug("Creating grid with minimal cell length ",
2562  min_cell_length);
2563  /* For the hyper-surface-crossing actions also unformed particles are
2564  * searched and therefore needed on the grid. */
2565  const bool include_unformed_particles = IC_switch_;
2566  const auto &grid =
2567  use_grid_ ? modus_.create_grid(ensembles_[i_ens], min_cell_length,
2568  dt, parameters_.coll_crit,
2569  include_unformed_particles)
2570  : modus_.create_grid(ensembles_[i_ens], min_cell_length,
2571  dt, parameters_.coll_crit,
2572  include_unformed_particles,
2574 
2575  const double gcell_vol = grid.cell_volume();
2576  /* (1.b) Iterate over cells and find actions. */
2577  grid.iterate_cells(
2578  [&](const ParticleList &search_list) {
2579  for (const auto &finder : action_finders_) {
2580  actions[i_ens].insert(finder->find_actions_in_cell(
2581  search_list, dt, gcell_vol, beam_momentum_));
2582  }
2583  },
2584  [&](const ParticleList &search_list,
2585  const ParticleList &neighbors_list) {
2586  for (const auto &finder : action_finders_) {
2587  actions[i_ens].insert(finder->find_actions_with_neighbors(
2588  search_list, neighbors_list, dt, beam_momentum_));
2589  }
2590  });
2591  }
2592  }
2593 
2594  /* \todo (optimizations) Adapt timestep size here */
2595 
2596  /* (2) Propagate from action to action until next output or timestep end */
2597  const double end_timestep_time = parameters_.labclock->next_time();
2598  while (next_output_time() < end_timestep_time) {
2599  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2600  run_time_evolution_timestepless(actions[i_ens], i_ens,
2601  next_output_time());
2602  }
2603  ++(*parameters_.outputclock);
2604 
2606  }
2607  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2608  run_time_evolution_timestepless(actions[i_ens], i_ens, end_timestep_time);
2609  }
2610 
2611  /* (3) Update potentials (if computed on the lattice) and
2612  * compute new momenta according to equations of motion */
2613  if (potentials_) {
2615  update_momenta(ensembles_, parameters_.labclock->timestep_duration(),
2616  *potentials_, FB_lat_.get(), FI3_lat_.get(), EM_lat_.get(),
2617  jmu_B_lat_.get());
2618  }
2619 
2620  /* (4) Expand universe if non-minkowskian metric; updates
2621  * positions and momenta according to the selected expansion */
2623  for (Particles &particles : ensembles_) {
2624  expand_space_time(&particles, parameters_, metric_);
2625  }
2626  }
2627 
2628  ++(*parameters_.labclock);
2629 
2630  /* (5) Check conservation laws.
2631  *
2632  * Check conservation of conserved quantities if potentials and string
2633  * fragmentation are off. If potentials are on then momentum is conserved
2634  * only in average. If string fragmentation is on, then energy and
2635  * momentum are only very roughly conserved in high-energy collisions. */
2638  std::string err_msg = conserved_initial_.report_deviations(ensembles_);
2639  if (!err_msg.empty()) {
2640  logg[LExperiment].error() << err_msg;
2641  throw std::runtime_error("Violation of conserved quantities!");
2642  }
2643  }
2644  }
2645 
2646  /* Increment once more the output clock in order to have it prepared for the
2647  * final_output() call. Once close to the end time, the while-loop above to
2648  * produce intermediate output is not entered as the next_output_time() is
2649  * never strictly smaller than end_timestep_time (they are usually equal).
2650  * Since in the final_output() function the current time of the output clock
2651  * is used to produce the output, this has to be incremented before producing
2652  * the final output and it makes sense to do it here.
2653  */
2654  ++(*parameters_.outputclock);
2655 
2656  if (pauli_blocker_) {
2657  logg[LExperiment].info(
2658  "Interactions: Pauli-blocked/performed = ", total_pauli_blocked_, "/",
2660  }
2661 }
const ExpansionProperties metric_
This struct contains information on the metric to be used.
Definition: experiment.h:604
std::vector< std::unique_ptr< ActionFinderInterface > > action_finders_
The Action finder objects.
Definition: experiment.h:423
double next_output_time() const
Shortcut for next output time.
Definition: experiment.h:349
const bool IC_dynamic_
This indicates if the IC is dynamic.
Definition: experiment.h:622
std::unique_ptr< GrandCanThermalizer > thermalizer_
Instance of class used for forced thermalization.
Definition: experiment.h:530
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:2690
void intermediate_output()
Intermediate output during an event.
Definition: experiment.h:2754
std::unique_ptr< RectangularLattice< std::pair< ThreeVector, ThreeVector > > > FI3_lat_
Lattices for the electric and magnetic component of the symmetry force.
Definition: experiment.h:482
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:341
const bool IC_switch_
This indicates whether the experiment will be used as initial condition for hydrodynamics.
Definition: experiment.h:619
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:478
std::unique_ptr< PauliBlocker > pauli_blocker_
An instance of PauliBlocker class that stores parameters needed for Pauli blocking calculations and c...
Definition: experiment.h:395
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:601
void update_potentials()
Recompute potentials on lattices if necessary.
Definition: experiment.h:2905
std::string report_deviations(const std::vector< Particles > &ensembles) const
Checks if the current particle list has still the same values and reports about differences.
constexpr int p
Proton.
void update_momenta(std::vector< Particles > &particles, double dt, const Potentials &pot, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *FB_lat, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *FI3_lat, RectangularLattice< std::pair< ThreeVector, ThreeVector >> *EM_lat, DensityLattice *jB_lat)
Updates the momenta of all particles at the current time step according to the equations of motion:
Definition: propagation.cc:111
void expand_space_time(Particles *particles, const ExperimentParameters &parameters, const ExpansionProperties &metric)
Modifies positions and momentum of all particles to account for space-time deformation.
Definition: propagation.cc:86
bool are_particles_identical_at_given_time(const ParticleData &p1, const ParticleData &p2, double time)
Utility function to compare two ParticleData instances with respect to their PDG code,...
void validate_and_adjust_particle_list(ParticleList &particle_list)
Validate a particle list adjusting each particle to be a valid SMASH particle.
Definition: experiment.cc:617
@ Largest
Make cells as large as possible.
ExpansionMode mode_
Type of metric used.
Definition: propagation.h:28
const CollisionCriterion coll_crit
Employed collision criterion.
bool strings_switch
This indicates whether string fragmentation is switched on.
int testparticles
Number of test-particles.

◆ do_final_decays()

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

Performs the final decays of an event.

Exceptions
runtime_errorif found actions cannot be performed

Definition at line 3006 of file experiment.h.

3006  {
3007  /* At end of time evolution: Force all resonances to decay. In order to handle
3008  * decay chains, we need to loop until no further actions occur. */
3009  bool actions_performed, decays_found;
3010  uint64_t interactions_old;
3011  do {
3012  decays_found = false;
3013  interactions_old = interactions_total_;
3014  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
3015  Actions actions;
3016 
3017  // Dileptons: shining of remaining resonances
3018  if (dilepton_finder_ != nullptr) {
3019  for (const auto &output : outputs_) {
3020  dilepton_finder_->shine_final(ensembles_[i_ens], output.get(), true);
3021  }
3022  }
3023  // Find actions.
3024  for (const auto &finder : action_finders_) {
3025  auto found_actions = finder->find_final_actions(ensembles_[i_ens]);
3026  if (!found_actions.empty()) {
3027  actions.insert(std::move(found_actions));
3028  decays_found = true;
3029  }
3030  }
3031  // Perform actions.
3032  while (!actions.is_empty()) {
3033  perform_action(*actions.pop(), i_ens, false);
3034  }
3035  }
3036  actions_performed = interactions_total_ > interactions_old;
3037  // Throw an error if actions were found but not performed
3038  if (decays_found && !actions_performed) {
3039  throw std::runtime_error("Final decays were found but not performed.");
3040  }
3041  // loop until no more decays occur
3042  } while (actions_performed);
3043 
3044  // Dileptons: shining of stable particles at the end
3045  if (dilepton_finder_ != nullptr) {
3046  for (const auto &output : outputs_) {
3047  for (Particles &particles : ensembles_) {
3048  dilepton_finder_->shine_final(particles, output.get(), false);
3049  }
3050  }
3051  }
3052 }
std::unique_ptr< DecayActionsFinderDilepton > dilepton_finder_
The Dilepton Action Finder.
Definition: experiment.h:426

◆ final_output()

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

Output at the end of an event.

Definition at line 3055 of file experiment.h.

3055  {
3056  /* make sure the experiment actually ran (note: we should compare this
3057  * to the start time, but we don't know that. Therefore, we check that
3058  * the time is positive, which should heuristically be the same). */
3059  double E_mean_field = 0.0;
3060  if (likely(parameters_.labclock > 0)) {
3061  const uint64_t wall_actions_this_interval =
3063  const uint64_t interactions_this_interval = interactions_total_ -
3065  wall_actions_this_interval;
3066  if (potentials_) {
3067  // using the lattice is necessary
3068  if ((jmu_B_lat_ != nullptr)) {
3070  EM_lat_.get(), parameters_);
3071  }
3072  }
3073  if (std::abs(parameters_.labclock->current_time() - end_time_) >
3074  really_small) {
3075  logg[LExperiment].warn()
3076  << "SMASH not propagated until configured end time. Current time = "
3077  << parameters_.labclock->current_time()
3078  << "fm. End time = " << end_time_ << "fm.";
3079  } else {
3081  ensembles_, interactions_this_interval, conserved_initial_,
3083  }
3084  int total_particles = 0;
3085  for (const Particles &particles : ensembles_) {
3086  total_particles += particles.size();
3087  }
3088  if (IC_switch_ && (total_particles == 0)) {
3089  const double initial_system_energy_plus_Pythia_violations =
3091  const double fraction_of_total_system_energy_removed =
3092  initial_system_energy_plus_Pythia_violations / total_energy_removed_;
3093  // Verify there is no more energy in the system if all particles were
3094  // removed when crossing the hypersurface
3095  if (std::fabs(fraction_of_total_system_energy_removed - 1.) >
3096  really_small) {
3097  throw std::runtime_error(
3098  "There is remaining energy in the system although all particles "
3099  "were removed.\n"
3100  "E_remain = " +
3101  std::to_string((initial_system_energy_plus_Pythia_violations -
3103  " [GeV]");
3104  } else {
3105  logg[LExperiment].info() << hline;
3106  logg[LExperiment].info()
3107  << "Time real: " << SystemClock::now() - time_start_;
3108  logg[LExperiment].info()
3109  << "Interactions before reaching hypersurface: "
3112  logg[LExperiment].info()
3113  << "Total number of particles removed on hypersurface: "
3115  }
3116  } else {
3117  const double precent_discarded =
3119  ? static_cast<double>(discarded_interactions_total_) * 100.0 /
3121  : 0.0;
3122  std::stringstream msg_discarded;
3123  msg_discarded
3124  << "Discarded interaction number: " << discarded_interactions_total_
3125  << " (" << precent_discarded
3126  << "% of the total interaction number including wall crossings)";
3127 
3128  logg[LExperiment].info() << hline;
3129  logg[LExperiment].info()
3130  << "Time real: " << SystemClock::now() - time_start_;
3131  logg[LExperiment].debug() << msg_discarded.str();
3132 
3134  precent_discarded > 1.0) {
3135  // The choosen threshold of 1% is a heuristical value
3136  logg[LExperiment].warn()
3137  << msg_discarded.str()
3138  << "\nThe number of discarded interactions is large, which means "
3139  "the assumption for the stochastic criterion of\n1 interaction "
3140  "per particle per timestep is probably violated. Consider "
3141  "reducing the timestep size.";
3142  }
3143 
3144  logg[LExperiment].info() << "Final interaction number: "
3146  }
3147 
3148  // Check if there are unformed particles
3149  int unformed_particles_count = 0;
3150  for (const Particles &particles : ensembles_) {
3151  for (const ParticleData &particle : particles) {
3152  if (particle.formation_time() > end_time_) {
3153  unformed_particles_count++;
3154  }
3155  }
3156  }
3157  if (unformed_particles_count > 0) {
3158  logg[LExperiment].warn(
3159  "End time might be too small. ", unformed_particles_count,
3160  " unformed particles were found at the end of the evolution.");
3161  }
3162  }
3163 
3164  // Keep track of how many ensembles had interactions
3166 
3167  for (const auto &output : outputs_) {
3168  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
3169  auto event_info = fill_event_info(
3170  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
3172  output->at_eventend(ensembles_[i_ens], {event_, i_ens}, event_info);
3173  }
3174  // For thermodynamic output
3175  output->at_eventend(ensembles_, event_);
3176 
3177  // For thermodynamic lattice output
3178  if (printout_rho_eckart_) {
3180  output->at_eventend(ThermodynamicQuantity::EckartDensity);
3181  }
3182  }
3183  if (printout_tmn_) {
3184  output->at_eventend(ThermodynamicQuantity::Tmn);
3185  }
3186  if (printout_tmn_landau_) {
3187  output->at_eventend(ThermodynamicQuantity::TmnLandau);
3188  }
3189  if (printout_v_landau_) {
3190  output->at_eventend(ThermodynamicQuantity::LandauVelocity);
3191  }
3192  if (printout_j_QBS_) {
3193  output->at_eventend(ThermodynamicQuantity::j_QBS);
3194  }
3195  }
3196 }
void count_nonempty_ensembles()
Counts the number of ensembles in wich interactions took place at the end of an event.
Definition: experiment.h:3199
double x0() const
Definition: fourvector.h:313
FourVector momentum() const
@ Stochastic
Stochastic Criteiron.
#define likely(x)
Tell the branch predictor that this expression is likely true.
Definition: macros.h:14
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37

◆ first_ensemble()

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

Provides external access to SMASH particles.

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

Definition at line 260 of file experiment.h.

260 { return &ensembles_[0]; }

◆ all_ensembles()

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

Getter for all ensembles.

Definition at line 262 of file experiment.h.

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

268 { 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 3234 of file experiment.h.

3234  {
3235  event_++;
3236 }

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

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

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

2665  {
2666  const double dt =
2667  propagate_straight_line(&particles, to_time, beam_momentum_);
2668  if (dilepton_finder_ != nullptr) {
2669  for (const auto &output : outputs_) {
2670  dilepton_finder_->shine(particles, output.get(), dt);
2671  }
2672  }
2673 }
double propagate_straight_line(Particles *particles, double to_time, const std::vector< FourVector > &beam_momentum)
Propagates the positions of all particles on a straight line to a given moment.
Definition: propagation.cc:44

◆ run_time_evolution_timestepless()

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

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

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

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

Definition at line 2690 of file experiment.h.

2691  {
2692  Particles &particles = ensembles_[i_ensemble];
2693  logg[LExperiment].debug(
2694  "Timestepless propagation: ", "Actions size = ", actions.size(),
2695  ", end time = ", end_time_propagation);
2696 
2697  // iterate over all actions
2698  while (!actions.is_empty()) {
2699  if (actions.earliest_time() > end_time_propagation) {
2700  break;
2701  }
2702  // get next action
2703  ActionPtr act = actions.pop();
2704  if (!act->is_valid(particles)) {
2706  logg[LExperiment].debug(~einhard::DRed(), "✘ ", act,
2707  " (discarded: invalid)");
2708  continue;
2709  }
2710  logg[LExperiment].debug(~einhard::Green(), "✔ ", act,
2711  ", action time = ", act->time_of_execution());
2712 
2713  /* (1) Propagate to the next action. */
2714  propagate_and_shine(act->time_of_execution(), particles);
2715 
2716  /* (2) Perform action.
2717  *
2718  * Update the positions of the incoming particles, because the information
2719  * in the action object will be outdated as the particles have been
2720  * propagated since the construction of the action. */
2721  act->update_incoming(particles);
2722  const bool performed = perform_action(*act, i_ensemble);
2723 
2724  /* No need to update actions for outgoing particles
2725  * if the action is not performed. */
2726  if (!performed) {
2727  continue;
2728  }
2729 
2730  /* (3) Update actions for newly-produced particles. */
2731 
2732  const double end_time_timestep = parameters_.labclock->next_time();
2733  // New actions are always search until the end of the current timestep
2734  const double time_left = end_time_timestep - act->time_of_execution();
2735  const ParticleList &outgoing_particles = act->outgoing_particles();
2736  // Grid cell volume set to zero, since there is no grid
2737  const double gcell_vol = 0.0;
2738  for (const auto &finder : action_finders_) {
2739  // Outgoing particles can still decay, cross walls...
2740  actions.insert(finder->find_actions_in_cell(outgoing_particles, time_left,
2741  gcell_vol, beam_momentum_));
2742  // ... and collide with other particles.
2743  actions.insert(finder->find_actions_with_surrounding_particles(
2744  outgoing_particles, particles, time_left, beam_momentum_));
2745  }
2746 
2748  }
2749 
2750  propagate_and_shine(end_time_propagation, particles);
2751 }
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:2664
void check_interactions_total(uint64_t interactions_total)
Make sure interactions_total can be represented as a 32-bit integer.
Definition: experiment.h:2682

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

2754  {
2755  const uint64_t wall_actions_this_interval =
2758  const uint64_t interactions_this_interval = interactions_total_ -
2760  wall_actions_this_interval;
2762  double E_mean_field = 0.0;
2765  double computational_frame_time = 0.0;
2766  if (potentials_) {
2767  // using the lattice is necessary
2768  if ((jmu_B_lat_ != nullptr)) {
2770  EM_lat_.get(), parameters_);
2771  /*
2772  * Mean field calculated in a box should remain approximately constant if
2773  * the system is in equilibrium, and so deviations from its original value
2774  * may signal a phase transition or other dynamical process. This
2775  * comparison only makes sense in the Box Modus, hence the condition.
2776  */
2777  if (modus_.is_box()) {
2778  double tmp = (E_mean_field - initial_mean_field_energy_) /
2779  (E_mean_field + initial_mean_field_energy_);
2780  /*
2781  * This is displayed when the system evolves away from its initial
2782  * configuration (which is when the total mean field energy in the box
2783  * deviates from its initial value).
2784  */
2785  if (std::abs(tmp) > 0.01) {
2786  logg[LExperiment].info()
2787  << "\n\n\n\t The mean field at t = "
2788  << parameters_.outputclock->current_time()
2789  << " [fm] differs from the mean field at t = 0:"
2790  << "\n\t\t initial_mean_field_energy_ = "
2791  << initial_mean_field_energy_ << " [GeV]"
2792  << "\n\t\t abs[(E_MF - E_MF(t=0))/(E_MF + E_MF(t=0))] = "
2793  << std::abs(tmp)
2794  << "\n\t\t E_MF/E_MF(t=0) = "
2795  << E_mean_field / initial_mean_field_energy_ << "\n\n";
2796  }
2797  }
2798  }
2799  }
2800 
2802  ensembles_, interactions_this_interval, conserved_initial_, time_start_,
2803  parameters_.outputclock->current_time(), E_mean_field,
2805  const LatticeUpdate lat_upd = LatticeUpdate::AtOutput;
2806 
2807  // save evolution data
2808  if (!(modus_.is_box() && parameters_.outputclock->current_time() <
2809  modus_.equilibration_time())) {
2810  for (const auto &output : outputs_) {
2811  if (output->is_dilepton_output() || output->is_photon_output() ||
2812  output->is_IC_output()) {
2813  continue;
2814  }
2815  for (int i_ens = 0; i_ens < parameters_.n_ensembles; i_ens++) {
2816  auto event_info = fill_event_info(
2817  ensembles_, E_mean_field, modus_.impact_parameter(), parameters_,
2819 
2820  output->at_intermediate_time(ensembles_[i_ens], parameters_.outputclock,
2821  density_param_, {event_, i_ens},
2822  event_info);
2823  computational_frame_time = event_info.current_time;
2824  }
2825  // For thermodynamic output
2826  output->at_intermediate_time(ensembles_, parameters_.outputclock,
2827  density_param_);
2828 
2829  // Thermodynamic output on the lattice versus time
2830  if (printout_rho_eckart_) {
2831  switch (dens_type_lattice_printout_) {
2832  case DensityType::Baryon:
2835  ensembles_, false);
2836  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2838  output->thermodynamics_lattice_output(*jmu_B_lat_,
2839  computational_frame_time);
2840  break;
2843  jmu_I3_lat_.get(), lat_upd, DensityType::BaryonicIsospin,
2844  density_param_, ensembles_, false);
2845  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2847  *jmu_I3_lat_);
2848  output->thermodynamics_lattice_output(*jmu_I3_lat_,
2849  computational_frame_time);
2850  break;
2851  case DensityType::None:
2852  break;
2853  default:
2856  density_param_, ensembles_, false);
2857  output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
2859  *jmu_custom_lat_);
2860  output->thermodynamics_lattice_output(*jmu_custom_lat_,
2861  computational_frame_time);
2862  }
2863  }
2867  ensembles_, false);
2868  if (printout_tmn_) {
2869  output->thermodynamics_output(ThermodynamicQuantity::Tmn,
2871  output->thermodynamics_lattice_output(
2872  ThermodynamicQuantity::Tmn, *Tmn_, computational_frame_time);
2873  }
2874  if (printout_tmn_landau_) {
2875  output->thermodynamics_output(ThermodynamicQuantity::TmnLandau,
2877  output->thermodynamics_lattice_output(
2879  computational_frame_time);
2880  }
2881  if (printout_v_landau_) {
2882  output->thermodynamics_output(ThermodynamicQuantity::LandauVelocity,
2884  output->thermodynamics_lattice_output(
2886  computational_frame_time);
2887  }
2888  }
2889  if (EM_lat_) {
2890  output->fields_output("Efield", "Bfield", *EM_lat_);
2891  }
2892  if (printout_j_QBS_) {
2893  output->thermodynamics_lattice_output(
2894  *j_QBS_lat_, computational_frame_time, ensembles_, density_param_);
2895  }
2896 
2897  if (thermalizer_) {
2898  output->thermodynamics_output(*thermalizer_);
2899  }
2900  }
2901  }
2902 }
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:644
LatticeUpdate
Enumerator option for lattice updates.
Definition: lattice.h:38

◆ update_potentials()

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

Recompute potentials on lattices if necessary.

Definition at line 2905 of file experiment.h.

2905  {
2906  if (potentials_) {
2907  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
2912  parameters_.labclock->timestep_duration(), true);
2913  }
2914  if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
2915  jmu_B_lat_ != nullptr) {
2920  parameters_.labclock->timestep_duration(), true);
2921  const size_t UBlattice_size = UB_lat_->size();
2922  for (size_t i = 0; i < UBlattice_size; i++) {
2923  auto jB = (*jmu_B_lat_)[i];
2924  const FourVector flow_four_velocity_B =
2925  std::abs(jB.rho()) > very_small_double ? jB.jmu_net() / jB.rho()
2926  : FourVector();
2927  double baryon_density = jB.rho();
2928  ThreeVector baryon_grad_j0 = jB.grad_j0();
2929  ThreeVector baryon_dvecj_dt = jB.dvecj_dt();
2930  ThreeVector baryon_curl_vecj = jB.curl_vecj();
2931  if (potentials_->use_skyrme()) {
2932  (*UB_lat_)[i] =
2933  flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
2934  (*FB_lat_)[i] =
2935  potentials_->skyrme_force(baryon_density, baryon_grad_j0,
2936  baryon_dvecj_dt, baryon_curl_vecj);
2937  }
2938  if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
2939  auto jI3 = (*jmu_I3_lat_)[i];
2940  const FourVector flow_four_velocity_I3 =
2941  std::abs(jI3.rho()) > very_small_double
2942  ? jI3.jmu_net() / jI3.rho()
2943  : FourVector();
2944  (*UI3_lat_)[i] = flow_four_velocity_I3 *
2945  potentials_->symmetry_pot(jI3.rho(), baryon_density);
2946  (*FI3_lat_)[i] = potentials_->symmetry_force(
2947  jI3.rho(), jI3.grad_j0(), jI3.dvecj_dt(), jI3.curl_vecj(),
2948  baryon_density, baryon_grad_j0, baryon_dvecj_dt,
2949  baryon_curl_vecj);
2950  }
2951  }
2952  }
2953  if (potentials_->use_coulomb()) {
2956  density_param_, ensembles_, true);
2957  for (size_t i = 0; i < EM_lat_->size(); i++) {
2958  ThreeVector electric_field = {0., 0., 0.};
2959  ThreeVector position = jmu_el_lat_->cell_center(i);
2960  jmu_el_lat_->integrate_volume(electric_field,
2962  potentials_->coulomb_r_cut(), position);
2963  ThreeVector magnetic_field = {0., 0., 0.};
2964  jmu_el_lat_->integrate_volume(magnetic_field,
2966  potentials_->coulomb_r_cut(), position);
2967  (*EM_lat_)[i] = std::make_pair(electric_field, magnetic_field);
2968  }
2969  } // if ((potentials_->use_skyrme() || ...
2970  if (potentials_->use_vdf() && jmu_B_lat_ != nullptr) {
2975  parameters_.labclock->timestep_duration(), true);
2978  fields_lat_.get(), old_fields_auxiliary_.get(),
2981  parameters_.labclock->timestep_duration());
2982  }
2983  const size_t UBlattice_size = UB_lat_->size();
2984  for (size_t i = 0; i < UBlattice_size; i++) {
2985  auto jB = (*jmu_B_lat_)[i];
2986  (*UB_lat_)[i] = potentials_->vdf_pot(jB.rho(), jB.jmu_net());
2989  (*FB_lat_)[i] = potentials_->vdf_force(
2990  jB.rho(), jB.drho_dxnu().x0(), jB.drho_dxnu().threevec(),
2991  jB.grad_rho_cross_vecj(), jB.jmu_net().x0(), jB.grad_j0(),
2992  jB.jmu_net().threevec(), jB.dvecj_dt(), jB.curl_vecj());
2993  break;
2995  auto Amu = (*fields_lat_)[i];
2996  (*FB_lat_)[i] = potentials_->vdf_force(
2997  Amu.grad_A0(), Amu.dvecA_dt(), Amu.curl_vecA());
2998  break;
2999  }
3000  } // for (size_t i = 0; i < UBlattice_size; i++)
3001  } // if potentials_->use_vdf()
3002  }
3003 }
std::unique_ptr< RectangularLattice< FourVector > > new_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0 + dt.
Definition: experiment.h:502
std::unique_ptr< RectangularLattice< std::array< FourVector, 4 > > > fields_four_gradient_auxiliary_
Auxiliary lattice for calculating the four-gradient of Amu.
Definition: experiment.h:505
std::unique_ptr< FieldsLattice > fields_lat_
Mean-field A^mu on the lattice.
Definition: experiment.h:447
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:465
std::unique_ptr< RectangularLattice< FourVector > > old_fields_auxiliary_
Auxiliary lattice for values of Amu at a time step t0.
Definition: experiment.h:500
std::unique_ptr< DensityLattice > jmu_el_lat_
Electric charge density on the lattice.
Definition: experiment.h:444
static ThreeVector B_field_integrand(ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
Integrand for calculating the magnetic field using the Biot-Savart formula.
Definition: potentials.h:338
static ThreeVector E_field_integrand(ThreeVector pos, DensityOnLattice &charge_density, ThreeVector point)
Integrand for calculating the electric field.
Definition: potentials.h:320
constexpr double very_small_double
A very small double, used to avoid division by zero.
Definition: constants.h:40
void update_fields_lattice(RectangularLattice< FieldsOnLattice > *fields_lat, RectangularLattice< FourVector > *old_fields, RectangularLattice< FourVector > *new_fields, RectangularLattice< std::array< FourVector, 4 >> *fields_four_grad_lattice, DensityLattice *jmu_B_lat, const LatticeUpdate fields_lat_update, const Potentials &potentials, const double time_step)
Updates the contents on the lattice of FieldsOnLattice type.
Definition: fields.cc:14
FieldDerivativesMode field_derivatives_mode
mode of calculating field derivatives

◆ compute_min_cell_length()

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

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

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

Definition at line 341 of file experiment.h.

341  {
344  }
345  return std::sqrt(4 * dt * dt + max_transverse_distance_sqr_);
346  }
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 349 of file experiment.h.

349  {
350  return parameters_.outputclock->next_time();
351  }

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

3199  {
3200  for (bool has_interaction : projectile_target_interact_) {
3201  if (has_interaction) {
3203  }
3204  }
3205 }
int nonempty_ensembles_
Number of ensembles containing an interaction.
Definition: experiment.h:576

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

3208  {
3210  return event_ >= nevents_;
3211  }
3214  return true;
3215  }
3216  if (event_ >= max_events_) {
3217  logg[LExperiment].warn()
3218  << "Maximum number of events (" << max_events_
3219  << ") exceeded. Stopping calculation. "
3220  << "The fraction of empty ensembles is "
3221  << (1.0 - static_cast<double>(nonempty_ensembles_) /
3223  << ". If this fraction is expected, try increasing the "
3224  "Maximum_Ensembles_Run.";
3225  return true;
3226  }
3227  return false;
3228  }
3229  throw std::runtime_error("Event counting option is invalid");
3230  return false;
3231 }
int minimum_nonempty_ensembles_
The number of ensembles, in which interactions take place, to be calculated.
Definition: experiment.h:562
int max_events_
Maximum number of events to be calculated in order obtain the desired number of non-empty events usin...
Definition: experiment.h:582
EventCounting event_counting_
The way in which the number of calculated events is specified.
Definition: experiment.h:570
int nevents_
Number of events.
Definition: experiment.h:552
@ 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 192 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 371 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 374 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 380 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 383 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 389 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 395 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 401 of file experiment.h.

◆ dilepton_output_

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

The Dilepton output.

Definition at line 404 of file experiment.h.

◆ photon_output_

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

The Photon output.

Definition at line 407 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 413 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 420 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 423 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 426 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 429 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 432 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 435 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 438 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 441 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 444 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 447 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 456 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 459 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 465 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 471 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 478 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 482 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 486 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 489 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 492 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 494 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 497 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 500 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 502 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 505 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 508 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 511 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 514 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 517 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 520 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 523 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 527 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 530 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 536 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 552 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 562 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 570 of file experiment.h.

◆ event_

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

Current event.

Definition at line 573 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 576 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 582 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 585 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 592 of file experiment.h.

◆ force_decays_

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

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

Definition at line 598 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 601 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 604 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 607 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 610 of file experiment.h.

◆ bremsstrahlung_switch_

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

This indicates whether bremsstrahlung is switched on.

Definition at line 613 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 619 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 622 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 625 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 631 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 640 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 646 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 649 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 652 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 658 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 664 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 670 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 676 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 682 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 688 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 694 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 699 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 704 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 707 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 710 of file experiment.h.


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