EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
parton.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file parton.F
1 *72*********************************************************************
2  SUBROUTINE parton(X,Q2,XPQ,XDPQ)
3 C *********************************************************************
4 C This subroutine gives the parton distributions (polarized
5 C and unpolarized) for the polarized leptoproduction. The
6 C parton distributions are optimized for polarized scattering.
7 C The following sets are included
8 C
9 C LST(15) = 101 Schaefer, Phys. Lett. B 208,2 (1988) 175
10 C for comparison with older PEPSI versions
11 C LST(15) = 102 free
12 C LST(15) = 103 free
13 C LST(15) = 104 Schaefer et al hep-ph/9505306
14 C using Glueck et al. Z. Phys. C 53 (1992) 127
15 C LST(15) = 105 Bartelski et al hep-ph/9502271 Set II(p,n)
16 C LST(15) = 106 Bartelski et al hep-ph/9502271 Set II(P,n)
17 C LST(15) = 107 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
18 C LST(15) = 108 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
19 C LST(15) = 109 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
20 C LST(15) = 110 Gehrmann et al hep-ph/9512406 Gluon A (LO)
21 C LST(15) = 111 Gehrmann et al hep-ph/9512406 Gluon B (LO)
22 C LST(15) = 112 Gehrmann et al hep-ph/9512406 Gluon C (LO)
23 C LST(15) = 113 Gehrmann et al hep-ph/9512406 Gluon A (LO) + (DGLAP)
24 C LST(15) = 114 Gehrmann et al hep-ph/9512406 Gluon B (LO) + (DGLAP)
25 C LST(15) = 115 Gehrmann et al hep-ph/9512406 Gluon C (LO) + (DGLAP)
26 C...............................................................................
27 C M. Glueck, E. Reya, M. Stratmann and W. Vogelsang,
28 C DO-TH 95/13, RAL-TR-95-042
29 C LST(15) = 116 'standard' scenario, next-to-leading order
30 C LST(15) = 117 'valence' scenario, next-to-leading order
31 C LST(15) = 118 'standard' scenario, leading order
32 C LST(15) = 119 'valence' scenario, leading order
33 C LST(15) = 120 Stanley J.Brodsky Nucl.Phys. B441(1995)
34 C LST(15) = 121 S.Keler & J.F.Owens Phys.Lett. B266(1991)
35 C & Phys.Rev. D19(1994)1199
36 C LST(15) = 124 D. de Florian et al., hep-ph/9711440 LO set 1
37 C LST(15) = 125 D. de Florian et al., hep-ph/9711440 LO set 2
38 C LST(15) = 126 D. de Florian et al., hep-ph/9711440 LO set 3
39 C LST(15) = 127 D. de Florian et al., hep-ph/9711440 NLO set 1
40 C LST(15) = 128 D. de Florian et al., hep-ph/9711440 NLO set 2
41 C LST(15) = 129 D. de Florian et al., hep-ph/9711440 NLO set 3
42 C LST(15) = 130 Fake sample :
43 C unpolarized Gehrmann et al hep-ph/9512406
44 C with Delta u(x) = 0.5 * u(x) and Delta d(x) = 0.
45 C LST(15) = 131 Fake sample :
46 C unpolarized Gehrmann et al hep-ph/9512406
47 C with Delta d(x) = 0.5 * d(x) and Delta u(x) = 0.
48 C LST(15) = 132 fit routine. (Is outside the official code.)
49 C
50 C CTEQ4 collaboration: http:
51 C UNPOL: Low Q2 parametrization is the only one used here
52 C POL: BOGUS, du=0.5* u(x) dd=-0.3*d(x) 0.0 all else
53 C LST(15) = 133 CTEQ4LQ
54 C
55 C MRS low Q2
56 C UNPOL MRS 95 low Q2
57 C POL: BOGUS, du=0.5* u(x) dd=-0.3*d(x) 0.0 all else
58 C LST(15)= 137 MRS low Q2
59 C LST(15) = 144 grsv2000 hep-ph/0011215 LO standard scenario
60 C LST(15) = 145 grsv2000 hep-ph/0011215 LO valence scenario
61 C LST(15) = 146 grsv2000 hep-ph/0011215 NLO standard scenario
62 C LST(15) = 147 grsv2000 hep-ph/0011215 NLO valence scenario
63 C..............................................................................
64 C unpolarized distributions
65 C
66 C LST(15) = 150 cteq5l LO
67 C LST(15) = 151 cteq5m NLO MSBAR
68 C LST(15) = 152 cteq5m1 NLO MSBAR (update)
69 C
70 C LST(15) = 161 mrs99 cor01 central gluon, a_s
71 C LST(15) = 162 mrs99 cor02 higher gluon
72 C LST(15) = 163 mrs99 cor03 lower gluon
73 C LST(15) = 164 mrs99 cor04 lower a_s
74 C LST(15) = 165 mrs99 cor05 higher a_s
75 C LST(15) = 166 mrs99 cor06 quarks up
76 C LST(15) = 167 mrs99 cor07 quarks down
77 C LST(15) = 168 mrs99 cor08 strange up
78 C LST(15) = 169 mrs99 cor09 strange down
79 C LST(15) = 170 mrs99 cor10 charm up
80 C LST(15) = 171 mrs99 cor11 charm down
81 C LST(15) = 172 mrs99 cor12 larger d/u
82 C
83 C LST(15) = 173 cteq6l LO
84 C LST(15) = 174 cteq6d DIS NLO
85 C LST(15) = 175 cteq6m NLO MSBAR
86 ***********************************************************************
87 
88 *
89 * to avoid variable conflictions, a second keep element is necessary
90 * with the same common block name (see LEPTO2)
91 *
92 
93  include "mc_set.inc"
94 
95  COMMON /leptou/ cut(14),lst(40),parl(30),
96  & xlp,ylp,w2lp,q2lp,ulp
97  REAL cut,parl,xlp,ylp,w2lp,q2lp,ulp
98  INTEGER lst
99  SAVE /leptou/
100 
101  INTEGER imxpdf
102  parameter(imxpdf=40)
103  COMMON /pepadm/cpdfnam(2,imxpdf),ipdfnam(2,imxpdf),
104  & iplst(10),cunpol,cpol
105  CHARACTER*256 cpdfnam,cunpol,cpol
106  INTEGER iplst,ipdfnam
107  SAVE /pepadm/
108 
109 **************************************************************
110 *
111 * IPLST(1) = 0 (default) : number of PDF warnings
112 * IPLST(2) = 11 (default) : unit -1- for pdf files
113 * IPLST(3) = 12 (default) : unit -2- for pdf files
114 *
115 **************************************************************
116 
117 
118  REAL x,q2,xpq(-6:6),xdpq(-6:6),ulmass
119  DOUBLE PRECISION dx,dq2
120  REAL xuv, xdv, xs, xg
121  REAL xub, xdb, xc, xbot
122  REAL xduv,xddv,xds,xdg
123  REAL xdub,xddb,xdc
124 
125 C...Variables for set 101
126  REAL valphau,vbetau,gamma
127  REAL valphad,vbetad
128  REAL vgamma,vc3
129  REAL vau0,vau1,vad0,vad1
130  REAL ffu0,ffu1,ffd0,ffd1,fg1p,fg1n
131 
132 C.. Variables for set 104
133  REAL vs,vn,vka,vsw
134  REAL valpha,vbeta
135  REAL vkleina,vkleinb
136  REAL va,vb,vc,vd,ve,vf
137  REAL va0,va10
138  REAL mnucl,www2
139 
140  REAL xd,xu,xdd,xdu,xdds,xdus
141 
142 C...Variables for sets 107-109
143  DOUBLE PRECISION mrsupv,mrsdnv,mrsdsea,mrsusea
144  DOUBLE PRECISION mrsstr,mrschm,mrsbot,mrsglu
145  INTEGER mode
146 
147 C...Variables for sets 122-129
148  DOUBLE PRECISION ddxuv,ddxdv,ddxdel,ddxudb,ddxsb,ddxgl,
149  & dxduv,dxddv,dxdqbar,dxdstr,dxdglu,
150  & ddxus,ddxds
151 C...Variables for sets 105 and 106
152  REAL fdelta,fm
153  REAL va1,va2,va3
154  REAL vc1,vc2
155 
156 C...Variables for set 110-112
157  REAL vetau ,vetad, vetaqb, vetag
158  REAL vnormu,vnormd,vnormqb,vnormg
159  REAL vau, vad, vaqb, vag
160  REAL vbu, vbd, vbqb, vbg
161  REAL vgu, vgd, vgqb, vgg
162  REAL vru, vrd, vrqb, vrg
163 
164 C...Variables for set 113-115
165  INTEGER iflag
166  DOUBLE PRECISION dxdg,dxdubar,dxddbar,dxds
167  REAL xupdb,xdmub
168 
169  REAL f2n,f2p
170 
171 C...Variables for sets 133-36
172  DOUBLE PRECISION ddxpq(-6:6)
173  INTEGER initcteq,cteqset
174  DOUBLE PRECISION ctq4pdf
175 
176 C...Variables for set 116-119
177  INTEGER iset
178 
179 
180 C...Variables for set 120-121
181  REAL y,x1,ss,z
182  REAL aj(4),bj(4),cj(4),rj(4),w(3),xw(3),ww(3),aaa(6,4)
183  REAL const0(6,4),const1(6,4),const2(6,4),const3(6,4),
184  & sumr(2),qg(4),qj(4)
185  REAL aq(4),bq(4),cq(4),dq(4),st1(4),st2(4),
186  & qp(4),qm(4),st3(4),st4(4)
187  DATA xw/-1.775e-02, -1.e-02, -2.2540e-03/,w/5.556e-03,
188  & 8.889e-03, 5.556e-03/,ww/0.3872983, 3.1329002e-02,
189  & 4.849505/,aj/ -0.2, -0.2, -0.3, -0.4/,bj/ 4.7, 3.7,
190  & 8., 8./,cj/ 0., 2.54, 0., 0./,rj/-0.9582284, 2.068491,
191  & -0.3953224, 14.54266/,aq/0.757, 3.784, 0.2897, 2./,
192  & bq/-0.645, -3.672, -0.2637, -1.25/,cq/3.23, 2.004,
193  & 1., 2./,dq/-3.118, -1.892, -0.9748, -1.25/,st1/3., 3.,
194  & 5., 4./,st2/4., 4., 6., 5./,st3/5., 5., 7., 6./,
195  & st4/6., 6., 8., 7./
196  DATA const0/-0.335 , 3.614 , 0.8673 , 0.0 , 0.0 , 0.0 ,
197  & -0.1612 , 4.667 , 0.0 , 0.0 , 0.0 , 0.0 ,
198  & -1.0 , 7.278 , 0.0 , 0.0 , 0.0 , 0.909,
199  & -1.0 , 5.304 , 0.0 , 0.0 , 0.0 , 3.017/
200  DATA const1/-0.1097 , .8395 ,-1.6637 , 1.1049, 0.0 , 0.0,
201  & -0.2092 , 0.7951,-1.0232 , 0.8616, 0.0 , 0.0,
202  & -0.3823,-0.7904,-1.6629,-0.0133,0.1211,-0.4023,
203  & -0.9342, 1.4654,-3.9141,9.0176,-5.9602,-4.7347/
204  DATA const2/-0.002442,-0.02186,0.342,-0.2369,0.0,0.0,
205  & 0.02657, 0.1081, 0.05799, 0.153,0.0,0.0,
206  & 0.02766, 0.8108,0.5719,0.5299,-0.1739,0.0063,
207  & 0.5454, -1.4292,2.8445,-10.426,7.515,3.3594/
208  DATA const3/ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
209  & 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
210  & 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 0.0 ,
211  & -0.1668,0.7569,-0.8411,4.0983,-2.7329,-0.9443/
212  DATA initcteq/ 1 /
213  DATA cteqset/10/
214 C.. Variables for helium amendment
215  REAL percs,percsp,percd
216  REAL admxs,admxsp,admxd
217  REAL facts,factsp,factd
218  REAL pari11
219 C
220 C.. Variables for set 144-147
221  DOUBLE PRECISION dxdd,dxdu,dummyg1p, dummyg1n
222  DOUBLE PRECISION dxdub,dxddb
223 
224 C.. Variables for set 150-152
225  INTEGER iparton
226  DOUBLE PRECISION ctq5pdf
227 
228 C.. Variable for set 173-175
229 C.. Use IPARTON from definition above
230  DOUBLE PRECISION ctq6pdf
231 
232  INTEGER i,j,k,icount,icoun1,icerr
233  DATA icount /0/,icoun1/0/,icerr/0/
234 
235  icount = icount + 1
236 
237 
238 
239 C Reset all partondensities
240  DO i=-6,6
241  xpq(i) = 0.
242  xdpq(i) = 0.
243  ENDDO
244 
245 C ***************************************************************************
246 C
247 C set 101 Schaefer, Phys. Lett. B 208,2 (1988) 175
248 C
249 C ***************************************************************************
250 
251 
252  IF(lst(15).EQ.101) THEN
253 
254 C unpolarized part
255 
256  valphau = 0.588
257  vbetau = 2.69
258  xuv = gamma(valphau)*gamma(vbetau+1.)/gamma(valphau+vbetau+1.)
259  xuv = 2./xuv*x**valphau*(1.-x)**vbetau
260 
261  valphad = 1.03
262  vbetad = 6.87
263  xdv = gamma(valphad)*gamma(vbetad+1.)/gamma(valphad+vbetad+1.)
264  xdv = 1./xdv*x**valphad*(1.-x)**vbetad
265 
266  vgamma = 14.6
267  vc3 = 0.0147
268  xub = vc3*(1.+vgamma)*(1.-x)**vgamma
269 
270  xdb = xub
271  xs = xub
272  xc = 0.02*(1.-x)**19
273 
274  xpq( 1) = xdv + xdb
275  xpq(-1) = xdb
276  xpq( 2) = xuv + xub
277  xpq(-2) = xub
278  xpq( 3) = xs
279  xpq(-3) = xs
280  xpq( 4) = xc
281  xpq(-4) = xc
282 
283 C polarized part
284 
285  vau0 = 0.12
286  vau1 = 0.3*vau0
287  vad0 = 7.0*vau0
288  vad1 = 7.0*vau1
289  valphau = 0.588
290  valphad = 1.03
291 
292  ffu0 = 1./(1.+ vau0*x**(-valphau)*(1.-x)*(1.-x))
293  ffu1 = 1./(1.+ vau1*x**(-valphau)*(1.-x)*(1.-x))
294  ffd0 = 1./(1.+ vad0*x**(-valphad)*(1.-x)*(1.-x))
295  ffd1 = 1./(1.+ vad1*x**(-valphad)*(1.-x)*(1.-x))
296  fg1p = 0.5*(4*xuv*ffu0/9.
297  & - xdv*(4*ffu0+4*ffu1-ffd0+2*ffd1)/27.)
298  fg1n = 0.5*( xuv*ffd0/9.
299  & - xdv*(-4*ffu0+8*ffu1+ffd0+ffd1)/27.)
300 
301  xduv = 6*(4*fg1p-fg1n)/5.
302 
303  vau0 = 0.12
304  vau1 = 0.3*vau0
305  vad0 = 7*vau0
306  vad1 = 7*vau1
307 
308  valphau = 0.588
309  valphad = 1.03
310 
311  ffu0 = 1./(1.+ vau0*x**(-valphau)*(1.-x)*(1.-x))
312  ffu1 = 1./(1.+ vau1*x**(-valphau)*(1.-x)*(1.-x))
313  ffd0 = 1./(1.+ vad0*x**(-valphad)*(1.-x)*(1.-x))
314  ffd1 = 1./(1.+ vad1*x**(-valphad)*(1.-x)*(1.-x))
315  fg1p = 0.5*(4*xuv*ffu0/9.
316  & - xdv*( 4*ffu0+4*ffu1-ffd0+2*ffd1)/27.)
317  fg1n = 0.5*(xuv*ffd0/9.
318  & - xdv*(-4*ffu0+8*ffu1+ffd0+ffd1)/27.)
319 
320  xddv = 6*(4*fg1n-fg1p)/5.
321  xdub = 0.
322  xddb = 0.
323 
324  xdpq( 1) = xddv + xddb
325  xdpq(-1) = xddb
326  xdpq( 2) = xduv + xdub
327  xdpq(-2) = xdub
328 
329  ELSEIF(lst(15).EQ.104) THEN
330 C ***************************************************************************
331 C
332 C set 104 Schaefer et al hep-ph/9505306
333 C
334 C ***************************************************************************
335 
336 C unpolarized part
337  IF (((q2.LT.0.25 ).OR.(q2.GT.1.e+8)).OR.
338  & (( x.LT.1.e-5).OR.( x.GT. 1. ))) THEN
339  icoun1 = icoun1 + 1
340  IF(icoun1.lt.10) THEN
341  WRITE(*,*) 'WARNING: (Q^2,x) exceeds allowed range'
342  ENDIF
343  IF(lst(22).eq.0.or.lst(25).eq.0) THEN
344  xpq(2) = 1.
345  ENDIF
346  RETURN
347  ENDIF
348 
349  mnucl = ulmass(2212)
350  www2 = mnucl*mnucl + q2*(1./x-1.)
351 
352  vs = log(log(q2/(0.232*0.232))/log(0.25/(0.232*0.232)))
353 
354  vn = 0.579 + 0.283*vs + 0.047*vs*vs
355  vka = 0.523 - 0.015*vs
356  va = 2.22 - 0.59*vs - 0.27*vs*vs
357  vb = 5.95 - 6.19*vs + 1.55*vs*vs
358  vd = 3.57 + 0.94*vs - 0.16*vs*vs
359  xdv = vn*x**vka*(1.+va*sqrt(x)+ vb*x)*(1.-x)**vd
360 
361  vn = 0.663 + 0.191*vs - 0.041*vs*vs + 0.031*vs*vs*vs
362  vka = 0.326
363  va = -1.97 + 6.74*vs - 1.96*vs*vs
364  vb = 24.4 -20.7*vs + 4.08*vs*vs
365  vd = 2.86 + 0.7*vs - 0.02*vs*vs
366  xuv = vn*x**vka*(1.+va*sqrt(x)+ vb*x)*(1.-x)**vd - xdv
367 
368  valpha = 1.396
369  vbeta = 1.331
370  vkleina = 0.412 - 0.171*vs
371  vkleinb = 0.566 - 0.496*vs
372  va = 0.363
373  vb = -1.196
374  vc = 1.029 + 1.785*vs - 0.459*vs*vs
375  vd = 4.696 + 2.109*vs
376  ve = 3.838 + 1.944*vs
377  vf = 2.845
378  xub = ( x**vkleina*(va+vb*x+vc*x*x)*(log(1./x))**vkleinb
379  & +vs**valpha*exp(-ve+sqrt(vf*vs**vbeta*log(1./x))))
380  & *(1.-x)**vd
381  xdb = xub
382 
383  vsw = 0.
384  valpha = 0.803
385  vbeta = 0.563
386  vkleina = 2.082 - 0.577*vs
387  va = -3.055 + 1.024*vs**0.67
388  vb = 27.4 - 20.0*vs**0.154
389  vd = 6.22
390  ve = 4.33 + 1.408*vs
391  vf = 8.27 - 0.437*vs
392 
393  IF (vs.GT.vsw) THEN
394  xs = (vs-vsw)**valpha/(log(1./x)**vkleina)
395  & *(1. + va*sqrt(x) + vb*x)*(1.-x)**vd
396  & *exp(-ve + sqrt(vf*vs**vbeta*log(1./x)))
397  ELSE
398  xs = 0
399  ENDIF
400 
401  vsw = 0.888
402  valpha = 1.01
403  vbeta = 0.37
404  vkleina = 0.
405  va = 0.
406  vb = 4.24 - 0.804*vs
407  vd = 3.46 + 1.076*vs
408  ve = 4.61 + 1.49*vs
409  vf = 2.555 + 1.961*vs
410 
411 *MM >> if w > 2* m_c = 3GeV
412  IF (www2.GT.9.AND.vs.GT.vsw) THEN
413  xc = (vs-vsw)**valpha/(log(1./x)**vkleina)
414  & *(1. + va*sqrt(x) + vb*x)*(1.-x)**vd
415  & *exp(-ve + sqrt(vf*vs**vbeta*log(1./x)))
416  ELSE
417  xc =0.
418  ENDIF
419 
420  vsw = 1.351
421  valpha = 1.00
422  vbeta = 0.51
423  vkleina = 0.
424  va = 0.
425  vb = 1.848
426  vd = 2.929 + 1.396*vs
427  ve = 4.71 + 1.514*vs
428  vf = 4.02 + 1.239*vs
429 
430 *MM >> if w > 2*m_b = 10 GeV
431  IF (www2.GT.100.AND.vs.GT.vsw) THEN
432  xbot = (vs-vsw)**valpha/(log(1./x)**vkleina)
433  & *(1. + va*sqrt(x) + vb*x)*(1.-x)**vd
434  & *exp(-ve + sqrt(vf*vs**vbeta*log(1./x)))
435  ELSE
436  xbot =0.
437  ENDIF
438 
439  valpha = 0.558
440  vbeta = 1.218
441  vkleina = 1.00 - 0.17*vs
442  vkleinb = 0.
443  va = 4.879*vs - 1.383*vs*vs
444  vb = 25.92 - 28.97*vs + 5.596*vs*vs
445  vc = -25.69 + 23.68*vs - 1.975*vs*vs
446  vd = 2.537 + 1.718*vs + 0.353*vs*vs
447  ve = 0.595 + 2.138*vs
448  vf = 4.066
449  xg = ( x**vkleina*(va+vb*x+vc*x*x)
450  & +vs**valpha*exp(-ve+sqrt(vf*vs**vbeta*log(1./x))))
451  & *(1-x)**vd
452 
453  xpq( 0) = xg
454  xpq( 1) = xdv + xdb
455  xpq(-1) = xdb
456  xpq( 2) = xuv + xub
457  xpq(-2) = xub
458  xpq( 3) = xs
459  xpq(-3) = xs
460  xpq( 4) = xc
461  xpq(-4) = xc
462  xpq( 5) = xbot
463  xpq(-5) = xbot
464 
465 C polarized part
466 
467  va0 = 0.225
468  va10 = 0.15
469  valphau = 0.326
470  valphad = 0.505
471 
472  ffu0 = 1./(1.+va0*x**(-valphau)*(1.-x)*(1.-x))
473  ffu1 = 1./(1.+va0*va10*x**(-valphau)*(1.-x)*(1.-x))
474  ffd0 = 1./(1.+va0*x**(-valphad)*(1.-x)*(1.-x))
475  ffd1 = 1./(1.+va0*va10*x**(-valphad)*(1.-x)*(1.-x))
476 
477  xduv = xuv*ffu0
478  & -xdv*(ffu0+ffu1)/3.
479  xddv = xdv*(ffd0-2*ffd1)/3.
480 
481  xdub = 0.
482  xddb = 0.
483  xdg = xuv*(1.-ffu0)/2.
484  & +xdv*(2*ffd1-ffd0+ffu1+ffu0-3.)/6.
485 
486  xdpq( 0) = xdg
487  xdpq( 1) = xddv + xddb
488  xdpq(-1) = xddb
489  xdpq( 2) = xduv + xdub
490  xdpq(-2) = xdub
491 
492  ELSEIF((lst(15).EQ.105).OR.(lst(15).EQ.106)) THEN
493 C***************************************************************************
494 C
495 C set 105 Bartelski et al hep-ph/9502271 Set II(p,n)
496 C set 106 Bartelski et al hep-ph/9502271 Set II(P,n)
497 C
498 C***************************************************************************
499 
500  xuv = 1.996*x**(-0.462)*(1.-x)**3.96
501  & *(1. - 0.39*sqrt(x) + 5.13*x)*x
502  xdv = 0.296*x**(-0.67)*(1.-x)**4.71
503  & *(1. + 5.03*sqrt(x) + 5.56*x)*x
504  fdelta = 0.099*x**(-0.462)*(1.-x)**9.27*(1.+25*x)
505  fm = 0.411*x**(-1.3 )*(1.-x)**9.27
506  & *(1. - 1.15*sqrt(x) + 15.6*x)
507  xub = 0.5*x*(0.392*fm - fdelta)
508  xdb = 0.5*x*(0.392*fm + fdelta)
509  xs = 0.5*x*0.196*fm
510  xc = 0.5*x*0.020*fm
511  xbot = 0.
512  xg = 0.775*x**(-1.3)*(1.-x)**5.3*(1.+5.2*x)*x
513 
514  xpq( 0) = xg
515  xpq( 1) = xdv + xdb
516  xpq(-1) = xdb
517  xpq( 2) = xuv + xub
518  xpq(-2) = xub
519  xpq( 3) = xs
520  xpq(-3) = xs
521  xpq( 4) = xc
522  xpq(-4) = xc
523 
524 C polarized part
525 
526  IF (lst(15).EQ.105) THEN
527  va1 = 0.924
528  va2 = -3.237
529  va3 = 11.3
530  xduv = x**(-0.462)*(1.-x)**3.96*(va1+va2*sqrt(x)+va3*x)*x
531 
532  va1 = -0.029
533  va2 = -0.163
534  va3 = -1.644
535  xddv = x**(-0.670)*(1.-x)**4.71*(va1+va2*sqrt(x)+va3*x)*x
536 
537  vc1 = 0.
538  vc2 = -0.993
539  vd = -0.015
540  fdelta = x**(-0.462)*(1.-x)**9.27*vd*(1.+25.*x)
541  fm = x**(-0.8 )*(1.-x)**9.27
542  & *(vc1 + vc2*sqrt(x))
543 
544  xdub = 0.5*x*(0.392*fm - fdelta)
545  xddb = 0.5*x*(0.392*fm + fdelta)
546  xds = 0.5*x*0.196*fm
547  xdc = 0.5*x*0.020*fm
548 
549  ELSEIF(lst(15).EQ.106) THEN
550 
551  va1 = 0.929
552  va2 = -3.005
553  va3 = 10.24
554  xduv = x**(-0.462)*(1.-x)**3.96*(va1+va2*sqrt(x)+va3*x)*x
555 
556  va1 = -0.066
557  va2 = -0.064
558  va3 = -1.644
559  xddv = x**(-0.670)*(1.-x)**4.71*(va1+va2*sqrt(x)+va3*x)*x
560 
561  vc1 = 0.
562  vc2 = -0.861
563  vd = -0.013
564  fdelta = x**(-0.462)*(1.-x)**9.27*vd*(1.+25.*x)
565  fm = x**(-0.8 )*(1.-x)**9.27
566  & *(vc1+ vc2*sqrt(x))
567  xdub = 0.5*x*(0.392*fm - fdelta)
568  xddb = 0.5*x*(0.392*fm + fdelta)
569  xds = 0.5*x*0.196*fm
570  xdc = 0.5*x*0.020*fm
571  ENDIF
572 
573  xdpq( 1) = xddv + xddb
574  xdpq(-1) = xddb
575  xdpq( 2) = xduv + xdub
576  xdpq(-2) = xdub
577  xdpq( 3) = xds
578  xdpq(-3) = xds
579  xdpq( 4) = xdc
580  xdpq(-4) = xdc
581  ELSEIF((lst(15).GE.107).AND.(lst(15).LE.109)) THEN
582 C ***************************************************************************
583 C
584 C LST(15) = 107 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
585 C LST(15) = 108 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
586 C LST(15) = 109 Gehrmann hep-ph/9512406 Gluon A (NLO) + DGLAP
587 C
588 C ***************************************************************************
589 
590  IF (lst(15).EQ.107) iflag = 0
591  IF (lst(15).EQ.108) iflag = 1
592  IF (lst(15).EQ.109) iflag = 2
593 
594  mode = 20 ! MSR(A prime) Parametrization
595 
596 
597 C unpolarized part
598 * MRS(A prime)
599 
600  dx=dble(x)
601  dq2=dble(q2)
602 
603  CALL mrseb(dx,dq2,mode,mrsupv,mrsdnv,mrsusea,
604  & mrsdsea,mrsstr,mrschm,mrsbot,mrsglu)
605 
606 
607  xpq( 0) = sngl(mrsglu)
608  xpq( 1) = sngl(mrsdnv + mrsdsea)
609  xpq(-1) = sngl( mrsdsea)
610  xpq( 2) = sngl(mrsupv + mrsusea)
611  xpq(-2) = sngl( mrsusea)
612  xpq( 3) = sngl(mrsstr)
613  xpq(-3) = sngl(mrsstr)
614  xpq( 4) = sngl(mrschm)
615  xpq(-4) = sngl(mrschm)
616  xpq( 5) = sngl(mrsbot)
617  xpq(-5) = sngl(mrsbot)
618 
619 C polarized part (evolved set of Gehrmann hep-ph/9512406 Gluon A,B,C (LO))
620  CALL polnlo(iflag,x,q2,dxduv,dxddv,dxdg,dxdubar,dxddbar,dxds)
621 
622  xdpq( 0) = sngl(dxdg)
623  xdpq( 1) = sngl(dxddv+dxddbar)
624  xdpq(-1) = sngl(dxddbar)
625  xdpq( 2) = sngl(dxduv+dxdubar)
626  xdpq(-2) = sngl(dxdubar)
627  xdpq( 3) = sngl(dxds)
628  xdpq(-3) = sngl(dxds)
629 
630 C ***************************************************************************
631 C
632 C LST(15) = 110 Gehrmann et al hep-ph/9512406 Gluon A (LO)
633 C LST(15) = 111 Gehrmann et al hep-ph/9512406 Gluon B (LO)
634 C LST(15) = 112 Gehrmann et al hep-ph/9512406 Gluon C (LO)
635 C
636 C LST(15) = 130 Fake sample :
637 C unpolarized Gehrmann et al hep-ph/9512406
638 C with Delta u(x) = 0.5 * u(x) and Delta d(x) = 0.
639 C LST(15) = 131 Fake sample :
640 C unpolarized Gehrmann et al hep-ph/9512406
641 C with Delta d(x) = 0.5 * d(x) and Delta u(x) = 0.
642 C
643 C ***************************************************************************
644 
645 
646 
647  ELSEIF((lst(15).GE.110).AND.(lst(15).LE.112) .or.
648  & (lst(15).GE.130).AND.(lst(15).LE.131) ) THEN
649 
650 C unpolarized part
651 
652  xuv = 3.221*x**0.564*(1.-x)**3.726
653  & *(1.-0.6889*x**0.200 + 2.254*x + 1.261*sqrt(x)*x)
654  xdv = 0.507*x**0.376*(1.-x)**4.476
655  & *(1.+1.615 *x**0.553 + 3.651*x + 1.299*sqrt(x)*x)
656  xupdb = (x**0.158*(0.738 - 0.981*x + 1.063*x*x)
657  & * (-log(x))**0.037
658  & + 0.00285*exp(sqrt(-4.010*log(x))))*(1.-x)**6.356
659  xdmub = 0.1067 * x**0.4036 * (1.+ 0.4680*x**0.8763 +
660  & x * 18.6760) * (1.- x)**8.6220
661  xs = 0.0034*(-log(x))**(-1.15)
662  & *(1.-2.392*sqrt(x) + 7.094*x)*(1.-x)**6.166
663  & *exp(sqrt(-6.719*log(x)))
664  xg = (x**0.731*(5.110 - 1.204*x - 1.911*x*x)
665  & * (-log(x))**(-0.4718)
666  & + 0.0527*exp(sqrt(-4.584*log(x))))*(1.-x)**5.566
667 
668  xpq( 0) = xg
669  xpq( 1) = xdv + (xupdb + xdmub)*0.5
670  xpq(-1) = (xupdb + xdmub)*0.5
671  xpq( 2) = xuv + (xupdb - xdmub)*0.5
672  xpq(-2) = (xupdb - xdmub)*0.5
673  xpq( 3) = xs
674  xpq(-3) = xs
675 
676 C polarized part
677 
678 C Gluon A
679  IF(lst(15).EQ.110) THEN
680  vetau = 0.823
681  vetad = -0.303
682  vetaqb = -0.0495
683  vetag = 1.9
684 
685  vau = 0.578
686  vad = 0.666
687  vaqb = 0.520
688  vag = 0.520
689 
690  vbu = 3.73
691  vbd = 4.73
692  vbqb = 15.06
693  vbg = 9.45
694 
695  vgu = 9.38
696  vgd = 10.46
697  vgqb = 2.30
698  vgg = 0.0
699 
700  vru = -4.26
701  vrd = -5.10
702  vrqb = -2.00
703  vrg = 0.00
704 
705  vnormu = (1.+ vgu*vau/(vau+vbu+1.))
706  & * gamma(vau)*gamma(vbu+1.)/gamma(vau+vbu+1.)
707  & + vru*gamma(vau+0.5)*gamma(vbu+1.)
708  & /gamma(vau+vbu+1.5)
709  vnormd = (1.+ vgd*vad/(vad+vbd+1.))
710  & * gamma(vad)*gamma(vbd+1.)/gamma(vad+vbd+1.)
711  & + vrd*gamma(vad+0.5)*gamma(vbd+1.)
712  & /gamma(vad+vbd+1.5)
713  vnormqb = (1.+ vgqb*vaqb/(vaqb+vbqb+1.))
714  & * gamma(vaqb)*gamma(vbqb+1.)/gamma(vaqb+vbqb+1.)
715  & + vrqb*gamma(vaqb+0.5)*gamma(vbqb+1.)
716  & /gamma(vaqb+vbqb+1.5)
717  vnormg = (1.+ vgg*vag/(vag+vbg+1.))
718  & * gamma(vag)*gamma(vbg+1.)/gamma(vag+vbg+1.)
719  & + vrg*gamma(vag+0.5)*gamma(vbg+1.)
720  & /gamma(vag+vbg+1.5)
721  vnormu = 1./vnormu
722  vnormd = 1./vnormd
723  vnormqb = 1./vnormqb
724  vnormg = 1./vnormg
725 
726  xduv = vetau*vnormu*x**vau*(1.-x)**vbu*(1.+vgu*x
727  & + vru*sqrt(x))
728  xddv = vetad*vnormd*x**vad*(1.-x)**vbd*(1.+vgd*x
729  & + vrd*sqrt(x))
730  xdub = vetaqb*vnormqb*x**vaqb*(1.-x)**vbqb*
731  & (1.+vgqb*x + vrqb*sqrt(x))
732  xddb = xdub
733  xds = xdub
734  xdg = vetag*vnormg*x**vag*(1.-x)**vbg*(1.+vgg*x
735  & + vrg*sqrt(x))
736 
737  xdpq( 0) = xdg
738  xdpq( 1) = xddv + xddb
739  xdpq(-1) = xddb
740  xdpq( 2) = xduv + xdub
741  xdpq(-2) = xdub
742  xdpq( 3) = xds
743  xdpq(-3) = xds
744 
745 C Gluon B
746  ELSE IF(lst(15).EQ.111) THEN
747 
748  vetau = 0.823
749  vetad = -0.303
750  vetaqb = -0.0495
751  vetag = 1.9
752 
753  vau = 0.585
754  vad = 0.662
755  vaqb = 0.524
756  vag = 0.524
757 
758  vbu = 3.73
759  vbd = 4.73
760  vbqb = 15.96
761  vbg = 6.87
762 
763  vgu = 9.31
764  vgd = 10.91
765  vgqb = 2.42
766  vgg = 1.
767 
768  vru = -4.28
769  vrd = -5.09
770  vrqb = -2.00
771  vrg = -2.00
772 
773  vnormu = (1.+ vgu*vau/(vau+vbu+1.))
774  & * gamma(vau)*gamma(vbu+1.)/gamma(vau+vbu+1.)
775  & + vru*gamma(vau+0.5)*gamma(vbu+1.)
776  & /gamma(vau+vbu+1.5)
777  vnormd = (1.+ vgd*vad/(vad+vbd+1.))
778  & * gamma(vad)*gamma(vbd+1.)/gamma(vad+vbd+1.)
779  & + vrd*gamma(vad+0.5)*gamma(vbd+1.)
780  & /gamma(vad+vbd+1.5)
781  vnormqb = (1.+ vgqb*vaqb/(vaqb+vbqb+1.))
782  & * gamma(vaqb)*gamma(vbqb+1.)/gamma(vaqb+vbqb+1.)
783  & + vrqb*gamma(vaqb+0.5)*gamma(vbqb+1.)
784  & /gamma(vaqb+vbqb+1.5)
785  vnormg = (1.+ vgg*vag/(vag+vbg+1.))
786  & * gamma(vag)*gamma(vbg+1.)/gamma(vag+vbg+1.)
787  & + vrg*gamma(vag+0.5)*gamma(vbg+1.)
788  & /gamma(vag+vbg+1.5)
789 
790  vnormu = 1./vnormu
791  vnormd = 1./vnormd
792  vnormqb = 1./vnormqb
793  vnormg = 1./vnormg
794 
795  xduv = vetau*vnormu*x**vau*(1.-x)**vbu
796  & *(1.+vgu*x + vru*sqrt(x))
797  xddv = vetad*vnormd*x**vad*(1.-x)**vbd
798  & *(1.+vgd*x + vrd*sqrt(x))
799  xdub = vetaqb*vnormqb*x**vaqb*(1.-x)**vbqb*
800  & (1.+vgqb*x + vrqb*sqrt(x))
801  xddb = xdub
802  xds = xdub
803  xdg = vetag*vnormg*x**vag*(1.-x)**vbg
804  & *(1.+vgg*x + vrg*sqrt(x))
805 
806  xdpq( 0) = xdg
807  xdpq( 1) = xddv + xddb
808  xdpq(-1) = xddb
809  xdpq( 2) = xduv + xdub
810  xdpq(-2) = xdub
811  xdpq( 3) = xds
812  xdpq(-3) = xds
813 
814 C Gluon C
815  ELSEIF(lst(15).EQ.112) THEN
816 
817  vetau = 0. 823
818  vetad = -0.303
819  vetaqb = -0.0495
820  vetag = 1.9
821 
822  vau = 0.582
823  vad = 0.660
824  vaqb = 0.456
825  vag = 0.456
826 
827  vbu = 3.73
828  vbd = 4.73
829  vbqb = 11.82
830  vbg = 8.72
831 
832  vgu = 9.50
833  vgd = 11.04
834  vgqb = 2.11
835  vgg = 0.0
836 
837  vru = -4.28
838  vrd = -5.06
839  vrqb = -1.95
840  vrg = -3.00
841 
842  vnormu = (1.+ vgu*vau/(vau+vbu+1.))
843  & * gamma(vau)*gamma(vbu+1.)/gamma(vau+vbu+1.)
844  & + vru*gamma(vau+0.5)*gamma(vbu+1.)
845  & /gamma(vau+vbu+1.5)
846  vnormd = (1.+ vgd*vad/(vad+vbd+1.))
847  & * gamma(vad)*gamma(vbd+1.)/gamma(vad+vbd+1.)
848  & + vrd*gamma(vad+0.5)*gamma(vbd+1.)
849  & /gamma(vad+vbd+1.5)
850  vnormqb = (1.+ vgqb*vaqb/(vaqb+vbqb+1.))
851  & * gamma(vaqb)*gamma(vbqb+1.)/gamma(vaqb+vbqb+1.)
852  & + vrqb*gamma(vaqb+0.5)*gamma(vbqb+1.)
853  & /gamma(vaqb+vbqb+1.5)
854  vnormg = (1.+ vgg*vag/(vag+vbg+1.))
855  & * gamma(vag)*gamma(vbg+1.)/gamma(vag+vbg+1.)
856  & + vrg*gamma(vag+0.5)*gamma(vbg+1.)
857  & /gamma(vag+vbg+1.5)
858 
859  vnormu = 1./vnormu
860  vnormd = 1./vnormd
861  vnormqb = 1./vnormqb
862  vnormg = 1./vnormg
863 
864  xduv = vetau*vnormu*x**vau*(1.-x)**vbu
865  & *(1.+vgu*x + vru*sqrt(x))
866  xddv = vetad*vnormd*x**vad*(1.-x)**vbd
867  & *(1.+vgd*x + vrd*sqrt(x))
868  xdub = vetaqb*vnormqb*x**vaqb*(1.-x)**vbqb*
869  & (1.+vgqb*x + vrqb*sqrt(x))
870  xddb = xdub
871  xds = xdub
872  xdg = vetag*vnormg*x**vag*(1.-x)**vbg
873  & *(1.+vgg*x + vrg*sqrt(x))
874 
875  xdpq( 0) = xdg
876  xdpq( 1) = xddv + xddb
877  xdpq(-1) = xddb
878  xdpq( 2) = xduv + xdub
879  xdpq(-2) = xdub
880  xdpq( 3) = xds
881  xdpq(-3) = xds
882 
883  ELSEIF(lst(15).EQ.130) THEN
884 C unpolarized Gehrmann et al hep-ph/9512406
885 C with Delta u(x) = 0.5 * u(x) and Delta d(x) = 0.
886  xdpq(2) = 0.5 * xpq(2)
887 
888  ELSEIF(lst(15).EQ.131) THEN
889 C unpolarized Gehrmann et al hep-ph/9512406
890 C with Delta d(x) = 0.5 * d(x) and Delta u(x) = 0.
891  xdpq(1) = 0.5 * xpq(1)
892 
893  ENDIF
894 
895 C***************************************************************************
896 C
897 C LST(15) = 113 Gehrmann et al hep-ph/9512406 Gluon A (LO) + (DGLAP)
898 C LST(15) = 114 Gehrmann et al hep-ph/9512406 Gluon B (LO) + (DGLAP)
899 C LST(15) = 115 Gehrmann et al hep-ph/9512406 Gluon C (LO) + (DGLAP)
900 C
901 C**************************************************************************
902 
903  ELSEIF ((lst(15).GE.113).AND.(lst(15).LE.115)) THEN
904 
905  IF (lst(15).EQ.113) iflag = 0
906  IF (lst(15).EQ.114) iflag = 1
907  IF (lst(15).EQ.115) iflag = 2
908 
909 C unpolarized part
910 * M. GLUECK, E.REYA, A.VOGT : *
911 * DO-TH 94/24 = DESY 94-206 *
912 * (TO APPEAR IN Z. PHYS. C) *
913 
914  dx=dble(x)
915  dq2=dble(q2)
916  CALL grv94lo(dx,dq2
917  & ,ddxuv,ddxdv,ddxdel,ddxudb,ddxsb,ddxgl)
918 
919  xpq( 0) = sngl(ddxgl)
920  xpq( 1) = sngl(ddxdv + (ddxudb + ddxdel)*0.5d0)
921  xpq(-1) = sngl( (ddxudb + ddxdel)*0.5d0)
922  xpq( 2) = sngl(ddxuv + (ddxudb - ddxdel)*0.5d0)
923  xpq(-2) = sngl( (ddxudb - ddxdel)*0.5d0)
924  xpq( 3) = sngl(ddxsb)
925  xpq(-3) = sngl(ddxsb)
926 
927 C polarized part (evolved set of Gehrmann hep-ph/9512406 Gluon A,B,C (LO))
928  CALL polpar(iflag,dx,dq2
929  & ,dxduv,dxddv,dxdg,dxdqbar,dxds)
930 
931  xdpq( 0) = sngl(dxdg)
932  xdpq( 1) = sngl(dxddv + dxdqbar)
933  xdpq(-1) = sngl(dxdqbar)
934  xdpq( 2) = sngl(dxduv + dxdqbar)
935  xdpq(-2) = sngl(dxdqbar)
936  xdpq( 3) = sngl(dxds)
937  xdpq(-3) = sngl(dxds)
938 
939 
940 C ***************************************************************************
941 C
942 C M. Glueck, E. Reya, M. Stratmann and W. Vogelsang,
943 C DO-TH 95/13, RAL-TR-95-042
944 C
945 C LST(15) = 116 'standard' scenario, next-to-leading order
946 C LST(15) = 117 'valence' scenario, next-to-leading order
947 C LST(15) = 118 'standard' scenario, leading order
948 C LST(15) = 119 'valence' scenario, leading order
949 C
950 C ***************************************************************************
951 
952  ELSEIF ((lst(15).GE.116).AND.(lst(15).LE.119)) THEN
953  IF (lst(15).EQ.116) iset = 1
954  IF (lst(15).EQ.117) iset = 2
955  IF (lst(15).EQ.118) iset = 3
956  IF (lst(15).EQ.119) iset = 4
957 
958 C unpolarized part
959 * M. GLUECK, E.REYA, A.VOGT : *
960 * DO-TH 94/24 = DESY 94-206 *
961 * (TO APPEAR IN Z. PHYS. C) *
962 
963 
964  dx=dble(x)
965  dq2=dble(q2)
966  CALL grv94lo(dx,dq2
967  & ,ddxuv,ddxdv,ddxdel,ddxudb,ddxsb,ddxgl)
968 
969  xpq( 0) = sngl(ddxgl)
970  xpq( 1) = sngl(ddxdv + (ddxudb + ddxdel)*0.5d0)
971  xpq(-1) = sngl( (ddxudb + ddxdel)*0.5d0)
972  xpq( 2) = sngl(ddxuv + (ddxudb - ddxdel)*0.5d0)
973  xpq(-2) = sngl( (ddxudb - ddxdel)*0.5d0)
974  xpq( 3) = sngl(ddxsb)
975  xpq(-3) = sngl(ddxsb)
976 
977 C Polarisierter Teil
978 
979  CALL parpol(iset,dx,dq2,
980  & dxduv, dxddv, dxdqbar, dxds, dxdg)
981 
982  xdpq( 0) = sngl(dxdg)
983  xdpq( 1) = sngl(dxddv) + sngl(dxdqbar)
984  xdpq(-1) = sngl(dxdqbar)
985  xdpq( 2) = sngl(dxduv) + sngl(dxdqbar)
986  xdpq(-2) = sngl(dxdqbar)
987  xdpq( 3) = sngl(dxds)
988  xdpq(-3) = sngl(dxds)
989 
990 C*************************************************************************
991 
992  ELSEIF(lst(15).EQ.120)THEN
993 
994 * This is polarization distribution of Stanley J.Brodsky et al.
995 * Nucl.Phys. B441(1995)
996  y=1.-x
997  x1=x**(-1.12)
998  DO i = 1,4
999  qp(i) = x1*(aq(i)*y**st1(i)+bq(i)*y**st2(i))
1000  qm(i) = x1*(cq(i)*y**st3(i)+dq(i)*y**st4(i))
1001  qg(i) = qp(i)+qm(i)
1002  qj(i) = qp(i)-qm(i)
1003  ENDDO
1004  DO i = 1,2
1005  qg(i) = qg(i)-2.*qg(3)
1006  qj(i) = qj(i)-2.*qj(3)
1007  ENDDO
1008  DO i = 1,4
1009  xpq(i-4+int(i/4.)*3) = qg(3)
1010  xdpq(i-4+int(i/4.)*3) = qj(3)
1011  ENDDO
1012 
1013  xpq(0) = qg(4)
1014  xpq(1) = qg(1)
1015  xpq(2) = qg(2)
1016  xdpq(0) = qj(4)
1017  xdpq(1) = qj(1)
1018  xdpq(2) = qj(2)
1019 
1020  ELSEIF(lst(15).EQ.121)THEN
1021 
1022 * This is polarization distribution of S.Keler & J.F.Owens
1023 * Phys.Lett. B266(1991)126 and Phys.Rev. D19(1994)1199
1024 
1025  IF(q2.LE.4.) q2 = 4.
1026  z=alog(alog(q2/ww(2))/ww(3))
1027 
1028  DO i = 1,6
1029  DO j = 1,4
1030  aaa(i,j) = const0(i,j) + const1(i,j)*z +
1031  & const2(i,j)*z*z + const3(i,j)*z*z*z
1032  ENDDO
1033  ENDDO
1034  DO i = 1,2
1035  sumr(i) = 0.
1036  DO 1 j = 1,3
1037  DO 1 k = 1,50
1038  x1=xw(j)+k/50.
1039  ss=1.+aaa(3,i)*x1+aaa(4,i)*x1*x1
1040  sumr(i)=sumr(i)+x1**aaa(1,i)*(1.-x1)*
1041  & *aaa(2,i)*ss*w(j)
1042  1 CONTINUE
1043  aaa(6,i)=(5.-2.*i)/sumr(i)
1044  ENDDO
1045  DO i = 1,4
1046  ss = 1.+aaa(3,i)*x+aaa(4,i)*x*x+aaa(5,i)*x*x*x
1047  qg(i) = x**aaa(1,i)*(1.-x)**aaa(2,i)*ss*aaa(6,i)
1048  ENDDO
1049  ss=qg(2)
1050  qg(3)=qg(3)/6.
1051  qg(2)=qg(1)-ss
1052  qg(1)=ss
1053  DO i = 1,4
1054  ss = 1. + cj(i)*x
1055  qj(i) = x**aj(i)*(1.-x)**bj(i)*ss*rj(i)
1056  ENDDO
1057 
1058  DO i = 1,4
1059  xpq(i-4+int(i/4.)*3) = qg(3)
1060  xdpq(i-4+int(i/4.)*3) = qj(3)
1061  ENDDO
1062  xpq(0) = qg(4)
1063  xpq(1) = qg(1)
1064  xpq(2) = qg(2)
1065  xdpq(0) = qj(4)
1066  xdpq(1) = qj(1)
1067  xdpq(2) = qj(2)
1068 C***********************************************************************
1069 C LST(15) = 124 D. de Florian et al., hep-ph/9711440 LO set 1
1070 C LST(15) = 125 D. de Florian et al., hep-ph/9711440 LO set 2
1071 C LST(15) = 126 D. de Florian et al., hep-ph/9711440 LO set 3
1072 C LST(15) = 127 D. de Florian et al., hep-ph/9711440 NLO set 1
1073 C LST(15) = 128 D. de Florian et al., hep-ph/9711440 NLO set 2
1074 C LST(15) = 129 D. de Florian et al., hep-ph/9711440 NLO set 3
1075 C
1076 C using GRV 1995 as unpolarized reference set
1077 C
1078 C***********************************************************************
1079  ELSEIF ((lst(15).GE.124).AND.(lst(15).LE.129)) THEN
1080 
1081 c IF(X.LE.0.0001) THEN
1082 c WRITE(*,*) 'Warning: x exceeds minum value, x set to 0.0001001'
1083 c X=0.0001001
1084 c ENDIF
1085 c IF(X.GE.0.9) THEN
1086 c WRITE(*,*) 'Warning: x exceeds maxinum value, x set to 0.89999'
1087 c X=0.89999
1088 c ENDIF
1089 c IF(Q2.LE.1.) THEN
1090 c WRITE(*,*)'Warning: Q2 exceeds minimum value, Q2 set to 1.0001'
1091 c Q2 = 1.0001
1092 c ENDIF
1093 c IF(Q2.GE.50000.0) THEN
1094 c WRITE(*,*) 'Warning: Q2 exceeds maximum value, Q2 set to 49999.9'
1095 c Q2 = 49999.9
1096 c ENDIF
1097 
1098  IF (lst(15).EQ.124) mode=1
1099  IF (lst(15).EQ.125) mode=2
1100  IF (lst(15).EQ.126) mode=3
1101  IF (lst(15).EQ.127) mode=4
1102  IF (lst(15).EQ.128) mode=5
1103  IF (lst(15).EQ.129) mode=6
1104 
1105 C unpolarized part
1106 * M. GLUECK, E.REYA, A.VOGT : *
1107 * DO-TH 94/24 = DESY 94-206 *
1108 * (TO APPEAR IN Z. PHYS. C) *
1109 
1110  dx=dble(x)
1111  dq2=dble(q2)
1112  IF ((lst(15).GE.124).AND.(lst(15).LE.126)) THEN
1113  CALL grv94lo(dx,dq2
1114  & ,ddxuv,ddxdv,ddxdel,ddxudb,ddxsb,ddxgl)
1115  ELSEIF ((lst(15).GE.127).AND.(lst(15).LE.129)) THEN
1116  CALL grv94ho(dx,dq2
1117  & ,ddxuv,ddxdv,ddxdel,ddxudb,ddxsb,ddxgl)
1118  ENDIF
1119 
1120  xpq( 0) = sngl(ddxgl)
1121  xpq( 1) = sngl(ddxdv + (ddxudb + ddxdel)*0.5d0)
1122  xpq(-1) = sngl( (ddxudb + ddxdel)*0.5d0)
1123  xpq( 2) = sngl(ddxuv + (ddxudb - ddxdel)*0.5d0)
1124  xpq(-2) = sngl( (ddxudb - ddxdel)*0.5d0)
1125  xpq( 3) = sngl(ddxsb)
1126  xpq(-3) = sngl(ddxsb)
1127 
1128 * polarized part
1129 
1130  CALL polfit(mode,dx,dq2,dxduv,dxddv,dxdqbar,dxdstr,dxdglu)
1131  xdpq( 0) = sngl(dxdglu)
1132  xdpq( 1) = sngl(dxddv) + sngl(dxdqbar)
1133  xdpq(-1) = sngl(dxdqbar)
1134  xdpq( 2) = sngl(dxduv) + sngl(dxdqbar)
1135  xdpq(-2) = sngl(dxdqbar)
1136  xdpq( 3) = sngl(dxdstr)
1137  xdpq(-3) = sngl(dxdstr)
1138  ELSEIF(lst(15).EQ.132)THEN
1139 *72*********************************************************************
1140  CALL fitparton(x,q2,xpq,xdpq)
1141 *72*********************************************************************
1142  ELSEIF ((lst(15).GE.133).AND.(lst(15).LE.136)) THEN
1143 
1144 C unpolarized part
1145 * "IMPROVED PARTON DISTRIBUTIONS FROM GLOBAL ANALYSIS OF
1146 * RECENT DEEP INELASTIC SCATTERING AND INCLUSIVE JET DATA"
1147 * By: H.L. Lai, J. Huston, S. Kuhlmann, F. Olness, J. Owens,
1148 * D. Soper W.K. Tung, H. Weerts
1149 * Phys. Rev. D55, 1280 (1997)
1150 * ( low Q2 set )
1151 *
1152 
1153 * from data statement above
1154 * CTEQSET=10 ! Low Q2
1155  dx=dble(x)
1156  dq2=dble(q2)
1157  IF(initcteq.eq.1) THEN
1158  initcteq=0
1159  CALL setctq4(cteqset)
1160  ENDIF
1161 
1162 * Note in CTEQ u and d indices are reversed WRT pepsi order:
1163 * Iparton is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5)
1164 * for (b, c, s, d, u, g, u_bar, ..., b_bar),
1165 
1166 
1167  ddxpq(-3) = ctq4pdf(-3,dx,dq2)
1168  ddxpq(-1) = ctq4pdf(-2,dx,dq2)
1169  ddxpq(-2) = ctq4pdf(-1,dx,dq2)
1170  ddxpq(0) = ctq4pdf(0,dx,dq2)
1171  ddxpq(2) = ctq4pdf(1,dx,dq2)
1172  ddxpq(1) = ctq4pdf(2,dx,dq2)
1173  ddxpq(3) = ctq4pdf(3,dx,dq2)
1174 
1175 * CTEQ routine double precision q(x). Make single xq(x)
1176  DO i=-3,3
1177  xpq( i) = x*sngl(ddxpq(i))
1178  ENDDO
1179 
1180  xdpq( 0) = 0.0
1181  xdpq( 1) = -0.3*xpq(1)
1182  xdpq(-1) = 0.0
1183  xdpq( 2) = 0.5*xpq(2)
1184  xdpq(-2) = 0.0
1185  xdpq( 3) = 0.0
1186  xdpq(-3) = 0.0
1187 
1188 * require F2 non-zero
1189  IF(xpq(2).EQ.0.0) xpq(2)=0.000001
1190 
1191  ELSEIF((lst(15).GE.137).AND.(lst(15).LE.138)) THEN
1192 
1193  mode = 10 ! MSR(A prime) Low Q2 Parametrization
1194 
1195 
1196 C unpolarized part
1197 * MRS(A prime)
1198 
1199  dx=dble(x)
1200  dq2=dble(q2)
1201 
1202  CALL mrseb(dx,dq2,mode,mrsupv,mrsdnv,mrsusea,
1203  & mrsdsea,mrsstr,mrschm,mrsbot,mrsglu)
1204 
1205 
1206  xpq( 0) = sngl(mrsglu)
1207  xpq( 1) = sngl(mrsdnv + mrsdsea)
1208  xpq(-1) = sngl( mrsdsea)
1209  xpq( 2) = sngl(mrsupv + mrsusea)
1210  xpq(-2) = sngl( mrsusea)
1211  xpq( 3) = sngl(mrsstr)
1212  xpq(-3) = sngl(mrsstr)
1213  xpq( 4) = sngl(mrschm)
1214  xpq(-4) = sngl(mrschm)
1215  xpq( 5) = sngl(mrsbot)
1216  xpq(-5) = sngl(mrsbot)
1217 
1218 
1219 
1220  xdpq( 0) = 0.0
1221  xdpq( 1) = -0.3*xpq(1)
1222  xdpq(-1) = 0.0
1223  xdpq( 2) = 0.5*xpq(2)
1224  xdpq(-2) = 0.0
1225  xdpq( 3) = 0.0
1226  xdpq(-3) = 0.0
1227 
1228 
1229 
1230 C***********************************************************************
1231 C LST(15) = 144 grsv2000 hep-ph/0011215 LO standard scenario
1232 C LST(15) = 145 grsv2000 hep-ph/0011215 LO valence scenario
1233 C LST(15) = 146 grsv2000 hep-ph/0011215 NLO standard scenario
1234 C LST(15) = 147 grsv2000 hep-ph/0011215 NLO valence scenario
1235 
1236 C
1237 C using GRV 1998 as unpolarized reference set
1238 C
1239 C***********************************************************************
1240  ELSEIF ((lst(15).GE.144).AND.(lst(15).LE.147)) THEN
1241 
1242  IF(x.LE.0.00001) THEN
1243 c$$$ WRITE(*,*) 'Warning: x exceeds minum value, x set to 0.0001001'
1244  x=0.0001001
1245  ENDIF
1246  IF(x.GE.1.0) THEN
1247 c$$$ WRITE(*,*) 'Warning: x exceeds maxinum value, x set to 0.99999'
1248  x=0.999999
1249  ENDIF
1250  IF(q2.LE.0.8) THEN
1251 c$$$ WRITE(*,*)'Warning: Q2 exceeds minimum value, Q2 set to 0.80001'
1252  q2 = 0.80001
1253  ENDIF
1254  IF(q2.GE.1.0e6) THEN
1255 c$$$ WRITE(*,*) 'Warning: Q2 exceeds maximum value, Q2 set to 999999.9'
1256  q2 = 999999.9
1257  ENDIF
1258 
1259 
1260 C unpolarized part
1261 * M. Glueck, E. Reya, A. Vogt : *
1262 * hep-ph/9806404 = DO-TH 98/07 = WUE-ITP-98-019 *
1263 * (To appear in Eur. Phys. J. C) *
1264 
1265  dx=dble(x)
1266  dq2=dble(q2)
1267  IF ((lst(15).GE.144).AND.(lst(15).LE.145)) THEN
1268  iset = 1 ! LO
1269  ELSEIF ((lst(15).GE.146).AND.(lst(15).LE.147)) THEN
1270  iset = 2 ! NLO
1271  ENDIF
1272  CALL grv98pa(iset, dx, dq2
1273  & ,ddxuv,ddxdv,ddxus, ddxds, ddxsb,ddxgl)
1274 
1275 
1276  xpq( 0) = sngl(ddxgl)
1277  xpq( 1) = sngl(ddxdv + ddxds)
1278  xpq(-1) = sngl(ddxds)
1279  xpq( 2) = sngl(ddxuv + ddxus)
1280  xpq(-2) = sngl(ddxus)
1281  xpq( 3) = sngl(ddxsb)
1282  xpq(-3) = sngl(ddxsb)
1283 
1284 * polarized part
1285 
1286 
1287  IF (lst(15).EQ.144) iset = 3
1288  IF (lst(15).EQ.145) iset = 4
1289  IF (lst(15).EQ.146) iset = 1
1290  IF (lst(15).EQ.147) iset = 2
1291 
1292 
1293 
1294  CALL parpolnew(iset, dx, dq2,
1295  & dxdu, dxdd, dxdub, dxddb, dxdstr, dxdglu,
1296  & dummyg1p, dummyg1n)
1297  xdpq( 0) = sngl(dxdglu)
1298  xdpq( 1) = sngl(dxdd)
1299  xdpq(-1) = sngl(dxddb)
1300  xdpq( 2) = sngl(dxdu)
1301  xdpq(-2) = sngl(dxdub)
1302  xdpq( 3) = sngl(dxdstr)
1303  xdpq(-3) = sngl(dxdstr)
1304 
1305 C*************************************************************************
1306 C LST(15) = 150 cteq5l LO
1307 C LST(15) = 151 cteq5m NLO MSBAR
1308 C LST(15) = 152 cteq5m1 NLO MSBAR (update)
1309 C*************************************************************************
1310 
1311  ELSEIF((lst(15).GE.150).AND.(lst(15).LE.152)) THEN
1312 
1313  IF(x.LE.0.00001) THEN
1314 c$$$ WRITE(*,*) 'Warning: x exceeds minum value, x set to 0.0000101'
1315  x=0.0000101
1316  ENDIF
1317  IF(x.GE.1.0) THEN
1318 c$$$ WRITE(*,*) 'Warning: x exceeds maxinum value, x set to 0.99999'
1319  x=0.999999
1320  ENDIF
1321  IF(q2.LE.1.0) THEN
1322 c$$$ WRITE(*,*)'Warning: Q2 exceeds minimum value, Q2 set to 1.00001'
1323  q2 = 1.00001
1324  ENDIF
1325  IF(q2.GE.1.0e4) THEN
1326 c$$$ WRITE(*,*) 'Warning: Q2 exceeds maximum value, Q2 set to 9999.9'
1327  q2 = 9999.9
1328  ENDIF
1329 
1330  xpq(-5) = x*sngl( ctq5pdf(-5, dble(x), dble(sqrt(q2))))
1331  xpq(-4) = x*sngl( ctq5pdf(-4, dble(x), dble(sqrt(q2))))
1332  xpq(-3) = x*sngl( ctq5pdf(-3, dble(x), dble(sqrt(q2))))
1333  xpq(-2) = x*sngl( ctq5pdf(-1, dble(x), dble(sqrt(q2))))
1334  xpq(-1) = x*sngl( ctq5pdf(-2, dble(x), dble(sqrt(q2))))
1335  xpq( 0) = x*sngl( ctq5pdf( 0, dble(x), dble(sqrt(q2))))
1336  xpq( 1) = x*sngl( ctq5pdf( 2, dble(x), dble(sqrt(q2))))
1337  xpq( 2) = x*sngl( ctq5pdf( 1, dble(x), dble(sqrt(q2))))
1338  xpq( 3) = x*sngl( ctq5pdf( 3, dble(x), dble(sqrt(q2))))
1339  xpq( 4) = x*sngl( ctq5pdf( 4, dble(x), dble(sqrt(q2))))
1340  xpq( 5) = x*sngl( ctq5pdf( 5, dble(x), dble(sqrt(q2))))
1341 
1342 C*****************************************************************************
1343 C LST(15) = 161 mrs99 cor01 central gluon, a_s
1344 C LST(15) = 162 mrs99 cor02 higher gluon
1345 C LST(15) = 163 mrs99 cor03 lower gluon
1346 C LST(15) = 164 mrs99 cor04 lower a_s
1347 C LST(15) = 165 mrs99 cor05 higher a_s
1348 C LST(15) = 166 mrs99 cor06 quarks up
1349 C LST(15) = 167 mrs99 cor07 quarks down
1350 C LST(15) = 168 mrs99 cor08 strange up
1351 C LST(15) = 169 mrs99 cor09 strange down
1352 C LST(15) = 170 mrs99 cor10 charm up
1353 C LST(15) = 171 mrs99 cor11 charm down
1354 C LST(15) = 172 mrs99 cor12 larger d/u
1355 C****************************************************************************
1356 
1357  ELSEIF((lst(15).GE.161).AND.(lst(15).LE.172)) THEN
1358 
1359  IF(x.LE.0.00001) THEN
1360 c$$$ WRITE(*,*) 'Warning: x exceeds minum value, x set to 0.0000101'
1361  x=0.0000101
1362  ENDIF
1363  IF(x.GE.1.0) THEN
1364 c$$$ WRITE(*,*) 'Warning: x exceeds maxinum value, x set to 0.99999'
1365  x=0.999999
1366  ENDIF
1367  IF(q2.LE.1.25) THEN
1368 c$$$ WRITE(*,*)'Warning: Q2 exceeds minimum value, Q2 set to 1.25001'
1369  q2 = 1.25001
1370  ENDIF
1371  IF(q2.GE.1.0e4) THEN
1372 c$$$ WRITE(*,*) 'Warning: Q2 exceeds maximum value, ',
1373 c$$$ & 'Q2 set to 9999999.9'
1374  q2 = 9999999.9
1375  ENDIF
1376 
1377  mode = lst(15)-160
1378  CALL mrs99(dble(x),dble(sqrt(q2)),mode,
1379  & mrsupv,mrsdnv,mrsusea,mrsdsea,mrsstr,mrschm,mrsbot,
1380  & mrsglu)
1381 
1382  xpq(-5) = sngl(mrsbot)
1383  xpq(-4) = sngl(mrschm)
1384  xpq(-3) = sngl(mrsstr)
1385  xpq(-2) = sngl(mrsusea)
1386  xpq(-1) = sngl(mrsdsea)
1387  xpq( 0) = sngl(mrsglu)
1388  xpq( 1) = sngl(mrsdsea)+sngl(mrsdnv)
1389  xpq( 2) = sngl(mrsusea)+sngl(mrsupv)
1390  xpq( 3) = sngl(mrsstr)
1391  xpq( 4) = sngl(mrschm)
1392  xpq( 5) = sngl(mrsbot)
1393 
1394 C Added By Joe Seele for the inclusion of the CTEQ 6 PDF's
1395 
1396  ELSEIF((lst(15).GE.173).AND.(lst(15).LE.175)) THEN
1397 
1398  IF(x.LE.0.00001) THEN
1399 c$$$ WRITE(*,*) 'Warning: x exceeds minum value, x set to 0.0000101'
1400  x=0.0000101
1401  ENDIF
1402  IF(x.GE.1.0) THEN
1403 c$$$ WRITE(*,*) 'Warning: x exceeds maxinum value, x set to 0.99999'
1404  x=0.999999
1405  ENDIF
1406  IF(q2.LE.1.0) THEN
1407 c$$$ WRITE(*,*)'Warning: Q2 exceeds minimum value, Q2 set to 1.00001'
1408  q2 = 1.00001
1409  ENDIF
1410  IF(q2.GE.1.0e4) THEN
1411 c$$$ WRITE(*,*) 'Warning: Q2 exceeds maximum value, Q2 set to 9999.9'
1412  q2 = 9999.9
1413  ENDIF
1414 
1415  xpq(-5) = x*sngl( ctq6pdf(-5, dble(x), dble(sqrt(q2))))
1416  xpq(-4) = x*sngl( ctq6pdf(-4, dble(x), dble(sqrt(q2))))
1417  xpq(-3) = x*sngl( ctq6pdf(-3, dble(x), dble(sqrt(q2))))
1418  xpq(-2) = x*sngl( ctq6pdf(-1, dble(x), dble(sqrt(q2))))
1419  xpq(-1) = x*sngl( ctq6pdf(-2, dble(x), dble(sqrt(q2))))
1420  xpq( 0) = x*sngl( ctq6pdf( 0, dble(x), dble(sqrt(q2))))
1421  xpq( 1) = x*sngl( ctq6pdf( 2, dble(x), dble(sqrt(q2))))
1422  xpq( 2) = x*sngl( ctq6pdf( 1, dble(x), dble(sqrt(q2))))
1423  xpq( 3) = x*sngl( ctq6pdf( 3, dble(x), dble(sqrt(q2))))
1424  xpq( 4) = x*sngl( ctq6pdf( 4, dble(x), dble(sqrt(q2))))
1425  xpq( 5) = x*sngl( ctq6pdf( 5, dble(x), dble(sqrt(q2))))
1426 
1427  IF(mcset_sfact.eq.2) THEN
1428  xpq( 3) = (x*sngl( ctq6pdf( 3, dble(x), dble(sqrt(q2)))))*
1429  1 ((x**(-0.3))*(1-x)*exp(-x/0.05))
1430  ENDIF
1431 
1432  IF(mcset_sbfact.eq.2) THEN
1433  xpq(-3) = (x*sngl( ctq6pdf(-3, dble(x), dble(sqrt(q2)))))*
1434  1 ((x**(-0.3))*(1-x)*exp(-x/0.05))
1435  ENDIF
1436 
1437 C*****************************************************************************
1438 
1439 
1440 
1441 
1442 C*************************************************************************
1443  ELSE
1444  WRITE(*,*) 'NOT KNOWN PARTON DENSITY FUNCTION: STOP',lst(15)
1445  stop
1446  ENDIF
1447 *72*********************************************************************
1448  DO i=-6,6
1449  IF(abs(xdpq(i)).gt.xpq(i)) THEN
1450  IF(icerr.lt.iplst(1).and.q2.gt.1.) THEN
1451  WRITE(*,*)'Non proper parton set : ',xdpq(i),xpq(i),i,x,q2
1452  icerr=icerr+1
1453  ENDIF
1454  IF(xdpq(i).gt.0) THEN
1455  xdpq(i)=xpq(i)
1456  ELSE
1457  xdpq(i)=-xpq(i)
1458  ENDIF
1459  ENDIF
1460  ENDDO
1461 *72*********************************************************************
1462 C This Amendment is to simulate HeliumIII
1463  IF(lst(39).NE.0) THEN
1464  f2p = 4.*(xpq(2)+xpq(-2)+xpq(4)+xpq(-4))
1465  & +xpq(1)+xpq(-1)+xpq(3)+xpq(-3)
1466  IF(lst(36).eq.0) THEN
1467  xdv=xpq(1)-xpq(-1)
1468  xuv=xpq(2)-xpq(-2)
1469  f2n = 4.*(xdv+2*xpq(-2)+xpq(4)+xpq(-4))
1470  & + xuv+2*xpq(-1)+xpq(3)+xpq(-3)
1471  ELSE
1472  f2n = 4.*(xpq(1)+xpq(-1)+xpq(4)+xpq(-4))+
1473  & xpq(2)+xpq(-2)+xpq(3)+xpq(-3)
1474  ENDIF
1475 
1476 C admitxurecoefficient for the unpolarized structure functions
1477  pari11 = 1./3.
1478 
1479  IF(lst(39).EQ.3) THEN
1480 C..........................
1481 C This is a real helium 3.
1482 C..........................
1483 C Percentage of neutron and proton to the helium3 asymmetry
1484  percs = 0.865
1485  percsp = 0.
1486  percd = -0.027
1487 C admixture coefficients for the spin states
1488  admxs = 1.0
1489  admxsp = 0.
1490  admxd = 0.
1491 C factors for the spin state (number of polarized nucleons/total number
1492 C of nucleons
1493  facts = f2n/(f2n+2*f2p)
1494  factsp = 0.
1495  factd = 2*f2p/(f2n+2*f2p)
1496  ELSEIF(lst(39).EQ.2) THEN
1497 C..........................
1498 C This is a real helium 3.
1499 C..........................
1500 C Percentage of the spin states
1501  percs = 0.903
1502  percsp = 0.014
1503  percd = -0.083
1504 C admixture coefficients for the spin states
1505  admxs = 1.0
1506  admxsp = 1./3.
1507  admxd = 1./3.
1508 C factors for the spin state (number of polarized nucleons/total number
1509 C of nucleons
1510  facts = f2n/(f2n+2*f2p)
1511  factsp = 1./3.
1512  factd = 1.
1513 C admitxurecoefficient for the unpolarized structure functions
1514  pari11 = 1./3.
1515  ELSEIF(lst(39).EQ.1) THEN
1516 C...........................
1517 C This is a ideal helium 3.
1518 C...........................
1519 C Percentage of the spin states
1520  percs = 1.0
1521  percsp = 0.0
1522  percd = 0.0
1523 C admixturecoefficients for the spin states
1524  admxs = 1.
1525  admxsp = 0.
1526  admxd = 0.
1527 C factors for the spin state (number of polarized nucleons/total number
1528 C of nucleons
1529  facts = f2n/(f2n+2*f2p)
1530  factsp = 0.
1531  factd = 0.
1532 C admitxurecoefficient for the unpolarized structure functions
1533  pari11 = 1./3.
1534  ENDIF
1535 
1536 C... For nuclear target, mix u- and d-valence distributions.
1537  IF(lst(36).EQ.0) THEN
1538 C... For nuclear target, mix u- and d-valence distributions.
1539  xpq(1)=(1.-pari11)*xdv+pari11*xuv+xpq(-1)
1540  xpq(2)=(1.-pari11)*xuv+pari11*xdv+xpq(-2)
1541 C Admixture in the polarized case
1542  xddv=xdpq(1)-xdpq(-1)
1543  xduv=xdpq(2)-xdpq(-2)
1544  xdpq(1)= ((1.-admxs )*xddv+admxs *xduv )
1545  & * facts * percs
1546  & +((1.-admxsp)*xddv+admxsp*xduv )
1547  & * factsp* percsp
1548  & +((1.-admxd )*xddv+admxd *xduv )
1549  & * factd * percd
1550  xdpq(2)= ((1.-admxs )*xduv+admxs *xddv )
1551  & * facts * percs
1552  & +((1.-admxsp)*xduv+admxsp*xddv )
1553  & * factsp* percsp
1554  & +((1.-admxd )*xduv+admxd *xddv )
1555  & * factd * percd
1556  DO i=-6,6
1557 C... For nuclear target, mix u- and d-valence distributions.
1558  IF(i.eq.1) THEN
1559 * XDPQ(-1) has to be diluted here by (FACTS*PERCS+FACTSP...
1560  xdpq(1)= xdpq(1) + xdpq(-1)
1561  ELSEIF(i.eq.2) THEN
1562 * XDPQ(-2) has to be diluted here by (FACTS*PERCS+FACTSP...
1563  xdpq(2)= xdpq(2) + xdpq(-2)
1564  ELSE
1565 * now for the sea quarks and gluon distributions
1566  xdpq(i)=xdpq(i)*
1567  & (facts*percs+factsp*percsp+factd*percd)
1568  ENDIF
1569  ENDDO
1570  ELSEIF(lst(36).EQ.1) THEN
1571  xd=xpq(1)
1572  xu=xpq(2)
1573  xpq(1) =(1.-pari11)*xd+pari11*xu
1574  xpq(2) =(1.-pari11)*xu+pari11*xd
1575  xd=xpq(-1)
1576  xu=xpq(-2)
1577  xpq(-1)=(1.-pari11)*xd+pari11*xu
1578  xpq(-2)=(1.-pari11)*xu+pari11*xd
1579 C Admixture in the polarized case
1580 *.............................................
1581 * FOR u and d QUARKS
1582 *.............................................
1583  xdd=xdpq(1)
1584  xdu=xdpq(2)
1585  xdds=xdpq(-1)
1586  xdus=xdpq(-2)
1587  DO i=-6,6
1588 C... For nuclear target, mix u- and d-valence distributions.
1589  IF(i.eq.1) THEN
1590  xdpq(1)= ((1.-admxs )*xdd+admxs *xdu )
1591  & * facts * percs
1592  & +((1.-admxsp)*xdd+admxsp*xdu )
1593  & * factsp* percsp
1594  & +((1.-admxd )*xdd+admxd *xdu )
1595  & * factd * percd
1596  ELSEIF(i.eq.2) THEN
1597  xdpq(2)= ((1.-admxs )*xdu+admxs *xdd )
1598  & * facts * percs
1599  & +((1.-admxsp)*xdu+admxsp*xdd )
1600  & * factsp* percsp
1601  & +((1.-admxd )*xdu+admxd *xdd )
1602  & * factd * percd
1603 * now for the sea quarks and gluon distributions
1604  ELSEIF(i.eq.-1) THEN
1605  xdpq(-1)=((1.-admxs)*xdds+admxs*xdus)
1606  & * facts * percs
1607  & +((1.-admxsp)*xdds+admxsp*xdus )
1608  & * factsp* percsp
1609  & +((1.-admxd )*xdds+admxd *xdus )
1610  & * factd * percd
1611  ELSEIF(i.eq.-2) THEN
1612  xdpq(-2)=((1.-admxs)*xdus+admxs*xdds )
1613  & * facts * percs
1614  & +((1.-admxsp)*xdus+admxsp*xdds )
1615  & * factsp* percsp
1616  & +((1.-admxd )*xdus+admxd *xdds )
1617  & * factd * percd
1618  ELSE
1619 * now for the sea quarks and gluon distributions
1620  xdpq(i)=xdpq(i)*
1621  & (facts*percs+factsp*percsp+factd*percd)
1622  ENDIF
1623  ENDDO
1624 C... For nuclear target, mix u- and d-valence distributions.
1625  ENDIF
1626  ENDIF
1627 
1628 * polarisation enhancement factor
1629  DO i=-6,6
1630  xdpq(i) = float(lst(38))* xdpq(i)
1631  ENDDO
1632 
1633  RETURN
1634  END