]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | |
2 | C********************************************************************* | |
3 | ||
4 | SUBROUTINE LUSHOW(IP1,IP2,QMAX) | |
5 | ||
6 | C...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 | ||
17 | C...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 | ||
58 | C...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 | ||
79 | C...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) | |
88 | C...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 | ||
116 | C...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 | ||
135 | C...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 | ||
158 | C...Boost interfering initial partons to rest frame | |
159 | C...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 | ||
199 | C...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 | ||
220 | C...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 | ||
242 | C...Position of aunt (sister to branching parton). | |
243 | C...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 | ||
272 | C...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 | ||
288 | C...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 | ||
315 | C...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 | ||
345 | C...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 | ||
368 | C...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 | ||
382 | C...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 | ||
409 | C...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 | ||
415 | C...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 | ||
422 | C...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 | ||
429 | C...Reset QCD probability for lepton. | |
430 | IF(KFL(1).GE.11.AND.KFL(1).LE.18) FBR=0. | |
431 | ||
432 | C...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 | ||
438 | C...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 | ||
455 | C...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 | ||
473 | C...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 | ||
484 | C...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 | ||
492 | C...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 | ||
498 | C...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 | ||
519 | C...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 | ||
538 | C...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 | ||
577 | C...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 | ||
592 | C...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 | ||
620 | C...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 | ||
638 | C...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 | ||
652 | C...Impose angular constraint in first branching from interference | |
653 | C...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 | ||
663 | C...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 | ||
680 | C...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 | ||
766 | C...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 | ||
802 | C...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 | ||
817 | C...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 | ||
837 | C...Find coefficient of azimuthal asymmetry due to soft gluon | |
838 | C...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 | ||
858 | C...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 | ||
884 | C...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 | ||
919 | C...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 | ||
950 | C...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 | ||
971 | C...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 | ||
982 | C...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 | ||
995 | C...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 | ||
1032 | C...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 | ||
1063 | C...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 | ||
1070 | C...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 |