LHAPDF/lhapdf5.3.1/wrapmrst.f

   1       subroutine MRSTevolve(x,Q,pdf)
   2       implicit real*8(a-h,o-z)
   3       include 'parmsetup.inc'
   4       character*16 name(nmxset)
   5       integer nmem(nmxset),ndef(nmxset),mmem
   6 c      integer member(nmxset)
   7       integer nset,iset
   8       common/NAME/name,nmem,ndef,mmem
   9       parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26,
  10      .nhess=30)
  11       real*8 pdf(-6:6)
  12       real*8 f1(nx,nq)
  13      .,f2(nx,nq)
  14      .,f3(nx,nq)
  15      .,f4(nx,nq)
  16      .,f5(nx,nq)
  17      .,f6(nx,nq)
  18      .,f7(nx,nq)
  19      .,f8(nx,nq)
  20      .,fc(nx,nqc),fb(nx,nqb)
  21       real*8 qq(nq),xx(nx),
  22      .cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4),
  23      .cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),
  24      .cc6(0:nhess,nx,nq,4,4),cc8(0:nhess,nx,nq,4,4),
  25      .ccc(0:nhess,nx,nqc,4,4),ccb(0:nhess,nx,nqb,4,4)
  26       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
  27       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
  28      .        1d-4,2d-4,4d-4,6d-4,8d-4,
  29      .        1d-3,2d-3,4d-3,6d-3,8d-3,
  30      .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
  31      .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
  32      .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
  33      .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
  34      .     .8d0,.9d0,1d0/
  35       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
  36      .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
  37      .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
  38      .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
  39      .        1.8d6,3.2d6,5.6d6,1d7/
  40       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
  41       save
  42
  43       xsave=x
  44       qsq = q*q
  45       q2save=qsq
  46
  47       xlog=dlog(x)
  48       qsqlog=dlog(qsq)
  49
  50       call getnset(iset)
  51 c      imem=member(iset)
  52       call getnmem(iset,imem)
  53       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
  54       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
  55       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
  56       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
  57       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
  58       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
  59
  60       chm=0.d0
  61       if(qsq.gt.emc2) then
  62       call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
  63       endif
  64
  65       bot=0.d0
  66       if(qsq.gt.emb2) then
  67       call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
  68       endif
  69
  70       pdf(0)  = glu
  71       pdf(1)  = dnv+dsea
  72       pdf(-1) = dsea
  73       pdf(2)  = upv+usea
  74       pdf(-2) = usea
  75       pdf(3)  = str
  76       pdf(-3) = str
  77       pdf(4)  = chm
  78       pdf(-4) = chm
  79       pdf(5)  = bot
  80       pdf(-5) = bot
  81       pdf(6)  = 0.0d0
  82       pdf(-6) = 0.0d0
  83
  84       x=xsave
  85       qsq=q2save
  86       return
  87 *
  88       entry MRSTread(nset)
  89       read(1,*)nmem(nset),ndef(nset)
  90 c      print *,nmem(nset),ndef(nset)
  91 c      do nm = 0,nmem-1
  92       do nm = 0,nmem(nset)
  93         do 20 n=1,nx-1
  94         do 20 m=1,nq
  95         read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
  96      .            f5(n,m),f7(n,m),f6(n,m),f8(n,m)
  97 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
  98   20  continue
  99 c      write(*,*)'PDF set ',nm,' first element ',f1(1,1)
 100       do 40 m=1,nq
 101       f1(nx,m)=0.d0
 102       f2(nx,m)=0.d0
 103       f3(nx,m)=0.d0
 104       f4(nx,m)=0.d0
 105       f5(nx,m)=0.d0
 106       f6(nx,m)=0.d0
 107       f7(nx,m)=0.d0
 108       f8(nx,m)=0.d0
 109   40  continue
 110
 111       do n=1,nx
 112       xxl(n)=dlog(xx(n))
 113       enddo
 114       do m=1,nq
 115       qql(m)=dlog(qq(m))
 116       enddo
 117
 118       call jeppe1(nm,nx,nq,xxl,qql,f1,cc1)
 119       call jeppe1(nm,nx,nq,xxl,qql,f2,cc2)
 120       call jeppe1(nm,nx,nq,xxl,qql,f3,cc3)
 121       call jeppe1(nm,nx,nq,xxl,qql,f4,cc4)
 122       call jeppe1(nm,nx,nq,xxl,qql,f6,cc6)
 123       call jeppe1(nm,nx,nq,xxl,qql,f8,cc8)
 124
 125       emc2=2.045
 126       emb2=18.5
 127
 128       do 44 m=1,nqc
 129       qqlc(m)=qql(m+nqc0)
 130       do 44 n=1,nx
 131       fc(n,m)=f5(n,m+nqc0)
 132    44 continue
 133       qqlc(1)=dlog(emc2)
 134       call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc)
 135
 136       do 45 m=1,nqb
 137       qqlb(m)=qql(m+nqb0)
 138       do 45 n=1,nx
 139       fb(n,m)=f7(n,m+nqb0)
 140    45 continue
 141       qqlb(1)=dlog(emb2)
 142       call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb)
 143
 144
 145       enddo
 146   50  format(8f10.5)
 147       return
 148 *
 149
 150       entry MRSTalfa(nflav,alfas,Qalfa)
 151         call getnset(iset)
 152 c        mem = member(iset)
 153 c       call setnmem(member(iset))
 154         call alphamrs(nflav,alfas,Qalfa)
 155       return
 156 *
 157       entry MRSTinit(Eorder,Q2fit)
 158       return
 159 *
 160       entry MRSTpdf(mem)
 161 c      if(mem.eq.0) mem=ndef
 162 c      imem = mem
 163         call getnset(iset)
 164 c       member(iset)=mem
 165         call setnmem(iset,mem)
 166
 167       return
 168 *
 169       end
 170
 171       subroutine jeppe1(imem,nx,my,xx,yy,ff,cc)
 172       implicit real*8(a-h,o-z)
 173       parameter(nnx=49,mmy=37,nhess=30)
 174       dimension xx(nx),yy(my),ff(nnx,mmy),
 175      xff1(nnx,mmy),ff2(nnx,mmy),
 176      xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
 177      xcl(16),cc(0:nhess,nx,my,4,4),iwt(16,16)
 178
 179       data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
 180      x            0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
 181      x            -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
 182      x            2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
 183      x            0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
 184      x            0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
 185      x            0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
 186      x            0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
 187      x            -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
 188      x            0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
 189      x            9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
 190      x            -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
 191      x            2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
 192      x            0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
 193      x            -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
 194      x            4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
 195
 196
 197       do 42 m=1,my
 198       dx=xx(2)-xx(1)
 199       ff1(1,m)=(ff(2,m)-ff(1,m))/dx
 200       dx=xx(nx)-xx(nx-1)
 201       ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
 202       do 41 n=2,nx-1
 203       ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
 204      xff(n+1,m))
 205    41 continue
 206    42 continue
 207
 208       do 44 n=1,nx
 209       dy=yy(2)-yy(1)
 210       ff2(n,1)=(ff(n,2)-ff(n,1))/dy
 211       dy=yy(my)-yy(my-1)
 212       ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
 213       do 43 m=2,my-1
 214       ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
 215      xff(n,m+1))
 216    43 continue
 217    44 continue
 218
 219       do 46 m=1,my
 220       dx=xx(2)-xx(1)
 221       ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
 222       dx=xx(nx)-xx(nx-1)
 223       ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
 224       do 45 n=2,nx-1
 225       ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
 226      xff2(n+1,m))
 227    45 continue
 228    46 continue
 229
 230       do 53 n=1,nx-1
 231       do 52 m=1,my-1
 232       d1=xx(n+1)-xx(n)
 233       d2=yy(m+1)-yy(m)
 234       d1d2=d1*d2
 235
 236       yy0(1)=ff(n,m)
 237       yy0(2)=ff(n+1,m)
 238       yy0(3)=ff(n+1,m+1)
 239       yy0(4)=ff(n,m+1)
 240
 241       yy1(1)=ff1(n,m)
 242       yy1(2)=ff1(n+1,m)
 243       yy1(3)=ff1(n+1,m+1)
 244       yy1(4)=ff1(n,m+1)
 245
 246       yy2(1)=ff2(n,m)
 247       yy2(2)=ff2(n+1,m)
 248       yy2(3)=ff2(n+1,m+1)
 249       yy2(4)=ff2(n,m+1)
 250
 251       yy12(1)=ff12(n,m)
 252       yy12(2)=ff12(n+1,m)
 253       yy12(3)=ff12(n+1,m+1)
 254       yy12(4)=ff12(n,m+1)
 255
 256       do 47 k=1,4
 257       z(k)=yy0(k)
 258       z(k+4)=yy1(k)*d1
 259       z(k+8)=yy2(k)*d2
 260       z(k+12)=yy12(k)*d1d2
 261    47 continue
 262
 263       do 49 l=1,16
 264       xxd=0.
 265       do 48 k=1,16
 266       xxd=xxd+iwt(k,l)*z(k)
 267    48 continue
 268       cl(l)=xxd
 269    49 continue
 270       l=0
 271       do 51 k=1,4
 272       do 50 j=1,4
 273       l=l+1
 274       cc(imem,n,m,k,j)=cl(l)
 275    50 continue
 276    51 continue
 277    52 continue
 278    53 continue
 279       return
 280       end
 281
 282       subroutine jeppe2(i,x,y,nx,my,xx,yy,cc,z)
 283 C--   G.W. 02/07/2007 Allow extrapolation to small x and large q.
 284       implicit real*8(a-h,o-z)
 285       parameter(nhess=30)
 286       dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)
 287
 288       n=locx(xx,nx,x)
 289       m=locx(yy,my,y)
 290
 291       if (n.gt.0.and.n.lt.nx.and.m.gt.0.and.m.lt.my) then
 292 C--   Do usual interpolation.
 293          t=(x-xx(n))/(xx(n+1)-xx(n))
 294          u=(y-yy(m))/(yy(m+1)-yy(m))
 295          z=0.d0
 296          do l=4,1,-1
 297             z=t*z+((cc(i,n,m,l,4)*u+cc(i,n,m,l,3))*u
 298      &           +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
 299          enddo
 300
 301       else if (n.eq.0.and.m.gt.0.and.m.lt.my) then
 302 C--   Extrapolate to small x.
 303          call jeppe3(i,xx(1),y,nx,my,xx,yy,cc,f0)
 304          call jeppe3(i,xx(2),y,nx,my,xx,yy,cc,f1)
 305          if (f0.gt.0.d0.and.f1.gt.0.d0) then
 306             z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
 307          else
 308             z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
 309          end if
 310
 311       else if (n.gt.0.and.m.eq.my) then
 312 C--   Extrapolate to large q.
 313          call jeppe3(i,x,yy(my),nx,my,xx,yy,cc,f0)
 314          call jeppe3(i,x,yy(my-1),nx,my,xx,yy,cc,f1)
 315          if (f0.gt.0.d0.and.f1.gt.0.d0) then
 316             z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
 317      &           (y-yy(my)))
 318          else
 319             z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 320          end if
 321
 322       else if (n.eq.0.and.m.eq.my) then
 323 C--   Extrapolate to small x AND large q.
 324          call jeppe3(i,xx(1),yy(my),nx,my,xx,yy,cc,f0)
 325          call jeppe3(i,xx(1),yy(my-1),nx,my,xx,yy,cc,f1)
 326          if (f0.gt.0.d0.and.f1.gt.0.d0) then
 327             z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
 328      &           (y-yy(my)))
 329          else
 330             z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 331          end if
 332          call jeppe3(i,xx(2),yy(my),nx,my,xx,yy,cc,f0)
 333          call jeppe3(i,xx(2),yy(my-1),nx,my,xx,yy,cc,f1)
 334          if (f0.gt.0.d0.and.f1.gt.0.d0) then
 335             z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
 336      &           (y-yy(my)))
 337          else
 338             z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 339          end if
 340          if (z0.gt.0.d0.and.z1.gt.0.d0) then
 341             z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
 342          else
 343             z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
 344          end if
 345
 346       else
 347 C--   Set parton distribution to zero otherwise.
 348          z = 0.d0
 349
 350       end if
 351
 352       return
 353       end
 354
 355 C--   G.W. 02/07/2007 Copy of the original jeppe2,
 356 C--   only used for extrapolation.
 357       subroutine jeppe3(i,x,y,nx,my,xx,yy,cc,z)
 358       implicit real*8(a-h,o-z)
 359       parameter(nhess=30)
 360       dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)
 361       n=locx(xx,nx,x)
 362       m=locx(yy,my,y)
 363       t=(x-xx(n))/(xx(n+1)-xx(n))
 364       u=(y-yy(m))/(yy(m+1)-yy(m))
 365       z=0.d0
 366       do l=4,1,-1
 367          z=t*z+((cc(i,n,m,l,4)*u+cc(i,n,m,l,3))*u
 368      &        +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
 369       enddo
 370       return
 371       end
 372
 373       integer function locx(xx,nx,x)
 374       implicit real*8(a-h,o-z)
 375       dimension xx(nx)
 376 c$$$      if(x.le.xx(1)) then
 377       if(x.eq.xx(1)) then ! G.W. 02/07/2007
 378       locx=1
 379       return
 380       endif
 381 c$$$      if(x.ge.xx(nx)) then
 382       if(x.eq.xx(nx)) then ! G.W. 02/07/2007
 383       locx=nx-1
 384       return
 385       endif
 386       ju=nx+1
 387       jl=0
 388     1 if((ju-jl).le.1) go to 2
 389       jm=(ju+jl)/2
 390       if(x.ge.xx(jm)) then
 391       jl=jm
 392       else
 393       ju=jm
 394       endif
 395       go to 1
 396     2 locx=jl
 397       return
 398       end
 399
 400       real*8 function  polderiv(x1,x2,x3,y1,y2,y3)
 401       implicit real*8(a-h,o-z)
 402       polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
 403      .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
 404       return
 405       end