Version: SMASH-1.5.1
smash::StringProcess Class Reference

#include <processstring.h>

String excitation processes used in SMASH.

Only one instance of this class should be created.

This class implements string excitation processes based on the UrQMD model Bass:1998ca, Bleicher:1999xi and subsequent fragmentation according to the LUND/PYTHIA fragmentation scheme Andersson:1983ia, Sjostrand:2014zea.

The class implements the following functionality:

  • given two colliding initial state particles it provides hadronic final state after single diffractive, double diffractive and non-diffractive string excitation
  • owns a Pythia8::SigmaTotal object to compute corresponding cross-sections
  • owns a Pythia object, that allows to fragment strings

Definition at line 46 of file processstring.h.

Collaboration diagram for smash::StringProcess:
[legend]

Public Member Functions

 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)
 Constructor, initializes pythia. More...
 
void common_setup_pythia (Pythia8::Pythia *pythia_in, double strange_supp, double diquark_supp, double stringz_a, double stringz_b, double string_sigma_T)
 Common setup of PYTHIA objects for soft and hard string routines. More...
 
void init_pythia_hadron_rndm ()
 Set PYTHIA random seeds to be desired values. More...
 
Pythia8::Pythia * get_ptr_pythia_parton ()
 Function to get the PYTHIA object for hard string routine. More...
 
std::array< double, 3 > cross_sections_diffractive (int pdg_a, int pdg_b, double sqrt_s)
 Interface to pythia_sigmatot_ to compute cross-sections of A+B-> different final states Schuler:1993wr. More...
 
void set_pmin_gluon_lightcone (double p_light_cone_min)
 set the minimum lightcone momentum scale carried by gluon. More...
 
void set_pow_fgluon (double betapow)
 lightcone momentum fraction of gluon is sampled according to probability distribution P(x) = 1/x * (1 - x)^{1 + pow_fgluon_beta_} in double-diffractive processes. More...
 
void set_pow_fquark (double alphapow, double betapow)
 lightcone momentum fraction of quark is sampled according to probability distribution \( P(x) = x^{pow_fquark_alpha_ - 1} * (1 - x)^{pow_fquark_beta_ - 1} \) in non-diffractive processes. More...
 
void set_sigma_qperp_ (double sigma_qperp)
 set the average amount of transverse momentum transfer sigma_qperp_. More...
 
void set_tension_string (double kappa_string)
 set the string tension, which is used in append_final_state. More...
 
void init (const ParticleList &incoming, double tcoll)
 initialization feed intial particles, time of collision and gamma factor of the center of mass. More...
 
void compute_incoming_lightcone_momenta ()
 compute the lightcone momenta of incoming particles where the longitudinal direction is set to be same as that of the three-momentum of particle A. More...
 
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)
 Determine string masses and directions in which strings are stretched. More...
 
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)
 Prepare kinematics of two strings, fragment them and append to final_state. More...
 
bool next_SDiff (bool is_AB_to_AX)
 Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns. More...
 
bool next_DDiff ()
 Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process, but lightcone momenta of gluons are sampled in the same was as the UrQMD model Bass:1998ca, Bleicher:1999xi. More...
 
bool next_NDiffSoft ()
 Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig. More...
 
bool next_NDiffHard ()
 Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions. More...
 
bool next_BBbarAnn ()
 Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca, Bleicher:1999xi does, it create two mesonic strings after annihilating one quark-antiquark pair. More...
 
void replace_constituent (Pythia8::Particle &particle, std::array< int, 5 > &excess_constituent)
 Convert a partonic PYTHIA partice into the desired species and update the excess of constituents. More...
 
void find_total_number_constituent (Pythia8::Event &event_intermediate, std::array< int, 5 > &nquark_total, std::array< int, 5 > &nantiq_total)
 Compute how many quarks and antiquarks we have in the system, and update the correspoing arrays with size 5. More...
 
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)
 Take total number of quarks and check if the system has enough constituents that need to be converted into other flavors. More...
 
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)
 Take total number of quarks and check if the system has enough constitents that need to be converted into other flavors. More...
 
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)
 Take the intermediate partonic state from PYTHIA event with mapped hadrons and convert constituents into the desired ones according to the excess of quarks and anti-quarks. More...
 
void compose_string_parton (bool find_forward_string, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
 Identify a set of partons, which are connected to form a color-neutral string, from a given PYTHIA event record. More...
 
void compose_string_junction (bool &find_forward_string, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
 Identify a set of partons and junction(s), which are connected to form a color-neutral string, from a given PYTHIA event record. More...
 
void find_junction_leg (bool sign_color, std::vector< int > &col, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
 Identify partons, which are associated with junction legs, from a given PYTHIA event record. More...
 
int get_index_forward (bool find_forward, int np_end, Pythia8::Event &event)
 Obtain index of the most forward or backward particle in a given PYTHIA event record. More...
 
ParticleList get_final_state ()
 a function to get the final state particle list which is called after the collision More...
 
int append_final_state (ParticleList &intermediate_particles, const FourVector &uString, const ThreeVector &evecLong)
 compute the formation time and fill the arrays with final-state particles as described in Andersson:1983ia. More...
 
int fragment_string (int idq1, int idq2, double mString, ThreeVector &evecLong, bool flip_string_ends, bool separate_fragment_baryon, ParticleList &intermediate_particles)
 perform string fragmentation to determine species and momenta of hadrons by implementing PYTHIA 8.2 Andersson:1983ia, Sjostrand:2014zea. More...
 

Static Public Member Functions

static void make_orthonormal_basis (ThreeVector &evec_polar, std::array< ThreeVector, 3 > &evec_basis)
 compute three orthonormal basis vectors from unit vector in the longitudinal direction More...
 
static void find_excess_constituent (PdgCode &pdg_actual, PdgCode &pdg_mapped, std::array< int, 5 > &excess_quark, std::array< int, 5 > &excess_antiq)
 Compare the valence quark contents of the actual and mapped hadrons and evaluate how many more constituents the actual hadron has compared to the mapped one. More...
 
static bool append_intermediate_list (int pdgid, FourVector momentum, ParticleList &intermediate_particles)
 append new particle from PYTHIA to a specific particle list More...
 
static void convert_KaonLS (int &pythia_id)
 convert Kaon-L or Kaon-S into K0 or Anti-K0 More...
 
static void quarks_from_diquark (int diquark, int &q1, int &q2, int &deg_spin)
 find two quarks from a diquark. More...
 
static int diquark_from_quarks (int q1, int q2)
 Construct diquark from two quarks. More...
 
static void make_string_ends (const PdgCode &pdgcode_in, int &idq1, int &idq2, double xi)
 make a random selection to determine partonic contents at the string ends. More...
 
static Pythia8::Vec4 set_Vec4 (double energy, const ThreeVector &mom)
 Easy setter of Pythia Vec4 from SMASH. More...
 
static FourVector reorient (Pythia8::Particle &particle, std::array< ThreeVector, 3 > &evec_basis)
 compute the four-momentum properly oriented in the lab frame. More...
 
static void assign_all_scaling_factors (int baryon_string, ParticleList &outgoing_particles, const ThreeVector &evecLong, double suppression_factor)
 Assign a cross section scaling factor to all outgoing particles. More...
 
static std::pair< int, int > find_leading (int nq1, int nq2, ParticleList &list)
 Find the leading string fragments. More...
 
static void assign_scaling_factor (int nquark, ParticleData &data, double suppression_factor)
 Assign a cross section scaling factor to the given particle. More...
 
static int pdg_map_for_pythia (PdgCode &pdg)
 Take pdg code and map onto particle specie which can be handled by PYTHIA. More...
 

Private Attributes

double PPosA_
 forward lightcone momentum p^{+} of incoming particle A in CM-frame [GeV] More...
 
double PPosB_
 forward lightcone momentum p^{+} of incoming particle B in CM-frame [GeV] More...
 
double PNegA_
 backward lightcone momentum p^{-} of incoming particle A in CM-frame [GeV] More...
 
double PNegB_
 backward lightcone momentum p^{-} of incoming particle B in CM-frame [GeV] More...
 
double massA_
 mass of incoming particle A [GeV] More...
 
double massB_
 mass of incoming particle B [GeV] More...
 
double sqrtsAB_
 sqrt of Mandelstam variable s of collision [GeV] More...
 
std::array< PdgCode, 2 > PDGcodes_
 PdgCodes of incoming particles. More...
 
std::array< FourVector, 2 > plab_
 momenta of incoming particles in the lab frame [GeV] More...
 
std::array< FourVector, 2 > pcom_
 momenta of incoming particles in the center of mass frame [GeV] More...
 
FourVector ucomAB_
 velocity four vector of the center of mass in the lab frame More...
 
ThreeVector vcomAB_
 velocity three vector of the center of mass in the lab frame More...
 
std::array< ThreeVector, 3 > evecBasisAB_
 Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of incoming particle A. More...
 
int NpartFinal_
 total number of final state particles More...
 
std::array< int, 2 > NpartString_
 number of particles fragmented from strings More...
 
double pmin_gluon_lightcone_
 the minimum lightcone momentum scale carried by a gluon [GeV] More...
 
double pow_fgluon_beta_
 parameter \(\beta\) for the gluon distribution function \( P(x) = x^{-1} (1 - x)^{1 + \beta} \) More...
 
double pow_fquark_alpha_
 parameter \(\alpha\) for the quark distribution function \( P(x) = x^{\alpha - 1} (1 - x)^{\beta - 1} \) More...
 
double pow_fquark_beta_
 parameter \(\beta\) for the quark distribution function \( P(x) = x^{\alpha - 1} (1 - x)^{\beta - 1} \) More...
 
double sigma_qperp_
 Transverse momentum spread of the excited strings. More...
 
double leading_frag_mean_
 mean value of the fragmentation function for the leading baryons More...
 
double leading_frag_width_
 width of the fragmentation function for the leading baryons More...
 
double kappa_tension_string_
 string tension [GeV/fm] More...
 
double additional_xsec_supp_
 additional cross-section suppression factor to take coherence effect into account. More...
 
double time_formation_const_
 constant proper time in the case of constant formation time [fm] More...
 
double soft_t_form_
 factor to be multiplied to formation times in soft strings More...
 
double time_collision_
 time of collision in the computational frame [fm] More...
 
bool use_yoyo_model_
 Whether to calculate the string formation times from the yoyo-model. More...
 
double prob_proton_to_d_uu_
 Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice. More...
 
double sqrt2_
 square root of 2 ( \(\sqrt{2}\)) More...
 
bool pythia_parton_initialized_ = false
 Remembers if Pythia is initialized or not. More...
 
ParticleList final_state_
 final state array which must be accessed after the collision More...
 
std::unique_ptr< Pythia8::Pythia > pythia_parton_
 PYTHIA object used in hard string routine. More...
 
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
 PYTHIA object used in fragmentation. More...
 
Pythia8::SigmaTotal pythia_sigmatot_
 An object to compute cross-sections. More...
 
Pythia8::Event event_intermediate_
 event record for intermediate partonic state in the hard string routine More...
 

Constructor & Destructor Documentation

◆ StringProcess()

smash::StringProcess::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 
)

Constructor, initializes pythia.

Should only be called once.

Parameters
[in]string_tensionvalue of kappa_tension_string_ [GeV/fm]
[in]time_formationvalue of time_formation_const_ [fm]
[in]gluon_betavalue of pow_fgluon_beta_
[in]gluon_pminvalue of pmin_gluon_lightcone_
[in]quark_alphavalue of pow_fquark_alpha_
[in]quark_betavalue of pow_fquark_beta_
[in]strange_suppstrangeness suppression factor (StringFlav:probStoUD) in fragmentation
[in]diquark_suppdiquark suppression factor (StringFlav:probQQtoQ) in fragmentation
[in]sigma_perpvalue of sigma_qperp_ [GeV]
[in]leading_frag_meanMean of the Gaussian used as fragmentation function for leading baryons.
[in]leading_frag_widthWidth of the Gaussian used as fragmentation function for leading baryons.
[in]stringz_aparameter (StringZ:aLund) for the fragmentation function
[in]stringz_bparameter (StringZ:bLund) for the fragmentation function [GeV^-2]
[in]string_sigma_Ttransverse momentum spread (StringPT:sigma) in fragmentation [GeV]
[in]factor_t_formto be multiplied to soft string formation times
[in]use_yoyo_modelWhether to calculate formation times from the yoyo-model.
[in]prob_proton_to_d_uuProbability of a nucleon to be split into the quark it contains once and a diquark another flavour.
See also
StringProcess::common_setup_pythia(Pythia8::Pythia *, double, double, double, double, double)
pythia8235/share/Pythia8/xmldoc/FlavourSelection.xml
pythia8235/share/Pythia8/xmldoc/Fragmentation.xml
pythia8235/share/Pythia8/xmldoc/MasterSwitches.xml
pythia8235/share/Pythia8/xmldoc/MultipartonInteractions.xml

Definition at line 19 of file processstring.cc.

28  : pmin_gluon_lightcone_(gluon_pmin),
29  pow_fgluon_beta_(gluon_beta),
30  pow_fquark_alpha_(quark_alpha),
31  pow_fquark_beta_(quark_beta),
32  sigma_qperp_(sigma_perp),
33  leading_frag_mean_(leading_frag_mean),
34  leading_frag_width_(leading_frag_width),
35  kappa_tension_string_(string_tension),
37  time_formation_const_(time_formation),
38  soft_t_form_(factor_t_form),
39  time_collision_(0.),
40  use_yoyo_model_(use_yoyo_model),
41  prob_proton_to_d_uu_(prob_proton_to_d_uu) {
42  // setup and initialize pythia for hard string process
43  pythia_parton_ = make_unique<Pythia8::Pythia>(PYTHIA_XML_DIR, false);
44  /* select only non-diffractive events
45  * diffractive ones are implemented in a separate routine */
46  pythia_parton_->readString("SoftQCD:nonDiffractive = on");
47  pythia_parton_->readString("MultipartonInteractions:pTmin = 1.5");
48  pythia_parton_->readString("HadronLevel:all = off");
49  common_setup_pythia(pythia_parton_.get(), strange_supp, diquark_supp,
50  stringz_a, stringz_b, string_sigma_T);
51 
52  // setup and initialize pythia for fragmentation
53  pythia_hadron_ = make_unique<Pythia8::Pythia>(PYTHIA_XML_DIR, false);
54  /* turn off all parton-level processes to implement only hadronization */
55  pythia_hadron_->readString("ProcessLevel:all = off");
56  common_setup_pythia(pythia_hadron_.get(), strange_supp, diquark_supp,
57  stringz_a, stringz_b, string_sigma_T);
58 
59  /* initialize PYTHIA */
60  pythia_hadron_->init();
61  pythia_sigmatot_.init(&pythia_hadron_->info, pythia_hadron_->settings,
62  &pythia_hadron_->particleData, &pythia_hadron_->rndm);
63  event_intermediate_.init("intermediate partons",
64  &pythia_hadron_->particleData);
65 
66  sqrt2_ = std::sqrt(2.);
67 
68  for (int imu = 0; imu < 3; imu++) {
69  evecBasisAB_[imu] = ThreeVector(0., 0., 0.);
70  }
71 
72  final_state_.clear();
73 }
double additional_xsec_supp_
additional cross-section suppression factor to take coherence effect into account.
double pmin_gluon_lightcone_
the minimum lightcone momentum scale carried by a gluon [GeV]
Definition: processstring.h:83
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
void common_setup_pythia(Pythia8::Pythia *pythia_in, double strange_supp, double diquark_supp, double stringz_a, double stringz_b, double string_sigma_T)
Common setup of PYTHIA objects for soft and hard string routines.
double prob_proton_to_d_uu_
Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice...
Pythia8::SigmaTotal pythia_sigmatot_
An object to compute cross-sections.
std::unique_ptr< Pythia8::Pythia > pythia_parton_
PYTHIA object used in hard string routine.
Pythia8::Event event_intermediate_
event record for intermediate partonic state in the hard string routine
double kappa_tension_string_
string tension [GeV/fm]
double pow_fquark_alpha_
parameter for the quark distribution function
Definition: processstring.h:93
bool use_yoyo_model_
Whether to calculate the string formation times from the yoyo-model.
ParticleList final_state_
final state array which must be accessed after the collision
double leading_frag_mean_
mean value of the fragmentation function for the leading baryons
double time_collision_
time of collision in the computational frame [fm]
double soft_t_form_
factor to be multiplied to formation times in soft strings
double leading_frag_width_
width of the fragmentation function for the leading baryons
double pow_fgluon_beta_
parameter for the gluon distribution function
Definition: processstring.h:88
double sqrt2_
square root of 2 ( )
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
double sigma_qperp_
Transverse momentum spread of the excited strings.
double time_formation_const_
constant proper time in the case of constant formation time [fm]
double pow_fquark_beta_
parameter for the quark distribution function
Definition: processstring.h:98
Here is the call graph for this function:

Member Function Documentation

◆ common_setup_pythia()

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

Common setup of PYTHIA objects for soft and hard string routines.

Parameters
[out]pythia_inpointer to the PYTHIA object
[in]strange_suppstrangeness suppression factor (StringFlav:probStoUD) in fragmentation
[in]diquark_suppdiquark suppression factor (StringFlav:probQQtoQ) in fragmentation
[in]stringz_aparameter (StringZ:aLund) for the fragmentation function
[in]stringz_bparameter (StringZ:bLund) for the fragmentation function
[in]string_sigma_Ttransverse momentum spread (StringPT:sigma) in fragmentation [GeV]
See also
pythia8235/share/Pythia8/xmldoc/FlavourSelection.xml
pythia8235/share/Pythia8/xmldoc/Fragmentation.xml

Definition at line 75 of file processstring.cc.

79  {
80  // choose parametrization for mass-dependent width
81  pythia_in->readString("ParticleData:modeBreitWigner = 4");
82  /* choose minimum transverse momentum scale
83  * involved in partonic interactions */
84  pythia_in->readString("MultipartonInteractions:pTmin = 1.5");
85  pythia_in->readString("MultipartonInteractions:nSample = 10000");
86  // transverse momentum spread in string fragmentation
87  pythia_in->readString("StringPT:sigma = " + std::to_string(string_sigma_T));
88  // diquark suppression factor in string fragmentation
89  pythia_in->readString("StringFlav:probQQtoQ = " +
90  std::to_string(diquark_supp));
91  // strangeness suppression factor in string fragmentation
92  pythia_in->readString("StringFlav:probStoUD = " +
93  std::to_string(strange_supp));
94  // parameters for the fragmentation function
95  pythia_in->readString("StringZ:aLund = " + std::to_string(stringz_a));
96  pythia_in->readString("StringZ:bLund = " + std::to_string(stringz_b));
97 
98  // manually set the parton distribution function
99  pythia_in->readString("PDF:pSet = 13");
100  pythia_in->readString("PDF:pSetB = 13");
101  pythia_in->readString("PDF:piSet = 1");
102  pythia_in->readString("PDF:piSetB = 1");
103  pythia_in->readString("Beams:idA = 2212");
104  pythia_in->readString("Beams:idB = 2212");
105  pythia_in->readString("Beams:eCM = 10.");
106 
107  // set PYTHIA random seed from outside
108  pythia_in->readString("Random:setSeed = on");
109  // suppress unnecessary output
110  pythia_in->readString("Print:quiet = on");
111  // No resonance decays, since the resonances will be handled by SMASH
112  pythia_in->readString("HadronLevel:Decay = off");
113  // set particle masses and widths in PYTHIA to be same with those in SMASH
114  for (auto &ptype : ParticleType::list_all()) {
115  int pdgid = ptype.pdgcode().get_decimal();
116  double mass_pole = ptype.mass();
117  double width_pole = ptype.width_at_pole();
118  // check if the particle species is in PYTHIA
119  if (pythia_in->particleData.isParticle(pdgid)) {
120  // set mass and width in PYTHIA
121  pythia_in->particleData.m0(pdgid, mass_pole);
122  pythia_in->particleData.mWidth(pdgid, width_pole);
123  } else if (pdgid == 310 || pdgid == 130) {
124  // set mass and width of Kaon-L and Kaon-S
125  pythia_in->particleData.m0(pdgid, kaon_mass);
126  pythia_in->particleData.mWidth(pdgid, 0.);
127  }
128  }
129 
130  // make energy-momentum conservation in PYTHIA more precise
131  pythia_in->readString("Check:epTolErr = 1e-6");
132  pythia_in->readString("Check:epTolWarn = 1e-8");
133 }
static const ParticleTypeList & list_all()
Definition: particletype.cc:55
constexpr double kaon_mass
Kaon mass in GeV.
Definition: constants.h:66
Here is the call graph for this function:
Here is the caller graph for this function:

◆ init_pythia_hadron_rndm()

void smash::StringProcess::init_pythia_hadron_rndm ( )
inline

Set PYTHIA random seeds to be desired values.

The value is recalculated such that it is allowed by PYTHIA.

See also
smash::maximum_rndm_seed_in_pythia

Definition at line 232 of file processstring.h.

232  {
233  const auto &log = logger<LogArea::Pythia>();
234  const int seed_new =
236 
237  pythia_hadron_->rndm.init(seed_new);
238  log.debug("pythia_hadron_ : rndm is initialized with seed ", seed_new);
239  }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
T uniform_int(T min, T max)
Definition: random.h:97
constexpr int maximum_rndm_seed_in_pythia
The maximum value of the random seed used in PYTHIA.
Definition: constants.h:113
Here is the call graph for this function:

◆ get_ptr_pythia_parton()

Pythia8::Pythia* smash::StringProcess::get_ptr_pythia_parton ( )
inline

Function to get the PYTHIA object for hard string routine.

Returns
pointer to the PYTHIA object used in hard string routine

Definition at line 247 of file processstring.h.

247 { return pythia_parton_.get(); }
std::unique_ptr< Pythia8::Pythia > pythia_parton_
PYTHIA object used in hard string routine.

◆ cross_sections_diffractive()

std::array<double, 3> smash::StringProcess::cross_sections_diffractive ( int  pdg_a,
int  pdg_b,
double  sqrt_s 
)
inline

Interface to pythia_sigmatot_ to compute cross-sections of A+B-> different final states Schuler:1993wr.

Parameters
[in]pdg_apdg code of incoming particle A
[in]pdg_bpdg code of incoming particle B
[in]sqrt_scollision energy in the center of mass frame [GeV]
Returns
array with single diffractive cross-sections AB->AX, AB->XB and double diffractive AB->XX.

Definition at line 258 of file processstring.h.

259  {
260  // This threshold magic is following Pythia. Todo(ryu): take care of this.
261  double sqrts_threshold = 2. * (1. + 1.0e-6);
262  /* In the case of mesons, the corresponding vector meson masses
263  * are used to evaluate the energy threshold. */
264  const int pdg_a_mod =
265  (std::abs(pdg_a) > 1000) ? pdg_a : 10 * (std::abs(pdg_a) / 10) + 3;
266  const int pdg_b_mod =
267  (std::abs(pdg_b) > 1000) ? pdg_b : 10 * (std::abs(pdg_b) / 10) + 3;
268  sqrts_threshold += pythia_hadron_->particleData.m0(pdg_a_mod) +
269  pythia_hadron_->particleData.m0(pdg_b_mod);
270  /* Constant cross-section for sub-processes below threshold equal to
271  * cross-section at the threshold. */
272  if (sqrt_s < sqrts_threshold) {
273  sqrt_s = sqrts_threshold;
274  }
275  pythia_sigmatot_.calc(pdg_a, pdg_b, sqrt_s);
276  return {pythia_sigmatot_.sigmaAX(), pythia_sigmatot_.sigmaXB(),
277  pythia_sigmatot_.sigmaXX()};
278  }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
Pythia8::SigmaTotal pythia_sigmatot_
An object to compute cross-sections.
Here is the caller graph for this function:

◆ set_pmin_gluon_lightcone()

void smash::StringProcess::set_pmin_gluon_lightcone ( double  p_light_cone_min)
inline

set the minimum lightcone momentum scale carried by gluon.

Todo:
The following set_ functions are replaced with constructor with arguments.

Must be cleaned up if necessary. This is relevant for the double-diffractive process. The minimum lightcone momentum fraction is set to be pmin_gluon_lightcone_/sqrtsAB.

Parameters
p_light_cone_mina value that we want to use for pmin_gluon_lightcone_.

Definition at line 294 of file processstring.h.

294  {
295  pmin_gluon_lightcone_ = p_light_cone_min;
296  }
double pmin_gluon_lightcone_
the minimum lightcone momentum scale carried by a gluon [GeV]
Definition: processstring.h:83

◆ set_pow_fgluon()

void smash::StringProcess::set_pow_fgluon ( double  betapow)
inline

lightcone momentum fraction of gluon is sampled according to probability distribution P(x) = 1/x * (1 - x)^{1 + pow_fgluon_beta_} in double-diffractive processes.

Parameters
betapowis a value that we want to use for pow_fgluon_beta_.

Definition at line 304 of file processstring.h.

304 { pow_fgluon_beta_ = betapow; }
double pow_fgluon_beta_
parameter for the gluon distribution function
Definition: processstring.h:88

◆ set_pow_fquark()

void smash::StringProcess::set_pow_fquark ( double  alphapow,
double  betapow 
)
inline

lightcone momentum fraction of quark is sampled according to probability distribution \( P(x) = x^{pow_fquark_alpha_ - 1} * (1 - x)^{pow_fquark_beta_ - 1} \) in non-diffractive processes.

Parameters
alphapowis a value that we want to use for pow_fquark_alpha_.
betapowis a value that we want to use for pow_fquark_beta_.

Definition at line 313 of file processstring.h.

313  {
314  pow_fquark_alpha_ = alphapow;
315  pow_fquark_beta_ = betapow;
316  }
double pow_fquark_alpha_
parameter for the quark distribution function
Definition: processstring.h:93
double pow_fquark_beta_
parameter for the quark distribution function
Definition: processstring.h:98

◆ set_sigma_qperp_()

void smash::StringProcess::set_sigma_qperp_ ( double  sigma_qperp)
inline

set the average amount of transverse momentum transfer sigma_qperp_.

Parameters
sigma_qperpis a value that we want to use for sigma_qperp_.

Definition at line 321 of file processstring.h.

321 { sigma_qperp_ = sigma_qperp; }
double sigma_qperp_
Transverse momentum spread of the excited strings.

◆ set_tension_string()

void smash::StringProcess::set_tension_string ( double  kappa_string)
inline

set the string tension, which is used in append_final_state.

Parameters
kappa_stringis a value that we want to use for string tension.

Definition at line 326 of file processstring.h.

326  {
327  kappa_tension_string_ = kappa_string;
328  }
double kappa_tension_string_
string tension [GeV/fm]

◆ init()

void smash::StringProcess::init ( const ParticleList &  incoming,
double  tcoll 
)

initialization feed intial particles, time of collision and gamma factor of the center of mass.

Parameters
[in]incomingis the list of initial state particles.
[in]tcollis time of collision.

Definition at line 221 of file processstring.cc.

221  {
222  PDGcodes_[0] = incoming[0].pdgcode();
223  PDGcodes_[1] = incoming[1].pdgcode();
224  massA_ = incoming[0].effective_mass();
225  massB_ = incoming[1].effective_mass();
226 
227  plab_[0] = incoming[0].momentum();
228  plab_[1] = incoming[1].momentum();
229 
230  // compute sqrts and velocity of the center of mass.
231  sqrtsAB_ = (plab_[0] + plab_[1]).abs();
232  ucomAB_ = (plab_[0] + plab_[1]) / sqrtsAB_;
234 
235  pcom_[0] = plab_[0].LorentzBoost(vcomAB_);
236  pcom_[1] = plab_[1].LorentzBoost(vcomAB_);
237 
238  const double pabscomAB = pCM(sqrtsAB_, massA_, massB_);
239  ThreeVector evec_coll = pcom_[0].threevec() / pabscomAB;
241 
243 
244  time_collision_ = tcoll;
245 }
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
Definition: fourvector.h:310
std::array< FourVector, 2 > plab_
momenta of incoming particles in the lab frame [GeV]
Definition: processstring.h:66
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
std::array< FourVector, 2 > pcom_
momenta of incoming particles in the center of mass frame [GeV]
Definition: processstring.h:68
double sqrtsAB_
sqrt of Mandelstam variable s of collision [GeV]
Definition: processstring.h:62
static void make_orthonormal_basis(ThreeVector &evec_polar, std::array< ThreeVector, 3 > &evec_basis)
compute three orthonormal basis vectors from unit vector in the longitudinal direction ...
void compute_incoming_lightcone_momenta()
compute the lightcone momenta of incoming particles where the longitudinal direction is set to be sam...
ThreeVector vcomAB_
velocity three vector of the center of mass in the lab frame
Definition: processstring.h:72
FourVector ucomAB_
velocity four vector of the center of mass in the lab frame
Definition: processstring.h:70
double time_collision_
time of collision in the computational frame [fm]
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
double massB_
mass of incoming particle B [GeV]
Definition: processstring.h:60
double massA_
mass of incoming particle A [GeV]
Definition: processstring.h:58
Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_orthonormal_basis()

void smash::StringProcess::make_orthonormal_basis ( ThreeVector evec_polar,
std::array< ThreeVector, 3 > &  evec_basis 
)
static

compute three orthonormal basis vectors from unit vector in the longitudinal direction

Parameters
[in]evec_polarunit three-vector in the longitudinal direction
[out]evec_basisorthonormal basis vectors of which evec_basis[0] is in the longitudinal direction while evec_basis[1] and evec_basis[2] span the transverse plane.

Definition at line 1639 of file processstring.cc.

1640  {
1641  if (std::abs(evec_polar.x3()) < (1. - 1.0e-8)) {
1642  double ex, ey, et;
1643  double theta, phi;
1644 
1645  // evec_basis[0] is set to be longitudinal direction
1646  evec_basis[0] = evec_polar;
1647 
1648  theta = std::acos(evec_basis[0].x3());
1649 
1650  ex = evec_basis[0].x1();
1651  ey = evec_basis[0].x2();
1652  et = std::sqrt(ex * ex + ey * ey);
1653  if (ey > 0.) {
1654  phi = std::acos(ex / et);
1655  } else {
1656  phi = -std::acos(ex / et);
1657  }
1658 
1659  /* The transverse plane is spanned
1660  * by evec_basis[1] and evec_basis[2]. */
1661  evec_basis[1].set_x1(cos(theta) * cos(phi));
1662  evec_basis[1].set_x2(cos(theta) * sin(phi));
1663  evec_basis[1].set_x3(-sin(theta));
1664 
1665  evec_basis[2].set_x1(-sin(phi));
1666  evec_basis[2].set_x2(cos(phi));
1667  evec_basis[2].set_x3(0.);
1668  } else {
1669  // if evec_polar is very close to the z axis
1670  if (evec_polar.x3() > 0.) {
1671  evec_basis[1] = ThreeVector(1., 0., 0.);
1672  evec_basis[2] = ThreeVector(0., 1., 0.);
1673  evec_basis[0] = ThreeVector(0., 0., 1.);
1674  } else {
1675  evec_basis[1] = ThreeVector(0., 1., 0.);
1676  evec_basis[2] = ThreeVector(1., 0., 0.);
1677  evec_basis[0] = ThreeVector(0., 0., -1.);
1678  }
1679  }
1680 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ compute_incoming_lightcone_momenta()

void smash::StringProcess::compute_incoming_lightcone_momenta ( )

compute the lightcone momenta of incoming particles where the longitudinal direction is set to be same as that of the three-momentum of particle A.

Definition at line 1682 of file processstring.cc.

1682  {
1683  PPosA_ = (pcom_[0].x0() + evecBasisAB_[0] * pcom_[0].threevec()) / sqrt2_;
1684  PNegA_ = (pcom_[0].x0() - evecBasisAB_[0] * pcom_[0].threevec()) / sqrt2_;
1685  PPosB_ = (pcom_[1].x0() + evecBasisAB_[0] * pcom_[1].threevec()) / sqrt2_;
1686  PNegB_ = (pcom_[1].x0() - evecBasisAB_[0] * pcom_[1].threevec()) / sqrt2_;
1687 }
double PPosA_
forward lightcone momentum p^{+} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:50
double PNegB_
backward lightcone momentum p^{-} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:56
std::array< FourVector, 2 > pcom_
momenta of incoming particles in the center of mass frame [GeV]
Definition: processstring.h:68
double PNegA_
backward lightcone momentum p^{-} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:54
double PPosB_
forward lightcone momentum p^{+} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:52
double sqrt2_
square root of 2 ( )
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
Here is the caller graph for this function:

◆ set_mass_and_direction_2strings()

bool smash::StringProcess::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 
)

Determine string masses and directions in which strings are stretched.

Parameters
[in]quarkspdg ids of string ends
[in]pstr_com4-momenta of strings in the C.o.m. frame [GeV]
[out]m_strmasses of strings [GeV]
[out]evec_strare directions in which strings are stretched.
Returns
whether masses are above the threshold

Definition at line 334 of file processstring.cc.

337  {
338  std::array<bool, 2> found_mass;
339  for (int i = 0; i < 2; i++) {
340  found_mass[i] = false;
341 
342  m_str[i] = pstr_com[i].sqr();
343  m_str[i] = (m_str[i] > 0.) ? std::sqrt(m_str[i]) : 0.;
344  const double threshold = pythia_hadron_->particleData.m0(quarks[i][0]) +
345  pythia_hadron_->particleData.m0(quarks[i][1]);
346  // string mass must be larger than threshold set by PYTHIA.
347  if (m_str[i] > threshold) {
348  found_mass[i] = true;
349  /* Determine direction in which string i is stretched.
350  * This is set to be same with the collision axis
351  * in the center of mass frame.
352  * Initial state partons inside incoming hadrons are
353  * moving along the collision axis,
354  * which is parallel to three momenta of incoming hadrons
355  * in the center of mass frame.
356  * Given that partons are assumed to be massless,
357  * their four momenta are null vectors and parallel to pnull.
358  * If we take unit three-vector of prs,
359  * which is pnull in the rest frame of string,
360  * it would be the direction in which string ends are moving. */
361  const ThreeVector mom = pcom_[i].threevec();
362  const FourVector pnull(mom.abs(), mom);
363  const FourVector prs = pnull.LorentzBoost(pstr_com[i].velocity());
364  evec_str[i] = prs.threevec() / prs.threevec().abs();
365  }
366  }
367 
368  return found_mass[0] && found_mass[1];
369 }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
std::array< FourVector, 2 > pcom_
momenta of incoming particles in the center of mass frame [GeV]
Definition: processstring.h:68
Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_final_state_2strings()

bool smash::StringProcess::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 
)

Prepare kinematics of two strings, fragment them and append to final_state.

Parameters
[in]quarkspdg ids of string ends
[in]pstr_com4-momenta of strings in the C.o.m. frame [GeV]
[in]m_strmasses of strings [GeV]
[out]evec_strare directions in which strings are stretched.
[in]flip_string_endsis whether or not we randomly switch string ends.
[in]separate_fragment_baryonis whether to fragment leading baryons (or anti-baryons) with a separate fragmentation function.
Returns
whether fragmentations and final state creation was successful

Definition at line 371 of file processstring.cc.

376  {
377  const std::array<FourVector, 2> ustr_com = {pstr_com[0] / m_str[0],
378  pstr_com[1] / m_str[1]};
379  for (int i = 0; i < 2; i++) {
380  ParticleList new_intermediate_particles;
381 
382  // determine direction in which string i is stretched.
383  ThreeVector evec = evec_str[i];
384  // perform fragmentation and add particles to final_state.
385  int nfrag = fragment_string(quarks[i][0], quarks[i][1], m_str[i], evec,
386  flip_string_ends, separate_fragment_baryon,
387  new_intermediate_particles);
388  if (nfrag <= 0) {
389  NpartString_[i] = 0;
390  return false;
391  }
392  NpartString_[i] =
393  append_final_state(new_intermediate_particles, ustr_com[i], evec);
394  assert(nfrag == NpartString_[i]);
395  }
396  if ((NpartString_[0] > 0) && (NpartString_[1] > 0)) {
398  return true;
399  }
400  return false;
401 }
std::array< int, 2 > NpartString_
number of particles fragmented from strings
Definition: processstring.h:81
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
int fragment_string(int idq1, int idq2, double mString, ThreeVector &evecLong, bool flip_string_ends, bool separate_fragment_baryon, ParticleList &intermediate_particles)
perform string fragmentation to determine species and momenta of hadrons by implementing PYTHIA 8...
int append_final_state(ParticleList &intermediate_particles, const FourVector &uString, const ThreeVector &evecLong)
compute the formation time and fill the arrays with final-state particles as described in Andersson:1...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ next_SDiff()

bool smash::StringProcess::next_SDiff ( bool  is_AB_to_AX)

Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns.

Parameters
[in]is_AB_to_AXspecifies which hadron to excite into a string. true : A + B -> A + X, false : A + B -> X + B
Returns
whether the process is successfully implemented.

Definition at line 250 of file processstring.cc.

250  {
251  NpartFinal_ = 0;
252  NpartString_[0] = 0;
253  NpartString_[1] = 0;
254  final_state_.clear();
255 
256  double massH = is_AB_to_AX ? massA_ : massB_;
257  double mstrMin = is_AB_to_AX ? massB_ : massA_;
258  double mstrMax = sqrtsAB_ - massH;
259 
260  int idqX1, idqX2;
261  double QTrn, QTrx, QTry;
262  double pabscomHX_sqr, massX;
263 
264  // decompose hadron into quarks
265  make_string_ends(is_AB_to_AX ? PDGcodes_[1] : PDGcodes_[0], idqX1, idqX2,
267  // string mass must be larger than threshold set by PYTHIA.
268  mstrMin = pythia_hadron_->particleData.m0(idqX1) +
269  pythia_hadron_->particleData.m0(idqX2);
270  // this threshold cannot be larger than maximum of allowed string mass.
271  if (mstrMin > mstrMax) {
272  return false;
273  }
274  // sample the transverse momentum transfer
275  QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
276  QTry = random::normal(0., sigma_qperp_ / sqrt2_);
277  QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
278  /* sample the string mass and
279  * evaluate the three-momenta of hadron and string. */
280  massX = random::power(-1.0, mstrMin, mstrMax);
281  pabscomHX_sqr = pCM_sqr(sqrtsAB_, massH, massX);
282  /* magnitude of the three momentum must be larger
283  * than the transverse momentum. */
284  const bool foundPabsX = pabscomHX_sqr > QTrn * QTrn;
285 
286  if (!foundPabsX) {
287  return false;
288  }
289  double sign_direction = is_AB_to_AX ? 1. : -1.;
290  // determine three momentum of the final state hadron
291  const ThreeVector cm_momentum =
292  sign_direction *
293  (evecBasisAB_[0] * std::sqrt(pabscomHX_sqr - QTrn * QTrn) +
294  evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry);
295  const FourVector pstrHcom(std::sqrt(pabscomHX_sqr + massH * massH),
296  cm_momentum);
297  const FourVector pstrXcom(std::sqrt(pabscomHX_sqr + massX * massX),
298  -cm_momentum);
299 
300  const FourVector ustrXcom = pstrXcom / massX;
301  /* determine direction in which the string is stretched.
302  * this is set to be same with the the collision axis
303  * in the center of mass frame. */
304  const ThreeVector threeMomentum =
305  is_AB_to_AX ? pcom_[1].threevec() : pcom_[0].threevec();
306  const FourVector pnull = FourVector(threeMomentum.abs(), threeMomentum);
307  const FourVector prs = pnull.LorentzBoost(ustrXcom.velocity());
308  ThreeVector evec = prs.threevec() / prs.threevec().abs();
309  // perform fragmentation and add particles to final_state.
310  ParticleList new_intermediate_particles;
311  bool separate_fragment_baryon = false;
312  int nfrag =
313  fragment_string(idqX1, idqX2, massX, evec, true, separate_fragment_baryon,
314  new_intermediate_particles);
315  if (nfrag < 1) {
316  NpartString_[0] = 0;
317  return false;
318  }
319  NpartString_[0] =
320  append_final_state(new_intermediate_particles, ustrXcom, evec);
321 
322  NpartString_[1] = 1;
323  PdgCode hadron_code = is_AB_to_AX ? PDGcodes_[0] : PDGcodes_[1];
324  ParticleData new_particle(ParticleType::find(hadron_code));
325  new_particle.set_4momentum(pstrHcom);
326  new_particle.set_cross_section_scaling_factor(1.);
327  new_particle.set_formation_time(time_collision_);
328  final_state_.push_back(new_particle);
329 
331  return true;
332 }
std::array< int, 2 > NpartString_
number of particles fragmented from strings
Definition: processstring.h:81
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
static void make_string_ends(const PdgCode &pdgcode_in, int &idq1, int &idq2, double xi)
make a random selection to determine partonic contents at the string ends.
int fragment_string(int idq1, int idq2, double mString, ThreeVector &evecLong, bool flip_string_ends, bool separate_fragment_baryon, ParticleList &intermediate_particles)
perform string fragmentation to determine species and momenta of hadrons by implementing PYTHIA 8...
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
int append_final_state(ParticleList &intermediate_particles, const FourVector &uString, const ThreeVector &evecLong)
compute the formation time and fill the arrays with final-state particles as described in Andersson:1...
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
double prob_proton_to_d_uu_
Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice...
std::array< FourVector, 2 > pcom_
momenta of incoming particles in the center of mass frame [GeV]
Definition: processstring.h:68
double sqrtsAB_
sqrt of Mandelstam variable s of collision [GeV]
Definition: processstring.h:62
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
ParticleList final_state_
final state array which must be accessed after the collision
double time_collision_
time of collision in the computational frame [fm]
T power(T n, T xMin, T xMax)
Draws a random number according to a power-law distribution ~ x^n.
Definition: random.h:200
double sqrt2_
square root of 2 ( )
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:247
double sigma_qperp_
Transverse momentum spread of the excited strings.
double massB_
mass of incoming particle B [GeV]
Definition: processstring.h:60
double massA_
mass of incoming particle A [GeV]
Definition: processstring.h:58
Here is the call graph for this function:
Here is the caller graph for this function:

◆ next_DDiff()

bool smash::StringProcess::next_DDiff ( )

Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process, but lightcone momenta of gluons are sampled in the same was as the UrQMD model Bass:1998ca, Bleicher:1999xi.

String masses are computed after pomeron exchange aquiring transverse momentum transfer.

Returns
whether the process is successfully implemented.

Definition at line 404 of file processstring.cc.

404  {
405  NpartFinal_ = 0;
406  NpartString_[0] = 0;
407  NpartString_[1] = 0;
408  final_state_.clear();
409 
410  std::array<std::array<int, 2>, 2> quarks;
411  std::array<FourVector, 2> pstr_com;
412  std::array<double, 2> m_str;
413  std::array<ThreeVector, 2> evec_str;
414  ThreeVector threeMomentum;
415 
416  // decompose hadron into quark (and diquark) contents
417  make_string_ends(PDGcodes_[0], quarks[0][0], quarks[0][1],
419  make_string_ends(PDGcodes_[1], quarks[1][0], quarks[1][1],
421  // sample the lightcone momentum fraction carried by gluons
422  const double xmin_gluon_fraction = pmin_gluon_lightcone_ / sqrtsAB_;
423  const double xfracA =
424  random::beta_a0(xmin_gluon_fraction, pow_fgluon_beta_ + 1.);
425  const double xfracB =
426  random::beta_a0(xmin_gluon_fraction, pow_fgluon_beta_ + 1.);
427  // sample the transverse momentum transfer
428  const double QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
429  const double QTry = random::normal(0., sigma_qperp_ / sqrt2_);
430  const double QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
431  // evaluate the lightcone momentum transfer
432  const double QPos = -QTrn * QTrn / (2. * xfracB * PNegB_);
433  const double QNeg = QTrn * QTrn / (2. * xfracA * PPosA_);
434  // compute four-momentum of string 1
435  threeMomentum = evecBasisAB_[0] * (PPosA_ + QPos - PNegA_ - QNeg) / sqrt2_ +
436  evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry;
437  pstr_com[0] =
438  FourVector((PPosA_ + QPos + PNegA_ + QNeg) / sqrt2_, threeMomentum);
439  // compute four-momentum of string 2
440  threeMomentum = evecBasisAB_[0] * (PPosB_ - QPos - PNegB_ + QNeg) / sqrt2_ -
441  evecBasisAB_[1] * QTrx - evecBasisAB_[2] * QTry;
442  pstr_com[1] =
443  FourVector((PPosB_ - QPos + PNegB_ - QNeg) / sqrt2_, threeMomentum);
444 
445  const bool found_masses =
446  set_mass_and_direction_2strings(quarks, pstr_com, m_str, evec_str);
447  if (!found_masses) {
448  return false;
449  }
450  const bool flip_string_ends = true;
451  const bool separate_fragment_baryon = false;
452  const bool success =
453  make_final_state_2strings(quarks, pstr_com, m_str, evec_str,
454  flip_string_ends, separate_fragment_baryon);
455  return success;
456 }
std::array< int, 2 > NpartString_
number of particles fragmented from strings
Definition: processstring.h:81
double PPosA_
forward lightcone momentum p^{+} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:50
T beta_a0(T xmin, T b)
Draws a random number from a beta-distribution with a = 0.
Definition: random.h:348
double pmin_gluon_lightcone_
the minimum lightcone momentum scale carried by a gluon [GeV]
Definition: processstring.h:83
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)
Prepare kinematics of two strings, fragment them and append to final_state.
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
static void make_string_ends(const PdgCode &pdgcode_in, int &idq1, int &idq2, double xi)
make a random selection to determine partonic contents at the string ends.
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
double prob_proton_to_d_uu_
Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice...
double PNegB_
backward lightcone momentum p^{-} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:56
double sqrtsAB_
sqrt of Mandelstam variable s of collision [GeV]
Definition: processstring.h:62
ParticleList final_state_
final state array which must be accessed after the collision
double PNegA_
backward lightcone momentum p^{-} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:54
double PPosB_
forward lightcone momentum p^{+} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:52
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)
Determine string masses and directions in which strings are stretched.
double pow_fgluon_beta_
parameter for the gluon distribution function
Definition: processstring.h:88
double sqrt2_
square root of 2 ( )
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:247
double sigma_qperp_
Transverse momentum spread of the excited strings.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ next_NDiffSoft()

bool smash::StringProcess::next_NDiffSoft ( )

Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig.

This involves a parton exchange in conjunction with momentum transfer. Probability distribution function of the lightcone momentum fraction carried by quark is based on the UrQMD model Bass:1998ca, Bleicher:1999xi.

Returns
whether the process is successfully implemented.
Exceptions
std::runtime_errorif incoming particles are neither mesonic nor baryonic

Definition at line 459 of file processstring.cc.

459  {
460  NpartFinal_ = 0;
461  NpartString_[0] = 0;
462  NpartString_[1] = 0;
463  final_state_.clear();
464 
465  std::array<std::array<int, 2>, 2> quarks;
466  std::array<FourVector, 2> pstr_com;
467  std::array<double, 2> m_str;
468  std::array<ThreeVector, 2> evec_str;
469 
470  // decompose hadron into quark (and diquark) contents
471  int idqA1, idqA2, idqB1, idqB2;
474 
475  const int bar_a = PDGcodes_[0].baryon_number(),
476  bar_b = PDGcodes_[1].baryon_number();
477  if (bar_a == 1 || // baryon-baryon, baryon-meson, baryon-antibaryon
478  (bar_a == 0 && bar_b == 1) || // meson-baryon
479  (bar_a == 0 && bar_b == 0)) { // meson-meson
480  quarks[0][0] = idqB1;
481  quarks[0][1] = idqA2;
482  quarks[1][0] = idqA1;
483  quarks[1][1] = idqB2;
484  } else if (((bar_a == 0) && (bar_b == -1)) || // meson-antibaryon
485  (bar_a == -1)) { // antibaryon-baryon, antibaryon-meson,
486  // antibaryon-antibaryon
487  quarks[0][0] = idqA1;
488  quarks[0][1] = idqB2;
489  quarks[1][0] = idqB1;
490  quarks[1][1] = idqA2;
491  } else {
492  std::stringstream ss;
493  ss << " StringProcess::next_NDiff : baryonA = " << bar_a
494  << ", baryonB = " << bar_b;
495  throw std::runtime_error(ss.str());
496  }
497  // sample the lightcone momentum fraction carried by quarks
498  const double xfracA = random::beta(pow_fquark_alpha_, pow_fquark_beta_);
499  const double xfracB = random::beta(pow_fquark_alpha_, pow_fquark_beta_);
500  // sample the transverse momentum transfer
501  const double QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
502  const double QTry = random::normal(0., sigma_qperp_ / sqrt2_);
503  const double QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
504  // evaluate the lightcone momentum transfer
505  const double QPos = -QTrn * QTrn / (2. * xfracB * PNegB_);
506  const double QNeg = QTrn * QTrn / (2. * xfracA * PPosA_);
507  const double dPPos = -xfracA * PPosA_ - QPos;
508  const double dPNeg = xfracB * PNegB_ - QNeg;
509  // compute four-momentum of string 1
510  ThreeVector threeMomentum =
511  evecBasisAB_[0] * (PPosA_ + dPPos - PNegA_ - dPNeg) / sqrt2_ +
512  evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry;
513  pstr_com[0] =
514  FourVector((PPosA_ + dPPos + PNegA_ + dPNeg) / sqrt2_, threeMomentum);
515  m_str[0] = pstr_com[0].sqr();
516  // compute four-momentum of string 2
517  threeMomentum = evecBasisAB_[0] * (PPosB_ - dPPos - PNegB_ + dPNeg) / sqrt2_ -
518  evecBasisAB_[1] * QTrx - evecBasisAB_[2] * QTry;
519  pstr_com[1] =
520  FourVector((PPosB_ - dPPos + PNegB_ - dPNeg) / sqrt2_, threeMomentum);
521 
522  const bool found_masses =
523  set_mass_and_direction_2strings(quarks, pstr_com, m_str, evec_str);
524  if (!found_masses) {
525  return false;
526  }
527  const bool flip_string_ends = false;
528  const bool separate_fragment_baryon = true;
529  const bool success =
530  make_final_state_2strings(quarks, pstr_com, m_str, evec_str,
531  flip_string_ends, separate_fragment_baryon);
532  return success;
533 }
std::array< int, 2 > NpartString_
number of particles fragmented from strings
Definition: processstring.h:81
double PPosA_
forward lightcone momentum p^{+} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:50
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)
Prepare kinematics of two strings, fragment them and append to final_state.
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
static void make_string_ends(const PdgCode &pdgcode_in, int &idq1, int &idq2, double xi)
make a random selection to determine partonic contents at the string ends.
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
double prob_proton_to_d_uu_
Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice...
double PNegB_
backward lightcone momentum p^{-} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:56
double pow_fquark_alpha_
parameter for the quark distribution function
Definition: processstring.h:93
ParticleList final_state_
final state array which must be accessed after the collision
double PNegA_
backward lightcone momentum p^{-} of incoming particle A in CM-frame [GeV]
Definition: processstring.h:54
double PPosB_
forward lightcone momentum p^{+} of incoming particle B in CM-frame [GeV]
Definition: processstring.h:52
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)
Determine string masses and directions in which strings are stretched.
double sqrt2_
square root of 2 ( )
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:247
double sigma_qperp_
Transverse momentum spread of the excited strings.
T beta(T a, T b)
Draws a random number from a beta-distribution, where probability density of is .
Definition: random.h:326
double pow_fquark_beta_
parameter for the quark distribution function
Definition: processstring.h:98
Here is the call graph for this function:
Here is the caller graph for this function:

◆ next_NDiffHard()

bool smash::StringProcess::next_NDiffHard ( )

Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions.

Returns
whether the process is successfully implemented.

Definition at line 536 of file processstring.cc.

536  {
537  const auto &log = logger<LogArea::Pythia>();
538  NpartFinal_ = 0;
539  final_state_.clear();
540 
541  log.debug("Hard non-diff. with ", PDGcodes_[0], " + ", PDGcodes_[1],
542  " at CM energy [GeV] ", sqrtsAB_);
543 
544  std::array<int, 2> pdg_for_pythia;
545  std::array<std::array<int, 5>, 2> excess_quark;
546  std::array<std::array<int, 5>, 2> excess_antiq;
547  for (int i = 0; i < 2; i++) {
548  for (int j = 0; j < 5; j++) {
549  excess_quark[i][j] = 0;
550  excess_antiq[i][j] = 0;
551  }
552 
553  // get PDG id used in PYTHIA event generation
554  pdg_for_pythia[i] = pdg_map_for_pythia(PDGcodes_[i]);
555  log.debug(" incoming particle ", i, " : ", PDGcodes_[i],
556  " is mapped onto ", pdg_for_pythia[i]);
557 
558  PdgCode pdgcode_for_pythia(std::to_string(pdg_for_pythia[i]));
559  /* evaluate how many more constituents incoming hadron has
560  * compared to the mapped one. */
561  find_excess_constituent(PDGcodes_[i], pdgcode_for_pythia, excess_quark[i],
562  excess_antiq[i]);
563  log.debug(" excess_quark[", i, "] = (", excess_quark[i][0], ", ",
564  excess_quark[i][1], ", ", excess_quark[i][2], ", ",
565  excess_quark[i][3], ", ", excess_quark[i][4], ")");
566  log.debug(" excess_antiq[", i, "] = (", excess_antiq[i][0], ", ",
567  excess_antiq[i][1], ", ", excess_antiq[i][2], ", ",
568  excess_antiq[i][3], ", ", excess_antiq[i][4], ")");
569  }
570 
571  int previous_idA = pythia_parton_->mode("Beams:idA"),
572  previous_idB = pythia_parton_->mode("Beams:idB");
573  double previous_eCM = pythia_parton_->parm("Beams:eCM");
574  // check if the initial state for PYTHIA remains same.
575  bool same_initial_state = previous_idA == pdg_for_pythia[0] &&
576  previous_idB == pdg_for_pythia[1] &&
577  std::abs(previous_eCM - sqrtsAB_) < really_small;
578 
579  /* Perform PYTHIA initialization if it was not previously initialized
580  * or the initial state changed. */
581  if (!pythia_parton_initialized_ || !same_initial_state) {
582  pythia_parton_->settings.mode("Beams:idA", pdg_for_pythia[0]);
583  pythia_parton_->settings.mode("Beams:idB", pdg_for_pythia[1]);
584  pythia_parton_->settings.parm("Beams:eCM", sqrtsAB_);
585 
587  log.debug("Pythia initialized with ", pdg_for_pythia[0], " + ",
588  pdg_for_pythia[1], " at CM energy [GeV] ", sqrtsAB_);
590  throw std::runtime_error("Pythia failed to initialize.");
591  }
592  }
593  /* Set the random seed of the Pythia random Number Generator.
594  * Pythia's random is controlled by SMASH in every single collision.
595  * In this way we ensure that the results are reproducible
596  * for every event if one knows SMASH random seed. */
597  const int seed_new = random::uniform_int(1, maximum_rndm_seed_in_pythia);
598  pythia_parton_->rndm.init(seed_new);
599  log.debug("pythia_parton_ : rndm is initialized with seed ", seed_new);
600 
601  // Short notation for Pythia event
602  Pythia8::Event &event_hadron = pythia_hadron_->event;
603  log.debug("Pythia hard event created");
604  bool final_state_success = pythia_parton_->next();
605  log.debug("Pythia final state computed, success = ", final_state_success);
606  if (!final_state_success) {
607  return false;
608  }
609 
610  ParticleList new_intermediate_particles;
611  ParticleList new_non_hadron_particles;
612 
613  Pythia8::Vec4 pSum = 0.;
614  event_intermediate_.reset();
615  /* Update the partonic intermediate state from PYTHIA output.
616  * Note that hadronization will be performed separately,
617  * after identification of strings and replacement of constituents. */
618  for (int i = 0; i < pythia_parton_->event.size(); i++) {
619  if (pythia_parton_->event[i].isFinal()) {
620  const int pdgid = pythia_parton_->event[i].id();
621  Pythia8::Vec4 pquark = pythia_parton_->event[i].p();
622  const double mass = pythia_parton_->particleData.m0(pdgid);
623 
624  const int status = pythia_parton_->event[i].status();
625  const int color = pythia_parton_->event[i].col();
626  const int anticolor = pythia_parton_->event[i].acol();
627 
628  pSum += pquark;
629  event_intermediate_.append(pdgid, status, color, anticolor, pquark, mass);
630  }
631  }
632  // add junctions to the intermediate state if there is any.
633  event_intermediate_.clearJunctions();
634  for (int i = 0; i < pythia_parton_->event.sizeJunction(); i++) {
635  const int kind = pythia_parton_->event.kindJunction(i);
636  std::array<int, 3> col;
637  for (int j = 0; j < 3; j++) {
638  col[j] = pythia_parton_->event.colJunction(i, j);
639  }
640  event_intermediate_.appendJunction(kind, col[0], col[1], col[2]);
641  }
642  /* The zeroth entry of event record is supposed to have the information
643  * on the whole system. Specify the total momentum and invariant mass. */
644  event_intermediate_[0].p(pSum);
645  event_intermediate_[0].m(pSum.mCalc());
646 
647  /* Replace quark constituents according to the excess of valence quarks
648  * and then rescale momenta of partons by constant factor
649  * to fulfill the energy-momentum conservation. */
650  bool correct_constituents =
651  restore_constituent(event_intermediate_, excess_quark, excess_antiq);
652  if (!correct_constituents) {
653  log.debug("failed to find correct partonic constituents.");
654  return false;
655  }
656 
657  int npart = event_intermediate_.size();
658  int ipart = 0;
659  while (ipart < npart) {
660  const int pdgid = event_intermediate_[ipart].id();
661  if (event_intermediate_[ipart].isFinal() &&
662  !event_intermediate_[ipart].isParton() &&
663  !pythia_parton_->particleData.isOctetHadron(pdgid)) {
664  log.debug("PDG ID from Pythia: ", pdgid);
665  FourVector momentum = reorient(event_intermediate_[ipart], evecBasisAB_);
666  log.debug("4-momentum from Pythia: ", momentum);
667  bool found_ptype =
668  append_intermediate_list(pdgid, momentum, new_non_hadron_particles);
669  if (!found_ptype) {
670  log.warn("PDG ID ", pdgid,
671  " does not exist in ParticleType - start over.");
672  final_state_success = false;
673  }
674  event_intermediate_.remove(ipart, ipart);
675  npart -= 1;
676  } else {
677  ipart += 1;
678  }
679  }
680 
681  bool hadronize_success = false;
682  bool find_forward_string = true;
683  log.debug("Hard non-diff: partonic process gives ",
684  event_intermediate_.size(), " partons.");
685  // identify and fragment strings until there is no parton left.
686  while (event_intermediate_.size() > 1) {
687  // dummy event to initialize the internal variables of PYTHIA.
688  pythia_hadron_->event.reset();
689  if (!pythia_hadron_->next()) {
690  log.debug(" Dummy event in hard string routine failed.");
691  hadronize_success = false;
692  break;
693  }
694 
695  if (event_intermediate_.sizeJunction() > 0) {
696  // identify string from a junction if there is any.
697  compose_string_junction(find_forward_string, event_intermediate_,
698  pythia_hadron_->event);
699  } else {
700  /* identify string from a most forward or backward parton.
701  * if there is no junction. */
702  compose_string_parton(find_forward_string, event_intermediate_,
703  pythia_hadron_->event);
704  }
705 
706  // fragment the (identified) string into hadrons.
707  hadronize_success = pythia_hadron_->forceHadronLevel(false);
708  log.debug("Pythia hadronized, success = ", hadronize_success);
709 
710  new_intermediate_particles.clear();
711  if (hadronize_success) {
712  for (int i = 0; i < event_hadron.size(); i++) {
713  if (event_hadron[i].isFinal()) {
714  int pythia_id = event_hadron[i].id();
715  log.debug("PDG ID from Pythia: ", pythia_id);
716  /* K_short and K_long need to be converted to K0
717  * since SMASH only knows K0 */
718  convert_KaonLS(pythia_id);
719 
720  /* evecBasisAB_[0] is a unit 3-vector in the collision axis,
721  * while evecBasisAB_[1] and evecBasisAB_[2] spans the transverse
722  * plane. Given that PYTHIA assumes z-direction to be the collision
723  * axis, pz from PYTHIA should be the momentum compoment in
724  * evecBasisAB_[0]. px and py are respectively the momentum components
725  * in two transverse directions evecBasisAB_[1] and evecBasisAB_[2].
726  */
727  FourVector momentum = reorient(event_hadron[i], evecBasisAB_);
728  log.debug("4-momentum from Pythia: ", momentum);
729  log.debug("appending the particle ", pythia_id,
730  " to the intermediate particle list.");
731  bool found_ptype = false;
732  if (event_hadron[i].isHadron()) {
733  found_ptype = append_intermediate_list(pythia_id, momentum,
734  new_intermediate_particles);
735  } else {
736  found_ptype = append_intermediate_list(pythia_id, momentum,
737  new_non_hadron_particles);
738  }
739  if (!found_ptype) {
740  log.warn("PDG ID ", pythia_id,
741  " does not exist in ParticleType - start over.");
742  hadronize_success = false;
743  }
744  }
745  }
746  }
747 
748  /* if hadronization is not successful,
749  * reset the event records, return false and then start over. */
750  if (!hadronize_success) {
751  break;
752  }
753 
754  FourVector uString = FourVector(1., 0., 0., 0.);
755  ThreeVector evec = find_forward_string ? evecBasisAB_[0] : -evecBasisAB_[0];
756  int nfrag = append_final_state(new_intermediate_particles, uString, evec);
757  NpartFinal_ += nfrag;
758 
759  find_forward_string = !find_forward_string;
760  }
761 
762  if (hadronize_success) {
763  // add the final state particles, which are not hadron.
764  for (ParticleData data : new_non_hadron_particles) {
765  data.set_cross_section_scaling_factor(1.);
766  data.set_formation_time(time_collision_);
767  final_state_.push_back(data);
768  }
769  } else {
770  final_state_.clear();
771  }
772 
773  return hadronize_success;
774 }
static FourVector reorient(Pythia8::Particle &particle, std::array< ThreeVector, 3 > &evec_basis)
compute the four-momentum properly oriented in the lab frame.
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:34
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
int append_final_state(ParticleList &intermediate_particles, const FourVector &uString, const ThreeVector &evecLong)
compute the formation time and fill the arrays with final-state particles as described in Andersson:1...
static int pdg_map_for_pythia(PdgCode &pdg)
Take pdg code and map onto particle specie which can be handled by PYTHIA.
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
std::unique_ptr< Pythia8::Pythia > pythia_parton_
PYTHIA object used in hard string routine.
double sqrtsAB_
sqrt of Mandelstam variable s of collision [GeV]
Definition: processstring.h:62
static bool append_intermediate_list(int pdgid, FourVector momentum, ParticleList &intermediate_particles)
append new particle from PYTHIA to a specific particle list
T uniform_int(T min, T max)
Definition: random.h:97
Pythia8::Event event_intermediate_
event record for intermediate partonic state in the hard string routine
void compose_string_parton(bool find_forward_string, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
Identify a set of partons, which are connected to form a color-neutral string, from a given PYTHIA ev...
ParticleList final_state_
final state array which must be accessed after the collision
bool pythia_parton_initialized_
Remembers if Pythia is initialized or not.
static void convert_KaonLS(int &pythia_id)
convert Kaon-L or Kaon-S into K0 or Anti-K0
void compose_string_junction(bool &find_forward_string, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
Identify a set of partons and junction(s), which are connected to form a color-neutral string...
double time_collision_
time of collision in the computational frame [fm]
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)
Take the intermediate partonic state from PYTHIA event with mapped hadrons and convert constituents i...
constexpr int maximum_rndm_seed_in_pythia
The maximum value of the random seed used in PYTHIA.
Definition: constants.h:113
std::array< ThreeVector, 3 > evecBasisAB_
Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of i...
Definition: processstring.h:77
static void find_excess_constituent(PdgCode &pdg_actual, PdgCode &pdg_mapped, std::array< int, 5 > &excess_quark, std::array< int, 5 > &excess_antiq)
Compare the valence quark contents of the actual and mapped hadrons and evaluate how many more consti...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ next_BBbarAnn()

bool smash::StringProcess::next_BBbarAnn ( )

Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca, Bleicher:1999xi does, it create two mesonic strings after annihilating one quark-antiquark pair.

Each string has mass equal to half of sqrts.

Returns
whether the process is successfully implemented.
Exceptions
std::invalid_argumentif incoming particles are not baryon-antibaryon pair

Definition at line 1527 of file processstring.cc.

1527  {
1528  const auto &log = logger<LogArea::Pythia>();
1529  const std::array<FourVector, 2> ustrcom = {FourVector(1., 0., 0., 0.),
1530  FourVector(1., 0., 0., 0.)};
1531 
1532  NpartFinal_ = 0;
1533  NpartString_[0] = 0;
1534  NpartString_[1] = 0;
1535  final_state_.clear();
1536 
1537  log.debug("Annihilation occurs between ", PDGcodes_[0], "+", PDGcodes_[1],
1538  " at CM energy [GeV] ", sqrtsAB_);
1539 
1540  // check if the initial state is baryon-antibaryon pair.
1541  PdgCode baryon = PDGcodes_[0], antibaryon = PDGcodes_[1];
1542  if (baryon.baryon_number() == -1) {
1543  std::swap(baryon, antibaryon);
1544  }
1545  if (baryon.baryon_number() != 1 || antibaryon.baryon_number() != -1) {
1546  throw std::invalid_argument("Expected baryon-antibaryon pair.");
1547  }
1548 
1549  // Count how many qqbar combinations are possible for each quark type
1550  constexpr int n_q_types = 5; // u, d, s, c, b
1551  std::vector<int> qcount_bar, qcount_antibar;
1552  std::vector<int> n_combinations;
1553  bool no_combinations = true;
1554  for (int i = 0; i < n_q_types; i++) {
1555  qcount_bar.push_back(baryon.net_quark_number(i + 1));
1556  qcount_antibar.push_back(-antibaryon.net_quark_number(i + 1));
1557  const int n_i = qcount_bar[i] * qcount_antibar[i];
1558  n_combinations.push_back(n_i);
1559  if (n_i > 0) {
1560  no_combinations = false;
1561  }
1562  }
1563 
1564  /* if it is a BBbar pair but there is no qqbar pair to annihilate,
1565  * nothing happens */
1566  if (no_combinations) {
1567  for (int i = 0; i < 2; i++) {
1568  NpartString_[i] = 1;
1569  ParticleData new_particle(ParticleType::find(PDGcodes_[i]));
1570  new_particle.set_4momentum(pcom_[i]);
1571  new_particle.set_cross_section_scaling_factor(1.);
1572  new_particle.set_formation_time(time_collision_);
1573  final_state_.push_back(new_particle);
1574  }
1576  return true;
1577  }
1578 
1579  // Select qqbar pair to annihilate and remove it away
1580  auto discrete_distr = random::discrete_dist<int>(n_combinations);
1581  const int q_annihilate = discrete_distr() + 1;
1582  qcount_bar[q_annihilate - 1]--;
1583  qcount_antibar[q_annihilate - 1]--;
1584 
1585  // Get the remaining quarks and antiquarks
1586  std::vector<int> remaining_quarks, remaining_antiquarks;
1587  for (int i = 0; i < n_q_types; i++) {
1588  for (int j = 0; j < qcount_bar[i]; j++) {
1589  remaining_quarks.push_back(i + 1);
1590  }
1591  for (int j = 0; j < qcount_antibar[i]; j++) {
1592  remaining_antiquarks.push_back(-(i + 1));
1593  }
1594  }
1595  assert(remaining_quarks.size() == 2);
1596  assert(remaining_antiquarks.size() == 2);
1597 
1598  const std::array<double, 2> mstr = {0.5 * sqrtsAB_, 0.5 * sqrtsAB_};
1599 
1600  // randomly select two quark-antiquark pairs
1601  if (random::uniform_int(0, 1) == 0) {
1602  std::swap(remaining_quarks[0], remaining_quarks[1]);
1603  }
1604  if (random::uniform_int(0, 1) == 0) {
1605  std::swap(remaining_antiquarks[0], remaining_antiquarks[1]);
1606  }
1607  // Make sure it satisfies kinematical threshold constraint
1608  bool kin_threshold_satisfied = true;
1609  for (int i = 0; i < 2; i++) {
1610  const double mstr_min =
1611  pythia_hadron_->particleData.m0(remaining_quarks[i]) +
1612  pythia_hadron_->particleData.m0(remaining_antiquarks[i]);
1613  if (mstr_min > mstr[i]) {
1614  kin_threshold_satisfied = false;
1615  }
1616  }
1617  if (!kin_threshold_satisfied) {
1618  return false;
1619  }
1620  // Fragment two strings
1621  for (int i = 0; i < 2; i++) {
1622  ParticleList new_intermediate_particles;
1623 
1624  ThreeVector evec = pcom_[i].threevec() / pcom_[i].threevec().abs();
1625  const int nfrag =
1626  fragment_string(remaining_quarks[i], remaining_antiquarks[i], mstr[i],
1627  evec, true, false, new_intermediate_particles);
1628  if (nfrag <= 0) {
1629  NpartString_[i] = 0;
1630  return false;
1631  }
1632  NpartString_[i] =
1633  append_final_state(new_intermediate_particles, ustrcom[i], evec);
1634  }
1636  return true;
1637 }
std::array< int, 2 > NpartString_
number of particles fragmented from strings
Definition: processstring.h:81
int NpartFinal_
total number of final state particles
Definition: processstring.h:79
int fragment_string(int idq1, int idq2, double mString, ThreeVector &evecLong, bool flip_string_ends, bool separate_fragment_baryon, ParticleList &intermediate_particles)
perform string fragmentation to determine species and momenta of hadrons by implementing PYTHIA 8...
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
int append_final_state(ParticleList &intermediate_particles, const FourVector &uString, const ThreeVector &evecLong)
compute the formation time and fill the arrays with final-state particles as described in Andersson:1...
std::array< PdgCode, 2 > PDGcodes_
PdgCodes of incoming particles.
Definition: processstring.h:64
std::array< FourVector, 2 > pcom_
momenta of incoming particles in the center of mass frame [GeV]
Definition: processstring.h:68
double sqrtsAB_
sqrt of Mandelstam variable s of collision [GeV]
Definition: processstring.h:62
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
T uniform_int(T min, T max)
Definition: random.h:97
ParticleList final_state_
final state array which must be accessed after the collision
double time_collision_
time of collision in the computational frame [fm]
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_excess_constituent()

void smash::StringProcess::find_excess_constituent ( PdgCode pdg_actual,
PdgCode pdg_mapped,
std::array< int, 5 > &  excess_quark,
std::array< int, 5 > &  excess_antiq 
)
static

Compare the valence quark contents of the actual and mapped hadrons and evaluate how many more constituents the actual hadron has compared to the mapped one.

excess_quark[i - 1] is how many more quarks with flavor i (PDG id i) pdg_actual has compared to pdg_mapped. excess_antiq[i - 1] is how many more antiquarks with flavor i (PDG id -i) pdg_actual has compared to pdg_mapped.

Parameters
[in]pdg_actualPDG code of actual incoming particle.
[in]pdg_mappedPDG code of mapped particles used in PYTHIA event generation.
[out]excess_quarkexcess of quarks.
[out]excess_antiqexcess of anti-quarks.

Definition at line 776 of file processstring.cc.

779  {
780  /* decompose PDG id of the actual hadron and mapped one
781  * to get the valence quark constituents */
782  std::array<int, 3> qcontent_actual = pdg_actual.quark_content();
783  std::array<int, 3> qcontent_mapped = pdg_mapped.quark_content();
784 
785  excess_quark = {0, 0, 0, 0, 0};
786  excess_antiq = {0, 0, 0, 0, 0};
787  for (int i = 0; i < 3; i++) {
788  if (qcontent_actual[i] > 0) { // quark content of the actual hadron
789  int j = qcontent_actual[i] - 1;
790  excess_quark[j] += 1;
791  }
792 
793  if (qcontent_mapped[i] > 0) { // quark content of the mapped hadron
794  int j = qcontent_mapped[i] - 1;
795  excess_quark[j] -= 1;
796  }
797 
798  if (qcontent_actual[i] < 0) { // antiquark content of the actual hadron
799  int j = std::abs(qcontent_actual[i]) - 1;
800  excess_antiq[j] += 1;
801  }
802 
803  if (qcontent_mapped[i] < 0) { // antiquark content of the mapped hadron
804  int j = std::abs(qcontent_mapped[i]) - 1;
805  excess_antiq[j] -= 1;
806  }
807  }
808 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ replace_constituent()

void smash::StringProcess::replace_constituent ( Pythia8::Particle &  particle,
std::array< int, 5 > &  excess_constituent 
)

Convert a partonic PYTHIA partice into the desired species and update the excess of constituents.

If the quark flavor i is converted into another flavor j, excess_constituent[i - 1] increases by 1 and excess_constituent[j - 1] decreases by 1. Note that this happens only if excess_constituent[i - 1] < 0 and excess_constituent[j - 1] > 0 (i.e., the incoming hadron has more constituents with flavor j and less constituents with flavor i, compared to the mapped hadron), so they get closer to 0 after the function call.

Parameters
[out]particlePYTHIA particle object to be converted.
[out]excess_constituentexcess in the number of quark constituents. If the particle has positive (negative) quark number, excess of quarks (anti-quarks) should be used.
See also
StringProcess::restore_constituent(Pythia8::Event &, std::array<std::array<int, 5>, 2> &, std::array<std::array<int, 5>, 2> &)

Definition at line 810 of file processstring.cc.

811  {
812  const auto &log = logger<LogArea::Pythia>();
813 
814  // If the particle is neither quark nor diquark, nothing to do.
815  if (!particle.isQuark() && !particle.isDiquark()) {
816  return;
817  }
818 
819  // If there is no excess of constituents, nothing to do.
820  const std::array<int, 5> excess_null = {0, 0, 0, 0, 0};
821  if (excess_constituent == excess_null) {
822  return;
823  }
824 
825  int nq = 0;
826  std::array<int, 2> pdgid = {0, 0};
827  int spin_deg = 0;
828  int pdgid_new = 0;
829  if (particle.isQuark()) {
830  nq = 1;
831  pdgid[0] = particle.id();
832  } else if (particle.isDiquark()) {
833  nq = 2;
834  quarks_from_diquark(particle.id(), pdgid[0], pdgid[1], spin_deg);
835  }
836 
837  for (int iq = 0; iq < nq; iq++) {
838  int jq = std::abs(pdgid[iq]) - 1;
839  int k_select = 0;
840  std::vector<int> k_found;
841  k_found.clear();
842  // check if the constituent needs to be converted.
843  if (excess_constituent[jq] < 0) {
844  for (int k = 0; k < 5; k++) {
845  // check which specie it can be converted into.
846  if (k != jq && excess_constituent[k] > 0) {
847  k_found.push_back(k);
848  }
849  }
850  }
851 
852  // make a random selection of specie and update the excess of constituent.
853  if (k_found.size() > 0) {
854  const int l =
855  random::uniform_int(0, static_cast<int>(k_found.size()) - 1);
856  k_select = k_found[l];
857  /* flavor jq + 1 is converted into k_select + 1
858  * and excess_constituent is updated. */
859  pdgid[iq] = pdgid[iq] > 0 ? k_select + 1 : -(k_select + 1);
860  excess_constituent[jq] += 1;
861  excess_constituent[k_select] -= 1;
862  }
863  }
864 
865  // determine PDG id of the converted parton.
866  if (particle.isQuark()) {
867  pdgid_new = pdgid[0];
868  } else if (particle.isDiquark()) {
869  if (std::abs(pdgid[0]) < std::abs(pdgid[1])) {
870  std::swap(pdgid[0], pdgid[1]);
871  }
872 
873  pdgid_new = std::abs(pdgid[0]) * 1000 + std::abs(pdgid[1]) * 100;
874  if (std::abs(pdgid[0]) == std::abs(pdgid[1])) {
875  pdgid_new += 3;
876  } else {
877  pdgid_new += spin_deg;
878  }
879 
880  if (particle.id() < 0) {
881  pdgid_new *= -1;
882  }
883  }
884  log.debug(" parton id = ", particle.id(), " is converted to ", pdgid_new);
885 
886  // update the constituent mass and energy.
887  Pythia8::Vec4 pquark = particle.p();
888  double mass_new = pythia_hadron_->particleData.m0(pdgid_new);
889  double e_new = std::sqrt(mass_new * mass_new + pquark.pAbs() * pquark.pAbs());
890  // update the particle object.
891  particle.id(pdgid_new);
892  particle.e(e_new);
893  particle.m(mass_new);
894 }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
static void quarks_from_diquark(int diquark, int &q1, int &q2, int &deg_spin)
find two quarks from a diquark.
T uniform_int(T min, T max)
Definition: random.h:97
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_total_number_constituent()

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

Compute how many quarks and antiquarks we have in the system, and update the correspoing arrays with size 5.

Note that elements of the array (0, 1, 2, 3, 4) correspond to d, u, s, c, b flavors.

Parameters
[in]event_intermediatePYTHIA partonic event record which contains output from PYTHIA (hard) event generation.
[out]nquark_totaltotal number of quarks in the system. This is computed based on event_intermediate.
[out]nantiq_totaltotal number of antiquarks in the system. This is computed based on event_intermediate.

Definition at line 896 of file processstring.cc.

898  {
899  for (int iflav = 0; iflav < 5; iflav++) {
900  nquark_total[iflav] = 0;
901  nantiq_total[iflav] = 0;
902  }
903 
904  for (int ip = 1; ip < event_intermediate.size(); ip++) {
905  if (!event_intermediate[ip].isFinal()) {
906  continue;
907  }
908  const int pdgid = event_intermediate[ip].id();
909  if (pdgid > 0) {
910  // quarks
911  for (int iflav = 0; iflav < 5; iflav++) {
912  nquark_total[iflav] +=
913  pythia_hadron_->particleData.nQuarksInCode(pdgid, iflav + 1);
914  }
915  } else {
916  // antiquarks
917  for (int iflav = 0; iflav < 5; iflav++) {
918  nantiq_total[iflav] += pythia_hadron_->particleData.nQuarksInCode(
919  std::abs(pdgid), iflav + 1);
920  }
921  }
922  }
923 }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
Here is the caller graph for this function:

◆ splitting_gluon_qqbar()

bool smash::StringProcess::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 
)

Take total number of quarks and check if the system has enough constituents that need to be converted into other flavors.

If that is not the case, a gluon is splitted into a quark-antiquark pair with desired flavor, so that their flavor can be changed afterwards. For example, if there is no antiquark in the system and we have excess_antiq = (1, -1, 0, 0, 0) (i.e., one ubar has to be converted into dbar), a gluon will be splitted into u-ubar pair.

Parameters
[out]event_intermediatePYTHIA partonic event record to be updated when a gluon happens to split into a qqbar pair.
[out]nquark_totaltotal number of quarks in the system. This is computed based on event_intermediate.
[out]nantiq_totaltotal number of antiquarks in the system. This is computed based on event_intermediate.
[in]sign_constituenttrue (false) if want to check quarks (antiquarks) and their excesses.
[in]excess_constituentexcess in the number of quark constituents. If sign_constituent is true (false), excess of quarks (anti-quarks) should be used.
Returns
false if there are not enough constituents and there is no gluon to split into desired quark-antiquark pair. Otherwise, it gives true.

Definition at line 925 of file processstring.cc.

928  {
929  const auto &log = logger<LogArea::Pythia>();
930 
931  Pythia8::Vec4 pSum = event_intermediate[0].p();
932 
933  /* compute total number of quark and antiquark constituents
934  * in the whole system. */
935  find_total_number_constituent(event_intermediate, nquark_total, nantiq_total);
936 
937  for (int iflav = 0; iflav < 5; iflav++) {
938  /* Find how many constituent will be in the system after
939  * changing the flavors.
940  * Note that nquark_total is number of constituent right after
941  * the pythia event (with mapped incoming hadrons), while the excess
942  * shows how many constituents we have more or less that nquark_total. */
943  int nquark_final =
944  excess_constituent[0][iflav] + excess_constituent[1][iflav];
945  if (sign_constituent) {
946  nquark_final += nquark_total[iflav];
947  } else {
948  nquark_final += nantiq_total[iflav];
949  }
950  /* Therefore, nquark_final should not be negative.
951  * negative nquark_final means that it will not be possible to
952  * find a constituent to change the flavor. */
953  bool enough_quark = nquark_final >= 0;
954  /* If that is the case, a gluon will be splitted into
955  * a quark-antiquark pair with the desired flavor. */
956  if (!enough_quark) {
957  log.debug(" not enough constituents with flavor ", iflav + 1,
958  " : try to split a gluon to qqbar.");
959  for (int ic = 0; ic < std::abs(nquark_final); ic++) {
960  /* Since each incoming hadron has its own count of the excess,
961  * it is necessary to find which one is problematic. */
962  int ih_mod = -1;
963  if (excess_constituent[0][iflav] < 0) {
964  ih_mod = 0;
965  } else {
966  ih_mod = 1;
967  }
968 
969  /* find the most forward or backward gluon
970  * depending on which incoming hadron is found to be an issue. */
971  int iforward = 1;
972  for (int ip = 2; ip < event_intermediate.size(); ip++) {
973  if (!event_intermediate[ip].isFinal() ||
974  !event_intermediate[ip].isGluon()) {
975  continue;
976  }
977 
978  const double y_gluon_current = event_intermediate[ip].y();
979  const double y_gluon_forward = event_intermediate[iforward].y();
980  if ((ih_mod == 0 && y_gluon_current > y_gluon_forward) ||
981  (ih_mod == 1 && y_gluon_current < y_gluon_forward)) {
982  iforward = ip;
983  }
984  }
985 
986  if (!event_intermediate[iforward].isGluon()) {
987  log.debug("There is no gluon to split into qqbar.");
988  return false;
989  }
990 
991  // four momentum of the original gluon
992  Pythia8::Vec4 pgluon = event_intermediate[iforward].p();
993 
994  const int pdgid = iflav + 1;
995  const double mass = pythia_hadron_->particleData.m0(pdgid);
996  const int status = event_intermediate[iforward].status();
997  /* color and anticolor indices.
998  * the color index of gluon goes to the quark, while
999  * the anticolor index goes to the antiquark */
1000  const int col = event_intermediate[iforward].col();
1001  const int acol = event_intermediate[iforward].acol();
1002 
1003  // three momenta of quark and antiquark
1004  std::array<double, 2> px_quark;
1005  std::array<double, 2> py_quark;
1006  std::array<double, 2> pz_quark;
1007  // energies of quark and antiquark
1008  std::array<double, 2> e_quark;
1009  // four momenta of quark and antiquark
1010  std::array<Pythia8::Vec4, 2> pquark;
1011  // transverse momentum scale of string fragmentation
1012  const double sigma_qt_frag = pythia_hadron_->parm("StringPT:sigma");
1013  // sample relative transverse momentum between quark and antiquark
1014  const double qx = random::normal(0., sigma_qt_frag / sqrt2_);
1015  const double qy = random::normal(0., sigma_qt_frag / sqrt2_);
1016  // setup kinematics
1017  for (int isign = 0; isign < 2; isign++) {
1018  /* As the simplest assumption, the three momentum of gluon
1019  * is equally distributed to quark and antiquark.
1020  * Then, they have a relative transverse momentum. */
1021  px_quark[isign] = 0.5 * pgluon.px() + (isign == 0 ? 1. : -1.) * qx;
1022  py_quark[isign] = 0.5 * pgluon.py() + (isign == 0 ? 1. : -1.) * qy;
1023  pz_quark[isign] = 0.5 * pgluon.pz();
1024  e_quark[isign] =
1025  std::sqrt(mass * mass + px_quark[isign] * px_quark[isign] +
1026  py_quark[isign] * py_quark[isign] +
1027  pz_quark[isign] * pz_quark[isign]);
1028  pquark[isign] = Pythia8::Vec4(px_quark[isign], py_quark[isign],
1029  pz_quark[isign], e_quark[isign]);
1030  }
1031 
1032  /* Total energy is not conserved at this point,
1033  * but this will be cured later. */
1034  pSum += pquark[0] + pquark[1] - pgluon;
1035  // add quark and antiquark to the event record
1036  event_intermediate.append(pdgid, status, col, 0, pquark[0], mass);
1037  event_intermediate.append(-pdgid, status, 0, acol, pquark[1], mass);
1038  // then remove the gluon from the record
1039  event_intermediate.remove(iforward, iforward);
1040 
1041  log.debug(" gluon at iforward = ", iforward, " is splitted into ",
1042  pdgid, ",", -pdgid, " qqbar pair.");
1043  /* Increase the total number of quarks and antiquarks by 1,
1044  * as we have extra ones from a gluon. */
1045  nquark_total[iflav] += 1;
1046  nantiq_total[iflav] += 1;
1047  }
1048  }
1049  }
1050 
1051  /* The zeroth entry of event record is supposed to have the information
1052  * on the whole system. Specify the total momentum and invariant mass. */
1053  event_intermediate[0].p(pSum);
1054  event_intermediate[0].m(pSum.mCalc());
1055 
1056  return true;
1057 }
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
void find_total_number_constituent(Pythia8::Event &event_intermediate, std::array< int, 5 > &nquark_total, std::array< int, 5 > &nantiq_total)
Compute how many quarks and antiquarks we have in the system, and update the correspoing arrays with ...
double sqrt2_
square root of 2 ( )
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:247
Here is the call graph for this function:
Here is the caller graph for this function:

◆ rearrange_excess()

void smash::StringProcess::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 
)

Take total number of quarks and check if the system has enough constitents that need to be converted into other flavors.

If that is not the case, excesses of quarks and antiquarks are modified such that the net quark number of each flavor is conserved. For example, if there is no antiquark in the system and we have excess_antiq = (1, -1, 0, 0, 0) (i.e., one ubar has to be converted into dbar), excess_antiq will be changed into (0, 0, 0, 0, 0) and (-1, 1, 0, 0, 0) will be added to excess_quark (i.e., one d quark has to be converted into u quark instead).

Number of quarks is checked if the first argument is the total number of quarks, and the second and third arguments are respectively excesses of quarks and antiquarks. Number of antiquarks is checked if the first argument is the total number of antiquarks, and the second and third arguments are respectively excesses of antiquarks and quarks.

Parameters
[in]nquark_totaltotal number of quarks (antiquarks) in the system.
[out]excess_quarkexcess of quarks (antiquarks) in incoming particles, compared to the mapped ones.
[out]excess_antiqexcess of anti-quarks (quarks) in incoming particles, compared to the mapped ones.

Definition at line 1059 of file processstring.cc.

1062  {
1063  const auto &log = logger<LogArea::Pythia>();
1064 
1065  for (int iflav = 0; iflav < 5; iflav++) {
1066  /* Find how many constituent will be in the system after
1067  * changing the flavors.
1068  * Note that nquark_total is number of constituent right after
1069  * the pythia event (with mapped incoming hadrons), while the excess
1070  * shows how many constituents we have more or less that nquark_total. */
1071  int nquark_final =
1072  nquark_total[iflav] + excess_quark[0][iflav] + excess_quark[1][iflav];
1073  /* Therefore, nquark_final should not be negative.
1074  * negative nquark_final means that it will not be possible to
1075  * find a constituent to change the flavor. */
1076  bool enough_quark = nquark_final >= 0;
1077  // If that is the case, excess of constituents will be modified
1078  if (!enough_quark) {
1079  log.debug(" not enough constituents with flavor ", iflav + 1,
1080  " : try to modify excess of constituents.");
1081  for (int ic = 0; ic < std::abs(nquark_final); ic++) {
1082  /* Since each incoming hadron has its own count of the excess,
1083  * it is necessary to find which one is problematic. */
1084  int ih_mod = -1;
1085  if (excess_quark[0][iflav] < 0) {
1086  ih_mod = 0;
1087  } else {
1088  ih_mod = 1;
1089  }
1090  /* Increase the excess of both quark and antiquark
1091  * with corresponding flavor (iflav + 1) by 1.
1092  * This is for conservation of the net quark number. */
1093  excess_quark[ih_mod][iflav] += 1;
1094  excess_antiq[ih_mod][iflav] += 1;
1095 
1096  /* Since incoming hadrons are mapped onto ones with
1097  * the same baryon number (or quark number),
1098  * summation of the excesses over all flavors should be zero.
1099  * Therefore, we need to find another flavor which has
1100  * a positive excess and subtract by 1. */
1101  for (int jflav = 0; jflav < 5; jflav++) {
1102  // another flavor with positive excess of constituents
1103  if (jflav != iflav && excess_quark[ih_mod][jflav] > 0) {
1104  /* Decrease the excess of both quark and antiquark
1105  * with corresponding flavor (jflav + 1) by 1. */
1106  excess_quark[ih_mod][jflav] -= 1;
1107  excess_antiq[ih_mod][jflav] -= 1;
1108  /* We only need to find one (another) flavor to subtract.
1109  * No more! */
1110  break;
1111  }
1112  }
1113  }
1114  }
1115  }
1116 }
Here is the caller graph for this function:

◆ restore_constituent()

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

Take the intermediate partonic state from PYTHIA event with mapped hadrons and convert constituents into the desired ones according to the excess of quarks and anti-quarks.

Quark (antiquark) flavor is changed and excess of quark (antiquark) is also updated by calling StringProcess::replace_constituent. Beginning with the most forward (or backward) constituent, conversion is done until the total net quark number of each flavor is same with that of incoming hadrons. (i.e., excess_quark minus excess_antiq of incoming hadrons becomes zero.)

However, note that there are some circumstances where this procedure is not directly carried out. For example, a proton-kaon(+) collision mapped onto a proton-pion(+) might be an issue if it involves d + dbar -> g g partonic interaction, given that we anticipate to change dbar to sbar. If such case occurs, we first try to split gluon into quark-antiquark pair with desired flavor. If there are not enough gluons to split, we try to modify the excesses of constituents such that the net quark number is conserved.

Parameters
[out]event_intermediatePYTHIA partonic event record to be updated according to the valence quark contents of incoming hadrons.
[out]excess_quarkexcess of quarks in incoming particles, compared to the mapped ones.
[out]excess_antiqexcess of anti-quarks in incoming particles, compared to the mapped ones.
See also
StringProcess::replace_constituent(Pythia8::Particle &, std::array<int, 5> &)
StringProcess::splitting_gluon_qqbar(Pythia8::Event &, std::array<int, 5> &, std::array<int, 5> &, bool, std::array<std::array<int, 5>, 2> &)
StringProcess::rearrange_excess(std::array<int, 5> &, std::array<std::array<int, 5>, 2> &, std::array<std::array<int, 5>, 2> &)

Definition at line 1118 of file processstring.cc.

1121  {
1122  const auto &log = logger<LogArea::Pythia>();
1123 
1124  Pythia8::Vec4 pSum = event_intermediate[0].p();
1125  const double energy_init = pSum.e();
1126  log.debug(" initial total energy [GeV] : ", energy_init);
1127 
1128  // Total number of quarks and antiquarks, respectively.
1129  std::array<int, 5> nquark_total;
1130  std::array<int, 5> nantiq_total;
1131 
1132  /* Split a gluon into qqbar if we do not have enough constituents
1133  * to be converted in the system. */
1134  bool split_for_quark = splitting_gluon_qqbar(
1135  event_intermediate, nquark_total, nantiq_total, true, excess_quark);
1136  bool split_for_antiq = splitting_gluon_qqbar(
1137  event_intermediate, nquark_total, nantiq_total, false, excess_antiq);
1138 
1139  /* Modify excess_quark and excess_antiq if we do not have enough constituents
1140  * to be converted in the system. */
1141  if (!split_for_quark || !split_for_antiq) {
1142  rearrange_excess(nquark_total, excess_quark, excess_antiq);
1143  rearrange_excess(nantiq_total, excess_antiq, excess_quark);
1144  }
1145 
1146  // Final check if there are enough constituents.
1147  for (int iflav = 0; iflav < 5; iflav++) {
1148  if (nquark_total[iflav] + excess_quark[0][iflav] + excess_quark[1][iflav] <
1149  0) {
1150  log.debug("Not enough quark constituents of flavor ", iflav + 1);
1151  return false;
1152  }
1153 
1154  if (nantiq_total[iflav] + excess_antiq[0][iflav] + excess_antiq[1][iflav] <
1155  0) {
1156  log.debug("Not enough antiquark constituents of flavor ", -(iflav + 1));
1157  return false;
1158  }
1159  }
1160 
1161  for (int ih = 0; ih < 2; ih++) {
1162  log.debug(" initial excess_quark[", ih, "] = (", excess_quark[ih][0], ", ",
1163  excess_quark[ih][1], ", ", excess_quark[ih][2], ", ",
1164  excess_quark[ih][3], ", ", excess_quark[ih][4], ")");
1165  log.debug(" initial excess_antiq[", ih, "] = (", excess_antiq[ih][0], ", ",
1166  excess_antiq[ih][1], ", ", excess_antiq[ih][2], ", ",
1167  excess_antiq[ih][3], ", ", excess_antiq[ih][4], ")");
1168  }
1169 
1170  bool recovered_quarks = false;
1171  while (!recovered_quarks) {
1172  /* Flavor conversion begins with the most forward and backward parton
1173  * respectively for incoming_particles_[0] and incoming_particles_[1]. */
1174  std::array<bool, 2> find_forward = {true, false};
1175  const std::array<int, 5> excess_null = {0, 0, 0, 0, 0};
1176  std::array<int, 5> excess_total = excess_null;
1177 
1178  for (int ih = 0; ih < 2; ih++) { // loop over incoming hadrons
1179  int nfrag = event_intermediate.size();
1180  for (int np_end = 0; np_end < nfrag - 1; np_end++) { // constituent loop
1181  /* select the np_end-th most forward or backward parton and
1182  * change its specie.
1183  * np_end = 0 corresponds to the most forward,
1184  * np_end = 1 corresponds to the second most forward and so on. */
1185  int iforward =
1186  get_index_forward(find_forward[ih], np_end, event_intermediate);
1187  pSum -= event_intermediate[iforward].p();
1188 
1189  if (event_intermediate[iforward].id() > 0) { // quark and diquark
1190  replace_constituent(event_intermediate[iforward], excess_quark[ih]);
1191  log.debug(" excess_quark[", ih, "] = (", excess_quark[ih][0], ", ",
1192  excess_quark[ih][1], ", ", excess_quark[ih][2], ", ",
1193  excess_quark[ih][3], ", ", excess_quark[ih][4], ")");
1194  } else { // antiquark and anti-diquark
1195  replace_constituent(event_intermediate[iforward], excess_antiq[ih]);
1196  log.debug(" excess_antiq[", ih, "] = (", excess_antiq[ih][0], ", ",
1197  excess_antiq[ih][1], ", ", excess_antiq[ih][2], ", ",
1198  excess_antiq[ih][3], ", ", excess_antiq[ih][4], ")");
1199  }
1200 
1201  const int pdgid = event_intermediate[iforward].id();
1202  Pythia8::Vec4 pquark = event_intermediate[iforward].p();
1203  const double mass = pythia_hadron_->particleData.m0(pdgid);
1204 
1205  const int status = event_intermediate[iforward].status();
1206  const int color = event_intermediate[iforward].col();
1207  const int anticolor = event_intermediate[iforward].acol();
1208 
1209  pSum += pquark;
1210  event_intermediate.append(pdgid, status, color, anticolor, pquark,
1211  mass);
1212 
1213  event_intermediate.remove(iforward, iforward);
1214  /* Now the last np_end + 1 entries in event_intermediate
1215  * are np_end + 1 most forward (or backward) partons. */
1216  }
1217 
1218  // Compute the excess of net quark numbers.
1219  for (int j = 0; j < 5; j++) {
1220  excess_total[j] += (excess_quark[ih][j] - excess_antiq[ih][j]);
1221  }
1222  }
1223 
1224  /* If there is no excess of net quark numbers,
1225  * quark content is considered to be correct. */
1226  recovered_quarks = excess_total == excess_null;
1227  }
1228  log.debug(" valence quark contents of hadons are recovered.");
1229 
1230  log.debug(" current total energy [GeV] : ", pSum.e());
1231  /* rescale momenta of all partons by a constant factor
1232  * to conserve the total energy. */
1233  while (true) {
1234  if (std::abs(pSum.e() - energy_init) < really_small * energy_init) {
1235  break;
1236  }
1237 
1238  double energy_current = pSum.e();
1239  double slope = 0.;
1240  for (int i = 1; i < event_intermediate.size(); i++) {
1241  slope += event_intermediate[i].pAbs2() / event_intermediate[i].e();
1242  }
1243 
1244  const double rescale_factor = 1. + (energy_init - energy_current) / slope;
1245  pSum = 0.;
1246  for (int i = 1; i < event_intermediate.size(); i++) {
1247  const double px = rescale_factor * event_intermediate[i].px();
1248  const double py = rescale_factor * event_intermediate[i].py();
1249  const double pz = rescale_factor * event_intermediate[i].pz();
1250  const double pabs = rescale_factor * event_intermediate[i].pAbs();
1251  const double mass = event_intermediate[i].m();
1252 
1253  event_intermediate[i].px(px);
1254  event_intermediate[i].py(py);
1255  event_intermediate[i].pz(pz);
1256  event_intermediate[i].e(std::sqrt(mass * mass + pabs * pabs));
1257  pSum += event_intermediate[i].p();
1258  }
1259  log.debug(" parton momenta are rescaled by factor of ", rescale_factor);
1260  }
1261 
1262  log.debug(" final total energy [GeV] : ", pSum.e());
1263  /* The zeroth entry of event record is supposed to have the information
1264  * on the whole system. Specify the total momentum and invariant mass. */
1265  event_intermediate[0].p(pSum);
1266  event_intermediate[0].m(pSum.mCalc());
1267 
1268  return true;
1269 }
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:34
void replace_constituent(Pythia8::Particle &particle, std::array< int, 5 > &excess_constituent)
Convert a partonic PYTHIA partice into the desired species and update the excess of constituents...
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
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)
Take total number of quarks and check if the system has enough constitents that need to be converted ...
int get_index_forward(bool find_forward, int np_end, Pythia8::Event &event)
Obtain index of the most forward or backward particle in a given PYTHIA event record.
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)
Take total number of quarks and check if the system has enough constituents that need to be converted...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ compose_string_parton()

void smash::StringProcess::compose_string_parton ( bool  find_forward_string,
Pythia8::Event &  event_intermediate,
Pythia8::Event &  event_hadronize 
)

Identify a set of partons, which are connected to form a color-neutral string, from a given PYTHIA event record.

All partons found are moved into a new event record for the further hadronization process. Note that col and acol of Pythia8::Particle contain information on the color flow. This function begins with the most forward (or backward) parton.

For example, quark (col = 1, acol = 0), gluon (col = 2, acol = 1) and antiquark (col = 0, acol = 2) correspond to a \( \bar{q} \, g \, q \) mesonic string. quark (col = 1, acol = 0) and diquark (col = 0, acol = 1) correspond to a \( qq \, q\) baryonic string.

Parameters
[in]find_forward_stringIf it is set to be true (false), it begins with forward (backward) parton.
[out]event_intermediatePYTHIA event record from which a string is identified. All partons found here are removed.
[out]event_hadronizePYTHIA event record to which partons in a string are added.

Definition at line 1271 of file processstring.cc.

1273  {
1274  const auto &log = logger<LogArea::Pythia>();
1275 
1276  Pythia8::Vec4 pSum = 0.;
1277  event_hadronize.reset();
1278 
1279  // select the most forward or backward parton.
1280  int iforward = get_index_forward(find_forward_string, 0, event_intermediate);
1281  log.debug("Hard non-diff: iforward = ", iforward, "(",
1282  event_intermediate[iforward].id(), ")");
1283 
1284  pSum += event_intermediate[iforward].p();
1285  event_hadronize.append(event_intermediate[iforward]);
1286 
1287  int col_to_find = event_intermediate[iforward].acol();
1288  int acol_to_find = event_intermediate[iforward].col();
1289  event_intermediate.remove(iforward, iforward);
1290  log.debug("Hard non-diff: event_intermediate reduces in size to ",
1291  event_intermediate.size());
1292 
1293  // trace color and anti-color indices and find corresponding partons.
1294  while (col_to_find != 0 || acol_to_find != 0) {
1295  log.debug(" col_to_find = ", col_to_find,
1296  ", acol_to_find = ", acol_to_find);
1297 
1298  int ifound = -1;
1299  for (int i = 1; i < event_intermediate.size(); i++) {
1300  const int pdgid = event_intermediate[i].id();
1301  bool found_col =
1302  col_to_find != 0 && col_to_find == event_intermediate[i].col();
1303  bool found_acol =
1304  acol_to_find != 0 && acol_to_find == event_intermediate[i].acol();
1305  if (found_col) {
1306  log.debug(" col_to_find ", col_to_find, " from i ", i, "(", pdgid,
1307  ") found");
1308  }
1309  if (found_acol) {
1310  log.debug(" acol_to_find ", acol_to_find, " from i ", i, "(", pdgid,
1311  ") found");
1312  }
1313 
1314  if (found_col && !found_acol) {
1315  ifound = i;
1316  col_to_find = event_intermediate[i].acol();
1317  break;
1318  } else if (!found_col && found_acol) {
1319  ifound = i;
1320  acol_to_find = event_intermediate[i].col();
1321  break;
1322  } else if (found_col && found_acol) {
1323  ifound = i;
1324  col_to_find = 0;
1325  acol_to_find = 0;
1326  break;
1327  }
1328  }
1329 
1330  if (ifound < 0) {
1331  event_intermediate.list();
1332  event_intermediate.listJunctions();
1333  event_hadronize.list();
1334  event_hadronize.listJunctions();
1335  if (col_to_find != 0) {
1336  log.error("No parton with col = ", col_to_find);
1337  }
1338  if (acol_to_find != 0) {
1339  log.error("No parton with acol = ", acol_to_find);
1340  }
1341  throw std::runtime_error("Hard string could not be identified.");
1342  } else {
1343  pSum += event_intermediate[ifound].p();
1344  // add a parton to the new event record.
1345  event_hadronize.append(event_intermediate[ifound]);
1346  // then remove from the original event record.
1347  event_intermediate.remove(ifound, ifound);
1348  log.debug("Hard non-diff: event_intermediate reduces in size to ",
1349  event_intermediate.size());
1350  }
1351  }
1352 
1353  /* The zeroth entry of event record is supposed to have the information
1354  * on the whole system. Specify the total momentum and invariant mass. */
1355  event_hadronize[0].p(pSum);
1356  event_hadronize[0].m(pSum.mCalc());
1357 }
int get_index_forward(bool find_forward, int np_end, Pythia8::Event &event)
Obtain index of the most forward or backward particle in a given PYTHIA event record.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ compose_string_junction()

void smash::StringProcess::compose_string_junction ( bool &  find_forward_string,
Pythia8::Event &  event_intermediate,
Pythia8::Event &  event_hadronize 
)

Identify a set of partons and junction(s), which are connected to form a color-neutral string, from a given PYTHIA event record.

All partons found are moved into a new event record for the further hadronization process. Junction topology in PYTHIA combines three quarks (antiquarks) to make a color-neutral baryonic (anti-baryonic) configuration. A junction (anti-junction) carries three color (anti-color) indices which are connected with quarks (antiquarks). This function begins with the first junction.

For example, if there is a kind-1 junction with legs (col = 1, 2 and 3), it will first look for three partons with color indices col = 1, 2 and 3 and trace color indices until each leg is `‘closed’' with quark. If there is no quark in the end, there should be an anti-junction and its legs are connected to partons with corresponding anti-colors.

Parameters
[out]find_forward_stringIf it is set to be true (false), it is a string in the forward (backward) direction.
[out]event_intermediatePYTHIA event record from which a string is identified. All partons and junction(s) found here are removed.
[out]event_hadronizePYTHIA event record to which partons in a string are added.
See also
StringProcess::find_junction_leg(bool, std::vector<int> &, Pythia8::Event &, Pythia8::Event &)

Definition at line 1359 of file processstring.cc.

1361  {
1362  const auto &log = logger<LogArea::Pythia>();
1363 
1364  event_hadronize.reset();
1365 
1366  /* Move the first junction to the event record for hadronization
1367  * and specify color or anti-color indices to be found.
1368  * If junction kind is an odd number, it connects three quarks
1369  * to make a color-neutral baryonic configuration.
1370  * Otherwise, it connects three antiquarks
1371  * to make a color-neutral anti-baryonic configuration. */
1372  const int kind = event_intermediate.kindJunction(0);
1373  bool sign_color = kind % 2 == 1;
1374  std::vector<int> col; // color or anti-color indices of the junction legs
1375  for (int j = 0; j < 3; j++) {
1376  col.push_back(event_intermediate.colJunction(0, j));
1377  }
1378  event_hadronize.appendJunction(kind, col[0], col[1], col[2]);
1379  event_intermediate.eraseJunction(0);
1380  log.debug("junction (", col[0], ", ", col[1], ", ", col[2], ") with kind ",
1381  kind, " will be handled.");
1382 
1383  bool found_string = false;
1384  while (!found_string) {
1385  // trace color or anti-color indices and find corresponding partons.
1386  find_junction_leg(sign_color, col, event_intermediate, event_hadronize);
1387  found_string = true;
1388  for (unsigned int j = 0; j < col.size(); j++) {
1389  found_string = found_string && col[j] == 0;
1390  }
1391  if (!found_string) {
1392  /* if there is any leg which is not closed with parton,
1393  * look over junctions and find connected ones. */
1394  log.debug(" still has leg(s) unfinished.");
1395  sign_color = !sign_color;
1396  std::vector<int> junction_to_move;
1397  for (int i = 0; i < event_intermediate.sizeJunction(); i++) {
1398  const int kind_new = event_intermediate.kindJunction(i);
1399  /* If the original junction is associated with positive baryon number,
1400  * it looks for anti-junctions whose legs are connected with
1401  * anti-quarks (anti-colors in general). */
1402  if (sign_color != (kind_new % 2 == 1)) {
1403  continue;
1404  }
1405 
1406  std::array<int, 3> col_new;
1407  for (int k = 0; k < 3; k++) {
1408  col_new[k] = event_intermediate.colJunction(i, k);
1409  }
1410 
1411  int n_legs_connected = 0;
1412  // loop over remaining legs
1413  for (unsigned int j = 0; j < col.size(); j++) {
1414  if (col[j] == 0) {
1415  continue;
1416  }
1417  for (int k = 0; k < 3; k++) {
1418  if (col[j] == col_new[k]) {
1419  n_legs_connected += 1;
1420  col[j] = 0;
1421  col_new[k] = 0;
1422  }
1423  }
1424  }
1425 
1426  // specify which junction is connected to the original one.
1427  if (n_legs_connected > 0) {
1428  for (int k = 0; k < 3; k++) {
1429  if (col_new[k] != 0) {
1430  col.push_back(col_new[k]);
1431  }
1432  }
1433  log.debug(" junction ", i, " (",
1434  event_intermediate.colJunction(i, 0), ", ",
1435  event_intermediate.colJunction(i, 1), ", ",
1436  event_intermediate.colJunction(i, 2), ") with kind ",
1437  kind_new, " will be added.");
1438  junction_to_move.push_back(i);
1439  }
1440  }
1441 
1442  /* If there is any connected junction,
1443  * move it to the event record for hadronization. */
1444  for (unsigned int i = 0; i < junction_to_move.size(); i++) {
1445  unsigned int imove = junction_to_move[i] - i;
1446  const int kind_add = event_intermediate.kindJunction(imove);
1447  std::array<int, 3> col_add;
1448  for (int k = 0; k < 3; k++) {
1449  col_add[k] = event_intermediate.colJunction(imove, k);
1450  }
1451  // add a junction to the new event record.
1452  event_hadronize.appendJunction(kind_add, col_add[0], col_add[1],
1453  col_add[2]);
1454  // then remove from the original event record.
1455  event_intermediate.eraseJunction(imove);
1456  }
1457  }
1458  }
1459 
1460  Pythia8::Vec4 pSum = event_hadronize[0].p();
1461  find_forward_string = pSum.pz() > 0.;
1462 }
void find_junction_leg(bool sign_color, std::vector< int > &col, Pythia8::Event &event_intermediate, Pythia8::Event &event_hadronize)
Identify partons, which are associated with junction legs, from a given PYTHIA event record...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_junction_leg()

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

Identify partons, which are associated with junction legs, from a given PYTHIA event record.

All partons found are moved into a new event record for the further hadronization process.

Parameters
[in]sign_colortrue (false) if the junction is associated with color (anti-color) indices, corresponding baryonic (anti-baryonic) string
[out]colset of color indices that need to be found. The value is set to be zero once the corresponding partons are found.
[out]event_intermediatePYTHIA event record from which a string is identified. All partons and junction(s) found here are removed.
[out]event_hadronizePYTHIA event record to which partons in a string are added.
See also
StringProcess::compose_string_junction(bool &, Pythia8::Event &, Pythia8::Event &)

Definition at line 1464 of file processstring.cc.

1466  {
1467  const auto &log = logger<LogArea::Pythia>();
1468 
1469  Pythia8::Vec4 pSum = event_hadronize[0].p();
1470  for (unsigned int j = 0; j < col.size(); j++) {
1471  if (col[j] == 0) {
1472  continue;
1473  }
1474  bool found_leg = false;
1475  while (!found_leg) {
1476  int ifound = -1;
1477  for (int i = 1; i < event_intermediate.size(); i++) {
1478  const int pdgid = event_intermediate[i].id();
1479  if (sign_color && col[j] == event_intermediate[i].col()) {
1480  log.debug(" col[", j, "] = ", col[j], " from i ", i, "(", pdgid,
1481  ") found");
1482  ifound = i;
1483  col[j] = event_intermediate[i].acol();
1484  break;
1485  } else if (!sign_color && col[j] == event_intermediate[i].acol()) {
1486  log.debug(" acol[", j, "] = ", col[j], " from i ", i, "(", pdgid,
1487  ") found");
1488  ifound = i;
1489  col[j] = event_intermediate[i].col();
1490  break;
1491  }
1492  }
1493 
1494  if (ifound < 0) {
1495  found_leg = true;
1496  if (event_intermediate.sizeJunction() == 0) {
1497  event_intermediate.list();
1498  event_intermediate.listJunctions();
1499  event_hadronize.list();
1500  event_hadronize.listJunctions();
1501  log.error("No parton with col = ", col[j],
1502  " connected with junction leg ", j);
1503  throw std::runtime_error("Hard string could not be identified.");
1504  }
1505  } else {
1506  pSum += event_intermediate[ifound].p();
1507  // add a parton to the new event record.
1508  event_hadronize.append(event_intermediate[ifound]);
1509  // then remove from the original event record.
1510  event_intermediate.remove(ifound, ifound);
1511  log.debug("Hard non-diff: event_intermediate reduces in size to ",
1512  event_intermediate.size());
1513  if (col[j] == 0) {
1514  found_leg = true;
1515  }
1516  }
1517  }
1518  }
1519 
1520  /* The zeroth entry of event record is supposed to have the information
1521  * on the whole system. Specify the total momentum and invariant mass. */
1522  event_hadronize[0].p(pSum);
1523  event_hadronize[0].m(pSum.mCalc());
1524 }
Here is the caller graph for this function:

◆ get_index_forward()

int smash::StringProcess::get_index_forward ( bool  find_forward,
int  np_end,
Pythia8::Event &  event 
)
inline

Obtain index of the most forward or backward particle in a given PYTHIA event record.

Parameters
[in]find_forwardif it looks for the most forward or backward particle.
[in]np_endnumber of the last particle entries to be excluded in lookup. In other words, it finds the most forward (or backward) particle among event[1, ... , event.size() - 1 - np_end].
[in]eventPYTHIA event record which contains particle entries. Note that event[0] is reserved for information on the entire system.
Returns
index of the selected particle, which is used to access the specific particle entry in the event record.

Definition at line 703 of file processstring.h.

704  {
705  int iforward = 1;
706  for (int ip = 2; ip < event.size() - np_end; ip++) {
707  const double y_quark_current = event[ip].y();
708  const double y_quark_forward = event[iforward].y();
709  if ((find_forward && y_quark_current > y_quark_forward) ||
710  (!find_forward && y_quark_current < y_quark_forward)) {
711  iforward = ip;
712  }
713  }
714  return iforward;
715  }
Here is the caller graph for this function:

◆ get_final_state()

ParticleList smash::StringProcess::get_final_state ( )
inline

a function to get the final state particle list which is called after the collision

Returns
ParticleList filled with the final state particles.

Definition at line 722 of file processstring.h.

722 { return final_state_; }
ParticleList final_state_
final state array which must be accessed after the collision
Here is the caller graph for this function:

◆ append_final_state()

int smash::StringProcess::append_final_state ( ParticleList &  intermediate_particles,
const FourVector uString,
const ThreeVector evecLong 
)

compute the formation time and fill the arrays with final-state particles as described in Andersson:1983ia.

Parameters
[out]intermediate_particleslist of fragmented particles to be appended
[in]uStringis velocity four vector of the string.
[in]evecLongis unit 3-vector in which string is stretched.
Returns
number of hadrons fragmented out of string.
Exceptions
std::invalid_argumentif fragmented particle is not hadron
std::invalid_argumentif string is neither mesonic nor baryonic

Definition at line 136 of file processstring.cc.

138  {
139  int nfrag = 0;
140  int bstring = 0;
141  double p_pos_tot = 0.0, p_neg_tot = 0.0;
142 
143  for (ParticleData &data : intermediate_particles) {
144  nfrag += 1;
145  bstring += data.pdgcode().baryon_number();
146 
147  const double pparallel = data.momentum().threevec() * evecLong;
148  p_pos_tot += (data.momentum().x0() + pparallel) / sqrt2_;
149  p_neg_tot += (data.momentum().x0() - pparallel) / sqrt2_;
150  }
151  assert(nfrag > 0);
152 
153  /* compute the cross section scaling factor for leading hadrons
154  * based on the number of valence quarks. */
155  assign_all_scaling_factors(bstring, intermediate_particles, evecLong,
157 
158  std::vector<double> xvertex_pos, xvertex_neg;
159  xvertex_pos.resize(nfrag + 1);
160  xvertex_neg.resize(nfrag + 1);
161  // x^{+} coordinates of the forward end
162  xvertex_pos[0] = p_pos_tot / kappa_tension_string_;
163  // x^{-} coordinates of the backward end
164  xvertex_neg[nfrag] = p_neg_tot / kappa_tension_string_;
165 
166  for (int i = 0; i < nfrag; i++) {
167  const double pparallel =
168  intermediate_particles[i].momentum().threevec() * evecLong;
169 
170  // recursively compute x^{+} coordinates of q-qbar formation vertex
171  xvertex_pos[i + 1] =
172  xvertex_pos[i] -
173  (intermediate_particles[i].momentum().x0() + pparallel) /
175 
176  // recursively compute x^{-} coordinates of q-qbar formation vertex
177  const int ineg = nfrag - i - 1;
178  xvertex_neg[ineg] =
179  xvertex_neg[ineg + 1] -
180  (intermediate_particles[ineg].momentum().x0() - pparallel) /
182  }
183 
184  // Velocity three-vector to perform Lorentz boost.
185  const ThreeVector vstring = uString.velocity();
186 
187  /* compute the formation times of hadrons
188  * from the lightcone coordinates of q-qbar formation vertices. */
189  for (int i = 0; i < nfrag; i++) {
190  // boost 4-momentum into the center of mass frame
191  FourVector momentum =
192  intermediate_particles[i].momentum().LorentzBoost(-vstring);
193  intermediate_particles[i].set_4momentum(momentum);
194 
195  if (use_yoyo_model_) {
196  // set the formation time and position in the rest frame of string
197  double t_prod = (xvertex_pos[i] + xvertex_neg[i + 1]) / sqrt2_;
198  FourVector fragment_position = FourVector(
199  t_prod, evecLong * (xvertex_pos[i] - xvertex_neg[i + 1]) / sqrt2_);
200  /* boost formation vertex into the center of mass frame
201  * and then into the lab frame */
202  fragment_position = fragment_position.LorentzBoost(-vstring);
203  fragment_position = fragment_position.LorentzBoost(-vcomAB_);
204  intermediate_particles[i].set_slow_formation_times(
206  soft_t_form_ * fragment_position.x0() + time_collision_);
207  } else {
208  ThreeVector v_calc = momentum.LorentzBoost(-vcomAB_).velocity();
209  double gamma_factor = 1.0 / std::sqrt(1 - (v_calc).sqr());
210  intermediate_particles[i].set_slow_formation_times(
212  time_formation_const_ * gamma_factor + time_collision_);
213  }
214 
215  final_state_.push_back(intermediate_particles[i]);
216  }
217 
218  return nfrag;
219 }
double additional_xsec_supp_
additional cross-section suppression factor to take coherence effect into account.
double kappa_tension_string_
string tension [GeV/fm]
static void assign_all_scaling_factors(int baryon_string, ParticleList &outgoing_particles, const ThreeVector &evecLong, double suppression_factor)
Assign a cross section scaling factor to all outgoing particles.
bool use_yoyo_model_
Whether to calculate the string formation times from the yoyo-model.
ThreeVector vcomAB_
velocity three vector of the center of mass in the lab frame
Definition: processstring.h:72
ParticleList final_state_
final state array which must be accessed after the collision
double time_collision_
time of collision in the computational frame [fm]
double soft_t_form_
factor to be multiplied to formation times in soft strings
double sqrt2_
square root of 2 ( )
double time_formation_const_
constant proper time in the case of constant formation time [fm]
Here is the call graph for this function:
Here is the caller graph for this function:

◆ append_intermediate_list()

static bool smash::StringProcess::append_intermediate_list ( int  pdgid,
FourVector  momentum,
ParticleList &  intermediate_particles 
)
inlinestatic

append new particle from PYTHIA to a specific particle list

Parameters
[in]pdgidPDG id of particle
[in]momentumfour-momentum of particle
[out]intermediate_particlesparticle list to which the new particle is added.
Returns
whether PDG id exists in ParticleType table.

Definition at line 748 of file processstring.h.

749  {
750  const std::string s = std::to_string(pdgid);
751  PdgCode pythia_code(s);
752  ParticleData new_particle(ParticleType::find(pythia_code));
753  new_particle.set_4momentum(momentum);
754  intermediate_particles.push_back(new_particle);
755  return true;
756  }
static const ParticleType & find(PdgCode pdgcode)
Returns the ParticleType object for the given pdgcode.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ convert_KaonLS()

static void smash::StringProcess::convert_KaonLS ( int &  pythia_id)
inlinestatic

convert Kaon-L or Kaon-S into K0 or Anti-K0

Parameters
[out]pythia_idis PDG id to be converted.

Definition at line 762 of file processstring.h.

762  {
763  if (pythia_id == 310 || pythia_id == 130) {
764  pythia_id = (random::uniform_int(0, 1) == 0) ? 311 : -311;
765  }
766  }
T uniform_int(T min, T max)
Definition: random.h:97
Here is the call graph for this function:
Here is the caller graph for this function:

◆ quarks_from_diquark()

void smash::StringProcess::quarks_from_diquark ( int  diquark,
int &  q1,
int &  q2,
int &  deg_spin 
)
static

find two quarks from a diquark.

Order does not matter.

Parameters
[in]diquarkPDG id of diquark
[out]q1PDG id of quark 1
[out]q2PDG id of quark 2
[out]deg_spinspin degeneracy

Definition at line 1689 of file processstring.cc.

1690  {
1691  // The 4-digit pdg id should be diquark.
1692  assert((std::abs(diquark) > 1000) && (std::abs(diquark) < 5510) &&
1693  (std::abs(diquark) % 100 < 10));
1694 
1695  // The fourth digit corresponds to the spin degeneracy.
1696  deg_spin = std::abs(diquark) % 10;
1697  // Diquark (anti-diquark) is decomposed into two quarks (antiquarks).
1698  const int sign_anti = diquark > 0 ? 1 : -1;
1699 
1700  // Obtain two quarks (or antiquarks) from the first and second digit.
1701  q1 = sign_anti * (std::abs(diquark) - (std::abs(diquark) % 1000)) / 1000;
1702  q2 = sign_anti * (std::abs(diquark) % 1000 - deg_spin) / 100;
1703 }
Here is the caller graph for this function:

◆ diquark_from_quarks()

int smash::StringProcess::diquark_from_quarks ( int  q1,
int  q2 
)
static

Construct diquark from two quarks.

Order does not matter.

Parameters
[in]q1PDG code of quark 1
[in]q2PDG code of quark 2
Returns
PDG code of diquark composed of q1 and q2

Definition at line 1705 of file processstring.cc.

1705  {
1706  assert((q1 > 0 && q2 > 0) || (q1 < 0 && q2 < 0));
1707  if (std::abs(q1) < std::abs(q2)) {
1708  std::swap(q1, q2);
1709  }
1710  int diquark = std::abs(q1 * 1000 + q2 * 100);
1711  /* Adding spin degeneracy = 2S+1. For identical quarks spin cannot be 0
1712  * because of Pauli exclusion principle, so spin 1 is assumed. Otherwise
1713  * S = 0 with probability 1/4 and S = 1 with probability 3/4. */
1714  diquark += (q1 != q2 && random::uniform_int(0, 3) == 0) ? 1 : 3;
1715  return (q1 < 0) ? -diquark : diquark;
1716 }
T uniform_int(T min, T max)
Definition: random.h:97
Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_string_ends()

void smash::StringProcess::make_string_ends ( const PdgCode pdgcode_in,
int &  idq1,
int &  idq2,
double  xi 
)
static

make a random selection to determine partonic contents at the string ends.

Parameters
[in]pdgcode_inis PdgCode of hadron which transforms into a string.
[out]idq1is PDG id of quark or anti-diquark.
[out]idq2is PDG id of anti-quark or diquark.
[in]xiprobability to split a nucleon into the quark it has only once and a diquark of another flavour.

Definition at line 1718 of file processstring.cc.

1719  {
1720  std::array<int, 3> quarks = pdg.quark_content();
1721  if (pdg.is_nucleon()) {
1722  // protons and neutrons treated seperately since single quarks is at a
1723  // different position in the PDG code
1724  if (pdg.charge() == 0) { // (anti)neutron
1725  if (random::uniform(0., 1.) < xi) {
1726  idq1 = quarks[0];
1727  idq2 = diquark_from_quarks(quarks[1], quarks[2]);
1728  } else {
1729  idq1 = quarks[1];
1730  idq2 = diquark_from_quarks(quarks[0], quarks[2]);
1731  }
1732  } else { // (anti)proton
1733  if (random::uniform(0., 1.) < xi) {
1734  idq1 = quarks[2];
1735  idq2 = diquark_from_quarks(quarks[0], quarks[1]);
1736  } else {
1737  idq1 = quarks[0];
1738  idq2 = diquark_from_quarks(quarks[1], quarks[2]);
1739  }
1740  }
1741  } else {
1742  if (pdg.is_meson()) {
1743  idq1 = quarks[1];
1744  idq2 = quarks[2];
1745  /* Some mesons with PDG id 11X are actually mixed state of uubar and
1746  * ddbar. have a random selection whether we have uubar or ddbar in this
1747  * case. */
1748  if (idq1 == 1 && idq2 == -1 && random::uniform_int(0, 1) == 0) {
1749  idq1 = 2;
1750  idq2 = -2;
1751  }
1752  } else {
1753  assert(pdg.is_baryon());
1754  // Get random quark to position 0
1755  std::swap(quarks[random::uniform_int(0, 2)], quarks[0]);
1756  idq1 = quarks[0];
1757  idq2 = diquark_from_quarks(quarks[1], quarks[2]);
1758  }
1759  }
1760  // Fulfil the convention: idq1 should be quark or anti-diquark
1761  if (idq1 < 0) {
1762  std::swap(idq1, idq2);
1763  }
1764 }
static int diquark_from_quarks(int q1, int q2)
Construct diquark from two quarks.
T uniform_int(T min, T max)
Definition: random.h:97
T uniform(T min, T max)
Definition: random.h:85
Here is the call graph for this function:
Here is the caller graph for this function:

◆ set_Vec4()

static Pythia8::Vec4 smash::StringProcess::set_Vec4 ( double  energy,
const ThreeVector mom 
)
inlinestatic

Easy setter of Pythia Vec4 from SMASH.

Parameters
[in]energytime component
[in]momspatial three-vector
Returns
Pythia Vec4 from energy and ThreeVector

Definition at line 800 of file processstring.h.

800  {
801  return Pythia8::Vec4(mom.x1(), mom.x2(), mom.x3(), energy);
802  }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ reorient()

static FourVector smash::StringProcess::reorient ( Pythia8::Particle &  particle,
std::array< ThreeVector, 3 > &  evec_basis 
)
inlinestatic

compute the four-momentum properly oriented in the lab frame.

While PYTHIA assumes that the collision axis is in z-direction, this is not necessarly the case in SMASH.

Parameters
[in]particleparticle object from PYTHIA event generation where z-axis is set to be the collision axis
[in]evec_basisthree basis vectors in the lab frame evec_basis[0] is unit vector in the collision axis and other two span the transverse plane
Returns
four-momentum of particle in the lab frame

Definition at line 815 of file processstring.h.

816  {
817  ThreeVector three_momentum = evec_basis[0] * particle.pz() +
818  evec_basis[1] * particle.px() +
819  evec_basis[2] * particle.py();
820  return FourVector(particle.e(), three_momentum);
821  }
Here is the caller graph for this function:

◆ fragment_string()

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

perform string fragmentation to determine species and momenta of hadrons by implementing PYTHIA 8.2 Andersson:1983ia, Sjostrand:2014zea.

Parameters
[in]idq1PDG id of quark or anti-diquark (carrying color index).
[in]idq2PDG id of diquark or anti-quark (carrying anti-color index).
[in]mStringthe string mass. [GeV]
[out]evecLongunit 3-vector specifying the direction of diquark or anti-diquark.
[in]flip_string_endsis whether or not we randomly switch string ends.
[in]separate_fragment_baryonis whether fragment leading baryon (or anti-baryon) with separate fragmentation function.
[out]intermediate_particleslist of fragmented particles
Returns
number of hadrons fragmented out of string.
Exceptions
std::runtime_errorif string mass is lower than threshold set by PYTHIA

Definition at line 1766 of file processstring.cc.

1769  {
1770  const auto &log = logger<LogArea::Pythia>();
1771  pythia_hadron_->event.reset();
1772  intermediate_particles.clear();
1773 
1774  log.debug("initial quark content for fragment_string : ", idq1, ", ", idq2);
1775  // PDG id of quark constituents of string ends
1776  std::array<int, 2> idqIn;
1777  idqIn[0] = idq1;
1778  idqIn[1] = idq2;
1779 
1780  int bstring = 0;
1781  // constituent masses of the string
1782  std::array<double, 2> m_const;
1783 
1784  for (int i = 0; i < 2; i++) {
1785  // evaluate total baryon number of the string times 3
1786  bstring += pythia_hadron_->particleData.baryonNumberType(idqIn[i]);
1787 
1788  m_const[i] = pythia_hadron_->particleData.m0(idqIn[i]);
1789  }
1790  log.debug("baryon number of string times 3 : ", bstring);
1791 
1792  if (flip_string_ends && random::uniform_int(0, 1) == 0) {
1793  /* in the case where we flip the string ends,
1794  * we need to flip the longitudinal unit vector itself
1795  * since it is set to be direction of diquark (anti-quark)
1796  * or anti-diquark. */
1797  evecLong = -evecLong;
1798  }
1799 
1800  if (m_const[0] + m_const[1] > mString) {
1801  throw std::runtime_error("String fragmentation: m1 + m2 > mString");
1802  }
1803  Pythia8::Vec4 pSum = 0.;
1804 
1805  int number_of_fragments = 0;
1806 
1807  if (separate_fragment_baryon && (std::abs(bstring) == 3) &&
1808  (mString > m_const[0] + m_const[1] + 1.)) {
1809  /* In the case of separate fragmentation function for leading baryons,
1810  * the leading baryon (antibaryon) is fragmented from the original string
1811  * with a gaussian fragmentation function.
1812  * It is then followed by fragmentation of the remaining mesonic string. */
1813  const double ssbar_supp = pythia_hadron_->parm("StringFlav:probStoUD");
1814  const double sigma_qt_frag = pythia_hadron_->parm("StringPT:sigma");
1815  std::array<ThreeVector, 3> evec_basis;
1816  make_orthonormal_basis(evecLong, evec_basis);
1817 
1818  /* transverse momentum (and magnitude) acquired by the qqbar pair
1819  * created in the middle of string */
1820  double QTrx, QTry, QTrn;
1821  /* lightcone momenta p^+ and p^- of the remaining (mesonic) string
1822  * p^{\pm} is defined as (E \pm p_{longitudinal}) / sqrts{2}. */
1823  double ppos_string_new, pneg_string_new;
1824  // transverse mass of the remaining (mesonic) string
1825  double mTrn_string_new;
1826  // transverse masses of quark constituents of string ends
1827  std::array<double, 2> m_trans;
1828 
1829  const int niter_max = 10000;
1830  int iiter = 0;
1831  bool found_leading_baryon = false;
1832  while (!found_leading_baryon) {
1833  int idnew_qqbar = 0;
1834  /* Determine flavor of the new qqbar pair
1835  * based on the strangeness suppression factor specifed in config. */
1836  if (random::uniform(0., 1. + ssbar_supp) < 1.) {
1837  if (random::uniform_int(0, 1) == 0) {
1838  idnew_qqbar = 1; // ddbar pair
1839  } else {
1840  idnew_qqbar = 2; // uubar pair
1841  }
1842  } else {
1843  idnew_qqbar = 3; // ssbar pair
1844  }
1845 
1846  int id_diquark = 0;
1847  if (bstring > 0) {
1848  /* q(idq1) + diquark(idq2) baryonic string fragments into
1849  * q(idq1) + qbar(-idnew_qqbar) mesonic string and
1850  * q(idnew_qqbar) + diquark(idq2) baryon. */
1851  id_diquark = std::abs(idq2);
1852  idqIn[1] = -idnew_qqbar;
1853  } else {
1854  /* anti-diquark(idq1) + qbar(idq2) anti-baryonic string fragments into
1855  * q(idnew_qqbar) + qbar(idq2) mesonic string and
1856  * anti-diquark(idq1) + qbar(-idnew_qqbar) antibaryon. */
1857  id_diquark = std::abs(idq1);
1858  idqIn[0] = idnew_qqbar;
1859  }
1860  log.debug("quark constituents for leading baryon : ", idnew_qqbar, ", ",
1861  id_diquark);
1862 
1863  // net quark number of d, u, s, c and b flavors
1864  std::array<int, 5> frag_net_q;
1865  /* Evaluate total net quark number of baryon (antibaryon)
1866  * from the valence quark constituents. */
1867  for (int iq = 0; iq < 5; iq++) {
1868  frag_net_q[iq] =
1869  (bstring > 0 ? 1 : -1) *
1870  (pythia_hadron_->particleData.nQuarksInCode(idnew_qqbar, iq + 1) +
1871  pythia_hadron_->particleData.nQuarksInCode(id_diquark, iq + 1));
1872  }
1873  const int frag_iso3 = frag_net_q[1] - frag_net_q[0];
1874  const int frag_strange = -frag_net_q[2];
1875  const int frag_charm = frag_net_q[3];
1876  const int frag_bottom = -frag_net_q[4];
1877  log.debug(" conserved charges of leading baryon : iso3 = ", frag_iso3,
1878  ", strangeness = ", frag_strange, ", charmness = ", frag_charm,
1879  ", bottomness = ", frag_bottom);
1880 
1881  std::vector<int> pdgid_possible;
1882  std::vector<double> weight_possible;
1883  std::vector<double> weight_summed;
1884  /* loop over hadronic species
1885  * Any hadron with the same valence quark contents is allowed and
1886  * the probability goes like spin degeneracy over mass. */
1887  for (auto &ptype : ParticleType::list_all()) {
1888  const int pdgid =
1889  (bstring > 0 ? 1 : -1) * std::abs(ptype.pdgcode().get_decimal());
1890  if ((pythia_hadron_->particleData.isParticle(pdgid)) &&
1891  (bstring == 3 * ptype.pdgcode().baryon_number()) &&
1892  (frag_iso3 == ptype.pdgcode().isospin3()) &&
1893  (frag_strange == ptype.pdgcode().strangeness()) &&
1894  (frag_charm == ptype.pdgcode().charmness()) &&
1895  (frag_bottom == ptype.pdgcode().bottomness())) {
1896  const int spin_degeneracy = ptype.pdgcode().spin_degeneracy();
1897  const double mass_pole = ptype.mass();
1898  const double weight =
1899  static_cast<double>(spin_degeneracy) / mass_pole;
1900  pdgid_possible.push_back(pdgid);
1901  weight_possible.push_back(weight);
1902 
1903  log.debug(" PDG id ", pdgid, " is possible with weight ", weight);
1904  }
1905  }
1906  const int n_possible = pdgid_possible.size();
1907  weight_summed.push_back(0.);
1908  for (int i = 0; i < n_possible; i++) {
1909  weight_summed.push_back(weight_summed[i] + weight_possible[i]);
1910  }
1911 
1912  /* Sample baryon (antibaryon) specie,
1913  * which is fragmented from the leading diquark (anti-diquark). */
1914  int pdgid_frag = 0;
1915  const double uspc = random::uniform(0., weight_summed[n_possible]);
1916  for (int i = 0; i < n_possible; i++) {
1917  if ((uspc >= weight_summed[i]) && (uspc < weight_summed[i + 1])) {
1918  pdgid_frag = pdgid_possible[i];
1919  log.debug("PDG id ", pdgid_frag, " selected for leading baryon.");
1920  break;
1921  }
1922  }
1923  // Sample mass.
1924  const double mass_frag = pythia_hadron_->particleData.mSel(pdgid_frag);
1925 
1926  // Sample transverse momentum carried by baryon (antibaryon).
1927  QTrx = random::normal(0., sigma_qt_frag / sqrt2_);
1928  QTry = random::normal(0., sigma_qt_frag / sqrt2_);
1929  log.debug("Transverse momentum (", QTrx, ", ", QTry,
1930  ") GeV selected for the leading baryon.");
1931  QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
1932  // transverse mass of the fragmented leading hadron
1933  const double mTrn_frag = std::sqrt(QTrn * QTrn + mass_frag * mass_frag);
1934 
1935  /* Sample lightcone momentum fraction carried by the baryon (antibaryon).
1936  * The corresponding fragmentation function has gaussian shape
1937  * i.e., exp( - (x - leading_frag_mean_)^2 / (2 leading_frag_width_^2) ).
1938  * This is done using the rejection method with an Lorentzian
1939  * envelop function, which is
1940  * 1 / (1 + (x - leading_frag_mean_)^2 / (4 leading_frag_width_^2) ). */
1941  double xfrac = 0.;
1942  bool xfrac_accepted = false;
1943  while (!xfrac_accepted) {
1944  const double angle =
1945  random::uniform(0., 1.) *
1946  (std::atan(0.5 * (1. - leading_frag_mean_) /
1948  std::atan(0.5 * leading_frag_mean_ / leading_frag_width_)) -
1949  std::atan(0.5 * leading_frag_mean_ / leading_frag_width_);
1950  // lightcone momentum fraction sampled from the envelop function
1951  xfrac = leading_frag_mean_ + 2. * leading_frag_width_ * std::tan(angle);
1952 
1953  /* The rejection method is applied to sample
1954  * according to the real fragmentation function. */
1955  const double xf_tmp =
1956  std::abs(xfrac - leading_frag_mean_) / leading_frag_width_;
1957  const double xf_env = (1. + really_small) / (1. + xf_tmp * xf_tmp / 4.);
1958  const double xf_pdf = std::exp(-xf_tmp * xf_tmp / 2.);
1959  if (random::uniform(0., xf_env) < xf_pdf && xfrac > 0. && xfrac < 1.) {
1960  xfrac_accepted = true;
1961  }
1962  }
1963 
1964  /* Determine kinematics of the fragmented baryon (antibaryon) and
1965  * remaining (mesonic) string.
1966  * It begins with the lightcone momentum, p^+ (ppos) and p^- (pneg). */
1967  const double ppos_frag = xfrac * mString / sqrt2_;
1968  const double pneg_frag = 0.5 * mTrn_frag * mTrn_frag / ppos_frag;
1969  ppos_string_new = mString / sqrt2_ - ppos_frag;
1970  pneg_string_new = mString / sqrt2_ - pneg_frag;
1971  mTrn_string_new =
1972  std::sqrt(std::max(0., 2. * ppos_string_new * pneg_string_new));
1973 
1974  /* Quark (antiquark) constituents of the remaining string are
1975  * different from those of the original string.
1976  * Therefore, the constituent masses have to be updated. */
1977  for (int i = 0; i < 2; i++) {
1978  m_const[i] = pythia_hadron_->particleData.m0(idqIn[i]);
1979  }
1980  if (bstring > 0) { // in the case of baryonic string
1981  /* Quark is coming from the original string
1982  * and therefore has zero transverse momentum. */
1983  m_trans[0] = m_const[0];
1984  /* Antiquark is coming from the newly produced qqbar pair
1985  * and therefore has transverse momentum, which is opposite to
1986  * that of the fragmented (leading) antibaryon. */
1987  m_trans[1] = std::sqrt(m_const[1] * m_const[1] + QTrn * QTrn);
1988  } else { // in the case of anti-baryonic string
1989  /* Quark is coming from the newly produced qqbar pair
1990  * and therefore has transverse momentum, which is opposite to
1991  * that of the fragmented (leading) baryon. */
1992  m_trans[0] = std::sqrt(m_const[0] * m_const[0] + QTrn * QTrn);
1993  /* Antiquark is coming from the original string
1994  * and therefore has zero transverse momentum. */
1995  m_trans[1] = m_const[1];
1996  }
1997 
1998  if (mTrn_string_new > m_trans[0] + m_trans[1]) {
1999  found_leading_baryon = true;
2000 
2001  FourVector mom_frag((ppos_frag + pneg_frag) / sqrt2_,
2002  evec_basis[0] * (ppos_frag - pneg_frag) / sqrt2_ +
2003  evec_basis[1] * QTrx + evec_basis[2] * QTry);
2004  log.debug("appending the leading baryon ", pdgid_frag,
2005  " to the intermediate particle list.");
2006  /* If the remaining string has enough transverse mass,
2007  * add fragmented baryon (antibaryon) to the intermediate list. */
2008  bool found_ptype = append_intermediate_list(pdgid_frag, mom_frag,
2009  intermediate_particles);
2010  if (!found_ptype) {
2011  log.error("PDG ID ", pdgid_frag, " should exist in ParticleType.");
2012  throw std::runtime_error("string fragmentation failed.");
2013  }
2014  number_of_fragments++;
2015  log.debug("proceed to the next step");
2016  }
2017 
2018  if (iiter == niter_max) {
2019  return 0;
2020  }
2021  iiter += 1;
2022  }
2023 
2024  // lightcone momentum p^+ of the quark constituents on the string ends
2025  std::array<double, 2> ppos_parton;
2026  // lightcone momentum p^- of the quark constituents on the string ends
2027  std::array<double, 2> pneg_parton;
2028 
2029  /* lightcone momenta of the string ends (quark and antiquark)
2030  * First, obtain ppos_parton[0] and ppos_parton[1]
2031  * (p^+ of quark and antiquark) by solving the following equations.
2032  * ppos_string_new = ppos_parton[0] + ppos_parton[1]
2033  * pneg_string_new = 0.5 * m_trans[0] * m_trans[0] / ppos_parton[0] +
2034  * 0.5 * m_trans[1] * m_trans[1] / ppos_parton[1] */
2035  const double pb_const =
2036  (mTrn_string_new * mTrn_string_new + m_trans[0] * m_trans[0] -
2037  m_trans[1] * m_trans[1]) /
2038  (4. * pneg_string_new);
2039  const double pc_const =
2040  0.5 * m_trans[0] * m_trans[0] * ppos_string_new / pneg_string_new;
2041  ppos_parton[0] = pb_const + (bstring > 0 ? -1. : 1.) *
2042  std::sqrt(pb_const * pb_const - pc_const);
2043  ppos_parton[1] = ppos_string_new - ppos_parton[0];
2044  /* Then, compute pneg_parton[0] and pneg_parton[1]
2045  * (p^- of quark and antiquark) from the dispersion relation.
2046  * 2 p^+ p^- = m_transverse^2 */
2047  for (int i = 0; i < 2; i++) {
2048  pneg_parton[i] = 0.5 * m_trans[i] * m_trans[i] / ppos_parton[i];
2049  }
2050 
2051  const int status = 1;
2052  int color, anticolor;
2053  ThreeVector three_mom;
2054  Pythia8::Vec4 pquark;
2055 
2056  // quark end of the remaining (mesonic) string
2057  color = 1;
2058  anticolor = 0;
2059  if (bstring > 0) { // in the case of baryonic string
2060  /* Quark is coming from the original string
2061  * and therefore has zero transverse momentum. */
2062  three_mom = evec_basis[0] * (ppos_parton[0] - pneg_parton[0]) / sqrt2_;
2063  } else { // in the case of anti-baryonic string
2064  /* Quark is coming from the newly produced qqbar pair
2065  * and therefore has transverse momentum, which is opposite to
2066  * that of the fragmented (leading) baryon. */
2067  three_mom = evec_basis[0] * (ppos_parton[0] - pneg_parton[0]) / sqrt2_ -
2068  evec_basis[1] * QTrx - evec_basis[2] * QTry;
2069  }
2070  pquark = set_Vec4((ppos_parton[0] + pneg_parton[0]) / sqrt2_, three_mom);
2071  pSum += pquark;
2072  pythia_hadron_->event.append(idqIn[0], status, color, anticolor, pquark,
2073  m_const[0]);
2074 
2075  // antiquark end of the remaining (mesonic) string
2076  color = 0;
2077  anticolor = 1;
2078  if (bstring > 0) { // in the case of baryonic string
2079  /* Antiquark is coming from the newly produced qqbar pair
2080  * and therefore has transverse momentum, which is opposite to
2081  * that of the fragmented (leading) antibaryon. */
2082  three_mom = evec_basis[0] * (ppos_parton[1] - pneg_parton[1]) / sqrt2_ -
2083  evec_basis[1] * QTrx - evec_basis[2] * QTry;
2084  } else { // in the case of anti-baryonic string
2085  /* Antiquark is coming from the original string
2086  * and therefore has zero transverse momentum. */
2087  three_mom = evec_basis[0] * (ppos_parton[1] - pneg_parton[1]) / sqrt2_;
2088  }
2089  pquark = set_Vec4((ppos_parton[1] + pneg_parton[1]) / sqrt2_, three_mom);
2090  pSum += pquark;
2091  pythia_hadron_->event.append(idqIn[1], status, color, anticolor, pquark,
2092  m_const[1]);
2093  } else {
2094  /* diquark (anti-quark) with PDG id idq2 is going in the direction of
2095  * evecLong.
2096  * quark with PDG id idq1 is going in the direction opposite to evecLong. */
2097  double sign_direction = 1.;
2098  if (bstring == -3) { // anti-baryonic string
2099  /* anti-diquark with PDG id idq1 is going in the direction of evecLong.
2100  * anti-quark with PDG id idq2 is going in the direction
2101  * opposite to evecLong. */
2102  sign_direction = -1;
2103  }
2104 
2105  // evaluate momenta of quarks
2106  const double pCMquark = pCM(mString, m_const[0], m_const[1]);
2107  const double E1 = std::sqrt(m_const[0] * m_const[0] + pCMquark * pCMquark);
2108  const double E2 = std::sqrt(m_const[1] * m_const[1] + pCMquark * pCMquark);
2109 
2110  ThreeVector direction = sign_direction * evecLong;
2111 
2112  // For status and (anti)color see \iref{Sjostrand:2007gs}.
2113  const int status1 = 1, color1 = 1, anticolor1 = 0;
2114  Pythia8::Vec4 pquark = set_Vec4(E1, -direction * pCMquark);
2115  pSum += pquark;
2116  pythia_hadron_->event.append(idqIn[0], status1, color1, anticolor1, pquark,
2117  m_const[0]);
2118 
2119  const int status2 = 1, color2 = 0, anticolor2 = 1;
2120  pquark = set_Vec4(E2, direction * pCMquark);
2121  pSum += pquark;
2122  pythia_hadron_->event.append(idqIn[1], status2, color2, anticolor2, pquark,
2123  m_const[1]);
2124  }
2125 
2126  log.debug("fragmenting a string with ", idqIn[0], ", ", idqIn[1]);
2127  // implement PYTHIA fragmentation
2128  pythia_hadron_->event[0].p(pSum);
2129  pythia_hadron_->event[0].m(pSum.mCalc());
2130  const bool successful_hadronization = pythia_hadron_->next();
2131  if (successful_hadronization) {
2132  for (int ipyth = 0; ipyth < pythia_hadron_->event.size(); ipyth++) {
2133  if (!pythia_hadron_->event[ipyth].isFinal()) {
2134  continue;
2135  }
2136  int pythia_id = pythia_hadron_->event[ipyth].id();
2137  /* K_short and K_long need are converted to K0
2138  * since SMASH only knows K0 */
2139  convert_KaonLS(pythia_id);
2140  FourVector momentum(
2141  pythia_hadron_->event[ipyth].e(), pythia_hadron_->event[ipyth].px(),
2142  pythia_hadron_->event[ipyth].py(), pythia_hadron_->event[ipyth].pz());
2143  log.debug("appending the fragmented hadron ", pythia_id,
2144  " to the intermediate particle list.");
2145  bool found_ptype =
2146  append_intermediate_list(pythia_id, momentum, intermediate_particles);
2147  if (!found_ptype) {
2148  log.warn("PDG ID ", pythia_id,
2149  " does not exist in ParticleType - start over.");
2150  intermediate_particles.clear();
2151  return 0;
2152  }
2153 
2154  number_of_fragments++;
2155  }
2156 
2157  return number_of_fragments;
2158  } else {
2159  return 0;
2160  }
2161 }
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:34
static Pythia8::Vec4 set_Vec4(double energy, const ThreeVector &mom)
Easy setter of Pythia Vec4 from SMASH.
std::unique_ptr< Pythia8::Pythia > pythia_hadron_
PYTHIA object used in fragmentation.
static bool append_intermediate_list(int pdgid, FourVector momentum, ParticleList &intermediate_particles)
append new particle from PYTHIA to a specific particle list
static const ParticleTypeList & list_all()
Definition: particletype.cc:55
static void make_orthonormal_basis(ThreeVector &evec_polar, std::array< ThreeVector, 3 > &evec_basis)
compute three orthonormal basis vectors from unit vector in the longitudinal direction ...
T uniform_int(T min, T max)
Definition: random.h:97
double leading_frag_mean_
mean value of the fragmentation function for the leading baryons
T uniform(T min, T max)
Definition: random.h:85
static void convert_KaonLS(int &pythia_id)
convert Kaon-L or Kaon-S into K0 or Anti-K0
double leading_frag_width_
width of the fragmentation function for the leading baryons
double sqrt2_
square root of 2 ( )
double normal(const T &mean, const T &sigma)
Returns a random number drawn from a normal distribution.
Definition: random.h:247
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
Here is the call graph for this function:
Here is the caller graph for this function:

◆ assign_all_scaling_factors()

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

Assign a cross section scaling factor to all outgoing particles.

The factor is only non-zero, when the outgoing particle carries a valence quark from the excited hadron. The assigned cross section scaling factor is equal to the number of the valence quarks from the fragmented hadron contained in the fragment divided by the total number of valence quarks of that fragment multiplied by a coherence factor

Parameters
[in]baryon_stringbaryon number of the string
[out]outgoing_particleslist of string fragments to which scaling factors are assigned
[in]evecLongdirection in which the string is stretched
[in]suppression_factoradditional coherence factor to be multiplied with scaling factor

Definition at line 2195 of file processstring.cc.

2198  {
2199  // Set each particle's cross section scaling factor to 0 first
2200  for (ParticleData &data : outgoing_particles) {
2201  data.set_cross_section_scaling_factor(0.0);
2202  }
2203  // sort outgoing particles according to the longitudinal velocity
2204  std::sort(outgoing_particles.begin(), outgoing_particles.end(),
2205  [&](ParticleData i, ParticleData j) {
2206  return i.momentum().velocity() * evecLong >
2207  j.momentum().velocity() * evecLong;
2208  });
2209  int nq1, nq2; // number of quarks at both ends of the string
2210  switch (baryon_string) {
2211  case 0:
2212  nq1 = -1;
2213  nq2 = 1;
2214  break;
2215  case 1:
2216  nq1 = 2;
2217  nq2 = 1;
2218  break;
2219  case -1:
2220  nq1 = -2;
2221  nq2 = -1;
2222  break;
2223  default:
2224  throw std::runtime_error("string is neither mesonic nor baryonic");
2225  }
2226  // Try to find nq1 on one string end and nq2 on the other string end and the
2227  // other way around. When the leading particles are close to the string ends,
2228  // the quarks are assumed to be distributed this way.
2229  std::pair<int, int> i = find_leading(nq1, nq2, outgoing_particles);
2230  std::pair<int, int> j = find_leading(nq2, nq1, outgoing_particles);
2231  if (baryon_string == 0 && i.second - i.first < j.second - j.first) {
2232  assign_scaling_factor(nq2, outgoing_particles[j.first], suppression_factor);
2233  assign_scaling_factor(nq1, outgoing_particles[j.second],
2234  suppression_factor);
2235  } else {
2236  assign_scaling_factor(nq1, outgoing_particles[i.first], suppression_factor);
2237  assign_scaling_factor(nq2, outgoing_particles[i.second],
2238  suppression_factor);
2239  }
2240 }
static void assign_scaling_factor(int nquark, ParticleData &data, double suppression_factor)
Assign a cross section scaling factor to the given particle.
static std::pair< int, int > find_leading(int nq1, int nq2, ParticleList &list)
Find the leading string fragments.
Here is the call graph for this function:
Here is the caller graph for this function:

◆ find_leading()

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

Find the leading string fragments.

Find the first particle, which can carry nq1, and the last particle, which can carry nq2 valence quarks and return their indices in the given list.

Parameters
[in]nq1number of valence quarks from excited hadron at forward end of the string
[in]nq2number of valance quarks from excited hadron at backward end of the string
[in]listlist of string fragments
Returns
indices of the leading hadrons in list

Definition at line 2178 of file processstring.cc.

2179  {
2180  assert(list.size() >= 2);
2181  int end = list.size() - 1;
2182  int i1, i2;
2183  for (i1 = 0;
2184  i1 <= end && !list[i1].pdgcode().contains_enough_valence_quarks(nq1);
2185  i1++) {
2186  }
2187  for (i2 = end;
2188  i2 >= 0 && !list[i2].pdgcode().contains_enough_valence_quarks(nq2);
2189  i2--) {
2190  }
2191  std::pair<int, int> indices(i1, i2);
2192  return indices;
2193 }
Here is the caller graph for this function:

◆ assign_scaling_factor()

void smash::StringProcess::assign_scaling_factor ( int  nquark,
ParticleData data,
double  suppression_factor 
)
static

Assign a cross section scaling factor to the given particle.

The scaling factor is the number of quarks from the excited hadron, that the fragment carries devided by the total number of quarks in this fragment multiplied by coherence factor.

Parameters
[in]nquarknumber of valence quarks from the excited hadron contained in the given string fragment data
[out]dataparticle to assign a scaling factor to
[in]suppression_factorcoherence factor to decrease scaling factor

Definition at line 2163 of file processstring.cc.

2164  {
2165  int nbaryon = data.pdgcode().baryon_number();
2166  if (nbaryon == 0) {
2167  // Mesons always get a scaling factor of 1/2 since there is never
2168  // a q-qbar pair at the end of a string so nquark is always 1
2169  data.set_cross_section_scaling_factor(0.5 * suppression_factor);
2170  } else if (data.is_baryon()) {
2171  // Leading baryons get a factor of 2/3 if they carry 2
2172  // and 1/3 if they carry 1 of the strings valence quarks
2173  data.set_cross_section_scaling_factor(suppression_factor * nquark /
2174  (3.0 * nbaryon));
2175  }
2176 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ pdg_map_for_pythia()

int smash::StringProcess::pdg_map_for_pythia ( PdgCode pdg)
static

Take pdg code and map onto particle specie which can be handled by PYTHIA.

Positively charged baryons are mapped onto proton and other baryons are mapped onto neutrons. Same rule applies for anti-baryons. Positively (negatively) charged mesons are mapped onto pi+ (pi-). Negatively and positively charged leptons are mapped respectivly onto electron and positron. Currently, we do not have cross sections for leptons and photons with high energy, so such collisions should not happen.

Parameters
[in]pdgPdgCode that will be mapped
Returns
mapped PDG id to be used in PYTHIA
Exceptions
std::runtime_errorif the incoming particle is neither hadron nor lepton.

Definition at line 2242 of file processstring.cc.

2242  {
2243  PdgCode pdg_mapped(0x0);
2244 
2245  if (pdg.baryon_number() == 1) { // baryon
2246  pdg_mapped = pdg.charge() > 0 ? PdgCode(pdg::p) : PdgCode(pdg::n);
2247  } else if (pdg.baryon_number() == -1) { // antibaryon
2248  pdg_mapped = pdg.charge() < 0 ? PdgCode(-pdg::p) : PdgCode(-pdg::n);
2249  } else if (pdg.is_hadron()) { // meson
2250  if (pdg.charge() >= 0) {
2251  pdg_mapped = PdgCode(pdg::pi_p);
2252  } else {
2253  pdg_mapped = PdgCode(pdg::pi_m);
2254  }
2255  } else if (pdg.is_lepton()) { // lepton
2256  pdg_mapped = pdg.charge() < 0 ? PdgCode(0x11) : PdgCode(-0x11);
2257  } else {
2258  throw std::runtime_error("StringProcess::pdg_map_for_pythia failed.");
2259  }
2260 
2261  return pdg_mapped.get_decimal();
2262 }
constexpr int pi_p
π⁺.
constexpr int p
Proton.
constexpr int pi_m
π⁻.
constexpr int n
Neutron.
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ PPosA_

double smash::StringProcess::PPosA_
private

forward lightcone momentum p^{+} of incoming particle A in CM-frame [GeV]

Definition at line 50 of file processstring.h.

◆ PPosB_

double smash::StringProcess::PPosB_
private

forward lightcone momentum p^{+} of incoming particle B in CM-frame [GeV]

Definition at line 52 of file processstring.h.

◆ PNegA_

double smash::StringProcess::PNegA_
private

backward lightcone momentum p^{-} of incoming particle A in CM-frame [GeV]

Definition at line 54 of file processstring.h.

◆ PNegB_

double smash::StringProcess::PNegB_
private

backward lightcone momentum p^{-} of incoming particle B in CM-frame [GeV]

Definition at line 56 of file processstring.h.

◆ massA_

double smash::StringProcess::massA_
private

mass of incoming particle A [GeV]

Definition at line 58 of file processstring.h.

◆ massB_

double smash::StringProcess::massB_
private

mass of incoming particle B [GeV]

Definition at line 60 of file processstring.h.

◆ sqrtsAB_

double smash::StringProcess::sqrtsAB_
private

sqrt of Mandelstam variable s of collision [GeV]

Definition at line 62 of file processstring.h.

◆ PDGcodes_

std::array<PdgCode, 2> smash::StringProcess::PDGcodes_
private

PdgCodes of incoming particles.

Definition at line 64 of file processstring.h.

◆ plab_

std::array<FourVector, 2> smash::StringProcess::plab_
private

momenta of incoming particles in the lab frame [GeV]

Definition at line 66 of file processstring.h.

◆ pcom_

std::array<FourVector, 2> smash::StringProcess::pcom_
private

momenta of incoming particles in the center of mass frame [GeV]

Definition at line 68 of file processstring.h.

◆ ucomAB_

FourVector smash::StringProcess::ucomAB_
private

velocity four vector of the center of mass in the lab frame

Definition at line 70 of file processstring.h.

◆ vcomAB_

ThreeVector smash::StringProcess::vcomAB_
private

velocity three vector of the center of mass in the lab frame

Definition at line 72 of file processstring.h.

◆ evecBasisAB_

std::array<ThreeVector, 3> smash::StringProcess::evecBasisAB_
private

Orthonormal basis vectors in the center of mass frame, where the 0th one is parallel to momentum of incoming particle A.

Definition at line 77 of file processstring.h.

◆ NpartFinal_

int smash::StringProcess::NpartFinal_
private

total number of final state particles

Definition at line 79 of file processstring.h.

◆ NpartString_

std::array<int, 2> smash::StringProcess::NpartString_
private

number of particles fragmented from strings

Definition at line 81 of file processstring.h.

◆ pmin_gluon_lightcone_

double smash::StringProcess::pmin_gluon_lightcone_
private

the minimum lightcone momentum scale carried by a gluon [GeV]

Definition at line 83 of file processstring.h.

◆ pow_fgluon_beta_

double smash::StringProcess::pow_fgluon_beta_
private

parameter \(\beta\) for the gluon distribution function \( P(x) = x^{-1} (1 - x)^{1 + \beta} \)

Definition at line 88 of file processstring.h.

◆ pow_fquark_alpha_

double smash::StringProcess::pow_fquark_alpha_
private

parameter \(\alpha\) for the quark distribution function \( P(x) = x^{\alpha - 1} (1 - x)^{\beta - 1} \)

Definition at line 93 of file processstring.h.

◆ pow_fquark_beta_

double smash::StringProcess::pow_fquark_beta_
private

parameter \(\beta\) for the quark distribution function \( P(x) = x^{\alpha - 1} (1 - x)^{\beta - 1} \)

Definition at line 98 of file processstring.h.

◆ sigma_qperp_

double smash::StringProcess::sigma_qperp_
private

Transverse momentum spread of the excited strings.

[GeV] Transverse momenta of strings are sampled according to gaussian distribution with width sigma_qperp_

Definition at line 104 of file processstring.h.

◆ leading_frag_mean_

double smash::StringProcess::leading_frag_mean_
private

mean value of the fragmentation function for the leading baryons

Definition at line 106 of file processstring.h.

◆ leading_frag_width_

double smash::StringProcess::leading_frag_width_
private

width of the fragmentation function for the leading baryons

Definition at line 108 of file processstring.h.

◆ kappa_tension_string_

double smash::StringProcess::kappa_tension_string_
private

string tension [GeV/fm]

Definition at line 110 of file processstring.h.

◆ additional_xsec_supp_

double smash::StringProcess::additional_xsec_supp_
private

additional cross-section suppression factor to take coherence effect into account.

Definition at line 115 of file processstring.h.

◆ time_formation_const_

double smash::StringProcess::time_formation_const_
private

constant proper time in the case of constant formation time [fm]

Definition at line 117 of file processstring.h.

◆ soft_t_form_

double smash::StringProcess::soft_t_form_
private

factor to be multiplied to formation times in soft strings

Definition at line 119 of file processstring.h.

◆ time_collision_

double smash::StringProcess::time_collision_
private

time of collision in the computational frame [fm]

Definition at line 121 of file processstring.h.

◆ use_yoyo_model_

bool smash::StringProcess::use_yoyo_model_
private

Whether to calculate the string formation times from the yoyo-model.

Definition at line 123 of file processstring.h.

◆ prob_proton_to_d_uu_

double smash::StringProcess::prob_proton_to_d_uu_
private

Probability of splitting a nucleon into the quark flavour it has only once and a diquark it has twice.

Definition at line 128 of file processstring.h.

◆ sqrt2_

double smash::StringProcess::sqrt2_
private

square root of 2 ( \(\sqrt{2}\))

Definition at line 130 of file processstring.h.

◆ pythia_parton_initialized_

bool smash::StringProcess::pythia_parton_initialized_ = false
private

Remembers if Pythia is initialized or not.

Definition at line 133 of file processstring.h.

◆ final_state_

ParticleList smash::StringProcess::final_state_
private

final state array which must be accessed after the collision

Definition at line 139 of file processstring.h.

◆ pythia_parton_

std::unique_ptr<Pythia8::Pythia> smash::StringProcess::pythia_parton_
private

PYTHIA object used in hard string routine.

Definition at line 142 of file processstring.h.

◆ pythia_hadron_

std::unique_ptr<Pythia8::Pythia> smash::StringProcess::pythia_hadron_
private

PYTHIA object used in fragmentation.

Definition at line 145 of file processstring.h.

◆ pythia_sigmatot_

Pythia8::SigmaTotal smash::StringProcess::pythia_sigmatot_
private

An object to compute cross-sections.

Definition at line 148 of file processstring.h.

◆ event_intermediate_

Pythia8::Event smash::StringProcess::event_intermediate_
private

event record for intermediate partonic state in the hard string routine

Definition at line 154 of file processstring.h.


The documentation for this class was generated from the following files: