oscaroutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2024
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/oscaroutput.h"
 
#include <filesystem>
#include <string>
 
#include "smash/action.h"
#include "smash/clock.h"
#include "smash/config.h"
#include "smash/forwarddeclarations.h"
 
namespace smash {
static constexpr int LHyperSurfaceCrossing = LogArea::HyperSurfaceCrossing::id;
 
template <OscarOutputFormat Format, int Contents>
OscarOutput<Format, Contents>::OscarOutput(
    const std::filesystem::path &path, const std::string &name,
    const std::vector<std::string> quantities)
    : OutputInterface(name),
      file_{path / (name + ((Format == ASCII) ? ".dat" : ".oscar") +
                    ((Format == OscarFormat1999) ? "1999" : "")),
            "w"},
      formatter_{Format == ASCII ? quantities
                 : (Format == OscarFormat2013)
                     ? OutputDefaultQuantities::oscar2013
                 : (Format == OscarFormat2013Extended)
                     ? OutputDefaultQuantities::oscar2013extended
                     : OutputDefaultQuantities::oscar1999} {
  if (Format != ASCII && !quantities.empty()) {
    throw std::logic_error(
        "Non-empty Quantities given alongside format other than ASCII.");
  }
  std::string format_name;
  if (Format == ASCII) {
    format_name = "ASCII";
  } else if (Format == OscarFormat2013) {
    format_name = "OSCAR2013";
  } else if (Format == OscarFormat2013Extended) {
    format_name = "OSCAR2013Extended";
  } else {
    format_name = "OSC1999A";
  }
  if (Format == ASCII || Format == OscarFormat2013 ||
      Format == OscarFormat2013Extended) {
    std::fprintf(file_.get(), "#!%s %s %s\n", format_name.c_str(), name.c_str(),
                 formatter_.quantities_line().c_str());
    std::fprintf(file_.get(), "# Units: %s\n", formatter_.unit_line().c_str());
    std::fprintf(file_.get(), "# %s\n", SMASH_VERSION);
  } else {
    const std::string &oscar_name =
        name == "particle_lists" ? "final_id_p_x" : name;
    // This is necessary because OSCAR1999A requires
    // this particular string for particle output.
 
    std::fprintf(file_.get(), "# %s\n# %s\n# %s\n", format_name.c_str(),
                 oscar_name.c_str(), SMASH_VERSION);
    std::fprintf(file_.get(), "# Block format:\n");
    if (oscar_name == "full_event_history") {
      std::fprintf(file_.get(),
                   "# nin nout density tot_weight part_weight proc_type\n");
    } else {
      std::fprintf(file_.get(), "# nin nout event_number ensemble_number\n");
    }
    std::fprintf(file_.get(), "# %s\n", formatter_.quantities_line().c_str());
    std::fprintf(
        file_.get(),
        "# End of event: 0 0 event_number ensemble_number impact_parameter\n");
    std::fprintf(file_.get(), "#\n");
  }
}
 
template <OscarOutputFormat Format, int Contents>
inline void OscarOutput<Format, Contents>::write(const Particles &particles) {
  for (const ParticleData &data : particles) {
    write_particledata(data);
  }
}
 
template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_eventstart(const Particles &particles,
                                                  const EventLabel &event_label,
                                                  const EventInfo &) {
  // We do not want the inital particle list or number to be printed in case of
  // IC output
  if (Contents & OscarAtEventstart && !(Contents & OscarParticlesIC)) {
    if (Format == ASCII || Format == OscarFormat2013 ||
        Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(), "# event %i ensemble %i in %zu\n",
                   event_label.event_number, event_label.ensemble_number,
                   particles.size());
    } else {
      /* OSCAR line prefix : initial particles; final particles; event id
       * First block of an event: initial = 0, final = number of particles
       */
      const size_t zero = 0;
      std::fprintf(file_.get(), "%zu %zu %i %i\n", zero, particles.size(),
                   event_label.event_number, event_label.ensemble_number);
    }
    write(particles);
  } else if (Contents & OscarParticlesIC) {
    if (Format == ASCII || Format == OscarFormat2013 ||
        Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(), "# event %i ensemble %i start\n",
                   event_label.event_number, event_label.ensemble_number);
    } else if (Format == OscarFormat1999) {
      const size_t zero = 0;
      std::fprintf(file_.get(), "%zu %zu %i %i\n", zero, zero,
                   event_label.event_number, event_label.ensemble_number);
    }
  }
}
 
template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_eventend(const Particles &particles,
                                                const EventLabel &event_label,
                                                const EventInfo &event) {
  if (Format == ASCII || Format == OscarFormat2013 ||
      Format == OscarFormat2013Extended) {
    if (Contents & OscarParticlesAtEventend ||
        (Contents & OscarParticlesAtEventendIfNotEmpty && !event.empty_event)) {
      std::fprintf(file_.get(), "# event %i ensemble %i out %zu\n",
                   event_label.event_number, event_label.ensemble_number,
                   particles.size());
      write(particles);
    }
    // Comment end of an event
    if (!(Contents & OscarParticlesIC)) {
      const char *empty_event_str = event.empty_event ? "no" : "yes";
      std::fprintf(file_.get(),
                   "# event %i ensemble %i end 0 impact %7.3f "
                   "scattering_projectile_target %s\n",
                   event_label.event_number, event_label.ensemble_number,
                   event.impact_parameter, empty_event_str);
    } else {
      std::fprintf(file_.get(), "# event %i ensemble %i end\n",
                   event_label.event_number, event_label.ensemble_number);
    }
  } else {
    /* OSCAR line prefix : initial particles; final particles; event id
     * Last block of an event: initial = number of particles, final = 0
     * Block ends with null interaction. */
    const size_t zero = 0;
    if (Contents & OscarParticlesAtEventend ||
        (Contents & OscarParticlesAtEventendIfNotEmpty && !event.empty_event)) {
      std::fprintf(file_.get(), "%zu %zu %i %i\n", particles.size(), zero,
                   event_label.event_number, event_label.ensemble_number);
      write(particles);
    }
    // Null interaction marks the end of an event
    std::fprintf(file_.get(), "%zu %zu %i %i %7.3f\n", zero, zero,
                 event_label.event_number, event_label.ensemble_number,
                 event.impact_parameter);
  }
  // Flush to disk
  std::fflush(file_.get());
 
  if (Contents & OscarParticlesIC) {
    // 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).");
    }
  }
}
 
template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_interaction(const Action &action,
                                                   const double density) {
  if (Contents & OscarInteractions) {
    if (Format == ASCII || Format == OscarFormat2013 ||
        Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(),
                   "# interaction in %zu out %zu rho %12.7f weight %12.7g"
                   " partial %12.7f type %5i\n",
                   action.incoming_particles().size(),
                   action.outgoing_particles().size(), density,
                   action.get_total_weight(), action.get_partial_weight(),
                   static_cast<int>(action.get_type()));
    } else {
      /* OSCAR line prefix : initial final
       * particle creation: 0 1
       * particle 2<->2 collision: 2 2
       * resonance formation: 2 1
       * resonance decay: 1 2
       * etc.*/
      std::fprintf(file_.get(), "%zu %zu %12.7f %12.7f %12.7f %5i\n",
                   action.incoming_particles().size(),
                   action.outgoing_particles().size(), density,
                   action.get_total_weight(), action.get_partial_weight(),
                   static_cast<int>(action.get_type()));
    }
    for (const auto &p : action.incoming_particles()) {
      write_particledata(p);
    }
    for (const auto &p : action.outgoing_particles()) {
      write_particledata(p);
    }
  } else if (Contents & OscarParticlesIC) {
    for (const auto &p : action.incoming_particles()) {
      write_particledata(p);
    }
  }
}
 
template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_intermediate_time(
    const Particles &particles, const std::unique_ptr<Clock> &,
    const DensityParameters &, const EventLabel &event_label,
    const EventInfo &) {
  if (Contents & OscarTimesteps) {
    if (Format == ASCII || Format == OscarFormat2013 ||
        Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(), "# event %i ensemble %i out %zu\n",
                   event_label.event_number, event_label.ensemble_number,
                   particles.size());
    } else {
      const size_t zero = 0;
      std::fprintf(file_.get(), "%zu %zu %i %i\n", particles.size(), zero,
                   event_label.event_number, event_label.ensemble_number);
    }
    write(particles);
  }
}
 
template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::write_particledata(
    const ParticleData &data) {
  std::fprintf(file_.get(), "%s\n", formatter_.data_line(data).c_str());
}
 
namespace {
template <int Contents>
std::unique_ptr<OutputInterface> create_select_format(
    bool modern_format, const std::filesystem::path &path,
    const OutputParameters &out_par, const std::string &name,
    const bool custom_format = false) {
  bool extended_format = (Contents & OscarInteractions) ? out_par.coll_extended
                                                        : out_par.part_extended;
  if (custom_format) {
    const auto &quantities = (Contents & OscarInteractions)
                                 ? out_par.quantities.at("Collisions")
                                 : out_par.quantities.at("Particles");
    return std::make_unique<OscarOutput<ASCII, Contents>>(path, name,
                                                          quantities);
  } else if (modern_format && extended_format) {
    return std::make_unique<OscarOutput<OscarFormat2013Extended, Contents>>(
        path, name);
  } else if (modern_format && !extended_format) {
    return std::make_unique<OscarOutput<OscarFormat2013, Contents>>(path, name);
  } else if (!modern_format && !extended_format) {
    return std::make_unique<OscarOutput<OscarFormat1999, Contents>>(path, name);
  } else {
    // Only remaining possibility: (!modern_format && extended_format)
    logg[LOutput].warn() << "Creating Oscar output: "
                         << "There is no extended Oscar1999 format.";
    return std::make_unique<OscarOutput<OscarFormat1999, Contents>>(path, name);
  }
}
}  // unnamed namespace
 
std::unique_ptr<OutputInterface> create_oscar_output(
    const std::string &format, const std::string &content,
    const std::filesystem::path &path, const OutputParameters &out_par) {
  if (format != "Oscar2013" && format != "Oscar1999" && format != "ASCII") {
    throw std::invalid_argument("Creating Oscar output: unknown format");
  }
  const bool modern_format = (format == "Oscar2013");
  const bool custom_format = (format == "ASCII");
  if (content == "Particles") {
    if (out_par.part_only_final == OutputOnlyFinal::Yes) {
      return create_select_format<OscarParticlesAtEventend>(
          modern_format, path, out_par, "particle_lists", custom_format);
    } else if (out_par.part_only_final == OutputOnlyFinal::IfNotEmpty) {
      return create_select_format<OscarParticlesAtEventendIfNotEmpty>(
          modern_format, path, out_par, "particle_lists", custom_format);
 
    } else {  // out_par.part_only_final == OutputOnlyFinal::No
      return create_select_format<OscarTimesteps | OscarAtEventstart |
                                  OscarParticlesAtEventend>(
          modern_format, path, out_par, "particle_lists", custom_format);
    }
  } else if (content == "Collisions") {
    if (out_par.coll_printstartend) {
      return create_select_format<OscarInteractions | OscarAtEventstart |
                                  OscarParticlesAtEventend>(
          modern_format, path, out_par, "full_event_history", custom_format);
    } else {
      return create_select_format<OscarInteractions>(
          modern_format, path, out_par, "full_event_history", custom_format);
    }
  } else if (content == "Dileptons") {
    if (modern_format && out_par.dil_extended) {
      return std::make_unique<
          OscarOutput<OscarFormat2013Extended, OscarInteractions>>(path,
                                                                   "Dileptons");
    } else if (modern_format && !out_par.dil_extended) {
      return std::make_unique<OscarOutput<OscarFormat2013, OscarInteractions>>(
          path, "Dileptons");
    } else if (!modern_format && !out_par.dil_extended) {
      return std::make_unique<OscarOutput<OscarFormat1999, OscarInteractions>>(
          path, "Dileptons");
    } else if (!custom_format && !modern_format && out_par.dil_extended) {
      logg[LOutput].warn()
          << "Creating Oscar output: "
          << "There is no extended Oscar1999 (dileptons) format.";
    }
  } else if (content == "Photons") {
    if (modern_format && !out_par.photons_extended) {
      return std::make_unique<OscarOutput<OscarFormat2013, OscarInteractions>>(
          path, "Photons");
    } else if (modern_format && out_par.photons_extended) {
      return std::make_unique<
          OscarOutput<OscarFormat2013Extended, OscarInteractions>>(path,
                                                                   "Photons");
    } else if (!modern_format && !out_par.photons_extended) {
      return std::make_unique<OscarOutput<OscarFormat1999, OscarInteractions>>(
          path, "Photons");
    } else if (!custom_format && !modern_format && out_par.photons_extended) {
      logg[LOutput].warn()
          << "Creating Oscar output: "
          << "There is no extended Oscar1999 (photons) format.";
    }
  } else if (content == "Initial_Conditions") {
    if (modern_format && !out_par.ic_extended) {
      return std::make_unique<
          OscarOutput<OscarFormat2013, OscarParticlesIC | OscarAtEventstart>>(
          path, "SMASH_IC");
    } else if (modern_format && out_par.ic_extended) {
      return std::make_unique<OscarOutput<
          OscarFormat2013Extended, OscarParticlesIC | OscarAtEventstart>>(
          path, "SMASH_IC");
    } else if (!modern_format && !out_par.ic_extended) {
      return std::make_unique<
          OscarOutput<OscarFormat1999, OscarParticlesIC | OscarAtEventstart>>(
          path, "SMASH_IC");
    } else if (!custom_format && !modern_format && out_par.ic_extended) {
      logg[LOutput].warn()
          << "Creating Oscar output: "
          << "There is no extended Oscar1999 (initial conditions) format.";
    }
  }
 
  throw std::invalid_argument("Create_oscar_output got unknown content.");
}
 
}  // namespace smash