← Back to branch summary

~maddevelopers/mg5amcnlo/WWW5_caching

« back to all changes in this revision

Viewing changes to users/mardelcourt/PROC_141512/PROC_141512/Source/PDF/mrst2002.f

  • Committer: John Doe
  • Date: 2013-03-25 20:27:02 UTC
  • Revision ID: john.doe@gmail.com-20130325202702-5sk3t1r8h33ca4p4
first clean version

Show diffs side-by-side

added added

removed removed

Lines of Context:
users/mardelcourt/PROC_141512/PROC_141512/Source/PDF/mrst2002.f
 
1
      subroutine mrst2002(x,q,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
 
2
C***************************************************************C
 
3
C                                                               C
 
4
C  This is a package for the new MRST 2002 updated NLO and      C
 
5
C  NNLO parton distributions.                                   C 
 
6
C  Reference: A.D. Martin, R.G. Roberts, W.J. Stirling and      C
 
7
C  R.S. Thorne, hep-ph/0211080                                  C
 
8
C                                                               C
 
9
C  There are 2 pdf sets corresponding to mode = 1, 2            C
 
10
C                                                               C
 
11
C  Mode=1 gives the NLO set with alpha_s(M_Z,NLO) = 0.1197      C  
 
12
C  This set reads a grid whose first number is 0.00949          C
 
13
C                                                               C
 
14
C  Mode=2 gives the NNLO set with alpha_s(M_Z,NNLO) = 0.1154    C
 
15
C  This set reads a grid whose first number is 0.00685          C
 
16
C                                                               C
 
17
C         Comments to : W.J.Stirling@durham.ac.uk               C
 
18
C                                                               C
 
19
C***************************************************************C
 
20
      implicit real*8(a-h,o-z)
 
21
      data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
 
22
      q2=q*q
 
23
      if(q2.lt.qsqmin.or.q2.gt.qsqmax) print 99,q2
 
24
      if(x.lt.xmin.or.x.gt.xmax)       print 98,x
 
25
          if(mode.eq.1) then
 
26
        call mrst_02_1(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu) 
 
27
      elseif(mode.eq.2) then
 
28
        call mrst_02_1(x,q2,upv,dnv,usea,dsea,str,chm,bot,glu) 
 
29
      endif 
 
30
  99  format('  WARNING:  Q^2 VALUE IS OUT OF RANGE   ','q2= ',e10.5)
 
31
  98  format('  WARNING:   X  VALUE IS OUT OF RANGE   ','x= ',e10.5)
 
32
      return
 
33
      end
 
34
 
 
35
      subroutine mrst_02_1(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
 
36
      implicit real*8(a-h,o-z)
 
37
c-fabio
 
38
      Character Tablefile*40
 
39
      data Tablefile/'mrst2002nlo.dat'/
 
40
      integer IU
 
41
      common/IU/IU
 
42
c
 
43
      parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
 
44
      real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
 
45
     .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
 
46
      real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
 
47
     .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
 
48
     .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
 
49
      real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
 
50
      data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
 
51
     .        1d-4,2d-4,4d-4,6d-4,8d-4,
 
52
     .        1d-3,2d-3,4d-3,6d-3,8d-3,
 
53
     .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
 
54
     .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
 
55
     .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
 
56
     .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
 
57
     .     .8d0,.9d0,1d0/
 
58
      data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
 
59
     .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
 
60
     .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
 
61
     .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
 
62
     .        1.8d6,3.2d6,5.6d6,1d7/
 
63
      data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
 
64
      data init/0/
 
65
      save
 
66
      xsave=x
 
67
      q2save=qsq
 
68
      if(init.ne.0) goto 10
 
69
c        open(unit=33,file='mrst2002nlo.dat',status='old')
 
70
        call OpenData(Tablefile)
 
71
        do 20 n=1,nx-1
 
72
        do 20 m=1,nq
 
73
        read(IU,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
 
74
     .            f5(n,m),f7(n,m),f6(n,m),f8(n,m)
 
75
c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
 
76
  20  continue
 
77
      do 40 m=1,nq
 
78
      f1(nx,m)=0.d0
 
79
      f2(nx,m)=0.d0
 
80
      f3(nx,m)=0.d0
 
81
      f4(nx,m)=0.d0
 
82
      f5(nx,m)=0.d0
 
83
      f6(nx,m)=0.d0
 
84
      f7(nx,m)=0.d0
 
85
      f8(nx,m)=0.d0
 
86
  40  continue
 
87
      do n=1,nx
 
88
      xxl(n)=dlog(xx(n))
 
89
      enddo
 
90
      do m=1,nq
 
91
      qql(m)=dlog(qq(m))
 
92
      enddo
 
93
 
 
94
      call jeppe1(nx,nq,xxl,qql,f1,cc1)
 
95
      call jeppe1(nx,nq,xxl,qql,f2,cc2)
 
96
      call jeppe1(nx,nq,xxl,qql,f3,cc3)
 
97
      call jeppe1(nx,nq,xxl,qql,f4,cc4)
 
98
      call jeppe1(nx,nq,xxl,qql,f6,cc6)
 
99
      call jeppe1(nx,nq,xxl,qql,f8,cc8)
 
100
 
 
101
      emc2=2.045
 
102
      emb2=18.5
 
103
 
 
104
      do 44 m=1,nqc
 
105
      qqlc(m)=qql(m+nqc0)
 
106
      do 44 n=1,nx
 
107
      fc(n,m)=f5(n,m+nqc0)
 
108
   44 continue
 
109
      qqlc(1)=dlog(emc2)
 
110
      call jeppe1(nx,nqc,xxl,qqlc,fc,ccc)
 
111
 
 
112
      do 45 m=1,nqb
 
113
      qqlb(m)=qql(m+nqb0)
 
114
      do 45 n=1,nx
 
115
      fb(n,m)=f7(n,m+nqb0)
 
116
   45 continue
 
117
      qqlb(1)=dlog(emb2)
 
118
      call jeppe1(nx,nqb,xxl,qqlb,fb,ccb)
 
119
 
 
120
 
 
121
      init=1
 
122
   10 continue
 
123
      
 
124
      xlog=dlog(x)
 
125
      qsqlog=dlog(qsq)
 
126
 
 
127
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
 
128
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
 
129
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
 
130
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
 
131
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
 
132
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
 
133
 
 
134
      chm=0.d0
 
135
      if(qsq.gt.emc2) then 
 
136
      call jeppe2(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
 
137
      endif
 
138
 
 
139
      bot=0.d0
 
140
      if(qsq.gt.emb2) then 
 
141
      call jeppe2(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
 
142
      endif
 
143
 
 
144
      x=xsave
 
145
      qsq=q2save
 
146
      return
 
147
   50 format(8f10.5)
 
148
      end
 
149
 
 
150
      subroutine mrst_02_2(x,qsq,upv,dnv,usea,dsea,str,chm,bot,glu)
 
151
      implicit real*8(a-h,o-z)
 
152
c-fabio
 
153
      Character Tablefile*40
 
154
      data Tablefile/'mrst2002nnlo.dat'/
 
155
      integer IU
 
156
      common/IU/IU
 
157
c
 
158
      parameter(nx=49,nq=37,np=8,nqc0=2,nqb0=11,nqc=35,nqb=26)
 
159
      real*8 f1(nx,nq),f2(nx,nq),f3(nx,nq),f4(nx,nq),f5(nx,nq),
 
160
     .f6(nx,nq),f7(nx,nq),f8(nx,nq),fc(nx,nqc),fb(nx,nqb)
 
161
      real*8 qq(nq),xx(nx),cc1(nx,nq,4,4),cc2(nx,nq,4,4),
 
162
     .cc3(nx,nq,4,4),cc4(nx,nq,4,4),cc6(nx,nq,4,4),cc8(nx,nq,4,4),
 
163
     .ccc(nx,nqc,4,4),ccb(nx,nqb,4,4)
 
164
      real*8 xxl(nx),qql(nq),qqlc(nqc),qqlb(nqb)
 
165
      data xx/1d-5,2d-5,4d-5,6d-5,8d-5,
 
166
     .        1d-4,2d-4,4d-4,6d-4,8d-4,
 
167
     .        1d-3,2d-3,4d-3,6d-3,8d-3,
 
168
     .        1d-2,1.4d-2,2d-2,3d-2,4d-2,6d-2,8d-2,
 
169
     .     .1d0,.125d0,.15d0,.175d0,.2d0,.225d0,.25d0,.275d0,
 
170
     .     .3d0,.325d0,.35d0,.375d0,.4d0,.425d0,.45d0,.475d0,
 
171
     .     .5d0,.525d0,.55d0,.575d0,.6d0,.65d0,.7d0,.75d0,
 
172
     .     .8d0,.9d0,1d0/
 
173
      data qq/1.25d0,1.5d0,2d0,2.5d0,3.2d0,4d0,5d0,6.4d0,8d0,1d1,
 
174
     .        1.2d1,1.8d1,2.6d1,4d1,6.4d1,1d2,
 
175
     .        1.6d2,2.4d2,4d2,6.4d2,1d3,1.8d3,3.2d3,5.6d3,1d4,
 
176
     .        1.8d4,3.2d4,5.6d4,1d5,1.8d5,3.2d5,5.6d5,1d6,
 
177
     .        1.8d6,3.2d6,5.6d6,1d7/
 
178
      data xmin,xmax,qsqmin,qsqmax/1d-5,1d0,1.25d0,1d7/
 
179
      data init/0/
 
180
      save
 
181
      xsave=x
 
182
      q2save=qsq
 
183
      if(init.ne.0) goto 10
 
184
c        open(unit=33,file='mrst2002nnlo.dat',status='old')
 
185
        call OpenData(Tablefile)
 
186
        do 20 n=1,nx-1
 
187
        do 20 m=1,nq
 
188
        read(IU,50)f1(n,m),f2(n,m),f3(n,m),f4(n,m),
 
189
     .            f5(n,m),f7(n,m),f6(n,m),f8(n,m)
 
190
c notation: 1=uval 2=val 3=glue 4=usea 5=chm 6=str 7=btm 8=dsea
 
191
  20  continue
 
192
      do 40 m=1,nq
 
193
      f1(nx,m)=0.d0
 
194
      f2(nx,m)=0.d0
 
195
      f3(nx,m)=0.d0
 
196
      f4(nx,m)=0.d0
 
197
      f5(nx,m)=0.d0
 
198
      f6(nx,m)=0.d0
 
199
      f7(nx,m)=0.d0
 
200
      f8(nx,m)=0.d0
 
201
  40  continue
 
202
      do n=1,nx
 
203
      xxl(n)=dlog(xx(n))
 
204
      enddo
 
205
      do m=1,nq
 
206
      qql(m)=dlog(qq(m))
 
207
      enddo
 
208
 
 
209
      call jeppe1(nx,nq,xxl,qql,f1,cc1)
 
210
      call jeppe1(nx,nq,xxl,qql,f2,cc2)
 
211
      call jeppe1(nx,nq,xxl,qql,f3,cc3)
 
212
      call jeppe1(nx,nq,xxl,qql,f4,cc4)
 
213
      call jeppe1(nx,nq,xxl,qql,f6,cc6)
 
214
      call jeppe1(nx,nq,xxl,qql,f8,cc8)
 
215
 
 
216
      emc2=2.045
 
217
      emb2=18.5
 
218
 
 
219
      do 44 m=1,nqc
 
220
      qqlc(m)=qql(m+nqc0)
 
221
      do 44 n=1,nx
 
222
      fc(n,m)=f5(n,m+nqc0)
 
223
   44 continue
 
224
      qqlc(1)=dlog(emc2)
 
225
      call jeppe1(nx,nqc,xxl,qqlc,fc,ccc)
 
226
 
 
227
      do 45 m=1,nqb
 
228
      qqlb(m)=qql(m+nqb0)
 
229
      do 45 n=1,nx
 
230
      fb(n,m)=f7(n,m+nqb0)
 
231
   45 continue
 
232
      qqlb(1)=dlog(emb2)
 
233
      call jeppe1(nx,nqb,xxl,qqlb,fb,ccb)
 
234
 
 
235
 
 
236
      init=1
 
237
   10 continue
 
238
      
 
239
      xlog=dlog(x)
 
240
      qsqlog=dlog(qsq)
 
241
 
 
242
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc1,upv)
 
243
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc2,dnv)
 
244
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc3,glu)
 
245
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc4,usea)
 
246
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc6,str)
 
247
      call jeppe2(xlog,qsqlog,nx,nq,xxl,qql,cc8,dsea)
 
248
 
 
249
      chm=0.d0
 
250
      if(qsq.gt.emc2) then 
 
251
      call jeppe2(xlog,qsqlog,nx,nqc,xxl,qqlc,ccc,chm)
 
252
      endif
 
253
 
 
254
      bot=0.d0
 
255
      if(qsq.gt.emb2) then 
 
256
      call jeppe2(xlog,qsqlog,nx,nqb,xxl,qqlb,ccb,bot)
 
257
      endif
 
258
 
 
259
      x=xsave
 
260
      qsq=q2save
 
261
      return
 
262
   50 format(8f10.5)
 
263
      end
 
264