Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pystrf.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pystrf.f
1 
2 C*********************************************************************
3 
4 C...PYSTRF
5 C...Handles the fragmentation of an arbitrary colour singlet
6 C...jet system according to the Lund string fragmentation model.
7 
8  SUBROUTINE pystrf(IP)
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...Commonblocks.
15  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
17  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18  SAVE /pyjets/,/pydat1/,/pydat2/
19 C...Local arrays. All MOPS variables ends with MO
20  dimension dps(5),kfl(3),pmq(3),px(3),py(3),gam(3),ie(2),pr(2),
21  &in(9),dhm(4),dhg(4),dp(5,5),irank(2),mju(4),iju(6),pju(5,5),
22  &tju(5),kfjh(2),njs(2),kfjs(2),pjs(4,5),mstu9t(8),paru9t(8),
23  &inmo(9),pm2qmo(2),xtmo(2),ejstr(2),ijuori(2),ibarrk(2),
24  &pbst(3,5),tjuold(5)
25 
26 C...Function: four-product of two vectors.
27  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)
28  dfour(i,j)=dp(i,4)*dp(j,4)-dp(i,1)*dp(j,1)-dp(i,2)*dp(j,2)-
29  &dp(i,3)*dp(j,3)
30 
31 C...Reset counters.
32  mstj(91)=0
33  nsav=n
34  mstu90=mstu(90)
35  np=0
36  kqsum=0
37  DO 100 j=1,5
38  dps(j)=0d0
39  100 CONTINUE
40  mju(1)=0
41  mju(2)=0
42  ntryfn=0
43  ijuori(1)=0
44  ijuori(2)=0
45 
46 C...Identify parton system.
47  i=ip-1
48  110 i=i+1
49  IF(i.GT.min(n,mstu(4)-mstu(32))) THEN
50  CALL pyerrm(12,'(PYSTRF:) failed to reconstruct jet system')
51  IF(mstu(21).GE.1) RETURN
52  ENDIF
53  IF(k(i,1).NE.1.AND.k(i,1).NE.2.AND.k(i,1).NE.41) goto 110
54  kc=pycomp(k(i,2))
55  IF(kc.EQ.0) goto 110
56  kq=kchg(kc,2)*isign(1,k(i,2))
57  IF(kq.EQ.0.AND.k(i,1).NE.41) goto 110
58  IF(n+5*np+11.GT.mstu(4)-mstu(32)-5) THEN
59  CALL pyerrm(11,'(PYSTRF:) no more memory left in PYJETS')
60  IF(mstu(21).GE.1) RETURN
61  ENDIF
62 
63 C...Take copy of partons to be considered. Check flavour sum.
64  np=np+1
65  DO 120 j=1,5
66  k(n+np,j)=k(i,j)
67  p(n+np,j)=p(i,j)
68  IF(j.NE.4) dps(j)=dps(j)+p(i,j)
69  120 CONTINUE
70  dps(4)=dps(4)+sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2+p(i,5)**2)
71  k(n+np,3)=i
72  IF(kq.NE.2) kqsum=kqsum+kq
73  IF(k(i,1).EQ.41) THEN
74  IF(mod(kqsum,2).EQ.0.AND.mju(1).EQ.0) THEN
75  mju(1)=n+np
76  ijuori(1)=i
77  ELSE
78  mju(2)=n+np
79  ijuori(2)=i
80  ENDIF
81  ENDIF
82  IF(k(i,1).EQ.2.OR.k(i,1).EQ.41) goto 110
83  IF(mod(kqsum,3).NE.0) THEN
84  CALL pyerrm(12,'(PYSTRF:) unphysical flavour combination')
85  IF(mstu(21).GE.1) RETURN
86  ENDIF
87  IF(mju(1).GT.0.OR.mju(2).GT.0) mstu(29)=1
88 
89 C...Boost copied system to CM frame (for better numerical precision).
90  IF(abs(dps(3)).LT.0.99d0*dps(4)) THEN
91  mbst=0
92  mstu(33)=1
93  CALL pyrobo(n+1,n+np,0d0,0d0,-dps(1)/dps(4),-dps(2)/dps(4),
94  & -dps(3)/dps(4))
95  ELSE
96  mbst=1
97  hhbz=sqrt(max(1d-6,dps(4)+dps(3))/max(1d-6,dps(4)-dps(3)))
98  DO 130 i=n+1,n+np
99  hhpmt=p(i,1)**2+p(i,2)**2+p(i,5)**2
100  IF(p(i,3).GT.0d0) THEN
101  hhpez=max(1d-10,(p(i,4)+p(i,3))/hhbz)
102  p(i,3)=0.5d0*(hhpez-hhpmt/hhpez)
103  p(i,4)=0.5d0*(hhpez+hhpmt/hhpez)
104  ELSE
105  hhpez=max(1d-10,(p(i,4)-p(i,3))*hhbz)
106  p(i,3)=-0.5d0*(hhpez-hhpmt/hhpez)
107  p(i,4)=0.5d0*(hhpez+hhpmt/hhpez)
108  ENDIF
109  130 CONTINUE
110  ENDIF
111 
112 C...Search for very nearby partons that may be recombined.
113  ntryr=0
114  ntrywr=0
115  paru12=paru(12)
116  paru13=paru(13)
117  mju(3)=mju(1)
118  mju(4)=mju(2)
119  nr=np
120  nrmin=2
121  IF(mju(1).GT.0) nrmin=nrmin+2
122  IF(mju(2).GT.0) nrmin=nrmin+2
123  140 IF(nr.GT.nrmin) THEN
124  pdrmin=2d0*paru12
125  DO 150 i=n+1,n+nr
126  IF(i.EQ.n+nr.AND.iabs(k(n+1,2)).NE.21) goto 150
127  i1=i+1
128  IF(i.EQ.n+nr) i1=n+1
129  IF(k(i,1).EQ.41.OR.k(i1,1).EQ.41) goto 150
130  IF(mju(1).NE.0.AND.i1.LT.mju(1).AND.iabs(k(i1,2)).NE.21)
131  & goto 150
132  IF(mju(2).NE.0.AND.i.GT.mju(2).AND.iabs(k(i,2)).NE.21)
133  & goto 150
134  pap=sqrt((p(i,1)**2+p(i,2)**2+p(i,3)**2)*(p(i1,1)**2+
135  & p(i1,2)**2+p(i1,3)**2))
136  pvp=p(i,1)*p(i1,1)+p(i,2)*p(i1,2)+p(i,3)*p(i1,3)
137  pdr=4d0*(pap-pvp)**2/max(1d-6,paru13**2*pap+2d0*(pap-pvp))
138  IF(pdr.LT.pdrmin) THEN
139  ir=i
140  pdrmin=pdr
141  ENDIF
142  150 CONTINUE
143 
144 C...Recombine very nearby partons to avoid machine precision problems.
145  IF(pdrmin.LT.paru12.AND.ir.EQ.n+nr) THEN
146  DO 160 j=1,4
147  p(n+1,j)=p(n+1,j)+p(n+nr,j)
148  160 CONTINUE
149  p(n+1,5)=sqrt(max(0d0,p(n+1,4)**2-p(n+1,1)**2-p(n+1,2)**2-
150  & p(n+1,3)**2))
151  nr=nr-1
152  goto 140
153  ELSEIF(pdrmin.LT.paru12) THEN
154  DO 170 j=1,4
155  p(ir,j)=p(ir,j)+p(ir+1,j)
156  170 CONTINUE
157  p(ir,5)=sqrt(max(0d0,p(ir,4)**2-p(ir,1)**2-p(ir,2)**2-
158  & p(ir,3)**2))
159  IF(mju(2).NE.0.AND.ir.GT.mju(2)) k(ir,2)=k(ir+1,2)
160  DO 190 i=ir+1,n+nr-1
161  k(i,1)=k(i+1,1)
162  k(i,2)=k(i+1,2)
163  DO 180 j=1,5
164  p(i,j)=p(i+1,j)
165  180 CONTINUE
166  190 CONTINUE
167  IF(ir.EQ.n+nr-1) k(ir,2)=k(n+nr,2)
168  nr=nr-1
169  IF(mju(1).GT.ir) mju(1)=mju(1)-1
170  IF(mju(2).GT.ir) mju(2)=mju(2)-1
171  goto 140
172  ENDIF
173  ENDIF
174  ntryr=ntryr+1
175 
176 C...Reset particle counter. Skip ahead if no junctions are present;
177 C...this is usually the case!
178  nrs=max(5*nr+11,np)
179  ntry=0
180  200 ntry=ntry+1
181  IF(ntry.GT.100.AND.ntryr.LE.8.AND.nr.GT.nrmin) THEN
182  paru12=4d0*paru12
183  paru13=2d0*paru13
184  goto 140
185  ELSEIF(ntry.GT.100.OR.ntryr.GT.100) THEN
186  CALL pyerrm(14,'(PYSTRF:) caught in infinite loop')
187  IF(mstu(21).GE.1) RETURN
188  ENDIF
189  i=n+nrs
190  mstu(90)=mstu90
191  IF(mju(1).EQ.0.AND.mju(2).EQ.0) goto 650
192  IF(mstj(12).GE.4) CALL pyerrm(29,'(PYSTRF:) sorry,'//
193  & ' junction strings not handled by MSTJ(12)>3 options')
194  DO 640 jt=1,2
195  njs(jt)=0
196  IF(mju(jt).EQ.0) goto 640
197  js=3-2*jt
198 
199 C++SKANDS
200 C...Find and sum up momentum on three sides of junction.
201 C...Begin with previous boost = zero.
202  ijrfit=0
203  DO 210 ix=1,3
204  tjuold(ix)=0d0
205  210 CONTINUE
206  tjuold(4)=1d0
207  220 iu=0
208 C...Beginning and end of string system in event record.
209  i1beg=n+1+(jt-1)*(nr-1)
210  i1end=n+nr+(jt-1)*(1-nr)
211 C...Look for junction string piece end points
212  DO 230 i1=i1beg,i1end,js
213  IF(k(i1,2).NE.21.AND.iu.LE.5.AND.ijrfit.EQ.0) THEN
214 C...Store junction string piece end points.
215 C 1-junction systems 2-junction systems
216 C IU : 1 2 3 4 1 2 3 4 5 6
217 C IJU(IU): q-g-g-q-g-g-j-g-q q-g-g-q-g-j-g-g-j-g-q-g-g-q
218  iu=iu+1
219  iju(iu)=i1
220  ENDIF
221 C...Sum over momenta, from junction outwards.
222  230 CONTINUE
223  DO 280 iu=1,3
224  pwt=0d0
225 C...Initialize junction drag and string piece 4-vectors.
226  DO 240 j=1,5
227  pbst(iu,j)=0d0
228  pju(iu,j)=0d0
229  240 CONTINUE
230 C...First two branches. Inwards out means opposite direction to JS.
231 C...(JS is 1 for JT=1, -1 for JT=2)
232  IF (iu.LT.3) THEN
233  i1a=iju(iu+1)-js
234  i1b=iju(iu)
235  idir=-js
236 C...Last branch (gq or gjgqgq). Direction now reversed.
237  ELSE
238  i1a=iju(iu)+js
239  i1b=i1end
240  idir=js
241  ENDIF
242  DO 270 i1=i1a,i1b,idir
243 C...Sum up momentum directions with exponential suppression
244 C...for use in finding junction rest frame below.
245  IF (k(i1,2).EQ.88) THEN
246 C...gjgqgq type system encountered. Use current PWT as start
247 C...for both strings.
248  pwtold=pwt
249  ELSE
250  IF (i1.EQ.iju(5)+idir) pwt=pwtold
251 C...Sum up string piece (boosted) 4-momenta.
252  DO 250 j=1,4
253  pju(iu,j)=pju(iu,j)+p(i1,j)
254  250 CONTINUE
255 C...Compute "junction drag" vectors from (boosted) 4-momenta (initial
256 C...boost is zero, see above). Skip parton if suppression factor large.
257  IF (pwt.GT.10d0) goto 270
258 C...Compute momentum in current frame:
259  tdp=tjuold(1)*p(i1,1)+tjuold(2)*p(i1,2)+tjuold(3)*p(i1,3)
260  bfc=tdp/(1d0+tjuold(4))+p(i1,4)
261  DO 260 j=1,3
262  ptmp=p(i1,j)+tjuold(j)*bfc
263  pbst(iu,j)=pbst(iu,j)+ptmp*exp(-pwt)
264  260 CONTINUE
265 C...Boosted energy
266  ptmp=tjuold(4)*p(i1,4)+tdp
267  pbst(iu,4)=pbst(iu,j)+ptmp*exp(-pwt)
268  pwt=pwt+ptmp/parj(48)
269  ENDIF
270  270 CONTINUE
271 C...Put |p| rather than m in 5th slot.
272  pbst(iu,5)=sqrt(pbst(iu,1)**2+pbst(iu,2)**2+pbst(iu,3)**2)
273  pju(iu,5)=sqrt(pju(iu,1)**2+pju(iu,2)**2+pju(iu,3)**2)
274  280 CONTINUE
275 
276 C...Calculate boost from present frame to next JRF candidate.
277  ijrfit=ijrfit+1
278  CALL pyjurf(pbst,tju)
279 
280 C...After some iterations do not take full step in new direction.
281  IF(ijrfit.GT.5) THEN
282  reduce=0.8d0**(ijrfit-5)
283  tju(1)=reduce*tju(1)
284  tju(2)=reduce*tju(2)
285  tju(3)=reduce*tju(3)
286  tju(4)=sqrt(1d0+tju(1)**2+tju(2)**2+tju(3)**2)
287  ENDIF
288 
289 C...Combine new boost (TJU) with old boost (TJUOLD)
290  tmp=tju(1)*tjuold(1)+tju(2)*tjuold(2)+tju(3)*tjuold(3)
291  DO 290 ix=1,3
292  tjuold(ix)=tju(ix)+tjuold(ix)*(tmp/(1d0+tjuold(4))+tju(4))
293  290 CONTINUE
294  tjuold(4)=sqrt(1d0+tjuold(1)**2+tjuold(2)**2+tjuold(3)**2)
295 
296 C...If last boost small, accept JRF, else iterate.
297 C...Also prevent possibility of infinite loop.
298  IF (abs((tju(4)-1d0)/tjuold(4)).GT.0.01d0.AND.
299  & ijrfit.LT.mstj(18)) THEN
300  goto 220
301  ELSEIF (ijrfit.GE.mstj(18)) THEN
302  CALL pyerrm(1,'(PYSTRF:) failed to converge on JRF')
303  ENDIF
304 
305 C...Now store total boost in TJU and change perception.
306 C...TJUOLD = boost vector from CM of string syst -> JRF. Henceforth,
307 C...TJU = junction motion vector in string CM, so the sign changes.
308  DO 300 j=1,3
309  tju(j)=-tjuold(j)
310  300 CONTINUE
311  tju(4)=sqrt(1d0+tju(1)**2+tju(2)**2+tju(3)**2)
312 
313 C--SKANDS
314 
315 C...Calculate string piece energies in junction rest frame.
316  DO 310 iu=1,3
317  pju(iu,5)=tju(4)*pju(iu,4)-tju(1)*pju(iu,1)-tju(2)*pju(iu,2)-
318  & tju(3)*pju(iu,3)
319  pbst(iu,5)=tju(4)*pbst(iu,4)-tju(1)*pbst(iu,1)-
320  & tju(2)*pbst(iu,2)-tju(3)*pbst(iu,3)
321  310 CONTINUE
322 
323 C...Start preparing for fragmentation of two strings from junction.
324  ista=i
325  ntryer=0
326  320 ntryer=ntryer+1
327  i=ista
328  DO 620 iu=1,2
329  ns=iabs(iju(iu+1)-iju(iu))
330 
331 C...Junction strings: find longitudinal string directions.
332  DO 350 is=1,ns
333  is1=iju(iu)+js*(is-1)
334  is2=iju(iu)+js*is
335  DO 330 j=1,5
336  dp(1,j)=0.5d0*p(is1,j)
337  IF(is.EQ.1) dp(1,j)=p(is1,j)
338  dp(2,j)=0.5d0*p(is2,j)
339  IF(is.EQ.ns) dp(2,j)=(-pbst(iu,j)+2d0*pbst(iu,5)*tju(j))*
340  & (pju(iu,5)/pbst(iu,5))
341  330 CONTINUE
342  IF(is.EQ.ns) dp(2,5)=sqrt(max(0d0,pju(iu,4)**2-
343  & pju(iu,1)**2-pju(iu,2)**2-pju(iu,3)**2))
344  dp(3,5)=dfour(1,1)
345  dp(4,5)=dfour(2,2)
346  dhkc=dfour(1,2)
347  IF(dp(3,5)+2d0*dhkc+dp(4,5).LE.0d0) THEN
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  dp(3,5)=0d0
351  dp(4,5)=0d0
352  dhkc=dfour(1,2)
353  ENDIF
354  dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
355  dhk1=0.5d0*((dp(4,5)+dhkc)/dhks-1d0)
356  dhk2=0.5d0*((dp(3,5)+dhkc)/dhks-1d0)
357  in1=n+nr+4*is-3
358  p(in1,5)=sqrt(dp(3,5)+2d0*dhkc+dp(4,5))
359  DO 340 j=1,4
360  p(in1,j)=(1d0+dhk1)*dp(1,j)-dhk2*dp(2,j)
361  p(in1+1,j)=(1d0+dhk2)*dp(2,j)-dhk1*dp(1,j)
362  340 CONTINUE
363  350 CONTINUE
364 
365 C...Junction strings: initialize flavour, momentum and starting pos.
366  isav=i
367  mstu91=mstu(90)
368  360 ntry=ntry+1
369  IF(ntry.GT.100.AND.ntryr.LE.8.AND.nr.GT.nrmin) THEN
370  paru12=4d0*paru12
371  paru13=2d0*paru13
372  goto 140
373  ELSEIF(ntry.GT.100) THEN
374  CALL pyerrm(14,'(PYSTRF:) caught in infinite loop')
375  IF(mstu(21).GE.1) RETURN
376  ENDIF
377  i=isav
378  mstu(90)=mstu91
379  irankj=0
380  ie(1)=k(n+1+(jt/2)*(np-1),3)
381  IF (mod(jt+iu,2).NE.0) THEN
382  ie(1)=k(iju(iu),3)
383  IF (np-nr.NE.0) THEN
384 C...If gluons have disappeared. Original IJU must be used.
385  it=ip
386  ne=1
387  370 it=it+1
388  IF (k(it,2).NE.21) THEN
389  ne=ne+1
390  ENDIF
391  IF (ne.EQ.iu+4*(jt-1)) THEN
392  ie(1)=it
393  ELSEIF (it.LE.ip+np) THEN
394  goto 370
395  ELSE
396  CALL pyerrm(14,'(PYSTRF:) '//
397  & 'Original IJU could not be reconstructed!')
398  ENDIF
399  ENDIF
400  ENDIF
401  in(4)=n+nr+1
402  in(5)=in(4)+1
403  in(6)=n+nr+4*ns+1
404  DO 390 jq=1,2
405  DO 380 in1=n+nr+2+jq,n+nr+4*ns-2+jq,4
406  p(in1,1)=2-jq
407  p(in1,2)=jq-1
408  p(in1,3)=1d0
409  380 CONTINUE
410  390 CONTINUE
411  kfl(1)=k(iju(iu),2)
412  px(1)=0d0
413  py(1)=0d0
414  gam(1)=0d0
415  DO 400 j=1,5
416  pju(iu+3,j)=0d0
417  400 CONTINUE
418 
419 C...Junction strings: find initial transverse directions.
420  DO 410 j=1,4
421  dp(1,j)=p(in(4),j)
422  dp(2,j)=p(in(4)+1,j)
423  dp(3,j)=0d0
424  dp(4,j)=0d0
425  410 CONTINUE
426  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
427  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
428  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
429  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
430  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
431  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1d0
432  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1d0
433  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1d0
434  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1d0
435  dhc12=dfour(1,2)
436  dhcx1=dfour(3,1)/dhc12
437  dhcx2=dfour(3,2)/dhc12
438  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
439  dhcy1=dfour(4,1)/dhc12
440  dhcy2=dfour(4,2)/dhc12
441  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
442  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
443  DO 420 j=1,4
444  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
445  p(in(6),j)=dp(3,j)
446  p(in(6)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
447  & dhcyx*dp(3,j))
448  420 CONTINUE
449 
450 C...Junction strings: produce new particle, origin.
451  430 i=i+1
452  IF(2*i-nsav.GE.mstu(4)-mstu(32)-5) THEN
453  CALL pyerrm(11,'(PYSTRF:) no more memory left in PYJETS')
454  IF(mstu(21).GE.1) RETURN
455  ENDIF
456  irankj=irankj+1
457  k(i,1)=1
458  k(i,3)=ie(1)
459  k(i,4)=0
460  k(i,5)=0
461 
462 C...Junction strings: generate flavour, hadron, pT, z and Gamma.
463  440 CALL pykfdi(kfl(1),0,kfl(3),k(i,2))
464  IF(k(i,2).EQ.0) goto 360
465  IF(irankj.EQ.1.AND.iabs(kfl(1)).LE.10.AND.
466  & iabs(kfl(3)).GT.10) THEN
467  IF(pyr(0).GT.parj(19)) goto 440
468  ENDIF
469  p(i,5)=pymass(k(i,2))
470  CALL pyptdi(kfl(1),px(3),py(3))
471  pr(1)=p(i,5)**2+(px(1)+px(3))**2+(py(1)+py(3))**2
472  CALL pyzdis(kfl(1),kfl(3),pr(1),z)
473  IF(iabs(kfl(1)).GE.4.AND.iabs(kfl(1)).LE.8.AND.
474  & mstu(90).LT.8) THEN
475  mstu(90)=mstu(90)+1
476  mstu(90+mstu(90))=i
477  paru(90+mstu(90))=z
478  ENDIF
479  gam(3)=(1d0-z)*(gam(1)+pr(1)/z)
480  DO 450 j=1,3
481  in(j)=in(3+j)
482  450 CONTINUE
483 
484 C...Junction strings: stepping within 'low' string region.
485  IF(in(1)+1.EQ.in(2).AND.z*p(in(1)+2,3)*p(in(2)+2,3)*
486  & p(in(1),5)**2.GE.pr(1)) THEN
487  p(in(1)+2,4)=z*p(in(1)+2,3)
488  p(in(2)+2,4)=pr(1)/(p(in(1)+2,4)*p(in(1),5)**2)
489  DO 460 j=1,4
490  p(i,j)=(px(1)+px(3))*p(in(3),j)+(py(1)+py(3))*p(in(3)+1,j)
491  460 CONTINUE
492  goto 560
493 C...Has used up energy of junction string, i.e. no more hadrons in it.
494  ELSEIF(in(1)+1.EQ.in(2).AND.in(1).EQ.n+nr+4*ns-3) THEN
495  DO 470 j=1,5
496  p(i,j)=0d0
497  470 CONTINUE
498  goto 600
499 C...Stepping from 'low' string region
500  ELSEIF(in(1)+1.EQ.in(2)) THEN
501  p(in(2)+2,4)=p(in(2)+2,3)
502  p(in(2)+2,1)=1d0
503  in(2)=in(2)+4
504  IF(in(2).GT.n+nr+4*ns) goto 360
505  IF(four(in(1),in(2)).LE.1d-2) THEN
506  p(in(1)+2,4)=p(in(1)+2,3)
507  p(in(1)+2,1)=0d0
508  in(1)=in(1)+4
509  ENDIF
510  ENDIF
511 
512 C...Junction strings: find new transverse directions.
513  480 IF(in(1).GT.n+nr+4*ns.OR.in(2).GT.n+nr+4*ns.OR.
514  & in(1).GT.in(2)) goto 360
515  IF(in(1).NE.in(4).OR.in(2).NE.in(5)) THEN
516  DO 490 j=1,4
517  dp(1,j)=p(in(1),j)
518  dp(2,j)=p(in(2),j)
519  dp(3,j)=0d0
520  dp(4,j)=0d0
521  490 CONTINUE
522  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
523  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
524  dhc12=dfour(1,2)
525  IF(dhc12.LE.1d-2) THEN
526  p(in(1)+2,4)=p(in(1)+2,3)
527  p(in(1)+2,1)=0d0
528  in(1)=in(1)+4
529  goto 480
530  ENDIF
531  in(3)=n+nr+4*ns+5
532  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
533  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
534  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
535  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1d0
536  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1d0
537  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1d0
538  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1d0
539  dhcx1=dfour(3,1)/dhc12
540  dhcx2=dfour(3,2)/dhc12
541  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
542  dhcy1=dfour(4,1)/dhc12
543  dhcy2=dfour(4,2)/dhc12
544  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
545  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
546  DO 500 j=1,4
547  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
548  p(in(3),j)=dp(3,j)
549  p(in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
550  & dhcyx*dp(3,j))
551  500 CONTINUE
552 C...Express pT with respect to new axes, if sensible.
553  pxp=-(px(3)*four(in(6),in(3))+py(3)*four(in(6)+1,in(3)))
554  pyp=-(px(3)*four(in(6),in(3)+1)+py(3)*four(in(6)+1,in(3)+1))
555  IF(abs(pxp**2+pyp**2-px(3)**2-py(3)**2).LT.0.01d0) THEN
556  px(3)=pxp
557  py(3)=pyp
558  ENDIF
559  ENDIF
560 
561 C...Junction strings: sum up known four-momentum, coefficients for m2.
562  DO 530 j=1,4
563  dhg(j)=0d0
564  p(i,j)=px(1)*p(in(6),j)+py(1)*p(in(6)+1,j)+px(3)*p(in(3),j)+
565  & py(3)*p(in(3)+1,j)
566  DO 510 in1=in(4),in(1)-4,4
567  p(i,j)=p(i,j)+p(in1+2,3)*p(in1,j)
568  510 CONTINUE
569  DO 520 in2=in(5),in(2)-4,4
570  p(i,j)=p(i,j)+p(in2+2,3)*p(in2,j)
571  520 CONTINUE
572  530 CONTINUE
573  dhm(1)=four(i,i)
574  dhm(2)=2d0*four(i,in(1))
575  dhm(3)=2d0*four(i,in(2))
576  dhm(4)=2d0*four(in(1),in(2))
577 
578 C...Junction strings: find coefficients for Gamma expression.
579  DO 550 in2=in(1)+1,in(2),4
580  DO 540 in1=in(1),in2-1,4
581  dhc=2d0*four(in1,in2)
582  dhg(1)=dhg(1)+p(in1+2,1)*p(in2+2,1)*dhc
583  IF(in1.EQ.in(1)) dhg(2)=dhg(2)-p(in2+2,1)*dhc
584  IF(in2.EQ.in(2)) dhg(3)=dhg(3)+p(in1+2,1)*dhc
585  IF(in1.EQ.in(1).AND.in2.EQ.in(2)) dhg(4)=dhg(4)-dhc
586  540 CONTINUE
587  550 CONTINUE
588 
589 C...Junction strings: solve (m2, Gamma) equation system for energies.
590  dhs1=dhm(3)*dhg(4)-dhm(4)*dhg(3)
591  IF(abs(dhs1).LT.1d-4) goto 360
592  dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(2)*dhg(3)-dhg(4)*
593  & (p(i,5)**2-dhm(1))+dhg(2)*dhm(3)
594  dhs3=dhm(2)*(gam(3)-dhg(1))-dhg(2)*(p(i,5)**2-dhm(1))
595  p(in(2)+2,4)=0.5d0*(sqrt(max(0d0,dhs2**2-4d0*dhs1*dhs3))/
596  & abs(dhs1)-dhs2/dhs1)
597  IF(dhm(2)+dhm(4)*p(in(2)+2,4).LE.0d0) goto 360
598  p(in(1)+2,4)=(p(i,5)**2-dhm(1)-dhm(3)*p(in(2)+2,4))/
599  & (dhm(2)+dhm(4)*p(in(2)+2,4))
600 
601 C...Junction strings: step to new region if necessary.
602  IF(p(in(2)+2,4).GT.p(in(2)+2,3)) THEN
603  p(in(2)+2,4)=p(in(2)+2,3)
604  p(in(2)+2,1)=1d0
605  in(2)=in(2)+4
606  IF(in(2).GT.n+nr+4*ns) goto 360
607  IF(four(in(1),in(2)).LE.1d-2) THEN
608  p(in(1)+2,4)=p(in(1)+2,3)
609  p(in(1)+2,1)=0d0
610  in(1)=in(1)+4
611  ENDIF
612  goto 480
613  ELSEIF(p(in(1)+2,4).GT.p(in(1)+2,3)) THEN
614  p(in(1)+2,4)=p(in(1)+2,3)
615  p(in(1)+2,1)=0d0
616  in(1)=in(1)+4
617  goto 480
618  ENDIF
619 
620 C...Junction strings: particle four-momentum, remainder, loop back.
621  560 DO 570 j=1,4
622  p(i,j)=p(i,j)+p(in(1)+2,4)*p(in(1),j)+
623  & p(in(2)+2,4)*p(in(2),j)
624  pju(iu+3,j)=pju(iu+3,j)+p(i,j)
625  570 CONTINUE
626  IF(p(i,4).LT.p(i,5)) goto 360
627  pju(iu+3,5)=tju(4)*pju(iu+3,4)-tju(1)*pju(iu+3,1)-
628  & tju(2)*pju(iu+3,2)-tju(3)*pju(iu+3,3)
629  IF(pju(iu+3,5).LT.pju(iu,5)) THEN
630  kfl(1)=-kfl(3)
631  px(1)=-px(3)
632  py(1)=-py(3)
633  gam(1)=gam(3)
634  IF(in(3).NE.in(6)) THEN
635  DO 580 j=1,4
636  p(in(6),j)=p(in(3),j)
637  p(in(6)+1,j)=p(in(3)+1,j)
638  580 CONTINUE
639  ENDIF
640  DO 590 jq=1,2
641  in(3+jq)=in(jq)
642  p(in(jq)+2,3)=p(in(jq)+2,3)-p(in(jq)+2,4)
643  p(in(jq)+2,1)=p(in(jq)+2,1)-(3-2*jq)*p(in(jq)+2,4)
644  590 CONTINUE
645  goto 430
646  ENDIF
647 
648 C...Junction strings: save quantities left after each string.
649  IF(iabs(kfl(1)).GT.10) goto 360
650  600 i=i-1
651  kfjh(iu)=kfl(1)
652  DO 610 j=1,4
653  pju(iu+3,j)=pju(iu+3,j)-p(i+1,j)
654  610 CONTINUE
655 
656 C...Junction strings: loopback if much unused energy in both strings.
657  pju(iu+3,5)=tju(4)*pju(iu+3,4)-tju(1)*pju(iu+3,1)-
658  & tju(2)*pju(iu+3,2)-tju(3)*pju(iu+3,3)
659  ejstr(iu)=pju(iu,5)-pju(iu+3,5)
660  620 CONTINUE
661  IF((min(ejstr(1),ejstr(2)).GT.parj(49).OR.
662  & ejstr(1).GT.parj(49)+pyr(0)*parj(50).OR.
663  & ejstr(2).GT.parj(49)+pyr(0)*parj(50))
664  & .AND.ntryer.LT.10) goto 320
665 
666 C...Junction strings: put together to new effective string endpoint.
667  njs(jt)=i-ista
668  kfls=2*int(pyr(0)+3d0*parj(4)/(1d0+3d0*parj(4)))+1
669  IF(kfjh(1).EQ.kfjh(2)) kfls=3
670  kfjs(jt)=isign(1000*max(iabs(kfjh(1)),iabs(kfjh(2)))+
671  & 100*min(iabs(kfjh(1)),iabs(kfjh(2)))+kfls,kfjh(1))
672  DO 630 j=1,4
673  pjs(jt,j)=pju(1,j)+pju(2,j)+p(mju(jt),j)
674  pjs(jt+2,j)=pju(4,j)+pju(5,j)
675  630 CONTINUE
676  pjs(jt,5)=sqrt(max(0d0,pjs(jt,4)**2-pjs(jt,1)**2-pjs(jt,2)**2-
677  & pjs(jt,3)**2))
678  pjs(jt+2,5)=0d0
679  640 CONTINUE
680 
681 C...Open versus closed strings. Choose breakup region for latter.
682  650 IF(mju(1).NE.0.AND.mju(2).NE.0) THEN
683  ns=mju(2)-mju(1)
684  nb=mju(1)-n
685  ELSEIF(mju(1).NE.0) THEN
686  ns=n+nr-mju(1)
687  nb=mju(1)-n
688  ELSEIF(mju(2).NE.0) THEN
689  ns=mju(2)-n
690  nb=1
691  ELSEIF(iabs(k(n+1,2)).NE.21) THEN
692  ns=nr-1
693  nb=1
694  ELSE
695  ns=nr+1
696  w2sum=0d0
697  DO 660 is=1,nr
698  p(n+nr+is,1)=0.5d0*four(n+is,n+is+1-nr*(is/nr))
699  w2sum=w2sum+p(n+nr+is,1)
700  660 CONTINUE
701  w2ran=pyr(0)*w2sum
702  nb=0
703  670 nb=nb+1
704  w2sum=w2sum-p(n+nr+nb,1)
705  IF(w2sum.GT.w2ran.AND.nb.LT.nr) goto 670
706  ENDIF
707 
708 C...Find longitudinal string directions (i.e. lightlike four-vectors).
709  DO 700 is=1,ns
710  is1=n+is+nb-1-nr*((is+nb-2)/nr)
711  is2=n+is+nb-nr*((is+nb-1)/nr)
712  DO 680 j=1,5
713  dp(1,j)=p(is1,j)
714  IF(iabs(k(is1,2)).EQ.21) dp(1,j)=0.5d0*dp(1,j)
715  IF(is1.EQ.mju(1)) dp(1,j)=pjs(1,j)-pjs(3,j)
716  dp(2,j)=p(is2,j)
717  IF(iabs(k(is2,2)).EQ.21) dp(2,j)=0.5d0*dp(2,j)
718  IF(is2.EQ.mju(2)) dp(2,j)=pjs(2,j)-pjs(4,j)
719  680 CONTINUE
720  IF(is1.EQ.mju(1)) dp(1,5)=sqrt(max(0d0,dp(1,4)**2-dp(1,1)**2-
721  & dp(1,2)**2-dp(1,3)**2))
722  IF(is2.EQ.mju(2)) dp(2,5)=sqrt(max(0d0,dp(2,4)**2-dp(2,1)**2-
723  & dp(2,2)**2-dp(2,3)**2))
724  dp(3,5)=dfour(1,1)
725  dp(4,5)=dfour(2,2)
726  dhkc=dfour(1,2)
727  IF(dp(3,5)+2d0*dhkc+dp(4,5).LE.0d0) goto 200
728  dhks=sqrt(dhkc**2-dp(3,5)*dp(4,5))
729  dhk1=0.5d0*((dp(4,5)+dhkc)/dhks-1d0)
730  dhk2=0.5d0*((dp(3,5)+dhkc)/dhks-1d0)
731  in1=n+nr+4*is-3
732  p(in1,5)=sqrt(dp(3,5)+2d0*dhkc+dp(4,5))
733  DO 690 j=1,4
734  p(in1,j)=(1d0+dhk1)*dp(1,j)-dhk2*dp(2,j)
735  p(in1+1,j)=(1d0+dhk2)*dp(2,j)-dhk1*dp(1,j)
736  690 CONTINUE
737  700 CONTINUE
738 
739 C...Begin initialization: sum up energy, set starting position.
740  isav=i
741  mstu91=mstu(90)
742  710 ntry=ntry+1
743  IF(ntry.GT.100.AND.ntryr.LE.8.AND.nr.GT.nrmin) THEN
744  paru12=4d0*paru12
745  paru13=2d0*paru13
746  goto 140
747  ELSEIF(ntry.GT.100) THEN
748  CALL pyerrm(14,'(PYSTRF:) caught in infinite loop')
749  IF(mstu(21).GE.1) RETURN
750  ENDIF
751  i=isav
752  mstu(90)=mstu91
753  DO 730 j=1,4
754  p(n+nrs,j)=0d0
755  DO 720 is=1,nr
756  p(n+nrs,j)=p(n+nrs,j)+p(n+is,j)
757  720 CONTINUE
758  730 CONTINUE
759  DO 750 jt=1,2
760  irank(jt)=0
761  IF(mju(jt).NE.0) irank(jt)=njs(jt)
762  IF(ns.GT.nr) irank(jt)=1
763  ibarrk(jt)=0
764  ie(jt)=k(n+1+(jt/2)*(np-1),3)
765  in(3*jt+1)=n+nr+1+4*(jt/2)*(ns-1)
766  in(3*jt+2)=in(3*jt+1)+1
767  in(3*jt+3)=n+nr+4*ns+2*jt-1
768  DO 740 in1=n+nr+2+jt,n+nr+4*ns-2+jt,4
769  p(in1,1)=2-jt
770  p(in1,2)=jt-1
771  p(in1,3)=1d0
772  740 CONTINUE
773  750 CONTINUE
774 
775 C.. MOPS variables and switches
776  nrvmo=0
777  xbmo=1d0
778  mstu(121)=0
779  mstu(122)=0
780 
781 C...Initialize flavour and pT variables for open string.
782  IF(ns.LT.nr) THEN
783  px(1)=0d0
784  py(1)=0d0
785  IF(ns.EQ.1.AND.mju(1)+mju(2).EQ.0) CALL pyptdi(0,px(1),py(1))
786  px(2)=-px(1)
787  py(2)=-py(1)
788  DO 760 jt=1,2
789  kfl(jt)=k(ie(jt),2)
790  IF(mju(jt).NE.0) kfl(jt)=kfjs(jt)
791  IF(mju(jt).NE.0.AND.iabs(kfl(jt)).GT.1000) ibarrk(jt)=1
792  mstj(93)=1
793  pmq(jt)=pymass(kfl(jt))
794  gam(jt)=0d0
795  760 CONTINUE
796 
797 C...Closed string: random initial breakup flavour, pT and vertex.
798  ELSE
799  kfl(3)=int(1d0+(2d0+parj(2))*pyr(0))*(-1)**int(pyr(0)+0.5d0)
800  ibmo=0
801  770 CALL pykfdi(kfl(3),0,kfl(1),kdump)
802 C.. Closed string: first vertex diq attempt => enforced second
803 C.. vertex diq
804  IF(iabs(kfl(1)).GT.10)THEN
805  ibmo=1
806  mstu(121)=0
807  goto 770
808  ENDIF
809  IF(ibmo.EQ.1) mstu(121)=-1
810  kfl(2)=-kfl(1)
811  CALL pyptdi(kfl(1),px(1),py(1))
812  px(2)=-px(1)
813  py(2)=-py(1)
814  pr3=min(25d0,0.1d0*p(n+nr+1,5)**2)
815  780 CALL pyzdis(kfl(1),kfl(2),pr3,z)
816  zr=pr3/(z*p(n+nr+1,5)**2)
817  IF(zr.GE.1d0) goto 780
818  DO 790 jt=1,2
819  mstj(93)=1
820  pmq(jt)=pymass(kfl(jt))
821  gam(jt)=pr3*(1d0-z)/z
822  in1=n+nr+3+4*(jt/2)*(ns-1)
823  p(in1,jt)=1d0-z
824  p(in1,3-jt)=jt-1
825  p(in1,3)=(2-jt)*(1d0-z)+(jt-1)*z
826  p(in1+1,jt)=zr
827  p(in1+1,3-jt)=2-jt
828  p(in1+1,3)=(2-jt)*(1d0-zr)+(jt-1)*zr
829  790 CONTINUE
830  ENDIF
831 C.. MOPS variables
832  DO 800 jt=1,2
833  xtmo(jt)=1d0
834  pm2qmo(jt)=pmq(jt)**2
835  IF(iabs(kfl(jt)).GT.10) pm2qmo(jt)=0d0
836  800 CONTINUE
837 
838 C...Find initial transverse directions (i.e. spacelike four-vectors).
839  DO 840 jt=1,2
840  IF(jt.EQ.1.OR.ns.EQ.nr-1.OR.mju(1)+mju(2).NE.0) THEN
841  in1=in(3*jt+1)
842  in3=in(3*jt+3)
843  DO 810 j=1,4
844  dp(1,j)=p(in1,j)
845  dp(2,j)=p(in1+1,j)
846  dp(3,j)=0d0
847  dp(4,j)=0d0
848  810 CONTINUE
849  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
850  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
851  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
852  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
853  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
854  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1d0
855  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1d0
856  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1d0
857  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1d0
858  dhc12=dfour(1,2)
859  dhcx1=dfour(3,1)/dhc12
860  dhcx2=dfour(3,2)/dhc12
861  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
862  dhcy1=dfour(4,1)/dhc12
863  dhcy2=dfour(4,2)/dhc12
864  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
865  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
866  DO 820 j=1,4
867  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
868  p(in3,j)=dp(3,j)
869  p(in3+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
870  & dhcyx*dp(3,j))
871  820 CONTINUE
872  ELSE
873  DO 830 j=1,4
874  p(in3+2,j)=p(in3,j)
875  p(in3+3,j)=p(in3+1,j)
876  830 CONTINUE
877  ENDIF
878  840 CONTINUE
879 
880 C...Remove energy used up in junction string fragmentation.
881  IF(mju(1)+mju(2).GT.0) THEN
882  DO 860 jt=1,2
883  IF(njs(jt).EQ.0) goto 860
884  DO 850 j=1,4
885  p(n+nrs,j)=p(n+nrs,j)-pjs(jt+2,j)
886  850 CONTINUE
887  860 CONTINUE
888  parjst=parj(33)
889  IF(mstj(11).EQ.2) parjst=parj(34)
890  wmin=parjst+pmq(1)+pmq(2)
891  wrem2=four(n+nrs,n+nrs)
892  IF(p(n+nrs,4).LT.0d0.OR.wrem2.LT.wmin**2) THEN
893  ntrywr=ntrywr+1
894  IF(mod(ntrywr,20).NE.0) ntryr=ntryr-1
895  goto 140
896  ENDIF
897  ENDIF
898 
899 C...Produce new particle: side, origin.
900  870 i=i+1
901  IF(2*i-nsav.GE.mstu(4)-mstu(32)-5) THEN
902  CALL pyerrm(11,'(PYSTRF:) no more memory left in PYJETS')
903  IF(mstu(21).GE.1) RETURN
904  ENDIF
905 C.. New side priority for popcorn systems
906  IF(mstu(121).LE.0)THEN
907  jt=1.5d0+pyr(0)
908  IF(iabs(kfl(3-jt)).GT.10) jt=3-jt
909  IF(iabs(kfl(3-jt)).GE.4.AND.iabs(kfl(3-jt)).LE.8) jt=3-jt
910  ENDIF
911  jr=3-jt
912  js=3-2*jt
913  irank(jt)=irank(jt)+1
914  k(i,1)=1
915  k(i,4)=0
916  k(i,5)=0
917 
918 C...Generate flavour, hadron and pT.
919  880 k(i,3)=ie(jt)
920  CALL pykfdi(kfl(jt),0,kfl(3),k(i,2))
921  IF(k(i,2).EQ.0) goto 710
922  mu90mo=mstu(90)
923  IF(mstu(121).EQ.-1) goto 910
924  IF(irank(jt).EQ.1.AND.iabs(kfl(jt)).LE.10.AND.
925  &iabs(kfl(3)).GT.10) THEN
926  IF(pyr(0).GT.parj(19)) goto 880
927  ENDIF
928  IF(ibarrk(jt).EQ.1.AND.mod(iabs(k(i,2)),10000).GT.1000)
929  &k(i,3)=ijuori(jt)
930  p(i,5)=pymass(k(i,2))
931  CALL pyptdi(kfl(jt),px(3),py(3))
932  pr(jt)=p(i,5)**2+(px(jt)+px(3))**2+(py(jt)+py(3))**2
933 
934 C...Final hadrons for small invariant mass.
935  mstj(93)=1
936  pmq(3)=pymass(kfl(3))
937  parjst=parj(33)
938  IF(mstj(11).EQ.2) parjst=parj(34)
939  wmin=parjst+pmq(1)+pmq(2)+parj(36)*pmq(3)
940  IF(iabs(kfl(jt)).GT.10.AND.iabs(kfl(3)).GT.10) wmin=
941  &wmin-0.5d0*parj(36)*pmq(3)
942  wrem2=four(n+nrs,n+nrs)
943  IF(wrem2.LT.0.10d0) goto 710
944  IF(wrem2.LT.max(wmin*(1d0+(2d0*pyr(0)-1d0)*parj(37)),
945  &parj(32)+pmq(1)+pmq(2))**2) goto 1080
946 
947 C...Choose z, which gives Gamma. Shift z for heavy flavours.
948  CALL pyzdis(kfl(jt),kfl(3),pr(jt),z)
949  IF(iabs(kfl(jt)).GE.4.AND.iabs(kfl(jt)).LE.8.AND.
950  &mstu(90).LT.8) THEN
951  mstu(90)=mstu(90)+1
952  mstu(90+mstu(90))=i
953  paru(90+mstu(90))=z
954  ENDIF
955  kfl1a=iabs(kfl(1))
956  kfl2a=iabs(kfl(2))
957  IF(max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
958  &mod(kfl2a/1000,10)).GE.4) THEN
959  pr(jr)=(pmq(jr)+pmq(3))**2+(px(jr)-px(3))**2+(py(jr)-py(3))**2
960  pw12=sqrt(max(0d0,(wrem2-pr(1)-pr(2))**2-4d0*pr(1)*pr(2)))
961  z=(wrem2+pr(jt)-pr(jr)+pw12*(2d0*z-1d0))/(2d0*wrem2)
962  pr(jr)=(pmq(jr)+parjst)**2+(px(jr)-px(3))**2+(py(jr)-py(3))**2
963  IF((1d0-z)*(wrem2-pr(jt)/z).LT.pr(jr)) goto 1080
964  ENDIF
965  gam(3)=(1d0-z)*(gam(jt)+pr(jt)/z)
966 
967 C.. MOPS baryon model modification
968  xtmo3=(1d0-z)*xtmo(jt)
969  IF(iabs(kfl(3)).LE.10) nrvmo=0
970  IF(iabs(kfl(3)).GT.10.AND.mstj(12).GE.4) THEN
971  gtstmo=1d0
972  ptstmo=1d0
973  rtstmo=pyr(0)
974  IF(iabs(kfl(jt)).LE.10)THEN
975  xbmo=min(xtmo3,1d0-(2d-10))
976  gbmo=gam(3)
977  pmmo=0d0
978  pgmo=gbmo+log(1d0-xbmo)*pm2qmo(jt)
979  gtstmo=1d0-parf(192)**pgmo
980  ELSE
981  IF(irank(jt).EQ.1) THEN
982  gbmo=gam(jt)
983  pmmo=0d0
984  xbmo=1d0
985  ENDIF
986  IF(xbmo.LT.1d0-(1d-10))THEN
987  pgnmo=gbmo*xtmo3/xbmo+pm2qmo(jt)*log(1d0-xtmo3)
988  gtstmo=(1d0-parf(192)**pgnmo)/(1d0-parf(192)**pgmo)
989  pgmo=pgnmo
990  ENDIF
991  IF(mstj(12).GE.5)THEN
992  pmnmo=sqrt((xbmo-xtmo3)*(gam(3)/xtmo3-gbmo/xbmo))
993  pmmo=pmmo+pmas(pycomp(k(i,2)),1)-pmas(pycomp(k(i,2)),3)
994  ptstmo=exp((pmmo-pmnmo)*parf(193))
995  pmmo=pmnmo
996  ENDIF
997  ENDIF
998 
999 C.. MOPS Accepting popcorn system hadron.
1000  IF(ptstmo*gtstmo.GT.rtstmo) THEN
1001  IF(irank(jt).EQ.1.OR.iabs(kfl(jt)).LE.10) THEN
1002  nrvmo=i-n-nr
1003  IF(i+nrvmo.GT.mstu(4)-mstu(32)-5) THEN
1004  CALL pyerrm(11,
1005  & '(PYSTRF:) no more memory left in PYJETS')
1006  IF(mstu(21).GE.1) RETURN
1007  ENDIF
1008  imo=i
1009  kflmo=kfl(jt)
1010  pmqmo=pmq(jt)
1011  pxmo=px(jt)
1012  pymo=py(jt)
1013  gammo=gam(jt)
1014  irmo=irank(jt)
1015  xmo=xtmo(jt)
1016  DO 900 j=1,9
1017  IF(j.LE.5) THEN
1018  DO 890 line=1,i-n-nr
1019  p(mstu(4)-mstu(32)-line,j)=p(n+nr+line,j)
1020  k(mstu(4)-mstu(32)-line,j)=k(n+nr+line,j)
1021  890 CONTINUE
1022  ENDIF
1023  inmo(j)=in(j)
1024  900 CONTINUE
1025  ENDIF
1026  ELSE
1027 C..Reject popcorn system, flag=-1 if enforcing new one
1028  mstu(121)=-1
1029  IF(ptstmo.GT.rtstmo) mstu(121)=-2
1030  ENDIF
1031  ENDIF
1032 
1033 
1034 C..Lift restoring string outside MOPS block
1035  910 IF(mstu(121).LT.0) THEN
1036  IF(mstu(121).EQ.-2) mstu(121)=0
1037  mstu(90)=mu90mo
1038  nrvmo=0
1039  IF(irank(jt).EQ.1.OR.iabs(kfl(jt)).LE.10) goto 880
1040  i=imo
1041  kfl(jt)=kflmo
1042  pmq(jt)=pmqmo
1043  px(jt)=pxmo
1044  py(jt)=pymo
1045  gam(jt)=gammo
1046  irank(jt)=irmo
1047  xtmo(jt)=xmo
1048  DO 930 j=1,9
1049  IF(j.LE.5) THEN
1050  DO 920 line=1,i-n-nr
1051  p(n+nr+line,j)=p(mstu(4)-mstu(32)-line,j)
1052  k(n+nr+line,j)=k(mstu(4)-mstu(32)-line,j)
1053  920 CONTINUE
1054  ENDIF
1055  in(j)=inmo(j)
1056  930 CONTINUE
1057  goto 880
1058  ENDIF
1059  xtmo(jt)=xtmo3
1060 C.. MOPS end of modification
1061 
1062  DO 940 j=1,3
1063  in(j)=in(3*jt+j)
1064  940 CONTINUE
1065 
1066 C...Stepping within or from 'low' string region easy.
1067  IF(in(1)+1.EQ.in(2).AND.z*p(in(1)+2,3)*p(in(2)+2,3)*
1068  &p(in(1),5)**2.GE.pr(jt)) THEN
1069  p(in(jt)+2,4)=z*p(in(jt)+2,3)
1070  p(in(jr)+2,4)=pr(jt)/(p(in(jt)+2,4)*p(in(1),5)**2)
1071  DO 950 j=1,4
1072  p(i,j)=(px(jt)+px(3))*p(in(3),j)+(py(jt)+py(3))*p(in(3)+1,j)
1073  950 CONTINUE
1074  goto 1040
1075  ELSEIF(in(1)+1.EQ.in(2)) THEN
1076  p(in(jr)+2,4)=p(in(jr)+2,3)
1077  p(in(jr)+2,jt)=1d0
1078  in(jr)=in(jr)+4*js
1079  IF(js*in(jr).GT.js*in(4*jr)) goto 710
1080  IF(four(in(1),in(2)).LE.1d-2) THEN
1081  p(in(jt)+2,4)=p(in(jt)+2,3)
1082  p(in(jt)+2,jt)=0d0
1083  in(jt)=in(jt)+4*js
1084  ENDIF
1085  ENDIF
1086 
1087 C...Find new transverse directions (i.e. spacelike string vectors).
1088  960 IF(js*in(1).GT.js*in(3*jr+1).OR.js*in(2).GT.js*in(3*jr+2).OR.
1089  &in(1).GT.in(2)) goto 710
1090  IF(in(1).NE.in(3*jt+1).OR.in(2).NE.in(3*jt+2)) THEN
1091  DO 970 j=1,4
1092  dp(1,j)=p(in(1),j)
1093  dp(2,j)=p(in(2),j)
1094  dp(3,j)=0d0
1095  dp(4,j)=0d0
1096  970 CONTINUE
1097  dp(1,4)=sqrt(dp(1,1)**2+dp(1,2)**2+dp(1,3)**2)
1098  dp(2,4)=sqrt(dp(2,1)**2+dp(2,2)**2+dp(2,3)**2)
1099  dhc12=dfour(1,2)
1100  IF(dhc12.LE.1d-2) THEN
1101  p(in(jt)+2,4)=p(in(jt)+2,3)
1102  p(in(jt)+2,jt)=0d0
1103  in(jt)=in(jt)+4*js
1104  goto 960
1105  ENDIF
1106  in(3)=n+nr+4*ns+5
1107  dp(5,1)=dp(1,1)/dp(1,4)-dp(2,1)/dp(2,4)
1108  dp(5,2)=dp(1,2)/dp(1,4)-dp(2,2)/dp(2,4)
1109  dp(5,3)=dp(1,3)/dp(1,4)-dp(2,3)/dp(2,4)
1110  IF(dp(5,1)**2.LE.dp(5,2)**2+dp(5,3)**2) dp(3,1)=1d0
1111  IF(dp(5,1)**2.GT.dp(5,2)**2+dp(5,3)**2) dp(3,3)=1d0
1112  IF(dp(5,2)**2.LE.dp(5,1)**2+dp(5,3)**2) dp(4,2)=1d0
1113  IF(dp(5,2)**2.GT.dp(5,1)**2+dp(5,3)**2) dp(4,3)=1d0
1114  dhcx1=dfour(3,1)/dhc12
1115  dhcx2=dfour(3,2)/dhc12
1116  dhcxx=1d0/sqrt(1d0+2d0*dhcx1*dhcx2*dhc12)
1117  dhcy1=dfour(4,1)/dhc12
1118  dhcy2=dfour(4,2)/dhc12
1119  dhcyx=dhcxx*(dhcx1*dhcy2+dhcx2*dhcy1)*dhc12
1120  dhcyy=1d0/sqrt(1d0+2d0*dhcy1*dhcy2*dhc12-dhcyx**2)
1121  DO 980 j=1,4
1122  dp(3,j)=dhcxx*(dp(3,j)-dhcx2*dp(1,j)-dhcx1*dp(2,j))
1123  p(in(3),j)=dp(3,j)
1124  p(in(3)+1,j)=dhcyy*(dp(4,j)-dhcy2*dp(1,j)-dhcy1*dp(2,j)-
1125  & dhcyx*dp(3,j))
1126  980 CONTINUE
1127 C...Express pT with respect to new axes, if sensible.
1128  pxp=-(px(3)*four(in(3*jt+3),in(3))+py(3)*
1129  & four(in(3*jt+3)+1,in(3)))
1130  pyp=-(px(3)*four(in(3*jt+3),in(3)+1)+py(3)*
1131  & four(in(3*jt+3)+1,in(3)+1))
1132  IF(abs(pxp**2+pyp**2-px(3)**2-py(3)**2).LT.0.01d0) THEN
1133  px(3)=pxp
1134  py(3)=pyp
1135  ENDIF
1136  ENDIF
1137 
1138 C...Sum up known four-momentum. Gives coefficients for m2 expression.
1139  DO 1010 j=1,4
1140  dhg(j)=0d0
1141  p(i,j)=px(jt)*p(in(3*jt+3),j)+py(jt)*p(in(3*jt+3)+1,j)+
1142  & px(3)*p(in(3),j)+py(3)*p(in(3)+1,j)
1143  DO 990 in1=in(3*jt+1),in(1)-4*js,4*js
1144  p(i,j)=p(i,j)+p(in1+2,3)*p(in1,j)
1145  990 CONTINUE
1146  DO 1000 in2=in(3*jt+2),in(2)-4*js,4*js
1147  p(i,j)=p(i,j)+p(in2+2,3)*p(in2,j)
1148  1000 CONTINUE
1149  1010 CONTINUE
1150  dhm(1)=four(i,i)
1151  dhm(2)=2d0*four(i,in(1))
1152  dhm(3)=2d0*four(i,in(2))
1153  dhm(4)=2d0*four(in(1),in(2))
1154 
1155 C...Find coefficients for Gamma expression.
1156  DO 1030 in2=in(1)+1,in(2),4
1157  DO 1020 in1=in(1),in2-1,4
1158  dhc=2d0*four(in1,in2)
1159  dhg(1)=dhg(1)+p(in1+2,jt)*p(in2+2,jt)*dhc
1160  IF(in1.EQ.in(1)) dhg(2)=dhg(2)-js*p(in2+2,jt)*dhc
1161  IF(in2.EQ.in(2)) dhg(3)=dhg(3)+js*p(in1+2,jt)*dhc
1162  IF(in1.EQ.in(1).AND.in2.EQ.in(2)) dhg(4)=dhg(4)-dhc
1163  1020 CONTINUE
1164  1030 CONTINUE
1165 
1166 C...Solve (m2, Gamma) equation system for energies taken.
1167  dhs1=dhm(jr+1)*dhg(4)-dhm(4)*dhg(jr+1)
1168  IF(abs(dhs1).LT.1d-4) goto 710
1169  dhs2=dhm(4)*(gam(3)-dhg(1))-dhm(jt+1)*dhg(jr+1)-dhg(4)*
1170  &(p(i,5)**2-dhm(1))+dhg(jt+1)*dhm(jr+1)
1171  dhs3=dhm(jt+1)*(gam(3)-dhg(1))-dhg(jt+1)*(p(i,5)**2-dhm(1))
1172  p(in(jr)+2,4)=0.5d0*(sqrt(max(0d0,dhs2**2-4d0*dhs1*dhs3))/
1173  &abs(dhs1)-dhs2/dhs1)
1174  IF(dhm(jt+1)+dhm(4)*p(in(jr)+2,4).LE.0d0) goto 710
1175  p(in(jt)+2,4)=(p(i,5)**2-dhm(1)-dhm(jr+1)*p(in(jr)+2,4))/
1176  &(dhm(jt+1)+dhm(4)*p(in(jr)+2,4))
1177 
1178 C...Step to new region if necessary.
1179  IF(p(in(jr)+2,4).GT.p(in(jr)+2,3)) THEN
1180  p(in(jr)+2,4)=p(in(jr)+2,3)
1181  p(in(jr)+2,jt)=1d0
1182  in(jr)=in(jr)+4*js
1183  IF(js*in(jr).GT.js*in(4*jr)) goto 710
1184  IF(four(in(1),in(2)).LE.1d-2) THEN
1185  p(in(jt)+2,4)=p(in(jt)+2,3)
1186  p(in(jt)+2,jt)=0d0
1187  in(jt)=in(jt)+4*js
1188  ENDIF
1189  goto 960
1190  ELSEIF(p(in(jt)+2,4).GT.p(in(jt)+2,3)) THEN
1191  p(in(jt)+2,4)=p(in(jt)+2,3)
1192  p(in(jt)+2,jt)=0d0
1193  in(jt)=in(jt)+4*js
1194  goto 960
1195  ENDIF
1196 
1197 C...Four-momentum of particle. Remaining quantities. Loop back.
1198  1040 DO 1050 j=1,4
1199  p(i,j)=p(i,j)+p(in(1)+2,4)*p(in(1),j)+p(in(2)+2,4)*p(in(2),j)
1200  p(n+nrs,j)=p(n+nrs,j)-p(i,j)
1201  1050 CONTINUE
1202  IF(p(in(1)+2,4).GT.1d0+paru(14).OR.p(in(1)+2,4).LT.-paru(14).OR.
1203  &p(in(2)+2,4).GT.1d0+paru(14).OR.p(in(2)+2,4).LT.-paru(14))
1204  &goto 200
1205  IF(p(i,4).LT.p(i,5)) goto 710
1206  kfl(jt)=-kfl(3)
1207  pmq(jt)=pmq(3)
1208  px(jt)=-px(3)
1209  py(jt)=-py(3)
1210  gam(jt)=gam(3)
1211  IF(in(3).NE.in(3*jt+3)) THEN
1212  DO 1060 j=1,4
1213  p(in(3*jt+3),j)=p(in(3),j)
1214  p(in(3*jt+3)+1,j)=p(in(3)+1,j)
1215  1060 CONTINUE
1216  ENDIF
1217  DO 1070 jq=1,2
1218  in(3*jt+jq)=in(jq)
1219  p(in(jq)+2,3)=p(in(jq)+2,3)-p(in(jq)+2,4)
1220  p(in(jq)+2,jt)=p(in(jq)+2,jt)-js*(3-2*jq)*p(in(jq)+2,4)
1221  1070 CONTINUE
1222  IF(ibarrk(jt).EQ.1.AND.mod(iabs(k(i,2)),10000).GT.1000)
1223  &ibarrk(jt)=0
1224  goto 870
1225 
1226 C...Final hadron: side, flavour, hadron, mass.
1227  1080 i=i+1
1228  k(i,1)=1
1229  k(i,3)=ie(jr)
1230  k(i,4)=0
1231  k(i,5)=0
1232  CALL pykfdi(kfl(jr),-kfl(3),kfldmp,k(i,2))
1233  IF(k(i,2).EQ.0) goto 710
1234  IF(ibarrk(jt).EQ.1.AND.mod(iabs(k(i-1,2)),10000).GT.1000)
1235  &ibarrk(jt)=0
1236  IF(ibarrk(jt).EQ.1.AND.mod(iabs(k(i,2)),10000).GT.1000)
1237  &k(i,3)=ijuori(jt)
1238  IF(ibarrk(jr).EQ.1.AND.mod(iabs(k(i,2)),10000).GT.1000)
1239  &k(i,3)=ijuori(jr)
1240  p(i,5)=pymass(k(i,2))
1241  pr(jr)=p(i,5)**2+(px(jr)-px(3))**2+(py(jr)-py(3))**2
1242 
1243 C...Final two hadrons: find common setup of four-vectors.
1244  jq=1
1245  IF(p(in(4)+2,3)*p(in(5)+2,3)*four(in(4),in(5)).LT.
1246  &p(in(7)+2,3)*p(in(8)+2,3)*four(in(7),in(8))) jq=2
1247  dhc12=four(in(3*jq+1),in(3*jq+2))
1248  dhr1=four(n+nrs,in(3*jq+2))/dhc12
1249  dhr2=four(n+nrs,in(3*jq+1))/dhc12
1250  IF(in(4).NE.in(7).OR.in(5).NE.in(8)) THEN
1251  px(3-jq)=-four(n+nrs,in(3*jq+3))-px(jq)
1252  py(3-jq)=-four(n+nrs,in(3*jq+3)+1)-py(jq)
1253  pr(3-jq)=p(i+(jt+jq-3)**2-1,5)**2+(px(3-jq)+(2*jq-3)*js*
1254  & px(3))**2+(py(3-jq)+(2*jq-3)*js*py(3))**2
1255  ENDIF
1256 
1257 C...Solve kinematics for final two hadrons, if possible.
1258  wrem2=2d0*dhr1*dhr2*dhc12
1259  fd=(sqrt(pr(1))+sqrt(pr(2)))/sqrt(wrem2)
1260  IF(mju(1)+mju(2).NE.0.AND.i.EQ.isav+2.AND.fd.GE.1d0) goto 200
1261  IF(fd.GE.1d0) goto 710
1262  fa=wrem2+pr(jt)-pr(jr)
1263  fb=sqrt(max(0d0,fa**2-4d0*wrem2*pr(jt)))
1264  prevcf=parj(42)
1265  IF(mstj(11).EQ.2) prevcf=parj(39)
1266  prev=1d0/(1d0+exp(min(50d0,prevcf*fb*parj(40))))
1267  fb=sign(fb,js*(pyr(0)-prev))
1268  kfl1a=iabs(kfl(1))
1269  kfl2a=iabs(kfl(2))
1270  IF(max(mod(kfl1a,10),mod(kfl1a/1000,10),mod(kfl2a,10),
1271  &mod(kfl2a/1000,10)).GE.6) fb=sign(sqrt(max(0d0,fa**2-
1272  &4d0*wrem2*pr(jt))),dble(js))
1273  DO 1090 j=1,4
1274  p(i-1,j)=(px(jt)+px(3))*p(in(3*jq+3),j)+(py(jt)+py(3))*
1275  & p(in(3*jq+3)+1,j)+0.5d0*(dhr1*(fa+fb)*p(in(3*jq+1),j)+
1276  & dhr2*(fa-fb)*p(in(3*jq+2),j))/wrem2
1277  p(i,j)=p(n+nrs,j)-p(i-1,j)
1278  1090 CONTINUE
1279  IF(p(i-1,4).LT.p(i-1,5).OR.p(i,4).LT.p(i,5)) goto 710
1280  dm2f1=p(i-1,4)**2-p(i-1,1)**2-p(i-1,2)**2-p(i-1,3)**2-p(i-1,5)**2
1281  dm2f2=p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2-p(i,5)**2
1282  IF(dm2f1.GT.1d-10*p(i-1,4)**2.OR.dm2f2.GT.1d-10*p(i,4)**2) THEN
1283  ntryfn=ntryfn+1
1284  IF(ntryfn.LT.100) goto 140
1285  CALL pyerrm(13,'(PYSTRF:) bad energies for final two hadrons')
1286  ENDIF
1287 
1288 C...Mark jets as fragmented and give daughter pointers.
1289  n=i-nrs+1
1290  DO 1100 i=nsav+1,nsav+np
1291  im=k(i,3)
1292  k(im,1)=k(im,1)+10
1293  IF(mstu(16).NE.2) THEN
1294  k(im,4)=nsav+1
1295  k(im,5)=nsav+1
1296  ELSE
1297  k(im,4)=nsav+2
1298  k(im,5)=n
1299  ENDIF
1300  1100 CONTINUE
1301 
1302 C...Document string system. Move up particles.
1303  nsav=nsav+1
1304  k(nsav,1)=11
1305  k(nsav,2)=92
1306  k(nsav,3)=ip
1307  k(nsav,4)=nsav+1
1308  k(nsav,5)=n
1309  DO 1110 j=1,4
1310  p(nsav,j)=dps(j)
1311  v(nsav,j)=v(ip,j)
1312  1110 CONTINUE
1313  p(nsav,5)=sqrt(max(0d0,dps(4)**2-dps(1)**2-dps(2)**2-dps(3)**2))
1314  v(nsav,5)=0d0
1315  DO 1130 i=nsav+1,n
1316  DO 1120 j=1,5
1317  k(i,j)=k(i+nrs-1,j)
1318  p(i,j)=p(i+nrs-1,j)
1319  v(i,j)=0d0
1320  1120 CONTINUE
1321  1130 CONTINUE
1322  mstu91=mstu(90)
1323  DO 1140 iz=mstu90+1,mstu91
1324  mstu9t(iz)=mstu(90+iz)-nrs+1-nsav+n
1325  paru9t(iz)=paru(90+iz)
1326  1140 CONTINUE
1327  mstu(90)=mstu90
1328 
1329 C...Order particles in rank along the chain. Update mother pointer.
1330  DO 1160 i=nsav+1,n
1331  DO 1150 j=1,5
1332  k(i-nsav+n,j)=k(i,j)
1333  p(i-nsav+n,j)=p(i,j)
1334  1150 CONTINUE
1335  1160 CONTINUE
1336  i1=nsav
1337  DO 1190 i=n+1,2*n-nsav
1338  IF(k(i,3).NE.ie(1).AND.k(i,3).NE.ijuori(1)) goto 1190
1339  i1=i1+1
1340  DO 1170 j=1,5
1341  k(i1,j)=k(i,j)
1342  p(i1,j)=p(i,j)
1343  1170 CONTINUE
1344  IF(mstu(16).NE.2) k(i1,3)=nsav
1345  DO 1180 iz=mstu90+1,mstu91
1346  IF(mstu9t(iz).EQ.i) THEN
1347  mstu(90)=mstu(90)+1
1348  mstu(90+mstu(90))=i1
1349  paru(90+mstu(90))=paru9t(iz)
1350  ENDIF
1351  1180 CONTINUE
1352  1190 CONTINUE
1353  DO 1220 i=2*n-nsav,n+1,-1
1354  IF(k(i,3).EQ.ie(1).OR.k(i,3).EQ.ijuori(1)) goto 1220
1355  i1=i1+1
1356  DO 1200 j=1,5
1357  k(i1,j)=k(i,j)
1358  p(i1,j)=p(i,j)
1359  1200 CONTINUE
1360  IF(mstu(16).NE.2) k(i1,3)=nsav
1361  DO 1210 iz=mstu90+1,mstu91
1362  IF(mstu9t(iz).EQ.i) THEN
1363  mstu(90)=mstu(90)+1
1364  mstu(90+mstu(90))=i1
1365  paru(90+mstu(90))=paru9t(iz)
1366  ENDIF
1367  1210 CONTINUE
1368  1220 CONTINUE
1369 
1370 C...Boost back particle system. Set production vertices.
1371  IF(mbst.EQ.0) THEN
1372  mstu(33)=1
1373  CALL pyrobo(nsav+1,n,0d0,0d0,dps(1)/dps(4),dps(2)/dps(4),
1374  & dps(3)/dps(4))
1375  ELSE
1376  DO 1230 i=nsav+1,n
1377  hhpmt=p(i,1)**2+p(i,2)**2+p(i,5)**2
1378  IF(p(i,3).GT.0d0) THEN
1379  hhpez=(p(i,4)+p(i,3))*hhbz
1380  p(i,3)=0.5d0*(hhpez-hhpmt/hhpez)
1381  p(i,4)=0.5d0*(hhpez+hhpmt/hhpez)
1382  ELSE
1383  hhpez=(p(i,4)-p(i,3))/hhbz
1384  p(i,3)=-0.5d0*(hhpez-hhpmt/hhpez)
1385  p(i,4)=0.5d0*(hhpez+hhpmt/hhpez)
1386  ENDIF
1387  1230 CONTINUE
1388  ENDIF
1389  DO 1250 i=nsav+1,n
1390  DO 1240 j=1,4
1391  v(i,j)=v(ip,j)
1392  1240 CONTINUE
1393  1250 CONTINUE
1394 
1395  RETURN
1396  END