Update to 6.5-10
[u/mrichter/AliRoot.git] / HERWIG / herwig6507.f
CommitLineData
31d78ebd 1C-----------------------------------------------------------------------
2C H E R W I G
3C
4C a Monte Carlo event generator for simulating
5C +---------------------------------------------------+
6C | Hadron Emission Reactions With Interfering Gluons |
7C +---------------------------------------------------+
8C I.G. Knowles(*), G. Marchesini(+), M.H.Seymour($,&) and B.R. Webber(#)
9C-----------------------------------------------------------------------
10C with Minimal Supersymmetric Standard Model Matrix Elements by
11C S. Moretti(") and K. Odagiri(^)
12C-----------------------------------------------------------------------
13C R parity violating Supersymmetric Decays and Matrix Elements by
14C P. Richardson(X)
15C-----------------------------------------------------------------------
16C matrix element corrections to top decay and Drell-Yan type processes
17C by G. Corcella(&)
18C-----------------------------------------------------------------------
19C Deep Inelastic Scattering and Heavy Flavour Electroproduction by
20C G. Abbiendi(@) and L. Stanco(%)
21C-----------------------------------------------------------------------
22C and Jet Photoproduction in Lepton-Hadron Collisions by J. Chyla(~)
23C-----------------------------------------------------------------------
24C(*) Department of Physics & Astronomy, University of Edinburgh
25C(+) Dipartimento di Fisica, Universita di Milano-Bicocca
26C($) Department of Physics & Astronomy, University of Manchester
27C(&) Theory Division, CERN
28C(#) Cavendish Laboratory, Cambridge
29C(") School of Physics & Astronomy, Southampton
30C(^) Academia Sinica, Taiwan
31C(X) Institute of Particle Physics Phenomenology, University of Durham
32C(@) Dipartimento di Fisica, Universita di Bologna
33C(%) Dipartimento di Fisica, Universita di Padova
34C(~) Institute of Physics, Prague
35C-----------------------------------------------------------------------
36C Version 6.507 - 8th March 2005
37C-----------------------------------------------------------------------
38C Main references:
39C
40C G.Corcella, I.G.Knowles, G.Marchesini, S.Moretti, K.Odagiri,
41C P.Richardson, M.H.Seymour and B.R.Webber, JHEP 0101 (2001) 010
42C
43C G.Marchesini, B.R.Webber, G.Abbiendi, I.G.Knowles, M.H.Seymour,
44C and L.Stanco, Computer Physics Communications 67 (1992) 465.
45C-----------------------------------------------------------------------
46C Please see the official HERWIG information page:
47C http://hepwww.rl.ac.uk/theory/seymour/herwig/
48C-----------------------------------------------------------------------
49CDECK ID>, CIRCEE.
50*CMZ :- -03/07/01 17.07.47 by Bryan Webber
51*-- Author : Bryan Webber
52C-----------------------------------------------------------------------
53 FUNCTION CIRCEE (X1, X2)
54C-----------------------------------------------------------------------
55C DUMMY FUNCTION: DELETE AND SET CIRCOP NON-ZERO
56C IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
57C-----------------------------------------------------------------------
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
64CDECK ID>, CIRCES.
65*CMZ :- -03/07/01 17.07.47 by Bryan Webber
66*-- Author : Bryan Webber
67C-----------------------------------------------------------------------
68 SUBROUTINE CIRCES (XX1M, XX2M, XROOTS, XACC, XVER, XREV, XCHAT)
69C-----------------------------------------------------------------------
70C DUMMY SUBROUTINE: DELETE AND SET CIRCOP NON-ZERO
71C IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
72C-----------------------------------------------------------------------
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
79CDECK ID>, CIRCGG.
80*CMZ :- -03/07/01 17.07.47 by Bryan Webber
81*-- Author : Bryan Webber
82C-----------------------------------------------------------------------
83 FUNCTION CIRCGG (X1, X2)
84C-----------------------------------------------------------------------
85C DUMMY FUNCTION: DELETE AND SET CIRCOP NON-ZERO
86C IN MAIN PROGRAM IF YOU USE CIRCE BEAM SPECTRUM PACKAGE
87C-----------------------------------------------------------------------
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
94CDECK ID>, DECADD.
95*CMZ :- -28/01/92 12.34.44 by Mike Seymour
96*-- Author : Luca Stanco
97C-----------------------------------------------------------------------
98 SUBROUTINE DECADD(LOGI)
99C-----------------------------------------------------------------------
100C DUMMY SUBROUTINE: DELETE AND SET BDECAY='CLEO'
101C IN MAIN PROGRAM IF YOU USE CLEO DECAY PACKAGE
102C-----------------------------------------------------------------------
103 LOGICAL LOGI
104 WRITE (6,10)
105 10 FORMAT(/10X,'DECADD CALLED BUT NOT LINKED')
106 STOP
107 END
108CDECK ID>, DEXAY.
109*CMZ :- -17/10/01 10.03.37 by Peter Richardson
110*-- Author : Peter Richardson
111C-----------------------------------------------------------------------
112 SUBROUTINE DEXAY(IMODE,POL)
113C-----------------------------------------------------------------------
114C DUMMY SUBROUTINE: DELETE AND SET TAUDEC='TAUOLA'
115C IN MAIN PROGRAM IF YOU USE TAUOLA DECAY PACKAGE
116C-----------------------------------------------------------------------
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
124CDECK ID>, EUDINI.
125*CMZ :- -28/01/92 12.34.44 by Mike Seymour
126*-- Author : Luca Stanco
127C-----------------------------------------------------------------------
128 SUBROUTINE EUDINI
129C-----------------------------------------------------------------------
130C DUMMY SUBROUTINE: DELETE AND SET BDECAY='EURO'
131C IN MAIN PROGRAM IF YOU USE EURODEC DECAY PACKAGE
132C-----------------------------------------------------------------------
133 WRITE (6,10)
134 10 FORMAT(/10X,'EUDINI CALLED BUT NOT LINKED')
135 STOP
136 END
137CDECK ID>, FILHEP.
138*CMZ :- -17/10/01 09:42:21 by Peter Richardson
139*-- Author : Martin W. Gruenewald
140C-----------------------------------------------------------------------
141 SUBROUTINE FILHEP(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
142C ----------------------------------------------------------------------
143C this subroutine fills one entry into the HEPEVT common
144C and updates the information for affected mother entries
145C used by TAUOLA
146C
147C written by Martin W. Gruenewald (91/01/28)
148C ----------------------------------------------------------------------
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)
156C
157C check address mode
158 IF (N.EQ.0) THEN
159C append mode
160 IHEP=NHEP+1
161 ELSE IF (N.GT.0) THEN
162C absolute position
163 IHEP=N
164 ELSE
165C relative position
166 IHEP=NHEP+N
167 END IF
168C check on IHEP
169 IF ((IHEP.LE.0).OR.(IHEP.GT.NMXHEP)) RETURN
170C 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)
182C 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
186C FLAG FOR PHOTOS...
187 QEDRAD(IHEP)=PHFLAG
188C update process:
189 DO IP=JMOHEP(1,IHEP),JMOHEP(2,IHEP)
190 IF(IP.GT.0)THEN
191C if there is a daughter at IHEP, mother entry at IP has decayed
192 IF(ISTHEP(IP).EQ.1)ISTHEP(IP)=2
193C 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
204CDECK ID>, FRAGMT.
205*CMZ :- -28/01/92 12.34.44 by Mike Seymour
206*-- Author : Luca Stanco
207C-----------------------------------------------------------------------
208 SUBROUTINE FRAGMT(I,J,K)
209C-----------------------------------------------------------------------
210C DUMMY SUBROUTINE: DELETE AND SET BDECAY='EURO'
211C IN MAIN PROGRAM IF YOU USE EURODEC DECAY PACKAGE
212C-----------------------------------------------------------------------
213 INTEGER I,J,K
214 WRITE (6,10)
215 10 FORMAT(/10X,'FRAGMT CALLED BUT NOT LINKED')
216 STOP
217 END
218CDECK ID>, HVCBVI.
219*CMZ :- -28/01/92 12.34.44 by Mike Seymour
220*-- Author : Mike Seymour
221C-----------------------------------------------------------------------
222 SUBROUTINE HVCBVI
223C-----------------------------------------------------------------------
224C DUMMY ROUTINE: DELETE IF YOU LINK TO BARYON NUMBER VIOLATN PACKAGE
225C-----------------------------------------------------------------------
226 WRITE (6,10)
227 10 FORMAT(/10X,'HVCBVI CALLED BUT NOT LINKED')
228 STOP
229 END
230CDECK ID>, HVHBVI.
231*CMZ :- -28/01/92 12.34.44 by Mike Seymour
232*-- Author : Mike Seymour
233C-----------------------------------------------------------------------
234 SUBROUTINE HVHBVI
235C-----------------------------------------------------------------------
236C DUMMY ROUTINE: DELETE IF YOU LINK TO BARYON NUMBER VIOLATN PACKAGE
237C-----------------------------------------------------------------------
238 WRITE (6,10)
239 10 FORMAT(/10X,'HERBVI CALLED BUT NOT LINKED')
240 STOP
241 END
242CDECK ID>, HWBAZF.
243*CMZ :- -26/04/91 11.11.54 by Bryan Webber
244*-- Author : Ian Knowles
245C-----------------------------------------------------------------------
246 SUBROUTINE HWBAZF(IPAR,JPAR,VEC1,VEC2,VEC3,VEC)
247C-----------------------------------------------------------------------
248C Azimuthal correlation functions for Collins' algorithm,
249C see I.G.Knowles, Comp. Phys. Comm. 58 (90) 271 for notation.
250C-----------------------------------------------------------------------
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
263C 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
272C 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
283C 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
292C 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
313CDECK ID>, HWBCON.
314*CMZ :- -11/10/01 12.01.52 by Peter Richardson
315*-- Author : Bryan Webber
316C-----------------------------------------------------------------------
317 SUBROUTINE HWBCON
318C-----------------------------------------------------------------------
319C MAKES COLOUR CONNECTIONS BETWEEN JETS
320C MODIFIED 12/10/97 BY BRW FOR SUSY PROCESSES
321C MODIFIED 11/01/01 BY PR FOR SPIN CORRELATIONS(PROBLEM WITH ORDER
322C OF DECAYS)
323C NEW VARAIBLE BACK TO ALLOW CODE TO SEARCH DOWN CHAIN
324C-----------------------------------------------------------------------
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)
336C---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
340C---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)
351C---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
354C---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)
358C--new bit to try and fix the problems for spin correlations
359C--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
372C---SPECIAL FOR GLUINO DECAYS
373 IF (ID.EQ.449) THEN
374 ID=IDHW(JC)
375C---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
380C---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
388C---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
396C---COULDNT FIND ONE
397 CALL HWWARN('HWBCON',101,*999)
398 5 JC=KC
399 ELSE
400C--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
408C--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
420C Wait for partner heavy quark to decay
421C RETURN
422C---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)
435C---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
442C---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
451C--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
459C---DIDN'T FIND PARTNER OF IHEP YET
460C CALL HWWARN('HWBCON',52,*20)
461 ENDIF
462 ENDIF
463 20 CONTINUE
464C---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
468C 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
473C 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
485CDECK ID>, HWBDED.
486*CMZ :- -22/04/96 13.54.08 by Mike Seymour
487*-- Author : Mike Seymour
488C-----------------------------------------------------------------------
489 SUBROUTINE HWBDED(IOPT)
490C FILL MISSING AREA OF DALITZ PLOT WITH 3-JET AND 2-JET+GAMMA EVENTS
491C IF (IOPT.EQ.1) SET UP EVENT RECORD
492C IF (IOPT.EQ.2) CLEAN UP EVENT RECORD AFTER SHOWERING
493C
494C********MODIFIED 13/11/00 BY BRW TO ALLOW MULTIPLE APPLICATION IN
495C*******SAME EVENT (FOR WW AND ZZ) N.B. NO CLEANUP CALLS FOR THESE!
496C-----------------------------------------------------------------------
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
509C---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
520C---ONLY APPLY THE CORRECTION TO HADRONIC DECAYS
521 IF (IDHW(JDAHEP(1,ICMF)).GT.12) GOTO 5
522C---GENERATE X1,X2 ACCORDING TO 1/((1-X1)*(1-X2))
523 100 CONTINUE
524C---CHOOSE X1
525 X(1)=1-(1-X1MAX)*((1-X1MIN)/(1-X1MAX))**HWRGEN(0)
526C---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)
533C---CALCULATE WEIGHT
534 W=2 * LOG((1-X1MIN)/(1-X1MAX))*LOG((1-X2MIN)/(1-X2MAX)) *
535 & (X(1)**2+X(2)**2)
536C---GENERATE UNWEIGHTED (X1,X2) PAIRS (EFFICIENCY IS ~50%)
537 IF (WMAX*HWRGEN(2).GT.W) GOTO 100
538C---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
544C---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)
550C---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)
557C---SWITCH OFF PHOTON EMISSION IN W DECAYS (THE M-E DOES NOT FACTORIZE)
558 IF (ICHRG(IDHW(ICMF)).NE.0) GAMFAC=0
559C---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
568C---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
582C---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
593C---STORE OLD MOMENTA
594 CALL HWVEQU(5,PHEP(1,JDAHEP(1,ICMF)),P1)
595 CALL HWVEQU(5,PHEP(1,JDAHEP(2,ICMF)),P2)
596C---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
606C---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
633C---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
642C---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)
663C---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))
670C---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))
675C---REORDER ENTRIES: IHEP=QUARK, JHEP=ANTI-QUARK, KHEP=EMITTED
676 IF (IHEP.EQ.LHEP) THEN
677 IHEP=JHEP
678 JHEP=LHEP
679 ENDIF
680C---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
689C---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
708C---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
757CDECK ID>, HWBDIS.
758*CMZ :- -17/05/94 09.33.08 by Mike Seymour
759*-- Author : Mike Seymour
760C-----------------------------------------------------------------------
761 SUBROUTINE HWBDIS(IOPT)
762C-----------------------------------------------------------------------
763C FILL MISSING AREA OF DIS PHASE-SPACE WITH 2+1-JET EVENTS
764C IF (IOPT.EQ.1) SET UP EVENT RECORD
765C IF (IOPT.EQ.2) CLEAN UP EVENT RECORD AFTER SHOWERING
766C-----------------------------------------------------------------------
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
782C---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)
798C---STORE OLD MOMENTA
799 CALL HWVEQU(5,P1,Q1)
800 CALL HWVEQU(5,P2,Q2)
801C---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)
816C---CHOOSE THE HADRONIC-PLANE CONFIGURATION, XP,ZP
817 IF (HWRGEN(0).LT.COMWGT) THEN
818C-----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
859C-----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
901C---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
924C---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
937C---DECIDE WHETHER TO MAKE AN EVENT HERE
938 IF (HWRGEN(4).GT.FAC+DIR) RETURN
939C---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)
961C---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
970C-----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
974C-----UNIFORMLY
975 PHI=(2*HWRGEN(6)-1)*PIFAC
976 ENDIF
977 ELSE
978C-----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
982C---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)
989C---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)
1005C---SPECIAL CASE FOR DIRECT PHOTON - COPY THE EXACT BEAM MOMENTUM
1006C---SHARE THE MISMATCH EQUALLY BETWEEN THE OUTGOING PARTONS
1007C---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))
1039C Decide which quark radiated and assign production vertices
1040 IF (BGF) THEN
1041C Boson-Gluon fusion case
1042 IF (1-ZP.LT.HWRGEN(0)) THEN
1043C 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
1049C 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
1056C QCD Compton case
1057 IF (1.LT.HWRGEN(0)*(1+(1-XP-ZP)**2+6*XP*(1-XP)*ZP*(1-ZP)))THEN
1058C 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
1064C 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
1071C---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
1101C---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
1115C---FACTORISATION SCALE
1116 EMSCA=SCALE
1117 EMIT=NEVHEP+NWGTS
1118 ELSEIF (IOPT.EQ.2) THEN
1119C---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
1228CDECK ID>, HWBDYP.
1229*CMZ :- -26/10/99 17.46.56 by Mike Seymour
1230*-- Author : Gennaro Corcella
1231C-----------------------------------------------------------------------
1232 SUBROUTINE HWBDYP(IOPT)
1233C MATRIX ELEMENT CORRECTIONS TO DRELL-YAN PROCESSES
1234C-----------------------------------------------------------------------
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
1254C-----CHOOSE WEIGHTS
1255 COMWGT1=0.1
1256 COMWGT2=0.55
1257C---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)
1264C-----SET THE VECTOR BOSON RAPIDITY
1265 Y=HALF*LOG((PHEP(4,ICMF)+PHEP(3,ICMF))/
1266 & (PHEP(4,ICMF)-PHEP(3,ICMF)))
1267C------SET PARTICLE IDENTIES
1268c------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)
1277C---STORE OLD MOMENTA
1278C------VECTOR BOSON MOMENTUM
1279 CALL HWVEQU(5,PHEP(1,ICMF),PBOS)
1280C----QUARK MOMENTUM
1281 CALL HWVEQU(5,PHEP(1,JMOHEP(1,ICMF)),P1)
1282C------ANTIQUARK MOMENTUM
1283 CALL HWVEQU(5,PHEP(1,JMOHEP(2,ICMF)),P2)
1284C-------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)
1287C------LEPTON MOMENTA IN THE BOSON REST FRAME
1288 CALL HWULOF(PHEP(1,ICMF),P2,P2N)
1289 CALL HWULOF(PHEP(1,ICMF),P3,P3N)
1290C------AZ=AZIMUTHAL ANGLE OF P3N
1291 AZ=ATAN2(P3N(2),P3N(1))
1292 CZ=COS(AZ)
1293 SZ=SIN(AZ)
1294C------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)
1300C------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)
1309C------ Q - QBAR ENERGY FRACTIONS (BORN PROCESS)
1310c---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))
1313C------ 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)
1316C-------CONSIDER Q(QBAR) IN THE INITIAL STATE
1317 RN=HWRGEN(9)
1318 IF (RN.LT.COMWGT1) THEN
1319C-------NO GLUON IN THE INITIAL STATE
1320 GLUIN=.FALSE.
1321C---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)
1328C---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)
1340C----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))
1343c---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
1346C-----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)
1349C------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))
1353C-------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
1359C--------GLUON IN THE INITIAL STATE
1360 GLUIN=.TRUE.
1361C---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)
1367C---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)
1378C--------INITIAL STATE GLUON COMING FROM HADRON 1
1379 IF (RN.LE.COMWGT2) THEN
1380 GP=.TRUE.
1381C--------ENERGY FRACTIONS and PDFs
1382c---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))
1385c---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
1393C-------INITIAL STATE GLUON COMING FROM HADRON 2
1394 GP=.FALSE.
1395C-------ENERGY FRACTIONS AND PDFs
1396c---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))
1399c---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)
1408C-------CHOOSE WHICH PARTON WILL EMIT
1409c---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
1414C-------FINAL WEIGHT FOR ALL THE CONSIDERED OPTIONS
1415 W1=GLUFAC*W*WW*ETA1*ETA2/(XI1*XI2)
1416 ENDIF
1417C--------ADD ONE MORE GLUON
1418 IF (W1.GT.HWRGEN(4)) THEN
1419 NTMP=NEVHEP+NWGTS
1420 ELSE
1421 RETURN
1422 ENDIF
1423C---------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
1445C----MOMENTA IN THE V-REST FRAME WITH NON EMITTER ALONG THE Z AXIS
1446C----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)
1500C---------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)
1516c------LEPTON MOMENTA IN THE BOSON REST FRAME, WITH THE DIRECTION
1517C------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)
1528C----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)
1558C------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)
1565C----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)
1587C-----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)
1594C--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)
1605C-----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
1609C---STATUS OF EMITTER/NON EMITTER
1610 ISTHEP(IHEP)=110+EMIT
1611 ISTHEP(JHEP)=110+NOEMIT
1612 ELSE
1613C-----GLUON COMING FROM THE 1ST HADRON
1614 IF (GP) THEN
1615 KHEP=CHEP-2
1616 ISTHEP(KHEP)=111
1617C----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
1630C-------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
1643C------GLUON COMING FROM THE HADRON 2
1644 IF (.NOT.GP) THEN
1645 KHEP=CHEP-1
1646 ISTHEP(KHEP)=112
1647C-------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
1660C-------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)
1680C---------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
1722C---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)
1776C--------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
1824C--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
1835CDECK ID>, HWBFIN.
1836*CMZ :- -26/04/91 10.18.56 by Bryan Webber
1837*-- Author : Bryan Webber
1838C-----------------------------------------------------------------------
1839 SUBROUTINE HWBFIN(IHEP)
1840C-----------------------------------------------------------------------
1841C DELETES INTERNAL LINES FROM SHOWER, MAKES COLOUR CONNECTION INDEX
1842C AND COPIES INTO /HEPEVT/ IN COLOUR ORDER.
1843C-----------------------------------------------------------------------
1844 INCLUDE 'HERWIG65.INC'
1845 INTEGER IHEP,ID,IJET,KHEP,IPAR,JPAR,NXPAR,IP,JP
1846 IF (IERROR.NE.0) RETURN
1847C---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))
1861C---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
1870C---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
1884C---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
1899C---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
1941CDECK ID>, HWBGEN.
1942*CMZ :- -14/10/99 18.04.56 by Mike Seymour
1943*-- Author : Bryan Webber
1944C-----------------------------------------------------------------------
1945 SUBROUTINE HWBGEN
1946C-----------------------------------------------------------------------
1947C BRANCHING GENERATOR WITH INTERFERING GLUONS
1948C HWBGEN EVOLVES QCD JETS ACCORDING TO THE METHOD OF
1949C G.MARCHESINI & B.R.WEBBER, NUCL. PHYS. B238(1984)1
1950C-----------------------------------------------------------------------
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
1959C---CHECK THAT EMSCA IS SET
1960 IF (EMSCA.LE.ZERO) CALL HWWARN('HWBGEN',200,*999)
1961 IF (HARDME) THEN
1962C---FORCE A BRANCH INTO THE `DEAD ZONE' IN E+E-
1963 JPR=IPROC/10
1964C**********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)
1966C**********END FIX
1967C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DIS
1968 IF (IPRO.EQ.90) CALL HWBDIS(1)
1969C---FORCE A BRANCH INTO THE `DEAD ZONE' IN DRELL-YAN PROCESSES
1970 IF (IPRO.EQ.13.OR.IPRO.EQ.14) CALL HWBDYP(1)
1971C---FORCE A BRANCH INTO THE `DEAD ZONE' IN TOP DECAYS
1972 CALL HWBTOP
1973 ENDIF
1974C---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
1999C---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
2011C---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)
2053C---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
2071C---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
2076C---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
2087C---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
2096C---ENTER PARTON JET IN /HEPEVT/
2097 CALL HWBFIN(IHEP)
2098 ELSE
2099C---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
2118C---COMBINE JETS
2119 ISTAT=20
2120 CALL HWBJCO
2121 ENDIF
2122 IF (.NOT.FROST) THEN
2123C---ATTACH SPECTATORS
2124 ISTAT=30
2125 CALL HWSSPC
2126 ENDIF
2127 IF (FROST) THEN
2128C---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
2133C---CONNECT COLOURS
2134 CALL HWBCON
2135 ISTAT=40
2136 LASHEP=NHEP
2137 IF (HARDME) THEN
2138C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN E+E-
2139 IF (IPROC/10.EQ.10) CALL HWBDED(2)
2140C---CLEAN UP IF THERE WAS A BRANCH IN THE `DEAD ZONE' IN DIS
2141 IF (IPRO.EQ.90) CALL HWBDIS(2)
2142C---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
2145C---IF THE CLEAN-UP OPERATION ADDED ANY PARTONS TO THE EVENT RECORD
2146C IT MIGHT NEED RESHOWERING
2147 IF (NHEP.GT.LASHEP) THEN
2148 LASHEP=NHEP
2149 GOTO 10
2150 ENDIF
2151 999 END
2152CDECK ID>, HWBGUP.
2153*CMZ :- -16/07/02 09.40.25 by Peter Richardson
2154*-- Author : Peter Richardson
2155C----------------------------------------------------------------------
2156 SUBROUTINE HWBGUP(ISTART,ICMF)
2157C----------------------------------------------------------------------
2158C Makes the colour connections and performs the parton shower
2159C for events read in from the GUPI (Generic User Process Interface)
2160C event common block
2161C----------------------------------------------------------------------
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)
2171C--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
2176C--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
2195C--now search for the partner
2196C--first search for the flavour partner if not looking for colour partner
2197C--search for the flavour partner of the particle
2198C--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.
2202C--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
2208C--antiflavour partner
2209 IF(JDAHEP(2,JLOC(I)).EQ.0) THEN
2210C--pair incoming particle with outgoing particle
2211C-- 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
2217C--pair incoming particle with incoming antiparticle
2218C-- 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
2225C--make the connection
2226 IF(FOUND) THEN
2227 JMOHEP(2,K) = JLOC(I)
2228 JDAHEP(2,JLOC(I)) = K
2229 ENDIF
2230 ENDIF
2231C--flavour partner
2232 IF(JMOHEP(2,JLOC(I)).EQ.0.AND.(.NOT.FOUND)) THEN
2233C--pair incoming antiparticle with outgoing antiparticle
2234C-- 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
2240C--pair incoming antiparticle with incoming particle
2241C-- 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
2248C--make the connection
2249 IF(FOUND) THEN
2250 JDAHEP(2,K) = JLOC(I)
2251 JMOHEP(2,JLOC(I)) = K
2252 ENDIF
2253 ENDIF
2254C--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
2263C--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
2275C--now the bit to find colour partners
2276 FOUND = .FALSE.
2277C--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
2291C--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
2327C--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
2332C--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
2340C--perform the shower
2341 30 CALL HWBGEN
2342 999 END
2343CDECK ID>, HWBJCO.
2344*CMZ :- -30/09/02 09.19.58 by Peter Richardson
2345*-- Author : Bryan Webber
2346C-----------------------------------------------------------------------
2347 SUBROUTINE HWBJCO
2348C-----------------------------------------------------------------------
2349C COMBINES JETS WITH REQUIRED KINEMATICS
2350C-----------------------------------------------------------------------
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
2370C---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
2391C---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
2398C---SPACELIKE JETS: FIND SPACELIKE PARTONS
2399 IF (NP.NE.2) CALL HWWARN('HWBJCO',103,*999)
2400C---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)
2406C---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
2420C---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
2427C---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
2432C---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
2442C---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
2458C---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
2482C---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
2493C---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
2499C---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)
2504C---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
2516C---calculate boost for struck parton
2517C PC is momentum of outgoing parton(s)
2518 IP2=JDAHEP(2,ICM)
2519 IF (.NOT.DISLOW) THEN
2520C---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)
2529C---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))
2554C---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
2570C 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
2596C--change to preserve either long mom or rapidity rather than long mom
2597C--by PR and BRW 30/9/02
2598 IF (PRESPL) THEN
2599C--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
2603C--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
2608C---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
2637C---NON-PARTON JET
2638 ISTHEP(IJT)=3
2639 ENDIF
2640 160 CONTINUE
2641 ENDIF
2642 ISTHEP(ICM)=120
2643 ELSE
2644C---TIMELIKE JETS
2645C---SPECIAL CASE: IF HARD PROCESS IS W/Z DECAY, PERFORM KINEMATIC
2646C 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
2656C 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)
2663C--- 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
2673C---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
2683C---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
2689C---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
2709C---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
2719C---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)
2734C--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)
2747C--End mod
2748 ELSE
2749 CALL HWVEQU(5,PC,PR)
2750 ENDIF
2751C---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
2754C---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
2769C---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)
2801C---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))
2817C---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))
2825C--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))
2829C--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
2834C---NON-PARTON JET
2835 ISTHEP(IHEP)=190
2836 ENDIF
2837 230 CONTINUE
2838 IF (ISTHEP(ICM).EQ.110) ISTHEP(ICM)=120
2839C---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)))
2847C--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)))
2853C---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
2870CDECK ID>, HWBMAS.
2871*CMZ :- -26/04/91 11.11.54 by Bryan Webber
2872*-- Author : Bryan Webber
2873C-----------------------------------------------------------------------
2874 SUBROUTINE HWBMAS
2875C-----------------------------------------------------------------------
2876C Passes backwards through a jet cascade calculating the masses
2877C and magnitudes of the longitudinal and transverse three momenta.
2878C Components given relative to direction of parent for a time-like
2879C vertex and with respect to z-axis for space-like vertices.
2880C
2881C On input PPAR(1-5,*) contains:
2882C (E*sqrt(Xi),Xi,3-mom (if external),E,M-sq (if external))
2883C
2884C On output PPAR(1-5,*) (if TMPAR(*)), containts:
2885C (P-trans,Xi or Xilast,P-long,E,M)
2886C-----------------------------------------------------------------------
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
2896C Find parent and partner of this branch
2897 IPAR=JMOPAR(1,JPAR)
2898 KPAR=JPAR+1
2899C Determine type of branching
2900 IF (TMPAR(IPAR)) THEN
2901C Time-like branching
2902C Compute mass of parent
2903 EXI=PPAR(1,JPAR)*PPAR(1,KPAR)
2904 PPAR(5,IPAR)=PPAR(5,JPAR)+PPAR(5,KPAR)+2.*EXI
2905C Compute three momentum of parent
2906 PISQ=PPAR(4,IPAR)*PPAR(4,IPAR)-PPAR(5,IPAR)
2907 PPAR(3,IPAR)=HWUSQR(PISQ)
2908C---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)
2934C---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
2959C 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
2968C Space-like branching
2969C Re-arrange such that JPAR is time-like
2970 IF (TMPAR(KPAR)) THEN
2971 KPAR=JPAR
2972 JPAR=JPAR+1
2973 ENDIF
2974C 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
2983C 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
2992CDECK ID>, HWBRAN.
2993*CMZ :- -14/10/99 18.04.56 by Mike Seymour
2994*-- Author : Bryan Webber & Mike Seymour
2995C-----------------------------------------------------------------------
2996 SUBROUTINE HWBRAN(KPAR)
2997C-----------------------------------------------------------------------
2998C BRANCHES TIMELIKE PARTON KPAR INTO TWO, PUTS PRODUCTS
2999C INTO NPAR+1 AND NPAR+2, AND INCREASES NPAR BY TWO
3000C-----------------------------------------------------------------------
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
3013C---SET SQRK(M,N) TO THE PROBABILITY THAT A GLUON WILL NOT PRODUCE A
3014C 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)
3038C--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
3050C--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
3058C--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
3062C--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
3070C--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
3084C--FIX AT HARD PROCESS SCALE IF POINTER NOT YET SET
3085C (CAN HAPPEN IN SUSY EVENTS)
3086 QMAX=EMSCA**2
3087 ENDIF
3088 ELSE
3089 QMAX=ENOW**2*PPAR(2,MPAR)
3090 ENDIF
3091C--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
3110C--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
3117C---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
3121C---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
3138C---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))
3157C---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
3178C--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
3192C--ENHANCE GLUON AND PHOTON EMISSION BY A FACTOR OF TWO IF THIS BRANCH
3193C 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
3196C--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)
3209C---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
3225C--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
3250C--BRANCHING HAS OCCURRED
3251 ZMIN=HWBVMC(ID2)/QNOW
3252 ZMAX=1.-ZMIN
3253 IF (ID.EQ.13) THEN
3254 IF (ID2.EQ.13) THEN
3255C--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))
3260C--CHOOSE Z1 DISTRIBUTED ON (ZMIN,ZMAX)
3261C 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
3268C--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
3276C--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
3283C--QUARK OR ANTIQUARK BRANCHING
3284 IF (ID2.EQ.13) THEN
3285C--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
3295C--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
3305C--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
3316C--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
3326C--BRANCHING REJECTED: REDUCE Q AND REPEAT
3327 QMAX=QNOW
3328 QLST=QNOW
3329 QNOW=-1.
3330 GOTO 5
3331 ENDIF
3332 ENDIF
3333C--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
3338C----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)))
3349C-----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)
3352C-----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
3360C---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)
3367C---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
3376C---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)
3383C---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
3411 IDPAR(NPAR)=ID2
3412 TMPAR(NPAR)=.TRUE.
3413 PPAR(1,NPAR)=QNOW*Z2
3414 PPAR(2,NPAR)=XI
3415 PPAR(4,NPAR)=ENOW*Z2
3416C---NEW MOTHER-DAUGHTER RELATIONS
3417 JDAPAR(1,KPAR)=MPAR
3418 JDAPAR(2,KPAR)=NPAR
3419 JMOPAR(1,MPAR)=KPAR
3420 JMOPAR(1,NPAR)=KPAR
3421C---NEW COLOUR CONNECTIONS
3422 JCOPAR(3,KPAR)=NPAR
3423 JCOPAR(4,KPAR)=MPAR
3424 JCOPAR(1,MPAR)=NPAR
3425 JCOPAR(2,MPAR)=KPAR
3426 JCOPAR(1,NPAR)=KPAR
3427 JCOPAR(2,NPAR)=MPAR
3428C
3429 ENDIF
3430 ENDIF
3431 IF (QNOW.LT.ZERO) THEN
3432C--BRANCHING STOPS
3433 IF (ID.EQ.IDPAR(2).AND.PPAR(5,2).GT.1D-6) THEN
3434 PPAR(5,KPAR)=PPAR(5,2)**2
3435 ELSE
3436 PPAR(5,KPAR)=RMASS(ID)**2
3437 ENDIF
3438 PMOM=PPAR(4,KPAR)**2-PPAR(5,KPAR)
3439 IF (PMOM.LT.-1E-6) CALL HWWARN('HWBRAN',104,*999)
3440 IF (PMOM.LT.ZERO) PMOM=ZERO
3441 PPAR(3,KPAR)=SQRT(PMOM)
3442 JDAPAR(1,KPAR)=0
3443 JDAPAR(2,KPAR)=0
3444 JCOPAR(3,KPAR)=0
3445 JCOPAR(4,KPAR)=0
3446 ENDIF
3447 999 END
3448CDECK ID>, HWBRCN.
3449*CMZ :- -31/03/00 17:54:05 by Peter Richardson
3450*-- Author : Peter Richardson
3451C-----------------------------------------------------------------------
3452 SUBROUTINE HWBRCN
3453C-----------------------------------------------------------------------
3454C SUBROUTINE TO REPLACE HWBCON IN RPARITY VIOLATING SUSY
3455C BASED ON HWBCON BY BRW
3456C-----------------------------------------------------------------------
3457 INCLUDE 'HERWIG65.INC'
3458 INTEGER IHEP,IST,ID,JC,JD,JHEP,IDP,IDM,IDP2,IDM2,
3459 & RHEP,IST2,ORG,ANTC,XHEP,IP,COLP
3460 LOGICAL BVVUSE,BVVHRD,BVDEC1,BVDEC2,COLRD,ACOLRD,BVHRD,BVHRD2,
3461 & BVDEC3
3462C--logical functions to decide if baryon number violating
3463C--BVDEC1 DELTAB=+1
3464 BVDEC1(IP) = ((IDHW(IP).GE.419.AND.IDHW(IP).LE.424).OR.
3465 & IDHW(IP).EQ.411.OR.IDHW(IP).EQ.412.OR.
3466 & IDHW(IP).EQ.449).AND.IDHW(JDAHEP(1,IP)).LE.6.
3467 & AND.IDHW(JDAHEP(1,IP)+1).LE.6.AND.
3468 & IDHW(JDAHEP(2,IP)).LE.6
3469C--BVDEC2 DELTAB=-1
3470 BVDEC2(IP) = ((IDHW(IP).GE.413.AND.IDHW(IP).LE.418).OR.
3471 & IDHW(IP).EQ.405.OR.IDHW(IP).EQ.406.OR.
3472 & IDHW(IP).EQ.449).AND.
3473 & IDHW(JDAHEP(1,IP)).GE.7.AND.IDHW(JDAHEP(1,IP)).LE.12.AND.
3474 & IDHW(JDAHEP(1,IP)+1).GE.7.AND.IDHW(JDAHEP(1,IP)+1).LE.12.AND.
3475 & IDHW(JDAHEP(2,IP)).GE.7.AND.IDHW(JDAHEP(2,IP)).LE.12
3476C--Neutralino and Chargino Decays
3477 BVDEC3(IP) = ((IDHW(IP).GE.450.AND.IDHW(IP).LE.457).AND.
3478 & (IDHW(JDAHEP(1,IP)).LE.12.AND.IDHW(JDAHEP(1,IP)+1).LE.12.
3479 & .AND.IDHW(JDAHEP(2,IP)).LE.12))
3480C--Now the hard vertices
3481 BVHRD(IP) = IDHW(IP).EQ.15.AND.IDHW(JMOHEP(1,IP)).LE.12.
3482 & AND.IDHW(JMOHEP(2,IP)).LE.12.AND.IDHW(JDAHEP(1,IP)).LE.12.
3483 & AND.IDHW(JDAHEP(2,IP)).GE.449.AND.IDHW(JDAHEP(2,IP)).LE.457
3484 BVHRD2(IP) = IDHW(IP).EQ.15.AND.IDHW(JMOHEP(1,IP)).LE.12.
3485 & AND.IDHW(JMOHEP(2,IP)).LE.12.AND.IDHW(JDAHEP(1,IP)).GE.198.
3486 & AND.IDHW(JDAHEP(1,IP)).LE.207.
3487 & AND.ABS(IDHEP(JDAHEP(2,IP))).GT.1000000
3488C--Those particles which are coloured
3489 COLRD(IP) = IP.LE.6.OR.IP.EQ.13.OR.IP.EQ.449.OR.
3490 & (IP.GE.401.AND.IP.LE.406).OR.(IP.GE.413.AND.IP.LE.418).OR.
3491 & (IP.GE.115.AND.IP.LE.120).OR.IP.EQ.59
3492C--Those particles which are anticoloured
3493 ACOLRD(IP) = (IP.GE.7.AND.IP.LE.12).OR.IP.EQ.13.OR.IP.EQ.449.OR.
3494 & (IP.GE.407.AND.IP.LE.412).OR.(IP.GE.419.AND.IP.LE.424).OR.
3495 & (IP.GE.109.AND.IP.LE.114).OR.IP.EQ.59
3496 IF (IERROR.NE.0) RETURN
3497C--Added 31/03/00 PR
3498 IF(NHEP.GT.NMXHEP) CALL HWWARN('HWBRCN',101,*999)
3499 COLP = 0
3500 IF(COLUPD.AND.HRDCOL(1,3).NE.0) THEN
3501 JD = 0
3502 DO IHEP = HRDCOL(1,3),HRDCOL(1,3)+4
3503 JD = JD+1
3504 IF(JD.NE.3) THEN
3505 JMOHEP(2,IHEP) = HRDCOL(1,JD)
3506 JDAHEP(2,IHEP) = HRDCOL(2,JD)
3507 ENDIF
3508 ENDDO
3509 COLUPD=.FALSE.
3510 DO IHEP=1,5
3511 DO JHEP=1,2
3512 HRDCOL(JHEP,IHEP)=0
3513 ENDDO
3514 ENDDO
3515 ELSEIF(COLUPD) THEN
3516 RETURN
3517 ENDIF
3518 DO 110 IHEP=1,NHEP
3519 IST=ISTHEP(IHEP)
3520 JD =0
3521 BVVUSE = .FALSE.
3522 BVVHRD = .FALSE.
3523C---LOOK FOR PARTONS WITHOUT COLOUR MOTHERS
3524 IF ((IST.LT.145.OR.IST.GT.152).AND.IST.NE.155) GOTO 110
3525 IF (JMOHEP(2,IHEP).EQ.0) THEN
3526C---FIND COLOUR-CONNECTED PARTON
3527 IF(IST.EQ.155.AND.ABS(IDHEP(IHEP)).EQ.6) THEN
3528 JC = JMOHEP(1,IHEP)
3529 ELSEIF(IST.EQ.155) THEN
3530 GOTO 110
3531 ELSE
3532 JC=JMOHEP(1,IHEP)
3533 ENDIF
3534 IF (IST.NE.152) JC=JMOHEP(1,JC)
3535C--Correction for BV
3536 IF(HRDCOL(1,1).NE.0) THEN
3537 IDP = IDHW(HRDCOL(1,1))
3538 IDP2 = 0
3539 ELSE
3540 IDP = 0
3541 IDP2 = 0
3542 ENDIF
3543 IDM = JMOHEP(1,JC)
3544 IF(BVDEC1(IDM).OR.BVDEC2(IDM)) THEN
3545 IF(IDHW(IDM).EQ.449.AND.JDAHEP(1,IDM).EQ.JC) THEN
3546 JC=JMOHEP(2,JC)
3547 ELSE
3548 JD = JMOHEP(2,JC)
3549 JC = IDM
3550 IF(JC.EQ.JD) JD= JDAHEP(2,JC-1)
3551 BVVUSE = .TRUE.
3552 ENDIF
3553C--NEW FOR BV HARD PROCESS
3554 ELSEIF(BVHRD(IDM)) THEN
3555 IF(IDHW(JDAHEP(2,JMOHEP(1,JC))).EQ.449) THEN
3556 JD = JMOHEP(2,JC)
3557 IDM2 = JDAHEP(2,HRDCOL(1,2))
3558 IF(JD.EQ.IDM2) JD = HRDCOL(1,1)
3559 IF(JC.EQ.JDAHEP(2,IDM2).AND.COLRD(IDHW(IHEP))) THEN
3560 JC = JMOHEP(2,JC)
3561 ELSEIF(JC.EQ.IDM2) THEN
3562 IF(JDAHEP(2,JMOHEP(2,JC)).EQ.JC) THEN
3563 JC = JMOHEP(2,JC)
3564 ELSE
3565 JMOHEP(2,IHEP)=JMOHEP(2,JC)
3566 GOTO 110
3567 ENDIF
3568 ELSE
3569 JC = HRDCOL(1,1)
3570 BVVUSE = .TRUE.
3571 BVVHRD = .TRUE.
3572 IF(ACOLRD(IDHW(IHEP))) JC = JD
3573 IF(JC.EQ.IDM2) GOTO 110
3574 ENDIF
3575 ELSE
3576 JC =JMOHEP(2,JC)
3577 BVVUSE = .TRUE.
3578 BVVHRD = .TRUE.
3579 ENDIF
3580 ELSEIF(BVHRD2(IDM)) THEN
3581 JD = JMOHEP(2,JC)
3582 IF(JC.EQ.JDAHEP(2,HRDCOL(1,2))) THEN
3583 JMOHEP(2,IHEP)=JMOHEP(2,JC)
3584 GOTO 110
3585 ENDIF
3586 IF(JD.EQ.JDAHEP(2,HRDCOL(1,2))) JD = HRDCOL(1,1)
3587 BVVUSE=.TRUE.
3588 BVVHRD = .TRUE.
3589 IF(JC.EQ.JDAHEP(2,HRDCOL(1,2))) THEN
3590 JC = JMOHEP(2,JC)
3591 ELSE
3592 JC = HRDCOL(1,1)
3593 ENDIF
3594 ELSE
3595 JC =JMOHEP(2,JC)
3596 ENDIF
3597 IF (JC.EQ.0) CALL HWWARN('HWBCON',51,*110)
3598C---FIND SPECTATOR WHEN JC IS DECAYED HEAVY QUARK OR SUSY PARTICLE
3599 IF (ISTHEP(JC).EQ.155) THEN
3600 IF (IDHEP(JMOHEP(1,JC)).EQ.94) THEN
3601C---DECAYED BEFORE HADRONIZING
3602 IF(BVVHRD) THEN
3603 JHEP = JC
3604 ELSEIF(BVVUSE) THEN
3605 JHEP=JDAHEP(2,JC-1)
3606 ELSE
3607 JHEP=JMOHEP(2,JC)
3608 ENDIF
3609 IF(JHEP.EQ.0.AND.ABS(IDHEP(JC)).EQ.6) THEN
3610 JHEP = JMOHEP(1,JMOHEP(1,JC))
3611 IF(BVDEC1(JMOHEP(1,JHEP)).OR.BVDEC2(JMOHEP(1,JHEP))) THEN
3612 JC = JHEP
3613 JHEP = JDAHEP(2,JC-1)
3614 ELSE
3615 JHEP = 0
3616 ENDIF
3617 ENDIF
3618 IF(BVVUSE.AND.ABS(IDHEP(JHEP)).GT.1000000.AND.
3619 & ISTHEP(JHEP).NE.155.OR.JHEP.EQ.0) GOTO 110
3620 ID=IDHW(JHEP)
3621 IF (ISTHEP(JHEP).EQ.155) THEN
3622C---SPECIAL FOR GLUINO DECAYS
3623 IF (ID.EQ.449) THEN
3624 ID=IDHW(JC)
3625 IF(BVVUSE) THEN
3626 ID=IDHW(IHEP)
3627 IF(ID.LE.6.OR.ID.EQ.13.OR.
3628 & (ID.GE.115.AND.ID.LE.120)) THEN
3629 ID = 7
3630 ELSE
3631 ID = 1
3632 ENDIF
3633 ENDIF
3634 CALL HWBRC1(JC,ID,JHEP,.TRUE.,*999)
3635 IF(BVVUSE.AND.JMOHEP(1,JC).EQ.JMOHEP(1,JD)) JC =JD
3636 ELSE
3637 JC=JDAHEP(2,JHEP)
3638 IF(COLRD(IDHW(IHEP)).AND.IDHW(JDAHEP(1,JHEP)).EQ.449)
3639 & JC=JDAHEP(1,JHEP)
3640 IF(BVVUSE.AND.JMOHEP(1,JC).EQ.JMOHEP(1,JD)) JC =JD
3641 ENDIF
3642 ELSE
3643 IF(BVVUSE) THEN
3644 IF(BVDEC2(JMOHEP(1,JHEP)).OR.JD.NE.JHEP.OR.
3645 & BVHRD(JMOHEP(1,JHEP)).OR.BVHRD2(JMOHEP(1,JHEP))) THEN
3646 JC = JD
3647 GOTO 100
3648 ELSE
3649 JMOHEP(2,IHEP)=JHEP
3650 ID = IDHW(JHEP)
3651 IF((ID.GE.7.AND.ID.LE.12).OR.
3652 & (ID.GE.109.AND.ID.LE.114)) JMOHEP(2,JHEP)=IHEP
3653 ENDIF
3654 ELSE
3655C--new for particles connected to BV
3656 IDM = JMOHEP(1,JHEP)
3657 IF(BVDEC1(IDM).OR.BVHRD(IDM).OR.BVHRD2(IDM)) THEN
3658 JC = JHEP
3659 IF(ABS(IDHEP(IHEP)).LT.1000000) GOTO 100
3660 JMOHEP(2,IHEP)=JHEP
3661 GOTO 110
3662 ENDIF
3663C--new for top's from BV
3664 ID = IDHW(JC)
3665 IDP = JMOHEP(1,JMOHEP(1,JMOHEP(1,JC)))
3666 IF((ID.EQ.6.AND.(BVDEC1(IDP))).
3667 & OR.(ID.EQ.12.AND.BVDEC2(IDP)).
3668 & OR.((ID.EQ.12.OR.ID.EQ.449).AND.BVHRD(IDP))) THEN
3669 JMOHEP(2,IHEP)=JHEP
3670 IF(JDAHEP(2,JHEP).EQ.JC) JDAHEP(2,JHEP)=IHEP
3671 ELSE
3672 IF((IDHW(IHEP).GE.7.AND.IDHW(IHEP).LE.12.
3673 & AND.IDHW(JHEP).GE.7.AND.IDHW(JHEP).LE.12).OR.
3674 & (IDHW(IHEP).LE.6.AND.IDHW(JHEP).LE.6)) THEN
3675 JMOHEP(2,IHEP)=JHEP
3676 ELSE
3677 JMOHEP(2,IHEP)=JHEP
3678 IF((COLRD(IDHW(IHEP)).AND.ACOLRD(IDHW(JHEP))).OR.
3679 & (.NOT.COLRD(IDHW(IHEP)).AND.
3680 & .NOT.ACOLRD(IDHW(JHEP)))) THEN
3681 IF(JDAHEP(2,JHEP).EQ.0) THEN
3682 JDAHEP(2,JHEP)=IHEP
3683 ELSEIF(JMOHEP(2,JDAHEP(2,JHEP)).NE.JHEP) THEN
3684 JDAHEP(2,JHEP)=IHEP
3685 ENDIF
3686 ELSE
3687 IF(JMOHEP(2,JHEP).EQ.JC) JMOHEP(2,JHEP)=IHEP
3688 ENDIF
3689 ENDIF
3690 ENDIF
3691 ENDIF
3692 GOTO 110
3693 ENDIF
3694 ELSE
3695 JC=JMOHEP(2,JC)
3696 ENDIF
3697 ENDIF
3698 100 CONTINUE
3699 IF(BVVUSE.AND.ABS(IDHEP(JC)).LT.1000000.AND.JC.NE.JD
3700 & .AND.JD.NE.0.AND.JD.NE.JMOHEP(1,JC)) JC = JD
3701 IF(BVVUSE.AND.ABS(IDHEP(JC)).GT.1000000) THEN
3702 IF(COLRD(IDHW(IHEP)).AND..NOT.BVVHRD) GOTO 110
3703 ENDIF
3704 IF(BVVUSE.AND.ISTHEP(JC).EQ.149) JC=JMOHEP(1,JMOHEP(1,JC))
3705C--SEARCH IN THE JET
3706 IF((ISTHEP(JC).GT.145.AND.ISTHEP(JC).LT.152).AND.
3707 & ISTHEP(IHEP).EQ.155) THEN
3708 JMOHEP(2,IHEP) = JC
3709 GOTO 110
3710 ENDIF
3711 CALL HWBRC2(COLP,IHEP,JC,.TRUE.,BVVUSE,BVVHRD)
3712 IF(COLP.NE.0) THEN
3713 JMOHEP(2,IHEP) = COLP
3714 IF(COLRD(IDHW(IHEP)).AND.ACOLRD(IDHW(COLP)).
3715 & AND.JDAHEP(2,COLP).EQ.0)
3716 & JDAHEP(2,COLP) = IHEP
3717 IF((IDHW(IHEP).GE.7.AND.IDHW(IHEP).LE.12).AND.
3718 & (IDHW(COLP).GE.7.AND.IDHW(COLP).LE.12)) THEN
3719 IF(JMOHEP(2,COLP).EQ.0) JMOHEP(2,COLP) = IHEP
3720 ENDIF
3721 ENDIF
3722 ENDIF
3723 110 CONTINUE
3724C---BREAK COLOUR CONNECTIONS WITH PHOTONS modified for Rslash
3725 IHEP=1
3726 130 IF (IHEP.LE.NHEP) THEN
3727 IF (IDHW(IHEP).EQ.59 .AND. ISTHEP(IHEP).EQ.149.AND.
3728 & (JMOHEP(2,IHEP).NE.IHEP.OR.JDAHEP(2,IHEP).NE.IHEP)) THEN
3729 IF(JMOHEP(2,IHEP).NE.0) THEN
3730 IF (JDAHEP(2,JMOHEP(2,IHEP)).EQ.IHEP)
3731 & JDAHEP(2,JMOHEP(2,IHEP))=JDAHEP(2,IHEP)
3732 ENDIF
3733 IF (JDAHEP(2,IHEP).NE.0) THEN
3734 IF (JMOHEP(2,JDAHEP(2,IHEP)).EQ.IHEP)
3735 & JMOHEP(2,JDAHEP(2,IHEP))=JMOHEP(2,IHEP)
3736 ENDIF
3737 DO RHEP=1,NHEP
3738 IST=ISTHEP(RHEP)
3739 IF((IST.GE.147.AND.IST.LE.149).AND.JDAHEP(2,RHEP).EQ.IHEP)
3740 & JDAHEP(2,RHEP)=JMOHEP(2,IHEP)
3741 ENDDO
3742 DO RHEP=1,NHEP
3743 IST=ISTHEP(RHEP)
3744 IF((IST.GE.147.AND.IST.LE.149).AND.JMOHEP(2,RHEP).EQ.IHEP)
3745 & JMOHEP(2,RHEP) = JDAHEP(2,IHEP)
3746 ENDDO
3747 JMOHEP(2,IHEP)=IHEP
3748 JDAHEP(2,IHEP)=IHEP
3749 ENDIF
3750 IHEP=IHEP+1
3751 GOTO 130
3752 ENDIF
3753C--Update the BV anticolour corrections
3754 DO 210 IHEP=1,NHEP+1
3755 IF(IHEP.EQ.1) GOTO 210
3756 IST2 = 0
3757 IF(IHEP.EQ.NHEP+1) THEN
3758 ANTC = HRDCOL(1,1)
3759 IF(ANTC.EQ.0.OR.(IDHW(JMOHEP(1,HRDCOL(1,2))).LE.6)) GOTO 210
3760 IST=155
3761 XHEP=HRDCOL(1,2)
3762 IF(ANTC.EQ.JDAHEP(2,XHEP)) ANTC=JDAHEP(1,JDAHEP(1,ANTC))
3763 IF(ANTC.NE.0.AND.JDAHEP(1,ANTC).NE.0) IST2=ISTHEP(ANTC)
3764 ELSE
3765 ANTC = JDAHEP(2,IHEP-1)
3766 IF(ANTC.NE.0) IST2=ISTHEP(ANTC)
3767 IST=ISTHEP(IHEP)
3768 IDM = IDHW(IHEP)
3769 XHEP=IHEP
3770 ENDIF
3771 JC = 0
3772 JHEP = 0
3773 JD = 0
3774 ORG = 0
3775 IF(IST.EQ.155.AND.IST2.EQ.155) THEN
3776 IDM = IDHW(XHEP)
3777 ORG = ANTC
3778 IF(BVDEC1(XHEP).OR.BVDEC2(XHEP).OR.BVHRD(XHEP).OR.
3779 & BVHRD2(XHEP)) THEN
3780 JC=ANTC
3781 ID = IDHW(JC)
3782 JHEP = JC
3783 IF(BVDEC1(JC).OR.BVDEC2(JC)) THEN
3784 IF(IHEP.EQ.(NHEP+1)) ANTC=JDAHEP(1,JC)
3785 GOTO 200
3786 ENDIF
3787 IF (ID.EQ.449) THEN
3788C--SPECIAL FOR GLUINO DECAYS
3789 ID=IDHW(XHEP)
3790 IF(IHEP.EQ.NHEP+1) ID = 407
3791 CALL HWBRC1(JC,ID,JHEP,.FALSE.,*999)
3792 ELSE
3793 IF(IDHW(JDAHEP(1,JHEP)).EQ.449) THEN
3794 JC=JDAHEP(1,JHEP)
3795 ELSE
3796 JC=JDAHEP(2,JHEP)
3797 ENDIF
3798 ENDIF
3799C--SEARCH IN JET
3800 CALL HWBRC2(COLP,XHEP,JC,.FALSE.,BVVUSE,.FALSE.)
3801 ANTC = COLP
3802 IF(IHEP.LE.NHEP.AND.ACOLRD(IDHW(IHEP)).AND.
3803 & COLRD(IDHW(COLP)).AND.JMOHEP(2,COLP).EQ.0) THEN
3804 JMOHEP(2,COLP) = IHEP
3805 ELSEIF(IHEP.LE.NHEP.AND.IDHW(IHEP).LE.6.AND.
3806 & IDHW(COLP).LE.6.AND.JDAHEP(2,COLP).EQ.0) THEN
3807 JDAHEP(2,COLP) = IHEP
3808 ELSEIF(IHEP.GT.NHEP.AND.
3809 & ((BVHRD(XHEP).AND.COLRD(JDAHEP(1,XHEP))).
3810 & OR.(BVHRD2(XHEP).AND.ACOLRD(JDAHEP(2,XHEP)))).AND.
3811 & ACOLRD(IDHW(COLP)).AND.JDAHEP(2,COLP).EQ.0) THEN
3812 JDAHEP(2,COLP) = IHEP
3813 ENDIF
3814 ENDIF
3815 ENDIF
3816 200 CONTINUE
3817 IF(IHEP.EQ.NHEP+1) THEN
3818 IF(HRDCOL(1,1).NE.ANTC.AND.ANTC.NE.0) THEN
3819 HRDCOL(1,1)=ANTC
3820 IF(JDAHEP(2,ANTC).EQ.IHEP) THEN
3821 IF(JDAHEP(2,JMOHEP(1,HRDCOL(1,2))).EQ.JDAHEP(2,HRDCOL(1,2)).
3822 & AND.JMOHEP(2,JDAHEP(2,HRDCOL(1,2))).EQ.JMOHEP(1,HRDCOL(1,2)))
3823 & THEN
3824 JDAHEP(2,ANTC) = JMOHEP(2,HRDCOL(1,2))
3825 ELSE
3826 JDAHEP(2,ANTC) = JMOHEP(1,HRDCOL(1,2))
3827 ENDIF
3828 ELSEIF(JMOHEP(2,ANTC).EQ.IHEP) THEN
3829 JMOHEP(2,ANTC) = JMOHEP(1,HRDCOL(1,2))
3830 ENDIF
3831 ENDIF
3832 ELSEIF(IHEP.NE.1) THEN
3833 IF(JDAHEP(2,IHEP-1).NE.ANTC.AND.ANTC.NE.0) JDAHEP(2,IHEP-1)=ANTC
3834 ENDIF
3835 210 CONTINUE
3836C--Update BV decaying particles connections
3837 DO 310 IHEP=1,NHEP+1
3838 IF(IHEP.EQ.1) GOTO 310
3839 IF(IHEP.EQ.NHEP+1) THEN
3840 ANTC=HRDCOL(1,1)
3841 IF(ANTC.EQ.0.OR.IDHW(JDAHEP(1,HRDCOL(1,2))).LE.6) GOTO 310
3842 IST=155
3843 XHEP=HRDCOL(1,2)
3844 IF(ANTC.EQ.JDAHEP(2,XHEP)) ANTC=JDAHEP(1,JDAHEP(1,ANTC))
3845 ELSE
3846 ANTC=JMOHEP(2,IHEP)
3847 IST=ISTHEP(IHEP)
3848 IDM = IDHW(IHEP)
3849 XHEP=IHEP
3850 ENDIF
3851 IST2 = 0
3852 JC = 0
3853 JD = 0
3854 IF(ANTC.NE.0.AND.IHEP.NE.NHEP+1) THEN
3855