1
subroutine decay_event(P5,wgt,ni,ic)
2
c********************************************************************
3
c Decay particle(s) in the event
5
c input: P(0,MaxParticles) :four momenta
7
c ni :number of particle in the event
9
c ic(6,MaxParticles):particle labels
11
c output:P(0,MaxParticles) :four momenta
13
c ic(6,MaxParticles):particle labels
15
c which particle has to be decayed is passed through
16
c a common block DEC_ID in decay.inc
18
c********************************************************************
24
integer ni, ic(7,MaxParticles)
25
double precision P5(0:4,MaxParticles),wgt
29
integer i,k,jj, jhel(MaxDecPart), idecay, idec(MaxParticles),j, im
30
integer JD(MaxDecPart),ICOL(2,MaxDecPart)
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
43
data iseed/1/ !no effect if already called
47
data s_unw,s_wei,s_ove/3*0d0/
51
include 'calc_values.inc'
53
c first map the momenta into objects 0:3
62
c-- First find out how many particles there are in the
65
call find_particle(ip,ni,ic,k,idec)
66
if (k.eq.0) return !no particle to decay
68
c-- Loop over particles to decay
72
c-- find the largest label present for the color flow
73
maxcol=max(ic(4,1),ic(5,1))
75
maxcol=max(maxcol,ic(4,i),ic(5,i))
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
83
pcms(i) = -p(i,1) - p(i,2)
86
call boostx(p(0,idecay),pcms,pd(0,1))
87
jhel(1) = ic(7,idecay) ! helicity of particle
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)
97
if(weight / r .gt. mxwgt) then
100
if (weight .gt. mxwgt) then
101
s_ove = s_ove-mxwgt+weight
106
c-- Multiply by the branching ratio
108
c-- Boost back to lab frame and put momentum into P(0,x)
112
do i=2,MaxDecPart !loop over decay products
113
call boostx(pd(0,i),pcms,p(0,i-1+ni))
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
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
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
132
c finally map the P back to 0:4
138
c p5(4,i)=pdgmass(ic(1,i))
142
p5(4,i)=pdgmass(ic(1,i))
149
SUBROUTINE EVENT(PD,JHEL,JD,ICOL,DWGT)
150
c**********************************************************
151
c Calculation of the decay event quantities
153
c input: pd(0:3,1) decay particle momentum
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
161
c**********************************************************
167
real*8 pd(0:3,MaxDecPart)
168
integer jhel(MaxDecPart),JD(MaxDecPart),ICOL(2,MaxDecPart)
174
include 'calc_values.inc'
178
real*8 pswgt,emmesq,fudge,factor
179
real*8 aa,aa1,aa2,aa3
180
integer isign,jl,jvl,i,j
190
data iseed/1/ !no effect if already called
207
c write(*,*) 'from event: imode,ip',imode,ip
209
If(abs(ip).eq.6) then
210
*-----------------------------------------------------
212
*-----------------------------------------------------
216
*------------------------
219
*------------------------
226
jd(2)=isign*5 !b or b~
227
jd(3)=isign*24 !w+ or w-
229
icol(1,2)=icol(1,1) ! color of the top
230
icol(2,2)=icol(2,1) !anti-color of the top
234
c-- color*(bose factor)*number of helicities
237
jhel(2) = get_hel(xran1(iseed),2)
238
jhel(3) = get_hel(xran1(iseed),3)
240
call phasespace(pd,pswgt)
241
if(pswgt.eq.0d0) return
243
if(isign.eq. 1) call emme_ffv(pd,jhel,emmesq)
244
if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)
246
dwgt=pswgt*emmesq*factor
247
c-- check that dwgt is a reasonable number
253
if(imode.ge.2.and.imode.le.8) then
254
*------------------------
257
*------------------------
267
icol(1,2)=icol(1,1) ! color of the top
268
icol(2,2)=icol(2,1) !anti-color of the top
270
if(isign.eq.1) then !t
277
jd(2)=isign*5 !b or b~
279
if(imode.ge.5) aa =dble(xran1(iseed))
280
if(imode.eq.8) aa1=dble(xran1(iseed))
282
c various cases imode=2,3,4,5,6,7,8
288
jd(jvl)=isign*12 !ve or ve~
289
jd(jl) =isign*(-11) !e+ or e-
291
elseif(imode.eq.3) then
292
*--------------------
294
*--------------------
295
jd(jvl)=isign*14 !vm or vm~
296
jd(jl) =isign*(-13) !mu+ or mu-
298
elseif(imode.eq.4) then
302
if(isign.eq.1) then !t
309
jd(jvl)=isign*16 !vt or vt~
310
jd(jl) =isign*(-15) !ta+ or ta-
312
elseif(imode.eq.5) then
314
* t -> b vl l+ (e+mu)
317
jd(jvl)=isign*12 !ve or ve~
318
jd(jl) =isign*(-11) !e+ or e-
320
jd(jvl)=isign*14 !vm or vm~
321
jd(jl) =isign*(-13) !mu+ or mu-
324
elseif(imode.eq.6) then
326
* t -> b vl l+ (e+mu+ta)
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-
335
if(isign.eq.1) then !t
342
jd(jvl)=isign*16 !vt or vt~
343
jd(jl) =isign*(-15) !ta+ or ta-
346
elseif(imode.eq.7) then
351
icol(1,3)=cindex !position 3 is always a particle
353
icol(1,4)=0 !position 4 is always a anti-particle
357
jd(jvl)=isign*2 !u or u~
358
jd(jl) =isign*(-1) !d~ or d
360
if(isign.eq.1) then !t
367
jd(jvl)=isign*4 !c or c~
368
jd(jl) =isign*(-3) !s~ or s
371
elseif(imode.eq.8) then
376
c first decide if W decays leptonically or hadronically
378
if(aa1.lt.3d0*br_w_lv) then !leptonic decay
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-
387
if(isign.eq.1) then !t
394
jd(jvl)=isign*16 !vt or vt~
395
jd(jl) =isign*(-15) !ta+ or ta-
398
else ! hadronic decay
401
icol(1,3)=cindex !position 3 is always a particle
403
icol(1,4)=0 !position 4 is always a anti-particle
407
jd(jvl)=isign*2 !u or u~
408
jd(jl) =isign*(-1) !d~ or d
410
if(isign.eq.1) then !t
417
jd(jvl)=isign*4 !c or c~
418
jd(jl) =isign*(-3) !s~ or s
422
endif !imode(from 2 to 8)
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
433
jhel(2) = get_hel(xran1(iseed),2)
434
jhel(3) = get_hel(xran1(iseed),2)
437
if(GV.eq.0d0) call error_trap("GV")
438
call phasespace(pd,pswgt)
439
if(pswgt.eq.0d0) return
441
if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq) !t
442
if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !t~
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
450
dwgt=dwgt*2d0/br_w_jj !jj=ud,cs
452
dwgt=dwgt*3d0/(3d0*br_w_lv) !l=e,mu,tau
455
c-- check that dwgt is a reasonable number
458
write(*,*) 'non-existing imode for top decays'
463
If(abs(ip).eq.15) then
464
*-----------------------------------------------------
466
*-----------------------------------------------------
470
*------------------------
472
*------------------------
482
if(isign.eq.1) then !tau-
489
jd(2)=isign*16 !vt or vt~
494
if(imode.eq.3) aa=dble(xran1(iseed))
500
jd(jvl)=isign*(-12) !ve~ or ve
501
jd(jl) =isign*(11) !e- or e+
503
elseif(imode.eq.2) then
504
*-----------------------
506
*-----------------------
507
jd(jvl)=isign*(-14) !vm~ or vm
508
jd(jl) =isign*( 13) !mu- or mu+
510
elseif(imode.eq.3) then
511
*----------------------
512
* tau -> vt l vl (e+mu)
513
*----------------------
515
jd(jvl)=isign*(-12) !ve~ or ve
516
jd(jl) =isign*(11) !e- or e+
518
jd(jvl)=isign*(-14) !vm~ or vm
519
jd(jl) =isign*( 13) !mu- or mu+
524
c-- color*(bose factor)*number of helicities
526
c-- helicities: by definition of 3rd and 4th particle
531
call phasespace(pd,pswgt)
532
if(pswgt.eq.0d0) return
534
if(isign.eq. 1) call emme_f3f(pd,jhel,emmesq) !tau-
535
if(isign.eq.-1) call emme_fx3f(pd,jhel,emmesq) !tau+
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
544
elseif(imode.eq.4) then
545
*------------------------
547
*------------------------
554
jd(2)=isign*16 !vt or vt~
555
jd(3)=isign*(-211) !pi- pi+
559
c-- color*(bose factor)*number of helicities
561
c-- fudge to normalize
562
fudge=lwidth*br_ta_pi/0.00373
565
jhel(2) = get_hel(xran1(iseed),2)
566
jhel(3) = 0 !not used
568
call phasespace(pd,pswgt)
569
if(pswgt.eq.0d0) return
571
if(isign.eq. 1) call emme_ffs (pd,jhel,emmesq)
572
if(isign.eq.-1) call emme_fxfs(pd,jhel,emmesq)
574
dwgt=pswgt*emmesq*factor
575
c-- check that dwgt is a reasonable number
579
elseIf(imode.eq.5) then
580
*------------------------
582
*------------------------
589
jd(2)=isign*16 !vt or vt~
590
jd(3)=isign*(-213) !rho- rho +
594
c-- color*(bose factor)*number of helicities
596
c-- fudge to normalize
597
fudge=lwidth*br_ta_ro/0.0182
600
jhel(2) = get_hel(xran1(iseed),2)
601
jhel(3) = get_hel(xran1(iseed),3)
603
call phasespace(pd,pswgt)
604
if(pswgt.eq.0d0) return
606
if(isign.eq. 1) call emme_ffv (pd,jhel,emmesq)
607
if(isign.eq.-1) call emme_fxfv(pd,jhel,emmesq)
609
dwgt=pswgt*emmesq*factor
610
c-- check that dwgt is a reasonable number
614
write(*,*) 'non-existing imode for tau decays'
620
*-----------------------------------------------------
622
*-----------------------------------------------------
631
c-- rnd number for flavour sums
632
aa =dble(xran1(iseed))
633
aa1=dble(xran1(iseed))
643
elseif(imode.eq.2) then
644
*--------------------
646
*--------------------
650
elseif(imode.eq.3) then
658
elseif(imode.eq.4) then
670
elseif(imode.eq.5) then
672
* z->e- e+,mu-mu+,ta-ta+
674
if(aa.lt.one/Three) then
677
elseif(aa.lt.two/three) then
686
elseif(imode.eq.6) then
687
*------------------------
689
*------------------------
695
if(aa.lt.one/Three) then
698
elseif(aa.lt.two/three) then
707
elseif(imode.eq.7) then
708
*------------------------
710
*------------------------
712
icol(1,2)=cindex !position 2 is always a particle
714
icol(1,3)=0 !position 3 is always a anti-particle
724
elseif(imode.eq.8) then
725
*------------------------
727
*------------------------
729
icol(1,2)=cindex !position 2 is always a particle
731
icol(1,3)=0 !position 3 is always a anti-particle
742
elseif(imode.eq.9) then
743
*------------------------
745
*------------------------
747
icol(1,2)=cindex !position 2 is always a particle
749
icol(1,3)=0 !position 3 is always a anti-particle
752
c-- probability of a uc or ds decay
754
xprob1=w_z_uu/(w_z_dd + w_z_uu)
756
if(aa1.lt.xprob1) then ! decay into ups
759
if(aa.lt. .5d0) then !u
773
else !decay into downs
775
if(aa.lt..5d0) then !d
789
endif !which decay: in ups or downs
791
elseif(imode.eq.10) then
792
*------------------------
794
*------------------------
796
icol(1,2)=cindex !position 2 is always a particle
798
icol(1,3)=0 !position 3 is always a anti-particle
801
c-- probability of a uc or ds decay
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)
806
if(aa1.lt.xprob1) then ! decay into ups
808
if(aa.lt.0.5d0) then !u
822
elseif(aa1.lt.xprob2) then ! decay into downs
824
if(aa.lt. .5d0) then !d
837
else ! decay into b's
847
write(*,*) 'non-existing imode for Z decays'
851
c-- color*(bose factor)*number of helicities
853
if(imode.ge.7) factor=factor*3d0 !quark color in jj
855
jhel(2) = get_hel(xran1(iseed),2)
856
jhel(3) = get_hel(xran1(iseed),2)
858
call phasespace(pd,pswgt)
859
if(pswgt.eq.0d0) return
861
call emme_vff(pd,jhel,emmesq)
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
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
872
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
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
887
c-- check that dwgt is a reasonable number
892
If(abs(ip).eq.24) then
893
*-----------------------------------------------------
895
*-----------------------------------------------------
905
if(isign.eq.1) then !W+
913
if(imode.ge.4) aa =dble(xran1(iseed))
914
if(imode.eq.7) aa1=dble(xran1(iseed))
918
*------------------------
920
*------------------------
921
jd(jvl)=isign*12 !ve or ve~
922
jd(jl) =isign*(-11) !e+ or e-
924
elseif(imode.eq.2) then
925
*--------------------
927
*--------------------
928
jd(jvl)=isign*14 !vm or vm~
929
jd(jl) =isign*(-13) !mu+ or mu-
931
elseif(imode.eq.3) then
935
if(isign.eq.1) then !t
942
jd(jvl)=isign*16 !vt or vt~
943
jd(jl) =isign*(-15) !ta+ or ta-
945
elseif(imode.eq.4) then
950
jd(jvl)=isign*12 !ve or ve~
951
jd(jl) =isign*(-11) !e+ or e-
953
jd(jvl)=isign*14 !vm or vm~
954
jd(jl) =isign*(-13) !mu+ or mu-
957
elseif(imode.eq.5) then
959
* w -> vl l+ (e+mu+ta)
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-
968
if(isign.eq.1) then !t
975
jd(jvl)=isign*16 !vt or vt~
976
jd(jl) =isign*(-15) !ta+ or ta-
979
elseif(imode.eq.6) then
984
icol(1,2)=cindex !position 2 is always a particle
986
icol(1,3)=0 !position 3 is always a anti-particle
990
jd(jvl)=isign*2 !u or u~
991
jd(jl) =isign*(-1) !d~ or d
993
if(isign.eq.1) then !t
1000
jd(jvl)=isign*4 !c or c~
1001
jd(jl) =isign*(-3) !s~ or s
1004
elseif(imode.eq.7) then
1005
*-------------------
1007
*-------------------
1009
c first decide if W decays leptonically or hadronically
1011
if(aa1.lt.3d0*br_w_lv) then !leptonic decay
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-
1020
if(isign.eq.1) then !t
1027
jd(jvl)=isign*16 !vt or vt~
1028
jd(jl) =isign*(-15) !ta+ or ta-
1031
hadro_dec=.true. !hadronic decay
1033
icol(1,2)=cindex !position 2 is always a particle
1035
icol(1,3)=0 !position 3 is always a anti-particle
1039
jd(jvl)=isign*2 !u or u~
1040
jd(jl) =isign*(-1) !d~ or d
1042
if(isign.eq.1) then !t
1049
jd(jvl)=isign*4 !c or c~
1050
jd(jl) =isign*(-3) !s~ or s
1055
c done with the modes: now I start the calculation
1058
c-- color*(bose factor)*number of helicities
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
1063
jhel(2) = get_hel(xran1(iseed),2)
1066
call phasespace(pd,pswgt)
1067
if(pswgt.eq.0d0) return
1069
call emme_vff(pd,jhel,emmesq)
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
1077
dwgt=dwgt*2d0/br_w_jj !jj=ud,cs
1079
dwgt=dwgt*3d0/(3d0*br_w_lv) !l=e,mu,tau
1083
c-- check that dwgt is a reasonable number
1084
call check_nan(dwgt)
1091
*-----------------------------------------------------
1093
*-----------------------------------------------------
1096
*------------------------
1098
*------------------------
1107
icol(1,2)=cindex !position 2 is always a particle
1109
icol(1,3)=0 !position 3 is always a anti-particle
1114
c-- color*(bose factor)*number of helicities
1117
jhel(2) = get_hel(xran1(iseed),2)
1120
call phasespace(pd,pswgt)
1121
if(pswgt.eq.0d0) return
1123
call emme_hff(pd,jhel,emmesq)
1125
dwgt=pswgt*emmesq*factor
1126
c-- check that dwgt is a reasonable number
1127
call check_nan(dwgt)
1130
elseif(imode.eq.2) then
1131
*------------------------
1133
*------------------------
1144
c-- color*(bose factor)*number of helicities
1147
jhel(2) = get_hel(xran1(iseed),2)
1150
call phasespace(pd,pswgt)
1151
if(pswgt.eq.0d0) return
1153
call emme_hff(pd,jhel,emmesq)
1155
dwgt=pswgt*emmesq*factor
1156
c-- check that dwgt is a reasonable number
1157
call check_nan(dwgt)
1161
elseif(imode.eq.3) then
1162
*------------------------
1164
*------------------------
1173
GXX(1)=ghtau(1)/lmass*0.105658389d0
1174
GXX(2)=ghtau(2)/lmass*0.105658389d0
1175
c-- color*(bose factor)*number of helicities
1178
jhel(2) = get_hel(xran1(iseed),2)
1181
call phasespace(pd,pswgt)
1182
if(pswgt.eq.0d0) return
1184
call emme_hff(pd,jhel,emmesq)
1186
dwgt=pswgt*emmesq*factor
1187
c-- check that dwgt is a reasonable number
1188
call check_nan(dwgt)
1191
elseif(imode.eq.4) then
1192
*------------------------
1194
*------------------------
1203
icol(1,2)=cindex !position 2 is always a particle
1205
icol(1,3)=0 !position 3 is always a anti-particle
1210
c-- color*(bose factor)*number of helicities
1213
jhel(2) = get_hel(xran1(iseed),2)
1214
jhel(3) = get_hel(xran1(iseed),2)
1216
call phasespace(pd,pswgt)
1217
if(pswgt.eq.0d0) return
1219
call emme_hff(pd,jhel,emmesq)
1221
dwgt=pswgt*emmesq*factor
1222
c-- check that dwgt is a reasonable number
1223
call check_nan(dwgt)
1226
elseif(imode.eq.5) then
1227
*------------------------
1229
*------------------------
1240
c-- color*(bose factor)*number of helicities
1243
jhel(2) = get_hel(xran1(iseed),2)
1244
jhel(3) = get_hel(xran1(iseed),2)
1246
call phasespace(pd,pswgt)
1247
if(pswgt.eq.0d0) return
1249
call emme_hff(pd,jhel,emmesq)
1251
dwgt=pswgt*emmesq*factor
1252
c-- check that dwgt is a reasonable number
1253
call check_nan(dwgt)
1256
elseif(imode.eq.6) then
1257
*------------------------
1259
*------------------------
1274
c-- color*(bose factor)*number of helicities
1276
c-- fudge factor for normalization
1277
factor=factor*hmass*2d0*pi*(SMBRG*SMWDTH)
1279
jhel(2) = get_hel(xran1(iseed),2)
1280
jhel(3) = get_hel(xran1(iseed),2)
1282
call phasespace(pd,pswgt)
1283
if(pswgt.eq.0d0) return
1285
call emme_hvv(pd,jhel,emmesq)
1287
dwgt=pswgt*emmesq*factor
1288
c-- check that dwgt is a reasonable number
1289
call check_nan(dwgt)
1292
elseif(imode.eq.7) then
1293
*------------------------
1295
*------------------------
1306
c-- color*(bose factor)*number of helicities
1308
c-- fudge factor for normalization
1309
factor=factor*hmass*16d0*pi*(SMBRGA*SMWDTH)
1311
jhel(2) = get_hel(xran1(iseed),2)
1312
jhel(3) = get_hel(xran1(iseed),2)
1314
call phasespace(pd,pswgt)
1315
if(pswgt.eq.0d0) return
1317
call emme_hvv(pd,jhel,emmesq)
1319
dwgt=pswgt*emmesq*factor
1320
c-- check that dwgt is a reasonable number
1321
call check_nan(dwgt)
1324
elseif(imode.eq.8) then
1325
*------------------------
1327
*------------------------
1337
c-- color*(bose factor)*number of helicities
1339
c-- fudge factor for normalization
1340
factor=factor*SMBRZGA*SMWDTH*
1341
&8d0*pi*hmass/(1d0-(zmass/hmass)**2)
1343
jhel(2) = get_hel(xran1(iseed),3) !-1,0,1
1344
jhel(3) = get_hel(xran1(iseed),2)
1346
call twomom (pd(0,1),pd(0,2),pd(0,3),pswgt)
1348
call phasespace(pd,pswgt)
1349
if(pswgt.eq.0d0) return
1351
call emme_hvv(pd,jhel,emmesq)
1353
dwgt=pswgt*emmesq*factor
1354
c-- check that dwgt is a reasonable number
1355
call check_nan(dwgt)
1358
elseif(imode.eq.9) then
1359
*------------------------
1361
*------------------------
1371
c-- color*(bose factor)*number of 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
1379
call phasespace(pd,pswgt)
1380
if(pswgt.eq.0d0) return
1382
call emme_hvv(pd,jhel,emmesq)
1384
dwgt=pswgt*emmesq*factor
1385
c-- check that dwgt is a reasonable number
1386
call check_nan(dwgt)
1390
elseif(imode.eq.10) then
1391
*--------------------------------
1392
* h -> w* w -> l vl l' vl' (l,l'=e,mu)
1393
*--------------------------------
1404
aa1=dble(xran1(iseed))
1405
aa2=dble(xran1(iseed))
1407
if(aa1.lt.half) then
1415
if(aa2.lt.half) then
1428
c-- color*(bose factor)*number of helicities
1436
if(GV.eq.0d0) call error_trap("GV")
1437
call phasespace(pd,pswgt)
1438
if(pswgt.eq.0d0) return
1440
call emme_h4f(pd,jhel,emmesq)
1442
dwgt=pswgt*emmesq*factor
1445
c-- check that dwgt is a reasonable number
1446
call check_nan(dwgt)
1449
elseif(imode.eq.11) then
1450
*--------------------------------
1451
* h -> w* w -> l vl l' vl' (l,l'=e,mu,ta)
1452
*--------------------------------
1463
aa1=dble(xran1(iseed))
1464
aa2=dble(xran1(iseed))
1466
if(aa1.lt.one/three) then
1469
elseif(aa1.lt.two/three) then
1478
if(aa2.lt.one/three) then
1481
elseif(aa2.lt.two/three) then
1496
c-- color*(bose factor)*number of helicities
1504
if(GV.eq.0d0) call error_trap("GV")
1505
call phasespace(pd,pswgt)
1506
if(pswgt.eq.0d0) return
1508
call emme_h4f(pd,jhel,emmesq)
1510
dwgt=pswgt*emmesq*factor
1513
c-- check that dwgt is a reasonable number
1514
call check_nan(dwgt)
1517
elseif(imode.eq.12) then
1518
*--------------------------------
1519
* h -> w* w -> j j l vl (jj=ud,cs;l=e,mu)
1520
*--------------------------------
1531
icol(1,2)=cindex !position 2 is always a particle
1533
icol(1,3)=0 !position 3 is always a anti-particle
1536
aa1=dble(xran1(iseed))
1537
aa2=dble(xran1(iseed))
1538
aa3=dble(xran1(iseed))
1540
if(aa1.lt.half) then
1542
if(aa2.lt.half) then
1551
if(aa3.lt.half) then
1561
if(aa2.lt.half) then
1570
if(aa3.lt.half) then
1585
c-- color*(bose factor)*number of helicities
1593
if(GV.eq.0d0) call error_trap("GV")
1594
call phasespace(pd,pswgt)
1595
if(pswgt.eq.0d0) return
1597
call emme_h4f(pd,jhel,emmesq)
1599
dwgt=pswgt*emmesq*factor
1602
c-- sum of W+ and W- possibilities
1604
c-- check that dwgt is a reasonable number
1605
call check_nan(dwgt)
1608
elseif(imode.eq.13) then
1609
*--------------------------------
1610
* h -> w* w -> j j l vl (jj=ud,cs;l=e,mu,ta)
1611
*--------------------------------
1622
icol(1,2)=cindex !position 2 is always a particle
1624
icol(1,3)=0 !position 3 is always a anti-particle
1627
aa1=dble(xran1(iseed))
1628
aa2=dble(xran1(iseed))
1629
aa3=dble(xran1(iseed))
1631
if(aa1.lt.half) then
1633
if(aa2.lt.half) then
1642
if(aa3.lt.one/three) then
1645
elseif(aa2.lt.two/three) then
1656
if(aa2.lt.half) then
1665
if(aa3.lt.one/three) then
1668
elseif(aa3.lt.two/three) then
1684
c-- color*(bose factor)*number of helicities
1692
if(GV.eq.0d0) call error_trap("GV")
1693
call phasespace(pd,pswgt)
1694
if(pswgt.eq.0d0) return
1696
call emme_h4f(pd,jhel,emmesq)
1698
dwgt=pswgt*emmesq*factor
1701
c-- sum of W+ and W- possibilities
1703
c-- check that dwgt is a reasonable number
1704
call check_nan(dwgt)
1707
elseif(imode.eq.14) then
1708
*------------------------
1710
*------------------------
1720
c-- color*(bose factor)*number of helicities
1723
jhel(2) = get_hel(xran1(iseed),3)
1724
jhel(3) = get_hel(xran1(iseed),3)
1726
call phasespace(pd,pswgt)
1727
if(pswgt.eq.0d0) return
1729
call emme_hvv(pd,jhel,emmesq)
1731
dwgt=pswgt*emmesq*factor
1732
c-- check that dwgt is a reasonable number
1733
call check_nan(dwgt)
1736
elseif(imode.eq.15) then
1737
*---------------------------------
1738
* h -> z* z -> l- l+ l-' l+'(l,l'=e,mu)
1739
*---------------------------------
1750
aa1=dble(xran1(iseed))
1751
aa2=dble(xran1(iseed))
1753
if(aa1.lt.half) then
1761
if(aa2.lt.half) then
1774
c-- color*(bose factor)*number of helicities
1775
factor=1d0*1d0*4d0/2d0
1778
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1781
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1784
if(GV.eq.0d0) call error_trap("GV")
1785
call phasespace(pd,pswgt)
1786
if(pswgt.eq.0d0) return
1788
call emme_h4f(pd,jhel,emmesq)
1790
dwgt=pswgt*emmesq*factor
1793
c-- check that dwgt is a reasonable number
1794
call check_nan(dwgt)
1798
elseif(imode.eq.16) then
1799
*---------------------------------
1800
* h -> z* z -> l- l+ l-' l+'(l,l'=e,mu,ta)
1801
*---------------------------------
1812
aa1=dble(xran1(iseed))
1813
aa2=dble(xran1(iseed))
1815
if(aa1.lt.one/three) then
1818
elseif(aa1.lt.two/three) then
1828
if(aa2.lt.one/three) then
1831
elseif(aa2.lt.two/three) then
1846
c-- color*(bose factor)*number of helicities
1847
factor=1d0*1d0*4d0/2d0
1850
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1853
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1856
if(GV.eq.0d0) call error_trap("GV")
1857
call phasespace(pd,pswgt)
1858
if(pswgt.eq.0d0) return
1860
call emme_h4f(pd,jhel,emmesq)
1862
dwgt=pswgt*emmesq*factor
1865
c-- check that dwgt is a reasonable number
1866
call check_nan(dwgt)
1869
elseif(imode.eq.17) then
1870
*---------------------------------
1871
* h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu )
1872
*---------------------------------
1883
icol(1,2)=cindex !position 2 is always a particle
1885
icol(1,3)=0 !position 3 is always a anti-particle
1894
aa =dble(xran1(iseed))
1895
aa1=dble(xran1(iseed))
1896
aa2=dble(xran1(iseed))
1898
c-- probability of a uc or ds decay
1900
xprob1=w_z_uu/(w_z_dd + w_z_uu)
1902
if(aa1.lt.xprob1) then ! decay into ups
1905
if(aa.lt. .5d0) then !u
1919
else !decay into downs
1921
if(aa.lt..5d0) then !d
1935
endif !which decay: in ups or downs
1938
if(aa2.lt.half) then
1946
c-- color*(bose factor)*number of helicities
1947
factor=3d0*1d0*4d0/2d0
1950
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
1953
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
1956
if(GV.eq.0d0) call error_trap("GV")
1957
call phasespace(pd,pswgt)
1958
if(pswgt.eq.0d0) return
1960
call emme_h4f(pd,jhel,emmesq)
1962
dwgt=pswgt*emmesq*factor
1965
if(multi_decay.eq.1) then
1966
dwgt=dwgt*2d0/xprob1 !j=u,c
1968
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
1970
c-- sum of two possibilities for the decay of a Z
1972
c-- check that dwgt is a reasonable number
1973
call check_nan(dwgt)
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
*---------------------------------
1990
icol(1,2)=cindex !position 2 is always a particle
1992
icol(1,3)=0 !position 3 is always a anti-particle
2001
aa =dble(xran1(iseed))
2002
aa1=dble(xran1(iseed))
2003
aa2=dble(xran1(iseed))
2005
xprob1=w_z_uu/(w_z_dd + w_z_uu)
2006
if(aa1.lt.xprob1) then ! decay into ups
2009
if(aa.lt. .5d0) then !u
2023
else !decay into downs
2025
if(aa.lt..5d0) then !d
2039
endif !which decay: in ups or downs
2042
if(aa2.lt.one/three) then
2045
elseif(aa2.lt.two/three) then
2056
c-- color*(bose factor)*number of helicities
2057
factor=3d0*1d0*4d0/2d0
2060
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2063
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2066
if(GV.eq.0d0) call error_trap("GV")
2067
call phasespace(pd,pswgt)
2068
if(pswgt.eq.0d0) return
2070
call emme_h4f(pd,jhel,emmesq)
2072
dwgt=pswgt*emmesq*factor
2075
if(multi_decay.eq.1) then
2076
dwgt=dwgt*2d0/xprob1 !j=u,c
2078
dwgt=dwgt*2d0/(1.-xprob1) !j=d,s
2080
c-- sum of two possibilities for the decay of a Z
2082
c-- check that dwgt is a reasonable number
2083
call check_nan(dwgt)
2086
elseif(imode.eq.19) then
2087
*---------------------------------
2088
* h -> z* z -> b b~ l- l+ (l=e,mu)
2089
*---------------------------------
2100
icol(1,2)=cindex !position 2 is always a particle
2102
icol(1,3)=0 !position 3 is always a anti-particle
2113
aa=dble(xran1(iseed))
2123
c-- color*(bose factor)*number of helicities
2124
factor=3d0*1d0*8d0/2d0
2126
jhel(2) = get_hel(xran1(iseed),2)
2127
jhel(3) = get_hel(xran1(iseed),2)
2128
jhel(4) = get_hel(xran1(iseed),2)
2131
if(GV.eq.0d0) call error_trap("GV")
2132
call phasespace(pd,pswgt)
2133
if(pswgt.eq.0d0) return
2135
call emme_h4f(pd,jhel,emmesq)
2137
dwgt=pswgt*emmesq*factor
2140
c-- sum of two possibilities for the decay of a Z
2142
c-- check that dwgt is a reasonable number
2143
call check_nan(dwgt)
2146
elseif(imode.eq.20) then
2147
*---------------------------------
2148
* h -> z* z -> b b~ l- l+ (l=e,mu,ta )
2149
*---------------------------------
2160
icol(1,2)=cindex !position 2 is always a particle
2162
icol(1,3)=0 !position 3 is always a anti-particle
2170
c-- color*(bose factor)*number of helicities
2175
aa=dble(xran1(iseed))
2176
if(aa.lt.one/three) then
2179
elseif(aa.lt.two/three) then
2189
c-- color*(bose factor)*number of helicities
2190
factor=3d0*1d0*8d0/2d0
2192
jhel(2) = get_hel(xran1(iseed),2)
2193
jhel(3) = get_hel(xran1(iseed),2)
2194
jhel(4) = get_hel(xran1(iseed),2)
2197
if(GV.eq.0d0) call error_trap("GV")
2198
call phasespace(pd,pswgt)
2199
if(pswgt.eq.0d0) return
2201
call emme_h4f(pd,jhel,emmesq)
2203
dwgt=pswgt*emmesq*factor
2206
c-- sum of two possibilities for the decay of a Z
2208
c-- check that dwgt is a reasonable number
2209
call check_nan(dwgt)
2212
elseif(imode.eq.21) then
2213
*---------------------------------
2214
* h -> z* z -> vl vl~ l-' l+'(l=e,mu,ta,l'=e,mu)
2215
*---------------------------------
2232
aa1=dble(xran1(iseed))
2233
aa2=dble(xran1(iseed))
2235
if(aa1.lt.one/three) then
2238
elseif(aa1.lt.two/three) then
2246
if(aa2.lt.half) then
2254
c-- color*(bose factor)*number of helicities
2255
factor=1d0*1d0*4d0/2d0
2258
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2261
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2264
if(GV.eq.0d0) call error_trap("GV")
2265
call phasespace(pd,pswgt)
2266
if(pswgt.eq.0d0) return
2268
call emme_h4f(pd,jhel,emmesq)
2270
dwgt=pswgt*emmesq*factor
2273
c-- sum of two possibilities for the decay of a Z
2275
c-- check that dwgt is a reasonable number
2276
call check_nan(dwgt)
2279
elseif(imode.eq.22) then
2280
*---------------------------------
2281
* h -> z* z -> vl vl~ l-' l+'(l,l'=e,mu,ta)
2282
*---------------------------------
2298
c-- color*(bose factor)*number of helicities
2301
aa1=dble(xran1(iseed))
2302
aa2=dble(xran1(iseed))
2304
if(aa1.lt.one/three) then
2307
elseif(aa1.lt.two/three) then
2315
if(aa2.lt.one/three) then
2318
elseif(aa2.lt.two/three) then
2329
c-- color*(bose factor)*number of helicities
2330
factor=1d0*1d0*4d0/2d0
2333
if(xran1(iseed) .gt. .5e0) jhel(2)=-1
2336
if(xran1(iseed) .gt. .5e0) jhel(4)=-1
2339
if(GV.eq.0d0) call error_trap("GV")
2340
call phasespace(pd,pswgt)
2341
if(pswgt.eq.0d0) return
2343
call emme_h4f(pd,jhel,emmesq)
2345
dwgt=pswgt*emmesq*factor
2348
c-- sum of two possibilities for the decay of a Z
2350
c-- check that dwgt is a reasonable number
2351
call check_nan(dwgt)
2358
write(*,*) 'from decay:end of decay reached'
2366
INTEGER FUNCTION GET_HEL(RND,ID)
2367
**********************************************************
2368
* GIVEN THE RND NUMBER, RETURNS
2371
**********************************************************
2384
if(rnd.gt.0.5e0) get_hel=1
2385
elseif(id.eq.3) then
2386
if(rnd.lt.0.333333e0)then
2388
elseif(rnd.lt.0.666666e0) then
2394
write(*,*) 'choice not implemented in gen_hel'
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-------------------------------------------------------
2411
PARAMETER (NDMX=50,MXDIM=10)
2419
INTEGER i,j,k,ia(MXDIM)
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
2434
DATA idum/1/ ! DOES NOT AFFECT PREVIOUS SETTING
2436
c-- start the loop over dimensions
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
2444
ia(j)=max(min(int(xn),NDMX),1)
2446
xo=xi(ia(j),j)-xi(ia(j)-1,j)
2447
rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
2452
x(j)=region(j)+rc*dx(j)
2462
SUBROUTINE ERROR_TRAP(STRING)
2463
C-------------------------------------------------------
2464
C RETURNS INFORMATION ABOUT SOME ERROR THAT OCCURED
2465
C-------------------------------------------------------
2475
IF(STRING(1:2).EQ."GV") THEN
2476
write(*,*) '*****************************************************'
2478
write(*,*) '* >>>>ERROR TRAP CALLED<<<< *'
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. *'
2484
write(*,*) '* PROGRAM STOPS HERE *'
2485
write(*,*) '*****************************************************'
2487
write(*,*) '*****************************************************'
2489
write(*,*) '* >>>>ERROR TRAP CALLED<<<< *'
2491
write(*,*) '* SOME ERROR OCCURRED *'
2493
write(*,*) '*****************************************************'
2502
SUBROUTINE CHECK_NAN(x)
2503
C-------------------------------------------------------
2504
C Check that x is real positive number
2505
C-------------------------------------------------------
2521
if(.not.(x.gt.0d0).and.x.ne.0) then
2524
write(*,*) 'Found total: ',n,' errors in points in PS'
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
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