hijing/src/lu4ent.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LU4ENT(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
0005
0006 C...Purpose: to store four partons or particles in their CM frame, with
0007 C...the first along the +z axis, the last in the xz plane with x > 0
0008 C...and the second having y < 0 and y > 0 with equal probability.
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
0016 C...Standard checks.
0017       MSTU(28)=0
0018       IF(MSTU(12).GE.1) CALL LULIST(0)
0019       IPA=MAX(1,IABS(IP))
0020       IF(IPA.GT.MSTU(4)-3) CALL LUERRM(21,
0021      &'(LU4ENT:) writing outside LUJETS momory')
0022       KC1=LUCOMP(KF1)
0023       KC2=LUCOMP(KF2)
0024       KC3=LUCOMP(KF3)
0025       KC4=LUCOMP(KF4)
0026       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0.OR.KC4.EQ.0) CALL LUERRM(12,
0027      &'(LU4ENT:) unknown flavour code')
0028
0029 C...Find masses. Reset K, P and V vectors.
0030       PM1=0.
0031       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
0032       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
0033       PM2=0.
0034       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
0035       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
0036       PM3=0.
0037       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
0038       IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
0039       PM4=0.
0040       IF(MSTU(10).EQ.1) PM4=P(IPA+3,5)
0041       IF(MSTU(10).GE.2) PM4=ULMASS(KF4)
0042       DO 100 I=IPA,IPA+3
0043       DO 100 J=1,5
0044       K(I,J)=0
0045       P(I,J)=0.
0046   100 V(I,J)=0.
0047
0048 C...Check flavours.
0049       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
0050       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
0051       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
0052       KQ4=KCHG(KC4,2)*ISIGN(1,KF4)
0053       IF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0.AND.KQ4.EQ.0) THEN
0054       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.KQ3.EQ.2.AND.(KQ1+KQ4.EQ.0.OR.
0055      &KQ1+KQ4.EQ.4)) THEN
0056       ELSEIF(KQ1.NE.0.AND.KQ1+KQ2.EQ.0.AND.KQ3.NE.0.AND.KQ3+KQ4.EQ.0.)
0057      &THEN
0058       ELSE
0059         CALL LUERRM(2,'(LU4ENT:) unphysical flavour combination')
0060       ENDIF
0061       K(IPA,2)=KF1
0062       K(IPA+1,2)=KF2
0063       K(IPA+2,2)=KF3
0064       K(IPA+3,2)=KF4
0065
0066 C...Store partons/particles in K vectors for normal case.
0067       IF(IP.GE.0) THEN
0068         K(IPA,1)=1
0069         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0.OR.KQ4.NE.0)) K(IPA,1)=2
0070         K(IPA+1,1)=1
0071         IF(KQ2.NE.0.AND.KQ1+KQ2.NE.0.AND.(KQ3.NE.0.OR.KQ4.NE.0))
0072      &  K(IPA+1,1)=2
0073         K(IPA+2,1)=1
0074         IF(KQ3.NE.0.AND.KQ4.NE.0) K(IPA+2,1)=2
0075         K(IPA+3,1)=1
0076
0077 C...Store partons for parton shower evolution from q-g-g-qbar or
0078 C...g-g-g-g event.
0079       ELSEIF(KQ1+KQ2.NE.0) THEN
0080         IF(KQ1.EQ.0.OR.KQ2.EQ.0.OR.KQ3.EQ.0.OR.KQ4.EQ.0) CALL LUERRM(2,
0081      &  '(LU4ENT:) requested flavours can not develop parton shower')
0082         K(IPA,1)=3
0083         K(IPA+1,1)=3
0084         K(IPA+2,1)=3
0085         K(IPA+3,1)=3
0086         KCS=4
0087         IF(KQ1.EQ.-1) KCS=5
0088         K(IPA,KCS)=MSTU(5)*(IPA+1)
0089         K(IPA,9-KCS)=MSTU(5)*(IPA+3)
0090         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
0091         K(IPA+1,9-KCS)=MSTU(5)*IPA
0092         K(IPA+2,KCS)=MSTU(5)*(IPA+3)
0093         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
0094         K(IPA+3,KCS)=MSTU(5)*IPA
0095         K(IPA+3,9-KCS)=MSTU(5)*(IPA+2)
0096
0097 C...Store partons for parton shower evolution from q-qbar-q-qbar event.
0098       ELSE
0099         IF(KQ1.EQ.0.OR.KQ2.EQ.0.OR.KQ3.EQ.0.OR.KQ4.EQ.0) CALL LUERRM(2,
0100      &  '(LU4ENT:) requested flavours can not develop parton shower')
0101         K(IPA,1)=3
0102         K(IPA+1,1)=3
0103         K(IPA+2,1)=3
0104         K(IPA+3,1)=3
0105         K(IPA,4)=MSTU(5)*(IPA+1)
0106         K(IPA,5)=K(IPA,4)
0107         K(IPA+1,4)=MSTU(5)*IPA
0108         K(IPA+1,5)=K(IPA+1,4)
0109         K(IPA+2,4)=MSTU(5)*(IPA+3)
0110         K(IPA+2,5)=K(IPA+2,4)
0111         K(IPA+3,4)=MSTU(5)*(IPA+2)
0112         K(IPA+3,5)=K(IPA+3,4)
0113       ENDIF
0114
0115 C...Check kinematics.
0116       MKERR=0
0117       IF(0.5*X1*PECM.LE.PM1.OR.0.5*X2*PECM.LE.PM2.OR.0.5*(2.-X1-X2-X4)*
0118      &PECM.LE.PM3.OR.0.5*X4*PECM.LE.PM4) MKERR=1
0119       PA1=SQRT(MAX(0.,(0.5*X1*PECM)**2-PM1**2))
0120       PA2=SQRT(MAX(0.,(0.5*X2*PECM)**2-PM2**2))
0121       PA3=SQRT(MAX(0.,(0.5*(2.-X1-X2-X4)*PECM)**2-PM3**2))
0122       PA4=SQRT(MAX(0.,(0.5*X4*PECM)**2-PM4**2))
0123       X24=X1+X2+X4-1.-X12-X14+(PM3**2-PM1**2-PM2**2-PM4**2)/PECM**2
0124       CTHE4=(X1*X4-2.*X14)*PECM**2/(4.*PA1*PA4)
0125       IF(ABS(CTHE4).GE.1.002) MKERR=1
0126       CTHE4=MAX(-1.,MIN(1.,CTHE4))
0127       STHE4=SQRT(1.-CTHE4**2)
0128       CTHE2=(X1*X2-2.*X12)*PECM**2/(4.*PA1*PA2)
0129       IF(ABS(CTHE2).GE.1.002) MKERR=1
0130       CTHE2=MAX(-1.,MIN(1.,CTHE2))
0131       STHE2=SQRT(1.-CTHE2**2)
0132       CPHI2=((X2*X4-2.*X24)*PECM**2-4.*PA2*CTHE2*PA4*CTHE4)/
0133      &(4.*PA2*STHE2*PA4*STHE4)
0134       IF(ABS(CPHI2).GE.1.05) MKERR=1
0135       CPHI2=MAX(-1.,MIN(1.,CPHI2))
0136       IF(MKERR.EQ.1) CALL LUERRM(13,
0137      &'(LU4ENT:) unphysical kinematical variable setup')
0138
0139 C...Store partons/particles in P vectors.
0140       P(IPA,3)=PA1
0141       P(IPA,4)=SQRT(PA1**2+PM1**2)
0142       P(IPA,5)=PM1
0143       P(IPA+3,1)=PA4*STHE4
0144       P(IPA+3,3)=PA4*CTHE4
0145       P(IPA+3,4)=SQRT(PA4**2+PM4**2)
0146       P(IPA+3,5)=PM4
0147       P(IPA+1,1)=PA2*STHE2*CPHI2
0148       P(IPA+1,2)=PA2*STHE2*SQRT(1.-CPHI2**2)*(-1.)**INT(RLU(0)+0.5)
0149       P(IPA+1,3)=PA2*CTHE2
0150       P(IPA+1,4)=SQRT(PA2**2+PM2**2)
0151       P(IPA+1,5)=PM2
0152       P(IPA+2,1)=-P(IPA+1,1)-P(IPA+3,1)
0153       P(IPA+2,2)=-P(IPA+1,2)
0154       P(IPA+2,3)=-P(IPA,3)-P(IPA+1,3)-P(IPA+3,3)
0155       P(IPA+2,4)=SQRT(P(IPA+2,1)**2+P(IPA+2,2)**2+P(IPA+2,3)**2+PM3**2)
0156       P(IPA+2,5)=PM3
0157
0158 C...Set N. Optionally fragment/decay.
0159       N=IPA+3
0160       IF(IP.EQ.0) CALL LUEXEC
0161
0162       RETURN
0163       END