LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf-5.9.1 / src / wrapmrstqed.f
1 ! -*- F90 -*-
2
3
4       subroutine MRSTqedevolve(x,Q,pdf,photon)
5       implicit real*8(a-h,o-z)
6       include 'parmsetup.inc'
7       character*16 name(nmxset)
8       integer nmem(nmxset),ndef(nmxset),mmem
9 !      integer member(nmxset)
10       integer nset,iset
11       common/NAME/name,nmem,ndef,mmem
12       parameter(nx=49,nq=37,np=9,nqc0=2,nqb0=11,nqc=35,nqb=26,          &
13      &nhess=30)
14       double precision gridx(nmxgridx),gridq(nmxgridq)
15       integer ngridx,ngridq,jx,jq
16       real*8 pdf(-6:6),photon,phot
17       real*8 f1(nx,nq)                                                  &
18      &,f2(nx,nq)                                                        &
19      &,f3(nx,nq)                                                        &
20      &,f4(nx,nq)                                                        &
21      &,f5(nx,nq)                                                        &
22      &,f6(nx,nq)                                                        &
23      &,f7(nx,nq)                                                        &
24      &,f8(nx,nq)                                                        &
25      &,f9(nx,nq)                                                        &
26      &,fc(nx,nqc),fb(nx,nqb)
27       real*8 qq(nq),xx(nx),                                             &
28      &cc1(0:nhess,nx,nq,4,4,nmxset),cc2(0:nhess,nx,nq,4,4,nmxset),      &
29      &cc3(0:nhess,nx,nq,4,4,nmxset),cc4(0:nhess,nx,nq,4,4,nmxset),      &
30      &cc6(0:nhess,nx,nq,4,4,nmxset),cc8(0:nhess,nx,nq,4,4,nmxset),      &
31      &cc9(0:nhess,nx,nq,4,4,nmxset),                                    &
32      &ccc(0:nhess,nx,nqc,4,4,nmxset),ccb(0:nhess,nx,nqb,4,4,nmxset)
33       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
34       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,                                 &
35      &              1d-4,2d-4,4d-4,6d-4,8d-4,                           &
36      &              1d-3,2d-3,4d-3,6d-3,8d-3,                           &
37      &              1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,               &
38      &           .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,     &
39      &           .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,     &
40      &           .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,        &
41      &           .8d0,.9d0,1d0/
42       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,       &
43      &        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,                          &
44      &        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,          &
45      &        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,              &
46      &        1.8d6,3.2d6,5.6d6,1d7/
47       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
48       save
49
50       xsave=x
51       qsq = q*q
52       q2save=qsq
53
54       xlog=dlog(x)
55       qsqlog=dlog(qsq)
56
57       call getnset(iset)
58 !      imem=member(iset)
59       call getnmem(iset,imem)
60       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1(0,1,1,1,1,iset),upv)
61       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2(0,1,1,1,1,iset),dnv)
62       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3(0,1,1,1,1,iset),glu)
63       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4(0,1,1,1,1,iset),usea)
64       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6(0,1,1,1,1,iset),str)
65       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8(0,1,1,1,1,iset),dsea)
66       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc9(0,1,1,1,1,iset),photon)
67
68       chm=0.d0
69       if(qsq.gt.emc2) then
70       call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc(0,1,1,1,1,iset),chm)
71       endif
72
73       bot=0.d0
74       if(qsq.gt.emb2) then
75       call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb(0,1,1,1,1,iset),bot)
76       endif
77
78       pdf(0)  = glu
79       pdf(1)  = dnv+dsea
80       pdf(-1) = dsea
81       pdf(2)  = upv+usea
82       pdf(-2) = usea
83       pdf(3)  = str
84       pdf(-3) = str
85       pdf(4)  = chm
86       pdf(-4) = chm
87       pdf(5)  = bot
88       pdf(-5) = bot
89       pdf(6) = 0.0d0
90       pdf(-6) = 0.0d0
91
92       x=xsave
93       qsq=q2save
94       return
95
96       entry MRSTqedgetgrid(nset,ngridx,ngridq,gridx,gridq)
97       do jx=1,nx
98           gridx(jx)=xx(jx)
99       enddo
100       do jq=1,nq
101           gridq(jq)=qq(jq)
102       enddo
103       ngridx=nx
104       ngridq=nq        
105       return
106 !
107       entry MRSTqedread(nset)
108       read(1,*)nmem(nset),ndef(nset)
109 !      print *,nmem(nset),ndef(nset)
110 !      do nm = 0,nmem-1
111       do nm = 0,nmem(nset)
112         do 20 n=1,nx-1
113         do 20 m=1,nq
114         read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),                      &
115      &                  f5(n,m),f7(n,m),f6(n,m),f8(n,m),f9(n,m)
116 ! notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea +9 - pho
117    20 continue
118 !      write(*,*)'PDF set ',nm,' first element ',f1(1,1)
119       do 40 m=1,nq
120       f1(nx,m)=0.d0
121       f2(nx,m)=0.d0
122       f3(nx,m)=0.d0
123       f4(nx,m)=0.d0
124       f5(nx,m)=0.d0
125       f6(nx,m)=0.d0
126       f7(nx,m)=0.d0
127       f8(nx,m)=0.d0
128       f9(nx,m)=0.d0
129    40 continue
130
131       do n=1,nx
132       xxl(n)=dlog(xx(n))
133       enddo
134       do m=1,nq
135       qql(m)=dlog(qq(m))
136       enddo
137
138       call jeppe1(nm,nx,nq,xxl,qql,f1,cc1(0,1,1,1,1,nset))
139       call jeppe1(nm,nx,nq,xxl,qql,f2,cc2(0,1,1,1,1,nset))
140       call jeppe1(nm,nx,nq,xxl,qql,f3,cc3(0,1,1,1,1,nset))
141       call jeppe1(nm,nx,nq,xxl,qql,f4,cc4(0,1,1,1,1,nset))
142       call jeppe1(nm,nx,nq,xxl,qql,f6,cc6(0,1,1,1,1,nset))
143       call jeppe1(nm,nx,nq,xxl,qql,f8,cc8(0,1,1,1,1,nset))
144       call jeppe1(nm,nx,nq,xxl,qql,f9,cc9(0,1,1,1,1,nset))
145
146       emc2=2.045
147       emb2=18.5
148
149       do 44 m=1,nqc
150       qqlc(m)=qql(m+nqc0)
151       do 44 n=1,nx
152       fc(n,m)=f5(n,m+nqc0)
153    44 continue
154       qqlc(1)=dlog(emc2)
155       call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc(0,1,1,1,1,nset))
156
157       do 45 m=1,nqb
158       qqlb(m)=qql(m+nqb0)
159       do 45 n=1,nx
160       fb(n,m)=f7(n,m+nqb0)
161    45 continue
162       qqlb(1)=dlog(emb2)
163       call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb(0,1,1,1,1,nset))
164
165
166       enddo
167    50 format(9f10.5)
168       return
169 !
170       entry MRSTqedalfa(nflav,alfas,Qalfa)
171         call getnset(iset)
172         call alphamrs(nflav,alfas,Qalfa)
173       return
174 !
175       entry MRSTqedinit(Eorder,Q2fit)
176       return
177 !
178       entry MRSTqedpdf(mem)
179         call getnset(iset)
180         call setnmem(iset,mem)
181
182       return
183 !
184       END