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:

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)
	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, FourVector >	get_potential_at_interaction_point () const
	Get the skyrme and asymmetry potential at the interaction point. More...
void	set_stochastic_pos_idx ()
	Setter function that stores a random incoming particle index latter used to determine the interaction point. More...

Protected Member Functions
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...
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 [59]. 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 [59]. 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 [59]. 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...
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
	)

Construct a ScatterActionMulti object.

Parameters

[in]	in_plist	List of reaction partners
[in]	time	Time at which the action is supposed to take place

Definition at line 25 of file scatteractionmulti.cc.

    : Action(in_plist, time), total_probability_(0.) {}

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 369 of file scatteractionmulti.cc.

                                              {
  logg[LScatterActionMulti].debug("Incoming particles for ScatterActionMulti: ",
                                  incoming_particles_);
 
  /* Decide for a particular final state. */
  const CollisionBranch* proc =
      choose_channel<CollisionBranch>(reaction_channels_, total_probability_);
  process_type_ = proc->get_type();
  outgoing_particles_ = proc->particle_list();
  partial_probability_ = proc->weight();
 
  logg[LScatterActionMulti].debug("Chosen channel: ", process_type_,
                                  outgoing_particles_);
 
  switch (process_type_) {
    case ProcessType::MultiParticleThreeMesonsToOne:
      /* n->1 annihilation */
      annihilation();
      break;
    case ProcessType::MultiParticleThreeToTwo:
    case ProcessType::MultiParticleFourToTwo:
    case ProcessType::MultiParticleFiveToTwo:
      /* n->2 scattering */
      n_to_two();
      break;
    default:
      throw InvalidScatterActionMulti(
          "ScatterActionMulti::generate_final_state: Invalid process type " +
          std::to_string(static_cast<int>(process_type_)) + " was requested.");
  }
 
  /* The production point of the new particles.  */
  FourVector middle_point = get_interaction_point();
 
  for (ParticleData& new_particle : outgoing_particles_) {
    // Boost to the computational frame
    new_particle.boost_momentum(
        -total_momentum_of_outgoing_particles().velocity());
    /* Set positions of the outgoing particles */
    new_particle.set_4position(middle_point);
  }
}

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 38 of file scatteractionmulti.cc.

                                                  {
  double xsec_scaling = 1.0;
  for (const ParticleData& in_part : incoming_particles_) {
    xsec_scaling *= in_part.xsec_scaling_factor();
  }
  return total_probability_ * xsec_scaling;
}

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 46 of file scatteractionmulti.cc.

                                                    {
  double xsec_scaling = 1.0;
  for (const ParticleData& in_part : incoming_particles_) {
    xsec_scaling *= in_part.xsec_scaling_factor();
  }
  return partial_probability_ * xsec_scaling;
}

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]	dt	timestep size
[in]	gcell_vol	gcell_vol grid cell volume
[in]	incl_multi	Which multi-particle reactions are enabled?

Definition at line 54 of file scatteractionmulti.cc.

                                                   {
  // 3 -> m
  if (incoming_particles_.size() == 3) {
    // 3 -> 1
    if (incl_multi[IncludedMultiParticleReactions::Meson_3to1] == 1) {
      if (three_different_pions(incoming_particles_[0], incoming_particles_[1],
                                incoming_particles_[2])) {
        // 3pi -> omega
        const ParticleTypePtr type_omega = ParticleType::try_find(0x223);
        if (type_omega) {
          add_reaction(std::make_unique<CollisionBranch>(
              *type_omega,
              probability_three_to_one(*type_omega, dt, gcell_vol,
                                       type_omega->spin_degeneracy()),
              ProcessType::MultiParticleThreeMesonsToOne));
        }
        // 3pi -> phi
        const ParticleTypePtr type_phi = ParticleType::try_find(0x333);
        if (type_phi) {
          add_reaction(std::make_unique<CollisionBranch>(
              *type_phi,
              probability_three_to_one(*type_phi, dt, gcell_vol,
                                       type_phi->spin_degeneracy()),
              ProcessType::MultiParticleThreeMesonsToOne));
        }
      } else if (two_pions_eta(incoming_particles_[0], incoming_particles_[1],
                               incoming_particles_[2])) {
        // eta2pi -> eta-prime
        const ParticleTypePtr type_eta_prime = ParticleType::try_find(0x331);
 
        int sym_factor_in = 1;
        if (incoming_particles_[0].type() == incoming_particles_[1].type() ||
            incoming_particles_[1].type() == incoming_particles_[2].type() ||
            incoming_particles_[2].type() == incoming_particles_[0].type()) {
          sym_factor_in = 2;  // 2 factorial
        }
 
        if (type_eta_prime) {
          add_reaction(std::make_unique<CollisionBranch>(
              *type_eta_prime,
              probability_three_to_one(
                  *type_eta_prime, dt, gcell_vol,
                  sym_factor_in * type_eta_prime->spin_degeneracy()),
              ProcessType::MultiParticleThreeMesonsToOne));
        }
      }
    }
    // 3 -> 2
    if (incl_multi[IncludedMultiParticleReactions::Deuteron_3to2] == 1) {
      const PdgCode pdg_a = incoming_particles_[0].pdgcode();
      const PdgCode pdg_b = incoming_particles_[1].pdgcode();
      const PdgCode pdg_c = incoming_particles_[2].pdgcode();
      const ParticleTypePtr type_deuteron =
          ParticleType::try_find(pdg::deuteron);
      const ParticleTypePtr type_anti_deuteron =
          ParticleType::try_find(pdg::antideuteron);
 
      const int spin_factor_inc = pdg_a.spin_degeneracy() *
                                  pdg_b.spin_degeneracy() *
                                  pdg_c.spin_degeneracy();
 
      if (type_deuteron && type_anti_deuteron) {
        // πpn → πd
        if ((pdg_a.is_pion() && pdg_b == pdg::p && pdg_c == pdg::n) ||
            (pdg_a.is_pion() && pdg_b == pdg::n && pdg_c == pdg::p) ||
            (pdg_a == pdg::p && pdg_b.is_pion() && pdg_c == pdg::n) ||
            (pdg_a == pdg::n && pdg_b.is_pion() && pdg_c == pdg::p) ||
            (pdg_a == pdg::p && pdg_b == pdg::n && pdg_c.is_pion()) ||
            (pdg_a == pdg::n && pdg_b == pdg::p && pdg_c.is_pion())) {
          // Get type of incoming π
          ParticleList::iterator it = std::find_if(
              incoming_particles_.begin(), incoming_particles_.end(),
              [](ParticleData x) { return x.is_pion(); });
          const ParticleType& type_pi = it->type();
 
          const double spin_degn =
              react_degen_factor(spin_factor_inc, type_pi.spin_degeneracy(),
                                 type_deuteron->spin_degeneracy());
 
          add_reaction(std::make_unique<CollisionBranch>(
              type_pi, *type_deuteron,
              probability_three_to_two(type_pi, *type_deuteron, dt, gcell_vol,
                                       spin_degn),
              ProcessType::MultiParticleThreeToTwo));
        }
 
        // πp̅n̅ → πd̅
        if ((pdg_a.is_pion() && pdg_b == -pdg::p && pdg_c == -pdg::n) ||
            (pdg_a.is_pion() && pdg_b == -pdg::n && pdg_c == -pdg::p) ||
            (pdg_a == -pdg::p && pdg_b.is_pion() && pdg_c == -pdg::n) ||
            (pdg_a == -pdg::n && pdg_b.is_pion() && pdg_c == -pdg::p) ||
            (pdg_a == -pdg::p && pdg_b == -pdg::n && pdg_c.is_pion()) ||
            (pdg_a == -pdg::n && pdg_b == -pdg::p && pdg_c.is_pion())) {
          // Get type of incoming π
          ParticleList::iterator it = std::find_if(
              incoming_particles_.begin(), incoming_particles_.end(),
              [](ParticleData x) { return x.is_pion(); });
          const ParticleType& type_pi = it->type();
 
          const double spin_degn =
              react_degen_factor(spin_factor_inc, type_pi.spin_degeneracy(),
                                 type_anti_deuteron->spin_degeneracy());
 
          add_reaction(std::make_unique<CollisionBranch>(
              type_pi, *type_anti_deuteron,
              probability_three_to_two(type_pi, *type_anti_deuteron, dt,
                                       gcell_vol, spin_degn),
              ProcessType::MultiParticleThreeToTwo));
        }
 
        // Nnp → Nd, N̅np → N̅d
        if ((pdg_a.is_nucleon() && pdg_b == pdg::p && pdg_c == pdg::n) ||
            (pdg_a.is_nucleon() && pdg_b == pdg::n && pdg_c == pdg::p) ||
            (pdg_a == pdg::p && pdg_b.is_nucleon() && pdg_c == pdg::n) ||
            (pdg_a == pdg::n && pdg_b.is_nucleon() && pdg_c == pdg::p) ||
            (pdg_a == pdg::p && pdg_b == pdg::n && pdg_c.is_nucleon()) ||
            (pdg_a == pdg::n && pdg_b == pdg::p && pdg_c.is_nucleon())) {
          int symmetry_factor = 1;  // already true for N̅np → N̅d case
 
          ParticleList::iterator it =
              std::find_if(incoming_particles_.begin(),
                           incoming_particles_.end(), [](ParticleData x) {
                             return x.pdgcode().antiparticle_sign() == -1;
                           });
          if (it == incoming_particles_.end()) {
            /* Meaning no anti-N found by find_if,
             * therefore not N̅np → N̅d, but Nnp → Nd. */
            symmetry_factor = 2;  // for Nnp → Nd (2 factorial)
            // It is already clear here that we have a double of two N
            if (pdg_a == pdg_b) {
              it = incoming_particles_.begin();
            } else {
              // If a and b are not the double, then c has to be part of it
              it = incoming_particles_.begin() + 2;
            }
          }
          const ParticleType& type_N = it->type();
 
          const double spin_degn =
              react_degen_factor(spin_factor_inc, type_N.spin_degeneracy(),
                                 type_deuteron->spin_degeneracy());
 
          add_reaction(std::make_unique<CollisionBranch>(
              type_N, *type_deuteron,
              probability_three_to_two(type_N, *type_deuteron, dt, gcell_vol,
                                       symmetry_factor * spin_degn),
              ProcessType::MultiParticleThreeToTwo));
        }
 
        // Np̅n̅ → Nd̅, N̅p̅n̅ → N̅d̅
        if ((pdg_a.is_nucleon() && pdg_b == -pdg::p && pdg_c == -pdg::n) ||
            (pdg_a.is_nucleon() && pdg_b == -pdg::n && pdg_c == -pdg::p) ||
            (pdg_a == -pdg::p && pdg_b.is_nucleon() && pdg_c == -pdg::n) ||
            (pdg_a == -pdg::n && pdg_b.is_nucleon() && pdg_c == -pdg::p) ||
            (pdg_a == -pdg::p && pdg_b == -pdg::n && pdg_c.is_nucleon()) ||
            (pdg_a == -pdg::n && pdg_b == -pdg::p && pdg_c.is_nucleon())) {
          int symmetry_factor = 1;  // already true for Np̅n̅ → Nd̅ case
 
          ParticleList::iterator it =
              std::find_if(incoming_particles_.begin(),
                           incoming_particles_.end(), [](ParticleData x) {
                             return x.pdgcode().antiparticle_sign() == 1;
                           });
          if (it == incoming_particles_.end()) {
            /* Meaning no N found by find_if,
             * therefore not Np̅n̅ → Nd̅, but N̅p̅n̅ → N̅d̅. */
            symmetry_factor = 2;  // for N̅p̅n̅ → N̅d̅ (2 factorial)
            // It is already clear here that we have a double of two N̅
            if (pdg_a == pdg_b) {
              it = incoming_particles_.begin();
            } else {
              // If a and b are not the double, then c has to be part of it
              it = incoming_particles_.begin() + 2;
            }
          }
          const ParticleType& type_N = it->type();
 
          const double spin_degn =
              react_degen_factor(spin_factor_inc, type_N.spin_degeneracy(),
                                 type_anti_deuteron->spin_degeneracy());
 
          add_reaction(std::make_unique<CollisionBranch>(
              type_N, *type_anti_deuteron,
              probability_three_to_two(type_N, *type_anti_deuteron, dt,
                                       gcell_vol, symmetry_factor * spin_degn),
              ProcessType::MultiParticleThreeToTwo));
        }
      }
    }
  }
  // 4 -> 2
  if (incoming_particles_.size() == 4 &&
      incl_multi[IncludedMultiParticleReactions::A3_Nuclei_4to2] == 1) {
    std::map<PdgCode, int> c;  // counts incoming PdgCodes
    int spin_factor_inc = 1;
    for (const ParticleData& data : incoming_particles_) {
      c[data.pdgcode()]++;
      spin_factor_inc *= data.pdgcode().spin_degeneracy();
    }
    // Nucleons, antinucleons, and pions can catalyze
    const int n_possible_catalysts_incoming =
        c[pdg::n] + c[pdg::p] + c[-pdg::p] + c[-pdg::n] + c[pdg::pi_p] +
        c[pdg::pi_z] + c[pdg::pi_m];
 
    for (PdgCode pdg_nucleus :
         {pdg::triton, pdg::antitriton, pdg::he3, pdg::antihe3,
          pdg::hypertriton, pdg::antihypertriton}) {
      const ParticleTypePtr type_nucleus = ParticleType::try_find(pdg_nucleus);
      // Nucleus can be formed if and only if:
      // 1) Incoming particles contain enough components (like p, n, Lambda)
      // 2) In (incoming particles - components) there is still a catalyst
      // This is including the situation like nnpp. Can be that t(nnp) is formed
      // and p is catalyst, can be that he-3(ppn) is formed and n is catalyst.
      // Both reactions should be added.
      const int n_nucleus_components_that_can_be_catalysts =
          pdg_nucleus.nucleus_p() + pdg_nucleus.nucleus_ap() +
          pdg_nucleus.nucleus_n() + pdg_nucleus.nucleus_an();
      const bool incoming_contain_nucleus_components =
          c[pdg::p] >= pdg_nucleus.nucleus_p() &&
          c[-pdg::p] >= pdg_nucleus.nucleus_ap() &&
          c[pdg::n] >= pdg_nucleus.nucleus_n() &&
          c[-pdg::n] >= pdg_nucleus.nucleus_an() &&
          c[pdg::Lambda] >= pdg_nucleus.nucleus_La() &&
          c[-pdg::Lambda] >= pdg_nucleus.nucleus_aLa();
      const bool can_form_nucleus =
          type_nucleus && incoming_contain_nucleus_components &&
          n_possible_catalysts_incoming -
                  n_nucleus_components_that_can_be_catalysts ==
              1;
 
      if (!can_form_nucleus) {
        continue;
      }
      // Find the catalyst
      std::map<PdgCode, int> catalyst_count = c;
      catalyst_count[pdg::p] -= pdg_nucleus.nucleus_p();
      catalyst_count[-pdg::p] -= pdg_nucleus.nucleus_ap();
      catalyst_count[pdg::n] -= pdg_nucleus.nucleus_n();
      catalyst_count[-pdg::n] -= pdg_nucleus.nucleus_an();
      catalyst_count[pdg::Lambda] -= pdg_nucleus.nucleus_La();
      catalyst_count[-pdg::Lambda] -= pdg_nucleus.nucleus_aLa();
      PdgCode pdg_catalyst = PdgCode::invalid();
      for (const auto i : catalyst_count) {
        if (i.second == 1) {
          pdg_catalyst = i.first;
          break;
        }
      }
      if (pdg_catalyst == PdgCode::invalid()) {
        logg[LScatterActionMulti].error("Something went wrong while forming",
                                        pdg_nucleus, " from ",
                                        incoming_particles_);
      }
      const ParticleTypePtr type_catalyst =
          ParticleType::try_find(pdg_catalyst);
      const double spin_degn =
          react_degen_factor(spin_factor_inc, type_catalyst->spin_degeneracy(),
                             type_nucleus->spin_degeneracy());
      double symmetry_factor = 1.0;
      for (const auto i : c) {
        symmetry_factor *= (i.second == 3)   ? 6.0  // 3!
                           : (i.second == 2) ? 2.0  // 2!
                                             : 1.0;
        if (i.second > 3 || i.second < 0) {
          logg[LScatterActionMulti].error("4<->2 error, incoming particles ",
                                          incoming_particles_);
        }
      }
 
      add_reaction(std::make_unique<CollisionBranch>(
          *type_catalyst, *type_nucleus,
          probability_four_to_two(*type_catalyst, *type_nucleus, dt, gcell_vol,
                                  symmetry_factor * spin_degn),
          ProcessType::MultiParticleFourToTwo));
    }
  }
  // 5 -> 2
  if (incoming_particles_.size() == 5) {
    if (incl_multi[IncludedMultiParticleReactions::NNbar_5to2] == 1) {
      if (all_incoming_particles_are_pions_have_zero_charge_only_one_piz()) {
        const int spin_factor_inc =
            incoming_particles_[0].pdgcode().spin_degeneracy() *
            incoming_particles_[1].pdgcode().spin_degeneracy() *
            incoming_particles_[2].pdgcode().spin_degeneracy() *
            incoming_particles_[3].pdgcode().spin_degeneracy() *
            incoming_particles_[4].pdgcode().spin_degeneracy();
        const double symmetry_factor = 4.0;  // 2! * 2!
 
        const ParticleTypePtr type_p = ParticleType::try_find(pdg::p);
        const ParticleTypePtr type_anti_p = ParticleType::try_find(-pdg::p);
        const ParticleTypePtr type_n = ParticleType::try_find(pdg::n);
        const ParticleTypePtr type_anti_n = ParticleType::try_find(-pdg::n);
 
        const double spin_degn =
            react_degen_factor(spin_factor_inc, type_p->spin_degeneracy(),
                               type_anti_p->spin_degeneracy());
 
        if (type_p && type_n) {
          const double prob = probability_five_to_two(
              type_p->mass(), dt, gcell_vol,
              symmetry_factor * spin_degn);  // same for ppbar and nnbar
          add_reaction(std::make_unique<CollisionBranch>(
              *type_p, *type_anti_p, prob,
              ProcessType::MultiParticleFiveToTwo));
          add_reaction(std::make_unique<CollisionBranch>(
              *type_n, *type_anti_n, prob,
              ProcessType::MultiParticleFiveToTwo));
        }
      }
    }
  }
}

reaction_channels()

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

inline

Get list of possible reaction channels.

Returns

list of possible reaction channels.

Definition at line 72 of file scatteractionmulti.h.

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

sqrts

center of mass energy of incoming particles (= mass of outgoing particle)

Returns

result of integral

Definition at line 412 of file scatteractionmulti.cc.

                                                                {
  const double m1 = incoming_particles_[0].type().mass();
  const double m2 = incoming_particles_[1].type().mass();
  const double m3 = incoming_particles_[2].type().mass();
 
  if (sqrts < m1 + m2 + m3) {
    return 0.0;
  }
  const double x1 = (m1 - m2) * (m1 - m2), x2 = (m1 + m2) * (m1 + m2),
               x3 = (sqrts - m3) * (sqrts - m3),
               x4 = (sqrts + m3) * (sqrts + m3);
  const double qmm = x3 - x1, qmp = x3 - x2, qpm = x4 - x1, qpp = x4 - x2;
  const double kappa = std::sqrt(qpm * qmp / (qpp * qmm));
  const double tmp = std::sqrt(qmm * qpp);
  const double c1 =
      4.0 * m1 * m2 * std::sqrt(qmm / qpp) * (x4 - m3 * sqrts + m1 * m2);
  const double c2 = 0.5 * (m1 * m1 + m2 * m2 + m3 * m3 + sqrts * sqrts) * tmp;
  const double c3 = 8 * m1 * m2 / tmp *
                    ((m1 * m1 + m2 * m2) * (m3 * m3 + sqrts * sqrts) -
                     2 * m1 * m1 * m2 * m2 - 2 * m3 * m3 * sqrts * sqrts);
  const double c4 =
      -8 * m1 * m2 / tmp * smash::pow_int(sqrts * sqrts - m3 * m3, 2);
  const double res =
      c1 * gsl_sf_ellint_Kcomp(kappa, GSL_PREC_DOUBLE) +
      c2 * gsl_sf_ellint_Ecomp(kappa, GSL_PREC_DOUBLE) +
      c3 * gsl_sf_ellint_Pcomp(kappa, -qmp / qmm, GSL_PREC_DOUBLE) +
      c4 * gsl_sf_ellint_Pcomp(kappa, -x1 * qmp / (x2 * qmm), GSL_PREC_DOUBLE);
  return res;
}

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_s

mandelstam s of reaction

Returns

phase space integral value for 4 bodies [GeV^4]

Definition at line 509 of file scatteractionmulti.cc.

                                                                       {
  int n_nucleons = 0, n_pions = 0, n_lambdas = 0;
  double sum_m = 0.0, prod_m = 1.0;
  for (const ParticleData& data : incoming_particles_) {
    const PdgCode pdg = data.type().pdgcode();
    n_nucleons += pdg.is_nucleon();  // including anti-nucleons
    n_pions += pdg.is_pion();
    n_lambdas += pdg.is_Lambda();  // including anti-Lambda
    sum_m += data.type().mass();
    prod_m *= data.type().mass();
  }
  const double x = 1.0 - sum_m / std::sqrt(man_s);
  const double x2 = x * x;
  const double x3 = x2 * x;
  double g = -1.0;
 
  if (n_nucleons == 3 && n_pions == 1) {  // NNNpi
    g = (1.0 + 0.862432 * x - 3.4853 * x2 + 1.70259 * x3) /
        (1.0 + 0.387376 * x - 1.34128 * x2 + 0.154489 * x3);
  } else if (n_nucleons == 4) {  // NNNN
    g = (1.0 - 1.72285 * x + 0.728331 * x2) /
        (1.0 - 0.967146 * x - 0.0103633 * x2);
  } else if (n_nucleons == 2 && n_lambdas == 1 && n_pions == 1) {  // LaNNpi
    g = (1.0 + 0.937064 * x - 3.56864 * x2 + 1.721 * x3) /
        (1.0 + 0.365202 * x - 1.2854 * x2 + 0.138444 * x3);
  } else if (n_nucleons == 3 && n_lambdas == 1) {  // LaNNN
    g = (1.0 + 0.882401 * x - 3.4074 * x2 + 1.62454 * x3) /
        (1.0 + 1.61741 * x - 2.12543 * x2 - 0.0902067 * x3);
  }
 
  if (g > 0.0) {
    return (std::sqrt(prod_m) * sum_m * sum_m * std::pow(x, 3.5) * g) /
           (840. * std::sqrt(2) * std::pow(M_PI, 4.0) * std::pow(1 - x, 4.0));
  } else {
    logg[LScatterActionMulti].error("parametrizaton_phi4: no parametrization ",
                                    "available for ", incoming_particles_);
    return 0.0;
  }
}

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]

out

out stream to be written to

Implements smash::Action.

Definition at line 656 of file scatteractionmulti.cc.

                                                                  {
  out << "MultiParticleScatter of " << incoming_particles_;
  if (outgoing_particles_.empty()) {
    out << " (not performed)";
  } else {
    out << " to " << outgoing_particles_;
  }
}

add_reaction()

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

private

Add a new reaction channel.

Parameters

[in]

p

Channel to be added.

Definition at line 29 of file scatteractionmulti.cc.

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

add_reactions()

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

private

Add several new reaction channels at once.

Parameters

[in]

pv

list of channels to be added.

Definition at line 33 of file scatteractionmulti.cc.

                                                             {
  add_processes<CollisionBranch>(std::move(pv), reaction_channels_,
                                 total_probability_);
}

annihilation()

void smash::ScatterActionMulti::annihilation ( )

private

Perform a n->1 annihilation process.

Exceptions

InvalidScatterActionMulti

Definition at line 580 of file scatteractionmulti.cc.

                                      {
  if (outgoing_particles_.size() != 1) {
    std::string s =
        "Annihilation: "
        "Incorrect number of particles in final state: ";
    s += std::to_string(outgoing_particles_.size()) + ".";
    throw InvalidScatterActionMulti(s);
  }
  // Set the momentum of the formed particle in its rest frame.
  outgoing_particles_[0].set_4momentum(
      total_momentum_of_outgoing_particles().abs(), 0., 0., 0.);
  // Make sure to assign formation times before boost to the computational frame
  assign_formation_time_to_outgoing_particles();
 
  logg[LScatterActionMulti].debug("Momentum of the new particle: ",
                                  outgoing_particles_[0].momentum());
}

n_to_two()

void smash::ScatterActionMulti::n_to_two ( )

private

Perform a n->2 process.

Definition at line 598 of file scatteractionmulti.cc.

                                  {
  sample_2body_phasespace();
  // Make sure to assign formation times before boost to the computational frame
  assign_formation_time_to_outgoing_particles();
  logg[LScatterActionMulti].debug(incoming_particles_.size(),
                                  "->2 scattering:", incoming_particles_,
                                  " -> ", outgoing_particles_);
}

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

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_out	type of outgoing particle
[in]	dt	timestep size
[in]	gcell_vol	grid cell volume
[in]	degen_sym_factor	degeneracy factor for reaction (including symmetry factors)

Returns

probabilty for 3-to-1 reaction

Definition at line 442 of file scatteractionmulti.cc.

                                      {
  const double e1 = incoming_particles_[0].momentum().x0();
  const double e2 = incoming_particles_[1].momentum().x0();
  const double e3 = incoming_particles_[2].momentum().x0();
  const double sqrts = sqrt_s();
 
  const double gamma_decay = type_out.get_partial_width(
      sqrts, {&incoming_particles_[0].type(), &incoming_particles_[1].type(),
              &incoming_particles_[2].type()});
 
  const double I_3 = calculate_I3(sqrts);
  const double ph_sp_3 =
      1. / (8 * M_PI * M_PI * M_PI) * 1. / (16 * sqrts * sqrts) * I_3;
 
  const double spec_f_val = type_out.spectral_function(sqrts);
 
  return dt / (gcell_vol * gcell_vol) * M_PI / (4. * e1 * e2 * e3) *
         gamma_decay / ph_sp_3 * spec_f_val * std::pow(hbarc, 5.0) *
         degen_sym_factor;
}

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

\[ 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_out1	type of outgoing particle 1
[in]	type_out2	type of outgoing particle 2
[in]	dt	timestep size
[in]	gcell_vol	grid cell volume
[in]	degen_sym_factor	degeneracy factor for reaction (including symmetry factors)

Returns

probabilty for 3-to-2 reaction

Definition at line 465 of file scatteractionmulti.cc.

                                                                 {
  const double e1 = incoming_particles_[0].momentum().x0();
  const double e2 = incoming_particles_[1].momentum().x0();
  const double e3 = incoming_particles_[2].momentum().x0();
  const double m4 = type_out1.mass();
  const double m5 = type_out2.mass();
 
  const double sqrts = sqrt_s();
  const double xs =
      CrossSections::two_to_three_xs(type_out1, type_out2, sqrts) / gev2_mb;
  const double lamb = lambda_tilde(sqrts * sqrts, m4 * m4, m5 * m5);
 
  const double I_3 = calculate_I3(sqrts);
  const double ph_sp_3 =
      1. / (8 * M_PI * M_PI * M_PI) * 1. / (16 * sqrts * sqrts) * I_3;
 
  return dt / (gcell_vol * gcell_vol) * 1. / (4. * e1 * e2 * e3) * lamb /
         (ph_sp_3 * 8 * M_PI * sqrts * sqrts) * xs * std::pow(hbarc, 5.0) *
         degen_sym_factor;
}

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

\[ 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_out1	type of outgoing particle 1
[in]	type_out2	type of outgoing particle 2
[in]	dt	timestep size
[in]	gcell_vol	grid cell volume
[in]	degen_sym_factor	degeneracy factor for reaction (including symmetry factors)

Returns

probabilty for 4-to-2 reaction

Definition at line 488 of file scatteractionmulti.cc.

                                                                 {
  const double e1 = incoming_particles_[0].momentum().x0();
  const double e2 = incoming_particles_[1].momentum().x0();
  const double e3 = incoming_particles_[2].momentum().x0();
  const double e4 = incoming_particles_[3].momentum().x0();
  const double m5 = type_out1.mass();
  const double m6 = type_out2.mass();
 
  const double man_s = sqrt_s() * sqrt_s();
  const double xs =
      CrossSections::two_to_four_xs(type_out1, type_out2, sqrt_s()) / gev2_mb;
  const double lamb = lambda_tilde(man_s, m5 * m5, m6 * m6);
  const double ph_sp_4 = parametrizaton_phi4(man_s);
 
  return dt / std::pow(gcell_vol, 3.0) * 1. / (16. * e1 * e2 * e3 * e4) * xs /
         (4. * M_PI * man_s) * lamb / ph_sp_4 * std::pow(hbarc, 8.0) *
         degen_sym_factor;
}

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_out	mass of outgoing particle types (assumes equal masses)
[in]	dt	timestep size
[in]	gcell_vol	grid cell volume
[in]	degen_sym_factor	degeneracy factor for reaction (including symmetry factors)

Returns

probabilty for 5-to-2 reaction

Definition at line 549 of file scatteractionmulti.cc.

                                         {
  const double e1 = incoming_particles_[0].momentum().x0();
  const double e2 = incoming_particles_[1].momentum().x0();
  const double e3 = incoming_particles_[2].momentum().x0();
  const double e4 = incoming_particles_[3].momentum().x0();
  const double e5 = incoming_particles_[4].momentum().x0();
 
  const double man_s = sqrt_s() * sqrt_s();
  const double lamb = lambda_tilde(man_s, mout * mout, mout * mout);
  const double ph_sp_5 = parametrizaton_phi5_pions(man_s);
 
  // Matching the NNbar anihilation cross section defintion for 2-to-5
  const double xs =
      std::max(0., ppbar_total(man_s) - ppbar_elastic(man_s)) / gev2_mb;
 
  return dt / std::pow(gcell_vol, 4.0) * 1. / (32. * e1 * e2 * e3 * e4 * e5) *
         xs / (4. * M_PI * man_s) * lamb / ph_sp_5 * std::pow(hbarc, 11.0) *
         degen_sym_factor;
}

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_s

mandelstam s of reaction

Returns

phase space value for 5 pions

Definition at line 571 of file scatteractionmulti.cc.

                                                                             {
  // see function documentation for parameter naming
  const double s_zero = 25 * pion_mass * pion_mass;
  const double fit_a = 2.1018e-10;
  const double fit_alpha = 1.982;
  return fit_a * std::pow(man_s - s_zero, 5.0) *
         std::pow(1 + man_s / s_zero, -fit_alpha);
}

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_inc	product of incoming spin degeneracy (denominator in above expression)
[in]	spin_degen_out1	degeneracy factor of outgoing particle 1
[in]	spin_degen_out2	degeneracy factor of outgoing particle 2

Returns

spin degeneracy factor

Definition at line 283 of file scatteractionmulti.h.

                                                             {
    return static_cast<double>(spin_degen_out1 * spin_degen_out2) /
           static_cast<double>(spin_factor_inc);
  }

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_a	data for first incoming particle
[in]	data_b	data for second incoming particle
[in]	data_c	data for third incoming particle

Returns

true if combination of π⁺,π⁻,π⁰

Definition at line 607 of file scatteractionmulti.cc.

                                      {
  // We want a combination of pi+, pi- and pi0
  const PdgCode pdg_a = data_a.pdgcode();
  const PdgCode pdg_b = data_b.pdgcode();
  const PdgCode pdg_c = data_c.pdgcode();
 
  return (pdg_a.is_pion() && pdg_b.is_pion() && pdg_c.is_pion()) &&
         (pdg_a != pdg_b && pdg_b != pdg_c && pdg_c != pdg_a);
}

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_a	data for first incoming particle
[in]	data_b	data for second incoming particle
[in]	data_c	data for third incoming particle

Returns

true if combination of π⁺,π⁻,η or π⁰,π⁰,η

Definition at line 619 of file scatteractionmulti.cc.

                                                                         {
  // We want a combination of pi0, pi0 and eta or pi+, pi- and eta
  const PdgCode pdg_a = data_a.pdgcode();
  const PdgCode pdg_b = data_b.pdgcode();
  const PdgCode pdg_c = data_c.pdgcode();
 
  return (pdg_a == pdg::pi_z && pdg_b == pdg::pi_z && pdg_c == pdg::eta) ||
         (pdg_a == pdg::pi_z && pdg_b == pdg::eta && pdg_c == pdg::pi_z) ||
         (pdg_a == pdg::eta && pdg_b == pdg::pi_z && pdg_c == pdg::pi_z) ||
 
         (pdg_a == pdg::eta && pdg_b == pdg::pi_m && pdg_c == pdg::pi_p) ||
         (pdg_a == pdg::eta && pdg_b == pdg::pi_p && pdg_c == pdg::pi_m) ||
         (pdg_a == pdg::pi_m && pdg_b == pdg::pi_p && pdg_c == pdg::eta) ||
         (pdg_a == pdg::pi_m && pdg_b == pdg::eta && pdg_c == pdg::pi_p) ||
         (pdg_a == pdg::pi_p && pdg_b == pdg::pi_m && pdg_c == pdg::eta) ||
         (pdg_a == pdg::pi_p && pdg_b == pdg::eta && pdg_c == pdg::pi_m);
}

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 639 of file scatteractionmulti.cc.

                                                                           {
  const bool all_inc_pi =
      all_of(incoming_particles_.begin(), incoming_particles_.end(),
             [](const ParticleData& data) { return data.is_pion(); });
  const int no_of_piz = std::count_if(
      incoming_particles_.begin(), incoming_particles_.end(),
      [](const ParticleData& data) { return data.pdgcode() == pdg::pi_z; });
 
  int total_state_charge = 0;
  for (const ParticleData& part : incoming_particles_) {
    total_state_charge += part.pdgcode().charge();
  }
 
  return (all_inc_pi && total_state_charge == 0 && no_of_piz == 1);
}

Member Data Documentation

total_probability_

double smash::ScatterActionMulti::total_probability_

private

Total probability of reaction.

Definition at line 322 of file scatteractionmulti.h.

partial_probability_

double smash::ScatterActionMulti::partial_probability_

private

Partial probability of the chosen outgoing channel.

Definition at line 325 of file scatteractionmulti.h.

reaction_channels_

CollisionBranchList smash::ScatterActionMulti::reaction_channels_

private

List of possible collisions.

Definition at line 328 of file scatteractionmulti.h.

---

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

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