4 *CMZ :- -26/04/91 11.11.55 by Bryan Webber
6 *-- Author : Bryan Webber & Luca Stanco
8 C-----------------------------------------------------------------------
10 SUBROUTINE HWEGAM(IHEP,ZMI,ZMA,WWA)
12 C-----------------------------------------------------------------------
14 C GENERATES A PHOTON IN WEIZSACKER-WILLIAMS (WWA=.TRUE.) OR
16 C ELSE EQUIVALENT PHOTON APPROX FROM INCOMING E+, E-, MU+ OR MU-
18 C-----------------------------------------------------------------------
20 INCLUDE 'HERWIG61.INC'
22 DOUBLE PRECISION HWR,HWRUNI,EGMIN,ZMIN,ZMAX,ZGAM,SS,ZMI,ZMA,
24 & PPL,PMI,QT2,Q2,QQMIN,QQMAX,S0,A,RPM(2)
26 INTEGER IHEP,IHADIS,HQ,I
34 IF (IERROR.NE.0) RETURN
36 IF (IHEP.LT.1.OR.IHEP.GT.2) CALL HWWARN('HWEGAM',500,*999)
48 IF (JDAHEP(1,IHADIS).NE.0) IHADIS=JDAHEP(1,IHADIS)
52 C---DEFINE LIMITS FOR GAMMA MOMENTUM FRACTION
54 IF (ZMI.LE.ZERO .OR. ZMA.GT.ONE) THEN
56 IF (IPRO.EQ.13.OR.IPRO.EQ.14) THEN
60 ELSEIF(IPRO.EQ.15.OR.IPRO.EQ.18.OR.IPRO.EQ.22.OR.IPRO.EQ.24.OR.
62 & IPRO.EQ.50.OR.IPRO.EQ.53.OR.IPRO.EQ.55)THEN
66 ELSEIF (IPRO.EQ.17.OR.IPRO.EQ.51) THEN
70 S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
72 ELSEIF (IPRO.EQ.16.OR.IPRO.EQ.19.OR.IPRO.EQ.95) THEN
74 S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
76 ELSEIF (IPRO.EQ.23) THEN
78 S0 = MAX(2*RMASS(1),RMASS(201)-GAMMAX*GAMH)**2
80 S0 = (PTMIN+SQRT(PTMIN**2+S0))**2
82 ELSEIF (IPRO.EQ.20) THEN
86 ELSEIF (IPRO.EQ.21) THEN
88 S0 = (PTMIN+SQRT(PTMIN**2+RMASS(198)**2))**2
92 ELSEIF(IPRO.EQ.30) THEN
94 S0 = 4.0D0*(PTMIN**2+RMMNSS**2)
96 ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
104 IF(HQ.GE.10.AND.HQ.LT.20) THEN
106 RPM(1) = ABS(RMASS(450))
108 IF(HQ.GT.10) RPM(1) = ABS(RMASS(449+MOD(HQ,10)))
110 ELSEIF(HQ.GE.20.AND.HQ.LT.30) THEN
112 RPM(1) = ABS(RMASS(454))
114 IF(HQ.GT.20) RPM(1) = ABS(RMASS(453+MOD(HQ,20)))
116 ELSEIF(HQ.EQ.30) THEN
120 ELSEIF(HQ.EQ.40) THEN
128 RPM(1) = MIN(RPM(1),RMASS(425+I))
134 RPM(1) = MIN(RMASS(405),RMASS(406))
140 ELSEIF(HQ.EQ.50) THEN
148 RPM(1) = MIN(RPM(1),RMASS(425+I))
154 RPM(2) = MIN(RPM(1),RMASS(433+2*I))
158 RPM(1) = MIN(RPM(1),RPM(2))
164 RPM(2) = MIN(RPM(2),RMASS(204+I))
176 RPM(1) = MIN(RPM(1),RMASS(401+I))
178 RPM(2) = MIN(RPM(2),RMASS(413+I))
182 RPM(1) = MIN(RPM(1),RPM(2))
188 RPM(2) = MIN(RPM(2),RMASS(204+I))
198 RPM(2) = MIN(RPM(2),RMASS(204+I))
202 ELSEIF(HQ.GE.60) THEN
212 S0 = RPM(1)+RPM(2)+TWO*(PTMIN**2+
214 & SQRT(RPM(1)*RPM(2)+PTMIN**2*(RPM(1)+RPM(2)+PTMIN**2)))
218 ELSEIF (IPRO.EQ.52) THEN
222 S0 = (PTMIN+SQRT(PTMIN**2+RMASS(HQ)**2))**2
224 ELSEIF (IPRO.EQ.60) THEN
234 IF (HQ.GT.6) HQ=2*HQ+107
236 IF (HQ.EQ.127) HQ=198
238 S0 = 4.D0*(PTMIN**2+RMASS(HQ)**2)
242 ELSEIF (IPRO.EQ.80) THEN
246 ELSEIF (IPRO.EQ.90) THEN
250 ELSEIF (IPRO.EQ.91.OR.IPRO.EQ.92) THEN
252 S0 = Q2MIN+4.D0*PTMIN**2
256 IF (HQ.GT.0) S0 = S0+4.D0*RMASS(HQ)**2
258 IF (IPRO.EQ.91) S0 = MAX(S0,EMMIN**2)
268 S0 = (SQRT(S0)+ABS(PHEP(5,IHADIS)))**2-PHEP(5,IHADIS)**2
270 S0 = MAX(S0,WHMIN**2)
272 ZMIN = S0 / (SS**2 - PHEP(5,IHEP)**2 - PHEP(5,IHADIS)**2)
278 C---UNKNOWN PROCESS: USE ENERGY CUTOFF, AND WARN USER
280 IF (FSTWGT) CALL HWWARN('HWEGAM',1,*999)
282 ZMIN = EGMIN / PHEP(4,IHEP)
296 C---APPLY USER DEFINED CUTS YWWMIN,YWWMAX AND INDIRECT LIMITS ON Z
300 ZMIN=MAX(ZMIN,YWWMIN,SQRT(Q2WWMN)/ABS(PHEP(3,IHEP)))
302 ZMAX=MIN(ZMAX,YWWMAX)
304 IF (ZMIN.GT.ZMAX) THEN
314 C---GENERATE GAMMA MOMENTUM FRACTION
318 10 IF (HWR().LT.A) THEN
320 ZGAM=(ZMIN/ZMAX)**HWR()*ZMAX
324 ZGAM=(ZMAX-ZMIN)*HWR()+ZMIN
328 GAMWT = GAMWT * .5*ALPHEM/PIFAC *
330 + (1+(1-ZGAM)**2)/(A/LOG(ZMAX/ZMIN)+(1-A)/(ZMAX-ZMIN)*ZGAM)
334 GAMWT = GAMWT * LOG((ONE-ZGAM)/ZGAM*(SS/PHEP(5,IHEP))**2)
338 C---Q2WWMN AND Q2WWMX ARE USER-DEFINED LIMITS IN THE Q**2 INTEGRATION
340 QQMAX=MIN(Q2WWMX,(ZGAM*PHEP(3,IHEP))**2)
342 QQMIN=MAX(Q2WWMN,(PHEP(5,IHEP)*ZGAM)**2/(1.-ZGAM))
344 IF (QQMIN.GT.QQMAX) CALL HWWARN('HWEGAM',50,*10)
346 Q2=EXP(HWRUNI(0,LOG(QQMIN),LOG(QQMAX)))
348 GAMWT = GAMWT * LOG(QQMAX/QQMIN)
352 IF (GAMWT.LT.ZERO) GAMWT=ZERO
376 C---FOR COLLINEAR KINEMATICS, ZGAM IS THE ENERGY FRACTION
378 PHEP(4,NHEP)=PHEP(4,IHEP)*ZGAM
380 PHEP(3,NHEP)=PHEP(3,IHEP)-SIGN(SQRT(
382 & (PHEP(4,IHEP)-PHEP(4,NHEP))**2-PHEP(5,IHEP)**2),PHEP(3,IHEP))
388 CALL HWUMAS(PHEP(1,NHEP))
392 C---FOR EXACT KINEMATICS, ZGAM IS TAKEN TO BE FRACTION OF (E+PZ)
394 PPL=ZGAM*(ABS(PHEP(3,IHEP))+PHEP(4,IHEP))
396 QT2=(ONE-ZGAM)*Q2-(ZGAM*PHEP(5,IHEP))**2
400 PHEP(5,NHEP)=-SQRT(Q2)
402 PHEP(4,NHEP)=(PPL+PMI)/TWO
404 PHEP(3,NHEP)=SIGN((PPL-PMI)/TWO,PHEP(3,IHEP))
406 CALL HWRAZM(SQRT(QT2),PHEP(1,NHEP),PHEP(2,NHEP))
410 C---UPDATE OVERALL CM FRAME
414 CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
416 CALL HWVSUM(4,PHEP(1,NHEP),PHEP(1,3),PHEP(1,3))
418 CALL HWUMAS(PHEP(1,3))
420 C---FILL OUTGOING LEPTON
424 IDHW(NHEP)=IDHW(IHEP)
428 IDHEP(NHEP)=IDHEP(IHEP)
440 CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP-1),PHEP(1,NHEP))
442 PHEP(5,NHEP)=PHEP(5,IHEP)
448 *CMZ :- -26/04/91 12.42.30 by Federico Carminati
450 *-- Author : Bryan Webber
452 C-----------------------------------------------------------------------
456 C-----------------------------------------------------------------------
458 C INITIALISES ELEMENTARY PROCESS
460 C-----------------------------------------------------------------------
462 INCLUDE 'HERWIG61.INC'
464 DOUBLE PRECISION HWRSET,DUMMY,SAFETY
468 PARAMETER (SAFETY=1.001)
472 C---NO OF WEIGHT GENERATED
476 C---ACCUMULATED WEIGHTS
480 C---ACCUMULATED WEIGHT-SQUARED
484 C---CURRENT MAX WEIGHT
488 C---LAST VALUE OF SCALE
492 C---NUMBER OF ERRORS REPORTED
496 C---NUMBER OF ERRORS UNREPORTED
500 C---FIND MAXIMUM EVENT WEIGHT IN CASES WHERE THIS IS REQUIRED
504 IF (WGTMAX.EQ.ZERO) THEN
510 WRITE(6,10) IPROC,IBRN,NBSH
512 10 FORMAT(/10X,'INITIAL SEARCH FOR MAX WEIGHT'//
514 & 10X,'PROCESS CODE IPROC = ',I11/
516 & 10X,'RANDOM NO. SEED 1 = ',I11/
518 & 10X,' SEED 2 = ',I11/
520 & 10X,'NUMBER OF SHOTS = ',I11)
532 20 FORMAT(/10X,'INITIAL SEARCH FINISHED')
534 IF (WBIGST*NWGTS.LT.SAFETY*WGTSUM)
536 & WGTMAX=SAFETY*WBIGST
550 WRITE(6,21) AVWGT,WGTMAX
552 21 FORMAT(/1P,10X,'INPUT EVT WEIGHT =',E12.4/
554 & 10X,'INPUT MAX WEIGHT =',E12.4)
560 C---RESET RANDOM NUMBER
570 *CMZ :- -01/04/99 19.55.17 by Mike Seymour
572 *-- Author : Mike Seymour
574 C-----------------------------------------------------------------------
576 SUBROUTINE HWEISR(IHEP)
578 C-----------------------------------------------------------------------
580 C GENERATES AN ISR PHOTON FROM INCOMING E+, E-, MU+ OR MU-
582 C-----------------------------------------------------------------------
584 INCLUDE 'HERWIG61.INC'
586 DOUBLE PRECISION HWR,QSQMAX,QSQMIN,A,B,B1,B2,B3,B4,B5,B6,B7,B8,
588 $ R,AA,T0,T1,C1,C2,T,Z(2),QSQ(2),PHI(2),C
596 C---IF ZMXISR IS ZERO, THERE CAN BE NO ISR
598 IF (ZMXISR.EQ.ZERO .OR. (IPRO.GT.3.AND.IPRO.NE.6)) RETURN
600 C---CHECK CONSISTENCY OF TMNISR AND ZMXISR
602 IF (ZMXISR**2.LT.TMNISR) CALL HWWARN('HWEISR',200,*999)
604 C---CALCULATE VIRTUALITY LIMITS
606 QSQMAX=4*PHEP(4,IHEP)**2
608 QSQMIN=PHEP(5,IHEP)**2
610 C---AND THEREFORE THE Z DEPENDENCE
614 B=A*(LOG(QSQMAX/QSQMIN)-1)
616 C---DECIDE HOW MUCH WEIGHT TO GIVE THE Z RESONANCE
620 IF (IPRO.EQ.1.OR.IPRO.EQ.6) THEN
624 ELSEIF (IPRO.EQ.2) THEN
628 ELSEIF (IPRO.EQ.3) THEN
638 T0=RMASS(200)**2/QSQMAX
640 T1=GAMZ*RMASS(200)/QSQMAX
652 C---GENERATE A T VALUE BETWEEN TMNISR AND 1 ACCORDING TO:
654 C ( b**2*log(zmxisr**2/t)/t + 2*b*(1-(1-zmxisr)**b)*((1-t)**(2*b-1)+1/t
656 C +(1-t0)**(2b-1)*aa*t1/((t-t0)**2+t1**2)) ) *theta(zmxisr**2-t)
658 C +( 2*b*(1-zmxisr)**b*((1-t)**(b-1)+1/t
660 C +(1-t0)**(b-1)*aa*t1/((t-t0)**2+t1**2)) ) *theta(zmxisr-t)
662 C +( (1-zmxisr)**(2*b) ) *delta(1-t)
666 B2=B1+2*(1-ZMXISR)**B*((1-TMNISR)**B-(1-ZMXISR)**B)
668 B3=B2+2*B*(1-ZMXISR)**B*LOG(ZMXISR/TMNISR)
670 B4=B3+2*B*(1-ZMXISR)**B*AA*(1-T0)**(B-1)
672 $ *(ATAN((ZMXISR-T0)/T1)-ATAN((TMNISR-T0)/T1))
674 B5=B4+(1-(1-ZMXISR)**B)*((1-TMNISR)**(2*B)-(1-ZMXISR**2)**(2*B))
676 B6=B5+2*B*(1-(1-ZMXISR)**B)*LOG(ZMXISR**2/TMNISR)
678 B7=B6+B**2*LOG(ZMXISR**2/TMNISR)**2/2
680 B8=B7+2*B*(1-(1-ZMXISR)**B)*AA*(1-T0)**(2*B-1)
682 $ *(ATAN((ZMXISR**2-T0)/T1)-ATAN((TMNISR-T0)/T1))
696 ELSEIF (R.LE.B4) THEN
704 T=1-(1-TMNISR)*(1-R*(1-((1-ZMXISR)/(1-TMNISR))**B))**(1/B)
706 ELSEIF (R.LE.B3) THEN
710 T=(TMNISR/ZMXISR)**R*ZMXISR
718 $ ATAN((ZMXISR-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
722 GAMWT=GAMWT*B8/(2*B*(1-ZMXISR)**B*((1-T)**(B-1)+1/T+
724 $ (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
728 IF (HWR().GT.HALF) Z(1)=T
742 $ (1-R*(1-((1-ZMXISR**2)/(1-TMNISR))**(2*B)))**(.5/B)
744 ELSEIF (R.LE.B6) THEN
748 T=(TMNISR/ZMXISR**2)**R*ZMXISR**2
750 ELSEIF (R.LE.B7) THEN
754 T=(TMNISR/ZMXISR**2)**SQRT(R)*ZMXISR**2
762 $ ATAN((ZMXISR**2-T0)/T1)*R+ATAN((TMNISR-T0)/T1)*(1-R))
766 GAMWT=GAMWT*B8/(B**2*LOG(ZMXISR**2/T)/T
768 $ + 2*B*(1-(1-ZMXISR)**B)*((1-T)**(2*B-1)+1/T+
770 $ (1-T0)**(B-1)*AA*T1/((T-T0)**2+T1**2)))
772 C---GENERATE A Z VALUE BETWEEN T/ZMXISR AND ZMXISR ACCORDING TO:
774 C 1/z+(1-z)**(b-1)+t/z**2*(1-t/z)**(b-1)
778 C2=C1+2/B*((1-T/ZMXISR)**B-(1-ZMXISR)**B)
786 Z(1)=(T/ZMXISR**2)**HWR()*ZMXISR
792 $ (1-HWR()*(1-((1-ZMXISR)/(1-T/ZMXISR))**B))**(1/B)
794 IF (2*R.LE.C2+C1) Z(1)=T/Z(1)
806 $ /(1/Z(1)+(1-Z(1))**(B-1)+T/Z(1)**2*(1-T/Z(1))**(B-1))
810 C---INCLUDE DISTRIBUTION FUNCTIONS
816 IF (Z(I).GT.ZMXISR) THEN
820 GAMWT=GAMWT*(1-ZMXISR)**B*EXP(3*B/4)*(1-B**2*PIFAC**2/12)
824 GAMWT=GAMWT*(B*(1-Z(I))**(B-1)*(1+Z(I)**2)/2
826 $ *EXP(B*Z(I)/2*(1+Z(I)/2))*(1-B**2*PIFAC**2/12)
828 $ +B**2/8*((1+Z(I))*((1+Z(I))**2+3*LOG(Z(I)))
830 $ -4*LOG(Z(I))/(1-Z(I))))
836 C---CHOOSE BOTH QSQ VALUES
840 IF (Z(I).GT.ZMXISR .OR. COLISR) THEN
848 C---ACCORDING TO 1/(QSQ+QSQMIN) FROM 0 TO (1-Z)*(T/(Z+T))*QSQMAX
850 20 QSQ(I)=(((1-Z(I))*(T/(Z(I)+T))
852 $ *QSQMAX/QSQMIN+1)**HWR()-1)*QSQMIN
854 C---AND REJECT TO QSQ/(QSQ+QSQMIN)**2
856 IF (HWR()*(QSQ(I)+QSQMIN).GT.QSQ(I)) GOTO 20
862 C---CHOOSE BOTH AZIMUTHS
868 C---USE S-HAT PRESCRIPTION TO MODIFY Z VALUES
872 IF ((1-Z(1))*QSQ(1).GT.(1-Z(2))*QSQ(2)) I=1
874 IF ((1-Z(2))*QSQ(2).GT.(1-Z(1))*QSQ(1)) I=2
880 Z(I)=Z(I)+QSQ(I)/QSQMAX
882 IF (QSQ(J).GT.ZERO) THEN
884 Z(J)=((QSQ(I)*QSQMAX+QSQ(J)*QSQMAX
886 $ -QSQ(I)*QSQ(J))/QSQMAX**2+T)/Z(I)
888 C=COS(PHI(1)-PHI(2))*SQRT(QSQ(1)*QSQ(2))/QSQMAX
890 Z(J)=Z(J)+(-2*C**2*(1-Z(I))+2*C*SQRT((1-Z(I))
892 $ *(C**2*(1-Z(I))+Z(I)**2*(1-Z(J)))))/Z(I)**2
898 ELSEIF (IHEP.EQ.2) THEN
900 C---EVERYTHING WAS GENERATED LAST TIME
904 C---ROUTINE CALLED UNEXPECTEDLY
906 CALL HWWARN('HWEISR',201,*999)
910 C---IF Z IS TOO LARGE THERE IS NO EMISSION
912 IF (Z(IHEP).GT.ZMXISR) RETURN
914 C---PUT NEW LEPTON IN EVENT RECORD
918 IDHW(NHEP)=IDHW(IHEP)
920 IDHEP(NHEP)=IDHEP(IHEP)
934 C---AND OUTGOING PHOTON
954 C---RECONSTRUCT PHOTON KINEMATICS (Z IS LIGHT-CONE MOMENTUM FRACTION)
956 PHEP(1,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*COS(PHI(IHEP))
958 PHEP(2,NHEP)=SQRT(QSQ(IHEP)*(1-Z(IHEP)))*SIN(PHI(IHEP))
960 PHEP(3,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)-QSQ(IHEP)/(4*PHEP(4,IHEP))
962 IF (IHEP.EQ.2) PHEP(3,NHEP)=-PHEP(3,NHEP)
964 PHEP(4,NHEP)=(1-Z(IHEP))*PHEP(4,IHEP)+QSQ(IHEP)/(4*PHEP(4,IHEP))
970 CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,NHEP),PHEP(1,NHEP-1))
972 CALL HWUMAS(PHEP(1,NHEP-1))
974 C---UPDATE OVERALL CM FRAME
976 JMOHEP(IHEP,3)=NHEP-1
978 CALL HWVDIF(4,PHEP(1,3),PHEP(1,IHEP),PHEP(1,3))
980 CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,3),PHEP(1,3))
982 CALL HWUMAS(PHEP(1,3))
988 *CMZ :- -26/04/91 11.11.55 by Bryan Webber
990 *-- Author : Bryan Webber
992 C-----------------------------------------------------------------------
996 C-----------------------------------------------------------------------
998 C SETS UP 2->1 (COLOUR SINGLET) HARD SUBPROCESS
1000 C-----------------------------------------------------------------------
1002 INCLUDE 'HERWIG61.INC'
1006 INTEGER ICMF,I,IBM,IHEP
1016 C---FIND BEAM AND TARGET
1018 IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1024 IDHEP(IHEP)=IDPDG(IDN(I))
1034 C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1036 IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1038 CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1040 IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1048 PHEP(5,IHEP)=RMASS(IDN(I))
1050 PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1052 PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1054 PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1060 PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1062 C---HARD CENTRE OF MASS
1066 IDHEP(ICMF)=IDPDG(IDCMF)
1070 CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1072 CALL HWUMAS(PHEP(1,ICMF))
1074 C---SET UP COLOUR STRUCTURE LABELS
1076 JMOHEP(2,NHEP+1)=NHEP+2
1078 JDAHEP(2,NHEP+1)=NHEP+2
1080 JMOHEP(2,NHEP+2)=NHEP+1
1082 JDAHEP(2,NHEP+2)=NHEP+1
1084 JDAHEP(1,NHEP+3)=NHEP+3
1086 JDAHEP(2,NHEP+3)=NHEP+3
1094 *CMZ :- -01/04/99 19.41.18 by Mike Seymour
1096 *-- Author : Bryan Webber
1098 C-----------------------------------------------------------------------
1102 C-----------------------------------------------------------------------
1104 C WHEN NEVHEP=0, CHOOSES X VALUES AND FINDS WEIGHT FOR PROCESS IPROC
1106 C OTHERWISE, CHOOSES AND LOADS ALL VARIABLES FOR HARD PROCESS
1108 C-----------------------------------------------------------------------
1110 INCLUDE 'HERWIG61.INC'
1112 DOUBLE PRECISION HWR
1116 IF (IERROR.NE.0) RETURN
1118 C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS NOT ACCEPTED
1122 C---FSTWGT IS .TRUE. DURING FIRST CALL TO HARD PROCESS ROUTINE
1126 C---FSTEVT IS .TRUE. THROUGHOUT THE FIRST EVENT
1130 C---SET COLOUR CORRECTION TO FALSE
1138 C---SET UP INITIAL STATE
1152 PHEP(5,NHEP)=RMASS(IPART1)
1164 IDHEP(NHEP)=IDPDG(IPART1)
1174 PHEP(3,NHEP)=-PBEAM2
1178 PHEP(5,NHEP)=RMASS(IPART2)
1190 IDHEP(NHEP)=IDPDG(IPART2)
1192 C---NEXT ENTRY IS OVERALL CM FRAME
1202 JMOHEP(1,NHEP)=NHEP-2
1204 JMOHEP(2,NHEP)=NHEP-1
1210 CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP-2),PHEP(1,NHEP))
1212 CALL HWUMAS(PHEP(1,NHEP))
1214 C Select a primary interaction point
1222 CALL HWVZRO(4,VTXPIP)
1226 CALL HWVEQU(3,VTXPIP,VHEP(1,NHEP))
1230 C---GENERATE PHOTONS (WEIZSACKER-WILLIAMS APPROX)
1232 C FOR HADRONIC PROCESSES WITH LEPTON BEAMS
1236 IF (IPRO.GT.10.AND.IPRO.LT.90) THEN
1238 IF (ABS(IDHEP(1)).EQ.11.OR.ABS(IDHEP(1)).EQ.13)
1240 & CALL HWEGAM(1,ZERO, ONE,.FALSE.)
1242 IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1244 & CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1246 ELSEIF (IPRO.GE.90) THEN
1248 IF (ABS(IDHEP(2)).EQ.11.OR.ABS(IDHEP(2)).EQ.13)
1250 & CALL HWEGAM(2,ZERO, ONE,.FALSE.)
1254 C---GENERATE ISR PHOTONS FOR LEPTONIC PROCESSES
1256 IF (IPRO.LT.10) THEN
1264 C---IF USER LIMITS WERE TOO TIGHT, MIGHT NOT BE ANY PHASE-SPACE
1266 IF (GAMWT.LE.ZERO) GOTO 30
1268 C---IF CMF HAS ACQUIRED A TRANSVERSE BOOST, OR USER REQUESTS IT ANYWAY,
1270 C BOOST EVENT RECORD BACK TO CMF
1272 IF (PHEP(1,3)**2+PHEP(2,3)**2.GT.ZERO .OR. USECMF) CALL HWUBST(1)
1274 C---ROUTINE LOOPS BACK TO HERE IF GENERATED WEIGHT WAS ACCEPTED
1278 C---IPRO=MOD(IPROC/100,100)
1282 IF (IPROC.LT.110.OR.IPROC.GE.120) THEN
1284 C--- E+E- -> Q-QBAR OR L-LBAR
1290 C--- E+E- -> Q-QBAR-GLUON
1296 ELSEIF (IPRO.EQ.2) THEN
1302 ELSEIF (IPRO.EQ.3) THEN
1308 ELSEIF (IPRO.EQ.4) THEN
1310 C---E+E- -> NUEB NUE H
1314 ELSEIF (IPRO.EQ.5 .AND. IPROC.LT.550) THEN
1316 C---EE -> EE GAMGAM -> EE FFBAR/WW
1320 ELSEIF (IPRO.EQ.5) THEN
1322 C---EE -> ENU GAMW -> ENU FF'BAR/WZ
1326 ELSEIF (IPRO.EQ.6) THEN
1332 ELSEIF (IPRO.EQ.13) THEN
1334 C---GAMMA/Z0/Z' DRELL-YAN PROCESS
1338 ELSEIF (IPRO.EQ.14) THEN
1340 C---W+/- PRODUCTION VIA DRELL-YAN PROCESS
1344 ELSEIF (IPRO.EQ.15) THEN
1346 C---QCD HARD 2->2 PROCESSES
1350 ELSEIF (IPRO.EQ.16) THEN
1352 C---HIGGS PRODUCTION VIA GLUON FUSION
1356 ELSEIF (IPRO.EQ.17) THEN
1358 C---QCD HEAVY FLAVOUR PRODUCTION
1362 ELSEIF (IPRO.EQ.18) THEN
1364 C---QCD DIRECT PHOTON + JET PRODUCTION
1368 ELSEIF (IPRO.EQ.19) THEN
1370 C---HIGGS PRODUCTION VIA W FUSION
1374 ELSEIF (IPRO.EQ.20) THEN
1376 C---TOP PRODUCTION FROM W EXCHANGE
1380 ELSEIF (IPRO.EQ.21) THEN
1382 C---VECTOR BOSON + JET PRODUCTION
1386 ELSEIF (IPRO.EQ.22) THEN
1388 C QCD direct photon pair production
1392 ELSEIF (IPRO.EQ.23) THEN
1394 C QCD Higgs plus jet production
1398 ELSEIF (IPRO.EQ.24) THEN
1400 C---COLOUR-SINGLET EXCHANGE
1404 ELSEIF (IPRO.EQ.30) THEN
1406 C---HADRON-HADRON SUSY PROCESSES
1410 ELSEIF(IPRO.EQ.40.OR.IPRO.EQ.41) THEN
1412 C---HADRON-HADRON R-PARITY VIOLATING SUSY PROCESSES
1416 ELSEIF (IPRO.EQ.50) THEN
1418 C Point-like photon two-jet production
1422 ELSEIF (IPRO.EQ.51) THEN
1424 C Point-like photon/QCD heavy flavour pair production
1428 ELSEIF (IPRO.EQ.52) THEN
1430 C Point-like photon/QCD heavy flavour single excitation
1434 ELSEIF (IPRO.EQ.53) THEN
1436 C Compton scattering of point-like photon and (anti)quark
1440 ELSEIF (IPRO.EQ.55) THEN
1442 C Point-like photon/higher twist meson production
1446 ELSEIF (IPRO.EQ.60) THEN
1448 C---QPM GAMMA-GAMMA-->QQBAR
1452 ELSEIF (IPRO.GE.70.AND.IPRO.LE.79) THEN
1454 C---BARYON-NUMBER VIOLATION, AND OTHER MULTI-W PRODUCTION PROCESSES
1458 ELSEIF (IPRO.EQ.80) THEN
1460 C---MINIMUM-BIAS: NO HARD SUBPROCESS
1466 ELSEIF (IPRO.EQ.90) THEN
1472 ELSEIF(IPRO.EQ.91) THEN
1474 C---BOSON - GLUON(QUARK) FUSION --> ANTIQUARK(GLUON) + QUARK
1478 ELSEIF(IPRO.EQ.92) THEN
1480 C---DEEP INELASTIC WITH EXTRA JET: OBSOLETE PROCESS
1484 40 FORMAT (1X,' IPROC=92** is no longer supported.'
1486 & /1X,' Please use IPROC=91** instead.')
1488 CALL HWWARN('HWEPRO',500,*999)
1490 ELSEIF(IPRO.EQ.95) THEN
1492 C---HIGGS PRODUCTION VIA W FUSION IN E P
1500 CALL HWWARN('HWEPRO',102,*999)
1506 IF (NOWGT) EVWGT=AVWGT
1514 C---IF AN EVENT IS CANCELLED BEFORE IT IS GENERATED, GIVE IT ZERO WEIGHT
1516 IF (IERROR.NE.0) THEN
1530 WSQSUM=WSQSUM+EVWGT**2
1532 IF (EVWGT.GT.WBIGST) THEN
1536 IF (NOWGT.AND.WBIGST.GT.WGTMAX) THEN
1538 IF (NEVHEP.NE.0) CALL HWWARN('HWEPRO',1,*999)
1546 ELSEIF (EVWGT.LT.ZERO) THEN
1548 IF (EVWGT.LT.-1.D-9) CALL HWWARN('HWEPRO',3,*999)
1554 IF (NEVHEP.NE.0) THEN
1556 C---LOW EFFICIENCY WARNINGS:
1558 C WARN AT 10*EFFMIN, STOP AT EFFMIN
1560 IF (10*EFFMIN*NWGTS.GT.NEVHEP) THEN
1562 IF (EFFMIN*NWGTS.GT.NEVHEP) CALL HWWARN('HWEPRO',200,*999)
1564 IF (EFFMIN.GT.ZERO) THEN
1566 IF (MOD(NWGTS,INT(10/EFFMIN)).EQ.0) THEN
1568 CALL HWWARN('HWEPRO',2,*999)
1580 GENEV=EVWGT.GT.WGTMAX*HWR()
1596 98 FORMAT(10X,' MAXIMUM WEIGHT =',1PG24.16)
1598 99 FORMAT(10X,'NEW MAXIMUM WEIGHT =',1PG24.16)
1604 *CMZ :- -26/04/91 11.11.55 by Bryan Webber
1606 *-- Author : Bryan Webber
1608 C-----------------------------------------------------------------------
1612 C-----------------------------------------------------------------------
1614 C SETS UP 2->2 HARD SUBPROCESS
1616 C-----------------------------------------------------------------------
1618 INCLUDE 'HERWIG61.INC'
1620 DOUBLE PRECISION HWUMBW,HWUPCM,PA,PCM
1622 INTEGER ICMF,IBM,I,J,K,IHEP,NTRY
1634 C---FIND BEAM AND TARGET
1636 IF (JDAHEP(1,I).NE.0) IBM=JDAHEP(1,I)
1642 IDHEP(IHEP)=IDPDG(IDN(I))
1652 C---SPECIAL - IF INCOMING PARTON IS INCOMING BEAM THEN COPY IT
1654 IF (XX(I).EQ.ONE.AND.IDHW(IBM).EQ.IDN(I)) THEN
1656 CALL HWVEQU(5,PHEP(1,IBM),PHEP(1,IHEP))
1658 IF (I.EQ.2) PHEP(3,IHEP)=-PHEP(3,IHEP)
1666 PHEP(5,IHEP)=RMASS(IDN(I))
1668 PA=XX(I)*(PHEP(4,IBM)+ABS(PHEP(3,IBM)))
1670 PHEP(4,IHEP)=0.5*(PA+PHEP(5,IHEP)**2/PA)
1672 PHEP(3,IHEP)=PA-PHEP(4,IHEP)
1678 PHEP(3,NHEP+2)=-PHEP(3,NHEP+2)
1680 C---HARD CENTRE OF MASS
1684 IDHEP(ICMF)=IDPDG(IDCMF)
1688 CALL HWVSUM(4,PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,ICMF))
1690 CALL HWUMAS(PHEP(1,ICMF))
1704 IDHEP(IHEP)=IDPDG(IDN(I))
1710 JDAHEP(I-2,ICMF)=IHEP
1712 20 PHEP(5,IHEP)=HWUMBW(IDN(I))
1714 PCM=HWUPCM(PHEP(5,NHEP+3),PHEP(5,NHEP+4),PHEP(5,NHEP+5))
1716 IF (PCM.LT.ZERO) THEN
1720 IF (NTRY.LE.NETRY) GO TO 19
1722 CALL HWWARN('HWETWO',103,*999)
1728 PHEP(4,IHEP)=SQRT(PCM**2+PHEP(5,IHEP)**2)
1730 PHEP(3,IHEP)=PCM*COSTH
1732 PHEP(1,IHEP)=SQRT((PCM+PHEP(3,IHEP))*(PCM-PHEP(3,IHEP)))
1734 CALL HWRAZM(PHEP(1,IHEP),PHEP(1,IHEP),PHEP(2,IHEP))
1736 CALL HWULOB(PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,IHEP))
1738 CALL HWVDIF(4,PHEP(1,NHEP+3),PHEP(1,IHEP),PHEP(1,NHEP+5))
1740 C---SET UP COLOUR STRUCTURE LABELS
1752 JMOHEP(2,NHEP+J)=NHEP+K
1754 30 JDAHEP(2,NHEP+K)=NHEP+J
1762 *CMZ :- -01/04/99 19.47.55 by Mike Seymour
1764 *-- Author : Ian Knowles
1766 C-----------------------------------------------------------------------
1770 C-----------------------------------------------------------------------
1772 C Four jet production in e^+e^- annihilation: qqbar+gg & qqbar+qqbar
1774 C IOP4JT controls the treatment of the colour flow interference term
1778 C IOP4JT(1)=0 neglect, =1 extreme 2341; =2 extreme 3421
1780 C qqbar-qqbar (identical quark flavour) case:
1782 C IOP4JT(2)=0 neglect, =1 extreme 4123; =2 extreme 2143
1786 C Matrix elements based on Ellis Ross & Terrano and Catani & Seymour
1790 C WARNING: Phase space factor inaccurate for JADE y_cut > 0.14.
1792 C-----------------------------------------------------------------------
1794 INCLUDE 'HERWIG61.INC'
1796 INTEGER LM,LP,IQK,I,J,IDMN,IDMX,ID1,ID2,IST(4)
1798 DOUBLE PRECISION HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,
1800 & HWH4J4,HWH4J5,HWH4J6,HWH4J7,QNOW,Q2NOW,QLST,SCUT,PSFAC,FACT,X,
1802 & COLA,COLB,COLC,CLF(7,6),P12,P13,P14,P23,P24,P34,FACTR,EP1,EP2,
1804 & EP3,EP4,GG1,GG2,GG12,GG3,GG13,GG23,GGINT,WTGG,QQ,QP,QQINT,QQ1,
1806 & QQ2,WTQQ,WTQP,HCS,WTAB,WTBA,WTOT,RCS
1810 LOGICAL INCLQG(6),INCLQQ(6,6),ORIENT
1812 EXTERNAL HWR,HWUALF,HWUAEM,HWULDO,HWH4J1,HWH4J2,HWH4J4,
1814 & HWH4J5,HWH4J6,HWH4J7
1816 SAVE HCS,QLST,WTQP,WTQQ,WTGG,FACTR,COLA,COLB,COLC,IDMN,IDMX,
1818 & CLF,GG1,GG2,GGINT,INCLQG,INCLQQ,LM,LP,QQ1,QQ2,QQINT,FACT,ORIENT,
1822 DATA QLST,IST/-1.,113,114,114,114/
1832 IF (NHEP+5.GT.NMXHEP) CALL HWWARN('HWH4JT',100,*999)
1836 IF (QNOW.NE.QLST) THEN
1844 C Calculate allowed fraction of Phase Space using parameterization
1848 PSFAC=(1.-6.*Y4JT)**5.50*(1.-173.3*Y4JT*(1.-247.3*Y4JT
1850 & *(1.+148.3*Y4JT*(1.+3.913*Y4JT))))
1852 & /(1.-8.352*Y4JT*(1.-1102.*Y4JT
1854 & *(1.+1603.*Y4JT*(1.+22.99*Y4JT))))
1858 PSFAC=(1.-6.*Y4JT)**4.62*(1.-44.72*Y4JT*(1.-176.0*Y4JT
1860 & *(1.+102.9*Y4JT*(1.-6.579*Y4JT))))
1862 & /(1.-3.392*Y4JT*(1.-946.5*Y4JT
1864 & *(1.+423.4*Y4JT*(1.-3.971*Y4JT))))
1868 FACT=GEV2NB*HWUAEM(Q2NOW)**2*CFFAC*FLOAT(NCOLO)*PSFAC
1874 COLB=CFFAC-HALF*CAFAC
1880 IF (JDAHEP(1,LM).NE.0) LM=JDAHEP(1,LM)
1884 IF (JDAHEP(1,LP).NE.0) LP=JDAHEP(1,LP)
1904 CALL HWUCFF(11,I,Q2NOW,CLF(1,I))
1906 IF (QNOW.GT.TWO*(RMASS(I)+RMASS(13))) THEN
1918 IF (QNOW.GT.TWO*(RMASS(I)+RMASS(J ))) THEN
1934 IF (MOD(IPROC/10,10).EQ.5) THEN
1946 C Generate phase space point and check it passes cuts
1948 CALL HWVEQU(5,PHEP(1,3),PHEP(1,NHEP+1))
1952 20 PHEP(5,NHEP+I)=0.
1954 30 CALL HWDFOR(PHEP(1,NHEP+1),PHEP(1,NHEP+2),PHEP(1,NHEP+3),
1956 & PHEP(1,NHEP+4),PHEP(1,NHEP+5))
1960 P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
1962 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+3),
1964 & PHEP(4,NHEP+3)/PHEP(4,NHEP+2))*P12
1968 P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
1970 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+4),
1972 & PHEP(4,NHEP+4)/PHEP(4,NHEP+2))*P13
1976 P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
1978 X=MIN(PHEP(4,NHEP+2)/PHEP(4,NHEP+5),
1980 & PHEP(4,NHEP+5)/PHEP(4,NHEP+2))*P14
1984 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
1986 X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+4),
1988 & PHEP(4,NHEP+4)/PHEP(4,NHEP+3))*P23
1992 P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
1994 X=MIN(PHEP(4,NHEP+3)/PHEP(4,NHEP+5),
1996 & PHEP(4,NHEP+5)/PHEP(4,NHEP+3))*P24
2000 P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2002 X=MIN(PHEP(4,NHEP+4)/PHEP(4,NHEP+5),
2004 & PHEP(4,NHEP+5)/PHEP(4,NHEP+4))*P34
2006 IF (X.GT.SCUT) GOTO 40
2020 P12=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+3))
2022 IF (P12.GT.SCUT) THEN
2024 P13=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+4))
2026 IF (P13.GT.SCUT) THEN
2028 P14=2*HWULDO(PHEP(1,NHEP+2),PHEP(1,NHEP+5))
2030 IF (P14.GT.SCUT) THEN
2032 P23=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+4))
2034 IF (P23.GT.SCUT) THEN
2036 P24=2*HWULDO(PHEP(1,NHEP+3),PHEP(1,NHEP+5))
2038 IF (P24.GT.SCUT) THEN
2040 P34=2*HWULDO(PHEP(1,NHEP+4),PHEP(1,NHEP+5))
2042 IF (P34.GT.SCUT) GOTO 40
2060 C Passed cuts: calculate contributions to Matrix Elements
2062 40 EMSCA=SQRT(MIN(P12,P13,P14,P23,P24,P34))
2064 FACTR=FACT*HWUALF(1,EMSCA)**2
2068 QF=HWULDO(PHEP(1,LP),PHEP(1,3))
2070 EF=Q2NOW/(2*SQRT(QF**2-HWULDO(PHEP(1,LP),PHEP(1,LP))*Q2NOW))
2076 E(I)=EF*PHEP(I,LP)+QF*PHEP(I,3)
2080 EP1=HWULDO(E,PHEP(1,NHEP+2))
2082 EP2=HWULDO(E,PHEP(1,NHEP+3))
2084 EP3=HWULDO(E,PHEP(1,NHEP+4))
2086 EP4=HWULDO(E,PHEP(1,NHEP+5))
2092 GG1=HWH4J1(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2094 & +HWH4J1(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2096 GG2=HWH4J1(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2098 & +HWH4J1(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2100 GG12=HWH4J2(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2102 & +HWH4J2(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2104 & +HWH4J2(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2106 & +HWH4J2(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2108 GG3=HWH4J4(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2110 & +HWH4J4(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2112 GG13=GG3+HWH4J5(P12,P13,P14,P23,P24,P34,EP1,EP2,EP3,EP4,ORIENT)
2114 & +HWH4J5(P12,P24,P23,P14,P13,P34,EP2,EP1,EP4,EP3,ORIENT)
2116 GG23=GG3+HWH4J5(P12,P14,P13,P24,P23,P34,EP1,EP2,EP4,EP3,ORIENT)
2118 & +HWH4J5(P12,P23,P24,P13,P14,P34,EP2,EP1,EP3,EP4,ORIENT)
2122 GG1 =COLA*(GG1 +GG13)
2124 GG2 =COLA*(GG2 +GG23)
2126 GGINT=COLB*(GG12-GG13-GG23)
2128 WTGG=FACTR*(GG1+GG2+GGINT)
2132 QP=HWH4J6(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2134 & +HWH4J6(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2136 & +HWH4J6(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2138 & +HWH4J6(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2140 QQ=HWH4J6(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2142 & +HWH4J6(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2144 & +HWH4J6(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2146 & +HWH4J6(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2148 QQINT=HWH4J7(P13,P12,P14,P23,P34,P24,EP1,EP3,EP2,EP4,ORIENT)
2150 & +HWH4J7(P24,P12,P23,P14,P34,P13,EP2,EP4,EP1,EP3,ORIENT)
2152 & +HWH4J7(P13,P23,P34,P12,P14,P24,EP3,EP1,EP2,EP4,ORIENT)
2154 & +HWH4J7(P24,P23,P12,P34,P14,P13,EP2,EP4,EP3,EP1,ORIENT)
2156 & +HWH4J7(P13,P14,P12,P34,P23,P24,EP1,EP3,EP4,EP2,ORIENT)
2158 & +HWH4J7(P24,P14,P34,P12,P23,P13,EP4,EP2,EP1,EP3,ORIENT)
2160 & +HWH4J7(P13,P34,P23,P14,P12,P24,EP3,EP1,EP4,EP2,ORIENT)
2162 & +HWH4J7(P24,P34,P14,P23,P12,P13,EP4,EP2,EP3,EP1,ORIENT)
2166 WTQP=FACTR*COLC*QP/TWO
2174 WTQQ=FACTR*(QQ1+QQ2+QQINT)/2
2184 IF (INCLQG(ID1)) THEN
2188 HCS=HCS+CLF(1,ID1)*WTGG
2190 IF (GENEV.AND.HCS.GT.RCS) THEN
2192 C Select colour flow
2198 IF (IOP4JT(1).EQ.1) THEN
2200 IF (GGINT.GE.ZERO) THEN
2206 WTBA=MAX(WTBA,WTBA+GGINT)
2210 ELSEIF (IOP4JT(1).EQ.2) THEN
2212 IF (GGINT.GE.ZERO) THEN
2218 WTAB=MAX(WTAB,WTAB+GGINT)
2222 ELSEIF (IOP4JT(1).NE.0) THEN
2224 CALL HWWARN('HWH4JT',101,*999)
2230 IF (WTAB.GT.HWR()*WTOT) THEN
2232 CALL HWHQCP( 13, 13,3142,91,*99)
2236 CALL HWHQCP( 13, 13,4123,92,*99)
2248 C Identical quark pairs
2250 IF (ID1.EQ.ID2.AND.INCLQQ(ID1,ID1)) THEN
2252 HCS=HCS+CLF(1,ID1)*WTQQ
2254 IF (GENEV.AND.HCS.GT.RCS) THEN
2256 C Select colour flow
2262 IF (IOP4JT(2).EQ.1) THEN
2264 IF (QQINT.GE.ZERO) THEN
2270 WTBA=MAX(WTBA,WTBA+QQINT)
2274 ELSEIF (IOP4JT(2).EQ.2) THEN
2276 IF (QQINT.GE.ZERO) THEN
2282 WTAB=MAX(WTAB,WTAB+QQINT)
2286 ELSEIF (IOP4JT(2).NE.0) THEN
2288 CALL HWWARN('HWH4JT',102,*999)
2294 IF (WTAB.GT.HWR()*WTOT) THEN
2296 CALL HWHQCP(ID1,ID1+6,4123,93,*99)
2300 CALL HWHQCP(ID1,ID1+6,2143,94,*99)
2306 C Unlike quark pairs
2308 ELSEIF (INCLQQ(ID1,ID2)) THEN
2310 HCS=HCS+(CLF(1,ID1)+CLF(1,ID2))*WTQP
2312 IF (GENEV.AND.HCS.GT.RCS) CALL HWHQCP(ID2,ID2+6,4123,95,*99)
2324 C Set up labels for selected final state
2352 IDHEP(J)=IDPDG(IDN(I))
2360 C And colour structure pointers
2366 JMOHEP(2,NHEP+1+I)=NHEP+1+J
2368 110 JDAHEP(2,NHEP+1+J)=NHEP+1+I
2376 *CMZ :- -01/04/99 19.47.55 by Mike Seymour
2378 *-- Author : Ian Knowles
2380 C-----------------------------------------------------------------------
2382 FUNCTION HWH4J1(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2384 C-----------------------------------------------------------------------
2386 C Evaluate `ERT' functions A, B, C, D, E; S12=(p1+p2)^2 etc.
2388 C-----------------------------------------------------------------------
2392 DOUBLE PRECISION HWH4J1,HWH4J2,HWH4J4,HWH4J5,HWH4J6,HWH4J7,
2394 & S12,S13,S14,S23,S24,S34,S123,S124,S134,S234,S,EP1,EP2,EP3,EP4,
2408 S=S12+S13+S14+S23+S24+S34
2410 HWH4J1=(S12*((S12+S14+S23+S34)**2+S13*(S12+S14-S24)+S24*(S12+S23))
2412 & +(S14*S23-S12*S34-S13*S24)*(S14+S23+S34)/2)
2414 & /(S13*S24*S134*S234)
2416 & +((S12+S24)*(S13+S34)-S14*S23)/(S13*S134**2)
2418 & +2*S23*(S-S13)/(S13*S134*S24) + S34/(2*S13*S24)
2424 & +4*((EP1*EP1*((S-S13)*(S23+S24)-S24*S34)
2426 & -EP1*EP2*(S12*(S123+S124)+(S+S12)*(S14+S23)+2*S14*S23
2428 & +S24*S134+S234*(S13+2*S234))
2430 & +EP1*EP3*(S*(S24-S12)+S12*S13+(S14+2*S234-S34)*S24)
2432 & -EP1*EP4*(S12*S124+S23*(S+S12+S14))
2434 & +EP2*EP2*((S-S24)*(S13+S14)+2*(S13+S34)*S234-S13*S34)
2436 & -EP2*EP3*((S+S23)*(S12+S14)+(S12+2*(S23+S234))*S234)
2438 & +EP2*EP4*(S12*(S24-S)+S13*(S+S23-S34)+2*(S13+S34-S234)*S234)
2440 & +EP3*EP3*(S14+2*S234)*S24
2442 & +EP3*EP4*(-S234*(2*(S12+S23)+S134)+S12*S34-S13*S24-S14*S23)
2444 & +EP4*EP4*S13*S23)*S134
2446 & +EP2*(EP1+EP3+EP4)*2*S14*S24*S234)/(S*S13*S24*S134**2*S234)
2456 C-----------------------------------------------------------------------
2458 ENTRY HWH4J2(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2460 C-----------------------------------------------------------------------
2470 S=S12+S13+S14+S23+S24+S34
2472 HWH4J2=(S12*S14*(S24+S34)+S24*(S12*(S14+S34)+S13*(S14-S24)))
2474 & /(S14*S23*S13*S134)
2476 & +S12*(S+S34)*S124/(S24*S234*S14*S134)
2478 & -(S13*(2*(S12+S24)+S23)+S14**2)/(S134*S13*S14)
2480 & +S12*S123*S124/(2*S13*S24*S14*S23)
2486 & +4*((EP1*EP1*(S12*S134*S234-4*S23*S24*S34)
2488 & +EP1*EP2*(2*(2*S13*S234+S14*S123)*S24-S12*S134*(S+S12+S34))
2490 & +EP1*EP3*(S12*(4*S24*S34-S134*(S12+S14-S24))
2492 & -4*(S13*S24-S14*S23)*S24)
2494 & +EP1*EP4*(4*(S13+S14)*S23*S24-S12*S134*(S12+S13-S23))
2496 & +EP2*EP2*(S12*S134-4*S13*S24)*S134
2498 & +EP2*EP3*(4*S13*(S12+S23+S24)*S24-S12*S134*(S12-S14+S24))
2500 & -EP2*EP4*(4*(S12*(S14+S134)+S13*(S134-S234))*S24
2502 & +S12*(S12-S13+S23)*S134)
2504 & -EP3*EP3*4*S12*S14*S24
2506 & -EP3*EP4*2*S12*(2*S14*S24+S12*S134))*S234
2508 & +(EP1*(EP1*(S23+S24)+EP2*(S134-2*S))
2510 & -(EP1+EP2)*(EP3+EP4)*S12+EP2*EP2*(S13+S14))*2*S14*S24*S123)
2512 & /(2*S*S13*S14*S234*S23*S24*S134)
2522 C-----------------------------------------------------------------------
2524 ENTRY HWH4J4(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2526 C-----------------------------------------------------------------------
2532 S=S12+S13+S14+S23+S24+S34
2534 HWH4J4=-(S12*(S34*(3*(S+S34)+S12)-S134*S234-2*(S13*S24+S14*S23))
2536 & +(S14*S23-S13*S24)*(S13-S14+S24-S23))/(2*S134*S234*S34**2)
2538 & -(S12*(S134**2/2+2*S13*S14+S34*(S13+S14-S34))
2540 & +S34*((S13+S14)*(S23+S24)+S14*S24+S13*S23)
2542 & +(S13*S24-S14*S23)*(S14-S13))/(S34*S134)**2
2548 & +4*((-EP1*EP1*2*(S23+S24)*S34
2550 & -EP1*EP2*(S13*(S23+3*S24)+S14*(3*S23+S24)-(4*S12-S34)*S34)
2552 & +EP1*EP3*((2*S12-S24)*S34-(S13-S14)*S24)
2554 & +EP1*EP4*((2*S12-S23)*S34+(S13-S14)*S23)
2556 & -EP2*EP2*2*(S13+S14)*S34
2558 & +EP2*EP3*(2*S12*S34-S14*(S23-S24+S34))
2560 & +EP2*EP4*(2*S12*S34+S13*(S23-S24-S34))
2562 & +EP3*EP3*2*S14*S24
2564 & +EP3*EP4*2*(S12*S34-S13*S24-S14*S23)
2566 & +EP4*EP4*2*S13*S23)/(S*S134*S234*S34**2)
2568 & +(EP1*EP2*(S134*(S134+2*S34)+4*(S13*S14-S34**2))
2570 & +EP2*EP3*2*(2*S13*S34+S14*(S13-S14+S34))
2572 & +EP2*EP4*2*(2*S14*S34-S13*(S13-S14-S34)))
2574 & /(S*(S134*S34)**2))
2584 C-----------------------------------------------------------------------
2586 ENTRY HWH4J5(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2588 C-----------------------------------------------------------------------
2598 S=S12+S13+S14+S23+S24+S34
2600 HWH4J5=(3*S12*S34**2-3*S13*S24*S34+3*S12*S24*S34+3*S14*S23*S34-
2602 $ S13*S24**2-S12*S23*S34+6*S12*S14*S34+2*S12*S13*S34-
2604 $ 2*S12**2*S34+S14*S23*S24-3*S13*S23*S24-2*S13*S14*S24+
2606 $ 4*S12*S14*S24+2*S12*S13*S24+3*S14*S23**2+2*S14**2*S23+
2608 $ 2*S14**2*S12+2*S12**2*S14+6*S12*S14*S23-2*S12*S13**2-
2610 $ 2*S12**2*S13)/(2*S13*S134*S234*S34)+
2612 $ (2*S12*S34**2-2*S13*S24*S34+S12*S24*S34+4*S13*S23*S34+
2614 $ 4*S12*S14*S34+2*S12*S13*S34+2*S12**2*S34-S13*S24**2+
2616 $ 3*S14*S23*S24+4*S13*S23*S24-2*S13*S14*S24+4*S12*S14*S24+
2618 $ 2*S12*S13*S24+2*S14*S23**2+4*S13*S23**2+2*S13*S14*S23+
2620 $ 2*S12*S14*S23+4*S12*S13*S23+2*S12*S14**2+4*S12**2*S13+
2622 $ 4*S12*S13*S14+2*S12**2*S14)/(2*S13*S134*S24*S34)-
2624 $ (S12*S34**2-2*S14*S24*S34-2*S13*S24*S34-S14*S23*S34+
2626 $ S13*S23*S34+S12*S14*S34+2*S12*S13*S34-2*S14**2*S24-
2628 $ 4*S13*S14*S24-4*S13**2*S24-S14**2*S23-S13**2*S23+
2630 $ S12*S13*S14-S12*S13**2)/(S13*S34*S134**2)
2636 & +EP1*EP1*((S13-S14+S23-3*S24)*S34+(S134+S14+2*S34)*S234)
2640 & +EP1*EP2*((2*(S12-S24)+S34)*S134-S14*(4*S12+S14+3*S23)
2642 & +S13*(S13+S23)+S24*S34 )*S24*S134
2644 & -EP1*EP2*(((2*S12*S134+S13*(2*(S12+S14+S23)-S24+S34)
2646 & +S14*(S14-S23)+(2*S14-S34)*S234)*S234)*S134
2648 & + 4*S13**2*S24*S234)
2650 & +EP1*EP3*(S12*(2*S13-S134)+S13*(S24+2*S234)+S14*(3*S24-S234)
2652 & +S34*(S234-3*S24))*S24*S134
2654 & +EP1*EP4*((S12*(S13-S14+3*S34)-S23*(S13+3*S14-S34))*S24
2656 & -(S12*(S13+S134+2*S34)+2*S13*S24
2658 & +(S13-2*S14)*S23)*S234)*S134
2660 & +EP2*EP2*(S13*((2*S13+S34)*S234+S24*(S134-2*S34))
2662 & +2*S14*S134*(S24+S234))*S134
2666 & -EP2*EP3*(((S12*(S13+2*S14-S34)+S14*(S+2*S23-S34))*S24
2668 & +(S12*(S13+S134)+(S13+S24+2*S234)*S14
2670 & +2*S13*(2*S23+S34))*S234)*S134
2672 & +4*S13**2*S24*S234)
2674 & +EP2*EP4*(((S12*(S13-2*S134)+S13*(S+2*S23-3*S34))*S24
2676 & -((S-3*S13+S23+2*S24)*S13+2*S12*S14
2678 & +2*S14*(S23+2*S24))*S234)*S134-4*S13**2*S24*S234)
2680 & +EP3*EP3*2*(S13*S234+S14*S24)*S24*S134
2682 & +EP3*EP4*(2*(S12*S34-S13*S24-S14*S23)*S24
2684 & -(S12*S134+2*S13*S23)*S234)*S134
2686 & +EP4*EP4*2*(S12*S234+S23*S24)*S13*S134
2688 HWH4J5=HWH4J5+4*SUM/(S*S234*S134**2*S13*S34*S24)
2698 C-----------------------------------------------------------------------
2700 ENTRY HWH4J6(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2702 C-----------------------------------------------------------------------
2712 S=S12+S13+S14+S23+S24+S34
2714 HWH4J6=(S23*(S123*S234-S*S23)+S12*(S123*S124-S*S12))/(S13*S123)**2
2716 & -(S12*S34*(S234-2*S23)+S14*S23*(S234-2*S34)
2718 & -S13*S24*(S234+S13))/(S13**2*S123*S134)
2724 & +4*(-EP1*EP1*2*S23*S34
2726 & +EP1*EP2*((S12-S23)*S34-S13*(S24-S34))
2728 & +(EP1*EP3+EP2*EP4)*2*(S12*S34-S13*S24+S14*S23)
2730 & -EP1*EP4*(S13*S24-(3*(S13+S14)+S34)*S23)
2732 & -(EP1+EP2+EP3)*EP4*2
2734 & *(S12*(S13+S23)+(S12+S13)*S23)*S134/S123
2736 & +EP2*EP2*S13*(S14+S34)
2738 & +EP2*EP3*(S13*(S14-S24)-(S12-S23)*S14)
2740 & -EP3*EP3*2*S12*S14
2742 & -EP3*EP4*(S13*S24-(3*(S13+S34)+S14)*S12)
2744 & +EP4*EP4*(S12+S23)*S13)/(S*S134*S123*S13**2)
2754 C-----------------------------------------------------------------------
2756 ENTRY HWH4J7(S12,S13,S14,S23,S24,S34,EP1,EP2,EP3,EP4,ORIENT)
2758 C-----------------------------------------------------------------------
2768 S=S12+S13+S14+S23+S24+S34
2770 HWH4J7=((S12*S34+S13*S24-S14*S23)*(S13+S14+S23+S24)-2*S12*S24*S34)
2772 & /(S13*S134*S23*S123)
2774 & -S12*(S12*S-S123*S124)/(S123**2*S13*S23)
2776 & -(S13+S14)*(S23+S24)*S34/(S13*S134*S23*S234)
2782 & +4*(+2*(EP1+EP2)*(S23*EP1-S13*EP2)*S34*S134
2784 & -EP1*EP2*2*S34**2*S123
2786 & +EP1*EP3*(S123*(S23+S24)*S34+2*S134*(S13*S24-S14*S23))
2788 & +EP1*EP4*(S123*(S23+S24)*S34+2*S12**2*S134*S234/S123
2790 & +2*S134*(S24*(S13-S12)-S23*(S12+S14)))
2792 & +EP2*EP3*(2*(S12*S34+S13*S24-S14*S23)*S134
2794 & +S123*(S13+S14)*S34)
2796 & +EP2*EP4*(S123*(S13+S14)*S34+2*S12**2*S234*S134/S123
2798 & -2*S134*(S12*S234-S13*S24+S14*S23))
2800 & -EP3*EP3*S12*(2*S24*S134+S123*S34)
2802 & +EP3*EP4*2*S12*(S134*(S23-S24)-S34*S123+S12*S134*S234/S123)
2804 & +EP4*EP4*S12*(2*S23*S134-S123*S34))
2806 & /(S*S13*S23*S123*S134*S234)