Change needed for G4
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapa02.f
1       subroutine A02evolve(xb,Q,fset)
2       implicit none
3       include 'parmsetup.inc'
4       integer npdf,npar,kschem,i,k,n,m,nxbb
5       integer nxb,nq,np,nvar,pdfmem,imem
6       parameter(nxb=99,nq=20,np=9,nvar=15)
7       character*16 name(nmxset)
8       integer nmem(nmxset),ndef(nmxset),mem
9       common/NAME/name,nmem,ndef,mem
10       integer nset
11       real*4 f(0:nvar,nxb,nq+1,0:np)
12       real*8 x1,xd,del,dels,delx,delx1,xlog1
13       real*8 pdfs(np),fset(-6:6),alfas
14       real*8 x,Q,xb,q2,xmin,xmax,qsq,qsqmin,qsqmax,a,b,ss
15       data xmin,xmax,qsqmin,qsqmax/1d-7,1d0,0.8d0,2d8/
16
17       save f,npdf,npar,pdfmem,dels,delx,x1,delx1,xlog1,nxbb
18
19       q2=Q*Q
20       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
21       if(xb.lt.xmin.or.xb.gt.xmax)       print 98,x
22   99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ')
23   98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ')
24
25       x=max(xb,xmin)
26       x=min(x,xmax)
27       qsq=max(q2,qsqmin)
28       qsq=min(qsq,qsqmax)
29
30       if (x.gt.x1) then
31         xd=(1-x1)**2-(1-x)**2
32         n=int(xd/delx1)+nxbb
33         a=xd/delx1-n+nxbb
34       else
35         xd=log(x)-xlog1
36         n=nxbb+int(xd/DELX)-1
37         a=xd/delx-n+nxbb
38       end if
39
40       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
41       m=int(ss/dels)+1
42       b=ss/dels-m+1
43
44       k=pdfmem
45
46       do i=1,npdf
47         if (n.gt.1.and.m.gt.1.and.n.ne.49) then 
48         pdfs(i)= f(k,n,m,i)*(1+a*b-a**2-b**2) 
49      +         + f(k,n+1,m+1,i)*a*b
50      +         + f(k,n+1,m,i)*a*(a-2*b+1)/2. 
51      +         + f(k,n,m+1,i)*b*(b-2*a+1)/2.
52      +         +  f(k,n-1,m,i)*a*(a-1)/2. 
53      +         + f(k,n,m-1,i)*b*(b-1)/2.
54         else
55         pdfs(i)= (1-a)*(1-b)*f(k,n,m,i) 
56      .         + (1-a)*b*f(k,n,m+1,i)
57      .         + a*(1-b)*f(k,n+1,m,i) 
58      .         + a*b*f(k,n+1,m+1,i)
59         end if
60       end do
61
62
63       fset(-6)=pdfs(9)
64       fset(-5)=pdfs(8)
65       fset(-4)=pdfs(7)
66       fset(-3)=pdfs(5)
67 c--reversed mrs 7/7/04 due
68       fset(-1)=pdfs(6)
69       fset(-2)=pdfs(4)
70       fset(0)=pdfs(3)
71 c--reversed mrs 7/7/04 due
72       fset(2)=pdfs(1)+pdfs(4)
73       fset(1)=pdfs(2)+pdfs(6)
74       fset(3)=pdfs(5)
75       fset(4)=pdfs(7)
76       fset(5)=pdfs(8)
77       fset(6)=pdfs(9)
78       return
79 *
80       entry A02alfa(alfas,Q)
81       q2=Q*Q
82       if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
83       qsq=max(q2,qsqmin)
84       qsq=min(qsq,qsqmax)
85
86       ss=log(log(qsq/0.04))-log(log(qsqmin/0.04))
87       m=int(ss/dels)+1
88       b=ss/dels-m+1
89
90       k=pdfmem
91       alfas=(1d0-b)*f(k,1,m,0)+b*f(k,1,m+1,0)
92       return
93 *
94       entry A02read(nset)
95 c following fix because members are 0-nvar
96 c      nmem = nvar + 1
97       nmem(nset) = nvar
98       ndef(nset) = 0
99       read(1,*) npdf
100       npar=nvar
101       do k=0,npar
102          do n=1,nxb-1
103             do m=1,nq
104                read(1,100) (f(k,n,m,i),i=0,npdf)
105             enddo
106          enddo
107       enddo
108   100 format (13f11.5)
109       return
110 *
111       entry A02init
112
113       dels=(log(log(qsqmax/0.04))-log(log(qsqmin/0.04)))/(nq-1)
114
115       nxbb=nxb/2
116       x1=0.3
117       xlog1=log(x1)
118       delx=(log(x1)-log(xmin))/(nxbb-1)
119       DELX1=(1-x1)**2/(nxbb+1)
120
121       do i=1,npdf
122          do m=1,nq
123             do k=1,npar
124                f(k,nxb,m,i)=0d0
125             end do
126          end do
127          do m=1,nq
128             do k=1,npar
129                f(k,nxb,m,0)=f(k,nxb-1,m,0)
130             end do
131          end do
132       end do
133       return
134 *
135       entry A02pdf(imem)
136       pdfmem=imem
137       if ((pdfmem.lt.0).or.(pdfmem.gt.npar)) then
138          write(*,*) 'A02 PDF set:'
139          write(*,*) 'PDF member out of range:'
140          write(*,*) 'member = ',pdfmem,'    member range = (0,',npar,')'
141          stop
142       endif
143       return
144       end