Removing obsolete constants.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / alphas.f
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
59 c      equivalence (Mc,Massc)
60 c      equivalence (Mb,Massb)
61 c      equivalence (Mt,Masst)
62
63       call setnset(nset)
64 c
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)
70 c           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)
102 c 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
160 c      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 
174 c
175 c       print *,'calling alphamrs with ',nflav,Q
176        call getnset(nset)
177 c       
178 c
179 c      alambda=0.323
180 c      alambda=0.3342
181 c -- 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
208 c - alambda found  -- save it !!!!
209 c - 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/
236 c  The value of Lambda required corresponds to nflav=4
237 c  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
240 c      print *,'calling mrslambda with ',nflav,Q,alambda
241       al2=alambda*alambda
242       q2=q*q
243       t=dlog(q2/al2)
244 c
245 c next line to be replaced by a getiord in the final version
246 c      iord=1
247 c
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    
266 c      if(QS.lt.QSDTT.and.nflav.eq.4) then
267 c        flav=3.
268 c       go to 11
269 c      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
298 c      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
402 c=====================================================================
403 c next is alphas routine from PDFLIB
404 c
405 c      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
410 c      CHARACTER*20 PARM(NCHDIM)
411 c      double precision 
412 c     +       VAL(NCHDIM)
413 c      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      
417 c      QCDL5=0.165d0
418       
419 c      print *,qcdl5,iord
420       LO = 1
421       if(iord.ne.0) LO = 2
422 c      if(LO.eq.1) then 
423 c         print *,'Calculating ALPHA_S at Leading Order',LO
424 c      else
425 c         print *,'Calculating ALPHA_S at Next to Leading Order',LO
426 c      endif
427
428       TMAS = 180.0d0
429     
430
431       ALPHAS2 = ZEROD
432 C.
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
454 C
455       Q   = SCALE
456       XLQ = TWOD *  LOG( Q/QCDL5 )
457       XLLQ =  LOG( XLQ )
458 c      KNF = NFL
459 c      NF = KNF
460      
461 c      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
471 c      ENDIF
472 c      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
493 c ========================================================================
494       SUBROUTINE CtLhAlphaNewSET(MC,MB,MT,Q0,ALPHA0,IORDER,IMODE)
495 c call to set quark masses for alpha_s, and choose lambda or its 
496 c 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
528 c set artificial quark masses, if necessary, in alpha_s to enforce 
529 c 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
537 c       print *,'CtLhAlphaNewSET: imode=',imode,' Nfl=',Nfl             !*** temporary ***
538 c       print *,'CtLhAlphaNewSET: mc,mb,mt=', mc, mb, mt                !*** temporary ***
539 c       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            
565 c external maximum number of flavors -- typically, nfmax=5 so 
566 c 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)
648 C calculate ALPHAS FROM RGE GIVEN AS0 AT QS0.
649 C
650 C IORD=1: LEADING ORDER DEFINED BY 
651 C                 Q*d(alpha)/d(Q) = c1*alpha**2 
652 C
653 C IORD=2,IMODE=1: QCD NUM CHOICE DEFINED BY 
654 C                 Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
655 C
656 C IORD=2,IMODE=2: AD HOC ALTERNATIVE DEFINED BY 
657 C                 Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
658 C
659 C IORD=2,IMODE=3: TRADITIONAL CTEQ CHOICE DEFINED BY 
660 C                 ALPHA = c3*(1 - c4*log(L)/L)/L, WHERE L=log((Q/lambda)**2)
661 C
662 C c1 = -beta0/(2*pi)     where beta0 =  11. -  (2./3.)*nf
663 C c2 = -beta1/(8*pi**2)  where beta1 = 102. - (38./3.)*nf
664 C
665 C c3 = -2/c1
666 C c4 = -2*c2/c1**2
667 C
668 C
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
689 C 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
699 C QCD NUM CHOICE: Q*d(alpha)/d(Q) = c1*alpha**2 + c2*alpha**3
700       IF(IMODE .EQ. 1) THEN
701
702 c use Newton's method to solve the equation, instead of the 
703 c 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  
724 20       CONTINUE
725          CtLhALPHAR = ALFA1
726          RETURN
727
728 C AD HOC ALTERNATIVE: Q*d(alpha)/d(Q) = c1*alpha**2 / (1 - (c2/c1)*alpha)
729       ELSEIF(IMODE .EQ. 2) THEN
730
731 c first get a good starting point, to be sure Newton's method doesn't go 
732 c 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
759 C 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
781 c 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
790 c 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)
806 C 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
823 c      call GetLam4M(nset,mem,alpha0)
824 c      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
831 c **********************************************************
832 c the input data file should probably specify a very large 
833 c value for xmt, because we never allow top quark as a parton
834 c in the PDF fitting, so there are never more than 5 active 
835 c flavors.  Presume it is most consistent therefore to 
836 c also keep top quark loop out of running of alpha.
837 c **********************************************************
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
851 c=================================================
852       subroutine getnset(nset)
853       integer iset,nset
854       save iset
855       nset = iset      
856 c      print *,'getting nset:',nset
857       return
858 c      
859       entry setnset(nset)
860       iset = nset
861 c      print *,'setting nset:',nset
862       return      
863       end
864 c=================================================
865       subroutine getnmem(nset,nmem)
866       include 'parmsetup.inc'
867       integer nmem,nset,member(nmxset)
868       save member
869       nmem = member(nset)      
870       return
871 c      
872       entry setnmem(nset,nmem)
873       member(nset) = nmem
874       return      
875       end
876 c=================================================