experiment.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2013-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#ifndef SRC_INCLUDE_EXPERIMENT_H_
#define SRC_INCLUDE_EXPERIMENT_H_

#include <limits>
#include <memory>
#include <string>
#include <vector>

#include "actionfinderfactory.h"
#include "actions.h"
#include "chrono.h"
#include "decayactionsfinder.h"
#include "decayactionsfinderdilepton.h"
#include "energymomentumtensor.h"
#include "fourvector.h"
#include "grandcan_thermalizer.h"
#include "grid.h"
#include "hypersurfacecrossingaction.h"
#include "outputparameters.h"
#include "pauliblocking.h"
#include "potential_globals.h"
#include "potentials.h"
#include "processstring.h"
#include "propagation.h"
#include "quantumnumbers.h"
#include "scatteractionphoton.h"
#include "scatteractionsfinder.h"
#include "thermalizationaction.h"
// Output
#include "binaryoutput.h"
#include "icoutput.h"
#include "oscaroutput.h"
#include "thermodynamicoutput.h"
#ifdef SMASH_USE_ROOT
#include "rootoutput.h"
#endif
#include "vtkoutput.h"
#include "wallcrossingaction.h"

namespace std {
template <typename T, typename Ratio>
static ostream &operator<<(ostream &out,
                           const chrono::duration<T, Ratio> &seconds) {
  using Seconds = chrono::duration<double>;
  using Minutes = chrono::duration<double, std::ratio<60>>;
  using Hours = chrono::duration<double, std::ratio<60 * 60>>;
  constexpr Minutes threshold_for_minutes{10};
  constexpr Hours threshold_for_hours{3};
  if (seconds < threshold_for_minutes) {
    return out << Seconds(seconds).count() << " [s]";
  }
  if (seconds < threshold_for_hours) {
    return out << Minutes(seconds).count() << " [min]";
  }
  return out << Hours(seconds).count() << " [h]";
}
}  // namespace std

namespace smash {

class ExperimentBase {
 public:
  ExperimentBase() = default;
  virtual ~ExperimentBase() = default;

  static std::unique_ptr<ExperimentBase> create(Configuration config,
                                                const bf::path &output_path);

  virtual void run() = 0;

  struct InvalidModusRequest : public std::invalid_argument {
    using std::invalid_argument::invalid_argument;
  };
};

template <typename Modus>
class Experiment;
template <typename Modus>
std::ostream &operator<<(std::ostream &out, const Experiment<Modus> &e);

template <typename Modus>
class Experiment : public ExperimentBase {
  friend class ExperimentBase;

 public:
  void run() override;

  explicit Experiment(Configuration config, const bf::path &output_path);

  void initialize_new_event();

  void run_time_evolution();

  void do_final_decays();

  void final_output(const int evt_num);

  Particles *particles() { return &particles_; }

  Modus *modus() { return &modus_; }

 private:
  template <typename Container>
  bool perform_action(Action &action,
                      const Container &particles_before_actions);
  void create_output(const std::string &format, const std::string &content,
                     const bf::path &output_path, const OutputParameters &par);

  void propagate_and_shine(double to_time);

  void run_time_evolution_timestepless(Actions &actions);

  void intermediate_output();

  void update_potentials();

  double compute_min_cell_length(double dt) const {
    return std::sqrt(4 * dt * dt + max_transverse_distance_sqr_);
  }

  double next_output_time() const {
    return parameters_.outputclock->next_time();
  }

  ExperimentParameters parameters_;

  DensityParameters density_param_;

  Modus modus_;

  Particles particles_;

  std::unique_ptr<Potentials> potentials_;

  std::unique_ptr<PauliBlocker> pauli_blocker_;

  OutputsList outputs_;

  OutputPtr dilepton_output_;

  OutputPtr photon_output_;

  std::vector<bool> nucleon_has_interacted_ = {};
  bool projectile_target_interact_ = false;

  std::vector<FourVector> beam_momentum_ = {};

  std::vector<std::unique_ptr<ActionFinderInterface>> action_finders_;

  std::unique_ptr<DecayActionsFinderDilepton> dilepton_finder_;

  std::unique_ptr<ActionFinderInterface> photon_finder_;

  int n_fractional_photons_ = 100;

  std::unique_ptr<DensityLattice> jmu_B_lat_;

  std::unique_ptr<DensityLattice> jmu_I3_lat_;

  std::unique_ptr<DensityLattice> jmu_custom_lat_;

  DensityType dens_type_lattice_printout_ = DensityType::None;

  std::unique_ptr<RectangularLattice<FourVector>> UB_lat_ = nullptr;

  std::unique_ptr<RectangularLattice<FourVector>> UI3_lat_ = nullptr;

  std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector>>>
      FB_lat_;

  std::unique_ptr<RectangularLattice<std::pair<ThreeVector, ThreeVector>>>
      FI3_lat_;

  std::unique_ptr<RectangularLattice<EnergyMomentumTensor>> Tmn_;

  bool printout_tmn_ = false;

  bool printout_tmn_landau_ = false;

  bool printout_v_landau_ = false;

  bool printout_lattice_td_ = false;

  std::unique_ptr<GrandCanThermalizer> thermalizer_;

  StringProcess *process_string_ptr_;

  const int nevents_;

  const double end_time_;

  const double delta_time_startup_;

  const bool force_decays_;

  const bool use_grid_;

  const ExpansionProperties metric_;

  const bool dileptons_switch_;

  const bool photons_switch_;

  const bool IC_output_switch_;

  const TimeStepMode time_step_mode_;

  double max_transverse_distance_sqr_ = std::numeric_limits<double>::max();

  QuantumNumbers conserved_initial_;

  SystemTimePoint time_start_ = SystemClock::now();

  DensityType dens_type_ = DensityType::None;

  uint64_t interactions_total_ = 0;

  uint64_t previous_interactions_total_ = 0;

  uint64_t wall_actions_total_ = 0;

  uint64_t previous_wall_actions_total_ = 0;

  uint64_t total_pauli_blocked_ = 0;

  uint64_t total_hypersurface_crossing_actions_ = 0;

  double total_energy_removed_ = 0.0;

  int64_t seed_ = -1;

  friend std::ostream &operator<<<>(std::ostream &out, const Experiment &e);
};

template <typename Modus>
std::ostream &operator<<(std::ostream &out, const Experiment<Modus> &e) {
  out << "End time: " << e.end_time_ << " fm/c\n";
  out << e.modus_;
  return out;
}

template <typename Modus>
void Experiment<Modus>::create_output(const std::string &format,
                                      const std::string &content,
                                      const bf::path &output_path,
                                      const OutputParameters &out_par) {
  const auto &log = logger<LogArea::Experiment>();
  log.info() << "Adding output " << content << " of format " << format
             << std::endl;

  if (format == "VTK" && content == "Particles") {
    outputs_.emplace_back(
        make_unique<VtkOutput>(output_path, content, out_par));
  } else if (format == "Root") {
#ifdef SMASH_USE_ROOT
    if (content == "Initial_Conditions") {
      outputs_.emplace_back(
          make_unique<RootOutput>(output_path, "SMASH_IC", out_par));
    } else {
      outputs_.emplace_back(
          make_unique<RootOutput>(output_path, content, out_par));
    }
#else
    log.error("Root output requested, but Root support not compiled in");
#endif
  } else if (format == "Binary") {
    if (content == "Collisions" || content == "Dileptons" ||
        content == "Photons") {
      outputs_.emplace_back(
          make_unique<BinaryOutputCollisions>(output_path, content, out_par));
    } else if (content == "Particles") {
      outputs_.emplace_back(
          make_unique<BinaryOutputParticles>(output_path, content, out_par));
    } else if (content == "Initial_Conditions") {
      outputs_.emplace_back(make_unique<BinaryOutputInitialConditions>(
          output_path, content, out_par));
    }
  } else if (format == "Oscar1999" || format == "Oscar2013") {
    outputs_.emplace_back(
        create_oscar_output(format, content, output_path, out_par));
  } else if (content == "Thermodynamics" && format == "ASCII") {
    outputs_.emplace_back(
        make_unique<ThermodynamicOutput>(output_path, content, out_par));
  } else if (content == "Thermodynamics" && format == "VTK") {
    printout_lattice_td_ = true;
    outputs_.emplace_back(
        make_unique<VtkOutput>(output_path, content, out_par));
  } else if (content == "Initial_Conditions" && format == "ASCII") {
    outputs_.emplace_back(
        make_unique<ICOutput>(output_path, "SMASH_IC", out_par));
  } else {
    log.error() << "Unknown combination of format (" << format
                << ") and content (" << content << "). Fix the config.";
  }
}

ExperimentParameters create_experiment_parameters(Configuration config);

template <typename Modus>
Experiment<Modus>::Experiment(Configuration config, const bf::path &output_path)
    : parameters_(create_experiment_parameters(config)),
      density_param_(DensityParameters(parameters_)),
      modus_(config["Modi"], parameters_),
      particles_(),
      nevents_(config.take({"General", "Nevents"})),
      end_time_(config.take({"General", "End_Time"})),
      delta_time_startup_(parameters_.labclock->timestep_duration()),
      force_decays_(
          config.take({"Collision_Term", "Force_Decays_At_End"}, true)),
      use_grid_(config.take({"General", "Use_Grid"}, true)),
      metric_(
          config.take({"General", "Metric_Type"}, ExpansionMode::NoExpansion),
          config.take({"General", "Expansion_Rate"}, 0.1)),
      dileptons_switch_(config.has_value({"Output", "Dileptons"})),
      photons_switch_(config.has_value({"Output", "Photons"})),
      IC_output_switch_(config.has_value({"Output", "Initial_Conditions"})),
      time_step_mode_(
          config.take({"General", "Time_Step_Mode"}, TimeStepMode::Fixed)) {
  const auto &log = logger<LogArea::Experiment>();
  log.info() << *this;

  // create finders
  if (dileptons_switch_) {
    dilepton_finder_ = make_unique<DecayActionsFinderDilepton>();
  }
  if (parameters_.photons_switch) {
    n_fractional_photons_ = config.take({"Output", "Photons", "Fractions"});
  }
  if (parameters_.two_to_one) {
    action_finders_.emplace_back(make_unique<DecayActionsFinder>());
  }
  bool no_coll = config.take({"Collision_Term", "No_Collisions"}, false);
  if ((parameters_.two_to_one || parameters_.included_2to2.any() ||
       parameters_.strings_switch) &&
      !no_coll) {
    auto scat_finder = make_unique<ScatterActionsFinder>(
        config, parameters_, nucleon_has_interacted_, modus_.total_N_number(),
        modus_.proj_N_number());
    max_transverse_distance_sqr_ =
        scat_finder->max_transverse_distance_sqr(parameters_.testparticles);
    process_string_ptr_ = scat_finder->get_process_string_ptr();
    action_finders_.emplace_back(std::move(scat_finder));
  } else {
    max_transverse_distance_sqr_ = maximum_cross_section / M_PI * fm2_mb;
    process_string_ptr_ = NULL;
  }
  const double modus_l = modus_.length();
  if (modus_l > 0.) {
    action_finders_.emplace_back(make_unique<WallCrossActionsFinder>(modus_l));
  }
  if (IC_output_switch_) {
    if (!modus_.is_collider()) {
      throw std::runtime_error(
          "Initial conditions can only be extracted in collider modus.");
    }
    double proper_time;
    if (config.has_value({"Output", "Initial_Conditions", "Proper_Time"})) {
      // Read in proper time from config
      proper_time =
          config.take({"Output", "Initial_Conditions", "Proper_Time"});
    } else {
      // Default proper time is the passing time of the two nuclei
      proper_time = modus_.nuclei_passing_time();
    }
    action_finders_.emplace_back(
        make_unique<HyperSurfaceCrossActionsFinder>(proper_time));
  }

  if (config.has_value({"Collision_Term", "Pauli_Blocking"})) {
    log.info() << "Pauli blocking is ON.";
    pauli_blocker_ = make_unique<PauliBlocker>(
        config["Collision_Term"]["Pauli_Blocking"], parameters_);
  }
  ParticleData::formation_power_ =
      config.take({"Collision_Term", "Power_Particle_Formation"}, 1.);

  // create outputs
  log.trace(source_location, " create OutputInterface objects");

  auto output_conf = config["Output"];
  dens_type_ = config.take({"Output", "Density_Type"}, DensityType::None);
  log.debug() << "Density type printed to headers: " << dens_type_;

  const OutputParameters output_parameters(std::move(output_conf));

  std::vector<std::string> output_contents = output_conf.list_upmost_nodes();
  for (const auto &content : output_contents) {
    auto this_output_conf = output_conf[content.c_str()];
    const std::vector<std::string> formats = this_output_conf.take({"Format"});
    if (output_path == "") {
      continue;
    }
    for (const auto &format : formats) {
      create_output(format, content, output_path, output_parameters);
    }
  }

  /* We can take away the Fermi motion flag, because the collider modus is
   * already initialized. We only need it when potentials are enabled, but we
   * always have to take it, otherwise SMASH will complain about unused
   * options. We have to provide a default value for modi other than Collider.
   */
  const FermiMotion motion =
      config.take({"Modi", "Collider", "Fermi_Motion"}, FermiMotion::Off);
  if (config.has_value({"Potentials"})) {
    if (time_step_mode_ == TimeStepMode::None) {
      log.error() << "Potentials only work with time steps!";
      throw std::invalid_argument("Can't use potentials without time steps!");
    }
    if (motion == FermiMotion::Frozen) {
      log.error() << "Potentials don't work with frozen Fermi momenta! "
                     "Use normal Fermi motion instead.";
      throw std::invalid_argument(
          "Can't use potentials "
          "with frozen Fermi momenta!");
    }
    log.info() << "Potentials are ON. Timestep is "
               << parameters_.labclock->timestep_duration();
    // potentials need testparticles and gaussian sigma from parameters_
    potentials_ = make_unique<Potentials>(config["Potentials"], parameters_);
  }

  // Create lattices
  if (config.has_value({"Lattice"})) {
    // Take lattice properties from config to assign them to all lattices
    const std::array<double, 3> l = config.take({"Lattice", "Sizes"});
    const std::array<int, 3> n = config.take({"Lattice", "Cell_Number"});
    const std::array<double, 3> origin = config.take({"Lattice", "Origin"});
    const bool periodic = config.take({"Lattice", "Periodic"});

    if (printout_lattice_td_) {
      dens_type_lattice_printout_ = output_parameters.td_dens_type;
      printout_tmn_ = output_parameters.td_tmn;
      printout_tmn_landau_ = output_parameters.td_tmn_landau;
      printout_v_landau_ = output_parameters.td_v_landau;
    }
    if (printout_tmn_ || printout_tmn_landau_ || printout_v_landau_) {
      Tmn_ = make_unique<RectangularLattice<EnergyMomentumTensor>>(
          l, n, origin, periodic, LatticeUpdate::AtOutput);
    }
    /* Create baryon and isospin density lattices regardless of config
       if potentials are on. This is because they allow to compute
       potentials faster */
    if (potentials_) {
      if (potentials_->use_skyrme()) {
        jmu_B_lat_ = make_unique<DensityLattice>(l, n, origin, periodic,
                                                 LatticeUpdate::EveryTimestep);
        UB_lat_ = make_unique<RectangularLattice<FourVector>>(
            l, n, origin, periodic, LatticeUpdate::EveryTimestep);
        FB_lat_ = make_unique<
            RectangularLattice<std::pair<ThreeVector, ThreeVector>>>(
            l, n, origin, periodic, LatticeUpdate::EveryTimestep);
      }
      if (potentials_->use_symmetry()) {
        jmu_I3_lat_ = make_unique<DensityLattice>(l, n, origin, periodic,
                                                  LatticeUpdate::EveryTimestep);
        UI3_lat_ = make_unique<RectangularLattice<FourVector>>(
            l, n, origin, periodic, LatticeUpdate::EveryTimestep);
        FI3_lat_ = make_unique<
            RectangularLattice<std::pair<ThreeVector, ThreeVector>>>(
            l, n, origin, periodic, LatticeUpdate::EveryTimestep);
      }
    } else {
      if (dens_type_lattice_printout_ == DensityType::Baryon) {
        jmu_B_lat_ = make_unique<DensityLattice>(l, n, origin, periodic,
                                                 LatticeUpdate::AtOutput);
      }
      if (dens_type_lattice_printout_ == DensityType::BaryonicIsospin) {
        jmu_I3_lat_ = make_unique<DensityLattice>(l, n, origin, periodic,
                                                  LatticeUpdate::AtOutput);
      }
    }
    if (dens_type_lattice_printout_ != DensityType::None &&
        dens_type_lattice_printout_ != DensityType::BaryonicIsospin &&
        dens_type_lattice_printout_ != DensityType::Baryon) {
      jmu_custom_lat_ = make_unique<DensityLattice>(l, n, origin, periodic,
                                                    LatticeUpdate::AtOutput);
    }
  } else if (printout_lattice_td_) {
    log.error(
        "If you want Thermodynamic VTK output, configure a lattice for it.");
  }

  // Store pointers to potential and lattice accessible for Action
  if (parameters_.potential_affect_threshold) {
    UB_lat_pointer = UB_lat_.get();
    UI3_lat_pointer = UI3_lat_.get();
    pot_pointer = potentials_.get();
  }

  // Create forced thermalizer
  if (config.has_value({"Forced_Thermalization"})) {
    Configuration &&th_conf = config["Forced_Thermalization"];
    thermalizer_ = modus_.create_grandcan_thermalizer(th_conf);
  }

  /* Take the seed setting only after the configuration was stored to a file
   * in smash.cc */
  seed_ = config.take({"General", "Randomseed"});
}

const std::string hline(67, '-');

template <typename Modus>
void Experiment<Modus>::initialize_new_event() {
  const auto &log = logger<LogArea::Experiment>();

  random::set_seed(seed_);
  log.info() << "random number seed: " << seed_;
  /* Set seed for the next event. It has to be positive, so it can be entered
   * in the config.
   *
   * We have to be careful about the minimal integer, whose absolute value
   * cannot be represented. */
  int64_t r = random::advance();
  while (r == INT64_MIN) {
    r = random::advance();
  }
  seed_ = std::abs(r);
  /* Set the random seed used in PYTHIA hadronization
   * to be same with the SMASH one.
   * In this way we ensure that the results are reproducible
   * for every event if one knows SMASH random seed. */
  if (process_string_ptr_ != NULL) {
    process_string_ptr_->init_pythia_hadron_rndm();
  }

  particles_.reset();

  // Sample particles according to the initial conditions
  double start_time = modus_.initial_conditions(&particles_, parameters_);
  /* For box modus make sure that particles are in the box. In principle, after
   * a correct initialization they should be, so this is just playing it safe.
   */
  modus_.impose_boundary_conditions(&particles_, outputs_);
  // Reset the simulation clock
  double timestep = delta_time_startup_;

  switch (time_step_mode_) {
    case TimeStepMode::Fixed:
      break;
    case TimeStepMode::None:
      timestep = end_time_ - start_time;
      // Take care of the box modus + timestepless propagation
      const double max_dt = modus_.max_timestep(max_transverse_distance_sqr_);
      if (max_dt > 0. && max_dt < timestep) {
        timestep = max_dt;
      }
      break;
  }
  std::unique_ptr<UniformClock> clock_for_this_event;
  if (modus_.is_list() && (timestep < 0.0)) {
    throw std::runtime_error(
        "Timestep for the given event is negative. \n"
        "This might happen if the formation times of the input particles are "
        "larger than the specified end time of the simulation.");
  }
  clock_for_this_event = make_unique<UniformClock>(start_time, timestep);
  parameters_.labclock = std::move(clock_for_this_event);

  // Reset the output clock
  parameters_.outputclock->reset(start_time, true);
  // remove time before starting time in case of custom output times.
  parameters_.outputclock->remove_times_in_past(start_time);

  log.debug("Lab clock: t_start = ", parameters_.labclock->current_time(),
            ", dt = ", parameters_.labclock->timestep_duration());

  /* Save the initial conserved quantum numbers and total momentum in
   * the system for conservation checks */
  conserved_initial_ = QuantumNumbers(particles_);
  wall_actions_total_ = 0;
  previous_wall_actions_total_ = 0;
  interactions_total_ = 0;
  previous_interactions_total_ = 0;
  total_pauli_blocked_ = 0;
  projectile_target_interact_ = false;
  total_hypersurface_crossing_actions_ = 0;
  total_energy_removed_ = 0.0;
  // Print output headers
  log.info() << hline;
  log.info() << "Time [fm]   Ediff [GeV]    Scatt.|Decays   "
                "Particles         Timing";
  log.info() << hline;
}

template <typename Modus>
template <typename Container>
bool Experiment<Modus>::perform_action(
    Action &action, const Container &particles_before_actions) {
  const auto &log = logger<LogArea::Experiment>();
  // Make sure to skip invalid and Pauli-blocked actions.
  if (!action.is_valid(particles_)) {
    log.debug(~einhard::DRed(), "✘ ", action, " (discarded: invalid)");
    return false;
  }
  action.generate_final_state();
  log.debug("Process Type is: ", action.get_type());
  if (pauli_blocker_ && action.is_pauli_blocked(particles_, *pauli_blocker_)) {
    total_pauli_blocked_++;
    return false;
  }
  if (modus_.is_collider()) {
    /* Mark incoming nucleons as interacted - now they are permitted
     * to collide with nucleons from their native nucleus */
    bool incoming_projectile = false;
    bool incoming_target = false;
    for (const auto &incoming : action.incoming_particles()) {
      assert(incoming.id() >= 0);
      if (incoming.id() < modus_.total_N_number()) {
        nucleon_has_interacted_[incoming.id()] = true;
      }
      if (incoming.id() < modus_.proj_N_number()) {
        incoming_projectile = true;
      }
      if (incoming.id() >= modus_.proj_N_number() &&
          incoming.id() < modus_.total_N_number()) {
        incoming_target = true;
      }
    }
    // Check whether particles from different nuclei interacted.
    if (incoming_projectile & incoming_target) {
      projectile_target_interact_ = true;
    }
  }
  /* Make sure to pick a non-zero integer, because 0 is reserved for "no
   * interaction yet". */
  const auto id_process = static_cast<uint32_t>(interactions_total_ + 1);
  action.perform(&particles_, id_process);
  interactions_total_++;
  if (action.get_type() == ProcessType::Wall) {
    wall_actions_total_++;
  }
  if (action.get_type() == ProcessType::HyperSurfaceCrossing) {
    total_hypersurface_crossing_actions_++;
    total_energy_removed_ += action.incoming_particles()[0].momentum().x0();
  }
  // Calculate Eckart rest frame density at the interaction point
  double rho = 0.0;
  if (dens_type_ != DensityType::None) {
    const FourVector r_interaction = action.get_interaction_point();
    constexpr bool compute_grad = false;
    const bool smearing = true;
    rho = std::get<0>(current_eckart(r_interaction.threevec(),
                                     particles_before_actions, density_param_,
                                     dens_type_, compute_grad, smearing));
  }
  for (const auto &output : outputs_) {
    if (!output->is_dilepton_output() && !output->is_photon_output()) {
      if (output->is_IC_output() &&
          action.get_type() == ProcessType::HyperSurfaceCrossing) {
        output->at_interaction(action, rho);
      } else if (!output->is_IC_output()) {
        output->at_interaction(action, rho);
      }
    }
  }

  // At every collision photons can be produced.
  // Note: We rely here on the lazy evaluation of the arguments to if.
  // It may happen that in a wall-crossing-action sqrt_s raises an exception.
  // Therefore we first have to check if the incoming particles can undergo
  // an em-interaction.
  if (photons_switch_ &&
      ScatterActionPhoton::is_photon_reaction(action.incoming_particles()) &&
      ScatterActionPhoton::is_kinematically_possible(
          action.sqrt_s(), action.incoming_particles())) {
    /* Time in the action constructor is relative to
     * current time of incoming */
    constexpr double action_time = 0.;
    ScatterActionPhoton photon_act(action.incoming_particles(), action_time,
                                   n_fractional_photons_,
                                   action.get_total_weight());

    photon_act.add_dummy_hadronic_process(action.get_total_weight());

    // Now add the actual photon reaction channel.
    photon_act.add_single_process();

    photon_act.perform_photons(outputs_);
  }
  log.debug(~einhard::Green(), "✔ ", action);
  return true;
}

std::string format_measurements(const Particles &particles,
                                uint64_t scatterings_this_interval,
                                const QuantumNumbers &conserved_initial,
                                SystemTimePoint time_start, double time);

template <typename Modus>
void Experiment<Modus>::run_time_evolution() {
  Actions actions;

  const auto &log = logger<LogArea::Experiment>();

  log.info() << format_measurements(particles_, 0u, conserved_initial_,
                                    time_start_,
                                    parameters_.labclock->current_time());

  while (parameters_.labclock->current_time() < end_time_) {
    const double t = parameters_.labclock->current_time();
    const double dt =
        std::min(parameters_.labclock->timestep_duration(), end_time_ - t);
    log.debug("Timestepless propagation for next ", dt, " fm/c.");

    // Perform forced thermalization if required
    if (thermalizer_ &&
        thermalizer_->is_time_to_thermalize(parameters_.labclock)) {
      const bool ignore_cells_under_treshold = true;
      thermalizer_->update_thermalizer_lattice(particles_, density_param_,
                                               ignore_cells_under_treshold);
      const double current_t = parameters_.labclock->current_time();
      thermalizer_->thermalize(particles_, current_t,
                               parameters_.testparticles);
      ThermalizationAction th_act(*thermalizer_, current_t);
      if (th_act.any_particles_thermalized()) {
        perform_action(th_act, particles_);
      }
    }

    if (particles_.size() > 0 && action_finders_.size() > 0) {
      /* (1.a) Create grid. */
      double min_cell_length = compute_min_cell_length(dt);
      log.debug("Creating grid with minimal cell length ", min_cell_length);
      const auto &grid =
          use_grid_ ? modus_.create_grid(particles_, min_cell_length, dt)
                    : modus_.create_grid(particles_, min_cell_length, dt,
                                         CellSizeStrategy::Largest);

      const double cell_vol = grid.cell_volume();

      /* (1.b) Iterate over cells and find actions. */
      grid.iterate_cells(
          [&](const ParticleList &search_list) {
            for (const auto &finder : action_finders_) {
              actions.insert(finder->find_actions_in_cell(
                  search_list, dt, cell_vol, beam_momentum_));
            }
          },
          [&](const ParticleList &search_list,
              const ParticleList &neighbors_list) {
            for (const auto &finder : action_finders_) {
              actions.insert(finder->find_actions_with_neighbors(
                  search_list, neighbors_list, dt, beam_momentum_));
            }
          });
    }

    /* \todo (optimizations) Adapt timestep size here */

    /* (2) Propagation from action to action until the end of timestep */
    run_time_evolution_timestepless(actions);

    /* (3) Update potentials (if computed on the lattice) and
     *     compute new momenta according to equations of motion */
    if (potentials_) {
      update_potentials();
      update_momenta(&particles_, parameters_.labclock->timestep_duration(),
                     *potentials_, FB_lat_.get(), FI3_lat_.get());
    }

    /* (4) Expand universe if non-minkowskian metric; updates
     *     positions and momenta according to the selected expansion */
    if (metric_.mode_ != ExpansionMode::NoExpansion) {
      expand_space_time(&particles_, parameters_, metric_);
    }

    ++(*parameters_.labclock);

    /* (5) Check conservation laws.
     *
     * Check conservation of conserved quantities if potentials and string
     * fragmentation are off.  If potentials are on then momentum is conserved
     * only in average.  If string fragmentation is on, then energy and
     * momentum are only very roughly conserved in high-energy collisions. */
    if (!potentials_ && !parameters_.strings_switch &&
        metric_.mode_ == ExpansionMode::NoExpansion && !IC_output_switch_) {
      std::string err_msg = conserved_initial_.report_deviations(particles_);
      if (!err_msg.empty()) {
        log.error() << err_msg;
        throw std::runtime_error("Violation of conserved quantities!");
      }
    }
  }

  if (pauli_blocker_) {
    log.info("Interactions: Pauli-blocked/performed = ", total_pauli_blocked_,
             "/", interactions_total_ - wall_actions_total_);
  }
}

template <typename Modus>
void Experiment<Modus>::propagate_and_shine(double to_time) {
  const double dt =
      propagate_straight_line(&particles_, to_time, beam_momentum_);
  if (dilepton_finder_ != nullptr) {
    for (const auto &output : outputs_) {
      dilepton_finder_->shine(particles_, output.get(), dt);
    }
  }
}

inline void check_interactions_total(uint64_t interactions_total) {
  constexpr uint64_t max_uint32 = std::numeric_limits<uint32_t>::max();
  if (interactions_total >= max_uint32) {
    throw std::runtime_error("Integer overflow in total interaction number!");
  }
}

template <typename Modus>
void Experiment<Modus>::run_time_evolution_timestepless(Actions &actions) {
  const auto &log = logger<LogArea::Experiment>();

  const double start_time = parameters_.labclock->current_time();
  const double end_time =
      std::min(parameters_.labclock->next_time(), end_time_);
  double time_left = end_time - start_time;
  log.debug("Timestepless propagation: ", "Actions size = ", actions.size(),
            ", start time = ", start_time, ", end time = ", end_time);

  // iterate over all actions
  while (!actions.is_empty()) {
    // get next action
    ActionPtr act = actions.pop();
    if (!act->is_valid(particles_)) {
      log.debug(~einhard::DRed(), "✘ ", act, " (discarded: invalid)");
      continue;
    }
    if (act->time_of_execution() > end_time) {
      log.error(act, " scheduled later than end time: t_action[fm/c] = ",
                act->time_of_execution(), ", t_end[fm/c] = ", end_time);
    }
    log.debug(~einhard::Green(), "✔ ", act);

    while (next_output_time() <= act->time_of_execution()) {
      log.debug("Propagating until output time: ", next_output_time());
      propagate_and_shine(next_output_time());
      ++(*parameters_.outputclock);
      intermediate_output();
    }

    /* (1) Propagate to the next action. */
    log.debug("Propagating until next action ", act,
              ", action time = ", act->time_of_execution());
    propagate_and_shine(act->time_of_execution());

    /* (2) Perform action.
     *
     * Update the positions of the incoming particles, because the information
     * in the action object will be outdated as the particles have been
     * propagated since the construction of the action. */
    act->update_incoming(particles_);
    const bool performed = perform_action(*act, particles_);

    /* No need to update actions for outgoing particles
     * if the action is not performed. */
    if (!performed) {
      continue;
    }

    /* (3) Update actions for newly-produced particles. */

    time_left = end_time - act->time_of_execution();
    const ParticleList &outgoing_particles = act->outgoing_particles();
    // Cell volume set to zero, since there is no grid
    const double cell_vol = 0.0;
    for (const auto &finder : action_finders_) {
      // Outgoing particles can still decay, cross walls...
      actions.insert(finder->find_actions_in_cell(outgoing_particles, time_left,
                                                  cell_vol, beam_momentum_));
      // ... and collide with other particles.
      actions.insert(finder->find_actions_with_surrounding_particles(
          outgoing_particles, particles_, time_left, beam_momentum_));
    }

    check_interactions_total(interactions_total_);
  }

  while (next_output_time() <= end_time) {
    log.debug("Propagating until output time: ", next_output_time());
    propagate_and_shine(next_output_time());
    ++(*parameters_.outputclock);
    // Avoid duplicating printout at event end time
    if (parameters_.outputclock->current_time() < end_time_) {
      intermediate_output();
    }
  }
  log.debug("Propagating to time ", end_time);
  propagate_and_shine(end_time);
}

template <typename Modus>
void Experiment<Modus>::intermediate_output() {
  const auto &log = logger<LogArea::Experiment>();
  const uint64_t wall_actions_this_interval =
      wall_actions_total_ - previous_wall_actions_total_;
  previous_wall_actions_total_ = wall_actions_total_;
  const uint64_t interactions_this_interval = interactions_total_ -
                                              previous_interactions_total_ -
                                              wall_actions_this_interval;
  previous_interactions_total_ = interactions_total_;
  log.info() << format_measurements(particles_, interactions_this_interval,
                                    conserved_initial_, time_start_,
                                    parameters_.outputclock->current_time());
  const LatticeUpdate lat_upd = LatticeUpdate::AtOutput;
  // save evolution data
  for (const auto &output : outputs_) {
    if (output->is_dilepton_output() || output->is_photon_output() ||
        output->is_IC_output()) {
      continue;
    }

    output->at_intermediate_time(particles_, parameters_.outputclock,
                                 density_param_);

    // Thermodynamic output on the lattice versus time
    switch (dens_type_lattice_printout_) {
      case DensityType::Baryon:
        update_lattice(jmu_B_lat_.get(), lat_upd, DensityType::Baryon,
                       density_param_, particles_, false);
        output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                      DensityType::Baryon, *jmu_B_lat_);
        break;
      case DensityType::BaryonicIsospin:
        update_lattice(jmu_I3_lat_.get(), lat_upd, DensityType::BaryonicIsospin,
                       density_param_, particles_, false);
        output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                      DensityType::BaryonicIsospin,
                                      *jmu_I3_lat_);
        break;
      case DensityType::None:
        break;
      default:
        update_lattice(jmu_custom_lat_.get(), lat_upd,
                       dens_type_lattice_printout_, density_param_, particles_,
                       false);
        output->thermodynamics_output(ThermodynamicQuantity::EckartDensity,
                                      dens_type_lattice_printout_,
                                      *jmu_custom_lat_);
    }
    if (printout_tmn_ || printout_tmn_landau_ || printout_v_landau_) {
      update_lattice(Tmn_.get(), lat_upd, dens_type_lattice_printout_,
                     density_param_, particles_);
      if (printout_tmn_) {
        output->thermodynamics_output(ThermodynamicQuantity::Tmn,
                                      dens_type_lattice_printout_, *Tmn_);
      }
      if (printout_tmn_landau_) {
        output->thermodynamics_output(ThermodynamicQuantity::TmnLandau,
                                      dens_type_lattice_printout_, *Tmn_);
      }
      if (printout_v_landau_) {
        output->thermodynamics_output(ThermodynamicQuantity::LandauVelocity,
                                      dens_type_lattice_printout_, *Tmn_);
      }
    }

    if (thermalizer_) {
      output->thermodynamics_output(*thermalizer_);
    }
  }
}

template <typename Modus>
void Experiment<Modus>::update_potentials() {
  if (potentials_) {
    if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
      update_lattice(jmu_I3_lat_.get(), LatticeUpdate::EveryTimestep,
                     DensityType::BaryonicIsospin, density_param_, particles_,
                     true);
    }
    if ((potentials_->use_skyrme() || potentials_->use_symmetry()) &&
        jmu_B_lat_ != nullptr) {
      update_lattice(jmu_B_lat_.get(), LatticeUpdate::EveryTimestep,
                     DensityType::Baryon, density_param_, particles_, true);
      const size_t UBlattice_size = UB_lat_->size();
      assert(UBlattice_size == UI3_lat_->size());
      for (size_t i = 0; i < UBlattice_size; i++) {
        auto jB = (*jmu_B_lat_)[i];
        const FourVector flow_four_velocity_B =
            std::abs(jB.density()) > really_small ? jB.jmu_net() / jB.density()
                                                  : FourVector();
        double baryon_density = jB.density();
        ThreeVector baryon_grad_rho = jB.grad_rho();
        ThreeVector baryon_dj_dt = jB.dj_dt();
        ThreeVector baryon_rot_j = jB.rot_j();
        if (potentials_->use_skyrme()) {
          (*UB_lat_)[i] =
              flow_four_velocity_B * potentials_->skyrme_pot(baryon_density);
          (*FB_lat_)[i] = potentials_->skyrme_force(
              baryon_density, baryon_grad_rho, baryon_dj_dt, baryon_rot_j);
        }
        if (potentials_->use_symmetry() && jmu_I3_lat_ != nullptr) {
          auto jI3 = (*jmu_I3_lat_)[i];
          const FourVector flow_four_velocity_I3 =
              std::abs(jI3.density()) > really_small
                  ? jI3.jmu_net() / jI3.density()
                  : FourVector();
          (*UI3_lat_)[i] =
              flow_four_velocity_I3 *
              potentials_->symmetry_pot(jI3.density(), baryon_density);
          (*FI3_lat_)[i] = potentials_->symmetry_force(
              jI3.density(), jI3.grad_rho(), jI3.dj_dt(), jI3.rot_j(),
              baryon_density, baryon_grad_rho, baryon_dj_dt, baryon_rot_j);
        }
      }
    }
  }
}

template <typename Modus>
void Experiment<Modus>::do_final_decays() {
  /* At end of time evolution: Force all resonances to decay. In order to handle
   * decay chains, we need to loop until no further actions occur. */
  uint64_t interactions_old;
  const auto particles_before_actions = particles_.copy_to_vector();
  do {
    Actions actions;

    interactions_old = interactions_total_;

    // Dileptons: shining of remaining resonances
    if (dilepton_finder_ != nullptr) {
      for (const auto &output : outputs_) {
        dilepton_finder_->shine_final(particles_, output.get(), true);
      }
    }
    // Find actions.
    for (const auto &finder : action_finders_) {
      actions.insert(finder->find_final_actions(particles_));
    }
    // Perform actions.
    while (!actions.is_empty()) {
      perform_action(*actions.pop(), particles_before_actions);
    }
    // loop until no more decays occur
  } while (interactions_total_ > interactions_old);

  // Dileptons: shining of stable particles at the end
  if (dilepton_finder_ != nullptr) {
    for (const auto &output : outputs_) {
      dilepton_finder_->shine_final(particles_, output.get(), false);
    }
  }
}

template <typename Modus>
void Experiment<Modus>::final_output(const int evt_num) {
  const auto &log = logger<LogArea::Experiment>();
  /* make sure the experiment actually ran (note: we should compare this
   * to the start time, but we don't know that. Therefore, we check that
   * the time is positive, which should heuristically be the same). */
  if (likely(parameters_.labclock > 0)) {
    const uint64_t wall_actions_this_interval =
        wall_actions_total_ - previous_wall_actions_total_;
    const uint64_t interactions_this_interval = interactions_total_ -
                                                previous_interactions_total_ -
                                                wall_actions_this_interval;
    log.info() << format_measurements(particles_, interactions_this_interval,
                                      conserved_initial_, time_start_,
                                      end_time_);
    if (IC_output_switch_ && (particles_.size() == 0)) {
      // Verify there is no more energy in the system if all particles were
      // removed when crossing the hypersurface
      const double remaining_energy =
          conserved_initial_.momentum().x0() - total_energy_removed_;
      if (remaining_energy > really_small) {
        throw std::runtime_error(
            "There is remaining energy in the system although all particles "
            "were removed.\n"
            "E_remain = " +
            std::to_string(remaining_energy) + " [GeV]");
      } else {
        log.info() << hline;
        log.info() << "Time real: " << SystemClock::now() - time_start_;
        log.info() << "Interactions before reaching hypersurface: "
                   << interactions_total_ - wall_actions_total_ -
                          total_hypersurface_crossing_actions_;
        log.info() << "Total number of particles removed on hypersurface: "
                   << total_hypersurface_crossing_actions_;
      }
    } else {
      log.info() << hline;
      log.info() << "Time real: " << SystemClock::now() - time_start_;
      log.info() << "Final interaction number: "
                 << interactions_total_ - wall_actions_total_;
    }

    // Check if there are unformed particles
    int unformed_particles_count = 0;
    for (const auto &particle : particles_) {
      if (particle.formation_time() > end_time_) {
        unformed_particles_count++;
      }
    }
    if (unformed_particles_count > 0) {
      log.warn("End time might be too small. ", unformed_particles_count,
               " unformed particles were found at the end of the evolution.");
    }
  }

  for (const auto &output : outputs_) {
    output->at_eventend(particles_, evt_num, modus_.impact_parameter(),
                        !projectile_target_interact_);
  }
}

template <typename Modus>
void Experiment<Modus>::run() {
  const auto &mainlog = logger<LogArea::Main>();
  for (int j = 0; j < nevents_; j++) {
    mainlog.info() << "Event " << j;

    // Sample initial particles, start clock, some printout and book-keeping
    initialize_new_event();
    /* In the ColliderModus, if the first collisions within the same nucleus are
     * forbidden, 'nucleon_has_interacted_', which records whether a nucleon has
     * collided with another nucleon, is initialized equal to false. If allowed,
     * 'nucleon_has_interacted' is initialized equal to true, which means these
     * incoming particles have experienced some fake scatterings, they can
     * therefore collide with each other later on since these collisions are not
     * "first" to them. */
    if (modus_.is_collider()) {
      if (!modus_.cll_in_nucleus()) {
        nucleon_has_interacted_.assign(modus_.total_N_number(), false);
      } else {
        nucleon_has_interacted_.assign(modus_.total_N_number(), true);
      }
    }
    /* In the ColliderModus, if Fermi motion is frozen, assign the beam momenta
     * to the nucleons in both the projectile and the target. */
    if (modus_.is_collider() && modus_.fermi_motion() == FermiMotion::Frozen) {
      for (int i = 0; i < modus_.total_N_number(); i++) {
        const auto mass_beam = particles_.copy_to_vector()[i].effective_mass();
        const auto v_beam = i < modus_.proj_N_number()
                                ? modus_.velocity_projectile()
                                : modus_.velocity_target();
        const auto gamma = 1.0 / std::sqrt(1.0 - v_beam * v_beam);
        beam_momentum_.emplace_back(FourVector(gamma * mass_beam, 0.0, 0.0,
                                               gamma * v_beam * mass_beam));
      }
    }

    // Output at event start
    for (const auto &output : outputs_) {
      output->at_eventstart(particles_, j);
    }

    run_time_evolution();

    if (force_decays_) {
      do_final_decays();
    }

    // Output at event end
    final_output(j);
  }
}

}  // namespace smash

#endif  // SRC_INCLUDE_EXPERIMENT_H_