« back to all changes in this revision
Viewing changes to MG4_DECAY/decay_couplings.f
- browse files at revision 36.21.85
- compare with another revision
- download diff
- download tarball
- view history from revision 36.21.85
- Committer: Chia-Hsien Shen
- Date: 2012-06-30 05:53:02 UTC
- mto: (36.21.97 decay_calculator_205) (238.3.11 1.6.0)
- mto: This revision was merged to the branch mainline in revision 249.
- Revision ID: anakin8120@gmail.com-20120630055302-5azmh3tx4vv2hlp1
move ./decay to ./mg5decay; resolve unit tests (n.b. __init__ does not check keys of input dictionaries, followed last revision)
- files added:
- DECAY
- mg5decay
- files removed:
- MG4_DECAY
- files modified:
- tests/input_files/full_sm_UFO/model.pkl
added
removed
MG4_DECAY/decay_couplings.f
1
c----------------------------------------------------------------------
2
C couplings.f
3
c----------------------------------------------------------------------
4
c This files takes the inputs of the standard model from a Les Houches
5
c file (param_card.dat) and calculates all the couplings that end up
6
c in the Feynman rules, i.e., in the HELAS routines.
7
c
8
c With the convention of the New Web Generation MadEvent in this
9
c part no non-trivial computation is made. The SM model secondary
10
c parameters have been all calculated by the SM calculator, SMCalc
11
c which produced param_card.dat.
12
c
13
c The only exception to the above rule is for alpha_S. In the case
14
c where a pdf is used in the initial state, then the values as(MZ)
15
c set in the les houches card is superseeded.
16
c Scales and pdf's are set in the run_card.dat.
17
c
18
c This file contains the following routines:
19
c
20
c- madgraph original call to set the parameters
21
c- lh_readin is called from here.
22
c It talks to the lh_readin through the common block values.
23
c subroutine setpara
24
c
25
c-routine to read the LH parameters in
26
c subroutine lh_readin
27
c
28
c-to set
29
c subroutine set_it(npara,ivalue,value,name,id,
30
c subroutine case_trap(string,length)
31
c subroutine no_spaces(buff,nchars)
32
c----------------------------------------------------------------------
33
34
35
subroutine setpara
36
c***********************************************************************
37
c This subroutine sets up the HELAS couplings of the STANDARD MODEL.
38
c***********************************************************************
39
implicit none
40
c
41
c local
42
c
43
integer i
44
real*8 dum
45
c
46
c common file with the couplings
47
c
48
include 'coupl.inc'
49
include 'decay.inc'
50
include 'calc_values.inc'
51
c
52
c local
53
c
54
double precision v
55
double precision ee, ee2, ez, ey, sw, cw, sc2
56
double precision gwne, gwud, lambda, lam4, xt, rew, rqcd
57
double precision alphas, alfa, alfaw, mfrun
58
external alphas, alfa, alfaw, mfrun
59
c
60
c Common to lh_readin and printout
61
c
62
double precision alpha, sin2w, gfermi, alfas
63
double precision mtMS,mbMS,mcMS,mtaMS!MSbar masses
64
double precision Vud,Vus !CKM matrix elements
65
common/values/ alpha,sin2w,gfermi,alfas,
66
& mtMS,mbMS,mcMS,mtaMS,
67
& Vud
68
c
69
c constants
70
c
71
double complex ci
72
parameter( ci = ( 0.0d0, 1.0d0 ) )
73
c
74
c alfas and run
75
c
76
include 'alfas.inc'
77
include 'run.inc'
78
c
79
c Auxiliary functions
80
c
81
REAL*8 XX,Y,Z
82
REAL*8 HFF,BETA,HEAVY
83
84
BETA(XX) = DSQRT(One-Four*XX)
85
HFF(XX,Y) = One/(Two*Four*PI)*XX*(BETA(Y))**3
86
HEAVY(XX,Y,Z)= ( One - Half*(y**2+z**2) - Half*(y**2-z**2)**2
87
& + Three*y*z*(XX**2 - One)/(XX**2 + One) )
88
& * sqrt( (One-y**2-z**2)**2 - Four * y**2 * z**2 )
89
c
90
c------------------------------------------
91
c Start calculating the couplings for HELAS
92
c------------------------------------------
93
c
94
call lh_readin
95
c
96
c Strong coupling
97
c
98
c As a rule we first check if a pdf has been chosen in the
99
c run_card.dat (which has been already read at this stage).
100
c If there pdfs in the initial state, then the alpha_s(MZ) used
101
c is set to the corresponding value.
102
103
if(scale.le.1d0) scale=zmass
104
105
if(lpp(1).ne.0.or.lpp(2).ne.0) then
106
if(asmz .le.0.01d0) asmz =0.118d0
107
else
108
asmz=alfas !value read from the param_card.dat
109
c nloop=2
110
endif
111
if(nloop.eq.0) nloop=1
112
113
G = DSQRT(4d0*PI*ALPHAS(SCALE)) ! use setting of the param_card.dat @ NLO
114
GG(1) = -G
115
GG(2) = -G
116
c
117
c auxiliary local values
118
c
119
cw = sqrt( One - sin2w )
120
ee2 = alpha * Fourpi
121
sw = sqrt( sin2w )
122
ee = sqrt( ee2 )
123
ez = ee/(sw*cw)
124
ey = ee*(sw/cw)
125
sc2 = sin2w*( One - sin2w )
126
v = Two*wmass*sw/ee ! the wmass is used to calculate v
127
lambda = hmass**2 / (Two * v**2)
128
c
129
c vector boson couplings
130
c
131
gw = ee/sw
132
gwwa = ee
133
gwwz = ee*cw/sw
134
c
135
c gauge & higgs boson coupling constants
136
c
137
gwwh = dcmplx( ee2/sin2w*Half*v, Zero )
138
gzzh = dcmplx( ee2/sc2*Half*v, Zero )
139
ghhh = dcmplx( -hmass**2/v*Three, Zero )
140
gwwhh = dcmplx( ee2/sin2w*Half, Zero )
141
gzzhh = dcmplx( ee2/sc2*Half, Zero)
142
ghhhh = ghhh/v
143
c
144
c fermion-fermion-vector couplings
145
c
146
gal(1) = dcmplx( ee , Zero )
147
gal(2) = dcmplx( ee , Zero )
148
gau(1) = dcmplx( -ee*Two/Three, Zero )
149
gau(2) = dcmplx( -ee*Two/Three, Zero )
150
gad(1) = dcmplx( ee/Three , Zero )
151
gad(2) = dcmplx( ee/Three , Zero )
152
153
gwf(1) = dcmplx( -ee/sqrt(Two*sin2w), Zero )
154
gwf(2) = dcmplx( Zero , Zero )
155
156
gzn(1) = dcmplx( -ez*Half , Zero )
157
gzn(2) = dcmplx( Zero , Zero )
158
gzl(1) = dcmplx( -ez*(-Half + sin2w) , Zero )
159
gzl(2) = dcmplx( -ey , Zero )
160
gzu(1) = dcmplx( -ez*( Half - sin2w*Two/Three), Zero )
161
gzu(2) = dcmplx( ey*Two/Three , Zero )
162
gzd(1) = dcmplx( -ez*(-Half + sin2w/Three) , Zero )
163
gzd(2) = dcmplx( -ey/Three , Zero )
164
c
165
c fermion-fermion-Higgs couplings (complex) hff(2)
166
c
167
c NOTE: the running mass is evaluated @ the same order
168
c nloop of alpha_s set by the PDF choice
169
c
170
171
if(mtMS.gt.1d0) then
172
ghtop(1) = dcmplx( -mtMS/v, Zero )
173
else
174
ghtop(1) = dcmplx( Zero,Zero)
175
endif
176
ghtop(2) = ghtop(1)
177
178
if(mbMS.gt.1d0) then
179
ghbot(1) = dcmplx( -mbMS/v, Zero )
180
else
181
ghbot(1) = dcmplx( Zero, Zero )
182
endif
183
ghbot(2) = ghbot(1)
184
185
if(mcMS.gt.1d0) then
186
ghcha(1) = dcmplx( -mcMS/v, Zero )
187
else
188
ghcha(1) = dcmplx( Zero, Zero )
189
endif
190
ghcha(2) = ghcha(1)
191
192
ghtau(1) = dcmplx( -mtaMS/v, Zero )
193
ghtau(2) = ghtau(1)
194
c
195
c CKM matrix:
196
c symmetric 3x3 matrix, Vud=Vcs, Vus=Vcd Vcb=Vub=0
197
c
198
c >>>>>>>>>>>>>>>***** NOTE****<<<<<<<<<<<<<<<<<<<<<<<<<
199
c these couplings matter only when interaction_CKM.dat
200
c is used to generate all the diagrams with off-diagonal
201
c couplings. The default of MadEvent is a diagonal
202
c CKM matrix.
203
204
Vus=DSQRT(1d0-Vud**2)
205
do i=1,2
206
gwfc(i) = gwf(i)*Vud
207
gwfs(i) = gwf(i)*Vus
208
gwfm(i) =-gwf(i)*Vus
209
enddo
210
211
c---------------------------------------------------------
212
c Set Photon Width to Zero, used by symmetry optimization
213
c---------------------------------------------------------
214
215
awidth = 0d0
216
c
217
c Z boson partial widths
218
c
219
decz = zmass / ( 24.0d0 * Pi )
220
w_z_nn = decz * ( gzn(1)**2 + gzn(2)**2 )
221
w_z_ll = decz * ( gzl(1)**2 + gzl(2)**2 )
222
decz = decz * Three
223
w_z_uu = decz * ( gzu(1)**2 + gzu(2)**2 )
224
w_z_dd = decz * ( gzd(1)**2 + gzd(2)**2 )
225
dum = dble( (gzl(2)+gzl(1))/(gzl(2)-gzl(1)) )
226
w_z_tau = w_z_ll * heavy( dum, lmass/zmass, lmass/zmass )
227
dum = dble( (gzu(2)+gzu(1))/(gzu(2)-gzu(1)) )
228
w_z_cc = w_z_uu * heavy( dum, cmass/zmass, cmass/zmass )
229
dum = dble( (gzd(2)+gzd(1))/(gzd(2)-gzd(1)) )
230
w_z_bb = w_z_dd * heavy( dum, bmass/zmass, bmass/zmass )
231
232
zwidth = Three*w_z_nn + Two*w_z_ll + w_z_tau
233
& + Two*w_z_dd + w_z_uu + w_z_cc + w_z_bb
234
235
c
236
c W boson partial widths
237
c
238
decw = wmass / ( 24.0d0 * Pi )
239
w_w_nl = decw * ( gwf(1)**2 + gwf(2)**2 )
240
dum = dble( (gwf(2)+gwf(1))/(gwf(2)-gwf(1)) )
241
w_w_tau = w_w_nl * heavy( dum, lmass/wmass, Zero )
242
w_w_ud = w_w_nl * Three
243
w_w_cs = w_w_ud * heavy( dum, cmass/wmass, Zero )
244
245
wwidth = Two*w_w_nl + w_w_tau + w_w_ud + w_w_cs
246
247
c
248
c top quark width
249
c
250
call topwid(tmass,wmass,bmass,wwidth,gw,twidth)
251
252
c
253
c tau width
254
c
255
lwidth = 2.36d-12 !tau width, PDG value
256
c
257
c LO withds of the Higgs into tt~,bb~,tau tau~.
258
c
259
if(hmass.gt.2d0*tmass) then
260
w_h_tt =3d0*cdabs(ghtop(1)**2)*hff(hmass,(tmass/hmass)**2)
261
else
262
w_h_tt =0d0
263
endif
264
265
w_h_bb =3d0*cdabs(ghbot(1)**2)*hff(hmass,(bmass/hmass)**2)
266
w_h_tau=cdabs(ghtau(1)**2)*hff(hmass,(lmass/hmass)**2)
267
w_h_bb =3d0*cdabs(ghbot(1)**2)*hff(hmass,(bmass/hmass)**2)
268
w_h_cc =3d0*cdabs(ghcha(1)**2)*hff(hmass,(cmass/hmass)**2)
269
w_h_tau=cdabs(ghtau(1)**2)*hff(hmass,(lmass/hmass)**2)
270
271
272
if(hmass.gt.2d0*wmass) then
273
w_h_ww=gfermi/8d0/pi/rt2*hmass**3*
274
& dsqrt(1-4d0*(wmass/hmass)**2)*
275
& (12d0*(wmass/hmass)**4 -4d0*(wmass/hmass)**2+1d0)
276
w_h_WLWL=gw**2/64d0/pi*hmass**3/wmass**2* !longitudinal W's
277
& dsqrt(1-4d0*(wmass/hmass)**2)*
278
& (1-2d0*(wmass/hmass)**2)**2
279
else
280
w_h_ww=0d0
281
w_h_WLWL=0d0
282
endif
283
284
if(hmass.gt.2d0*zmass) then
285
w_h_zz=gfermi/16d0/pi/rt2*hmass**3*
286
& dsqrt(1-4d0*(zmass/hmass)**2)*
287
& (12d0*(zmass/hmass)**4 -4d0*(zmass/hmass)**2+1d0)
288
w_h_ZLZL=gw**2/128d0/pi*hmass**3/wmass**2* !longitudinal Z's
289
& dsqrt(1-4d0*(zmass/hmass)**2)*
290
& (1-2d0*(zmass/hmass)**2)**2
291
else
292
w_h_zz=0d0
293
w_h_zLzL=0d0
294
endif
295
296
c--------------------------
297
c start interface to HDECAY
298
c--------------------------
299
c do not change things here unless you exactly know what you are doing
300
c
301
ihiggs = 0
302
nnlo = 0
303
ipole = 0
304
305
ionsh = 0
306
ionwz = 0
307
iofsusy = 1
308
309
hdals = asmz
310
c-- do not set masses to zero here
311
ams = 0.190d0 !strange pole mass
312
313
if(cmass.gt.0d0) then
314
amc = cmass
315
else
316
amc = 1.42d0 !charm pole mass
317
endif
318
319
if(bmass.gt.0d0) then
320
amb = bmass
321
else
322
amb = 4.7d0 !bottom pole mass
323
endif
324
325
if(tmass.gt.0d0) then
326
amt = tmass
327
else
328
amt = 174.3d0 !top pole mass
329
endif
330
331
if(lmass.gt.0d0) then
332
almass = lmass
333
else
334
almass = 1.777d0 !tau mass
335
endif
336
337
ammuon = 0.105658389d0 !muon mass
338
339
alph = 137.0359895D0 !alpha_EM
340
gf = gfermi
341
342
if(wwidth.gt.0d0) then
343
gamw = wwidth
344
else
345
gamw=2.12d0
346
endif
347
348
if(zwidth.gt.0d0) then
349
gamz = zwidth
350
else
351
gamz=2.495d0
352
endif
353
354
amz = zmass
355
amw = wmass
356
357
hdmhbeg = hmass
358
hdmhend = hmass
359
360
c
361
c this calculates the branching ratios of the Higgs
362
c at the best of our knowledge
363
c
364
call hdecay
365
hwidth = SMWDTH
366
367
c------------------------
368
c end interface to HDECAY
369
c------------------------
370
371
372
c----------------------------
373
c end subroutine coupsm
374
c----------------------------
375
376
377
return
378
end
379
380
381
subroutine lh_readin
382
c----------------------------------------------------------------------
383
c Read the parameters from the lh file
384
c
385
c 1. Input values for the EW sector
386
c 2. Higgs mass and width
387
c 3. Fermion masses (pole and MSbar) and widths
388
c----------------------------------------------------------------------
389
implicit none
390
c
391
c parameters
392
c
393
integer maxpara
394
parameter (maxpara=1000)
395
c
396
c local
397
c
398
integer npara,l1,l2
399
character*132 buff
400
real*8 real_value
401
real*8 value(maxpara)
402
integer ivalue(maxpara),n
403
character*20 name(maxpara),bn,dumstring
404
logical block_found,done,fopened
405
integer i,name_length,idum
406
407
c
408
c block info
409
c
410
character*20 block_name
411
c
412
c Common
413
c
414
include 'coupl.inc'
415
include 'decay.inc'
416
c
417
c Common to lh_readin and printout
418
c
419
double precision alpha, sin2w, gfermi, alfas
420
double precision mtMS,mbMS,mcMS,mtaMS!MSbar masses
421
double precision Vud,Vus !CKM matrix elements
422
common/values/ alpha,sin2w,gfermi,alfas,
423
& mtMS,mbMS,mcMS,mtaMS,
424
& Vud
425
c
426
c----------
427
c start
428
c----------
429
c
430
n=0
431
rewind(lunp)
432
done=.false.
433
434
do while(.not.done)
435
block_found=.false.
436
c
437
c looks for the blocks or for decay
438
c
439
do while(.not.block_found)
440
read(lunp,'(a132)',end=99,err=99) buff
441
c-- change to lower case
442
call case_trap(buff,20)
443
if(buff(1:5).eq.'block') then
444
if(index(buff,"#").ne.0) l1=index(buff,"#")-1
445
block_name=buff(6:min(l1,26))
446
call no_spaces(block_name,name_length)
447
c write(*,*) block_name(1:name_length)
448
block_found=.true.
449
elseif(buff(1:5).eq.'decay') then
450
n=n+1
451
l1=30
452
if(index(buff,"#").ne.0) l1=index(buff,"#")-1 ! ignore comments
453
read(buff(6:l1),*) ivalue(n),value(n)
454
name(n)="decay"
455
endif
456
end do
457
c
458
c
459
460
if(block_found) then
461
do while(.true.)
462
read(lunp,'(a132)',end=99,err=99) buff
463
call case_trap(buff,20)
464
if(buff(1:1).eq.'b'.or.buff(1:1).eq.'d') then
465
backspace lunp
466
exit
467
endif
468
if(buff(1:1).ne.'#') then !if it not a comment
469
n=n+1
470
l1=30
471
if(index(buff,"#").ne.0) l1=index(buff,"#")-1 ! ignore comments
472
c
473
c WARNING:... not all blocks have the same sintax!! You need to change it
474
c depending on the block you are reading
475
c
476
477
if(block_name(1:5).eq."mgckm") then
478
read(buff(1:l1),*) ivalue(n),idum,value(n)
479
elseif (block_name(1:6).eq."spinfo".or.
480
&block_name(1:6).eq."dcinfo") then
481
read(buff(1:l1),*) ivalue(n),dumstring
482
else
483
read(buff(1:l1),*) ivalue(n),value(n)
484
endif
485
name(n)=block_name(1:name_length)
486
c write(*,"(1x,i2,2x,e16.8,1x,a)")
487
c & ivalue(n),value(n),name(n)
488
endif
489
end do ! do while in the block
490
else
491
done=.true.
492
endif
493
end do ! do while the entire file
494
495
496
99 continue
497
498
bn="sminputs"
499
call set_it(n,ivalue,value,name,1,bn,alpha,128.9d0)
500
alpha=1d0/alpha
501
call set_it(n,ivalue,value,name,2,bn,gfermi,0.1166d-4)
502
call set_it(n,ivalue,value,name,3,bn,alfas,0.119d0)
503
call set_it(n,ivalue,value,name,4,bn,zmass,91.188d0)
504
call set_it(n,ivalue,value,name,6,bn,tmass,174.3d0)
505
call set_it(n,ivalue,value,name,7,bn,lmass,1.777d0)
506
bn="mgsmparam"
507
call set_it(n,ivalue,value,name,1,bn,sin2w,0.2312d0)
508
call set_it(n,ivalue,value,name,2,bn,wmass,80.419d0)
509
bn="mgyukawa"
510
call set_it(n,ivalue,value,name,4,bn,mcMS,1.25d0)
511
call set_it(n,ivalue,value,name,5,bn,mbMS,4.2d0)
512
call set_it(n,ivalue,value,name,6,bn,mtMS,174d0)
513
call set_it(n,ivalue,value,name,15,bn,mtaMS,1.777d0)
514
bn="mgckm"
515
call set_it(n,ivalue,value,name,1,bn,vud,1d0)
516
bn="mass"
517
call set_it(n,ivalue,value,name,4,bn,cmass,1.4d0)
518
call set_it(n,ivalue,value,name,5,bn,bmass,4.7d0)
519
call set_it(n,ivalue,value,name,6,bn,tmass,tmass*1d0)
520
call set_it(n,ivalue,value,name,15,bn,lmass,lmass*1d0)
521
call set_it(n,ivalue,value,name,25,bn,hmass,120d0)
522
call set_it(n,ivalue,value,name,23,bn,zmass,zmass*1d0)
523
call set_it(n,ivalue,value,name,24,bn,wmass,wmass*1d0)
524
525
c bn="decay"
526
c call set_it(n,ivalue,value,name,6,bn,twidth,1.5083d0)
527
c call set_it(n,ivalue,value,name,25,bn,hwidth,0.0037d0)
528
c call set_it(n,ivalue,value,name,23,bn,zwidth,2.441d0)
529
c call set_it(n,ivalue,value,name,24,bn,wwidth,2.0476d0)
530
531
write(*,*)
532
write(*,*) ' >>> Widths in param_card.dat are ignored <<<'
533
write(*,*)
534
535
return
536
end
537
538
539
subroutine set_it(npara,ivalue,value,name,id,
540
& block_name,var,def_value)
541
c----------------------------------------------------------------------------------
542
c finds the parameter value in block_name and associate var to it.
543
c If it is not found a default is given.
544
c----------------------------------------------------------------------------------
545
implicit none
546
547
c
548
c parameters
549
c
550
integer maxpara
551
parameter (maxpara=100)
552
c
553
c arguments
554
c
555
integer npara,ivalue(maxpara),id
556
character*20 block_name,name(maxpara)
557
real*8 var,def_value,value(maxpara)
558
c
559
c local
560
c
561
logical found
562
integer i
563
c
564
c start
565
c
566
found=.false.
567
do i=1,npara
568
found = (id.eq.ivalue(i)).and.(name(i).eq.block_name)
569
if(found) then
570
var=value(i)
571
exit
572
endif
573
enddo
574
575
if (.not.found) then
576
c write (*,*) "Warning: parameter not found"
577
c write (*,*) " setting it to default value ",def_value
578
var=def_value
579
endif
580
return
581
582
end
583
584
585
subroutine case_trap(string,length)
586
c**********************************************************
587
c change string to lowercase if the input is not
588
c**********************************************************
589
implicit none
590
c
591
c ARGUMENT
592
c
593
character*(*) string
594
integer length
595
c
596
c LOCAL
597
c
598
integer i,k
599
600
do i=1,length
601
k=ichar(string(i:i))
602
if(k.ge.65.and.k.le.90) then !upper case A-Z
603
k=ichar(string(i:i))+32
604
string(i:i)=char(k)
605
endif
606
enddo
607
608
return
609
end
610
611
612
subroutine no_spaces(buff,nchars)
613
c**********************************************************************
614
c Given buff a buffer of words separated by spaces
615
c returns it where all space are moved to the right
616
c returns also the length of the single word.
617
c maxlength is the length of the buffer
618
c**********************************************************************
619
implicit none
620
c
621
c Constants
622
c
623
integer maxline
624
parameter (maxline=20)
625
character*1 null
626
parameter (null=' ')
627
c
628
c Arguments
629
c
630
character*(maxline) buff
631
integer nchars,maxlength
632
c
633
c Local
634
c
635
integer i,j
636
character*(maxline) temp
637
c-----
638
c Begin Code
639
c-----
640
nchars=0
641
c write (*,*) "buff=",buff(1:maxlength)
642
do i=1,maxline
643
if(buff(i:i).ne.null) then
644
nchars=nchars+1
645
temp(nchars:nchars)=buff(i:i)
646
endif
647
c write(*,*) i,":",buff(1:maxlength),":",temp(1:nchars),":"
648
enddo
649
buff=temp
650
end
651
652
653
SUBROUTINE TOPWID(RMT,RMW,RMB,RGW,GW,RGT)
654
c*************************************************************************
655
c THE TOTAL WEAK DECAY WIDTH OF THE TOP QUARK, INCLUDING
656
c THE EFFECTS OF BOTTOM MASS AND, IF IGW=1, A FINITE W WIDTH.
657
c From James Stirling 6-10-94
658
c
659
c RMT=TOP MASS
660
c RMW=W MASS
661
c RMB=B MASS
662
c RGW=W WIDTH
663
c GW =WEAK COUPLING
664
c
665
c RGT=TOP WIDTH
666
c
667
c*************************************************************************
668
IMPLICIT COMPLEX*16(A-H,O-Z)
669
REAL*8 RMT,RMB,RMW,XW,XB,RGW,RGT,GW
670
*
671
PI=4.*DATAN(1.D0)
672
XGW=dcmplx(GW/2d0/dsqrt(2d0))
673
*
674
XB=RMB/RMT
675
XW=RMW/RMT
676
*
677
RM=XB**2
678
OM=1.+RM-DCMPLX(RMW**2,RMW*RGW)/RMT**2
679
Y1=OM+CDSQRT(OM*OM-4.*RM)
680
Y0=OM-CDSQRT(OM*OM-4.*RM)
681
Z1=2.
682
Z0=2.*CDSQRT(RM)
683
*
684
D0=(-Y0**8+3.*Y0**7*RM+3.*Y0**7-8.*Y0**6*RM-12.*Y0**5*RM**
685
. 2-12.*Y0**5*RM+96.*Y0**4*RM**2-48.*Y0**3*RM**3-48.*Y0**3*
686
. RM**2-128.*Y0**2*RM**3+192.*Y0*RM**4+192.*Y0*RM**3-256.*
687
. RM**4)/(24.*Y0**4*(Y1-Y0))
688
D1=(-Y1**8+3.*Y1**7*RM+3.*Y1**7-8.*Y1**6*RM-12.*Y1**5*RM**
689
. 2-12.*Y1**5*RM+96.*Y1**4*RM**2-48.*Y1**3*RM**3-48.*Y1**3*
690
. RM**2-128.*Y1**2*RM**3+192.*Y1*RM**4+192.*Y1*RM**3-256.*
691
. RM**4)/(24.*Y1**4*(Y1-Y0))
692
A4=(32.*RM**4*(Y1-Y0))/(3.*Y1*Y0*(Y1-Y0))
693
A3=(8.*RM**3*(-3.*Y1**2*Y0*RM-3.*Y1**2*Y0+4.*Y1**2*RM+3.*
694
. Y1*Y0**2*RM+3.*Y1*Y0**2-4.*Y0**2*RM))/(3.*Y1**2*Y0**2*(Y1
695
. -Y0))
696
A2=(8.*RM**3*(2.*Y1**3*Y0**2-3.*Y1**3*Y0*RM-3.*Y1**3*Y0+4.
697
. *Y1**3*RM-2.*Y1**2*Y0**3+3.*Y1*Y0**3*RM+3.*Y1*Y0**3-4.*Y0
698
. **3*RM))/(3.*Y1**3*Y0**3*(Y1-Y0))
699
A1=(2.*RM**2*(3.*Y1**4*Y0**3*RM+3.*Y1**4*Y0**3+8.*Y1**4*Y0
700
. **2*RM-12.*Y1**4*Y0*RM**2-12.*Y1**4*Y0*RM+16.*Y1**4*RM**2
701
. -3.*Y1**3*Y0**4*RM-3.*Y1**3*Y0**4-8.*Y1**2*Y0**4*RM+12.*
702
. Y1*Y0**4*RM**2+12.*Y1*Y0**4*RM-16.*Y0**4*RM**2))/(3.*Y1**
703
. 4*Y0**4*(Y1-Y0))
704
B0=(Y1**3-3.*Y1**2*RM-3.*Y1**2+8.*Y1*RM-Y0**3+3.*Y0**2*RM+
705
. 3.*Y0**2-8.*Y0*RM)/(24.*(Y1-Y0))
706
B1=(Y1+Y0-3.*RM-3.)/24.
707
B2=1./24.
708
*
709
RINT=D0*CDLOG((Z1-Y0)/(Z0-Y0))
710
. -D1*CDLOG((Y1-Z1)/(Y1-Z0))
711
. -A4/3.*(1./Z1**3-1./Z0**3)
712
. -A3/2.*(1./Z1**2-1./Z0**2)
713
. -A2 *(1./Z1 -1./Z0 )
714
. +A1*CDLOG(Z1/Z0)
715
. +B0 *(Z1 -Z0 )
716
. +B1/2.*(Z1**2-Z0**2)
717
. +B2/3.*(Z1**3-Z0**3)
718
*
719
XGW4=XGW**4
720
*
721
* TOTAL WIDTH INCLUDES FLAVOUR & COLOUR FACTORS
722
RGT=RMT**3/(RMW*RGW)*XGW4/(8.*PI**3)*DIMAG(RINT)
723
RGT=9.*RGT
724
RETURN
725
END
Loggerhead is a web-based interface for Breezy
Version: 2.0.1