]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.2.2/alphas.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
More shadow variable warnings fixed
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / alphas.f
CommitLineData
3c5d1739 1 real*8 function alphasPDF(Q)
2 implicit none
3 integer nset
4 real*8 Q,a
5 call getnset(nset)
6 call evolveAs(nset,Q,a)
7 alphasPDF=a
8 return
9 end
10*
11 real*8 function alphasPDFM(nset,Q)
12 implicit none
13 integer nset
14 real*8 Q,a
15 call evolveAs(nset,Q,a)
16 alphasPDFM=a
17 return
18 end
19*
20 subroutine GetQmass(nnf,mass)
21 implicit none
22 integer nset
23 real*8 Q,alphas,alfas0,scale0,mass
24 integer nnf,nf,order
25 nset = 1
26 call GetQmassM(nset,nnf,mass)
27 return
28 entry GetOrderAs(order)
29 nset = 1
30 call GetOrderAsM(nset,order)
31 return
32 end
33
34
35 subroutine evolveAs(nset,Q,alphas)
36 implicit none
37 include 'parmsetup.inc'
38 integer iset,imem
39 common/SET/iset,imem
40 integer nset
41 character*16 s1,s2,s3
42 real*8 Q,alphas,alfas0,scale0,Mc,Mb,Mt
43 real*8 Q0(nmxset)
44 real*8 AlfasQ(nmxset)
45 real*8 b0,b1,b2,L,As,mass
46 double precision CtLhALPI,CtLhAlphaNew
47 parameter (b0=1.2202,b1=0.4897,b2=0.1913)
48 integer n
49 integer order
50 integer EvlOrd(nmxset)
51 integer parm(nmxset)
52 integer Etype(nmxset)
53 integer Method(nmxset)
54 integer nnf,naf,nf(nmxset)
55 real*8 cmass,bmass,tmass
56 common/masses_LHA/cmass,bmass,tmass
57 save EvlOrd,parm,Etype,alfasQ,Q0,Method
58
59c equivalence (Mc,Massc)
60c equivalence (Mb,Massb)
61c equivalence (Mt,Masst)
62
63 call setnset(nset)
64c
65 if (method(nset).eq.0) then
66 if ((Etype(nset).eq.1).or.(Etype(nset).eq.2)) then
67 L=log(Q/Q0(nset))
68 As=alfasQ(nset)
69 if (Etype(nset).eq.2) call GetParmPDF(nset,parm,As)
70c print *,Etype,parm,As
71 if (EvlOrd(nset).eq.0) L=b0*L
72 if (EvlOrd(nset).eq.1) L=(b0+As*b1)*L
73 if (EvlOrd(nset).eq.2) L=(b0+As*b1+As**2*b2)*L
74 . -0.5*As**2*b0*b1/2d0*L**2
75 alphas=As/(1.0+As*L)
76 endif
77 endif
78 if (method(nset).eq.1) then
79 call alfasevolve(nset,alphas,Q)
80 elseif (method(nset).eq.2) then
81 call alphamrs(5,alphas,Q)
82 elseif (method(nset).eq.3) then
83 alphas = 3.1415926535898d0*CtLhALPI(Q)
84 elseif (method(nset).eq.4) then
85 alphas = CtLhAlphaNew(Q)
86 elseif (method(nset).eq.5) then
87 call alphamrs(3,alphas,Q)
88 elseif (method(nset).eq.6) then
89 call alphamrs(4,alphas,Q)
90 endif
91 return
92*
93 entry GetQmassM(nset,nnf,mass)
94 n=abs(nnf)
95 mass=0d0
96 if (n.eq.4) mass=cmass
97 if (n.eq.5) mass=bmass
98 if (n.eq.6) mass=tmass
99 return
100
101 entry GetNactive(naf,q)
102c compute nnfn = number of active flavors at scale q.
103 n = 3
104 if(q .ge. cmass) n = 4
105 if(q .ge. bmass) n = n + 1
106 if(q .ge. tmass) n = n + 1
107 naf = n
108 return
109
110 entry GetAlfas(nset,alfas0,scale0)
111 scale0=Q0(nset)
112 alfas0=alfasQ(nset)
113 if (Etype(nset).eq.2) call GetParmPDF(nset,parm(nset),alfas0)
114 return
115*
116 entry GetOrderAsM(nset,order)
117 order=EvlOrd(nset)
118 return
119*
120 entry InitAlphasPDF(nset)
121 Etype(nset)=-1
122 EvlOrd(nset)=-1
123 read(1,*) s1,s2,s3
124 if (index(s2,'lo').eq.1) EvlOrd(nset)=0
125 if (index(s2,'nlo').eq.1) EvlOrd(nset)=1
126 if (index(s2,'nnlo').eq.1) EvlOrd(nset)=2
127 if (EvlOrd(nset).lt.0) then
128 write(*,*) 'File description error:'
129 write(*,*) 'Unknown alpha_s evolution order ',s2
130 stop
131 endif
132 if (index(s1,'Fixed').eq.1) then
133 Etype(nset)=1
134 parm(nset)=-1
135 read(1,*) alfasQ(nset),Q0(nset),cmass,bmass,tmass
136 endif
137 if (index(s1,'Variable').eq.1) then
138 Etype(nset)=2
139 alfasQ(nset)=0d0
140 read(1,*) parm(nset),Q0(nset),cmass,bmass,tmass
141 endif
142 if (Etype(nset).lt.0) then
143 write(*,*) 'File description error:'
144 write(*,*) 'Unknown alpha_s evolution method ',s1
145 stop
146 endif
147 Method(nset)=-1
148 if (index(s3,'Internal').eq.1) Method(nset)=0
149 if (index(s3,'EvolCode').eq.1) Method(nset)=1
150 if (index(s3,'MRSTalfa').eq.1) Method(nset)=2
151 if (index(s3,'CTEQalfa').eq.1) Method(nset)=3
152 if (index(s3,'CTEQABalfa').eq.1) Method(nset)=4
153 if (index(s3,'MRST3alfa').eq.1) Method(nset)=5
154 if (index(s3,'MRST4alfa').eq.1) Method(nset)=6
155 if (Method(nset).lt.0) then
156 write(*,*) 'File description error:'
157 write(*,*) 'Unknown alpha_s method ',s3
158 stop
159 endif
160c print *,'alpha_s evolution method = ',method
161 return
162 end
163
164 subroutine alphamrs(nflav,alpha,Q)
165 IMPLICIT real*8 (a-h,o-z)
166 common/SET/iset,imem
167 common/masses_LHA/cmass,bmass,tmass
168 dimension parms(31)
169 DATA PI/3.14159/
170 DATA TOL/.0005/
171 data memold/-1/
172 integer nset
173 save alambda,memold
174c
175c print *,'calling alphamrs with ',nflav,Q
176 call getnset(nset)
177c
178c
179c alambda=0.323
180c alambda=0.3342
181c -- find alambda corresponding to alphas given in first parameter
182 call getnmem(nset,imem)
183 if(imem.ne.memold) then
184 call listPDF(nset,imem,parms)
185 alphas = parms(1)
186 print *,alphas
187 memold=imem
188 qz2=8315.
189 qz = dsqrt(qz2)
190 alambda = 0.3000
191 astep = 0.010
192 tol2 = 0.0000001
193 idir = +1
194 10 continue
195 alambda = alambda + idir*astep
196 call mrslambda(nflav,alpha,qz,alambda)
197
198 if(idir*(alphas-alpha).gt.0.0) then
199 go to 20
200 else
201 astep = 0.5*astep
202 idir = -1*idir
203 go to 20
204 endif
205 20 continue
206 if(abs(alpha-alphas).gt.tol2) go to 10
207 endif
208c - alambda found -- save it !!!!
209c - next call mrslambda to get alphas at q with the correct alambda
210 call mrslambda(nflav,alpha,q,alambda)
211 RETURN
212 END
213*
214 subroutine rgras(alpha,Q2)
215 IMPLICIT real*8 (a-h,o-z)
216 real*8 mc,mb,mt
217 include 'parmsetup.inc'
218 character*16 name(nmxset)
219 integer nmem(nmxset),ndef(nmxset),mem
220 common/NAME/name,nmem,ndef,mem
221 q=dsqrt(q2)
222 call getnset(nset)
223 nflav=5
224 if(name(nset).eq.'QCDNUM_MRST3') nflav=3
225 if(name(nset).eq.'QCDNUM_MRST4') nflav=4
226 call alphamrs(nflav,alpha,Q)
227 return
228 end
229*
230* new vewrsion of mrslambda 13/5/2004 - includes lo, nlo, and nnlo with flags
231*
232 subroutine mrslambda(nflav,alpha,Q,alambda)
233 IMPLICIT REAL*8(A-H,O-Z)
234 DATA PI/3.14159/
235 DATA TOL/.0005/
236c The value of Lambda required corresponds to nflav=4
237c iord=0 gives leading order, iord=1 gives NLO, iord=2 gives NNLO
238 qsdt=8.18 !! This is the value of 4m_c^2
239 qsct=74.0 !! This is the value of 4m_b^2
240c print *,'calling mrslambda with ',nflav,Q,alambda
241 al2=alambda*alambda
242 q2=q*q
243 t=dlog(q2/al2)
244c
245c next line to be replaced by a getiord in the final version
246c iord=1
247c
248 call getnset(nset)
249 call GetOrderAsM(nset,iord)
250 ITH=0
251 TT=T
252 qsdtt=qsdt/4.
253 qsctt=qsct/4.
254 AL=ALAMBDA
255 AL2=AL*AL
256 FLAV=4.
257 if(nflav.eq.3) flav=3.
258 QS=AL2*dEXP(T)
259
260 if(qs.lt.0.5d0) then !! running stops below 0.5
261 qs=0.5d0
262 t=dlog(qs/al2)
263 tt=t
264 endif
265
266c if(QS.lt.QSDTT.and.nflav.eq.4) then
267c flav=3.
268c go to 11
269c endif
270
271
272 IF(QS.gt.QSCTT.and.nflav.gt.4) GO TO 12
273 IF(QS.lt.QSDTT.and.nflav.gt.3) GO TO 312
274 11 CONTINUE
275 B0=11-2.*FLAV/3.
276 X1=4.*PI/B0
277 IF(IORD.eq.0) then
278 ALPHA=X1/T
279 ELSE
280 if(iord.gt.1) then
281 alpha=qwikalf(t,iord,flav)
282 go to 51
283 endif
284 B1=102.-38.*FLAV/3.
285 X2=B1/B0**2
286 AS2=X1/T*(1.-X2*dLOG(T)/T)
287 5 AS=AS2
288 F=-T+X1/AS-X2*dLOG(X1/AS+X2)
289 FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
290 AS2=AS-F/FP
291 DEL=ABS(F/FP/AS)
292 IF((DEL-TOL).GT.0.) go to 5
293 ALPHA=AS2
294 51 continue
295 ENDIF
296 IF(ITH.EQ.0) RETURN
297 GO TO (13,14,15) ITH
298c GO TO 5
299 12 ITH=1
300 T=dLOG(QSCTT/AL2)
301 GO TO 11
302 13 ALFQC4=ALPHA
303 FLAV=5.
304 ITH=2
305
306 GO TO 11
307 14 ALFQC5=ALPHA
308 ITH=3
309 T=TT
310 GO TO 11
311 15 ALFQS5=ALPHA
312 ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
313 ALPHA=1./ALFINV
314 RETURN
315
316 311 CONTINUE
317 B0=11-2.*FLAV/3.
318 X1=4.*PI/B0
319 IF(IORD.eq.0) then
320 ALPHA=X1/T
321
322 ELSE
323 if(iord.gt.1) then
324 alpha=qwikalf(t,iord,flav)
325 go to 351
326 endif
327 B1=102.-38.*FLAV/3.
328 X2=B1/B0**2
329 AS2=X1/T*(1.-X2*dLOG(T)/T)
330 35 AS=AS2
331 F=-T+X1/AS-X2*dLOG(X1/AS+X2)
332 FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
333 AS2=AS-F/FP
334 DEL=ABS(F/FP/AS)
335 IF((DEL-TOL).GT.0.) go to 35
336 ALPHA=AS2
337
338 351 continue
339 endif
340 IF(ITH.EQ.0) RETURN
341 GO TO (313,314,315) ITH
342 312 ITH=1
343 T=dLOG(QSDTT/AL2)
344 GO TO 311
345 313 ALFQC4=ALPHA
346 FLAV=3.
347 ITH=2
348 GO TO 311
349 314 ALFQC3=ALPHA
350 ITH=3
351 T=TT
352 GO TO 311
353 315 ALFQS3=ALPHA
354 ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
355 ALPHA=1./ALFINV
356 RETURN
357 END
358
359
360
361 double precision function qwikalf(t,iord,flav)
362 implicit real*8(a-h,o-z)
363 dimension z3(6),z4(6),z5(6),zz3(6),zz4(6),zz5(6)
364 data z3/ -.161667E+01,0.954244E+01,
365 .0.768623E+01,0.101523E+00,-.360127E-02,0.457867E-04/
366 data z4/ -.172239E+01,0.831185E+01,
367 .0.721463E+01,0.835531E-01,-.285436E-02,0.349129E-04/
368 data z5/ -.872190E+00,0.572816E+01,
369 .0.716119E+01,0.195884E-01,-.300199E-03,0.151741E-05/
370 data zz3/-.155611E+02,0.168406E+02,
371 .0.603014E+01,0.257682E+00,-.970217E-02,0.127628E-03/
372 data zz4/-.106762E+02,0.118497E+02,0.664964E+01,
373 .0.112996E+00,-.317551E-02,0.302434E-04/
374 data zz5/-.531860E+01,0.708503E+01,0.698352E+01,
375 .0.274170E-01,-.426894E-03,0.217591E-05/
376
377 data pi/3.14159/
378 nfm2=flav-2.
379 x=dsqrt(t)
380 x2=x*x
381 x3=x*x2
382 x4=x*x3
383 x5=x*x4
384 go to (1,2) iord
385 1 go to (3,4,5) nfm2
386 3 y=z3(1)+z3(2)*x+z3(3)*x2+z3(4)*x3+z3(5)*x4+z3(6)*x5
387 go to 10
388 4 y=z4(1)+z4(2)*x+z4(3)*x2+z4(4)*x3+z4(5)*x4+z4(6)*x5
389 go to 10
390 5 y=z5(1)+z5(2)*x+z5(3)*x2+z5(4)*x3+z5(5)*x4+z5(6)*x5
391 go to 10
392 2 go to (6,7,8) nfm2
393 6 y=zz3(1)+zz3(2)*x+zz3(3)*x2+zz3(4)*x3+zz3(5)*x4+zz3(6)*x5
394 go to 10
395 7 y=zz4(1)+zz4(2)*x+zz4(3)*x2+zz4(4)*x3+zz4(5)*x4+zz4(6)*x5
396 go to 10
397 8 y=zz5(1)+zz5(2)*x+zz5(3)*x2+zz5(4)*x3+zz5(5)*x4+zz5(6)*x5
398 go to 10
399 10 qwikalf=4.*pi/y
400 return
401 end
402c=====================================================================
403c next is alphas routine from PDFLIB
404c
405c DOUBLE PRECISION FUNCTION ALPHAS2(SCALE,qcdl5)
406 subroutine aspdflib(alphas2,SCALE,iord,qcdl5)
407 implicit real*8 (a-h,o-z)
408 double precision
409 + NF,KNF
410c CHARACTER*20 PARM(NCHDIM)
411c double precision
412c + VAL(NCHDIM)
413c DATA XMC/1.5D0/,XMB/4.75D0/,XMT/100.D0/
414 DATA XMC/1.43D0/,XMB/4.30D0/,XMT/100.D0/
415 DATA ZEROD/0.D0/,PONED/0.001D0/,ONED/1.D0/,TWOD/2.D0/
416
417c QCDL5=0.165d0
418
419c print *,qcdl5,iord
420 LO = 1
421 if(iord.ne.0) LO = 2
422c if(LO.eq.1) then
423c print *,'Calculating ALPHA_S at Leading Order',LO
424c else
425c print *,'Calculating ALPHA_S at Next to Leading Order',LO
426c endif
427
428 TMAS = 180.0d0
429
430
431 ALPHAS2 = ZEROD
432C.
433 PI=4.0D0*ATAN(ONED)
434 B6 = (33.D0-2.D0*6.D0)/PI/12.D0
435 BP6 = (153.D0 - 19.D0*6.D0) / PI / TWOD / (33.D0 - 2.D0*6.D0)
436 B5 = (33.D0-2.D0*5.D0)/PI/12.D0
437 BP5 = (153.D0 - 19.D0*5.D0) / PI / TWOD / (33.D0 - 2.D0*5.D0)
438 B4 = (33.D0-2.D0*4.D0)/PI/12.D0
439 BP4 = (153.D0 - 19.D0*4.D0) / PI / TWOD / (33.D0 - 2.D0*4.D0)
440 B3 = (33.D0-2.D0*3.D0)/PI/12.D0
441 BP3 = (153.D0 - 19.D0*3.D0) / PI / TWOD / (33.D0 - 2.D0*3.D0)
442 XLC = TWOD * LOG( XMC/QCDL5)
443 XLB = TWOD * LOG( XMB/QCDL5)
444 XLT = TWOD * LOG( XMT/QCDL5 * TMAS/XMT)
445 XLLC = LOG( XLC)
446 XLLB = LOG( XLB)
447 XLLT = LOG( XLT)
448 C65 = ONED/( ONED/(B5 * XLT) - XLLT*BP5/(B5 * XLT)**2 )
449 + - ONED/( ONED/(B6 * XLT) - XLLT*BP6/(B6 * XLT)**2 )
450 C45 = ONED/( ONED/(B5 * XLB) - XLLB*BP5/(B5 * XLB)**2 )
451 + - ONED/( ONED/(B4 * XLB) - XLLB*BP4/(B4 * XLB)**2 )
452 C35 = ONED/( ONED/(B4 * XLC) - XLLC*BP4/(B4 * XLC)**2 )
453 + - ONED/( ONED/(B3 * XLC) - XLLC*BP3/(B3 * XLC)**2 ) + C45
454C
455 Q = SCALE
456 XLQ = TWOD * LOG( Q/QCDL5 )
457 XLLQ = LOG( XLQ )
458c KNF = NFL
459c NF = KNF
460
461c IF ( NF .LT. ZEROD) THEN
462 IF ( Q .GT. XMT * TMAS/XMT) THEN
463 NF = 6.D0
464 ELSEIF ( Q .GT. XMB ) THEN
465 NF = 5.D0
466 ELSEIF ( Q .GT. XMC ) THEN
467 NF = 4.D0
468 ELSE
469 NF = 3.D0
470 ENDIF
471c ENDIF
472c print *,nf
473 IF(NF .GT. 6.D0) NF = 6.D0
474 IF ( NF .EQ. 6.D0 ) THEN
475 ALF = ONED/(ONED/(ONED/(B6*XLQ)- BP6/(B6*XLQ)**2*XLLQ) + C65)
476 IF (LO.EQ.1) ALF = ONED/B6/XLQ
477 ELSEIF ( NF .EQ. 5.D0 ) THEN
478 ALF = ONED/(B5 * XLQ) - BP5/(B5 * XLQ)**2 * XLLQ
479 IF (LO.EQ.1) ALF = ONED/B5/XLQ
480 ELSEIF ( NF .EQ. 4.D0 ) THEN
481 ALF = ONED/(ONED/(ONED/(B4*XLQ)- BP4/(B4*XLQ)**2*XLLQ) + C45)
482 IF (LO.EQ.1) ALF = ONED/B4/XLQ
483 ELSEIF ( NF .EQ. 3.D0 ) THEN
484 ALF = ONED/(ONED/(ONED/(B3*XLQ)- BP3/(B3*XLQ)**2*XLLQ) + C35)
485 IF (LO.EQ.1) ALF = ONED/B3/XLQ
486 ELSE
487 WRITE(*,*)'Error in Alphas2'
488 STOP
489 ENDIF
490 ALPHAS2 = ALF
491 RETURN
492 END
493c ========================================================================
494 SUBROUTINE CtLhAlphaNewSET(MC,MB,MT,Q0,ALPHA0,IORDER,IMODE)
495c call to set quark masses for alpha_s, and choose lambda or its
496c equivalent to make alpha_s take the value alpha0 at scale Q0.
497
498 IMPLICIT NONE
499 DOUBLE PRECISION MC, MB, MT, Q0, ALPHA0
500 INTEGER IORDER, IMODE
501
502 DOUBLE PRECISION UDSCBT, QQ0, AALPHA0
503 INTEGER IIORDER, IIMODE
504 COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
505 COMMON /ALSCALE/ QQ0, AALPHA0
506
507 double precision Adummy
508 integer Nfl, Idummy
509 common / QCDtable / Adummy, Nfl, Idummy
510
511 if((imode .lt. 1) .or. (imode .gt. 3)) then
512 print *,'CtLhAlphaNewSET: fatal imode=',imode
513 stop
514 endif
515
516 IIMODE = IMODE
517 QQ0 = Q0
518 AALPHA0 = ALPHA0
519 IIORDER = IORDER
520
521 UDSCBT(1) = .005D0
522 UDSCBT(2) = .010D0
523 UDSCBT(3) = .300D0
524 UDSCBT(4) = MC
525 UDSCBT(5) = MB
526 UDSCBT(6) = MT
527
528c set artificial quark masses, if necessary, in alpha_s to enforce
529c the requested maximum number of flavors...
530 if(Nfl .le. 5) UDSCBT(6) = 6.d99
531 if(Nfl .le. 4) UDSCBT(5) = 5.d99
532 if(Nfl .le. 3) UDSCBT(4) = 4.d99
533 if(Nfl .le. 2) UDSCBT(3) = 3.d99
534 if(Nfl .le. 1) UDSCBT(2) = 2.d99
535 if(Nfl .le. 0) UDSCBT(1) = 1.d99
536
537c print *,'CtLhAlphaNewSET: imode=',imode,' Nfl=',Nfl !*** temporary ***
538c print *,'CtLhAlphaNewSET: mc,mb,mt=', mc, mb, mt !*** temporary ***
539c print *,'CtLhAlphaNewSET: Q0, alpha0=',Q0, alpha0 !*** temporary ***
540
541
542 RETURN
543 END
544
545 FUNCTION CtLhAlphaNew(Q)
546 IMPLICIT NONE
547 DOUBLE PRECISION CtLhAlphaNew, Q
548
549 DOUBLE PRECISION Q2, CtLhQALPHAS, UDSCBT
550 integer nf,ier
551
552 INTEGER IIMODE, IIORDER
553 COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
554
555 if((iimode .ge. 1) .and. (iimode .le. 3)) then
556 Q2=Q*Q
557 nf=6
558 if (Q2.lt. UDSCBT(6)**2) nf=5
559 if (Q2.lt. UDSCBT(5)**2) nf=4
560 if (Q2.lt. UDSCBT(4)**2) nf=3
561 if (Q2.lt. UDSCBT(3)**2) nf=2
562 if (Q2.lt. UDSCBT(2)**2) nf=1
563 if (Q2.lt. UDSCBT(1)**2) nf=0
564
565c external maximum number of flavors -- typically, nfmax=5 so
566c top quark is not a parton even at large Q...
567
568 CtLhAlphaNew=CtLhQALPHAS(Q2,nf,ier)
569 if(ier .ne. 0) then
570 print *,'warning in CtLhAlphaNew, Q=',Q,' nf=',nf,
571 & ' ier=',ier,' CtLhAlphaNew=',CtLhAlphaNew
572 endif
573
574 else
575 print *,'CtLhAlphaNew: undefined mode=',iimode
576 stop
577 endif
578
579
580 return
581 end
582
583 FUNCTION CtLhQALPHAS(QQ2,NF,IERR)
584 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
585
586 COMMON /ALSCALE/ Q0, ALPHA0
587 COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
588
589 Q02 = Q0**2
590 ALP0 = ALPHA0
591 IOR = IIORDER
592
593 NFFF = NF
594 CtLhQALPHAS = CtLhA0TOA1(QQ2,Q02,ALP0,IOR,NFFF,IERR)
595
596 RETURN
597 END
598
599 FUNCTION CtLhA0TOA1(QSU,QS0,AS0,IORD,NFF,IERR)
600 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
601
602 COMMON /QMASSES/ UDSCBT(6), IIMODE, IIORDER
603
604 QS1 = QSU
605
606 QMU0 = SQRT(QS0)
607 QMU1 = SQRT(QS1)
608
609 DO I = 1, 6
610 IF(QMU0.GE.UDSCBT(I)) NF0 = I
611 IF(QMU1.GE.UDSCBT(I)) NF1 = I
612 ENDDO
613
614 IF(NF1.LT.NF0) THEN
615 IST = -1
616 JST = 0
617 ELSE
618 IST = 1
619 JST = 1
620 ENDIF
621
622 ALFA0 = AS0
623 Q00 = QS0
624
625 IERR = 0
626 DO 50 NF = NF0,NF1,IST
627
628 IF(NF.NE.NF1) THEN
629 Q21 = UDSCBT(NF+JST)**2
630 ELSE
631 Q21 = QS1
632 ENDIF
633 IMODE = IIMODE
634 ALFA1 = CtLhALPHAR(Q21,Q00,ALFA0,NF,IORD,IMODE,JERR)
635 IERR = IERR + JERR !IERR is sum of all errors
636 ALFA0 = ALFA1
637 Q00 = Q21
638
639 50 CONTINUE
640
641 CtLhA0TOA1 = ALFA0
642 NFF = NF1
643
644 RETURN
645 END
646
647 FUNCTION CtLhALPHAR(QSQ,QS0,AS0,NF,IORD,IMODE,IERR)
648C calculate ALPHAS FROM RGE GIVEN AS0 AT QS0.
649C
650C IORD=1: LEADING ORDER DEFINED BY
651C Q*d(alpha)/d(Q) = c1*alpha**2
652C
653C IORD=2,IMODE=1: QCD NUM CHOICE DEFINED BY
654C Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
655C
656C IORD=2,IMODE=2: AD HOC ALTERNATIVE DEFINED BY
657C Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
658C
659C IORD=2,IMODE=3: TRADITIONAL CTEQ CHOICE DEFINED BY
660C ALPHA = c3*(1 - c4*log(L)/L)/L, WHERE L=log((Q/lambda)**2)
661C
662C c1 = -beta0/(2*pi) where beta0 = 11. - (2./3.)*nf
663C c2 = -beta1/(8*pi**2) where beta1 = 102. - (38./3.)*nf
664C
665C c3 = -2/c1
666C c4 = -2*c2/c1**2
667C
668C
669 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
670
671
672 DATA PI / 3.14159265358979d0 /
673
674 BET0 = 11.d0 -(2*NF)/3.d0
675 BET1 = 102.d0 - (38*NF)/3.d0
676 B0 = BET0/(4.d0*PI)
677 B1 = BET1/(4.d0*PI*BET0)
678 IERR = 0
679
680 TERM0 = 1.d0/AS0+B0*LOG(QSQ/QS0)
681 IF(TERM0.LE.0.) THEN
682 CtLhALPHAR = 100.
683 IERR = 1
684 PRINT *,'CtLhALPHAR WARNING: RETURN 100.'
685 RETURN
686 ENDIF
687 ALFA0 = 1.d0/TERM0
688
689C ORDER=1 IS LEADING ORDER, WHICH IS SAME FOR ALL IMODE.
690 IF(IORD.EQ.1) THEN
691 CtLhALPHAR = ALFA0
692 RETURN
693 ELSEIF(IORD.NE.2) THEN
694 PRINT *,'FATAL ERROR: UNDEFINED ORDER IN CtLhALPHAR'
695 STOP
696 ENDIF
697
698
699C QCD NUM CHOICE: Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
700 IF(IMODE .EQ. 1) THEN
701
702c use Newton's method to solve the equation, instead of the
703c simple iterative method used in qcdnum (jcp 9/01)
704 DO ITER = 1, 20
705 ARG = (1.d0/ALFA0+B1)/(1.d0/AS0+B1)
706 IF(ARG.LE.0.) THEN
707 CtLhALPHAR = 10.
708 IERR = 1
709 PRINT *,'CtLhALPHAR WARNING: RETURN 10.'
710 RETURN
711 ENDIF
712
713 TERM = TERM0 + B1*LOG(ARG) - 1.d0/ALFA0
714 ALFA1 = ALFA0/(1.d0 + ALFA0*(1.d0 + B1*ALFA0)*TERM)
715
716 IF(ABS(ALFA1-ALFA0).LT.1.E-12) GOTO 20
717 ALFA0 = ALFA1
718 ENDDO
719
720 CtLhALPHAR = 10.
721 IERR = 1
722 RETURN
723
72420 CONTINUE
725 CtLhALPHAR = ALFA1
726 RETURN
727
728C AD HOC ALTERNATIVE: Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
729 ELSEIF(IMODE .EQ. 2) THEN
730
731c first get a good starting point, to be sure Newton's method doesn't go
732c to the wrong root.
733 BEST = 9.d99
734 DO ITRY = 0, 20
735 IF(ITRY.NE.0) ALFA0 = (1./B1)*(ITRY-0.5D0)/20.D0
736 F = -1.d0/ALFA0 + TERM0 - B1*LOG(ALFA0/AS0)
737 IF(ABS(F) .LT. BEST) THEN
738 BEST = ABS(F)
739 ALBST = ALFA0
740 ENDIF
741 ENDDO
742
743 ALFA0 = ALBST
744 DO ITER=1, 20
745 F = -1.d0/ALFA0 + TERM0 - B1*LOG(ALFA0/AS0)
746 ALFA1 = ALFA0/(1.d0 + ALFA0*F/(1.d0 - B1*ALFA0))
747 IF(ABS(ALFA1-ALFA0) .LT. 1.E-12) GOTO 30
748 ALFA0 = ALFA1
749 ENDDO
750
751 CtLhALPHAR = 10.
752 IERR = 1
753 RETURN
754
755 30 CONTINUE
756 CtLhALPHAR = ALFA1
757 RETURN
758
759C TRADITIONAL CTEQ CHOICE: ALPHA = c3*(1 - c4*log(L)/L)/L, WHERE L=log((Q/lambda)**2)
760 ELSEIF(IMODE .EQ. 3) THEN
761
762 Z = -LOG(B0*AS0)
763 TMP = BET1/BET0**2
764
765 DO ITER = 1, 20
766 F = EXP(Z) - (1.D0 - TMP*Z*EXP(-Z))/(B0*AS0)
767 FPRI = EXP(Z) + TMP*(1.D0-Z)*EXP(-Z)/(B0*AS0)
768 ZNEW = Z - F/FPRI
769 IF(ABS(Z-ZNEW) .LT. 1.E-10) GOTO 40
770 Z = ZNEW
771 ENDDO
772
773 CtLhALPHAR = 10.
774 IERR = 1
775 RETURN
776
777 40 CONTINUE
778 XLAMSQ = QS0 * EXP(-EXP(ZNEW))
779 XL = LOG(QSQ/XLAMSQ)
780
781c return a fixed value if no solution...
782 IF(XL .LE. 0.D0) THEN
783 CtLhALPHAR = 10.D0
784 IER = 1
785 RETURN
786 ENDIF
787
788 CtLhALPHAR = (1.d0 - TMP*LOG(XL)/XL)/(B0*XL)
789
790c place a cutoff if comes out very large...
791 if(CtLhALPHAR .gt. 10.d0) then
792 CtLhALPHAR = 10.d0
793 ier = 1
794 endif
795
796 RETURN
797
798 ELSE
799 PRINT *,'FATAL UNDEFINED IMODE=',IMODE
800 STOP
801 ENDIF
802
803 END
804
805 FUNCTION CtLhALPInew (AMU)
806C Returns effective g**2/(4pi**2) = alpha/pi.
807 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
808 data pi / 3.14159265358979d0 /
809
810 q = amu
811 alpha = CtLhAlphaNew(q)
812 CtLhalpinew = alpha/pi
813
814 RETURN
815 END
816
817 subroutine CTEQ6NewAlpha(nset,mem)
818 IMPLICIT DOUBLE PRECISION (A-H, O-Z)
819 common / QCDtable / Alambda, Nfl, Iorder
820 dimension parms(2)
821
822 q0 = 91.188d0
823c call GetLam4M(nset,mem,alpha0)
824c call GetLam5M(nset,mem,ximode)
825 call listPDF(nset,mem,parms)
826 alpha0 = parms(1)
827 ximode = parms(2)
828 imode = int(ximode)
829 iiorder = iorder
830
831c **********************************************************
832c the input data file should probably specify a very large
833c value for xmt, because we never allow top quark as a parton
834c in the PDF fitting, so there are never more than 5 active
835c flavors. Presume it is most consistent therefore to
836c also keep top quark loop out of running of alpha.
837c **********************************************************
838
839 call GetQmass(4,cmass)
840 call GetQmass(5,bmass)
841 call GetQmass(6,tmass)
842
843 write(6,*) 'CTEQ6NewAlpha: mc, mb, mt=',
844 & cmass, bmass, tmass
845
846 call CtLhAlphaNewSET(XMC,XMB,XMT,Q0,ALPHA0,IIORDER,IMODE)
847
848 RETURN
849 END
850
851c=================================================
852 subroutine getnset(nset)
853 integer iset,nset
854 save iset
855 nset = iset
856c print *,'getting nset:',nset
857 return
858c
859 entry setnset(nset)
860 iset = nset
861c print *,'setting nset:',nset
862 return
863 end
864c=================================================
865 subroutine getnmem(nset,nmem)
866 include 'parmsetup.inc'
867 integer nmem,nset,member(nmxset)
868 save member
869 nmem = member(nset)
870 return
871c
872 entry setnmem(nset,nmem)
873 member(nset) = nmem
874 return
875 end
876c=================================================