action.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2021
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_ACTION_H_
#define SRC_INCLUDE_SMASH_ACTION_H_
 
#include <stdexcept>
#include <utility>
#include <vector>
 
#include "lattice.h"
#include "particles.h"
#include "pauliblocking.h"
#include "potentials.h"
#include "processbranch.h"
#include "random.h"
 
namespace smash {
static constexpr int LAction = LogArea::Action::id;
 
class Action {
 public:
  Action(const ParticleList &in_part, double time)
      : incoming_particles_(in_part),
        time_of_execution_(time + in_part[0].position().x0()) {}
 
  Action(const ParticleData &in_part, const ParticleData &out_part, double time,
         ProcessType type)
      : incoming_particles_({in_part}),
        outgoing_particles_({out_part}),
        time_of_execution_(time + in_part.position().x0()),
        process_type_(type) {}
 
  Action(const ParticleList &in_part, const ParticleList &out_part,
         double absolute_execution_time, ProcessType type)
      : incoming_particles_(std::move(in_part)),
        outgoing_particles_(std::move(out_part)),
        time_of_execution_(absolute_execution_time),
        process_type_(type) {}
 
  Action(const Action &) = delete;
 
  virtual ~Action();
 
  bool operator<(const Action &rhs) const {
    return time_of_execution_ < rhs.time_of_execution_;
  }
 
  virtual double get_total_weight() const = 0;
 
  virtual double get_partial_weight() const = 0;
 
  virtual ProcessType get_type() const { return process_type_; }
 
  template <typename Branch>
  void add_process(ProcessBranchPtr<Branch> &p,
                   ProcessBranchList<Branch> &subprocesses,
                   double &total_weight) {
    if (p->weight() > 0) {
      total_weight += p->weight();
      subprocesses.emplace_back(std::move(p));
    }
  }
 
  template <typename Branch>
  void add_processes(ProcessBranchList<Branch> pv,
                     ProcessBranchList<Branch> &subprocesses,
                     double &total_weight) {
    subprocesses.reserve(subprocesses.size() + pv.size());
    for (auto &proc : pv) {
      if (proc->weight() > 0) {
        total_weight += proc->weight();
        subprocesses.emplace_back(std::move(proc));
      }
    }
  }
 
  virtual void generate_final_state() = 0;
 
  virtual void perform(Particles *particles, uint32_t id_process);
 
  bool is_valid(const Particles &particles) const;
 
  bool is_pauli_blocked(const std::vector<Particles> &ensembles,
                        const PauliBlocker &p_bl) const;
 
  const ParticleList &incoming_particles() const;
 
  void update_incoming(const Particles &particles);
 
  const ParticleList &outgoing_particles() const { return outgoing_particles_; }
 
  double time_of_execution() const { return time_of_execution_; }
 
  virtual void check_conservation(const uint32_t id_process) const;
 
  double sqrt_s() const { return total_momentum().abs(); }
 
  FourVector total_momentum_of_outgoing_particles() const;
 
  FourVector get_interaction_point() const;
 
  std::pair<FourVector, FourVector> get_potential_at_interaction_point() const;
 
  void set_stochastic_pos_idx() {
    const int max_inc_idx = incoming_particles_.size() - 1;
    stochastic_position_idx_ = random::uniform_int(0, max_inc_idx);
  }
 
  static double lambda_tilde(double a, double b, double c) {
    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;
  }
 
  class InvalidResonanceFormation : public std::invalid_argument {
    using std::invalid_argument::invalid_argument;
  };
 
  static void sample_manybody_phasespace_impl(
      double sqrts, const std::vector<double> &m,
      std::vector<FourVector> &sampled_momenta);
 
 protected:
  ParticleList incoming_particles_;
 
  ParticleList outgoing_particles_;
 
  const double time_of_execution_;
 
  ProcessType process_type_;
 
  double box_length_ = -1.0;
 
  int stochastic_position_idx_ = -1;
 
  FourVector total_momentum() const {
    FourVector mom(0.0, 0.0, 0.0, 0.0);
    for (const auto &p : incoming_particles_) {
      mom += p.momentum();
    }
    return mom;
  }
 
  template <typename Branch>
  const Branch *choose_channel(const ProcessBranchList<Branch> &subprocesses,
                               double total_weight) {
    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();
  }
 
  virtual std::pair<double, double> sample_masses(
      double kinetic_energy_cm) const;
 
  virtual void sample_angles(std::pair<double, double> masses,
                             double kinetic_energy_cm);
 
  void sample_2body_phasespace();
 
  virtual void sample_manybody_phasespace();
 
  void assign_formation_time_to_outgoing_particles();
 
  virtual void format_debug_output(std::ostream &out) const = 0;
 
  friend std::ostream &operator<<(std::ostream &out, const Action &action) {
    action.format_debug_output(out);
    return out;
  }
 
 private:
  const ParticleType &type_of_pout(const ParticleData &p_out) const {
    return p_out.type();
  }
  const ParticleType &type_of_pout(const ParticleTypePtr &p_out) const {
    return *p_out;
  }
};
 
inline std::vector<ActionPtr> &operator+=(std::vector<ActionPtr> &lhs,
                                          std::vector<ActionPtr> &&rhs) {
  if (lhs.size() == 0) {
    lhs = std::move(rhs);
  } else {
    lhs.insert(lhs.end(), std::make_move_iterator(rhs.begin()),
               std::make_move_iterator(rhs.end()));
  }
  return lhs;
}
 
inline std::ostream &operator<<(std::ostream &out, const ActionPtr &action) {
  return out << *action;
}
 
std::ostream &operator<<(std::ostream &out, const ActionList &actions);
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_ACTION_H_