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