1 *CMZ : 17/07/98 15.44.32 by Federico Carminati
3 C*********************************************************************
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)
11 COMMON /LUJETS/ N,K(200000,5),P(200000,5),V(200000,5)
14 COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
17 COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
20 DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
21 &KFLO(2),PXO(2),PYO(2),WO(2)
23 C...Reset counters. Identify parton system and take copy. Check flavour.
31 IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
32 CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system')
33 IF(MSTU(21).GE.1) RETURN
35 IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
38 KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
41 IF(KQ.NE.2) KQSUM=KQSUM+KQ
45 120 DPS(J)=DPS(J)+P(I,J)
47 IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
48 &K(I+1,1).EQ.2)) GOTO 110
49 IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
50 CALL LUERRM(12,'(LUINDF:) unphysical flavour combination')
51 IF(MSTU(21).GE.1) RETURN
54 C...Boost copied system to CM frame. Find CM energy and sum flavours.
57 CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4),
58 & -DPS(2)/DPS(4),-DPS(3)/DPS(4))
63 DO 140 I=NSAV+1,NSAV+NJET
67 NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
68 ELSEIF(KFA.GT.1000) THEN
71 IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
72 IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
76 C...Loop over attempts made. Reset counters.
81 CALL LUERRM(14,'(LUINDF:) caught in infinite loop')
82 IF(MSTU(21).GE.1) RETURN
89 C...Loop over jets to be fragmented.
90 DO 230 IP1=NSAV+1,NSAV+NJET
94 C...Initial flavour and momentum values. Jet along +z axis.
96 IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
98 WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
100 C...Initial values for quark or diquark jet.
101 170 IF(IABS(K(IP1,2)).NE.21) THEN
104 CALL LUPTDI(0,PXO(1),PYO(1))
107 C...Initial values for gluon treated like random quark jet.
108 ELSEIF(MSTJ(2).LE.2) THEN
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))
115 C...Initial values for gluon treated like quark-antiquark jet pair,
116 C...sharing energy according to Altarelli-Parisi splitting function.
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)
122 CALL LUPTDI(0,PXO(1),PYO(1))
125 WO(1)=WF*RLU(0)**(1./3.)
129 C...Initial values for rank, flavour, pT and W+.
138 C...New hadron. Generate flavour and hadron species.
140 IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
141 CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETS')
142 IF(MSTU(21).GE.1) RETURN
149 200 CALL LUKFDI(KFL1,0,KFL2,K(I,2))
150 IF(K(I,2).EQ.0) GOTO 180
151 IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
152 &IABS(KFL2).GT.10) THEN
153 IF(RLU(0).GT.PARJ(19)) GOTO 200
156 C...Find hadron mass. Generate four-momentum.
157 P(I,5)=ULMASS(K(I,2))
158 CALL LUPTDI(KFL1,PX2,PY2)
161 PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
162 CALL LUZDIS(KFL1,KFL2,PR,Z)
163 P(I,3)=0.5*(Z*W-PR/(Z*W))
164 P(I,4)=0.5*(Z*W+PR/(Z*W))
165 IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
166 &P(I,3).LE.0.001) THEN
167 IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
173 C...Remaining flavour and momentum.
181 C...Check if pL acceptable. Go back for new hadron if enough energy.
182 IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) I=I-1
183 IF(W.GT.PARJ(31)) GOTO 190
185 IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
186 IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
188 C...Rotate jet to new direction.
189 THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
190 PHI=ULANGL(P(IP1,1),P(IP1,2))
192 CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
193 K(K(IP1,3),4)=NSAV1+1
196 C...End of jet generation loop. Skip conservation in some cases.
198 IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470
199 IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
201 C...Subtract off produced hadron flavours, finished if zero.
202 DO 240 I=NSAV+NJET+1,N
204 KFLA=MOD(KFA/1000,10)
208 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
209 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
211 IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
212 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
213 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
216 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
217 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
218 IF(NREQ.EQ.0) GOTO 320
220 C...Take away flavour of low-momentum particles until enough freedom.
224 DO 260 I=NSAV+NJET+1,N
225 P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
226 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
227 260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
228 IF(IREM.EQ.0) GOTO 150
231 KFLA=MOD(KFA/1000,10)
234 IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
235 IF(K(IREM,1).EQ.8) GOTO 250
237 ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
238 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
239 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
241 IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
242 IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
243 IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
246 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
247 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
248 IF(NREQ.GT.NREM) GOTO 250
249 DO 270 I=NSAV+NJET+1,N
250 270 IF(K(I,1).EQ.8) K(I,1)=1
252 C...Find combination of existing and new flavours for hadron.
254 IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
255 IF(NREQ.LT.NREM) NFET=1
256 IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
258 IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0)
259 KFLF(J)=ISIGN(1,NFL(1))
260 IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
261 290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
262 IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
264 IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
265 &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).
266 <.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
267 IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0))
268 IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
269 IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1))
270 IF(NFET.LE.2) KFLF(3)=0
271 IF(KFLF(3).NE.0) THEN
272 KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
273 & 100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
274 IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.)
275 & KFLFC=KFLFC+ISIGN(2,KFLFC)
279 CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF)
281 DO 300 J=1,MAX(2,NFET)
282 300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
284 C...Store hadron at random among free positions.
285 NPOS=MIN(1+INT(RLU(0)*NREM),NREM)
286 DO 310 I=NSAV+NJET+1,N
287 IF(K(I,1).EQ.7) NPOS=NPOS-1
288 IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
291 P(I,5)=ULMASS(K(I,2))
292 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
295 NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
296 &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
297 IF(NREM.GT.0) GOTO 280
299 C...Compensate for missing momentum in global scheme (3 options).
300 320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
303 DO 330 I=NSAV+NJET+1,N
304 330 PSI(J)=PSI(J)+P(I,J)
305 PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
307 DO 340 I=NSAV+NJET+1,N
308 IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
309 IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
310 & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
311 340 IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
312 DO 360 I=NSAV+NJET+1,N
313 IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
314 IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
315 & PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
316 IF(MOD(MSTJ(3),5).EQ.3) PW=1.
318 350 P(I,J)=P(I,J)-PSI(J)*PW/PWS
319 360 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
321 C...Compensate for missing momentum withing each jet separately.
322 ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
327 DO 390 I=NSAV+NJET+1,N
331 PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
332 & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
334 380 P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
335 P(IR2,4)=P(IR2,4)+P(I,4)
336 390 P(IR2,5)=P(IR2,5)+PLS
339 400 IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
340 DO 420 I=NSAV+NJET+1,N
343 PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
344 & (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
346 410 P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
348 420 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
351 C...Scale momenta for energy conservation.
352 IF(MOD(MSTJ(3),5).NE.0) THEN
356 DO 430 I=NSAV+NJET+1,N
359 430 PQS=PQS+P(I,5)**2/P(I,4)
360 IF(PMS.GE.PECM) GOTO 150
363 PFAC=(PECM-PQS)/(PES-PQS)
366 DO 460 I=NSAV+NJET+1,N
368 450 P(I,J)=PFAC*P(I,J)
369 P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
371 460 PQS=PQS+P(I,5)**2/P(I,4)
372 IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440
375 C...Origin of produced particles and parton daughter pointers.
376 470 DO 480 I=NSAV+NJET+1,N
377 IF(MSTU(16).NE.2) K(I,3)=NSAV+1
378 480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
379 DO 490 I=NSAV+1,NSAV+NJET
382 IF(MSTU(16).NE.2) THEN
386 K(I1,4)=K(I1,4)-NJET+1
387 K(I1,5)=K(I1,5)-NJET+1
388 IF(K(I1,5).LT.K(I1,4)) THEN
395 C...Document independent fragmentation system. Remove copy of jets.
404 500 V(NSAV,J)=V(IP,J)
405 P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
411 510 V(I-NJET+1,J)=V(I,J)
414 C...Boost back particle system. Set production vertices.
415 IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
416 &DPS(2)/DPS(4),DPS(3)/DPS(4))