CDECK ID>, HWHIGZ. *CMZ :- -02/05/91 11.18.44 by Federico Carminati *-- Author : Mike Seymour C----------------------------------------------------------------------- SUBROUTINE HWHIGZ C----------------------------------------------------------------------- C HIGGS PRODUCTION VIA THE BJORKEN PROCESS: E+E- --> Z(*) --> Z(*)H C WHERE ONE OR BOTH OF THE Zs IS OFF-SHELL C USES ALGORITHM OF BERENDS AND KLEISS: NUCL.PHYS. B260(1985)32 C C MEAN EVWGT = CROSS-SECTION (IN NB) * HIGGS BRANCHING FRACTION C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWUAEM,HWHIGY,HWRUNI,HWR,HWULDO,EMZ,CVE,CAE, & POL1,POL2,CE1,CE2,CE3,PMAX,EMZ2,S,B,FACTR,EMH,EMFAC,EMH2,A,XP, & CV,CA,BRHIGQ,BR,X1,X2,FAC1,FAC2,XPP,XPPSQ,COEF,X,XSQ,PROB,C1,C2, & CHIGG,PTHETA,SHIGG,C3,PHIMAX,CPHI,SPHI,C2PHI,S2PHI,PCM,ELST INTEGER IDEC,I,NLOOP,ICMF,IHIG,IZED,IFER,IANT,ID1,ID2,IN1,IN2 EXTERNAL HWUAEM,HWHIGY,HWRUNI,HWR,HWULDO SAVE CVE,CAE,CE1,CE2,CE3,PMAX,EMZ2,S,EMH,B,FACTR,A,EMH2 EQUIVALENCE (EMZ,RMASS(200)) DATA ELST/0/ C---SET UP CONSTANTS IN1=1 IF (JDAHEP(1,IN1).NE.0) IN1=JDAHEP(1,IN1) IN2=2 IF (JDAHEP(1,IN2).NE.0) IN2=JDAHEP(1,IN2) IF (FSTWGT.OR.ELST.NE.PHEP(5,3)) THEN ELST=PHEP(5,3) CVE=VFCH(11,1) CAE=AFCH(11,1) POL1=1.-EPOLN(3)*PPOLN(3) POL2=EPOLN(3)-PPOLN(3) CE1=(POL1*(CVE**2+CAE**2)+POL2*2.*CVE*CAE) CE2=(POL1*2.*CVE*CAE+POL2*(CVE**2+CAE**2)) IF ((IDHW(IN1).GT.IDHW(IN2).AND.PHEP(3,IN1).LT.ZERO).OR. & (IDHW(IN2).GT.IDHW(IN1).AND.PHEP(3,IN2).LT.ZERO)) CE2=-CE2 IF (TPOL) CE3=(CVE**2-CAE**2) PMAX=4 EMZ2=EMZ**2 S=PHEP(5,3)**2 B=EMZ*GAMZ/S FACTR=GEV2NB*CE1*(HWUAEM(RMASS(201)**2)*ENHANC(11))**2 & /(12.*S*SWEIN*(1.-SWEIN))*B/((1-EMZ2/S)**2+B**2) ENDIF IF (.NOT.GENEV) THEN C---CHOOSE HIGGS MASS, AND CALCULATE EVENT WEIGHT EVWGT=0D0 CALL HWHIGM(EMH,EMFAC) IF (EMH.LE.ZERO .OR. EMH.GT.PHEP(5,3)) RETURN EMSCA=EMH EMH2=EMH**2 A=4*EMH2/S XP=1+(EMH2-EMZ2)/S EVWGT=FACTR*HWHIGY(A,B,XP)*EMFAC C---INCLUDE BRANCHING RATIO OF HIGGS IDEC=MOD(IPROC,100) IF (IDEC.GT.0.AND.IDEC.LE.12) EVWGT=EVWGT*BRHIG(IDEC) IF (IDEC.EQ.0) THEN BRHIGQ=0 DO 10 I=1,6 10 BRHIGQ=BRHIGQ+BRHIG(I) EVWGT=EVWGT*BRHIGQ ENDIF C Add Z branching fractions CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,0) EVWGT=EVWGT*BR IF (IDEC.EQ.10) THEN CALL HWDBOZ(198,ID1,ID2,CV,CA,BR,1) CALL HWDBOZ(199,ID1,ID2,CV,CA,BR,1) EVWGT=EVWGT*BR ELSEIF (IDEC.EQ.11) THEN CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1) CALL HWDBOZ(200,ID1,ID2,CV,CA,BR,1) EVWGT=EVWGT*BR ENDIF ELSE C---GENERATE EVENT ICMF=NHEP+1 IHIG=NHEP+2 IZED=NHEP+3 IFER=NHEP+4 IANT=NHEP+5 CALL HWVEQU(5,PHEP(1,3),PHEP(1,ICMF)) NHEP=NHEP+5 C---CHOOSE ENERGY FRACTION OF HIGGS X1=SQRT(A) X2=1+0.25*A XP=1+(EMH2-EMZ2)/S FAC1=ATAN((X1-XP)/B) FAC2=ATAN((X2-XP)/B) XPP=MIN(X2,MAX(X1+B,XP)) XPPSQ=XPP**2 NLOOP=0 COEF=1./((12+2*A-12*XPP+XPPSQ)*SQRT(XPPSQ-A)) 20 NLOOP=NLOOP+1 IF (NLOOP.GT.NBTRY) CALL HWWARN('HWHIGZ',101,*999) X=XP+B*TAN(HWRUNI(1,FAC1,FAC2)) XSQ=X**2 PROB=COEF*((12+2*A-12*X+XSQ)*SQRT(XSQ-A)) IF (PROB.GT.PMAX) THEN PMAX=1.1*PROB CALL HWWARN('HWHIGZ',1,*999) WRITE (6,21) PMAX 21 FORMAT(7X,'NEW HWHIGZ MAX WEIGHT =',F8.4) ENDIF IF (PROB.LT.PMAX*HWR()) GOTO 20 C Choose Z decay mode CALL HWDBOZ(200,IDHW(IFER),IDHW(IANT),CV,CA,BR,0) C1=CE1*(CV**2+CA**2) C2=CE2*2.*CV*CA C---CHOOSE HIGGS DIRECTION C First polar angle NLOOP=0 COEF=(XSQ-A)/(8.*(1.-X)+XSQ+A) 30 NLOOP=NLOOP+1 IF (NLOOP.GT.NBTRY) CALL HWWARN('HWHIGZ',102,*999) CHIGG=HWRUNI(2,-ONE, ONE) PTHETA=1-COEF*CHIGG**2 IF (PTHETA.LT.HWR()) GOTO 30 SHIGG=SQRT(1-CHIGG**2) C Now azimuthal angle IF (TPOL) THEN C3=CE3*(CV*2+CA**2) COEF=COEF*SHIGG**2*C3/C1 PHIMAX=PTHETA+ABS(COEF) 40 CALL HWRAZM(ONE,CPHI,SPHI) C2PHI=2.*CPHI**2-1. S2PHI=2.*CPHI*SPHI PROB=PTHETA-COEF*(C2PHI*COSS+S2PHI*SINS) IF (PROB.LT.HWR()*PHIMAX) GOTO 40 ELSE CALL HWRAZM(ONE,CPHI,SPHI) ENDIF C Construct Higgs and Z momenta PHEP(5,IHIG)=EMH PHEP(4,IHIG)=X*PHEP(5,ICMF)/2 PCM=SQRT(PHEP(4,IHIG)**2-EMH2) PHEP(3,IHIG)=CHIGG*PCM PHEP(1,IHIG)=SHIGG*PCM*CPHI PHEP(2,IHIG)=SHIGG*PCM*SPHI CALL HWVDIF(4,PHEP(1,ICMF),PHEP(1,IHIG),PHEP(1,IZED)) CALL HWUMAS(PHEP(1,IZED)) C Choose orientation of Z decay NLOOP=0 COEF=2.*(C1+ABS(C2))*HWULDO(PHEP(1,IN1),PHEP(1,IZED)) & *HWULDO(PHEP(1,IN2),PHEP(1,IZED))/S IF (TPOL) COEF=COEF*(C1+ABS(C2)+ABS(C3))/(C1+ABS(C2)) PCM=PHEP(5,IZED)/2 PHEP(5,IFER)=0 PHEP(5,IANT)=0 50 NLOOP=NLOOP+1 IF (NLOOP.GT.NBTRY) CALL HWWARN('HWHIGZ',103,*999) CALL HWDTWO(PHEP(1,IZED),PHEP(1,IFER),PHEP(1,IANT), & PCM,TWO,.TRUE.) PROB=C1*(PHEP(4,IFER)*PHEP(4,IANT)-PHEP(3,IFER)*PHEP(3,IANT)) & +C2*(PHEP(4,IFER)*PHEP(3,IANT)-PHEP(3,IFER)*PHEP(4,IANT)) IF (TPOL) PROB=PROB+C3* & (COSS*(PHEP(1,IFER)*PHEP(1,IANT)-PHEP(2,IFER)*PHEP(2,IANT)) & +SINS*(PHEP(1,IFER)*PHEP(2,IANT)+PHEP(2,IFER)*PHEP(1,IANT))) IF (PROB.LT.HWR()*COEF) GOTO 50 C---SET UP STATUS CODES, ISTHEP(ICMF)=120 ISTHEP(IHIG)=190 ISTHEP(IZED)=195 ISTHEP(IFER)=113 ISTHEP(IANT)=114 C---COLOR CONNECTIONS, JMOHEP(1,ICMF)=1 JMOHEP(2,ICMF)=2 JDAHEP(1,ICMF)=IHIG JDAHEP(2,ICMF)=IZED JMOHEP(1,IHIG)=ICMF JMOHEP(1,IZED)=ICMF JMOHEP(1,IFER)=IZED JMOHEP(1,IANT)=IZED JMOHEP(2,IFER)=IANT JMOHEP(2,IANT)=IFER JDAHEP(1,IZED)=IFER JDAHEP(2,IZED)=IANT JDAHEP(2,IFER)=IANT JDAHEP(2,IANT)=IFER C---IDENTITY CODES IDHW(ICMF)=200 IDHW(IHIG)=201 IDHW(IZED)=200 IDHEP(ICMF)=IDPDG(IDHW(ICMF)) IDHEP(IHIG)=IDPDG(IDHW(IHIG)) IDHEP(IZED)=IDPDG(IDHW(IZED)) IDHEP(IFER)=IDPDG(IDHW(IFER)) IDHEP(IANT)=IDPDG(IDHW(IANT)) ENDIF 999 END CDECK ID>, HWHPH2. *CMZ :- -12/01/93 10.12.43 by Bryan Webber *-- Author : Ian Knowles C----------------------------------------------------------------------- SUBROUTINE HWHPH2 C----------------------------------------------------------------------- C QQD direct photon pair production: mean EVWGT = sigma in nb C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWRUNI,HWUALF,HWHPPB,EPS,RCS,ET,EJ,KK,KK2, & YJ1INF,YJ1SUP,Z1,YJ2INF,YJ2SUP,Z2,FACT,FACTR,RS,S,T,U,CSTU,TQSQ, & DSTU,HCS INTEGER ID,ID1,ID2 EXTERNAL HWR,HWRUNI,HWUALF,HWHPPB SAVE HCS,CSTU,DSTU,FACT PARAMETER (EPS=1.D-9) IF (GENEV) THEN RCS=HCS*HWR() ELSE EVWGT=0. CALL HWRPOW(ET,EJ) KK=ET/PHEP(5,3) KK2=KK**2 IF (KK.GE.ONE) RETURN YJ1INF=MAX( YJMIN , LOG((1.-SQRT(1.-KK2))/KK) ) YJ1SUP=MIN( YJMAX , LOG((1.+SQRT(1.-KK2))/KK) ) IF (YJ1INF.GE.YJ1SUP) RETURN Z1=EXP(HWRUNI(1,YJ1INF,YJ1SUP)) YJ2INF=MAX( YJMIN , -LOG(2./KK-1./Z1) ) YJ2SUP=MIN( YJMAX , LOG(2./KK-Z1) ) IF (YJ2INF.GE.YJ2SUP) RETURN Z2=EXP(HWRUNI(2,YJ2INF,YJ2SUP)) XX(1)=0.5*(Z1+Z2)*KK IF (XX(1).GE.ONE) RETURN XX(2)=XX(1)/(Z1*Z2) IF (XX(2).GE.ONE) RETURN COSTH=(Z1-Z2)/(Z1+Z2) S=XX(1)*XX(2)*PHEP(5,3)**2 RS=0.5*SQRT(S) T=-0.5*S*(1.-COSTH) U=-S-T EMSCA=SQRT(2.*S*T*U/(S*S+T*T+U*U)) FACT=GEV2NB*PIFAC*0.5*ET*EJ*(YJ1SUP-YJ1INF)*(YJ2SUP-YJ2INF) & *(ALPHEM/S)**2 CALL HWSGEN(.FALSE.) CSTU=2.*(U/T+T/U)/CAFAC IF (DISF(13,1).GT.EPS.AND.DISF(13,2).GT.EPS) THEN TQSQ=0. DO 10 ID=1,6 10 IF (RMASS(ID).LT.RS) TQSQ=TQSQ+QFCH(ID)**2 DSTU=DISF(13,1)*DISF(13,2)*FACT*HWHPPB(S,T,U) & /64.*(HWUALF(1,EMSCA)*TQSQ/PIFAC)**2 ENDIF ENDIF HCS=0. DO 30 ID=1,6 FACTR=FACT*CSTU*QFCH(ID)**4 C q+qbar ---> gamma+gamma ID1=ID ID2=ID+6 IF (DISF(ID1,1).LT.EPS.OR.DISF(ID2,2).LT.EPS) GOTO 20 HCS=HCS+FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(59,59,2134,61,*99) C qbar+q ---> gamma+gamma 20 ID1=ID+6 ID2=ID IF (DISF(ID1,1).LT.EPS.OR.DISF(ID2,2).LT.EPS) GOTO 30 HCS=HCS+FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(59,59,2134,62,*99) 30 CONTINUE C g+g ---> gamma+gamma ID1=13 ID2=13 HCS=HCS+DSTU IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(59,59,2134,63,*99) EVWGT=HCS RETURN C Generate event 99 IDN(1)=ID1 IDN(2)=ID2 IDCMF=15 CALL HWETWO 999 END CDECK ID>, HWHPHO. *CMZ :- -26/04/91 14.55.45 by Federico Carminati *-- Author : Bryan Webber C----------------------------------------------------------------------- SUBROUTINE HWHPHO C----------------------------------------------------------------------- C QCD DIRECT PHOTON + JET PRODUCTION: MEAN EVWGT = SIGMA IN NB C----------------------------------------------------------------------- INCLUDE 'HERWIG61.INC' DOUBLE PRECISION HWR,HWRUNI,HWUALF,HWHPPB,EPS,RCS,ET,EJ,KK,KK2, & YJ1INF,YJ1SUP,Z1,YJ2INF,YJ2SUP,Z2,FACT,FACTR,FACTF,RS,S,T,U,CF, & AF,CSTU,CTSU,CUST,DSTU,HCS,TQCH INTEGER ID,ID1,ID2 EXTERNAL HWR,HWRUNI,HWUALF,HWHPPB SAVE HCS PARAMETER (EPS=1.D-9) IF (GENEV) THEN RCS=HCS*HWR() ELSE EVWGT=0. CALL HWRPOW(ET,EJ) KK=ET/PHEP(5,3) KK2=KK**2 IF (KK.GE.ONE) RETURN YJ1INF=MAX( YJMIN , LOG((1.-SQRT(1.-KK2))/KK) ) YJ1SUP=MIN( YJMAX , LOG((1.+SQRT(1.-KK2))/KK) ) IF (YJ1INF.GE.YJ1SUP) RETURN Z1=EXP(HWRUNI(1,YJ1INF,YJ1SUP)) YJ2INF=MAX( YJMIN , -LOG(2./KK-1./Z1) ) YJ2SUP=MIN( YJMAX , LOG(2./KK-Z1) ) IF (YJ2INF.GE.YJ2SUP) RETURN Z2=EXP(HWRUNI(2,YJ2INF,YJ2SUP)) XX(1)=0.5*(Z1+Z2)*KK IF (XX(1).GE.ONE) RETURN XX(2)=XX(1)/(Z1*Z2) IF (XX(2).GE.ONE) RETURN COSTH=(Z1-Z2)/(Z1+Z2) S=XX(1)*XX(2)*PHEP(5,3)**2 RS=0.5*SQRT(S) T=-0.5*S*(1.-COSTH) U=-S-T C---SET EMSCA TO HARD PROCESS SCALE (APPROX ET-JET) EMSCA=SQRT(2.*S*T*U/(S*S+T*T+U*U)) FACT=GEV2NB*PIFAC*0.5*ET*EJ*ALPHEM & *HWUALF(1,EMSCA)*(YJ1SUP-YJ1INF)*(YJ2SUP-YJ2INF)/S**2 CALL HWSGEN(.FALSE.) C CF=2.*CFFAC/CAFAC AF=-1./CAFAC CSTU=CF*(U/T+T/U) CTSU=AF*(U/S+S/U) CUST=AF*(T/S+S/T) IF (DISF(13,1).GT.EPS.AND.DISF(13,2).GT.EPS) THEN TQCH=0. DO 10 ID=1,6 10 IF (RMASS(ID).LT.RS) TQCH=TQCH+QFCH(ID) DSTU=DISF(13,1)*DISF(13,2)*FACT*HWHPPB(S,T,U) & *5./768.*(HWUALF(1,EMSCA)*TQCH/PIFAC)**2 ENDIF ENDIF C HCS=0. DO 30 ID=1,6 FACTR=FACT*QFCH(ID)**2 C---QUARK FIRST ID1=ID IF (DISF(ID1,1).LT.EPS) GOTO 20 ID2=ID1+6 HCS=HCS+CSTU*FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP( 13, 59,2314,41,*9) ID2=13 HCS=HCS+CTSU*FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1, 59,3124,42,*9) C---QBAR FIRST 20 ID1=ID+6 IF (DISF(ID1,1).LT.EPS) GOTO 30 ID2=ID HCS=HCS+CSTU*FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP( 13, 59,3124,43,*9) ID2=13 HCS=HCS+CTSU*FACTR*DISF(ID1,1)*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID1, 59,2314,44,*9) 30 CONTINUE C---GLUON FIRST ID1=13 FACTF=FACT*CUST*DISF(ID1,1) DO 50 ID=1,6 FACTR=FACTF*QFCH(ID)**2 ID2=ID IF (DISF(ID2,2).LT.EPS) GOTO 40 HCS=HCS+FACTR*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID2, 59,2314,45,*9) 40 ID2=ID+6 IF (DISF(ID2,2).LT.EPS) GOTO 50 HCS=HCS+FACTR*DISF(ID2,2) IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID2, 59,3124,46,*9) 50 CONTINUE C g+g ---> g+gamma ID2=13 HCS=HCS+DSTU IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP( 13, 59,2314,47,*9) EVWGT=HCS RETURN C---GENERATE EVENT 9 IDN(1)=ID1 IDN(2)=ID2 IDCMF=15 CALL HWETWO 999 END CDECK ID>, HWHPPB. *CMZ :- -12/01/93 10.12.43 by Bryan Webber *-- Author : Ian Knowles C----------------------------------------------------------------------- FUNCTION HWHPPB(S,T,U) C----------------------------------------------------------------------- C Quark box diagram contribution to photon/gluon scattering C Internal quark mass neglected: m_q << U,T,S C----------------------------------------------------------------------- IMPLICIT NONE DOUBLE PRECISION HWHPPB,S,T,U,S2,T2,U2,PI2,ALNTU,ALNST,ALNSU PI2=ACOS(-1.D0)**2 S2=S**2 T2=T**2 U2=U**2 ALNTU=LOG(T/U) ALNST=LOG(-S/T) ALNSU=ALNST+ALNTU HWHPPB=5.*4. & +((2.*S2+2.*(U2-T2)*ALNTU+(T2+U2)*(ALNTU**2+PI2))/S2)**2 & +((2.*U2+2.*(T2-S2)*ALNST+(T2+S2)* ALNST**2 )/U2)**2 & +((2.*T2+2.*(U2-S2)*ALNSU+(U2+S2)* ALNSU**2 )/T2)**2 & +4.*PI2*(((T2-S2+(T2+S2)*ALNST)/U2)**2 & +((U2-S2+(U2+S2)*ALNSU)/T2)**2) END