kinematics.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2015-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */

#ifndef SRC_INCLUDE_KINEMATICS_H_
#define SRC_INCLUDE_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_Ekin(double e_kin, double m_P, double m_T) {
  return m_P * m_P + m_T * m_T + 2 * m_T * (m_P + e_kin);
}

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);
}

}  // namespace smash

#endif  // SRC_INCLUDE_KINEMATICS_H_