pythiaeRHIC/pythia/py3ent.f

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