Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lustrf.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lustrf.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lustrf(IP)
5 C...Purpose: to handle the fragmentation of an arbitrary colour singlet
6 C...jet system according to the Lund string fragmentation model.
7  IMPLICIT DOUBLE PRECISION(d)
8  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  dimension dps(5),kfl(3),pmq(3),px(3),py(3),gam(3),ie(2),pr(2),
15  &in(9),dhm(4),dhg(4),dp(5,5),irank(2),mju(4),iju(3),pju(5,5),
16  &tju(5),kfjh(2),njs(2),kfjs(2),pjs(4,5)
17 
18 C...Function: four-product of two vectors.
19  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)
20  dfour(i,j)=dp(i,4)*dp(j,4)-dp(i,1)*dp(j,1)-dp(i,2)*dp(j,2)-
21  &dp(i,3)*dp(j,3)
22 
23 C...Reset counters. Identify parton system.
24  mstj(91)=0
25  nsav=n
26  np=0
27  kqsum=0
28  DO 100 j=1,5
29  100 dps(j)=0.
30  mju(1)=0
31  mju(2)=0
32  i=ip-1
33  110 i=i+1
34  IF(i.GT.min(n,mstu(4)-mstu(32))) THEN
35  CALL luerrm(12,'(LUSTRF:) failed to reconstruct jet system')
36  IF(mstu(21).GE.1) RETURN
37  ENDIF
38  IF(k(i,1).NE.1.AND.k(i,1).NE.2.AND.k(i,1).NE.41) goto 110
39  kc=lucomp(k(i,2))
40  IF(kc.EQ.0) goto 110
41  kq=kchg(kc,2)*isign(1,k(i,2))
42  IF(kq.EQ.0) goto 110
43  IF(n+5*np+11.GT.mstu(4)-mstu(32)-5) THEN
44  CALL luerrm(11,'(LUSTRF:) no more memory left in LUJETS')
45  IF(mstu(21).GE.1) RETURN
46  ENDIF
47 
48 C...Take copy of partons to be considered. Check flavour sum.
49  np=np+1
50  DO 120 j=1,5
51  k(n+np,j)=k(i,j)
52  p(n+np,j)=p(i,j)
53  120 dps(j)=dps(j)+p(i,j)
54  k(n+np,3)=i
55  IF(p(n+np,4)**2.LT.p(n+np,1)**2+p(n+np,2)**2+p(n+np,3)**2) THEN
56  p(n+np,4)=sqrt(p(n+np,1)**2+p(n+np,2)**2+p(n+np,3)**2+
57  & p(n+np,5)**2)
58  dps(4)=dps(4)+max(0.,p(n+np,4)-p(i,4))
59  ENDIF
60  IF(kq.NE.2) kqsum=kqsum+kq
61  IF(k(i,1).EQ.41) THEN
62  kqsum=kqsum+2*kq
63  IF(kqsum.EQ.kq) mju(1)=n+np
64  IF(kqsum.NE.kq) mju(2)=n+np
65  ENDIF
66  IF(k(i,1).EQ.2.OR.k(i,1).EQ.41) goto 110
67  IF(kqsum.NE.0) THEN
68  CALL luerrm(12,'(LUSTRF:) unphysical flavour combination')
69  IF(mstu(21).GE.1) RETURN
70  ENDIF
71 
72 C...Boost copied system to CM frame (for better numerical precision).
73  CALL ludbrb(n+1,n+np,0.,0.,-dps(1)/dps(4),-dps(2)/dps(4),
74  &-dps(3)/dps(4))
75 
76 C...Search for very nearby partons that may be recombined.
77  ntryr=0
78  paru12=paru(12)
79  paru13=paru(13)
80  mju(3)=mju(1)
81  mju(4)=mju(2)
82  nr=np
83  130 IF(nr.GE.3) THEN
84  pdrmin=2.*paru12
85  DO 140 i=n+1,n+nr
86  IF(i.EQ.n+nr.AND.iabs(k(n+1,2)).NE.21) goto 140
87  i1=i+1
88  IF(i.EQ.n+nr) i1=n+1
89  IF(k(i,1).EQ.41.OR.k(i1,1).EQ.41) goto 140
90  IF(mju(1).NE.0.AND.i1.LT.mju(1).AND.iabs(k(i1,2)).NE.21)
91  & goto 140
92  IF(mju(2).NE.0.AND.i.GT.mju(2).AND.iabs(k(i,2)).NE.21) goto 140
93  pap=sqrt((p(i,1)**2+p(i,2)**2+p(i,3)**2)*(p(i1,1)**2+
94  & p(i1,2)**2+p(i1,3)**2))
95  pvp=p(i,1)*p(i1,1)+p(i,2)*p(i1,2)+p(i,3)*p(i1,3)
96  pdr=4.*(pap-pvp)**2/(paru13**2*pap+2.*(pap-pvp))
97  IF(pdr.LT.pdrmin) THEN
98  ir=i
99  pdrmin=pdr
100  ENDIF
101  140 CONTINUE
102 
103 C...Recombine very nearby partons to avoid machine precision problems.
104  IF(pdrmin.LT.paru12.AND.ir.EQ.n+nr) THEN
105  DO 150 j=1,4
106  150 p(n+1,j)=p(n+1,j)+p(n+nr,j)
107  p(n+1,5)=sqrt(max(0.,p(n+1,4)**2-p(n+1,1)**2-p(n+1,2)**2-
108  & p(n+1,3)**2))
109  nr=nr-1
110  goto 130
111  ELSEIF(pdrmin.LT.paru12) THEN
112  DO 160 j=1,4
113  160 p(ir,j)=p(ir,j)+p(ir+1,j)
114  p(ir,5)=sqrt(max(0.,p(ir,4)**2-p(ir,1)**2-p(ir,2)**2-
115  & p(ir,3)**2))
116  DO 170 i=ir+1,n+nr-1
117  k(i,2)=k(i+1,2)
118  DO 170 j=1,5
119  170 p(i,j)=p(i+1,j)
120  IF(ir.EQ.n+nr-1) k(ir,2)=k(n+nr,2)
121  nr=nr-1
122  IF(mju(1).GT.ir) mju(1)=mju(1)-1
123  IF(mju(2).GT.ir) mju(2)=mju(2)-1
124  goto 130
125  ENDIF
126  ENDIF
127  ntryr=ntryr+1
128 
129 C...Reset particle counter. Skip ahead if no junctions are present;
130 C...this is usually the case!
131  nrs=max(5*nr+11,np)
132  ntry=0
133  180 ntry=ntry+1
134  IF(ntry.GT.100.AND.ntryr.LE.4) THEN
135  paru12=4.*paru12
136  paru13=2.*paru13
137  goto 130
138  ELSEIF(ntry.GT.100) THEN
139  CALL luerrm(14,'(LUSTRF:) caught in infinite loop')
140  IF(mstu(21).GE.1) RETURN
141  ENDIF
142  i=n+nrs
143  IF(mju(1).EQ.0.AND.mju(2).EQ.0) goto 500
144  DO 490 jt=1,2
145  njs(jt)=0
146  IF(mju(jt).EQ.0) goto 490
147  js=3-2*jt
148 
149 C...Find and sum up momentum on three sides of junction. Check flavours.
150  DO 190 iu=1,3
151  iju(iu)=0
152  DO 190 j=1,5
153  190 pju(iu,j)=0.
154  iu=0
155  DO 200 i1=n+1+(jt-1)*(nr-1),n+nr+(jt-1)*(1-nr),js
156  IF(k(i1,2).NE.21.AND.iu.LE.2) THEN
157  iu=iu+1
158  iju(iu)=i1
159  ENDIF
160  DO 200 j=1,4
161  200 pju(iu,j)=pju(iu,j)+p(i1,j)
162  DO 210 iu=1,3
163  210 pju(iu,5)=sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
164  IF(k(iju(3),2)/100.NE.10*k(iju(1),2)+k(iju(2),2).AND.
165  &k(iju(3),2)/100.NE.10*k(iju(2),2)+k(iju(1),2)) THEN
166  CALL luerrm(12,'(LUSTRF:) unphysical flavour combination')
167  IF(mstu(21).GE.1) RETURN
168  ENDIF
169 
170 C...Calculate (approximate) boost to rest frame of junction.
171  t12=(pju(1,1)*pju(2,1)+pju(1,2)*pju(2,2)+pju(1,3)*pju(2,3))/
172  &(pju(1,5)*pju(2,5))
173  t13=(pju(1,1)*pju(3,1)+pju(1,2)*pju(3,2)+pju(1,3)*pju(3,3))/
174  &(pju(1,5)*pju(3,5))
175  t23=(pju(2,1)*pju(3,1)+pju(2,2)*pju(3,2)+pju(2,3)*pju(3,3))/
176  &(pju(2,5)*pju(3,5))
177  t11=sqrt((2./3.)*(1.-t12)*(1.-t13)/(1.-t23))
178  t22=sqrt((2./3.)*(1.-t12)*(1.-t23)/(1.-t13))
179  tsq=sqrt((2.*t11*t22+t12-1.)*(1.+t12))
180  t1f=(tsq-t22*(1.+t12))/(1.-t12**2)
181  t2f=(tsq-t11*(1.+t12))/(1.-t12**2)
182  DO 220 j=1,3
183  220 tju(j)=-(t1f*pju(1,j)/pju(1,5)+t2f*pju(2,j)/pju(2,5))
184  tju(4)=sqrt(1.+tju(1)**2+tju(2)**2+tju(3)**2)
185  DO 230 iu=1,3
186  230 pju(iu,5)=tju(4)*pju(iu,4)-tju(1)*pju(iu,1)-tju(2)*pju(iu,2)-
187  &tju(3)*pju(iu,3)
188 
189 C...Put junction at rest if motion could give inconsistencies.
190  IF(pju(1,5)+pju(2,5).GT.pju(1,4)+pju(2,4)) THEN
191  DO 240 j=1,3
192  240 tju(j)=0.
193  tju(4)=1.
194  pju(1,5)=pju(1,4)
195  pju(2,5)=pju(2,4)
196  pju(3,5)=pju(3,4)
197  ENDIF
198 
199 C...Start preparing for fragmentation of two strings from junction.
200  ista=i
201  DO 470 iu=1,2
202  ns=iju(iu+1)-iju(iu)
203 
204 C...Junction strings: find longitudinal string directions.
205  DO 260 is=1,ns
206  is1=iju(iu)+is-1
207  is2=iju(iu)+is
208  DO 250 j=1,5
209  dp(1,j)=0.5*p(is1,j)
210  IF(is.EQ.1) dp(1,j)=p(is1,j)
211  dp(2,j)=0.5*p(is2,j)
212  250 IF(is.EQ.ns) dp(2,j)=-pju(iu,j)
213  IF(is.EQ.ns) dp(2,4)=sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
214  IF(is.EQ.ns) dp(2,5)=0.
215  dp(3,5)=dfour(1,1)
216  dp(4,5)=dfour(2,2)
217  dhkc=dfour(1,2)
218  IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.) THEN
219  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
220  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
221  dp(3,5)=0d0
222  dp(4,5)=0d0
223  dhkc=dfour(1,2)
224  ENDIF
225  dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
226  dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
227  dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
228  in1=n+nr+4*is-3
229  p(in1,5)=sqrt(dp(3,5)+2.*dhkc+dp(4,5))
230  DO 260 j=1,4
231  p(in1,j)=(1.+dhk1)*dp(1,j)-dhk2*dp(2,j)
232  260 p(in1+1,j)=(1.+dhk2)*dp(2,j)-dhk1*dp(1,j)
233 
234 C...Junction strings: initialize flavour, momentum and starting pos.
235  isav=i
236  270 ntry=ntry+1
237  IF(ntry.GT.100.AND.ntryr.LE.4) THEN
238  paru12=4.*paru12
239  paru13=2.*paru13
240  goto 130
241  ELSEIF(ntry.GT.100) THEN
242  CALL luerrm(14,'(LUSTRF:) caught in infinite loop')
243  IF(mstu(21).GE.1) RETURN
244  ENDIF
245  i=isav
246  irankj=0
247  ie(1)=k(n+1+(jt/2)*(np-1),3)
248  in(4)=n+nr+1
249  in(5)=in(4)+1
250  in(6)=n+nr+4*ns+1
251  DO 280 jq=1,2
252  DO 280 in1=n+nr+2+jq,n+nr+4*ns-2+jq,4
253  p(in1,1)=2-jq
254  p(in1,2)=jq-1
255  280 p(in1,3)=1.
256  kfl(1)=k(iju(iu),2)
257  px(1)=0.
258  py(1)=0.
259  gam(1)=0.
260  DO 290 j=1,5
261  290 pju(iu+3,j)=0.
262 
263 C...Junction strings: find initial transverse directions.
264  DO 300 j=1,4
265  dp(1,j)=p(in(4),j)
266  dp(2,j)=p(in(4)+1,j)
267  dp(3,j)=0.
268  300 dp(4,j)=0.
269  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
270  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
271  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
272  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
273  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
274  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
275  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
276  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
277  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
278  dhc12=dfour(1,2)
279  dhcx1=dfour(3,1)/dhc12
280  dhcx2=dfour(3,2)/dhc12
281  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
282  dhcy1=dfour(4,1)/dhc12
283  dhcy2=dfour(4,2)/dhc12
284  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
285  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
286  DO 310 j=1,4
287  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
288  p(in(6),j)=dp(3,j)
289  310 p(in(6)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
290  &dhcyx*dp(3,j))
291 
292 C...Junction strings: produce new particle, origin.
293  320 i=i+1
294  IF(2*i-nsav.GE.mstu(4)-mstu(32)-5) THEN
295  CALL luerrm(11,'(LUSTRF:) no more memory left in LUJETS')
296  IF(mstu(21).GE.1) RETURN
297  ENDIF
298  irankj=irankj+1
299  k(i,1)=1
300  k(i,3)=ie(1)
301  k(i,4)=0
302  k(i,5)=0
303 
304 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
305  330 CALL lukfdi(kfl(1),0,kfl(3),k(i,2))
306  IF(k(i,2).EQ.0) goto 270
307  IF(mstj(12).GE.3.AND.irankj.EQ.1.AND.iabs(kfl(1)).LE.10.AND.
308  &iabs(kfl(3)).GT.10) THEN
309  IF(rlu(0).GT.parj(19)) goto 330
310  ENDIF
311  p(i,5)=ulmass(k(i,2))
312  CALL luptdi(kfl(1),px(3),py(3))
313  pr(1)=p(i,5)**2+(px(1)+px(3))**2+(py(1)+py(3))**2
314  CALL luzdis(kfl(1),kfl(3),pr(1),z)
315  gam(3)=(1.-z)*(gam(1)+pr(1)/z)
316  DO 340 j=1,3
317  340 in(j)=in(3+j)
318 
319 C...Junction strings: stepping within or from 'low' string region easy.
320  IF(in(1)+1.EQ.in(2).AND.z*p(in(1)+2,3)*p(in(2)+2,3)*
321  &p(in(1),5)**2.GE.pr(1)) THEN
322  p(in(1)+2,4)=z*p(in(1)+2,3)
323  p(in(2)+2,4)=pr(1)/(p(in(1)+2,4)*p(in(1),5)**2)
324  DO 350 j=1,4
325  350 p(i,j)=(px(1)+px(3))*p(in(3),j)+(py(1)+py(3))*p(in(3)+1,j)
326  goto 420
327  ELSEIF(in(1)+1.EQ.in(2)) THEN
328  p(in(2)+2,4)=p(in(2)+2,3)
329  p(in(2)+2,1)=1.
330  in(2)=in(2)+4
331  IF(in(2).GT.n+nr+4*ns) goto 270
332  IF(four(in(1),in(2)).LE.1e-2) THEN
333  p(in(1)+2,4)=p(in(1)+2,3)
334  p(in(1)+2,1)=0.
335  in(1)=in(1)+4
336  ENDIF
337  ENDIF
338 
339 C...Junction strings: find new transverse directions.
340  360 IF(in(1).GT.n+nr+4*ns.OR.in(2).GT.n+nr+4*ns.OR.
341  &in(1).GT.in(2)) goto 270
342  IF(in(1).NE.in(4).OR.in(2).NE.in(5)) THEN
343  DO 370 j=1,4
344  dp(1,j)=p(in(1),j)
345  dp(2,j)=p(in(2),j)
346  dp(3,j)=0.
347  370 dp(4,j)=0.
348  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
349  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
350  dhc12=dfour(1,2)
351  IF(dhc12.LE.1e-2) THEN
352  p(in(1)+2,4)=p(in(1)+2,3)
353  p(in(1)+2,1)=0.
354  in(1)=in(1)+4
355  goto 360
356  ENDIF
357  in(3)=n+nr+4*ns+5
358  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
359  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
360  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
361  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
362  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
363  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
364  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
365  dhcx1=dfour(3,1)/dhc12
366  dhcx2=dfour(3,2)/dhc12
367  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
368  dhcy1=dfour(4,1)/dhc12
369  dhcy2=dfour(4,2)/dhc12
370  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
371  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
372  DO 380 j=1,4
373  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
374  p(in(3),j)=dp(3,j)
375  380 p(in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
376  & dhcyx*dp(3,j))
377 C...Express pT with respect to new axes, if sensible.
378  pxp=-(px(3)*four(in(6),in(3))+py(3)*four(in(6)+1,in(3)))
379  pyp=-(px(3)*four(in(6),in(3)+1)+py(3)*four(in(6)+1,in(3)+1))
380  IF(abs(pxp**2+pyp**2-px(3)**2-py(3)**2).LT.0.01) THEN
381  px(3)=pxp
382  py(3)=pyp
383  ENDIF
384  ENDIF
385 
386 C...Junction strings: sum up known four-momentum, coefficients for m2.
387  DO 400 j=1,4
388  dhg(j)=0.
389  p(i,j)=px(1)*p(in(6),j)+py(1)*p(in(6)+1,j)+px(3)*p(in(3),j)+
390  &py(3)*p(in(3)+1,j)
391  DO 390 in1=in(4),in(1)-4,4
392  390 p(i,j)=p(i,j)+p(in1+2,3)*p(in1,j)
393  DO 400 in2=in(5),in(2)-4,4
394  400 p(i,j)=p(i,j)+p(in2+2,3)*p(in2,j)
395  dhm(1)=four(i,i)
396  dhm(2)=2.*four(i,in(1))
397  dhm(3)=2.*four(i,in(2))
398  dhm(4)=2.*four(in(1),in(2))
399 
400 C...Junction strings: find coefficients for Gamma expression.
401  DO 410 in2=in(1)+1,in(2),4
402  DO 410 in1=in(1),in2-1,4
403  dhc=2.*four(in1,in2)
404  dhg(1)=dhg(1)+p(in1+2,1)*p(in2+2,1)*dhc
405  IF(in1.EQ.in(1)) dhg(2)=dhg(2)-p(in2+2,1)*dhc
406  IF(in2.EQ.in(2)) dhg(3)=dhg(3)+p(in1+2,1)*dhc
407  410 IF(in1.EQ.in(1).AND.in2.EQ.in(2)) dhg(4)=dhg(4)-dhc
408 
409 C...Junction strings: solve (m2, Gamma) equation system for energies.
410  dhs1=dhm(3)*dhg(4)-dhm(4)*dhg(3)
411  IF(abs(dhs1).LT.1e-4) goto 270
412  dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(2)*dhg(3)-dhg(4)*
413  &(p(i,5)**2-dhm(1))+dhg(2)*dhm(3)
414  dhs3=dhm(2)*(gam(3)-dhg(1))-dhg(2)*(p(i,5)**2-dhm(1))
415  p(in(2)+2,4)=0.5*(sqrt(max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
416  &dhs2/dhs1)
417  IF(dhm(2)+dhm(4)*p(in(2)+2,4).LE.0.) goto 270
418  p(in(1)+2,4)=(p(i,5)**2-dhm(1)-dhm(3)*p(in(2)+2,4))/
419  &(dhm(2)+dhm(4)*p(in(2)+2,4))
420 
421 C...Junction strings: step to new region if necessary.
422  IF(p(in(2)+2,4).GT.p(in(2)+2,3)) THEN
423  p(in(2)+2,4)=p(in(2)+2,3)
424  p(in(2)+2,1)=1.
425  in(2)=in(2)+4
426  IF(in(2).GT.n+nr+4*ns) goto 270
427  IF(four(in(1),in(2)).LE.1e-2) THEN
428  p(in(1)+2,4)=p(in(1)+2,3)
429  p(in(1)+2,1)=0.
430  in(1)=in(1)+4
431  ENDIF
432  goto 360
433  ELSEIF(p(in(1)+2,4).GT.p(in(1)+2,3)) THEN
434  p(in(1)+2,4)=p(in(1)+2,3)
435  p(in(1)+2,1)=0.
436  in(1)=in(1)+js
437  goto 710
438  ENDIF
439 
440 C...Junction strings: particle four-momentum, remainder, loop back.
441  420 DO 430 j=1,4
442  p(i,j)=p(i,j)+p(in(1)+2,4)*p(in(1),j)+p(in(2)+2,4)*p(in(2),j)
443  430 pju(iu+3,j)=pju(iu+3,j)+p(i,j)
444  IF(p(i,4).LE.0.) goto 270
445  pju(iu+3,5)=tju(4)*pju(iu+3,4)-tju(1)*pju(iu+3,1)-
446  &tju(2)*pju(iu+3,2)-tju(3)*pju(iu+3,3)
447  IF(pju(iu+3,5).LT.pju(iu,5)) THEN
448  kfl(1)=-kfl(3)
449  px(1)=-px(3)
450  py(1)=-py(3)
451  gam(1)=gam(3)
452  IF(in(3).NE.in(6)) THEN
453  DO 440 j=1,4
454  p(in(6),j)=p(in(3),j)
455  440 p(in(6)+1,j)=p(in(3)+1,j)
456  ENDIF
457  DO 450 jq=1,2
458  in(3+jq)=in(jq)
459  p(in(jq)+2,3)=p(in(jq)+2,3)-p(in(jq)+2,4)
460  450 p(in(jq)+2,1)=p(in(jq)+2,1)-(3-2*jq)*p(in(jq)+2,4)
461  goto 320
462  ENDIF
463 
464 C...Junction strings: save quantities left after each string.
465  IF(iabs(kfl(1)).GT.10) goto 270
466  i=i-1
467  kfjh(iu)=kfl(1)
468  DO 460 j=1,4
469  460 pju(iu+3,j)=pju(iu+3,j)-p(i+1,j)
470  470 CONTINUE
471 
472 C...Junction strings: put together to new effective string endpoint.
473  njs(jt)=i-ista
474  kfjs(jt)=k(k(mju(jt+2),3),2)
475  kfls=2*int(rlu(0)+3.*parj(4)/(1.+3.*parj(4)))+1
476  IF(kfjh(1).EQ.kfjh(2)) kfls=3
477  IF(ista.NE.i) kfjs(jt)=isign(1000*max(iabs(kfjh(1)),
478  &iabs(kfjh(2)))+100*min(iabs(kfjh(1)),iabs(kfjh(2)))+
479  &kfls,kfjh(1))
480  DO 480 j=1,4
481  pjs(jt,j)=pju(1,j)+pju(2,j)+p(mju(jt),j)
482  480 pjs(jt+2,j)=pju(4,j)+pju(5,j)
483  pjs(jt,5)=sqrt(max(0.,pjs(jt,4)**2-pjs(jt,1)**2-pjs(jt,2)**2-
484  &pjs(jt,3)**2))
485  490 CONTINUE
486 
487 C...Open versus closed strings. Choose breakup region for latter.
488  500 IF(mju(1).NE.0.AND.mju(2).NE.0) THEN
489  ns=mju(2)-mju(1)
490  nb=mju(1)-n
491  ELSEIF(mju(1).NE.0) THEN
492  ns=n+nr-mju(1)
493  nb=mju(1)-n
494  ELSEIF(mju(2).NE.0) THEN
495  ns=mju(2)-n
496  nb=1
497  ELSEIF(iabs(k(n+1,2)).NE.21) THEN
498  ns=nr-1
499  nb=1
500  ELSE
501  ns=nr+1
502  w2sum=0.
503  DO 510 is=1,nr
504  p(n+nr+is,1)=0.5*four(n+is,n+is+1-nr*(is/nr))
505  510 w2sum=w2sum+p(n+nr+is,1)
506  w2ran=rlu(0)*w2sum
507  nb=0
508  520 nb=nb+1
509  w2sum=w2sum-p(n+nr+nb,1)
510  IF(w2sum.GT.w2ran.AND.nb.LT.nr) goto 520
511  ENDIF
512 
513 C...Find longitudinal string directions (i.e. lightlike four-vectors).
514  DO 540 is=1,ns
515  is1=n+is+nb-1-nr*((is+nb-2)/nr)
516  is2=n+is+nb-nr*((is+nb-1)/nr)
517  DO 530 j=1,5
518  dp(1,j)=p(is1,j)
519  IF(iabs(k(is1,2)).EQ.21) dp(1,j)=0.5*dp(1,j)
520  IF(is1.EQ.mju(1)) dp(1,j)=pjs(1,j)-pjs(3,j)
521  dp(2,j)=p(is2,j)
522  IF(iabs(k(is2,2)).EQ.21) dp(2,j)=0.5*dp(2,j)
523  530 IF(is2.EQ.mju(2)) dp(2,j)=pjs(2,j)-pjs(4,j)
524  dp(3,5)=dfour(1,1)
525  dp(4,5)=dfour(2,2)
526  dhkc=dfour(1,2)
527  IF(dp(3,5)+2.*dhkc+dp(4,5).LE.0.) THEN
528  dp(3,5)=dp(1,5)**2
529  dp(4,5)=dp(2,5)**2
530  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2+dp(1,5)**2)
531  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2+dp(2,5)**2)
532  dhkc=dfour(1,2)
533  ENDIF
534  dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
535  dhk1=0.5*((dp(4,5)+dhkc)/dhks-1.)
536  dhk2=0.5*((dp(3,5)+dhkc)/dhks-1.)
537  in1=n+nr+4*is-3
538  p(in1,5)=sqrt(dp(3,5)+2.*dhkc+dp(4,5))
539  DO 540 j=1,4
540  p(in1,j)=(1.+dhk1)*dp(1,j)-dhk2*dp(2,j)
541  540 p(in1+1,j)=(1.+dhk2)*dp(2,j)-dhk1*dp(1,j)
542 
543 C...Begin initialization: sum up energy, set starting position.
544  isav=i
545  550 ntry=ntry+1
546  IF(ntry.GT.100.AND.ntryr.LE.4) THEN
547  paru12=4.*paru12
548  paru13=2.*paru13
549  goto 130
550  ELSEIF(ntry.GT.100) THEN
551  CALL luerrm(14,'(LUSTRF:) caught in infinite loop')
552  IF(mstu(21).GE.1) RETURN
553  ENDIF
554  i=isav
555  DO 560 j=1,4
556  p(n+nrs,j)=0.
557  DO 560 is=1,nr
558  560 p(n+nrs,j)=p(n+nrs,j)+p(n+is,j)
559  DO 570 jt=1,2
560  irank(jt)=0
561  IF(mju(jt).NE.0) irank(jt)=njs(jt)
562  IF(ns.GT.nr) irank(jt)=1
563  ie(jt)=k(n+1+(jt/2)*(np-1),3)
564  in(3*jt+1)=n+nr+1+4*(jt/2)*(ns-1)
565  in(3*jt+2)=in(3*jt+1)+1
566  in(3*jt+3)=n+nr+4*ns+2*jt-1
567  DO 570 in1=n+nr+2+jt,n+nr+4*ns-2+jt,4
568  p(in1,1)=2-jt
569  p(in1,2)=jt-1
570  570 p(in1,3)=1.
571 
572 C...Initialize flavour and pT variables for open string.
573  IF(ns.LT.nr) THEN
574  px(1)=0.
575  py(1)=0.
576  IF(ns.EQ.1.AND.mju(1)+mju(2).EQ.0) CALL luptdi(0,px(1),py(1))
577  px(2)=-px(1)
578  py(2)=-py(1)
579  DO 580 jt=1,2
580  kfl(jt)=k(ie(jt),2)
581  IF(mju(jt).NE.0) kfl(jt)=kfjs(jt)
582  mstj(93)=1
583  pmq(jt)=ulmass(kfl(jt))
584  580 gam(jt)=0.
585 
586 C...Closed string: random initial breakup flavour, pT and vertex.
587  ELSE
588  kfl(3)=int(1.+(2.+parj(2))*rlu(0))*(-1)**int(rlu(0)+0.5)
589  CALL lukfdi(kfl(3),0,kfl(1),kdump)
590  kfl(2)=-kfl(1)
591  IF(iabs(kfl(1)).GT.10.AND.rlu(0).GT.0.5) THEN
592  kfl(2)=-(kfl(1)+isign(10000,kfl(1)))
593  ELSEIF(iabs(kfl(1)).GT.10) THEN
594  kfl(1)=-(kfl(2)+isign(10000,kfl(2)))
595  ENDIF
596  CALL luptdi(kfl(1),px(1),py(1))
597  px(2)=-px(1)
598  py(2)=-py(1)
599  pr3=min(25.,0.1*p(n+nr+1,5)**2)
600  590 CALL luzdis(kfl(1),kfl(2),pr3,z)
601  zr=pr3/(z*p(n+nr+1,5)**2)
602  IF(zr.GE.1.) goto 590
603  DO 600 jt=1,2
604  mstj(93)=1
605  pmq(jt)=ulmass(kfl(jt))
606  gam(jt)=pr3*(1.-z)/z
607  in1=n+nr+3+4*(jt/2)*(ns-1)
608  p(in1,jt)=1.-z
609  p(in1,3-jt)=jt-1
610  p(in1,3)=(2-jt)*(1.-z)+(jt-1)*z
611  p(in1+1,jt)=zr
612  p(in1+1,3-jt)=2-jt
613  600 p(in1+1,3)=(2-jt)*(1.-zr)+(jt-1)*zr
614  ENDIF
615 
616 C...Find initial transverse directions (i.e. spacelike four-vectors).
617  DO 640 jt=1,2
618  IF(jt.EQ.1.OR.ns.EQ.nr-1) THEN
619  in1=in(3*jt+1)
620  in3=in(3*jt+3)
621  DO 610 j=1,4
622  dp(1,j)=p(in1,j)
623  dp(2,j)=p(in1+1,j)
624  dp(3,j)=0.
625  610 dp(4,j)=0.
626  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
627  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
628  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
629  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
630  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
631  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
632  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
633  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
634  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
635  dhc12=dfour(1,2)
636  dhcx1=dfour(3,1)/dhc12
637  dhcx2=dfour(3,2)/dhc12
638  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
639  dhcy1=dfour(4,1)/dhc12
640  dhcy2=dfour(4,2)/dhc12
641  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
642  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
643  DO 620 j=1,4
644  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
645  p(in3,j)=dp(3,j)
646  620 p(in3+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
647  & dhcyx*dp(3,j))
648  ELSE
649  DO 630 j=1,4
650  p(in3+2,j)=p(in3,j)
651  630 p(in3+3,j)=p(in3+1,j)
652  ENDIF
653  640 CONTINUE
654 
655 C...Remove energy used up in junction string fragmentation.
656  IF(mju(1)+mju(2).GT.0) THEN
657  DO 660 jt=1,2
658  IF(njs(jt).EQ.0) goto 660
659  DO 650 j=1,4
660  650 p(n+nrs,j)=p(n+nrs,j)-pjs(jt+2,j)
661  660 CONTINUE
662  ENDIF
663 
664 C...Produce new particle: side, origin.
665  670 i=i+1
666  IF(2*i-nsav.GE.mstu(4)-mstu(32)-5) THEN
667  CALL luerrm(11,'(LUSTRF:) no more memory left in LUJETS')
668  IF(mstu(21).GE.1) RETURN
669  ENDIF
670  jt=1.5+rlu(0)
671  IF(iabs(kfl(3-jt)).GT.10) jt=3-jt
672  jr=3-jt
673  js=3-2*jt
674  irank(jt)=irank(jt)+1
675  k(i,1)=1
676  k(i,3)=ie(jt)
677  k(i,4)=0
678  k(i,5)=0
679 
680 C...Generate flavour, hadron and pT.
681  680 CALL lukfdi(kfl(jt),0,kfl(3),k(i,2))
682  IF(k(i,2).EQ.0) goto 550
683  IF(mstj(12).GE.3.AND.irank(jt).EQ.1.AND.iabs(kfl(jt)).LE.10.AND.
684  &iabs(kfl(3)).GT.10) THEN
685  IF(rlu(0).GT.parj(19)) goto 680
686  ENDIF
687  p(i,5)=ulmass(k(i,2))
688  CALL luptdi(kfl(jt),px(3),py(3))
689  pr(jt)=p(i,5)**2+(px(jt)+px(3))**2+(py(jt)+py(3))**2
690 
691 C...Final hadrons for small invariant mass.
692  mstj(93)=1
693  pmq(3)=ulmass(kfl(3))
694  wmin=parj(32+mstj(11))+pmq(1)+pmq(2)+parj(36)*pmq(3)
695  IF(iabs(kfl(jt)).GT.10.AND.iabs(kfl(3)).GT.10) wmin=
696  &wmin-0.5*parj(36)*pmq(3)
697  wrem2=four(n+nrs,n+nrs)
698  IF(wrem2.LT.0.10) goto 550
699  IF(wrem2.LT.max(wmin*(1.+(2.*rlu(0)-1.)*parj(37)),
700  &parj(32)+pmq(1)+pmq(2))**2) goto 810
701 
702 C...Choose z, which gives Gamma. Shift z for heavy flavours.
703  CALL luzdis(kfl(jt),kfl(3),pr(jt),z)
704  kfl1a=iabs(kfl(1))
705  kfl2a=iabs(kfl(2))
706  IF(max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
707  &mod(kfl2a/1000,10)).GE.4) THEN
708  pr(jr)=(pmq(jr)+pmq(3))**2+(px(jr)-px(3))**2+(py(jr)-py(3))**2
709  pw12=sqrt(max(0.,(wrem2-pr(1)-pr(2))**2-4.*pr(1)*pr(2)))
710  z=(wrem2+pr(jt)-pr(jr)+pw12*(2.*z-1.))/(2.*wrem2)
711  pr(jr)=(pmq(jr)+parj(32+mstj(11)))**2+(px(jr)-px(3))**2+
712  & (py(jr)-py(3))**2
713  IF((1.-z)*(wrem2-pr(jt)/z).LT.pr(jr)) goto 810
714  ENDIF
715  gam(3)=(1.-z)*(gam(jt)+pr(jt)/z)
716  DO 690 j=1,3
717  690 in(j)=in(3*jt+j)
718 
719 C...Stepping within or from 'low' string region easy.
720  IF(in(1)+1.EQ.in(2).AND.z*p(in(1)+2,3)*p(in(2)+2,3)*
721  &p(in(1),5)**2.GE.pr(jt)) THEN
722  p(in(jt)+2,4)=z*p(in(jt)+2,3)
723  p(in(jr)+2,4)=pr(jt)/(p(in(jt)+2,4)*p(in(1),5)**2)
724  DO 700 j=1,4
725  700 p(i,j)=(px(jt)+px(3))*p(in(3),j)+(py(jt)+py(3))*p(in(3)+1,j)
726  goto 770
727  ELSEIF(in(1)+1.EQ.in(2)) THEN
728  p(in(jr)+2,4)=p(in(jr)+2,3)
729  p(in(jr)+2,jt)=1.
730  in(jr)=in(jr)+4*js
731  IF(js*in(jr).GT.js*in(4*jr)) goto 550
732  IF(four(in(1),in(2)).LE.1e-2) THEN
733  p(in(jt)+2,4)=p(in(jt)+2,3)
734  p(in(jt)+2,jt)=0.
735  in(jt)=in(jt)+4*js
736  ENDIF
737  ENDIF
738 
739 C...Find new transverse directions (i.e. spacelike string vectors).
740  710 IF(js*in(1).GT.js*in(3*jr+1).OR.js*in(2).GT.js*in(3*jr+2).OR.
741  &in(1).GT.in(2)) goto 550
742  IF(in(1).NE.in(3*jt+1).OR.in(2).NE.in(3*jt+2)) THEN
743  DO 720 j=1,4
744  dp(1,j)=p(in(1),j)
745  dp(2,j)=p(in(2),j)
746  dp(3,j)=0.
747  720 dp(4,j)=0.
748  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
749  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
750  dhc12=dfour(1,2)
751  IF(dhc12.LE.1e-2) THEN
752  p(in(jt)+2,4)=p(in(jt)+2,3)
753  p(in(jt)+2,jt)=0.
754  in(jt)=in(jt)+4*js
755  goto 710
756  ENDIF
757  in(3)=n+nr+4*ns+5
758  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
759  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
760  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
761  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1.
762  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1.
763  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1.
764  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1.
765  dhcx1=dfour(3,1)/dhc12
766  dhcx2=dfour(3,2)/dhc12
767  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
768  dhcy1=dfour(4,1)/dhc12
769  dhcy2=dfour(4,2)/dhc12
770  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
771  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
772  DO 730 j=1,4
773  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
774  p(in(3),j)=dp(3,j)
775  730 p(in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
776  & dhcyx*dp(3,j))
777 C...Express pT with respect to new axes, if sensible.
778  pxp=-(px(3)*four(in(3*jt+3),in(3))+py(3)*
779  & four(in(3*jt+3)+1,in(3)))
780  pyp=-(px(3)*four(in(3*jt+3),in(3)+1)+py(3)*
781  & four(in(3*jt+3)+1,in(3)+1))
782  IF(abs(pxp**2+pyp**2-px(3)**2-py(3)**2).LT.0.01) THEN
783  px(3)=pxp
784  py(3)=pyp
785  ENDIF
786  ENDIF
787 
788 C...Sum up known four-momentum. Gives coefficients for m2 expression.
789  DO 750 j=1,4
790  dhg(j)=0.
791  p(i,j)=px(jt)*p(in(3*jt+3),j)+py(jt)*p(in(3*jt+3)+1,j)+
792  &px(3)*p(in(3),j)+py(3)*p(in(3)+1,j)
793  DO 740 in1=in(3*jt+1),in(1)-4*js,4*js
794  740 p(i,j)=p(i,j)+p(in1+2,3)*p(in1,j)
795  DO 750 in2=in(3*jt+2),in(2)-4*js,4*js
796  750 p(i,j)=p(i,j)+p(in2+2,3)*p(in2,j)
797  dhm(1)=four(i,i)
798  dhm(2)=2.*four(i,in(1))
799  dhm(3)=2.*four(i,in(2))
800  dhm(4)=2.*four(in(1),in(2))
801 
802 C...Find coefficients for Gamma expression.
803  DO 760 in2=in(1)+1,in(2),4
804  DO 760 in1=in(1),in2-1,4
805  dhc=2.*four(in1,in2)
806  dhg(1)=dhg(1)+p(in1+2,jt)*p(in2+2,jt)*dhc
807  IF(in1.EQ.in(1)) dhg(2)=dhg(2)-js*p(in2+2,jt)*dhc
808  IF(in2.EQ.in(2)) dhg(3)=dhg(3)+js*p(in1+2,jt)*dhc
809  760 IF(in1.EQ.in(1).AND.in2.EQ.in(2)) dhg(4)=dhg(4)-dhc
810 
811 C...Solve (m2, Gamma) equation system for energies taken.
812  dhs1=dhm(jr+1)*dhg(4)-dhm(4)*dhg(jr+1)
813  IF(abs(dhs1).LT.1e-4) goto 550
814  dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(jt+1)*dhg(jr+1)-dhg(4)*
815  &(p(i,5)**2-dhm(1))+dhg(jt+1)*dhm(jr+1)
816  dhs3=dhm(jt+1)*(gam(3)-dhg(1))-dhg(jt+1)*(p(i,5)**2-dhm(1))
817  p(in(jr)+2,4)=0.5*(sqrt(max(0d0,dhs2**2-4.*dhs1*dhs3))/abs(dhs1)-
818  &dhs2/dhs1)
819  IF(dhm(jt+1)+dhm(4)*p(in(jr)+2,4).LE.0.) goto 550
820  p(in(jt)+2,4)=(p(i,5)**2-dhm(1)-dhm(jr+1)*p(in(jr)+2,4))/
821  &(dhm(jt+1)+dhm(4)*p(in(jr)+2,4))
822 
823 C...Step to new region if necessary.
824  IF(p(in(jr)+2,4).GT.p(in(jr)+2,3)) THEN
825  p(in(jr)+2,4)=p(in(jr)+2,3)
826  p(in(jr)+2,jt)=1.
827  in(jr)=in(jr)+4*js
828  IF(js*in(jr).GT.js*in(4*jr)) goto 550
829  IF(four(in(1),in(2)).LE.1e-2) THEN
830  p(in(jt)+2,4)=p(in(jt)+2,3)
831  p(in(jt)+2,jt)=0.
832  in(jt)=in(jt)+4*js
833  ENDIF
834  goto 710
835  ELSEIF(p(in(jt)+2,4).GT.p(in(jt)+2,3)) THEN
836  p(in(jt)+2,4)=p(in(jt)+2,3)
837  p(in(jt)+2,jt)=0.
838  in(jt)=in(jt)+4*js
839  goto 710
840  ENDIF
841 
842 C...Four-momentum of particle. Remaining quantities. Loop back.
843  770 DO 780 j=1,4
844  p(i,j)=p(i,j)+p(in(1)+2,4)*p(in(1),j)+p(in(2)+2,4)*p(in(2),j)
845  780 p(n+nrs,j)=p(n+nrs,j)-p(i,j)
846  IF(p(i,4).LE.0.) goto 550
847  kfl(jt)=-kfl(3)
848  pmq(jt)=pmq(3)
849  px(jt)=-px(3)
850  py(jt)=-py(3)
851  gam(jt)=gam(3)
852  IF(in(3).NE.in(3*jt+3)) THEN
853  DO 790 j=1,4
854  p(in(3*jt+3),j)=p(in(3),j)
855  790 p(in(3*jt+3)+1,j)=p(in(3)+1,j)
856  ENDIF
857  DO 800 jq=1,2
858  in(3*jt+jq)=in(jq)
859  p(in(jq)+2,3)=p(in(jq)+2,3)-p(in(jq)+2,4)
860  800 p(in(jq)+2,jt)=p(in(jq)+2,jt)-js*(3-2*jq)*p(in(jq)+2,4)
861  goto 670
862 
863 C...Final hadron: side, flavour, hadron, mass.
864  810 i=i+1
865  k(i,1)=1
866  k(i,3)=ie(jr)
867  k(i,4)=0
868  k(i,5)=0
869  CALL lukfdi(kfl(jr),-kfl(3),kfldmp,k(i,2))
870  IF(k(i,2).EQ.0) goto 550
871  p(i,5)=ulmass(k(i,2))
872  pr(jr)=p(i,5)**2+(px(jr)-px(3))**2+(py(jr)-py(3))**2
873 
874 C...Final two hadrons: find common setup of four-vectors.
875  jq=1
876  IF(p(in(4)+2,3)*p(in(5)+2,3)*four(in(4),in(5)).LT.p(in(7),3)*
877  &p(in(8),3)*four(in(7),in(8))) jq=2
878  dhc12=four(in(3*jq+1),in(3*jq+2))
879  dhr1=four(n+nrs,in(3*jq+2))/dhc12
880  dhr2=four(n+nrs,in(3*jq+1))/dhc12
881  IF(in(4).NE.in(7).OR.in(5).NE.in(8)) THEN
882  px(3-jq)=-four(n+nrs,in(3*jq+3))-px(jq)
883  py(3-jq)=-four(n+nrs,in(3*jq+3)+1)-py(jq)
884  pr(3-jq)=p(i+(jt+jq-3)**2-1,5)**2+(px(3-jq)+(2*jq-3)*js*
885  & px(3))**2+(py(3-jq)+(2*jq-3)*js*py(3))**2
886  ENDIF
887 
888 C...Solve kinematics for final two hadrons, if possible.
889  wrem2=wrem2+(px(1)+px(2))**2+(py(1)+py(2))**2
890  fd=(sqrt(pr(1))+sqrt(pr(2)))/sqrt(wrem2)
891  IF(mju(1)+mju(2).NE.0.AND.i.EQ.isav+2.AND.fd.GE.1.) goto 180
892  IF(fd.GE.1.) goto 550
893  fa=wrem2+pr(jt)-pr(jr)
894  IF(mstj(11).EQ.2) prev=0.5*fd**parj(37+mstj(11))
895  IF(mstj(11).NE.2) prev=0.5*exp(max(-100.,log(fd)*
896  &parj(37+mstj(11))*(pr(1)+pr(2))**2))
897  fb=sign(sqrt(max(0.,fa**2-4.*wrem2*pr(jt))),js*(rlu(0)-prev))
898  kfl1a=iabs(kfl(1))
899  kfl2a=iabs(kfl(2))
900  IF(max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
901  &mod(kfl2a/1000,10)).GE.6) fb=sign(sqrt(max(0.,fa**2-
902  &4.*wrem2*pr(jt))),float(js))
903  DO 820 j=1,4
904  p(i-1,j)=(px(jt)+px(3))*p(in(3*jq+3),j)+(py(jt)+py(3))*
905  &p(in(3*jq+3)+1,j)+0.5*(dhr1*(fa+fb)*p(in(3*jq+1),j)+
906  &dhr2*(fa-fb)*p(in(3*jq+2),j))/wrem2
907  820 p(i,j)=p(n+nrs,j)-p(i-1,j)
908 
909 C...Mark jets as fragmented and give daughter pointers.
910  n=i-nrs+1
911  DO 830 i=nsav+1,nsav+np
912  im=k(i,3)
913  k(im,1)=k(im,1)+10
914  IF(mstu(16).NE.2) THEN
915  k(im,4)=nsav+1
916  k(im,5)=nsav+1
917  ELSE
918  k(im,4)=nsav+2
919  k(im,5)=n
920  ENDIF
921  830 CONTINUE
922 
923 C...Document string system. Move up particles.
924  nsav=nsav+1
925  k(nsav,1)=11
926  k(nsav,2)=92
927  k(nsav,3)=ip
928  k(nsav,4)=nsav+1
929  k(nsav,5)=n
930  DO 840 j=1,4
931  p(nsav,j)=dps(j)
932  840 v(nsav,j)=v(ip,j)
933  p(nsav,5)=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
934  v(nsav,5)=0.
935  DO 850 i=nsav+1,n
936  DO 850 j=1,5
937  k(i,j)=k(i+nrs-1,j)
938  p(i,j)=p(i+nrs-1,j)
939  850 v(i,j)=0.
940 
941 C...Order particles in rank along the chain. Update mother pointer.
942  DO 860 i=nsav+1,n
943  DO 860 j=1,5
944  k(i-nsav+n,j)=k(i,j)
945  860 p(i-nsav+n,j)=p(i,j)
946  i1=nsav
947  DO 880 i=n+1,2*n-nsav
948  IF(k(i,3).NE.ie(1)) goto 880
949  i1=i1+1
950  DO 870 j=1,5
951  k(i1,j)=k(i,j)
952  870 p(i1,j)=p(i,j)
953  IF(mstu(16).NE.2) k(i1,3)=nsav
954  880 CONTINUE
955  DO 900 i=2*n-nsav,n+1,-1
956  IF(k(i,3).EQ.ie(1)) goto 900
957  i1=i1+1
958  DO 890 j=1,5
959  k(i1,j)=k(i,j)
960  890 p(i1,j)=p(i,j)
961  IF(mstu(16).NE.2) k(i1,3)=nsav
962  900 CONTINUE
963 
964 C...Boost back particle system. Set production vertices.
965  CALL ludbrb(nsav+1,n,0.,0.,dps(1)/dps(4),dps(2)/dps(4),
966  &dps(3)/dps(4))
967  DO 910 i=nsav+1,n
968  DO 910 j=1,4
969  910 v(i,j)=v(ip,j)
970 
971  RETURN
972  END