New functions and constructors added and some other fixes.
[u/mrichter/AliRoot.git] / PYTHIA / jetset / 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(4000,5),P(4000,5),V(4000,5) 
9       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) 
10       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) 
11       SAVE /LUJETS/,/LUDAT1/,/LUDAT2/ 
12       DIMENSION DPS(5),KFL(3),PMQ(3),PX(3),PY(3),GAM(3),IE(2),PR(2), 
13      &IN(9),DHM(4),DHG(4),DP(5,5),IRANK(2),MJU(4),IJU(3),PJU(5,5), 
14      &TJU(5),KFJH(2),NJS(2),KFJS(2),PJS(4,5),MSTU9T(8),PARU9T(8) 
15  
16 C...Function: four-product of two vectors. 
17       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) 
18       DFOUR(I,J)=DP(I,4)*DP(J,4)-DP(I,1)*DP(J,1)-DP(I,2)*DP(J,2)- 
19      &DP(I,3)*DP(J,3) 
20  
21 C...Reset counters. Identify parton system. 
22       MSTJ(91)=0 
23       NSAV=N 
24       MSTU90=MSTU(90) 
25       NP=0 
26       KQSUM=0 
27       DO 100 J=1,5 
28       DPS(J)=0D0 
29   100 CONTINUE 
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       IF(J.NE.4) DPS(J)=DPS(J)+P(I,J) 
54   120 CONTINUE 
55       DPS(4)=DPS(4)+SQRT(DBLE(P(I,1))**2+DBLE(P(I,2))**2+ 
56      &DBLE(P(I,3))**2+DBLE(P(I,5))**2) 
57       K(N+NP,3)=I 
58       IF(KQ.NE.2) KQSUM=KQSUM+KQ 
59       IF(K(I,1).EQ.41) THEN 
60         KQSUM=KQSUM+2*KQ 
61         IF(KQSUM.EQ.KQ) MJU(1)=N+NP 
62         IF(KQSUM.NE.KQ) MJU(2)=N+NP 
63       ENDIF 
64       IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110 
65       IF(KQSUM.NE.0) THEN 
66         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination') 
67         IF(MSTU(21).GE.1) RETURN 
68       ENDIF 
69  
70 C...Boost copied system to CM frame (for better numerical precision). 
71       IF(ABS(DPS(3)).LT.0.99D0*DPS(4)) THEN 
72         MBST=0 
73         MSTU(33)=1 
74         CALL LUDBRB(N+1,N+NP,0.,0.,-DPS(1)/DPS(4),-DPS(2)/DPS(4), 
75      &  -DPS(3)/DPS(4)) 
76       ELSE 
77         MBST=1 
78         HHBZ=SQRT(MAX(1D-6,DPS(4)+DPS(3))/MAX(1D-6,DPS(4)-DPS(3))) 
79         DO 130 I=N+1,N+NP 
80         HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2 
81         IF(P(I,3).GT.0.) THEN 
82           HHPEZ=(P(I,4)+P(I,3))/HHBZ 
83           P(I,3)=0.5*(HHPEZ-HHPMT/HHPEZ) 
84           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ) 
85         ELSE 
86           HHPEZ=(P(I,4)-P(I,3))*HHBZ 
87           P(I,3)=-0.5*(HHPEZ-HHPMT/HHPEZ) 
88           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ) 
89         ENDIF 
90   130   CONTINUE 
91       ENDIF 
92  
93 C...Search for very nearby partons that may be recombined. 
94       NTRYR=0 
95       PARU12=PARU(12) 
96       PARU13=PARU(13) 
97       MJU(3)=MJU(1) 
98       MJU(4)=MJU(2) 
99       NR=NP 
100   140 IF(NR.GE.3) THEN 
101         PDRMIN=2.*PARU12 
102         DO 150 I=N+1,N+NR 
103         IF(I.EQ.N+NR.AND.IABS(K(N+1,2)).NE.21) GOTO 150 
104         I1=I+1 
105         IF(I.EQ.N+NR) I1=N+1 
106         IF(K(I,1).EQ.41.OR.K(I1,1).EQ.41) GOTO 150 
107         IF(MJU(1).NE.0.AND.I1.LT.MJU(1).AND.IABS(K(I1,2)).NE.21) 
108      &  GOTO 150 
109         IF(MJU(2).NE.0.AND.I.GT.MJU(2).AND.IABS(K(I,2)).NE.21) GOTO 150 
110         PAP=SQRT((P(I,1)**2+P(I,2)**2+P(I,3)**2)*(P(I1,1)**2+ 
111      &  P(I1,2)**2+P(I1,3)**2)) 
112         PVP=P(I,1)*P(I1,1)+P(I,2)*P(I1,2)+P(I,3)*P(I1,3) 
113         PDR=4.*(PAP-PVP)**2/MAX(1E-6,PARU13**2*PAP+2.*(PAP-PVP)) 
114         IF(PDR.LT.PDRMIN) THEN 
115           IR=I 
116           PDRMIN=PDR 
117         ENDIF 
118   150   CONTINUE 
119  
120 C...Recombine very nearby partons to avoid machine precision problems. 
121         IF(PDRMIN.LT.PARU12.AND.IR.EQ.N+NR) THEN 
122           DO 160 J=1,4 
123           P(N+1,J)=P(N+1,J)+P(N+NR,J) 
124   160     CONTINUE 
125           P(N+1,5)=SQRT(MAX(0.,P(N+1,4)**2-P(N+1,1)**2-P(N+1,2)**2- 
126      &    P(N+1,3)**2)) 
127           NR=NR-1 
128           GOTO 140 
129         ELSEIF(PDRMIN.LT.PARU12) THEN 
130           DO 170 J=1,4 
131           P(IR,J)=P(IR,J)+P(IR+1,J) 
132   170     CONTINUE 
133           P(IR,5)=SQRT(MAX(0.,P(IR,4)**2-P(IR,1)**2-P(IR,2)**2- 
134      &    P(IR,3)**2)) 
135           DO 190 I=IR+1,N+NR-1 
136           K(I,2)=K(I+1,2) 
137           DO 180 J=1,5 
138           P(I,J)=P(I+1,J) 
139   180     CONTINUE 
140   190     CONTINUE 
141           IF(IR.EQ.N+NR-1) K(IR,2)=K(N+NR,2) 
142           NR=NR-1 
143           IF(MJU(1).GT.IR) MJU(1)=MJU(1)-1 
144           IF(MJU(2).GT.IR) MJU(2)=MJU(2)-1 
145           GOTO 140 
146         ENDIF 
147       ENDIF 
148       NTRYR=NTRYR+1 
149  
150 C...Reset particle counter. Skip ahead if no junctions are present; 
151 C...this is usually the case! 
152       NRS=MAX(5*NR+11,NP) 
153       NTRY=0 
154   200 NTRY=NTRY+1 
155       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN 
156         PARU12=4.*PARU12 
157         PARU13=2.*PARU13 
158         GOTO 140 
159       ELSEIF(NTRY.GT.100) THEN 
160         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
161         IF(MSTU(21).GE.1) RETURN 
162       ENDIF 
163       I=N+NRS 
164       MSTU(90)=MSTU90 
165       IF(MJU(1).EQ.0.AND.MJU(2).EQ.0) GOTO 580 
166       DO 570 JT=1,2 
167       NJS(JT)=0 
168       IF(MJU(JT).EQ.0) GOTO 570 
169       JS=3-2*JT 
170  
171 C...Find and sum up momentum on three sides of junction. Check flavours. 
172       DO 220 IU=1,3 
173       IJU(IU)=0 
174       DO 210 J=1,5 
175       PJU(IU,J)=0. 
176   210 CONTINUE 
177   220 CONTINUE 
178       IU=0 
179       DO 240 I1=N+1+(JT-1)*(NR-1),N+NR+(JT-1)*(1-NR),JS 
180       IF(K(I1,2).NE.21.AND.IU.LE.2) THEN 
181         IU=IU+1 
182         IJU(IU)=I1 
183       ENDIF 
184       DO 230 J=1,4 
185       PJU(IU,J)=PJU(IU,J)+P(I1,J) 
186   230 CONTINUE 
187   240 CONTINUE 
188       DO 250 IU=1,3 
189       PJU(IU,5)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2) 
190   250 CONTINUE 
191       IF(K(IJU(3),2)/100.NE.10*K(IJU(1),2)+K(IJU(2),2).AND. 
192      &K(IJU(3),2)/100.NE.10*K(IJU(2),2)+K(IJU(1),2)) THEN 
193         CALL LUERRM(12,'(LUSTRF:) unphysical flavour combination') 
194         IF(MSTU(21).GE.1) RETURN 
195       ENDIF 
196  
197 C...Calculate (approximate) boost to rest frame of junction. 
198       T12=(PJU(1,1)*PJU(2,1)+PJU(1,2)*PJU(2,2)+PJU(1,3)*PJU(2,3))/ 
199      &(PJU(1,5)*PJU(2,5)) 
200       T13=(PJU(1,1)*PJU(3,1)+PJU(1,2)*PJU(3,2)+PJU(1,3)*PJU(3,3))/ 
201      &(PJU(1,5)*PJU(3,5)) 
202       T23=(PJU(2,1)*PJU(3,1)+PJU(2,2)*PJU(3,2)+PJU(2,3)*PJU(3,3))/ 
203      &(PJU(2,5)*PJU(3,5)) 
204       T11=SQRT((2./3.)*(1.-T12)*(1.-T13)/(1.-T23)) 
205       T22=SQRT((2./3.)*(1.-T12)*(1.-T23)/(1.-T13)) 
206       TSQ=SQRT((2.*T11*T22+T12-1.)*(1.+T12)) 
207       T1F=(TSQ-T22*(1.+T12))/(1.-T12**2) 
208       T2F=(TSQ-T11*(1.+T12))/(1.-T12**2) 
209       DO 260 J=1,3 
210       TJU(J)=-(T1F*PJU(1,J)/PJU(1,5)+T2F*PJU(2,J)/PJU(2,5)) 
211   260 CONTINUE 
212       TJU(4)=SQRT(1.+TJU(1)**2+TJU(2)**2+TJU(3)**2) 
213       DO 270 IU=1,3 
214       PJU(IU,5)=TJU(4)*PJU(IU,4)-TJU(1)*PJU(IU,1)-TJU(2)*PJU(IU,2)- 
215      &TJU(3)*PJU(IU,3) 
216   270 CONTINUE 
217  
218 C...Put junction at rest if motion could give inconsistencies. 
219       IF(PJU(1,5)+PJU(2,5).GT.PJU(1,4)+PJU(2,4)) THEN 
220         DO 280 J=1,3 
221         TJU(J)=0. 
222   280   CONTINUE 
223         TJU(4)=1. 
224         PJU(1,5)=PJU(1,4) 
225         PJU(2,5)=PJU(2,4) 
226         PJU(3,5)=PJU(3,4) 
227       ENDIF 
228  
229 C...Start preparing for fragmentation of two strings from junction. 
230       ISTA=I 
231       DO 550 IU=1,2 
232       NS=IJU(IU+1)-IJU(IU) 
233  
234 C...Junction strings: find longitudinal string directions. 
235       DO 310 IS=1,NS 
236       IS1=IJU(IU)+IS-1 
237       IS2=IJU(IU)+IS 
238       DO 290 J=1,5 
239       DP(1,J)=0.5*P(IS1,J) 
240       IF(IS.EQ.1) DP(1,J)=P(IS1,J) 
241       DP(2,J)=0.5*P(IS2,J) 
242       IF(IS.EQ.NS) DP(2,J)=-PJU(IU,J) 
243   290 CONTINUE 
244       IF(IS.EQ.NS) DP(2,4)=SQRT(PJU(IU,1)**2+PJU(IU,2)**2+PJU(IU,3)**2) 
245       IF(IS.EQ.NS) DP(2,5)=0. 
246       DP(3,5)=DFOUR(1,1) 
247       DP(4,5)=DFOUR(2,2) 
248       DHKC=DFOUR(1,2) 
249       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN 
250         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) 
251         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) 
252         DP(3,5)=0D0 
253         DP(4,5)=0D0 
254         DHKC=DFOUR(1,2) 
255       ENDIF 
256       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5)) 
257       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.) 
258       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.) 
259       IN1=N+NR+4*IS-3 
260       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5)) 
261       DO 300 J=1,4 
262       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J) 
263       P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J) 
264   300 CONTINUE 
265   310 CONTINUE 
266  
267 C...Junction strings: initialize flavour, momentum and starting pos. 
268       ISAV=I 
269       MSTU91=MSTU(90) 
270   320 NTRY=NTRY+1 
271       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN 
272         PARU12=4.*PARU12 
273         PARU13=2.*PARU13 
274         GOTO 140 
275       ELSEIF(NTRY.GT.100) THEN 
276         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
277         IF(MSTU(21).GE.1) RETURN 
278       ENDIF 
279       I=ISAV 
280       MSTU(90)=MSTU91 
281       IRANKJ=0 
282       IE(1)=K(N+1+(JT/2)*(NP-1),3) 
283       IN(4)=N+NR+1 
284       IN(5)=IN(4)+1 
285       IN(6)=N+NR+4*NS+1 
286       DO 340 JQ=1,2 
287       DO 330 IN1=N+NR+2+JQ,N+NR+4*NS-2+JQ,4 
288       P(IN1,1)=2-JQ 
289       P(IN1,2)=JQ-1 
290       P(IN1,3)=1. 
291   330 CONTINUE 
292   340 CONTINUE 
293       KFL(1)=K(IJU(IU),2) 
294       PX(1)=0. 
295       PY(1)=0. 
296       GAM(1)=0. 
297       DO 350 J=1,5 
298       PJU(IU+3,J)=0. 
299   350 CONTINUE 
300  
301 C...Junction strings: find initial transverse directions. 
302       DO 360 J=1,4 
303       DP(1,J)=P(IN(4),J) 
304       DP(2,J)=P(IN(4)+1,J) 
305       DP(3,J)=0. 
306       DP(4,J)=0. 
307   360 CONTINUE 
308       DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) 
309       DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) 
310       DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
311       DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
312       DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
313       IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1. 
314       IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1. 
315       IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1. 
316       IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1. 
317       DHC12=DFOUR(1,2) 
318       DHCX1=DFOUR(3,1)/DHC12 
319       DHCX2=DFOUR(3,2)/DHC12 
320       DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12) 
321       DHCY1=DFOUR(4,1)/DHC12 
322       DHCY2=DFOUR(4,2)/DHC12 
323       DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
324       DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2) 
325       DO 370 J=1,4 
326       DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
327       P(IN(6),J)=DP(3,J) 
328       P(IN(6)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)- 
329      &DHCYX*DP(3,J)) 
330   370 CONTINUE 
331  
332 C...Junction strings: produce new particle, origin. 
333   380 I=I+1 
334       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN 
335         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS') 
336         IF(MSTU(21).GE.1) RETURN 
337       ENDIF 
338       IRANKJ=IRANKJ+1 
339       K(I,1)=1 
340       K(I,3)=IE(1) 
341       K(I,4)=0 
342       K(I,5)=0 
343  
344 C...Junction strings: generate flavour, hadron, pT, z and Gamma. 
345   390 CALL LUKFDI(KFL(1),0,KFL(3),K(I,2)) 
346       IF(K(I,2).EQ.0) GOTO 320 
347       IF(MSTJ(12).GE.3.AND.IRANKJ.EQ.1.AND.IABS(KFL(1)).LE.10.AND. 
348      &IABS(KFL(3)).GT.10) THEN 
349         IF(RLU(0).GT.PARJ(19)) GOTO 390 
350       ENDIF 
351       P(I,5)=ULMASS(K(I,2)) 
352       CALL LUPTDI(KFL(1),PX(3),PY(3)) 
353       PR(1)=P(I,5)**2+(PX(1)+PX(3))**2+(PY(1)+PY(3))**2 
354       CALL LUZDIS(KFL(1),KFL(3),PR(1),Z) 
355       IF(IABS(KFL(1)).GE.4.AND.IABS(KFL(1)).LE.8.AND. 
356      &MSTU(90).LT.8) THEN 
357         MSTU(90)=MSTU(90)+1 
358         MSTU(90+MSTU(90))=I 
359         PARU(90+MSTU(90))=Z 
360       ENDIF 
361       GAM(3)=(1.-Z)*(GAM(1)+PR(1)/Z) 
362       DO 400 J=1,3 
363       IN(J)=IN(3+J) 
364   400 CONTINUE 
365  
366 C...Junction strings: stepping within or from 'low' string region easy. 
367       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)* 
368      &P(IN(1),5)**2.GE.PR(1)) THEN 
369         P(IN(1)+2,4)=Z*P(IN(1)+2,3) 
370         P(IN(2)+2,4)=PR(1)/(P(IN(1)+2,4)*P(IN(1),5)**2) 
371         DO 410 J=1,4 
372         P(I,J)=(PX(1)+PX(3))*P(IN(3),J)+(PY(1)+PY(3))*P(IN(3)+1,J) 
373   410   CONTINUE 
374         GOTO 500 
375       ELSEIF(IN(1)+1.EQ.IN(2)) THEN 
376         P(IN(2)+2,4)=P(IN(2)+2,3) 
377         P(IN(2)+2,1)=1. 
378         IN(2)=IN(2)+4 
379         IF(IN(2).GT.N+NR+4*NS) GOTO 320 
380         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN 
381           P(IN(1)+2,4)=P(IN(1)+2,3) 
382           P(IN(1)+2,1)=0. 
383           IN(1)=IN(1)+4 
384         ENDIF 
385       ENDIF 
386  
387 C...Junction strings: find new transverse directions. 
388   420 IF(IN(1).GT.N+NR+4*NS.OR.IN(2).GT.N+NR+4*NS.OR. 
389      &IN(1).GT.IN(2)) GOTO 320 
390       IF(IN(1).NE.IN(4).OR.IN(2).NE.IN(5)) THEN 
391         DO 430 J=1,4 
392         DP(1,J)=P(IN(1),J) 
393         DP(2,J)=P(IN(2),J) 
394         DP(3,J)=0. 
395         DP(4,J)=0. 
396   430   CONTINUE 
397         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) 
398         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) 
399         DHC12=DFOUR(1,2) 
400         IF(DHC12.LE.1E-2) THEN 
401           P(IN(1)+2,4)=P(IN(1)+2,3) 
402           P(IN(1)+2,1)=0. 
403           IN(1)=IN(1)+4 
404           GOTO 420 
405         ENDIF 
406         IN(3)=N+NR+4*NS+5 
407         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
408         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
409         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
410         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1. 
411         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1. 
412         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1. 
413         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1. 
414         DHCX1=DFOUR(3,1)/DHC12 
415         DHCX2=DFOUR(3,2)/DHC12 
416         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12) 
417         DHCY1=DFOUR(4,1)/DHC12 
418         DHCY2=DFOUR(4,2)/DHC12 
419         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
420         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2) 
421         DO 440 J=1,4 
422         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
423         P(IN(3),J)=DP(3,J) 
424         P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)- 
425      &  DHCYX*DP(3,J)) 
426   440   CONTINUE 
427 C...Express pT with respect to new axes, if sensible. 
428         PXP=-(PX(3)*FOUR(IN(6),IN(3))+PY(3)*FOUR(IN(6)+1,IN(3))) 
429         PYP=-(PX(3)*FOUR(IN(6),IN(3)+1)+PY(3)*FOUR(IN(6)+1,IN(3)+1)) 
430         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN 
431           PX(3)=PXP 
432           PY(3)=PYP 
433         ENDIF 
434       ENDIF 
435  
436 C...Junction strings: sum up known four-momentum, coefficients for m2. 
437       DO 470 J=1,4 
438       DHG(J)=0. 
439       P(I,J)=PX(1)*P(IN(6),J)+PY(1)*P(IN(6)+1,J)+PX(3)*P(IN(3),J)+ 
440      &PY(3)*P(IN(3)+1,J) 
441       DO 450 IN1=IN(4),IN(1)-4,4 
442       P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J) 
443   450 CONTINUE 
444       DO 460 IN2=IN(5),IN(2)-4,4 
445       P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J) 
446   460 CONTINUE 
447   470 CONTINUE 
448       DHM(1)=FOUR(I,I) 
449       DHM(2)=2.*FOUR(I,IN(1)) 
450       DHM(3)=2.*FOUR(I,IN(2)) 
451       DHM(4)=2.*FOUR(IN(1),IN(2)) 
452  
453 C...Junction strings: find coefficients for Gamma expression. 
454       DO 490 IN2=IN(1)+1,IN(2),4 
455       DO 480 IN1=IN(1),IN2-1,4 
456       DHC=2.*FOUR(IN1,IN2) 
457       DHG(1)=DHG(1)+P(IN1+2,1)*P(IN2+2,1)*DHC 
458       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-P(IN2+2,1)*DHC 
459       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+P(IN1+2,1)*DHC 
460       IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC 
461   480 CONTINUE 
462   490 CONTINUE 
463  
464 C...Junction strings: solve (m2, Gamma) equation system for energies. 
465       DHS1=DHM(3)*DHG(4)-DHM(4)*DHG(3) 
466       IF(ABS(DHS1).LT.1E-4) GOTO 320 
467       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(2)*DHG(3)-DHG(4)* 
468      &(P(I,5)**2-DHM(1))+DHG(2)*DHM(3) 
469       DHS3=DHM(2)*(GAM(3)-DHG(1))-DHG(2)*(P(I,5)**2-DHM(1)) 
470       P(IN(2)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)- 
471      &DHS2/DHS1) 
472       IF(DHM(2)+DHM(4)*P(IN(2)+2,4).LE.0.) GOTO 320 
473       P(IN(1)+2,4)=(P(I,5)**2-DHM(1)-DHM(3)*P(IN(2)+2,4))/ 
474      &(DHM(2)+DHM(4)*P(IN(2)+2,4)) 
475  
476 C...Junction strings: step to new region if necessary. 
477       IF(P(IN(2)+2,4).GT.P(IN(2)+2,3)) THEN 
478         P(IN(2)+2,4)=P(IN(2)+2,3) 
479         P(IN(2)+2,1)=1. 
480         IN(2)=IN(2)+4 
481         IF(IN(2).GT.N+NR+4*NS) GOTO 320 
482         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN 
483           P(IN(1)+2,4)=P(IN(1)+2,3) 
484           P(IN(1)+2,1)=0. 
485           IN(1)=IN(1)+4 
486         ENDIF 
487         GOTO 420 
488       ELSEIF(P(IN(1)+2,4).GT.P(IN(1)+2,3)) THEN 
489         P(IN(1)+2,4)=P(IN(1)+2,3) 
490         P(IN(1)+2,1)=0. 
491         IN(1)=IN(1)+JS 
492         GOTO 820 
493       ENDIF 
494  
495 C...Junction strings: particle four-momentum, remainder, loop back. 
496   500 DO 510 J=1,4 
497       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J) 
498       PJU(IU+3,J)=PJU(IU+3,J)+P(I,J) 
499   510 CONTINUE 
500       IF(P(I,4).LT.P(I,5)) GOTO 320 
501       PJU(IU+3,5)=TJU(4)*PJU(IU+3,4)-TJU(1)*PJU(IU+3,1)- 
502      &TJU(2)*PJU(IU+3,2)-TJU(3)*PJU(IU+3,3) 
503       IF(PJU(IU+3,5).LT.PJU(IU,5)) THEN 
504         KFL(1)=-KFL(3) 
505         PX(1)=-PX(3) 
506         PY(1)=-PY(3) 
507         GAM(1)=GAM(3) 
508         IF(IN(3).NE.IN(6)) THEN 
509           DO 520 J=1,4 
510           P(IN(6),J)=P(IN(3),J) 
511           P(IN(6)+1,J)=P(IN(3)+1,J) 
512   520     CONTINUE 
513         ENDIF 
514         DO 530 JQ=1,2 
515         IN(3+JQ)=IN(JQ) 
516         P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4) 
517         P(IN(JQ)+2,1)=P(IN(JQ)+2,1)-(3-2*JQ)*P(IN(JQ)+2,4) 
518   530   CONTINUE 
519         GOTO 380 
520       ENDIF 
521  
522 C...Junction strings: save quantities left after each string. 
523       IF(IABS(KFL(1)).GT.10) GOTO 320 
524       I=I-1 
525       KFJH(IU)=KFL(1) 
526       DO 540 J=1,4 
527       PJU(IU+3,J)=PJU(IU+3,J)-P(I+1,J) 
528   540 CONTINUE 
529   550 CONTINUE 
530  
531 C...Junction strings: put together to new effective string endpoint. 
532       NJS(JT)=I-ISTA 
533       KFJS(JT)=K(K(MJU(JT+2),3),2) 
534       KFLS=2*INT(RLU(0)+3.*PARJ(4)/(1.+3.*PARJ(4)))+1 
535       IF(KFJH(1).EQ.KFJH(2)) KFLS=3 
536       IF(ISTA.NE.I) KFJS(JT)=ISIGN(1000*MAX(IABS(KFJH(1)), 
537      &IABS(KFJH(2)))+100*MIN(IABS(KFJH(1)),IABS(KFJH(2)))+ 
538      &KFLS,KFJH(1)) 
539       DO 560 J=1,4 
540       PJS(JT,J)=PJU(1,J)+PJU(2,J)+P(MJU(JT),J) 
541       PJS(JT+2,J)=PJU(4,J)+PJU(5,J) 
542   560 CONTINUE 
543       PJS(JT,5)=SQRT(MAX(0.,PJS(JT,4)**2-PJS(JT,1)**2-PJS(JT,2)**2- 
544      &PJS(JT,3)**2)) 
545   570 CONTINUE 
546  
547 C...Open versus closed strings. Choose breakup region for latter. 
548   580 IF(MJU(1).NE.0.AND.MJU(2).NE.0) THEN 
549         NS=MJU(2)-MJU(1) 
550         NB=MJU(1)-N 
551       ELSEIF(MJU(1).NE.0) THEN 
552         NS=N+NR-MJU(1) 
553         NB=MJU(1)-N 
554       ELSEIF(MJU(2).NE.0) THEN 
555         NS=MJU(2)-N 
556         NB=1 
557       ELSEIF(IABS(K(N+1,2)).NE.21) THEN 
558         NS=NR-1 
559         NB=1 
560       ELSE 
561         NS=NR+1 
562         W2SUM=0. 
563         DO 590 IS=1,NR 
564         P(N+NR+IS,1)=0.5*FOUR(N+IS,N+IS+1-NR*(IS/NR)) 
565         W2SUM=W2SUM+P(N+NR+IS,1) 
566   590   CONTINUE 
567         W2RAN=RLU(0)*W2SUM 
568         NB=0 
569   600   NB=NB+1 
570         W2SUM=W2SUM-P(N+NR+NB,1) 
571         IF(W2SUM.GT.W2RAN.AND.NB.LT.NR) GOTO 600 
572       ENDIF 
573  
574 C...Find longitudinal string directions (i.e. lightlike four-vectors). 
575       DO 630 IS=1,NS 
576       IS1=N+IS+NB-1-NR*((IS+NB-2)/NR) 
577       IS2=N+IS+NB-NR*((IS+NB-1)/NR) 
578       DO 610 J=1,5 
579       DP(1,J)=P(IS1,J) 
580       IF(IABS(K(IS1,2)).EQ.21) DP(1,J)=0.5*DP(1,J) 
581       IF(IS1.EQ.MJU(1)) DP(1,J)=PJS(1,J)-PJS(3,J) 
582       DP(2,J)=P(IS2,J) 
583       IF(IABS(K(IS2,2)).EQ.21) DP(2,J)=0.5*DP(2,J) 
584       IF(IS2.EQ.MJU(2)) DP(2,J)=PJS(2,J)-PJS(4,J) 
585   610 CONTINUE 
586       DP(3,5)=DFOUR(1,1) 
587       DP(4,5)=DFOUR(2,2) 
588       DHKC=DFOUR(1,2) 
589       IF(DP(3,5)+2.*DHKC+DP(4,5).LE.0.) THEN 
590         DP(3,5)=DP(1,5)**2 
591         DP(4,5)=DP(2,5)**2 
592         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2+DP(1,5)**2) 
593         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2+DP(2,5)**2) 
594         DHKC=DFOUR(1,2) 
595       ENDIF 
596       DHKS=SQRT(DHKC**2-DP(3,5)*DP(4,5)) 
597       DHK1=0.5*((DP(4,5)+DHKC)/DHKS-1.) 
598       DHK2=0.5*((DP(3,5)+DHKC)/DHKS-1.) 
599       IN1=N+NR+4*IS-3 
600       P(IN1,5)=SQRT(DP(3,5)+2.*DHKC+DP(4,5)) 
601       DO 620 J=1,4 
602       P(IN1,J)=(1.+DHK1)*DP(1,J)-DHK2*DP(2,J) 
603       P(IN1+1,J)=(1.+DHK2)*DP(2,J)-DHK1*DP(1,J) 
604   620 CONTINUE 
605   630 CONTINUE 
606  
607 C...Begin initialization: sum up energy, set starting position. 
608       ISAV=I 
609       MSTU91=MSTU(90) 
610   640 NTRY=NTRY+1 
611       IF(NTRY.GT.100.AND.NTRYR.LE.4) THEN 
612         PARU12=4.*PARU12 
613         PARU13=2.*PARU13 
614         GOTO 140 
615       ELSEIF(NTRY.GT.100) THEN 
616         CALL LUERRM(14,'(LUSTRF:) caught in infinite loop') 
617         IF(MSTU(21).GE.1) RETURN 
618       ENDIF 
619       I=ISAV 
620       MSTU(90)=MSTU91 
621       DO 660 J=1,4 
622       P(N+NRS,J)=0. 
623       DO 650 IS=1,NR 
624       P(N+NRS,J)=P(N+NRS,J)+P(N+IS,J) 
625   650 CONTINUE 
626   660 CONTINUE 
627       DO 680 JT=1,2 
628       IRANK(JT)=0 
629       IF(MJU(JT).NE.0) IRANK(JT)=NJS(JT) 
630       IF(NS.GT.NR) IRANK(JT)=1 
631       IE(JT)=K(N+1+(JT/2)*(NP-1),3) 
632       IN(3*JT+1)=N+NR+1+4*(JT/2)*(NS-1) 
633       IN(3*JT+2)=IN(3*JT+1)+1 
634       IN(3*JT+3)=N+NR+4*NS+2*JT-1 
635       DO 670 IN1=N+NR+2+JT,N+NR+4*NS-2+JT,4 
636       P(IN1,1)=2-JT 
637       P(IN1,2)=JT-1 
638       P(IN1,3)=1. 
639   670 CONTINUE 
640   680 CONTINUE 
641  
642 C...Initialize flavour and pT variables for open string. 
643       IF(NS.LT.NR) THEN 
644         PX(1)=0. 
645         PY(1)=0. 
646         IF(NS.EQ.1.AND.MJU(1)+MJU(2).EQ.0) CALL LUPTDI(0,PX(1),PY(1)) 
647         PX(2)=-PX(1) 
648         PY(2)=-PY(1) 
649         DO 690 JT=1,2 
650         KFL(JT)=K(IE(JT),2) 
651         IF(MJU(JT).NE.0) KFL(JT)=KFJS(JT) 
652         MSTJ(93)=1 
653         PMQ(JT)=ULMASS(KFL(JT)) 
654         GAM(JT)=0. 
655   690   CONTINUE 
656  
657 C...Closed string: random initial breakup flavour, pT and vertex. 
658       ELSE 
659         KFL(3)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5) 
660         CALL LUKFDI(KFL(3),0,KFL(1),KDUMP) 
661         KFL(2)=-KFL(1) 
662         IF(IABS(KFL(1)).GT.10.AND.RLU(0).GT.0.5) THEN 
663           KFL(2)=-(KFL(1)+ISIGN(10000,KFL(1))) 
664         ELSEIF(IABS(KFL(1)).GT.10) THEN 
665           KFL(1)=-(KFL(2)+ISIGN(10000,KFL(2))) 
666         ENDIF 
667         CALL LUPTDI(KFL(1),PX(1),PY(1)) 
668         PX(2)=-PX(1) 
669         PY(2)=-PY(1) 
670         PR3=MIN(25.,0.1*P(N+NR+1,5)**2) 
671   700   CALL LUZDIS(KFL(1),KFL(2),PR3,Z) 
672         ZR=PR3/(Z*P(N+NR+1,5)**2) 
673         IF(ZR.GE.1.) GOTO 700 
674         DO 710 JT=1,2 
675         MSTJ(93)=1 
676         PMQ(JT)=ULMASS(KFL(JT)) 
677         GAM(JT)=PR3*(1.-Z)/Z 
678         IN1=N+NR+3+4*(JT/2)*(NS-1) 
679         P(IN1,JT)=1.-Z 
680         P(IN1,3-JT)=JT-1 
681         P(IN1,3)=(2-JT)*(1.-Z)+(JT-1)*Z 
682         P(IN1+1,JT)=ZR 
683         P(IN1+1,3-JT)=2-JT 
684         P(IN1+1,3)=(2-JT)*(1.-ZR)+(JT-1)*ZR 
685   710   CONTINUE 
686       ENDIF 
687  
688 C...Find initial transverse directions (i.e. spacelike four-vectors). 
689       DO 750 JT=1,2 
690       IF(JT.EQ.1.OR.NS.EQ.NR-1) THEN 
691         IN1=IN(3*JT+1) 
692         IN3=IN(3*JT+3) 
693         DO 720 J=1,4 
694         DP(1,J)=P(IN1,J) 
695         DP(2,J)=P(IN1+1,J) 
696         DP(3,J)=0. 
697         DP(4,J)=0. 
698   720   CONTINUE 
699         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) 
700         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) 
701         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
702         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
703         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
704         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1. 
705         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1. 
706         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1. 
707         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1. 
708         DHC12=DFOUR(1,2) 
709         DHCX1=DFOUR(3,1)/DHC12 
710         DHCX2=DFOUR(3,2)/DHC12 
711         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12) 
712         DHCY1=DFOUR(4,1)/DHC12 
713         DHCY2=DFOUR(4,2)/DHC12 
714         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
715         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2) 
716         DO 730 J=1,4 
717         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
718         P(IN3,J)=DP(3,J) 
719         P(IN3+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)- 
720      &  DHCYX*DP(3,J)) 
721   730   CONTINUE 
722       ELSE 
723         DO 740 J=1,4 
724         P(IN3+2,J)=P(IN3,J) 
725         P(IN3+3,J)=P(IN3+1,J) 
726   740   CONTINUE 
727       ENDIF 
728   750 CONTINUE 
729  
730 C...Remove energy used up in junction string fragmentation. 
731       IF(MJU(1)+MJU(2).GT.0) THEN 
732         DO 770 JT=1,2 
733         IF(NJS(JT).EQ.0) GOTO 770 
734         DO 760 J=1,4 
735         P(N+NRS,J)=P(N+NRS,J)-PJS(JT+2,J) 
736   760   CONTINUE 
737   770   CONTINUE 
738       ENDIF 
739  
740 C...Produce new particle: side, origin. 
741   780 I=I+1 
742       IF(2*I-NSAV.GE.MSTU(4)-MSTU(32)-5) THEN 
743         CALL LUERRM(11,'(LUSTRF:) no more memory left in LUJETS') 
744         IF(MSTU(21).GE.1) RETURN 
745       ENDIF 
746       JT=1.5+RLU(0) 
747       IF(IABS(KFL(3-JT)).GT.10) JT=3-JT 
748       IF(IABS(KFL(3-JT)).GE.4.AND.IABS(KFL(3-JT)).LE.8) JT=3-JT 
749       JR=3-JT 
750       JS=3-2*JT 
751       IRANK(JT)=IRANK(JT)+1 
752       K(I,1)=1 
753       K(I,3)=IE(JT) 
754       K(I,4)=0 
755       K(I,5)=0 
756  
757 C...Generate flavour, hadron and pT. 
758   790 CALL LUKFDI(KFL(JT),0,KFL(3),K(I,2)) 
759       IF(K(I,2).EQ.0) GOTO 640 
760       IF(MSTJ(12).GE.3.AND.IRANK(JT).EQ.1.AND.IABS(KFL(JT)).LE.10.AND. 
761      &IABS(KFL(3)).GT.10) THEN 
762         IF(RLU(0).GT.PARJ(19)) GOTO 790 
763       ENDIF 
764       P(I,5)=ULMASS(K(I,2)) 
765       CALL LUPTDI(KFL(JT),PX(3),PY(3)) 
766       PR(JT)=P(I,5)**2+(PX(JT)+PX(3))**2+(PY(JT)+PY(3))**2 
767  
768 C...Final hadrons for small invariant mass. 
769       MSTJ(93)=1 
770       PMQ(3)=ULMASS(KFL(3)) 
771       PARJST=PARJ(33) 
772       IF(MSTJ(11).EQ.2) PARJST=PARJ(34) 
773       WMIN=PARJST+PMQ(1)+PMQ(2)+PARJ(36)*PMQ(3) 
774       IF(IABS(KFL(JT)).GT.10.AND.IABS(KFL(3)).GT.10) WMIN= 
775      &WMIN-0.5*PARJ(36)*PMQ(3) 
776       WREM2=FOUR(N+NRS,N+NRS) 
777       IF(WREM2.LT.0.10) GOTO 640 
778       IF(WREM2.LT.MAX(WMIN*(1.+(2.*RLU(0)-1.)*PARJ(37)), 
779      &PARJ(32)+PMQ(1)+PMQ(2))**2) GOTO 940 
780  
781 C...Choose z, which gives Gamma. Shift z for heavy flavours. 
782       CALL LUZDIS(KFL(JT),KFL(3),PR(JT),Z) 
783       IF(IABS(KFL(JT)).GE.4.AND.IABS(KFL(JT)).LE.8.AND. 
784      &MSTU(90).LT.8) THEN 
785         MSTU(90)=MSTU(90)+1 
786         MSTU(90+MSTU(90))=I 
787         PARU(90+MSTU(90))=Z 
788       ENDIF 
789       KFL1A=IABS(KFL(1)) 
790       KFL2A=IABS(KFL(2)) 
791       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10), 
792      &MOD(KFL2A/1000,10)).GE.4) THEN 
793         PR(JR)=(PMQ(JR)+PMQ(3))**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2 
794         PW12=SQRT(MAX(0.,(WREM2-PR(1)-PR(2))**2-4.*PR(1)*PR(2))) 
795         Z=(WREM2+PR(JT)-PR(JR)+PW12*(2.*Z-1.))/(2.*WREM2) 
796         PR(JR)=(PMQ(JR)+PARJST)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2 
797         IF((1.-Z)*(WREM2-PR(JT)/Z).LT.PR(JR)) GOTO 940 
798       ENDIF 
799       GAM(3)=(1.-Z)*(GAM(JT)+PR(JT)/Z) 
800       DO 800 J=1,3 
801       IN(J)=IN(3*JT+J) 
802   800 CONTINUE 
803  
804 C...Stepping within or from 'low' string region easy. 
805       IF(IN(1)+1.EQ.IN(2).AND.Z*P(IN(1)+2,3)*P(IN(2)+2,3)* 
806      &P(IN(1),5)**2.GE.PR(JT)) THEN 
807         P(IN(JT)+2,4)=Z*P(IN(JT)+2,3) 
808         P(IN(JR)+2,4)=PR(JT)/(P(IN(JT)+2,4)*P(IN(1),5)**2) 
809         DO 810 J=1,4 
810         P(I,J)=(PX(JT)+PX(3))*P(IN(3),J)+(PY(JT)+PY(3))*P(IN(3)+1,J) 
811   810   CONTINUE 
812         GOTO 900 
813       ELSEIF(IN(1)+1.EQ.IN(2)) THEN 
814         P(IN(JR)+2,4)=P(IN(JR)+2,3) 
815         P(IN(JR)+2,JT)=1. 
816         IN(JR)=IN(JR)+4*JS 
817         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 640 
818         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN 
819           P(IN(JT)+2,4)=P(IN(JT)+2,3) 
820           P(IN(JT)+2,JT)=0. 
821           IN(JT)=IN(JT)+4*JS 
822         ENDIF 
823       ENDIF 
824  
825 C...Find new transverse directions (i.e. spacelike string vectors). 
826   820 IF(JS*IN(1).GT.JS*IN(3*JR+1).OR.JS*IN(2).GT.JS*IN(3*JR+2).OR. 
827      &IN(1).GT.IN(2)) GOTO 640 
828       IF(IN(1).NE.IN(3*JT+1).OR.IN(2).NE.IN(3*JT+2)) THEN 
829         DO 830 J=1,4 
830         DP(1,J)=P(IN(1),J) 
831         DP(2,J)=P(IN(2),J) 
832         DP(3,J)=0. 
833         DP(4,J)=0. 
834   830   CONTINUE 
835         DP(1,4)=SQRT(DP(1,1)**2+DP(1,2)**2+DP(1,3)**2) 
836         DP(2,4)=SQRT(DP(2,1)**2+DP(2,2)**2+DP(2,3)**2) 
837         DHC12=DFOUR(1,2) 
838         IF(DHC12.LE.1E-2) THEN 
839           P(IN(JT)+2,4)=P(IN(JT)+2,3) 
840           P(IN(JT)+2,JT)=0. 
841           IN(JT)=IN(JT)+4*JS 
842           GOTO 820 
843         ENDIF 
844         IN(3)=N+NR+4*NS+5 
845         DP(5,1)=DP(1,1)/DP(1,4)-DP(2,1)/DP(2,4) 
846         DP(5,2)=DP(1,2)/DP(1,4)-DP(2,2)/DP(2,4) 
847         DP(5,3)=DP(1,3)/DP(1,4)-DP(2,3)/DP(2,4) 
848         IF(DP(5,1)**2.LE.DP(5,2)**2+DP(5,3)**2) DP(3,1)=1. 
849         IF(DP(5,1)**2.GT.DP(5,2)**2+DP(5,3)**2) DP(3,3)=1. 
850         IF(DP(5,2)**2.LE.DP(5,1)**2+DP(5,3)**2) DP(4,2)=1. 
851         IF(DP(5,2)**2.GT.DP(5,1)**2+DP(5,3)**2) DP(4,3)=1. 
852         DHCX1=DFOUR(3,1)/DHC12 
853         DHCX2=DFOUR(3,2)/DHC12 
854         DHCXX=1D0/SQRT(1D0+2D0*DHCX1*DHCX2*DHC12) 
855         DHCY1=DFOUR(4,1)/DHC12 
856         DHCY2=DFOUR(4,2)/DHC12 
857         DHCYX=DHCXX*(DHCX1*DHCY2+DHCX2*DHCY1)*DHC12 
858         DHCYY=1D0/SQRT(1D0+2D0*DHCY1*DHCY2*DHC12-DHCYX**2) 
859         DO 840 J=1,4 
860         DP(3,J)=DHCXX*(DP(3,J)-DHCX2*DP(1,J)-DHCX1*DP(2,J)) 
861         P(IN(3),J)=DP(3,J) 
862         P(IN(3)+1,J)=DHCYY*(DP(4,J)-DHCY2*DP(1,J)-DHCY1*DP(2,J)- 
863      &  DHCYX*DP(3,J)) 
864   840   CONTINUE 
865 C...Express pT with respect to new axes, if sensible. 
866         PXP=-(PX(3)*FOUR(IN(3*JT+3),IN(3))+PY(3)* 
867      &  FOUR(IN(3*JT+3)+1,IN(3))) 
868         PYP=-(PX(3)*FOUR(IN(3*JT+3),IN(3)+1)+PY(3)* 
869      &  FOUR(IN(3*JT+3)+1,IN(3)+1)) 
870         IF(ABS(PXP**2+PYP**2-PX(3)**2-PY(3)**2).LT.0.01) THEN 
871           PX(3)=PXP 
872           PY(3)=PYP 
873         ENDIF 
874       ENDIF 
875  
876 C...Sum up known four-momentum. Gives coefficients for m2 expression. 
877       DO 870 J=1,4 
878       DHG(J)=0. 
879       P(I,J)=PX(JT)*P(IN(3*JT+3),J)+PY(JT)*P(IN(3*JT+3)+1,J)+ 
880      &PX(3)*P(IN(3),J)+PY(3)*P(IN(3)+1,J) 
881       DO 850 IN1=IN(3*JT+1),IN(1)-4*JS,4*JS 
882       P(I,J)=P(I,J)+P(IN1+2,3)*P(IN1,J) 
883   850 CONTINUE 
884       DO 860 IN2=IN(3*JT+2),IN(2)-4*JS,4*JS 
885       P(I,J)=P(I,J)+P(IN2+2,3)*P(IN2,J) 
886   860 CONTINUE 
887   870 CONTINUE 
888       DHM(1)=FOUR(I,I) 
889       DHM(2)=2.*FOUR(I,IN(1)) 
890       DHM(3)=2.*FOUR(I,IN(2)) 
891       DHM(4)=2.*FOUR(IN(1),IN(2)) 
892  
893 C...Find coefficients for Gamma expression. 
894       DO 890 IN2=IN(1)+1,IN(2),4 
895       DO 880 IN1=IN(1),IN2-1,4 
896       DHC=2.*FOUR(IN1,IN2) 
897       DHG(1)=DHG(1)+P(IN1+2,JT)*P(IN2+2,JT)*DHC 
898       IF(IN1.EQ.IN(1)) DHG(2)=DHG(2)-JS*P(IN2+2,JT)*DHC 
899       IF(IN2.EQ.IN(2)) DHG(3)=DHG(3)+JS*P(IN1+2,JT)*DHC 
900       IF(IN1.EQ.IN(1).AND.IN2.EQ.IN(2)) DHG(4)=DHG(4)-DHC 
901   880 CONTINUE 
902   890 CONTINUE 
903  
904 C...Solve (m2, Gamma) equation system for energies taken. 
905       DHS1=DHM(JR+1)*DHG(4)-DHM(4)*DHG(JR+1) 
906       IF(ABS(DHS1).LT.1E-4) GOTO 640 
907       DHS2=DHM(4)*(GAM(3)-DHG(1))-DHM(JT+1)*DHG(JR+1)-DHG(4)* 
908      &(P(I,5)**2-DHM(1))+DHG(JT+1)*DHM(JR+1) 
909       DHS3=DHM(JT+1)*(GAM(3)-DHG(1))-DHG(JT+1)*(P(I,5)**2-DHM(1)) 
910       P(IN(JR)+2,4)=0.5*(SQRT(MAX(0D0,DHS2**2-4.*DHS1*DHS3))/ABS(DHS1)- 
911      &DHS2/DHS1) 
912       IF(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4).LE.0.) GOTO 640 
913       P(IN(JT)+2,4)=(P(I,5)**2-DHM(1)-DHM(JR+1)*P(IN(JR)+2,4))/ 
914      &(DHM(JT+1)+DHM(4)*P(IN(JR)+2,4)) 
915  
916 C...Step to new region if necessary. 
917       IF(P(IN(JR)+2,4).GT.P(IN(JR)+2,3)) THEN 
918         P(IN(JR)+2,4)=P(IN(JR)+2,3) 
919         P(IN(JR)+2,JT)=1. 
920         IN(JR)=IN(JR)+4*JS 
921         IF(JS*IN(JR).GT.JS*IN(4*JR)) GOTO 640 
922         IF(FOUR(IN(1),IN(2)).LE.1E-2) THEN 
923           P(IN(JT)+2,4)=P(IN(JT)+2,3) 
924           P(IN(JT)+2,JT)=0. 
925           IN(JT)=IN(JT)+4*JS 
926         ENDIF 
927         GOTO 820 
928       ELSEIF(P(IN(JT)+2,4).GT.P(IN(JT)+2,3)) THEN 
929         P(IN(JT)+2,4)=P(IN(JT)+2,3) 
930         P(IN(JT)+2,JT)=0. 
931         IN(JT)=IN(JT)+4*JS 
932         GOTO 820 
933       ENDIF 
934  
935 C...Four-momentum of particle. Remaining quantities. Loop back. 
936   900 DO 910 J=1,4 
937       P(I,J)=P(I,J)+P(IN(1)+2,4)*P(IN(1),J)+P(IN(2)+2,4)*P(IN(2),J) 
938       P(N+NRS,J)=P(N+NRS,J)-P(I,J) 
939   910 CONTINUE 
940       IF(P(I,4).LT.P(I,5)) GOTO 640 
941       KFL(JT)=-KFL(3) 
942       PMQ(JT)=PMQ(3) 
943       PX(JT)=-PX(3) 
944       PY(JT)=-PY(3) 
945       GAM(JT)=GAM(3) 
946       IF(IN(3).NE.IN(3*JT+3)) THEN 
947         DO 920 J=1,4 
948         P(IN(3*JT+3),J)=P(IN(3),J) 
949         P(IN(3*JT+3)+1,J)=P(IN(3)+1,J) 
950   920   CONTINUE 
951       ENDIF 
952       DO 930 JQ=1,2 
953       IN(3*JT+JQ)=IN(JQ) 
954       P(IN(JQ)+2,3)=P(IN(JQ)+2,3)-P(IN(JQ)+2,4) 
955       P(IN(JQ)+2,JT)=P(IN(JQ)+2,JT)-JS*(3-2*JQ)*P(IN(JQ)+2,4) 
956   930 CONTINUE 
957       GOTO 780 
958  
959 C...Final hadron: side, flavour, hadron, mass. 
960   940 I=I+1 
961       K(I,1)=1 
962       K(I,3)=IE(JR) 
963       K(I,4)=0 
964       K(I,5)=0 
965       CALL LUKFDI(KFL(JR),-KFL(3),KFLDMP,K(I,2)) 
966       IF(K(I,2).EQ.0) GOTO 640 
967       P(I,5)=ULMASS(K(I,2)) 
968       PR(JR)=P(I,5)**2+(PX(JR)-PX(3))**2+(PY(JR)-PY(3))**2 
969  
970 C...Final two hadrons: find common setup of four-vectors. 
971       JQ=1 
972       IF(P(IN(4)+2,3)*P(IN(5)+2,3)*FOUR(IN(4),IN(5)).LT.P(IN(7),3)* 
973      &P(IN(8),3)*FOUR(IN(7),IN(8))) JQ=2 
974       DHC12=FOUR(IN(3*JQ+1),IN(3*JQ+2)) 
975       DHR1=FOUR(N+NRS,IN(3*JQ+2))/DHC12 
976       DHR2=FOUR(N+NRS,IN(3*JQ+1))/DHC12 
977       IF(IN(4).NE.IN(7).OR.IN(5).NE.IN(8)) THEN 
978         PX(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3))-PX(JQ) 
979         PY(3-JQ)=-FOUR(N+NRS,IN(3*JQ+3)+1)-PY(JQ) 
980         PR(3-JQ)=P(I+(JT+JQ-3)**2-1,5)**2+(PX(3-JQ)+(2*JQ-3)*JS* 
981      &  PX(3))**2+(PY(3-JQ)+(2*JQ-3)*JS*PY(3))**2 
982       ENDIF 
983  
984 C...Solve kinematics for final two hadrons, if possible. 
985       WREM2=WREM2+(PX(1)+PX(2))**2+(PY(1)+PY(2))**2 
986       FD=(SQRT(PR(1))+SQRT(PR(2)))/SQRT(WREM2) 
987       IF(MJU(1)+MJU(2).NE.0.AND.I.EQ.ISAV+2.AND.FD.GE.1.) GOTO 200 
988       IF(FD.GE.1.) GOTO 640 
989       FA=WREM2+PR(JT)-PR(JR) 
990       IF(MSTJ(11).NE.2) PREV=0.5*EXP(MAX(-50.,LOG(FD)*PARJ(38)* 
991      &(PR(1)+PR(2))**2)) 
992       IF(MSTJ(11).EQ.2) PREV=0.5*FD**PARJ(39) 
993       FB=SIGN(SQRT(MAX(0.,FA**2-4.*WREM2*PR(JT))),JS*(RLU(0)-PREV)) 
994       KFL1A=IABS(KFL(1)) 
995       KFL2A=IABS(KFL(2)) 
996       IF(MAX(MOD(KFL1A,10),MOD(KFL1A/1000,10),MOD(KFL2A,10), 
997      &MOD(KFL2A/1000,10)).GE.6) FB=SIGN(SQRT(MAX(0.,FA**2- 
998      &4.*WREM2*PR(JT))),FLOAT(JS)) 
999       DO 950 J=1,4 
1000       P(I-1,J)=(PX(JT)+PX(3))*P(IN(3*JQ+3),J)+(PY(JT)+PY(3))* 
1001      &P(IN(3*JQ+3)+1,J)+0.5*(DHR1*(FA+FB)*P(IN(3*JQ+1),J)+ 
1002      &DHR2*(FA-FB)*P(IN(3*JQ+2),J))/WREM2 
1003       P(I,J)=P(N+NRS,J)-P(I-1,J) 
1004   950 CONTINUE 
1005       IF(P(I-1,4).LT.P(I-1,5).OR.P(I,4).LT.P(I,5)) GOTO 640 
1006  
1007 C...Mark jets as fragmented and give daughter pointers. 
1008       N=I-NRS+1 
1009       DO 960 I=NSAV+1,NSAV+NP 
1010       IM=K(I,3) 
1011       K(IM,1)=K(IM,1)+10 
1012       IF(MSTU(16).NE.2) THEN 
1013         K(IM,4)=NSAV+1 
1014         K(IM,5)=NSAV+1 
1015       ELSE 
1016         K(IM,4)=NSAV+2 
1017         K(IM,5)=N 
1018       ENDIF 
1019   960 CONTINUE 
1020  
1021 C...Document string system. Move up particles. 
1022       NSAV=NSAV+1 
1023       K(NSAV,1)=11 
1024       K(NSAV,2)=92 
1025       K(NSAV,3)=IP 
1026       K(NSAV,4)=NSAV+1 
1027       K(NSAV,5)=N 
1028       DO 970 J=1,4 
1029       P(NSAV,J)=DPS(J) 
1030       V(NSAV,J)=V(IP,J) 
1031   970 CONTINUE 
1032       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2)) 
1033       V(NSAV,5)=0. 
1034       DO 990 I=NSAV+1,N 
1035       DO 980 J=1,5 
1036       K(I,J)=K(I+NRS-1,J) 
1037       P(I,J)=P(I+NRS-1,J) 
1038       V(I,J)=0. 
1039   980 CONTINUE 
1040   990 CONTINUE 
1041       MSTU91=MSTU(90) 
1042       DO 1000 IZ=MSTU90+1,MSTU91 
1043       MSTU9T(IZ)=MSTU(90+IZ)-NRS+1-NSAV+N 
1044       PARU9T(IZ)=PARU(90+IZ) 
1045  1000 CONTINUE 
1046       MSTU(90)=MSTU90 
1047  
1048 C...Order particles in rank along the chain. Update mother pointer. 
1049       DO 1020 I=NSAV+1,N 
1050       DO 1010 J=1,5 
1051       K(I-NSAV+N,J)=K(I,J) 
1052       P(I-NSAV+N,J)=P(I,J) 
1053  1010 CONTINUE 
1054  1020 CONTINUE 
1055       I1=NSAV 
1056       DO 1050 I=N+1,2*N-NSAV 
1057       IF(K(I,3).NE.IE(1)) GOTO 1050 
1058       I1=I1+1 
1059       DO 1030 J=1,5 
1060       K(I1,J)=K(I,J) 
1061       P(I1,J)=P(I,J) 
1062  1030 CONTINUE 
1063       IF(MSTU(16).NE.2) K(I1,3)=NSAV 
1064       DO 1040 IZ=MSTU90+1,MSTU91 
1065       IF(MSTU9T(IZ).EQ.I) THEN 
1066         MSTU(90)=MSTU(90)+1 
1067         MSTU(90+MSTU(90))=I1 
1068         PARU(90+MSTU(90))=PARU9T(IZ) 
1069       ENDIF 
1070  1040 CONTINUE 
1071  1050 CONTINUE 
1072       DO 1080 I=2*N-NSAV,N+1,-1 
1073       IF(K(I,3).EQ.IE(1)) GOTO 1080 
1074       I1=I1+1 
1075       DO 1060 J=1,5 
1076       K(I1,J)=K(I,J) 
1077       P(I1,J)=P(I,J) 
1078  1060 CONTINUE 
1079       IF(MSTU(16).NE.2) K(I1,3)=NSAV 
1080       DO 1070 IZ=MSTU90+1,MSTU91 
1081       IF(MSTU9T(IZ).EQ.I) THEN 
1082         MSTU(90)=MSTU(90)+1 
1083         MSTU(90+MSTU(90))=I1 
1084         PARU(90+MSTU(90))=PARU9T(IZ) 
1085       ENDIF 
1086  1070 CONTINUE 
1087  1080 CONTINUE 
1088  
1089 C...Boost back particle system. Set production vertices. 
1090       IF(MBST.EQ.0) THEN 
1091         MSTU(33)=1 
1092         CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),DPS(2)/DPS(4), 
1093      &  DPS(3)/DPS(4)) 
1094       ELSE 
1095         DO 1090 I=NSAV+1,N 
1096         HHPMT=P(I,1)**2+P(I,2)**2+P(I,5)**2 
1097         IF(P(I,3).GT.0.) THEN 
1098           HHPEZ=(P(I,4)+P(I,3))*HHBZ 
1099           P(I,3)=0.5*(HHPEZ-HHPMT/HHPEZ) 
1100           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ) 
1101         ELSE 
1102           HHPEZ=(P(I,4)-P(I,3))/HHBZ 
1103           P(I,3)=-0.5*(HHPEZ-HHPMT/HHPEZ) 
1104           P(I,4)=0.5*(HHPEZ+HHPMT/HHPEZ) 
1105         ENDIF 
1106  1090   CONTINUE 
1107       ENDIF 
1108       DO 1110 I=NSAV+1,N 
1109       DO 1100 J=1,4 
1110       V(I,J)=V(IP,J) 
1111  1100 CONTINUE 
1112  1110 CONTINUE 
1113  
1114       RETURN 
1115       END