]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | |
2 | C********************************************************************* | |
3 | ||
4 | SUBROUTINE LUINDF(IP) | |
5 | ||
6 | C...Purpose: to handle the fragmentation of a jet system (or a single | |
7 | C...jet) according to independent fragmentation models. | |
8 | IMPLICIT DOUBLE PRECISION(D) | |
9 | COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) | |
10 | COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) | |
11 | COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4) | |
12 | SAVE /LUJETS/,/LUDAT1/,/LUDAT2/ | |
13 | DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3), | |
14 | &KFLO(2),PXO(2),PYO(2),WO(2) | |
15 | ||
16 | C...Reset counters. Identify parton system and take copy. Check flavour. | |
17 | NSAV=N | |
18 | MSTU90=MSTU(90) | |
19 | NJET=0 | |
20 | KQSUM=0 | |
21 | DO 100 J=1,5 | |
22 | DPS(J)=0. | |
23 | 100 CONTINUE | |
24 | I=IP-1 | |
25 | 110 I=I+1 | |
26 | IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN | |
27 | CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system') | |
28 | IF(MSTU(21).GE.1) RETURN | |
29 | ENDIF | |
30 | IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110 | |
31 | KC=LUCOMP(K(I,2)) | |
32 | IF(KC.EQ.0) GOTO 110 | |
33 | KQ=KCHG(KC,2)*ISIGN(1,K(I,2)) | |
34 | IF(KQ.EQ.0) GOTO 110 | |
35 | NJET=NJET+1 | |
36 | IF(KQ.NE.2) KQSUM=KQSUM+KQ | |
37 | DO 120 J=1,5 | |
38 | K(NSAV+NJET,J)=K(I,J) | |
39 | P(NSAV+NJET,J)=P(I,J) | |
40 | DPS(J)=DPS(J)+P(I,J) | |
41 | 120 CONTINUE | |
42 | K(NSAV+NJET,3)=I | |
43 | IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND. | |
44 | &K(I+1,1).EQ.2)) GOTO 110 | |
45 | IF(NJET.NE.1.AND.KQSUM.NE.0) THEN | |
46 | CALL LUERRM(12,'(LUINDF:) unphysical flavour combination') | |
47 | IF(MSTU(21).GE.1) RETURN | |
48 | ENDIF | |
49 | ||
50 | C...Boost copied system to CM frame. Find CM energy and sum flavours. | |
51 | IF(NJET.NE.1) THEN | |
52 | MSTU(33)=1 | |
53 | CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4), | |
54 | & -DPS(2)/DPS(4),-DPS(3)/DPS(4)) | |
55 | ENDIF | |
56 | PECM=0. | |
57 | DO 130 J=1,3 | |
58 | NFI(J)=0 | |
59 | 130 CONTINUE | |
60 | DO 140 I=NSAV+1,NSAV+NJET | |
61 | PECM=PECM+P(I,4) | |
62 | KFA=IABS(K(I,2)) | |
63 | IF(KFA.LE.3) THEN | |
64 | NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2)) | |
65 | ELSEIF(KFA.GT.1000) THEN | |
66 | KFLA=MOD(KFA/1000,10) | |
67 | KFLB=MOD(KFA/100,10) | |
68 | IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2)) | |
69 | IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2)) | |
70 | ENDIF | |
71 | 140 CONTINUE | |
72 | ||
73 | C...Loop over attempts made. Reset counters. | |
74 | NTRY=0 | |
75 | 150 NTRY=NTRY+1 | |
76 | IF(NTRY.GT.200) THEN | |
77 | CALL LUERRM(14,'(LUINDF:) caught in infinite loop') | |
78 | IF(MSTU(21).GE.1) RETURN | |
79 | ENDIF | |
80 | N=NSAV+NJET | |
81 | MSTU(90)=MSTU90 | |
82 | DO 160 J=1,3 | |
83 | NFL(J)=NFI(J) | |
84 | IFET(J)=0 | |
85 | KFLF(J)=0 | |
86 | 160 CONTINUE | |
87 | ||
88 | C...Loop over jets to be fragmented. | |
89 | DO 230 IP1=NSAV+1,NSAV+NJET | |
90 | MSTJ(91)=0 | |
91 | NSAV1=N | |
92 | MSTU91=MSTU(90) | |
93 | ||
94 | C...Initial flavour and momentum values. Jet along +z axis. | |
95 | KFLH=IABS(K(IP1,2)) | |
96 | IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10) | |
97 | KFLO(2)=0 | |
98 | WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2) | |
99 | ||
100 | C...Initial values for quark or diquark jet. | |
101 | 170 IF(IABS(K(IP1,2)).NE.21) THEN | |
102 | NSTR=1 | |
103 | KFLO(1)=K(IP1,2) | |
104 | CALL LUPTDI(0,PXO(1),PYO(1)) | |
105 | WO(1)=WF | |
106 | ||
107 | C...Initial values for gluon treated like random quark jet. | |
108 | ELSEIF(MSTJ(2).LE.2) THEN | |
109 | NSTR=1 | |
110 | IF(MSTJ(2).EQ.2) MSTJ(91)=1 | |
111 | KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5) | |
112 | CALL LUPTDI(0,PXO(1),PYO(1)) | |
113 | WO(1)=WF | |
114 | ||
115 | C...Initial values for gluon treated like quark-antiquark jet pair, | |
116 | C...sharing energy according to Altarelli-Parisi splitting function. | |
117 | ELSE | |
118 | NSTR=2 | |
119 | IF(MSTJ(2).EQ.4) MSTJ(91)=1 | |
120 | KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5) | |
121 | KFLO(2)=-KFLO(1) | |
122 | CALL LUPTDI(0,PXO(1),PYO(1)) | |
123 | PXO(2)=-PXO(1) | |
124 | PYO(2)=-PYO(1) | |
125 | WO(1)=WF*RLU(0)**(1./3.) | |
126 | WO(2)=WF-WO(1) | |
127 | ENDIF | |
128 | ||
129 | C...Initial values for rank, flavour, pT and W+. | |
130 | DO 220 ISTR=1,NSTR | |
131 | 180 I=N | |
132 | MSTU(90)=MSTU91 | |
133 | IRANK=0 | |
134 | KFL1=KFLO(ISTR) | |
135 | PX1=PXO(ISTR) | |
136 | PY1=PYO(ISTR) | |
137 | W=WO(ISTR) | |
138 | ||
139 | C...New hadron. Generate flavour and hadron species. | |
140 | 190 I=I+1 | |
141 | IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN | |
142 | CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETS') | |
143 | IF(MSTU(21).GE.1) RETURN | |
144 | ENDIF | |
145 | IRANK=IRANK+1 | |
146 | K(I,1)=1 | |
147 | K(I,3)=IP1 | |
148 | K(I,4)=0 | |
149 | K(I,5)=0 | |
150 | 200 CALL LUKFDI(KFL1,0,KFL2,K(I,2)) | |
151 | IF(K(I,2).EQ.0) GOTO 180 | |
152 | IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND. | |
153 | &IABS(KFL2).GT.10) THEN | |
154 | IF(RLU(0).GT.PARJ(19)) GOTO 200 | |
155 | ENDIF | |
156 | ||
157 | C...Find hadron mass. Generate four-momentum. | |
158 | P(I,5)=ULMASS(K(I,2)) | |
159 | CALL LUPTDI(KFL1,PX2,PY2) | |
160 | P(I,1)=PX1+PX2 | |
161 | P(I,2)=PY1+PY2 | |
162 | PR=P(I,5)**2+P(I,1)**2+P(I,2)**2 | |
163 | CALL LUZDIS(KFL1,KFL2,PR,Z) | |
164 | MZSAV=0 | |
165 | IF(IABS(KFL1).GE.4.AND.IABS(KFL1).LE.8.AND.MSTU(90).LT.8) THEN | |
166 | MZSAV=1 | |
167 | MSTU(90)=MSTU(90)+1 | |
168 | MSTU(90+MSTU(90))=I | |
169 | PARU(90+MSTU(90))=Z | |
170 | ENDIF | |
171 | P(I,3)=0.5*(Z*W-PR/MAX(1E-4,Z*W)) | |
172 | P(I,4)=0.5*(Z*W+PR/MAX(1E-4,Z*W)) | |
173 | IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND. | |
174 | &P(I,3).LE.0.001) THEN | |
175 | IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180 | |
176 | P(I,3)=0.0001 | |
177 | P(I,4)=SQRT(PR) | |
178 | Z=P(I,4)/W | |
179 | ENDIF | |
180 | ||
181 | C...Remaining flavour and momentum. | |
182 | KFL1=-KFL2 | |
183 | PX1=-PX2 | |
184 | PY1=-PY2 | |
185 | W=(1.-Z)*W | |
186 | DO 210 J=1,5 | |
187 | V(I,J)=0. | |
188 | 210 CONTINUE | |
189 | ||
190 | C...Check if pL acceptable. Go back for new hadron if enough energy. | |
191 | IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) THEN | |
192 | I=I-1 | |
193 | IF(MZSAV.EQ.1) MSTU(90)=MSTU(90)-1 | |
194 | ENDIF | |
195 | IF(W.GT.PARJ(31)) GOTO 190 | |
196 | N=I | |
197 | 220 CONTINUE | |
198 | IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32) | |
199 | IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170 | |
200 | ||
201 | C...Rotate jet to new direction. | |
202 | THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2)) | |
203 | PHI=ULANGL(P(IP1,1),P(IP1,2)) | |
204 | MSTU(33)=1 | |
205 | CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0) | |
206 | K(K(IP1,3),4)=NSAV1+1 | |
207 | K(K(IP1,3),5)=N | |
208 | ||
209 | C...End of jet generation loop. Skip conservation in some cases. | |
210 | 230 CONTINUE | |
211 | IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 490 | |
212 | IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150 | |
213 | ||
214 | C...Subtract off produced hadron flavours, finished if zero. | |
215 | DO 240 I=NSAV+NJET+1,N | |
216 | KFA=IABS(K(I,2)) | |
217 | KFLA=MOD(KFA/1000,10) | |
218 | KFLB=MOD(KFA/100,10) | |
219 | KFLC=MOD(KFA/10,10) | |
220 | IF(KFLA.EQ.0) THEN | |
221 | IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB | |
222 | IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB | |
223 | ELSE | |
224 | IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2)) | |
225 | IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2)) | |
226 | IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2)) | |
227 | ENDIF | |
228 | 240 CONTINUE | |
229 | NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ | |
230 | &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3 | |
231 | IF(NREQ.EQ.0) GOTO 320 | |
232 | ||
233 | C...Take away flavour of low-momentum particles until enough freedom. | |
234 | NREM=0 | |
235 | 250 IREM=0 | |
236 | P2MIN=PECM**2 | |
237 | DO 260 I=NSAV+NJET+1,N | |
238 | P2=P(I,1)**2+P(I,2)**2+P(I,3)**2 | |
239 | IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I | |
240 | IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2 | |
241 | 260 CONTINUE | |
242 | IF(IREM.EQ.0) GOTO 150 | |
243 | K(IREM,1)=7 | |
244 | KFA=IABS(K(IREM,2)) | |
245 | KFLA=MOD(KFA/1000,10) | |
246 | KFLB=MOD(KFA/100,10) | |
247 | KFLC=MOD(KFA/10,10) | |
248 | IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8 | |
249 | IF(K(IREM,1).EQ.8) GOTO 250 | |
250 | IF(KFLA.EQ.0) THEN | |
251 | ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB | |
252 | IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN | |
253 | IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN | |
254 | ELSE | |
255 | IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2)) | |
256 | IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2)) | |
257 | IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2)) | |
258 | ENDIF | |
259 | NREM=NREM+1 | |
260 | NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ | |
261 | &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3 | |
262 | IF(NREQ.GT.NREM) GOTO 250 | |
263 | DO 270 I=NSAV+NJET+1,N | |
264 | IF(K(I,1).EQ.8) K(I,1)=1 | |
265 | 270 CONTINUE | |
266 | ||
267 | C...Find combination of existing and new flavours for hadron. | |
268 | 280 NFET=2 | |
269 | IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3 | |
270 | IF(NREQ.LT.NREM) NFET=1 | |
271 | IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0 | |
272 | DO 290 J=1,NFET | |
273 | IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0) | |
274 | KFLF(J)=ISIGN(1,NFL(1)) | |
275 | IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2)) | |
276 | IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3)) | |
277 | 290 CONTINUE | |
278 | IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0)) | |
279 | &GOTO 280 | |
280 | IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR. | |
281 | &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3) | |
282 | &.LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280 | |
283 | IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0)) | |
284 | IF(NFET.EQ.0) KFLF(2)=-KFLF(1) | |
285 | IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1)) | |
286 | IF(NFET.LE.2) KFLF(3)=0 | |
287 | IF(KFLF(3).NE.0) THEN | |
288 | KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+ | |
289 | & 100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1)) | |
290 | IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.) | |
291 | & KFLFC=KFLFC+ISIGN(2,KFLFC) | |
292 | ELSE | |
293 | KFLFC=KFLF(1) | |
294 | ENDIF | |
295 | CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF) | |
296 | IF(KF.EQ.0) GOTO 280 | |
297 | DO 300 J=1,MAX(2,NFET) | |
298 | NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J)) | |
299 | 300 CONTINUE | |
300 | ||
301 | C...Store hadron at random among free positions. | |
302 | NPOS=MIN(1+INT(RLU(0)*NREM),NREM) | |
303 | DO 310 I=NSAV+NJET+1,N | |
304 | IF(K(I,1).EQ.7) NPOS=NPOS-1 | |
305 | IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310 | |
306 | K(I,1)=1 | |
307 | K(I,2)=KF | |
308 | P(I,5)=ULMASS(K(I,2)) | |
309 | P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2) | |
310 | 310 CONTINUE | |
311 | NREM=NREM-1 | |
312 | NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+ | |
313 | &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3 | |
314 | IF(NREM.GT.0) GOTO 280 | |
315 | ||
316 | C...Compensate for missing momentum in global scheme (3 options). | |
317 | 320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN | |
318 | DO 340 J=1,3 | |
319 | PSI(J)=0. | |
320 | DO 330 I=NSAV+NJET+1,N | |
321 | PSI(J)=PSI(J)+P(I,J) | |
322 | 330 CONTINUE | |
323 | 340 CONTINUE | |
324 | PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2 | |
325 | PWS=0. | |
326 | DO 350 I=NSAV+NJET+1,N | |
327 | IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4) | |
328 | IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+ | |
329 | & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4)) | |
330 | IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1. | |
331 | 350 CONTINUE | |
332 | DO 370 I=NSAV+NJET+1,N | |
333 | IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4) | |
334 | IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+ | |
335 | & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4)) | |
336 | IF(MOD(MSTJ(3),5).EQ.3) PW=1. | |
337 | DO 360 J=1,3 | |
338 | P(I,J)=P(I,J)-PSI(J)*PW/PWS | |
339 | 360 CONTINUE | |
340 | P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2) | |
341 | 370 CONTINUE | |
342 | ||
343 | C...Compensate for missing momentum withing each jet separately. | |
344 | ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN | |
345 | DO 390 I=N+1,N+NJET | |
346 | K(I,1)=0 | |
347 | DO 380 J=1,5 | |
348 | P(I,J)=0. | |
349 | 380 CONTINUE | |
350 | 390 CONTINUE | |
351 | DO 410 I=NSAV+NJET+1,N | |
352 | IR1=K(I,3) | |
353 | IR2=N+IR1-NSAV | |
354 | K(IR2,1)=K(IR2,1)+1 | |
355 | PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/ | |
356 | & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2) | |
357 | DO 400 J=1,3 | |
358 | P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J) | |
359 | 400 CONTINUE | |
360 | P(IR2,4)=P(IR2,4)+P(I,4) | |
361 | P(IR2,5)=P(IR2,5)+PLS | |
362 | 410 CONTINUE | |
363 | PSS=0. | |
364 | DO 420 I=N+1,N+NJET | |
365 | IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2)) | |
366 | 420 CONTINUE | |
367 | DO 440 I=NSAV+NJET+1,N | |
368 | IR1=K(I,3) | |
369 | IR2=N+IR1-NSAV | |
370 | PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/ | |
371 | & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2) | |
372 | DO 430 J=1,3 | |
373 | P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS* | |
374 | & P(IR1,J) | |
375 | 430 CONTINUE | |
376 | P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2) | |
377 | 440 CONTINUE | |
378 | ENDIF | |
379 | ||
380 | C...Scale momenta for energy conservation. | |
381 | IF(MOD(MSTJ(3),5).NE.0) THEN | |
382 | PMS=0. | |
383 | PES=0. | |
384 | PQS=0. | |
385 | DO 450 I=NSAV+NJET+1,N | |
386 | PMS=PMS+P(I,5) | |
387 | PES=PES+P(I,4) | |
388 | PQS=PQS+P(I,5)**2/P(I,4) | |
389 | 450 CONTINUE | |
390 | IF(PMS.GE.PECM) GOTO 150 | |
391 | NECO=0 | |
392 | 460 NECO=NECO+1 | |
393 | PFAC=(PECM-PQS)/(PES-PQS) | |
394 | PES=0. | |
395 | PQS=0. | |
396 | DO 480 I=NSAV+NJET+1,N | |
397 | DO 470 J=1,3 | |
398 | P(I,J)=PFAC*P(I,J) | |
399 | 470 CONTINUE | |
400 | P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2) | |
401 | PES=PES+P(I,4) | |
402 | PQS=PQS+P(I,5)**2/P(I,4) | |
403 | 480 CONTINUE | |
404 | IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 460 | |
405 | ENDIF | |
406 | ||
407 | C...Origin of produced particles and parton daughter pointers. | |
408 | 490 DO 500 I=NSAV+NJET+1,N | |
409 | IF(MSTU(16).NE.2) K(I,3)=NSAV+1 | |
410 | IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3) | |
411 | 500 CONTINUE | |
412 | DO 510 I=NSAV+1,NSAV+NJET | |
413 | I1=K(I,3) | |
414 | K(I1,1)=K(I1,1)+10 | |
415 | IF(MSTU(16).NE.2) THEN | |
416 | K(I1,4)=NSAV+1 | |
417 | K(I1,5)=NSAV+1 | |
418 | ELSE | |
419 | K(I1,4)=K(I1,4)-NJET+1 | |
420 | K(I1,5)=K(I1,5)-NJET+1 | |
421 | IF(K(I1,5).LT.K(I1,4)) THEN | |
422 | K(I1,4)=0 | |
423 | K(I1,5)=0 | |
424 | ENDIF | |
425 | ENDIF | |
426 | 510 CONTINUE | |
427 | ||
428 | C...Document independent fragmentation system. Remove copy of jets. | |
429 | NSAV=NSAV+1 | |
430 | K(NSAV,1)=11 | |
431 | K(NSAV,2)=93 | |
432 | K(NSAV,3)=IP | |
433 | K(NSAV,4)=NSAV+1 | |
434 | K(NSAV,5)=N-NJET+1 | |
435 | DO 520 J=1,4 | |
436 | P(NSAV,J)=DPS(J) | |
437 | V(NSAV,J)=V(IP,J) | |
438 | 520 CONTINUE | |
439 | P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2)) | |
440 | V(NSAV,5)=0. | |
441 | DO 540 I=NSAV+NJET,N | |
442 | DO 530 J=1,5 | |
443 | K(I-NJET+1,J)=K(I,J) | |
444 | P(I-NJET+1,J)=P(I,J) | |
445 | V(I-NJET+1,J)=V(I,J) | |
446 | 530 CONTINUE | |
447 | 540 CONTINUE | |
448 | N=N-NJET+1 | |
449 | DO 550 IZ=MSTU90+1,MSTU(90) | |
450 | MSTU(90+IZ)=MSTU(90+IZ)-NJET+1 | |
451 | 550 CONTINUE | |
452 | ||
453 | C...Boost back particle system. Set production vertices. | |
454 | IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4), | |
455 | &DPS(2)/DPS(4),DPS(3)/DPS(4)) | |
456 | DO 570 I=NSAV+1,N | |
457 | DO 560 J=1,4 | |
458 | V(I,J)=V(IP,J) | |
459 | 560 CONTINUE | |
460 | 570 CONTINUE | |
461 | ||
462 | RETURN | |
463 | END |