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:

Collaboration diagram for smash::ScatterActionMulti:

Classes
class	InvalidScatterActionMulti

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 Particles &particles, const PauliBlocker &p_bl) const
	Check if the action is Pauli-blocked. More...
const ParticleList &	incoming_particles () const
	Get the list of particles that go into the action. More...
void	update_incoming (const Particles &particles)
	Update the incoming particles that are stored in this action to the state they have in the global particle list. More...
const ParticleList &	outgoing_particles () const
	Get the list of particles that resulted from the action. More...
double	time_of_execution () const
	Get the time at which the action is supposed to be performed. More...
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
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...
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...
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 (e.g. 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 (similar to Cassing:2001ds [14]). 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...

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

                                              {
  logg[LScatterActionMulti].debug("Incoming particles: ", 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;
    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 32 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 40 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 48 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));
        }
      }
    }
  }
}

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.

{ 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 413 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 23 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 27 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 355 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 373 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:

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 (e.g.

Xu:2004mz [53] (Sec.IIB)).

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 309 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 (similar to Cassing:2001ds [14]).

\[ P_{3 \rightarrow 2} = \frac{1}{4E_1E_2E_3} \frac{\Delta t}{\Delta^3 x} \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 332 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:

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 282 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 200 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 381 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 393 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:

Member Data Documentation

total_probability_

double smash::ScatterActionMulti::total_probability_

private

Total probability of reaction.

Definition at line 233 of file scatteractionmulti.h.

partial_probability_

double smash::ScatterActionMulti::partial_probability_

private

Partial probability of the chosen outgoing channel.

Definition at line 236 of file scatteractionmulti.h.

reaction_channels_

CollisionBranchList smash::ScatterActionMulti::reaction_channels_

private

List of possible collisions.

Definition at line 239 of file scatteractionmulti.h.

---

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

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