pdgcode.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2020
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_PDGCODE_H_
#define SRC_INCLUDE_SMASH_PDGCODE_H_
 
#include <algorithm>
#include <array>
#include <cassert>
#include <cstdlib>
#include <iosfwd>
#include <sstream>
#include <stdexcept>
#include <string>
 
#include <iostream>
 
#include "pdgcode_constants.h"
 
namespace smash {
 
class PdgCode {
 public:
  struct InvalidPdgCode : public std::invalid_argument {
    using std::invalid_argument::invalid_argument;
  };
 
  /****************************************************************************
   *                                                                          *
   * First, the constructors                                                  *
   *                                                                          *
   ****************************************************************************/
 
  PdgCode() : dump_(0x0) {}
  explicit PdgCode(const std::string& codestring) {
    set_from_string(codestring);
  }
 
  PdgCode(std::int32_t codenumber) : dump_(0x0) {  // NOLINT(runtime/explicit)
    digits_.antiparticle_ = false;
    if (codenumber < 0) {
      digits_.antiparticle_ = true;
      codenumber = -codenumber;
    }
    set_fields(codenumber);
  }
  explicit PdgCode(const std::uint32_t abscode) : dump_(0x0) {
    // use the first bit for the antiparticle_ boolean.
    digits_.antiparticle_ = ((abscode & 0x80000000u) != 0);
    set_fields(abscode);
  }
 
  /****************************************************************************
   *                                                                          *
   * test function and export functions                                       *
   *                                                                          *
   ****************************************************************************/
 
  inline int test_code() const {
    int fail = 0;
    if (digits_.n_ > 9) {
      fail |= 1 << 6;
    }
    if (digits_.n_R_ > 9) {
      fail |= 1 << 5;
    }
    if (digits_.n_L_ > 9) {
      fail |= 1 << 4;
    }
    if (digits_.n_q1_ > 6) {
      fail |= 1 << 3;
    }
    if (digits_.n_q2_ > 6) {
      fail |= 1 << 2;
    }
    if (digits_.n_q3_ > 6) {
      fail |= 1 << 1;
    }
    if (digits_.n_J_ > 15) {
      fail |= 1;
    }
    return fail;
  }
 
  void check() const {
    // n_J must be odd for mesons and even for baryons (and cannot be zero)
    if (is_hadron()) {
      if (baryon_number() == 0) {
        // mesons: special cases K0_L=0x130 and K0_S=0x310
        if ((digits_.n_J_ % 2 == 0) && dump() != 0x130 && dump() != 0x310) {
          throw InvalidPdgCode("Invalid PDG code " + string() +
                               " (meson with even n_J)");
        }
      } else {
        if ((digits_.n_J_ % 2 != 0) || digits_.n_J_ == 0) {
          throw InvalidPdgCode("Invalid PDG code " + string() +
                               " (baryon with odd n_J)");
        }
      }
    } else {
      if (digits_.n_J_ == 0 && dump() != 0x0) {
        throw InvalidPdgCode("Invalid PDG code " + string() + " (n_J==0)");
      }
    }
    /* The antiparticle flag only makes sense for particle types
     * that have an antiparticle. */
    if (digits_.antiparticle_ && !has_antiparticle()) {
      throw InvalidPdgCode("Invalid PDG code " + string() +
                           " (cannot be negative)");
    }
  }
 
  inline std::uint32_t dump() const {
    // this cuts the three unused bits.
    return (dump_ & 0x8fffffff);
  }
 
  inline std::int32_t code() const { return antiparticle_sign() * ucode(); }
 
  inline std::string string() const {
    std::stringstream ss;
    ss << get_decimal();
    return ss.str();
  }
 
  PdgCode get_antiparticle() const {
    PdgCode result = *this;
    result.digits_.antiparticle_ = !digits_.antiparticle_;
    return result;
  }
 
  static PdgCode from_decimal(const int pdgcode_decimal) {
    // Nucleus and special codes with 2J+1 > 9
    if (std::abs(pdgcode_decimal) > 1E7) {
      return PdgCode(std::to_string(pdgcode_decimal));
    }
    int a = pdgcode_decimal;
    int hex_pdg = 0, tmp = 1;
    while (a) {
      hex_pdg += (a % 10) * tmp;
      tmp *= 16;
      a = a / 10;
    }
    return PdgCode(hex_pdg);
  }
 
  /****************************************************************************
   *                                                                          *
   * accessors of various properties                                          *
   *                                                                          *
   ****************************************************************************/
 
  inline bool is_nucleus() const {
    assert(digits_.is_nucleus_ == nucleus_.is_nucleus_);
    return nucleus_.is_nucleus_;
  }
 
  inline bool is_hadron() const {
    return (digits_.n_q3_ != 0 && digits_.n_q2_ != 0 && !is_nucleus());
  }
 
  inline bool is_lepton() const {
    return (digits_.n_q1_ == 0 && digits_.n_q2_ == 0 && digits_.n_q3_ == 1 &&
            !is_nucleus());
  }
 
  inline int baryon_number() const {
    if (is_nucleus()) {
      return static_cast<int>(nucleus_.A_) * antiparticle_sign();
    }
    if (!is_hadron() || digits_.n_q1_ == 0) {
      return 0;
    }
    return antiparticle_sign();
  }
  inline bool is_baryon() const { return is_hadron() && digits_.n_q1_ != 0; }
 
  inline bool is_meson() const { return is_hadron() && digits_.n_q1_ == 0; }
 
  inline bool is_nucleon() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::p || abs_code == pdg::n);
  }
 
  inline bool is_proton() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::p);
  }
 
  inline bool is_neutron() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::n);
  }
 
  inline bool is_Nstar1535() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::N1535_p || abs_code == pdg::N1535_z);
  }
 
  inline bool is_Delta() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::Delta_pp || abs_code == pdg::Delta_p ||
            abs_code == pdg::Delta_z || abs_code == pdg::Delta_m);
  }
 
  inline bool is_hyperon() const { return is_hadron() && digits_.n_q1_ == 3; }
 
  inline bool is_Omega() const {
    return is_hyperon() && digits_.n_q2_ == 3 && digits_.n_q3_ == 3;
  }
 
  inline bool is_Xi() const {
    return is_hyperon() && digits_.n_q2_ == 3 && digits_.n_q3_ != 3;
  }
 
  inline bool is_Lambda() const {
    return is_hyperon() && digits_.n_q2_ == 1 && digits_.n_q3_ == 2;
  }
 
  inline bool is_Sigma() const {
    return is_hyperon() && digits_.n_q2_ != 3 && !is_Lambda();
  }
 
  inline bool is_kaon() const {
    const auto abs_code = std::abs(code());
    return (abs_code == pdg::K_p) || (abs_code == pdg::K_z);
  }
 
  inline bool is_pion() const {
    const auto c = code();
    return (c == pdg::pi_z) || (c == pdg::pi_p) || (c == pdg::pi_m);
  }
 
  inline bool is_rho() const {
    const auto c = code();
    return (c == pdg::rho_z) || (c == pdg::rho_p) || (c == pdg::rho_m);
  }
 
  inline bool is_deuteron() const {
    const int dec = get_decimal();
    return is_nucleus() && (dec == pdg::decimal_d || dec == pdg::decimal_antid);
  }
 
  bool has_antiparticle() const {
    if (is_nucleus()) {
      return true;
    }
    if (is_hadron()) {
      return (baryon_number() != 0) || (digits_.n_q2_ != digits_.n_q3_);
    } else {
      return digits_.n_q3_ == 1;  // leptons!
    }
  }
 
  inline int isospin3() const {
    /* net_quark_number(2) is the number of u quarks,
     * net_quark_number(1) is the number of d quarks. */
    return net_quark_number(2) - net_quark_number(1);
  }
 
  inline double frac_strange() const {
    /* The quarkonium state has 0 net strangeness
     *  but there are actually 2 strange quarks out of 2 total */
    if (is_hadron() && digits_.n_q3_ == 3 && digits_.n_q2_ == 3) {
      return 1.;
    } else {
      // For all other cases, there isn't both a strange and anti-strange
      if (is_baryon()) {
        return std::abs(strangeness()) / 3.;
      } else if (is_meson()) {
        return std::abs(strangeness()) / 2.;
      } else {
        /* If not baryon or meson, this should be 0, as AQM does not
         * extend to non-hadrons */
        return 0.;
      }
    }
  }
 
  inline int strangeness() const { return -net_quark_number(3); }
 
  inline int charmness() const { return +net_quark_number(4); }
 
  inline int bottomness() const { return -net_quark_number(5); }
 
  int charge() const {
    if (is_hadron() || is_nucleus()) {
      // Q will accumulate 3*charge (please excuse the upper case. I
      // want to distinguish this from q which might be interpreted as
      // shorthand for "quark".)
      int Q = 0;
      /* This loops over d,u,s,c,b,t quarks (the latter can be safely ignored,
       * but I don't think this will be a bottle neck. */
      for (int i = 1; i < 7; i++) {
        /* u,c,t quarks have charge = 2/3 e, while d,s,b quarks have -1/3 e.
         * The antiparticle sign is already in net_quark_number. */
        Q += (i % 2 == 0 ? 2 : -1) * net_quark_number(i);
      }
      return Q / 3;
    }
    /* non-hadron:
     * Leptons: 11, 13, 15 are e, μ, τ and have a charge -1, while
     *          12, 14, 16 are the neutrinos that have no charge. */
    if (digits_.n_q3_ == 1) {
      return -1 * (digits_.n_J_ % 2) * antiparticle_sign();
    }
    /* Bosons: 24 is the W+, all else is uncharged.
     * we ignore the first digits so that this also finds strange gauge
     * boson "resonances" (in particular, \f$\tilde \chi_1^+\f$ with PDG
     * Code 1000024). */
    if ((dump_ & 0x0000ffff) == 0x24) {
      return antiparticle_sign();
    }
    // default (this includes all other Bosons) is 0.
    return 0;
  }
 
  inline unsigned int spin() const {
    if (is_nucleus()) {
      /* Currently the only nucleus I care about is deutron,
       * which has spin one. */
      return 2;
    }
 
    if (is_hadron()) {
      if (digits_.n_J_ == 0) {
        return 0;  // special cases: K0_L=0x130 & K0_S=0x310
      } else {
        return digits_.n_J_ - 1;
      }
    }
    /* this assumes that we only have white particles (no single
     * quarks): Electroweak fermions have 11-17, so the
     * second-to-last-digit is the spin. The same for the Bosons: they
     * have 21-29 and 2spin = 2 (this fails for the Higgs). */
    return digits_.n_q3_;
  }
  inline unsigned int spin_degeneracy() const {
    if (is_hadron() && digits_.n_J_ > 0) {
      return digits_.n_J_;
    }
    return spin() + 1;
  }
  inline int antiparticle_sign() const {
    return (digits_.antiparticle_ ? -1 : +1);
  }
  inline std::int32_t quarks() const {
    if (!is_hadron() || is_nucleus()) {
      return 0;
    }
    return chunks_.quarks_;
  }
 
  std::array<int, 3> quark_content() const {
    std::array<int, 3> result = {static_cast<int>(digits_.n_q1_),
                                 static_cast<int>(digits_.n_q2_),
                                 static_cast<int>(digits_.n_q3_)};
    if (is_hadron()) {
      // Antibaryons
      if (digits_.n_q1_ != 0 && digits_.antiparticle_) {
        for (size_t i = 0; i < 3; i++) {
          result[i] = -result[i];
        }
      }
      // Mesons
      if (digits_.n_q1_ == 0) {
        // Own antiparticle
        if (digits_.n_q2_ == digits_.n_q3_) {
          result[2] = -result[2];
        } else {
          // Like pi-
          if (digits_.antiparticle_) {
            result[1] = -result[1];
            // Like pi+
          } else {
            result[2] = -result[2];
          }
        }
        // add extra minus sign according to the pdg convention
        if (digits_.n_q2_ != digits_.n_q3_ && digits_.n_q2_ % 2 == 1) {
          for (int i = 1; i <= 2; i++) {
            result[i] = -result[i];
          }
        }
      }
    } else {
      result = {0, 0, 0};
    }
    return result;
  }
 
  bool contains_enough_valence_quarks(int valence_quarks_required) const;
 
  /****************************************************************************
   *                                                                          *
   * operations with more than one PDG Code                                   *
   *                                                                          *
   ****************************************************************************/
 
  inline bool operator<(const PdgCode rhs) const {
    return dump_ < rhs.dump_;
    /* the complex thing to do here is to calculate:
     *   code() < rhs.code()
     * but for getting a total order that's overkill. The uint32_t value in
     * dump_ works just fine. */
  }
 
  inline bool operator==(const PdgCode rhs) const { return dump_ == rhs.dump_; }
 
  inline bool operator!=(const PdgCode rhs) const { return !(*this == rhs); }
 
  inline bool is_antiparticle_of(const PdgCode rhs) const {
    return code() == -rhs.code();
  }
 
  friend std::istream& operator>>(std::istream& is, PdgCode& code);
 
  static PdgCode invalid() { return PdgCode(0x0); }
 
  int32_t get_decimal() const {
    if (is_nucleus()) {
      // ±10LZZZAAAI
      return antiparticle_sign() *
             (nucleus_.I_ + 10 * nucleus_.A_ + 10000 * nucleus_.Z_ +
              10000000 * nucleus_.n_Lambda_ + 1000000000);
    }
    int n_J_1 = 0;
    int n_J_2 = digits_.n_J_;
    if (n_J_2 > 9) {
      n_J_1 = n_J_2 - 9;
      n_J_2 = 9;
    }
    return antiparticle_sign() *
           (n_J_2 + digits_.n_q3_ * 10 + digits_.n_q2_ * 100 +
            digits_.n_q1_ * 1000 + digits_.n_L_ * 10000 +
            digits_.n_R_ * 100000 + digits_.n_ * 1000000 + n_J_1 * 10000000);
  }
 
  void deexcite() {
    if (!is_nucleus()) {
      chunks_.excitation_ = 0;
    } else {
      nucleus_.I_ = 0;
    }
  }
 
  int net_quark_number(const int quark) const;
 
 private:
/* amend this line with something that identifies your compiler if its
 * bit field order is like in the gnu c compiler for 64 bit
 * architectures (if you are unsure, try one and check the pdgcode
 * test). */
#if defined(__GNUC__) && (defined(__x86_64__) || defined(__i386__)) || \
    defined(DOXYGEN)
#define SMASH_BITFIELD_ORDER_ 1
/* put your compiler here if the bit field order is reversed w.r.t. gnu
 * c compiler for 64 bit. */
#elif defined(__OTHER_COMPILER__)
#define SMASH_BITFIELD_ORDER_ 2
#else
#error "Please determine the correct bit-field order for your target/compiler"
#endif
 
  union {
    struct {
#if SMASH_BITFIELD_ORDER_ == 1
      std::uint32_t n_J_ : 4;
      std::uint32_t n_q3_ : 4;
      std::uint32_t n_q2_ : 4;
      std::uint32_t n_q1_ : 4;
      std::uint32_t n_L_ : 4;
      std::uint32_t n_R_ : 4;
      std::uint32_t n_ : 4, : 2;
      bool is_nucleus_ : 1;
      bool antiparticle_ : 1;
#else  // reverse ordering
      bool antiparticle_ : 1;
      bool is_nucleus : 1, : 2;
      std::uint32_t n_ : 4;
      std::uint32_t n_R_ : 4;
      std::uint32_t n_L_ : 4;
      std::uint32_t n_q1_ : 4;
      std::uint32_t n_q2_ : 4;
      std::uint32_t n_q3_ : 4;
      std::uint32_t n_J_ : 4;
#endif
    } digits_;
    std::uint32_t dump_;
    struct {
#if SMASH_BITFIELD_ORDER_ == 1
      std::uint32_t : 4;
      std::uint32_t quarks_ : 12;
      std::uint32_t excitation_ : 12, : 4;
#else  
      std::uint32_t : 4, excitation_ : 12;
      std::uint32_t quarks_ : 12, : 4;
#endif
    } chunks_;
    struct {
#if SMASH_BITFIELD_ORDER_ == 1
      std::uint32_t n_Lambda_ : 6;
      std::uint32_t Z_ : 10;
      std::uint32_t A_ : 10;
      std::uint32_t I_ : 4;
      bool is_nucleus_ : 1;
      bool antiparticle_ : 1;
#else  // Reverse ordering
      bool antiparticle_ : 1;
      bool is_nucleus_ : 1;
      std::uint32_t I_ : 4;
      std::uint32_t A_ : 10;
      std::uint32_t Z_ : 10;
      std::uint32_t n_Lambda_ : 6;
#endif
    } nucleus_;
  };
 
  inline std::uint32_t ucode() const { return (dump_ & 0x0fffffff); }
 
  inline std::uint32_t get_digit_from_char(const char inp) const {
    // Decimal digit
    if (48 <= inp && inp <= 57) {
      return inp - 48;
    }
    // Hexdecimal digit, uppercase
    if (65 <= inp && inp <= 70) {
      return inp - 65 + 10;
    }
    // Hexdecimal digit, lowercase
    if (97 <= inp && inp <= 102) {
      return inp - 97 + 10;
    }
    throw InvalidPdgCode("PdgCode: Invalid character " + std::string(&inp, 1) +
                         " found.\n");
  }
 
  inline void set_from_string(const std::string& codestring) {
    dump_ = 0;
    // Implicit with the above: digits_.antiparticle_ = false;
    digits_.n_ = digits_.n_R_ = digits_.n_L_ = digits_.n_q1_ = digits_.n_q2_ =
        digits_.n_q3_ = digits_.n_J_ = digits_.is_nucleus_ = 0;
    size_t length = codestring.size();
    if (length < 1) {
      throw InvalidPdgCode("Empty string does not contain PDG Code\n");
    }
    int c = 0;
    /* Look at current character; if it is a + or minus sign, read it
     * and advance to next char. */
    if (codestring[c] == '-') {
      digits_.antiparticle_ = true;
      ++c;
    } else if (codestring[c] == '+') {
      digits_.antiparticle_ = false;
      ++c;
    }
    // Save if the first character was a sign:
    unsigned int sign = c;
 
    // Nucleus
    if (length == 10 + sign) {
      nucleus_.is_nucleus_ = true;
      if (codestring[c] != '1' || codestring[c + 1] != '0') {
        throw InvalidPdgCode("Pdg code of nucleus \"" + codestring +
                             "\" should start with 10\n");
      }
      c += 2;
      // ±10LZZZAAAI is the standard for nuclei
      std::array<int, 8> digits;
      for (int i = 0; i < 8; i++) {
        digits[i] = get_digit_from_char(codestring[c + i]);
      }
      nucleus_.n_Lambda_ = digits[0];
      nucleus_.Z_ = 100 * digits[1] + 10 * digits[2] + digits[3];
      nucleus_.A_ = 100 * digits[4] + 10 * digits[5] + digits[6];
      nucleus_.I_ = digits[7];
      return;
    }
 
    // Codestring shouldn't be longer than 8 + sign, except for nuclei
    if (length > 8 + sign) {
      throw InvalidPdgCode("String \"" + codestring +
                           "\" too long for PDG Code\n");
    }
    /* Please note that in what follows, we actually need c++, not ++c.
     * first digit is used for n_J if the last digit is not enough. */
    if (length > 7 + sign) {
      digits_.n_J_ += get_digit_from_char(codestring[c++]);
    }
    // Codestring has 7 digits? 7th from last goes in n_.
    if (length > 6 + sign) {
      digits_.n_ = get_digit_from_char(codestring[c++]);
    }
    // It has 6 or 7 digits? 6th from last is n_R_.
    if (length > 5 + sign) {
      digits_.n_R_ = get_digit_from_char(codestring[c++]);
    }
    // 5th from last is n_L_.
    if (length > 4 + sign) {
      digits_.n_L_ = get_digit_from_char(codestring[c++]);
    }
    // 4th from last is n_q1_.
    if (length > 3 + sign) {
      digits_.n_q1_ = get_digit_from_char(codestring[c++]);
      if (digits_.n_q1_ > 6) {
        throw InvalidPdgCode("Invalid PDG code " + codestring + " (n_q1>6)");
      }
    }
    // 3rd from last is n_q2_.
    if (length > 2 + sign) {
      digits_.n_q2_ = get_digit_from_char(codestring[c++]);
      if (digits_.n_q2_ > 6) {
        throw InvalidPdgCode("Invalid PDG code " + codestring + " (n_q2>6)");
      }
    }
    // Next to last is n_q3_.
    if (length > 1 + sign) {
      digits_.n_q3_ = get_digit_from_char(codestring[c++]);
      if (digits_.n_q3_ > 6) {
        throw InvalidPdgCode("Invalid PDG code " + codestring + " (n_q3>6)");
      }
    }
    // Last digit is the spin degeneracy.
    if (length > sign) {
      digits_.n_J_ += get_digit_from_char(codestring[c++]);
    } else {
      throw InvalidPdgCode(
          "String \"" + codestring +
          "\" only consists of a sign, that is no valid PDG Code\n");
    }
    check();
  }
 
  inline void set_fields(std::uint32_t abscode) {
    /* "dump_ =" overwrites antiparticle_, but this needs to have been set
     * already, so we carry it around the assignment. */
    bool ap = digits_.antiparticle_;
    dump_ = abscode & 0x0fffffff;
    digits_.antiparticle_ = ap;
    int test = test_code();
    if (test > 0) {
      throw InvalidPdgCode("Invalid digits " + std::to_string(test) +
                           " in PDG Code " + string());
    }
    check();
  }
};
 
static_assert(sizeof(PdgCode) == 4, "should fit into 32 bit integer");
 
std::istream& operator>>(std::istream& is, PdgCode& code);
std::ostream& operator<<(std::ostream& is, const PdgCode& code);
 
inline bool is_dilepton(const PdgCode pdg1, const PdgCode pdg2) {
  const auto c1 = pdg1.code();
  const auto c2 = pdg2.code();
  const auto min = std::min(c1, c2);
  const auto max = std::max(c1, c2);
  return (max == 0x11 && min == -0x11) || (max == 0x13 && min == -0x13);
}
 
inline bool has_lepton_pair(const PdgCode pdg1, const PdgCode pdg2,
                            const PdgCode pdg3) {
  return is_dilepton(pdg1, pdg2) || is_dilepton(pdg1, pdg3) ||
         is_dilepton(pdg2, pdg3);
}
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_PDGCODE_H_