hijing/src/lu3ent.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LU3ENT(IP,KF1,KF2,KF3,PECM,X1,X3)
0005
0006 C...Purpose: to store three partons or particles in their CM frame,
0007 C...with the first along the +z axis and the third in the (x,z)
0008 C...plane with x > 0.
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)-2) CALL LUERRM(21,
0021      &'(LU3ENT:) writing outside LUJETS memory')
0022       KC1=LUCOMP(KF1)
0023       KC2=LUCOMP(KF2)
0024       KC3=LUCOMP(KF3)
0025       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0) CALL LUERRM(12,
0026      &'(LU3ENT:) unknown flavour code')
0027
0028 C...Find masses. Reset K, P and V vectors.
0029       PM1=0.
0030       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
0031       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
0032       PM2=0.
0033       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
0034       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
0035       PM3=0.
0036       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
0037       IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
0038       DO 100 I=IPA,IPA+2
0039       DO 100 J=1,5
0040       K(I,J)=0
0041       P(I,J)=0.
0042   100 V(I,J)=0.
0043
0044 C...Check flavours.
0045       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
0046       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
0047       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
0048       IF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0) THEN
0049       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.(KQ1+KQ3.EQ.0.OR.KQ1+KQ3.EQ.4))
0050      &THEN
0051       ELSE
0052         CALL LUERRM(2,'(LU3ENT:) unphysical flavour combination')
0053       ENDIF
0054       K(IPA,2)=KF1
0055       K(IPA+1,2)=KF2
0056       K(IPA+2,2)=KF3
0057
0058 C...Store partons/particles in K vectors for normal case.
0059       IF(IP.GE.0) THEN
0060         K(IPA,1)=1
0061         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0)) K(IPA,1)=2
0062         K(IPA+1,1)=1
0063         IF(KQ2.NE.0.AND.KQ3.NE.0) K(IPA+1,1)=2
0064         K(IPA+2,1)=1
0065
0066 C...Store partons in K vectors for parton shower evolution.
0067       ELSE
0068         IF(KQ1.EQ.0.OR.KQ2.EQ.0.OR.KQ3.EQ.0) CALL LUERRM(2,
0069      &  '(LU3ENT:) requested flavours can not develop parton shower')
0070         K(IPA,1)=3
0071         K(IPA+1,1)=3
0072         K(IPA+2,1)=3
0073         KCS=4
0074         IF(KQ1.EQ.-1) KCS=5
0075         K(IPA,KCS)=MSTU(5)*(IPA+1)
0076         K(IPA,9-KCS)=MSTU(5)*(IPA+2)
0077         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
0078         K(IPA+1,9-KCS)=MSTU(5)*IPA
0079         K(IPA+2,KCS)=MSTU(5)*IPA
0080         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
0081       ENDIF
0082
0083 C...Check kinematics.
0084       MKERR=0
0085       IF(0.5*X1*PECM.LE.PM1.OR.0.5*(2.-X1-X3)*PECM.LE.PM2.OR.
0086      &0.5*X3*PECM.LE.PM3) MKERR=1
0087       PA1=SQRT(MAX(0.,(0.5*X1*PECM)**2-PM1**2))
0088       PA2=SQRT(MAX(0.,(0.5*(2.-X1-X3)*PECM)**2-PM2**2))
0089       PA3=SQRT(MAX(0.,(0.5*X3*PECM)**2-PM3**2))
0090       CTHE2=(PA3**2-PA1**2-PA2**2)/(2.*PA1*PA2)
0091       CTHE3=(PA2**2-PA1**2-PA3**2)/(2.*PA1*PA3)
0092       IF(ABS(CTHE2).GE.1.001.OR.ABS(CTHE3).GE.1.001) MKERR=1
0093       CTHE3=MAX(-1.,MIN(1.,CTHE3))
0094       IF(MKERR.NE.0) CALL LUERRM(13,
0095      &'(LU3ENT:) unphysical kinematical variable setup')
0096
0097 C...Store partons/particles in P vectors.
0098       P(IPA,3)=PA1
0099       P(IPA,4)=SQRT(PA1**2+PM1**2)
0100       P(IPA,5)=PM1
0101       P(IPA+2,1)=PA3*SQRT(1.-CTHE3**2)
0102       P(IPA+2,3)=PA3*CTHE3
0103       P(IPA+2,4)=SQRT(PA3**2+PM3**2)
0104       P(IPA+2,5)=PM3
0105       P(IPA+1,1)=-P(IPA+2,1)
0106       P(IPA+1,3)=-P(IPA,3)-P(IPA+2,3)
0107       P(IPA+1,4)=SQRT(P(IPA+1,1)**2+P(IPA+1,3)**2+PM2**2)
0108       P(IPA+1,5)=PM2
0109
0110 C...Set N. Optionally fragment/decay.
0111       N=IPA+2
0112       IF(IP.EQ.0) CALL LUEXEC
0113
0114       RETURN
0115       END