particletype.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2012-2019
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
#ifndef SRC_INCLUDE_PARTICLETYPE_H_
#define SRC_INCLUDE_PARTICLETYPE_H_

#include <cassert>
#include <string>
#include <utility>
#include <vector>

#include "forwarddeclarations.h"
#include "macros.h"
#include "pdgcode.h"

namespace smash {

enum class Parity {
  Pos,
  Neg
};

enum class WhichDecaymodes {
  All,
  Hadronic,
  Dileptons
};

inline Parity operator-(Parity p) {
  switch (p) {
    case Parity::Pos:
      return Parity::Neg;
    case Parity::Neg:
      return Parity::Pos;
  }
  // This is unreachable and should be optimized away.
  // It is required to silence a compiler warning.
  throw std::runtime_error("unreachable");
}

inline Parity operator*(Parity x, Parity y) {
  if (x == y) {
    return Parity::Pos;
  } else {
    return Parity::Neg;
  }
}

inline void operator*=(Parity &x, Parity y) {
  if (x == y) {
    x = Parity::Pos;
  } else {
    x = Parity::Neg;
  }
}

class ParticleType {
 public:
  // If this is changed, make sure to update the userguide in
  // `include/configuration.h`.
  static constexpr double width_cutoff = 1e-5;

  ParticleType(std::string n, double m, double w, Parity p, PdgCode id);

  ParticleType(const ParticleType &) = delete;
  ParticleType &operator=(const ParticleType &) = delete;

  ParticleType(ParticleType &&) = default;
  ParticleType &operator=(ParticleType &&) = default;

  const DecayModes &decay_modes() const;

  const std::string &name() const { return name_; }

  double mass() const { return mass_; }

  double mass_sqr() const { return mass_ * mass_; }

  double width_at_pole() const { return width_; }

  Parity parity() const { return parity_; }

  PdgCode pdgcode() const { return pdgcode_; }

  bool has_antiparticle() const { return pdgcode_.has_antiparticle(); }

  ParticleTypePtr get_antiparticle() const;

  int antiparticle_sign() const { return pdgcode_.antiparticle_sign(); }

  int isospin() const;

  int isospin3() const { return I3_; }

  double isospin3_rel() const {
    unsigned int I = isospin();
    return (I == 0) ? 0 : static_cast<double>(isospin3()) / I;
  }

  IsoParticleType *iso_multiplet() const { return iso_multiplet_; }

  int32_t charge() const { return charge_; }

  unsigned int spin() const { return pdgcode_.spin(); }

  bool is_hadron() const { return pdgcode_.is_hadron(); }

  bool is_lepton() const { return pdgcode_.is_lepton(); }

  bool is_baryon() const { return pdgcode_.is_baryon(); }

  bool is_meson() const { return pdgcode_.is_meson(); }

  int baryon_number() const { return pdgcode_.baryon_number(); }

  int strangeness() const { return pdgcode_.strangeness(); }

  bool is_nucleon() const { return pdgcode_.is_nucleon(); }

  bool is_Delta() const { return pdgcode_.is_Delta(); }

  bool is_rho() const { return pdgcode_.is_rho(); }

  inline bool is_Nstar() const {
    return is_baryon() && isospin() == 1 && !pdgcode_.is_nucleon() &&
           pdgcode_.strangeness() == 0 && pdgcode_.charmness() == 0;
  }

  bool is_Nstar1535() const { return pdgcode_.is_Nstar1535(); }

  inline bool is_Deltastar() const {
    return is_baryon() && isospin() == 3 && !pdgcode_.is_Delta() &&
           pdgcode_.strangeness() == 0 && pdgcode_.charmness() == 0;
  }

  inline bool is_stable() const { return width_ < width_cutoff; }

  inline bool is_nucleus() const { return pdgcode_.is_nucleus(); }

  inline bool is_deuteron() const { return pdgcode_.is_deuteron(); }

  inline bool is_dprime() const {
    return is_nucleus() && std::abs(pdgcode_.get_decimal()) == 1000010021;
  }

  double min_mass_kinematic() const;

  double min_mass_spectral() const;

  double partial_width(const double m, const DecayBranch *mode) const;

  double total_width(const double m) const;

  bool wanted_decaymode(const DecayType &t, WhichDecaymodes wh) const;

  DecayBranchList get_partial_widths(const FourVector p, const ThreeVector x,
                                     WhichDecaymodes wh) const;

  double get_partial_width(const double m, const ParticleType &t_a,
                           const ParticleType &t_b) const;

  double get_partial_in_width(const double m, const ParticleData &p_a,
                              const ParticleData &p_b) const;

  double spectral_function(double m) const;

  double spectral_function_no_norm(double m) const;

  double spectral_function_const_width(double m) const;

  double spectral_function_simple(double m) const;

  double sample_resonance_mass(const double mass_stable,
                               const double cms_energy, int L = 0) const;

  std::pair<double, double> sample_resonance_masses(const ParticleType &t2,
                                                    const double cms_energy,
                                                    int L = 0) const;

  void dump_width_and_spectral_function() const;

  static const ParticleTypeList &list_all();

  static ParticleTypePtrList &list_nucleons();
  static ParticleTypePtrList &list_anti_nucleons();
  static ParticleTypePtrList &list_Deltas();
  static ParticleTypePtrList &list_anti_Deltas();
  static ParticleTypePtrList &list_baryon_resonances();
  static ParticleTypePtrList &list_light_nuclei();

  static const ParticleTypePtr try_find(PdgCode pdgcode);

  static const ParticleType &find(PdgCode pdgcode);

  struct PdgNotFoundFailure : public std::runtime_error {
    using std::runtime_error::runtime_error;
  };

  static bool exists(PdgCode pdgcode);

  static bool exists(const std::string &name);

  static void create_type_list(const std::string &particles);

  bool operator==(const ParticleType &rhs) const {
    return pdgcode() == rhs.pdgcode();
  }
  bool operator!=(const ParticleType &rhs) const {
    return pdgcode() != rhs.pdgcode();
  }
  bool operator<(const ParticleType &rhs) const {
    return pdgcode() < rhs.pdgcode();
  }

  static void check_consistency();

  ParticleTypePtr operator&() const;

  struct LoadFailure : public std::runtime_error {
    using std::runtime_error::runtime_error;
  };

 private:
  std::string name_;
  double mass_;
  double width_;
  Parity parity_;
  PdgCode pdgcode_;
  mutable double min_mass_kinematic_;
  mutable double min_mass_spectral_;
  mutable double norm_factor_ = -1.;
  int32_t charge_;
  mutable int isospin_;
  int I3_;

  IsoParticleType *iso_multiplet_ = nullptr;

  mutable double max_factor1_ = 1.;
  mutable double max_factor2_ = 1.;

  friend std::ostream &operator<<(std::ostream &out, const ParticleType &type);
};

class ParticleTypePtr {
 public:
  const ParticleType &operator*() const { return lookup(); }

  const ParticleType *operator->() const {
    // this requires std::addressof because &lookup() would call
    // ParticleType::operator& and return ParticleTypePtr again
    return std::addressof(lookup());
  }

  ParticleTypePtr() = default;

  bool operator==(const ParticleTypePtr &rhs) const {
    return index_ == rhs.index_;
  }

  bool operator!=(const ParticleTypePtr &rhs) const {
    return index_ != rhs.index_;
  }
  bool operator<(const ParticleTypePtr &rhs) const {
    return index_ < rhs.index_;
  }

  operator bool() const { return index_ != 0xffff; }

 private:
  friend ParticleTypePtr ParticleType::operator&() const;

  explicit ParticleTypePtr(std::uint16_t i) : index_(i) {}

  const ParticleType &lookup() const {
    assert(index_ != 0xffff);
    return ParticleType::list_all()[index_];
  }

  std::uint16_t index_ = 0xffff;
};

inline ParticleTypePtr ParticleType::get_antiparticle() const {
  assert(has_antiparticle());
  return &find(pdgcode_.get_antiparticle());
}

}  // namespace smash

#endif  // SRC_INCLUDE_PARTICLETYPE_H_