]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hipyset1_35/pymult_hijing.F
- Reset TProcessID count after each event
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pymult_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE PYMULT_HIJING(MMUL)
6
7C...Initializes treatment of multiple interactions, selects kinematics
8C...of hardest interaction if low-pT physics included in run, and
9C...generates all non-hardest interactions.
10#include "lujets_hijing.inc"
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 "pyint5_hijing.inc"
19 DIMENSION NMUL(20),SIGM(20),KSTR(500,2)
20 SAVE XT2,XT2FAC,XC2,XTS,IRBIN,RBIN,NMUL,SIGM
21
22C...Initialization of multiple interaction treatment.
23 IF(MMUL.EQ.1) THEN
24 IF(MSTP(122).GE.1) WRITE(MSTU(11),1000) MSTP(82)
25 ISUB=96
26 MINT(1)=96
27 VINT(63)=0.
28 VINT(64)=0.
29 VINT(143)=1.
30 VINT(144)=1.
31
32C...Loop over phase space points: xT2 choice in 20 bins.
33 100 SIGSUM=0.
34 DO 120 IXT2=1,20
35 NMUL(IXT2)=MSTP(83)
36 SIGM(IXT2)=0.
37 DO 110 ITRY=1,MSTP(83)
38 RSCA=0.05*((21-IXT2)-RLU_HIJING(0))
39 XT2=VINT(149)*(1.+VINT(149))/(VINT(149)+RSCA)-VINT(149)
40 XT2=MAX(0.01*VINT(149),XT2)
41 VINT(25)=XT2
42
43C...Choose tau and y*. Calculate cos(theta-hat).
44 IF(RLU_HIJING(0).LE.COEF(ISUB,1)) THEN
45 TAUP=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU_HIJING(0)
46 TAU=XT2*(1.+TAUP)**2/(4.*TAUP)
47 ELSE
48 TAU=XT2*(1.+TAN(RLU_HIJING(0)*ATAN(SQRT(1./XT2-1.)))**2)
49 ENDIF
50 VINT(21)=TAU
51 CALL PYKLIM_HIJING(2)
52 RYST=RLU_HIJING(0)
53 MYST=1
54 IF(RYST.GT.COEF(ISUB,7)) MYST=2
55 IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3
56 CALL PYKMAP_HIJING(2,MYST,RLU_HIJING(0))
57 VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU_HIJING(0))
58
59C...Calculate differential cross-section.
60 VINT(71)=0.5*VINT(1)*SQRT(XT2)
61 CALL PYSIGH_HIJING(NCHN,SIGS)
62 110 SIGM(IXT2)=SIGM(IXT2)+SIGS
63 120 SIGSUM=SIGSUM+SIGM(IXT2)
64 SIGSUM=SIGSUM/(20.*MSTP(83))
65
66C...Reject result if sigma(parton-parton) is smaller than hadronic one.
67 IF(SIGSUM.LT.1.1*VINT(106)) THEN
68 IF(MSTP(122).GE.1) WRITE(MSTU(11),1100) PARP(82),SIGSUM
69 PARP(82)=0.9*PARP(82)
70 VINT(149)=4.*PARP(82)**2/VINT(2)
71 GOTO 100
72 ENDIF
73 IF(MSTP(122).GE.1) WRITE(MSTU(11),1200) PARP(82), SIGSUM
74
75C...Start iteration to find k factor.
76 YKE=SIGSUM/VINT(106)
77 SO=0.5
78 XI=0.
79 YI=0.
80 XK=0.5
81 IIT=0
82 130 IF(IIT.EQ.0) THEN
83 XK=2.*XK
84 ELSEIF(IIT.EQ.1) THEN
85 XK=0.5*XK
86 ELSE
87 XK=XI+(YKE-YI)*(XF-XI)/(YF-YI)
88 ENDIF
89
90C...Evaluate overlap integrals.
91 IF(MSTP(82).EQ.2) THEN
92 SP=0.5*PARU(1)*(1.-EXP(-XK))
93 SOP=SP/PARU(1)
94 ELSE
95 IF(MSTP(82).EQ.3) DELTAB=0.02
96 IF(MSTP(82).EQ.4) DELTAB=MIN(0.01,0.05*PARP(84))
97 SP=0.
98 SOP=0.
99 B=-0.5*DELTAB
100 140 B=B+DELTAB
101 IF(MSTP(82).EQ.3) THEN
102 OV=EXP(-B**2)/PARU(2)
103 ELSE
104 CQ2=PARP(84)**2
105 OV=((1.-PARP(83))**2*EXP(-MIN(100.,B**2))+2.*PARP(83)*
106 & (1.-PARP(83))*2./(1.+CQ2)*EXP(-MIN(100.,B**2*2./(1.+CQ2)))+
107 & PARP(83)**2/CQ2*EXP(-MIN(100.,B**2/CQ2)))/PARU(2)
108 ENDIF
109 PACC=1.-EXP(-MIN(100.,PARU(1)*XK*OV))
110 SP=SP+PARU(2)*B*DELTAB*PACC
111 SOP=SOP+PARU(2)*B*DELTAB*OV*PACC
112 IF(B.LT.1..OR.B*PACC.GT.1E-6) GOTO 140
113 ENDIF
114 YK=PARU(1)*XK*SO/SP
115
116C...Continue iteration until convergence.
117 IF(YK.LT.YKE) THEN
118 XI=XK
119 YI=YK
120 IF(IIT.EQ.1) IIT=2
121 ELSE
122 XF=XK
123 YF=YK
124 IF(IIT.EQ.0) IIT=1
125 ENDIF
126 IF(ABS(YK-YKE).GE.1E-5*YKE) GOTO 130
127
128C...Store some results for subsequent use.
129 VINT(145)=SIGSUM
130 VINT(146)=SOP/SO
131 VINT(147)=SOP/SP
132
133C...Initialize iteration in xT2 for hardest interaction.
134 ELSEIF(MMUL.EQ.2) THEN
135 IF(MSTP(82).LE.0) THEN
136 ELSEIF(MSTP(82).EQ.1) THEN
137 XT2=1.
138 XT2FAC=XSEC(96,1)/VINT(106)*VINT(149)/(1.-VINT(149))
139 ELSEIF(MSTP(82).EQ.2) THEN
140 XT2=1.
141 XT2FAC=VINT(146)*XSEC(96,1)/VINT(106)*VINT(149)*(1.+VINT(149))
142 ELSE
143 XC2=4.*CKIN(3)**2/VINT(2)
144 IF(CKIN(3).LE.CKIN(5).OR.MINT(82).GE.2) XC2=0.
145 ENDIF
146
147 ELSEIF(MMUL.EQ.3) THEN
148C...Low-pT or multiple interactions (first semihard interaction):
149C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
150C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
151 ISUB=MINT(1)
152 IF(MSTP(82).LE.0) THEN
153 XT2=0.
154 ELSEIF(MSTP(82).EQ.1) THEN
155 XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(RLU_HIJING(0)))
156 ELSEIF(MSTP(82).EQ.2) THEN
157 IF(XT2.LT.1..AND.EXP(-XT2FAC*XT2/(VINT(149)*(XT2+
158 & VINT(149)))).GT.RLU_HIJING(0)) XT2=1.
159 IF(XT2.GE.1.) THEN
160 XT2=(1.+VINT(149))*XT2FAC/(XT2FAC-(1.+VINT(149))*LOG(1.-
161 & RLU_HIJING(0)*(1.-EXP(-XT2FAC/(VINT(149)*(1.+VINT(149)
162 $ ))))))-VINT(149)
163 ELSE
164 XT2=-XT2FAC/LOG(EXP(-XT2FAC/(XT2+VINT(149)))+RLU_HIJING(0)*
165 & (EXP(-XT2FAC/VINT(149))-EXP(-XT2FAC/(XT2+VINT(149)))))-
166 & VINT(149)
167 ENDIF
168 XT2=MAX(0.01*VINT(149),XT2)
169 ELSE
170 XT2=(XC2+VINT(149))*(1.+VINT(149))/(1.+VINT(149)-
171 & RLU_HIJING(0)*(1.-XC2))-VINT(149)
172 XT2=MAX(0.01*VINT(149),XT2)
173 ENDIF
174 VINT(25)=XT2
175
176C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
177 IF(MSTP(82).LE.1.AND.XT2.LT.VINT(149)) THEN
178 IF(MINT(82).EQ.1) NGEN(0,1)=NGEN(0,1)-1
179 IF(MINT(82).EQ.1) NGEN(ISUB,1)=NGEN(ISUB,1)-1
180 ISUB=95
181 MINT(1)=ISUB
182 VINT(21)=0.01*VINT(149)
183 VINT(22)=0.
184 VINT(23)=0.
185 VINT(25)=0.01*VINT(149)
186
187 ELSE
188C...Multiple interactions (first semihard interaction).
189C...Choose tau and y*. Calculate cos(theta-hat).
190 IF(RLU_HIJING(0).LE.COEF(ISUB,1)) THEN
191 TAUP=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU_HIJING(0)
192 TAU=XT2*(1.+TAUP)**2/(4.*TAUP)
193 ELSE
194 TAU=XT2*(1.+TAN(RLU_HIJING(0)*ATAN(SQRT(1./XT2-1.)))**2)
195 ENDIF
196 VINT(21)=TAU
197 CALL PYKLIM_HIJING(2)
198 RYST=RLU_HIJING(0)
199 MYST=1
200 IF(RYST.GT.COEF(ISUB,7)) MYST=2
201 IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3
202 CALL PYKMAP_HIJING(2,MYST,RLU_HIJING(0))
203 VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU_HIJING(0))
204 ENDIF
205 VINT(71)=0.5*VINT(1)*SQRT(VINT(25))
206
207C...Store results of cross-section calculation.
208 ELSEIF(MMUL.EQ.4) THEN
209 ISUB=MINT(1)
210 XTS=VINT(25)
211 IF(ISET(ISUB).EQ.1) XTS=VINT(21)
212 IF(ISET(ISUB).EQ.2) XTS=(4.*VINT(48)+2.*VINT(63)+2.*VINT(64))/
213 & VINT(2)
214 IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) XTS=VINT(26)
215 RBIN=MAX(0.000001,MIN(0.999999,XTS*(1.+VINT(149))/
216 & (XTS+VINT(149))))
217 IRBIN=INT(1.+20.*RBIN)
218 IF(ISUB.EQ.96) NMUL(IRBIN)=NMUL(IRBIN)+1
219 IF(ISUB.EQ.96) SIGM(IRBIN)=SIGM(IRBIN)+VINT(153)
220
221C...Choose impact parameter.
222 ELSEIF(MMUL.EQ.5) THEN
223 IF(MSTP(82).EQ.3) THEN
224 VINT(148)=RLU_HIJING(0)/(PARU(2)*VINT(147))
225 ELSE
226 RTYPE=RLU_HIJING(0)
227 CQ2=PARP(84)**2
228 IF(RTYPE.LT.(1.-PARP(83))**2) THEN
229 B2=-LOG(RLU_HIJING(0))
230 ELSEIF(RTYPE.LT.1.-PARP(83)**2) THEN
231 B2=-0.5*(1.+CQ2)*LOG(RLU_HIJING(0))
232 ELSE
233 B2=-CQ2*LOG(RLU_HIJING(0))
234 ENDIF
235 VINT(148)=((1.-PARP(83))**2*EXP(-MIN(100.,B2))+2.*PARP(83)*
236 & (1.-PARP(83))*2./(1.+CQ2)*EXP(-MIN(100.,B2*2./(1.+CQ2)))+
237 & PARP(83)**2/CQ2*EXP(-MIN(100.,B2/CQ2)))/(PARU(2)*VINT(147))
238 ENDIF
239
240C...Multiple interactions (variable impact parameter) : reject with
241C...probability exp(-overlap*cross-section above pT/normalization).
242 RNCOR=(IRBIN-20.*RBIN)*NMUL(IRBIN)
243 SIGCOR=(IRBIN-20.*RBIN)*SIGM(IRBIN)
244 DO 150 IBIN=IRBIN+1,20
245 RNCOR=RNCOR+NMUL(IBIN)
246 150 SIGCOR=SIGCOR+SIGM(IBIN)
247 SIGABV=(SIGCOR/RNCOR)*VINT(149)*(1.-XTS)/(XTS+VINT(149))
248 VINT(150)=EXP(-MIN(100.,VINT(146)*VINT(148)*SIGABV/VINT(106)))
249
250C...Generate additional multiple semihard interactions.
251 ELSEIF(MMUL.EQ.6) THEN
252
253C...Reconstruct strings in hard scattering.
254 ISUB=MINT(1)
255 NMAX=MINT(84)+4
256 IF(ISET(ISUB).EQ.1) NMAX=MINT(84)+2
257 NSTR=0
258 DO 170 I=MINT(84)+1,NMAX
259 KCS=KCHG(LUCOMP_HIJING(K(I,2)),2)*ISIGN(1,K(I,2))
260 IF(KCS.EQ.0) GOTO 170
261 DO 160 J=1,4
262 IF(KCS.EQ.1.AND.(J.EQ.2.OR.J.EQ.4)) GOTO 160
263 IF(KCS.EQ.-1.AND.(J.EQ.1.OR.J.EQ.3)) GOTO 160
264 IF(J.LE.2) THEN
265 IST=MOD(K(I,J+3)/MSTU(5),MSTU(5))
266 ELSE
267 IST=MOD(K(I,J+1),MSTU(5))
268 ENDIF
269 IF(IST.LT.MINT(84).OR.IST.GT.I) GOTO 160
270 IF(KCHG(LUCOMP_HIJING(K(IST,2)),2).EQ.0) GOTO 160
271 NSTR=NSTR+1
272 IF(J.EQ.1.OR.J.EQ.4) THEN
273 KSTR(NSTR,1)=I
274 KSTR(NSTR,2)=IST
275 ELSE
276 KSTR(NSTR,1)=IST
277 KSTR(NSTR,2)=I
278 ENDIF
279 160 CONTINUE
280 170 CONTINUE
281
282C...Set up starting values for iteration in xT2.
283 XT2=VINT(25)
284 IF(ISET(ISUB).EQ.1) XT2=VINT(21)
285 IF(ISET(ISUB).EQ.2) XT2=(4.*VINT(48)+2.*VINT(63)+2.*VINT(64))/
286 & VINT(2)
287 IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) XT2=VINT(26)
288 ISUB=96
289 MINT(1)=96
290 IF(MSTP(82).LE.1) THEN
291 XT2FAC=XSEC(ISUB,1)*VINT(149)/((1.-VINT(149))*VINT(106))
292 ELSE
293 XT2FAC=VINT(146)*VINT(148)*XSEC(ISUB,1)/VINT(106)*
294 & VINT(149)*(1.+VINT(149))
295 ENDIF
296 VINT(63)=0.
297 VINT(64)=0.
298 VINT(151)=0.
299 VINT(152)=0.
300 VINT(143)=1.-VINT(141)
301 VINT(144)=1.-VINT(142)
302
303C...Iterate downwards in xT2.
304 180 IF(MSTP(82).LE.1) THEN
305 XT2=XT2FAC*XT2/(XT2FAC-XT2*LOG(RLU_HIJING(0)))
306 IF(XT2.LT.VINT(149)) GOTO 220
307 ELSE
308 IF(XT2.LE.0.01*VINT(149)) GOTO 220
309 XT2=XT2FAC*(XT2+VINT(149))/(XT2FAC-(XT2+VINT(149))*
310 & LOG(RLU_HIJING(0)))-VINT(149)
311 IF(XT2.LE.0.) GOTO 220
312 XT2=MAX(0.01*VINT(149),XT2)
313 ENDIF
314 VINT(25)=XT2
315
316C...Choose tau and y*. Calculate cos(theta-hat).
317 IF(RLU_HIJING(0).LE.COEF(ISUB,1)) THEN
318 TAUP=(2.*(1.+SQRT(1.-XT2))/XT2-1.)**RLU_HIJING(0)
319 TAU=XT2*(1.+TAUP)**2/(4.*TAUP)
320 ELSE
321 TAU=XT2*(1.+TAN(RLU_HIJING(0)*ATAN(SQRT(1./XT2-1.)))**2)
322 ENDIF
323 VINT(21)=TAU
324 CALL PYKLIM_HIJING(2)
325 RYST=RLU_HIJING(0)
326 MYST=1
327 IF(RYST.GT.COEF(ISUB,7)) MYST=2
328 IF(RYST.GT.COEF(ISUB,7)+COEF(ISUB,8)) MYST=3
329 CALL PYKMAP_HIJING(2,MYST,RLU_HIJING(0))
330 VINT(23)=SQRT(MAX(0.,1.-XT2/TAU))*(-1)**INT(1.5+RLU_HIJING(0))
331
332C...Check that x not used up. Accept or reject kinematical variables.
333 X1M=SQRT(TAU)*EXP(VINT(22))
334 X2M=SQRT(TAU)*EXP(-VINT(22))
335 IF(VINT(143)-X1M.LT.0.01.OR.VINT(144)-X2M.LT.0.01) GOTO 180
336 VINT(71)=0.5*VINT(1)*SQRT(XT2)
337 CALL PYSIGH_HIJING(NCHN,SIGS)
338 IF(SIGS.LT.XSEC(ISUB,1)*RLU_HIJING(0)) GOTO 180
339
340C...Reset K, P and V vectors. Select some variables.
341 DO 190 I=N+1,N+2
342 DO 190 J=1,5
343 K(I,J)=0
344 P(I,J)=0.
345 190 V(I,J)=0.
346 RFLAV=RLU_HIJING(0)
347 PT=0.5*VINT(1)*SQRT(XT2)
348 PHI=PARU(2)*RLU_HIJING(0)
349 CTH=VINT(23)
350
351C...Add first parton to event record.
352 K(N+1,1)=3
353 K(N+1,2)=21
354 IF(RFLAV.GE.MAX(PARP(85),PARP(86))) K(N+1,2)=
355 & 1+INT((2.+PARJ(2))*RLU_HIJING(0))
356 P(N+1,1)=PT*COS(PHI)
357 P(N+1,2)=PT*SIN(PHI)
358 P(N+1,3)=0.25*VINT(1)*(VINT(41)*(1.+CTH)-VINT(42)*(1.-CTH))
359 P(N+1,4)=0.25*VINT(1)*(VINT(41)*(1.+CTH)+VINT(42)*(1.-CTH))
360 P(N+1,5)=0.
361
362C...Add second parton to event record.
363 K(N+2,1)=3
364 K(N+2,2)=21
365 IF(K(N+1,2).NE.21) K(N+2,2)=-K(N+1,2)
366 P(N+2,1)=-P(N+1,1)
367 P(N+2,2)=-P(N+1,2)
368 P(N+2,3)=0.25*VINT(1)*(VINT(41)*(1.-CTH)-VINT(42)*(1.+CTH))
369 P(N+2,4)=0.25*VINT(1)*(VINT(41)*(1.-CTH)+VINT(42)*(1.+CTH))
370 P(N+2,5)=0.
371
372 IF(RFLAV.LT.PARP(85).AND.NSTR.GE.1) THEN
373C....Choose relevant string pieces to place gluons on.
374 DO 210 I=N+1,N+2
375 DMIN=1E8
376 DO 200 ISTR=1,NSTR
377 I1=KSTR(ISTR,1)
378 I2=KSTR(ISTR,2)
379 DIST=(P(I,4)*P(I1,4)-P(I,1)*P(I1,1)-P(I,2)*P(I1,2)-
380 & P(I,3)*P(I1,3))*(P(I,4)*P(I2,4)-P(I,1)*P(I2,1)-
381 & P(I,2)*P(I2,2)-P(I,3)*P(I2,3))/MAX(1.,P(I1,4)*P(I2,4)-
382 & P(I1,1)*P(I2,1)-P(I1,2)*P(I2,2)-P(I1,3)*P(I2,3))
383 IF(ISTR.EQ.1.OR.DIST.LT.DMIN) THEN
384 DMIN=DIST
385 IST1=I1
386 IST2=I2
387 ISTM=ISTR
388 ENDIF
389 200 CONTINUE
390
391C....Colour flow adjustments, new string pieces.
392 IF(K(IST1,4)/MSTU(5).EQ.IST2) K(IST1,4)=MSTU(5)*I+
393 & MOD(K(IST1,4),MSTU(5))
394 IF(MOD(K(IST1,5),MSTU(5)).EQ.IST2) K(IST1,5)=
395 & MSTU(5)*(K(IST1,5)/MSTU(5))+I
396 K(I,5)=MSTU(5)*IST1
397 K(I,4)=MSTU(5)*IST2
398 IF(K(IST2,5)/MSTU(5).EQ.IST1) K(IST2,5)=MSTU(5)*I+
399 & MOD(K(IST2,5),MSTU(5))
400 IF(MOD(K(IST2,4),MSTU(5)).EQ.IST1) K(IST2,4)=
401 & MSTU(5)*(K(IST2,4)/MSTU(5))+I
402 KSTR(ISTM,2)=I
403 KSTR(NSTR+1,1)=I
404 KSTR(NSTR+1,2)=IST2
405 210 NSTR=NSTR+1
406
407C...String drawing and colour flow for gluon loop.
408 ELSEIF(K(N+1,2).EQ.21) THEN
409 K(N+1,4)=MSTU(5)*(N+2)
410 K(N+1,5)=MSTU(5)*(N+2)
411 K(N+2,4)=MSTU(5)*(N+1)
412 K(N+2,5)=MSTU(5)*(N+1)
413 KSTR(NSTR+1,1)=N+1
414 KSTR(NSTR+1,2)=N+2
415 KSTR(NSTR+2,1)=N+2
416 KSTR(NSTR+2,2)=N+1
417 NSTR=NSTR+2
418
419C...String drawing and colour flow for q-qbar pair.
420 ELSE
421 K(N+1,4)=MSTU(5)*(N+2)
422 K(N+2,5)=MSTU(5)*(N+1)
423 KSTR(NSTR+1,1)=N+1
424 KSTR(NSTR+1,2)=N+2
425 NSTR=NSTR+1
426 ENDIF
427
428C...Update remaining energy; iterate.
429 N=N+2
430 IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
431 CALL LUERRM_HIJING(11
432 $ ,'(PYMULT_HIJING:) no more memory left in LUJETS_HIJING'
433 $ )
434 IF(MSTU(21).GE.1) RETURN
435 ENDIF
436 MINT(31)=MINT(31)+1
437 VINT(151)=VINT(151)+VINT(41)
438 VINT(152)=VINT(152)+VINT(42)
439 VINT(143)=VINT(143)-VINT(41)
440 VINT(144)=VINT(144)-VINT(42)
441 IF(MINT(31).LT.240) GOTO 180
442 220 CONTINUE
443 ENDIF
444
445C...Format statements for printout.
446 1000 FORMAT(/1X
447 $ ,'****** PYMULT_HIJING: initialization of multiple inter'
448 $ ,'actions for MSTP(82) =',I2,' ******')
449 1100 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
450 &E9.2,' mb: rejected')
451 1200 FORMAT(8X,'pT0 =',F5.2,' GeV gives sigma(parton-parton) =',1P,
452 &E9.2,' mb: accepted')
453
454 RETURN
455 END