oscaroutput.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
#include "smash/oscaroutput.h"

#include <string>

#include <boost/filesystem.hpp>

#include "smash/action.h"
#include "smash/clock.h"
#include "smash/config.h"
#include "smash/configuration.h"
#include "smash/cxx14compat.h"
#include "smash/forwarddeclarations.h"
#include "smash/particles.h"

namespace smash {

template <OscarOutputFormat Format, int Contents>
OscarOutput<Format, Contents>::OscarOutput(const bf::path &path,
                                           const std::string &name)
    : OutputInterface(name),
      file_{path /
                (name + ".oscar" + ((Format == OscarFormat1999) ? "1999" : "")),
            "w"} {
  if (Format == OscarFormat2013) {
    std::fprintf(file_.get(),
                 "#!OSCAR2013 %s t x y z mass "
                 "p0 px py pz pdg ID charge\n",
                 name.c_str());
    std::fprintf(file_.get(),
                 "# Units: fm fm fm fm "
                 "GeV GeV GeV GeV GeV none none e\n");
    std::fprintf(file_.get(), "# %s\n", VERSION_MAJOR);
  } else if (Format == OscarFormat2013Extended) {
    std::fprintf(file_.get(),
                 "#!OSCAR2013Extended %s t x y z mass p0 px py pz"
                 " pdg ID charge ncoll form_time xsecfac proc_id_origin"
                 " proc_type_origin time_last_coll pdg_mother1 pdg_mother2\n",
                 name.c_str());
    std::fprintf(file_.get(),
                 "# Units: fm fm fm fm GeV GeV GeV GeV GeV"
                 " none none e none fm none none none fm none none\n");
    std::fprintf(file_.get(), "# %s\n", VERSION_MAJOR);
  } else {
    const std::string &oscar_name =
        name == "particle_lists" ? "final_id_p_x" : name;
    // This is necessary because OSCAR199A requires
    // this particular string for particle output.

    std::fprintf(file_.get(), "# OSC1999A\n# %s\n# %s\n", oscar_name.c_str(),
                 VERSION_MAJOR);
    std::fprintf(file_.get(), "# Block format:\n");
    std::fprintf(file_.get(), "# nin nout event_number\n");
    std::fprintf(file_.get(), "# id pdg 0 px py pz p0 mass x y z t\n");
    std::fprintf(file_.get(),
                 "# End of event: 0 0 event_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 int event_number) {
  current_event_ = event_number;
  if (Contents & OscarAtEventstart) {
    if (Format == OscarFormat2013 || Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(), "# event %i in %zu\n", event_number + 1,
                   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\n", zero, particles.size(),
                   event_number + 1);
    }
    if (!(Contents == OscarParticlesIC)) {
      // We do not want the inital particle list to be printed in case of IC
      // output
      write(particles);
    }
  }
}

template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_eventend(const Particles &particles,
                                                const int event_number,
                                                double impact_parameter,
                                                bool empty_event) {
  if (Format == OscarFormat2013 || Format == OscarFormat2013Extended) {
    if (Contents & OscarParticlesAtEventend) {
      std::fprintf(file_.get(), "# event %i out %zu\n", event_number + 1,
                   particles.size());
      write(particles);
    }
    // Comment end of an event
    const char *empty_event_str = empty_event ? "yes" : "no";
    std::fprintf(file_.get(), "# event %i end 0 impact %7.3f empty %s\n",
                 event_number + 1, impact_parameter, empty_event_str);
  } 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) {
      std::fprintf(file_.get(), "%zu %zu %i\n", particles.size(), zero,
                   event_number + 1);
      write(particles);
    }
    // Null interaction marks the end of an event
    std::fprintf(file_.get(), "%zu %zu %i %7.3f\n", zero, zero,
                 event_number + 1, impact_parameter);
  }
  // Flush to disk
  std::fflush(file_.get());

  if (Contents == OscarParticlesIC) {
    const auto &log = logger<LogArea::HyperSurfaceCrossing>();
    // 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).");
    }
  }
}

template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::at_interaction(const Action &action,
                                                   const double density) {
  if (Contents & OscarInteractions) {
    if (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 &) {
  if (Contents & OscarTimesteps) {
    if (Format == OscarFormat2013 || Format == OscarFormat2013Extended) {
      std::fprintf(file_.get(), "# event %i out %zu\n", current_event_ + 1,
                   particles.size());
    } else {
      const size_t zero = 0;
      std::fprintf(file_.get(), "%zu %zu %i\n", particles.size(), zero,
                   current_event_ + 1);
    }
    write(particles);
  }
}

template <OscarOutputFormat Format, int Contents>
void OscarOutput<Format, Contents>::write_particledata(
    const ParticleData &data) {
  const FourVector pos = data.position();
  const FourVector mom = data.momentum();
  if (Format == OscarFormat2013) {
    std::fprintf(file_.get(), "%g %g %g %g %g %.9g %.9g %.9g %.9g %s %i %i\n",
                 pos.x0(), pos.x1(), pos.x2(), pos.x3(), data.effective_mass(),
                 mom.x0(), mom.x1(), mom.x2(), mom.x3(),
                 data.pdgcode().string().c_str(), data.id(),
                 data.type().charge());
  } else if (Format == OscarFormat2013Extended) {
    const auto h = data.get_history();
    std::fprintf(
        file_.get(),
        "%g %g %g %g %g %.9g %.9g %.9g"
        " %.9g %s %i %i %i %g %g %i %i %g %s %s\n",
        pos.x0(), pos.x1(), pos.x2(), pos.x3(), data.effective_mass(), mom.x0(),
        mom.x1(), mom.x2(), mom.x3(), data.pdgcode().string().c_str(),
        data.id(), data.type().charge(), h.collisions_per_particle,
        data.formation_time(), data.xsec_scaling_factor(), h.id_process,
        static_cast<int>(h.process_type), h.time_last_collision,
        h.p1.string().c_str(), h.p2.string().c_str());
  } else {
    std::fprintf(file_.get(), "%i %s %i %g %g %g %g %g %g %g %g %g\n",
                 data.id(), data.pdgcode().string().c_str(), 0, mom.x1(),
                 mom.x2(), mom.x3(), mom.x0(), data.effective_mass(), pos.x1(),
                 pos.x2(), pos.x3(), pos.x0());
  }
}

namespace {
template <int Contents>
std::unique_ptr<OutputInterface> create_select_format(
    bool modern_format, const bf::path &path, const OutputParameters &out_par,
    const std::string &name) {
  const auto &log = logger<LogArea::Output>();
  bool extended_format = (Contents & OscarInteractions) ? out_par.coll_extended
                                                        : out_par.part_extended;
  if (modern_format && extended_format) {
    return make_unique<OscarOutput<OscarFormat2013Extended, Contents>>(path,
                                                                       name);
  } else if (modern_format && !extended_format) {
    return make_unique<OscarOutput<OscarFormat2013, Contents>>(path, name);
  } else if (!modern_format && !extended_format) {
    return make_unique<OscarOutput<OscarFormat1999, Contents>>(path, name);
  } else {
    // Only remaining possibility: (!modern_format && extended_format)
    log.warn() << "Creating Oscar output: "
               << "There is no extended Oscar1999 format.";
    return 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 bf::path &path,
    const OutputParameters &out_par) {
  const auto &log = logger<LogArea::Output>();
  if (format != "Oscar2013" && format != "Oscar1999") {
    throw std::invalid_argument("Creating Oscar output: unknown format");
  }
  const bool modern_format = (format == "Oscar2013");
  if (content == "Particles") {
    if (out_par.part_only_final) {
      return create_select_format<OscarParticlesAtEventend>(
          modern_format, path, out_par, "particle_lists");
    } else {
      return create_select_format<OscarTimesteps | OscarAtEventstart |
                                  OscarParticlesAtEventend>(
          modern_format, path, out_par, "particle_lists");
    }
  } else if (content == "Collisions") {
    if (out_par.coll_printstartend) {
      return create_select_format<OscarInteractions | OscarAtEventstart |
                                  OscarParticlesAtEventend>(
          modern_format, path, out_par, "full_event_history");
    } else {
      return create_select_format<OscarInteractions>(
          modern_format, path, out_par, "full_event_history");
    }
  } else if (content == "Dileptons") {
    if (modern_format && out_par.dil_extended) {
      return make_unique<
          OscarOutput<OscarFormat2013Extended, OscarInteractions>>(path,
                                                                   "Dileptons");
    } else if (modern_format && !out_par.dil_extended) {
      return make_unique<OscarOutput<OscarFormat2013, OscarInteractions>>(
          path, "Dileptons");
    } else if (!modern_format && !out_par.dil_extended) {
      return make_unique<OscarOutput<OscarFormat1999, OscarInteractions>>(
          path, "Dileptons");
    } else if (!modern_format && out_par.dil_extended) {
      log.warn() << "Creating Oscar output: "
                 << "There is no extended Oscar1999 (dileptons) format.";
    }
  } else if (content == "Photons") {
    if (modern_format && !out_par.photons_extended) {
      return make_unique<OscarOutput<OscarFormat2013, OscarInteractions>>(
          path, "Photons");
    } else if (modern_format && out_par.photons_extended) {
      return make_unique<
          OscarOutput<OscarFormat2013Extended, OscarInteractions>>(path,
                                                                   "Photons");
    } else if (!modern_format && !out_par.photons_extended) {
      return make_unique<OscarOutput<OscarFormat1999, OscarInteractions>>(
          path, "Photons");
    } else if (!modern_format && out_par.photons_extended) {
      log.warn() << "Creating Oscar output: "
                 << "There is no extended Oscar1999 (photons) format.";
    }
  } else if (content == "Initial_Conditions") {
    if (modern_format && !out_par.ic_extended) {
      return make_unique<OscarOutput<OscarFormat2013, OscarParticlesIC>>(
          path, "SMASH_IC");
    } else if (modern_format && out_par.ic_extended) {
      return make_unique<
          OscarOutput<OscarFormat2013Extended, OscarParticlesIC>>(path,
                                                                  "SMASH_IC");
    } else if (!modern_format && !out_par.ic_extended) {
      return make_unique<OscarOutput<OscarFormat1999, OscarParticlesIC>>(
          path, "SMASH_IC");
    } else if (!modern_format && out_par.ic_extended) {
      log.warn()
          << "Creating Oscar output: "
          << "There is no extended Oscar1999 (initial conditions) format.";
    }
  }

  throw std::invalid_argument("Create_oscar_output got unknown content.");
}

}  // namespace smash