~maddevelopers/mg5amcnlo/2.9.4 : revision 36.21.85

1

program decay

2

c***********************************************************************

3

c 1. reads events from file *

4

c 2. asks which particle to decay *

5

c 3. asks the decay mode *

6

c 4. decays particle *

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***********************************************************************

35

implicit none

36

include 'decay.inc'

37

c

38

c parameters

39

c

40

real*8 smallmass

41

parameter (smallmass=1d0)

42

c

43

c Local

44

c

45

c

46

integer nexternal, ic(7,MaxParticles),idecay(MaxParticles)

47

double precision P(0:4,MaxParticles),wgt

48

integer i,k,iten

49

integer nevent,nfound(0:MaxParticles)

50

data nevent/0/

51

c

52

real*8 sum(0:maxievent),maxwgt(0:maxievent)

53

character*50 dec_name

54

character*3 name

55

character*132 buff

56

character*140 buff2

57

integer iseed

58

real*8 aa

59

real*8 scale,aqcd,aqed

60

integer ievent,eventnr,eventnumber

61

integer maxpup

62

parameter(maxpup=100)

63

integer idbmup(2),pdfgup(2),pdfsup(2),idwtup,nprup,lprup

64

double precision ebmup(2),xsecup,xerrup,xmaxup

65

C

66

logical done,firsttime,fopenend,lwrite,newbanner

67

LOGICAL DEBUG

68

data firsttime/.true./

69

data done/.false./

70

DATA DEBUG/.TRUE./

71

c

72

c external

73

c

74

real*8 dot

75

real xran1

76

c

77

c Global

78

c

79

include 'coupl.inc'

80

data id /201*' '/

81

c

82

character*70 infile,outfile

83

integer run_mode

84

common/mode/run_mode,infile

85

c

86

integer iunit

87

common/to_lh/iunit

88

89

integer sizeievent

90

common/ievents/sizeievent

91

c

92

c Date

93

c

94

c CHARACTER*20 DATE

95

c EXTERNAL DATE_Y2KBUG

96

97

c-----

98

c Begin Code

99

c-----

100

c

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 *'

108

write(*,*) '* *'

109

write(*,*) '* 27-July-2006 *'

110

write(*,*) '*****************************************************'

111

c

112

c run mode: 0 calculates partial widths

113

c 1 decays event in file

114

c

115

write(*,*)

116

write(*,*) ' Input run mode:'

117

write(*,*) ' ---------------'

118

write(*,*)

119

write(*,*) ' 0 = calculates decay widths'

120

write(*,*) ' 1 = decay events in file'

121

write(*,*)

122

123

read(*,'(i1)') run_mode

124

125

c

126

c initialize rnd number for decay

127

c

128

call rnd_init(lunrnd,iseed)

129

c

130

c

131

c open scratch file which will contain the param_card.dat

132

c

133

open(lunp,status='scratch')

134

c

135

if(run_mode.eq.0) then

136

c

137

c calculates the decay width at LO

138

c

139

write(*,*) '*****************************************************'

140

write(*,*) '* run_mode=0 => calculting decay widths *'

141

write(*,*) '*****************************************************'

142

write(*,*)

143

write(*,*) ' Input parameter file (e.g. ../Cards/param_card.dat)'

144

write(*,*) ' ---------------------------------------------------'

145

write(*,*)

146

read(*,'(a)') infile

147

open(lunr,file=infile,status='old',err=102)

148

write(*,*)

149

write(*,*) '*****************************************************'

150

write(*,*) '* Using file: *'

151

write(*,'(1x,a1,a25,a1)') ' ',infile,' '

152

write(*,*) '* for the input params. *'

153

write(*,*) '* *'

154

write(*,*) '* >>>>>>Total widths are recalculated here<<<<< *'

155

write(*,*) '*****************************************************'

156

write(*,*)

157

158

goto 103

159

102 write(*,*) 'file', infile ,' not found : stopping here'

160

stop

161

103 continue

162

163

c

164

c write the param information into the scratch file

165

c

166

done=.false.

167

do while(.not.done)

168

read(lunr,'(a132)',err=101,end=101) buff

169

write(lunp,'(a132)') buff

170

enddo

171

c

172

101 continue

173

174

elseif(run_mode.eq.1) then

175

176

write(*,*) '*****************************************************'

177

write(*,*) '* run_mode=1 => decaying events *'

178

write(*,*) '* *'

179

write(*,*) '* Using the param_card.dat in the banner for *'

180

write(*,*) '* the input params. *'

181

write(*,*) '* *'

182

write(*,*) '* >>>>>>Total widths are recalculated here<<<<< *'

183

write(*,*) '*****************************************************'

184

c

185

write(*,*) 'input event file: (e.g. events.lhe)'

186

read(*,'(a)') infile

187

c

188

write(*,*) 'name for output file: (e.g. dec-events.lhe)'

189

read(*,'(a)') outfile

190

c

191

open(lunr,file=infile,status='old')

192

open(lunw,status='scratch')

193

open(luni,file=outfile,status='unknown')

194

c

195

c copy old banner into new banner and the param_card.dat into the scratch file

196

c

197

done=.false.

198

lwrite=.false.

199

newbanner=.false.

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.

207

c

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

215

enddo

216

c

217

rewind(lunr)

218

c

219

endif

220

221

call setpara

222

call printout

223

call set_mass

224

call set_id

225

c

226

c

227

c

228

write(*,*) '*****************************************************'

229

aa=xran1(iseed)

230

write(*,'(a22,f8.6,a24)') ' * first rnd number= ',aa ,

231

& ' *'

232

write(*,*) '*****************************************************'

233

c

234

c Input the particle to be decayed

235

c

236

write(*,*)

237

write(*,*) ' Implemented decays are for:'

238

write(*,*) ' ---------------------------'

239

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)

246

c

247

c identify the particle ip

248

c

249

ip=0

250

i=-25

251

do while(i.le.25.and.ip.eq.0)

252

if(name.eq.id(i)) ip=i

253

i=i+1

254

enddo

255

c

256

if(ip.eq.0) goto 99 ! particle not existing

257

c

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.

261

c

262

call decay_modes(dec_name)

263

if(imode.eq.0) goto 99 ! decay not implemented

264

c

265

c setting up branching ratio, widths and so on..

266

c

267

call init

268

c

269

c write out decay information

270

c

271

write(*,*)

272

write(*,*) '----------------------------------------------'

273

write(*,*) ' Decay information '

274

write(*,*) '----------------------------------------------'

275

write(*,*) ' '

276

write(*,*) ' particle name : ',id(ip)

277

write(*,*) ' decay mode : ',dec_name

278

c

279

c calculate the MC partial width for normalization purposes

280

c

281

call tot_decay

282

c

283

if(run_mode.eq.0) goto 98 !If only decay width is needed, I'm done

284

c

285

c first check that the particle to be decayed

286

c is present in the events and gather some info

287

c on the events

288

c

289

290

write(*,*) ' '

291

write(*,*) ' Wait....Reading In Event File '

292

write(*,*) ' '

293

done=.false.

294

nevent=0

295

sizeievent=0

296

do i=0,maxievent

297

sum(i)=0d0

298

maxwgt(i)=-1d0

299

enddo

300

do i=0,MaxParticles

301

nfound(i)=0

302

enddo

303

do while (.not. done)

304

call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)

305

if(.not.done) then

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))

311

nevent=nevent+1

312

endif

313

enddo

314

do i=1, sizeievent

315

sum(0)=sum(0)+sum(i)

316

maxwgt(0)=maxwgt(0)+maxwgt(i)

317

enddo

318

319

320

rewind(lunr)

321

322

write(*,*) '----------------------------------------------'

323

write(*,*) ' Input events information '

324

write(*,*) '----------------------------------------------'

325

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

331

write(*,*)

332

333

if(nfound(0).eq.nevent) then

334

write(*,*) ' There are no events with particle ',

335

& id(ip),' to decay !!'

336

stop

337

else

338

do i=0,MaxParticles

339

if(nfound(i).gt.0) then

340

write(*,'(1x,a14,i2,a1,a2,a13,i10)')

341

& ' Events with ',i,' ',id(ip),' :',nfound(i)

342

endif

343

enddo

344

endif

345

346

347

c

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.

353

c

354

write(*,*) ' '

355

write(*,*) '----------------------------------------------'

356

write(*,*) ' Decay events in file : in progress '

357

write(*,*) '----------------------------------------------'

358

write(*,*) ' '

359

done=.false.

360

i = 0

361

do while (.not. done)

362

call read_event(lunr,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)

363

if (.not. done) then

364

i=i+1

365

call decay_event(P,wgt,nexternal,ic)

366

call write_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)

367

endif

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'

371

endif

372

enddo

373

write(*,*) ' 100% done'

374

c

375

c decayed event information

376

c

377

c-- statistics of the decay_event routine

378

c

379

write(*,*)

380

write(*,*) '----------------------------------------------'

381

write(*,*) ' Statistics of the unweighting '

382

write(*,*) '----------------------------------------------'

383

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'

388

IF(DEBUG) then

389

write(*,*) ' Events over = ' ,n_ove

390

write(*,*) ' Integral over = ' ,s_ove/real(n_wei), ' GeV'

391

ENDIF

392

write(*,*) ' '

393

c

394

write(*,*) ' '

395

write(*,*) ' Wait....Writing Out Decayed Event file '

396

write(*,*) ' '

397

rewind(lunw)

398

done=.false.

399

nevent=0

400

do i=0,sizeievent

401

sum(i)=0d0

402

maxwgt(i)=-1d0

403

enddo

404

do while (.not. done)

405

call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)

406

if(.not.done) then

407

eventnumber=eventnr(ievent)

408

sum(eventnumber)=sum(eventnumber)+wgt

409

maxwgt(eventnumber) = max(wgt,maxwgt(eventnumber))

410

nevent=nevent+1

411

endif

412

enddo

413

sum(0)=0d0

414

maxwgt(0)=-1d0

415

do i=1, sizeievent

416

sum(0)=sum(0)+sum(i)

417

maxwgt(0)=maxwgt(0)+maxwgt(i)

418

enddo

419

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

423

iseed=iseed+1

424

close(lunrnd)

425

300 continue

426

cend RF

427

428

429

write(*,*) '----------------------------------------------'

430

write(*,*) ' Output events information '

431

write(*,*) '----------------------------------------------'

432

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

438

write(*,*) ' '

439

c

440

if(newbanner) then

441

write(luni,'(a)') "<MGDecayInfo>"

442

else

443

write(luni,'(a)') "#********************************************************************"

444

write(luni,'(a)') '# particle decayed '

445

write(luni,'(a)') "#********************************************************************"

446

endif

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

455

if(newbanner) then

456

write(luni,'(a)') '</MGDecayInfo>'

457

write(luni,'(a)') '</header>'

458

else

459

write(luni,'(a)') '#********************************************************************'

460

write(luni,'(a)') '-->'

461

endif

462

463

464

C Write out compulsory init info

465

write(luni,'(a)') '<init>'

466

467

rewind(lunr)

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

475

do i=1,sizeievent

476

read(lunr,*) xsecup,xerrup,xmaxup,lprup

477

write(luni,91) sum(i),xerrup*sum(i)/xsecup,maxwgt(i),lprup

478

enddo

479

enddo

480

write(luni,'(a)') '</init>'

481

482

90 FORMAT(2i6,2e19.11,2i2,2i6,i2,i3)

483

91 FORMAT(3e19.11,i4)

484

485

rewind(lunw)

486

i=0

487

done = .false.

488

do while (.not. done)

489

call read_event(lunw,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2,done)

490

491

if (.not. done) then

492

call write_event(luni,P,wgt,nexternal,ic,ievent,scale,aqcd,aqed,buff2)

493

494

endif

495

i=i+1

496

enddo

497

498

499

98 continue

500

501

write(luni,'(a)')'</LesHouchesEvents>'

502

503

c

504

c close datafiles

505

c

506

close(luni)

507

close(lunr)

508

close(lunw)

509

510

stop

511

512

99 write(*,*) 'error'

513

stop

514

515

end

516

517

518

519

integer function eventnr(number)

520

c*************************************************************************

521

c function that should be given a number and returns

522

c the number:

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.

531

c 4 etc.

532

c It updates also a common block with the total number

533

c of different values used as the argument of this function

534

c

535

c before the first call to this function sizeievent should

536

c should be set to 0.

537

c*************************************************************************

538

include "decay.inc"

539

integer i,number,test(maxievent)

540

logical exist

541

542

integer sizeievent

543

common/ievents/sizeievent

544

545

save test

546

547

exist = .false.

548

if (sizeievent.lt.1)goto 188

549

do i=1, sizeievent

550

if (number.eq.test(i))then

551

exist=.true.

552

eventnr=i

553

endif

554

enddo

555

556

188 continue

557

if (.not.exist)then

558

sizeievent=sizeievent+1

559

eventnr=sizeievent

560

test(sizeievent)=number

561

endif

562

563

if (sizeievent.ge.maxievent)then

564

write (*,*)

565

& 'Error, too many different event numbers in event file'

566

return

567

endif

568

569

return

570

end

571

572

573

574

575

576

577

subroutine set_id

578

c****************************************

579

c to set the particles ids

580

*****************************************

581

implicit none

582

include 'decay.inc'

583

c

584

c

585

c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~

586

id(1) ='d'

587

id(2) ='u'

588

id(3) ='s'

589

id(4) ='c'

590

id(5) ='b'

591

id(6) ='t'

592

id(-1)='d~'

593

id(-2)='u~'

594

id(-3)='s~'

595

id(-4)='c~'

596

id(-5)='b~'

597

id(-6)='t~'

598

c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~

599

id(11) ='e-'

600

id(12) ='ve'

601

id(13) ='mu-'

602

id(14) ='vm'

603

id(15) ='ta-'

604

id(16) ='vt'

605

id(-11)='e+'

606

id(-12)='ve~'

607

id(-13)='mu+'

608

id(-14)='vm~'

609

id(-15)='ta+'

610

id(-16)='vt~'

611

c-Bosons: g a z w+ w- h

612

id(21) ='g'

613

id(22) ='a'

614

id(23) ='z'

615

id(24) ='w+'

616

id(-24) ='w-'

617

id(25) ='h'

618

return

619

end

620

621

622

subroutine set_mass

623

c****************************************

624

c to set the particles ids

625

*****************************************

626

implicit none

627

include 'decay.inc'

628

include 'coupl.inc'

629

c

630

c

631

c-Quarks: d u s c b t d~ u~ s~ c~ b~ t~

632

pdgmass(1) =0d0

633

pdgmass(2) =0d0

634

pdgmass(3) =0d0

635

pdgmass(4) =cmass

636

pdgmass(5) =bmass

637

pdgmass(6) =tmass

638

pdgmass(-1)=0d0

639

pdgmass(-2)=0d0

640

pdgmass(-3)=0d0

641

pdgmass(-4)=cmass

642

pdgmass(-5)=bmass

643

pdgmass(-6)=tmass

644

c-Leptons: e- mu- ta- ve vm vt e+ mu+ ta+ ve~ vm~ vt~

645

pdgmass(11) =0d0

646

pdgmass(12) =0d0

647

pdgmass(13) =0d0

648

pdgmass(14) =0d0

649

pdgmass(15) =lmass

650

pdgmass(16) =0d0

651

pdgmass(-11)=0d0

652

pdgmass(-12)=0d0

653

pdgmass(-13)=0d0

654

pdgmass(-14)=0d0

655

pdgmass(-15)=lmass

656

pdgmass(-16)=0d0

657

c-Bosons: g a z w+ w- h

658

pdgmass(21) =0d0

659

pdgmass(22) =0d0

660

pdgmass(23) =zmass

661

pdgmass(24) =wmass

662

pdgmass(-24) =wmass

663

pdgmass(25) =hmass

664

return

665

end

666

667

subroutine decay_modes(dec_name)

668

c*********************************************

669

c to set the particles ids

670

c*********************************************

671

implicit none

672

include 'decay.inc'

673

c

674

c argument

675

c

676

character*50 dec_name

677

c

678

c common

679

c

680

include 'coupl.inc'

681

c

682

c local

683

c

684

character*50 string(30)

685

integer i

686

c----------

687

c START

688

c----------

689

690

write(*,*) 'particle to decay :',id(ip)

691

write(*,*)

692

write(*,*) 'Implemented decay modes:'

693

write(*,*) '------------------------'

694

c

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)'

704

705

write(*,'(i2,2x,a40)') (i,string(i),i=1,8)

706

write(*,*) ' '

707

write(*,*) 'your choice is:'

708

read (*,*) imode

709

if(.not.(imode.ge.1.and.imode.le.8)) then

710

write(*,*) 'choice not implemented'

711

imode=0

712

return

713

endif

714

if(imode.eq.1) then !number of particles in the decay

715

nd=2

716

else

717

nd=3

718

endif

719

dec_name=string(imode)

720

m1=tmass

721

return

722

endif

723

c

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)

734

write(*,*) ' '

735

write(*,*) 'your choice is:'

736

read (*,*) imode

737

if(.not.(imode.ge.1.and.imode.le.8)) then

738

write(*,*) 'choice not implemented'

739

imode=0

740

return

741

endif

742

if(imode.eq.1) then !number of particles in the decay

743

nd=2

744

else

745

nd=3

746

endif

747

m1=tmass

748

dec_name=string(imode)

749

return

750

endif

751

c

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)

759

write(*,*) ' '

760

write(*,*) 'your choice is:'

761

read (*,*) imode

762

if(.not.(imode.ge.1.and.imode.le.5)) then

763

write(*,*) 'choice not implemented'

764

imode=0

765

return

766

endif

767

if(imode.le.3) then !number of particles in the decay

768

nd=3

769

else

770

nd=2

771

endif

772

m1=lmass

773

dec_name=string(imode)

774

return

775

endif

776

c

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)

784

write(*,*) ' '

785

write(*,*) 'your choice is:'

786

read (*,*) imode

787

if(.not.(imode.ge.1.and.imode.le.5)) then

788

write(*,*) 'choice not implemented'

789

imode=0

790

endif

791

if(imode.le.3) then !number of particles in the decay

792

nd=3

793

else

794

nd=2

795

endif

796

m1=lmass

797

dec_name=string(imode)

798

return

799

endif

800

c

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'

814

815

write(*,'(i2,2x,a40)') (i,string(i),i=1,10)

816

write(*,*) ' '

817

write(*,*) 'your choice is:'

818

read (*,*) imode

819

if(.not.(imode.ge.1.and.imode.le.12)) then

820

write(*,*) 'choice not implemented'

821

imode=0

822

return

823

endif

824

m1=zmass

825

nd=2 !number of particles in the decay

826

dec_name=string(imode)

827

return

828

endif

829

c

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)'

838

839

write(*,'(i2,2x,a40)') (i,string(i),i=1,7)

840

write(*,*) 'your choice is:'

841

write(*,*) ' '

842

read (*,*) imode

843

if(.not.(imode.ge.1.and.imode.le.7)) then

844

write(*,*) 'choice not implemented'

845

imode=0

846

endif

847

nd=2 !number of particles in the decay

848

m1=wmass

849

dec_name=string(imode)

850

return

851

endif

852

c

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:'

863

write(*,*) ' '

864

read (*,*) imode

865

if(.not.(imode.ge.1.and.imode.le.7)) then

866

write(*,*) 'choice not implemented'

867

imode=0

868

return

869

endif

870

m1=wmass

871

nd=2 !number of particles in the decay

872

dec_name=string(imode)

873

return

874

endif

875

c

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) "

899

900

write(*,'(i2,2x,a50)') (i,string(i),i=1,22)

901

write(*,*) 'your choice is:'

902

write(*,*) ' '

903

read (*,*) imode

904

if(.not.(imode.ge.1.and.imode.le.22)) then

905

write(*,*) 'choice not implemented'

906

imode=0

907

return

908

endif

909

c-- number of decay products

910

ND=4

911

m1=hmass

912

if(imode.le.9.or.imode.eq.14) ND=2

913

dec_name=string(imode)

914

return

915

endif

916

917

write(*,*) 'no decay mode implemented for particle ',id(ip)

918

919

920

return

921

end

922

923

924

double precision function dot(p1,p2)

925

C****************************************************************************

926

C 4-Vector Dot product

927

C****************************************************************************

928

implicit none

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)

931

end

932

933

934

subroutine init

935

c********************************************************

936

c sets up : the branching ratios (best=HDECAY or PDG)

937

c LO partial widths

938

c ND: the number of decay products

939

c********************************************************

940

implicit none

941

c

942

c Include

943

c

944

include 'decay.inc'

945

include 'coupl.inc'

946

include 'calc_values.inc'

947

c

948

c-----

949

cstart

950

c-----

951

c

952

953

c

954

c PDG2002 values for the Branching ratios:

955

c

956

957

c-- tau

958

br_ta_lv =0.175d0

959

br_ta_pi =0.111d0

960

br_ta_ro =0.254d0

961

c-- w

962

br_w_lv =0.1068d0

963

br_w_jj =1d0-3d0*br_w_lv

964

c-- z

965

br_z_ll =0.0336d0

966

br_z_vv =0.2000d0

967

br_z_cc =0.1176d0

968

br_z_bb =0.1514d0

969

br_z_jj =0.6991d0-br_z_cc-br_z_bb

970

971

write(*,*)

972

write(*,*) '----------------------------------------------'

973

write(*,*) ' Relevant branching ratios '

974

write(*,*) '----------------------------------------------'

975

write(*,*)

976

977

c-------------------------------------------------------------------

978

if(abs(ip).eq.6) then !top

979

980

981

write(*,16) 'BR (t -> b w ) = ', 1d0

982

bratio=1d0

983

if(imode.eq.2) then

984

write(*,16) 'BR (w -> e ve ) = ', br_w_lv

985

bratio=br_w_lv

986

elseif(imode.eq.3) then

987

write(*,16) 'BR (w -> mu vm ) = ', br_w_lv

988

bratio=br_w_lv

989

elseif(imode.eq.4) then

990

write(*,16) 'BR (w -> tau vt ) = ', br_w_lv

991

bratio=br_w_lv

992

elseif(imode.eq.5) then

993

write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'

994

bratio=2d0*br_w_lv

995

elseif(imode.eq.6) then

996

write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'

997

bratio=3d0*br_w_lv

998

elseif(imode.eq.7) then

999

write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'

1000

bratio=br_w_jj

1001

elseif(imode.eq.8) then

1002

write(*,15) 'BR (w ->anything) = ',1d0

1003

bratio=1d0

1004

endif

1005

1006

If(imode.eq.1) then !t -> b w+

1007

calc_width=twidth

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

1018

calc_width=twidth

1019

endif

1020

1021

endif

1022

c-------------------------------------------------------------------

1023

if(abs(ip).eq.15) then !tau

1024

if(imode.eq.1) then

1025

write(*,16) 'BR (ta -> vt e ve) = ', br_ta_lv

1026

bratio=br_ta_lv

1027

elseif(imode.eq.2) then

1028

write(*,16) 'BR (ta -> vt mu vm) = ', br_ta_lv

1029

bratio=br_ta_lv

1030

elseif(imode.eq.3) then

1031

write(*,15) 'BR (ta -> vt l vl) = ', 2d0*br_ta_lv,' (l=e,mu)'

1032

bratio=2d0*br_ta_lv

1033

elseif(imode.eq.4) then

1034

write(*,16) 'BR (ta -> vt pi ) = ', br_ta_pi

1035

bratio=br_ta_pi

1036

elseif(imode.eq.5) then

1037

write(*,16) 'BR (ta -> vt rho ) = ', br_ta_ro

1038

bratio=br_ta_ro

1039

endif

1040

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

1049

endif

1050

1051

endif

1052

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

1057

bratio=br_z_ll

1058

elseif(imode.eq.2) then

1059

write(*,16) 'BR (z -> mu- mu+) = ', br_z_ll

1060

bratio=br_z_ll

1061

elseif(imode.eq.3) then

1062

write(*,16) 'BR (z -> ta- ta+) = ', br_z_ll

1063

bratio=br_z_ll

1064

elseif(imode.eq.4) then

1065

write(*,15) 'BR (z -> l- l+ ) = ', 2d0*br_z_ll,' (l=e,mu)'

1066

bratio=2d0*br_z_ll

1067

elseif(imode.eq.5) then

1068

write(*,15) 'BR (z -> l- l+ ) = ', 3d0*br_z_ll,' (l=e,mu,ta)'

1069

bratio=3d0*br_z_ll

1070

elseif(imode.eq.6) then

1071

write(*,15) 'BR (z -> vl vl~ ) = ', br_z_vv,' (l=e,mu,ta)'

1072

bratio=br_z_vv

1073

elseif(imode.eq.7) then

1074

write(*,16) 'BR (z -> b b~ ) = ', br_z_bb

1075

bratio=br_z_bb

1076

elseif(imode.eq.8)then

1077

write(*,16) 'BR (z -> c c~ ) = ', br_z_cc

1078

bratio=br_z_cc

1079

elseif(imode.eq.9)then

1080

write(*,15) 'BR (z -> j j~ ) = ', br_z_jj+br_z_cc,

1081

& ' (j=u+d+c+s)'

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,

1090

& ' (j=l+j)'

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. ) = ',

1094

& 1d0,' (j=l+v+j)'

1095

bratio=1d0

1096

endif

1097

1098

1099

If(imode.eq.1.or.imode.eq.2) then ! z->e- e+

1100

calc_width=w_z_ll

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+

1104

calc_width=w_z_tau

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~

1110

calc_width=w_z_bb

1111

elseif(imode.eq.8) then ! z->c c~

1112

calc_width=w_z_cc

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

1120

calc_width=zwidth

1121

endif

1122

endif

1123

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

1128

bratio=br_w_lv

1129

elseif(imode.eq.2) then

1130

write(*,16) 'BR (w -> mu vm ) = ', br_w_lv

1131

bratio=br_w_lv

1132

elseif(imode.eq.3) then

1133

write(*,16) 'BR (w -> tau vt ) = ', br_w_lv

1134

bratio=br_w_lv

1135

elseif(imode.eq.4) then

1136

write(*,15) 'BR (w -> l vl ) = ', 2d0*br_w_lv,' (l=e,mu)'

1137

bratio=2d0*br_w_lv

1138

elseif(imode.eq.5) then

1139

write(*,15) 'BR (w -> l vl ) = ', 3d0*br_w_lv,' (l=e,mu,ta)'

1140

bratio=3d0*br_w_lv

1141

elseif(imode.eq.6) then

1142

write(*,15) 'BR (w -> j j ) = ', br_w_jj ,' (jj=ud+cs)'

1143

bratio=br_w_jj

1144

elseif(imode.eq.7) then

1145

write(*,15) 'BR (w ->anything) = 1 (=e+mu+ta+ud+cs)'

1146

bratio=1d0

1147

endif

1148

1149

If(imode.eq.1.or.imode.eq.2) then !w-> l ve

1150

calc_width=w_w_nl

1151

elseif(imode.eq.3) then ! w -> ta vt

1152

calc_width=w_w_tau

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

1160

calc_width=wwidth

1161

endif

1162

1163

endif

1164

1165

c-------------------------------------------------------------------

1166

if(ip.eq.25) then ! higgs

1167

1168

write(*,15) 'Higgs mass = ', hmass, ' GeV'

1169

write(*,15) 'Higgs tot width = ', SMWDTH, ' GeV'

1170

if(imode.eq.1) then

1171

write(*,16) 'BR (h -> b b~ ) = ', SMBRB

1172

bratio=SMBRB

1173

elseif(imode.eq.2) then

1174

write(*,16) 'BR (h -> ta+ ta-) = ', SMBRL

1175

bratio=SMBRL

1176

elseif(imode.eq.3) then

1177

write(*,16) 'BR (h -> mu+ mu-) = ', SMBRM

1178

bratio=SMBRM

1179

elseif(imode.eq.4) then

1180

write(*,16) 'BR (h -> c c~ ) = ', SMBRC

1181

bratio=SMBRC

1182

elseif(imode.eq.5) then

1183

write(*,16) 'BR (h -> t t~ ) = ', SMBRT

1184

bratio=SMBRT

1185

elseif(imode.eq.6) then

1186

write(*,16) 'BR (h -> g g ) = ', SMBRG

1187

bratio=SMBRG

1188

elseif(imode.eq.7) then

1189

write(*,16) 'BR (h -> a a ) = ', SMBRGA

1190

bratio=SMBRGA

1191

elseif(imode.eq.8) then

1192

write(*,16) 'BR (h -> z a ) = ', SMBRZGA

1193

bratio=SMBRZGA

1194

elseif(imode.ge.9.and.imode.le.13) then

1195

write(*,16) 'BR (h -> w+ w- ) = ', SMBRW

1196

bratio=SMBRW

1197

elseif(imode.ge.14.and.imode.le.22) then

1198

write(*,16) 'BR (h -> z z ) = ', SMBRZ

1199

bratio=SMBRZ

1200

endif

1201

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

1216

endif

1217

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,

1226

& ' (j=u+d+c+s)'

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,

1231

& ' (j=u+d+c+s)'

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

1250

endif

1251

1252

1253

write(*,16) 'Best tot Br Ratio = ', bratio

1254

1255

If(imode.eq.1) then ! h->b b~

1256

calc_width=w_h_bb

1257

elseif(imode.eq.2) then ! h->ta- ta+

1258

calc_width=w_h_tau

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~

1262

calc_width=w_h_cc

1263

elseif(imode.eq.5) then ! h-> t t~

1264

if(hmass.le.2d0*tmass) then

1265

write(*,*)

1266

write(*,*) 'Error: this decay mode assumes m_h>2 m_t'

1267

write(*,*)

1268

stop

1269

endif

1270

calc_width=w_h_tt

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

1277

write(*,*)

1278

write(*,*) 'Error: this decay mode assumes m_h>m_Z'

1279

write(*,*)

1280

stop

1281

endif

1282

calc_width= SMBRZGA*SMWDTH

1283

elseif(imode.eq.9) then ! h-> w w

1284

if(hmass.le.2d0*wmass) then

1285

write(*,*)

1286

write(*,*) 'Error: this decay mode assumes m_h>2*m_W'

1287

write(*,*)

1288

stop

1289

endif

1290

calc_width=w_h_ww

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

1302

write(*,*)

1303

write(*,*) 'Error: this decay mode assumes m_h>2*m_Z'

1304

write(*,*)

1305

stop

1306

endif

1307

calc_width=w_h_zz

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

1332

endif

1333

1334

calc_br=calc_width/SMWDTH

1335

write(*,16) 'Calc Br Ratio = ', calc_br

1336

1337

endif

1338

1339

1340

1341

1342

1343

1344

15 format( 3x,a,f9.5,a )

1345

16 format( 3x,a,f9.5 )

1346

1347

1348

1349

return

1350

end

1351

1352

1353

subroutine tot_decay

1354

c**********************************************************

1355

c vegas evaluation of the width

1356

c**********************************************************

1357

implicit none

1358

include 'decay.inc'

1359

c

1360

c LOCAL

1361

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

1365

LOGICAL WRITEOUT

1366

DATA WRITEOUT /.TRUE./

1367

c

1368

c EXTERNAL

1369

c

1370

EXTERNAL FXN

1371

c

1372

c-----

1373

c Begin

1374

c-----

1375

c

1376

c

1377

c write out that calculation is going on

1378

c

1379

write(*,*)

1380

write(*,*) '----------------------------------------------'

1381

write(*,*) ' MC Evaluation of the Partial Width '

1382

write(*,*) '----------------------------------------------'

1383

c

1384

c-MC

1385

c

1386

ndim =3*ND-4 !number of dimensions of the phase space

1387

c

1388

c-warm up

1389

c

1390

init =0

1391

ncall =100000

1392

itmx =5

1393

nprn =0

1394

do i=1,ndim

1395

region(i)=0E0

1396

region(ndim+i)=1E0

1397

enddo

1398

call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)

1399

c

1400

c-final run

1401

c

1402

mxwgt =0d0 !reset the maximum weight to zero

1403

ncall =maxpoints

1404

itmx =1 !only one large iteration

1405

nprn =0

1406

init =1

1407

call vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,chi2a)

1408

c

1409

c-result

1410

c

1411

MC_width=tgral

1412

mxwgt=mxwgt*real(ncall)

1413

eff=MC_width/mxwgt*100e0

1414

C

1415

C write results on the screen

1416

C

1417

write(*,*) ' '

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)

1423

c

1424

c write the result on a log file when WRITEOUT=.true.

1425

c

1426

IF(WRITEOUT) THEN

1427

open(lunout,file='all_decay.dat',status='unknown')

1428

call toend(lunout)

1429

diff=0d0

1430

if(calc_width.gt.0d0) then

1431

diff=(calc_width-MC_width)/(calc_width+MC_width)*100d0

1432

endif

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

1436

close(lunout)

1437

ENDIF

1438

1439

1440

return

1441

end

1442

1443

1444

real*4 function fxn(xx,wgt)

1445

c**********************************************************

1446

c vegas evaluation of the width

1447

c**********************************************************

1448

implicit none

1449

include 'decay.inc'

1450

c

1451

c arguments

1452

c

1453

real*4 xx(mxdim),wgt

1454

c

1455

c Local

1456

c

1457

integer IDS(MaxDecPart),ICOL(2,MaxDecPart)

1458

integer jhel(MaxDecPart)

1459

double precision pd(0:3,MaxDecPart),px(0:3),weight

1460

integer i

1461

c

1462

c External

1463

c

1464

real xran1

1465

integer iseed

1466

data iseed/1/ !this value is irrelevant

1467

integer get_hel

1468

c

1469

c Global

1470

c

1471

include 'coupl.inc'

1472

include 'calc_values.inc'

1473

c

1474

c-----

1475

c Begin

1476

c-----

1477

c-- momentum-random if you want to check Lorentz Invariance

1478

PD(0,1)=M1

1479

PD(1,1)=0d0

1480

PD(2,1)=0d0

1481

PD(3,1)=0d0

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)

1486

1487

c-- helicity

1488

if(abs(ip).eq.24.or.ip.eq.23) then

1489

jhel(1) = get_hel(xran1(iseed),3)

1490

else

1491

jhel(1) = get_hel(xran1(iseed),2)

1492

endif

1493

c-- pass information to the common block

1494

do i=1,mxdim

1495

x(i)=xx(i)

1496

enddo

1497

C

1498

CALL EVENT(PD,JHEL,IDS,ICOL,WEIGHT)

1499

C

1500

c write(*,*) 'from fxn: weight=',weight

1501

mxwgt = max(weight*wgt,mxwgt) !max weight = wgt * fxn

1502

fxn = REAL(weight)

1503

C

1504

return

1505

end

1506

1507

1508

subroutine toend(iunit)

1509

c**********************************************************

1510

c read a file until the end

1511

c**********************************************************

1512

ios = 0

1513

dowhile(ios.eq.0)

1514

read(unit=iunit,fmt='(1x)',iostat=ios)

1515

enddo

1516

end

1517

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***********************************************************

1525

implicit none

1526

include 'decay.inc'

1527

c

1528

c Arguments

1529

c

1530

integer n,ipp,nexternal,ic(7,MaxParticles),idecay(MaxParticles)

1531

c

1532

c Local

1533

c

1534

integer i

1535

1536

c-----

1537

c Begin Code

1538

c-----

1539

1540

n=0

1541

do i=1,nexternal

1542

idecay(i)=0

1543

if (ic(1,i) .eq. ipp .and. ic(6,i) .eq. 1 ) then

1544

n=n+1

1545

idecay(n)=i

1546

endif

1547

enddo

1548

1549

return

1550

end

1551

1552

1553

1554

subroutine rnd_init(iunit,iseed)

1555

c***********************************************************

1556

c initialize rnd number generator xran1 by reading iseed.dat

1557

c***********************************************************

1558

implicit none

1559

c

1560

c Arguments

1561

c

1562

integer iunit,iseed

1563

c

1564

c External

1565

c

1566

real xran1

1567

c

1568

c Local

1569

c

1570

real*8 aa

1571

character*50 adummy

1572

logical done

1573

character*(132) buff

1574

c-----

1575

c Begin Code

1576

c-----

1577

1578

c

1579

c calculating decay rates only

1580

c

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 ,' *'

1588

1589

rewind(iunit)

1590

write(iunit,fmt='(i10)') iseed-1

1591

close(iunit)

1592

write(*,*) '*****************************************************'

1593

goto 201

1594

c- if iseed.dat does not exist then set it to -1.

1595

200 iseed=-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(*,*) '*****************************************************'

1602

201 continue

1603

close(iunit)

1604

return

1605

1606

return

1607

99 write (*,*) 'error'

1608

1609

return

1610

end

1611

1612

1613