Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyptmi.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyptmi.f
1 
2 C*********************************************************************
3 
4 C...PYPTMI
5 C...Handles the generation of additional interactions in the new
6 C...multiple interactions framework.
7 C...MODE=-1 : Initalize MI from scratch.
8 C...MODE= 0 : Generate trial interaction. Start at PT2NOW, solve
9 C... Sudakov for PT2, abort if below PT2CUT.
10 C...MODE= 1 : Accept interaction at PT2NOW and store variables.
11 C...MODE= 2 : Decide sea/val/cmp for kicked-out quark at PT2NOW
12 C...PT2NOW : Starting (max) PT2 scale for evolution.
13 C...PT2CUT : Lower limit for evolution.
14 C...PT2 : Result of evolution. Generated PT2 for trial interaction.
15 C...IFAIL : Status return code.
16 C... = 0: All is well.
17 C... < 0: Phase space exhausted, generation to be terminated.
18 C... > 0: Additional interaction vetoed, but continue evolution.
19 
20  SUBROUTINE pyptmi(MODE,PT2NOW,PT2CUT,PT2,IFAIL)
21 C...Double precision and integer declarations.
22  IMPLICIT DOUBLE PRECISION(a-h, o-z)
23  IMPLICIT INTEGER(i-n)
24  INTEGER pyk,pychge,pycomp
25 C...Parameter statement for maximum size of showers.
26  parameter(maxnur=1000)
27 C...Commonblocks.
28  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
29  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
30  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
31  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
32  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
33  common/pypars/mstp(200),parp(200),msti(200),pari(200)
34  common/pyint1/mint(400),vint(400)
35  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
36  common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
37  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
38  common/pyint7/sigt(0:6,0:6,0:5)
39  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
40  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
41  & xmi(2,240),pt2mi(240),imisep(0:240)
42  common/pyismx/mimx,jsmx,kflamx,kflcmx,kfbeam(2),nisgen(2,240),
43  & pt2mx,pt2amx,zmx,rm2cmx,q2bmx,phimx
44  common/pyctag/nct,mct(4000,2)
45 C...Local arrays and saved variables.
46  dimension wdtp(0:400),wdte(0:400,0:5),xpq(-25:25)
47 
48  SAVE /pypart/,/pyjets/,/pydat1/,/pydat2/,/pydat3/,/pypars/,
49  & /pyint1/,/pyint2/,/pyint3/,/pyint5/,/pyint7/,/pyintm/,
50  & /pyismx/,/pyctag/
51  SAVE xt2fac,sigs
52 
53  ifail=0
54 C...Set MI subprocess = QCD 2 -> 2.
55  isub=96
56 
57 C----------------------------------------------------------------------
58 C...MODE=-1: Initialize from scratch
59  IF (mode.EQ.-1) THEN
60 C...Initialize PT2 array.
61  pt2mi(1)=vint(54)
62 C...Initialize list of incoming beams and partons from two sides.
63  DO 110 js=1,2
64  DO 100 mi=1,240
65  imi(js,mi,1)=0
66  imi(js,mi,2)=0
67  100 CONTINUE
68  nmi(js)=1
69  imi(js,1,1)=mint(84)+js
70  imi(js,1,2)=0
71  xmi(js,1)=vint(40+js)
72 C...Rescale x values to fractions of photon energy.
73  IF(mint(18+js).EQ.1) xmi(js,1)=vint(40+js)/vint(154+js)
74 C...Hard reset: hard interaction initiators motherless by definition.
75  k(mint(84)+js,3)=2+js
76  k(mint(84)+js,4)=mod(k(mint(84)+js,4),mstu(5))
77  k(mint(84)+js,5)=mod(k(mint(84)+js,5),mstu(5))
78  110 CONTINUE
79  imisep(0)=mint(84)
80  imisep(1)=n
81  IF (mod(mstp(81),10).GE.1) THEN
82  IF(mstp(82).LE.1) THEN
83  sigrat=xsec(isub,1)/max(1d-10,vint(315)*vint(316)*sigt(0,0
84  & ,5))
85  IF(mint(141).NE.0.OR.mint(142).NE.0) sigrat=sigrat*
86  & vint(317)/(vint(318)*vint(320))
87  xt2fac=sigrat*vint(149)/(1d0-vint(149))
88  ELSE
89  xt2fac=vint(146)*vint(148)*xsec(isub,1)/
90  & max(1d-10,sigt(0,0,5))*vint(149)*(1d0+vint(149))
91  ENDIF
92  ENDIF
93 C...Zero entries relating to scatterings beyond the first.
94  DO 120 mi=2,240
95  imi(1,mi,1)=0
96  imi(2,mi,1)=0
97  imi(1,mi,2)=0
98  imi(2,mi,2)=0
99  imisep(mi)=imisep(1)
100  pt2mi(mi)=0d0
101  xmi(1,mi)=0d0
102  xmi(2,mi)=0d0
103  120 CONTINUE
104 C...Initialize factors for PDF reshaping.
105  DO 140 js=1,2
106  kfbeam(js)=mint(10+js)
107  IF(mint(18+js).EQ.1) kfbeam(js)=22
108  kfabm=iabs(kfbeam(js))
109  kfsbm=isign(1,kfbeam(js))
110 
111 C...Zero flavour content of incoming beam particle.
112  kfival(js,1)=0
113  kfival(js,2)=0
114  kfival(js,3)=0
115 C... Flavour content of baryon.
116  IF(kfabm.GT.1000) THEN
117  kfival(js,1)=kfsbm*mod(kfabm/1000,10)
118  kfival(js,2)=kfsbm*mod(kfabm/100,10)
119  kfival(js,3)=kfsbm*mod(kfabm/10,10)
120 C... Flavour content of pi+-, K+-.
121  ELSEIF(kfabm.EQ.211) THEN
122  kfival(js,1)=kfsbm*2
123  kfival(js,2)=-kfsbm
124  ELSEIF(kfabm.EQ.321) THEN
125  kfival(js,1)=-kfsbm*3
126  kfival(js,2)=kfsbm*2
127 C... Flavour content of pi0, gamma, K0S, K0L not defined yet.
128  ENDIF
129 
130 C...Zero initial valence and companion content.
131  DO 130 ifl=-6,6
132  nvc(js,ifl)=0
133  130 CONTINUE
134  140 CONTINUE
135 C...Set up colour line tags starting from hard interaction initiators.
136  nct=0
137 C...Reset colour tag array and colour processing flags.
138  DO 150 i=imisep(0)+1,n
139  mct(i,1)=0
140  mct(i,2)=0
141  k(i,4)=mod(k(i,4),mstu(5)**2)
142  k(i,5)=mod(k(i,5),mstu(5)**2)
143  150 CONTINUE
144 C... Consider each side in turn.
145  DO 170 js=1,2
146  i1=imi(js,1,1)
147  i2=imi(3-js,1,1)
148  DO 160 jcs=4,5
149  IF (k(i1,2).NE.21.AND.(9-2*jcs).NE.isign(1,k(i1,2)))
150  & goto 160
151  IF (k(i1,jcs)/mstu(5)**2.NE.0) goto 160
152  kcs=jcs
153  CALL pycttr(i1,kcs,i2)
154  IF(mint(51).NE.0) RETURN
155  160 CONTINUE
156  170 CONTINUE
157 
158 C...Range checking for companion quark pdf large-x param.
159  IF (mstp(87).LT.0) THEN
160  CALL pyerrm(19,'(PYPTMI:) MSTP(87) out of range. Forced'//
161  & ' MSTP(87)=0')
162  mstp(87)=0
163  ELSEIF (mstp(87).GT.4) THEN
164  CALL pyerrm(19,'(PYPTMI:) MSTP(87) out of range. Forced'//
165  & ' MSTP(87)=4')
166  mstp(87)=4
167  ENDIF
168 
169 C----------------------------------------------------------------------
170 C...MODE=0: Generate trial interaction. Return codes:
171 C...IFAIL < 0: Phase space exhausted, generation to be terminated.
172 C...IFAIL = 0: Additional interaction generated at PT2.
173 C...IFAIL > 0: Additional interaction vetoed, but continue evolution.
174  ELSEIF (mode.EQ.0) THEN
175 C...Abolute MI max scale = VINT(62)
176  xt2=4d0*min(pt2now,vint(62))/vint(2)
177  180 IF(mstp(82).LE.1) THEN
178  xt2=xt2fac*xt2/(xt2fac-xt2*log(pyr(0)))
179  IF(xt2.LT.vint(149)) ifail=-2
180  ELSE
181  IF(xt2.LE.0.01001d0*vint(149)) THEN
182  ifail=-3
183  ELSE
184  xt2=xt2fac*(xt2+vint(149))/(xt2fac-(xt2+vint(149))*
185  & log(pyr(0)))-vint(149)
186  ENDIF
187  ENDIF
188 C...Also exit if below lower limit or if higher trial branching
189 C...already found.
190  pt2=0.25d0*vint(2)*xt2
191  IF (pt2.LE.pt2cut) ifail=-4
192  IF (pt2.LE.pt2mx) ifail=-5
193  IF (ifail.NE.0) THEN
194  pt2=0d0
195  RETURN
196  ENDIF
197  IF(mstp(82).GE.2) pt2=max(0.25d0*vint(2)*0.01d0*vint(149),pt2)
198  vint(25)=4d0*pt2/vint(2)
199  xt2=vint(25)
200 
201 C...Choose tau and y*. Calculate cos(theta-hat).
202  IF(pyr(0).LE.coef(isub,1)) THEN
203  taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**pyr(0)
204  tau=xt2*(1d0+taut)**2/(4d0*taut)
205  ELSE
206  tau=xt2*(1d0+tan(pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
207  ENDIF
208  vint(21)=tau
209 C...New: require shat > 1.
210  IF(tau*vint(2).LT.1d0) goto 180
211  CALL pyklim(2)
212  ryst=pyr(0)
213  myst=1
214  IF(ryst.GT.coef(isub,8)) myst=2
215  IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
216  CALL pykmap(2,myst,pyr(0))
217  vint(23)=sqrt(max(0d0,1d0-xt2/tau))*(-1)**int(1.5d0+pyr(0))
218 
219 C...Check that x not used up. Accept or reject kinematical variables.
220  x1m=sqrt(tau)*exp(vint(22))
221  x2m=sqrt(tau)*exp(-vint(22))
222  IF(vint(143)-x1m.LT.0.01d0.OR.vint(144)-x2m.LT.0.01d0) goto 180
223  vint(71)=0.5d0*vint(1)*sqrt(xt2)
224  CALL pysigh(nchn,sigs)
225  IF(mint(141).NE.0.OR.mint(142).NE.0) sigs=sigs*vint(320)
226  IF(sigs.LT.xsec(isub,1)*pyr(0)) goto 180
227  IF(mint(141).NE.0.OR.mint(142).NE.0) sigs=sigs/vint(320)
228 
229 C...Save if highest PT so far.
230  IF (pt2.GT.pt2mx) THEN
231  jsmx=0
232  mimx=mint(31)+1
233  pt2mx=pt2
234  ENDIF
235 
236 C----------------------------------------------------------------------
237 C...MODE=1: Generate and save accepted scattering.
238  ELSEIF (mode.EQ.1) THEN
239  pt2=pt2now
240 C...Reset K, P, V, and MCT vectors.
241  DO 200 i=n+1,n+4
242  DO 190 j=1,5
243  k(i,j)=0
244  p(i,j)=0d0
245  v(i,j)=0d0
246  190 CONTINUE
247  mct(i,1)=0
248  mct(i,2)=0
249  200 CONTINUE
250 
251  ntry=0
252 C...Choose flavour of reacting partons (and subprocess).
253  210 ntry=ntry+1
254  IF (ntry.GT.50) THEN
255  CALL pyerrm(9,'(PYPTMI:) Unable to generate additional '
256  & //'interaction. Giving up!')
257  mint(51)=1
258  RETURN
259  ENDIF
260  rsigs=sigs*pyr(0)
261  DO 220 ichn=1,nchn
262  kfl1=isig(ichn,1)
263  kfl2=isig(ichn,2)
264  iconmi=isig(ichn,3)
265  rsigs=rsigs-sigh(ichn)
266  IF(rsigs.LE.0d0) goto 230
267  220 CONTINUE
268 
269 C...Reassign to appropriate process codes.
270  230 isubmi=iconmi/10
271  iconmi=mod(iconmi,10)
272 
273 C...Choose new quark flavour for annihilation graphs
274  IF(isubmi.EQ.12.OR.isubmi.EQ.53) THEN
275  sh=vint(21)*vint(2)
276  CALL pywidt(21,sh,wdtp,wdte)
277  240 rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*pyr(0)
278  DO 250 i=1,mdcy(21,3)
279  kflf=kfdp(i+mdcy(21,2)-1,1)
280  rkfl=rkfl-(wdte(i,1)+wdte(i,2)+wdte(i,4))
281  IF(rkfl.LE.0d0) goto 260
282  250 CONTINUE
283  260 IF(isubmi.EQ.53.AND.iconmi.LE.2) THEN
284  IF(kflf.GE.4) goto 240
285  ELSEIF(isubmi.EQ.53.AND.iconmi.LE.4) THEN
286  kflf=4
287  iconmi=iconmi-2
288  ELSEIF(isubmi.EQ.53) THEN
289  kflf=5
290  iconmi=iconmi-4
291  ENDIF
292  ENDIF
293 
294 C...Final state flavours and colour flow: default values
295  js=1
296  kfl3=kfl1
297  kfl4=kfl2
298  kcc=20
299  kcs=isign(1,kfl1)
300 
301  IF(isubmi.EQ.11) THEN
302 C...f + f' -> f + f' (g exchange); th = (p(f)-p(f))**2
303  kcc=iconmi
304  IF(kfl1*kfl2.LT.0) kcc=kcc+2
305 
306  ELSEIF(isubmi.EQ.12) THEN
307 C...f + fbar -> f' + fbar'; th = (p(f)-p(f'))**2
308  kfl3=isign(kflf,kfl1)
309  kfl4=-kfl3
310  kcc=4
311 
312  ELSEIF(isubmi.EQ.13) THEN
313 C...f + fbar -> g + g; th arbitrary
314  kfl3=21
315  kfl4=21
316  kcc=iconmi+4
317 
318  ELSEIF(isubmi.EQ.28) THEN
319 C...f + g -> f + g; th = (p(f)-p(f))**2
320  IF(kfl1.EQ.21) js=2
321  kcc=iconmi+6
322  IF(kfl1.EQ.21) kcc=kcc+2
323  IF(kfl1.NE.21) kcs=isign(1,kfl1)
324  IF(kfl2.NE.21) kcs=isign(1,kfl2)
325 
326  ELSEIF(isubmi.EQ.53) THEN
327 C...g + g -> f + fbar; th arbitrary
328  kcs=(-1)**int(1.5d0+pyr(0))
329  kfl3=isign(kflf,kcs)
330  kfl4=-kfl3
331  kcc=iconmi+10
332 
333  ELSEIF(isubmi.EQ.68) THEN
334 C...g + g -> g + g; th arbitrary
335  kcc=iconmi+12
336  kcs=(-1)**int(1.5d0+pyr(0))
337  ENDIF
338 
339 C...Check that massive sea quarks have non-zero phase space for g -> Q Q
340  IF (iabs(kfl3).EQ.4.OR.iabs(kfl4).EQ.4.OR.iabs(kfl3).EQ.5
341  & .OR.iabs(kfl4).EQ.5) THEN
342  rmmax2=max(pmas(pycomp(kfl3),1),pmas(pycomp(kfl4),1))**2
343  IF (pt2.LE.1.05*rmmax2) THEN
344  IF (ntry.EQ.1) CALL pyerrm(9,'(PYPTMI:) Heavy quarks'
345  & //' created below threshold. Rejected.')
346  goto 210
347  ENDIF
348  ENDIF
349 
350 C...Store flavours of scattering.
351  mint(13)=kfl1
352  mint(14)=kfl2
353  mint(15)=kfl1
354  mint(16)=kfl2
355  mint(21)=kfl3
356  mint(22)=kfl4
357 
358 C...Set flavours and mothers of scattering partons.
359  k(n+1,1)=14
360  k(n+2,1)=14
361  k(n+3,1)=3
362  k(n+4,1)=3
363  k(n+1,2)=kfl1
364  k(n+2,2)=kfl2
365  k(n+3,2)=kfl3
366  k(n+4,2)=kfl4
367  k(n+1,3)=mint(83)+1
368  k(n+2,3)=mint(83)+2
369  k(n+3,3)=n+1
370  k(n+4,3)=n+2
371 
372 C...Store colour connection indices.
373  DO 270 j=1,2
374  jc=j
375  IF(kcs.EQ.-1) jc=3-j
376  IF(icol(kcc,1,jc).NE.0) k(n+1,j+3)=n+icol(kcc,1,jc)
377  IF(icol(kcc,2,jc).NE.0) k(n+2,j+3)=n+icol(kcc,2,jc)
378  IF(icol(kcc,3,jc).NE.0) k(n+3,j+3)=mstu(5)*(n+icol(kcc,3,jc))
379  IF(icol(kcc,4,jc).NE.0) k(n+4,j+3)=mstu(5)*(n+icol(kcc,4,jc))
380  270 CONTINUE
381 
382 C...Store incoming and outgoing partons in their CM-frame.
383  shr=sqrt(vint(21))*vint(1)
384  p(n+1,3)=0.5d0*shr
385  p(n+1,4)=0.5d0*shr
386  p(n+2,3)=-0.5d0*shr
387  p(n+2,4)=0.5d0*shr
388  p(n+3,5)=pymass(k(n+3,2))
389  p(n+4,5)=pymass(k(n+4,2))
390  IF(p(n+3,5)+p(n+4,5).GE.shr) THEN
391  ifail=1
392  RETURN
393  ENDIF
394  p(n+3,4)=0.5d0*(shr+(p(n+3,5)**2-p(n+4,5)**2)/shr)
395  p(n+3,3)=sqrt(max(0d0,p(n+3,4)**2-p(n+3,5)**2))
396  p(n+4,4)=shr-p(n+3,4)
397  p(n+4,3)=-p(n+3,3)
398 
399 C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4)
400  phi=paru(2)*pyr(0)
401  CALL pyrobo(n+3,n+4,acos(vint(23)),phi,0d0,0d0,0d0)
402 
403 C...Global statistics.
404  mint(351)=mint(351)+1
405  vint(351)=vint(351)+sqrt(p(n+3,1)**2+p(n+3,2)**2)
406  IF (mint(351).EQ.1) vint(356)=sqrt(p(n+3,1)**2+p(n+3,2)**2)
407 
408 C...Keep track of loose colour ends and information on scattering.
409  mint(31)=mint(31)+1
410  mint(36)=mint(31)
411  pt2mi(mint(36))=pt2
412  imisep(mint(31))=n+4
413  DO 280 js=1,2
414  imi(js,mint(31),1)=n+js
415  imi(js,mint(31),2)=0
416  xmi(js,mint(31))=vint(40+js)
417  nmi(js)=nmi(js)+1
418 C...Update cumulative counters
419  vint(142+js)=vint(142+js)-vint(40+js)
420  vint(150+js)=vint(150+js)+vint(40+js)
421  280 CONTINUE
422 
423 C...Add to list of final state partons
424  ipart(npart+1)=n+3
425  ipart(npart+2)=n+4
426  ptpart(npart+1)=sqrt(pt2)
427  ptpart(npart+2)=sqrt(pt2)
428  npart=npart+2
429 
430 C...Initialize ISR
431  nisgen(1,mint(31))=0
432  nisgen(2,mint(31))=0
433 
434 C...Update ER
435  n=n+4
436  IF(n.GT.mstu(4)-mstu(32)-10) THEN
437  CALL pyerrm(11,'(PYMIGN:) no more memory left in PYJETS')
438  mint(51)=1
439  RETURN
440  ENDIF
441 
442 C...Finally, assign colour tags to new partons
443  DO 300 js=1,2
444  i1=imi(js,mint(31),1)
445  i2=imi(3-js,mint(31),1)
446  DO 290 jcs=4,5
447  IF (k(i1,2).NE.21.AND.(9-2*jcs).NE.isign(1,k(i1,2)))
448  & goto 290
449  IF (k(i1,jcs)/mstu(5)**2.NE.0) goto 290
450  kcs=jcs
451  CALL pycttr(i1,kcs,i2)
452  IF(mint(51).NE.0) RETURN
453  290 CONTINUE
454  300 CONTINUE
455 
456 C----------------------------------------------------------------------
457 C...MODE=2: Decide whether quarks in last scattering were valence,
458 C...companion, or sea.
459  ELSEIF (mode.EQ.2) THEN
460  js=mint(30)
461  mi=mint(36)
462  pt2=pt2now
463  kfsbm=isign(1,mint(10+js))
464  ifl=k(imi(js,mi,1),2)
465  imi(js,mi,2)=0
466  IF (iabs(ifl).GE.6) THEN
467  IF (iabs(ifl).EQ.6) THEN
468  CALL pyerrm(29,'(PYPTMI:) top in initial state!')
469  ENDIF
470  RETURN
471  ENDIF
472 C...Get PDFs at X(rescaled) and PT2 of the current initiator.
473 C...(Do not include the parton itself in the X rescaling.)
474  x=xmi(js,mi)
475  xrsc=x/(vint(142+js)+x)
476 C...Note: XPSVC = x*pdf.
477  mint(30)=js
478  CALL pypdfu(kfbeam(js),xrsc,pt2,xpq)
479  sea=xpsvc(ifl,-1)
480  val=xpsvc(ifl,0)
481  cmp=0d0
482  DO 310 ivc=1,nvc(js,ifl)
483  cmp=cmp+xpsvc(ifl,ivc)
484  310 CONTINUE
485 
486 C...Decide (Extra factor x cancels in the dvision).
487  320 rvcs=pyr(0)*(sea+val+cmp)
488  ivnow=1
489  330 IF (rvcs.LE.val.AND.ivnow.GE.1) THEN
490 C...Safety check that valence present; pi0/gamma/K0S/K0L special cases.
491  ivnow=0
492  IF(kfival(js,1).EQ.ifl) ivnow=ivnow+1
493  IF(kfival(js,2).EQ.ifl) ivnow=ivnow+1
494  IF(kfival(js,3).EQ.ifl) ivnow=ivnow+1
495  IF(kfival(js,1).EQ.0) THEN
496  IF(kfbeam(js).EQ.111.AND.iabs(ifl).LE.2) ivnow=1
497  IF(kfbeam(js).EQ.22.AND.iabs(ifl).LE.5) ivnow=1
498  IF((kfbeam(js).EQ.130.OR.kfbeam(js).EQ.310).AND.
499  & (iabs(ifl).EQ.1.OR.iabs(ifl).EQ.3)) ivnow=1
500  ELSE
501 C...Count down valence remaining. Do not count current scattering.
502  DO 340 i1=1,nmi(js)
503  IF (i1.EQ.mint(36)) goto 340
504  IF (k(imi(js,i1,1),2).EQ.ifl.AND.imi(js,i1,2).EQ.0)
505  & ivnow=ivnow-1
506  340 CONTINUE
507  ENDIF
508  IF(ivnow.EQ.0) goto 330
509 C...Mark valence.
510  imi(js,mi,2)=0
511 C...Sets valence content of gamma, pi0, K0S, K0L if not done.
512  IF(kfival(js,1).EQ.0) THEN
513  IF(kfbeam(js).EQ.111.OR.kfbeam(js).EQ.22) THEN
514  kfival(js,1)=ifl
515  kfival(js,2)=-ifl
516  ELSEIF(kfbeam(js).EQ.130.OR.kfbeam(js).EQ.310) THEN
517  kfival(js,1)=ifl
518  IF(iabs(ifl).EQ.1) kfival(js,2)=isign(3,-ifl)
519  IF(iabs(ifl).NE.1) kfival(js,2)=isign(1,-ifl)
520  ENDIF
521  ENDIF
522 
523  ELSEIF (rvcs.LE.val+sea) THEN
524 C...If sea, add opposite sign companion parton. Store X and I.
525  nvc(js,-ifl)=nvc(js,-ifl)+1
526  xassoc(js,-ifl,nvc(js,-ifl))=xmi(js,mi)
527 C...Set pointer to companion
528  imi(js,mi,2)=-nvc(js,-ifl)
529 
530  ELSE
531 C...If companion, decide which one.
532  IF (nvc(js,ifl).EQ.0) THEN
533  cmp=0d0
534  CALL pyerrm(9,'(PYPTMI:) No cmp quark, but pdf != 0!')
535  goto 320
536  ENDIF
537  cmpsum=val+sea
538  isel=0
539  350 isel=isel+1
540  cmpsum=cmpsum+xpsvc(ifl,isel)
541  IF (rvcs.GT.cmpsum.AND.isel.LT.nvc(js,ifl)) goto 350
542 C...Find original sea (anti-)quark. Do not consider current scattering.
543  iassoc=0
544  DO 360 i1=1,nmi(js)
545  IF (i1.EQ.mint(36)) goto 360
546  IF (k(imi(js,i1,1),2).NE.-ifl) goto 360
547  IF (-imi(js,i1,2).EQ.isel) THEN
548  imi(js,mi,2)=imi(js,i1,1)
549  imi(js,i1,2)=imi(js,mi,1)
550  ENDIF
551  360 CONTINUE
552 C...Mark companion "out-kicked".
553  xassoc(js,ifl,isel)=-xassoc(js,ifl,isel)
554  ENDIF
555 
556  ENDIF
557  RETURN
558  END