tabulation.h

Go to the documentation of this file.

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

#ifndef SRC_INCLUDE_TABULATION_H_
#define SRC_INCLUDE_TABULATION_H_

#include <functional>
#include <map>
#include <memory>
#include <vector>

#include "forwarddeclarations.h"
#include "integrate.h"
#include "kinematics.h"
#include "particletype.h"

namespace smash {

enum class Extrapolation {
  Zero = 0,
  Const = 1,
  Linear = 2,
};

class Tabulation {
 public:
  Tabulation(double x_min, double range, int num,
             std::function<double(double)> f);

  double get_value_step(double x) const;

  double get_value_linear(
      double x, Extrapolation extrapolation = Extrapolation::Linear) const;

 protected:
  std::vector<double> values_;

  const double x_min_;

  const double x_max_;

  const 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 std::unique_ptr<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 make_unique<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 std::unique_ptr<Tabulation> spectral_integral_unstable(
    Integrator2dCuhre& 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 make_unique<Tabulation>(m1_min + m2_min, range, 100, [&](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_TABULATION_H_