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