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