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