]>
Commit | Line | Data |
---|---|---|
e74335a4 | 1 | * $Id$ |
2 | ||
3 | C********************************************************************* | |
4 | ||
5 | SUBROUTINE PYRAND_HIJING | |
6 | ||
7 | C...Generates quantities characterizing the high-pT scattering at the | |
8 | C...parton level according to the matrix elements. Chooses incoming, | |
9 | C...reacting partons, their momentum fractions and one of the possible | |
10 | C...subprocesses. | |
11 | #include "ludat1_hijing.inc" | |
12 | #include "ludat2_hijing.inc" | |
13 | #include "pysubs_hijing.inc" | |
14 | #include "pypars_hijing.inc" | |
15 | #include "pyint1_hijing.inc" | |
16 | #include "pyint2_hijing.inc" | |
17 | #include "pyint3_hijing.inc" | |
18 | #include "pyint4_hijing.inc" | |
19 | #include "pyint5_hijing.inc" | |
20 | ||
21 | C...Initial values, specifically for (first) semihard interaction. | |
22 | MINT(17)=0 | |
23 | MINT(18)=0 | |
24 | VINT(143)=1. | |
25 | VINT(144)=1. | |
26 | IF(MSUB(95).EQ.1.OR.MINT(82).GE.2) CALL PYMULT_HIJING(2) | |
27 | ISUB=0 | |
28 | 100 MINT(51)=0 | |
29 | ||
30 | C...Choice of process type - first event of overlay. | |
31 | IF(MINT(82).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) THEN | |
32 | RSUB=XSEC(0,1)*RLU_HIJING(0) | |
33 | DO 110 I=1,200 | |
34 | IF(MSUB(I).NE.1) GOTO 110 | |
35 | ISUB=I | |
36 | RSUB=RSUB-XSEC(I,1) | |
37 | IF(RSUB.LE.0.) GOTO 120 | |
38 | 110 CONTINUE | |
39 | 120 IF(ISUB.EQ.95) ISUB=96 | |
40 | ||
41 | C...Choice of inclusive process type - overlayed events. | |
42 | ELSEIF(MINT(82).GE.2.AND.ISUB.EQ.0) THEN | |
43 | RSUB=VINT(131)*RLU_HIJING(0) | |
44 | ISUB=96 | |
45 | IF(RSUB.GT.VINT(106)) ISUB=93 | |
46 | IF(RSUB.GT.VINT(106)+VINT(104)) ISUB=92 | |
47 | IF(RSUB.GT.VINT(106)+VINT(104)+VINT(103)) ISUB=91 | |
48 | ENDIF | |
49 | IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)+1 | |
50 | IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)+1 | |
51 | MINT(1)=ISUB | |
52 | ||
53 | C...Find resonances (explicit or implicit in cross-section). | |
54 | MINT(72)=0 | |
55 | KFR1=0 | |
56 | IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN | |
57 | KFR1=KFPR(ISUB,1) | |
58 | ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN | |
59 | KFR1=25 | |
60 | ENDIF | |
61 | IF(KFR1.NE.0) THEN | |
62 | TAUR1=PMAS(KFR1,1)**2/VINT(2) | |
63 | GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2) | |
64 | MINT(72)=1 | |
65 | MINT(73)=KFR1 | |
66 | VINT(73)=TAUR1 | |
67 | VINT(74)=GAMR1 | |
68 | ENDIF | |
69 | IF(ISUB.EQ.141) THEN | |
70 | KFR2=23 | |
71 | TAUR2=PMAS(KFR2,1)**2/VINT(2) | |
72 | GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2) | |
73 | MINT(72)=2 | |
74 | MINT(74)=KFR2 | |
75 | VINT(75)=TAUR2 | |
76 | VINT(76)=GAMR2 | |
77 | ENDIF | |
78 | ||
79 | C...Find product masses and minimum pT of process, | |
80 | C...optionally with broadening according to a truncated Breit-Wigner. | |
81 | VINT(63)=0. | |
82 | VINT(64)=0. | |
83 | MINT(71)=0 | |
84 | VINT(71)=CKIN(3) | |
85 | IF(MINT(82).GE.2) VINT(71)=0. | |
86 | IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN | |
87 | DO 130 I=1,2 | |
88 | IF(KFPR(ISUB,I).EQ.0) THEN | |
89 | ELSEIF(MSTP(42).LE.0) THEN | |
90 | VINT(62+I)=PMAS(KFPR(ISUB,I),1)**2 | |
91 | ELSE | |
92 | VINT(62+I)=ULMASS_HIJING(KFPR(ISUB,I))**2 | |
93 | ENDIF | |
94 | 130 CONTINUE | |
95 | IF(MIN(VINT(63),VINT(64)).LT.CKIN(6)**2) MINT(71)=1 | |
96 | IF(MINT(71).EQ.1) VINT(71)=MAX(CKIN(3),CKIN(5)) | |
97 | ENDIF | |
98 | ||
99 | IF(ISET(ISUB).EQ.0) THEN | |
100 | C...Double or single diffractive, or elastic scattering: | |
101 | C...choose m^2 according to 1/m^2 (diffractive), constant (elastic) | |
102 | IS=INT(1.5+RLU_HIJING(0)) | |
103 | VINT(63)=VINT(3)**2 | |
104 | VINT(64)=VINT(4)**2 | |
105 | IF(ISUB.EQ.92.OR.ISUB.EQ.93) VINT(62+IS)=PARP(111)**2 | |
106 | IF(ISUB.EQ.93) VINT(65-IS)=PARP(111)**2 | |
107 | SH=VINT(2) | |
108 | SQM1=VINT(3)**2 | |
109 | SQM2=VINT(4)**2 | |
110 | SQM3=VINT(63) | |
111 | SQM4=VINT(64) | |
112 | SQLA12=(SH-SQM1-SQM2)**2-4.*SQM1*SQM2 | |
113 | SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4 | |
114 | THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH | |
115 | THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH | |
116 | THL=0.5*(THTER1-THTER2) | |
117 | THU=0.5*(THTER1+THTER2) | |
118 | THM=MIN(MAX(THL,PARP(101)),THU) | |
119 | JTMAX=0 | |
120 | IF(ISUB.EQ.92.OR.ISUB.EQ.93) JTMAX=ISUB-91 | |
121 | DO 140 JT=1,JTMAX | |
122 | MINT(13+3*JT-IS*(2*JT-3))=1 | |
123 | SQMMIN=VINT(59+3*JT-IS*(2*JT-3)) | |
124 | SQMI=VINT(8-3*JT+IS*(2*JT-3))**2 | |
125 | SQMJ=VINT(3*JT-1-IS*(2*JT-3))**2 | |
126 | SQMF=VINT(68-3*JT+IS*(2*JT-3)) | |
127 | SQUA=0.5*SH/SQMI*((1.+(SQMI-SQMJ)/SH)*THM+SQMI-SQMF- | |
128 | & SQMJ**2/SH+(SQMI+SQMJ)*SQMF/SH+(SQMI-SQMJ)**2/SH**2*SQMF) | |
129 | QUAR=SH/SQMI*(THM*(THM+SH-SQMI-SQMJ-SQMF*(1.-(SQMI-SQMJ)/SH))+ | |
130 | & SQMI*SQMJ-SQMJ*SQMF*(1.+(SQMI-SQMJ-SQMF)/SH)) | |
131 | SQMMAX=SQUA+SQRT(MAX(0.,SQUA**2-QUAR)) | |
132 | IF(ABS(QUAR/SQUA**2).LT.1.E-06) SQMMAX=0.5*QUAR/SQUA | |
133 | SQMMAX=MIN(SQMMAX,(VINT(1)-SQRT(SQMF))**2) | |
134 | VINT(59+3*JT-IS*(2*JT-3))=SQMMIN*(SQMMAX/SQMMIN)**RLU_HIJING(0) | |
135 | 140 CONTINUE | |
136 | C...Choose t-hat according to exp(B*t-hat+C*t-hat^2). | |
137 | SQM3=VINT(63) | |
138 | SQM4=VINT(64) | |
139 | SQLA34=(SH-SQM3-SQM4)**2-4.*SQM3*SQM4 | |
140 | THTER1=SQM1+SQM2+SQM3+SQM4-(SQM1-SQM2)*(SQM3-SQM4)/SH-SH | |
141 | THTER2=SQRT(MAX(0.,SQLA12))*SQRT(MAX(0.,SQLA34))/SH | |
142 | THL=0.5*(THTER1-THTER2) | |
143 | THU=0.5*(THTER1+THTER2) | |
144 | B=VINT(121) | |
145 | C=VINT(122) | |
146 | IF(ISUB.EQ.92.OR.ISUB.EQ.93) THEN | |
147 | B=0.5*B | |
148 | C=0.5*C | |
149 | ENDIF | |
150 | THM=MIN(MAX(THL,PARP(101)),THU) | |
151 | EXPTH=0. | |
152 | THARG=B*(THM-THU) | |
153 | IF(THARG.GT.-20.) EXPTH=EXP(THARG) | |
154 | 150 TH=THU+LOG(EXPTH+(1.-EXPTH)*RLU_HIJING(0))/B | |
155 | TH=MAX(THM,MIN(THU,TH)) | |
156 | RATLOG=MIN((B+C*(TH+THM))*(TH-THM),(B+C*(TH+THU))*(TH-THU)) | |
157 | IF(RATLOG.LT.LOG(RLU_HIJING(0))) GOTO 150 | |
158 | VINT(21)=1. | |
159 | VINT(22)=0. | |
160 | VINT(23)=MIN(1.,MAX(-1.,(2.*TH-THTER1)/THTER2)) | |
161 | ||
162 | C...Note: in the following, by In is meant the integral over the | |
163 | C...quantity multiplying coefficient cn. | |
164 | C...Choose tau according to h1(tau)/tau, where | |
165 | C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) + | |
166 | C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) + | |
167 | C...I0/I4*c4*1/(tau+tau_R') + | |
168 | C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and | |
169 | C...c0 + c1 + c2 + c3 + c4 + c5 = 1 | |
170 | ELSEIF(ISET(ISUB).GE.1.AND.ISET(ISUB).LE.4) THEN | |
171 | CALL PYKLIM_HIJING(1) | |
172 | IF(MINT(51).NE.0) GOTO 100 | |
173 | RTAU=RLU_HIJING(0) | |
174 | MTAU=1 | |
175 | IF(RTAU.GT.COEF(ISUB,1)) MTAU=2 | |
176 | IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)) MTAU=3 | |
177 | IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)) MTAU=4 | |
178 | IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)) | |
179 | & MTAU=5 | |
180 | IF(RTAU.GT.COEF(ISUB,1)+COEF(ISUB,2)+COEF(ISUB,3)+COEF(ISUB,4)+ | |
181 | & COEF(ISUB,5)) MTAU=6 | |
182 | CALL PYKMAP_HIJING(1,MTAU,RLU_HIJING(0)) | |
183 | ||
184 | C...2 -> 3, 4 processes: | |
185 | C...Choose tau' according to h4(tau,tau')/tau', where | |
186 | C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and | |
187 | C...c0 + c1 = 1. | |
188 | IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN | |
189 | CALL PYKLIM_HIJING(4) | |
190 | IF(MINT(51).NE.0) GOTO 100 | |
191 | RTAUP=RLU_HIJING(0) | |
192 | MTAUP=1 | |
193 | IF(RTAUP.GT.COEF(ISUB,15)) MTAUP=2 | |
194 | CALL PYKMAP_HIJING(4,MTAUP,RLU_HIJING(0)) | |
195 | ENDIF | |
196 | ||
197 | C...Choose y* according to h2(y*), where | |
198 | C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) + | |
199 | C...I0/I3*c3*1/cosh(y*), I0 = y*max-y*min, and c1 + c2 + c3 = 1. | |
200 | CALL PYKLIM_HIJING(2) | |
201 | IF(MINT(51).NE.0) GOTO 100 | |
202 | RYST=RLU_HIJING(0) | |
203 | MYST=1 | |
204 | IF(RYST.GT.COEF(ISUB,7)) MYST=2 | |
205 | IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3 | |
206 | CALL PYKMAP_HIJING(2,MYST,RLU_HIJING(0)) | |
207 | ||
208 | C...2 -> 2 processes: | |
209 | C...Choose cos(theta-hat) (cth) according to h3(cth), where | |
210 | C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) + | |
211 | C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2, | |
212 | C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products), | |
213 | C...and c0 + c1 + c2 + c3 + c4 = 1. | |
214 | CALL PYKLIM_HIJING(3) | |
215 | IF(MINT(51).NE.0) GOTO 100 | |
216 | IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN | |
217 | RCTH=RLU_HIJING(0) | |
218 | MCTH=1 | |
219 | IF(RCTH.GT.COEF(ISUB,10)) MCTH=2 | |
220 | IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)) MCTH=3 | |
221 | IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)) MCTH=4 | |
222 | IF(RCTH.GT.COEF(ISUB,10)+COEF(ISUB,11)+COEF(ISUB,12)+ | |
223 | & COEF(ISUB,13)) MCTH=5 | |
224 | CALL PYKMAP_HIJING(3,MCTH,RLU_HIJING(0)) | |
225 | ENDIF | |
226 | ||
227 | C...Low-pT or multiple interactions (first semihard interaction). | |
228 | ELSEIF(ISET(ISUB).EQ.5) THEN | |
229 | CALL PYMULT_HIJING(3) | |
230 | ISUB=MINT(1) | |
231 | ENDIF | |
232 | ||
233 | C...Choose azimuthal angle. | |
234 | VINT(24)=PARU(2)*RLU_HIJING(0) | |
235 | ||
236 | C...Check against user cuts on kinematics at parton level. | |
237 | MINT(51)=0 | |
238 | IF(ISUB.LE.90.OR.ISUB.GT.100) CALL PYKLIM_HIJING(0) | |
239 | IF(MINT(51).NE.0) GOTO 100 | |
240 | IF(MINT(82).EQ.1.AND.MSTP(141).GE.1) THEN | |
241 | MCUT=0 | |
242 | IF(MSUB(91)+MSUB(92)+MSUB(93)+MSUB(94)+MSUB(95).EQ.0) | |
243 | & CALL PYKCUT_HIJING(MCUT) | |
244 | IF(MCUT.NE.0) GOTO 100 | |
245 | ENDIF | |
246 | ||
247 | C...Calculate differential cross-section for different subprocesses. | |
248 | CALL PYSIGH_HIJING(NCHN,SIGS) | |
249 | ||
250 | C...Calculations for Monte Carlo estimate of all cross-sections. | |
251 | IF(MINT(82).EQ.1.AND.ISUB.LE.90.OR.ISUB.GE.96) THEN | |
252 | XSEC(ISUB,2)=XSEC(ISUB,2)+SIGS | |
253 | ELSEIF(MINT(82).EQ.1) THEN | |
254 | XSEC(ISUB,2)=XSEC(ISUB,2)+XSEC(ISUB,1) | |
255 | ENDIF | |
256 | ||
257 | C...Multiple interactions: store results of cross-section calculation. | |
258 | IF(MINT(43).EQ.4.AND.MSTP(82).GE.3) THEN | |
259 | VINT(153)=SIGS | |
260 | CALL PYMULT_HIJING(4) | |
261 | ENDIF | |
262 | ||
263 | C...Weighting using estimate of maximum of differential cross-section. | |
264 | VIOL=SIGS/XSEC(ISUB,1) | |
265 | IF(VIOL.LT.RLU_HIJING(0)) GOTO 100 | |
266 | ||
267 | C...Check for possible violation of estimated maximum of differential | |
268 | C...cross-section used in weighting. | |
269 | IF(MSTP(123).LE.0) THEN | |
270 | IF(VIOL.GT.1.) THEN | |
271 | WRITE(MSTU(11),1000) VIOL,NGEN(0,3)+1 | |
272 | WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26) | |
273 | STOP | |
274 | ENDIF | |
275 | ELSEIF(MSTP(123).EQ.1) THEN | |
276 | IF(VIOL.GT.VINT(108)) THEN | |
277 | VINT(108)=VIOL | |
278 | C IF(VIOL.GT.1.) THEN | |
279 | C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1 | |
280 | C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23), | |
281 | C & VINT(26) | |
282 | C ENDIF | |
283 | ENDIF | |
284 | ELSEIF(VIOL.GT.VINT(108)) THEN | |
285 | VINT(108)=VIOL | |
286 | IF(VIOL.GT.1.) THEN | |
287 | XDIF=XSEC(ISUB,1)*(VIOL-1.) | |
288 | XSEC(ISUB,1)=XSEC(ISUB,1)+XDIF | |
289 | IF(MSUB(ISUB).EQ.1.AND.(ISUB.LE.90.OR.ISUB.GT.96)) | |
290 | & XSEC(0,1)=XSEC(0,1)+XDIF | |
291 | C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1 | |
292 | C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26) | |
293 | C IF(ISUB.LE.9) THEN | |
294 | C WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1) | |
295 | C ELSEIF(ISUB.LE.99) THEN | |
296 | C WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1) | |
297 | C ELSE | |
298 | C WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1) | |
299 | C ENDIF | |
300 | VINT(108)=1. | |
301 | ENDIF | |
302 | ENDIF | |
303 | ||
304 | C...Multiple interactions: choose impact parameter. | |
305 | VINT(148)=1. | |
306 | IF(MINT(43).EQ.4.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3) | |
307 | &THEN | |
308 | CALL PYMULT_HIJING(5) | |
309 | IF(VINT(150).LT.RLU_HIJING(0)) GOTO 100 | |
310 | ENDIF | |
311 | IF(MINT(82).EQ.1.AND.MSUB(95).EQ.1) THEN | |
312 | IF(ISUB.LE.90.OR.ISUB.GE.95) NGEN(95,1)=NGEN(95,1)+1 | |
313 | IF(ISUB.LE.90.OR.ISUB.GE.96) NGEN(96,2)=NGEN(96,2)+1 | |
314 | ENDIF | |
315 | IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1 | |
316 | ||
317 | C...Choose flavour of reacting partons (and subprocess). | |
318 | RSIGS=SIGS*RLU_HIJING(0) | |
319 | QT2=VINT(48) | |
320 | RQQBAR=PARP(87)*(1.-(QT2/(QT2+(PARP(88)*PARP(82))**2))**2) | |
321 | IF(ISUB.NE.95.AND.(ISUB.NE.96.OR.MSTP(82).LE.1.OR. | |
322 | &RLU_HIJING(0).GT.RQQBAR)) THEN | |
323 | DO 190 ICHN=1,NCHN | |
324 | KFL1=ISIG(ICHN,1) | |
325 | KFL2=ISIG(ICHN,2) | |
326 | MINT(2)=ISIG(ICHN,3) | |
327 | RSIGS=RSIGS-SIGH(ICHN) | |
328 | IF(RSIGS.LE.0.) GOTO 210 | |
329 | 190 CONTINUE | |
330 | ||
331 | C...Multiple interactions: choose qqbar preferentially at small pT. | |
332 | ELSEIF(ISUB.EQ.96) THEN | |
333 | CALL PYSPLI_HIJING(MINT(11),21,KFL1,KFLDUM) | |
334 | CALL PYSPLI_HIJING(MINT(12),21,KFL2,KFLDUM) | |
335 | MINT(1)=11 | |
336 | MINT(2)=1 | |
337 | IF(KFL1.EQ.KFL2.AND.RLU_HIJING(0).LT.0.5) MINT(2)=2 | |
338 | ||
339 | C...Low-pT: choose string drawing configuration. | |
340 | ELSE | |
341 | KFL1=21 | |
342 | KFL2=21 | |
343 | RSIGS=6.*RLU_HIJING(0) | |
344 | MINT(2)=1 | |
345 | IF(RSIGS.GT.1.) MINT(2)=2 | |
346 | IF(RSIGS.GT.2.) MINT(2)=3 | |
347 | ENDIF | |
348 | ||
349 | C...Reassign QCD process. Partons before initial state radiation. | |
350 | 210 IF(MINT(2).GT.10) THEN | |
351 | MINT(1)=MINT(2)/10 | |
352 | MINT(2)=MOD(MINT(2),10) | |
353 | ENDIF | |
354 | MINT(15)=KFL1 | |
355 | MINT(16)=KFL2 | |
356 | MINT(13)=MINT(15) | |
357 | MINT(14)=MINT(16) | |
358 | VINT(141)=VINT(41) | |
359 | VINT(142)=VINT(42) | |
360 | ||
361 | C...Format statements for differential cross-section maximum violations. | |
362 | 1000 FORMAT(1X,'Error: maximum violated by',1P,E11.3,1X, | |
363 | &'in event',1X,I7,'.'/1X,'Execution stopped!') | |
364 | 1100 FORMAT(1X,'ISUB = ',I3,'; Point of violation:'/1X,'tau =',1P, | |
365 | &E11.3,', y* =',E11.3,', cthe = ',0P,F11.7,', tau'' =',1P,E11.3) | |
366 | 1200 FORMAT(1X,'Warning: maximum violated by',1P,E11.3,1X, | |
367 | &'in event',1X,I7) | |
368 | 1300 FORMAT(1X,'XSEC(',I1,',1) increased to',1P,E11.3) | |
369 | 1400 FORMAT(1X,'XSEC(',I2,',1) increased to',1P,E11.3) | |
370 | 1500 FORMAT(1X,'XSEC(',I3,',1) increased to',1P,E11.3) | |
371 | ||
372 | RETURN | |
373 | END |