kinematics.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2015-2021
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
 
#ifndef SRC_INCLUDE_SMASH_KINEMATICS_H_
#define SRC_INCLUDE_SMASH_KINEMATICS_H_
 
#include <array>
#include <cassert>
#include <sstream>
 
#include "constants.h"
 
namespace smash {
 
inline double center_of_velocity_v(double s, double ma, double mb) {
  const double m_sum = ma + mb;
  const double m_dif = ma - mb;
  return std::sqrt((s - m_sum * m_sum) / (s - m_dif * m_dif));
}
 
inline double fixed_target_projectile_v(double s, double ma, double mb) {
  const double inv_gamma = 2 * ma * mb / (s - ma * ma - mb * mb);
  return std::sqrt(1.0 - inv_gamma * inv_gamma);
}
 
template <typename T>
T pCM_sqr_from_s(const T s, const T mass_a, const T mass_b) noexcept {
  const auto mass_a_sqr = mass_a * mass_a;
  const auto x = s + mass_a_sqr - mass_b * mass_b;
  return x * x * (T(0.25) / s) - mass_a_sqr;
}
 
template <typename T>
T pCM_from_s(const T s, const T mass_a, const T mass_b) noexcept {
  const auto psqr = pCM_sqr_from_s(s, mass_a, mass_b);
  return psqr > T(0.) ? std::sqrt(psqr) : T(0.);
}
 
template <typename T>
T pCM(const T sqrts, const T mass_a, const T mass_b) noexcept {
  return pCM_from_s(sqrts * sqrts, mass_a, mass_b);
}
 
template <typename T>
T pCM_sqr(const T sqrts, const T mass_a, const T mass_b) noexcept {
  return pCM_sqr_from_s(sqrts * sqrts, mass_a, mass_b);
}
 
template <typename T>
std::array<T, 2> get_t_range(const T sqrts, const T m1, const T m2, const T m3,
                             const T m4) {
  const T p_i = pCM(sqrts, m1, m2);  // initial-state CM momentum
  const T p_f = pCM(sqrts, m3, m4);  // final-state CM momentum
  const T sqrt_t0 = (m1 * m1 - m2 * m2 - m3 * m3 + m4 * m4) / (2. * sqrts);
  const T t0 = sqrt_t0 * sqrt_t0;
  const T t_min = t0 - (p_i - p_f) * (p_i - p_f);
  const T t_max = t0 - (p_i + p_f) * (p_i + p_f);
  assert(t_min >= t_max);
  return {t_min, t_max};
}
 
static inline void check_energy(double mandelstam_s, double m_sum) {
  if (mandelstam_s < m_sum * m_sum) {
    std::stringstream err;
    err << "plab_from_s: s too small: " << mandelstam_s << " < "
        << m_sum * m_sum;
    throw std::runtime_error(err.str());
  }
}
 
static inline void check_radicand(double mandelstam_s, double radicand) {
  if (radicand < 0) {
    std::stringstream err;
    err << "plab_from_s: negative radicand: " << mandelstam_s;
    throw std::runtime_error(err.str());
  }
}
 
inline double plab_from_s(double mandelstam_s, double mass) {
  const double radicand = mandelstam_s * (mandelstam_s - 4 * mass * mass);
#ifndef NDEBUG
  const double m_sum = 2 * mass;
  check_energy(mandelstam_s, m_sum);
  check_radicand(mandelstam_s, radicand);
#endif
  return std::sqrt(radicand) / (2 * mass);
}
 
inline double plab_from_s(double mandelstam_s) {
  return plab_from_s(mandelstam_s, nucleon_mass);
}
 
inline double plab_from_s(double mandelstam_s, double m_projectile,
                          double m_target) {
  const double m_sum = m_projectile + m_target;
  const double m_diff = m_projectile - m_target;
  const double radicand =
      (mandelstam_s - m_sum * m_sum) * (mandelstam_s - m_diff * m_diff);
/* This is equivalent to:
 * const double radicand
 *     = (mandelstam_s - m_a_sq - m_b_sq) * (mandelstam_s - m_a_sq - m_b_sq)
 *     - 4 * m_a_sq * m_b_sq; */
#ifndef NDEBUG
  check_energy(mandelstam_s, m_sum);
  check_radicand(mandelstam_s, radicand);
#endif
  return std::sqrt(radicand) / (2 * m_target);
}
 
inline double s_from_Etot(double e_tot, double m_P, double m_T) {
  return m_P * m_P + m_T * m_T + 2 * m_T * e_tot;
}
inline double s_from_Etot(double e_tot_p, double e_tot_t, double m_p,
                          double m_t) {
  double pz_p = std::sqrt(e_tot_p * e_tot_p - m_p * m_p);
  double pz_t = std::sqrt(e_tot_t * e_tot_t - m_t * m_t);
  return std::pow(e_tot_p + e_tot_t, 2) - std::pow(pz_p - pz_t, 2);
}
inline double s_from_Ekin(double e_kin, double m_P, double m_T) {
  return s_from_Etot(e_kin + m_P, m_P, m_T);
}
inline double s_from_Ekin(double e_kin_p, double e_kin_t, double m_p,
                          double m_t) {
  return s_from_Etot(e_kin_p + m_t, e_kin_t + m_t, m_p, m_t);
}
inline double s_from_plab(double plab, double m_P, double m_T) {
  return m_P * m_P + m_T * m_T + 2 * m_T * std::sqrt(m_P * m_P + plab * plab);
}
inline double s_from_plab(double plab_p, double plab_t, double m_p,
                          double m_t) {
  return s_from_Etot(std::sqrt(m_p * m_p + plab_p * plab_p),
                     std::sqrt(m_t * m_t + plab_t * plab_t), plab_p, plab_t);
}
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_KINEMATICS_H_