#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:

[legend]

Collaboration diagram for smash::ScatterActionMulti:

[legend]

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

Public Member Functions inherited from smash::Action
	Action (const ParticleList &in_part, double time)
	Construct an action object with incoming particles and relative time. More...

	Action (const ParticleData &in_part, const ParticleData &out_part, double time, ProcessType type)
	Construct an action object with the incoming particles, relative time, and the already known outgoing particles and type of the process. More...

	Action (const ParticleList &in_part, const ParticleList &out_part, double absolute_execution_time, ProcessType type)
	Construct an action object with the incoming particles, absolute time, and the already known outgoing particles and type of the process. More...

	Action (const Action &)=delete
	Copying is disabled. Use pointers or create a new Action. More...

virtual	~Action ()
	Virtual Destructor. More...

bool	operator< (const Action &rhs) const
	Determine whether one action takes place before another in time. More...

virtual ProcessType	get_type () const
	Get the process type. More...

template<typename Branch >
void	add_process (ProcessBranchPtr< Branch > &p, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add a new subprocess. More...

template<typename Branch >
void	add_processes (ProcessBranchList< Branch > pv, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add several new subprocesses at once. More...

virtual void	perform (Particles *particles, uint32_t id_process)
	Actually perform the action, e.g. More...

bool	is_valid (const Particles &particles) const
	Check whether the action still applies. More...

bool	is_pauli_blocked (const std::vector< Particles > &ensembles, const PauliBlocker &p_bl) const
	Check if the action is Pauli-blocked. More...

const ParticleList &	incoming_particles () const
	Get the list of particles that go into the action. More...

void	update_incoming (const Particles &particles)
	Update the incoming particles that are stored in this action to the state they have in the global particle list. More...

const ParticleList &	outgoing_particles () const
	Get the list of particles that resulted from the action. More...

double	time_of_execution () const
	Get the time at which the action is supposed to be performed. More...

virtual void	check_conservation (const uint32_t id_process) const
	Check various conservation laws. More...

double	sqrt_s () const
	Determine the total energy in the center-of-mass frame [GeV]. More...

FourVector	total_momentum_of_outgoing_particles () const
	Calculate the total kinetic momentum of the outgoing particles. More...

FourVector	get_interaction_point () const
	Get the interaction point. More...

std::pair< FourVector, 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_3body_phasespace ()
	Sample the full 3-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...

virtual void	sample_5body_phasespace ()
	Sample the full 5-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	three_to_two ()
	Perform a 3->2 process. More...

void	five_to_two ()
	Perform a 5->2 process. More...

double	probability_three_to_one (const ParticleType &type_out, double dt, const double gcell_vol, const int degen_factor=1) const
	Calculate the probability for a 3m-to-1 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [49]. More...

double	probability_three_to_two (const ParticleType &type_out1, const ParticleType &type_out2, double dt, const double gcell_vol, const double degen_factor=1.0) const
	Calculate the probability for a 3-to-2 reaction according to the stochastic collision criterion as given in Staudenmaier:2021lrg [49]. More...

double	probability_five_to_two (const double m_out, double dt, const double gcell_vol, const double degen_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 [20]. More...

double	parametrizaton_phi5_pions (const double man_s) const
	Calculate the parametrized 5-pion phase space. More...

double	calculate_I3 (const double sqrts) const
	Calculate the integration necessary for the three-body 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 3->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...

Protected Attributes inherited from smash::Action
ParticleList	incoming_particles_
	List with data of incoming particles. More...

ParticleList	outgoing_particles_
	Initially this stores only the PDG codes of final-state particles. More...

const double	time_of_execution_
	Time at which the action is supposed to be performed (absolute time in the lab frame in fm/c). More...

ProcessType	process_type_
	type of process More...

double	box_length_ = -1.0
	Box length: needed to determine coordinates of collision correctly in case of collision through the wall. More...

int	stochastic_position_idx_ = -1
	This stores a randomly-chosen index to an incoming particle. More...

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

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

smash::Action::Action

Action(const ParticleList &in_part, double time)

Construct an action object with incoming particles and relative time.

Definition: action.h:44

smash::ScatterActionMulti::total_probability_

double total_probability_

Total probability of reaction.

Definition: scatteractionmulti.h:280

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 278 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:
       /* 3->2 scattering */
       three_to_two();
       break;
     case ProcessType::MultiParticleFiveToTwo:
       /* 5->2 scattering */
       five_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);
   }
 }

Here is the call graph for this function:

◆ 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 33 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 41 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 49 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(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(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(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(PdgCode::from_decimal(pdg::decimal_d));
       const ParticleTypePtr type_anti_deuteron =
           ParticleType::try_find(PdgCode::from_decimal(pdg::decimal_antid));
  
       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(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(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(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(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));
         }
       }
     }
   }
   // 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(make_unique<CollisionBranch>(
               *type_p, *type_anti_p, prob,
               ProcessType::MultiParticleFiveToTwo));
           add_reaction(make_unique<CollisionBranch>(
               *type_n, *type_anti_n, prob,
               ProcessType::MultiParticleFiveToTwo));
         }
       }
     }
   }
 }

Here is the call graph for this function:

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

72 { return reaction_channels_; }

◆ 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 510 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 24 of file scatteractionmulti.cc.

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

Here is the caller graph for this function:

◆ 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 28 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 427 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());
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ three_to_two()

void smash::ScatterActionMulti::three_to_two ( )

private

Perform a 3->2 process.

Definition at line 445 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("3->2 scattering:", incoming_particles_,
                                   " -> ", outgoing_particles_);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ five_to_two()

void smash::ScatterActionMulti::five_to_two ( )

private

Perform a 5->2 process.

Definition at line 453 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("5->2 scattering:", incoming_particles_,
                                   " -> ", outgoing_particles_);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ probability_three_to_one()

double smash::ScatterActionMulti::probability_three_to_one	(	const ParticleType &	type_out,
		double	dt,
		const double	gcell_vol,
		const int	degen_factor = `1`
	)		const

private

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

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_factor	degeneracy factor for reaction (including symmetry factors)

Returns: probabilty for 3-to-1 reaction

Definition at line 350 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_factor;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

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

\[ 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 1
[in]	dt	timestep size
[in]	gcell_vol	grid cell volume
[in]	degen_factor	degeneracy factor for reaction (including symmetry factors)

Returns: probabilty for 3-to-2 reaction

Definition at line 373 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_factor;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ probability_five_to_two()

double smash::ScatterActionMulti::probability_five_to_two	(	const double	m_out,
		double	dt,
		const double	gcell_vol,
		const double	degen_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 [20].

\[ 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_factor	degeneracy factor for reaction (including symmetry factors)

Returns: probabilty for 5-to-2 reaction

Definition at line 396 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_factor;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ 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 418 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);
 }

Here is the caller graph for this function:

◆ calculate_I3()

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

private

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

                                                                 {
   static Integrator integrate;
   const double m1 = incoming_particles_[0].effective_mass();
   const double m2 = incoming_particles_[1].effective_mass();
   const double m3 = incoming_particles_[2].effective_mass();
   const double lower_bound = (m1 + m2) * (m1 + m2);
   const double upper_bound = (sqrts - m3) * (sqrts - m3);
   const auto result = integrate(lower_bound, upper_bound, [&](double m12_sqr) {
     const double m12 = std::sqrt(m12_sqr);
     const double e2_star = (m12_sqr - m1 * m1 + m2 * m2) / (2 * m12);
     const double e3_star = (sqrts * sqrts - m12_sqr - m3 * m3) / (2 * m12);
     const double m23_sqr_min =
         (e2_star + e3_star) * (e2_star + e3_star) -
         std::pow(std::sqrt(e2_star * e2_star - m2 * m2) +
                      std::sqrt(e3_star * e3_star - m3 * m3),
                  2.0);
     const double m23_sqr_max =
         (e2_star + e3_star) * (e2_star + e3_star) -
         std::pow(std::sqrt(e2_star * e2_star - m2 * m2) -
                      std::sqrt(e3_star * e3_star - m3 * m3),
                  2.0);
     return m23_sqr_max - m23_sqr_min;
   });
  
   return result;
 }

Here is the caller graph for this function:

◆ 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 3->2 reaction.

\[D_{spin} = \frac{(2J_{out1}+1)(2J_{out2}+1)} {(2J_{in1}+1)(2J_{in2}+1)(2J_{in3}+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 241 of file scatteractionmulti.h.

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

Here is the caller graph for this function:

◆ 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 461 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ 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 473 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ 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 493 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);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

◆ total_probability_

double smash::ScatterActionMulti::total_probability_

private

Total probability of reaction.

Definition at line 280 of file scatteractionmulti.h.

◆ partial_probability_

double smash::ScatterActionMulti::partial_probability_

private

Partial probability of the chosen outgoing channel.

Definition at line 283 of file scatteractionmulti.h.

◆ reaction_channels_

CollisionBranchList smash::ScatterActionMulti::reaction_channels_

private

List of possible collisions.

Definition at line 286 of file scatteractionmulti.h.

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

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

Classes

Public Member Functions

Protected Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Constructor & Destructor Documentation

◆ ScatterActionMulti()

Member Function Documentation

◆ generate_final_state()

◆ get_total_weight()

◆ get_partial_weight()

◆ add_possible_reactions()

◆ reaction_channels()

◆ format_debug_output()

◆ add_reaction()

◆ add_reactions()

◆ annihilation()

◆ three_to_two()

◆ five_to_two()

◆ probability_three_to_one()

◆ probability_three_to_two()

◆ probability_five_to_two()

◆ parametrizaton_phi5_pions()

◆ calculate_I3()

◆ react_degen_factor()

◆ three_different_pions()

◆ two_pions_eta()

◆ all_incoming_particles_are_pions_have_zero_charge_only_one_piz()

Member Data Documentation

◆ total_probability_

◆ partial_probability_

◆ reaction_channels_