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