smash::ScatterAction Class Reference

#include <scatteraction.h>

ScatterAction is a special action which takes two incoming particles and performs a scattering, producing one or more final-state particles.

Definition at line 30 of file scatteraction.h.

Inheritance diagram for smash::ScatterAction:

Classes
class	InvalidScatterAction
	Thrown when ScatterAction is called to perform with unknown ProcessType. More...

Public Member Functions
	ScatterAction (const ParticleData &in_part1, const ParticleData &in_part2, double time, bool isotropic=false, double string_formation_time=1.0, double box_length=-1.0, bool is_total_parametrized=false)
	Construct a ScatterAction object. More...
void	add_collision (CollisionBranchPtr p)
	Add a new collision channel. More...
void	add_collisions (CollisionBranchList pv)
	Add several new collision channels at once. More...
double	transverse_distance_sqr () const
	Calculate the transverse distance of the two incoming particles in their local rest frame. More...
double	cov_transverse_distance_sqr () const
	Calculate the transverse distance of the two incoming particles in their local rest frame written in a covariant form. More...
double	mandelstam_s () const
	Determine the Mandelstam s variable,. More...
double	relative_velocity () const
	Get the relative velocity of the two incoming particles. More...
void	generate_final_state () override
	Generate the final-state of the scattering process. More...
double	get_total_weight () const override
	Get the total cross section of scattering particles. More...
double	get_partial_weight () const override
	Get the partial cross section of the chosen channel. More...
void	sample_angles (std::pair< double, double > masses, double kinetic_energy_cm) override
	Sample final-state angles in a 2->2 collision (possibly anisotropic). More...
void	add_all_scatterings (const ScatterActionsFinderParameters &finder_parameters)
	Add all possible scattering subprocesses for this action object. More...
void	set_parametrized_total_cross_section (const ScatterActionsFinderParameters &finder_parameters)
	Given the incoming particles, assigns the correct parametrization of the total cross section. More...
const CollisionBranchList &	collision_channels ()
	Get list of possible collision channels. More...
void	set_string_interface (StringProcess *str_proc)
	Set the StringProcess object to be used. More...
virtual double	cross_section () const
	Get the total cross section of the scattering particles, either from a parametrization, or from the sum of partials. More...
Public Member Functions inherited from smash::Action
	Action (const ParticleList &in_part, double time)
	Construct an action object with incoming particles and relative time. More...
	Action (const ParticleData &in_part, const ParticleData &out_part, double time, ProcessType type)
	Construct an action object with the incoming particles, relative time, and the already known outgoing particles and type of the process. More...
	Action (const ParticleList &in_part, const ParticleList &out_part, double absolute_execution_time, ProcessType type)
	Construct an action object with the incoming particles, absolute time, and the already known outgoing particles and type of the process. More...
	Action (const Action &)=delete
	Copying is disabled. Use pointers or create a new Action. More...
virtual	~Action ()
	Virtual Destructor. More...
bool	operator< (const Action &rhs) const
	Determine whether one action takes place before another in time. More...
virtual ProcessType	get_type () const
	Get the process type. More...
template<typename Branch >
void	add_process (ProcessBranchPtr< Branch > &p, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add a new subprocess. More...
template<typename Branch >
void	add_processes (ProcessBranchList< Branch > pv, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add several new subprocesses at once. More...
virtual double	perform (Particles *particles, uint32_t id_process)
	Actually perform the action, e.g. More...
bool	is_valid (const Particles &particles) const
	Check whether the action still applies. More...
bool	is_pauli_blocked (const std::vector< Particles > &ensembles, const PauliBlocker &p_bl) const
	Check if the action is Pauli-blocked. More...
const ParticleList &	incoming_particles () const
	Get the list of particles that go into the action. More...
void	update_incoming (const Particles &particles)
	Update the incoming particles that are stored in this action to the state they have in the global particle list. More...
const ParticleList &	outgoing_particles () const
	Get the list of particles that resulted from the action. More...
double	time_of_execution () const
	Get the time at which the action is supposed to be performed. More...
virtual double	check_conservation (const uint32_t id_process) const
	Check various conservation laws. More...
double	sqrt_s () const
	Determine the total energy in the center-of-mass frame [GeV]. More...
FourVector	total_momentum_of_outgoing_particles () const
	Calculate the total kinetic momentum of the outgoing particles. More...
FourVector	get_interaction_point () const
	Get the interaction point. More...
std::pair< FourVector, FourVector >	get_potential_at_interaction_point () const
	Get the skyrme and asymmetry potential at the interaction point. More...
void	set_stochastic_pos_idx ()
	Setter function that stores a random incoming particle index latter used to determine the interaction point. More...

Protected Member Functions
double	cm_momentum () const
	Get the momentum of the center of mass of the incoming particles in the calculation frame. More...
double	cm_momentum_squared () const
	Get the squared momentum of the center of mass of the incoming particles in the calculation frame. More...
ThreeVector	beta_cm () const
	Get the velocity of the center of mass of the scattering/incoming particles in the calculation frame. More...
double	gamma_cm () const
	Get the gamma factor corresponding to a boost to the center of mass frame of the colliding particles. More...
void	elastic_scattering ()
	Perform an elastic two-body scattering, i.e. just exchange momentum. More...
void	inelastic_scattering ()
	Perform an inelastic two-body scattering, i.e. new particles are formed. More...
void	two_to_many_scattering ()
	Perform an inelastic two-to-many-body scattering (more than 2) More...
void	create_string_final_state ()
	Creates the final states for string-processes after they are performed. More...
void	string_excitation ()
	Todo(ryu): document better - it is not really UrQMD-based, isn't it? Perform the UrQMD-based string excitation and decay. More...
void	format_debug_output (std::ostream &out) const override
	Writes information about this scatter action to the out stream. More...
Protected Member Functions inherited from smash::Action
FourVector	total_momentum () const
	Sum of 4-momenta of incoming particles. More...
template<typename Branch >
const Branch *	choose_channel (const ProcessBranchList< Branch > &subprocesses, double total_weight)
	Decide for a particular final-state channel via Monte-Carlo and return it as a ProcessBranch. More...
virtual std::pair< double, double >	sample_masses (double kinetic_energy_cm) const
	Sample final-state masses in general X->2 processes (thus also fixing the absolute c.o.m. More...
void	sample_2body_phasespace ()
	Sample the full 2-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...
virtual void	sample_manybody_phasespace ()
	Sample the full n-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...
void	assign_formation_time_to_outgoing_particles ()
	Assign the formation time to the outgoing particles. More...

Protected Attributes
CollisionBranchList	collision_channels_
	List of possible collisions. More...
double	sum_of_partial_cross_sections_
	Current sum of partial hadronic cross sections. More...
double	partial_cross_section_
	Partial cross-section to the chosen outgoing channel. More...
bool	isotropic_ = false
	Do this collision isotropically? More...
double	string_formation_time_ = 1.0
	Time fragments take to be fully formed in hard string excitation. More...
Protected Attributes inherited from smash::Action
ParticleList	incoming_particles_
	List with data of incoming particles. More...
ParticleList	outgoing_particles_
	Initially this stores only the PDG codes of final-state particles. More...
const double	time_of_execution_
	Time at which the action is supposed to be performed (absolute time in the lab frame in fm). More...
ProcessType	process_type_
	type of process More...
double	box_length_ = -1.0
	Box length: needed to determine coordinates of collision correctly in case of collision through the wall. More...
int	stochastic_position_idx_ = -1
	This stores a randomly-chosen index to an incoming particle. More...

Private Member Functions
bool	is_elastic () const
	Check if the scattering is elastic. More...
void	resonance_formation ()
	Perform a 2->1 resonance-formation process. More...
void	rescale_outgoing_branches ()
	Loop over the possible branches and rescales their weight according to the desired total cross section. More...
ParticleTypePtr	try_find_pseudoresonance (const PseudoResonance method, const StringTransitionParameters &transition) const
	Try to find a pseudo-resonance that can be created from the incoming particles using a given method. More...

Private Attributes
StringProcess *	string_process_ = nullptr
	Pointer to interface class for strings. More...
bool	is_total_parametrized_ = false
	Whether the total cross section is parametrized. More...
std::optional< double >	parametrized_total_cross_section_ = std::nullopt
	If cross section is parametrized, store the value. More...
bool	were_processes_added_ = false
	Lock for calling add_all_scatterings only once. More...

Static Private Attributes
static std::set< std::set< ParticleTypePtr > >	warned_no_rescaling_available {}
	Warn about zero cross section only once per particle type pair. More...

Additional Inherited Members
Static Public Member Functions inherited from smash::Action
static double	lambda_tilde (double a, double b, double c)
	Little helper function that calculates the lambda function (sometimes written with a tilde to better distinguish it) that appears e.g. More...
static void	sample_manybody_phasespace_impl (double sqrts, const std::vector< double > &m, std::vector< FourVector > &sampled_momenta)
	Implementation of the full n-body phase-space sampling (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...

Constructor & Destructor Documentation

ScatterAction()

smash::ScatterAction::ScatterAction	(	const ParticleData &	in_part1,
		const ParticleData &	in_part2,
		double	time,
		bool	isotropic = false,
		double	string_formation_time = 1.0,
		double	box_length = -1.0,
		bool	is_total_parametrized = false
	)

Construct a ScatterAction object.

Parameters

[in]	in_part1	first scattering partner
[in]	in_part2	second scattering partner
[in]	time	Time at which the action is supposed to take place
[in]	isotropic	if true, do the collision isotropically
[in]	string_formation_time	Time string fragments take to form
[in]	box_length	Passing box length to determine coordinate of the collision, in case it happened through the wall in a box. If negative, then there is no wrapping.
[in]	is_total_parametrized	Whether the total cross section used for collision finding is parametrized

Definition at line 28 of file scatteraction.cc.

    : Action({in_part_a, in_part_b}, time),
      sum_of_partial_cross_sections_(0.),
      isotropic_(isotropic),
      string_formation_time_(string_formation_time),
      is_total_parametrized_(is_total_parametrized) {
  box_length_ = box_length;
  if (is_total_parametrized_) {
    parametrized_total_cross_section_ = NAN;
  }
}

Member Function Documentation

add_collision()

void smash::ScatterAction::add_collision

(

CollisionBranchPtr

p

)

Add a new collision channel.

Parameters

[in]

p

Channel to be added.

Definition at line 43 of file scatteraction.cc.

                                                      {
  add_process<CollisionBranch>(p, collision_channels_,
                               sum_of_partial_cross_sections_);
}

add_collisions()

void smash::ScatterAction::add_collisions

(

CollisionBranchList

pv

)

Add several new collision channels at once.

Parameters

[in]

pv

list of channels to be added.

Definition at line 48 of file scatteraction.cc.

                                                         {
  add_processes<CollisionBranch>(std::move(pv), collision_channels_,
                                 sum_of_partial_cross_sections_);
}

transverse_distance_sqr()

double smash::ScatterAction::transverse_distance_sqr

(

)

const

Calculate the transverse distance of the two incoming particles in their local rest frame.

According to UrQMD criterion, Bass:1998ca [6] eq. (3.27):

• position of particle a: \(\mathbf{x}_a\)
• position of particle b: \(\mathbf{x}_b\)
• momentum of particle a: \(\mathbf{p}_a\)
• momentum of particle b: \(\mathbf{p}_b\)

\[ d^2_\mathrm{coll} = (\mathbf{x}_a - \mathbf{x}_b)^2 - \frac{\bigl[(\mathbf{x}_a - \mathbf{x}_b) \cdot (\mathbf{p}_a - \mathbf{p}_b)\bigr]^2 } {(\mathbf{p}_a - \mathbf{p}_b)^2} \]

Returns

squared distance \(d^2_\mathrm{coll}\).

Definition at line 344 of file scatteraction.cc.

                                                    {
  // local copy of particles (since we need to boost them)
  ParticleData p_a = incoming_particles_[0];
  ParticleData p_b = incoming_particles_[1];
  /* Boost particles to center-of-momentum frame. */
  const ThreeVector velocity = beta_cm();
  p_a.boost(velocity);
  p_b.boost(velocity);
  const ThreeVector pos_diff =
      p_a.position().threevec() - p_b.position().threevec();
  const ThreeVector mom_diff =
      p_a.momentum().threevec() - p_b.momentum().threevec();
 
  logg[LScatterAction].debug("Particle ", incoming_particles_,
                             " position difference [fm]: ", pos_diff,
                             ", momentum difference [GeV]: ", mom_diff);
 
  const double dp2 = mom_diff.sqr();
  const double dr2 = pos_diff.sqr();
  /* Zero momentum leads to infite distance. */
  if (dp2 < really_small) {
    return dr2;
  }
  const double dpdr = pos_diff * mom_diff;
 
  /* UrQMD squared distance criterion, in the center of momentum frame:
   * position of particle a: x_a
   * position of particle b: x_b
   * momentum of particle a: p_a
   * momentum of particle b: p_b
   * d^2_{coll} = (x_a - x_b)^2 - ((x_a - x_b) . (p_a - p_b))^2 / (p_a - p_b)^2
   */
  const double result = dr2 - dpdr * dpdr / dp2;
  return result > 0.0 ? result : 0.0;
}

cov_transverse_distance_sqr()

double smash::ScatterAction::cov_transverse_distance_sqr

(

)

const

Calculate the transverse distance of the two incoming particles in their local rest frame written in a covariant form.

Equivalent to the UrQMD transverse distance. See Hirano:2012yy [27] (5.6)-(5.11).

Returns

squared distance \(d^2_\mathrm{coll}\).

Definition at line 380 of file scatteraction.cc.

                                                        {
  // local copy of particles (since we need to boost them)
  ParticleData p_a = incoming_particles_[0];
  ParticleData p_b = incoming_particles_[1];
 
  const FourVector delta_x = p_a.position() - p_b.position();
  const double mom_diff_sqr =
      (p_a.momentum().threevec() - p_b.momentum().threevec()).sqr();
  const double x_sqr = delta_x.sqr();
 
  if (mom_diff_sqr < really_small) {
    return -x_sqr;
  }
 
  const double p_a_sqr = p_a.momentum().sqr();
  const double p_b_sqr = p_b.momentum().sqr();
  const double p_a_dot_x = p_a.momentum().Dot(delta_x);
  const double p_b_dot_x = p_b.momentum().Dot(delta_x);
  const double p_a_dot_p_b = p_a.momentum().Dot(p_b.momentum());
 
  const double b_sqr =
      -x_sqr -
      (p_a_sqr * std::pow(p_b_dot_x, 2) + p_b_sqr * std::pow(p_a_dot_x, 2) -
       2 * p_a_dot_p_b * p_a_dot_x * p_b_dot_x) /
          (std::pow(p_a_dot_p_b, 2) - p_a_sqr * p_b_sqr);
  return b_sqr > 0.0 ? b_sqr : 0.0;
}

mandelstam_s()

double smash::ScatterAction::mandelstam_s

(

)

const

Determine the Mandelstam s variable,.

\[s = (p_a + p_b)^2\]

Equal to the square of CMS energy.

Returns

Mandelstam s

Definition at line 321 of file scatteraction.cc.

{ return total_momentum().sqr(); }

relative_velocity()

double smash::ScatterAction::relative_velocity

(

)

const

Get the relative velocity of the two incoming particles.

For a defintion see e.g. Seifert:2017oyb [53], eq. (5)

Returns

relative velocity.

Definition at line 335 of file scatteraction.cc.

                                              {
  const double m1 = incoming_particles()[0].effective_mass();
  const double m2 = incoming_particles()[1].effective_mass();
  const double m_s = mandelstam_s();
  const double lamb = lambda_tilde(m_s, m1 * m1, m2 * m2);
  return std::sqrt(lamb) / (2. * incoming_particles()[0].momentum().x0() *
                            incoming_particles()[1].momentum().x0());
}

generate_final_state()

void smash::ScatterAction::generate_final_state ( )

overridevirtual

Generate the final-state of the scattering process.

Performs either elastic or inelastic scattering.

Exceptions

InvalidScatterAction

Implements smash::Action.

Reimplemented in smash::ScatterActionPhoton.

Definition at line 53 of file scatteraction.cc.

                                         {
  logg[LScatterAction].debug("Incoming particles: ", incoming_particles_);
 
  const CollisionBranch *proc = choose_channel<CollisionBranch>(
      collision_channels_, is_total_parametrized_
                               ? *parametrized_total_cross_section_
                               : sum_of_partial_cross_sections_);
  process_type_ = proc->get_type();
  outgoing_particles_ = proc->particle_list();
  partial_cross_section_ = proc->weight();
 
  logg[LScatterAction].debug("Chosen channel: ", process_type_,
                             outgoing_particles_);
 
  /* The production point of the new particles.  */
  FourVector middle_point = get_interaction_point();
 
  switch (process_type_) {
    case ProcessType::Elastic:
      /* 2->2 elastic scattering */
      elastic_scattering();
      break;
    case ProcessType::TwoToOne:
      /* resonance formation */
      resonance_formation();
      break;
    case ProcessType::TwoToTwo:
      /* 2->2 inelastic scattering */
      /* Sample the particle momenta in CM system. */
      inelastic_scattering();
      break;
    case ProcessType::TwoToThree:
    case ProcessType::TwoToFour:
    case ProcessType::TwoToFive:
      /* 2->m scattering */
      two_to_many_scattering();
      break;
    case ProcessType::StringSoftSingleDiffractiveAX:
    case ProcessType::StringSoftSingleDiffractiveXB:
    case ProcessType::StringSoftDoubleDiffractive:
    case ProcessType::StringSoftAnnihilation:
    case ProcessType::StringSoftNonDiffractive:
    case ProcessType::StringHard:
      string_excitation();
      break;
    default:
      throw InvalidScatterAction(
          "ScatterAction::generate_final_state: Invalid process type " +
          std::to_string(static_cast<int>(process_type_)) + " was requested. " +
          "(PDGcode1=" + incoming_particles_[0].pdgcode().string() +
          ", PDGcode2=" + incoming_particles_[1].pdgcode().string() + ")");
  }
 
  for (ParticleData &new_particle : outgoing_particles_) {
    // Boost to the computational frame
    new_particle.boost_momentum(
        -total_momentum_of_outgoing_particles().velocity());
    /* Set positions of the outgoing particles */
    if (proc->get_type() != ProcessType::Elastic) {
      new_particle.set_4position(middle_point);
    }
  }
}

get_total_weight()

double smash::ScatterAction::get_total_weight ( ) const

overridevirtual

Get the total cross section of scattering particles.

Returns

total cross section.

Implements smash::Action.

Reimplemented in smash::ScatterActionPhoton.

Definition at line 302 of file scatteraction.cc.

                                             {
  return sum_of_partial_cross_sections_ *
         incoming_particles_[0].xsec_scaling_factor() *
         incoming_particles_[1].xsec_scaling_factor();
}

get_partial_weight()

double smash::ScatterAction::get_partial_weight ( ) const

overridevirtual

Get the partial cross section of the chosen channel.

Returns

partial cross section.

Implements smash::Action.

Definition at line 308 of file scatteraction.cc.

                                               {
  return partial_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
         incoming_particles_[1].xsec_scaling_factor();
}

sample_angles()

void smash::ScatterAction::sample_angles ( std::pair< double, double >  masses, double  kinetic_energy_cm  )

overridevirtual

Sample final-state angles in a 2->2 collision (possibly anisotropic).

NN → NN: Choose angular distribution according to Cugnon parametrization, see Cugnon:1996kh [19].

NN → NΔ: Sample scattering angles in center-of-mass frame from an anisotropic angular distribution, using the same distribution as for elastic pp scattering, as suggested in Cugnon:1996kh [19].

NN → NR: Fit to HADES data, see Agakishiev:2014wqa [2].

Reimplemented from smash::Action.

Definition at line 464 of file scatteraction.cc.

                                                            {
  if (is_string_soft_process(process_type_) ||
      (process_type_ == ProcessType::StringHard)) {
    // We potentially have more than two particles, so the following angular
    // distributions don't work. Instead we just keep the angular
    // distributions generated by string fragmentation.
    return;
  }
  assert(outgoing_particles_.size() == 2);
 
  // NN scattering is anisotropic currently
  const bool nn_scattering = incoming_particles_[0].type().is_nucleon() &&
                             incoming_particles_[1].type().is_nucleon();
  /* Elastic process is anisotropic and
   * the angular distribution is based on the NN elastic scattering. */
  const bool el_scattering = process_type_ == ProcessType::Elastic;
 
  const double mass_in_a = incoming_particles_[0].effective_mass();
  const double mass_in_b = incoming_particles_[1].effective_mass();
 
  ParticleData *p_a = &outgoing_particles_[0];
  ParticleData *p_b = &outgoing_particles_[1];
 
  const double mass_a = masses.first;
  const double mass_b = masses.second;
 
  const std::array<double, 2> t_range = get_t_range<double>(
      kinetic_energy_cm, mass_in_a, mass_in_b, mass_a, mass_b);
  Angles phitheta;
  if (el_scattering && !isotropic_) {
    double mandelstam_s_new = 0.;
    if (nn_scattering) {
      mandelstam_s_new = mandelstam_s();
    } else {
      /* In the case of elastic collisions other than NN collisions,
       * there is an ambiguity on how to get the lab-frame momentum (plab),
       * since the incoming particles can have different masses.
       * Right now, we first obtain the center-of-mass momentum
       * of the collision (pcom_now).
       * Then, the lab-frame momentum is evaluated from the mandelstam s,
       * which yields the original center-of-mass momentum
       * when nucleon mass is assumed. */
      const double pcm_now = pCM_from_s(mandelstam_s(), mass_in_a, mass_in_b);
      mandelstam_s_new =
          4. * std::sqrt(pcm_now * pcm_now + nucleon_mass * nucleon_mass);
    }
    double bb, a, plab = plab_from_s(mandelstam_s_new);
    if (nn_scattering &&
        p_a->pdgcode().antiparticle_sign() ==
            p_b->pdgcode().antiparticle_sign() &&
        std::abs(p_a->type().charge() + p_b->type().charge()) == 1) {
      // proton-neutron and antiproton-antineutron
      bb = std::max(Cugnon_bnp(plab), really_small);
      a = (plab < 0.8) ? 1. : 0.64 / (plab * plab);
    } else {
      /* all others including pp, nn and AQM elastic processes
       * This is applied for all particle pairs, which are allowed to
       * interact elastically. */
      bb = std::max(Cugnon_bpp(plab), really_small);
      a = 1.;
    }
    double t = random::expo(bb, t_range[0], t_range[1]);
    if (random::canonical() > 1. / (1. + a)) {
      t = t_range[0] + t_range[1] - t;
    }
    // determine scattering angles in center-of-mass frame
    phitheta = Angles(2. * M_PI * random::canonical(),
                      1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
  } else if (nn_scattering && p_a->pdgcode().is_Delta() &&
             p_b->pdgcode().is_nucleon() &&
             p_a->pdgcode().antiparticle_sign() ==
                 p_b->pdgcode().antiparticle_sign() &&
             !isotropic_) {
    const double plab = plab_from_s(mandelstam_s());
    const double bb = std::max(Cugnon_bpp(plab), really_small);
    double t = random::expo(bb, t_range[0], t_range[1]);
    if (random::canonical() > 0.5) {
      t = t_range[0] + t_range[1] - t;  // symmetrize
    }
    phitheta = Angles(2. * M_PI * random::canonical(),
                      1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
  } else if (nn_scattering && p_b->pdgcode().is_nucleon() && !isotropic_ &&
             (p_a->type().is_Nstar() || p_a->type().is_Deltastar())) {
    const std::array<double, 4> p{1.46434, 5.80311, -6.89358, 1.94302};
    const double a = p[0] + mass_a * (p[1] + mass_a * (p[2] + mass_a * p[3]));
    /*  If the resonance is so heavy that the index "a" exceeds 30,
     *  the power function turns out to be too sharp. Take t directly to be
     *  t_0 in such a case. */
    double t = t_range[0];
    if (a < 30) {
      t = random::power(-a, t_range[0], t_range[1]);
    }
    if (random::canonical() > 0.5) {
      t = t_range[0] + t_range[1] - t;  // symmetrize
    }
    phitheta = Angles(2. * M_PI * random::canonical(),
                      1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
  } else {
    /* isotropic angular distribution */
    phitheta.distribute_isotropically();
  }
 
  ThreeVector pscatt = phitheta.threevec();
  // 3-momentum of first incoming particle in center-of-mass frame
  ThreeVector pcm =
      incoming_particles_[0].momentum().lorentz_boost(beta_cm()).threevec();
  pscatt.rotate_z_axis_to(pcm);
 
  // final-state CM momentum
  const double p_f = pCM(kinetic_energy_cm, mass_a, mass_b);
  if (!(p_f > 0.0)) {
    logg[LScatterAction].warn("Particle: ", p_a->pdgcode(),
                              " radial momentum: ", p_f);
    logg[LScatterAction].warn("Etot: ", kinetic_energy_cm, " m_a: ", mass_a,
                              " m_b: ", mass_b);
  }
  p_a->set_4momentum(mass_a, pscatt * p_f);
  p_b->set_4momentum(mass_b, -pscatt * p_f);
 
  /* Debug message is printed before boost, so that p_a and p_b are
   * the momenta in the center of mass frame and thus opposite to
   * each other.*/
  logg[LScatterAction].debug("p_a: ", *p_a, "\np_b: ", *p_b);
}

add_all_scatterings()

void smash::ScatterAction::add_all_scatterings

(

const ScatterActionsFinderParameters &

finder_parameters

)

Add all possible scattering subprocesses for this action object.

This can only be called once per ScatterAction instance.

Parameters

[in]

finder_parameters

parameters for collision finding.

Definition at line 117 of file scatteraction.cc.

                                                             {
  if (were_processes_added_) {
    logg[LScatterAction].fatal() << "Trying to add processes again.";
    throw std::logic_error(
        "add_all_scatterings should be called only once per ScatterAction "
        "instance");
  } else {
    were_processes_added_ = true;
  }
  CrossSections xs(incoming_particles_, sqrt_s(),
                   get_potential_at_interaction_point());
  CollisionBranchList processes =
      xs.generate_collision_list(finder_parameters, string_process_);
 
  /* Add various subprocesses.*/
  add_collisions(std::move(processes));
 
  /* If the string processes are not triggered by a probability, then they
   * always happen as long as the parametrized total cross section is larger
   * than the sum of the cross sections of the non-string processes, and the
   * square root s exceeds the threshold by at least 0.9 GeV. The cross section
   * of the string processes are counted by taking the difference between the
   * parametrized total and the sum of the non-strings. */
  if (!finder_parameters.strings_with_probability &&
      xs.string_probability(finder_parameters) > 0) {
    const double xs_diff =
        xs.high_energy(finder_parameters) - sum_of_partial_cross_sections_;
    if (xs_diff > 0.) {
      add_collisions(
          xs.string_excitation(xs_diff, string_process_, finder_parameters));
    }
  }
 
  ParticleTypePtr pseudoresonance =
      try_find_pseudoresonance(finder_parameters.pseudoresonance_method,
                               finder_parameters.transition_high_energy);
  if (pseudoresonance && finder_parameters.two_to_one) {
    const double xs_total = is_total_parametrized_
                                ? *parametrized_total_cross_section_
                                : xs.high_energy(finder_parameters);
    const double xs_gap = xs_total - sum_of_partial_cross_sections_;
    // The pseudo-resonance is only created if there is a (positive) cross
    // section gap
    if (xs_gap > really_small) {
      auto pseudoresonance_branch = std::make_unique<CollisionBranch>(
          *pseudoresonance, xs_gap, ProcessType::TwoToOne);
      add_collision(std::move(pseudoresonance_branch));
      logg[LScatterAction].debug()
          << "Pseudoresonance between " << incoming_particles_[0].type().name()
          << " and " << incoming_particles_[1].type().name() << "is"
          << pseudoresonance->name() << " with cross section " << xs_gap
          << "mb.";
    }
  }
  were_processes_added_ = true;
  // Rescale the branches so that their sum matches the parametrization
  if (is_total_parametrized_) {
    rescale_outgoing_branches();
  }
}

set_parametrized_total_cross_section()

void smash::ScatterAction::set_parametrized_total_cross_section

(

const ScatterActionsFinderParameters &

finder_parameters

)

Given the incoming particles, assigns the correct parametrization of the total cross section.

Parameters

[in]

finder_parameters

Parameters for collision finding.

Definition at line 285 of file scatteraction.cc.

                                                             {
  CrossSections xs(incoming_particles_, sqrt_s(),
                   get_potential_at_interaction_point());
 
  if (is_total_parametrized_) {
    parametrized_total_cross_section_ =
        xs.parametrized_total(finder_parameters);
  } else {
    logg[LScatterAction].fatal()
        << "Trying to parametrize total cross section when it shouldn't be.";
    throw std::logic_error(
        "This function can only be called on ScatterAction objects with "
        "parametrized cross section.");
  }
}

collision_channels()

const CollisionBranchList& smash::ScatterAction::collision_channels ( )

inline

Get list of possible collision channels.

Returns

list of possible collision channels.

Definition at line 165 of file scatteraction.h.

                                                  {
    return collision_channels_;
  }

set_string_interface()

void smash::ScatterAction::set_string_interface ( StringProcess *  str_proc )

inline

Set the StringProcess object to be used.

The StringProcess object is used to handle string excitation and to generate final state particles.

Parameters

[in]

str_proc

String process object to be used.

Definition at line 185 of file scatteraction.h.

                                                     {
    string_process_ = str_proc;
  }

cross_section()

virtual double smash::ScatterAction::cross_section ( ) const

inlinevirtual

Get the total cross section of the scattering particles, either from a parametrization, or from the sum of partials.

Returns

total cross section.

Definition at line 195 of file scatteraction.h.

                                       {
    if (is_total_parametrized_) {
      return *parametrized_total_cross_section_;
    }
    return sum_of_partial_cross_sections_;
  }

cm_momentum()

double smash::ScatterAction::cm_momentum ( ) const

protected

Get the momentum of the center of mass of the incoming particles in the calculation frame.

Returns

center of mass momentum.

Definition at line 323 of file scatteraction.cc.

                                        {
  const double m1 = incoming_particles_[0].effective_mass();
  const double m2 = incoming_particles_[1].effective_mass();
  return pCM(sqrt_s(), m1, m2);
}

cm_momentum_squared()

double smash::ScatterAction::cm_momentum_squared ( ) const

protected

Get the squared momentum of the center of mass of the incoming particles in the calculation frame.

Returns

center of mass momentum squared.

Definition at line 329 of file scatteraction.cc.

                                                {
  const double m1 = incoming_particles_[0].effective_mass();
  const double m2 = incoming_particles_[1].effective_mass();
  return pCM_sqr(sqrt_s(), m1, m2);
}

beta_cm()

ThreeVector smash::ScatterAction::beta_cm ( ) const

protected

Get the velocity of the center of mass of the scattering/incoming particles in the calculation frame.

Note: Do not use this function to boost the outgoing particles. Use total_momentum_of_outgoing_particles(), which corrects for the effect of potentials on intial and final state.

Returns

boost velocity between center of mass and calculation frame.

Definition at line 313 of file scatteraction.cc.

                                         {
  return total_momentum().velocity();
}

gamma_cm()

double smash::ScatterAction::gamma_cm ( ) const

protected

Get the gamma factor corresponding to a boost to the center of mass frame of the colliding particles.

Returns

gamma factor.

Definition at line 317 of file scatteraction.cc.

                                     {
  return (1. / std::sqrt(1.0 - beta_cm().sqr()));
}

elastic_scattering()

void smash::ScatterAction::elastic_scattering ( )

protected

Perform an elastic two-body scattering, i.e. just exchange momentum.

Definition at line 597 of file scatteraction.cc.

                                       {
  // copy initial particles into final state
  outgoing_particles_[0] = incoming_particles_[0];
  outgoing_particles_[1] = incoming_particles_[1];
  // resample momenta
  sample_angles({outgoing_particles_[0].effective_mass(),
                 outgoing_particles_[1].effective_mass()},
                sqrt_s());
}

inelastic_scattering()

void smash::ScatterAction::inelastic_scattering ( )

protected

Perform an inelastic two-body scattering, i.e. new particles are formed.

Definition at line 607 of file scatteraction.cc.

                                         {
  // create new particles
  sample_2body_phasespace();
  assign_formation_time_to_outgoing_particles();
}

two_to_many_scattering()

void smash::ScatterAction::two_to_many_scattering ( )

protected

Perform an inelastic two-to-many-body scattering (more than 2)

Definition at line 613 of file scatteraction.cc.

                                           {
  sample_manybody_phasespace();
  assign_formation_time_to_outgoing_particles();
  logg[LScatterAction].debug("2->", outgoing_particles_.size(),
                             " scattering:", incoming_particles_, " -> ",
                             outgoing_particles_);
}

create_string_final_state()

void smash::ScatterAction::create_string_final_state ( )

protected

Creates the final states for string-processes after they are performed.

Definition at line 643 of file scatteraction.cc.

                                              {
  outgoing_particles_ = string_process_->get_final_state();
  assign_formation_time_to_outgoing_particles();
  /* Check momentum difference for debugging */
  FourVector out_mom;
  for (ParticleData data : outgoing_particles_) {
    out_mom += data.momentum();
  }
  logg[LPythia].debug("Incoming momenta string:", total_momentum());
  logg[LPythia].debug("Outgoing momenta string:", out_mom);
}

string_excitation()

void smash::ScatterAction::string_excitation ( )

protected

Todo(ryu): document better - it is not really UrQMD-based, isn't it? Perform the UrQMD-based string excitation and decay.

Definition at line 659 of file scatteraction.cc.

                                      {
  assert(incoming_particles_.size() == 2);
  // Disable floating point exception trap for Pythia
  {
    DisableFloatTraps guard;
    /* initialize the string_process_ object for this particular collision */
    string_process_->init(incoming_particles_, time_of_execution_);
    /* implement collision */
    bool success = false;
    int ntry = 0;
    const int ntry_max = 10000;
    while (!success && ntry < ntry_max) {
      ntry++;
      switch (process_type_) {
        case ProcessType::StringSoftSingleDiffractiveAX:
          /* single diffractive to A+X */
          success = string_process_->next_SDiff(true);
          break;
        case ProcessType::StringSoftSingleDiffractiveXB:
          /* single diffractive to X+B */
          success = string_process_->next_SDiff(false);
          break;
        case ProcessType::StringSoftDoubleDiffractive:
          /* double diffractive */
          success = string_process_->next_DDiff();
          break;
        case ProcessType::StringSoftNonDiffractive:
          /* soft non-diffractive */
          success = string_process_->next_NDiffSoft();
          break;
        case ProcessType::StringSoftAnnihilation:
          /* soft BBbar 2 mesonic annihilation */
          success = string_process_->next_BBbarAnn();
          break;
        case ProcessType::StringHard:
          success = string_process_->next_NDiffHard();
          break;
        default:
          logg[LPythia].error("Unknown string process required.");
          success = false;
      }
    }
    if (ntry == ntry_max) {
      /* If pythia fails to form a string, it is usually because the energy
       * is not large enough. In this case, annihilation is then enforced. If
       * this process still does not not produce any results, it defaults to
       * an elastic collision. */
      bool success_newtry = false;
 
      /* Check if the initial state is a baryon-antibaryon state.*/
      PdgCode part1 = incoming_particles_[0].pdgcode(),
              part2 = incoming_particles_[1].pdgcode();
      bool is_BBbar_Pair = (part1.baryon_number() != 0) &&
                           (part1.baryon_number() == -part2.baryon_number());
 
      /* Decide on the new process .*/
      if (is_BBbar_Pair) {
        process_type_ = ProcessType::StringSoftAnnihilation;
      } else {
        process_type_ = ProcessType::StringSoftDoubleDiffractive;
      }
      /* Perform the new process*/
      int ntry_new = 0;
      while (!success_newtry && ntry_new < ntry_max) {
        ntry_new++;
        if (is_BBbar_Pair) {
          success_newtry = string_process_->next_BBbarAnn();
        } else {
          success_newtry = string_process_->next_DDiff();
        }
      }
 
      if (success_newtry) {
        create_string_final_state();
      }
 
      if (!success_newtry) {
        /* If annihilation fails:
         * Particles are normally added after process selection for
         * strings, outgoing_particles is still uninitialized, and memory
         * needs to be allocated. We also shift the process_type_ to elastic
         * so that sample_angles does a proper treatment. */
        outgoing_particles_.reserve(2);
        outgoing_particles_.push_back(ParticleData{incoming_particles_[0]});
        outgoing_particles_.push_back(ParticleData{incoming_particles_[1]});
        process_type_ = ProcessType::FailedString;
        elastic_scattering();
      }
    } else {
      create_string_final_state();
    }
  }
}

is_elastic()

bool smash::ScatterAction::is_elastic ( ) const

private

Check if the scattering is elastic.

Returns

whether the scattering is elastic.

resonance_formation()

void smash::ScatterAction::resonance_formation ( )

private

Perform a 2->1 resonance-formation process.

Exceptions

InvalidResonanceFormation

Definition at line 621 of file scatteraction.cc.

                                        {
  if (outgoing_particles_.size() != 1) {
    std::string s =
        "resonance_formation: "
        "Incorrect number of particles in final state: ";
    s += std::to_string(outgoing_particles_.size()) + " (";
    s += incoming_particles_[0].pdgcode().string() + " + ";
    s += incoming_particles_[1].pdgcode().string() + ")";
    throw InvalidResonanceFormation(s);
  }
  // Set the momentum of the formed resonance in its rest frame.
  outgoing_particles_[0].set_4momentum(
      total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
  assign_formation_time_to_outgoing_particles();
  /* this momentum is evaluated in the computational frame. */
  logg[LScatterAction].debug("Momentum of the new particle: ",
                             outgoing_particles_[0].momentum());
}

rescale_outgoing_branches()

void smash::ScatterAction::rescale_outgoing_branches ( )

private

Loop over the possible branches and rescales their weight according to the desired total cross section.

In case the current sum of partials is close to 0, a warning is issued as this would not happen in an usual run, and an elastic process is added to match the total.

Definition at line 179 of file scatteraction.cc.

                                              {
  if (!were_processes_added_) {
    logg[LScatterAction].fatal()
        << "Trying to rescale branches before adding processes.";
    throw std::logic_error(
        "This function can only be called after having added processes.");
  }
  if (sum_of_partial_cross_sections_ < really_small) {
    const ParticleTypePtr type_a = &incoming_particles_[0].type();
    const ParticleTypePtr type_b = &incoming_particles_[1].type();
    // This is a std::set instead of std::pair because the order of particles
    // does not matter here
    const std::set<ParticleTypePtr> pair{type_a, type_b};
    if (!warned_no_rescaling_available.count(pair)) {
      logg[LScatterAction].warn()
          << "Total cross section between " << type_a->name() << " and "
          << type_b->name() << "is roughly zero at sqrt(s) = " << sqrt_s()
          << " GeV, and no rescaling to match the parametrized value will be "
             "done.\nAn elastic process will be added, instead, to match the "
             "total cross section.\nFor this pair of particles, this warning "
             "will be subsequently suppressed.";
      warned_no_rescaling_available.insert(pair);
    }
    auto elastic_branch = std::make_unique<CollisionBranch>(
        *type_a, *type_b, *parametrized_total_cross_section_,
        ProcessType::Elastic);
    add_collision(std::move(elastic_branch));
  } else {
    const double reweight =
        *parametrized_total_cross_section_ / sum_of_partial_cross_sections_;
    logg[LScatterAction].debug("Reweighting ", sum_of_partial_cross_sections_,
                               " to ", *parametrized_total_cross_section_);
    for (auto &proc : collision_channels_) {
      proc->set_weight(proc->weight() * reweight);
    }
  }
}

try_find_pseudoresonance()

ParticleTypePtr smash::ScatterAction::try_find_pseudoresonance ( const PseudoResonance  method, const StringTransitionParameters &  transition  ) const

private

Try to find a pseudo-resonance that can be created from the incoming particles using a given method.

Parameters

[in]	method	used to select the pseudo-resonance among possible candidates. See user guide description for more information.
[in]	transition	parameters for the string transition region, which are also used to determine when a pseudo-resonance can be created.

Returns

the appropriate pseudo-resonance, if there is any, or an invalid pointer otherwise.

Definition at line 217 of file scatteraction.cc.

                                                        {
  const ParticleTypePtr type_a = &incoming_particles_[0].type();
  const ParticleTypePtr type_b = &incoming_particles_[1].type();
  const double desired_mass = sqrt_s();
 
  double string_offset = 0.5 * transition.sqrts_add_lower;
  const bool nucleon_and_pion = (type_a->is_nucleon() && type_b->is_pion()) ||
                                (type_a->is_pion() && type_b->is_nucleon());
  const bool two_nucleons = type_a->is_nucleon() && type_b->is_nucleon();
  const bool two_pions = type_a->is_pion() && type_b->is_pion();
  const bool nucleon_and_kaon = (type_a->is_nucleon() && type_b->is_kaon()) ||
                                (type_a->is_kaon() && type_b->is_nucleon());
  if (nucleon_and_pion) {
    string_offset = transition.sqrts_range_Npi.first - pion_mass - nucleon_mass;
  } else if (two_nucleons) {
    string_offset = transition.sqrts_range_NN.first - 2 * nucleon_mass;
  } else if (two_pions) {
    string_offset = transition.pipi_offset;
  } else if (nucleon_and_kaon) {
    string_offset = transition.KN_offset;
  }
  /*
   * Artificial cutoff to create a pseudo-resonance only close to the string
   * transition, where data or first-principle models are unhelpful
   */
  if (desired_mass < incoming_particles_[0].effective_mass() +
                         incoming_particles_[1].effective_mass() +
                         string_offset) {
    return {};
  }
 
  if (method == PseudoResonance::None) {
    return {};
  } else if (method == PseudoResonance::LargestFromUnstable ||
             method == PseudoResonance::ClosestFromUnstable) {
    if (type_a->is_stable() && type_b->is_stable()) {
      return {};
    }
  }
 
  // If this list is empty, there are no possible pseudo-resonances.
  ParticleTypePtrList list = list_possible_resonances(type_a, type_b);
  if (std::empty(list)) {
    return {};
  }
 
  if (method == PseudoResonance::Largest ||
      method == PseudoResonance::LargestFromUnstable) {
    auto largest = *std::max_element(list.begin(), list.end(),
                                     [](ParticleTypePtr a, ParticleTypePtr b) {
                                       return a->mass() < b->mass();
                                     });
    return largest;
  } else if (method == PseudoResonance::Closest ||
             method == PseudoResonance::ClosestFromUnstable) {
    auto comparison = [&desired_mass](ParticleTypePtr a, ParticleTypePtr b) {
      return std::abs(a->mass() - desired_mass) <
             std::abs(b->mass() - desired_mass);
    };
    auto closest = *std::min_element(list.begin(), list.end(), comparison);
    return closest;
  } else {
    throw std::logic_error("Unknown method for selecting pseudoresonance.");
  }
}

Member Data Documentation

collision_channels_

CollisionBranchList smash::ScatterAction::collision_channels_

protected

List of possible collisions.

Definition at line 261 of file scatteraction.h.

sum_of_partial_cross_sections_

double smash::ScatterAction::sum_of_partial_cross_sections_

protected

Current sum of partial hadronic cross sections.

Definition at line 264 of file scatteraction.h.

partial_cross_section_

double smash::ScatterAction::partial_cross_section_

protected

Partial cross-section to the chosen outgoing channel.

Definition at line 267 of file scatteraction.h.

isotropic_

bool smash::ScatterAction::isotropic_ = false

protected

Do this collision isotropically?

Definition at line 270 of file scatteraction.h.

string_formation_time_

double smash::ScatterAction::string_formation_time_ = 1.0

protected

Time fragments take to be fully formed in hard string excitation.

Definition at line 273 of file scatteraction.h.

string_process_

StringProcess* smash::ScatterAction::string_process_ = nullptr

private

Pointer to interface class for strings.

Definition at line 313 of file scatteraction.h.

is_total_parametrized_

bool smash::ScatterAction::is_total_parametrized_ = false

private

Whether the total cross section is parametrized.

Definition at line 316 of file scatteraction.h.

parametrized_total_cross_section_

std::optional<double> smash::ScatterAction::parametrized_total_cross_section_ = std::nullopt

private

If cross section is parametrized, store the value.

Definition at line 319 of file scatteraction.h.

were_processes_added_

bool smash::ScatterAction::were_processes_added_ = false

private

Lock for calling add_all_scatterings only once.

Definition at line 322 of file scatteraction.h.

warned_no_rescaling_available

std::set<std::set<ParticleTypePtr> > smash::ScatterAction::warned_no_rescaling_available {}

inlinestaticprivate

Warn about zero cross section only once per particle type pair.

Definition at line 326 of file scatteraction.h.

---

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

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