/* This file calculates the parton distribution functions. */ /**************************************************************************/ /* */ /* IMPORTANT : all dimensional quantities in this program */ /* are in units GeV */ /* */ /**************************************************************************/ #include #include #include static int param_par_s; static double scale_sqr_s, n_s; /*************************************************************************/ /* */ /* UNPOLARIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* the function up_valence() gives the unpolarized valence up quark distribution as function of x for different parametrizations */ double up_valence(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double n1, n2, eud, gud, norm, test; n1 = 0.26; n2 = 3.82; eud = 14.4; gud = 16.99; norm = 9.790758; /* test = 1./0.26 * ( 3.82 * pow(x,n1) * pow(1. - x, n2 - 1.) * (1. + eud * pow(x,0.5) + gud * x) - pow(1. - x, n2) * ( eud/2. * pow(x,n1 - 0.5) + gud * pow(x,n1) ) ); */ return 3./norm * 1./pow(x, 1. - n1) * pow(1. - x, n2) * (1. + eud * pow(x,0.5) + gud * x) - down_valence(param, x, scale_sqr); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { double n1, n2, eud, gud, norm, test; n1 = 0.45; n2 = 3.91; eud = 2.46; gud = 3.32; norm = 1.823211; return 3./norm * 1./pow(x, 1. - n1) * pow(1. - x, n2) * (1. + eud * pow(x,0.5) + gud * x) - down_valence(param, x, scale_sqr); break; } case 3: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double n1, n2, eud, gud, norm; n1 = 0.26; n2 = 3.82; eud = 14.4; gud = 16.99; norm = 9.790758; return 3./norm * 1./pow(x, 1. - n1) * pow(1. - x, n2) * (1. + eud * pow(x,0.5) + gud * x) - down_valence(param, x, scale_sqr); break; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return upv/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return upv/x; break; } case 6: /* MRST98 (central gluon) at Q^2 = 1 GeV^2 */ { double A_u, eta_1, eta_2, eps_u, gamma_u; A_u = 0.6051; eta_1 = 0.4089; eta_2 = 3.395; eps_u = 2.078; gamma_u = 14.56; return A_u * pow(x, -1. + eta_1) * pow(1. - x, eta_2) * (1. + eps_u * sqrt(x) + gamma_u * x); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { double A_u, eta_1, eta_2, eps_u, gamma_u; A_u = 0.7763; eta_1 = 0.4398; eta_2 = 3.427; eps_u = 1.152; gamma_u = 12.36; return A_u * pow(x, -1. + eta_1) * pow(1. - x, eta_2) * (1. + eps_u * sqrt(x) + gamma_u * x); break; } case 8: /* CTEQ4L (LO) */ { int iparton1, iparton2; /* iset = 3; setctq4_(&iset);*/ iparton1=1; iparton2=-1; return CTQ4PDF(iparton1, x, Q) - CTQ4PDF(iparton2, x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq4_(&iset);*/ iparton1=1; iparton2=-1; return CTQ4PDF(iparton1, x, Q) - CTQ4PDF(iparton2, x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq5_(&iset);*/ iparton1=1; iparton2=-1; return CTQ5PDF(iparton1, x, Q) - CTQ5PDF(iparton2, x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return upv/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return upv/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return upv/x; break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return upv/x; break; } case 15 : /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { double A_u, eta_1, eta_2, eps_u, gamma_u; A_u = 0.262; eta_1 = 0.31; eta_2 = 3.50; eps_u = 3.83; gamma_u = 37.65; return A_u * pow(x, -1. + eta_1) * pow(1. - x, eta_2) * (1. + eps_u * sqrt(x) + gamma_u * x); break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return upv/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq6_(&iset);*/ iparton1=1; iparton2=-1; return CTQ6PDF(iparton1, x, Q) - CTQ6PDF(iparton2, x, Q); break; } } } double norm_up_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_up_valence, 0., 1., parton_npoints); } double int_up_valence(double x) { return up_valence(param_par_s, x, scale_sqr_s); } double mom_up_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_up_valence, 0., 1., parton_npoints); } double int_mom_up_valence(double x) { return x * up_valence(param_par_s, x, scale_sqr_s); } double moment_up_valence(int param, int parton_npoints, double scale_sqr, double n) { param_par_s = param; scale_sqr_s = scale_sqr; n_s = n; return int_gauss7(int_moment_up_valence, 0., 1., parton_npoints); } double int_moment_up_valence(double x) { return pow(x,n_s) * up_valence(param_par_s, x, scale_sqr_s); } /* The function down_valence() gives the unpolarized valence down quark distribution as function of x for different parametrizations */ double down_valence(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double n3, n4, ed, gd, norm, test; n3 = 0.78; n4 = 4.57; ed = -0.87; gd = 0.82; norm = 0.263275; /* test = 1./0.78 * ( 4.57 * pow(x,n3) * pow(1. - x, n4 - 1.) * (1. + ed * pow(x,0.5) + gd * x) + pow(1. - x, n4) * (- ed/2. * pow(x,n3 - 0.5) - gd * pow(x,n3) ) ); */ return 1./norm * 1./pow(x, 1. - n3) * pow(1. - x, n4) * (1. + ed * pow(x,0.5) + gd * x); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { double n3, n4, ed, gd, norm, test; n3 = 0.35; n4 = 4.66; ed = 11.4; gd = 3.0; norm = 4.599439; return 1./norm * 1./pow(x, 1. - n3) * pow(1. - x, n4) * (1. + ed * pow(x,0.5) + gd * x); break; } case 3: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double n3, n4, ed, gd, norm; n3 = 0.78; n4 = 4.57; ed = -0.87; gd = 0.82; norm = 0.263275; return 1./norm * 1./pow(x, 1. - n3) * pow(1. - x, n4) * (1. + ed * pow(x,0.5) + gd * x); break; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return dnv/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dnv/x; break; } case 6: /* MRST98 (central gluon) at Q^2 = 1 GeV^2 */ { double A_d, eta_3, eta_4, eps_d, gamma_d; A_d = 0.05811; eta_3 = 0.2882; eta_4 = 3.874; eps_d = 34.69; gamma_d = 28.96; return A_d * pow(x, -1. + eta_3) * pow(1. - x, eta_4) * (1. + eps_d * sqrt(x) + gamma_d * x); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { double A_d, eta_3, eta_4, eps_d, gamma_d; A_d = 0.06015; eta_3 = 0.2694; eta_4 = 3.941; eps_d = 27.96; gamma_d = 38.35; return A_d * pow(x, -1. + eta_3) * pow(1. - x, eta_4) * (1. + eps_d * sqrt(x) + gamma_d * x); break; } case 8: /* CTEQ4L (LO) */ { int iparton1, iparton2; /* iset = 3; setctq4_(&iset); */ iparton1=2; iparton2=-2; return CTQ4PDF(iparton1, x, Q) - CTQ4PDF(iparton2, x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq4_(&iset); */ iparton1=2; iparton2=-2; return CTQ4PDF(iparton1, x, Q) - CTQ4PDF(iparton2, x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq5_(&iset); */ iparton1=2; iparton2=-2; return CTQ5PDF(iparton1, x, Q) - CTQ5PDF(iparton2, x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dnv/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dnv/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dnv/x; break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dnv/x; break; } case 15 : /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { double A_d, eta_3, eta_4, eps_d, gamma_d; A_d = 0.061; eta_3 = 0.35; eta_4 = 4.03; eps_d = 49.05; gamma_d = 8.65; return A_d * pow(x, -1. + eta_3) * pow(1. - x, eta_4) * (1. + eps_d * sqrt(x) + gamma_d * x); break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dnv/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton1, iparton2; /* iset = 1; setctq6_(&iset); */ iparton1=2; iparton2=-2; return CTQ6PDF(iparton1, x, Q) - CTQ6PDF(iparton2, x, Q); break; } } } double norm_down_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_down_valence, 0., 1., parton_npoints); } double int_down_valence(double x) { return down_valence(param_par_s, x, scale_sqr_s); } double mom_down_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_down_valence, 0., 1., parton_npoints); } double int_mom_down_valence(double x) { return x * down_valence(param_par_s, x, scale_sqr_s); } double moment_down_valence(int param, int parton_npoints, double scale_sqr, double n) { param_par_s = param; scale_sqr_s = scale_sqr; n_s = n; return int_gauss7(int_moment_down_valence, 0., 1., parton_npoints); } double int_moment_down_valence(double x) { return pow(x,n_s) * down_valence(param_par_s, x, scale_sqr_s); } /* The function up_bar() gives the u_bar quark distribution (remark total up quark sea is obtained as u_sea + u_bar = 2 * u_bar) as function of x for different parametrizations */ double up_bar(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 1./5. * sea_distr(1, x, scale_sqr); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { return .5 * ( 2./5. * sea_distr(2, x, scale_sqr) - diff_distr(2, x, scale_sqr) ); break; } case 3: { return 0.; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return usea/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return usea/x; break; } case 6: /* MRST98 at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(6, x, scale_sqr) - diff_distr(6, x, scale_sqr) ); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(6, x, scale_sqr) - diff_distr(6, x, scale_sqr) ); break; } case 8: /* CTEQ4L (L0) */ { int iparton; /* iset = 3; setctq4_(&iset); */ iparton =-1; return CTQ4PDF(iparton , x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton; /* iset = 1; setctq4_(&iset); */ iparton =-1; return CTQ4PDF(iparton , x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton; /* iset = 1; setctq5_(&iset); */ iparton =-1; return CTQ5PDF(iparton , x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return usea/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return usea/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return usea/x; break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return usea/x; break; } case 15 : /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(15, x, scale_sqr) - diff_distr(15, x, scale_sqr) ); break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return usea/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton; /* iset = 1; setctq6_(&iset); */ iparton =-1; return CTQ6PDF(iparton , x, Q); break; } } } double mom_up_bar(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_up_bar, 0., 1., parton_npoints); } double int_mom_up_bar(double x) { return 2. * x * up_bar(param_par_s, x, scale_sqr_s); } /* The function down_bar() gives the d_bar quark distribution (remark total down quark sea is obtained as d_sea + d_bar = 2 * d_bar) as function of x for different parametrizations */ double down_bar(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 1./5. * sea_distr(1, x, scale_sqr); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { return .5 * ( 2./5. * sea_distr(2, x, scale_sqr) + diff_distr(2, x, scale_sqr) ); break; } case 3: { return 0.; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return dsea/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dsea/x; break; } case 6: /* MRST98 at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(6, x, scale_sqr) + diff_distr(6, x, scale_sqr) ); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(6, x, scale_sqr) + diff_distr(6, x, scale_sqr) ); break; } case 8: /* CTEQ4L (LO) */ { int iparton; /* iset = 3; setctq4_(&iset); */ iparton =-2; return CTQ4PDF(iparton , x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton; /* iset = 1; setctq4_(&iset); */ iparton =-2; return CTQ4PDF(iparton , x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton; /* iset = 1; setctq5_(&iset); */ iparton =-2; return CTQ5PDF(iparton, x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dsea/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dsea/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return dsea/x; break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dsea/x; break; } case 15: /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { return .5 * ( 2./5. * sea_distr(15, x, scale_sqr) + diff_distr(15, x, scale_sqr) ); break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return dsea/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton; /* iset = 1; setctq6_(&iset); */ iparton =-2; return CTQ6PDF(iparton, x, Q); break; } } } double mom_down_bar(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_down_bar, 0., 1., parton_npoints); } double int_mom_down_bar(double x) { return 2. * x * down_bar(param_par_s, x, scale_sqr_s); } /* The function strange_sea() gives the strange sea quark distribution as function of x for different parametrizations */ double strange_sea(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 1./2. * 1./5. * sea_distr(1, x, scale_sqr); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { return 1./2. * 1./5. * sea_distr(2, x, scale_sqr); break; } case 3: { return 0.; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return str/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return str/x; break; } case 6: /* MRST98 at Q^2 = 1 GeV^2 */ { return 1./2. * 1./5. * sea_distr(6, x, scale_sqr); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { return 1./2. * 1./5. * sea_distr(6, x, scale_sqr); break; } case 8: /* CTEQ4L (LO) */ { int iparton; /* iset = 3; setctq4_(&iset); */ iparton =3; return CTQ4PDF(iparton , x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton; /* iset = 1; setctq4_(&iset); */ iparton =3; return CTQ4PDF(iparton , x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton; /* iset = 1; setctq5_(&iset); */ iparton =3; return CTQ5PDF(iparton , x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return str/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return str/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return str/x; break; } case 14: /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return str/x; break; } case 15: /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { return 1./2. * 1./5. * sea_distr(15, x, scale_sqr); break; } case 16: /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return str/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton; /* iset = 1; setctq6_(&iset); */ iparton =3; return CTQ6PDF(iparton , x, Q); break; } } } double mom_strange_sea(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_strange_sea, 0., 1., parton_npoints); } double int_mom_strange_sea(double x) { return 2. * x * strange_sea(param_par_s, x, scale_sqr_s); } /* The function sea_distr() gives the TOTAL sea quark distribution as function of x for different parametrizations */ double sea_distr(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double ns, es, gs, As; ns = 10.; es = -2.21; gs = 6.22; As = 1.87; return As * 1./x * pow(1. - x, ns) * (1. + es * pow(x,0.5) + gs * x); break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { double ns, es, gs, As; ns = 10.; es = -2.68; gs = 7.38; As = 1.93; return As * 1./x * pow(1. - x, ns) * (1. + es * pow(x,0.5) + gs * x); break; } case 3: { return 0.; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; /*Return q_sea(x) (u_sea(x)+ubar_sea(x)+d_sea(x)+dbar_sea(x)u+ s_sea(x)+sbar_sea(x)=2.*(u_bar(x)+d_bar(x)+s_bar(x))) :*/ MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return (2. * usea + 2. * dsea + 2. * str)/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (2. * usea + 2. * dsea + 2. * str)/x; break; } case 6: /* MRST98 at Q^2 = 1 GeV^2 */ { double lams, ns, es, gs, As; lams = 0.2712; ns = 7.808; es = 2.283; gs = 20.69; As = 0.2004; return As * pow(x, -1. - lams) * pow(1. - x, ns) * (1. + es * pow(x,0.5) + gs * x); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { double lams, ns, es, gs, As; lams = 0.2819; ns = 8.212; es = 3.725; gs = 21.80; As = 0.1786; return As * pow(x, -1. - lams) * pow(1. - x, ns) * (1. + es * pow(x,0.5) + gs * x); break; } case 8: /* CTEQ4L (LO) */ { int iparton1, iparton2, iparton3; /* iset = 3; setctq4_(&iset); */ iparton1=-1; iparton2=-2; iparton3=3; /*Return q_sea(x) (u_sea(x)+ubar_sea(x)+d_sea(x)+dbar_sea(x)u+ s_sea(x)+sbar_sea(x)=2.*(u_bar(x)+d_bar(x)+s_bar(x))) :*/ return 2.*( CTQ4PDF(iparton1,x,Q) + CTQ4PDF(iparton2,x,Q)+ + CTQ4PDF(iparton3,x,Q)); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton1, iparton2, iparton3; /* iset = 1; setctq4_(&iset); */ iparton1=-1; iparton2=-2; iparton3=3; /*Return q_sea(x) (u_sea(x)+ubar_sea(x)+d_sea(x)+dbar_sea(x)u+ s_sea(x)+sbar_sea(x)=2.*(u_bar(x)+d_bar(x)+s_bar(x))) :*/ return 2.*( CTQ4PDF(iparton1,x,Q) + CTQ4PDF(iparton2,x,Q)+ + CTQ4PDF(iparton3,x,Q)); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton1, iparton2, iparton3; /* iset = 1; setctq5_(&iset); */ iparton1=-1; iparton2=-2; iparton3=3; /*Return q_sea(x) (u_sea(x)+ubar_sea(x)+d_sea(x)+dbar_sea(x)u+ s_sea(x)+sbar_sea(x)=2.*(u_bar(x)+d_bar(x)+s_bar(x))) :*/ return 2.*( CTQ5PDF(iparton1,x,Q) + CTQ5PDF(iparton2,x,Q)+ + CTQ5PDF(iparton3,x,Q)); break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (2. * usea + 2. * dsea + 2. * str)/x; break; } case 15: /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { double lams, ns, es, gs, As; lams = 0.12; ns = 7.66; es = -1.34; gs = 7.40; As = 0.759; return As * pow(x, -1. - lams) * pow(1. - x, ns) * (1. + es * pow(x,0.5) + gs * x); break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (2. * usea + 2. * dsea + 2. * str)/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton1, iparton2, iparton3; /* iset = 1; setctq6_(&iset); */ iparton1=-1; iparton2=-2; iparton3=3; /*Return q_sea(x) (u_sea(x)+ubar_sea(x)+d_sea(x)+dbar_sea(x)u+ s_sea(x)+sbar_sea(x)=2.*(u_bar(x)+d_bar(x)+s_bar(x))) :*/ return 2.*( CTQ6PDF(iparton1,x,Q) + CTQ6PDF(iparton2,x,Q)+ + CTQ6PDF(iparton3,x,Q)); break; } } } /* The function diff_distr() gives the difference of the down_bar - up_bar quark distribution as function of x for different parametrizations */ double diff_distr(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 0.; break; } case 2: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, D_o */ { double ns, nd, Ad; ns = 10.; nd = 0.45; Ad = 0.163; return Ad * 1./x * pow(x,nd) * pow(1. - x, ns); break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (dsea - usea)/x; break; } case 6: /* MRST98 at Q^2 = 1 GeV^2 */ { double ns, nd, Ad, gd, dd; ns = 7.808; nd = 1.183; Ad = 1.290; gd = 9.987; dd = -33.34; return Ad * 1./x * pow(x,nd) * pow(1. - x, ns + 2.) * (1. + gd * x + dd * pow(x, 2.)); break; } case 7: /* MRST98 (higher gluon, k_perp_av = 0 GeV) at Q^2 = 1 GeV^2 */ { double ns, nd, Ad, gd, dd; ns = 8.212; nd = 1.157; Ad = 1.260; gd = 9.778; dd = -30.83; return Ad * 1./x * pow(x,nd) * pow(1. - x, ns + 2.) * (1. + gd * x + dd * pow(x, 2.)); break; } case 14: /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (dsea - usea)/x; break; } case 15: /* MRST02 NNLO at Q^2 = 1 GeV^2 */ { double ns, nd, Ad, gd, dd; ns = 9.66; nd = 1.24; Ad = 1.432; gd = 9.86; dd = -29.04; return Ad * 1./x * pow(x,nd) * pow(1. - x, ns) * (1. + gd * x + dd * pow(x, 2.)); break; } case 16: /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return (dsea - usea)/x; break; } } } double norm_dbar_min_ubar(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_norm_dbar_min_ubar, 0., 1., parton_npoints); } double int_norm_dbar_min_ubar(double x) { return diff_distr(param_par_s, x, scale_sqr_s); } /* the function pion_valence() gives the PION unpolarized valence quark distribution as function of x for different parametrizations */ double pion_valence(int param, double x, double scale_sqr) { double Q; Q = sqrt(scale_sqr); switch(param) { case 100: /* parametrization of Gluck, Reya, Schienbein, hep-ph 9903288 */ { double n1, n2, eud, gud, norm, test; return 1.391 /pow(x,.447) * pow(1. - x, .426); break; } } } double norm_pion_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_pion_valence, 0., 1., parton_npoints); } double int_pion_valence(double x) { return pion_valence(param_par_s, x, scale_sqr_s); } double mom_pion_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_mom_pion_valence, 0., 1., parton_npoints); } double int_mom_pion_valence(double x) { return x * pion_valence(param_par_s, x, scale_sqr_s); } /*************************************************************************/ /* The function F2_distr() gives the parton distribution function F2 for DI electron proton scattering as function of x for different parametrizations */ double F2_distr(int param, double x, double scale_sqr) { switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 1./9. * x * ( 4. * up_valence(1, x, scale_sqr) + down_valence(1, x, scale_sqr) + 11./5. * sea_distr(1, x, scale_sqr) ); break; } case 3: /* valence quark contribution parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { return 1./9. * x * ( 4. * up_valence(1, x, scale_sqr) + down_valence(1, x, scale_sqr) ); break; } } } /*************************************************************************/ /* The function gluon_distr() gives the gluon distribution as function of x for different parametrizations */ double gluon_distr(int param, double x, double scale_sqr) { double Q; Q=sqrt(scale_sqr); switch(param) { case 1: /* parametrization of Martin, Roberts, Stirling, PRD 47 (1993) 867, table II, S_o */ { double ng, Ag; ng = 5.1; Ag = 2.72; return Ag * 1./x * pow(1. - x, ng); break; } case 4: /* MRS95G */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int imode; imode=2; MRSEB(x,Q,imode,upv,dnv,usea,dsea,str,chm,bot,glu); return glu/x; break; } case 5: /* MRST98 */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRS98(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return glu/x; break; } case 8: /* CTEQ4L (LO) */ { int iparton; /* iset = 3; setctq4_(&iset); */ iparton =0; return CTQ4PDF(iparton , x, Q); break; } case 9: /* CTEQ4M (NLO MSbar) */ { int iparton; /* iset = 1; setctq4_(&iset); */ iparton =0; return CTQ4PDF(iparton , x, Q); break; } case 10: /* CTEQ5M (NLO MSbar) */ { int iparton; /* iset = 1; setctq5_(&iset); */ iparton =0; return CTQ5PDF(iparton , x, Q); break; } case 11: /* GRV (NLO MSbar) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 2; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return glu/x; break; } case 12: /* GRV (NLO DIS) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 3; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return glu/x; break; } case 13: /* GRV (LO) */ { double upv, dnv, usea, dsea, str, glu; int iset, iset2; iset = 1; iset2 = 0; GRV98PA(iset,x,scale_sqr,upv,dnv,usea,dsea,str,glu,iset2); return glu/x; break; } case 14 : /* MRST01 NLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRST2001(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return glu/x; break; } case 16 : /* MRST02 NNLO */ { double upv, dnv, usea, dsea, str, chm, bot, glu; int mode; mode=1; MRSTNNLO(x,Q,mode,upv,dnv,usea,dsea,str,chm,bot,glu); return glu/x; break; } case 17: /* CTEQ6M (NLO MSbar) */ { int iparton; /* iset = 1; setctq6_(&iset); */ iparton =0; return CTQ6PDF(iparton , x, Q); break; } } } /*************************************************************************/ /* */ /* POLARIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* The function pol_up_valence() gives the polarized valence up quark distribution as function of x for different parametrizations */ double pol_up_valence(int param, double x, double scale_sqr) { switch(param) { case 1 : /* SU(6) model by Gordon, Gosthtasbpour and Ramsey */ { return spin_dilution(param,x) * ( up_valence(param,x,scale_sqr) - 2./3. * down_valence(param,x,scale_sqr) ); break; } case 6 : /* model by Leader, Sidorov and Stamenov, using MRST98 unpolarized distribution at Q^2 = 1 GeV^2 */ { double eta_u, A_u, a_u; eta_u = 0.918; a_u = 0.250; A_u = 0.882; return eta_u * A_u * pow(x, a_u) * up_valence(6,x,1.); break; } } } double norm_pol_up_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_pol_up_valence, 0., 1., parton_npoints); } double int_pol_up_valence(double x) { return pol_up_valence(param_par_s, x, scale_sqr_s); } /* The function pol_down_valence() gives the polarized valence down quark distribution as function of x for different parametrizations */ double pol_down_valence(int param, double x, double scale_sqr) { switch(param) { case 1 : /* SU(6) model by Gordon, Gosthtasbpour and Ramsey */ { return spin_dilution(param,x) * (- 1./3. * down_valence(param,x,scale_sqr) ); break; } case 6 : /* model by Leader, Sidorov and Stamenov, using MRST98 unpolarized distribution at Q^2 = 1 GeV^2 */ { double eta_d, A_d, a_d; eta_d = -0.339; a_d = 0.231; A_d = 1.768; return eta_d * A_d * pow(x, a_d) * down_valence(6,x,1.); break; } } } double norm_pol_down_valence(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_pol_down_valence, 0., 1., parton_npoints); } double int_pol_down_valence(double x) { return pol_down_valence(param_par_s, x, scale_sqr_s); } /* The function spin_dilution() gives the spin dilution function as function of x for different unpolarized quark parametrizations (-> fit of H_0) */ double spin_dilution(int param, double x) { double h0; switch(param) { case 1 : /* fit of BSR using MRS S_o distributions */ { h0 = 0.055; break; } case 2 : /* fit of BSR using MRS D_o distributions */ { h0 = 0.075; break; } case 5 : /* fit of BSR using MRST98 distributions at 2. GeV^2 */ { h0 = 0.085; break; } case 9 : /* fit of BSR using CTEQ4L distributions at 2.6 GeV^2 */ { h0 = 0.075; break; } } return 1./(1. + h0 * pow(1. - x, 2.)/sqrt(x)); } double bjorken_sumrule(int param, int parton_npoints, double scale_sqr) { double norm; param_par_s = param; scale_sqr_s = scale_sqr; norm = 1./6. * 1.267; /* g_A/g_V = 1.267 from PDG98 */ return 1./norm * int_gauss7(int_bjorken_sumrule, 0., 1., parton_npoints); } double int_bjorken_sumrule(double x) { return 1./6. * ( pol_up_valence(param_par_s, x, scale_sqr_s) - pol_down_valence(param_par_s, x, scale_sqr_s)); } /*************************************************************************/ /* The function pol_strange_bar() gives the polarized strange_bar quark distribution as function of x for different parametrizations */ double pol_strange_bar(int param, double x, double scale_sqr) { switch(param) { case 1 : { return 0.; break; } case 6 : /* model by Leader, Sidorov and Stamenov, using MRST98 unpolarized distribution at Q^2 = 1 GeV^2 */ { double etas, As, as; etas = -0.054; as = 0.576; As = 1.6478; return etas * As * pow(x, as) * sea_distr(6,x,1.); break; } } } /* The function pol_up_bar() gives the polarized up_bar quark distribution as function of x for different parametrizations */ double pol_up_bar(int param, double x, double scale_sqr) { return pol_strange_bar(param, x, scale_sqr); } /* The function pol_down_bar() gives the polarized down_bar quark distribution as function of x for different parametrizations */ double pol_down_bar(int param, double x, double scale_sqr) { return pol_strange_bar(param, x, scale_sqr); } double norm_pol_strange_bar(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_pol_strange_bar, 0., 1., parton_npoints); } double int_pol_strange_bar(double x) { return pol_strange_bar(param_par_s, x, scale_sqr_s); } /*************************************************************************/ double pol_strange_sea(int param, double x, double scale_sqr) { return -1./3. * x * strange_sea(param, x, scale_sqr); } double norm_pol_strange_sea(int param, int parton_npoints, double scale_sqr) { param_par_s = param; scale_sqr_s = scale_sqr; return int_gauss7(int_pol_strange_sea, 0., 1., parton_npoints); } double int_pol_strange_sea(double x) { return pol_strange_sea(param_par_s, x, scale_sqr_s); }