* Vadim Guzey, Theory group, Jefferson Lab * Email: vguzey@jlab.org * January, 2009 * Double precision FORTRAN code for the *PROTON CFFs* * [PDFs H AND E convoluted with * Leading Order hard coefficient functions] (H_script and E_script) * and summed with quark charges squared * (LO expressions entering DVCS observables) * in the DUAL parameterization of GPDs * V. Guzey and T. Teckentrup, Phys.Rev.D74 (2006) 054027 [hep-ph/0607099] * V. Guzey, T. Teckentrup, e-Print: arXiv:0810.3899 * IMPORTANT: Used non-factorizable, Regge-motivated model for the t-dependence * (see the paper for details). * For given values of xB (0.01 < xB < 0.5), Q**2 (1 < Q**2 < 10 GeV**2), * t (0.005 < |t| < 1 GeV**2) and for * -0.6 < J_u < 0.6 and -1 < J_d < 1 * the code gives the imaginary and real parts of H_script and E_script * HOW TO USE THE CODE * The code performs 3D linear interpolations using grids supplied by the * data files Hscript_dual_guidal.dat and Escript_dual_guidal_Ju**_Jd**.dat. *(There 49 data files for the CFF Hscript corresponding to 49 different * combinations of (Ju,Jd)) * Note that only the following values of (Ju,Jd) are allowed * Ju=(-0.6,-0.4,-0.2,0.0,0.2,0.4,0.6) * Jd=(-1.0,-0.6,-0.2,0.0,0.2,0.6,0.1) * * The data files should be stored in the same directory with the code. * An example of the use is given below: * The main subroutine Dual is called with given Ju (aJu), * Bjorken x (xB), Q**2 (Q2) and t. The subroutine outputs * Hscript_Im (the imaginary part of H_script), etc. implicit double precision (a-h,o-z) * Example aJu=0.0 aJd=0.0 xB=0.1 Q2=2.5 t=-0.1 * Calling the main subroutine call Dual_Proton(xB,Q2,t,aJu,aJd, >Hscript_Im_p,Hscript_Re_p,Escript_Im_p,Escript_Re_p) write(6,66)xB,Q2,t,Hscript_Im_p,Hscript_Re_p, >Escript_Im_p,Escript_Re_p 66 format(7f10.4) stop end C Main subroutine for CFFs of GPDs H and E C Input: xB = Bjorken x C Q2= Q**2 [GeV**2] C t = 4-momentum transfer [GeV**2] C (aJu,aJd) = (Ju,Jd) = proton total angular momentum fraction carries by (u,d) quarks C Output: Hscript_Im_p = imaginary part of the CFF H C Escript_Re_p = real part of the CFF H C Escript_Im_p = imaginary part of the CFF E C Escript_Re_p = real part of the CFF E subroutine Dual_Proton(xB,Q2,t,aJu,aJd, >Hscript_Im_p,Hscript_Re_p,Escript_Im_p,Escript_Re_p) implicit double precision (a-h,o-z) Parameter (N_E_data_files=49) Character FileName1(N_E_data_files)*40,FileName2*40 * Arrays of variables from the read data files dimension xB_array(27),Q2_array(19),t_array(13), >Hscript_Re_array(13,27,19),Hscript_Im_array(13,27,19) * Temporary arrays used for 3D interpolation dimension Q2tmp_Im(19),Q2tmp_Re(19),xBtmp_Im(27), >xBtmp_Re(27),ttmp_Im(13),ttmp_Re(13) * The same for the GPD E dimension Q2E_array(10),Escript_Re_array(13,27,10), >Escript_Im_array(13,27,10) dimension Q2Etmp_Im(10),Q2Etmp_Re(10),xBEtmp_Im(27), >xBEtmp_Re(27),tEtmp_Im(13),tEtmp_Re(13) data FileName1 /'Escript_dual_guidal_Ju06_Jd1_2009.dat', >'Escript_dual_guidal_Ju06_Jd06_2009.dat', >'Escript_dual_guidal_Ju06_Jd02_2009.dat', >'Escript_dual_guidal_Ju06_Jd00_2009.dat', >'Escript_dual_guidal_Ju06_Jdm02_2009.dat', >'Escript_dual_guidal_Ju06_Jdm06_2009.dat', >'Escript_dual_guidal_Ju06_Jdm1_2009.dat', >'Escript_dual_guidal_Ju04_Jd1_2009.dat', >'Escript_dual_guidal_Ju04_Jd06_2009.dat', >'Escript_dual_guidal_Ju04_Jd02_2009.dat', >'Escript_dual_guidal_Ju04_Jd00_2009.dat', >'Escript_dual_guidal_Ju04_Jdm02_2009.dat', >'Escript_dual_guidal_Ju04_Jdm06_2009.dat', >'Escript_dual_guidal_Ju04_Jdm1_2009.dat', >'Escript_dual_guidal_Ju02_Jd1_2009.dat', >'Escript_dual_guidal_Ju02_Jd06_2009.dat', >'Escript_dual_guidal_Ju02_Jd02_2009.dat', >'Escript_dual_guidal_Ju02_Jd00_2009.dat', >'Escript_dual_guidal_Ju02_Jdm02_2009.dat', >'Escript_dual_guidal_Ju02_Jdm06_2009.dat', >'Escript_dual_guidal_Ju02_Jdm1_2009.dat', >'Escript_dual_guidal_Ju00_Jd1_2009.dat', >'Escript_dual_guidal_Ju00_Jd06_2009.dat', >'Escript_dual_guidal_Ju00_Jd02_2009.dat', >'Escript_dual_guidal_Ju00_Jd00_2009.dat', >'Escript_dual_guidal_Ju00_Jdm02_2009.dat', >'Escript_dual_guidal_Ju00_Jdm06_2009.dat', >'Escript_dual_guidal_Ju00_Jdm1_2009.dat', >'Escript_dual_guidal_Jum02_Jd1_2009.dat', >'Escript_dual_guidal_Jum02_Jd06_2009.dat', >'Escript_dual_guidal_Jum02_Jd02_2009.dat', >'Escript_dual_guidal_Jum02_Jd00_2009.dat', >'Escript_dual_guidal_Jum02_Jdm02_2009.dat', >'Escript_dual_guidal_Jum02_Jdm06_2009.dat', >'Escript_dual_guidal_Jum02_Jdm1_2009.dat', >'Escript_dual_guidal_Jum04_Jd1_2009.dat', >'Escript_dual_guidal_Jum04_Jd06_2009.dat', >'Escript_dual_guidal_Jum04_Jd02_2009.dat', >'Escript_dual_guidal_Jum04_Jd00_2009.dat', >'Escript_dual_guidal_Jum04_Jdm02_2009.dat', >'Escript_dual_guidal_Jum04_Jdm06_2009.dat', >'Escript_dual_guidal_Jum04_Jdm1_2009.dat', >'Escript_dual_guidal_Jum06_Jd1_2009.dat', >'Escript_dual_guidal_Jum06_Jd06_2009.dat', >'Escript_dual_guidal_Jum06_Jd02_2009.dat', >'Escript_dual_guidal_Jum06_Jd00_2009.dat', >'Escript_dual_guidal_Jum06_Jdm02_2009.dat', >'Escript_dual_guidal_Jum06_Jdm06_2009.dat', >'Escript_dual_guidal_Jum06_Jdm1_2009.dat'/ common /Dual_Data/xB_array,Q2_array,t_array, >Hscript_Re_array,Hscript_Im_array,Q2E_array, >Escript_Re_array,Escript_Im_array C Reading the data file for CFF H open(1,file='Hscript_dual_guidal_2009.dat',status='unknown') Nt=13 Nx=27 Nq=19 do it=1,Nt do ix=1,Nx do iq=1,Nq read(1,*)xB_array(ix),Q2_array(iq),t_array(it), >Hscript_Im_array(it,ix,iq),Hscript_Re_array(it,ix,iq) enddo enddo enddo close(1) C Determining the appropriate data file for the given (J_u,J_d) C Go through 49 possibilities if (aJu.eq.0.6.and.aJd.eq.1.0) then FileName2=FileName1(1) else if (aJu.eq.0.6.and.aJd.eq.0.6) then FileName2=FileName1(2) else if (aJu.eq.0.6.and.aJd.eq.0.2) then FileName2=FileName1(3) else if (aJu.eq.0.6.and.aJd.eq.0.0) then FileName2=FileName1(4) else if (aJu.eq.0.6.and.aJd.eq.-0.2) then FileName2=FileName1(5) else if (aJu.eq.0.6.and.aJd.eq.-0.6) then FileName2=FileName1(6) else if (aJu.eq.0.6.and.aJd.eq.-1.0) then FileName2=FileName1(7) else if (aJu.eq.0.4.and.aJd.eq.1.0) then FileName2=FileName1(8) else if (aJu.eq.0.4.and.aJd.eq.0.6) then FileName2=FileName1(9) else if (aJu.eq.0.4.and.aJd.eq.0.2) then FileName2=FileName1(10) else if (aJu.eq.0.4.and.aJd.eq.0.0) then FileName2=FileName1(11) else if (aJu.eq.0.4.and.aJd.eq.-0.2) then FileName2=FileName1(12) else if (aJu.eq.0.4.and.aJd.eq.-0.6) then FileName2=FileName1(13) else if (aJu.eq.0.4.and.aJd.eq.-1.0) then FileName2=FileName1(14) else if (aJu.eq.0.2.and.aJd.eq.1.0) then FileName2=FileName1(15) else if (aJu.eq.0.2.and.aJd.eq.0.6) then FileName2=FileName1(16) else if (aJu.eq.0.2.and.aJd.eq.0.2) then FileName2=FileName1(17) else if (aJu.eq.0.2.and.aJd.eq.0.0) then FileName2=FileName1(18) else if (aJu.eq.0.2.and.aJd.eq.-0.2) then FileName2=FileName1(19) else if (aJu.eq.0.2.and.aJd.eq.-0.6) then FileName2=FileName1(20) else if (aJu.eq.0.2.and.aJd.eq.-1.0) then FileName2=FileName1(21) else if (aJu.eq.0.0.and.aJd.eq.1.0) then FileName2=FileName1(22) else if (aJu.eq.0.0.and.aJd.eq.0.6) then FileName2=FileName1(23) else if (aJu.eq.0.0.and.aJd.eq.0.2) then FileName2=FileName1(24) else if (aJu.eq.0.0.and.aJd.eq.0.0) then FileName2=FileName1(25) else if (aJu.eq.0.0.and.aJd.eq.-0.2) then FileName2=FileName1(26) else if (aJu.eq.0.0.and.aJd.eq.-0.6) then FileName2=FileName1(27) else if (aJu.eq.0.0.and.aJd.eq.-1.0) then FileName2=FileName1(28) else if (aJu.eq.-0.2.and.aJd.eq.1.0) then FileName2=FileName1(29) else if (aJu.eq.-0.2.and.aJd.eq.0.6) then FileName2=FileName1(30) else if (aJu.eq.-0.2.and.aJd.eq.0.2) then FileName2=FileName1(31) else if (aJu.eq.-0.2.and.aJd.eq.0.0) then FileName2=FileName1(32) else if (aJu.eq.-0.2.and.aJd.eq.-0.2) then FileName2=FileName1(33) else if (aJu.eq.-0.2.and.aJd.eq.-0.6) then FileName2=FileName1(34) else if (aJu.eq.-0.2.and.aJd.eq.-1.0) then FileName2=FileName1(35) else if (aJu.eq.-0.4.and.aJd.eq.1.0) then FileName2=FileName1(36) else if (aJu.eq.-0.4.and.aJd.eq.0.6) then FileName2=FileName1(37) else if (aJu.eq.-0.4.and.aJd.eq.0.2) then FileName2=FileName1(38) else if (aJu.eq.-0.4.and.aJd.eq.0.0) then FileName2=FileName1(39) else if (aJu.eq.-0.4.and.aJd.eq.-0.2) then FileName2=FileName1(40) else if (aJu.eq.-0.4.and.aJd.eq.-0.6) then FileName2=FileName1(41) else if (aJu.eq.-0.4.and.aJd.eq.-1.0) then FileName2=FileName1(42) else if (aJu.eq.-0.6.and.aJd.eq.1.0) then FileName2=FileName1(43) else if (aJu.eq.-0.6.and.aJd.eq.0.6) then FileName2=FileName1(44) else if (aJu.eq.-0.6.and.aJd.eq.0.2) then FileName2=FileName1(45) else if (aJu.eq.-0.6.and.aJd.eq.0.0) then FileName2=FileName1(46) else if (aJu.eq.-0.6.and.aJd.eq.-0.2) then FileName2=FileName1(47) else if (aJu.eq.-0.6.and.aJd.eq.-0.6) then FileName2=FileName1(48) else if (aJu.eq.-0.6.and.aJd.eq.-1.0) then FileName2=FileName1(49) else write(6,*) 'Ju and Jd out or range' endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif endif c write(6,*) 'J_u and J_d out of range' endif open(2,file=FileName2,status='unknown') Nt=13 Nx=27 NqE=10 do it=1,Nt do ix=1,Nx do iq=1,NqE read(2,*)a,Q2E_array(iq),a, >Escript_Im_array(it,ix,iq),Escript_Re_array(it,ix,iq) enddo enddo enddo close(2) * 3D interpolation done separately in all three variables do it=1,Nt do ix=1,Nx do iq=1,Nq Q2tmp_Im(iq)=Hscript_Im_array(it,ix,iq) Q2tmp_Re(iq)=Hscript_Re_array(it,ix,iq) enddo do iq=1,NqE Q2Etmp_Im(iq)=Escript_Im_array(it,ix,iq) Q2Etmp_Re(iq)=Escript_Re_array(it,ix,iq) enddo call linear(Nq,Q2_array,Q2tmp_Im,Q2,xBtmp_Im(ix)) call linear(Nq,Q2_array,Q2tmp_Re,Q2,xBtmp_Re(ix)) call linear(NqE,Q2E_array,Q2Etmp_Im,Q2,xBEtmp_Im(ix)) call linear(NqE,Q2E_array,Q2Etmp_Re,Q2,xBEtmp_Re(ix)) enddo call linear(Nx,xB_array,xBtmp_Im,xB,ttmp_Im(it)) call linear(Nx,xB_array,xBtmp_Re,xB,ttmp_Re(it)) call linear(Nx,xB_array,xBEtmp_Im,xB,tEtmp_Im(it)) call linear(Nx,xB_array,xBEtmp_Re,xB,tEtmp_Re(it)) enddo call linear(Nt,t_array,ttmp_Im,t,Hscript_Im_p) call linear(Nt,t_array,ttmp_Re,t,Hscript_Re_p) call linear(Nt,t_array,tEtmp_Im,t,Escript_Im_p) call linear(Nt,t_array,tEtmp_Re,t,Escript_Re_p) return end * 1D linear interpolation subroutine linear(n,xa,ya,x,y) implicit double precision (a-h,o-z) integer n,ns,i dimension xa(n),ya(n) C just in case there is some rounding error y=0.0 C Later added by T. Teckentrup ns=0 dif=abs(x-xa(1)) do i=1,n-1 C Finding the closest index i dift=abs(x-xa(i)) if (dift.le.dif) then ns=i dif=dift endif enddo C x could be greater or smaller than xa(i) if (x.gt.xa(ns)) then y1=ya(ns) y2=ya(ns+1) x1=xa(ns) x2=xa(ns+1) y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2) else y1=ya(ns-1) y2=ya(ns) x1=xa(ns-1) x2=xa(ns) y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2) endif return end