← Back to branch summary

~mg5core1/mg5amcnlo/2.6.4

« back to all changes in this revision

Viewing changes to MG4_DECAY/hdecay.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:
MG4_DECAY/hdecay.f
1
 
C         Last modification on March 10th 2001 by M.S.
2
 
C ==================================================================
3
 
C ================= PROGRAM HDECAY: COMMENTS =======================
4
 
C ==================================================================
5
 
C
6
 
C                       ***************
7
 
C                       * VERSION 2.0 *
8
 
C                       ***************
9
 
C
10
 
C
11
 
C  This program calculates the total decay widths and the branching 
12
 
C  ratios of the C Standard Model Higgs boson (HSM) as well as those 
13
 
C  of the neutral (HL= the light CP-even, HH= the heavy CP-even, HA= 
14
 
C  the pseudoscalar) and the charged (HC) Higgs bosons of the Minimal
15
 
C  Supersymmetric extension of the Standard Model (MSSM). It includes:
16
 
C
17
 
C - All the decay channels which are kinematically allowed and which
18
 
C   have branching ratios larger than 10**(-4). 
19
 
C
20
 
C - All QCD corrections to the fermionic and gluonic decay modes.
21
 
C   Most of these corrections are mapped into running masses in a
22
 
C   consistent way with some freedom for including high order terms. 
23
 
C
24
 
C - Below--threshold three--body decays with off--shell top quarks
25
 
C   or ONE off-shell gauge boson, as well as some decays with one
26
 
C   off-shell Higgs boson in the MSSM. 
27
 
C
28
 
C - Double off-shell decays: HSM,HL,HH --> W*W*,Z*Z* -->4 fermions,
29
 
C   which could be important for Higgs masses close to MW or MZ.
30
 
C
31
 
C - In the MSSM, the radiative corrections with full squark mixing and 
32
 
C   uses the RG improved values of Higgs masses and couplings with the 
33
 
C   main NLO corrections implemented (based on M.Carena, M. Quiros and
34
 
C   C.E.M. Wagner, Nucl. Phys. B461 (1996) 407, hep-ph/9508343). 
35
 
C
36
 
C - In the MSSM, all the decays into CHARGINOS, NEUTRALINOS, SLEPTONS 
37
 
C   and SQUARKS (with mixing in the stop and sbottom sectors). 
38
 
C
39
 
C - Chargino, slepton and squark loops in the 2 photon decays and squark
40
 
C   loops in the gluonic decays (including QCD corrections). 
41
 
C
42
 
C  ===================================================================
43
 
C  This program has been written by A.Djouadi, J.Kalinowski and M.Spira.
44
 
C  For details on how to use the program see: Comp. Phys. Commun. 108
45
 
C  (1998) 56, hep-ph/9704448. For any question, comment, suggestion or
46
 
C  complaint, please contact us at:
47
 
C          djouadi@lpm.univ-montp2.fr
48
 
C          kalino@fuw.edu.pl
49
 
C          Michael.Spira@cern.ch
50
 
 
51
 
 
52
 
C ================ IT USES AS INPUT PARAMETERS:
53
 
C
54
 
C   IHIGGS: =0: CALCULATE BRANCHING RATIOS OF SM HIGGS BOSON
55
 
C           =1: CALCULATE BRANCHING RATIOS OF MSSM h BOSON
56
 
C           =2: CALCULATE BRANCHING RATIOS OF MSSM H BOSON
57
 
C           =3: CALCULATE BRANCHING RATIOS OF MSSM A BOSON
58
 
C           =4: CALCULATE BRANCHING RATIOS OF MSSM H+ BOSON
59
 
C           =5: CALCULATE BRANCHING RATIOS OF ALL MSSM HIGGS BOSONS
60
 
C
61
 
C TGBET:    TAN(BETA) FOR MSSM
62
 
C MABEG:    START VALUE OF M_A FOR MSSM AND M_H FOR SM
63
 
C MAEND:    END VALUE OF M_A FOR MSSM AND M_H FOR SM
64
 
C NMA:      NUMBER OF ITERATIONS FOR M_A
65
 
C ALS(MZ):  VALUE FOR ALPHA_S(M_Z)
66
 
C MSBAR(1): MSBAR MASS OF STRANGE QUARK AT SCALE Q=1 GEV
67
 
C MC:       CHARM POLE MASS
68
 
C MB:       BOTTOM POLE MASS
69
 
C MT:       TOP POLE MASS
70
 
C MTAU:     TAU MASS
71
 
C MMUON:    MUON MASS
72
 
C ALPH:     INVERSE QED COUPLING
73
 
C GF:       FERMI CONSTANT
74
 
C GAMW:     W WIDTH
75
 
C GAMZ:     Z WIDTH
76
 
C MZ:       Z MASS
77
 
C MW:       W MASS
78
 
C VUS:      CKM PARAMETER V_US
79
 
C VCB:      CKM PARAMETER V_CB
80
 
C VUB/VCB:  RATIO V_UB/V_CB
81
 
C 1ST AND 2ND GENERATION:
82
 
C MSL1:      SUSY BREAKING MASS PARAMETERS OF LEFT HANDED SLEPTONS 
83
 
C MER1:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED SLEPTONS 
84
 
C MQL1:      SUSY BREAKING MASS PARAMETERS OF LEFT HANDED SUPS
85
 
C MUR1:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED SUPS
86
 
C MDR1:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED SDOWNS 
87
 
C 3RD GENERATION:
88
 
C MSL:      SUSY BREAKING MASS PARAMETERS OF LEFT HANDED STAUS 
89
 
C MER:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED STAUS 
90
 
C MSQ:      SUSY BREAKING MASS PARAMETERS OF LEFT HANDED STOPS
91
 
C MUR:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED STOPS
92
 
C MDR:      SUSY BREAKING MASS PARAMETERS OF RIGHT HANDED SBOTTOMS 
93
 
C AL:       STAU TRILINEAR SOFT BREAKING TERMS 
94
 
C AU:       STOP TRILINEAR SOFT BREAKING TERMS.
95
 
C AD:       SBOTTOM TRILINEAR SOFT BREAKING TERMS.
96
 
C MU:       SUSY HIGGS MASS PARAMETER
97
 
C M2:       GAUGINO MASS PARAMETER. 
98
 
C
99
 
C NNLO (M): =0: USE O(ALPHA_S) FORMULA FOR POLE MASS --> MSBAR MASS
100
 
C           =1: USE O(ALPHA_S**2) FORMULA FOR POLE MASS --> MSBAR MASS
101
 
C
102
 
C ON-SHELL: =0: INCLUDE OFF_SHELL DECAYS H,A --> T*T*, A --> Z*H,
103
 
C               H --> W*H+,Z*A, H+ --> W*A, W*H, T*B
104
 
C           =1: EXCLUDE THE OFF-SHELL DECAYS ABOVE
105
 
C
106
 
C ON-SH-WZ: =0: INCLUDE DOUBLE OFF-SHELL PAIR DECAYS PHI --> W*W*,Z*Z*
107
 
C           =1: INCLUDE ONLY SINGLE OFF-SHELL DECAYS PHI --> W*W,Z*Z
108
 
C
109
 
C IPOLE:    =0 COMPUTES RUNNING HIGGS MASSES (FASTER) 
110
 
C           =1 COMPUTES POLE HIGGS MASSES 
111
 
C
112
 
C OFF-SUSY: =0: INCLUDE DECAYS (AND LOOPS) INTO SUPERSYMMETRIC PARTICLES
113
 
C           =1: EXCLUDE DECAYS (AND LOOPS) INTO SUPERSYMMETRIC PARTICLES
114
 
C
115
 
C INIDEC:   =0: PRINT OUT SUMS OF CHARGINO/NEUTRALINO/SFERMION DECAYS
116
 
C           =1: PRINT OUT INDIVIDUAL CHARGINO/NEUTRALINO/SFERMION DECAYS
117
 
C
118
 
C NF-GG:    NUMBER OF LIGHT FLAVORS INCLUDED IN THE GLUONIC DECAYS 
119
 
C            PHI --> GG* --> GQQ (3,4 OR 5)
120
 
C           
121
 
C =======================================================================
122
 
C ============== BEGINNING OF THE MAIN PROGRAM ==========================
123
 
C =======================================================================
124
 
C
125
 
C      PROGRAM HDECAY
126
 
      SUBROUTINE HDECAY
127
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
128
 
      PARAMETER(K=6,NI=87,NSA=85,NSB=86,NLA=88,NLB=89,NHA=90,NHB=91,
129
 
     .          NHC=92,NAA=93,NAB=94,NCA=95,NCB=96,NRA=97,NRB=98,
130
 
     .          NSUSYL=81,NSUSYA=82,NSUSYH=83,NSUSYC=84,NPAR=80,
131
 
     .          NSUSYLA=79,NSUSYLB=78,NSUSYLC=77,NSUSYLD=76,NSUSYLE=75,
132
 
     .          NSUSYLF=59,NSUSYHF=58,
133
 
     .          NSUSYHA=74,NSUSYHB=73,NSUSYHC=72,NSUSYHD=71,NSUSYHE=70,
134
 
     .          NSUSYAA=69,NSUSYAB=68,NSUSYAC=67,NSUSYAD=66,NSUSYAE=65,
135
 
     .          NSUSYCA=64,NSUSYCB=63,NSUSYCC=62,NSUSYCD=61,NSUSYCE=60)
136
 
      DIMENSION GMN(4),XMN(4),GMC(2),GMST(2),GMSB(2),GMSL(2),
137
 
     .          GMSU(2),GMSD(2),GMSE(2),GMSN(2)
138
 
      DIMENSION HLBRSC(2,2),HLBRSN(4,4),HHBRSC(2,2),HHBRSN(4,4),
139
 
     .          HABRSC(2,2),HABRSN(4,4),HCBRSU(2,4),
140
 
     .          HHBRST(2,2),HHBRSB(2,2),HCBRSTB(2,2) 
141
 
      DIMENSION AC1(2,2),AC2(2,2),AC3(2,2),
142
 
     .          AN1(4,4),AN2(4,4),AN3(4,4),
143
 
     .          ACNL(2,4),ACNR(2,4)
144
 
      DIMENSION GLTT(2,2),GLBB(2,2),GHTT(2,2),GHBB(2,2),GCTB(2,2),
145
 
     .          GLEE(2,2),GHEE(2,2),GCEN(2,2)
146
 
      COMMON/MASSES/AMS,AMC,AMB,AMT
147
 
      COMMON/STRANGE/AMSB
148
 
      COMMON/PARAM/GF,ALPH,AMTAU,AMMUON,AMZ,AMW
149
 
      COMMON/CKMPAR/VUS,VCB,VUB
150
 
      COMMON/HMASS/AMSM,AMA,AML,AMH,AMCH,AMAR
151
 
      COMMON/BREAK/AMEL,AMER,AMSQ,AMUR,AMDR,AL,AU,AD,AMU,AM2
152
 
      COMMON/SFER1ST/AMQL1,AMUR1,AMDR1,AMEL1,AMER1
153
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW,GAMZ
154
 
      COMMON/COUP/GAT,GAB,GLT,GLB,GHT,GHB,GZAH,GZAL,
155
 
     .            GHHH,GLLL,GHLL,GLHH,GHAA,GLAA,GLVV,GHVV,
156
 
     .            GLPM,GHPM,B,A
157
 
      COMMON/ALS/XLAMBDA,AMC0,AMB0,AMT0,N0
158
 
      COMMON/FLAG/IHIGGS,NNLO,IPOLE
159
 
      COMMON/ONSHELL/IONSH,IONWZ,IOFSUSY
160
 
      COMMON/OLDFASH/NFGG
161
 
      COMMON/WIDTHSM/SMBRB,SMBRL,SMBRM,SMBRS,SMBRC,SMBRT,SMBRG,SMBRGA,
162
 
     .               SMBRZGA,SMBRW,SMBRZ,SMWDTH
163
 
      COMMON/WIDTHA/ABRB,ABRL,ABRM,ABRS,ABRC,ABRT,ABRG,ABRGA,ABRZGA,
164
 
     .              ABRZ,AWDTH
165
 
      COMMON/WIDTHHL/HLBRB,HLBRL,HLBRM,HLBRS,HLBRC,HLBRT,HLBRG,HLBRGA,
166
 
     .               HLBRZGA,HLBRW,HLBRZ,HLBRA,HLBRAZ,HLBRHW,HLWDTH
167
 
      COMMON/WIDTHHH/HHBRB,HHBRL,HHBRM,HHBRS,HHBRC,HHBRT,HHBRG,HHBRGA,
168
 
     .               HHBRZGA,HHBRW,HHBRZ,HHBRH,HHBRA,HHBRAZ,HHBRHW,
169
 
     .               HHWDTH
170
 
      COMMON/WIDTHHC/HCBRB,HCBRL,HCBRM,HCBRBU,HCBRS,HCBRC,HCBRT,HCBRW,
171
 
     .               HCBRA,HCWDTH
172
 
      COMMON/WISUSY/HLBRSC,HLBRSN,HHBRSC,HHBRSN,HABRSC,HABRSN,HCBRSU,
173
 
     .              HLBRCHT,HHBRCHT,HABRCHT,HLBRNET,HHBRNET,HABRNET,
174
 
     .              HCBRCNT,HLBRSL,HHBRSL,HCBRSL,HABRSL,HABRST,HABRSB,
175
 
     .              HHBRSQ,HHBRST,HHBRSB,HHBRSQT,HCBRSQ,HCBRSTB,
176
 
     .              HCBRSQT,HLBRSQ,HLBRSQT
177
 
      COMMON/WISFER/BHLSLNL,BHLSLEL,BHLSLER,BHLSQUL,BHLSQUR,BHLSQDL,
178
 
     .              BHLSQDR,BHLST(2,2),BHLSB(2,2),BHLSTAU(2,2),
179
 
     .              BHHSLNL,BHHSLEL,BHHSLER,BHHSQUL,BHHSQUR,BHHSQDL,
180
 
     .              BHHSQDR,BHHST(2,2),BHHSB(2,2),BHHSTAU(2,2),
181
 
     .              BHASTAU,BHASB,BHAST,
182
 
     .              BHCSL00,BHCSL11,BHCSL21,BHCSQ,BHCSTB(2,2)
183
 
      COMMON/SMASS/GMN,XMN,GMC,GMST,GMSB,GMSL,GMSU,GMSD,GMSE,GMSN 
184
 
c
185
 
      double precision  hdals, hdmhbeg, hdmhend, tgbet
186
 
      common /hdparms/  hdals, hdmhbeg, hdmhend, tgbet
187
 
c
188
 
      PI = 4*DATAN(1D0)
189
 
 
190
 
      ALSMZ=HDALS
191
 
      AMABEG=hdmhbeg
192
 
      AMAEND=hdmhend
193
 
 
194
 
c      OPEN(NI,FILE='hdecay.in')
195
 
c      OPEN(NPAR,FILE='br.input')
196
 
 
197
 
 
198
 
c      READ(NI,101)IHIGGS
199
 
c      READ(NI,100)TGBET
200
 
c      READ(NI,100)AMABEG
201
 
c      READ(NI,100)AMAEND
202
 
c      READ(NI,101)NMA
203
 
c      READ(NI,100)ALSMZ
204
 
c      READ(NI,100)AMS
205
 
c      READ(NI,100)AMC
206
 
c      READ(NI,100)AMB
207
 
c      READ(NI,100)AMT
208
 
c      READ(NI,100)AMTAU
209
 
c      READ(NI,100)AMMUON
210
 
c      READ(NI,100)ALPH
211
 
c      READ(NI,100)GF
212
 
c      READ(NI,100)GAMW
213
 
c      READ(NI,100)GAMZ
214
 
c      READ(NI,100)AMZ
215
 
c      READ(NI,100)AMW
216
 
c      READ(NI,100)VUS
217
 
c      READ(NI,100)VCB
218
 
c      READ(NI,100)RVUB
219
 
c      READ(NI,100)AMU
220
 
c      READ(NI,100)AM2
221
 
c      READ(NI,100)AMEL1
222
 
c      READ(NI,100)AMER1
223
 
c      READ(NI,100)AMQL1
224
 
c      READ(NI,100)AMUR1
225
 
c      READ(NI,100)AMDR1
226
 
c      READ(NI,100)AMEL
227
 
c      READ(NI,100)AMER
228
 
c      READ(NI,100)AMSQ
229
 
c      READ(NI,100)AMUR
230
 
c      READ(NI,100)AMDR
231
 
c      READ(NI,100)AL
232
 
c      READ(NI,100)AU
233
 
c      READ(NI,100)AD
234
 
c      READ(NI,101)NNLO
235
 
c      READ(NI,101)IONSH
236
 
c      READ(NI,101)IONWZ
237
 
c      READ(NI,101)IPOLE
238
 
c      READ(NI,101)IOFSUSY
239
 
c      READ(NI,101)INDIDEC
240
 
      TGBET    = 10.D0
241
 
      NMA      = 1
242
 
      VUS      = 0.2205D0
243
 
      VCB      = 0.04D0
244
 
      RVUB     = 0.08D0
245
 
      AMU       = -200.D0
246
 
      AM2       = 200.D0
247
 
      AMEL1     = 500.D0
248
 
      AMER1     = 500.D0
249
 
      AMQL1     = 500.D0
250
 
      AMUR1     = 500.D0
251
 
      AMDR1     = 500.D0
252
 
      AMEL      = 500.D0
253
 
      AMER      = 500.D0
254
 
      AMSQ      = 1000.D0
255
 
      AMUR      = 1000.D0
256
 
      AMDR      = 1000.D0
257
 
      AL        = 1000.D0
258
 
      AU        = 1000.D0
259
 
      AD        = 1000.D0
260
 
      NFGG      = 4
261
 
       
262
 
 
263
 
      VUB=RVUB*VCB
264
 
      ALPH=1.D0/ALPH
265
 
      AMSB = AMS
266
 
 
267
 
 
268
 
 
269
 
C--CHECK NFGG
270
 
      IF(NFGG.GT.5.OR.NFGG.LT.3)THEN
271
 
C       WRITE(6,*)'NF-GG NOT VALID. TAKING THE DEFAULT NF-GG = 3....'
272
 
       NFGG = 3
273
 
      ENDIF
274
 
 
275
 
100   FORMAT(10X,G30.20)
276
 
101   FORMAT(10X,I30)
277
 
 
278
 
      AMC0=AMC
279
 
      AMB0=AMB
280
 
      AMT0=AMT
281
 
      ACC=1.D-8
282
 
      NLOOP=2
283
 
      XLAMBDA=XITLA(NLOOP,ALSMZ,ACC)
284
 
      N0=5
285
 
      CALL ALSINI(ACC)
286
 
  
287
 
C--INITIALIZE COEFFICIENTS FOR POLYLOGARITHMS
288
 
      NBER = 18
289
 
      CALL BERNINI(NBER)
290
 
 
291
 
 
292
 
      DO 9999 II=1,NMA
293
 
       IF(NMA.NE.1)THEN
294
 
        AMAR = AMABEG + (AMAEND-AMABEG)/(NMA-1D0)*(II-1D0)
295
 
       ELSE
296
 
        AMAR = AMABEG
297
 
       ENDIF
298
 
       AMSM = AMAR
299
 
       AMA = AMAR
300
 
       
301
 
***********************************************
302
 
 
303
 
      CALL HDEC(TGBET)
304
 
      
305
 
 9999 continue
306
 
 
307
 
      RETURN
308
 
      END
309
 
 
310
 
C =====================================================================
311
 
C =========== BEGINNING OF THE SUBROUTINE FOR THE DECAYS ==============
312
 
C !!!!!!!!!!!!!! Any change below this line is at your own risk!!!!!!!!
313
 
C =====================================================================
314
 
 
315
 
      SUBROUTINE HDEC(TGBET)
316
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
317
 
      DOUBLE PRECISION LAMB
318
 
      DIMENSION XX(4),YY(4)
319
 
      DIMENSION AMCHAR(2),AMNEUT(4),XMNEUT(4),
320
 
     .          AC1(2,2),AC2(2,2),AC3(2,2),
321
 
     .          AN1(4,4),AN2(4,4),AN3(4,4),
322
 
     .          ACNL(2,4),ACNR(2,4),
323
 
     .          AMST(2),AMSB(2),AMSL(2),
324
 
     .          AMSU(2),AMSD(2),AMSE(2),AMSN(2),
325
 
     .          GLTT(2,2),GLBB(2,2),GLEE(2,2),
326
 
     .          GHTT(2,2),GHBB(2,2),GHEE(2,2),
327
 
     .          GCTB(2,2),GCEN(2,2)
328
 
      DIMENSION GMST(2),GMSB(2),GMSL(2),GMSU(2),GMSD(2),GMSE(2),
329
 
     .          GMSN(2)
330
 
      DIMENSION HLBRSC(2,2),HLBRSN(4,4),HHBRSC(2,2),
331
 
     .          HHBRSN(4,4),HABRSC(2,2),HABRSN(4,4),HCBRSU(2,4),
332
 
     .          HHBRST(2,2),HHBRSB(2,2),HCBRSTB(2,2) 
333
 
      DIMENSION WHLCH(2,2),WHLNE(4,4),WHHCH(2,2),WHHNE(4,4),
334
 
     .          WHACH(2,2),WHANE(4,4),WHCCN(2,4),
335
 
     .          WHHST(2,2),WHHSB(2,2),WHHSTAU(2,2),WHCSTB(2,2), 
336
 
     .          WHLST(2,2),WHLSB(2,2),WHLSTAU(2,2)
337
 
      COMPLEX*16 CF,CG,CI1,CI2,CA,CB,CTT,CTB,CTC,CTW,CLT,CLB,CLW,
338
 
     .           CAT,CAB,CAC,CAW,CAH,CTH,CLH,CX1,CX2,CAX1,CAX2,CTL,CAL,
339
 
     .           CSL,CSQ,CSB1,CSB2,CST1,CST2,CSL1,CSL2,
340
 
     .           CXL,CXQ,CXB1,CXB2,CXT1,CXT2,CXL1,CXL2
341
 
      COMPLEX*16 CSEL,CSER,CSUL,CSUR,CSDL,CSDR,
342
 
     .           CXEL,CXER,CXUL,CXUR,CXDL,CXDR
343
 
      COMMON/HMASS/AMSM,AMA,AML,AMH,AMCH,AMAR
344
 
      COMMON/CHIMASS/AMCHI
345
 
      COMMON/MASSES/AMS,AMC,AMB,AMT
346
 
      COMMON/ALS/XLAMBDA,AMC0,AMB0,AMT0,N0
347
 
      COMMON/PARAM/GF,ALPH,AMTAU,AMMUON,AMZ,AMW
348
 
      COMMON/CKMPAR/VUS,VCB,VUB
349
 
      COMMON/BREAK/AMEL,AMER,AMSQ,AMUR,AMDR,AL,AU,AD,AMU,AM2
350
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW,GAMZ
351
 
      COMMON/ONSHELL/IONSH,IONWZ,IOFSUSY
352
 
      COMMON/OLDFASH/NFGG
353
 
      COMMON/FLAG/IHIGGS,NNLO,IPOLE
354
 
      COMMON/WIDTHSM/SMBRB,SMBRL,SMBRM,SMBRS,SMBRC,SMBRT,SMBRG,SMBRGA,
355
 
     .               SMBRZGA,SMBRW,SMBRZ,SMWDTH
356
 
      COMMON/WIDTHA/ABRB,ABRL,ABRM,ABRS,ABRC,ABRT,ABRG,ABRGA,ABRZGA,
357
 
     .              ABRZ,AWDTH
358
 
      COMMON/WIDTHHL/HLBRB,HLBRL,HLBRM,HLBRS,HLBRC,HLBRT,HLBRG,HLBRGA,
359
 
     .               HLBRZGA,HLBRW,HLBRZ,HLBRA,HLBRAZ,HLBRHW,HLWDTH
360
 
      COMMON/WIDTHHH/HHBRB,HHBRL,HHBRM,HHBRS,HHBRC,HHBRT,HHBRG,HHBRGA,
361
 
     .               HHBRZGA,HHBRW,HHBRZ,HHBRH,HHBRA,HHBRAZ,HHBRHW,
362
 
     .               HHWDTH
363
 
      COMMON/WIDTHHC/HCBRB,HCBRL,HCBRM,HCBRBU,HCBRS,HCBRC,HCBRT,HCBRW,
364
 
     .               HCBRA,HCWDTH
365
 
      COMMON/WISUSY/HLBRSC,HLBRSN,HHBRSC,HHBRSN,HABRSC,HABRSN,HCBRSU,
366
 
     .              HLBRCHT,HHBRCHT,HABRCHT,HLBRNET,HHBRNET,HABRNET,
367
 
     .              HCBRCNT,HLBRSL,HHBRSL,HCBRSL,HABRSL,HABRST,HABRSB,
368
 
     .              HHBRSQ,HHBRST,HHBRSB,HHBRSQT,HCBRSQ,HCBRSTB,
369
 
     .              HCBRSQT,HLBRSQ,HLBRSQT
370
 
      COMMON/WISFER/BHLSLNL,BHLSLEL,BHLSLER,BHLSQUL,BHLSQUR,BHLSQDL,
371
 
     .              BHLSQDR,BHLST(2,2),BHLSB(2,2),BHLSTAU(2,2),
372
 
     .              BHHSLNL,BHHSLEL,BHHSLER,BHHSQUL,BHHSQUR,BHHSQDL,
373
 
     .              BHHSQDR,BHHST(2,2),BHHSB(2,2),BHHSTAU(2,2),
374
 
     .              BHASTAU,BHASB,BHAST,
375
 
     .              BHCSL00,BHCSL11,BHCSL21,BHCSQ,BHCSTB(2,2)
376
 
      COMMON/SMASS/AMNEUT,XMNEUT,AMCHAR,AMST,AMSB,AMSL,
377
 
     .              AMSU,AMSD,AMSE,AMSN 
378
 
      COMMON/COUP/GAT,GAB,GLT,GLB,GHT,GHB,GZAH,GZAL,
379
 
     .            GHHH,GLLL,GHLL,GLHH,GHAA,GLAA,GLVV,GHVV,
380
 
     .            GLPM,GHPM,B,A
381
 
      HVV(X,Y)= GF/(4.D0*PI*DSQRT(2.D0))*X**3/2.D0*BETA(Y)
382
 
     .            *(1.D0-4.D0*Y+12.D0*Y**2)
383
 
      AFF(X,Y)= GF/(4*PI*DSQRT(2.D0))*X**3*Y*(BETA(Y))
384
 
      HFF(X,Y)= GF/(4*PI*DSQRT(2.D0))*X**3*Y*(BETA(Y))**3
385
 
      CFF(Z,TB,X,Y)= GF/(4*PI*DSQRT(2.D0))*Z**3*LAMB(X,Y)
386
 
     .              *((1.D0-X-Y)*(X*TB**2+Y/TB**2)-4.D0*X*Y)
387
 
      HV(V)=3.D0*(1.D0-8.D0*V+20.D0*V**2)/DSQRT((4.D0*V-1.D0))
388
 
     .      *DACOS((3.D0*V-1.D0)/2.D0/DSQRT(V**3))
389
 
     .      -(1.D0-V)*(47.D0/2.D0*V-13.D0/2.D0+1.D0/V)
390
 
     .      -3.D0/2.D0*(1.D0-6.D0*V+4.D0*V**2)*DLOG(V)
391
 
      HVH(X,Y)=0.25D0*( (1-X)*(-2+4*X-2*X**2+9*Y+9*X*Y-6*Y**2)
392
 
     .        /(3*Y)-2*(1-X-X**2+X**3-3*Y-2*X*Y-3*X**2*Y+3*Y**2
393
 
     .        +3*X*Y**2-Y**3)*(-PI/2- DATAN((1-2*X+X**2-Y-X*Y)/
394
 
     .         ((1-X)*DSQRT(-1.D0+2*X+2*Y-(X-Y)**2))))/DSQRT(-1.D0
395
 
     .         +2*X-(X-Y)**2+2*Y)-(1+X**2-2*Y-2*X*Y+Y**2)*DLOG(X))
396
 
      QCD0(X) = (1+X**2)*(4*SP((1-X)/(1+X)) + 2*SP((X-1)/(X+1))
397
 
     .        - 3*DLOG((1+X)/(1-X))*DLOG(2/(1+X))
398
 
     .        - 2*DLOG((1+X)/(1-X))*DLOG(X))
399
 
     .        - 3*X*DLOG(4/(1-X**2)) - 4*X*DLOG(X)
400
 
      HQCDM(X)=QCD0(X)/X+(3+34*X**2-13*X**4)/16/X**3*DLOG((1+X)/(1-X))
401
 
     .        + 3.D0/8/X**2*(7*X**2-1)
402
 
      AQCDM(X)=QCD0(X)/X + (19+2*X**2+3*X**4)/16/X*DLOG((1+X)/(1-X))
403
 
     .        + 3.D0/8*(7-X**2)
404
 
      HQCD(X)=(4.D0/3*HQCDM(BETA(X))
405
 
     .        +2*(4.D0/3-DLOG(X))*(1-10*X)/(1-4*X))*ASH/PI
406
 
     .       + (29.14671D0 + RATCOUP*(1.570D0 - 2*DLOG(HIGTOP)/3
407
 
     .                                     + DLOG(X)**2/9))*(ASH/PI)**2
408
 
     .       + (164.14D0 - 25.77D0*5 + 0.259D0*5**2)*(ASH/PI)**3
409
 
      AQCD(X)=(4.D0/3*AQCDM(BETA(X))
410
 
     .        +2*(4.D0/3-DLOG(X))*(1-6*X)/(1-4*X))*ASH/PI
411
 
     .       + (29.14671D0 + RATCOUP*(23/6.D0 - DLOG(HIGTOP)
412
 
     .                                     + DLOG(X)**2/6))*(ASH/PI)**2
413
 
     .       + (164.14D0 - 25.77D0*5 + 0.259D0*5**2)*(ASH/PI)**3
414
 
      QCDH(X)=1.D0+HQCD(X)
415
 
      QCDA(X)=1.D0+AQCD(X)
416
 
      TQCDH(X)=1.D0+4.D0/3*HQCDM(BETA(X))*ASH/PI
417
 
      TQCDA(X)=1.D0+4.D0/3*AQCDM(BETA(X))*ASH/PI
418
 
      QCDC(X,Y)=1.D0+4/3.D0*ASH/PI*(9/4.D0 + (3-2*X+2*Y)/4*DLOG(X/Y)
419
 
     .         +((1.5D0-X-Y)*LAMB(X,Y)**2+5*X*Y)/2/LAMB(X,Y)
420
 
     .         /(1-X-Y)*DLOG(XI(X,Y)*XI(Y,X))
421
 
     .         + BIJ(X,Y))
422
 
     .         + ASH/PI*(2*(4/3.D0-DLOG(X))
423
 
     .         - (X*2*(4/3.D0-DLOG(X)) + Y*2*(4/3.D0-DLOG(Y)))/(1-X-Y)
424
 
     .         - (X*2*(4/3.D0-DLOG(X))*(1-X+Y)
425
 
     .           +Y*2*(4/3.D0-DLOG(Y))*(1+X-Y))/LAMB(X,Y)**2)
426
 
      QCDCI(X,Y)=1.D0+4/3.D0*ASH/PI*(3 + (Y-X)/2*DLOG(X/Y)
427
 
     .         +(2*(1-X-Y)+LAMB(X,Y)**2)/2/LAMB(X,Y)
428
 
     .         *DLOG(XI(X,Y)*XI(Y,X))
429
 
     .         + BIJ(X,Y))
430
 
     .         + ASH/PI*(2*(4/3.D0-DLOG(X)) + 2*(4/3.D0-DLOG(Y))
431
 
     .         - (X*2*(4/3.D0-DLOG(X))*(1-X+Y)
432
 
     .           +Y*2*(4/3.D0-DLOG(Y))*(1+X-Y))/LAMB(X,Y)**2)
433
 
      QCDCM(X,Y)=1.D0+4/3.D0*ASH/PI*(9/4.D0 + (3-2*X+2*Y)/4*DLOG(X/Y)
434
 
     .         +((1.5D0-X-Y)*LAMB(X,Y)**2+5*X*Y)/2/LAMB(X,Y)
435
 
     .         /(1-X-Y)*DLOG(4*X*Y/(1-X-Y+LAMB(X,Y))**2)
436
 
     .         + BIJ(X,Y))
437
 
      QCDCMI(X,Y)=1.D0+4/3.D0*ASH/PI*(3 + (Y-X)/2*DLOG(X/Y)
438
 
     .         +(2*(1-X-Y)*LAMB(X,Y)**2)/2/LAMB(X,Y)
439
 
     .         *DLOG(4*X*Y/(1-X-Y+LAMB(X,Y))**2)
440
 
     .         + BIJ(X,Y))
441
 
      CQCD(Z,TB,X,Y)= GF/(4*PI*DSQRT(2.D0))*Z**3*LAMB(X,Y)
442
 
     .              *((1.D0-X-Y)*(X*TB**2*QCDC(X,Y)
443
 
     .                           +Y/TB**2*QCDC(Y,X))
444
 
     .               -4.D0*X*Y*QCDCI(X,Y))
445
 
      CQCDM(Z,TB,X,Y)= GF/(4*PI*DSQRT(2.D0))*Z**3*LAMB(X,Y)
446
 
     .              *((1.D0-X-Y)*(X*TB**2*QCDCM(X,Y)
447
 
     .                           +Y/TB**2*QCDCM(Y,X))
448
 
     .               -4.D0*X*Y*QCDCMI(X,Y))
449
 
      ELW(AMH,AMF,QF,ACF)=ALPH/PI*3.D0/2*QF**2
450
 
     .                              *(3.D0/2-DLOG(AMH**2/AMF**2))
451
 
     .      +GF/8/DSQRT(2.D0)/PI**2*(ACF*AMT**2
452
 
     .        +AMW**2*(3*DLOG(CS)/SS-5)+AMZ**2*(0.5D0
453
 
     .          -3*(1-4*SS*DABS(QF))**2))
454
 
      CF(CA) = -CDLOG(-(1+CDSQRT(1-CA))/(1-CDSQRT(1-CA)))**2/4
455
 
      CG(CA) = CDSQRT(1-CA)/2*CDLOG(-(1+CDSQRT(1-CA))/(1-CDSQRT(1-CA)))
456
 
      CI1(CA,CB) = CA*CB/2/(CA-CB)
457
 
     .           + CA**2*CB**2/2/(CA-CB)**2*(CF(CA)-CF(CB))
458
 
     .           + CA**2*CB/(CA-CB)**2*(CG(CA)-CG(CB))
459
 
      CI2(CA,CB) = -CA*CB/2/(CA-CB)*(CF(CA)-CF(CB))
460
 
      HGGQCD(ASG,NF)=1.D0+ASG/PI*(95.D0/4.D0-NF*7.D0/6.D0)
461
 
      SGGQCD(ASG)=ASG/PI*17.D0/6.D0
462
 
      AGGQCD(ASG,NF)=1.D0+ASG/PI*(97.D0/4.D0-NF*7.D0/6.D0)
463
 
      HFFSELF(AMH)=1.D0+GF*AMH**2/16.D0/PI**2/DSQRT(2.D0)*2.117203D0
464
 
     .            -(GF*AMH**2/16.D0/PI**2/DSQRT(2.D0))**2*32.6567D0
465
 
      HVVSELF(AMH)=1.D0+GF*AMH**2/16.D0/PI**2/DSQRT(2.D0)*2.800952D0
466
 
     .            +(GF*AMH**2/16.D0/PI**2/DSQRT(2.D0))**2*62.0308D0
467
 
 
468
 
      PI=4D0*DATAN(1D0)
469
 
      SS=1.D0-(AMW/AMZ)**2
470
 
      CS=1.D0-SS
471
 
 
472
 
C--DECOUPLING THE TOP QUARK FROM ALPHAS
473
 
      AMT0=3.D8
474
 
 
475
 
C--TOP QUARK DECAY WIDTH
476
 
      GAMT0 = GF*AMT**3/8/DSQRT(2D0)/PI*(1-AMW**2/AMT**2)**2
477
 
     .                                 *(1+2*AMW**2/AMT**2)
478
 
      IF(IHIGGS.NE.0.AND.AMT.GT.AMCH+AMB)THEN
479
 
       GAMT1 = GF*AMT**3/8/DSQRT(2D0)/PI*(1-AMCH**2/AMT**2)**2
480
 
     .        *((AMB/AMT)**2*TGBET**2 + 1/TGBET**2)
481
 
      ELSE
482
 
       GAMT1 = 0
483
 
      ENDIF
484
 
      GAMT1 = GAMT0+GAMT1
485
 
 
486
 
      IF(IHIGGS.EQ.0)THEN
487
 
 
488
 
C        =========================================================
489
 
C                              SM HIGGS DECAYS
490
 
C        =========================================================
491
 
      AMXX=AMH
492
 
      AMH=AMSM
493
 
C     =============  RUNNING MASSES 
494
 
      RMS = RUNM(AMH,3)
495
 
      RMC = RUNM(AMH,4)
496
 
      RMB = RUNM(AMH,5)
497
 
      RMT = RUNM(AMH,6)
498
 
      RATCOUP = 1
499
 
      HIGTOP = AMH**2/AMT**2
500
 
 
501
 
      ASH=ALPHAS_HD(AMH,2)
502
 
      AMC0=1.D8
503
 
      AMB0=2.D8
504
 
      AS3=ALPHAS_HD(AMH,2)
505
 
      AMC0=AMC
506
 
      AS4=ALPHAS_HD(AMH,2)
507
 
      AMB0=AMB
508
 
C     AMT0=AMT
509
 
C     =============== PARTIAL WIDTHS 
510
 
C  H ---> G G
511
 
C
512
 
       EPS=1.D-8
513
 
       NFEXT = 3
514
 
       ASG = AS3
515
 
       CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
516
 
       CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
517
 
       CAT = 2*CTT*(1+(1-CTT)*CF(CTT))
518
 
       CAB = 2*CTB*(1+(1-CTB)*CF(CTB))
519
 
       FQCD=HGGQCD(ASG,NFEXT)
520
 
       XFAC = CDABS(CAT+CAB)**2*FQCD
521
 
       HGG=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8
522
 
 
523
 
C  H ---> G G* ---> G CC   TO BE ADDED TO H ---> CC
524
 
       NFEXT = 4
525
 
       ASG = AS4
526
 
       FQCD=HGGQCD(ASG,NFEXT)
527
 
       XFAC = CDABS(CAT+CAB)**2*FQCD
528
 
       DCC=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8 - HGG
529
 
C  H ---> G G* ---> G BB   TO BE ADDED TO H ---> BB
530
 
       NFEXT = 5
531
 
       ASG = ASH
532
 
       FQCD=HGGQCD(ASG,NFEXT)
533
 
       XFAC = CDABS(CAT+CAB)**2*FQCD
534
 
       DBB=HVV(AMH,0.D0)*(ASG/PI)**2*XFAC/8 - HGG - DCC
535
 
 
536
 
      IF(NFGG.EQ.5)THEN
537
 
       HGG = HGG + DBB + DCC
538
 
       DBB = 0
539
 
       DCC = 0
540
 
      ELSEIF(NFGG.EQ.4)THEN
541
 
       HGG = HGG + DCC
542
 
       DCC = 0
543
 
      ENDIF
544
 
 
545
 
C  H ---> MU MU
546
 
      IF(AMH.LE.2*AMMUON) THEN
547
 
       HMM = 0
548
 
      ELSE
549
 
      HMM=HFF(AMH,(AMMUON/AMH)**2)
550
 
     .    *(1+ELW(AMH,AMMUON,-1.D0,7.D0))
551
 
     .    *HFFSELF(AMH)
552
 
      ENDIF
553
 
C  H ---> TAU TAU
554
 
      IF(AMH.LE.2*AMTAU) THEN
555
 
       HLL = 0
556
 
      ELSE
557
 
      HLL=HFF(AMH,(AMTAU/AMH)**2)
558
 
     .    *(1+ELW(AMH,AMTAU,-1.D0,7.D0))
559
 
     .    *HFFSELF(AMH)
560
 
      ENDIF
561
 
C  H --> SS
562
 
      IF(AMH.LE.2*AMS) THEN
563
 
       HSS = 0
564
 
      ELSE
565
 
       HS2=3.D0*HFF(AMH,(RMS/AMH)**2)
566
 
     .    *QCDH(RMS**2/AMH**2)
567
 
     .    *(1+ELW(AMH,RMS,-1.D0/3.D0,7.D0))
568
 
     .    *HFFSELF(AMH)
569
 
       IF(HS2.LT.0.D0) HS2 = 0
570
 
       HS1=3.D0*HFF(AMH,(AMS/AMH)**2)
571
 
     .    *TQCDH(AMS**2/AMH**2)
572
 
     .    *HFFSELF(AMH)
573
 
       RAT = 2*AMS/AMH
574
 
       HSS = QQINT(RAT,HS1,HS2)
575
 
      ENDIF
576
 
C  H --> CC
577
 
      IF(AMH.LE.2*AMC) THEN
578
 
       HCC = 0
579
 
      ELSE
580
 
       HC2=3.D0*HFF(AMH,(RMC/AMH)**2)
581
 
     .    *QCDH(RMC**2/AMH**2)
582
 
     .    *(1+ELW(AMH,RMC,2.D0/3.D0,7.D0))
583
 
     .    *HFFSELF(AMH)
584
 
     .   + DCC
585
 
       IF(HC2.LT.0.D0) HC2 = 0
586
 
       HC1=3.D0*HFF(AMH,(AMC/AMH)**2)
587
 
     .    *TQCDH(AMC**2/AMH**2)
588
 
     .    *HFFSELF(AMH)
589
 
       RAT = 2*AMC/AMH
590
 
       HCC = QQINT(RAT,HC1,HC2)
591
 
      ENDIF
592
 
C  H --> BB :
593
 
      IF(AMH.LE.2*AMB) THEN
594
 
       HBB = 0
595
 
      ELSE
596
 
       HB2=3.D0*HFF(AMH,(RMB/AMH)**2)
597
 
     .    *QCDH(RMB**2/AMH**2)
598
 
     .    *(1+ELW(AMH,RMB,-1.D0/3.D0,1.D0))
599
 
     .    *HFFSELF(AMH)
600
 
     .   + DBB
601
 
       IF(HB2.LT.0.D0) HB2 = 0
602
 
       HB1=3.D0*HFF(AMH,(AMB/AMH)**2)
603
 
     .    *TQCDH(AMB**2/AMH**2)
604
 
     .    *HFFSELF(AMH)
605
 
       RAT = 2*AMB/AMH
606
 
       HBB = QQINT(RAT,HB1,HB2)
607
 
      ENDIF
608
 
C  H ---> TT
609
 
      RATCOUP = 0
610
 
      IF(IONSH.EQ.0)THEN
611
 
       DLD=3D0
612
 
       DLU=5D0
613
 
       XM1 = 2D0*AMT-DLD
614
 
       XM2 = 2D0*AMT+DLU
615
 
       IF (AMH.LE.AMT+AMW+AMB) THEN
616
 
       HTT=0.D0
617
 
       ELSEIF (AMH.LE.XM1) THEN
618
 
        FACTT=6.D0*GF**2*AMH**3*AMT**2/2.D0/128.D0/PI**3
619
 
        CALL HTOTTS(AMH,AMT,AMB,AMW,HTTS)
620
 
        HTT=FACTT*HTTS
621
 
       ELSEIF (AMH.LE.XM2) THEN
622
 
        XX(1) = XM1-1D0
623
 
        XX(2) = XM1
624
 
        XX(3) = XM2
625
 
        XX(4) = XM2+1D0
626
 
        FACTT=6.D0*GF**2*XX(1)**3*AMT**2/2.D0/128.D0/PI**3
627
 
        CALL HTOTTS(XX(1),AMT,AMB,AMW,HTTS)
628
 
        YY(1)=FACTT*HTTS
629
 
        FACTT=6.D0*GF**2*XX(2)**3*AMT**2/2.D0/128.D0/PI**3
630
 
        CALL HTOTTS(XX(2),AMT,AMB,AMW,HTTS)
631
 
        YY(2)=FACTT*HTTS
632
 
        XMT = RUNM(XX(3),6)
633
 
        XY2=3.D0*HFF(XX(3),(XMT/XX(3))**2)
634
 
     .    *QCDH(XMT**2/XX(3)**2)
635
 
     .    *HFFSELF(XX(3))
636
 
        IF(XY2.LT.0.D0) XY2 = 0
637
 
        XY1=3.D0*HFF(XX(3),(AMT/XX(3))**2)
638
 
     .    *TQCDH(AMT**2/XX(3)**2)
639
 
     .    *HFFSELF(XX(3))
640
 
        RAT = 2*AMT/XX(3)
641
 
        YY(3) = QQINT(RAT,XY1,XY2)
642
 
        XMT = RUNM(XX(4),6)
643
 
        XY2=3.D0*HFF(XX(4),(XMT/XX(4))**2)
644
 
     .    *QCDH(XMT**2/XX(4)**2)
645
 
     .    *HFFSELF(XX(4))
646
 
        IF(XY2.LT.0.D0) XY2 = 0
647
 
        XY1=3.D0*HFF(XX(4),(AMT/XX(4))**2)
648
 
     .    *TQCDH(AMT**2/XX(4)**2)
649
 
     .    *HFFSELF(XX(4))
650
 
        RAT = 2*AMT/XX(4)
651
 
        YY(4) = QQINT(RAT,XY1,XY2)
652
 
        HTT = FINT(AMH,XX,YY)
653
 
       ELSE
654
 
        HT2=3.D0*HFF(AMH,(RMT/AMH)**2)
655
 
     .    *QCDH(RMT**2/AMH**2)
656
 
     .    *HFFSELF(AMH)
657
 
        IF(HT2.LT.0.D0) HT2 = 0
658
 
        HT1=3.D0*HFF(AMH,(AMT/AMH)**2)
659
 
     .    *TQCDH(AMT**2/AMH**2)
660
 
     .    *HFFSELF(AMH)
661
 
        RAT = 2*AMT/AMH
662
 
        HTT = QQINT(RAT,HT1,HT2)
663
 
       ENDIF
664
 
      ELSE
665
 
       IF (AMH.LE.2.D0*AMT) THEN
666
 
        HTT=0.D0
667
 
       ELSE
668
 
        HT2=3.D0*HFF(AMH,(RMT/AMH)**2)
669
 
     .    *QCDH(RMT**2/AMH**2)
670
 
     .    *HFFSELF(AMH)
671
 
        IF(HT2.LT.0.D0) HT2 = 0
672
 
        HT1=3.D0*HFF(AMH,(AMT/AMH)**2)
673
 
     .    *TQCDH(AMT**2/AMH**2)
674
 
     .    *HFFSELF(AMH)
675
 
        RAT = 2*AMT/AMH
676
 
        HTT = QQINT(RAT,HT1,HT2)
677
 
       ENDIF
678
 
      ENDIF
679
 
C  H ---> GAMMA GAMMA
680
 
       EPS=1.D-8
681
 
       CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
682
 
       CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
683
 
       CTC = 4*AMC**2/AMH**2*DCMPLX(1D0,-EPS)
684
 
       CTL = 4*AMTAU**2/AMH**2*DCMPLX(1D0,-EPS)
685
 
       CTW = 4*AMW**2/AMH**2*DCMPLX(1D0,-EPS)
686
 
       CAW = -(2+3*CTW+3*CTW*(2-CTW)*CF(CTW))
687
 
       CAT = 4/3D0 * 2*CTT*(1+(1-CTT)*CF(CTT))
688
 
       CAB = 1/3D0 * 2*CTB*(1+(1-CTB)*CF(CTB))
689
 
       CAC = 4/3D0 * 2*CTC*(1+(1-CTC)*CF(CTC))
690
 
       CAL =         2*CTL*(1+(1-CTL)*CF(CTL))
691
 
       XFAC = CDABS(CAT+CAB+CAC+CAL+CAW)**2
692
 
       HGA=HVV(AMH,0.D0)*(ALPH/PI)**2/16.D0*XFAC
693
 
C  H ---> Z GAMMA
694
 
      IF(AMH.LE.AMZ)THEN
695
 
       HZGA=0
696
 
      ELSE
697
 
       EPS=1.D-8
698
 
       TS = SS/CS
699
 
       FT = -3*2D0/3*(1-4*2D0/3*SS)/DSQRT(SS*CS)
700
 
       FB = 3*1D0/3*(-1+4*1D0/3*SS)/DSQRT(SS*CS)
701
 
       CTT = 4*AMT**2/AMH**2*DCMPLX(1D0,-EPS)
702
 
       CTB = 4*AMB**2/AMH**2*DCMPLX(1D0,-EPS)
703
 
       CTW = 4*AMW**2/AMH**2*DCMPLX(1D0,-EPS)
704
 
       CLT = 4*AMT**2/AMZ**2*DCMPLX(1D0,-EPS)
705
 
       CLB = 4*AMB**2/AMZ**2*DCMPLX(1D0,-EPS)
706
 
       CLW = 4*AMW**2/AMZ**2*DCMPLX(1D0,-EPS)
707
 
       CAT = FT*(CI1(CTT,CLT) - CI2(CTT,CLT))
708
 
       CAB = FB*(CI1(CTB,CLB) - CI2(CTB,CLB))
709
 
       CAW = -1/DSQRT(TS)*(4*(3-TS)*CI2(CTW,CLW)
710
 
     .     + ((1+2/CTW)*TS - (5+2/CTW))*CI1(CTW,CLW))
711
 
       XFAC = CDABS(CAT+CAB+CAW)**2
712
 
       ACOUP = DSQRT(2D0)*GF*AMZ**2*SS*CS/PI**2
713
 
       HZGA = GF/(4.D0*PI*DSQRT(2.D0))*AMH**3*(ALPH/PI)*ACOUP/16.D0
714
 
     .        *XFAC*(1-AMZ**2/AMH**2)**3
715
 
      ENDIF
716
 
C  H ---> W W
717
 
      IF(IONWZ.EQ.0)THEN
718
 
       DLD=2D0
719
 
       DLU=2D0
720
 
       XM1 = 2D0*AMW-DLD
721
 
       XM2 = 2D0*AMW+DLU
722
 
       IF (AMH.LE.XM1) THEN
723
 
        CALL HTOVV(AMH,AMW,GAMW,HTWW)
724
 
        HWW = 3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/AMH**3*HTWW
725
 
       ELSEIF (AMH.LE.XM2) THEN
726
 
        XX(1) = XM1-1D0
727
 
        XX(2) = XM1
728
 
        XX(3) = XM2
729
 
        XX(4) = XM2+1D0
730
 
        CALL HTOVV(XX(1),AMW,GAMW,HTWW)
731
 
        YY(1)=3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/XX(1)**3*HTWW
732
 
        CALL HTOVV(XX(2),AMW,GAMW,HTWW)
733
 
        YY(2)=3D0/2D0*GF*AMW**4/DSQRT(2D0)/PI/XX(2)**3*HTWW
734
 
        YY(3)=HVV(XX(3),AMW**2/XX(3)**2)
735
 
     .       *HVVSELF(XX(3))
736
 
        YY(4)=HVV(XX(4),AMW**2/XX(4)**2)
737
 
     .       *HVVSELF(XX(4))
738
 
        HWW = FINT(AMH,XX,YY)
739
 
       ELSE
740
 
        HWW=HVV(AMH,AMW**2/AMH**2)
741
 
     .     *HVVSELF(AMH)
742
 
       ENDIF
743
 
      ELSE
744
 
       DLD=2D0
745
 
       DLU=2D0
746
 
       XM1 = 2D0*AMW-DLD
747
 
       XM2 = 2D0*AMW+DLU
748
 
      IF (AMH.LE.AMW) THEN
749
 
       HWW=0
750
 
      ELSE IF (AMH.LE.XM1) THEN
751
 
       CWW=3.D0*GF**2*AMW**4/16.D0/PI**3
752
 
       HWW=HV(AMW**2/AMH**2)*CWW*AMH
753
 
      ELSE IF (AMH.LT.XM2) THEN
754
 
       CWW=3.D0*GF**2*AMW**4/16.D0/PI**3
755
 
       XX(1) = XM1-1D0
756
 
       XX(2) = XM1
757
 
       XX(3) = XM2
758
 
       XX(4) = XM2+1D0
759
 
       YY(1)=HV(AMW**2/XX(1)**2)*CWW*XX(1)
760
 
       YY(2)=HV(AMW**2/XX(2)**2)*CWW*XX(2)
761
 
       YY(3)=HVV(XX(3),AMW**2/XX(3)**2)
762
 
     .      *HVVSELF(XX(3))
763
 
       YY(4)=HVV(XX(4),AMW**2/XX(4)**2)
764
 
     .      *HVVSELF(XX(4))
765
 
       HWW = FINT(AMH,XX,YY)
766
 
      ELSE
767
 
       HWW=HVV(AMH,AMW**2/AMH**2)
768
 
     .     *HVVSELF(AMH)
769
 
      ENDIF
770
 
      ENDIF
771
 
C  H ---> Z Z
772
 
      IF(IONWZ.EQ.0)THEN
773
 
       DLD=2D0
774
 
       DLU=2D0
775
 
       XM1 = 2D0*AMZ-DLD
776
 
       XM2 = 2D0*AMZ+DLU
777
 
       IF (AMH.LE.XM1) THEN
778
 
        CALL HTOVV(AMH,AMZ,GAMZ,HTZZ)
779
 
        HZZ = 3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/AMH**3*HTZZ
780
 
       ELSEIF (AMH.LE.XM2) THEN
781
 
        XX(1) = XM1-1D0
782
 
        XX(2) = XM1
783
 
        XX(3) = XM2
784
 
        XX(4) = XM2+1D0
785
 
        CALL HTOVV(XX(1),AMZ,GAMZ,HTZZ)
786
 
        YY(1)=3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/XX(1)**3*HTZZ
787
 
        CALL HTOVV(XX(2),AMZ,GAMZ,HTZZ)
788
 
        YY(2)=3D0/4D0*GF*AMZ**4/DSQRT(2D0)/PI/XX(2)**3*HTZZ
789
 
        YY(3)=HVV(XX(3),AMZ**2/XX(3)**2)/2
790
 
     .       *HVVSELF(XX(3))
791
 
        YY(4)=HVV(XX(4),AMZ**2/XX(4)**2)/2
792
 
     .       *HVVSELF(XX(4))
793
 
        HZZ = FINT(AMH,XX,YY)
794
 
       ELSE
795
 
        HZZ=HVV(AMH,AMZ**2/AMH**2)/2.D0
796
 
     .      *HVVSELF(AMH)
797
 
       ENDIF
798
 
      ELSE
799
 
      DLD=2D0
800
 
      DLU=2D0
801
 
      XM1 = 2D0*AMZ-DLD
802
 
      XM2 = 2D0*AMZ+DLU
803
 
      IF (AMH.LE.AMZ) THEN
804
 
      HZZ=0
805
 
      ELSE IF (AMH.LE.XM1) THEN
806
 
      CZZ=3.D0*GF**2*AMZ**4/192.D0/PI**3*(7-40/3.D0*SS+160/9.D0*SS**2)
807
 
      HZZ=HV(AMZ**2/AMH**2)*CZZ*AMH
808
 
      ELSE IF (AMH.LT.XM2) THEN
809
 
      CZZ=3.D0*GF**2*AMZ**4/192.D0/PI**3*(7-40/3.D0*SS+160/9.D0*SS**2)
810
 
      XX(1) = XM1-1D0
811
 
      XX(2) = XM1
812
 
      XX(3) = XM2
813
 
      XX(4) = XM2+1D0
814
 
      YY(1)=HV(AMZ**2/XX(1)**2)*CZZ*XX(1)
815
 
      YY(2)=HV(AMZ**2/XX(2)**2)*CZZ*XX(2)
816
 
      YY(3)=HVV(XX(3),AMZ**2/XX(3)**2)/2
817
 
     .      *HVVSELF(XX(3))
818
 
      YY(4)=HVV(XX(4),AMZ**2/XX(4)**2)/2
819
 
     .      *HVVSELF(XX(4))
820
 
      HZZ = FINT(AMH,XX,YY)
821
 
      ELSE
822
 
      HZZ=HVV(AMH,AMZ**2/AMH**2)/2.D0
823
 
     .   *HVVSELF(AMH)
824
 
      ENDIF
825
 
      ENDIF
826
 
C
827
 
C    ==========  TOTAL WIDTH AND BRANCHING RATIOS 
828
 
C
829
 
      WTOT=HLL+HMM+HSS+HCC+HBB+HTT+HGG+HGA+HZGA+HWW+HZZ
830
 
      SMBRT=HTT/WTOT
831
 
      SMBRB=HBB/WTOT
832
 
      SMBRL=HLL/WTOT
833
 
      SMBRM=HMM/WTOT
834
 
      SMBRC=HCC/WTOT
835
 
      SMBRS=HSS/WTOT
836
 
      SMBRG=HGG/WTOT
837
 
      SMBRGA=HGA/WTOT
838
 
      SMBRZGA=HZGA/WTOT
839
 
      SMBRW=HWW/WTOT
840
 
      SMBRZ=HZZ/WTOT
841
 
      SMWDTH=WTOT
842
 
 
843
 
      AMH=AMXX
844
 
 
845
 
      endif
846
 
 
847
 
      RETURN
848
 
      END
849
 
 
850
 
      DOUBLE PRECISION FUNCTION BIJ(X,Y)
851
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
852
 
      DOUBLE PRECISION LAMB
853
 
      BIJ = (1-X-Y)/LAMB(X,Y)*(
854
 
     .          4*SP(XI(X,Y)*XI(Y,X))
855
 
     .        - 2*SP(-XI(X,Y)) - 2*SP(-XI(Y,X))
856
 
     .        + 2*DLOG(XI(X,Y)*XI(Y,X))*DLOG(1-XI(X,Y)*XI(Y,X))
857
 
     .        - DLOG(XI(X,Y))*DLOG(1+XI(X,Y))
858
 
     .        - DLOG(XI(Y,X))*DLOG(1+XI(Y,X))
859
 
     .          )
860
 
     .        -4*(DLOG(1-XI(X,Y)*XI(Y,X))
861
 
     .      +XI(X,Y)*XI(Y,X)/(1-XI(X,Y)*XI(Y,X))*DLOG(XI(X,Y)*XI(Y,X)))
862
 
     .        + (LAMB(X,Y)+X-Y)/LAMB(X,Y)*(DLOG(1+XI(X,Y))
863
 
     .                 - XI(X,Y)/(1+XI(X,Y))*DLOG(XI(X,Y)))
864
 
     .        + (LAMB(X,Y)-X+Y)/LAMB(X,Y)*(DLOG(1+XI(Y,X))
865
 
     .                 - XI(Y,X)/(1+XI(Y,X))*DLOG(XI(Y,X)))
866
 
      RETURN
867
 
      END
868
 
 
869
 
      DOUBLE PRECISION FUNCTION BETA(X)
870
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
871
 
      BETA=DSQRT(1.D0-4.D0*X)
872
 
      RETURN
873
 
      END
874
 
 
875
 
      DOUBLE PRECISION FUNCTION LAMB(X,Y)
876
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
877
 
      LAMB=DSQRT((1.D0-X-Y)**2-4.D0*X*Y)
878
 
      RETURN
879
 
      END
880
 
 
881
 
      DOUBLE PRECISION FUNCTION XI(X,Y)
882
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
883
 
      DOUBLE PRECISION LAMB
884
 
      XI = 2*X/(1-X-Y+LAMB(X,Y))
885
 
      RETURN
886
 
      END
887
 
 
888
 
 
889
 
C ===================== THE FUNCTION F0 ===============
890
 
      COMPLEX*16 FUNCTION F0(M1,M2,QSQ)
891
 
      IMPLICIT REAL*8 (A-H,M,O-Z)
892
 
      COMPLEX*16 CD,CR,CQ2,IEPS,CBET,CXX
893
 
      M1SQ = M1*M1
894
 
      M2SQ = M2*M2
895
 
      AQSQ = DABS(QSQ)
896
 
      IEPS = DCMPLX(1.D0,1.D-12)
897
 
      CQ2 = QSQ*IEPS
898
 
      CD = (M1SQ-M2SQ)/CQ2
899
 
      CR = CDSQRT((1+CD)**2 - 4*M1SQ/CQ2)
900
 
      IF(QSQ.EQ.0.D0) THEN
901
 
       F0 = 0.D0
902
 
      ELSE
903
 
       IF(M1.EQ.M2) THEN
904
 
        F0 = -2.D0 + CR*CDLOG(-(1+CR)/(1-CR))
905
 
       ELSE
906
 
        CBET = CDSQRT(1-4*M1*M2/(CQ2 - (M1-M2)**2))
907
 
        CXX = (CBET-1)/(CBET+1)
908
 
        F0 = -1 + ((QSQ+M2SQ-M1SQ)/2/QSQ - M2SQ/(M2SQ-M1SQ))
909
 
     .                                           *DLOG(M2SQ/M1SQ)
910
 
     .     - (QSQ-(M1-M2)**2)/QSQ*CBET*CDLOG(CXX)
911
 
       ENDIF
912
 
      ENDIF
913
 
      RETURN
914
 
      END
915
 
 
916
 
C     ************************************************************
917
 
C     SUBROUTINE FOR HSM ---> V*V* ---> 4F
918
 
C     ************************************************************
919
 
      SUBROUTINE HTOVV(AMH,AMV,GAMV,HTVV)
920
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
921
 
      COMMON/VVOFF/AMH1,AMV1,GAMV1
922
 
      COMMON/PREC/IP
923
 
      EXTERNAL FTOVV1
924
 
      IP=20
925
 
      AMH1=AMH
926
 
      AMV1=AMV
927
 
      GAMV1=GAMV
928
 
      DLT=1D0/IP
929
 
      SUM=0D0
930
 
      DO 1 I=1,IP
931
 
       UU=DLT*I
932
 
       DD=UU-DLT
933
 
       CALL QGAUS1(FTOVV1,DD,UU,RES)
934
 
       SUM=SUM+RES
935
 
1     CONTINUE
936
 
      HTVV=SUM
937
 
      RETURN
938
 
      END
939
 
 
940
 
      DOUBLE PRECISION FUNCTION FTOVV1(XX)
941
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
942
 
      COMMON/FIRST/X1
943
 
      COMMON/PREC/IP
944
 
      EXTERNAL FTOVV2
945
 
      X1=XX
946
 
      DLT=1D0/IP
947
 
      SUM=0D0
948
 
      DO 1 I=1,IP
949
 
       UU=DLT*I
950
 
       DD=UU-DLT
951
 
       CALL QGAUS2(FTOVV2,DD,UU,RES)
952
 
       SUM=SUM+RES
953
 
1     CONTINUE
954
 
      FTOVV1=SUM
955
 
      RETURN
956
 
      END
957
 
 
958
 
      DOUBLE PRECISION FUNCTION FTOVV2(XX)
959
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
960
 
      DIMENSION YY(2)
961
 
      COMMON/FIRST/X1
962
 
      YY(1)=X1
963
 
      YY(2)=XX
964
 
      FTOVV2=FTOVV(YY)
965
 
      RETURN
966
 
      END
967
 
 
968
 
      DOUBLE PRECISION FUNCTION FTOVV(XX)
969
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
970
 
      DOUBLE PRECISION LAMB
971
 
      DIMENSION XX(2)
972
 
      COMMON/VVOFF/AMH,AMW,GAMW
973
 
      LAMB(X,Y)=DSQRT((1.D0-X-Y)**2-4.D0*X*Y)
974
 
      PI=4D0*DATAN(1D0)
975
 
      ICASE = 1
976
 
      IF(ICASE.EQ.0)THEN
977
 
       YY = AMH**2
978
 
       Y1 = DATAN((YY-AMW**2)/AMW/GAMW)
979
 
       Y2 = -DATAN((AMW**2)/AMW/GAMW)
980
 
       DJAC = Y1-Y2
981
 
       T1 = TAN(Y1*XX(1)+Y2*(1.D0-XX(1)))
982
 
       SP = AMW**2 + AMW*GAMW*T1
983
 
       YY = (AMH-DSQRT(SP))**2
984
 
       Y1 = DATAN((YY-AMW**2)/AMW/GAMW)
985
 
       Y2 = -DATAN((AMW**2)/AMW/GAMW)
986
 
       DJAC = DJAC*(Y1-Y2)
987
 
       T2 = TAN(Y1*XX(2)+Y2*(1.D0-XX(2)))
988
 
       SM = AMW**2 + AMW*GAMW*T2
989
 
       AM2=AMH**2
990
 
       GAM = AM2*LAMB(SP/AM2,SM/AM2)*(1+LAMB(SP/AM2,SM/AM2)**2*AMH**4
991
 
     .                               /SP/SM/12)
992
 
       PRO1 = SP/AMW**2
993
 
       PRO2 = SM/AMW**2
994
 
       FTOVV = PRO1*PRO2*GAM*DJAC/PI**2
995
 
      ELSE
996
 
       SP = AMH**2*XX(1)
997
 
       SM = (AMH-DSQRT(SP))**2*XX(2)
998
 
       DJAC = AMH**2*(AMH-DSQRT(SP))**2/PI**2
999
 
       AM2=AMH**2
1000
 
       GAM = AM2*LAMB(SP/AM2,SM/AM2)*(1+LAMB(SP/AM2,SM/AM2)**2*AMH**4
1001
 
     .                               /SP/SM/12)
1002
 
       PRO1 = SP*GAMW/AMW/((SP-AMW**2)**2+AMW**2*GAMW**2)
1003
 
       PRO2 = SM*GAMW/AMW/((SM-AMW**2)**2+AMW**2*GAMW**2)
1004
 
       FTOVV = PRO1*PRO2*GAM*DJAC
1005
 
      ENDIF
1006
 
      RETURN
1007
 
      END
1008
 
 
1009
 
C     ************************************************************
1010
 
C     SUBROUTINE FOR HSM ---> TT* ---> TBW
1011
 
C     ************************************************************
1012
 
      SUBROUTINE HTOTTS(AMH,AMT,AMB,AMW,HTTS)
1013
 
      IMPLICIT REAL*8(A-Z)
1014
 
      INTEGER IP,K
1015
 
      COMMON/PREC1/IP
1016
 
      EXTERNAL FUNSTT1
1017
 
      COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
1018
 
      COMMON/TOP0/AMH0,AMT0,AMB0,AMW0
1019
 
      AMH0=AMH
1020
 
      AMT0=AMT
1021
 
      AMB0=AMB
1022
 
      AMW0=AMW
1023
 
      IP=5
1024
 
      M1=AMB
1025
 
      M2=AMT
1026
 
      M3=AMW
1027
 
C     FIRST INTEGRATE OVER X2, i.e. (1+3) SYSTEM
1028
 
C        CHECK WHETHER ENOUGH PHASE SPACE
1029
 
      MASTOT=M1+M2+M3
1030
 
      IF(MASTOT.GE.AMH) GOTO 12
1031
 
      ECM=AMH
1032
 
      S=ECM**2
1033
 
      U1=(ECM-M2)**2
1034
 
      D1=(M1+M3)**2
1035
 
      U=(S-D1+M2**2)/s
1036
 
      D=(S-U1+M2**2)/s
1037
 
      DEL=(U-D)/IP
1038
 
      U=D+DEL
1039
 
      XSEC=0.D0
1040
 
      DO K=1,IP
1041
 
      CALL QGAUS1(FUNSTT1,D,U,SS)
1042
 
      D=U
1043
 
      U=D+DEL
1044
 
      XSEC=XSEC+SS
1045
 
      ENDDO
1046
 
      HTTS=XSEC
1047
 
12    CONTINUE
1048
 
      RETURN
1049
 
      END
1050
 
 
1051
 
      DOUBLE PRECISION FUNCTION FUNSTT1(XL)
1052
 
      IMPLICIT REAL*8(A-Z)
1053
 
      INTEGER IP,I
1054
 
      COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
1055
 
      COMMON/PREC1/IP
1056
 
      EXTERNAL FUNSTT2
1057
 
      X2=XL
1058
 
      S13=S-S*X2+M2**2
1059
 
      TEM=2.D0*DSQRT(S13)
1060
 
      E2S=(S-S13-M2**2)/TEM
1061
 
      E3S=(S13+M3**2-M1**2)/TEM
1062
 
C     SECOND INTEGRAL OVER X1, i.e. (2+3) SYSTEM
1063
 
      U1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)-DSQRT(E3S**2-M3**2))**2
1064
 
      D1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)+DSQRT(E3S**2-M3**2))**2
1065
 
      U=(S-D1+M1**2)/s
1066
 
      D=(S-U1+M1**2)/s
1067
 
      DEL=(U-D)/IP
1068
 
      FUNSTT1=0.d0
1069
 
      U=D+DEL
1070
 
      DO I=1,IP
1071
 
      CALL QGAUS2(FUNSTT2,D,U,SS)
1072
 
      FUNSTT1=FUNSTT1+SS
1073
 
      D=U
1074
 
      U=D+DEL
1075
 
      ENDDO
1076
 
      RETURN
1077
 
      END
1078
 
 
1079
 
      DOUBLE PRECISION FUNCTION FUNSTT2(XK)
1080
 
      IMPLICIT REAL*8(A-Z)
1081
 
      COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
1082
 
      X1=XK
1083
 
      CALL ELEMSTT(SS)
1084
 
      FUNSTT2=SS
1085
 
      RETURN
1086
 
      END
1087
 
 
1088
 
      SUBROUTINE ELEMSTT(RES)
1089
 
      IMPLICIT REAL*8(A-Z)
1090
 
      COMMON/IKSY0/X1,X2,M1,M2,M3,ECM,S
1091
 
      COMMON/TOP0/AMH,AMT,AMB,AMW
1092
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW0,GAMZ0
1093
 
      GAMT=GAMT0**2*AMT**2/AMH**4
1094
 
      GAMW=GAMW0**2*AMW**2/AMH**4
1095
 
      W=AMW**2/AMH**2
1096
 
      T=AMT**2/AMH**2
1097
 
      Y1=1-X2
1098
 
      Y2=1-X1
1099
 
      X0=2.D0-X1-X2
1100
 
      W1=(1.D0-X2)
1101
 
      W3=(1.-X1-X2)
1102
 
      W11=1.D0/((1.D0-X2)**2+GAMT)
1103
 
      W33=1.D0/(W3**2+GAMW**2)
1104
 
      W13=W1*W3*W11*W33
1105
 
 
1106
 
      R11=4*T*W-16.*T*W*Y1-4.*T*Y2*Y1+8.*T*Y1+32.*T*W**2-20
1107
 
     . .*T*Y1**2+8.*W*Y2*Y1+4.*W*Y1**2-4.*Y2*Y1**2-16.*T**2*W-
1108
 
     .  32.*T**2*Y1+4.*T**2-16.*T**3-8.*W**2+4.*Y1**2-4.*Y1**3
1109
 
      R33=-4.*T*W+4.*T*W*Y2-2.*T*W*Y2*Y1+4.*T*W*Y1+T*W*Y2**2-
1110
 
     .  3.*T*W*Y1**2+2.*T*Y2*Y1-3.*T*Y2*Y1**2+4.*T*W**2-4.*T*W**3
1111
 
     .  +T*Y2**2-3.*T*Y2**2*Y1-T*Y2**3+T*Y1**2-T*Y1**3+4.*T**2
1112
 
     .  *W-4.*T**2*W*Y2-4.*T**2*W*Y1-2.*T**2*Y2*Y1-4.*T**2*W**2-
1113
 
     .  T**2*Y2**2-T**2*Y1**2+4.*W**2*Y2*Y1-8.*W**3*Y2-8.*W**3*Y1
1114
 
     .  +4.*W**3+8.*W**4
1115
 
      R13=8.*W-24.*T*W+16.*T*W*Y1 -4.*T*Y2+16.*T*Y2*Y1-4.*T*
1116
 
     .  Y1+16.*T*W**2+4.*T*Y2**2+12.*T*Y1**2-8.*W*Y2-12.*W*Y2*Y1
1117
 
     .  -8.*W*Y1+4.*W*Y1**2-4.*Y2*Y1+8.*Y2*Y1**2+16.*T**2*W+8.
1118
 
     .  *T**2*Y2+8.*T**2*Y1+16.*W**2*Y2+24.*W**2*Y1+4.*Y2**2*Y1-
1119
 
     .  32.*W**3-4.*Y1**2+4.*Y1**3
1120
 
      RES=R11*W11+4.D0*R33*W33/T-2.D0*R13*W13
1121
 
      RETURN
1122
 
      END
1123
 
 
1124
 
C     **************************************************
1125
 
C     SUBROUTINE FOR A -> TT* -> TBW
1126
 
C     **************************************************
1127
 
 
1128
 
      SUBROUTINE ATOTT(AMA,AMT,AMB,AMW,AMCH,ATT0)
1129
 
      IMPLICIT REAL*8(A-Z)
1130
 
      INTEGER IP,K
1131
 
      COMMON/PREC1/IP
1132
 
      EXTERNAL FUNATT1
1133
 
      COMMON/IKSY1/X1,X2,M1,M2,M3,ECM,S
1134
 
      COMMON/TOP1/AMA1,AMT1,AMB1,AMW1,AMCH1
1135
 
      AMA1=AMA
1136
 
      AMT1=AMT
1137
 
      AMB1=AMB
1138
 
      AMW1=AMW
1139
 
      AMCH1=AMCH
1140
 
      IP=5
1141
 
      M1=AMB
1142
 
      M2=AMT
1143
 
      M3=AMW
1144
 
C        FIRST INTEGRATE OVER X2, i.e. (1+3) SYSTEM
1145
 
C        CHECK WHETHER ENOUGH PHASE SPACE
1146
 
      MASTOT=M1+M2+M3
1147
 
      IF(MASTOT.GE.AMA) GOTO 12
1148
 
      ECM=AMA
1149
 
      S=ECM**2
1150
 
      U1=(ECM-M2)**2
1151
 
      D1=(M1+M3)**2
1152
 
      U=(S-D1+M2**2)/s
1153
 
      D=(S-U1+M2**2)/s
1154
 
      DEL=(U-D)/IP
1155
 
      U=D+DEL
1156
 
      XSEC=0.D0
1157
 
      DO K=1,IP
1158
 
      CALL QGAUS1(FUNATT1,D,U,SS)
1159
 
      D=U
1160
 
      U=D+DEL
1161
 
      XSEC=XSEC+SS
1162
 
      ENDDO
1163
 
      ATT0=XSEC
1164
 
 12   CONTINUE
1165
 
      RETURN
1166
 
      END
1167
 
 
1168
 
      DOUBLE PRECISION FUNCTION FUNATT1(XL)
1169
 
      IMPLICIT REAL*8(A-Z)
1170
 
      INTEGER IP,I
1171
 
      COMMON/IKSY1/X1,X2,M1,M2,M3,ECM,S
1172
 
      COMMON/PREC1/IP
1173
 
      EXTERNAL FUNATT2
1174
 
      X2=XL
1175
 
      S13=S-S*X2+M2**2
1176
 
      TEM=2.D0*DSQRT(S13)
1177
 
      E2S=(S-S13-M2**2)/TEM
1178
 
      E3S=(S13+M3**2-M1**2)/TEM
1179
 
C     SECOND INTEGRAL OVER X1, i.e. (2+3) SYSTEM
1180
 
      U1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)-DSQRT(E3S**2-M3**2))**2
1181
 
      D1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)+DSQRT(E3S**2-M3**2))**2
1182
 
      U=(S-D1+M1**2)/s
1183
 
      D=(S-U1+M1**2)/s
1184
 
      DEL=(U-D)/IP
1185
 
      FUNATT1=0.d0
1186
 
      U=D+DEL
1187
 
      DO I=1,IP
1188
 
      CALL QGAUS2(FUNATT2,D,U,SS)
1189
 
      FUNATT1=FUNATT1+SS
1190
 
      D=U
1191
 
      U=D+DEL
1192
 
      ENDDO
1193
 
      RETURN
1194
 
      END
1195
 
 
1196
 
      DOUBLE PRECISION FUNCTION FUNATT2(XK)
1197
 
      IMPLICIT REAL*8(A-Z)
1198
 
      COMMON/IKSY1/X1,X2,M1,M2,M3,ECM,S
1199
 
      X1=XK
1200
 
      CALL ELEMATT(SS)
1201
 
      FUNATT2=SS
1202
 
      RETURN
1203
 
      END
1204
 
 
1205
 
      SUBROUTINE ELEMATT(RES)
1206
 
      IMPLICIT REAL*8(A-Z)
1207
 
      COMMON/IKSY1/X1,X2,M1,M2,M3,ECM,S
1208
 
      COMMON/TOP1/AMA,AMT,AMB,AMW,AMCH
1209
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW,GAMZ
1210
 
      GAMT=GAMT1**2*AMT**2/AMA**4
1211
 
      GAMC=GAMC0**2*AMCH**2/AMA**4
1212
 
      CH=AMCH**2/AMA**2
1213
 
      W=AMW**2/AMA**2
1214
 
      T=AMT**2/AMA**2
1215
 
      Y1=1-X1
1216
 
      Y2=1-X2
1217
 
      X0=2.D0-X1-X2
1218
 
      W1=(1.D0-x2)
1219
 
      W2=(1.D0-X0+W-CH)
1220
 
      W22=1.D0/ ((1.D0-X0+W-CH)**2+GAMC)
1221
 
      W11=1.D0/((1.D0-X2)**2+GAMT)
1222
 
      W12=W1*W2*W11*W22
1223
 
      R11=4.D0*T*W-4.D0*T*Y1*Y2+8.D0*T*Y2-4.D0*T*Y2**2+8.D0*W*Y1*Y2+4.D0
1224
 
     .  *W*Y2**2-4.D0*Y1*Y2**2+4.D0*T**2-8.D0*W**2+4.D0*Y2**2-4.D0*Y2**3
1225
 
      R22=-16.D0*W+16.D0*T*W-8.D0*T*Y1*Y2-4.D0*T*Y1**2-4.D0*T*Y2**2+16.
1226
 
     .D0*W*Y1+8.D0*W*Y1*Y2+16.D0*W*Y2+4.D0*W*Y1**2+4.D0*W*Y2**2+8.D0*Y1*
1227
 
     . Y2-12.D0*Y1*Y2**2-12.D0*Y1**2*Y2-16.D0*W**2+4.D0*Y1**2-4.D0*Y1**3
1228
 
     . +4.D0*Y2**2-4.D0*Y2**3
1229
 
      R12=16.D0*W-16.D0*T*W-8.D0*T*Y1+16.D0*T*Y1*Y2-8.D0*T*Y2+8.D0*T*Y1
1230
 
     . **2+8.D0*T*Y2**2-16.D0*W*Y1-8.D0*W*Y1*Y2-16.D0*W*Y2-8.D0*W*Y2**2-
1231
 
     . 8.D0*Y1*Y2+16.D0*Y1*Y2**2+8.D0*Y1**2*Y2+16.D0*W**2-8.D0*Y2**2
1232
 
     . +8.D0*Y2**3
1233
 
      RES=R11*W11+R22*W22+R12*W12
1234
 
      RETURN
1235
 
      END
1236
 
 
1237
 
C     ************************************************************
1238
 
C     SUBROUTINE FOR H ---> TT* ---> TBW
1239
 
C     ************************************************************
1240
 
      SUBROUTINE HTOTT(AMH,AMT,AMB,AMW,AMCH,TB,GHT,GAT,GHVV,HTT0)
1241
 
      IMPLICIT REAL*8(A-Z)
1242
 
      INTEGER IP,K
1243
 
      COMMON/PREC1/IP
1244
 
      EXTERNAL FUNHTT1
1245
 
      COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
1246
 
      COMMON/TOP2/AMH2,AMT2,AMB2,AMW2,AMCH2,TB2,GHT2,GAT2,GHVV2
1247
 
      AMH2=AMH
1248
 
      AMT2=AMT
1249
 
      AMB2=AMB
1250
 
      AMW2=AMW
1251
 
      AMCH2=AMCH
1252
 
      TB2=TB
1253
 
      GHT2=GHT
1254
 
      GAT2=GAT
1255
 
      GHVV2=GHVV
1256
 
      IP=5
1257
 
      M1=AMB
1258
 
      M2=AMT
1259
 
      M3=AMW
1260
 
C     FIRST INTEGRATE OVER X2, i.e. (1+3) SYSTEM
1261
 
C        CHECK WHETHER ENOUGH PHASE SPACE
1262
 
      MASTOT=M1+M2+M3
1263
 
      IF(MASTOT.GE.AMH) GOTO 12
1264
 
      ECM=AMH
1265
 
      S=ECM**2
1266
 
      U1=(ECM-M2)**2
1267
 
      D1=(M1+M3)**2
1268
 
      U=(S-D1+M2**2)/s
1269
 
      D=(S-U1+M2**2)/s
1270
 
      DEL=(U-D)/IP
1271
 
      U=D+DEL
1272
 
      XSEC=0.D0
1273
 
      DO K=1,IP
1274
 
      CALL QGAUS1(FUNHTT1,D,U,SS)
1275
 
      D=U
1276
 
      U=D+DEL
1277
 
      XSEC=XSEC+SS
1278
 
      ENDDO
1279
 
      HTT0=XSEC
1280
 
 12   CONTINUE
1281
 
      RETURN
1282
 
      END
1283
 
 
1284
 
      DOUBLE PRECISION FUNCTION FUNHTT1(XL)
1285
 
      IMPLICIT REAL*8(A-Z)
1286
 
      INTEGER IP,I
1287
 
      COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
1288
 
      COMMON/PREC1/IP
1289
 
      EXTERNAL FUNHTT2
1290
 
      X2=XL
1291
 
      S13=S-S*X2+M2**2
1292
 
      TEM=2.D0*DSQRT(S13)
1293
 
      E2S=(S-S13-M2**2)/TEM
1294
 
      E3S=(S13+M3**2-M1**2)/TEM
1295
 
C     SECOND INTEGRAL OVER X1, i.e. (2+3) SYSTEM
1296
 
      U1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)-DSQRT(E3S**2-M3**2))**2
1297
 
      D1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)+DSQRT(E3S**2-M3**2))**2
1298
 
      U=(S-D1+M1**2)/s
1299
 
      D=(S-U1+M1**2)/s
1300
 
      DEL=(U-D)/IP
1301
 
      FUNHTT1=0.d0
1302
 
      U=D+DEL
1303
 
      DO I=1,IP
1304
 
      CALL QGAUS2(FUNHTT2,D,U,SS)
1305
 
      FUNHTT1=FUNHTT1+SS
1306
 
      D=U
1307
 
      U=D+DEL
1308
 
      ENDDO
1309
 
      RETURN
1310
 
      END
1311
 
 
1312
 
      DOUBLE PRECISION FUNCTION FUNHTT2(XK)
1313
 
      IMPLICIT REAL*8(A-Z)
1314
 
      COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
1315
 
      X1=XK
1316
 
      CALL ELEMHTT(SS)
1317
 
      FUNHTT2=SS
1318
 
      RETURN
1319
 
      END
1320
 
 
1321
 
      SUBROUTINE ELEMHTT(RES)
1322
 
      IMPLICIT REAL*8(A-Z)
1323
 
      COMMON/IKSY2/X1,X2,M1,M2,M3,ECM,S
1324
 
      COMMON/TOP2/AMH,AMT,AMB,AMW,AMCH,TB,GHT,GAT,GHVV
1325
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW0,GAMZ0
1326
 
      GAMT=GAMT1**2*AMT**2/AMH**4
1327
 
      GAMC=GAMC0**2*AMCH**2/AMH**4
1328
 
      GAMW=GAMW0**2*AMW**2/AMH**4
1329
 
      CH=AMCH**2/AMH**2
1330
 
      W=AMW**2/AMH**2
1331
 
      T=AMT**2/AMH**2
1332
 
      Y1=1-X2
1333
 
      Y2=1-X1
1334
 
      X0=2.D0-X1-X2
1335
 
      W1=(1.D0-X2)
1336
 
      W2=(1.D0-X0+W-CH)
1337
 
      W3=-(1.-X1-X2)
1338
 
      W22=1.D0/ ((1.D0-X0+W-CH)**2+GAMC)
1339
 
      W11=1.D0/((1.D0-X2)**2+GAMT)
1340
 
      W33=1.D0/(W3**2+GAMW**2)
1341
 
      W12=W1*W2*W11*W22
1342
 
      W13=W1*W3*W11*W33
1343
 
      W23=W2*W3*W22*W33
1344
 
 
1345
 
      R11=4*T*W-16.*T*W*Y1-4.*T*Y2*Y1+8.*T*Y1+32.*T*W**2-20
1346
 
     . .*T*Y1**2+8.*W*Y2*Y1+4.*W*Y1**2-4.*Y2*Y1**2-16.*T**2*W-
1347
 
     .  32.*T**2*Y1+4.*T**2-16.*T**3-8.*W**2+4.*Y1**2-4.*Y1**3
1348
 
      R22=-16.*W+16.*T*W-8.*T*Y2*Y1-4.*T*Y2**2-4.*T*Y1**2+16
1349
 
     .  .*W*Y2 + 8.*W*Y2*Y1 + 16.*W*Y1 + 4.*W*Y2**2 + 4.*W*Y1**2+8.*Y2*
1350
 
     .  Y1-12.*Y2*Y1**2-12.*Y2**2*Y1-16.*W**2+4.*Y2**2-4.*Y2**3
1351
 
     .  +4.*Y1**2-4.*Y1**3
1352
 
      R33=-4.*T*W+4.*T*W*Y2-2.*T*W*Y2*Y1+4.*T*W*Y1+T*W*Y2**2-
1353
 
     .  3.*T*W*Y1**2+2.*T*Y2*Y1-3.*T*Y2*Y1**2+4.*T*W**2-4.*T*W**3
1354
 
     .  +T*Y2**2-3.*T*Y2**2*Y1-T*Y2**3+T*Y1**2-T*Y1**3+4.*T**2
1355
 
     .  *W-4.*T**2*W*Y2-4.*T**2*W*Y1-2.*T**2*Y2*Y1-4.*T**2*W**2-
1356
 
     .  T**2*Y2**2-T**2*Y1**2+4.*W**2*Y2*Y1-8.*W**3*Y2-8.*W**3*Y1
1357
 
     .  +4.*W**3+8.*W**4
1358
 
      R12=-16.*W+48.*T*W-16.*T*W*Y2+16.*T*W*Y1+8.*T*Y2-32.*T
1359
 
     .  *Y2*Y1+8.*T*Y1-8.*T*Y2**2 - 24.*T*Y1**2+16.*W*Y2+8.*W*Y2*
1360
 
     .  Y1+16.*W*Y1+8.*W*Y1**2+8.*Y2*Y1-16.*Y2*Y1**2-16.*T**2*Y2
1361
 
     .  -16.*T**2*Y1-8.*Y2**2*Y1-16.*W**2+8.*Y1**2-8.*Y1**3
1362
 
      R13=8.*W-24.*T*W+16.*T*W*Y1 -4.*T*Y2+16.*T*Y2*Y1-4.*T*
1363
 
     .  Y1+16.*T*W**2+4.*T*Y2**2+12.*T*Y1**2-8.*W*Y2-12.*W*Y2*Y1
1364
 
     .  -8.*W*Y1+4.*W*Y1**2-4.*Y2*Y1+8.*Y2*Y1**2+16.*T**2*W+8.
1365
 
     .  *T**2*Y2+8.*T**2*Y1+16.*W**2*Y2+24.*W**2*Y1+4.*Y2**2*Y1-
1366
 
     .  32.*W**3-4.*Y1**2+4.*Y1**3
1367
 
      R23=16.*W-16.*T*W+8.*T*W*Y2+8.*T*W*Y1+8.*T*Y2*Y1+4.*T*
1368
 
     .  Y2**2+4.*T*Y1**2-16.*W*Y2-16.*W*Y1-4.*W*Y2**2+4.*W*Y1**2
1369
 
     .  -8.*Y2*Y1+12.*Y2*Y1**2+8.*W**2*Y2-8.*W**2*Y1+12.*Y2**2*
1370
 
     .  Y1-4.*Y2**2+4.*Y2**3-4.*Y1**2+4.*Y1**3
1371
 
      GLVV=DSQRT(1.D0-GHVV**2)
1372
 
      RES=GHT**2*R11*W11+GLVV**2*GAT**2*R22*W22+
1373
 
     .    4.D0*GHVV**2*R33*W33/T+2.D0*GHT*GLVV*GAT*R12*W12+
1374
 
     .    2.D0*GHT*GHVV*R13*W13+2.D0*GHVV*GLVV*GAT*R23*W23
1375
 
      RETURN
1376
 
      END
1377
 
 
1378
 
C     ************************************************************
1379
 
C     SUBROUTINE FOR H+ ---> BT* ---> BBW
1380
 
C     ************************************************************
1381
 
      SUBROUTINE CTOTT(AMCH,AMT,AMB,AMW,CTT0)
1382
 
      IMPLICIT REAL*8(A-Z)
1383
 
      INTEGER IP,K
1384
 
      COMMON/PREC1/IP
1385
 
      EXTERNAL FUNCTT1
1386
 
      COMMON/IKSY3/X1,X2,M1,M2,M3,ECM,S
1387
 
      COMMON/TOP3/AMH3,AMT3,AMB3,AMW3
1388
 
      AMH3=AMCH
1389
 
      AMT3=AMT
1390
 
      AMB3=AMB
1391
 
      AMW3=AMW
1392
 
      IP=5
1393
 
      M1=AMB
1394
 
      M2=AMB
1395
 
      M3=AMW
1396
 
C     FIRST INTEGRATE OVER X2, i.e. (1+3) SYSTEM
1397
 
C        CHECK WHETHER ENOUGH PHASE SPACE
1398
 
      MASTOT=M1+M2+M3
1399
 
      IF(MASTOT.GE.AMCH) GOTO 12
1400
 
      ECM=AMCH
1401
 
      S=ECM**2
1402
 
      U1=(ECM-M2)**2
1403
 
      D1=(M1+M3)**2
1404
 
      U=(S-D1+M2**2)/s
1405
 
      D=(S-U1+M2**2)/s
1406
 
      DEL=(U-D)/IP
1407
 
      U=D+DEL
1408
 
      XSEC=0.D0
1409
 
      DO K=1,IP
1410
 
      CALL QGAUS1(FUNCTT1,D,U,SS)
1411
 
      D=U
1412
 
      U=D+DEL
1413
 
      XSEC=XSEC+SS
1414
 
      ENDDO
1415
 
      CTT0=XSEC
1416
 
12    CONTINUE
1417
 
      RETURN
1418
 
      END
1419
 
 
1420
 
      DOUBLE PRECISION FUNCTION FUNCTT1(XL)
1421
 
      IMPLICIT REAL*8(A-Z)
1422
 
      INTEGER IP,I
1423
 
      COMMON/IKSY3/X1,X2,M1,M2,M3,ECM,S
1424
 
      COMMON/PREC1/IP
1425
 
      EXTERNAL FUNCTT2
1426
 
      X2=XL
1427
 
      S13=S-S*X2+M2**2
1428
 
      TEM=2.D0*DSQRT(S13)
1429
 
      E2S=(S-S13-M2**2)/TEM
1430
 
      E3S=(S13+M3**2-M1**2)/TEM
1431
 
C     SECOND INTEGRAL OVER X1, i.e. (2+3) SYSTEM
1432
 
      U1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)-DSQRT(E3S**2-M3**2))**2
1433
 
      D1=(E2S+E3S)**2-(DSQRT(E2S**2-M2**2)+DSQRT(E3S**2-M3**2))**2
1434
 
      U=(S-D1+M1**2)/s
1435
 
      D=(S-U1+M1**2)/s
1436
 
      DEL=(U-D)/IP
1437
 
      FUNCTT1=0.d0
1438
 
      U=D+DEL
1439
 
      DO I=1,IP
1440
 
      CALL QGAUS2(FUNCTT2,D,U,SS)
1441
 
      FUNCTT1=FUNCTT1+SS
1442
 
      D=U
1443
 
      U=D+DEL
1444
 
      ENDDO
1445
 
      RETURN
1446
 
      END
1447
 
 
1448
 
      DOUBLE PRECISION FUNCTION FUNCTT2(XK)
1449
 
      IMPLICIT REAL*8(A-Z)
1450
 
      COMMON/IKSY3/X1,X2,M1,M2,M3,ECM,S
1451
 
      X1=XK
1452
 
      CALL ELEMCTT(SS)
1453
 
      FUNCTT2=SS
1454
 
      RETURN
1455
 
      END
1456
 
 
1457
 
      SUBROUTINE ELEMCTT(RES)
1458
 
      IMPLICIT REAL*8(A-Z)
1459
 
      COMMON/IKSY3/X1,X2,M1,M2,M3,ECM,S
1460
 
      COMMON/TOP3/AMCH,AMT,AMB,AMW
1461
 
      COMMON/WZWDTH/GAMC0,GAMT0,GAMT1,GAMW,GAMZ
1462
 
      GAMT=GAMT1**2*AMT**2/AMCH**4
1463
 
      W=AMW**2/AMCH**2
1464
 
      T=AMT**2/AMCH**2
1465
 
      B=AMB**2/AMCH**2
1466
 
      RES=((1.D0-X1-W)*(1.D0-X2-W)+W*(X1+X2-1.D0+W))/
1467
 
     .   ((1.D0-X2+B-T)**2+GAMT)
1468
 
      RETURN
1469
 
      END
1470
 
 
1471
 
C   *****************  INTEGRATION ROUTINE ***********************
1472
 
C    Returns SS as integral of FUNC from A to B, by 10-point Gauss-
1473
 
C    Legendre integration
1474
 
      SUBROUTINE QGAUS1(FUNC,A,B,SS)
1475
 
      IMPLICIT REAL*8(A-Z)
1476
 
      INTEGER J
1477
 
      DIMENSION X(5),W(5)
1478
 
      EXTERNAL FUNC
1479
 
      DATA X/.1488743389D0,.4333953941D0,.6794095682D0
1480
 
     .  ,.8650633666D0,.9739065285D0/
1481
 
      DATA W/.2955242247D0,.2692667193D0,.2190863625D0
1482
 
     .  ,.1494513491D0,.0666713443D0/
1483
 
      XM=0.5D0*(B+A)
1484
 
      XR=0.5D0*(B-A)
1485
 
      SS=0.D0
1486
 
      DO 11 J=1,5
1487
 
        DX=XR*X(J)
1488
 
        SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
1489
 
11    CONTINUE
1490
 
      SS=XR*SS
1491
 
      RETURN
1492
 
      END
1493
 
 
1494
 
C     Returns SS as integral of FUNC from A to B, by 10-point Gauss-
1495
 
C      Legendre integration
1496
 
      SUBROUTINE QGAUS2(FUNC,A,B,SS)
1497
 
      IMPLICIT REAL*8(A-Z)
1498
 
      INTEGER J
1499
 
      DIMENSION X(5),W(5)
1500
 
      EXTERNAL FUNC
1501
 
      DATA X/.1488743389D0,.4333953941D0,.6794095682D0
1502
 
     .  ,.8650633666D0,.9739065285D0/
1503
 
      DATA W/.2955242247D0,.2692667193D0,.2190863625D0
1504
 
     .  ,.1494513491D0,.0666713443D0/
1505
 
      XM=0.5D0*(B+A)
1506
 
      XR=0.5D0*(B-A)
1507
 
      SS=0.D0
1508
 
      DO 11 J=1,5
1509
 
        DX=XR*X(J)
1510
 
        SS=SS+W(J)*(FUNC(XM+DX)+FUNC(XM-DX))
1511
 
11    CONTINUE
1512
 
      SS=XR*SS
1513
 
      RETURN
1514
 
      END
1515
 
 
1516
 
 
1517
 
C ******************************************************************
1518
 
 
1519
 
C      DOUBLE PRECISION FUNCTION RUNP(Q,NF)
1520
 
C      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1521
 
C      COMMON/RUN/XMSB,XMHAT,XKFAC
1522
 
C      RUNP = RUNM(Q,NF)
1523
 
C      RUNP = RUNM(Q/2.D0,NF)*XKFAC
1524
 
C      RETURN
1525
 
C      END
1526
 
 
1527
 
      DOUBLE PRECISION FUNCTION RUNM(Q,NF)
1528
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1529
 
      PARAMETER (NN=6)
1530
 
      PARAMETER (ZETA3 = 1.202056903159594D0)
1531
 
      DIMENSION AM(NN),YMSB(NN)
1532
 
      COMMON/ALS/XLAMBDA,AMCA,AMBA,AMTA,N0A
1533
 
      COMMON/MASSES/AMS,AMC,AMB,AMT
1534
 
      COMMON/STRANGE/AMSB
1535
 
      COMMON/RUN/XMSB,XMHAT,XKFAC
1536
 
      COMMON/FLAG/IHIGGS,NNLO,IPOLE
1537
 
      SAVE ISTRANGE
1538
 
      B0(NF)=(33.D0-2.D0*NF)/12D0
1539
 
      B1(NF) = (102D0-38D0/3D0*NF)/16D0
1540
 
      B2(NF) = (2857D0/2D0-5033D0/18D0*NF+325D0/54D0*NF**2)/64D0
1541
 
      G0(NF) = 1D0
1542
 
      G1(NF) = (202D0/3D0-20D0/9D0*NF)/16D0
1543
 
      G2(NF) = (1249D0-(2216D0/27D0+160D0/3D0*ZETA3)*NF
1544
 
     .       - 140D0/81D0*NF**2)/64D0
1545
 
      C1(NF) = G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2
1546
 
      C2(NF) = ((G1(NF)/B0(NF) - B1(NF)*G0(NF)/B0(NF)**2)**2
1547
 
     .       + G2(NF)/B0(NF) + B1(NF)**2*G0(NF)/B0(NF)**3
1548
 
     .       - B1(NF)*G1(NF)/B0(NF)**2 - B2(NF)*G0(NF)/B0(NF)**2)/2D0
1549
 
      TRAN(X,XK)=1D0+4D0/3D0*ALPHAS_HD(X,2)/PI+XK*(ALPHAS_HD(X,2)/PI)**2
1550
 
      CQ(X,NF)=(2D0*B0(NF)*X)**(G0(NF)/B0(NF))
1551
 
     .            *(1D0+C1(NF)*X+C2(NF)*X**2)
1552
 
      DATA ISTRANGE/0/
1553
 
      PI=4D0*DATAN(1D0)
1554
 
      ACC = 1.D-8
1555
 
      AM(1) = 0
1556
 
      AM(2) = 0
1557
 
C--------------------------------------------
1558
 
      IMSBAR = 0
1559
 
      IF(IMSBAR.EQ.1)THEN
1560
 
       IF(ISTRANGE.EQ.0)THEN
1561
 
C--STRANGE POLE MASS FROM MSBAR-MASS AT 1 GEV
1562
 
        AMSD = XLAMBDA
1563
 
        AMSU = 1.D8
1564
 
123     AMS  = (AMSU+AMSD)/2
1565
 
        AM(3) = AMS
1566
 
        XMSB = AMS/CQ(ALPHAS_HD(AMS,2)/PI,3)
1567
 
     .            *CQ(ALPHAS_HD(1.D0,2)/PI,3)/TRAN(AMS,0D0)
1568
 
        DD = (XMSB-AMSB)/AMSB
1569
 
        IF(DABS(DD).GE.ACC)THEN
1570
 
         IF(DD.LE.0.D0)THEN
1571
 
          AMSD = AM(3)
1572
 
         ELSE
1573
 
          AMSU = AM(3)
1574
 
         ENDIF
1575
 
         GOTO 123
1576
 
        ENDIF
1577
 
        ISTRANGE=1
1578
 
       ENDIF
1579
 
       AM(3) = AMSB
1580
 
      ELSE
1581
 
       AMS=AMSB
1582
 
       AM(3) = AMS
1583
 
      ENDIF
1584
 
C--------------------------------------------
1585
 
      AM(3) = AMSB
1586
 
      AM(4) = AMC
1587
 
      AM(5) = AMB
1588
 
      AM(6) = AMT
1589
 
      XK = 16.11D0
1590
 
      DO 1 I=1,NF-1
1591
 
       XK = XK - 1.04D0*(1.D0-AM(I)/AM(NF))
1592
 
1     CONTINUE
1593
 
      IF(NF.GE.4)THEN
1594
 
       XMSB = AM(NF)/TRAN(AM(NF),0D0)
1595
 
       XMHAT = XMSB/CQ(ALPHAS_HD(AM(NF),2)/PI,NF)
1596
 
      ELSE
1597
 
       XMSB = 0
1598
 
       XMHAT = 0
1599
 
      ENDIF
1600
 
      YMSB(3) = AMSB
1601
 
      IF(NF.EQ.3)THEN
1602
 
       YMSB(4) = YMSB(3)*CQ(ALPHAS_HD(AM(4),2)/PI,3)/
1603
 
     .                   CQ(ALPHAS_HD(1.D0,2)/PI,3)
1604
 
       YMSB(5) = YMSB(4)*CQ(ALPHAS_HD(AM(5),2)/PI,4)/
1605
 
     .                   CQ(ALPHAS_HD(AM(4),2)/PI,4)
1606
 
       YMSB(6) = YMSB(5)*CQ(ALPHAS_HD(AM(6),2)/PI,5)/
1607
 
     .                   CQ(ALPHAS_HD(AM(5),2)/PI,5)
1608
 
      ELSEIF(NF.EQ.4)THEN
1609
 
       YMSB(4) = XMSB
1610
 
       YMSB(5) = YMSB(4)*CQ(ALPHAS_HD(AM(5),2)/PI,4)/
1611
 
     .                   CQ(ALPHAS_HD(AM(4),2)/PI,4)
1612
 
       YMSB(6) = YMSB(5)*CQ(ALPHAS_HD(AM(6),2)/PI,5)/
1613
 
     .                   CQ(ALPHAS_HD(AM(5),2)/PI,5)
1614
 
      ELSEIF(NF.EQ.5)THEN
1615
 
       YMSB(5) = XMSB
1616
 
       YMSB(4) = YMSB(5)*CQ(ALPHAS_HD(AM(4),2)/PI,4)/
1617
 
     .                   CQ(ALPHAS_HD(AM(5),2)/PI,4)
1618
 
       YMSB(6) = YMSB(5)*CQ(ALPHAS_HD(AM(6),2)/PI,5)/
1619
 
     .                   CQ(ALPHAS_HD(AM(5),2)/PI,5)
1620
 
      ELSEIF(NF.EQ.6)THEN
1621
 
       YMSB(6) = XMSB
1622
 
       YMSB(5) = YMSB(6)*CQ(ALPHAS_HD(AM(5),2)/PI,5)/
1623
 
     .                   CQ(ALPHAS_HD(AM(6),2)/PI,5)
1624
 
       YMSB(4) = YMSB(5)*CQ(ALPHAS_HD(AM(4),2)/PI,4)/
1625
 
     .                   CQ(ALPHAS_HD(AM(5),2)/PI,4)
1626
 
      ENDIF
1627
 
      IF(Q.LT.AMC)THEN
1628
 
       N0=3
1629
 
       Q0 = 1.D0
1630
 
      ELSEIF(Q.LE.AMB)THEN
1631
 
       N0=4
1632
 
       Q0 = AMC
1633
 
      ELSEIF(Q.LE.AMT)THEN
1634
 
       N0=5
1635
 
       Q0 = AMB
1636
 
      ELSE
1637
 
       N0=6
1638
 
       Q0 = AMT
1639
 
      ENDIF
1640
 
      IF(NNLO.EQ.1.AND.NF.GT.3)THEN
1641
 
       XKFAC = TRAN(AM(NF),0D0)/TRAN(AM(NF),XK)
1642
 
      ELSE
1643
 
       XKFAC = 1D0
1644
 
      ENDIF
1645
 
      RUNM = YMSB(N0)*CQ(ALPHAS_HD(Q,2)/PI,N0)/
1646
 
     .               CQ(ALPHAS_HD(Q0,2)/PI,N0)
1647
 
     .       * XKFAC
1648
 
      RETURN
1649
 
      END
1650
 
 
1651
 
      DOUBLE PRECISION FUNCTION ALPHAS_HD(Q,N)
1652
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1653
 
      DIMENSION XLB(6)
1654
 
      COMMON/ALSLAM/XLB1(6),XLB2(6)
1655
 
      COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
1656
 
      B0(NF)=33.D0-2.D0*NF
1657
 
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1658
 
      ALS1(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
1659
 
      ALS2(NF,X)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB(NF)**2))
1660
 
     .          *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB(NF)**2))
1661
 
     .           /DLOG(X**2/XLB(NF)**2))
1662
 
      PI=4.D0*DATAN(1.D0)
1663
 
      IF(N.EQ.1)THEN
1664
 
       DO 1 I=1,6
1665
 
        XLB(I)=XLB1(I)
1666
 
1      CONTINUE
1667
 
      ELSE
1668
 
       DO 2 I=1,6
1669
 
        XLB(I)=XLB2(I)
1670
 
2      CONTINUE
1671
 
      ENDIF
1672
 
      IF(Q.LT.AMC)THEN
1673
 
       NF=3
1674
 
      ELSEIF(Q.LE.AMB)THEN
1675
 
       NF=4
1676
 
      ELSEIF(Q.LE.AMT)THEN
1677
 
       NF=5
1678
 
      ELSE
1679
 
       NF=6
1680
 
      ENDIF
1681
 
      IF(N.EQ.1)THEN
1682
 
        ALPHAS_HD=ALS1(NF,Q)
1683
 
      ELSE
1684
 
        ALPHAS_HD=ALS2(NF,Q)
1685
 
      ENDIF
1686
 
      RETURN
1687
 
      END
1688
 
 
1689
 
      SUBROUTINE ALSINI(ACC)
1690
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1691
 
      DIMENSION XLB(6)
1692
 
      COMMON/ALSLAM/XLB1(6),XLB2(6)
1693
 
      COMMON/ALS/XLAMBDA,AMC,AMB,AMT,N0
1694
 
      PI=4.D0*DATAN(1.D0)
1695
 
      XLB1(1)=0D0
1696
 
      XLB1(2)=0D0
1697
 
      XLB2(1)=0D0
1698
 
      XLB2(2)=0D0
1699
 
      IF(N0.EQ.3)THEN
1700
 
       XLB(3)=XLAMBDA
1701
 
       XLB(4)=XLB(3)*(XLB(3)/AMC)**(2.D0/25.D0)
1702
 
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1703
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1704
 
      ELSEIF(N0.EQ.4)THEN
1705
 
       XLB(4)=XLAMBDA
1706
 
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1707
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1708
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1709
 
      ELSEIF(N0.EQ.5)THEN
1710
 
       XLB(5)=XLAMBDA
1711
 
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1712
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1713
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1714
 
      ELSEIF(N0.EQ.6)THEN
1715
 
       XLB(6)=XLAMBDA
1716
 
       XLB(5)=XLB(6)*(XLB(6)/AMT)**(-2.D0/23.D0)
1717
 
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1718
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1719
 
      ENDIF
1720
 
      DO 1 I=1,6
1721
 
       XLB1(I)=XLB(I)
1722
 
1     CONTINUE
1723
 
      IF(N0.EQ.3)THEN
1724
 
       XLB(3)=XLAMBDA
1725
 
       XLB(4)=XLB(3)*(XLB(3)/AMC)**(2.D0/25.D0)
1726
 
     .             *(2.D0*DLOG(AMC/XLB(3)))**(-107.D0/1875.D0)
1727
 
       XLB(4)=XITER(AMC,XLB(3),3,XLB(4),4,ACC)
1728
 
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1729
 
     .             *(2.D0*DLOG(AMB/XLB(4)))**(-963.D0/13225.D0)
1730
 
       XLB(5)=XITER(AMB,XLB(4),4,XLB(5),5,ACC)
1731
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1732
 
     .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1733
 
       XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1734
 
      ELSEIF(N0.EQ.4)THEN
1735
 
       XLB(4)=XLAMBDA
1736
 
       XLB(5)=XLB(4)*(XLB(4)/AMB)**(2.D0/23.D0)
1737
 
     .             *(2.D0*DLOG(AMB/XLB(4)))**(-963.D0/13225.D0)
1738
 
       XLB(5)=XITER(AMB,XLB(4),4,XLB(5),5,ACC)
1739
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1740
 
     .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1741
 
       XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1742
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1743
 
     .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1744
 
       XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1745
 
      ELSEIF(N0.EQ.5)THEN
1746
 
       XLB(5)=XLAMBDA
1747
 
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1748
 
     .             *(2.D0*DLOG(AMB/XLB(5)))**(963.D0/14375.D0)
1749
 
       XLB(4)=XITER(AMB,XLB(5),5,XLB(4),4,ACC)
1750
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1751
 
     .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1752
 
       XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1753
 
       XLB(6)=XLB(5)*(XLB(5)/AMT)**(2.D0/21.D0)
1754
 
     .            *(2.D0*DLOG(AMT/XLB(5)))**(-321.D0/3381.D0)
1755
 
       XLB(6)=XITER(AMT,XLB(5),5,XLB(6),6,ACC)
1756
 
      ELSEIF(N0.EQ.6)THEN
1757
 
       XLB(6)=XLAMBDA
1758
 
       XLB(5)=XLB(6)*(XLB(6)/AMT)**(-2.D0/23.D0)
1759
 
     .            *(2.D0*DLOG(AMT/XLB(6)))**(321.D0/3703.D0)
1760
 
       XLB(5)=XITER(AMT,XLB(6),6,XLB(5),5,ACC)
1761
 
       XLB(4)=XLB(5)*(XLB(5)/AMB)**(-2.D0/25.D0)
1762
 
     .             *(2.D0*DLOG(AMB/XLB(5)))**(963.D0/14375.D0)
1763
 
       XLB(4)=XITER(AMB,XLB(5),5,XLB(4),4,ACC)
1764
 
       XLB(3)=XLB(4)*(XLB(4)/AMC)**(-2.D0/27.D0)
1765
 
     .             *(2.D0*DLOG(AMC/XLB(4)))**(107.D0/2025.D0)
1766
 
       XLB(3)=XITER(AMC,XLB(4),4,XLB(3),3,ACC)
1767
 
      ENDIF
1768
 
      DO 2 I=1,6
1769
 
       XLB2(I)=XLB(I)
1770
 
2     CONTINUE
1771
 
      RETURN
1772
 
      END
1773
 
 
1774
 
      DOUBLE PRECISION FUNCTION XITER(Q,XLB1,NF1,XLB,NF2,ACC)
1775
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1776
 
      B0(NF)=33.D0-2.D0*NF
1777
 
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1778
 
      ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
1779
 
     .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
1780
 
     .              /DLOG(X**2/XLB**2))
1781
 
      AA(NF)=12D0*PI/B0(NF)
1782
 
      BB(NF)=B1(NF)/AA(NF)
1783
 
      XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
1784
 
      PI=4.D0*DATAN(1.D0)
1785
 
      XLB2=XLB
1786
 
      II=0
1787
 
1     II=II+1
1788
 
      X=DLOG(Q**2/XLB2**2)
1789
 
      ALP=ALS2(NF1,Q,XLB1)
1790
 
      A=AA(NF2)/ALP
1791
 
      B=BB(NF2)*ALP
1792
 
      XX=XIT(A,B,X)
1793
 
      XLB2=Q*DEXP(-XX/2.D0)
1794
 
      Y1=ALS2(NF1,Q,XLB1)
1795
 
      Y2=ALS2(NF2,Q,XLB2)
1796
 
      DY=DABS(Y2-Y1)/Y1
1797
 
      IF(DY.GE.ACC) GOTO 1
1798
 
      XITER=XLB2
1799
 
      RETURN
1800
 
      END
1801
 
 
1802
 
      DOUBLE PRECISION FUNCTION FINT(Z,XX,YY)
1803
 
C--ONE-DIMENSIONAL CUBIC INTERPOLATION
1804
 
C--Z  = WANTED POINT
1805
 
C--XX = ARRAY OF 4 DISCRETE X-VALUES AROUND Z
1806
 
C--YY = ARRAY OF 4 DISCRETE FUNCTION-VALUES AROUND Z
1807
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1808
 
      DIMENSION XX(4),YY(4)
1809
 
      X = DLOG(Z)
1810
 
      X0=DLOG(XX(1))
1811
 
      X1=DLOG(XX(2))
1812
 
      X2=DLOG(XX(3))
1813
 
      X3=DLOG(XX(4))
1814
 
      Y0=DLOG(YY(1))
1815
 
      Y1=DLOG(YY(2))
1816
 
      Y2=DLOG(YY(3))
1817
 
      Y3=DLOG(YY(4))
1818
 
      A0=(X-X1)*(X-X2)*(X-X3)/(X0-X1)/(X0-X2)/(X0-X3)
1819
 
      A1=(X-X0)*(X-X2)*(X-X3)/(X1-X0)/(X1-X2)/(X1-X3)
1820
 
      A2=(X-X0)*(X-X1)*(X-X3)/(X2-X0)/(X2-X1)/(X2-X3)
1821
 
      A3=(X-X0)*(X-X1)*(X-X2)/(X3-X0)/(X3-X1)/(X3-X2)
1822
 
      FINT=DEXP(A0*Y0+A1*Y1+A2*Y2+A3*Y3)
1823
 
      RETURN
1824
 
      END
1825
 
 
1826
 
      DOUBLE PRECISION FUNCTION SP(X)
1827
 
C--REAL DILOGARITHM (SPENCE-FUNCTION)
1828
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1829
 
      COMPLEX*16 CX,LI2
1830
 
      CX = DCMPLX(X,0.D0)
1831
 
      SP = DREAL(LI2(CX))
1832
 
      RETURN
1833
 
      END
1834
 
 
1835
 
      COMPLEX*16 FUNCTION LI2(X)
1836
 
C--COMPLEX DILOGARITHM (SPENCE-FUNCTION)
1837
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1838
 
      COMPLEX*16 X,Y,CLI2
1839
 
      COMMON/CONST/ZETA2,ZETA3
1840
 
      ZERO=1.D-16
1841
 
      XR=DREAL(X)
1842
 
      XI=DIMAG(X)
1843
 
      R2=XR*XR+XI*XI
1844
 
      LI2=0
1845
 
      IF(R2.LE.ZERO)THEN
1846
 
        LI2=X
1847
 
        RETURN
1848
 
      ENDIF
1849
 
      RR=XR/R2
1850
 
      IF(R2.EQ.1.D0.AND.XI.EQ.0.D0)THEN
1851
 
        IF(XR.EQ.1.D0)THEN
1852
 
          LI2=DCMPLX(ZETA2)
1853
 
        ELSE
1854
 
          LI2=-DCMPLX(ZETA2/2.D0)
1855
 
        ENDIF
1856
 
        RETURN
1857
 
      ELSEIF(R2.GT.1.D0.AND.RR.GT.0.5D0)THEN
1858
 
        Y=(X-1.D0)/X
1859
 
        LI2=CLI2(Y)+ZETA2-CDLOG(X)*CDLOG(1.D0-X)+0.5D0*CDLOG(X)**2
1860
 
        RETURN
1861
 
      ELSEIF(R2.GT.1.D0.AND.RR.LE.0.5D0)THEN
1862
 
        Y=1.D0/X
1863
 
        LI2=-CLI2(Y)-ZETA2-0.5D0*CDLOG(-X)**2
1864
 
        RETURN
1865
 
      ELSEIF(R2.LE.1.D0.AND.XR.GT.0.5D0)THEN
1866
 
        Y=1.D0-X
1867
 
        LI2=-CLI2(Y)+ZETA2-CDLOG(X)*CDLOG(1.D0-X)
1868
 
       RETURN
1869
 
      ELSEIF(R2.LE.1.D0.AND.XR.LE.0.5D0)THEN
1870
 
        Y=X
1871
 
        LI2=CLI2(Y)
1872
 
        RETURN
1873
 
      ENDIF
1874
 
      END
1875
 
 
1876
 
      COMPLEX*16 FUNCTION CLI2(X)
1877
 
C--TAYLOR-EXPANSION FOR COMPLEX DILOGARITHM (SPENCE-FUNCTION)
1878
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1879
 
      COMPLEX*16 X,Z
1880
 
      COMMON/BERNOULLI/B2(18),B12(18),B3(18)
1881
 
      COMMON/POLY/NBER
1882
 
      N=NBER-1
1883
 
      Z=-CDLOG(1.D0-X)
1884
 
      CLI2=B2(NBER)
1885
 
      DO 111 I=N,1,-1
1886
 
        CLI2=Z*CLI2+B2(I)
1887
 
111   CONTINUE
1888
 
      CLI2=Z**2*CLI2+Z
1889
 
      RETURN
1890
 
      END
1891
 
 
1892
 
      DOUBLE PRECISION FUNCTION FACULT(N)
1893
 
C--DOUBLE PRECISION VERSION OF FACULTY
1894
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1895
 
      FACULT=1.D0
1896
 
      IF(N.EQ.0)RETURN
1897
 
      DO 999 I=1,N
1898
 
        FACULT=FACULT*DFLOAT(I)
1899
 
999   CONTINUE
1900
 
      RETURN
1901
 
      END
1902
 
 
1903
 
      SUBROUTINE BERNINI(N)
1904
 
C--INITIALIZATION OF COEFFICIENTS FOR POLYLOGARITHMS
1905
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1906
 
      DIMENSION B(18),PB(19)
1907
 
      COMMON/BERNOULLI/B2(18),B12(18),B3(18)
1908
 
      COMMON/CONST/ZETA2,ZETA3
1909
 
      COMMON/POLY/NBER
1910
 
 
1911
 
      NBER=N
1912
 
      PI=4.D0*DATAN(1.D0)
1913
 
 
1914
 
      B(1)=-1.D0/2.D0
1915
 
      B(2)=1.D0/6.D0
1916
 
      B(3)=0.D0
1917
 
      B(4)=-1.D0/30.D0
1918
 
      B(5)=0.D0
1919
 
      B(6)=1.D0/42.D0
1920
 
      B(7)=0.D0
1921
 
      B(8)=-1.D0/30.D0
1922
 
      B(9)=0.D0
1923
 
      B(10)=5.D0/66.D0
1924
 
      B(11)=0.D0
1925
 
      B(12)=-691.D0/2730.D0
1926
 
      B(13)=0.D0
1927
 
      B(14)=7.D0/6.D0
1928
 
      B(15)=0.D0
1929
 
      B(16)=-3617.D0/510.D0
1930
 
      B(17)=0.D0
1931
 
      B(18)=43867.D0/798.D0
1932
 
      ZETA2=PI**2/6.D0
1933
 
      ZETA3=1.202056903159594D0
1934
 
 
1935
 
      DO 995 I=1,18
1936
 
        B2(I)=B(I)/FACULT(I+1)
1937
 
        B12(I)=DFLOAT(I+1)/FACULT(I+2)*B(I)/2.D0
1938
 
        PB(I+1)=B(I)
1939
 
        B3(I)=0.D0
1940
 
995   CONTINUE
1941
 
      PB(1)=1.D0
1942
 
      DO 996 I=1,18
1943
 
      DO 996 J=0,I
1944
 
        B3(I)=B3(I)+PB(J+1)*PB(I-J+1)/FACULT(I-J)/FACULT(J+1)
1945
 
     .                                            /DFLOAT(I+1)
1946
 
996   CONTINUE
1947
 
 
1948
 
      RETURN
1949
 
      END
1950
 
 
1951
 
      DOUBLE PRECISION FUNCTION QQINT(RAT,H1,H2)
1952
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1953
 
      N = 2
1954
 
      QQINT = RAT**N * H1 + (1-RAT**N) * H2
1955
 
      RETURN
1956
 
      END
1957
 
 
1958
 
      DOUBLE PRECISION FUNCTION XITLA(NO,ALP,ACC)
1959
 
C--ITERATION ROUTINE TO DETERMINE IMPROVED LAMBDA'S
1960
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1961
 
      COMMON/PARAM/GF,ALPH,AMTAU,AMMUON,AMZ,AMW
1962
 
      B0(NF)=33.D0-2.D0*NF
1963
 
      B1(NF)=6.D0*(153.D0-19.D0*NF)/B0(NF)**2
1964
 
      ALS2(NF,X,XLB)=12.D0*PI/(B0(NF)*DLOG(X**2/XLB**2))
1965
 
     .              *(1.D0-B1(NF)*DLOG(DLOG(X**2/XLB**2))
1966
 
     .              /DLOG(X**2/XLB**2))
1967
 
      AA(NF)=12D0*PI/B0(NF)
1968
 
      BB(NF)=B1(NF)/AA(NF)
1969
 
      XIT(A,B,X)=A/2.D0*(1D0+DSQRT(1D0-4D0*B*DLOG(X)))
1970
 
      PI=4.D0*DATAN(1.D0)
1971
 
      NF=5
1972
 
      Q=AMZ
1973
 
      XLB=Q*DEXP(-AA(NF)/ALP/2.D0)
1974
 
      IF(NO.EQ.1)GOTO 111
1975
 
      II=0
1976
 
1     II=II+1
1977
 
      X=DLOG(Q**2/XLB**2)
1978
 
      A=AA(NF)/ALP
1979
 
      B=BB(NF)*ALP
1980
 
      XX=XIT(A,B,X)
1981
 
      XLB=Q*DEXP(-XX/2.D0)
1982
 
      Y1=ALP
1983
 
      Y2=ALS2(NF,Q,XLB)
1984
 
      DY=DABS(Y2-Y1)/Y1
1985
 
      IF(DY.GE.ACC) GOTO 1
1986
 
111   XITLA=XLB
1987
 
      RETURN
1988
 
      END
1989
 
 
1990
 
 
1991