/* This file calculates for a given quark flavour the GENERALIZED parton distribution functions */ /**************************************************************************/ /* */ /* IMPORTANT : all dimensional quantities in this program */ /* are in units GeV */ /* */ /**************************************************************************/ #include #define eps_doubleint 0.001 /* default 0.001 */ #define npoints_doubleint 50 /* default 50, when b is an integer in the profile functions */ static int flavour_s5, param_s5; static double x_s5, xi_s5, t_s5, u_s5, zeta_s5, k_perp_s5, b_s5, b_profile_val_s5, b_profile_sea_s5, scale_sqr_s5; static double brho_s; void stupid(double b_rho) { brho_s=b_rho; } /*************************************************************************/ /* */ /* UNPOLARIZED TWIST-2 GENERALIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* The function H_up_distr() returns the off-forward up quark distribution Hup(x,xi,t) of the proton as function of x in the range -1 < x < 1 */ double H_up_distr(int param, double x, double xi, double t, int factt, double b_profile_val, double b_profile_sea, double scale_sqr) { int tdep_sea; double f1p, f1n, sea_ff; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); tdep_sea = 1; switch(tdep_sea) { case 1 : /* same t-dependence as for valence part */ { sea_ff = f1p + 0.5 * f1n; break; } case 2 : /* different t-dependence for sea part */ { double b_slope = 5.; /* in GeV^-2 */ sea_ff = exp(b_slope * t / 2.); break; } } if ((param >= 1) && (param < 15)) /* factorized model q(x) * F_1(t) with valence + sea quark distribution for q(x) */ { if(x > 0.) { return ( up_valence(param, x, scale_sqr) * (f1p + 0.5 * f1n) + up_bar(param, x, scale_sqr) * sea_ff ); /* return up_valence(param, x, scale_sqr) * (f1p + 0.5 * f1n); */ } else { return - up_bar(param, -x, scale_sqr) * sea_ff; /* return 0.; */ } } if ((param >= 21) && (param < 35)) /* xi dependent distribution based on quark distribution ansatz for the double distributions */ { switch(factt) { case 1: /* factorized t-dependence */ { return offforward_distr(1,param-20,x,xi, b_profile_val,b_profile_sea,scale_sqr) * (f1p + 0.5 * f1n); break; } case 2: /* unfactorized t-dependence : regge ansatz */ { return offforward_distr_regge(1,param-20,x,xi,t, b_profile_val,b_profile_sea, scale_sqr); break; } case 3: /* Experimental exponential */ { return offforward_distr(1,param-20,x,xi, b_profile_val,b_profile_sea,scale_sqr) * exp(t*brho_s); break; } } } } /* The function H_down_distr() returns the off-forward down quark distribution Hdown(x,xi,t) of the proton as function of x in the range -1 < x < 1 */ double H_down_distr(int param, double x, double xi, double t, int factt, double b_profile_val, double b_profile_sea, double scale_sqr) { int tdep_sea; double f1p, f1n, sea_ff; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); tdep_sea = 1; switch(tdep_sea) { case 1 : /* same t-dependence as for valence part */ { sea_ff = f1p + 2. * f1n; break; } case 2 : /* different t-dependence for sea part */ { double b_slope = 5.; /* in GeV^-2 */ sea_ff = exp(b_slope * t / 2.); break; } } if ((param >= 1) && (param < 15)) /* factorized model q(x) * F_1(t) with valence + sea quark distribution for q(x) */ { if(x > 0.) { return ( down_valence(param,x,scale_sqr) * (f1p + 2. * f1n) + down_bar(param,x,scale_sqr) * sea_ff ); /* return down_valence(param,x,scale_sqr) * (f1p + 2. * f1n); */ } else { return - down_bar(param,-x,scale_sqr) * sea_ff; /* return 0.; */ } } if ((param >= 21) && (param < 35)) /* xi dependent distribution based on quark distribution ansatz for the double distributions */ { switch(factt) { case 1: /* factorized t-dependence */ { return offforward_distr(2,param-20,x,xi, b_profile_val,b_profile_sea,scale_sqr) * (f1p + 2. * f1n); break; } case 2: /* unfactorized t-dependence : regge ansatz */ { return offforward_distr_regge(2,param-20,x,xi,t, b_profile_val,b_profile_sea, scale_sqr); break; } case 3: /* Experimental exponential */ { return offforward_distr(1,param-20,x,xi, b_profile_val,b_profile_sea,scale_sqr) * exp(t*brho_s); break; } } } } /*************************************************************************/ /*************************************************************************/ /* The function E_up_distr() returns the off-forward up quark distribution Eup(x,xi,t) of the proton as function of x in the range -1 < x < 1. */ double E_up_distr(int param_e, int e_model, double x, double xi, double t, double scale_sqr, double Ju, double momu, double momuval) { if ((param_e >= 1) && (param_e < 15)) /* factorized models q(x) * F_2(t) with valence distribution for q(x) */ { double qdistr, formfactor, f2p, f2n; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); switch(e_model) { case 1 : /* valence u-quark x distribution */ { if(x > 0.) return up_valence(param_e,x,scale_sqr) * (f2p + 0.5 * f2n); else return 0.; break; } case 2 : { if(x > 0.) { return ( up_valence(param_e,x,scale_sqr) + up_bar(param_e,x,scale_sqr)) * (f2p + 0.5 * f2n); } else { return - up_bar(param_e,-x,scale_sqr) * (f2p + 0.5 * f2n); } break; } } } if ((param_e >= 21) && (param_e < 35)) /* xi dependent distribution from double distributions */ { double f2p, f2n, eval, vmpart, Au, kappau, Bu; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); kappau = 2. * form_factor_F2(1,0.) + form_factor_F2(2,0.); eval = offforward_valence_distr(1,param_e-20,x,xi,1.,scale_sqr); switch(e_model) { case 1 : /* only valence u-quark distribution */ { return kappau/2. * eval * f2p / form_factor_F2(1,0.); break; } case 2 : /* valence u-quark distribution + VM contribution */ { Au = (2. * Ju - momu) / momuval; Bu = kappau - 2. * Au; if(fabs(x) <= xi) vmpart = 3./4. * 1/xi * (1 - pow(x/xi, 2.)); else vmpart = 0.; return ( Au * eval + Bu * vmpart ) * f2p/form_factor_F2(1,0.); break; } case 3 : /* VM contribution only */ { Au = (2. * Ju - momu) / momuval; Bu = kappau - 2. * Au; if(fabs(x) <= xi) vmpart = 3./4. * 1/xi * (1 - pow(x/xi, 2.)); else vmpart = 0.; return ( Bu * vmpart ) * f2p/form_factor_F2(1,0.); break; } } } } /* The function E_down_distr() returns the off-forward down quark distribution Edown(x,xi,t) of the proton as function of x in the range -1 < x < 1. */ double E_down_distr(int param_e, int e_model, double x, double xi, double t, double scale_sqr, double Jd, double momd, double momdval) { if ((param_e >= 1) && (param_e < 15)) /* factorized model q(x) * F_2(t) with valence + sea quark distribution for q(x) */ { double qdistr, formfactor, f2p, f2n; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); switch(e_model) { case 1 : /* valence d-quark distribution */ { if(x > 0.) return down_valence(param_e,x,scale_sqr) * (f2p + 2. * f2n); else return 0.; break; } case 2 : { if(x > 0.) { return ( down_valence(param_e,x,scale_sqr) + down_bar(param_e,x,scale_sqr)) * (f2p + 2. * f2n); } else { return - down_bar(param_e,-x,scale_sqr) * (f2p + 2. * f2n); } break; } } } if ((param_e >= 21) && (param_e < 35)) /* xi dependent distribution from double distributions */ { double f2p, f2n, eval, vmpart, result, Ad, kappad, Bd; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); kappad = form_factor_F2(1,0.) + 2. * form_factor_F2(2,0.); eval = offforward_valence_distr(2,param_e-20,x,xi,1.,scale_sqr); switch(e_model) { case 1 : /* only valence d-quark distribution */ { result = kappad * eval * f2p / form_factor_F2(1,0.); break; } case 2 : /* valence d-quark distribution + VM contribution */ { Ad = (2. * Jd - momd) / momdval; Bd = kappad - Ad; if(fabs(x) <= xi) vmpart = 3./4. * 1/xi * (1 - pow(x/xi, 2.)); else vmpart = 0.; result = ( Ad * eval + Bd * vmpart ) * f2p/form_factor_F2(1,0.); break; } case 3 : /* VM contribution only */ { Ad = (2. * Jd - momd) / momdval; Bd = kappad - Ad; if(fabs(x) <= xi) vmpart = 3./4. * 1/xi * (1 - pow(x/xi, 2.)); else vmpart = 0.; result = ( Bd * vmpart ) * f2p/form_factor_F2(1,0.); break; } } return result; } } /*************************************************************************/ /*************************************************************************/ /* the function offforward_distr() gives the ofpd for the quark of a given flavour as function of x in the range -1 < x < 1. */ double offforward_distr(int flavour, int param, double x, double xi, double b_profile_val, double b_profile_sea, double scale_sqr) { int dd; double ofpd; dd = 1; /* default = 1 */ switch(dd) { case 1 : /* parametrization through symmetric Double Distributions (DD in beta and alpha) with variable profile function (parameter b) */ { x_s5 = x; xi_s5 = xi; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; b_profile_sea_s5 = b_profile_sea; scale_sqr_s5 = scale_sqr; if(x >= xi) { ofpd = int_gauss5(int_symm_double_distr, - (1. - x)/(1. + xi), (1. - x)/(1. - xi), npoints_doubleint); } else { if( (-xi < x) && (x < xi) ) { ofpd = int_gauss5(int_symm_double_distr, - (1. - x)/(1. + xi), x/xi - eps_doubleint, npoints_doubleint) - int_gauss5(int_symm_double_antidistr,-(1. + x)/(1. + xi), -x/xi - eps_doubleint, npoints_doubleint); } else { ofpd = - int_gauss5(int_symm_double_antidistr, - (1. + x)/(1. + xi), (1. + x)/(1. - xi), npoints_doubleint); } } break; } case 2 : /* parametrization through NFPD and Double Distributions (DD in x and y) with asymptotic profile function (fixed b = 1) */ { double X, zeta; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x >= xi) { ofpd = 1./(1. + xi) * nonforward_distr(flavour,param,X,zeta,scale_sqr); } else { if( (-xi < x) && (x < xi) ) { ofpd = 1./(1. + xi) * ( nonforward_distr(flavour,param,X,zeta,scale_sqr) - nonforward_antidistr(flavour,param,zeta - X,zeta,scale_sqr) ); } else { ofpd = 1./(1. + xi) * (-1.) * nonforward_antidistr(flavour,param,zeta - X,zeta,scale_sqr); } } break; } } return ofpd; } double offforward_valence_distr(int flavour, int param, double x, double xi, double b_profile_val, double scale_sqr) { int dd; double ofpd; dd = 1; /* default = 1 */ switch(dd) { case 1 : /* parametrization through symmetric Double Distributions (DD in beta and alpha) with variable profile function (parameter b) */ { x_s5 = x; xi_s5 = xi; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; scale_sqr_s5 = scale_sqr; if(x >= xi) { ofpd = int_gauss5(int_symm_double_valence_distr, - (1. - x)/(1. + xi), (1. - x)/(1. - xi), 100); } else { if( (-xi < x) && (x < xi) ) { ofpd = int_gauss5(int_symm_double_valence_distr, - (1. - x)/(1. + xi), x/xi - 0.00001, 100); } else ofpd = 0.; } break; } case 2 : /* parametrization through NFPD and Double Distributions (DD in x and y) with asymptotic profile function (fixed b = 1) */ { double X, zeta; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x > -xi) { ofpd = 1./(1. + xi) * nonforward_valence_distr(flavour,param,X,zeta,scale_sqr); } else { ofpd = 0.; } break; } } return ofpd; } /* the function offforward_distr_regge() gives the gpd for the quark of a given flavour as function of x in the range -1 < x < 1 using a non-factorizable Regge-type ansatz in t. */ double offforward_distr_regge(int flavour, int param, double x, double xi, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { double gpd; /* parametrization through symmetric Double Distributions (DD in beta and alpha) with variable profile function (parameter b) */ x_s5 = x; xi_s5 = xi; t_s5 = t; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; b_profile_sea_s5 = b_profile_sea; scale_sqr_s5 = scale_sqr; if(x >= xi) { gpd = int_gauss5(int_symm_double_distr_regge, - (1. - x)/(1. + xi), (1. - x)/(1. - xi), npoints_doubleint); } else { if( (-xi < x) && (x < xi) ) { gpd = int_gauss5(int_symm_double_distr_regge, -(1. - x)/(1. + xi), x/xi - eps_doubleint, npoints_doubleint) - int_gauss5(int_symm_double_antidistr_regge, -(1. + x)/(1. + xi), -x/xi - eps_doubleint, npoints_doubleint); } else { gpd = - int_gauss5(int_symm_double_antidistr_regge, - (1. + x)/(1. + xi), (1. + x)/(1. - xi), npoints_doubleint); } } return gpd; } double offforward_valence_distr_regge(int flavour, int param, double x, double xi, double t, double b_profile_val,double scale_sqr) { double gpd; /* parametrization through symmetric Double Distributions (DD in beta and alpha) with variable profile function (parameter b) */ x_s5 = x; xi_s5 = xi; t_s5 = t; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; scale_sqr_s5 = scale_sqr; if(x >= xi) { gpd = int_gauss5(int_symm_double_valence_distr_regge, - (1. - x)/(1. + xi), (1. - x)/(1. - xi), npoints_doubleint); } else { if( (-xi < x) && (x < xi) ) { gpd = int_gauss5(int_symm_double_valence_distr_regge, - (1. - x)/(1. + xi), x/xi - 0.00001, npoints_doubleint); } else { gpd = 0.0; } } return gpd; } /*************************************************************************/ /*************************************************************************/ /* The functions symm_double_distr() are the symmetric double distributions (in arguments beta and alpha) */ double int_symm_double_distr(double alpha) { return symm_double_distr(flavour_s5,param_s5, x_s5 - xi_s5 * alpha, alpha, b_profile_val_s5,b_profile_sea_s5,scale_sqr_s5); } double symm_double_distr(int flavour, int param, double beta, double alpha, double b_profile_val, double b_profile_sea, double scale_sqr) { double symm_dd; if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in antidistr", alpha, beta); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr) + symm_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr); break; } case 2: /* down quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr) + symm_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr); break; } } return symm_dd; } double int_symm_double_antidistr(double alpha) { return symm_double_antidistr(flavour_s5,param_s5, - x_s5 - xi_s5 * alpha, alpha, b_profile_sea_s5, scale_sqr_s5); } double symm_double_antidistr(int flavour, int param, double beta, double alpha, double b_profile_sea, double scale_sqr) { double symm_dd; if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in antidistr", alpha, beta); switch(flavour) { case 1 : /* up quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr); break; } case 2 : /* down quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr); break; } } return symm_dd; } double int_symm_double_valence_distr(double alpha) { return symm_double_valence_distr(flavour_s5, param_s5, x_s5 - xi_s5 * alpha, alpha, b_profile_val_s5, scale_sqr_s5); } double symm_double_valence_distr(int flavour, int param, double beta, double alpha, double b_profile_val, double scale_sqr) { double symm_dd; if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in valence distr", alpha, beta); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr); break; } case 2: /* down quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr); break; } } return symm_dd; } double symm_profile_function(double beta, double alpha, double b_profile) { double norm; norm = 1./pow(2., 2.* b_profile + 1.) * gamma_function(2.* b_profile + 2.) / pow(gamma_function(b_profile + 1.), 2.); return norm * pow(pow(1. - fabs(beta), 2.) - pow(alpha, 2.),b_profile) / pow(1. - fabs(beta), 2.* b_profile + 1.); } /* The functions symm_double_distr_regge() are the symmetric double distributions in the Regge type ansatz (in arguments beta and alpha) */ double int_symm_double_distr_regge(double alpha) { return symm_double_distr_regge(flavour_s5,param_s5, x_s5 - xi_s5 * alpha, alpha, t_s5, b_profile_val_s5,b_profile_sea_s5, scale_sqr_s5); } double symm_double_distr_regge(int flavour, int param, double beta, double alpha, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { double symm_dd, alphap; alphap = 0.8; /* Regge trajectory slope */ if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in valence distr", alpha, beta); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = ( symm_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr) + symm_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr) ) * pow(1./fabs(beta), alphap * t); break; } case 2: /* down quark distribution */ { symm_dd = ( symm_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr) + symm_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr) ) * pow(1./fabs(beta), alphap * t); break; } } return symm_dd; } double int_symm_double_antidistr_regge(double alpha) { return symm_double_antidistr_regge(flavour_s5,param_s5, - x_s5 - xi_s5 * alpha, alpha, t_s5, b_profile_sea_s5, scale_sqr_s5); } double symm_double_antidistr_regge(int flavour, int param, double beta, double alpha, double t, double b_profile_sea, double scale_sqr) { double symm_dd, alphap; alphap = 0.8; /* Regge trajectory slope */ if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in valence distr", alpha, beta); switch(flavour) { case 1 : /* up quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr) * pow(1./fabs(beta), alphap * t); break; } case 2 : /* down quark distribution */ { symm_dd = symm_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr) * pow(1./fabs(beta), alphap * t); break; } } return symm_dd; } double int_symm_double_valence_distr_regge(double alpha) { return symm_double_valence_distr_regge(flavour_s5,param_s5, x_s5 - xi_s5 * alpha,alpha,t_s5, b_profile_val_s5, scale_sqr_s5); } double symm_double_valence_distr_regge(int flavour, int param, double beta, double alpha, double t, double b_profile_val, double scale_sqr) { double symm_dd, alphap; alphap = 0.8; /* Regge trajectory slope */ if(beta <= 0.) printf("\n alpha = %f, beta = %f, argument 0 or negative in valence distr", alpha, beta); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = ( symm_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr) ) * pow(1./fabs(beta), alphap * t); break; } case 2: /* down quark distribution */ { symm_dd = ( symm_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr) ) * pow(1./fabs(beta), alphap * t); break; } } return symm_dd; } /*************************************************************************/ /*************************************************************************/ double nonforward_distr(int flavour, int param, double x, double zeta, double scale_sqr) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; scale_sqr_s5 = scale_sqr; if (x >= zeta) return int_gauss5(int_double_distr, 0., (1. - x)/(1. - zeta), npoints_doubleint); else return int_gauss5(int_double_distr, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_antidistr(int flavour, int param, double x, double zeta, double scale_sqr) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; scale_sqr_s5 = scale_sqr; if (x >= zeta) return int_gauss5(int_double_antidistr,0.,(1. - x)/(1. - zeta), npoints_doubleint); else return int_gauss5(int_double_antidistr, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_valence_distr(int flavour, int param, double x, double zeta, double scale_sqr) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; scale_sqr_s5 = scale_sqr; if (x >= zeta) return int_gauss5(int_double_valence_distr,0.,(1. - x)/(1. - zeta), npoints_doubleint); else { return int_gauss5(int_double_valence_distr, 0., x/zeta - eps_doubleint, npoints_doubleint); } } /*************************************************************************/ /*************************************************************************/ double int_double_distr(double y) { return double_distr(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, scale_sqr_s5); } double double_distr(int flavour, int param, double x, double y, double scale_sqr) { double qdistr; if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr) + up_bar(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr) + down_bar(param, x, scale_sqr); break; } } return profile_function(x, y) * qdistr; } double int_double_antidistr(double y) { return double_antidistr(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, scale_sqr_s5); } double double_antidistr(int flavour, int param, double x, double y, double scale_sqr) { double qdistr; if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1 : /* up quark distribution */ { qdistr = up_bar(param, x, scale_sqr); break; } case 2 : /* down quark distribution */ { qdistr = down_bar(param, x, scale_sqr); break; } } return profile_function(x, y) * qdistr; } double int_double_valence_distr(double y) { return double_valence_distr(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, scale_sqr_s5); } double double_valence_distr(int flavour, int param, double x, double y, double scale_sqr) { double qdistr; if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr); break; } } return profile_function(x, y) * qdistr; } double profile_function(double x, double y) { int profile; double asympt; asympt = 6. * y * (1. - x - y) / pow(1. - x,3.); profile = 1; switch(profile) { case 1 : /* asymptotic profile */ { return asympt; break; } case 2 : /* asymptotic profile + 2nd order Gegenbauer term (for test only!) */ { double a2; a2 = -0.58; return asympt * (1. + a2 * gegenbauer_1p5(2, 2.*y/(1. - x) - 1.)); break; } case 3 : /* b-profile, with b = 2 */ { double normb, b; b = 2.; normb = 30.; return normb * pow(y * (1. - x - y), b)/pow(1. - x, 2.*b + 1.); break; } } } /*************************************************************************/ /* */ /* D-TERM CONTRIBUTION TO THE ISOSCALAR GPDs H and E */ /* */ /*************************************************************************/ /* D-term contribution to GPD of each flavor. The argument alpha runs from -1 -> 1. To obtain D-term to isoscalar GPD : multiply by N_f = 3. */ double spd_dterm(double alpha) { double norm, d3, d5, term1, term3, term5, N_flavors; term1 = gegenbauer_1p5(1, alpha); term3 = gegenbauer_1p5(3, alpha); term5 = gegenbauer_1p5(5, alpha); N_flavors = 3.; norm = 4.0; d3 = 0.3; d5 = 0.1; return - 1./N_flavors * norm * (1. - pow(alpha, 2.)) * (term1 + d3 * term3 + d5 * term5); } /*************************************************************************/ /* */ /* UNPOLARIZED TWIST-3 GENERALIZED PARTON DISTRIBUTIONS */ /* IN WANDZURA-WILCZEK APPROXIMATION */ /* */ /*************************************************************************/ /* The function H_up_twist3ww_distr() returns the part of the off-forward twist-3 WW up quark distribution of the proton which can be expressed as a derivative of the twist-2 GPD. Function of u in the range -1 < u < 1 */ double H_up_twist3ww_distr(int param, double u, double xi, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { int tdep_sea; double f1p, f1n, sea_ff; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); tdep_sea = 1; switch(tdep_sea) { case 1 : /* same t-dependence as for valence part */ { sea_ff = f1p + 0.5 * f1n; break; } case 2 : /* different t-dependence for sea part */ { double b_slope = 5.; /* in GeV^-2 */ sea_ff = exp(b_slope * t / 2.); break; } } if ((param >= 21) && (param < 35)) /* xi dependent twist-3 distribution based on quark distribution ansatz for the double distributions */ { return offforward_twist3ww_distr(1,param-20,u,xi, b_profile_val,b_profile_sea,scale_sqr) * (f1p + 0.5 * f1n); } else return 0.; } /* The function H_down_twist3ww_distr() returns the off-forward twist-3 WW down quark distribution of the proton which can be expressed as a derivative of the twist-2 GPD. Function of u in the range -1 < u < 1 */ double H_down_twist3ww_distr(int param, double u, double xi, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { int tdep_sea; double f1p, f1n, sea_ff; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); tdep_sea = 1; switch(tdep_sea) { case 1 : /* same t-dependence as for valence part */ { sea_ff = f1p + 2. * f1n; break; } case 2 : /* different t-dependence for sea part */ { double b_slope = 5.; /* in GeV^-2 */ sea_ff = exp(b_slope * t / 2.); break; } } if ((param >= 21) && (param < 35)) /* xi dependent twist-3 distribution based on quark distribution ansatz for the double distributions */ { return offforward_twist3ww_distr(2,param-20,u,xi, b_profile_val,b_profile_sea,scale_sqr) * (f1p + 2. * f1n); } else return 0.; } /*************************************************************************/ /*************************************************************************/ /* The function E_up_twist3ww_distr() returns the off-forward twist-3 WW up quark distribution of the proton which can be expressed as a derivative of the twist-2 GPD. Function of u in the range -1 < u < 1. */ double E_up_twist3ww_distr(int param_e, double u, double xi, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { if ((param_e >= 21) && (param_e < 35)) /* xi dependent twist-3 distribution based on quark distribution ansatz for the double distributions */ { double f2p, f2n, result; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); result = offforward_twist3ww_distr(1,param_e - 20,u,xi, b_profile_val,b_profile_sea,scale_sqr) * (f2p + 0.5 * f2n); return result; } } /* The function E_down_twist3ww_distr() returns the off-forward twist-3 WW down quark distribution of the proton which can be expressed as a derivative of the twist-2 GPD. Function of u in the range -1 < u < 1. */ double E_down_twist3ww_distr(int param_e, double u, double xi, double t, double b_profile_val, double b_profile_sea, double scale_sqr) { if ((param_e >= 21) && (param_e < 35)) /* xi dependent twist-3 distribution based on quark distribution ansatz for the double distributions */ { double f2p, f2n, result; f2p = form_factor_F2(1,-t); f2n = form_factor_F2(2,-t); result = offforward_twist3ww_distr(2,param_e - 20,u,xi, b_profile_val,b_profile_sea,scale_sqr) * (f2p + 2. * f2n); return result; } } /*************************************************************************/ /*************************************************************************/ /* the function offforward_twist3ww_distr() gives the part of the twist-3 WW GPD which can be expressed in terms of derivatives of the twist-2 GPDs. Function for the quark of a given flavor and of u in the range -1 < u < 1. Parametrization through symmetric Double Distributions (DD in beta and alpha) with variable profile function (parameter b). */ double offforward_twist3ww_distr(int flavour, int param, double u, double xi, double b_profile_val, double b_profile_sea, double scale_sqr) { double ofpd; u_s5 = u; xi_s5 = xi; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; b_profile_sea_s5 = b_profile_sea; scale_sqr_s5 = scale_sqr; if(u >= xi) { ofpd = int_gauss5(int_symm_deriv_double_distr, - (1. - u)/(1. + xi), (1. - u)/(1. - xi), npoints_doubleint); } else { if( (-xi < u) && (u < xi) ) { ofpd = int_gauss5(int_symm_deriv_double_distr,-(1. - u)/(1. + xi), u/xi - eps_doubleint, npoints_doubleint) - int_gauss5(int_symm_deriv_double_antidistr,-(1. + u)/(1. + xi), -u/xi - eps_doubleint, npoints_doubleint); } else { ofpd = - int_gauss5(int_symm_deriv_double_antidistr, - (1. + u)/(1. + xi), (1. + u)/(1. - xi), npoints_doubleint); } } return ofpd; } double offforward_twist3ww_valence_distr(int flavour, int param, double u, double xi, double b_profile_val,double scale_sqr) { double ofpd; u_s5 = u; xi_s5 = xi; flavour_s5 = flavour; param_s5 = param; b_profile_val_s5 = b_profile_val; scale_sqr_s5 = scale_sqr; if(u >= xi) { ofpd = int_gauss5(int_symm_deriv_double_valence_distr, - (1. - u)/(1. + xi), (1. - u)/(1. - xi), npoints_doubleint); } else { if( (-xi < u) && (u < xi) ) { ofpd = int_gauss5(int_symm_deriv_double_valence_distr, - (1. - u)/(1. + xi), u/xi - eps_doubleint, npoints_doubleint); } else ofpd = 0.; } return ofpd; } /*************************************************************************/ /*************************************************************************/ /* The functions symm_deriv_double_distr() are the symmetric double distributions with a derivative to alpha (in arguments beta and alpha) */ double int_symm_deriv_double_distr(double alpha) { return (u_s5 - xi_s5 * alpha) / xi_s5 * symm_deriv_double_distr(flavour_s5,param_s5,u_s5 - xi_s5 * alpha, alpha,b_profile_val_s5, b_profile_sea_s5, scale_sqr_s5); } double symm_deriv_double_distr(int flavour,int param,double beta,double alpha, double b_profile_val, double b_profile_sea, double scale_sqr) { double symm_dd; if(beta <= 0) printf("\n argument 0 or negative!!!"); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr) + symm_deriv_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr); break; } case 2: /* down quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr) + symm_deriv_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr); break; } } return symm_dd; } double int_symm_deriv_double_antidistr(double alpha) { return (- u_s5 - xi_s5 * alpha) / xi_s5 * symm_deriv_double_antidistr(flavour_s5,param_s5, - u_s5 - xi_s5 * alpha, alpha, b_profile_sea_s5, scale_sqr_s5); } double symm_deriv_double_antidistr(int flavour, int param, double beta, double alpha, double b_profile_sea, double scale_sqr) { double symm_dd; if(beta <= 0) printf("\n argument 0 or negative!!!"); switch(flavour) { case 1 : /* up quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_sea) * up_bar(param, beta, scale_sqr); break; } case 2 : /* down quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_sea) * down_bar(param, beta, scale_sqr); break; } } return symm_dd; } double int_symm_deriv_double_valence_distr(double alpha) { return (u_s5 - xi_s5 * alpha) / xi_s5 * symm_deriv_double_valence_distr(flavour_s5, param_s5, u_s5 - xi_s5 * alpha, alpha, b_profile_val_s5, scale_sqr_s5); } double symm_deriv_double_valence_distr(int flavour, int param, double beta, double alpha, double b_profile_val, double scale_sqr) { double symm_dd; if(beta <= 0) printf("\n argument 0 or negative!!!"); switch(flavour) { case 1: /* up quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_val) * up_valence(param, beta, scale_sqr); break; } case 2: /* down quark distribution */ { symm_dd = symm_deriv_profile_function(beta, alpha, b_profile_val) * down_valence(param, beta, scale_sqr); break; } } return symm_dd; } double symm_deriv_profile_function(double beta, double alpha, double b_profile) { double norm; norm = 1./pow(2., 2.* b_profile + 1.) * gamma_function(2.* b_profile + 2.) / pow(gamma_function(b_profile + 1.), 2.); return norm * (- 2. * b_profile * alpha) * pow(pow(1. - fabs(beta), 2.) - pow(alpha, 2.), b_profile - 1.) / pow(1. - fabs(beta), 2.* b_profile + 1.); } /*************************************************************************/ /* */ /* K_PERP DEPENDENT UNPOLARIZED GENERALIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* The function H_up_distr_kperp() returns the k_perp dependent GPD for the up quark in the proton : Hup(x,k_perp,xi,t) */ double H_up_distr_kperp(int param,double x,double k_perp, double xi, double t) { double f1p, f1n, result; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); result = offforward_distr_kperp(1,param-40,x,k_perp,xi) * (f1p + 0.5 * f1n); return result; } /* The function H_down_distr_kperp() returns the k_perp dependent GPD for the down quark in the proton : Hdown(x,k_perp,xi,t) */ double H_down_distr_kperp(int param,double x,double k_perp,double xi,double t) { double f1p, f1n, result; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); result = offforward_distr_kperp(2,param-40,x,k_perp,xi) * (f1p + 2. * f1n); return result; } /*************************************************************************/ /*************************************************************************/ double offforward_distr_kperp(int flavour, int param, double x, double k_perp, double xi) { double X, zeta, ofpd; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(xi > 0.01) /* for numerical reasons */ { if( (x > 0.) && (fabs(x - xi) < (2. * eps_doubleint)) ) x = xi; if( (x < 0.) && (fabs(x + xi) < (2. * eps_doubleint)) ) x = -xi; } if(x >= xi ) { ofpd = 1./(1. + xi) * nonforward_distr_kperp(flavour, param, X, k_perp, zeta); } else { if( (-xi < x) && (x < xi) ) { ofpd = 1./(1. + xi) * ( nonforward_distr_kperp(flavour,param,X,k_perp,zeta) - nonforward_antidistr_kperp(flavour,param,zeta - X,k_perp,zeta) ); } else { ofpd = 1./(1. + xi) * (-1.) * nonforward_antidistr_kperp(flavour,param,zeta - X,k_perp,zeta); } } return ofpd; } double offforward_valence_distr_kperp(int flavour, int param, double x, double k_perp, double xi) { double X, zeta, ofpd; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x > -xi) { ofpd = 1./(1. + xi) * nonforward_valence_distr_kperp(flavour, param, X, k_perp, zeta); } else { ofpd = 0.; } return ofpd; } /*************************************************************************/ /*************************************************************************/ double nonforward_distr_kperp(int flavour, int param, double x, double k_perp, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; k_perp_s5 = k_perp; param_s5=param; if (x >= zeta) return int_gauss5(int_double_distr_kperp, 0., (1. - x)/(1. - zeta),25); else return int_gauss5(int_double_distr_kperp, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_antidistr_kperp(int flavour, int param, double x, double k_perp, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; k_perp_s5 = k_perp; param_s5=param; if (x >= zeta) return int_gauss5(int_double_antidistr_kperp,0.,(1. - x)/(1. - zeta),25); else return int_gauss5(int_double_antidistr_kperp, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_valence_distr_kperp(int flavour, int param, double x, double k_perp, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; k_perp_s5 = k_perp; param_s5=param; if (x >= zeta) return int_gauss5(int_double_valence_distr_kperp, 0.,(1. - x)/(1. - zeta),25); else return int_gauss5(int_double_valence_distr_kperp, 0., x/zeta - eps_doubleint, 25); } /*************************************************************************/ /*************************************************************************/ double int_double_distr_kperp(double y) { return double_distr_kperp(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, k_perp_s5); } double double_distr_kperp(int flavour, int param, double x, double y, double k_perp) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr) + up_bar(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr) + down_bar(param, x, scale_sqr); break; } } gaussian = 4. * pow(PI,2.)/( SIGMA_NU * y * (1. - x - y) ) * exp(- pow(k_perp,2.)/(2. * SIGMA_NU * y * (1. - x - y)) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } double int_double_antidistr_kperp(double y) { return double_antidistr_kperp(flavour_s5, param_s5, x_s5 - zeta_s5 * y,y,k_perp_s5); } double double_antidistr_kperp(int flavour, int param, double x, double y, double k_perp) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ switch(flavour) { case 1 : /* up quark distribution */ { qdistr = up_bar(param, x, scale_sqr); break; } case 2 : /* down quark distribution */ { qdistr = down_bar(param, x, scale_sqr); break; } } gaussian = 4. * pow(PI,2.)/( SIGMA_NU * y * (1. - x - y) ) * exp(- pow(k_perp,2.)/(2. * SIGMA_NU * y * (1. - x - y)) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } double int_double_valence_distr_kperp(double y) { return double_valence_distr_kperp(flavour_s5, param_s5, x_s5 - zeta_s5 * y,y,k_perp_s5); } double double_valence_distr_kperp(int flavour, int param, double x, double y, double k_perp) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr); break; } } gaussian = 4. * pow(PI,2.)/( SIGMA_NU * y * (1. - x - y) ) * exp(- pow(k_perp,2.)/(2. * SIGMA_NU * y * (1. - x - y)) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } /*****************************************************************************/ /* */ /* UNPOLARIZED GENERALIZED PARTON DISTRIBUTIONS IN IMPACT PARAMETER (b) SPACE*/ /* */ /*****************************************************************************/ double H_up_distr_bspace(int param, double x, double b, double xi, double t) { double f1p, f1n, result; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); result = offforward_distr_bspace(1,param-40,x,b,xi) * (f1p + 0.5 * f1n); return result; } double H_down_distr_bspace(int param, double x, double b, double xi, double t) { double f1p, f1n, result; f1p = form_factor_F1(1,-t); f1n = form_factor_F1(2,-t); result = offforward_distr_bspace(2,param-40,x,b,xi) * (f1p + 2. * f1n); return result; } /*************************************************************************/ /*************************************************************************/ double offforward_distr_bspace(int flavour, int param, double x, double b, double xi) { double X, zeta, ofpd; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x >= xi) { ofpd = 1./(1. + xi) * nonforward_distr_bspace(flavour,param,X,b,zeta); } else { if( (-xi < x) && (x < xi) ) { ofpd = 1./(1. + xi) * ( nonforward_distr_bspace(flavour,param,X,b,zeta) - nonforward_antidistr_bspace(flavour,param,zeta - X,b,zeta) ); } else { ofpd = 1./(1. + xi) * (-1.) * nonforward_antidistr_bspace(flavour,param,zeta - X, b, zeta); } } return ofpd; } double offforward_valence_distr_bspace(int flavour, int param, double x, double b, double xi) { double X, zeta, ofpd; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x > -xi) { ofpd = 1./(1. + xi) * nonforward_valence_distr_bspace(flavour, param, X, b, zeta); } else { ofpd = 0.; } return ofpd; } /*************************************************************************/ /*************************************************************************/ double nonforward_distr_bspace(int flavour, int param, double x, double b, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5=param; b_s5 = b; if (x >= zeta) return int_gauss5(int_double_distr_bspace,0., (1. - x)/(1. - zeta),25); else return int_gauss5(int_double_distr_bspace, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_antidistr_bspace(int flavour, int param, double x,double b,double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5=param; b_s5 = b; if (x >= zeta) return int_gauss5(int_double_antidistr_bspace, 0.,(1. - x)/(1. - zeta),25); else return int_gauss5(int_double_antidistr_bspace, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_valence_distr_bspace(int flavour, int param, double x, double b, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5=param; b_s5 = b; if (x >= zeta) return int_gauss5(int_double_valence_distr_bspace, 0.,(1. - x)/(1. - zeta),25); else return int_gauss5(int_double_valence_distr_bspace, 0., x/zeta - eps_doubleint, 25); } /*************************************************************************/ /*************************************************************************/ double int_double_distr_bspace(double y) { return double_distr_bspace(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, b_s5); } double double_distr_bspace(int flavour,int param, double x, double y, double b) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr) + up_bar(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr) + down_bar(param, x, scale_sqr); break; } } gaussian = 2.*PI * exp(- 0.5 * SIGMA_NU * pow(b,2.) * y * (1. - x - y) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } double int_double_antidistr_bspace(double y) { return double_antidistr_bspace(flavour_s5, param_s5, x_s5 - zeta_s5 * y,y,b_s5); } double double_antidistr_bspace(int flavour, int param, double x, double y, double b) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ switch(flavour) { case 1 : /* up quark distribution */ { qdistr = up_bar(param, x, scale_sqr); break; } case 2 : /* down quark distribution */ { qdistr = down_bar(param, x, scale_sqr); break; } } gaussian = 2.*PI * exp(- 0.5 * SIGMA_NU * pow(b,2.) * y * (1. - x - y) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } double int_double_valence_distr_bspace(double y) { return double_valence_distr_bspace(flavour_s5, param_s5, x_s5 - zeta_s5 * y,y,b_s5); } double double_valence_distr_bspace(int flavour, int param, double x, double y, double b) { double qdistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* up quark distribution */ { qdistr = up_valence(param, x, scale_sqr); break; } case 2: /* down quark distribution */ { qdistr = down_valence(param, x, scale_sqr); break; } } gaussian = 2.*PI * exp(- 0.5 * SIGMA_NU * pow(b,2.) * y * (1. - x - y) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qdistr/norm * gaussian; } /*************************************************************************/ /* */ /* POLARIZED GENERALIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* The function Hhel_up_distr() returns the off-forward up quark polarized distribution Htilde(x,xi,t) of the proton as function of x in the range -1 < x < 1. */ double Hhel_up_distr(int param, double x, double xi, double t, double scale_sqr) { double formfactor; if((param == 1) || (param == 6)) /* factorized model qhel(x) * G_A(t) with only polarized valence quark distribution for qhel(x) */ { double qheldistr; if ( x <= 0.) return 0.; else { qheldistr = pol_up_valence(param, x, scale_sqr); formfactor = 1/g_AXIAL * form_factor_G_A(-t); return qheldistr * formfactor; } } if ((param == 21) || (param == 26)) /* xi dependent distribution based on polarized quark distribution ansatz for the double distributions */ { double qheldistr; qheldistr = offforward_heldistr(1, param - 20, x, xi, scale_sqr); formfactor = 1/g_AXIAL * form_factor_G_A(-t); return qheldistr * formfactor; } } /* The function Hhel_down_distr() returns the off-forward down quark polarized distribution Htilde(x,xi,t) of the proton as function of x in the range -1 < x < 1. */ double Hhel_down_distr(int param, double x, double xi, double t, double scale_sqr) { double formfactor; if ((param == 1) || (param == 6)) /* factorized model qhel(x) * G_A(t) with only polarized valence quark distribution for qhel(x) */ { double qheldistr; if ( x <= 0.) return 0.; else { qheldistr = pol_down_valence(param, x, scale_sqr); formfactor = 1/g_AXIAL * form_factor_G_A(-t); return qheldistr * formfactor; } } if ((param == 21) || (param == 26)) /* xi dependent distribution based on polarized quark distribution ansatz for the double distributions */ { double qheldistr; qheldistr = offforward_heldistr(2, param - 20, x, xi, scale_sqr); formfactor = 1/g_AXIAL * form_factor_G_A(-t); return qheldistr * formfactor; } } double Hhel_strange_distr(int param, double x, double xi, double t, double scale_sqr) { double formfactor; if ((param == 1) || (param == 6)) /* factorized model qhel(x) * G_A(t) with polarized quark distribution for qhel(x) */ { double qheldistr; if ( x <= 0.) return 0.; else return 0.; } if ((param == 21) || (param == 26)) /* xi dependent distribution based on polarized quark distribution ansatz for the double distributions */ { double qheldistr; qheldistr = offforward_heldistr(3, param - 20, x, xi, scale_sqr); formfactor = 1/g_AXIAL * form_factor_G_A(-t); return qheldistr * formfactor; } } /* The function offforward_heldistr() returns the OFF-forward polarized parton distribution (at t = 0) for a given quark flavour in the proton. */ double offforward_heldistr(int flavour, int param, double x, double xi, double scale_sqr) { double X, zeta, ofpd; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if(x >= xi) { ofpd = 1./(1. + xi) * nonforward_heldistr(flavour,param,X,zeta,scale_sqr); } else { if( (-xi < x) && (x < xi) ) { ofpd = 1./(1. + xi) * ( nonforward_heldistr(flavour,param,X,zeta,scale_sqr) + nonforward_helantidistr(flavour,param,zeta - X,zeta,scale_sqr) ); } else { ofpd = 1./(1. + xi) * nonforward_helantidistr(flavour,param,zeta - X,zeta,scale_sqr); } } return ofpd; } /*************************************************************************/ /*************************************************************************/ double nonforward_heldistr(int flavour, int param, double x, double zeta, double scale_sqr) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; scale_sqr_s5 = scale_sqr; if (x >= zeta) return int_gauss5(int_double_heldistr, 0., (1. - x)/(1. - zeta),40); else return int_gauss5(int_double_heldistr, 0., x/zeta - eps_doubleint, npoints_doubleint); } double nonforward_helantidistr(int flavour, int param, double x, double zeta, double scale_sqr) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; scale_sqr_s5 = scale_sqr; if (x >= zeta) return int_gauss5(int_double_helantidistr,0.,(1. - x)/(1. - zeta),40); else return int_gauss5(int_double_helantidistr, 0., x/zeta - eps_doubleint, npoints_doubleint); } /*************************************************************************/ /*************************************************************************/ double int_double_heldistr(double y) { return double_heldistr(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, scale_sqr_s5); } double double_heldistr(int flavour, int param, double x, double y, double scale_sqr) { double qheldistr, norm; if(x == 0) printf("\n argument 0 in double hel distribution!"); switch(flavour) { case 1: /* polarized up quark helicity distribution */ { qheldistr = pol_up_valence(param, x, scale_sqr) + pol_up_bar(param, x, scale_sqr); break; } case 2: /* polarized down quark distribution */ { qheldistr = pol_down_valence(param, x, scale_sqr) + pol_down_bar(param, x, scale_sqr); break; } case 3: /* polarized strange quark distribution */ { qheldistr = pol_strange_bar(param, x, scale_sqr); break; } } norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qheldistr/norm; } double int_double_helantidistr(double y) { return double_helantidistr(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, scale_sqr_s5); } double double_helantidistr(int flavour, int param, double x, double y, double scale_sqr) { double qheldistr, norm; if(x == 0) printf("\n argument 0 in double hel distribution!"); switch(flavour) { case 1: /* polarized up quark helicity distribution */ { qheldistr = pol_up_bar(param, x, scale_sqr); break; } case 2: /* polarized down quark distribution */ { qheldistr = pol_down_bar(param, x, scale_sqr); break; } case 3: /* polarized strange quark distribution */ { qheldistr = pol_strange_bar(param, x, scale_sqr); break; } } norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qheldistr/norm; } /*************************************************************************/ /* */ /* K_PERP DEPENDENT POLARIZED GENERALIZED PARTON DISTRIBUTIONS */ /* */ /*************************************************************************/ /* The function Hhel_up_distr_kperp() returns the k_perp dependent helicity OFPD for the up quark in the proton : Hhelup(x,k_perp,xi,t) */ double Hhel_up_distr_kperp(int param, double x, double k_perp, double xi, double t) { double formfactor; formfactor = 1/g_AXIAL * form_factor_G_A(-t); return formfactor * offforward_heldistr_kperp(1,param-40, x, k_perp, xi); } /* The function Hhel_down_distr_kperp() returns the k_perp dependent helicity OFPD for the down quark in the proton : Hheldown(x,k_perp,xi,t) */ double Hhel_down_distr_kperp(int param, double x, double k_perp, double xi, double t) { double formfactor; formfactor = 1/g_AXIAL * form_factor_G_A(-t); return formfactor * offforward_heldistr_kperp(2,param-40, x, k_perp, xi); } /*************************************************************************/ /*************************************************************************/ /* The function offforward_heldistr_kperp() returns the k_perp dependent OFF-forward valence polarized parton distribution (at t=0) for a given quark flavour in the proton. */ double offforward_heldistr_kperp(int flavour, int param, double x, double k_perp, double xi) { double X, zeta; X = (x + xi)/(1. + xi); zeta = 2. * xi / (1. + xi); if (x <= (-xi)) return 0.; else return 1./(1. + xi) * nonforward_heldistr_kperp(flavour, param, X, k_perp, zeta); } /*************************************************************************/ /*************************************************************************/ double nonforward_heldistr_kperp(int flavour, int param, double x, double k_perp, double zeta) { x_s5 = x; zeta_s5 = zeta; flavour_s5 = flavour; param_s5 = param; k_perp_s5 = k_perp; if (x >= zeta) return int_gauss5(int_double_heldistr_kperp,0.,(1. - x)/(1. - zeta),20); else return int_gauss5(int_double_heldistr_kperp, 0., x/zeta, 20); } /*************************************************************************/ /*************************************************************************/ double int_double_heldistr_kperp(double y) { return double_heldistr_kperp(flavour_s5, param_s5, x_s5 - zeta_s5 * y, y, k_perp_s5); } double double_heldistr_kperp(int flavour, int param, double x, double y, double k_perp) { double qheldistr, norm, gaussian, scale_sqr; scale_sqr = 2.; /* fixed scale */ if(x == 0) printf("\n argument 0!!!"); switch(flavour) { case 1: /* valence up quark helicity distribution */ { qheldistr = pol_up_valence(param, x, scale_sqr); break; } case 2: /* valence down quark distribution */ { qheldistr = pol_down_valence(param, x, scale_sqr); break; } } gaussian = 4. * pow(PI,2.)/( SIGMA_NU * y * (1. - x - y) ) * exp(- pow(k_perp,2.)/(2. * SIGMA_NU * y * (1. - x - y)) ); norm = pow(1. - x,3.); return 6. * y * (1. - x - y) * qheldistr/norm * gaussian; }