]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/pythia/pysspa.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Addition of CPV library as a separate detector from PHOS
[u/mrichter/AliRoot.git] / PYTHIA / pythia / pysspa.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE PYSSPA(IPU1,IPU2)
5
6C...Generates spacelike parton showers.
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 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
12 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
13 COMMON/PYINT1/MINT(400),VINT(400)
14 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
15 COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
16 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
17 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/
18 DIMENSION KFLS(4),IS(2),XS(2),ZS(2),Q2S(2),TEVCSV(2),TEVESV(2),
19 &XFS(2,-25:25),XFA(-25:25),XFB(-25:25),XFN(-25:25),WTAPC(-25:25),
20 &WTAPE(-25:25),WTSF(-25:25),THE2(2),ALAM(2),DQ2(3),DPC(3),DPD(4),
21 &DPB(4),ROBO(5),MORE(2),KFBEAM(2),Q2MNCS(2),KCFI(2),NFIS(2),
22 &THEFIS(2,2),ISFI(2)
23 DATA IS/2*0/
24
25C...Read out basic information; set global Q^2 scale.
26 IPUS1=IPU1
27 IPUS2=IPU2
28 ISUB=MINT(1)
29 Q2MX=VINT(56)
30 IF(ISET(ISUB).EQ.2) Q2MX=PARP(67)*VINT(56)
31
32C...Initialize QCD evolution and check phase space.
33 Q2MNC=PARP(62)**2
34 Q2MNCS(1)=Q2MNC
35 IF(MSTP(66).EQ.1.AND.MINT(107).EQ.3)
36 &Q2MNCS(1)=MAX(Q2MNC,VINT(283))
37 Q2MNCS(2)=Q2MNC
38 IF(MSTP(66).EQ.1.AND.MINT(108).EQ.3)
39 &Q2MNCS(2)=MAX(Q2MNC,VINT(284))
40 MCEV=0
41 XEC0=2.*PARP(65)/VINT(1)
42 ALAMS=PARU(112)
43 PARU(112)=PARP(61)
44 FQ2C=1.
45 TCMX=0.
46 IF(MINT(47).GE.2.AND.(MINT(47).NE.5.OR.MSTP(12).GE.1)) THEN
47 MCEV=1
48 IF(MSTP(64).EQ.1) FQ2C=PARP(63)
49 IF(MSTP(64).EQ.2) FQ2C=PARP(64)
50 TCMX=LOG(FQ2C*Q2MX/PARP(61)**2)
51 IF(Q2MX.LT.MAX(Q2MNC,2.*PARP(61)**2).OR.TCMX.LT.0.2)
52 & MCEV=0
53 ENDIF
54
55C...Initialize QED evolution and check phase space.
56 Q2MNE=PARP(68)**2
57 MEEV=0
58 XEE=1E-6
59 SPME=PMAS(11,1)**2
60 TEMX=0.
61 FWTE=10.
62 IF(MINT(45).EQ.3.OR.MINT(46).EQ.3) THEN
63 MEEV=1
64 TEMX=LOG(Q2MX/SPME)
65 IF(Q2MX.LE.Q2MNE.OR.TEMX.LT.0.2) MEEV=0
66 ENDIF
67 IF(MCEV.EQ.0.AND.MEEV.EQ.0) RETURN
68
69C...Initial values: flavours, momenta, virtualities.
70 NS=N
71 100 N=NS
72 DO 120 JT=1,2
73 MORE(JT)=1
74 KFBEAM(JT)=MINT(10+JT)
75 IF(MINT(18+JT).EQ.1)KFBEAM(JT)=22
76 KFLS(JT)=MINT(14+JT)
77 KFLS(JT+2)=KFLS(JT)
78 XS(JT)=VINT(40+JT)
79 IF(MINT(18+JT).EQ.1) XS(JT)=VINT(40+JT)/VINT(154+JT)
80 ZS(JT)=1.
81 Q2S(JT)=Q2MX
82 TEVCSV(JT)=TCMX
83 ALAM(JT)=PARP(61)
84 THE2(JT)=100.
85 TEVESV(JT)=TEMX
86 DO 110 KFL=-25,25
87 XFS(JT,KFL)=XSFX(JT,KFL)
88 110 CONTINUE
89 120 CONTINUE
90 DSH=VINT(44)
91 IF(ISET(ISUB).GE.3.AND.ISET(ISUB).LE.5) DSH=VINT(26)*VINT(2)
92
93C...Find if interference with final state partons.
94 MFIS=0
95 IF(MSTP(67).GE.1.AND.MSTP(67).LE.3) MFIS=MSTP(67)
96 IF(MFIS.NE.0) THEN
97 DO 140 I=1,2
98 KCFI(I)=0
99 KCA=LUCOMP(IABS(KFLS(I)))
100 IF(KCA.NE.0) KCFI(I)=KCHG(KCA,2)*ISIGN(1,KFLS(I))
101 NFIS(I)=0
102 IF(KCFI(I).NE.0) THEN
103 IF(I.EQ.1) IPFS=IPUS1
104 IF(I.EQ.2) IPFS=IPUS2
105 DO 130 J=1,2
106 ICSI=MOD(K(IPFS,3+J),MSTU(5))
107 IF(ICSI.GT.0.AND.ICSI.NE.IPUS1.AND.ICSI.NE.IPUS2.AND.
108 & (KCFI(I).EQ.(-1)**(J+1).OR.KCFI(I).EQ.2)) THEN
109 NFIS(I)=NFIS(I)+1
110 THEFIS(I,NFIS(I))=ULANGL(P(ICSI,3),SQRT(P(ICSI,1)**2+
111 & P(ICSI,2)**2))
112 IF(I.EQ.2) THEFIS(I,NFIS(I))=PARU(1)-THEFIS(I,NFIS(I))
113 ENDIF
114 130 CONTINUE
115 ENDIF
116 140 CONTINUE
117 IF(NFIS(1)+NFIS(2).EQ.0) MFIS=0
118 ENDIF
119
120C...Pick up leg with highest virtuality.
121 150 N=N+1
122 JT=1
123 IF(N.GT.NS+1.AND.Q2S(2).GT.Q2S(1)) JT=2
124 IF(MORE(JT).EQ.0) JT=3-JT
125 KFLB=KFLS(JT)
126 XB=XS(JT)
127 DO 160 KFL=-25,25
128 XFB(KFL)=XFS(JT,KFL)
129 160 CONTINUE
130 DSHR=2D0*SQRT(DSH)
131 DSHZ=DSH/DBLE(ZS(JT))
132
133C...Check if allowed to branch.
134 MCEV=0
135 IF(IABS(KFLB).LE.10.OR.KFLB.EQ.21) THEN
136 MCEV=1
137 XEC=MAX(XEC0,XB*(1./(1.-PARP(66))-1.))
138 IF(XB.GE.1.-2.*XEC) MCEV=0
139 ENDIF
140 MEEV=0
141 IF(MINT(44+JT).EQ.3) THEN
142 MEEV=1
143 IF(XB.GE.1.-2.*XEE) MEEV=0
144 IF((IABS(KFLB).LE.10.OR.KFLB.EQ.21).AND.XB.GE.1.-2.*XEC) MEEV=0
145C***Currently kill QED shower for resolved photoproduction.
146 IF(MINT(18+JT).EQ.1) MEEV=0
147C***Currently kill shower for W inside electron.
148 IF(IABS(KFLB).EQ.24) THEN
149 MCEV=0
150 MEEV=0
151 ENDIF
152 ENDIF
153 IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
154 Q2B=0.
155 GOTO 250
156 ENDIF
157
158C...Maximum Q2 with or without Q2 ordering. Effective Lambda and n_f.
159 Q2B=Q2S(JT)
160 TEVCB=TEVCSV(JT)
161 TEVEB=TEVESV(JT)
162 IF(MSTP(62).LE.1) THEN
163 IF(ZS(JT).GT.0.99999) THEN
164 Q2B=Q2S(JT)
165 ELSE
166 Q2B=0.5*(1./ZS(JT)+1.)*Q2S(JT)+0.5*(1./ZS(JT)-1.)*(Q2S(3-JT)-
167 & SNGL(DSH)+SQRT((SNGL(DSH)+Q2S(1)+Q2S(2))**2+8.*Q2S(1)*Q2S(2)*
168 & ZS(JT)/(1.-ZS(JT))))
169 ENDIF
170 IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
171 IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
172 ENDIF
173 IF(MCEV.EQ.1) THEN
174 ALSDUM=ULALPS(FQ2C*Q2B)
175 TEVCB=TEVCB+2.*LOG(ALAM(JT)/PARU(117))
176 ALAM(JT)=PARU(117)
177 B0=(33.-2.*MSTU(118))/6.
178 ENDIF
179 TEVCBS=TEVCB
180 TEVEBS=TEVEB
181
182C...Select side for interference with final state partons.
183 IF(MFIS.GE.1.AND.N.LE.NS+2) THEN
184 IFI=N-NS
185 ISFI(IFI)=0
186 IF(IABS(KCFI(IFI)).EQ.1.AND.NFIS(IFI).EQ.1) THEN
187 ISFI(IFI)=1
188 ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.1) THEN
189 IF(RLU(0).GT.0.5) ISFI(IFI)=1
190 ELSEIF(KCFI(IFI).EQ.2.AND.NFIS(IFI).EQ.2) THEN
191 ISFI(IFI)=1
192 IF(RLU(0).GT.0.5) ISFI(IFI)=2
193 ENDIF
194 ENDIF
195
196C...Calculate Altarelli-Parisi weights.
197 DO 170 KFL=-25,25
198 WTAPC(KFL)=0.
199 WTAPE(KFL)=0.
200 WTSF(KFL)=0.
201 170 CONTINUE
202C...q -> q, g -> q.
203 IF(IABS(KFLB).LE.10) THEN
204 WTAPC(KFLB)=(8./3.)*LOG((1.-XEC-XB)*(XB+XEC)/(XEC*(1.-XEC)))
205 WTAPC(21)=0.5*(XB/(XB+XEC)-XB/(1.-XEC))
206C...f -> f, gamma -> f.
207 ELSEIF(IABS(KFLB).LE.20) THEN
208 WTAPF1=LOG((1.-XEE-XB)*(XB+XEE)/(XEE*(1.-XEE)))
209 WTAPF2=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))
210 WTAPE(KFLB)=2.*(WTAPF1+WTAPF2)
211 IF(MSTP(12).GE.1) WTAPE(22)=XB/(XB+XEE)-XB/(1.-XEE)
212C...f -> g, g -> g.
213 ELSEIF(KFLB.EQ.21) THEN
214 WTAPQ=(16./3.)*(SQRT((1.-XEC)/XB)-SQRT((XB+XEC)/XB))
215 DO 180 KFL=1,MSTP(58)
216 WTAPC(KFL)=WTAPQ
217 WTAPC(-KFL)=WTAPQ
218 180 CONTINUE
219 WTAPC(21)=6.*LOG((1.-XEC-XB)/XEC)
220C...f -> gamma, W+, W-.
221 ELSEIF(KFLB.EQ.22) THEN
222 WTAPF=LOG((1.-XEE-XB)*(1.-XEE)/(XEE*(XB+XEE)))/XB
223 WTAPE(11)=WTAPF
224 WTAPE(-11)=WTAPF
225 ELSEIF(KFLB.EQ.24) THEN
226 WTAPE(-11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
227 & (XEE*(XB+XEE)))/XB
228 ELSEIF(KFLB.EQ.-24) THEN
229 WTAPE(11)=1./(4.*PARU(102))*LOG((1.-XEE-XB)*(1.-XEE)/
230 & (XEE*(XB+XEE)))/XB
231 ENDIF
232
233C...Calculate structure function weights and sum.
234 NTRY=0
235 190 NTRY=NTRY+1
236 IF(NTRY.GT.500) THEN
237 MINT(51)=1
238 RETURN
239 ENDIF
240 WTSUMC=0.
241 WTSUME=0.
242 XFBO=MAX(1E-10,XFB(KFLB))
243 DO 200 KFL=-25,25
244 WTSF(KFL)=XFB(KFL)/XFBO
245 WTSUMC=WTSUMC+WTAPC(KFL)*WTSF(KFL)
246 WTSUME=WTSUME+WTAPE(KFL)*WTSF(KFL)
247 200 CONTINUE
248 WTSUMC=MAX(0.0001,WTSUMC)
249 WTSUME=MAX(0.0001/FWTE,WTSUME)
250
251C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
252 NTRY2=0
253 210 NTRY2=NTRY2+1
254 IF(NTRY2.GT.500) THEN
255 MINT(51)=1
256 RETURN
257 ENDIF
258 IF(MCEV.EQ.1) THEN
259 IF(MSTP(64).LE.0) THEN
260 TEVCB=TEVCB+LOG(RLU(0))*PARU(2)/(PARU(111)*WTSUMC)
261 ELSEIF(MSTP(64).EQ.1) THEN
262 TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/WTSUMC))
263 ELSE
264 TEVCB=TEVCB*EXP(MAX(-50.,LOG(RLU(0))*B0/(5.*WTSUMC)))
265 ENDIF
266 ENDIF
267 IF(MEEV.EQ.1) THEN
268 TEVEB=TEVEB*EXP(MAX(-50.,LOG(RLU(0))*PARU(2)/
269 & (PARU(101)*FWTE*WTSUME*TEMX)))
270 ENDIF
271
272C...Translate t into Q2 scale; choose between QCD and QED evolution.
273 220 IF(MCEV.EQ.1) Q2CB=ALAM(JT)**2*EXP(MAX(-50.,TEVCB))/FQ2C
274 IF(MEEV.EQ.1) Q2EB=SPME*EXP(MAX(-50.,TEVEB))
275 MCE=0
276 IF(MCEV.EQ.0.AND.MEEV.EQ.0) THEN
277 ELSEIF(MCEV.EQ.1.AND.MEEV.EQ.0) THEN
278 IF(Q2CB.GT.Q2MNCS(JT)) MCE=1
279 ELSEIF(MCEV.EQ.0.AND.MEEV.EQ.1) THEN
280 IF(Q2EB.GT.Q2MNE) MCE=2
281 ELSEIF(Q2MNCS(JT).GT.Q2MNE) THEN
282 MCE=1
283 IF(Q2EB.GT.Q2CB.OR.Q2CB.LE.Q2MNCS(JT)) MCE=2
284 IF(MCE.EQ.2.AND.Q2EB.LE.Q2MNE) MCE=0
285 ELSE
286 MCE=2
287 IF(Q2CB.GT.Q2EB.OR.Q2EB.LE.Q2MNE) MCE=1
288 IF(MCE.EQ.1.AND.Q2CB.LE.Q2MNCS(JT)) MCE=0
289 ENDIF
290
291C...Evolution possibly ended. Update t values.
292 IF(MCE.EQ.0) THEN
293 Q2B=0.
294 GOTO 250
295 ELSEIF(MCE.EQ.1) THEN
296 Q2B=Q2CB
297 Q2REF=FQ2C*Q2B
298 IF(MEEV.EQ.1) TEVEB=LOG(Q2B/SPME)
299 ELSE
300 Q2B=Q2EB
301 Q2REF=Q2B
302 IF(MCEV.EQ.1) TEVCB=LOG(FQ2C*Q2B/ALAM(JT)**2)
303 ENDIF
304
305C...Select flavour for branching parton.
306 IF(MCE.EQ.1) WTRAN=RLU(0)*WTSUMC
307 IF(MCE.EQ.2) WTRAN=RLU(0)*WTSUME
308 KFLA=-25
309 230 KFLA=KFLA+1
310 IF(MCE.EQ.1) WTRAN=WTRAN-WTAPC(KFLA)*WTSF(KFLA)
311 IF(MCE.EQ.2) WTRAN=WTRAN-WTAPE(KFLA)*WTSF(KFLA)
312 IF(KFLA.LE.24.AND.WTRAN.GT.0.) GOTO 230
313 IF(KFLA.EQ.25) THEN
314 Q2B=0.
315 GOTO 250
316 ENDIF
317
318C...Choose z value and corrective weight.
319 WTZ=0.
320C...q -> q + g.
321 IF(IABS(KFLA).LE.10.AND.IABS(KFLB).LE.10) THEN
322 Z=1.-((1.-XB-XEC)/(1.-XEC))*
323 & (XEC*(1.-XEC)/((XB+XEC)*(1.-XB-XEC)))**RLU(0)
324 WTZ=0.5*(1.+Z**2)
325C...q -> g + q.
326 ELSEIF(IABS(KFLA).LE.10.AND.KFLB.EQ.21) THEN
327 Z=XB/(SQRT(XB+XEC)+RLU(0)*(SQRT(1.-XEC)-SQRT(XB+XEC)))**2
328 WTZ=0.5*(1.+(1.-Z)**2)*SQRT(Z)
329C...f -> f + gamma.
330 ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).LE.20) THEN
331 IF(WTAPF1.GT.RLU(0)*(WTAPF1+WTAPF2)) THEN
332 Z=1.-((1.-XB-XEE)/(1.-XEE))*
333 & (XEE*(1.-XEE)/((XB+XEE)*(1.-XB-XEE)))**RLU(0)
334 ELSE
335 Z=XB+XB*(XEE/(1.-XEE))*
336 & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
337 ENDIF
338 WTZ=0.5*(1.+Z**2)*(Z-XB)/(1.-XB)
339C...f -> gamma + f.
340 ELSEIF(IABS(KFLA).LE.20.AND.KFLB.EQ.22) THEN
341 Z=XB+XB*(XEE/(1.-XEE))*
342 & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
343 WTZ=0.5*(1.+(1.-Z)**2)*XB*(Z-XB)/Z
344C...f -> W+- + f'.
345 ELSEIF(IABS(KFLA).LE.20.AND.IABS(KFLB).EQ.24) THEN
346 Z=XB+XB*(XEE/(1.-XEE))*
347 & ((1.-XB-XEE)*(1.-XEE)/(XEE*(XB+XEE)))**RLU(0)
348 WTZ=0.5*(1.+(1.-Z)**2)*(XB*(Z-XB)/Z)*(Q2B/(Q2B+PMAS(24,1)**2))
349C...g -> q + q~.
350 ELSEIF(KFLA.EQ.21.AND.IABS(KFLB).LE.10) THEN
351 Z=XB/(1.-XEC)+RLU(0)*(XB/(XB+XEC)-XB/(1.-XEC))
352 WTZ=1.-2.*Z*(1.-Z)
353C...g -> g + g.
354 ELSEIF(KFLA.EQ.21.AND.KFLB.EQ.21) THEN
355 Z=1./(1.+((1.-XEC-XB)/XB)*(XEC/(1.-XEC-XB))**RLU(0))
356 WTZ=(1.-Z*(1.-Z))**2
357C...gamma -> f + f~.
358 ELSEIF(KFLA.EQ.22.AND.IABS(KFLB).LE.20) THEN
359 Z=XB/(1.-XEE)+RLU(0)*(XB/(XB+XEE)-XB/(1.-XEE))
360 WTZ=1.-2.*Z*(1.-Z)
361 ENDIF
362 IF(MCE.EQ.2) WTZ=(WTZ/FWTE)*(TEVEB/TEMX)
363
364C...Option with resummation of soft gluon emission as effective z shift.
365 IF(MCE.EQ.1) THEN
366 IF(MSTP(65).GE.1) THEN
367 RSOFT=6.
368 IF(KFLB.NE.21) RSOFT=8./3.
369 Z=Z*(TEVCB/TEVCSV(JT))**(RSOFT*XEC/((XB+XEC)*B0))
370 IF(Z.LE.XB) GOTO 210
371 ENDIF
372
373C...Option with alpha_s(k_T^2): demand k_T^2 > cutoff, reweight.
374 IF(MSTP(64).GE.2) THEN
375 IF((1.-Z)*Q2B.LT.Q2MNCS(JT)) GOTO 210
376 ALPRAT=TEVCB/(TEVCB+LOG(1.-Z))
377 IF(ALPRAT.LT.5.*RLU(0)) GOTO 210
378 IF(ALPRAT.GT.5.) WTZ=WTZ*ALPRAT/5.
379 ENDIF
380
381C...Impose angular constraint in first branching from interference
382C...with final state partons.
383 IF(MFIS.GE.1.AND.N.LE.NS+2.AND.NTRY2.LT.200) THEN
384 THE2D=(4.*Q2B)/(DSH*(1.-Z))
385 IF(N.EQ.NS+1.AND.ISFI(1).GE.1) THEN
386 IF(THE2D.GT.THEFIS(1,ISFI(1))**2) GOTO 210
387 ELSEIF(N.EQ.NS+2.AND.ISFI(2).GE.1) THEN
388 IF(THE2D.GT.THEFIS(2,ISFI(2))**2) GOTO 210
389 ENDIF
390 ENDIF
391
392C...Option with angular ordering requirement.
393 IF(MSTP(62).GE.3.AND.NTRY2.LT.200) THEN
394 THE2T=(4.*Z**2*Q2B)/(VINT(2)*(1.-Z)*XB**2)
395 IF(THE2T.GT.THE2(JT)) GOTO 210
396 ENDIF
397 ENDIF
398
399C...Weighting with new structure functions.
400 MINT(105)=MINT(102+JT)
401 MINT(109)=MINT(106+JT)
402 IF(MSTP(57).LE.1) THEN
403 CALL PYSTFU(KFBEAM(JT),XB,Q2REF,XFN)
404 ELSE
405 CALL PYSTFL(KFBEAM(JT),XB,Q2REF,XFN)
406 ENDIF
407 XFBN=XFN(KFLB)
408 IF(XFBN.LT.1E-20) THEN
409 IF(KFLA.EQ.KFLB) THEN
410 TEVCB=TEVCBS
411 TEVEB=TEVEBS
412 WTAPC(KFLB)=0.
413 WTAPE(KFLB)=0.
414 GOTO 190
415 ELSEIF(MCE.EQ.1.AND.TEVCBS-TEVCB.GT.0.2) THEN
416 TEVCB=0.5*(TEVCBS+TEVCB)
417 GOTO 220
418 ELSEIF(MCE.EQ.2.AND.TEVEBS-TEVEB.GT.0.2) THEN
419 TEVEB=0.5*(TEVEBS+TEVEB)
420 GOTO 220
421 ELSE
422 XFBN=1E-10
423 XFN(KFLB)=XFBN
424 ENDIF
425 ENDIF
426 DO 240 KFL=-25,25
427 XFB(KFL)=XFN(KFL)
428 240 CONTINUE
429 XA=XB/Z
430 IF(MSTP(57).LE.1) THEN
431 CALL PYSTFU(KFBEAM(JT),XA,Q2REF,XFA)
432 ELSE
433 CALL PYSTFL(KFBEAM(JT),XA,Q2REF,XFA)
434 ENDIF
435 XFAN=XFA(KFLA)
436 IF(XFAN.LT.1E-20) GOTO 190
437 WTSFA=WTSF(KFLA)
438 IF(WTZ*XFAN/XFBN.LT.RLU(0)*WTSFA) GOTO 190
439
440C...Define two hard scatterers in their CM-frame.
441 250 IF(N.EQ.NS+2) THEN
442 DQ2(JT)=Q2B
443 DPLCM=SQRT((DSH+DQ2(1)+DQ2(2))**2-4D0*DQ2(1)*DQ2(2))/DSHR
444 DO 270 JR=1,2
445 I=NS+JR
446 IF(JR.EQ.1) IPO=IPUS1
447 IF(JR.EQ.2) IPO=IPUS2
448 DO 260 J=1,5
449 K(I,J)=0
450 P(I,J)=0.
451 V(I,J)=0.
452 260 CONTINUE
453 K(I,1)=14
454 K(I,2)=KFLS(JR+2)
455 K(I,4)=IPO
456 K(I,5)=IPO
457 P(I,3)=DPLCM*(-1)**(JR+1)
458 P(I,4)=(DSH+DQ2(3-JR)-DQ2(JR))/DSHR
459 P(I,5)=-SQRT(SNGL(DQ2(JR)))
460 K(IPO,1)=14
461 K(IPO,3)=I
462 K(IPO,4)=MOD(K(IPO,4),MSTU(5))+MSTU(5)*I
463 K(IPO,5)=MOD(K(IPO,5),MSTU(5))+MSTU(5)*I
464 270 CONTINUE
465
466C...Find maximum allowed mass of timelike parton.
467 ELSEIF(N.GT.NS+2) THEN
468 JR=3-JT
469 DQ2(3)=Q2B
470 DPC(1)=P(IS(1),4)
471 DPC(2)=P(IS(2),4)
472 DPC(3)=0.5*(ABS(P(IS(1),3))+ABS(P(IS(2),3)))
473 DPD(1)=DSH+DQ2(JR)+DQ2(JT)
474 DPD(2)=DSHZ+DQ2(JR)+DQ2(3)
475 DPD(3)=SQRT(DPD(1)**2-4D0*DQ2(JR)*DQ2(JT))
476 DPD(4)=SQRT(DPD(2)**2-4D0*DQ2(JR)*DQ2(3))
477 IKIN=0
478 IF(Q2S(JR).GE.0.25*Q2MNC.AND.DPD(1)-DPD(3).GE.
479 & 1D-10*DPD(1)) IKIN=1
480 IF(IKIN.EQ.0) DMSMA=(DQ2(JT)/DBLE(ZS(JT))-DQ2(3))*(DSH/
481 & (DSH+DQ2(JT))-DSH/(DSHZ+DQ2(3)))
482 IF(IKIN.EQ.1) DMSMA=(DPD(1)*DPD(2)-DPD(3)*DPD(4))/(2.*
483 & DQ2(JR))-DQ2(JT)-DQ2(3)
484
485C...Generate timelike parton shower (if required).
486 IT=N
487 DO 280 J=1,5
488 K(IT,J)=0
489 P(IT,J)=0.
490 V(IT,J)=0.
491 280 CONTINUE
492 K(IT,1)=3
493C...f -> f + g (gamma).
494 IF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).LE.20) THEN
495 K(IT,2)=21
496 IF(IABS(KFLB).GE.11) K(IT,2)=22
497C...f -> g (gamma, W+-) + f.
498 ELSEIF(IABS(KFLB).LE.20.AND.IABS(KFLS(JT+2)).GT.20) THEN
499 K(IT,2)=KFLB
500 IF(KFLS(JT+2).EQ.24) THEN
501 K(IT,2)=-12
502 ELSEIF(KFLS(JT+2).EQ.-24) THEN
503 K(IT,2)=12
504 ENDIF
505C...g (gamma) -> f + f~, g + g.
506 ELSE
507 K(IT,2)=-KFLS(JT+2)
508 IF(KFLS(JT+2).GT.20) K(IT,2)=KFLS(JT+2)
509 ENDIF
510 P(IT,5)=ULMASS(K(IT,2))
511 IF(SNGL(DMSMA).LE.P(IT,5)**2) GOTO 100
512 IF(MSTP(63).GE.1.AND.MCE.EQ.1) THEN
513 MSTJ48=MSTJ(48)
514 PARJ85=PARJ(85)
515 P(IT,4)=(DSHZ-DSH-P(IT,5)**2)/DSHR
516 P(IT,3)=SQRT(P(IT,4)**2-P(IT,5)**2)
517 IF(MSTP(63).EQ.1) THEN
518 Q2TIM=DMSMA
519 ELSEIF(MSTP(63).EQ.2) THEN
520 Q2TIM=MIN(SNGL(DMSMA),PARP(71)*Q2S(JT))
521 ELSE
522 Q2TIM=DMSMA
523 MSTJ(48)=1
524 IF(IKIN.EQ.0) DPT2=DMSMA*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
525 IF(IKIN.EQ.1) DPT2=DMSMA*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
526 & DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)))/(4.*DSH*DPC(3)**2)
527 PARJ(85)=SQRT(MAX(0.,SNGL(DPT2)))*
528 & (1./P(IT,4)+1./P(IS(JT),4))
529 ENDIF
530 CALL LUSHOW(IT,0,SQRT(Q2TIM))
531 MSTJ(48)=MSTJ48
532 PARJ(85)=PARJ85
533 IF(N.GE.IT+1) P(IT,5)=P(IT+1,5)
534 ENDIF
535
536C...Reconstruct kinematics of branching: timelike parton shower.
537 DMS=P(IT,5)**2
538 IF(IKIN.EQ.0) DPT2=(DMSMA-DMS)*(DSHZ+DQ2(3))/(DSH+DQ2(JT))
539 IF(IKIN.EQ.1) DPT2=(DMSMA-DMS)*(0.5*DPD(1)*DPD(2)+0.5*DPD(3)*
540 & DPD(4)-DQ2(JR)*(DQ2(JT)+DQ2(3)+DMS))/(4.*DSH*DPC(3)**2)
541 IF(DPT2.LT.0.) GOTO 100
542 DPB(1)=(0.5*DPD(2)-DPC(JR)*(DSHZ+DQ2(JR)-DQ2(JT)-DMS)/
543 & DSHR)/DPC(3)-DPC(3)
544 P(IT,1)=SQRT(SNGL(DPT2))
545 P(IT,3)=DPB(1)*(-1)**(JT+1)
546 P(IT,4)=SQRT(DPT2+DPB(1)**2+DMS)
547 IF(N.GE.IT+1) THEN
548 DPB(1)=SQRT(DPB(1)**2+DPT2)
549 DPB(2)=SQRT(DPB(1)**2+DMS)
550 DPB(3)=P(IT+1,3)
551 DPB(4)=SQRT(DPB(3)**2+DMS)
552 DBEZ=(DPB(4)*DPB(1)-DPB(3)*DPB(2))/(DPB(4)*DPB(2)-DPB(3)*
553 & DPB(1))
554 CALL LUDBRB(IT+1,N,0.,0.,0D0,0D0,DBEZ)
555 THE=ULANGL(P(IT,3),P(IT,1))
556 CALL LUDBRB(IT+1,N,THE,0.,0D0,0D0,0D0)
557 ENDIF
558
559C...Reconstruct kinematics of branching: spacelike parton.
560 DO 290 J=1,5
561 K(N+1,J)=0
562 P(N+1,J)=0.
563 V(N+1,J)=0.
564 290 CONTINUE
565 K(N+1,1)=14
566 K(N+1,2)=KFLB
567 P(N+1,1)=P(IT,1)
568 P(N+1,3)=P(IT,3)+P(IS(JT),3)
569 P(N+1,4)=P(IT,4)+P(IS(JT),4)
570 P(N+1,5)=-SQRT(SNGL(DQ2(3)))
571
572C...Define colour flow of branching.
573 K(IS(JT),3)=N+1
574 K(IT,3)=N+1
575 IM1=N+1
576 IM2=N+1
577C...f -> f + gamma (Z, W).
578 IF(IABS(K(IT,2)).GE.22) THEN
579 K(IT,1)=1
580 ID1=IS(JT)
581 ID2=IS(JT)
582C...f -> gamma (Z, W) + f.
583 ELSEIF(IABS(K(IS(JT),2)).GE.22) THEN
584 ID1=IT
585 ID2=IT
586C...gamma -> q + q~, g + g.
587 ELSEIF(K(N+1,2).EQ.22) THEN
588 ID1=IS(JT)
589 ID2=IT
590 IM1=ID2
591 IM2=ID1
592C...q -> q + g.
593 ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21.AND.K(IT,2).EQ.21) THEN
594 ID1=IT
595 ID2=IS(JT)
596C...q -> g + q.
597 ELSEIF(K(N+1,2).GT.0.AND.K(N+1,2).NE.21) THEN
598 ID1=IS(JT)
599 ID2=IT
600C...q~ -> q~ + g.
601 ELSEIF(K(N+1,2).LT.0.AND.K(IT,2).EQ.21) THEN
602 ID1=IS(JT)
603 ID2=IT
604C...q~ -> g + q~.
605 ELSEIF(K(N+1,2).LT.0) THEN
606 ID1=IT
607 ID2=IS(JT)
608C...g -> g + g; g -> q + q~.
609 ELSEIF((K(IT,2).EQ.21.AND.RLU(0).GT.0.5).OR.K(IT,2).LT.0) THEN
610 ID1=IS(JT)
611 ID2=IT
612 ELSE
613 ID1=IT
614 ID2=IS(JT)
615 ENDIF
616 IF(IM1.EQ.N+1) K(IM1,4)=K(IM1,4)+ID1
617 IF(IM2.EQ.N+1) K(IM2,5)=K(IM2,5)+ID2
618 K(ID1,4)=K(ID1,4)+MSTU(5)*IM1
619 K(ID2,5)=K(ID2,5)+MSTU(5)*IM2
620 IF(ID1.NE.ID2) THEN
621 K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
622 K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
623 ENDIF
624 N=N+1
625
626C...Boost to new CM-frame.
627 DBSVX=DBLE((P(N,1)+P(IS(JR),1))/(P(N,4)+P(IS(JR),4)))
628 DBSVZ=DBLE((P(N,3)+P(IS(JR),3))/(P(N,4)+P(IS(JR),4)))
629 IF(DBSVX**2+DBSVZ**2.GE.1D0) GOTO 100
630 CALL LUDBRB(NS+1,N,0.,0.,-DBSVX,0D0,-DBSVZ)
631 IR=N+(JT-1)*(IS(1)-N)
632 CALL LUDBRB(NS+1,N,-ULANGL(P(IR,3),P(IR,1)),PARU(2)*RLU(0),
633 & 0D0,0D0,0D0)
634 ENDIF
635
636C...Update kinematics variables.
637 IS(JT)=N
638 DQ2(JT)=Q2B
639 IF(MSTP(62).GE.3) THE2(JT)=THE2T
640 DSH=DSHZ
641
642C...Save quantities; loop back.
643 Q2S(JT)=Q2B
644 IF((MCEV.EQ.1.AND.Q2B.GE.0.25*Q2MNC).OR.
645 &(MEEV.EQ.1.AND.Q2B.GE.Q2MNE)) THEN
646 KFLS(JT+2)=KFLS(JT)
647 KFLS(JT)=KFLA
648 XS(JT)=XA
649 ZS(JT)=Z
650 DO 300 KFL=-25,25
651 XFS(JT,KFL)=XFA(KFL)
652 300 CONTINUE
653 TEVCSV(JT)=TEVCB
654 TEVESV(JT)=TEVEB
655 ELSE
656 MORE(JT)=0
657 IF(JT.EQ.1) IPU1=N
658 IF(JT.EQ.2) IPU2=N
659 ENDIF
660 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
661 CALL LUERRM(11,'(PYSSPA:) no more memory left in LUJETS')
662 IF(MSTU(21).GE.1) N=NS
663 IF(MSTU(21).GE.1) RETURN
664 ENDIF
665 IF(MORE(1).EQ.1.OR.MORE(2).EQ.1) GOTO 150
666
667C...Boost hard scattering partons to frame of shower initiators.
668 DO 310 J=1,3
669 ROBO(J+2)=(P(NS+1,J)+P(NS+2,J))/(P(NS+1,4)+P(NS+2,4))
670 310 CONTINUE
671 K(N+2,1)=1
672 DO 320 J=1,5
673 P(N+2,J)=P(NS+1,J)
674 320 CONTINUE
675 ROBOT=ROBO(3)**2+ROBO(4)**2+ROBO(5)**2
676 IF(ROBOT.GE.0.999999) THEN
677 ROBOT=1.00001*SQRT(ROBOT)
678 ROBO(3)=ROBO(3)/ROBOT
679 ROBO(4)=ROBO(4)/ROBOT
680 ROBO(5)=ROBO(5)/ROBOT
681 ENDIF
682 CALL LUDBRB(N+2,N+2,0.,0.,-DBLE(ROBO(3)),-DBLE(ROBO(4)),
683 &-DBLE(ROBO(5)))
684 ROBO(2)=ULANGL(P(N+2,1),P(N+2,2))
685 ROBO(1)=ULANGL(P(N+2,3),SQRT(P(N+2,1)**2+P(N+2,2)**2))
686 CALL LUDBRB(MINT(83)+5,NS,ROBO(1),ROBO(2),DBLE(ROBO(3)),
687 &DBLE(ROBO(4)),DBLE(ROBO(5)))
688
689C...Store user information. Reset Lambda value.
690 K(IPU1,3)=MINT(83)+3
691 K(IPU2,3)=MINT(83)+4
692 DO 330 JT=1,2
693 MINT(12+JT)=KFLS(JT)
694 VINT(140+JT)=XS(JT)
695 IF(MINT(18+JT).EQ.1) VINT(140+JT)=VINT(154+JT)*XS(JT)
696 330 CONTINUE
697 PARU(112)=ALAMS
698
699 RETURN
700 END