LHAPDF/lhapdf5.5.1/src/wrapmrst2006-lite.f

   1 ! -*- F90 -*-
   2
   3
   4 subroutine MRST2006evolve(x,Q,pdf)
   5   implicit none
   6   ! implicit real*8 (a-h,o-z)
   7   include 'parmsetup.inc'
   8   character*16 name(nmxset)
   9   character*80 line
  10   character*512 setpath
  11   integer i,f,nhess,nx,nq,np,nqc0,nqb0,nqc,nqb,n,m,io
  12   double precision x,q,xmin,xmax,qsqmin,qsqmax,emc2,emb2,eps,dummy,qsq, &
  13        xlog,qsqlog,res,ExtrapLHAPDF,InterpLHAPDF
  14   double precision emc2minuseps,emb2minuseps
  15   double precision xsave,q2save
  16   double precision pdf(-6:6),alfas,Qalfa,Q2fit,eorder
  17   double precision upv,dnv,glu,usea,str,dsea,sbar,chm,bot
  18   integer iset,nflav,mem,nset,nmem(nmxset),ndef(nmxset),mmem,imem
  19   common/NAME/name,nmem,ndef,mmem
  20   parameter(nx=64,nq=48,np=9,nqc0=4,nqb0=14,nqc=nq-nqc0,nqb=nq-nqb0)
  21   parameter(xmin=1d-6,xmax=1d0,qsqmin=1d0,qsqmax=1d9,eps=1d-6)
  22   parameter(nhess=0,emc2=1.43d0**2,emb2=4.3d0**2) ! MRST 2006 NNLO
  23   parameter(emc2minuseps=emc2-eps,emb2minuseps=emb2-eps)
  24   double precision f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq), &
  25        f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb),f9(nx,nq)
  26   double precision qq(nq),xx(nx),cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4), &
  27        cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),cc6(0:nhess,nx,nq,4,4), &
  28        cc8(0:nhess,nx,nq,4,4),cc9(0:nhess,nx,nq,4,4),ccc(0:nhess,nx,nqc,4,4), &
  29        ccb(0:nhess,nx,nqb,4,4)
  30   double precision xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
  31   data xx/1d-6,2d-6,4d-6,6d-6,8d-6,1d-5,2d-5,4d-5,6d-5,8d-5,1d-4,2d-4,4d-4,6d-4, &
  32        8d-4,1d-3,2d-3,4d-3,6d-3,8d-3,1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2, &
  33        .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,.3d0,.325d0,.35d0,.375d0, &
  34        .4d0,.425d0,.45d0,.475d0,.5d0,.525d0,.55d0,.575d0,.6d0,.625d0,.65d0,.675d0, &
  35        .7d0,.725d0,.75d0,.775d0,.8d0,.825d0,.85d0,.875d0,.9d0,.925d0,.95d0,.975d0,1d0/
  36   data qq/1.d0,1.25d0,1.5d0,emc2minuseps,emc2,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,1.2d1, &
  37        emb2minuseps,emb2,2.6d1,4d1,6.4d1,1d2,1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3, &
  38        1d4,1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,1.8d6,3.2d6,5.6d6,1d7,1.8d7,3.2d7, &
  39        5.6d7,1d8,1.8d8,3.2d8,5.6d8,1d9/
  40   save
  41
  42   xsave=x
  43   qsq=q*q
  44   q2save=qsq
  45
  46   if (qsq.gt.qq(nqc0).and.qsq.lt.qq(nqc0+1)) qsq = qq(nqc0)
  47   if (qsq.gt.qq(nqb0).and.qsq.lt.qq(nqb0+1)) qsq = qq(nqb0)
  48
  49   xlog=log10(x)
  50   qsqlog=log10(qsq)
  51
  52   call getnset(iset)
  53   call getnmem(iset,imem)
  54
  55   upv =  0.d0
  56   dnv =  0.d0
  57   glu =  0.d0
  58   usea = 0.d0
  59   str =  0.d0
  60   dsea = 0.d0
  61   sbar = 0.d0
  62   chm =  0.d0
  63   bot =  0.d0
  64
  65   if (x.le.0.d0.or.x.gt.xmax.or.qsq.lt.qsqmin) then
  66      print *,"Error in GetOnePDF: x,qsq = ",x,qsq
  67   else if (x.lt.xmin.or.qsq.gt.qsqmax) then ! extrapolate
  68      print *, "Warning in GetOnePDF, extrapolating: f = ",f,", x = ",x,", q = ",q
  69      upv =  ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc1)
  70      dnv =  ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc2)
  71      glu =  ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc3)
  72      usea = ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc4)
  73      str =  ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc6)
  74      dsea = ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc8)
  75      sbar = ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc9)
  76      if (qsq.ge.emc2) then ! chm
  77         chm = ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc)
  78      endif
  79      if (qsq.ge.emb2) then ! bot
  80         bot = ExtrapLHAPDF(imem,nhess,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb)
  81      endif
  82   else
  83      upv =  InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc1)
  84      dnv =  InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc2)
  85      glu =  InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc3)
  86      usea = InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc4)
  87      str =  InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc6)
  88      dsea = InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc8)
  89      sbar = InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nq,xxl,qql,cc9)
  90      if (qsq.ge.emc2) then ! chm
  91         chm = InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc)
  92      endif
  93      if (qsq.ge.emb2) then ! bot
  94         bot = InterpLHAPDF(imem,nhess,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb)
  95      endif
  96
  97   end if
  98
  99   pdf(0)  = glu
 100   pdf(1)  = dnv+dsea
 101   pdf(-1) = dsea
 102   pdf(2)  = upv+usea
 103   pdf(-2) = usea
 104   pdf(3)  = str
 105   pdf(-3) = sbar
 106   pdf(4)  = chm
 107   pdf(-4) = chm
 108   pdf(5)  = bot
 109   pdf(-5) = bot
 110   pdf(6)  = 0.0d0
 111   pdf(-6) = 0.0d0
 112
 113   x=xsave
 114   qsq=q2save
 115   return
 116
 117   !----------------------------------------------------------------------
 118
 119   entry MRST2006read(nset)
 120   ! - dummy read in to get to End: (stream 1 is still open)
 121   read(1,*)nmem(nset),ndef(nset)
 122
 123   do i=0,nmem(nset)
 124      do n=1,nx-1
 125         do m=1,nq
 126            read(1,'(a)'),line
 127         enddo
 128      enddo
 129   enddo
 130
 131
 132   return
 133   !----------------------------------------------------------------------
 134
 135   entry MRST2006alfa(nflav,alfas,Qalfa)
 136   ! print *,'MRST2006alfa',nflav,alfas,Qalfa
 137   call getnset(iset)
 138   call alphamrs(nflav,alfas,Qalfa)
 139   ! print *,'MRST2006alfa',nflav,alfas,Qalfa
 140   return
 141
 142   !----------------------------------------------------------------------
 143
 144   entry MRST2006init(Eorder,Q2fit)
 145   return
 146
 147   !----------------------------------------------------------------------
 148
 149   entry MRST2006pdf(mem)
 150   call getnset(iset)
 151   call setnmem(iset,mem)
 152
 153   ! have to reopen stream 1
 154   call getsetpath(setpath)
 155   open(1,file=setpath(1:len_trim(setpath)))
 156
 157   line = ''
 158   do while (line(2:11).ne.'Evolution:')
 159      read(1,'(a)'),line
 160   enddo
 161   read(1,'(a)'),line
 162   read(1,'(a)'),line
 163
 164   read(1,*)nmem(iset),ndef(iset)
 165
 166   ! - dummy read up to the member requested
 167   do i=0,mem-1
 168      do n=1,nx-1
 169         do m=1,nq
 170            read(1,'(a)')line
 171         enddo
 172      enddo
 173   enddo
 174
 175
 176   !- read in the data of the requested member
 177   do  n=1,nx-1
 178      do  m=1,nq
 179         ! print *,i,n,m
 180         read(1,50) f1(n,m),f2(n,m),f3(n,m),f4(n,m),f5(n,m),f7(n,m),f6(n,m),f8(n,m),f9(n,m)
 181         ! write(6,50) f1(n,m),f2(n,m),f3(n,m),f4(n,m),f5(n,m),f7(n,m),f6(n,m),f8(n,m),f9(n,m)
 182         !--   Notation:1=upv 2=dnv 3=glu 4=usea 5=chm 6=str 7=bot 8=dsea 9=sbar
 183      enddo
 184   enddo
 185   ! write(*,*)'PDF set ',i,' first element ',f1(1,1)
 186   do m=1,nq
 187      f1(nx,m)=0d0
 188      f2(nx,m)=0d0
 189      f3(nx,m)=0d0
 190      f4(nx,m)=0d0
 191      f5(nx,m)=0d0
 192      f6(nx,m)=0d0
 193      f7(nx,m)=0d0
 194      f8(nx,m)=0d0
 195      f9(nx,m)=0d0
 196   enddo
 197
 198   do n=1,nx
 199      xxl(n)=log10(xx(n))
 200   enddo
 201   do m=1,nq
 202      qql(m)=log10(qq(m))
 203   enddo
 204
 205   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f1,cc1)
 206   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f2,cc2)
 207   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f3,cc3)
 208   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f4,cc4)
 209   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f6,cc6)
 210   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f8,cc8)
 211   call InitialLHAPDF(i,nhess,nx,nq,nqc0,nqb0,xxl,qql,f9,cc9)
 212
 213   do m=1,nqc
 214      qqlc(m)=qql(m+nqc0)
 215      do  n=1,nx
 216         fc(n,m)=f5(n,m+nqc0)
 217      enddo
 218   enddo
 219   call InitialLHAPDF(i,nhess,nx,nqc,nqc0-nqc0,nqb0-nqc0,xxl,qqlc,fc,ccc)
 220
 221   do m=1,nqb
 222      qqlb(m)=qql(m+nqb0)
 223      do n=1,nx
 224         fb(n,m)=f7(n,m+nqb0)
 225      enddo
 226   enddo
 227   call InitialLHAPDF(i,nhess,nx,nqb,nqc0-nqb0,nqb0-nqb0,xxl,qqlb,fb,ccb)
 228
 229   close(1)
 230
 231 50 format(9e13.5)
 232
 233 end subroutine MRST2006evolve
 234
 235
 236 !----------------------------------------------------------------------
 237
 238
 239 subroutine InitialLHAPDF(i,nhess,nx,my,myc0,myb0,xx,yy,ff,cc)
 240   implicit none
 241   integer nhess,i,nx,my,myc0,myb0,j,k,l,m,n
 242   double precision xx(nx),yy(my),ff(nx,my),ff1(nx,my),ff2(nx,my),ff12(nx,my), &
 243        yy0(4),yy1(4),yy2(4),yy12(4),z(16),cl(16),cc(0:nhess,nx,my,4,4), &
 244        iwt(16,16),dx,dy,polderiv,d1,d2,d1d2,xxd
 245   data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, &
 246        0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0, &
 247        -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0, &
 248        2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0, &
 249        0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0, &
 250        0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0, &
 251        0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1, &
 252        0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1, &
 253        -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0, &
 254        0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0, &
 255        9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2, &
 256        -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2, &
 257        2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0, &
 258        0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0, &
 259        -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1, &
 260        4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
 261
 262   do m=1,my
 263      dx=xx(2)-xx(1)
 264      ff1(1,m)=(ff(2,m)-ff(1,m))/dx
 265      dx=xx(nx)-xx(nx-1)
 266      ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
 267      do n=2,nx-1
 268         ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),ff(n+1,m))
 269      enddo
 270   enddo
 271
 272   !--   Calculate the derivatives at qsq=emc2-eps,emc2,emb2-eps,emb2
 273   !--   in a similar way as at the endpoints qsqmin and qsqmax.
 274   do n=1,nx
 275      do m = 1, my
 276         if (m.eq.1.or.m.eq.myc0+1.or.m.eq.myb0+1) then
 277            dy=yy(m+1)-yy(m)
 278            ff2(n,m)=(ff(n,m+1)-ff(n,m))/dy
 279         else if (m.eq.my.or.m.eq.myc0.or.m.eq.myb0) then
 280            dy=yy(m)-yy(m-1)
 281            ff2(n,m)=(ff(n,m)-ff(n,m-1))/dy
 282         else
 283            ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),ff(n,m+1))
 284         end if
 285      end do
 286   end do
 287
 288   do m=1,my
 289      dx=xx(2)-xx(1)
 290      ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
 291      dx=xx(nx)-xx(nx-1)
 292      ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
 293      do n=2,nx-1
 294         ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),ff2(n+1,m))
 295      enddo
 296   enddo
 297
 298   do n=1,nx-1
 299      do m=1,my-1
 300         d1=xx(n+1)-xx(n)
 301         d2=yy(m+1)-yy(m)
 302         d1d2=d1*d2
 303
 304         yy0(1)=ff(n,m)
 305         yy0(2)=ff(n+1,m)
 306         yy0(3)=ff(n+1,m+1)
 307         yy0(4)=ff(n,m+1)
 308
 309         yy1(1)=ff1(n,m)
 310         yy1(2)=ff1(n+1,m)
 311         yy1(3)=ff1(n+1,m+1)
 312         yy1(4)=ff1(n,m+1)
 313
 314         yy2(1)=ff2(n,m)
 315         yy2(2)=ff2(n+1,m)
 316         yy2(3)=ff2(n+1,m+1)
 317         yy2(4)=ff2(n,m+1)
 318
 319         yy12(1)=ff12(n,m)
 320         yy12(2)=ff12(n+1,m)
 321         yy12(3)=ff12(n+1,m+1)
 322         yy12(4)=ff12(n,m+1)
 323
 324         do k=1,4
 325            z(k)=yy0(k)
 326            z(k+4)=yy1(k)*d1
 327            z(k+8)=yy2(k)*d2
 328            z(k+12)=yy12(k)*d1d2
 329         enddo
 330
 331         do l=1,16
 332            xxd=0.d0
 333            do k=1,16
 334               xxd=xxd+iwt(k,l)*z(k)
 335            enddo
 336            cl(l)=xxd
 337         enddo
 338         l=0
 339         do k=1,4
 340            do j=1,4
 341               l=l+1
 342               cc(i,n,m,k,j)=cl(l)
 343            enddo
 344         enddo
 345      enddo
 346   enddo
 347   return
 348 end subroutine InitialLHAPDF
 349
 350
 351 !----------------------------------------------------------------------
 352
 353
 354 double precision function InterpLHAPDF(i,nhess,x,y,nx,my,xx,yy,cc)
 355   implicit none
 356   integer i,nx,my,nhess,locx,l,m,n
 357   double precision xx(nx),yy(my),cc(0:nhess,nx,my,4,4),x,y,z,t,u
 358
 359   n=locx(xx,nx,x)
 360   m=locx(yy,my,y)
 361
 362   t=(x-xx(n))/(xx(n+1)-xx(n))
 363   u=(y-yy(m))/(yy(m+1)-yy(m))
 364
 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 + cc(i,n,m,l,2))*u + cc(i,n,m,l,1)
 368   enddo
 369
 370   InterpLHAPDF = z
 371
 372   return
 373 end function InterpLHAPDF
 374
 375
 376 !----------------------------------------------------------------------
 377
 378
 379 double precision function ExtrapLHAPDF(i,nhess,x,y,nx,my,xx,yy,cc)
 380   implicit none
 381   integer i,nx,my,nhess,locx,n,m
 382   double precision xx(nx),yy(my),cc(0:nhess,nx,my,4,4),x,y,z,f0,f1,z0,z1,InterpLHAPDF
 383
 384   n=locx(xx,nx,x)           ! 0: below xmin, nx: above xmax
 385   m=locx(yy,my,y)           ! 0: below qsqmin, my: above qsqmax
 386
 387   !--   If extrapolation in small x only:
 388   if (n.eq.0.and.m.gt.0.and.m.lt.my) then
 389      f0 = InterpLHAPDF(i,nhess,xx(1),y,nx,my,xx,yy,cc)
 390      f1 = InterpLHAPDF(i,nhess,xx(2),y,nx,my,xx,yy,cc)
 391      if (f0.gt.0.d0.and.f1.gt.0.d0) then
 392         z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
 393      else
 394         z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
 395      end if
 396      !--   If extrapolation into large q only:
 397   else if (n.gt.0.and.m.eq.my) then
 398      f0 = InterpLHAPDF(i,nhess,x,yy(my),nx,my,xx,yy,cc)
 399      f1 = InterpLHAPDF(i,nhess,x,yy(my-1),nx,my,xx,yy,cc)
 400      if (f0.gt.0.d0.and.f1.gt.0.d0) then
 401         z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*(y-yy(my)))
 402      else
 403         z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 404      end if
 405      !--   If extrapolation into large q AND small x:
 406   else if (n.eq.0.and.m.eq.my) then
 407      f0 = InterpLHAPDF(i,nhess,xx(1),yy(my),nx,my,xx,yy,cc)
 408      f1 = InterpLHAPDF(i,nhess,xx(1),yy(my-1),nx,my,xx,yy,cc)
 409      if (f0.gt.0.d0.and.f1.gt.0.d0) then
 410         z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*(y-yy(my)))
 411      else
 412         z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 413      end if
 414      f0 = InterpLHAPDF(i,nhess,xx(2),yy(my),nx,my,xx,yy,cc)
 415      f1 = InterpLHAPDF(i,nhess,xx(2),yy(my-1),nx,my,xx,yy,cc)
 416      if (f0.gt.0.d0.and.f1.gt.0.d0) then
 417         z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*(y-yy(my)))
 418      else
 419         z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
 420      end if
 421      if (z0.gt.0.d0.and.z1.gt.0.d0) then
 422         z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
 423      else
 424         z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
 425      end if
 426   else
 427      print *,"Error in ExtrapLHAPDF"
 428      stop
 429   end if
 430
 431   ExtrapLHAPDF = z
 432
 433   return
 434 end function ExtrapLHAPDF
 435
 436
 437 !----------------------------------------------------------------------