experiment.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2012-2020
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/experiment.h"
 
#include <cstdint>
 
#include "smash/boxmodus.h"
#include "smash/collidermodus.h"
#include "smash/cxx14compat.h"
#include "smash/listmodus.h"
#include "smash/spheremodus.h"
 
namespace smash {
 
/* ExperimentBase carries everything that is needed for the evolution */
ExperimentPtr ExperimentBase::create(Configuration config,
                                     const bf::path &output_path) {
  logg[LExperiment].trace() << source_location;
  const std::string modus_chooser = config.read({"General", "Modus"});
  logg[LExperiment].debug() << "Modus for this calculation: " << modus_chooser;
 
  if (modus_chooser == "Box") {
    return make_unique<Experiment<BoxModus>>(config, output_path);
  } else if (modus_chooser == "List") {
    return make_unique<Experiment<ListModus>>(config, output_path);
  } else if (modus_chooser == "Collider") {
    return make_unique<Experiment<ColliderModus>>(config, output_path);
  } else if (modus_chooser == "Sphere") {
    return make_unique<Experiment<SphereModus>>(config, output_path);
  } else {
    throw InvalidModusRequest("Invalid Modus (" + modus_chooser +
                              ") requested from ExperimentBase::create.");
  }
}
 
ExperimentParameters create_experiment_parameters(Configuration config) {
  logg[LExperiment].trace() << source_location;
 
  const int ntest = config.take({"General", "Testparticles"}, 1);
  if (ntest <= 0) {
    throw std::invalid_argument("Testparticle number should be positive!");
  }
 
  const std::string modus_chooser = config.take({"General", "Modus"});
  // remove config maps of unused Modi
  config["Modi"].remove_all_but(modus_chooser);
 
  double box_length = -1.0;
  if (config.has_value({"Modi", "Box", "Length"})) {
    box_length = config.read({"Modi", "Box", "Length"});
  }
 
  /* If this Delta_Time option is absent (this can be for timestepless mode)
   * just assign 1.0 fm/c, reasonable value will be set at event initialization
   */
  const double dt = config.take({"General", "Delta_Time"}, 1.);
  const double t_end = config.read({"General", "End_Time"});
 
  // Enforce a small time step, if the box modus is used
  if (box_length > 0.0 && dt > box_length / 10.0) {
    throw std::invalid_argument(
        "Please decrease the timestep size. "
        "A value of (dt < l_box / 10) is recommended in the boxmodus.");
  }
 
  // define output clock
  std::unique_ptr<Clock> output_clock = nullptr;
  if (config.has_value({"Output", "Output_Times"})) {
    if (config.has_value({"Output", "Output_Interval"})) {
      throw std::invalid_argument(
          "Please specify either Output_Interval or Output_Times");
    }
    std::vector<double> output_times = config.take({"Output", "Output_Times"});
    // Add an output time larger than the end time so that the next time is
    // always defined during the time evolution
    output_times.push_back(t_end + 1.);
    output_clock = make_unique<CustomClock>(output_times);
  } else {
    const double output_dt = config.take({"Output", "Output_Interval"}, t_end);
    output_clock = make_unique<UniformClock>(0.0, output_dt);
  }
 
  // Add proper error messages if photons are not configured properly.
  // 1) Missing Photon config section.
  if (config["Output"].has_value({"Photons"}) &&
      (!config.has_value({"Collision_Term", "Photons"}))) {
    throw std::invalid_argument(
        "Photon output is enabled although photon production is disabled. "
        "Photon production can be configured in the \"Photon\" subsection "
        "of the \"Collision_Term\".");
  }
 
  // 2) Missing Photon output section.
  bool missing_output_2to2 = false;
  bool missing_output_brems = false;
  if (!(config["Output"].has_value({"Photons"}))) {
    if (config.has_value({"Collision_Term", "Photons", "2to2_Scatterings"})) {
      missing_output_2to2 =
          config.read({"Collision_Term", "Photons", "2to2_Scatterings"});
    }
    if (config.has_value({"Collision_Term", "Photons", "Bremsstrahlung"})) {
      missing_output_brems =
          config.read({"Collision_Term", "Photons", "Bremsstrahlung"});
    }
 
    if (missing_output_2to2 || missing_output_brems) {
      throw std::invalid_argument(
          "Photon output is disabled although photon production is enabled. "
          "Please enable the photon output.");
    }
  }
 
  // Add proper error messages if dileptons are not configured properly.
  // 1) Missing Dilepton config section.
  if (config["Output"].has_value({"Dileptons"}) &&
      (!config.has_value({"Collision_Term", "Dileptons"}))) {
    throw std::invalid_argument(
        "Dilepton output is enabled although dilepton production is disabled. "
        "Dilepton production can be configured in the \"Dileptons\" subsection "
        "of the \"Collision_Term\".");
  }
 
  // 2) Missing Dilepton output section.
  bool missing_output_decays = false;
  if (!(config["Output"].has_value({"Dileptons"}))) {
    if (config.has_value({"Collision_Term", "Dileptons", "Decays"})) {
      missing_output_decays =
          config.read({"Collision_Term", "Dileptons", "Decays"});
    }
 
    if (missing_output_decays) {
      throw std::invalid_argument(
          "Dilepton output is disabled although dilepton production is "
          "enabled. "
          "Please enable the dilepton output.");
    }
  }
 
  auto config_coll = config["Collision_Term"];
  /* Elastic collisions between the nucleons with the square root s
   * below low_snn_cut are excluded. */
  const double low_snn_cut =
      config_coll.take({"Elastic_NN_Cutoff_Sqrts"}, 1.98);
  const auto proton = ParticleType::try_find(pdg::p);
  const auto pion = ParticleType::try_find(pdg::pi_z);
  if (proton && pion &&
      low_snn_cut > proton->mass() + proton->mass() + pion->mass()) {
    logg[LExperiment].warn("The cut-off should be below the threshold energy",
                           " of the process: NN to NNpi");
  }
  const bool potential_affect_threshold =
      config.take({"Lattice", "Potentials_Affect_Thresholds"}, false);
  const double scale_xs = config_coll.take({"Cross_Section_Scaling"}, 1.0);
  const double maximum_cross_section_default =
      ParticleType::exists("d'") ? 2000.0 : 200.0;
  double maximum_cross_section = config_coll.take(
      {"Maximum_Cross_Section"}, maximum_cross_section_default);
  maximum_cross_section *= scale_xs;
  return {
      make_unique<UniformClock>(0.0, dt),
      std::move(output_clock),
      ntest,
      config.take({"General", "Gaussian_Sigma"}, 1.),
      config.take({"General", "Gauss_Cutoff_In_Sigma"}, 4.),
      config_coll.take({"Collision_Criterion"}, CollisionCriterion::Covariant),
      config_coll.take({"Two_to_One"}, true),
      config_coll.take({"Included_2to2"}, ReactionsBitSet().set()),
      config_coll.take({"Multi_Particle_Reactions"},
                       MultiParticleReactionsBitSet().reset()),
      config_coll.take({"Strings"}, modus_chooser != "Box"),
      config_coll.take({"Use_AQM"}, true),
      config_coll.take({"Resonance_Lifetime_Modifier"}, 1.),
      config_coll.take({"Strings_with_Probability"}, true),
      config_coll.take({"NNbar_Treatment"}, NNbarTreatment::Strings),
      low_snn_cut,
      potential_affect_threshold,
      box_length,
      maximum_cross_section,
      scale_xs,
      config_coll.take({"Additional_Elastic_Cross_Section"}, 0.0)};
}
 
std::string format_measurements(const Particles &particles,
                                uint64_t scatterings_this_interval,
                                const QuantumNumbers &conserved_initial,
                                SystemTimePoint time_start, double time,
                                double E_mean_field,
                                double E_mean_field_initial) {
  const SystemTimeSpan elapsed_seconds = SystemClock::now() - time_start;
 
  const QuantumNumbers current_values(particles);
  const QuantumNumbers difference = current_values - conserved_initial;
 
  // Make sure there are no FPEs in case of IC output, were there will
  // eventually be no more particles in the system
  const double current_energy =
      (particles.size() > 0) ? current_values.momentum().x0() : 0.0;
  const double energy_per_part =
      (particles.size() > 0)
          ? (current_energy + E_mean_field) / particles.size()
          : 0.0;
 
  std::ostringstream ss;
  // clang-format off
  ss << field<7, 3> << time
    // total kinetic energy in the system
     << field<11, 3> << current_energy
    // total mean field energy in the system
     << field<11, 3> << E_mean_field
    // total energy in the system
     << field<12, 3> << current_energy + E_mean_field
    // total energy per particle in the system
     << field<12, 6> << energy_per_part;
    // change in total energy per particle (unless IC output is enabled)
    if (particles.size() == 0) {
     ss << field<13, 6> << "N/A";
    } else {
     ss << field<13, 6> << (difference.momentum().x0()
                            + E_mean_field - E_mean_field_initial)
                            / particles.size();
    }
    ss << field<14, 3> << scatterings_this_interval
     << field<10, 3> << particles.size()
     << field<9, 3> << elapsed_seconds;
  // clang-format on
  return ss.str();
}
 
double calculate_mean_field_energy(
    const Potentials &potentials,
    RectangularLattice<smash::DensityOnLattice> &jmu_B_lat,
    const ExperimentParameters &parameters) {
  // basic parameters and variables
  const double V_cell = (jmu_B_lat.cell_sizes())[0] *
                        (jmu_B_lat.cell_sizes())[1] *
                        (jmu_B_lat.cell_sizes())[2];
 
  double E_mean_field = 0.0;
  double density_mean = 0.0;
  double density_variance = 0.0;
 
  /*
   * We anticipate having other options, like the vector DFT potentials, in the
   * future, hence we include checking which potentials are used.
   */
  if (potentials.use_skyrme()) {
    /*
     * Calculating the symmetry energy contribution to the total mean field
     * energy in the system is not implemented at this time.
     */
    if (potentials.use_symmetry() &&
        parameters.outputclock->current_time() == 0.0) {
      logg[LExperiment].warn()
          << "Note:"
          << "\nSymmetry energy is not included in the mean field calculation."
          << "\n\n";
    }
 
    /*
     * Skyrme potential parameters:
     * C1GeV are the Skyrme coefficients converted to GeV,
     * b1 are the powers of the baryon number density entering the expression
     * for the energy density of the system. Note that these exponents are
     * larger by 1 than those for the energy of a particle (which are used in
     * Potentials class). The formula for a total mean field energy due to a
     * Skyrme potential is E_MF = \sum_i (C_i/b_i) ( n_B^b_i )/( n_0^(b_i - 1) )
     * where nB is the local rest frame baryon number density and n_0 is the
     * saturation density. Then the single particle potential follows from
     * V = d E_MF / d n_B .
     */
    double C1GeV = (potentials.skyrme_a()) / 1000.0;
    double C2GeV = (potentials.skyrme_b()) / 1000.0;
    double b1 = 2.0;
    double b2 = (potentials.skyrme_tau()) + 1.0;
 
    /*
     * Note: calculating the mean field only works if lattice is used.
     * We iterate over the nodes of the baryon density lattice to sum their
     * contributions to the total mean field.
     */
    int number_of_nodes = 0;
    double lattice_mean_field_total = 0.0;
 
    for (auto &node : jmu_B_lat) {
      number_of_nodes++;
      // the rest frame density
      double nB = node.density();
      // the computational frame density
      const double j0 = node.jmu_net().x0();
 
      const double abs_nB = std::abs(nB);
      if ((abs_nB < really_small) || (std::abs(j0) < really_small)) {
        continue;
      }
      density_mean += j0;
      density_variance += j0 * j0;
 
      /*
       * The mean-field energy for the Skyrme potential. Note: this expression
       * is only exact in the rest frame, and is expected to significantly
       * deviate from the correct value for systems that are considerably
       * relativistic. Note: symmetry energy is not taken into the account.
       *
       * TODO: Add symmetry energy.
       */
      double mean_field_contribution_1 = (C1GeV / b1) * std::pow(abs_nB, b1) /
                                         std::pow(nuclear_density, b1 - 1);
      double mean_field_contribution_2 = (C2GeV / b2) * std::pow(abs_nB, b2) /
                                         std::pow(nuclear_density, b2 - 1);
 
      lattice_mean_field_total +=
          V_cell * (mean_field_contribution_1 + mean_field_contribution_2);
    }
 
    // logging statistical properties of the density calculation
    density_mean = density_mean / number_of_nodes;
    density_variance = density_variance / number_of_nodes;
    double density_scaled_variance =
        std::sqrt(density_variance - density_mean * density_mean) /
        density_mean;
    logg[LExperiment].debug() << "\t\t\t\t\t";
    logg[LExperiment].debug()
        << "\n\t\t\t\t\t            density mean = " << density_mean;
    logg[LExperiment].debug()
        << "\n\t\t\t\t\t density scaled variance = " << density_scaled_variance;
    logg[LExperiment].debug()
        << "\n\t\t\t\t\t        total mean_field = "
        << lattice_mean_field_total * parameters.testparticles << "\n";
 
    /*
     * E_mean_field is multiplied by the number of testparticles because the
     * total kinetic energy tracked is that of all particles, including
     * testparticles, and so this is more consistent with the general paradigm.
     */
    E_mean_field = lattice_mean_field_total * parameters.testparticles;
  }
 
  return E_mean_field;
}
 
EventInfo fill_event_info(const Particles &particles, double E_mean_field,
                          double modus_impact_parameter,
                          const ExperimentParameters &parameters,
                          bool projectile_target_interact) {
  const QuantumNumbers current_values(particles);
  const double E_kinetic_total = current_values.momentum().x0();
  const double E_total = E_kinetic_total + E_mean_field;
 
  EventInfo event_info{modus_impact_parameter,
                       parameters.box_length,
                       parameters.outputclock->current_time(),
                       E_kinetic_total,
                       E_mean_field,
                       E_total,
                       parameters.testparticles,
                       !projectile_target_interact};
  return event_info;
}
 
}  // namespace smash