Version update.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapmrst.f
1       subroutine MRSTevolve(x,Q,pdf)
2       implicit real*8(a-h,o-z)
3       include 'parmsetup.inc'
4       character*16 name(nmxset)
5       integer nmem(nmxset),ndef(nmxset),mmem
6 c      integer member(nmxset)
7       integer nset,iset
8       common/NAME/name,nmem,ndef,mmem
9       parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26,
10      .nhess=30)
11       real*8 pdf(-6:6)
12       real*8 f1(nx,nq)
13      .,f2(nx,nq)
14      .,f3(nx,nq)
15      .,f4(nx,nq)
16      .,f5(nx,nq)
17      .,f6(nx,nq)
18      .,f7(nx,nq)
19      .,f8(nx,nq)
20      .,fc(nx,nqc),fb(nx,nqb)
21       real*8 qq(nq),xx(nx),
22      .cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4),
23      .cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),
24      .cc6(0:nhess,nx,nq,4,4),cc8(0:nhess,nx,nq,4,4),
25      .ccc(0:nhess,nx,nqc,4,4),ccb(0:nhess,nx,nqb,4,4)
26       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
27       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
28      .        1d-4,2d-4,4d-4,6d-4,8d-4,
29      .        1d-3,2d-3,4d-3,6d-3,8d-3,
30      .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
31      .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
32      .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
33      .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
34      .     .8d0,.9d0,1d0/
35       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
36      .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
37      .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
38      .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
39      .        1.8d6,3.2d6,5.6d6,1d7/
40       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
41       save
42       
43       xsave=x
44       qsq = q*q
45       q2save=qsq
46
47       xlog=dlog(x)
48       qsqlog=dlog(qsq)
49
50       call getnset(iset)
51 c      imem=member(iset)
52       call getnmem(iset,imem)
53       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
54       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
55       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
56       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
57       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
58       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
59
60       chm=0.d0
61       if(qsq.gt.emc2) then 
62       call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
63       endif
64
65       bot=0.d0
66       if(qsq.gt.emb2) then 
67       call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
68       endif
69
70       pdf(0)  = glu
71       pdf(1)  = dnv+dsea
72       pdf(-1) = dsea
73       pdf(2)  = upv+usea
74       pdf(-2) = usea
75       pdf(3)  = str
76       pdf(-3) = str
77       pdf(4)  = chm
78       pdf(-4) = chm
79       pdf(5)  = bot
80       pdf(-5) = bot
81       pdf(6)  = 0.0d0
82       pdf(-6) = 0.0d0
83       
84       x=xsave
85       qsq=q2save
86       return
87 *
88       entry MRSTread(nset)
89       read(1,*)nmem(nset),ndef(nset)
90 c      print *,nmem(nset),ndef(nset)
91 c      do nm = 0,nmem-1
92       do nm = 0,nmem(nset)
93         do 20 n=1,nx-1
94         do 20 m=1,nq
95         read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
96      .            f5(n,m),f7(n,m),f6(n,m),f8(n,m)
97 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
98   20  continue
99 c      write(*,*)'PDF set ',nm,' first element ',f1(1,1)
100       do 40 m=1,nq
101       f1(nx,m)=0.d0
102       f2(nx,m)=0.d0
103       f3(nx,m)=0.d0
104       f4(nx,m)=0.d0
105       f5(nx,m)=0.d0
106       f6(nx,m)=0.d0
107       f7(nx,m)=0.d0
108       f8(nx,m)=0.d0
109   40  continue
110
111       do n=1,nx
112       xxl(n)=dlog(xx(n))
113       enddo
114       do m=1,nq
115       qql(m)=dlog(qq(m))
116       enddo
117
118       call jeppe1(nm,nx,nq,xxl,qql,f1,cc1)
119       call jeppe1(nm,nx,nq,xxl,qql,f2,cc2)
120       call jeppe1(nm,nx,nq,xxl,qql,f3,cc3)
121       call jeppe1(nm,nx,nq,xxl,qql,f4,cc4)
122       call jeppe1(nm,nx,nq,xxl,qql,f6,cc6)
123       call jeppe1(nm,nx,nq,xxl,qql,f8,cc8)
124
125       emc2=2.045
126       emb2=18.5
127
128       do 44 m=1,nqc
129       qqlc(m)=qql(m+nqc0)
130       do 44 n=1,nx
131       fc(n,m)=f5(n,m+nqc0)
132    44 continue
133       qqlc(1)=dlog(emc2)
134       call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc)
135
136       do 45 m=1,nqb
137       qqlb(m)=qql(m+nqb0)
138       do 45 n=1,nx
139       fb(n,m)=f7(n,m+nqb0)
140    45 continue
141       qqlb(1)=dlog(emb2)
142       call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb)
143
144
145       enddo
146   50  format(8f10.5)
147       return
148 *    
149
150       entry MRSTalfa(nflav,alfas,Qalfa)
151         call getnset(iset)
152 c        mem = member(iset)
153 c       call setnmem(member(iset))
154         call alphamrs(nflav,alfas,Qalfa)
155       return
156 *
157       entry MRSTinit(Eorder,Q2fit)
158       return
159 *
160       entry MRSTpdf(mem)
161 c      if(mem.eq.0) mem=ndef
162 c      imem = mem
163         call getnset(iset)
164 c       member(iset)=mem
165         call setnmem(iset,mem)
166
167       return
168 *
169       end
170
171       subroutine jeppe1(imem,nx,my,xx,yy,ff,cc)
172       implicit real*8(a-h,o-z)
173       parameter(nnx=49,mmy=37,nhess=30)
174       dimension xx(nx),yy(my),ff(nnx,mmy),
175      xff1(nnx,mmy),ff2(nnx,mmy),
176      xff12(nnx,mmy),yy0(4),yy1(4),yy2(4),yy12(4),z(16),wt(16,16),
177      xcl(16),cc(0:nhess,nx,my,4,4),iwt(16,16)
178
179       data iwt/1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
180      x            0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
181      x            -3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0,
182      x            2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0,
183      x            0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
184      x            0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
185      x            0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1,
186      x            0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1,
187      x            -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,
188      x            0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0,
189      x            9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2,
190      x            -6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2,
191      x            2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
192      x            0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0,
193      x            -6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1,
194      x            4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1/
195
196
197       do 42 m=1,my
198       dx=xx(2)-xx(1)
199       ff1(1,m)=(ff(2,m)-ff(1,m))/dx
200       dx=xx(nx)-xx(nx-1)
201       ff1(nx,m)=(ff(nx,m)-ff(nx-1,m))/dx
202       do 41 n=2,nx-1
203       ff1(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff(n-1,m),ff(n,m),
204      xff(n+1,m))
205    41 continue
206    42 continue
207
208       do 44 n=1,nx
209       dy=yy(2)-yy(1)
210       ff2(n,1)=(ff(n,2)-ff(n,1))/dy
211       dy=yy(my)-yy(my-1)
212       ff2(n,my)=(ff(n,my)-ff(n,my-1))/dy
213       do 43 m=2,my-1
214       ff2(n,m)=polderiv(yy(m-1),yy(m),yy(m+1),ff(n,m-1),ff(n,m),
215      xff(n,m+1))
216    43 continue
217    44 continue
218
219       do 46 m=1,my
220       dx=xx(2)-xx(1)
221       ff12(1,m)=(ff2(2,m)-ff2(1,m))/dx
222       dx=xx(nx)-xx(nx-1)
223       ff12(nx,m)=(ff2(nx,m)-ff2(nx-1,m))/dx
224       do 45 n=2,nx-1
225       ff12(n,m)=polderiv(xx(n-1),xx(n),xx(n+1),ff2(n-1,m),ff2(n,m),
226      xff2(n+1,m))
227    45 continue
228    46 continue
229
230       do 53 n=1,nx-1
231       do 52 m=1,my-1
232       d1=xx(n+1)-xx(n)
233       d2=yy(m+1)-yy(m)
234       d1d2=d1*d2
235
236       yy0(1)=ff(n,m)
237       yy0(2)=ff(n+1,m)
238       yy0(3)=ff(n+1,m+1)
239       yy0(4)=ff(n,m+1)
240
241       yy1(1)=ff1(n,m)
242       yy1(2)=ff1(n+1,m)
243       yy1(3)=ff1(n+1,m+1)
244       yy1(4)=ff1(n,m+1)
245
246       yy2(1)=ff2(n,m)
247       yy2(2)=ff2(n+1,m)
248       yy2(3)=ff2(n+1,m+1)
249       yy2(4)=ff2(n,m+1)
250
251       yy12(1)=ff12(n,m)
252       yy12(2)=ff12(n+1,m)
253       yy12(3)=ff12(n+1,m+1)
254       yy12(4)=ff12(n,m+1)
255
256       do 47 k=1,4
257       z(k)=yy0(k)
258       z(k+4)=yy1(k)*d1
259       z(k+8)=yy2(k)*d2
260       z(k+12)=yy12(k)*d1d2
261    47 continue
262
263       do 49 l=1,16
264       xxd=0.
265       do 48 k=1,16
266       xxd=xxd+iwt(k,l)*z(k)
267    48 continue
268       cl(l)=xxd
269    49 continue
270       l=0
271       do 51 k=1,4
272       do 50 j=1,4
273       l=l+1
274       cc(imem,n,m,k,j)=cl(l)
275    50 continue
276    51 continue
277    52 continue
278    53 continue
279       return
280       end
281
282       subroutine jeppe2(i,x,y,nx,my,xx,yy,cc,z)
283 C--   G.W. 02/07/2007 Allow extrapolation to small x and large q.
284       implicit real*8(a-h,o-z)
285       parameter(nhess=30)
286       dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)      
287       
288       n=locx(xx,nx,x)
289       m=locx(yy,my,y)
290       
291       if (n.gt.0.and.n.lt.nx.and.m.gt.0.and.m.lt.my) then
292 C--   Do usual interpolation.
293          t=(x-xx(n))/(xx(n+1)-xx(n))
294          u=(y-yy(m))/(yy(m+1)-yy(m))
295          z=0.d0
296          do l=4,1,-1
297             z=t*z+((cc(i,n,m,l,4)*u+cc(i,n,m,l,3))*u
298      &           +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
299          enddo
300          
301       else if (n.eq.0.and.m.gt.0.and.m.lt.my) then
302 C--   Extrapolate to small x.
303          call jeppe3(i,xx(1),y,nx,my,xx,yy,cc,f0)
304          call jeppe3(i,xx(2),y,nx,my,xx,yy,cc,f1)
305          if (f0.gt.0.d0.and.f1.gt.0.d0) then
306             z = exp(log(f0)+(log(f1)-log(f0))/(xx(2)-xx(1))*(x-xx(1)))
307          else
308             z = f0+(f1-f0)/(xx(2)-xx(1))*(x-xx(1))
309          end if
310          
311       else if (n.gt.0.and.m.eq.my) then
312 C--   Extrapolate to large q.
313          call jeppe3(i,x,yy(my),nx,my,xx,yy,cc,f0)
314          call jeppe3(i,x,yy(my-1),nx,my,xx,yy,cc,f1)
315          if (f0.gt.0.d0.and.f1.gt.0.d0) then
316             z = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
317      &           (y-yy(my)))
318          else
319             z = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
320          end if
321          
322       else if (n.eq.0.and.m.eq.my) then
323 C--   Extrapolate to small x AND large q.
324          call jeppe3(i,xx(1),yy(my),nx,my,xx,yy,cc,f0)
325          call jeppe3(i,xx(1),yy(my-1),nx,my,xx,yy,cc,f1)
326          if (f0.gt.0.d0.and.f1.gt.0.d0) then
327             z0 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
328      &           (y-yy(my)))
329          else
330             z0 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
331          end if
332          call jeppe3(i,xx(2),yy(my),nx,my,xx,yy,cc,f0)
333          call jeppe3(i,xx(2),yy(my-1),nx,my,xx,yy,cc,f1)
334          if (f0.gt.0.d0.and.f1.gt.0.d0) then
335             z1 = exp(log(f0)+(log(f0)-log(f1))/(yy(my)-yy(my-1))*
336      &           (y-yy(my)))
337          else
338             z1 = f0+(f0-f1)/(yy(my)-yy(my-1))*(y-yy(my))
339          end if
340          if (z0.gt.0.d0.and.z1.gt.0.d0) then
341             z = exp(log(z0)+(log(z1)-log(z0))/(xx(2)-xx(1))*(x-xx(1)))
342          else
343             z = z0+(z1-z0)/(xx(2)-xx(1))*(x-xx(1))
344          end if
345
346       else
347 C--   Set parton distribution to zero otherwise.
348          z = 0.d0
349
350       end if
351       
352       return
353       end
354
355 C--   G.W. 02/07/2007 Copy of the original jeppe2,
356 C--   only used for extrapolation.
357       subroutine jeppe3(i,x,y,nx,my,xx,yy,cc,z)
358       implicit real*8(a-h,o-z)
359       parameter(nhess=30)
360       dimension xx(nx),yy(my),cc(0:nhess,nx,my,4,4)      
361       n=locx(xx,nx,x)
362       m=locx(yy,my,y)
363       t=(x-xx(n))/(xx(n+1)-xx(n))
364       u=(y-yy(m))/(yy(m+1)-yy(m))
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
368      &        +cc(i,n,m,l,2))*u+cc(i,n,m,l,1)
369       enddo
370       return
371       end
372       
373       integer function locx(xx,nx,x)
374       implicit real*8(a-h,o-z)
375       dimension xx(nx)
376 c$$$      if(x.le.xx(1)) then
377       if(x.eq.xx(1)) then ! G.W. 02/07/2007
378       locx=1
379       return
380       endif
381 c$$$      if(x.ge.xx(nx)) then 
382       if(x.eq.xx(nx)) then ! G.W. 02/07/2007
383       locx=nx-1  
384       return
385       endif
386       ju=nx+1
387       jl=0
388     1 if((ju-jl).le.1) go to 2
389       jm=(ju+jl)/2
390       if(x.ge.xx(jm)) then
391       jl=jm
392       else
393       ju=jm
394       endif
395       go to 1
396     2 locx=jl
397       return
398       end
399
400       real*8 function  polderiv(x1,x2,x3,y1,y2,y3)
401       implicit real*8(a-h,o-z)
402       polderiv=(x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*
403      .(y2-y3))+x2*x2*(y1-y3)+x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3))
404       return
405       end