action.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#ifndef SRC_INCLUDE_ACTION_H_
#define SRC_INCLUDE_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 {

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 Particles &particles,
                        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_; }

  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;

  class InvalidResonanceFormation : public std::invalid_argument {
    using std::invalid_argument::invalid_argument;
  };

 protected:
  ParticleList incoming_particles_;

  ParticleList outgoing_particles_;

  const double time_of_execution_;

  ProcessType process_type_;

  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) {
    const auto &log = logger<LogArea::Action>();
    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. */
    log.fatal(source_location,
              "Problem in choose_channel: ", subprocesses.size(), " ",
              weight_sum, " ", total_weight, " ",
              //          random_weight, "\n", *this);
              random_weight, "\n");
    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 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_ACTION_H_