Version: SMASH-1.5
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 31 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)
 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...
 
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, double low_snn_cut, bool strings_switch, bool use_AQM, bool strings_with_probability, NNbarTreatment nnbar_treatment)
 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 Particles &particles, 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...
 
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...
 

Protected Member Functions

double mandelstam_s () const
 Determine the Mandelstam s variable,. More...
 
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 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 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...
 

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

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

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 
)

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

Definition at line 32 of file scatteraction.cc.

35  : Action({in_part_a, in_part_b}, time),
37  isotropic_(isotropic),
38  string_formation_time_(string_formation_time) {}
Action(const ParticleList &in_part, double time)
Construct an action object with incoming particles and relative time.
Definition: action.h:43
double string_formation_time_
Time fragments take to be fully formed in hard string excitation.
double total_cross_section_
Total hadronic cross section.
bool isotropic_
Do this collision isotropically?

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 }
double total_cross_section_
Total hadronic cross section.
constexpr int p
Proton.
CollisionBranchList collision_channels_
List of possible collisions.
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 }
double total_cross_section_
Total hadronic cross section.
CollisionBranchList collision_channels_
List of possible collisions.
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 (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 168 of file scatteraction.cc.

168  {
169  // local copy of particles (since we need to boost them)
170  ParticleData p_a = incoming_particles_[0];
171  ParticleData p_b = incoming_particles_[1];
172  /* Boost particles to center-of-momentum frame. */
173  const ThreeVector velocity = beta_cm();
174  p_a.boost(velocity);
175  p_b.boost(velocity);
176  const ThreeVector pos_diff =
177  p_a.position().threevec() - p_b.position().threevec();
178  const ThreeVector mom_diff =
179  p_a.momentum().threevec() - p_b.momentum().threevec();
180 
181  const auto &log = logger<LogArea::ScatterAction>();
182  log.debug("Particle ", incoming_particles_,
183  " position difference [fm]: ", pos_diff,
184  ", momentum difference [GeV]: ", mom_diff);
185 
186  const double dp2 = mom_diff.sqr();
187  const double dr2 = pos_diff.sqr();
188  /* Zero momentum leads to infite distance. */
189  if (dp2 < really_small) {
190  return dr2;
191  }
192  const double dpdr = pos_diff * mom_diff;
193 
202  return dr2 - dpdr * dpdr / dp2;
203 }
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:34
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering particles in the calculation frame...
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
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  const auto &log = logger<LogArea::ScatterAction>();
51 
52  log.debug("Incoming particles: ", incoming_particles_);
53 
54  /* Decide for a particular final state. */
55  const CollisionBranch *proc = choose_channel<CollisionBranch>(
57  process_type_ = proc->get_type();
58  outgoing_particles_ = proc->particle_list();
59  partial_cross_section_ = proc->weight();
60 
61  log.debug("Chosen channel: ", process_type_, outgoing_particles_);
62 
63  /* The production point of the new particles. */
64  FourVector middle_point = get_interaction_point();
65 
66  switch (process_type_) {
68  /* 2->2 elastic scattering */
70  break;
72  /* resonance formation */
74  break;
76  /* 2->2 inelastic scattering */
77  /* Sample the particle momenta in CM system. */
79  break;
87  break;
88  default:
89  throw InvalidScatterAction(
90  "ScatterAction::generate_final_state: Invalid process type " +
91  std::to_string(static_cast<int>(process_type_)) + " was requested. " +
92  "(PDGcode1=" + incoming_particles_[0].pdgcode().string() +
93  ", PDGcode2=" + incoming_particles_[1].pdgcode().string() + ")");
94  }
95 
96  for (ParticleData &new_particle : outgoing_particles_) {
97  // Boost to the computational frame
98  new_particle.boost_momentum(
100  /* Set positions of the outgoing particles */
101  if (proc->get_type() != ProcessType::Elastic) {
102  new_particle.set_4position(middle_point);
103  }
104  }
105 }
double partial_cross_section_
Partial cross-section to the chosen outgoing channel.
double diffractive. Two strings are formed, one from A and one from B.
FourVector get_interaction_point() const
Get the interaction point.
Definition: action.cc:68
ProcessType process_type_
type of process
Definition: action.h:320
a special case of baryon-antibaryon annihilation.
2->2 inelastic scattering
void resonance_formation()
Perform a 2->1 resonance-formation process.
void string_excitation()
Todo(ryu): document better - it is not really UrQMD-based, isn&#39;t it? Perform the UrQMD-based string e...
double total_cross_section_
Total hadronic cross section.
elastic scattering: particles remain the same, only momenta change
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
hard string process involving 2->2 QCD process by PYTHIA.
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:123
void elastic_scattering()
Perform an elastic two-body scattering, i.e. just exchange momentum.
resonance formation (2->1)
CollisionBranchList collision_channels_
List of possible collisions.
non-diffractive. Two strings are formed both have ends in A and B.
(41-45) soft string excitations.
void inelastic_scattering()
Perform an inelastic two-body scattering, i.e. new particles are formed.
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 136 of file scatteraction.cc.

136  {
137  return total_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
138  incoming_particles_[1].xsec_scaling_factor();
139 }
double total_cross_section_
Total hadronic cross section.
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303

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

141  {
142  return partial_cross_section_ * incoming_particles_[0].xsec_scaling_factor() *
143  incoming_particles_[1].xsec_scaling_factor();
144 }
double partial_cross_section_
Partial cross-section to the chosen outgoing channel.
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303

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

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.
NN → NR: Fit to HADES data, see Agakishiev:2014wqa.

Reimplemented from smash::Action.

Definition at line 249 of file scatteraction.cc.

250  {
253  // We potentially have more than two particles, so the following angular
254  // distributions don't work. Instead we just keep the angular
255  // distributions generated by string fragmentation.
256  return;
257  }
258  assert(outgoing_particles_.size() == 2);
259  const auto &log = logger<LogArea::ScatterAction>();
260 
261  // NN scattering is anisotropic currently
262  const bool nn_scattering = incoming_particles_[0].type().is_nucleon() &&
263  incoming_particles_[1].type().is_nucleon();
264  /* Elastic process is anisotropic and
265  * the angular distribution is based on the NN elastic scattering. */
266  const bool el_scattering = process_type_ == ProcessType::Elastic;
267 
268  const double mass_in_a = incoming_particles_[0].effective_mass();
269  const double mass_in_b = incoming_particles_[1].effective_mass();
270 
271  ParticleData *p_a = &outgoing_particles_[0];
272  ParticleData *p_b = &outgoing_particles_[1];
273 
274  const double mass_a = masses.first;
275  const double mass_b = masses.second;
276 
277  const std::array<double, 2> t_range = get_t_range<double>(
278  kinetic_energy_cm, mass_in_a, mass_in_b, mass_a, mass_b);
279  Angles phitheta;
280  if (el_scattering && !isotropic_) {
284  double mandelstam_s_new = 0.;
285  if (nn_scattering) {
286  mandelstam_s_new = mandelstam_s();
287  } else {
288  /* In the case of elastic collisions other than NN collisions,
289  * there is an ambiguity on how to get the lab-frame momentum (plab),
290  * since the incoming particles can have different masses.
291  * Right now, we first obtain the center-of-mass momentum
292  * of the collision (pcom_now).
293  * Then, the lab-frame momentum is evaluated from the mandelstam s,
294  * which yields the original center-of-mass momentum
295  * when nucleon mass is assumed. */
296  const double pcm_now = pCM_from_s(mandelstam_s(), mass_in_a, mass_in_b);
297  mandelstam_s_new =
298  4. * std::sqrt(pcm_now * pcm_now + nucleon_mass * nucleon_mass);
299  }
300  double bb, a, plab = plab_from_s(mandelstam_s_new);
301  if (nn_scattering &&
302  p_a->pdgcode().antiparticle_sign() ==
303  p_b->pdgcode().antiparticle_sign() &&
304  std::abs(p_a->type().charge() + p_b->type().charge()) == 1) {
305  // proton-neutron and antiproton-antineutron
306  bb = std::max(Cugnon_bnp(plab), really_small);
307  a = (plab < 0.8) ? 1. : 0.64 / (plab * plab);
308  } else {
309  /* all others including pp, nn and AQM elastic processes
310  * This is applied for all particle pairs, which are allowed to
311  * interact elastically. */
312  bb = std::max(Cugnon_bpp(plab), really_small);
313  a = 1.;
314  }
315  double t = random::expo(bb, t_range[0], t_range[1]);
316  if (random::canonical() > 1. / (1. + a)) {
317  t = t_range[0] + t_range[1] - t;
318  }
319  // determine scattering angles in center-of-mass frame
320  phitheta = Angles(2. * M_PI * random::canonical(),
321  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
322  } else if (nn_scattering && p_a->pdgcode().is_Delta() &&
323  p_b->pdgcode().is_nucleon() &&
324  p_a->pdgcode().antiparticle_sign() ==
325  p_b->pdgcode().antiparticle_sign() &&
326  !isotropic_) {
330  const double plab = plab_from_s(mandelstam_s());
331  const double bb = std::max(Cugnon_bpp(plab), really_small);
332  double t = random::expo(bb, t_range[0], t_range[1]);
333  if (random::canonical() > 0.5) {
334  t = t_range[0] + t_range[1] - t; // symmetrize
335  }
336  phitheta = Angles(2. * M_PI * random::canonical(),
337  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
338  } else if (nn_scattering && p_b->pdgcode().is_nucleon() && !isotropic_ &&
339  (p_a->type().is_Nstar() || p_a->type().is_Deltastar())) {
341  const std::array<double, 4> p{1.46434, 5.80311, -6.89358, 1.94302};
342  const double a = p[0] + mass_a * (p[1] + mass_a * (p[2] + mass_a * p[3]));
343  /* If the resonance is so heavy that the index "a" exceeds 30,
344  * the power function turns out to be too sharp. Take t directly to be
345  * t_0 in such a case. */
346  double t = t_range[0];
347  if (a < 30) {
348  t = random::power(-a, t_range[0], t_range[1]);
349  }
350  if (random::canonical() > 0.5) {
351  t = t_range[0] + t_range[1] - t; // symmetrize
352  }
353  phitheta = Angles(2. * M_PI * random::canonical(),
354  1. - 2. * (t - t_range[0]) / (t_range[1] - t_range[0]));
355  } else {
356  /* isotropic angular distribution */
357  phitheta.distribute_isotropically();
358  }
359 
360  ThreeVector pscatt = phitheta.threevec();
361  // 3-momentum of first incoming particle in center-of-mass frame
362  ThreeVector pcm =
363  incoming_particles_[0].momentum().LorentzBoost(beta_cm()).threevec();
364  pscatt.rotate_z_axis_to(pcm);
365 
366  // final-state CM momentum
367  const double p_f = pCM(kinetic_energy_cm, mass_a, mass_b);
368  if (!(p_f > 0.0)) {
369  log.warn("Particle: ", p_a->pdgcode(), " radial momentum: ", p_f);
370  log.warn("Etot: ", kinetic_energy_cm, " m_a: ", mass_a, " m_b: ", mass_b);
371  }
372  p_a->set_4momentum(mass_a, pscatt * p_f);
373  p_b->set_4momentum(mass_b, -pscatt * p_f);
374 
375  /* Debug message is printed before boost, so that p_a and p_b are
376  * the momenta in the center of mass frame and thus opposite to
377  * each other.*/
378  log.debug("p_a: ", *p_a, "\np_b: ", *p_b);
379 }
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:163
constexpr double really_small
Numerical error tolerance.
Definition: constants.h:34
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
ProcessType process_type_
type of process
Definition: action.h:320
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering particles in the calculation frame...
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:66
constexpr double nucleon_mass
Nucleon mass in GeV.
Definition: constants.h:52
T canonical()
Definition: random.h:110
static double Cugnon_bnp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic np...
double plab_from_s(double mandelstam_s, double mass)
Convert Mandelstam-s to p_lab in a fixed-target collision.
Definition: kinematics.h:157
double mandelstam_s() const
Determine the Mandelstam s variable,.
elastic scattering: particles remain the same, only momenta change
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
hard string process involving 2->2 QCD process by PYTHIA.
bool isotropic_
Do this collision isotropically?
static double Cugnon_bpp(double plab)
Computes the B coefficients from the Cugnon parametrization of the angular distribution in elastic pp...
T power(T n, T xMin, T xMax)
Draws a random number according to a power-law distribution ~ x^n.
Definition: random.h:200
constexpr int p
Proton.
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
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,
double  low_snn_cut,
bool  strings_switch,
bool  use_AQM,
bool  strings_with_probability,
NNbarTreatment  nnbar_treatment 
)

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]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

Definition at line 107 of file scatteraction.cc.

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

132  {
133  return collision_channels_;
134  }
CollisionBranchList collision_channels_
List of possible collisions.

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

152  {
153  string_process_ = str_proc;
154  }
StringProcess * string_process_
Pointer to interface class for strings.

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

161 { return total_cross_section_; }
double total_cross_section_
Total hadronic cross section.
Here is the caller graph for this function:

◆ mandelstam_s()

double smash::ScatterAction::mandelstam_s ( ) const
protected

Determine the Mandelstam s variable,.

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

Equal to the square of CMS energy.

Returns
Mandelstam s

Definition at line 154 of file scatteraction.cc.

154 { return total_momentum().sqr(); }
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:323
double sqr() const
calculate the square of the vector (which is a scalar)
Definition: fourvector.h:437
Here is the call graph for this function:
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 156 of file scatteraction.cc.

156  {
157  const double m1 = incoming_particles_[0].effective_mass();
158  const double m2 = incoming_particles_[1].effective_mass();
159  return pCM(sqrt_s(), m1, m2);
160 }
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:265
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 162 of file scatteraction.cc.

162  {
163  const double m1 = incoming_particles_[0].effective_mass();
164  const double m2 = incoming_particles_[1].effective_mass();
165  return pCM_sqr(sqrt_s(), m1, m2);
166 }
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:91
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:265
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 particles in the calculation frame.

Returns
boost velocity between center of mass and calculation frame.

Definition at line 146 of file scatteraction.cc.

146  {
147  return total_momentum().velocity();
148 }
ThreeVector velocity() const
Get the velocity (3-vector divided by zero component).
Definition: fourvector.h:310
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:323
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 150 of file scatteraction.cc.

150  {
151  return (1. / std::sqrt(1.0 - beta_cm().sqr()));
152 }
ThreeVector beta_cm() const
Get the velocity of the center of mass of the scattering particles in the calculation frame...
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 381 of file scatteraction.cc.

381  {
382  // copy initial particles into final state
385  // resample momenta
386  sample_angles({outgoing_particles_[0].effective_mass(),
387  outgoing_particles_[1].effective_mass()},
388  sqrt_s());
389 }
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
void sample_angles(std::pair< double, double > masses, double kinetic_energy_cm) override
Sample final-state angles in a 2->2 collision (possibly anisotropic).
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:265
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 391 of file scatteraction.cc.

391  {
392  // create new particles
394  /* Set the formation time of the 2 particles to the larger formation time of
395  * the incoming particles, if it is larger than the execution time;
396  * execution time is otherwise taken to be the formation time */
397  const double t0 = incoming_particles_[0].formation_time();
398  const double t1 = incoming_particles_[1].formation_time();
399 
400  const size_t index_tmax = (t0 > t1) ? 0 : 1;
401  const double form_time_begin =
402  incoming_particles_[index_tmax].begin_formation_time();
403  const double sc =
404  incoming_particles_[index_tmax].initial_xsec_scaling_factor();
405  if (t0 > time_of_execution_ || t1 > time_of_execution_) {
406  /* The newly produced particles are supposed to continue forming exactly
407  * like the latest forming ingoing particle. Therefore the details on the
408  * formation are adopted. The initial cross section scaling factor of the
409  * incoming particles is considered to also be the scaling factor of the
410  * newly produced outgoing particles.
411  */
412  outgoing_particles_[0].set_slow_formation_times(form_time_begin,
413  std::max(t0, t1));
414  outgoing_particles_[1].set_slow_formation_times(form_time_begin,
415  std::max(t0, t1));
416  outgoing_particles_[0].set_cross_section_scaling_factor(sc);
417  outgoing_particles_[1].set_cross_section_scaling_factor(sc);
418  } else {
419  outgoing_particles_[0].set_formation_time(time_of_execution_);
420  outgoing_particles_[1].set_formation_time(time_of_execution_);
421  }
422 }
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
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:208
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:317
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 465 of file scatteraction.cc.

465  {
466  assert(incoming_particles_.size() == 2);
467  const auto &log = logger<LogArea::Pythia>();
468  // Disable floating point exception trap for Pythia
469  {
470  DisableFloatTraps guard;
471  /* initialize the string_process_ object for this particular collision */
473  /* implement collision */
474  bool success = false;
475  int ntry = 0;
476  const int ntry_max = 10000;
477  while (!success && ntry < ntry_max) {
478  ntry++;
479  switch (process_type_) {
481  /* single diffractive to A+X */
482  success = string_process_->next_SDiff(true);
483  break;
485  /* single diffractive to X+B */
486  success = string_process_->next_SDiff(false);
487  break;
489  /* double diffractive */
490  success = string_process_->next_DDiff();
491  break;
493  /* soft non-diffractive */
494  success = string_process_->next_NDiffSoft();
495  break;
497  /* soft BBbar 2 mesonic annihilation */
498  success = string_process_->next_BBbarAnn();
499  break;
501  success = string_process_->next_NDiffHard();
502  break;
503  default:
504  log.error("Unknown string process required.");
505  success = false;
506  }
507  }
508  if (ntry == ntry_max) {
509  /* If pythia fails to form a string, it is usually because the energy
510  * is not large enough. In this case, an elastic scattering happens.
511  *
512  * Since particles are normally added after process selection for
513  * strings, outgoing_particles is still uninitialized, and memory
514  * needs to be allocated. We also shift the process_type_ to elastic
515  * so that sample_angles does a proper treatment. */
516  outgoing_particles_.reserve(2);
517  outgoing_particles_.push_back(ParticleData{incoming_particles_[0]});
518  outgoing_particles_.push_back(ParticleData{incoming_particles_[1]});
521  } else {
523  /* If the incoming particles already were unformed, the formation
524  * times and cross section scaling factors need to be adjusted */
525  const double tform_in = std::max(incoming_particles_[0].formation_time(),
526  incoming_particles_[1].formation_time());
527  if (tform_in > time_of_execution_) {
528  const double fin =
529  (incoming_particles_[0].formation_time() >
530  incoming_particles_[1].formation_time())
531  ? incoming_particles_[0].initial_xsec_scaling_factor()
532  : incoming_particles_[1].initial_xsec_scaling_factor();
533  for (size_t i = 0; i < outgoing_particles_.size(); i++) {
534  const double tform_out = outgoing_particles_[i].formation_time();
535  const double fout =
536  outgoing_particles_[i].initial_xsec_scaling_factor();
537  /* The new cross section scaling factor will be the product of the
538  * cross section scaling factor of the ingoing particles and of the
539  * outgoing ones (since the outgoing ones are also string fragments
540  * and thus take time to form). */
541  outgoing_particles_[i].set_cross_section_scaling_factor(fin * fout);
542  /* If the unformed incoming particles' formation time is larger than
543  * the current outgoing particle's formation time, then the latter
544  * is overwritten by the former*/
545  if (tform_in > tform_out) {
546  outgoing_particles_[i].set_slow_formation_times(time_of_execution_,
547  tform_in);
548  }
549  }
550  }
551  /* Check momentum difference for debugging */
552  FourVector out_mom;
553  for (ParticleData data : outgoing_particles_) {
554  out_mom += data.momentum();
555  }
556  log.debug("Incoming momenta string:", total_momentum());
557  log.debug("Outgoing momenta string:", out_mom);
558  }
559  }
560 }
bool next_NDiffSoft()
Soft Non-diffractive process is modelled in accordance with dual-topological approach Capella:1978ig...
double diffractive. Two strings are formed, one from A and one from B.
ProcessType process_type_
type of process
Definition: action.h:320
a special case of baryon-antibaryon annihilation.
StringProcess * string_process_
Pointer to interface class for strings.
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:323
elastic scattering: particles remain the same, only momenta change
void init(const ParticleList &incoming, double tcoll)
initialization feed intial particles, time of collision and gamma factor of the center of mass...
bool next_SDiff(bool is_AB_to_AX)
Single-diffractive process is based on single pomeron exchange described in Ingelman:1984ns.
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
bool next_DDiff()
Double-diffractive process ( A + B -> X + X ) is similar to the single-diffractive process...
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
hard string process involving 2->2 QCD process by PYTHIA.
void elastic_scattering()
Perform an elastic two-body scattering, i.e. just exchange momentum.
ParticleList get_final_state()
a function to get the final state particle list which is called after the collision ...
bool next_BBbarAnn()
Baryon-antibaryon annihilation process Based on what UrQMD Bass:1998ca, Bleicher:1999xi does...
non-diffractive. Two strings are formed both have ends in A and B.
bool next_NDiffHard()
Hard Non-diffractive process is based on PYTHIA 8 with partonic showers and interactions.
(41-45) soft string excitations.
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:317
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 424 of file scatteraction.cc.

424  {
425  const auto &log = logger<LogArea::ScatterAction>();
426 
427  if (outgoing_particles_.size() != 1) {
428  std::string s =
429  "resonance_formation: "
430  "Incorrect number of particles in final state: ";
431  s += std::to_string(outgoing_particles_.size()) + " (";
432  s += incoming_particles_[0].pdgcode().string() + " + ";
433  s += incoming_particles_[1].pdgcode().string() + ")";
434  throw InvalidResonanceFormation(s);
435  }
436  // Set the momentum of the formed resonance in its rest frame.
437  outgoing_particles_[0].set_4momentum(
438  total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
439  /* Set the formation time of the resonance to the larger formation time of the
440  * incoming particles, if it is larger than the execution time; execution time
441  * is otherwise taken to be the formation time */
442  const double t0 = incoming_particles_[0].formation_time();
443  const double t1 = incoming_particles_[1].formation_time();
444 
445  const size_t index_tmax = (t0 > t1) ? 0 : 1;
446  const double begin_form_time =
447  incoming_particles_[index_tmax].begin_formation_time();
448  const double sc =
449  incoming_particles_[index_tmax].initial_xsec_scaling_factor();
450  if (t0 > time_of_execution_ || t1 > time_of_execution_) {
451  outgoing_particles_[0].set_slow_formation_times(begin_form_time,
452  std::max(t0, t1));
453  outgoing_particles_[0].set_cross_section_scaling_factor(sc);
454  } else {
455  outgoing_particles_[0].set_formation_time(time_of_execution_);
456  }
457  /* this momentum is evaluated in the computational frame. */
458  log.debug("Momentum of the new particle: ",
459  outgoing_particles_[0].momentum());
460 }
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:311
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:303
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:123
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:317
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 221 of file scatteraction.h.

◆ total_cross_section_

double smash::ScatterAction::total_cross_section_
protected

Total hadronic cross section.

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

◆ isotropic_

bool smash::ScatterAction::isotropic_ = false
protected

Do this collision isotropically?

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

◆ string_process_

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

Pointer to interface class for strings.

Definition at line 250 of file scatteraction.h.


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