Update to 6.5-10
[u/mrichter/AliRoot.git] / HERWIG / herwig6510.f
1 C  HERWIG---AliRoot/HERWIG
2 C-----------------------------------------------------------------------
3 C                           H E R W I G
4 C
5 C            a Monte Carlo event generator for simulating
6 C        +---------------------------------------------------+
7 C        | Hadron Emission Reactions With Interfering Gluons |
8 C        +---------------------------------------------------+
9 C I.G. Knowles(*), G. Marchesini(+), M.H.Seymour($,&) and B.R. Webber(#)
10 C-----------------------------------------------------------------------
11 C with Minimal Supersymmetric Standard Model Matrix Elements by
12 C                  S. Moretti(") and K. Odagiri(^)
13 C-----------------------------------------------------------------------
14 C R parity violating Supersymmetric Decays and Matrix Elements by
15 C                          P. Richardson(X)
16 C-----------------------------------------------------------------------
17 C matrix element corrections to top decay and Drell-Yan type processes
18 C                         by G. Corcella(&)
19 C-----------------------------------------------------------------------
20 C Deep Inelastic Scattering and Heavy Flavour Electroproduction by
21 C                  G. Abbiendi(@) and L. Stanco(%)
22 C-----------------------------------------------------------------------
23 C and Jet Photoproduction in Lepton-Hadron Collisions by J. Chyla(~)
24 C-----------------------------------------------------------------------
25 C(*)  Department of Physics & Astronomy, University of Edinburgh
26 C(+)  Dipartimento di Fisica, Universita di Milano-Bicocca
27 C($)  School of Physics & Astronomy, University of Manchester
28 C(&)  Theory Physics Group, CERN
29 C(#)  Cavendish Laboratory, Cambridge
30 C(")  School of Physics & Astronomy, Southampton
31 C(^)  Academia Sinica, Taiwan
32 C(X)  Institute of Particle Physics Phenomenology, University of Durham
33 C(@)  Dipartimento di Fisica, Universita di Bologna
34 C(%)  Dipartimento di Fisica, Universita di Padova
35 C(~)  Institute of Physics, Prague
36 C-----------------------------------------------------------------------
37 C                  Version 6.510 - 31st October 2005
38 C-----------------------------------------------------------------------
39 C Main references:
40 C
41 C    G.Corcella, I.G.Knowles, G.Marchesini, S.Moretti, K.Odagiri,
42 C    P.Richardson, M.H.Seymour and B.R.Webber, JHEP 0101 (2001) 010
43 C
44 C    G.Marchesini,  B.R.Webber,  G.Abbiendi,  I.G.Knowles,  M.H.Seymour,
45 C    and L.Stanco, Computer Physics Communications 67 (1992) 465.
46 C-----------------------------------------------------------------------
47 C Please see the official HERWIG information page:
48 C    http://hepwww.rl.ac.uk/theory/seymour/herwig/
49 C-----------------------------------------------------------------------
50 CDECK  ID>, CIRCEE.
51 *CMZ :-        -03/07/01  17.07.47  by  Bryan Webber
52 *-- Author :    Bryan Webber
53 C-----------------------------------------------------------------------
54       FUNCTION CIRCEE (X1, X2)
55 C-----------------------------------------------------------------------
56 C     DUMMY FUNCTION: DELETE AND SET CIRCOP NON-ZERO
57 C     IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
58 C-----------------------------------------------------------------------
59       IMPLICIT NONE
60       DOUBLE PRECISION CIRCEE, X1, X2
61       WRITE (6,10)
62    10 FORMAT(/10X,'CIRCEE CALLED BUT NOT LINKED')
63       CIRCEE = 0.0D0
64       STOP
65       END
66 CDECK  ID>, CIRCES.
67 *CMZ :-        -03/07/01  17.07.47  by  Bryan Webber
68 *-- Author :    Bryan Webber
69 C-----------------------------------------------------------------------
70       SUBROUTINE CIRCES (XX1M, XX2M, XROOTS, XACC, XVER, XREV, XCHAT)
71 C-----------------------------------------------------------------------
72 C     DUMMY SUBROUTINE: DELETE AND SET CIRCOP NON-ZERO
73 C     IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
74 C-----------------------------------------------------------------------
75       IMPLICIT NONE
76       DOUBLE PRECISION XX1M, XX2M, XROOTS
77       INTEGER XACC, XVER, XREV, XCHAT
78       WRITE (6,10)
79    10 FORMAT(/10X,'CIRCES CALLED BUT NOT LINKED')
80       STOP
81       END
82 CDECK  ID>, CIRCGG.
83 *CMZ :-        -03/07/01  17.07.47  by  Bryan Webber
84 *-- Author :    Bryan Webber
85 C-----------------------------------------------------------------------
86       FUNCTION CIRCGG (X1, X2)
87 C-----------------------------------------------------------------------
88 C     DUMMY FUNCTION: DELETE AND SET CIRCOP NON-ZERO
89 C     IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
90 C-----------------------------------------------------------------------
91       IMPLICIT NONE
92       DOUBLE PRECISION CIRCGG, X1, X2
93       WRITE (6,10)
94    10 FORMAT(/10X,'CIRCGG CALLED BUT NOT LINKED')
95       CIRCGG = 0.0D0
96       STOP
97       END
98 CDECK  ID>, DECADD.
99 *CMZ :-        -28/01/92  12.34.44  by  Mike Seymour
100 *-- Author :    Luca Stanco
101 C-----------------------------------------------------------------------
102       SUBROUTINE DECADD(LOGI)
103 C-----------------------------------------------------------------------
104 C     DUMMY SUBROUTINE: DELETE AND SET BDECAY='CLEO'
105 C     IN MAIN PROGRAM IF YOU USE CLEO DECAY PACKAGE
106 C-----------------------------------------------------------------------
107       IMPLICIT NONE
108       LOGICAL LOGI
109       WRITE (6,10)
110    10 FORMAT(/10X,'DECADD CALLED BUT NOT LINKED')
111       STOP
112       END
113 CDECK  ID>, DEXAY.
114 *CMZ :-        -17/10/01  10.03.37  by  Peter Richardson
115 *-- Author :    Peter Richardson
116 C-----------------------------------------------------------------------
117       SUBROUTINE DEXAY(IMODE,POL)
118 C-----------------------------------------------------------------------
119 C     DUMMY SUBROUTINE: DELETE AND SET TAUDEC='TAUOLA'
120 C     IN MAIN PROGRAM IF YOU USE TAUOLA DECAY PACKAGE
121 C-----------------------------------------------------------------------
122       IMPLICIT NONE
123       INTEGER IMODE
124       REAL POL(4)
125       WRITE (6,10)
126    10 FORMAT(/10X,'DEXAY CALLED BUT NOT LINKED')
127       STOP
128       END
129 CDECK  ID>, EUDINI.
130 *CMZ :-        -28/01/92  12.34.44  by  Mike Seymour
131 *-- Author :    Luca Stanco
132 C-----------------------------------------------------------------------
133       SUBROUTINE EUDINI
134 C-----------------------------------------------------------------------
135 C     DUMMY SUBROUTINE: DELETE AND SET BDECAY='EURO'
136 C     IN MAIN PROGRAM IF YOU USE EURODEC DECAY PACKAGE
137 C-----------------------------------------------------------------------
138       IMPLICIT NONE
139       WRITE (6,10)
140    10 FORMAT(/10X,'EUDINI CALLED BUT NOT LINKED')
141       STOP
142       END
143 CDECK  ID>, FILHEP.
144 *CMZ :-        -17/10/01  09:42:21  by  Peter Richardson
145 *-- Author :    Martin W. Gruenewald
146 C-----------------------------------------------------------------------
147       SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
148 C ----------------------------------------------------------------------
149 C this subroutine fills one entry into the HEPEVT common
150 C and updates the information for affected mother entries
151 C used by TAUOLA
152 C
153 C written by Martin W. Gruenewald (91/01/28)
154 C ----------------------------------------------------------------------
155       INCLUDE 'HERWIG65.INC'
156       LOGICAL QEDRAD
157       COMMON /PHORAD/ QEDRAD(NMXHEP)
158       INTEGER N,IHEP,IST,ID,JMO1,JMO2,JDA1,JDA2,I,IP
159       REAL PINV
160       LOGICAL PHFLAG
161       REAL*4 P4(4)
162 C
163 C check address mode
164       IF (N.EQ.0) THEN
165 C append mode
166         IHEP=NHEP+1
167       ELSE IF (N.GT.0) THEN
168 C absolute position
169         IHEP=N
170       ELSE
171 C relative position
172         IHEP=NHEP+N
173       END IF
174 C check on IHEP
175       IF ((IHEP.LE.0).OR.(IHEP.GT.NMXHEP)) RETURN
176 C add entry
177       NHEP=IHEP
178       ISTHEP(IHEP)=IST
179       IDHEP(IHEP)=ID
180       JMOHEP(1,IHEP)=JMO1
181       IF(JMO1.LT.0)JMOHEP(1,IHEP)=JMOHEP(1,IHEP)+IHEP
182       JMOHEP(2,IHEP)=JMO2
183       IF(JMO2.LT.0)JMOHEP(2,IHEP)=JMOHEP(2,IHEP)+IHEP
184       JDAHEP(1,IHEP)=JDA1
185       JDAHEP(2,IHEP)=JDA2
186       DO I=1,4
187         PHEP(I,IHEP)=P4(I)
188 C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
189         VHEP(I,IHEP)=0.0
190       END DO
191       PHEP(5,IHEP)=PINV
192 C FLAG FOR PHOTOS...
193       QEDRAD(IHEP)=PHFLAG
194 C update process:
195       DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
196         IF(IP.GT.0)THEN
197 C if there is a daughter at IHEP, mother entry at IP has decayed
198           IF(ISTHEP(IP).EQ.1)ISTHEP(IP)=2
199 C and daughter pointers of mother entry must be updated
200           IF(JDAHEP(1,IP).EQ.0)THEN
201             JDAHEP(1,IP)=IHEP
202             JDAHEP(2,IP)=IHEP
203           ELSE
204             JDAHEP(2,IP)=MAX(IHEP,JDAHEP(2,IP))
205           END IF
206         END IF
207       END DO
208       END
209 CDECK  ID>, FRAGMT.
210 *CMZ :-        -28/01/92  12.34.44  by  Mike Seymour
211 *-- Author :    Luca Stanco
212 C-----------------------------------------------------------------------
213       SUBROUTINE FRAGMT(I,J,K)
214 C-----------------------------------------------------------------------
215 C     DUMMY SUBROUTINE: DELETE AND SET BDECAY='EURO'
216 C     IN MAIN PROGRAM IF YOU USE EURODEC DECAY PACKAGE
217 C-----------------------------------------------------------------------
218       IMPLICIT NONE
219       INTEGER I,J,K
220       WRITE (6,10)
221    10 FORMAT(/10X,'FRAGMT CALLED BUT NOT LINKED')
222       STOP
223       END
224 CDECK  ID>, HVCBVI.
225 *CMZ :-        -28/01/92  12.34.44  by  Mike Seymour
226 *-- Author :    Mike Seymour
227 C-----------------------------------------------------------------------
228       SUBROUTINE HVCBVI
229 C-----------------------------------------------------------------------
230 C     DUMMY ROUTINE: DELETE IF YOU LINK TO BARYON NUMBER VIOLATN PACKAGE
231 C-----------------------------------------------------------------------
232       IMPLICIT NONE
233       WRITE (6,10)
234    10 FORMAT(/10X,'HVCBVI CALLED BUT NOT LINKED')
235       STOP
236       END
237 CDECK  ID>, HVHBVI.
238 *CMZ :-        -28/01/92  12.34.44  by  Mike Seymour
239 *-- Author :    Mike Seymour
240 C-----------------------------------------------------------------------
241       SUBROUTINE HVHBVI
242 C-----------------------------------------------------------------------
243 C     DUMMY ROUTINE: DELETE IF YOU LINK TO BARYON NUMBER VIOLATN PACKAGE
244 C-----------------------------------------------------------------------
245       IMPLICIT NONE
246       WRITE (6,10)
247    10 FORMAT(/10X,'HERBVI CALLED BUT NOT LINKED')
248       STOP
249       END
250 CDECK  ID>, HWBAZF.
251 *CMZ :-        -26/04/91  11.11.54  by  Bryan Webber
252 *-- Author :    Ian Knowles
253 C-----------------------------------------------------------------------
254       SUBROUTINE HWBAZF(IPAR,JPAR,VEC1,VEC2,VEC3,VEC)
255 C-----------------------------------------------------------------------
256 C     Azimuthal correlation functions for Collins' algorithm,
257 C     see I.G.Knowles, Comp. Phys. Comm. 58 (90) 271 for notation.
258 C-----------------------------------------------------------------------
259       INCLUDE 'HERWIG65.INC'
260       DOUBLE PRECISION Z1,Z2,DOT12,DOT23,DOT31,TR,FN(7),VEC1(2),VEC2(2),
261      & VEC3(2),VEC(2)
262       INTEGER IPAR,JPAR
263       LOGICAL GLUI,GLUJ
264       IF (.NOT.AZSPIN) RETURN
265       Z1=PPAR(4,JPAR)/PPAR(4,IPAR)
266       Z2=1.-Z1
267       GLUI=IDPAR(IPAR).EQ.13
268       GLUJ=IDPAR(JPAR).EQ.13
269       IF (GLUI) THEN
270          IF (GLUJ) THEN
271 C           Branching: g--->gg
272             FN(2)=Z2/Z1
273             FN(3)=1./FN(2)
274             FN(4)=Z1*Z2
275             FN(1)=FN(2)+FN(3)+FN(4)
276             FN(5)=FN(2)+2.*Z1
277             FN(6)=FN(3)+2.*Z2
278             FN(7)=FN(4)-2.
279          ELSE
280 C           Branching: g--->qqbar
281             FN(1)=(Z1*Z1+Z2*Z2)/2.
282             FN(2)=0.
283             FN(3)=0.
284             FN(4)=-Z1*Z2
285             FN(5)=-(2.*Z1-1.)/2.
286             FN(6)=-FN(5)
287             FN(7)=FN(1)
288          ENDIF
289       ELSE
290          IF (GLUJ) THEN
291 C           Branching: q--->gq
292             FN(1)=(1.+Z2*Z2)/(2.*Z1)
293             FN(2)=Z2/Z1
294             FN(3)=0.
295             FN(4)=0.
296             FN(5)=FN(1)
297             FN(6)=(1.+Z2)/2.
298             FN(7)=-FN(6)
299          ELSE
300 C           Branching: q--->qg
301             FN(1)=(1.+Z1*Z1)/(2.*Z2)
302             FN(2)=0.
303             FN(3)=Z1/Z2
304             FN(4)=0.
305             FN(5)=(1.+Z1)/2.
306             FN(6)=FN(1)
307             FN(7)=-FN(5)
308          ENDIF
309       ENDIF
310       DOT12=VEC1(1)*VEC2(1)+VEC1(2)*VEC2(2)
311       DOT23=VEC2(1)*VEC3(1)+VEC2(2)*VEC3(2)
312       DOT31=VEC3(1)*VEC1(1)+VEC3(2)*VEC1(2)
313       TR=1./(FN(1)+FN(2)*DOT23+FN(3)*DOT31+FN(4)*DOT12)
314       VEC(1)=((FN(2)+FN(5)*DOT23)*VEC1(1)
315      &       +(FN(3)+FN(6)*DOT31)*VEC2(1)
316      &       +(FN(4)+FN(7)*DOT12)*VEC3(1))*TR
317       VEC(2)=((FN(2)+FN(5)*DOT23)*VEC1(2)
318      &       +(FN(3)+FN(6)*DOT31)*VEC2(2)
319      &       +(FN(4)+FN(7)*DOT12)*VEC3(2))*TR
320       END
321 CDECK  ID>, HWBCON.
322 *CMZ :-        -11/10/01  12.01.52  by  Peter Richardson
323 *-- Author :    Bryan Webber
324 C-----------------------------------------------------------------------
325       SUBROUTINE HWBCON
326 C-----------------------------------------------------------------------
327 C     MAKES COLOUR CONNECTIONS BETWEEN JETS
328 C     MODIFIED 12/10/97 BY BRW FOR SUSY PROCESSES
329 C     MODIFIED 11/01/01 BY PR  FOR SPIN CORRELATIONS(PROBLEM WITH ORDER
330 C                                                    OF DECAYS)
331 C     NEW VARAIBLE BACK TO ALLOW CODE TO SEARCH DOWN CHAIN
332 C-----------------------------------------------------------------------
333       INCLUDE 'HERWIG65.INC'
334       INTEGER IHEP,IST,ID,JC,KC,JD,JHEP,LHEP,ID2,NTRY,KHEP
335       LOGICAL BACK
336       IF (IERROR.NE.0) RETURN
337       IF(.NOT.RPARTY) THEN
338         CALL HWBRCN
339         RETURN
340       ENDIF
341       DO 20 IHEP=1,NHEP
342       BACK = .FALSE.
343       IST=ISTHEP(IHEP)
344 C---LOOK FOR PARTONS WITHOUT COLOUR MOTHERS
345       IF (IST.LT.145.OR.IST.GT.152) GOTO 20
346  51   IF (JMOHEP(2,IHEP).EQ.0.OR.BACK.OR.
347      &     ISTHEP(JMOHEP(2,IHEP)).EQ.155) THEN
348 C---FIND COLOUR-CONNECTED PARTON
349         IF(BACK) GOTO 52
350         IF(JMOHEP(2,IHEP).EQ.0) THEN
351           JC=JMOHEP(1,IHEP)
352           IF (IST.NE.152) JC=JMOHEP(1,JC)
353           JC =JMOHEP(2,JC)
354         ELSE
355           JC = JMOHEP(2,IHEP)
356           JHEP = JC
357         ENDIF
358         IF (JC.EQ.0) THEN
359           CALL HWWARN('HWBCON',51)
360           GOTO 20
361         ENDIF
362 C---FIND SPECTATOR WHEN JC IS DECAYED HEAVY QUARK OR SUSY PARTICLE
363  52       IF (ISTHEP(JC).EQ.155.OR.BACK) THEN
364           IF (IDHEP(JMOHEP(1,JC)).EQ.94.OR.BACK) THEN
365 C---DECAYED BEFORE HADRONIZING
366             IF(BACK.OR.(JMOHEP(2,IHEP).NE.0.AND.
367      &                  ISTHEP(JMOHEP(2,IHEP)).EQ.155)) GOTO 53
368             JHEP=JMOHEP(2,JC)
369 C--new bit to try and fix the problems for spin correlations
370 C--move one step further up the tree and hope this helps
371             IF (JHEP.EQ.0) THEN
372               NTRY = 0
373  1            NTRY = NTRY+1
374               JC   = JMOHEP(1,JC)
375               JHEP = JMOHEP(2,JC)
376               IF(JHEP.NE.0.AND.ISTHEP(JHEP).EQ.155)
377      &             JHEP = JMOHEP(2,JHEP)
378               IF(JHEP.EQ.0.AND.NTRY.LT.NHEP) GOTO 1
379               IF(NHEP.EQ.NTRY) GOTO 20
380             ENDIF
381  53         ID=IDHW(JHEP)
382             IF (ISTHEP(JHEP).EQ.155) THEN
383 C---SPECIAL FOR GLUINO DECAYS
384               IF (ID.EQ.449) THEN
385                 ID=IDHW(JC)
386 C---N.B. WILL NEED MODS WHEN SUSY PARTICLES CAN SHOWER
387                 IF (ID.EQ.449.OR.ID.EQ.13.OR.
388      &             (ID.GE.401.AND.ID.LE.406).OR.
389      &             (ID.GE.413.AND.ID.LE.418).OR.
390      &             ID.LE.6.OR.(ID.GE.115.AND.ID.LE.120)) THEN
391 C---LOOK FOR ANTI(S)QUARK OR GLUON
392                   DO KC=JDAHEP(1,JHEP),JDAHEP(2,JHEP)
393                     ID=IDHW(KC)
394                     IF ((ID.GE.  7.AND.ID.LE. 13).OR.
395      &                  (ID.GE.407.AND.ID.LE.412).OR.
396      &                  (ID.GE.419.AND.ID.LE.424)) GOTO 5
397                   ENDDO
398                 ELSE
399 C---LOOK FOR (S)QUARK OR GLUON
400                   DO KC=JDAHEP(1,JHEP),JDAHEP(2,JHEP)
401                     ID=IDHW(KC)
402                     IF (ID.LE.  6.OR. ID.EQ. 13.OR.
403      &                 (ID.GE.401.AND.ID.LE.406).OR.
404      &                 (ID.GE.413.AND.ID.LE.418)) GOTO 5
405                   ENDDO
406                 ENDIF
407 C---COULDNT FIND ONE
408                 CALL HWWARN('HWBCON',101)
409                 GOTO 999
410     5           JC=KC
411               ELSE
412 C--PR MOD 30/6/99 should fix HWCFOR 104 errors
413                 ID2 = IDHW(IHEP)
414                 IF(IDHW(JDAHEP(1,JHEP)).EQ.449.AND.
415      &             (ID2.LE.6.OR.(ID2.GE.115.AND.ID2.LE.120).OR.
416      &             (ID2.GE.401.AND.ID2.LE.406).OR.ID2.EQ.13.OR.
417      &             (ID2.GE.413.AND.ID2.LE.418).OR.ID2.EQ.449)) THEN
418                   JC = JDAHEP(1,JHEP)
419                 ELSE
420 C--modifcation for top ME correction (modified for additional photon radiation)
421                   IF(IDHW(JHEP).EQ.6) THEN
422                     JC = JDAHEP(1,JHEP)+1
423                   ELSE
424                     JC = JDAHEP(1,JHEP)+1
425                     IF(IDHW(JDAHEP(1,JHEP)+2).EQ.13) JC=JC+1
426                   ENDIF
427                 ENDIF
428               ENDIF
429             ELSEIF (ID.EQ.6.OR.ID.EQ.12.OR.
430      &      (ID.GE.209.AND.ID.LE.218).OR.
431      &      (ID.GE.401.AND.ID.LE.424).OR.ID.EQ.449) THEN
432 C Wait for partner heavy quark to decay
433 C              RETURN
434 C---N.B. MAY BE A PROBLEM HERE
435               GOTO 20
436             ELSE
437               JMOHEP(2,IHEP)=JHEP
438               JDAHEP(2,JHEP)=IHEP
439               GOTO 20
440             ENDIF
441           ELSE
442             JC=JMOHEP(2,JC)
443           ENDIF
444         ENDIF
445         JC=JDAHEP(1,JC)
446         JD=JDAHEP(2,JC)
447 C---SEARCH IN CORRESPONDING JET
448         IF (JD.LT.JC) JD=JC
449         LHEP=0
450         DO 10 JHEP=JC,JD
451         IF (ISTHEP(JHEP).LT.145.OR.ISTHEP(JHEP).GT.152) GOTO 10
452         IF (JDAHEP(2,JHEP).EQ.IHEP) LHEP=JHEP
453         IF (JDAHEP(2,JHEP).NE.0) GOTO 10
454 C---JOIN IHEP AND JHEP
455         ID=IDHW(JHEP)
456         JMOHEP(2,IHEP)=JHEP
457         JDAHEP(2,JHEP)=IHEP
458         GOTO 20
459    10   CONTINUE
460         IF (LHEP.NE.0) THEN
461           JMOHEP(2,IHEP)=LHEP
462         ELSE
463 C--search down the tree
464           DO 50 KHEP=JC,JD
465           IF(ISTHEP(KHEP).EQ.3.AND.ISTHEP(JDAHEP(1,KHEP)).EQ.155) THEN
466             JHEP = JDAHEP(1,KHEP)
467             BACK = .TRUE.
468             GOTO 51
469           ENDIF
470  50       CONTINUE
471 C---DIDN'T FIND PARTNER OF IHEP YET
472 C          CALL HWWARN('HWBCON',52)
473 C          GOTO 20
474         ENDIF
475       ENDIF
476   20  CONTINUE
477 C---BREAK COLOUR CONNECTIONS WITH PHOTONS
478       IHEP=1
479   30  IF (IHEP.LE.NHEP) THEN
480         IF (IDHW(IHEP).EQ.59 .AND. ISTHEP(IHEP).EQ.149) THEN
481 C  BRW FIX 13/03/99
482           IF (JMOHEP(2,IHEP).NE.0) THEN
483             IF (JDAHEP(2,JMOHEP(2,IHEP)).EQ.IHEP)
484      &        JDAHEP(2,JMOHEP(2,IHEP))=JDAHEP(2,IHEP)
485           ENDIF
486 C  END FIX
487           IF (JDAHEP(2,IHEP).NE.0) THEN
488             IF (JMOHEP(2,JDAHEP(2,IHEP)).EQ.IHEP)
489      &        JMOHEP(2,JDAHEP(2,IHEP))=JMOHEP(2,IHEP)
490           ENDIF
491           JMOHEP(2,IHEP)=IHEP
492           JDAHEP(2,IHEP)=IHEP
493         ENDIF
494         IHEP=IHEP+1
495         GOTO 30
496       ENDIF
497  999  RETURN
498       END
499 CDECK  ID>, HWBDED.
500 *CMZ :-        -22/04/96  13.54.08  by  Mike Seymour
501 *-- Author :    Mike Seymour
502 C-----------------------------------------------------------------------
503       SUBROUTINE HWBDED(IOPT)
504 C     FILL MISSING AREA OF DALITZ PLOT WITH 3-JET AND 2-JET+GAMMA EVENTS
505 C     IF (IOPT.EQ.1) SET UP EVENT RECORD
506 C     IF (IOPT.EQ.2) CLEAN UP EVENT RECORD AFTER SHOWERING
507 C
508 C********MODIFIED 13/11/00 BY BRW TO ALLOW MULTIPLE APPLICATION IN
509 C*******SAME EVENT (FOR WW AND ZZ) N.B. NO CLEANUP CALLS FOR THESE!
510 C-----------------------------------------------------------------------
511       INCLUDE 'HERWIG65.INC'
512       DOUBLE PRECISION HWBVMC,HWRGEN,HWUALF,HWUSQR,X(3),W,WMAX,WSUM,
513      & X1MIN,X1MAX,X2MIN,X2MAX,QSCALE,GAMFAC,GLUFAC,R(3,3),CS,SN,M(3),
514      & E(3),LAMBDA,A,B,C,PTSQ,EM,P1(5),P2(5),PVRT(4),EPS,MASDEP
515       INTEGER ID,ID3,EMIT,NOEMIT,IEVT,IHEP,JHEP,KHEP,ICMF,IOPT,IEDT(3),
516      & I,NDEL,LHEP,IP,JP,KP,IDUN
517       EXTERNAL HWBVMC,HWRGEN,HWUALF,HWUSQR
518       SAVE X,WMAX,P1,P2
519       SAVE WSUM,     X1MIN,X1MAX,EMIT,ICMF,IEVT
520       DATA WSUM,WMAX,X1MIN,X1MAX,EMIT,ICMF,IEVT
521      & /0.994651D0,1.84096D0,0.0D0,0.773459D0,3*0.0D0/
522       LAMBDA(A,B,C)=(A**2+B**2+C**2-2*A*B-2*B*C-2*C*A)/(4*A)
523       IF (IOPT.EQ.1) THEN
524 C---FIND AN UNTREATED CMF
525         IF (IEVT.EQ.NEVHEP+NWGTS) RETURN
526         IEVT=0
527         ICMF=0
528  5      IDUN=ICMF
529         DO 10 IHEP=IDUN+1,NHEP
530  10       IF (ICMF.EQ.IDUN .AND. ISTHEP(IHEP).EQ.110 .AND.
531      &    JDAHEP(2,IHEP).EQ.JDAHEP(1,IHEP)+1) ICMF=IHEP
532         IF (ICMF.EQ.IDUN) RETURN
533         EM=PHEP(5,ICMF)
534         IF (EM.LT.2*HWBVMC(1)) GOTO 5
535 C---ONLY APPLY THE CORRECTION TO HADRONIC DECAYS
536         IF (IDHW(JDAHEP(1,ICMF)).GT.12) GOTO 5
537 C---GENERATE X1,X2 ACCORDING TO 1/((1-X1)*(1-X2))
538  100    CONTINUE
539 C---CHOOSE X1
540         X(1)=1-(1-X1MAX)*((1-X1MIN)/(1-X1MAX))**HWRGEN(0)
541 C---CHOOSE X2
542         X2MIN=MAX(X(1),1-X(1))
543         X2MAX=(4*X(1)-3+2*DREAL(  DCMPLX(  X(1)**3+135*(X(1)-1)**3,
544      &    3*HWUSQR(3*(128*X(1)**4-368*X(1)**3+405*X(1)**2-216*X(1)+54))*
545      &    (X(1)-1)  )**(1./3)  ))/3
546         IF (X2MAX.GE.ONE.OR.X2MIN.GE.ONE.OR.X2MAX.LE.X2MIN) GOTO 100
547         X(2)=1-(1-X2MAX)*((1-X2MIN)/(1-X2MAX))**HWRGEN(1)
548 C---CALCULATE WEIGHT
549         W=2 * LOG((1-X1MIN)/(1-X1MAX))*LOG((1-X2MIN)/(1-X2MAX)) *
550      &    (X(1)**2+X(2)**2)
551 C---GENERATE UNWEIGHTED (X1,X2) PAIRS (EFFICIENCY IS ~50%)
552         IF (WMAX*HWRGEN(2).GT.W) GOTO 100
553 C---SYMMETRIZE X1,X2
554         X(3)=2-X(1)-X(2)
555         IF (HWRGEN(5).GT.HALF) THEN
556           X(1)=X(2)
557           X(2)=2-X(3)-X(1)
558         ENDIF
559 C---CHOOSE WHICH PARTON WILL EMIT
560         EMIT=1
561         IF (HWRGEN(6).LT.X(1)**2/(X(1)**2+X(2)**2)) EMIT=2
562         NOEMIT=3-EMIT
563         IHEP=JDAHEP(  EMIT,ICMF)
564         JHEP=JDAHEP(NOEMIT,ICMF)
565 C---PREFACTORS FOR GAMMA AND GLUON CASES
566         QSCALE=HWUSQR((1-X(1))*(1-X(2))*(1-X(3)))*EM/X(NOEMIT)
567         ID=IDHW(JDAHEP(1,ICMF))
568         GAMFAC=ALPFAC*ALPHEM*ICHRG(ID)**2/(18*PIFAC)
569         GLUFAC=0
570         IF (QSCALE.GT.HWBVMC(13))
571      &    GLUFAC=CFFAC/(2*PIFAC)*HWUALF(1,QSCALE)
572 C---SWITCH OFF PHOTON EMISSION IN W DECAYS (THE M-E DOES NOT FACTORIZE)
573         IF (ICHRG(IDHW(ICMF)).NE.0) GAMFAC=0
574 C---IN FRACTION FAC*WSUM OF EVENTS ADD A GAMMA/GLUON
575         IF     (GAMFAC*WSUM .GT. HWRGEN(3)) THEN
576           ID3=59
577         ELSEIF (GLUFAC*WSUM .GT. HWRGEN(4)) THEN
578           ID3=13
579         ELSE
580           EMIT=0
581           GOTO 5
582         ENDIF
583 C---CHECK INFRA-RED CUT-OFF FOR GAMMA/GLUON
584         M(EMIT)=PHEP(5,IHEP)+VQCUT
585         M(NOEMIT)=PHEP(5,JHEP)+VQCUT
586         M(3)=HWBVMC(ID3)
587         E(1)=HALF*EM*(X(1)+(M(1)**2-M(2)**2-M(3)**2)/EM**2)
588         E(2)=HALF*EM*(X(2)+(M(2)**2-M(3)**2-M(1)**2)/EM**2)
589         E(3)=EM-E(1)-E(2)
590         PTSQ=-LAMBDA(E(NOEMIT)**2-M(NOEMIT)**2,E(3)**2-M(3)**2,
591      &    E(EMIT)**2-M(EMIT)**2)
592         IF (PTSQ.LE.ZERO .OR.
593      $       E(1).LE.M(1).OR.E(2).LE.M(2).OR.E(3).LE.M(3)) THEN
594           EMIT=0
595           GOTO 5
596         ENDIF
597 C---CALCULATE MASS-DEPENDENT SUPRESSION
598         IF (MOD(IPROC,10).GT.0) THEN
599           EPS=(RMASS(ID)/EM)**2
600           MASDEP=X(1)**2+X(2)**2
601      $         -4*EPS*X(3)-2*EPS*((1-X(2))/(1-X(1))+(1-X(1))/(1-X(2)))
602      $         -4*EPS**2*X(3)**2/((1-X(1))*(1-X(2)))
603           IF (MASDEP.LT.HWRGEN(7)*(X(1)**2+X(2)**2)) THEN
604             EMIT=0
605             GOTO 5
606           ENDIF
607         ENDIF
608 C---STORE OLD MOMENTA
609         CALL HWVEQU(5,PHEP(1,JDAHEP(1,ICMF)),P1)
610         CALL HWVEQU(5,PHEP(1,JDAHEP(2,ICMF)),P2)
611 C---GET THE NON-EMITTING PARTON'S CMF DIRECTION
612         CALL HWULOF(PHEP(1,ICMF),PHEP(1,JHEP),PHEP(1,JHEP))
613         CALL HWRAZM(ONE,CS,SN)
614         CALL HWUROT(PHEP(1,JHEP),CS,SN,R)
615         M(EMIT)=PHEP(5,IHEP)
616         M(NOEMIT)=PHEP(5,JHEP)
617         M(3)=RMASS(ID3)
618         KHEP=JDAHEP(2,ICMF)
619         LHEP=KHEP+1
620         IF (NHEP.GT.KHEP) THEN
621 C---MOVE UP REST OF EVENT
622            DO IP=NHEP,LHEP,-1
623               JP=IP+1
624               ISTHEP(JP)= ISTHEP(IP)
625               IDHW(JP)=IDHW(IP)
626               IDHEP(JP)=IDHEP(IP)
627               KP=JMOHEP(1,IP)
628               IF (KP.GT.KHEP) THEN
629                  KP=KP+1
630               ELSE
631                  IF (JDAHEP(1,KP).EQ.IP) JDAHEP(1,KP)=JP
632                  IF (JDAHEP(2,KP).EQ.IP) JDAHEP(2,KP)=JP
633               ENDIF
634               JMOHEP(1,JP)=KP
635               KP=JMOHEP(2,IP)
636               IF (KP.GT.KHEP) KP=KP+1
637               JMOHEP(2,JP)=KP
638               KP=JDAHEP(1,IP)
639               IF (KP.GT.KHEP) KP=KP+1
640               JDAHEP(1,JP)=KP
641               KP=JDAHEP(2,IP)
642               IF (KP.GT.KHEP) KP=KP+1
643               JDAHEP(2,JP)=KP
644               CALL HWVEQU(5,PHEP(1,IP),PHEP(1,JP))
645               CALL HWVEQU(4,VHEP(1,IP),VHEP(1,JP))
646            ENDDO
647         ENDIF
648 C---REORDER ENTRIES: IHEP=EMITTER, JHEP=NON-EMITTER, KHEP=EMITTED
649         NHEP=NHEP+1
650         IF (IDHW(IHEP).LT.IDHW(JHEP)) THEN
651           IHEP=JDAHEP(1,ICMF)
652           JHEP=LHEP
653         ELSE
654           IHEP=LHEP
655           JHEP=JDAHEP(1,ICMF)
656         ENDIF
657 C---SET UP MOMENTA
658         PHEP(5,JHEP)=M(NOEMIT)
659         PHEP(5,IHEP)=M(EMIT)
660         PHEP(5,KHEP)=M(3)
661         PHEP(4,JHEP)=HALF*EM*(X(NOEMIT)+
662      &                  (M(NOEMIT)**2-M(EMIT)**2-M(3)**2)/EM**2)
663         PHEP(4,IHEP)=HALF*EM*(X(EMIT)+
664      &                  (M(EMIT)**2-M(NOEMIT)**2-M(3)**2)/EM**2)
665         PHEP(4,KHEP)=EM-PHEP(4,IHEP)-PHEP(4,JHEP)
666         PHEP(3,JHEP)=HWUSQR(PHEP(4,JHEP)**2-PHEP(5,JHEP)**2)
667         PHEP(3,IHEP)=( (PHEP(4,KHEP)**2-PHEP(5,KHEP)**2) -
668      &    (PHEP(4,IHEP)**2-PHEP(5,IHEP)**2) -
669      &    (PHEP(3,JHEP)**2) )*HALF/PHEP(3,JHEP)
670         PHEP(3,KHEP)=-PHEP(3,IHEP)-PHEP(3,JHEP)
671         PHEP(2,JHEP)=0
672         PHEP(2,IHEP)=0
673         PHEP(2,KHEP)=0
674         PHEP(1,JHEP)=0
675         PHEP(1,IHEP)=HWUSQR(PHEP(4,IHEP)**2-
676      &    PHEP(3,IHEP)**2-PHEP(5,IHEP)**2)
677         PHEP(1,KHEP)=-PHEP(1,IHEP)
678 C---ORIENT IN CMF, THEN BOOST TO LAB
679         CALL HWUROB(R,PHEP(1,IHEP),PHEP(1,IHEP))
680         CALL HWUROB(R,PHEP(1,JHEP),PHEP(1,JHEP))
681         CALL HWUROB(R,PHEP(1,KHEP),PHEP(1,KHEP))
682         CALL HWULOB(PHEP(1,ICMF),PHEP(1,IHEP),PHEP(1,IHEP))
683         CALL HWULOB(PHEP(1,ICMF),PHEP(1,JHEP),PHEP(1,JHEP))
684         CALL HWULOB(PHEP(1,ICMF),PHEP(1,KHEP),PHEP(1,KHEP))
685 C---CALCULATE PRODUCTION VERTICES
686         CALL HWVZRO(4,VHEP(1,JHEP))
687         CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,KHEP),PVRT)
688         CALL HWUDKL(ID,PVRT,VHEP(1,KHEP))
689         CALL HWVEQU(4,VHEP(1,KHEP),VHEP(1,IHEP))
690 C---REORDER ENTRIES: IHEP=QUARK, JHEP=ANTI-QUARK, KHEP=EMITTED
691         IF (IHEP.EQ.LHEP) THEN
692           IHEP=JHEP
693           JHEP=LHEP
694         ENDIF
695 C---STATUS, ID AND POINTERS
696         ISTHEP(JHEP)=114
697         IDHW(JHEP)=IDHW(KHEP)
698         IDHEP(JHEP)=IDHEP(KHEP)
699         IDHW(KHEP)=ID3
700         IDHEP(KHEP)=IDPDG(ID3)
701         JDAHEP(2,ICMF)=JHEP
702         JMOHEP(1,JHEP)=ICMF
703         JDAHEP(1,JHEP)=0
704 C---COLOUR CONNECTIONS AND GLUON POLARIZATION
705         JMOHEP(2,JHEP)=IHEP
706         JDAHEP(2,IHEP)=JHEP
707         IF (ID3.EQ.13) THEN
708           JMOHEP(2,IHEP)=KHEP
709           JMOHEP(2,KHEP)=JHEP
710           JDAHEP(2,JHEP)=KHEP
711           JDAHEP(2,KHEP)=IHEP
712           GPOLN=((1-X(1))**2+(1-X(2))**2)/(4*(1-X(3)))
713           GPOLN=1/(1+GPOLN)
714         ELSE
715           JMOHEP(2,IHEP)=JHEP
716           JMOHEP(2,KHEP)=KHEP
717           JDAHEP(2,JHEP)=IHEP
718           JDAHEP(2,KHEP)=KHEP
719         ENDIF
720         IEVT=NEVHEP+NWGTS
721         GOTO 5
722       ELSEIF (IOPT.EQ.2) THEN
723 C---MAKE THREE-JET EVENTS FROM THE `DEAD-ZONE' LOOK LIKE TWO-JET EVENTS
724         IF (EMIT.EQ.0.OR.IEVT.NE.NEVHEP+NWGTS) THEN
725           RETURN
726         ELSEIF (EMIT.EQ.1) THEN
727           IHEP=JDAHEP(1,JDAHEP(1,ICMF)+1)
728           JHEP=JDAHEP(1,JDAHEP(1,ICMF))
729         ELSE
730           IHEP=JDAHEP(1,JDAHEP(2,ICMF))
731           JHEP=JDAHEP(1,JDAHEP(1,ICMF)+1)
732           JDAHEP(1,JDAHEP(2,ICMF))=JHEP
733           IDHW(JHEP)=IDHW(IHEP)
734           IF (ISTHEP(IHEP+1).EQ.100 .AND. ISTHEP(JHEP+1).EQ.100)
735      &      CALL HWVEQU(5,PHEP(1,IHEP+1),PHEP(1,JHEP+1))
736         ENDIF
737         JMOHEP(2,JDAHEP(1,ICMF))=JDAHEP(2,ICMF)
738         JDAHEP(2,JDAHEP(1,ICMF))=JDAHEP(2,ICMF)
739         JMOHEP(2,JDAHEP(2,ICMF))=JDAHEP(1,ICMF)
740         JDAHEP(2,JDAHEP(2,ICMF))=JDAHEP(1,ICMF)
741         CALL HWVEQU(5,P1,PHEP(1,JDAHEP(1,ICMF)))
742         CALL HWVEQU(5,P2,PHEP(1,JDAHEP(2,ICMF)))
743         CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,JHEP),PHEP(1,JHEP))
744         CALL HWUMAS(PHEP(1,JHEP))
745         JDAHEP(2,JHEP)=JDAHEP(2,IHEP)
746         IEDT(1)=JDAHEP(1,ICMF)+1
747         IEDT(2)=IHEP
748         IEDT(3)=IHEP+1
749         NDEL=3
750         IF (ISTHEP(IHEP+1).NE.100) NDEL=2
751         CALL HWUEDT(NDEL,IEDT)
752         DO 410 I=1,2
753           IHEP=JDAHEP(1,JDAHEP(I,ICMF))
754           JMOHEP(1,IHEP)=JDAHEP(I,ICMF)
755           IF (ISTHEP(IHEP+1).EQ.100) THEN
756             JMOHEP(1,IHEP+1)=JMOHEP(1,IHEP)
757             JMOHEP(2,IHEP+1)=JMOHEP(2,JMOHEP(1,IHEP))
758           ENDIF
759           DO 400 JHEP=JDAHEP(1,IHEP),JDAHEP(2,IHEP)
760             JMOHEP(1,JHEP)=IHEP
761  400      CONTINUE
762           CALL HWVZRO(4,VHEP(1,JDAHEP(I,ICMF)))
763           CALL HWVZRO(4,VHEP(1,IHEP))
764           IF (ISTHEP(IHEP+1).EQ.100) CALL HWVZRO(4,VHEP(1,IHEP+1))
765  410    CONTINUE
766         EMIT=0
767         IEVT=0
768       ELSE
769         CALL HWWARN('HWBDED',500)
770       ENDIF
771       END
772 CDECK  ID>, HWBDIS.
773 *CMZ :-        -17/05/94  09.33.08  by  Mike Seymour
774 *-- Author :    Mike Seymour
775 C-----------------------------------------------------------------------
776       SUBROUTINE HWBDIS(IOPT)
777 C-----------------------------------------------------------------------
778 C     FILL MISSING AREA OF DIS PHASE-SPACE WITH 2+1-JET EVENTS
779 C     IF (IOPT.EQ.1) SET UP EVENT RECORD
780 C     IF (IOPT.EQ.2) CLEAN UP EVENT RECORD AFTER SHOWERING
781 C-----------------------------------------------------------------------
782       INCLUDE 'HERWIG65.INC'
783       DOUBLE PRECISION HWRGEN,HWBVMC,HWUALF,HWULDO,P1(5),P2(5),P3(5),
784      & PCMF(5),L(5),R(3,3),Q,XBJ,RN,XPMIN,XPMAX,XP,ZPMIN,ZPMAX,ZP,FAC,
785      & X1,X2,XTSQ,XT,PTSQ,SIN1,SIN2,W1,W2,CFAC,PDFOLD(13),PDFNEW(13),
786      & PHI,SCALE,Q1(5),Q2(5),DIR1,DIR2,DIR,PM(5),POLD,PNEW,COMINT,
787      & BGFINT,COMWGT,C1,C2,CM,B1,B2,BM,PVRT(4)
788       INTEGER IOPT,EMIT,ICMF,IHEP,JHEP,IIN,IOUT,ILEP,IHAD,ID,IDNEW,
789      & IEDT(3),NDEL,NTRY,ITEMP
790       LOGICAL BGF
791       EXTERNAL HWRGEN,HWBVMC,HWUALF,HWULDO
792       SAVE BGF,IIN,IOUT,ICMF,ID,Q1,Q2,XP,XBJ
793       SAVE EMIT,COMINT,BGFINT,COMWGT,C1,C2,CM,B1,B2,BM
794       DATA EMIT,COMINT,BGFINT,COMWGT/0D0,3.9827D0,1.2462D0,0.3D0/
795       DATA C1,C2,CM,B1,B2,BM/0.56D0,0.20D0,10D0,0.667D0,0.167D0,3D0/
796       IF (IERROR.NE.0) RETURN
797       IF (IOPT.EQ.1) THEN
798 C---FIND AN UNTREATED CMF
799         IF (EMIT.EQ.NEVHEP+NWGTS) RETURN
800         ICMF=0
801         DO 10 IHEP=1,NHEP
802  10       IF (ICMF.EQ.0 .AND. ISTHEP(IHEP).EQ.110 .AND.
803      &    JDAHEP(2,IHEP).EQ.JDAHEP(1,IHEP)+1) ICMF=IHEP
804         IF (ICMF.EQ.0) RETURN
805         IIN=JMOHEP(2,ICMF)
806         IOUT=JDAHEP(2,ICMF)
807         ILEP=JMOHEP(1,ICMF)
808         CALL HWVEQU(5,PHEP(1,IIN),P1)
809         CALL HWVEQU(5,PHEP(1,IOUT),P2)
810         CALL HWVEQU(5,PHEP(1,ILEP),L)
811         IHAD=2
812         IF (JDAHEP(1,IHAD).NE.0) IHAD=JDAHEP(1,IHAD)
813         ID=IDHW(IIN)
814 C---STORE OLD MOMENTA
815         CALL HWVEQU(5,P1,Q1)
816         CALL HWVEQU(5,P2,Q2)
817 C---BOOST AND ROTATE THE MOMENTA TO THE BREIT FRAME
818         CALL HWVDIF(4,P2,P1,PCMF)
819         CALL HWUMAS(PCMF)
820         CALL HWVEQU(5,PHEP(1,IHAD),PM)
821         Q=-PCMF(5)
822         XBJ=HALF*Q**2/HWULDO(PM,PCMF)
823         CALL HWVSCA(4,HALF/XBJ,PCMF,PCMF)
824         CALL HWVSUM(4,PM,PCMF,PCMF)
825         CALL HWUMAS(PCMF)
826         CALL HWULOF(PCMF,L,L)
827         CALL HWULOF(PCMF,PM,PM)
828         CALL HWUROT(PM,ONE,ZERO,R)
829         CALL HWUROF(R,L,L)
830         PHI=ATAN2(L(2),L(1))
831         CALL HWUROT(PM,COS(PHI),SIN(PHI),R)
832 C---CHOOSE THE HADRONIC-PLANE CONFIGURATION, XP,ZP
833         IF (HWRGEN(0).LT.COMWGT) THEN
834 C-----CONSIDER GENERATING A QCD COMPTON EVENT
835           BGF=.FALSE.
836           P3(5)=RMASS(13)
837  100      RN=HWRGEN(1)
838           IF (RN.LT.C1) THEN
839             ZP=HWRGEN(2)
840             XPMAX=MIN(ZP,1-ZP)
841             XP=HWRGEN(3)*XPMAX
842             FAC=1/C1*2*XPMAX/((1-XP)*(1-ZP))*
843      $           (1+(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
844             IF (HWRGEN(4).LT.HALF) THEN
845               ZPMAX=ZP
846               ZP=XP
847               XP=ZPMAX
848             ENDIF
849           ELSEIF (RN.LT.C1+C2) THEN
850             XPMAX=0.83
851             XP=XPMAX*HWRGEN(2)
852             ZPMIN=MAX(XP,1-XP)
853             ZPMAX=1-2./3.*XP*(1+DREAL( DCMPLX(10-45*XP+18*XP**2,3*SQRT(
854      $         3*(9+66*XP-93*XP**2+12*XP**3-8*XP**4+24*XP**5-8*XP**6)))
855      $         **(1./3.) * DCMPLX(0.5D0,0.86602540378444D0) ))
856             ZP=1-((1-ZPMIN)/(1-ZPMAX))**HWRGEN(4)*(1-ZPMAX)
857             FAC=1/C2*XPMAX*LOG((1-ZPMIN)/(1-ZPMAX))/(1-XP)*
858      $           (1+(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
859           ELSE
860             ZPMAX=0.85
861             ZP=ZPMAX*HWRGEN(2)
862             XPMIN=MAX(ZP,1-ZP)
863             XPMAX=(1+4*ZP*(1-ZP))/(1+6*ZP*(1-ZP))
864             XP=1-((1-XPMIN)/(1-XPMAX))**HWRGEN(4)*(1-XPMAX)
865             FAC=1/(1-C1-C2)*ZPMAX*LOG((1-XPMIN)/(1-XPMAX))/(1-ZP)*
866      $           (1+(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
867           ENDIF
868           XPMAX=(1+4*ZP*(1-ZP))/(1+6*ZP*(1-ZP))
869           ZPMAX=1-2./3.*XP*(1+DREAL( DCMPLX(10-45*XP+18*XP**2,3*SQRT(
870      $         3*(9+66*XP-93*XP**2+12*XP**3-8*XP**4+24*XP**5-8*XP**6)))
871      $         **(1./3.) * DCMPLX(0.5D0,0.86602540378444D0) ))
872           IF (XP.GT.XPMAX.OR.ZP.GT.ZPMAX.OR.CM*HWRGEN(4).GT.FAC)
873      $         GOTO 100
874         ELSE
875 C-----CONSIDER GENERATING A BGF EVENT
876           BGF=.TRUE.
877           P3(5)=P1(5)
878           P1(5)=RMASS(13)
879  110      RN=HWRGEN(1)
880           IF (RN.LT.B1) THEN
881             ZP=HWRGEN(2)
882             XPMAX=MIN(ZP,1-ZP)
883             XP=HWRGEN(3)*XPMAX
884             FAC=1/B1*2*XPMAX/(1-ZP)*
885      $           ((  XP+ZP-2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP
886      $           +(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
887             IF (HWRGEN(4).LT.HALF) XP=1-XP
888           ELSEIF (RN.LT.B1+B2) THEN
889             XPMAX=0.83
890             XP=XPMAX*HWRGEN(2)
891             ZPMIN=MAX(XP,1-XP)
892             ZPMAX=1-2./3.*XP*(1+DREAL( DCMPLX(10-45*XP+18*XP**2,3*SQRT(
893      $         3*(9+66*XP-93*XP**2+12*XP**3-8*XP**4+24*XP**5-8*XP**6)))
894      $         **(1./3.) * DCMPLX(0.5D0,0.86602540378444D0) ))
895             ZP=1-((1-ZPMIN)/(1-ZPMAX))**HWRGEN(4)*(1-ZPMAX)
896             FAC=1/B2*XPMAX*LOG((1-ZPMIN)/(1-ZPMAX))*
897      $           ((  XP+ZP-2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP
898      $           +(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
899           ELSE
900             XPMAX=0.83
901             XP=XPMAX*HWRGEN(2)
902             ZPMAX=MIN(XP,1-XP)
903             ZPMIN=2./3.*XP*(1+DREAL( DCMPLX(10-45*XP+18*XP**2,3*SQRT(
904      $         3*(9+66*XP-93*XP**2+12*XP**3-8*XP**4+24*XP**5-8*XP**6)))
905      $         **(1./3.) * DCMPLX(0.5D0,0.86602540378444D0) ))
906             ZP=(ZPMAX-ZPMIN)*HWRGEN(4)+ZPMIN
907             FAC=1/(1-B1-B2)*XPMAX*(ZPMAX-ZPMIN)/(1-ZP)*
908      $           ((  XP+ZP-2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP
909      $           +(1-XP-ZP+2*XP*ZP)**2+2*(1-XP)*(1-ZP)*XP*ZP)
910           ENDIF
911           ZPMAX=1-2./3.*XP*(1+DREAL( DCMPLX(10-45*XP+18*XP**2,3*SQRT(
912      $         3*(9+66*XP-93*XP**2+12*XP**3-8*XP**4+24*XP**5-8*XP**6)))
913      $         **(1./3.) * DCMPLX(0.5D0,0.86602540378444D0) ))
914           IF (ZP.GT.ZPMAX.OR.ZP.LT.ONE-ZPMAX.OR.BM*HWRGEN(4).GT.FAC)
915      $         GOTO 110
916         ENDIF
917 C---CALCULATE THE ADDITIONAL FACTORS IN THE WEIGHT
918         IF (BGF) THEN
919           IDNEW=13
920           CFAC=1./2
921           FAC=BGFINT/(1-COMWGT)
922         ELSE
923           IDNEW=ID
924           CFAC=4./3
925           FAC=COMINT/COMWGT
926         ENDIF
927         SCALE=Q*SQRT((1-XP)*(1-ZP)*ZP/XP+1)
928         ITEMP=ISTAT
929         ISTAT=7
930         CALL HWSFUN(XBJ,Q,IDHW(IHAD),NSTRU,PDFOLD,2)
931         ISTAT=ITEMP
932         IF (PDFOLD(ID).LE.ZERO) THEN
933           CALL HWWARN('HWBDIS',100)
934           GOTO 999
935         ENDIF
936         IF (XP.GT.XBJ) THEN
937           CALL HWSFUN(XBJ/XP,SCALE,IDHW(IHAD),NSTRU,PDFNEW,2)
938           FAC=CFAC/(2*PIFAC) * HWUALF(1,SCALE) * FAC *
939      $         PDFNEW(IDNEW)/PDFOLD(ID)
940         ELSE
941           FAC=0
942         ENDIF
943 C---FOR PHOTON BEAMS, INCLUDE DIRECT PHOTON COUPLING
944         IF (IDHW(IHAD).EQ.59) THEN
945           ZPMIN=2./3.*XBJ*(1+DREAL( DCMPLX(10-45*XBJ+18*XBJ**2,3*SQRT(
946      $         3*(9+66*XBJ-93*XBJ**2+12*XBJ**3-8*XBJ**4+24*XBJ**5
947      $         -8*XBJ**6)))**(1./3.)*DCMPLX(0.5D0,0.86602540378444D0) ))
948           ZPMAX=1-ZPMIN
949           DIR1=(XBJ**2+(1-XBJ)**2)*(LOG(ZPMAX/ZPMIN)-(ZPMAX-ZPMIN))
950           DIR2=4*XBJ*(1-XBJ)*(ZPMAX-ZPMIN)
951           DIR=QFCH(MOD(ID-1,6)+1)**2*ALPHEM/(2*PIFAC*PDFOLD(ID))*XBJ
952      $         *(DIR1+DIR2)
953         ELSE
954           DIR=0
955         ENDIF
956 C---DECIDE WHETHER TO MAKE AN EVENT HERE
957         IF (HWRGEN(4).GT.FAC+DIR) RETURN
958 C---FOR DIRECT COUPLING, CHOOSE ZP VALUE
959         IF ((FAC+DIR)*HWRGEN(8).GT.FAC) THEN
960           IF ((DIR1+DIR2)*HWRGEN(9).LT.DIR1) THEN
961             NTRY=0
962  120        NTRY=NTRY+2
963             ZP=1-(ZPMAX/ZPMIN)**HWRGEN(NTRY+1)*ZPMIN
964             IF ((ZPMIN**2+(1-ZPMIN)**2)*HWRGEN(NTRY).GT.ZP**2+(1-ZP)**2)
965      $           GOTO 120
966           ELSE
967             ZP=SQRT((ZPMAX-ZPMIN)*HWRGEN(10)+ZPMIN**2)
968           ENDIF
969           XP=XBJ
970           BGF=.TRUE.
971           P3(5)=P2(5)
972           P1(5)=0
973         ENDIF
974         X1=1-   ZP /XP
975         X2=1-(1-ZP)/XP
976         XTSQ=4*(1-XP)*(1-ZP)*ZP/XP
977         XT=SQRT(XTSQ)
978         SIN1=XT/SQRT(X1**2+XTSQ)
979         SIN2=XT/SQRT(X2**2+XTSQ)
980 C---CHOOSE THE AZIMUTH BETWEEN THE TWO PLANES
981         IF (BGF) THEN
982           W1=XP**2*(X1**2+1.5*XTSQ)
983         ELSE
984           W1=1
985         ENDIF
986         W2=XP**2*(X2**2+1.5*XTSQ)
987         IF (HWRGEN(5)*(W1+W2).GT.W2) THEN
988           IF (BGF) THEN
989 C-----WEIGHTED BY (1+SIN1*COS(PHI))**2
990  200        PHI=(2*HWRGEN(6)-1)*PIFAC
991             IF (HWRGEN(7)*(1+SIN1)**2.GT.(1+SIN1*COS(PHI))**2) GOTO 200
992           ELSE
993 C-----UNIFORMLY
994             PHI=(2*HWRGEN(6)-1)*PIFAC
995           ENDIF
996         ELSE
997 C-----WEIGHTED BY (1-SIN2*COS(PHI))**2
998  210      PHI=(2*HWRGEN(6)-1)*PIFAC
999           IF (HWRGEN(7)*(1+SIN2)**2.GT.(1-SIN2*COS(PHI))**2) GOTO 210
1000         ENDIF
1001 C---RECONSTRUCT MOMENTA AND BOOST BACK TO LAB
1002         P1(1)=0
1003         P1(2)=0
1004         P1(3)=HALF*Q/XP
1005         P1(4)=SQRT(P1(3)**2+P1(5)**2)
1006         PTSQ=((ZP*Q*(P1(4)+P1(3)-Q)-P2(5)**2)*(P1(4)-P1(3)+(1-ZP)*Q)
1007      $       -P3(5)**2*ZP*Q)/(P1(4)-P1(3)+Q)
1008 C---CHECK INFRARED CUTOFF FOR THIS PARTON TYPE
1009         IF (PTSQ.LT.MAX(HWBVMC(ID),HWBVMC(IDHW(IOUT)))**2) RETURN
1010         P2(1)=SQRT(PTSQ)*COS(PHI)
1011         P2(2)=SQRT(PTSQ)*SIN(PHI)
1012         P2(3)=-0.5*(ZP*Q-(PTSQ+P2(5)**2)/(ZP*Q))
1013         P2(4)= 0.5*(ZP*Q+(PTSQ+P2(5)**2)/(ZP*Q))
1014         P3(1)=P1(1)-P2(1)
1015         P3(2)=P1(2)-P2(2)
1016         P3(3)=P1(3)-P2(3)-Q
1017         P3(4)=P1(4)-P2(4)
1018         CALL HWUROB(R,P1,P1)
1019         CALL HWUROB(R,P2,P2)
1020         CALL HWUROB(R,P3,P3)
1021         CALL HWULOB(PCMF,P1,P1)
1022         CALL HWULOB(PCMF,P2,P2)
1023         CALL HWULOB(PCMF,P3,P3)
1024 C---SPECIAL CASE FOR DIRECT PHOTON - COPY THE EXACT BEAM MOMENTUM
1025 C---SHARE THE MISMATCH EQUALLY BETWEEN THE OUTGOING PARTONS
1026 C---AND PUT THEM BACK ON SHELL
1027         IF (XP.EQ.XBJ) THEN
1028           CALL HWVDIF(4,PHEP(1,IHAD),P1,PM)
1029           CALL HWVSCA(4,HALF,PM,PM)
1030           CALL HWVSUM(4,PM,P2,P2)
1031           CALL HWVSUM(4,PM,P3,P3)
1032           CALL HWUMAS(P2)
1033           CALL HWUMAS(P3)
1034           CALL HWVEQU(5,PHEP(1,IHAD),P1)
1035           CALL HWVSUM(4,P2,P3,PCMF)
1036           CALL HWUMAS(PCMF)
1037           POLD=HWULDO(P2,PCMF)**2/PCMF(5)**2-SIGN(P2(5)**2,P2(5))
1038           PNEW=PCMF(5)**2/4-RMASS(ID)**2
1039           IF (PCMF(5).LE.ZERO.OR.POLD.LE.ZERO.OR.PNEW.LE.ZERO) RETURN
1040           CALL HWVSCA(4,SQRT(PNEW/POLD),P2,P2)
1041           CALL HWVSCA(4,HALF-HWULDO(P2,PCMF)/PCMF(5)**2,PCMF,PM)
1042           CALL HWVSUM(4,PM,P2,P2)
1043           CALL HWUMAS(P2)
1044           CALL HWVDIF(4,PCMF,P2,P3)
1045           CALL HWUMAS(P3)
1046         ENDIF
1047         NHEP=NHEP+1
1048         CALL HWVEQU(5,P1,PHEP(1,IIN))
1049         IF (BGF.AND.ID.GT.6.OR..NOT.BGF.AND.ID.LT.7) THEN
1050           CALL HWVEQU(5,P2,PHEP(1,IOUT))
1051           CALL HWVEQU(5,P3,PHEP(1,NHEP))
1052         ELSE
1053           CALL HWVEQU(5,P3,PHEP(1,IOUT))
1054           CALL HWVEQU(5,P2,PHEP(1,NHEP))
1055         ENDIF
1056         CALL HWVSUM(4,PHEP(1,ILEP),PHEP(1,IIN),PHEP(1,ICMF))
1057         CALL HWUMAS(PHEP(1,ICMF))
1058 C Decide which quark radiated and assign production vertices
1059         IF (BGF) THEN
1060 C Boson-Gluon fusion case
1061           IF (1-ZP.LT.HWRGEN(0)) THEN
1062 C Gluon splitting to quark
1063             CALL HWVZRO(4,VHEP(1,NHEP-1))
1064             CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
1065             CALL HWUDKL(ID,PVRT,VHEP(1,NHEP))
1066             CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
1067           ELSE
1068 C Gluon splitting to antiquark
1069             CALL HWVZRO(4,VHEP(1,NHEP))
1070             CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP-1),PVRT)
1071             CALL HWUDKL(ID,PVRT,VHEP(1,NHEP-1))
1072             CALL HWVEQU(4,VHEP(1,NHEP-1),VHEP(1,NHEP-4))
1073           ENDIF
1074         ELSE
1075 C QCD Compton case
1076           IF (1.LT.HWRGEN(0)*(1+(1-XP-ZP)**2+6*XP*(1-XP)*ZP*(1-ZP)))THEN
1077 C Incoming quark radiated the gluon
1078             CALL HWVZRO(4,VHEP(1,NHEP-1))
1079             CALL HWVDIF(4,PHEP(1,NHEP-4),PHEP(1,NHEP),PVRT)
1080             CALL HWUDKL(ID,PVRT,VHEP(1,NHEP))
1081             CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-4))
1082           ELSE
1083 C Outgoing quark radiated the gluon
1084             CALL HWVZRO(4,VHEP(1,NHEP-4))
1085             CALL HWVSUM(4,PHEP(1,NHEP-1),PHEP(1,NHEP),PVRT)
1086             CALL HWUDKL(ID,PVRT,VHEP(1,NHEP))
1087             CALL HWVEQU(4,VHEP(1,NHEP),VHEP(1,NHEP-1))
1088           ENDIF
1089         ENDIF
1090 C---STATUS, ID AND POINTERS
1091         ISTHEP(NHEP)=114
1092         IF (BGF) THEN
1093           IF (XP.EQ.XBJ) THEN
1094             IDHW(IIN)=59
1095             IDHEP(IIN)=IDPDG(59)
1096           ELSE
1097             IDHW(IIN)=13
1098             IDHEP(IIN)=IDPDG(13)
1099           ENDIF
1100           IF (ID.LT.7) THEN
1101             IDHW(NHEP)=IDHW(IOUT)
1102             IDHEP(NHEP)=IDHEP(IOUT)
1103             IDHW(IOUT)=MOD(ID,6)+6
1104             IDHEP(IOUT)=IDPDG(IDHW(IOUT))
1105           ELSE
1106             IDHW(NHEP)=MOD(ID,6)
1107             IDHEP(NHEP)=IDPDG(IDHW(NHEP))
1108           ENDIF
1109         ELSEIF (ID.LT.7) THEN
1110           IDHW(NHEP)=13
1111           IDHEP(NHEP)=IDPDG(13)
1112         ELSE
1113           IDHW(NHEP)=IDHW(IOUT)
1114           IDHEP(NHEP)=IDHEP(IOUT)
1115           IDHW(IOUT)=13
1116           IDHEP(IOUT)=IDPDG(13)
1117         ENDIF
1118         JDAHEP(2,ICMF)=NHEP
1119         JMOHEP(1,NHEP)=ICMF
1120 C---COLOUR CONNECTIONS
1121         IF (XP.EQ.XBJ) THEN
1122           JMOHEP(2,IIN)=IIN
1123           JDAHEP(2,IIN)=IIN
1124           JMOHEP(2,IOUT)=NHEP
1125           JDAHEP(2,IOUT)=NHEP
1126           JMOHEP(2,NHEP)=IOUT
1127           JDAHEP(2,NHEP)=IOUT
1128         ELSE
1129           JDAHEP(2,IIN)=NHEP
1130           JDAHEP(2,NHEP)=IOUT
1131           JMOHEP(2,IOUT)=NHEP
1132           JMOHEP(2,NHEP)=IIN
1133         ENDIF
1134 C---FACTORISATION SCALE
1135         EMSCA=SCALE
1136         EMIT=NEVHEP+NWGTS
1137       ELSEIF (IOPT.EQ.2) THEN
1138 C---MAKE TWO-JET EVENTS LOOK LIKE ONE-JET EVENTS
1139         IF (EMIT.NE.NEVHEP+NWGTS .OR. XP.EQ.XBJ) RETURN
1140         IF (.NOT.BGF) THEN
1141           CALL HWVEQU(5,Q1,PHEP(1,IIN))
1142           CALL HWVEQU(5,Q2,PHEP(1,IOUT))
1143           JMOHEP(2,IIN)=IOUT
1144           JDAHEP(2,IIN)=IOUT
1145           JMOHEP(2,IOUT)=IIN
1146           JDAHEP(2,IOUT)=IIN
1147           JDAHEP(2,ICMF)=IOUT
1148           IHEP=JDAHEP(1,IOUT)
1149           JHEP=JDAHEP(1,IOUT+1)
1150           CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,JHEP),PHEP(1,IHEP))
1151           CALL HWUMAS(PHEP(1,IHEP))
1152           JDAHEP(2,IHEP)=JDAHEP(2,JHEP)
1153           IEDT(1)=IOUT+1
1154           IEDT(2)=JHEP
1155           IEDT(3)=JHEP+1
1156           NDEL=3
1157           IF (ISTHEP(JHEP+1).NE.100) NDEL=2
1158           IHEP=JDAHEP(1,IOUT)
1159           JMOHEP(1,IHEP)=IOUT
1160           IF (ISTHEP(IHEP+1).EQ.100) THEN
1161             JMOHEP(1,IHEP+1)=IOUT
1162             JMOHEP(2,IHEP+1)=IIN
1163           ENDIF
1164           DO 300 JHEP=JDAHEP(1,IHEP),JDAHEP(2,IHEP)
1165             JMOHEP(1,JHEP)=IHEP
1166  300      CONTINUE
1167           IF (IDHW(IOUT).EQ.13) IDHW(IOUT)=IDHW(IOUT+1)
1168           IDHEP(IOUT)=IDPDG(IDHW(IOUT))
1169           IDHW(IHEP)=IDHW(IOUT)
1170           CALL HWUEDT(NDEL,IEDT)
1171         ELSEIF (ID.LT.7) THEN
1172           CALL HWVEQU(5,Q1,PHEP(1,IIN))
1173           CALL HWVEQU(5,Q2,PHEP(1,IOUT+1))
1174           JMOHEP(2,IIN)=IOUT+1
1175           JDAHEP(2,IIN)=IOUT+1
1176           JMOHEP(2,IOUT+1)=IIN
1177           JDAHEP(2,IOUT+1)=IIN
1178           JDAHEP(2,ICMF)=IOUT+1
1179           IHEP=JDAHEP(1,IIN)
1180           JHEP=JDAHEP(1,IOUT)
1181           CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,JHEP),PHEP(1,IHEP))
1182           CALL HWUMAS(PHEP(1,IHEP))
1183           CALL HWVDIF(4,PHEP(1,ICMF),PHEP(1,JHEP),PHEP(1,ICMF))
1184           CALL HWUMAS(PHEP(1,ICMF))
1185           CALL HWUEMV(JDAHEP(2,JHEP)-JDAHEP(1,JHEP)+1,
1186      $         JDAHEP(1,JHEP),JDAHEP(2,IHEP))
1187           JHEP=JDAHEP(1,IOUT)
1188           JDAHEP(2,IHEP)=JDAHEP(2,JHEP)
1189           IEDT(1)=IOUT
1190           IEDT(2)=JHEP
1191           IEDT(3)=JHEP+1
1192           NDEL=3
1193           IF (ISTHEP(JHEP+1).NE.100) NDEL=2
1194           CALL HWUEDT(NDEL,IEDT)
1195           IHEP=JDAHEP(1,IIN)
1196           DO 400 JHEP=JDAHEP(1,IHEP),JDAHEP(2,IHEP)
1197             JMOHEP(1,JHEP)=IHEP
1198  400      CONTINUE
1199           IDHW(IIN)=ID
1200           IDHEP(IIN)=IDPDG(ID)
1201           IDHW(IHEP)=ID
1202         ELSE
1203           CALL HWVEQU(5,Q1,PHEP(1,IIN))
1204           CALL HWVEQU(5,Q2,PHEP(1,IOUT))
1205           JMOHEP(2,IIN)=IOUT
1206           JDAHEP(2,IIN)=IOUT
1207           JMOHEP(2,IOUT)=IIN
1208           JDAHEP(2,IOUT)=IIN
1209           JDAHEP(2,ICMF)=IOUT
1210           IHEP=JDAHEP(1,IIN)
1211           JHEP=JDAHEP(1,IOUT+1)
1212           CALL HWVDIF(4,PHEP(1,IHEP),PHEP(1,JHEP),PHEP(1,IHEP))
1213           CALL HWUMAS(PHEP(1,IHEP))
1214           CALL HWVDIF(4,PHEP(1,ICMF),PHEP(1,JHEP),PHEP(1,ICMF))
1215           CALL HWUMAS(PHEP(1,ICMF))
1216           CALL HWUEMV(JDAHEP(2,JHEP)-JDAHEP(1,JHEP)+1,
1217      $         JDAHEP(1,JHEP),JDAHEP(1,IHEP)-1)
1218           JHEP=JDAHEP(1,IOUT+1)
1219           JDAHEP(1,IHEP)=JDAHEP(1,JHEP)
1220           IEDT(1)=IOUT+1
1221           IEDT(2)=JHEP
1222           IEDT(3)=JHEP+1
1223           NDEL=3
1224           IF (ISTHEP(JHEP+1).NE.100.OR.JHEP.EQ.NHEP) NDEL=2
1225           CALL HWUEDT(NDEL,IEDT)
1226           IHEP=JDAHEP(1,IIN)
1227           DO 500 JHEP=JDAHEP(1,IHEP),JDAHEP(2,IHEP)
1228             JMOHEP(1,JHEP)=IHEP
1229  500      CONTINUE
1230           IDHW(IIN)=ID
1231           IDHEP(IIN)=IDPDG(ID)
1232           IDHW(IHEP)=ID
1233         ENDIF
1234         CALL HWVZRO(4,VHEP(1,IIN))
1235         CALL HWVZRO(4,VHEP(1,JDAHEP(1,IIN)))
1236         IF (ISTHEP(JDAHEP(1,IIN)+1).EQ.100)
1237      $       CALL HWVZRO(4,VHEP(1,JDAHEP(1,IIN)+1))
1238         CALL HWVZRO(4,VHEP(1,IOUT))
1239         CALL HWVZRO(4,VHEP(1,JDAHEP(1,IOUT)))
1240         IF (ISTHEP(JDAHEP(1,IOUT)+1).EQ.100)
1241      $       CALL HWVZRO(4,VHEP(1,JDAHEP(1,IOUT)+1))
1242         EMIT=0
1243       ELSE
1244         CALL HWWARN('HWBDIS',500)
1245       ENDIF
1246  999  RETURN
1247       END
1248 CDECK  ID>, HWBDYP.
1249 *CMZ :-        -26/10/99  17.46.56  by  Mike Seymour
1250 *-- Author :    Gennaro Corcella
1251 C-----------------------------------------------------------------------
1252       SUBROUTINE HWBDYP(IOPT)
1253 C     MATRIX ELEMENT CORRECTIONS TO DRELL-YAN PROCESSES
1254 C-----------------------------------------------------------------------
1255       INCLUDE 'HERWIG65.INC'
1256       DOUBLE PRECISION HWBVMC,HWRGEN,HWUALF,HWUSQR,PMODK,AZ,CZ,
1257      & T,U,S,EM,TMIN,TMAX,PMOD2,GLUFAC,SMIN,SMAX,SZ,TEST,
1258      & JAC,M(3),W1,W,PMOD3,SCAPR,CPHI,SPHI,SCALE,XI1,XI2,
1259      & PDFOLD1(13),PDFOLD2(13),PDFNEW1(13),PDFNEW2(13),ETA1,ETA2,Y,
1260      & COMWGT1,COMWGT2,WW,COS3,MODP,RN,BETA1,SIN3,R3(3,3),CTH,STH,M1,
1261      & M2,M3,GAMMA1,R5(3,3),CW,SW,R4(3,3),SCALE1,X1,X2,X3,MM,
1262      & PHAD1(5),PHAD2(5),P1(5),P2(5),P3(5),P4(5),PF(5),PV(5),PK(5),
1263      & PR(5),PNE(5),PE(5),PP1(5),PP2(5),PZ(5),PS(5),PD(5),P2N(5),
1264      & PBOS(5),PLAB(5),PTOT(5),P3N(5),SVNTN
1265       LOGICAL GLUIN,GP
1266       INTEGER EMIT,NOEMIT,IHEP,JHEP,KHEP,ICMF,IOPT,CHEP,
1267      & ID2,ID1,K,ID4,ID5,IDBOS,IHAD1,IHAD2,NTMP
1268       EXTERNAL HWBVMC,HWRGEN,HWUALF,HWUSQR
1269       SAVE PS,PF,ICMF,ID4,ID5
1270       SAVE EMIT,NTMP
1271       DATA EMIT,NTMP/2*0/
1272       IF (IOPT.EQ.1) THEN
1273         EMIT=0
1274         NTMP=0
1275 C-----CHOOSE WEIGHTS
1276         COMWGT1=0.1
1277         COMWGT2=0.55
1278 C---FIND AN UNTREATED CMF
1279         ICMF=0
1280         DO 10 IHEP=1,NHEP
1281  10     IF (ICMF.EQ.0 .AND. ISTHEP(IHEP).EQ.110.AND.
1282      &         JDAHEP(2,IHEP).EQ.JDAHEP(1,IHEP)+1) ICMF=IHEP
1283         IF (ICMF.EQ.0) RETURN
1284         EM=PHEP(5,ICMF)
1285 C-----SET THE VECTOR BOSON RAPIDITY
1286         Y=HALF*LOG((PHEP(4,ICMF)+PHEP(3,ICMF))/
1287      &       (PHEP(4,ICMF)-PHEP(3,ICMF)))
1288 C------SET PARTICLE IDENTIES
1289 c------ID1=QUARK, ID2=ANTIQUARK, IDBOS=VECTOR BOSON, ID4-5 BOSON DECAY
1290         IDBOS=IDHW(ICMF)
1291         ID1=IDHW(JMOHEP(1,ICMF))
1292         ID2=IDHW(JMOHEP(2,ICMF))
1293         ID4=IDHW(JDAHEP(1,ICMF))
1294         ID5=IDHW(JDAHEP(2,ICMF))
1295         M1=RMASS(ID1)
1296         M2=RMASS(ID2)
1297         M3=RMASS(13)
1298 C---STORE OLD MOMENTA
1299 C------VECTOR BOSON MOMENTUM
1300         CALL HWVEQU(5,PHEP(1,ICMF),PBOS)
1301 C----QUARK MOMENTUM
1302         CALL HWVEQU(5,PHEP(1,JMOHEP(1,ICMF)),P1)
1303 C------ANTIQUARK MOMENTUM
1304         CALL HWVEQU(5,PHEP(1,JMOHEP(2,ICMF)),P2)
1305 C-------VECTOR DECAY (LEPTON) PRODUCT MOMENTA
1306         CALL HWVEQU(5,PHEP(1,JDAHEP(1,ICMF)),P3)
1307         CALL HWVEQU(5,PHEP(1,JDAHEP(2,ICMF)),P4)
1308 C------LEPTON MOMENTA IN THE BOSON REST FRAME
1309         CALL HWULOF(PHEP(1,ICMF),P2,P2N)
1310         CALL HWULOF(PHEP(1,ICMF),P3,P3N)
1311 C------AZ=AZIMUTHAL ANGLE OF P3N
1312         AZ=ATAN2(P3N(2),P3N(1))
1313         CZ=COS(AZ)
1314         SZ=SIN(AZ)
1315 C------PHI=ANGLE BETWEEN P2N AND P3N
1316         SCAPR=P2N(1)*P3N(1)+P2N(2)*P3N(2)+P2N(3)*P3N(3)
1317         PMOD2=SQRT(P2N(1)**2+P2N(2)**2+P2N(3)**2)
1318         PMOD3=SQRT(P3N(1)**2+P3N(2)**2+P3N(3)**2)
1319         CPHI=SCAPR/(PMOD3*PMOD2)
1320         SPHI=SQRT(1-CPHI**2)
1321 C------HADRON MOMENTA
1322         IHAD1=1
1323         IHAD2=2
1324         IF (JDAHEP(1,IHAD1).NE.0) IHAD1=JDAHEP(1,IHAD1)
1325         IF (JDAHEP(1,IHAD2).NE.0) IHAD2=JDAHEP(1,IHAD2)
1326         CALL HWVEQU(5,PHEP(1,IHAD1),PHAD1)
1327         CALL HWVEQU(5,PHEP(1,IHAD2),PHAD2)
1328         CALL HWVSUM(4,PHAD1,PHAD2,PTOT)
1329         CALL HWUMAS(PTOT)
1330 C------ Q - QBAR ENERGY FRACTIONS (BORN PROCESS)
1331 c---minorimprovement---mhs---4/8/04---include mass effects correctly
1332         ETA1=(P1(4)+P1(3))/(PHAD1(4)+PHAD1(3))
1333         ETA2=(P2(4)-P2(3))/(PHAD2(4)-PHAD2(3))
1334 C------ PDFs FOR THE BORN PROCESS
1335         CALL HWSFUN(ETA1,EM,IDHW(IHAD1),NSTRU,PDFOLD1,1)
1336         CALL HWSFUN(ETA2,EM,IDHW(IHAD2),NSTRU,PDFOLD2,2)
1337 C-------CONSIDER Q(QBAR) IN THE INITIAL STATE
1338         RN=HWRGEN(9)
1339         IF (RN.LT.COMWGT1) THEN
1340 C-------NO GLUON IN THE INITIAL STATE
1341           GLUIN=.FALSE.
1342 C---CHOOSE S ACCORDING TO 1/S**2
1343           SVNTN=17
1344           SMIN=HALF*EM**2*(7-SQRT(SVNTN))
1345           SMAX=PTOT(5)**2
1346           IF (SMAX.LE.SMIN) RETURN
1347           S=SMIN*SMAX/(SMIN+HWRGEN(0)*(SMAX-SMIN))
1348           JAC=S**2*(1/SMIN-1/SMAX)
1349 C---CHOOSE T ACCORDING TO (S-EM**2)/(T*U)=1/T+1/U
1350           TMAX=-HALF*EM**2*(3-HWUSQR(1+8*EM**2/S))
1351           TMIN=EM**2-S-TMAX
1352           IF (TMAX.LE.TMIN) RETURN
1353           T=TMAX*(TMIN/TMAX)**HWRGEN(1)
1354           IF (HWRGEN(2).GT.HALF) T=EM**2-S-T
1355           U=EM**2-S-T
1356           JAC=JAC*2*T*U/(S-EM**2)*LOG(TMIN/TMAX)
1357           SCALE=SQRT(U*T/S)
1358           SCALE1=SQRT(U*T/S+EM**2)
1359           GLUFAC=0
1360           IF (SCALE1.GT.HWBVMC(13)) GLUFAC=HWUALF(1,SCALE1)/(2*PIFAC)
1361 C----Q-QBAR ENERGY FRACTIONS FOR Q QBAR-> VG
1362           XI1=(HALF/PHAD1(4))*EXP(Y)*SQRT(S*(S+T)/(S+U))
1363           XI2=S/(4*XI1*PHAD1(4)*PHAD2(4))
1364 c---minorimprovement---mhs---4/8/04---apply infrared cutoff for large x
1365           IF ((1-XI1)*SCALE.LT.HWBVMC(ID1)) RETURN
1366           IF ((1-XI2)*SCALE.LT.HWBVMC(ID2)) RETURN
1367 C-----PDFs WITH AN EMITTED GLUON
1368           CALL HWSFUN(XI1,SCALE,IDHW(IHAD1),NSTRU,PDFNEW1,1)
1369           CALL HWSFUN(XI2,SCALE,IDHW(IHAD2),NSTRU,PDFNEW2,2)
1370 C------CALCULATE WEIGHT
1371           W=JAC*((EM**2-T)**2+(EM**2-U)**2)/(S**2*T*U)
1372           W1=(GLUFAC/COMWGT1)*W*PDFNEW1(ID1)*PDFNEW2(ID2)/(PDFOLD1(ID1)*
1373      &         PDFOLD2(ID2))*(CFFAC*ETA1*ETA2/(XI1*XI2))
1374 C-------CHOOSE WHICH PARTON WILL EMIT
1375           EMIT=1
1376           IF (HWRGEN(6).LT.(EM**2-U)**2/((EM**2-U)**2+(EM**2-T)**2))
1377      &         EMIT=2
1378           NOEMIT=3-EMIT
1379         ELSE
1380 C--------GLUON IN THE INITIAL STATE
1381           GLUIN=.TRUE.
1382 C---CHOOSE S ACCORDING TO 1/S**2
1383           SMIN=EM**2
1384           SMAX=PTOT(5)**2
1385           IF (SMAX.LE.SMIN) RETURN
1386           S=SMIN*SMAX/(SMIN+HWRGEN(0)*(SMAX-SMIN))
1387           JAC=S**2*(1/SMIN-1/SMAX)
1388 C---CHOOSE T ACCORDING TO 1/T
1389           TMAX=-HALF*EM**2*(3-HWUSQR(1+8*EM**2/S))
1390           TMIN=EM**2-S
1391           IF (TMAX.LE.TMIN) RETURN
1392           T=TMAX*(TMIN/TMAX)**HWRGEN(1)
1393           JAC=JAC*T*LOG(TMAX/TMIN)
1394           U=EM**2-S-T
1395           SCALE=SQRT(U*T/S)
1396           SCALE1=SQRT(U*T/S+EM**2)
1397           GLUFAC=0
1398           IF (SCALE1.GT.HWBVMC(13)) GLUFAC=HWUALF(1,SCALE1)/(2*PIFAC)
1399 C--------INITIAL STATE GLUON COMING FROM HADRON 1
1400           IF (RN.LE.COMWGT2) THEN
1401             GP=.TRUE.
1402 C--------ENERGY FRACTIONS and PDFs
1403 c---bug fix---mhs---4/8/04---swap u and t in mtm frac definitions
1404             XI1=(HALF/PHAD1(4))*EXP(Y)*SQRT(S*(S+T)/(S+U))
1405             XI2=S/(4*XI1*PHAD1(4)*PHAD2(4))
1406 c---minorimprovement---mhs---4/8/04---apply infrared cutoff for large x
1407             IF ((1-XI1)*SCALE.LT.HWBVMC(13)) RETURN
1408             IF ((1-XI2)*SCALE.LT.HWBVMC(ID2)) RETURN
1409             CALL HWSFUN(XI1,SCALE,IDHW(IHAD1),NSTRU,PDFNEW1,1)
1410             CALL HWSFUN(XI2,SCALE,IDHW(IHAD2),NSTRU,PDFNEW2,2)
1411             WW=PDFNEW1(13)*PDFNEW2(ID2)/((COMWGT2-COMWGT1)*
1412      &           PDFOLD1(ID1)*PDFOLD2(ID2))
1413           ELSE
1414 C-------INITIAL STATE GLUON COMING FROM HADRON 2
1415             GP=.FALSE.
1416 C-------ENERGY FRACTIONS AND PDFs
1417 c---bug fix---mhs---4/8/04---swap u and t in mtm frac definitions
1418             XI1=(HALF/PHAD1(4))*EXP(Y)*SQRT(S*(S+U)/(S+T))
1419             XI2=S/(4*XI1*PHAD1(4)*PHAD2(4))
1420 c---minorimprovement---mhs---4/8/04---apply infrared cutoff for large x
1421             IF ((1-XI1)*SCALE.LT.HWBVMC(ID1)) RETURN
1422             IF ((1-XI2)*SCALE.LT.HWBVMC(13)) RETURN
1423             CALL HWSFUN(XI1,SCALE,IDHW(IHAD1),NSTRU,PDFNEW1,1)
1424             CALL HWSFUN(XI2,SCALE,IDHW(IHAD2),NSTRU,PDFNEW2,2)
1425             WW=PDFNEW1(ID1)*PDFNEW2(13)/((1-COMWGT2)*
1426      &           PDFOLD1(ID1)*PDFOLD2(ID2))
1427           ENDIF
1428           W=-HALF*JAC*((EM**2-T)**2+(EM**2-S)**2)/(S**3*T)
1429 C-------CHOOSE WHICH PARTON WILL EMIT
1430 c---bug fix---mhs---4/8/04---swap emitter and nonemitter
1431           EMIT=2
1432           IF (HWRGEN(10).LT.(EM**2-S)**2/((EM**2-S)**2+(EM**2-T)**2))
1433      &         EMIT=1
1434           NOEMIT=3-EMIT
1435 C-------FINAL WEIGHT FOR ALL THE CONSIDERED OPTIONS
1436           W1=GLUFAC*W*WW*ETA1*ETA2/(XI1*XI2)
1437         ENDIF
1438 C--------ADD ONE MORE GLUON
1439         IF (W1.GT.HWRGEN(4)) THEN
1440           NTMP=NEVHEP+NWGTS
1441         ELSE
1442           RETURN
1443         ENDIF
1444 C---------INCLUDE MASSES
1445         S=S+M1**2+M2**2+M3**2
1446         IF (.NOT.GLUIN) THEN
1447           TEST=((S+M1**2-M2**2)*(S+M3**2-EM**2)-2*S*(M1**2+M3**2-T))**2
1448      $         -((S-M1**2-M2**2)**2-4*M1**2*M2**2)*
1449      $         ((S-M3**2-EM**2)**2-4*M3**2*EM**2)
1450         ELSEIF (GP) THEN
1451           TEST=((S+M3**2-M2**2)*(S+M1**2-EM**2)-2*S*(M3**2+M1**2-T))**2
1452      $         -((S-M3**2-M2**2)**2-4*M3**2*M2**2)*
1453      $         ((S-M1**2-EM**2)**2-4*M1**2*EM**2)
1454         ELSE
1455           TEST=((S+M3**2-M1**2)*(S+M2**2-EM**2)-2*S*(M3**2+M2**2-T))**2
1456      $         -((S-M3**2-M1**2)**2-4*M3**2*M1**2)*
1457      $         ((S-M2**2-EM**2)**2-4*M2**2*EM**2)
1458         ENDIF
1459         IF (TEST.GE.0) THEN
1460           EMIT=0
1461           RETURN
1462         ENDIF
1463         M(1)=M1
1464         M(2)=M2
1465         M(3)=M3
1466 C----MOMENTA IN THE V-REST FRAME WITH NON EMITTER ALONG THE Z AXIS
1467 C----V=BOSON,K=GLUON,E=EMITTER,NE=NON-EMITTER
1468         PV(1)=0
1469         PV(2)=0
1470         PV(3)=0
1471         PV(4)=EM
1472         PV(5)=EM
1473         PNE(2)=0
1474         PNE(1)=0
1475         IF (.NOT.GLUIN) THEN
1476           PK(4)=(S-M(3)**2-EM**2)/(2*EM)
1477           PMODK=SQRT(PK(4)**2-M(3)**2)
1478           IF (EMIT.EQ.1) THEN
1479             MM=M(1)
1480             X1=T
1481             X2=U
1482             X3=-1
1483           ELSE
1484             MM=M(2)
1485             X1=U
1486             X2=T
1487             X3=+1
1488           ENDIF
1489           PNE(4)=(EM**2+MM**2-X1)/(2*EM)
1490           PNE(3)=X3*SQRT(PNE(4)**2-MM**2)
1491           COS3=HALF*(X2-MM**2-M(3)**2+2*PNE(4)*PK(4))/(PNE(3)*PMODK)
1492         ELSE
1493           PK(4)=(EM**2+M(3)**2-U)/(2*EM)
1494           PMODK=SQRT(PK(4)**2-M(3)**2)
1495           IF (EMIT.EQ.1) THEN
1496             IF (GP) THEN
1497               MM=M(1)
1498               X3=+1
1499             ELSE
1500               MM=M(2)
1501               X3=-1
1502             ENDIF
1503             PNE(4)=(S-MM**2-EM**2)/(2*EM)
1504             PNE(3)=X3*SQRT(PNE(4)**2-MM**2)
1505             COS3=HALF*(T-MM**2-M(3)**2+2*PNE(4)*PK(4))/(PNE(3)*PMODK)
1506           ELSE
1507             IF (GP) THEN
1508               MM=M(2)
1509               X3=-1
1510             ELSE
1511               MM=M(1)
1512               X3=+1
1513             ENDIF
1514             PNE(4)=(EM**2+MM**2-T)/(2*EM)
1515             PNE(3)=X3*SQRT(PNE(4)**2-MM**2)
1516             COS3=HALF*(MM**2+M(3)**2-S+2*PNE(4)*PK(4))/(PNE(3)*PMODK)
1517           ENDIF
1518         ENDIF
1519         CALL HWUMAS(PNE)
1520         SIN3=SQRT(1-COS3**2)
1521 C---------DEFINE A RANDOM ROTATION AROUND THE Z-AXIS
1522         CALL HWRAZM(PMODK*SIN3,PK(1),PK(2))
1523         PK(3)=PMODK*COS3
1524         CALL HWUMAS(PK)
1525         DO K=1,4
1526           IF (.NOT.GLUIN) THEN
1527             PE(K)=PV(K)+PK(K)-PNE(K)
1528           ELSE
1529             IF (EMIT.EQ.1) THEN
1530               PE(K)=PV(K)+PNE(K)-PK(K)
1531             ELSE
1532               PE(K)=PNE(K)+PK(K)-PV(K)
1533             ENDIF
1534           ENDIF
1535         ENDDO
1536         CALL HWUMAS(PE)
1537 c------LEPTON MOMENTA IN THE BOSON REST FRAME, WITH THE DIRECTION
1538 C------TAKEN FROM THE BORN PROCESS
1539         PS(5)=P3(5)
1540         PS(4)=(EM**2+P3(5)**2-P4(5)**2)/(2*EM)
1541         PS(3)=-SQRT(PS(4)**2-P3(5)**2)*CPHI
1542         PS(2)=SQRT(PS(4)**2-P3(5)**2)*SPHI*SZ
1543         PS(1)=SQRT(PS(4)**2-P3(5)**2)*SPHI*CZ
1544         PF(5)=P4(5)
1545         PF(4)=(EM**2+P4(5)**2-P3(5)**2)/(2*EM)
1546         PF(3)=-PS(3)
1547         PF(2)=-PS(2)
1548         PF(1)=-PS(1)
1549 C----FIND A STATIONARY VECTOR PLAB IN THE LAB FRAME
1550         IF (.NOT.GLUIN) THEN
1551           IF (EMIT.EQ.1) THEN
1552             CALL HWVEQU(5,PE,PP1)
1553             CALL HWVEQU(5,PNE,PP2)
1554           ELSE
1555             CALL HWVEQU(5,PNE,PP1)
1556             CALL HWVEQU(5,PE,PP2)
1557           ENDIF
1558         ELSE
1559           IF (GP) THEN
1560             CALL HWVEQU(5,PK,PP1)
1561             IF (EMIT.EQ.1) THEN
1562               CALL HWVEQU(5,PE,PP2)
1563             ELSE
1564               CALL HWVEQU(5,PNE,PP2)
1565             ENDIF
1566           ELSE
1567             CALL HWVEQU(5,PK,PP2)
1568             IF (EMIT.EQ.1) THEN
1569               CALL HWVEQU(5,PE,PP1)
1570             ELSE
1571               CALL HWVEQU(5,PNE,PP1)
1572             ENDIF
1573           ENDIF
1574         ENDIF
1575         CALL HWVSCA(4,1/XI1,PP1,PP1)
1576         CALL HWVSCA(4,1/XI2,PP2,PP2)
1577         CALL HWVSUM(4,PP1,PP2,PLAB)
1578         CALL HWUMAS(PLAB)
1579 C------BOOST TO PLAB REST FRAME
1580         CALL HWULOF(PLAB,PE,PE)
1581         CALL HWULOF(PLAB,PNE,PNE)
1582         CALL HWULOF(PLAB,PK,PK)
1583         CALL HWULOF(PLAB,PS,PS)
1584         CALL HWULOF(PLAB,PF,PF)
1585         CALL HWULOF(PLAB,PV,PV)
1586 C----PUT THE INITIAL PARTON BELONGING TO HADRON 1 ON THE Z-AXIS
1587         IF (.NOT.GLUIN) THEN
1588           IF (EMIT.EQ.1) THEN
1589             CALL HWVEQU(5,PE,PZ)
1590           ELSE
1591             CALL HWVEQU(5,PNE,PZ)
1592           ENDIF
1593         ELSE
1594           IF (GP) THEN
1595             CALL HWVEQU(5,PK,PZ)
1596           ELSE
1597             IF (EMIT.EQ.1) THEN
1598               CALL HWVEQU(5,PE,PZ)
1599             ELSE
1600               CALL HWVEQU(5,PNE,PZ)
1601             ENDIF
1602           ENDIF
1603         ENDIF
1604         MODP=SQRT(PZ(1)**2+PZ(2)**2)
1605         CTH=PZ(1)/MODP
1606         STH=PZ(2)/MODP
1607         CALL HWUROT(PZ,CTH,STH,R3)
1608 C-----ROTATE EVERYTHING BY R3
1609         CALL HWUROF(R3,PE,PE)
1610         CALL HWUROF(R3,PNE,PNE)
1611         CALL HWUROF(R3,PV,PV)
1612         CALL HWUROF(R3,PK,PK)
1613         CALL HWUROF(R3,PS,PS)
1614         CALL HWUROF(R3,PF,PF)
1615 C--REORDER ENTRIES:--IHEP=EMITTER,JHEP=NON-EMITTER,KHEP=EMITTED
1616         IF (.NOT.GLUIN) THEN
1617           IHEP=JMOHEP(EMIT,ICMF)
1618           JHEP=JMOHEP(NOEMIT,ICMF)
1619         ENDIF
1620         CHEP=ICMF
1621         IDHW(CHEP)=15
1622         IDHEP(CHEP)=IDPDG(15)
1623         ICMF=ICMF+1
1624         IDHW(ICMF)=IDBOS
1625         IDHEP(ICMF)=IDPDG(IDBOS)
1626 C-----NO GLUON IN THE INITIAL STATE: JUST ADD IT AFTER THE VECTOR BOSON
1627         IF (.NOT.GLUIN) THEN
1628           KHEP=ICMF+1
1629           ISTHEP(KHEP)=114
1630 C---STATUS OF EMITTER/NON EMITTER
1631           ISTHEP(IHEP)=110+EMIT
1632           ISTHEP(JHEP)=110+NOEMIT
1633         ELSE
1634 C-----GLUON COMING FROM THE 1ST HADRON
1635           IF (GP) THEN
1636             KHEP=CHEP-2
1637             ISTHEP(KHEP)=111
1638 C----EMIT=1
1639             IF (EMIT.EQ.1) THEN
1640               IHEP=KHEP+1
1641               ISTHEP(IHEP)=112
1642               JHEP=ICMF+1
1643               ISTHEP(JHEP)=114
1644               IDHW(IHEP)=ID2
1645               IF (ID1.LE.6) THEN
1646                 IDHW(JHEP)=ID1+6
1647               ELSE
1648                 IDHW(JHEP)=ID1-6
1649               ENDIF
1650             ELSE
1651 C-------EMIT=2
1652               JHEP=KHEP+1
1653               ISTHEP(JHEP)=112
1654               IDHW(JHEP)=ID2
1655               IHEP=ICMF+1
1656               ISTHEP(IHEP)=114
1657               IF (ID1.LE.6) THEN
1658                 IDHW(IHEP)=ID1+6
1659               ELSE
1660                 IDHW(IHEP)=ID1-6
1661               ENDIF
1662             ENDIF
1663           ENDIF
1664 C------GLUON COMING FROM THE HADRON 2
1665           IF (.NOT.GP) THEN
1666             KHEP=CHEP-1
1667             ISTHEP(KHEP)=112
1668 C-------EMIT=1
1669             IF (EMIT.EQ.1) THEN
1670               IHEP=KHEP-1
1671               ISTHEP(IHEP)=111
1672               IDHW(IHEP)=ID1
1673               JHEP=ICMF+1
1674               ISTHEP(JHEP)=114
1675               IF (ID2.LE.6) THEN
1676                 IDHW(JHEP)=ID2+6
1677               ELSE
1678                 IDHW(JHEP)=ID2-6
1679               ENDIF
1680             ELSE
1681 C-------EMIT=2
1682               JHEP=KHEP-1
1683               ISTHEP(JHEP)=111
1684               IDHW(JHEP)=ID1
1685               IHEP=ICMF+1
1686               ISTHEP(IHEP)=114
1687               IF (ID2.LE.6) THEN
1688                 IDHW(IHEP)=ID2+6
1689               ELSE
1690                 IDHW(IHEP)=ID2-6
1691               ENDIF
1692             ENDIF
1693           ENDIF
1694         ENDIF
1695         IDHEP(IHEP)=IDPDG(IDHW(IHEP))
1696         IDHEP(JHEP)=IDPDG(IDHW(JHEP))
1697         ISTHEP(ICMF)=113
1698         ISTHEP(CHEP)=110
1699         IDHW(KHEP)=13
1700         IDHEP(KHEP)=IDPDG(13)
1701 C---------DEFINE MOMENTA IN THE LAB FRAME
1702         CALL HWVEQU(5,PV,PHEP(1,ICMF))
1703         CALL HWVEQU(5,PK,PHEP(1,KHEP))
1704         CALL HWVEQU(5,PNE,PHEP(1,JHEP))
1705         CALL HWVEQU(5,PE,PHEP(1,IHEP))
1706         IF (.NOT.GLUIN) THEN
1707           CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,JHEP),PHEP(1,CHEP))
1708         ELSE
1709           IF (EMIT.EQ.1) THEN
1710             CALL HWVSUM(4,PHEP(1,IHEP),PHEP(1,KHEP),PHEP(1,CHEP))
1711           ELSE
1712             CALL HWVSUM(4,PHEP(1,KHEP),PHEP(1,JHEP),PHEP(1,CHEP))
1713           ENDIF
1714         ENDIF
1715         CALL HWUMAS(PHEP(1,CHEP))
1716         IF (.NOT.GLUIN) THEN
1717           JMOHEP(1,JHEP)=CHEP
1718           JMOHEP(1,IHEP)=CHEP
1719           JDAHEP(1,JHEP)=CHEP
1720           JDAHEP(1,IHEP)=CHEP
1721           JMOHEP(1,KHEP)=CHEP
1722           JDAHEP(1,KHEP)=0
1723           JMOHEP(1,ICMF)=CHEP
1724           JMOHEP(2,ICMF)=ICMF
1725           JDAHEP(1,ICMF)=0
1726           JDAHEP(2,ICMF)=ICMF
1727         ENDIF
1728         IF (GLUIN) THEN
1729           JMOHEP(2,ICMF)=ICMF
1730           JDAHEP(2,ICMF)=ICMF
1731           JMOHEP(1,KHEP)=CHEP
1732           JDAHEP(1,KHEP)=CHEP
1733           JMOHEP(1,IHEP)=CHEP
1734           JMOHEP(1,JHEP)=CHEP
1735           IF (EMIT.EQ.1) THEN
1736             JDAHEP(1,IHEP)=CHEP
1737             JDAHEP(1,JHEP)=0
1738           ELSE
1739             JDAHEP(1,JHEP)=CHEP
1740             JDAHEP(1,IHEP)=0
1741           ENDIF
1742         ENDIF
1743 C---COLOUR CONNECTIONS
1744         IF (.NOT.GLUIN) THEN
1745           IF (IDHW(IHEP).LT.IDHW(JHEP)) THEN
1746             JMOHEP(2,KHEP)=IHEP
1747             JDAHEP(2,KHEP)=JHEP
1748             JMOHEP(2,IHEP)=JHEP
1749             JDAHEP(2,IHEP)=KHEP
1750             JDAHEP(2,JHEP)=IHEP
1751             JMOHEP(2,JHEP)=KHEP
1752           ELSE
1753             JMOHEP(2,KHEP)=JHEP
1754             JDAHEP(2,KHEP)=IHEP
1755             JMOHEP(2,JHEP)=IHEP
1756             JDAHEP(2,JHEP)=KHEP
1757             JDAHEP(2,IHEP)=JHEP
1758             JMOHEP(2,IHEP)=KHEP
1759           ENDIF
1760         ENDIF
1761         IF (GLUIN) THEN
1762           IF (EMIT.EQ.1) THEN
1763             IF (IDHEP(IHEP).GT.0) THEN
1764               JMOHEP(2,IHEP)=JHEP
1765               JDAHEP(2,IHEP)=KHEP
1766               JMOHEP(2,JHEP)=KHEP
1767               JDAHEP(2,JHEP)=IHEP
1768               JMOHEP(2,KHEP)=IHEP
1769               JDAHEP(2,KHEP)=JHEP
1770             ELSE
1771               JMOHEP(2,IHEP)=KHEP
1772               JDAHEP(2,IHEP)=JHEP
1773               JMOHEP(2,JHEP)=IHEP
1774               JDAHEP(2,JHEP)=KHEP
1775               JMOHEP(2,KHEP)=JHEP
1776               JDAHEP(2,KHEP)=IHEP
1777             ENDIF
1778           ELSE
1779             IF (IDHEP(JHEP).GT.0) THEN
1780               JMOHEP(2,JHEP)=IHEP
1781               JDAHEP(2,JHEP)=KHEP
1782               JMOHEP(2,IHEP)=KHEP
1783               JDAHEP(2,IHEP)=JHEP
1784               JMOHEP(2,KHEP)=JHEP
1785               JDAHEP(2,KHEP)=IHEP
1786             ELSE
1787               JMOHEP(2,JHEP)=KHEP
1788               JDAHEP(2,JHEP)=IHEP
1789               JMOHEP(2,IHEP)=JHEP
1790               JDAHEP(2,IHEP)=KHEP
1791               JMOHEP(2,KHEP)=IHEP
1792               JDAHEP(2,KHEP)=JHEP
1793             ENDIF
1794           ENDIF
1795         ENDIF
1796         EMSCA=SQRT(EM**2+PHEP(1,ICMF)**2+PHEP(2,ICMF)**2)
1797 C--------SET STATUS AND LEPTON MOMENTA AFTER THE PARTON SHOWER
1798       ELSEIF (IOPT.EQ.2) THEN
1799         IF (EMIT.EQ.0.OR.NEVHEP+NWGTS.NE.NTMP) RETURN
1800         ISTHEP(JDAHEP(1,ICMF))=195
1801         IDHW(NHEP+1)=ID4
1802         IDHW(NHEP+2)=ID5
1803         IDHEP(NHEP+1)=IDPDG(ID4)
1804         IDHEP(NHEP+2)=IDPDG(ID5)
1805         ISTHEP(NHEP+1)=113
1806         ISTHEP(NHEP+2)=114
1807         CW=PHEP(3,ICMF)/SQRT(PHEP(1,ICMF)**2+PHEP(2,ICMF)**2+
1808      &       PHEP(3,ICMF)**2)
1809         SW=SQRT(1-CW**2)
1810         CALL HWUROT(PHEP(1,ICMF),CW,SW,R4)
1811         CALL HWUROF(R4,PHEP(1,ICMF),PR)
1812         PR(4)=PHEP(4,ICMF)
1813         CALL HWUMAS(PR)
1814         CALL HWUROF(R4,PS,PS)
1815         CALL HWUROF(R4,PF,PF)
1816         CALL HWUMAS(PS)
1817         CALL HWUMAS(PF)
1818         CALL HWUROT(PHEP(1,JDAHEP(1,ICMF)),CW,SW,R5)
1819         CALL HWUROF(R5,PHEP(1,JDAHEP(1,ICMF)),PD)
1820         PD(4)=PHEP(4,JDAHEP(1,ICMF))
1821         CALL HWUMAS(PD)
1822         BETA1=(PR(4)*PR(3)-SQRT(PR(4)**2*PD(3)**2-PR(3)**2*PD(3)**2+
1823      &       PD(3)**4))/(PD(3)**2+PR(4)**2)
1824         GAMMA1=1/SQRT(1-BETA1**2)
1825         PHEP(4,NHEP+1)=GAMMA1*PS(4)-BETA1*GAMMA1*PS(3)
1826         PHEP(3,NHEP+1)=-BETA1*GAMMA1*PS(4)+GAMMA1*PS(3)
1827         PHEP(4,NHEP+2)=GAMMA1*PF(4)-BETA1*GAMMA1*PF(3)
1828         PHEP(3,NHEP+2)=-BETA1*GAMMA1*PF(4)+GAMMA1*PF(3)
1829         PHEP(1,NHEP+1)=PS(1)
1830         PHEP(2,NHEP+1)=PS(2)
1831         PHEP(1,NHEP+2)=PF(1)
1832         PHEP(2,NHEP+2)=PF(2)
1833         CALL HWUMAS(PHEP(1,NHEP+1))
1834         CALL HWUMAS(PHEP(1,NHEP+2))
1835         CALL HWUROB(R5,PHEP(1,NHEP+1),PHEP(1,NHEP+1))
1836         CALL HWUROB(R5,PHEP(1,NHEP+2),PHEP(1,NHEP+2))
1837         JDAHEP(1,JDAHEP(1,ICMF))=NHEP+1
1838         JDAHEP(2,JDAHEP(1,ICMF))=NHEP+2
1839         JMOHEP(1,NHEP+1)=JDAHEP(1,ICMF)
1840         JMOHEP(1,NHEP+2)=JDAHEP(1,ICMF)
1841         JMOHEP(2,NHEP+1)=NHEP+2
1842         JDAHEP(2,NHEP+1)=NHEP+2
1843         JMOHEP(2,NHEP+2)=NHEP+1
1844         JDAHEP(2,NHEP+2)=NHEP+1
1845 C--special for spin correlations(relabel in spin common block)
1846         IF(SYSPIN.AND.NSPN.NE.0) THEN
1847           IDSPN(2) = NHEP+1
1848           IDSPN(3) = NHEP+2
1849           ISNHEP(NHEP+1) = 2
1850           ISNHEP(NHEP+2) = 3
1851         ENDIF
1852         NHEP=NHEP+2
1853         EMIT=0
1854       ENDIF
1855       END
1856 CDECK  ID>, HWBFIN.
1857 *CMZ :-        -26/04/91  10.18.56  by  Bryan Webber
1858 *-- Author :    Bryan Webber
1859 C-----------------------------------------------------------------------
1860       SUBROUTINE HWBFIN(IHEP)
1861 C-----------------------------------------------------------------------
1862 C     DELETES INTERNAL LINES FROM SHOWER, MAKES COLOUR CONNECTION INDEX
1863 C     AND COPIES INTO /HEPEVT/ IN COLOUR ORDER.
1864 C-----------------------------------------------------------------------
1865       INCLUDE 'HERWIG65.INC'
1866       INTEGER IHEP,ID,IJET,KHEP,IPAR,JPAR,NXPAR,IP,JP
1867       IF (IERROR.NE.0) RETURN
1868 C---SAVE VIRTUAL PARTON DATA
1869       NHEP=NHEP+1
1870       IF(NHEP.GT.NMXHEP) THEN
1871         CALL HWWARN('HWBFIN',100)
1872         GOTO 999
1873       ENDIF
1874       ID=IDPAR(2)
1875       IDHW(NHEP)=ID
1876       IDHEP(NHEP)=IDPDG(ID)
1877       ISTHEP(NHEP)=ISTHEP(IHEP)+20
1878       JMOHEP(1,NHEP)=IHEP
1879       JMOHEP(2,NHEP)=JMOHEP(1,IHEP)
1880       JDAHEP(1,IHEP)=NHEP
1881       JDAHEP(1,NHEP)=0
1882       JDAHEP(2,NHEP)=0
1883       CALL HWVEQU(5,PPAR(1,2),PHEP(1,NHEP))
1884       CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP))
1885 C---FINISHED FOR SPECTATOR OR NON-PARTON JETS
1886       IF (ISTHEP(NHEP).GT.136) RETURN
1887       IF (ID.GT.13.AND.ID.LT.209 .AND. ID.NE.59) RETURN
1888       IF (ID.GT.220.AND.ABS(IDPDG(ID)).LT.1000000) RETURN
1889       IF (ID.GT.424.AND.ID.NE.449) RETURN
1890       IF (.NOT.TMPAR(2).AND.ID.EQ.59) RETURN
1891       IDHEP(NHEP)=94
1892       IJET=NHEP
1893       IF (NPAR.GT.2) THEN
1894 C---SAVE CONE DATA
1895         NHEP=NHEP+1
1896         IF(NHEP.GT.NMXHEP) THEN
1897           CALL HWWARN('HWBFIN',101)
1898           GOTO 999
1899         ENDIF
1900         IDHW(NHEP)=IDPAR(1)
1901         IDHEP(NHEP)=0
1902         ISTHEP(NHEP)=100
1903         JMOHEP(1,NHEP)=IHEP
1904         JMOHEP(2,NHEP)=JCOPAR(1,1)
1905         JDAHEP(1,NHEP)=0
1906         JDAHEP(2,NHEP)=0
1907         CALL HWVEQU(5,PPAR,PHEP(1,NHEP))
1908         CALL HWVEQU(4,VPAR(1,2),VHEP(1,NHEP))
1909       ENDIF
1910       KHEP=NHEP
1911 C---START WITH ANTICOLOUR DAUGHTER OF HARDEST PARTON
1912       IPAR=2
1913       JPAR=JCOPAR(4,IPAR)
1914       NXPAR=NPAR/2
1915       DO 20 IP=1,NXPAR
1916       DO 10 JP=1,NXPAR
1917       IF (JPAR.EQ.0) GOTO 15
1918       IF (JCOPAR(2,JPAR).EQ.IPAR) THEN
1919         IPAR=JPAR
1920         JPAR=JCOPAR(4,IPAR)
1921       ELSE
1922         IPAR=JPAR
1923         JPAR=JCOPAR(1,IPAR)
1924       ENDIF
1925    10 CONTINUE
1926 C---COULDN'T FIND COLOUR PARTNER
1927       CALL HWWARN('HWBFIN',1)
1928    15 JPAR=JCOPAR(1,IPAR)
1929       KHEP=KHEP+1
1930       IF(KHEP.GT.NMXHEP) THEN
1931         CALL HWWARN('HWBFIN',102)
1932         GOTO 999
1933       ENDIF
1934       ID=IDPAR(IPAR)
1935       IF (TMPAR(IPAR)) THEN
1936         IF (ID.LT.14) THEN
1937           ISTHEP(KHEP)=139
1938         ELSEIF (ID.EQ.59) THEN
1939           ISTHEP(KHEP)=139
1940         ELSEIF (ID.LT.109) THEN
1941           ISTHEP(KHEP)=130
1942         ELSEIF (ID.LT.120) THEN
1943           ISTHEP(KHEP)=139
1944         ELSEIF (ABS(IDPDG(ID)).LT.1000000) THEN
1945           ISTHEP(KHEP)=130
1946         ELSEIF (ID.LT.425) THEN
1947           ISTHEP(KHEP)=139
1948         ELSEIF (ID.EQ.449) THEN
1949           ISTHEP(KHEP)=139
1950         ELSE
1951           ISTHEP(KHEP)=130
1952         ENDIF
1953       ELSE
1954         ISTHEP(KHEP)=ISTHEP(IHEP)+24
1955       ENDIF
1956       IDHW(KHEP)=ID
1957       IDHEP(KHEP)=IDPDG(ID)
1958       CALL HWVEQU(5,PPAR(1,IPAR),PHEP(1,KHEP))
1959       CALL HWVEQU(4,VPAR(1,IPAR),VHEP(1,KHEP))
1960       JMOHEP(1,KHEP)=IJET
1961       JMOHEP(2,KHEP)=KHEP+1
1962       JDAHEP(1,KHEP)=0
1963       JDAHEP(2,KHEP)=KHEP-1
1964    20 CONTINUE
1965       JMOHEP(2,KHEP)=0
1966       JDAHEP(2,NHEP+1)=0
1967       JDAHEP(1,IJET)=NHEP+1
1968       JDAHEP(2,IJET)=KHEP
1969       NHEP=KHEP
1970  999  RETURN
1971       END
1972 CDECK  ID>, HWBGEN.
1973 *CMZ :-        -14/10/99  18.04.56  by  Mike Seymour
1974 *-- Author :    Bryan Webber
1975 C-----------------------------------------------------------------------
1976       SUBROUTINE HWBGEN
1977 C-----------------------------------------------------------------------
1978 C     BRANCHING GENERATOR WITH INTERFERING GLUONS
1979 C     HWBGEN EVOLVES QCD JETS ACCORDING TO THE METHOD OF
1980 C     G.MARCHESINI & B.R.WEBBER, NUCL. PHYS. B238(1984)1
1981 C-----------------------------------------------------------------------
1982       INCLUDE 'HERWIG65.INC'
1983       DOUBLE PRECISION HWULDO,HWRGAU,EINHEP,ERTXI,RTXI,XF
1984       INTEGER NTRY,LASHEP,IHEP,NRHEP,ID,IST,JHEP,KPAR,I,J,IRHEP(NMXJET),
1985      & IRST(NMXJET),JPR
1986       LOGICAL HWRLOG
1987       EXTERNAL HWULDO,HWRGAU
1988       IF (IERROR.NE.0) RETURN
1989       IF (IPRO.EQ.80) RETURN
1990 C---CHECK THAT EMSCA IS SET
1991       IF (EMSCA.LE.ZERO) CALL HWWARN('HWBGEN',200)
1992       IF (HARDME) THEN
1993 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN E+E-
1994         JPR=IPROC/10
1995 C**********13/11/00 BRW FIX TO ALLOW ALSO WW AND ZZ
1996         IF (JPR.EQ.10.OR.JPR.EQ.20.OR.JPR.EQ.25) CALL HWBDED(1)
1997 C**********END FIX
1998 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DIS
1999         IF (IPRO.EQ.90) CALL HWBDIS(1)
2000 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DRELL-YAN PROCESSES
2001         IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(1)
2002 C---FORCE A BRANCH INTO THE `DEAD ZONE' IN TOP DECAYS
2003         CALL HWBTOP
2004       ENDIF
2005 C---GENERATE INTRINSIC PT ONCE AND FOR ALL
2006       DO 5 JNHAD=1,2
2007         IF (PTRMS.NE.0.) THEN
2008           PTINT(1,JNHAD)=HWRGAU(1,ZERO,PXRMS)
2009           PTINT(2,JNHAD)=HWRGAU(2,ZERO,PXRMS)
2010           PTINT(3,JNHAD)=PTINT(1,JNHAD)**2+PTINT(2,JNHAD)**2
2011         ELSE
2012           CALL HWVZRO(3,PTINT(1,JNHAD))
2013         ENDIF
2014  5    CONTINUE
2015       NTRY=0
2016       LASHEP=NHEP
2017  10   NTRY=NTRY+1
2018       IF (NTRY.GT.NETRY) THEN
2019         CALL HWWARN('HWBGEN',ISLENT*100)
2020         GOTO 999
2021       ENDIF
2022       NRHEP=0
2023       NHEP=LASHEP
2024       FROST=.FALSE.
2025       DO 100 IHEP=1,LASHEP
2026       IST=ISTHEP(IHEP)
2027       IF (IST.GE.111.AND.IST.LE.115) THEN
2028        NRHEP=NRHEP+1
2029        IRHEP(NRHEP)=IHEP
2030        IRST(NRHEP)=IST
2031        ID=IDHW(IHEP)
2032        IF (IST.NE.115) THEN
2033 C---FOUND A PARTON TO EVOLVE
2034         NEVPAR=IHEP
2035         NPAR=2
2036         IDPAR(1)=17
2037         IDPAR(2)=ID
2038         TMPAR(1)=.TRUE.
2039         PPAR(2,1)=0.
2040         PPAR(4,1)=1.
2041         DO 15 J=1,2
2042         DO 15 I=1,2
2043         JMOPAR(I,J)=0
2044  15     JCOPAR(I,J)=0
2045 C---SET UP EVOLUTION SCALE AND FRAME
2046         JHEP=JMOHEP(2,IHEP)
2047         IF (ID.EQ.13) THEN
2048           IF (HWRLOG(HALF)) JHEP=JDAHEP(2,IHEP)
2049         ELSEIF (IST.GT.112) THEN
2050           IF ((ID.GT.6.AND.ID.LT.13).OR.
2051      &        (ID.GT.214.AND.ID.LT.221)) JHEP=JDAHEP(2,IHEP)
2052         ELSE
2053           IF (ID.LT.7.OR.(ID.GT.208.AND.ID.LT.215)) JHEP=JDAHEP(2,IHEP)
2054         ENDIF
2055         IF (JHEP.LE.0.OR.JHEP.GT.NHEP) THEN
2056           CALL HWWARN('HWBGEN',1)
2057           JHEP=IHEP
2058         ENDIF
2059         JCOPAR(1,1)=JHEP
2060         EINHEP=PHEP(4,IHEP)
2061         ERTXI=HWULDO(PHEP(1,IHEP),PHEP(1,JHEP))
2062         IF (ERTXI.LT.ZERO) ERTXI=0.
2063         IF (IST.LE.112.AND.IHEP.EQ.JHEP) ERTXI=0.
2064         IF (ISTHEP(JHEP).EQ.155) THEN
2065           ERTXI=ERTXI/PHEP(5,JHEP)
2066           RTXI=1.
2067         ELSE
2068           ERTXI=SQRT(ERTXI)
2069           RTXI=ERTXI/EINHEP
2070         ENDIF
2071         IF (RTXI.EQ.ZERO) THEN
2072           XF=1.
2073           PPAR(1,1)=0.
2074           PPAR(3,1)=1.
2075           PPAR(1,2)=EINHEP
2076           PPAR(2,2)=0.
2077           PPAR(4,2)=EINHEP
2078         ELSE
2079           XF=1./RTXI
2080           PPAR(1,1)=1.
2081           PPAR(3,1)=0.
2082           PPAR(1,2)=ERTXI
2083           PPAR(2,2)=1.
2084           PPAR(4,2)=ERTXI
2085         ENDIF
2086         IF (PPAR(4,2).LT.PHEP(5,IHEP)) PPAR(4,2)=PHEP(5,IHEP)
2087 C---STORE MASS
2088         PPAR(5,2)=PHEP(5,IHEP)
2089         CALL HWVZRO(4,VPAR(1,1))
2090         CALL HWVZRO(4,VPAR(1,2))
2091         IF (IST.GT.112) THEN
2092           TMPAR(2)=.TRUE.
2093           INHAD=0
2094           JNHAD=0
2095           XFACT=0.
2096         ELSE
2097           TMPAR(2)=.FALSE.
2098           JNHAD=IST-110
2099           INHAD=JNHAD
2100           IF (JDAHEP(1,JNHAD).NE.0) INHAD=JDAHEP(1,JNHAD)
2101           XFACT=XF/PHEP(4,INHAD)
2102           ANOMSC(1,JNHAD)=ZERO
2103           ANOMSC(2,JNHAD)=ZERO
2104         ENDIF
2105 C---FOR QUARKS IN A COLOUR SINGLET, ALLOW SOFT MATRIX-ELEMENT CORRECTION
2106         HARDST=PPAR(4,2)
2107         IF (SOFTME.AND.IDHW(IHEP).LT.13.AND.
2108      $       ((JMOHEP(2,JHEP).EQ.IHEP.AND.JDAHEP(2,JHEP).EQ.IHEP).OR.
2109      $       ISTHEP(JHEP).EQ.155)) HARDST=0
2110 C---CREATE BRANCHES AND COMPUTE ENERGIES
2111         DO 20 KPAR=2,NMXPAR
2112         IF (TMPAR(KPAR)) THEN
2113           CALL HWBRAN(KPAR)
2114         ELSE
2115           CALL HWSBRN(KPAR)
2116         ENDIF
2117         IF (IERROR.NE.0) RETURN
2118         IF (FROST) GOTO 100
2119         IF (KPAR.EQ.NPAR) GOTO 30
2120  20     CONTINUE
2121 C---COMPUTE MASSES AND 3-MOMENTA
2122  30     CONTINUE
2123         CALL HWBMAS
2124         IF (AZSPIN) CALL HWBSPN
2125         IF (TMPAR(2)) THEN
2126            CALL HWBTIM(2,1)
2127         ELSE
2128            CALL HWBSPA
2129         ENDIF
2130 C---ENTER PARTON JET IN /HEPEVT/
2131         CALL HWBFIN(IHEP)
2132        ELSE
2133 C---COPY SPECTATOR
2134         NHEP=NHEP+1
2135         IF (ID.GT.120.AND.ID.LT.133 .OR. ID.GE.198.AND.ID.LE.201) THEN
2136           ISTHEP(NHEP)=190
2137         ELSE
2138           ISTHEP(NHEP)=152
2139         ENDIF
2140         IDHW(NHEP)=ID
2141         IDHEP(NHEP)=IDPDG(ID)
2142         JMOHEP(1,NHEP)=IHEP
2143         JMOHEP(2,NHEP)=0
2144         JDAHEP(2,NHEP)=0
2145         JDAHEP(1,IHEP)=NHEP
2146         CALL HWVEQU(5,PHEP(1,IHEP),PHEP(1,NHEP))
2147        ENDIF
2148        ISTHEP(IHEP)=ISTHEP(IHEP)+10
2149       ENDIF
2150  100  CONTINUE
2151       IF (.NOT.FROST) THEN
2152 C---COMBINE JETS
2153         ISTAT=20
2154         CALL HWBJCO
2155       ENDIF
2156       IF (.NOT.FROST) THEN
2157 C---ATTACH SPECTATORS
2158         ISTAT=30
2159         CALL HWSSPC
2160       ENDIF
2161       IF (FROST) THEN
2162 C---BAD JET: RESTORE PARTONS AND RE-EVOLVE
2163          DO 120 I=1,NRHEP
2164  120     ISTHEP(IRHEP(I))=IRST(I)
2165          GOTO 10
2166       ENDIF
2167 C---CONNECT COLOURS
2168       CALL HWBCON
2169       ISTAT=40
2170       LASHEP=NHEP
2171       IF (HARDME) THEN
2172 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN E+E-
2173         IF (IPROC/10.EQ.10) CALL HWBDED(2)
2174 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DIS
2175         IF (IPRO.EQ.90) CALL HWBDIS(2)
2176 C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DRELL-YAN PROC
2177         IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(2)
2178       ENDIF
2179 C---IF THE CLEAN-UP OPERATION ADDED ANY PARTONS TO THE EVENT RECORD
2180 C   IT MIGHT NEED RESHOWERING
2181       IF (NHEP.GT.LASHEP) THEN
2182         LASHEP=NHEP
2183         GOTO 10
2184       ENDIF
2185  999  RETURN
2186       END
2187 CDECK  ID>, HWBGUP.
2188 *CMZ :-        -16/07/02  09.40.25  by  Peter Richardson
2189 *-- Author :    Peter Richardson
2190 C----------------------------------------------------------------------
2191       SUBROUTINE HWBGUP(ISTART,ICMF)
2192 C----------------------------------------------------------------------
2193 C     Makes the colour connections and performs the parton shower
2194 C     for events read in from the GUPI (Generic User Process Interface)
2195 C     event common block
2196 C----------------------------------------------------------------------
2197       INCLUDE 'HERWIG65.INC'
2198       INTEGER MAXNUP
2199       PARAMETER (MAXNUP=500)
2200       INTEGER NUP,IDPRUP,IDUP,ISTUP,MOTHUP,ICOLUP
2201       DOUBLE PRECISION XWGTUP,SCALUP,AQEDUP,AQCDUP,PUP,VTIMUP,SPINUP
2202       COMMON/HEPEUP/NUP,IDPRUP,XWGTUP,SCALUP,AQEDUP,AQCDUP,
2203      &              IDUP(MAXNUP),ISTUP(MAXNUP),MOTHUP(2,MAXNUP),
2204      &              ICOLUP(2,MAXNUP),PUP(5,MAXNUP),VTIMUP(MAXNUP),
2205      &              SPINUP(MAXNUP)
2206 C--Local variables
2207       INTEGER ISTART,ICMF,J,K,I,JCOL,ICOL
2208       LOGICAL FOUND
2209       COMMON /HWGUP/ILOC(NMXHEP),JLOC(MAXNUP)
2210       INTEGER ILOC,JLOC
2211 C--now we need to do the colour connections
2212  20   ISTART = ISTART+1
2213       IF(ISTART.GT.NHEP) GOTO 30
2214       IF(ISTART.EQ.ICMF) ISTART = ISTART+1
2215       IF(JMOHEP(2,ISTART).NE.0.AND.JDAHEP(2,ISTART).NE.0) GOTO 20
2216       K = ISTART
2217       J = ILOC(K)
2218       IF(ICOLUP(1,J).NE.0) THEN
2219         JCOL = 1
2220         ICOL = ICOLUP(1,J)
2221       ELSE
2222         JCOL = 2
2223         ICOL = ICOLUP(2,J)
2224       ENDIF
2225       IF(ICOL.EQ.0) THEN
2226         JMOHEP(2,K) = K
2227         JDAHEP(2,K) = K
2228         GOTO 20
2229       ENDIF
2230 C--now search for the partner
2231 C--first search for the flavour partner if not looking for colour partner
2232 C--search for the flavour partner of the particle
2233 C--this must be set or HERWIG won't work
2234  10   IF(JDAHEP(2,K).NE.0.AND.JMOHEP(2,K).NE.0) GOTO 20
2235       IF(ICOL.EQ.0) THEN
2236         FOUND = .FALSE.
2237 C--look for unpaired particle
2238         DO 15 I=1,NUP
2239           IF(JLOC(I).EQ.0) GOTO 15
2240           IF(IDUP(I).EQ.21.OR.IDUP(I).EQ.9) GOTO 15
2241           IF(JLOC(I).EQ.ISTART) GOTO 15
2242           IF(ICOLUP(1,I).EQ.0.AND.ICOLUP(2,I).EQ.0) GOTO 15
2243 C--antiflavour partner
2244           IF(JDAHEP(2,JLOC(I)).EQ.0) THEN
2245 C--pair incoming     particle with outgoing     particle
2246 C-- or  outgoing antiparticle with outgoing     particle
2247             IF(ISTUP(I).GT.0.AND.IDUP(I).GT.0.AND.
2248      &         ((IDUP(J).GT.0.AND.ISTUP(J).EQ.-1).OR.
2249      &          (IDUP(J).LT.0.AND.ISTUP(J).GT.0 )))  THEN
2250               FOUND = .TRUE.
2251               JCOL = 1
2252 C--pair incoming     particle with incoming antiparticle
2253 C-- or  outgoing antiparticle with incoming antiparticle
2254             ELSEIF(IDUP(I).LT.0.AND.ISTUP(I).EQ.-1.AND.
2255      &             ((IDUP(J).GT.0.AND.ISTUP(J).EQ.-1).OR.
2256      &              (IDUP(J).LT.0.AND.ISTUP(J).GT.0 ))) THEN
2257               FOUND = .TRUE.
2258               JCOL = 2
2259             ENDIF
2260 C--make the connection
2261             IF(FOUND) THEN
2262               JMOHEP(2,K)       = JLOC(I)
2263               JDAHEP(2,JLOC(I)) = K
2264             ENDIF
2265           ENDIF
2266 C--flavour partner
2267           IF(JMOHEP(2,JLOC(I)).EQ.0.AND.(.NOT.FOUND)) THEN
2268 C--pair incoming antiparticle with outgoing antiparticle
2269 C-- or  outgoing     particle with outgoing antiparticle
2270             IF(IDUP(I).LT.0.AND.ISTUP(I).GT.0.AND.
2271      &         ((IDUP(J).LT.0.AND.ISTUP(J).EQ.-1).OR.
2272      &          (IDUP(J).GT.0.AND.ISTUP(J).GT.0 ))) THEN
2273               FOUND = .TRUE.
2274               JCOL = 2
2275 C--pair incoming antiparticle with incoming     particle
2276 C-- or  outgoing     particle with incoming     particle
2277             ELSEIF(IDUP(I).GT.0.AND.ISTUP(I).EQ.-1.AND.
2278      &             ((IDUP(J).LT.0.AND.ISTUP(J).EQ.-1).OR.
2279      &              (IDUP(J).GT.0.AND.ISTUP(J).GT.0 ))) THEN
2280               FOUND = .TRUE.
2281               JCOL = 1
2282             ENDIF
2283 C--make the connection
2284             IF(FOUND) THEN
2285               JDAHEP(2,K) = JLOC(I)
2286               JMOHEP(2,JLOC(I)) = K
2287             ENDIF
2288           ENDIF
2289 C--set up the search for the next partner
2290           IF(FOUND) THEN
2291             FOUND = .FALSE.
2292             ICOL = ICOLUP(JCOL,I)
2293             K = JLOC(I)
2294             J = I
2295             GOTO 10
2296           ENDIF
2297  15     CONTINUE
2298 C--if no other choice then connect to the first particle in the loop
2299         IF(JDAHEP(2,K).EQ.0.AND.JMOHEP(2,ISTART).EQ.0) THEN
2300            JDAHEP(2,K) = ISTART
2301            JMOHEP(2,ISTART) = K
2302         ELSEIF(JDAHEP(2,ISTART).EQ.0.AND.JMOHEP(2,K).EQ.0) THEN
2303            JMOHEP(2,K) = ISTART
2304            JDAHEP(2,ISTART) = K
2305         ELSE
2306           CALL HWWARN('HWBGUP',100)
2307           GOTO 999
2308         ENDIF
2309         GOTO 20
2310       ENDIF
2311 C--now the bit to find colour partners
2312       FOUND = .FALSE.
2313 C--special for particle from a decaying coloured particle
2314       IF(MOTHUP(1,J).NE.0) THEN
2315         IF(ISTUP(MOTHUP(1,J)).EQ.2.OR.ISTUP(MOTHUP(1,J)).EQ.3) THEN
2316           IF(IDUP(J).LT.0.AND.ICOL.EQ.ICOLUP(2,MOTHUP(1,J))) THEN
2317             JDAHEP(2,K) = JLOC(MOTHUP(1,J))
2318             JMOHEP(2,K) = JLOC(MOTHUP(1,J))
2319             GOTO 20
2320           ELSEIF(IDUP(J).GT.0.AND.ICOL.EQ.ICOLUP(1,MOTHUP(1,J))) THEN
2321             JDAHEP(2,K) = JLOC(MOTHUP(1,J))
2322             JMOHEP(2,K) = JLOC(MOTHUP(1,J))
2323             GOTO 20
2324           ENDIF
2325         ENDIF
2326       ENDIF
2327 C--search for the partner
2328       DO I=1,NUP
2329         IF(ICOLUP(1,I).EQ.ICOL.AND.I.NE.J) THEN
2330           IF((JCOL.EQ.1.AND.ISTUP(J).EQ.-1.AND.ISTUP(I).GT.0).OR.
2331      &       (JCOL.EQ.2.AND.ISTUP(J).GT.0.AND.ISTUP(I).GE.0)) THEN
2332             JDAHEP(2,K)       = JLOC(I)
2333             JMOHEP(2,JLOC(I)) = K
2334             FOUND = .TRUE.
2335           ELSEIF((JCOL.EQ.1.AND.ISTUP(J).GT.0.AND.ISTUP(I).EQ.-1).OR.
2336      &          (JCOL.EQ.2.AND.ISTUP(J).EQ.-1.AND.ISTUP(I).EQ.-1)) THEN
2337             JMOHEP(2,K)       = JLOC(I)
2338             JDAHEP(2,JLOC(I)) = K
2339             FOUND = .TRUE.
2340           ENDIF
2341           IF(FOUND) JCOL = 2
2342         ELSEIF(ICOLUP(2,I).EQ.ICOL.AND.I.NE.J) THEN
2343           IF((JCOL.EQ.1.AND.ISTUP(J).EQ.-1.AND.ISTUP(I).EQ.-1).OR.
2344      &       (JCOL.EQ.2.AND.ISTUP(J).GT.0.AND.ISTUP(I).EQ.-1)) THEN
2345             JDAHEP(2,K) = JLOC(I)
2346             JMOHEP(2,JLOC(I)) = K
2347             FOUND = .TRUE.
2348           ELSEIF((JCOL.EQ.1.AND.ISTUP(J).GE.0.AND.ISTUP(I).GE.0).OR.
2349      &           (JCOL.EQ.2.AND.ISTUP(J).EQ.-1.AND.ISTUP(I).GE.0)) THEN
2350             JMOHEP(2,K) = JLOC(I)
2351             JDAHEP(2,JLOC(I)) = K
2352             FOUND = .TRUE.
2353           ENDIF
2354           IF(FOUND) JCOL = 1
2355         ENDIF
2356         IF(FOUND) THEN
2357           K = JLOC(I)
2358           J = I
2359           ICOL = ICOLUP(JCOL,I)
2360           GOTO 10
2361         ENDIF
2362       ENDDO
2363 C--special for self connected gluons
2364       IF(IDUP(J).EQ.21.OR.IDUP(J).EQ.9.AND.
2365      &     ICOLUP(1,J).EQ.ICOLUP(2,J)) THEN
2366         JMOHEP(2,K) = K
2367         JDAHEP(2,K) = K
2368 C--options for self connected gluons
2369         IF(LHGLSF) THEN
2370           CALL HWWARN('HWBGUP',1)
2371         ELSE
2372           CALL HWWARN('HWBGUP',101)
2373           GOTO 999
2374         ENDIF
2375         GOTO 20
2376       ENDIF
2377 C--perform the shower
2378  30   CALL HWBGEN
2379  999  RETURN
2380       END
2381 CDECK  ID>, HWBJCO.
2382 *CMZ :-        -30/09/02  09.19.58  by  Peter Richardson
2383 *-- Author :    Bryan Webber
2384 C-----------------------------------------------------------------------
2385       SUBROUTINE HWBJCO
2386 C-----------------------------------------------------------------------
2387 C     COMBINES JETS WITH REQUIRED KINEMATICS
2388 C-----------------------------------------------------------------------
2389       INCLUDE 'HERWIG65.INC'
2390       DOUBLE PRECISION HWULDO,EPS,PTX,PTY,PF,PTINF,PTCON,CN,CP,SP,PP0,
2391      & PM0,ET0,DET,ECM,EMJ,EMP,EMS,DMS,ES,DPF,ALF,AL(2),ET(2),PP(2),
2392      & PT(3),PA(5),PB(5),PC(5),PQ(5),PR(5),PS(5),RR(3,3),RS(3,3),ETC,
2393      & PJ(NMXJET),PM(NMXJET),PBR(5),RBR(3,3),DISP(4),PLAB(5)
2394       INTEGER LJET,IJ1,IST,IP,ICM,IP1,IP2,NP,IHEP,MHEP,JP,KP,LP,KHEP,
2395      & JHEP,NE,IJT,IEND(2),IJET(NMXJET),IPAR(NMXJET)
2396       LOGICAL AZCOR,JETRAD,DISPRO,DISLOW
2397       EXTERNAL HWULDO
2398       PARAMETER (EPS=1.D-4)
2399       IF (IERROR.NE.0) RETURN
2400       AZCOR=AZSOFT.OR.AZSPIN
2401       LJET=131
2402   10  IJET(1)=1
2403   20  IJ1=IJET(1)
2404       DO 40 IHEP=IJ1,NHEP
2405       IST=ISTHEP(IHEP)
2406       IF (IST.EQ.137.OR.IST.EQ.138) IST=133
2407       IF (IST.EQ.LJET) THEN
2408 C---FOUND AN UNBOOSTED JET - FIND PARTNERS
2409         IP=JMOHEP(1,IHEP)
2410         ICM=JMOHEP(1,IP)
2411         DISPRO=IPRO/10.EQ.9.AND.IDHW(ICM).EQ.15
2412         DISLOW=DISPRO.AND.JDAHEP(1,ICM).EQ.JDAHEP(2,ICM)-1
2413         IF (IST.EQ.131) THEN
2414           IP1=JMOHEP(1,ICM)
2415           IP2=JMOHEP(2,ICM)
2416         ELSE
2417           IP1=JDAHEP(1,ICM)
2418           IP2=JDAHEP(2,ICM)
2419         ENDIF
2420         IF (IP1.NE.IP) THEN
2421           CALL HWWARN('HWBJCO',100)
2422           GOTO 999
2423         ENDIF
2424         NP=0
2425         DO 30 JHEP=IP1,IP2
2426         NP=NP+1
2427         IPAR(NP)=JHEP
2428   30    IJET(NP)=JDAHEP(1,JHEP)
2429         GOTO 50
2430       ENDIF
2431   40  CONTINUE
2432 C---NO MORE JETS?
2433       IF (LJET.EQ.131) THEN
2434         LJET=133
2435         GOTO 10
2436       ENDIF
2437       RETURN
2438   50  IF (LJET.EQ.131) THEN
2439 C---SPACELIKE JETS: FIND SPACELIKE PARTONS
2440         IF (NP.NE.2) THEN
2441           CALL HWWARN('HWBJCO',103)
2442           GOTO 999
2443         ENDIF
2444 C---special for DIS: FIND BOOST AND ROTATION FROM LAB TO BREIT FRAME
2445         IF (DISPRO.AND.BREIT) THEN
2446           IP=2
2447           IF (JDAHEP(1,IP).NE.0) IP=JDAHEP(1,IP)
2448           CALL HWVDIF(4,PHEP(1,JMOHEP(1,ICM)),PHEP(1,JDAHEP(1,ICM)),PB)
2449           CALL HWUMAS(PB)
2450 C---IF Q**2<10**-2, SOMETHING MUST HAVE ALREADY GONE WRONG
2451           IF (PB(5)**2.LT.1.D-2) THEN
2452             CALL HWWARN('HWBJCO',102)
2453             GOTO 999
2454           ENDIF
2455           CALL HWVSCA(4,PB(5)**2/HWULDO(PHEP(1,IP),PB),PHEP(1,IP),PBR)
2456           CALL HWVSUM(4,PB,PBR,PBR)
2457           CALL HWUMAS(PBR)
2458           CALL HWULOF(PBR,PB,PB)
2459           CALL HWUROT(PB,ONE,ZERO,RBR)
2460         ENDIF
2461         PTX=0.
2462         PTY=0.
2463         PF=1.D0
2464         DO 90 IP=1,2
2465         MHEP=IJET(IP)
2466         IF (JDAHEP(1,MHEP).EQ.0) THEN
2467 C---SPECIAL FOR NON-PARTON JETS
2468           IHEP=MHEP
2469           GOTO 70
2470         ELSE
2471           IST=134+IP
2472           DO 60 IHEP=MHEP,NHEP
2473   60      IF (ISTHEP(IHEP).EQ.IST) GOTO 70
2474 C---COULDN'T FIND SPACELIKE PARTON
2475           CALL HWWARN('HWBJCO',101)
2476           GOTO 999
2477         ENDIF
2478   70    CALL HWVSCA(3,PF,PHEP(1,IHEP),PS)
2479         IF (PTINT(3,IP).GT.ZERO) THEN
2480 C---ADD INTRINSIC PT
2481           PT(1)=PTINT(1,IP)
2482           PT(2)=PTINT(2,IP)
2483           PT(3)=0.
2484           CALL HWUROT(PS, ONE,ZERO,RS)
2485           CALL HWUROB(RS,PT,PT)
2486           CALL HWVSUM(3,PS,PT,PS)
2487         ENDIF
2488         JP=IJET(IP)+1
2489         IF (AZCOR.AND.JP.LE.NHEP.AND.IDHW(JP).EQ.17) THEN
2490 C---ALIGN CONE WITH INTERFERING PARTON
2491           CALL HWUROT(PS, ONE,ZERO,RS)
2492           CALL HWUROF(RS,PHEP(1,JP),PR)
2493           PTCON=PR(1)**2+PR(2)**2
2494           KP=JMOHEP(2,JP)
2495           IF (KP.EQ.0) THEN
2496             CALL HWWARN('HWBJCO',1)
2497             PTINF=0.
2498           ELSE
2499             CALL HWVEQU(4,PHEP(1,KP),PB)
2500             IF (DISPRO.AND.BREIT) THEN
2501               CALL HWULOF(PBR,PB,PB)
2502               CALL HWUROF(RBR,PB,PB)
2503             ENDIF
2504             PTINF=PB(1)**2+PB(2)**2
2505             IF (PTINF.LT.EPS) THEN
2506 C---COLLINEAR JETS: ALIGN CONES
2507               KP=JDAHEP(1,KP)+1
2508 C---BUG FIX BY MHS 17/03/05: RETURNED TO VERSION 6.500!
2509               IF (ISTHEP(KP).EQ.100.AND.ISTHEP(KP-1).GE.141
2510      $                             .AND.ISTHEP(KP-1).LE.144) THEN
2511 C---END FIX
2512                 CALL HWVEQU(4,PHEP(1,KP),PB)
2513                 IF (DISPRO.AND.BREIT) THEN
2514                   CALL HWULOF(PBR,PB,PB)
2515                   CALL HWUROF(RBR,PB,PB)
2516                 ENDIF
2517                 PTINF=PB(1)**2+PB(2)**2
2518               ELSE
2519                 PTINF=0.
2520               ENDIF
2521             ENDIF
2522           ENDIF
2523           IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN
2524             CN=1./SQRT(PTINF*PTCON)
2525             CP=CN*(PR(1)*PB(1)+PR(2)*PB(2))
2526             SP=CN*(PR(1)*PB(2)-PR(2)*PB(1))
2527           ELSE
2528             CALL HWRAZM( ONE,CP,SP)
2529           ENDIF
2530         ELSE
2531           CALL HWRAZM( ONE,CP,SP)
2532         ENDIF
2533 C---ROTATE SO SPACELIKE IS ALONG AXIS (APART FROM INTRINSIC PT)
2534         CALL HWUROT(PS,CP,SP,RS)
2535         IHEP=IJET(IP)
2536         KHEP=JDAHEP(2,IHEP)
2537         IF (KHEP.LT.IHEP) KHEP=IHEP
2538         IEND(IP)=KHEP
2539         DO 80 JHEP=IHEP,KHEP
2540         CALL HWUROF(RS,PHEP(1,JHEP),PHEP(1,JHEP))
2541   80    CALL HWUROF(RS,VHEP(1,JHEP),VHEP(1,JHEP))
2542         PP(IP)=PHEP(4,IHEP)+PF*PHEP(3,IHEP)
2543         ET(IP)=PHEP(1,IHEP)**2+PHEP(2,IHEP)**2-PHEP(5,IHEP)**2
2544 C---REDEFINE HARD CM
2545         PTX=PTX+PHEP(1,IHEP)
2546         PTY=PTY+PHEP(2,IHEP)
2547   90    PF=-PF
2548         PHEP(1,ICM)=PTX
2549         PHEP(2,ICM)=PTY
2550 C---special for DIS: keep lepton momenta fixed
2551         IF (DISPRO) THEN
2552           IP1=JMOHEP(1,ICM)
2553           IP2=JDAHEP(1,ICM)
2554           IJT=IJET(1)
2555 C---IJT will be used to store lepton momentum transfer
2556           CALL HWVDIF(4,PHEP(1,IP1),PHEP(1,IP2),PHEP(1,IJT))
2557           CALL HWUMAS(PHEP(1,IJT))
2558           IF (IDHEP(IP1).EQ.IDHEP(IP2)) THEN
2559             IDHW(IJT)=200
2560           ELSEIF (IDHEP(IP1).LT.IDHEP(IP2)) THEN
2561             IDHW(IJT)=199
2562           ELSE
2563             IDHW(IJT)=198
2564           ENDIF
2565           IDHEP(IJT)=IDPDG(IDHW(IJT))
2566           ISTHEP(IJT)=3
2567 C---calculate boost for struck parton
2568 C   PC is momentum of outgoing parton(s)
2569           IP2=JDAHEP(2,ICM)
2570           IF (.NOT.DISLOW) THEN
2571 C---FOR heavy QQbar PQ and PC are old and new QQbar momenta
2572             CALL HWVSUM(4,PHEP(1,IP2-1),PHEP(1,IP2),PQ)
2573             CALL HWUMAS(PQ)
2574             PC(5)=PQ(5)
2575           ELSE
2576             PC(5)=PHEP(5,JDAHEP(1,IP2))
2577           ENDIF
2578           CALL HWVSUM(2,PHEP(1,IJT),PHEP(1,IJET(2)),PC)
2579           ET(1)=ET(2)
2580 C---USE BREIT FRAME BOSON MOMENTUM IF NECESSARY
2581           IF (BREIT) THEN
2582             ET(2)=ET(1)+PC(5)**2+PHEP(5,IJET(2))**2
2583             PM0=PHEP(5,IJT)
2584             PP0=-PM0
2585           ELSE
2586             ET(2)=PC(1)**2+PC(2)**2+PC(5)**2
2587             PP0=PHEP(4,IJT)+PHEP(3,IJT)
2588             PM0=PHEP(4,IJT)-PHEP(3,IJT)
2589           ENDIF
2590           ET0=(PP0*PM0)+ET(1)-ET(2)
2591           DET=ET0**2-4.*(PP0*PM0)*ET(1)
2592           IF (DET.LT.ZERO) THEN
2593             FROST=.TRUE.
2594             RETURN
2595           ENDIF
2596           ALF=(SQRT(DET)-ET0)/(2.*PP0*PP(2))
2597           PB(1)=0.
2598           PB(2)=0.
2599           PB(5)=2.D0
2600           PB(3)=ALF-(1./ALF)
2601           PB(4)=ALF+(1./ALF)
2602           DO 100 IHEP=IJET(2),IEND(2)
2603           CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
2604           CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
2605 C---BOOST FROM BREIT FRAME IF NECESSARY
2606           IF (BREIT) THEN
2607             CALL HWUROB(RBR,PHEP(1,IHEP),PHEP(1,IHEP))
2608             CALL HWULOB(PBR,PHEP(1,IHEP),PHEP(1,IHEP))
2609             CALL HWUROB(RBR,VHEP(1,IHEP),VHEP(1,IHEP))
2610             CALL HWULB4(PBR,VHEP(1,IHEP),VHEP(1,IHEP))
2611           ENDIF
2612   100     ISTHEP(IHEP)=ISTHEP(IHEP)+10
2613           CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJET(2)),DISP)
2614           DO 110 IHEP=IJET(2),IEND(2)
2615   110     CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
2616           IF (IEND(2).GT.IJET(2)+1) ISTHEP(IJET(2)+1)=100
2617           CALL HWVSUM(4,PHEP(1,IJT),PHEP(1,IJET(2)),PC)
2618           CALL HWVSUM(4,PHEP(1,IP1),PHEP(1,IJET(2)),PHEP(1,ICM))
2619           CALL HWUMAS(PHEP(1,ICM))
2620         ELSEIF (IPRO/10.EQ.5) THEN
2621 C Special to preserve photon momentum
2622            ETC=PTX**2+PTY**2+PHEP(5,ICM)**2
2623            ET0=ETC+ET(1)-ET(2)
2624            DET=ET0**2-4.*ETC*ET(1)
2625            IF (DET.LT.ZERO) THEN
2626               FROST=.TRUE.
2627               RETURN
2628            ENDIF
2629            ALF=(SQRT(DET)+ET0-2.*ET(1))/(2.*PP(1)*PP(2))
2630            PB(1)=0.
2631            PB(2)=0.
2632            PB(3)=ALF-1./ALF
2633            PB(4)=ALF+1./ALF
2634            PB(5)=2.
2635            IJT=IJET(2)
2636            DO 120 IHEP=IJT,IEND(2)
2637            CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
2638            CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
2639   120      ISTHEP(IHEP)=ISTHEP(IHEP)+10
2640            CALL HWVDIF(4,VHEP(1,IPAR(2)),VHEP(1,IJT),DISP)
2641            DO 130 IHEP=IJT,IEND(2)
2642   130      CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
2643            IF (IEND(2).GT.IJT+1) ISTHEP(IJT+1)=100
2644            ISTHEP(IJET(1))=ISTHEP(IJET(1))+10
2645            CALL HWVSUM(2,PHEP(3,IPAR(1)),PHEP(3,IJT),PHEP(3,ICM))
2646         ELSE
2647 C--change to preserve either long mom or rapidity rather than long mom
2648 C--by PR and BRW 30/9/02
2649           IF (PRESPL) THEN
2650 C--PRESERVE LONG MOM OF CMF
2651             PHEP(4,ICM)=
2652      &            SQRT(PTX**2+PTY**2+PHEP(3,ICM)**2+PHEP(5,ICM)**2)
2653           ELSE
2654 C--PRESERVE RAPIDITY OF CMF
2655             DET=SQRT(ONE+(PTX**2+PTY**2)/(PHEP(4,ICM)**2
2656      &                -PHEP(3,ICM)**2))
2657             CALL HWVSCA(2,DET,PHEP(3,ICM),PHEP(3,ICM))
2658           ENDIF
2659 C---NOW BOOST TO REQUIRED Q**2 AND X-F
2660           PP0=PHEP(4,ICM)+PHEP(3,ICM)
2661           PM0=PHEP(4,ICM)-PHEP(3,ICM)
2662           ET0=(PP0*PM0)+ET(1)-ET(2)
2663           DET=ET0**2-4.*(PP0*PM0)*ET(1)
2664           IF (DET.LT.ZERO) THEN
2665             FROST=.TRUE.
2666             RETURN
2667           ENDIF
2668           DET=SQRT(DET)+ET0
2669           AL(1)= 2.*PM0*PP(1)/DET
2670           AL(2)=(PM0/PP(2))*(1.-2.*ET(1)/DET)
2671           PB(1)=0.
2672           PB(2)=0.
2673           PB(5)=2.
2674           DO 160 IP=1,2
2675           PB(3)=AL(IP)-(1./AL(IP))
2676           PB(4)=AL(IP)+(1./AL(IP))
2677           IJT=IJET(IP)
2678           DO 140 IHEP=IJT,IEND(IP)
2679           CALL HWULOF(PB,PHEP(1,IHEP),PHEP(1,IHEP))
2680           CALL HWULF4(PB,VHEP(1,IHEP),VHEP(1,IHEP))
2681   140     ISTHEP(IHEP)=ISTHEP(IHEP)+10
2682           CALL HWVDIF(4,VHEP(1,IPAR(IP)),VHEP(1,IJT),DISP)
2683           DO 150 IHEP=IJT,IEND(IP)
2684   150     CALL HWVSUM(4,DISP,VHEP(1,IHEP),VHEP(1,IHEP))
2685           IF (IEND(IP).GT.IJT+1) THEN
2686             ISTHEP(IJT+1)=100
2687           ELSEIF (IEND(IP).EQ.IJT) THEN
2688 C---NON-PARTON JET
2689             ISTHEP(IJT)=3
2690           ENDIF
2691   160     CONTINUE
2692         ENDIF
2693         ISTHEP(ICM)=120
2694       ELSE
2695 C---TIMELIKE JETS
2696 C---SPECIAL CASE: IF HARD PROCESS IS W/Z DECAY, PERFORM KINEMATIC
2697 C   RECONSTRUCTION IN ITS REST FRAME INSTEAD OF THE LAB FRAME
2698         IF (IDHW(ICM).GE.198.AND.IDHW(ICM).LE.200.AND.WZRFR) THEN
2699           CALL HWVEQU(5,PHEP(1,ICM),PLAB)
2700           CALL HWULOF(PLAB,PHEP(1,ICM),PHEP(1,ICM))
2701           CALL HWULF4(PLAB,VHEP(1,ICM),VHEP(1,ICM))
2702           DO 165 IP=1,NP
2703             CALL HWULOF(PLAB,PHEP(1,IPAR(IP)),PHEP(1,IPAR(IP)))
2704             CALL HWULF4(PLAB,VHEP(1,IPAR(IP)),VHEP(1,IPAR(IP)))
2705  165      CONTINUE
2706         ENDIF
2707 C   special for DIS: preserve outgoing lepton momentum
2708         IF (DISPRO) THEN
2709           CALL HWVEQU(5,PHEP(1,IPAR(1)),PHEP(1,IJET(1)))
2710           ISTHEP(IJET(1))=1
2711           LP=2
2712         ELSE
2713           CALL HWVEQU(5,PHEP(1,ICM),PC)
2714 C--- PQ AND PC ARE OLD AND NEW PARTON CM
2715           CALL HWVSUM(4,PHEP(1,IPAR(1)),PHEP(1,IPAR(2)),PQ)
2716           PQ(5)=PHEP(5,ICM)
2717           IF (NP.GT.2) THEN
2718             DO 170 KP=3,NP
2719   170       CALL HWVSUM(4,PHEP(1,IPAR(KP)),PQ,PQ)
2720           ENDIF
2721           LP=1
2722         ENDIF
2723         IF (.NOT.DISLOW) THEN
2724 C---FIND JET CM MOMENTA
2725           ECM=PQ(5)
2726           EMS=0.
2727           JETRAD=.FALSE.
2728           DO 180 KP=LP,NP
2729           EMJ=PHEP(5,IJET(KP))
2730           EMP=PHEP(5,IPAR(KP))
2731           JETRAD=JETRAD.OR.EMJ.NE.EMP
2732           EMS=EMS+EMJ
2733           PM(KP)= EMJ**2
2734 C---N.B. ROUNDING ERRORS HERE AT HIGH ENERGIES
2735           PJ(KP)=(HWULDO(PHEP(1,IPAR(KP)),PQ)/ECM)**2-EMP**2
2736           IF (PJ(KP).LE.ZERO) THEN
2737             CALL HWWARN('HWBJCO',104)
2738             GOTO 999
2739           ENDIF
2740   180     CONTINUE
2741           PF=1.
2742           IF (JETRAD) THEN
2743 C---JETS DID RADIATE
2744             IF (EMS.GE.ECM) THEN
2745               FROST=.TRUE.
2746               GOTO 240
2747             ENDIF
2748             DO 200 NE=1,NETRY
2749             EMS=-ECM
2750             DMS=0.
2751             DO 190 KP=LP,NP
2752             ES=SQRT(PF*PJ(KP)+PM(KP))
2753             EMS=EMS+ES
2754   190       DMS=DMS+PJ(KP)/ES
2755             DPF=2.*EMS/DMS
2756             IF (DPF.GT.PF) DPF=0.9*PF
2757             PF=PF-DPF
2758   200       IF (ABS(DPF).LT.EPS) GOTO 210
2759             CALL HWWARN('HWBJCO',105)
2760             GOTO 999
2761           ENDIF
2762   210     CONTINUE
2763         ENDIF
2764 C---BOOST PC AND PQ TO BREIT FRAME IF NECESSARY
2765         IF (DISPRO.AND.BREIT) THEN
2766           CALL HWULOF(PBR,PC,PC)
2767           CALL HWUROF(RBR,PC,PC)
2768           IF (.NOT.DISLOW) THEN
2769             CALL HWULOF(PBR,PQ,PQ)
2770             CALL HWUROF(RBR,PQ,PQ)
2771           ENDIF
2772         ENDIF
2773         DO 230 IP=LP,NP
2774 C---FIND CM ROTATION FOR JET IP
2775         IF (.NOT.DISLOW) THEN
2776           CALL HWVEQU(4,PHEP(1,IPAR(IP)),PR)
2777           IF (DISPRO.AND.BREIT) THEN
2778             CALL HWULOF(PBR,PR,PR)
2779             CALL HWUROF(RBR,PR,PR)
2780           ENDIF
2781 C--Modified by MHS 17/08/05 to do unboost in 2 stages (trans,long)
2782           PA(1)=PQ(1)
2783           PA(2)=PQ(2)
2784           PA(3)=ZERO
2785           PA(5)=SQRT(PQ(3)**2+PQ(5)**2)
2786           PA(4)=PQ(4)
2787           CALL HWULOF(PA,PR,PR)
2788           PA(1)=ZERO
2789           PA(2)=ZERO
2790           PA(3)=PQ(3)
2791           PA(4)=PA(5)
2792           PA(5)=PQ(5)
2793           CALL HWULOF(PA,PR,PR)
2794 C--End mod
2795           CALL HWUROT(PR, ONE,ZERO,RR)
2796           PR(1)=ZERO
2797           PR(2)=ZERO
2798           PR(3)=SQRT(PF*PJ(IP))
2799           PR(4)=SQRT(PF*PJ(IP)+PM(IP))
2800           PR(5)=PHEP(5,IJET(IP))
2801           CALL HWUROB(RR,PR,PR)
2802 C--Modified by BRW 25/10/02 to do boost in 2 stages (long,trans)
2803           PA(1)=ZERO
2804           PA(2)=ZERO
2805           PA(3)=PC(3)
2806           PA(5)=PC(5)
2807           PA(4)=SQRT(PA(3)**2+PA(5)**2)
2808           CALL HWULOB(PA,PR,PR)
2809           PA(1)=PC(1)
2810           PA(2)=PC(2)
2811           PA(3)=ZERO
2812           PA(5)=PA(4)
2813           PA(4)=PC(4)
2814           CALL HWULOB(PA,PR,PR)
2815 C--End mod
2816         ELSE
2817           CALL HWVEQU(5,PC,PR)
2818         ENDIF
2819 C---NOW PR IS LAB/BREIT MOMENTUM OF JET IP
2820         KP=IJET(IP)+1
2821         IF (AZCOR.AND.KP.LE.NHEP.AND.IDHW(KP).EQ.17) THEN
2822 C---ALIGN CONE WITH INTERFERING PARTON
2823           CALL HWUROT(PR, ONE,ZERO,RS)
2824           JP=JMOHEP(2,KP)
2825           IF (JP.EQ.0) THEN
2826             CALL HWWARN('HWBJCO',2)
2827             PTINF=0.
2828           ELSE
2829             CALL HWVEQU(4,PHEP(1,JP),PS)
2830             IF (DISPRO.AND.BREIT) THEN
2831               CALL HWULOF(PBR,PS,PS)
2832               CALL HWUROF(RBR,PS,PS)
2833             ENDIF
2834             CALL HWUROF(RS,PS,PS)
2835             PTINF=PS(1)**2+PS(2)**2
2836             IF (PTINF.LT.EPS) THEN
2837 C---COLLINEAR JETS: ALIGN CONES
2838               JP=JDAHEP(1,JP)+1
2839 C---BUG FIX BY MHS 17/03/05: RETURNED TO VERSION 6.500!
2840               IF (ISTHEP(JP).EQ.100.AND.ISTHEP(JP-1).GE.141
2841      $                             .AND.ISTHEP(JP-1).LE.144) THEN
2842 C---END FIX
2843                 CALL HWVEQU(4,PHEP(1,JP),PS)
2844                 IF (DISPRO.AND.BREIT) THEN
2845                   CALL HWULOF(PBR,PS,PS)
2846                   CALL HWUROF(RBR,PS,PS)
2847                 ENDIF
2848                 CALL HWUROF(RS,PS,PS)
2849                 PTINF=PS(1)**2+PS(2)**2
2850               ELSE
2851                 PTINF=0.
2852               ENDIF
2853             ENDIF
2854           ENDIF
2855           CALL HWVEQU(4,PHEP(1,KP),PB)
2856           IF (DISPRO.AND.BREIT) THEN
2857             CALL HWULOF(PBR,PB,PB)
2858             CALL HWUROF(RBR,PB,PB)
2859           ENDIF
2860           PTCON=PB(1)**2+PB(2)**2
2861           IF (PTCON.NE.ZERO.AND.PTINF.NE.ZERO) THEN
2862             CN=1./SQRT(PTINF*PTCON)
2863             CP=CN*(PS(1)*PB(1)+PS(2)*PB(2))
2864             SP=CN*(PS(1)*PB(2)-PS(2)*PB(1))
2865           ELSE
2866             CALL HWRAZM( ONE,CP,SP)
2867           ENDIF
2868         ELSE
2869           CALL HWRAZM( ONE,CP,SP)
2870         ENDIF
2871         CALL HWUROT(PR,CP,SP,RS)
2872 C---FIND BOOST FOR JET IP
2873         ALF=(PHEP(3,IJET(IP))+PHEP(4,IJET(IP)))/
2874      &      (PR(4)+SQRT((PR(4)+PR(5))*(PR(4)-PR(5))))
2875         PB(1)=0.
2876         PB(2)=0.
2877         PB(3)=ALF-(1./ALF)
2878         PB(4)=ALF+(1./ALF)
2879         PB(5)=2.
2880         IHEP=IJET(IP)
2881         KHEP=JDAHEP(2,IHEP)
2882         IF (KHEP.LT.IHEP) KHEP=IHEP
2883         DO 220 JHEP=IHEP,KHEP
2884         CALL HWULOF(PB,PHEP(1,JHEP),PHEP(1,JHEP))
2885         CALL HWUROB(RS,PHEP(1,JHEP),PHEP(1,JHEP))
2886         CALL HWULF4(PB,VHEP(1,JHEP),VHEP(1,JHEP))
2887         CALL HWUROB(RS,VHEP(1,JHEP),VHEP(1,JHEP))
2888 C---BOOST FROM BREIT FRAME IF NECESSARY
2889         IF (DISPRO.AND.BREIT) THEN
2890           CALL HWUROB(RBR,PHEP(1,JHEP),PHEP(1,JHEP))
2891           CALL HWULOB(PBR,PHEP(1,JHEP),PHEP(1,JHEP))
2892           CALL HWUROB(RBR,VHEP(1,JHEP),VHEP(1,JHEP))
2893           CALL HWULB4(PBR,VHEP(1,JHEP),VHEP(1,JHEP))
2894         ENDIF
2895         CALL HWVSUM(4,VHEP(1,JHEP),VHEP(1,IPAR(IP)),VHEP(1,JHEP))
2896 C--MHS FIX 07/03/05 FOR VERTEX POSITION OF LONG LIVED NON-PARTON JETS
2897         IF (KHEP.EQ.IHEP.AND.(IDHW(JHEP).GE.121.AND.IDHW(JHEP).LE.132
2898      $       .OR.IDHW(JHEP).EQ.59))
2899      $       CALL HWVSUM(4,VTXPIP,VHEP(1,JHEP),VHEP(1,JHEP))
2900 C--END FIX
2901   220   ISTHEP(JHEP)=ISTHEP(JHEP)+10
2902         IF (KHEP.GT.IHEP+1) THEN
2903           ISTHEP(IHEP+1)=100
2904         ELSEIF (KHEP.EQ.IHEP) THEN
2905 C---NON-PARTON JET
2906           ISTHEP(IHEP)=190
2907         ENDIF
2908   230   CONTINUE
2909         IF (ISTHEP(ICM).EQ.110) ISTHEP(ICM)=120
2910 C---SPECIAL CASE: FOR W/Z DECAY BOOST BACK TO THE LAB FRAME
2911  240    IF (IDHW(ICM).GE.198.AND.IDHW(ICM).LE.200.AND.WZRFR) THEN
2912           CALL HWULOB(PLAB,PHEP(1,ICM),PHEP(1,ICM))
2913           CALL HWULB4(PLAB,VHEP(1,ICM),VHEP(1,ICM))
2914           DO 260 IP=1,NP
2915             CALL HWULOB(PLAB,PHEP(1,IPAR(IP)),PHEP(1,IPAR(IP)))
2916             CALL HWULB4(PLAB,VHEP(1,IPAR(IP)),VHEP(1,IPAR(IP)))
2917             CALL HWULOB(PLAB,PHEP(1,IJET(IP)),PHEP(1,IJET(IP)))
2918 C--MHS FIX 07/03/05 - DO NOT REBOOST PRIMARY VERTEX
2919             IF (ISTHEP(IJET(IP)).EQ.190)
2920      $           CALL HWVDIF(4,VHEP(1,IJET(IP)),VTXPIP,VHEP(1,IJET(IP)))
2921             CALL HWULB4(PLAB,VHEP(1,IJET(IP)),VHEP(1,IJET(IP)))
2922             IF (ISTHEP(IJET(IP)).EQ.190)
2923      $           CALL HWVSUM(4,VHEP(1,IJET(IP)),VTXPIP,VHEP(1,IJET(IP)))
2924 C---END FIX
2925             IF (JDAHEP(1,IJET(IP)).GT.0) THEN
2926               IF (JDAHEP(2,IJET(IP)).GT.JDAHEP(1,IJET(IP))) THEN
2927                 CALL HWULOB(PLAB,PHEP(1,IJET(IP)+1),PHEP(1,IJET(IP)+1))
2928                 CALL HWULB4(PLAB,VHEP(1,IJET(IP)+1),VHEP(1,IJET(IP)+1))
2929               ENDIF
2930               DO 250 IHEP=JDAHEP(1,IJET(IP)),JDAHEP(2,IJET(IP))
2931                 CALL HWULOB(PLAB,PHEP(1,IHEP),PHEP(1,IHEP))
2932                 CALL HWULB4(PLAB,VHEP(1,IHEP),VHEP(1,IHEP))
2933  250          CONTINUE
2934             ENDIF
2935  260      CONTINUE
2936         ENDIF
2937         IF (FROST) RETURN
2938       ENDIF
2939       GOTO 20
2940  999  RETURN
2941       END
2942 CDECK  ID>, HWBMAS.
2943 *CMZ :-        -26/04/91  11.11.54  by  Bryan Webber
2944 *-- Author :    Bryan Webber
2945 C-----------------------------------------------------------------------
2946       SUBROUTINE HWBMAS
2947 C-----------------------------------------------------------------------
2948 C     Passes  backwards through a  jet cascade  calculating the masses
2949 C     and magnitudes of the longitudinal and transverse three momenta.
2950 C     Components given relative to direction of parent for a time-like
2951 C     vertex and with respect to z-axis for space-like vertices.
2952 C
2953 C     On input PPAR(1-5,*) contains:
2954 C     (E*sqrt(Xi),Xi,3-mom (if external),E,M-sq (if external))
2955 C
2956 C     On output PPAR(1-5,*) (if TMPAR(*)), containts:
2957 C     (P-trans,Xi or Xilast,P-long,E,M)
2958 C-----------------------------------------------------------------------
2959       INCLUDE 'HERWIG65.INC'
2960       DOUBLE PRECISION HWUSQR,EXI,PISQ,PJPK,EJEK,PTSQ,Z,ZMIN,ZMAX,
2961      $     EMI,EMJ,EMK,C,NQ,HWBVMC,RHO,POLD,PNEW,EOLD,ENEW,A,B
2962       INTEGER IPAR,JPAR,KPAR,MPAR,I,J,K
2963       EXTERNAL HWUSQR
2964       IF (IERROR.NE.0) RETURN
2965       IF (NPAR.GT.2) THEN
2966         DO 30 MPAR=NPAR-1,3,-2
2967          JPAR=MPAR
2968 C Find parent and partner of this branch
2969          IPAR=JMOPAR(1,JPAR)
2970          KPAR=JPAR+1
2971 C Determine type of branching
2972          IF (TMPAR(IPAR)) THEN
2973 C Time-like branching
2974 C           Compute mass of parent
2975             EXI=PPAR(1,JPAR)*PPAR(1,KPAR)
2976             PPAR(5,IPAR)=PPAR(5,JPAR)+PPAR(5,KPAR)+2.*EXI
2977 C           Compute three momentum of parent
2978             PISQ=PPAR(4,IPAR)*PPAR(4,IPAR)-PPAR(5,IPAR)
2979             PPAR(3,IPAR)=HWUSQR(PISQ)
2980 C---SPECIAL FOR G-->QQBAR: READJUST ANGULAR DISTRIBUTION
2981             IF (IDPAR(IPAR).EQ.13 .AND. IDPAR(JPAR).LT.13) THEN
2982               Z=PPAR(4,JPAR)/PPAR(4,IPAR)
2983               ZMIN=HWBVMC(IDPAR(JPAR))/PPAR(1,JPAR)*Z
2984               RHO=(Z*(3-Z*(3-2*Z))-ZMIN*(3-ZMIN*(3-2*ZMIN)))
2985      $             /(2*(1-2*ZMIN)*(1-ZMIN*(1-ZMIN)))
2986               NQ=PPAR(3,IPAR)*(PPAR(3,IPAR)+PPAR(4,IPAR))
2987               EMI=PPAR(5,IPAR)
2988               EMJ=PPAR(5,JPAR)
2989               EMK=PPAR(5,KPAR)
2990               ZMIN=MAX((EMI+EMJ-EMK)/(2*(EMI+NQ)),
2991      $      (EMI+EMJ-EMK-SQRT(ABS((EMI-EMJ-EMK)**2-4*EMJ*EMK)))/(2*EMI))
2992               ZMAX=1-MAX((EMI-EMJ+EMK)/(2*(EMI+NQ)),
2993      $      (EMI-EMJ+EMK-SQRT(ABS((EMI-EMJ-EMK)**2-4*EMJ*EMK)))/(2*EMI))
2994               C=2*RMASS(IDPAR(JPAR))**2/EMI
2995               Z=(4*ZMIN*(1.5*(1+C-ZMIN)+ZMIN**2)*(1-RHO)
2996      $          +4*ZMAX*(1.5*(1+C-ZMAX)+ZMAX**2)*RHO-2-3*C)/(1+2*C)**1.5
2997               Z=SQRT(1+2*C)*SINH(LOG(Z+SQRT(Z**2+1))/3)+0.5
2998               Z=(Z*NQ+(EMI+EMJ-EMK)/2)/(NQ+EMI)
2999               PPAR(4,JPAR)=Z*PPAR(4,IPAR)
3000               PPAR(4,KPAR)=PPAR(4,IPAR)-PPAR(4,JPAR)
3001               PPAR(3,JPAR)=HWUSQR(PPAR(4,JPAR)**2-EMJ)
3002               PPAR(3,KPAR)=HWUSQR(PPAR(4,KPAR)**2-EMK)
3003               PPAR(2,JPAR)=EXI/(PPAR(4,JPAR)*PPAR(4,KPAR))
3004               IF(JDAPAR(2,JPAR).NE.0)PPAR(2,JDAPAR(2,JPAR))=PPAR(2,JPAR)
3005               IF(JDAPAR(2,KPAR).NE.0)PPAR(2,JDAPAR(2,KPAR))=PPAR(2,JPAR)
3006 C---FIND DESCENDENTS OF THIS SPLITTING AND READJUST THEIR MOMENTA TOO
3007               DO 20 J=JPAR+2,NPAR-1,2
3008                 I=J
3009  10             I=JMOPAR(1,I)
3010                 IF (I.GT.IPAR) GOTO 10
3011                 IF (I.EQ.IPAR) THEN
3012                   I=JMOPAR(1,J)
3013                   K=J+1
3014                   POLD=PPAR(3,J)+PPAR(3,K)
3015                   EOLD=PPAR(4,J)+PPAR(4,K)
3016                   PNEW=HWUSQR(PPAR(4,I)**2-PPAR(5,I))
3017                   ENEW=PPAR(4,I)
3018                   A=(ENEW*EOLD-PNEW*POLD)/PPAR(5,I)
3019                   B=(PNEW*EOLD-ENEW*POLD)/PPAR(5,I)
3020                   PPAR(3,J)=A*PPAR(3,J)+B*PPAR(4,J)
3021                   PPAR(4,J)=(PPAR(4,J)+B*PPAR(3,J))/A
3022                   PPAR(3,K)=PNEW-PPAR(3,J)
3023                   PPAR(4,K)=ENEW-PPAR(4,J)
3024                   PPAR(2,J)=1-(PPAR(3,J)*PPAR(3,K)+PPAR(1,J)*PPAR(1,K))
3025      $                 /(PPAR(4,J)*PPAR(4,K))
3026                   IF (JDAPAR(2,J).NE.0) PPAR(2,JDAPAR(2,J))=PPAR(2,J)
3027                   IF (JDAPAR(2,K).NE.0) PPAR(2,JDAPAR(2,K))=PPAR(2,J)
3028                 ENDIF
3029  20           CONTINUE
3030             ENDIF
3031 C           Compute daughter' transverse and longitudinal momenta
3032             PJPK=PPAR(3,JPAR)*PPAR(3,KPAR)
3033             EJEK=PPAR(4,JPAR)*PPAR(4,KPAR)-EXI
3034             PTSQ=(PJPK+EJEK)*(PJPK-EJEK)/PISQ
3035             PPAR(1,JPAR)=HWUSQR(PTSQ)
3036             PPAR(3,JPAR)=HWUSQR(PPAR(3,JPAR)*PPAR(3,JPAR)-PTSQ)
3037             PPAR(1,KPAR)=-PPAR(1,JPAR)
3038             PPAR(3,KPAR)= PPAR(3,IPAR)-PPAR(3,JPAR)
3039          ELSE
3040 C Space-like branching
3041 C           Re-arrange such that JPAR is time-like
3042             IF (TMPAR(KPAR)) THEN
3043                KPAR=JPAR
3044                JPAR=JPAR+1
3045             ENDIF
3046 C           Compute time-like branch
3047             PTSQ=(2.-PPAR(2,JPAR))*PPAR(1,JPAR)*PPAR(1,JPAR)
3048      &          -PPAR(5,JPAR)
3049             PPAR(1,JPAR)=HWUSQR(PTSQ)
3050             PPAR(3,JPAR)=(1.-PPAR(2,JPAR))*PPAR(4,JPAR)
3051             PPAR(3,IPAR)=PPAR(3,KPAR)-PPAR(3,JPAR)
3052             PPAR(5,IPAR)=0.
3053             PPAR(1,KPAR)=0.
3054          ENDIF
3055 C Reset Xi to Xilast
3056          PPAR(2,KPAR)=PPAR(2,IPAR)
3057  30    CONTINUE
3058       ENDIF
3059       DO 40 IPAR=2,NPAR
3060  40   PPAR(5,IPAR)=HWUSQR(PPAR(5,IPAR))
3061       PPAR(1,2)=0.
3062       PPAR(2,2)=0.
3063       END
3064 CDECK  ID>, HWBRAN.
3065 *CMZ :-        -14/10/99  18.04.56  by  Mike Seymour
3066 *-- Author :    Bryan Webber & Mike Seymour
3067 C-----------------------------------------------------------------------
3068       SUBROUTINE HWBRAN(KPAR)
3069 C-----------------------------------------------------------------------
3070 C     BRANCHES TIMELIKE PARTON KPAR INTO TWO, PUTS PRODUCTS
3071 C     INTO NPAR+1 AND NPAR+2, AND INCREASES NPAR BY TWO
3072 C-----------------------------------------------------------------------
3073       INCLUDE 'HERWIG65.INC'
3074       DOUBLE PRECISION HWBVMC,HWRGEN,HWUALF,HWUTAB,HWRUNI,HWULDO,PMOM,
3075      & QNOW,QLST,QKTHR,RN,QQBAR,DQQ,QGTHR,SNOW,QSUD,ZMIN,ZMAX,ZRAT,WMIN,
3076      & QLAM,Z1,Z2,ETEST,ZTEST,ENOW,XI,XIPREV,EPREV,QMAX,QGAM,SLST,SFNL,
3077      & TARG,ALF,BETA0(3:6),BETAP(3:6),SQRK(4:6,5),REJFAC,Z,X1,X2,OTHXI,
3078      & OTHZ,X3,FF,AW,XCUT,CC,JJ,HWUSQR
3079       INTEGER HWRINT,KPAR,ID,JD,IS,NTRY,N,ID1,ID2,MPAR,ISUD(13),IHEP,
3080      & JHEP,M,NF,NN,IREJ,NREJ,ITOP
3081       EXTERNAL HWBVMC,HWRGEN,HWUALF,HWUTAB,HWRUNI,HWULDO,HWRINT,HWUSQR
3082       SAVE BETA0,BETAP,SQRK
3083       SAVE ISUD
3084       DATA ISUD,BETA0/2,2,3,4,5,6,2,2,3,4,5,6,1,4*ZERO/
3085       IF (IERROR.NE.0) RETURN
3086 C---SET SQRK(M,N) TO THE PROBABILITY THAT A GLUON WILL NOT PRODUCE A
3087 C   QUARK-ANTIQUARK PAIR BETWEEN SCALES RMASS(M) AND 2*HWBVMC(N)
3088       IF (SUDORD.NE.1.AND.BETA0(3).EQ.ZERO) THEN
3089         DO 100 M=3,6
3090           BETA0(M)=(11.*CAFAC-2.*M)*0.5
3091  100      BETAP(M)=(17.*CAFAC**2-(5.*CAFAC+3.*CFFAC)*M)
3092      &            /BETA0(M)*0.25/PIFAC
3093         DO 120 N=1,5
3094           DO 110 M=4,6
3095             IF (M.LE.N) THEN
3096               SQRK(M,N)=ONE
3097             ELSEIF (M.EQ.4.OR.M.EQ.N+1) THEN
3098               NF=M
3099               IF (2*HWBVMC(N).GT.RMASS(M)) NF=M+1
3100               SQRK(M,N)=((BETAP(NF-1)+1/HWUALF(1,2*HWBVMC(N)))/
3101      $             (BETAP(NF-1)+1/HWUALF(1,RMASS(M))))**(1/BETA0(NF-1))
3102             ELSE
3103               SQRK(M,N)=SQRK(M-1,N)*
3104      $             ((BETAP(M-1)+1/HWUALF(1,RMASS(M-1)))/
3105      $             (BETAP(M-1)+1/HWUALF(1,RMASS(M))))**(1/BETA0(M-1))
3106             ENDIF
3107  110      CONTINUE
3108  120    CONTINUE
3109       ENDIF
3110       ID=IDPAR(KPAR)
3111 C--TEST FOR PARTON TYPE
3112       IF (ID.LE.13) THEN
3113         JD=ID
3114         IS=ISUD(ID)
3115       ELSEIF (ID.GE.209.AND.ID.LE.220) THEN
3116         JD=ID-208
3117         IS=7
3118       ELSE
3119         IS=0
3120       END IF
3121       QNOW=-1.
3122       IF (IS.NE.0) THEN
3123 C--TIMELIKE PARTON BRANCHING
3124         ENOW=PPAR(4,KPAR)
3125         XIPREV=PPAR(2,KPAR)
3126         IF (JMOPAR(1,KPAR).EQ.0) THEN
3127           EPREV=PPAR(4,KPAR)
3128         ELSE
3129           EPREV=PPAR(4,JMOPAR(1,KPAR))
3130         ENDIF
3131 C--IF THIS IS CHARGED & PHOTONS ARE ALLOWED, ANGLES MIGHT NOT BE ORDERED
3132         QMAX=0
3133         QLST=PPAR(1,KPAR)
3134         IF (ICHRG(ID).NE.0 .AND. VPCUT.LT.PPAR(1,2)) THEN
3135 C--LOOK FOR A PREVIOUS G->QQBAR, IF ANY
3136           MPAR=KPAR
3137  1        IF (JMOPAR(1,MPAR).NE.0) THEN
3138             IF (IDPAR(JMOPAR(1,MPAR)).EQ.ID) THEN
3139               MPAR=JMOPAR(1,MPAR)
3140               GOTO 1
3141             ENDIF
3142           ENDIF
3143 C--IF CLIMBED TO THE TOP OF THE LIST, FIND QED INTERFERENCE PARTNER
3144           IF (MPAR.EQ.2) THEN
3145             JHEP=0
3146             IF (ID.LT.7) THEN
3147               IHEP=JDAHEP(2,JCOPAR(1,1))
3148               IF (IHEP.GT.0) JHEP=JDAHEP(2,IHEP)
3149             ELSE
3150               IHEP=JMOHEP(2,JCOPAR(1,1))
3151               IF (IHEP.GT.0) JHEP=JMOHEP(2,IHEP)
3152             ENDIF
3153             IF (IHEP.GT.0.AND.JHEP.GT.0) THEN
3154                QMAX=HWULDO(PHEP(1,IHEP),PHEP(1,JHEP))
3155      &              *(ENOW/PPAR(4,2))**2
3156             ELSE
3157 C--FIX AT HARD PROCESS SCALE IF POINTER NOT YET SET
3158 C  (CAN HAPPEN IN SUSY EVENTS)
3159                QMAX=EMSCA**2
3160             ENDIF
3161           ELSE
3162             QMAX=ENOW**2*PPAR(2,MPAR)
3163           ENDIF
3164 C--IF PREVIOUS BRANCHING WAS Q->QGAMMA, LOOK FOR A QCD BRANCHING
3165           MPAR=KPAR
3166  2        IF (JMOPAR(1,MPAR).NE.0) THEN
3167             IF (IDPAR(JDAPAR(1,JMOPAR(1,MPAR))).EQ.59 .OR.
3168      &        IDPAR(JDAPAR(2,JMOPAR(1,MPAR))).EQ.59) THEN
3169               MPAR=JMOPAR(1,MPAR)
3170               GOTO 2
3171             ENDIF
3172           ENDIF
3173           QLST=ENOW**2*PPAR(2,MPAR)
3174           QMAX=SQRT(MAX(ZERO,MIN(
3175      &         QMAX , EPREV**2*XIPREV , ENOW**2*XIPREV*(2-XIPREV))))
3176           QLST=SQRT(MIN(
3177      &         QLST , EPREV**2*XIPREV , ENOW**2*XIPREV*(2-XIPREV)))
3178         ENDIF
3179         NTRY=0
3180     5   NTRY=NTRY+1
3181         IF (NTRY.GT.NBTRY) THEN
3182           CALL HWWARN('HWBRAN',100)
3183           GOTO 999
3184         ENDIF
3185         IF (ID.EQ.13) THEN
3186 C--GLUON -> QUARK+ANTIQUARK OPTION
3187           IF (QLST.GT.QCDL3) THEN
3188             DO 8 N=1,NFLAV
3189             QKTHR=2.*HWBVMC(N)
3190             IF (QLST.GT.QKTHR) THEN
3191               RN=HWRGEN(N)
3192               IF (SUDORD.NE.1) THEN
3193 C---FIND IN WHICH FLAVOUR INTERVAL THE UPPER LIMIT LIES
3194                 NF=3
3195                 DO 200 M=MAX(3,N),NFLAV
3196  200              IF (QLST.GT.RMASS(M)) NF=M
3197 C---CALCULATE THE FORM FACTOR
3198                 IF (NF.EQ.MAX(3,N)) THEN
3199                   SFNL=((BETAP(NF)+1/HWUALF(1,QKTHR))/
3200      $                 (BETAP(NF)+1/HWUALF(1,QLST)))**(1/BETA0(NF))
3201                   SLST=SFNL
3202                 ELSE
3203                   SFNL=((BETAP(NF)+1/HWUALF(1,RMASS(NF)))/
3204      $                 (BETAP(NF)+1/HWUALF(1,QLST)))**(1/BETA0(NF))
3205                   SLST=SFNL*SQRK(NF,N)
3206                 ENDIF
3207               ENDIF
3208               IF (RN.GT.1.E-3) THEN
3209                 QQBAR=QCDL3*(QLST/QCDL3)**(RN**BETAF)
3210               ELSE
3211                 QQBAR=QCDL3
3212               ENDIF
3213               IF (SUDORD.NE.1) THEN
3214 C---FIND IN WHICH FLAVOUR INTERVAL THE SOLUTION LIES
3215                 IF (RN.GE.SFNL) THEN
3216                   NN=NF
3217                 ELSEIF (RN.GE.SLST) THEN
3218                   NN=MAX(3,N)
3219                   DO 210 M=MAX(3,N)+1,NF-1
3220  210                IF (RN.GE.SLST/SQRK(M,N)) NN=M
3221                 ELSE
3222                   NN=0
3223                   QQBAR=QCDL3
3224                 ENDIF
3225                 IF (NN.GT.0) THEN
3226                   IF (NN.EQ.NF) THEN
3227                     TARG=HWUALF(1,QLST)
3228                   ELSE
3229                     TARG=HWUALF(1,RMASS(NN+1))
3230                     RN=RN/SLST*SQRK(NN+1,N)
3231                   ENDIF
3232                   TARG=1/((BETAP(NN)+1/TARG)*RN**BETA0(NN)-BETAP(NN))
3233 C---NOW SOLVE HWUALF(1,QQBAR)=TARG FOR QQBAR ITERATIVELY
3234  7                QQBAR=MAX(QQBAR,HALF*QKTHR)
3235                   ALF=HWUALF(1,QQBAR)
3236                   IF (ABS(ALF-TARG).GT.ACCUR) THEN
3237                     NTRY=NTRY+1
3238                     IF (NTRY.GT.NBTRY) THEN
3239                       CALL HWWARN('HWBRAN',101)
3240                       GOTO 999
3241                     ENDIF
3242                     QQBAR=QQBAR*(1+3*PIFAC*(ALF-TARG)
3243      $                   /(BETA0(NN)*ALF**2*(1+BETAP(NN)*ALF)))
3244                     GOTO 7
3245                   ENDIF
3246                 ENDIF
3247               ENDIF
3248               IF (QQBAR.GT.QNOW.AND.QQBAR.GT.QKTHR) THEN
3249                 QNOW=QQBAR
3250                 ID2=N
3251               ENDIF
3252             ELSE
3253               GOTO 9
3254             ENDIF
3255     8       CONTINUE
3256           ENDIF
3257 C--GLUON->DIQUARKS OPTION
3258     9     IF (QLST.LT.QDIQK) THEN
3259             IF (PDIQK.NE.ZERO) THEN
3260               RN=HWRGEN(0)
3261               DQQ=QLST*EXP(-RN/PDIQK)
3262               IF (DQQ.GT.QNOW) THEN
3263                 IF (DQQ.GT.2.*RMASS(115)) THEN
3264                   QNOW=DQQ
3265                   ID2=115
3266                 ENDIF
3267               ENDIF
3268             ENDIF
3269           ENDIF
3270         ENDIF
3271 C--ENHANCE GLUON AND PHOTON EMISSION BY A FACTOR OF TWO IF THIS BRANCH
3272 C  IS CAPABLE OF BEING THE HARDEST SO FAR
3273         NREJ=1
3274         IF (TMPAR(2).AND.0.25*MAX(QLST,QMAX).GT.HARDST) NREJ=2
3275 C--BRANCHING ID->ID+GLUON
3276         QGTHR=HWBVMC(ID)+HWBVMC(13)
3277         IF (QLST.GT.QGTHR) THEN
3278          DO 300 IREJ=1,NREJ
3279           RN=HWRGEN(1)
3280           SLST=HWUTAB(SUD(1,IS),QEV(1,IS),NQEV,QLST,INTER)
3281           IF (RN.EQ.ZERO) THEN
3282             SNOW=2.
3283           ELSE
3284             SNOW=SLST/RN
3285           ENDIF
3286           IF (SNOW.LT.ONE) THEN
3287             QSUD=HWUTAB(QEV(1,IS),SUD(1,IS),NQEV,SNOW,INTER)
3288 C---IF FORM FACTOR DID NOT GET INVERTED CORRECTLY TRY LINEAR INSTEAD
3289             IF (QSUD.GT.QLST) THEN
3290               SNOW=HWUTAB(SUD(1,IS),QEV(1,IS),NQEV,QLST,1)/RN
3291               QSUD=HWUTAB(QEV(1,IS),SUD(1,IS),NQEV,SNOW,1)
3292               IF (QSUD.GT.QLST) THEN
3293                 CALL HWWARN('HWBRAN',1)
3294                 QSUD=-1
3295               ENDIF
3296             ENDIF
3297             IF (QSUD.GT.QGTHR.AND.QSUD.GT.QNOW) THEN
3298               ID2=13
3299               QNOW=QSUD
3300             ENDIF
3301           ENDIF
3302  300     CONTINUE
3303         ENDIF
3304 C--BRANCHING ID->ID+PHOTON
3305         IF (ICHRG(ID).NE.0) THEN
3306           QGTHR=MAX(HWBVMC(ID)+HWBVMC(59),HWBVMC(59)*EXP(0.75))
3307           IF (QMAX.GT.QGTHR) THEN
3308            DO 400 IREJ=1,NREJ
3309             RN=HWRGEN(2)
3310             IF (RN.EQ.ZERO) THEN
3311               QGAM=0
3312             ELSE
3313               QGAM=(LOG(QMAX/HWBVMC(59))-0.75)**2
3314      &            +PIFAC*9/(ICHRG(ID)**2*ALPFAC*ALPHEM)*LOG(RN)
3315               IF (QGAM.GT.ZERO) THEN
3316                 QGAM=HWBVMC(59)*EXP(0.75+SQRT(QGAM))
3317               ELSE
3318                 QGAM=0
3319               ENDIF
3320             ENDIF
3321             IF (QGAM.GT.QGTHR.AND.QGAM.GT.QNOW) THEN
3322               ID2=59
3323               QNOW=QGAM
3324             ENDIF
3325  400       CONTINUE
3326           ENDIF
3327         ENDIF
3328         IF (QNOW.GT.ZERO) THEN
3329 C--BRANCHING HAS OCCURRED
3330           ZMIN=HWBVMC(ID2)/QNOW
3331           ZMAX=1.-ZMIN
3332           IF (ID.EQ.13) THEN
3333             IF (ID2.EQ.13) THEN
3334 C--GLUON -> GLUON + GLUON
3335               ID1=13
3336               WMIN=ZMIN*ZMAX
3337               ETEST=(1.-WMIN)**2*HWUALF(5-SUDORD*2,QNOW*WMIN)
3338               ZRAT=(ZMAX*(1-ZMIN))/(ZMIN*(1-ZMAX))
3339 C--CHOOSE Z1 DISTRIBUTED ON (ZMIN,ZMAX)
3340 C  ACCORDING TO GLUON BRANCHING FUNCTION
3341    10         Z1=ZMAX/(ZMAX+(1-ZMAX)*ZRAT**HWRGEN(0))
3342               Z2=1.-Z1
3343               ZTEST=(1.-(Z1*Z2))**2*HWUALF(5-SUDORD*2,QNOW*(Z1*Z2))
3344               IF (ZTEST.LT.ETEST*HWRGEN(1)) GOTO 10
3345               Z=Z1
3346             ELSEIF (ID2.NE.115) THEN
3347 C--GLUON -> QUARKS
3348               ID1=ID2+6
3349               ETEST=ZMIN**2+ZMAX**2
3350    20         Z1=HWRUNI(0,ZMIN,ZMAX)
3351               Z2=1.-Z1
3352               ZTEST=Z1*Z1+Z2*Z2
3353               IF (ZTEST.LT.ETEST*HWRGEN(0)) GOTO 20
3354             ELSE
3355 C--GLUON -> DIQUARKS
3356               ID2=HWRINT(115,117)
3357               ID1=ID2-6
3358               Z1=HWRUNI(0,ZMIN,ZMAX)
3359               Z2=1.-Z1
3360             ENDIF
3361           ELSE
3362 C--QUARK OR ANTIQUARK BRANCHING
3363             IF (ID2.EQ.13) THEN
3364 C--TO GLUON
3365               ZMAX=1.-HWBVMC(ID)/QNOW
3366               WMIN=MIN(ZMIN*(1.-ZMIN),ZMAX*(1.-ZMAX))
3367               ETEST=(1.+ZMAX**2)*HWUALF(5-SUDORD*2,QNOW*WMIN)
3368               ZRAT=ZMAX/ZMIN
3369    30         Z1=ZMIN*ZRAT**HWRGEN(0)
3370               Z2=1.-Z1
3371               ZTEST=(1.+Z2*Z2)*HWUALF(5-SUDORD*2,QNOW*Z1*Z2)
3372               IF (ZTEST.LT.ETEST*HWRGEN(1)) GOTO 30
3373             ELSE
3374 C--TO PHOTON
3375               ZMIN=  HWBVMC(59)/QNOW
3376               ZMAX=1-HWBVMC(ID)/QNOW
3377               ZRAT=ZMAX/ZMIN
3378               ETEST=1+(1-ZMIN)**2
3379    40         Z1=ZMIN*ZRAT**HWRGEN(0)
3380               Z2=1-Z1
3381               ZTEST=1+Z2*Z2
3382               IF (ZTEST.LT.ETEST*HWRGEN(1)) GOTO 40
3383             ENDIF
3384 C--QUARKS EMIT ON LOWER SIDE, ANTIQUARKS ON UPPER SIDE
3385             Z=Z1
3386             IF (JD.LE.6) THEN
3387               Z1=Z2
3388               Z2=1.-Z2
3389               ID1=ID
3390             ELSE
3391               ID1=ID2
3392               ID2=ID
3393             ENDIF
3394           ENDIF
3395 C--UPDATE THIS BRANCH AND CREATE NEW BRANCHES
3396           XI=(QNOW/ENOW)**2
3397           IF (ID1.NE.59.AND.ID2.NE.59) THEN
3398             IF (ID.EQ.13.AND.ID1.NE.13) THEN
3399               QLAM=QNOW
3400             ELSE
3401               QLAM=QNOW*Z1*Z2
3402             ENDIF
3403             IF (SUDORD.EQ.1.AND.HWUALF(2,QLAM).LT.HWRGEN(0) .OR.
3404      &           (2.-XI)*(QNOW*Z1*Z2)**2.GT.EMSCA**2) THEN
3405 C--BRANCHING REJECTED: REDUCE Q AND REPEAT
3406                 QMAX=QNOW
3407                 QLST=QNOW
3408                 QNOW=-1.
3409                 GOTO 5
3410             ENDIF
3411           ENDIF
3412 C--IF THIS IS HARDEST EMISSION SO FAR, APPLY MATRIX-ELEMENT CORRECTION
3413           IF (ID.NE.13.OR.ID1.EQ.13) THEN
3414             QLAM=QNOW*Z1*Z2
3415             REJFAC=1
3416             IF (TMPAR(2).AND.QLAM.GT.HARDST) THEN
3417 C----SOFT MATRIX-ELEMENT CORRECTION TO TOP DECAYS
3418               ITOP=JCOPAR(1,1)
3419               IF (ISTHEP(ITOP).EQ.155.AND.(IDHW(ITOP).EQ.6
3420      $             .OR.IDHW(ITOP).EQ.12)) THEN
3421                 AW=(PHEP(5,JDAHEP(1,ITOP))/PHEP(5,ITOP))**2
3422                 FF=0.5*(1-AW)*(1-2*AW+1/AW)