Version: SMASH-3.3
Classes | Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
smash::ScatterAction Class Reference
Action Classes

#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:
smash::Action smash::BremsstrahlungAction smash::ScatterActionPhoton

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, const SpinInteractionType spin_interaction_type=SpinInteractionType::Off)
 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...
 
void assign_unpolarized_spin_vector_to_outgoing_particles ()
 Assign an unpolarized spin vector to all outgoing particles. 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 spin_interaction ()
 Perform spin interaction in binary interactions. More...
 
void string_spin_interaction ()
 Perform spin interaction in string excitations. 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...
 
virtual 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...
 
SpinInteractionType spin_interaction_type_ = SpinInteractionType::Off
 What kind of spin interaction to use. 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,
const SpinInteractionType  spin_interaction_type = SpinInteractionType::Off 
)

Construct a ScatterAction object.

Parameters
[in]in_part1first scattering partner
[in]in_part2second scattering partner
[in]timeTime at which the action is supposed to take place
[in]isotropicif true, do the collision isotropically
[in]string_formation_timeTime string fragments take to form
[in]box_lengthPassing 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_parametrizedWhether the total cross section used for collision finding is parametrized
[in]spin_interaction_typeWhich type of spin interaction to use

Definition at line 28 of file scatteraction.cc.

35  : Action({in_part_a, in_part_b}, time),
36  sum_of_partial_cross_sections_(0.),
37  isotropic_(isotropic),
38  string_formation_time_(string_formation_time),
39  is_total_parametrized_(is_total_parametrized),
40  spin_interaction_type_(spin_interaction_type) {
41  box_length_ = box_length;
42  if (is_total_parametrized_) {
43  parametrized_total_cross_section_ = smash_NaN<double>;
44  }
45 }
smash::Action::Action
Action(const ParticleList &in_part, double time)
Construct an action object with incoming particles and relative time.
Definition: action.h:44
smash::Action::box_length_
double box_length_
Box length: needed to determine coordinates of collision correctly in case of collision through the w...
Definition: action.h:388
smash::ScatterAction::spin_interaction_type_
SpinInteractionType spin_interaction_type_
What kind of spin interaction to use.
Definition: scatteraction.h:338
smash::ScatterAction::isotropic_
bool isotropic_
Do this collision isotropically?
Definition: scatteraction.h:286
smash::ScatterAction::string_formation_time_
double string_formation_time_
Time fragments take to be fully formed in hard string excitation.
Definition: scatteraction.h:289
smash::ScatterAction::parametrized_total_cross_section_
std::optional< double > parametrized_total_cross_section_
If cross section is parametrized, store the value.
Definition: scatteraction.h:335
smash::ScatterAction::is_total_parametrized_
bool is_total_parametrized_
Whether the total cross section is parametrized.
Definition: scatteraction.h:332
smash::ScatterAction::sum_of_partial_cross_sections_
double sum_of_partial_cross_sections_
Current sum of partial hadronic cross sections.
Definition: scatteraction.h:280

Member Function Documentation

◆ add_collision()

void smash::ScatterAction::add_collision ( CollisionBranchPtr  p)

Add a new collision channel.

Parameters
[in]pChannel to be added.

Definition at line 47 of file scatteraction.cc.

47  {
48  add_process<CollisionBranch>(p, collision_channels_,
49  sum_of_partial_cross_sections_);
50 }
smash::ScatterAction::collision_channels_
CollisionBranchList collision_channels_
List of possible collisions.
Definition: scatteraction.h:277
smash::pdg::p
constexpr int p
Proton.
Definition: pdgcode_constants.h:28

◆ add_collisions()

void smash::ScatterAction::add_collisions ( CollisionBranchList  pv)

Add several new collision channels at once.

Parameters
[in]pvlist of channels to be added.

Definition at line 52 of file scatteraction.cc.

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

◆ 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 366 of file scatteraction.cc.

366  {
367  // local copy of particles (since we need to boost them)
368  ParticleData p_a = incoming_particles_[0];
369  ParticleData p_b = incoming_particles_[1];
370  /* Boost particles to center-of-momentum frame. */
371  const ThreeVector velocity = beta_cm();
372  p_a.boost(velocity);
373  p_b.boost(velocity);
374  const ThreeVector pos_diff =
375  p_a.position().threevec() - p_b.position().threevec();
376  const ThreeVector mom_diff =
377  p_a.momentum().threevec() - p_b.momentum().threevec();
378 
379  logg[LScatterAction].debug("Particle ", incoming_particles_,
380  " position difference [fm]: ", pos_diff,
381  ", momentum difference [GeV]: ", mom_diff);
382 
383  const double dp2 = mom_diff.sqr();
384  const double dr2 = pos_diff.sqr();
385  /* Zero momentum leads to infite distance. */
386  if (dp2 < really_small) {
387  return dr2;
388  }
389  const double dpdr = pos_diff * mom_diff;
390 
391  /* UrQMD squared distance criterion, in the center of momentum frame:
392  * position of particle a: x_a
393  * position of particle b: x_b
394  * momentum of particle a: p_a
395  * momentum of particle b: p_b
396  * d^2_{coll} = (x_a - x_b)^2 - ((x_a - x_b) . (p_a - p_b))^2 / (p_a - p_b)^2
397  */
398  const double result = dr2 - dpdr * dpdr / dp2;
399  return result > 0.0 ? result : 0.0;
400 }
smash::Action::incoming_particles_
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:364
smash::ScatterAction::beta_cm
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering/incoming particles in the calculation frame.
Definition: scatteraction.cc:335
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:40
smash::really_small
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:41
smash::LScatterAction
static constexpr int LScatterAction
Definition: bremsstrahlungaction.cc:17

◆ 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 [28] (5.6)-(5.11).

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

Definition at line 402 of file scatteraction.cc.

402  {
403  // local copy of particles (since we need to boost them)
404  ParticleData p_a = incoming_particles_[0];
405  ParticleData p_b = incoming_particles_[1];
406 
407  const FourVector delta_x = p_a.position() - p_b.position();
408  const double mom_diff_sqr =
409  (p_a.momentum().threevec() - p_b.momentum().threevec()).sqr();
410  const double x_sqr = delta_x.sqr();
411 
412  if (mom_diff_sqr < really_small) {
413  return -x_sqr;
414  }
415 
416  const double p_a_sqr = p_a.momentum().sqr();
417  const double p_b_sqr = p_b.momentum().sqr();
418  const double p_a_dot_x = p_a.momentum().Dot(delta_x);
419  const double p_b_dot_x = p_b.momentum().Dot(delta_x);
420  const double p_a_dot_p_b = p_a.momentum().Dot(p_b.momentum());
421 
422  const double b_sqr =
423  -x_sqr -
424  (p_a_sqr * std::pow(p_b_dot_x, 2) + p_b_sqr * std::pow(p_a_dot_x, 2) -
425  2 * p_a_dot_p_b * p_a_dot_x * p_b_dot_x) /
426  (std::pow(p_a_dot_p_b, 2) - p_a_sqr * p_b_sqr);
427  return b_sqr > 0.0 ? b_sqr : 0.0;
428 }

◆ 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 343 of file scatteraction.cc.

343 { return total_momentum().sqr(); }
smash::Action::total_momentum
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:398
smash::FourVector::sqr
double sqr() const
calculate the square of the vector (which is a scalar)
Definition: fourvector.h:460

◆ 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 [55], eq. (5)

Returns
relative velocity.

Definition at line 357 of file scatteraction.cc.

357  {
358  const double m1 = incoming_particles()[0].effective_mass();
359  const double m2 = incoming_particles()[1].effective_mass();
360  const double m_s = mandelstam_s();
361  const double lamb = lambda_tilde(m_s, m1 * m1, m2 * m2);
362  return std::sqrt(lamb) / (2. * incoming_particles()[0].momentum().x0() *
363  incoming_particles()[1].momentum().x0());
364 }
smash::Action::incoming_particles
const ParticleList & incoming_particles() const
Get the list of particles that go into the action.
Definition: action.cc:58
smash::Action::lambda_tilde
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 ...
Definition: action.h:315
smash::ScatterAction::mandelstam_s
double mandelstam_s() const
Determine the Mandelstam s variable,.
Definition: scatteraction.cc:343

◆ 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 57 of file scatteraction.cc.

57  {
58  logg[LScatterAction].debug("Incoming particles: ", incoming_particles_);
59 
60  const CollisionBranch *proc = choose_channel<CollisionBranch>(
61  collision_channels_, is_total_parametrized_
62  ? *parametrized_total_cross_section_
63  : sum_of_partial_cross_sections_);
64  process_type_ = proc->get_type();
65  outgoing_particles_ = proc->particle_list();
66  partial_cross_section_ = proc->weight();
67 
68  logg[LScatterAction].debug("Chosen channel: ", process_type_,
69  outgoing_particles_);
70 
71  /* The production point of the new particles. */
72  FourVector middle_point = get_interaction_point();
73 
74  switch (process_type_) {
75  case ProcessType::Elastic:
76  /* 2->2 elastic scattering */
77  elastic_scattering();
78  spin_interaction();
79  break;
80  case ProcessType::TwoToOne:
81  /* resonance formation */
82  resonance_formation();
83  break;
84  case ProcessType::TwoToTwo:
85  /* 2->2 inelastic scattering */
86  /* Sample the particle momenta in CM system. */
87  inelastic_scattering();
88  spin_interaction();
89  break;
90  case ProcessType::TwoToThree:
91  case ProcessType::TwoToFour:
92  case ProcessType::TwoToFive:
93  /* 2->m scattering */
94  two_to_many_scattering();
95  break;
96  case ProcessType::StringSoftSingleDiffractiveAX:
97  case ProcessType::StringSoftSingleDiffractiveXB:
98  case ProcessType::StringSoftDoubleDiffractive:
99  case ProcessType::StringSoftAnnihilation:
100  case ProcessType::StringSoftNonDiffractive:
101  case ProcessType::StringHard:
102  string_excitation();
103  string_spin_interaction();
104  break;
105  default:
106  throw InvalidScatterAction(
107  "ScatterAction::generate_final_state: Invalid process type " +
108  std::to_string(static_cast<int>(process_type_)) + " was requested. " +
109  "(PDGcode1=" + incoming_particles_[0].pdgcode().string() +
110  ", PDGcode2=" + incoming_particles_[1].pdgcode().string() + ")");
111  }
112 
113  const bool core_in_incoming =
114  std::any_of(incoming_particles_.begin(), incoming_particles_.end(),
115  [](const ParticleData &p) { return p.is_core(); });
116  for (ParticleData &new_particle : outgoing_particles_) {
117  // Boost to the computational frame
118  new_particle.boost_momentum(
119  -total_momentum_of_outgoing_particles().velocity());
120  /* Set positions of the outgoing particles */
121  if (proc->get_type() != ProcessType::Elastic) {
122  new_particle.set_4position(middle_point);
123  if (core_in_incoming) {
124  new_particle.fluidize();
125  }
126  if (new_particle.type().pdgcode().is_heavy_flavor()) {
127  // Particle weight is the product of incoming weights
128  const double perturbative_weight = std::accumulate(
129  incoming_particles_.begin(), incoming_particles_.end(), 1.0,
130  [](double w, const ParticleData &p) {
131  return w * p.perturbative_weight();
132  });
133  new_particle.set_perturbative_weight(perturbative_weight);
134  }
135  }
136  }
137 }
smash::Action::total_momentum_of_outgoing_particles
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:160
smash::Action::outgoing_particles_
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:372
smash::Action::get_interaction_point
FourVector get_interaction_point() const
Get the interaction point.
Definition: action.cc:68
smash::Action::process_type_
ProcessType process_type_
type of process
Definition: action.h:381
smash::ScatterAction::string_spin_interaction
void string_spin_interaction()
Perform spin interaction in string excitations.
Definition: scatteraction.cc:864
smash::ScatterAction::resonance_formation
void resonance_formation()
Perform a 2->1 resonance-formation process.
Definition: scatteraction.cc:649
smash::ScatterAction::partial_cross_section_
double partial_cross_section_
Partial cross-section to the chosen outgoing channel.
Definition: scatteraction.h:283
smash::ScatterAction::two_to_many_scattering
void two_to_many_scattering()
Perform an inelastic two-to-many-body scattering (more than 2)
Definition: scatteraction.cc:638
smash::ScatterAction::string_excitation
void string_excitation()
Todo(ryu): document better - it is not really UrQMD-based, isn't it? Perform the UrQMD-based string e...
Definition: scatteraction.cc:690
smash::ScatterAction::elastic_scattering
void elastic_scattering()
Perform an elastic two-body scattering, i.e. just exchange momentum.
Definition: scatteraction.cc:619
smash::ScatterAction::inelastic_scattering
void inelastic_scattering()
Perform an inelastic two-body scattering, i.e. new particles are formed.
Definition: scatteraction.cc:629
smash::ScatterAction::spin_interaction
void spin_interaction()
Perform spin interaction in binary interactions.
Definition: scatteraction.cc:795
smash::ProcessType::TwoToOne
@ TwoToOne
See here for a short description.
smash::ProcessType::StringSoftDoubleDiffractive
@ StringSoftDoubleDiffractive
See here for a short description.
smash::ProcessType::TwoToFive
@ TwoToFive
See here for a short description.
smash::ProcessType::StringSoftSingleDiffractiveXB
@ StringSoftSingleDiffractiveXB
See here for a short description.
smash::ProcessType::TwoToTwo
@ TwoToTwo
See here for a short description.
smash::ProcessType::Elastic
@ Elastic
See here for a short description.
smash::ProcessType::TwoToFour
@ TwoToFour
See here for a short description.
smash::ProcessType::StringSoftAnnihilation
@ StringSoftAnnihilation
See here for a short description.
smash::ProcessType::StringSoftNonDiffractive
@ StringSoftNonDiffractive
See here for a short description.
smash::ProcessType::StringSoftSingleDiffractiveAX
@ StringSoftSingleDiffractiveAX
See here for a short description.
smash::ProcessType::StringHard
@ StringHard
See here for a short description.
smash::ProcessType::TwoToThree
@ TwoToThree
See here for a short description.
smash::to_string
std::string to_string(ThermodynamicQuantity quantity)
Convert a ThermodynamicQuantity enum value to its corresponding string.
Definition: stringify.cc:26

◆ 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 324 of file scatteraction.cc.

324  {
325  return sum_of_partial_cross_sections_ *
326  incoming_particles_[0].xsec_scaling_factor() *
327  incoming_particles_[1].xsec_scaling_factor();
328 }

◆ 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 330 of file scatteraction.cc.

330  {
331  return partial_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
332  incoming_particles_[1].xsec_scaling_factor();
333 }

◆ 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 486 of file scatteraction.cc.

487  {
488  if (is_string_soft_process(process_type_) ||
489  (process_type_ == ProcessType::StringHard)) {
490  // We potentially have more than two particles, so the following angular
491  // distributions don't work. Instead we just keep the angular
492  // distributions generated by string fragmentation.
493  return;
494  }
495  assert(outgoing_particles_.size() == 2);
496 
497  // NN scattering is anisotropic currently
498  const bool nn_scattering = incoming_particles_[0].type().is_nucleon() &&
499  incoming_particles_[1].type().is_nucleon();
500  /* Elastic process is anisotropic and
501  * the angular distribution is based on the NN elastic scattering. */
502  const bool el_scattering = process_type_ == ProcessType::Elastic;
503 
504  const double mass_in_a = incoming_particles_[0].effective_mass();
505  const double mass_in_b = incoming_particles_[1].effective_mass();
506 
507  ParticleData *p_a = &outgoing_particles_[0];
508  ParticleData *p_b = &outgoing_particles_[1];
509 
510  const double mass_a = masses.first;
511  const double mass_b = masses.second;
512 
513  const std::array<double, 2> t_range = get_t_range<double>(
514  kinetic_energy_cm, mass_in_a, mass_in_b, mass_a, mass_b);
515  Angles phitheta;
516  if (el_scattering && !isotropic_) {
520  double mandelstam_s_new = 0.;
521  if (nn_scattering) {
522  mandelstam_s_new = mandelstam_s();
523  } else {
524  /* In the case of elastic collisions other than NN collisions,
525  * there is an ambiguity on how to get the lab-frame momentum (plab),
526  * since the incoming particles can have different masses.
527  * Right now, we first obtain the center-of-mass momentum
528  * of the collision (pcom_now).
529  * Then, the lab-frame momentum is evaluated from the mandelstam s,
530  * which yields the original center-of-mass momentum
531  * when nucleon mass is assumed. */
532  const double pcm_now = pCM_from_s(mandelstam_s(), mass_in_a, mass_in_b);
533  mandelstam_s_new =
534  4. * std::sqrt(pcm_now * pcm_now + nucleon_mass * nucleon_mass);
535  }
536  double bb, a, plab = plab_from_s(mandelstam_s_new);
537  if (nn_scattering &&
538  p_a->pdgcode().antiparticle_sign() ==
539  p_b->pdgcode().antiparticle_sign() &&
540  std::abs(p_a->type().charge() + p_b->type().charge()) == 1) {
541  // proton-neutron and antiproton-antineutron
542  bb = std::max(Cugnon_bnp(plab), really_small);
543  a = (plab < 0.8) ? 1. : 0.64 / (plab * plab);
544  } else {
545  /* all others including pp, nn and AQM elastic processes
546  * This is applied for all particle pairs, which are allowed to
547  * interact elastically. */
548  bb = std::max(Cugnon_bpp(plab), really_small);
549  a = 1.;
550  }
551  double t = random::expo(bb, t_range[0], t_range[1]);
552  if (random::canonical() > 1. / (1. + a)) {
553  t = t_range[0] + t_range[1] - t;
554  }
555  // determine scattering angles in center-of-mass frame
556  phitheta = Angles(2. * M_PI * random::canonical(),
557  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
558  } else if (nn_scattering && p_a->pdgcode().is_Delta() &&
559  p_b->pdgcode().is_nucleon() &&
560  p_a->pdgcode().antiparticle_sign() ==
561  p_b->pdgcode().antiparticle_sign() &&
562  !isotropic_) {
566  const double plab = plab_from_s(mandelstam_s());
567  const double bb = std::max(Cugnon_bpp(plab), really_small);
568  double t = random::expo(bb, t_range[0], t_range[1]);
569  if (random::canonical() > 0.5) {
570  t = t_range[0] + t_range[1] - t; // symmetrize
571  }
572  phitheta = Angles(2. * M_PI * random::canonical(),
573  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
574  } else if (nn_scattering && p_b->pdgcode().is_nucleon() && !isotropic_ &&
575  (p_a->type().is_Nstar() || p_a->type().is_Deltastar())) {
577  const std::array<double, 4> p{1.46434, 5.80311, -6.89358, 1.94302};
578  const double a = p[0] + mass_a * (p[1] + mass_a * (p[2] + mass_a * p[3]));
579  /* If the resonance is so heavy that the index "a" exceeds 30,
580  * the power function turns out to be too sharp. Take t directly to be
581  * t_0 in such a case. */
582  double t = t_range[0];
583  if (a < 30) {
584  t = random::power(-a, t_range[0], t_range[1]);
585  }
586  if (random::canonical() > 0.5) {
587  t = t_range[0] + t_range[1] - t; // symmetrize
588  }
589  phitheta = Angles(2. * M_PI * random::canonical(),
590  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
591  } else {
592  /* isotropic angular distribution */
593  phitheta.distribute_isotropically();
594  }
595 
596  ThreeVector pscatt = phitheta.threevec();
597  // 3-momentum of first incoming particle in center-of-mass frame
598  ThreeVector pcm =
599  incoming_particles_[0].momentum().lorentz_boost(beta_cm()).threevec();
600  pscatt.rotate_z_axis_to(pcm);
601 
602  // final-state CM momentum
603  const double p_f = pCM(kinetic_energy_cm, mass_a, mass_b);
604  if (!(p_f > 0.0)) {
605  logg[LScatterAction].warn("Particle: ", p_a->pdgcode(),
606  " radial momentum: ", p_f);
607  logg[LScatterAction].warn("Etot: ", kinetic_energy_cm, " m_a: ", mass_a,
608  " m_b: ", mass_b);
609  }
610  p_a->set_4momentum(mass_a, pscatt * p_f);
611  p_b->set_4momentum(mass_b, -pscatt * p_f);
612 
613  /* Debug message is printed before boost, so that p_a and p_b are
614  * the momenta in the center of mass frame and thus opposite to
615  * each other.*/
616  logg[LScatterAction].debug("p_a: ", *p_a, "\np_b: ", *p_b);
617 }
smash::random::power
T power(T n, T xMin, T xMax)
Draws a random number according to a power-law distribution ~ x^n.
Definition: random.h:203
smash::random::expo
T expo(T A, T x1, T x2)
Draws a random number x from an exponential distribution exp(A*x), where A is assumed to be positive,...
Definition: random.h:166
smash::random::canonical
T canonical()
Definition: random.h:113
smash::plab_from_s
double plab_from_s(double mandelstam_s, double mass)
Convert Mandelstam-s to p_lab in a fixed-target collision.
Definition: kinematics.h:157
smash::Cugnon_bnp
static double Cugnon_bnp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic np...
Definition: scatteraction.cc:474
smash::pCM
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
smash::nucleon_mass
constexpr double nucleon_mass
Nucleon mass in GeV.
Definition: constants.h:62
smash::pCM_from_s
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:66
smash::is_string_soft_process
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
Definition: processbranch.cc:18
smash::Cugnon_bpp
static double Cugnon_bpp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic pp...
Definition: scatteraction.cc:455

◆ 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_parametersparameters for collision finding.

Definition at line 139 of file scatteraction.cc.

140  {
141  if (were_processes_added_) {
142  logg[LScatterAction].fatal() << "Trying to add processes again.";
143  throw std::logic_error(
144  "add_all_scatterings should be called only once per ScatterAction "
145  "instance");
146  } else {
147  were_processes_added_ = true;
148  }
149  CrossSections xs(incoming_particles_, sqrt_s(),
150  get_potential_at_interaction_point());
151  CollisionBranchList processes =
152  xs.generate_collision_list(finder_parameters, string_process_);
153 
154  /* Add various subprocesses.*/
155  add_collisions(std::move(processes));
156 
157  /* If the string processes are not triggered by a probability, then they
158  * always happen as long as the parametrized total cross section is larger
159  * than the sum of the cross sections of the non-string processes, and the
160  * square root s exceeds the threshold by at least 0.9 GeV. The cross section
161  * of the string processes are counted by taking the difference between the
162  * parametrized total and the sum of the non-strings. */
163  if (!finder_parameters.strings_with_probability &&
164  xs.string_probability(finder_parameters) > 0) {
165  const double xs_diff =
166  xs.high_energy(finder_parameters) - sum_of_partial_cross_sections_;
167  if (xs_diff > 0.) {
168  add_collisions(
169  xs.string_excitation(xs_diff, string_process_, finder_parameters));
170  }
171  }
172 
173  ParticleTypePtr pseudoresonance =
174  try_find_pseudoresonance(finder_parameters.pseudoresonance_method,
175  finder_parameters.transition_high_energy);
176  if (pseudoresonance && finder_parameters.two_to_one) {
177  const double xs_total = is_total_parametrized_
178  ? *parametrized_total_cross_section_
179  : xs.high_energy(finder_parameters);
180  const double xs_gap = xs_total - sum_of_partial_cross_sections_;
181  // The pseudo-resonance is only created if there is a (positive) cross
182  // section gap
183  if (xs_gap > really_small) {
184  auto pseudoresonance_branch = std::make_unique<CollisionBranch>(
185  *pseudoresonance, xs_gap, ProcessType::TwoToOne);
186  add_collision(std::move(pseudoresonance_branch));
187  logg[LScatterAction].debug()
188  << "Pseudoresonance between " << incoming_particles_[0].type().name()
189  << " and " << incoming_particles_[1].type().name() << "is"
190  << pseudoresonance->name() << " with cross section " << xs_gap
191  << "mb.";
192  }
193  }
194  were_processes_added_ = true;
195  // Rescale the branches so that their sum matches the parametrization
196  if (is_total_parametrized_) {
197  rescale_outgoing_branches();
198  }
199 }
smash::Action::get_potential_at_interaction_point
std::pair< FourVector, FourVector > get_potential_at_interaction_point() const
Get the skyrme and asymmetry potential at the interaction point.
Definition: action.cc:112
smash::Action::sqrt_s
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:271
smash::ScatterAction::rescale_outgoing_branches
void rescale_outgoing_branches()
Loop over the possible branches and rescales their weight according to the desired total cross sectio...
Definition: scatteraction.cc:201
smash::ScatterAction::add_collision
void add_collision(CollisionBranchPtr p)
Add a new collision channel.
Definition: scatteraction.cc:47
smash::ScatterAction::string_process_
StringProcess * string_process_
Pointer to interface class for strings.
Definition: scatteraction.h:329
smash::ScatterAction::add_collisions
void add_collisions(CollisionBranchList pv)
Add several new collision channels at once.
Definition: scatteraction.cc:52
smash::ScatterAction::were_processes_added_
bool were_processes_added_
Lock for calling add_all_scatterings only once.
Definition: scatteraction.h:341
smash::ScatterAction::try_find_pseudoresonance
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.
Definition: scatteraction.cc:239

◆ 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_parametersParameters for collision finding.

Definition at line 307 of file scatteraction.cc.

308  {
309  CrossSections xs(incoming_particles_, sqrt_s(),
310  get_potential_at_interaction_point());
311 
312  if (is_total_parametrized_) {
313  parametrized_total_cross_section_ =
314  xs.parametrized_total(finder_parameters);
315  } else {
316  logg[LScatterAction].fatal()
317  << "Trying to parametrize total cross section when it shouldn't be.";
318  throw std::logic_error(
319  "This function can only be called on ScatterAction objects with "
320  "parametrized cross section.");
321  }
322 }

◆ collision_channels()

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

Get list of possible collision channels.

Returns
list of possible collision channels.

Definition at line 167 of file scatteraction.h.

167  {
168  return collision_channels_;
169  }

◆ 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_procString process object to be used.

Definition at line 187 of file scatteraction.h.

187  {
188  string_process_ = str_proc;
189  }

◆ 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 197 of file scatteraction.h.

197  {
198  if (is_total_parametrized_) {
199  return *parametrized_total_cross_section_;
200  }
201  return sum_of_partial_cross_sections_;
202  }

◆ 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 345 of file scatteraction.cc.

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

◆ 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 351 of file scatteraction.cc.

351  {
352  const double m1 = incoming_particles_[0].effective_mass();
353  const double m2 = incoming_particles_[1].effective_mass();
354  return pCM_sqr(sqrt_s(), m1, m2);
355 }
smash::pCM_sqr
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91

◆ 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 335 of file scatteraction.cc.

335  {
336  return total_momentum().velocity();
337 }
smash::FourVector::velocity
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
Definition: fourvector.h:333

◆ 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 339 of file scatteraction.cc.

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

◆ elastic_scattering()

void smash::ScatterAction::elastic_scattering ( )
protected

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

Definition at line 619 of file scatteraction.cc.

619  {
620  // copy initial particles into final state
621  outgoing_particles_[0] = incoming_particles_[0];
622  outgoing_particles_[1] = incoming_particles_[1];
623  // resample momenta
624  sample_angles({outgoing_particles_[0].effective_mass(),
625  outgoing_particles_[1].effective_mass()},
626  sqrt_s());
627 }
smash::ScatterAction::sample_angles
void sample_angles(std::pair< double, double > masses, double kinetic_energy_cm) override
Sample final-state angles in a 2->2 collision (possibly anisotropic).
Definition: scatteraction.cc:486

◆ inelastic_scattering()

void smash::ScatterAction::inelastic_scattering ( )
protected

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

Definition at line 629 of file scatteraction.cc.

629  {
630  // create new particles
631  sample_2body_phasespace();
632  assign_formation_time_to_outgoing_particles();
633  if (spin_interaction_type_ != SpinInteractionType::Off) {
634  assign_unpolarized_spin_vector_to_outgoing_particles();
635  }
636 }
smash::Action::sample_2body_phasespace
virtual void sample_2body_phasespace()
Sample the full 2-body phase-space (masses, momenta, angles) in the center-of-mass frame for the fina...
Definition: action.cc:305
smash::Action::assign_formation_time_to_outgoing_particles
void assign_formation_time_to_outgoing_particles()
Assign the formation time to the outgoing particles.
Definition: action.cc:191
smash::Action::assign_unpolarized_spin_vector_to_outgoing_particles
void assign_unpolarized_spin_vector_to_outgoing_particles()
Assign an unpolarized spin vector to all outgoing particles.
Definition: action.cc:478
SpinInteractionType::Off
@ Off
No spin interactions.

◆ 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 638 of file scatteraction.cc.

638  {
639  sample_manybody_phasespace();
640  assign_formation_time_to_outgoing_particles();
641  if (spin_interaction_type_ != SpinInteractionType::Off) {
642  assign_unpolarized_spin_vector_to_outgoing_particles();
643  }
644  logg[LScatterAction].debug("2->", outgoing_particles_.size(),
645  " scattering:", incoming_particles_, " -> ",
646  outgoing_particles_);
647 }
smash::Action::sample_manybody_phasespace
virtual void sample_manybody_phasespace()
Sample the full n-body phase-space (masses, momenta, angles) in the center-of-mass frame for the fina...
Definition: action.cc:449

◆ 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 674 of file scatteraction.cc.

674  {
675  outgoing_particles_ = string_process_->get_final_state();
676  assign_formation_time_to_outgoing_particles();
677  /* Check momentum difference for debugging */
678  FourVector out_mom;
679  for (ParticleData data : outgoing_particles_) {
680  out_mom += data.momentum();
681  }
682  logg[LPythia].debug("Incoming momenta string:", total_momentum());
683  logg[LPythia].debug("Outgoing momenta string:", out_mom);
684 }
smash::StringProcess::get_final_state
ParticleList get_final_state()
a function to get the final state particle list which is called after the collision
Definition: stringprocess.h:783
smash::LPythia
static constexpr int LPythia
Definition: stringprocess.h:26

◆ 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 690 of file scatteraction.cc.

690  {
691  assert(incoming_particles_.size() == 2);
692  // Disable floating point exception trap for Pythia
693  {
694  DisableFloatTraps guard;
695  /* initialize the string_process_ object for this particular collision */
696  string_process_->init(incoming_particles_, time_of_execution_);
697  /* implement collision */
698  bool success = false;
699  int ntry = 0;
700  const int ntry_max = 10000;
701  while (!success && ntry < ntry_max) {
702  ntry++;
703  switch (process_type_) {
704  case ProcessType::StringSoftSingleDiffractiveAX:
705  /* single diffractive to A+X */
706  success = string_process_->next_SDiff(true);
707  break;
708  case ProcessType::StringSoftSingleDiffractiveXB:
709  /* single diffractive to X+B */
710  success = string_process_->next_SDiff(false);
711  break;
712  case ProcessType::StringSoftDoubleDiffractive:
713  /* double diffractive */
714  success = string_process_->next_DDiff();
715  break;
716  case ProcessType::StringSoftNonDiffractive:
717  /* soft non-diffractive */
718  success = string_process_->next_NDiffSoft();
719  break;
720  case ProcessType::StringSoftAnnihilation:
721  /* soft BBbar 2 mesonic annihilation */
722  success = string_process_->next_BBbarAnn();
723  break;
724  case ProcessType::StringHard:
725  success = string_process_->next_NDiffHard();
726  break;
727  default:
728  logg[LPythia].error("Unknown string process required.");
729  success = false;
730  }
731  }
732  if (ntry == ntry_max) {
733  /* If pythia fails to form a string, it is usually because the energy
734  * is not large enough. In this case, annihilation is then enforced. If
735  * this process still does not not produce any results, it defaults to
736  * an elastic collision. */
737  bool success_newtry = false;
738 
739  /* Check if the initial state is a baryon-antibaryon state.*/
740  PdgCode part1 = incoming_particles_[0].pdgcode(),
741  part2 = incoming_particles_[1].pdgcode();
742  bool is_BBbar_Pair = (part1.baryon_number() != 0) &&
743  (part1.baryon_number() == -part2.baryon_number());
744 
745  /* Decide on the new process .*/
746  if (is_BBbar_Pair) {
747  process_type_ = ProcessType::StringSoftAnnihilation;
748  } else {
749  process_type_ = ProcessType::StringSoftDoubleDiffractive;
750  }
751  /* Perform the new process*/
752  int ntry_new = 0;
753  while (!success_newtry && ntry_new < ntry_max) {
754  ntry_new++;
755  if (is_BBbar_Pair) {
756  success_newtry = string_process_->next_BBbarAnn();
757  } else {
758  success_newtry = string_process_->next_DDiff();
759  }
760  }
761 
762  if (success_newtry) {
763  create_string_final_state();
764  }
765 
766  if (!success_newtry) {
767  /* If annihilation fails:
768  * Particles are normally added after process selection for
769  * strings, outgoing_particles is still uninitialized, and memory
770  * needs to be allocated. We also shift the process_type_ to elastic
771  * so that sample_angles does a proper treatment. */
772  outgoing_particles_.reserve(2);
773  outgoing_particles_.push_back(ParticleData{incoming_particles_[0]});
774  outgoing_particles_.push_back(ParticleData{incoming_particles_[1]});
775  process_type_ = ProcessType::FailedString;
776  elastic_scattering();
777  }
778  } else {
779  create_string_final_state();
780  }
781  }
782 }
smash::Action::time_of_execution_
const double time_of_execution_
Time at which the action is supposed to be performed (absolute time in the lab frame in fm).
Definition: action.h:378
smash::ScatterAction::create_string_final_state
void create_string_final_state()
Creates the final states for string-processes after they are performed.
Definition: scatteraction.cc:674
smash::StringProcess::next_SDiff
bool next_SDiff(bool is_AB_to_AX)
Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns .
Definition: stringprocess.cc:242
smash::StringProcess::next_NDiffSoft
bool next_NDiffSoft()
Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig ...
Definition: stringprocess.cc:451
smash::StringProcess::next_DDiff
bool next_DDiff()
Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process,...
Definition: stringprocess.cc:396
smash::StringProcess::next_BBbarAnn
bool next_BBbarAnn()
Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca , Bleicher:1999xi does,...
Definition: stringprocess.cc:1527
smash::StringProcess::init
void init(const ParticleList &incoming, double tcoll)
initialization feed intial particles, time of collision and gamma factor of the center of mass.
Definition: stringprocess.cc:213
smash::StringProcess::next_NDiffHard
bool next_NDiffHard()
Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions.
Definition: stringprocess.cc:528
smash::ProcessType::FailedString
@ FailedString
See here for a short description.

◆ spin_interaction()

void smash::ScatterAction::spin_interaction ( )
protected

Perform spin interaction in binary interactions.

At the moment, we include a spin-flip in the y component of elastic scatterings if enabled.

Definition at line 795 of file scatteraction.cc.

795  {
796  if (spin_interaction_type_ != SpinInteractionType::Off) {
797  /* 2->2 elastic scattering */
798  if (process_type_ == ProcessType::Elastic) {
799  // Final boost to the outgoing particle momenta
800  boost_spin_vectors_after_elastic_scattering(outgoing_particles_[0],
801  outgoing_particles_[1]);
802  }
803 
804  /* 2->1 resonance formation */
805  if (process_type_ == ProcessType::TwoToOne) {
806  /*
807  * @brief Λ+π → Σ* resonance formation with Λ–spin bookkeeping.
808  *
809  * We do not simulate a direct inelastic Λ+π scattering; instead we form a
810  * Σ* resonance and let it decay later. To preserve Λ polarization per
811  * arXiv:2404.15890v2, we treat the Σ* spin vector as a proxy for the
812  * would-be outgoing Λ spin: at formation, we set the Σ* spin to the
813  * incoming Λ spin and apply a possible spin flip according to the
814  * Λ–flip/non-flip fractions extracted from the paper. On Σ* → Λ+π decay,
815  * the Σ* spin vector is copied to the Λ, thus transporting Λ polarization
816  * through the resonance stage.
817  */
818  // Identify if the outgoing resonance is a Σ*
819  if (outgoing_particles_[0].is_sigmastar()) {
820  // Check that one of the incoming particles is a Λ and the other a π
821  const bool has_lambda = incoming_particles_[0].pdgcode().is_Lambda() ||
822  incoming_particles_[1].pdgcode().is_Lambda();
823  const bool has_pion = incoming_particles_[0].is_pion() ||
824  incoming_particles_[1].is_pion();
825  if (has_lambda && has_pion) {
826  auto &lambda = (incoming_particles_[0].pdgcode().is_Lambda())
827  ? incoming_particles_[0]
828  : incoming_particles_[1];
829  auto &sigma_star = outgoing_particles_[0];
830 
831  // Perform spin flip with probability of 2/9
832  int random_int = random::uniform_int(1, 9);
833  FourVector final_spin_vector = lambda.spin_vector();
834 
835  if (random_int <= 7) {
836  // No spin flip
837  final_spin_vector = final_spin_vector.lorentz_boost(
838  outgoing_particles_[0].velocity());
839  outgoing_particles_[0].set_spin_vector(final_spin_vector);
840  } else {
841  // Spin flip in Lambda rest frame
842  ThreeVector lambda_velocity = lambda.velocity();
843  final_spin_vector =
844  final_spin_vector.lorentz_boost(lambda_velocity);
845 
846  // Flip the spatial spin vector components
847  final_spin_vector[1] = -final_spin_vector[1];
848  final_spin_vector[2] = -final_spin_vector[2];
849  final_spin_vector[3] = -final_spin_vector[3];
850 
851  // Boost back to computational frame and to Sigma* frame
852  final_spin_vector =
853  final_spin_vector.lorentz_boost(-lambda_velocity);
854  final_spin_vector =
855  final_spin_vector.lorentz_boost(sigma_star.velocity());
856  sigma_star.set_spin_vector(final_spin_vector);
857  }
858  }
859  }
860  }
861  }
862 }
smash::random::uniform_int
T uniform_int(T min, T max)
Definition: random.h:100
smash::boost_spin_vectors_after_elastic_scattering
static void boost_spin_vectors_after_elastic_scattering(ParticleData &outgoing_particle_a, ParticleData &outgoing_particle_b)
Definition: scatteraction.cc:784

◆ string_spin_interaction()

void smash::ScatterAction::string_spin_interaction ( )
protected

Perform spin interaction in string excitations.

At the moment, we assign unpolarized spin vectors to the outgoing particles, unless the process is single diffractive, in which case we copy the spin vector of the elastically scattered particle to the corresponding outgoing particle.

Definition at line 864 of file scatteraction.cc.

864  {
865  // Check if spin interaction is disabled
866  if (spin_interaction_type_ == SpinInteractionType::Off) {
867  return;
868  }
869  const bool is_AB_to_AX =
870  (process_type_ == ProcessType::StringSoftSingleDiffractiveAX);
871  const bool is_AB_to_XB =
872  (process_type_ == ProcessType::StringSoftSingleDiffractiveXB);
873 
874  /* This logic relies on the assumption that the surviving hadron is
875  * always appended as the final element in the outgoing particle list.
876  * This ordering is guaranteed by StringProcess::next_SDiff(bool
877  * is_AB_to_AX). If that implementation changes, the behavior here must
878  * be re-evaluated. */
879  if (is_AB_to_AX || is_AB_to_XB) {
880  const std::size_t idx_hadron_in = is_AB_to_AX ? 0 : 1;
881 
882  // Boost spin vector of surviving hadron to outgoing frame
883  const FourVector final_spin_vector =
884  incoming_particles_[idx_hadron_in].spin_vector().lorentz_boost(
885  outgoing_particles_.back().velocity());
886 
887  outgoing_particles_.back().set_spin_vector(final_spin_vector);
888 
889  /* Set unpolarized spin vector for all newly created particles (all but
890  * the last one) */
891  for (auto it = outgoing_particles_.begin();
892  it != outgoing_particles_.end() - 1; ++it) {
893  it->set_unpolarized_spin_vector();
894  }
895  } else {
896  for (auto &particle : outgoing_particles_) {
897  particle.set_unpolarized_spin_vector();
898  }
899  }
900 }

◆ 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 649 of file scatteraction.cc.

649  {
650  if (outgoing_particles_.size() != 1) {
651  std::string s =
652  "resonance_formation: "
653  "Incorrect number of particles in final state: ";
654  s += std::to_string(outgoing_particles_.size()) + " (";
655  s += incoming_particles_[0].pdgcode().string() + " + ";
656  s += incoming_particles_[1].pdgcode().string() + ")";
657  throw InvalidResonanceFormation(s);
658  }
659  // Set the momentum of the formed resonance in its rest frame.
660  outgoing_particles_[0].set_4momentum(
661  total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
662  assign_formation_time_to_outgoing_particles();
663  if (spin_interaction_type_ != SpinInteractionType::Off) {
664  assign_unpolarized_spin_vector_to_outgoing_particles();
665  }
666  /* this momentum is evaluated in the computational frame. */
667  logg[LScatterAction].debug("Momentum of the new particle: ",
668  outgoing_particles_[0].momentum());
669 }

◆ 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 201 of file scatteraction.cc.

201  {
202  if (!were_processes_added_) {
203  logg[LScatterAction].fatal()
204  << "Trying to rescale branches before adding processes.";
205  throw std::logic_error(
206  "This function can only be called after having added processes.");
207  }
208  if (sum_of_partial_cross_sections_ < really_small) {
209  const ParticleTypePtr type_a = &incoming_particles_[0].type();
210  const ParticleTypePtr type_b = &incoming_particles_[1].type();
211  // This is a std::set instead of std::pair because the order of particles
212  // does not matter here
213  const std::set<ParticleTypePtr> pair{type_a, type_b};
214  if (!warned_no_rescaling_available.count(pair)) {
215  logg[LScatterAction].warn()
216  << "Total cross section between " << type_a->name() << " and "
217  << type_b->name() << "is roughly zero at sqrt(s) = " << sqrt_s()
218  << " GeV, and no rescaling to match the parametrized value will be "
219  "done.\nAn elastic process will be added, instead, to match the "
220  "total cross section.\nFor this pair of particles, this warning "
221  "will be subsequently suppressed.";
222  warned_no_rescaling_available.insert(pair);
223  }
224  auto elastic_branch = std::make_unique<CollisionBranch>(
225  *type_a, *type_b, *parametrized_total_cross_section_,
226  ProcessType::Elastic);
227  add_collision(std::move(elastic_branch));
228  } else {
229  const double reweight =
230  *parametrized_total_cross_section_ / sum_of_partial_cross_sections_;
231  logg[LScatterAction].debug("Reweighting ", sum_of_partial_cross_sections_,
232  " to ", *parametrized_total_cross_section_);
233  for (auto &proc : collision_channels_) {
234  proc->set_weight(proc->weight() * reweight);
235  }
236  }
237 }
smash::ScatterAction::warned_no_rescaling_available
static std::set< std::set< ParticleTypePtr > > warned_no_rescaling_available
Warn about zero cross section only once per particle type pair.
Definition: scatteraction.h:345

◆ 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]methodused to select the pseudo-resonance among possible candidates. See user guide description for more information.
[in]transitionparameters 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 239 of file scatteraction.cc.

241  {
242  const ParticleTypePtr type_a = &incoming_particles_[0].type();
243  const ParticleTypePtr type_b = &incoming_particles_[1].type();
244  const double desired_mass = sqrt_s();
245 
246  double string_offset = 0.5 * transition.sqrts_add_lower;
247  const bool nucleon_and_pion = (type_a->is_nucleon() && type_b->is_pion()) ||
248  (type_a->is_pion() && type_b->is_nucleon());
249  const bool two_nucleons = type_a->is_nucleon() && type_b->is_nucleon();
250  const bool two_pions = type_a->is_pion() && type_b->is_pion();
251  const bool nucleon_and_kaon = (type_a->is_nucleon() && type_b->is_kaon()) ||
252  (type_a->is_kaon() && type_b->is_nucleon());
253  if (nucleon_and_pion) {
254  string_offset = transition.sqrts_range_Npi.first - pion_mass - nucleon_mass;
255  } else if (two_nucleons) {
256  string_offset = transition.sqrts_range_NN.first - 2 * nucleon_mass;
257  } else if (two_pions) {
258  string_offset = transition.pipi_offset;
259  } else if (nucleon_and_kaon) {
260  string_offset = transition.KN_offset;
261  }
262  /*
263  * Artificial cutoff to create a pseudo-resonance only close to the string
264  * transition, where data or first-principle models are unhelpful
265  */
266  if (desired_mass < incoming_particles_[0].effective_mass() +
267  incoming_particles_[1].effective_mass() +
268  string_offset) {
269  return {};
270  }
271 
272  if (method == PseudoResonance::None) {
273  return {};
274  } else if (method == PseudoResonance::LargestFromUnstable ||
275  method == PseudoResonance::ClosestFromUnstable) {
276  if (type_a->is_stable() && type_b->is_stable()) {
277  return {};
278  }
279  }
280 
281  // If this list is empty, there are no possible pseudo-resonances.
282  ParticleTypePtrList list = list_possible_resonances(type_a, type_b);
283  if (std::empty(list)) {
284  return {};
285  }
286 
287  if (method == PseudoResonance::Largest ||
288  method == PseudoResonance::LargestFromUnstable) {
289  auto largest = *std::max_element(list.begin(), list.end(),
290  [](ParticleTypePtr a, ParticleTypePtr b) {
291  return a->mass() < b->mass();
292  });
293  return largest;
294  } else if (method == PseudoResonance::Closest ||
295  method == PseudoResonance::ClosestFromUnstable) {
296  auto comparison = [&desired_mass](ParticleTypePtr a, ParticleTypePtr b) {
297  return std::abs(a->mass() - desired_mass) <
298  std::abs(b->mass() - desired_mass);
299  };
300  auto closest = *std::min_element(list.begin(), list.end(), comparison);
301  return closest;
302  } else {
303  throw std::logic_error("Unknown method for selecting pseudoresonance.");
304  }
305 }
PseudoResonance::Closest
@ Closest
Resonance with the pole mass closest from the invariant mass of incoming particles for all processes.
PseudoResonance::ClosestFromUnstable
@ ClosestFromUnstable
Closest resonance for a given mass from processes with at least one resonance in the incoming particl...
PseudoResonance::None
@ None
No pseudo-resonance is created.
PseudoResonance::LargestFromUnstable
@ LargestFromUnstable
Heaviest possible resonance from processes with at least one resonance in the incoming particles.
PseudoResonance::Largest
@ Largest
Resonance of largest mass for all processes.
smash::list_possible_resonances
ParticleTypePtrList list_possible_resonances(const ParticleTypePtr type_a, const ParticleTypePtr type_b)
Lists the possible resonances that decay into two particles.
Definition: particletype.cc:776
smash::pion_mass
constexpr double pion_mass
Pion mass in GeV.
Definition: constants.h:69

Member Data Documentation

◆ collision_channels_

CollisionBranchList smash::ScatterAction::collision_channels_
protected

List of possible collisions.

Definition at line 277 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 280 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 283 of file scatteraction.h.

◆ isotropic_

bool smash::ScatterAction::isotropic_ = false
protected

Do this collision isotropically?

Definition at line 286 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 289 of file scatteraction.h.

◆ string_process_

StringProcess* smash::ScatterAction::string_process_ = nullptr
private

Pointer to interface class for strings.

Definition at line 329 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 332 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 335 of file scatteraction.h.

◆ spin_interaction_type_

SpinInteractionType smash::ScatterAction::spin_interaction_type_ = SpinInteractionType::Off
private

What kind of spin interaction to use.

Definition at line 338 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 341 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 345 of file scatteraction.h.

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