Reading QA thresholds from external file II
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapmrstqed.f
1       subroutine MRSTqedevolve(x,Q,pdf,photon)
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=9,nqc0=2,nqb0=11,nqc=35,nqb=26,
10      .nhess=30)
11       real*8 pdf(-6:6),photon,phot
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      .,f9(nx,nq)
21      .,fc(nx,nqc),fb(nx,nqb)
22       real*8 qq(nq),xx(nx),
23      .cc1(0:nhess,nx,nq,4,4),cc2(0:nhess,nx,nq,4,4),
24      .cc3(0:nhess,nx,nq,4,4),cc4(0:nhess,nx,nq,4,4),
25      .cc6(0:nhess,nx,nq,4,4),cc8(0:nhess,nx,nq,4,4),
26      .cc9(0:nhess,nx,nq,4,4),
27      .ccc(0:nhess,nx,nqc,4,4),ccb(0:nhess,nx,nqb,4,4)
28       real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
29       data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
30      .        1d-4,2d-4,4d-4,6d-4,8d-4,
31      .        1d-3,2d-3,4d-3,6d-3,8d-3,
32      .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
33      .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
34      .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
35      .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
36      .     .8d0,.9d0,1d0/
37       data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
38      .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
39      .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
40      .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
41      .        1.8d6,3.2d6,5.6d6,1d7/
42       data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
43       save
44       
45       xsave=x
46       qsq = q*q
47       q2save=qsq
48
49       xlog=dlog(x)
50       qsqlog=dlog(qsq)
51
52       call getnset(iset)
53 c      imem=member(iset)
54       call getnmem(iset,imem)
55       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
56       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
57       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
58       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
59       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
60       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
61       call jeppe2(imem,xlog,qsqlog,nx,nq,xxl,qql,cc9,photon)
62
63       chm=0.d0
64       if(qsq.gt.emc2) then 
65       call jeppe2(imem,xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
66       endif
67
68       bot=0.d0
69       if(qsq.gt.emb2) then 
70       call jeppe2(imem,xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
71       endif
72
73       pdf(0)  = glu
74       pdf(1)  = dnv+dsea
75       pdf(-1) = dsea
76       pdf(2)  = upv+usea
77       pdf(-2) = usea
78       pdf(3)  = str
79       pdf(-3) = str
80       pdf(4)  = chm
81       pdf(-4) = chm
82       pdf(5)  = bot
83       pdf(-5) = bot
84       pdf(6) = 0.0d0
85       pdf(-6) = 0.0d0
86       
87       x=xsave
88       qsq=q2save
89       return
90 *
91       entry MRSTqedread(nset)
92       read(1,*)nmem(nset),ndef(nset)
93 c      print *,nmem(nset),ndef(nset)
94 c      do nm = 0,nmem-1
95       do nm = 0,nmem(nset)
96         do 20 n=1,nx-1
97         do 20 m=1,nq
98         read(1,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
99      .            f5(n,m),f7(n,m),f6(n,m),f8(n,m),f9(n,m)
100 c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea +9 - photon
101   20  continue
102 c      write(*,*)'PDF set ',nm,' first element ',f1(1,1)
103       do 40 m=1,nq
104       f1(nx,m)=0.d0
105       f2(nx,m)=0.d0
106       f3(nx,m)=0.d0
107       f4(nx,m)=0.d0
108       f5(nx,m)=0.d0
109       f6(nx,m)=0.d0
110       f7(nx,m)=0.d0
111       f8(nx,m)=0.d0
112       f9(nx,m)=0.d0
113   40  continue
114
115       do n=1,nx
116       xxl(n)=dlog(xx(n))
117       enddo
118       do m=1,nq
119       qql(m)=dlog(qq(m))
120       enddo
121
122       call jeppe1(nm,nx,nq,xxl,qql,f1,cc1)
123       call jeppe1(nm,nx,nq,xxl,qql,f2,cc2)
124       call jeppe1(nm,nx,nq,xxl,qql,f3,cc3)
125       call jeppe1(nm,nx,nq,xxl,qql,f4,cc4)
126       call jeppe1(nm,nx,nq,xxl,qql,f6,cc6)
127       call jeppe1(nm,nx,nq,xxl,qql,f8,cc8)
128       call jeppe1(nm,nx,nq,xxl,qql,f9,cc9)
129
130       emc2=2.045
131       emb2=18.5
132
133       do 44 m=1,nqc
134       qqlc(m)=qql(m+nqc0)
135       do 44 n=1,nx
136       fc(n,m)=f5(n,m+nqc0)
137    44 continue
138       qqlc(1)=dlog(emc2)
139       call jeppe1(nm,nx,nqc,xxl,qqlc,fc,ccc)
140
141       do 45 m=1,nqb
142       qqlb(m)=qql(m+nqb0)
143       do 45 n=1,nx
144       fb(n,m)=f7(n,m+nqb0)
145    45 continue
146       qqlb(1)=dlog(emb2)
147       call jeppe1(nm,nx,nqb,xxl,qqlb,fb,ccb)
148
149
150       enddo
151   50  format(9f10.5)
152       return
153 *    
154       end