]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/jetset/lushow.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New functions and constructors added and some other fixes.
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lushow.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LUSHOW(IP1,IP2,QMAX)
5
6C...Purpose: to generate timelike parton showers from given partons.
7 IMPLICIT DOUBLE PRECISION(D)
8 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
9 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
10 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
11 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
12 DIMENSION PMTH(5,50),PS(5),PMA(4),PMSD(4),IEP(4),IPA(4),
13 &KFLA(4),KFLD(4),KFL(4),ITRY(4),ISI(4),ISL(4),DP(4),DPT(5,4),
14 &KSH(0:40),KCII(2),NIIS(2),IIIS(2,2),THEIIS(2,2),PHIIIS(2,2),
15 &ISII(2)
16
17C...Initialization of cutoff masses etc.
18 IF(MSTJ(41).LE.0.OR.(MSTJ(41).EQ.1.AND.QMAX.LE.PARJ(82)).OR.
19 &QMAX.LE.MIN(PARJ(82),PARJ(83))) RETURN
20 DO 100 IFL=0,40
21 KSH(IFL)=0
22 100 CONTINUE
23 KSH(21)=1
24 PMTH(1,21)=ULMASS(21)
25 PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25*PARJ(82)**2)
26 PMTH(3,21)=2.*PMTH(2,21)
27 PMTH(4,21)=PMTH(3,21)
28 PMTH(5,21)=PMTH(3,21)
29 PMTH(1,22)=ULMASS(22)
30 PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25*PARJ(83)**2)
31 PMTH(3,22)=2.*PMTH(2,22)
32 PMTH(4,22)=PMTH(3,22)
33 PMTH(5,22)=PMTH(3,22)
34 PMQTH1=PARJ(82)
35 IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83))
36 PMQTH2=PMTH(2,21)
37 IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
38 DO 110 IFL=1,8
39 KSH(IFL)=1
40 PMTH(1,IFL)=ULMASS(IFL)
41 PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25*PMQTH1**2)
42 PMTH(3,IFL)=PMTH(2,IFL)+PMQTH2
43 PMTH(4,IFL)=SQRT(PMTH(1,IFL)**2+0.25*PARJ(82)**2)+PMTH(2,21)
44 PMTH(5,IFL)=SQRT(PMTH(1,IFL)**2+0.25*PARJ(83)**2)+PMTH(2,22)
45 110 CONTINUE
46 DO 120 IFL=11,17,2
47 IF(MSTJ(41).GE.2) KSH(IFL)=1
48 PMTH(1,IFL)=ULMASS(IFL)
49 PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25*PARJ(83)**2)
50 PMTH(3,IFL)=PMTH(2,IFL)+PMTH(2,22)
51 PMTH(4,IFL)=PMTH(3,IFL)
52 PMTH(5,IFL)=PMTH(3,IFL)
53 120 CONTINUE
54 PT2MIN=MAX(0.5*PARJ(82),1.1*PARJ(81))**2
55 ALAMS=PARJ(81)**2
56 ALFM=LOG(PT2MIN/ALAMS)
57
58C...Store positions of shower initiating partons.
59 IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
60 NPA=1
61 IPA(1)=IP1
62 ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
63 &MSTU(32))) THEN
64 NPA=2
65 IPA(1)=IP1
66 IPA(2)=IP2
67 ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0
68 &.AND.IP2.GE.-3) THEN
69 NPA=IABS(IP2)
70 DO 130 I=1,NPA
71 IPA(I)=IP1+I-1
72 130 CONTINUE
73 ELSE
74 CALL LUERRM(12,
75 & '(LUSHOW:) failed to reconstruct showering system')
76 IF(MSTU(21).GE.1) RETURN
77 ENDIF
78
79C...Check on phase space available for emission.
80 IREJ=0
81 DO 140 J=1,5
82 PS(J)=0.
83 140 CONTINUE
84 PM=0.
85 DO 160 I=1,NPA
86 KFLA(I)=IABS(K(IPA(I),2))
87 PMA(I)=P(IPA(I),5)
88C...Special cutoff masses for t, l, h with variable masses.
89 IFLA=KFLA(I)
90 IF(KFLA(I).GE.6.AND.KFLA(I).LE.8) THEN
91 IFLA=37+KFLA(I)+ISIGN(2,K(IPA(I),2))
92 PMTH(1,IFLA)=PMA(I)
93 PMTH(2,IFLA)=SQRT(PMTH(1,IFLA)**2+0.25*PMQTH1**2)
94 PMTH(3,IFLA)=PMTH(2,IFLA)+PMQTH2
95 PMTH(4,IFLA)=SQRT(PMTH(1,IFLA)**2+0.25*PARJ(82)**2)+PMTH(2,21)
96 PMTH(5,IFLA)=SQRT(PMTH(1,IFLA)**2+0.25*PARJ(83)**2)+PMTH(2,22)
97 ENDIF
98 IF(KFLA(I).LE.40) THEN
99 IF(KSH(KFLA(I)).EQ.1) PMA(I)=PMTH(3,IFLA)
100 ENDIF
101 PM=PM+PMA(I)
102 IF(KFLA(I).GT.40) THEN
103 IREJ=IREJ+1
104 ELSE
105 IF(KSH(KFLA(I)).EQ.0.OR.PMA(I).GT.QMAX) IREJ=IREJ+1
106 ENDIF
107 DO 150 J=1,4
108 PS(J)=PS(J)+P(IPA(I),J)
109 150 CONTINUE
110 160 CONTINUE
111 IF(IREJ.EQ.NPA) RETURN
112 PS(5)=SQRT(MAX(0.,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
113 IF(NPA.EQ.1) PS(5)=PS(4)
114 IF(PS(5).LE.PM+PMQTH1) RETURN
115
116C...Check if 3-jet matrix elements to be used.
117 M3JC=0
118 IF(NPA.EQ.2.AND.MSTJ(47).GE.1) THEN
119 IF(KFLA(1).GE.1.AND.KFLA(1).LE.8.AND.KFLA(2).GE.1.AND.
120 & KFLA(2).LE.8) M3JC=1
121 IF((KFLA(1).EQ.11.OR.KFLA(1).EQ.13.OR.KFLA(1).EQ.15.OR.
122 & KFLA(1).EQ.17).AND.KFLA(2).EQ.KFLA(1)) M3JC=1
123 IF((KFLA(1).EQ.11.OR.KFLA(1).EQ.13.OR.KFLA(1).EQ.15.OR.
124 & KFLA(1).EQ.17).AND.KFLA(2).EQ.KFLA(1)+1) M3JC=1
125 IF((KFLA(1).EQ.12.OR.KFLA(1).EQ.14.OR.KFLA(1).EQ.16.OR.
126 & KFLA(1).EQ.18).AND.KFLA(2).EQ.KFLA(1)-1) M3JC=1
127 IF(MSTJ(47).EQ.2.OR.MSTJ(47).EQ.4) M3JC=1
128 M3JCM=0
129 IF(M3JC.EQ.1.AND.MSTJ(47).GE.3.AND.KFLA(1).EQ.KFLA(2)) THEN
130 M3JCM=1
131 QME=(2.*PMTH(1,KFLA(1))/PS(5))**2
132 ENDIF
133 ENDIF
134
135C...Find if interference with initial state partons.
136 MIIS=0
137 IF(MSTJ(50).GE.1.AND.MSTJ(50).LE.3.AND.NPA.EQ.2) MIIS=MSTJ(50)
138 IF(MIIS.NE.0) THEN
139 DO 180 I=1,2
140 KCII(I)=0
141 KCA=LUCOMP(KFLA(I))
142 IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2))
143 NIIS(I)=0
144 IF(KCII(I).NE.0) THEN
145 DO 170 J=1,2
146 ICSI=MOD(K(IPA(I),3+J)/MSTU(5),MSTU(5))
147 IF(ICSI.GT.0.AND.ICSI.NE.IPA(1).AND.ICSI.NE.IPA(2).AND.
148 & (KCII(I).EQ.(-1)**(J+1).OR.KCII(I).EQ.2)) THEN
149 NIIS(I)=NIIS(I)+1
150 IIIS(I,NIIS(I))=ICSI
151 ENDIF
152 170 CONTINUE
153 ENDIF
154 180 CONTINUE
155 IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0
156 ENDIF
157
158C...Boost interfering initial partons to rest frame
159C...and reconstruct their polar and azimuthal angles.
160 IF(MIIS.NE.0) THEN
161 DO 200 I=1,2
162 DO 190 J=1,5
163 K(N+I,J)=K(IPA(I),J)
164 P(N+I,J)=P(IPA(I),J)
165 V(N+I,J)=0.
166 190 CONTINUE
167 200 CONTINUE
168 DO 220 I=3,2+NIIS(1)
169 DO 210 J=1,5
170 K(N+I,J)=K(IIIS(1,I-2),J)
171 P(N+I,J)=P(IIIS(1,I-2),J)
172 V(N+I,J)=0.
173 210 CONTINUE
174 220 CONTINUE
175 DO 240 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
176 DO 230 J=1,5
177 K(N+I,J)=K(IIIS(2,I-2-NIIS(1)),J)
178 P(N+I,J)=P(IIIS(2,I-2-NIIS(1)),J)
179 V(N+I,J)=0.
180 230 CONTINUE
181 240 CONTINUE
182 CALL LUDBRB(N+1,N+2+NIIS(1)+NIIS(2),0.,0.,-DBLE(PS(1)/PS(4)),
183 & -DBLE(PS(2)/PS(4)),-DBLE(PS(3)/PS(4)))
184 PHI=ULANGL(P(N+1,1),P(N+1,2))
185 CALL LUDBRB(N+1,N+2+NIIS(1)+NIIS(2),0.,-PHI,0D0,0D0,0D0)
186 THE=ULANGL(P(N+1,3),P(N+1,1))
187 CALL LUDBRB(N+1,N+2+NIIS(1)+NIIS(2),-THE,0.,0D0,0D0,0D0)
188 DO 250 I=3,2+NIIS(1)
189 THEIIS(1,I-2)=ULANGL(P(N+I,3),SQRT(P(N+I,1)**2+P(N+I,2)**2))
190 PHIIIS(1,I-2)=ULANGL(P(N+I,1),P(N+I,2))
191 250 CONTINUE
192 DO 260 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
193 THEIIS(2,I-2-NIIS(1))=PARU(1)-ULANGL(P(N+I,3),
194 & SQRT(P(N+I,1)**2+P(N+I,2)**2))
195 PHIIIS(2,I-2-NIIS(1))=ULANGL(P(N+I,1),P(N+I,2))
196 260 CONTINUE
197 ENDIF
198
199C...Define imagined single initiator of shower for parton system.
200 NS=N
201 IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
202 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
203 IF(MSTU(21).GE.1) RETURN
204 ENDIF
205 IF(NPA.GE.2) THEN
206 K(N+1,1)=11
207 K(N+1,2)=21
208 K(N+1,3)=0
209 K(N+1,4)=0
210 K(N+1,5)=0
211 P(N+1,1)=0.
212 P(N+1,2)=0.
213 P(N+1,3)=0.
214 P(N+1,4)=PS(5)
215 P(N+1,5)=PS(5)
216 V(N+1,5)=PS(5)**2
217 N=N+1
218 ENDIF
219
220C...Loop over partons that may branch.
221 NEP=NPA
222 IM=NS
223 IF(NPA.EQ.1) IM=NS-1
224 270 IM=IM+1
225 IF(N.GT.NS) THEN
226 IF(IM.GT.N) GOTO 510
227 KFLM=IABS(K(IM,2))
228 IF(KFLM.GT.40) GOTO 270
229 IF(KSH(KFLM).EQ.0) GOTO 270
230 IFLM=KFLM
231 IF(KFLM.GE.6.AND.KFLM.LE.8) IFLM=37+KFLM+ISIGN(2,K(IM,2))
232 IF(P(IM,5).LT.PMTH(2,IFLM)) GOTO 270
233 IGM=K(IM,3)
234 ELSE
235 IGM=-1
236 ENDIF
237 IF(N+NEP.GT.MSTU(4)-MSTU(32)-5) THEN
238 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
239 IF(MSTU(21).GE.1) RETURN
240 ENDIF
241
242C...Position of aunt (sister to branching parton).
243C...Origin and flavour of daughters.
244 IAU=0
245 IF(IGM.GT.0) THEN
246 IF(K(IM-1,3).EQ.IGM) IAU=IM-1
247 IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
248 ENDIF
249 IF(IGM.GE.0) THEN
250 K(IM,4)=N+1
251 DO 280 I=1,NEP
252 K(N+I,3)=IM
253 280 CONTINUE
254 ELSE
255 K(N+1,3)=IPA(1)
256 ENDIF
257 IF(IGM.LE.0) THEN
258 DO 290 I=1,NEP
259 K(N+I,2)=K(IPA(I),2)
260 290 CONTINUE
261 ELSEIF(KFLM.NE.21) THEN
262 K(N+1,2)=K(IM,2)
263 K(N+2,2)=K(IM,5)
264 ELSEIF(K(IM,5).EQ.21) THEN
265 K(N+1,2)=21
266 K(N+2,2)=21
267 ELSE
268 K(N+1,2)=K(IM,5)
269 K(N+2,2)=-K(IM,5)
270 ENDIF
271
272C...Reset flags on daughers and tries made.
273 DO 300 IP=1,NEP
274 K(N+IP,1)=3
275 K(N+IP,4)=0
276 K(N+IP,5)=0
277 KFLD(IP)=IABS(K(N+IP,2))
278 IF(KCHG(LUCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1
279 ITRY(IP)=0
280 ISL(IP)=0
281 ISI(IP)=0
282 IF(KFLD(IP).LE.40) THEN
283 IF(KSH(KFLD(IP)).EQ.1) ISI(IP)=1
284 ENDIF
285 300 CONTINUE
286 ISLM=0
287
288C...Maximum virtuality of daughters.
289 IF(IGM.LE.0) THEN
290 DO 310 I=1,NPA
291 IF(NPA.GE.3) P(N+I,4)=(PS(4)*P(IPA(I),4)-PS(1)*P(IPA(I),1)-
292 & PS(2)*P(IPA(I),2)-PS(3)*P(IPA(I),3))/PS(5)
293 P(N+I,5)=MIN(QMAX,PS(5))
294 IF(NPA.GE.3) P(N+I,5)=MIN(P(N+I,5),P(N+I,4))
295 IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
296 310 CONTINUE
297 ELSE
298 IF(MSTJ(43).LE.2) PEM=V(IM,2)
299 IF(MSTJ(43).GE.3) PEM=P(IM,4)
300 P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
301 P(N+2,5)=MIN(P(IM,5),(1.-V(IM,1))*PEM)
302 IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
303 ENDIF
304 DO 320 I=1,NEP
305 PMSD(I)=P(N+I,5)
306 IF(ISI(I).EQ.1) THEN
307 IFLD=KFLD(I)
308 IF(KFLD(I).GE.6.AND.KFLD(I).LE.8) IFLD=37+KFLD(I)+
309 & ISIGN(2,K(N+I,2))
310 IF(P(N+I,5).LE.PMTH(3,IFLD)) P(N+I,5)=PMTH(1,IFLD)
311 ENDIF
312 V(N+I,5)=P(N+I,5)**2
313 320 CONTINUE
314
315C...Choose one of the daughters for evolution.
316 330 INUM=0
317 IF(NEP.EQ.1) INUM=1
318 DO 340 I=1,NEP
319 IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
320 340 CONTINUE
321 DO 350 I=1,NEP
322 IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
323 IFLD=KFLD(I)
324 IF(KFLD(I).GE.6.AND.KFLD(I).LE.8) IFLD=37+KFLD(I)+
325 & ISIGN(2,K(N+I,2))
326 IF(P(N+I,5).GE.PMTH(2,IFLD)) INUM=I
327 ENDIF
328 350 CONTINUE
329 IF(INUM.EQ.0) THEN
330 RMAX=0.
331 DO 360 I=1,NEP
332 IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQTH2) THEN
333 RPM=P(N+I,5)/PMSD(I)
334 IFLD=KFLD(I)
335 IF(KFLD(I).GE.6.AND.KFLD(I).LE.8) IFLD=37+KFLD(I)+
336 & ISIGN(2,K(N+I,2))
337 IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IFLD)) THEN
338 RMAX=RPM
339 INUM=I
340 ENDIF
341 ENDIF
342 360 CONTINUE
343 ENDIF
344
345C...Store information on choice of evolving daughter.
346 INUM=MAX(1,INUM)
347 IEP(1)=N+INUM
348 DO 370 I=2,NEP
349 IEP(I)=IEP(I-1)+1
350 IF(IEP(I).GT.N+NEP) IEP(I)=N+1
351 370 CONTINUE
352 DO 380 I=1,NEP
353 KFL(I)=IABS(K(IEP(I),2))
354 380 CONTINUE
355 ITRY(INUM)=ITRY(INUM)+1
356 IF(ITRY(INUM).GT.200) THEN
357 CALL LUERRM(14,'(LUSHOW:) caught in infinite loop')
358 IF(MSTU(21).GE.1) RETURN
359 ENDIF
360 Z=0.5
361 IF(KFL(1).GT.40) GOTO 430
362 IF(KSH(KFL(1)).EQ.0) GOTO 430
363 IFL=KFL(1)
364 IF(KFL(1).GE.6.AND.KFL(1).LE.8) IFL=37+KFL(1)+
365 &ISIGN(2,K(IEP(1),2))
366 IF(P(IEP(1),5).LT.PMTH(2,IFL)) GOTO 430
367
368C...Select side for interference with initial state partons.
369 IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN
370 III=IEP(1)-NS-1
371 ISII(III)=0
372 IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN
373 ISII(III)=1
374 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.1) THEN
375 IF(RLU(0).GT.0.5) ISII(III)=1
376 ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.2) THEN
377 ISII(III)=1
378 IF(RLU(0).GT.0.5) ISII(III)=2
379 ENDIF
380 ENDIF
381
382C...Calculate allowed z range.
383 IF(NEP.EQ.1) THEN
384 PMED=PS(4)
385 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
386 PMED=P(IM,5)
387 ELSE
388 IF(INUM.EQ.1) PMED=V(IM,1)*PEM
389 IF(INUM.EQ.2) PMED=(1.-V(IM,1))*PEM
390 ENDIF
391 IF(MOD(MSTJ(43),2).EQ.1) THEN
392 ZC=PMTH(2,21)/PMED
393 ZCE=PMTH(2,22)/PMED
394 ELSE
395 ZC=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,21)/PMED)**2)))
396 IF(ZC.LT.1E-4) ZC=(PMTH(2,21)/PMED)**2
397 ZCE=0.5*(1.-SQRT(MAX(0.,1.-(2.*PMTH(2,22)/PMED)**2)))
398 IF(ZCE.LT.1E-4) ZCE=(PMTH(2,22)/PMED)**2
399 ENDIF
400 ZC=MIN(ZC,0.491)
401 ZCE=MIN(ZCE,0.491)
402 IF((MSTJ(41).EQ.1.AND.ZC.GT.0.49).OR.(MSTJ(41).GE.2.AND.
403 &MIN(ZC,ZCE).GT.0.49)) THEN
404 P(IEP(1),5)=PMTH(1,IFL)
405 V(IEP(1),5)=P(IEP(1),5)**2
406 GOTO 430
407 ENDIF
408
409C...Integral of Altarelli-Parisi z kernel for QCD.
410 IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
411 FBR=6.*LOG((1.-ZC)/ZC)+MSTJ(45)*(0.5-ZC)
412 ELSEIF(MSTJ(49).EQ.0) THEN
413 FBR=(8./3.)*LOG((1.-ZC)/ZC)
414
415C...Integral of Altarelli-Parisi z kernel for scalar gluon.
416 ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
417 FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1.-2.*ZC)
418 ELSEIF(MSTJ(49).EQ.1) THEN
419 FBR=(1.-2.*ZC)/3.
420 IF(IGM.EQ.0.AND.M3JC.EQ.1) FBR=4.*FBR
421
422C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
423 ELSEIF(KFL(1).EQ.21) THEN
424 FBR=6.*MSTJ(45)*(0.5-ZC)
425 ELSE
426 FBR=2.*LOG((1.-ZC)/ZC)
427 ENDIF
428
429C...Reset QCD probability for lepton.
430 IF(KFL(1).GE.11.AND.KFL(1).LE.18) FBR=0.
431
432C...Integral of Altarelli-Parisi kernel for photon emission.
433 IF(MSTJ(41).GE.2.AND.KFL(1).GE.1.AND.KFL(1).LE.18) THEN
434 FBRE=(KCHG(KFL(1),1)/3.)**2*2.*LOG((1.-ZCE)/ZCE)
435 IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE
436 ENDIF
437
438C...Inner veto algorithm starts. Find maximum mass for evolution.
439 390 PMS=V(IEP(1),5)
440 IF(IGM.GE.0) THEN
441 PM2=0.
442 DO 400 I=2,NEP
443 PM=P(IEP(I),5)
444 IF(KFL(I).LE.40) THEN
445 IFLI=KFL(I)
446 IF(KFL(I).GE.6.AND.KFL(I).LE.8) IFLI=37+KFL(I)+
447 & ISIGN(2,K(IEP(I),2))
448 IF(KSH(KFL(I)).EQ.1) PM=PMTH(2,IFLI)
449 ENDIF
450 PM2=PM2+PM
451 400 CONTINUE
452 PMS=MIN(PMS,(P(IM,5)-PM2)**2)
453 ENDIF
454
455C...Select mass for daughter in QCD evolution.
456 B0=27./6.
457 DO 410 IFF=4,MSTJ(45)
458 IF(PMS.GT.4.*PMTH(2,IFF)**2) B0=(33.-2.*IFF)/6.
459 410 CONTINUE
460 IF(FBR.LT.1E-3) THEN
461 PMSQCD=0.
462 ELSEIF(MSTJ(44).LE.0) THEN
463 PMSQCD=PMS*EXP(MAX(-50.,LOG(RLU(0))*PARU(2)/(PARU(111)*FBR)))
464 ELSEIF(MSTJ(44).EQ.1) THEN
465 PMSQCD=4.*ALAMS*(0.25*PMS/ALAMS)**(RLU(0)**(B0/FBR))
466 ELSE
467 PMSQCD=PMS*EXP(MAX(-50.,ALFM*B0*LOG(RLU(0))/FBR))
468 ENDIF
469 IF(ZC.GT.0.49.OR.PMSQCD.LE.PMTH(4,IFL)**2) PMSQCD=PMTH(2,IFL)**2
470 V(IEP(1),5)=PMSQCD
471 MCE=1
472
473C...Select mass for daughter in QED evolution.
474 IF(MSTJ(41).GE.2.AND.KFL(1).GE.1.AND.KFL(1).LE.18) THEN
475 PMSQED=PMS*EXP(MAX(-50.,LOG(RLU(0))*PARU(2)/(PARU(101)*FBRE)))
476 IF(ZCE.GT.0.49.OR.PMSQED.LE.PMTH(5,IFL)**2) PMSQED=
477 & PMTH(2,IFL)**2
478 IF(PMSQED.GT.PMSQCD) THEN
479 V(IEP(1),5)=PMSQED
480 MCE=2
481 ENDIF
482 ENDIF
483
484C...Check whether daughter mass below cutoff.
485 P(IEP(1),5)=SQRT(V(IEP(1),5))
486 IF(P(IEP(1),5).LE.PMTH(3,IFL)) THEN
487 P(IEP(1),5)=PMTH(1,IFL)
488 V(IEP(1),5)=P(IEP(1),5)**2
489 GOTO 430
490 ENDIF
491
492C...Select z value of branching: q -> qgamma.
493 IF(MCE.EQ.2) THEN
494 Z=1.-(1.-ZCE)*(ZCE/(1.-ZCE))**RLU(0)
495 IF(1.+Z**2.LT.2.*RLU(0)) GOTO 390
496 K(IEP(1),5)=22
497
498C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
499 ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
500 Z=1.-(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
501 IF(1.+Z**2.LT.2.*RLU(0)) GOTO 390
502 K(IEP(1),5)=21
503 ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*(0.5-ZC).LT.RLU(0)*FBR) THEN
504 Z=(1.-ZC)*(ZC/(1.-ZC))**RLU(0)
505 IF(RLU(0).GT.0.5) Z=1.-Z
506 IF((1.-Z*(1.-Z))**2.LT.RLU(0)) GOTO 390
507 K(IEP(1),5)=21
508 ELSEIF(MSTJ(49).NE.1) THEN
509 Z=ZC+(1.-2.*ZC)*RLU(0)
510 IF(Z**2+(1.-Z)**2.LT.RLU(0)) GOTO 390
511 KFLB=1+INT(MSTJ(45)*RLU(0))
512 PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
513 IF(PMQ.GE.1.) GOTO 390
514 PMQ0=4.*PMTH(2,21)**2/V(IEP(1),5)
515 IF(MOD(MSTJ(43),2).EQ.0.AND.(1.+0.5*PMQ)*SQRT(1.-PMQ).LT.
516 & RLU(0)*(1.+0.5*PMQ0)*SQRT(1.-PMQ0)) GOTO 390
517 K(IEP(1),5)=KFLB
518
519C...Ditto for scalar gluon model.
520 ELSEIF(KFL(1).NE.21) THEN
521 Z=1.-SQRT(ZC**2+RLU(0)*(1.-2.*ZC))
522 K(IEP(1),5)=21
523 ELSEIF(RLU(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
524 Z=ZC+(1.-2.*ZC)*RLU(0)
525 K(IEP(1),5)=21
526 ELSE
527 Z=ZC+(1.-2.*ZC)*RLU(0)
528 KFLB=1+INT(MSTJ(45)*RLU(0))
529 PMQ=4.*PMTH(2,KFLB)**2/V(IEP(1),5)
530 IF(PMQ.GE.1.) GOTO 390
531 K(IEP(1),5)=KFLB
532 ENDIF
533 IF(MCE.EQ.1.AND.MSTJ(44).GE.2) THEN
534 IF(Z*(1.-Z)*V(IEP(1),5).LT.PT2MIN) GOTO 390
535 IF(ALFM/LOG(V(IEP(1),5)*Z*(1.-Z)/ALAMS).LT.RLU(0)) GOTO 390
536 ENDIF
537
538C...Check if z consistent with chosen m.
539 IF(KFL(1).EQ.21) THEN
540 KFLGD1=IABS(K(IEP(1),5))
541 KFLGD2=KFLGD1
542 ELSE
543 KFLGD1=KFL(1)
544 KFLGD2=IABS(K(IEP(1),5))
545 ENDIF
546 IF(NEP.EQ.1) THEN
547 PED=PS(4)
548 ELSEIF(NEP.GE.3) THEN
549 PED=P(IEP(1),4)
550 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
551 PED=0.5*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
552 ELSE
553 IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
554 IF(IEP(1).EQ.N+2) PED=(1.-V(IM,1))*PEM
555 ENDIF
556 IF(MOD(MSTJ(43),2).EQ.1) THEN
557 IFLGD1=KFLGD1
558 IF(KFLGD1.GE.6.AND.KFLGD1.LE.8) IFLGD1=IFL
559 PMQTH3=0.5*PARJ(82)
560 IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
561 PMQ1=(PMTH(1,IFLGD1)**2+PMQTH3**2)/V(IEP(1),5)
562 PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(IEP(1),5)
563 ZD=SQRT(MAX(0.,(1.-V(IEP(1),5)/PED**2)*((1.-PMQ1-PMQ2)**2-
564 & 4.*PMQ1*PMQ2)))
565 ZH=1.+PMQ1-PMQ2
566 ELSE
567 ZD=SQRT(MAX(0.,1.-V(IEP(1),5)/PED**2))
568 ZH=1.
569 ENDIF
570 ZL=0.5*(ZH-ZD)
571 ZU=0.5*(ZH+ZD)
572 IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 390
573 IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*
574 &(1.-ZU)))
575 IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
576
577C...Width suppression for q -> q + g.
578 IF(MSTJ(40).NE.0.AND.KFL(1).NE.21) THEN
579 IF(IGM.EQ.0) THEN
580 EGLU=0.5*PS(5)*(1.-Z)*(1.+V(IEP(1),5)/V(NS+1,5))
581 ELSE
582 EGLU=PMED*(1.-Z)
583 ENDIF
584 CHI=PARJ(89)**2/(PARJ(89)**2+EGLU**2)
585 IF(MSTJ(40).EQ.1) THEN
586 IF(CHI.LT.RLU(0)) GOTO 390
587 ELSEIF(MSTJ(40).EQ.2) THEN
588 IF(1.-CHI.LT.RLU(0)) GOTO 390
589 ENDIF
590 ENDIF
591
592C...Three-jet matrix element correction.
593 IF(IGM.EQ.0.AND.M3JC.EQ.1) THEN
594 X1=Z*(1.+V(IEP(1),5)/V(NS+1,5))
595 X2=1.-V(IEP(1),5)/V(NS+1,5)
596 X3=(1.-X1)+(1.-X2)
597 IF(MCE.EQ.2) THEN
598 KI1=K(IPA(INUM),2)
599 KI2=K(IPA(3-INUM),2)
600 QF1=KCHG(IABS(KI1),1)*ISIGN(1,KI1)/3.
601 QF2=KCHG(IABS(KI2),1)*ISIGN(1,KI2)/3.
602 WSHOW=QF1**2*(1.-X1)/X3*(1.+(X1/(2.-X2))**2)+
603 & QF2**2*(1.-X2)/X3*(1.+(X2/(2.-X1))**2)
604 WME=(QF1*(1.-X1)/X3-QF2*(1.-X2)/X3)**2*(X1**2+X2**2)
605 ELSEIF(MSTJ(49).NE.1) THEN
606 WSHOW=1.+(1.-X1)/X3*(X1/(2.-X2))**2+
607 & (1.-X2)/X3*(X2/(2.-X1))**2
608 WME=X1**2+X2**2
609 IF(M3JCM.EQ.1) WME=WME-QME*X3-0.5*QME**2-
610 & (0.5*QME+0.25*QME**2)*((1.-X2)/MAX(1E-7,1.-X1)+
611 & (1.-X1)/MAX(1E-7,1.-X2))
612 ELSE
613 WSHOW=4.*X3*((1.-X1)/(2.-X2)**2+(1.-X2)/(2.-X1)**2)
614 WME=X3**2
615 IF(MSTJ(102).GE.2) WME=X3**2-2.*(1.+X3)*(1.-X1)*(1.-X2)*
616 & PARJ(171)
617 ENDIF
618 IF(WME.LT.RLU(0)*WSHOW) GOTO 390
619
620C...Impose angular ordering by rejection of nonordered emission.
621 ELSEIF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2) THEN
622 MAOM=1
623 ZM=V(IM,1)
624 IF(IEP(1).EQ.N+2) ZM=1.-V(IM,1)
625 THE2ID=Z*(1.-Z)*(ZM*P(IM,4))**2/V(IEP(1),5)
626 IAOM=IM
627 420 IF(K(IAOM,5).EQ.22) THEN
628 IAOM=K(IAOM,3)
629 IF(K(IAOM,3).LE.NS) MAOM=0
630 IF(MAOM.EQ.1) GOTO 420
631 ENDIF
632 IF(MAOM.EQ.1) THEN
633 THE2IM=V(IAOM,1)*(1.-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
634 IF(THE2ID.LT.THE2IM) GOTO 390
635 ENDIF
636 ENDIF
637
638C...Impose user-defined maximum angle at first branching.
639 IF(MSTJ(48).EQ.1) THEN
640 IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
641 THE2ID=Z*(1.-Z)*PS(4)**2/V(IEP(1),5)
642 IF(THE2ID.LT.1./PARJ(85)**2) GOTO 390
643 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
644 THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
645 IF(THE2ID.LT.1./PARJ(85)**2) GOTO 390
646 ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
647 THE2ID=Z*(1.-Z)*(0.5*P(IM,4))**2/V(IEP(1),5)
648 IF(THE2ID.LT.1./PARJ(86)**2) GOTO 390
649 ENDIF
650 ENDIF
651
652C...Impose angular constraint in first branching from interference
653C...with initial state partons.
654 IF(MIIS.GE.2.AND.IEP(1).LE.NS+3) THEN
655 THE2D=MAX((1.-Z)/Z,Z/(1.-Z))*V(IEP(1),5)/(0.5*P(IM,4))**2
656 IF(IEP(1).EQ.NS+2.AND.ISII(1).GE.1) THEN
657 IF(THE2D.GT.THEIIS(1,ISII(1))**2) GOTO 390
658 ELSEIF(IEP(1).EQ.NS+3.AND.ISII(2).GE.1) THEN
659 IF(THE2D.GT.THEIIS(2,ISII(2))**2) GOTO 390
660 ENDIF
661 ENDIF
662
663C...End of inner veto algorithm. Check if only one leg evolved so far.
664 430 V(IEP(1),1)=Z
665 ISL(1)=0
666 ISL(2)=0
667 IF(NEP.EQ.1) GOTO 460
668 IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 330
669 DO 440 I=1,NEP
670 IF(ITRY(I).EQ.0.AND.KFLD(I).LE.40) THEN
671 IF(KSH(KFLD(I)).EQ.1) THEN
672 IFLD=KFLD(I)
673 IF(KFLD(I).GE.6.AND.KFLD(I).LE.8) IFLD=37+KFLD(I)+
674 & ISIGN(2,K(N+I,2))
675 IF(P(N+I,5).GE.PMTH(2,IFLD)) GOTO 330
676 ENDIF
677 ENDIF
678 440 CONTINUE
679
680C...Check if chosen multiplet m1,m2,z1,z2 is physical.
681 IF(NEP.EQ.3) THEN
682 PA1S=(P(N+1,4)+P(N+1,5))*(P(N+1,4)-P(N+1,5))
683 PA2S=(P(N+2,4)+P(N+2,5))*(P(N+2,4)-P(N+2,5))
684 PA3S=(P(N+3,4)+P(N+3,5))*(P(N+3,4)-P(N+3,5))
685 PTS=0.25*(2.*PA1S*PA2S+2.*PA1S*PA3S+2.*PA2S*PA3S-
686 & PA1S**2-PA2S**2-PA3S**2)/PA1S
687 IF(PTS.LE.0.) GOTO 330
688 ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
689 DO 450 I1=N+1,N+2
690 KFLDA=IABS(K(I1,2))
691 IF(KFLDA.GT.40) GOTO 450
692 IF(KSH(KFLDA).EQ.0) GOTO 450
693 IFLDA=KFLDA
694 IF(KFLDA.GE.6.AND.KFLDA.LE.8) IFLDA=37+KFLDA+
695 & ISIGN(2,K(I1,2))
696 IF(P(I1,5).LT.PMTH(2,IFLDA)) GOTO 450
697 IF(KFLDA.EQ.21) THEN
698 KFLGD1=IABS(K(I1,5))
699 KFLGD2=KFLGD1
700 ELSE
701 KFLGD1=KFLDA
702 KFLGD2=IABS(K(I1,5))
703 ENDIF
704 I2=2*N+3-I1
705 IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
706 PED=0.5*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
707 ELSE
708 IF(I1.EQ.N+1) ZM=V(IM,1)
709 IF(I1.EQ.N+2) ZM=1.-V(IM,1)
710 PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
711 & 4.*V(N+1,5)*V(N+2,5))
712 PED=PEM*(0.5*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/V(IM,5)
713 ENDIF
714 IF(MOD(MSTJ(43),2).EQ.1) THEN
715 PMQTH3=0.5*PARJ(82)
716 IF(KFLGD2.EQ.22) PMQTH3=0.5*PARJ(83)
717 IFLGD1=KFLGD1
718 IF(KFLGD1.GE.6.AND.KFLGD1.LE.8) IFLGD1=IFLDA
719 PMQ1=(PMTH(1,IFLGD1)**2+PMQTH3**2)/V(I1,5)
720 PMQ2=(PMTH(1,KFLGD2)**2+PMQTH3**2)/V(I1,5)
721 ZD=SQRT(MAX(0.,(1.-V(I1,5)/PED**2)*((1.-PMQ1-PMQ2)**2-
722 & 4.*PMQ1*PMQ2)))
723 ZH=1.+PMQ1-PMQ2
724 ELSE
725 ZD=SQRT(MAX(0.,1.-V(I1,5)/PED**2))
726 ZH=1.
727 ENDIF
728 ZL=0.5*(ZH-ZD)
729 ZU=0.5*(ZH+ZD)
730 IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(1)=1
731 IF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU)) ISL(2)=1
732 IF(KFLDA.EQ.21) V(I1,4)=LOG(ZU*(1.-ZL)/MAX(1E-20,ZL*(1.-ZU)))
733 IF(KFLDA.NE.21) V(I1,4)=LOG((1.-ZL)/MAX(1E-10,1.-ZU))
734 450 CONTINUE
735 IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
736 ISL(3-ISLM)=0
737 ISLM=3-ISLM
738 ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
739 ZDR1=MAX(0.,V(N+1,3)/MAX(1E-6,V(N+1,4))-1.)
740 ZDR2=MAX(0.,V(N+2,3)/MAX(1E-6,V(N+2,4))-1.)
741 IF(ZDR2.GT.RLU(0)*(ZDR1+ZDR2)) ISL(1)=0
742 IF(ISL(1).EQ.1) ISL(2)=0
743 IF(ISL(1).EQ.0) ISLM=1
744 IF(ISL(2).EQ.0) ISLM=2
745 ENDIF
746 IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 330
747 ENDIF
748 IFLD1=KFLD(1)
749 IF(KFLD(1).GE.6.AND.KFLD(1).LE.8) IFLD1=37+KFLD(1)+
750 &ISIGN(2,K(N+1,2))
751 IFLD2=KFLD(2)
752 IF(KFLD(2).GE.6.AND.KFLD(2).LE.8) IFLD2=37+KFLD(2)+
753 &ISIGN(2,K(N+2,2))
754 IF(IGM.GT.0.AND.MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
755 &PMTH(2,IFLD1).OR.P(N+2,5).GE.PMTH(2,IFLD2))) THEN
756 PMQ1=V(N+1,5)/V(IM,5)
757 PMQ2=V(N+2,5)/V(IM,5)
758 ZD=SQRT(MAX(0.,(1.-V(IM,5)/PEM**2)*((1.-PMQ1-PMQ2)**2-
759 & 4.*PMQ1*PMQ2)))
760 ZH=1.+PMQ1-PMQ2
761 ZL=0.5*(ZH-ZD)
762 ZU=0.5*(ZH+ZD)
763 IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 330
764 ENDIF
765
766C...Accepted branch. Construct four-momentum for initial partons.
767 460 MAZIP=0
768 MAZIC=0
769 IF(NEP.EQ.1) THEN
770 P(N+1,1)=0.
771 P(N+1,2)=0.
772 P(N+1,3)=SQRT(MAX(0.,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
773 & P(N+1,5))))
774 P(N+1,4)=P(IPA(1),4)
775 V(N+1,2)=P(N+1,4)
776 ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
777 PED1=0.5*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
778 P(N+1,1)=0.
779 P(N+1,2)=0.
780 P(N+1,3)=SQRT(MAX(0.,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
781 P(N+1,4)=PED1
782 P(N+2,1)=0.
783 P(N+2,2)=0.
784 P(N+2,3)=-P(N+1,3)
785 P(N+2,4)=P(IM,5)-PED1
786 V(N+1,2)=P(N+1,4)
787 V(N+2,2)=P(N+2,4)
788 ELSEIF(NEP.EQ.3) THEN
789 P(N+1,1)=0.
790 P(N+1,2)=0.
791 P(N+1,3)=SQRT(MAX(0.,PA1S))
792 P(N+2,1)=SQRT(PTS)
793 P(N+2,2)=0.
794 P(N+2,3)=0.5*(PA3S-PA2S-PA1S)/P(N+1,3)
795 P(N+3,1)=-P(N+2,1)
796 P(N+3,2)=0.
797 P(N+3,3)=-(P(N+1,3)+P(N+2,3))
798 V(N+1,2)=P(N+1,4)
799 V(N+2,2)=P(N+2,4)
800 V(N+3,2)=P(N+3,4)
801
802C...Construct transverse momentum for ordinary branching in shower.
803 ELSE
804 ZM=V(IM,1)
805 PZM=SQRT(MAX(0.,(PEM+P(IM,5))*(PEM-P(IM,5))))
806 PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4.*V(N+1,5)*V(N+2,5)
807 IF(PZM.LE.0.) THEN
808 PTS=0.
809 ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
810 PTS=(PEM**2*(ZM*(1.-ZM)*V(IM,5)-(1.-ZM)*V(N+1,5)-
811 & ZM*V(N+2,5))-0.25*PMLS)/PZM**2
812 ELSE
813 PTS=PMLS*(ZM*(1.-ZM)*PEM**2/V(IM,5)-0.25)/PZM**2
814 ENDIF
815 PT=SQRT(MAX(0.,PTS))
816
817C...Find coefficient of azimuthal asymmetry due to gluon polarization.
818 HAZIP=0.
819 IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21.
820 & AND.IAU.NE.0) THEN
821 IF(K(IGM,3).NE.0) MAZIP=1
822 ZAU=V(IGM,1)
823 IF(IAU.EQ.IM+1) ZAU=1.-V(IGM,1)
824 IF(MAZIP.EQ.0) ZAU=0.
825 IF(K(IGM,2).NE.21) THEN
826 HAZIP=2.*ZAU/(1.+ZAU**2)
827 ELSE
828 HAZIP=(ZAU/(1.-ZAU*(1.-ZAU)))**2
829 ENDIF
830 IF(K(N+1,2).NE.21) THEN
831 HAZIP=HAZIP*(-2.*ZM*(1.-ZM))/(1.-2.*ZM*(1.-ZM))
832 ELSE
833 HAZIP=HAZIP*(ZM*(1.-ZM)/(1.-ZM*(1.-ZM)))**2
834 ENDIF
835 ENDIF
836
837C...Find coefficient of azimuthal asymmetry due to soft gluon
838C...interference.
839 HAZIC=0.
840 IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.
841 & K(N+2,2).EQ.21).AND.IAU.NE.0) THEN
842 IF(K(IGM,3).NE.0) MAZIC=N+1
843 IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
844 IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
845 & ZM.GT.0.5) MAZIC=N+2
846 IF(K(IAU,2).EQ.22) MAZIC=0
847 ZS=ZM
848 IF(MAZIC.EQ.N+2) ZS=1.-ZM
849 ZGM=V(IGM,1)
850 IF(IAU.EQ.IM-1) ZGM=1.-V(IGM,1)
851 IF(MAZIC.EQ.0) ZGM=1.
852 IF(MAZIC.NE.0) HAZIC=(P(IM,5)/P(IGM,5))*
853 & SQRT((1.-ZS)*(1.-ZGM)/(ZS*ZGM))
854 HAZIC=MIN(0.95,HAZIC)
855 ENDIF
856 ENDIF
857
858C...Construct kinematics for ordinary branching in shower.
859 470 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
860 IF(MOD(MSTJ(43),2).EQ.1) THEN
861 P(N+1,4)=PEM*V(IM,1)
862 ELSE
863 P(N+1,4)=PEM*(0.5*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
864 & SQRT(PMLS)*ZM)/V(IM,5)
865 ENDIF
866 PHI=PARU(2)*RLU(0)
867 P(N+1,1)=PT*COS(PHI)
868 P(N+1,2)=PT*SIN(PHI)
869 IF(PZM.GT.0.) THEN
870 P(N+1,3)=0.5*(V(N+2,5)-V(N+1,5)-V(IM,5)+2.*PEM*P(N+1,4))/PZM
871 ELSE
872 P(N+1,3)=0.
873 ENDIF
874 P(N+2,1)=-P(N+1,1)
875 P(N+2,2)=-P(N+1,2)
876 P(N+2,3)=PZM-P(N+1,3)
877 P(N+2,4)=PEM-P(N+1,4)
878 IF(MSTJ(43).LE.2) THEN
879 V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
880 V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
881 ENDIF
882 ENDIF
883
884C...Rotate and boost daughters.
885 IF(IGM.GT.0) THEN
886 IF(MSTJ(43).LE.2) THEN
887 BEX=P(IGM,1)/P(IGM,4)
888 BEY=P(IGM,2)/P(IGM,4)
889 BEZ=P(IGM,3)/P(IGM,4)
890 GA=P(IGM,4)/P(IGM,5)
891 GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1.+GA)-
892 & P(IM,4))
893 ELSE
894 BEX=0.
895 BEY=0.
896 BEZ=0.
897 GA=1.
898 GABEP=0.
899 ENDIF
900 THE=ULANGL(P(IM,3)+GABEP*BEZ,SQRT((P(IM,1)+GABEP*BEX)**2+
901 & (P(IM,2)+GABEP*BEY)**2))
902 PHI=ULANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
903 DO 480 I=N+1,N+2
904 DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
905 & SIN(THE)*COS(PHI)*P(I,3)
906 DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
907 & SIN(THE)*SIN(PHI)*P(I,3)
908 DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
909 DP(4)=P(I,4)
910 DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
911 DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
912 P(I,1)=DP(1)+DGABP*BEX
913 P(I,2)=DP(2)+DGABP*BEY
914 P(I,3)=DP(3)+DGABP*BEZ
915 P(I,4)=GA*(DP(4)+DBP)
916 480 CONTINUE
917 ENDIF
918
919C...Weight with azimuthal distribution, if required.
920 IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
921 DO 490 J=1,3
922 DPT(1,J)=P(IM,J)
923 DPT(2,J)=P(IAU,J)
924 DPT(3,J)=P(N+1,J)
925 490 CONTINUE
926 DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
927 DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
928 DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
929 DO 500 J=1,3
930 DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/DPMM
931 DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/DPMM
932 500 CONTINUE
933 DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
934 DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
935 IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1*PARJ(82)) THEN
936 CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
937 & DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
938 IF(MAZIP.NE.0) THEN
939 IF(1.+HAZIP*(2.*CAD**2-1.).LT.RLU(0)*(1.+ABS(HAZIP)))
940 & GOTO 470
941 ENDIF
942 IF(MAZIC.NE.0) THEN
943 IF(MAZIC.EQ.N+2) CAD=-CAD
944 IF((1.-HAZIC)*(1.-HAZIC*CAD)/(1.+HAZIC**2-2.*HAZIC*CAD)
945 & .LT.RLU(0)) GOTO 470
946 ENDIF
947 ENDIF
948 ENDIF
949
950C...Azimuthal anisotropy due to interference with initial state partons.
951 IF(MOD(MIIS,2).EQ.1.AND.IGM.EQ.NS+1.AND.(K(N+1,2).EQ.21.OR.
952 &K(N+2,2).EQ.21)) THEN
953 III=IM-NS-1
954 IF(ISII(III).GE.1) THEN
955 IAZIID=N+1
956 IF(K(N+1,2).NE.21) IAZIID=N+2
957 IF(K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
958 & P(N+1,4).GT.P(N+2,4)) IAZIID=N+2
959 THEIID=ULANGL(P(IAZIID,3),SQRT(P(IAZIID,1)**2+P(IAZIID,2)**2))
960 IF(III.EQ.2) THEIID=PARU(1)-THEIID
961 PHIIID=ULANGL(P(IAZIID,1),P(IAZIID,2))
962 HAZII=MIN(0.95,THEIID/THEIIS(III,ISII(III)))
963 CAD=COS(PHIIID-PHIIIS(III,ISII(III)))
964 PHIREL=ABS(PHIIID-PHIIIS(III,ISII(III)))
965 IF(PHIREL.GT.PARU(1)) PHIREL=PARU(2)-PHIREL
966 IF((1.-HAZII)*(1.-HAZII*CAD)/(1.+HAZII**2-2.*HAZII*CAD)
967 & .LT.RLU(0)) GOTO 470
968 ENDIF
969 ENDIF
970
971C...Continue loop over partons that may branch, until none left.
972 IF(IGM.GE.0) K(IM,1)=14
973 N=N+NEP
974 NEP=2
975 IF(N.GT.MSTU(4)-MSTU(32)-5) THEN
976 CALL LUERRM(11,'(LUSHOW:) no more memory left in LUJETS')
977 IF(MSTU(21).GE.1) N=NS
978 IF(MSTU(21).GE.1) RETURN
979 ENDIF
980 GOTO 270
981
982C...Set information on imagined shower initiator.
983 510 IF(NPA.GE.2) THEN
984 K(NS+1,1)=11
985 K(NS+1,2)=94
986 K(NS+1,3)=IP1
987 IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
988 K(NS+1,4)=NS+2
989 K(NS+1,5)=NS+1+NPA
990 IIM=1
991 ELSE
992 IIM=0
993 ENDIF
994
995C...Reconstruct string drawing information.
996 DO 520 I=NS+1+IIM,N
997 IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
998 K(I,1)=1
999 ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND.
1000 &IABS(K(I,2)).LE.18) THEN
1001 K(I,1)=1
1002 ELSEIF(K(I,1).LE.10) THEN
1003 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
1004 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
1005 ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
1006 ID1=MOD(K(I,4),MSTU(5))
1007 IF(K(I,2).GE.1.AND.K(I,2).LE.8) ID1=MOD(K(I,4),MSTU(5))+1
1008 ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
1009 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
1010 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
1011 K(ID1,4)=K(ID1,4)+MSTU(5)*I
1012 K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
1013 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
1014 K(ID2,5)=K(ID2,5)+MSTU(5)*I
1015 ELSE
1016 ID1=MOD(K(I,4),MSTU(5))
1017 ID2=ID1+1
1018 K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
1019 K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
1020 IF(IABS(K(I,2)).LE.10.OR.K(ID1,1).GE.11) THEN
1021 K(ID1,4)=K(ID1,4)+MSTU(5)*I
1022 K(ID1,5)=K(ID1,5)+MSTU(5)*I
1023 ELSE
1024 K(ID1,4)=0
1025 K(ID1,5)=0
1026 ENDIF
1027 K(ID2,4)=0
1028 K(ID2,5)=0
1029 ENDIF
1030 520 CONTINUE
1031
1032C...Transformation from CM frame.
1033 IF(NPA.GE.2) THEN
1034 BEX=PS(1)/PS(4)
1035 BEY=PS(2)/PS(4)
1036 BEZ=PS(3)/PS(4)
1037 GA=PS(4)/PS(5)
1038 GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
1039 & /(1.+GA)-P(IPA(1),4))
1040 ELSE
1041 BEX=0.
1042 BEY=0.
1043 BEZ=0.
1044 GABEP=0.
1045 ENDIF
1046 THE=ULANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
1047 &+GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
1048 PHI=ULANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
1049 IF(NPA.EQ.3) THEN
1050 CHI=ULANGL(COS(THE)*COS(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(THE)*
1051 & SIN(PHI)*(P(IPA(2),2)+GABEP*BEY)-SIN(THE)*(P(IPA(2),3)+GABEP*
1052 & BEZ),-SIN(PHI)*(P(IPA(2),1)+GABEP*BEX)+COS(PHI)*(P(IPA(2),2)+
1053 & GABEP*BEY))
1054 MSTU(33)=1
1055 CALL LUDBRB(NS+1,N,0.,CHI,0D0,0D0,0D0)
1056 ENDIF
1057 DBEX=DBLE(BEX)
1058 DBEY=DBLE(BEY)
1059 DBEZ=DBLE(BEZ)
1060 MSTU(33)=1
1061 CALL LUDBRB(NS+1,N,THE,PHI,DBEX,DBEY,DBEZ)
1062
1063C...Decay vertex of shower.
1064 DO 540 I=NS+1,N
1065 DO 530 J=1,5
1066 V(I,J)=V(IP1,J)
1067 530 CONTINUE
1068 540 CONTINUE
1069
1070C...Delete trivial shower, else connect initiators.
1071 IF(N.EQ.NS+NPA+IIM) THEN
1072 N=NS
1073 ELSE
1074 DO 550 IP=1,NPA
1075 K(IPA(IP),1)=14
1076 K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
1077 K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
1078 K(NS+IIM+IP,3)=IPA(IP)
1079 IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
1080 IF(K(NS+IIM+IP,1).NE.1) THEN
1081 K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
1082 K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
1083 ENDIF
1084 550 CONTINUE
1085 ENDIF
1086
1087 RETURN
1088 END