#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_tension	value of kappa_tension_string_ [GeV/fm]
[in]	time_formation	value of time_formation_const_ [fm]
[in]	gluon_beta	value of pow_fgluon_beta_
[in]	gluon_pmin	value of pmin_gluon_lightcone_
[in]	quark_alpha	value of pow_fquark_alpha_
[in]	quark_beta	value of pow_fquark_beta_
[in]	strange_supp	strangeness suppression factor (StringFlav:probStoUD) in fragmentation
[in]	diquark_supp	diquark suppression factor (StringFlav:probQQtoQ) in fragmentation
[in]	sigma_perp	value of sigma_qperp_ [GeV]
[in]	leading_frag_mean	Mean of the Gaussian used as fragmentation function for leading baryons.
[in]	leading_frag_width	Width of the Gaussian used as fragmentation function for leading baryons.
[in]	stringz_a	parameter (StringZ:aLund) for the fragmentation function
[in]	stringz_b	parameter (StringZ:bLund) for the fragmentation function [GeV^-2]
[in]	string_sigma_T	transverse momentum spread (StringPT:sigma) in fragmentation [GeV]
[in]	factor_t_form	to be multiplied to soft string formation times
[in]	use_yoyo_model	Whether to calculate formation times from the yoyo-model.
[in]	prob_proton_to_d_uu	Probability 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.

     : pmin_gluon_lightcone_(gluon_pmin),
       pow_fgluon_beta_(gluon_beta),
       pow_fquark_alpha_(quark_alpha),
       pow_fquark_beta_(quark_beta),
       sigma_qperp_(sigma_perp),
       leading_frag_mean_(leading_frag_mean),
       leading_frag_width_(leading_frag_width),
       kappa_tension_string_(string_tension),
       additional_xsec_supp_(0.7),
       time_formation_const_(time_formation),
       soft_t_form_(factor_t_form),
       time_collision_(0.),
       use_yoyo_model_(use_yoyo_model),
       prob_proton_to_d_uu_(prob_proton_to_d_uu) {
   // setup and initialize pythia for hard string process
   pythia_parton_ = make_unique<Pythia8::Pythia>(PYTHIA_XML_DIR, false);
   /* select only non-diffractive events
    * diffractive ones are implemented in a separate routine */
   pythia_parton_->readString("SoftQCD:nonDiffractive = on");
   pythia_parton_->readString("MultipartonInteractions:pTmin = 1.5");
   pythia_parton_->readString("HadronLevel:all = off");
   common_setup_pythia(pythia_parton_.get(), strange_supp, diquark_supp,
                       stringz_a, stringz_b, string_sigma_T);
 
   // setup and initialize pythia for fragmentation
   pythia_hadron_ = make_unique<Pythia8::Pythia>(PYTHIA_XML_DIR, false);
   /* turn off all parton-level processes to implement only hadronization */
   pythia_hadron_->readString("ProcessLevel:all = off");
   common_setup_pythia(pythia_hadron_.get(), strange_supp, diquark_supp,
                       stringz_a, stringz_b, string_sigma_T);
 
   /* initialize PYTHIA */
   pythia_hadron_->init();
   pythia_sigmatot_.init(&pythia_hadron_->info, pythia_hadron_->settings,
                         &pythia_hadron_->particleData, &pythia_hadron_->rndm);
   event_intermediate_.init("intermediate partons",
                            &pythia_hadron_->particleData);
 
   sqrt2_ = std::sqrt(2.);
 
   for (int imu = 0; imu < 3; imu++) {
     evecBasisAB_[imu] = ThreeVector(0., 0., 0.);
   }
 
   final_state_.clear();
 }

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_in	pointer to the PYTHIA object
[in]	strange_supp	strangeness suppression factor (StringFlav:probStoUD) in fragmentation
[in]	diquark_supp	diquark suppression factor (StringFlav:probQQtoQ) in fragmentation
[in]	stringz_a	parameter (StringZ:aLund) for the fragmentation function
[in]	stringz_b	parameter (StringZ:bLund) for the fragmentation function
[in]	string_sigma_T	transverse 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.

                                                                {
   // choose parametrization for mass-dependent width
   pythia_in->readString("ParticleData:modeBreitWigner = 4");
   /* choose minimum transverse momentum scale
    * involved in partonic interactions */
   pythia_in->readString("MultipartonInteractions:pTmin = 1.5");
   pythia_in->readString("MultipartonInteractions:nSample = 10000");
   // transverse momentum spread in string fragmentation
   pythia_in->readString("StringPT:sigma = " + std::to_string(string_sigma_T));
   // diquark suppression factor in string fragmentation
   pythia_in->readString("StringFlav:probQQtoQ = " +
                         std::to_string(diquark_supp));
   // strangeness suppression factor in string fragmentation
   pythia_in->readString("StringFlav:probStoUD = " +
                         std::to_string(strange_supp));
   // parameters for the fragmentation function
   pythia_in->readString("StringZ:aLund = " + std::to_string(stringz_a));
   pythia_in->readString("StringZ:bLund = " + std::to_string(stringz_b));
 
   // manually set the parton distribution function
   pythia_in->readString("PDF:pSet = 13");
   pythia_in->readString("PDF:pSetB = 13");
   pythia_in->readString("PDF:piSet = 1");
   pythia_in->readString("PDF:piSetB = 1");
   pythia_in->readString("Beams:idA = 2212");
   pythia_in->readString("Beams:idB = 2212");
   pythia_in->readString("Beams:eCM = 10.");
 
   // set PYTHIA random seed from outside
   pythia_in->readString("Random:setSeed = on");
   // suppress unnecessary output
   pythia_in->readString("Print:quiet = on");
   // No resonance decays, since the resonances will be handled by SMASH
   pythia_in->readString("HadronLevel:Decay = off");
   // set particle masses and widths in PYTHIA to be same with those in SMASH
   for (auto &ptype : ParticleType::list_all()) {
     int pdgid = ptype.pdgcode().get_decimal();
     double mass_pole = ptype.mass();
     double width_pole = ptype.width_at_pole();
     // check if the particle species is in PYTHIA
     if (pythia_in->particleData.isParticle(pdgid)) {
       // set mass and width in PYTHIA
       pythia_in->particleData.m0(pdgid, mass_pole);
       pythia_in->particleData.mWidth(pdgid, width_pole);
     } else if (pdgid == 310 || pdgid == 130) {
       // set mass and width of Kaon-L and Kaon-S
       pythia_in->particleData.m0(pdgid, kaon_mass);
       pythia_in->particleData.mWidth(pdgid, 0.);
     }
   }
 
   // make energy-momentum conservation in PYTHIA more precise
   pythia_in->readString("Check:epTolErr = 1e-6");
   pythia_in->readString("Check:epTolWarn = 1e-8");
 }

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.

                                  {
     const auto &log = logger<LogArea::Pythia>();
     const int seed_new =
         random::uniform_int(1, maximum_rndm_seed_in_pythia);
 
     pythia_hadron_->rndm.init(seed_new);
     log.debug("pythia_hadron_ : rndm is initialized with seed ", seed_new);
   }

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(); }

smash::StringProcess::pythia_parton_

std::unique_ptr< Pythia8::Pythia > pythia_parton_

PYTHIA object used in hard string routine.

Definition: processstring.h:142

◆ 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_a	pdg code of incoming particle A
[in]	pdg_b	pdg code of incoming particle B
[in]	sqrt_s	collision 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.

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

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_min a value that we want to use for pmin_gluon_lightcone_.

Definition at line 294 of file processstring.h.

                                                          {
     pmin_gluon_lightcone_ = p_light_cone_min;
   }

◆ 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

betapow is a value that we want to use for pow_fgluon_beta_.

Definition at line 304 of file processstring.h.

304 { pow_fgluon_beta_ = betapow; }

smash::StringProcess::pow_fgluon_beta_

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

alphapow	is a value that we want to use for pow_fquark_alpha_.
betapow	is a value that we want to use for pow_fquark_beta_.

Definition at line 313 of file processstring.h.

                                                        {
     pow_fquark_alpha_ = alphapow;
     pow_fquark_beta_ = betapow;
   }

◆ 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_qperp is a value that we want to use for sigma_qperp_.

Definition at line 321 of file processstring.h.

321 { sigma_qperp_ = sigma_qperp; }

smash::StringProcess::sigma_qperp_

double sigma_qperp_

Transverse momentum spread of the excited strings.

Definition: processstring.h:104

◆ 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_string is a value that we want to use for string tension.

Definition at line 326 of file processstring.h.

                                                {
     kappa_tension_string_ = kappa_string;
   }

◆ 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]	incoming	is the list of initial state particles.
[in]	tcoll	is time of collision.

Definition at line 220 of file processstring.cc.

                                                                    {
   PDGcodes_[0] = incoming[0].pdgcode();
   PDGcodes_[1] = incoming[1].pdgcode();
   massA_ = incoming[0].effective_mass();
   massB_ = incoming[1].effective_mass();
 
   plab_[0] = incoming[0].momentum();
   plab_[1] = incoming[1].momentum();
 
   // compute sqrts and velocity of the center of mass.
   sqrtsAB_ = (plab_[0] + plab_[1]).abs();
   ucomAB_ = (plab_[0] + plab_[1]) / sqrtsAB_;
   vcomAB_ = ucomAB_.velocity();
 
   pcom_[0] = plab_[0].LorentzBoost(vcomAB_);
   pcom_[1] = plab_[1].LorentzBoost(vcomAB_);
 
   const double pabscomAB = pCM(sqrtsAB_, massA_, massB_);
   ThreeVector evec_coll = pcom_[0].threevec() / pabscomAB;
   make_orthonormal_basis(evec_coll, evecBasisAB_);
 
   compute_incoming_lightcone_momenta();
 
   time_collision_ = tcoll;
 }

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_polar	unit three-vector in the longitudinal direction
[out]	evec_basis	orthonormal 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 1638 of file processstring.cc.

                                                                    {
   if (std::abs(evec_polar.x3()) < (1. - 1.0e-8)) {
     double ex, ey, et;
     double theta, phi;
 
     // evec_basis[0] is set to be longitudinal direction
     evec_basis[0] = evec_polar;
 
     theta = std::acos(evec_basis[0].x3());
 
     ex = evec_basis[0].x1();
     ey = evec_basis[0].x2();
     et = std::sqrt(ex * ex + ey * ey);
     if (ey > 0.) {
       phi = std::acos(ex / et);
     } else {
       phi = -std::acos(ex / et);
     }
 
     /* The transverse plane is spanned
      * by evec_basis[1] and evec_basis[2]. */
     evec_basis[1].set_x1(cos(theta) * cos(phi));
     evec_basis[1].set_x2(cos(theta) * sin(phi));
     evec_basis[1].set_x3(-sin(theta));
 
     evec_basis[2].set_x1(-sin(phi));
     evec_basis[2].set_x2(cos(phi));
     evec_basis[2].set_x3(0.);
   } else {
     // if evec_polar is very close to the z axis
     if (evec_polar.x3() > 0.) {
       evec_basis[1] = ThreeVector(1., 0., 0.);
       evec_basis[2] = ThreeVector(0., 1., 0.);
       evec_basis[0] = ThreeVector(0., 0., 1.);
     } else {
       evec_basis[1] = ThreeVector(0., 1., 0.);
       evec_basis[2] = ThreeVector(1., 0., 0.);
       evec_basis[0] = ThreeVector(0., 0., -1.);
     }
   }
 }

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 1681 of file processstring.cc.

                                                        {
   PPosA_ = (pcom_[0].x0() + evecBasisAB_[0] * pcom_[0].threevec()) / sqrt2_;
   PNegA_ = (pcom_[0].x0() - evecBasisAB_[0] * pcom_[0].threevec()) / sqrt2_;
   PPosB_ = (pcom_[1].x0() + evecBasisAB_[0] * pcom_[1].threevec()) / sqrt2_;
   PNegB_ = (pcom_[1].x0() - evecBasisAB_[0] * pcom_[1].threevec()) / sqrt2_;
 }

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]	quarks	pdg ids of string ends
[in]	pstr_com	4-momenta of strings in the C.o.m. frame [GeV]
[out]	m_str	masses of strings [GeV]
[out]	evec_str	are directions in which strings are stretched.

Returns: whether masses are above the threshold

Definition at line 333 of file processstring.cc.

                                         {
   std::array<bool, 2> found_mass;
   for (int i = 0; i < 2; i++) {
     found_mass[i] = false;
 
     m_str[i] = pstr_com[i].sqr();
     m_str[i] = (m_str[i] > 0.) ? std::sqrt(m_str[i]) : 0.;
     const double threshold = pythia_hadron_->particleData.m0(quarks[i][0]) +
                              pythia_hadron_->particleData.m0(quarks[i][1]);
     // string mass must be larger than threshold set by PYTHIA.
     if (m_str[i] > threshold) {
       found_mass[i] = true;
       /* Determine direction in which string i is stretched.
        * This is set to be same with the collision axis
        * in the center of mass frame.
        * Initial state partons inside incoming hadrons are
        * moving along the collision axis,
        * which is parallel to three momenta of incoming hadrons
        * in the center of mass frame.
        * Given that partons are assumed to be massless,
        * their four momenta are null vectors and parallel to pnull.
        * If we take unit three-vector of prs,
        * which is pnull in the rest frame of string,
        * it would be the direction in which string ends are moving. */
       const ThreeVector mom = pcom_[i].threevec();
       const FourVector pnull(mom.abs(), mom);
       const FourVector prs = pnull.LorentzBoost(pstr_com[i].velocity());
       evec_str[i] = prs.threevec() / prs.threevec().abs();
     }
   }
 
   return found_mass[0] && found_mass[1];
 }

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]	quarks	pdg ids of string ends
[in]	pstr_com	4-momenta of strings in the C.o.m. frame [GeV]
[in]	m_str	masses of strings [GeV]
[out]	evec_str	are directions in which strings are stretched.
[in]	flip_string_ends	is whether or not we randomly switch string ends.
[in]	separate_fragment_baryon	is 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 370 of file processstring.cc.

                                    {
   const std::array<FourVector, 2> ustr_com = {pstr_com[0] / m_str[0],
                                               pstr_com[1] / m_str[1]};
   for (int i = 0; i < 2; i++) {
     ParticleList new_intermediate_particles;
 
     // determine direction in which string i is stretched.
     ThreeVector evec = evec_str[i];
     // perform fragmentation and add particles to final_state.
     int nfrag = fragment_string(quarks[i][0], quarks[i][1], m_str[i], evec,
                                 flip_string_ends, separate_fragment_baryon,
                                 new_intermediate_particles);
     if (nfrag <= 0) {
       NpartString_[i] = 0;
       return false;
     }
     NpartString_[i] =
         append_final_state(new_intermediate_particles, ustr_com[i], evec);
     assert(nfrag == NpartString_[i]);
   }
   if ((NpartString_[0] > 0) && (NpartString_[1] > 0)) {
     NpartFinal_ = NpartString_[0] + NpartString_[1];
     return true;
   }
   return false;
 }

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_AX specifies 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 249 of file processstring.cc.

                                                {
   NpartFinal_ = 0;
   NpartString_[0] = 0;
   NpartString_[1] = 0;
   final_state_.clear();
 
   double massH = is_AB_to_AX ? massA_ : massB_;
   double mstrMin = is_AB_to_AX ? massB_ : massA_;
   double mstrMax = sqrtsAB_ - massH;
 
   int idqX1, idqX2;
   double QTrn, QTrx, QTry;
   double pabscomHX_sqr, massX;
 
   // decompose hadron into quarks
   make_string_ends(is_AB_to_AX ? PDGcodes_[1] : PDGcodes_[0], idqX1, idqX2,
                    prob_proton_to_d_uu_);
   // string mass must be larger than threshold set by PYTHIA.
   mstrMin = pythia_hadron_->particleData.m0(idqX1) +
             pythia_hadron_->particleData.m0(idqX2);
   // this threshold cannot be larger than maximum of allowed string mass.
   if (mstrMin > mstrMax) {
     return false;
   }
   // sample the transverse momentum transfer
   QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
   QTry = random::normal(0., sigma_qperp_ / sqrt2_);
   QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
   /* sample the string mass and
    * evaluate the three-momenta of hadron and string. */
   massX = random::power(-1.0, mstrMin, mstrMax);
   pabscomHX_sqr = pCM_sqr(sqrtsAB_, massH, massX);
   /* magnitude of the three momentum must be larger
    * than the transverse momentum. */
   const bool foundPabsX = pabscomHX_sqr > QTrn * QTrn;
 
   if (!foundPabsX) {
     return false;
   }
   double sign_direction = is_AB_to_AX ? 1. : -1.;
   // determine three momentum of the final state hadron
   const ThreeVector cm_momentum =
       sign_direction *
       (evecBasisAB_[0] * std::sqrt(pabscomHX_sqr - QTrn * QTrn) +
        evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry);
   const FourVector pstrHcom(std::sqrt(pabscomHX_sqr + massH * massH),
                             cm_momentum);
   const FourVector pstrXcom(std::sqrt(pabscomHX_sqr + massX * massX),
                             -cm_momentum);
 
   const FourVector ustrXcom = pstrXcom / massX;
   /* determine direction in which the string is stretched.
    * this is set to be same with the the collision axis
    * in the center of mass frame. */
   const ThreeVector threeMomentum =
       is_AB_to_AX ? pcom_[1].threevec() : pcom_[0].threevec();
   const FourVector pnull = FourVector(threeMomentum.abs(), threeMomentum);
   const FourVector prs = pnull.LorentzBoost(ustrXcom.velocity());
   ThreeVector evec = prs.threevec() / prs.threevec().abs();
   // perform fragmentation and add particles to final_state.
   ParticleList new_intermediate_particles;
   bool separate_fragment_baryon = false;
   int nfrag =
       fragment_string(idqX1, idqX2, massX, evec, true, separate_fragment_baryon,
                       new_intermediate_particles);
   if (nfrag < 1) {
     NpartString_[0] = 0;
     return false;
   }
   NpartString_[0] =
       append_final_state(new_intermediate_particles, ustrXcom, evec);
 
   NpartString_[1] = 1;
   PdgCode hadron_code = is_AB_to_AX ? PDGcodes_[0] : PDGcodes_[1];
   ParticleData new_particle(ParticleType::find(hadron_code));
   new_particle.set_4momentum(pstrHcom);
   new_particle.set_cross_section_scaling_factor(1.);
   new_particle.set_formation_time(time_collision_);
   final_state_.push_back(new_particle);
 
   NpartFinal_ = NpartString_[0] + NpartString_[1];
   return true;
 }

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 403 of file processstring.cc.

                                {
   NpartFinal_ = 0;
   NpartString_[0] = 0;
   NpartString_[1] = 0;
   final_state_.clear();
 
   std::array<std::array<int, 2>, 2> quarks;
   std::array<FourVector, 2> pstr_com;
   std::array<double, 2> m_str;
   std::array<ThreeVector, 2> evec_str;
   ThreeVector threeMomentum;
 
   // decompose hadron into quark (and diquark) contents
   make_string_ends(PDGcodes_[0], quarks[0][0], quarks[0][1],
                    prob_proton_to_d_uu_);
   make_string_ends(PDGcodes_[1], quarks[1][0], quarks[1][1],
                    prob_proton_to_d_uu_);
   // sample the lightcone momentum fraction carried by gluons
   const double xmin_gluon_fraction = pmin_gluon_lightcone_ / sqrtsAB_;
   const double xfracA =
       random::beta_a0(xmin_gluon_fraction, pow_fgluon_beta_ + 1.);
   const double xfracB =
       random::beta_a0(xmin_gluon_fraction, pow_fgluon_beta_ + 1.);
   // sample the transverse momentum transfer
   const double QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
   const double QTry = random::normal(0., sigma_qperp_ / sqrt2_);
   const double QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
   // evaluate the lightcone momentum transfer
   const double QPos = -QTrn * QTrn / (2. * xfracB * PNegB_);
   const double QNeg = QTrn * QTrn / (2. * xfracA * PPosA_);
   // compute four-momentum of string 1
   threeMomentum = evecBasisAB_[0] * (PPosA_ + QPos - PNegA_ - QNeg) / sqrt2_ +
                   evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry;
   pstr_com[0] =
       FourVector((PPosA_ + QPos + PNegA_ + QNeg) / sqrt2_, threeMomentum);
   // compute four-momentum of string 2
   threeMomentum = evecBasisAB_[0] * (PPosB_ - QPos - PNegB_ + QNeg) / sqrt2_ -
                   evecBasisAB_[1] * QTrx - evecBasisAB_[2] * QTry;
   pstr_com[1] =
       FourVector((PPosB_ - QPos + PNegB_ - QNeg) / sqrt2_, threeMomentum);
 
   const bool found_masses =
       set_mass_and_direction_2strings(quarks, pstr_com, m_str, evec_str);
   if (!found_masses) {
     return false;
   }
   const bool flip_string_ends = true;
   const bool separate_fragment_baryon = false;
   const bool success =
       make_final_state_2strings(quarks, pstr_com, m_str, evec_str,
                                 flip_string_ends, separate_fragment_baryon);
   return success;
 }

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_error if incoming particles are neither mesonic nor baryonic

Definition at line 458 of file processstring.cc.

                                    {
   NpartFinal_ = 0;
   NpartString_[0] = 0;
   NpartString_[1] = 0;
   final_state_.clear();
 
   std::array<std::array<int, 2>, 2> quarks;
   std::array<FourVector, 2> pstr_com;
   std::array<double, 2> m_str;
   std::array<ThreeVector, 2> evec_str;
 
   // decompose hadron into quark (and diquark) contents
   int idqA1, idqA2, idqB1, idqB2;
   make_string_ends(PDGcodes_[0], idqA1, idqA2, prob_proton_to_d_uu_);
   make_string_ends(PDGcodes_[1], idqB1, idqB2, prob_proton_to_d_uu_);
 
   const int bar_a = PDGcodes_[0].baryon_number(),
             bar_b = PDGcodes_[1].baryon_number();
   if (bar_a == 1 ||  // baryon-baryon, baryon-meson, baryon-antibaryon
       (bar_a == 0 && bar_b == 1) ||  // meson-baryon
       (bar_a == 0 && bar_b == 0)) {  // meson-meson
     quarks[0][0] = idqB1;
     quarks[0][1] = idqA2;
     quarks[1][0] = idqA1;
     quarks[1][1] = idqB2;
   } else if (((bar_a == 0) && (bar_b == -1)) ||  // meson-antibaryon
              (bar_a == -1)) {  // antibaryon-baryon, antibaryon-meson,
                                // antibaryon-antibaryon
     quarks[0][0] = idqA1;
     quarks[0][1] = idqB2;
     quarks[1][0] = idqB1;
     quarks[1][1] = idqA2;
   } else {
     std::stringstream ss;
     ss << "  StringProcess::next_NDiff : baryonA = " << bar_a
        << ", baryonB = " << bar_b;
     throw std::runtime_error(ss.str());
   }
   // sample the lightcone momentum fraction carried by quarks
   const double xfracA = random::beta(pow_fquark_alpha_, pow_fquark_beta_);
   const double xfracB = random::beta(pow_fquark_alpha_, pow_fquark_beta_);
   // sample the transverse momentum transfer
   const double QTrx = random::normal(0., sigma_qperp_ / sqrt2_);
   const double QTry = random::normal(0., sigma_qperp_ / sqrt2_);
   const double QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
   // evaluate the lightcone momentum transfer
   const double QPos = -QTrn * QTrn / (2. * xfracB * PNegB_);
   const double QNeg = QTrn * QTrn / (2. * xfracA * PPosA_);
   const double dPPos = -xfracA * PPosA_ - QPos;
   const double dPNeg = xfracB * PNegB_ - QNeg;
   // compute four-momentum of string 1
   ThreeVector threeMomentum =
       evecBasisAB_[0] * (PPosA_ + dPPos - PNegA_ - dPNeg) / sqrt2_ +
       evecBasisAB_[1] * QTrx + evecBasisAB_[2] * QTry;
   pstr_com[0] =
       FourVector((PPosA_ + dPPos + PNegA_ + dPNeg) / sqrt2_, threeMomentum);
   m_str[0] = pstr_com[0].sqr();
   // compute four-momentum of string 2
   threeMomentum = evecBasisAB_[0] * (PPosB_ - dPPos - PNegB_ + dPNeg) / sqrt2_ -
                   evecBasisAB_[1] * QTrx - evecBasisAB_[2] * QTry;
   pstr_com[1] =
       FourVector((PPosB_ - dPPos + PNegB_ - dPNeg) / sqrt2_, threeMomentum);
 
   const bool found_masses =
       set_mass_and_direction_2strings(quarks, pstr_com, m_str, evec_str);
   if (!found_masses) {
     return false;
   }
   const bool flip_string_ends = false;
   const bool separate_fragment_baryon = true;
   const bool success =
       make_final_state_2strings(quarks, pstr_com, m_str, evec_str,
                                 flip_string_ends, separate_fragment_baryon);
   return success;
 }

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 535 of file processstring.cc.

                                    {
   const auto &log = logger<LogArea::Pythia>();
   NpartFinal_ = 0;
   final_state_.clear();
 
   log.debug("Hard non-diff. with ", PDGcodes_[0], " + ", PDGcodes_[1],
             " at CM energy [GeV] ", sqrtsAB_);
 
   std::array<int, 2> pdg_for_pythia;
   std::array<std::array<int, 5>, 2> excess_quark;
   std::array<std::array<int, 5>, 2> excess_antiq;
   for (int i = 0; i < 2; i++) {
     for (int j = 0; j < 5; j++) {
       excess_quark[i][j] = 0;
       excess_antiq[i][j] = 0;
     }
 
     // get PDG id used in PYTHIA event generation
     pdg_for_pythia[i] = pdg_map_for_pythia(PDGcodes_[i]);
     log.debug("  incoming particle ", i, " : ", PDGcodes_[i],
               " is mapped onto ", pdg_for_pythia[i]);
 
     PdgCode pdgcode_for_pythia(std::to_string(pdg_for_pythia[i]));
     /* evaluate how many more constituents incoming hadron has
      * compared to the mapped one. */
     find_excess_constituent(PDGcodes_[i], pdgcode_for_pythia, excess_quark[i],
                             excess_antiq[i]);
     log.debug("    excess_quark[", i, "] = (", excess_quark[i][0], ", ",
               excess_quark[i][1], ", ", excess_quark[i][2], ", ",
               excess_quark[i][3], ", ", excess_quark[i][4], ")");
     log.debug("    excess_antiq[", i, "] = (", excess_antiq[i][0], ", ",
               excess_antiq[i][1], ", ", excess_antiq[i][2], ", ",
               excess_antiq[i][3], ", ", excess_antiq[i][4], ")");
   }
 
   int previous_idA = pythia_parton_->mode("Beams:idA"),
       previous_idB = pythia_parton_->mode("Beams:idB");
   double previous_eCM = pythia_parton_->parm("Beams:eCM");
   // check if the initial state for PYTHIA remains same.
   bool same_initial_state = previous_idA == pdg_for_pythia[0] &&
                             previous_idB == pdg_for_pythia[1] &&
                             std::abs(previous_eCM - sqrtsAB_) < really_small;
 
   /* Perform PYTHIA initialization if it was not previously initialized
    * or the initial state changed. */
   if (!pythia_parton_initialized_ || !same_initial_state) {
     pythia_parton_->settings.mode("Beams:idA", pdg_for_pythia[0]);
     pythia_parton_->settings.mode("Beams:idB", pdg_for_pythia[1]);
     pythia_parton_->settings.parm("Beams:eCM", sqrtsAB_);
 
     pythia_parton_initialized_ = pythia_parton_->init();
     log.debug("Pythia initialized with ", pdg_for_pythia[0], " + ",
               pdg_for_pythia[1], " at CM energy [GeV] ", sqrtsAB_);
     if (!pythia_parton_initialized_) {
       throw std::runtime_error("Pythia failed to initialize.");
     }
   }
   /* Set the random seed of the Pythia random Number Generator.
    * Pythia's random is controlled by SMASH in every single collision.
    * In this way we ensure that the results are reproducible
    * for every event if one knows SMASH random seed. */
   const int seed_new = random::uniform_int(1, maximum_rndm_seed_in_pythia);
   pythia_parton_->rndm.init(seed_new);
   log.debug("pythia_parton_ : rndm is initialized with seed ", seed_new);
 
   // Short notation for Pythia event
   Pythia8::Event &event_hadron = pythia_hadron_->event;
   log.debug("Pythia hard event created");
   bool final_state_success = pythia_parton_->next();
   log.debug("Pythia final state computed, success = ", final_state_success);
   if (!final_state_success) {
     return false;
   }
 
   ParticleList new_intermediate_particles;
   ParticleList new_non_hadron_particles;
 
   Pythia8::Vec4 pSum = 0.;
   event_intermediate_.reset();
   /* Update the partonic intermediate state from PYTHIA output.
    * Note that hadronization will be performed separately,
    * after identification of strings and replacement of constituents. */
   for (int i = 0; i < pythia_parton_->event.size(); i++) {
     if (pythia_parton_->event[i].isFinal()) {
       const int pdgid = pythia_parton_->event[i].id();
       Pythia8::Vec4 pquark = pythia_parton_->event[i].p();
       const double mass = pythia_parton_->particleData.m0(pdgid);
 
       const int status = pythia_parton_->event[i].status();
       const int color = pythia_parton_->event[i].col();
       const int anticolor = pythia_parton_->event[i].acol();
 
       pSum += pquark;
       event_intermediate_.append(pdgid, status, color, anticolor, pquark, mass);
     }
   }
   // add junctions to the intermediate state if there is any.
   event_intermediate_.clearJunctions();
   for (int i = 0; i < pythia_parton_->event.sizeJunction(); i++) {
     const int kind = pythia_parton_->event.kindJunction(i);
     std::array<int, 3> col;
     for (int j = 0; j < 3; j++) {
       col[j] = pythia_parton_->event.colJunction(i, j);
     }
     event_intermediate_.appendJunction(kind, col[0], col[1], col[2]);
   }
   /* The zeroth entry of event record is supposed to have the information
    * on the whole system. Specify the total momentum and invariant mass. */
   event_intermediate_[0].p(pSum);
   event_intermediate_[0].m(pSum.mCalc());
 
   /* Replace quark constituents according to the excess of valence quarks
    * and then rescale momenta of partons by constant factor
    * to fulfill the energy-momentum conservation. */
   bool correct_constituents =
       restore_constituent(event_intermediate_, excess_quark, excess_antiq);
   if (!correct_constituents) {
     log.debug("failed to find correct partonic constituents.");
     return false;
   }
 
   int npart = event_intermediate_.size();
   int ipart = 0;
   while (ipart < npart) {
     const int pdgid = event_intermediate_[ipart].id();
     if (event_intermediate_[ipart].isFinal() &&
         !event_intermediate_[ipart].isParton() &&
         !pythia_parton_->particleData.isOctetHadron(pdgid)) {
       log.debug("PDG ID from Pythia: ", pdgid);
       FourVector momentum = reorient(event_intermediate_[ipart], evecBasisAB_);
       log.debug("4-momentum from Pythia: ", momentum);
       bool found_ptype =
           append_intermediate_list(pdgid, momentum, new_non_hadron_particles);
       if (!found_ptype) {
         log.warn("PDG ID ", pdgid,
                  " does not exist in ParticleType - start over.");
         final_state_success = false;
       }
       event_intermediate_.remove(ipart, ipart);
       npart -= 1;
     } else {
       ipart += 1;
     }
   }
 
   bool hadronize_success = false;
   bool find_forward_string = true;
   log.debug("Hard non-diff: partonic process gives ",
             event_intermediate_.size(), " partons.");
   // identify and fragment strings until there is no parton left.
   while (event_intermediate_.size() > 1) {
     // dummy event to initialize the internal variables of PYTHIA.
     pythia_hadron_->event.reset();
     if (!pythia_hadron_->next()) {
       log.debug("  Dummy event in hard string routine failed.");
       hadronize_success = false;
       break;
     }
 
     if (event_intermediate_.sizeJunction() > 0) {
       // identify string from a junction if there is any.
       compose_string_junction(find_forward_string, event_intermediate_,
                               pythia_hadron_->event);
     } else {
       /* identify string from a most forward or backward parton.
        * if there is no junction. */
       compose_string_parton(find_forward_string, event_intermediate_,
                             pythia_hadron_->event);
     }
 
     // fragment the (identified) string into hadrons.
     hadronize_success = pythia_hadron_->forceHadronLevel(false);
     log.debug("Pythia hadronized, success = ", hadronize_success);
 
     new_intermediate_particles.clear();
     if (hadronize_success) {
       for (int i = 0; i < event_hadron.size(); i++) {
         if (event_hadron[i].isFinal()) {
           int pythia_id = event_hadron[i].id();
           log.debug("PDG ID from Pythia: ", pythia_id);
           /* K_short and K_long need to be converted to K0
            * since SMASH only knows K0 */
           convert_KaonLS(pythia_id);
 
           /* evecBasisAB_[0] is a unit 3-vector in the collision axis,
            * while evecBasisAB_[1] and evecBasisAB_[2] spans the transverse
            * plane. Given that PYTHIA assumes z-direction to be the collision
            * axis, pz from PYTHIA should be the momentum compoment in
            * evecBasisAB_[0]. px and py are respectively the momentum components
            * in two transverse directions evecBasisAB_[1] and evecBasisAB_[2].
            */
           FourVector momentum = reorient(event_hadron[i], evecBasisAB_);
           log.debug("4-momentum from Pythia: ", momentum);
           log.debug("appending the particle ", pythia_id,
                     " to the intermediate particle list.");
           bool found_ptype = false;
           if (event_hadron[i].isHadron()) {
             found_ptype = append_intermediate_list(pythia_id, momentum,
                                                    new_intermediate_particles);
           } else {
             found_ptype = append_intermediate_list(pythia_id, momentum,
                                                    new_non_hadron_particles);
           }
           if (!found_ptype) {
             log.warn("PDG ID ", pythia_id,
                      " does not exist in ParticleType - start over.");
             hadronize_success = false;
           }
         }
       }
     }
 
     /* if hadronization is not successful,
      * reset the event records, return false and then start over. */
     if (!hadronize_success) {
       break;
     }
 
     FourVector uString = FourVector(1., 0., 0., 0.);
     ThreeVector evec = find_forward_string ? evecBasisAB_[0] : -evecBasisAB_[0];
     int nfrag = append_final_state(new_intermediate_particles, uString, evec);
     NpartFinal_ += nfrag;
 
     find_forward_string = !find_forward_string;
   }
 
   if (hadronize_success) {
     // add the final state particles, which are not hadron.
     for (ParticleData data : new_non_hadron_particles) {
       data.set_cross_section_scaling_factor(1.);
       data.set_formation_time(time_collision_);
       final_state_.push_back(data);
     }
   } else {
     final_state_.clear();
   }
 
   return hadronize_success;
 }

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_argument if incoming particles are not baryon-antibaryon pair

Definition at line 1526 of file processstring.cc.

                                   {
   const auto &log = logger<LogArea::Pythia>();
   const std::array<FourVector, 2> ustrcom = {FourVector(1., 0., 0., 0.),
                                              FourVector(1., 0., 0., 0.)};
 
   NpartFinal_ = 0;
   NpartString_[0] = 0;
   NpartString_[1] = 0;
   final_state_.clear();
 
   log.debug("Annihilation occurs between ", PDGcodes_[0], "+", PDGcodes_[1],
             " at CM energy [GeV] ", sqrtsAB_);
 
   // check if the initial state is baryon-antibaryon pair.
   PdgCode baryon = PDGcodes_[0], antibaryon = PDGcodes_[1];
   if (baryon.baryon_number() == -1) {
     std::swap(baryon, antibaryon);
   }
   if (baryon.baryon_number() != 1 || antibaryon.baryon_number() != -1) {
     throw std::invalid_argument("Expected baryon-antibaryon pair.");
   }
 
   // Count how many qqbar combinations are possible for each quark type
   constexpr int n_q_types = 5;  // u, d, s, c, b
   std::vector<int> qcount_bar, qcount_antibar;
   std::vector<int> n_combinations;
   bool no_combinations = true;
   for (int i = 0; i < n_q_types; i++) {
     qcount_bar.push_back(baryon.net_quark_number(i + 1));
     qcount_antibar.push_back(-antibaryon.net_quark_number(i + 1));
     const int n_i = qcount_bar[i] * qcount_antibar[i];
     n_combinations.push_back(n_i);
     if (n_i > 0) {
       no_combinations = false;
     }
   }
 
   /* if it is a BBbar pair but there is no qqbar pair to annihilate,
    * nothing happens */
   if (no_combinations) {
     for (int i = 0; i < 2; i++) {
       NpartString_[i] = 1;
       ParticleData new_particle(ParticleType::find(PDGcodes_[i]));
       new_particle.set_4momentum(pcom_[i]);
       new_particle.set_cross_section_scaling_factor(1.);
       new_particle.set_formation_time(time_collision_);
       final_state_.push_back(new_particle);
     }
     NpartFinal_ = NpartString_[0] + NpartString_[1];
     return true;
   }
 
   // Select qqbar pair to annihilate and remove it away
   auto discrete_distr = random::discrete_dist<int>(n_combinations);
   const int q_annihilate = discrete_distr() + 1;
   qcount_bar[q_annihilate - 1]--;
   qcount_antibar[q_annihilate - 1]--;
 
   // Get the remaining quarks and antiquarks
   std::vector<int> remaining_quarks, remaining_antiquarks;
   for (int i = 0; i < n_q_types; i++) {
     for (int j = 0; j < qcount_bar[i]; j++) {
       remaining_quarks.push_back(i + 1);
     }
     for (int j = 0; j < qcount_antibar[i]; j++) {
       remaining_antiquarks.push_back(-(i + 1));
     }
   }
   assert(remaining_quarks.size() == 2);
   assert(remaining_antiquarks.size() == 2);
 
   const std::array<double, 2> mstr = {0.5 * sqrtsAB_, 0.5 * sqrtsAB_};
 
   // randomly select two quark-antiquark pairs
   if (random::uniform_int(0, 1) == 0) {
     std::swap(remaining_quarks[0], remaining_quarks[1]);
   }
   if (random::uniform_int(0, 1) == 0) {
     std::swap(remaining_antiquarks[0], remaining_antiquarks[1]);
   }
   // Make sure it satisfies kinematical threshold constraint
   bool kin_threshold_satisfied = true;
   for (int i = 0; i < 2; i++) {
     const double mstr_min =
         pythia_hadron_->particleData.m0(remaining_quarks[i]) +
         pythia_hadron_->particleData.m0(remaining_antiquarks[i]);
     if (mstr_min > mstr[i]) {
       kin_threshold_satisfied = false;
     }
   }
   if (!kin_threshold_satisfied) {
     return false;
   }
   // Fragment two strings
   for (int i = 0; i < 2; i++) {
     ParticleList new_intermediate_particles;
 
     ThreeVector evec = pcom_[i].threevec() / pcom_[i].threevec().abs();
     const int nfrag =
         fragment_string(remaining_quarks[i], remaining_antiquarks[i], mstr[i],
                         evec, true, false, new_intermediate_particles);
     if (nfrag <= 0) {
       NpartString_[i] = 0;
       return false;
     }
     NpartString_[i] =
         append_final_state(new_intermediate_particles, ustrcom[i], evec);
   }
   NpartFinal_ = NpartString_[0] + NpartString_[1];
   return true;
 }

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_actual	PDG code of actual incoming particle.
[in]	pdg_mapped	PDG code of mapped particles used in PYTHIA event generation.
[out]	excess_quark	excess of quarks.
[out]	excess_antiq	excess of anti-quarks.

Definition at line 775 of file processstring.cc.

                                                                             {
   /* decompose PDG id of the actual hadron and mapped one
    * to get the valence quark constituents */
   std::array<int, 3> qcontent_actual = pdg_actual.quark_content();
   std::array<int, 3> qcontent_mapped = pdg_mapped.quark_content();
 
   excess_quark = {0, 0, 0, 0, 0};
   excess_antiq = {0, 0, 0, 0, 0};
   for (int i = 0; i < 3; i++) {
     if (qcontent_actual[i] > 0) {  // quark content of the actual hadron
       int j = qcontent_actual[i] - 1;
       excess_quark[j] += 1;
     }
 
     if (qcontent_mapped[i] > 0) {  // quark content of the mapped hadron
       int j = qcontent_mapped[i] - 1;
       excess_quark[j] -= 1;
     }
 
     if (qcontent_actual[i] < 0) {  // antiquark content of the actual hadron
       int j = std::abs(qcontent_actual[i]) - 1;
       excess_antiq[j] += 1;
     }
 
     if (qcontent_mapped[i] < 0) {  // antiquark content of the mapped hadron
       int j = std::abs(qcontent_mapped[i]) - 1;
       excess_antiq[j] -= 1;
     }
   }
 }

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]	particle	PYTHIA particle object to be converted.
[out]	excess_constituent	excess 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 809 of file processstring.cc.

                                                                      {
   const auto &log = logger<LogArea::Pythia>();
 
   // If the particle is neither quark nor diquark, nothing to do.
   if (!particle.isQuark() && !particle.isDiquark()) {
     return;
   }
 
   // If there is no excess of constituents, nothing to do.
   const std::array<int, 5> excess_null = {0, 0, 0, 0, 0};
   if (excess_constituent == excess_null) {
     return;
   }
 
   int nq = 0;
   std::array<int, 2> pdgid = {0, 0};
   int spin_deg = 0;
   int pdgid_new = 0;
   if (particle.isQuark()) {
     nq = 1;
     pdgid[0] = particle.id();
   } else if (particle.isDiquark()) {
     nq = 2;
     quarks_from_diquark(particle.id(), pdgid[0], pdgid[1], spin_deg);
   }
 
   for (int iq = 0; iq < nq; iq++) {
     int jq = std::abs(pdgid[iq]) - 1;
     int k_select = 0;
     std::vector<int> k_found;
     k_found.clear();
     // check if the constituent needs to be converted.
     if (excess_constituent[jq] < 0) {
       for (int k = 0; k < 5; k++) {
         // check which specie it can be converted into.
         if (k != jq && excess_constituent[k] > 0) {
           k_found.push_back(k);
         }
       }
     }
 
     // make a random selection of specie and update the excess of constituent.
     if (k_found.size() > 0) {
       const int l =
           random::uniform_int(0, static_cast<int>(k_found.size()) - 1);
       k_select = k_found[l];
       /* flavor jq + 1 is converted into k_select + 1
        * and excess_constituent is updated. */
       pdgid[iq] = pdgid[iq] > 0 ? k_select + 1 : -(k_select + 1);
       excess_constituent[jq] += 1;
       excess_constituent[k_select] -= 1;
     }
   }
 
   // determine PDG id of the converted parton.
   if (particle.isQuark()) {
     pdgid_new = pdgid[0];
   } else if (particle.isDiquark()) {
     if (std::abs(pdgid[0]) < std::abs(pdgid[1])) {
       std::swap(pdgid[0], pdgid[1]);
     }
 
     pdgid_new = std::abs(pdgid[0]) * 1000 + std::abs(pdgid[1]) * 100;
     if (std::abs(pdgid[0]) == std::abs(pdgid[1])) {
       pdgid_new += 3;
     } else {
       pdgid_new += spin_deg;
     }
 
     if (particle.id() < 0) {
       pdgid_new *= -1;
     }
   }
   log.debug("  parton id = ", particle.id(), " is converted to ", pdgid_new);
 
   // update the constituent mass and energy.
   Pythia8::Vec4 pquark = particle.p();
   double mass_new = pythia_hadron_->particleData.m0(pdgid_new);
   double e_new = std::sqrt(mass_new * mass_new + pquark.pAbs() * pquark.pAbs());
   // update the particle object.
   particle.id(pdgid_new);
   particle.e(e_new);
   particle.m(mass_new);
 }

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_intermediate	PYTHIA partonic event record which contains output from PYTHIA (hard) event generation.
[out]	nquark_total	total number of quarks in the system. This is computed based on event_intermediate.
[out]	nantiq_total	total number of antiquarks in the system. This is computed based on event_intermediate.

Definition at line 895 of file processstring.cc.

                                     {
   for (int iflav = 0; iflav < 5; iflav++) {
     nquark_total[iflav] = 0;
     nantiq_total[iflav] = 0;
   }
 
   for (int ip = 1; ip < event_intermediate.size(); ip++) {
     if (!event_intermediate[ip].isFinal()) {
       continue;
     }
     const int pdgid = event_intermediate[ip].id();
     if (pdgid > 0) {
       // quarks
       for (int iflav = 0; iflav < 5; iflav++) {
         nquark_total[iflav] +=
             pythia_hadron_->particleData.nQuarksInCode(pdgid, iflav + 1);
       }
     } else {
       // antiquarks
       for (int iflav = 0; iflav < 5; iflav++) {
         nantiq_total[iflav] += pythia_hadron_->particleData.nQuarksInCode(
             std::abs(pdgid), iflav + 1);
       }
     }
   }
 }

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_intermediate	PYTHIA partonic event record to be updated when a gluon happens to split into a qqbar pair.
[out]	nquark_total	total number of quarks in the system. This is computed based on event_intermediate.
[out]	nantiq_total	total number of antiquarks in the system. This is computed based on event_intermediate.
[in]	sign_constituent	true (false) if want to check quarks (antiquarks) and their excesses.
[in]	excess_constituent	excess 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 924 of file processstring.cc.

                                                        {
   const auto &log = logger<LogArea::Pythia>();
 
   Pythia8::Vec4 pSum = event_intermediate[0].p();
 
   /* compute total number of quark and antiquark constituents
    * in the whole system. */
   find_total_number_constituent(event_intermediate, nquark_total, nantiq_total);
 
   for (int iflav = 0; iflav < 5; iflav++) {
     /* Find how many constituent will be in the system after
      * changing the flavors.
      * Note that nquark_total is number of constituent right after
      * the pythia event (with mapped incoming hadrons), while the excess
      * shows how many constituents we have more or less that nquark_total. */
     int nquark_final =
         excess_constituent[0][iflav] + excess_constituent[1][iflav];
     if (sign_constituent) {
       nquark_final += nquark_total[iflav];
     } else {
       nquark_final += nantiq_total[iflav];
     }
     /* Therefore, nquark_final should not be negative.
      * negative nquark_final means that it will not be possible to
      * find a constituent to change the flavor. */
     bool enough_quark = nquark_final >= 0;
     /* If that is the case, a gluon will be splitted into
      * a quark-antiquark pair with the desired flavor. */
     if (!enough_quark) {
       log.debug("  not enough constituents with flavor ", iflav + 1,
                 " : try to split a gluon to qqbar.");
       for (int ic = 0; ic < std::abs(nquark_final); ic++) {
         /* Since each incoming hadron has its own count of the excess,
          * it is necessary to find which one is problematic. */
         int ih_mod = -1;
         if (excess_constituent[0][iflav] < 0) {
           ih_mod = 0;
         } else {
           ih_mod = 1;
         }
 
         /* find the most forward or backward gluon
          * depending on which incoming hadron is found to be an issue. */
         int iforward = 1;
         for (int ip = 2; ip < event_intermediate.size(); ip++) {
           if (!event_intermediate[ip].isFinal() ||
               !event_intermediate[ip].isGluon()) {
             continue;
           }
 
           const double y_gluon_current = event_intermediate[ip].y();
           const double y_gluon_forward = event_intermediate[iforward].y();
           if ((ih_mod == 0 && y_gluon_current > y_gluon_forward) ||
               (ih_mod == 1 && y_gluon_current < y_gluon_forward)) {
             iforward = ip;
           }
         }
 
         if (!event_intermediate[iforward].isGluon()) {
           log.debug("There is no gluon to split into qqbar.");
           return false;
         }
 
         // four momentum of the original gluon
         Pythia8::Vec4 pgluon = event_intermediate[iforward].p();
 
         const int pdgid = iflav + 1;
         const double mass = pythia_hadron_->particleData.m0(pdgid);
         const int status = event_intermediate[iforward].status();
         /* color and anticolor indices.
          * the color index of gluon goes to the quark, while
          * the anticolor index goes to the antiquark */
         const int col = event_intermediate[iforward].col();
         const int acol = event_intermediate[iforward].acol();
 
         // three momenta of quark and antiquark
         std::array<double, 2> px_quark;
         std::array<double, 2> py_quark;
         std::array<double, 2> pz_quark;
         // energies of quark and antiquark
         std::array<double, 2> e_quark;
         // four momenta of quark and antiquark
         std::array<Pythia8::Vec4, 2> pquark;
         // transverse momentum scale of string fragmentation
         const double sigma_qt_frag = pythia_hadron_->parm("StringPT:sigma");
         // sample relative transverse momentum between quark and antiquark
         const double qx = random::normal(0., sigma_qt_frag / sqrt2_);
         const double qy = random::normal(0., sigma_qt_frag / sqrt2_);
         // setup kinematics
         for (int isign = 0; isign < 2; isign++) {
           /* As the simplest assumption, the three momentum of gluon
            * is equally distributed to quark and antiquark.
            * Then, they have a relative transverse momentum. */
           px_quark[isign] = 0.5 * pgluon.px() + (isign == 0 ? 1. : -1.) * qx;
           py_quark[isign] = 0.5 * pgluon.py() + (isign == 0 ? 1. : -1.) * qy;
           pz_quark[isign] = 0.5 * pgluon.pz();
           e_quark[isign] =
               std::sqrt(mass * mass + px_quark[isign] * px_quark[isign] +
                         py_quark[isign] * py_quark[isign] +
                         pz_quark[isign] * pz_quark[isign]);
           pquark[isign] = Pythia8::Vec4(px_quark[isign], py_quark[isign],
                                         pz_quark[isign], e_quark[isign]);
         }
 
         /* Total energy is not conserved at this point,
          * but this will be cured later. */
         pSum += pquark[0] + pquark[1] - pgluon;
         // add quark and antiquark to the event record
         event_intermediate.append(pdgid, status, col, 0, pquark[0], mass);
         event_intermediate.append(-pdgid, status, 0, acol, pquark[1], mass);
         // then remove the gluon from the record
         event_intermediate.remove(iforward, iforward);
 
         log.debug("  gluon at iforward = ", iforward, " is splitted into ",
                   pdgid, ",", -pdgid, " qqbar pair.");
         /* Increase the total number of quarks and antiquarks by 1,
          * as we have extra ones from a gluon. */
         nquark_total[iflav] += 1;
         nantiq_total[iflav] += 1;
       }
     }
   }
 
   /* The zeroth entry of event record is supposed to have the information
    * on the whole system. Specify the total momentum and invariant mass. */
   event_intermediate[0].p(pSum);
   event_intermediate[0].m(pSum.mCalc());
 
   return true;
 }

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_total	total number of quarks (antiquarks) in the system.
[out]	excess_quark	excess of quarks (antiquarks) in incoming particles, compared to the mapped ones.
[out]	excess_antiq	excess of anti-quarks (quarks) in incoming particles, compared to the mapped ones.

Definition at line 1058 of file processstring.cc.

                                                  {
   const auto &log = logger<LogArea::Pythia>();
 
   for (int iflav = 0; iflav < 5; iflav++) {
     /* Find how many constituent will be in the system after
      * changing the flavors.
      * Note that nquark_total is number of constituent right after
      * the pythia event (with mapped incoming hadrons), while the excess
      * shows how many constituents we have more or less that nquark_total. */
     int nquark_final =
         nquark_total[iflav] + excess_quark[0][iflav] + excess_quark[1][iflav];
     /* Therefore, nquark_final should not be negative.
      * negative nquark_final means that it will not be possible to
      * find a constituent to change the flavor. */
     bool enough_quark = nquark_final >= 0;
     // If that is the case, excess of constituents will be modified
     if (!enough_quark) {
       log.debug("  not enough constituents with flavor ", iflav + 1,
                 " : try to modify excess of constituents.");
       for (int ic = 0; ic < std::abs(nquark_final); ic++) {
         /* Since each incoming hadron has its own count of the excess,
          * it is necessary to find which one is problematic. */
         int ih_mod = -1;
         if (excess_quark[0][iflav] < 0) {
           ih_mod = 0;
         } else {
           ih_mod = 1;
         }
         /* Increase the excess of both quark and antiquark
          * with corresponding flavor (iflav + 1) by 1.
          * This is for conservation of the net quark number. */
         excess_quark[ih_mod][iflav] += 1;
         excess_antiq[ih_mod][iflav] += 1;
 
         /* Since incoming hadrons are mapped onto ones with
          * the same baryon number (or quark number),
          * summation of the excesses over all flavors should be zero.
          * Therefore, we need to find another flavor which has
          * a positive excess and subtract by 1. */
         for (int jflav = 0; jflav < 5; jflav++) {
           // another flavor with positive excess of constituents
           if (jflav != iflav && excess_quark[ih_mod][jflav] > 0) {
             /* Decrease the excess of both quark and antiquark
              * with corresponding flavor (jflav + 1) by 1. */
             excess_quark[ih_mod][jflav] -= 1;
             excess_antiq[ih_mod][jflav] -= 1;
             /* We only need to find one (another) flavor to subtract.
              * No more! */
             break;
           }
         }
       }
     }
   }
 }

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_intermediate	PYTHIA partonic event record to be updated according to the valence quark contents of incoming hadrons.
[out]	excess_quark	excess of quarks in incoming particles, compared to the mapped ones.
[out]	excess_antiq	excess 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 1117 of file processstring.cc.

                                                  {
   const auto &log = logger<LogArea::Pythia>();
 
   Pythia8::Vec4 pSum = event_intermediate[0].p();
   const double energy_init = pSum.e();
   log.debug("  initial total energy [GeV] : ", energy_init);
 
   // Total number of quarks and antiquarks, respectively.
   std::array<int, 5> nquark_total;
   std::array<int, 5> nantiq_total;
 
   /* Split a gluon into qqbar if we do not have enough constituents
    * to be converted in the system. */
   bool split_for_quark = splitting_gluon_qqbar(
       event_intermediate, nquark_total, nantiq_total, true, excess_quark);
   bool split_for_antiq = splitting_gluon_qqbar(
       event_intermediate, nquark_total, nantiq_total, false, excess_antiq);
 
   /* Modify excess_quark and excess_antiq if we do not have enough constituents
    * to be converted in the system. */
   if (!split_for_quark || !split_for_antiq) {
     rearrange_excess(nquark_total, excess_quark, excess_antiq);
     rearrange_excess(nantiq_total, excess_antiq, excess_quark);
   }
 
   // Final check if there are enough constituents.
   for (int iflav = 0; iflav < 5; iflav++) {
     if (nquark_total[iflav] + excess_quark[0][iflav] + excess_quark[1][iflav] <
         0) {
       log.debug("Not enough quark constituents of flavor ", iflav + 1);
       return false;
     }
 
     if (nantiq_total[iflav] + excess_antiq[0][iflav] + excess_antiq[1][iflav] <
         0) {
       log.debug("Not enough antiquark constituents of flavor ", -(iflav + 1));
       return false;
     }
   }
 
   for (int ih = 0; ih < 2; ih++) {
     log.debug("  initial excess_quark[", ih, "] = (", excess_quark[ih][0], ", ",
               excess_quark[ih][1], ", ", excess_quark[ih][2], ", ",
               excess_quark[ih][3], ", ", excess_quark[ih][4], ")");
     log.debug("  initial excess_antiq[", ih, "] = (", excess_antiq[ih][0], ", ",
               excess_antiq[ih][1], ", ", excess_antiq[ih][2], ", ",
               excess_antiq[ih][3], ", ", excess_antiq[ih][4], ")");
   }
 
   bool recovered_quarks = false;
   while (!recovered_quarks) {
     /* Flavor conversion begins with the most forward and backward parton
      * respectively for incoming_particles_[0] and incoming_particles_[1]. */
     std::array<bool, 2> find_forward = {true, false};
     const std::array<int, 5> excess_null = {0, 0, 0, 0, 0};
     std::array<int, 5> excess_total = excess_null;
 
     for (int ih = 0; ih < 2; ih++) {  // loop over incoming hadrons
       int nfrag = event_intermediate.size();
       for (int np_end = 0; np_end < nfrag - 1; np_end++) {  // constituent loop
         /* select the np_end-th most forward or backward parton and
          * change its specie.
          * np_end = 0 corresponds to the most forward,
          * np_end = 1 corresponds to the second most forward and so on. */
         int iforward =
             get_index_forward(find_forward[ih], np_end, event_intermediate);
         pSum -= event_intermediate[iforward].p();
 
         if (event_intermediate[iforward].id() > 0) {  // quark and diquark
           replace_constituent(event_intermediate[iforward], excess_quark[ih]);
           log.debug("    excess_quark[", ih, "] = (", excess_quark[ih][0], ", ",
                     excess_quark[ih][1], ", ", excess_quark[ih][2], ", ",
                     excess_quark[ih][3], ", ", excess_quark[ih][4], ")");
         } else {  // antiquark and anti-diquark
           replace_constituent(event_intermediate[iforward], excess_antiq[ih]);
           log.debug("    excess_antiq[", ih, "] = (", excess_antiq[ih][0], ", ",
                     excess_antiq[ih][1], ", ", excess_antiq[ih][2], ", ",
                     excess_antiq[ih][3], ", ", excess_antiq[ih][4], ")");
         }
 
         const int pdgid = event_intermediate[iforward].id();
         Pythia8::Vec4 pquark = event_intermediate[iforward].p();
         const double mass = pythia_hadron_->particleData.m0(pdgid);
 
         const int status = event_intermediate[iforward].status();
         const int color = event_intermediate[iforward].col();
         const int anticolor = event_intermediate[iforward].acol();
 
         pSum += pquark;
         event_intermediate.append(pdgid, status, color, anticolor, pquark,
                                   mass);
 
         event_intermediate.remove(iforward, iforward);
         /* Now the last np_end + 1 entries in event_intermediate
          * are np_end + 1 most forward (or backward) partons. */
       }
 
       // Compute the excess of net quark numbers.
       for (int j = 0; j < 5; j++) {
         excess_total[j] += (excess_quark[ih][j] - excess_antiq[ih][j]);
       }
     }
 
     /* If there is no excess of net quark numbers,
      * quark content is considered to be correct. */
     recovered_quarks = excess_total == excess_null;
   }
   log.debug("  valence quark contents of hadons are recovered.");
 
   log.debug("  current total energy [GeV] : ", pSum.e());
   /* rescale momenta of all partons by a constant factor
    * to conserve the total energy. */
   while (true) {
     if (std::abs(pSum.e() - energy_init) < really_small * energy_init) {
       break;
     }
 
     double energy_current = pSum.e();
     double slope = 0.;
     for (int i = 1; i < event_intermediate.size(); i++) {
       slope += event_intermediate[i].pAbs2() / event_intermediate[i].e();
     }
 
     const double rescale_factor = 1. + (energy_init - energy_current) / slope;
     pSum = 0.;
     for (int i = 1; i < event_intermediate.size(); i++) {
       const double px = rescale_factor * event_intermediate[i].px();
       const double py = rescale_factor * event_intermediate[i].py();
       const double pz = rescale_factor * event_intermediate[i].pz();
       const double pabs = rescale_factor * event_intermediate[i].pAbs();
       const double mass = event_intermediate[i].m();
 
       event_intermediate[i].px(px);
       event_intermediate[i].py(py);
       event_intermediate[i].pz(pz);
       event_intermediate[i].e(std::sqrt(mass * mass + pabs * pabs));
       pSum += event_intermediate[i].p();
     }
     log.debug("  parton momenta are rescaled by factor of ", rescale_factor);
   }
 
   log.debug("  final total energy [GeV] : ", pSum.e());
   /* The zeroth entry of event record is supposed to have the information
    * on the whole system. Specify the total momentum and invariant mass. */
   event_intermediate[0].p(pSum);
   event_intermediate[0].m(pSum.mCalc());
 
   return true;
 }

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_string	If it is set to be true (false), it begins with forward (backward) parton.
[out]	event_intermediate	PYTHIA event record from which a string is identified. All partons found here are removed.
[out]	event_hadronize	PYTHIA event record to which partons in a string are added.

Definition at line 1270 of file processstring.cc.

                                                                          {
   const auto &log = logger<LogArea::Pythia>();
 
   Pythia8::Vec4 pSum = 0.;
   event_hadronize.reset();
 
   // select the most forward or backward parton.
   int iforward = get_index_forward(find_forward_string, 0, event_intermediate);
   log.debug("Hard non-diff: iforward = ", iforward, "(",
             event_intermediate[iforward].id(), ")");
 
   pSum += event_intermediate[iforward].p();
   event_hadronize.append(event_intermediate[iforward]);
 
   int col_to_find = event_intermediate[iforward].acol();
   int acol_to_find = event_intermediate[iforward].col();
   event_intermediate.remove(iforward, iforward);
   log.debug("Hard non-diff: event_intermediate reduces in size to ",
             event_intermediate.size());
 
   // trace color and anti-color indices and find corresponding partons.
   while (col_to_find != 0 || acol_to_find != 0) {
     log.debug("  col_to_find = ", col_to_find,
               ", acol_to_find = ", acol_to_find);
 
     int ifound = -1;
     for (int i = 1; i < event_intermediate.size(); i++) {
       const int pdgid = event_intermediate[i].id();
       bool found_col =
           col_to_find != 0 && col_to_find == event_intermediate[i].col();
       bool found_acol =
           acol_to_find != 0 && acol_to_find == event_intermediate[i].acol();
       if (found_col) {
         log.debug("  col_to_find ", col_to_find, " from i ", i, "(", pdgid,
                   ") found");
       }
       if (found_acol) {
         log.debug("  acol_to_find ", acol_to_find, " from i ", i, "(", pdgid,
                   ") found");
       }
 
       if (found_col && !found_acol) {
         ifound = i;
         col_to_find = event_intermediate[i].acol();
         break;
       } else if (!found_col && found_acol) {
         ifound = i;
         acol_to_find = event_intermediate[i].col();
         break;
       } else if (found_col && found_acol) {
         ifound = i;
         col_to_find = 0;
         acol_to_find = 0;
         break;
       }
     }
 
     if (ifound < 0) {
       event_intermediate.list();
       event_intermediate.listJunctions();
       event_hadronize.list();
       event_hadronize.listJunctions();
       if (col_to_find != 0) {
         log.error("No parton with col = ", col_to_find);
       }
       if (acol_to_find != 0) {
         log.error("No parton with acol = ", acol_to_find);
       }
       throw std::runtime_error("Hard string could not be identified.");
     } else {
       pSum += event_intermediate[ifound].p();
       // add a parton to the new event record.
       event_hadronize.append(event_intermediate[ifound]);
       // then remove from the original event record.
       event_intermediate.remove(ifound, ifound);
       log.debug("Hard non-diff: event_intermediate reduces in size to ",
                 event_intermediate.size());
     }
   }
 
   /* The zeroth entry of event record is supposed to have the information
    * on the whole system. Specify the total momentum and invariant mass. */
   event_hadronize[0].p(pSum);
   event_hadronize[0].m(pSum.mCalc());
 }

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_string	If it is set to be true (false), it is a string in the forward (backward) direction.
[out]	event_intermediate	PYTHIA event record from which a string is identified. All partons and junction(s) found here are removed.
[out]	event_hadronize	PYTHIA 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 1358 of file processstring.cc.

                                                                            {
   const auto &log = logger<LogArea::Pythia>();
 
   event_hadronize.reset();
 
   /* Move the first junction to the event record for hadronization
    * and specify color or anti-color indices to be found.
    * If junction kind is an odd number, it connects three quarks
    * to make a color-neutral baryonic configuration.
    * Otherwise, it connects three antiquarks
    * to make a color-neutral anti-baryonic configuration. */
   const int kind = event_intermediate.kindJunction(0);
   bool sign_color = kind % 2 == 1;
   std::vector<int> col;  // color or anti-color indices of the junction legs
   for (int j = 0; j < 3; j++) {
     col.push_back(event_intermediate.colJunction(0, j));
   }
   event_hadronize.appendJunction(kind, col[0], col[1], col[2]);
   event_intermediate.eraseJunction(0);
   log.debug("junction (", col[0], ", ", col[1], ", ", col[2], ") with kind ",
             kind, " will be handled.");
 
   bool found_string = false;
   while (!found_string) {
     // trace color or anti-color indices and find corresponding partons.
     find_junction_leg(sign_color, col, event_intermediate, event_hadronize);
     found_string = true;
     for (unsigned int j = 0; j < col.size(); j++) {
       found_string = found_string && col[j] == 0;
     }
     if (!found_string) {
       /* if there is any leg which is not closed with parton,
        * look over junctions and find connected ones. */
       log.debug("  still has leg(s) unfinished.");
       sign_color = !sign_color;
       std::vector<int> junction_to_move;
       for (int i = 0; i < event_intermediate.sizeJunction(); i++) {
         const int kind_new = event_intermediate.kindJunction(i);
         /* If the original junction is associated with positive baryon number,
          * it looks for anti-junctions whose legs are connected with
          * anti-quarks (anti-colors in general). */
         if (sign_color != (kind_new % 2 == 1)) {
           continue;
         }
 
         std::array<int, 3> col_new;
         for (int k = 0; k < 3; k++) {
           col_new[k] = event_intermediate.colJunction(i, k);
         }
 
         int n_legs_connected = 0;
         // loop over remaining legs
         for (unsigned int j = 0; j < col.size(); j++) {
           if (col[j] == 0) {
             continue;
           }
           for (int k = 0; k < 3; k++) {
             if (col[j] == col_new[k]) {
               n_legs_connected += 1;
               col[j] = 0;
               col_new[k] = 0;
             }
           }
         }
 
         // specify which junction is connected to the original one.
         if (n_legs_connected > 0) {
           for (int k = 0; k < 3; k++) {
             if (col_new[k] != 0) {
               col.push_back(col_new[k]);
             }
           }
           log.debug("  junction ", i, " (",
                     event_intermediate.colJunction(i, 0), ", ",
                     event_intermediate.colJunction(i, 1), ", ",
                     event_intermediate.colJunction(i, 2), ") with kind ",
                     kind_new, " will be added.");
           junction_to_move.push_back(i);
         }
       }
 
       /* If there is any connected junction,
        * move it to the event record for hadronization. */
       for (unsigned int i = 0; i < junction_to_move.size(); i++) {
         unsigned int imove = junction_to_move[i] - i;
         const int kind_add = event_intermediate.kindJunction(imove);
         std::array<int, 3> col_add;
         for (int k = 0; k < 3; k++) {
           col_add[k] = event_intermediate.colJunction(imove, k);
         }
         // add a junction to the new event record.
         event_hadronize.appendJunction(kind_add, col_add[0], col_add[1],
                                        col_add[2]);
         // then remove from the original event record.
         event_intermediate.eraseJunction(imove);
       }
     }
   }
 
   Pythia8::Vec4 pSum = event_hadronize[0].p();
   find_forward_string = pSum.pz() > 0.;
 }

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_color	true (false) if the junction is associated with color (anti-color) indices, corresponding baryonic (anti-baryonic) string
[out]	col	set of color indices that need to be found. The value is set to be zero once the corresponding partons are found.
[out]	event_intermediate	PYTHIA event record from which a string is identified. All partons and junction(s) found here are removed.
[out]	event_hadronize	PYTHIA event record to which partons in a string are added.

See also: StringProcess::compose_string_junction(bool &, Pythia8::Event &, Pythia8::Event &)

Definition at line 1463 of file processstring.cc.

                                                                      {
   const auto &log = logger<LogArea::Pythia>();
 
   Pythia8::Vec4 pSum = event_hadronize[0].p();
   for (unsigned int j = 0; j < col.size(); j++) {
     if (col[j] == 0) {
       continue;
     }
     bool found_leg = false;
     while (!found_leg) {
       int ifound = -1;
       for (int i = 1; i < event_intermediate.size(); i++) {
         const int pdgid = event_intermediate[i].id();
         if (sign_color && col[j] == event_intermediate[i].col()) {
           log.debug("  col[", j, "] = ", col[j], " from i ", i, "(", pdgid,
                     ") found");
           ifound = i;
           col[j] = event_intermediate[i].acol();
           break;
         } else if (!sign_color && col[j] == event_intermediate[i].acol()) {
           log.debug("  acol[", j, "] = ", col[j], " from i ", i, "(", pdgid,
                     ") found");
           ifound = i;
           col[j] = event_intermediate[i].col();
           break;
         }
       }
 
       if (ifound < 0) {
         found_leg = true;
         if (event_intermediate.sizeJunction() == 0) {
           event_intermediate.list();
           event_intermediate.listJunctions();
           event_hadronize.list();
           event_hadronize.listJunctions();
           log.error("No parton with col = ", col[j],
                     " connected with junction leg ", j);
           throw std::runtime_error("Hard string could not be identified.");
         }
       } else {
         pSum += event_intermediate[ifound].p();
         // add a parton to the new event record.
         event_hadronize.append(event_intermediate[ifound]);
         // then remove from the original event record.
         event_intermediate.remove(ifound, ifound);
         log.debug("Hard non-diff: event_intermediate reduces in size to ",
                   event_intermediate.size());
         if (col[j] == 0) {
           found_leg = true;
         }
       }
     }
   }
 
   /* The zeroth entry of event record is supposed to have the information
    * on the whole system. Specify the total momentum and invariant mass. */
   event_hadronize[0].p(pSum);
   event_hadronize[0].m(pSum.mCalc());
 }

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_forward	if it looks for the most forward or backward particle.
[in]	np_end	number 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]	event	PYTHIA 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.

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

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_; }

smash::StringProcess::final_state_

ParticleList final_state_

final state array which must be accessed after the collision

Definition: processstring.h:139

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_particles	list of fragmented particles to be appended
[in]	uString	is velocity four vector of the string.
[in]	evecLong	is unit 3-vector in which string is stretched.

Returns: number of hadrons fragmented out of string.

Exceptions

std::invalid_argument	if fragmented particle is not hadron
std::invalid_argument	if string is neither mesonic nor baryonic

Definition at line 136 of file processstring.cc.

                                                                    {
   int nfrag = 0;
   int bstring = 0;
   double p_pos_tot = 0.0, p_neg_tot = 0.0;
 
   for (ParticleData &data : intermediate_particles) {
     nfrag += 1;
     bstring += data.pdgcode().baryon_number();
 
     const double pparallel = data.momentum().threevec() * evecLong;
     p_pos_tot += (data.momentum().x0() + pparallel) / sqrt2_;
     p_neg_tot += (data.momentum().x0() - pparallel) / sqrt2_;
   }
   assert(nfrag > 0);
 
   /* compute the cross section scaling factor for leading hadrons
    * based on the number of valence quarks. */
   assign_all_scaling_factors(bstring, intermediate_particles, evecLong,
                              additional_xsec_supp_);
 
   std::vector<double> xvertex_pos, xvertex_neg;
   xvertex_pos.resize(nfrag + 1);
   xvertex_neg.resize(nfrag + 1);
   // x^{+} coordinates of the forward end
   xvertex_pos[0] = p_pos_tot / kappa_tension_string_;
   // x^{-} coordinates of the backward end
   xvertex_neg[nfrag] = p_neg_tot / kappa_tension_string_;
 
   for (int i = 0; i < nfrag; i++) {
     const double pparallel =
         intermediate_particles[i].momentum().threevec() * evecLong;
 
     // recursively compute x^{+} coordinates of q-qbar formation vertex
     xvertex_pos[i + 1] =
         xvertex_pos[i] -
         (intermediate_particles[i].momentum().x0() + pparallel) /
             (kappa_tension_string_ * sqrt2_);
 
     // recursively compute x^{-} coordinates of q-qbar formation vertex
     const int ineg = nfrag - i - 1;
     xvertex_neg[ineg] =
         xvertex_neg[ineg + 1] -
         (intermediate_particles[ineg].momentum().x0() - pparallel) /
             (kappa_tension_string_ * sqrt2_);
   }
 
   // Velocity three-vector to perform Lorentz boost.
   const ThreeVector vstring = uString.velocity();
 
   /* compute the formation times of hadrons
    * from the lightcone coordinates of q-qbar formation vertices. */
   for (int i = 0; i < nfrag; i++) {
     // boost 4-momentum into the center of mass frame
     FourVector momentum =
         intermediate_particles[i].momentum().LorentzBoost(-vstring);
     intermediate_particles[i].set_4momentum(momentum);
 
     if (use_yoyo_model_) {
       // set the formation time and position in the rest frame of string
       double t_prod = (xvertex_pos[i] + xvertex_neg[i + 1]) / sqrt2_;
       FourVector fragment_position = FourVector(
           t_prod, evecLong * (xvertex_pos[i] - xvertex_neg[i + 1]) / sqrt2_);
       /* boost formation vertex into the center of mass frame
        * and then into the lab frame */
       fragment_position = fragment_position.LorentzBoost(-vstring);
       fragment_position = fragment_position.LorentzBoost(-vcomAB_);
       intermediate_particles[i].set_formation_time(
           soft_t_form_ * fragment_position.x0() + time_collision_);
     } else {
       ThreeVector v_calc = momentum.LorentzBoost(-vcomAB_).velocity();
       double gamma_factor = 1.0 / std::sqrt(1 - (v_calc).sqr());
       intermediate_particles[i].set_slow_formation_times(
           time_collision_,
           time_formation_const_ * gamma_factor + time_collision_);
     }
 
     final_state_.push_back(intermediate_particles[i]);
   }
 
   return nfrag;
 }

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]	pdgid	PDG id of particle
[in]	momentum	four-momentum of particle
[out]	intermediate_particles	particle list to which the new particle is added.

Returns: whether PDG id exists in ParticleType table.

Definition at line 748 of file processstring.h.

                                                                              {
     const std::string s = std::to_string(pdgid);
     PdgCode pythia_code(s);
     ParticleData new_particle(ParticleType::find(pythia_code));
     new_particle.set_4momentum(momentum);
     intermediate_particles.push_back(new_particle);
     return true;
   }

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_id is PDG id to be converted.

Definition at line 762 of file processstring.h.

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

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]	diquark	PDG id of diquark
[out]	q1	PDG id of quark 1
[out]	q2	PDG id of quark 2
[out]	deg_spin	spin degeneracy

Definition at line 1688 of file processstring.cc.

                                                        {
   // The 4-digit pdg id should be diquark.
   assert((std::abs(diquark) > 1000) && (std::abs(diquark) < 5510) &&
          (std::abs(diquark) % 100 < 10));
 
   // The fourth digit corresponds to the spin degeneracy.
   deg_spin = std::abs(diquark) % 10;
   // Diquark (anti-diquark) is decomposed into two quarks (antiquarks).
   const int sign_anti = diquark > 0 ? 1 : -1;
 
   // Obtain two quarks (or antiquarks) from the first and second digit.
   q1 = sign_anti * (std::abs(diquark) - (std::abs(diquark) % 1000)) / 1000;
   q2 = sign_anti * (std::abs(diquark) % 1000 - deg_spin) / 100;
 }

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]	q1	PDG code of quark 1
[in]	q2	PDG code of quark 2

Returns: PDG code of diquark composed of q1 and q2

Definition at line 1704 of file processstring.cc.

                                                      {
   assert((q1 > 0 && q2 > 0) || (q1 < 0 && q2 < 0));
   if (std::abs(q1) < std::abs(q2)) {
     std::swap(q1, q2);
   }
   int diquark = std::abs(q1 * 1000 + q2 * 100);
   /* Adding spin degeneracy = 2S+1. For identical quarks spin cannot be 0
    * because of Pauli exclusion principle, so spin 1 is assumed. Otherwise
    * S = 0 with probability 1/4 and S = 1 with probability 3/4. */
   diquark += (q1 != q2 && random::uniform_int(0, 3) == 0) ? 1 : 3;
   return (q1 < 0) ? -diquark : diquark;
 }

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_in	is PdgCode of hadron which transforms into a string.
[out]	idq1	is PDG id of quark or anti-diquark.
[out]	idq2	is PDG id of anti-quark or diquark.
[in]	xi	probability to split a nucleon into the quark it has only once and a diquark of another flavour.

Definition at line 1717 of file processstring.cc.

                                                 {
   std::array<int, 3> quarks = pdg.quark_content();
   if (pdg.is_nucleon()) {
     // protons and neutrons treated seperately since single quarks is at a
     // different position in the PDG code
     if (pdg.charge() == 0) {  // (anti)neutron
       if (random::uniform(0., 1.) < xi) {
         idq1 = quarks[0];
         idq2 = diquark_from_quarks(quarks[1], quarks[2]);
       } else {
         idq1 = quarks[1];
         idq2 = diquark_from_quarks(quarks[0], quarks[2]);
       }
     } else {  // (anti)proton
       if (random::uniform(0., 1.) < xi) {
         idq1 = quarks[2];
         idq2 = diquark_from_quarks(quarks[0], quarks[1]);
       } else {
         idq1 = quarks[0];
         idq2 = diquark_from_quarks(quarks[1], quarks[2]);
       }
     }
   } else {
     if (pdg.is_meson()) {
       idq1 = quarks[1];
       idq2 = quarks[2];
       /* Some mesons with PDG id 11X are actually mixed state of uubar and
        * ddbar. have a random selection whether we have uubar or ddbar in this
        * case. */
       if (idq1 == 1 && idq2 == -1 && random::uniform_int(0, 1) == 0) {
         idq1 = 2;
         idq2 = -2;
       }
     } else {
       assert(pdg.is_baryon());
       // Get random quark to position 0
       std::swap(quarks[random::uniform_int(0, 2)], quarks[0]);
       idq1 = quarks[0];
       idq2 = diquark_from_quarks(quarks[1], quarks[2]);
     }
   }
   // Fulfil the convention: idq1 should be quark or anti-diquark
   if (idq1 < 0) {
     std::swap(idq1, idq2);
   }
 }

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]	energy	time component
[in]	mom	spatial three-vector

Returns: Pythia Vec4 from energy and ThreeVector

Definition at line 800 of file processstring.h.

                                                                      {
     return Pythia8::Vec4(mom.x1(), mom.x2(), mom.x3(), energy);
   }

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]	particle	particle object from PYTHIA event generation where z-axis is set to be the collision axis
[in]	evec_basis	three 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.

                                                                    {
     ThreeVector three_momentum = evec_basis[0] * particle.pz() +
                                  evec_basis[1] * particle.px() +
                                  evec_basis[2] * particle.py();
     return FourVector(particle.e(), three_momentum);
   }

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]	idq1	PDG id of quark or anti-diquark (carrying color index).
[in]	idq2	PDG id of diquark or anti-quark (carrying anti-color index).
[in]	mString	the string mass. [GeV]
[out]	evecLong	unit 3-vector specifying the direction of diquark or anti-diquark.
[in]	flip_string_ends	is whether or not we randomly switch string ends.
[in]	separate_fragment_baryon	is whether fragment leading baryon (or anti-baryon) with separate fragmentation function.
[out]	intermediate_particles	list of fragmented particles

Returns: number of hadrons fragmented out of string.

Exceptions

std::runtime_error if string mass is lower than threshold set by PYTHIA

Definition at line 1765 of file processstring.cc.

                                                                          {
   const auto &log = logger<LogArea::Pythia>();
   pythia_hadron_->event.reset();
   intermediate_particles.clear();
 
   log.debug("initial quark content for fragment_string : ", idq1, ", ", idq2);
   // PDG id of quark constituents of string ends
   std::array<int, 2> idqIn;
   idqIn[0] = idq1;
   idqIn[1] = idq2;
 
   int bstring = 0;
   // constituent masses of the string
   std::array<double, 2> m_const;
 
   for (int i = 0; i < 2; i++) {
     // evaluate total baryon number of the string times 3
     bstring += pythia_hadron_->particleData.baryonNumberType(idqIn[i]);
 
     m_const[i] = pythia_hadron_->particleData.m0(idqIn[i]);
   }
   log.debug("baryon number of string times 3 : ", bstring);
 
   if (flip_string_ends && random::uniform_int(0, 1) == 0) {
     /* in the case where we flip the string ends,
      * we need to flip the longitudinal unit vector itself
      * since it is set to be direction of diquark (anti-quark)
      * or anti-diquark. */
     evecLong = -evecLong;
   }
 
   if (m_const[0] + m_const[1] > mString) {
     throw std::runtime_error("String fragmentation: m1 + m2 > mString");
   }
   Pythia8::Vec4 pSum = 0.;
 
   int number_of_fragments = 0;
 
   if (separate_fragment_baryon && (std::abs(bstring) == 3) &&
       (mString > m_const[0] + m_const[1] + 1.)) {
     /* In the case of separate fragmentation function for leading baryons,
      * the leading baryon (antibaryon) is fragmented from the original string
      * with a gaussian fragmentation function.
      * It is then followed by fragmentation of the remaining mesonic string. */
     const double ssbar_supp = pythia_hadron_->parm("StringFlav:probStoUD");
     const double sigma_qt_frag = pythia_hadron_->parm("StringPT:sigma");
     std::array<ThreeVector, 3> evec_basis;
     make_orthonormal_basis(evecLong, evec_basis);
 
     /* transverse momentum (and magnitude) acquired by the qqbar pair
      * created in the middle of string */
     double QTrx, QTry, QTrn;
     /* lightcone momenta p^+ and p^- of the remaining (mesonic) string
      * p^{\pm} is defined as (E \pm p_{longitudinal}) / sqrts{2}. */
     double ppos_string_new, pneg_string_new;
     // transverse mass of the remaining (mesonic) string
     double mTrn_string_new;
     // transverse masses of quark constituents of string ends
     std::array<double, 2> m_trans;
 
     const int niter_max = 10000;
     int iiter = 0;
     bool found_leading_baryon = false;
     while (!found_leading_baryon) {
       int idnew_qqbar = 0;
       /* Determine flavor of the new qqbar pair
        * based on the strangeness suppression factor specifed in config. */
       if (random::uniform(0., 1. + ssbar_supp) < 1.) {
         if (random::uniform_int(0, 1) == 0) {
           idnew_qqbar = 1;  // ddbar pair
         } else {
           idnew_qqbar = 2;  // uubar pair
         }
       } else {
         idnew_qqbar = 3;  // ssbar pair
       }
 
       int id_diquark = 0;
       if (bstring > 0) {
         /* q(idq1) + diquark(idq2) baryonic string fragments into
          * q(idq1) + qbar(-idnew_qqbar) mesonic string and
          * q(idnew_qqbar) + diquark(idq2) baryon. */
         id_diquark = std::abs(idq2);
         idqIn[1] = -idnew_qqbar;
       } else {
         /* anti-diquark(idq1) + qbar(idq2) anti-baryonic string fragments into
          * q(idnew_qqbar) + qbar(idq2) mesonic string and
          * anti-diquark(idq1) + qbar(-idnew_qqbar) antibaryon. */
         id_diquark = std::abs(idq1);
         idqIn[0] = idnew_qqbar;
       }
       log.debug("quark constituents for leading baryon : ", idnew_qqbar, ", ",
                 id_diquark);
 
       // net quark number of d, u, s, c and b flavors
       std::array<int, 5> frag_net_q;
       /* Evaluate total net quark number of baryon (antibaryon)
        * from the valence quark constituents. */
       for (int iq = 0; iq < 5; iq++) {
         frag_net_q[iq] =
             (bstring > 0 ? 1 : -1) *
             (pythia_hadron_->particleData.nQuarksInCode(idnew_qqbar, iq + 1) +
              pythia_hadron_->particleData.nQuarksInCode(id_diquark, iq + 1));
       }
       const int frag_iso3 = frag_net_q[1] - frag_net_q[0];
       const int frag_strange = -frag_net_q[2];
       const int frag_charm = frag_net_q[3];
       const int frag_bottom = -frag_net_q[4];
       log.debug("  conserved charges of leading baryon : iso3 = ", frag_iso3,
                 ", strangeness = ", frag_strange, ", charmness = ", frag_charm,
                 ", bottomness = ", frag_bottom);
 
       std::vector<int> pdgid_possible;
       std::vector<double> weight_possible;
       std::vector<double> weight_summed;
       /* loop over hadronic species
        * Any hadron with the same valence quark contents is allowed and
        * the probability goes like spin degeneracy over mass. */
       for (auto &ptype : ParticleType::list_all()) {
         const int pdgid =
             (bstring > 0 ? 1 : -1) * std::abs(ptype.pdgcode().get_decimal());
         if ((pythia_hadron_->particleData.isParticle(pdgid)) &&
             (bstring == 3 * ptype.pdgcode().baryon_number()) &&
             (frag_iso3 == ptype.pdgcode().isospin3()) &&
             (frag_strange == ptype.pdgcode().strangeness()) &&
             (frag_charm == ptype.pdgcode().charmness()) &&
             (frag_bottom == ptype.pdgcode().bottomness())) {
           const int spin_degeneracy = ptype.pdgcode().spin_degeneracy();
           const double mass_pole = ptype.mass();
           const double weight =
               static_cast<double>(spin_degeneracy) / mass_pole;
           pdgid_possible.push_back(pdgid);
           weight_possible.push_back(weight);
 
           log.debug("  PDG id ", pdgid, " is possible with weight ", weight);
         }
       }
       const int n_possible = pdgid_possible.size();
       weight_summed.push_back(0.);
       for (int i = 0; i < n_possible; i++) {
         weight_summed.push_back(weight_summed[i] + weight_possible[i]);
       }
 
       /* Sample baryon (antibaryon) specie,
        * which is fragmented from the leading diquark (anti-diquark). */
       int pdgid_frag = 0;
       const double uspc = random::uniform(0., weight_summed[n_possible]);
       for (int i = 0; i < n_possible; i++) {
         if ((uspc >= weight_summed[i]) && (uspc < weight_summed[i + 1])) {
           pdgid_frag = pdgid_possible[i];
           log.debug("PDG id ", pdgid_frag, " selected for leading baryon.");
           break;
         }
       }
       // Sample mass.
       const double mass_frag = pythia_hadron_->particleData.mSel(pdgid_frag);
 
       // Sample transverse momentum carried by baryon (antibaryon).
       QTrx = random::normal(0., sigma_qt_frag / sqrt2_);
       QTry = random::normal(0., sigma_qt_frag / sqrt2_);
       log.debug("Transverse momentum (", QTrx, ", ", QTry,
                 ") GeV selected for the leading baryon.");
       QTrn = std::sqrt(QTrx * QTrx + QTry * QTry);
       // transverse mass of the fragmented leading hadron
       const double mTrn_frag = std::sqrt(QTrn * QTrn + mass_frag * mass_frag);
 
       /* Sample lightcone momentum fraction carried by the baryon (antibaryon).
        * The corresponding fragmentation function has gaussian shape
        * i.e., exp( - (x - leading_frag_mean_)^2 / (2 leading_frag_width_^2) ).
        * This is done using the rejection method with an Lorentzian
        * envelop function, which is
        * 1 / (1 + (x - leading_frag_mean_)^2 / (4 leading_frag_width_^2) ). */
       double xfrac = 0.;
       bool xfrac_accepted = false;
       while (!xfrac_accepted) {
         const double angle =
             random::uniform(0., 1.) *
                 (std::atan(0.5 * (1. - leading_frag_mean_) /
                            leading_frag_width_) +
                  std::atan(0.5 * leading_frag_mean_ / leading_frag_width_)) -
             std::atan(0.5 * leading_frag_mean_ / leading_frag_width_);
         // lightcone momentum fraction sampled from the envelop function
         xfrac = leading_frag_mean_ + 2. * leading_frag_width_ * std::tan(angle);
 
         /* The rejection method is applied to sample
          * according to the real fragmentation function. */
         const double xf_tmp =
             std::abs(xfrac - leading_frag_mean_) / leading_frag_width_;
         const double xf_env = (1. + really_small) / (1. + xf_tmp * xf_tmp / 4.);
         const double xf_pdf = std::exp(-xf_tmp * xf_tmp / 2.);
         if (random::uniform(0., xf_env) < xf_pdf && xfrac > 0. && xfrac < 1.) {
           xfrac_accepted = true;
         }
       }
 
       /* Determine kinematics of the fragmented baryon (antibaryon) and
        * remaining (mesonic) string.
        * It begins with the lightcone momentum, p^+ (ppos) and p^- (pneg). */
       const double ppos_frag = xfrac * mString / sqrt2_;
       const double pneg_frag = 0.5 * mTrn_frag * mTrn_frag / ppos_frag;
       ppos_string_new = mString / sqrt2_ - ppos_frag;
       pneg_string_new = mString / sqrt2_ - pneg_frag;
       mTrn_string_new =
           std::sqrt(std::max(0., 2. * ppos_string_new * pneg_string_new));
 
       /* Quark (antiquark) constituents of the remaining string are
        * different from those of the original string.
        * Therefore, the constituent masses have to be updated. */
       for (int i = 0; i < 2; i++) {
         m_const[i] = pythia_hadron_->particleData.m0(idqIn[i]);
       }
       if (bstring > 0) {  // in the case of baryonic string
         /* Quark is coming from the original string
          * and therefore has zero transverse momentum. */
         m_trans[0] = m_const[0];
         /* Antiquark is coming from the newly produced qqbar pair
          * and therefore has transverse momentum, which is opposite to
          * that of the fragmented (leading) antibaryon. */
         m_trans[1] = std::sqrt(m_const[1] * m_const[1] + QTrn * QTrn);
       } else {  // in the case of anti-baryonic string
         /* Quark is coming from the newly produced qqbar pair
          * and therefore has transverse momentum, which is opposite to
          * that of the fragmented (leading) baryon. */
         m_trans[0] = std::sqrt(m_const[0] * m_const[0] + QTrn * QTrn);
         /* Antiquark is coming from the original string
          * and therefore has zero transverse momentum. */
         m_trans[1] = m_const[1];
       }
 
       if (mTrn_string_new > m_trans[0] + m_trans[1]) {
         found_leading_baryon = true;
 
         FourVector mom_frag((ppos_frag + pneg_frag) / sqrt2_,
                             evec_basis[0] * (ppos_frag - pneg_frag) / sqrt2_ +
                                 evec_basis[1] * QTrx + evec_basis[2] * QTry);
         log.debug("appending the leading baryon ", pdgid_frag,
                   " to the intermediate particle list.");
         /* If the remaining string has enough transverse mass,
          * add fragmented baryon (antibaryon) to the intermediate list. */
         bool found_ptype = append_intermediate_list(pdgid_frag, mom_frag,
                                                     intermediate_particles);
         if (!found_ptype) {
           log.error("PDG ID ", pdgid_frag, " should exist in ParticleType.");
           throw std::runtime_error("string fragmentation failed.");
         }
         number_of_fragments++;
         log.debug("proceed to the next step");
       }
 
       if (iiter == niter_max) {
         return 0;
       }
       iiter += 1;
     }
 
     // lightcone momentum p^+ of the quark constituents on the string ends
     std::array<double, 2> ppos_parton;
     // lightcone momentum p^- of the quark constituents on the string ends
     std::array<double, 2> pneg_parton;
 
     /* lightcone momenta of the string ends (quark and antiquark)
      * First, obtain ppos_parton[0] and ppos_parton[1]
      * (p^+ of quark and antiquark) by solving the following equations.
      * ppos_string_new = ppos_parton[0] + ppos_parton[1]
      * pneg_string_new = 0.5 * m_trans[0] * m_trans[0] / ppos_parton[0] +
      *                   0.5 * m_trans[1] * m_trans[1] / ppos_parton[1] */
     const double pb_const =
         (mTrn_string_new * mTrn_string_new + m_trans[0] * m_trans[0] -
          m_trans[1] * m_trans[1]) /
         (4. * pneg_string_new);
     const double pc_const =
         0.5 * m_trans[0] * m_trans[0] * ppos_string_new / pneg_string_new;
     ppos_parton[0] = pb_const + (bstring > 0 ? -1. : 1.) *
                                     std::sqrt(pb_const * pb_const - pc_const);
     ppos_parton[1] = ppos_string_new - ppos_parton[0];
     /* Then, compute pneg_parton[0] and pneg_parton[1]
      * (p^- of quark and antiquark) from the dispersion relation.
      * 2 p^+ p^- = m_transverse^2 */
     for (int i = 0; i < 2; i++) {
       pneg_parton[i] = 0.5 * m_trans[i] * m_trans[i] / ppos_parton[i];
     }
 
     const int status = 1;
     int color, anticolor;
     ThreeVector three_mom;
     Pythia8::Vec4 pquark;
 
     // quark end of the remaining (mesonic) string
     color = 1;
     anticolor = 0;
     if (bstring > 0) {  // in the case of baryonic string
       /* Quark is coming from the original string
        * and therefore has zero transverse momentum. */
       three_mom = evec_basis[0] * (ppos_parton[0] - pneg_parton[0]) / sqrt2_;
     } else {  // in the case of anti-baryonic string
       /* Quark is coming from the newly produced qqbar pair
        * and therefore has transverse momentum, which is opposite to
        * that of the fragmented (leading) baryon. */
       three_mom = evec_basis[0] * (ppos_parton[0] - pneg_parton[0]) / sqrt2_ -
                   evec_basis[1] * QTrx - evec_basis[2] * QTry;
     }
     pquark = set_Vec4((ppos_parton[0] + pneg_parton[0]) / sqrt2_, three_mom);
     pSum += pquark;
     pythia_hadron_->event.append(idqIn[0], status, color, anticolor, pquark,
                                  m_const[0]);
 
     // antiquark end of the remaining (mesonic) string
     color = 0;
     anticolor = 1;
     if (bstring > 0) {  // in the case of baryonic string
       /* Antiquark is coming from the newly produced qqbar pair
        * and therefore has transverse momentum, which is opposite to
        * that of the fragmented (leading) antibaryon. */
       three_mom = evec_basis[0] * (ppos_parton[1] - pneg_parton[1]) / sqrt2_ -
                   evec_basis[1] * QTrx - evec_basis[2] * QTry;
     } else {  // in the case of anti-baryonic string
       /* Antiquark is coming from the original string
        * and therefore has zero transverse momentum. */
       three_mom = evec_basis[0] * (ppos_parton[1] - pneg_parton[1]) / sqrt2_;
     }
     pquark = set_Vec4((ppos_parton[1] + pneg_parton[1]) / sqrt2_, three_mom);
     pSum += pquark;
     pythia_hadron_->event.append(idqIn[1], status, color, anticolor, pquark,
                                  m_const[1]);
   } else {
     /* diquark (anti-quark) with PDG id idq2 is going in the direction of
      * evecLong.
      * quark with PDG id idq1 is going in the direction opposite to evecLong. */
     double sign_direction = 1.;
     if (bstring == -3) {  // anti-baryonic string
       /* anti-diquark with PDG id idq1 is going in the direction of evecLong.
        * anti-quark with PDG id idq2 is going in the direction
        * opposite to evecLong. */
       sign_direction = -1;
     }
 
     // evaluate momenta of quarks
     const double pCMquark = pCM(mString, m_const[0], m_const[1]);
     const double E1 = std::sqrt(m_const[0] * m_const[0] + pCMquark * pCMquark);
     const double E2 = std::sqrt(m_const[1] * m_const[1] + pCMquark * pCMquark);
 
     ThreeVector direction = sign_direction * evecLong;
 
     // For status and (anti)color see \iref{Sjostrand:2007gs}.
     const int status1 = 1, color1 = 1, anticolor1 = 0;
     Pythia8::Vec4 pquark = set_Vec4(E1, -direction * pCMquark);
     pSum += pquark;
     pythia_hadron_->event.append(idqIn[0], status1, color1, anticolor1, pquark,
                                  m_const[0]);
 
     const int status2 = 1, color2 = 0, anticolor2 = 1;
     pquark = set_Vec4(E2, direction * pCMquark);
     pSum += pquark;
     pythia_hadron_->event.append(idqIn[1], status2, color2, anticolor2, pquark,
                                  m_const[1]);
   }
 
   log.debug("fragmenting a string with ", idqIn[0], ", ", idqIn[1]);
   // implement PYTHIA fragmentation
   pythia_hadron_->event[0].p(pSum);
   pythia_hadron_->event[0].m(pSum.mCalc());
   const bool successful_hadronization = pythia_hadron_->next();
   if (successful_hadronization) {
     for (int ipyth = 0; ipyth < pythia_hadron_->event.size(); ipyth++) {
       if (!pythia_hadron_->event[ipyth].isFinal()) {
         continue;
       }
       int pythia_id = pythia_hadron_->event[ipyth].id();
       /* K_short and K_long need are converted to K0
        * since SMASH only knows K0 */
       convert_KaonLS(pythia_id);
       FourVector momentum(
           pythia_hadron_->event[ipyth].e(), pythia_hadron_->event[ipyth].px(),
           pythia_hadron_->event[ipyth].py(), pythia_hadron_->event[ipyth].pz());
       log.debug("appending the fragmented hadron ", pythia_id,
                 " to the intermediate particle list.");
       bool found_ptype =
           append_intermediate_list(pythia_id, momentum, intermediate_particles);
       if (!found_ptype) {
         log.warn("PDG ID ", pythia_id,
                  " does not exist in ParticleType - start over.");
         intermediate_particles.clear();
         return 0;
       }
 
       number_of_fragments++;
     }
 
     return number_of_fragments;
   } else {
     return 0;
   }
 }

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_string	baryon number of the string
[out]	outgoing_particles	list of string fragments to which scaling factors are assigned
[in]	evecLong	direction in which the string is stretched
[in]	suppression_factor	additional coherence factor to be multiplied with scaling factor

Definition at line 2194 of file processstring.cc.

                                                                           {
   // Set each particle's cross section scaling factor to 0 first
   for (ParticleData &data : outgoing_particles) {
     data.set_cross_section_scaling_factor(0.0);
   }
   // sort outgoing particles according to the longitudinal velocity
   std::sort(outgoing_particles.begin(), outgoing_particles.end(),
             [&](ParticleData i, ParticleData j) {
               return i.momentum().velocity() * evecLong >
                      j.momentum().velocity() * evecLong;
             });
   int nq1, nq2;  // number of quarks at both ends of the string
   switch (baryon_string) {
     case 0:
       nq1 = -1;
       nq2 = 1;
       break;
     case 1:
       nq1 = 2;
       nq2 = 1;
       break;
     case -1:
       nq1 = -2;
       nq2 = -1;
       break;
     default:
       throw std::runtime_error("string is neither mesonic nor baryonic");
   }
   // Try to find nq1 on one string end and nq2 on the other string end and the
   // other way around. When the leading particles are close to the string ends,
   // the quarks are assumed to be distributed this way.
   std::pair<int, int> i = find_leading(nq1, nq2, outgoing_particles);
   std::pair<int, int> j = find_leading(nq2, nq1, outgoing_particles);
   if (baryon_string == 0 && i.second - i.first < j.second - j.first) {
     assign_scaling_factor(nq2, outgoing_particles[j.first], suppression_factor);
     assign_scaling_factor(nq1, outgoing_particles[j.second],
                           suppression_factor);
   } else {
     assign_scaling_factor(nq1, outgoing_particles[i.first], suppression_factor);
     assign_scaling_factor(nq2, outgoing_particles[i.second],
                           suppression_factor);
   }
 }

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]	nq1	number of valence quarks from excited hadron at forward end of the string
[in]	nq2	number of valance quarks from excited hadron at backward end of the string
[in]	list	list of string fragments

Returns: indices of the leading hadrons in list

Definition at line 2177 of file processstring.cc.

                                                                     {
   assert(list.size() >= 2);
   int end = list.size() - 1;
   int i1, i2;
   for (i1 = 0;
        i1 <= end && !list[i1].pdgcode().contains_enough_valence_quarks(nq1);
        i1++) {
   }
   for (i2 = end;
        i2 >= 0 && !list[i2].pdgcode().contains_enough_valence_quarks(nq2);
        i2--) {
   }
   std::pair<int, int> indices(i1, i2);
   return indices;
 }

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]	nquark	number of valence quarks from the excited hadron contained in the given string fragment `data`
[out]	data	particle to assign a scaling factor to
[in]	suppression_factor	coherence factor to decrease scaling factor

Definition at line 2162 of file processstring.cc.

                                                                      {
   int nbaryon = data.pdgcode().baryon_number();
   if (nbaryon == 0) {
     // Mesons always get a scaling factor of 1/2 since there is never
     // a q-qbar pair at the end of a string so nquark is always 1
     data.set_cross_section_scaling_factor(0.5 * suppression_factor);
   } else if (data.is_baryon()) {
     // Leading baryons get a factor of 2/3 if they carry 2
     // and 1/3 if they carry 1 of the strings valence quarks
     data.set_cross_section_scaling_factor(suppression_factor * nquark /
                                           (3.0 * nbaryon));
   }
 }

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] pdg PdgCode that will be mapped

Returns: mapped PDG id to be used in PYTHIA

Exceptions

std::runtime_error if the incoming particle is neither hadron nor lepton.

Definition at line 2241 of file processstring.cc.

                                                   {
   PdgCode pdg_mapped(0x0);
 
   if (pdg.baryon_number() == 1) {  // baryon
     pdg_mapped = pdg.charge() > 0 ? PdgCode(pdg::p) : PdgCode(pdg::n);
   } else if (pdg.baryon_number() == -1) {  // antibaryon
     pdg_mapped = pdg.charge() < 0 ? PdgCode(-pdg::p) : PdgCode(-pdg::n);
   } else if (pdg.is_hadron()) {  // meson
     if (pdg.charge() >= 0) {
       pdg_mapped = PdgCode(pdg::pi_p);
     } else {
       pdg_mapped = PdgCode(pdg::pi_m);
     }
   } else if (pdg.is_lepton()) {  // lepton
     pdg_mapped = pdg.charge() < 0 ? PdgCode(0x11) : PdgCode(-0x11);
   } else {
     throw std::runtime_error("StringProcess::pdg_map_for_pythia failed.");
   }
 
   return pdg_mapped.get_decimal();
 }

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:

src/include/smash/processstring.h
src/processstring.cc

Public Member Functions

Static Public Member Functions

Private Attributes

Constructor & Destructor Documentation

◆ StringProcess()

Member Function Documentation

◆ common_setup_pythia()

◆ init_pythia_hadron_rndm()

◆ get_ptr_pythia_parton()

◆ cross_sections_diffractive()

◆ set_pmin_gluon_lightcone()

◆ set_pow_fgluon()

◆ set_pow_fquark()

◆ set_sigma_qperp_()

◆ set_tension_string()

◆ init()

◆ make_orthonormal_basis()

◆ compute_incoming_lightcone_momenta()

◆ set_mass_and_direction_2strings()

◆ make_final_state_2strings()

◆ next_SDiff()

◆ next_DDiff()

◆ next_NDiffSoft()

◆ next_NDiffHard()

◆ next_BBbarAnn()

◆ find_excess_constituent()

◆ replace_constituent()

◆ find_total_number_constituent()

◆ splitting_gluon_qqbar()

◆ rearrange_excess()

◆ restore_constituent()

◆ compose_string_parton()

◆ compose_string_junction()

◆ find_junction_leg()

◆ get_index_forward()

◆ get_final_state()

◆ append_final_state()

◆ append_intermediate_list()

◆ convert_KaonLS()

◆ quarks_from_diquark()

◆ diquark_from_quarks()

◆ make_string_ends()

◆ set_Vec4()

◆ reorient()

◆ fragment_string()

◆ assign_all_scaling_factors()

◆ find_leading()

◆ assign_scaling_factor()

◆ pdg_map_for_pythia()

Member Data Documentation

◆ PPosA_

◆ PPosB_

◆ PNegA_

◆ PNegB_

◆ massA_

◆ massB_

◆ sqrtsAB_

◆ PDGcodes_

◆ plab_

◆ pcom_

◆ ucomAB_

◆ vcomAB_

◆ evecBasisAB_

◆ NpartFinal_

◆ NpartString_

◆ pmin_gluon_lightcone_

◆ pow_fgluon_beta_

◆ pow_fquark_alpha_

◆ pow_fquark_beta_

◆ sigma_qperp_

◆ leading_frag_mean_

◆ leading_frag_width_

◆ kappa_tension_string_

◆ additional_xsec_supp_

◆ time_formation_const_

◆ soft_t_form_

◆ time_collision_

◆ use_yoyo_model_

◆ prob_proton_to_d_uu_

◆ sqrt2_