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