]>
Commit | Line | Data |
---|---|---|
e74335a4 | 1 | * $Id$ |
2 | ||
3 | C********************************************************************* | |
4 | ||
5 | SUBROUTINE LUSTRF_HIJING(IP) | |
6 | C...Purpose: to handle the fragmentation of an arbitrary colour singlet | |
7 | C...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 | ||
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 | 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 | ||
48 | C...Take copy of partons to be considered. Check flavour sum. | |
49 | NP=NP+1 | |
50 | DO 120 J=1,5 | |
51 | K(N+NP,J)=K(I,J) | |
52 | P(N+NP,J)=P(I,J) | |
53 | 120 DPS(J)=DPS(J)+P(I,J) | |
54 | K(N+NP,3)=I | |
55 | IF(P(N+NP,4)**2.LT.P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2) THEN | |
56 | P(N+NP,4)=SQRT(P(N+NP,1)**2+P(N+NP,2)**2+P(N+NP,3)**2+ | |
57 | & P(N+NP,5)**2) | |
58 | DPS(4)=DPS(4)+MAX(0.,P(N+NP,4)-P(I,4)) | |
59 | ENDIF | |
60 | IF(KQ.NE.2) KQSUM=KQSUM+KQ | |
61 | IF(K(I,1).EQ.41) THEN | |
62 | KQSUM=KQSUM+2*KQ | |
63 | IF(KQSUM.EQ.KQ) MJU(1)=N+NP | |
64 | IF(KQSUM.NE.KQ) MJU(2)=N+NP | |
65 | ENDIF | |
66 | IF(K(I,1).EQ.2.OR.K(I,1).EQ.41) GOTO 110 | |
67 | IF(KQSUM.NE.0) THEN | |
68 | CALL LUERRM_HIJING(12 | |
69 | $ ,'(LUSTRF_HIJING:) unphysical flavour combination') | |
70 | IF(MSTU(21).GE.1) RETURN | |
71 | ENDIF | |
72 | ||
73 | C...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 | ||
77 | C...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 | ||
104 | C...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 | ||
130 | C...Reset particle counter. Skip ahead if no junctions are present; | |
131 | C...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 | ||
151 | C...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 | ||
173 | C...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 | ||
192 | C...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 | ||
202 | C...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 | ||
207 | C...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 | ||
237 | C...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 | ||
267 | C...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 | ||
296 | C...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 | ||
309 | C...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 | ||
324 | C...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 | ||
344 | C...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)) | |
382 | C...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 | ||
391 | C...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 | ||
405 | C...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 | ||
414 | C...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 | ||
426 | C...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 | ||
445 | C...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 | ||
469 | C...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 | ||
477 | C...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 | ||
492 | C...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 | ||
518 | C...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 | ||
548 | C...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 | ||
578 | C...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 | ||
593 | C...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 | ||
624 | C...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 | ||
663 | C...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 | ||
672 | C...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 | ||
689 | C...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 | ||
700 | C...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 | ||
711 | C...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 | ||
728 | C...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 | ||
748 | C...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)) | |
786 | C...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 | ||
797 | C...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 | ||
811 | C...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 | ||
820 | C...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 | ||
832 | C...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 | ||
851 | C...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 | ||
872 | C...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 | ||
883 | C...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 | ||
897 | C...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 | ||
919 | C...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 | ||
933 | C...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 | ||
951 | C...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 | ||
974 | C...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 |