Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:24 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023 Duane Byer
0003 
0004 #include "Prokudin.h"
0005 
0006 #include "../interp/Interpolate.h"
0007 
0008 #include <mstwpdf.h>
0009 
0010 #include <array>
0011 #include <cmath>
0012 #include <fstream>
0013 #include <sstream>
0014 #include <string>
0015 
0016 #define SF_SET_DIR "grids"
0017 #define PROKUDIN_DIR "prokudin"
0018 
0019 // The following parameters come from appendix A in [2].
0020 
0021 namespace {
0022 
0023 unsigned const NUM_FLAVORS = 6;
0024 
0025 // We only have data files for proton structure function.
0026 double const M  = 0.93827208816;
0027 double const mh = 0.13957039;
0028 double const PI = 3.1415926535897932384626433;
0029 double const E  = 2.7182818284590452353602875;
0030 
0031 // Parameters for PDF.
0032 double const F1_MEAN_K_PERP_SQ = 0.25;
0033 double const G1_MEAN_K_PERP_SQ = 0.76 * F1_MEAN_K_PERP_SQ;
0034 double const D1_MEAN_P_PERP_SQ = 0.2;
0035 
0036 // Parameters for transversity.
0037 double const H1_ALPHA = 1.11;
0038 double const H1_BETA = 3.64;
0039 double const H1_N[6] = {
0040     0.46, -1.00, 0.,
0041     0., 0., 0.,
0042 };
0043 double const H1_MEAN_K_PERP_SQ = 0.25;
0044 // Reference [2] is not consistent about the choice of `mean_hT`, so here we
0045 // choose something in between the two options.
0046 double const HT_MEAN_K_PERP_SQ = 0.15;
0047 
0048 // Parameters for Collins.
0049 double const COLLINS_M_SQ = 1.50;
0050 double const COLLINS_M = std::sqrt(COLLINS_M_SQ);
0051 double const COLLINS_GAMMA = 1.06;
0052 double const COLLINS_DELTA = 0.07;
0053 double const COLLINS_N_FAV = 0.49;
0054 double const COLLINS_N_DISFAV = -1.;
0055 double const COLLINS_MEAN_P_PERP_SQ = (D1_MEAN_P_PERP_SQ * COLLINS_M_SQ)
0056     / (D1_MEAN_P_PERP_SQ + COLLINS_M_SQ);
0057 
0058 // Parameters for Sivers.
0059 double const SIVERS_M_1_SQ = 0.19;
0060 double const SIVERS_M_1 = std::sqrt(SIVERS_M_1_SQ);
0061 double const SIVERS_ALPHA[6] = {
0062     0.35, 0.44, 0.,
0063     0., 0., 0.,
0064 };
0065 double const SIVERS_BETA[6] = {
0066     2.6, 0.90, 0.,
0067     0., 0., 0.,
0068 };
0069 double const SIVERS_N[6] = {
0070     0.40, -0.97, 0.,
0071     0., 0., 0.,
0072 };
0073 double const SIVERS_MEAN_K_PERP_SQ = (F1_MEAN_K_PERP_SQ * SIVERS_M_1_SQ)
0074     / (F1_MEAN_K_PERP_SQ + SIVERS_M_1_SQ);
0075 
0076 // Parameters for Boer-Mulders.
0077 double const BM_M_1_SQ = 0.34;
0078 double const BM_M_1 = std::sqrt(BM_M_1_SQ);
0079 double const BM_ALPHA[6] = {
0080     0.73, 1.08, 0.79,
0081     0.79, 0.79, 0.79,
0082 };
0083 double const BM_BETA = 3.46;
0084 double const BM_A[6] = {
0085     0.35, -0.90, 0.24,
0086     0.04, -0.40, -1.,
0087 };
0088 double const BM_LAMBDA[6] = {
0089     2.1, -1.111, 0.,
0090     0., 0., 0.,
0091 };
0092 double const BM_MEAN_K_PERP_SQ = (F1_MEAN_K_PERP_SQ * BM_M_1_SQ)
0093     / (F1_MEAN_K_PERP_SQ + BM_M_1_SQ);
0094 
0095 // Parameters for pretzelosity.
0096 double const PRETZ_M_TT_SQ = 0.18;
0097 double const PRETZ_M_TT = std::sqrt(PRETZ_M_TT_SQ);
0098 double const PRETZ_ALPHA = 2.5;
0099 double const PRETZ_BETA = 2.;
0100 double const PRETZ_N[6] = {
0101     1., -1., 0.,
0102     0., 0., 0.,
0103 };
0104 
0105 double sq(double x) {
0106     return x * x;
0107 }
0108 double const PRETZ_MEAN_K_PERP_SQ = (F1_MEAN_K_PERP_SQ * PRETZ_M_TT_SQ)
0109     / (F1_MEAN_K_PERP_SQ + PRETZ_M_TT_SQ);
0110 
0111 double G(double ph_t_sq, double l) {
0112     // Equation [2.5.2].
0113     return std::exp(-ph_t_sq / l) / (PI * l);
0114 }
0115 double lambda(double z, double mean_kperp_sq, double mean_pperp_sq) {
0116     return sq(z) * mean_kperp_sq + mean_pperp_sq;
0117 }
0118 
0119 struct DataFileNotFoundError : public std::runtime_error {
0120     DataFileNotFoundError(char const* name) :
0121         std::runtime_error(std::string("Couldn't load data file '") + name + "'") {
0122     }
0123 };
0124 struct DataFileParseError : public std::runtime_error {
0125     DataFileParseError(char const* name) :
0126         std::runtime_error(std::string("Couldn't parse data file '") + name + "'") {
0127     }
0128 };
0129 
0130 // Finds a grid file.
0131 std::istream& find_file(std::ifstream& fin, char const* file_name) {
0132     fin.open(std::string("../share/" SF_SET_DIR "/" PROKUDIN_DIR "/") + file_name);
0133     if (fin) {
0134         return fin;
0135     }
0136     fin.open(std::string(SF_SET_DIR "/" PROKUDIN_DIR "/") + file_name);
0137     if (fin) {
0138         return fin;
0139     }
0140     fin.open(std::string(PROKUDIN_DIR "/") + file_name);
0141     if (fin) {
0142         return fin;
0143     }
0144     fin.open(file_name);
0145     if (fin) {
0146         return fin;
0147     }
0148     throw DataFileNotFoundError(file_name);
0149 }
0150 
0151 // Convenience method for loading grid data from a file. Automatically searches
0152 // in several possible directories for the file.
0153 template<typename T, std::size_t N, std::size_t K>
0154 std::array<Grid<T, N>, K> load_grids(char const* file_name) {
0155     std::ifstream in;
0156     find_file(in, file_name);
0157     std::vector<std::array<T, N + K> > data;
0158     std::string line;
0159     while (std::getline(in, line)) {
0160         std::stringstream line_in(line);
0161         std::array<T, N + K> next;
0162         for (std::size_t idx = 0; idx < N + K; ++idx) {
0163             line_in >> next[idx];
0164         }
0165         data.push_back(next);
0166         if (!line_in) {
0167             throw DataFileParseError(file_name);
0168         }
0169     }
0170     // By default, assume the grids are only accurate to single precision.
0171     return read_grids<T, N, K>(data);
0172 }
0173 
0174 double charge(unsigned fl) {
0175     switch (fl) {
0176     // Up.
0177     case 0:
0178         return 2./3.;
0179     // Down.
0180     case 1:
0181     // Strange.
0182     case 2:
0183         return -1./3.;
0184     // Up bar.
0185     case 3:
0186         return -2./3.;
0187     // Down bar.
0188     case 4:
0189     // Strange bar.
0190     case 5:
0191         return 1./3.;
0192     default:
0193         return 0.;
0194     }
0195 }
0196 
0197 struct ProkudinImpl {
0198     // PDF are calculated with the MSTWPDF library.
0199     std::ifstream file_pdf;
0200     mstw::c_mstwpdf pdf;
0201 
0202     // Load data files from WWSIDIS repository for the TMDs and FFs.
0203     std::array<Grid<double, 2>, 6> data_D1_pi_plus;
0204     std::array<Grid<double, 2>, 6> data_D1_pi_minus;
0205     std::array<Grid<double, 2>, 6> data_g1;
0206     std::array<Grid<double, 2>, 6> data_xgT;
0207     std::array<Grid<double, 2>, 2> data_xh1LperpM1;
0208     // Soffer bound.
0209     std::array<Grid<double, 2>, 6> data_sb;
0210 
0211     // Fragmentation functions.
0212     CubicView<double, 2> interp_D1_pi_plus[6];
0213     CubicView<double, 2> interp_D1_pi_minus[6];
0214 
0215     // Transverse momentum distributions.
0216     CubicView<double, 2> interp_g1[6];
0217     CubicView<double, 2> interp_xgT[6];
0218     CubicView<double, 2> interp_xh1LperpM1[2];
0219     CubicView<double, 2> interp_sb[6];
0220 
0221     ProkudinImpl() :
0222             file_pdf(),
0223             pdf(find_file(file_pdf, "mstw2008lo.00.dat"), false, true),
0224             data_D1_pi_plus(
0225                 load_grids<double, 2, 6>("fragmentationpiplus.dat")),
0226             data_D1_pi_minus(
0227                 load_grids<double, 2, 6>("fragmentationpiminus.dat")),
0228             data_g1(load_grids<double, 2, 6>("g1.dat")),
0229             data_xgT(load_grids<double, 2, 6>("gT_u_d_ubar_dbar_s_sbar.dat")),
0230             data_xh1LperpM1(load_grids<double, 2, 2>("xh1Lperp_u_d.dat")),
0231             data_sb(load_grids<double, 2, 6>("SofferBound.dat")),
0232             interp_D1_pi_plus{
0233                 CubicView<double, 2>(data_D1_pi_plus[0]),
0234                 CubicView<double, 2>(data_D1_pi_plus[1]),
0235                 CubicView<double, 2>(data_D1_pi_plus[2]),
0236                 CubicView<double, 2>(data_D1_pi_plus[3]),
0237                 CubicView<double, 2>(data_D1_pi_plus[4]),
0238                 CubicView<double, 2>(data_D1_pi_plus[5]),
0239             },
0240             interp_D1_pi_minus{
0241                 CubicView<double, 2>(data_D1_pi_minus[0]),
0242                 CubicView<double, 2>(data_D1_pi_minus[1]),
0243                 CubicView<double, 2>(data_D1_pi_minus[2]),
0244                 CubicView<double, 2>(data_D1_pi_minus[3]),
0245                 CubicView<double, 2>(data_D1_pi_minus[4]),
0246                 CubicView<double, 2>(data_D1_pi_minus[5]),
0247             },
0248             interp_g1{
0249                 CubicView<double, 2>(data_g1[0]),
0250                 CubicView<double, 2>(data_g1[1]),
0251                 CubicView<double, 2>(data_g1[2]),
0252                 CubicView<double, 2>(data_g1[3]),
0253                 CubicView<double, 2>(data_g1[4]),
0254                 CubicView<double, 2>(data_g1[5]),
0255             },
0256             interp_xgT{
0257                 CubicView<double, 2>(data_xgT[0]),
0258                 CubicView<double, 2>(data_xgT[1]),
0259                 CubicView<double, 2>(data_xgT[4]),
0260                 CubicView<double, 2>(data_xgT[2]),
0261                 CubicView<double, 2>(data_xgT[3]),
0262                 CubicView<double, 2>(data_xgT[5]),
0263             },
0264             interp_xh1LperpM1{
0265                 CubicView<double, 2>(data_xh1LperpM1[0]),
0266                 CubicView<double, 2>(data_xh1LperpM1[1]),
0267             },
0268             interp_sb{
0269                 CubicView<double, 2>(data_sb[0]),
0270                 CubicView<double, 2>(data_sb[1]),
0271                 CubicView<double, 2>(data_sb[2]),
0272                 CubicView<double, 2>(data_sb[3]),
0273                 CubicView<double, 2>(data_sb[4]),
0274                 CubicView<double, 2>(data_sb[5]),
0275             } {
0276         file_pdf.close();
0277     }
0278 };
0279 
0280 }
0281 
0282 struct ProkudinSfSet::Impl {
0283     ProkudinImpl impl;
0284 };
0285 
0286 ProkudinSfSet::ProkudinSfSet(ProkudinSfSet&& other) noexcept :
0287         _impl(nullptr) {
0288     std::swap(_impl, other._impl);
0289 }
0290 ProkudinSfSet& ProkudinSfSet::operator=(ProkudinSfSet&& other) noexcept {
0291     std::swap(_impl, other._impl);
0292     return *this;
0293 }
0294 
0295 ProkudinSfSet::ProkudinSfSet() : SfSet() {
0296     _impl = new Impl();
0297 }
0298 
0299 ProkudinSfSet::~ProkudinSfSet() {
0300     if (_impl != nullptr) {
0301         delete _impl;
0302     }
0303 }
0304 
0305 double ProkudinSfSet::F_UUT(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0306     // Equation [2.5.1a].
0307     double result = 0.;
0308     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0309         result += sq(charge(fl))*xf1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0310     }
0311     double l = lambda(z, F1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0312     return G(ph_t_sq, l)*result;
0313 }
0314 double ProkudinSfSet::F_UU_cos_phih(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0315     // Equation [2.7.9a].
0316     double Q = std::sqrt(Q_sq);
0317     double ph_t = std::sqrt(ph_t_sq);
0318     double result = 0.;
0319     // Uses a WW-type approximation to rewrite in terms of `xf1`.
0320     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0321         result += sq(charge(fl))*xf1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0322     }
0323     double l = lambda(z, F1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0324     return -2.*F1_MEAN_K_PERP_SQ/Q*ph_t*(z/l)*G(ph_t_sq, l)*result;
0325 }
0326 double ProkudinSfSet::F_UU_cos_2phih(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0327     // Equation [2.5.9a].
0328     double result = 0.;
0329     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0330         result += sq(charge(fl))*xh1perpM1(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0331     }
0332     double l = lambda(z, BM_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0333     return 4.*M*mh*ph_t_sq*sq(z/l)*G(ph_t_sq, l)*result;
0334 }
0335 
0336 double ProkudinSfSet::F_UL_sin_phih(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0337     // Equation [2.7.6a].
0338     double Q = std::sqrt(Q_sq);
0339     double ph_t = std::sqrt(ph_t_sq);
0340     double result = 0.;
0341     // Use WW-type approximation to rewrite in terms of `xh1LperpM1`.
0342     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0343         result += sq(charge(fl))*xh1LperpM1(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0344     }
0345     // Approximate width with `H1_MEAN_K_PERP_SQ`.
0346     double l = lambda(z, H1_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0347     return -8.*M*mh*z*ph_t/(Q*l)*G(ph_t_sq, l)*result;
0348 }
0349 double ProkudinSfSet::F_UL_sin_2phih(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0350     // Equation [2.6.2a].
0351     double result = 0.;
0352     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0353         result += sq(charge(fl))*xh1LperpM1(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0354     }
0355     // Approximate width with `H1_MEAN_K_PERP_SQ`.
0356     double l = lambda(z, H1_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0357     return 4.*M*mh*ph_t_sq*sq(z/l)*G(ph_t_sq, l)*result;
0358 }
0359 
0360 double ProkudinSfSet::F_UTT_sin_phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0361     // Equation [2.5.7a].
0362     double ph_t = std::sqrt(ph_t_sq);
0363     double result = 0.;
0364     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0365         result += sq(charge(fl))*xf1TperpM1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0366     }
0367     double l = lambda(z, SIVERS_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0368     return -2.*M*z*ph_t/l*G(ph_t_sq, l)*result;
0369 }
0370 double ProkudinSfSet::F_UT_sin_2phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0371     // Equation [2.7.8a].
0372     double Q = std::sqrt(Q_sq);
0373     double result_1 = 0.;
0374     // Use WW-type approximation to rewrite in terms of `xf1TperpM1`.
0375     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0376         result_1 += sq(charge(fl))*xf1TperpM1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0377     }
0378     double result_2 = 0.;
0379     // Use WW-type approximation to rewrite in terms of `h1TperpM2`.
0380     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0381         result_2 += sq(charge(fl))*xh1TperpM2(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0382     }
0383     // Approximate width with `SIVERS_MEAN_K_PERP_SQ`.
0384     double l_1 = lambda(z, SIVERS_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0385     // Approximate width with `PRETZ_MEAN_K_PERP_SQ`.
0386     double l_2 = lambda(z, PRETZ_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0387     // The paragraph following [2.7.8a] has a mistake in the WW-type
0388     // approximation linking `h1TM1 + h1TperpM1` with `h1TperpM2`, due to a
0389     // missing factor of 2.
0390     return 2.*M*ph_t_sq/Q*(
0391         SIVERS_MEAN_K_PERP_SQ*sq(z/l_1)*G(ph_t_sq, l_1)*result_1
0392         - 2.*M*mh*sq(z/l_2)*G(ph_t_sq, l_2)*result_2);
0393 }
0394 double ProkudinSfSet::F_UT_sin_3phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0395     // Equation [2.5.10a].
0396     double ph_t = std::sqrt(ph_t_sq);
0397     double result = 0.;
0398     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0399         result += sq(charge(fl))*xh1TperpM2(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0400     }
0401     double l = lambda(z, PRETZ_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0402     return 2.*sq(M)*mh*std::pow(z*ph_t/l, 3)*G(ph_t_sq, l)*result;
0403 }
0404 double ProkudinSfSet::F_UT_sin_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0405     // Equation [2.7.7a].
0406     double Q = std::sqrt(Q_sq);
0407     double result = 0.;
0408     // WW-type approximation used here (see [2] for details).
0409     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0410         result += sq(charge(fl))*xh1M1(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0411     }
0412     double l = lambda(z, PRETZ_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0413     return 8.*sq(M)*mh*sq(z)/(Q*l)*(1. - ph_t_sq/l)*G(ph_t_sq, l)*result;
0414 }
0415 double ProkudinSfSet::F_UT_sin_phih_p_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0416     // Equation [2.5.8a].
0417     double ph_t = std::sqrt(ph_t_sq);
0418     double result = 0.;
0419     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0420         result += sq(charge(fl))*xh1(fl, x, Q_sq)*H1perpM1(h, fl, z, Q_sq);
0421     }
0422     double l = lambda(z, H1_MEAN_K_PERP_SQ, COLLINS_MEAN_P_PERP_SQ);
0423     return 2.*mh*z*ph_t/l*G(ph_t_sq, l)*result;
0424 }
0425 
0426 double ProkudinSfSet::F_LL(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0427     // Equation [2.5.5a].
0428     double result = 0.;
0429     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0430         result += sq(charge(fl))*xg1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0431     }
0432     double l = lambda(z, G1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0433     return G(ph_t_sq, l)*result;
0434 }
0435 double ProkudinSfSet::F_LL_cos_phih(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0436     // Equation [2.7.5a].
0437     double Q = std::sqrt(Q_sq);
0438     double ph_t = std::sqrt(ph_t_sq);
0439     double result = 0.;
0440     // Uses a WW-type approximation to rewrite in terms of `xg1`.
0441     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0442         result += sq(charge(fl))*xg1(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0443     }
0444     // Approximate width with `G1_MEAN_K_PERP_SQ`.
0445     double l = lambda(z, G1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0446     return -2.*G1_MEAN_K_PERP_SQ*z*ph_t/(Q*l)*G(ph_t_sq, l)*result;
0447 }
0448 
0449 double ProkudinSfSet::F_LT_cos_phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0450     // Equation [2.6.1a].
0451     double ph_t = std::sqrt(ph_t_sq);
0452     double result = 0.;
0453     // Uses a WW-type approximation to rewrite in terms of `xgT`.
0454     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0455         result += sq(charge(fl))*xgT(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0456     }
0457     // Approximate width with `G1_MEAN_K_PERP_SQ`.
0458     double l = lambda(z, G1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0459     return 2.*M*x*z*ph_t/l*G(ph_t_sq, l)*result;
0460 }
0461 double ProkudinSfSet::F_LT_cos_2phih_m_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0462     // Equation [2.7.4a].
0463     double Q = std::sqrt(Q_sq);
0464     double result = 0.;
0465     // Uses a WW-type approximation to rewrite in terms of `xgT`.
0466     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0467         result += sq(charge(fl))*xgT(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0468     }
0469     // Approximate width with `G1_MEAN_K_PERP_SQ`.
0470     double l = lambda(z, G1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0471     return -2.*G1_MEAN_K_PERP_SQ*M*x*ph_t_sq*sq(z/l)/Q*G(ph_t_sq, l)*result;
0472 }
0473 double ProkudinSfSet::F_LT_cos_phis(Hadron h, double x, double z, double Q_sq, double ph_t_sq) const {
0474     // Equation [2.7.2a].
0475     double Q = std::sqrt(Q_sq);
0476     double result = 0.;
0477     for (unsigned fl = 0; fl < NUM_FLAVORS; ++fl) {
0478         result += sq(charge(fl))*xgT(fl, x, Q_sq)*D1(h, fl, z, Q_sq);
0479     }
0480     // Approximate width with `G1_MEAN_K_PERP_SQ`.
0481     double l = lambda(z, G1_MEAN_K_PERP_SQ, D1_MEAN_P_PERP_SQ);
0482     return -2.*M*x/Q*G(ph_t_sq, l)*result;
0483 }
0484 
0485 // Fragmentation functions.
0486 double ProkudinSfSet::D1(Hadron h, unsigned fl, double z, double Q_sq) const {
0487     switch (h) {
0488     case Hadron::PI_P:
0489         return _impl->impl.interp_D1_pi_plus[fl]({ z, Q_sq });
0490     case Hadron::PI_M:
0491         return _impl->impl.interp_D1_pi_minus[fl]({ z, Q_sq });
0492     default:
0493         return 0.;
0494     }
0495 }
0496 double ProkudinSfSet::H1perpM1(Hadron h, unsigned fl, double z, double Q_sq) const {
0497     double collins_coeff = 0.;
0498     if (h == Hadron::PI_P) {
0499         if (fl == 0 || fl == 4) {
0500             collins_coeff = COLLINS_N_FAV;
0501         } else if (fl == 1 || fl == 3) {
0502             collins_coeff = COLLINS_N_DISFAV;
0503         }
0504     } else if (h == Hadron::PI_M) {
0505         if (fl == 1 || fl == 3) {
0506             collins_coeff = COLLINS_N_FAV;
0507         } else if (fl == 0 || fl == 4) {
0508             collins_coeff = COLLINS_N_DISFAV;
0509         }
0510     } else {
0511         return 0.;
0512     }
0513     return std::sqrt(E/2.)/(z*mh*COLLINS_M)
0514         *sq(COLLINS_MEAN_P_PERP_SQ)/D1_MEAN_P_PERP_SQ
0515         *collins_coeff
0516         *std::pow(z, COLLINS_GAMMA)*std::pow(1. - z, COLLINS_DELTA)
0517         *std::pow(COLLINS_GAMMA + COLLINS_DELTA, COLLINS_GAMMA + COLLINS_DELTA)
0518         *std::pow(COLLINS_GAMMA, -COLLINS_GAMMA)
0519         *std::pow(COLLINS_DELTA, -COLLINS_DELTA)
0520         *D1(h, fl, z, Q_sq);
0521 }
0522 
0523 // Parton distribution functions.
0524 double ProkudinSfSet::xf1(unsigned fl, double x, double Q_sq) const {
0525     double Q = std::sqrt(Q_sq);
0526     switch (fl) {
0527     case 0:
0528         return _impl->impl.pdf.parton(8, x, Q) + _impl->impl.pdf.parton(-2, x, Q);
0529     case 1:
0530         return _impl->impl.pdf.parton(7, x, Q) + _impl->impl.pdf.parton(-1, x, Q);
0531     case 2:
0532         return _impl->impl.pdf.parton(3, x, Q);
0533     case 3:
0534         return _impl->impl.pdf.parton(-2, x, Q);
0535     case 4:
0536         return _impl->impl.pdf.parton(-1, x, Q);
0537     case 5:
0538         return _impl->impl.pdf.parton(-3, x, Q);
0539     default:
0540         return 0.;
0541     }
0542 }
0543 
0544 // Transverse momentum distributions.
0545 double ProkudinSfSet::xf1TperpM1(unsigned fl, double x, double Q_sq) const {
0546     // Equation [2.A.4].
0547     return -std::sqrt(E/2.)/(M*SIVERS_M_1)
0548         *sq(SIVERS_MEAN_K_PERP_SQ)/F1_MEAN_K_PERP_SQ
0549         *SIVERS_N[fl]
0550         *std::pow(x, SIVERS_ALPHA[fl])*std::pow(1. - x, SIVERS_BETA[fl])
0551         *std::pow(SIVERS_ALPHA[fl] + SIVERS_BETA[fl],
0552             SIVERS_ALPHA[fl] + SIVERS_BETA[fl])
0553         *std::pow(SIVERS_ALPHA[fl], -SIVERS_ALPHA[fl])
0554         *std::pow(SIVERS_BETA[fl], -SIVERS_BETA[fl])
0555         *xf1(fl, x, Q_sq);
0556 }
0557 double ProkudinSfSet::xg1(unsigned fl, double x, double Q_sq) const {
0558     return x*_impl->impl.interp_g1[fl]({ x, Q_sq });
0559 }
0560 double ProkudinSfSet::xgT(unsigned fl, double x, double Q_sq) const {
0561     return _impl->impl.interp_xgT[fl]({ x, Q_sq });
0562 }
0563 double ProkudinSfSet::xh1(unsigned fl, double x, double Q_sq) const {
0564     // Use the Soffer bound to get an upper limit on transversity (Equation
0565     // [2.A.7]).
0566     return x*H1_N[fl]
0567         *std::pow(x, H1_ALPHA)*std::pow(1. - x, H1_BETA)
0568         *std::pow(H1_ALPHA + H1_BETA, H1_ALPHA + H1_BETA)
0569         *std::pow(H1_ALPHA, -H1_ALPHA)
0570         *std::pow(H1_BETA, -H1_BETA)
0571         *_impl->impl.interp_sb[fl]({ x, Q_sq });
0572 }
0573 double ProkudinSfSet::xh1M1(unsigned fl, double x, double Q_sq) const {
0574     return H1_MEAN_K_PERP_SQ/(2.*sq(M))*xh1(fl, x, Q_sq);
0575 }
0576 double ProkudinSfSet::xh1LperpM1(unsigned fl, double x, double Q_sq) const {
0577     // Data only exists for up and down quarks.
0578     if (!(fl == 0 || fl == 1)) {
0579         return 0.;
0580     } else {
0581         return _impl->impl.interp_xh1LperpM1[fl]({ x, Q_sq });
0582     }
0583 }
0584 double ProkudinSfSet::xh1TperpM2(unsigned fl, double x, double Q_sq) const {
0585     // Equation [2.A.24].
0586     return E/(2.*sq(M)*PRETZ_M_TT_SQ)
0587         *std::pow(PRETZ_MEAN_K_PERP_SQ, 3)/F1_MEAN_K_PERP_SQ
0588         *PRETZ_N[fl]
0589         *std::pow(x, PRETZ_ALPHA)*std::pow(1. - x, PRETZ_BETA)
0590         *std::pow(PRETZ_ALPHA + PRETZ_BETA, PRETZ_ALPHA + PRETZ_BETA)
0591         *std::pow(PRETZ_ALPHA, -PRETZ_ALPHA)
0592         *std::pow(PRETZ_BETA, -PRETZ_BETA)
0593         *(xf1(fl, x, Q_sq) - xg1(fl, x, Q_sq));
0594 }
0595 double ProkudinSfSet::xh1perpM1(unsigned fl, double x, double Q_sq) const {
0596     // Equation [2.A.18].
0597     return -std::sqrt(E/2.)/(M*BM_M_1)
0598         *sq(BM_MEAN_K_PERP_SQ)/F1_MEAN_K_PERP_SQ
0599         *BM_LAMBDA[fl] * BM_A[fl]
0600         *std::pow(x, BM_ALPHA[fl])*std::pow(1. - x, BM_BETA)
0601         *std::pow(BM_ALPHA[fl] + BM_BETA, BM_ALPHA[fl] + BM_BETA)
0602         *std::pow(BM_ALPHA[fl], -BM_ALPHA[fl])
0603         *std::pow(BM_BETA, -BM_BETA)
0604         *xf1(fl, x, Q_sq);
0605 }
0606
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!