* Vadim Guzey, Theory group, Jefferson Lab * Email: vguzey@jlab.org * January, 2009 * Double precision FORTRAN code for the *NEUTRON 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_Neutron(xB,Q2,t,aJu,aJd, >Hscript_Im_n,Hscript_Re_n,Escript_Im_n,Escript_Re_n) write(6,66)xB,Q2,t,Hscript_Im_n,Hscript_Re_n, >Escript_Im_n,Escript_Re_n 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_n = imaginary part of the CFF H C Escript_Re_n = real part of the CFF H C Escript_Im_n = imaginary part of the CFF E C Escript_Re_n = real part of the CFF E subroutine Dual_Neutron(xB,Q2,t,aJu,aJd, >Hscript_Im_n,Hscript_Re_n,Escript_Im_n,Escript_Re_n) implicit double precision (a-h,o-z) Parameter (N_E_data_files=49) Character FileName1_n(N_E_data_files)*46, >FileName2_n*46 * Arrays of variables from the read data files dimension xB_array_n(27),Q2_array_n(19),t_array_n(13), >Hscript_Re_array_n(13,27,19),Hscript_Im_array_n(13,27,19) * Temporary arrays used for 3D interpolation dimension Q2tmp_Im_n(19),Q2tmp_Re_n(19),xBtmp_Im_n(27), >xBtmp_Re_n(27),ttmp_Im_n(13),ttmp_Re_n(13) * The same for the GPD E dimension Q2E_array_n(10),Escript_Re_array_n(13,27,10), >Escript_Im_array_n(13,27,10) dimension Q2Etmp_Im_n(10),Q2Etmp_Re_n(10),xBEtmp_Im_n(27), >xBEtmp_Re_n(27),tEtmp_Im_n(13),tEtmp_Re_n(13) data FileName1_n /'Escript_dual_guidal_Ju06_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Ju06_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Ju04_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Ju02_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Ju00_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Jum02_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Jum04_Jdm1_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jd1_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jd06_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jd02_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jd00_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jdm02_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jdm06_neutron_2009.dat', >'Escript_dual_guidal_Jum06_Jdm1_neutron_2009.dat'/ common /Dual_Data/xB_array_n,Q2_array_n,t_array_n, >Hscript_Re_array_n,Hscript_Im_array_n,Q2E_array_n, >Escript_Re_array_n,Escript_Im_array_n C Reading the data file for CFF H open(11,file='Hscript_dual_guidal_neutron_2009.dat', >status='unknown') Nt=13 Nx=27 Nq=19 do it=1,Nt do ix=1,Nx do iq=1,Nq read(11,*)xB_array_n(ix),Q2_array_n(iq),t_array_n(it), >Hscript_Im_array_n(it,ix,iq),Hscript_Re_array_n(it,ix,iq) enddo enddo enddo close(11) 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_n=FileName1_n(1) else if (aJu.eq.0.6.and.aJd.eq.0.6) then FileName2_n=FileName1_n(2) else if (aJu.eq.0.6.and.aJd.eq.0.2) then FileName2_n=FileName1_n(3) else if (aJu.eq.0.6.and.aJd.eq.0.0) then FileName2_n=FileName1_n(4) else if (aJu.eq.0.6.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(5) else if (aJu.eq.0.6.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(6) else if (aJu.eq.0.6.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(7) else if (aJu.eq.0.4.and.aJd.eq.1.0) then FileName2_n=FileName1_n(8) else if (aJu.eq.0.4.and.aJd.eq.0.6) then FileName2_n=FileName1_n(9) else if (aJu.eq.0.4.and.aJd.eq.0.2) then FileName2_n=FileName1_n(10) else if (aJu.eq.0.4.and.aJd.eq.0.0) then FileName2_n=FileName1_n(11) else if (aJu.eq.0.4.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(12) else if (aJu.eq.0.4.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(13) else if (aJu.eq.0.4.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(14) else if (aJu.eq.0.2.and.aJd.eq.1.0) then FileName2_n=FileName1_n(15) else if (aJu.eq.0.2.and.aJd.eq.0.6) then FileName2_n=FileName1_n(16) else if (aJu.eq.0.2.and.aJd.eq.0.2) then FileName2_n=FileName1_n(17) else if (aJu.eq.0.2.and.aJd.eq.0.0) then FileName2_n=FileName1_n(18) else if (aJu.eq.0.2.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(19) else if (aJu.eq.0.2.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(20) else if (aJu.eq.0.2.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(21) else if (aJu.eq.0.0.and.aJd.eq.1.0) then FileName2_n=FileName1_n(22) else if (aJu.eq.0.0.and.aJd.eq.0.6) then FileName2_n=FileName1_n(23) else if (aJu.eq.0.0.and.aJd.eq.0.2) then FileName2_n=FileName1_n(24) else if (aJu.eq.0.0.and.aJd.eq.0.0) then FileName2_n=FileName1_n(25) else if (aJu.eq.0.0.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(26) else if (aJu.eq.0.0.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(27) else if (aJu.eq.0.0.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(28) else if (aJu.eq.-0.2.and.aJd.eq.1.0) then FileName2_n=FileName1_n(29) else if (aJu.eq.-0.2.and.aJd.eq.0.6) then FileName2_n=FileName1_n(30) else if (aJu.eq.-0.2.and.aJd.eq.0.2) then FileName2_n=FileName1_n(31) else if (aJu.eq.-0.2.and.aJd.eq.0.0) then FileName2_n=FileName1_n(32) else if (aJu.eq.-0.2.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(33) else if (aJu.eq.-0.2.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(34) else if (aJu.eq.-0.2.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(35) else if (aJu.eq.-0.4.and.aJd.eq.1.0) then FileName2_n=FileName1_n(36) else if (aJu.eq.-0.4.and.aJd.eq.0.6) then FileName2_n=FileName1_n(37) else if (aJu.eq.-0.4.and.aJd.eq.0.2) then FileName2_n=FileName1_n(38) else if (aJu.eq.-0.4.and.aJd.eq.0.0) then FileName2_n=FileName1_n(39) else if (aJu.eq.-0.4.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(40) else if (aJu.eq.-0.4.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(41) else if (aJu.eq.-0.4.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(42) else if (aJu.eq.-0.6.and.aJd.eq.1.0) then FileName2_n=FileName1_n(43) else if (aJu.eq.-0.6.and.aJd.eq.0.6) then FileName2_n=FileName1_n(44) else if (aJu.eq.-0.6.and.aJd.eq.0.2) then FileName2_n=FileName1_n(45) else if (aJu.eq.-0.6.and.aJd.eq.0.0) then FileName2_n=FileName1_n(46) else if (aJu.eq.-0.6.and.aJd.eq.-0.2) then FileName2_n=FileName1_n(47) else if (aJu.eq.-0.6.and.aJd.eq.-0.6) then FileName2_n=FileName1_n(48) else if (aJu.eq.-0.6.and.aJd.eq.-1.0) then FileName2_n=FileName1_n(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 endif open(12,file=FileName2_n,status='unknown') Nt=13 Nx=27 NqE=10 do it=1,Nt do ix=1,Nx do iq=1,NqE read(12,*)a,Q2E_array_n(iq),a, >Escript_Im_array_n(it,ix,iq),Escript_Re_array_n(it,ix,iq) enddo enddo enddo close(12) * 3D interpolation done separately in all three variables do it=1,Nt do ix=1,Nx do iq=1,Nq Q2tmp_Im_n(iq)=Hscript_Im_array_n(it,ix,iq) Q2tmp_Re_n(iq)=Hscript_Re_array_n(it,ix,iq) enddo do iq=1,NqE Q2Etmp_Im_n(iq)=Escript_Im_array_n(it,ix,iq) Q2Etmp_Re_n(iq)=Escript_Re_array_n(it,ix,iq) enddo call linear(Nq,Q2_array_n,Q2tmp_Im_n,Q2,xBtmp_Im_n(ix)) call linear(Nq,Q2_array_n,Q2tmp_Re_n,Q2,xBtmp_Re_n(ix)) call linear(NqE,Q2E_array_n,Q2Etmp_Im_n,Q2,xBEtmp_Im_n(ix)) call linear(NqE,Q2E_array_n,Q2Etmp_Re_n,Q2,xBEtmp_Re_n(ix)) enddo call linear(Nx,xB_array_n,xBtmp_Im_n,xB,ttmp_Im_n(it)) call linear(Nx,xB_array_n,xBtmp_Re_n,xB,ttmp_Re_n(it)) call linear(Nx,xB_array_n,xBEtmp_Im_n,xB,tEtmp_Im_n(it)) call linear(Nx,xB_array_n,xBEtmp_Re_n,xB,tEtmp_Re_n(it)) enddo call linear(Nt,t_array_n,ttmp_Im_n,t,Hscript_Im_n) call linear(Nt,t_array_n,ttmp_Re_n,t,Hscript_Re_n) call linear(Nt,t_array_n,tEtmp_Im_n,t,Escript_Im_n) call linear(Nt,t_array_n,tEtmp_Re_n,t,Escript_Re_n) 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