LHAPDF/lhapdf-5.9.1/src/wrapheragrid.f

   1       subroutine HERAGRIDevolve(xin,qin,pdf)
   2       implicit real*8 (a-h,o-z)
   3       include 'parmsetup.inc'
   4       double precision parm(nopmax)
   5       character*16 name(nmxset)
   6       integer nmem(nmxset),ndef(nmxset),mmem
   7       common/NAME/name,nmem,ndef,mmem
   8       double precision gridx(nmxgridx),gridq(nmxgridq)
   9       integer ngridx,ngridq,jx,jq
  10       CHARACTER*80 LINE
  11       dimension pdf(-6:6)
  12       integer init,set,i,j,k,l,nset,iset
  13       parameter(nhess=35)
  14       double precision fgrid(0:nhess,161,161,0:7),grid(322)
  15       double precision agrid(161),alfas,Qalfa
  16       common/fgridz/fgrid
  17       double precision up,dn,usea,dsea,str,chm,bot,glu
  18       double precision heragrid
  19       double precision qq(5),yntmp(5)
  20       real*8 mc,mc2,mb,mb2,mt,mt2,mz,mz2,alfa0,scale0
  21       COMMON/QCCONS/                                                    &
  22      &PI,PROTON,EUTRON,UCLEON,UDSCBT(6),AAM2H,BBM2H,AAM2L,BBM2L,        &
  23      &AAAR2,BBBR2,FL_FAC,CBMSTF(4:7),CHARGE(4:7),                       &
  24      &C1S3,C2S3,C4S3,C5S3,C8S3,C11S3,C14S3,C16S3,C20S3,C22S3,C28S3,     &
  25      &C38S3,C40S3,C44S3,C52S3,C136S3,C11S6,C2S9,C4S9,C10S9,C14S9,C16S9, &
  26      &C40S9,C44S9,C62S9,C112S9,C182S9,C11S12,C35S18,C61S12,C215S1,      &
  27      &C29S12,CPI2S3,CPIA,CPIB,CPIC,CPID,CPIE,CPIF,CCA,CCF,CTF,CATF,CFTF
  28       save
  29       x=xin
  30       q2=qin*qin
  31       call getnset(iset)
  32       up =   HERAGRID(x,Q2,grid,1)
  33       dn =   HERAGRID(x,Q2,grid,2)
  34       usea = HERAGRID(x,Q2,grid,3)
  35       dsea = HERAGRID(x,Q2,grid,4)
  36       if(name(iset)(1:9).eq.'HERAGRID1') then
  37         up = up + usea
  38         dn = dn + dsea
  39       endif
  40       str =  HERAGRID(x,Q2,grid,5)
  41       chm =  HERAGRID(x,Q2,grid,6)
  42       bot =  HERAGRID(x,Q2,grid,7)
  43       glu =  HERAGRID(x,Q2,grid,0)
  44       pdf(-6) = 0.0d0
  45        pdf(6) = 0.0d0
  46       pdf(-5) = bot
  47        pdf(5) = bot
  48       pdf(-4) = chm
  49        pdf(4) = chm
  50       pdf(-3) = str
  51        pdf(3) = str
  52       pdf(-2) = usea
  53        pdf(2) = up
  54       pdf(-1) = dsea
  55        pdf(1) = dn
  56        pdf(0) = glu
  57       return
  58 !
  59       entry HERAGRIDgetgrid(nset,ngridx,ngridq,gridx,gridq)
  60
  61       ngridx=161
  62       do jx=1,ngridx
  63           gridx(jx)=grid(jx)
  64       enddo
  65       ngridq=161
  66       do jq=1,ngridq
  67           gridq(jq)=grid(jq+161)
  68       enddo
  69       return
  70 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  71       entry HERAGRIDread(nset)
  72 !
  73       call getnmem(nset,imem)
  74
  75       read(1,*)nmem(nset),ndef(nset)
  76
  77       do iq=1,23
  78          read(1,*) (grid(161+(iq-1)*7+ii),ii=1,7)
  79       enddo
  80
  81       do jx=1,23
  82          read(1,*) (grid((jx-1)*7+ii),ii=1,7)
  83       enddo
  84
  85     !read in to alphas grid
  86       do iq=1,23
  87          read(1,*) (agrid((iq-1)*7+ii),ii=1,7)
  88       enddo
  89
  90
  91       do ns=0,nmem(nset)
  92         do k = 0,7
  93           do IQ=1,161
  94             do JX=1,23
  95               read(1,*)(fgrid(ns,(jx-1)*7+ii,iq,k),ii=1,7)
  96             enddo
  97           enddo
  98         enddo
  99       enddo
 100       return
 101 !
 102 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 103       entry HERAGRIDalfa(alfas,Qalfa)
 104
 105        call getnset(iset)
 106        q2 = Qalfa*Qalfa
 107        nf=6
 108        if (Q2.lt.mt2) nf=5
 109        if (Q2.lt.mb2) nf=4
 110        if (Q2.lt.mc2) nf=3
 111        call GetOrderAsM(iset,iord)
 112        iord=iord+1
 113        call listPDF(iset,imem,parm)
 114        alfa0=parm(1)
 115        alfas=A0TOA1(q2,mz2,alfa0,iord,nf,ierr)
 116
 117       return
 118 !
 119 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 120       entry HERAGRIDinit(Eorder,Q2fit)
 121
 122       call getnset(iset)
 123       mz = 91.187d0
 124       mz2=mz*mz
 125       call getQmassM(iset,4,mc)
 126       call getQmassM(iset,5,mb)
 127       call getQmassM(iset,6,mt)
 128       mc2=mc*mc
 129       mb2=mb*mb
 130       mt2=mt*mt
 131       UDSCBT(1) = 0.005
 132       UDSCBT(2) = 0.01
 133       UDSCBT(3) = 0.3
 134       UDSCBT(4) = mc
 135       UDSCBT(5) = mb
 136       UDSCBT(6) = mt
 137
 138       return
 139 !
 140 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 141       entry HERAGRIDpdf(mem)
 142       imem = mem
 143       call getnset(iset)
 144       call setnmem(iset,imem)
 145       return
 146 !
 147  1000 format(5e13.5)
 148       end
 149 !
 150 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 151 ! attempt at polynomial (order 4) interpolation based on  polint
 152       double precision function HERAGRID(x,Q2,grid,n)
 153       implicit real*8(a-h,o-z)
 154       integer iset,imem,nhess
 155       parameter(nhess=35)
 156       double precision grid(322),x,Q2
 157       double precision fgrid(0:nhess,161,161,0:7)
 158       double precision ya(5,5),yntmp(5),ymtmp(5)
 159       common/fgridz/fgrid
 160
 161       npt = 5
 162       call getnset(iset)
 163       call getnmem(iset,imem)
 164
 165 ! find the x bins around x
 166       nxlow = -1
 167       nxhi = 162
 168       do while (nxhi-nxlow.gt.1)
 169          nxmid = (nxlow+nxhi)/2
 170          if(x.ge.grid(nxmid)) then
 171             nxlow = nxmid
 172          else
 173             nxhi = nxmid
 174          endif
 175       enddo
 176
 177 ! find the q2 bins around q2
 178       nqlow = -1
 179       nqhi = 162
 180       do while (nqhi-nqlow.gt.1)
 181          nqmid = (nqlow+nqhi)/2
 182          if(q2.ge.grid(161+nqmid)) then
 183             nqlow = nqmid
 184          else
 185             nqhi = nqmid
 186          endif
 187       enddo
 188
 189 ! fill the temp 4x4 funtion array (allowing for endpoints and extrapolation)
 190         if(nxlow.le.0) nxlow=1
 191         if(nxlow.ge.161) nxlow=160
 192         if(nxlow.eq.1) then
 193            nxbot = 1
 194         elseif(nxlow.eq.2) then
 195            nxbot = 2
 196         else
 197            if(nxlow.eq.160) then
 198               nxbot = 4
 199            else
 200               nxbot = 3
 201            endif
 202         endif
 203         if(nqlow.le.0) nqlow=1
 204         if(nqlow.ge.161) nqlow=160
 205         if(nqlow.eq.1) then
 206            nqbot = 1
 207         elseif(nqlow.eq.2) then
 208            nqbot = 2
 209         else
 210            if(nqlow.eq.160) then
 211               nqbot = 4
 212            else
 213               nqbot = 3
 214            endif
 215         endif
 216
 217         do nx=1,5
 218           do nq=1,5
 219               ya(nx,nq) = fgrid(imem,nxlow+nx-nxbot,nqlow+nq-nqbot,n)
 220           enddo
 221         enddo
 222
 223       do j=1,5
 224         do k=1,5
 225           yntmp(k)=ya(j,k)
 226         enddo
 227         call herapolint(grid(161+nqlow-nqbot+1),yntmp,npt,q2,ymtmp(j),dy)
 228       enddo
 229       call herapolint(grid(nxlow-nxbot+1),ymtmp,npt,x,y,dy)
 230
 231       heragrid=y
 232
 233       return
 234       end
 235 !=========================================================
 236          SUBROUTINE HERAPOLINT (XA,YA,N,X,Y,DY)
 237
 238       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
 239 !                                        Adapted from "Numerical Recipes"
 240       PARAMETER (NMAX=10)
 241       DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
 242       NS=1
 243       DIF=ABS(X-XA(1))
 244       DO 11 I=1,N
 245         DIFT=ABS(X-XA(I))
 246         IF (DIFT.LT.DIF) THEN
 247           NS=I
 248           DIF=DIFT
 249         ENDIF
 250         C(I)=YA(I)
 251         D(I)=YA(I)
 252 11    CONTINUE
 253       Y=YA(NS)
 254       NS=NS-1
 255       DO 13 M=1,N-1
 256         DO 12 I=1,N-M
 257           HO=XA(I)-X
 258           HP=XA(I+M)-X
 259           W=C(I+1)-D(I)
 260           DEN=HO-HP
 261           IF(DEN.EQ.0.) RETURN
 262           DEN=W/DEN
 263           D(I)=HP*DEN
 264           C(I)=HO*DEN
 265 12      CONTINUE
 266         IF (2*NS.LT.N-M)THEN
 267           DY=C(NS+1)
 268         ELSE
 269           DY=D(NS)
 270           NS=NS-1
 271         ENDIF
 272         Y=Y+DY
 273 13    CONTINUE
 274       RETURN
 275       END
 276 !=========================================================