stringprocess.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2017-2020,2022
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_STRINGPROCESS_H_
#define SRC_INCLUDE_SMASH_STRINGPROCESS_H_
 
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>
 
#include "Pythia8/Pythia.h"
 
#include "constants.h"
#include "logging.h"
#include "particledata.h"
 
namespace smash {
static constexpr int LPythia = LogArea::Pythia::id;
 
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 stringz_a_leading_;
  double stringz_b_leading_;
  double stringz_a_produce_;
  double stringz_b_produce_;
  double strange_supp_;
  double diquark_supp_;
  double popcorn_rate_;
  double string_sigma_T_;
  double kappa_tension_string_;
  double additional_xsec_supp_;
  double time_formation_const_;
  double soft_t_form_;
  double time_collision_;
  bool mass_dependent_formation_times_;
  double prob_proton_to_d_uu_;
 
  bool separate_fragment_baryon_;
 
  bool use_monash_tune_;
 
  ParticleList final_state_;
 
  typedef std::map<std::pair<int, int>, std::unique_ptr<Pythia8::Pythia>>
      pythia_map;
 
  pythia_map hard_map_;
 
  std::unique_ptr<Pythia8::Pythia> pythia_hadron_;
 
  Pythia8::SigmaTotal pythia_sigmatot_;
 
  Pythia8::StringFlav pythia_stringflav_;
 
  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 stringz_a_leading,
                double stringz_b_leading, double stringz_a,
                double stringz_b,  double string_sigma_T,
                double factor_t_form,
                bool mass_dependent_formation_times,
                double prob_proton_to_d_uu,
                bool separate_fragment_baryon, double popcorn_rate,
                bool use_monash_tune);
 
  void common_setup_pythia(Pythia8::Pythia *pythia_in, double strange_supp,
                           double diquark_supp, double popcorn_rate,
                           double stringz_a, double stringz_b,
                           double string_sigma_T);
 
  void init_pythia_hadron_rndm() {
    const int seed_new =
        random::uniform_int(1, maximum_rndm_seed_in_pythia);
 
    pythia_hadron_->rndm.init(seed_new);
    logg[LPythia].debug("pythia_hadron_ : rndm is initialized with seed ",
                        seed_new);
  }
 
  // clang-format on
 
  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_; }
 
  void clear_final_state() { final_state_.clear(); }
 
  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);
    ParticleTypePtr new_type = ParticleType::try_find(pythia_code);
    if (new_type) {
      ParticleData new_particle(ParticleType::find(pythia_code));
      new_particle.set_4momentum(momentum);
      intermediate_particles.push_back(new_particle);
      return true;
    } else {
      // if the particle does not exist in SMASH the pythia event is rerun
      return false;
    }
  }
 
  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);
 
  int fragment_off_hadron(bool from_forward,
                          bool separate_fragment_baryon,
                          std::array<ThreeVector, 3> &evec_basis,
                          double &ppos_string, double &pneg_string,
                          double &QTrx_string_pos, double &QTrx_string_neg,
                          double &QTry_string_pos, double &QTry_string_neg,
                          Pythia8::FlavContainer &flav_string_pos,
                          Pythia8::FlavContainer &flav_string_neg,
                          std::vector<int> &pdgid_frag,
                          std::vector<FourVector> &momentum_frag);
 
  int get_hadrontype_from_quark(int idq1, int idq2);
 
  int get_resonance_from_quark(int idq1, int idq2, double mass);
 
  bool make_lightcone_final_two(
      bool separate_fragment_hadron,
      double ppos_string, double pneg_string,
      double mTrn_had_forward, double mTrn_had_backward,
      double &ppos_had_forward, double &ppos_had_backward,
      double &pneg_had_forward, double &pneg_had_backward);
 
  static double sample_zLund(double a, double b, double mTrn);
 
  bool remake_kinematics_fragments(
      Pythia8::Event &event_fragments, std::array<ThreeVector, 3> &evec_basis,
      double ppos_string, double pneg_string,
      double QTrx_string, double QTry_string,
      double QTrx_add_pos, double QTry_add_pos,
      double QTrx_add_neg, double QTry_add_neg);
 
  void shift_rapidity_event(
      Pythia8::Event &event, std::array<ThreeVector, 3> &evec_basis,
      double factor_yrapid, double diff_yrapid) {
    Pythia8::Vec4 pvec_string_now = Pythia8::Vec4(0., 0., 0., 0.);
    // loop over all particles in the record
    for (int ipyth = 1; ipyth < event.size(); ipyth++) {
      if (!event[ipyth].isFinal()) {
        continue;
      }
 
      FourVector p_frag = FourVector(
          event[ipyth].e(), event[ipyth].px(),
          event[ipyth].py(), event[ipyth].pz());
      const double E_frag = p_frag.x0();
      const double pl_frag = p_frag.threevec() * evec_basis[0];
      const double ppos_frag = (E_frag + pl_frag) * M_SQRT1_2;
      const double pneg_frag = (E_frag - pl_frag) * M_SQRT1_2;
      const double mTrn_frag = std::sqrt(2. * ppos_frag * pneg_frag);
      // evaluate the old momentum rapidity
      const double y_frag = 0.5 * std::log(ppos_frag / pneg_frag);
 
      // obtain the new momentum rapidity
      const double y_new_frag = factor_yrapid * y_frag + diff_yrapid;
      // compute the new four momentum
      const double E_new_frag = mTrn_frag * std::cosh(y_new_frag);
      const double pl_new_frag = mTrn_frag * std::sinh(y_new_frag);
      ThreeVector mom_new_frag =
          p_frag.threevec() + (pl_new_frag - pl_frag) * evec_basis[0];
      Pythia8::Vec4 pvec_new_frag = set_Vec4(E_new_frag, mom_new_frag);
      event[ipyth].p(pvec_new_frag);
      pvec_string_now += pvec_new_frag;
    }
    event[0].p(pvec_string_now);
    event[0].m(pvec_string_now.mCalc());
  }
 
  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);
 
 
  double getPPosA() { return PPosA_;}
 
  double getPNegA() { return PNegA_;}
 
  double getPPosB() { return PPosB_;}
 
  double getPnegB() { return PNegB_;}
 
  double get_massA() { return massA_;}
 
  double get_massB() { return massB_;}
 
  double get_sqrts() { return sqrtsAB_;}
 
  std::array<PdgCode, 2> get_PDGs() { return PDGcodes_;}
 
  std::array<FourVector, 2> get_plab() { return plab_;}
 
  std::array<FourVector, 2> get_pcom() { return pcom_;}
 
  FourVector get_ucom() { return ucomAB_;}
 
  ThreeVector get_vcom() { return vcomAB_;}
 
  double get_tcoll() { return time_collision_;}
 
  // clang-format on
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_STRINGPROCESS_H_