icoutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2019-2020,2022-2024
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/icoutput.h"
 
#include <filesystem>
#include <fstream>
 
#include "smash/action.h"
 
namespace smash {
static constexpr int LHyperSurfaceCrossing = LogArea::HyperSurfaceCrossing::id;
 
ICOutput::ICOutput(const std::filesystem::path &path, const std::string &name,
                   const OutputParameters &out_par)
    : OutputInterface(name),
      file_{path / "SMASH_IC.dat", "w"},
      out_par_(out_par) {
  std::fprintf(
      file_.get(),
      "# %s initial conditions: hypersurface of constant proper time\n",
      SMASH_VERSION);
  std::fprintf(file_.get(),
               "# tau x y eta mt px py Rap pdg charge "
               "baryon_number strangeness\n");
  std::fprintf(file_.get(),
               "# fm fm fm none GeV GeV GeV none none e "
               "none none\n");
}
 
ICOutput::~ICOutput() {}
 
void ICOutput::at_eventstart(const Particles &, const EventLabel &event_label,
                             const EventInfo &event) {
  if (event.n_ensembles != 1) {
    throw std::logic_error(
        "ICOutput shouldn't be used with multiple parallel ensembles.");
  }
  std::fprintf(file_.get(), "# event %i ensemble %i start\n",
               event_label.event_number, event_label.ensemble_number);
}
 
void ICOutput::at_eventend(const Particles &particles,
                           const EventLabel &event_label,
                           const EventInfo &event) {
  if (event.n_ensembles != 1) {
    throw std::logic_error(
        "ICOutput shouldn't be used with multiple parallel ensembles.");
  }
  std::fprintf(file_.get(), "# event %i ensemble %i end\n",
               event_label.event_number, event_label.ensemble_number);
 
  // If the runtime is too short some particles might not yet have
  // reached the hypersurface. Warning is printed.
  if (particles.size() != 0 && !event.impose_kinematic_cut_for_SMASH_IC) {
    logg[LHyperSurfaceCrossing].warn(
        "End time might be too small for initial conditions output. "
        "Hypersurface has not yet been crossed by ",
        particles.size(), " particle(s).");
  }
}
 
void ICOutput::at_intermediate_time(const Particles &,
                                    const std::unique_ptr<Clock> &,
                                    const DensityParameters &,
                                    const EventLabel &, const EventInfo &) {
  // Dummy, but virtual function needs to be declared.
}
 
void ICOutput::at_interaction(const Action &action, const double) {
  assert(action.get_type() == ProcessType::HyperSurfaceCrossing);
  assert(action.incoming_particles().size() == 1);
 
  ParticleData particle = action.incoming_particles()[0];
 
  // transverse mass
  const double m_trans =
      std::sqrt(particle.effective_mass() * particle.effective_mass() +
                particle.momentum()[1] * particle.momentum()[1] +
                particle.momentum()[2] * particle.momentum()[2]);
  // momentum space rapidity
  const double rapidity =
      0.5 * std::log((particle.momentum()[0] + particle.momentum()[3]) /
                     (particle.momentum()[0] - particle.momentum()[3]));
 
  // Determine if particle is spectator:
  // Fulfilled if particle is initial nucleon, aka has no prior interactions
  bool is_spectator = particle.get_history().collisions_per_particle == 0;
 
  // write particle data excluding spectators
  if (!is_spectator) {
    std::fprintf(file_.get(), "%g %g %g %g %g %g %g %g %s %i %i %i \n",
                 particle.position().tau(), particle.position()[1],
                 particle.position()[2], particle.position().eta(), m_trans,
                 particle.momentum()[1], particle.momentum()[2], rapidity,
                 particle.pdgcode().string().c_str(), particle.type().charge(),
                 particle.type().baryon_number(),
                 particle.type().strangeness());
  }
 
  if (IC_proper_time_ < 0.0) {
    // First particle that is removed, overwrite negative default
    IC_proper_time_ = particle.position().tau();
  } else {
    // Verify that all other particles have the same proper time
    const double next_proper_time = particle.position().tau();
    if (!((next_proper_time - IC_proper_time_) < really_small))
      throw std::runtime_error(
          "Hypersurface proper time changed during evolution.");
  }
}
 
}  // namespace smash