Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pythia_hijing.F
CommitLineData
e74335a4 1* $Id$
2
3C*********************************************************************
4
5 SUBROUTINE PYTHIA_HIJING
6
7C...Administers the generation of a high-pt event via calls to a number
8C...of subroutines; also computes cross-sections.
9#include "lujets_hijing.inc"
10#include "ludat1_hijing.inc"
11#include "ludat2_hijing.inc"
12#include "pysubs_hijing.inc"
13#include "pypars_hijing.inc"
14#include "pyint1_hijing.inc"
15#include "pyint2_hijing.inc"
16#include "pyint5_hijing.inc"
17
18C...Loop over desired number of overlayed events (normally 1).
19 MINT(7)=0
20 MINT(8)=0
21 NOVL=1
22 IF(MSTP(131).NE.0) CALL PYOVLY_HIJING(2)
23 IF(MSTP(131).NE.0) NOVL=MINT(81)
24 MINT(83)=0
25 MINT(84)=MSTP(126)
26 MSTU(70)=0
27 DO 190 IOVL=1,NOVL
28 IF(MINT(84)+100.GE.MSTU(4)) THEN
29 CALL LUERRM_HIJING(11,
30 & '(PYTHIA_HIJING:) no more space in LUJETS_HIJING '/
31 $ /'for overlayed events')
32 IF(MSTU(21).GE.1) GOTO 200
33 ENDIF
34 MINT(82)=IOVL
35
36C...Generate variables of hard scattering.
37 100 CONTINUE
38 IF(IOVL.EQ.1) NGEN(0,2)=NGEN(0,2)+1
39 MINT(31)=0
40 MINT(51)=0
41 CALL PYRAND_HIJING
42 ISUB=MINT(1)
43 IF(IOVL.EQ.1) THEN
44 NGEN(ISUB,2)=NGEN(ISUB,2)+1
45
46C...Store information on hard interaction.
47 DO 110 J=1,200
48 MSTI(J)=0
49 110 PARI(J)=0.
50 MSTI(1)=MINT(1)
51 MSTI(2)=MINT(2)
52 MSTI(11)=MINT(11)
53 MSTI(12)=MINT(12)
54 MSTI(15)=MINT(15)
55 MSTI(16)=MINT(16)
56 MSTI(17)=MINT(17)
57 MSTI(18)=MINT(18)
58 PARI(11)=VINT(1)
59 PARI(12)=VINT(2)
60 IF(ISUB.NE.95) THEN
61 DO 120 J=13,22
62 120 PARI(J)=VINT(30+J)
63 PARI(33)=VINT(41)
64 PARI(34)=VINT(42)
65 PARI(35)=PARI(33)-PARI(34)
66 PARI(36)=VINT(21)
67 PARI(37)=VINT(22)
68 PARI(38)=VINT(26)
69 PARI(41)=VINT(23)
70 ENDIF
71 ENDIF
72
73 IF(MSTP(111).EQ.-1) GOTO 160
74 IF(ISUB.LE.90.OR.ISUB.GE.95) THEN
75C...Hard scattering (including low-pT):
76C...reconstruct kinematics and colour flow of hard scattering.
77 CALL PYSCAT_HIJING
78 IF(MINT(51).EQ.1) GOTO 100
79
80C...Showering of initial state partons (optional).
81 IPU1=MINT(84)+1
82 IPU2=MINT(84)+2
83 IF(MSTP(61).GE.1.AND.MINT(43).NE.1.AND.ISUB.NE.95)
84 & CALL PYSSPA_HIJING(IPU1,IPU2)
85 NSAV1=N
86
87C...Multiple interactions.
88 IF(MSTP(81).GE.1.AND.MINT(43).EQ.4.AND.ISUB.NE.95)
89 & CALL PYMULT_HIJING(6)
90 MINT(1)=ISUB
91 NSAV2=N
92
93C...Hadron remnants and primordial kT.
94 CALL PYREMN_HIJING(IPU1,IPU2)
95 IF(MINT(51).EQ.1) GOTO 100
96 NSAV3=N
97
98C...Showering of final state partons (optional).
99 IPU3=MINT(84)+3
100 IPU4=MINT(84)+4
101 IF(MSTP(71).GE.1.AND.ISUB.NE.95.AND.K(IPU3,1).GT.0.AND.
102 & K(IPU3,1).LE.10.AND.K(IPU4,1).GT.0.AND.K(IPU4,1).LE.10) THEN
103 QMAX=SQRT(PARP(71)*VINT(52))
104 IF(ISUB.EQ.5) QMAX=SQRT(PMAS(23,1)**2)
105 IF(ISUB.EQ.8) QMAX=SQRT(PMAS(24,1)**2)
106 CALL LUSHOW_HIJING(IPU3,IPU4,QMAX)
107 ENDIF
108
109C...Sum up transverse and longitudinal momenta.
110 IF(IOVL.EQ.1) THEN
111 PARI(65)=2.*PARI(17)
112 DO 130 I=MSTP(126)+1,N
113 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 130
114 PT=SQRT(P(I,1)**2+P(I,2)**2)
115 PARI(69)=PARI(69)+PT
116 IF(I.LE.NSAV1.OR.I.GT.NSAV3) PARI(66)=PARI(66)+PT
117 IF(I.GT.NSAV1.AND.I.LE.NSAV2) PARI(68)=PARI(68)+PT
118 130 CONTINUE
119 PARI(67)=PARI(68)
120 PARI(71)=VINT(151)
121 PARI(72)=VINT(152)
122 PARI(73)=VINT(151)
123 PARI(74)=VINT(152)
124 ENDIF
125
126C...Decay of final state resonances.
127 IF(MSTP(41).GE.1.AND.ISUB.NE.95) CALL PYRESD_HIJING
128
129 ELSE
130C...Diffractive and elastic scattering.
131 CALL PYDIFF_HIJING
132 IF(IOVL.EQ.1) THEN
133 PARI(65)=2.*PARI(17)
134 PARI(66)=PARI(65)
135 PARI(69)=PARI(65)
136 ENDIF
137 ENDIF
138
139C...Recalculate energies from momenta and masses (if desired).
140 IF(MSTP(113).GE.1) THEN
141 DO 140 I=MINT(83)+1,N
142 140 IF(K(I,1).GT.0.AND.K(I,1).LE.10) P(I,4)=SQRT(P(I,1)**2+
143 & P(I,2)**2+P(I,3)**2+P(I,5)**2)
144 ENDIF
145
146C...Rearrange partons along strings, check invariant mass cuts.
147 MSTU(28)=0
148 CALL LUPREP_HIJING(MINT(84)+1)
149 IF(MSTP(112).EQ.1.AND.MSTU(28).EQ.3) GOTO 100
150 IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) THEN
151 DO 150 I=MINT(84)+1,N
152 IF(K(I,2).NE.94) GOTO 150
153 K(I+1,3)=MOD(K(I+1,4)/MSTU(5),MSTU(5))
154 K(I+2,3)=MOD(K(I+2,4)/MSTU(5),MSTU(5))
155 150 CONTINUE
156 CALL LUEDIT_HIJING(12)
157 CALL LUEDIT_HIJING(14)
158 IF(MSTP(125).EQ.0) CALL LUEDIT_HIJING(15)
159 IF(MSTP(125).EQ.0) MINT(4)=0
160 ENDIF
161
162C...Introduce separators between sections in LULIST_HIJING event listing.
163 IF(IOVL.EQ.1.AND.MSTP(125).LE.0) THEN
164 MSTU(70)=1
165 MSTU(71)=N
166 ELSEIF(IOVL.EQ.1) THEN
167 MSTU(70)=3
168 MSTU(71)=2
169 MSTU(72)=MINT(4)
170 MSTU(73)=N
171 ENDIF
172
173C...Perform hadronization (if desired).
174 IF(MSTP(111).GE.1) CALL LUEXEC_HIJING
175 IF(MSTP(125).EQ.0.OR.MSTP(125).EQ.1) CALL LUEDIT_HIJING(14)
176
177C...Calculate Monte Carlo estimates of cross-sections.
178 160 IF(IOVL.EQ.1) THEN
179 IF(MSTP(111).NE.-1) NGEN(ISUB,3)=NGEN(ISUB,3)+1
180 NGEN(0,3)=NGEN(0,3)+1
181 XSEC(0,3)=0.
182 DO 170 I=1,200
183 IF(I.EQ.96) THEN
184 XSEC(I,3)=0.
185 ELSEIF(MSUB(95).EQ.1.AND.(I.EQ.11.OR.I.EQ.12.OR.I.EQ.13.OR.
186 & I.EQ.28.OR.I.EQ.53.OR.I.EQ.68)) THEN
187 XSEC(I,3)=XSEC(96,2)*NGEN(I,3)/MAX(1.,FLOAT(NGEN(96,1))*
188 & FLOAT(NGEN(96,2)))
189 ELSEIF(NGEN(I,1).EQ.0) THEN
190 XSEC(I,3)=0.
191 ELSEIF(NGEN(I,2).EQ.0) THEN
192 XSEC(I,3)=XSEC(I,2)*NGEN(0,3)/(FLOAT(NGEN(I,1))*
193 & FLOAT(NGEN(0,2)))
194 ELSE
195 XSEC(I,3)=XSEC(I,2)*NGEN(I,3)/(FLOAT(NGEN(I,1))*
196 & FLOAT(NGEN(I,2)))
197 ENDIF
198 170 XSEC(0,3)=XSEC(0,3)+XSEC(I,3)
199 IF(MSUB(95).EQ.1) THEN
200 NGENS=NGEN(91,3)+NGEN(92,3)+NGEN(93,3)+NGEN(94,3)+NGEN(95,3)
201 XSECS=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+XSEC(95,3)
202 XMAXS=XSEC(95,1)
203 IF(MSUB(91).EQ.1) XMAXS=XMAXS+XSEC(91,1)
204 IF(MSUB(92).EQ.1) XMAXS=XMAXS+XSEC(92,1)
205 IF(MSUB(93).EQ.1) XMAXS=XMAXS+XSEC(93,1)
206 IF(MSUB(94).EQ.1) XMAXS=XMAXS+XSEC(94,1)
207 FAC=1.
208 IF(NGENS.LT.NGEN(0,3)) FAC=(XMAXS-XSECS)/(XSEC(0,3)-XSECS)
209 XSEC(11,3)=FAC*XSEC(11,3)
210 XSEC(12,3)=FAC*XSEC(12,3)
211 XSEC(13,3)=FAC*XSEC(13,3)
212 XSEC(28,3)=FAC*XSEC(28,3)
213 XSEC(53,3)=FAC*XSEC(53,3)
214 XSEC(68,3)=FAC*XSEC(68,3)
215 XSEC(0,3)=XSEC(91,3)+XSEC(92,3)+XSEC(93,3)+XSEC(94,3)+
216 & XSEC(95,1)
217 ENDIF
218
219C...Store final information.
220 MINT(5)=MINT(5)+1
221 MSTI(3)=MINT(3)
222 MSTI(4)=MINT(4)
223 MSTI(5)=MINT(5)
224 MSTI(6)=MINT(6)
225 MSTI(7)=MINT(7)
226 MSTI(8)=MINT(8)
227 MSTI(13)=MINT(13)
228 MSTI(14)=MINT(14)
229 MSTI(21)=MINT(21)
230 MSTI(22)=MINT(22)
231 MSTI(23)=MINT(23)
232 MSTI(24)=MINT(24)
233 MSTI(25)=MINT(25)
234 MSTI(26)=MINT(26)
235 MSTI(31)=MINT(31)
236 PARI(1)=XSEC(0,3)
237 PARI(2)=XSEC(0,3)/MINT(5)
238 PARI(31)=VINT(141)
239 PARI(32)=VINT(142)
240 IF(ISUB.NE.95.AND.MINT(7)*MINT(8).NE.0) THEN
241 PARI(42)=2.*VINT(47)/VINT(1)
242 DO 180 IS=7,8
243 PARI(36+IS)=P(MINT(IS),3)/VINT(1)
244 PARI(38+IS)=P(MINT(IS),4)/VINT(1)
245 I=MINT(IS)
246 PR=MAX(1E-20,P(I,5)**2+P(I,1)**2+P(I,2)**2)
247 PARI(40+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
248 & SQRT(PR),1E20)),P(I,3))
249 PR=MAX(1E-20,P(I,1)**2+P(I,2)**2)
250 PARI(42+IS)=SIGN(LOG(MIN((SQRT(PR+P(I,3)**2)+ABS(P(I,3)))/
251 & SQRT(PR),1E20)),P(I,3))
252 PARI(44+IS)=P(I,3)/SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
253 PARI(46+IS)=ULANGL_HIJING(P(I,3),SQRT(P(I,1)**2+P(I,2)**2))
254 PARI(48+IS)=ULANGL_HIJING(P(I,1),P(I,2))
255 180 CONTINUE
256 ENDIF
257 PARI(61)=VINT(148)
258 IF(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.3) THEN
259 MSTU(161)=MINT(21)
260 MSTU(162)=0
261 ELSE
262 MSTU(161)=MINT(21)
263 MSTU(162)=MINT(22)
264 ENDIF
265 ENDIF
266
267C...Prepare to go to next overlayed event.
268 MSTI(41)=IOVL
269 IF(IOVL.GE.2.AND.IOVL.LE.10) MSTI(40+IOVL)=ISUB
270 IF(MSTU(70).LT.10) THEN
271 MSTU(70)=MSTU(70)+1
272 MSTU(70+MSTU(70))=N
273 ENDIF
274 MINT(83)=N
275 MINT(84)=N+MSTP(126)
276 190 CONTINUE
277
278C...Information on overlayed events.
279 IF(MSTP(131).EQ.1.AND.MSTP(133).GE.1) THEN
280 PARI(91)=VINT(132)
281 PARI(92)=VINT(133)
282 PARI(93)=VINT(134)
283 IF(MSTP(133).EQ.2) PARI(93)=PARI(93)*XSEC(0,3)/VINT(131)
284 ENDIF
285
286C...Transform to the desired coordinate frame.
287 200 CALL PYFRAM_HIJING(MSTP(124))
288
289 RETURN
290 END