]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/jetset/lustrf.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Mostly minor style modifications to be ready for cloning with EMCAL
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lustrf.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LUSTRF(IP)
5C...Purpose: to handle the fragmentation of an arbitrary colour singlet
6C...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
16C...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
21C...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
48C...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
70C...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
93C...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
120C...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
150C...Reset particle counter. Skip ahead if no junctions are present;
151C...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
171C...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
197C...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
218C...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
229C...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
234C...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
267C...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
301C...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
332C...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
344C...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
366C...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
387C...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
427C...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
436C...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
453C...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
464C...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
476C...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
495C...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
522C...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
531C...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
547C...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
574C...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
607C...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
642C...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
657C...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
688C...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
730C...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
740C...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
757C...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
768C...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
781C...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
804C...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
825C...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
865C...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
876C...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
893C...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
904C...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
916C...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
935C...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
959C...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
970C...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
984C...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
1007C...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
1021C...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
1048C...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
1089C...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