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