~maddevelopers/mg5amcnlo/3.0.1

« back to all changes in this revision

Viewing changes to MG4_DECAY/alfas_functions.f

move ./decay to ./mg5decay; resolve unit tests (n.b. __init__ does not check keys of input dictionaries, followed last revision)

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
C
2
 
C-----------------------------------------------------------------------------
3
 
C
4
 
      double precision function alfa(alfa0,qsq )
5
 
C
6
 
C-----------------------------------------------------------------------------
7
 
C
8
 
C       This function returns the 1-loop value of alpha.
9
 
C
10
 
C       INPUT: 
11
 
C               qsq   = Q^2
12
 
C
13
 
C-----------------------------------------------------------------------------
14
 
C
15
 
      implicit none
16
 
      double precision  qsq,alfa0
17
 
c
18
 
c constants
19
 
c
20
 
      double precision  One, Three, Pi,zmass
21
 
      parameter( One = 1.0d0, Three = 3.0d0 )
22
 
      parameter( Pi = 3.14159265358979323846d0 )
23
 
      parameter( zmass = 91.188d0 )
24
 
cc
25
 
      alfa = alfa0 / ( 1.0d0 - alfa0*dlog( qsq/zmass**2 ) /Three /Pi )
26
 
ccc
27
 
      return
28
 
      end
29
 
 
30
 
C
31
 
C-----------------------------------------------------------------------------
32
 
C
33
 
      double precision function alfaw(alfaw0,qsq,nh )
34
 
C
35
 
C-----------------------------------------------------------------------------
36
 
C
37
 
C       This function returns the 1-loop value of alpha_w.
38
 
C
39
 
C       INPUT: 
40
 
C               qsq = Q^2
41
 
C               nh  = # of Higgs doublets
42
 
C
43
 
C-----------------------------------------------------------------------------
44
 
C
45
 
      implicit none
46
 
      double precision  qsq, alphaw, dum,alfaw0
47
 
      integer  nh, nq
48
 
c
49
 
c         include
50
 
c
51
 
          
52
 
c
53
 
c constants
54
 
c
55
 
      double precision  Two, Four, Pi, Twpi, zmass,tmass
56
 
      parameter( Two = 2.0d0, Four = 4.0d0 )
57
 
      parameter( Pi = 3.14159265358979323846d0 )
58
 
      parameter( Twpi = 3.0d0*Four*Pi )
59
 
      parameter( zmass = 91.188d0,tmass=174d0 )
60
 
cc
61
 
      if ( qsq.ge.tmass**2 ) then
62
 
         nq = 6
63
 
      else
64
 
         nq = 5
65
 
      end if
66
 
      dum = ( 22.0d0 - Four*nq - nh/Two ) / Twpi
67
 
      alfaw = alfaw0 / ( 1.0d0 + dum*alfaw0*dlog( qsq/zmass**2 ) )
68
 
ccc
69
 
      return
70
 
      end
71
 
 
72
 
C-----------------------------------------------------------------------------
73
 
C
74
 
      DOUBLE PRECISION FUNCTION ALPHAS(Q)
75
 
c
76
 
c     Evaluation of strong coupling constant alpha_S
77
 
c     Author: R.K. Ellis
78
 
c
79
 
c     q -- scale at which alpha_s is to be evaluated
80
 
c
81
 
c-- common block alfas.inc
82
 
c     asmz -- value of alpha_s at the mass of the Z-boson
83
 
c     nloop -- the number of loops (1,2, or 3) at which beta 
84
 
c
85
 
c     function is evaluated to determine running.
86
 
c     the values of the cmass and the bmass should be set
87
 
c     in common block qmass.
88
 
C-----------------------------------------------------------------------------
89
 
 
90
 
      IMPLICIT NONE
91
 
c
92
 
      include 'alfas.inc'
93
 
      DOUBLE PRECISION Q,T,AMZ0,AMB,AMC
94
 
      DOUBLE PRECISION AS_OUT
95
 
      INTEGER NLOOP0,NF3,NF4,NF5
96
 
      PARAMETER(NF5=5,NF4=4,NF3=3)
97
 
C
98
 
      REAL*8       CMASS,BMASS
99
 
      COMMON/QMASS/CMASS,BMASS
100
 
      DATA CMASS,BMASS/1.42D0,4.7D0/  ! HEAVY QUARK MASSES FOR THRESHOLDS
101
 
C
102
 
      REAL*8 ZMASS
103
 
      DATA ZMASS/91.188D0/
104
 
C
105
 
      SAVE AMZ0,NLOOP0,AMB,AMC
106
 
      DATA AMZ0,NLOOP0/0D0,0/
107
 
      IF (Q .LE. 0D0) THEN 
108
 
         WRITE(6,*) 'q .le. 0 in alphas'
109
 
         WRITE(6,*) 'q= ',Q
110
 
         STOP
111
 
      ENDIF
112
 
      IF (asmz .LE. 0D0) THEN 
113
 
         WRITE(6,*) 'asmz .le. 0 in alphas',asmz
114
 
c         WRITE(6,*) 'continue with asmz=0.1185'
115
 
         STOP
116
 
         asmz=0.1185D0
117
 
      ENDIF
118
 
      IF (CMASS .LE. 0.3D0) THEN 
119
 
         WRITE(6,*) 'cmass .le. 0.3GeV in alphas',CMASS
120
 
c         WRITE(6,*) 'continue with cmass=1.5GeV?'
121
 
         STOP
122
 
         CMASS=1.42D0
123
 
      ENDIF
124
 
      IF (BMASS .LE. 0D0) THEN 
125
 
         WRITE(6,*) 'bmass .le. 0 in alphas',BMASS
126
 
         WRITE(6,*) 'COMMON/QMASS/CMASS,BMASS'
127
 
         STOP
128
 
         BMASS=4.7D0
129
 
      ENDIF
130
 
c--- establish value of coupling at b- and c-mass and save
131
 
      IF ((asmz .NE. AMZ0) .OR. (NLOOP .NE. NLOOP0)) THEN
132
 
         AMZ0=asmz
133
 
         NLOOP0=NLOOP
134
 
         T=2D0*DLOG(BMASS/ZMASS)
135
 
         CALL NEWTON1(T,asmz,AMB,NLOOP,NF5)
136
 
         T=2D0*DLOG(CMASS/BMASS)
137
 
         CALL NEWTON1(T,AMB,AMC,NLOOP,NF4)
138
 
      ENDIF
139
 
 
140
 
c--- evaluate strong coupling at scale q
141
 
      IF (Q  .LT. BMASS) THEN
142
 
           IF (Q  .LT. CMASS) THEN
143
 
             T=2D0*DLOG(Q/CMASS)
144
 
             CALL NEWTON1(T,AMC,AS_OUT,NLOOP,NF3)
145
 
           ELSE
146
 
             T=2D0*DLOG(Q/BMASS)
147
 
             CALL NEWTON1(T,AMB,AS_OUT,NLOOP,NF4)
148
 
           ENDIF
149
 
      ELSE
150
 
      T=2D0*DLOG(Q/ZMASS)
151
 
      CALL NEWTON1(T,asmz,AS_OUT,NLOOP,NF5)
152
 
      ENDIF
153
 
      ALPHAS=AS_OUT
154
 
      RETURN
155
 
      END
156
 
 
157
 
 
158
 
      SUBROUTINE NEWTON1(T,A_IN,A_OUT,NLOOP,NF)
159
 
C     Author: R.K. Ellis
160
 
 
161
 
c---  calculate a_out using nloop beta-function evolution 
162
 
c---  with nf flavours, given starting value as-in
163
 
c---  given as_in and logarithmic separation between 
164
 
c---  input scale and output scale t.
165
 
c---  Evolution is performed using Newton's method,
166
 
c---  with a precision given by tol.
167
 
 
168
 
      IMPLICIT NONE
169
 
      INTEGER NLOOP,NF
170
 
      REAL*8 T,A_IN,A_OUT,AS,TOL,F2,F3,F,FP,DELTA
171
 
      REAL*8 B0(3:5),C1(3:5),C2(3:5),DEL(3:5)
172
 
      PARAMETER(TOL=5.D-4)
173
 
C---     B0=(11.-2.*NF/3.)/4./PI
174
 
      DATA B0/0.716197243913527D0,0.66314559621623D0,0.61009394851893D0/
175
 
C---     C1=(102.D0-38.D0/3.D0*NF)/4.D0/PI/(11.D0-2.D0/3.D0*NF)
176
 
      DATA C1/.565884242104515D0,0.49019722472304D0,0.40134724779695D0/
177
 
C---     C2=(2857.D0/2.D0-5033*NF/18.D0+325*NF**2/54)
178
 
C---     /16.D0/PI**2/(11.D0-2.D0/3.D0*NF)
179
 
      DATA C2/0.453013579178645D0,0.30879037953664D0,0.14942733137107D0/
180
 
C---     DEL=SQRT(4*C2-C1**2)
181
 
      DATA DEL/1.22140465909230D0,0.99743079911360D0,0.66077962451190D0/
182
 
      F2(AS)=1D0/AS+C1(NF)*LOG((C1(NF)*AS)/(1D0+C1(NF)*AS))
183
 
      F3(AS)=1D0/AS+0.5D0*C1(NF)
184
 
     & *LOG((C2(NF)*AS**2)/(1D0+C1(NF)*AS+C2(NF)*AS**2))
185
 
     & -(C1(NF)**2-2D0*C2(NF))/DEL(NF)
186
 
     & *ATAN((2D0*C2(NF)*AS+C1(NF))/DEL(NF))
187
 
 
188
 
           
189
 
      A_OUT=A_IN/(1D0+A_IN*B0(NF)*T)
190
 
      IF (NLOOP .EQ. 1) RETURN
191
 
      A_OUT=A_IN/(1D0+B0(NF)*A_IN*T+C1(NF)*A_IN*LOG(1D0+A_IN*B0(NF)*T))
192
 
      IF (A_OUT .LT. 0D0) AS=0.3D0
193
 
 30   AS=A_OUT
194
 
 
195
 
      IF (NLOOP .EQ. 2) THEN
196
 
      F=B0(NF)*T+F2(A_IN)-F2(AS)
197
 
      FP=1D0/(AS**2*(1D0+C1(NF)*AS))
198
 
      ENDIF
199
 
      IF (NLOOP .EQ. 3) THEN
200
 
      F=B0(NF)*T+F3(A_IN)-F3(AS)
201
 
      FP=1D0/(AS**2*(1D0+C1(NF)*AS+C2(NF)*AS**2))
202
 
      ENDIF
203
 
      A_OUT=AS-F/FP
204
 
      DELTA=ABS(F/FP/AS)
205
 
      IF (DELTA .GT. TOL) GO TO 30
206
 
      RETURN
207
 
      END
208
 
 
209
 
 
210
 
C-----------------------------------------------------------------------------
211
 
C
212
 
      double precision function mfrun(mf,scale,asmz,nloop)
213
 
C
214
 
C-----------------------------------------------------------------------------
215
 
C
216
 
C       This function returns the 2-loop value of a MSbar fermion mass
217
 
C       at a given scale.
218
 
C
219
 
C       INPUT: mf    = MSbar mass of fermion at MSbar fermion mass scale 
220
 
C              scale = scale at which the running mass is evaluated
221
 
C              asmz  = AS(MZ) : this is passed to alphas(scale,asmz,nloop)
222
 
C              nloop = # of loops in the evolution
223
 
C       
224
 
C
225
 
C
226
 
C       EXTERNAL:      double precision alphas(scale,asmz,nloop)
227
 
C                      
228
 
C-----------------------------------------------------------------------------
229
 
C
230
 
      implicit none
231
 
C
232
 
C     ARGUMENTS
233
 
C
234
 
      double precision  mf,scale,asmz
235
 
      integer           nloop
236
 
C
237
 
C     LOCAL
238
 
C
239
 
      double precision  beta0, beta1,gamma0,gamma1
240
 
      double precision  A1,as,asmf,l2
241
 
      integer  nf
242
 
C
243
 
C     EXTERNAL
244
 
C
245
 
      double precision  alphas
246
 
      external          alphas
247
 
c
248
 
c     CONSTANTS
249
 
c
250
 
      double precision  One, Two, Three, Pi
251
 
      parameter( One = 1.0d0, Two = 2.0d0, Three = 3.0d0 )
252
 
      parameter( Pi = 3.14159265358979323846d0) 
253
 
      double precision tmass
254
 
      parameter(tmass=174d0)
255
 
cc
256
 
C
257
 
C
258
 
      if ( mf.gt.tmass ) then
259
 
         nf = 6
260
 
      else
261
 
         nf = 5
262
 
      end if
263
 
 
264
 
      beta0 = ( 11.0d0 - Two/Three *nf )/4d0
265
 
      beta1 = ( 102d0  - 38d0/Three*nf )/16d0
266
 
      gamma0= 1d0
267
 
      gamma1= ( 202d0/3d0  - 20d0/9d0*nf )/16d0
268
 
      A1    = -beta1*gamma0/beta0**2+gamma1/beta0
269
 
      as    = alphas(scale)
270
 
      asmf  = alphas(mf)
271
 
      l2    = (1+ A1*as/Pi)/(1+ A1*asmf/Pi)
272
 
      
273
 
      
274
 
      mfrun = mf * (as/asmf)**(gamma0/beta0)
275
 
 
276
 
      if(nloop.eq.2) mfrun =mfrun*l2
277
 
ccc
278
 
      return
279
 
      end
280