* V. Guzey, July 2010, JLab * Nuclear diffractive parton distributions in the leading twist theory of nuclear shadowing * L. Frankfurt, V. Guzey and M. Strikman, Phys. Rept. (in preparation) * Suggested abbreviation: FGS10 * For the nucleus of Pb-208 * FOR ANY GIVEN VALUE OF POMERON LIGHT-CONE FRACTION X_P (10^{-4} < X < 0.1), * LIGHT-CONE FRACTION BETA (0.01 < BETA 1), AND Q^2 (4 < Q^2 < 16,000) GEV^2, * THE CODE GIVES THE RATIOS OF NUCLEAR TO FREE * NUCLEON NEXT-TO-LEADING ORDER (NLO) DIFFRACTIVE PDFS, * f_{j/A}^{D(3)}(x_Pom,beta,Q2)/[A*f_{j/N}^{D(3)}(x_Pom,beta,Q2)] * AND ALSO THE RATIO OF THE NUCLEAR TO NUCLEON NLO DIFFRACTIVE STRUCTURE FUNCTIONS * F_{2A}^{D(3)}(x_Pom,beta,Q2)/[A*F_{2N}^{D(3)}(x_Pom,beta,Q2)] * Note the flavor structure: ubar=dbar=sbar, charm is generated through the evolution * **** HOW TO USE THIS CODE *** * THE MAIN SOUBROUTINE IS CALLED '''LT_diffraction(xpom,beta,Q2,ubar,glue,f2,cbar)' * WHEN THE USER SPECIFIES THE VALUES OF xpom, beta, and Q2, * THE SUBROUTINE OUTPUTS THE FOLLOWING RATIOS: * * U_BAR(NUCLEUS)/[A*U_BAR(proton)] .... ubar (anti u-quark ratios) * Note that * U_BAR(NUCLEUS)/[A*U_BAR(proton)]=D_BAR(NUCLEUS)/[A*D_BAR(proton)] * =S_BAR(NUCLEUS)/[A*S_BAR(proton)] * GLUE(NUCLEUS)/[A*GLUE(proton)] .... glue (gluon ratios) * F_2^{D(3)}(NUCLEUS)/[A* _2^{D(3)}(proton)] .... f2 (F_2 ratios) * C_BAR(NUCLEUS)/[A*C_BAR(proton)] .... cbar (anti c-quark ratios) * **** WHAT DOES THE CODE DO? * THE CODE READS THE DATA FILE *( QCDEvolution_diffraction_XXp_modelX_Jul10.dat) AND * INTERPOLATES BETWEEN THE xpom, beta, and Q^2 GRID POINTS OF THE DATA FILE * (subroutine polin3) USING THE LINEAR INTERPOLATION. * THE DATA FILE CONTAINS 28 GRID POINTS IN XPOM, 19 POINTS IN BETA, * AND 7 GRID POINTS IN Q^2 (ACTUALLY, IN LOG(Q^2)/LOG(4.)) implicit double precision (a-h,o-z) dimension xpom_array(28),beta_array(19) data xpom_array /0.0001,0.0002,0.0003,0.0004,0.0005, >0.0006,0.0007,0.0008,0.0009,0.001,0.002,0.003,0.004, >0.005,0.006,0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05, >0.06,0.07,0.08,0.09,0.1/ data beta_array/0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08, >0.09,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0/ Q2=4.0 xpom=0.001 c do j=1,28 c xpom=xpom_array(j) do i=1,19 beta=beta_array(i) call LT_diffraction(xpom,beta,Q2,ubar,glue,f2,cbar) write(6,*) beta,ubar,glue,f2 7 format(9f8.5) enddo stop end subroutine LT_diffraction(xpom,beta,Q2,ubar,glue,f2,cbar) implicit double precision (a-h,o-z) integer iread dimension x1a(36),x2a(7),ya1(36,7,28),ya2(36,7,28), >ya3(36,7,28),ya4(36,7,28),res(4) common/fgs03/x1a,x2a,ya1,ya2,ya3,ya4 data iread/0/ if (iread.eq.0) then open(1,file= >'QCDEvolution_diffraction_pb208p_model2_Jul10.dat', >status='unknown') do k=1,28 do j=1,7 read(1,*)x2a(j) do i=1,36 read(1,*)x1a(i),ya1(i,j,k),ya2(i,j,k),ya3(i,j,k),ya4(i,j,k) enddo enddo enddo close(1) endif iread=1 call polin3(xpom,beta,q2,res) ubar=res(1) glue=res(2) f2=res(3) cbar=res(4) return end subroutine polin3(xpom,x1,x2,res) implicit double precision (a-h,o-z) integer j,k,iparton,NMAX,MMAX,Nxpom dimension x1a(36),x2a(7), > B(100),C(100),D(100),res(4) dimension yntmp1(7),ymtmp1(36),yntmp2(7), >ymtmp2(36),yntmp3(7),ymtmp3(36),yntmp4(7),ymtmp4(36) dimension ya1(36,7,28),ya2(36,7,28),ya3(36,7,28),ya4(36,7,28) dimension res1a(28),res2a(28),res3a(28),res4a(28) dimension xpom_array(28) common/fgs03/x1a,x2a,ya1,ya2,ya3,ya4 data xpom_array /0.0001,0.0002,0.0003,0.0004,0.0005, >0.0006,0.0007, >0.0008,0.0009,0.001,0.002,0.003,0.004,0.005,0.006, >0.007,0.008,0.009,0.01,0.02,0.03,0.04,0.05,0.06,0.07, >0.08,0.09,0.1 / parameter(NMAX=7,MMAX=36,Nxpom=28) x2a(1)=1. x2a(2)=2. x2a(3)=3. x2a(4)=4. x2a(5)=5. x2a(6)=6. x2a(7)=7. do ixpom=1,28 do j=1,36 do k=1,7 yntmp1(k)=ya1(j,k,ixpom) yntmp2(k)=ya2(j,k,ixpom) yntmp3(k)=ya3(j,k,ixpom) yntmp4(k)=ya4(j,k,ixpom) enddo x2d=Log(x2)/Log(4.) call linear(NMAX,x2a,yntmp1,x2d,ymtmp1(j)) call linear(NMAX,x2a,yntmp2,x2d,ymtmp2(j)) call linear(NMAX,x2a,yntmp3,x2d,ymtmp3(j)) call linear(NMAX,x2a,yntmp4,x2d,ymtmp4(j)) enddo call linear(MMAX,x1a,ymtmp1,x1,res1a(ixpom)) call linear(MMAX,x1a,ymtmp2,x1,res2a(ixpom)) call linear(MMAX,x1a,ymtmp3,x1,res3a(ixpom)) call linear(MMAX,x1a,ymtmp4,x1,res4a(ixpom)) enddo call linear(Nxpom,xpom_array,res1a,xpom,res(1)) call linear(Nxpom,xpom_array,res2a,xpom,res(2)) call linear(Nxpom,xpom_array,res3a,xpom,res(3)) call linear(Nxpom,xpom_array,res4a,xpom,res(4)) return end C Linear interpolation C Vadim Guzey, 06/10/2003 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