]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/src/hwdcle.f
Change in CookdEdx() by Prashant
[u/mrichter/AliRoot.git] / HERWIG / src / hwdcle.f
CommitLineData
3820ca8e 1
2CDECK ID>, HWDCLE.
3
4*CMZ :- -28/01/92 12.34.44 by Mike Seymour
5
6*-- Author : Luca Stanco
7
8C-----------------------------------------------------------------------
9
10 SUBROUTINE HWDCLE(IHEP)
11
12C-----------------------------------------------------------------------
13
14C INTERFACE TO QQ-CLEO MONTE CARLO (LS 11/12/91)
15
16C-----------------------------------------------------------------------
17
18 INCLUDE 'HERWIG61.INC'
19
20 INTEGER IHEP,IIHEP,NHEPHF,QQLMAT
21
22 LOGICAL QQLERR
23
24 CHARACTER*8 NAME
25
26 EXTERNAL QQLMAT
27
28C---QQ-CLEO COMMON'S
29
30C*** MCPARS.INC
31
32 INTEGER MCTRK, NTRKS, MCVRTX, NVTXS, MCHANS, MCDTRS, MPOLQQ
33
34 INTEGER MCNUM, MCSTBL, MCSTAB, MCTLQQ, MDECQQ
35
36 INTEGER MHLPRB, MHLLST, MHLANG, MCPLST, MFDECA
37
38 PARAMETER (MCTRK = 512)
39
40 PARAMETER (NTRKS = MCTRK)
41
42 PARAMETER (MCVRTX = 256)
43
44 PARAMETER (NVTXS = MCVRTX)
45
46 PARAMETER (MCHANS = 4000)
47
48 PARAMETER (MCDTRS = 8000)
49
50 PARAMETER (MPOLQQ = 300)
51
52 PARAMETER (MCNUM = 500)
53
54 PARAMETER (MCSTBL = 40)
55
56 PARAMETER (MCSTAB = 512)
57
58 PARAMETER (MCTLQQ = 100)
59
60 PARAMETER (MDECQQ = 300)
61
62 PARAMETER (MHLPRB = 500)
63
64 PARAMETER (MHLLST = 1000)
65
66 PARAMETER (MHLANG = 500)
67
68 PARAMETER (MCPLST = 200)
69
70 PARAMETER (MFDECA = 5)
71
72C*** MCPROP.INC
73
74 REAL AMASS, CHARGE, CTAU, SPIN, RWIDTH, RMASMN, RMASMX
75
76 REAL RMIXPP, RCPMIX
77
78 INTEGER NPMNQQ, NPMXQQ, IDMC, INVMC, LPARTY, CPARTY
79
80 INTEGER IMIXPP, ICPMIX
81
82 COMMON/MCMAS1/
83
84 * NPMNQQ, NPMXQQ,
85
86 * AMASS(-20:MCNUM), CHARGE(-20:MCNUM), CTAU(-20:MCNUM),
87
88 * IDMC(-20:MCNUM), SPIN(-20:MCNUM),
89
90 * RWIDTH(-20:MCNUM), RMASMN(-20:MCNUM), RMASMX(-20:MCNUM),
91
92 * LPARTY(-20:MCNUM), CPARTY(-20:MCNUM),
93
94 * IMIXPP(-20:MCNUM), RMIXPP(-20:MCNUM),
95
96 * ICPMIX(-20:MCNUM), RCPMIX(-20:MCNUM),
97
98 * INVMC(0:MCSTBL)
99
100C
101
102 INTEGER NPOLQQ, IPOLQQ
103
104 COMMON/MCPOL1/
105
106 * NPOLQQ, IPOLQQ(5,MPOLQQ)
107
108C
109
110 CHARACTER QNAME*10, PNAME*10
111
112 COMMON/MCNAMS/
113
114 * QNAME(37), PNAME(-20:MCNUM)
115
116C
117
118C*** MCCOMS.INC
119
120 INTEGER NCTLQQ, NDECQQ, IVRSQQ, IORGQQ, IRS1QQ
121
122 INTEGER IEVTQQ, IRUNQQ, IBMRAD
123
124 INTEGER NTRKMC, QQNTRK, NSTBMC, NSTBQQ, NCHGMC, NCHGQQ
125
126 INTEGER IRANQQ, IRANMC, IRANCC, IRS2QQ
127
128 INTEGER IPFTQQ, IPCDQQ, IPRNTV, ITYPEV, IDECSV, IDAUTV
129
130 INTEGER ISTBMC, NDAUTV
131
132 INTEGER IVPROD, IVDECA
133
134 REAL BFLDQQ
135
136 REAL ENERQQ, BEAMQQ, BMPSQQ, BMNGQQ, EWIDQQ, BWPSQQ, BWNGQQ
137
138 REAL BPOSQQ, BSIZQQ
139
140 REAL ECM, P4CMQQ, P4PHQQ, ENERNW, BEAMNW, BEAMP, BEAMN
141
142 REAL PSAV, P4QQ, HELCQQ
143
144 CHARACTER DATEQQ*20, TIMEQQ*20, FOUTQQ*80, FCTLQQ*80, FDECQQ*80
145
146 CHARACTER FGEOQQ*80
147
148 CHARACTER CCTLQQ*80, CDECQQ*80
149
150C
151
152 COMMON/MCCM1A/
153
154 * NCTLQQ, NDECQQ, IVRSQQ, IORGQQ, IRS1QQ(3), BFLDQQ,
155
156 * ENERQQ, BEAMQQ, BMPSQQ, BMNGQQ, EWIDQQ, BWPSQQ, BWNGQQ,
157
158 * BPOSQQ(3), BSIZQQ(3),
159
160 * IEVTQQ, IRUNQQ,
161
162 * IBMRAD, ECM, P4CMQQ(4), P4PHQQ(4),
163
164 * ENERNW, BEAMNW, BEAMP, BEAMN,
165
166 * NTRKMC, QQNTRK, NSTBMC, NSTBQQ, NCHGMC, NCHGQQ,
167
168 * IRANQQ(2), IRANMC(2), IRANCC(2), IRS2QQ(5),
169
170 * IPFTQQ(MCTRK), IPCDQQ(MCTRK), IPRNTV(MCTRK), ITYPEV(MCTRK,2),
171
172 * IDECSV(MCTRK), IDAUTV(MCTRK), ISTBMC(MCTRK), NDAUTV(MCTRK),
173
174 * IVPROD(MCTRK), IVDECA(MCTRK),
175
176 * PSAV(MCTRK,4), HELCQQ(MCTRK), P4QQ(4,MCTRK)
177
178C
179
180 COMMON/MCCM1B/
181
182 * DATEQQ, TIMEQQ, FOUTQQ, FCTLQQ, FDECQQ, FGEOQQ,
183
184 * CCTLQQ(MCTLQQ), CDECQQ(MDECQQ)
185
186 INTEGER IDSTBL
187
188 COMMON/MCCM1C/
189
190 * IDSTBL(MCSTAB)
191
192C
193
194 INTEGER IFINAL(MCTRK), IFINSV(MCSTAB), NFINAL
195
196 EQUIVALENCE (IFINAL,ISTBMC), (IFINSV,IDSTBL), (NFINAL,NSTBMC)
197
198C
199
200 INTEGER NVRTX, ITRKIN, NTRKOU, ITRKOU, IVKODE
201
202 REAL XVTX, TVTX, RVTX
203
204 COMMON/MCCM2/
205
206 * NVRTX, XVTX(MCVRTX,3), TVTX(MCVRTX), RVTX(MCVRTX),
207
208 * ITRKIN(MCVRTX), NTRKOU(MCVRTX), ITRKOU(MCVRTX),
209
210 * IVKODE(MCVRTX)
211
212C*** MCGEN.INC
213
214 INTEGER QQIST,QQIFR,QQN,QQK,QQMESO,QQNC,QQKC,QQLASTN
215
216 REAL QQPUD,QQPS1,QQSIGM,QQMAS,QQPAR,QQCMIX,QQCND,QQBSPI,QQBSYM,QQP
217
218 REAL QQPC,QQCZF
219
220C
221
222 COMMON/DATA1/QQIST,QQIFR,QQPUD,QQPS1,QQSIGM,QQMAS(15),QQPAR(25)
223
224 COMMON/DATA2/QQCZF(15),QQMESO(36),QQCMIX(6,2)
225
226 COMMON/DATA3/QQCND(3)
227
228 COMMON/DATA5/QQBSPI(5),QQBSYM(3)
229
230 COMMON/JET/QQN,QQK(250,2),QQP(250,5),QQNC,QQKC(10),QQPC(10,4),
231
232 * QQLASTN
233
234C---
235
236 IF(FSTEVT) THEN
237
238C---INITIALIZE QQ-CLEO
239
240 CALL QQINIT(QQLERR)
241
242 IF(QQLERR) CALL HWWARN('HWDEUR',500,*999)
243
244 ENDIF
245
246C---CONSTRUCT THE HADRON FOR QQ-CLEO
247
248C NOTE: THE IDPDG CODE IS PROVIDED THROUGH THE QQLMAT ROUTINE
249
250C FROM THE CLEO PACKAGE (QQ-CLEO <--> IDPDG CODE TRANSFORMATION)
251
252 QQN=1
253
254 IDHEP(IHEP)=IDPDG(IDHW(IHEP))
255
256 QQK(1,1)=0
257
258 QQK(1,2)=QQLMAT(IDHEP(IHEP),1)
259
260 QQP(1,1)=PHEP(1,IHEP)
261
262 QQP(1,2)=PHEP(2,IHEP)
263
264 QQP(1,3)=PHEP(3,IHEP)
265
266 QQP(1,5)=AMASS(QQK(1,2))
267
268 QQP(1,4)=SQRT(QQP(1,5)**2+QQP(1,1)**2+QQP(1,2)**2+QQP(1,3)**2)
269
270C---LET QQ-CLEO DO THE JOB
271
272 QQNTRK=0
273
274 NVRTX=0
275
276 CALL DECADD(.FALSE.)
277
278C---UPDATE THE HERWIG TABLE : LOOP OVER QQN-CLEO FINAL PARTICLES
279
280 DO 40 IIHEP=1,QQN
281
282 NHEP=NHEP+1
283
284 ISTHEP(NHEP)=198
285
286 IF(ITYPEV(IIHEP,2).GE.0) ISTHEP(NHEP)=1
287
288 IDHEP(NHEP)=QQLMAT(ITYPEV(IIHEP,1),2)
289
290 CALL HWUIDT(1,IDHEP(NHEP),IDHW(NHEP),NAME)
291
292 IF(IIHEP.EQ.1) THEN
293
294 ISTHEP(IHEP)=199
295
296 JDAHEP(1,IHEP)=NHEP
297
298 JDAHEP(2,IHEP)=NHEP
299
300 ISTHEP(NHEP)=199
301
302 NHEPHF=NHEP
303
304 JMOHEP(1,NHEP)=IHEP
305
306 JMOHEP(2,NHEP)=IHEP
307
308 ELSE
309
310 JMOHEP(1,NHEP)=IPRNTV(IIHEP)+NHEPHF-1
311
312 JMOHEP(2,NHEP)=NHEPHF
313
314 ENDIF
315
316 JDAHEP(1,NHEP)=0
317
318 JDAHEP(2,NHEP)=0
319
320 IF(NDAUTV(IIHEP).GT.0) THEN
321
322 JDAHEP(1,NHEP)=IDAUTV(IIHEP)+NHEPHF-1
323
324 JDAHEP(2,NHEP)=JDAHEP(1,NHEP)+NDAUTV(IIHEP)-1
325
326 ENDIF
327
328 PHEP(1,NHEP)=QQP(IIHEP,1)
329
330 PHEP(2,NHEP)=QQP(IIHEP,2)
331
332 PHEP(3,NHEP)=QQP(IIHEP,3)
333
334 PHEP(4,NHEP)=QQP(IIHEP,4)
335
336 PHEP(5,NHEP)=QQP(IIHEP,5)
337
338 VHEP(1,NHEP)=XVTX(IVPROD(IIHEP),1)
339
340 VHEP(2,NHEP)=XVTX(IVPROD(IIHEP),2)
341
342 VHEP(3,NHEP)=XVTX(IVPROD(IIHEP),3)
343
344 VHEP(4,NHEP)=0.
345
346 40 CONTINUE
347
348 999 END
349
350CDECK ID>, HWDEUR.
351
352*CMZ :- -28/01/92 12.34.44 by Mike Seymour
353
354*-- Author : Luca Stanco
355
356C-----------------------------------------------------------------------
357
358 SUBROUTINE HWDEUR(IHEP)
359
360C-----------------------------------------------------------------------
361
362C INTERFACE TO EURODEC PACKAGE (LS 10/29/91)
363
364C-----------------------------------------------------------------------
365
366 INCLUDE 'HERWIG61.INC'
367
368 INTEGER IHEP,IIHEP,NHEPHF,IEUPDG,IPDGEU
369
370 CHARACTER*8 NAME
371
372C---EURODEC COMMON'S : INITIAL INPUT
373
374 INTEGER EULUN0,EULUN1,EULUN2,EURUN,EUEVNT
375
376 CHARACTER*4 EUDATD,EUTIT
377
378 REAL AMINIE(12),EUWEI
379
380 COMMON/INPOUT/EULUN0,EULUN1,EULUN2
381
382 COMMON/FILNAM/EUDATD,EUTIT
383
384 COMMON/HVYINI/AMINIE
385
386 COMMON/RUNINF/EURUN,EUEVNT,EUWEI
387
388C---EURODEC WORKING COMMON'S
389
390 INTEGER NPMAX,NTMAX
391
392 PARAMETER (NPMAX=18,NTMAX=2000)
393
394 INTEGER EUNP,EUIP(NPMAX),EUPHEL(NPMAX),EUTEIL,EUINDX(NTMAX),
395
396 & EUORIG(NTMAX),EUDCAY(NTMAX),EUTHEL(NTMAX)
397
398 REAL EUAPM(NPMAX),EUPCM(5,NPMAX),EUPVTX(3,NPMAX),EUPTEI(5,NTMAX),
399
400 & EUSECV(3,NTMAX)
401
402 COMMON/MOMGEN/EUNP,EUIP,EUAPM,EUPCM,EUPHEL,EUPVTX
403
404 COMMON/RESULT/EUTEIL,EUPTEI,EUINDX,EUORIG,EUDCAY,EUTHEL,EUSECV
405
406C---EURODEC COMMON'S FOR DECAY PROPERTIES
407
408 INTEGER NGMAX,NCMAX
409
410 PARAMETER (NGMAX=400,NCMAX=9000)
411
412 INTEGER EUNPA,EUIPC(NGMAX),EUIPDG(NGMAX),EUIDP(NGMAX),
413
414 & EUCONV(NCMAX)
415
416 REAL EUPM(NGMAX),EUPLT(NGMAX)
417
418 COMMON/PCTABL/EUNPA,EUIPC,EUIPDG,EUPM,EUPLT,EUIDP
419
420 COMMON/CONVRT/EUCONV
421
422C---
423
424 IF(FSTEVT) THEN
425
426C---CHANGE HERE THE DEFAULT VALUES OF EURODEC COMMON'S
427
428C
429
430C---INITIALIZE EURODEC COMMON'S
431
432CC CALL EUDCIN
433
434C---INITIALIZE EURODEC
435
436 CALL EUDINI
437
438 ENDIF
439
440C---CONSTRUCT THE HADRON FOR EURODEC FROM ID1,ID2
441
442 EUNP=1
443
444 IDHEP(IHEP)=IDPDG(IDHW(IHEP))
445
446 EUIP(1)=IPDGEU(IDHEP(IHEP))
447
448 EUAPM(1)=EUPM(EUCONV(IABS(EUIP(1))))
449
450 EUPCM(1,1)=PHEP(1,IHEP)
451
452 EUPCM(2,1)=PHEP(2,IHEP)
453
454 EUPCM(3,1)=PHEP(3,IHEP)
455
456 EUPCM(5,1)=SQRT(PHEP(1,IHEP)**2+PHEP(2,IHEP)**2+PHEP(3,IHEP)**2)
457
458 EUPCM(4,1)=SQRT(EUPCM(5,1)**2+EUAPM(1)**2)
459
460C NOT POLARIZED HADRONS
461
462 EUPHEL(1)=0
463
464C HADRONS START FROM PRIMARY VERTEX
465
466 EUPVTX(1,1)=0.
467
468 EUPVTX(2,1)=0.
469
470 EUPVTX(3,1)=0.
471
472C---LET EURODEC DO THE JOB
473
474 EUTEIL=0
475
476 CALL FRAGMT(1,1,0)
477
478C---UPDATE THE HERWIG TABLE : LOOP OVER N-EURODEC FINAL PARTICLES
479
480 DO 40 IIHEP=1,EUTEIL
481
482 NHEP=NHEP+1
483
484 ISTHEP(NHEP)=198
485
486 IF(EUDCAY(IIHEP).EQ.0) ISTHEP(NHEP)=1
487
488 IDHEP(NHEP)=IEUPDG(EUINDX(IIHEP))
489
490 CALL HWUIDT(1,IDHEP(NHEP),IDHW(NHEP),NAME)
491
492 IF(IIHEP.EQ.1) THEN
493
494 ISTHEP(IHEP)=199
495
496 JDAHEP(1,IHEP)=NHEP
497
498 JDAHEP(2,IHEP)=NHEP
499
500 ISTHEP(NHEP)=199
501
502 NHEPHF=NHEP
503
504 JMOHEP(1,NHEP)=IHEP
505
506 JMOHEP(2,NHEP)=IHEP
507
508 JDAHEP(1,NHEP)=EUDCAY(IIHEP)/10000+NHEPHF-1
509
510 JDAHEP(2,NHEP)=MOD(EUDCAY(IIHEP),10000)+NHEPHF-1
511
512 ELSE
513
514 JMOHEP(1,NHEP)=MOD(EUORIG(IIHEP),10000)+NHEPHF-1
515
516 JMOHEP(2,NHEP)=NHEPHF
517
518 JDAHEP(1,NHEP)=EUDCAY(IIHEP)/10000+NHEPHF-1
519
520 JDAHEP(2,NHEP)=MOD(EUDCAY(IIHEP),10000)+NHEPHF-1
521
522 ENDIF
523
524 PHEP(1,NHEP)=EUPTEI(1,IIHEP)
525
526 PHEP(2,NHEP)=EUPTEI(2,IIHEP)
527
528 PHEP(3,NHEP)=EUPTEI(3,IIHEP)
529
530 PHEP(4,NHEP)=EUPTEI(4,IIHEP)
531
532 PHEP(5,NHEP)=EUPTEI(5,IIHEP)
533
534 VHEP(1,NHEP)=EUSECV(1,IIHEP)
535
536 VHEP(2,NHEP)=EUSECV(2,IIHEP)
537
538 VHEP(3,NHEP)=EUSECV(3,IIHEP)
539
540 VHEP(4,NHEP)=0.
541
542 IF (IIHEP.GT.NTMAX) CALL HWWARN('HWDEUR',99,*999)
543
544 40 CONTINUE
545
546 999 END
547
548CDECK ID>, HWDFOR.
549
550*CMZ :- -01/04/99 19.52.44 by Mike Seymour
551
552*-- Author : Ian Knowles
553
554C-----------------------------------------------------------------------
555
556 SUBROUTINE HWDFOR(P0,P1,P2,P3,P4)
557
558C-----------------------------------------------------------------------
559
560C Generates 4-body decay 0->1+2+3+4 using pure phase space
561
562C-----------------------------------------------------------------------
563
564 IMPLICIT NONE
565
566 DOUBLE PRECISION HWR,P0(5),P1(5),P2(5),P3(5),P4(5),B,C,AA,BB,
567
568 & CC,DD,EE,TT,S1,RS1,FF,S2,PP,QQ,RR,P1CM,P234(5),P2CM,P34(5),P3CM
569
570 DOUBLE PRECISION TWO
571
572 PARAMETER (TWO=2.D0)
573
574 EXTERNAL HWR
575
576 B=P0(5)-P1(5)
577
578 C=P2(5)+P3(5)+P4(5)
579
580 IF (B.LT.C) CALL HWWARN('HWDFOR',100,*999)
581
582 AA=(P0(5)+P1(5))**2
583
584 BB=B**2
585
586 CC=C**2
587
588 DD=(P3(5)+P4(5))**2
589
590 EE=(P3(5)-P4(5))**2
591
592 TT=(B-C)*P0(5)**7/16
593
594C Select squared masses S1 and S2 of 234 and 34 subsystems
595
596 10 S1=BB+HWR()*(CC-BB)
597
598 RS1=SQRT(S1)
599
600 FF=(RS1-P2(5))**2
601
602 S2=DD+HWR()*(FF-DD)
603
604 PP=(AA-S1)*(BB-S1)
605
606 QQ=((RS1+P2(5))**2-S2)*(FF-S2)/S1
607
608 RR=(S2-DD)*(S2-EE)/S2
609
610 IF (PP*QQ*RR*(FF-DD)**2.LT.TT*S1*S2*HWR()**2) GOTO 10
611
612C Do two body decays: 0-->1+234, 234-->2+34 and 34-->3+4
613
614 P1CM=SQRT(PP/4)/P0(5)
615
616 P234(5)=RS1
617
618 P2CM=SQRT(QQ/4)
619
620 P34(5)=SQRT(S2)
621
622 P3CM=SQRT(RR/4)
623
624 CALL HWDTWO(P0 ,P1,P234,P1CM,TWO,.TRUE.)
625
626 CALL HWDTWO(P234,P2,P34 ,P2CM,TWO,.TRUE.)
627
628 CALL HWDTWO(P34 ,P3,P4 ,P3CM,TWO,.TRUE.)
629
630 999 END
631
632CDECK ID>, HWDFIV.
633
634*CMZ :- -01/04/99 19.52.44 by Mike Seymour
635
636*-- Author : Ian Knowles
637
638C-----------------------------------------------------------------------
639
640 SUBROUTINE HWDFIV(P0,P1,P2,P3,P4,P5)
641
642C-----------------------------------------------------------------------
643
644C Generates 5-body decay 0->1+2+3+4+5 using pure phase space
645
646C-----------------------------------------------------------------------
647
648 IMPLICIT NONE
649
650 DOUBLE PRECISION HWR,P0(5),P1(5),P2(5),P3(5),P4(5),P5(5),B,C,
651
652 & AA,BB,CC,DD,EE,FF,TT,S1,RS1,GG,S2,RS2,HH,S3,PP,QQ,RR,SS,P1CM,
653
654 & P2345(5),P2CM,P345(5),P3CM,P45(5),P4CM
655
656 DOUBLE PRECISION TWO
657
658 PARAMETER (TWO=2.D0)
659
660 EXTERNAL HWR
661
662 B=P0(5)-P1(5)
663
664 C=P2(5)+P3(5)+P4(5)+P5(5)
665
666 IF (B.LT.C) CALL HWWARN('HWDFIV',100,*999)
667
668 AA=(P0(5)+P1(5))**2
669
670 BB=B**2
671
672 CC=C**2
673
674 DD=(P3(5)+P4(5)+P5(5))**2
675
676 EE=(P4(5)+P5(5))**2
677
678 FF=(P4(5)-P5(5))**2
679
680 TT=(B-C)*P0(5)**11/729
681
682C Select squared masses S1, S2 and S3 of 2345, 345 and 45 subsystems
683
684 10 S1=BB+HWR()*(CC-BB)
685
686 RS1=SQRT(S1)
687
688 GG=(RS1-P2(5))**2
689
690 S2=DD+HWR()*(GG-DD)
691
692 RS2=SQRT(S2)
693
694 HH=(RS2-P3(5))**2
695
696 S3=EE+HWR()*(HH-EE)
697
698 PP=(AA-S1)*(BB-S1)
699
700 QQ=((RS1+P2(5))**2-S2)*(GG-S2)/S1
701
702 RR=((RS2+P3(5))**2-S3)*(HH-S3)/S2
703
704 SS=(S3-EE)*(S3-FF)/S3
705
706 IF (PP*QQ*RR*SS*((GG-DD)*(HH-EE))**2.LT.TT*S1*S2*S3*HWR()**2)
707
708 & GOTO 10
709
710C Do two body decays: 0-->1+2345, 2345-->2+345, 345-->3+45 and 45-->4+5
711
712 P1CM=SQRT(PP/4)/P0(5)
713
714 P2345(5)=RS1
715
716 P2CM=SQRT(QQ/4)
717
718 P345(5)=RS2
719
720 P3CM=SQRT(RR/4)
721
722 P45(5)=SQRT(S3)
723
724 P4CM=SQRT(SS/4)
725
726 CALL HWDTWO(P0 ,P1,P2345,P1CM,TWO,.TRUE.)
727
728 CALL HWDTWO(P2345,P2,P345 ,P2CM,TWO,.TRUE.)
729
730 CALL HWDTWO(P345 ,P3,P45 ,P3CM,TWO,.TRUE.)
731
732 CALL HWDTWO(P45 ,P4,P5 ,P4CM,TWO,.TRUE.)
733
734 999 END
735
736CDECK ID>, HWDHAD.
737
738*CMZ :- -26/04/91 14.01.26 by Federico Carminati
739
740*-- Author : Ian Knowles, Bryan Webber & Mike Seymour
741
742C-----------------------------------------------------------------------
743
744 SUBROUTINE HWDHAD
745
746C-----------------------------------------------------------------------
747
748C GENERATES DECAYS OF UNSTABLE HADRONS AND LEPTONS
749
750C-----------------------------------------------------------------------
751
752 INCLUDE 'HERWIG61.INC'
753
754 DOUBLE PRECISION HWR,HWULDO,RN,BF,COSANG,RSUM,DIST(4),VERTX(4),
755
756 & PMIX,WTMX,WTMX2,XS,DOT1,DOT2,HWDPWT,HWDWWT,XXX,YYY
757
758 INTEGER IHEP,ID,MHEP,IDM,I,IDS,IM,MO,IPDG
759
760 LOGICAL STABLE
761
762 EXTERNAL HWR,HWDPWT,HWDWWT,HWULDO
763
764 IF (IERROR.NE.0) RETURN
765
766 DO 100 IHEP=1,NMXHEP
767
768 IF (IHEP.GT.NHEP) THEN
769
770 ISTAT=90
771
772 RETURN
773
774 ELSEIF (ISTHEP(IHEP).EQ.120 .AND.
775
776 & JDAHEP(1,IHEP).EQ.IHEP.AND.JDAHEP(2,IHEP).EQ.IHEP) THEN
777
778C---COPY COLOUR SINGLET CMF
779
780 NHEP=NHEP+1
781
782 IF (NHEP.GT.NMXHEP) CALL HWWARN('HWDHAD',100,*999)
783
784 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
785
786 CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
787
788 IDHW(NHEP)=IDHW(IHEP)
789
790 IDHEP(NHEP)=IDHEP(IHEP)
791
792 ISTHEP(NHEP)=190
793
794 JMOHEP(1,NHEP)=IHEP
795
796 JMOHEP(2,NHEP)=NHEP
797
798 JDAHEP(2,NHEP)=NHEP
799
800 JDAHEP(1,IHEP)=NHEP
801
802 JDAHEP(2,IHEP)=NHEP
803
804 ELSEIF (ISTHEP(IHEP).GE.190.AND.ISTHEP(IHEP).LE.193) THEN
805
806C---FIRST CHECK FOR STABILITY
807
808 ID=IDHW(IHEP)
809
810 IF (RSTAB(ID)) THEN
811
812 ISTHEP(IHEP)=1
813
814 JDAHEP(1,IHEP)=0
815
816 JDAHEP(2,IHEP)=0
817
818C---SPECIAL FOR GAUGE BOSON DECAY
819
820 IF (ID.GE.198.AND.ID.LE.200) CALL HWDBOS(IHEP)
821
822C---SPECIAL FOR HIGGS BOSON DECAY
823
824 IF (ID.EQ.201) CALL HWDHIG(ZERO)
825
826 ELSE
827
828C---UNSTABLE.
829
830C Calculate position of decay vertex
831
832 IF (DKLTM(ID).EQ.ZERO) THEN
833
834 CALL HWVEQU(4,VHEP(1,IHEP),VERTX)
835
836 MHEP=IHEP
837
838 IDM=ID
839
840 ELSE
841
842 CALL HWUDKL(ID,PHEP(1,IHEP),DIST)
843
844 CALL HWVSUM(4,VHEP(1,IHEP),DIST,VERTX)
845
846 IF (MAXDKL) THEN
847
848 CALL HWDXLM(VERTX,STABLE)
849
850 IF (STABLE) THEN
851
852 ISTHEP(IHEP)=1
853
854 JDAHEP(1,IHEP)=0
855
856 JDAHEP(2,IHEP)=0
857
858 GOTO 100
859
860 ENDIF
861
862 ENDIF
863
864 IF (MIXING.AND.(ID.EQ.221.OR.ID.EQ.223.OR.
865
866 & ID.EQ.245.OR.ID.EQ.247)) THEN
867
868C Select flavour of decaying b-meson allowing for flavour oscillation
869
870 IDS=MOD(ID,3)
871
872 XXX=XMRCT(IDS)*DIST(4)/PHEP(4,IHEP)
873
874 YYY=YMRCT(IDS)*DIST(4)/PHEP(4,IHEP)
875
876 IF (ABS(YYY).LT.10) THEN
877
878 PMIX=HALF*(ONE-COS(XXX)/COSH(YYY))
879
880 ELSE
881
882 PMIX=HALF
883
884 ENDIF
885
886 IF (HWR().LE.PMIX) THEN
887
888 IF (ID.LE.223) THEN
889
890 IDM=ID+24
891
892 ELSE
893
894 IDM=ID-24
895
896 ENDIF
897
898 ELSE
899
900 IDM=ID
901
902 ENDIF
903
904C Introduce a decaying neutral b-meson
905
906 IF (NHEP+1.GT.NMXHEP) CALL HWWARN('HWDHAD',101,*999)
907
908 MHEP=NHEP+1
909
910 ISTHEP(MHEP)=ISTHEP(IHEP)
911
912 ISTHEP(IHEP)=200
913
914 JDAHEP(1,IHEP)=MHEP
915
916 JDAHEP(2,IHEP)=MHEP
917
918 IDHW(MHEP)=IDM
919
920 IDHEP(MHEP)=IDPDG(IDM)
921
922 JMOHEP(1,MHEP)=IHEP
923
924 JMOHEP(2,MHEP)=JMOHEP(2,IHEP)
925
926 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,MHEP))
927
928 CALL HWVEQU(4,VERTX,VHEP(1,MHEP))
929
930 NHEP=NHEP+1
931
932 ELSE
933
934 MHEP=IHEP
935
936 IDM=ID
937
938 ENDIF
939
940 ENDIF
941
942C Use CLEO/EURODEC packages for b-hadrons if requested
943
944 IF ((IDM.GE.221.AND.IDM.LE.231).OR.
945
946 & (IDM.GE.245.AND.IDM.LE.254)) THEN
947
948 IF (BDECAY.EQ.'CLEO') THEN
949
950 CALL HWDCLE(MHEP)
951
952 GOTO 100
953
954 ELSEIF (BDECAY.EQ.'EURO') THEN
955
956 CALL HWDEUR(MHEP)
957
958 GOTO 100
959
960 ENDIF
961
962 ENDIF
963
964C Choose decay mode
965
966 ISTHEP(MHEP)=ISTHEP(MHEP)+5
967
968 RN=HWR()
969
970 BF=0.
971
972 IM=LSTRT(IDM)
973
974 DO 10 I=1,NMODES(IDM)
975
976 BF=BF+BRFRAC(IM)
977
978 IF (BF.GE.RN) GOTO 20
979
980 10 IM=LNEXT(IM)
981
982 CALL HWWARN('HWDHAD',50,*20)
983
984 20 IF ((IDKPRD(1,IM).GE.1.AND.IDKPRD(1,IM).LE.13).OR.
985
986 & (IDKPRD(3,IM).GE.1.AND.IDKPRD(3,IM).LE.13)) THEN
987
988C Partonic decay of a heavy-(b,c)-hadron, store details
989
990 NQDK=NQDK+1
991
992 IF (NQDK.GT.NMXQDK) CALL HWWARN('HWDHAD',102,*999)
993
994 LOCQ(NQDK)=MHEP
995
996 IMQDK(NQDK)=IM
997
998 CALL HWVEQU(4,VERTX,VTXQDK(1,NQDK))
999
1000 GOTO 100
1001
1002 ELSE
1003
1004C Exclusive decay, add decay products to event record
1005
1006 IF (NHEP+NPRODS(IM).GT.NMXHEP)
1007
1008 & CALL HWWARN('HWDHAD',103,*999)
1009
1010 JDAHEP(1,MHEP)=NHEP+1
1011
1012 DO 30 I=1,NPRODS(IM)
1013
1014 NHEP=NHEP+1
1015
1016 IDHW(NHEP)=IDKPRD(I,IM)
1017
1018 IDHEP(NHEP)=IDPDG(IDKPRD(I,IM))
1019
1020 ISTHEP(NHEP)=193
1021
1022 JMOHEP(1,NHEP)=MHEP
1023
1024 JMOHEP(2,NHEP)=JMOHEP(2,MHEP)
1025
1026 PHEP(5,NHEP)=RMASS(IDKPRD(I,IM))
1027
1028 30 CALL HWVEQU(4,VERTX,VHEP(1,NHEP))
1029
1030 JDAHEP(2,MHEP)=NHEP
1031
1032 ENDIF
1033
1034C Next choose momenta:
1035
1036 IF (NPRODS(IM).EQ.1) THEN
1037
1038C 1-body decay: K0(BR) --> K0S,K0L
1039
1040 CALL HWVEQU(4,PHEP(1,MHEP),PHEP(1,NHEP))
1041
1042 ELSEIF (NPRODS(IM).EQ.2) THEN
1043
1044C 2-body decay
1045
1046C---SPECIAL TREATMENT OF POLARIZED MESONS
1047
1048 COSANG=TWO
1049
1050 IF (ID.EQ.IDHW(JMOHEP(1,MHEP))) THEN
1051
1052 MO=JMOHEP(1,MHEP)
1053
1054 RSUM=0
1055
1056 DO 40 I=1,3
1057
1058 40 RSUM=RSUM+RHOHEP(I,MO)
1059
1060 IF (RSUM.GT.ZERO) THEN
1061
1062 RSUM=RSUM*HWR()
1063
1064 IF (RSUM.LT.RHOHEP(1,MO)) THEN
1065
1066C---(1+COSANG)**2
1067
1068 COSANG=MAX(HWR(),HWR(),HWR())*TWO-ONE
1069
1070 ELSEIF (RSUM.LT.RHOHEP(1,MO)+RHOHEP(2,MO)) THEN
1071
1072C---1-COSANG**2
1073
1074 COSANG=2*COS((ACOS(HWR()*TWO-ONE)+PIFAC)/THREE)
1075
1076 ELSE
1077
1078C---(1-COSANG)**2
1079
1080 COSANG=MIN(HWR(),HWR(),HWR())*TWO-ONE
1081
1082 ENDIF
1083
1084 ENDIF
1085
1086 ENDIF
1087
1088 CALL HWDTWO(PHEP(1,MHEP),PHEP(1,NHEP-1),
1089
1090 & PHEP(1,NHEP),CMMOM(IM),COSANG,.FALSE.)
1091
1092 ELSEIF (NPRODS(IM).EQ.3) THEN
1093
1094C 3-body decay
1095
1096 IF (NME(IM).EQ.100) THEN
1097
1098C Use free massless (V-A)*(V-A) Matrix Element
1099
1100 CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-1),PHEP(1,NHEP-2),
1101
1102 & PHEP(1,NHEP),HWDWWT)
1103
1104 ELSEIF (NME(IM).EQ.101) THEN
1105
1106C Use bound massless (V-A)*(V-A) Matrix Element
1107
1108 WTMX=((PHEP(5,MHEP)-PHEP(5,NHEP))
1109
1110 & *(PHEP(5,MHEP)+PHEP(5,NHEP))
1111
1112 & +(PHEP(5,NHEP-1)-PHEP(5,NHEP-2))
1113
1114 & *(PHEP(5,NHEP-1)+PHEP(5,NHEP-2)))/TWO
1115
1116 WTMX2=WTMX**2
1117
1118 IPDG=ABS(IDHEP(MHEP))
1119
1120 XS=ONE-MAX(RMASS(MOD(IPDG/1000,10)),
1121
1122 & RMASS(MOD(IPDG/100,10)),RMASS(MOD(IPDG/10,10)))
1123
1124 & /(RMASS(MOD(IPDG/1000,10))+RMASS(MOD(IPDG/100,10))
1125
1126 & +RMASS(MOD(IPDG/10,10)))
1127
1128 50 CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-1),PHEP(1,NHEP-2),
1129
1130 & PHEP(1,NHEP),HWDWWT)
1131
1132 DOT1=HWULDO(PHEP(1,MHEP),PHEP(1,NHEP-1))
1133
1134 DOT2=HWULDO(PHEP(1,MHEP),PHEP(1,NHEP-2))
1135
1136 IF (DOT1*(WTMX-DOT1-XS*DOT2).LT.HWR()*WTMX2) GOTO 50
1137
1138 ELSE
1139
1140 CALL HWDTHR(PHEP(1,MHEP),PHEP(1,NHEP-2),PHEP(1,NHEP-1),
1141
1142 & PHEP(1,NHEP),HWDPWT)
1143
1144 ENDIF
1145
1146 ELSEIF (NPRODS(IM).EQ.4) THEN
1147
1148C 4-body decay
1149
1150 CALL HWDFOR(PHEP(1,MHEP ),PHEP(1,NHEP-3),PHEP(1,NHEP-2),
1151
1152 & PHEP(1,NHEP-1),PHEP(1,NHEP))
1153
1154 ELSEIF (NPRODS(IM).EQ.5) THEN
1155
1156C 5-body decay
1157
1158 CALL HWDFIV(PHEP(1,MHEP ),PHEP(1,NHEP-4),PHEP(1,NHEP-3),
1159
1160 & PHEP(1,NHEP-2),PHEP(1,NHEP-1),PHEP(1,NHEP))
1161
1162 ELSE
1163
1164 CALL HWWARN('HWDHAD',104,*999)
1165
1166 ENDIF
1167
1168 ENDIF
1169
1170 ENDIF
1171
1172 100 CONTINUE
1173
1174C---MAY HAVE OVERFLOWED /HEPEVT/
1175
1176 CALL HWWARN('HWDHAD',105,*999)
1177
1178 999 END
1179
1180CDECK ID>, HWDHGC.
1181
1182*CMZ :- -26/04/91 11.11.55 by Bryan Webber
1183
1184*-- Author : Mike Seymour
1185
1186C-----------------------------------------------------------------------
1187
1188 SUBROUTINE HWDHGC(TAU,FNREAL,FNIMAG)
1189
1190C-----------------------------------------------------------------------
1191
1192C CALCULATE THE COMPLEX FUNCTION F OF HHG eq 2.18
1193
1194C FOR USE IN H-->GAMMGAMM DECAYS
1195
1196C-----------------------------------------------------------------------
1197
1198 INCLUDE 'HERWIG61.INC'
1199
1200 DOUBLE PRECISION TAU,FNREAL,FNIMAG,FNLOG,FNSQR
1201
1202 IF (TAU.GT.ONE) THEN
1203
1204 FNREAL=(ASIN(1/SQRT(TAU)))**2
1205
1206 FNIMAG=0
1207
1208 ELSEIF (TAU.LT.ONE) THEN
1209
1210 FNSQR=SQRT(1-TAU)
1211
1212 FNLOG=LOG((1+FNSQR)/(1-FNSQR))
1213
1214 FNREAL=-0.25 * (FNLOG**2 - PIFAC**2)
1215
1216 FNIMAG= 0.5 * PIFAC*FNLOG
1217
1218 ELSE
1219
1220 FNREAL=0.25*PIFAC**2
1221
1222 FNIMAG=0
1223
1224 ENDIF
1225
1226 END