hijing/src/luindf.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUINDF(IP)
0005
0006 C...Purpose: to handle the fragmentation of a jet system (or a single
0007 C...jet) according to independent fragmentation models.
0008       IMPLICIT DOUBLE PRECISION(D)
0009       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0010       SAVE /LUJETS/
0011       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0012       SAVE /LUDAT1/
0013       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0014       SAVE /LUDAT2/
0015       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
0016      &KFLO(2),PXO(2),PYO(2),WO(2)
0017
0018 C...Reset counters. Identify parton system and take copy. Check flavour.
0019       NSAV=N
0020       NJET=0
0021       KQSUM=0
0022       DO 100 J=1,5
0023   100 DPS(J)=0.
0024       I=IP-1
0025   110 I=I+1
0026       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
0027         CALL LUERRM(12,'(LUINDF:) failed to reconstruct jet system')
0028         IF(MSTU(21).GE.1) RETURN
0029       ENDIF
0030       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
0031       KC=LUCOMP(K(I,2))
0032       IF(KC.EQ.0) GOTO 110
0033       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
0034       IF(KQ.EQ.0) GOTO 110
0035       NJET=NJET+1
0036       IF(KQ.NE.2) KQSUM=KQSUM+KQ
0037       DO 120 J=1,5
0038       K(NSAV+NJET,J)=K(I,J)
0039       P(NSAV+NJET,J)=P(I,J)
0040   120 DPS(J)=DPS(J)+P(I,J)
0041       K(NSAV+NJET,3)=I
0042       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
0043      &K(I+1,1).EQ.2)) GOTO 110
0044       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
0045         CALL LUERRM(12,'(LUINDF:) unphysical flavour combination')
0046         IF(MSTU(21).GE.1) RETURN
0047       ENDIF
0048
0049 C...Boost copied system to CM frame. Find CM energy and sum flavours.
0050       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)/DPS(4),
0051      &-DPS(2)/DPS(4),-DPS(3)/DPS(4))
0052       PECM=0.
0053       DO 130 J=1,3
0054   130 NFI(J)=0
0055       DO 140 I=NSAV+1,NSAV+NJET
0056       PECM=PECM+P(I,4)
0057       KFA=IABS(K(I,2))
0058       IF(KFA.LE.3) THEN
0059         NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
0060       ELSEIF(KFA.GT.1000) THEN
0061         KFLA=MOD(KFA/1000,10)
0062         KFLB=MOD(KFA/100,10)
0063         IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
0064         IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
0065       ENDIF
0066   140 CONTINUE
0067
0068 C...Loop over attempts made. Reset counters.
0069       NTRY=0
0070   150 NTRY=NTRY+1
0071       N=NSAV+NJET
0072       IF(NTRY.GT.200) THEN
0073         CALL LUERRM(14,'(LUINDF:) caught in infinite loop')
0074         IF(MSTU(21).GE.1) RETURN
0075       ENDIF
0076       DO 160 J=1,3
0077       NFL(J)=NFI(J)
0078       IFET(J)=0
0079   160 KFLF(J)=0
0080
0081 C...Loop over jets to be fragmented.
0082       DO 230 IP1=NSAV+1,NSAV+NJET
0083       MSTJ(91)=0
0084       NSAV1=N
0085
0086 C...Initial flavour and momentum values. Jet along +z axis.
0087       KFLH=IABS(K(IP1,2))
0088       IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
0089       KFLO(2)=0
0090       WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
0091
0092 C...Initial values for quark or diquark jet.
0093   170 IF(IABS(K(IP1,2)).NE.21) THEN
0094         NSTR=1
0095         KFLO(1)=K(IP1,2)
0096         CALL LUPTDI(0,PXO(1),PYO(1))
0097         WO(1)=WF
0098
0099 C...Initial values for gluon treated like random quark jet.
0100       ELSEIF(MSTJ(2).LE.2) THEN
0101         NSTR=1
0102         IF(MSTJ(2).EQ.2) MSTJ(91)=1
0103         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
0104         CALL LUPTDI(0,PXO(1),PYO(1))
0105         WO(1)=WF
0106
0107 C...Initial values for gluon treated like quark-antiquark jet pair,
0108 C...sharing energy according to Altarelli-Parisi splitting function.
0109       ELSE
0110         NSTR=2
0111         IF(MSTJ(2).EQ.4) MSTJ(91)=1
0112         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU(0))*(-1)**INT(RLU(0)+0.5)
0113         KFLO(2)=-KFLO(1)
0114         CALL LUPTDI(0,PXO(1),PYO(1))
0115         PXO(2)=-PXO(1)
0116         PYO(2)=-PYO(1)
0117         WO(1)=WF*RLU(0)**(1./3.)
0118         WO(2)=WF-WO(1)
0119       ENDIF
0120
0121 C...Initial values for rank, flavour, pT and W+.
0122       DO 220 ISTR=1,NSTR
0123   180 I=N
0124       IRANK=0
0125       KFL1=KFLO(ISTR)
0126       PX1=PXO(ISTR)
0127       PY1=PYO(ISTR)
0128       W=WO(ISTR)
0129
0130 C...New hadron. Generate flavour and hadron species.
0131   190 I=I+1
0132       IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
0133         CALL LUERRM(11,'(LUINDF:) no more memory left in LUJETS')
0134         IF(MSTU(21).GE.1) RETURN
0135       ENDIF
0136       IRANK=IRANK+1
0137       K(I,1)=1
0138       K(I,3)=IP1
0139       K(I,4)=0
0140       K(I,5)=0
0141   200 CALL LUKFDI(KFL1,0,KFL2,K(I,2))
0142       IF(K(I,2).EQ.0) GOTO 180
0143       IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
0144      &IABS(KFL2).GT.10) THEN
0145         IF(RLU(0).GT.PARJ(19)) GOTO 200
0146       ENDIF
0147
0148 C...Find hadron mass. Generate four-momentum.
0149       P(I,5)=ULMASS(K(I,2))
0150       CALL LUPTDI(KFL1,PX2,PY2)
0151       P(I,1)=PX1+PX2
0152       P(I,2)=PY1+PY2
0153       PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
0154       CALL LUZDIS(KFL1,KFL2,PR,Z)
0155       P(I,3)=0.5*(Z*W-PR/(Z*W))
0156       P(I,4)=0.5*(Z*W+PR/(Z*W))
0157       IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
0158      &P(I,3).LE.0.001) THEN
0159         IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
0160         P(I,3)=0.0001
0161         P(I,4)=SQRT(PR)
0162         Z=P(I,4)/W
0163       ENDIF
0164
0165 C...Remaining flavour and momentum.
0166       KFL1=-KFL2
0167       PX1=-PX2
0168       PY1=-PY2
0169       W=(1.-Z)*W
0170       DO 210 J=1,5
0171   210 V(I,J)=0.
0172
0173 C...Check if pL acceptable. Go back for new hadron if enough energy.
0174       IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) I=I-1
0175       IF(W.GT.PARJ(31)) GOTO 190
0176   220 N=I
0177       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
0178       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
0179
0180 C...Rotate jet to new direction.
0181       THE=ULANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
0182       PHI=ULANGL(P(IP1,1),P(IP1,2))
0183       CALL LUDBRB(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
0184       K(K(IP1,3),4)=NSAV1+1
0185       K(K(IP1,3),5)=N
0186
0187 C...End of jet generation loop. Skip conservation in some cases.
0188   230 CONTINUE
0189       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470
0190       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
0191
0192 C...Subtract off produced hadron flavours, finished if zero.
0193       DO 240 I=NSAV+NJET+1,N
0194       KFA=IABS(K(I,2))
0195       KFLA=MOD(KFA/1000,10)
0196       KFLB=MOD(KFA/100,10)
0197       KFLC=MOD(KFA/10,10)
0198       IF(KFLA.EQ.0) THEN
0199         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
0200         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
0201       ELSE
0202         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
0203         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
0204         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
0205       ENDIF
0206   240 CONTINUE
0207       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0208      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0209       IF(NREQ.EQ.0) GOTO 320
0210
0211 C...Take away flavour of low-momentum particles until enough freedom.
0212       NREM=0
0213   250 IREM=0
0214       P2MIN=PECM**2
0215       DO 260 I=NSAV+NJET+1,N
0216       P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
0217       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
0218   260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
0219       IF(IREM.EQ.0) GOTO 150
0220       K(IREM,1)=7
0221       KFA=IABS(K(IREM,2))
0222       KFLA=MOD(KFA/1000,10)
0223       KFLB=MOD(KFA/100,10)
0224       KFLC=MOD(KFA/10,10)
0225       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
0226       IF(K(IREM,1).EQ.8) GOTO 250
0227       IF(KFLA.EQ.0) THEN
0228         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
0229         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
0230         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
0231       ELSE
0232         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
0233         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
0234         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
0235       ENDIF
0236       NREM=NREM+1
0237       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0238      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0239       IF(NREQ.GT.NREM) GOTO 250
0240       DO 270 I=NSAV+NJET+1,N
0241   270 IF(K(I,1).EQ.8) K(I,1)=1
0242
0243 C...Find combination of existing and new flavours for hadron.
0244   280 NFET=2
0245       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
0246       IF(NREQ.LT.NREM) NFET=1
0247       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
0248       DO 290 J=1,NFET
0249       IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU(0)
0250       KFLF(J)=ISIGN(1,NFL(1))
0251       IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
0252   290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
0253       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
0254      &GOTO 280
0255       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
0256      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).
0257      &LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
0258       IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU(0))
0259       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
0260       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU(0)),-KFLF(1))
0261       IF(NFET.LE.2) KFLF(3)=0
0262       IF(KFLF(3).NE.0) THEN
0263         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
0264      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
0265         IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU(0).GT.1.)
0266      &  KFLFC=KFLFC+ISIGN(2,KFLFC)
0267       ELSE
0268         KFLFC=KFLF(1)
0269       ENDIF
0270       CALL LUKFDI(KFLFC,KFLF(2),KFLDMP,KF)
0271       IF(KF.EQ.0) GOTO 280
0272       DO 300 J=1,MAX(2,NFET)
0273   300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
0274
0275 C...Store hadron at random among free positions.
0276       NPOS=MIN(1+INT(RLU(0)*NREM),NREM)
0277       DO 310 I=NSAV+NJET+1,N
0278       IF(K(I,1).EQ.7) NPOS=NPOS-1
0279       IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
0280       K(I,1)=1
0281       K(I,2)=KF
0282       P(I,5)=ULMASS(K(I,2))
0283       P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0284   310 CONTINUE
0285       NREM=NREM-1
0286       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
0287      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
0288       IF(NREM.GT.0) GOTO 280
0289
0290 C...Compensate for missing momentum in global scheme (3 options).
0291   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
0292         DO 330 J=1,3
0293         PSI(J)=0.
0294         DO 330 I=NSAV+NJET+1,N
0295   330   PSI(J)=PSI(J)+P(I,J)
0296         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
0297         PWS=0.
0298         DO 340 I=NSAV+NJET+1,N
0299         IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
0300         IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
0301      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
0302   340   IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
0303         DO 360 I=NSAV+NJET+1,N
0304         IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
0305         IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
0306      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
0307         IF(MOD(MSTJ(3),5).EQ.3) PW=1.
0308         DO 350 J=1,3
0309   350   P(I,J)=P(I,J)-PSI(J)*PW/PWS
0310   360   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0311
0312 C...Compensate for missing momentum withing each jet separately.
0313       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
0314         DO 370 I=N+1,N+NJET
0315         K(I,1)=0
0316         DO 370 J=1,5
0317   370   P(I,J)=0.
0318         DO 390 I=NSAV+NJET+1,N
0319         IR1=K(I,3)
0320         IR2=N+IR1-NSAV
0321         K(IR2,1)=K(IR2,1)+1
0322         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
0323      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
0324         DO 380 J=1,3
0325   380   P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
0326         P(IR2,4)=P(IR2,4)+P(I,4)
0327   390   P(IR2,5)=P(IR2,5)+PLS
0328         PSS=0.
0329         DO 400 I=N+1,N+NJET
0330   400   IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
0331         DO 420 I=NSAV+NJET+1,N
0332         IR1=K(I,3)
0333         IR2=N+IR1-NSAV
0334         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
0335      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
0336         DO 410 J=1,3
0337   410   P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
0338      &  P(IR1,J)
0339   420   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0340       ENDIF
0341
0342 C...Scale momenta for energy conservation.
0343       IF(MOD(MSTJ(3),5).NE.0) THEN
0344         PMS=0.
0345         PES=0.
0346         PQS=0.
0347         DO 430 I=NSAV+NJET+1,N
0348         PMS=PMS+P(I,5)
0349         PES=PES+P(I,4)
0350   430   PQS=PQS+P(I,5)**2/P(I,4)
0351         IF(PMS.GE.PECM) GOTO 150
0352         NECO=0
0353   440   NECO=NECO+1
0354         PFAC=(PECM-PQS)/(PES-PQS)
0355         PES=0.
0356         PQS=0.
0357         DO 460 I=NSAV+NJET+1,N
0358         DO 450 J=1,3
0359   450   P(I,J)=PFAC*P(I,J)
0360         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
0361         PES=PES+P(I,4)
0362   460   PQS=PQS+P(I,5)**2/P(I,4)
0363         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440
0364       ENDIF
0365
0366 C...Origin of produced particles and parton daughter pointers.
0367   470 DO 480 I=NSAV+NJET+1,N
0368       IF(MSTU(16).NE.2) K(I,3)=NSAV+1
0369   480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
0370       DO 490 I=NSAV+1,NSAV+NJET
0371       I1=K(I,3)
0372       K(I1,1)=K(I1,1)+10
0373       IF(MSTU(16).NE.2) THEN
0374         K(I1,4)=NSAV+1
0375         K(I1,5)=NSAV+1
0376       ELSE
0377         K(I1,4)=K(I1,4)-NJET+1
0378         K(I1,5)=K(I1,5)-NJET+1
0379         IF(K(I1,5).LT.K(I1,4)) THEN
0380           K(I1,4)=0
0381           K(I1,5)=0
0382         ENDIF
0383       ENDIF
0384   490 CONTINUE
0385
0386 C...Document independent fragmentation system. Remove copy of jets.
0387       NSAV=NSAV+1
0388       K(NSAV,1)=11
0389       K(NSAV,2)=93
0390       K(NSAV,3)=IP
0391       K(NSAV,4)=NSAV+1
0392       K(NSAV,5)=N-NJET+1
0393       DO 500 J=1,4
0394       P(NSAV,J)=DPS(J)
0395   500 V(NSAV,J)=V(IP,J)
0396       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
0397       V(NSAV,5)=0.
0398       DO 510 I=NSAV+NJET,N
0399       DO 510 J=1,5
0400       K(I-NJET+1,J)=K(I,J)
0401       P(I-NJET+1,J)=P(I,J)
0402   510 V(I-NJET+1,J)=V(I,J)
0403       N=N-NJET+1
0404
0405 C...Boost back particle system. Set production vertices.
0406       IF(NJET.NE.1) CALL LUDBRB(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
0407      &DPS(2)/DPS(4),DPS(3)/DPS(4))
0408       DO 520 I=NSAV+1,N
0409       DO 520 J=1,4
0410   520 V(I,J)=V(IP,J)
0411
0412       RETURN
0413       END