Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymult.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymult.f
1 
2 C*********************************************************************
3 
4 C...PYMULT
5 C...Initializes treatment of multiple interactions, selects kinematics
6 C...of hardest interaction if low-pT physics included in run, and
7 C...generates all non-hardest interactions.
8 
9  SUBROUTINE pymult(MMUL)
10 
11 C...Double precision and integer declarations.
12  IMPLICIT DOUBLE PRECISION(a-h, o-z)
13  IMPLICIT INTEGER(i-n)
14  INTEGER pyk,pychge,pycomp
15 C...Commonblocks.
16  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
17  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
18  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
19  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
20  common/pypars/mstp(200),parp(200),msti(200),pari(200)
21  common/pyint1/mint(400),vint(400)
22  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
23  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
24  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
25  common/pyint7/sigt(0:6,0:6,0:5)
26  SAVE /pyjets/,/pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/,
27  &/pyint2/,/pyint3/,/pyint5/,/pyint7/
28 C...Local arrays and saved variables.
29  dimension nmul(20),sigm(20),kstr(500,2),vintsv(80)
30  SAVE xt2,xt2fac,xc2,xts,irbin,rbin,nmul,sigm,p83a,p83b,p83c,
31  &cq2i,cq2r,pik,bdiv,b,plowb,phighb,pallb,s4a,s4b,s4c,powip,
32  &rpwip,b2rpdv,b2rpmx,bavg,vnt145,vnt146,vnt147
33 
34 C...Initialization of multiple interaction treatment.
35  IF(mmul.EQ.1) THEN
36  IF(mstp(122).GE.1) WRITE(mstu(11),5000) mstp(82)
37  isub=96
38  mint(1)=96
39  vint(63)=0d0
40  vint(64)=0d0
41  vint(143)=1d0
42  vint(144)=1d0
43 
44 C...Loop over phase space points: xT2 choice in 20 bins.
45  100 sigsum=0d0
46  DO 120 ixt2=1,20
47  nmul(ixt2)=mstp(83)
48  sigm(ixt2)=0d0
49  DO 110 itry=1,mstp(83)
50  rsca=0.05d0*((21-ixt2)-pyr(0))
51  xt2=vint(149)*(1d0+vint(149))/(vint(149)+rsca)-vint(149)
52  xt2=max(0.01d0*vint(149),xt2)
53  vint(25)=xt2
54 
55 C...Choose tau and y*. Calculate cos(theta-hat).
56  IF(pyr(0).LE.coef(isub,1)) THEN
57  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
58  tau=xt2*(1d0+taut)**2/(4d0*taut)
59  ELSE
60  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
61  ENDIF
62  vint(21)=tau
63  CALL pyklim(2)
64  ryst=pyr(0)
65  myst=1
66  IF(ryst.GT.coef(isub,8)) myst=2
67  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
68  CALL pykmap(2,myst,pyr(0))
69  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
70 
71 C...Calculate differential cross-section.
72  vint(71)=0.5d0*vint(1)*sqrt(xt2)
73  CALL pysigh(nchn,sigs)
74  sigm(ixt2)=sigm(ixt2)+sigs
75  110 CONTINUE
76  sigsum=sigsum+sigm(ixt2)
77  120 CONTINUE
78  sigsum=sigsum/(20d0*mstp(83))
79 
80 C...Reject result if sigma(parton-parton) is smaller than hadronic one.
81  IF(sigsum.LT.1.1d0*sigt(0,0,5)) THEN
82  IF(mstp(122).GE.1) WRITE(mstu(11),5100)
83  & parp(82)*(vint(1)/parp(89))**parp(90),sigsum
84  parp(82)=0.9d0*parp(82)
85  vint(149)=4d0*(parp(82)*(vint(1)/parp(89))**parp(90))**2/
86  & vint(2)
87  goto 100
88  ENDIF
89  IF(mstp(122).GE.1) WRITE(mstu(11),5200)
90  & parp(82)*(vint(1)/parp(89))**parp(90), sigsum
91 
92 C...Start iteration to find k factor.
93  yke=sigsum/max(1d-10,sigt(0,0,5))
94  p83a=(1d0-parp(83))**2
95  p83b=2d0*parp(83)*(1d0-parp(83))
96  p83c=parp(83)**2
97  cq2i=1d0/parp(84)**2
98  cq2r=2d0/(1d0+parp(84)**2)
99  so=0.5d0
100  xi=0d0
101  yi=0d0
102  xf=0d0
103  yf=0d0
104  xk=0.5d0
105  iit=0
106  130 IF(iit.EQ.0) THEN
107  xk=2d0*xk
108  ELSEIF(iit.EQ.1) THEN
109  xk=0.5d0*xk
110  ELSE
111  xk=xi+(yke-yi)*(xf-xi)/(yf-yi)
112  ENDIF
113 
114 C...Evaluate overlap integrals. Find where to divide the b range.
115  IF(mstp(82).EQ.2) THEN
116  sp=0.5d0*paru(1)*(1d0-exp(-xk))
117  sop=sp/paru(1)
118  ELSE
119  IF(mstp(82).EQ.3) THEN
120  deltab=0.02d0
121  ELSEIF(mstp(82).EQ.4) THEN
122  deltab=min(0.01d0,0.05d0*parp(84))
123  ELSE
124  powip=max(0.4d0,parp(83))
125  rpwip=2d0/powip-1d0
126  deltab=max(0.02d0,0.02d0*(2d0/powip)**(1d0/powip))
127  so=0d0
128  ENDIF
129  sp=0d0
130  sop=0d0
131  bsp=0d0
132  sohigh=0d0
133  ibdiv=0
134  b=-0.5d0*deltab
135  140 b=b+deltab
136  IF(mstp(82).EQ.3) THEN
137  ov=exp(-b**2)/paru(2)
138  ELSEIF(mstp(82).EQ.4) THEN
139  ov=(p83a*exp(-min(50d0,b**2))+
140  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
141  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
142  ELSE
143  ov=exp(-b**powip)/paru(2)
144  so=so+paru(2)*b*deltab*ov
145  ENDIF
146  IF(ibdiv.EQ.1) sohigh=sohigh+paru(2)*b*deltab*ov
147  pacc=1d0-exp(-min(50d0,paru(1)*xk*ov))
148  sp=sp+paru(2)*b*deltab*pacc
149  sop=sop+paru(2)*b*deltab*ov*pacc
150  bsp=bsp+b*paru(2)*b*deltab*pacc
151  IF(ibdiv.EQ.0.AND.paru(1)*xk*ov.LT.1d0) THEN
152  ibdiv=1
153  bdiv=b+0.5d0*deltab
154  ENDIF
155  IF(b.LT.1d0.OR.b*pacc.GT.1d-6) goto 140
156  ENDIF
157  yk=paru(1)*xk*so/sp
158 
159 C...Continue iteration until convergence.
160  IF(yk.LT.yke) THEN
161  xi=xk
162  yi=yk
163  IF(iit.EQ.1) iit=2
164  ELSE
165  xf=xk
166  yf=yk
167  IF(iit.EQ.0) iit=1
168  ENDIF
169  IF(abs(yk-yke).GE.1d-5*yke) goto 130
170 
171 C...Store some results for subsequent use.
172  bavg=bsp/sp
173  vint(145)=sigsum
174  vint(146)=sop/so
175  vint(147)=sop/sp
176  vnt145=vint(145)
177  vnt146=vint(146)
178  vnt147=vint(147)
179 C...PIK = PARU(1)*XK = (VINT(146)/VINT(147))*sigma_jet/sigma_nondiffr.
180  pik=(vnt146/vnt147)*yke
181 
182 C...Find relative weight for low and high impact parameter.
183  plowb=paru(1)*bdiv**2
184  IF(mstp(82).EQ.3) THEN
185  phighb=pik*0.5*exp(-bdiv**2)
186  ELSEIF(mstp(82).EQ.4) THEN
187  s4a=p83a*exp(-bdiv**2)
188  s4b=p83b*exp(-bdiv**2*cq2r)
189  s4c=p83c*exp(-bdiv**2*cq2i)
190  phighb=pik*0.5*(s4a+s4b+s4c)
191  ELSEIF(parp(83).GE.1.999d0) THEN
192  phighb=pik*sohigh
193  b2rpdv=bdiv**powip
194  ELSE
195  phighb=pik*sohigh
196  b2rpdv=bdiv**powip
197  b2rpmx=max(2d0*rpwip,b2rpdv)
198  ENDIF
199  pallb=plowb+phighb
200 
201 C...Initialize iteration in xT2 for hardest interaction.
202  ELSEIF(mmul.EQ.2) THEN
203  vint(145)=vnt145
204  vint(146)=vnt146
205  vint(147)=vnt147
206  IF(mstp(82).LE.0) THEN
207  ELSEIF(mstp(82).EQ.1) THEN
208  xt2=1d0
209  sigrat=xsec(96,1)/max(1d-10,vint(315)*vint(316)*sigt(0,0,5))
210  IF(mint(141).NE.0.OR.mint(142).NE.0) sigrat=sigrat*
211  & vint(317)/(vint(318)*vint(320))
212  xt2fac=sigrat*vint(149)/(1d0-vint(149))
213  ELSEIF(mstp(82).EQ.2) THEN
214  xt2=1d0
215  xt2fac=vnt146*xsec(96,1)/max(1d-10,sigt(0,0,5))*
216  & vint(149)*(1d0+vint(149))
217  ELSE
218  xc2=4d0*ckin(3)**2/vint(2)
219  IF(ckin(3).LE.ckin(5).OR.mint(82).GE.2) xc2=0d0
220  ENDIF
221 
222 C...Select impact parameter for hardest interaction.
223  IF(mstp(82).LE.2) RETURN
224  142 IF(pyr(0)*pallb.LT.plowb) THEN
225 C...Treatment in low b region.
226  mint(39)=1
227  b=bdiv*sqrt(pyr(0))
228  IF(mstp(82).EQ.3) THEN
229  ov=exp(-b**2)/paru(2)
230  ELSEIF(mstp(82).EQ.4) THEN
231  ov=(p83a*exp(-min(50d0,b**2))+
232  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
233  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
234  ELSE
235  ov=exp(-b**powip)/paru(2)
236  ENDIF
237  vint(148)=ov/vnt147
238  pacc=1d0-exp(-min(50d0,pik*ov))
239  xt2=1d0
240  xt2fac=vnt146*vint(148)*xsec(96,1)/max(1d-10,sigt(0,0,5))*
241  & vint(149)*(1d0+vint(149))
242  ELSE
243 C...Treatment in high b region.
244  mint(39)=2
245  IF(mstp(82).EQ.3) THEN
246  b=sqrt(bdiv**2-log(pyr(0)))
247  ov=exp(-b**2)/paru(2)
248  ELSEIF(mstp(82).EQ.4) THEN
249  s4rndm=pyr(0)*(s4a+s4b+s4c)
250  IF(s4rndm.LT.s4a) THEN
251  b=sqrt(bdiv**2-log(pyr(0)))
252  ELSEIF(s4rndm.LT.s4a+s4b) THEN
253  b=sqrt(bdiv**2-log(pyr(0))/cq2r)
254  ELSE
255  b=sqrt(bdiv**2-log(pyr(0))/cq2i)
256  ENDIF
257  ov=(p83a*exp(-min(50d0,b**2))+
258  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
259  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
260  ELSEIF(parp(83).GE.1.999d0) THEN
261  144 b2rpw=b2rpdv-log(pyr(0))
262  accip=(b2rpw/b2rpdv)**rpwip
263  IF(accip.LT.pyr(0)) goto 144
264  ov=exp(-b2rpw)/paru(2)
265  b=b2rpw**(1d0/powip)
266  ELSE
267  146 b2rpw=b2rpdv-2d0*log(pyr(0))
268  accip=(b2rpw/b2rpmx)**rpwip*exp(-0.5d0*(b2rpw-b2rpmx))
269  IF(accip.LT.pyr(0)) goto 146
270  ov=exp(-b2rpw)/paru(2)
271  b=b2rpw**(1d0/powip)
272  ENDIF
273  vint(148)=ov/vnt147
274  pacc=(1d0-exp(-min(50d0,pik*ov)))/(pik*ov)
275  ENDIF
276  IF(pacc.LT.pyr(0)) goto 142
277  vint(139)=b/bavg
278 
279  ELSEIF(mmul.EQ.3) THEN
280 C...Low-pT or multiple interactions (first semihard interaction):
281 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
282 C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
283  isub=mint(1)
284  vint(145)=vnt145
285  vint(146)=vnt146
286  vint(147)=vnt147
287  IF(mstp(82).LE.0) THEN
288  xt2=0d0
289  ELSEIF(mstp(82).EQ.1) THEN
290  xt2=xt2fac*xt2/(xt2fac-xt2*log(pyr(0)))
291 C...Use with "Sudakov" for low b values when impact parameter dependence.
292  ELSEIF(mstp(82).EQ.2.OR.mint(39).EQ.1) THEN
293  IF(xt2.LT.1d0.AND.exp(-xt2fac*xt2/(vint(149)*(xt2+
294  & vint(149)))).GT.pyr(0)) xt2=1d0
295  IF(xt2.GE.1d0) THEN
296  xt2=(1d0+vint(149))*xt2fac/(xt2fac-(1d0+vint(149))*log(1d0-
297  & pyr(0)*(1d0-exp(-xt2fac/(vint(149)*(1d0+vint(149)))))))-
298  & vint(149)
299  ELSE
300  xt2=-xt2fac/log(exp(-xt2fac/(xt2+vint(149)))+pyr(0)*
301  & (exp(-xt2fac/vint(149))-exp(-xt2fac/(xt2+vint(149)))))-
302  & vint(149)
303  ENDIF
304  xt2=max(0.01d0*vint(149),xt2)
305 C...Use without "Sudakov" for high b values when impact parameter dep.
306  ELSE
307  xt2=(xc2+vint(149))*(1d0+vint(149))/(1d0+vint(149)-
308  & pyr(0)*(1d0-xc2))-vint(149)
309  xt2=max(0.01d0*vint(149),xt2)
310  ENDIF
311  vint(25)=xt2
312 
313 C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
314  IF(mstp(82).LE.1.AND.xt2.LT.vint(149)) THEN
315  IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)-mint(143)
316  IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)-mint(143)
317  isub=95
318  mint(1)=isub
319  vint(21)=0.01d0*vint(149)
320  vint(22)=0d0
321  vint(23)=0d0
322  vint(25)=0.01d0*vint(149)
323 
324  ELSE
325 C...Multiple interactions (first semihard interaction).
326 C...Choose tau and y*. Calculate cos(theta-hat).
327  IF(pyr(0).LE.coef(isub,1)) THEN
328  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
329  tau=xt2*(1d0+taut)**2/(4d0*taut)
330  ELSE
331  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
332  ENDIF
333  vint(21)=tau
334  CALL pyklim(2)
335  ryst=pyr(0)
336  myst=1
337  IF(ryst.GT.coef(isub,8)) myst=2
338  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
339  CALL pykmap(2,myst,pyr(0))
340  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
341  ENDIF
342  vint(71)=0.5d0*vint(1)*sqrt(vint(25))
343 
344 C...Store results of cross-section calculation.
345  ELSEIF(mmul.EQ.4) THEN
346  isub=mint(1)
347  vint(145)=vnt145
348  vint(146)=vnt146
349  vint(147)=vnt147
350  xts=vint(25)
351  IF(iset(isub).EQ.1) xts=vint(21)
352  IF(iset(isub).EQ.2)
353  & xts=(4d0*vint(48)+2d0*vint(63)+2d0*vint(64))/vint(2)
354  IF(iset(isub).GE.3.AND.iset(isub).LE.5) xts=vint(26)
355  rbin=max(0.000001d0,min(0.999999d0,xts*(1d0+vint(149))/
356  & (xts+vint(149))))
357  irbin=int(1d0+20d0*rbin)
358  IF(isub.EQ.96.AND.mstp(171).EQ.0) THEN
359  nmul(irbin)=nmul(irbin)+1
360  sigm(irbin)=sigm(irbin)+vint(153)
361  ENDIF
362 
363 C...Choose impact parameter if not already done.
364  ELSEIF(mmul.EQ.5) THEN
365  isub=mint(1)
366  vint(145)=vnt145
367  vint(146)=vnt146
368  vint(147)=vnt147
369  150 IF(mint(39).GT.0) THEN
370  ELSEIF(mstp(82).EQ.3) THEN
371  expb2=pyr(0)
372  b2=-log(pyr(0))
373  vint(148)=expb2/(paru(2)*vnt147)
374  vint(139)=sqrt(b2)/bavg
375  ELSEIF(mstp(82).EQ.4) THEN
376  rtype=pyr(0)
377  IF(rtype.LT.p83a) THEN
378  b2=-log(pyr(0))
379  ELSEIF(rtype.LT.p83a+p83b) THEN
380  b2=-log(pyr(0))/cq2r
381  ELSE
382  b2=-log(pyr(0))/cq2i
383  ENDIF
384  vint(148)=(p83a*exp(-min(50d0,b2))+
385  & p83b*cq2r*exp(-min(50d0,b2*cq2r))+
386  & p83c*cq2i*exp(-min(50d0,b2*cq2i)))/(paru(2)*vnt147)
387  vint(139)=sqrt(b2)/bavg
388  ELSEIF(parp(83).GE.1.999d0) THEN
389  powip=max(2d0,parp(83))
390  rpwip=2d0/powip-1d0
391  prob1=powip/(2d0*exp(-1d0)+powip)
392  160 IF(pyr(0).LT.prob1) THEN
393  b2rpw=pyr(0)**(0.5d0*powip)
394  accip=exp(-b2rpw)
395  ELSE
396  b2rpw=1d0-log(pyr(0))
397  accip=b2rpw**rpwip
398  ENDIF
399  IF(accip.LT.pyr(0)) goto 160
400  vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
401  vint(139)=b2rpw**(1d0/powip)/bavg
402  ELSE
403  powip=max(0.4d0,parp(83))
404  rpwip=2d0/powip-1d0
405  prob1=rpwip/(rpwip+2d0**rpwip*exp(-rpwip))
406  170 IF(pyr(0).LT.prob1) THEN
407  b2rpw=2d0*rpwip*pyr(0)
408  accip=(b2rpw/rpwip)**rpwip*exp(rpwip-b2rpw)
409  ELSE
410  b2rpw=2d0*(rpwip-log(pyr(0)))
411  accip=(0.5d0*b2rpw/rpwip)**rpwip*exp(rpwip-0.5d0*b2rpw)
412  ENDIF
413  IF(accip.lt .pyr(0)) goto 170
414  vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
415  vint(139)=b2rpw**(1d0/powip)/bavg
416  ENDIF
417 
418 C...Multiple interactions (variable impact parameter) : reject with
419 C...probability exp(-overlap*cross-section above pT/normalization).
420 C...Does not apply to low-b region, where "Sudakov" already included.
421  vint(150)=1d0
422  IF(mint(39).NE.1) THEN
423  rncor=(irbin-20d0*rbin)*nmul(irbin)
424  sigcor=(irbin-20d0*rbin)*sigm(irbin)
425  DO 180 ibin=irbin+1,20
426  rncor=rncor+nmul(ibin)
427  sigcor=sigcor+sigm(ibin)
428  180 CONTINUE
429  sigabv=(sigcor/rncor)*vint(149)*(1d0-xts)/(xts+vint(149))
430  IF(mstp(171).EQ.1) sigabv=sigabv*vint(2)/vint(289)
431  vint(150)=exp(-min(50d0,vnt146*vint(148)*
432  & sigabv/max(1d-10,sigt(0,0,5))))
433  ENDIF
434  IF(mstp(86).EQ.3.OR.(mstp(86).EQ.2.AND.isub.NE.11.AND.
435  & isub.NE.12.AND.isub.NE.13.AND.isub.NE.28.AND.isub.NE.53
436  & .AND.isub.NE.68.AND.isub.NE.95.AND.isub.NE.96)) THEN
437  IF(vint(150).LT.pyr(0)) goto 150
438  vint(150)=1d0
439  ENDIF
440 
441 C...Generate additional multiple semihard interactions.
442  ELSEIF(mmul.EQ.6) THEN
443  isubsv=mint(1)
444  vint(145)=vnt145
445  vint(146)=vnt146
446  vint(147)=vnt147
447  DO 190 j=11,80
448  vintsv(j)=vint(j)
449  190 CONTINUE
450  isub=96
451  mint(1)=96
452  vint(151)=0d0
453  vint(152)=0d0
454 
455 C...Reconstruct strings in hard scattering.
456  nmax=mint(84)+4
457  IF(iset(isubsv).EQ.1) nmax=mint(84)+2
458  IF(iset(isubsv).EQ.11) nmax=mint(84)+2+mint(3)
459  nstr=0
460  DO 210 i=mint(84)+1,nmax
461  kcs=kchg(pycomp(k(i,2)),2)*isign(1,k(i,2))
462  IF(kcs.EQ.0) goto 210
463  DO 200 j=1,4
464  IF(kcs.EQ.1.AND.(j.EQ.2.OR.j.EQ.4)) goto 200
465  IF(kcs.EQ.-1.AND.(j.EQ.1.OR.j.EQ.3)) goto 200
466  IF(j.LE.2) THEN
467  ist=mod(k(i,j+3)/mstu(5),mstu(5))
468  ELSE
469  ist=mod(k(i,j+1),mstu(5))
470  ENDIF
471  IF(ist.LT.mint(84).OR.ist.GT.i) goto 200
472  IF(kchg(pycomp(k(ist,2)),2).EQ.0) goto 200
473  nstr=nstr+1
474  IF(j.EQ.1.OR.j.EQ.4) THEN
475  kstr(nstr,1)=i
476  kstr(nstr,2)=ist
477  ELSE
478  kstr(nstr,1)=ist
479  kstr(nstr,2)=i
480  ENDIF
481  200 CONTINUE
482  210 CONTINUE
483 
484 C...Set up starting values for iteration in xT2.
485  xt2=4d0*vint(62)/vint(2)
486  IF(mstp(82).LE.1) THEN
487  sigrat=xsec(isub,1)/max(1d-10,vint(315)*vint(316)*sigt(0,0,5))
488  IF(mint(141).NE.0.OR.mint(142).NE.0) sigrat=sigrat*
489  & vint(317)/(vint(318)*vint(320))
490  xt2fac=sigrat*vint(149)/(1d0-vint(149))
491  ELSE
492  xt2fac=vnt146*vint(148)*xsec(isub,1)/
493  & max(1d-10,sigt(0,0,5))*vint(149)*(1d0+vint(149))
494  ENDIF
495  vint(63)=0d0
496  vint(64)=0d0
497  vint(143)=1d0-vint(141)
498  vint(144)=1d0-vint(142)
499 
500 C...Iterate downwards in xT2.
501  220 IF(mstp(82).LE.1) THEN
502  xt2=xt2fac*xt2/(xt2fac-xt2*log(pyr(0)))
503  IF(xt2.LT.vint(149)) goto 270
504  ELSE
505  IF(xt2.LE.0.01001d0*vint(149)) goto 270
506  xt2=xt2fac*(xt2+vint(149))/(xt2fac-(xt2+vint(149))*
507  & log(pyr(0)))-vint(149)
508  IF(xt2.LE.0d0) goto 270
509  xt2=max(0.01d0*vint(149),xt2)
510  ENDIF
511  vint(25)=xt2
512 
513 C...Choose tau and y*. Calculate cos(theta-hat).
514  IF(pyr(0).LE.coef(isub,1)) THEN
515  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
516  tau=xt2*(1d0+taut)**2/(4d0*taut)
517  ELSE
518  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
519  ENDIF
520  vint(21)=tau
521  CALL pyklim(2)
522  ryst=pyr(0)
523  myst=1
524  IF(ryst.GT.coef(isub,8)) myst=2
525  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
526  CALL pykmap(2,myst,pyr(0))
527  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
528 
529 C...Check that x not used up. Accept or reject kinematical variables.
530  x1m=sqrt(tau)*exp(vint(22))
531  x2m=sqrt(tau)*exp(-vint(22))
532  IF(vint(143)-x1m.LT.0.01d0.OR.vint(144)-x2m.LT.0.01d0) goto 220
533  vint(71)=0.5d0*vint(1)*sqrt(xt2)
534  CALL pysigh(nchn,sigs)
535  IF(mint(141).NE.0.OR.mint(142).NE.0) sigs=sigs*vint(320)
536  IF(sigs.LT.xsec(isub,1)*pyr(0)) goto 220
537 
538 C...Reset K, P and V vectors. Select some variables.
539  DO 240 i=n+1,n+2
540  DO 230 j=1,5
541  k(i,j)=0
542  p(i,j)=0d0
543  v(i,j)=0d0
544  230 CONTINUE
545  240 CONTINUE
546  rflav=pyr(0)
547  pt=0.5d0*vint(1)*sqrt(xt2)
548  phi=paru(2)*pyr(0)
549  cth=vint(23)
550 
551 C...Add first parton to event record.
552  k(n+1,1)=3
553  k(n+1,2)=21
554  IF(rflav.GE.max(parp(85),parp(86))) k(n+1,2)=
555  & 1+int((2d0+parj(2))*pyr(0))
556  p(n+1,1)=pt*cos(phi)
557  p(n+1,2)=pt*sin(phi)
558  p(n+1,3)=0.25d0*vint(1)*(vint(41)*(1d0+cth)-vint(42)*(1d0-cth))
559  p(n+1,4)=0.25d0*vint(1)*(vint(41)*(1d0+cth)+vint(42)*(1d0-cth))
560  p(n+1,5)=0d0
561 
562 C...Add second parton to event record.
563  k(n+2,1)=3
564  k(n+2,2)=21
565  IF(k(n+1,2).NE.21) k(n+2,2)=-k(n+1,2)
566  p(n+2,1)=-p(n+1,1)
567  p(n+2,2)=-p(n+1,2)
568  p(n+2,3)=0.25d0*vint(1)*(vint(41)*(1d0-cth)-vint(42)*(1d0+cth))
569  p(n+2,4)=0.25d0*vint(1)*(vint(41)*(1d0-cth)+vint(42)*(1d0+cth))
570  p(n+2,5)=0d0
571 
572  IF(rflav.LT.parp(85).AND.nstr.GE.1) THEN
573 C....Choose relevant string pieces to place gluons on.
574  DO 260 i=n+1,n+2
575  dmin=1d8
576  DO 250 istr=1,nstr
577  i1=kstr(istr,1)
578  i2=kstr(istr,2)
579  dist=(p(i,4)*p(i1,4)-p(i,1)*p(i1,1)-p(i,2)*p(i1,2)-
580  & p(i,3)*p(i1,3))*(p(i,4)*p(i2,4)-p(i,1)*p(i2,1)-
581  & p(i,2)*p(i2,2)-p(i,3)*p(i2,3))/max(1d0,p(i1,4)*p(i2,4)-
582  & p(i1,1)*p(i2,1)-p(i1,2)*p(i2,2)-p(i1,3)*p(i2,3))
583  IF(istr.EQ.1.OR.dist.LT.dmin) THEN
584  dmin=dist
585  ist1=i1
586  ist2=i2
587  istm=istr
588  ENDIF
589  250 CONTINUE
590 
591 C....Colour flow adjustments, new string pieces.
592  IF(k(ist1,4)/mstu(5).EQ.ist2) k(ist1,4)=mstu(5)*i+
593  & mod(k(ist1,4),mstu(5))
594  IF(mod(k(ist1,5),mstu(5)).EQ.ist2) k(ist1,5)=
595  & mstu(5)*(k(ist1,5)/mstu(5))+i
596  k(i,5)=mstu(5)*ist1
597  k(i,4)=mstu(5)*ist2
598  IF(k(ist2,5)/mstu(5).EQ.ist1) k(ist2,5)=mstu(5)*i+
599  & mod(k(ist2,5),mstu(5))
600  IF(mod(k(ist2,4),mstu(5)).EQ.ist1) k(ist2,4)=
601  & mstu(5)*(k(ist2,4)/mstu(5))+i
602  kstr(istm,2)=i
603  kstr(nstr+1,1)=i
604  kstr(nstr+1,2)=ist2
605  nstr=nstr+1
606  260 CONTINUE
607 
608 C...String drawing and colour flow for gluon loop.
609  ELSEIF(k(n+1,2).EQ.21) THEN
610  k(n+1,4)=mstu(5)*(n+2)
611  k(n+1,5)=mstu(5)*(n+2)
612  k(n+2,4)=mstu(5)*(n+1)
613  k(n+2,5)=mstu(5)*(n+1)
614  kstr(nstr+1,1)=n+1
615  kstr(nstr+1,2)=n+2
616  kstr(nstr+2,1)=n+2
617  kstr(nstr+2,2)=n+1
618  nstr=nstr+2
619 
620 C...String drawing and colour flow for qqbar pair.
621  ELSE
622  k(n+1,4)=mstu(5)*(n+2)
623  k(n+2,5)=mstu(5)*(n+1)
624  kstr(nstr+1,1)=n+1
625  kstr(nstr+1,2)=n+2
626  nstr=nstr+1
627  ENDIF
628 
629 C...Global statistics.
630  mint(351)=mint(351)+1
631  vint(351)=vint(351)+pt
632  IF (mint(351).EQ.1) vint(356)=pt
633 
634 C...Update remaining energy; iterate.
635  n=n+2
636  IF(n.GT.mstu(4)-mstu(32)-10) THEN
637  CALL pyerrm(11,'(PYMULT:) no more memory left in PYJETS')
638  mint(51)=1
639  RETURN
640  ENDIF
641  mint(31)=mint(31)+1
642  vint(151)=vint(151)+vint(41)
643  vint(152)=vint(152)+vint(42)
644  vint(143)=vint(143)-vint(41)
645  vint(144)=vint(144)-vint(42)
646 C...Allow FSR for UE
647  IF(mstp(152).EQ.1) CALL pyshow(n-1,n,sqrt(parp(71))*pt)
648  IF(mint(31).LT.240) goto 220
649  270 CONTINUE
650  mint(1)=isubsv
651  DO 280 j=11,80
652  vint(j)=vintsv(j)
653  280 CONTINUE
654  ENDIF
655 
656 C...Format statements for printout.
657  5000 FORMAT(/1x,'****** PYMULT: initialization of multiple inter',
658  &'actions for MSTP(82) =',i2,' ******')
659  5100 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,
660  &d9.2,' mb: rejected')
661  5200 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,
662  &d9.2,' mb: accepted')
663 
664  RETURN
665  END