« back to all changes in this revision
Viewing changes to MG4_DECAY/decay_event.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_event.f
1
subroutine decay_event(P5,wgt,ni,ic)
2
c********************************************************************
3
c Decay particle(s) in the event
4
c
5
c input: P(0,MaxParticles) :four momenta
6
c wgt :event weight
7
c ni :number of particle in the event
8
c before the decay
9
c ic(6,MaxParticles):particle labels
10
c
11
c output:P(0,MaxParticles) :four momenta
12
c wgt :event weight
13
c ic(6,MaxParticles):particle labels
14
c
15
c which particle has to be decayed is passed through
16
c a common block DEC_ID in decay.inc
17
c
18
c********************************************************************
19
implicit none
20
include 'decay.inc'
21
c
22
c Arguments
23
c
24
integer ni, ic(7,MaxParticles)
25
double precision P5(0:4,MaxParticles),wgt
26
c
27
c Local
28
c
29
integer i,k,jj, jhel(MaxDecPart), idecay, idec(MaxParticles),j, im
30
integer JD(MaxDecPart),ICOL(2,MaxDecPart)
31
integer maxcol
32
double precision pcms(0:3), pd(0:3,MaxDecPart),pswgt
33
double precision p(0:3,MaxParticles)
34
double precision totwidth, dwgt, goal_wgt
35
double precision weight,r
36
logical done
37
real wgt_vegas
38
c
39
c External
40
c
41
real xran1
42
integer iseed
43
data iseed/1/ !no effect if already called
44
c
45
c Global
46
c
47
data s_unw,s_wei,s_ove/3*0d0/
48
data n_wei,n_ove/2*0/
49
c
50
include 'coupl.inc'
51
include 'calc_values.inc'
52
c
53
c first map the momenta into objects 0:3
54
c
55
do i=1,MaxParticles
56
do j=0,3
57
p(j,i)=p5(j,i)
58
enddo
59
enddo
60
61
c
62
c-- First find out how many particles there are in the
63
c event to decay
64
c
65
call find_particle(ip,ni,ic,k,idec)
66
if (k.eq.0) return !no particle to decay
67
c
68
c-- Loop over particles to decay
69
c
70
do jj=1,k
71
idecay=idec(jj)
72
c-- find the largest label present for the color flow
73
maxcol=max(ic(4,1),ic(5,1))
74
do i=2,ni
75
maxcol=max(maxcol,ic(4,i),ic(5,i))
76
enddo
77
cindex=maxcol+1 !available color index
78
c-- setup color of the decay particle
79
icol(1,1) = ic(4,idecay) ! color of particle
80
icol(2,1) = ic(5,idecay) ! anti-color of particle
81
c-- Boost to the momentum of particle to the CMS frame
82
do i=0,3
83
pcms(i) = -p(i,1) - p(i,2)
84
enddo
85
pcms(0)=-pcms(0)
86
call boostx(p(0,idecay),pcms,pd(0,1))
87
jhel(1) = ic(7,idecay) ! helicity of particle
88
c-- Unweigthed decay
89
done=.false.
90
do while (.not. done)
91
CALL GET_X(X,WGT_VEGAS)
92
CALL EVENT(PD(0,1),JHEL(1),JD(1),ICOL,WEIGHT)
93
weight=weight*real(wgt_vegas)
94
n_wei=n_wei+1
95
s_wei=s_wei+weight
96
r = xran1(iseed)
97
if(weight / r .gt. mxwgt) then
98
done=.true.
99
s_unw=s_unw+mxwgt
100
if (weight .gt. mxwgt) then
101
s_ove = s_ove-mxwgt+weight
102
n_ove = n_ove + 1
103
endif
104
endif
105
enddo
106
c-- Multiply by the branching ratio
107
wgt = wgt * bratio
108
c-- Boost back to lab frame and put momentum into P(0,x)
109
do i=1,3
110
pcms(i)=-pcms(i)
111
enddo
112
do i=2,MaxDecPart !loop over decay products
113
call boostx(pd(0,i),pcms,p(0,i-1+ni))
114
enddo
115
c-- Set the new particle information
116
c-- The info on the decayed particles are kept intact but
117
c-- status changed to decayed
118
do i=1,ND
119
ic(1,ni+i)=jd(i+1) !particle's IDs
120
ic(2,ni+i)=idecay !Mother info
121
ic(3,ni+i)=idecay !Mother info
122
ic(4,ni+i)=icol(1,i+1) ! color info
123
ic(5,ni+i)=icol(2,i+1) !anti-color info
124
ic(6,ni+i)=+1 !Final Particle
125
ic(7,ni+i)=jhel(1+i) !Helicity info
126
enddo
127
ni = ni + nd !Keep intermediate state decaying particle
128
ic(6,idecay)=+2 !This is decayed particle
129
c-- end loop over particles to decay
130
enddo
131
c
132
c finally map the P back to 0:4
133
c
134
do i=1,ni
135
do j=0,3
136
p5(j,i)=p(j,i)
137
enddo
138
c p5(4,i)=pdgmass(ic(1,i))
139
enddo
140
141
do i=ni-nd,ni
142
p5(4,i)=pdgmass(ic(1,i))
143
enddo
144
145
return
146
end
147
148
149
SUBROUTINE EVENT(PD,JHEL,JD,ICOL,DWGT)
150
c**********************************************************
151
c Calculation of the decay event quantities
152
c
153
c input: pd(0:3,1) decay particle momentum
154
c
155
c ouput: pd(0:3,MaxDecPart) momenta
156
c jhel(MaxDecPart) helicities
157
c jd(MaxDecPart) particle id's
158
c icol(2,MaxDecPart) color labels
159
c DWGT pswgt*emmesq*factor
160
c
161
c**********************************************************
162
implicit none
163
include 'decay.inc'
164
c
165
c Arguments
166
c
167
real*8 pd(0:3,MaxDecPart)
168
integer jhel(MaxDecPart),JD(MaxDecPart),ICOL(2,MaxDecPart)
169
real*8 dwgt
170
c
171
c Global
172
c
173
include 'coupl.inc'
174
include 'calc_values.inc'
175
c
176
c local
177
c
178
real*8 pswgt,emmesq,fudge,factor
179
real*8 aa,aa1,aa2,aa3
180
integer isign,jl,jvl,i,j
181
logical hadro_dec
182
integer multi_decay
183
real*8 xprob1,xprob2
184
c
185
c EXTERNAL
186
c
187
integer get_hel
188
real xran1
189
integer iseed
190
data iseed/1/ !no effect if already called
191
192
C------
193
C START
194
C------
195
196
DWGT =0d0
197
emmesq =0d0
198
isign =ip/abs(ip)
199
aa =0d0
200
aa1 =0d0
201
aa2 =0d0
202
do i=2,MaxDecPart
203
icol(1,i)=0
204
icol(2,i)=0
205
enddo
206
207
c write(*,*) 'from event: imode,ip',imode,ip
208
209
If(abs(ip).eq.6) then
210
*-----------------------------------------------------
211
* top decays
212
*-----------------------------------------------------
213
214
215
If(imode.eq.1) then
216
*------------------------
217
* t -> b w+
218
* t~ -> b~ w-
219
*------------------------
220
221
c-- masses
222
M1=tmass
223
M2=bmass
224
M3=wmass
225
c-- id's
226
jd(2)=isign*5 !b or b~
227
jd(3)=isign*24 !w+ or w-
228
c-- color
229
icol(1,2)=icol(1,1) ! color of the top
230
icol(2,2)=icol(2,1) !anti-color of the top
231
c-- couplings
232
GXX(1)=gwf(1)
233
GXX(2)=gwf(2)
234
c-- color*(bose factor)*number of helicities
235
factor=1d0*1d0*6d0
236
c-- helicities
237
jhel(2) = get_hel(xran1(iseed),2)
238
jhel(3) = get_hel(xran1(iseed),3)
239
c-- phase space
240
call phasespace(pd,pswgt)
241
if(pswgt.eq.0d0) return
242
c-- matrix element
243
if(isign.eq. 1) call emme_ffv(pd,jhel,emmesq)
244
if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)
245
c-- weight
246
dwgt=pswgt*emmesq*factor
247
c-- check that dwgt is a reasonable number
248
call check_nan(dwgt)
249
250
return
251
endif !end imode=1
252
253
if(imode.ge.2.and.imode.le.8) then
254
*------------------------
255
* t -> b vl l+ ,jj
256
* t~ -> b~ l vl~ ,jj
257
*------------------------
258
259
c-- masses
260
M1=tmass
261
M2=bmass
262
M3=0d0
263
M4=0d0
264
MV=WMASS
265
GV=WWIDTH
266
c-- color
267
icol(1,2)=icol(1,1) ! color of the top
268
icol(2,2)=icol(2,1) !anti-color of the top
269
c-- id's
270
if(isign.eq.1) then !t
271
jl=4
272
jvl=3
273
else !t~
274
jl=3
275
jvl=4
276
endif
277
jd(2)=isign*5 !b or b~
278
c-- rnd number
279
if(imode.ge.5) aa =dble(xran1(iseed))
280
if(imode.eq.8) aa1=dble(xran1(iseed))
281
c
282
c various cases imode=2,3,4,5,6,7,8
283
c
284
if(imode.eq.2) then
285
*-------------------
286
* t -> b ve e+
287
*-------------------
288
jd(jvl)=isign*12 !ve or ve~
289
jd(jl) =isign*(-11) !e+ or e-
290
291
elseif(imode.eq.3) then
292
*--------------------
293
* t -> b vm mu+
294
*--------------------
295
jd(jvl)=isign*14 !vm or vm~
296
jd(jl) =isign*(-13) !mu+ or mu-
297
298
elseif(imode.eq.4) then
299
*-------------------
300
* t -> b vt ta+
301
*-------------------
302
if(isign.eq.1) then !t
303
M3=0D0
304
M4=lmass
305
else !t~
306
M3=lmass
307
M4=0d0
308
endif
309
jd(jvl)=isign*16 !vt or vt~
310
jd(jl) =isign*(-15) !ta+ or ta-
311
312
elseif(imode.eq.5) then
313
*-------------------
314
* t -> b vl l+ (e+mu)
315
*-------------------
316
if(aa.lt.Half) then
317
jd(jvl)=isign*12 !ve or ve~
318
jd(jl) =isign*(-11) !e+ or e-
319
else
320
jd(jvl)=isign*14 !vm or vm~
321
jd(jl) =isign*(-13) !mu+ or mu-
322
endif
323
324
elseif(imode.eq.6) then
325
*-------------------
326
* t -> b vl l+ (e+mu+ta)
327
*-------------------
328
if(aa.lt.one/Three) then
329
jd(jvl)=isign*12 !ve or ve~
330
jd(jl) =isign*(-11) !e+ or e-
331
elseif(aa.lt.two/three) then
332
jd(jvl)=isign*14 !vm or vm~
333
jd(jl) =isign*(-13) !mu+ or mu-
334
else
335
if(isign.eq.1) then !t
336
M3=0D0
337
M4=lmass
338
else !t~
339
M3=lmass
340
M4=0d0
341
endif
342
jd(jvl)=isign*16 !vt or vt~
343
jd(jl) =isign*(-15) !ta+ or ta-
344
endif
345
346
elseif(imode.eq.7) then
347
*-------------------
348
* t -> b j j (ud,cs)
349
*-------------------
350
c-- color
351
icol(1,3)=cindex !position 3 is always a particle
352
icol(2,3)=0
353
icol(1,4)=0 !position 4 is always a anti-particle
354
icol(2,4)=cindex
355
356
if(aa.lt..5d0) then
357
jd(jvl)=isign*2 !u or u~
358
jd(jl) =isign*(-1) !d~ or d
359
else
360
if(isign.eq.1) then !t
361
M3=cmass
362
M4=0d0
363
else !t~
364
M3=0d0
365
M4=cmass
366
endif
367
jd(jvl)=isign*4 !c or c~
368
jd(jl) =isign*(-3) !s~ or s
369
endif
370
371
elseif(imode.eq.8) then
372
*-------------------
373
* t -> b anything
374
*-------------------
375
c
376
c first decide if W decays leptonically or hadronically
377
c
378
if(aa1.lt.3d0*br_w_lv) then !leptonic decay
379
hadro_dec=.false.
380
if(aa.lt.one/Three) then
381
jd(jvl)=isign*12 !ve or ve~
382
jd(jl) =isign*(-11) !e+ or e-
383
elseif(aa.lt.two/three) then
384
jd(jvl)=isign*14 !vm or vm~
385
jd(jl) =isign*(-13) !mu+ or mu-
386
else
387
if(isign.eq.1) then !t
388
M3=0D0
389
M4=lmass
390
else !t~
391
M3=lmass
392
M4=0d0
393
endif
394
jd(jvl)=isign*16 !vt or vt~
395
jd(jl) =isign*(-15) !ta+ or ta-
396
endif
397
398
else ! hadronic decay
399
hadro_dec=.true.
400
c-- color
401
icol(1,3)=cindex !position 3 is always a particle
402
icol(2,3)=0
403
icol(1,4)=0 !position 4 is always a anti-particle
404
icol(2,4)=cindex
405
406
if(aa.lt..5d0) then
407
jd(jvl)=isign*2 !u or u~
408
jd(jl) =isign*(-1) !d~ or d
409
else
410
if(isign.eq.1) then !t
411
M3=cmass
412
M4=0d0
413
else !t~
414
M3=0d0
415
M4=cmass
416
endif
417
jd(jvl)=isign*4 !c or c~
418
jd(jl) =isign*(-3) !s~ or s
419
endif
420
endif
421
422
endif !imode(from 2 to 8)
423
endif !imode
424
425
c-- couplings
426
GXX(1)=gwf(1)
427
GXX(2)=gwf(2)
428
c-- color*(bose factor)*number of helicities
429
factor=(3d0/3d0)*1d0*4d0
430
if(imode.eq.7) factor=factor*3d0 !quark color in jj
431
if(imode.eq.8.and.hadro_dec) factor=factor*3d0 !quark color in jj
432
c-- helicities
433
jhel(2) = get_hel(xran1(iseed),2)
434
jhel(3) = get_hel(xran1(iseed),2)
435
jhel(4) = -jhel(3)
436
c-- phase space
437
if(GV.eq.0d0) call error_trap("GV")
438
call phasespace(pd,pswgt)
439
if(pswgt.eq.0d0) return
440
c-- matrix element
441
if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq) !t
442
if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !t~
443
c-- weight
444
dwgt=pswgt*emmesq*factor
445
if(imode.eq.5) dwgt=dwgt*2d0 !l=e,mu
446
if(imode.eq.6) dwgt=dwgt*3d0 !l=e,mu,tau
447
if(imode.eq.7) dwgt=dwgt*2d0 !jj=ud,cs
448
if(imode.eq.8) then
449
if(hadro_dec) then
450
dwgt=dwgt*2d0/br_w_jj !jj=ud,cs
451
else
452
dwgt=dwgt*3d0/(3d0*br_w_lv) !l=e,mu,tau
453
endif
454
endif
455
c-- check that dwgt is a reasonable number
456
call check_nan(dwgt)
457
return
458
write(*,*) 'non-existing imode for top decays'
459
endif ! top decays
460
461
462
463
If(abs(ip).eq.15) then
464
*-----------------------------------------------------
465
* tau decays
466
*-----------------------------------------------------
467
468
469
if(imode.le.3) then
470
*------------------------
471
* ta -> vt vl l
472
*------------------------
473
474
c-- masses
475
M1=lmass
476
M2=0d0
477
M3=0d0
478
M4=0d0
479
MV=WMASS
480
GV=WWIDTH
481
c-- id's
482
if(isign.eq.1) then !tau-
483
jl=3
484
jvl=4
485
else !tau+
486
jl=4
487
jvl=3
488
endif
489
jd(2)=isign*16 !vt or vt~
490
c-- couplings
491
GXX(1)=gwf(1)
492
GXX(2)=gwf(2)
493
c-- rnd number
494
if(imode.eq.3) aa=dble(xran1(iseed))
495
496
if (imode.eq.1) then
497
*-------------------
498
* ta -> vt e ve
499
*-------------------
500
jd(jvl)=isign*(-12) !ve~ or ve
501
jd(jl) =isign*(11) !e- or e+
502
503
elseif(imode.eq.2) then
504
*-----------------------
505
* tau -> vt mu vm
506
*-----------------------
507
jd(jvl)=isign*(-14) !vm~ or vm
508
jd(jl) =isign*( 13) !mu- or mu+
509
510
elseif(imode.eq.3) then
511
*----------------------
512
* tau -> vt l vl (e+mu)
513
*----------------------
514
if(aa.lt.Half) then
515
jd(jvl)=isign*(-12) !ve~ or ve
516
jd(jl) =isign*(11) !e- or e+
517
else
518
jd(jvl)=isign*(-14) !vm~ or vm
519
jd(jl) =isign*( 13) !mu- or mu+
520
endif
521
522
endif
523
524
c-- color*(bose factor)*number of helicities
525
factor=1d0*1d0*1d0
526
c-- helicities: by definition of 3rd and 4th particle
527
jhel(2) = -isign
528
jhel(3) = -1
529
jhel(4) = 1
530
c-- phase space
531
call phasespace(pd,pswgt)
532
if(pswgt.eq.0d0) return
533
c-- matrix element
534
if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq) !tau-
535
if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !tau+
536
c-- weight
537
dwgt=pswgt*emmesq*factor
538
if(imode.eq.3) dwgt=dwgt*2d0 !l=e,mu
539
c-- check that dwgt is a reasonable number
540
call check_nan(dwgt)
541
return
542
543
544
elseif(imode.eq.4) then
545
*------------------------
546
* tau -> vt pi
547
*------------------------
548
549
c-- masses
550
M1=lmass
551
M2=0d0
552
M3=0.1395d0
553
c-- id's
554
jd(2)=isign*16 !vt or vt~
555
jd(3)=isign*(-211) !pi- pi+
556
c-- couplings
557
GXX(1)=gwf(1)
558
GXX(2)=gwf(2)
559
c-- color*(bose factor)*number of helicities
560
factor=1d0*1d0*2d0
561
c-- fudge to normalize
562
fudge=lwidth*br_ta_pi/0.00373
563
factor=factor*fudge
564
c-- helicities
565
jhel(2) = get_hel(xran1(iseed),2)
566
jhel(3) = 0 !not used
567
c-- phase space
568
call phasespace(pd,pswgt)
569
if(pswgt.eq.0d0) return
570
c-- matrix element
571
if(isign.eq. 1) call emme_ffs (pd,jhel,emmesq)
572
if(isign.eq.-1) call emme_fxfs(pd,jhel,emmesq)
573
c-- weight
574
dwgt=pswgt*emmesq*factor
575
c-- check that dwgt is a reasonable number
576
call check_nan(dwgt)
577
return
578
579
elseIf(imode.eq.5) then
580
*------------------------
581
* tau -> vt rho
582
*------------------------
583
584
c-- masses
585
M1=lmass
586
M2=0d0
587
M3=0.770d0
588
c-- id's
589
jd(2)=isign*16 !vt or vt~
590
jd(3)=isign*(-213) !rho- rho +
591
c-- couplings
592
GXX(1)=gwf(1)
593
GXX(2)=gwf(2)
594
c-- color*(bose factor)*number of helicities
595
factor=1d0*1d0*6d0
596
c-- fudge to normalize
597
fudge=lwidth*br_ta_ro/0.0182
598
factor=factor*fudge
599
c-- helicities
600
jhel(2) = get_hel(xran1(iseed),2)
601
jhel(3) = get_hel(xran1(iseed),3)
602
c-- phase space
603
call phasespace(pd,pswgt)
604
if(pswgt.eq.0d0) return
605
c-- matrix element
606
if(isign.eq. 1) call emme_ffv (pd,jhel,emmesq)
607
if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)
608
c-- weight
609
dwgt=pswgt*emmesq*factor
610
c-- check that dwgt is a reasonable number
611
call check_nan(dwgt)
612
return
613
endif !imodes
614
write(*,*) 'non-existing imode for tau decays'
615
endif !tau
616
617
618
619
If(ip.eq.23) then
620
*-----------------------------------------------------
621
* Z decays
622
*-----------------------------------------------------
623
624
c-- masses
625
M1=zmass
626
M2=0D0
627
M3=0D0
628
c-- couplings
629
GXX(1)=gzl(1)
630
GXX(2)=gzl(2)
631
c-- rnd number for flavour sums
632
aa =dble(xran1(iseed))
633
aa1=dble(xran1(iseed))
634
635
636
if (imode.eq.1) then
637
*-------------------
638
* z->e- e+
639
*-------------------
640
jd(2) = 11
641
jd(3) =-11
642
643
elseif(imode.eq.2) then
644
*--------------------
645
* z->mu- mu+
646
*--------------------
647
jd(2) = 13
648
jd(3) =-13
649
650
elseif(imode.eq.3) then
651
*-------------------
652
* z->ta- ta+
653
*-------------------
654
M2=lmass
655
jd(2) = 15
656
jd(3) =-15
657
658
elseif(imode.eq.4) then
659
*-------------------
660
* z->e- e+,mu-mu+
661
*-------------------
662
if(aa.lt.Half) then
663
jd(2) = 11
664
jd(3) =-11
665
else
666
jd(2) = 13
667
jd(3) =-13
668
endif
669
670
elseif(imode.eq.5) then
671
*-------------------
672
* z->e- e+,mu-mu+,ta-ta+
673
*-------------------
674
if(aa.lt.one/Three) then
675
jd(2) = 11
676
jd(3) =-11
677
elseif(aa.lt.two/three) then
678
jd(2) = 13
679
jd(3) =-13
680
else
681
M2=lmass
682
jd(2) = 15
683
jd(3) =-15
684
endif
685
686
elseif(imode.eq.6) then
687
*------------------------
688
* z->vl vl~
689
*------------------------
690
691
c-- couplings
692
GXX(1)=gzn(1)
693
GXX(2)=gzn(2)
694
695
if(aa.lt.one/Three) then
696
jd(2) = 12
697
jd(3) =-12
698
elseif(aa.lt.two/three) then
699
jd(2) = 14
700
jd(3) =-14
701
else
702
jd(2) = 16
703
jd(3) =-16
704
endif
705
706
707
elseif(imode.eq.7) then
708
*------------------------
709
* z-> b b~
710
*------------------------
711
c-- color
712
icol(1,2)=cindex !position 2 is always a particle
713
icol(2,2)=0
714
icol(1,3)=0 !position 3 is always a anti-particle
715
icol(2,3)=cindex
716
717
c-- couplings
718
GXX(1)=gzd(1)
719
GXX(2)=gzd(2)
720
M2=bmass
721
jd(2) = 5
722
jd(3) =-5
723
724
elseif(imode.eq.8) then
725
*------------------------
726
* z-> c c~
727
*------------------------
728
c-- color
729
icol(1,2)=cindex !position 2 is always a particle
730
icol(2,2)=0
731
icol(1,3)=0 !position 3 is always a anti-particle
732
icol(2,3)=cindex
733
734
c-- couplings
735
GXX(1)=gzu(1)
736
GXX(2)=gzu(2)
737
M2=cmass
738
jd(2) = 4
739
jd(3) =-4
740
741
742
elseif(imode.eq.9) then
743
*------------------------
744
* z-> u,d,c,s
745
*------------------------
746
c-- color
747
icol(1,2)=cindex !position 2 is always a particle
748
icol(2,2)=0
749
icol(1,3)=0 !position 3 is always a anti-particle
750
icol(2,3)=cindex
751
752
c-- probability of a uc or ds decay
753
754
xprob1=w_z_uu/(w_z_dd + w_z_uu)
755
756
if(aa1.lt.xprob1) then ! decay into ups
757
multi_decay=1
758
759
if(aa.lt. .5d0) then !u
760
M2=0d0
761
jd(2) = 2
762
jd(3) =-2
763
GXX(1)=gzu(1)
764
GXX(2)=gzu(2)
765
else !c
766
jd(2) = 4
767
jd(3) =-4
768
M2=cmass
769
GXX(1)=gzu(1)
770
GXX(2)=gzu(2)
771
endif
772
773
else !decay into downs
774
multi_decay=2
775
if(aa.lt..5d0) then !d
776
M2=0d0
777
jd(2) = 1
778
jd(3) =-1
779
GXX(1)=gzd(1)
780
GXX(2)=gzd(2)
781
else !s
782
jd(2) = 3
783
jd(3) =-3
784
M2=0d0
785
GXX(1)=gzd(1)
786
GXX(2)=gzd(2)
787
endif
788
789
endif !which decay: in ups or downs
790
791
elseif(imode.eq.10) then
792
*------------------------
793
* z-> u,d,c,s,b
794
*------------------------
795
c-- color
796
icol(1,2)=cindex !position 2 is always a particle
797
icol(2,2)=0
798
icol(1,3)=0 !position 3 is always a anti-particle
799
icol(2,3)=cindex
800
801
c-- probability of a uc or ds decay
802
803
xprob1=2d0*w_z_uu/(2d0*w_z_dd + 2d0*w_z_uu+w_z_bb)
804
xprob2=xprob1+2d0*w_z_dd/(2d0*w_z_dd + 2d0*w_z_uu+w_z_bb)
805
806
if(aa1.lt.xprob1) then ! decay into ups
807
multi_decay=1
808
if(aa.lt.0.5d0) then !u
809
jd(2) = 2
810
jd(3) =-2
811
M2=0d0
812
GXX(1)=gzu(1)
813
GXX(2)=gzu(2)
814
else !c
815
jd(2) = 4
816
jd(3) =-4
817
M2=cmass
818
GXX(1)=gzu(1)
819
GXX(2)=gzu(2)
820
endif
821
822
elseif(aa1.lt.xprob2) then ! decay into downs
823
multi_decay=2
824
if(aa.lt. .5d0) then !d
825
M2=0d0
826
jd(2) = 1
827
jd(3) =-1
828
GXX(1)=gzd(1)
829
GXX(2)=gzd(2)
830
else !s
831
jd(2) = 3
832
jd(3) =-3
833
M2=0d0
834
GXX(1)=gzd(1)
835
GXX(2)=gzd(2)
836
endif
837
else ! decay into b's
838
multi_decay=3
839
jd(2) = 5
840
jd(3) =-5
841
M2=bmass
842
GXX(1)=gzd(1)
843
GXX(2)=gzd(2)
844
endif
845
846
else
847
write(*,*) 'non-existing imode for Z decays'
848
endif ! imode
849
850
M3=M2
851
c-- color*(bose factor)*number of helicities
852
factor=1d0*1d0*4d0
853
if(imode.ge.7) factor=factor*3d0 !quark color in jj
854
c-- helicities
855
jhel(2) = get_hel(xran1(iseed),2)
856
jhel(3) = get_hel(xran1(iseed),2)
857
c-- phase space
858
call phasespace(pd,pswgt)
859
if(pswgt.eq.0d0) return
860
c-- matrix element
861
call emme_vff(pd,jhel,emmesq)
862
c-- weight
863
dwgt=pswgt*emmesq*factor
864
if(imode.eq.4) dwgt=dwgt*2d0 !l=e,mu
865
if(imode.eq.5) dwgt=dwgt*3d0 !l=e,mu,tau
866
if(imode.eq.6) dwgt=dwgt*3d0 !vl=ve,vm,vt
867
c
868
if(imode.eq.9) then !j=u,d,c,s
869
if(multi_decay.eq.1) then
870
dwgt=dwgt*2d0/xprob1 !j=u,c
871
else
872
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
873
endif
874
endif
875
c
876
if(imode.eq.10) then !j=u,d,c,s,b
877
if(multi_decay.eq.1) then
878
dwgt=dwgt*2d0/xprob1 !j=u,c
879
elseif(multi_decay.eq.2) then
880
dwgt=dwgt*2d0/(xprob2-xprob1) !j=d,s
881
elseif(multi_decay.eq.3) then
882
dwgt=dwgt/(1.-xprob2) !j=b
883
endif
884
endif
885
886
return
887
c-- check that dwgt is a reasonable number
888
call check_nan(dwgt)
889
endif ! Z
890
891
892
If(abs(ip).eq.24) then
893
*-----------------------------------------------------
894
* W decays
895
*-----------------------------------------------------
896
897
c-- masses
898
M1=wmass
899
M2=0d0
900
M3=0d0
901
c-- couplings
902
GXX(1)=gwf(1)
903
GXX(2)=gwf(2)
904
c-- id's
905
if(isign.eq.1) then !W+
906
jl=3
907
jvl=2
908
else !W-
909
jl=2
910
jvl=3
911
endif
912
c-- rnd number
913
if(imode.ge.4) aa =dble(xran1(iseed))
914
if(imode.eq.7) aa1=dble(xran1(iseed))
915
916
917
if(imode.eq.1) then
918
*------------------------
919
* w-> e ve
920
*------------------------
921
jd(jvl)=isign*12 !ve or ve~
922
jd(jl) =isign*(-11) !e+ or e-
923
924
elseif(imode.eq.2) then
925
*--------------------
926
* w -> mu vm
927
*--------------------
928
jd(jvl)=isign*14 !vm or vm~
929
jd(jl) =isign*(-13) !mu+ or mu-
930
931
elseif(imode.eq.3) then
932
*-------------------
933
* w -> ta vt
934
*-------------------
935
if(isign.eq.1) then !t
936
M2=0D0
937
M3=lmass
938
else !t~
939
M2=lmass
940
M3=0d0
941
endif
942
jd(jvl)=isign*16 !vt or vt~
943
jd(jl) =isign*(-15) !ta+ or ta-
944
945
elseif(imode.eq.4) then
946
*-------------------
947
* w -> vl l+ (e+mu)
948
*-------------------
949
if(aa.lt.Half) then
950
jd(jvl)=isign*12 !ve or ve~
951
jd(jl) =isign*(-11) !e+ or e-
952
else
953
jd(jvl)=isign*14 !vm or vm~
954
jd(jl) =isign*(-13) !mu+ or mu-
955
endif
956
957
elseif(imode.eq.5) then
958
*-------------------
959
* w -> vl l+ (e+mu+ta)
960
*-------------------
961
if(aa.lt.one/Three) then
962
jd(jvl)=isign*12 !ve or ve~
963
jd(jl) =isign*(-11) !e+ or e-
964
elseif(aa.lt.two/three) then
965
jd(jvl)=isign*14 !vm or vm~
966
jd(jl) =isign*(-13) !mu+ or mu-
967
else
968
if(isign.eq.1) then !t
969
M2=0D0
970
M3=lmass
971
else !t~
972
M2=lmass
973
M3=0d0
974
endif
975
jd(jvl)=isign*16 !vt or vt~
976
jd(jl) =isign*(-15) !ta+ or ta-
977
endif
978
979
elseif(imode.eq.6) then
980
*-------------------
981
* w -> j j (ud,cs)
982
*-------------------
983
c-- color
984
icol(1,2)=cindex !position 2 is always a particle
985
icol(2,2)=0
986
icol(1,3)=0 !position 3 is always a anti-particle
987
icol(2,3)=cindex
988
989
if(aa.lt..5d0) then
990
jd(jvl)=isign*2 !u or u~
991
jd(jl) =isign*(-1) !d~ or d
992
else
993
if(isign.eq.1) then !t
994
M2=cmass
995
M3=0d0
996
else !t~
997
M2=0d0
998
M3=cmass
999
endif
1000
jd(jvl)=isign*4 !c or c~
1001
jd(jl) =isign*(-3) !s~ or s
1002
endif
1003
1004
elseif(imode.eq.7) then
1005
*-------------------
1006
* w -> anything
1007
*-------------------
1008
c
1009
c first decide if W decays leptonically or hadronically
1010
c
1011
if(aa1.lt.3d0*br_w_lv) then !leptonic decay
1012
hadro_dec=.false.
1013
if(aa.lt.one/Three) then
1014
jd(jvl)=isign*12 !ve or ve~
1015
jd(jl) =isign*(-11) !e+ or e-
1016
elseif(aa.lt.two/three) then
1017
jd(jvl)=isign*14 !vm or vm~
1018
jd(jl) =isign*(-13) !mu+ or mu-
1019
else
1020
if(isign.eq.1) then !t
1021
M2=0D0
1022
M3=lmass
1023
else !t~
1024
M2=lmass
1025
M3=0d0
1026
endif
1027
jd(jvl)=isign*16 !vt or vt~
1028
jd(jl) =isign*(-15) !ta+ or ta-
1029
endif
1030
else
1031
hadro_dec=.true. !hadronic decay
1032
c-- color
1033
icol(1,2)=cindex !position 2 is always a particle
1034
icol(2,2)=0
1035
icol(1,3)=0 !position 3 is always a anti-particle
1036
icol(2,3)=cindex
1037
1038
if(aa.lt..5d0) then
1039
jd(jvl)=isign*2 !u or u~
1040
jd(jl) =isign*(-1) !d~ or d
1041
else
1042
if(isign.eq.1) then !t
1043
M2=cmass
1044
M3=0d0
1045
else !t~
1046
M2=0d0
1047
M3=cmass
1048
endif
1049
jd(jvl)=isign*4 !c or c~
1050
jd(jl) =isign*(-3) !s~ or s
1051
endif
1052
endif
1053
endif
1054
c
1055
c done with the modes: now I start the calculation
1056
c
1057
1058
c-- color*(bose factor)*number of helicities
1059
factor=1d0*1d0*2d0
1060
if(imode.eq.6) factor=factor*3d0 !quark color in jj
1061
if(imode.eq.7.and.hadro_dec) factor=factor*3d0 !quark color in jj
1062
c-- helicities
1063
jhel(2) = get_hel(xran1(iseed),2)
1064
jhel(3) = -jhel(2)
1065
c-- phase space
1066
call phasespace(pd,pswgt)
1067
if(pswgt.eq.0d0) return
1068
c-- matrix element
1069
call emme_vff(pd,jhel,emmesq)
1070
c-- weight
1071
dwgt=pswgt*emmesq*factor
1072
if(imode.eq.4) dwgt=dwgt*2d0 !l=e,mu
1073
if(imode.eq.5) dwgt=dwgt*3d0 !l=e,mu,tau
1074
if(imode.eq.6) dwgt=dwgt*2d0 !j=ud+cs
1075
if(imode.eq.7) then
1076
if(hadro_dec) then
1077
dwgt=dwgt*2d0/br_w_jj !jj=ud,cs
1078
else
1079
dwgt=dwgt*3d0/(3d0*br_w_lv) !l=e,mu,tau
1080
endif
1081
endif
1082
1083
c-- check that dwgt is a reasonable number
1084
call check_nan(dwgt)
1085
1086
return
1087
endif ! w
1088
1089
1090
If(ip.eq.25) then
1091
*-----------------------------------------------------
1092
* Higgs decays
1093
*-----------------------------------------------------
1094
1095
If(imode.eq.1) then
1096
*------------------------
1097
* h->b b~
1098
*------------------------
1099
c-- masses
1100
M1=hmass
1101
M2=bmass
1102
M3=bmass
1103
c-- id's
1104
jd(2) = 5
1105
jd(3) =-5
1106
c-- color
1107
icol(1,2)=cindex !position 2 is always a particle
1108
icol(2,2)=0
1109
icol(1,3)=0 !position 3 is always a anti-particle
1110
icol(2,3)=cindex
1111
c-- couplings
1112
GXX(1)=ghbot(1)
1113
GXX(2)=ghbot(2)
1114
c-- color*(bose factor)*number of helicities
1115
factor=3d0*1d0*2d0
1116
c-- helicities
1117
jhel(2) = get_hel(xran1(iseed),2)
1118
jhel(3) = jhel(2)
1119
c-- phase space
1120
call phasespace(pd,pswgt)
1121
if(pswgt.eq.0d0) return
1122
c-- matrix element
1123
call emme_hff(pd,jhel,emmesq)
1124
c-- weight
1125
dwgt=pswgt*emmesq*factor
1126
c-- check that dwgt is a reasonable number
1127
call check_nan(dwgt)
1128
return
1129
1130
elseif(imode.eq.2) then
1131
*------------------------
1132
* h->ta- ta+
1133
*------------------------
1134
c-- masses
1135
M1=hmass
1136
M2=lmass
1137
M3=lmass
1138
c-- id's
1139
jd(2) = 15
1140
jd(3) =-15
1141
c-- couplings
1142
GXX(1)=ghtau(1)
1143
GXX(2)=ghtau(2)
1144
c-- color*(bose factor)*number of helicities
1145
factor=1d0*1d0*2d0
1146
c-- helicities
1147
jhel(2) = get_hel(xran1(iseed),2)
1148
jhel(3) = jhel(2)
1149
c-- phase space
1150
call phasespace(pd,pswgt)
1151
if(pswgt.eq.0d0) return
1152
c-- matrix element
1153
call emme_hff(pd,jhel,emmesq)
1154
c-- weight
1155
dwgt=pswgt*emmesq*factor
1156
c-- check that dwgt is a reasonable number
1157
call check_nan(dwgt)
1158
1159
return
1160
1161
elseif(imode.eq.3) then
1162
*------------------------
1163
* h->mu- mu+
1164
*------------------------
1165
c-- masses
1166
M1=hmass
1167
M2=0d0
1168
M3=0d0
1169
c-- id's
1170
jd(2) = 13
1171
jd(3) =-13
1172
c-- couplings
1173
GXX(1)=ghtau(1)/lmass*0.105658389d0
1174
GXX(2)=ghtau(2)/lmass*0.105658389d0
1175
c-- color*(bose factor)*number of helicities
1176
factor=1d0*1d0*2d0
1177
c-- helicities
1178
jhel(2) = get_hel(xran1(iseed),2)
1179
jhel(3) = jhel(2)
1180
c-- phase space
1181
call phasespace(pd,pswgt)
1182
if(pswgt.eq.0d0) return
1183
c-- matrix element
1184
call emme_hff(pd,jhel,emmesq)
1185
c-- weight
1186
dwgt=pswgt*emmesq*factor
1187
c-- check that dwgt is a reasonable number
1188
call check_nan(dwgt)
1189
return
1190
1191
elseif(imode.eq.4) then
1192
*------------------------
1193
* h-> c c~
1194
*------------------------
1195
c-- masses
1196
M1=hmass
1197
M2=cmass
1198
M3=cmass
1199
c-- id's
1200
jd(2) = 4
1201
jd(3) =-4
1202
c-- color
1203
icol(1,2)=cindex !position 2 is always a particle
1204
icol(2,2)=0
1205
icol(1,3)=0 !position 3 is always a anti-particle
1206
icol(2,3)=cindex
1207
c-- couplings
1208
GXX(1)=ghcha(1)
1209
GXX(2)=ghcha(2)
1210
c-- color*(bose factor)*number of helicities
1211
factor=3d0*1d0*4d0
1212
c-- helicities
1213
jhel(2) = get_hel(xran1(iseed),2)
1214
jhel(3) = get_hel(xran1(iseed),2)
1215
c-- phase space
1216
call phasespace(pd,pswgt)
1217
if(pswgt.eq.0d0) return
1218
c-- matrix element
1219
call emme_hff(pd,jhel,emmesq)
1220
c-- weight
1221
dwgt=pswgt*emmesq*factor
1222
c-- check that dwgt is a reasonable number
1223
call check_nan(dwgt)
1224
return
1225
1226
elseif(imode.eq.5) then
1227
*------------------------
1228
* h-> t t~
1229
*------------------------
1230
c-- masses
1231
M1=hmass
1232
M2=tmass
1233
M3=tmass
1234
c-- id's
1235
jd(2) = 6
1236
jd(3) =-6
1237
c-- couplings
1238
GXX(1)=ghtop(1)
1239
GXX(2)=ghtop(2)
1240
c-- color*(bose factor)*number of helicities
1241
factor=3d0*1d0*4d0
1242
c-- helicities
1243
jhel(2) = get_hel(xran1(iseed),2)
1244
jhel(3) = get_hel(xran1(iseed),2)
1245
c-- phase space
1246
call phasespace(pd,pswgt)
1247
if(pswgt.eq.0d0) return
1248
c-- matrix element
1249
call emme_hff(pd,jhel,emmesq)
1250
c-- weight
1251
dwgt=pswgt*emmesq*factor
1252
c-- check that dwgt is a reasonable number
1253
call check_nan(dwgt)
1254
return
1255
1256
elseif(imode.eq.6) then
1257
*------------------------
1258
* h-> g g
1259
*------------------------
1260
c-- masses
1261
M1=hmass
1262
M2=0d0
1263
M3=0d0
1264
c-- id's
1265
jd(2) = 21
1266
jd(3) = 21
1267
c-- color
1268
icol(1,2)=cindex
1269
icol(2,2)=cindex+1
1270
icol(1,3)=cindex+1
1271
icol(2,3)=cindex
1272
c-- couplings
1273
GX=cmplx(1d0)
1274
c-- color*(bose factor)*number of helicities
1275
factor=8d0*.5d0*4d0
1276
c-- fudge factor for normalization
1277
factor=factor*hmass*2d0*pi*(SMBRG*SMWDTH)
1278
c-- helicities
1279
jhel(2) = get_hel(xran1(iseed),2)
1280
jhel(3) = get_hel(xran1(iseed),2)
1281
c-- phase space
1282
call phasespace(pd,pswgt)
1283
if(pswgt.eq.0d0) return
1284
c-- matrix element
1285
call emme_hvv(pd,jhel,emmesq)
1286
c-- weight
1287
dwgt=pswgt*emmesq*factor
1288
c-- check that dwgt is a reasonable number
1289
call check_nan(dwgt)
1290
return
1291
1292
elseif(imode.eq.7) then
1293
*------------------------
1294
* h-> a a
1295
*------------------------
1296
1297
c-- masses
1298
M1=hmass
1299
M2=0d0
1300
M3=0d0
1301
c-- id's
1302
jd(2) = 22
1303
jd(3) = 22
1304
c-- couplings
1305
GX=cmplx(1d0)
1306
c-- color*(bose factor)*number of helicities
1307
factor=1d0*.5d0*4d0
1308
c-- fudge factor for normalization
1309
factor=factor*hmass*16d0*pi*(SMBRGA*SMWDTH)
1310
c-- helicities
1311
jhel(2) = get_hel(xran1(iseed),2)
1312
jhel(3) = get_hel(xran1(iseed),2)
1313
c-- phase space
1314
call phasespace(pd,pswgt)
1315
if(pswgt.eq.0d0) return
1316
c-- matrix element
1317
call emme_hvv(pd,jhel,emmesq)
1318
c-- weight
1319
dwgt=pswgt*emmesq*factor
1320
c-- check that dwgt is a reasonable number
1321
call check_nan(dwgt)
1322
return
1323
1324
elseif(imode.eq.8) then
1325
*------------------------
1326
* h-> z a
1327
*------------------------
1328
c-- masses
1329
M1=hmass
1330
M2=zmass
1331
M3=0d0
1332
c-- id's
1333
jd(2) = 23
1334
jd(3) = 22
1335
c-- couplings
1336
GX=cmplx(1d0)
1337
c-- color*(bose factor)*number of helicities
1338
factor=1d0*1d0*6d0
1339
c-- fudge factor for normalization
1340
factor=factor*SMBRZGA*SMWDTH*
1341
&8d0*pi*hmass/(1d0-(zmass/hmass)**2)
1342
c-- helicities
1343
jhel(2) = get_hel(xran1(iseed),3) !-1,0,1
1344
jhel(3) = get_hel(xran1(iseed),2)
1345
c-- phase space
1346
call twomom (pd(0,1),pd(0,2),pd(0,3),pswgt)
1347
c-- phase space
1348
call phasespace(pd,pswgt)
1349
if(pswgt.eq.0d0) return
1350
c-- matrix element
1351
call emme_hvv(pd,jhel,emmesq)
1352
c-- weight
1353
dwgt=pswgt*emmesq*factor
1354
c-- check that dwgt is a reasonable number
1355
call check_nan(dwgt)
1356
return
1357
1358
elseif(imode.eq.9) then
1359
*------------------------
1360
* h-> w w
1361
*------------------------
1362
c-- masses
1363
M1=hmass
1364
M2=wmass
1365
M3=wmass
1366
c-- id's
1367
jd(2) = 24
1368
jd(3) =-24
1369
c-- couplings
1370
GX=gwwh
1371
c-- color*(bose factor)*number of helicities
1372
factor=1d0*1d0*9d0
1373
c-- helicities
1374
jhel(2) = get_hel(xran1(iseed),3)
1375
jhel(3) = get_hel(xran1(iseed),3)
1376
jhel(2) = INT(3e0*xran1(iseed)) - 1
1377
jhel(3) = INT(3e0*xran1(iseed)) - 1
1378
c-- phase space
1379
call phasespace(pd,pswgt)
1380
if(pswgt.eq.0d0) return
1381
c-- matrix element
1382
call emme_hvv(pd,jhel,emmesq)
1383
c-- weight
1384
dwgt=pswgt*emmesq*factor
1385
c-- check that dwgt is a reasonable number
1386
call check_nan(dwgt)
1387
return
1388
1389
1390
elseif(imode.eq.10) then
1391
*--------------------------------
1392
* h -> w* w -> l vl l' vl' (l,l'=e,mu)
1393
*--------------------------------
1394
1395
c-- masses
1396
M1=hmass
1397
M2=0d0
1398
M3=0d0
1399
M4=0d0
1400
M5=0d0
1401
MV=WMASS
1402
GV=WWIDTH
1403
c-- id's
1404
aa1=dble(xran1(iseed))
1405
aa2=dble(xran1(iseed))
1406
c-- W*-
1407
if(aa1.lt.half) then
1408
jd(2) = 11 !e-
1409
jd(3) = -12 !ve~
1410
else
1411
jd(2) = 13 !mu-
1412
jd(3) = -14 !vm~
1413
endif
1414
c-- W*+
1415
if(aa2.lt.half) then
1416
jd(4) = 12 !ve
1417
jd(5) = -11 !e+
1418
else
1419
jd(4) = 14 !vm
1420
jd(5) = -13 !mu+
1421
endif
1422
c-- couplings
1423
GX =gwwh
1424
GXX(1) =gwf(1)
1425
GXX(2) =gwf(2)
1426
GXX1(1)=gwf(1)
1427
GXX1(2)=gwf(2)
1428
c-- color*(bose factor)*number of helicities
1429
factor=1d0*1d0
1430
c-- helicities
1431
jhel(2)=-1
1432
jhel(3)=1
1433
jhel(4)=-1
1434
jhel(5)=1
1435
c-- phase space
1436
if(GV.eq.0d0) call error_trap("GV")
1437
call phasespace(pd,pswgt)
1438
if(pswgt.eq.0d0) return
1439
c-- matrix element
1440
call emme_h4f(pd,jhel,emmesq)
1441
c-- weight
1442
dwgt=pswgt*emmesq*factor
1443
c-- sum of flavours
1444
dwgt=dwgt*4d0
1445
c-- check that dwgt is a reasonable number
1446
call check_nan(dwgt)
1447
return
1448
1449
elseif(imode.eq.11) then
1450
*--------------------------------
1451
* h -> w* w -> l vl l' vl' (l,l'=e,mu,ta)
1452
*--------------------------------
1453
1454
c-- masses
1455
M1=hmass
1456
M2=0d0
1457
M3=0d0
1458
M4=0d0
1459
M5=0d0
1460
MV=WMASS
1461
GV=WWIDTH
1462
c-- id's
1463
aa1=dble(xran1(iseed))
1464
aa2=dble(xran1(iseed))
1465
c-- W*-
1466
if(aa1.lt.one/three) then
1467
jd(2) = 11 !e-
1468
jd(3) = -12 !ve~
1469
elseif(aa1.lt.two/three) then
1470
jd(2) = 13 !mu-
1471
jd(3) = -14 !vm~
1472
else
1473
M2 = lmass
1474
jd(2) = 15 !ta-
1475
jd(3) = -16 !vt~
1476
endif
1477
c-- W*+
1478
if(aa2.lt.one/three) then
1479
jd(4) = 12 !ve
1480
jd(5) = -11 !e+
1481
elseif(aa2.lt.two/three) then
1482
jd(4) = 14 !vm
1483
jd(5) = -13 !mu+
1484
else
1485
M5 = lmass
1486
jd(4) = 16 !vt
1487
jd(5) = -15 !ta+
1488
endif
1489
1490
c-- couplings
1491
GX =gwwh
1492
GXX(1) =gwf(1)
1493
GXX(2) =gwf(2)
1494
GXX1(1)=gwf(1)
1495
GXX1(2)=gwf(2)
1496
c-- color*(bose factor)*number of helicities
1497
factor=1d0*1d0
1498
c-- helicities
1499
jhel(2)=-1
1500
jhel(3)=1
1501
jhel(4)=-1
1502
jhel(5)=1
1503
c-- phase space
1504
if(GV.eq.0d0) call error_trap("GV")
1505
call phasespace(pd,pswgt)
1506
if(pswgt.eq.0d0) return
1507
c-- matrix element
1508
call emme_h4f(pd,jhel,emmesq)
1509
c-- weight
1510
dwgt=pswgt*emmesq*factor
1511
c-- sum of flavours
1512
dwgt=dwgt*9d0
1513
c-- check that dwgt is a reasonable number
1514
call check_nan(dwgt)
1515
return
1516
1517
elseif(imode.eq.12) then
1518
*--------------------------------
1519
* h -> w* w -> j j l vl (jj=ud,cs;l=e,mu)
1520
*--------------------------------
1521
1522
c-- masses
1523
M1=hmass
1524
M2=0d0
1525
M3=0d0
1526
M4=0d0
1527
M5=0d0
1528
MV=WMASS
1529
GV=WWIDTH
1530
c-- color
1531
icol(1,2)=cindex !position 2 is always a particle
1532
icol(2,2)=0
1533
icol(1,3)=0 !position 3 is always a anti-particle
1534
icol(2,3)=cindex
1535
c-- id's
1536
aa1=dble(xran1(iseed))
1537
aa2=dble(xran1(iseed))
1538
aa3=dble(xran1(iseed))
1539
1540
if(aa1.lt.half) then
1541
c-- W*-
1542
if(aa2.lt.half) then
1543
jd(2) = 1 !d
1544
jd(3) = -2 !u~
1545
else
1546
M3 = cmass
1547
jd(2) = 3 !s
1548
jd(3) = -4 !c~
1549
endif
1550
c-- W*+
1551
if(aa3.lt.half) then
1552
jd(4) = 12 !ve
1553
jd(5) = -11 !e+
1554
else
1555
jd(4) = 14 !vm
1556
jd(5) = -13 !mu+
1557
endif
1558
1559
else
1560
c-- W*+
1561
if(aa2.lt.half) then
1562
jd(2) = 2 !u
1563
jd(3) = -1 !d~
1564
else
1565
M2 = cmass
1566
jd(2) = 4 !c
1567
jd(3) = -3 !s~
1568
endif
1569
c-- W*-
1570
if(aa3.lt.half) then
1571
jd(4) = 11 !e-
1572
jd(5) = -12 !ve~
1573
else
1574
jd(4) = 13 !mu-
1575
jd(5) = -14 !vm~
1576
endif
1577
endif
1578
1579
c-- couplings
1580
GX =gwwh
1581
GXX(1) =gwf(1)
1582
GXX(2) =gwf(2)
1583
GXX1(1)=gwf(1)
1584
GXX1(2)=gwf(2)
1585
c-- color*(bose factor)*number of helicities
1586
factor=3d0*1d0
1587
c-- helicities
1588
jhel(2)=-1
1589
jhel(3)=1
1590
jhel(4)=-1
1591
jhel(5)=1
1592
c-- phase space
1593
if(GV.eq.0d0) call error_trap("GV")
1594
call phasespace(pd,pswgt)
1595
if(pswgt.eq.0d0) return
1596
c-- matrix element
1597
call emme_h4f(pd,jhel,emmesq)
1598
c-- weight
1599
dwgt=pswgt*emmesq*factor
1600
c-- sum of flavours
1601
dwgt=dwgt*4d0
1602
c-- sum of W+ and W- possibilities
1603
dwgt=dwgt*2d0
1604
c-- check that dwgt is a reasonable number
1605
call check_nan(dwgt)
1606
return
1607
1608
elseif(imode.eq.13) then
1609
*--------------------------------
1610
* h -> w* w -> j j l vl (jj=ud,cs;l=e,mu,ta)
1611
*--------------------------------
1612
1613
c-- masses
1614
M1=hmass
1615
M2=0d0
1616
M3=0d0
1617
M4=0d0
1618
M5=0d0
1619
MV=WMASS
1620
GV=WWIDTH
1621
c-- color
1622
icol(1,2)=cindex !position 2 is always a particle
1623
icol(2,2)=0
1624
icol(1,3)=0 !position 3 is always a anti-particle
1625
icol(2,3)=cindex
1626
c-- id's
1627
aa1=dble(xran1(iseed))
1628
aa2=dble(xran1(iseed))
1629
aa3=dble(xran1(iseed))
1630
1631
if(aa1.lt.half) then
1632
c-- W*-
1633
if(aa2.lt.half) then
1634
jd(2) = 1 !d
1635
jd(3) = -2 !u~
1636
else
1637
M3 = cmass
1638
jd(2) = 3 !s
1639
jd(3) = -4 !c~
1640
endif
1641
c-- W*+
1642
if(aa3.lt.one/three) then
1643
jd(4) = 12 !ve
1644
jd(5) = -11 !e+
1645
elseif(aa2.lt.two/three) then
1646
jd(4) = 14 !vm
1647
jd(5) = -13 !mu+
1648
else
1649
M5 = lmass
1650
jd(4) = 16 !vt
1651
jd(5) = -15 !ta+
1652
endif
1653
1654
else
1655
c-- W*+
1656
if(aa2.lt.half) then
1657
jd(2) = 2 !u
1658
jd(3) = -1 !d~
1659
else
1660
M2 = cmass
1661
jd(2) = 4 !c
1662
jd(3) = -3 !s~
1663
endif
1664
c-- W*-
1665
if(aa3.lt.one/three) then
1666
jd(4) = 11 !e-
1667
jd(5) = -12 !ve~
1668
elseif(aa3.lt.two/three) then
1669
jd(4) = 13 !mu-
1670
jd(5) = -14 !vm~
1671
else
1672
M4 = lmass
1673
jd(4) = 15 !ta-
1674
jd(5) = -16 !vt~
1675
endif
1676
endif
1677
1678
c-- couplings
1679
GX =gwwh
1680
GXX(1) =gwf(1)
1681
GXX(2) =gwf(2)
1682
GXX1(1)=gwf(1)
1683
GXX1(2)=gwf(2)
1684
c-- color*(bose factor)*number of helicities
1685
factor=3d0*1d0
1686
c-- helicities
1687
jhel(2)=-1
1688
jhel(3)=1
1689
jhel(4)=-1
1690
jhel(5)=1
1691
c-- phase space
1692
if(GV.eq.0d0) call error_trap("GV")
1693
call phasespace(pd,pswgt)
1694
if(pswgt.eq.0d0) return
1695
c-- matrix element
1696
call emme_h4f(pd,jhel,emmesq)
1697
c-- weight
1698
dwgt=pswgt*emmesq*factor
1699
c-- sum of flavours
1700
dwgt=dwgt*6d0
1701
c-- sum of W+ and W- possibilities
1702
dwgt=dwgt*2d0
1703
c-- check that dwgt is a reasonable number
1704
call check_nan(dwgt)
1705
return
1706
1707
elseif(imode.eq.14) then
1708
*------------------------
1709
* h-> z z
1710
*------------------------
1711
c-- masses
1712
M1=hmass
1713
M2=zmass
1714
M3=zmass
1715
c-- id's
1716
jd(2) = 23
1717
jd(3) = 23
1718
c-- couplings
1719
GX=gzzh
1720
c-- color*(bose factor)*number of helicities
1721
factor=1d0*.5d0*9d0
1722
c-- helicities
1723
jhel(2) = get_hel(xran1(iseed),3)
1724
jhel(3) = get_hel(xran1(iseed),3)
1725
c-- phase space
1726
call phasespace(pd,pswgt)
1727
if(pswgt.eq.0d0) return
1728
c-- matrix element
1729
call emme_hvv(pd,jhel,emmesq)
1730
c-- weight
1731
dwgt=pswgt*emmesq*factor
1732
c-- check that dwgt is a reasonable number
1733
call check_nan(dwgt)
1734
return
1735
1736
elseif(imode.eq.15) then
1737
*---------------------------------
1738
* h -> z* z -> l- l+ l-' l+'(l,l'=e,mu)
1739
*---------------------------------
1740
1741
c-- masses
1742
M1=hmass
1743
M2=0d0
1744
M3=0d0
1745
M4=0d0
1746
M5=0d0
1747
MV=ZMASS
1748
GV=ZWIDTH
1749
c-- id's
1750
aa1=dble(xran1(iseed))
1751
aa2=dble(xran1(iseed))
1752
c Z1
1753
if(aa1.lt.half) then
1754
jd(2) = 11 !e-
1755
jd(3) = -11 !e+
1756
else
1757
jd(2) = 13 !mu-
1758
jd(3) = -13 !mu+
1759
endif
1760
c Z2
1761
if(aa2.lt.half) then
1762
jd(4) = 11 !e-
1763
jd(5) = -11 !e+
1764
else
1765
jd(4) = 13 !mu-
1766
jd(5) = -13 !mu+
1767
endif
1768
c-- couplings
1769
GX =gzzh
1770
GXX(1)=gzl(1)
1771
GXX(2)=gzl(2)
1772
GXX1(1)=gzl(1)
1773
GXX1(2)=gzl(2)
1774
c-- color*(bose factor)*number of helicities
1775
factor=1d0*1d0*4d0/2d0
1776
c-- helicities
1777
jhel(2)=+1
1778
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1779
jhel(3)=-jhel(2)
1780
jhel(4)=+1
1781
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1782
jhel(5)=-jhel(4)
1783
c-- phase space
1784
if(GV.eq.0d0) call error_trap("GV")
1785
call phasespace(pd,pswgt)
1786
if(pswgt.eq.0d0) return
1787
c-- matrix element
1788
call emme_h4f(pd,jhel,emmesq)
1789
c-- weight
1790
dwgt=pswgt*emmesq*factor
1791
c-- sum of flavours
1792
dwgt=dwgt*4d0
1793
c-- check that dwgt is a reasonable number
1794
call check_nan(dwgt)
1795
return
1796
1797
1798
elseif(imode.eq.16) then
1799
*---------------------------------
1800
* h -> z* z -> l- l+ l-' l+'(l,l'=e,mu,ta)
1801
*---------------------------------
1802
1803
c-- masses
1804
M1=hmass
1805
M2=0d0
1806
M3=0d0
1807
M4=0d0
1808
M5=0d0
1809
MV=ZMASS
1810
GV=ZWIDTH
1811
c-- id's
1812
aa1=dble(xran1(iseed))
1813
aa2=dble(xran1(iseed))
1814
c Z
1815
if(aa1.lt.one/three) then
1816
jd(2) = 11 !e-
1817
jd(3) = -11 !e+
1818
elseif(aa1.lt.two/three) then
1819
jd(2) = 13 !mu-
1820
jd(3) = -13 !mu+
1821
else
1822
M2 = lmass
1823
M3 = lmass
1824
jd(2) = 15 !ta-
1825
jd(3) = -15 !ta+
1826
endif
1827
c Z
1828
if(aa2.lt.one/three) then
1829
jd(4) = 11 !e-
1830
jd(5) = -11 !e+
1831
elseif(aa2.lt.two/three) then
1832
jd(4) = 13 !mu-
1833
jd(5) = -13 !mu+
1834
else
1835
M4 = lmass
1836
M5 = lmass
1837
jd(4) = 15 !ta-
1838
jd(5) = -15 !ta+
1839
endif
1840
c-- couplings
1841
GX =gzzh
1842
GXX(1)=gzl(1)
1843
GXX(2)=gzl(2)
1844
GXX1(1)=gzl(1)
1845
GXX1(2)=gzl(2)
1846
c-- color*(bose factor)*number of helicities
1847
factor=1d0*1d0*4d0/2d0
1848
c-- helicities
1849
jhel(2)=+1
1850
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1851
jhel(3)=-jhel(2)
1852
jhel(4)=+1
1853
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1854
jhel(5)=-jhel(4)
1855
c-- phase space
1856
if(GV.eq.0d0) call error_trap("GV")
1857
call phasespace(pd,pswgt)
1858
if(pswgt.eq.0d0) return
1859
c-- matrix element
1860
call emme_h4f(pd,jhel,emmesq)
1861
c-- weight
1862
dwgt=pswgt*emmesq*factor
1863
c-- sum of flavours
1864
dwgt=dwgt*9d0
1865
c-- check that dwgt is a reasonable number
1866
call check_nan(dwgt)
1867
return
1868
1869
elseif(imode.eq.17) then
1870
*---------------------------------
1871
* h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu )
1872
*---------------------------------
1873
1874
c-- masses
1875
M1=hmass
1876
M2=0d0
1877
M3=0d0
1878
M4=0d0
1879
M5=0d0
1880
MV=ZMASS
1881
GV=ZWIDTH
1882
c-- color
1883
icol(1,2)=cindex !position 2 is always a particle
1884
icol(2,2)=0
1885
icol(1,3)=0 !position 3 is always a anti-particle
1886
icol(2,3)=cindex
1887
c-- couplings
1888
GX =gzzh
1889
GXX(1)=gzl(1)
1890
GXX(2)=gzl(2)
1891
GXX1(1)=gzl(1)
1892
GXX1(2)=gzl(2)
1893
c-- id's
1894
aa =dble(xran1(iseed))
1895
aa1=dble(xran1(iseed))
1896
aa2=dble(xran1(iseed))
1897
c Z1
1898
c-- probability of a uc or ds decay
1899
1900
xprob1=w_z_uu/(w_z_dd + w_z_uu)
1901
1902
if(aa1.lt.xprob1) then ! decay into ups
1903
multi_decay=1
1904
1905
if(aa.lt. .5d0) then !u
1906
M2=0d0
1907
jd(2) = 2
1908
jd(3) =-2
1909
GXX(1)=gzu(1)
1910
GXX(2)=gzu(2)
1911
else !c
1912
jd(2) = 4
1913
jd(3) =-4
1914
M2=cmass
1915
GXX(1)=gzu(1)
1916
GXX(2)=gzu(2)
1917
endif
1918
1919
else !decay into downs
1920
multi_decay=2
1921
if(aa.lt..5d0) then !d
1922
M2=0d0
1923
jd(2) = 1
1924
jd(3) =-1
1925
GXX(1)=gzd(1)
1926
GXX(2)=gzd(2)
1927
else !s
1928
jd(2) = 3
1929
jd(3) =-3
1930
M2=0d0
1931
GXX(1)=gzd(1)
1932
GXX(2)=gzd(2)
1933
endif
1934
1935
endif !which decay: in ups or downs
1936
M3=M2
1937
c Z2
1938
if(aa2.lt.half) then
1939
jd(4) = 11 !e-
1940
jd(5) = -11 !e+
1941
else
1942
jd(4) = 13 !mu-
1943
jd(5) = -13 !mu+
1944
endif
1945
1946
c-- color*(bose factor)*number of helicities
1947
factor=3d0*1d0*4d0/2d0
1948
c-- helicities
1949
jhel(2)=+1
1950
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1951
jhel(3)=-jhel(2)
1952
jhel(4)=+1
1953
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1954
jhel(5)=-jhel(4)
1955
c-- phase space
1956
if(GV.eq.0d0) call error_trap("GV")
1957
call phasespace(pd,pswgt)
1958
if(pswgt.eq.0d0) return
1959
c-- matrix element
1960
call emme_h4f(pd,jhel,emmesq)
1961
c-- weight
1962
dwgt=pswgt*emmesq*factor
1963
c-- sum of flavours
1964
dwgt=dwgt*2d0
1965
if(multi_decay.eq.1) then
1966
dwgt=dwgt*2d0/xprob1 !j=u,c
1967
else
1968
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
1969
endif
1970
c-- sum of two possibilities for the decay of a Z
1971
dwgt=dwgt*2d0
1972
c-- check that dwgt is a reasonable number
1973
call check_nan(dwgt)
1974
return
1975
1976
elseif(imode.eq.18) then
1977
*---------------------------------
1978
* h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu,ta)
1979
*---------------------------------
1980
1981
c-- masses
1982
M1=hmass
1983
M2=0d0
1984
M3=0d0
1985
M4=0d0
1986
M5=0d0
1987
MV=ZMASS
1988
GV=ZWIDTH
1989
c-- color
1990
icol(1,2)=cindex !position 2 is always a particle
1991
icol(2,2)=0
1992
icol(1,3)=0 !position 3 is always a anti-particle
1993
icol(2,3)=cindex
1994
c-- couplings
1995
GX =gzzh
1996
GXX(1)=gzl(1)
1997
GXX(2)=gzl(2)
1998
GXX1(1)=gzl(1)
1999
GXX1(2)=gzl(2)
2000
c-- id's
2001
aa =dble(xran1(iseed))
2002
aa1=dble(xran1(iseed))
2003
aa2=dble(xran1(iseed))
2004
c Z1
2005
xprob1=w_z_uu/(w_z_dd + w_z_uu)
2006
if(aa1.lt.xprob1) then ! decay into ups
2007
multi_decay=1
2008
2009
if(aa.lt. .5d0) then !u
2010
M2=0d0
2011
jd(2) = 2
2012
jd(3) =-2
2013
GXX(1)=gzu(1)
2014
GXX(2)=gzu(2)
2015
else !c
2016
jd(2) = 4
2017
jd(3) =-4
2018
M2=cmass
2019
GXX(1)=gzu(1)
2020
GXX(2)=gzu(2)
2021
endif
2022
2023
else !decay into downs
2024
multi_decay=2
2025
if(aa.lt..5d0) then !d
2026
M2=0d0
2027
jd(2) = 1
2028
jd(3) =-1
2029
GXX(1)=gzd(1)
2030
GXX(2)=gzd(2)
2031
else !s
2032
jd(2) = 3
2033
jd(3) =-3
2034
M2=0d0
2035
GXX(1)=gzd(1)
2036
GXX(2)=gzd(2)
2037
endif
2038
2039
endif !which decay: in ups or downs
2040
M3=M2
2041
c Z2
2042
if(aa2.lt.one/three) then
2043
jd(4) = 11 !e-
2044
jd(5) = -11 !e+
2045
elseif(aa2.lt.two/three) then
2046
jd(4) = 13 !mu-
2047
jd(5) = -13 !mu+
2048
else
2049
M4 = lmass
2050
M5 = lmass
2051
jd(4) = 15 !ta-
2052
jd(5) = -15 !ta+
2053
endif
2054
2055
2056
c-- color*(bose factor)*number of helicities
2057
factor=3d0*1d0*4d0/2d0
2058
c-- helicities
2059
jhel(2)=+1
2060
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2061
jhel(3)=-jhel(2)
2062
jhel(4)=+1
2063
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2064
jhel(5)=-jhel(4)
2065
c-- phase space
2066
if(GV.eq.0d0) call error_trap("GV")
2067
call phasespace(pd,pswgt)
2068
if(pswgt.eq.0d0) return
2069
c-- matrix element
2070
call emme_h4f(pd,jhel,emmesq)
2071
c-- weight
2072
dwgt=pswgt*emmesq*factor
2073
c-- sum of flavours
2074
dwgt=dwgt*3d0
2075
if(multi_decay.eq.1) then
2076
dwgt=dwgt*2d0/xprob1 !j=u,c
2077
else
2078
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
2079
endif
2080
c-- sum of two possibilities for the decay of a Z
2081
dwgt=dwgt*2d0
2082
c-- check that dwgt is a reasonable number
2083
call check_nan(dwgt)
2084
return
2085
2086
elseif(imode.eq.19) then
2087
*---------------------------------
2088
* h -> z* z -> b b~ l- l+ (l=e,mu)
2089
*---------------------------------
2090
2091
c-- masses
2092
M1=hmass
2093
M2=bmass
2094
M3=bmass
2095
M4=0d0
2096
M5=0d0
2097
MV=ZMASS
2098
GV=ZWIDTH
2099
c-- color
2100
icol(1,2)=cindex !position 2 is always a particle
2101
icol(2,2)=0
2102
icol(1,3)=0 !position 3 is always a anti-particle
2103
icol(2,3)=cindex
2104
c-- couplings
2105
GX =gzzh
2106
GXX(1) =gzd(1)
2107
GXX(2) =gzd(2)
2108
GXX1(1)=gzl(1)
2109
GXX1(2)=gzl(2)
2110
c-- id's
2111
jd(2) = 5
2112
jd(3) =-5
2113
aa=dble(xran1(iseed))
2114
c Z2
2115
if(aa.lt.half) then
2116
jd(4) = 11 !e-
2117
jd(5) = -11 !e+
2118
else
2119
jd(4) = 13 !mu-
2120
jd(5) = -13 !mu+
2121
endif
2122
2123
c-- color*(bose factor)*number of helicities
2124
factor=3d0*1d0*8d0/2d0
2125
c-- helicities
2126
jhel(2) = get_hel(xran1(iseed),2)
2127
jhel(3) = get_hel(xran1(iseed),2)
2128
jhel(4) = get_hel(xran1(iseed),2)
2129
jhel(5) = -jhel(4)
2130
c-- phase space
2131
if(GV.eq.0d0) call error_trap("GV")
2132
call phasespace(pd,pswgt)
2133
if(pswgt.eq.0d0) return
2134
c-- matrix element
2135
call emme_h4f(pd,jhel,emmesq)
2136
c-- weight
2137
dwgt=pswgt*emmesq*factor
2138
c-- sum of flavours
2139
dwgt=dwgt*2d0
2140
c-- sum of two possibilities for the decay of a Z
2141
dwgt=dwgt*2d0
2142
c-- check that dwgt is a reasonable number
2143
call check_nan(dwgt)
2144
return
2145
2146
elseif(imode.eq.20) then
2147
*---------------------------------
2148
* h -> z* z -> b b~ l- l+ (l=e,mu,ta )
2149
*---------------------------------
2150
2151
c-- masses
2152
M1=hmass
2153
M2=bmass
2154
M3=bmass
2155
M4=0d0
2156
M5=0d0
2157
MV=ZMASS
2158
GV=ZWIDTH
2159
c-- color
2160
icol(1,2)=cindex !position 2 is always a particle
2161
icol(2,2)=0
2162
icol(1,3)=0 !position 3 is always a anti-particle
2163
icol(2,3)=cindex
2164
c-- couplings
2165
GX =gzzh
2166
GXX(1) =gzd(1)
2167
GXX(2) =gzd(2)
2168
GXX1(1)=gzl(1)
2169
GXX1(2)=gzl(2)
2170
c-- color*(bose factor)*number of helicities
2171
factor=3d0*1d0*16d0
2172
c-- id's
2173
jd(2) = 5
2174
jd(3) =-5
2175
aa=dble(xran1(iseed))
2176
if(aa.lt.one/three) then
2177
jd(4) = 11 !e-
2178
jd(5) = -11 !e+
2179
elseif(aa.lt.two/three) then
2180
jd(4) = 13 !mu-
2181
jd(5) = -13 !mu+
2182
else
2183
M4 = lmass
2184
M5 = lmass
2185
jd(4) = 15 !ta-
2186
jd(5) = -15 !ta+
2187
endif
2188
2189
c-- color*(bose factor)*number of helicities
2190
factor=3d0*1d0*8d0/2d0
2191
c-- helicities
2192
jhel(2) = get_hel(xran1(iseed),2)
2193
jhel(3) = get_hel(xran1(iseed),2)
2194
jhel(4) = get_hel(xran1(iseed),2)
2195
jhel(5) = -jhel(4)
2196
c-- phase space
2197
if(GV.eq.0d0) call error_trap("GV")
2198
call phasespace(pd,pswgt)
2199
if(pswgt.eq.0d0) return
2200
c-- matrix element
2201
call emme_h4f(pd,jhel,emmesq)
2202
c-- weight
2203
dwgt=pswgt*emmesq*factor
2204
c-- sum of flavours
2205
dwgt=dwgt*3d0
2206
c-- sum of two possibilities for the decay of a Z
2207
dwgt=dwgt*2d0
2208
c-- check that dwgt is a reasonable number
2209
call check_nan(dwgt)
2210
return
2211
2212
elseif(imode.eq.21) then
2213
*---------------------------------
2214
* h -> z* z -> vl vl~ l-' l+'(l=e,mu,ta,l'=e,mu)
2215
*---------------------------------
2216
2217
c-- masses
2218
M1=hmass
2219
M2=0d0
2220
M3=0d0
2221
M4=0d0
2222
M5=0d0
2223
MV=ZMASS
2224
GV=ZWIDTH
2225
c-- couplings
2226
GX =gzzh
2227
GXX(1)=gzn(1)
2228
GXX(2)=gzn(2)
2229
GXX1(1)=gzl(1)
2230
GXX1(2)=gzl(2)
2231
c-- id's
2232
aa1=dble(xran1(iseed))
2233
aa2=dble(xran1(iseed))
2234
c Z
2235
if(aa1.lt.one/three) then
2236
jd(2) = 12 !ve
2237
jd(3) = -12 !ve~
2238
elseif(aa1.lt.two/three) then
2239
jd(2) = 14 !vm
2240
jd(3) = -14 !vm~
2241
else
2242
jd(2) = 16 !vt
2243
jd(3) = -16 !vt~
2244
endif
2245
c Z2
2246
if(aa2.lt.half) then
2247
jd(4) = 11 !e-
2248
jd(5) = -11 !e+
2249
else
2250
jd(4) = 13 !mu-
2251
jd(5) = -13 !mu+
2252
endif
2253
2254
c-- color*(bose factor)*number of helicities
2255
factor=1d0*1d0*4d0/2d0
2256
c-- helicities
2257
jhel(2)=+1
2258
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2259
jhel(3)=-jhel(2)
2260
jhel(4)=+1
2261
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2262
jhel(5)=-jhel(4)
2263
c-- phase space
2264
if(GV.eq.0d0) call error_trap("GV")
2265
call phasespace(pd,pswgt)
2266
if(pswgt.eq.0d0) return
2267
c-- matrix element
2268
call emme_h4f(pd,jhel,emmesq)
2269
c-- weight
2270
dwgt=pswgt*emmesq*factor
2271
c-- sum of flavours
2272
dwgt=dwgt*6d0
2273
c-- sum of two possibilities for the decay of a Z
2274
dwgt=dwgt*2d0
2275
c-- check that dwgt is a reasonable number
2276
call check_nan(dwgt)
2277
return
2278
2279
elseif(imode.eq.22) then
2280
*---------------------------------
2281
* h -> z* z -> vl vl~ l-' l+'(l,l'=e,mu,ta)
2282
*---------------------------------
2283
2284
c-- masses
2285
M1=hmass
2286
M2=0d0
2287
M3=0d0
2288
M4=0d0
2289
M5=0d0
2290
MV=ZMASS
2291
GV=ZWIDTH
2292
c-- couplings
2293
GX =gzzh
2294
GXX(1) =gzn(1)
2295
GXX(2) =gzn(2)
2296
GXX1(1)=gzl(1)
2297
GXX1(2)=gzl(2)
2298
c-- color*(bose factor)*number of helicities
2299
factor=1d0*1d0*16d0
2300
c-- id's
2301
aa1=dble(xran1(iseed))
2302
aa2=dble(xran1(iseed))
2303
c Z
2304
if(aa1.lt.one/three) then
2305
jd(2) = 12 !ve
2306
jd(3) = -12 !ve~
2307
elseif(aa1.lt.two/three) then
2308
jd(2) = 14 !vm
2309
jd(3) = -14 !vm~
2310
else
2311
jd(2) = 16 !vt
2312
jd(3) = -16 !vt~
2313
endif
2314
c Z
2315
if(aa2.lt.one/three) then
2316
jd(4) = 11 !e-
2317
jd(5) = -11 !e+
2318
elseif(aa2.lt.two/three) then
2319
jd(4) = 13 !mu-
2320
jd(5) = -13 !mu+
2321
else
2322
M4 = lmass
2323
M5 = lmass
2324
jd(4) = 15 !ta-
2325
jd(5) = -15 !ta+
2326
endif
2327
2328
2329
c-- color*(bose factor)*number of helicities
2330
factor=1d0*1d0*4d0/2d0
2331
c-- helicities
2332
jhel(2)=+1
2333
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2334
jhel(3)=-jhel(2)
2335
jhel(4)=+1
2336
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2337
jhel(5)=-jhel(4)
2338
c-- phase space
2339
if(GV.eq.0d0) call error_trap("GV")
2340
call phasespace(pd,pswgt)
2341
if(pswgt.eq.0d0) return
2342
c-- matrix element
2343
call emme_h4f(pd,jhel,emmesq)
2344
c-- weight
2345
dwgt=pswgt*emmesq*factor
2346
c-- sum of flavours
2347
dwgt=dwgt*9d0
2348
c-- sum of two possibilities for the decay of a Z
2349
dwgt=dwgt*2d0
2350
c-- check that dwgt is a reasonable number
2351
call check_nan(dwgt)
2352
return
2353
2354
endif !imode
2355
2356
endif !higgs
2357
2358
write(*,*) 'from decay:end of decay reached'
2359
2360
return
2361
end
2362
2363
2364
2365
2366
INTEGER FUNCTION GET_HEL(RND,ID)
2367
**********************************************************
2368
* GIVEN THE RND NUMBER, RETURNS
2369
* +1, -1 WHEN ID=2
2370
* +1,0,-1 WHEN ID=3
2371
**********************************************************
2372
implicit none
2373
c
2374
c arguments
2375
c
2376
real*4 rnd
2377
integer id
2378
c----------
2379
c Begin
2380
c----------
2381
2382
if(id.eq.2) then
2383
get_hel=-1
2384
if(rnd.gt.0.5e0) get_hel=1
2385
elseif(id.eq.3) then
2386
if(rnd.lt.0.333333e0)then
2387
get_hel=-1
2388
elseif(rnd.lt.0.666666e0) then
2389
get_hel=0
2390
else
2391
get_hel=1
2392
endif
2393
else
2394
write(*,*) 'choice not implemented in gen_hel'
2395
endif
2396
return
2397
end
2398
2399
2400
SUBROUTINE GET_X(X,WGT)
2401
C-------------------------------------------------------
2402
C PROVIDES THE X(NDIM) AND THE WGT OBTAINED FROM
2403
C A PREVIOUS RUN OF VEGAS
2404
C-------------------------------------------------------
2405
C
2406
IMPLICIT NONE
2407
C
2408
C PARAMETERS
2409
C
2410
INTEGER NDMX,MXDIM
2411
PARAMETER (NDMX=50,MXDIM=10)
2412
C
2413
C ARGUMENTS
2414
C
2415
REAL*4 X(MXDIM),WGT
2416
C
2417
C LOCAL
2418
C
2419
INTEGER i,j,k,ia(MXDIM)
2420
REAL rc,xo,xn,rndn
2421
C
2422
C GLOBAL
2423
C
2424
INTEGER NG,NDIM
2425
REAL region(2*MXDIM),xi(NDMX,MXDIM),xnd,dx(MXDIM)
2426
COMMON /VEGAS_PAR1/NG,NDIM
2427
COMMON /VEGAS_PAR2/region,xi,xnd,dx
2428
C
2429
C EXTERNAL
2430
C
2431
REAL*4 XRAN1
2432
EXTERNAL XRAN1
2433
INTEGER IDUM
2434
DATA idum/1/ ! DOES NOT AFFECT PREVIOUS SETTING
2435
c
2436
c-- start the loop over dimensions
2437
c
2438
wgt=1.
2439
do 17 j=1,ndim
2440
c fax: avoid random numbers exactly equal to 0. or 1.
2441
303 rndn=xran1(idum)
2442
if(rndn.eq.1e0.or.rndn.eq.0e0) goto 303
2443
xn=rndn*xnd+1.
2444
ia(j)=max(min(int(xn),NDMX),1)
2445
if(ia(j).gt.1)then
2446
xo=xi(ia(j),j)-xi(ia(j)-1,j)
2447
rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
2448
else
2449
xo=xi(ia(j),j)
2450
rc=(xn-ia(j))*xo
2451
endif
2452
x(j)=region(j)+rc*dx(j)
2453
wgt=wgt*xo*xnd
2454
17 continue
2455
2456
RETURN
2457
END
2458
2459
2460
2461
2462
SUBROUTINE ERROR_TRAP(STRING)
2463
C-------------------------------------------------------
2464
C RETURNS INFORMATION ABOUT SOME ERROR THAT OCCURED
2465
C-------------------------------------------------------
2466
C
2467
IMPLICIT NONE
2468
C
2469
C ARGUMENTS
2470
C
2471
CHARACTER*20 STRING
2472
C
2473
C BEGIN
2474
C
2475
IF(STRING(1:2).EQ."GV") THEN
2476
write(*,*) '*****************************************************'
2477
write(*,*) '* *'
2478
write(*,*) '* >>>>ERROR TRAP CALLED<<<< *'
2479
write(*,*) '* *'
2480
write(*,*) '* the width of the vector boson(s) entering *'
2481
write(*,*) '* the selected decay should be set >0 in *'
2482
write(*,*) '* the banner of the event file or in setpara.f. *'
2483
write(*,*) '* *'
2484
write(*,*) '* PROGRAM STOPS HERE *'
2485
write(*,*) '*****************************************************'
2486
ELSE
2487
write(*,*) '*****************************************************'
2488
write(*,*) '* *'
2489
write(*,*) '* >>>>ERROR TRAP CALLED<<<< *'
2490
write(*,*) '* *'
2491
write(*,*) '* SOME ERROR OCCURRED *'
2492
write(*,*) '* *'
2493
write(*,*) '*****************************************************'
2494
ENDIF
2495
2496
STOP
2497
2498
RETURN
2499
END
2500
2501
2502
SUBROUTINE CHECK_NAN(x)
2503
C-------------------------------------------------------
2504
C Check that x is real positive number
2505
C-------------------------------------------------------
2506
C
2507
IMPLICIT NONE
2508
C
2509
C ARGUMENTS
2510
C
2511
real*8 x
2512
c
2513
c LOCAL
2514
c
2515
integer n
2516
data n/0/
2517
SAVE
2518
C
2519
C BEGIN
2520
C
2521
if(.not.(x.gt.0d0).and.x.ne.0) then
2522
n=n+1
2523
x=0d0
2524
write(*,*) 'Found total: ',n,' errors in points in PS'
2525
2526
c open(unit=20,file='decay_error.log',status="old",err=100)
2527
c write(unit=20,fmt='(a50,1x,i5)') 'error in one point in PS',n
2528
c close(20)
2529
c return
2530
c 100 open(unit=20,file='decay_error.log',status="new")
2531
c write(unit=20,fmt='(a50,1x,i5)') 'error in one point in PS',n
2532
c close(20)
2533
2534
endif
2535
2536
return
2537
end
Loggerhead is a web-based interface for Breezy
Version: 2.0.1