icoutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2019-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#include "smash/icoutput.h"

#include <fstream>

#include <boost/filesystem.hpp>

#include "smash/action.h"

namespace smash {

ICOutput::ICOutput(const bf::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",
      VERSION_MAJOR);
  std::fprintf(file_.get(), "# tau x y eta mt px py Rap ID charge\n");
  std::fprintf(file_.get(), "# fm fm fm none GeV GeV GeV none none e\n");
}

ICOutput::~ICOutput() {}

void ICOutput::at_eventstart(const Particles &, const int event_number) {
  std::fprintf(file_.get(), "# event %i start\n", event_number + 1);
}

void ICOutput::at_eventend(const Particles &particles, const int event_number,
                           double, bool) {
  const auto &log = logger<LogArea::HyperSurfaceCrossing>();
  std::fprintf(file_.get(), "# event %i end\n", event_number + 1);

  // If the runtime is too short some particles might not yet have
  // reached the hypersurface. Warning is printed.
  if (particles.size() != 0) {
    log.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 &) {
  // 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 = sqrt(particle.type().mass() * particle.type().mass() +
                              particle.momentum()[1] * particle.momentum()[1] +
                              particle.momentum()[2] * particle.momentum()[2]);
  // momentum space rapidity
  const double rapidity =
      0.5 * log((particle.momentum()[0] + particle.momentum()[3]) /
                (particle.momentum()[0] - particle.momentum()[3]));

  // write particle data
  std::fprintf(file_.get(), "%g %g %g %g %g %g %g %g %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.id(), particle.type().charge());

  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