PHOS/shaker/lu3ent.f

   1 *CMZ :          17/07/98  15.44.31  by  Federico Carminati
   2 *-- Author :
   3 C*********************************************************************
   4
   5       SUBROUTINE LU3ENT(IP,KF1,KF2,KF3,PECM,X1,X3)
   6
   7 C...Purpose: to store three partons or particles in their CM frame,
   8 C...with the first along the +z axis and the third in the (x,z)
   9 C...plane with x > 0.
  10 *KEEP,LUJETS.
  11       COMMON /LUJETS/ N,K(200000,5),P(200000,5),V(200000,5)
  12       SAVE /LUJETS/
  13 *KEEP,LUDAT1.
  14       COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  15       SAVE /LUDAT1/
  16 *KEEP,LUDAT2.
  17       COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  18       SAVE /LUDAT2/
  19 *KEND.
  20
  21 C...Standard checks.
  22       MSTU(28)=0
  23       IF(MSTU(12).GE.1) CALL LULIST(0)
  24       IPA=MAX(1,IABS(IP))
  25       IF(IPA.GT.MSTU(4)-2) CALL LUERRM(21,
  26      &'(LU3ENT:) writing outside LUJETS memory')
  27       KC1=LUCOMP(KF1)
  28       KC2=LUCOMP(KF2)
  29       KC3=LUCOMP(KF3)
  30       IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0) CALL LUERRM(12,
  31      &'(LU3ENT:) unknown flavour code')
  32
  33 C...Find masses. Reset K, P and V vectors.
  34       PM1=0.
  35       IF(MSTU(10).EQ.1) PM1=P(IPA,5)
  36       IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
  37       PM2=0.
  38       IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
  39       IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
  40       PM3=0.
  41       IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
  42       IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
  43       DO 100 I=IPA,IPA+2
  44       DO 100 J=1,5
  45       K(I,J)=0
  46       P(I,J)=0.
  47   100 V(I,J)=0.
  48
  49 C...Check flavours.
  50       KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
  51       KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
  52       KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
  53       IF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0) THEN
  54       ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.(KQ1+KQ3.EQ.0.OR.KQ1+KQ3.EQ.4))
  55      &THEN
  56       ELSE
  57         CALL LUERRM(2,'(LU3ENT:) unphysical flavour combination')
  58       ENDIF
  59       K(IPA,2)=KF1
  60       K(IPA+1,2)=KF2
  61       K(IPA+2,2)=KF3
  62
  63 C...Store partons/particles in K vectors for normal case.
  64       IF(IP.GE.0) THEN
  65         K(IPA,1)=1
  66         IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0)) K(IPA,1)=2
  67         K(IPA+1,1)=1
  68         IF(KQ2.NE.0.AND.KQ3.NE.0) K(IPA+1,1)=2
  69         K(IPA+2,1)=1
  70
  71 C...Store partons in K vectors for parton shower evolution.
  72       ELSE
  73         IF(KQ1.EQ.0.OR.KQ2.EQ.0.OR.KQ3.EQ.0) CALL LUERRM(2,
  74      &  '(LU3ENT:) requested flavours can not develop parton shower')
  75         K(IPA,1)=3
  76         K(IPA+1,1)=3
  77         K(IPA+2,1)=3
  78         KCS=4
  79         IF(KQ1.EQ.-1) KCS=5
  80         K(IPA,KCS)=MSTU(5)*(IPA+1)
  81         K(IPA,9-KCS)=MSTU(5)*(IPA+2)
  82         K(IPA+1,KCS)=MSTU(5)*(IPA+2)
  83         K(IPA+1,9-KCS)=MSTU(5)*IPA
  84         K(IPA+2,KCS)=MSTU(5)*IPA
  85         K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
  86       ENDIF
  87
  88 C...Check kinematics.
  89       MKERR=0
  90       IF(0.5*X1*PECM.LE.PM1.OR.0.5*(2.-X1-X3)*PECM.LE.PM2.OR.
  91      &0.5*X3*PECM.LE.PM3) MKERR=1
  92       PA1=SQRT(MAX(0.,(0.5*X1*PECM)**2-PM1**2))
  93       PA2=SQRT(MAX(0.,(0.5*(2.-X1-X3)*PECM)**2-PM2**2))
  94       PA3=SQRT(MAX(0.,(0.5*X3*PECM)**2-PM3**2))
  95       CTHE2=(PA3**2-PA1**2-PA2**2)/(2.*PA1*PA2)
  96       CTHE3=(PA2**2-PA1**2-PA3**2)/(2.*PA1*PA3)
  97       IF(ABS(CTHE2).GE.1.001.OR.ABS(CTHE3).GE.1.001) MKERR=1
  98       CTHE3=MAX(-1.,MIN(1.,CTHE3))
  99       IF(MKERR.NE.0) CALL LUERRM(13,
 100      &'(LU3ENT:) unphysical kinematical variable setup')
 101
 102 C...Store partons/particles in P vectors.
 103       P(IPA,3)=PA1
 104       P(IPA,4)=SQRT(PA1**2+PM1**2)
 105       P(IPA,5)=PM1
 106       P(IPA+2,1)=PA3*SQRT(1.-CTHE3**2)
 107       P(IPA+2,3)=PA3*CTHE3
 108       P(IPA+2,4)=SQRT(PA3**2+PM3**2)
 109       P(IPA+2,5)=PM3
 110       P(IPA+1,1)=-P(IPA+2,1)
 111       P(IPA+1,3)=-P(IPA,3)-P(IPA+2,3)
 112       P(IPA+1,4)=SQRT(P(IPA+1,1)**2+P(IPA+1,3)**2+PM2**2)
 113       P(IPA+1,5)=PM2
 114
 115 C...Set N. Optionally fragment/decay.
 116       N=IPA+2
 117       IF(IP.EQ.0) CALL LUEXEC
 118
 119       RETURN
 120       END