]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/src/hwdhig.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding rule violations corrected.
[u/mrichter/AliRoot.git] / HERWIG / src / hwdhig.f
CommitLineData
3820ca8e 1
2CDECK ID>, HWDHIG.
3
4*CMZ :- -24/04/92 14.23.44 by Mike Seymour
5
6*-- Author : Mike Seymour
7
8C-----------------------------------------------------------------------
9
10 SUBROUTINE HWDHIG(GAMINP)
11
12C-----------------------------------------------------------------------
13
14C HIGGS DECAY ROUTINE
15
16C A) FOR GAMinp=0 FIND AND DECAY HIGGS
17
18C B) FOR GAMinp>0 CALCULATE TOTAL HIGGS WIDTH
19
20C FOR EMH=GAMINP. STORE RESULT IN GAMINP.
21
22C-----------------------------------------------------------------------
23
24 INCLUDE 'HERWIG61.INC'
25
26 DOUBLE PRECISION HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,GAMINP,EMH,
27
28 & EMF,COLFAC,ENF,K1,K0,BET0,BET1,GAM0,GAM1,SCLOG,CFAC,XF,EM,GAMLIM,
29
30 & GAM,XW,EMW,XZ,EMZ,YW,YZ,EMI,TAUT,TAUW,WIDHIG,VECDEC,EMB,GAMB,
31
32 & TMIN,TMAX1,EM1,TMAX2,EM2,X1,X2,PROB,PCM,SUMR,SUMI,TAUTR,TAUTI,
33
34 & TAUWR,TAUWI,GFACTR
35
36 INTEGER HWRINT,IHIG,I,IFERM,NLOOK,I1,I2,IPART,IMODE,IDEC,MMAX
37
38 LOGICAL HWRLOG
39
40 EXTERNAL HWDHGF,HWR,HWRUNI,HWUSQR,HWUPCM,HWRINT,HWRLOG
41
42 SAVE GAM,EM,VECDEC
43
44 PARAMETER (NLOOK=100)
45
46 DIMENSION VECDEC(2,0:NLOOK)
47
48 EQUIVALENCE (EMW,RMASS(198)),(EMZ,RMASS(200))
49
50 DATA GAMLIM,GAM,EM/10D0,2*0D0/
51
52C---IF DECAY, FIND HIGGS (HWDHAD WILL HAVE GIVEN IT STATUS=1)
53
54 IF (GAMINP.EQ.ZERO) THEN
55
56 IHIG=0
57
58 DO 10 I=1,NHEP
59
60 10 IF (IHIG.EQ.0.AND.IDHW(I).EQ.201.AND.ISTHEP(I).EQ.1) IHIG=I
61
62 IF (IHIG.EQ.0) CALL HWWARN('HWDHIG',101,*999)
63
64 EMH=PHEP(5,IHIG)
65
66 IF (EMH.LE.ZERO) CALL HWWARN('HWDHIG',102,*999)
67
68 EMSCA=EMH
69
70 ELSE
71
72 EMH=GAMINP
73
74 IF (EMH.LE.ZERO) THEN
75
76 GAMINP=0
77
78 RETURN
79
80 ENDIF
81
82 ENDIF
83
84C---CALCULATE BRANCHING FRACTIONS
85
86C---FERMIONS
87
88C---NLL CORRECTION TO QUARK DECAY RATE (HHG eq 2.6-9)
89
90 ENF=0
91
92 DO 1 I=1,6
93
94 1 IF (2*RMASS(I).LT.EMH) ENF=ENF+1
95
96 K1=5/PIFAC**2
97
98 K0=3/(4*PIFAC**2)
99
100 BET0=(11*CAFAC-2*ENF)/3
101
102 BET1=(34*CAFAC**2-(10*CAFAC+6*CFFAC)*ENF)/3
103
104 GAM0=-8
105
106 GAM1=-404./3+40*ENF/9
107
108 SCLOG=LOG(EMH**2/QCDLAM**2)
109
110 CFAC=1 + ( K1/K0 - 2*GAM0 + GAM0*BET1/BET0**2*LOG(SCLOG)
111
112 & + (GAM0*BET1-GAM1*BET0)/BET0**2) / (BET0*SCLOG)
113
114 DO 100 IFERM=1,9
115
116 IF (IFERM.LE.6) THEN
117
118 EMF=RMASS(IFERM)
119
120 XF=(EMF/EMH)**2
121
122 COLFAC=FLOAT(NCOLO)
123
124 IF (EMF.GT.QCDLAM)
125
126 & EMF=EMF*(LOG(EMH/QCDLAM)/LOG(EMF/QCDLAM))**(GAM0/(2*BET0))
127
128 ELSE
129
130 EMF=RMASS(107+IFERM*2)
131
132 XF=(EMF/EMH)**2
133
134 COLFAC=1
135
136 CFAC=1
137
138 ENDIF
139
140 IF (FOUR*XF.LT.ONE) THEN
141
142 GFACTR=ALPHEM/(8.*SWEIN*EMW**2)
143
144 BRHIG(IFERM)=COLFAC*GFACTR*EMH*EMF**2 * (1-4*XF)**1.5 * CFAC
145
146 ELSE
147
148 BRHIG(IFERM)=0
149
150 ENDIF
151
152 100 CONTINUE
153
154C---W*W*/Z*Z*
155
156 IF (ABS(EM-EMH).GE.GAMLIM*GAM) THEN
157
158C---OFF EDGE OF LOOK-UP TABLE
159
160 XW=(EMW/EMH)**2
161
162 XZ=(EMZ/EMH)**2
163
164 YW=EMW*GAMW/EMH**2
165
166 YZ=EMZ*GAMZ/EMH**2
167
168 BRHIG(10)=.50*GFACTR * EMH**3 * HWDHGF(XW,YW)
169
170 BRHIG(11)=.25*GFACTR * EMH**3 * HWDHGF(XZ,YZ)
171
172 ELSE
173
174C---LOOK IT UP
175
176 EMI=((EMH-EM)/(GAM*GAMLIM)+1)*NLOOK/2.0
177
178 I1=INT(EMI)
179
180 I2=INT(EMI+1)
181
182 BRHIG(10)=.50*GFACTR * EMH**3 * ( VECDEC(1,I1)*(I2-EMI) +
183
184 & VECDEC(1,I2)*(EMI-I1) )
185
186 BRHIG(11)=.25*GFACTR * EMH**3 * ( VECDEC(2,I1)*(I2-EMI) +
187
188 & VECDEC(2,I2)*(EMI-I1) )
189
190 ENDIF
191
192C---GAMMAGAMMA
193
194 TAUT=(2*RMASS(6)/EMH)**2
195
196 TAUW=(2*EMW/EMH)**2
197
198 CALL HWDHGC(TAUT,TAUTR,TAUTI)
199
200 CALL HWDHGC(TAUW,TAUWR,TAUWI)
201
202 SUMR=4./3*( - 2*TAUT*( 1 + (1-TAUT)*TAUTR ) ) * ENHANC(6)
203
204 & +(2 + 3*TAUW*( 1 + (2-TAUW)*TAUWR ) ) * ENHANC(10)
205
206 SUMI=4./3*( - 2*TAUT*( (1-TAUT)*TAUTI ) ) * ENHANC(6)
207
208 & +( 3*TAUW*( (2-TAUW)*TAUWI ) ) * ENHANC(10)
209
210 BRHIG(12)=GFACTR*.03125*(ALPHEM/PIFAC)**2
211
212 & *EMH**3 * (SUMR**2 + SUMI**2)
213
214 WIDHIG=0
215
216 DO 200 IPART=1, 12
217
218 IF (IPART.LT.12) BRHIG(IPART)=BRHIG(IPART)*ENHANC(IPART)**2
219
220 200 WIDHIG=WIDHIG+BRHIG(IPART)
221
222 IF (WIDHIG.EQ.ZERO) CALL HWWARN('HWDHIG',103,*999)
223
224 DO 300 IPART=1, 12
225
226 300 BRHIG(IPART)=BRHIG(IPART)/WIDHIG
227
228 IF (EM.NE.RMASS(201)) THEN
229
230C---SET UP W*W*/Z*Z* LOOKUP TABLES
231
232 EM=EMH
233
234 GAM=WIDHIG
235
236 GAMLIM=MAX(GAMLIM,GAMMAX)
237
238 DO 400 I=0,NLOOK
239
240 EMH=(I*2.0/NLOOK-1)*GAM*GAMLIM+EM
241
242 XW=(EMW/EMH)**2
243
244 XZ=(EMZ/EMH)**2
245
246 YW=EMW*GAMW/EMH**2
247
248 YZ=EMZ*GAMZ/EMH**2
249
250 VECDEC(1,I)=HWDHGF(XW,YW)
251
252 VECDEC(2,I)=HWDHGF(XZ,YZ)
253
254 400 CONTINUE
255
256 EMH=EM
257
258 ENDIF
259
260 IF (GAMINP.GT.ZERO) THEN
261
262 GAMINP=WIDHIG
263
264 RETURN
265
266 ENDIF
267
268C---SEE IF USER SPECIFIED A DECAY MODE
269
270 IMODE=MOD(IPROC,100)
271
272C---IF NOT, CHOOSE ONE
273
274 IF (IMODE.LT.1.OR.IMODE.GT.12) THEN
275
276 MMAX=12
277
278 IF (IMODE.LT.1) MMAX=6
279
280 500 IMODE=HWRINT(1,MMAX)
281
282 IF (BRHIG(IMODE).LT.HWR()) GOTO 500
283
284 ENDIF
285
286C---SEE IF SPECIFIED DECAY IS POSSIBLE
287
288 IF (BRHIG(IMODE).EQ.ZERO) CALL HWWARN('HWDHIG',104,*999)
289
290 IF (IMODE.LE.6) THEN
291
292 IDEC=IMODE
293
294 ELSEIF (IMODE.LE.9) THEN
295
296 IDEC=107+IMODE*2
297
298 ELSEIF (IMODE.EQ.10) THEN
299
300 IDEC=198
301
302 ELSEIF (IMODE.EQ.11) THEN
303
304 IDEC=200
305
306 ELSEIF (IMODE.EQ.12) THEN
307
308 IDEC=59
309
310 ENDIF
311
312C---STATUS, IDs AND POINTERS
313
314 ISTHEP(IHIG)=195
315
316 DO 600 I=1,2
317
318 ISTHEP(NHEP+I)=193
319
320 IDHW(NHEP+I)=IDEC
321
322 IDHEP(NHEP+I)=IDPDG(IDEC)
323
324 JDAHEP(I,IHIG)=NHEP+I
325
326 JMOHEP(1,NHEP+I)=IHIG
327
328 JMOHEP(2,NHEP+I)=NHEP+(3-I)
329
330 JDAHEP(2,NHEP+I)=NHEP+(3-I)
331
332 PHEP(5,NHEP+I)=RMASS(IDEC)
333
334 IDEC=IDEC+6
335
336 IF (IDEC.EQ.204) IDEC=199
337
338 IF (IDEC.EQ.206) IDEC=200
339
340 IF (IDEC.EQ. 65) IDEC= 59
341
342 600 CONTINUE
343
344C---ALLOW W/Z TO BE OFF-SHELL
345
346 IF (IMODE.EQ.10.OR.IMODE.EQ.11) THEN
347
348 IF (IMODE.EQ.10) THEN
349
350 EMB=EMW
351
352 GAMB=GAMW
353
354 ELSE
355
356 EMB=EMZ
357
358 GAMB=GAMZ
359
360 ENDIF
361
362C---STANDARD MASS DISTRIBUTION
363
364 700 TMIN=ATAN(-EMB/GAMB)
365
366 TMAX1=ATAN((EMH**2/EMB-EMB)/GAMB)
367
368 EM1=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX1))+EMB))
369
370 TMAX2=ATAN(((EMH-EM1)**2/EMB-EMB)/GAMB)
371
372 EM2=HWUSQR(EMB*(GAMB*TAN(HWRUNI(0,TMIN,TMAX2))+EMB))
373
374 X1=(EM1/EMH)**2
375
376 X2=(EM2/EMH)**2
377
378C---CORRECT MASS DISTRIBUTION
379
380 PROB=HWUSQR(1+X1**2+X2**2-2*X1-2*X2-2*X1*X2)
381
382 & * ((X1+X2-1)**2 + 8*X1*X2)
383
384 IF (.NOT.HWRLOG(PROB)) GOTO 700
385
386C---CALCULATE SPIN DENSITY MATRIX
387
388 RHOHEP(1,NHEP+1)=4*X1*X2 / (8*X1*X2 + (X1+X2-1)**2)
389
390 RHOHEP(2,NHEP+1)=(X1+X2-1)**2 / (8*X1*X2 + (X1+X2-1)**2)
391
392 RHOHEP(3,NHEP+1)=RHOHEP(1,NHEP+1)
393
394C---SYMMETRIZE DISTRIBUTIONS IN PARTICLES 1,2
395
396 IF (HWRLOG(HALF)) THEN
397
398 PHEP(5,NHEP+1)=EM1
399
400 PHEP(5,NHEP+2)=EM2
401
402 ELSE
403
404 PHEP(5,NHEP+1)=EM2
405
406 PHEP(5,NHEP+2)=EM1
407
408 ENDIF
409
410 ENDIF
411
412C---DO DECAY
413
414 PCM=HWUPCM(EMH,PHEP(5,NHEP+1),PHEP(5,NHEP+2))
415
416 IF (PCM.LT.ZERO) CALL HWWARN('HWDHIG',105,*999)
417
418 CALL HWDTWO(PHEP(1,IHIG),PHEP(1,NHEP+1),PHEP(1,NHEP+2),
419
420 & PCM,TWO,.TRUE.)
421
422 NHEP=NHEP+2
423
424C---IF QUARK DECAY, HADRONIZE
425
426 IF (IMODE.LE.6) THEN
427
428 ISTHEP(NHEP-1)=113
429
430 ISTHEP(NHEP)=114
431
432 CALL HWBGEN
433
434 CALL HWDHOB
435
436 CALL HWCFOR
437
438 CALL HWCDEC
439
440 ENDIF
441
442 999 END
443
444CDECK ID>, HWDHOB.
445
446*CMZ :- -20/10/99 09:46:43 by Peter Richardson
447
448*-- Author : Ian Knowles & Bryan Webber
449
450C-----------------------------------------------------------------------
451
452 SUBROUTINE HWDHOB
453
454C-----------------------------------------------------------------------
455
456C Performs decays of heavy objects (heavy quarks & SUSY particles)
457
458C MODIFIED TO INCLUDE R-PARITY VIOLATING SUSY PR 9/4/99
459
460C-----------------------------------------------------------------------
461
462 INCLUDE 'HERWIG61.INC'
463
464 DOUBLE PRECISION HWUMBW,HWUPCM,HWR,SDKM,RN,BF,PCM,
465
466 & EMMX,EMWSQ,GMWSQ,EMLIM,PW(5),EMTST,HWDPWT,HWDWWT,HWULDO,PDW(5,3)
467
468 INTEGER IST(3),IHEP,IS,ID,IM,I,JHEP,KHEP,LHEP,MHEP,NPR,ISM,JCM,
469
470 & MTRY,NTRY,IDM,IDM2,THEP,CLSAVE(2),WHEP,RHEP
471
472 LOGICAL FOUND
473
474 EXTERNAL HWR,HWDPWT,HWDWWT
475
476 DATA IST/113,114,114/
477
478 IF (IERROR.NE.0) RETURN
479
480 10 FOUND=.FALSE.
481
482 CLSAVE(1) = 0
483
484 CLSAVE(2) = 0
485
486 DO 60 IHEP=1,NMXHEP
487
488 IS=ISTHEP(IHEP)
489
490 ID=IDHW(IHEP)
491
492 IF (.NOT.RSTAB(ID).AND.(ID.EQ.6.OR.ID.EQ.12.OR.
493
494 & (ID.GE.203.AND.ID.LE.218).OR.ABS(IDPDG(ID)).GT.1000000).AND.
495
496 & (IS.EQ.190.OR.(IS.GE.147.AND.IS.LE.151))) THEN
497
498 FOUND=.TRUE.
499
500 IF(.NOT.RPARTY) THEN
501
502 NHEP = NHEP+1
503
504 ISTHEP(NHEP) = 3
505
506 IDHW(NHEP) = 20
507
508 IDHEP(NHEP) = 0
509
510 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
511
512 CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
513
514 JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
515
516 JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
517
518 JDAHEP(1,NHEP)=JDAHEP(1,IHEP)
519
520 JDAHEP(2,NHEP)=JDAHEP(2,IHEP)
521
522 ENDIF
523
524C Make a copy of decaying object
525
526 NHEP=NHEP+1
527
528 ISTHEP(NHEP)=155
529
530 IDHW(NHEP)=IDHW(IHEP)
531
532 IDHEP(NHEP)=IDHEP(IHEP)
533
534 CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
535
536 CALL HWVEQU(4,VHEP(1,IHEP),VHEP(1,NHEP))
537
538 JMOHEP(1,NHEP)=JMOHEP(1,IHEP)
539
540 JMOHEP(2,NHEP)=JMOHEP(2,IHEP)
541
542 MTRY=0
543
544 15 MTRY=MTRY+1
545
546C Select decay mode
547
548 RN=HWR()
549
550 BF=0.
551
552 IM=LSTRT(ID)
553
554 DO 20 I=1,NMODES(ID)
555
556 BF=BF+BRFRAC(IM)
557
558 IF (BF.GE.RN) GOTO 30
559
560 20 IM=LNEXT(IM)
561
562 CALL HWWARN('HWDHOB',50,*30)
563
564 30 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHOB',100,*999)
565
566 NPR=NPRODS(IM)
567
568 JDAHEP(1,NHEP)=NHEP+1
569
570 JDAHEP(2,NHEP)=NHEP+NPR
571
572C Reset colour pointers (if set)
573
574 JHEP=JMOHEP(2,IHEP)
575
576 IF (JHEP.GT.0) THEN
577
578 IF (JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
579
580 IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
581
582 & .AND.ABS(IDHEP(JHEP)).GT.1000000
583
584 & .AND.JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1) = NHEP
585
586 ENDIF
587
588 JHEP=JDAHEP(2,IHEP)
589
590 IF (JHEP.GT.0) THEN
591
592 IF (JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
593
594 IF(.NOT.RPARTY.AND.ISTHEP(JHEP).EQ.155
595
596 & .AND.ABS(IDHEP(JHEP)).GT.1000000
597
598 & .AND.JMOHEP(2,JHEP-1).EQ.IHEP) JMOHEP(2,JHEP-1) = NHEP
599
600 ENDIF
601
602C--Reset colour pointers if baryon number violated
603
604 IF(.NOT.RPARTY) THEN
605
606 DO JHEP=1,NHEP
607
608 IF(ISTHEP(JHEP).EQ.155
609
610 & .AND.ABS(IDHEP(JHEP)).GT.1000000.AND.
611
612 & JDAHEP(2,JHEP-1).EQ.IHEP) JDAHEP(2,JHEP-1)= NHEP
613
614 IF(JDAHEP(2,JHEP).EQ.IHEP) JDAHEP(2,JHEP)=NHEP
615
616 IF(JMOHEP(2,JHEP).EQ.IHEP) JMOHEP(2,JHEP)=NHEP
617
618 ENDDO
619
620 IF(HRDCOL(1,1).EQ.IHEP) HRDCOL(1,1)=NHEP
621
622 ENDIF
623
624C Relabel original track
625
626 ISTHEP(IHEP)=3
627
628 JMOHEP(2,IHEP)=JMOHEP(1,IHEP)
629
630 JDAHEP(1,IHEP)=NHEP
631
632 JDAHEP(2,IHEP)=NHEP
633
634C Label decay products and choose masses
635
636 LHEP=NHEP
637
638 MHEP=LHEP+1
639
640 NTRY=0
641
642 35 NTRY=NTRY+1
643
644 SDKM=PHEP(5,NHEP)
645
646 DO 40 I=1,NPR
647
648 NHEP=NHEP+1
649
650 IDHW(NHEP)=IDKPRD(I,IM)
651
652 IDHEP(NHEP)=IDPDG(IDKPRD(I,IM))
653
654 ISTHEP(NHEP)=IST(I)
655
656 JMOHEP(1,NHEP)=LHEP
657
658 JDAHEP(1,NHEP)=0
659
660 PHEP(5,NHEP)=HWUMBW(IDKPRD(I,IM))
661
662 40 SDKM=SDKM-PHEP(5,NHEP)
663
664 IF (SDKM.LT.ZERO) THEN
665
666 NHEP=NHEP-NPR
667
668 IF (NTRY.LE.NETRY) GO TO 35
669
670 CALL HWWARN('HWDHOB',1,*45)
671
672 45 IF (MTRY.LE.NETRY) GO TO 15
673
674 CALL HWWARN('HWDHOB',101,*999)
675
676 ENDIF
677
678C Assign production vertices to decay products
679
680 CALL HWUDKL(ID,PHEP(1,IHEP),VHEP(1,MHEP))
681
682 CALL HWVSUM(4,VHEP(1,IHEP),VHEP(1,MHEP),VHEP(1,MHEP))
683
684 CALL HWVEQU(4,VHEP(1,MHEP),VHEP(1,NHEP))
685
686 IF (NPR.EQ.2) THEN
687
688C Two body decay: LHEP -> MHEP + NHEP
689
690 PCM=HWUPCM(PHEP(5,IHEP),PHEP(5,MHEP),PHEP(5,NHEP))
691
692 CALL HWDTWO(PHEP(1,IHEP),PHEP(1,MHEP),
693
694 & PHEP(1,NHEP),PCM,TWO,.FALSE.)
695
696 ELSEIF (NPR.EQ.3) THEN
697
698C Three body decay: LHEP -> KHEP + MHEP + NHEP
699
700 KHEP=MHEP
701
702 MHEP=MHEP+1
703
704C Provisional colour self-connection of KHEP
705
706 JMOHEP(2,KHEP)=KHEP
707
708 JDAHEP(2,KHEP)=KHEP
709
710 IF (NME(IM).EQ.100) THEN
711
712C Generate decay momenta using full (V-A)*(V-A) matrix element
713
714 EMMX=PHEP(5,IHEP)-PHEP(5,NHEP)
715
716 EMWSQ=RMASS(198)**2
717
718 GMWSQ=(RMASS(198)*GAMW)**2
719
720 EMLIM=GMWSQ
721
722 IF (EMMX.LT.RMASS(198)) EMLIM=EMLIM+(EMWSQ-EMMX**2)**2
723
724 50 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
725
726 & PHEP(1,KHEP),PHEP(1,NHEP),HWDWWT)
727
728 CALL HWVSUM(4,PHEP(1,KHEP),PHEP(1,MHEP),PW)
729
730 PW(5)=HWULDO(PW,PW)
731
732 EMTST=(EMWSQ-PW(5))**2
733
734 IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 50
735
736 PW(5)=SQRT(PW(5))
737
738C Assign production vertices to 1 and 2
739
740 CALL HWUDKL(198,PW,VHEP(1,KHEP))
741
742 CALL HWVSUM(4,VHEP(1,NHEP),VHEP(1,KHEP),VHEP(1,KHEP))
743
744 ELSEIF(NME(IM).EQ.300) THEN
745
746C Generate momenta using 3-body RPV matrix element
747
748 CALL HWDRME(LHEP,KHEP)
749
750 ELSE
751
752C Three body phase space decay
753
754 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,MHEP),
755
756 & PHEP(1,KHEP),PHEP(1,NHEP),HWDPWT)
757
758 ENDIF
759
760 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
761
762 ELSEIF(NPR.EQ.4) THEN
763
764C Four body decay: LHEP -> KHEP + RHEP + MHEP + NHEP
765
766 KHEP = MHEP
767
768 RHEP = MHEP+1
769
770 MHEP = MHEP+2
771
772C Provisional colour connections of KHEP and RHEP
773
774 JMOHEP(2,KHEP)=RHEP
775
776 JDAHEP(2,KHEP)=RHEP
777
778 JMOHEP(2,RHEP)=KHEP
779
780 JDAHEP(2,RHEP)=KHEP
781
782C Four body phase space decay
783
784 CALL HWDFOR(PHEP(1,IHEP),PHEP(1,KHEP),PHEP(1,RHEP),
785
786 & PHEP(1,MHEP),PHEP(1,NHEP))
787
788 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,RHEP))
789
790 CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,MHEP))
791
792 ELSE
793
794 CALL HWWARN('HWDHOB',102,*999)
795
796 ENDIF
797
798C Colour connections
799
800 IF (ID.EQ.6.OR.ID.EQ.12.OR.(ID.GE.209.AND.ID.LE.212)
801
802 & .OR.(ID.GE.215.AND.ID.LE.218)) THEN
803
804 IF (NPR.EQ.3.AND.NME(IM).EQ.100) THEN
805
806C usual heavy quark decay
807
808 JMOHEP(2,KHEP)=MHEP
809
810 JDAHEP(2,KHEP)=MHEP
811
812 JMOHEP(2,MHEP)=KHEP
813
814 JDAHEP(2,MHEP)=KHEP
815
816 JMOHEP(2,NHEP)=LHEP
817
818 JDAHEP(2,NHEP)=LHEP
819
820 ELSEIF (ABS(IDHEP(MHEP)).EQ.37) THEN
821
822C heavy quark to charged Higgs
823
824 JMOHEP(2,MHEP)=MHEP
825
826 JDAHEP(2,MHEP)=MHEP
827
828 JMOHEP(2,NHEP)=LHEP
829
830 JDAHEP(2,NHEP)=LHEP
831
832 ELSEIF (ABS(IDHEP(NHEP)).EQ.37) THEN
833
834 JMOHEP(2,MHEP)=LHEP
835
836 JDAHEP(2,MHEP)=LHEP
837
838 JMOHEP(2,NHEP)=NHEP
839
840 JDAHEP(2,NHEP)=NHEP
841
842 ELSE
843
844 CALL HWWARN('HWDHOB',103,*999)
845
846 ENDIF
847
848 ELSE
849
850 IF(.NOT.RPARTY.AND.
851
852 & ((NPR.EQ.2.AND.ID.GE.401.AND.ID.LT.448.AND.
853
854 & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132)
855
856 & .OR.(NPR.EQ.3.AND.ID.GE.449.AND.ID.LE.457.AND.
857
858 & IDHW(MHEP).LE.132.AND.IDHW(NHEP).LE.132.AND.
859
860 & IDHW(MHEP-1).LE.132))) THEN
861
862C R-parity violating SUSY decays
863
864 IF(NPR.EQ.2) THEN
865
866C--Rparity slepton colour connections
867
868 IF(ID.GE.425.AND.ID.LE.448) THEN
869
870 IF(IDHW(MHEP).GT.12) THEN
871
872 JMOHEP(2,MHEP) = MHEP
873
874 JDAHEP(2,MHEP) = MHEP
875
876 JMOHEP(2,NHEP) = NHEP
877
878 JDAHEP(2,NHEP) = NHEP
879
880 ELSE
881
882 JMOHEP(2,MHEP) = NHEP
883
884 JDAHEP(2,MHEP) = NHEP
885
886 JMOHEP(2,NHEP) = MHEP
887
888 JDAHEP(2,NHEP) = MHEP
889
890 ENDIF
891
892C--Rparity squark colour connections
893
894 ELSE
895
896 IF(IDHEP(LHEP).GT.0) THEN
897
898C--LQD decay colour connections
899
900 IF(IDHW(MHEP).GT.12) THEN
901
902 JMOHEP(2,MHEP) = MHEP
903
904 JDAHEP(2,MHEP) = MHEP
905
906 JMOHEP(2,NHEP) = LHEP
907
908 JDAHEP(2,NHEP) = LHEP
909
910 ELSE
911
912C--UDD decay colour connections
913
914 HVFCEN = .TRUE.
915
916 CALL HWDRCL(LHEP,MHEP,CLSAVE)
917
918 ENDIF
919
920 ELSE
921
922C--Antisquark connections
923
924 IF(IDHW(MHEP).GT.12) THEN
925
926 JMOHEP(2,MHEP) = MHEP
927
928 JDAHEP(2,MHEP) = MHEP
929
930 JMOHEP(2,NHEP) = LHEP
931
932 JDAHEP(2,NHEP) = LHEP
933
934 ELSE
935
936 HVFCEN = .TRUE.
937
938 CALL HWDRCL(LHEP,MHEP,CLSAVE)
939
940 ENDIF
941
942 ENDIF
943
944 ENDIF
945
946 ELSE
947
948 IF(ID.GE.450.AND.ID.LE.457) THEN
949
950C--Rparity Neutralino/Chargino colour connection
951
952 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
953
954 & AND.IDHW(NHEP).LE.12) THEN
955
956 HVFCEN = .TRUE.
957
958 CALL HWDRCL(LHEP,MHEP,CLSAVE)
959
960 ELSE
961
962 JMOHEP(2,MHEP) = NHEP
963
964 JDAHEP(2,MHEP) = NHEP
965
966 JMOHEP(2,NHEP) = MHEP
967
968 JDAHEP(2,NHEP) = MHEP
969
970 ENDIF
971
972C--Rparity gluino colour connections
973
974 ELSEIF(ID.EQ.449) THEN
975
976 IF(IDHW(MHEP-1).LE.12.AND.IDHW(MHEP).LE.12.
977
978 & AND.IDHW(NHEP).LE.12) THEN
979
980 HVFCEN = .TRUE.
981
982 CALL HWDRCL(LHEP,MHEP,CLSAVE)
983
984C--Now the lepton number violating decay
985
986 ELSE
987
988 IF(IDHW(MHEP).LE.6) THEN
989
990 JMOHEP(2,MHEP) = LHEP
991
992 JDAHEP(2,MHEP) = NHEP
993
994 JMOHEP(2,NHEP) = MHEP
995
996 JDAHEP(2,NHEP) = LHEP
997
998 ELSE
999
1000 JMOHEP(2,MHEP) = NHEP
1001
1002 JDAHEP(2,MHEP) = LHEP
1003
1004 JMOHEP(2,NHEP) = LHEP
1005
1006 JDAHEP(2,NHEP) = MHEP
1007
1008 ENDIF
1009
1010 ENDIF
1011
1012 ELSE
1013
1014 CALL HWWARN('HWDHOB',104,*999)
1015
1016 ENDIF
1017
1018 ENDIF
1019
1020 ELSE
1021
1022C Normal SUSY decays
1023
1024 IF (ID.LE.448.AND.ID.GT.207) THEN
1025
1026C Squark (or slepton)
1027
1028 IF (IDHW(MHEP).EQ.449) THEN
1029
1030 IF (IDHEP(LHEP).GT.0) THEN
1031
1032 JMOHEP(2,MHEP)=LHEP
1033
1034 JDAHEP(2,MHEP)=NHEP
1035
1036 JMOHEP(2,NHEP)=MHEP
1037
1038 JDAHEP(2,NHEP)=LHEP
1039
1040 ELSE
1041
1042 JMOHEP(2,MHEP)=NHEP
1043
1044 JDAHEP(2,MHEP)=LHEP
1045
1046 JMOHEP(2,NHEP)=LHEP
1047
1048 JDAHEP(2,NHEP)=MHEP
1049
1050 ENDIF
1051
1052 ELSE
1053
1054 IF(NPR.EQ.3.AND.IDHW(MHEP).LE.12) THEN
1055
1056 JMOHEP(2,MHEP)=NHEP
1057
1058 JDAHEP(2,MHEP)=NHEP
1059
1060 JMOHEP(2,NHEP)=MHEP
1061
1062 JDAHEP(2,NHEP)=MHEP
1063
1064 ELSE
1065
1066 JMOHEP(2,MHEP)=MHEP
1067
1068 JDAHEP(2,MHEP)=MHEP
1069
1070 JMOHEP(2,NHEP)=LHEP
1071
1072 JDAHEP(2,NHEP)=LHEP
1073
1074 ENDIF
1075
1076 ENDIF
1077
1078 ELSEIF (ID.EQ.449) THEN
1079
1080C Gluino
1081
1082 IF (IDHW(NHEP).EQ.13) THEN
1083
1084 JMOHEP(2,MHEP)=MHEP
1085
1086 JDAHEP(2,MHEP)=MHEP
1087
1088 JMOHEP(2,NHEP)=LHEP
1089
1090 JDAHEP(2,NHEP)=LHEP
1091
1092 ELSEIF (IDHEP(MHEP).GT.0) THEN
1093
1094 JMOHEP(2,MHEP)=LHEP
1095
1096 JDAHEP(2,MHEP)=NHEP
1097
1098 JMOHEP(2,NHEP)=MHEP
1099
1100 JDAHEP(2,NHEP)=LHEP
1101
1102 ELSE
1103
1104 JMOHEP(2,MHEP)=NHEP
1105
1106 JDAHEP(2,MHEP)=LHEP
1107
1108 JMOHEP(2,NHEP)=LHEP
1109
1110 JDAHEP(2,NHEP)=MHEP
1111
1112 ENDIF
1113
1114 ELSE
1115
1116C Gaugino or Higgs
1117
1118 JMOHEP(2,MHEP)=NHEP
1119
1120 JDAHEP(2,MHEP)=NHEP
1121
1122 JMOHEP(2,NHEP)=MHEP
1123
1124 JDAHEP(2,NHEP)=MHEP
1125
1126 ENDIF
1127
1128 ENDIF
1129
1130 ENDIF
1131
1132C---SPECIAL CASE FOR THREE-BODY TOP DECAYS:
1133
1134C RELABEL THEM AS TWO TWO-BODY DECAYS FOR PARTON SHOWERING
1135
1136 IF ((ID.EQ.6.OR.ID.EQ.12).AND.NPR.EQ.3.AND.NME(IM).EQ.100) THEN
1137
1138C---STORE W DECAY PRODUCTS
1139
1140 CALL HWVEQU(10,PHEP(1,KHEP),PDW)
1141
1142C---BOOST THEM INTO W REST FRAME
1143
1144 CALL HWULOF(PW,PDW(1,1),PDW(1,3))
1145
1146C---REPLACE THEM BY W
1147
1148 CALL HWVEQU(5,PW,PHEP(1,KHEP))
1149
1150 WHEP=KHEP
1151
1152 IDHW(KHEP)=198
1153
1154 IF (ID.EQ.12) IDHW(KHEP)=199
1155
1156 IDHEP(KHEP)=IDPDG(IDHW(KHEP))
1157
1158 JMOHEP(2,KHEP)=KHEP
1159
1160 JDAHEP(2,KHEP)=KHEP
1161
1162 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,KHEP))
1163
1164C---AND MOVE B UP
1165
1166 CALL HWVEQU(5,PHEP(1,NHEP),PHEP(1,MHEP))
1167
1168 IDHW(MHEP)=IDHW(NHEP)
1169
1170 IDHEP(MHEP)=IDHEP(NHEP)
1171
1172 JDAHEP(2,LHEP)=MHEP
1173
1174 JMOHEP(2,MHEP)=JMOHEP(2,NHEP)
1175
1176 JDAHEP(2,MHEP)=JDAHEP(2,NHEP)
1177
1178 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,MHEP))
1179
1180 NHEP=MHEP
1181
1182C---DO PARTON SHOWER
1183
1184 EMSCA=PHEP(5,IHEP)
1185
1186 CALL HWBGEN
1187
1188 IF (IERROR.NE.0) RETURN
1189
1190C---FIND BOOSTED W MOMENTUM
1191
1192 NTRY=0
1193
1194 41 NTRY=NTRY+1
1195
1196 IF (NTRY.GT.NHEP.OR.WHEP.LE.0.OR.WHEP.GT.NHEP)
1197
1198 $ CALL HWWARN('HWDHOB',101,*999)
1199
1200 WHEP=JDAHEP(1,WHEP)
1201
1202 IF (ISTHEP(WHEP).NE.190) GOTO 41
1203
1204C---AND HENCE ITS CHILDRENS MOMENTA
1205
1206 CALL HWULOB(PHEP(1,WHEP),PDW(1,3),PHEP(1,NHEP+1))
1207
1208 CALL HWVDIF(4,PHEP(1,WHEP),PHEP(1,NHEP+1),PHEP(1,NHEP+2))
1209
1210 PHEP(5,NHEP+2)=PDW(5,2)
1211
1212C---LABEL THEM
1213
1214 ISTHEP(WHEP)=195
1215
1216 DO 51 I=1,2
1217
1218 IDHW(NHEP+I)=IDKPRD(I,IM)
1219
1220 IDHEP(NHEP+I)=IDPDG(IDHW(NHEP+I))
1221
1222 ISTHEP(NHEP+I)=112+I
1223
1224 JDAHEP(I,WHEP)=NHEP+I
1225
1226 JMOHEP(1,NHEP+I)=WHEP
1227
1228 JMOHEP(2,NHEP+I)=NHEP+3-I
1229
1230 JDAHEP(2,NHEP+I)=NHEP+3-I
1231
1232 51 CONTINUE
1233
1234 NHEP=NHEP+2
1235
1236C---ASSIGN PRODUCTION VERTICES TO 1 AND 2
1237
1238 CALL HWUDKL(198,PW,VHEP(1,NHEP))
1239
1240 CALL HWVSUM(4,VHEP(1,WHEP),VHEP(1,NHEP),VHEP(1,NHEP))
1241
1242 CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
1243
1244C---DO PARTON SHOWERS
1245
1246 EMSCA=PW(5)
1247
1248 CALL HWBGEN
1249
1250 IF (IERROR.NE.0) RETURN
1251
1252 ELSE
1253
1254C Do parton showers
1255
1256 EMSCA=PHEP(5,IHEP)
1257
1258 CALL HWBGEN
1259
1260 IF (IERROR.NE.0) RETURN
1261
1262 ENDIF
1263
1264 ENDIF
1265
1266C--New to correct colour connections in Rslash
1267
1268 IF(CLSAVE(1).NE.0) THEN
1269
1270 THEP = MHEP+1
1271
1272 ID = IDHW(CLSAVE(1))
1273
1274 IDM = IDHW(JMOHEP(1,CLSAVE(1)))
1275
1276 IDM2 = IDHW(LHEP)
1277
1278 IF(IDM.EQ.15) ID=IDHW(JMOHEP(1,JMOHEP(1,CLSAVE(1))))
1279
1280 IF((ID.LE.6.AND.((IDM.GE.419.AND.IDM.LE.424).OR.IDM.EQ.411.OR.
1281
1282 & IDM.EQ.412).
1283
1284 & AND.((IDM2.GE.413.AND.IDM2.LE.418)
1285
1286 & .OR.IDM2.EQ.449).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1287
1288 & .OR.(ID.LE.6.AND.IDM.EQ.449.AND.
1289
1290 & (((IDM2.GE.413.AND.IDM2.LE.418).OR.IDM2.EQ.405.OR.IDM2.EQ.406)
1291
1292 & .OR.IDM2.EQ.449)).OR.
1293
1294 & (IDM.EQ.15.AND.ID.LE.12.AND.ID.GE.7.AND.((IDM2.GE.413.AND.
1295
1296 & IDM2.LE.418).OR.IDM2.EQ.449.OR.IDM2.
1297
1298 & EQ.405.OR.IDM2.EQ.406))) THEN
1299
1300 IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1301
1302 IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1303
1304 & JMOHEP(2,CLSAVE(2)) = THEP
1305
1306 JDAHEP(2,MHEP) = CLSAVE(1)
1307
1308 JDAHEP(2,THEP) = CLSAVE(2)
1309
1310 ELSE
1311
1312 IF(IDHW(CLSAVE(1)).NE.13.AND.IDHW(CLSAVE(1)).NE.449)
1313
1314 & JMOHEP(2,CLSAVE(2)) = MHEP
1315
1316 JDAHEP(2,MHEP) = CLSAVE(2)
1317
1318 JDAHEP(2,THEP) = CLSAVE(1)
1319
1320 ENDIF
1321
1322 ELSEIF((ID.GT.6.AND.ID.LE.12.
1323
1324 & AND.((IDM.GE.413.AND.IDM.LE.418).OR.IDM.EQ.405.OR.
1325
1326 & IDM.EQ.406).AND.
1327
1328 & ((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1329
1330 & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1331
1332 & (ID.GT.6.AND.ID.LE.12.AND.IDM.EQ.449.
1333
1334 & AND.((IDM2.GE.419.AND.IDM2.LE.424).OR.IDM2.EQ.449.OR.
1335
1336 & IDM2.EQ.411.OR.IDM2.EQ.412)).OR.
1337
1338 & (IDM.EQ.15.AND.ID.LE.6.AND.((IDM2.GE.419.AND.
1339
1340 & IDM2.LE.424).OR.IDM2.EQ.449.OR.IDM2.EQ.411.OR.
1341
1342 & IDM2.EQ.412))) THEN
1343
1344 IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
1345
1346 JDAHEP(2,CLSAVE(2))=THEP
1347
1348 JMOHEP(2,MHEP)=CLSAVE(1)
1349
1350 JMOHEP(2,THEP)=CLSAVE(2)
1351
1352 ELSE
1353
1354 JDAHEP(2,CLSAVE(2))=MHEP
1355
1356 JMOHEP(2,MHEP)=CLSAVE(2)
1357
1358 JMOHEP(2,THEP)=CLSAVE(1)
1359
1360 ENDIF
1361
1362 ENDIF
1363
1364 COLUPD = .FALSE.
1365
1366 CALL HWBCON
1367
1368 ENDIF
1369
1370 IF (IHEP.EQ.NHEP) GOTO 70
1371
1372 60 CONTINUE
1373
1374 70 IF (FOUND) THEN
1375
1376C Fix any SUSY colour disconnections
1377
1378 DO 80 IHEP=1,NHEP
1379
1380 IF (ISTHEP(IHEP).GE.147.AND.ISTHEP(IHEP).LE.151
1381
1382 & .AND.JDAHEP(2,IHEP).EQ.0) THEN
1383
1384 IM=JMOHEP(1,IHEP)
1385
1386C Chase connection back through SUSY decays
1387
1388 75 IM=JMOHEP(1,IM)
1389
1390 ISM=ISTHEP(IM)
1391
1392 IF (ISM.EQ.120) GOTO 80
1393
1394 IF (ISM.NE.123.AND.ISM.NE.124.AND.ISM.NE.155) GOTO 75
1395
1396C Look for unclustered parton to connect
1397
1398 DO JHEP=1,NHEP
1399
1400 IF (ISTHEP(JHEP).GE.147.AND.ISTHEP(JHEP).LE.151) THEN
1401
1402 JCM=JMOHEP(2,JHEP)
1403
1404 IF (JCM.EQ.IM) THEN
1405
1406C Found it: connect
1407
1408 JMOHEP(2,JHEP)=IHEP
1409
1410 JDAHEP(2,IHEP)=JHEP
1411
1412 GOTO 80
1413
1414 ENDIF
1415
1416 ENDIF
1417
1418 ENDDO
1419
1420C Not found: need to go further back
1421
1422 GOTO 75
1423
1424 ENDIF
1425
1426 80 CONTINUE
1427
1428C Go back to check for further heavy decay products
1429
1430 GOTO 10
1431
1432 ENDIF
1433
1434 999 END
1435
1436CDECK ID>, HWDHVY.
1437
1438*CMZ :- -26/04/91 12.19.24 by Federico Carminati
1439
1440*-- Author : Ian Knowles & Bryan Webber
1441
1442C-----------------------------------------------------------------------
1443
1444 SUBROUTINE HWDHVY
1445
1446C-----------------------------------------------------------------------
1447
1448C Performs partonic decays of hadrons containing heavy quark(s):
1449
1450C either, meson/baryon spectator model weak decays;
1451
1452C or, quarkonia -> 2-gluons, q-qbar, 3-gluons, or 2-gluons + photon.
1453
1454C-----------------------------------------------------------------------
1455
1456 INCLUDE 'HERWIG61.INC'
1457
1458 DOUBLE PRECISION HWULDO,HWR,XS,XB,EMWSQ,GMWSQ,EMLIM,PW(4),
1459
1460 & EMTST,X1,X2,X3,TEST,HWDWWT,HWDPWT
1461
1462 INTEGER IST(3),I,IHEP,IM,ID,IDQ,IQ,IS,J
1463
1464 EXTERNAL HWR,HWDWWT,HWDPWT,HWULDO
1465
1466 DATA IST/113,114,114/
1467
1468 IF (IERROR.NE.0) RETURN
1469
1470 DO 100 I=1,NMXQDK
1471
1472 IF (I.GT.NQDK) THEN
1473
1474 NQDK=0
1475
1476 RETURN
1477
1478 ENDIF
1479
1480 IHEP=LOCQ(I)
1481
1482 IF (ISTHEP(IHEP).EQ.199) GOTO 100
1483
1484 IM=IMQDK(I)
1485
1486 IF (NHEP+NPRODS(IM).GT.NMXHEP) CALL HWWARN('HWDHVY',100,*999)
1487
1488 IF (IDKPRD(4,IM).NE.0) THEN
1489
1490C Weak decay of meson or baryon
1491
1492C Idenitify decaying heavy quark and spectator
1493
1494 ID=IDHW(IHEP)
1495
1496 IF (ID.EQ.136.OR.ID.EQ.140.OR.ID.EQ.144.OR.
1497
1498 & ID.EQ.150.OR.ID.EQ.155.OR.ID.EQ.158.OR.ID.EQ.161.OR.
1499
1500 & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.11)) THEN
1501
1502C c hadron or c decay of B_c+
1503
1504 IDQ=4
1505
1506 IQ=NHEP+1
1507
1508 IS=NHEP+2
1509
1510 ELSEIF (ID.EQ.171.OR.ID.EQ.175.OR.ID.EQ.179.OR.
1511
1512 & ID.EQ.185.OR.ID.EQ.190.OR.ID.EQ.194.OR.ID.EQ.196.OR.
1513
1514 & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.5)) THEN
1515
1516C cbar hadron or cbar decay of B_c-
1517
1518 IDQ=10
1519
1520 IS=NHEP+1
1521
1522 IQ=NHEP+2
1523
1524 ELSEIF ((ID.GE.221.AND.ID.LE.229).OR.
1525
1526 & (ID.EQ.230.AND.IDKPRD(4,IM).EQ.10)) THEN
1527
1528C b hadron or b decay of B_c-
1529
1530 IDQ=5
1531
1532 IQ=NHEP+1
1533
1534 IS=NHEP+2
1535
1536 ELSEIF ((ID.GE.245.AND.ID.LE.253).OR.
1537
1538 & (ID.EQ.254.AND.IDKPRD(4,IM).EQ.4)) THEN
1539
1540C bbar hadron or bbar decay of B_c+
1541
1542 IDQ=11
1543
1544 IS=NHEP+1
1545
1546 IQ=NHEP+2
1547
1548 ELSE
1549
1550C Decay not recognized
1551
1552 CALL HWWARN('HWDHVY',101,*999)
1553
1554 ENDIF
1555
1556C Label constituents
1557
1558 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWDHVY',102,*999)
1559
1560 ISTHEP(IHEP)=199
1561
1562 JDAHEP(1,IHEP)=NHEP+1
1563
1564 JDAHEP(2,IHEP)=NHEP+2
1565
1566 IDHW(IQ)=IDQ
1567
1568 IDHW(IS)=IDKPRD(4,IM)
1569
1570 IDHEP(IQ)=IDPDG(IDQ)
1571
1572 IDHEP(IS)=IDPDG(IDKPRD(4,IM))
1573
1574 ISTHEP(IQ)=155
1575
1576 ISTHEP(IS)=115
1577
1578 JMOHEP(1,IQ)=IHEP
1579
1580 JMOHEP(2,IQ)=IS
1581
1582 JDAHEP(1,IQ)=NHEP+3
1583
1584 JDAHEP(2,IQ)=NHEP+5
1585
1586 JMOHEP(1,IS)=IHEP
1587
1588 JMOHEP(2,IS)=NHEP+5
1589
1590 JDAHEP(1,IS)=0
1591
1592 JDAHEP(2,IS)=NHEP+5
1593
1594 NHEP=NHEP+2
1595
1596C and weak decay product jets
1597
1598 DO 10 J=1,3
1599
1600 NHEP=NHEP+1
1601
1602 IDHW(NHEP)=IDKPRD(J,IM)
1603
1604 IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1605
1606 ISTHEP(NHEP)=IST(J)
1607
1608 JMOHEP(1,NHEP)=IQ
1609
1610 JDAHEP(1,NHEP)=0
1611
1612 10 PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1613
1614 JMOHEP(2,NHEP-2)=NHEP-1
1615
1616 JDAHEP(2,NHEP-2)=NHEP-1
1617
1618 JMOHEP(2,NHEP-1)=NHEP-2
1619
1620 JDAHEP(2,NHEP-1)=NHEP-2
1621
1622 JMOHEP(2,NHEP )=IQ
1623
1624 JDAHEP(2,NHEP )=IQ
1625
1626C Share momenta in ratio of masses, preserving specator mass
1627
1628 XS=RMASS(IDHW(IS))/PHEP(5,IHEP)
1629
1630 XB=ONE-XS
1631
1632 CALL HWVSCA(5,XB,PHEP(1,IHEP),PHEP(1,IQ))
1633
1634 CALL HWVSCA(5,XS,PHEP(1,IHEP),PHEP(1,IS))
1635
1636 IF (NME(IM).EQ.100) THEN
1637
1638C Generate decay momenta using full (V-A)*(V-A) matrix element
1639
1640 EMWSQ=RMASS(198)**2
1641
1642 GMWSQ=(RMASS(198)*GAMW)**2
1643
1644 EMLIM=GMWSQ+(EMWSQ-(PHEP(5,IQ)-PHEP(5,NHEP))**2)**2
1645
1646 20 CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-1),
1647
1648 & PHEP(1,NHEP-2),PHEP(1,NHEP),HWDWWT)
1649
1650 CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1651
1652 EMTST=(HWULDO(PW,PW)-EMWSQ)**2
1653
1654 IF ((EMTST+GMWSQ)*HWR().GT.EMLIM) GOTO 20
1655
1656 ELSE
1657
1658C Use phase space
1659
1660 CALL HWDTHR(PHEP(1,IQ ),PHEP(1,NHEP-2),
1661
1662 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1663
1664 CALL HWVSUM(4,PHEP(1,NHEP-2),PHEP(1,NHEP-1),PW)
1665
1666 ENDIF
1667
1668C Set up production vertices
1669
1670 CALL HWVZRO(4,VHEP(1,IQ))
1671
1672 CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,IS))
1673
1674 CALL HWVEQU(4,VHEP(1,IQ),VHEP(1,NHEP))
1675
1676 CALL HWUDKL(198,PW,VHEP(1,NHEP-2))
1677
1678 CALL HWVSUM(4,VHEP(1,IQ),VHEP(1,NHEP-2),VHEP(1,NHEP-2))
1679
1680 CALL HWVEQU(4,VHEP(1,NHEP-2),VHEP(1,NHEP-1))
1681
1682 EMSCA=PHEP(5,IQ)
1683
1684 ELSE
1685
1686C Quarkonium decay
1687
1688C Label products
1689
1690 ISTHEP(IHEP)=199
1691
1692 JDAHEP(1,IHEP)=NHEP+1
1693
1694 DO 30 J=1,NPRODS(IM)
1695
1696 NHEP=NHEP+1
1697
1698 IDHW(NHEP)=IDKPRD(J,IM)
1699
1700 IDHEP(NHEP)=IDPDG(IDKPRD(J,IM))
1701
1702 ISTHEP(NHEP)=IST(J)
1703
1704 JMOHEP(1,NHEP)=IHEP
1705
1706 JDAHEP(1,NHEP)=0
1707
1708 PHEP(5,NHEP)=RMASS(IDKPRD(J,IM))
1709
1710 30 CALL HWVZRO(4,VHEP(1,NHEP))
1711
1712 JDAHEP(2,IHEP)=NHEP
1713
1714C Establish colour connections and select momentum configuration
1715
1716 IF (NPRODS(IM).EQ.3) THEN
1717
1718 IF (IDKPRD(3,IM).EQ.13) THEN
1719
1720C 3-gluon decay
1721
1722 JMOHEP(2,NHEP-2)=NHEP
1723
1724 JMOHEP(2,NHEP-1)=NHEP-2
1725
1726 JMOHEP(2,NHEP )=NHEP-1
1727
1728 JDAHEP(2,NHEP-2)=NHEP-1
1729
1730 JDAHEP(2,NHEP-1)=NHEP
1731
1732 JDAHEP(2,NHEP )=NHEP-2
1733
1734 ELSE
1735
1736C or 2-gluon + photon decay
1737
1738 JMOHEP(2,NHEP-2)=NHEP-1
1739
1740 JMOHEP(2,NHEP-1)=NHEP-2
1741
1742 JMOHEP(2,NHEP )=NHEP
1743
1744 JDAHEP(2,NHEP-2)=NHEP-1
1745
1746 JDAHEP(2,NHEP-1)=NHEP-2
1747
1748 JDAHEP(2,NHEP )=NHEP
1749
1750 ENDIF
1751
1752 IF (NME(IM).EQ.130) THEN
1753
1754C Use Ore & Powell orthopositronium matrix element
1755
1756 40 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1757
1758 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1759
1760 X1=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-2))/PHEP(5,IHEP)**2
1761
1762 X2=TWO*HWULDO(PHEP(1,IHEP),PHEP(1,NHEP-1))/PHEP(5,IHEP)**2
1763
1764 X3=TWO-X1-X2
1765
1766 TEST=((X1*(ONE-X1))**2+(X2*(ONE-X2))**2+(X3*(ONE-X3))**2)
1767
1768 & /(X1*X2*X3)**2
1769
1770 IF (TEST.LT.TWO*HWR()) GOTO 40
1771
1772 ELSE
1773
1774C Use phase space
1775
1776 CALL HWDTHR(PHEP(1,IHEP),PHEP(1,NHEP-2),
1777
1778 & PHEP(1,NHEP-1),PHEP(1,NHEP),HWDPWT)
1779
1780 ENDIF
1781
1782 ELSE
1783
1784C Parapositronium 2-gluon or q-qbar decay
1785
1786 JMOHEP(2,NHEP-1)=NHEP
1787
1788 JMOHEP(2,NHEP )=NHEP-1
1789
1790 JDAHEP(2,NHEP-1)=NHEP
1791
1792 JDAHEP(2,NHEP )=NHEP-1
1793
1794 CALL HWDTWO(PHEP(1,IHEP),PHEP(1,NHEP-1),
1795
1796 & PHEP(1,NHEP),CMMOM(IM),TWO,.FALSE.)
1797
1798 ENDIF
1799
1800 EMSCA=PHEP(5,IHEP)
1801
1802 ENDIF
1803
1804C Process this new hard scatter
1805
1806 CALL HWVEQU(4,VTXQDK(1,I),VTXPIP)
1807
1808 CALL HWBGEN
1809
1810 CALL HWCFOR
1811
1812 CALL HWCDEC
1813
1814 CALL HWDHAD
1815
1816 100 CONTINUE
1817
1818 NQDK=0
1819
1820 999 END
1821
1822CDECK ID>, HWDRCL.
1823
1824*CMZ :- -20/07/99 10:56:12 by Peter Richardson
1825
1826*-- Author : Peter Richardson
1827
1828C-----------------------------------------------------------------------
1829
1830 SUBROUTINE HWDRCL(IHEP,MHEP,CLSAVE)
1831
1832C-----------------------------------------------------------------------
1833
1834C Sets the colour connections in Baryon number violating decays
1835
1836C-----------------------------------------------------------------------
1837
1838 INCLUDE 'HERWIG61.INC'
1839
1840 INTEGER IHEP,MHEP,ID,ID2,IDM2,IDM3,COLCON(2,2,3),FLACON(2,3),JHEP,
1841
1842 & DECAY,COLANT,KHEP,IDM,IDMB,IDMB2,IDMB3,IDMB4,QHEP,IDM4,
1843
1844 & CLSAVE(2),XHEP,I,HWRINT,THEP
1845
1846 LOGICAL CONBV
1847
1848C--Colour connections for the decays
1849
1850 DATA COLCON/-1,1,-1,-2,-2,1,-3,-1,-1,1,-2,-1/
1851
1852 DATA FLACON/1,-1,1,-1,-1,0/
1853
1854C--identify the decay
1855
1856 IF(IERROR.NE.0) RETURN
1857
1858 ID = IDHW(IHEP)
1859
1860 ID2 = IDHW(MHEP)
1861
1862 IF(ID.GE.450.AND.ID.LE.457) THEN
1863
1864 DECAY = 1
1865
1866 ELSEIF(ID.EQ.449) THEN
1867
1868 DECAY = 2
1869
1870 ELSEIF((ID.GE.411.AND.ID.LE.424).OR.ID.EQ.405.OR.ID.EQ.406) THEN
1871
1872 DECAY = 3
1873
1874 ELSE
1875
1876C--UNKNOWN DECAY
1877
1878 CALL HWWARN('HWDRCL',100,*999)
1879
1880 ENDIF
1881
1882 COLANT = 1
1883
1884C--identify the colour partner
1885
1886 IF(DECAY.GT.1.AND.ID2.LE.6) THEN
1887
1888C--colour partner
1889
1890 COLANT = 2
1891
1892 KHEP = JDAHEP(2,IHEP-1)
1893
1894 ELSEIF(DECAY.GT.1.AND.ID2.GE.7) THEN
1895
1896C--anticolour partner
1897
1898 COLANT = 3
1899
1900 KHEP = JMOHEP(2,IHEP)
1901
1902 ELSE
1903
1904 KHEP=IHEP
1905
1906 ENDIF
1907
1908 IDM = IDHW(JMOHEP(1,KHEP))
1909
1910 IF(ABS(IDPDG(IDM)).GT.1000000.OR.IDM.EQ.15) THEN
1911
1912 IDM2 = IDHW(JDAHEP(1,JMOHEP(1,KHEP)))
1913
1914 IDM3 = IDHW(JDAHEP(2,JMOHEP(1,KHEP)))
1915
1916 IDM4 = IDHW(JDAHEP(2,JMOHEP(1,KHEP))-1)
1917
1918 QHEP = JMOHEP(1,KHEP)
1919
1920 IDMB = IDHW(JMOHEP(1,QHEP))
1921
1922 IDMB2 = IDHW(JMOHEP(2,QHEP))
1923
1924 IDMB3 = IDHW(JDAHEP(1,QHEP))
1925
1926 IDMB4 = IDHW(JDAHEP(2,QHEP))
1927
1928 ENDIF
1929
1930C--Now decide if the colour partner decayed via BV
1931
1932 IF(COLANT.EQ.2.AND.((((IDM.GE.413.AND.IDM.LE.418).OR.
1933
1934 & IDM.EQ.449.OR.IDM.EQ.405.OR.IDM.EQ.406).AND.
1935
1936 & (IDM2.GE.7.AND.IDM2.LE.12.AND.
1937
1938 & IDM3.GE.7.AND.IDM3.LE.12.AND.
1939
1940 & IDM4.GE.7.AND.IDM4.LE.12)).OR.
1941
1942 & (IDM.EQ.15.AND.IDMB.LE.6.AND.IDMB2.LE.6.AND.
1943
1944 & ((IDMB3.GE.7.AND.IDMB4.GE.12.AND.IDMB4.EQ.449).OR.
1945
1946 & (IDMB3.GE.198.AND.IDMB3.LE.207.AND.
1947
1948 & ABS(IDPDG(IDMB4)).GT.1000000))))) THEN
1949
1950 CONBV = .TRUE.
1951
1952 COLUPD = .TRUE.
1953
1954 HVFCEN = .FALSE.
1955
1956 XHEP = JMOHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1957
1958 ELSEIF(COLANT.EQ.3.AND.((((IDM.GE.419.AND.IDM.LE.424).OR.
1959
1960 & IDM.EQ.449.OR.IDM.EQ.411.OR.IDM.EQ.412).AND.
1961
1962 & (IDM2.LE.6.AND.IDM3.LE.6.AND.IDM4.LE.6)).OR.
1963
1964 & (IDM.EQ.15.AND.IDMB.GE.7.AND.IDMB.LE.12.AND.
1965
1966 & IDMB2.GE.7.AND.IDMB2.LE.12.AND.((IDMB3.LE.6.AND.
1967
1968 & IDMB4.EQ.449).OR.(ABS(IDPDG(IDMB4)).GT.1000000
1969
1970 & .AND.IDMB3.GE.198.AND.IDMB3.LE.207))))) THEN
1971
1972 CONBV = .TRUE.
1973
1974 COLUPD = .TRUE.
1975
1976 HVFCEN = .FALSE.
1977
1978 XHEP = JDAHEP(2,JDAHEP(2,JMOHEP(1,KHEP)))
1979
1980 ELSE
1981
1982 CONBV = .FALSE.
1983
1984 COLUPD = .FALSE.
1985
1986 XHEP = 0
1987
1988 ENDIF
1989
1990 IF(CONBV) THEN
1991
1992 IF(IDM.NE.15) THEN
1993
1994 CLSAVE(1) = JDAHEP(2,JMOHEP(1,KHEP))-1
1995
1996 CLSAVE(2) = CLSAVE(1)+1
1997
1998 ELSE
1999
2000 IF(IDMB4.EQ.449) THEN
2001
2002 DO I=1,2
2003
2004 CLSAVE(I) = JMOHEP(I,JMOHEP(1,KHEP))
2005
2006 IF(CLSAVE(I).EQ.XHEP) CLSAVE(I)=JDAHEP(1,JMOHEP(1,KHEP))
2007
2008 ENDDO
2009
2010 ELSE
2011
2012 CLSAVE(1) = JMOHEP(1,JMOHEP(1,KHEP))
2013
2014 CLSAVE(2) = JMOHEP(2,JMOHEP(1,KHEP))
2015
2016 ENDIF
2017
2018 ENDIF
2019
2020 ELSE
2021
2022 CLSAVE(1)=0
2023
2024 CLSAVE(2)=0
2025
2026 ENDIF
2027
2028C--Now set the colours for angular ordering
2029
2030 THEP = MHEP-1
2031
2032 IF(DECAY.EQ.1) THEN
2033
2034 IF(ID2.LE.6) THEN
2035
2036 JMOHEP(2,THEP) = THEP+HWRINT(1,2)
2037
2038 JDAHEP(2,THEP) = THEP
2039
2040 ELSE
2041
2042 JMOHEP(2,THEP) = THEP
2043
2044 JDAHEP(2,THEP) = THEP+HWRINT(1,2)
2045
2046 ENDIF
2047
2048 ELSEIF(DECAY.EQ.2) THEN
2049
2050 IF(ID2.LE.6) THEN
2051
2052 JMOHEP(2,THEP) = IHEP
2053
2054 JDAHEP(2,THEP) = THEP
2055
2056 ELSE
2057
2058 JMOHEP(2,THEP) = THEP
2059
2060 JDAHEP(2,THEP) = IHEP
2061
2062 ENDIF
2063
2064 ENDIF
2065
2066C--Colour of the second two
2067
2068 DO JHEP=1,2
2069
2070 IF(ID2.LE.6) THEN
2071
2072 JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2073
2074 & COLCON(HWRINT(1,2),JHEP,DECAY)
2075
2076 JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2077
2078 ELSE
2079
2080 JDAHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+
2081
2082 & COLCON(HWRINT(1,2),JHEP,DECAY)
2083
2084 JMOHEP(2,MHEP+JHEP-1) = MHEP+JHEP-1+FLACON(JHEP,DECAY)
2085
2086 ENDIF
2087
2088 ENDDO
2089
2090C--Now set the colours of the colour partner
2091
2092 IF(DECAY.GT.1.AND..NOT.CONBV) THEN
2093
2094 IF(ID2.LE.6) JMOHEP(2,KHEP) = MHEP+HWRINT(0,1)
2095
2096 IF(ID2.GE.7) JDAHEP(2,KHEP) = MHEP+HWRINT(0,1)
2097
2098 ELSEIF(CONBV) THEN
2099
2100 IF(ID2.GT.6) THEN
2101
2102 JMOHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2103
2104 IF(JMOHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2105
2106 JMOHEP(2,CLSAVE(2)) = MHEP+1
2107
2108 ELSE
2109
2110 JMOHEP(2,CLSAVE(2)) = MHEP
2111
2112 ENDIF
2113
2114 ELSE
2115
2116 JDAHEP(2,CLSAVE(1)) = MHEP+HWRINT(0,1)
2117
2118 IF(JDAHEP(2,CLSAVE(1)).EQ.MHEP) THEN
2119
2120 JDAHEP(2,CLSAVE(2)) = MHEP+1
2121
2122 ELSE
2123
2124 JDAHEP(2,CLSAVE(2)) = MHEP
2125
2126 ENDIF
2127
2128 ENDIF
2129
2130 ENDIF
2131
2132 999 END
2133
2134CDECK ID>, HWDRME.
2135
2136*CMZ :- -20/07/99 10:56:12 by Peter Richardson
2137
2138*-- Author : Peter Richardson
2139
2140C-----------------------------------------------------------------------
2141
2142 SUBROUTINE HWDRME(LHEP,MHEP)
2143
2144C-----------------------------------------------------------------------
2145
2146C SUBROUTINE TO IMPLEMENT ALL RPARITY DECAY MATRIX ELEMENTS
2147
2148C-----------------------------------------------------------------------
2149
2150 INCLUDE 'HERWIG61.INC'
2151
2152 DOUBLE PRECISION SM(6),SW(6),HWULDO,INFCOL,AM, M12SQ,M23SQ,MSGN,
2153
2154 & M13SQ,A(6),B(6),SWEAK,MW,DECMOM(5),TEST(4),EPS,
2155
2156 & M12SQT(6),M23SQT(6),M13SQT(6),LIMIT,M(4),RAND,
2157
2158 & MC(2),MX2(6),MX(6),HWDPWT,HWR,HWDRM1,LAMD(3)
2159
2160 EXTERNAL HWDRM1,HWULDO,HWDPWT,HWR
2161
2162 INTEGER K,SN(3),LHEP,CSP,I,SB(3),J,ND,LTRY,MHEP,NSP,ID(3),IG,
2163
2164 & IDHWTP,IDHPTP,MTRY
2165
2166 PARAMETER(EPS=1D-20)
2167
2168 IF(IERROR.NE.0) RETURN
2169
2170C--Electroweak parameters, etc
2171
2172 SWEAK = SQRT(SWEIN)
2173
2174 MW = RMASS(198)
2175
2176 M(4) = PHEP(5,LHEP)
2177
2178 IG = IDHW(LHEP)
2179
2180C--Find the masses of the final state and zero parameters
2181
2182 DO K=1,3
2183
2184 ID(K) = IDHW(MHEP+K-1)
2185
2186 IF(ID(K).LE.12) THEN
2187
2188 SN(K)=ID(K)
2189
2190 ELSE
2191
2192 SN(K)=ID(K)-120
2193
2194 ENDIF
2195
2196 IF(SN(K).GT.6) SN(K)=SN(K)-6
2197
2198 M(K) = PHEP(5,LHEP+K)
2199
2200 SB(K)=SN(K)
2201
2202 LAMD(K) = ZERO
2203
2204 ENDDO
2205
2206 DO J=1,6
2207
2208 MX2(J) = ZERO
2209
2210 MX(J) = ZERO
2211
2212 M13SQT(J) = ZERO
2213
2214 M23SQT(J) = ZERO
2215
2216 M12SQT(J) = ZERO
2217
2218 ENDDO
2219
2220C--Evaluate the coefficents for the mode we want
2221
2222 IF(IG.GE.450.AND.IG.LE.453) THEN
2223
2224C--NEUTRALINO
2225
2226 NSP = IG-449
2227
2228 AM = RMASS(IG)
2229
2230 MSGN = ZSGNSS(NSP)
2231
2232 MC(1) = ZMIXSS(NSP,3)/(2*MW*COSB*SWEAK)
2233
2234 MC(2) = ZMIXSS(NSP,4)/(2*MW*SINB*SWEAK)
2235
2236C--Calculate the combinations of couplings needed
2237
2238 IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2239
2240C--first for the UDD modes
2241
2242 DO J=1,2
2243
2244 A(J) = M(1)*MC(2)*QMIXSS(SN(1),2,J)
2245
2246 & +SLFCH(SN(1),NSP)*QMIXSS(SN(1),1,J)
2247
2248 B(J) = MSGN*(M(1)*MC(2)*QMIXSS(SN(1),1,J)
2249
2250 & +SRFCH(SN(1),NSP)*QMIXSS(SN(1),2,J))
2251
2252 MX2(J) = QMIXSS(SN(1),2,J)
2253
2254 A(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2255
2256 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2257
2258 B(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2259
2260 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2261
2262 MX2(J+2) = QMIXSS(SN(2),2,J)
2263
2264 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2265
2266 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2267
2268 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2269
2270 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2271
2272 MX2(J+2) = QMIXSS(SN(3),2,J)
2273
2274 ENDDO
2275
2276 DO K=1,3
2277
2278 SN(K) = SN(K)+400
2279
2280 SB(K) = SB(K)+412
2281
2282 ENDDO
2283
2284 ELSEIF(ID(1).GE.121.AND.ID(2).GE.121.AND.ID(3).GE.121) THEN
2285
2286C--Now for the LLE modes
2287
2288 DO J=1,2
2289
2290 A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2291
2292 & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2293
2294 B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2295
2296 & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(2),1,J)
2297
2298 MX2(J)= LMIXSS(SN(1),1,J)
2299
2300 A(J+2) = ZERO
2301
2302 B(J+2) = SLFCH(10+SN(2),NSP)*LMIXSS(SN(2),1,J)
2303
2304 MX2(J+2) = LMIXSS(SN(2),1,J)
2305
2306 A(J+4) = M(3)*MC(1)*LMIXSS(SN(3),2,J)
2307
2308 & +SLFCH(10+SN(3),NSP)*LMIXSS(SN(3),1,J)
2309
2310 B(J+4) = MSGN*(M(3)*MC(1)*LMIXSS(SN(3),1,J)
2311
2312 & +SRFCH(10+SN(3),NSP)*LMIXSS(SN(3),2,J))
2313
2314 MX2(4+J) = LMIXSS(SN(3),2,J)
2315
2316 ENDDO
2317
2318 DO J=1,3
2319
2320 SN(J) = SN(J) + 424
2321
2322 SB(J) = SB(J) + 436
2323
2324 ENDDO
2325
2326 ELSE
2327
2328C--Now for both types of LQD modes
2329
2330 IF(MOD(SN(1),2).EQ.0) THEN
2331
2332C--First the neutrino,down,antidown mode
2333
2334 DO J=1,2
2335
2336 A(J) = ZERO
2337
2338 B(J) = SLFCH(10+SN(1),NSP)*
2339
2340 & LMIXSS(SN(1),1,J)
2341
2342 MX2(J) = LMIXSS(SN(1),1,J)
2343
2344 A(J+2) = MSGN*(M(2)*MC(1)*QMIXSS(SN(2),1,J)
2345
2346 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2347
2348 B(J+2) = M(2)*MC(1)*QMIXSS(SN(2),2,J)
2349
2350 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2351
2352 MX2(2+J) = QMIXSS(SN(2),1,J)
2353
2354 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2355
2356 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2357
2358 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2359
2360 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2361
2362 MX2(J+4) = QMIXSS(SN(3),2,J)
2363
2364 ENDDO
2365
2366 ELSE
2367
2368C--Now the charged lepton, antiup,down modes
2369
2370 DO J=1,2
2371
2372 A(J) = MSGN*(M(1)*MC(1)*LMIXSS(SN(1),1,J)
2373
2374 & +SRFCH(10+SN(1),NSP)*LMIXSS(SN(1),2,J))
2375
2376 B(J) = M(1)*MC(1)*LMIXSS(SN(1),2,J)
2377
2378 & +SLFCH(10+SN(1),NSP)*LMIXSS(SN(1),1,J)
2379
2380 MX2(J) = LMIXSS(SN(1),1,J)
2381
2382 A(J+2) =MSGN*(M(2)*MC(2)*QMIXSS(SN(2),1,J)
2383
2384 & +SRFCH(SN(2),NSP)*QMIXSS(SN(2),2,J))
2385
2386 B(J+2) = M(2)*MC(2)*QMIXSS(SN(2),2,J)
2387
2388 & +SLFCH(SN(2),NSP)*QMIXSS(SN(2),1,J)
2389
2390 MX2(2+J) = QMIXSS(SN(2),1,J)
2391
2392 A(J+4) = M(3)*MC(1)*QMIXSS(SN(3),2,J)
2393
2394 & +SLFCH(SN(3),NSP)*QMIXSS(SN(3),1,J)
2395
2396 B(J+4) = MSGN*(M(3)*MC(1)*QMIXSS(SN(3),1,J)
2397
2398 & +SRFCH(SN(3),NSP)*QMIXSS(SN(3),2,J))
2399
2400 MX2(J+4) = QMIXSS(SN(3),2,J)
2401
2402 ENDDO
2403
2404 ENDIF
2405
2406 SN(1) = SN(1) + 424
2407
2408 SB(1) = SB(1) + 436
2409
2410 DO J=2,3
2411
2412 SN(J) = SN(J) + 400
2413
2414 SB(J) = SB(J) + 412
2415
2416 ENDDO
2417
2418 ENDIF
2419
2420 DO K=1,3
2421
2422 SM(2*K-1) = RMASS(SN(K))
2423
2424 SM(2*K) = RMASS(SB(K))
2425
2426 SW(2*K-1) = HBAR/RLTIM(SN(K))
2427
2428 SW(2*K) = HBAR/RLTIM(SB(K))
2429
2430 ENDDO
2431
2432 ND = 3
2433
2434 DO K=1,3
2435
2436 LAMD(K) = ONE
2437
2438 ENDDO
2439
2440 INFCOL = ONE
2441
2442 ELSEIF(IG.EQ.449) THEN
2443
2444C--GLUINO
2445
2446C--First obtian the masses and widths needed
2447
2448 AM = RMASS(IG)
2449
2450 ND = 3
2451
2452C--Calculate the combinations of couplings needed
2453
2454 IF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2455
2456C--first for the UDD modes
2457
2458 INFCOL = -0.5D0
2459
2460C--Couplings
2461
2462 DO I=1,3
2463
2464 DO J=1,2
2465
2466 A(2*I-2+J) = -QMIXSS(SN(I),1,J)
2467
2468 B(2*I-2+J) = QMIXSS(SN(I),2,J)
2469
2470 MX2(2*I-2+J) = QMIXSS(SN(I),2,J)
2471
2472 ENDDO
2473
2474 SN(I) = SN(I)+400
2475
2476 SB(I) = SB(I)+412
2477
2478 ENDDO
2479
2480 ELSE
2481
2482 INFCOL = ONE
2483
2484C--Now for both types of LQD modes
2485
2486 IF(MOD(SN(1),2).EQ.0) THEN
2487
2488C--First the neutrino,down,antidown mode
2489
2490 DO J=1,2
2491
2492 A(J) = ZERO
2493
2494 B(J) = ZERO
2495
2496 MX2(J) = ZERO
2497
2498 A(J+2) = QMIXSS(SN(2),2,J)
2499
2500 B(J+2) = -QMIXSS(SN(2),1,J)
2501
2502 MX2(J+2) = QMIXSS(SN(2),1,J)
2503
2504 A(J+4) = -QMIXSS(SN(3),1,J)
2505
2506 B(J+4) = QMIXSS(SN(3),2,J)
2507
2508 MX2(4+J) = QMIXSS(SN(3),2,J)
2509
2510 ENDDO
2511
2512 ELSEIF(MOD(SN(1),2).EQ.1) THEN
2513
2514C--Now the charged lepton, antiup,down modes
2515
2516 DO J=1,2
2517
2518 A(J) = ZERO
2519
2520 B(J) = ZERO
2521
2522 MX2(J) = ZERO
2523
2524 A(J+2) = QMIXSS(SN(2),2,J)
2525
2526 B(J+2) = -QMIXSS(SN(2),1,J)
2527
2528 MX2(J+2) = QMIXSS(SN(2),1,J)
2529
2530 A(J+4) = -QMIXSS(SN(3),1,J)
2531
2532 B(J+4) = QMIXSS(SN(3),2,J)
2533
2534 MX2(J+4) = QMIXSS(SN(3),2,J)
2535
2536 ENDDO
2537
2538 ENDIF
2539
2540 SN(1) = SN(1) + 424
2541
2542 SB(1) = SB(1) + 436
2543
2544 DO K=2,3
2545
2546 SN(K) = SN(K) + 400
2547
2548 SB(K) = SB(K) + 412
2549
2550 ENDDO
2551
2552 ENDIF
2553
2554 DO K=1,3
2555
2556 SM(2*K-1) = RMASS(SN(K))
2557
2558 SM(2*K) = RMASS(SB(K))
2559
2560 SW(2*K-1) = HBAR/RLTIM(SN(K))
2561
2562 SW(2*K) = HBAR/RLTIM(SB(K))
2563
2564 ENDDO
2565
2566 DO K=1,3
2567
2568 LAMD(K) = ONE
2569
2570 ENDDO
2571
2572 ELSEIF(IG.GE.454.AND.IG.LE.457) THEN
2573
2574C--CHARGINO
2575
2576 CSP = IG-453
2577
2578 IF(CSP.GT.2) CSP = CSP-2
2579
2580 AM = RMASS(IG)
2581
2582 INFCOL = -ONE
2583
2584 MSGN = WSGNSS(CSP)
2585
2586 MC(1) = ONE/(SQRT(2.0D0)*MW*COSB)
2587
2588 MC(2) = ONE/(SQRT(2.0D0)*MW*SINB)
2589
2590C--Calculate the combinations of the couplings needed
2591
2592 IF(ID(1).GT.120.AND.ID(2).GT.120.AND.ID(3).GT.120) THEN
2593
2594C--first for the LLE modes, three modes
2595
2596 IF(MOD(SN(1),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2597
2598C--the one diagram mode nubar,positron, nu
2599
2600 DO J=1,2
2601
2602 A(J+4) = LMIXSS(SN(3)-1,1,J)*WMXUSS(CSP,1)
2603
2604 & -RMASS(SN(3)+119)*MC(1)*LMIXSS(SN(3)-1,2,J)*WMXUSS(CSP,2)
2605
2606 B(J+4) = ZERO
2607
2608 MX2(J+4) = LMIXSS(SN(3)-1,2,J)
2609
2610 ENDDO
2611
2612 ND = 1
2613
2614 SN(3) = SN(3)+423
2615
2616 SB(3) = SB(3)+435
2617
2618 ELSEIF(MOD(SN(1),2).EQ.0.AND.MOD(SN(2),2).EQ.0) THEN
2619
2620C--the first two diagram mode nu, nu, positron
2621
2622 DO J=1,2
2623
2624 A(J) = ZERO
2625
2626 B(J) = LMIXSS(SN(1)-1,1,J)*WMXUSS(CSP,1)
2627
2628 & -RMASS(SN(1)+119)*MC(1)*LMIXSS(SN(1)-1,2,J)*WMXUSS(CSP,2)
2629
2630 A(J+2) = ZERO
2631
2632 B(J+2) = LMIXSS(SN(2)-1,1,J)*WMXUSS(CSP,1)
2633
2634 & -RMASS(SN(2)+119)*MC(1)*LMIXSS(SN(2)-1,2,J)*WMXUSS(CSP,2)
2635
2636 MX2(J) = LMIXSS(SN(1)-1,1,J)
2637
2638 MX2(J+2) = LMIXSS(SN(2)-1,1,J)
2639
2640 ENDDO
2641
2642 ND = 2
2643
2644 DO J=1,2
2645
2646 SN(J) = SN(J)+423
2647
2648 SB(J) = SB(J)+435
2649
2650 ENDDO
2651
2652 ELSE
2653
2654C--the second two diagram mode positron, positron, electron
2655
2656 DO J=1,2
2657
2658 A(J) = -M(1)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(1)+1,1,J)
2659
2660 B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(1)+1,1,J)
2661
2662 A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*LMIXSS(SN(2)+1,1,J)
2663
2664 B(J+2) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2665
2666 MX2(J) = LMIXSS(SN(1)+1,1,J)
2667
2668 MX2(J+2) = LMIXSS(SN(2)+1,1,J)
2669
2670 ENDDO
2671
2672 DO J=1,2
2673
2674 SN(J) = SN(J)+425
2675
2676 SB(J) = SB(J)+437
2677
2678 ENDDO
2679
2680 ND = 2
2681
2682 ENDIF
2683
2684 DO K=1,3
2685
2686 LAMD(K) = ONE
2687
2688 ENDDO
2689
2690 ELSEIF(ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
2691
2692C--now for the UDD
2693
2694 IF(MOD(SN(1),2).EQ.0) THEN
2695
2696C--two diagram mode
2697
2698 LAMD(1) = LAMDA3(SN(2)/2,SN(1)/2,(SN(3)+1)/2)
2699
2700 LAMD(2) = LAMDA3(SN(1)/2,SN(2)/2,(SN(3)+1)/2)
2701
2702 DO J=1,2
2703
2704 A(J) = WMXUSS(CSP,1)*QMIXSS(SN(1)-1,1,J)
2705
2706 & -RMASS(SN(1)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(1)-1,2,J)
2707
2708 B(J) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(1)-1,1,J)
2709
2710 A(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2711
2712 & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2713
2714 B(J+2) = -MSGN*M(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)-1,1,J)
2715
2716 MX2(J) = QMIXSS(SN(1)-1,2,J)
2717
2718 MX2(J+2) = QMIXSS(SN(2)-1,2,J)
2719
2720 ENDDO
2721
2722 DO J=1,2
2723
2724 SN(J) = SN(J) + 399
2725
2726 SB(J) = SB(J) + 411
2727
2728 ENDDO
2729
2730 ND = 2
2731
2732 ELSE
2733
2734C--three diagram mode
2735
2736 LAMD(1) = LAMDA3((SN(1)+1)/2,(SN(2)+1)/2,(SN(3)+1)/2)
2737
2738 LAMD(2) = LAMDA3((SN(2)+1)/2,(SN(1)+1)/2,(SN(3)+1)/2)
2739
2740 LAMD(3) = LAMDA3((SN(3)+1)/2,(SN(2)+1)/2,(SN(1)+1)/2)
2741
2742 DO I=1,3
2743
2744 DO J=1,2
2745
2746 A(J+2*I-2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(I)+1,1,J)
2747
2748 & -RMASS(SN(I)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(I)+1,2,J))
2749
2750 B(J+2*I-2) = -M(I)*MC(1)*WMXUSS(CSP,2)
2751
2752 & *QMIXSS(SN(I)+1,1,J)
2753
2754 MX2(J+2*I-2) = QMIXSS(SN(I)+1,2,J)
2755
2756 ENDDO
2757
2758 SN(I) = SN(I) + 401
2759
2760 SB(I) = SB(I) + 413
2761
2762 ENDDO
2763
2764 ND = 3
2765
2766 ENDIF
2767
2768 ELSE
2769
2770C--now for the LQD modes
2771
2772 IF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.0) THEN
2773
2774C--first one diagram mode nubar, dbar, up
2775
2776 DO J=1,2
2777
2778 A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2779
2780 & QMIXSS(SN(3)-1,1,J)
2781
2782 B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2783
2784 & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2785
2786 MX2(J+4) = QMIXSS(SN(3)-1,2,J)
2787
2788 ENDDO
2789
2790 SN(3) = SN(3) + 399
2791
2792 SB(3) = SB(3) + 411
2793
2794 ND = 1
2795
2796 ELSEIF(MOD(SN(2),2).EQ.0.AND.MOD(SN(3),2).EQ.0) THEN
2797
2798C--second one diagram mode positron, ubar, up
2799
2800 DO J=1,2
2801
2802 A(J+4) = -MSGN*M(3)*WMXVSS(CSP,2)*MC(2)*
2803
2804 & QMIXSS(SN(3)-1,1,J)
2805
2806 B(J+4) = WMXUSS(CSP,1)*QMIXSS(SN(3)-1,1,J)
2807
2808 & -RMASS(SN(3)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(3)-1,2,1)
2809
2810 MX2(J+4) = QMIXSS(SN(3)-1,2,J)
2811
2812 ENDDO
2813
2814 SN(3) = SN(3) + 399
2815
2816 SB(3) = SB(3) + 411
2817
2818 ND = 1
2819
2820 ELSEIF(MOD(SN(2),2).EQ.1.AND.MOD(SN(3),2).EQ.1) THEN
2821
2822C--first two diagram mode positron, dbar, down
2823
2824 DO J=1,2
2825
2826 A(J) = -M(1)*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)+1,1,J)
2827
2828 B(J) = MSGN*WMXVSS(CSP,1)*LMIXSS(SN(2)+1,1,J)
2829
2830 A(J+2) = -M(2)*WMXUSS(CSP,2)*MC(1)*QMIXSS(SN(2)+1,1,J)
2831
2832 B(J+2) = MSGN*(WMXVSS(CSP,1)*QMIXSS(SN(2)+1,1,J)
2833
2834 & -RMASS(SN(2)+1)*MC(2)*WMXVSS(CSP,2)*QMIXSS(SN(2)+1,2,J))
2835
2836 MX2(J) = LMIXSS(SN(1)+1,1,J)
2837
2838 MX2(J+2) = QMIXSS(SN(2)+1,1,J)
2839
2840 ENDDO
2841
2842 SN(1) = SN(1) + 425
2843
2844 SB(1) = SB(1) + 437
2845
2846 SN(2) = SN(2) + 401
2847
2848 SB(2) = SB(2) + 413
2849
2850 ND = 2
2851
2852 ELSE
2853
2854C--second two diagram mode nu, up, dbar
2855
2856 DO J=1,2
2857
2858 A(J) = ZERO
2859
2860 B(J) = WMXUSS(CSP,1)*LMIXSS(SN(1)-1,1,J)
2861
2862 & -RMASS(119+SN(1))*MC(1)*WMXUSS(CSP,2)*LMIXSS(SN(1)-1,2,J)
2863
2864 A(J+2) = -MSGN*M(2)*MC(2)*WMXVSS(CSP,2)*
2865
2866 & QMIXSS(SN(2)-1,1,J)
2867
2868 B(J+2) = WMXUSS(CSP,1)*QMIXSS(SN(2)-1,1,J)
2869
2870 & -RMASS(SN(2)-1)*MC(1)*WMXUSS(CSP,2)*QMIXSS(SN(2)-1,2,J)
2871
2872 MX2(J) = LMIXSS(SN(1)-1,1,J)
2873
2874 MX2(J+2) = QMIXSS(SN(2)-1,1,J)
2875
2876 ENDDO
2877
2878 SN(1) = SN(1) + 423
2879
2880 SB(1) = SB(1) + 435
2881
2882 SN(2) = SN(2) + 399
2883
2884 SB(2) = SB(2) + 411
2885
2886 ND = 2
2887
2888 ENDIF
2889
2890 DO K=1,3
2891
2892 LAMD(K) = ONE
2893
2894 ENDDO
2895
2896 ENDIF
2897
2898 IF(ND.EQ.1) THEN
2899
2900 DO K=1,2
2901
2902 SM(2*K-1) = 0.0D0
2903
2904 SM(2*K) = 0.0D0
2905
2906 SW(2*K-1) = 0.0D0
2907
2908 SW(2*K) = 0.0D0
2909
2910 ENDDO
2911
2912 SM(5) = RMASS(SN(3))
2913
2914 SM(6) = RMASS(SB(3))
2915
2916 SW(5) = HBAR/RLTIM(SN(3))
2917
2918 SW(6) = HBAR/RLTIM(SB(3))
2919
2920 ELSE
2921
2922 DO K=1,2
2923
2924 SM(2*K-1) = RMASS(SN(K))
2925
2926 SM(2*K) = RMASS(SB(K))
2927
2928 SW(2*K-1) = HBAR/RLTIM(SN(K))
2929
2930 SW(2*K) = HBAR/RLTIM(SB(K))
2931
2932 SM(4+K) = ZERO
2933
2934 SW(4+K) = ZERO
2935
2936 ENDDO
2937
2938 ENDIF
2939
2940 ELSE
2941
2942C--UNKNOWN
2943
2944 CALL HWWARN('HWDRME',500,*999)
2945
2946 ENDIF
2947
2948C--Set mixing to zero if diagram not available
2949
2950 IF((AM.LT.(ABS(SM(1))+M(1)).OR.ABS(SM(1)).LT.(M(2)+M(3)))
2951
2952 & .AND.ABS(MX2(1)).GT.ZERO.AND.ND.NE.1) MX(1) = MX2(1)*LAMD(1)
2953
2954 IF((AM.LT.(ABS(SM(2))+M(1)).OR.ABS(SM(2)).LT.(M(2)+M(3)))
2955
2956 & .AND.ABS(MX2(2)).GT.ZERO.AND.ND.NE.1) MX(2) = MX2(2)*LAMD(1)
2957
2958 IF((AM.LT.(ABS(SM(3))+M(2)).OR.ABS(SM(3)).LT.(M(1)+M(3)))
2959
2960 & .AND.ABS(MX2(3)).GT.ZERO.AND.ND.NE.1) MX(3) = MX2(3)*LAMD(2)
2961
2962 IF((AM.LT.(ABS(SM(4))+M(2)).OR.ABS(SM(4)).LT.(M(1)+M(3)))
2963
2964 & .AND.ABS(MX2(4)).GT.ZERO.AND.ND.NE.1) MX(4) = MX2(4)*LAMD(2)
2965
2966 IF((AM.LT.(ABS(SM(5))+M(3)).OR.ABS(SM(5)).LT.(M(1)+M(2)))
2967
2968 & .AND.ABS(MX2(5)).GT.ZERO.AND.ND.NE.2) MX(5) = MX2(5)*LAMD(3)
2969
2970 IF((AM.LT.(ABS(SM(6))+M(3)).OR.ABS(SM(6)).LT.(M(1)+M(2)))
2971
2972 & .AND.ABS(MX2(6)).GT.ZERO.AND.ND.NE.2) MX(6) = MX2(6)*LAMD(3)
2973
2974C--Calculate the limiting points
2975
2976 DO J=1,2
2977
2978 IF(ND.NE.1) THEN
2979
2980 IF(ABS(MX(J)).GT.EPS) CALL HWDRM5(M23SQT(J),M13SQT(J),
2981
2982 & M12SQT(J),A(J),B(J),M(2),M(3),M(1),M(4),SM(J),SW(J))
2983
2984 IF(ABS(MX(J+2)).GT.EPS) CALL HWDRM5(M13SQT(2+J),M23SQT(2+J),
2985
2986 & M12SQT(2+J),A(2+J),B(2+J),M(1),M(3),M(2),M(4),SM(2+J),SW(2+J))
2987
2988 ENDIF
2989
2990 IF(ND.NE.2) THEN
2991
2992 IF(ABS(MX(J+4)).GT.EPS) CALL HWDRM5(M12SQT(4+J),M23SQT(4+J),
2993
2994 & M13SQT(4+J),A(4+J),B(4+J),M(1),M(2),M(3),M(4),SM(4+J),SW(4+J))
2995
2996 ENDIF
2997
2998 ENDDO
2999
3000C--Now evaluate the limit using these points
3001
3002 LIMIT = ZERO
3003
3004 DO 100 I=1,6
3005
3006 IF(ABS(MX(I)).LT.EPS) GOTO 100
3007
3008 LIMIT = LIMIT+HWDRM1(TEST,M12SQT(I),M13SQT(I),M23SQT(I),A,B,MX,
3009
3010 & M,SM,SW,INFCOL,AM,0,ND)
3011
3012 100 CONTINUE
3013
3014 LIMIT = TWO*LIMIT
3015
3016C--Now evaluate at a random point
3017
3018 MTRY = 0
3019
3020 25 MTRY = MTRY+1
3021
3022 LTRY = 0
3023
3024 35 LTRY = LTRY+1
3025
3026 CALL HWDTHR(PHEP(1,LHEP),PHEP(1,MHEP),
3027
3028 & PHEP(1,MHEP+1),PHEP(1,MHEP+2),HWDPWT)
3029
3030C--Now calculate the m12sq etc for the actual point
3031
3032 M12SQ = M(1)**2+M(2)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+1))
3033
3034 M13SQ = M(1)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP),PHEP(1,MHEP+2))
3035
3036 M23SQ = M(2)**2+M(3)**2+2*HWULDO(PHEP(1,MHEP+1),PHEP(1,MHEP+2))
3037
3038C--Now calulate the matrix element
3039
3040 TEST(4) = HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,
3041
3042 & M,SM,SW,INFCOL,AM,1,ND)
3043
3044C--Now test the value againest the limit
3045
3046 RAND = HWR()*LIMIT
3047
3048 IF(TEST(4).GT.LIMIT) THEN
3049
3050 LIMIT = 1.1D0*TEST(4)
3051
3052 CALL HWWARN('HWDRME',51,*150)
3053
3054 ENDIF
3055
3056 150 IF(TEST(4).LT.RAND.AND.LTRY.LT.NETRY) THEN
3057
3058 GOTO 35
3059
3060 ELSEIF(LTRY.GE.NETRY) THEN
3061
3062 IF(MTRY.LE.NETRY) THEN
3063
3064 LIMIT = LIMIT*0.9D0
3065
3066 CALL HWWARN('HWDRME',52,*25)
3067
3068 ELSE
3069
3070 CALL HWWARN('HWDRME',100,*999)
3071
3072 ENDIF
3073
3074 ENDIF
3075
3076C--Reorder the particles in gluino decay to get angular ordering right
3077
3078 IF(IG.EQ.449.AND.ID(1).LE.12.AND.ID(2).LE.12.AND.ID(3).LE.12) THEN
3079
3080 DO LTRY=1,3
3081
3082 IF(TEST(LTRY).GT.RAND) THEN
3083
3084 IF(LTRY.EQ.2) THEN
3085
3086 IDHWTP = IDHW(MHEP)
3087
3088 IDHW(MHEP) = IDHW(MHEP+1)
3089
3090 IDHW(MHEP+1) = IDHWTP
3091
3092 IDHPTP = IDHEP(MHEP)
3093
3094 IDHEP(MHEP) = IDHEP(MHEP+1)
3095
3096 IDHEP(MHEP+1) = IDHPTP
3097
3098 CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3099
3100 CALL HWVEQU(5,PHEP(1,MHEP+1),PHEP(1,MHEP))
3101
3102 CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+1))
3103
3104 ELSEIF(LTRY.EQ.3) THEN
3105
3106 IDHWTP = IDHW(MHEP)
3107
3108 IDHW(MHEP) = IDHW(MHEP+2)
3109
3110 IDHW(MHEP+2) = IDHWTP
3111
3112 IDHPTP = IDHEP(MHEP)
3113
3114 IDHEP(MHEP) = IDHEP(MHEP+2)
3115
3116 IDHEP(MHEP+2) = IDHPTP
3117
3118 DO I=1,5
3119
3120 CALL HWVEQU(5,PHEP(1,MHEP),DECMOM)
3121
3122 CALL HWVEQU(5,PHEP(1,MHEP+2),PHEP(1,MHEP))
3123
3124 CALL HWVEQU(5,DECMOM,PHEP(1,MHEP+2))
3125
3126 ENDDO
3127
3128 ENDIF
3129
3130 GOTO 52
3131
3132 ENDIF
3133
3134 RAND=RAND-TEST(LTRY)
3135
3136 ENDDO
3137
3138 ENDIF
3139
3140 52 CONTINUE
3141
3142 999 END
3143
3144CDECK ID>, HWDRM1.
3145
3146*CMZ :- -20/07/99 10:56:12 by Peter Richardson
3147
3148*-- Author : Peter Richardson
3149
3150C-----------------------------------------------------------------------
3151
3152 FUNCTION HWDRM1(TEST,M12SQ,M13SQ,M23SQ,A,B,MX,M,SM,SW
3153
3154 & ,INFCOL,AM,LM,ND)
3155
3156C-----------------------------------------------------------------------
3157
3158C FUNCTION TO GIVE THE R-PARITY VIOLATING MATRIX ELEMENT AT A GIVEN
3159
3160C PHASE-SPACE POINT
3161
3162C-----------------------------------------------------------------------
3163
3164 IMPLICIT NONE
3165
3166 DOUBLE PRECISION M12SQ,M13SQ,M23SQ,MX(6),A(6),B(6),SM(6),SW(6),
3167
3168 & INFCOL,AM,TERM(21),TEST(4),PLN,NPLN,ZERO,
3169
3170 & M(4),HWDRM1,HWDRM2,HWDRM3,HWDRM4
3171
3172 PARAMETER (ZERO=0)
3173
3174 EXTERNAL HWDRM2,HWDRM3,HWDRM4
3175
3176 INTEGER LM,K,ND
3177
3178C--Zero the array
3179
3180 DO K=1,21
3181
3182 TERM(K) = 0.0D0
3183
3184 ENDDO
3185
3186 HWDRM1 = 0.0D0
3187
3188C--The amplitude
3189
3190 IF(ABS(MX(1)).GT.ZERO.AND.ND.NE.1) THEN
3191
3192 TERM(1) = MX(1)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(1),
3193
3194 & SW(1),A(1),B(1))
3195
3196 IF(ABS(MX(2)).GT.ZERO) TERM(7)= MX(1)*MX(2)*HWDRM3(M23SQ,M(2),
3197
3198 & M(3),M(1),M(4),SM(1),SM(2),SW(1),SW(2),A(1),A(2),B(1),B(2))
3199
3200 IF(ABS(MX(3)).GT.ZERO) TERM(10)=-MX(1)*MX(3)*HWDRM4(M13SQ,M23SQ,
3201
3202 & M(1),M(3),M(2),M(4),SM(3),SM(1),SW(3),SW(1),A(1),A(3),B(1),B(3))
3203
3204 IF(ABS(MX(4)).GT.ZERO) TERM(11)=-MX(1)*MX(4)*HWDRM4(M13SQ,M23SQ,
3205
3206 & M(1),M(3),M(2),M(4),SM(4),SM(1),SW(4),SW(1),A(1),A(4),B(1),B(4))
3207
3208 IF(ABS(MX(5)).GT.ZERO) TERM(12)=-MX(1)*MX(5)*HWDRM4(M23SQ,M12SQ,
3209
3210 & M(3),M(2),M(1),M(4),SM(1),SM(5),SW(1),SW(5),A(5),A(1),B(5),B(1))
3211
3212 IF(ABS(MX(6)).GT.ZERO) TERM(13)=-MX(1)*MX(6)*HWDRM4(M23SQ,M12SQ,
3213
3214 & M(3),M(2),M(1),M(4),SM(1),SM(6),SW(1),SW(6),A(6),A(1),B(6),B(1))
3215
3216 ENDIF
3217
3218 IF(ABS(MX(2)).GT.ZERO.AND.ND.NE.1) THEN
3219
3220 TERM(2) = MX(2)**2*HWDRM2(M23SQ,M(2),M(3),M(1),M(4),SM(2),
3221
3222 & SW(2),A(2),B(2))
3223
3224 IF(ABS(MX(3)).GT.ZERO) TERM(14)=-MX(2)*MX(3)*HWDRM4(M13SQ,M23SQ,
3225
3226 & M(1),M(3),M(2),M(4),SM(3),SM(2),SW(3),SW(2),A(2),A(3),B(2),B(3))
3227
3228 IF(ABS(MX(4)).GT.ZERO) TERM(15)=-MX(2)*MX(4)*HWDRM4(M13SQ,M23SQ,
3229
3230 & M(1),M(3),M(2),M(4),SM(4),SM(2),SW(4),SW(2),A(2),A(4),B(2),B(4))
3231
3232 IF(ABS(MX(5)).GT.ZERO) TERM(16)=-MX(2)*MX(5)*HWDRM4(M23SQ,M12SQ,
3233
3234 & M(3),M(2),M(1),M(4),SM(2),SM(5),SW(2),SW(5),A(5),A(2),B(5),B(2))
3235
3236 IF(ABS(MX(6)).GT.ZERO) TERM(17)=-MX(2)*MX(6)*HWDRM4(M23SQ,M12SQ,
3237
3238 & M(3),M(2),M(1),M(4),SM(2),SM(6),SW(2),SW(6),A(6),A(2),B(6),B(2))
3239
3240 ENDIF
3241
3242 IF(ABS(MX(3)).GT.ZERO.AND.ND.NE.1) THEN
3243
3244 TERM(3) = MX(3)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(3),
3245
3246 & SW(3),A(3),B(3))
3247
3248 IF(ABS(MX(4)).GT.ZERO) TERM(8)= MX(3)*MX(4)*HWDRM3(M13SQ,M(1),
3249
3250 & M(3),M(2),M(4),SM(3),SM(4),SW(3),SW(4),A(3),A(4),B(3),B(4))
3251
3252 IF(ABS(MX(5)).GT.ZERO) TERM(18)=-MX(3)*MX(5)*HWDRM4(M12SQ,M13SQ,
3253
3254 & M(2),M(1),M(3),M(4),SM(5),SM(3),SW(5),SW(3),A(3),A(5),B(3),B(5))
3255
3256 IF(ABS(MX(6)).GT.ZERO) TERM(19)=-MX(3)*MX(6)*HWDRM4(M12SQ,M13SQ,
3257
3258 & M(2),M(1),M(3),M(4),SM(6),SM(3),SW(6),SW(3),A(3),A(6),B(3),B(6))
3259
3260 ENDIF
3261
3262 IF(ABS(MX(4)).GT.ZERO.AND.ND.NE.1) THEN
3263
3264 TERM(4) = MX(4)**2*HWDRM2(M13SQ,M(1),M(3),M(2),M(4),SM(4),
3265
3266 & SW(4),A(4),B(4))
3267
3268 IF(ABS(MX(5)).GT.ZERO) TERM(20)=-MX(4)*MX(5)*HWDRM4(M12SQ,M13SQ,
3269
3270 & M(2),M(1),M(3),M(4),SM(5),SM(4),SW(5),SW(4),A(4),A(5),B(4),B(5))
3271
3272 IF(ABS(MX(6)).GT.ZERO) TERM(21)=-MX(4)*MX(6)*HWDRM4(M12SQ,M13SQ,
3273
3274 & M(2),M(1),M(3),M(4),SM(6),SM(4),SW(6),SW(4),A(4),A(6),B(4),B(6))
3275
3276 ENDIF
3277
3278 IF(ABS(MX(5)).GT.ZERO.AND.ND.NE.2) THEN
3279
3280 TERM(5) = MX(5)**2*HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(5),
3281
3282 & SW(5),A(5),B(5))
3283
3284 IF(ABS(MX(6)).GT.ZERO) TERM(9)= MX(5)*MX(6)*HWDRM3(M12SQ,M(1),
3285
3286 & M(2),M(3),M(4),SM(5),SM(6),SW(5),SW(6),A(5),A(6),B(5),B(6))
3287
3288 ENDIF
3289
3290 IF(ABS(MX(6)).GT.ZERO.AND.ND.NE.2) TERM(6) = MX(6)**2*
3291
3292 & HWDRM2(M12SQ,M(1),M(2),M(3),M(4),SM(6),SW(6),A(6),B(6))
3293
3294 DO K=10,21
3295
3296 TERM(K)=TERM(K)*INFCOL
3297
3298 ENDDO
3299
3300C--Add them up
3301
3302 DO K=1,21
3303
3304 HWDRM1 = HWDRM1+TERM(K)
3305
3306 ENDDO
3307
3308C--Different colour flows for the gluino
3309
3310 IF(LM.NE.0) THEN
3311
3312 NPLN = 0.0D0
3313
3314 PLN = 0.0D0
3315
3316 DO K=1,9
3317
3318 PLN = PLN+TERM(K)
3319
3320 ENDDO
3321
3322 DO K=10,21
3323
3324 NPLN= NPLN+TERM(K)
3325
3326 ENDDO
3327
3328 DO K=1,3
3329
3330 TEST(K) = (TERM(2*K-1)+TERM(2*K)+TERM(6+K))*(1+NPLN/PLN)
3331
3332 ENDDO
3333
3334 ELSE
3335
3336 DO K=1,3
3337
3338 TEST(K) = 0.0D0
3339
3340 ENDDO
3341
3342 ENDIF
3343
3344 IF(TEST(4).LT.ZERO) CALL HWWARN('HWDRM1',50,*999)
3345
3346 999 END
3347
3348CDECK ID>, HWDRM2.
3349
3350*CMZ :- -20/07/99 10:56:12 by Peter Richardson
3351
3352*-- Author : Peter Richardson
3353
3354C-----------------------------------------------------------------------
3355
3356 FUNCTION HWDRM2(X,MA,MB,MC,MD,MR1,GAM1,A,B)
3357
3358C-----------------------------------------------------------------------
3359
3360C Function to compute the matrix element squared part of a 3-body
3361
3362C R-parity decay
3363
3364C-----------------------------------------------------------------------
3365
3366 IMPLICIT NONE
3367
3368 DOUBLE PRECISION X,MA,MB,MC,MD,A,B,HWDRM2,MR1,GAM1
3369
3370 HWDRM2 = (X - MA**2 - MB**2)*(4*A*B*MC*MD +
3371
3372 & (A**2 + B**2)*(-X + MC**2 + MD**2))/
3373
3374 & ((X-MR1**2)**2+GAM1**2*MR1**2)
3375
3376 END