CMM C Modified the calculation of sigma to give it at any Q,W. C Need to interpolate the amplitudes in Q and fixed x_bj C C adjust nx and/or nq as necessary. Note: nq should stay preferably at 10! C C program main subroutine sig(qsqmin,qsqmax,wsq1,qbar) C C qsqmin = Q^2_min C qsqmax = Q^2_max C wsq1 = average W C qbar = average Q^2 C implicit none external readin, bour integer i,j integer iw,nw integer iq,mq integer k1,ncut integer nx,nq,ntemp,ntempq parameter(nx=20,nq=10) real*8 yy(nx), qq(nq),qq1(nq) real*8 wsq,w,wsq1 real*8 wmin, wmax real*8 qsqmin, qsqmax common /kins/ yy,qq C arrays and common blocks for the amplitudes of different flavours real*8 sigma, asumsq(nx,nq),tempx(nx),tempq(nq) common /amps2/ asumsq real*8 sig3 real*8 alpha, hcnb, bour, pi real*8 qbar,xbar,qbar1 real*8 logxb,res,err parameter (hcnb = 3.8937929D5) parameter (alpha = 7.2973525D-3) parameter (pi = 3.14159237D0) integer ordn C Set LO (1) or NLO (2) do k1 = 1,2 C read-in data from amplitude files CMM Read in data, depending whether LO (ordn=1) or NLO (ordn=2) ordn = k1 call readin(ordn) if (qsqmin.lt.qq(1).or.qsqmax.gt.qq(nq)) then print *,"Either Q^2_min too small or Q^2_max too large!" stop endif C open output file if (ordn.eq.1) then open(unit =51,file='sigwlo.dat',status='unknown') open(unit =52,file='sigloq.dat',status='unknown') else open(unit =51,file='sigwnlo.dat',status='unknown') open(unit =52,file='signloq.dat',status='unknown') endif C C sig vs. W for fixed Q^2 C C C depending on Q^2_min, determine number of points in Q^2 to be left out C in interpolation. C ncut = 0 do j=1,nq if (qq(j).lt.qsqmin) then ncut = ncut + 1 else goto 500 endif enddo 500 ntemp = nx C C if nq-ncut > 10, determine most advatageous region in Q^2 for C interpolation C if ((nq-ncut).gt.10) then ntempq = 10+ncut else ntempq = nq endif do i = 1,nx do j = ncut+1,ntempq tempq(j-ncut) = log(asumsq(i,j)) qq1(j-ncut) = qq(j) enddo call ratint(qq1,tempq,ntempq-ncut,qbar,res,err) tempx(i) = res enddo do iw = nx,1,-1 w = qbar*(1D0-yy(iw))/yy(iw) wsq = w**2 sig3 = pi*alpha**2*hcnb*dexp(tempx(iw)) > /(wsq*bour(qbar)) write(51,102) ordn, dsqrt(w), qbar, yy(iw), sig3 enddo CMM Fixed W plots vs. qsq ntemp = nx do i = 1,nx qbar1 = wsq1**2*yy(i)/(1D0-yy(i)) if (qbar1.ge.qsqmin.and.qbar1.le.qsqmax) then do j = ncut+1,ntempq tempq(j-ncut) = log(asumsq(i,j)) enddo call ratint(qq1,tempq,ntempq-ncut,qbar1,res,err) wsq = wsq1**2 sig3 = pi*alpha**2*hcnb*dexp(res)/(wsq**2*bour(qbar1)) write(52,102) ordn, wsq1, qbar1, yy(i), sig3 else goto 400 endif 400 enddo enddo 102 FORMAT(I1,1X,4(F12.5,1X)) 103 FORMAT(1X,11(F12.5,1X)) close(51) close(52) stop end CMM Input data files which contain re and imag parts of DVCS amplitudes CMM for the appropriate flavour and order CMM delta(j=1) CMM Q(i) ReA(j,i) ImA(j.i) CMM ... CMM delta(j=nx) CMM ... CMM Q(nq) ReA(nx,nq) ImA(nx.nq) subroutine readin(ordn) implicit none integer ordn integer nx,nq parameter(nx=20,nq=10) real*8 yy(nx),qq(nq) real*8 dum1,dum2 integer i,j real*8 x,q common /kins/ yy,qq real*8 reu(nx,nq),imu(nx,nq),red(nx,nq),imd(nx,nq) real*8 res(nx,nq),ims(nx,nq),reg(nx,nq),img(nx,nq) real*8 rsum(nx,nq), isum(nx,nq) real*8 asumsq(nx,nq) common /amps2/ asumsq CMM Read in data, depending whether LO (ordn=1) or NLO (ordn=2) C C del and q loops C if (ordn.eq.1) then open(unit =11,file='luamp.dat',status='unknown') do j = 1, nx read(11,*) x yy(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='ldamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lsamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='lgamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 reg(j,i)= dum1 img(j,i) = dum2 enddo enddo close(unit=11) else if (ordn.eq.2) then open(unit =11,file='nlouamp.dat',status='unknown') do j = 1, nx read(11,*) x yy(j) = x do i = 1,nq read(11,*) q,dum1,dum2 qq(i) = q reu(j,i) = dum1 imu(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='nlodamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 red(j,i) = dum1 imd(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='nlosamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 res(j,i)= dum1 ims(j,i) = dum2 enddo enddo close(unit=11) open(unit =11,file='nlogamp.dat',status='unknown') do j = 1, nx read(11,*) x do i = 1,nq read(11,*) q,dum1,dum2 reg(j,i) = dum1 img(j,i) = dum2 enddo enddo close(unit=11) else stop endif do j = 1, nx do i = 1,nq rsum(j,i) = reu(j,i) + red(j,i) + res(j,i) + reg(j,i) isum(j,i) = imu(j,i) + imd(j,i) + ims(j,i) + img(j,i) asumsq(j,i) = rsum(j,i)**2 + isum(j,i)**2 enddo enddo return 104 FORMAT(11(E12.4,1X)) end C C Simple model for B slope of t-dep.! One can put here whatever one wants! C real*8 function bour(q) implicit none real*8 q bour = 8D0*(1D0-0.15D0*log(q/2)) return end SUBROUTINE RATINT(XA,YA,N,X,Y,DY) IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (NMAX=11,TINY=1.E-25,MXX=1050) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 HH=ABS(X-XA(1)) DO 11 I=1,N H=ABS(X-XA(I)) IF (H.EQ.0.)THEN Y=YA(I) DY=0.0 RETURN ELSE IF (H.LT.HH) THEN NS=I HH=H ENDIF C(I)=YA(I) D(I)=YA(I)+TINY 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M W=C(I+1)-D(I) H=XA(I+M)-X T=(XA(I)-X)*D(I)/H DD=T-C(I+1) IF(DD.EQ.0.) GOTO 12 DD=W/DD D(I)=C(I+1)*DD C(I)=T*DD 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END