Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lu4ent.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lu4ent.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lu4ent(IP,KF1,KF2,KF3,KF4,PECM,X1,X2,X4,X12,X14)
5 
6 C...Purpose: to store four partons or particles in their CM frame, with
7 C...the first along the +z axis, the last in the xz plane with x > 0
8 C...and the second having y < 0 and y > 0 with equal probability.
9  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
10  SAVE /lujets/
11  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
12  SAVE /ludat1/
13  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
14  SAVE /ludat2/
15 
16 C...Standard checks.
17  mstu(28)=0
18  IF(mstu(12).GE.1) CALL lulist(0)
19  ipa=max(1,iabs(ip))
20  IF(ipa.GT.mstu(4)-3) CALL luerrm(21,
21  &'(LU4ENT:) writing outside LUJETS momory')
22  kc1=lucomp(kf1)
23  kc2=lucomp(kf2)
24  kc3=lucomp(kf3)
25  kc4=lucomp(kf4)
26  IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0.OR.kc4.EQ.0) CALL luerrm(12,
27  &'(LU4ENT:) unknown flavour code')
28 
29 C...Find masses. Reset K, P and V vectors.
30  pm1=0.
31  IF(mstu(10).EQ.1) pm1=p(ipa,5)
32  IF(mstu(10).GE.2) pm1=ulmass(kf1)
33  pm2=0.
34  IF(mstu(10).EQ.1) pm2=p(ipa+1,5)
35  IF(mstu(10).GE.2) pm2=ulmass(kf2)
36  pm3=0.
37  IF(mstu(10).EQ.1) pm3=p(ipa+2,5)
38  IF(mstu(10).GE.2) pm3=ulmass(kf3)
39  pm4=0.
40  IF(mstu(10).EQ.1) pm4=p(ipa+3,5)
41  IF(mstu(10).GE.2) pm4=ulmass(kf4)
42  DO 100 i=ipa,ipa+3
43  DO 100 j=1,5
44  k(i,j)=0
45  p(i,j)=0.
46  100 v(i,j)=0.
47 
48 C...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(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0.AND.kq4.EQ.0) THEN
54  ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.kq3.EQ.2.AND.(kq1+kq4.EQ.0.OR.
55  &kq1+kq4.EQ.4)) THEN
56  ELSEIF(kq1.NE.0.AND.kq1+kq2.EQ.0.AND.kq3.NE.0.AND.kq3+kq4.EQ.0.)
57  &THEN
58  ELSE
59  CALL luerrm(2,'(LU4ENT:) unphysical flavour combination')
60  ENDIF
61  k(ipa,2)=kf1
62  k(ipa+1,2)=kf2
63  k(ipa+2,2)=kf3
64  k(ipa+3,2)=kf4
65 
66 C...Store partons/particles in K vectors for normal case.
67  IF(ip.GE.0) THEN
68  k(ipa,1)=1
69  IF(kq1.NE.0.AND.(kq2.NE.0.OR.kq3.NE.0.OR.kq4.NE.0)) k(ipa,1)=2
70  k(ipa+1,1)=1
71  IF(kq2.NE.0.AND.kq1+kq2.NE.0.AND.(kq3.NE.0.OR.kq4.NE.0))
72  & k(ipa+1,1)=2
73  k(ipa+2,1)=1
74  IF(kq3.NE.0.AND.kq4.NE.0) k(ipa+2,1)=2
75  k(ipa+3,1)=1
76 
77 C...Store partons for parton shower evolution from q-g-g-qbar or
78 C...g-g-g-g event.
79  ELSEIF(kq1+kq2.NE.0) THEN
80  IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0.OR.kq4.EQ.0) CALL luerrm(2,
81  & '(LU4ENT:) requested flavours can not develop parton shower')
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 
97 C...Store partons for parton shower evolution from q-qbar-q-qbar event.
98  ELSE
99  IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0.OR.kq4.EQ.0) CALL luerrm(2,
100  & '(LU4ENT:) requested flavours can not develop parton shower')
101  k(ipa,1)=3
102  k(ipa+1,1)=3
103  k(ipa+2,1)=3
104  k(ipa+3,1)=3
105  k(ipa,4)=mstu(5)*(ipa+1)
106  k(ipa,5)=k(ipa,4)
107  k(ipa+1,4)=mstu(5)*ipa
108  k(ipa+1,5)=k(ipa+1,4)
109  k(ipa+2,4)=mstu(5)*(ipa+3)
110  k(ipa+2,5)=k(ipa+2,4)
111  k(ipa+3,4)=mstu(5)*(ipa+2)
112  k(ipa+3,5)=k(ipa+3,4)
113  ENDIF
114 
115 C...Check kinematics.
116  mkerr=0
117  IF(0.5*x1*pecm.LE.pm1.OR.0.5*x2*pecm.LE.pm2.OR.0.5*(2.-x1-x2-x4)*
118  &pecm.LE.pm3.OR.0.5*x4*pecm.LE.pm4) mkerr=1
119  pa1=sqrt(max(0.,(0.5*x1*pecm)**2-pm1**2))
120  pa2=sqrt(max(0.,(0.5*x2*pecm)**2-pm2**2))
121  pa3=sqrt(max(0.,(0.5*(2.-x1-x2-x4)*pecm)**2-pm3**2))
122  pa4=sqrt(max(0.,(0.5*x4*pecm)**2-pm4**2))
123  x24=x1+x2+x4-1.-x12-x14+(pm3**2-pm1**2-pm2**2-pm4**2)/pecm**2
124  cthe4=(x1*x4-2.*x14)*pecm**2/(4.*pa1*pa4)
125  IF(abs(cthe4).GE.1.002) mkerr=1
126  cthe4=max(-1.,min(1.,cthe4))
127  sthe4=sqrt(1.-cthe4**2)
128  cthe2=(x1*x2-2.*x12)*pecm**2/(4.*pa1*pa2)
129  IF(abs(cthe2).GE.1.002) mkerr=1
130  cthe2=max(-1.,min(1.,cthe2))
131  sthe2=sqrt(1.-cthe2**2)
132  cphi2=((x2*x4-2.*x24)*pecm**2-4.*pa2*cthe2*pa4*cthe4)/
133  &(4.*pa2*sthe2*pa4*sthe4)
134  IF(abs(cphi2).GE.1.05) mkerr=1
135  cphi2=max(-1.,min(1.,cphi2))
136  IF(mkerr.EQ.1) CALL luerrm(13,
137  &'(LU4ENT:) unphysical kinematical variable setup')
138 
139 C...Store partons/particles in P vectors.
140  p(ipa,3)=pa1
141  p(ipa,4)=sqrt(pa1**2+pm1**2)
142  p(ipa,5)=pm1
143  p(ipa+3,1)=pa4*sthe4
144  p(ipa+3,3)=pa4*cthe4
145  p(ipa+3,4)=sqrt(pa4**2+pm4**2)
146  p(ipa+3,5)=pm4
147  p(ipa+1,1)=pa2*sthe2*cphi2
148  p(ipa+1,2)=pa2*sthe2*sqrt(1.-cphi2**2)*(-1.)**int(rlu(0)+0.5)
149  p(ipa+1,3)=pa2*cthe2
150  p(ipa+1,4)=sqrt(pa2**2+pm2**2)
151  p(ipa+1,5)=pm2
152  p(ipa+2,1)=-p(ipa+1,1)-p(ipa+3,1)
153  p(ipa+2,2)=-p(ipa+1,2)
154  p(ipa+2,3)=-p(ipa,3)-p(ipa+1,3)-p(ipa+3,3)
155  p(ipa+2,4)=sqrt(p(ipa+2,1)**2+p(ipa+2,2)**2+p(ipa+2,3)**2+pm3**2)
156  p(ipa+2,5)=pm3
157 
158 C...Set N. Optionally fragment/decay.
159  n=ipa+3
160  IF(ip.EQ.0) CALL luexec
161 
162  RETURN
163  END