tabulation.h

Go to the documentation of this file.

/*
 *    Copyright (c) 2015-2020
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 */
 
#ifndef SRC_INCLUDE_SMASH_TABULATION_H_
#define SRC_INCLUDE_SMASH_TABULATION_H_
 
#include <fstream>
#include <functional>
#include <map>
#include <memory>
#include <vector>
 
#include "forwarddeclarations.h"
#include "integrate.h"
#include "kinematics.h"
#include "particletype.h"
#include "sha256.h"
 
namespace smash {
 
enum class Extrapolation {
  Zero = 0,
  Const = 1,
  Linear = 2,
};
 
class Tabulation {
 public:
  Tabulation() : values_({}), x_min_(0.0), x_max_(0.0), inv_dx_(0.0) {}
 
  Tabulation(double x_min, double range, size_t num,
             std::function<double(double)> f);
 
  bool is_empty() const { return values_.empty(); }
 
  static Tabulation from_file(std::ifstream& stream, sha256::Hash hash);
 
  double get_value_step(double x) const;
 
  double get_value_linear(
      double x, Extrapolation extrapolation = Extrapolation::Linear) const;
 
  void write(std::ofstream& stream, sha256::Hash hash) const;
 
 protected:
  std::vector<double> values_;
 
  double x_min_;
 
  double x_max_;
 
  double inv_dx_;
};
 
inline double spec_func_integrand_1res(double resonance_mass, double sqrts,
                                       double stable_mass,
                                       const ParticleType& type) {
  if (sqrts <= stable_mass + resonance_mass) {
    return 0.;
  }
 
  /* Integrand is the spectral function weighted by the CM momentum of the
   * final state. */
  return type.spectral_function(resonance_mass) *
         pCM(sqrts, stable_mass, resonance_mass);
}
 
inline double spec_func_integrand_2res(double sqrts, double res_mass_1,
                                       double res_mass_2,
                                       const ParticleType& t1,
                                       const ParticleType& t2) {
  if (sqrts <= res_mass_1 + res_mass_2) {
    return 0.;
  }
 
  /* Integrand is the product of the spectral function weighted by the
   * CM momentum of the final state. */
  return t1.spectral_function(res_mass_1) * t2.spectral_function(res_mass_2) *
         pCM(sqrts, res_mass_1, res_mass_2);
}
 
inline Tabulation spectral_integral_semistable(Integrator& integrate,
                                               const ParticleType& resonance,
                                               const ParticleType& stable,
                                               double range) {
  const double m_min = resonance.min_mass_kinematic();
  const double m_stable = stable.mass();
  return Tabulation(m_min + m_stable, range, 100, [&](double srts) {
    return integrate(m_min, srts - m_stable, [&](double m) {
      return spec_func_integrand_1res(m, srts, m_stable, resonance);
    });
  });
}
 
inline Tabulation spectral_integral_unstable(Integrator2d& integrate2d,
                                             const ParticleType& res1,
                                             const ParticleType& res2,
                                             double range) {
  const double m1_min = res1.min_mass_kinematic();
  const double m2_min = res2.min_mass_kinematic();
  return Tabulation(m1_min + m2_min, range, 125, [&](double srts) {
    const double m1_max = srts - m2_min;
    const double m2_max = srts - m1_min;
    return integrate2d(
        m1_min, m1_max, m2_min, m2_max, [&](double m1, double m2) {
          return spec_func_integrand_2res(srts, m1, m2, res1, res2);
        });
  });
}
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_TABULATION_H_