Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymirm.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymirm.f
1 
2 C*********************************************************************
3 
4 C...PYMIRM
5 C...Picks primordial kT and shares longitudinal momentum among
6 C...beam remnants.
7 
8  SUBROUTINE pymirm
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...The event record
15  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16 C...Parameters
17  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
18  common/pypars/mstp(200),parp(200),msti(200),pari(200)
19  common/pyint1/mint(400),vint(400)
20 C...The common block of colour tags.
21  common/pyctag/nct,mct(4000,2)
22 C...The common block of dangling ends
23  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
24  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
25  & xmi(2,240),pt2mi(240),imisep(0:240)
26  SAVE /pyjets/,/pydat1/,/pypars/,/pyint1/,/pyintm/,/pyctag/
27 C...Local variables
28  dimension w(0:2,0:2),vb(3),nnxt(2),ivalq(2),icomq(2)
29 C...W(I,J)| J=0 | 1 | 2 |
30 C... I=0 | Wrem**2 | W+ | W- |
31 C... 1 | W1**2 | W1+ | W1- |
32 C... 2 | W2**2 | W2+ | W2- |
33 C...4-product
34  four(i,j)=p(i,4)*p(j,4)-p(i,1)*p(j,1)-p(i,2)*p(j,2)-p(i,3)*p(j,3)
35 C...Tentative parametrization of <kT> as a function of Q.
36  sigpt(q)=max(parj(21),2.1d0*q/(7d0+q))
37 C SIGPT(Q)=MAX(0.36D0,4D0*SQRT(Q)/(10D0+SQRT(Q))
38 C SIGPT(Q)=MAX(PARJ(21),3D0*SQRT(Q)/(5D0+SQRT(Q))
39  getpt(q,sigma)=min(sigma*sqrt(-log(pyr(0))),parp(93))
40 C...Lambda kinematic function.
41  flam(a,b,c)=a**2+b**2+c**2-2d0*(a*b+b*c+c*a)
42 
43 C...Beginning and end of beam remnant partons
44  nout=mint(53)
45  isub=mint(1)
46 
47 C...Loopback point if kinematic choices gives impossible configuration.
48  ntry=0
49  100 ntry=ntry+1
50 
51 C...Assign kT values on each side separately.
52  DO 180 js=1,2
53 
54 C...First zero all kT on this side. Skip if no kT to generate.
55  DO 110 im=1,nmi(js)
56  p(imi(js,im,1),1)=0d0
57  p(imi(js,im,1),2)=0d0
58  110 CONTINUE
59  IF(mstp(91).LE.0) goto 180
60 
61 C...Now assign kT to each (non-collapsed) parton in IMI.
62  DO 170 im=1,nmi(js)
63  i=imi(js,im,1)
64 C...Select kT according to truncated gaussian or 1/kt6 tails.
65 C...For first interaction, either use rms width = PARP(91) or fitted.
66  IF (im.EQ.1) THEN
67  sigma=parp(91)
68  IF (mstp(91).GE.11.AND.mstp(91).LE.20) THEN
69  q=sqrt(pt2mi(im))
70  sigma=sigpt(q)
71  ENDIF
72  ELSE
73 C...For subsequent interactions and BR partons use fragmentation width.
74  sigma=parj(21)
75  ENDIF
76  phi=paru(2)*pyr(0)
77  pt=0d0
78  IF(ntry.LE.100) THEN
79  111 IF (mstp(91).EQ.1.OR.mstp(91).EQ.11) THEN
80  pt=getpt(q,sigma)
81  ptx=pt*cos(phi)
82  pty=pt*sin(phi)
83  ELSEIF (mstp(91).EQ.2) THEN
84  CALL pyerrm(11,'(PYMIRM:) Sorry, MSTP(91)=2 not '//
85  & 'available, using MSTP(91)=1.')
86  CALL pygive('MSTP(91)=1')
87  goto 111
88  ELSEIF(mstp(91).EQ.3.OR.mstp(91).EQ.13) THEN
89 C...Use distribution with kt**6 tails, rms width = PARP(91).
90  eps=sqrt(3d0/2d0)*sigma
91 C...Generate PTX and PTY separately, each propto 1/KT**6
92  DO 119 ixy=1,2
93 C...Decide which interval to try
94  112 p12=1d0/(1d0+27d0/40d0*sigma**6/eps**6)
95  IF (pyr(0).LT.p12) THEN
96 C...Use flat approx with accept/reject up to EPS.
97  pt=pyr(0)*eps
98  wt=(3d0/2d0*sigma**2/(pt**2+3d0/2d0*sigma**2))**3
99  IF (pyr(0).GT.wt) goto 112
100  ELSE
101 C...Above EPS, use 1/kt**6 approx with accept/reject.
102  pt=eps/(pyr(0)**(1d0/5d0))
103  wt=pt**6/(pt**2+3d0/2d0*sigma**2)**3
104  IF (pyr(0).GT.wt) goto 112
105  ENDIF
106  msign=1
107  IF (pyr(0).GT.0.5d0) msign=-1
108  IF (ixy.EQ.1) ptx=msign*pt
109  IF (ixy.EQ.2) pty=msign*pt
110  119 CONTINUE
111  ELSEIF (mstp(91).EQ.4.OR.mstp(91).EQ.14) THEN
112  ptx=sigma*(sqrt(6d0)*pyr(0)-sqrt(3d0/2d0))
113  pty=sigma*(sqrt(6d0)*pyr(0)-sqrt(3d0/2d0))
114  ENDIF
115 C...Adjust final PT. Impose upper cutoff, or zero for soft evts.
116  pt=sqrt(ptx**2+pty**2)
117  wt=1d0
118  IF (pt.GT.parp(93)) wt=sqrt(parp(93)/pt)
119  IF(isub.EQ.95.AND.im.EQ.1) wt=0d0
120  ptx=ptx*wt
121  pty=pty*wt
122  pt=sqrt(ptx**2+pty**2)
123  ENDIF
124 
125  p(i,1)=p(i,1)+ptx
126  p(i,2)=p(i,2)+pty
127 
128 C...Compensation kicks, with varying degree of local anticorrelations.
129  mcorr=mstp(90)
130  IF (mcorr.EQ.0.OR.isub.EQ.95) THEN
131  ptcx=-ptx/(nmi(js)-1)
132  ptcy=-pty/(nmi(js)-1)
133  IF(isub.EQ.95) THEN
134  ptcx=-ptx/(nmi(js)-2)
135  ptcy=-pty/(nmi(js)-2)
136  ENDIF
137  DO 120 imc=1,nmi(js)
138  IF (imc.EQ.im) goto 120
139  IF(isub.EQ.95.AND.imc.EQ.1) goto 120
140  p(imi(js,imc,1),1)=p(imi(js,imc,1),1)+ptcx
141  p(imi(js,imc,1),2)=p(imi(js,imc,1),2)+ptcy
142  120 CONTINUE
143  ELSEIF (mcorr.GE.1) THEN
144  DO 140 msid=4,5
145  nnxt(msid-3)=0
146 C...Count up # of neighbours on either side
147  imo=i
148  130 imo=k(imo,msid)/mstu(5)
149  IF (imo.EQ.0) goto 140
150  nnxt(msid-3)=nnxt(msid-3)+1
151 C...Stop at quarks and junctions
152  IF (mcorr.EQ.1.AND.k(imo,2).EQ.21) goto 130
153  140 CONTINUE
154 C...How should compensation be shared when unequal numbers on the
155 C...two sides? 50/50 regardless? N1:N2? Assume latter for now.
156  nsum=nnxt(1)+nnxt(2)
157  t1=0
158  DO 160 msid=4,5
159 C...Total momentum to be compensated on this side
160  IF (nnxt(msid-3).EQ.0) goto 160
161  ptcx=-(nnxt(msid-3)*ptx)/nsum
162  ptcy=-(nnxt(msid-3)*pty)/nsum
163 C...RS: compensation supression factor as we go out from parton I.
164 C...Hardcoded behaviour RS=0.5, i.e. 1/2**n falloff,
165 C...since (for now) MSTP(90) provides enough variability.
166  rs=0.5d0
167  fac=(1d0-rs)/(rs*(1-rs**nnxt(msid-3)))
168  imo=i
169  150 ida=imo
170  imo=k(imo,msid)/mstu(5)
171  IF (imo.EQ.0) goto 160
172  fac=fac*rs
173  IF (k(imo,2).NE.88) THEN
174  p(imo,1)=p(imo,1)+fac*ptcx
175  p(imo,2)=p(imo,2)+fac*ptcy
176  IF (mcorr.EQ.1.AND.k(imo,2).EQ.21) goto 150
177 C...If we reach junction, divide out the kT that would have been
178 C...assigned to the junction on each of its other legs.
179  ELSE
180  l1=mod(k(imo,4),mstu(5))
181  l2=k(imo,5)/mstu(5)
182  l3=mod(k(imo,5),mstu(5))
183  p(l1,1)=p(l1,1)+0.5d0*fac*ptcx
184  p(l1,2)=p(l1,2)+0.5d0*fac*ptcy
185  p(l2,1)=p(l2,1)+0.5d0*fac*ptcx
186  p(l2,2)=p(l2,2)+0.5d0*fac*ptcy
187  p(l3,1)=p(l3,1)+0.5d0*fac*ptcx
188  p(l3,2)=p(l3,2)+0.5d0*fac*ptcy
189  p(ida,1)=p(ida,1)-0.5d0*fac*ptcx
190  p(ida,2)=p(ida,2)-0.5d0*fac*ptcy
191  ENDIF
192 
193  160 CONTINUE
194  ENDIF
195  170 CONTINUE
196 C...End assignment of kT values to initiators and remnants.
197  180 CONTINUE
198 
199 C...Check kinematics constraints for non-BR partons.
200  DO 190 im=1,mint(31)
201  shat=xmi(1,im)*xmi(2,im)*vint(2)
202  pt1=sqrt(p(imi(1,im,1),1)**2+p(imi(1,im,1),2)**2)
203  pt2=sqrt(p(imi(2,im,1),1)**2+p(imi(2,im,1),2)**2)
204  pt1pt2=p(imi(1,im,1),1)*p(imi(2,im,1),1)
205  & +p(imi(1,im,1),2)*p(imi(2,im,1),2)
206  IF (shat.LT.2d0*(pt1*pt2-pt1pt2).AND.ntry.LE.100) THEN
207  IF(ntry.GE.100) THEN
208 C...Kill this event and start another.
209  CALL pyerrm(11,
210  & '(PYMIRM:) No consistent (x,kT) sets found')
211  mint(51)=1
212  RETURN
213  ENDIF
214  goto 100
215  ENDIF
216  190 CONTINUE
217 
218 C...Calculate W+ and W- available for combined remnant system.
219  w(0,1)=vint(1)
220  w(0,2)=vint(1)
221  DO 200 im=1,mint(31)
222  pt2 = (p(imi(1,im,1),1)+p(imi(2,im,1),1))**2
223  & +(p(imi(1,im,1),2)+p(imi(2,im,1),2))**2
224  st=xmi(1,im)*xmi(2,im)*vint(2)+pt2
225  w(0,1)=w(0,1)-sqrt(xmi(1,im)/xmi(2,im)*st)
226  w(0,2)=w(0,2)-sqrt(xmi(2,im)/xmi(1,im)*st)
227  200 CONTINUE
228 C...Also store Wrem**2 = W+ * W-
229  w(0,0)=w(0,1)*w(0,2)
230 
231  IF (w(0,0).LT.0d0.AND.ntry.LE.100) THEN
232  IF(ntry.GE.100) THEN
233 C...Kill this event and start another.
234  CALL pyerrm(11,
235  & '(PYMIRM:) Negative beam remnant mass squared unavoidable')
236  mint(51)=1
237  RETURN
238  ENDIF
239  goto 100
240  ENDIF
241 
242 C...Assign unscaled x values to partons/hadrons in each of the
243 C...beam remnants and calculate unscaled W+ and W- from them.
244  ntryx=0
245  210 ntryx=ntryx+1
246  DO 280 js=1,2
247  w(js,1)=0d0
248  w(js,2)=0d0
249  DO 270 im=mint(31)+1,nmi(js)
250  i=imi(js,im,1)
251  kf=k(i,2)
252  kfa=iabs(kf)
253  icomp=imi(js,im,2)
254 
255 C...Skip collapsed gluons and junctions. Reset.
256  IF (kfa.EQ.21.AND.k(i,1).EQ.14) goto 270
257  IF (kfa.EQ.88) goto 270
258  x=0d0
259  ivalq(1)=0
260  ivalq(2)=0
261  icomq(1)=0
262  icomq(2)=0
263 
264 C...If gluon then only beam remnant, so takes all.
265  IF(kfa.EQ.21) THEN
266  x=1d0
267 C...If valence quark then use parametrized valence distribution.
268  ELSEIF(kfa.LE.6.AND.icomp.EQ.0) THEN
269  ivalq(1)=kf
270 C...If companion quark then derive from companion x.
271  ELSEIF(kfa.LE.6) THEN
272  icomq(1)=icomp
273 C...If valence diquark then use two parametrized valence distributions.
274  ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
275  & icomp.EQ.0) THEN
276  ivalq(1)=isign(kfa/1000,kf)
277  ivalq(2)=isign(mod(kfa/100,10),kf)
278 C...If valence+sea diquark then combine valence + companion choices.
279  ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0.AND.
280  & icomp.LT.mstu(5)) THEN
281  IF(kfa/1000.EQ.iabs(k(icomp,2))) THEN
282  ivalq(1)=isign(mod(kfa/100,10),kf)
283  ELSE
284  ivalq(1)=isign(kfa/1000,kf)
285  ENDIF
286  icomq(1)=icomp
287 C...Extra code: workaround for diquark made out of two sea
288 C...quarks, but where not (yet) ICOMP > MSTU(5).
289  DO 220 im1=1,mint(31)
290  IF(imi(js,im1,2).EQ.i.AND.imi(js,im1,1).NE.icomp) THEN
291  icomq(2)=imi(js,im1,1)
292  ivalq(1)=0
293  ENDIF
294  220 CONTINUE
295 C...If sea diquark then sum of two derived from companion x.
296  ELSEIF(kfa.GT.1000.AND.mod(kfa/10,10).EQ.0) THEN
297  icomq(1)=mod(icomp,mstu(5))
298  icomq(2)=icomp/mstu(5)
299 C...If meson or baryon then use fragmentation function.
300 C...Somewhat arbitrary split into old and new flavour, but OK normally.
301  ELSE
302  kfl3=mod(kfa/10,10)
303  IF(mod(kfa/1000,10).EQ.0) THEN
304  kfl1=mod(kfa/100,10)
305  ELSE
306  kfl1=mod(kfa,10000)-10*kfl3-1
307  IF(mod(kfa/1000,10).EQ.mod(kfa/100,10).AND.
308  & mod(kfa,10).EQ.2) kfl1=kfl1+2
309  ENDIF
310  pr=p(i,5)**2+p(i,1)**2+p(i,2)**2
311  CALL pyzdis(kfl1,kfl3,pr,x)
312  ENDIF
313 
314  DO 260 iq=1,2
315 C...Calculation of x of valence quark: assume form (1-x)^a/sqrt(x),
316 C...where a=3.5 for u in proton, =2 for d in proton and =0.8 for meson.
317 C...In other baryons combine u and d from proton appropriately.
318  IF(ivalq(iq).NE.0) THEN
319  nval=0
320  IF(kfival(js,1).EQ.ivalq(iq)) nval=nval+1
321  IF(kfival(js,2).EQ.ivalq(iq)) nval=nval+1
322  IF(kfival(js,3).EQ.ivalq(iq)) nval=nval+1
323 C...Meson.
324  IF(kfival(js,3).EQ.0) THEN
325  mdu=0
326 C...Baryon with three identical quarks: mix u and d forms.
327  ELSEIF(nval.EQ.3) THEN
328  mdu=int(pyr(0)+5d0/3d0)
329 C...Baryon, one of two identical quarks: u form.
330  ELSEIF(nval.EQ.2) THEN
331  mdu=2
332 C...Baryon with two identical quarks, but not the one picked: d form.
333  ELSEIF(kfival(js,1).EQ.kfival(js,2).OR.kfival(js,2).EQ.
334  & kfival(js,3).OR.kfival(js,1).EQ.kfival(js,3)) THEN
335  mdu=1
336 C...Baryon with three nonidentical quarks: mix u and d forms.
337  ELSE
338  mdu=int(pyr(0)+5d0/3d0)
339  ENDIF
340  xpow=0.8d0
341  IF(mdu.EQ.1) xpow=3.5d0
342  IF(mdu.EQ.2) xpow=2d0
343  230 xx=pyr(0)**2
344  IF((1d0-xx)**xpow.LT.pyr(0)) goto 230
345  x=x+xx
346  ENDIF
347 
348 C...Calculation of x of companion quark.
349  IF(icomq(iq).NE.0) THEN
350  xcomp=1d-4
351  DO 240 im1=1,mint(31)
352  IF(imi(js,im1,1).EQ.icomq(iq)) xcomp=xmi(js,im1)
353  240 CONTINUE
354  npow=max(0,min(4,mstp(87)))
355  250 xx=xcomp*(1d0/(1d0-pyr(0)*(1d0-xcomp))-1d0)
356  corr=((1d0-xcomp-xx)/(1d0-xcomp))**npow*
357  & (xcomp**2+xx**2)/(xcomp+xx)**2
358  IF(corr.LT.pyr(0)) goto 250
359  x=x+xx
360  ENDIF
361  260 CONTINUE
362 
363 C...Optionally enchance x of composite systems (e.g. diquarks)
364  IF (kfa.GT.100) x=parp(79)*x
365 
366 C...Store x. Also calculate light cone energies of each system.
367  xmi(js,im)=x
368  w(js,js)=w(js,js)+x
369  w(js,3-js)=w(js,3-js)+(p(i,5)**2+p(i,1)**2+p(i,2)**2)/x
370  270 CONTINUE
371  w(js,js)=w(js,js)*w(0,js)
372  w(js,3-js)=w(js,3-js)/w(0,js)
373  w(js,0)=w(js,1)*w(js,2)
374  280 CONTINUE
375 
376 C...Check W1 W2 < Wrem (can be done before rescaling, since W
377 C...insensitive to global rescalings of the BR x values).
378  IF (sqrt(w(1,0))+sqrt(w(2,0)).GT.sqrt(w(0,0)).AND.ntryx.LE.100)
379  & THEN
380  goto 210
381  ELSEIF (ntryx.GT.100.AND.ntry.LE.100) THEN
382  goto 100
383  ELSEIF (ntryx.GT.100) THEN
384  CALL pyerrm(11,'(PYMIRM:) No consistent (x,kT) sets found')
385  mint(57)=mint(57)+1
386  mint(51)=1
387  RETURN
388  ENDIF
389 
390 C...Compute x rescaling factors
391  comtrm=w(0,0)+sqrt(flam(w(0,0),w(1,0),w(2,0)))
392  r1=(comtrm+w(1,0)-w(2,0))/(2d0*w(1,1)*w(0,2))
393  r2=(comtrm+w(2,0)-w(1,0))/(2d0*w(2,2)*w(0,1))
394 
395  IF (r1.LT.0.OR.r2.LT.0) THEN
396  CALL pyerrm(19,'(PYMIRM:) negative rescaling factors !')
397  mint(57)=mint(57)+1
398  mint(51)=1
399  ENDIF
400 
401 C...Rescale W(1,*) and W(2,*) (not really necessary, but consistent).
402  w(1,1)=w(1,1)*r1
403  w(1,2)=w(1,2)/r1
404  w(2,1)=w(2,1)/r2
405  w(2,2)=w(2,2)*r2
406 
407 C...Rescale BR x values.
408  DO 290 im=mint(31)+1,max(nmi(1),nmi(2))
409  xmi(1,im)=xmi(1,im)*r1
410  xmi(2,im)=xmi(2,im)*r2
411  290 CONTINUE
412 
413 C...Now we have a consistent set of x and kT values.
414 C...First set up the initiators and their daughters correctly.
415  DO 300 im=1,mint(31)
416  i1=imi(1,im,1)
417  i2=imi(2,im,1)
418  st=xmi(1,im)*xmi(2,im)*vint(2)+(p(i1,1)+p(i2,1))**2+
419  & (p(i1,2)+p(i2,2))**2
420  pt12=p(i1,1)**2+p(i1,2)**2
421  pt22=p(i2,1)**2+p(i2,2)**2
422 C...p_z
423  p(i1,3)=sqrt(flam(st,pt12,pt22)/(4d0*st))
424  p(i2,3)=-p(i1,3)
425 C...Energies (masses should be zero at this stage)
426  p(i1,4)=sqrt(pt12+p(i1,3)**2)
427  p(i2,4)=sqrt(pt22+p(i2,3)**2)
428 
429 C...Transverse 12 system initiator velocity:
430  vb(1)=(p(i1,1)+p(i2,1))/sqrt(st)
431  vb(2)=(p(i1,2)+p(i2,2))/sqrt(st)
432 C...Boost to overall initiator system rest frame
433  CALL pyrobo(i1,i1,0d0,0d0,-vb(1),-vb(2),0d0)
434  CALL pyrobo(i2,i2,0d0,0d0,-vb(1),-vb(2),0d0)
435 C...Compute phi,theta coordinates of I1 and rotate z axis.
436  phi=pyangl(p(i1,1),p(i1,2))
437  the=pyangl(p(i1,3),sqrt(p(i1,1)**2+p(i1,2)**2))
438  CALL pyrobo(i1,i1,0d0,-phi,0d0,0d0,0d0)
439  CALL pyrobo(i2,i2,0d0,-phi,0d0,0d0,0d0)
440  CALL pyrobo(i1,i1,-the,0d0,0d0,0d0,0d0)
441  CALL pyrobo(i2,i2,-the,0d0,0d0,0d0,0d0)
442 
443 C...Now boost initiators + daughters back to LAB system
444 C...(also update documentation lines for MI = 1.)
445  vb(3)=(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
446  imin=imisep(im-1)+1
447  IF (im.EQ.1) imin=mint(83)+5
448  imax=imisep(im)
449  CALL pyrobo(imin,imax,the,phi,vb(1),vb(2),0d0)
450  CALL pyrobo(imin,imax,0d0,0d0,0d0,0d0,vb(3))
451 
452  300 CONTINUE
453 
454 
455 C...For the beam remnant partons/hadrons, we only need to set pz and E.
456  DO 320 js=1,2
457  DO 310 im=mint(31)+1,nmi(js)
458  i=imi(js,im,1)
459 C...Skip collapsed gluons and junctions.
460  IF (k(i,2).EQ.21.AND.k(i,1).EQ.14) goto 310
461  IF (kfa.EQ.88) goto 310
462  rmt2=p(i,5)**2+p(i,1)**2+p(i,2)**2
463  p(i,4)=0.5d0*(xmi(js,im)*w(0,js)+rmt2/(xmi(js,im)*w(0,js)))
464  p(i,3)=0.5d0*(xmi(js,im)*w(0,js)-rmt2/(xmi(js,im)*w(0,js)))
465  IF (js.EQ.2) p(i,3)=-p(i,3)
466  310 CONTINUE
467  320 CONTINUE
468 
469 
470 C...Documentation lines
471  DO 340 js=1,2
472  in=mint(83)+js+2
473  io=imi(js,1,1)
474  k(in,1)=21
475  k(in,2)=k(io,2)
476  k(in,3)=mint(83)+js
477  k(in,4)=0
478  k(in,5)=0
479  DO 330 j=1,5
480  p(in,j)=p(io,j)
481  v(in,j)=v(io,j)
482  330 CONTINUE
483  mct(in,1)=mct(io,1)
484  mct(in,2)=mct(io,2)
485  340 CONTINUE
486 
487 C...Final state colour reconnections.
488  IF (mstp(95).NE.1.OR.mint(31).LE.1) goto 380
489 
490 C...Number of colour tags for which a recoupling will be tried.
491  ntot=nct
492 C...Number of recouplings to try
493  mint(34)=0
494  nrecp=0
495  niter=0
496  350 nrecp=mint(34)
497  niter=niter+1
498  iiter=0
499  360 iiter=iiter+1
500  IF (iiter.LE.parp(78)*ntot) THEN
501 C...Select two colour tags at random
502 C...NB: jj strings do not have colour tags assigned to them,
503 C...thus they are as yet not affected by anything done here.
504  jct=pyr(0)*nct+1
505  kct=mod(int(jct+pyr(0)*nct),nct)+1
506  ij1=0
507  ij2=0
508  ik1=0
509  ik2=0
510 C...Find final state partons with this (anti)colour
511  DO 370 i=mint(84)+1,n
512  IF (k(i,1).EQ.3) THEN
513  IF (mct(i,1).EQ.jct) ij1=i
514  IF (mct(i,2).EQ.jct) ij2=i
515  IF (mct(i,1).EQ.kct) ik1=i
516  IF (mct(i,2).EQ.kct) ik2=i
517  ENDIF
518  370 CONTINUE
519 C...Only consider recouplings not involving junctions for now.
520  IF (ij1.EQ.0.OR.ij2.EQ.0.OR.ik1.EQ.0.OR.ik2.EQ.0) goto 360
521 
522  rlo=2d0*four(ij1,ij2)*2d0*four(ik1,ik2)
523  rln=2d0*four(ij1,ik2)*2d0*four(ik1,ij2)
524  IF (rln.LT.rlo.AND.mct(ij2,1).NE.kct.AND.mct(ik2,1).NE.jct) THEN
525  mct(ij2,2)=kct
526  mct(ik2,2)=jct
527 C...Count up number of reconnections
528  mint(34)=mint(34)+1
529  ENDIF
530  IF (mint(34).LE.1000) THEN
531  goto 360
532  ELSE
533  CALL pyerrm(4,'(PYMIRM:) caught in infinite loop')
534  goto 380
535  ENDIF
536  ENDIF
537  IF (nrecp.LT.mint(34)) goto 350
538 
539 C...Signal PYPREP to use /PYCTAG/ information rather than K(I,KCS).
540  380 mint(33)=1
541 
542  RETURN
543  END