]>
Commit | Line | Data |
---|---|---|
4e9e3152 | 1 | subroutine CTEQ65evolve(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 | common/NAME/name,nmem,ndef,mmem | |
7 | real*8 pdf(-6:6) | |
8 | integer nset | |
9 | Character Line*80 | |
10 | PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 40) | |
11 | PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX) | |
12 | Common | |
13 | > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX) | |
14 | > / CtqPar2 / Nx, Nt, NfMx | |
15 | > / XQrange / Qini, Qmax, Xmin | |
16 | > / QCDtable / Alambda, Nfl, Iorder | |
17 | > / Masstbl / Amass(6) | |
18 | ||
19 | common/masses_LHA/cMass,bMass,tMass | |
20 | ||
21 | data pi / 3.141592653589793d0 / | |
22 | save | |
23 | c | |
24 | call getnset(iset) | |
25 | call getnmem(iset,imem) | |
26 | ||
27 | U = X * CtLhCtq65Pdf(imem,1,X,Q) | |
28 | D = X * CtLhCtq65Pdf(imem,2,X,Q) | |
29 | USEA = X * CtLhCtq65Pdf(imem,-1,X,Q) | |
30 | DSEA = X * CtLhCtq65Pdf(imem,-2,X,Q) | |
31 | STR = X * CtLhCtq65Pdf(imem,3,X,Q) | |
32 | CHM = X * CtLhCtq65Pdf(imem,-4,X,Q) | |
33 | BOT = X * CtLhCtq65Pdf(imem,-5,X,Q) | |
34 | GLU = X * CtLhCtq65Pdf(imem,0,X,Q) | |
35 | c | |
36 | pdf(0) = glu | |
37 | pdf(1) = d | |
38 | pdf(-1) = dsea | |
39 | pdf(2) = u | |
40 | pdf(-2) = usea | |
41 | pdf(3) = str | |
42 | pdf(-3) = str | |
43 | pdf(4) = chm | |
44 | pdf(-4) = chm | |
45 | pdf(5) = bot | |
46 | pdf(-5) = bot | |
47 | pdf(6) = 0.0d0 | |
48 | pdf(-6) = 0.0d0 | |
49 | return | |
50 | * | |
51 | entry CTEQ65cevolve(x,Q,pdf) !like cteq65evolve, but allows c.ne.cbar and s.ne.sbar | |
52 | c | |
53 | call getnset(iset) | |
54 | call getnmem(iset,imem) | |
55 | ||
56 | U = X * CtLhCtq65Pdf(imem,1,X,Q) | |
57 | D = X * CtLhCtq65Pdf(imem,2,X,Q) | |
58 | USEA = X * CtLhCtq65Pdf(imem,-1,X,Q) | |
59 | DSEA = X * CtLhCtq65Pdf(imem,-2,X,Q) | |
60 | STR = X * CtLhCtq65Pdf(imem,3,X,Q) | |
61 | SBAR = X * CtLhCtq65Pdf(imem,-3,X,Q) | |
62 | CHM = X * CtLhCtq65Pdf(imem,4,X,Q) | |
63 | CBAR = X * CtLhCtq65Pdf(imem,-4,X,Q) | |
64 | BOT = X * CtLhCtq65Pdf(imem,5,X,Q) | |
65 | GLU = X * CtLhCtq65Pdf(imem,0,X,Q) | |
66 | c | |
67 | pdf(0) = glu | |
68 | pdf(1) = d | |
69 | pdf(-1) = dsea | |
70 | pdf(2) = u | |
71 | pdf(-2) = usea | |
72 | pdf(3) = str | |
73 | pdf(-3) = sbar | |
74 | pdf(4) = chm | |
75 | pdf(-4) = cbar | |
76 | pdf(5) = bot | |
77 | pdf(-5) = bot | |
78 | pdf(6) = 0.0d0 | |
79 | pdf(-6) = 0.0d0 | |
80 | ||
81 | return | |
82 | * | |
83 | entry CTEQ65read(nset) | |
84 | call CtLhbldat1 | |
85 | call CtLhbldat2 | |
86 | ||
87 | call LHct65set | |
88 | ||
89 | read(1,*)nmem(nset),ndef(nset) !*** nmem+1=number of members; ndef is not used for anything *** | |
90 | if(nmem(nset) .gt. nhess) then | |
91 | print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess | |
92 | stop | |
93 | endif | |
94 | ||
95 | MxVal = 3 | |
96 | Read (1, '(A)') Line | |
97 | Read (1, '(A)') Line | |
98 | Read (1, *) Dr, Fl, Al, (Amass(I),I=1,6) | |
99 | Iorder = Nint(Dr) | |
100 | Nfl = Nint(Fl) | |
101 | Alambda = Al | |
102 | ||
103 | cMass = Amass(4) | |
104 | bMass = Amass(5) | |
105 | tMass = Amass(6) | |
106 | ||
107 | Read (1, '(A)') Line | |
108 | C This is the .pds (WKT) format | |
109 | Read (1, *) N0, N0, N0, NfMx, N0, N0 | |
110 | Read (1, '(A)') Line | |
111 | Read (1, *) NX, NT, N0, N0, N0 | |
112 | Read (1, '(A)') (Line,I=1,4) | |
113 | Read (1, *) QINI, QMAX, (aa,TV(I), I =0, NT) | |
114 | ||
115 | Read (1, '(A)') Line | |
116 | Read (1, *) XMIN, aa, (XV(I), I =1, NX) | |
117 | XV(0)=0D0 | |
118 | C | |
119 | C Since quark = anti-quark for nfl>2 at this stage, | |
120 | C we Read out only the non-redundent data points | |
121 | C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) | |
122 | ||
123 | ||
124 | Nblk = (NX+1) * (NT+1) | |
125 | Npts = Nblk * (NfMx+1+MxVal) | |
126 | ||
127 | do ihess = 0,nmem(nset) !*** new version: allows nmem < nhess *** | |
128 | ||
129 | Read (1, '(A)') Line | |
130 | Read (1, '(A)') Line | |
131 | Read (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts) | |
132 | ||
133 | enddo | |
134 | return | |
135 | * | |
136 | entry CTEQ65cread(nset) !like CTEQ65read, but c.ne.cbar allowed | |
137 | call CtLhbldat1 | |
138 | call CtLhbldat2 | |
139 | call LHct65set | |
140 | ||
141 | read(1,*)nmem(nset),ndef(nset) !*** nmem+1=number of members; ndef is not used for anything *** | |
142 | if(nmem(nset) .gt. nhess) then | |
143 | print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess | |
144 | stop | |
145 | endif | |
146 | ||
147 | MxVal = 4 !one more species free | |
148 | Read (1, '(A)') Line | |
149 | Read (1, '(A)') Line | |
150 | Read (1, *) Dr, Fl, Al, (Amass(I),I=1,6) | |
151 | Iorder = Nint(Dr) | |
152 | Nfl = Nint(Fl) | |
153 | Alambda = Al | |
154 | ||
155 | cMass = Amass(4) | |
156 | bMass = Amass(5) | |
157 | tMass = Amass(6) | |
158 | ||
159 | Read (1, '(A)') Line | |
160 | C This is the .pds (WKT) format | |
161 | Read (1, *) N0, N0, N0, NfMx, N0, N0 | |
162 | Read (1, '(A)') Line | |
163 | Read (1, *) NX, NT, N0, N0, N0 | |
164 | Read (1, '(A)') (Line,I=1,4) | |
165 | Read (1, *) QINI, QMAX, (aa,TV(I), I =0, NT) | |
166 | ||
167 | Read (1, '(A)') Line | |
168 | Read (1, *) XMIN, aa, (XV(I), I =1, NX) | |
169 | XV(0)=0D0 | |
170 | C | |
171 | C we Read only the non-redundent data points | |
172 | ||
173 | Nblk = (NX+1) * (NT+1) | |
174 | Npts = Nblk * (NfMx+1+MxVal) | |
175 | ||
176 | do ihess = 0,nmem(nset) !*** new version: allows nmem < nhess *** | |
177 | ||
178 | Read (1, '(A)') Line | |
179 | Read (1, '(A)') Line | |
180 | Read (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts) | |
181 | ||
182 | enddo | |
183 | return | |
184 | * | |
185 | entry CTEQ65alfa(alfas,Qalfa) | |
186 | alfas = pi*CtLhALPI(Qalfa) | |
187 | return | |
188 | * | |
189 | entry CTEQ65init(Eorder,Q2fit) | |
190 | return | |
191 | * | |
192 | entry CTEQ65pdf(mem) | |
193 | c imem = mem | |
194 | call getnset(iset) | |
195 | call setnmem(iset,mem) | |
196 | return | |
197 | * | |
198 | end | |
199 | ||
200 | subroutine LHct65set | |
201 | Implicit Double Precision (A-H,O-Z) | |
202 | common /ctq65co/ xlast, qlast, nxsave | |
203 | nxsave = -1000 | |
204 | xlast = -2. | |
205 | qlast = -2. | |
206 | return | |
207 | end | |
208 | ||
209 | c=========================================================================== | |
210 | Function CtLhPartonX65 (iset,IPRTN, XX, QQ) | |
211 | c Given the parton distribution function in the array U in | |
212 | c COMMON / PEVLDT / , this routine interpolates to find | |
213 | c the parton distribution at an arbitray point in x and q. | |
214 | c | |
215 | Implicit Double Precision (A-H,O-Z) | |
216 | ||
217 | PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 40) | |
218 | PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX) | |
219 | ||
220 | Common | |
221 | > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX) | |
222 | > / CtqPar2 / Nx, Nt, NfMx | |
223 | > / XQrange / Qini, Qmax, Xmin | |
224 | common /ctq65co/ xlast,qlast, nxsave | |
225 | ||
226 | parameter(nqvec = 4) | |
227 | ||
228 | Dimension fvec(4), fij(4) | |
229 | Dimension xvpow(0:mxx) | |
230 | Data OneP / 1.00001 / | |
231 | Data xpow / 0.3d0 / !**** choice of interpolation variable | |
232 | Save xvpow | |
233 | ||
234 | save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3 | |
235 | save const1 , const2, const3, const4, const5, const6 | |
236 | save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet | |
237 | ||
238 | c store the powers used for interpolation on first call... | |
239 | if(nx .ne. nxsave) then | |
240 | nxsave = nx | |
241 | xvpow(0) = 0.D0 | |
242 | do i = 1, nx | |
243 | xvpow(i) = xv(i)**xpow | |
244 | enddo | |
245 | endif | |
246 | ||
247 | X = XX | |
248 | Q = QQ | |
249 | ||
250 | if((x.lt.xmin).or.(x.gt.1.d0)) print 98,x | |
251 | 98 format(' WARNING: X=',e12.5,' OUT OF RANGE') | |
252 | if((q.lt.qini).or.(q.gt.qmax)) print 99,q | |
253 | 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE') | |
254 | ||
255 | c skip the initialization in x if same as in the previous call. | |
256 | if(x .eq. xlast) goto 100 | |
257 | xlast = x | |
258 | ||
259 | c ------------- find lower end of interval containing x, i.e., | |
260 | c get jx such that xv(jx) .le. x .le. xv(jx+1)... | |
261 | JLx = -1 | |
262 | JU = Nx+1 | |
263 | 11 If (JU-JLx .GT. 1) Then | |
264 | JM = (JU+JLx) / 2 | |
265 | If (X .Ge. XV(JM)) Then | |
266 | JLx = JM | |
267 | Else | |
268 | JU = JM | |
269 | Endif | |
270 | Goto 11 | |
271 | Endif | |
272 | C Ix 0 1 2 Jx JLx Nx-2 Nx | |
273 | C |---|---|---|...|---|-x-|---|...|---|---| | |
274 | C x 0 Xmin x 1 | |
275 | C | |
276 | If (JLx .LE. -1) Then | |
277 | Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX65 x=',x | |
278 | Stop | |
279 | ElseIf (JLx .Eq. 0) Then | |
280 | Jx = 0 | |
281 | Elseif (JLx .LE. Nx-2) Then | |
282 | ||
283 | C For interior points, keep x in the middle, as shown above | |
284 | Jx = JLx - 1 | |
285 | Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then | |
286 | ||
287 | C We tolerate a slight over-shoot of one (OneP=1.00001), | |
288 | C perhaps due to roundoff or whatever, but not more than that. | |
289 | C Keep at least 4 points >= Jx | |
290 | Jx = JLx - 2 | |
291 | Else | |
292 | Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX65 x=',x | |
293 | Stop | |
294 | Endif | |
295 | C ---------- Note: JLx uniquely identifies the x-bin; Jx does not. | |
296 | ||
297 | C This is the variable to be interpolated in | |
298 | ss = x**xpow | |
299 | ||
300 | If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then | |
301 | ||
302 | c initiation work for "interior bins": store the lattice points in s... | |
303 | svec1 = xvpow(jx) | |
304 | svec2 = xvpow(jx+1) | |
305 | svec3 = xvpow(jx+2) | |
306 | svec4 = xvpow(jx+3) | |
307 | ||
308 | s12 = svec1 - svec2 | |
309 | s13 = svec1 - svec3 | |
310 | s23 = svec2 - svec3 | |
311 | s24 = svec2 - svec4 | |
312 | s34 = svec3 - svec4 | |
313 | ||
314 | sy2 = ss - svec2 | |
315 | sy3 = ss - svec3 | |
316 | ||
317 | c constants needed for interpolating in s at fixed t lattice points... | |
318 | const1 = s13/s23 | |
319 | const2 = s12/s23 | |
320 | const3 = s34/s23 | |
321 | const4 = s24/s23 | |
322 | s1213 = s12 + s13 | |
323 | s2434 = s24 + s34 | |
324 | sdet = s12*s34 - s1213*s2434 | |
325 | tmp = sy2*sy3/sdet | |
326 | const5 = (s34*sy2-s2434*sy3)*tmp/s12 | |
327 | const6 = (s1213*sy2-s12*sy3)*tmp/s34 | |
328 | ||
329 | EndIf | |
330 | ||
331 | 100 continue | |
332 | ||
333 | c skip the initialization in q if same as in the previous call. | |
334 | if(q .eq. qlast) goto 110 | |
335 | qlast = q | |
336 | ||
337 | tt = log(log(Q/Al)) | |
338 | ||
339 | c --------------Now find lower end of interval containing Q, i.e., | |
340 | c get jq such that qv(jq) .le. q .le. qv(jq+1)... | |
341 | JLq = -1 | |
342 | JU = NT+1 | |
343 | 12 If (JU-JLq .GT. 1) Then | |
344 | JM = (JU+JLq) / 2 | |
345 | If (tt .GE. TV(JM)) Then | |
346 | JLq = JM | |
347 | Else | |
348 | JU = JM | |
349 | Endif | |
350 | Goto 12 | |
351 | Endif | |
352 | ||
353 | If (JLq .LE. 0) Then | |
354 | Jq = 0 | |
355 | Elseif (JLq .LE. Nt-2) Then | |
356 | C keep q in the middle, as shown above | |
357 | Jq = JLq - 1 | |
358 | Else | |
359 | C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq. | |
360 | Jq = Nt - 3 | |
361 | ||
362 | Endif | |
363 | C This is the interpolation variable in Q | |
364 | ||
365 | If (JLq.GE.1 .and. JLq.LE.Nt-2) Then | |
366 | c store the lattice points in t... | |
367 | tvec1 = Tv(jq) | |
368 | tvec2 = Tv(jq+1) | |
369 | tvec3 = Tv(jq+2) | |
370 | tvec4 = Tv(jq+3) | |
371 | ||
372 | t12 = tvec1 - tvec2 | |
373 | t13 = tvec1 - tvec3 | |
374 | t23 = tvec2 - tvec3 | |
375 | t24 = tvec2 - tvec4 | |
376 | t34 = tvec3 - tvec4 | |
377 | ||
378 | ty2 = tt - tvec2 | |
379 | ty3 = tt - tvec3 | |
380 | ||
381 | tmp1 = t12 + t13 | |
382 | tmp2 = t24 + t34 | |
383 | ||
384 | tdet = t12*t34 - tmp1*tmp2 | |
385 | ||
386 | EndIf | |
387 | ||
388 | 110 continue | |
389 | ||
390 | c get the pdf function values at the lattice points... | |
391 | c In this code, we store 10 flavors: u,ubar,d,dbar,s,sbar,c,cbar,b=bbar,g | |
392 | c hence Iprtn=5 (b) is obtained from -5 (bbar) | |
393 | ||
394 | If (Iprtn .GE. 5) Then | |
395 | Ip = - Iprtn | |
396 | Else | |
397 | Ip = Iprtn | |
398 | EndIf | |
399 | jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1 | |
400 | ||
401 | Do it = 1, nqvec | |
402 | ||
403 | J1 = jtmp + it*(NX+1) | |
404 | ||
405 | If (Jx .Eq. 0) Then | |
406 | C For the first 4 x points, interpolate x^2*f(x,Q) | |
407 | C This applies to the two lowest bins JLx = 0, 1 | |
408 | C We cannot put the JLx.eq.1 bin into the "interior" section | |
409 | C (as we do for q), since Upd(J1) is undefined. | |
410 | fij(1) = 0 | |
411 | fij(2) = Upd(iset,J1+1) * XV(1)**2 | |
412 | fij(3) = Upd(iset,J1+2) * XV(2)**2 | |
413 | fij(4) = Upd(iset,J1+3) * XV(3)**2 | |
414 | C | |
415 | C Use CtLhPolint which allows x to be anywhere w.r.t. the grid | |
416 | ||
417 | Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx) | |
418 | ||
419 | If (x .GT. 0D0) fvec(it) = Fx / x**2 | |
420 | C Pdf is undefined for x.eq.0 | |
421 | ElseIf (JLx .Eq. Nx-1) Then | |
422 | C This is the highest x bin: | |
423 | ||
424 | c** fix allow 4 consecutive elements with iset... mrw 19.9.2005 | |
425 | fij(1) = Upd(iset,j1) | |
426 | fij(2) = Upd(iset,j1+1) | |
427 | fij(3) = Upd(iset,j1+2) | |
428 | fij(4) = Upd(iset,j1+3) | |
429 | Call CtLhPolint4 (XVpow(Nx-3), Fij(1), 4, ss, Fx, Dfx) | |
430 | ||
431 | fvec(it) = Fx | |
432 | ||
433 | Else | |
434 | ||
435 | C for all interior points, use Jon's in-line function | |
436 | C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2) | |
437 | c (This is cubic spline interpolation, as used by cteq; it was | |
438 | c changed to polint in previous Durham releases (jcp).) | |
439 | sf2 = Upd(iset,J1+1) | |
440 | sf3 = Upd(iset,J1+2) | |
441 | ||
442 | Fvec(it) = (const5*(Upd(iset,J1) | |
443 | & - sf2*const1 + sf3*const2) | |
444 | & + const6*(Upd(iset,J1+3) | |
445 | & + sf2*const3 - sf3*const4) | |
446 | & + sf2*sy3 - sf3*sy2) / s23 | |
447 | ||
448 | Endif | |
449 | ||
450 | enddo | |
451 | C We now have the four values Fvec(1:4) | |
452 | c interpolate in t... | |
453 | ||
454 | If (JLq .LE. 0) Then | |
455 | C 1st Q-bin, as well as extrapolation to lower Q | |
456 | Call CtLhPolint4(TV(0), Fvec(1), 4, tt, ff, Dfq) | |
457 | ||
458 | ElseIf (JLq .GE. Nt-1) Then | |
459 | C Last Q-bin, as well as extrapolation to higher Q | |
460 | Call CtLhPolint4(TV(Nt-3), Fvec(1), 4, tt, ff, Dfq) | |
461 | Else | |
462 | C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2) | |
463 | C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for | |
464 | C the full range QV(0:Nt) (in contrast to XV) | |
465 | tf2 = fvec(2) | |
466 | tf3 = fvec(3) | |
467 | ||
468 | g1 = ( tf2*t13 - tf3*t12) / t23 | |
469 | g4 = (-tf2*t34 + tf3*t24) / t23 | |
470 | ||
471 | h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12 | |
472 | & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34) | |
473 | ||
474 | ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23 | |
475 | EndIf | |
476 | ||
477 | CtLhPartonX65 = ff | |
478 | ||
479 | Return | |
480 | End | |
481 | c=========================================================================== | |
482 | Function CtLhCtq65Pdf (iset,Iparton, X, Q) | |
483 | Implicit Double Precision (A-H,O-Z) | |
484 | Logical Warn | |
485 | Common | |
486 | > / CtqPar2 / Nx, Nt, NfMx | |
487 | > / QCDtable / Alambda, Nfl, Iorder | |
488 | ||
489 | Data Warn /.true./ | |
490 | save Warn | |
491 | ||
492 | If (X .lt. 0D0 .or. X .gt. 1D0) Then | |
493 | Print *, 'X out of range in CtLhCtq65Pdf: ', X | |
494 | Stop | |
495 | Endif | |
496 | If (Q .lt. Alambda) Then | |
497 | Print *, 'Q out of range in CtLhCtq65Pdf: ', Q | |
498 | Stop | |
499 | Endif | |
500 | ||
501 | c added to force pdf = 0.0 at x=1.0 exactly - mrw | |
502 | if(x .eq. 1.0d0) then | |
503 | CtLhCtq65Pdf = 0.0d0 | |
504 | return | |
505 | endif | |
506 | c | |
507 | If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then | |
508 | If (Warn) Then | |
509 | C put a warning for calling extra flavor. | |
510 | Warn = .false. | |
511 | Print *, 'Warning: Iparton out of range in CtLhCtq65Pdf: ' | |
512 | > , Iparton | |
513 | Endif | |
514 | CtLhCtq65Pdf = 0D0 | |
515 | Return | |
516 | Endif | |
517 | ||
518 | CtLhCtq65Pdf = CtLhPartonX65 (iset,Iparton, X, Q) | |
519 | if(CtLhCtq65Pdf.lt.0.D0) CtLhCtq65Pdf = 0.D0 | |
520 | ||
521 | Return | |
522 | ||
523 | C ******************** | |
524 | End |