Version: SMASH-2.2
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:
[legend]
Collaboration diagram for smash::ScatterAction:
[legend]

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)
 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 (double elastic_parameter, bool two_to_one, ReactionsBitSet included_2to2, MultiParticleReactionsBitSet included_multi, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment, double scale_xs, double additional_el_xs)
 Add all possible scattering subprocesses for this action object. 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. 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 void 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 void 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, FourVectorget_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 total_cross_section_
 Total hadronic cross section. 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/c). 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...
 

Private Attributes

StringProcessstring_process_ = nullptr
 Pointer to interface class for strings. 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 
)

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.

Definition at line 29 of file scatteraction.cc.

33  : Action({in_part_a, in_part_b}, time),
35  isotropic_(isotropic),
36  string_formation_time_(string_formation_time) {
37  box_length_ = box_length;
38 }
Action(const ParticleList &in_part, double time)
Construct an action object with incoming particles and relative time.
Definition: action.h:44
double box_length_
Box length: needed to determine coordinates of collision correctly in case of collision through the w...
Definition: action.h:364
bool isotropic_
Do this collision isotropically?
double string_formation_time_
Time fragments take to be fully formed in hard string excitation.
double total_cross_section_
Total hadronic cross section.

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

40  {
41  add_process<CollisionBranch>(p, collision_channels_, total_cross_section_);
42 }
CollisionBranchList collision_channels_
List of possible collisions.
constexpr int p
Proton.
Here is the caller graph for this function:

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

44  {
45  add_processes<CollisionBranch>(std::move(pv), collision_channels_,
47 }
Here is the caller graph for this function:

◆ 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 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_\mathrm{coll} = (\vec{x_a} - \vec{x_b})^2 - \frac{((\vec{x_a} - \vec{x_b}) \cdot (\vec{p_a} - \vec{p_b}))^2 } {(\vec{p_a} - \vec{p_b})^2}\]

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

UrQMD squared distance criterion: Bass:1998ca [4] (3.27): in 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

Definition at line 184 of file scatteraction.cc.

184  {
185  // local copy of particles (since we need to boost them)
186  ParticleData p_a = incoming_particles_[0];
187  ParticleData p_b = incoming_particles_[1];
188  /* Boost particles to center-of-momentum frame. */
189  const ThreeVector velocity = beta_cm();
190  p_a.boost(velocity);
191  p_b.boost(velocity);
192  const ThreeVector pos_diff =
193  p_a.position().threevec() - p_b.position().threevec();
194  const ThreeVector mom_diff =
195  p_a.momentum().threevec() - p_b.momentum().threevec();
196 
197  logg[LScatterAction].debug("Particle ", incoming_particles_,
198  " position difference [fm]: ", pos_diff,
199  ", momentum difference [GeV]: ", mom_diff);
200 
201  const double dp2 = mom_diff.sqr();
202  const double dr2 = pos_diff.sqr();
203  /* Zero momentum leads to infite distance. */
204  if (dp2 < really_small) {
205  return dr2;
206  }
207  const double dpdr = pos_diff * mom_diff;
208 
217  const double result = dr2 - dpdr * dpdr / dp2;
218  return result > 0.0 ? result : 0.0;
219 }
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:340
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering/incoming particles in the calculation frame.
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:37
static constexpr int LScatterAction
Here is the call graph for this function:

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

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

Definition at line 221 of file scatteraction.cc.

221  {
222  // local copy of particles (since we need to boost them)
223  ParticleData p_a = incoming_particles_[0];
224  ParticleData p_b = incoming_particles_[1];
225 
226  const FourVector delta_x = p_a.position() - p_b.position();
227  const double mom_diff_sqr =
228  (p_a.momentum().threevec() - p_b.momentum().threevec()).sqr();
229  const double x_sqr = delta_x.sqr();
230 
231  if (mom_diff_sqr < really_small) {
232  return -x_sqr;
233  }
234 
235  const double p_a_sqr = p_a.momentum().sqr();
236  const double p_b_sqr = p_b.momentum().sqr();
237  const double p_a_dot_x = p_a.momentum().Dot(delta_x);
238  const double p_b_dot_x = p_b.momentum().Dot(delta_x);
239  const double p_a_dot_p_b = p_a.momentum().Dot(p_b.momentum());
240 
241  const double b_sqr =
242  -x_sqr -
243  (p_a_sqr * std::pow(p_b_dot_x, 2) + p_b_sqr * std::pow(p_a_dot_x, 2) -
244  2 * p_a_dot_p_b * p_a_dot_x * p_b_dot_x) /
245  (std::pow(p_a_dot_p_b, 2) - p_a_sqr * p_b_sqr);
246  return b_sqr > 0.0 ? b_sqr : 0.0;
247 }
Here is the call graph for this function:

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

161 { return total_momentum().sqr(); }
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:374
double sqr() const
calculate the square of the vector (which is a scalar)
Definition: fourvector.h:455
Here is the call graph for this function:
Here is the caller graph for this function:

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

Returns
relative velocity.

Definition at line 175 of file scatteraction.cc.

175  {
176  const double m1 = incoming_particles()[0].effective_mass();
177  const double m2 = incoming_particles()[1].effective_mass();
178  const double m_s = mandelstam_s();
179  const double lamb = lambda_tilde(m_s, m1 * m1, m2 * m2);
180  return std::sqrt(lamb) / (2. * incoming_particles()[0].momentum().x0() *
181  incoming_particles()[1].momentum().x0());
182 }
const ParticleList & incoming_particles() const
Get the list of particles that go into the action.
Definition: action.cc:57
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:310
double mandelstam_s() const
Determine the Mandelstam s variable,.
Here is the call graph for this function:

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

49  {
50  logg[LScatterAction].debug("Incoming particles: ", incoming_particles_);
51 
52  /* Decide for a particular final state. */
53  const CollisionBranch *proc = choose_channel<CollisionBranch>(
55  process_type_ = proc->get_type();
56  outgoing_particles_ = proc->particle_list();
57  partial_cross_section_ = proc->weight();
58 
59  logg[LScatterAction].debug("Chosen channel: ", process_type_,
61 
62  /* The production point of the new particles. */
63  FourVector middle_point = get_interaction_point();
64 
65  switch (process_type_) {
67  /* 2->2 elastic scattering */
69  break;
71  /* resonance formation */
73  break;
75  /* 2->2 inelastic scattering */
76  /* Sample the particle momenta in CM system. */
78  break;
82  /* 2->m scattering */
84  break;
92  break;
93  default:
94  throw InvalidScatterAction(
95  "ScatterAction::generate_final_state: Invalid process type " +
96  std::to_string(static_cast<int>(process_type_)) + " was requested. " +
97  "(PDGcode1=" + incoming_particles_[0].pdgcode().string() +
98  ", PDGcode2=" + incoming_particles_[1].pdgcode().string() + ")");
99  }
100 
101  for (ParticleData &new_particle : outgoing_particles_) {
102  // Boost to the computational frame
103  new_particle.boost_momentum(
105  /* Set positions of the outgoing particles */
106  if (proc->get_type() != ProcessType::Elastic) {
107  new_particle.set_4position(middle_point);
108  }
109  }
110 }
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:155
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:348
FourVector get_interaction_point() const
Get the interaction point.
Definition: action.cc:67
ProcessType process_type_
type of process
Definition: action.h:357
void resonance_formation()
Perform a 2->1 resonance-formation process.
double partial_cross_section_
Partial cross-section to the chosen outgoing channel.
void two_to_many_scattering()
Perform an inelastic two-to-many-body scattering (more than 2)
void string_excitation()
Todo(ryu): document better - it is not really UrQMD-based, isn't it? Perform the UrQMD-based string e...
void elastic_scattering()
Perform an elastic two-body scattering, i.e. just exchange momentum.
void inelastic_scattering()
Perform an inelastic two-body scattering, i.e. new particles are formed.
@ TwoToOne
resonance formation (2->1)
@ StringSoftDoubleDiffractive
double diffractive. Two strings are formed, one from A and one from B.
@ TwoToFive
2->5 scattering
@ StringSoftSingleDiffractiveXB
single diffractive AB->XB.
@ TwoToTwo
2->2 inelastic scattering
@ Elastic
elastic scattering: particles remain the same, only momenta change
@ TwoToFour
2->4 scattering
@ StringSoftAnnihilation
a special case of baryon-antibaryon annihilation.
@ StringSoftNonDiffractive
non-diffractive. Two strings are formed both have ends in A and B.
@ StringSoftSingleDiffractiveAX
(41-45) soft string excitations.
@ StringHard
hard string process involving 2->2 QCD process by PYTHIA.
@ TwoToThree
2->3 scattering
Here is the call graph for this function:

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

143  {
144  return total_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
145  incoming_particles_[1].xsec_scaling_factor();
146 }

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

148  {
149  return partial_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
150  incoming_particles_[1].xsec_scaling_factor();
151 }

◆ 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 [15].

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

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

Reimplemented from smash::Action.

Definition at line 293 of file scatteraction.cc.

294  {
297  // We potentially have more than two particles, so the following angular
298  // distributions don't work. Instead we just keep the angular
299  // distributions generated by string fragmentation.
300  return;
301  }
302  assert(outgoing_particles_.size() == 2);
303 
304  // NN scattering is anisotropic currently
305  const bool nn_scattering = incoming_particles_[0].type().is_nucleon() &&
306  incoming_particles_[1].type().is_nucleon();
307  /* Elastic process is anisotropic and
308  * the angular distribution is based on the NN elastic scattering. */
309  const bool el_scattering = process_type_ == ProcessType::Elastic;
310 
311  const double mass_in_a = incoming_particles_[0].effective_mass();
312  const double mass_in_b = incoming_particles_[1].effective_mass();
313 
314  ParticleData *p_a = &outgoing_particles_[0];
315  ParticleData *p_b = &outgoing_particles_[1];
316 
317  const double mass_a = masses.first;
318  const double mass_b = masses.second;
319 
320  const std::array<double, 2> t_range = get_t_range<double>(
321  kinetic_energy_cm, mass_in_a, mass_in_b, mass_a, mass_b);
322  Angles phitheta;
323  if (el_scattering && !isotropic_) {
327  double mandelstam_s_new = 0.;
328  if (nn_scattering) {
329  mandelstam_s_new = mandelstam_s();
330  } else {
331  /* In the case of elastic collisions other than NN collisions,
332  * there is an ambiguity on how to get the lab-frame momentum (plab),
333  * since the incoming particles can have different masses.
334  * Right now, we first obtain the center-of-mass momentum
335  * of the collision (pcom_now).
336  * Then, the lab-frame momentum is evaluated from the mandelstam s,
337  * which yields the original center-of-mass momentum
338  * when nucleon mass is assumed. */
339  const double pcm_now = pCM_from_s(mandelstam_s(), mass_in_a, mass_in_b);
340  mandelstam_s_new =
341  4. * std::sqrt(pcm_now * pcm_now + nucleon_mass * nucleon_mass);
342  }
343  double bb, a, plab = plab_from_s(mandelstam_s_new);
344  if (nn_scattering &&
345  p_a->pdgcode().antiparticle_sign() ==
346  p_b->pdgcode().antiparticle_sign() &&
347  std::abs(p_a->type().charge() + p_b->type().charge()) == 1) {
348  // proton-neutron and antiproton-antineutron
349  bb = std::max(Cugnon_bnp(plab), really_small);
350  a = (plab < 0.8) ? 1. : 0.64 / (plab * plab);
351  } else {
352  /* all others including pp, nn and AQM elastic processes
353  * This is applied for all particle pairs, which are allowed to
354  * interact elastically. */
355  bb = std::max(Cugnon_bpp(plab), really_small);
356  a = 1.;
357  }
358  double t = random::expo(bb, t_range[0], t_range[1]);
359  if (random::canonical() > 1. / (1. + a)) {
360  t = t_range[0] + t_range[1] - t;
361  }
362  // determine scattering angles in center-of-mass frame
363  phitheta = Angles(2. * M_PI * random::canonical(),
364  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
365  } else if (nn_scattering && p_a->pdgcode().is_Delta() &&
366  p_b->pdgcode().is_nucleon() &&
367  p_a->pdgcode().antiparticle_sign() ==
368  p_b->pdgcode().antiparticle_sign() &&
369  !isotropic_) {
373  const double plab = plab_from_s(mandelstam_s());
374  const double bb = std::max(Cugnon_bpp(plab), really_small);
375  double t = random::expo(bb, t_range[0], t_range[1]);
376  if (random::canonical() > 0.5) {
377  t = t_range[0] + t_range[1] - t; // symmetrize
378  }
379  phitheta = Angles(2. * M_PI * random::canonical(),
380  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
381  } else if (nn_scattering && p_b->pdgcode().is_nucleon() && !isotropic_ &&
382  (p_a->type().is_Nstar() || p_a->type().is_Deltastar())) {
384  const std::array<double, 4> p{1.46434, 5.80311, -6.89358, 1.94302};
385  const double a = p[0] + mass_a * (p[1] + mass_a * (p[2] + mass_a * p[3]));
386  /* If the resonance is so heavy that the index "a" exceeds 30,
387  * the power function turns out to be too sharp. Take t directly to be
388  * t_0 in such a case. */
389  double t = t_range[0];
390  if (a < 30) {
391  t = random::power(-a, t_range[0], t_range[1]);
392  }
393  if (random::canonical() > 0.5) {
394  t = t_range[0] + t_range[1] - t; // symmetrize
395  }
396  phitheta = Angles(2. * M_PI * random::canonical(),
397  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
398  } else {
399  /* isotropic angular distribution */
400  phitheta.distribute_isotropically();
401  }
402 
403  ThreeVector pscatt = phitheta.threevec();
404  // 3-momentum of first incoming particle in center-of-mass frame
405  ThreeVector pcm =
406  incoming_particles_[0].momentum().lorentz_boost(beta_cm()).threevec();
407  pscatt.rotate_z_axis_to(pcm);
408 
409  // final-state CM momentum
410  const double p_f = pCM(kinetic_energy_cm, mass_a, mass_b);
411  if (!(p_f > 0.0)) {
412  logg[LScatterAction].warn("Particle: ", p_a->pdgcode(),
413  " radial momentum: ", p_f);
414  logg[LScatterAction].warn("Etot: ", kinetic_energy_cm, " m_a: ", mass_a,
415  " m_b: ", mass_b);
416  }
417  p_a->set_4momentum(mass_a, pscatt * p_f);
418  p_b->set_4momentum(mass_b, -pscatt * p_f);
419 
420  /* Debug message is printed before boost, so that p_a and p_b are
421  * the momenta in the center of mass frame and thus opposite to
422  * each other.*/
423  logg[LScatterAction].debug("p_a: ", *p_a, "\np_b: ", *p_b);
424 }
T power(T n, T xMin, T xMax)
Draws a random number according to a power-law distribution ~ x^n.
Definition: random.h:203
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
T canonical()
Definition: random.h:113
double plab_from_s(double mandelstam_s, double mass)
Convert Mandelstam-s to p_lab in a fixed-target collision.
Definition: kinematics.h:157
static double Cugnon_bnp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic np...
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
constexpr double nucleon_mass
Nucleon mass in GeV.
Definition: constants.h:58
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:66
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
static double Cugnon_bpp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic pp...
Here is the call graph for this function:
Here is the caller graph for this function:

◆ add_all_scatterings()

void smash::ScatterAction::add_all_scatterings ( double  elastic_parameter,
bool  two_to_one,
ReactionsBitSet  included_2to2,
MultiParticleReactionsBitSet  included_multi,
double  low_snn_cut,
bool  strings_switch,
bool  use_AQM,
bool  strings_with_probability,
NNbarTreatment  nnbar_treatment,
double  scale_xs,
double  additional_el_xs 
)

Add all possible scattering subprocesses for this action object.

Parameters
[in]elastic_parameterIf non-zero, given global elastic cross section.
[in]two_to_one2->1 reactions enabled?
[in]included_2to2Which 2->2 reactions are enabled?
[in]included_multiWhich multi-particle reactions are enabled?
[in]low_snn_cutElastic collisions with CME below are forbidden.
[in]strings_switchAre string processes enabled?
[in]use_AQMuse elastic cross sections via AQM?
[in]strings_with_probabilityAre string processes triggered according to a probability?
[in]nnbar_treatmentNNbar treatment through resonance, strings or none
[in]scale_xsFactor by which all (partial) cross sections are scaled
[in]additional_el_xsAdditional constant elastic cross section

Definition at line 112 of file scatteraction.cc.

116  {
117  CrossSections xs(incoming_particles_, sqrt_s(),
119  CollisionBranchList processes = xs.generate_collision_list(
120  elastic_parameter, two_to_one, included_2to2, included_multi, low_snn_cut,
121  strings_switch, use_AQM, strings_with_probability, nnbar_treatment,
122  string_process_, scale_xs, additional_el_xs);
123 
124  /* Add various subprocesses.*/
125  add_collisions(std::move(processes));
126 
127  /* If the string processes are not triggered by a probability, then they
128  * always happen as long as the parametrized total cross section is larger
129  * than the sum of the cross sections of the non-string processes, and the
130  * square root s exceeds the threshold by at least 0.9 GeV. The cross section
131  * of the string processes are counted by taking the difference between the
132  * parametrized total and the sum of the non-strings. */
133  if (!strings_with_probability &&
134  xs.string_probability(strings_switch, strings_with_probability, use_AQM,
135  nnbar_treatment == NNbarTreatment::Strings) == 1.) {
136  const double xs_diff = xs.high_energy() - cross_section();
137  if (xs_diff > 0.) {
138  add_collisions(xs.string_excitation(xs_diff, string_process_, use_AQM));
139  }
140  }
141 }
std::pair< FourVector, FourVector > get_potential_at_interaction_point() const
Get the skyrme and asymmetry potential at the interaction point.
Definition: action.cc:111
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:266
StringProcess * string_process_
Pointer to interface class for strings.
void add_collisions(CollisionBranchList pv)
Add several new collision channels at once.
virtual double cross_section() const
Get the total cross section of the scattering particles.
@ Strings
Use string fragmentation.
Here is the call graph for this function:

◆ 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.

165  {
166  return collision_channels_;
167  }

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

185  {
186  string_process_ = str_proc;
187  }

◆ cross_section()

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

Get the total cross section of the scattering particles.

Returns
total cross section.

Definition at line 194 of file scatteraction.h.

194 { return total_cross_section_; }
Here is the caller graph for this function:

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

163  {
164  const double m1 = incoming_particles_[0].effective_mass();
165  const double m2 = incoming_particles_[1].effective_mass();
166  return pCM(sqrt_s(), m1, m2);
167 }
Here is the call graph for this function:
Here is the caller graph for this function:

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

169  {
170  const double m1 = incoming_particles_[0].effective_mass();
171  const double m2 = incoming_particles_[1].effective_mass();
172  return pCM_sqr(sqrt_s(), m1, m2);
173 }
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91
Here is the call graph for this function:

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

153  {
154  return total_momentum().velocity();
155 }
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
Definition: fourvector.h:328
Here is the call graph for this function:
Here is the caller graph for this function:

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

157  {
158  return (1. / std::sqrt(1.0 - beta_cm().sqr()));
159 }
Here is the call graph for this function:

◆ elastic_scattering()

void smash::ScatterAction::elastic_scattering ( )
protected

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

Definition at line 426 of file scatteraction.cc.

426  {
427  // copy initial particles into final state
430  // resample momenta
431  sample_angles({outgoing_particles_[0].effective_mass(),
432  outgoing_particles_[1].effective_mass()},
433  sqrt_s());
434 }
void sample_angles(std::pair< double, double > masses, double kinetic_energy_cm) override
Sample final-state angles in a 2->2 collision (possibly anisotropic).
Here is the call graph for this function:
Here is the caller graph for this function:

◆ inelastic_scattering()

void smash::ScatterAction::inelastic_scattering ( )
protected

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

Definition at line 436 of file scatteraction.cc.

436  {
437  // create new particles
440 }
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:300
void assign_formation_time_to_outgoing_particles()
Assign the formation time to the outgoing particles.
Definition: action.cc:186
Here is the call graph for this function:
Here is the caller graph for this function:

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

442  {
445  logg[LScatterAction].debug("2->", outgoing_particles_.size(),
446  " scattering:", incoming_particles_, " -> ",
448 }
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:443
Here is the call graph for this function:
Here is the caller graph for this function:

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

472  {
475  /* Check momentum difference for debugging */
476  FourVector out_mom;
477  for (ParticleData data : outgoing_particles_) {
478  out_mom += data.momentum();
479  }
480  logg[LPythia].debug("Incoming momenta string:", total_momentum());
481  logg[LPythia].debug("Outgoing momenta string:", out_mom);
482 }
ParticleList get_final_state()
a function to get the final state particle list which is called after the collision
static constexpr int LPythia
Definition: stringprocess.h:26
Here is the call graph for this function:
Here is the caller graph for this function:

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

488  {
489  assert(incoming_particles_.size() == 2);
490  // Disable floating point exception trap for Pythia
491  {
492  DisableFloatTraps guard;
493  /* initialize the string_process_ object for this particular collision */
495  /* implement collision */
496  bool success = false;
497  int ntry = 0;
498  const int ntry_max = 10000;
499  while (!success && ntry < ntry_max) {
500  ntry++;
501  switch (process_type_) {
503  /* single diffractive to A+X */
504  success = string_process_->next_SDiff(true);
505  break;
507  /* single diffractive to X+B */
508  success = string_process_->next_SDiff(false);
509  break;
511  /* double diffractive */
512  success = string_process_->next_DDiff();
513  break;
515  /* soft non-diffractive */
516  success = string_process_->next_NDiffSoft();
517  break;
519  /* soft BBbar 2 mesonic annihilation */
520  success = string_process_->next_BBbarAnn();
521  break;
523  success = string_process_->next_NDiffHard();
524  break;
525  default:
526  logg[LPythia].error("Unknown string process required.");
527  success = false;
528  }
529  }
530  if (ntry == ntry_max) {
531  /* If pythia fails to form a string, it is usually because the energy
532  * is not large enough. In this case, annihilation is then enforced. If
533  * this process still does not not produce any results, it defaults to
534  * an elastic collision. */
535  bool success_newtry = false;
536 
537  /* Check if the initial state is a baryon-antibaryon state.*/
538  PdgCode part1 = incoming_particles_[0].pdgcode(),
539  part2 = incoming_particles_[1].pdgcode();
540  bool is_BBbar_Pair = (part1.baryon_number() != 0) &&
541  (part1.baryon_number() == -part2.baryon_number());
542 
543  /* Decide on the new process .*/
544  if (is_BBbar_Pair) {
546  } else {
548  }
549  /* Perform the new process*/
550  int ntry_new = 0;
551  while (!success_newtry && ntry_new < ntry_max) {
552  ntry_new++;
553  if (is_BBbar_Pair) {
554  success_newtry = string_process_->next_BBbarAnn();
555  } else {
556  success_newtry = string_process_->next_DDiff();
557  }
558  }
559 
560  if (success_newtry) {
562  }
563 
564  if (!success_newtry) {
565  /* If annihilation fails:
566  * Particles are normally added after process selection for
567  * strings, outgoing_particles is still uninitialized, and memory
568  * needs to be allocated. We also shift the process_type_ to elastic
569  * so that sample_angles does a proper treatment. */
570  outgoing_particles_.reserve(2);
571  outgoing_particles_.push_back(ParticleData{incoming_particles_[0]});
572  outgoing_particles_.push_back(ParticleData{incoming_particles_[1]});
575  }
576  } else {
578  }
579  }
580 }
const double time_of_execution_
Time at which the action is supposed to be performed (absolute time in the lab frame in fm/c).
Definition: action.h:354
void create_string_final_state()
Creates the final states for string-processes after they are performed.
bool next_SDiff(bool is_AB_to_AX)
Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns .
bool next_NDiffSoft()
Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig ...
bool next_DDiff()
Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process,...
bool next_BBbarAnn()
Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca , Bleicher:1999xi does,...
void init(const ParticleList &incoming, double tcoll)
initialization feed intial particles, time of collision and gamma factor of the center of mass.
bool next_NDiffHard()
Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions.
@ FailedString
Soft String NNbar annihilation process can fail by lack of energy.
Here is the call graph for this function:
Here is the caller graph for this function:

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

450  {
451  if (outgoing_particles_.size() != 1) {
452  std::string s =
453  "resonance_formation: "
454  "Incorrect number of particles in final state: ";
455  s += std::to_string(outgoing_particles_.size()) + " (";
456  s += incoming_particles_[0].pdgcode().string() + " + ";
457  s += incoming_particles_[1].pdgcode().string() + ")";
458  throw InvalidResonanceFormation(s);
459  }
460  // Set the momentum of the formed resonance in its rest frame.
461  outgoing_particles_[0].set_4momentum(
462  total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
464  /* this momentum is evaluated in the computational frame. */
465  logg[LScatterAction].debug("Momentum of the new particle: ",
466  outgoing_particles_[0].momentum());
467 }
Here is the call graph for this function:
Here is the caller graph for this function:

Member Data Documentation

◆ collision_channels_

CollisionBranchList smash::ScatterAction::collision_channels_
protected

List of possible collisions.

Definition at line 255 of file scatteraction.h.

◆ total_cross_section_

double smash::ScatterAction::total_cross_section_
protected

Total hadronic cross section.

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

◆ isotropic_

bool smash::ScatterAction::isotropic_ = false
protected

Do this collision isotropically?

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

◆ string_process_

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

Pointer to interface class for strings.

Definition at line 284 of file scatteraction.h.


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