Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pymihk.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pymihk.f
1 
2 C*********************************************************************
3 
4 C...PYMIHK
5 C...Finds left-behind remnant flavour content and hooks up
6 C...the colour flow between the hard scattering and remnants
7 
8  SUBROUTINE pymihk
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/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
19  common/pypars/mstp(200),parp(200),msti(200),pari(200)
20  common/pyint1/mint(400),vint(400)
21 C...The common block of dangling ends
22  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
23  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
24  & xmi(2,240),pt2mi(240),imisep(0:240)
25  SAVE /pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/,/pyintm/
26 C...Local variables
27  parameter(nersiz=4000)
28  COMMON /pycbls/mco(nersiz,2),ncc,jcco(nersiz,2),jccn(nersiz,2)
29  & ,maccpt
30  COMMON /pyctag/nct,mct(nersiz,2)
31  SAVE /pycbls/,/pyctag/
32  dimension jst(2,3),iv(2,3),idq(3),nvsum(2),nbrtot(2),ng(2)
33  & ,itjunc(2),mout(2),insr(1000,3),istr(6),ymi(240)
34  DATA nerrpr/0/
35  SAVE nerrpr
36  four(i,j)=p(i,4)*p(j,4)-p(i,3)*p(j,3)-p(i,2)*p(j,2)-p(i,1)*p(j,1)
37 
38 C...Set up error checkers
39  iboost=0
40 
41 C...Initialize colour arrays: MCO (Original) and MCT (New)
42  DO 110 i=mint(84)+1,nersiz
43  DO 100 jc=1,2
44  mct(i,jc)=0
45  mco(i,jc)=0
46  100 CONTINUE
47 C...Also zero colour tracing information, if existed.
48  IF (i.LE.n) THEN
49  k(i,4)=mod(k(i,4),mstu(5)**2)
50  k(i,5)=mod(k(i,5),mstu(5)**2)
51  ENDIF
52  110 CONTINUE
53 
54 C...Initialize colour tag collapse arrays:
55 C...JCCO (Original) and JCCN (New).
56  DO 130 mg=mint(84)+1,nersiz
57  DO 120 jc=1,2
58  jcco(mg,jc)=0
59  jccn(mg,jc)=0
60  120 CONTINUE
61  130 CONTINUE
62 
63 C...Zero gluon insertion array
64  DO 150 im=1,1000
65  DO 140 j=1,3
66  insr(im,j)=0
67  140 CONTINUE
68  150 CONTINUE
69 
70 C...Compute hard scattering system rapidities
71  IF (mstp(89).EQ.1) THEN
72  DO 160 im=1,240
73  IF (im.LE.mint(31)) THEN
74  ymi(im)=log(xmi(1,im)/xmi(2,im))
75  ELSE
76 C...Set (unsigned) rapidity = 100 for beam remnant systems.
77  ymi(im)=100d0
78  ENDIF
79  160 CONTINUE
80  ENDIF
81 
82 C...Treat each side separately
83  DO 290 js=1,2
84 
85 C...Initialize side.
86  ng(js)=0
87  jv=0
88  kfs=isign(1,mint(10+js))
89 
90 C...Set valence content of pi0, gamma, K0S, K0L if not yet done.
91  IF(kfival(js,1).EQ.0) THEN
92  IF(mint(10+js).EQ.111) THEN
93  kfival(js,1)=int(1.5d0+pyr(0))
94  kfival(js,2)=-kfival(js,1)
95  ELSEIF(mint(10+js).EQ.22) THEN
96  pyrkf=pyr(0)
97  kfival(js,1)=1
98  IF(pyrkf.GT.0.1d0) kfival(js,1)=2
99  IF(pyrkf.GT.0.5d0) kfival(js,1)=3
100  IF(pyrkf.GT.0.6d0) kfival(js,1)=4
101  kfival(js,2)=-kfival(js,1)
102  ELSEIF(mint(10+js).EQ.130.OR.mint(10+js).EQ.310) THEN
103  IF(pyr(0).GT.0.5d0) THEN
104  kfival(js,1)=1
105  kfival(js,2)=-3
106  ELSE
107  kfival(js,1)=3
108  kfival(js,2)=-1
109  ENDIF
110  ENDIF
111  ENDIF
112 
113 C...Initialize beam remnant sea and valence content flavour by flavour.
114  nvsum(js)=0
115  nbrtot(js)=0
116  DO 210 jfa=1,6
117 C...Count up original number of JFA valence quarks and antiquarks.
118  nvalq=0
119  nvalqb=0
120  nsea=0
121  DO 170 j=1,3
122  IF(kfival(js,j).EQ.jfa) nvalq=nvalq+1
123  IF(kfival(js,j).EQ.-jfa) nvalqb=nvalqb+1
124  170 CONTINUE
125  nvsum(js)=nvsum(js)+nvalq+nvalqb
126 C...Subtract kicked out valence and determine sea from flavour cons.
127  DO 180 im=1,nmi(js)
128  ifl = k(imi(js,im,1),2)
129  ifa = iabs(ifl)
130  ifs = isign(1,ifl)
131  IF (ifl.EQ.jfa.AND.imi(js,im,2).EQ.0) THEN
132 C...Subtract K.O. valence quark from remainder.
133  nvalq=nvalq-1
134  jv=nvsum(js)-nvalq-nvalqb
135  iv(js,jv)=imi(js,im,1)
136  ELSEIF (ifl.EQ.-jfa.AND.imi(js,im,2).EQ.0) THEN
137 C...Subtract K.O. valence antiquark from remainder.
138  nvalqb=nvalqb-1
139  jv=nvsum(js)-nvalq-nvalqb
140  iv(js,jv)=imi(js,im,1)
141  ELSEIF (ifa.EQ.jfa) THEN
142 C...Outside sea without companion: add opposite sea flavour inside.
143  IF (imi(js,im,2).LT.0) nsea=nsea-ifs
144  ENDIF
145  180 CONTINUE
146 C...Check if space left in PYJETS for additional BR flavours
147  nflsum=iabs(nsea)+nvalq+nvalqb
148  nbrtot(js)=nbrtot(js)+nflsum
149  IF (n+nflsum+1.GT.mstu(4)) THEN
150  CALL pyerrm(11,'(PYMIHK:) no more memory left in PYJETS')
151  mint(51)=1
152  RETURN
153  ENDIF
154 C...Add required val+sea content to beam remnant.
155  IF (nflsum.GT.0) THEN
156  DO 200 ia=1,nflsum
157 C...Insert beam remnant quark as p.t. symbolic parton in ER.
158  n=n+1
159  DO 190 ix=1,5
160  k(n,ix)=0
161  p(n,ix)=0d0
162  v(n,ix)=0d0
163  190 CONTINUE
164  k(n,1)=3
165  k(n,2)=isign(jfa,nsea)
166  IF (ia.LE.nvalq) k(n,2)=jfa
167  IF (ia.GT.nvalq.AND.ia.LE.nvalq+nvalqb) k(n,2)=-jfa
168  k(n,3)=mint(83)+js
169 C...Also update NMI, IMI, and IV arrays.
170  nmi(js)=nmi(js)+1
171  imi(js,nmi(js),1)=n
172  imi(js,nmi(js),2)=-1
173  IF (ia.LE.nvalq+nvalqb) THEN
174  imi(js,nmi(js),2)=0
175  jv=jv+1
176  iv(js,jv)=imi(js,nmi(js),1)
177  ENDIF
178  200 CONTINUE
179  ENDIF
180  210 CONTINUE
181 
182  im=0
183  220 im=im+1
184  IF (im.LE.nmi(js)) THEN
185  IF (k(imi(js,im,1),2).EQ.21) THEN
186  ng(js)=ng(js)+1
187 C...Add fictitious parent gluons for companion pairs.
188  ELSEIF (imi(js,im,2).NE.0.AND.k(imi(js,im,1),2).GT.0) THEN
189 C...Randomly assign companions to sea quarks which have none.
190  IF (imi(js,im,2).LT.0) THEN
191  imc=pyr(0)*nmi(js)
192  230 imc=mod(imc,nmi(js))+1
193  IF (k(imi(js,imc,1),2).NE.-k(imi(js,im,1),2)) goto 230
194  IF (imi(js,imc,2).GE.0) goto 230
195  imi(js, im,2) = imi(js,imc,1)
196  imi(js,imc,2) = imi(js, im,1)
197  ENDIF
198 C...Add fictitious parent gluon
199  n=n+1
200  DO 240 ix=1,5
201  k(n,ix)=0
202  p(n,ix)=0d0
203  v(n,ix)=0d0
204  240 CONTINUE
205  k(n,1)=14
206  k(n,2)=21
207  k(n,3)=mint(83)+js
208 C...Set gluon (anti-)colour daughter pointers
209  k(n,4)=imi(js, im,1)
210  k(n,5)=imi(js, im,2)
211 C...Set quark (anti-)colour parent pointers
212  k(imi(js, im,2),5)=k(imi(js, im,2),5)+mstu(5)*n
213  k(imi(js, im,1),4)=k(imi(js, im,1),4)+mstu(5)*n
214 C...Add gluon to IMI
215  nmi(js)=nmi(js)+1
216  imi(js,nmi(js),1)=n
217  imi(js,nmi(js),2)=0
218  ENDIF
219  goto 220
220  ENDIF
221 
222 C...If incoming (anti-)baryon, insert inside (anti-)junction.
223 C...Set up initial v-v-j-v configuration. Otherwise set up
224 C...mesonic v-vbar configuration
225  IF (iabs(mint(10+js)).GT.1000) THEN
226 C...Determine junction type (1: B=1 2: B=-1)
227  itjunc(js) = (3-kfs)/2
228 C...Insert junction.
229  n=n+1
230  DO 250 ix=1,5
231  k(n,ix)=0
232  p(n,ix)=0d0
233  v(n,ix)=0d0
234  250 CONTINUE
235 C...Set special junction codes:
236  k(n,1)=42
237  k(n,2)=88
238 C...Set parent to side.
239  k(n,3)=mint(83)+js
240  k(n,4)=itjunc(js)*mstu(5)
241  k(n,5)=0
242 C...Connect valence quarks to junction.
243  mout(js)=0
244  manti=itjunc(js)-1
245 C...Set (anti)colour mother = junction.
246  DO 260 jv=1,3
247  k(iv(js,jv),4+manti)=mod(k(iv(js,jv),4+manti),mstu(5))
248  & +mstu(5)*n
249 C...Keep track of partons adjacent to junction:
250  jst(js,jv)=iv(js,jv)
251  260 CONTINUE
252  ELSE
253 C...Mesons: set up initial q-qbar topology
254  itjunc(js)=0
255  IF (k(iv(js,1),2).GT.0) THEN
256  iq=iv(js,1)
257  iqbar=iv(js,2)
258  ELSE
259  iq=iv(js,2)
260  iqbar=iv(js,1)
261  ENDIF
262  iv(js,3)=0
263  jst(js,1)=iq
264  jst(js,2)=iqbar
265  jst(js,3)=0
266  k(iq,4)=mod(k(iq,4),mstu(5))+mstu(5)*iqbar
267  k(iqbar,5)=mod(k(iqbar,5),mstu(5))+mstu(5)*iq
268 C...Special for mesons. Insert gluon if BR empty.
269  IF (nbrtot(js).EQ.0) THEN
270  n=n+1
271  DO 270 ix=1,5
272  k(n,ix)=0
273  p(n,ix)=0d0
274  v(n,ix)=0d0
275  270 CONTINUE
276  k(n,1)=3
277  k(n,2)=21
278  k(n,3)=mint(83)+js
279  k(n,4)=0
280  k(n,5)=0
281  nbrtot(js)=1
282  ng(js)=ng(js)+1
283 C...Add gluon to IMI
284  nmi(js)=nmi(js)+1
285  imi(js,nmi(js),1)=n
286  imi(js,nmi(js),2)=0
287  ENDIF
288  mout(js)=0
289  ENDIF
290 
291 C...Count up number of valence quarks outside BR.
292  DO 280 jv=1,3
293  IF (jst(js,jv).LE.mint(53).AND.jst(js,jv).GT.0)
294  & mout(js)=mout(js)+1
295  280 CONTINUE
296 
297  290 CONTINUE
298 
299 C...Now both sides have been prepared in an initial vvjv (baryonic) or
300 C...v(g)vbar (mesonic) configuration.
301 
302 C...Create colour line tags starting from initiators.
303  nct=0
304  DO 320 im=1,mint(31)
305 C...Consider each side in turn.
306  DO 310 js=1,2
307  i1=imi(js,im,1)
308  i2=imi(3-js,im,1)
309  DO 300 jcs=4,5
310  IF (k(i1,2).NE.21.AND.(9-2*jcs).NE.isign(1,k(i1,2)))
311  & goto 300
312  IF (k(i1,jcs)/mstu(5)**2.NE.0) goto 300
313 
314  kcs=jcs
315  CALL pycttr(i1,kcs,i2)
316  IF(mint(51).NE.0) RETURN
317 
318  300 CONTINUE
319  310 CONTINUE
320  320 CONTINUE
321 
322  DO 340 js=1,2
323 C...Create colour tags for beam remnant partons.
324  DO 330 im=mint(31)+1,nmi(js)
325  ip=imi(js,im,1)
326  IF (k(ip,2).NE.21) THEN
327  jc=(3-isign(1,k(ip,2)))/2
328  IF (mct(ip,jc).EQ.0) THEN
329  nct=nct+1
330  mct(ip,jc)=nct
331  ENDIF
332  ELSE
333 C...Gluons
334  icd=k(ip,4)
335  iad=k(ip,5)
336  IF (icd.NE.0) THEN
337 C...Fictituous gluons just inherit from their quark daughters.
338  icc=mct(icd,1)
339  iac=mct(iad,2)
340  ELSE
341 C...Real beam remnant gluons get their own colours
342  icc=nct+1
343  iac=nct+2
344  nct=nct+2
345  ENDIF
346  mct(ip,1)=icc
347  mct(ip,2)=iac
348  ENDIF
349  330 CONTINUE
350  340 CONTINUE
351 
352 C...Create colour tags for colour lines which are detached from the
353 C...initial state.
354 
355  DO 360 mqgst=1,2
356  DO 350 i=mint(84)+1,n
357 
358 C...Look for coloured string endpoint, or (later) leftover gluon.
359  IF (k(i,1).NE.3) goto 350
360  kc=pycomp(k(i,2))
361  IF(kc.EQ.0) goto 350
362  kq=kchg(kc,2)
363  IF(kq.EQ.0.OR.(mqgst.EQ.1.AND.kq.EQ.2)) goto 350
364 
365 C...Pick up loose string end with no previous tag.
366  kcs=4
367  IF(kq*isign(1,k(i,2)).LT.0) kcs=5
368  IF(mct(i,kcs-3).NE.0) goto 350
369 
370  CALL pycttr(i,kcs,i)
371  IF(mint(51).NE.0) RETURN
372 
373  350 CONTINUE
374  360 CONTINUE
375 
376 C...Store original colour tags
377  DO 370 i=mint(84)+1,n
378  mco(i,1)=mct(i,1)
379  mco(i,2)=mct(i,2)
380  370 CONTINUE
381 
382 C...Iteratively add gluons to already existing string pieces, enforcing
383 C...various possible orderings, and rejecting insertions that would give
384 C...rise to singlet gluons.
385 C...<kappa tau> normalization.
386  rm0=1.5d0
387  mretry=0
388  parp80=parp(80)
389 
390 C...Set up simplified kinematics.
391 C...Boost hard interaction systems.
392  iboost=iboost+1
393  DO 380 im=1,mint(31)
394  beta=(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
395  CALL pyrobo(imisep(im-1)+1,imisep(im),0d0,0d0,0d0,0d0,beta)
396  380 CONTINUE
397 C...Assign preliminary beam remnant momenta.
398  DO 390 i=mint(53)+1,n
399  js=k(i,3)
400  p(i,1)=0d0
401  p(i,2)=0d0
402  IF (k(i,2).NE.88) THEN
403  p(i,4)=0.5d0*vint(142+js)*vint(1)/max(1,nmi(js)-mint(31))
404  p(i,3)=p(i,4)
405  IF (js.EQ.2) p(i,3)=-p(i,3)
406  ELSE
407 C...Junctions are wildcards for the present.
408  p(i,4)=0d0
409  p(i,3)=0d0
410  ENDIF
411  390 CONTINUE
412 
413 C...Reset colour processing information.
414  400 DO 410 i=mint(84)+1,n
415  k(i,4)=mod(k(i,4),mstu(5)**2)
416  k(i,5)=mod(k(i,5),mstu(5)**2)
417  410 CONTINUE
418 
419  ncc=0
420  DO 430 js=1,2
421 C...If meson, without gluon in BR, collapse q-qbar colour tags:
422  IF (itjunc(js).EQ.0) THEN
423  jc1=mct(jst(js,1),1)
424  jc2=mct(jst(js,2),2)
425  ncc=ncc+1
426  jcco(ncc,1)=max(jc1,jc2)
427  jcco(ncc,2)=min(jc1,jc2)
428 C...Collapse colour tags in event record
429  DO 420 i=mint(84)+1,n
430  IF (mct(i,1).EQ.jcco(ncc,1)) mct(i,1)=jcco(ncc,2)
431  IF (mct(i,2).EQ.jcco(ncc,1)) mct(i,2)=jcco(ncc,2)
432  420 CONTINUE
433  ENDIF
434  430 CONTINUE
435 
436  440 js=1
437  IF (pyr(0).GT.0.5d0.OR.ng(1).EQ.0) js=2
438  IF (ng(js).GT.0) THEN
439  nopt=0
440  rlopt=1d9
441 C...Start at random gluon (optimizes speed for random attachments)
442  nmgl=0
443  imgl=pyr(0)*nmi(js)+1
444  450 imgl=mod(imgl,nmi(js))+1
445  nmgl=nmgl+1
446 C...Only loop through NMI once (with upper limit to save time)
447  IF (nmgl.LE.nmi(js).AND.nopt.LE.3) THEN
448  igl = imi(js,imgl,1)
449 C...If not gluon or if already connected, try next.
450  IF (k(igl,2).NE.21.OR.k(igl,4)/mstu(5).NE.0
451  & .OR.k(igl,5)/mstu(5).NE.0) goto 450
452 C...Now loop through all possible insertions of this gluon.
453  nmp1=0
454  imp1=pyr(0)*nmi(js)+1
455  460 imp1=mod(imp1,nmi(js))+1
456  nmp1=nmp1+1
457  IF (imp1.EQ.imgl) goto 460
458 C...Only loop through NMI once (with upper limit to save time).
459  IF (nmp1.LE.nmi(js).AND.nopt.LE.3) THEN
460  ip1 = imi(js,imp1,1)
461 C...Try both colour mother and colour anti-mother.
462 C...Randomly select which one to try first.
463  nanti=0
464  manti=pyr(0)*2
465  470 manti=mod(manti+1,2)
466  nanti=nanti+1
467  IF (nanti.LE.2) THEN
468  ip2 =mod(k(ip1,4+manti)/mstu(5),mstu(5))
469 C...Reject if no appropriate mother (or if mother is fictitious
470 C...parent gluon.)
471  IF (ip2.LE.0) goto 470
472  IF (k(ip2,2).EQ.21.AND.ip2.GT.mint(53)) goto 470
473 C...Also reject if this link has already been tried.
474  IF (k(ip1,4+manti)/mstu(5)**2.EQ.2) goto 470
475  IF (k(ip2,5-manti)/mstu(5)**2.EQ.2) goto 470
476 C...Set flag to indicate that this link has now been tried for this
477 C...gluon. IP2 may be junction, which has several mothers.
478  k(ip1,4+manti)=k(ip1,4+manti)+2*mstu(5)**2
479  IF (k(ip2,2).NE.88) THEN
480  k(ip2,5-manti)=k(ip2,5-manti)+2*mstu(5)**2
481  ENDIF
482 
483 C...JCG1: Original colour tag of gluon on IP1 side
484 C...JCG2: Original colour tag of gluon on IP2 side
485 C...JCP1: Original colour tag of IP1 on gluon side
486 C...JCP2: Original colour tag of IP2 on gluon side.
487  jcg1=mco(igl,2-manti)
488  jcg2=mco(igl,1+manti)
489  jcp1=mco(ip1,1+manti)
490  jcp2=mco(ip2,2-manti)
491 
492  CALL pymihg(jcp1,jcg1,jcp2,jcg2)
493 C...Reject gluon attachments that give rise to singlet gluons.
494  IF (maccpt.EQ.0) goto 470
495 
496 C...Update colours
497  jcg1=mct(igl,2-manti)
498  jcg2=mct(igl,1+manti)
499  jcp1=mct(ip1,1+manti)
500  jcp2=mct(ip2,2-manti)
501 
502 C...Select whether to accept this insertion
503  IF (mstp(89).EQ.0) THEN
504 C...Random insertions: no measure.
505  rl=1d0
506 C...For random ordering, we want to suppress beam remnant breakups
507 C...already at this point.
508  IF (ip1.GT.mint(53).AND.ip2.GT.mint(53)
509  & .AND.mout(js).NE.0.AND.pyr(0).GT.parp80) THEN
510  nmp1=0
511  nmgl=0
512  goto 470
513  ENDIF
514  ELSEIF (mstp(89).EQ.1) THEN
515 C...Rapidity ordering:
516 C...YGL = Rapidity of gluon.
517  ygl=ymi(imgl)
518 C...If fictitious gluon
519  IF (ygl.EQ.100d0) THEN
520  ygl=(3-2*js)*100d0
521  ida1=mod(k(igl,4),mstu(5))
522  ida2=mod(k(igl,5),mstu(5))
523  DO 480 imt=1,nmi(js)
524 C...Select (arbitrarily) the most central daughter.
525  IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2)
526  & THEN
527  IF (abs(ygl).GT.abs(ymi(imt))) ygl=ymi(imt)
528  ENDIF
529  480 CONTINUE
530  ENDIF
531 C...YP1 = Rapidity IP1
532  yp1=ymi(imp1)
533 C...If fictitious gluon
534  IF (yp1.EQ.100d0) THEN
535  yp1=(3-2*js)*yp1
536  ida1=mod(k(ip1,4),mstu(5))
537  ida2=mod(k(ip1,5),mstu(5))
538  DO 490 imt=1,nmi(js)
539 C...Select (arbitrarily) the most central daughter.
540  IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2)
541  & THEN
542  IF (abs(yp1).GT.abs(ymi(imt))) yp1=ymi(imt)
543  ENDIF
544  490 CONTINUE
545  ENDIF
546 C...YP2 = Rapidity of mother system
547  IF (k(ip2,2).NE.88) THEN
548  DO 500 imt=1,nmi(js)
549  IF (imi(js,imt,1).EQ.ip2) yp2=ymi(imt)
550  500 CONTINUE
551 C...If fictitious gluon
552  IF (yp2.EQ.100d0) THEN
553  yp2=(3-2*js)*yp2
554  ida1=mod(k(ip2,4),mstu(5))
555  ida2=mod(k(ip2,5),mstu(5))
556  DO 510 imt=1,nmi(js)
557 C...Select (arbitrarily) the most central daughter.
558  IF (imi(js,imt,1).EQ.ida1.OR.imi(js,imt,1).EQ.ida2
559  & ) THEN
560  IF (abs(yp2).GT.abs(ymi(imt))) yp2=ymi(imt)
561  ENDIF
562  510 CONTINUE
563  ENDIF
564 C...Assign (arbitrarily) 100D0 to junction also
565  ELSE
566  yp2=(3-2*js)*100d0
567  ENDIF
568  rl=abs(ygl-yp1)+abs(ygl-yp2)
569  ELSEIF (mstp(89).EQ.2) THEN
570 C...Lambda ordering:
571 C...Compute lambda measure for this insertion.
572  rl=1d0
573  DO 520 ist=1,6
574  istr(ist)=0
575  520 CONTINUE
576 C...If IP2 is junction, not caught below.
577  IF (jcp2.EQ.0) THEN
578  itju=mod(k(ip2,4)/mstu(5),mstu(5))
579 C...Anti-junction is colour endpoint et vv., always on JCG2.
580  istr(5-itju)=ip2
581  ENDIF
582  DO 530 i=mint(84)+1,n
583  IF (k(i,1).LT.10) THEN
584 C...The new string pieces
585  IF (mct(i,1).EQ.jcg1) istr(1)=i
586  IF (mct(i,2).EQ.jcg1) istr(2)=i
587  IF (mct(i,1).EQ.jcg2) istr(3)=i
588  IF (mct(i,2).EQ.jcg2) istr(4)=i
589  ENDIF
590  530 CONTINUE
591 C...Also identify junctions as string endpoints.
592  DO 540 i=mint(84)+1,n
593  icmo=mod(k(i,4)/mstu(5),mstu(5))
594  iamo=mod(k(i,5)/mstu(5),mstu(5))
595 C...Find partons adjacent to junctions.
596  IF (k(icmo,1).EQ.42.AND.mct(i,1).EQ.jcg1.AND.istr(2)
597  & .EQ.0) istr(2) = icmo
598  IF (k(iamo,1).EQ.42.AND.mct(i,2).EQ.jcg1.AND.istr(1)
599  & .EQ.0) istr(1) = iamo
600  IF (k(icmo,1).EQ.42.AND.mct(i,1).EQ.jcg2.AND.istr(4)
601  & .EQ.0) istr(4) = icmo
602  IF (k(iamo,1).EQ.42.AND.mct(i,2).EQ.jcg2.AND.istr(3)
603  & .EQ.0) istr(3) = iamo
604  540 CONTINUE
605 C...The old string piece
606  istr(5)=istr(1+2*manti)
607  istr(6)=istr(4-2*manti)
608  rl=max(1d0,four(istr(1),istr(2)))*max(1d0,four(istr(3)
609  & ,istr(4)))/max(1d0,four(istr(5),istr(6)))
610  rl=log(rl)
611  ENDIF
612 C...Allow some breadth to speed things up.
613  IF (abs(1d0-rl/rlopt).LT.0.05d0) THEN
614  nopt=nopt+1
615  ELSEIF (rl.GT.rlopt) THEN
616  goto 470
617  ELSE
618  nopt=1
619  rlopt=rl
620  ENDIF
621 C...INSR(NOPT,1)=Gluon colour mother
622 C...INSR(NOPT,2)=Gluon
623 C...INSR(NOPT,3)=Gluon anticolour mother
624  IF (nopt.GT.1000) goto 470
625  insr(nopt,1+2*manti)=ip2
626  insr(nopt,2)=igl
627  insr(nopt,3-2*manti)=ip1
628  IF (mstp(89).GT.0.OR.nopt.EQ.0) goto 470
629  ENDIF
630  IF (mstp(89).GT.0.OR.nopt.EQ.0) goto 460
631  ENDIF
632 C...Reset link test information.
633  DO 550 i=mint(84)+1,n
634  k(i,4)=mod(k(i,4),mstu(5)**2)
635  k(i,5)=mod(k(i,5),mstu(5)**2)
636  550 CONTINUE
637  IF (mstp(89).GT.0.OR.nopt.EQ.0) goto 450
638  ENDIF
639 C...Now we have a list of best gluon insertions, none of which cause
640 C...singlets to arise. If list is empty, try again a few times. Note:
641 C...this should never happen if we have a meson with a gluon inserted
642 C...in the beam remnant, since that breaks up the colour line.
643  IF (nopt.EQ.0) THEN
644 C...Abandon BR-g-BR suppression for retries. This is not serious, it
645 C...just means we happened to start with trying a bad sequence.
646  parp80=1d0
647  IF (mretry.LE.10.AND.(itjunc(1).NE.0.OR.jst(1,3).EQ.0).and
648  & .(itjunc(2).NE.0.OR.jst(2,3).EQ.0)) THEN
649  mretry=mretry+1
650  DO 590 js=1,2
651  IF (itjunc(js).NE.0) THEN
652  jst(js,1)=iv(js,1)
653  jst(js,2)=iv(js,2)
654  jst(js,3)=iv(js,3)
655 C...Reset valence quark parent pointers
656  DO 560 i=mint(53)+1,n
657  IF (k(i,2).EQ.88.AND.k(i,3).EQ.js) iju=i
658  560 CONTINUE
659  manti=itjunc(js)-1
660 C...Set (anti)colour mother = junction.
661  DO 570 jv=1,3
662  k(iv(js,jv),4+manti)=mod(k(iv(js,jv),4+manti),mstu(5))
663  & +mstu(5)*iju
664  570 CONTINUE
665  ELSE
666 C...Same for mesons. JST unchanged, so needn't be restored.
667  iq=jst(js,1)
668  iqbar=jst(js,2)
669  k(iq,4)=mod(k(iq,4),mstu(5))+mstu(5)*iqbar
670  k(iqbar,5)=mod(k(iqbar,5),mstu(5))+mstu(5)*iq
671  ENDIF
672 C...Also reset gluon parent pointers.
673  ng(js)=0
674  DO 580 im=1,nmi(js)
675  i=imi(js,im,1)
676  IF (k(i,2).EQ.21) THEN
677  k(i,4)=mod(k(i,4),mstu(5))
678  k(i,5)=mod(k(i,5),mstu(5))
679  ng(js)=ng(js)+1
680  ENDIF
681  580 CONTINUE
682  590 CONTINUE
683 C...Reset colour tags
684  DO 600 i=mint(84)+1,n
685  mct(i,1)=mco(i,1)
686  mct(i,2)=mco(i,2)
687  600 CONTINUE
688  goto 400
689  ELSE
690  IF(nerrpr.LT.5) THEN
691  nerrpr=nerrpr+1
692  CALL pylist(4)
693  CALL pyerrm(19,'(PYMIHK:) No physical colour flow found!')
694  WRITE(mstu(11),*) 'NG:', ng,' MOUT:', mout(js)
695  ENDIF
696 C...Kill event and start another.
697  mint(51)=1
698  RETURN
699  ENDIF
700  ELSE
701 C...Select between insertions, suppressing insertions wholly in the BR.
702  iin=pyr(0)*nopt+1
703  610 iin=mod(iin,nopt)+1
704  IF (insr(iin,1).GT.mint(53).AND.insr(iin,3).GT.mint(53)
705  & .AND.mout(js).NE.0.AND.pyr(0).GT.parp80) goto 610
706  ENDIF
707 
708 C...Now we know which gluon to insert where. Colour tags in JCCO and
709 C...colour connection information should be updated, NG(JS) should be
710 C...counted down, and a new loop performed if there are still gluons
711 C...left on any side.
712  icm=insr(iin,1)
713  iacm=insr(iin,3)
714  igl=insr(iin,2)
715 C...JCG : Original gluon colour tag
716 C...JCAG: Original gluon anticolour tag.
717 C...JCM : Original anticolour tag of gluon colour mother
718 C...JACM: Original colour tag of gluon anticolour mother
719  jcg=mco(igl,1)
720  jcm=mco(icm,2)
721  jacg=mco(igl,2)
722  jacm=mco(iacm,1)
723 
724  CALL pymihg(jacm,jacg,jcm,jcg)
725  IF (maccpt.EQ.0) THEN
726  IF(nerrpr.LT.5) THEN
727  nerrpr=nerrpr+1
728  CALL pylist(4)
729  CALL pyerrm(11,'(PYMIHK:) Unphysical colour flow!')
730  WRITE(mstu(11),*) 'attaching', igl,' between', icm, iacm
731  ENDIF
732 C...Kill event and start another.
733  mint(51)=1
734  RETURN
735  ELSE
736 C...If everything went fine, store new JCCN in JCCO.
737  ncc=ncc+1
738  DO 620 icc=1,ncc
739  jcco(icc,1)=jccn(icc,1)
740  jcco(icc,2)=jccn(icc,2)
741  620 CONTINUE
742  ENDIF
743 
744 C...One gluon attached is counted as equivalent to one end outside.
745  mout(js)=1
746 C...Set IGL colour mother = ICM.
747  k(igl,4)=mod(k(igl,4),mstu(5))+mstu(5)*icm
748 C...Set ICM anticolour mother = IGL colour.
749  IF (k(icm,2).NE.88) THEN
750  k(icm,5)=mod(k(icm,5),mstu(5))+mstu(5)*igl
751  ELSE
752 C...If ICM is junction, just update JST array for now.
753  DO 630 msj=1,3
754  IF (jst(js,msj).EQ.iacm) jst(js,msj)=igl
755  630 CONTINUE
756  ENDIF
757 C...Set IGL anticolour mother = IACM.
758  k(igl,5)=mod(k(igl,5),mstu(5))+mstu(5)*iacm
759 C...Set IACM anticolour mother = IGL anticolour.
760  IF (k(iacm,2).NE.88) THEN
761  k(iacm,4)=mod(k(iacm,4),mstu(5))+mstu(5)*igl
762  ELSE
763 C...If IACM is junction, just update JST array for now.
764  DO 640 msj=1,3
765  IF (jst(js,msj).EQ.icm) jst(js,msj)=igl
766  640 CONTINUE
767  ENDIF
768 C...Count down # unconnected gluons.
769  ng(js)=ng(js)-1
770  ENDIF
771  IF (ng(1).GT.0.OR.ng(2).GT.0) goto 440
772 
773  DO 840 js=1,2
774 C...Collapse fictitious gluons.
775  DO 670 igl=mint(53)+1,n
776  IF (k(igl,2).EQ.21.AND.k(igl,3).EQ.mint(83)+js.AND.
777  & k(igl,1).EQ.14) THEN
778  icm=k(igl,4)/mstu(5)
779  iam=k(igl,5)/mstu(5)
780  icd=mod(k(igl,4),mstu(5))
781  iad=mod(k(igl,5),mstu(5))
782 C...Set gluon daughters pointing to gluon mothers
783  k(iad,5)=mod(k(iad,5),mstu(5))+mstu(5)*iam
784  k(icd,4)=mod(k(icd,4),mstu(5))+mstu(5)*icm
785 C...Set gluon mothers pointing to gluon daughters.
786  IF (k(icm,2).NE.88) THEN
787  k(icm,5)=mod(k(icm,5),mstu(5))+mstu(5)*icd
788  ELSE
789 C...Special case: mother=junction. Just update JST array for now.
790  DO 650 msj=1,3
791  IF (jst(js,msj).EQ.igl) jst(js,msj)=icd
792  650 CONTINUE
793  ENDIF
794  IF (k(iam,2).NE.88) THEN
795  k(iam,4)=mod(k(iam,4),mstu(5))+mstu(5)*iad
796  ELSE
797  DO 660 msj=1,3
798  IF (jst(js,msj).EQ.igl) jst(js,msj)=iad
799  660 CONTINUE
800  ENDIF
801  ENDIF
802  670 CONTINUE
803 
804 C...Erase collapsed gluons from NMI and IMI (but keep them in ER)
805  im=nmi(js)+1
806  680 im=im-1
807  IF (im.GT.mint(31).AND.k(imi(js,im,1),2).NE.21) goto 680
808  IF (im.GT.mint(31)) THEN
809  nmi(js)=nmi(js)-1
810  DO 690 imr=im,nmi(js)
811  imi(js,imr,1)=imi(js,imr+1,1)
812  imi(js,imr,2)=imi(js,imr+1,2)
813  690 CONTINUE
814  goto 680
815  ENDIF
816 
817 C...Finally, connect junction.
818  IF (itjunc(js).NE.0) THEN
819  DO 700 i=mint(53)+1,n
820  IF (k(i,2).EQ.88.AND.k(i,3).EQ.mint(83)+js) iju=i
821  700 CONTINUE
822 C...NBRJQ counts # of jq, NBRVQ # of jv, inside BR.
823  nbrjq =0
824  nbrvq =0
825  DO 720 msj=1,3
826  idq(msj)=0
827 C...Find jq with no glue inbetween inside beam remnant.
828  IF (jst(js,msj).GT.mint(53).AND.iabs(k(jst(js,msj),2)).LE.5)
829  & THEN
830  nbrjq=nbrjq+1
831 C...Set IDQ = -I if q non-valence and = +I if q valence.
832  idq(nbrjq)=-jst(js,msj)
833  DO 710 jv=1,3
834  IF (iv(js,jv).EQ.jst(js,msj)) THEN
835  idq(nbrjq)=jst(js,msj)
836  nbrvq=nbrvq+1
837  ENDIF
838  710 CONTINUE
839  ENDIF
840  i12=mod(msj+1,2)
841  i45=5
842  IF (msj.EQ.3) i45=4
843  k(iju,i45)=k(iju,i45)+(mstu(5)**i12)*jst(js,msj)
844  720 CONTINUE
845 
846 C...Check if diquark can be formed.
847  IF ((mstp(88).GE.0.AND.nbrvq.GE.2).OR.(nbrjq.GE.2.AND.mstp(88)
848  & .GE.1)) THEN
849 C...If there is less than 2 valence quarks connected to junction
850 C...and MSTP(88)>1, use random non-valence quarks to fill up.
851  IF (nbrvq.LE.1) THEN
852  ndiq=nbrvq
853  730 jflip=nbrjq*pyr(0)+1
854  IF (idq(jflip).LT.0) THEN
855  idq(jflip)=-idq(jflip)
856  ndiq=ndiq+1
857  ENDIF
858  IF (ndiq.LE.1) goto 730
859  ENDIF
860 C...Place selected quarks first in IDQ, ordered in flavour.
861  DO 740 jdq=1,3
862  IF (idq(jdq).LE.0) THEN
863  itemp1 = idq(jdq)
864  idq(jdq)= idq(3)
865  idq(3) = -itemp1
866  IF (iabs(k(idq(1),2)).LT.iabs(k(idq(2),2))) THEN
867  itemp1 = idq(1)
868  idq(1) = idq(2)
869  idq(2) = itemp1
870  ENDIF
871  ENDIF
872  740 CONTINUE
873 C...Choose diquark spin.
874  IF (nbrvq.EQ.2) THEN
875 C...If the selected quarks are both valence, we may use SU(6) rules
876 C...to figure out which spin the diquark has, by a subdivision of the
877 C...original beam hadron into the selected diquark system plus a kicked
878 C...out quark, IKO.
879  jko=6
880  DO 760 jdq=1,2
881  DO 750 jv=1,3
882  IF (idq(jdq).EQ.iv(js,jv)) jko=jko-jv
883  750 CONTINUE
884  760 CONTINUE
885  iko=iv(js,jko)
886  CALL pyspli(mint(10+js),k(iko,2),kfdum,kfdq)
887  ELSE
888 C...If one or more of the selected quarks are not valence, we cannot use
889 C...SU(6) subdivisions of the original beam hadron. Instead, with the
890 C...flavours of the diquark already selected, we assume for now
891 C...50:50 spin-1:spin-0 (where spin-0 possible).
892  kfdq=1000*k(idq(1),2)+100*k(idq(2),2)
893  is=3
894  IF (k(idq(1),2).NE.k(idq(2),2).AND.
895  & (1d0+3d0*parj(4))*pyr(0).LT.1d0) is=1
896  kfdq=kfdq+isign(is,kfdq)
897  ENDIF
898 
899 C...Collapse diquark-j-quark system to baryon, if allowed and possible.
900 C...Note: third quark can per definition not also be valence,
901 C...therefore we can only do this if we are allowed to use sea quarks.
902  770 IF (idq(3).NE.0.AND.mstp(88).GE.2) THEN
903  ntry=0
904  780 ntry=ntry+1
905  CALL pykfdi(kfdq,k(iabs(idq(3)),2),kfdum,kfbar)
906  IF (kfbar.EQ.0.AND.ntry.LE.100) THEN
907  goto 780
908  ELSEIF(ntry.GT.100) THEN
909 C...If no baryon can be found, give up and form diquark.
910  idq(3)=0
911  goto 770
912  ELSE
913 C...Replace junction by baryon.
914  k(iju,1)=1
915  k(iju,2)=kfbar
916  k(iju,3)=mint(83)+js
917  k(iju,4)=0
918  k(iju,5)=0
919  p(iju,5)=pymass(kfbar)
920  DO 790 msj=1,3
921 C...Prepare removal of participating quarks from ER.
922  k(jst(js,msj),1)=-1
923  790 CONTINUE
924  ENDIF
925  ELSE
926 C...If collapse to baryon not possible or not allowed, replace junction
927 C...by diquark. This way, collapsed gluons that were pointing at the
928 C...junction will now point (correctly) at diquark.
929  manti=itjunc(js)-1
930  k(iju,1)=3
931  k(iju,2)=kfdq
932  k(iju,3)=mint(83)+js
933  k(iju,4)=0
934  k(iju,5)=0
935  DO 800 msj=1,3
936  ip=jst(js,msj)
937  IF (ip.NE.idq(1).AND.ip.NE.idq(2)) THEN
938  k(iju,4+manti)=0
939  k(iju,5-manti)=ip*mstu(5)
940  k(ip,4+manti)=mod(k(ip,4+manti),mstu(5))+
941  & mstu(5)*iju
942  mct(iju,2-manti)=mct(ip,1+manti)
943  ELSE
944 C...Prepare removal of participating quarks from ER.
945  k(ip,1)=-1
946  ENDIF
947  800 CONTINUE
948  ENDIF
949 
950 C...Update so ER pointers to collapsed quarks
951 C...now go to collapsed object.
952  DO 820 i=mint(84)+1,n
953  IF ((k(i,3).EQ.mint(83)+js.OR.k(i,3).EQ.mint(83)+2+js).and
954  & .k(i,1).GT.0) THEN
955  DO 810 isid=4,5
956  imo=k(i,isid)/mstu(5)
957  ida=mod(k(i,isid),mstu(5))
958  IF (imo.GT.0) THEN
959  IF (k(imo,1).EQ.-1) imo=iju
960  ENDIF
961  IF (ida.GT.0) THEN
962  IF (k(ida,1).EQ.-1) ida=iju
963  ENDIF
964  k(i,isid)=ida+mstu(5)*imo
965  810 CONTINUE
966  ENDIF
967  820 CONTINUE
968  ENDIF
969  ENDIF
970 
971 C...Finally, if beam remnant is empty, insert a gluon in beam remnant.
972 C...(this only happens for baryons, where we want to force the gluon
973 C...to sit next to the junction. Mesons handled above.)
974  IF (nbrtot(js).EQ.0) THEN
975  n=n+1
976  DO 830 ix=1,5
977  k(n,ix)=0
978  p(n,ix)=0d0
979  v(n,ix)=0d0
980  830 CONTINUE
981  igl=n
982  k(igl,1)=3
983  k(igl,2)=21
984  k(igl,3)=mint(83)+js
985  IF (itjunc(js).NE.0) THEN
986 C...Incoming baryons. Pick random leg in JST (NVSUM = 3 for baryons)
987  jleg=pyr(0)*nvsum(js)+1
988  i1=jst(js,jleg)
989  jst(js,jleg)=igl
990  jct=mct(i1,itjunc(js))
991  mct(igl,3-itjunc(js))=jct
992  nct=nct+1
993  mct(igl,itjunc(js))=nct
994  manti=itjunc(js)-1
995  ELSE
996 C...Meson. Should not happen.
997  CALL pyerrm(19,'(PYMIHK:) Empty meson beam remnant')
998  IF(nerrpr.LT.5) THEN
999  WRITE(mstu(11),*) 'This should not have been possible!'
1000  CALL pylist(4)
1001  nerrpr=nerrpr+1
1002  ENDIF
1003  mint(51)=1
1004  RETURN
1005  ENDIF
1006  i2=mod(k(i1,4+manti)/mstu(5),mstu(5))
1007  k(i1,4+manti)=mod(k(i1,4+manti),mstu(5))+mstu(5)*igl
1008  k(igl,5-manti)=mod(k(igl,5-manti),mstu(5))+mstu(5)*i1
1009  k(igl,4+manti)=mod(k(igl,4+manti),mstu(5))+mstu(5)*i2
1010  IF (k(i2,2).NE.88) THEN
1011  k(i2,5-manti)=mod(k(i2,5-manti),mstu(5))+mstu(5)*igl
1012  ELSE
1013  IF (mod(k(i2,4),mstu(5)).EQ.i1) THEN
1014  k(i2,4)=(k(i2,4)/mstu(5))*mstu(5)+igl
1015  ELSEIF(mod(k(i2,5)/mstu(5),mstu(5)).EQ.i1) THEN
1016  k(i2,5)=mod(k(i2,5),mstu(5))+mstu(5)*igl
1017  ELSE
1018  k(i2,5)=(k(i2,5)/mstu(5))*mstu(5)+igl
1019  ENDIF
1020  ENDIF
1021  ENDIF
1022  840 CONTINUE
1023 
1024 C...Remove collapsed quarks and junctions from ER and update IMI.
1025  CALL pyedit(11)
1026 
1027 C...Also update beam remnant part of IMI.
1028  nmi(1)=mint(31)
1029  nmi(2)=mint(31)
1030  DO 850 i=mint(53)+1,n
1031  IF (k(i,1).LE.0) goto 850
1032 C...Restore BR quark/diquark/baryon pointers in IMI.
1033  IF ((k(i,2).NE.21.OR.k(i,1).NE.14).AND.k(i,2).NE.88) THEN
1034  js=k(i,3)-mint(83)
1035  nmi(js)=nmi(js)+1
1036  imi(js,nmi(js),1)=i
1037  imi(js,nmi(js),2)=0
1038  ENDIF
1039  850 CONTINUE
1040 
1041 C...Restore companion information from collapsed gluons.
1042  DO 870 i=mint(53)+1,n
1043  IF (k(i,2).EQ.21.AND.k(i,1).EQ.14) THEN
1044  js=k(i,3)-mint(83)
1045  jcd=mod(k(i,4),mstu(5))
1046  jad=mod(k(i,5),mstu(5))
1047  DO 860 im=1,nmi(js)
1048  IF (imi(js,im,1).EQ.jcd) imc=im
1049  IF (imi(js,im,1).EQ.jad) ima=im
1050  860 CONTINUE
1051  imi(js,imc,2)=imi(js,ima,1)
1052  imi(js,ima,2)=imi(js,imc,1)
1053  ENDIF
1054  870 CONTINUE
1055 
1056 C...Renumber colour lines (since some have disappeared)
1057  jct=0
1058  jcd=0
1059  880 jct=jct+1
1060  mfound=0
1061  i=mint(84)
1062  890 i=i+1
1063  IF (i.EQ.n+1) THEN
1064  IF (mfound.EQ.0) jcd=jcd+1
1065  ELSEIF (mct(i,1).EQ.jct.AND.k(i,1).GE.1) THEN
1066  mct(i,1)=jct-jcd
1067  mfound=1
1068  ELSEIF (mct(i,2).EQ.jct.AND.k(i,1).GE.1) THEN
1069  mct(i,2)=jct-jcd
1070  mfound=1
1071  ENDIF
1072  IF (i.LE.n) goto 890
1073  IF (jct.LT.nct) goto 880
1074  nct=jct-jcd
1075 
1076 C...Reset hard interaction subsystems to their CM frames.
1077  IF (iboost.EQ.1) THEN
1078  DO 900 im=1,mint(31)
1079  beta=-(xmi(1,im)-xmi(2,im))/(xmi(1,im)+xmi(2,im))
1080  CALL pyrobo(imisep(im-1)+1,imisep(im),0d0,0d0,0d0,0d0,beta)
1081  900 CONTINUE
1082 C...Zero beam remnant longitudinal momenta and energies
1083  DO 910 i=mint(53)+1,n
1084  p(i,3)=0d0
1085  p(i,4)=0d0
1086  910 CONTINUE
1087  ELSE
1088  CALL pyerrm(9
1089  & ,'(PYMIHK:) Inconsistent kinematics. Too many boosts.')
1090 C...Kill event and start another.
1091  mint(51)=1
1092  RETURN
1093  ENDIF
1094 
1095  9999 RETURN
1096  END