3 C*********************************************************************
5 SUBROUTINE PYRAND_HIJING
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
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"
21 C...Initial values, specifically for (first) semihard interaction.
26 IF(MSUB(95).EQ.1.OR.MINT(82).GE.2) CALL PYMULT_HIJING(2)
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)
34 IF(MSUB(I).NE.1) GOTO 110
37 IF(RSUB.LE.0.) GOTO 120
39 120 IF(ISUB.EQ.95) ISUB=96
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)
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
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
53 C...Find resonances (explicit or implicit in cross-section).
56 IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
58 ELSEIF(ISUB.GE.71.AND.ISUB.LE.77) THEN
62 TAUR1=PMAS(KFR1,1)**2/VINT(2)
63 GAMR1=PMAS(KFR1,1)*PMAS(KFR1,2)/VINT(2)
71 TAUR2=PMAS(KFR2,1)**2/VINT(2)
72 GAMR2=PMAS(KFR2,1)*PMAS(KFR2,2)/VINT(2)
79 C...Find product masses and minimum pT of process,
80 C...optionally with broadening according to a truncated Breit-Wigner.
85 IF(MINT(82).GE.2) VINT(71)=0.
86 IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN
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
92 VINT(62+I)=ULMASS_HIJING(KFPR(ISUB,I))**2
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))
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))
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
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)
120 IF(ISUB.EQ.92.OR.ISUB.EQ.93) JTMAX=ISUB-91
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)
136 C...Choose t-hat according to exp(B*t-hat+C*t-hat^2).
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)
146 IF(ISUB.EQ.92.OR.ISUB.EQ.93) THEN
150 THM=MIN(MAX(THL,PARP(101)),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
160 VINT(23)=MIN(1.,MAX(-1.,(2.*TH-THTER1)/THTER2))
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
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))
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))
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
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
193 IF(RTAUP.GT.COEF(ISUB,15)) MTAUP=2
194 CALL PYKMAP_HIJING(4,MTAUP,RLU_HIJING(0))
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
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))
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
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))
227 C...Low-pT or multiple interactions (first semihard interaction).
228 ELSEIF(ISET(ISUB).EQ.5) THEN
229 CALL PYMULT_HIJING(3)
233 C...Choose azimuthal angle.
234 VINT(24)=PARU(2)*RLU_HIJING(0)
236 C...Check against user cuts on kinematics at parton level.
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
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
247 C...Calculate differential cross-section for different subprocesses.
248 CALL PYSIGH_HIJING(NCHN,SIGS)
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)
257 C...Multiple interactions: store results of cross-section calculation.
258 IF(MINT(43).EQ.4.AND.MSTP(82).GE.3) THEN
260 CALL PYMULT_HIJING(4)
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
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
271 WRITE(MSTU(11),1000) VIOL,NGEN(0,3)+1
272 WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
275 ELSEIF(MSTP(123).EQ.1) THEN
276 IF(VIOL.GT.VINT(108)) THEN
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),
284 ELSEIF(VIOL.GT.VINT(108)) 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)
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)
298 C WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)
304 C...Multiple interactions: choose impact parameter.
306 IF(MINT(43).EQ.4.AND.(ISUB.LE.90.OR.ISUB.GE.96).AND.MSTP(82).GE.3)
308 CALL PYMULT_HIJING(5)
309 IF(VINT(150).LT.RLU_HIJING(0)) GOTO 100
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
315 IF(ISUB.LE.90.OR.ISUB.GE.96) MINT(31)=MINT(31)+1
317 C...Choose flavour of reacting partons (and subprocess).
318 RSIGS=SIGS*RLU_HIJING(0)
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
327 RSIGS=RSIGS-SIGH(ICHN)
328 IF(RSIGS.LE.0.) GOTO 210
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)
337 IF(KFL1.EQ.KFL2.AND.RLU_HIJING(0).LT.0.5) MINT(2)=2
339 C...Low-pT: choose string drawing configuration.
343 RSIGS=6.*RLU_HIJING(0)
345 IF(RSIGS.GT.1.) MINT(2)=2
346 IF(RSIGS.GT.2.) MINT(2)=3
349 C...Reassign QCD process. Partons before initial state radiation.
350 210 IF(MINT(2).GT.10) THEN
352 MINT(2)=MOD(MINT(2),10)
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,
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)