Version: SMASH-3.3
smash::ScatterActionMulti Class Reference

#include <scatteractionmulti.h>

ScatterActionMulti is a special action which takes any number of incoming particles and performs a scattering with the use of the stochastic criterion, producing one or more final-state particles.

Definition at line 23 of file scatteractionmulti.h.

Inheritance diagram for smash::ScatterActionMulti:
smash::Action

Classes

class  InvalidScatterActionMulti
 Thrown when ScatterActionMulti is called to perform with unknown combination of incoming and outgoing number of particles or unknown process type. More...
 

Public Member Functions

 ScatterActionMulti (const ParticleList &in_plist, double time, const SpinInteractionType spin_interaction_type=SpinInteractionType::Off)
 Construct a ScatterActionMulti object. More...
 
void generate_final_state () override
 Generate the final-state of the multi-particle scattering process. More...
 
double get_total_weight () const override
 Get the total probability for the reaction (scaled with the cross section scaling factors of the incoming particles). More...
 
double get_partial_weight () const override
 Get the partial probability for the chosen channel (scaled with the cross section scaling factors of the incoming particles). More...
 
void add_possible_reactions (double dt, const double gcell_vol, const MultiParticleReactionsBitSet incl_multi)
 Add all possible multi-particle reactions for the given incoming particles. More...
 
const CollisionBranchList & reaction_channels ()
 Get list of possible reaction channels. More...
 
double calculate_I3 (const double sqrts) const
 Calculate the integration necessary for the three-body phase space. More...
 
double parametrizaton_phi4 (const double man_s) const
 Calculate the parametrized 4-body relativistic phase space integral. 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, 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...
 
void assign_unpolarized_spin_vector_to_outgoing_particles ()
 Assign an unpolarized spin vector to all outgoing particles. More...
 

Protected Member Functions

void format_debug_output (std::ostream &out) const override
 Writes information about this 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_angles (std::pair< double, double > masses, double kinetic_energy_cm)
 Sample final-state momenta in general X->2 processes (here: using an isotropical angular distribution). 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...
 

Private Member Functions

void add_reaction (CollisionBranchPtr p)
 Add a new reaction channel. More...
 
void add_reactions (CollisionBranchList pv)
 Add several new reaction channels at once. More...
 
void annihilation ()
 Perform a n->1 annihilation process. More...
 
void n_to_two ()
 Perform a n->2 process. More...
 
double probability_three_to_one (const ParticleType &type_out, double dt, const double gcell_vol, const int degen_sym_factor=1) const
 Calculate the probability for a 3m-to-1 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61]. More...
 
double probability_three_to_two (const ParticleType &type_out1, const ParticleType &type_out2, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
 Calculate the probability for a 3-to-2 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61]. More...
 
double probability_four_to_two (const ParticleType &type_out1, const ParticleType &type_out2, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
 Calculate the probability for a 4-to-2 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61]. More...
 
double probability_five_to_two (const double m_out, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
 Calculate the probability for a 5-to-2 reaction according to the stochastic collision criterion as given in Garcia-Montero:2021haa [24]. More...
 
double parametrizaton_phi5_pions (const double man_s) const
 Calculate the parametrized 5-pion phase space. More...
 
double react_degen_factor (const int spin_factor_inc, const int spin_degen_out1, const int spin_degen_out2) const
 Determine the spin degeneracy factor ( \(D_{spin}\)) for the N->2 reaction. More...
 
bool three_different_pions (const ParticleData &data_a, const ParticleData &data_b, const ParticleData &data_c) const
 Check wether the three incoming particles are π⁺,π⁻,π⁰ in any order. More...
 
bool two_pions_eta (const ParticleData &data_a, const ParticleData &data_b, const ParticleData &data_c) const
 Check wether the three incoming particles are π⁺,π⁻,η or π⁰,π⁰,η in any order. More...
 
bool all_incoming_particles_are_pions_have_zero_charge_only_one_piz () const
 Check if 5 incoming particles match intial pion state for 5-to-2, which is pi+ pi- pi+ pi- pi0 in order to match the NNbar resonance treatment. More...
 

Private Attributes

double total_probability_
 Total probability of reaction. More...
 
double partial_probability_
 Partial probability of the chosen outgoing channel. More...
 
SpinInteractionType spin_interaction_type_
 Spin interaction type. More...
 
CollisionBranchList reaction_channels_
 List of possible collisions. 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...
 
- 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...
 

Constructor & Destructor Documentation

◆ ScatterActionMulti()

smash::ScatterActionMulti::ScatterActionMulti ( const ParticleList &  in_plist,
double  time,
const SpinInteractionType  spin_interaction_type = SpinInteractionType::Off 
)

Construct a ScatterActionMulti object.

Parameters
[in]in_plistList of reaction partners
[in]timeTime at which the action is supposed to take place
[in]spin_interaction_typeType of spin interaction

Definition at line 25 of file scatteractionmulti.cc.

28  : Action(in_plist, time),
30  spin_interaction_type_(spin_interaction_type) {}
Action(const ParticleList &in_part, double time)
Construct an action object with incoming particles and relative time.
Definition: action.h:44
double total_probability_
Total probability of reaction.
SpinInteractionType spin_interaction_type_
Spin interaction type.

Member Function Documentation

◆ generate_final_state()

void smash::ScatterActionMulti::generate_final_state ( )
overridevirtual

Generate the final-state of the multi-particle scattering process.

Assign position and momenta to outgoing particles.

Exceptions
InvalidScatterActionMulti

Implements smash::Action.

Definition at line 372 of file scatteractionmulti.cc.

372  {
373  logg[LScatterActionMulti].debug("Incoming particles for ScatterActionMulti: ",
375 
376  /* Decide for a particular final state. */
377  const CollisionBranch* proc =
378  choose_channel<CollisionBranch>(reaction_channels_, total_probability_);
379  process_type_ = proc->get_type();
380  outgoing_particles_ = proc->particle_list();
381  partial_probability_ = proc->weight();
382 
383  logg[LScatterActionMulti].debug("Chosen channel: ", process_type_,
385 
386  switch (process_type_) {
388  /* n->1 annihilation */
389  annihilation();
390  break;
394  /* n->2 scattering */
395  n_to_two();
396  break;
397  default:
398  throw InvalidScatterActionMulti(
399  "ScatterActionMulti::generate_final_state: Invalid process type " +
400  std::to_string(static_cast<int>(process_type_)) + " was requested.");
401  }
402 
403  /* The production point of the new particles. */
404  FourVector middle_point = get_interaction_point();
405 
406  for (ParticleData& new_particle : outgoing_particles_) {
407  // Boost to the computational frame
408  new_particle.boost_momentum(
410  /* Set positions of the outgoing particles */
411  new_particle.set_4position(middle_point);
412  }
413 }
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:160
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:372
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:364
FourVector get_interaction_point() const
Get the interaction point.
Definition: action.cc:68
ProcessType process_type_
type of process
Definition: action.h:381
CollisionBranchList reaction_channels_
List of possible collisions.
double partial_probability_
Partial probability of the chosen outgoing channel.
void n_to_two()
Perform a n->2 process.
void annihilation()
Perform a n->1 annihilation process.
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:40
static constexpr int LScatterActionMulti
@ MultiParticleThreeToTwo
See here for a short description.
@ MultiParticleThreeMesonsToOne
See here for a short description.
@ MultiParticleFourToTwo
See here for a short description.
@ MultiParticleFiveToTwo
See here for a short description.
std::string to_string(ThermodynamicQuantity quantity)
Convert a ThermodynamicQuantity enum value to its corresponding string.
Definition: stringify.cc:26

◆ get_total_weight()

double smash::ScatterActionMulti::get_total_weight ( ) const
overridevirtual

Get the total probability for the reaction (scaled with the cross section scaling factors of the incoming particles).

Returns
total probability.

Implements smash::Action.

Definition at line 41 of file scatteractionmulti.cc.

41  {
42  double xsec_scaling = 1.0;
43  for (const ParticleData& in_part : incoming_particles_) {
44  xsec_scaling *= in_part.xsec_scaling_factor();
45  }
46  return total_probability_ * xsec_scaling;
47 }

◆ get_partial_weight()

double smash::ScatterActionMulti::get_partial_weight ( ) const
overridevirtual

Get the partial probability for the chosen channel (scaled with the cross section scaling factors of the incoming particles).

Returns
partial probability.

Implements smash::Action.

Definition at line 49 of file scatteractionmulti.cc.

49  {
50  double xsec_scaling = 1.0;
51  for (const ParticleData& in_part : incoming_particles_) {
52  xsec_scaling *= in_part.xsec_scaling_factor();
53  }
54  return partial_probability_ * xsec_scaling;
55 }

◆ add_possible_reactions()

void smash::ScatterActionMulti::add_possible_reactions ( double  dt,
const double  gcell_vol,
const MultiParticleReactionsBitSet  incl_multi 
)

Add all possible multi-particle reactions for the given incoming particles.

Parameters
[in]dttimestep size
[in]gcell_volgcell_vol grid cell volume
[in]incl_multiWhich multi-particle reactions are enabled?

Definition at line 57 of file scatteractionmulti.cc.

59  {
60  // 3 -> m
61  if (incoming_particles_.size() == 3) {
62  // 3 -> 1
63  if (incl_multi[IncludedMultiParticleReactions::Meson_3to1] == 1) {
65  incoming_particles_[2])) {
66  // 3pi -> omega
67  const ParticleTypePtr type_omega = ParticleType::try_find(0x223);
68  if (type_omega) {
69  add_reaction(std::make_unique<CollisionBranch>(
70  *type_omega,
71  probability_three_to_one(*type_omega, dt, gcell_vol,
72  type_omega->spin_degeneracy()),
74  }
75  // 3pi -> phi
76  const ParticleTypePtr type_phi = ParticleType::try_find(0x333);
77  if (type_phi) {
78  add_reaction(std::make_unique<CollisionBranch>(
79  *type_phi,
80  probability_three_to_one(*type_phi, dt, gcell_vol,
81  type_phi->spin_degeneracy()),
83  }
85  incoming_particles_[2])) {
86  // eta2pi -> eta-prime
87  const ParticleTypePtr type_eta_prime = ParticleType::try_find(0x331);
88 
89  int sym_factor_in = 1;
90  if (incoming_particles_[0].type() == incoming_particles_[1].type() ||
91  incoming_particles_[1].type() == incoming_particles_[2].type() ||
92  incoming_particles_[2].type() == incoming_particles_[0].type()) {
93  sym_factor_in = 2; // 2 factorial
94  }
95 
96  if (type_eta_prime) {
97  add_reaction(std::make_unique<CollisionBranch>(
98  *type_eta_prime,
100  *type_eta_prime, dt, gcell_vol,
101  sym_factor_in * type_eta_prime->spin_degeneracy()),
103  }
104  }
105  }
106  // 3 -> 2
107  if (incl_multi[IncludedMultiParticleReactions::Deuteron_3to2] == 1) {
108  const PdgCode pdg_a = incoming_particles_[0].pdgcode();
109  const PdgCode pdg_b = incoming_particles_[1].pdgcode();
110  const PdgCode pdg_c = incoming_particles_[2].pdgcode();
111  const ParticleTypePtr type_deuteron =
113  const ParticleTypePtr type_anti_deuteron =
115 
116  const int spin_factor_inc = pdg_a.spin_degeneracy() *
117  pdg_b.spin_degeneracy() *
118  pdg_c.spin_degeneracy();
119 
120  if (type_deuteron && type_anti_deuteron) {
121  // πpn → πd
122  if ((pdg_a.is_pion() && pdg_b == pdg::p && pdg_c == pdg::n) ||
123  (pdg_a.is_pion() && pdg_b == pdg::n && pdg_c == pdg::p) ||
124  (pdg_a == pdg::p && pdg_b.is_pion() && pdg_c == pdg::n) ||
125  (pdg_a == pdg::n && pdg_b.is_pion() && pdg_c == pdg::p) ||
126  (pdg_a == pdg::p && pdg_b == pdg::n && pdg_c.is_pion()) ||
127  (pdg_a == pdg::n && pdg_b == pdg::p && pdg_c.is_pion())) {
128  // Get type of incoming π
129  ParticleList::iterator it = std::find_if(
131  [](ParticleData x) { return x.is_pion(); });
132  const ParticleType& type_pi = it->type();
133 
134  const double spin_degn =
135  react_degen_factor(spin_factor_inc, type_pi.spin_degeneracy(),
136  type_deuteron->spin_degeneracy());
137 
138  add_reaction(std::make_unique<CollisionBranch>(
139  type_pi, *type_deuteron,
140  probability_three_to_two(type_pi, *type_deuteron, dt, gcell_vol,
141  spin_degn),
143  }
144 
145  // πp̅n̅ → πd̅
146  if ((pdg_a.is_pion() && pdg_b == -pdg::p && pdg_c == -pdg::n) ||
147  (pdg_a.is_pion() && pdg_b == -pdg::n && pdg_c == -pdg::p) ||
148  (pdg_a == -pdg::p && pdg_b.is_pion() && pdg_c == -pdg::n) ||
149  (pdg_a == -pdg::n && pdg_b.is_pion() && pdg_c == -pdg::p) ||
150  (pdg_a == -pdg::p && pdg_b == -pdg::n && pdg_c.is_pion()) ||
151  (pdg_a == -pdg::n && pdg_b == -pdg::p && pdg_c.is_pion())) {
152  // Get type of incoming π
153  ParticleList::iterator it = std::find_if(
155  [](ParticleData x) { return x.is_pion(); });
156  const ParticleType& type_pi = it->type();
157 
158  const double spin_degn =
159  react_degen_factor(spin_factor_inc, type_pi.spin_degeneracy(),
160  type_anti_deuteron->spin_degeneracy());
161 
162  add_reaction(std::make_unique<CollisionBranch>(
163  type_pi, *type_anti_deuteron,
164  probability_three_to_two(type_pi, *type_anti_deuteron, dt,
165  gcell_vol, spin_degn),
167  }
168 
169  // Nnp → Nd, N̅np → N̅d
170  if ((pdg_a.is_nucleon() && pdg_b == pdg::p && pdg_c == pdg::n) ||
171  (pdg_a.is_nucleon() && pdg_b == pdg::n && pdg_c == pdg::p) ||
172  (pdg_a == pdg::p && pdg_b.is_nucleon() && pdg_c == pdg::n) ||
173  (pdg_a == pdg::n && pdg_b.is_nucleon() && pdg_c == pdg::p) ||
174  (pdg_a == pdg::p && pdg_b == pdg::n && pdg_c.is_nucleon()) ||
175  (pdg_a == pdg::n && pdg_b == pdg::p && pdg_c.is_nucleon())) {
176  int symmetry_factor = 1; // already true for N̅np → N̅d case
177 
178  ParticleList::iterator it =
179  std::find_if(incoming_particles_.begin(),
180  incoming_particles_.end(), [](ParticleData x) {
181  return x.pdgcode().antiparticle_sign() == -1;
182  });
183  if (it == incoming_particles_.end()) {
184  /* Meaning no anti-N found by find_if,
185  * therefore not N̅np → N̅d, but Nnp → Nd. */
186  symmetry_factor = 2; // for Nnp → Nd (2 factorial)
187  // It is already clear here that we have a double of two N
188  if (pdg_a == pdg_b) {
189  it = incoming_particles_.begin();
190  } else {
191  // If a and b are not the double, then c has to be part of it
192  it = incoming_particles_.begin() + 2;
193  }
194  }
195  const ParticleType& type_N = it->type();
196 
197  const double spin_degn =
198  react_degen_factor(spin_factor_inc, type_N.spin_degeneracy(),
199  type_deuteron->spin_degeneracy());
200 
201  add_reaction(std::make_unique<CollisionBranch>(
202  type_N, *type_deuteron,
203  probability_three_to_two(type_N, *type_deuteron, dt, gcell_vol,
204  symmetry_factor * spin_degn),
206  }
207 
208  // Np̅n̅ → Nd̅, N̅p̅n̅ → N̅d̅
209  if ((pdg_a.is_nucleon() && pdg_b == -pdg::p && pdg_c == -pdg::n) ||
210  (pdg_a.is_nucleon() && pdg_b == -pdg::n && pdg_c == -pdg::p) ||
211  (pdg_a == -pdg::p && pdg_b.is_nucleon() && pdg_c == -pdg::n) ||
212  (pdg_a == -pdg::n && pdg_b.is_nucleon() && pdg_c == -pdg::p) ||
213  (pdg_a == -pdg::p && pdg_b == -pdg::n && pdg_c.is_nucleon()) ||
214  (pdg_a == -pdg::n && pdg_b == -pdg::p && pdg_c.is_nucleon())) {
215  int symmetry_factor = 1; // already true for Np̅n̅ → Nd̅ case
216 
217  ParticleList::iterator it =
218  std::find_if(incoming_particles_.begin(),
219  incoming_particles_.end(), [](ParticleData x) {
220  return x.pdgcode().antiparticle_sign() == 1;
221  });
222  if (it == incoming_particles_.end()) {
223  /* Meaning no N found by find_if,
224  * therefore not Np̅n̅ → Nd̅, but N̅p̅n̅ → N̅d̅. */
225  symmetry_factor = 2; // for N̅p̅n̅ → N̅d̅ (2 factorial)
226  // It is already clear here that we have a double of two N̅
227  if (pdg_a == pdg_b) {
228  it = incoming_particles_.begin();
229  } else {
230  // If a and b are not the double, then c has to be part of it
231  it = incoming_particles_.begin() + 2;
232  }
233  }
234  const ParticleType& type_N = it->type();
235 
236  const double spin_degn =
237  react_degen_factor(spin_factor_inc, type_N.spin_degeneracy(),
238  type_anti_deuteron->spin_degeneracy());
239 
240  add_reaction(std::make_unique<CollisionBranch>(
241  type_N, *type_anti_deuteron,
242  probability_three_to_two(type_N, *type_anti_deuteron, dt,
243  gcell_vol, symmetry_factor * spin_degn),
245  }
246  }
247  }
248  }
249  // 4 -> 2
250  if (incoming_particles_.size() == 4 &&
252  std::map<PdgCode, int> c; // counts incoming PdgCodes
253  int spin_factor_inc = 1;
254  for (const ParticleData& data : incoming_particles_) {
255  c[data.pdgcode()]++;
256  spin_factor_inc *= data.pdgcode().spin_degeneracy();
257  }
258  // Nucleons, antinucleons, and pions can catalyze
259  const int n_possible_catalysts_incoming =
260  c[pdg::n] + c[pdg::p] + c[-pdg::p] + c[-pdg::n] + c[pdg::pi_p] +
261  c[pdg::pi_z] + c[pdg::pi_m];
262 
263  for (PdgCode pdg_nucleus :
266  const ParticleTypePtr type_nucleus = ParticleType::try_find(pdg_nucleus);
267  // Nucleus can be formed if and only if:
268  // 1) Incoming particles contain enough components (like p, n, Lambda)
269  // 2) In (incoming particles - components) there is still a catalyst
270  // This is including the situation like nnpp. Can be that t(nnp) is formed
271  // and p is catalyst, can be that he-3(ppn) is formed and n is catalyst.
272  // Both reactions should be added.
273  const int n_nucleus_components_that_can_be_catalysts =
274  pdg_nucleus.nucleus_p() + pdg_nucleus.nucleus_ap() +
275  pdg_nucleus.nucleus_n() + pdg_nucleus.nucleus_an();
276  const bool incoming_contain_nucleus_components =
277  c[pdg::p] >= pdg_nucleus.nucleus_p() &&
278  c[-pdg::p] >= pdg_nucleus.nucleus_ap() &&
279  c[pdg::n] >= pdg_nucleus.nucleus_n() &&
280  c[-pdg::n] >= pdg_nucleus.nucleus_an() &&
281  c[pdg::Lambda] >= pdg_nucleus.nucleus_La() &&
282  c[-pdg::Lambda] >= pdg_nucleus.nucleus_aLa();
283  const bool can_form_nucleus =
284  type_nucleus && incoming_contain_nucleus_components &&
285  n_possible_catalysts_incoming -
286  n_nucleus_components_that_can_be_catalysts ==
287  1;
288 
289  if (!can_form_nucleus) {
290  continue;
291  }
292  // Find the catalyst
293  std::map<PdgCode, int> catalyst_count = c;
294  catalyst_count[pdg::p] -= pdg_nucleus.nucleus_p();
295  catalyst_count[-pdg::p] -= pdg_nucleus.nucleus_ap();
296  catalyst_count[pdg::n] -= pdg_nucleus.nucleus_n();
297  catalyst_count[-pdg::n] -= pdg_nucleus.nucleus_an();
298  catalyst_count[pdg::Lambda] -= pdg_nucleus.nucleus_La();
299  catalyst_count[-pdg::Lambda] -= pdg_nucleus.nucleus_aLa();
300  PdgCode pdg_catalyst = PdgCode::invalid();
301  for (const auto i : catalyst_count) {
302  if (i.second == 1) {
303  pdg_catalyst = i.first;
304  break;
305  }
306  }
307  if (pdg_catalyst == PdgCode::invalid()) {
308  logg[LScatterActionMulti].error("Something went wrong while forming",
309  pdg_nucleus, " from ",
311  }
312  const ParticleTypePtr type_catalyst =
313  ParticleType::try_find(pdg_catalyst);
314  const double spin_degn =
315  react_degen_factor(spin_factor_inc, type_catalyst->spin_degeneracy(),
316  type_nucleus->spin_degeneracy());
317  double symmetry_factor = 1.0;
318  for (const auto i : c) {
319  symmetry_factor *= (i.second == 3) ? 6.0 // 3!
320  : (i.second == 2) ? 2.0 // 2!
321  : 1.0;
322  if (i.second > 3 || i.second < 0) {
323  logg[LScatterActionMulti].error("4<->2 error, incoming particles ",
325  }
326  }
327 
328  add_reaction(std::make_unique<CollisionBranch>(
329  *type_catalyst, *type_nucleus,
330  probability_four_to_two(*type_catalyst, *type_nucleus, dt, gcell_vol,
331  symmetry_factor * spin_degn),
333  }
334  }
335  // 5 -> 2
336  if (incoming_particles_.size() == 5) {
337  if (incl_multi[IncludedMultiParticleReactions::NNbar_5to2] == 1) {
339  const int spin_factor_inc =
340  incoming_particles_[0].pdgcode().spin_degeneracy() *
341  incoming_particles_[1].pdgcode().spin_degeneracy() *
342  incoming_particles_[2].pdgcode().spin_degeneracy() *
343  incoming_particles_[3].pdgcode().spin_degeneracy() *
344  incoming_particles_[4].pdgcode().spin_degeneracy();
345  const double symmetry_factor = 4.0; // 2! * 2!
346 
347  const ParticleTypePtr type_p = ParticleType::try_find(pdg::p);
348  const ParticleTypePtr type_anti_p = ParticleType::try_find(-pdg::p);
349  const ParticleTypePtr type_n = ParticleType::try_find(pdg::n);
350  const ParticleTypePtr type_anti_n = ParticleType::try_find(-pdg::n);
351 
352  const double spin_degn =
353  react_degen_factor(spin_factor_inc, type_p->spin_degeneracy(),
354  type_anti_p->spin_degeneracy());
355 
356  if (type_p && type_n) {
357  const double prob = probability_five_to_two(
358  type_p->mass(), dt, gcell_vol,
359  symmetry_factor * spin_degn); // same for ppbar and nnbar
360  add_reaction(std::make_unique<CollisionBranch>(
361  *type_p, *type_anti_p, prob,
363  add_reaction(std::make_unique<CollisionBranch>(
364  *type_n, *type_anti_n, prob,
366  }
367  }
368  }
369  }
370 }
static const ParticleTypePtr try_find(PdgCode pdgcode)
Returns the ParticleTypePtr for the given pdgcode.
Definition: particletype.cc:89
static PdgCode invalid()
PdgCode 0x0 is guaranteed not to be valid by the PDG standard, but it passes all tests here,...
Definition: pdgcode.h:826
double probability_four_to_two(const ParticleType &type_out1, const ParticleType &type_out2, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
Calculate the probability for a 4-to-2 reaction according to the stochastic collision criterion as gi...
bool two_pions_eta(const ParticleData &data_a, const ParticleData &data_b, const ParticleData &data_c) const
Check wether the three incoming particles are π⁺,π⁻,η or π⁰,π⁰,η in any order.
double probability_three_to_one(const ParticleType &type_out, double dt, const double gcell_vol, const int degen_sym_factor=1) const
Calculate the probability for a 3m-to-1 reaction according to the stochastic collision criterion as g...
bool all_incoming_particles_are_pions_have_zero_charge_only_one_piz() const
Check if 5 incoming particles match intial pion state for 5-to-2, which is pi+ pi- pi+ pi- pi0 in ord...
double react_degen_factor(const int spin_factor_inc, const int spin_degen_out1, const int spin_degen_out2) const
Determine the spin degeneracy factor ( ) for the N->2 reaction.
double probability_five_to_two(const double m_out, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
Calculate the probability for a 5-to-2 reaction according to the stochastic collision criterion as gi...
double probability_three_to_two(const ParticleType &type_out1, const ParticleType &type_out2, double dt, const double gcell_vol, const double degen_sym_factor=1.0) const
Calculate the probability for a 3-to-2 reaction according to the stochastic collision criterion as gi...
bool three_different_pions(const ParticleData &data_a, const ParticleData &data_b, const ParticleData &data_c) const
Check wether the three incoming particles are π⁺,π⁻,π⁰ in any order.
void add_reaction(CollisionBranchPtr p)
Add a new reaction channel.
@ NNbar_5to2
@ A3_Nuclei_4to2
@ Deuteron_3to2
@ Meson_3to1
constexpr int pi_p
π⁺.
constexpr int64_t antideuteron
Anti-deuteron in decimal digits.
constexpr int64_t triton
Triton.
constexpr int64_t antihe3
Anti-He-3.
constexpr int p
Proton.
constexpr int64_t hypertriton
Hypertriton.
constexpr int64_t antihypertriton
Anti-Hypertriton.
constexpr int64_t antitriton
Anti-triton.
constexpr int pi_z
π⁰.
constexpr int n
Neutron.
constexpr int64_t deuteron
Deuteron.
constexpr int Lambda
Λ.
constexpr int pi_m
π⁻.
constexpr int64_t he3
He-3.

◆ reaction_channels()

const CollisionBranchList& smash::ScatterActionMulti::reaction_channels ( )
inline

Get list of possible reaction channels.

Returns
list of possible reaction channels.

Definition at line 75 of file scatteractionmulti.h.

75 { return reaction_channels_; }

◆ calculate_I3()

double smash::ScatterActionMulti::calculate_I3 ( const double  sqrts) const

Calculate the integration necessary for the three-body phase space.

The defintion for the integral is given by

\[I_3 = \int dm^2_{23}dm^2_{12} = \int^{(M-m_3)^2}_{(m_1+m_2)^2}[m^2_{23, max}- m^2_{23, min}]dm^2_{12}\]

see PDG book (chapter Kinematics) for defintions of variables. The numbered masses reference the incoming particles and \(M\) the mass of the outgoing particles in this case, since we are looking at the backreaction to the 1-to-3 decay.

Parameters
[in]sqrtscenter of mass energy of incoming particles (= mass of outgoing particle)
Returns
result of integral

Definition at line 415 of file scatteractionmulti.cc.

415  {
416  const double m1 = incoming_particles_[0].type().mass();
417  const double m2 = incoming_particles_[1].type().mass();
418  const double m3 = incoming_particles_[2].type().mass();
419 
420  if (sqrts < m1 + m2 + m3) {
421  return 0.0;
422  }
423  const double x1 = (m1 - m2) * (m1 - m2), x2 = (m1 + m2) * (m1 + m2),
424  x3 = (sqrts - m3) * (sqrts - m3),
425  x4 = (sqrts + m3) * (sqrts + m3);
426  const double qmm = x3 - x1, qmp = x3 - x2, qpm = x4 - x1, qpp = x4 - x2;
427  const double kappa = std::sqrt(qpm * qmp / (qpp * qmm));
428  const double tmp = std::sqrt(qmm * qpp);
429  const double c1 =
430  4.0 * m1 * m2 * std::sqrt(qmm / qpp) * (x4 - m3 * sqrts + m1 * m2);
431  const double c2 = 0.5 * (m1 * m1 + m2 * m2 + m3 * m3 + sqrts * sqrts) * tmp;
432  const double c3 = 8 * m1 * m2 / tmp *
433  ((m1 * m1 + m2 * m2) * (m3 * m3 + sqrts * sqrts) -
434  2 * m1 * m1 * m2 * m2 - 2 * m3 * m3 * sqrts * sqrts);
435  const double c4 =
436  -8 * m1 * m2 / tmp * smash::pow_int(sqrts * sqrts - m3 * m3, 2);
437  const double res =
438  c1 * gsl_sf_ellint_Kcomp(kappa, GSL_PREC_DOUBLE) +
439  c2 * gsl_sf_ellint_Ecomp(kappa, GSL_PREC_DOUBLE) +
440  c3 * gsl_sf_ellint_Pcomp(kappa, -qmp / qmm, GSL_PREC_DOUBLE) +
441  c4 * gsl_sf_ellint_Pcomp(kappa, -x1 * qmp / (x2 * qmm), GSL_PREC_DOUBLE);
442  return res;
443 }
constexpr T pow_int(const T base, unsigned const exponent)
Efficient template for calculating integer powers using squaring.
Definition: pow.h:23

◆ parametrizaton_phi4()

double smash::ScatterActionMulti::parametrizaton_phi4 ( const double  man_s) const

Calculate the parametrized 4-body relativistic phase space integral.

The 4-body phase space is a n = 4 case of general integral over

\[d\Phi_n = (2\pi)^4 \prod_{i=1}^n \frac{d^3p_i}{2E_i (2\pi)^3} \times \delta(E_{tot} - \sum E_i) \delta^{(3)}(p_{tot} - \sum p_i)\]

. This is a Lorentz-invariant quantity, so the result of the integration depends only on

\[ s = E_{tot}^2 - p_{tot}^2\]

and masses of the particles

\[ m_i^2 = E_i^2 - p_i^2 \]

. The dimension in general case is

\[\mathrm{GeV}^{2n-4}\]

. The

\[\hbar = c = 1 \]

convention is used here, so when this integral is used in the acceptance probability for collision, one has to restore

\[\hbar\]

to obtain correct dimensionless probability. More on phase space integrals can be found, for example, in CERN-68-15 report. For developers, I (oliiny) have compiled a document with properties and parametrizations of many-body phase space integrals here: github.com/smash-transport/smash-devel/files/7791360/n_body_relativistic_phase_space.pdf

Parameters
[in]man_smandelstam s of reaction
Returns
phase space integral value for 4 bodies [GeV^4]

Definition at line 512 of file scatteractionmulti.cc.

512  {
513  int n_nucleons = 0, n_pions = 0, n_lambdas = 0;
514  double sum_m = 0.0, prod_m = 1.0;
515  for (const ParticleData& data : incoming_particles_) {
516  const PdgCode pdg = data.type().pdgcode();
517  n_nucleons += pdg.is_nucleon(); // including anti-nucleons
518  n_pions += pdg.is_pion();
519  n_lambdas += pdg.is_Lambda(); // including anti-Lambda
520  sum_m += data.type().mass();
521  prod_m *= data.type().mass();
522  }
523  const double x = 1.0 - sum_m / std::sqrt(man_s);
524  const double x2 = x * x;
525  const double x3 = x2 * x;
526  double g = -1.0;
527 
528  if (n_nucleons == 3 && n_pions == 1) { // NNNpi
529  g = (1.0 + 0.862432 * x - 3.4853 * x2 + 1.70259 * x3) /
530  (1.0 + 0.387376 * x - 1.34128 * x2 + 0.154489 * x3);
531  } else if (n_nucleons == 4) { // NNNN
532  g = (1.0 - 1.72285 * x + 0.728331 * x2) /
533  (1.0 - 0.967146 * x - 0.0103633 * x2);
534  } else if (n_nucleons == 2 && n_lambdas == 1 && n_pions == 1) { // LaNNpi
535  g = (1.0 + 0.937064 * x - 3.56864 * x2 + 1.721 * x3) /
536  (1.0 + 0.365202 * x - 1.2854 * x2 + 0.138444 * x3);
537  } else if (n_nucleons == 3 && n_lambdas == 1) { // LaNNN
538  g = (1.0 + 0.882401 * x - 3.4074 * x2 + 1.62454 * x3) /
539  (1.0 + 1.61741 * x - 2.12543 * x2 - 0.0902067 * x3);
540  }
541 
542  if (g > 0.0) {
543  return (std::sqrt(prod_m) * sum_m * sum_m * std::pow(x, 3.5) * g) /
544  (840. * std::sqrt(2) * std::pow(M_PI, 4.0) * std::pow(1 - x, 4.0));
545  } else {
546  logg[LScatterActionMulti].error("parametrizaton_phi4: no parametrization ",
547  "available for ", incoming_particles_);
548  return 0.0;
549  }
550 }

◆ format_debug_output()

void smash::ScatterActionMulti::format_debug_output ( std::ostream &  out) const
overrideprotectedvirtual

Writes information about this action to the out stream.

Parameters
[out]outout stream to be written to

Implements smash::Action.

Definition at line 665 of file scatteractionmulti.cc.

665  {
666  out << "MultiParticleScatter of " << incoming_particles_;
667  if (outgoing_particles_.empty()) {
668  out << " (not performed)";
669  } else {
670  out << " to " << outgoing_particles_;
671  }
672 }

◆ add_reaction()

void smash::ScatterActionMulti::add_reaction ( CollisionBranchPtr  p)
private

Add a new reaction channel.

Parameters
[in]pChannel to be added.

Definition at line 32 of file scatteractionmulti.cc.

32  {
33  add_process<CollisionBranch>(p, reaction_channels_, total_probability_);
34 }

◆ add_reactions()

void smash::ScatterActionMulti::add_reactions ( CollisionBranchList  pv)
private

Add several new reaction channels at once.

Parameters
[in]pvlist of channels to be added.

Definition at line 36 of file scatteractionmulti.cc.

36  {
37  add_processes<CollisionBranch>(std::move(pv), reaction_channels_,
39 }

◆ annihilation()

void smash::ScatterActionMulti::annihilation ( )
private

Perform a n->1 annihilation process.

Exceptions
InvalidScatterActionMulti

Definition at line 583 of file scatteractionmulti.cc.

583  {
584  if (outgoing_particles_.size() != 1) {
585  std::string s =
586  "Annihilation: "
587  "Incorrect number of particles in final state: ";
588  s += std::to_string(outgoing_particles_.size()) + ".";
589  throw InvalidScatterActionMulti(s);
590  }
591  // Set the momentum of the formed particle in its rest frame.
592  outgoing_particles_[0].set_4momentum(
593  total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
594  // Make sure to assign formation times before boost to the computational frame
598  }
599 
600  logg[LScatterActionMulti].debug("Momentum of the new particle: ",
601  outgoing_particles_[0].momentum());
602 }
void assign_formation_time_to_outgoing_particles()
Assign the formation time to the outgoing particles.
Definition: action.cc:191
void assign_unpolarized_spin_vector_to_outgoing_particles()
Assign an unpolarized spin vector to all outgoing particles.
Definition: action.cc:478
@ Off
No spin interactions.

◆ n_to_two()

void smash::ScatterActionMulti::n_to_two ( )
private

Perform a n->2 process.

Definition at line 604 of file scatteractionmulti.cc.

604  {
606  // Make sure to assign formation times before boost to the computational frame
610  }
612  "->2 scattering:", incoming_particles_,
613  " -> ", outgoing_particles_);
614 }
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

◆ probability_three_to_one()

double smash::ScatterActionMulti::probability_three_to_one ( const ParticleType type_out,
double  dt,
const double  gcell_vol,
const int  degen_sym_factor = 1 
) const
private

Calculate the probability for a 3m-to-1 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61].

The formula for the probablilty is not taken from a reference, but derived following the same idea as specified e.g. in the paper above.

\[ P_{3\rightarrow 1} = \frac{\Delta t}{(\Delta^3x)^2} \frac{\pi}{4E_1E_2E_3}\frac{\Gamma_{1\rightarrow3}}{\Phi_3} \mathcal{A}(\sqrt{s}),\]

where \(\Phi_3\) represents the 3-body phase space:

\[\Phi_3 = \frac{1}{(2\pi)^3)}\frac{1}{16M^2}I_3.\]

The defintion for the integral \(I_3\) is given in the documentation of calculate_I3(). Degeneracy and symmetry factors are neglected in the formula, since they are treated as input for the function.

Parameters
[in]type_outtype of outgoing particle
[in]dttimestep size
[in]gcell_volgrid cell volume
[in]degen_sym_factordegeneracy factor for reaction (including symmetry factors)
Returns
probabilty for 3-to-1 reaction

Definition at line 445 of file scatteractionmulti.cc.

447  {
448  const double e1 = incoming_particles_[0].momentum().x0();
449  const double e2 = incoming_particles_[1].momentum().x0();
450  const double e3 = incoming_particles_[2].momentum().x0();
451  const double sqrts = sqrt_s();
452 
453  const double gamma_decay = type_out.get_partial_width(
454  sqrts, {&incoming_particles_[0].type(), &incoming_particles_[1].type(),
455  &incoming_particles_[2].type()});
456 
457  const double I_3 = calculate_I3(sqrts);
458  const double ph_sp_3 =
459  1. / (8 * M_PI * M_PI * M_PI) * 1. / (16 * sqrts * sqrts) * I_3;
460 
461  const double spec_f_val = type_out.spectral_function(sqrts);
462 
463  return dt / (gcell_vol * gcell_vol) * M_PI / (4. * e1 * e2 * e3) *
464  gamma_decay / ph_sp_3 * spec_f_val * std::pow(hbarc, 5.0) *
465  degen_sym_factor;
466 }
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:271
double calculate_I3(const double sqrts) const
Calculate the integration necessary for the three-body phase space.
constexpr double hbarc
GeV <-> fm conversion factor.
Definition: constants.h:29

◆ probability_three_to_two()

double smash::ScatterActionMulti::probability_three_to_two ( const ParticleType type_out1,
const ParticleType type_out2,
double  dt,
const double  gcell_vol,
const double  degen_sym_factor = 1.0 
) const
private

Calculate the probability for a 3-to-2 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61].

\[ P_{3 \rightarrow 2} = \frac{1}{4E_1E_2E_3} \frac{\Delta t}{(\Delta^3 x)^2} \frac{\tilde{\lambda}}{\Phi_38\pi s}\sigma_{2 \rightarrow 3},\]

where \(\Phi_3\) represents the 3-body phase space:

\[\Phi_3 = \frac{1}{(2\pi)^3)}\frac{1}{16M^2}I_3.\]

The defintion for the integral \(I_3\) is given in the documentation of calculate_I3(). Degeneracy and symmetry factors are neglected in the formula, since they are treated as input for the function.

Parameters
[in]type_out1type of outgoing particle 1
[in]type_out2type of outgoing particle 2
[in]dttimestep size
[in]gcell_volgrid cell volume
[in]degen_sym_factordegeneracy factor for reaction (including symmetry factors)
Returns
probabilty for 3-to-2 reaction

Definition at line 468 of file scatteractionmulti.cc.

470  {
471  const double e1 = incoming_particles_[0].momentum().x0();
472  const double e2 = incoming_particles_[1].momentum().x0();
473  const double e3 = incoming_particles_[2].momentum().x0();
474  const double m4 = type_out1.mass();
475  const double m5 = type_out2.mass();
476 
477  const double sqrts = sqrt_s();
478  const double xs =
479  CrossSections::two_to_three_xs(type_out1, type_out2, sqrts) / gev2_mb;
480  const double lamb = lambda_tilde(sqrts * sqrts, m4 * m4, m5 * m5);
481 
482  const double I_3 = calculate_I3(sqrts);
483  const double ph_sp_3 =
484  1. / (8 * M_PI * M_PI * M_PI) * 1. / (16 * sqrts * sqrts) * I_3;
485 
486  return dt / (gcell_vol * gcell_vol) * 1. / (4. * e1 * e2 * e3) * lamb /
487  (ph_sp_3 * 8 * M_PI * sqrts * sqrts) * xs * std::pow(hbarc, 5.0) *
488  degen_sym_factor;
489 }
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
static double two_to_three_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->3 cross section for the scattering of the given particle types.
constexpr double gev2_mb
GeV^-2 <-> mb conversion factor.
Definition: constants.h:35

◆ probability_four_to_two()

double smash::ScatterActionMulti::probability_four_to_two ( const ParticleType type_out1,
const ParticleType type_out2,
double  dt,
const double  gcell_vol,
const double  degen_sym_factor = 1.0 
) const
private

Calculate the probability for a 4-to-2 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [61].

\[ P_{4 \rightarrow 2} = \frac{1}{16E_1E_2E_3E_4} \frac{\Delta t}{(\Delta^3 x)^3} \frac{\tilde{\lambda}}{\Phi_44\pi s} \sigma_{2 \rightarrow 4},\]

where \(\Phi_4\) represents the 4-body phase space. Degeneracy and symmetry factors are neglected in the formula, since they are treated as input for the function.

Parameters
[in]type_out1type of outgoing particle 1
[in]type_out2type of outgoing particle 2
[in]dttimestep size
[in]gcell_volgrid cell volume
[in]degen_sym_factordegeneracy factor for reaction (including symmetry factors)
Returns
probabilty for 4-to-2 reaction

Definition at line 491 of file scatteractionmulti.cc.

493  {
494  const double e1 = incoming_particles_[0].momentum().x0();
495  const double e2 = incoming_particles_[1].momentum().x0();
496  const double e3 = incoming_particles_[2].momentum().x0();
497  const double e4 = incoming_particles_[3].momentum().x0();
498  const double m5 = type_out1.mass();
499  const double m6 = type_out2.mass();
500 
501  const double man_s = sqrt_s() * sqrt_s();
502  const double xs =
503  CrossSections::two_to_four_xs(type_out1, type_out2, sqrt_s()) / gev2_mb;
504  const double lamb = lambda_tilde(man_s, m5 * m5, m6 * m6);
505  const double ph_sp_4 = parametrizaton_phi4(man_s);
506 
507  return dt / std::pow(gcell_vol, 3.0) * 1. / (16. * e1 * e2 * e3 * e4) * xs /
508  (4. * M_PI * man_s) * lamb / ph_sp_4 * std::pow(hbarc, 8.0) *
509  degen_sym_factor;
510 }
static double two_to_four_xs(const ParticleType &type_in1, const ParticleType &type_in2, double sqrts)
Determine 2->4 cross section for the scattering of the given particle types.
double parametrizaton_phi4(const double man_s) const
Calculate the parametrized 4-body relativistic phase space integral.

◆ probability_five_to_two()

double smash::ScatterActionMulti::probability_five_to_two ( const double  m_out,
double  dt,
const double  gcell_vol,
const double  degen_sym_factor = 1.0 
) const
private

Calculate the probability for a 5-to-2 reaction according to the stochastic collision criterion as given in Garcia-Montero:2021haa [24].

\[ P_{5 \rightarrow 2} = \frac{1}{32E_1E_2E_3E_4E_5} \frac{\Delta t}{(\Delta^3 x)^4} \frac{\tilde{\lambda}}{\Phi_54\pi s}\sigma_{2 \rightarrow 5},\]

where the defintion is given without the necessary symmetry and spin degeneracy factors, which are input to the function and \(\Phi_5\) represents the 5-body phase space, which is paramaterized for the relevent 5 pion state here, see documentation of parametrizaton_phi5_pions().

Parameters
[in]m_outmass of outgoing particle types (assumes equal masses)
[in]dttimestep size
[in]gcell_volgrid cell volume
[in]degen_sym_factordegeneracy factor for reaction (including symmetry factors)
Returns
probabilty for 5-to-2 reaction

Definition at line 552 of file scatteractionmulti.cc.

554  {
555  const double e1 = incoming_particles_[0].momentum().x0();
556  const double e2 = incoming_particles_[1].momentum().x0();
557  const double e3 = incoming_particles_[2].momentum().x0();
558  const double e4 = incoming_particles_[3].momentum().x0();
559  const double e5 = incoming_particles_[4].momentum().x0();
560 
561  const double man_s = sqrt_s() * sqrt_s();
562  const double lamb = lambda_tilde(man_s, mout * mout, mout * mout);
563  const double ph_sp_5 = parametrizaton_phi5_pions(man_s);
564 
565  // Matching the NNbar anihilation cross section defintion for 2-to-5
566  const double xs =
567  std::max(0., ppbar_total(man_s) - ppbar_elastic(man_s)) / gev2_mb;
568 
569  return dt / std::pow(gcell_vol, 4.0) * 1. / (32. * e1 * e2 * e3 * e4 * e5) *
570  xs / (4. * M_PI * man_s) * lamb / ph_sp_5 * std::pow(hbarc, 11.0) *
571  degen_sym_factor;
572 }
double parametrizaton_phi5_pions(const double man_s) const
Calculate the parametrized 5-pion phase space.
double ppbar_total(double mandelstam_s)
ppbar total cross section parametrization Source: Bass:1998ca
double ppbar_elastic(double mandelstam_s)
ppbar elastic cross section parametrization Source: Bass:1998ca

◆ parametrizaton_phi5_pions()

double smash::ScatterActionMulti::parametrizaton_phi5_pions ( const double  man_s) const
private

Calculate the parametrized 5-pion phase space.

The defintion for the paramterization is given by

\[\Phi_5^{param.} = A(s-s_0)^5(1+\frac{s}{s_0})^{-\alpha}\]

with \(s_0 = 25 m_{\pi}^2\). \(A\) and \(\alpha\) are fitted to reproduce the phase space distribution.

Parameters
[in]man_smandelstam s of reaction
Returns
phase space value for 5 pions

Definition at line 574 of file scatteractionmulti.cc.

574  {
575  // see function documentation for parameter naming
576  const double s_zero = 25 * pion_mass * pion_mass;
577  const double fit_a = 2.1018e-10;
578  const double fit_alpha = 1.982;
579  return fit_a * std::pow(man_s - s_zero, 5.0) *
580  std::pow(1 + man_s / s_zero, -fit_alpha);
581 }
constexpr double pion_mass
Pion mass in GeV.
Definition: constants.h:69

◆ react_degen_factor()

double smash::ScatterActionMulti::react_degen_factor ( const int  spin_factor_inc,
const int  spin_degen_out1,
const int  spin_degen_out2 
) const
inlineprivate

Determine the spin degeneracy factor ( \(D_{spin}\)) for the N->2 reaction.

\[D_{spin} = \frac{(2J_{out1}+1)(2J_{out2}+1)} {(2J_{in1}+1)(2J_{in2}+1)(2J_{in3}+1)...(2J_{inN}+1)}\]

Parameters
[in]spin_factor_incproduct of incoming spin degeneracy (denominator in above expression)
[in]spin_degen_out1degeneracy factor of outgoing particle 1
[in]spin_degen_out2degeneracy factor of outgoing particle 2
Returns
spin degeneracy factor

Definition at line 286 of file scatteractionmulti.h.

288  {
289  return static_cast<double>(spin_degen_out1 * spin_degen_out2) /
290  static_cast<double>(spin_factor_inc);
291  }

◆ three_different_pions()

bool smash::ScatterActionMulti::three_different_pions ( const ParticleData data_a,
const ParticleData data_b,
const ParticleData data_c 
) const
private

Check wether the three incoming particles are π⁺,π⁻,π⁰ in any order.

Wrapper for unwieldy if statment.

Parameters
[in]data_adata for first incoming particle
[in]data_bdata for second incoming particle
[in]data_cdata for third incoming particle
Returns
true if combination of π⁺,π⁻,π⁰

Definition at line 616 of file scatteractionmulti.cc.

618  {
619  // We want a combination of pi+, pi- and pi0
620  const PdgCode pdg_a = data_a.pdgcode();
621  const PdgCode pdg_b = data_b.pdgcode();
622  const PdgCode pdg_c = data_c.pdgcode();
623 
624  return (pdg_a.is_pion() && pdg_b.is_pion() && pdg_c.is_pion()) &&
625  (pdg_a != pdg_b && pdg_b != pdg_c && pdg_c != pdg_a);
626 }

◆ two_pions_eta()

bool smash::ScatterActionMulti::two_pions_eta ( const ParticleData data_a,
const ParticleData data_b,
const ParticleData data_c 
) const
private

Check wether the three incoming particles are π⁺,π⁻,η or π⁰,π⁰,η in any order.

Wrapper for unwieldy if statment.

Parameters
[in]data_adata for first incoming particle
[in]data_bdata for second incoming particle
[in]data_cdata for third incoming particle
Returns
true if combination of π⁺,π⁻,η or π⁰,π⁰,η

Definition at line 628 of file scatteractionmulti.cc.

630  {
631  // We want a combination of pi0, pi0 and eta or pi+, pi- and eta
632  const PdgCode pdg_a = data_a.pdgcode();
633  const PdgCode pdg_b = data_b.pdgcode();
634  const PdgCode pdg_c = data_c.pdgcode();
635 
636  return (pdg_a == pdg::pi_z && pdg_b == pdg::pi_z && pdg_c == pdg::eta) ||
637  (pdg_a == pdg::pi_z && pdg_b == pdg::eta && pdg_c == pdg::pi_z) ||
638  (pdg_a == pdg::eta && pdg_b == pdg::pi_z && pdg_c == pdg::pi_z) ||
639 
640  (pdg_a == pdg::eta && pdg_b == pdg::pi_m && pdg_c == pdg::pi_p) ||
641  (pdg_a == pdg::eta && pdg_b == pdg::pi_p && pdg_c == pdg::pi_m) ||
642  (pdg_a == pdg::pi_m && pdg_b == pdg::pi_p && pdg_c == pdg::eta) ||
643  (pdg_a == pdg::pi_m && pdg_b == pdg::eta && pdg_c == pdg::pi_p) ||
644  (pdg_a == pdg::pi_p && pdg_b == pdg::pi_m && pdg_c == pdg::eta) ||
645  (pdg_a == pdg::pi_p && pdg_b == pdg::eta && pdg_c == pdg::pi_m);
646 }
constexpr int eta
η.

◆ all_incoming_particles_are_pions_have_zero_charge_only_one_piz()

bool smash::ScatterActionMulti::all_incoming_particles_are_pions_have_zero_charge_only_one_piz ( ) const
private

Check if 5 incoming particles match intial pion state for 5-to-2, which is pi+ pi- pi+ pi- pi0 in order to match the NNbar resonance treatment.

Definition at line 648 of file scatteractionmulti.cc.

649  {
650  const bool all_inc_pi =
652  [](const ParticleData& data) { return data.is_pion(); });
653  const int no_of_piz = std::count_if(
655  [](const ParticleData& data) { return data.pdgcode() == pdg::pi_z; });
656 
657  int total_state_charge = 0;
658  for (const ParticleData& part : incoming_particles_) {
659  total_state_charge += part.pdgcode().charge();
660  }
661 
662  return (all_inc_pi && total_state_charge == 0 && no_of_piz == 1);
663 }
bool all_of(Container &&c, UnaryPredicate &&p)
Convenience wrapper for std::all_of that operates on a complete container.
Definition: algorithms.h:80

Member Data Documentation

◆ total_probability_

double smash::ScatterActionMulti::total_probability_
private

Total probability of reaction.

Definition at line 325 of file scatteractionmulti.h.

◆ partial_probability_

double smash::ScatterActionMulti::partial_probability_
private

Partial probability of the chosen outgoing channel.

Definition at line 328 of file scatteractionmulti.h.

◆ spin_interaction_type_

SpinInteractionType smash::ScatterActionMulti::spin_interaction_type_
private

Spin interaction type.

Definition at line 331 of file scatteractionmulti.h.

◆ reaction_channels_

CollisionBranchList smash::ScatterActionMulti::reaction_channels_
private

List of possible collisions.

Definition at line 334 of file scatteractionmulti.h.


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