smash::Action Class Reference

#include <action.h>

Action is the base class for a generic process that takes a number of incoming particles and transforms them into any number of outgoing particles.

Currently such an action can be either a decay, a two-body collision, a wallcrossing or a thermalization. (see derived classes).

Definition at line 35 of file action.h.

Inheritance diagram for smash::Action:

Classes
class	InvalidResonanceFormation
	Thrown for example when ScatterAction is called to perform with a wrong number of final-state particles or when the energy is too low to produce the resonance. More...

Public Member Functions
	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 double	get_total_weight () const =0
	Return the total weight value, which is mainly used for the weight output entry. More...
virtual double	get_partial_weight () const =0
	Return the specific weight for the chosen outgoing channel, which is mainly used for the partial weight output entry. 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	generate_final_state ()=0
	Generate the final state for this action. More...
virtual double	perform (Particles *particles, uint32_t id_process)
	Actually perform the action, e.g. More...
bool	is_valid (const Particles &particles) const
	Check whether the action still applies. More...
bool	is_pauli_blocked (const std::vector< Particles > &ensembles, const PauliBlocker &p_bl) const
	Check if the action is Pauli-blocked. More...
const ParticleList &	incoming_particles () const
	Get the list of particles that go into the action. More...
void	update_incoming (const Particles &particles)
	Update the incoming particles that are stored in this action to the state they have in the global particle list. More...
const ParticleList &	outgoing_particles () const
	Get the list of particles that resulted from the action. More...
double	time_of_execution () const
	Get the time at which the action is supposed to be performed. More...
virtual double	check_conservation (const uint32_t id_process) const
	Check various conservation laws. More...
double	sqrt_s () const
	Determine the total energy in the center-of-mass frame [GeV]. More...
FourVector	total_momentum_of_outgoing_particles () const
	Calculate the total kinetic momentum of the outgoing particles. More...
FourVector	get_interaction_point () const
	Get the interaction point. More...
std::pair< FourVector, FourVector >	get_potential_at_interaction_point () const
	Get the skyrme and asymmetry potential at the interaction point. More...
void	set_stochastic_pos_idx ()
	Setter function that stores a random incoming particle index latter used to determine the interaction point. More...

Static Public Member Functions
static double	lambda_tilde (double a, double b, double c)
	Little helper function that calculates the lambda function (sometimes written with a tilde to better distinguish it) that appears e.g. More...
static void	sample_manybody_phasespace_impl (double sqrts, const std::vector< double > &m, std::vector< FourVector > &sampled_momenta)
	Implementation of the full n-body phase-space sampling (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...

Protected Member Functions
FourVector	total_momentum () const
	Sum of 4-momenta of incoming particles. More...
template<typename Branch >
const Branch *	choose_channel (const ProcessBranchList< Branch > &subprocesses, double total_weight)
	Decide for a particular final-state channel via Monte-Carlo and return it as a ProcessBranch. More...
virtual std::pair< double, double >	sample_masses (double kinetic_energy_cm) const
	Sample final-state masses in general X->2 processes (thus also fixing the absolute c.o.m. More...
virtual void	sample_angles (std::pair< double, double > masses, double kinetic_energy_cm)
	Sample final-state momenta in general X->2 processes (here: using an isotropical angular distribution). More...
void	sample_2body_phasespace ()
	Sample the full 2-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...
virtual void	sample_manybody_phasespace ()
	Sample the full n-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...
void	assign_formation_time_to_outgoing_particles ()
	Assign the formation time to the outgoing particles. More...
virtual void	format_debug_output (std::ostream &out) const =0
	Writes information about this action to the out stream. More...

Protected Attributes
ParticleList	incoming_particles_
	List with data of incoming particles. More...
ParticleList	outgoing_particles_
	Initially this stores only the PDG codes of final-state particles. More...
const double	time_of_execution_
	Time at which the action is supposed to be performed (absolute time in the lab frame in fm). More...
ProcessType	process_type_
	type of process More...
double	box_length_ = -1.0
	Box length: needed to determine coordinates of collision correctly in case of collision through the wall. More...
int	stochastic_position_idx_ = -1
	This stores a randomly-chosen index to an incoming particle. More...

Private Member Functions
const ParticleType &	type_of_pout (const ParticleData &p_out) const
	Get the type of a given particle. More...
const ParticleType &	type_of_pout (const ParticleTypePtr &p_out) const
	Get the particle type for given pointer to a particle type. More...

Friends
std::ostream &	operator<< (std::ostream &out, const Action &action)
	Dispatches formatting to the virtual Action::format_debug_output function. More...

Constructor & Destructor Documentation

Action() [1/4]

smash::Action::Action ( const ParticleList &  in_part, double  time  )

inline

Construct an action object with incoming particles and relative time.

Parameters

[in]	in_part	list of incoming particles
[in]	time	time at which the action is supposed to take place (relative to the current time of the incoming particles)

Definition at line 44 of file action.h.

      : incoming_particles_(in_part),
        time_of_execution_(time + in_part[0].position().x0()) {}

Action() [2/4]

smash::Action::Action ( const ParticleData &  in_part, const ParticleData &  out_part, double  time, ProcessType  type  )

inline

Construct an action object with the incoming particles, relative time, and the already known outgoing particles and type of the process.

Parameters

[in]	in_part	list of incoming particles
[in]	out_part	list of outgoing particles
[in]	time	time at which the action is supposed to take place (relative to the current time of the incoming particles)
[in]	type	type of the interaction

Definition at line 58 of file action.h.

      : incoming_particles_({in_part}),
        outgoing_particles_({out_part}),
        time_of_execution_(time + in_part.position().x0()),
        process_type_(type) {}

Action() [3/4]

smash::Action::Action ( const ParticleList &  in_part, const ParticleList &  out_part, double  absolute_execution_time, ProcessType  type  )

inline

Construct an action object with the incoming particles, absolute time, and the already known outgoing particles and type of the process.

Parameters

[in]	in_part	list of incoming particles
[in]	out_part	list of outgoing particles
[in]	absolute_execution_time	absolute time at which the action is supposed to take place
[in]	type	type of the interaction

Definition at line 75 of file action.h.

      : incoming_particles_(std::move(in_part)),
        outgoing_particles_(std::move(out_part)),
        time_of_execution_(absolute_execution_time),
        process_type_(type) {}

Action() [4/4]

smash::Action::Action ( const Action &  )

delete

Copying is disabled. Use pointers or create a new Action.

~Action()

smash::Action::~Action ( )

virtualdefault

Virtual Destructor.

Destructor.

The declaration of the destructor is necessary to make it virtual.

Member Function Documentation

operator<()

bool smash::Action::operator< ( const Action &  rhs ) const

inline

Determine whether one action takes place before another in time.

Returns

if the first argument action takes place before the other

Definition at line 96 of file action.h.

                                          {
    return time_of_execution_ < rhs.time_of_execution_;
  }

get_total_weight()

virtual double smash::Action::get_total_weight ( ) const

pure virtual

Return the total weight value, which is mainly used for the weight output entry.

It has different meanings depending of the type of action. It is the total cross section in case of a ScatterAction, the total decay width in case of a DecayAction and the shining weight in case of a DecayActionDilepton.

Prefer to use a more specific function. If there is no weight for the action type, 0 should be returned.

Returns

total cross section, decay width or shining weight

Implemented in smash::WallcrossingAction, smash::ScatterActionPhoton, smash::ScatterActionMulti, smash::ScatterAction, smash::HypersurfacecrossingAction, smash::DecayActionDilepton, smash::DecayAction, smash::BremsstrahlungAction, smash::ThermalizationAction, and smash::FreeforallAction.

get_partial_weight()

virtual double smash::Action::get_partial_weight ( ) const

pure virtual

Return the specific weight for the chosen outgoing channel, which is mainly used for the partial weight output entry.

For scatterings it will be the partial cross section, for decays (including dilepton decays) the partial decay width.

If there is no weight for the action type, 0 should be returned.

Returns

specific weight for the chosen output channel.

Implemented in smash::WallcrossingAction, smash::ScatterActionMulti, smash::ScatterAction, smash::HypersurfacecrossingAction, smash::DecayAction, smash::ThermalizationAction, and smash::FreeforallAction.

get_type()

virtual ProcessType smash::Action::get_type ( ) const

inlinevirtual

Get the process type.

Returns

type of the process

Definition at line 131 of file action.h.

{ return process_type_; }

add_process()

template<typename Branch >

void smash::Action::add_process ( ProcessBranchPtr< Branch > &  p, ProcessBranchList< Branch > &  subprocesses, double &  total_weight  )

inline

Add a new subprocess.

Parameters

[in]	p	process to be added
[out]	subprocesses	processes, where p is added to
[out]	total_weight	summed weights of all the subprocesses

Definition at line 141 of file action.h.

                                         {
    if (p->weight() > 0) {
      total_weight += p->weight();
      subprocesses.emplace_back(std::move(p));
    }
  }

add_processes()

template<typename Branch >

void smash::Action::add_processes ( ProcessBranchList< Branch >  pv, ProcessBranchList< Branch > &  subprocesses, double &  total_weight  )

inline

Add several new subprocesses at once.

Parameters

[in]	pv	processes list to be added
[out]	subprocesses	processes, where pv are added to
[out]	total_weight	summed weights of all the subprocesses

Definition at line 158 of file action.h.

                                           {
    subprocesses.reserve(subprocesses.size() + pv.size());
    for (auto &proc : pv) {
      if (proc->weight() > 0) {
        total_weight += proc->weight();
        subprocesses.emplace_back(std::move(proc));
      }
    }
  }

generate_final_state()

virtual void smash::Action::generate_final_state ( )

pure virtual

Generate the final state for this action.

This function selects a subprocess by Monte-Carlo decision and sets up the final-state particles in phase space.

Implemented in smash::WallcrossingAction, smash::ScatterActionPhoton, smash::ScatterActionMulti, smash::ScatterAction, smash::HypersurfacecrossingAction, smash::DecayAction, smash::BremsstrahlungAction, smash::ThermalizationAction, and smash::FreeforallAction.

perform()

double smash::Action::perform ( Particles *  particles, uint32_t  id_process  )

virtual

Actually perform the action, e.g.

carry out a decay or scattering by updating the particle list.

This function removes the initial-state particles from the particle list and then inserts the final-state particles. It does not do any sanity checks, but assumes that is_valid has been called to determine if the action is still valid.

Parameters

[in]	id_process	unique id of the performed process
[out]	particles	particle list that is updated

Returns

the amount of energy violated in Pythia processes (if any)

Note that you are required to increase id_process before the next call, such that you get unique numbers.

Definition at line 128 of file action.cc.

                                                                {
  assert(id_process != 0);
  double energy_violation = 0.;
  for (ParticleData &p : outgoing_particles_) {
    // store the history info
    if (process_type_ != ProcessType::Wall) {
      p.set_history(p.get_history().collisions_per_particle + 1, id_process,
                    process_type_, time_of_execution_, incoming_particles_);
    }
  }
 
  /* For elastic collisions and box wall crossings it is not necessary to remove
   * particles from the list and insert new ones, it is enough to update their
   * properties. */
  particles->update(incoming_particles_, outgoing_particles_,
                    (process_type_ != ProcessType::Elastic) &&
                        (process_type_ != ProcessType::Wall));
 
  logg[LAction].debug("Particle map now has ", particles->size(), " elements.");
 
  /* Check the conservation laws if the modifications of the total kinetic
   * energy of the outgoing particles by the mean field potentials are not
   * taken into account. */
  if (UB_lat_pointer == nullptr && UI3_lat_pointer == nullptr) {
    energy_violation = check_conservation(id_process);
  }
  return energy_violation;
}

is_valid()

bool smash::Action::is_valid

(

const Particles &

particles

)

const

Check whether the action still applies.

It can happen that a different action removed the incoming_particles from the set of existing particles in the experiment, or that the particle has scattered elastically in the meantime. In this case the Action doesn't apply anymore and should be discarded.

Parameters

[in]

particles

current particle list

Returns

true, if action still applies; false otherwise

Definition at line 29 of file action.cc.

                                                      {
  return std::all_of(
      incoming_particles_.begin(), incoming_particles_.end(),
      [&particles](const ParticleData &p) { return particles.is_valid(p); });
}

is_pauli_blocked()

bool smash::Action::is_pauli_blocked	(	const std::vector< Particles > &	ensembles,
		const PauliBlocker &	p_bl
	)		const

Check if the action is Pauli-blocked.

If there are baryons in the final state then blocking probability is \( 1 - \Pi (1-f_i) \), where the product is taken by all fermions in the final state and \( f_i \) denotes the phase-space density at the position of i-th final-state fermion.

Parameters

[in]	ensembles	current particle list, all ensembles
[in]	p_bl	PauliBlocker that stores the configurations concerning Pauli-blocking.

Returns

true, if the action is Pauli-blocked, false otherwise

Definition at line 35 of file action.cc.

                                                              {
  // Wall-crossing actions should never be blocked: currently
  // if the action is blocked, a particle continues to propagate in a straight
  // line. This would simply bring it out of the box.
  if (process_type_ == ProcessType::Wall) {
    return false;
  }
  for (const auto &p : outgoing_particles_) {
    if (p.is_baryon()) {
      const auto f =
          p_bl.phasespace_dens(p.position().threevec(), p.momentum().threevec(),
                               ensembles, p.pdgcode(), incoming_particles_);
      if (f > random::uniform(0., 1.)) {
        logg[LPauliBlocking].debug("Action ", *this,
                                   " is pauli-blocked with f = ", f);
        return true;
      }
    }
  }
  return false;
}

incoming_particles()

const ParticleList & smash::Action::incoming_particles

(

)

const

Get the list of particles that go into the action.

Returns

a list of incoming particles

Definition at line 58 of file action.cc.

                                                     {
  return incoming_particles_;
}

update_incoming()

void smash::Action::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.

Parameters

[in]

particles

current particle list

Definition at line 62 of file action.cc.

                                                       {
  for (auto &p : incoming_particles_) {
    p = particles.lookup(p);
  }
}

outgoing_particles()

const ParticleList& smash::Action::outgoing_particles ( ) const

inline

Get the list of particles that resulted from the action.

Returns

list of outgoing particles

Definition at line 247 of file action.h.

{ return outgoing_particles_; }

time_of_execution()

double smash::Action::time_of_execution ( ) const

inline

Get the time at which the action is supposed to be performed.

Returns

absolute time in the calculation frame in fm

Definition at line 254 of file action.h.

{ return time_of_execution_; }

check_conservation()

double smash::Action::check_conservation ( const uint32_t  id_process ) const

virtual

Check various conservation laws.

Parameters

[in]

id_process

process id only used for debugging output

Returns

the amount of energy conservation violated by Pythia processes (if any)

Reimplemented in smash::HypersurfacecrossingAction.

Definition at line 471 of file action.cc.

                                                                 {
  QuantumNumbers before(incoming_particles_);
  QuantumNumbers after(outgoing_particles_);
  double energy_violation = 0.;
  if (before != after) {
    std::stringstream particle_names;
    for (const auto &p : incoming_particles_) {
      particle_names << p.type().name();
    }
    particle_names << " vs. ";
    for (const auto &p : outgoing_particles_) {
      particle_names << p.type().name();
    }
    particle_names << "\n";
    std::string err_msg = before.report_deviations(after);
    /* Pythia does not conserve energy and momentum at high energy, so we just
     * print the warning and continue. */
    if ((is_string_soft_process(process_type_)) ||
        (process_type_ == ProcessType::StringHard)) {
      logg[LAction].warn() << "Conservation law violations due to Pythia\n"
                           << particle_names.str() << err_msg;
      energy_violation = after.momentum()[0] - before.momentum()[0];
      return energy_violation;
    }
    /* We allow decay of particles stable under the strong interaction to decay
     * at the end, so just warn about such a "weak" process violating
     * conservation laws */
    if (process_type_ == ProcessType::Decay &&
        incoming_particles_[0].type().is_stable()) {
      logg[LAction].warn()
          << "Conservation law violations of strong interaction in weak or "
             "e.m. decay\n"
          << particle_names.str() << err_msg;
      return energy_violation;
    }
    /* If particles are added or removed, it is not surprising that conservation
     * laws are potentially violated. Do not warn the user but print some
     * information for debug */
    if (process_type_ == ProcessType::Freeforall) {
      logg[LAction].debug()
          << "Conservation law violation, but we want it (Freeforall Action).\n"
          << particle_names.str() << err_msg;
      return energy_violation;
    }
    logg[LAction].error() << "Conservation law violations detected\n"
                          << particle_names.str() << err_msg;
    if (id_process == ID_PROCESS_PHOTON) {
      throw std::runtime_error("Conservation laws violated in photon process");
    } else {
      throw std::runtime_error("Conservation laws violated in process " +
                               std::to_string(id_process));
    }
  }
  return energy_violation;
}

sqrt_s()

double smash::Action::sqrt_s ( ) const

inline

Determine the total energy in the center-of-mass frame [GeV].

Returns

\( \sqrt{s}\) of incoming particles

Definition at line 271 of file action.h.

{ return total_momentum().abs(); }

total_momentum_of_outgoing_particles()

FourVector smash::Action::total_momentum_of_outgoing_particles

(

)

const

Calculate the total kinetic momentum of the outgoing particles.

Use this to determine the momemtum and boost of the outgoing particles by calcluating the total momentum of the incoming particles and correcting it for the effect of potentials. This function is used when the species of the outgoing particles are already determined.

Returns

total kinetic momentum of the outgoing particles [GeV]

Definition at line 157 of file action.cc.

                                                              {
  const auto potentials = get_potential_at_interaction_point();
  /* scale_B returns the difference of the total force scales of the skyrme
   * potential between the initial and final states. */
  double scale_B = 0.0;
  /* scale_I3 returns the difference of the total force scales of the symmetry
   * potential between the initial and final states. */
  double scale_I3 = 0.0;
  for (const auto &p_in : incoming_particles_) {
    // Get the force scale of the incoming particle.
    const auto scale =
        ((pot_pointer != nullptr) ? pot_pointer->force_scale(p_in.type())
                                  : std::make_pair(0.0, 0));
    scale_B += scale.first;
    scale_I3 += scale.second * p_in.type().isospin3_rel();
  }
  for (const auto &p_out : outgoing_particles_) {
    // Get the force scale of the outgoing particle.
    const auto scale = ((pot_pointer != nullptr)
                            ? pot_pointer->force_scale(type_of_pout(p_out))
                            : std::make_pair(0.0, 0));
    scale_B -= scale.first;
    scale_I3 -= scale.second * type_of_pout(p_out).isospin3_rel();
  }
  /* Rescale to get the potential difference between the
   * initial and final state, and thus get the total momentum
   * of the outgoing particles*/
  return total_momentum() + potentials.first * scale_B +
         potentials.second * scale_I3;
}

get_interaction_point()

FourVector smash::Action::get_interaction_point

(

)

const

Get the interaction point.

Returns

four vector of interaction point

Definition at line 68 of file action.cc.

                                               {
  // Estimate for the interaction point in the calculational frame
  ThreeVector interaction_point = ThreeVector(0., 0., 0.);
  std::vector<ThreeVector> propagated_positions;
  for (const auto &part : incoming_particles_) {
    ThreeVector propagated_position =
        part.position().threevec() +
        part.velocity() * (time_of_execution_ - part.position().x0());
    propagated_positions.push_back(propagated_position);
    interaction_point += propagated_position;
  }
  interaction_point /= incoming_particles_.size();
  /*
   * In case of periodic boundaries interaction point is not necessarily
   * (x1 + x2)/2. Consider only one dimension, e.g. x, the rest are analogous.
   * Instead of x, there can be x + k * L, where k is any integer and L
   * is period.Interaction point is either. Therefore, interaction point is
   * (x1 + k * L + x2 + m * L) / 2  = (x1 + x2) / 2 + n * L / 2. We need
   * this interaction point to be with [0, L], so n can be {-1, 0, 1}.
   * Which n to choose? Our guiding principle is that n should be such that
   * interaction point is closest to interacting particles.
   */
  if (box_length_ > 0 && stochastic_position_idx_ < 0) {
    assert(incoming_particles_.size() == 2);
    const ThreeVector r = propagated_positions[0] - propagated_positions[1];
    for (int i = 0; i < 3; i++) {
      const double d = std::abs(r[i]);
      if (d > 0.5 * box_length_) {
        if (interaction_point[i] >= 0.5 * box_length_) {
          interaction_point[i] -= 0.5 * box_length_;
        } else {
          interaction_point[i] += 0.5 * box_length_;
        }
      }
    }
  }
  /* In case of scatterings via the stochastic criterion, use postion of random
   * incoming particle to prevent density hotspots in grid cell centers. */
  if (stochastic_position_idx_ >= 0) {
    return incoming_particles_[stochastic_position_idx_].position();
  }
  return FourVector(time_of_execution_, interaction_point);
}

get_potential_at_interaction_point()

std::pair< FourVector, FourVector > smash::Action::get_potential_at_interaction_point

(

)

const

Get the skyrme and asymmetry potential at the interaction point.

Returns

skyrme and asymmetry potential [GeV]

Definition at line 112 of file action.cc.

          {
  const ThreeVector r = get_interaction_point().threevec();
  FourVector UB = FourVector();
  FourVector UI3 = FourVector();
  /* Check:
   * Lattice is turned on. */
  if (UB_lat_pointer != nullptr) {
    UB_lat_pointer->value_at(r, UB);
  }
  if (UI3_lat_pointer != nullptr) {
    UI3_lat_pointer->value_at(r, UI3);
  }
  return std::make_pair(UB, UI3);
}

set_stochastic_pos_idx()

void smash::Action::set_stochastic_pos_idx ( )

inline

Setter function that stores a random incoming particle index latter used to determine the interaction point.

Definition at line 303 of file action.h.

                                {
    const int max_inc_idx = incoming_particles_.size() - 1;
    stochastic_position_idx_ = random::uniform_int(0, max_inc_idx);
  }

lambda_tilde()

static double smash::Action::lambda_tilde ( double  a, double  b, double  c  )

inlinestatic

Little helper function that calculates the lambda function (sometimes written with a tilde to better distinguish it) that appears e.g.

in the relative velocity or 3-to-2 probability calculation, where it is used with a=s, b=m1^2 and c=m2^2. Defintion found e.g. in Seifert:2017oyb [50], eq. (5).

Definition at line 315 of file action.h.

                                                           {
    const double res = (a - b - c) * (a - b - c) - 4. * b * c;
    if (res < 0.0) {
      // floating point precision problem
      return 0.0;
    }
    return res;
  }

sample_manybody_phasespace_impl()

void smash::Action::sample_manybody_phasespace_impl ( double  sqrts, const std::vector< double > &  m, std::vector< FourVector > &  sampled_momenta  )

static

Implementation of the full n-body phase-space sampling (masses, momenta, angles) in the center-of-mass frame for the final state particles.

Function is static for convenient testing.

Using the M-method from CERN-68-15 report, paragraph 9.6 1) Generate invariant masses M12, M123, M1234, etc from distribution dM12 x dM123 x dM1234 x ... This is not trivial because of the integration limits. Here the idea is to change variables to T12 = M12 - (m1 + m2), T123 = M123 - (m1 + m2 + m3), etc. Then we need to generate uniform T such that 0 <= T12 <= T123 <= T1234 <= ... <= sqrts - sum (m_i). For the latter there is a trick: generate values uniformly in [0, sqrts - sum (m_i)] and then sort the values. 2) accept or reject this combination of invariant masses with weight proportional to R2(sqrt, M_{n-1}, m_n) x R2(M_{n-1}, M_{n-2}, m_{n-1}) x ... x R2(M2, m1, m2) x (prod M_i). Maximum weight is estmated heuristically, here I'm using an idea by Scott Pratt that maximum is close to T12 = T123 = T1234 = ... = (sqrts - sum (m_i)) / (n - 1)

Definition at line 313 of file action.cc.

                                            {
  const size_t n = m.size();
  assert(n > 1);
  sampled_momenta.resize(n);
 
  // Arrange a convenient vector of m1, m1 + m2, m1 + m2 + m3, ...
  std::vector<double> msum(n);
  std::partial_sum(m.begin(), m.end(), msum.begin());
  const double msum_all = msum[n - 1];
  int rejection_counter = -1;
  if (sqrts <= msum_all) {
    logg[LAction].error() << "sample_manybody_phasespace_impl: "
                          << "Can't sample when sqrts = " << sqrts
                          << " < msum = " << msum_all;
  }
 
  double w, r01;
  std::vector<double> Minv(n);
 
  double weight_sqr_max = 1;
  const double Ekin_share = (sqrts - msum_all) / (n - 1);
  for (size_t i = 1; i < n; i++) {
    // This maximum estimate idea is due Scott Pratt: maximum should be
    // roughly at equal kinetic energies
    weight_sqr_max *= pCM_sqr(i * Ekin_share + msum[i],
                              (i - 1) * Ekin_share + msum[i - 1], m[i]);
  }
  // Maximum estimate is rough and can be wrong. We multiply it by additional
  // factor to be on the safer side.
  const double safety_factor = 1.1 + (n - 2) * 0.2;
  weight_sqr_max *= (safety_factor * safety_factor);
  bool first_warning = true;
 
  do {
    // Generate invariant masses of 1, 12, 123, 1243, etc.
    // Minv = {m1, M12, M123, ..., M123n-1, sqrts}
    Minv[0] = 0.0;
    Minv[n - 1] = sqrts - msum_all;
    for (size_t i = 1; i < n - 1; i++) {
      Minv[i] = random::uniform(0.0, sqrts - msum_all);
    }
    std::sort(Minv.begin(), Minv.end());
    for (size_t i = 0; i < n; i++) {
      Minv[i] += msum[i];
    }
 
    double weight_sqr = 1;
    for (size_t i = 1; i < n; i++) {
      weight_sqr *= pCM_sqr(Minv[i], Minv[i - 1], m[i]);
    }
 
    rejection_counter++;
    r01 = random::canonical();
    w = weight_sqr / weight_sqr_max;
    if (w > 1.0) {
      logg[LAction].warn()
          << "sample_manybody_phasespace_impl: alarm, weight > 1, w^2 = " << w
          << ". Increase safety factor." << std::endl;
    }
    if (rejection_counter > 20 && first_warning) {
      logg[LAction].warn() << "sample_manybody_phasespace_impl: "
                           << "likely hanging, way too many rejections,"
                           << " n = " << n << ", sqrts = " << sqrts
                           << ", msum = " << msum_all;
      first_warning = false;
    }
  } while (w < r01 * r01);
 
  // Boost particles to the right frame
  std::vector<ThreeVector> beta(n);
  for (size_t i = n - 1; i > 0; i--) {
    const double pcm = pCM(Minv[i], Minv[i - 1], m[i]);
    Angles phitheta;
    phitheta.distribute_isotropically();
    const ThreeVector isotropic_unitvector = phitheta.threevec();
    sampled_momenta[i] = FourVector(std::sqrt(m[i] * m[i] + pcm * pcm),
                                    pcm * isotropic_unitvector);
    if (i >= 2) {
      beta[i - 2] = pcm * isotropic_unitvector /
                    std::sqrt(pcm * pcm + Minv[i - 1] * Minv[i - 1]);
    }
    if (i == 1) {
      sampled_momenta[0] = FourVector(std::sqrt(m[0] * m[0] + pcm * pcm),
                                      -pcm * isotropic_unitvector);
    }
  }
 
  for (size_t i = 0; i < n - 2; i++) {
    // After each boost except the last one the sum of 3-momenta should be 0
    FourVector ptot = FourVector(0.0, 0.0, 0.0, 0.0);
    for (size_t j = 0; j <= i + 1; j++) {
      ptot += sampled_momenta[j];
    }
    logg[LAction].debug() << "Total momentum of 0.." << i + 1 << " = "
                          << ptot.threevec() << " and should be (0, 0, 0). "
                          << std::endl;
 
    // Boost the first i+1 particles to the next CM frame
    for (size_t j = 0; j <= i + 1; j++) {
      sampled_momenta[j] = sampled_momenta[j].lorentz_boost(beta[i]);
    }
  }
 
  FourVector ptot_all = FourVector(0.0, 0.0, 0.0, 0.0);
  for (size_t j = 0; j < n; j++) {
    ptot_all += sampled_momenta[j];
  }
  logg[LAction].debug() << "Total 4-momentum = " << ptot_all << ", should be ("
                        << sqrts << ", 0, 0, 0)" << std::endl;
}

total_momentum()

FourVector smash::Action::total_momentum ( ) const

inlineprotected

Sum of 4-momenta of incoming particles.

Definition at line 379 of file action.h.

                                    {
    FourVector mom(0.0, 0.0, 0.0, 0.0);
    for (const auto &p : incoming_particles_) {
      mom += p.momentum();
    }
    return mom;
  }

choose_channel()

template<typename Branch >

const Branch* smash::Action::choose_channel ( const ProcessBranchList< Branch > &  subprocesses, double  total_weight  )

inlineprotected

Decide for a particular final-state channel via Monte-Carlo and return it as a ProcessBranch.

Template Parameters

Branch

Type of processbranch

Parameters

[in]	subprocesses	list of possible processes
[in]	total_weight	summed weight of all processes

Returns

ProcessBranch that is sampled

Definition at line 397 of file action.h.

                                                    {
    double random_weight = random::uniform(0., total_weight);
    double weight_sum = 0.;
    /* Loop through all subprocesses and select one by Monte Carlo, based on
     * their weights.  */
    for (const auto &proc : subprocesses) {
      weight_sum += proc->weight();
      if (random_weight <= weight_sum) {
        /* Return the full process information. */
        return proc.get();
      }
    }
    /* Should never get here. */
    logg[LAction].fatal(SMASH_SOURCE_LOCATION,
                        "Problem in choose_channel: ", subprocesses.size(), " ",
                        weight_sum, " ", total_weight, " ",
                        //          random_weight, "\n", *this);
                        random_weight, "\n");
    std::abort();
  }

sample_masses()

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

protectedvirtual

Sample final-state masses in general X->2 processes (thus also fixing the absolute c.o.m.

momentum).

Parameters

[in]

kinetic_energy_cm

total kinetic energy of the outgoing particles in their center of mass frame [GeV]

Exceptions

InvalidResonanceFormation

Returns

masses of final state particles

Reimplemented in smash::DecayAction.

Definition at line 250 of file action.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() +
                                 incoming_particles_[1].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);
  } else if (!t_b.is_stable() && t_a.is_stable()) {
    masses.second = t_b.sample_resonance_mass(t_a.mass(), kinetic_energy_cm);
  } 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);
  }
  return masses;
}

sample_angles()

void smash::Action::sample_angles ( std::pair< double, double >  masses, double  kinetic_energy_cm  )

protectedvirtual

Sample final-state momenta in general X->2 processes (here: using an isotropical angular distribution).

Parameters

[in]	kinetic_energy_cm	total kinetic energy of the outgoing particles in their center of mass frame [GeV]
[in]	masses	masses of each of the final state particles

Reimplemented in smash::ScatterAction.

Definition at line 279 of file action.cc.

                                                           {
  ParticleData *p_a = &outgoing_particles_[0];
  ParticleData *p_b = &outgoing_particles_[1];
 
  const double pcm = pCM(kinetic_energy_cm, masses.first, masses.second);
  if (!(pcm > 0.0)) {
    logg[LAction].warn("Particle: ", p_a->pdgcode(), " radial momentum: ", pcm);
    logg[LAction].warn("Ektot: ", kinetic_energy_cm, " m_a: ", masses.first,
                       " m_b: ", masses.second);
  }
  /* Here we assume an isotropic angular distribution. */
  Angles phitheta;
  phitheta.distribute_isotropically();
 
  p_a->set_4momentum(masses.first, phitheta.threevec() * pcm);
  p_b->set_4momentum(masses.second, -phitheta.threevec() * pcm);
  /* Debug message is printed before boost, so that p_a and p_b are
   * the momenta in the center of mass frame and thus opposite to
   * each other.*/
  logg[LAction].debug("p_a: ", *p_a, "\np_b: ", *p_b);
}

sample_2body_phasespace()

void smash::Action::sample_2body_phasespace ( )

protected

Sample the full 2-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles.

Definition at line 302 of file action.cc.

                                     {
  /* This function only operates on 2-particle final states. */
  assert(outgoing_particles_.size() == 2);
  const FourVector p_tot = total_momentum_of_outgoing_particles();
  const double cm_kin_energy = p_tot.abs();
  // first sample the masses
  const std::pair<double, double> masses = sample_masses(cm_kin_energy);
  // after the masses are fixed (and thus also pcm), sample the angles
  sample_angles(masses, cm_kin_energy);
}

sample_manybody_phasespace()

void smash::Action::sample_manybody_phasespace ( )

protectedvirtual

Sample the full n-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles.

Exceptions

std::invalid_argument

if one outgoing particle is a resonance

Reimplemented in smash::DecayActionDilepton.

Definition at line 442 of file action.cc.

                                        {
  const size_t n = outgoing_particles_.size();
  if (n < 3) {
    throw std::invalid_argument(
        "sample_manybody_phasespace: number of outgoing particles should be 3 "
        "or more");
  }
  bool all_stable = true;
  for (size_t i = 0; i < n; i++) {
    all_stable = all_stable && outgoing_particles_[i].type().is_stable();
  }
  if (!all_stable) {
    throw std::invalid_argument(
        "sample_manybody_phasespace: Found resonance in to be sampled outgoing "
        "particles, but assumes stable particles.");
  }
 
  std::vector<double> m(n);
  for (size_t i = 0; i < n; i++) {
    m[i] = outgoing_particles_[i].type().mass();
  }
  std::vector<FourVector> p(n);
 
  sample_manybody_phasespace_impl(sqrt_s(), m, p);
  for (size_t i = 0; i < n; i++) {
    outgoing_particles_[i].set_4momentum(p[i]);
  }
}

assign_formation_time_to_outgoing_particles()

void smash::Action::assign_formation_time_to_outgoing_particles ( )

protected

Assign the formation time to the outgoing particles.

The formation time is set to the largest formation time of the incoming particles, if it is larger than the execution time. The newly produced particles are supposed to continue forming exactly like the latest forming ingoing particle. Therefore the details on the formation are adopted. The initial cross section scaling factor of the incoming particles is considered to also be the scaling factor of the newly produced outgoing particles. If the formation time is smaller than the exectution time, the execution time is taken to be the formation time.

Note: Make sure to assign the formation times before boosting the outgoing particles to the computational frame.

Definition at line 188 of file action.cc.

                                                         {
  /* Find incoming particle with largest formation time i.e. the last formed
   * incoming particle. If all particles form at the same time, take the one
   * with the lowest cross section scaling factor */
  ParticleList::iterator last_formed_in_part;
  bool all_incoming_same_formation_time =
      std::all_of(incoming_particles_.begin() + 1, incoming_particles_.end(),
                  [&](const ParticleData &data_comp) {
                    return std::abs(incoming_particles_[0].formation_time() -
                                    data_comp.formation_time()) < really_small;
                  });
  if (all_incoming_same_formation_time) {
    last_formed_in_part =
        std::min_element(incoming_particles_.begin(), incoming_particles_.end(),
                         [](const ParticleData &a, const ParticleData &b) {
                           return a.initial_xsec_scaling_factor() <
                                  b.initial_xsec_scaling_factor();
                         });
  } else {
    last_formed_in_part =
        std::max_element(incoming_particles_.begin(), incoming_particles_.end(),
                         [](const ParticleData &a, const ParticleData &b) {
                           return a.formation_time() < b.formation_time();
                         });
  }
 
  const double form_time_begin = last_formed_in_part->begin_formation_time();
  const double sc = last_formed_in_part->initial_xsec_scaling_factor();
 
  if (last_formed_in_part->formation_time() > time_of_execution_) {
    for (ParticleData &new_particle : outgoing_particles_) {
      if (new_particle.initial_xsec_scaling_factor() < 1.0) {
        /* The new cross section scaling factor will be the product of the
         * cross section scaling factor of the ingoing particles and of the
         * outgoing ones (since the outgoing ones are also string fragments
         * and thus take time to form). */
        double sc_out = new_particle.initial_xsec_scaling_factor();
        new_particle.set_cross_section_scaling_factor(sc * sc_out);
        if (last_formed_in_part->formation_time() >
            new_particle.formation_time()) {
          /* If the unformed incoming particles' formation time is larger than
           * the current outgoing particle's formation time, then the latter
           * is overwritten by the former*/
          new_particle.set_slow_formation_times(
              time_of_execution_, last_formed_in_part->formation_time());
        }
      } else {
        // not a string product
        new_particle.set_slow_formation_times(
            form_time_begin, last_formed_in_part->formation_time());
        new_particle.set_cross_section_scaling_factor(sc);
      }
    }
  } else {
    for (ParticleData &new_particle : outgoing_particles_) {
      if (new_particle.initial_xsec_scaling_factor() == 1.0) {
        new_particle.set_formation_time(time_of_execution_);
      }
    }
  }
}

type_of_pout() [1/2]

const ParticleType& smash::Action::type_of_pout ( const ParticleData &  p_out ) const

inlineprivate

Get the type of a given particle.

Parameters

[in]

p_out

particle of which the type will be returned

Returns

type of given particle

Definition at line 499 of file action.h.

                                                                    {
    return p_out.type();
  }

type_of_pout() [2/2]

const ParticleType& smash::Action::type_of_pout ( const ParticleTypePtr &  p_out ) const

inlineprivate

Get the particle type for given pointer to a particle type.

Helper function for total_momentum_of_outgoing_particles

Parameters

[in]

p_out

pointer to a particle type

Returns

particle type

Definition at line 510 of file action.h.

                                                                       {
    return *p_out;
  }

Member Data Documentation

incoming_particles_

ParticleList smash::Action::incoming_particles_

protected

List with data of incoming particles.

Definition at line 345 of file action.h.

outgoing_particles_

ParticleList smash::Action::outgoing_particles_

protected

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

After perform was called it contains the complete particle data of the outgoing particles.

Definition at line 353 of file action.h.

time_of_execution_

const double smash::Action::time_of_execution_

protected

Time at which the action is supposed to be performed (absolute time in the lab frame in fm).

Definition at line 359 of file action.h.

process_type_

ProcessType smash::Action::process_type_

protected

type of process

Definition at line 362 of file action.h.

box_length_

double smash::Action::box_length_ = -1.0

protected

Box length: needed to determine coordinates of collision correctly in case of collision through the wall.

Ignored if negative.

Definition at line 369 of file action.h.

stochastic_position_idx_

int smash::Action::stochastic_position_idx_ = -1

protected

This stores a randomly-chosen index to an incoming particle.

If non-negative, the the interaction point equals the postion of the chosen particle (index). This is done for the stochastic criterion.

Definition at line 376 of file action.h.

---

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

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