Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymign.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymign.f
1 
2 C*********************************************************************
3 
4 C...PYMIGN
5 C...Initializes treatment of new multiple interactions scenario,
6 C...selects kinematics of hardest interaction if low-pT physics
7 C...included in run, and generates all non-hardest interactions.
8 
9  SUBROUTINE pymign(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  EXTERNAL pyalps
16  DOUBLE PRECISION pyalps
17 C...Commonblocks.
18  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
19  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
20  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
21  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
22  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
23  common/pypars/mstp(200),parp(200),msti(200),pari(200)
24  common/pyint1/mint(400),vint(400)
25  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
26  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
27  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
28  common/pyint7/sigt(0:6,0:6,0:5)
29  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
30  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
31  & xmi(2,240),pt2mi(240),imisep(0:240)
32  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,
33  &/pyint1/,/pyint2/,/pyint3/,/pyint5/,/pyint7/,/pyintm/
34 C...Local arrays and saved variables.
35  dimension nmul(20),sigm(20),kstr(500,2),vintsv(80),
36  &wdtp(0:400),wdte(0:400,0:5),xpq(-25:25),ksav(4,5),psav(4,5)
37  SAVE xt2,xt2fac,xc2,xts,irbin,rbin,nmul,sigm,p83a,p83b,p83c,
38  &cq2i,cq2r,pik,bdiv,b,plowb,phighb,pallb,s4a,s4b,s4c,powip,
39  &rpwip,b2rpdv,b2rpmx,bavg,vnt145,vnt146,vnt147
40 
41 C...Initialization of multiple interaction treatment.
42  IF(mmul.EQ.1) THEN
43  IF(mstp(122).GE.1) WRITE(mstu(11),5000) mstp(82)
44  isub=96
45  mint(1)=96
46  vint(63)=0d0
47  vint(64)=0d0
48  vint(143)=1d0
49  vint(144)=1d0
50 
51 C...Loop over phase space points: xT2 choice in 20 bins.
52  100 sigsum=0d0
53  DO 120 ixt2=1,20
54  nmul(ixt2)=mstp(83)
55  sigm(ixt2)=0d0
56  DO 110 itry=1,mstp(83)
57  rsca=0.05d0*((21-ixt2)-pyr(0))
58  xt2=vint(149)*(1d0+vint(149))/(vint(149)+rsca)-vint(149)
59  xt2=max(0.01d0*vint(149),xt2)
60  vint(25)=xt2
61 
62 C...Choose tau and y*. Calculate cos(theta-hat).
63  IF(pyr(0).LE.coef(isub,1)) THEN
64  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
65  tau=xt2*(1d0+taut)**2/(4d0*taut)
66  ELSE
67  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
68  ENDIF
69  vint(21)=tau
70  CALL pyklim(2)
71  ryst=pyr(0)
72  myst=1
73  IF(ryst.GT.coef(isub,8)) myst=2
74  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
75  CALL pykmap(2,myst,pyr(0))
76  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
77 
78 C...Calculate differential cross-section.
79  vint(71)=0.5d0*vint(1)*sqrt(xt2)
80  CALL pysigh(nchn,sigs)
81  sigm(ixt2)=sigm(ixt2)+sigs
82  110 CONTINUE
83  sigsum=sigsum+sigm(ixt2)
84  120 CONTINUE
85  sigsum=sigsum/(20d0*mstp(83))
86 
87 C...Reject result if sigma(parton-parton) is smaller than hadronic one.
88  IF(sigsum.LT.1.1d0*sigt(0,0,5)) THEN
89  IF(mstp(122).GE.1) WRITE(mstu(11),5100)
90  & parp(82)*(vint(1)/parp(89))**parp(90),sigsum
91  parp(82)=0.9d0*parp(82)
92  vint(149)=4d0*(parp(82)*(vint(1)/parp(89))**parp(90))**2/
93  & vint(2)
94  goto 100
95  ENDIF
96  IF(mstp(122).GE.1) WRITE(mstu(11),5200)
97  & parp(82)*(vint(1)/parp(89))**parp(90), sigsum
98 
99 C...Start iteration to find k factor.
100  yke=sigsum/max(1d-10,sigt(0,0,5))
101  p83a=(1d0-parp(83))**2
102  p83b=2d0*parp(83)*(1d0-parp(83))
103  p83c=parp(83)**2
104  cq2i=1d0/parp(84)**2
105  cq2r=2d0/(1d0+parp(84)**2)
106  so=0.5d0
107  xi=0d0
108  yi=0d0
109  xf=0d0
110  yf=0d0
111  xk=0.5d0
112  iit=0
113  130 IF(iit.EQ.0) THEN
114  xk=2d0*xk
115  ELSEIF(iit.EQ.1) THEN
116  xk=0.5d0*xk
117  ELSE
118  xk=xi+(yke-yi)*(xf-xi)/(yf-yi)
119  ENDIF
120 
121 C...Evaluate overlap integrals. Find where to divide the b range.
122  IF(mstp(82).EQ.2) THEN
123  sp=0.5d0*paru(1)*(1d0-exp(-xk))
124  sop=sp/paru(1)
125  ELSE
126  IF(mstp(82).EQ.3) THEN
127  deltab=0.02d0
128  ELSEIF(mstp(82).EQ.4) THEN
129  deltab=min(0.01d0,0.05d0*parp(84))
130  ELSE
131  powip=max(0.4d0,parp(83))
132  rpwip=2d0/powip-1d0
133  deltab=max(0.02d0,0.02d0*(2d0/powip)**(1d0/powip))
134  so=0d0
135  ENDIF
136  sp=0d0
137  sop=0d0
138  bsp=0d0
139  sohigh=0d0
140  ibdiv=0
141  b=-0.5d0*deltab
142  140 b=b+deltab
143  IF(mstp(82).EQ.3) THEN
144  ov=exp(-b**2)/paru(2)
145  ELSEIF(mstp(82).EQ.4) THEN
146  ov=(p83a*exp(-min(50d0,b**2))+
147  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
148  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
149  ELSE
150  ov=exp(-b**powip)/paru(2)
151  so=so+paru(2)*b*deltab*ov
152  ENDIF
153  IF(ibdiv.EQ.1) sohigh=sohigh+paru(2)*b*deltab*ov
154  pacc=1d0-exp(-min(50d0,paru(1)*xk*ov))
155  sp=sp+paru(2)*b*deltab*pacc
156  sop=sop+paru(2)*b*deltab*ov*pacc
157  bsp=bsp+b*paru(2)*b*deltab*pacc
158  IF(ibdiv.EQ.0.AND.paru(1)*xk*ov.LT.1d0) THEN
159  ibdiv=1
160  bdiv=b+0.5d0*deltab
161  ENDIF
162  IF(b.LT.1d0.OR.b*pacc.GT.1d-6) goto 140
163  ENDIF
164  yk=paru(1)*xk*so/sp
165 
166 C...Continue iteration until convergence.
167  IF(yk.LT.yke) THEN
168  xi=xk
169  yi=yk
170  IF(iit.EQ.1) iit=2
171  ELSE
172  xf=xk
173  yf=yk
174  IF(iit.EQ.0) iit=1
175  ENDIF
176  IF(abs(yk-yke).GE.1d-5*yke) goto 130
177 
178 C...Store some results for subsequent use.
179  bavg=bsp/sp
180  vint(145)=sigsum
181  vint(146)=sop/so
182  vint(147)=sop/sp
183  vnt145=vint(145)
184  vnt146=vint(146)
185  vnt147=vint(147)
186 C...PIK = PARU(1)*XK = (VINT(146)/VINT(147))*sigma_jet/sigma_nondiffr.
187  pik=(vnt146/vnt147)*yke
188 
189 C...Find relative weight for low and high impact parameter..
190  plowb=paru(1)*bdiv**2
191  IF(mstp(82).EQ.3) THEN
192  phighb=pik*0.5*exp(-bdiv**2)
193  ELSEIF(mstp(82).EQ.4) THEN
194  s4a=p83a*exp(-bdiv**2)
195  s4b=p83b*exp(-bdiv**2*cq2r)
196  s4c=p83c*exp(-bdiv**2*cq2i)
197  phighb=pik*0.5*(s4a+s4b+s4c)
198  ELSEIF(parp(83).GE.1.999d0) THEN
199  phighb=pik*sohigh
200  b2rpdv=bdiv**powip
201  ELSE
202  phighb=pik*sohigh
203  b2rpdv=bdiv**powip
204  b2rpmx=max(2d0*rpwip,b2rpdv)
205  ENDIF
206  pallb=plowb+phighb
207 
208 C...Initialize iteration in xT2 for hardest interaction.
209  ELSEIF(mmul.EQ.2) THEN
210  vint(145)=vnt145
211  vint(146)=vnt146
212  vint(147)=vnt147
213  IF(mstp(82).LE.0) THEN
214  ELSEIF(mstp(82).EQ.1) THEN
215  xt2=1d0
216  sigrat=xsec(96,1)/max(1d-10,vint(315)*vint(316)*sigt(0,0,5))
217  IF(mint(141).NE.0.OR.mint(142).NE.0) sigrat=sigrat*
218  & vint(317)/(vint(318)*vint(320))
219  xt2fac=sigrat*vint(149)/(1d0-vint(149))
220  ELSEIF(mstp(82).EQ.2) THEN
221  xt2=1d0
222  xt2fac=vnt146*xsec(96,1)/max(1d-10,sigt(0,0,5))*
223  & vint(149)*(1d0+vint(149))
224  ELSE
225  xc2=4d0*ckin(3)**2/vint(2)
226  IF(ckin(3).LE.ckin(5).OR.mint(82).GE.2) xc2=0d0
227  ENDIF
228 
229 C...Select impact parameter for hardest interaction.
230  IF(mstp(82).LE.2) RETURN
231  142 IF(pyr(0)*pallb.LT.plowb) THEN
232 C...Treatment in low b region.
233  mint(39)=1
234  b=bdiv*sqrt(pyr(0))
235  IF(mstp(82).EQ.3) THEN
236  ov=exp(-b**2)/paru(2)
237  ELSEIF(mstp(82).EQ.4) THEN
238  ov=(p83a*exp(-min(50d0,b**2))+
239  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
240  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
241  ELSE
242  ov=exp(-b**powip)/paru(2)
243  ENDIF
244  vint(148)=ov/vnt147
245  pacc=1d0-exp(-min(50d0,pik*ov))
246  xt2=1d0
247  xt2fac=vnt146*vint(148)*xsec(96,1)/max(1d-10,sigt(0,0,5))*
248  & vint(149)*(1d0+vint(149))
249  ELSE
250 C...Treatment in high b region.
251  mint(39)=2
252  IF(mstp(82).EQ.3) THEN
253  b=sqrt(bdiv**2-log(pyr(0)))
254  ov=exp(-b**2)/paru(2)
255  ELSEIF(mstp(82).EQ.4) THEN
256  s4rndm=pyr(0)*(s4a+s4b+s4c)
257  IF(s4rndm.LT.s4a) THEN
258  b=sqrt(bdiv**2-log(pyr(0)))
259  ELSEIF(s4rndm.LT.s4a+s4b) THEN
260  b=sqrt(bdiv**2-log(pyr(0))/cq2r)
261  ELSE
262  b=sqrt(bdiv**2-log(pyr(0))/cq2i)
263  ENDIF
264  ov=(p83a*exp(-min(50d0,b**2))+
265  & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
266  & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
267  ELSEIF(parp(83).GE.1.999d0) THEN
268  144 b2rpw=b2rpdv-log(pyr(0))
269  accip=(b2rpw/b2rpdv)**rpwip
270  IF(accip.LT.pyr(0)) goto 144
271  ov=exp(-b2rpw)/paru(2)
272  b=b2rpw**(1d0/powip)
273  ELSE
274  146 b2rpw=b2rpdv-2d0*log(pyr(0))
275  accip=(b2rpw/b2rpmx)**rpwip*exp(-0.5d0*(b2rpw-b2rpmx))
276  IF(accip.LT.pyr(0)) goto 146
277  ov=exp(-b2rpw)/paru(2)
278  b=b2rpw**(1d0/powip)
279  ENDIF
280  vint(148)=ov/vnt147
281  pacc=(1d0-exp(-min(50d0,pik*ov)))/(pik*ov)
282  ENDIF
283  IF(pacc.LT.pyr(0)) goto 142
284  vint(139)=b/bavg
285 
286  ELSEIF(mmul.EQ.3) THEN
287 C...Low-pT or multiple interactions (first semihard interaction):
288 C...choose xT2 according to dpT2/pT2**2*exp(-(sigma above pT2)/norm)
289 C...or (MSTP(82)>=2) dpT2/(pT2+pT0**2)**2*exp(-....).
290  isub=mint(1)
291  vint(145)=vnt145
292  vint(146)=vnt146
293  vint(147)=vnt147
294  IF(mstp(82).LE.0) THEN
295  xt2=0d0
296  ELSEIF(mstp(82).EQ.1) THEN
297  xt2=xt2fac*xt2/(xt2fac-xt2*log(pyr(0)))
298 C...Use with "Sudakov" for low b values when impact parameter dependence.
299  ELSEIF(mstp(82).EQ.2.OR.mint(39).EQ.1) THEN
300  IF(xt2.LT.1d0.AND.exp(-xt2fac*xt2/(vint(149)*(xt2+
301  & vint(149)))).GT.pyr(0)) xt2=1d0
302  IF(xt2.GE.1d0) THEN
303  xt2=(1d0+vint(149))*xt2fac/(xt2fac-(1d0+vint(149))*log(1d0-
304  & pyr(0)*(1d0-exp(-xt2fac/(vint(149)*(1d0+vint(149)))))))-
305  & vint(149)
306  ELSE
307  xt2=-xt2fac/log(exp(-xt2fac/(xt2+vint(149)))+pyr(0)*
308  & (exp(-xt2fac/vint(149))-exp(-xt2fac/(xt2+vint(149)))))-
309  & vint(149)
310  ENDIF
311  xt2=max(0.01d0*vint(149),xt2)
312 C...Use without "Sudakov" for high b values when impact parameter dep.
313  ELSE
314  xt2=(xc2+vint(149))*(1d0+vint(149))/(1d0+vint(149)-
315  & pyr(0)*(1d0-xc2))-vint(149)
316  xt2=max(0.01d0*vint(149),xt2)
317  ENDIF
318  vint(25)=xt2
319 
320 C...Low-pT: choose xT2, tau, y* and cos(theta-hat) fixed.
321  IF(mstp(82).LE.1.AND.xt2.LT.vint(149)) THEN
322  IF(mint(82).EQ.1) ngen(0,1)=ngen(0,1)-mint(143)
323  IF(mint(82).EQ.1) ngen(isub,1)=ngen(isub,1)-mint(143)
324  isub=95
325  mint(1)=isub
326  vint(21)=1d-12*vint(149)
327  vint(22)=0d0
328  vint(23)=0d0
329  vint(25)=1d-12*vint(149)
330 
331  ELSE
332 C...Multiple interactions (first semihard interaction).
333 C...Choose tau and y*. Calculate cos(theta-hat).
334  IF(pyr(0).LE.coef(isub,1)) THEN
335  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
336  tau=xt2*(1d0+taut)**2/(4d0*taut)
337  ELSE
338  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
339  ENDIF
340  vint(21)=tau
341  CALL pyklim(2)
342  ryst=pyr(0)
343  myst=1
344  IF(ryst.GT.coef(isub,8)) myst=2
345  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
346  CALL pykmap(2,myst,pyr(0))
347  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
348  ENDIF
349  vint(71)=0.5d0*vint(1)*sqrt(vint(25))
350 
351 C...Store results of cross-section calculation.
352  ELSEIF(mmul.EQ.4) THEN
353  isub=mint(1)
354  vint(145)=vnt145
355  vint(146)=vnt146
356  vint(147)=vnt147
357  xts=vint(25)
358  IF(iset(isub).EQ.1) xts=vint(21)
359  IF(iset(isub).EQ.2)
360  & xts=(4d0*vint(48)+2d0*vint(63)+2d0*vint(64))/vint(2)
361  IF(iset(isub).GE.3.AND.iset(isub).LE.5) xts=vint(26)
362  rbin=max(0.000001d0,min(0.999999d0,xts*(1d0+vint(149))/
363  & (xts+vint(149))))
364  irbin=int(1d0+20d0*rbin)
365  IF(isub.EQ.96.AND.mstp(171).EQ.0) THEN
366  nmul(irbin)=nmul(irbin)+1
367  sigm(irbin)=sigm(irbin)+vint(153)
368  ENDIF
369 
370 C...Choose impact parameter if not already done.
371  ELSEIF(mmul.EQ.5) THEN
372  isub=mint(1)
373  vint(145)=vnt145
374  vint(146)=vnt146
375  vint(147)=vnt147
376  150 IF(mint(39).GT.0) THEN
377  ELSEIF(mstp(82).EQ.3) THEN
378  expb2=pyr(0)
379  b2=-log(pyr(0))
380  vint(148)=expb2/(paru(2)*vnt147)
381  vint(139)=sqrt(b2)/bavg
382  ELSEIF(mstp(82).EQ.4) THEN
383  rtype=pyr(0)
384  IF(rtype.LT.p83a) THEN
385  b2=-log(pyr(0))
386  ELSEIF(rtype.LT.p83a+p83b) THEN
387  b2=-log(pyr(0))/cq2r
388  ELSE
389  b2=-log(pyr(0))/cq2i
390  ENDIF
391  vint(148)=(p83a*exp(-min(50d0,b2))+
392  & p83b*cq2r*exp(-min(50d0,b2*cq2r))+
393  & p83c*cq2i*exp(-min(50d0,b2*cq2i)))/(paru(2)*vnt147)
394  vint(139)=sqrt(b2)/bavg
395  ELSEIF(parp(83).GE.1.999d0) THEN
396  powip=max(2d0,parp(83))
397  rpwip=2d0/powip-1d0
398  prob1=powip/(2d0*exp(-1d0)+powip)
399  160 IF(pyr(0).LT.prob1) THEN
400  b2rpw=pyr(0)**(0.5d0*powip)
401  accip=exp(-b2rpw)
402  ELSE
403  b2rpw=1d0-log(pyr(0))
404  accip=b2rpw**rpwip
405  ENDIF
406  IF(accip.LT.pyr(0)) goto 160
407  vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
408  vint(139)=b2rpw**(1d0/powip)/bavg
409  ELSE
410  powip=max(0.4d0,parp(83))
411  rpwip=2d0/powip-1d0
412  prob1=rpwip/(rpwip+2d0**rpwip*exp(-rpwip))
413  170 IF(pyr(0).LT.prob1) THEN
414  b2rpw=2d0*rpwip*pyr(0)
415  accip=(b2rpw/rpwip)**rpwip*exp(rpwip-b2rpw)
416  ELSE
417  b2rpw=2d0*(rpwip-log(pyr(0)))
418  accip=(0.5d0*b2rpw/rpwip)**rpwip*exp(rpwip-0.5d0*b2rpw)
419  ENDIF
420  IF(accip.lt .pyr(0)) goto 170
421  vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
422  vint(139)=b2rpw**(1d0/powip)/bavg
423  ENDIF
424 
425 C...Multiple interactions (variable impact parameter) : reject with
426 C...probability exp(-overlap*cross-section above pT/normalization).
427 C...Does not apply to low-b region, where "Sudakov" already included.
428  vint(150)=1d0
429  IF(mint(39).NE.1) THEN
430  rncor=(irbin-20d0*rbin)*nmul(irbin)
431  sigcor=(irbin-20d0*rbin)*sigm(irbin)
432  DO 180 ibin=irbin+1,20
433  rncor=rncor+nmul(ibin)
434  sigcor=sigcor+sigm(ibin)
435  180 CONTINUE
436  sigabv=(sigcor/rncor)*vint(149)*(1d0-xts)/(xts+vint(149))
437  IF(mstp(171).EQ.1) sigabv=sigabv*vint(2)/vint(289)
438  vint(150)=exp(-min(50d0,vnt146*vint(148)*
439  & sigabv/max(1d-10,sigt(0,0,5))))
440  ENDIF
441  IF(mstp(86).EQ.3.OR.(mstp(86).EQ.2.AND.isub.NE.11.AND.
442  & isub.NE.12.AND.isub.NE.13.AND.isub.NE.28.AND.isub.NE.53
443  & .AND.isub.NE.68.AND.isub.NE.95.AND.isub.NE.96)) THEN
444  IF(vint(150).LT.pyr(0)) goto 150
445  vint(150)=1d0
446  ENDIF
447 
448 C...Generate additional multiple semihard interactions.
449  ELSEIF(mmul.EQ.6) THEN
450 
451 C...Save data for hardest initeraction, to be restored.
452  isubsv=mint(1)
453  vint(145)=vnt145
454  vint(146)=vnt146
455  vint(147)=vnt147
456  m13sv=mint(13)
457  m14sv=mint(14)
458  m15sv=mint(15)
459  m16sv=mint(16)
460  m21sv=mint(21)
461  m22sv=mint(22)
462  DO 190 j=11,80
463  vintsv(j)=vint(j)
464  190 CONTINUE
465  v141sv=vint(141)
466  v142sv=vint(142)
467 
468 C...Store data on hardest interaction.
469  xmi(1,1)=vint(141)
470  xmi(2,1)=vint(142)
471  pt2mi(1)=vint(54)
472  imisep(0)=mint(84)
473  imisep(1)=n
474 
475 C...Change process to generate; sum of x values so far.
476  isub=96
477  mint(1)=96
478  vint(143)=1d0-vint(141)
479  vint(144)=1d0-vint(142)
480  vint(151)=0d0
481  vint(152)=0d0
482 
483 C...Initialize factors for PDF reshaping.
484  DO 230 js=1,2
485  kfbeam=mint(10+js)
486  kfabm=iabs(kfbeam)
487  kfsbm=isign(1,kfbeam)
488 
489 C...Zero flavour content of incoming beam particle.
490  kfival(js,1)=0
491  kfival(js,2)=0
492  kfival(js,3)=0
493 C...Flavour content of baryon.
494  IF(kfabm.GT.1000) THEN
495  kfival(js,1)=kfsbm*mod(kfabm/1000,10)
496  kfival(js,2)=kfsbm*mod(kfabm/100,10)
497  kfival(js,3)=kfsbm*mod(kfabm/10,10)
498 C...Flavour content of pi+-, K+-.
499  ELSEIF(kfabm.EQ.211) THEN
500  kfival(js,1)=kfsbm*2
501  kfival(js,2)=-kfsbm
502  ELSEIF(kfabm.EQ.321) THEN
503  kfival(js,1)=-kfsbm*3
504  kfival(js,2)=kfsbm*2
505 C...Flavour content of pi0, gamma, K0S, K0L not defined yet.
506  ENDIF
507 
508 C...Zero initial valence and companion content.
509  DO 200 ifl=-6,6
510  nvc(js,ifl)=0
511  200 CONTINUE
512 
513 C...Initiate listing of all incoming partons from two sides.
514  nmi(js)=0
515  DO 210 i=mint(84)+1,n
516  IF(k(i,3).EQ.mint(83)+2+js) THEN
517  imi(js,1,1)=i
518  imi(js,1,2)=0
519  ENDIF
520  210 CONTINUE
521 
522 C...Decide whether quarks in hard scattering were valence or sea.
523  ifl=k(imi(js,1,1),2)
524  IF (iabs(ifl).GT.6) goto 230
525 
526 C...Get PDFs at X and Q2 of the parton shower initiator for the
527 C...hard scattering.
528  x=vint(140+js)
529  IF(mstp(61).GE.1) THEN
530  q2=parp(62)**2
531  ELSE
532  q2=vint(54)
533  ENDIF
534 C...Note: XPSVC = x*pdf.
535  mint(30)=js
536  CALL pypdfu(kfbeam,x,q2,xpq)
537  sea=xpsvc(ifl,-1)
538  val=xpsvc(ifl,0)
539 
540 C...Decide (Extra factor x cancels in the division).
541  rvcs=pyr(0)*(sea+val)
542  ivnow=1
543  220 IF (rvcs.LE.val.AND.ivnow.GE.1) THEN
544 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
545  ivnow=0
546  IF(kfival(js,1).EQ.ifl) ivnow=ivnow+1
547  IF(kfival(js,2).EQ.ifl) ivnow=ivnow+1
548  IF(kfival(js,3).EQ.ifl) ivnow=ivnow+1
549  IF(kfival(js,1).EQ.0) THEN
550  IF(kfbeam.EQ.111.AND.iabs(ifl).LE.2) ivnow=1
551  IF(kfbeam.EQ.22.AND.iabs(ifl).LE.5) ivnow=1
552  IF((kfbeam.EQ.130.OR.kfbeam.EQ.310).AND.
553  & (iabs(ifl).EQ.1.OR.iabs(ifl).EQ.3)) ivnow=1
554  ENDIF
555  IF(ivnow.EQ.0) goto 220
556 C...Mark valence.
557  imi(js,1,2)=0
558 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
559  IF(kfival(js,1).EQ.0) THEN
560  IF(kfbeam.EQ.111.OR.kfbeam.EQ.22) THEN
561  kfival(js,1)=ifl
562  kfival(js,2)=-ifl
563  ELSEIF(kfbeam.EQ.130.OR.kfbeam.EQ.310) THEN
564  kfival(js,1)=ifl
565  IF(iabs(ifl).EQ.1) kfival(js,2)=isign(3,-ifl)
566  IF(iabs(ifl).NE.1) kfival(js,2)=isign(1,-ifl)
567  ENDIF
568  ENDIF
569 
570 C...If sea, add opposite sign companion parton. Store X and I.
571  ELSE
572  nvc(js,-ifl)=nvc(js,-ifl)+1
573  xassoc(js,-ifl,nvc(js,-ifl))=x
574 C...Set pointer to companion
575  imi(js,1,2)=-nvc(js,-ifl)
576  ENDIF
577  230 CONTINUE
578 
579 C...Update counter number of multiple interactions.
580  nmi(1)=1
581  nmi(2)=1
582 
583 C...Set up starting values for iteration in xT2.
584  IF(mstp(86).EQ.3.OR.(mstp(86).EQ.2.AND.isubsv.NE.11.AND.
585  & isubsv.NE.12.AND.isubsv.NE.13.AND.isubsv.NE.28.AND.
586  & isubsv.NE.53.AND.isubsv.NE.68.AND.isubsv.NE.95.AND.
587  & isubsv.NE.96)) THEN
588  xt2=(1d0-vint(141))*(1d0-vint(142))
589  ELSE
590  xt2=vint(25)
591  IF(iset(isubsv).EQ.1) xt2=vint(21)
592  IF(iset(isubsv).EQ.2)
593  & xt2=(4d0*vint(48)+2d0*vint(63)+2d0*vint(64))/vint(2)
594  IF(iset(isubsv).GE.3.AND.iset(isubsv).LE.5) xt2=vint(26)
595  ENDIF
596  IF(mstp(82).LE.1) THEN
597  sigrat=xsec(isub,1)/max(1d-10,vint(315)*vint(316)*sigt(0,0,5))
598  IF(mint(141).NE.0.OR.mint(142).NE.0) sigrat=sigrat*
599  & vint(317)/(vint(318)*vint(320))
600  xt2fac=sigrat*vint(149)/(1d0-vint(149))
601  ELSE
602  xt2fac=vnt146*vint(148)*xsec(isub,1)/
603  & max(1d-10,sigt(0,0,5))*vint(149)*(1d0+vint(149))
604  ENDIF
605  vint(63)=0d0
606  vint(64)=0d0
607 
608 C...Iterate downwards in xT2.
609  240 IF((mint(35).EQ.2.AND.mstp(81).EQ.10).OR.isubsv.EQ.95) THEN
610  xt2=0d0
611  goto 440
612  ELSEIF(mstp(82).LE.1) THEN
613  xt2=xt2fac*xt2/(xt2fac-xt2*log(pyr(0)))
614  IF(xt2.LT.vint(149)) goto 440
615  ELSE
616  IF(xt2.LE.0.01001d0*vint(149)) goto 440
617  xt2=xt2fac*(xt2+vint(149))/(xt2fac-(xt2+vint(149))*
618  & log(pyr(0)))-vint(149)
619  IF(xt2.LE.0d0) goto 440
620  xt2=max(0.01d0*vint(149),xt2)
621  ENDIF
622  vint(25)=xt2
623 
624 C...Choose tau and y*. Calculate cos(theta-hat).
625  IF(pyr(0).LE.coef(isub,1)) THEN
626  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
627  tau=xt2*(1d0+taut)**2/(4d0*taut)
628  ELSE
629  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
630  ENDIF
631  vint(21)=tau
632 C...New: require shat > 1.
633  IF(tau*vint(2).LT.1d0) goto 240
634  CALL pyklim(2)
635  ryst=pyr(0)
636  myst=1
637  IF(ryst.GT.coef(isub,8)) myst=2
638  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
639  CALL pykmap(2,myst,pyr(0))
640  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
641 
642 C...Check that x not used up. Accept or reject kinematical variables.
643  x1m=sqrt(tau)*exp(vint(22))
644  x2m=sqrt(tau)*exp(-vint(22))
645  IF(vint(143)-x1m.LT.0.01d0.OR.vint(144)-x2m.LT.0.01d0) goto 240
646  vint(71)=0.5d0*vint(1)*sqrt(xt2)
647  CALL pysigh(nchn,sigs)
648  IF(mint(141).NE.0.OR.mint(142).NE.0) sigs=sigs*vint(320)
649  IF(sigs.LT.xsec(isub,1)*pyr(0)) goto 240
650  IF(mint(141).NE.0.OR.mint(142).NE.0) sigs=sigs/vint(320)
651 
652 C...Reset K, P and V vectors.
653  DO 260 i=n+1,n+4
654  DO 250 j=1,5
655  k(i,j)=0
656  p(i,j)=0d0
657  v(i,j)=0d0
658  250 CONTINUE
659  260 CONTINUE
660  pt=0.5d0*vint(1)*sqrt(xt2)
661 
662 C...Choose flavour of reacting partons (and subprocess).
663  rsigs=sigs*pyr(0)
664  DO 270 ichn=1,nchn
665  kfl1=isig(ichn,1)
666  kfl2=isig(ichn,2)
667  iconmi=isig(ichn,3)
668  rsigs=rsigs-sigh(ichn)
669  IF(rsigs.LE.0d0) goto 280
670  270 CONTINUE
671 
672 C...Reassign to appropriate process codes.
673  280 isubmi=iconmi/10
674  iconmi=mod(iconmi,10)
675 
676 C...Choose new quark flavour for annihilation graphs
677  IF(isubmi.EQ.12.OR.isubmi.EQ.53) THEN
678  sh=tau*vint(2)
679  CALL pywidt(21,sh,wdtp,wdte)
680  290 rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*pyr(0)
681  DO 300 i=1,mdcy(21,3)
682  kflf=kfdp(i+mdcy(21,2)-1,1)
683  rkfl=rkfl-(wdte(i,1)+wdte(i,2)+wdte(i,4))
684  IF(rkfl.LE.0d0) goto 310
685  300 CONTINUE
686  310 IF(isubmi.EQ.53.AND.iconmi.LE.2) THEN
687  IF(kflf.GE.4) goto 290
688  ELSEIF(isubmi.EQ.53.AND.iconmi.LE.4) THEN
689  kflf=4
690  iconmi=iconmi-2
691  ELSEIF(isubmi.EQ.53) THEN
692  kflf=5
693  iconmi=iconmi-4
694  ENDIF
695  ENDIF
696 
697 C...Final state flavours and colour flow: default values
698  js=1
699  kfl3=kfl1
700  kfl4=kfl2
701  kcc=20
702  kcs=isign(1,kfl1)
703 
704  IF(isubmi.EQ.11) THEN
705 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2
706  kcc=iconmi
707  IF(kfl1*kfl2.LT.0) kcc=kcc+2
708 
709  ELSEIF(isubmi.EQ.12) THEN
710 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2
711  kfl3=isign(kflf,kfl1)
712  kfl4=-kfl3
713  kcc=4
714 
715  ELSEIF(isubmi.EQ.13) THEN
716 C...f + fbar -> g + g; th arbitrary
717  kfl3=21
718  kfl4=21
719  kcc=iconmi+4
720 
721  ELSEIF(isubmi.EQ.28) THEN
722 C...f + g -> f + g; th = (p(f)-p(f))**2
723  IF(kfl1.EQ.21) js=2
724  kcc=iconmi+6
725  IF(kfl1.EQ.21) kcc=kcc+2
726  IF(kfl1.NE.21) kcs=isign(1,kfl1)
727  IF(kfl2.NE.21) kcs=isign(1,kfl2)
728 
729  ELSEIF(isubmi.EQ.53) THEN
730 C...g + g -> f + fbar; th arbitrary
731  kcs=(-1)**int(1.5d0+pyr(0))
732  kfl3=isign(kflf,kcs)
733  kfl4=-kfl3
734  kcc=iconmi+10
735 
736  ELSEIF(isubmi.EQ.68) THEN
737 C...g + g -> g + g; th arbitrary
738  kcc=iconmi+12
739  kcs=(-1)**int(1.5d0+pyr(0))
740  ENDIF
741 
742 C...Store flavours of scattering.
743  mint(13)=kfl1
744  mint(14)=kfl2
745  mint(15)=kfl1
746  mint(16)=kfl2
747  mint(21)=kfl3
748  mint(22)=kfl4
749 
750 C...Set flavours and mothers of scattering partons.
751  k(n+1,1)=14
752  k(n+2,1)=14
753  k(n+3,1)=3
754  k(n+4,1)=3
755  k(n+1,2)=kfl1
756  k(n+2,2)=kfl2
757  k(n+3,2)=kfl3
758  k(n+4,2)=kfl4
759  k(n+1,3)=mint(83)+1
760  k(n+2,3)=mint(83)+2
761  k(n+3,3)=n+1
762  k(n+4,3)=n+2
763 
764 C...Store colour connection indices.
765  DO 320 j=1,2
766  jc=j
767  IF(kcs.EQ.-1) jc=3-j
768  IF(icol(kcc,1,jc).NE.0) k(n+1,j+3)=n+icol(kcc,1,jc)
769  IF(icol(kcc,2,jc).NE.0) k(n+2,j+3)=n+icol(kcc,2,jc)
770  IF(icol(kcc,3,jc).NE.0) k(n+3,j+3)=mstu(5)*(n+icol(kcc,3,jc))
771  IF(icol(kcc,4,jc).NE.0) k(n+4,j+3)=mstu(5)*(n+icol(kcc,4,jc))
772  320 CONTINUE
773 
774 C...Store incoming and outgoing partons in their CM-frame.
775  shr=sqrt(tau)*vint(1)
776  p(n+1,3)=0.5d0*shr
777  p(n+1,4)=0.5d0*shr
778  p(n+2,3)=-0.5d0*shr
779  p(n+2,4)=0.5d0*shr
780  p(n+3,5)=pymass(k(n+3,2))
781  p(n+4,5)=pymass(k(n+4,2))
782  IF(p(n+3,5)+p(n+4,5).GE.shr) goto 240
783  p(n+3,4)=0.5d0*(shr+(p(n+3,5)**2-p(n+4,5)**2)/shr)
784  p(n+3,3)=sqrt(max(0d0,p(n+3,4)**2-p(n+3,5)**2))
785  p(n+4,4)=shr-p(n+3,4)
786  p(n+4,3)=-p(n+3,3)
787 
788 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4)
789  phi=paru(2)*pyr(0)
790  CALL pyrobo(n+3,n+4,acos(vint(23)),phi,0d0,0d0,0d0)
791 
792 C...Set up default values before showers.
793  mint(31)=mint(31)+1
794  ipu1=n+1
795  ipu2=n+2
796  ipu3=n+3
797  ipu4=n+4
798  vint(141)=vint(41)
799  vint(142)=vint(42)
800  n=n+4
801 
802 C...Showering of initial state partons (optional).
803 C...Note: no showering of final state partons here; it comes later.
804  IF(mstp(84).GE.1.AND.mstp(61).GE.1) THEN
805  mint(51)=0
806  alamsv=parj(81)
807  parj(81)=parp(72)
808  nsav=n
809  DO 340 i=1,4
810  DO 330 j=1,5
811  ksav(i,j)=k(n-4+i,j)
812  psav(i,j)=p(n-4+i,j)
813  330 CONTINUE
814  340 CONTINUE
815  CALL pysspa(ipu1,ipu2)
816  parj(81)=alamsv
817 C...If shower failed then restore to situation before shower.
818  IF(mint(51).GE.1) THEN
819  n=nsav
820  DO 360 i=1,4
821  DO 350 j=1,5
822  k(n-4+i,j)=ksav(i,j)
823  p(n-4+i,j)=psav(i,j)
824  350 CONTINUE
825  360 CONTINUE
826  ipu1=n-3
827  ipu2=n-2
828  vint(141)=vint(41)
829  vint(142)=vint(42)
830  ENDIF
831  ENDIF
832 
833 C...Keep track of loose colour ends and information on scattering.
834  370 imi(1,mint(31),1)=ipu1
835  imi(2,mint(31),1)=ipu2
836  imi(1,mint(31),2)=0
837  imi(2,mint(31),2)=0
838  xmi(1,mint(31))=vint(141)
839  xmi(2,mint(31))=vint(142)
840  pt2mi(mint(31))=vint(54)
841  imisep(mint(31))=n
842 
843 C...Decide whether quarks in last scattering were valence, companion or
844 C...sea.
845  DO 430 js=1,2
846  kfbeam=mint(10+js)
847  kfsbm=isign(1,mint(10+js))
848  ifl=k(imi(js,mint(31),1),2)
849  imi(js,mint(31),2)=0
850  IF (iabs(ifl).GT.6) goto 430
851 
852 C...Get PDFs at X and Q2 of the parton shower initiator for the
853 C...last scattering. At this point VINT(143:144) do not yet
854 C...include the scattered x values VINT(141:142).
855  x=vint(140+js)/vint(142+js)
856  IF(mstp(84).GE.1.AND.mstp(61).GE.1) THEN
857  q2=parp(62)**2
858  ELSE
859  q2=vint(54)
860  ENDIF
861 C...Note: XPSVC = x*pdf.
862  mint(30)=js
863  CALL pypdfu(kfbeam,x,q2,xpq)
864  sea=xpsvc(ifl,-1)
865  val=xpsvc(ifl,0)
866  cmp=0d0
867  DO 380 ivc=1,nvc(js,ifl)
868  cmp=cmp+xpsvc(ifl,ivc)
869  380 CONTINUE
870 
871 C...Decide (Extra factor x cancels in the dvision).
872  rvcs=pyr(0)*(sea+val+cmp)
873  ivnow=1
874  390 IF (rvcs.LE.val.AND.ivnow.GE.1) THEN
875 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
876  ivnow=0
877  IF(kfival(js,1).EQ.ifl) ivnow=ivnow+1
878  IF(kfival(js,2).EQ.ifl) ivnow=ivnow+1
879  IF(kfival(js,3).EQ.ifl) ivnow=ivnow+1
880  IF(kfival(js,1).EQ.0) THEN
881  IF(kfbeam.EQ.111.AND.iabs(ifl).LE.2) ivnow=1
882  IF(kfbeam.EQ.22.AND.iabs(ifl).LE.5) ivnow=1
883  IF((kfbeam.EQ.130.OR.kfbeam.EQ.310).AND.
884  & (iabs(ifl).EQ.1.OR.iabs(ifl).EQ.3)) ivnow=1
885  ELSE
886  DO 400 i1=1,nmi(js)
887  IF (k(imi(js,i1,1),2).EQ.ifl.AND.imi(js,i1,2).EQ.0)
888  & ivnow=ivnow-1
889  400 CONTINUE
890  ENDIF
891  IF(ivnow.EQ.0) goto 390
892 C...Mark valence.
893  imi(js,mint(31),2)=0
894 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
895  IF(kfival(js,1).EQ.0) THEN
896  IF(kfbeam.EQ.111.OR.kfbeam.EQ.22) THEN
897  kfival(js,1)=ifl
898  kfival(js,2)=-ifl
899  ELSEIF(kfbeam.EQ.130.OR.kfbeam.EQ.310) THEN
900  kfival(js,1)=ifl
901  IF(iabs(ifl).EQ.1) kfival(js,2)=isign(3,-ifl)
902  IF(iabs(ifl).NE.1) kfival(js,2)=isign(1,-ifl)
903  ENDIF
904  ENDIF
905 
906  ELSEIF (rvcs.LE.val+sea.OR.nvc(js,ifl).EQ.0) THEN
907 C...If sea, add opposite sign companion parton. Store X and I.
908  nvc(js,-ifl)=nvc(js,-ifl)+1
909  xassoc(js,-ifl,nvc(js,-ifl))=x
910 C...Set pointer to companion
911  imi(js,mint(31),2)=-nvc(js,-ifl)
912  ELSE
913 C...If companion, decide which one.
914  cmpsum=val+sea
915  isel=0
916  410 isel=isel+1
917  cmpsum=cmpsum+xpsvc(ifl,isel)
918  IF (rvcs.GT.cmpsum.AND.isel.LT.nvc(js,ifl)) goto 410
919 C...Find original sea (anti-)quark:
920  iassoc=0
921  DO 420 i1=1,nmi(js)
922  IF (k(imi(js,i1,1),2).NE.-ifl) goto 420
923  IF (-imi(js,i1,2).EQ.isel) THEN
924  imi(js,mint(31),2)=imi(js,i1,1)
925  imi(js,i1,2)=imi(js,mint(31),1)
926  ENDIF
927  420 CONTINUE
928 C...Change X to what associated companion had, so that the correct
929 C...amount of momentum can be subtracted from the companion sum below.
930  x=xassoc(js,ifl,isel)
931 C...Mark companion read.
932  xassoc(js,ifl,isel)=0d0
933  ENDIF
934  430 CONTINUE
935 
936 C...Global statistics.
937  mint(351)=mint(351)+1
938  vint(351)=vint(351)+pt
939  IF (mint(351).EQ.1) vint(356)=pt
940 
941 C...Update remaining energy and other counters.
942  IF(n.GT.mstu(4)-mstu(32)-10) THEN
943  CALL pyerrm(11,'(PYMIGN:) no more memory left in PYJETS')
944  mint(51)=1
945  RETURN
946  ENDIF
947  nmi(1)=nmi(1)+1
948  nmi(2)=nmi(2)+1
949  vint(151)=vint(151)+vint(41)
950  vint(152)=vint(152)+vint(42)
951  vint(143)=vint(143)-vint(141)
952  vint(144)=vint(144)-vint(142)
953 
954 C...Iterate, with more interactions allowed.
955  IF(mint(31).LT.240) goto 240
956  440 CONTINUE
957 
958 C...Restore saved quantities for hardest interaction.
959  mint(1)=isubsv
960  mint(13)=m13sv
961  mint(14)=m14sv
962  mint(15)=m15sv
963  mint(16)=m16sv
964  mint(21)=m21sv
965  mint(22)=m22sv
966  DO 450 j=11,80
967  vint(j)=vintsv(j)
968  450 CONTINUE
969  vint(141)=v141sv
970  vint(142)=v142sv
971 
972  ENDIF
973 
974 C...Format statements for printout.
975  5000 FORMAT(/1x,'****** PYMIGN: initialization of multiple inter',
976  &'actions for MSTP(82) =',i2,' ******')
977  5100 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,
978  &d9.2,' mb: rejected')
979  5200 FORMAT(8x,'pT0 =',f5.2,' GeV gives sigma(parton-parton) =',1p,
980  &d9.2,' mb: accepted')
981 
982  RETURN
983  END