src/sfset/Pavia.cxx

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Duane Byer
0003
0004 #include "Pavia.h"
0005
0006 #include "../interp/Interpolate.h"
0007
0008 #include <fstream>
0009 #include <stdexcept>
0010 #include <sstream>
0011 #include <string>
0012
0013 #define SF_SET_DIR "grids"
0014 #define PAVIA_DIR "pavia"
0015
0016 namespace {
0017
0018 struct DataFileNotFoundError : public std::runtime_error {
0019     DataFileNotFoundError(char const* name) :
0020         std::runtime_error(std::string("Couldn't load data file '") + name + "'") {
0021     }
0022 };
0023 struct DataFileParseError : public std::runtime_error {
0024     DataFileParseError(char const* name) :
0025         std::runtime_error(std::string("Couldn't parse data file '") + name + "'") {
0026     }
0027 };
0028
0029 std::istream& find_file(std::ifstream& fin, char const* file_name) {
0030     fin.open(std::string("../share/" SF_SET_DIR "/" PAVIA_DIR "/") + file_name);
0031     if (fin) {
0032         return fin;
0033     }
0034     fin.open(std::string(SF_SET_DIR "/" PAVIA_DIR "/") + file_name);
0035     if (fin) {
0036         return fin;
0037     }
0038     fin.open(std::string(PAVIA_DIR "/") + file_name);
0039     if (fin) {
0040         return fin;
0041     }
0042     fin.open(file_name);
0043     if (fin) {
0044         return fin;
0045     }
0046     throw DataFileNotFoundError(file_name);
0047 }
0048
0049 std::array<Grid<double, 4>, 2> load_grids(char const* file_name) {
0050     std::ifstream in;
0051     find_file(in, file_name);
0052     std::vector<std::array<double, 6> > data;
0053     std::string line;
0054     std::getline(in, line);
0055     while (std::getline(in, line)) {
0056         std::string blank;
0057         std::stringstream line_in(line);
0058         std::array<double, 6> next;
0059         line_in >> blank >> blank;
0060         line_in >> next[0] >> next[1] >> next[2] >> next[3];
0061         line_in >> blank >> blank >> blank >> blank >> blank;
0062         line_in >> next[4] >> next[5];
0063         data.push_back(next);
0064         if (!line_in) {
0065             throw DataFileParseError(file_name);
0066         }
0067     }
0068     // By default, assume the grids are only accurate to single precision.
0069     return read_grids<double, 4, 2>(data);
0070 }
0071
0072 }
0073
0074 struct PaviaSfSet::Impl {
0075     std::array<Grid<double, 4>, 2> grid_pip;
0076     std::array<Grid<double, 4>, 2> grid_pim;
0077     CubicView<double, 4> interp_pip_uu;
0078     CubicView<double, 4> interp_pip_ut_sivers;
0079     CubicView<double, 4> interp_pim_uu;
0080     CubicView<double, 4> interp_pim_ut_sivers;
0081     Impl() :
0082         grid_pip(load_grids("grid_Pip_2.txt")),
0083         grid_pim(load_grids("grid_Pim_2.txt")),
0084         interp_pip_uu(grid_pip[0]),
0085         interp_pip_ut_sivers(grid_pip[1]),
0086         interp_pim_uu(grid_pim[0]),
0087         interp_pim_ut_sivers(grid_pim[1]) {
0088     }
0089 };
0090
0091 PaviaSfSet::PaviaSfSet(PaviaSfSet&& other) noexcept :
0092         _impl(nullptr) {
0093     std::swap(_impl, other._impl);
0094 }
0095 PaviaSfSet& PaviaSfSet::operator=(PaviaSfSet&& other) noexcept {
0096     std::swap(_impl, other._impl);
0097     return *this;
0098 }
0099
0100 PaviaSfSet::PaviaSfSet() : SfSet() {
0101     _impl = new Impl();
0102 }
0103
0104 PaviaSfSet::~PaviaSfSet() {
0105     if (_impl != nullptr) {
0106         delete _impl;
0107     }
0108 }
0109
0110 double PaviaSfSet::F_UUT(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0111     double Q = std::sqrt(Q_sq);
0112     double qTtoQ = std::sqrt(ph_t_sq) / (z * Q);
0113     bool interp_x = false;
0114     bool interp_z = false;
0115     double old_x;
0116     double old_z;
0117     if (x < 0.02) {
0118         interp_x = true;
0119         old_x = x;
0120         x = 0.02;
0121     }
0122     if (z < 0.1) {
0123         interp_z = true;
0124         old_z = z;
0125         z = 0.1;
0126     }
0127     double result;
0128     if (h == Hadron::PI_P) {
0129         result = _impl->interp_pip_uu({ Q, x, z, qTtoQ });
0130     } else if (h == Hadron::PI_M) {
0131         result = _impl->interp_pim_uu({ Q, x, z, qTtoQ });
0132     } else {
0133         return std::numeric_limits<double>::quiet_NaN();
0134     }
0135     if (interp_x) {
0136         result *= old_x / 0.02;
0137     }
0138     if (interp_z) {
0139         result *= old_z / 0.1;
0140     }
0141     return result;
0142 }
0143
0144 double PaviaSfSet::F_UTT_sin_phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0145     double Q = std::sqrt(Q_sq);
0146     double qTtoQ = std::sqrt(ph_t_sq) / (z * Q);
0147     bool interp_x = false;
0148     bool interp_z = false;
0149     double old_x;
0150     double old_z;
0151     if (x < 0.02) {
0152         interp_x = true;
0153         old_x = x;
0154         x = 0.02;
0155     }
0156     if (z < 0.1) {
0157         interp_z = true;
0158         old_z = z;
0159         z = 0.1;
0160     }
0161     double result;
0162     if (h == Hadron::PI_P) {
0163         result = _impl->interp_pip_ut_sivers({ Q, x, z, qTtoQ });
0164     } else if (h == Hadron::PI_M) {
0165         result = _impl->interp_pim_ut_sivers({ Q, x, z, qTtoQ });
0166     } else {
0167         return std::numeric_limits<double>::quiet_NaN();
0168     }
0169     if (interp_x) {
0170         result *= old_x / 0.02;
0171     }
0172     if (interp_z) {
0173         result *= old_z / 0.1;
0174     }
0175     return result;
0176 }
0177