4 *CMZ :- -02/05/91 10.58.29 by Federico Carminati
6 *-- Author : Zoltan Kunszt, modified by Bryan Webber
8 C-----------------------------------------------------------------------
12 C-----------------------------------------------------------------------
14 C E+E- -> W+W-/Z0Z0 (BASED ON ZOLTAN KUNSZT'S PROGRAM)
16 C-----------------------------------------------------------------------
18 INCLUDE 'HERWIG61.INC'
22 DOUBLE PRECISION HWUAEM,HWR,HWUPCM,ETOT,STOT,FLUXW,GAMM,GIMM,
24 & WM2,WXMIN,WX1MAX,WX2MAX,FJAC1,FJAC2,WX1,WX2,WMM1,WMM2,XXM,W2BO,
26 & PST,WEIGHT,TOTSIG,WMASS,WWIDTH,ELST,CV,CA,BR,XMASS,PLAB,PRW,PCM,
28 & AMPWW(4),CCC,HELSUM,HELCTY,BRZED(12),BRTOT,CPFAC,CPALL,RLL(12),
32 INTEGER IB,IBOS,I,ID1,ID2,NTRY,IDP(10),IDBOS(2),J1,J2,IPRC,ILST,
34 & IDZOLT(16),MAP(12),NEWHEP
38 EXTERNAL HWUAEM,HWR,HWUPCM
40 SAVE IDP,STOT,FLUXW,GAMM,GIMM,WM2,WXMIN,WX1MAX,FJAC1,ELST,ILST,
44 COMMON/HWHEWP/XMASS(10),PLAB(5,10),PRW(5,2),PCM(5,10)
46 COMMON/HWHEWQ/ZH(7,7),ZCH(7,7),ZD(7,7)
48 COMMON/HWHEWR/CPFAC(12,12,8),CPALL(8)
52 DATA IDZOLT/4,3,8,7,12,11,4*0,2,1,6,5,10,9/
54 DATA MAP/12,11,2,1,14,13,4,3,16,15,6,5/
56 IF (IERROR.NE.0) RETURN
58 EISBM1=IDHW(1).LT.IDHW(2)
70 CALL HWVEQU(5,PRW(1,IB),PHEP(1,IBOS))
72 IF (EISBM1) PHEP(3,IBOS)=-PHEP(3,IBOS)
74 CALL HWVZRO(4,VHEP(1,IBOS))
76 CALL HWUDKL(IDBOS(IB),PHEP(1,IBOS),DIST)
78 CALL HWVSUM(4,VHEP(1,IBOS),DIST,DIST)
82 IDHEP(IBOS)=IDPDG(IDBOS(IB))
92 CALL HWVEQU(5,PLAB(1,2*IB+I),PHEP(1,NHEP+I))
94 IF (EISBM1) PHEP(3,NHEP+I)=-PHEP(3,NHEP+I)
96 CALL HWVEQU(4,DIST,VHEP(1,NHEP+I))
98 C---STATUS, IDs AND POINTERS
102 IDHW(NHEP+I)=IDP(2*IB+I)
104 IDHEP(NHEP+I)=IDPDG(IDP(2*IB+I))
106 JDAHEP(I,IBOS)=NHEP+I
108 JMOHEP(1,NHEP+I)=IBOS
110 JMOHEP(2,NHEP+I)=JMOHEP(1,IBOS)
116 JMOHEP(2,NHEP)=NHEP-1
118 JDAHEP(2,NHEP)=NHEP-1
120 JMOHEP(2,NHEP-1)=NHEP
122 JDAHEP(2,NHEP-1)=NHEP
134 IF (ETOT.NE.ELST .OR. IPRC.NE.ILST) THEN
138 FLUXW=GEV2NB*.125*(HWUAEM(STOT)/PIFAC)**4/STOT
150 ELSEIF (IPRC.EQ.50) THEN
160 C---LOAD FERMION COUPLINGS TO Z
164 RLL(I)=VFCH(MAP(I),1)+AFCH(MAP(I),1)
166 RRL(I)=VFCH(MAP(I),1)-AFCH(MAP(I),1)
184 IF (MOD(J1-1,4).GE.2) CCC=CCC*CAFAC
186 IF (MOD(J2-1,4).GE.2) CCC=CCC*CAFAC
188 CPFAC(J1,J2,1)=CCC*(RLL(2)**2*RLL(J1)*RLL(J2))**2
190 CPFAC(J1,J2,2)=CCC*(RLL(2)**2*RRL(J1)*RLL(J2))**2
192 CPFAC(J1,J2,3)=CCC*(RLL(2)**2*RLL(J1)*RRL(J2))**2
194 CPFAC(J1,J2,4)=CCC*(RLL(2)**2*RRL(J1)*RRL(J2))**2
196 CPFAC(J1,J2,5)=CCC*(RRL(2)**2*RLL(J1)*RLL(J2))**2
198 CPFAC(J1,J2,6)=CCC*(RRL(2)**2*RRL(J1)*RLL(J2))**2
200 CPFAC(J1,J2,7)=CCC*(RRL(2)**2*RLL(J1)*RRL(J2))**2
202 CPFAC(J1,J2,8)=CCC*(RRL(2)**2*RRL(J1)*RRL(J2))**2
206 IF (J1.EQ.1.AND.J2.EQ.1) CPALL(I)=0
208 CPALL(I)=CPALL(I)+CPFAC(J1,J2,I)
210 BRZED(J1)=BRZED(J1)+CPFAC(J1,J2,I)
212 BRTOT=BRTOT+CPFAC(J1,J2,I)
222 70 BRZED(I)=BRZED(I)/BRTOT
226 CALL HWWARN('HWHEWW',500,*999)
236 WXMIN=ATAN(-WMASS/WWIDTH)
238 WX1MAX=ATAN((STOT-WM2)*GIMM)
252 WX1=WXMIN+FJAC1*HWR()
254 WMM1=GAMM*TAN(WX1)+WM2
258 WX2MAX=ATAN(((ETOT-XMASS(1))**2-WM2)*GIMM)
262 WX2=WXMIN+FJAC2*HWR()
264 WMM2=GAMM*TAN(WX2)+WM2
268 IF (HWRLOG(HALF))THEN
278 C---CTMAX=ANGULAR CUT ON COS W-ANGLE
280 CALL HWHEW0(1,ETOT,XMASS(1),PRW(1,1),W2BO,CTMAX)
282 IF (W2BO.EQ.ZERO) RETURN
284 C---FOR ZZ EVENTS, FORCE BOSE STATISTICS, BY KILLING EVENTS WITH COS1<0
288 IF (PRW(3,1).LT.ZERO) RETURN
290 C---AND THEN SYMMETRIZE (THIS PROCEDURE VASTLY IMPROVES EFFICIENCY)
292 IF (HWRLOG(HALF)) THEN
312 C---LET THE W BOSONS DECAY
320 CALL HWDBOZ(IDBOS(IB),ID1,ID2,CV,CA,BR,1)
322 PST=HWUPCM(XMASS(IB),RMASS(ID1),RMASS(ID2))
324 IF (PST.LT.ZERO) THEN
326 CALL HWDBOZ(IDBOS(IB),ID1,ID2,CV,CA,BR,2)
328 IF (NTRY.LE.NBTRY) GOTO 80
330 C CALL HWWARN('HWHEWW',1,*999)
342 PLAB(5,2*IB+1)=RMASS(ID1)
344 PLAB(5,2*IB+2)=RMASS(ID2)
346 CALL HWDTWO(PRW(1,IB),PLAB(1,2*IB+1),PLAB(1,2*IB+2),
352 WEIGHT=FLUXW*W2BO*FJAC1*FJAC2*(0.5D0*PIFAC*GIMM)**2
356 CALL HWHEW2(6,PCM(1,1),ZH,ZCH,ZD)
360 CALL HWHEW3(5,6,3,4,1,2,AMPWW)
362 TOTSIG=9.*AMPWW(1)+6.*(AMPWW(2)+AMPWW(3))+4.*AMPWW(4)
364 EVWGT=TOTSIG*WEIGHT*BR
368 ID1=IDZOLT(IDPDG(IDP(3)))
370 ID2=IDZOLT(IDPDG(IDP(5)))
372 CALL HWHEW5(5,6,3,4,1,2,HELSUM,HELCTY,ID1,ID2)
374 EVWGT=HELCTY*WEIGHT*BR/(BRZED(ID1)*BRZED(ID2))
384 *CMZ :- -18/05/99 14.55.44 by Kosuke Odagiri
386 *-- Author : Bryan Webber
388 C-----------------------------------------------------------------------
392 C-----------------------------------------------------------------------
394 C QCD HEAVY FLAVOUR PRODUCTION: MEAN EVWGT = SIGMA IN NB
396 C-----------------------------------------------------------------------
398 INCLUDE 'HERWIG61.INC'
400 DOUBLE PRECISION HWR,HWRUNI,HWUALF,EPS,RCS,Z1,Z2,ET,EJ,
402 & QM2,QPE,FACTR,S,T,U,ST,TU,US,TUS,UST,EN,RN,AF,ASTU,
404 & AUST,CF,CN,CS,CSTU,CSUT,CTSU,CTUS,HCS,UT,SU,GT,DIST,KK,KK2,
406 & YJ1INF,YJ1SUP,YJ2INF,YJ2SUP
408 INTEGER IQ1,IQ2,ID1,ID2
412 EXTERNAL HWR,HWRUNI,HWUALF
414 SAVE HCS,ASTU,AUST,CSTU,CSUT,CTSU,CTUS,S,T,TU,U,US
416 PARAMETER (EPS=1.D-9)
432 IF (KK.GE.ONE) RETURN
434 YJ1INF = MAX( YJMIN, LOG((ONE-SQRT(ONE-KK2))/KK) )
436 YJ1SUP = MIN( YJMAX, LOG((ONE+SQRT(ONE-KK2))/KK) )
438 IF (YJ1INF.GE.YJ1SUP) RETURN
440 Z1=EXP(HWRUNI(1,YJ1INF,YJ1SUP))
442 YJ2INF = MAX( YJMIN, -LOG(TWO/KK-ONE/Z1) )
444 YJ2SUP = MIN( YJMAX, LOG(TWO/KK-Z1) )
446 IF (YJ2INF.GE.YJ2SUP) RETURN
448 Z2=EXP(HWRUNI(2,YJ2INF,YJ2SUP))
450 XX(1)=HALF*(Z1+Z2)*KK
452 IF (XX(1).GE.ONE) RETURN
456 IF (XX(2).GE.ONE) RETURN
458 S=XX(1)*XX(2)*PHEP(5,3)**2
466 IF (QPE.LE.ZERO) RETURN
468 COSTH=HALF*ET*(Z1-Z2)/SQRT(Z1*Z2*QPE)
470 IF (ABS(COSTH).GT.ONE) RETURN
472 C---REDEFINE S, T, U AS P1.P2, -P1.P3, -P1.P4
476 T=-HALF*(1.+Z2/Z1)*(HALF*ET)**2
480 C---SET EMSCA TO HEAVY HARD PROCESS SCALE
482 EMSCA=SQRT(4.*S*T*U/(S*S+T*T+U*U))
484 FACTR = GEV2NB*.125*PIFAC*EJ*ET*(HWUALF(1,EMSCA)/S)**2
486 & *(YJ1SUP-YJ1INF)*(YJ2SUP-YJ2INF)
514 ASTU=AF*(1.-2.*UST+QM2/T)
516 AUST=AF*(1.-2.*TUS+QM2/S)
522 C-----------------------------------------------------------------------
524 C---Heavy flavour colour decomposition modifications below (KO)
526 C-----------------------------------------------------------------------
528 CS=(TU+UT-CN/TUS)*(HALF-TUS+QM2/S-QM2**2/U/T/TWO)
530 CSTU=CF*CS/(ONE+TU**2)
532 CSUT=CF*CS/(ONE+UT**2)
534 CS=(SU+US-CN/UST)*(HALF-UST+QM2/T-QM2**2/U/S/TWO)
536 CTSU=-FACTR*CS/(ONE+SU**2)
538 CTUS=-FACTR*CS/(ONE+US**2)
540 C-----------------------------------------------------------------------
542 C CS=HALF/TU-QM2/T-HALF*(QM2/T)**2
544 C CSTU=CF*(CS- US**2-QM2/S - CN*(CS+QM2*QM2/(S*T)))
546 C CS=HALF*TU-QM2/U-HALF*(QM2/U)**2
548 C CSUT=CF*(CS-1./ST**2-QM2/S - CN*(CS+QM2*QM2/(S*U)))
550 C CS=HALF*US-QM2/S-HALF*(QM2/S)**2
552 C CTSU=-FACTR*(CS-1./TU**2-QM2/T - CN*(CS+QM2*QM2/(S*T)))
554 C CS=HALF/US-QM2/U-HALF*(QM2/U)**2
556 C CTUS=-FACTR*(CS- ST**2-QM2/T - CN*(CS+QM2*QM2/(T*U)))
558 C-----------------------------------------------------------------------
570 IF (DISF(ID1,1).LT.EPS) GOTO 6
572 HQ1=ID1.EQ.IQ1.OR.ID1.EQ.IQ2
576 IF (DISF(ID2,2).LT.EPS) GOTO 5
578 HQ2=ID2.EQ.IQ1.OR.ID2.EQ.IQ2
580 DIST=DISF(ID1,1)*DISF(ID2,2)
584 C---PROCESSES INVOLVING HEAVY CONSTITUENT
586 C N.B. NEGLECT CASE THAT BOTH ARE HEAVY
588 IF (HQ1.AND.HQ2) GOTO 5
598 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3421, 3,*9)
600 ELSEIF (ID2.NE.13) THEN
604 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3142, 9,*9)
610 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3142,10,*9)
614 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3421,11,*9)
618 ELSEIF (ID1.NE.13) THEN
626 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,2413,17,*9)
628 ELSEIF (ID2.NE.13) THEN
632 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,4312,20,*9)
638 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,2413,21,*9)
642 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,4312,22,*9)
654 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,2413,23,*9)
658 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3421,24,*9)
660 ELSEIF (ID2.LT.13) THEN
664 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,3142,25,*9)
668 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1,ID2,4312,26,*9)
674 ELSEIF (ID2.NE.13.AND.ID2.EQ.ID1+6) THEN
676 C---LIGHT Q-QBAR ANNIHILATION
680 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(IQ1,IQ2,2413, 4,*9)
682 ELSEIF (ID1.NE.13.AND.ID1.EQ.ID2+6) THEN
684 C---LIGHT QBAR-Q ANNIHILATION
688 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(IQ2,IQ1,3142,12,*9)
690 ELSEIF (ID1.EQ.13.AND.ID2.EQ.13) THEN
696 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(IQ1,IQ2,2413,27,*9)
700 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(IQ1,IQ2,4123,28,*9)
724 C Calculate coefficients for constructing spin density matrices
726 IF (IHPRO.EQ.7 .OR.IHPRO.EQ.8 .OR.
728 & IHPRO.EQ.15.OR.IHPRO.EQ.16) THEN
730 C qqbar-->gg or qbarq-->gg
748 ELSEIF (IHPRO.EQ.10.OR.IHPRO.EQ.11.OR.
750 & IHPRO.EQ.21.OR.IHPRO.EQ.22.OR.
752 & IHPRO.EQ.23.OR.IHPRO.EQ.24.OR.
754 & IHPRO.EQ.25.OR.IHPRO.EQ.26) THEN
756 C qg-->qg or qbarg-->qbarg or gq-->gq or gqbar-->gqbar
774 ELSEIF (IHPRO.EQ.27.OR.IHPRO.EQ.28) THEN
794 ELSEIF (IHPRO.EQ.29.OR.IHPRO.EQ.30.OR.
808 GCOEF(1)=GT*GT-GCOEF(2)-GCOEF(3)-GCOEF(4)
810 GCOEF(5)=GT*(GT-2.*S*S)-GCOEF(2)
812 GCOEF(6)=GT*(GT-2.*T*T)-GCOEF(3)
814 GCOEF(7)=GT*(GT-2.*U*U)-GCOEF(4)
828 *CMZ :- -23/08/94 13.22.29 by Mike Seymour
830 *-- Author : Ulrich Baur & Nigel Glover, adapted by Ian Knowles
832 C-----------------------------------------------------------------------
834 FUNCTION HWHIG1(S,T,U,EH2,EQ2,I,J,K,I1,J1,K1)
836 C-----------------------------------------------------------------------
838 C Basic matrix elements for Higgs + jet production; used in HWHIGA
840 C-----------------------------------------------------------------------
844 DOUBLE COMPLEX HWHIG1,HWHIG2,HWHIG5,BI(4),CI(7),DI(3)
846 DOUBLE PRECISION S,T,U,EH2,EQ2,S1,T1,U1,ONE,TWO,FOUR,HALF
848 INTEGER I,J,K,I1,J1,K1
850 COMMON/CINTS/BI,CI,DI
852 PARAMETER (ONE =1.D0, TWO =2.D0, FOUR =4.D0, HALF =0.5D0)
854 C-----------------------------------------------------------------------
856 C +++ helicity amplitude for: g+g --> g+H
858 C-----------------------------------------------------------------------
866 HWHIG1=EQ2*FOUR*DSQRT(TWO*S*T*U)*(
868 & -FOUR*(ONE/(U*T)+ONE/(U*U1)+ONE/(T*T1))
870 & -FOUR*((TWO*S+T)*BI(K)/U1**2+(TWO*S+U)*BI(J)/T1**2)/S
872 & -(S-FOUR*EQ2)*(S1*CI(I1)+(U-S)*CI(J1)+(T-S)*CI(K1))/(S*T*U)
874 & -8.D0*EQ2*(CI(J1)/(T*T1)+CI(K1)/(U*U1))
876 & +HALF*(S-FOUR*EQ2)*(S*T*DI(K)+U*S*DI(J)-U*T*DI(I))/(S*T*U)
880 & -TWO*(U*CI(K)+T*CI(J)+U1*CI(K1)+T1*CI(J1)-U*T*DI(I))/S**2 )
884 C-----------------------------------------------------------------------
886 ENTRY HWHIG2(S,T,U,EH2,EQ2,I,J,K,I1,J1,K1)
888 C-----------------------------------------------------------------------
890 C ++- helicity amplitude for: g+g --> g+H
892 C-----------------------------------------------------------------------
900 HWHIG2=EQ2*FOUR*DSQRT(TWO*S*T*U)*(FOUR*EH2
902 & +(EH2-FOUR*EQ2)*(S1*CI(4)+T1*CI(5)+U1*CI(6))
904 & -HALF*(EH2-FOUR*EQ2)*(S*T*DI(3)+U*S*DI(2)+U*T*DI(1)) )/(S*T*U)
908 C-----------------------------------------------------------------------
910 ENTRY HWHIG5(S,T,U,EH2,EQ2,I,J,K,I1,J1,K1)
912 C-----------------------------------------------------------------------
914 C Amplitude for: q+qbar --> g+H
916 C-----------------------------------------------------------------------
918 HWHIG5=DCMPLX(TWO+TWO*S/(S-EH2))*BI(I)+DCMPLX(FOUR*EQ2-U-T)*CI(K)