2
c***********************************************************************
3
c 1. reads events from file *
4
c 2. asks which particle to decay *
5
c 3. asks the decay mode *
7
c 5. write out decayed events. *
8
c----------------------------------------------------------------------*
9
c FIRST VERSION 16-May-2003 *
10
c----------------------------------------------------------------------*
11
c LAST UPDATE 26-Sep-2003: *
12
c - W->anything and t->b anything decays added. *
13
c - rnd number generator seed changes in sequential runs. *
14
c----------------------------------------------------------------------*
15
c LAST UPDATE 07-Nov-2003: *
16
c - exactly turns unweighted evts in unweighted evts. *
17
c - decaying identical particles in one event fixed. *
18
c - bug in the Z->jets quark fractions corrected. *
19
c - error trap routines added *
20
c----------------------------------------------------------------------*
21
c LAST UPDATE 25-Feb-2004: *
22
c - H->ZZ>leptons : wrong extra factor of two deleted *
23
c----------------------------------------------------------------------*
24
c LAST UPDATE 12-Dec-2004: *
25
c - Merging it to the official web template *
26
c----------------------------------------------------------------------*
27
c LAST UPDATE June-10-2005: *
28
c - Updated for use with new rw_events.f format *
29
c LAST UPDATE August-18-2005: *
30
c - Fixed a bug in the previous update *
31
c----------------------------------------------------------------------*
32
c LAST UPDATE July-27-2006: *
33
c - Make it compliant with MadGraph v 4.0 *
34
c***********************************************************************
41
parameter (smallmass=1d0)
46
integer nexternal, ic(7,MaxParticles),idecay(MaxParticles)
47
double precision P(0:4,MaxParticles),wgt
49
integer nevent,nfound(0:MaxParticles)
52
real*8 sum(0:maxievent),maxwgt(0:maxievent)
59
real*8 scale,aqcd,aqed
60
integer ievent,eventnr,eventnumber
63
integer idbmup(2),pdfgup(2),pdfsup(2),idwtup,nprup,lprup
64
double precision ebmup(2),xsecup,xerrup,xmaxup
66
logical done,firsttime,fopenend,lwrite,newbanner
68
data firsttime/.true./
82
character*70 infile,outfile
84
common/mode/run_mode,infile
90
common/ievents/sizeievent
95
c EXTERNAL DATE_Y2KBUG
101
write(*,*) '*****************************************************'
102
write(*,*) '* DECAY *'
103
write(*,*) '* a MadEvent program *'
104
write(*,*) '* for decaying unstable particles *'
105
write(*,*) '* in the Standar Model *'
106
write(*,*) '* --------------------------------- *'
107
write(*,*) '* version compliant with MG_ME_V4.0 *'
109
write(*,*) '* 27-July-2006 *'
110
write(*,*) '*****************************************************'
112
c run mode: 0 calculates partial widths
113
c 1 decays event in file
116
write(*,*) ' Input run mode:'
117
write(*,*) ' ---------------'
119
write(*,*) ' 0 = calculates decay widths'
120
write(*,*) ' 1 = decay events in file'
123
read(*,'(i1)') run_mode
126
c initialize rnd number for decay
128
call rnd_init(lunrnd,iseed)
131
c open scratch file which will contain the param_card.dat
133
open(lunp,status='scratch')
135
if(run_mode.eq.0) then
137
c calculates the decay width at LO
139
write(*,*) '*****************************************************'
140
write(*,*) '* run_mode=0 => calculting decay widths *'
141
write(*,*) '*****************************************************'
143
write(*,*) ' Input parameter file (e.g. ../Cards/param_card.dat)'
144
write(*,*) ' ---------------------------------------------------'
147
open(lunr,file=infile,status='old',err=102)
149
write(*,*) '*****************************************************'
150
write(*,*) '* Using file: *'
151
write(*,'(1x,a1,a25,a1)') ' ',infile,' '
152
write(*,*) '* for the input params. *'
154
write(*,*) '* >>>>>>Total widths are recalculated here<<<<< *'
155
write(*,*) '*****************************************************'
159
102 write(*,*) 'file', infile ,' not found : stopping here'
164
c write the param information into the scratch file
168
read(lunr,'(a132)',err=101,end=101) buff
169
write(lunp,'(a132)') buff
174
elseif(run_mode.eq.1) then
176
write(*,*) '*****************************************************'
177
write(*,*) '* run_mode=1 => decaying events *'
179
write(*,*) '* Using the param_card.dat in the banner for *'
180
write(*,*) '* the input params. *'
182
write(*,*) '* >>>>>>Total widths are recalculated here<<<<< *'
183
write(*,*) '*****************************************************'
185
write(*,*) 'input event file: (e.g. events.lhe)'
188
write(*,*) 'name for output file: (e.g. dec-events.lhe)'
189
read(*,'(a)') outfile
191
open(lunr,file=infile,status='old')
192
open(lunw,status='scratch')
193
open(luni,file=outfile,status='unknown')
195
c copy old banner into new banner and the param_card.dat into the scratch file
200
read(lunr,'(a132)',err=99) buff
201
do while(index(buff,'<init>') .eq. 0)
202
if(index(buff,"<header>") .ne. 0) newbanner=.true.
203
if (index(buff,"<slha>") .ne. 0 .or.
204
$ index(buff,"Begin param_card.dat") .ne. 0) lwrite=.true.
205
if (index(buff,"</slha>") .ne. 0 .or.
206
$ index(buff,"End param_card.dat") .ne. 0) lwrite=.false.
208
if(index(buff,'</header>') .eq. 0 .and.
209
$ (newbanner .or. index(buff,'-->') .eq. 0 ))
210
$ write(luni,'(a)') buff(1:len_trim(buff))
211
if(lwrite) write(lunp,'(a)') buff(1:len_trim(buff))
212
c if(lwrite) write(*,'(a50)') 'found in param_card: ',buff
213
c if(.not.lwrite) write(*,'(a50)') 'found in banner: ',buff
214
read(lunr,'(a132)',err=99) buff
228
write(*,*) '*****************************************************'
230
write(*,'(a22,f8.6,a24)') ' * first rnd number= ',aa ,
232
write(*,*) '*****************************************************'
234
c Input the particle to be decayed
237
write(*,*) ' Implemented decays are for:'
238
write(*,*) ' ---------------------------'
240
write(*,*) ' Leptons: ta- ta+'
241
write(*,*) ' Quarks : t t~'
242
write(*,*) ' Bosons : z w+ w- h'
243
write(*,*) ' Input particle to be decayed (e.g. t~):'
244
read(*,'(a3)',err=99) name
245
call case_trap(name,3)
247
c identify the particle ip
251
do while(i.le.25.and.ip.eq.0)
252
if(name.eq.id(i)) ip=i
256
if(ip.eq.0) goto 99 ! particle not existing
258
c find out whether the decay requested is implemented.
259
c In case it is, printout the possibilities and
260
c ask the user to choose.
262
call decay_modes(dec_name)
263
if(imode.eq.0) goto 99 ! decay not implemented
265
c setting up branching ratio, widths and so on..
269
c write out decay information
272
write(*,*) '----------------------------------------------'
273
write(*,*) ' Decay information '
274
write(*,*) '----------------------------------------------'
276
write(*,*) ' particle name : ',id(ip)
277
write(*,*) ' decay mode : ',dec_name
279
c calculate the MC partial width for normalization purposes
283
if(run_mode.eq.0) goto 98 !If only decay width is needed, I'm done
285
c first check that the particle to be decayed
286
c is present in the events and gather some info
291
write(*,*) ' Wait....Reading In Event File '
303
do while (.not. done)
304
call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
306
call find_particle(ip,nexternal,ic,k,idecay)
307
nfound(k)=nfound(k)+1
308
eventnumber=eventnr(ievent)
309
sum(eventnumber)=sum(eventnumber)+wgt
310
maxwgt(eventnumber) = max(wgt,maxwgt(eventnumber))
316
maxwgt(0)=maxwgt(0)+maxwgt(i)
322
write(*,*) '----------------------------------------------'
323
write(*,*) ' Input events information '
324
write(*,*) '----------------------------------------------'
326
write(*,*) ' Input Event file : ',infile
327
write(*,*) ' Number of Events : ',nevent
328
write(*,*) ' Integrated weight (pb) : ',sum(0)
329
write(*,*) ' Max wgt : ',maxwgt(0)
330
write(*,*) ' Average wgt : ',sum(0)/nevent
333
if(nfound(0).eq.nevent) then
334
write(*,*) ' There are no events with particle ',
335
& id(ip),' to decay !!'
339
if(nfound(i).gt.0) then
340
write(*,'(1x,a14,i2,a1,a2,a13,i10)')
341
& ' Events with ',i,' ',id(ip),' :',nfound(i)
348
c read each event, decay it, and write it out
349
c note that it is NOT assumed that all the events
350
c have the same number of particles in the final
351
c state or the same type. Each event is treated
352
c independently from the others.
355
write(*,*) '----------------------------------------------'
356
write(*,*) ' Decay events in file : in progress '
357
write(*,*) '----------------------------------------------'
361
do while (.not. done)
362
call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
365
call decay_event(P,wgt,nexternal,ic)
366
call write_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)
368
iten=int(real(nevent)/10e0)+1
369
if(real(i)/real(iten).eq.real(i/iten)) then
370
write(*,'(10x,i3,a)')int(real(i)/real(nevent)*100),' % done'
373
write(*,*) ' 100% done'
375
c decayed event information
377
c-- statistics of the decay_event routine
380
write(*,*) '----------------------------------------------'
381
write(*,*) ' Statistics of the unweighting '
382
write(*,*) '----------------------------------------------'
384
write(*,*) ' Events weighted = ' ,n_wei
385
write(*,*) ' Events unweighted = ' ,nevent
386
write(*,*) ' Width (wgt) = ' ,s_wei/real(n_wei), ' GeV'
387
write(*,*) ' Width (unwgt) = ' ,s_unw/real(n_wei), ' GeV'
389
write(*,*) ' Events over = ' ,n_ove
390
write(*,*) ' Integral over = ' ,s_ove/real(n_wei), ' GeV'
395
write(*,*) ' Wait....Writing Out Decayed Event file '
404
do while (.not. done)
405
call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
407
eventnumber=eventnr(ievent)
408
sum(eventnumber)=sum(eventnumber)+wgt
409
maxwgt(eventnumber) = max(wgt,maxwgt(eventnumber))
417
maxwgt(0)=maxwgt(0)+maxwgt(i)
420
c RF: To write the correct iseed in the banner, read it again form file iseed.dat
421
open(lunrnd,file='./iseed.dat',status="old",err=300)
422
read(lunrnd,err=300,end=300,fmt='(i10)') iseed
429
write(*,*) '----------------------------------------------'
430
write(*,*) ' Output events information '
431
write(*,*) '----------------------------------------------'
433
write(*,*) ' Output Event file : ',outfile
434
write(*,*) ' Number of Events : ',nevent
435
write(*,*) ' Integrated weight (pb) : ',sum(0)
436
write(*,*) ' Max wgt : ',maxwgt(0)
437
write(*,*) ' Average wgt : ',sum(0)/nevent
441
write(luni,'(a)') "<MGDecayInfo>"
443
write(luni,'(a)') "#********************************************************************"
444
write(luni,'(a)') '# particle decayed '
445
write(luni,'(a)') "#********************************************************************"
447
write(luni,'(a,a3)') '# particle name : ',id(ip)
448
write(luni,'(a,a)') '# decay mode : ',dec_name
449
write(luni,'(a,e10.5)') '# MC partial width: ',MC_width
450
write(luni,'(a,i10 )') '# Rnd seed : ',iseed
451
write(luni,'(a,i10 )') '# Number of Events : ',nevent
452
write(luni,'(a,e10.5)') '# Integrated weight (pb): ',sum(0)
453
write(luni,'(a,e10.5)') '# Max wgt : ',maxwgt(0)
454
write(luni,'(a,e10.5)') '# Average wgt : ',sum(0)/nevent
456
write(luni,'(a)') '</MGDecayInfo>'
457
write(luni,'(a)') '</header>'
459
write(luni,'(a)') '#********************************************************************'
460
write(luni,'(a)') '-->'
464
C Write out compulsory init info
465
write(luni,'(a)') '<init>'
468
do while(index(buff,'</init') .eq. 0)
469
read(lunr,'(a132)',err=99) buff
470
if (index(buff,'<init') .eq. 0) cycle
471
read(lunr,*) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
472
$ (pdfsup(i),i=1,2),idwtup,nprup
473
write(luni,90) (idbmup(i),i=1,2),(ebmup(i),i=1,2),(pdfgup(i),i=1,2),
474
$ (pdfsup(i),i=1,2),idwtup,nprup
476
read(lunr,*) xsecup,xerrup,xmaxup,lprup
477
write(luni,91) sum(i),xerrup*sum(i)/xsecup,maxwgt(i),lprup
480
write(luni,'(a)') '</init>'
482
90 FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)
483
91 FORMAT(3e19.11,i4)
488
do while (.not. done)
489
call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)
492
call write_event(luni,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)
501
write(luni,'(a)')'</LesHouchesEvents>'
512
99 write(*,*) 'error'
519
integer function eventnr(number)
520
c*************************************************************************
521
c function that should be given a number and returns
523
c 1 for the first number (i.e. the first time this
524
c function gets executed)
525
c 2 for the second number if the second number is
526
c diffent from the first, otherwise it returns 1
527
c 3 for the third number, if the third number is
528
c different from the first or second, otherwise
529
c it returns 1, or 2 depending if the third number
530
c is equal to the first or the second.
532
c It updates also a common block with the total number
533
c of different values used as the argument of this function
535
c before the first call to this function sizeievent should
536
c should be set to 0.
537
c*************************************************************************
539
integer i,number,test(maxievent)
543
common/ievents/sizeievent
548
if (sizeievent.lt.1)goto 188
550
if (number.eq.test(i))then
558
sizeievent=sizeievent+1
560
test(sizeievent)=number
563
if (sizeievent.ge.maxievent)then
565
& 'Error, too many different event numbers in event file'
578
c****************************************
579
c to set the particles ids
580
*****************************************
585
c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~
598
c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~
611
c-Bosons: g a z w+ w- h
623
c****************************************
624
c to set the particles ids
625
*****************************************
631
c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~
644
c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~
657
c-Bosons: g a z w+ w- h
667
subroutine decay_modes(dec_name)
668
c*********************************************
669
c to set the particles ids
670
c*********************************************
676
character*50 dec_name
684
character*50 string(30)
690
write(*,*) 'particle to decay :',id(ip)
692
write(*,*) 'Implemented decay modes:'
693
write(*,*) '------------------------'
695
if(ip.eq.6) then !top
696
string(1) = ' t -> b w+ '
697
string(2) = ' t -> b ve e+'
698
string(3) = ' t -> b vm mu+'
699
string(4) = ' t -> b vt ta+'
700
string(5) = ' t -> b vl l+ (e+mu)'
701
string(6) = ' t -> b vl l+ (e+mu+ta)'
702
string(7) = ' t -> b j j (ud+cs)'
703
string(8) = ' t -> b anything (e+mu+ta+ud+cs)'
705
write(*,'(i2,2x,a40)') (i,string(i),i=1,8)
707
write(*,*) 'your choice is:'
709
if(.not.(imode.ge.1.and.imode.le.8)) then
710
write(*,*) 'choice not implemented'
714
if(imode.eq.1) then !number of particles in the decay
719
dec_name=string(imode)
724
if(ip.eq.-6) then !anti-top
725
string(1) = ' t~ -> b~ w- '
726
string(2) = ' t~ -> b~ e- ve~'
727
string(3) = ' t~ -> b~ mu- vm~'
728
string(4) = ' t~ -> b~ ta- vt~'
729
string(5) = ' t~ -> b~ l- vl~ (e+mu)'
730
string(6) = ' t~ -> b~ l- vl~ (e+mu+ta)'
731
string(7) = ' t~ -> b~ j j (ud+cs)'
732
string(8) = ' t~ -> b~ anything (e+mu+ta+ud+cs)'
733
write(*,'(i2,2x,a40)') (i,string(i),i=1,8)
735
write(*,*) 'your choice is:'
737
if(.not.(imode.ge.1.and.imode.le.8)) then
738
write(*,*) 'choice not implemented'
742
if(imode.eq.1) then !number of particles in the decay
748
dec_name=string(imode)
752
if(ip.eq.15) then ! tau-
753
string(1) = ' ta- -> vt e- ve~ '
754
string(2) = ' ta- -> vt mu- vm~ '
755
string(3) = ' ta- -> vt l- vl~ (e+mu) '
756
string(4) = ' ta- -> vt pi- '
757
string(5) = ' ta- -> vt rho(770)-'
758
write(*,'(i2,2x,a30)') (i,string(i),i=1,5)
760
write(*,*) 'your choice is:'
762
if(.not.(imode.ge.1.and.imode.le.5)) then
763
write(*,*) 'choice not implemented'
767
if(imode.le.3) then !number of particles in the decay
773
dec_name=string(imode)
777
if(ip.eq.-15) then ! tau+
778
string(1) = ' ta+ -> vt~ ve e+ '
779
string(2) = ' ta+ -> vt~ vm mu+ '
780
string(3) = ' ta+ -> vt~ vl l+ (e+mu)'
781
string(4) = ' ta+ -> vt~ pi+ '
782
string(5) = ' ta+ -> vt~ rho(770)+'
783
write(*,'(i2,2x,a30)') (i,string(i),i=1,5)
785
write(*,*) 'your choice is:'
787
if(.not.(imode.ge.1.and.imode.le.5)) then
788
write(*,*) 'choice not implemented'
791
if(imode.le.3) then !number of particles in the decay
797
dec_name=string(imode)
801
if(ip.eq.23) then ! z
802
string(1) = ' z -> e- e+ '
803
string(2) = ' z -> mu- mu+ '
804
string(3) = ' z -> ta- ta+ '
805
string(4) = ' z -> l- l+ (e+mu)'
806
string(5) = ' z -> l- l+ (e+mu+ta )'
807
string(6) = ' z -> vl vl~ (ve+vm+vt)'
808
string(7) = ' z -> b b~ '
809
string(8) = ' z -> c c~ '
810
string(9) = ' z -> j j~ (u+d+c+s) '
811
string(10)= ' z -> j j~ (u+d+c+s+b)'
812
string(11)= ' z -> visible (e+mu+ta+u+d+c+s+b)'
813
string(12)= ' z -> anything'
815
write(*,'(i2,2x,a40)') (i,string(i),i=1,10)
817
write(*,*) 'your choice is:'
819
if(.not.(imode.ge.1.and.imode.le.12)) then
820
write(*,*) 'choice not implemented'
825
nd=2 !number of particles in the decay
826
dec_name=string(imode)
830
if(ip.eq.24) then ! w+
831
string(1) = ' w+ -> e+ ve'
832
string(2) = ' w+ -> mu+ vm'
833
string(3) = ' w+ -> ta+ vt '
834
string(4) = ' w+ -> l+ vl (e+mu)'
835
string(5) = ' w+ -> l+ vl (e+mu+ta)'
836
string(6) = ' w+ -> j j (ud+cs)'
837
string(7) = ' w+ -> anything (e+mu+ta+ud+cs)'
839
write(*,'(i2,2x,a40)') (i,string(i),i=1,7)
840
write(*,*) 'your choice is:'
843
if(.not.(imode.ge.1.and.imode.le.7)) then
844
write(*,*) 'choice not implemented'
847
nd=2 !number of particles in the decay
849
dec_name=string(imode)
853
if(ip.eq.-24) then ! w-
854
string(1) = ' w- -> e- ve~'
855
string(2) = ' w- -> mu- vm~'
856
string(3) = ' w- -> ta- vt~ '
857
string(4) = ' w- -> l- vl~ (e+mu)'
858
string(5) = ' w- -> l- vl~ (e+mu+ta)'
859
string(6) = ' w- -> j j (ud+cs)'
860
string(7) = ' w- -> anything (e+mu+ta+ud+cs)'
861
write(*,'(i2,2x,a40)') (i,string(i),i=1,7)
862
write(*,*) 'your choice is:'
865
if(.not.(imode.ge.1.and.imode.le.7)) then
866
write(*,*) 'choice not implemented'
871
nd=2 !number of particles in the decay
872
dec_name=string(imode)
876
if(ip.eq.25) then ! h
877
string(1) =' h -> b b~ '
878
string(2) =' h -> ta- ta+'
879
string(3) =' h -> mu- mu+'
880
string(4) =' h -> c c~ '
881
string(5) =' h -> t t~ (when m_h>2*m_t)'
882
string(6) =' h -> g g '
883
string(7) =' h -> a a '
884
string(8) =' h -> z a (when m_h> m_z)'
885
string(9) =' h -> w+ w- (when m_h>2*m_w)'
886
string(10)=" h -> w* w -> l vl l' vl' (l,l'=e,mu)"
887
string(11)=" h -> w* w -> l vl l' vl' (l,l'=e,mu,ta)"
888
string(12)=" h -> w* w -> j j l vl (jj=ud,cs;l=e,mu)"
889
string(13)=" h -> w* w -> j j l vl (jj=ud,cs;l=e,mu,ta)"
890
string(14)=' h -> z z (when m_h>2*m_z)'
891
string(15)=" h -> z* z -> l- l+ l-' l+'(l,l'=e,mu)"
892
string(16)=" h -> z* z -> l- l+ l-' l+'(l,l'=e,mu,ta)"
893
string(17)=" h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu )"
894
string(18)=" h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu,ta)"
895
string(19)=" h -> z* z -> b b~ l- l+ (l=e,mu )"
896
string(20)=" h -> z* z -> b b~ l- l+ (l=e,mu,ta )"
897
string(21)=" h -> z* z -> vl vl~ l-' l+'(l=e,mu,ta;l'=e,mu) "
898
string(22)=" h -> z* z -> vl vl~ l-' l+'(l,l'=e,mu,ta) "
900
write(*,'(i2,2x,a50)') (i,string(i),i=1,22)
901
write(*,*) 'your choice is:'
904
if(.not.(imode.ge.1.and.imode.le.22)) then
905
write(*,*) 'choice not implemented'
909
c-- number of decay products
912
if(imode.le.9.or.imode.eq.14) ND=2
913
dec_name=string(imode)
917
write(*,*) 'no decay mode implemented for particle ',id(ip)
924
double precision function dot(p1,p2)
925
C****************************************************************************
926
C 4-Vector Dot product
927
C****************************************************************************
929
double precision p1(0:3),p2(0:3)
930
dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
935
c********************************************************
936
c sets up : the branching ratios (best=HDECAY or PDG)
938
c ND: the number of decay products
939
c********************************************************
946
include 'calc_values.inc'
954
c PDG2002 values for the Branching ratios:
963
br_w_jj =1d0-3d0*br_w_lv
969
br_z_jj =0.6991d0-br_z_cc-br_z_bb
972
write(*,*) '----------------------------------------------'
973
write(*,*) ' Relevant branching ratios '
974
write(*,*) '----------------------------------------------'
977
c-------------------------------------------------------------------
978
if(abs(ip).eq.6) then !top
981
write(*,16) 'BR (t -> b w ) = ', 1d0
984
write(*,16) 'BR (w -> e ve ) = ', br_w_lv
986
elseif(imode.eq.3) then
987
write(*,16) 'BR (w -> mu vm ) = ', br_w_lv
989
elseif(imode.eq.4) then
990
write(*,16) 'BR (w -> tau vt ) = ', br_w_lv
992
elseif(imode.eq.5) then
993
write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'
995
elseif(imode.eq.6) then
996
write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'
998
elseif(imode.eq.7) then
999
write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'
1001
elseif(imode.eq.8) then
1002
write(*,15) 'BR (w ->anything) = ',1d0
1006
If(imode.eq.1) then !t -> b w+
1008
elseif(imode.eq.2.or.imode.eq.3.or.imode.eq.5) then ! t -> b vl l+
1009
calc_width=twidth*(w_w_nl/wwidth)
1010
if(imode.eq.5) calc_width=calc_width*2d0
1011
elseif(imode.eq.4) then ! t -> b vt tau+
1012
calc_width=twidth*(w_w_tau/wwidth)
1013
elseif(imode.eq.6) then ! t -> b vl l+
1014
calc_width=twidth*((2d0*w_w_nl+w_w_tau)/wwidth)
1015
elseif(imode.eq.7) then ! t -> b j j
1016
calc_width=twidth*(w_w_ud+w_w_cs)/wwidth
1017
elseif(imode.eq.8) then ! t -> b anything
1022
c-------------------------------------------------------------------
1023
if(abs(ip).eq.15) then !tau
1025
write(*,16) 'BR (ta -> vt e ve) = ', br_ta_lv
1027
elseif(imode.eq.2) then
1028
write(*,16) 'BR (ta -> vt mu vm) = ', br_ta_lv
1030
elseif(imode.eq.3) then
1031
write(*,15) 'BR (ta -> vt l vl) = ', 2d0*br_ta_lv,' (l=e,mu)'
1033
elseif(imode.eq.4) then
1034
write(*,16) 'BR (ta -> vt pi ) = ', br_ta_pi
1036
elseif(imode.eq.5) then
1037
write(*,16) 'BR (ta -> vt rho ) = ', br_ta_ro
1041
if(imode.eq.1.or.imode.eq.2) then !ta -> vt vl l
1042
calc_width=lwidth*br_ta_lv
1043
elseIf(imode.eq.3) then
1044
calc_width=2d0*lwidth*br_ta_lv
1045
elseIf(imode.eq.5) then ! tau -> vt rho
1046
calc_width=lwidth*br_ta_ro
1047
elseif(imode.eq.4) then ! tau -> vt pi
1048
calc_width=lwidth*br_ta_pi
1053
c-------------------------------------------------------------------
1054
if(ip.eq.23) then ! z
1055
if (imode.eq.1) then
1056
write(*,16) 'BR (z -> e- e+ ) = ', br_z_ll
1058
elseif(imode.eq.2) then
1059
write(*,16) 'BR (z -> mu- mu+) = ', br_z_ll
1061
elseif(imode.eq.3) then
1062
write(*,16) 'BR (z -> ta- ta+) = ', br_z_ll
1064
elseif(imode.eq.4) then
1065
write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'
1067
elseif(imode.eq.5) then
1068
write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'
1070
elseif(imode.eq.6) then
1071
write(*,15) 'BR (z -> vl vl~ ) = ', br_z_vv,' (l=e,mu,ta)'
1073
elseif(imode.eq.7) then
1074
write(*,16) 'BR (z -> b b~ ) = ', br_z_bb
1076
elseif(imode.eq.8)then
1077
write(*,16) 'BR (z -> c c~ ) = ', br_z_cc
1079
elseif(imode.eq.9)then
1080
write(*,15) 'BR (z -> j j~ ) = ', br_z_jj+br_z_cc,
1082
bratio=br_z_jj+br_z_cc
1083
elseif(imode.eq.10)then
1084
write(*,15) 'BR (z -> j j~ ) = ',
1085
& br_z_jj+br_z_bb+br_z_cc,' (j=u+d+c+s+b)'
1086
bratio=br_z_jj+br_z_bb+br_z_cc
1087
elseif(imode.eq.11)then
1088
write(*,15) 'BR (z -> visible) = ',
1089
& br_z_jj+br_z_bb+br_z_cc+3d0*br_z_ll,
1091
bratio=br_z_jj+br_z_bb+br_z_cc+3d0*br_z_ll
1092
elseif(imode.eq.12)then
1093
write(*,15) 'BR (z -> anyth. ) = ',
1099
If(imode.eq.1.or.imode.eq.2) then ! z->e- e+
1101
elseif(imode.eq.4) then ! z->e- e+,mu-mu+,ta-ta+
1102
calc_width=2d0*w_z_ll
1103
elseif(imode.eq.3) then ! z->ta- ta+
1105
elseif(imode.eq.5) then ! z->e- e+,mu-mu+,ta-ta+
1106
calc_width=2d0*w_z_ll+w_z_tau
1107
elseif(imode.eq.6) then ! z-> vl vl~
1108
calc_width=w_z_nn*3d0
1109
elseif(imode.eq.7) then ! z->b b~
1111
elseif(imode.eq.8) then ! z->c c~
1113
elseif(imode.eq.9) then ! z->jj (u+d+c+s)
1114
calc_width=Two*w_z_dd + w_z_uu + w_z_cc
1115
elseif(imode.eq.10) then !* z->jj (u+d+c+s+b)
1116
calc_width=Two*w_z_dd + w_z_uu + w_z_cc + w_z_bb
1117
elseif(imode.eq.11) then !* z->any (l+u+d+c+s+b)
1118
calc_width=zwidth-w_z_nn*3d0
1119
elseif(imode.eq.12) then !* total width
1124
c-------------------------------------------------------------------
1125
if(abs(ip).eq.24) then ! w
1126
if (imode.eq.1) then
1127
write(*,16) 'BR (w -> e ve ) = ', br_w_lv
1129
elseif(imode.eq.2) then
1130
write(*,16) 'BR (w -> mu vm ) = ', br_w_lv
1132
elseif(imode.eq.3) then
1133
write(*,16) 'BR (w -> tau vt ) = ', br_w_lv
1135
elseif(imode.eq.4) then
1136
write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'
1138
elseif(imode.eq.5) then
1139
write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'
1141
elseif(imode.eq.6) then
1142
write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'
1144
elseif(imode.eq.7) then
1145
write(*,15) 'BR (w ->anything) = 1 (=e+mu+ta+ud+cs)'
1149
If(imode.eq.1.or.imode.eq.2) then !w-> l ve
1151
elseif(imode.eq.3) then ! w -> ta vt
1153
elseif(imode.eq.4) then ! w -> l vl (e,mu)
1154
calc_width=2d0*w_w_nl
1155
elseif(imode.eq.5) then ! w -> l vl (e,mu,tau)
1156
calc_width=2d0*w_w_nl+w_w_tau
1157
elseif(imode.eq.6) then ! w -> ud+cs
1158
calc_width=w_w_ud+w_w_cs
1159
elseif(imode.eq.7) then ! w -> anything
1165
c-------------------------------------------------------------------
1166
if(ip.eq.25) then ! higgs
1168
write(*,15) 'Higgs mass = ', hmass, ' GeV'
1169
write(*,15) 'Higgs tot width = ', SMWDTH, ' GeV'
1171
write(*,16) 'BR (h -> b b~ ) = ', SMBRB
1173
elseif(imode.eq.2) then
1174
write(*,16) 'BR (h -> ta+ ta-) = ', SMBRL
1176
elseif(imode.eq.3) then
1177
write(*,16) 'BR (h -> mu+ mu-) = ', SMBRM
1179
elseif(imode.eq.4) then
1180
write(*,16) 'BR (h -> c c~ ) = ', SMBRC
1182
elseif(imode.eq.5) then
1183
write(*,16) 'BR (h -> t t~ ) = ', SMBRT
1185
elseif(imode.eq.6) then
1186
write(*,16) 'BR (h -> g g ) = ', SMBRG
1188
elseif(imode.eq.7) then
1189
write(*,16) 'BR (h -> a a ) = ', SMBRGA
1191
elseif(imode.eq.8) then
1192
write(*,16) 'BR (h -> z a ) = ', SMBRZGA
1194
elseif(imode.ge.9.and.imode.le.13) then
1195
write(*,16) 'BR (h -> w+ w- ) = ', SMBRW
1197
elseif(imode.ge.14.and.imode.le.22) then
1198
write(*,16) 'BR (h -> z z ) = ', SMBRZ
1202
if(imode.eq.10) then
1203
write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'
1204
bratio=bratio*(2d0*br_w_lv)**2
1205
elseif(imode.eq.11) then
1206
write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'
1207
bratio=bratio*(3d0*br_w_lv)**2
1208
elseif(imode.eq.12) then
1209
write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'
1210
write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'
1211
bratio=bratio*2d0*br_w_jj*2d0*br_w_lv
1212
elseif(imode.eq.13) then
1213
write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'
1214
write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'
1215
bratio=bratio*2d0*br_w_jj*3d0*br_w_lv
1218
if(imode.eq.15) then
1219
write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'
1220
bratio=bratio*(2d0*br_z_ll)**2
1221
elseif(imode.eq.16) then
1222
write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'
1223
bratio=bratio*(3d0*br_z_ll)**2
1224
elseif(imode.eq.17) then
1225
write(*,15) 'BR (z -> j j~ ) = ', br_z_jj+br_z_cc,
1227
write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'
1228
bratio=bratio*(br_z_jj+br_z_cc)*2d0*br_z_ll*2d0
1229
elseif(imode.eq.18) then
1230
write(*,15) 'BR (z -> j j~ ) = ', br_z_jj+br_z_cc,
1232
write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'
1233
bratio=bratio*(br_z_jj+br_z_cc)*3d0*br_z_ll*2d0
1234
elseif(imode.eq.19) then
1235
write(*,16) 'BR (z -> b b~ ) = ', br_z_bb
1236
write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'
1237
bratio=bratio*br_z_bb*2d0*br_z_ll*2d0
1238
elseif(imode.eq.20) then
1239
write(*,16) 'BR (z -> b b~ ) = ', br_z_bb
1240
write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'
1241
bratio=bratio*br_z_bb*3d0*br_z_ll*2d0
1242
elseif(imode.eq.21) then
1243
write(*,15) 'BR (z -> vl vl~ ) = ', br_z_vv,' (l=e,mu,ta)'
1244
write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'
1245
bratio=bratio*br_z_vv*2d0*br_z_ll*2d0
1246
elseif(imode.eq.22) then
1247
write(*,15) 'BR (z -> vl vl~ ) = ', br_z_vv,' (l=e,mu,ta)'
1248
write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'
1249
bratio=bratio*br_z_vv*3d0*br_z_ll*2d0
1253
write(*,16) 'Best tot Br Ratio = ', bratio
1255
If(imode.eq.1) then ! h->b b~
1257
elseif(imode.eq.2) then ! h->ta- ta+
1259
elseif(imode.eq.3) then ! h->mu- mu+
1260
calc_width=w_h_tau/(lmass/0.105658389d0)**2
1261
elseif(imode.eq.4) then ! h-> c c~
1263
elseif(imode.eq.5) then ! h-> t t~
1264
if(hmass.le.2d0*tmass) then
1266
write(*,*) 'Error: this decay mode assumes m_h>2 m_t'
1271
elseif(imode.eq.6) then ! h-> g g
1272
calc_width=SMBRG*SMWDTH
1273
elseif(imode.eq.7) then ! h-> a a
1274
calc_width= SMBRGA*SMWDTH
1275
elseif(imode.eq.8) then ! h-> z a
1276
if(hmass.le.zmass) then
1278
write(*,*) 'Error: this decay mode assumes m_h>m_Z'
1282
calc_width= SMBRZGA*SMWDTH
1283
elseif(imode.eq.9) then ! h-> w w
1284
if(hmass.le.2d0*wmass) then
1286
write(*,*) 'Error: this decay mode assumes m_h>2*m_W'
1291
elseif(imode.eq.10) then ! h -> w* w -> e ve mu vmu
1292
calc_width= SMBRW*SMWDTH*(2d0*w_w_nl/wwidth)**2
1293
elseif(imode.eq.11) then ! h -> w* w -> l vl l' vl' (l,l'=e,mu,ta)
1294
calc_width= SMBRW*SMWDTH*((2d0*w_w_nl+w_w_tau)/wwidth)**2
1295
elseif(imode.eq.12) then ! h -> w* w -> j j l vl (jj=ud,cs;l=e,mu)
1296
calc_width= SMBRW*SMWDTH*2d0*(2d0*w_w_nl*(w_w_ud+w_w_cs))/wwidth**2
1297
elseif(imode.eq.13) then !h -> w* w -> j j l vl (jj=ud,cs;l=e,mu,ta)
1298
calc_width= SMBRW*SMWDTH*2d0*
1299
& ((2d0*w_w_nl+w_w_tau)*(w_w_ud+w_w_cs))/wwidth**2
1300
elseif(imode.eq.14) then ! h-> z z
1301
if(hmass.le.2d0*zmass) then
1303
write(*,*) 'Error: this decay mode assumes m_h>2*m_Z'
1308
elseif(imode.eq.15) then ! h -> z* z -> l- l+ l-' l+'(l,l'=e,mu)
1309
calc_width= SMBRZ*SMWDTH*(2d0*w_z_ll/zwidth)**2
1310
elseif(imode.eq.16) then ! h -> z* z -> l- l+ l-' l+'(l,l'=e,mu,ta)
1311
calc_width= SMBRZ*SMWDTH*
1312
& ((2d0*w_z_ll+w_z_tau)/zwidth)**2
1313
elseif(imode.eq.17) then ! h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu )
1314
calc_width= SMBRZ*SMWDTH*2d0*
1315
& (2d0*w_z_ll*(Two*w_z_dd + w_z_uu + w_z_cc))/zwidth**2
1316
elseif(imode.eq.18) then ! h -> z* z -> j j~ l- l+ (j=u,d,c,s;l=e,mu,ta)
1317
calc_width= SMBRZ*SMWDTH*2d0*
1318
& ((2d0*w_z_ll+w_z_tau)*
1319
& (Two*w_z_dd + w_z_uu + w_z_cc))/zwidth**2
1320
elseif(imode.eq.19) then ! h -> z* z -> b b~ l- l+ (l=e,mu )
1321
calc_width=SMBRZ*SMWDTH*2d0*
1322
& w_z_ll*w_z_bb*2d0/zwidth**2
1323
elseif(imode.eq.20) then ! h -> z* z -> b b~ l- l+ (l=e,mu,ta )
1324
calc_width=SMBRZ*SMWDTH*2d0*
1325
& (2d0*w_z_ll+w_z_tau)*w_z_bb/zwidth**2
1326
elseif(imode.eq.21) then ! h -> z* z -> vl vl~ l-' l+'(l=e,mu,ta;l'=e,mu)
1327
calc_width=SMBRZ*SMWDTH*2d0*
1328
& 3d0*w_z_nn*2d0*w_z_ll/zwidth**2
1329
elseif(imode.eq.22) then ! h -> z* z -> vl vl~ l-' l+'(l,l'=e,mu,ta)
1330
calc_width=SMBRZ*SMWDTH*2d0*
1331
& 3d0*w_z_nn*(2d0*w_z_ll+w_z_tau)/zwidth**2
1334
calc_br=calc_width/SMWDTH
1335
write(*,16) 'Calc Br Ratio = ', calc_br
1344
15 format( 3x,a,f9.5,a )
1345
16 format( 3x,a,f9.5 )
1353
subroutine tot_decay
1354
c**********************************************************
1355
c vegas evaluation of the width
1356
c**********************************************************
1362
INTEGER init,itmx,ncall,nprn,ndim,i
1363
REAL*4 region(2*mxdim),fxn,tgral,chi2a,sd
1364
real*8 diff,eff,ave_wei
1366
DATA WRITEOUT /.TRUE./
1377
c write out that calculation is going on
1380
write(*,*) '----------------------------------------------'
1381
write(*,*) ' MC Evaluation of the Partial Width '
1382
write(*,*) '----------------------------------------------'
1386
ndim =3*ND-4 !number of dimensions of the phase space
1398
call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)
1402
mxwgt =0d0 !reset the maximum weight to zero
1404
itmx =1 !only one large iteration
1407
call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)
1412
mxwgt=mxwgt*real(ncall)
1413
eff=MC_width/mxwgt*100e0
1415
C write results on the screen
1418
write(*,*) ' LO partial width = ' ,calc_width, ' GeV'
1419
write(*,*) ' MC partial width = ' ,MC_width,' +-',sd,' GeV'
1420
write(*,*) ' Max. wgt = ' ,mxwgt
1421
write(*,1) ' Unwgt. Eff.(%) = ' ,eff
1422
1 format(1x,a21,1x,f4.1)
1424
c write the result on a log file when WRITEOUT=.true.
1427
open(lunout,file='all_decay.dat',status='unknown')
1430
if(calc_width.gt.0d0) then
1431
diff=(calc_width-MC_width)/(calc_width+MC_width)*100d0
1433
write(lunout,'(a3,2x,i2,2x,e12.5,2x,e10.5,
1434
& 1x,f5.2,3x,e10.5,3x,f5.1)')
1435
& id(ip),imode,calc_width,MC_width,diff,mxwgt,eff
1444
real*4 function fxn(xx,wgt)
1445
c**********************************************************
1446
c vegas evaluation of the width
1447
c**********************************************************
1453
real*4 xx(mxdim),wgt
1457
integer IDS(MaxDecPart),ICOL(2,MaxDecPart)
1458
integer jhel(MaxDecPart)
1459
double precision pd(0:3,MaxDecPart),px(0:3),weight
1466
data iseed/1/ !this value is irrelevant
1472
include 'calc_values.inc'
1477
c-- momentum-random if you want to check Lorentz Invariance
1482
PD(1,1)=xran1(iseed)*10d0
1483
PD(2,1)=xran1(iseed)*10d0
1484
PD(3,1)=xran1(iseed)*10d0
1485
PD(0,1)=dsqrt(M1**2+PD(1,1)**2+PD(2,1)**2+PD(3,1)**2)
1488
if(abs(ip).eq.24.or.ip.eq.23) then
1489
jhel(1) = get_hel(xran1(iseed),3)
1491
jhel(1) = get_hel(xran1(iseed),2)
1493
c-- pass information to the common block
1498
CALL EVENT(PD,JHEL,IDS,ICOL,WEIGHT)
1500
c write(*,*) 'from fxn: weight=',weight
1501
mxwgt = max(weight*wgt,mxwgt) !max weight = wgt * fxn
1508
subroutine toend(iunit)
1509
c**********************************************************
1510
c read a file until the end
1511
c**********************************************************
1514
read(unit=iunit,fmt='(1x)',iostat=ios)
1518
subroutine find_particle(ipp,nexternal,ic,n,idecay)
1519
c***********************************************************
1520
c looks for particle ip in the event ic and
1521
c returns the n instances found which
1522
c are not decayed yet. idecay(i) gives the position
1523
c of the i-th particle to be decayed.
1524
c***********************************************************
1530
integer n,ipp,nexternal,ic(7,MaxParticles),idecay(MaxParticles)
1543
if (ic(1,i) .eq. ipp .and. ic(6,i) .eq. 1 ) then
1554
subroutine rnd_init(iunit,iseed)
1555
c***********************************************************
1556
c initialize rnd number generator xran1 by reading iseed.dat
1557
c***********************************************************
1573
character*(132) buff
1579
c calculating decay rates only
1581
c- try to open the iseed.dat file
1582
open(iunit,file='./iseed.dat',status="old",err=200)
1583
read(iunit,err=200,end=200,fmt='(i10)') iseed
1584
write(*,*) '*****************************************************'
1585
write(*,*) '* reading seed from iseed.dat *'
1586
write(*,'(a22,i6,a26)')
1587
& ' * rnd number seed = ',iseed ,' *'
1590
write(iunit,fmt='(i10)') iseed-1
1592
write(*,*) '*****************************************************'
1594
c- if iseed.dat does not exist then set it to -1.
1596
open(iunit,file='./iseed.dat',status="new")
1597
write(iunit,fmt='(i10)') iseed-1
1598
write(*,*) '*****************************************************'
1599
write(*,*) '* no iseed.dat => iseed=-1 *'
1600
write(*,*) '* iseed.dat now written *'
1601
write(*,*) '*****************************************************'
1607
99 write (*,*) 'error'