parametrizations.cc

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2013-2026
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#include "smash/parametrizations.h"
 
#include <cmath>
#include <initializer_list>
#include <iostream>
#include <memory>
#include <set>
#include <vector>
 
#include "smash/average.h"
#include "smash/clebschgordan.h"
#include "smash/constants.h"
#include "smash/interpolation.h"
#include "smash/kinematics.h"
#include "smash/lowess.h"
#include "smash/parametrizations_data.h"
#include "smash/pow.h"
 
namespace smash {
 
bool parametrization_exists(const PdgCode& pdg_a, const PdgCode& pdg_b) {
  const bool two_nucleons = pdg_a.is_nucleon() && pdg_b.is_nucleon();
  const bool nucleon_and_kaon = (pdg_a.is_nucleon() && pdg_b.is_kaon()) ||
                                (pdg_a.is_kaon() && pdg_b.is_nucleon());
  const bool nucleon_and_pion = (pdg_a.is_nucleon() && pdg_b.is_pion()) ||
                                (pdg_a.is_pion() && pdg_b.is_nucleon());
  const bool two_pions = pdg_a.is_pion() && pdg_b.is_pion();
  return two_nucleons || nucleon_and_kaon || nucleon_and_pion || two_pions;
}
 
double xs_high_energy(double mandelstam_s, bool is_opposite_charge, double ma,
                      double mb, double P, double R1, double R2) {
  const double M = 2.1206;
  const double H = 0.272;
  const double eta1 = 0.4473;
  const double eta2 = 0.5486;
  const double s_sab = mandelstam_s / (ma + mb + M) / (ma + mb + M);
  double xs =
      H * std::log(s_sab) * std::log(s_sab) + P + R1 * std::pow(s_sab, -eta1);
  xs = is_opposite_charge ? xs + R2 * std::pow(s_sab, -eta2)
                          : xs - R2 * std::pow(s_sab, -eta2);
  return xs;
}
 
double pp_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, false, 0.939, 0.939, 34.41, 13.07, 7.394);
}
 
double ppbar_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, true, 0.939, 0.939, 34.41, 13.07, 7.394);
}
 
double np_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, false, 0.939, 0.939, 34.41, 12.52, 6.66);
}
 
double npbar_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, true, 0.939, 0.939, 34.41, 12.52, 6.66);
}
 
double piplusp_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, false, 0.939, 0.138, 18.75, 9.56, 1.767);
}
 
double piminusp_high_energy(double mandelstam_s) {
  return xs_high_energy(mandelstam_s, true, 0.939, 0.138, 18.75, 9.56, 1.767);
}
 
double xs_ppbar_annihilation(double mandelstam_s) {
  const double xs_ref = 120.;
  const double s_ref = 4. * nucleon_mass * nucleon_mass;
  const double constant_a = 0.05;
  const double constant_b = 0.6;
  const double factor = constant_a * constant_a * s_ref /
                            ((mandelstam_s - s_ref) * (mandelstam_s - s_ref) +
                             constant_a * constant_a * s_ref) +
                        constant_b;
  return xs_ref * (s_ref / mandelstam_s) * factor;
}
 
double xs_string_hard(double mandelstam_s, double xs_0, double e_0,
                      double lambda_pow) {
  const double sqrts = std::sqrt(mandelstam_s);
  if (sqrts < e_0) {
    return 0.;
  } else {
    double xs = xs_0 * std::pow(std::log(sqrts / e_0), lambda_pow);
    return xs;
  }
}
 
double NN_string_hard(double mandelstam_s) {
  return xs_string_hard(mandelstam_s, 0.087, 4.1, 3.8);
}
 
double Npi_string_hard(double mandelstam_s) {
  return xs_string_hard(mandelstam_s, 0.042, 3.5, 4.2);
}
 
double pipi_string_hard(double mandelstam_s) {
  return xs_string_hard(mandelstam_s, 0.013, 2.3, 4.7);
}
 
double pipluspiminus_total(double sqrts) {
  if (pipluspiminus_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIPLUSPIMINUS_TOT_SQRTS, PIPLUSPIMINUS_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.01, 10);
    pipluspiminus_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double last = *(PIPLUSPIMINUS_TOT_SQRTS.end() - 1);
  if (sqrts < last)
    return (*pipluspiminus_total_interpolation)(sqrts);
  else
    return pipi_string_hard(sqrts * sqrts);
}
 
double pizeropizero_total(double sqrts) {
  if (pizeropizero_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIZEROPIZERO_TOT_SQRTS, PIZEROPIZERO_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.01, 10);
    pizeropizero_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double last = *(PIZEROPIZERO_TOT_SQRTS.end() - 1);
  if (sqrts < last)
    return (*pizeropizero_total_interpolation)(sqrts);
  else
    return pipi_string_hard(sqrts * sqrts);
}
 
double piplusp_total(double sqrts) {
  if (piplusp_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIPLUSP_TOT_SQRTS, PIPLUSP_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.01, 10);
    piplusp_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double last = *(PIPLUSP_TOT_SQRTS.end() - 1);
  if (sqrts < last)
    return (*piplusp_total_interpolation)(sqrts);
  else
    return piplusp_high_energy(sqrts * sqrts);
}
 
/* pi+ p elastic cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
static double piplusp_elastic_pdg(double mandelstam_s) {
  if (piplusp_elastic_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIPLUSP_ELASTIC_P_LAB, PIPLUSP_ELASTIC_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.1, 5);
    piplusp_elastic_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  return (*piplusp_elastic_interpolation)(p_lab);
}
 
double piplusp_elastic_AQM(double mandelstam_s, double m1, double m2) {
  const double p_lab =
      plab_from_s_heavier_particle_at_rest(mandelstam_s, m1, m2);
  if (p_lab < 3.05) {  // the plab from which the param starts to explode
    return 7.5;        // this will be scaled down by 2/3 for meson-meson
  } else {
    const auto logp = std::log(p_lab);
    return 11.4 * std::pow(p_lab, -0.4) + 0.079 * logp * logp;
  }
}
 
double piplusp_elastic(double mandelstam_s) {
  double sigma;
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  if (mandelstam_s < 2.25) {
    sigma = really_small;
  } else if (mandelstam_s > 4.84) {
    const auto logp = std::log(p_lab);
    sigma = 11.4 * std::pow(p_lab, -0.4) + 0.079 * logp * logp;
  } else {
    sigma = piplusp_elastic_pdg(mandelstam_s);
  }
 
  // The elastic contributions from decays still need to be subtracted.
  if (piplusp_elastic_res_interpolation == nullptr) {
    std::vector<double> x = PIPLUSP_RES_SQRTS;
    for (auto& i : x) {
      i = i * i;
    }
    std::vector<double> y = PIPLUSP_RES_SIG;
    piplusp_elastic_res_interpolation = std::make_unique<InterpolateDataSpline>(
        x, y, ExtrapolationType::Constant);
  }
  sigma -= (*piplusp_elastic_res_interpolation)(mandelstam_s);
  if (sigma < 0) {
    sigma = really_small;
  }
  return sigma;
}
 
/* pi+ p to Sigma+ K+ cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
double piplusp_sigmapluskplus_pdg(double mandelstam_s) {
  if (piplusp_sigmapluskplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(PIPLUSP_SIGMAPLUSKPLUS_P_LAB,
                                                PIPLUSP_SIGMAPLUSKPLUS_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.2, 5);
    piplusp_sigmapluskplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  /* If p_lab is beyond the upper bound of the linear interpolation,
   * InterpolateDataLinear will return the value at the upper bound if
   * ExtrapolationType::Constant is set and this is what we want here. */
  return (*piplusp_sigmapluskplus_interpolation)(p_lab);
}
 
double piminusp_total(double sqrts) {
  if (piminusp_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIMINUSP_TOT_SQRTS, PIMINUSP_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.01, 6);
    piminusp_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double last = *(PIMINUSP_TOT_SQRTS.end() - 1);
  if (sqrts < last)
    return (*piminusp_total_interpolation)(sqrts);
  else
    return piminusp_high_energy(sqrts * sqrts);
}
 
/* pi- p elastic cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
static double piminusp_elastic_pdg(double mandelstam_s) {
  if (piminusp_elastic_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIMINUSP_ELASTIC_P_LAB, PIMINUSP_ELASTIC_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.2, 6);
    piminusp_elastic_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  return (*piminusp_elastic_interpolation)(p_lab);
}
 
double piminusp_elastic(double mandelstam_s) {
  double sigma;
  if (mandelstam_s < 1.69) {
    sigma = really_small;
  } else if (mandelstam_s > 4.84) {
    const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
    /* std::log(0) = -∞ and std::pow(0, negative) = +∞ both raise FE_DIVBYZERO.
     * Handle this case explicitly to preserve the mathematical result without
     * raising a floating-point exception. */
    if (p_lab == 0.0) {
      sigma = std::numeric_limits<double>::infinity();
    } else {
      const auto logp = std::log(p_lab);
      sigma = 1.76 + 11.2 * std::pow(p_lab, -0.64) + 0.043 * logp * logp;
    }
  } else {
    sigma = piminusp_elastic_pdg(mandelstam_s);
  }
  /* Tune down the elastic cross section when sqrt s is between 1.8 GeV
   * and 1.97 GeV so that the total cross section can fit the data. The
   * scaling factor is chosen so that the it's equal to one and its
   * derivate vanishes at the both ends. The minimum scaling factor in this
   * region is 0.88-0.12=0.76. */
  if (mandelstam_s > 3.24 && mandelstam_s < 3.8809) {
    sigma *= (0.12 * std::cos(2 * M_PI * (std::sqrt(mandelstam_s) - 1.8) /
                              (1.97 - 1.8)) +
              0.88);
  }
  // The elastic contributions from decays still need to be subtracted.
  if (piminusp_elastic_res_interpolation == nullptr) {
    std::vector<double> x = PIMINUSP_RES_SQRTS;
    for (auto& i : x) {
      i = i * i;
    }
    std::vector<double> y = PIMINUSP_RES_SIG;
    auto [dedup_x, dedup_y] = dedup_avg(x, y);
    piminusp_elastic_res_interpolation =
        std::make_unique<InterpolateDataSpline>(dedup_x, dedup_y,
                                                ExtrapolationType::Constant);
  }
  sigma -= (*piminusp_elastic_res_interpolation)(mandelstam_s);
  if (sigma < 0) {
    sigma = really_small;
  }
  return sigma;
}
 
/* pi- p -> Lambda K0 cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
double piminusp_lambdak0_pdg(double mandelstam_s) {
  if (piminusp_lambdak0_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(PIMINUSP_LAMBDAK0_P_LAB, PIMINUSP_LAMBDAK0_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.2, 6);
    piminusp_lambdak0_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  return (*piminusp_lambdak0_interpolation)(p_lab);
}
 
/* pi- p -> Sigma- K+ cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
double piminusp_sigmaminuskplus_pdg(double mandelstam_s) {
  if (piminusp_sigmaminuskplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(PIMINUSP_SIGMAMINUSKPLUS_P_LAB,
                                                PIMINUSP_SIGMAMINUSKPLUS_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.2, 6);
    piminusp_sigmaminuskplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, pion_mass, nucleon_mass);
  return (*piminusp_sigmaminuskplus_interpolation)(p_lab);
}
 
/* pi- p -> Sigma0 K0 cross section parametrization, resonance contribution.
 *
 * The data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one sqrts value, the corresponding cross sections
 * are averaged. */
double piminusp_sigma0k0_res(double mandelstam_s) {
  if (piminusp_sigma0k0_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(PIMINUSP_SIGMA0K0_RES_SQRTS,
                                                PIMINUSP_SIGMA0K0_RES_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.2, 6);
    piminusp_sigma0k0_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double sqrts = std::sqrt(mandelstam_s);
  return (*piminusp_sigma0k0_interpolation)(sqrts);
}
 
double pp_elastic(double mandelstam_s) {
  const double p_lab = plab_from_s(mandelstam_s);
  if (p_lab < 0.435) {
    return 5.12 * nucleon_mass /
               (mandelstam_s - 4 * nucleon_mass * nucleon_mass) +
           1.67;
  } else if (p_lab < 0.8) {
    return 23.5 + 1000 * pow_int(p_lab - 0.7, 4);
  } else if (p_lab < 2.0) {
    return 1250 / (p_lab + 50) - 4 * (p_lab - 1.3) * (p_lab - 1.3);
  } else if (p_lab < 2.776) {
    return 77 / (p_lab + 1.5);
  } else {
    const auto logp = std::log(p_lab);
    return 11.9 + 26.9 * std::pow(p_lab, -1.21) + 0.169 * logp * logp -
           1.85 * logp;
  }
}
 
double pp_elastic_high_energy(double mandelstam_s, double m1, double m2) {
  const double p_lab =
      plab_from_s_heavier_particle_at_rest(mandelstam_s, m1, m2);
  /* std::log(0) = -∞ and std::pow(0, negative) = +∞ both raise FE_DIVBYZERO.
   * Handle this case explicitly to preserve the mathematical result without
   * raising a floating-point exception. */
  if (p_lab == 0.0) {
    return std::numeric_limits<double>::infinity();
  } else {
    const auto logp = std::log(p_lab);
    return 11.9 + 26.9 * std::pow(p_lab, -1.21) + 0.169 * logp * logp -
           1.85 * logp;
  }
}
 
double pp_total(double mandelstam_s) {
  const double p_lab = plab_from_s(mandelstam_s);
  if (p_lab < 0.4) {
    return 34 * std::pow(p_lab / 0.4, -2.104);
  } else if (p_lab < 0.8) {
    return 23.5 + 1000 * pow_int(p_lab - 0.7, 4);
  } else if (p_lab < 1.5) {
    return 23.5 + 24.6 / (1 + std::exp(-(p_lab - 1.2) / 0.1));
  } else if (p_lab < 5.0) {
    return 41 + 60 * (p_lab - 0.9) * std::exp(-1.2 * p_lab);
  } else {
    const auto logp = std::log(p_lab);
    return 48.0 + 0.522 * logp * logp - 4.51 * logp;
  }
}
 
double np_elastic(double mandelstam_s) {
  const double p_lab = plab_from_s(mandelstam_s);
  if (p_lab < 0.525) {
    return 17.05 * nucleon_mass /
               (mandelstam_s - 4 * nucleon_mass * nucleon_mass) -
           6.83;
  } else if (p_lab < 0.8) {
    return 33 + 196 * std::pow(std::abs(p_lab - 0.95), 2.5);
  } else if (p_lab < 2.0) {
    return 31 / std::sqrt(p_lab);
  } else if (p_lab < 2.776) {
    return 77 / (p_lab + 1.5);
  } else {
    const auto logp = std::log(p_lab);
    return 11.9 + 26.9 * std::pow(p_lab, -1.21) + 0.169 * logp * logp -
           1.85 * logp;
  }
}
 
double np_total(double mandelstam_s) {
  const double p_lab = plab_from_s(mandelstam_s);
  const auto logp = std::log(p_lab);
  if (p_lab < 0.4) {
    return 6.3555 * std::pow(p_lab, -3.2481) * std::exp(-0.377 * logp * logp);
  } else if (p_lab < 1.0) {
    return 33 + 196 * std::pow(std::abs(p_lab - 0.95), 2.5);
  } else if (p_lab < 2.0) {
    return 24.2 + 8.9 * p_lab;
  } else if (p_lab < 5.0) {
    return 42;
  } else {
    return 48.0 + 0.522 * logp * logp - 4.51 * logp;
  }
}
 
double ppbar_elastic(double mandelstam_s) {
  if (mandelstam_s < 4 * nucleon_mass * nucleon_mass) {
    // Needed, since called directly from p_52
    return 0.0;
  }
  const double p_lab = plab_from_s(mandelstam_s);
  if (p_lab < 0.3) {
    return 78.6;
  } else if (p_lab < 5.0) {
    return 31.6 + 18.3 / p_lab - 1.1 / (p_lab * p_lab) - 3.8 * p_lab;
  } else {
    const auto logp = std::log(p_lab);
    return 10.2 + 52.7 * std::pow(p_lab, -1.16) + 0.125 * logp * logp -
           1.28 * logp;
  }
}
 
double ppbar_total(double mandelstam_s) {
  if (mandelstam_s < 4 * nucleon_mass * nucleon_mass) {
    // Needed, since called directly from p_52
    return 0.0;
  }
  const double p_lab = plab_from_s(mandelstam_s);
  if (p_lab < 0.3) {
    return 271.6 * std::exp(-1.1 * p_lab * p_lab);
  } else if (p_lab < 5.0) {
    return 75.0 + 43.1 / p_lab + 2.6 / (p_lab * p_lab) - 3.9 * p_lab;
  } else {
    const auto logp = std::log(p_lab);
    return 38.4 + 77.6 * std::pow(p_lab, -0.64) + 0.26 * logp * logp -
           1.2 * logp;
  }
}
 
double deuteron_pion_elastic(double mandelstam_s) {
  const double tmp = std::sqrt(mandelstam_s) - 2.172;
  return 4.0 + 0.27 / (tmp * tmp + 0.065 * 0.065);
}
 
double deuteron_nucleon_elastic(double mandelstam_s) {
  const double excess = (mandelstam_s - 7.93);
  const double excess_sqr = excess * excess;
  return 2500.0 * std::exp(-excess_sqr / 0.003) +
         600.0 * std::exp(-excess_sqr / 0.1) + 10.0;
}
 
double deuteron_pion_inelastic(double pion_kinetic_energy) {
  const double x = pion_kinetic_energy;
  return x * (4.3 + 10.0 * x) / ((x - 0.16) * (x - 0.16) + 0.007);
}
 
double deuteron_nucleon_inelastic(double N_kinetic_energy) {
  const double x = N_kinetic_energy;
  return x * (1.0 + 50 * x) / (x * x + 0.01) +
         4 * x / ((x - 0.008) * (x - 0.008) + 0.0004);
}
 
double deuteron_antinucleon_inelastic(double aN_kinetic_energy) {
  return 55.0 / (aN_kinetic_energy + 0.17);
}
 
double kplusp_total(double mandelstam_s) {
  if (kplusp_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KPLUSP_TOT_PLAB, KPLUSP_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.1, 5);
    kplusp_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kplusp_total_interpolation)(p_lab);
}
 
double kplusn_total(double mandelstam_s) {
  if (kplusn_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KPLUSN_TOT_PLAB, KPLUSN_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.05, 5);
    kplusn_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kplusn_total_interpolation)(p_lab);
}
 
double kminusp_total(double mandelstam_s) {
  if (kminusp_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KMINUSP_TOT_PLAB, KMINUSP_TOT_SIG);
    // Parametrization data KMINUSP_TOT_PLAB, KMINUSP_TOT_SIG is pre-smoothed
    dedup_y = smooth(dedup_x, dedup_y, 0.01, 5);
    kminusp_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kminusp_total_interpolation)(p_lab);
}
 
double kminusn_total(double mandelstam_s) {
  if (kminusn_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KMINUSN_TOT_PLAB, KMINUSN_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.05, 5);
    kminusn_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kminusn_total_interpolation)(p_lab);
}
 
double kplusp_elastic_background(double mandelstam_s) {
  constexpr double a0 = 10.508;  // mb
  constexpr double a1 = -3.716;  // mb/GeV
  constexpr double a2 = 1.845;   // mb/GeV^2
  constexpr double a3 = -0.764;  // GeV^-1
  constexpr double a4 = 0.508;   // GeV^-2
 
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  const double p_lab2 = p_lab * p_lab;
 
  return (a0 + a1 * p_lab + a2 * p_lab2) / (1 + a3 * p_lab + a4 * p_lab2);
}
 
double kplusn_elastic_background(double mandelstam_s) {
  return 0.25 * kplusp_elastic_background(mandelstam_s);
}
 
double kplusn_k0p(double mandelstam_s) {
  return 0.25 * kplusp_elastic_background(mandelstam_s);
}
 
/* K- p elastic cross section parametrization, PDG data.
 *
 * The PDG data is smoothed using the LOWESS algorithm. If more than one
 * cross section was given for one p_lab value, the corresponding cross sections
 * are averaged. */
static double kminusp_elastic_pdg(double mandelstam_s) {
  if (kminusp_elastic_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KMINUSP_ELASTIC_P_LAB, KMINUSP_ELASTIC_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.1, 5);
    kminusp_elastic_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kminusp_elastic_interpolation)(p_lab);
}
 
double kminusp_elastic_background(double mandelstam_s) {
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  double sigma;
  if (std::sqrt(mandelstam_s) < 1.68) {
    /* The parametrization here also works for anti-K0 n, Lambda pi0,
     * Sigma+ pi-, Sigma- pi+, Sigma0 pi0 with different parameters a0, a1, a2.
     *
     * The values of the parameters are *not* taken from the source above,
     * they come from a fit to PDG data. */
    constexpr double a0 = 186.03567644;  // mb GeV^2
    constexpr double a1 = 0.22002795;    // Gev
    constexpr double a2 = 0.64907116;
 
    /* In sigma a ratio p_i/p_f is omitted as both are set to p_lab. Keeping it
     * is unnecessary and it would trigger a FPE if p_lab were zero. */
    const double ratio = a1 * a1 / (a1 * a1 + p_lab * p_lab);
    sigma = a0 / mandelstam_s * std::pow(ratio, a2);
  } else {
    sigma = kminusp_elastic_pdg(mandelstam_s);
  }
  // The elastic contributions from decays still need to be subtracted.
  if (kminusp_elastic_res_interpolation == nullptr) {
    std::vector<double> x = KMINUSP_RES_SQRTS;
    for (auto& i : x) {
      i = plab_from_s(i * i, kaon_mass, nucleon_mass);
    }
    std::vector<double> y = KMINUSP_RES_SIG;
    kminusp_elastic_res_interpolation = std::make_unique<InterpolateDataSpline>(
        x, y, ExtrapolationType::Constant);
  }
  const auto old_sigma = sigma;
  sigma -= (*kminusp_elastic_res_interpolation)(p_lab);
  if (sigma < 0) {
    std::cout << "NEGATIVE SIGMA: sigma=" << sigma
              << ", sqrt(s)=" << std::sqrt(mandelstam_s)
              << ", sig_el_exp=" << old_sigma
              << ", sig_el_res=" << (*kminusp_elastic_res_interpolation)(p_lab)
              << std::endl;
  }
  assert(sigma >= 0);
  return sigma;
}
 
double kminusn_elastic_background(double) { return 4.0; }
 
double k0p_elastic_background(double mandelstam_s) {
  // by isospin symmetry
  return kplusn_elastic_background(mandelstam_s);
}
 
double k0n_elastic_background(double mandelstam_s) {
  // by isospin symmetry
  return kplusp_elastic_background(mandelstam_s);
}
 
double kbar0p_elastic_background(double mandelstam_s) {
  // by isospin symmetry
  return kminusn_elastic_background(mandelstam_s);
}
 
double kbar0n_elastic_background(double mandelstam_s) {
  // by isospin symmetry
  return kminusp_elastic_background(mandelstam_s);
}
 
double kplusp_inelastic_background(double mandelstam_s) {
  if (kplusp_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KPLUSP_TOT_PLAB, KPLUSP_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.1, 5);
    kplusp_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kplusp_total_interpolation)(p_lab)-kplusp_elastic_background(
      mandelstam_s);
}
 
double kplusn_inelastic_background(double mandelstam_s) {
  if (kplusn_total_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(KPLUSN_TOT_PLAB, KPLUSN_TOT_SIG);
    dedup_y = smooth(dedup_x, dedup_y, 0.05, 5);
    kplusn_total_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
  return (*kplusn_total_interpolation)(p_lab)-kplusn_elastic_background(
             mandelstam_s) -
         kplusn_k0p(mandelstam_s);
}
 
/**
 * Calculate and store isospin ratios for K N -> K Delta reactions.
 *
 * See the documentation of `KaonNucleonRatios` for details.
 *
 * \param[in] ratios An empty map where the ratios for K N -> K Delta
 *                      reactions are stored.
 */
static void initialize(std::unordered_map<std::pair<uint64_t, uint64_t>, double,
                                          pair_hash>& ratios) {
  const auto& type_p = ParticleType::find(pdg::p);
  const auto& type_n = ParticleType::find(pdg::n);
  const auto& type_K_p = ParticleType::find(pdg::K_p);
  const auto& type_K_z = ParticleType::find(pdg::K_z);
  const auto& type_Delta_pp = ParticleType::find(pdg::Delta_pp);
  const auto& type_Delta_p = ParticleType::find(pdg::Delta_p);
  const auto& type_Delta_z = ParticleType::find(pdg::Delta_z);
  const auto& type_Delta_m = ParticleType::find(pdg::Delta_m);
 
  /* Store the isospin ratio of the given reaction relative to all other
   * possible isospin-symmetric reactions. */
  auto add_to_ratios = [&](const ParticleType& a, const ParticleType& b,
                           const ParticleType& c, const ParticleType& d,
                           double weight_numerator, double weight_other) {
    assert(weight_numerator + weight_other != 0);
    const auto key =
        std::make_pair(pack(a.pdgcode().code(), b.pdgcode().code()),
                       pack(c.pdgcode().code(), d.pdgcode().code()));
    const double ratio = weight_numerator / (weight_numerator + weight_other);
    ratios[key] = ratio;
  };
 
  /* All inelastic channels are K N -> K Delta -> K pi N or charge exchange,
   * with identical cross section, weighted by the isospin factor.
   *
   * For charge exchange, the isospin factors are 1,
   * so they are excluded here. */
  {
    const auto weight1 = isospin_clebsch_gordan_sqr_2to2(
        type_p, type_K_p, type_K_z, type_Delta_pp);
    const auto weight2 = isospin_clebsch_gordan_sqr_2to2(
        type_p, type_K_p, type_K_p, type_Delta_p);
 
    add_to_ratios(type_p, type_K_p, type_K_z, type_Delta_pp, weight1, weight2);
    add_to_ratios(type_p, type_K_p, type_K_p, type_Delta_p, weight2, weight1);
  }
  {
    const auto weight1 = isospin_clebsch_gordan_sqr_2to2(
        type_n, type_K_p, type_K_z, type_Delta_p);
    const auto weight2 = isospin_clebsch_gordan_sqr_2to2(
        type_n, type_K_p, type_K_p, type_Delta_z);
 
    add_to_ratios(type_n, type_K_p, type_K_z, type_Delta_p, weight1, weight2);
    add_to_ratios(type_n, type_K_p, type_K_p, type_Delta_z, weight2, weight1);
  }
  /* K+ and K0 have the same mass and spin, their cross sections are assumed to
   * only differ in isospin factors. */
  {
    const auto weight1 = isospin_clebsch_gordan_sqr_2to2(
        type_p, type_K_z, type_K_z, type_Delta_p);
    const auto weight2 = isospin_clebsch_gordan_sqr_2to2(
        type_p, type_K_z, type_K_p, type_Delta_z);
 
    add_to_ratios(type_p, type_K_z, type_K_z, type_Delta_p, weight1, weight2);
    add_to_ratios(type_p, type_K_z, type_K_p, type_Delta_z, weight2, weight1);
  }
  {
    const auto weight1 = isospin_clebsch_gordan_sqr_2to2(
        type_n, type_K_z, type_K_z, type_Delta_z);
    const auto weight2 = isospin_clebsch_gordan_sqr_2to2(
        type_n, type_K_z, type_K_p, type_Delta_m);
 
    add_to_ratios(type_n, type_K_z, type_K_z, type_Delta_z, weight1, weight2);
    add_to_ratios(type_n, type_K_z, type_K_p, type_Delta_m, weight2, weight1);
  }
}
 
double KaonNucleonRatios::get_ratio(const ParticleType& a,
                                    const ParticleType& b,
                                    const ParticleType& c,
                                    const ParticleType& d) const {
  /* If this method is called with anti-nucleons, flip all particles to
   * anti-particles; the ratio is equal */
  int flip = 0;
  for (const auto& p : {&a, &b, &c, &d}) {
    if (p->is_nucleon()) {
      if (flip == 0) {
        flip = p->antiparticle_sign();
      } else {
        assert(p->antiparticle_sign() == flip);
      }
    }
  }
  const auto key = std::make_pair(
      pack(a.pdgcode().code() * flip, b.pdgcode().code() * flip),
      pack(c.pdgcode().code() * flip, d.pdgcode().code() * flip));
  if (ratios_.empty()) {
    initialize(ratios_);
  }
  return ratios_.at(key);
}
 
/*thread_local (see commit 897d0b8)*/ KaonNucleonRatios kaon_nucleon_ratios;
 
double kminusp_kbar0n(double mandelstam_s) {
  constexpr double a0 = 100;   // mb GeV^2
  constexpr double a1 = 0.15;  // GeV
  constexpr unsigned a2 = 2;
 
  const double p_lab = plab_from_s(mandelstam_s, kaon_mass, nucleon_mass);
 
  /* In this expression a ratio p_i/p_f is omitted as both are set to p_lab.
   * Keeping it is unnecessary and it would trigger a FPE if p_lab were zero. */
  return a0 / mandelstam_s * pow_int(a1 * a1 / (a1 * a1 + p_lab * p_lab), a2);
}
 
double kminusp_piminussigmaplus(double sqrts) {
  return 0.0788265 / ((sqrts - 1.38841) * (sqrts - 1.38841));
}
 
double kminusp_piplussigmaminus(double sqrts) {
  return 0.0196741 / ((sqrts - 1.42318) * (sqrts - 1.42318));
}
 
double kminusp_pi0sigma0(double sqrts) {
  return 0.0403364 / ((sqrts - 1.39830305) * (sqrts - 1.39830305));
}
 
double kminusp_pi0lambda(double sqrts) {
  return 0.05932562 / ((sqrts - 1.38786692) * (sqrts - 1.38786692));
}
 
double kminusn_piminussigma0(double sqrts) {
  return kminusp_piminussigmaplus(sqrts) + kminusp_piplussigmaminus(sqrts) -
         2. * kminusp_pi0sigma0(sqrts);
}
 
double kminusn_piminuslambda(double sqrts) {
  return 2. * kminusp_pi0lambda(sqrts);
}
 
// All K+ p and K+ n channels are forbidden by isospin.
 
double lambdalambda_ximinusp(double sqrts_sqrts0, double p_N, double p_lambda) {
  assert(p_lambda != 0);
  assert(sqrts_sqrts0 >= 0);
  return 37.15 / 2 * p_N / p_lambda * std::pow(sqrts_sqrts0, -0.16);
}
 
double lambdalambda_xi0n(double sqrts_sqrts0, double p_N, double p_lambda) {
  return lambdalambda_ximinusp(sqrts_sqrts0, p_N, p_lambda);
}
 
double lambdasigmaplus_xi0p(double sqrts_sqrts0) {
  assert(sqrts_sqrts0 >= 0);
  return 24.3781 * std::pow(sqrts_sqrts0, -0.479);
}
 
double lambdasigmaminus_ximinusn(double sqrts_sqrts0) {
  return lambdasigmaplus_xi0p(sqrts_sqrts0);
}
 
double lambdasigma0_ximinusp(double sqrts_sqrts0) {
  assert(sqrts_sqrts0 >= 0);
  if (sqrts_sqrts0 < 0.03336) {
    return 6.475 * std::pow(sqrts_sqrts0, -0.4167);
  } else {
    return 14.5054 * std::pow(sqrts_sqrts0, -0.1795);
  }
}
 
double lambdasigma0_xi0n(double sqrts_sqrts0) {
  return lambdasigma0_ximinusp(sqrts_sqrts0);
}
 
double sigma0sigma0_ximinusp(double sqrts_sqrts0) {
  assert(sqrts_sqrts0 >= 0);
  if (sqrts_sqrts0 < 0.09047) {
    return 5.625 * std::pow(sqrts_sqrts0, -0.318);
  } else {
    return 4.174 * std::pow(sqrts_sqrts0, -0.4421);
  }
}
 
double sigma0sigma0_xi0n(double sqrts_sqrts0) {
  return sigma0sigma0_ximinusp(sqrts_sqrts0);
}
 
double sigmaplussigmaminus_xi0p(double sqrts_sqrts0) {
  return 4 * sigma0sigma0_ximinusp(sqrts_sqrts0);
}
 
double sigma0sigmaminus_ximinusn(double sqrts_sqrts0) {
  return 4 * sigma0sigma0_ximinusp(sqrts_sqrts0);
}
 
double sigmaplussigmaminus_ximinusp(double sqrts_sqrts0) {
  return 14.194 * std::pow(sqrts_sqrts0, -0.442);
}
 
double sigmaplussigmaminus_xi0n(double sqrts_sqrts0) {
  return sigmaplussigmaminus_ximinusp(sqrts_sqrts0);
}
 
std::optional<double> Dzeropiplus_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzeropiplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DZEROPIPLUS_ELASTIC_SIG);
      Dzeropiplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzeropiplus_elastic_interpolation)(sqrts);
  }
}
 
double Dzeropiplus_Dpluspizero(double sqrts) {
  if (Dzeropiplus_Dpluspizero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DPI_SQRTS, DZEROPIPLUS_DPLUSPIZERO_SIG);
    Dzeropiplus_Dpluspizero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dzeropiplus_Dpluspizero_interpolation)(sqrts);
}
 
std::optional<double> Dzeropiminus_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzeropiminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DZEROPIMINUS_ELASTIC_SIG);
      Dzeropiminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzeropiminus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dzeropizero_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzeropizero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DZEROPIZERO_ELASTIC_SIG);
      Dzeropizero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzeropizero_elastic_interpolation)(sqrts);
  }
}
 
double Dzeropizero_Dpluspiminus(double sqrts) {
  if (Dzeropizero_Dpluspiminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DPI_SQRTS, DZEROPIZERO_DPLUSPIMINUS_SIG);
    Dzeropizero_Dpluspiminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dzeropizero_Dpluspiminus_interpolation)(sqrts);
}
 
std::optional<double> Dpluspiplus_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dpluspiplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DPLUSPIPLUS_ELASTIC_SIG);
      Dpluspiplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dpluspiplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dpluspiminus_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dpluspiminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DPLUSPIMINUS_ELASTIC_SIG);
      Dpluspiminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dpluspiminus_elastic_interpolation)(sqrts);
  }
}
 
double Dpluspiminus_Dzeropizero(double sqrts) {
  if (Dpluspiminus_Dzeropizero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DPI_SQRTS, DPLUSPIMINUS_DZEROPIZERO_SIG);
    Dpluspiminus_Dzeropizero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dpluspiminus_Dzeropizero_interpolation)(sqrts);
}
 
std::optional<double> Dpluspizero_elastic(double sqrts) {
  if (sqrts > *(DPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dpluspizero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DPI_SQRTS, DPLUSPIZERO_ELASTIC_SIG);
      Dpluspizero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dpluspizero_elastic_interpolation)(sqrts);
  }
}
 
double Dpluspizero_Dzeropiplus(double sqrts) {
  if (Dpluspizero_Dzeropiplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DPI_SQRTS, DPLUSPIZERO_DZEROPIPLUS_SIG);
    Dpluspizero_Dzeropiplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dpluspizero_Dzeropiplus_interpolation)(sqrts);
}
 
std::optional<double> Dpluseta_elastic(double sqrts) {
  if (sqrts > *(DETA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dpluseta_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DETA_SQRTS, DPLUSETA_ELASTIC_SIG);
      Dpluseta_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dpluseta_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dzeroeta_elastic(double sqrts) {
  if (sqrts > *(DETA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzeroeta_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DETA_SQRTS, DZEROETA_ELASTIC_SIG);
      Dzeroeta_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzeroeta_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusKplus_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusKplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DPLUSKPLUS_ELASTIC_SIG);
      DplusKplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusKplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusKzero_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusKzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DPLUSKZERO_ELASTIC_SIG);
      DplusKzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusKzero_elastic_interpolation)(sqrts);
  }
}
 
double DplusKzero_DzeroKplus(double sqrts) {
  if (DplusKzero_DzeroKplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DKAON_SQRTS, DPLUSKZERO_DZEROKPLUS_SIG);
    DplusKzero_DzeroKplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DplusKzero_DzeroKplus_interpolation)(sqrts);
}
 
double DzeroKplus_DplusKzero(double sqrts) {
  if (DzeroKplus_DplusKzero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DKAON_SQRTS, DZEROKPLUS_DPLUSKZERO_SIG);
    DzeroKplus_DplusKzero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DzeroKplus_DplusKzero_interpolation)(sqrts);
}
 
std::optional<double> DzeroKplus_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroKplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DZEROKPLUS_ELASTIC_SIG);
      DzeroKplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroKplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DzeroKzero_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroKzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DZEROKZERO_ELASTIC_SIG);
      DzeroKzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroKzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusKbarzero_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusKbarzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DPLUSKBARZERO_ELASTIC_SIG);
      DplusKbarzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusKbarzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusKminus_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusKminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DPLUSKMINUS_ELASTIC_SIG);
      DplusKminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusKminus_elastic_interpolation)(sqrts);
  }
}
 
double DplusKminus_DzeroKbarzero(double sqrts) {
  if (DplusKminus_DzeroKbarzero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DKAON_SQRTS, DPLUSKMINUS_DZEROKBARZERO_SIG);
    DplusKminus_DzeroKbarzero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DplusKminus_DzeroKbarzero_interpolation)(sqrts);
}
 
double DzeroKbarzero_DplusKminus(double sqrts) {
  if (DzeroKbarzero_DplusKminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DKAON_SQRTS, DZEROKBARZERO_DPLUSKMINUS_SIG);
    DzeroKbarzero_DplusKminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DzeroKbarzero_DplusKminus_interpolation)(sqrts);
}
 
std::optional<double> DzeroKbarzero_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroKbarzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DZEROKBARZERO_ELASTIC_SIG);
      DzeroKbarzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroKbarzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DzeroKminus_elastic(double sqrts) {
  if (sqrts > *(DKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroKminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DKAON_SQRTS, DZEROKMINUS_ELASTIC_SIG);
      DzeroKminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroKminus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dstarpluspiplus_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarpluspiplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARPLUSPIPLUS_ELASTIC_SIG);
      Dstarpluspiplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarpluspiplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dstarpluspiminus_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarpluspiminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARPLUSPIMINUS_ELASTIC_SIG);
      Dstarpluspiminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarpluspiminus_elastic_interpolation)(sqrts);
  }
}
 
double Dstarpluspiminus_Dstarzeropizero(double sqrts) {
  if (Dstarpluspiminus_Dstarzeropizero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARPI_SQRTS, DSTARPLUSPIMINUS_DSTARZEROPIZERO_SIG);
    Dstarpluspiminus_Dstarzeropizero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dstarpluspiminus_Dstarzeropizero_interpolation)(sqrts);
}
 
std::optional<double> Dstarpluspizero_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarpluspizero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARPLUSPIZERO_ELASTIC_SIG);
      Dstarpluspizero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarpluspizero_elastic_interpolation)(sqrts);
  }
}
 
double Dstarpluspizero_Dstarzeropiplus(double sqrts) {
  if (Dstarpluspizero_Dstarzeropiplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARPI_SQRTS, DSTARPLUSPIZERO_DSTARZEROPIPLUS_SIG);
    Dstarpluspizero_Dstarzeropiplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dstarpluspizero_Dstarzeropiplus_interpolation)(sqrts);
}
 
double Dstarzeropiplus_Dstarpluspizero(double sqrts) {
  if (Dstarzeropiplus_Dstarpluspizero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARPI_SQRTS, DSTARZEROPIPLUS_DSTARPLUSPIZERO_SIG);
    Dstarzeropiplus_Dstarpluspizero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dstarzeropiplus_Dstarpluspizero_interpolation)(sqrts);
}
 
std::optional<double> Dstarzeropiplus_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarzeropiplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARZEROPIPLUS_ELASTIC_SIG);
      Dstarzeropiplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarzeropiplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dstarzeropiminus_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarzeropiminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARZEROPIMINUS_ELASTIC_SIG);
      Dstarzeropiminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarzeropiminus_elastic_interpolation)(sqrts);
  }
}
 
double Dstarzeropizero_Dstarpluspiminus(double sqrts) {
  if (Dstarzeropizero_Dstarpluspiminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARPI_SQRTS, DSTARZEROPIZERO_DSTARPLUSPIMINUS_SIG);
    Dstarzeropizero_Dstarpluspiminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dstarzeropizero_Dstarpluspiminus_interpolation)(sqrts);
}
 
std::optional<double> Dstarzeropizero_elastic(double sqrts) {
  if (sqrts > *(DSTARPI_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarzeropizero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARPI_SQRTS, DSTARZEROPIZERO_ELASTIC_SIG);
      Dstarzeropizero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarzeropizero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dstarpluseta_elastic(double sqrts) {
  if (sqrts > *(DSTARETA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarpluseta_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARETA_SQRTS, DSTARPLUSETA_ELASTIC_SIG);
      Dstarpluseta_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarpluseta_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dstarzeroeta_elastic(double sqrts) {
  if (sqrts > *(DSTARETA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dstarzeroeta_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARETA_SQRTS, DSTARZEROETA_ELASTIC_SIG);
      Dstarzeroeta_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dstarzeroeta_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarplusKplus_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarplusKplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARPLUSKPLUS_ELASTIC_SIG);
      DstarplusKplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarplusKplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarplusKzero_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarplusKzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARPLUSKZERO_ELASTIC_SIG);
      DstarplusKzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarplusKzero_elastic_interpolation)(sqrts);
  }
}
 
double DstarplusKzero_DstarzeroKplus(double sqrts) {
  if (DstarplusKzero_DstarzeroKplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARKAON_SQRTS, DSTARPLUSKZERO_DSTARZEROKPLUS_SIG);
    DstarplusKzero_DstarzeroKplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DstarplusKzero_DstarzeroKplus_interpolation)(sqrts);
}
 
double DstarzeroKplus_DstarplusKzero(double sqrts) {
  if (DstarzeroKplus_DstarplusKzero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DSTARKAON_SQRTS, DSTARZEROKPLUS_DSTARPLUSKZERO_SIG);
    DstarzeroKplus_DstarplusKzero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DstarzeroKplus_DstarplusKzero_interpolation)(sqrts);
}
 
std::optional<double> DstarzeroKplus_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarzeroKplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARZEROKPLUS_ELASTIC_SIG);
      DstarzeroKplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarzeroKplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarzeroKzero_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarzeroKzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARZEROKZERO_ELASTIC_SIG);
      DstarzeroKzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarzeroKzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarplusKbarzero_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarplusKbarzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARPLUSKBARZERO_ELASTIC_SIG);
      DstarplusKbarzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarplusKbarzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarplusKminus_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarplusKminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARPLUSKMINUS_ELASTIC_SIG);
      DstarplusKminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarplusKminus_elastic_interpolation)(sqrts);
  }
}
 
double DstarplusKminus_DstarzeroKbarzero(double sqrts) {
  if (DstarplusKminus_DstarzeroKbarzero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DSTARKAON_SQRTS, DSTARPLUSKMINUS_DSTARZEROKBARZERO_SIG);
    DstarplusKminus_DstarzeroKbarzero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DstarplusKminus_DstarzeroKbarzero_interpolation)(sqrts);
}
 
double DstarzeroKbarzero_DstarplusKminus(double sqrts) {
  if (DstarzeroKbarzero_DstarplusKminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DSTARKAON_SQRTS, DSTARZEROKBARZERO_DSTARPLUSKMINUS_SIG);
    DstarzeroKbarzero_DstarplusKminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DstarzeroKbarzero_DstarplusKminus_interpolation)(sqrts);
}
 
std::optional<double> DstarzeroKbarzero_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarzeroKbarzero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARZEROKBARZERO_ELASTIC_SIG);
      DstarzeroKbarzero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarzeroKbarzero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DstarzeroKminus_elastic(double sqrts) {
  if (sqrts > *(DSTARKAON_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DstarzeroKminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DSTARKAON_SQRTS, DSTARZEROKMINUS_ELASTIC_SIG);
      DstarzeroKminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DstarzeroKminus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dplusn_elastic(double sqrts) {
  if (sqrts > *(DN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dplusn_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DPLUSN_ELASTIC_SIG);
      Dplusn_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dplusn_elastic_interpolation)(sqrts);
  }
}
 
double Dplusn_Dzerop(double sqrts) {
  if (Dplusn_Dzerop_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DPLUSN_DZEROP_SIG);
    Dplusn_Dzerop_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dplusn_Dzerop_interpolation)(sqrts);
}
 
std::optional<double> Dplusp_elastic(double sqrts) {
  if (sqrts > *(DN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dplusp_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DPLUSP_ELASTIC_SIG);
      Dplusp_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dplusp_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dzeron_elastic(double sqrts) {
  if (sqrts > *(DN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzeron_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DZERON_ELASTIC_SIG);
      Dzeron_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzeron_elastic_interpolation)(sqrts);
  }
}
 
double Dzerop_Dplusn(double sqrts) {
  if (Dzerop_Dplusn_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DZEROP_DPLUSN_SIG);
    Dzerop_Dplusn_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dzerop_Dplusn_interpolation)(sqrts);
}
 
std::optional<double> Dzerop_elastic(double sqrts) {
  if (sqrts > *(DN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dzerop_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] = dedup_avg<double>(DN_SQRTS, DZEROP_ELASTIC_SIG);
      Dzerop_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dzerop_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dminusn_elastic(double sqrts) {
  if (sqrts > *(DBARN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dminusn_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARN_SQRTS, DMINUSN_ELASTIC_SIG);
      Dminusn_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dminusn_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dminusp_elastic(double sqrts) {
  if (sqrts > *(DBARN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dminusp_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARN_SQRTS, DMINUSP_ELASTIC_SIG);
      Dminusp_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dminusp_elastic_interpolation)(sqrts);
  }
}
 
double Dminusp_Dbarzeron(double sqrts) {
  if (Dminusp_Dbarzeron_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DBARN_SQRTS, DMINUSP_DBARZERON_SIG);
    Dminusp_Dbarzeron_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dminusp_Dbarzeron_interpolation)(sqrts);
}
 
double Dbarzeron_Dminusp(double sqrts) {
  if (Dbarzeron_Dminusp_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DBARN_SQRTS, DBARZERON_DMINUSP_SIG);
    Dbarzeron_Dminusp_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*Dbarzeron_Dminusp_interpolation)(sqrts);
}
 
std::optional<double> Dbarzeron_elastic(double sqrts) {
  if (sqrts > *(DBARN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dbarzeron_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARN_SQRTS, DBARZERON_ELASTIC_SIG);
      Dbarzeron_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dbarzeron_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> Dbarzerop_elastic(double sqrts) {
  if (sqrts > *(DBARN_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (Dbarzerop_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARN_SQRTS, DBARZEROP_ELASTIC_SIG);
      Dbarzerop_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*Dbarzerop_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusDeltaplus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusDeltaplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAPLUS_ELASTIC_SIG);
      DplusDeltaplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusDeltaplus_elastic_interpolation)(sqrts);
  }
}
 
double DplusDeltaplus_DzeroDeltaplusplus(double sqrts) {
  if (DplusDeltaplus_DzeroDeltaplusplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAPLUS_DZERODELTAPLUSPLUS_SIG);
    DplusDeltaplus_DzeroDeltaplusplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DplusDeltaplus_DzeroDeltaplusplus_interpolation)(sqrts);
}
 
std::optional<double> DplusDeltaplusplus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusDeltaplusplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAPLUSPLUS_ELASTIC_SIG);
      DplusDeltaplusplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusDeltaplusplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DplusDeltaminus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusDeltaminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAMINUS_ELASTIC_SIG);
      DplusDeltaminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusDeltaminus_elastic_interpolation)(sqrts);
  }
}
 
double DplusDeltaminus_DzeroDeltazero(double sqrts) {
  if (DplusDeltaminus_DzeroDeltazero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAMINUS_DZERODELTAZERO_SIG);
    DplusDeltaminus_DzeroDeltazero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DplusDeltaminus_DzeroDeltazero_interpolation)(sqrts);
}
 
std::optional<double> DplusDeltazero_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DplusDeltazero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAZERO_ELASTIC_SIG);
      DplusDeltazero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DplusDeltazero_elastic_interpolation)(sqrts);
  }
}
 
double DplusDeltazero_DzeroDeltaplus(double sqrts) {
  if (DplusDeltazero_DzeroDeltaplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DPLUSDELTAZERO_DZERODELTAPLUS_SIG);
    DplusDeltazero_DzeroDeltaplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DplusDeltazero_DzeroDeltaplus_interpolation)(sqrts);
}
 
double DzeroDeltaplus_DplusDeltazero(double sqrts) {
  if (DzeroDeltaplus_DplusDeltazero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DZERODELTAPLUS_DPLUSDELTAZERO_SIG);
    DzeroDeltaplus_DplusDeltazero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DzeroDeltaplus_DplusDeltazero_interpolation)(sqrts);
}
 
std::optional<double> DzeroDeltaplus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroDeltaplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DZERODELTAPLUS_ELASTIC_SIG);
      DzeroDeltaplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroDeltaplus_elastic_interpolation)(sqrts);
  }
}
 
double DzeroDeltaplusplus_DplusDeltaplus(double sqrts) {
  if (DzeroDeltaplusplus_DplusDeltaplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DZERODELTAPLUSPLUS_DPLUSDELTAPLUS_SIG);
    DzeroDeltaplusplus_DplusDeltaplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DzeroDeltaplusplus_DplusDeltaplus_interpolation)(sqrts);
}
 
std::optional<double> DzeroDeltaplusplus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroDeltaplusplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DZERODELTAPLUSPLUS_ELASTIC_SIG);
      DzeroDeltaplusplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroDeltaplusplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DzeroDeltaminus_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroDeltaminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DZERODELTAMINUS_ELASTIC_SIG);
      DzeroDeltaminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroDeltaminus_elastic_interpolation)(sqrts);
  }
}
 
double DzeroDeltazero_DplusDeltaminus(double sqrts) {
  if (DzeroDeltazero_DplusDeltaminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] =
        dedup_avg<double>(DDELTA_SQRTS, DZERODELTAZERO_DPLUSDELTAMINUS_SIG);
    DzeroDeltazero_DplusDeltaminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DzeroDeltazero_DplusDeltaminus_interpolation)(sqrts);
}
 
std::optional<double> DzeroDeltazero_elastic(double sqrts) {
  if (sqrts > *(DDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DzeroDeltazero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DDELTA_SQRTS, DZERODELTAZERO_ELASTIC_SIG);
      DzeroDeltazero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DzeroDeltazero_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DminusDeltaplus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DminusDeltaplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DMINUSDELTAPLUS_ELASTIC_SIG);
      DminusDeltaplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DminusDeltaplus_elastic_interpolation)(sqrts);
  }
}
 
double DminusDeltaplus_DbarzeroDeltazero(double sqrts) {
  if (DminusDeltaplus_DbarzeroDeltazero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DMINUSDELTAPLUS_DBARZERODELTAZERO_SIG);
    DminusDeltaplus_DbarzeroDeltazero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DminusDeltaplus_DbarzeroDeltazero_interpolation)(sqrts);
}
 
std::optional<double> DminusDeltaplusplus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DminusDeltaplusplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DMINUSDELTAPLUSPLUS_ELASTIC_SIG);
      DminusDeltaplusplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DminusDeltaplusplus_elastic_interpolation)(sqrts);
  }
}
 
double DminusDeltaplusplus_DbarzeroDeltaplus(double sqrts) {
  if (DminusDeltaplusplus_DbarzeroDeltaplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DMINUSDELTAPLUSPLUS_DBARZERODELTAPLUS_SIG);
    DminusDeltaplusplus_DbarzeroDeltaplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DminusDeltaplusplus_DbarzeroDeltaplus_interpolation)(sqrts);
}
 
std::optional<double> DminusDeltaminus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DminusDeltaminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DMINUSDELTAMINUS_ELASTIC_SIG);
      DminusDeltaminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DminusDeltaminus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DminusDeltazero_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DminusDeltazero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DMINUSDELTAZERO_ELASTIC_SIG);
      DminusDeltazero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DminusDeltazero_elastic_interpolation)(sqrts);
  }
}
 
double DminusDeltazero_DbarzeroDeltaminus(double sqrts) {
  if (DminusDeltazero_DbarzeroDeltaminus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DMINUSDELTAZERO_DBARZERODELTAMINUS_SIG);
    DminusDeltazero_DbarzeroDeltaminus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DminusDeltazero_DbarzeroDeltaminus_interpolation)(sqrts);
}
 
double DbarzeroDeltaplus_DminusDeltaplusplus(double sqrts) {
  if (DbarzeroDeltaplus_DminusDeltaplusplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DBARZERODELTAPLUS_DMINUSDELTAPLUSPLUS_SIG);
    DbarzeroDeltaplus_DminusDeltaplusplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DbarzeroDeltaplus_DminusDeltaplusplus_interpolation)(sqrts);
}
 
std::optional<double> DbarzeroDeltaplus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DbarzeroDeltaplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DBARZERODELTAPLUS_ELASTIC_SIG);
      DbarzeroDeltaplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DbarzeroDeltaplus_elastic_interpolation)(sqrts);
  }
}
 
std::optional<double> DbarzeroDeltaplusplus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DbarzeroDeltaplusplus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DBARZERODELTAPLUSPLUS_ELASTIC_SIG);
      DbarzeroDeltaplusplus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DbarzeroDeltaplusplus_elastic_interpolation)(sqrts);
  }
}
 
double DbarzeroDeltaminus_DminusDeltazero(double sqrts) {
  if (DbarzeroDeltaminus_DminusDeltazero_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DBARZERODELTAMINUS_DMINUSDELTAZERO_SIG);
    DbarzeroDeltaminus_DminusDeltazero_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DbarzeroDeltaminus_DminusDeltazero_interpolation)(sqrts);
}
 
std::optional<double> DbarzeroDeltaminus_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DbarzeroDeltaminus_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DBARZERODELTAMINUS_ELASTIC_SIG);
      DbarzeroDeltaminus_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DbarzeroDeltaminus_elastic_interpolation)(sqrts);
  }
}
 
double DbarzeroDeltazero_DminusDeltaplus(double sqrts) {
  if (DbarzeroDeltazero_DminusDeltaplus_interpolation == nullptr) {
    auto [dedup_x, dedup_y] = dedup_avg<double>(
        DBARDELTA_SQRTS, DBARZERODELTAZERO_DMINUSDELTAPLUS_SIG);
    DbarzeroDeltazero_DminusDeltaplus_interpolation =
        std::make_unique<InterpolateDataLinear<double>>(
            dedup_x, dedup_y, ExtrapolationType::Constant);
  }
  return (*DbarzeroDeltazero_DminusDeltaplus_interpolation)(sqrts);
}
 
std::optional<double> DbarzeroDeltazero_elastic(double sqrts) {
  if (sqrts > *(DBARDELTA_SQRTS.end() - 1)) {
    return std::nullopt;
  } else {
    if (DbarzeroDeltazero_elastic_interpolation == nullptr) {
      auto [dedup_x, dedup_y] =
          dedup_avg<double>(DBARDELTA_SQRTS, DBARZERODELTAZERO_ELASTIC_SIG);
      DbarzeroDeltazero_elastic_interpolation =
          std::make_unique<InterpolateDataLinear<double>>(
              dedup_x, dedup_y, ExtrapolationType::Zero);
    }
    return (*DbarzeroDeltazero_elastic_interpolation)(sqrts);
  }
}
 
}  // namespace smash