This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lu4ent.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LU4ENT(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
5
6C...Purpose: to store four partons or particles in their CM frame, with
7C...the first along the +z axis, the last in the xz plane with x > 0
8C...and the second having y < 0 and y > 0 with equal probability.
9 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
10 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
11 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
12 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/
13
14C...Standard checks.
15 MSTU(28)=0
16 IF(MSTU(12).GE.1) CALL LULIST(0)
17 IPA=MAX(1,IABS(IP))
18 IF(IPA.GT.MSTU(4)-3) CALL LUERRM(21,
19 &'(LU4ENT:) writing outside LUJETS momory')
20 KC1=LUCOMP(KF1)
21 KC2=LUCOMP(KF2)
22 KC3=LUCOMP(KF3)
23 KC4=LUCOMP(KF4)
24 IF(KC1.EQ.0.OR.KC2.EQ.0.OR.KC3.EQ.0.OR.KC4.EQ.0) CALL LUERRM(12,
25 &'(LU4ENT:) unknown flavour code')
26
27C...Find masses. Reset K, P and V vectors.
28 PM1=0.
29 IF(MSTU(10).EQ.1) PM1=P(IPA,5)
30 IF(MSTU(10).GE.2) PM1=ULMASS(KF1)
31 PM2=0.
32 IF(MSTU(10).EQ.1) PM2=P(IPA+1,5)
33 IF(MSTU(10).GE.2) PM2=ULMASS(KF2)
34 PM3=0.
35 IF(MSTU(10).EQ.1) PM3=P(IPA+2,5)
36 IF(MSTU(10).GE.2) PM3=ULMASS(KF3)
37 PM4=0.
38 IF(MSTU(10).EQ.1) PM4=P(IPA+3,5)
39 IF(MSTU(10).GE.2) PM4=ULMASS(KF4)
40 DO 110 I=IPA,IPA+3
41 DO 100 J=1,5
42 K(I,J)=0
43 P(I,J)=0.
44 V(I,J)=0.
45 100 CONTINUE
46 110 CONTINUE
47
48C...Check flavours.
49 KQ1=KCHG(KC1,2)*ISIGN(1,KF1)
50 KQ2=KCHG(KC2,2)*ISIGN(1,KF2)
51 KQ3=KCHG(KC3,2)*ISIGN(1,KF3)
52 KQ4=KCHG(KC4,2)*ISIGN(1,KF4)
53 IF(MSTU(19).EQ.1) THEN
54 MSTU(19)=0
55 ELSEIF(KQ1.EQ.0.AND.KQ2.EQ.0.AND.KQ3.EQ.0.AND.KQ4.EQ.0) THEN
56 ELSEIF(KQ1.NE.0.AND.KQ2.EQ.2.AND.KQ3.EQ.2.AND.(KQ1+KQ4.EQ.0.OR.
57 &KQ1+KQ4.EQ.4)) THEN
58 ELSEIF(KQ1.NE.0.AND.KQ1+KQ2.EQ.0.AND.KQ3.NE.0.AND.KQ3+KQ4.EQ.0.)
59 &THEN
60 ELSE
61 CALL LUERRM(2,'(LU4ENT:) unphysical flavour combination')
62 ENDIF
63 K(IPA,2)=KF1
64 K(IPA+1,2)=KF2
65 K(IPA+2,2)=KF3
66 K(IPA+3,2)=KF4
67
68C...Store partons/particles in K vectors for normal case.
69 IF(IP.GE.0) THEN
70 K(IPA,1)=1
71 IF(KQ1.NE.0.AND.(KQ2.NE.0.OR.KQ3.NE.0.OR.KQ4.NE.0)) K(IPA,1)=2
72 K(IPA+1,1)=1
73 IF(KQ2.NE.0.AND.KQ1+KQ2.NE.0.AND.(KQ3.NE.0.OR.KQ4.NE.0))
74 & K(IPA+1,1)=2
75 K(IPA+2,1)=1
76 IF(KQ3.NE.0.AND.KQ4.NE.0) K(IPA+2,1)=2
77 K(IPA+3,1)=1
78
79C...Store partons for parton shower evolution from q-g-g-qbar or
80C...g-g-g-g event.
81 ELSEIF(KQ1+KQ2.NE.0) THEN
82 K(IPA,1)=3
83 K(IPA+1,1)=3
84 K(IPA+2,1)=3
85 K(IPA+3,1)=3
86 KCS=4
87 IF(KQ1.EQ.-1) KCS=5
88 K(IPA,KCS)=MSTU(5)*(IPA+1)
89 K(IPA,9-KCS)=MSTU(5)*(IPA+3)
90 K(IPA+1,KCS)=MSTU(5)*(IPA+2)
91 K(IPA+1,9-KCS)=MSTU(5)*IPA
92 K(IPA+2,KCS)=MSTU(5)*(IPA+3)
93 K(IPA+2,9-KCS)=MSTU(5)*(IPA+1)
94 K(IPA+3,KCS)=MSTU(5)*IPA
95 K(IPA+3,9-KCS)=MSTU(5)*(IPA+2)
96
97C...Store partons for parton shower evolution from q-qbar-q-qbar event.
98 ELSE
99 K(IPA,1)=3
100 K(IPA+1,1)=3
101 K(IPA+2,1)=3
102 K(IPA+3,1)=3
103 K(IPA,4)=MSTU(5)*(IPA+1)
104 K(IPA,5)=K(IPA,4)
105 K(IPA+1,4)=MSTU(5)*IPA
106 K(IPA+1,5)=K(IPA+1,4)
107 K(IPA+2,4)=MSTU(5)*(IPA+3)
108 K(IPA+2,5)=K(IPA+2,4)
109 K(IPA+3,4)=MSTU(5)*(IPA+2)
110 K(IPA+3,5)=K(IPA+3,4)
111 ENDIF
112
113C...Check kinematics.
114 MKERR=0
115 IF(0.5*X1*PECM.LE.PM1.OR.0.5*X2*PECM.LE.PM2.OR.0.5*(2.-X1-X2-X4)*
116 &PECM.LE.PM3.OR.0.5*X4*PECM.LE.PM4) MKERR=1
117 PA1=SQRT(MAX(1E-10,(0.5*X1*PECM)**2-PM1**2))
118 PA2=SQRT(MAX(1E-10,(0.5*X2*PECM)**2-PM2**2))
119 PA4=SQRT(MAX(1E-10,(0.5*X4*PECM)**2-PM4**2))
120 X24=X1+X2+X4-1.-X12-X14+(PM3**2-PM1**2-PM2**2-PM4**2)/PECM**2
121 CTHE4=(X1*X4-2.*X14)*PECM**2/(4.*PA1*PA4)
122 IF(ABS(CTHE4).GE.1.002) MKERR=1
123 CTHE4=MAX(-1.,MIN(1.,CTHE4))
124 STHE4=SQRT(1.-CTHE4**2)
125 CTHE2=(X1*X2-2.*X12)*PECM**2/(4.*PA1*PA2)
126 IF(ABS(CTHE2).GE.1.002) MKERR=1
127 CTHE2=MAX(-1.,MIN(1.,CTHE2))
128 STHE2=SQRT(1.-CTHE2**2)
129 CPHI2=((X2*X4-2.*X24)*PECM**2-4.*PA2*CTHE2*PA4*CTHE4)/
130 &MAX(1E-8*PECM**2,4.*PA2*STHE2*PA4*STHE4)
131 IF(ABS(CPHI2).GE.1.05) MKERR=1
132 CPHI2=MAX(-1.,MIN(1.,CPHI2))
133 IF(MKERR.EQ.1) CALL LUERRM(13,
134 &'(LU4ENT:) unphysical kinematical variable setup')
135
136C...Store partons/particles in P vectors.
137 P(IPA,3)=PA1
138 P(IPA,4)=SQRT(PA1**2+PM1**2)
139 P(IPA,5)=PM1
140 P(IPA+3,1)=PA4*STHE4
141 P(IPA+3,3)=PA4*CTHE4
142 P(IPA+3,4)=SQRT(PA4**2+PM4**2)
143 P(IPA+3,5)=PM4
144 P(IPA+1,1)=PA2*STHE2*CPHI2
145 P(IPA+1,2)=PA2*STHE2*SQRT(1.-CPHI2**2)*(-1.)**INT(RLU(0)+0.5)
146 P(IPA+1,3)=PA2*CTHE2
147 P(IPA+1,4)=SQRT(PA2**2+PM2**2)
148 P(IPA+1,5)=PM2
149 P(IPA+2,1)=-P(IPA+1,1)-P(IPA+3,1)
150 P(IPA+2,2)=-P(IPA+1,2)
151 P(IPA+2,3)=-P(IPA,3)-P(IPA+1,3)-P(IPA+3,3)
152 P(IPA+2,4)=SQRT(P(IPA+2,1)**2+P(IPA+2,2)**2+P(IPA+2,3)**2+PM3**2)
153 P(IPA+2,5)=PM3
154
155C...Set N. Optionally fragment/decay.
156 N=IPA+3
157 IF(IP.EQ.0) CALL LUEXEC
158
159 RETURN
160 END