jamimpl.cpp

#include "config.hpp"
#include "jamimpl.hpp"
#include <cstdlib>
#include <vector>
#include <string>
#include <memory>

#ifdef USE_LIBJAM2
# include <sstream>
# include <fstream>
# include <Pythia8/Settings.h>
# include <libjam2.hpp>
#endif

#ifdef USE_LIBJAM1
# include "ParticleSample.hpp"
# include "libjam1.hpp"
#endif

namespace idt {
namespace runjam {

#ifdef USE_LIBJAM2
  static libjam2::irunner* runner_instance = nullptr;
#endif

  int getParticleStableCode(int kf) {
#ifdef USE_LIBJAM2
    if (runner_instance)
      return runner_instance->get_particle_stable_code(kf);
#endif

#if defined(USE_LIBJAM1)
    return libjam1::determineStableCode(kf);
#else
    switch (std::abs(kf)) {
    case 111: case 211: case 221: // pions
    case 311: case 321: // kaons
    case 2212: case 2112: // nucleons
    case 3112: case 3122: case 3212: case 3222: // Lambda and Sigma
    case 3312: case 3322: case 3334: // Xi and Omega
      return 1;
    default:
      return 2;
    }
#endif
  }

#if defined(USE_LIBJAM2) || defined(USE_LIBJAM1)
  double getParticleMass(int kf) {
# ifdef USE_LIBJAM2
    if (runner_instance)
      return runner_instance->get_particle_mass(kf);
# endif

# ifdef USE_LIBJAM1
    return libjam1::jamMass(kf);
# else
    std::cerr << "FATAL: libjam2::runner not yet initialized." << std::endl;
    std::exit(3);
# endif
  }
#endif

  class NullRunner: public IJamRunner {
    std::string message;

    void abort() const {
      std::cerr << message << std::endl;
      std::exit(3);
    }

  public:
    NullRunner(const char* message): message(message) {}

    void initialize(runjam_context const& ctx, std::string const& cascadeMode) override {
      (void) ctx;
      (void) cascadeMode;
      this->abort();
    }
    void finalize() override {
      this->abort();
    }
    void adjust_mass(ParticleSampleBase& psample) const override {
      (void) psample;
      this->abort();
    }
    void do_decay(ParticleSampleBase& psamp, std::vector<Particle>& final_state) override {
      (void) psamp;
      (void) final_state;
      this->abort();
    }
    void do_cascade(ParticleSampleBase& psamp, std::vector<Particle>& final_state, int iev = 0) override {
      (void) psamp;
      (void) final_state;
      (void) iev;
      this->abort();
    }
    double get_average_collision_number() const override {
      this->abort();
      return 0.0;
    }
  };

#ifdef USE_LIBJAM2

  class Jam2Runner: public IJamRunner {
    std::unique_ptr<libjam2::irunner> m_runner;
    double m_oversample = 1.0;

  public:
    void initialize(runjam_context const& ctx, std::string const& cascadeMode) override {
      std::string const fconfig = ctx.get_config<std::string>("runjam_jam2_config", "/dev/null");
      m_runner = libjam2::create_runner(ctx, fconfig);
      m_oversample = ctx.get_config("runjam_oversampling_factor", 1.0);
      if (m_oversample != std::floor(m_oversample)) {
        std::cerr << "jam2: 'runjam_oversampling_factor=" << m_oversample << "' needs to be an integer." << std::endl;
        std::exit(2);
      }

      auto const settings = m_runner->settings();

      settings->flag("Random:setSeed", true);
      settings->mode("Random:seed", ctx.seed());
      settings->mode("Cascade:model", cascadeMode == "decay" ? 0 : 3);
      settings->mode("Cascade:overSample", (int) m_oversample);

      std::string const cachedir = ctx.cachedir();
      settings->word("Cascade:dataFile1", cachedir + "/BWintjam2a.dat");
      settings->word("Cascade:dataFile2", cachedir + "/BWintjam2b.dat");
      settings->word("Cascade:bwFileName1", cachedir + "/BWintjam2a.dat");
      settings->word("Cascade:bwFileName2", cachedir + "/BWintjam2b.dat");

      // weak decay
      {
        bool const weakdecay = ctx.get_config("runjam_switch_weak_decay", false);
        std::string const setting = weakdecay ? ":mayDecay = on" : ":mayDecay = off";

        settings->readString("111" + setting); // pi0
        settings->readString("311" + setting); // k0
        settings->readString("-311" + setting); // ak0
        settings->readString("310" + setting); // k0_S
        settings->readString("130" + setting); // k0_L

        settings->readString("411" + setting);  // D+
        settings->readString("421" + setting);  // D0
        settings->readString("221" + setting);  // eta
        settings->readString("331" + setting);  // eta'
        settings->readString("441" + setting);  // xeta_c
        settings->readString("310" + setting);
        settings->readString("431" + setting);
        settings->readString("511" + setting);
        settings->readString("521" + setting);
        settings->readString("531" + setting);
        settings->readString("3122" + setting);
        settings->readString("3112" + setting);
        settings->readString("3212" + setting);
        settings->readString("3222" + setting);
        settings->readString("3312" + setting);
        settings->readString("3322" + setting);
        settings->readString("3334" + setting);
      }

      // phi decay
      settings->flag("333:mayDecay", ctx.get_config("runjam_phi_decays", true));

      m_runner->initialize();
      if (runner_instance == nullptr)
        runner_instance = m_runner.get();
    }

    void finalize() override {
      m_runner.reset();
    }

    void adjust_mass(ParticleSampleBase& psample) const override {
      for (Particle& part: psample) {
        if (part.pdg == 0) continue;
        double const px = part.px;
        double const py = part.py;
        double const pz = part.pz;
        double const e = part.e;
        if (e < 0.0) {
          double const m = m_runner->get_particle_mass(part.pdg);
          part.mass = m;
          part.e = std::sqrt(px * px + py * py + pz * pz + m * m);
        } else {
          part.mass = std::sqrt(e * e - (px * px + py * py + pz * pz));
        }
      }
    }

    void do_cascade(ParticleSampleBase& psamp, std::vector<Particle>& final_state, int iev) override {
      if (psamp.getOverSamplingFactor() != m_oversample) {
        std::cerr << "ParticleSample: mismatching oversampling factor between the initial condition and JAM2 initialization." << std::endl;
        std::exit(1);
      }

      (void) iev;
      std::vector<Particle> initial_state(psamp.begin(), psamp.end());
      m_runner->run(initial_state, final_state);
    }

    void do_decay(ParticleSampleBase& psamp, std::vector<Particle>& final_state) override {
      // settings に指定している。
      this->do_cascade(psamp, final_state, 0);
    }

    double get_average_collision_number() const override {
      return m_runner->get_event_collision_number();
    }
  };

  std::unique_ptr<IJamRunner> create_jam2_runner() {
    return std::make_unique<Jam2Runner>();
  }

#else

  std::unique_ptr<IJamRunner> create_jam2_runner() {
    return std::make_unique<NullRunner>("FATAL: runjam: compiled without the JAM2 support.");
  }

#endif


#ifdef USE_LIBJAM1
  static void storeParticlesInJam1(Particle const* begin, Particle const* end) {
    int nv = 0;
    int nbary = 0;
    int nmeson = 0;

    for (; begin != end; ++begin) {
      Particle const& particle = *begin;
      int const kf = particle.pdg;
      if (kf == 0) continue;

      int const kc = libjam1::jamComp(kf);      // internal particle code.
      int const ks = libjam1::determineStableCode(kf);
      int const ibary = libjam1::getBaryonNumber(kc, kf);  // baryon number
      if (ibary == 0)
        nmeson++;
      else
        nbary++;
      nv++;

      //...Zero the vector.
      libjam1::jamZero(nv);
      libjam1::setK(1,  nv, ks);
      libjam1::setK(2,  nv, kf);
      libjam1::setK(3,  nv, 0);
      libjam1::setK(4,  nv, 0);
      libjam1::setK(5,  nv, -1);
      libjam1::setK(6,  nv, 0);
      libjam1::setK(7,  nv, 1);
      libjam1::setK(8,  nv, 1);
      libjam1::setK(9,  nv, ibary);
      libjam1::setK(10, nv, 0);
      libjam1::setK(11, nv, 0);

      double const x  = particle.x;
      double const y  = particle.y;
      double const z  = particle.z;
      double const t  = particle.t;
      double const px = particle.px  ;
      double const py = particle.py  ;
      double const pz = particle.pz  ;
      double const pe = particle.e   ;
      double const pm = particle.mass;

      libjam1::setP(1, nv, px);
      libjam1::setP(2, nv, py);
      libjam1::setP(3, nv, pz);
      libjam1::setP(4, nv, pe);
      libjam1::setP(5, nv, pm);

      libjam1::setR(1, nv, x);
      libjam1::setR(2, nv, y);
      libjam1::setR(3, nv, z);
      libjam1::setR(4, nv, t);
      libjam1::setR(5, nv, t); // formation time
      libjam1::setV(1, nv, x); // vertex
      libjam1::setV(2, nv, y); // vertex
      libjam1::setV(3, nv, z); // vertex
      libjam1::setV(4, nv, t); // vertex

      // Set resonance decay time.
      double decayTime = 1.0e+35;
      if (ks == 2) decayTime = t + libjam1::jamDecayTime(1, kf, kc, ks, pm, pe);
      libjam1::setV(5, nv, decayTime);
    }

    libjam1::setNV(nv);          // set total number of particles
    libjam1::setNBARY(nbary);    // set total number of baryons
    libjam1::setNMESON(nmeson);  // set total number of mesons
  }

  static int loadParticlesFromJam1(std::vector<Particle>& list) {
    list.clear();
    int const nv = libjam1::getNV();
    int const ntest = libjam1::getMSTC(5);
    for (int i = 1; i <= nv; i++) {
      //if(libjam1::getK(1,i) > 10) continue;
      list.emplace_back();
      Particle& part = list.back();
      part.pdg  = libjam1::getK(2, i);
      part.mass = libjam1::getP(5, i);
      part.px = libjam1::getP(1, i);
      part.py = libjam1::getP(2, i);
      part.pz = libjam1::getP(3, i);
      part.e  = libjam1::getP(4, i);
      part.x  = libjam1::getR(1, i);
      part.y  = libjam1::getR(2, i);
      part.z  = libjam1::getR(3, i);
      part.t  = libjam1::getR(4, i);
    }
    return ntest;
  }

  class Jam1Runner: public IJamRunner {
  public:
    void initialize(runjam_context const& ctx, std::string const& cascadeMode) override {
      (void) cascadeMode;

      int         const cfg_nevent     = ctx.nevent(1);
      std::string const cfg_outdir     = ctx.outdir();
      int         const cfg_jamseed    = ctx.get_config("runjam_jamseed", ctx.seed());
      bool        const cfg_phi_decays = ctx.get_config("runjam_phi_decays", true);
      bool        const cfg_weakdecay  = ctx.get_config("runjam_switch_weak_decay", false);

      // Initialize JAM
      libjam1::setMSTC(  1, cfg_jamseed); // int seed = 1921;
      libjam1::setMSTC(  2, cfg_nevent);  // number of event
      libjam1::setPARC(  6, 5.0); // scale of display
      libjam1::setMSTC(  8,   0); // job mode
      libjam1::setMSTC( 16,   0); // display on/off
      //libjam1::setMSTC( 38,   6); // IO number for jamlist
      //libjam1::setMSTC( 39,   0); //no output fname(4)
      libjam1::setMSTC( 54,   0); // avoid first coll inside the same nucleus off

      // Default settings
      //libjam1::setMSTC(156,   0); // analysis of collision distribution
      //libjam1::setMSTC(161,   0); // no analysis from jam internal subr
      //libjam1::setMSTC(162,   0); // Output collision history
      //libjam1::setMSTC(165,   0); //

      // Additional settings
      libjam1::setMSTC(  4, 100); // user defined frame
      libjam1::setPARC(  7, 1.0); // Output time interval (fm/c)
      //libjam1::setMSTC( 41,   0); // 0:no resonance decay after simulation
      //libjam1::setMSTC( 61,   0); // isotropic resonance decay option
      libjam1::setMSTC( 81,   0); // 1:hard scattering on/off
      libjam1::setMSTC(156,   1); // analysis of collision distribution
      libjam1::setMSTC(161,   0); // no analysis from jam internal subr
      libjam1::setMSTC(162,   1); // Output collision history
      libjam1::setMSTC(165,   1); //

      // JAM*.DAT files
      if (cfg_outdir.size() > 0) {
        std::string const dir = cfg_outdir + "/";
        libjam1::setFNAME(2, dir + libjam1::getFNAME(2));
        libjam1::setFNAME(3, dir + libjam1::getFNAME(3));
        libjam1::setFNAME(4, dir + libjam1::getFNAME(4));
        libjam1::setFNAME(8, dir);
      }

      //libjam1::setMDCY(libjam1::jamComp(111) , 1, 0); // no pi0 decay
      //libjam1::setMDCY(libjam1::jamComp(3122), 1, 1); // Lambda decay
      //libjam1::setMDCY(libjam1::jamComp(3222), 1, 1); // Sigma- decay
      //libjam1::setMDCY(libjam1::jamComp(3212), 1, 1); // Sigma0 decay
      //libjam1::setMDCY(libjam1::jamComp(3112), 1, 1); // Sigma+ decay
      if (!cfg_phi_decays)
        libjam1::setMDCY(libjam1::jamComp(333), 1, 0); // no phi decay

      if (cfg_weakdecay) libjam1::setMSTC(42, 0); // allow weak decay

      // JAMINIT
      std::string frame = "user"; // comp. frame in this case, user defined
      double dt = 100.0;          // collision time(fm/c)
      int nstep = 1;              // time step (i.e. no time step)
      double bmin = 0.0;         // minimum impact parameter (dummy)
      double bmax = 0.0;         // maximum impact parameter (dummy)
      libjam1::setPARD(16, 10.0); // user defined frame.
      libjam1::jamInit(cfg_nevent, bmin, bmax, dt, nstep, frame.c_str(), "p ", "p ", "2gev");
    }

    void finalize() override {
      libjam1::jamFin();
    }

    void adjust_mass(ParticleSampleBase& psample) const override {
      for (Particle& part: psample) {
        if (part.pdg == 0) continue;
        double const px = part.px;
        double const py = part.py;
        double const pz = part.pz;
        double const e = part.e;
        if (e < 0.0) {
          double const m = runjam::getParticleMass(part.pdg);
          part.mass = m;
          part.e = std::sqrt(px * px + py * py + pz * pz + m * m);
        } else {
          part.mass = std::sqrt(e * e - (px * px + py * py + pz * pz));
        }
      }
    }

    void do_decay(ParticleSampleBase& psamp, std::vector<Particle>& final_state) override {
      double const ntest = psamp.getOverSamplingFactor();
      libjam1::setMSTC(5, ntest);
      storeParticlesInJam1(psamp.begin(), psamp.end());
      libjam1::finalResonanceDecay();
      loadParticlesFromJam1(final_state);
    }

    void do_cascade(ParticleSampleBase& psamp, std::vector<Particle>& final_state, int iev = 0) override {
      double const ntest = psamp.getOverSamplingFactor();
      libjam1::setMSTC(5, ntest);
      storeParticlesInJam1(psamp.begin(), psamp.end());
      libjam1::jamEvt(iev);
      loadParticlesFromJam1(final_state);
    }

    double get_average_collision_number() const override {
      int const ntest = libjam1::getMSTC(5);
      return (libjam1::getMSTD(41) + libjam1::getMSTD(42)) / ntest;
    }
  };

  std::unique_ptr<IJamRunner> create_jam1_runner() {
    return std::make_unique<Jam1Runner>();
  }

#else

  std::unique_ptr<IJamRunner> create_jam1_runner() {
    return std::make_unique<NullRunner>("FATAL: runjam: compiled without the JAM1 support.");
  }

#endif

  std::unique_ptr<IJamRunner> create_runner(runjam_context const& ctx) {
#if defined(USE_LIBJAM1) || defined(USE_LIBJAM2)
    int const libjam_version = ctx.get_config("runjam_jam_version", DEFAULT_JAM_VERSION);
    switch (libjam_version) {
    case 1:
#ifdef USE_LIBJAM1
      return runjam::create_jam1_runner();
#else
      std::cerr << "runjam was compiled without the support for JAM1" << std::endl;
      break;
#endif
    case 2:
#ifdef USE_LIBJAM2
      return runjam::create_jam2_runner();
#else
      std::cerr << "runjam was compiled without the support for JAM2" << std::endl;
      break;
#endif
    default:
      std::cerr << "runjam: invalid value runjam_jam_version='" << libjam_version << "'" << std::endl;
      break;
    }
#else
    std::cerr << "  runjam was compiled without the support for JAM" << std::endl;
#endif
    return std::unique_ptr<IJamRunner>();
  }

  std::vector<std::string> libjam_versions() {
    std::vector<std::string> ret;
#ifdef USE_LIBJAM1
    ret.push_back("libjam1");
#endif
#ifdef USE_LIBJAM2
    ret.push_back(libjam2::version_string());
#endif
    return ret;
  }
}
}