LHAPDF/lhapdf5.5.1/src/wrapEVLCTEQ.f

   1 ! -*- F90 -*-
   2
   3
   4       subroutine EVLCTEQevolve(x,Q,f)
   5       implicit none
   6       include 'parmsetup.inc'
   7       character*16 name(nmxset)
   8       integer nmem(nmxset),ndef(nmxset),mmem
   9       common/NAME/name,nmem,ndef,mmem
  10       real*8 f(-6:6)
  11       real*8 x,Q,xx,qq,up,dn,ubar,dbar,gluon,sbar,cbar,bbar,tbar,Blam
  12       real*8 xmin,qmax,anx,anq,Q2fit,qini,Ahdn,ordi,anf,mass,pi,alfas
  13       real*8 Q0,alfas0
  14       real*8 CtLhpardis,CtLhALPI,parmflavor,CtLhastolam
  15       integer nflavor,Eorder,j,nx,nq
  16       integer nset
  17       data pi / 3.141592653589793d0 /
  18       external parmflavor
  19       save nflavor,xmin,qmax,anx,anq
  20
  21 !
  22         xx = x
  23         qq = Q
  24         up   = xx*CtLhpardis(1,xx,qq)
  25         dn   = xx*CtLhpardis(2,xx,qq)
  26         ubar = xx*CtLhpardis(-1,xx,qq)
  27         dbar = xx*CtLhpardis(-2,xx,qq)
  28         gluon= xx*CtLhpardis( 0,xx,qq)
  29
  30         if(nflavor .ge. 3) then
  31            sbar = xx*CtLhpardis(-3,xx,qq)
  32         else
  33            sbar = 0.d0
  34         endif
  35
  36         if(nflavor .ge. 4) then
  37            cbar = xx*CtLhpardis(-4,xx,qq)
  38         else
  39            cbar = 0.d0
  40         endif
  41
  42         if(nflavor .ge. 5) then
  43            bbar = xx*CtLhpardis(-5,xx,qq)
  44         else
  45            bbar = 0.d0
  46         endif
  47
  48         if(nflavor .eq. 6) then
  49            tbar = xx*CtLhpardis(-6,xx,qq)
  50         else
  51            tbar = 0.d0
  52         endif
  53
  54         f(6) = tbar
  55         f(5) = bbar
  56         f(4) = cbar
  57         f(3) = sbar
  58         f(2) = up
  59         f(1) = dn
  60         f(0) = gluon
  61         f(-1) = dbar
  62         f(-2) = ubar
  63         f(-3) = sbar
  64         f(-4) = cbar
  65         f(-5) = bbar
  66         f(-6) = tbar
  67 !
  68         return
  69 !
  70       entry EVLCTEQread(nset)
  71       call CtLhbldat1
  72       call CtLhbldat2
  73       read(1,*) xmin,qmax,nx,nq
  74       anx=nx
  75       anq=nq-1
  76       return
  77 !
  78       entry EVLCTEQalfa(alfas,Q)
  79       alfas = pi*CtLhALPI(Q)
  80       return
  81 !
  82       entry EVLCTEQinit(nset,Eorder,Q2fit)
  83 !
  84       Ahdn=1d0
  85       call CtLhParPdf(1,'IHDN',Ahdn,j)
  86       call CtLhParPdf(1,'QMAX',qmax,j)
  87       call CtLhParPdf(1,'XMIN',xmin,j)
  88       call GetOrderAsM(nset,j)
  89       ordi=j+1.0
  90       call CtLhParQcd(1,'ORDR',ordi,j)
  91       call GetThresholdM(nset,4,mass)
  92       call CtLhParQcd(1,'M4',mass,j)
  93       call GetThresholdM(nset,5,mass)
  94       call CtLhParQcd(1,'M5',mass,j)
  95       mass=180d0
  96       call CtLhParQcd(1,'M6',mass,j)
  97       call GetNfM(nset,j)
  98       nflavor=j
  99       Anf=nflavor
 100       call CtLhParQcd(1,'NFL',Anf,j)
 101       qini=sqrt(Q2fit)
 102       call CtLhParPdf(1,'QINI',qini,j)
 103       call CtLhEvlpar(1,'NFMX', Anf, j)
 104       ordi=Eorder+1.0
 105       call CtLhParPdf(1,'IKNL',ordi,j)
 106       return
 107 !
 108       entry EVLCTEQpdf(nset)
 109       call GetOrderAsM(nset,j)
 110       call GetAlfas(nset,alfas0,Q0)
 111       Blam=CtLhastolam(alfas0,Q0,j+1,nflavor)
 112       call CtLhSetLam (nflavor,Blam,j+1)
 113       call CtLhParPdf(1,'NX',anx,j)
 114       call CtLhParPdf(1,'NT',anq,j)
 115       call CtLhEvolve (parmflavor,j)
 116       if (j .ne. 0) then
 117          write(*,*) 'EVLCTEQ Evolve Error code :',j
 118          stop
 119       endif
 120       return
 121 !
 122       END
 123 !
 124       double precision function parmflavor(i,x)
 125       implicit none
 126       real*8 x,f(-6:6)
 127       integer i,i0
 128       integer iset
 129 !
 130       call getnset(iset)
 131       call parmPDF(iset,x,f)
 132       i0=i
 133       if (i.eq.-2) i0=-1
 134       if (i.eq.-1) i0=-2
 135       if (i.eq.1) i0=2
 136       if (i.eq.2) i0=1
 137       parmflavor=f(i0)/x
 138       return
 139 !
 140       END
 141
 142       function CtLhastolam(as,q,nloop,nf)
 143       implicit double precision (a-h,o-z)
 144       data pi / 3.141592653589793d0 /
 145       xlp = nloop-1d0
 146       b  = (33.0-2.0*nf)/pi/12.0
 147       bp = (153.0 - 19.0*nf) / pi / 2.0 / (33.0 - 2.0*nf) * xlp
 148       t  = 1.0/b/as
 149 !-----------------------------------------------------------
 150 ! Solve the equation
 151 !
 152     1 xlt = log(t)
 153       ot = t
 154 !-----------------------------------------------------------
 155 ! Solve the equation
 156 ! Value and Derivative of alfa with respect to t
 157 !
 158       as0  = 1/b/t - bp*xlt/(b*t)**2
 159       as1  = - 1/b/t**2 -bp/b**2*(1-2*xlt)/t**3
 160       t  = (as-as0)/as1 + t
 161       if(abs(ot-t)/ot.gt..00001)goto 1
 162       CtLhastolam = q/exp(t/2.0)
 163       return
 164       END