Version: SMASH-1.8
smash::Action Class Referenceabstract

#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:
[legend]
Collaboration diagram for smash::Action:
[legend]

Classes

class  InvalidResonanceFormation
 

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 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, FourVectorget_potential_at_interaction_point () const
 Get the skyrme and asymmetry potential at the interaction point. 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_3body_phasespace ()
 Sample the full 3-body phase-space (masses, momenta, angles) in the center-of-mass frame for the final state particles. More...
 
virtual void format_debug_output (std::ostream &out) const =0
 

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/c). More...
 
ProcessType process_type_
 type of process More...
 

Private Member Functions

const ParticleTypetype_of_pout (const ParticleData &p_out) const
 Get the type of a given particle. More...
 
const ParticleTypetype_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)
 

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_partlist of incoming particles
[in]timetime 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.

45  : incoming_particles_(in_part),
46  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_partlist of incoming particles
[in]out_partlist of outgoing particles
[in]timetime at which the action is supposed to take place (relative to the current time of the incoming particles)
[in]typetype of the interaction

Definition at line 58 of file action.h.

60  : incoming_particles_({in_part}),
61  outgoing_particles_({out_part}),
62  time_of_execution_(time + in_part.position().x0()),
63  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_partlist of incoming particles
[in]out_partlist of outgoing particles
[in]absolute_execution_timeabsolute time at which the action is supposed to take place
[in]typetype of the interaction

Definition at line 75 of file action.h.

77  : incoming_particles_(std::move(in_part)),
78  outgoing_particles_(std::move(out_part)),
79  time_of_execution_(absolute_execution_time),
80  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.

96  {
97  return time_of_execution_ < rhs.time_of_execution_;
98  }

◆ 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::ScatterAction, smash::BremsstrahlungAction, smash::DecayAction, smash::ScatterActionPhoton, smash::DecayActionDilepton, smash::HypersurfacecrossingAction, smash::WallcrossingAction, and smash::ThermalizationAction.

Here is the caller graph for this function:

◆ 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::ScatterAction, smash::DecayAction, smash::HypersurfacecrossingAction, smash::WallcrossingAction, and smash::ThermalizationAction.

Here is the caller graph for this function:

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

131 { return process_type_; }
Here is the caller graph for this function:

◆ 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]pprocess to be added
[out]subprocessesprocesses, where p is added to
[out]total_weightsummed weights of all the subprocesses

Definition at line 141 of file action.h.

143  {
144  if (p->weight() > 0) {
145  total_weight += p->weight();
146  subprocesses.emplace_back(std::move(p));
147  }
148  }

◆ 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]pvprocesses list to be added
[out]subprocessesprocesses, where pv are added to
[out]total_weightsummed weights of all the subprocesses

Definition at line 158 of file action.h.

160  {
161  subprocesses.reserve(subprocesses.size() + pv.size());
162  for (auto &proc : pv) {
163  if (proc->weight() > 0) {
164  total_weight += proc->weight();
165  subprocesses.emplace_back(std::move(proc));
166  }
167  }
168  }

◆ 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::ScatterAction, smash::ScatterActionPhoton, smash::DecayAction, smash::BremsstrahlungAction, smash::HypersurfacecrossingAction, smash::WallcrossingAction, and smash::ThermalizationAction.

◆ perform()

void 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_processunique id of the performed process
[out]particlesparticle list that is updated

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

Definition at line 95 of file action.cc.

95  {
96  assert(id_process != 0);
97 
98  for (ParticleData &p : outgoing_particles_) {
99  // store the history info
101  p.set_history(p.get_history().collisions_per_particle + 1, id_process,
103  }
104  }
105 
106  /* For elastic collisions and box wall crossings it is not necessary to remove
107  * particles from the list and insert new ones, it is enough to update their
108  * properties. */
109  particles->update(incoming_particles_, outgoing_particles_,
112 
113  logg[LAction].debug("Particle map now has ", particles->size(), " elements.");
114 
115  /* Check the conservation laws if the modifications of the total kinetic
116  * energy of the outgoing particles by the mean field potentials are not
117  * taken into account. */
118  if (UB_lat_pointer == nullptr && UI3_lat_pointer == nullptr) {
119  check_conservation(id_process);
120  }
121 }
Here is the call graph for this function:

◆ 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]particlescurrent particle list
Returns
true, if action still applies; false otherwise

Definition at line 30 of file action.cc.

30  {
31  return std::all_of(
33  [&particles](const ParticleData &p) { return particles.is_valid(p); });
34 }
Here is the call graph for this function:

◆ is_pauli_blocked()

bool smash::Action::is_pauli_blocked ( const Particles particles,
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]particlescurrent particle list
[in]p_blPauliBlocker that stores the configurations concerning Pauli-blocking.
Returns
true, if the action is Pauli-blocked, false otherwise

Definition at line 36 of file action.cc.

37  {
38  // Wall-crossing actions should never be blocked: currently
39  // if the action is blocked, a particle continues to propagate in a straight
40  // line. This would simply bring it out of the box.
42  return false;
43  }
44  for (const auto &p : outgoing_particles_) {
45  if (p.is_baryon()) {
46  const auto f =
47  p_bl.phasespace_dens(p.position().threevec(), p.momentum().threevec(),
48  particles, p.pdgcode(), incoming_particles_);
49  if (f > random::uniform(0., 1.)) {
50  logg[LPauliBlocking].debug("Action ", *this,
51  " is pauli-blocked with f = ", f);
52  return true;
53  }
54  }
55  }
56  return false;
57 }
Here is the call graph for this function:

◆ 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 59 of file action.cc.

59  {
60  return incoming_particles_;
61 }
Here is the caller graph for this function:

◆ 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]particlescurrent particle list

Definition at line 63 of file action.cc.

63  {
64  for (auto &p : incoming_particles_) {
65  p = particles.lookup(p);
66  }
67 }
Here is the call graph for this function:

◆ 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 245 of file action.h.

245 { return outgoing_particles_; }
Here is the caller graph for this function:

◆ 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/c

Definition at line 252 of file action.h.

252 { return time_of_execution_; }

◆ check_conservation()

void smash::Action::check_conservation ( const uint32_t  id_process) const
virtual

Check various conservation laws.

Parameters
[in]id_processprocess id only used for debugging output

Reimplemented in smash::HypersurfacecrossingAction.

Definition at line 246 of file action.cc.

246  {
247  QuantumNumbers before(incoming_particles_);
248  QuantumNumbers after(outgoing_particles_);
249  if (before != after) {
250  std::stringstream particle_names;
251  for (const auto &p : incoming_particles_) {
252  particle_names << p.type().name();
253  }
254  particle_names << " vs. ";
255  for (const auto &p : outgoing_particles_) {
256  particle_names << p.type().name();
257  }
258  particle_names << "\n";
259  std::string err_msg = before.report_deviations(after);
260  logg[LAction].error() << particle_names.str() << err_msg;
261  /* Pythia does not conserve energy and momentum at high energy, so we just
262  * print the error and continue. */
265  return;
266  }
267  if (id_process == ID_PROCESS_PHOTON) {
268  throw std::runtime_error("Conservation laws violated in photon process");
269  } else {
270  throw std::runtime_error("Conservation laws violated in process " +
271  std::to_string(id_process));
272  }
273  }
274 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 266 of file action.h.

266 { return total_momentum().abs(); }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ total_momentum_of_outgoing_particles()

FourVector smash::Action::total_momentum_of_outgoing_particles ( ) const

Calculate the total kinetic momentum of the outgoing particles.

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 123 of file action.cc.

123  {
124  const auto potentials = get_potential_at_interaction_point();
125  /* scale_B returns the difference of the total force scales of the skyrme
126  * potential between the initial and final states. */
127  double scale_B = 0.0;
128  /* scale_I3 returns the difference of the total force scales of the symmetry
129  * potential between the initial and final states. */
130  double scale_I3 = 0.0;
131  for (const auto &p_in : incoming_particles_) {
132  // Get the force scale of the incoming particle.
133  const auto scale =
134  ((pot_pointer != nullptr) ? pot_pointer->force_scale(p_in.type())
135  : std::make_pair(0.0, 0));
136  scale_B += scale.first;
137  scale_I3 += scale.second * p_in.type().isospin3_rel();
138  }
139  for (const auto &p_out : outgoing_particles_) {
140  // Get the force scale of the outgoing particle.
141  const auto scale = ((pot_pointer != nullptr)
143  : std::make_pair(0.0, 0));
144  scale_B -= scale.first;
145  scale_I3 -= scale.second * type_of_pout(p_out).isospin3_rel();
146  }
147  /* Rescale to get the potential difference between the
148  * initial and final state, and thus get the total momentum
149  * of the outgoing particles*/
150  return total_momentum() + potentials.first * scale_B +
151  potentials.second * scale_I3;
152 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ get_interaction_point()

FourVector smash::Action::get_interaction_point ( ) const

Get the interaction point.

Returns
four vector of interaction point

Definition at line 69 of file action.cc.

69  {
70  // Estimate for the interaction point in the calculational frame
71  FourVector interaction_point = FourVector(0., 0., 0., 0.);
72  for (const auto &part : incoming_particles_) {
73  interaction_point += part.position();
74  }
75  interaction_point /= incoming_particles_.size();
76  return interaction_point;
77 }
Here is the caller graph for this function:

◆ 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 79 of file action.cc.

80  {
81  const ThreeVector r = get_interaction_point().threevec();
82  FourVector UB = FourVector();
83  FourVector UI3 = FourVector();
84  /* Check:
85  * Lattice is turned on. */
86  if (UB_lat_pointer != nullptr) {
87  UB_lat_pointer->value_at(r, UB);
88  }
89  if (UI3_lat_pointer != nullptr) {
90  UI3_lat_pointer->value_at(r, UI3);
91  }
92  return std::make_pair(UB, UI3);
93 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ total_momentum()

FourVector smash::Action::total_momentum ( ) const
inlineprotected

Sum of 4-momenta of incoming particles.

Definition at line 324 of file action.h.

324  {
325  FourVector mom(0.0, 0.0, 0.0, 0.0);
326  for (const auto &p : incoming_particles_) {
327  mom += p.momentum();
328  }
329  return mom;
330  }
Here is the caller graph for this function:

◆ 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
BranchType of processbranch
Parameters
[in]subprocesseslist of possible processes
[in]total_weightsummed weight of all processes
Returns
ProcessBranch that is sampled

Definition at line 342 of file action.h.

343  {
344  double random_weight = random::uniform(0., total_weight);
345  double weight_sum = 0.;
346  /* Loop through all subprocesses and select one by Monte Carlo, based on
347  * their weights. */
348  for (const auto &proc : subprocesses) {
349  weight_sum += proc->weight();
350  if (random_weight <= weight_sum) {
351  /* Return the full process information. */
352  return proc.get();
353  }
354  }
355  /* Should never get here. */
357  "Problem in choose_channel: ", subprocesses.size(), " ",
358  weight_sum, " ", total_weight, " ",
359  // random_weight, "\n", *this);
360  random_weight, "\n");
361  std::abort();
362  }
Here is the call graph for this function:

◆ 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_cmtotal 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 154 of file action.cc.

155  {
156  const ParticleType &t_a = outgoing_particles_[0].type();
157  const ParticleType &t_b = outgoing_particles_[1].type();
158  // start with pole masses
159  std::pair<double, double> masses = {t_a.mass(), t_b.mass()};
160 
161  if (kinetic_energy_cm < t_a.min_mass_kinematic() + t_b.min_mass_kinematic()) {
162  const std::string reaction = incoming_particles_[0].type().name() +
163  incoming_particles_[1].type().name() + "→" +
164  t_a.name() + t_b.name();
165  throw InvalidResonanceFormation(
166  reaction + ": not enough energy, " + std::to_string(kinetic_energy_cm) +
167  " < " + std::to_string(t_a.min_mass_kinematic()) + " + " +
168  std::to_string(t_b.min_mass_kinematic()));
169  }
170 
171  /* If one of the particles is a resonance, sample its mass. */
172  if (!t_a.is_stable() && t_b.is_stable()) {
173  masses.first = t_a.sample_resonance_mass(t_b.mass(), kinetic_energy_cm);
174  } else if (!t_b.is_stable() && t_a.is_stable()) {
175  masses.second = t_b.sample_resonance_mass(t_a.mass(), kinetic_energy_cm);
176  } else if (!t_a.is_stable() && !t_b.is_stable()) {
177  // two resonances in final state
178  masses = t_a.sample_resonance_masses(t_b, kinetic_energy_cm);
179  }
180  return masses;
181 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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_cmtotal kinetic energy of the outgoing particles in their center of mass frame [GeV]
[in]massesmasses of each of the final state particles

Reimplemented in smash::ScatterAction.

Definition at line 183 of file action.cc.

184  {
185  ParticleData *p_a = &outgoing_particles_[0];
186  ParticleData *p_b = &outgoing_particles_[1];
187 
188  const double pcm = pCM(kinetic_energy_cm, masses.first, masses.second);
189  if (!(pcm > 0.0)) {
190  logg[LAction].warn("Particle: ", p_a->pdgcode(), " radial momentum: ", pcm);
191  logg[LAction].warn("Ektot: ", kinetic_energy_cm, " m_a: ", masses.first,
192  " m_b: ", masses.second);
193  }
194  /* Here we assume an isotropic angular distribution. */
195  Angles phitheta;
196  phitheta.distribute_isotropically();
197 
198  p_a->set_4momentum(masses.first, phitheta.threevec() * pcm);
199  p_b->set_4momentum(masses.second, -phitheta.threevec() * pcm);
200  /* Debug message is printed before boost, so that p_a and p_b are
201  * the momenta in the center of mass frame and thus opposite to
202  * each other.*/
203  logg[LAction].debug("p_a: ", *p_a, "\np_b: ", *p_b);
204 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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 206 of file action.cc.

206  {
207  /* This function only operates on 2-particle final states. */
208  assert(outgoing_particles_.size() == 2);
209  const FourVector p_tot = total_momentum_of_outgoing_particles();
210  const double cm_kin_energy = p_tot.abs();
211  // first sample the masses
212  const std::pair<double, double> masses = sample_masses(cm_kin_energy);
213  // after the masses are fixed (and thus also pcm), sample the angles
214  sample_angles(masses, cm_kin_energy);
215 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ sample_3body_phasespace()

void smash::Action::sample_3body_phasespace ( )
protectedvirtual

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

Reimplemented in smash::BremsstrahlungAction, and smash::DecayActionDilepton.

Definition at line 217 of file action.cc.

217  {
218  assert(outgoing_particles_.size() == 3);
219  const double m_a = outgoing_particles_[0].type().mass(),
220  m_b = outgoing_particles_[1].type().mass(),
221  m_c = outgoing_particles_[2].type().mass();
222  const double sqrts = sqrt_s();
223 
224  // sample mab from pCM(sqrt, mab, mc) pCM (mab, ma, mb) <= sqrts^2/4
225  double mab, r, probability, pcm_ab, pcm;
226  do {
227  mab = random::uniform(m_a + m_b, sqrts - m_c);
228  r = random::canonical();
229  pcm = pCM(sqrts, mab, m_c);
230  pcm_ab = pCM(mab, m_a, m_b);
231  probability = pcm * pcm_ab * 4 / (sqrts * sqrts);
232  } while (r > probability);
233  Angles phitheta;
234  phitheta.distribute_isotropically();
235  outgoing_particles_[2].set_4momentum(m_c, pcm * phitheta.threevec());
236  const ThreeVector beta_cm =
237  pcm * phitheta.threevec() / std::sqrt(pcm * pcm + mab * mab);
238 
239  phitheta.distribute_isotropically();
240  outgoing_particles_[0].set_4momentum(m_a, pcm_ab * phitheta.threevec());
241  outgoing_particles_[1].set_4momentum(m_b, -pcm_ab * phitheta.threevec());
242  outgoing_particles_[0].boost_momentum(beta_cm);
243  outgoing_particles_[1].boost_momentum(beta_cm);
244 }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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_outparticle of which the type will be returned
Returns
type of given particle

Definition at line 425 of file action.h.

425  {
426  return p_out.type();
427  }
Here is the call graph for this function:
Here is the caller graph for this function:

◆ 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_outpointer to a particle type
Returns
particle type

Definition at line 436 of file action.h.

436  {
437  return *p_out;
438  }

Member Data Documentation

◆ incoming_particles_

ParticleList smash::Action::incoming_particles_
protected

List with data of incoming particles.

Definition at line 304 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 312 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/c).

Definition at line 318 of file action.h.

◆ process_type_

ProcessType smash::Action::process_type_
protected

type of process

Definition at line 321 of file action.h.


The documentation for this class was generated from the following files:
smash::Action::sample_angles
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...
Definition: action.cc:183
smash::Action::incoming_particles_
ParticleList incoming_particles_
List with data of incoming particles.
Definition: action.h:304
smash::ProcessType::StringHard
hard string process involving 2->2 QCD process by PYTHIA.
smash::Action::total_momentum
FourVector total_momentum() const
Sum of 4-momenta of incoming particles.
Definition: action.h:324
smash::Action::time_of_execution_
const double time_of_execution_
Time at which the action is supposed to be performed (absolute time in the lab frame in fm/c).
Definition: action.h:318
smash::LPauliBlocking
static constexpr int LPauliBlocking
Definition: action.cc:28
smash::Potentials::force_scale
static std::pair< double, int > force_scale(const ParticleType &data)
Evaluates the scaling factor of the forces acting on the particles.
Definition: potentials.cc:164
smash::Action::check_conservation
virtual void check_conservation(const uint32_t id_process) const
Check various conservation laws.
Definition: action.cc:246
smash::is_string_soft_process
bool is_string_soft_process(ProcessType p)
Check if a given process type is a soft string excitation.
Definition: processbranch.cc:18
smash::logg
std::array< einhard::Logger<>, std::tuple_size< LogArea::AreaTuple >::value > logg
An array that stores all pre-configured Logger objects.
Definition: logging.cc:39
smash::all_of
bool all_of(Container &&c, UnaryPredicate &&p)
Convenience wrapper for std::all_of that operates on a complete container.
Definition: algorithms.h:80
smash::pCM
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept
Definition: kinematics.h:79
smash::UI3_lat_pointer
RectangularLattice< FourVector > * UI3_lat_pointer
Pointer to the symmmetry potential on the lattice.
Definition: potential_globals.cc:16
source_location
#define source_location
Hackery that is required to output the location in the source code where the log statement occurs.
Definition: logging.h:240
smash::ProcessType::Wall
box wall crossing
smash::UB_lat_pointer
RectangularLattice< FourVector > * UB_lat_pointer
Pointer to the skyrme potential on the lattice.
Definition: potential_globals.cc:15
smash::Action::type_of_pout
const ParticleType & type_of_pout(const ParticleData &p_out) const
Get the type of a given particle.
Definition: action.h:425
smash::ParticleType::isospin3_rel
double isospin3_rel() const
Definition: particletype.h:179
smash::Action::process_type_
ProcessType process_type_
type of process
Definition: action.h:321
smash::LAction
static constexpr int LAction
Definition: action.h:25
smash::Action::get_potential_at_interaction_point
std::pair< FourVector, FourVector > get_potential_at_interaction_point() const
Get the skyrme and asymmetry potential at the interaction point.
Definition: action.cc:79
smash::Action::total_momentum_of_outgoing_particles
FourVector total_momentum_of_outgoing_particles() const
Calculate the total kinetic momentum of the outgoing particles.
Definition: action.cc:123
smash::Action::outgoing_particles_
ParticleList outgoing_particles_
Initially this stores only the PDG codes of final-state particles.
Definition: action.h:312
smash::FourVector::abs
double abs() const
calculate the lorentz invariant absolute value
Definition: fourvector.h:454
smash::Action::sqrt_s
double sqrt_s() const
Determine the total energy in the center-of-mass frame [GeV].
Definition: action.h:266
smash::ProcessType::Elastic
elastic scattering: particles remain the same, only momenta change
smash::Action::get_interaction_point
FourVector get_interaction_point() const
Get the interaction point.
Definition: action.cc:69
smash::pdg::p
constexpr int p
Proton.
Definition: pdgcode_constants.h:28
smash::random::uniform
T uniform(T min, T max)
Definition: random.h:88
smash::ID_PROCESS_PHOTON
constexpr std::uint32_t ID_PROCESS_PHOTON
Process ID for any photon process.
Definition: constants.h:131
smash::pot_pointer
Potentials * pot_pointer
Pointer to a Potential class.
Definition: potential_globals.cc:17
smash::Action::sample_masses
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....
Definition: action.cc:154
smash::random::canonical
T canonical()
Definition: random.h:113
smash::FourVector::threevec
ThreeVector threevec() const
Definition: fourvector.h:319