nucleus.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2014-2022
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#ifndef SRC_INCLUDE_SMASH_NUCLEUS_H_
#define SRC_INCLUDE_SMASH_NUCLEUS_H_
 
#include <map>
#include <stdexcept>
#include <vector>
 
#include "configuration.h"
#include "constants.h"
#include "forwarddeclarations.h"
#include "fourvector.h"
#include "particledata.h"
#include "threevector.h"
 
namespace smash {
 
class Nucleus {
 public:
  Nucleus() = default;
 
  Nucleus(Configuration &config, int nTest);
 
  Nucleus(const std::map<PdgCode, int> &particle_list, int nTest);
  virtual ~Nucleus() = default;
 
  double mass() const;
 
  virtual ThreeVector distribute_nucleon();
 
  double woods_saxon(double x);
 
  virtual void arrange_nucleons();
 
  virtual void set_parameters_automatic();
 
  virtual void set_parameters_from_config(Configuration &config);
 
  virtual void generate_fermi_momenta();
 
  void boost(double beta_scalar);
 
  void fill_from_list(const std::map<PdgCode, int> &particle_list,
                      int testparticles);
 
  void shift(double z_offset, double x_offset, double simulation_time);
 
  virtual void rotate() {}
 
  void copy_particles(Particles *particles);
 
  inline size_t size() const { return particles_.size(); }
 
  inline size_t number_of_particles() const {
    size_t nop = particles_.size() / testparticles_;
    /* If size() is not a multiple of testparticles_, this will throw an
     * error. */
    if (nop * testparticles_ != particles_.size()) {
      throw TestparticleConfusion(
          "Number of test particles and test particles"
          "per particle are incompatible.");
    }
    return nop;
  }
 
  inline size_t number_of_protons() const {
    size_t proton_counter = 0;
    /* If n_protons is not a multiple of testparticles_, this will throw an
     * error. */
    for (auto &particle : particles_) {
      if (particle.type().pdgcode() == pdg::p) {
        proton_counter++;
      }
    }
 
    size_t n_protons = proton_counter / testparticles_;
 
    if (n_protons * testparticles_ != proton_counter) {
      throw TestparticleConfusion(
          "Number of test protons and test particles"
          "per proton are incompatible.");
    }
 
    return n_protons;
  }
 
  FourVector center() const;
 
  void set_label(BelongsTo label) {
    for (ParticleData &data : particles_) {
      data.set_belongs_to(label);
    }
  }
 
  void align_center() {
    FourVector centerpoint = center();
    for (auto p = particles_.begin(); p != particles_.end(); ++p) {
      p->set_4position(p->position() - centerpoint);
    }
  }
 
  // This function as well as nucleon_density_unnormalized could in principle
  // be defined without the second argument
  virtual double nucleon_density(double r, double) const;
  virtual double nucleon_density_unnormalized(double r, double) const;
  virtual double calculate_saturation_density() const;
  virtual void set_saturation_density(double density) {
    saturation_density_ = density;
  }
 
  struct TestparticleConfusion : public std::length_error {
    using std::length_error::length_error;
  };
 
 private:
  double diffusiveness_;
  double nuclear_radius_;
  double proton_radius_ = 1.2;
  size_t testparticles_ = 1;
 
 protected:
  std::vector<ParticleData> particles_;
 
  // Needed as public member for inheritance to deformed nuclei
  double saturation_density_ = nuclear_density;
 
  void random_euler_angles();
 
  double euler_phi_;
  double euler_theta_;
  double euler_psi_;
 
 public:
  inline std::vector<ParticleData>::iterator begin() {
    return particles_.begin();
  }
  inline std::vector<ParticleData>::iterator end() { return particles_.end(); }
  inline std::vector<ParticleData>::const_iterator cbegin() const {
    return particles_.cbegin();
  }
  inline std::vector<ParticleData>::const_iterator cend() const {
    return particles_.cend();
  }
  inline void set_diffusiveness(double diffuse) { diffusiveness_ = diffuse; }
  inline double get_diffusiveness() const { return diffusiveness_; }
  inline double get_saturation_density() const { return saturation_density_; }
  inline double default_nuclear_radius() {
    int A = number_of_particles();
 
    if (A <= 16) {
      // radius: rough guess for all nuclei not listed explicitly with A <= 16
      return (proton_radius_ * std::cbrt(A));
    } else {
      // radius taken from \iref{Rybczynski:2013yba}
      return (1.12 * std::pow(A, 1.0 / 3.0) - 0.86 * std::pow(A, -1.0 / 3.0));
    }
  }
  inline void set_nuclear_radius(double rad) { nuclear_radius_ = rad; }
  inline double get_nuclear_radius() const { return nuclear_radius_; }
  friend std::ostream &operator<<(std::ostream &, const Nucleus &);
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_NUCLEUS_H_