smash::DecayAction Class Reference

#include <decayaction.h>

DecayAction is a special action which takes one single particle in the initial state and makes it decay into a number of daughter particles (currently two or three).

Definition at line 25 of file decayaction.h.

Inheritance diagram for smash::DecayAction:

Collaboration diagram for smash::DecayAction:

Classes
class	InvalidDecay
	Thrown when DecayAction is called to perform with 0 or more than 2 entries in outgoing_particles. More...

Public Member Functions
	DecayAction (const ParticleData &p, double time)
	Construct a DecayAction from a particle p. More...
void	add_decays (DecayBranchList pv)
	Add several new decays at once. More...
void	add_decay (DecayBranchPtr p)
	Add one new decay. More...
void	generate_final_state () override
	Generate the final state of the decay process. More...
std::pair< double, double >	sample_masses (double kinetic_energy_cm) const override
	Sample the masses of the final particles. More...
double	get_total_weight () const override
	Return the total width of the decay process. More...
double	get_partial_weight () const override
	Get partial width of chosen channel. More...
double	total_width () const
	Get total decay width. More...
Public Member Functions inherited from smash::Action
	Action (const ParticleList &in_part, double time)
	Construct an action object with incoming particles and relative time. More...
	Action (const ParticleData &in_part, const ParticleData &out_part, double time, ProcessType type)
	Construct an action object with the incoming particles, relative time, and the already known outgoing particles and type of the process. More...
	Action (const ParticleList &in_part, const ParticleList &out_part, double absolute_execution_time, ProcessType type)
	Construct an action object with the incoming particles, absolute time, and the already known outgoing particles and type of the process. More...
	Action (const Action &)=delete
	Copying is disabled. Use pointers or create a new Action. More...
virtual	~Action ()
	Virtual Destructor. More...
bool	operator< (const Action &rhs) const
	Determine whether one action takes place before another in time. More...
virtual ProcessType	get_type () const
	Get the process type. More...
template<typename Branch >
void	add_process (ProcessBranchPtr< Branch > &p, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add a new subprocess. More...
template<typename Branch >
void	add_processes (ProcessBranchList< Branch > pv, ProcessBranchList< Branch > &subprocesses, double &total_weight)
	Add several new subprocesses at once. More...
virtual void	perform (Particles *particles, uint32_t id_process)
	Actually perform the action, e.g. More...
bool	is_valid (const Particles &particles) const
	Check whether the action still applies. More...
bool	is_pauli_blocked (const Particles &particles, const PauliBlocker &p_bl) const
	Check if the action is Pauli-blocked. More...
const ParticleList &	incoming_particles () const
	Get the list of particles that go into the action. More...
void	update_incoming (const Particles &particles)
	Update the incoming particles that are stored in this action to the state they have in the global particle list. More...
const ParticleList &	outgoing_particles () const
	Get the list of particles that resulted from the action. More...
double	time_of_execution () const
	Get the time at which the action is supposed to be performed. More...
void	check_conservation (const uint32_t id_process) const
	Check various conservation laws. More...
double	sqrt_s () const
	Determine the total energy in the center-of-mass frame [GeV]. More...
FourVector	total_momentum_of_outgoing_particles () const
	Calculate the total kinetic momentum of the outgoing particles. More...
FourVector	get_interaction_point () const
	Get the interaction point. More...
std::pair< FourVector, FourVector >	get_potential_at_interaction_point () const
	Get the skyrme and asymmetry potential at the interaction point. More...

Protected Member Functions
void	format_debug_output (std::ostream &out) const override
	Writes information about this decay action to the out stream. More...
virtual void	one_to_three ()
	Kinematics of a 1-to-3 decay process. 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 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...

Protected Attributes
DecayBranchList	decay_channels_
	List of possible decays. More...
double	total_width_
	total decay width More...
double	partial_width_
	partial decay width to the chosen outgoing channel More...
int	L_ = 0
	Angular momentum of the decay. 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...

Constructor & Destructor Documentation

DecayAction()

smash::DecayAction::DecayAction	(	const ParticleData &	p,
		double	time
	)

Construct a DecayAction from a particle p.

It does not initialize the list of possible decay processes. You need to call add_processes after construction.

Parameters

[in]	p	The particle that should decay if the action is performed.
[in]	time	Time at which the action is supposed to take place

Definition at line 21 of file decayaction.cc.

    : Action({p}, time), total_width_(0.) {}

Member Function Documentation

add_decays()

void smash::DecayAction::add_decays

(

DecayBranchList

pv

)

Add several new decays at once.

Parameters

[in]

pv

List of decays to be added.

Definition at line 24 of file decayaction.cc.

                                               {
  add_processes<DecayBranch>(std::move(pv), decay_channels_, total_width_);
}

add_decay()

void smash::DecayAction::add_decay

(

DecayBranchPtr

p

)

Add one new decay.

Parameters

[in]

p

Decay to be added.

Definition at line 28 of file decayaction.cc.

                                            {
  add_process<DecayBranch>(p, decay_channels_, total_width_);
}

Here is the caller graph for this function:

generate_final_state()

void smash::DecayAction::generate_final_state ( )

overridevirtual

Generate the final state of the decay process.

Performs a decay of one particle to two or three particles.

Exceptions

InvalidDecay

Implements smash::Action.

Definition at line 149 of file decayaction.cc.

                                       {
  const auto &log = logger<LogArea::DecayModes>();
  log.debug("Process: Resonance decay. ");
  /* Execute a decay process for the selected particle.
   *
   * randomly select one of the decay modes of the particle
   * according to their relative weights. Then decay the particle
   * by calling function one_to_two or one_to_three.
   */
  const DecayBranch *proc =
      choose_channel<DecayBranch>(decay_channels_, total_width_);
  outgoing_particles_ = proc->particle_list();
  // set positions of the outgoing particles
  for (auto &p : outgoing_particles_) {
    p.set_4position(incoming_particles_[0].position());
  }
  process_type_ = proc->get_type();
  L_ = proc->angular_momentum();
  partial_width_ = proc->weight();

  switch (outgoing_particles_.size()) {
    case 2:
      sample_2body_phasespace();
      break;
    case 3:
      one_to_three();
      break;
    default:
      throw InvalidDecay(
          "DecayAction::perform: Only 1->2 or 1->3 processes are supported. "
          "Decay from 1->" +
          std::to_string(outgoing_particles_.size()) +
          " was requested. (PDGcode=" +
          incoming_particles_[0].pdgcode().string() + ", mass=" +
          std::to_string(incoming_particles_[0].effective_mass()) + ")");
  }

  // Set formation time.
  for (auto &p : outgoing_particles_) {
    log.debug("particle momenta in lrf ", p);
    // assuming decaying particles are always fully formed
    p.set_formation_time(time_of_execution_);
    // Boost to the computational frame
    p.boost_momentum(-total_momentum_of_outgoing_particles().velocity());
    log.debug("particle momenta in comp ", p);
  }
}

Here is the call graph for this function:

Here is the caller graph for this function:

sample_masses()

std::pair< double, double > smash::DecayAction::sample_masses ( double  kinetic_energy_cm ) const

overridevirtual

Sample the masses of the final particles.

Returns

Pair of sampled masses of particle 1 and 2

Reimplemented from smash::Action.

Definition at line 199 of file decayaction.cc.

                                    {
  const ParticleType &t_a = outgoing_particles_[0].type();
  const ParticleType &t_b = outgoing_particles_[1].type();

  // start with pole masses
  std::pair<double, double> masses = {t_a.mass(), t_b.mass()};

  if (kinetic_energy_cm < t_a.min_mass_kinematic() + t_b.min_mass_kinematic()) {
    const std::string reaction =
        incoming_particles_[0].type().name() + "→" + t_a.name() + t_b.name();
    throw InvalidResonanceFormation(
        reaction + ": not enough energy, " + std::to_string(kinetic_energy_cm) +
        " < " + std::to_string(t_a.min_mass_kinematic()) + " + " +
        std::to_string(t_b.min_mass_kinematic()));
  }

  // If one of the particles is a resonance, sample its mass.
  if (!t_a.is_stable() && t_b.is_stable()) {
    masses.first = t_a.sample_resonance_mass(t_b.mass(), kinetic_energy_cm, L_);
  } else if (!t_b.is_stable() && t_a.is_stable()) {
    masses.second =
        t_b.sample_resonance_mass(t_a.mass(), kinetic_energy_cm, L_);
  } else if (!t_a.is_stable() && !t_b.is_stable()) {
    // two resonances in final state
    masses = t_a.sample_resonance_masses(t_b, kinetic_energy_cm, L_);
  }

  return masses;
}

Here is the call graph for this function:

get_total_weight()

double smash::DecayAction::get_total_weight ( ) const

inlineoverridevirtual

Return the total width of the decay process.

Implements smash::Action.

Reimplemented in smash::DecayActionDilepton.

Definition at line 66 of file decayaction.h.

{ return total_width_; }

get_partial_weight()

double smash::DecayAction::get_partial_weight ( ) const

inlineoverridevirtual

Get partial width of chosen channel.

Returns

Partial width of chosen channel

Implements smash::Action.

Definition at line 72 of file decayaction.h.

{ return partial_width_; }

total_width()

double smash::DecayAction::total_width ( ) const

inline

Get total decay width.

Returns

Total width of decay

Definition at line 78 of file decayaction.h.

{ return total_width_; }

one_to_three()

void smash::DecayAction::one_to_three ( )

protectedvirtual

Kinematics of a 1-to-3 decay process.

Sample the masses and momenta of the decay products in the center-of-momentum frame.

Reimplemented in smash::DecayActionDilepton.

Definition at line 32 of file decayaction.cc.

                               {
  const auto &log = logger<LogArea::DecayModes>();
  ParticleData &outgoing_a = outgoing_particles_[0];
  ParticleData &outgoing_b = outgoing_particles_[1];
  ParticleData &outgoing_c = outgoing_particles_[2];
  const ParticleType &outgoing_a_type = outgoing_a.type();
  const ParticleType &outgoing_b_type = outgoing_b.type();
  const ParticleType &outgoing_c_type = outgoing_c.type();

  log.debug("Note: Doing 1->3 decay!");

  const double mass_a = outgoing_a_type.mass();
  const double mass_b = outgoing_b_type.mass();
  const double mass_c = outgoing_c_type.mass();
  const double mass_resonance = incoming_particles_[0].effective_mass();

  // mandelstam-s limits for pairs ab and bc
  const double s_ab_max = (mass_resonance - mass_c) * (mass_resonance - mass_c);
  const double s_ab_min = (mass_a + mass_b) * (mass_a + mass_b);
  const double s_bc_max = (mass_resonance - mass_a) * (mass_resonance - mass_a);
  const double s_bc_min = (mass_b + mass_c) * (mass_b + mass_c);

  log.debug("s_ab limits: ", s_ab_min, " ", s_ab_max);
  log.debug("s_bc limits: ", s_bc_min, " ", s_bc_max);

  /* randomly pick values for s_ab and s_bc
   * until the pair is within the Dalitz plot */
  double dalitz_bc_max = 0.0, dalitz_bc_min = 1.0;
  double s_ab = 0.0, s_bc = 0.5;
  while (s_bc > dalitz_bc_max || s_bc < dalitz_bc_min) {
    s_ab = random::uniform(s_ab_min, s_ab_max);
    s_bc = random::uniform(s_bc_min, s_bc_max);
    const double e_b_rest =
        (s_ab - mass_a * mass_a + mass_b * mass_b) / (2 * std::sqrt(s_ab));
    const double e_c_rest =
        (mass_resonance * mass_resonance - s_ab - mass_c * mass_c) /
        (2 * std::sqrt(s_ab));
    dalitz_bc_max = (e_b_rest + e_c_rest) * (e_b_rest + e_c_rest) -
                    (std::sqrt(e_b_rest * e_b_rest - mass_b * mass_b) -
                     std::sqrt(e_c_rest * e_c_rest - mass_c * mass_c)) *
                        (std::sqrt(e_b_rest * e_b_rest - mass_b * mass_b) -
                         std::sqrt(e_c_rest * e_c_rest - mass_c * mass_c));
    dalitz_bc_min = (e_b_rest + e_c_rest) * (e_b_rest + e_c_rest) -
                    (std::sqrt(e_b_rest * e_b_rest - mass_b * mass_b) +
                     std::sqrt(e_c_rest * e_c_rest - mass_c * mass_c)) *
                        (std::sqrt(e_b_rest * e_b_rest - mass_b * mass_b) +
                         std::sqrt(e_c_rest * e_c_rest - mass_c * mass_c));
  }

  log.debug("s_ab: ", s_ab, " s_bc: ", s_bc, " min: ", dalitz_bc_min,
            " max: ", dalitz_bc_max);

  // Compute energy and momentum magnitude
  const double energy_a =
      (mass_resonance * mass_resonance + mass_a * mass_a - s_bc) /
      (2 * mass_resonance);
  const double energy_c =
      (mass_resonance * mass_resonance + mass_c * mass_c - s_ab) /
      (2 * mass_resonance);
  const double energy_b =
      (s_ab + s_bc - mass_a * mass_a - mass_c * mass_c) / (2 * mass_resonance);
  const double momentum_a = std::sqrt(energy_a * energy_a - mass_a * mass_a);
  const double momentum_c = std::sqrt(energy_c * energy_c - mass_c * mass_c);
  const double momentum_b = std::sqrt(energy_b * energy_b - mass_b * mass_b);

  const double total_energy = sqrt_s();
  if (std::abs(energy_a + energy_b + energy_c - total_energy) > really_small) {
    log.warn("1->3: Ea + Eb + Ec: ", energy_a + energy_b + energy_c,
             " Total E: ", total_energy);
  }
  log.debug("Calculating the angles...");

  // momentum_a direction is random
  Angles phitheta;
  phitheta.distribute_isotropically();
  // This is the angle of the plane of the three decay particles
  outgoing_a.set_4momentum(mass_a, phitheta.threevec() * momentum_a);

  // Angle between a and b
  double theta_ab = std::acos(
      (energy_a * energy_b - 0.5 * (s_ab - mass_a * mass_a - mass_b * mass_b)) /
      (momentum_a * momentum_b));
  log.debug("theta_ab: ", theta_ab, " Ea: ", energy_a, " Eb: ", energy_b,
            " sab: ", s_ab, " pa: ", momentum_a, " pb: ", momentum_b);
  bool phi_has_changed = phitheta.add_to_theta(theta_ab);
  outgoing_b.set_4momentum(mass_b, phitheta.threevec() * momentum_b);

  // Angle between b and c
  double theta_bc = std::acos(
      (energy_b * energy_c - 0.5 * (s_bc - mass_b * mass_b - mass_c * mass_c)) /
      (momentum_b * momentum_c));
  log.debug("theta_bc: ", theta_bc, " Eb: ", energy_b, " Ec: ", energy_c,
            " sbc: ", s_bc, " pb: ", momentum_b, " pc: ", momentum_c);
  /* pass information on whether phi has changed during the last adding
   * on to add_to_theta: */
  phitheta.add_to_theta(theta_bc, phi_has_changed);
  outgoing_c.set_4momentum(mass_c, phitheta.threevec() * momentum_c);

  // Momentum check
  FourVector ptot =
      outgoing_a.momentum() + outgoing_b.momentum() + outgoing_c.momentum();

  if (std::abs(ptot.x0() - total_energy) > really_small) {
    log.warn("1->3 energy not conserved! Before: ", total_energy,
             " After: ", ptot.x0());
  }
  if (std::abs(ptot.x1()) > really_small ||
      std::abs(ptot.x2()) > really_small ||
      std::abs(ptot.x3()) > really_small) {
    log.warn("1->3 momentum check failed. Total momentum: ", ptot.threevec());
  }

  log.debug("outgoing_a: ", outgoing_a.momentum(),
            "\noutgoing_b: ", outgoing_b.momentum(),
            "\noutgoing_c: ", outgoing_c.momentum());
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

decay_channels_

DecayBranchList smash::DecayAction::decay_channels_

protected

List of possible decays.

Definition at line 97 of file decayaction.h.

total_width_

double smash::DecayAction::total_width_

protected

total decay width

Definition at line 100 of file decayaction.h.

partial_width_

double smash::DecayAction::partial_width_

protected

partial decay width to the chosen outgoing channel

Definition at line 103 of file decayaction.h.

L_

int smash::DecayAction::L_ = 0

protected

Angular momentum of the decay.

Definition at line 106 of file decayaction.h.

---

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

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