]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/pythia/pyrand.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
added copyright, updated class description
[u/mrichter/AliRoot.git] / PYTHIA / pythia / pyrand.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE PYRAND
5
6C...Generates quantities characterizing the high-pT scattering at the
7C...parton level according to the matrix elements. Chooses incoming,
8C...reacting partons, their momentum fractions and one of the possible
9C...subprocesses.
10 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
11 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
12 COMMON/PYSUBS/MSEL,MSUB(200),KFIN(2,-40:40),CKIN(200)
13 COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
14 COMMON/PYINT1/MINT(400),VINT(400)
15 COMMON/PYINT2/ISET(200),KFPR(200,2),COEF(200,20),ICOL(40,4,2)
16 COMMON/PYINT3/XSFX(2,-40:40),ISIG(1000,3),SIGH(1000)
17 COMMON/PYINT4/WIDP(21:40,0:40),WIDE(21:40,0:40),WIDS(21:40,3)
18 COMMON/PYINT5/NGEN(0:200,3),XSEC(0:200,3)
19 COMMON/PYINT7/SIGT(0:6,0:6,0:5)
20 COMMON/PYINT9/DXSEC(0:200)
21 DOUBLE PRECISION DXSEC
22 COMMON/PYUPPR/NUP,KUP(20,7),PUP(20,5),NFUP,IFUP(10,2),Q2UP(0:10)
23 SAVE /LUDAT1/,/LUDAT2/
24 SAVE /PYSUBS/,/PYPARS/,/PYINT1/,/PYINT2/,/PYINT3/,/PYINT4/,
25 &/PYINT5/,/PYINT7/,/PYINT9/,/PYUPPR/
26 DIMENSION XPQ(-25:25),PMM(2),PDIF(4),BHAD(4)
27
28C...Parameters and data used in elastic/diffractive treatment.
29 DATA EPS/0.0808/, ALP/0.25/, CRES/2./, PMRC/1.062/, SMP/0.880/
30 DATA BHAD/2.3,1.4,1.4,0.23/
31
32C...Initial values, specifically for (first) semihard interaction.
33 MINT(10)=0
34 MINT(17)=0
35 MINT(18)=0
36 VINT(143)=1.
37 VINT(144)=1.
38 MFAIL=0
39 IF(MSTP(171).EQ.1.AND.MSTP(172).EQ.2) MFAIL=1
40 ISUB=0
41 LOOP=0
42 100 LOOP=LOOP+1
43 MINT(51)=0
44
45C...Choice of process type - first event of pileup.
46 IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN
47
48C...For gamma-p or gamma-gamma first pick between alternatives.
49 IF(MINT(121).GT.1) CALL PYSAVE(4,IGA)
50 MINT(122)=IGA
51
52C...For gamma + gamma with different nature, flip at random.
53 IF(MINT(11).EQ.22.AND.MINT(12).EQ.22.AND.MINT(123).GE.4.AND.
54 & RLU(0).GT.0.5) THEN
55 MINTSV=MINT(41)
56 MINT(41)=MINT(42)
57 MINT(42)=MINTSV
58 MINTSV=MINT(45)
59 MINT(45)=MINT(46)
60 MINT(46)=MINTSV
61 MINTSV=MINT(107)
62 MINT(107)=MINT(108)
63 MINT(108)=MINTSV
64 IF(MINT(47).EQ.2.OR.MINT(47).EQ.3) MINT(47)=5-MINT(47)
65 ENDIF
66
67C...Pick process type.
68 RSUB=XSEC(0,1)*RLU(0)
69 DO 110 I=1,200
70 IF(MSUB(I).NE.1) GOTO 110
71 ISUB=I
72 RSUB=RSUB-XSEC(I,1)
73 IF(RSUB.LE.0.) GOTO 120
74 110 CONTINUE
75 120 IF(ISUB.EQ.95) ISUB=96
76 IF(ISUB.EQ.96) CALL PYMULT(2)
77
78C...Choice of inclusive process type - pileup events.
79 ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN
80 RSUB=VINT(131)*RLU(0)
81 ISUB=96
82 IF(RSUB.GT.SIGT(0,0,5)) ISUB=94
83 IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)) ISUB=93
84 IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)) ISUB=92
85 IF(RSUB.GT.SIGT(0,0,5)+SIGT(0,0,4)+SIGT(0,0,3)+SIGT(0,0,2))
86 & ISUB=91
87 IF(ISUB.EQ.96) CALL PYMULT(2)
88 ENDIF
89 IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1
90 IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1
91 IF(ISUB.EQ.96.AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
92 &NGEN(97,1)=NGEN(97,1)+1
93 MINT(1)=ISUB
94 ISTSB=ISET(ISUB)
95
96C...Find resonances (explicit or implicit in cross-section).
97 MINT(72)=0
98 KFR1=0
99 IF(ISTSB.EQ.1.OR.ISTSB.EQ.3.OR.ISTSB.EQ.5) THEN
100 KFR1=KFPR(ISUB,1)
101 ELSEIF(ISUB.EQ.24.OR.ISUB.EQ.25.OR.ISUB.EQ.110.OR.ISUB.EQ.165.OR.
102 &ISUB.EQ.171.OR.ISUB.EQ.176) THEN
103 KFR1=23
104 ELSEIF(ISUB.EQ.23.OR.ISUB.EQ.26.OR.ISUB.EQ.166.OR.ISUB.EQ.172.OR.
105 &ISUB.EQ.177) THEN
106 KFR1=24
107 ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
108 KFR1=25
109 IF(MSTP(46).EQ.5) THEN
110 KFR1=30
111 PMAS(30,1)=PARP(45)
112 PMAS(30,2)=PARP(45)**3/(96.*PARU(1)*PARP(47)**2)
113 ENDIF
114 ENDIF
115 CKMX=CKIN(2)
116 IF(CKMX.LE.0.) CKMX=VINT(1)
117 IF(KFR1.NE.0) THEN
118 IF(CKIN(1).GT.PMAS(KFR1,1)+20.*PMAS(KFR1,2).OR.
119 & CKMX.LT.PMAS(KFR1,1)-20.*PMAS(KFR1,2)) KFR1=0
120 ENDIF
121 IF(KFR1.NE.0) THEN
122 TAUR1=PMAS(KFR1,1)**2/VINT(2)
123 GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
124 MINT(72)=1
125 MINT(73)=KFR1
126 VINT(73)=TAUR1
127 VINT(74)=GAMR1
128 ENDIF
129 IF(ISUB.EQ.141) THEN
130 KFR2=23
131 TAUR2=PMAS(KFR2,1)**2/VINT(2)
132 GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
133 IF(CKIN(1).GT.PMAS(KFR2,1)+20.*PMAS(KFR2,2).OR.
134 & CKMX.LT.PMAS(KFR2,1)-20.*PMAS(KFR2,2)) KFR2=0
135 IF(KFR2.NE.0.AND.KFR1.NE.0) THEN
136 MINT(72)=2
137 MINT(74)=KFR2
138 VINT(75)=TAUR2
139 VINT(76)=GAMR2
140 ELSEIF(KFR2.NE.0) THEN
141 KFR1=KFR2
142 TAUR1=TAUR2
143 GAMR1=GAMR2
144 MINT(72)=1
145 MINT(73)=KFR1
146 VINT(73)=TAUR1
147 VINT(74)=GAMR1
148 ENDIF
149 ENDIF
150
151C...Find product masses and minimum pT of process,
152C...optionally with broadening according to a truncated Breit-Wigner.
153 VINT(63)=0.
154 VINT(64)=0.
155 MINT(71)=0
156 VINT(71)=CKIN(3)
157 IF(MINT(82).GE.2) VINT(71)=0.
158 VINT(80)=1.
159 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4) THEN
160 NBW=0
161 DO 130 I=1,2
162 IF(KFPR(ISUB,I).EQ.0) THEN
163 ELSEIF(MSTP(42).LE.0.OR.PMAS(LUCOMP(KFPR(ISUB,I)),2).LT.
164 & PARP(41)) THEN
165 VINT(62+I)=PMAS(LUCOMP(KFPR(ISUB,I)),1)**2
166 ELSE
167 NBW=NBW+1
168 ENDIF
169 130 CONTINUE
170 IF(NBW.GE.1) THEN
171 CALL PYOFSH(4,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
172 IF(MINT(51).EQ.1) THEN
173 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
174 IF(MFAIL.EQ.1) THEN
175 MSTI(61)=1
176 RETURN
177 ENDIF
178 GOTO 100
179 ENDIF
180 VINT(63)=PQM3**2
181 VINT(64)=PQM4**2
182 ENDIF
183 IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1
184 IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5))
185 ELSEIF(ISTSB.EQ.6) THEN
186 CALL PYOFSH(6,0,KFPR(ISUB,1),KFPR(ISUB,2),0.,PQM3,PQM4)
187 IF(MINT(51).EQ.1) THEN
188 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
189 IF(MFAIL.EQ.1) THEN
190 MSTI(61)=1
191 RETURN
192 ENDIF
193 GOTO 100
194 ENDIF
195 VINT(63)=PQM3**2
196 VINT(64)=PQM4**2
197 ENDIF
198
199C...Prepare for additional variable choices in 2 -> 3.
200 IF(ISTSB.EQ.5) THEN
201 VINT(201)=0.
202 IF(KFPR(ISUB,2).GT.0) VINT(201)=PMAS(KFPR(ISUB,2),1)
203 VINT(206)=VINT(201)
204 VINT(204)=PMAS(23,1)
205 IF(ISUB.EQ.124) VINT(204)=PMAS(24,1)
206 IF(ISUB.EQ.121.OR.ISUB.EQ.122.OR.ISUB.EQ.181.OR.ISUB.EQ.182.OR.
207 & ISUB.EQ.186.OR.ISUB.EQ.187) VINT(204)=VINT(201)
208 VINT(209)=VINT(204)
209 ENDIF
210
211C...Select incoming VDM particle (rho/omega/phi/J/psi).
212 IF(ISTSB.NE.0.AND.(MINT(101).GE.2.OR.MINT(102).GE.2).AND.
213 &(MINT(123).EQ.2.OR.MINT(123).EQ.5.OR.MINT(123).EQ.7)) THEN
214 VRN=RLU(0)*SIGT(0,0,5)
215 IF(MINT(101).LE.1) THEN
216 I1MN=0
217 I1MX=0
218 ELSE
219 I1MN=1
220 I1MX=MINT(101)
221 ENDIF
222 IF(MINT(102).LE.1) THEN
223 I2MN=0
224 I2MX=0
225 ELSE
226 I2MN=1
227 I2MX=MINT(102)
228 ENDIF
229 DO 150 I1=I1MN,I1MX
230 KFV1=110*I1+3
231 DO 140 I2=I2MN,I2MX
232 KFV2=110*I2+3
233 VRN=VRN-SIGT(I1,I2,5)
234 IF(VRN.LE.0.) GOTO 160
235 140 CONTINUE
236 150 CONTINUE
237 160 IF(MINT(101).GE.2) MINT(103)=KFV1
238 IF(MINT(102).GE.2) MINT(104)=KFV2
239 ENDIF
240
241 IF(ISTSB.EQ.0) THEN
242C...Elastic scattering or single or double diffractive scattering.
243
244C...Select incoming particle (rho/omega/phi/J/psi for VDM) and mass.
245 MINT(103)=MINT(11)
246 MINT(104)=MINT(12)
247 PMM(1)=VINT(3)
248 PMM(2)=VINT(4)
249 IF(MINT(101).GE.2.OR.MINT(102).GE.2) THEN
250 JJ=ISUB-90
251 VRN=RLU(0)*SIGT(0,0,JJ)
252 IF(MINT(101).LE.1) THEN
253 I1MN=0
254 I1MX=0
255 ELSE
256 I1MN=1
257 I1MX=MINT(101)
258 ENDIF
259 IF(MINT(102).LE.1) THEN
260 I2MN=0
261 I2MX=0
262 ELSE
263 I2MN=1
264 I2MX=MINT(102)
265 ENDIF
266 DO 180 I1=I1MN,I1MX
267 KFV1=110*I1+3
268 DO 170 I2=I2MN,I2MX
269 KFV2=110*I2+3
270 VRN=VRN-SIGT(I1,I2,JJ)
271 IF(VRN.LE.0.) GOTO 190
272 170 CONTINUE
273 180 CONTINUE
274 190 IF(MINT(101).GE.2) THEN
275 MINT(103)=KFV1
276 PMM(1)=ULMASS(KFV1)
277 ENDIF
278 IF(MINT(102).GE.2) THEN
279 MINT(104)=KFV2
280 PMM(2)=ULMASS(KFV2)
281 ENDIF
282 ENDIF
283
284C...Side/sides of diffractive system.
285 MINT(17)=0
286 MINT(18)=0
287 IF(ISUB.EQ.92.OR.ISUB.EQ.94) MINT(17)=1
288 IF(ISUB.EQ.93.OR.ISUB.EQ.94) MINT(18)=1
289
290C...Find masses of particles and minimal masses of diffractive states.
291 DO 200 JT=1,2
292 PDIF(JT)=PMM(JT)
293 VINT(66+JT)=PDIF(JT)
294 IF(MINT(16+JT).EQ.1) PDIF(JT)=PDIF(JT)+PARP(102)
295 200 CONTINUE
296 SH=VINT(2)
297 SQM1=PMM(1)**2
298 SQM2=PMM(2)**2
299 SQM3=PDIF(1)**2
300 SQM4=PDIF(2)**2
301 SMRES1=(PMM(1)+PMRC)**2
302 SMRES2=(PMM(2)+PMRC)**2
303
304C...Find elastic slope and lower limit diffractive slope.
305 IHA=MAX(2,IABS(MINT(103))/110)
306 IF(IHA.GE.5) IHA=1
307 IHB=MAX(2,IABS(MINT(104))/110)
308 IF(IHB.GE.5) IHB=1
309 IF(ISUB.EQ.91) THEN
310 BMN=2.*BHAD(IHA)+2.*BHAD(IHB)+4.*SH**EPS-4.2
311 ELSEIF(ISUB.EQ.92) THEN
312 BMN=MAX(2.,2.*BHAD(IHB))
313 ELSEIF(ISUB.EQ.93) THEN
314 BMN=MAX(2.,2.*BHAD(IHA))
315 ELSEIF(ISUB.EQ.94) THEN
316 BMN=2.*ALP*4.
317 ENDIF
318
319C...Determine maximum possible t range and coefficient of generation.
320 SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2
321 SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
322 THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
323 THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
324 THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
325 & (SQM1*SQM4-SQM2*SQM3)/SH
326 THL=-0.5*(THA+THB)
327 THU=THC/THL
328 THRND=EXP(MAX(-50.,BMN*(THL-THU)))-1.
329
330C...Select diffractive mass/masses according to dm^2/m^2.
331 210 DO 220 JT=1,2
332 IF(MINT(16+JT).EQ.0) THEN
333 PDIF(2+JT)=PDIF(JT)
334 ELSE
335 PMMIN=PDIF(JT)
336 PMMAX=MAX(VINT(2+JT),VINT(1)-PDIF(3-JT))
337 PDIF(2+JT)=PMMIN*(PMMAX/PMMIN)**RLU(0)
338 ENDIF
339 220 CONTINUE
340 SQM3=PDIF(3)**2
341 SQM4=PDIF(4)**2
342
343C..Additional mass factors, including resonance enhancement.
344 IF(PDIF(3)+PDIF(4).GE.VINT(1)) GOTO 210
345 IF(ISUB.EQ.92) THEN
346 FSD=(1.-SQM3/SH)*(1.+CRES*SMRES1/(SMRES1+SQM3))
347 IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
348 ELSEIF(ISUB.EQ.93) THEN
349 FSD=(1.-SQM4/SH)*(1.+CRES*SMRES2/(SMRES2+SQM4))
350 IF(FSD.LT.RLU(0)*(1.+CRES)) GOTO 210
351 ELSEIF(ISUB.EQ.94) THEN
352 FDD=(1.-(PDIF(3)+PDIF(4))**2/SH)*(SH*SMP/(SH*SMP+SQM3*SQM4))*
353 & (1.+CRES*SMRES1/(SMRES1+SQM3))*(1.+CRES*SMRES2/(SMRES2+SQM4))
354 IF(FDD.LT.RLU(0)*(1.+CRES)**2) GOTO 210
355 ENDIF
356
357C...Select t according to exp(Bmn*t) and correct to right slope.
358 TH=THU+LOG(1.+THRND*RLU(0))/BMN
359 IF(ISUB.GE.92) THEN
360 IF(ISUB.EQ.92) THEN
361 BADD=2.*ALP*LOG(SH/SQM3)
362 IF(BHAD(IHB).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHB)-2.)
363 ELSEIF(ISUB.EQ.93) THEN
364 BADD=2.*ALP*LOG(SH/SQM4)
365 IF(BHAD(IHA).LT.1.) BADD=MAX(0.,BADD+2.*BHAD(IHA)-2.)
366 ELSEIF(ISUB.EQ.94) THEN
367 BADD=2.*ALP*(LOG(EXP(4.)+SH/(ALP*SQM3*SQM4))-4.)
368 ENDIF
369 IF(EXP(MAX(-50.,BADD*(TH-THU))).LT.RLU(0)) GOTO 210
370 ENDIF
371
372C...Check whether m^2 and t choices are consistent.
373 SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4
374 THA=SH-(SQM1+SQM2+SQM3+SQM4)+(SQM1-SQM2)*(SQM3-SQM4)/SH
375 THB=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH
376 IF(THB.LE.1E-8) GOTO 210
377 THC=(SQM3-SQM1)*(SQM4-SQM2)+(SQM1+SQM4-SQM2-SQM3)*
378 & (SQM1*SQM4-SQM2*SQM3)/SH
379 THLM=-0.5*(THA+THB)
380 THUM=THC/THLM
381 IF(TH.LT.THLM.OR.TH.GT.THUM) GOTO 210
382
383C...Information to output.
384 VINT(21)=1.
385 VINT(22)=0.
386 VINT(23)=MIN(1.,MAX(-1.,(THA+2.*TH)/THB))
387 VINT(45)=TH
388 VINT(59)=2.*SQRT(MAX(0.,-(THC+THA*TH+TH**2)))/THB
389 VINT(63)=PDIF(3)**2
390 VINT(64)=PDIF(4)**2
391
392C...Note: in the following, by In is meant the integral over the
393C...quantity multiplying coefficient cn.
394C...Choose tau according to h1(tau)/tau, where
395C...h1(tau) = c1 + I1/I2*c2*1/tau + I1/I3*c3*1/(tau+tau_R) +
396C...I1/I4*c4*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
397C...I1/I5*c5*1/(tau+tau_R') +
398C...I1/I6*c6*tau/((s*tau-m'^2)^2+(m'*Gamma')^2) +
399C...I1/I7*c7*tau/(1.-tau), and
400C...c1 + c2 + c3 + c4 + c5 + c6 + c7 = 1.
401 ELSEIF(ISTSB.GE.1.AND.ISTSB.LE.6) THEN
402 CALL PYKLIM(1)
403 IF(MINT(51).NE.0) THEN
404 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
405 IF(MFAIL.EQ.1) THEN
406 MSTI(61)=1
407 RETURN
408 ENDIF
409 GOTO 100
410 ENDIF
411 RTAU=RLU(0)
412 MTAU=1
413 IF(RTAU.GT.COEF(ISUB,1)) MTAU=2
414 IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3
415 IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4
416 IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4))
417 & MTAU=5
418 IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
419 & COEF(ISUB,5)) MTAU=6
420 IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+
421 & COEF(ISUB,5)+COEF(ISUB,6)) MTAU=7
422 CALL PYKMAP(1,MTAU,RLU(0))
423
424C...2 -> 3, 4 processes:
425C...Choose tau' according to h4(tau,tau')/tau', where
426C...h4(tau,tau') = c1 + I1/I2*c2*(1 - tau/tau')^3/tau' +
427C...I1/I3*c3*1/(1 - tau'), and c1 + c2 + c3 = 1.
428 IF(ISTSB.GE.3.AND.ISTSB.LE.5) THEN
429 CALL PYKLIM(4)
430 IF(MINT(51).NE.0) THEN
431 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
432 IF(MFAIL.EQ.1) THEN
433 MSTI(61)=1
434 RETURN
435 ENDIF
436 GOTO 100
437 ENDIF
438 RTAUP=RLU(0)
439 MTAUP=1
440 IF(RTAUP.GT.COEF(ISUB,18)) MTAUP=2
441 IF(RTAUP.GT.COEF(ISUB,18)+COEF(ISUB,19)) MTAUP=3
442 CALL PYKMAP(4,MTAUP,RLU(0))
443 ENDIF
444
445C...Choose y* according to h2(y*), where
446C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
447C...I0/I3*c3*1/cosh(y*) + I0/I4*c4*1/(1-exp(y*-y*max)) +
448C...I0/I5*c5*1/(1-exp(-y*-y*min)), I0 = y*max-y*min,
449C...and c1 + c2 + c3 + c4 + c5 = 1.
450 CALL PYKLIM(2)
451 IF(MINT(51).NE.0) THEN
452 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
453 IF(MFAIL.EQ.1) THEN
454 MSTI(61)=1
455 RETURN
456 ENDIF
457 GOTO 100
458 ENDIF
459 RYST=RLU(0)
460 MYST=1
461 IF(RYST.GT.COEF(ISUB,8)) MYST=2
462 IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)) MYST=3
463 IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)) MYST=4
464 IF(RYST.GT.COEF(ISUB,8)+COEF(ISUB,9)+COEF(ISUB,10)+
465 & COEF(ISUB,11)) MYST=5
466 CALL PYKMAP(2,MYST,RLU(0))
467
468C...2 -> 2 processes:
469C...Choose cos(theta-hat) (cth) according to h3(cth), where
470C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
471C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
472C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
473C...and c0 + c1 + c2 + c3 + c4 = 1.
474 CALL PYKLIM(3)
475 IF(MINT(51).NE.0) THEN
476 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
477 IF(MFAIL.EQ.1) THEN
478 MSTI(61)=1
479 RETURN
480 ENDIF
481 GOTO 100
482 ENDIF
483 IF(ISTSB.EQ.2.OR.ISTSB.EQ.4.OR.ISTSB.EQ.6) THEN
484 RCTH=RLU(0)
485 MCTH=1
486 IF(RCTH.GT.COEF(ISUB,13)) MCTH=2
487 IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)) MCTH=3
488 IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)) MCTH=4
489 IF(RCTH.GT.COEF(ISUB,13)+COEF(ISUB,14)+COEF(ISUB,15)+
490 & COEF(ISUB,16)) MCTH=5
491 CALL PYKMAP(3,MCTH,RLU(0))
492 ENDIF
493
494C...2 -> 3 : select pT1, phi1, pT2, phi2, y3 for 3 outgoing.
495 IF(ISTSB.EQ.5) THEN
496 CALL PYKMAP(5,0,0.)
497 IF(MINT(51).NE.0) THEN
498 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
499 IF(MFAIL.EQ.1) THEN
500 MSTI(61)=1
501 RETURN
502 ENDIF
503 GOTO 100
504 ENDIF
505 ENDIF
506
507C...Low-pT or multiple interactions (first semihard interaction).
508 ELSEIF(ISTSB.EQ.9) THEN
509 CALL PYMULT(3)
510 ISUB=MINT(1)
511
512C...Generate user-defined process: kinematics plus weight.
513 ELSEIF(ISTSB.EQ.11) THEN
514 MSTI(51)=0
515 CALL PYUPEV(ISUB,SIGS)
516 IF(NUP.LE.0) THEN
517 MINT(51)=2
518 MSTI(51)=1
519 IF(MINT(82).EQ.1) THEN
520 NGEN(0,1)=NGEN(0,1)-1
521 NGEN(0,2)=NGEN(0,2)-1
522 NGEN(ISUB,1)=NGEN(ISUB,1)-1
523 ENDIF
524 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
525 RETURN
526 ENDIF
527
528C...Construct 'trivial' kinematical variables needed.
529 KFL1=KUP(1,2)
530 KFL2=KUP(2,2)
531 VINT(41)=2.*PUP(1,4)/VINT(1)
532 VINT(42)=2.*PUP(2,4)/VINT(1)
533 VINT(21)=VINT(41)*VINT(42)
534 VINT(22)=0.5*LOG(VINT(41)/VINT(42))
535 VINT(44)=VINT(21)*VINT(2)
536 VINT(43)=SQRT(MAX(0.,VINT(44)))
537 VINT(56)=Q2UP(0)
538 VINT(55)=SQRT(MAX(0.,VINT(56)))
539
540C...Construct other kinematical variables needed (approximately).
541 VINT(23)=0.
542 VINT(26)=VINT(21)
543 VINT(45)=-0.5*VINT(44)
544 VINT(46)=-0.5*VINT(44)
545 VINT(49)=VINT(43)
546 VINT(50)=VINT(44)
547 VINT(51)=VINT(55)
548 VINT(52)=VINT(56)
549 VINT(53)=VINT(55)
550 VINT(54)=VINT(56)
551 VINT(25)=0.
552 VINT(48)=0.
553 DO 230 IUP=3,NUP
554 IF(KUP(IUP,1).EQ.1) VINT(25)=VINT(25)+2.*(PUP(IUP,5)**2+
555 & PUP(IUP,1)**2+PUP(IUP,2)**2)/VINT(1)
556 IF(KUP(IUP,1).EQ.1) VINT(48)=VINT(48)+0.5*(PUP(IUP,1)**2+
557 & PUP(IUP,2)**2)
558 230 CONTINUE
559 VINT(47)=SQRT(VINT(48))
560
561C...Calculate structure function weights.
562 IF(MINT(47).GE.2) THEN
563 DO 250 I=3-MIN(2,MINT(45)),MIN(2,MINT(46))
564 MINT(105)=MINT(102+I)
565 MINT(109)=MINT(106+I)
566 IF(MSTP(57).LE.1) THEN
567 CALL PYSTFU(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
568 ELSE
569 CALL PYSTFL(MINT(10+I),VINT(40+I),Q2UP(0),XPQ)
570 ENDIF
571 DO 240 KFL=-25,25
572 XSFX(I,KFL)=XPQ(KFL)
573 240 CONTINUE
574 250 CONTINUE
575 ENDIF
576 ENDIF
577
578C...Choose azimuthal angle.
579 VINT(24)=PARU(2)*RLU(0)
580
581C...Check against user cuts on kinematics at parton level.
582 MINT(51)=0
583 IF((ISUB.LE.90.OR.ISUB.GT.100).AND.ISTSB.LE.10) CALL PYKLIM(0)
584 IF(MINT(51).NE.0) THEN
585 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
586 IF(MFAIL.EQ.1) THEN
587 MSTI(61)=1
588 RETURN
589 ENDIF
590 GOTO 100
591 ENDIF
592 IF(MINT(82).EQ.1.AND.MSTP(141).GE.1.AND.ISTSB.LE.10) THEN
593 MCUT=0
594 IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0)
595 & CALL PYKCUT(MCUT)
596 IF(MCUT.NE.0) THEN
597 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
598 IF(MFAIL.EQ.1) THEN
599 MSTI(61)=1
600 RETURN
601 ENDIF
602 GOTO 100
603 ENDIF
604 ENDIF
605
606C...Calculate differential cross-section for different subprocesses.
607 IF(ISTSB.LE.10) CALL PYSIGH(NCHN,SIGS)
608 SIGSOR=SIGS
609 SIGLPT=SIGT(0,0,5)
610
611C...Multiply cross-section by user-defined weights.
612 IF(MSTP(173).EQ.1) THEN
613 SIGS=PARP(173)*SIGS
614 DO 260 ICHN=1,NCHN
615 SIGH(ICHN)=PARP(173)*SIGH(ICHN)
616 260 CONTINUE
617 SIGLPT=PARP(173)*SIGLPT
618 ENDIF
619 WTXS=1.
620 SIGSWT=SIGS
621 VINT(99)=1.
622 VINT(100)=1.
623 IF(MINT(82).EQ.1.AND.MSTP(142).GE.1) THEN
624 IF(ISUB.NE.96.AND.MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+
625 & MSUB(95).EQ.0) CALL PYEVWT(WTXS)
626 SIGSWT=WTXS*SIGS
627 VINT(99)=WTXS
628 IF(MSTP(142).EQ.1) VINT(100)=1./WTXS
629 ENDIF
630
631C...Calculations for Monte Carlo estimate of all cross-sections.
632 IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN
633 IF(MSTP(142).LE.1) THEN
634 XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
635 DXSEC(ISUB)=DXSEC(ISUB)+SIGS
636 ELSE
637 XSEC(ISUB,2)=XSEC(ISUB,2)+SIGSWT
638 DXSEC(ISUB)=DXSEC(ISUB)+SIGSWT
639 ENDIF
640 ELSEIF(MINT(82).EQ.1) THEN
641 XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS
642 DXSEC(ISUB)=DXSEC(ISUB)+SIGS
643 ENDIF
644 IF((ISUB.EQ.95.OR.ISUB.EQ.96).AND.LOOP.EQ.1.AND.MINT(82).EQ.1)
645 &THEN
646 XSEC(97,2)=XSEC(97,2)+SIGLPT
647 DXSEC(97)=DXSEC(97)+SIGLPT
648 ENDIF
649
650C...Multiple interactions: store results of cross-section calculation.
651 IF(MINT(50).EQ.1.AND.MSTP(82).GE.3) THEN
652 VINT(153)=SIGSOR
653 CALL PYMULT(4)
654 ENDIF
655
656C...Check that weight not negative.
657 VIOL=SIGSWT/XSEC(ISUB,1)
658 IF(ISUB.EQ.96.AND.MSTP(173).EQ.1) VIOL=VIOL/PARP(174)
659 IF(MSTP(123).LE.0) THEN
660 IF(VIOL.LT.-1E-3) THEN
661 WRITE(MSTU(11),5000) VIOL,NGEN(0,3)+1
662 WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
663 STOP
664 ENDIF
665 ELSE
666 IF(VIOL.LT.MIN(-1E-3,VINT(109))) THEN
667 VINT(109)=VIOL
668 WRITE(MSTU(11),5200) VIOL,NGEN(0,3)+1
669 WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
670 ENDIF
671 ENDIF
672
673C...Weighting using estimate of maximum of differential cross-section.
674 IF(MFAIL.EQ.0) THEN
675 IF(VIOL.LT.RLU(0)) THEN
676 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
677 GOTO 100
678 ENDIF
679 ELSEIF(ISUB.NE.95.AND.ISUB.NE.96) THEN
680 IF(VIOL.LT.RLU(0)) THEN
681 MSTI(61)=1
682 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
683 RETURN
684 ENDIF
685 ELSE
686 RATND=SIGLPT/XSEC(95,1)
687 IF(LOOP.EQ.1.AND.RATND.LT.RLU(0)) THEN
688 MSTI(61)=1
689 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
690 RETURN
691 ENDIF
692 VIOL=VIOL/RATND
693 IF(VIOL.LT.RLU(0)) THEN
694 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
695 GOTO 100
696 ENDIF
697 ENDIF
698
699C...Check for possible violation of estimated maximum of differential
700C...cross-section used in weighting.
701 IF(MSTP(123).LE.0) THEN
702 IF(VIOL.GT.1.) THEN
703 WRITE(MSTU(11),5300) VIOL,NGEN(0,3)+1
704 WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
705 STOP
706 ENDIF
707 ELSEIF(MSTP(123).EQ.1) THEN
708 IF(VIOL.GT.VINT(108)) THEN
709 VINT(108)=VIOL
710 IF(VIOL.GT.1.) THEN
711 MINT(10)=1
712 WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
713 WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),
714 & VINT(26)
715 ENDIF
716 ENDIF
717 ELSEIF(VIOL.GT.VINT(108)) THEN
718 VINT(108)=VIOL
719 IF(VIOL.GT.1.) THEN
720 MINT(10)=1
721 XDIF=XSEC(ISUB,1)*(VIOL-1.)
722 XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF
723 IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96))
724 & XSEC(0,1)=XSEC(0,1)+XDIF
725 WRITE(MSTU(11),5400) VIOL,NGEN(0,3)+1
726 WRITE(MSTU(11),5100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
727 IF(ISUB.LE.9) THEN
728 WRITE(MSTU(11),5500) ISUB,XSEC(ISUB,1)
729 ELSEIF(ISUB.LE.99) THEN
730 WRITE(MSTU(11),5600) ISUB,XSEC(ISUB,1)
731 ELSE
732 WRITE(MSTU(11),5700) ISUB,XSEC(ISUB,1)
733 ENDIF
734 VINT(108)=1.
735 ENDIF
736 ENDIF
737
738C...Multiple interactions: choose impact parameter.
739 VINT(148)=1.
740 IF(MINT(50).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
741 &THEN
742 CALL PYMULT(5)
743 IF(VINT(150).LT.RLU(0)) THEN
744 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
745 IF(MFAIL.EQ.1) THEN
746 MSTI(61)=1
747 RETURN
748 ENDIF
749 GOTO 100
750 ENDIF
751 ENDIF
752 IF(MINT(82).EQ.1) NGEN(0,2)=NGEN(0,2)+1
753 IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN
754 IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1
755 IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1
756 ENDIF
757 IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
758
759C...Choose flavour of reacting partons (and subprocess).
760 IF(ISTSB.GE.11) GOTO 280
761 RSIGS=SIGS*RLU(0)
762 QT2=VINT(48)
763 RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2)
764 IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR.
765 &RLU(0).GT.RQQBAR)) THEN
766 DO 270 ICHN=1,NCHN
767 KFL1=ISIG(ICHN,1)
768 KFL2=ISIG(ICHN,2)
769 MINT(2)=ISIG(ICHN,3)
770 RSIGS=RSIGS-SIGH(ICHN)
771 IF(RSIGS.LE.0.) GOTO 280
772 270 CONTINUE
773
774C...Multiple interactions: choose qq~ preferentially at small pT.
775 ELSEIF(ISUB.EQ.96) THEN
776 MINT(105)=MINT(103)
777 MINT(109)=MINT(107)
778 CALL PYSPLI(MINT(11),21,KFL1,KFLDUM)
779 MINT(105)=MINT(104)
780 MINT(109)=MINT(108)
781 CALL PYSPLI(MINT(12),21,KFL2,KFLDUM)
782 MINT(1)=11
783 MINT(2)=1
784 IF(KFL1.EQ.KFL2.AND.RLU(0).LT.0.5) MINT(2)=2
785
786C...Low-pT: choose string drawing configuration.
787 ELSE
788 KFL1=21
789 KFL2=21
790 RSIGS=6.*RLU(0)
791 MINT(2)=1
792 IF(RSIGS.GT.1.) MINT(2)=2
793 IF(RSIGS.GT.2.) MINT(2)=3
794 ENDIF
795
796C...Reassign QCD process. Partons before initial state radiation.
797 280 IF(MINT(2).GT.10) THEN
798 MINT(1)=MINT(2)/10
799 MINT(2)=MOD(MINT(2),10)
800 ENDIF
801 IF(MINT(82).EQ.1.AND.MSTP(111).GE.0) NGEN(MINT(1),2)=
802 &NGEN(MINT(1),2)+1
803 MINT(15)=KFL1
804 MINT(16)=KFL2
805 MINT(13)=MINT(15)
806 MINT(14)=MINT(16)
807 VINT(141)=VINT(41)
808 VINT(142)=VINT(42)
809 VINT(151)=0.
810 VINT(152)=0.
811
812C...Calculate x value of photon for parton inside photon inside e.
813 DO 310 JT=1,2
814 MINT(18+JT)=0
815 VINT(154+JT)=0.
816 MSPLI=0
817 IF(JT.EQ.1.AND.MINT(43).LE.2) MSPLI=1
818 IF(JT.EQ.2.AND.MOD(MINT(43),2).EQ.1) MSPLI=1
819 IF(IABS(MINT(14+JT)).LE.8.OR.MINT(14+JT).EQ.21) MSPLI=MSPLI+1
820 IF(MSPLI.EQ.2) THEN
821 KFLH=MINT(14+JT)
822 XHRD=VINT(140+JT)
823 Q2HRD=VINT(54)
824 MINT(105)=MINT(102+JT)
825 MINT(109)=MINT(106+JT)
826 IF(MSTP(57).LE.1) THEN
827 CALL PYSTFU(22,XHRD,Q2HRD,XPQ)
828 ELSE
829 CALL PYSTFL(22,XHRD,Q2HRD,XPQ)
830 ENDIF
831 WTMX=4.*XPQ(KFLH)
832 IF(MSTP(13).EQ.2) THEN
833 Q2PMS=Q2HRD/PMAS(11,1)**2
834 WTMX=WTMX*LOG(MAX(2.,Q2PMS*(1.-XHRD)/XHRD**2))
835 ENDIF
836 290 XE=XHRD**RLU(0)
837 XG=MIN(0.999999,XHRD/XE)
838 IF(MSTP(57).LE.1) THEN
839 CALL PYSTFU(22,XG,Q2HRD,XPQ)
840 ELSE
841 CALL PYSTFL(22,XG,Q2HRD,XPQ)
842 ENDIF
843 WT=(1.+(1.-XE)**2)*XPQ(KFLH)
844 IF(MSTP(13).EQ.2) WT=WT*LOG(MAX(2.,Q2PMS*(1.-XE)/XE**2))
845 IF(WT.LT.RLU(0)*WTMX) GOTO 290
846 MINT(18+JT)=1
847 VINT(154+JT)=XE
848 DO 300 KFLS=-25,25
849 XSFX(JT,KFLS)=XPQ(KFLS)
850 300 CONTINUE
851 ENDIF
852 310 CONTINUE
853
854C...Pick scale where photon is resolved.
855 IF(MINT(107).EQ.3) VINT(283)=PARP(15)**2*
856 &(VINT(54)/PARP(15)**2)**RLU(0)
857 IF(MINT(108).EQ.3) VINT(284)=PARP(15)**2*
858 &(VINT(54)/PARP(15)**2)**RLU(0)
859 IF(MINT(121).GT.1) CALL PYSAVE(2,IGA)
860
861C...Format statements for differential cross-section maximum violations.
862 5000 FORMAT(1X,'Error: negative cross-section fraction',1P,E11.3,1X,
863 &'in event',1X,I7,'.'/1X,'Execution stopped!')
864 5100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P,
865 &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3)
866 5200 FORMAT(1X,'Warning: negative cross-section fraction',1P,E11.3,1X,
867 &'in event',1X,I7)
868 5300 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X,
869 &'in event',1X,I7,'.'/1X,'Execution stopped!')
870 5400 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X,
871 &'in event',1X,I7)
872 5500 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3)
873 5600 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3)
874 5700 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3)
875
876 RETURN
877 END