]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/pyrand_hijing.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
- Reset TProcessID count after each event
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pyrand_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE PYRAND_HIJING
6
7C...Generates quantities characterizing the high-pT scattering at the
8C...parton level according to the matrix elements. Chooses incoming,
9C...reacting partons, their momentum fractions and one of the possible
10C...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
21C...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
30C...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
41C...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
53C...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
79C...Find product masses and minimum pT of process,
80C...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
100C...Double or single diffractive, or elastic scattering:
101C...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
136C...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
162C...Note: in the following, by In is meant the integral over the
163C...quantity multiplying coefficient cn.
164C...Choose tau according to h1(tau)/tau, where
165C...h1(tau) = c0 + I0/I1*c1*1/tau + I0/I2*c2*1/(tau+tau_R) +
166C...I0/I3*c3*tau/((s*tau-m^2)^2+(m*Gamma)^2) +
167C...I0/I4*c4*1/(tau+tau_R') +
168C...I0/I5*c5*tau/((s*tau-m'^2)^2+(m'*Gamma')^2), and
169C...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
184C...2 -> 3, 4 processes:
185C...Choose tau' according to h4(tau,tau')/tau', where
186C...h4(tau,tau') = c0 + I0/I1*c1*(1 - tau/tau')^3/tau', and
187C...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
197C...Choose y* according to h2(y*), where
198C...h2(y*) = I0/I1*c1*(y*-y*min) + I0/I2*c2*(y*max-y*) +
199C...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
208C...2 -> 2 processes:
209C...Choose cos(theta-hat) (cth) according to h3(cth), where
210C...h3(cth) = c0 + I0/I1*c1*1/(A - cth) + I0/I2*c2*1/(A + cth) +
211C...I0/I3*c3*1/(A - cth)^2 + I0/I4*c4*1/(A + cth)^2,
212C...A = 1 + 2*(m3*m4/sh)^2 (= 1 for massless products),
213C...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
227C...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
233C...Choose azimuthal angle.
234 VINT(24)=PARU(2)*RLU_HIJING(0)
235
236C...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
247C...Calculate differential cross-section for different subprocesses.
248 CALL PYSIGH_HIJING(NCHN,SIGS)
249
250C...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
257C...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
263C...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
267C...Check for possible violation of estimated maximum of differential
268C...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
278C IF(VIOL.GT.1.) THEN
279C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
280C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),
281C & VINT(26)
282C 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
291C WRITE(MSTU(11),1200) VIOL,NGEN(0,3)+1
292C WRITE(MSTU(11),1100) ISUB,VINT(21),VINT(22),VINT(23),VINT(26)
293C IF(ISUB.LE.9) THEN
294C WRITE(MSTU(11),1300) ISUB,XSEC(ISUB,1)
295C ELSEIF(ISUB.LE.99) THEN
296C WRITE(MSTU(11),1400) ISUB,XSEC(ISUB,1)
297C ELSE
298C WRITE(MSTU(11),1500) ISUB,XSEC(ISUB,1)
299C ENDIF
300 VINT(108)=1.
301 ENDIF
302 ENDIF
303
304C...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
317C...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
331C...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
339C...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
349C...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
361C...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