processstring.h

Go to the documentation of this file.

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

#ifndef SRC_INCLUDE_PROCESSSTRING_H_
#define SRC_INCLUDE_PROCESSSTRING_H_

#include <memory>
#include <string>
#include <utility>
#include <vector>

#include "Pythia8/Pythia.h"

#include "constants.h"
#include "logging.h"
#include "particledata.h"

namespace smash {

// \todo Sangwook: make file (processstring) and class (StringProcess) name
// consistent

class StringProcess {
 private:
  // The following 4 variables are in the center of mass frame
  double PPosA_;
  double PPosB_;
  double PNegA_;
  double PNegB_;
  double massA_;
  double massB_;
  double sqrtsAB_;
  std::array<PdgCode, 2> PDGcodes_;
  std::array<FourVector, 2> plab_;
  std::array<FourVector, 2> pcom_;
  FourVector ucomAB_;
  ThreeVector vcomAB_;
  std::array<ThreeVector, 3> evecBasisAB_;
  int NpartFinal_;
  std::array<int, 2> NpartString_;
  double pmin_gluon_lightcone_;
  double pow_fgluon_beta_;
  double pow_fquark_alpha_;
  double pow_fquark_beta_;
  double sigma_qperp_;
  double leading_frag_mean_;
  double leading_frag_width_;
  double kappa_tension_string_;
  double additional_xsec_supp_;
  double time_formation_const_;
  double soft_t_form_;
  double time_collision_;
  bool use_yoyo_model_;
  double prob_proton_to_d_uu_;
  double sqrt2_;

  bool pythia_parton_initialized_ = false;

  ParticleList final_state_;

  std::unique_ptr<Pythia8::Pythia> pythia_parton_;

  std::unique_ptr<Pythia8::Pythia> pythia_hadron_;

  Pythia8::SigmaTotal pythia_sigmatot_;

  Pythia8::Event event_intermediate_;

 public:
  // clang-format off

  StringProcess(double string_tension, double time_formation,
                double gluon_beta, double gluon_pmin,
                double quark_alpha, double quark_beta,
                double strange_supp, double diquark_supp,
                double sigma_perp, double leading_frag_mean,
                double leading_frag_width, double stringz_a,
                double stringz_b,  double string_sigma_T,
                double factor_t_form, bool use_yoyo_model,
                double prob_proton_to_d_uu);

  void common_setup_pythia(Pythia8::Pythia *pythia_in, double strange_supp,
                           double diquark_supp, double stringz_a,
                           double stringz_b, double string_sigma_T);

  void init_pythia_hadron_rndm() {
    const auto &log = logger<LogArea::Pythia>();
    const int seed_new =
        random::uniform_int(1, maximum_rndm_seed_in_pythia);

    pythia_hadron_->rndm.init(seed_new);
    log.debug("pythia_hadron_ : rndm is initialized with seed ", seed_new);
  }

  // clang-format on

  Pythia8::Pythia *get_ptr_pythia_parton() { return pythia_parton_.get(); }

  std::array<double, 3> cross_sections_diffractive(int pdg_a, int pdg_b,
                                                   double sqrt_s) {
    // This threshold magic is following Pythia. Todo(ryu): take care of this.
    double sqrts_threshold = 2. * (1. + 1.0e-6);
    /* In the case of mesons, the corresponding vector meson masses
     * are used to evaluate the energy threshold. */
    const int pdg_a_mod =
        (std::abs(pdg_a) > 1000) ? pdg_a : 10 * (std::abs(pdg_a) / 10) + 3;
    const int pdg_b_mod =
        (std::abs(pdg_b) > 1000) ? pdg_b : 10 * (std::abs(pdg_b) / 10) + 3;
    sqrts_threshold += pythia_hadron_->particleData.m0(pdg_a_mod) +
                       pythia_hadron_->particleData.m0(pdg_b_mod);
    /* Constant cross-section for sub-processes below threshold equal to
     * cross-section at the threshold. */
    if (sqrt_s < sqrts_threshold) {
      sqrt_s = sqrts_threshold;
    }
    pythia_sigmatot_.calc(pdg_a, pdg_b, sqrt_s);
    return {pythia_sigmatot_.sigmaAX(), pythia_sigmatot_.sigmaXB(),
            pythia_sigmatot_.sigmaXX()};
  }

  void set_pmin_gluon_lightcone(double p_light_cone_min) {
    pmin_gluon_lightcone_ = p_light_cone_min;
  }
  void set_pow_fgluon(double betapow) { pow_fgluon_beta_ = betapow; }
  void set_pow_fquark(double alphapow, double betapow) {
    pow_fquark_alpha_ = alphapow;
    pow_fquark_beta_ = betapow;
  }
  void set_sigma_qperp_(double sigma_qperp) { sigma_qperp_ = sigma_qperp; }
  void set_tension_string(double kappa_string) {
    kappa_tension_string_ = kappa_string;
  }

  // clang-format off

  void init(const ParticleList &incoming, double tcoll);
  static void make_orthonormal_basis(ThreeVector &evec_polar,
                                     std::array<ThreeVector, 3> &evec_basis);
  void compute_incoming_lightcone_momenta();
  bool set_mass_and_direction_2strings(
      const std::array<std::array<int, 2>, 2> &quarks,
      const std::array<FourVector, 2> &pstr_com,
      std::array<double, 2> &m_str,
      std::array<ThreeVector, 2> &evec_str);
  bool make_final_state_2strings(
      const std::array<std::array<int, 2>, 2> &quarks,
      const std::array<FourVector, 2> &pstr_com,
      const std::array<double, 2> &m_str,
      const std::array<ThreeVector, 2> &evec_str,
      bool flip_string_ends, bool separate_fragment_baryon);

  bool next_SDiff(bool is_AB_to_AX);
  bool next_DDiff();
  bool next_NDiffSoft();
  bool next_NDiffHard();
  bool next_BBbarAnn();

  static void find_excess_constituent(PdgCode &pdg_actual, PdgCode &pdg_mapped,
                                      std::array<int, 5> &excess_quark,
                                      std::array<int, 5> &excess_antiq);
  void replace_constituent(Pythia8::Particle &particle,
                           std::array<int, 5> &excess_constituent);

  void find_total_number_constituent(Pythia8::Event &event_intermediate,
                                     std::array<int, 5> &nquark_total,
                                     std::array<int, 5> &nantiq_total);

  bool splitting_gluon_qqbar(Pythia8::Event &event_intermediate,
      std::array<int, 5> &nquark_total, std::array<int, 5> &nantiq_total,
      bool sign_constituent,
      std::array<std::array<int, 5>, 2> &excess_constituent);

  void rearrange_excess(std::array<int, 5> &nquark_total,
                        std::array<std::array<int, 5>, 2> &excess_quark,
                        std::array<std::array<int, 5>, 2> &excess_antiq);

  bool restore_constituent(Pythia8::Event &event_intermediate,
                           std::array<std::array<int, 5>, 2> &excess_quark,
                           std::array<std::array<int, 5>, 2> &excess_antiq);

  void compose_string_parton(bool find_forward_string,
                             Pythia8::Event &event_intermediate,
                             Pythia8::Event &event_hadronize);
  void compose_string_junction(bool &find_forward_string,
                               Pythia8::Event &event_intermediate,
                               Pythia8::Event &event_hadronize);

  void find_junction_leg(bool sign_color, std::vector<int> &col,
                         Pythia8::Event &event_intermediate,
                         Pythia8::Event &event_hadronize);

  int get_index_forward(bool find_forward, int np_end,
                        Pythia8::Event &event) {
    int iforward = 1;
    for (int ip = 2; ip < event.size() - np_end; ip++) {
      const double y_quark_current = event[ip].y();
      const double y_quark_forward = event[iforward].y();
      if ((find_forward && y_quark_current > y_quark_forward) ||
          (!find_forward && y_quark_current < y_quark_forward)) {
        iforward = ip;
      }
    }
    return iforward;
  }

  ParticleList get_final_state() { return final_state_; }

  int append_final_state(ParticleList &intermediate_particles,
                         const FourVector &uString,
                         const ThreeVector &evecLong);

  static bool append_intermediate_list(int pdgid, FourVector momentum,
                                       ParticleList &intermediate_particles) {
    const std::string s = std::to_string(pdgid);
    PdgCode pythia_code(s);
    ParticleData new_particle(ParticleType::find(pythia_code));
    new_particle.set_4momentum(momentum);
    intermediate_particles.push_back(new_particle);
    return true;
  }

  static void convert_KaonLS(int &pythia_id) {
    if (pythia_id == 310 || pythia_id == 130) {
      pythia_id = (random::uniform_int(0, 1) == 0) ? 311 : -311;
    }
  }

  static void quarks_from_diquark(int diquark, int &q1, int &q2, int &deg_spin);

  static int diquark_from_quarks(int q1, int q2);

  static void make_string_ends(const PdgCode &pdgcode_in, int &idq1, int &idq2, double xi);

  static Pythia8::Vec4 set_Vec4(double energy, const ThreeVector &mom) {
    return Pythia8::Vec4(mom.x1(), mom.x2(), mom.x3(), energy);
  }

  static FourVector reorient(Pythia8::Particle &particle,
                             std::array<ThreeVector, 3> &evec_basis) {
    ThreeVector three_momentum = evec_basis[0] * particle.pz() +
                                 evec_basis[1] * particle.px() +
                                 evec_basis[2] * particle.py();
    return FourVector(particle.e(), three_momentum);
  }

  int fragment_string(int idq1, int idq2, double mString,
                      ThreeVector &evecLong, bool flip_string_ends,
                      bool separate_fragment_baryon,
                      ParticleList &intermediate_particles);

  static void assign_all_scaling_factors(int baryon_string,
                                         ParticleList& outgoing_particles,
                                         const ThreeVector &evecLong,
                                         double suppression_factor);

  static std::pair<int, int> find_leading(int nq1, int nq2, ParticleList& list);

  static void assign_scaling_factor(int nquark, ParticleData& data,
                                    double suppression_factor);

  static int pdg_map_for_pythia(PdgCode &pdg);

  // clang-format on
};

}  // namespace smash

#endif  // SRC_INCLUDE_PROCESSSTRING_H_