smash::StringProcess Class Reference

#include <stringprocess.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 [8], Bleicher:1999xi [13] and subsequent fragmentation according to the LUND/PYTHIA fragmentation scheme Andersson:1983ia [5], Sjostrand:2014zea [60], Bierlich:2022pfr [11].

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 stringprocess.h.

Collaboration diagram for smash::StringProcess:

Public Types
enum class	LeadingStatus : int {   LeadingQuark = 202 , LeadingDiquark = 63 , NonLeadingParton = 23 , FromLeadingQuark = 203 ,   FromLeadingDiquark = 204 }
	PYTHIA status codes used to track leading/valence ancestry. More...

Public Member Functions
	StringProcess (Configuration &config)
	Constructor, initializes PYTHIA. More...
void	common_setup_pythia (Pythia8::Pythia *pythia_in, double strange_supp, double diquark_supp, double popcorn_rate, double stringz_a, double stringz_b, double string_sigma_T)
	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...
void	tag_leading_hadrons (Pythia8::Event &event)
	Tag leading hadrons in a hadronized string. More...
bool	next (ProcessType type)
	Generate the next string process for a given process type. More...
bool	string_above_threshold (const Pythia8::Event &event)
	Check whether all strings in a PYTHIA event are above fragmentation threshold. More...
double	estimate_string_threshold (int p_left, int p_right)
	Estimate the minimum invariant mass required for a string to fragment. More...
void	set_color_by_type (Pythia8::Particle &p, int color)
	Set the color or anticolor index of a particle according to its type. More...
std::optional< ParticleList >	hadronize (const Pythia8::Event &string_evt)
	Hadronize a single partonic string configuration using Pythia8 and convert the produced hadrons into SMASH particles. 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 [56]. 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 form_intermediate_particles. More...
void	set_mpi_initialization_sqrts (double sqrts)
	Set the center-of-mass energy used to initialize MPI-capable Pythia objects. 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...
void	replace_constituent (Pythia8::Particle &particle, std::array< int, 5 > &excess_constituent)
	Convert a partonic PYTHIA particle 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 constituents 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...
void	form_intermediate_particles (ParticleList &intermediate_particles, const FourVector &pString, const ThreeVector &evecLong, double additional_xsec_supp=1.0, bool find_and_scale_leading=true)
	Set formation times and cross-section scaling factors for fragmented hadrons as described in Andersson:1983ia [5]. More...
double	getPPosA ()
double	getPNegA ()
double	getPPosB ()
double	getPnegB ()
double	get_massA ()
double	get_massB ()
double	get_sqrts ()
std::array< PdgCode, 2 >	get_PDGs ()
std::array< FourVector, 2 >	get_plab ()
std::array< FourVector, 2 >	get_pcom ()
FourVector	get_ucom ()
ThreeVector	get_vcom ()
double	get_tcoll ()
ParticleList	get_final_state ()

Static Public Member Functions
static FourVector	make_smash_4vec (const Pythia8::Vec4 &p)
	Convert a PYTHIA four-vector into a SMASH four-vector. More...
static Pythia8::Vec4	make_pythia_4vec (const FourVector &p)
	Convert a SMASH four-vector into a PYTHIA four-vector. More...
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 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...

Public Attributes
Pythia8::RotBstMatrix	to_cm_
	Rotation/boost matrix to transform particles to the center-of-mass frame. More...

Private Types
typedef std::map< std::pair< int, int >, std::unique_ptr< Pythia8::Pythia > >	pythia_map
	Map containing PYTHIA objects for hard string routines. More...

Private Member Functions
std::vector< bool >	compute_beam_valence_flags (Pythia8::Pythia &pythia)
	Compute flags identifying beam valence partons (quarks or diquarks) that act as leading partons after the initial interaction. More...
int	leading_hadron_status_from_endpoint (const Pythia8::Particle &end)
	Determine the custom leading-hadron status code from a string endpoint. More...
bool	is_leading_parton (const Pythia8::Particle &p)
	Check whether a particle is tagged as a leading parton. More...
bool	is_leading_from_quark (const Pythia8::Particle &p)
	Check whether a particle is tagged as originating from a leading quark. More...
bool	is_leading_from_diquark (const Pythia8::Particle &p)
	Check whether a particle is tagged as originating from a leading diquark. More...
bool	is_leading (const Pythia8::Particle &p)
	Check whether a particle is tagged as originating from a leading endpoint. More...
bool	append_string (const Pythia8::Vec4 &p_str, const std::array< int, 2 > &ends, int color_tag, bool use_projectile_axis, bool random_flip_of_endpoints=false)
	Append a single two-endpoint string as an independent PYTHIA event. More...
bool	next_SDiff (bool is_AB_to_AX)
	Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns [32]. 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 [8], Bleicher:1999xi [13]. More...
bool	next_NDiffSoft ()
	Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig [18]. More...
bool	next_Hard (ProcessType type)
	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 [8], Bleicher:1999xi [13] does, it create two mesonic strings after annihilating one quark-antiquark pair. 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...
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	stringz_a_leading_
	parameter (StringZ:aLund) for the fragmentation function of leading baryon in soft non-diffractive string processes More...
double	stringz_b_leading_
	parameter (StringZ:bLund) for the fragmentation function of leading baryon in soft non-diffractive string processes More...
double	stringz_a_produce_
	parameter (StringZ:aLund) for the fragmentation function of other (produced) hadrons in soft non-diffractive string processes More...
double	stringz_b_produce_
	parameter (StringZ:bLund) for the fragmentation function of other (produced) hadrons in soft non-diffractive string processes More...
double	strange_supp_
	strange quark suppression factor More...
double	diquark_supp_
	diquark suppression factor More...
double	popcorn_rate_
	popcorn rate More...
double	damp_popcorn_
	damp popcorn meson from diquark remnant endpoint rate More...
double	string_sigma_T_
	transverse momentum spread in string fragmentation More...
double	kappa_tension_string_
	string tension [GeV/fm] 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	mass_dependent_formation_times_
	Whether the formation time should depend on the mass of the fragment according to Andersson:1983ia [5] eq. 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...
bool	separate_fragment_baryon_
	Whether to use a separate fragmentation function for leading baryons. More...
bool	use_monash_tune_
	Whether to use the monash tune Skands:2014pea [61] for all string processes. More...
double	additional_xsec_supp_
	additional cross-section suppression factor to take coherence effect into account. More...
std::optional< double >	mpi_initialization_sqrts_
	Optional center-of-mass energy used to initialize MPI-capable Pythia objects. More...
std::vector< std::string >	pythia_settings_
	Additional Pythia 8 settings passed to each internal Pythia instance. More...
ParticleList	final_state_
	final state array which must be accessed after the collision More...
pythia_map	hard_map_
	Map object to contain the different pythia objects. 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::StringFlav	pythia_stringflav_
	An object for the flavor selection in string fragmentation in the case of separate fragmentation function for leading baryon. More...
Pythia8::Event	event_intermediate_
	event record for intermediate partonic state in the hard string routine More...
std::vector< Pythia8::Event >	string_parton_events_
	PYTHIA event records containing string partons to be hadronized. More...

Member Typedef Documentation

pythia_map

typedef std::map<std::pair<int, int>, std::unique_ptr<Pythia8::Pythia> > smash::StringProcess::pythia_map

private

Map containing PYTHIA objects for hard string routines.

Particle IDs are used as the keys to obtain the respective object. This was introduced to reduce the amount of Pythia init() calls.

Definition at line 288 of file stringprocess.h.

Member Enumeration Documentation

LeadingStatus

enum smash::StringProcess::LeadingStatus : int

strong

PYTHIA status codes used to track leading/valence ancestry.

PYTHIA uses integer status codes to classify particles in the event record. SMASH uses a combination of standard PYTHIA codes and custom internal codes to tag leading partons and hadrons originating from leading quark/diquark endpoints.

Custom values are chosen to avoid collisions with commonly used PYTHIA internal codes. Standard PYTHIA values are used deliberately where they activate specific hadronization machinery.

Enumerator
LeadingQuark	Custom status assigned to leading (valence) quarks.
LeadingDiquark	Standard PYTHIA beam-remnant status used for leading diquarks. Leading diquarks are assigned status 63 so that PYTHIA treats them as beam remnants during hadronization. This allows PYTHIA's beam remnant machinery to be used for leading baryon production, including popcorn suppression and the optional hard-remnant-baryon treatment. In particular, when BeamRemnants:hardRemnantBaryon is enabled, PYTHIA replaces the standard Lund symmetric fragmentation function for the leading baryon by a dedicated remnant-baryon fragmentation function controlled by BeamRemnants:aRemnantBaryon and BeamRemnants:bRemnantBaryon. This produces harder leading baryons than ordinary string fragmentation. Therefore status 63 is intentionally used here rather than a custom status code.
NonLeadingParton	Standard PYTHIA status for non-leading partons that should be hadronized. The PYTHIA manual recommends status 23 for particles provided as input to standalone hadronization. The precise value is not physically important here, but using the recommended status keeps the event record conventional. See: https://pythia.org/latest-manual/HadronLevelStandalone.html
FromLeadingQuark	Custom status assigned to hadrons containing a leading quark.
FromLeadingDiquark	Custom status assigned to hadrons containing a leading diquark.

Definition at line 454 of file stringprocess.h.

                           : int {
    /// Custom status assigned to leading (valence) quarks.
    LeadingQuark = 202,
 
    /// Standard PYTHIA beam-remnant status used for leading diquarks.
    ///
    /// Leading diquarks are assigned status 63 so that PYTHIA treats them
    /// as beam remnants during hadronization. This allows PYTHIA's beam
    /// remnant machinery to be used for leading baryon production, including
    /// popcorn suppression and the optional hard-remnant-baryon treatment.
    ///
    /// In particular, when BeamRemnants:hardRemnantBaryon is enabled,
    /// PYTHIA replaces the standard Lund symmetric fragmentation function
    /// for the leading baryon by a dedicated remnant-baryon fragmentation
    /// function controlled by BeamRemnants:aRemnantBaryon and
    /// BeamRemnants:bRemnantBaryon. This produces harder leading baryons
    /// than ordinary string fragmentation.
    ///
    /// Therefore status 63 is intentionally used here rather than a custom
    /// status code.
    LeadingDiquark = 63,
 
    /// Standard PYTHIA status for non-leading partons that should be
    /// hadronized.
    ///
    /// The PYTHIA manual recommends status 23 for particles provided as
    /// input to standalone hadronization. The precise value is not
    /// physically important here, but using the recommended status keeps
    /// the event record conventional.
    ///
    /// See:
    /// https://pythia.org/latest-manual/HadronLevelStandalone.html
    NonLeadingParton = 23,
 
    /// Custom status assigned to hadrons containing a leading quark.
    FromLeadingQuark = 203,
 
    /// Custom status assigned to hadrons containing a leading diquark.
    FromLeadingDiquark = 204,
  };

Constructor & Destructor Documentation

StringProcess()

smash::StringProcess::StringProcess ( Configuration &  config )

explicit

Constructor, initializes PYTHIA.

Should only be called once.

All parameters are taken from config via Configuration::take.

Parameters

[in,out]

config

SMASH configuration object (keys are consumed)

Definition at line 29 of file stringprocess.cc.

    : pmin_gluon_lightcone_(
          config.take(InputKeys::collTerm_stringParam_gluonPMin)),
      pow_fgluon_beta_(config.take(InputKeys::collTerm_stringParam_gluonBeta)),
      pow_fquark_alpha_(
          config.take(InputKeys::collTerm_stringParam_quarkAlpha)),
      pow_fquark_beta_(config.take(InputKeys::collTerm_stringParam_quarkBeta)),
      sigma_qperp_(config.take(InputKeys::collTerm_stringParam_sigmaPerp)),
      stringz_a_leading_(
          config.take(InputKeys::collTerm_stringParam_stringZALeading)),
      stringz_b_leading_(
          config.take(InputKeys::collTerm_stringParam_stringZBLeading)),
      stringz_a_produce_(config.take(InputKeys::collTerm_stringParam_stringZA)),
      stringz_b_produce_(config.take(InputKeys::collTerm_stringParam_stringZB)),
      strange_supp_(
          config.take(InputKeys::collTerm_stringParam_strangeSuppression)),
      diquark_supp_(
          config.take(InputKeys::collTerm_stringParam_diquarkSuppression)),
      popcorn_rate_(config.take(InputKeys::collTerm_stringParam_popcornRate)),
      damp_popcorn_(config.take(InputKeys::collTerm_stringParam_dampPopcorn)),
      string_sigma_T_(
          config.take(InputKeys::collTerm_stringParam_stringSigmaT)),
      kappa_tension_string_(
          config.take(InputKeys::collTerm_stringParam_stringTension)),
      time_formation_const_(
          config.take(InputKeys::collTerm_stringParam_formationTime)),
      soft_t_form_(config.take(InputKeys::collTerm_stringParam_formTimeFactor)),
      mass_dependent_formation_times_(config.take(
          InputKeys::collTerm_stringParam_mDependentFormationTimes)),
      prob_proton_to_d_uu_(
          config.take(InputKeys::collTerm_stringParam_probabilityPToDUU)),
      separate_fragment_baryon_(
          config.take(InputKeys::collTerm_stringParam_separateFragmentBaryon)),
      use_monash_tune_(
          config.take(InputKeys::collTerm_stringParam_useMonashTune)),
      additional_xsec_supp_(
          config.take(InputKeys::collTerm_stringParam_unformedXsecSuppression)),
      pythia_settings_(
          config.take(InputKeys::collTerm_stringParam_pythiaSettings)) {
  // setup and initialize pythia for fragmentation
  pythia_hadron_ = std::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_,
                      popcorn_rate_, stringz_a_produce_, stringz_b_produce_,
                      string_sigma_T_);
  pythia_hadron_->init();
  /*
   * The const_cast<type>() function is used to obtain the reference of the
   * PrivateInfo object in the pythia_hadron_.
   * This cast is needed since Pythia 8.302 which included a major architecture
   * change. The Info object of a Pythia object is now private, only a const
   * reference can be obtained.
   * In order to reference the PrivateInfo object during initialization, we
   * cast the const reference to obtain the stored address.
   */
  pythia_sigmatot_.initInfoPtr(
      const_cast<Pythia8::Info&>(pythia_hadron_->info));
  pythia_sigmatot_.init();
 
  pythia_stringflav_.initInfoPtr(
      const_cast<Pythia8::Info&>(pythia_hadron_->info));
  pythia_stringflav_.init();
 
  event_intermediate_.init("intermediate partons",
                           &pythia_hadron_->particleData);
 
  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

compute_beam_valence_flags()

std::vector< bool > smash::StringProcess::compute_beam_valence_flags ( Pythia8::Pythia &  pythia )

private

Compute flags identifying beam valence partons (quarks or diquarks) that act as leading partons after the initial interaction.

The returned vector is aligned with the current Pythia event record (pythia.event): each entry is true if the corresponding particle is a valence parton originating from one of the incoming beam particles, i.e. a quark or diquark that carries the beam quantum numbers and should be considered "leading".

This typically includes valence quarks extracted during the initial scattering, as well as surviving diquark remnants of the beam.

Note

The size of the returned vector equals pythia.event.size().

Parameters

[in,out]

pythia

Pythia instance containing the event to inspect.

Returns

Per-particle flags for beam valence (leading) partons.

Definition at line 775 of file stringprocess.cc.

                           {
  const Pythia8::Event& event = pythia.event;
 
  auto find_final_copy = [&](int iPos) -> int {
    if (iPos <= 0 || iPos >= event.size())
      return -1;
    if (event[iPos].isFinal())
      return iPos;
 
    const int id = event[iPos].id();
    const auto ds = event[iPos].daughterListRecursive();
    for (int j : ds) {
      if (j > 0 && j < event.size() && event[j].isFinal() &&
          event[j].id() == id)
        return j;
    }
    return -1;
  };
 
  auto tag_from_beam = [&](const Pythia8::BeamParticle& beam,
                           std::vector<bool>& isValenceFinal) -> void {
    for (int i = 0; i < beam.size(); ++i) {
      if (!beam[i].isValence()) {
        continue;
      }
 
      const int j = find_final_copy(beam[i].iPos());
 
      if (j >= 0 && (event[j].isQuark() || event[j].isDiquark())) {
        isValenceFinal[j] = true;
      }
    }
  };
  auto tag_from_unresolved_diff_systems =
      [&](std::vector<bool>& isValenceFinal) -> void {
    for (int i = 1; i < event.size(); ++i) {
      if (event[i].statusAbs() != 15) {
        continue;
      }
 
      std::vector<int> endpoint_daughters;
      bool resolved_system = false;
 
      for (int d : event[i].daughterList()) {
        if (d <= 0 || d >= event.size()) {
          continue;
        }
 
        if (event[d].isQuark() || event[d].isDiquark()) {
          endpoint_daughters.push_back(d);
          continue;
        }
 
        if (!event[d].isGluon()) {
          resolved_system = true;
        }
      }
 
      if (resolved_system || endpoint_daughters.empty()) {
        continue;
      }
 
      bool already_tagged = false;
      for (int j : endpoint_daughters) {
        if (isValenceFinal[j]) {
          already_tagged = true;
          break;
        }
      }
 
      if (already_tagged) {
        continue;
      }
 
      for (int j : endpoint_daughters) {
        isValenceFinal[j] = true;
      }
    }
  };
 
  std::vector<bool> isValenceFinal(event.size(), false);
  tag_from_beam(pythia.beamA, isValenceFinal);
  tag_from_beam(pythia.beamB, isValenceFinal);
 
  tag_from_unresolved_diff_systems(isValenceFinal);
 
  return isValenceFinal;
}

Here is the caller graph for this function:

leading_hadron_status_from_endpoint()

int smash::StringProcess::leading_hadron_status_from_endpoint ( const Pythia8::Particle &  end )

inlineprivate

Determine the custom leading-hadron status code from a string endpoint.

For a string endpoint particle, return the appropriate SMASH leading status: diquark endpoints map to FROM_Leading_Diquark, otherwise FROM_Leading_Quark.

Parameters

[in]

end

String endpoint particle (quark or diquark).

Returns

Integer status code corresponding to the appropriate LeadingStatus.

Definition at line 222 of file stringprocess.h.

                                                                             {
    return static_cast<int>(end.isDiquark() ? LeadingStatus::FromLeadingDiquark
                                            : LeadingStatus::FromLeadingQuark);
  }

Here is the caller graph for this function:

is_leading_parton()

bool smash::StringProcess::is_leading_parton ( const Pythia8::Particle &  p )

inlineprivate

Check whether a particle is tagged as a leading parton.

Parameters

[in]

p

Pythia particle.

Returns

True if p has statusAbs() equal to LeadingStatus::Leading_Quark or LeadingStatus::Leading_Diquark.

Definition at line 234 of file stringprocess.h.

                                                          {
    return p.statusAbs() == static_cast<int>(LeadingStatus::LeadingQuark) ||
           p.statusAbs() == static_cast<int>(LeadingStatus::LeadingDiquark);
  }

Here is the caller graph for this function:

is_leading_from_quark()

bool smash::StringProcess::is_leading_from_quark ( const Pythia8::Particle &  p )

inlineprivate

Check whether a particle is tagged as originating from a leading quark.

Parameters

[in]

p

Pythia particle.

Returns

True if p has statusAbs() equal to LeadingStatus::FROM_Leading_Quark.

Definition at line 246 of file stringprocess.h.

                                                              {
    return p.statusAbs() == static_cast<int>(LeadingStatus::FromLeadingQuark);
  }

Here is the caller graph for this function:

is_leading_from_diquark()

bool smash::StringProcess::is_leading_from_diquark ( const Pythia8::Particle &  p )

inlineprivate

Check whether a particle is tagged as originating from a leading diquark.

Parameters

[in]

p

Pythia particle.

Returns

True if p has statusAbs() equal to LeadingStatus::FROM_Leading_Diquark.

Definition at line 257 of file stringprocess.h.

                                                                {
    return p.statusAbs() == static_cast<int>(LeadingStatus::FromLeadingDiquark);
  }

Here is the caller graph for this function:

is_leading()

bool smash::StringProcess::is_leading ( const Pythia8::Particle &  p )

inlineprivate

Check whether a particle is tagged as originating from a leading endpoint.

This is the union of is_leading_from_quark() and is_leading_from_diquark().

Parameters

[in]

p

Pythia particle.

Returns

True if p is tagged as leading-from-quark or leading-from-diquark.

Definition at line 272 of file stringprocess.h.

                                                   {
    return is_leading_from_quark(p) || is_leading_from_diquark(p);
  }

Here is the call graph for this function:

Here is the caller graph for this function:

append_string()

bool smash::StringProcess::append_string ( const Pythia8::Vec4 &  p_str, const std::array< int, 2 > &  ends, int  color_tag, bool  use_projectile_axis, bool  random_flip_of_endpoints = false  )

private

Append a single two-endpoint string as an independent PYTHIA event.

A minimal partonic event is constructed containing a string with the specified endpoint flavors and total four-momentum. The endpoint momenta are determined from two-body kinematics in the string rest frame and are aligned with either the projectile or target beam direction before being boosted back to the collision frame.

Strings whose invariant mass is below the estimated fragmentation threshold are rejected.

Parameters

[in]	p_str	Four-momentum of the string.
[in]	ends	PDG ids of the two string endpoints.
[in]	color_tag	Color tag used to connect the string endpoints.
[in]	use_projectile_axis	Whether to align the string with the projectile beam direction. If false, the target beam direction is used.
[in]	random_flip_of_endpoints	Whether to randomly reverse the order of the two string endpoints before assigning endpoint momenta.

Returns

True if a valid string event was created and appended, false otherwise.

Definition at line 433 of file stringprocess.cc.

                                                                 {
  const bool flip_endpoints =
      random_flip_of_endpoints && random::canonical() > 0.5;
 
  const int id1 = flip_endpoints ? ends[1] : ends[0];
  const int id2 = flip_endpoints ? ends[0] : ends[1];
  const auto& pd = pythia_hadron_->particleData;
 
  const double m_string = p_str.mCalc();
 
  if (m_string < estimate_string_threshold(id1, id2))
    return false;
 
  const double m1 = pd.m0(id1);
  const double m2 = pd.m0(id2);
 
  const double p_cm = pCM(m_string, m1, m2);
  if (p_cm < really_small)
    return false;
 
  const double E1 = std::sqrt(p_cm * p_cm + m1 * m1);
  const double E2 = std::sqrt(p_cm * p_cm + m2 * m2);
 
  const int ibeam = use_projectile_axis ? 0 : 1;
 
  const ThreeVector mom = pcom_[ibeam].threevec();
  Pythia8::Vec4 beam_axis_cm(mom.x1(), mom.x2(), mom.x3(), mom.abs());
 
  Pythia8::RotBstMatrix to_string_rest;
  to_string_rest.bstback(p_str);
 
  Pythia8::Vec4 beam_axis_rest = to_string_rest * beam_axis_cm;
 
  const double norm = std::sqrt(beam_axis_rest.px() * beam_axis_rest.px() +
                                beam_axis_rest.py() * beam_axis_rest.py() +
                                beam_axis_rest.pz() * beam_axis_rest.pz());
 
  if (norm < really_small)
    return false;
 
  const double nx = beam_axis_rest.px() / norm;
  const double ny = beam_axis_rest.py() / norm;
  const double nz = beam_axis_rest.pz() / norm;
 
  const Pythia8::Vec4 p1_rest(nx * p_cm, ny * p_cm, nz * p_cm, E1);
  const Pythia8::Vec4 p2_rest(-nx * p_cm, -ny * p_cm, -nz * p_cm, E2);
 
  Pythia8::RotBstMatrix boost;
  boost.bst(p_str);
 
  const Pythia8::Vec4 p1_cm = boost * p1_rest;
  const Pythia8::Vec4 p2_cm = boost * p2_rest;
 
  // Create event for this string.
  string_parton_events_.emplace_back();
  Pythia8::Event& evt = string_parton_events_.back();
  evt.init("Soft SMASH string", &pythia_hadron_->particleData);
 
  evt.append(90, -11, 0, 0, p1_cm + p2_cm, (p1_cm + p2_cm).mCalc());
 
  const int status_quark = static_cast<int>(LeadingStatus::LeadingQuark);
  const int status_diquark = static_cast<int>(LeadingStatus::LeadingDiquark);
 
  const bool id1_is_quark = pythia_hadron_->particleData.isQuark(id1);
  const bool id2_is_quark = pythia_hadron_->particleData.isQuark(id2);
 
  const int i1 = evt.append(id1, id1_is_quark ? status_quark : status_diquark,
                            0, 0, p1_cm, m1);
  const int i2 = evt.append(id2, id2_is_quark ? status_quark : status_diquark,
                            0, 0, p2_cm, m2);
 
  set_color_by_type(evt[i1], color_tag);
  set_color_by_type(evt[i2], color_tag);
  evt.rotbst(to_cm_.inverse());
  return true;
}

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 )

private

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

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 513 of file stringprocess.cc.

                                               {
  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_ * M_SQRT1_2);
  QTry = random::normal(0., sigma_qperp_ * M_SQRT1_2);
  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;
  }
 
  // determine three momentum of the final state hadron
  double sign_direction = is_AB_to_AX ? 1. : -1.;
 
  // longitudinal component (along +z by convention)
  const double pL = std::sqrt(pabscomHX_sqr - QTrn * QTrn);
 
  // momentum of the outgoing hadron H in the CM frame
  const ThreeVector cm_momentum(sign_direction * QTrx, sign_direction * QTry,
                                sign_direction * pL);
 
  const FourVector pstrHcom(std::sqrt(pabscomHX_sqr + massH * massH),
                            cm_momentum);
  const FourVector pstrXcom(std::sqrt(pabscomHX_sqr + massX * massX),
                            -cm_momentum);
 
  if (!append_string(make_pythia_4vec(pstrXcom), {idqX2, idqX1}, 101,
                     !is_AB_to_AX)) {
    return false;
  }
  PdgCode hadron_code = is_AB_to_AX ? PDGcodes_[0] : PDGcodes_[1];
  ParticleData new_particle(ParticleType::find(hadron_code));
 
  auto had_mom = make_pythia_4vec(pstrHcom);
  had_mom.rotbst(to_cm_.inverse());
 
  new_particle.set_4momentum(make_smash_4vec(had_mom));
  new_particle.set_cross_section_scaling_factor(1.);
  new_particle.set_formation_time(time_collision_);
  final_state_.push_back(new_particle);
 
  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 ( )

private

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 [8], Bleicher:1999xi [13].

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

Returns

whether the process is successfully implemented.

Definition at line 581 of file stringprocess.cc.

                               {
  std::array<std::array<int, 2>, 2> quarks;
 
  // decompose hadron into quark (and diquark) contents
  make_string_ends(PDGcodes_[0], quarks[0][1], quarks[0][0],
                   prob_proton_to_d_uu_);
  make_string_ends(PDGcodes_[1], quarks[1][1], quarks[1][0],
                   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_ * M_SQRT1_2);
  const double QTry = random::normal(0., sigma_qperp_ * M_SQRT1_2);
  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 pz1 = ((PPosA_ + QPos) - (PNegA_ + QNeg)) * M_SQRT1_2;
  const double E1 = ((PPosA_ + QPos) + (PNegA_ + QNeg)) * M_SQRT1_2;
  Pythia8::Vec4 p_str1(QTrx, QTry, pz1, E1);
 
  const double pz2 = ((PPosB_ - QPos) - (PNegB_ - QNeg)) * M_SQRT1_2;
  const double E2 = ((PPosB_ - QPos) + (PNegB_ - QNeg)) * M_SQRT1_2;
  Pythia8::Vec4 p_str2(-QTrx, -QTry, pz2, E2);
  if (random::uniform_int(0, 1) == 0) {
    std::swap(quarks[0][0], quarks[0][1]);
  }
  if (random::uniform_int(0, 1) == 0) {
    std::swap(quarks[1][0], quarks[1][1]);
  }
  return append_string(p_str1, quarks[0], 101, true) &&
         append_string(p_str2, quarks[1], 102, false);
}

Here is the call graph for this function:

Here is the caller graph for this function:

next_NDiffSoft()

bool smash::StringProcess::next_NDiffSoft ( )

private

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

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 [8], Bleicher:1999xi [13].

Returns

whether the process is successfully implemented.

Exceptions

std::runtime_error

if incoming particles are neither mesonic nor baryonic

Definition at line 715 of file stringprocess.cc.

                                   {
  // ids of the two ends of each string
  std::array<int, 2> endsA;  // string built with p_strA
  std::array<int, 2> endsB;  // string built with p_strB
 
  // 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();
  const int bar_b = PDGcodes_[1].baryon_number();
 
  if (bar_a == 1 ||                  // baryon-*
      (bar_a == 0 && bar_b == 1) ||  // meson-baryon
      (bar_a == 0 && bar_b == 0)) {  // meson-meson
    endsA = {idqA2, idqB1};
    endsB = {idqB2, idqA1};
  } else if ((bar_a == 0 && bar_b == -1) ||  // meson-antibaryon
             (bar_a == -1)) {                // antibaryon-*
    endsA = {idqB2, idqA1};
    endsB = {idqA2, idqB1};
  } else {
    std::stringstream ss;
    ss << "StringProcess::next_NDiffSoft: baryonA = " << bar_a
       << ", baryonB = " << bar_b;
    throw std::runtime_error(ss.str());
  }
 
  // CM frame with A along +z and B along -z; PPos/PNeg are lightcone components
  // in this frame.
 
  // 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 qx = random::normal(0., sigma_qperp_ * M_SQRT1_2);
  const double qy = random::normal(0., sigma_qperp_ * M_SQRT1_2);
  const double qT2 = qx * qx + qy * qy;
 
  // evaluate the lightcone momentum transfer
  const double QPos = -qT2 / (2. * xfracB * PNegB_);
  const double QNeg = qT2 / (2. * xfracA * PPosA_);
  const double dPPos = -xfracA * PPosA_ - QPos;
  const double dPNeg = xfracB * PNegB_ - QNeg;
 
  // string from hadron A side (with +qT)
  const double pzA = ((PPosA_ + dPPos) - (PNegA_ + dPNeg)) * M_SQRT1_2;
  const double EA = ((PPosA_ + dPPos) + (PNegA_ + dPNeg)) * M_SQRT1_2;
  Pythia8::Vec4 p_strA(qx, qy, pzA, EA);
 
  // string from hadron B side (with -qT)
  const double pzB = ((PPosB_ - dPPos) - (PNegB_ - dPNeg)) * M_SQRT1_2;
  const double EB = ((PPosB_ - dPPos) + (PNegB_ - dPNeg)) * M_SQRT1_2;
  Pythia8::Vec4 p_strB(-qx, -qy, pzB, EB);
  return append_string(p_strA, endsA, 101, true) &&
         append_string(p_strB, endsB, 102, false);
}

Here is the call graph for this function:

Here is the caller graph for this function:

next_Hard()

bool smash::StringProcess::next_Hard ( ProcessType  type )

private

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

Returns

whether the process is successfully implemented.

Definition at line 993 of file stringprocess.cc.

                                              {
  logg[LPythia].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]);
    logg[LPythia].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]);
    logg[LPythia].debug("    excess_quark[", i, "] = (", excess_quark[i][0],
                        ", ", excess_quark[i][1], ", ", excess_quark[i][2],
                        ", ", excess_quark[i][3], ", ", excess_quark[i][4],
                        ")");
    logg[LPythia].debug("    excess_antiq[", i, "] = (", excess_antiq[i][0],
                        ", ", excess_antiq[i][1], ", ", excess_antiq[i][2],
                        ", ", excess_antiq[i][3], ", ", excess_antiq[i][4],
                        ")");
  }
 
  std::pair<int, int> idAB{pdg_for_pythia[0], pdg_for_pythia[1]};
 
  // If an entry for the calculated particle IDs does not exist, create one
  // and initialize it accordingly
  if (hard_map_.count(idAB) == 0) {
    hard_map_[idAB] = std::make_unique<Pythia8::Pythia>(PYTHIA_XML_DIR, false);
    hard_map_[idAB]->readString("SoftQCD:nonDiffractive = on");
    hard_map_[idAB]->readString("SoftQCD:singleDiffractiveXB = on");
    hard_map_[idAB]->readString("SoftQCD:singleDiffractiveAX = on");
    hard_map_[idAB]->readString("SoftQCD:doubleDiffractive = on");
    hard_map_[idAB]->readString("HadronLevel:all = off");
    common_setup_pythia(hard_map_[idAB].get(), strange_supp_, diquark_supp_,
                        popcorn_rate_, stringz_a_produce_, stringz_b_produce_,
                        string_sigma_T_);
    hard_map_[idAB]->settings.flag("Beams:allowVariableEnergy", true);
    hard_map_[idAB]->settings.mode("Beams:idA", idAB.first);
    hard_map_[idAB]->settings.mode("Beams:idB", idAB.second);
    hard_map_[idAB]->settings.parm(
        "Beams:eCM", mpi_initialization_sqrts_.value_or(sqrtsAB_));
    logg[LPythia].debug("Pythia object initialized with ", pdg_for_pythia[0],
                        " + ", pdg_for_pythia[1], " at CM energy [GeV] ",
                        mpi_initialization_sqrts_.value_or(sqrtsAB_));
    if (!hard_map_[idAB]->init()) {
      throw std::runtime_error("Pythia failed to initialize.");
    }
  }
  const int seed_new = random::uniform_int(1, maximum_rndm_seed_in_pythia);
  hard_map_[idAB]->rndm.init(seed_new);
  logg[LPythia].debug("hard_map_[", idAB.first, "][", idAB.second,
                      "] : rndm is initialized with seed ", seed_new);
 
  // Change the energy using the Pythia 8.302+ feature
  hard_map_[idAB]->setKinematics(sqrtsAB_);
 
  bool final_state_success = false;
  /* Hard-process codes follow the Pythia convention:
   * https://www.pythia.org/latest-manual/QCDSoftProcesses.html
   *
   * Note: Pythia effectively classifies processes by the last digit
   * of the code (e.g. 101 -> 1, 102 -> 2, ...).
   *
   * The terminology "soft" vs "hard" is somewhat misleading:
   * Pythia's SoftQCD corresponds to a semi-perturbative (semi-hard)
   * model valid down to low pT, while HardQCD represents purely
   * perturbative processes and is not applicable over the full pT range.
   */
 
  switch (type) {
    case ProcessType::StringHardNonDiffractive:
      final_state_success = hard_map_[idAB]->next(1);
      break;
    case ProcessType::StringHardDoubleDiffractive:
      final_state_success = hard_map_[idAB]->next(5);
      break;
    case ProcessType::StringHardSingleDiffractiveAX:
      final_state_success = hard_map_[idAB]->next(4);
      break;
    case ProcessType::StringHardSingleDiffractiveXB:
      final_state_success = hard_map_[idAB]->next(3);
      break;
    default:
      logg[LPythia].error("Unknown string process required.");
      final_state_success = false;
      break;
  }
 
  logg[LPythia].debug("Pythia final state computed, success = ",
                      final_state_success);
  if (!final_state_success) {
    return false;
  }
 
  auto valence_tags = compute_beam_valence_flags(*hard_map_[idAB]);
 
  for (auto& p : hard_map_[idAB]->event) {
    if (!p.isQuark() && !p.isDiquark() && !p.isGluon()) {
      continue;
    }
 
    if (valence_tags[p.index()] && (p.isQuark() || p.isDiquark())) {
      p.statusCode(static_cast<int>(p.isDiquark()
                                        ? LeadingStatus::LeadingDiquark
                                        : LeadingStatus::LeadingQuark));
    } else {
      p.statusCode(static_cast<int>(LeadingStatus::NonLeadingParton));
    }
  }
  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 < hard_map_[idAB]->event.size(); i++) {
    if (hard_map_[idAB]->event[i].isFinal()) {
      const int pdgid = hard_map_[idAB]->event[i].id();
      Pythia8::Vec4 pquark = hard_map_[idAB]->event[i].p();
      const double mass = pquark.mCalc();
      const int status = hard_map_[idAB]->event[i].status();
      const int color = hard_map_[idAB]->event[i].col();
      const int anticolor = hard_map_[idAB]->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 < hard_map_[idAB]->event.sizeJunction(); i++) {
    const int kind = hard_map_[idAB]->event.kindJunction(i);
    std::array<int, 3> col;
    for (int j = 0; j < 3; j++) {
      col[j] = hard_map_[idAB]->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());
 
  const int elastic_side =
      type == ProcessType::StringHardSingleDiffractiveAX   ? 0
      : type == ProcessType::StringHardSingleDiffractiveXB ? 1
                                                           : -1;
  if (elastic_side >= 0 &&
      PDGcodes_[elastic_side].get_decimal() != pdg_for_pythia[elastic_side]) {
    const int mapped_id = pdg_for_pythia[elastic_side];
    const int actual_id = PDGcodes_[elastic_side].get_decimal();
 
    int elastic_index = -1;
 
    for (int i = 1; i < event_intermediate_.size(); ++i) {
      const auto& p = event_intermediate_[i];
 
      if (p.status() == 14 && p.id() == mapped_id) {
        elastic_index = i;
        break;
      }
    }
 
    if (elastic_index < 0) {
      logg[LPythia].warn("Could not find elastic mapped hadron ", mapped_id,
                         " to replace by ", actual_id);
      return false;
    }
 
    auto& elastic = event_intermediate_[elastic_index];
 
    std::vector<int> recoil_indices;
    Pythia8::Vec4 p_recoil_old;
 
    for (int i = 1; i < event_intermediate_.size(); ++i) {
      if (i == elastic_index)
        continue;
 
      const auto& p = event_intermediate_[i];
 
      recoil_indices.push_back(i);
      p_recoil_old += p.p();
    }
 
    if (recoil_indices.empty()) {
      event_intermediate_.list();
      logg[LPythia].warn(
          "Could not find recoil system for elastic mapped hadron ", mapped_id);
      return false;
    }
 
    Pythia8::Vec4 p_elastic_new = elastic.p();
    Pythia8::Vec4 p_recoil_new = p_recoil_old;
 
    const double m_elastic_new =
        ParticleType::find(PdgCode::from_decimal(actual_id)).mass();
 
    const double m_recoil_system = p_recoil_old.mCalc();
 
    if (!Pythia8::pShift(p_elastic_new, p_recoil_new, m_elastic_new,
                         m_recoil_system)) {
      logg[LPythia].warn(
          "Could not shift momenta when remapping elastic hadron from ",
          mapped_id, " to ", actual_id, ".");
      return false;
    }
 
    Pythia8::RotBstMatrix recoil_bst;
    recoil_bst.bstback(p_recoil_old);
    recoil_bst.bst(p_recoil_new);
 
    elastic.id(actual_id);
    elastic.m(m_elastic_new);
    elastic.p(p_elastic_new);
 
    for (const int idx : recoil_indices) {
      auto& recoiler = event_intermediate_[idx];
 
      Pythia8::Vec4 p_new = recoiler.p();
      p_new.rotbst(recoil_bst);
      recoiler.p(p_new);
    }
    excess_quark[elastic_side] = {0, 0, 0, 0, 0};
    excess_antiq[elastic_side] = {0, 0, 0, 0, 0};
  }
  /* 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) {
    logg[LPythia].debug("failed to find correct partonic constituents.");
    return false;
  }
  event_intermediate_.rotbst(to_cm_.inverse());
  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() &&
        !hard_map_[idAB]->particleData.isOctetHadron(pdgid)) {
      logg[LPythia].debug("PDG ID from Pythia: ", pdgid);
      Pythia8::Vec4 momentum_pythia = event_intermediate_[ipart].p();
      FourVector momentum = make_smash_4vec(momentum_pythia);
 
      logg[LPythia].debug("4-momentum from Pythia: ", momentum);
      bool found_ptype =
          append_intermediate_list(pdgid, momentum, new_non_hadron_particles);
      if (!found_ptype) {
        logg[LPythia].warn("PDG ID ", pdgid,
                           " does not exist in ParticleType - start over.");
        return false;
      }
      event_intermediate_.remove(ipart, ipart);
      npart -= 1;
    } else {
      ipart += 1;
    }
  }
 
  for (ParticleData& non_hadron : new_non_hadron_particles) {
    non_hadron.set_cross_section_scaling_factor(1.);
    non_hadron.set_formation_time(time_collision_);
    final_state_.push_back(non_hadron);
  }
  logg[LPythia].debug("Hard non-diff: partonic process gives ",
                      event_intermediate_.size(), " partons.");
  bool find_forward_string = true;
  while (event_intermediate_.size() > 1) {
    Pythia8::Event string_event;
    string_event.init("Hard string", &pythia_hadron_->particleData);
    if (event_intermediate_.sizeJunction() > 0) {
      // identify string from a junction if there is any.
      compose_string_junction(find_forward_string, event_intermediate_,
                              string_event);
    } else {
      /* identify string from a most forward or backward parton.
       * if there is no junction. */
      compose_string_parton(find_forward_string, event_intermediate_,
                            string_event);
    }
    if (!string_above_threshold(string_event)) {
      return false;
    }
    string_parton_events_.push_back(string_event);
    find_forward_string = !find_forward_string;
  }
 
  return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

next_BBbarAnn()

bool smash::StringProcess::next_BBbarAnn ( )

private

Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca [8], Bleicher:1999xi [13] 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 620 of file stringprocess.cc.

                                  {
  logg[LPythia].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++) {
      ParticleData new_particle(ParticleType::find(PDGcodes_[i]));
      Pythia8::Vec4 pcom_pyth = make_pythia_4vec(pcom_[i]);
      pcom_pyth.rotbst(to_cm_.inverse());
      new_particle.set_4momentum(make_smash_4vec(pcom_pyth));
      new_particle.set_cross_section_scaling_factor(1.);
      new_particle.set_formation_time(time_collision_);
      final_state_.push_back(new_particle);
    }
    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;
  }
 
  for (int i = 0; i < 2; i++) {
    if (!append_string(make_pythia_4vec(pcom_[i]),
                       {remaining_quarks[i], remaining_antiquarks[i]}, 101 + i,
                       i == 0, true)) {
      return false;
    }
  }
  return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

common_setup_pythia()

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

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]	popcorn_rate	parameter (StringFlav:popcornRate) to determine the production rate of popcorn mesons from the diquark end of a string.
[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

pythia8302/share/Pythia8/xmldoc/FlavourSelection.xml

pythia8302/share/Pythia8/xmldoc/Fragmentation.xml

Definition at line 102 of file stringprocess.cc.

                                                               {
  // choose parametrization for mass-dependent width
  pythia_in->readString("ParticleData:modeBreitWigner = 4");
 
  // Global Lund fragmentation
  pythia_in->readString("StringZ:aLund = " + std::to_string(stringz_a));
  pythia_in->readString("StringZ:bLund = " + std::to_string(stringz_b));
  pythia_in->readString("BeamRemnants:dampPopcorn = " +
                        std::to_string(damp_popcorn_));
  if (separate_fragment_baryon_) {
    pythia_in->readString("BeamRemnants:hardRemnantBaryon = on");
    pythia_in->readString("BeamRemnants:aRemnantBaryon = " +
                          std::to_string(stringz_a_leading_));
    pythia_in->readString("BeamRemnants:bRemnantBaryon = " +
                          std::to_string(stringz_b_leading_));
  }
  // 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));
  pythia_in->readString("StringFlav:popcornRate = " +
                        std::to_string(popcorn_rate));
 
  // 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");
 
  if (use_monash_tune_) {
    pythia_in->readString("Tune:ee = 7");
    pythia_in->readString("Tune:pp = 14");
  }
 
  for (const auto& setting : pythia_settings_) {
    if (!pythia_in->readString(setting)) {
      throw std::runtime_error(
          "Failed to apply Pythia setting \"" + setting +
          "\". Please check that it is a valid Pythia 8 setting.");
    }
  }
}

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 435 of file stringprocess.h.

                                 {
    const int seed_new = random::uniform_int(1, maximum_rndm_seed_in_pythia);
 
    pythia_hadron_->rndm.init(seed_new);
    logg[LPythia].debug("pythia_hadron_ : rndm is initialized with seed ",
                        seed_new);
  }

Here is the call graph for this function:

tag_leading_hadrons()

void smash::StringProcess::tag_leading_hadrons

(

Pythia8::Event &

event

)

Tag leading hadrons in a hadronized string.

After string fragmentation, this function identifies the hadrons that originate from the leading string endpoints and assigns them a custom status code (see LeadingStatus).

The procedure is:

• identify the two string endpoint partons (quarks or diquarks),
• boost the event into the string rest frame,
• sort all final-state hadrons by longitudinal momentum,
• associate each endpoint with the hadron closest to its end of the string,
• assign a leading-hadron status based on the flavor content of the endpoint parton.

Diquark endpoints are processed before quark endpoints to ensure that leading baryons are preferentially associated with beam-remnant diquarks. A hadron can only be tagged once.

For diquark endpoints, only baryons with matching baryon-number sign are considered. Quark endpoints may tag either mesons or baryons, but baryons must again have a compatible baryon-number sign.

The event is temporarily transformed into the string rest frame during the identification procedure and restored to its original frame before returning.

Parameters

[in,out]

event

Hadronized Pythia event. The event is modified in-place by assigning custom status codes to identified leading hadrons.

Definition at line 865 of file stringprocess.cc.

                                                           {
  // Find the two string endpoint partons.
  int idx1 = -1;
  int idx2 = -1;
  const Pythia8::ParticleData& pd = pythia_hadron_->particleData;
 
  for (const auto& part : event) {
    if (part.isQuark() || part.isDiquark()) {
      if (idx1 < 0) {
        idx1 = part.index();
      } else {
        idx2 = part.index();
        break;
      }
    }
  }
 
  if (idx1 < 0 || idx2 < 0) {
    return;
  }
 
  // Work in the string rest frame. This makes the leading hadrons the ones
  // closest to the two string endpoints in longitudinal momentum.
  Pythia8::RotBstMatrix toRest;
  toRest.toCMframe(event[idx1].p(), event[idx2].p());
  event.rotbst(toRest);
 
  // Collect all final-state hadrons and order them from backward to forward.
  std::vector<int> hadrons;
  hadrons.reserve(event.size());
 
  for (int i = 0; i < event.size(); ++i) {
    if (event[i].isFinal() && event[i].isHadron()) {
      hadrons.push_back(i);
    }
  }
 
  std::sort(hadrons.begin(), hadrons.end(),
            [&](int a, int b) { return event[a].pz() < event[b].pz(); });
 
  // Identify the forward and backward string endpoints in the rest frame.
  const int forward_endpoint =
      event[idx1].pz() > event[idx2].pz() ? idx1 : idx2;
  const int backward_endpoint = forward_endpoint == idx1 ? idx2 : idx1;
 
  auto matches_endpoint = [&](int endpoint, int hadron) {
    const auto& end = event[endpoint];
    const auto& h = event[hadron];
 
    // Do not allow the same hadron to be tagged by both endpoints.
    if (is_leading(h)) {
      return false;
    }
 
    // A diquark endpoint should tag a leading baryon with the same
    // baryon-number sign.
    if (end.isDiquark()) {
      return pd.isBaryon(h.id()) && h.id() * end.id() > 0;
    }
 
    // If a quark endpoint tags a baryon, require a compatible baryon sign.
    if (pd.isBaryon(h.id())) {
      return h.id() * end.id() > 0;
    }
 
    // Mesons are compatible with quark endpoints.
    return true;
  };
 
  auto find_edge_hadron = [&](int endpoint, bool forward) -> int {
    // Search from the edge corresponding to this endpoint: largest pz for
    // the forward endpoint, smallest pz for the backward endpoint.
    if (forward) {
      for (auto it = hadrons.rbegin(); it != hadrons.rend(); ++it) {
        if (matches_endpoint(endpoint, *it)) {
          return *it;
        }
      }
    } else {
      for (int h : hadrons) {
        if (matches_endpoint(endpoint, h)) {
          return h;
        }
      }
    }
 
    return -1;
  };
 
  auto tag_endpoint = [&](int endpoint, bool forward) {
    // Only original leading partons should define leading hadrons.
    if (!is_leading_parton(event[endpoint])) {
      return;
    }
 
    const int hadron = find_edge_hadron(endpoint, forward);
 
    if (hadron >= 0) {
      event[hadron].status(
          leading_hadron_status_from_endpoint(event[endpoint]));
    }
  };
 
  auto tag_if_diquark = [&](int endpoint, bool forward) {
    if (event[endpoint].isDiquark()) {
      tag_endpoint(endpoint, forward);
    }
  };
 
  auto tag_if_quark = [&](int endpoint, bool forward) {
    if (!event[endpoint].isDiquark()) {
      tag_endpoint(endpoint, forward);
    }
  };
 
  // Tag diquark endpoints first. This lets baryonic leading hadrons claim the
  // appropriate edge hadrons before quark endpoints are considered.
  tag_if_diquark(backward_endpoint, false);
  tag_if_diquark(forward_endpoint, true);
 
  // Tag the remaining quark endpoints.
  tag_if_quark(backward_endpoint, false);
  tag_if_quark(forward_endpoint, true);
 
  // Restore the original event frame.
  event.rotbst(toRest.inverse());
}

Here is the call graph for this function:

Here is the caller graph for this function:

make_smash_4vec()

static FourVector smash::StringProcess::make_smash_4vec ( const Pythia8::Vec4 &  p )

inlinestatic

Convert a PYTHIA four-vector into a SMASH four-vector.

Parameters

[in]

p

PYTHIA four-vector.

Returns

Corresponding SMASH four-vector.

Definition at line 535 of file stringprocess.h.

                                                          {
    return FourVector(p.e(), p.px(), p.py(), p.pz());
  }

Here is the caller graph for this function:

make_pythia_4vec()

static Pythia8::Vec4 smash::StringProcess::make_pythia_4vec ( const FourVector &  p )

inlinestatic

Convert a SMASH four-vector into a PYTHIA four-vector.

Parameters

[in]

p

SMASH four-vector.

Returns

Corresponding PYTHIA four-vector.

Definition at line 545 of file stringprocess.h.

                                                           {
    return Pythia8::Vec4(p.x1(), p.x2(), p.x3(), p.x0());
  }

Here is the caller graph for this function:

next()

bool smash::StringProcess::next

(

ProcessType

type

)

Generate the next string process for a given process type.

Parameters

[in]

type

Type of string process to generate.

Returns

Whether the process was successfully generated.

Definition at line 259 of file stringprocess.cc.

                                         {
  string_parton_events_.clear();
  final_state_.clear();
 
  bool parton_level_success = false;
 
  switch (type) {
    case ProcessType::StringSoftNonDiffractive:
      parton_level_success = next_NDiffSoft();
      break;
    case ProcessType::StringSoftSingleDiffractiveAX:
      parton_level_success = next_SDiff(true);
      break;
    case ProcessType::StringSoftSingleDiffractiveXB:
      parton_level_success = next_SDiff(false);
      break;
    case ProcessType::StringSoftDoubleDiffractive:
      parton_level_success = next_DDiff();
      break;
    case ProcessType::StringSoftAnnihilation:
      parton_level_success = next_BBbarAnn();
      break;
    case ProcessType::StringHardNonDiffractive:
    case ProcessType::StringHardDoubleDiffractive:
    case ProcessType::StringHardSingleDiffractiveAX:
    case ProcessType::StringHardSingleDiffractiveXB:
      parton_level_success = next_Hard(type);
      break;
    default:
      logg[LPythia].error("Unknown string process required.");
      return false;
  }
 
  if (!parton_level_success) {
    return false;
  }
 
  for (const Pythia8::Event& string_event : string_parton_events_) {
    auto hadrons = hadronize(string_event);
 
    if (!hadrons) {
      final_state_.clear();
      return false;
    }
 
    final_state_.insert(final_state_.end(), hadrons->begin(), hadrons->end());
  }
 
  return true;
}

Here is the call graph for this function:

Here is the caller graph for this function:

string_above_threshold()

bool smash::StringProcess::string_above_threshold

(

const Pythia8::Event &

event

)

Check whether all strings in a PYTHIA event are above fragmentation threshold.

Parameters

[in]

event

PYTHIA event record containing the string partons.

Returns

True if all strings are above threshold.

Definition at line 1350 of file stringprocess.cc.

                                                                    {
  const double string_mass = event[0].mCalc();
  const auto& particle_data = pythia_hadron_->particleData;
 
  std::vector<int> endpoints;
  endpoints.reserve(event.size());
 
  for (const Pythia8::Particle& p : event) {
    if (!p.isFinal() || p.isGluon() || !p.isParton()) {
      continue;
    }
 
    const int id = p.id();
 
    if (!particle_data.isQuark(id) && !particle_data.isDiquark(id)) {
      logg[LPythia].error("Colored non-quark/diquark in threshold check: id=",
                          id, ", col=", p.col(), ", acol=", p.acol());
      return false;
    }
 
    endpoints.push_back(id);
  }
 
  // Closed gluon string.
  if (endpoints.size() < 2) {
    return string_mass > pion_mass * 2.0;
  }
 
  // Ordinary open string.
  if (endpoints.size() == 2) {
    return string_mass > estimate_string_threshold(endpoints[0], endpoints[1]);
  }
 
  // Baryonic/antibaryonic endpoint: combine two same-sign quarks.
  if (endpoints.size() == 3) {
    for (std::size_t i = 0; i < endpoints.size(); ++i) {
      for (std::size_t j = i + 1; j < endpoints.size(); ++j) {
        const int id_i = endpoints[i];
        const int id_j = endpoints[j];
 
        if (!particle_data.isQuark(id_i) || !particle_data.isQuark(id_j)) {
          continue;
        }
 
        if (id_i * id_j < 0) {
          continue;
        }
 
        const int diquark_id = diquark_from_quarks(id_i, id_j);
 
        const std::size_t k = 3 - i - j;
        const int remaining_id = endpoints[k];
 
        return string_mass >
               estimate_string_threshold(diquark_id, remaining_id);
      }
    }
 
    logg[LPythia].error(
        "Cannot reduce 3 string endpoints to quark-diquark "
        "threshold: ",
        endpoints[0], ", ", endpoints[1], ", ", endpoints[2]);
    return false;
  }
 
  logg[LPythia].error(
      "Unexpected number of string endpoints in threshold check: ",
      endpoints.size());
  return false;
}

Here is the call graph for this function:

Here is the caller graph for this function:

estimate_string_threshold()

double smash::StringProcess::estimate_string_threshold	(	int	p_left,
		int	p_right
	)

Estimate the minimum invariant mass required for a string to fragment.

Parameters

[in]	p_left	Index of the left string endpoint in the PYTHIA event.
[in]	p_right	Index of the right string endpoint in the PYTHIA event.

Returns

Estimated string threshold mass [GeV].

Definition at line 1297 of file stringprocess.cc.

                                                                       {
  constexpr double fragmentation_margin = 2.0 * pion_mass;  // GeV
 
  auto& particle_data = pythia_hadron_->particleData;
 
  const bool left_is_quark = particle_data.isQuark(left_endpoint_id);
  const bool right_is_quark = particle_data.isQuark(right_endpoint_id);
  const bool left_is_diquark = particle_data.isDiquark(left_endpoint_id);
  const bool right_is_diquark = particle_data.isDiquark(right_endpoint_id);
 
  if (!(left_is_quark || left_is_diquark) ||
      !(right_is_quark || right_is_diquark)) {
    throw std::invalid_argument(
        "String threshold estimate requires quark or diquark endpoints.");
  }
 
  if (left_is_diquark && right_is_diquark) {
    const int lightest_left_baryon =
        pythia_stringflav_.combineToLightest(left_endpoint_id, 2);  // u
    const int lightest_right_baryon =
        pythia_stringflav_.combineToLightest(right_endpoint_id, -2);  // anti-u
 
    return particle_data.m0(lightest_left_baryon) +
           particle_data.m0(lightest_right_baryon) + fragmentation_margin;
  }
 
  const int lightest_hadron =
      pythia_stringflav_.combineToLightest(left_endpoint_id, right_endpoint_id);
 
  const PdgCode pdg = PdgCode::from_decimal(lightest_hadron);
  return ParticleType::find(pdg).mass() + fragmentation_margin;
}

Here is the call graph for this function:

Here is the caller graph for this function:

set_color_by_type()

void smash::StringProcess::set_color_by_type	(	Pythia8::Particle &	p,
		int	color
	)

Set the color or anticolor index of a particle according to its type.

Quarks receive a color index, antiquarks receive an anticolor index, and diquarks are treated according to their PYTHIA color-flow convention.

Parameters

[out]	p	Particle whose color information is modified.
[in]	color	Color tag to assign.

Definition at line 1331 of file stringprocess.cc.

                                                                   {
  const int id = p.id();
  const bool is_anti = (id < 0);
 
  const auto& pd = pythia_hadron_->particleData;
  const bool is_quark = pd.isQuark(id);
  const bool is_diquark = pd.isDiquark(id);
 
  if ((is_quark && !is_anti) || (is_diquark && is_anti)) {
    // quark or antidiquark: color flows out
    p.cols(color, 0);
  } else if ((is_quark && is_anti) || (is_diquark && !is_anti)) {
    // antiquark or diquark: color flows in
    p.cols(0, color);
  } else {
    throw std::runtime_error("Cannot set color: not a quark or diquark (id=" +
                             std::to_string(id) + ")");
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

hadronize()

std::optional< ParticleList > smash::StringProcess::hadronize

(

const Pythia8::Event &

string_evt

)

Hadronize a single partonic string configuration using Pythia8 and convert the produced hadrons into SMASH particles.

The input event is interpreted as a single color-singlet string system. Fragmentation is performed in the string rest frame. The produced hadrons are converted to SMASH ParticleData objects, assigned formation times and cross-section scaling factors, and finally boosted out of the string rest frame.

Leading hadrons originating from valence quark or diquark endpoints may be identified and assigned reduced cross sections according to the leading-hadron prescription.

Parameters

[in]

string_evt

Partonic string event to hadronize.

Returns

List of fragmented hadrons on success.

std::nullopt if fragmentation or particle conversion fails.

Definition at line 310 of file stringprocess.cc.

                                    {
  ParticleList intermediate_particles;
  // usually 0 is "system", so need at least indices 1 and 2
  if (string_evt.size() < 3) {
    logg[LPythia].error("String event too small to hadronize.");
    return std::nullopt;
  }
 
  bool has_leading_parton = false;
 
  const bool has_junction = string_evt.sizeJunction() > 0;
 
  for (const Pythia8::Particle& particle : string_evt) {
    if (is_leading_parton(particle)) {
      has_leading_parton = true;
      break;
    }
  }
 
  const bool has_leading_and_is_open = has_leading_parton && !has_junction;
 
  pythia_hadron_->event.reset();
  pythia_hadron_->event[0].p(string_evt[0].p());
  pythia_hadron_->event[0].m(string_evt[0].p().mCalc());
  for (int i = 1; i < string_evt.size(); ++i) {
    pythia_hadron_->event.append(string_evt[i]);
  }
 
  pythia_hadron_->event.clearJunctions();
  for (int i = 0; i < string_evt.sizeJunction(); ++i) {
    pythia_hadron_->event.appendJunction(
        string_evt.kindJunction(i), string_evt.colJunction(i, 0),
        string_evt.colJunction(i, 1), string_evt.colJunction(i, 2));
  }
 
  const Pythia8::Vec4 p_str = string_evt[0].p();
 
  Pythia8::RotBstMatrix to_string_rest;
  to_string_rest.bstback(p_str);
 
  pythia_hadron_->event.rotbst(to_string_rest);
 
  if (!pythia_hadron_->forceHadronLevel(false)) {
    pythia_hadron_->event.list();
    logg[LPythia].error("Pythia fragmentation failed for one string.");
    return std::nullopt;
  }
 
  if (has_leading_and_is_open) {
    tag_leading_hadrons(pythia_hadron_->event);
  }
 
  const FourVector pString_smash(
      p_str.e(), ThreeVector(p_str.px(), p_str.py(), p_str.pz()));
 
  const double m2 = pString_smash.sqr();
 
  if (m2 <= 0.0) {
    logg[LPythia].error("String has non-positive invariant mass.");
    return std::nullopt;
  }
 
  // evecLong: unit vector along string momentum
  // fallback if |p| ~ 0
  ThreeVector evecLong(p_str.px(), p_str.py(), p_str.pz());
 
  const double pabs = evecLong.abs();
 
  if (pabs > 1e-12) {
    evecLong = evecLong / pabs;
  } else {
    evecLong = ThreeVector(0.0, 0.0, 1.0);
  }
 
  // Build intermediate SMASH particles from final hadrons
  intermediate_particles.reserve(pythia_hadron_->event.size());
  for (int i = 0; i < pythia_hadron_->event.size(); ++i) {
    const auto& particle = pythia_hadron_->event[i];
 
    if (!(particle.isFinal() && particle.isHadron())) {
      continue;
    }
 
    const FourVector mom_smash = make_smash_4vec(particle.p());
 
    int particle_id = particle.id();
    convert_KaonLS(particle_id);
 
    if (!append_intermediate_list(particle_id, mom_smash,
                                  intermediate_particles)) {
      logg[LPythia].error("Unknown hadron in SMASH during hadronization: PDG=",
                          particle.id());
 
      return std::nullopt;
    }
 
    if (is_leading_from_diquark(particle)) {
      intermediate_particles.back().set_cross_section_scaling_factor(
          2.0 / 3.0 * additional_xsec_supp_);
 
    } else if (is_leading_from_quark(particle)) {
      intermediate_particles.back().set_cross_section_scaling_factor(
          (pythia_hadron_->particleData.isBaryon(particle.id()) ? 1.0 / 3.0
                                                                : 0.5) *
          additional_xsec_supp_);
    } else {
      intermediate_particles.back().set_cross_section_scaling_factor(0.0);
    }
  }
 
  if (intermediate_particles.empty()) {
    logg[LPythia].error("Hadronization produced no final hadrons.");
    return std::nullopt;
  }
 
  const bool should_assign_scaling = has_junction && has_leading_parton;
 
  form_intermediate_particles(intermediate_particles, pString_smash, evecLong,
                              additional_xsec_supp_, should_assign_scaling);
 
  return intermediate_particles;
}

Here is the call graph for this function:

Here is the caller graph for this function:

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 [56].

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 618 of file stringprocess.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 654 of file stringprocess.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 664 of file stringprocess.h.

{ pow_fgluon_beta_ = betapow; }

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 673 of file stringprocess.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 681 of file stringprocess.h.

{ sigma_qperp_ = sigma_qperp; }

set_tension_string()

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

inline

set the string tension, which is used in form_intermediate_particles.

Parameters

kappa_string

is a value that we want to use for string tension.

Definition at line 686 of file stringprocess.h.

                                               {
    kappa_tension_string_ = kappa_string;
  }

set_mpi_initialization_sqrts()

void smash::StringProcess::set_mpi_initialization_sqrts ( double  sqrts )

inline

Set the center-of-mass energy used to initialize MPI-capable Pythia objects.

If no value is set, the initialization energy is determined from the incoming hadrons.

Parameters

sqrts

Center-of-mass energy [GeV] used for MPI initialization.

Definition at line 699 of file stringprocess.h.

                                                  {
    mpi_initialization_sqrts_ = sqrts;
  }

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 226 of file stringprocess.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();
 
  sqrtsAB_ = (plab_[0] + plab_[1]).abs();
  ucomAB_ = (plab_[0] + plab_[1]) / sqrtsAB_;
  vcomAB_ = ucomAB_.velocity();
 
  const Pythia8::Vec4 pA_lab = make_pythia_4vec(plab_[0]);
  const Pythia8::Vec4 pB_lab = make_pythia_4vec(plab_[1]);
 
  to_cm_.reset();
  to_cm_.toCMframe(pA_lab, pB_lab);
 
  const Pythia8::Vec4 pA_cm = to_cm_ * pA_lab;
  const Pythia8::Vec4 pB_cm = to_cm_ * pB_lab;
 
  pcom_[0] = make_smash_4vec(pA_cm);
  pcom_[1] = make_smash_4vec(pB_cm);
  ThreeVector evec_polar(pA_cm.px(), pA_cm.py(), pA_cm.pz());
  evec_polar /= std::sqrt(evec_polar.sqr());
 
  make_orthonormal_basis(evec_polar, 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 2169 of file stringprocess.cc.

                                                                   {
  assert(std::fabs(evec_polar.sqr() - 1.) < really_small);
 
  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(std::cos(theta) * std::cos(phi));
    evec_basis[1].set_x2(std::cos(theta) * std::sin(phi));
    evec_basis[1].set_x3(-std::sin(theta));
 
    evec_basis[2].set_x1(-std::sin(phi));
    evec_basis[2].set_x2(std::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.);
    }
  }
 
  assert(std::fabs(evec_basis[1] * evec_basis[2]) < really_small);
  assert(std::fabs(evec_basis[2] * evec_basis[0]) < really_small);
  assert(std::fabs(evec_basis[0] * evec_basis[1]) < really_small);
}

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 2218 of file stringprocess.cc.

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

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 1420 of file stringprocess.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 particle 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 1454 of file stringprocess.cc.

                                                                     {
  // 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;
    }
  }
  logg[LPythia].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 1539 of file stringprocess.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 1568 of file stringprocess.cc.

                                                       {
  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) {
      logg[LPythia].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()) {
          logg[LPythia].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 * M_SQRT1_2);
        const double qy = random::normal(0., sigma_qt_frag * M_SQRT1_2);
        // 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);
 
        logg[LPythia].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 constituents 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 1701 of file stringprocess.cc.

                                                 {
  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) {
      logg[LPythia].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 1758 of file stringprocess.cc.

                                                 {
  Pythia8::Vec4 pSum = event_intermediate[0].p();
  const double energy_init = pSum.e();
  logg[LPythia].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) {
      logg[LPythia].debug("Not enough quark constituents of flavor ",
                          iflav + 1);
      return false;
    }
 
    if (nantiq_total[iflav] + excess_antiq[0][iflav] + excess_antiq[1][iflav] <
        0) {
      logg[LPythia].debug("Not enough antiquark constituents of flavor ",
                          -(iflav + 1));
      return false;
    }
  }
 
  for (int ih = 0; ih < 2; ih++) {
    logg[LPythia].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], ")");
    logg[LPythia].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]);
          logg[LPythia].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]);
          logg[LPythia].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 = pquark.mCalc();
 
        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;
  }
  logg[LPythia].debug("  valence quark contents of hadons are recovered.");
 
  logg[LPythia].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) <=
        std::abs(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();
    }
    logg[LPythia].debug("  parton momenta are rescaled by factor of ",
                        rescale_factor);
  }
 
  logg[LPythia].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 1917 of file stringprocess.cc.

                                                                         {
  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);
  logg[LPythia].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);
  logg[LPythia].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) {
    logg[LPythia].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) {
        logg[LPythia].debug("  col_to_find ", col_to_find, " from i ", i, "(",
                            pdgid, ") found");
      }
      if (found_acol) {
        logg[LPythia].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) {
        logg[LPythia].error("No parton with col = ", col_to_find);
      }
      if (acol_to_find != 0) {
        logg[LPythia].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);
      logg[LPythia].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 2005 of file stringprocess.cc.

                                                                           {
  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);
  logg[LPythia].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. */
      logg[LPythia].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]);
            }
          }
          logg[LPythia].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 2108 of file stringprocess.cc.

                                                                     {
  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()) {
          logg[LPythia].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()) {
          logg[LPythia].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();
          logg[LPythia].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);
        logg[LPythia].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 991 of file stringprocess.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:

form_intermediate_particles()

void smash::StringProcess::form_intermediate_particles	(	ParticleList &	intermediate_particles,
		const FourVector &	pString,
		const ThreeVector &	evecLong,
		double	additional_xsec_supp = 1.0,
		bool	find_and_scale_leading = true
	)

Set formation times and cross-section scaling factors for fragmented hadrons as described in Andersson:1983ia [5].

The input particles are expected to be in the string rest frame. This function modifies them in place: particle momenta and production vertices are boosted out of the string rest frame, formation times are assigned, and optionally leading-hadron cross-section scaling factors are recomputed.

Parameters

[in,out]	intermediate_particles	Fragmented hadrons to process.
[in]	pString	Four-momentum of the string.
[in]	evecLong	Unit 3-vector along which the string is stretched.
[in]	additional_xsec_supp	Additional multiplicative factor applied to cross-section scaling, e.g. coherence or medium effects.
[in]	find_and_scale_leading	If true, identify leading hadrons from valence quark or diquark endpoints and assign their cross-section scaling factors according to the leading-hadron prescription.

Precondition

intermediate_particles is not empty.

intermediate_particles contains only hadrons.

Definition at line 184 of file stringprocess.cc.

                                 {
  assert(intermediate_particles.size() > 0);
 
  int bstring = 0;
  for (const ParticleData& data : intermediate_particles) {
    bstring += data.pdgcode().baryon_number();
  }
 
  const ThreeVector vstring = pString.velocity();
  if (find_and_scale_leading) {
    assign_all_scaling_factors(bstring, intermediate_particles, evecLong,
                               additional_xsec_supp);
  }
 
  for (ParticleData& particle : intermediate_particles) {
    const FourVector mom_in_string_restframe = particle.momentum();
    const ThreeVector velocity_in_string_restframe =
        mom_in_string_restframe.velocity();
    const double gamma_string = 1.0 / particle.inverse_gamma();
 
    const FourVector p_com = mom_in_string_restframe.lorentz_boost(-vstring);
    particle.set_4momentum(p_com);
    const double tau_prod =
        mass_dependent_formation_times_
            ? M_SQRT2 * particle.effective_mass() / kappa_tension_string_
            : time_formation_const_;
 
    const double t_prod_string = tau_prod * gamma_string;
 
    FourVector fragment_position(t_prod_string,
                                 t_prod_string * velocity_in_string_restframe);
 
    fragment_position = fragment_position.lorentz_boost(-vstring);
 
    particle.set_slow_formation_times(
        time_collision_,
        soft_t_form_ * fragment_position.x0() + time_collision_);
  }
}

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 1039 of file stringprocess.h.

                                                                             {
    const std::string s = std::to_string(pdgid);
    PdgCode pythia_code(s);
    ParticleTypePtr new_type = ParticleType::try_find(pythia_code);
    if (new_type) {
      ParticleData new_particle(ParticleType::find(pythia_code));
      new_particle.set_4momentum(momentum);
      intermediate_particles.push_back(new_particle);
      return true;
    } else {
      // if the particle does not exist in SMASH the pythia event is rerun
      return false;
    }
  }

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 1059 of file stringprocess.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 2225 of file stringprocess.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 2241 of file stringprocess.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 2254 of file stringprocess.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:

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 2334 of file stringprocess.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 valence 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 2317 of file stringprocess.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 2302 of file stringprocess.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 2381 of file stringprocess.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:

getPPosA()

double smash::StringProcess::getPPosA ( )

inline

Returns

forward lightcone momentum incoming particle A in CM-frame [GeV]

See also

PPosA_

Definition at line 1167 of file stringprocess.h.

{ return PPosA_; }

getPNegA()

double smash::StringProcess::getPNegA ( )

inline

Returns

backward lightcone momentum incoming particle Af in CM-frame [GeV]

See also

PNegA_

Definition at line 1173 of file stringprocess.h.

{ return PNegA_; }

getPPosB()

double smash::StringProcess::getPPosB ( )

inline

Returns

forward lightcone momentum incoming particle B in CM-frame [GeV]

See also

PPosB_

Definition at line 1179 of file stringprocess.h.

{ return PPosB_; }

getPnegB()

double smash::StringProcess::getPnegB ( )

inline

Returns

backward lightcone momentum incoming particle B in CM-frame [GeV]

See also

PNegB_

Definition at line 1185 of file stringprocess.h.

{ return PNegB_; }

get_massA()

double smash::StringProcess::get_massA ( )

inline

Returns

mass of incoming particle A [GeV]

See also

massA_

Definition at line 1191 of file stringprocess.h.

{ return massA_; }

get_massB()

double smash::StringProcess::get_massB ( )

inline

Returns

mass of incoming particle B [GeV]

See also

massB_

Definition at line 1197 of file stringprocess.h.

{ return massB_; }

get_sqrts()

double smash::StringProcess::get_sqrts ( )

inline

Returns

sqrt of mandelstam s [GeV]

See also

sqrtsAB_

Definition at line 1203 of file stringprocess.h.

{ return sqrtsAB_; }

get_PDGs()

std::array<PdgCode, 2> smash::StringProcess::get_PDGs ( )

inline

Returns

array with PDG Codes of incoming particles

See also

PDGcodes_

Definition at line 1209 of file stringprocess.h.

{ return PDGcodes_; }

get_plab()

std::array<FourVector, 2> smash::StringProcess::get_plab ( )

inline

Returns

momenta of incoming particles in lab frame [GeV]

See also

plab_

Definition at line 1215 of file stringprocess.h.

{ return plab_; }

get_pcom()

std::array<FourVector, 2> smash::StringProcess::get_pcom ( )

inline

Returns

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

See also

pcom_

Definition at line 1221 of file stringprocess.h.

{ return pcom_; }

get_ucom()

FourVector smash::StringProcess::get_ucom ( )

inline

Returns

velocity four vector of the COM in the lab frame

See also

ucomAB_

Definition at line 1227 of file stringprocess.h.

{ return ucomAB_; }

get_vcom()

ThreeVector smash::StringProcess::get_vcom ( )

inline

Returns

velocity three vector of the COM in the lab frame

See also

vcomAB_

Definition at line 1233 of file stringprocess.h.

{ return vcomAB_; }

get_tcoll()

double smash::StringProcess::get_tcoll ( )

inline

Returns

collision time

See also

time_collision_

Definition at line 1239 of file stringprocess.h.

{ return time_collision_; }

get_final_state()

ParticleList smash::StringProcess::get_final_state ( )

inline

Returns

final state

See also

final_state_

Definition at line 1245 of file stringprocess.h.

{ return final_state_; }

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 51 of file stringprocess.h.

PPosB_

double smash::StringProcess::PPosB_

private

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

Definition at line 54 of file stringprocess.h.

PNegA_

double smash::StringProcess::PNegA_

private

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

Definition at line 57 of file stringprocess.h.

PNegB_

double smash::StringProcess::PNegB_

private

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

Definition at line 60 of file stringprocess.h.

massA_

double smash::StringProcess::massA_

private

mass of incoming particle A [GeV]

Definition at line 62 of file stringprocess.h.

massB_

double smash::StringProcess::massB_

private

mass of incoming particle B [GeV]

Definition at line 64 of file stringprocess.h.

sqrtsAB_

double smash::StringProcess::sqrtsAB_

private

sqrt of Mandelstam variable s of collision [GeV]

Definition at line 66 of file stringprocess.h.

PDGcodes_

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

private

PdgCodes of incoming particles.

Definition at line 68 of file stringprocess.h.

plab_

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

private

momenta of incoming particles in the lab frame [GeV]

Definition at line 70 of file stringprocess.h.

pcom_

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

private

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

Definition at line 72 of file stringprocess.h.

ucomAB_

FourVector smash::StringProcess::ucomAB_

private

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

Definition at line 74 of file stringprocess.h.

vcomAB_

ThreeVector smash::StringProcess::vcomAB_

private

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

Definition at line 76 of file stringprocess.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 81 of file stringprocess.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 stringprocess.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 stringprocess.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 stringprocess.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 stringprocess.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 stringprocess.h.

stringz_a_leading_

double smash::StringProcess::stringz_a_leading_

private

parameter (StringZ:aLund) for the fragmentation function of leading baryon in soft non-diffractive string processes

Definition at line 109 of file stringprocess.h.

stringz_b_leading_

double smash::StringProcess::stringz_b_leading_

private

parameter (StringZ:bLund) for the fragmentation function of leading baryon in soft non-diffractive string processes

Definition at line 114 of file stringprocess.h.

stringz_a_produce_

double smash::StringProcess::stringz_a_produce_

private

parameter (StringZ:aLund) for the fragmentation function of other (produced) hadrons in soft non-diffractive string processes

Definition at line 119 of file stringprocess.h.

stringz_b_produce_

double smash::StringProcess::stringz_b_produce_

private

parameter (StringZ:bLund) for the fragmentation function of other (produced) hadrons in soft non-diffractive string processes

Definition at line 124 of file stringprocess.h.

strange_supp_

double smash::StringProcess::strange_supp_

private

strange quark suppression factor

Definition at line 126 of file stringprocess.h.

diquark_supp_

double smash::StringProcess::diquark_supp_

private

diquark suppression factor

Definition at line 128 of file stringprocess.h.

popcorn_rate_

double smash::StringProcess::popcorn_rate_

private

popcorn rate

Definition at line 130 of file stringprocess.h.

damp_popcorn_

double smash::StringProcess::damp_popcorn_

private

damp popcorn meson from diquark remnant endpoint rate

Definition at line 132 of file stringprocess.h.

string_sigma_T_

double smash::StringProcess::string_sigma_T_

private

transverse momentum spread in string fragmentation

Definition at line 134 of file stringprocess.h.

kappa_tension_string_

double smash::StringProcess::kappa_tension_string_

private

string tension [GeV/fm]

Definition at line 136 of file stringprocess.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 138 of file stringprocess.h.

soft_t_form_

double smash::StringProcess::soft_t_form_

private

factor to be multiplied to formation times in soft strings

Definition at line 140 of file stringprocess.h.

time_collision_

double smash::StringProcess::time_collision_

private

time of collision in the computational frame [fm]

Definition at line 142 of file stringprocess.h.

mass_dependent_formation_times_

bool smash::StringProcess::mass_dependent_formation_times_

private

Whether the formation time should depend on the mass of the fragment according to Andersson:1983ia [5] eq.

2.45:

\( \tau = \sqrt{2}\frac{m}{\kappa} \)

The formation time and position is not calculated directly using the yoyo model because the spacetime rapidity where a string fragment forms is not equal to the fragment's momentum space rapidity. This cannot be easily combined with possible interactions before the formation time.

Definition at line 154 of file stringprocess.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 159 of file stringprocess.h.

separate_fragment_baryon_

bool smash::StringProcess::separate_fragment_baryon_

private

Whether to use a separate fragmentation function for leading baryons.

Definition at line 162 of file stringprocess.h.

use_monash_tune_

bool smash::StringProcess::use_monash_tune_

private

Whether to use the monash tune Skands:2014pea [61] for all string processes.

Definition at line 168 of file stringprocess.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 174 of file stringprocess.h.

mpi_initialization_sqrts_

std::optional<double> smash::StringProcess::mpi_initialization_sqrts_

private

Optional center-of-mass energy used to initialize MPI-capable Pythia objects.

If unset, the initialization energy is determined from the incoming hadrons.

Definition at line 181 of file stringprocess.h.

pythia_settings_

std::vector<std::string> smash::StringProcess::pythia_settings_

private

Additional Pythia 8 settings passed to each internal Pythia instance.

These settings are applied after the corresponding SMASH string parameters and therefore override them if both configure the same Pythia setting.

Definition at line 190 of file stringprocess.h.

final_state_

ParticleList smash::StringProcess::final_state_

private

final state array which must be accessed after the collision

Definition at line 280 of file stringprocess.h.

hard_map_

pythia_map smash::StringProcess::hard_map_

private

Map object to contain the different pythia objects.

Definition at line 291 of file stringprocess.h.

pythia_hadron_

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

private

PYTHIA object used in fragmentation.

Definition at line 294 of file stringprocess.h.

pythia_sigmatot_

Pythia8::SigmaTotal smash::StringProcess::pythia_sigmatot_

private

An object to compute cross-sections.

Definition at line 297 of file stringprocess.h.

pythia_stringflav_

Pythia8::StringFlav smash::StringProcess::pythia_stringflav_

private

An object for the flavor selection in string fragmentation in the case of separate fragmentation function for leading baryon.

Definition at line 303 of file stringprocess.h.

event_intermediate_

Pythia8::Event smash::StringProcess::event_intermediate_

private

event record for intermediate partonic state in the hard string routine

Definition at line 309 of file stringprocess.h.

string_parton_events_

std::vector<Pythia8::Event> smash::StringProcess::string_parton_events_

private

PYTHIA event records containing string partons to be hadronized.

Definition at line 340 of file stringprocess.h.

to_cm_

Pythia8::RotBstMatrix smash::StringProcess::to_cm_

Rotation/boost matrix to transform particles to the center-of-mass frame.

Definition at line 558 of file stringprocess.h.

---

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

• src/include/smash/stringprocess.h
• src/stringprocess.cc