quantumnumbers.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_QUANTUMNUMBERS_H_
#define SRC_INCLUDE_SMASH_QUANTUMNUMBERS_H_
 
#include <string>
#include <vector>
 
#include "particles.h"
 
namespace smash {
 
class QuantumNumbers {
 public:
  QuantumNumbers()
      : momentum_(0., 0., 0., 0.),
        charge_(0),
        isospin3_(0),
        strangeness_(0),
        charmness_(0),
        bottomness_(0),
        baryon_number_(0) {}
 
  QuantumNumbers(const FourVector& m, const int q, const int i3, const int s,
                 const int c, const int b, const int B)
      : momentum_(m),
        charge_(q),
        isospin3_(i3),
        strangeness_(s),
        charmness_(c),
        bottomness_(b),
        baryon_number_(B) {}
 
  explicit QuantumNumbers(const Particles& particles) : QuantumNumbers() {
    for (const ParticleData& data : particles) {
      add_values(data);
    }
  }
 
  explicit QuantumNumbers(const std::vector<Particles>& ensembles)
      : QuantumNumbers() {
    for (const Particles& particles : ensembles) {
      for (const ParticleData& data : particles) {
        add_values(data);
      }
    }
  }
 
  explicit QuantumNumbers(const ParticleList& part) : QuantumNumbers() {
    for (const auto& p : part) {
      add_values(p);
    }
  }
 
  void add_values(const ParticleData& p) {
    momentum_ += p.momentum();
    charge_ += p.pdgcode().charge();
    isospin3_ += p.pdgcode().isospin3();
    strangeness_ += p.pdgcode().strangeness();
    charmness_ += p.pdgcode().charmness();
    bottomness_ += p.pdgcode().bottomness();
    baryon_number_ += p.pdgcode().baryon_number();
  }
 
  FourVector momentum() const { return momentum_; }
 
  int charge() const { return charge_; }
  int isospin3() const { return isospin3_; }
  int strangeness() const { return strangeness_; }
  int charmness() const { return charmness_; }
  int bottomness() const { return bottomness_; }
  int baryon_number() const { return baryon_number_; }
 
  bool operator==(const QuantumNumbers& rhs) const {
    // clang-format off
    // invariant mass of the current momentum
    const double m0_invariant = std::sqrt(std::fabs(momentum_.sqr()));
    // invariant mass of the momentum to be compared
    const double m1_invariant = std::sqrt(std::fabs(rhs.momentum_.sqr()));
    // average invariant mass
    const double avg_m_inv = 0.5 * (m0_invariant + m1_invariant);
    /* precision used to check relative error
     * 0.1 is multiplied to make the relative error check kicks in
     * for collision energies beyond 10 GeV. */
    const double precision_high = 0.1 * small_number;
    // difference in the four-momentum
    FourVector momentum_diff = momentum_ - rhs.momentum_;
    /* Check whether momentum is conserved.
     * Even if the absolute difference in each component is larger than
     * small_number, momentum is still considered to be conserved when
     * the deviation is small enough compared to the given invariant mass. */
    bool momentum_conserved = momentum_ == rhs.momentum_ ||
        (std::fabs(momentum_diff.x0()) < precision_high * avg_m_inv &&
         std::fabs(momentum_diff.x1()) < precision_high * avg_m_inv &&
         std::fabs(momentum_diff.x2()) < precision_high * avg_m_inv &&
         std::fabs(momentum_diff.x3()) < precision_high * avg_m_inv);
    return (momentum_conserved &&
            charge_ == rhs.charge_ &&
            isospin3_ == rhs.isospin3_ &&
            strangeness_ == rhs.strangeness_ &&
            charmness_ == rhs.charmness_ &&
            bottomness_ == rhs.bottomness_ &&
            baryon_number_ == rhs.baryon_number_);
    // clang-format on
  }
  bool operator!=(const QuantumNumbers& rhs) const { return !(*this == rhs); }
 
  QuantumNumbers operator-(const QuantumNumbers& rhs) const {
    return {momentum_ - rhs.momentum_,          charge_ - rhs.charge_,
            isospin3_ - rhs.isospin3_,          strangeness_ - rhs.strangeness_,
            charmness_ - rhs.charmness_,        bottomness_ - rhs.bottomness_,
            baryon_number_ - rhs.baryon_number_};
  }
 
  std::string report_deviations(const std::vector<Particles>& ensembles) const {
    QuantumNumbers current_values(ensembles);
    return report_deviations(current_values);
  }
 
  std::string report_deviations(const QuantumNumbers& rhs) const;
 
 private:
  FourVector momentum_;
 
  int charge_;
  int isospin3_;
  int strangeness_;
  int charmness_;
  int bottomness_;
  int baryon_number_;
};
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_QUANTUMNUMBERS_H_