pythiaeRHIC/pythia/py4ent.f

0001
0002 C*********************************************************************
0003
0004 C...PY4ENT
0005 C...Stores four partons or particles in their CM frame, with
0006 C...the first along the +z axis, the last in the xz plane with x > 0
0007 C...and the second having y < 0 and y > 0 with equal probability.
0008
0009       SUBROUTINE PY4ENT(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
0010
0011 C...Double precision and integer declarations.
0012       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0013       IMPLICIT INTEGER(I-N)
0014       INTEGER PYK,PYCHGE,PYCOMP
0015 C...Commonblocks.
0016       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0017       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0018       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0019       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
0020
0021 C...Standard checks.
0022       MSTU(28)=0
0023       IF(MSTU(12).NE.12345) CALL PYLIST(0)
0024       IPA=MAX(1,IABS(IP))
0025       IF(IPA.GT.MSTU(4)-3) CALL PYERRM(21,
0026      &'(PY4ENT:) writing outside PYJETS momory')
0027       KC1=PYCOMP(KF1)
0028       KC2=PYCOMP(KF2)
0029       KC3=PYCOMP(KF3)
0030       KC4=PYCOMP(KF4)
0031       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0.OR.KC4.EQ.0) CALL PYERRM(12,
0032      &'(PY4ENT:) unknown flavour code')
0033
0034 C...Find masses. Reset K, P and V vectors.
0035       PM1=0D0
0036       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
0037       IF(MSTU(10).GE.2) PM1=PYMASS(KF1)
0038       PM2=0D0
0039       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
0040       IF(MSTU(10).GE.2) PM2=PYMASS(KF2)
0041       PM3=0D0
0042       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
0043       IF(MSTU(10).GE.2) PM3=PYMASS(KF3)
0044       PM4=0D0
0045       IF(MSTU(10).EQ.1) PM4=P(IPA+3,5)
0046       IF(MSTU(10).GE.2) PM4=PYMASS(KF4)
0047       DO 110 I=IPA,IPA+3
0048         DO 100 J=1,5
0049           K(I,J)=0
0050           P(I,J)=0D0
0051           V(I,J)=0D0
0052   100   CONTINUE
0053   110 CONTINUE
0054
0055 C...Check flavours.
0056       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
0057       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
0058       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
0059       KQ4=KCHG(KC4,2)*ISIGN(1,KF4)
0060       IF(MSTU(19).EQ.1) THEN
0061         MSTU(19)=0
0062       ELSEIF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0.AND.KQ4.EQ.0) THEN
0063       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.KQ3.EQ.2.AND.(KQ1+KQ4.EQ.0.OR.
0064      &  KQ1+KQ4.EQ.4)) THEN
0065       ELSEIF(KQ1.NE.0.AND.KQ1+KQ2.EQ.0.AND.KQ3.NE.0.AND.KQ3+KQ4.EQ.0D0)
0066      &  THEN
0067       ELSE
0068         CALL PYERRM(2,'(PY4ENT:) unphysical flavour combination')
0069       ENDIF
0070       K(IPA,2)=KF1
0071       K(IPA+1,2)=KF2
0072       K(IPA+2,2)=KF3
0073       K(IPA+3,2)=KF4
0074
0075 C...Store partons/particles in K vectors for normal case.
0076       IF(IP.GE.0) THEN
0077         K(IPA,1)=1
0078         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0.OR.KQ4.NE.0)) K(IPA,1)=2
0079         K(IPA+1,1)=1
0080         IF(KQ2.NE.0.AND.KQ1+KQ2.NE.0.AND.(KQ3.NE.0.OR.KQ4.NE.0))
0081      &  K(IPA+1,1)=2
0082         K(IPA+2,1)=1
0083         IF(KQ3.NE.0.AND.KQ4.NE.0) K(IPA+2,1)=2
0084         K(IPA+3,1)=1
0085
0086 C...Store partons for parton shower evolution from q-g-g-qbar or
0087 C...g-g-g-g event.
0088       ELSEIF(KQ1+KQ2.NE.0) THEN
0089         K(IPA,1)=3
0090         K(IPA+1,1)=3
0091         K(IPA+2,1)=3
0092         K(IPA+3,1)=3
0093         KCS=4
0094         IF(KQ1.EQ.-1) KCS=5
0095         K(IPA,KCS)=MSTU(5)*(IPA+1)
0096         K(IPA,9-KCS)=MSTU(5)*(IPA+3)
0097         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
0098         K(IPA+1,9-KCS)=MSTU(5)*IPA
0099         K(IPA+2,KCS)=MSTU(5)*(IPA+3)
0100         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
0101         K(IPA+3,KCS)=MSTU(5)*IPA
0102         K(IPA+3,9-KCS)=MSTU(5)*(IPA+2)
0103
0104 C...Store partons for parton shower evolution from q-qbar-q-qbar event.
0105       ELSE
0106         K(IPA,1)=3
0107         K(IPA+1,1)=3
0108         K(IPA+2,1)=3
0109         K(IPA+3,1)=3
0110         K(IPA,4)=MSTU(5)*(IPA+1)
0111         K(IPA,5)=K(IPA,4)
0112         K(IPA+1,4)=MSTU(5)*IPA
0113         K(IPA+1,5)=K(IPA+1,4)
0114         K(IPA+2,4)=MSTU(5)*(IPA+3)
0115         K(IPA+2,5)=K(IPA+2,4)
0116         K(IPA+3,4)=MSTU(5)*(IPA+2)
0117         K(IPA+3,5)=K(IPA+3,4)
0118       ENDIF
0119
0120 C...Check kinematics.
0121       MKERR=0
0122       IF(0.5D0*X1*PECM.LE.PM1.OR.0.5D0*X2*PECM.LE.PM2.OR.
0123      &0.5D0*(2D0-X1-X2-X4)*PECM.LE.PM3.OR.0.5D0*X4*PECM.LE.PM4)
0124      &MKERR=1
0125       PA1=SQRT(MAX(1D-10,(0.5D0*X1*PECM)**2-PM1**2))
0126       PA2=SQRT(MAX(1D-10,(0.5D0*X2*PECM)**2-PM2**2))
0127       PA4=SQRT(MAX(1D-10,(0.5D0*X4*PECM)**2-PM4**2))
0128       X24=X1+X2+X4-1D0-X12-X14+(PM3**2-PM1**2-PM2**2-PM4**2)/PECM**2
0129       CTHE4=(X1*X4-2D0*X14)*PECM**2/(4D0*PA1*PA4)
0130       IF(ABS(CTHE4).GE.1.002D0) MKERR=1
0131       CTHE4=MAX(-1D0,MIN(1D0,CTHE4))
0132       STHE4=SQRT(1D0-CTHE4**2)
0133       CTHE2=(X1*X2-2D0*X12)*PECM**2/(4D0*PA1*PA2)
0134       IF(ABS(CTHE2).GE.1.002D0) MKERR=1
0135       CTHE2=MAX(-1D0,MIN(1D0,CTHE2))
0136       STHE2=SQRT(1D0-CTHE2**2)
0137       CPHI2=((X2*X4-2D0*X24)*PECM**2-4D0*PA2*CTHE2*PA4*CTHE4)/
0138      &MAX(1D-8*PECM**2,4D0*PA2*STHE2*PA4*STHE4)
0139       IF(ABS(CPHI2).GE.1.05D0) MKERR=1
0140       CPHI2=MAX(-1D0,MIN(1D0,CPHI2))
0141       IF(MKERR.EQ.1) CALL PYERRM(13,
0142      &'(PY4ENT:) unphysical kinematical variable setup')
0143
0144 C...Store partons/particles in P vectors.
0145       P(IPA,3)=PA1
0146       P(IPA,4)=SQRT(PA1**2+PM1**2)
0147       P(IPA,5)=PM1
0148       P(IPA+3,1)=PA4*STHE4
0149       P(IPA+3,3)=PA4*CTHE4
0150       P(IPA+3,4)=SQRT(PA4**2+PM4**2)
0151       P(IPA+3,5)=PM4
0152       P(IPA+1,1)=PA2*STHE2*CPHI2
0153       P(IPA+1,2)=PA2*STHE2*SQRT(1D0-CPHI2**2)*(-1D0)**INT(PYR(0)+0.5D0)
0154       P(IPA+1,3)=PA2*CTHE2
0155       P(IPA+1,4)=SQRT(PA2**2+PM2**2)
0156       P(IPA+1,5)=PM2
0157       P(IPA+2,1)=-P(IPA+1,1)-P(IPA+3,1)
0158       P(IPA+2,2)=-P(IPA+1,2)
0159       P(IPA+2,3)=-P(IPA,3)-P(IPA+1,3)-P(IPA+3,3)
0160       P(IPA+2,4)=SQRT(P(IPA+2,1)**2+P(IPA+2,2)**2+P(IPA+2,3)**2+PM3**2)
0161       P(IPA+2,5)=PM3
0162
0163 C...Set N. Optionally fragment/decay.
0164       N=IPA+3
0165       IF(IP.EQ.0) CALL PYEXEC
0166
0167       RETURN
0168       END