Merging the VirtualMC branch to the main development branch (HEAD)
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pymult_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE PYMULT_HIJING(MMUL)   
6     
7 C...Initializes treatment of multiple interactions, selects kinematics  
8 C...of hardest interaction if low-pT physics included in run, and   
9 C...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     
22 C...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     
32 C...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     
43 C...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     
59 C...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     
66 C...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     
75 C...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     
90 C...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     
116 C...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     
128 C...Store some results for subsequent use.  
129         VINT(145)=SIGSUM    
130         VINT(146)=SOP/SO    
131         VINT(147)=SOP/SP    
132     
133 C...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    
148 C...Low-pT or multiple interactions (first semihard interaction):   
149 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)    
150 C...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     
176 C...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    
188 C...Multiple interactions (first semihard interaction). 
189 C...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     
207 C...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     
221 C...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     
240 C...Multiple interactions (variable impact parameter) : reject with 
241 C...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     
250 C...Generate additional multiple semihard interactions. 
251       ELSEIF(MMUL.EQ.6) THEN    
252     
253 C...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     
282 C...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     
303 C...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     
316 C...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     
332 C...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     
340 C...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     
351 C...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     
362 C...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    
373 C....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     
391 C....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     
407 C...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     
419 C...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     
428 C...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     
445 C...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