/**************************************************************************/ /**************************************************************************/ /* */ /* This program calculates the hard (high Q^2, low t) */ /* pion electroproduction observables on a nucleon */ /* e + N -> e + pi + N */ /* */ /* All dimensional quantities in this program are in units GeV. */ /* */ /**************************************************************************/ /**************************************************************************/ /*************************************************************/ /* */ /* Authors : Marc Vanderhaeghen */ /* Michel Guidal */ /* */ /* First run : 13/10/1997 */ /* */ /* Last update : 04/11/1999 */ /* */ /*************************************************************/ /************************************/ /* GENERAL STRUCTURE OF THE PROGRAM */ /************************************/ /**************************************************************************/ /* */ /* -> main() : all input/output of the program */ /* -> e + N -> e + pi + N : CROSS SECTIONS and ASYMMETRIES */ /* -> e + N -> e + pi + N : RESPONSE FUNCTIONS */ /* -> gamma^* + N -> pi + N AMPLITUDES */ /* */ /**************************************************************************/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include /**************************************************************************/ /**************************************************************************/ #define PI_M 0.138 /* 0.138 */ /* isospin averaged pion mass in GeV */ #define g_PINN 13.13 /* 13.13 */ //#define g_PINN 14.5 /* JML gpiNN value*/ /* piNN coupling from VPI analysis '98 */ /* g**2/(4 pi) = 13.72 */ /**************************************************************************/ /**************************************************************************/ double dsigma_dt(int reaction, int param, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t); double int_dsigma_dt(double phi); double dsigma_dtdphi(int reaction, int param, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t); double dsigma_dt_anal(int reaction, int param, double el_hel, double epsilon, double Q_sqr, double s, double t); double sigma_long_tot(int reaction, int param, double Q_sqr, double s, double tmin); double int_sigma_long_tot(double t); double dsigma_lab_fivefold(int reaction, int param, double k_in_mom_lab, double Q_sqr, double s, double mes_th_lab, double Phi); double dsigma_dqsqr_dy_dt(int reaction, int param, double sqrt_sep,double Q_sqr, double x_b, double t); double perp_spin_asymm(int reaction,int param,int wavef, double x_b, double t); double extrapol_spin_asymm(int reaction, int param, int wavef, double x_b, double t); /**************************************************************************/ double response_R_T(int reaction, int param, double Q_sqr, double s, double t); double response_R_L(int reaction, int param, double Q_sqr, double s, double t); double response_R_TT(int reaction, int param, double Q_sqr,double s,double t); double response_R_TL(int reaction, int param, double Q_sqr,double s,double t); double response_R_TLp(int reaction, int param, double Q_sqr,double s,double t); double response_R_L_ortho(int reaction, int param, double Q_sqr,double s,double t); dcomplex J_hadron(int reaction, int param, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in); /**************************************************************************/ dcomplex t_gap_ppi_hard(int reaction, int param, int hardsoft, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in); dcomplex Jlong_pilong_prod_1gluon(int reaction, int param, double sp_out, double sp_in, vec p_in, vec p_out, fvec q4_in); double intd_axial_part(double x); double intc_axial_part(double x); double int_soft_integral(double x); /**************************************************************************/ double Hhel_p_distr(int reaction, int param, double x, double xi, double t); double nonforward_p_heldistr(int reaction, int param, double x, double zeta, double t); /* ********************************************************************* */ /* ********************************************************************* */ FILE *fp, *fp1, *fp2, *fp3, *fp4, *fp5; /* DO NOT USE static _s5 variables as they are used in subroutine ofpd.c */ static int reaction_s, param_s, wavefunction_s, hardsoft_s,flavour_s,HplusE_s, deriv_s,norm_s, choice_scale_s, parpipole, reaction_s2, param_s2, flavour_s2, choice_alpha_s2, reaction_s3, param_s3, reaction_s4, param_s4; static long npoints_mc; static double el_hel_s, phi_s, epsilon_s, Q_sqr_s, s_s, xi_s, t_s, x_s, zeta_s, b_s, Q_sqr_s2, s_s2, xi_s2,t_s2,zeta_s2,z_s2,k_perp_s2,scale_sqr_s2, Q_sqr_s3, s_s3, xi_s3, t_s3, x_s3, b_s3, tmin_s3,tmax_s3, Q_sqr_s4, s_s4, xi_s4, t_s4, mperp_s4, xi_perp_s4, zeta_s4, scale_sqr_s4; main() { int reaction, mechanism, param, observable, dummy, dummy2, dummy2b, dummy3; double a, b, c, amin, bmin, cmin, Q_sqr, Q_sqr2, Q_sqr_min, Q_sqr_max, Q_sqr_step, x_b, k_in_mom_lab, ga_en_in_lab, ga_mom_in_lab, epsilon, s, t, W, tmin, tmin_test, tmin_approx, phi, mt, mt_min, mt_max, mt_step, mt2, pi_th_lab, pi_mom_lab, pi_mom_lab_test, pi_en_lab, ga_mom_cm, ga_en_cm, pi_mom_cm, pi_en_cm, result1, result2, result3, result4, test; char filed[80], fileqcd[80]="_pqcd", filehardsoft[80], filescale[80]; /* ********************************************************************* */ /* ********************************************************************* */ /* */ /* Choose the reaction you want to calculate : */ /* -1. (e, e' p pi_minus) */ /* 0. (e, e' p pi_zero) */ /* 1. (e, e' n pi_plus) */ /* */ /* ********************************************************************* */ /* ********************************************************************* */ printf("\n Give the isospin channel you want to calculate?"); printf("\n -1 = pi^- production"); printf("\n 0 = pi^0 production"); printf("\n 1 = pi^+ production"); printf("\n"); scanf("%d", &reaction); printf("\n Give the contribution you want to calculate?"); printf("\n 1 : perturbative 1-gluon exchange"); printf("\n 2 : soft overlap contribution"); printf("\n 3 : sum of 1 and 2 (still under construction)"); printf("\n"); scanf("%d", &hardsoft_s); if(hardsoft_s != 2) { printf("\nDo you want to evaluate the charged pion pole contribution?"); printf("\n1 = Axial vector contribution and NO pion pole"); printf("\n2 = ONLY pion pole contribution (using L.O. PQCD pion EM FF)"); printf("\n3 = Axial vector contribution + pion pole (using L.O. PQCD pion EM FF)"); printf("\n4 = ONLY pion pole contribution (using power corrected pion EM FF)"); printf("\n5 = Axial vector contribution + pion pole (using power corrected pion EM FF)"); printf("\n6 = ONLY pion pole contribution (using phenomenological pion EM FF)"); printf("\n7 = Axial vector contribution + pion pole (using phenomenological pion EM FF)"); printf("\n"); scanf("%d",&parpipole); } if((parpipole != 2) && (parpipole != 4) && (parpipole != 6)) { printf("\n Choose model for polarized OFPD"); printf("\n 1 : factorized model based on MRS S_o distribution"); printf("\n 6 : factorized model based on MRST98 distribution at Q^2 = 1 GeV^2 (Leader et al.)"); printf("\n 21 : xi dependence based on MRS S_o distribution "); printf("\n 26 : xi dependence based on MRST98 distribution at Q^2 = 1 GeV^2 (Leader et al.)"); printf("\n"); scanf("%d", ¶m); } if((parpipole != 4) && (parpipole != 6)) { printf("\n Choice for the strong coupling alpha_s :"); printf("\n 1 = constant coupling 0.56 at scale around 1 GeV^2"); printf("\n 2 = running coupling (LAMBDA_QCD = 0.2 GeV)"); printf("\n 3 = running coupling which freezes to a maximum value of 0.5"); printf("\n 4 = analytical model for the strong coupling which is IR finite"); printf("\n"); scanf("%d", &choice_alpha_s2); } switch(hardsoft_s) { case 1: { char filetemp[80]=".dat"; char filetemp2[80]="_pipole_asymff.dat"; char filetemp3[80]="_pipole_ff.dat"; char filetemp4[80]="_pipole_phenomff.dat"; if(parpipole == 1) strcpy(filehardsoft,filetemp); if((parpipole == 2) || (parpipole == 3)) strcpy(filehardsoft,filetemp2); if((parpipole == 4) || (parpipole == 5)) strcpy(filehardsoft,filetemp3); if((parpipole == 6) || (parpipole == 7)) strcpy(filehardsoft,filetemp4); break; } case 2: { char filetemp[80]="_soft.dat"; strcpy(filehardsoft,filetemp); break; } case 3: { char filetemp[80]="_sum.dat"; char filetemp2[80]="_sum_pipole_asymff.dat"; char filetemp3[80]="_sum_pipole_ff.dat"; char filetemp4[80]="_sum_pipole_phenomff.dat"; if(parpipole == 1) strcpy(filehardsoft,filetemp); if((parpipole == 2) || (parpipole == 3)) strcpy(filehardsoft,filetemp2); if((parpipole == 4) || (parpipole == 5)) strcpy(filehardsoft,filetemp3); if((parpipole == 6) || (parpipole == 7)) strcpy(filehardsoft,filetemp4); break; } } if((parpipole != 2) && (parpipole != 4) && (parpipole != 6)) { switch(param) { case 1 : { char fileout[80]="_pol_mrss0"; strcpy(filed,strcat(fileqcd, strcat(fileout,filehardsoft)) ); break; } case 6 : { char fileout[80]="_pol_mrst98"; strcpy(filed,strcat(fileqcd, strcat(fileout,filehardsoft)) ); break; } case 21 : { char fileout[80]="_xidep_pol_mrss0"; strcpy(filed,strcat(fileqcd, strcat(fileout,filehardsoft)) ); break; } case 26 : { char fileout[80]="_xidep_pol_mrst98"; strcpy(filed,strcat(fileqcd, strcat(fileout,filehardsoft)) ); break; } } } else { strcpy(filed,strcat(fileqcd,filehardsoft) ); } printf("\n Choose one of the following menu items :"); printf("\n 1 : DIFFERENTIAL CROSS SECTION as function of t"); printf("\n 2 : DIFFERENTIAL CROSS SECTION as function of W"); printf("\n 3 : DIFFERENTIAL CROSS SECTION as function of x_B"); printf("\n 4 : DIFFERENTIAL CROSS SECTION as function of Q^2"); printf("\n 5 : TOTAL CROSS SECTION"); printf("\n 6 : 5-FOLD DIFFERENTIAL CROSS SECTION in LAB"); printf("\n 7 : 3-FOLD DIFFERENTIAL CROSS SECTION (Q^2, y, t) for collider kinematics"); printf("\n : plot as function of t"); printf("\n 9 : SPIN ASYMMETRY as function of x_B at fixed t"); printf("\n 10 : SPIN ASYMMETRY as function of t at fixed x_B"); printf("\n 11 : extrapolation of pion FF from SPIN ASYMMETRY at fixed x_B"); printf("\n"); scanf("%d", &observable); switch(observable) { case 1 : /* differential cross section */ { char fileb[80], filexb[80], filel1[80], filel[80]; char filec[80]="_x0p", filez[80]; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_qsqr"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_qsqr"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_qsqr"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); dummy = (int) Q_sqr; if (Q_sqr>9.) sprintf(fileb,"%2i",dummy); else sprintf(fileb,"%1i",dummy); printf("\n\n Give the value of x_B (e.g. 0.3)"); printf("\n"); scanf("%lf", &x_b); dummy2 = (int) (10. * x_b); sprintf(filexb,"%1i",dummy2); strcpy(filel1,strcat(filel,fileb)); strcpy(filez,strcat(filel1,strcat(strcat(filec,filexb),filed))); ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); pi_mom_cm = sqrt( ( pow(s, 2.) - 2. * s * (pow(NU_M,2.) + pow(PI_M,2.)) + pow(pow(NU_M,2.) - pow(PI_M,2.),2.) ) / (4. * s)); pi_en_cm = sqrt( pow(pi_mom_cm, 2.) + pow(PI_M, 2.) ); tmin = - Q_sqr + pow(PI_M,2.) - 2. * (ga_en_cm * pi_en_cm - ga_mom_cm * pi_mom_cm); printf("\n\n Give the lowest value of -t (in GeV^2, e.g. tmin)\n"); scanf("%lf", &mt_min); if(mt_min==0) mt_min=-tmin; printf("\n\n Give the highest value of -t (in GeV^2, e.g. 1.5)\n"); scanf("%lf", &mt_max); if(mt_max==0) mt_max=1.5; printf("\n Give the step in t (in GeV^2, e.g. 0.1)\n"); scanf("%lf", &mt_step); for(mt = mt_min; mt <= mt_max; mt += mt_step) { if(mt<-tmin){ printf("The value chosen for t %lf is below the tmin calculated: %lf",mt,-tmin); exit(0); } if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); t = -mt; result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%e \t %e \t %e\n", mt, -tmin, result1); fclose(fp1); } /* for(pi_th_lab = 0.; pi_th_lab <= 5.; pi_th_lab += .25) { a = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab * cos(pi_th_lab * PI/180.), 2.); b = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab * cos(pi_th_lab * PI/180.); c = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab * cos(pi_th_lab * PI/180.)); if ((b*b + a*c) > 0) { pi_mom_lab = 1./a * ( b + sqrt(pow(b,2.) + a * c)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %e \t %e\n", pi_th_lab, -t, result1); } fclose(fp1); } */ break; } case 2 : /* differential cross section as function of W */ { double it; char fileq[80]="_qsqr"; char filea[80], fileb[80], filel[80]; char filec[80]="_pqcd.dat"; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_t"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_t"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_t"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n\n Give the value of -t in GeV^2 (e.g. 1.0, 0.0 => tmin)"); printf("\n"); scanf("%lf", &it); t=it; dummy = (int) Q_sqr; if (Q_sqr>9.) sprintf(fileb,"%2i",dummy); else sprintf(fileb,"%1i",dummy); dummy2 = (int) t; if (t>9.) sprintf(filea,"%2i",dummy2); else sprintf(filea,"%1i",dummy2); if ( (fp1 = fopen(strcat(strcat(strcat(filel,filea), strcat(fileq,fileb)),filed),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(W = 2.; W <= 25.; W += .5) { s = pow(W,2.); ga_en_in_lab = (s - pow(NU_M,2.) + Q_sqr) / (2. * NU_M); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); amin = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab,2.); bmin = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab; cmin = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin_test = - 2. * NU_M * (ga_en_in_lab - pi_mom_lab_test); if ((bmin*bmin + amin*cmin) > 0) { pi_mom_lab = 1./amin * ( bmin + sqrt(pow(bmin,2.) + amin * cmin)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); } result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %e \t %e\n", W, -t, result1); } fclose(fp1); break; } case 3 : /* differential cross section as function of x_B */ { double mt; char fileq[80]="_qsqr"; char filet1[80], filet2[80], fileb[80], filel[80]; char filec[80]="_pqcd.dat"; char filez[80]; char filetp[80]="p"; char file00[80]="00"; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_xb_t"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_xb_t"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_xb_t"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n\n Give the value of -t in GeV^2 (e.g. 1.0, 0.0 => tmin)"); printf("\n"); scanf("%lf", &mt); dummy = (int) Q_sqr; if (Q_sqr>9.) sprintf(fileb,"%2i",dummy); else sprintf(fileb,"%1i",dummy); dummy2 = (int) (mt); sprintf(filet1,"%1i",dummy2); mt2 = 100. * ( mt - dummy2); if(mt2 > 0.) { dummy2b = (int) (mt2); sprintf(filet2,"%2i",dummy2b); } else strcpy(filet2,file00); strcpy(filez, strcat(strcat(strcat(filel,strcat(filet1,strcat(filetp,filet2))),strcat(fileq,fileb)),filed) ); for(x_b = 0.01; x_b <= .95; x_b += .01) { if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); s = pow(NU_M,2.) + Q_sqr * (1 - x_b)/x_b; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); pi_mom_cm = sqrt( ( pow(s, 2.) - 2. * s * (pow(NU_M,2.) + pow(PI_M,2.)) + pow(pow(NU_M,2.) - pow(PI_M,2.),2.) ) / (4. * s)); pi_en_cm = sqrt( pow(pi_mom_cm, 2.) + pow(PI_M, 2.) ); tmin = - Q_sqr + pow(PI_M,2.) - 2. * (ga_en_cm * pi_en_cm - ga_mom_cm * pi_mom_cm); tmin_test = pow((- Q_sqr - pow(PI_M,2.) - pow(NU_M,2.) + pow(NU_M,2.))/(2. * sqrt(s)),2.) - pow(ga_mom_cm - pi_mom_cm, 2.); tmin_approx = - pow(NU_M,2.) * pow(x_b,2.)/(1. - x_b); /* printf("\nx_b = %f \t tmin = %e \t tmin_approx = %e", x_b, tmin, tmin_approx); */ if( mt == 0.) { t = tmin; result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %e \t %e\n", x_b, -t, result1); } else { t = -mt; if( t < tmin) { result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %e \t %e\n", x_b, -t, result1); } } fclose(fp1); } break; } case 4 : /* differential cross section as function of Q^2 */ { int xwoption; int choice_flux; char fileb[80], filez[80]; char filetmin[80]="_tmin"; char filet[80]="_mt"; char filetph[80]="_nophase_mt"; char filetp[80]="p"; char filet1[80], filet2[80]; char filel[80]; printf("\n\n fixed x_B or W ?"); printf("\n 1 = dsigma/dt at fixed x_B"); printf("\n 2 = dsigma/domega at fixed W"); printf("\n"); scanf("%d", &xwoption); switch(reaction) { case -1 : { char fileout1[]="en_pimp_long_xb0p"; char fileout2[]="en_pimp_long_w"; if(xwoption == 1) strcpy(filel,fileout1); if(xwoption == 2) strcpy(filel,fileout2); break; } case 0 : { char fileout1[]="ep_piop_long_xb0p"; char fileout2[]="ep_piop_long_w"; if(xwoption == 1) strcpy(filel,fileout1); if(xwoption == 2) strcpy(filel,fileout2); break; } case 1 : { char fileout1[]="ep_pipn_long_xb0p"; char fileout2[]="ep_pipn_long_w"; if(xwoption == 1) strcpy(filel,fileout1); if(xwoption == 2) strcpy(filel,fileout2); break; } } switch(xwoption) { case 1 : /* fixed x_B */ { printf("\n\n Give the value of x_B (e.g. 0.3)"); printf("\n"); scanf("%lf", &x_b); dummy = (int) (10. * x_b); sprintf(fileb,"%1i",dummy); printf("\n Give the value of -t (if 0. then t = t_min)"); printf("\n"); scanf("%lf", &mt); if(mt == 0.) { strcpy(filez, strcat(strcat(filel,fileb), strcat(filetmin,filed) ) ); printf("\n Give the lowest value of Q_sqr (in GeV^2, e.g. 2.)\n"); scanf("%lf", &Q_sqr_min); printf("\n Give the highest value of Q_sqr (in GeV^2, e.g. 25.)\n"); scanf("%lf", &Q_sqr_max); printf("\n Give the step in Q_sqr (in GeV^2, e.g. 1.)\n"); scanf("%lf", &Q_sqr_step); for(Q_sqr = Q_sqr_min; Q_sqr <= Q_sqr_max; Q_sqr += Q_sqr_step) { if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; amin = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab, 2.); bmin = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab; cmin = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin_test = - 2.* NU_M * (ga_en_in_lab - pi_mom_lab_test); if ((bmin*bmin + amin*cmin) > 0) { pi_mom_lab = 1./amin * ( bmin + sqrt(pow(bmin,2.) + amin * cmin)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); } result1 = response_R_L(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %e \t %e \t %e\n", Q_sqr, -t, result1, pow(Q_sqr,3.) * result1); fclose(fp1); } } else { dummy2 = (int) (mt); sprintf(filet1,"%1i",dummy2); mt2 = 100. * ( mt - dummy2); dummy2b = (int) (mt2); sprintf(filet2,"%2i",dummy2b); if ( (fp1 = fopen(strcat(strcat(filel,fileb), strcat(strcat(filet,strcat(strcat(filet1,filetp),filet2)),filed) ),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); /*cynthia 27.03.2004: add Q2 range option for t!=tmin*/ printf("\n Give the lowest value of Q_sqr (in GeV^2, e.g. 2.)\n"); scanf("%lf", &Q_sqr_min); printf("\n Give the highest value of Q_sqr (in GeV^2, e.g. 25.)\n"); scanf("%lf", &Q_sqr_max); printf("\n Give the step in Q_sqr (in GeV^2, e.g. 1.)\n"); scanf("%lf", &Q_sqr_step); for(Q_sqr = Q_sqr_min; Q_sqr <= Q_sqr_max; Q_sqr += Q_sqr_step) { ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; result1 = response_R_L(reaction, param, Q_sqr, s, -mt); printf("%f \t %e \t %e \t %e\n", Q_sqr, mt, result1, pow(Q_sqr,3.) * result1); fprintf(fp1,"%f \t %e \t %e \t %e\n", Q_sqr, mt, result1, pow(Q_sqr,3.) * result1); } fclose(fp1); } break; } case 2 : /* fixed W */ { double w, conv; printf("\n\n Give the value of W in GeV (e.g. 2.65)"); printf("\n"); scanf("%lf", &w); dummy = (int) (100. * w); sprintf(fileb,"%3i",dummy); if ( (fp1 = fopen(strcat(filel, strcat(fileb,filed)),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); // for(Q_sqr = 2.; Q_sqr <= 25.; Q_sqr += .5) //TH - add extendend range in Q2 here for(Q_sqr = 0.35; Q_sqr <= 3.8; Q_sqr += .1) //for(Q_sqr = 0.35; Q_sqr <= 0.35; Q_sqr += .1) { s = pow(w,2.); x_b = Q_sqr/(Q_sqr+s-pow(NU_M,2.)); ga_en_in_lab=Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); amin = pow(NU_M +ga_en_in_lab,2.) - pow(ga_mom_in_lab,2.); bmin = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab; cmin = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin_test = - 2.*NU_M * (ga_en_in_lab - pi_mom_lab_test); if ((bmin*bmin + amin*cmin) > 0) { pi_mom_lab = 1./amin * ( bmin + sqrt(pow(bmin,2.) + amin * cmin)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); printf("\n\n TH - print tmin"); printf("%f \t %f \n",tmin_test,t); } /* ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); amin = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab, 2.); bmin = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab; cmin = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin_test = - 2. * NU_M * (ga_en_in_lab - pi_mom_lab_test); if ((bmin*bmin + amin*cmin) > 0) { pi_mom_lab = 1./amin * ( bmin + sqrt(pow(bmin,2.) + amin * cmin)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); }*/ ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); pi_mom_cm = sqrt(( pow(s - pow(NU_M,2.) - pow(PI_M,2.), 2.) - pow(2. * NU_M * PI_M, 2.) ) / (4. * s) ); result1 = response_R_L(reaction, param, Q_sqr, s, t); conv = 2. * pi_mom_cm * ga_mom_cm / (2. * PI); fprintf(fp1,"%f \t %e \t %e\n", Q_sqr, -t, result1 * conv); } fclose(fp1); break; } } break; } case 5 : /* total cross section as function of W */ /*cyn 26/04/2004: add option */ { char fileb[80], filel[80]; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_tot_qsqr"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_tot_qsqr"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_tot_qsqr"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); dummy = (int) Q_sqr; if (Q_sqr>9.) sprintf(fileb,"%2i",dummy); else sprintf(fileb,"%1i",dummy); strcat(filel,fileb); printf("\n\n Give the value of x_B (e.g. 0.3)"); printf("\n"); scanf("%lf", &x_b); dummy = (int) (10. * x_b); sprintf(fileb,"xb%1i",dummy); strcat(filel,fileb); strcat(filel,filehardsoft); if ( (fp1 = fopen(filel,"w")) == NULL) fprintf(stderr,"\nUnable to open file"); ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; amin = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab, 2.); bmin = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab; cmin = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin_test = - 2.* NU_M * (ga_en_in_lab - pi_mom_lab_test); if ((bmin*bmin + amin*cmin) > 0) { pi_mom_lab = 1./amin * ( bmin + sqrt(pow(bmin,2.) + amin * cmin)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); } result1 = sigma_long_tot(reaction, param, Q_sqr, s, t); fprintf(fp1,"%f \t %f \t %e \t %e\n", Q_sqr, x_b, t, result1); printf("%f \t %f \n",tmin_test,t); /*for(W = 2.; W <= 25.; W += .5) { s = pow(W,2.); ga_en_in_lab = (s - pow(NU_M,2.) + Q_sqr) / (2. * NU_M); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); pi_mom_lab = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab); tmin = - 2. * NU_M * (ga_en_in_lab - pi_mom_lab); result1 = sigma_long_tot(reaction, param, Q_sqr, s, tmin); fprintf(fp1,"%f \t %e \t %e\n", W, tmin, result1); } */ fclose(fp1); break; } case 6 : /* fivefold differential cross section in LAB */ { char fileq[80]="gev_qsqr"; char filep[80]="p"; char filez[80]="0"; char filex[80]="_x0p"; char filea[80], fileb[80], filebb[80], filec[80]; char filel[80]; double eprime,y; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_5fold_"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_5fold_"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_5fold_"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of beam energy in GeV (e.g. 200.)"); printf("\n"); scanf("%lf", &k_in_mom_lab); printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n\n Give the value of x_B (e.g. 0.3)"); printf("\n"); scanf("%lf", &x_b); dummy = (int) k_in_mom_lab; if (k_in_mom_lab>9.) sprintf(filea,"%2i",dummy); else sprintf(filea,"%1i",dummy); dummy2 = (int) Q_sqr; if (Q_sqr>10.) { sprintf(fileb,"%2i",dummy2); strcpy(filebb,filez); } else { dummy2 = (int) (Q_sqr); sprintf(fileb,"%1i",dummy2); Q_sqr2 = 10. * ( Q_sqr - dummy2); dummy2b = (int) (Q_sqr2); sprintf(filebb,"%1i",dummy2b); } dummy3 = (int) (x_b * 10.); sprintf(filec,"%1i",dummy3); ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; eprime=k_in_mom_lab-ga_en_in_lab; y=ga_en_in_lab/k_in_mom_lab; if ( (fp1 = fopen(strcat(strcat(strcat(filel,filea), strcat(fileq,strcat(fileb,strcat(filep,filebb)))), strcat(strcat(filex,filec), filed) ),"w")) == NULL) fprintf(stderr,"\nUnable to open file"); for(pi_th_lab = -25.; pi_th_lab <= 25.; pi_th_lab += 2.) { a = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab * cos(pi_th_lab * PI/180.), 2.); b = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab * cos(pi_th_lab * PI/180.); c = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); pi_mom_lab_test = .5 * (s - pow(NU_M,2.)) / (ga_en_in_lab + NU_M - ga_mom_in_lab * cos(pi_th_lab * PI/180.)); if ((b*b + a*c) > 0) { pi_mom_lab = 1./a * ( b + sqrt(pow(b,2.) + a * c)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); result1 = dsigma_lab_fivefold(reaction, param, k_in_mom_lab, Q_sqr, s, pi_th_lab * PI/180., 0.); result2=result1*2*PI*NU_M*ga_en_in_lab/eprime/x_b/(s-pow(NU_M,2.)); result3=result1*2*PI/2./eprime/k_in_mom_lab*ga_en_in_lab/x_b; result4=result3*Q_sqr/y; fprintf(fp1,"%f \t %e \t %e \t %e\n", pi_th_lab, -t, result1, result3); printf("%f \t %e \t %e \t %e \t %e \t %e\n", pi_th_lab, -t, result1, result2, result3, result4); } } fclose(fp1); break; } case 7 : /* 3-fold diff. cross section for collider kinematics */ { double sqrt_sep, y; char fileq[80]="gev_qsqr"; char filep[80]="p"; char filez[80]="0"; char filex[80]="_x0p"; char filea[80], fileb[80], filebb[80], filec[80]; char filel[80]; char filetot[80]; switch(reaction) { case -1 : { char fileout[]="en_pimp_long_3fold_"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_long_3fold_"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_long_3fold_"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of the collider energy sqrt(s_ep) in GeV (e.g. 40.)"); printf("\n"); scanf("%lf", &sqrt_sep); printf("\n\n Give the value of Q^2 in GeV^2 (e.g. 5.0)"); printf("\n"); scanf("%lf", &Q_sqr); printf("\n\n Give the value of x_B (e.g. 0.2)"); printf("\n"); scanf("%lf", &x_b); y = Q_sqr / ( x_b * (pow(sqrt_sep,2.) - pow(NU_M,2.)) ); k_in_mom_lab = (pow(sqrt_sep,2.) - pow(NU_M,2.))/(2.* NU_M); epsilon = ( 1. - y - Q_sqr/(4. * pow(k_in_mom_lab,2.)) ) / ( 1. - y + .5 * pow(y,2.) + Q_sqr/(4. * pow(k_in_mom_lab,2.)) ); printf("\n E_e_lab = %e \t y = %e \t eps = %e \n", k_in_mom_lab, y, epsilon); dummy = (int) sqrt_sep; if (sqrt_sep>10.) sprintf(filea,"%2i",dummy); else sprintf(filea,"%1i",dummy); dummy2 = (int) Q_sqr; if (Q_sqr>10.) { sprintf(fileb,"%2i",dummy2); strcpy(filebb,filez); } else { dummy2 = (int) (Q_sqr); sprintf(fileb,"%1i",dummy2); Q_sqr2 = 10. * ( Q_sqr - dummy2); dummy2b = (int) (Q_sqr2); sprintf(filebb,"%1i",dummy2b); } dummy3 = (int) (x_b * 10.); sprintf(filec,"%1i",dummy3); ga_en_in_lab = Q_sqr / (2. * NU_M * x_b); ga_mom_in_lab = sqrt(Q_sqr + pow(ga_en_in_lab,2.)); s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; strcpy(filetot,strcat(strcat(strcat(filel,filea),strcat(fileq,strcat(fileb,strcat(filep,filebb)))),strcat(strcat(filex,filec),filed) )); for(pi_th_lab = 0.; pi_th_lab <= 5.; pi_th_lab += .25) { if ( (fp1 = fopen(filetot,"a")) == NULL) fprintf(stderr,"\n Unable to open file"); a = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab * cos(pi_th_lab * PI/180.), 2.); b = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab * cos(pi_th_lab * PI/180.); c = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); if ((b*b + a*c) > 0) { pi_mom_lab = 1./a * ( b + sqrt(pow(b,2.) + a * c)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); result1 = dsigma_dqsqr_dy_dt(reaction, param, sqrt_sep, Q_sqr, x_b, t); fprintf(fp1,"%f \t %e \t %e\n", pi_th_lab, -t,result1); printf("\n%f \t %e \t %e", pi_th_lab, -t, result1); fclose(fp1); } } break; } case 9 : /* spin asymmetry as function of x_B at fixed t */ { int wavef; double mt; char filea[80],filea2[80],fileb[80],filel[80], filewf[80]; char filec[80]="_pqcd.dat"; char filez[80]; char fileop[80]="0p"; printf("\n\n Give the choice for the pion distribution amplitude"); printf("\n 1 = asymptotic DA"); printf("\n 2 = CZ pion DA"); printf("\n"); scanf("%d", &wavef); switch(wavef) { case 1 : { char fileout[]="_asymptotic"; strcpy(filewf,fileout); break; } case 2 : { char fileout[]="_cz"; strcpy(filewf,fileout); break; } } switch(reaction) { case -1 : { char fileout[]="en_pimp_spin_asymm_t"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_spin_asymm_t"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_spin_asymm_t"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of -t in GeV^2"); printf("\n"); scanf("%lf", &mt); t = -mt; dummy2 = (int) mt; if (mt>9.) sprintf(filea,"%2i",dummy2); else { dummy3 = (int) (10. * mt); sprintf(filea2,"%1i",dummy3); strcpy(filea,strcat(fileop,filea2)); } strcpy(filez, strcat(strcat(filel,filea),strcat(filewf,filed)) ); for(x_b = 0.05; x_b <= .5; x_b += .01) { if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); result1 = perp_spin_asymm(reaction, param, wavef, x_b, t); Q_sqr=10.; s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; result2 = response_R_L_ortho(reaction, param, Q_sqr, s, t); result3=response_R_L(reaction, param, Q_sqr, s, t); printf("%f \t %e \t %e \t %e \t %e \t %e\n", x_b, result2, result3, result1,result2/result3,result2/result3*2./PI); fprintf(fp1,"%f \t %e \t %e\n", x_b, t, result1); fclose(fp1); } break; } case 10 : /* spin asymmetry at fixed t as function of x_B */ { int wavef; char filea[80],filea2[80],filexb[80],filel[80], filewf[80]; char filec[80]="_pqcd.dat"; char filez[80]; printf("\n\n Give the choice for the pion distribution amplitude"); printf("\n 1 = asymptotic DA"); printf("\n 2 = CZ pion DA"); printf("\n"); scanf("%d", &wavef); switch(wavef) { case 1 : { char fileout[]="_asymptotic"; strcpy(filewf,fileout); break; } case 2 : { char fileout[]="_cz"; strcpy(filewf,fileout); break; } } switch(reaction) { case -1 : { char fileout[]="en_pimp_spin_asymm_xb0p"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_spin_asymm_xb0p"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_spin_asymm_xb0p"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of x_B (e.g. 0.3)"); printf("\n"); scanf("%lf", &x_b); dummy2 = (int) (10. * x_b); sprintf(filexb,"%1i",dummy2); /* Q_sqr = 10.; s = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); pi_mom_cm = sqrt( ( pow(s, 2.) - 2. * s * (pow(NU_M,2.) + pow(PI_M,2.)) + pow(pow(NU_M,2.) - pow(PI_M,2.),2.) ) / (4. * s)); pi_en_cm = sqrt( pow(pi_mom_cm, 2.) + pow(PI_M, 2.) ); tmin = - Q_sqr + pow(PI_M,2.) - 2. * (ga_en_cm * pi_en_cm - ga_mom_cm * pi_mom_cm); */ tmin_approx = - pow(NU_M,2.) * pow(x_b,2.)/(1. - x_b); tmin = tmin_approx; strcpy(filez, strcat(strcat(filel,filexb),strcat(filewf,filed)) ); for(mt = -tmin; mt <= 1.; mt += .01) { if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); t = -mt; result1 = perp_spin_asymm(reaction, param, wavef, x_b, t); fprintf(fp1,"%f \t %e \t %e\n", x_b, mt, result1); fclose(fp1); } break; } case 11 : /* extrapolation of pion FF from spin asymmetry at fixed x_B */ { int wavef; double mt; char filea[80],filea2[80],fileb[80],filel[80], filewf[80]; char filec[80]="_pqcd.dat"; char filez[80]; char fileop[80]="0p"; char fileopo[80]="0p0"; printf("\n\n Give the choice for the pion distribution amplitude"); printf("\n 1 = asymptotic DA"); printf("\n 2 = CZ pion DA"); printf("\n"); scanf("%d", &wavef); switch(wavef) { case 1 : { char fileout[]="_asymptotic"; strcpy(filewf,fileout); break; } case 2 : { char fileout[]="_cz"; strcpy(filewf,fileout); break; } } switch(reaction) { case -1 : { char fileout[]="en_pimp_spin_asymm_xb"; strcpy(filel,fileout); break; } case 0 : { char fileout[]="ep_piop_spin_asymm_xb"; strcpy(filel,fileout); break; } case 1 : { char fileout[]="ep_pipn_spin_asymm_xb"; strcpy(filel,fileout); break; } } printf("\n\n Give the value of x_b"); printf("\n"); scanf("%lf", &x_b); dummy2 = (int) x_b; dummy3 = (int) (100. * x_b); if(dummy3 < 10) { sprintf(filea2,"%1i",dummy3); strcpy(filea,strcat(fileopo,filea2)); } else { sprintf(filea2,"%2i",dummy3); strcpy(filea,strcat(fileop,filea2)); } strcpy(filez, strcat(strcat(filel,filea),strcat(filewf,filed)) ); for(mt = 0.0025; mt <= .3; mt += .0025) { if ( (fp1 = fopen(filez,"a")) == NULL) fprintf(stderr,"\nUnable to open file"); result1 = extrapol_spin_asymm(reaction,param,wavef,x_b,-mt); fprintf(fp1,"%f \t %e \n", mt, result1); fclose(fp1); } break; } } } /*************************************************************************/ /*************************************************************************/ /* */ /* e + N -> e + pi + N : CROSS SECTIONS and ASYMMETRIES */ /* */ /*************************************************************************/ /*************************************************************************/ /* The function dsigma_dt() gives the e + N -> e + pi + N differential cross section dsigma/dt integrated over Phi as function of electron helicity : el_hel, epsilon, Q_sqr, s and t. The units of the cross section are (nanobarns / GeV^2). */ double dsigma_dt(int reaction, int param, double el_hel, double phi_min, double phi_max, double epsilon, double Q_sqr, double s, double t) { reaction_s = reaction; param_s2 = param; el_hel_s = el_hel; epsilon_s = epsilon; Q_sqr_s = Q_sqr; s_s = s; t_s = t; return int1d(int_dsigma_dt, phi_min, phi_max, 4); } double int_dsigma_dt(double phi) { return dsigma_dtdphi(reaction_s, param_s2, el_hel_s, phi, epsilon_s, Q_sqr_s, s_s, t_s); } /* The function dsigma_dtdphi() gives the e + N -> e + pi + N differential cross section dsigma/dt dPhi as function of electron helicity : el_hel, out-of-plane angle Phi, epsilon, Q_sqr, s and t. The units of the cross section are (nanobarns / GeV^2). */ double dsigma_dtdphi(int reaction, int param, double el_hel, double phi, double epsilon, double Q_sqr, double s, double t) { double unit_conv; /* unit_conv = 1. / (2. * PI) * 1. / (16. * PI * pow(s - pow(NU_M,2.),2.)) * pow(0.197, 2.) * pow(10., 7.); */ unit_conv = 1. / (2. * PI); return unit_conv * ( response_R_T(reaction, param, Q_sqr, s, t) + epsilon * response_R_L(reaction, param, Q_sqr, s, t) - epsilon * cos(2. * phi) * response_R_TT(reaction, param, Q_sqr, s, t) + sqrt(2. * epsilon * (1. + epsilon)) * cos(phi) * response_R_TL(reaction, param, Q_sqr, s, t) + el_hel * sqrt(2. * epsilon * (1. - epsilon)) * sin(phi) * response_R_TLp(reaction, param, Q_sqr, s, t) ); } /* The function dsigma_dt_anal() gives the e + N -> e + pi + N differential cross section dsigma/dt as function of electron helicity : el_hel, epsilon, Q_sqr, s and t. The units of the cross section are (nanobarns / GeV^2). */ double dsigma_dt_anal(int reaction, int param, double el_hel, double epsilon, double Q_sqr, double s, double t) { return response_R_T(reaction, param, Q_sqr, s, t) + epsilon * response_R_L(reaction, param, Q_sqr, s, t); } double sigma_long_tot(int reaction, int param, double Q_sqr, double s, double tmin) { double result; reaction_s2 = reaction; param_s2 = param; Q_sqr_s2 = Q_sqr; s_s2 = s; result = int_gauss(int_sigma_long_tot, tmin - 4., tmin, 100); return result; } double int_sigma_long_tot(double t) { return response_R_L(reaction_s2, param_s2, Q_sqr_s2, s_s2, t); } /*************************************************************************/ /* The function dsigma_lab_fivefold() gives the e + N -> e + pi + N fivevold differential cross section. The units of the cross section are (nanobarns / (GeV sr^2). */ double dsigma_lab_fivefold(int reaction, int param, double k_in_mom_lab, double Q_sqr, double s, double mes_th_lab, double Phi) { int variable; double virtual_photon_flux, jacobian, x_b, k_out_mom_lab, k_out_th_lab, ga_en_in_lab, ga_mom_in_lab, epsilon, ga_mom_in_cm, mes_mom_out_cm, a, b, c, pi_mom_lab, pi_en_lab, t; x_b = Q_sqr / (s + Q_sqr - pow(NU_M,2.)); k_out_mom_lab = k_in_mom_lab - Q_sqr / (2. * NU_M * x_b); k_out_th_lab = 2. * asin(sqrt(Q_sqr / (4. * k_in_mom_lab * k_out_mom_lab))); ga_en_in_lab = k_in_mom_lab - k_out_mom_lab; ga_mom_in_lab = sqrt(pow(k_in_mom_lab, 2.) + pow(k_out_mom_lab, 2.) - 2. * k_in_mom_lab * k_out_mom_lab * cos(k_out_th_lab) ); epsilon = 1. / (1. + 2. * pow(ga_mom_in_lab, 2.) / Q_sqr * pow(tan(0.5 * k_out_th_lab), 2.) ); virtual_photon_flux = pow(ELEC,2.)/(4. * PI) * 1./(2. * pow(PI,2.)) * k_out_mom_lab / k_in_mom_lab * (s - pow(NU_M,2.))/(2. * NU_M) * 1./Q_sqr * 1./(1. - epsilon); ga_mom_in_cm = NU_M/sqrt(s) * ga_mom_in_lab; mes_mom_out_cm = (s - pow(NU_M,2.))/(2. * sqrt(s)); a = pow(NU_M + ga_en_in_lab, 2.) - pow(ga_mom_in_lab * cos(mes_th_lab), 2.); b = 0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)) * ga_mom_in_lab * cos(mes_th_lab); c = pow(0.5 * (s - pow(NU_M,2.) + pow(PI_M, 2.)), 2.) - pow(PI_M * (NU_M + ga_en_in_lab), 2.); if ((b*b + a*c) > 0) { pi_mom_lab = 1./a * ( b + sqrt(pow(b,2.) + a * c)); pi_en_lab = sqrt(pow(pi_mom_lab,2.) + pow(PI_M,2.)); t = - 2. * NU_M * (ga_en_in_lab - pi_en_lab); jacobian = s / pow(ga_en_in_lab + NU_M - ga_mom_in_lab * cos(mes_th_lab),2.); variable = 2; switch(variable) { case 1: /* dsigma/dk_e domega_e domega^M_Lab */ { return virtual_photon_flux * 2. * ga_mom_in_cm * mes_mom_out_cm * jacobian * 1./(2. * PI) * epsilon * response_R_L(reaction, param, Q_sqr, s, t); break; } case 2: /* dsigma/dk_e domega_e dt */ { return virtual_photon_flux * epsilon * response_R_L(reaction, param, Q_sqr, s, t); break; } } } else return 0.; } /* The function dsigma_dqsqr_dy_dt() gives the e + N -> e + pi + N 3-fold differential cross section. The units of the cross section are ( nanobarns/GeV^4 ). */ double dsigma_dqsqr_dy_dt(int reaction, int param, double sqrt_sep, double Q_sqr, double x_b, double t) { int variable; double virtual_photon_flux, W_sqr, y, k_in_mom_lab, epsilon; W_sqr = pow(NU_M,2.) + Q_sqr * (1. - x_b)/x_b; y = Q_sqr / ( x_b * (pow(sqrt_sep,2.) - pow(NU_M,2.)) ); k_in_mom_lab = (pow(sqrt_sep,2.) - pow(NU_M,2.)) / (2. * NU_M); epsilon = ( 1. - y - Q_sqr/(4. * pow(k_in_mom_lab,2.)) ) / ( 1. - y + .5 * pow(y,2.) + Q_sqr/(4. * pow(k_in_mom_lab,2.)) ); virtual_photon_flux = pow(ELEC,2.)/(4. * PI) * 1./(2. * PI) * 1./Q_sqr * (1. - x_b) * y / (1. - epsilon); return virtual_photon_flux * epsilon * response_R_L(reaction, param, Q_sqr, W_sqr, t); } /*************************************************************************/ /*************************************************************************/ /* the function perp_spin_asymm() returns the spin asymmetry for a proton target polarized perpendicular to the reaction plane. Definition by Frankfurt, Pobylitsa, Polyakov and Strikman, PRD 69, 014010 (1999). */ double perp_spin_asymm(int reaction,int param,int wavef, double x_b, double t) { dcomplex direct, crossed; double delta_perp_sqr, delta_perp, ff_ps, xi, re_a, im_a, sqr_a, piwavef; delta_perp_sqr = -t * 4. * (1. - x_b) / pow(2. - x_b,2.) - pow(NU_M,2.) * 4. * pow(x_b,2.) / pow(2. - x_b,2.); if(delta_perp_sqr > 0.) { delta_perp = sqrt(delta_perp_sqr); ff_ps = form_factor_G_P(t); xi = x_b / (2. - x_b); reaction_s4 = reaction; param_s4 = param; xi_s4 = xi; t_s4 = t; direct = Complex(int_gauss4(intd_axial_part, 0., 1., 30) + Hhel_p_distr(reaction,param,xi,xi,t) * log((1. - xi)/xi), - PI * Hhel_p_distr(reaction,param,xi, xi, t) ); crossed = Complex(int_gauss4(intc_axial_part, 0., 1., 30), 0.0); re_a = 6. * (direct.re + crossed.re); im_a = 6. * direct.im; sqr_a = pow(re_a, 2.) + pow(im_a, 2.); if(wavef == 1) piwavef = 1.; /* 1 = asymptotic distribution amplitude (default) */ if(wavef == 2) piwavef = 5./3.; /* 5/3 = CZ pion distribution amplitude */ return delta_perp / (PI * NU_M) * 18. * piwavef * ff_ps * im_a / ( sqr_a * (1. - pow(xi,2.)) - 81. * pow(piwavef,2.) * t / pow(2. * NU_M,2.) * pow(ff_ps,2.) - 18. * xi * piwavef * ff_ps * re_a ); } else return 0.; } double extrapol_spin_asymm(int reaction, int param, int wavef, double x_b, double t) { dcomplex direct, crossed; double delta_perp_sqr, delta_perp, ff_ps, xi, re_a, im_a, sqr_a, piwavef; delta_perp_sqr = -t * 4. * (1. - x_b) / pow(2. - x_b,2.) - pow(NU_M,2.) * 4. * pow(x_b,2.) / pow(2. - x_b,2.); if(delta_perp_sqr > 0.) { ff_ps = form_factor_G_P(t); xi = x_b / (2. - x_b); reaction_s4 = reaction; param_s4 = param; xi_s4 = xi; t_s4 = t; direct = Complex(int_gauss4(intd_axial_part, 0., 1., 30) + Hhel_p_distr(reaction,param,xi,xi,t) * log((1. - xi)/xi), - PI * Hhel_p_distr(reaction,param,xi, xi, t) ); crossed = Complex(int_gauss4(intc_axial_part, 0., 1., 30), 0.0); re_a = 6. * (direct.re + crossed.re); im_a = 6. * direct.im; sqr_a = pow(re_a, 2.) + pow(im_a, 2.); if(wavef == 1) piwavef = 1.; /* 1 = asymptotic distribution amplitude (default) */ if(wavef == 2) piwavef = 5./3.; /* 5/3 = CZ pion distribution amplitude */ return -2./9. * 1./(pow(PI_M,2.) * g_AXIAL) * (-t + pow(PI_M,2.)) * ( sqr_a * (1. - pow(xi,2.)) - 81. * pow(piwavef,2.) * t / pow(2. * NU_M,2.) * pow(ff_ps,2.) - 18. * xi * piwavef * ff_ps * re_a ) / (18. * piwavef * ff_ps); } else return 0.; } /*************************************************************************/ /*************************************************************************/ /* */ /* e + N -> e + pi + N : RESPONSE FUNCTIONS */ /* */ /*************************************************************************/ /*************************************************************************/ /* The function response_R_T() gives the response function R_T for the virtual photon + N -> pi + N process for pion and the kinematic variables Q_sqr, s, and t. */ double response_R_T(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result_plus = 0., result_minus = 0., unit_conv, cross; dcomplex contr_plus, contr_minus; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, -1, sp_in); result_plus += pow(Cabs(contr_plus),2.); result_minus += pow(Cabs(contr_minus),2.); } } choice = 1; switch(choice) { case 1 : /* for the T cross section in units (nanobarns / GeV^2) */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.25 * (result_plus + result_minus); break; } case 2 : /* for the R_T response function */ { return 0.25 * (result_plus + result_minus); break; } } } /* The function response_R_L() gives the response function R_L for the virtual photon + N -> pi + N process for pion and the kinematic variables Q_sqr, s, and t. */ double response_R_L(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result_zero = 0., resultspin = 0., unit_conv, cross, Lambda; dcomplex contr_zero; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_zero = J_hadron(reaction, param, Q_sqr, s, t, sp_out, 0, sp_in); result_zero += pow(Cabs(contr_zero),2.); } } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the L cross section in units (nanobarns / GeV^2) */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 * result_zero; break; } case 2 : /* for the L cross section in units (nanobarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); /* printf("\n%e \t %e", Q_sqr, (s - pow(NU_M,2.))/Lambda ); */ unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 * result_zero; break; } case 3 : /* for the R_L response function */ { return 0.5 * result_zero; break; } } } /* The function response_R_TT() gives the response function R_TT for the virtual photon + N -> pi + N process for pion and the kinematic variables Q_sqr, s, and t. */ double response_R_TT(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross; dcomplex contr; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr = Cmul(J_hadron(reaction, param, Q_sqr, s, t, sp_out, +1, sp_in), Conjg(J_hadron(reaction, param, Q_sqr, s, t, sp_out, -1, sp_in) ) ); result += contr.re; } } choice = 1; switch(choice) { case 1 : /* for the TT cross section in units (nanobarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 * result; break; } case 2 : /* for the R_TT response function */ { return 0.5 * result; break; } } } /* The function response_R_TL() gives the response function R_TL for the virtual photon + N -> pi + N process for pion and the kinematic variables Q_sqr, s, and t. */ double response_R_TL(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross; dcomplex contr_plus, contr_minus, contr_zero; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, -1, sp_in); contr_zero = J_hadron(reaction, param, Q_sqr, s, t, sp_out, 0, sp_in); result += Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).re; } } choice = 1; switch(choice) { case 1 : /* for the TL cross section in units (nanobarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the R_TL response function */ { return 0.5 / sqrt(2.) * result; break; } } } /* The function response_R_TLp() gives the response function R_TL' for the virtual photon + N -> pi + N process for pion and the kinematic variables Q_sqr, s, and t. */ double response_R_TLp(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, sp_in, result = 0., unit_conv, cross; dcomplex contr_plus, contr_minus, contr_zero; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { for(sp_in = -0.5; sp_in <= 0.5; sp_in += 1.0) { contr_plus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, +1, sp_in); contr_minus = J_hadron(reaction, param, Q_sqr, s, t, sp_out, -1, sp_in); contr_zero = J_hadron(reaction, param, Q_sqr, s, t, sp_out, 0, sp_in); result += Cmul(Csub(contr_plus, contr_minus), Conjg(contr_zero) ).im; } } choice = 1; switch(choice) { case 1 : /* for the TL' cross section in units (nanobarns / GeV^2). */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 / sqrt(2.) * result; break; } case 2 : /* for the R_TL' response function */ { return 0.5 / sqrt(2.) * result; break; } } } /* The function response_R_L_ortho() gives the response function R_L_ortho for the virtual photon + N -> pi + N process for kaon (1 = pi^+ n, 2 = pi^- p, 3 = pi^0 p) and the kinematic variables Q_sqr, s, and t. */ double response_R_L_ortho(int reaction, int param, double Q_sqr, double s, double t) { int choice; double sp_out, result = 0., unit_conv, cross, Lambda; dcomplex contr; for(sp_out = -0.5; sp_out <= 0.5; sp_out += 1.0) { contr = Cmul(J_hadron(reaction, param, Q_sqr, s, t, sp_out, 0, .5), Conjg(J_hadron(reaction, param, Q_sqr, s, t, sp_out, 0, -.5) ) ); result += 2.*contr.im; } choice = 2; /* default = 2 */ switch(choice) { case 1 : /* for the L cross section in units (nanobarns / GeV^2) */ { unit_conv = 1. / (16. * PI * pow(s - pow(NU_M, 2.),2.)) * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 * result; break; } case 2 : /* for the L cross section in units (nanobarns / GeV^2) (remark choice of virtual photon flux factor to compare to experimentally extracted reduced cross sections) */ { Lambda = sqrt( pow(s,2.) + pow(Q_sqr,2.) + pow(NU_M,4.) + 2. * s * Q_sqr - 2. * s * pow(NU_M,2.) + 2. * pow(NU_M,2.) * Q_sqr ); /* printf("\n%e \t %e", Q_sqr, (s - pow(NU_M,2.))/Lambda ); */ unit_conv = 1. / (16. * PI) * 1./(s - pow(NU_M, 2.)) * 1./Lambda * pow(0.197, 2.) * pow(10., 7.); return unit_conv * 0.5 * result; break; } case 3 : /* for the R_L response function */ { return 0.5 * result; break; } } } /*************************************************************************/ /*************************************************************************/ /* The relativistic hadronic matrixelement considered is given by J_hadron(). */ dcomplex J_hadron(int reaction, int param, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in) { return t_gap_ppi_hard(reaction, param, hardsoft_s, Q_sqr, s, t, sp_out, pol_ph, sp_in); } /*************************************************************************/ /* */ /* HARD PI ELECTROPRODUCTION AMPLITUDE */ /* */ /*************************************************************************/ dcomplex t_gap_ppi_hard(int reaction, int param, int hardsoft, double Q_sqr, double s, double t, double sp_out, int pol_ph, double sp_in) { int mu; double ga_mom_cm, ga_en_cm, mes_mom_cm, mes_en, mes_th, E_in, E_out, norm_in, norm_out, cos_th, xbj, zeta; dcomplex hard_piampl, soft_piampl, lorentz_soft, ampl, result; vec k_ga, p_mes, p_in, p_out; fvec k4_ga, p4_mes, pol4_ph; spinor spinor_p_in, spinor_p_out; ga_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) + Q_sqr, 2.) + pow(2. * NU_M, 2.) * Q_sqr ) / (4. * s) ); ga_en_cm = sqrt(s) - sqrt(pow(ga_mom_cm,2.) + pow(NU_M,2.)); k_ga = Vpolar(ga_mom_cm, 0., 0.); k4_ga = V4ini(Complex(ga_en_cm,0.), k_ga); mes_mom_cm = sqrt( ( pow(s - pow(NU_M,2.) - pow(PI_M,2.), 2.) - pow(2. * NU_M * PI_M, 2.) ) / (4. * s) ); mes_en = sqrt( pow(mes_mom_cm, 2.) + pow(PI_M, 2.) ); cos_th = (t + Q_sqr - pow(PI_M,2.) + 2. * ga_en_cm * mes_en) / (2. * ga_mom_cm * mes_mom_cm); if (cos_th > 1.) { /* printf("\n Problem cos_th > 1 : %e, for W = %f and Q^2 = %f", cos_th, sqrt(s), Q_sqr); */ cos_th = 1.; } mes_th = acos( cos_th ); p_mes = Vpolar(mes_mom_cm, mes_th, 0.); p4_mes = V4ini(Complex(mes_en,0.), p_mes); p_in = VRmul(k_ga, - 1.0); p_out = VRmul(p_mes,- 1.0); E_in = sqrt( pow( Vnorm(k_ga), 2.) + pow(NU_M,2.) ); E_out = sqrt( pow( Vnorm(p_mes), 2.) + pow(NU_M,2.) ); norm_in = sqrt(E_in + NU_M); norm_out = sqrt(E_out + NU_M); spinor_p_in = spinor_pos(p_in, NU_M, sp_in); spinor_p_out = spinor_pos(p_out, NU_M, sp_out); xbj = Q_sqr / (s - pow(NU_M,2.) + Q_sqr); zeta = xbj; if( (hardsoft == 1) || (hardsoft == 3) ) { hard_piampl = Jlong_pilong_prod_1gluon(reaction,param,sp_out,sp_in, p_in, p_out, k4_ga); } if( (hardsoft == 2) || (hardsoft == 3) ) { reaction_s4 = reaction; param_s4 = param; Q_sqr_s4 = Q_sqr; zeta_s4 = zeta; t_s4 = t; lorentz_soft = CRmul(1./sqrt(2.), spinleft_mat_spinright(spinor_p_out, SMmul(SMadd(dirac_gamma(0), dirac_gamma(3)), dirac_gamma5() ), spinor_p_in ) ); soft_piampl = CRmul(1./sqrt(2.) * sqrt(2.) * F_PI * norm_in * norm_out * ELEC * int_gauss4(int_soft_integral,0.,zeta,30), lorentz_soft ); } switch(hardsoft) { case 1 : /* 1-gluon exchange amplitude */ { ampl = hard_piampl; break; } case 2 : /* soft overlap amplitude */ { ampl = soft_piampl; break; } case 3 : /* 1-gluon exchange amplitude + soft overlap amplitude */ { printf("\n"); Cecho(hard_piampl); Cecho(soft_piampl); ampl = Cadd(hard_piampl, soft_piampl); break; } } return ampl; } dcomplex Jlong_pilong_prod_1gluon(int reaction, int param, double sp_out, double sp_in, vec p_in, vec p_out, fvec q4_in) { int wavefunction; dcomplex direct, crossed, lorentz_struct, piprod_axial, lorentz_ps_struct, piprod_ps; double p_in_en, p_out_en, Q_sqr, norm_in, norm_out, scale_sqr, pdotq, xbj, Pbdotq, xip, xi, delta_sqr, Mb_sqr, pi_charge, pqcd_piff, ff_hard, ff_soft, pi_cutoffsqr; fvec p4_in, p4_out, p4_av, delta, p4_sud, n4_sud; spinor spinor_p_in, spinor_p_out; p_in_en = sqrt(pow(Vnorm(p_in), 2.) + pow(NU_M, 2.)); p_out_en = sqrt(pow(Vnorm(p_out), 2.) + pow(NU_M, 2.)); p4_in = V4ini(Complex(p_in_en, 0.0), p_in); p4_out = V4ini(Complex(p_out_en, 0.0), p_out); p4_av = V4Rmul(V4add(p4_in,p4_out), 0.5); delta = V4sub(p4_out, p4_in); Q_sqr = CRmul(-1., V4mul(q4_in, q4_in)).re; pdotq = V4mul(p4_in,q4_in).re; xbj = Q_sqr / (2. * pdotq); delta_sqr = V4mul(delta,delta).re; Mb_sqr = pow(NU_M,2.) - 0.25 * delta_sqr; Pbdotq = V4mul(p4_av,q4_in).re; /* xip = 0.5 * Pbdotq / Mb_sqr * (-1. + sqrt(1. + Q_sqr * Mb_sqr / pow(Pbdotq,2.))); xi = xip * ( Q_sqr - delta_sqr + pow(PI_M,2.) ) / (Q_sqr + pow(2. * xip,2.) * Mb_sqr); printf("\n %e\t%e\t%e\t%e",xbj, 2.*xi, 2.*xip,xbj/(1. - xbj/2.)); */ xip = .5 * xbj/(1. - xbj/2.); p4_sud = V4Rmul(V4add(p4_av, V4Rmul(q4_in, - Mb_sqr * 2.*xip / Q_sqr) ), 1./( 1. + Mb_sqr * pow(2. * xip,2.) / Q_sqr ) ); n4_sud = V4Rmul(V4add(V4Rmul(p4_av, 2. * xip), q4_in), (4.* xip/Q_sqr)/( 1. + Mb_sqr * pow(2. * xip,2.) / Q_sqr ) ); norm_in = sqrt(p_in_en + NU_M); norm_out = sqrt(p_out_en + NU_M); spinor_p_in = spinor_pos(p_in, NU_M, sp_in); spinor_p_out = spinor_pos(p_out, NU_M, sp_out); if((parpipole != 2) && (parpipole != 4) && (parpipole != 6)) { lorentz_struct = CRmul(- 0.5 * 1./sqrt(Q_sqr), spinleft_mat_spinright(spinor_p_out, SMmul(fvec_slash(n4_sud), dirac_gamma5() ), spinor_p_in) ); reaction_s4 = reaction; param_s4 = param; xi_s4 = xip; t_s4 = delta_sqr; Q_sqr_s4 = Q_sqr; direct = Complex(int_gauss4(intd_axial_part, 0., 1., 30) + Hhel_p_distr(reaction,param,xip,xip,delta_sqr) * log((1. - xip)/xip), - PI * Hhel_p_distr(reaction,param,xip, xip, delta_sqr) ); crossed = Complex(int_gauss4(intc_axial_part, 0., 1., 30), 0.0); //TH - note F_PI is the pion decay constant - value used is 0.0924, defined in ../../h/ piprod_axial = Cmul(CRmul(4. * PI * alpha_strong(choice_alpha_s2,Q_sqr) * ELEC * 4./9. * (3. * sqrt(2.) * F_PI) * norm_in * norm_out, lorentz_struct ), Cadd(direct,crossed) ); } else piprod_axial = Complex(0.,0.); if((parpipole > 1) && (reaction != 0)) { if((parpipole == 2) || (parpipole == 3)) { //TH - test more reasonable parameterization here (12/5/06) pqcd_piff = 16.* PI * alpha_strong(choice_alpha_s2,Q_sqr) * pow(F_PI,2.) / Q_sqr; //pqcd_piff = 1./(1. + 1.86*Q_sqr + 0.12*pow(Q_sqr,2.)); } if((parpipole == 4) || (parpipole == 5)) { ff_hard = pion_emff(4, Q_sqr); ff_soft = pion_emff_soft_overlap(Q_sqr); pqcd_piff = ff_hard + ff_soft; } if((parpipole == 6) || (parpipole == 7)) { //TH (11/17/06) Lpi2=0.462 is charge radius from Bebek (r=0.7fm) //pi_cutoffsqr = 0.462; //TH (11/17/06) Lpi2=0.540 is the charge radius from Amendolia (r=0.661fm) pi_cutoffsqr = 0.540; pqcd_piff = 1./(1. + Q_sqr/pi_cutoffsqr); // TH (11/7/06) - more reasonable parameterization for Fpi - and interestingly does // not make much difference //pqcd_piff = 1./(1. + 1.77*Q_sqr + 0.12*pow(Q_sqr,2.)); //pi_cutoffsqr = Q_sqr*pqcd_piff/(1-pqcd_piff); //TH (11/17/06) - new parameterization including Fpi-1 //pqcd_piff = 1./(1. + 1.86*Q_sqr + 0.12*pow(Q_sqr,2.)); //pi_cutoffsqr = Q_sqr*pqcd_piff/(1-pqcd_piff); } lorentz_ps_struct = spinleft_mat_spinright(spinor_p_out, dirac_gamma5(), spinor_p_in); if(reaction == -1) pi_charge = -1.; if(reaction == +1) pi_charge = +1.; /*VGG1*/ /*piprod_ps = CRmul(pi_charge * ELEC * sqrt(Q_sqr) * pqcd_piff * sqrt(2.) * g_PINN / (- delta_sqr + pow(PI_M,2.)) * norm_in * norm_out, lorentz_ps_struct ); */ /*VGG2*/ /*cynthia 25/03/2004 new t parametrisation for E~*/ piprod_ps = CRmul(pi_charge * ELEC * sqrt(Q_sqr) * pqcd_piff * sqrt(2.) * g_PINN / (- delta_sqr + pow(PI_M,2.)) * norm_in * norm_out*exp(delta_sqr), lorentz_ps_struct ); } else piprod_ps = Complex(0.,0.); return Cadd(piprod_axial, piprod_ps); } double intd_axial_part(double x) { return ( Hhel_p_distr(reaction_s4, param_s4, x, xi_s4, t_s4) - Hhel_p_distr(reaction_s4, param_s4, xi_s4, xi_s4, t_s4) ) / (x - xi_s4); } double intc_axial_part(double x) { double ofpd; if(reaction_s4 == 0) ofpd = Hhel_p_distr(reaction_s4, param_s4, x, xi_s4, t_s4); if((reaction_s4 == -1) || (reaction_s4 == 1)) ofpd = - Hhel_p_distr(reaction_s4, param_s4, -x, xi_s4, t_s4); return ofpd / (x + xi_s4); } /*************************************************************************/ double int_soft_integral(double x) { double sigma_pi, sigma_nu; sigma_pi = 0.67; sigma_nu = sigma_pi; return nonforward_p_heldistr(reaction_s4, param_s4, x, zeta_s4, t_s4) * 8. * pow(PI,2.) /(sigma_pi + sigma_nu) * exp(- Q_sqr_s4/(2. * (sigma_pi + sigma_nu)) * (zeta_s4 - x)/x ); } /*************************************************************************/ /* */ /* OFF-FORWARD PARTON DISTRIBUTIONS */ /* that appear in PION electroproduction */ /* */ /*************************************************************************/ /* The function Hhel_p_distr() returns the off-forward valence polarized parton distribution Htilde(x,xi,t) of the proton that appears in pi production with 0 < x < 1, where Htilde(x,xi,t) is as defined by Ji (PRL 78 (1997) 610. */ double Hhel_p_distr(int reaction, int param, double x, double xi, double t) { double scale_sqr = 1.; /* this option is not active for polarized distr */ switch(reaction) { case -1 : { return (-1./3. * ( Hhel_up_distr(param,x,xi,t,scale_sqr) - Hhel_down_distr(param,x,xi,t,scale_sqr) ) -2./3. * ( Hhel_up_distr(param,-x,xi,t,scale_sqr) - Hhel_down_distr(param,-x,xi,t,scale_sqr) ) ); break; } case 0 : { return 1./sqrt(2.) * ( 2./3. * ( Hhel_up_distr(param,x,xi,t,scale_sqr) - Hhel_up_distr(param,-x,xi,t,scale_sqr) ) + 1./3.* ( Hhel_down_distr(param,x,xi,t,scale_sqr) - Hhel_down_distr(param,-x,xi,t,scale_sqr) ) ); break; } case 1 : { return (-2./3. * ( Hhel_up_distr(param,x,xi,t,scale_sqr) - Hhel_down_distr(param,x,xi,t,scale_sqr) ) -1./3. * ( Hhel_up_distr(param,-x,xi,t,scale_sqr) - Hhel_down_distr(param,-x,xi,t,scale_sqr) ) ); break; } } } /* The function nonforward_p_heldistr() returns the NON-forward valence polarized parton distribution of the proton that appears in pi production. */ double nonforward_p_heldistr(int reaction, int param, double x, double zeta, double t) { double scale_sqr = 1.; /* this option is not active for polarized distr */ double formfactor; formfactor = 1/g_AXIAL * form_factor_G_A(-t); switch(reaction) { case -1 : { return ( nonforward_heldistr(1,param,x,zeta,scale_sqr) - nonforward_heldistr(2,param,x,zeta,scale_sqr) ) * formfactor; break; } case 0 : { return 1./sqrt(2.) * ( 2./3. * nonforward_heldistr(1,param,x,zeta,scale_sqr) + 1./3. * nonforward_heldistr(2,param,x,zeta,scale_sqr) ) * formfactor; break; } case 1 : { return (- nonforward_heldistr(1,param,x,zeta,scale_sqr) + nonforward_heldistr(2,param,x,zeta,scale_sqr) ) * formfactor; break; } } }