Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lu3ent.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lu3ent.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lu3ent(IP,KF1,KF2,KF3,PECM,X1,X3)
5 
6 C...Purpose: to store three partons or particles in their CM frame,
7 C...with the first along the +z axis and the third in the (x,z)
8 C...plane with x > 0.
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)-2) CALL luerrm(21,
21  &'(LU3ENT:) writing outside LUJETS memory')
22  kc1=lucomp(kf1)
23  kc2=lucomp(kf2)
24  kc3=lucomp(kf3)
25  IF(kc1.EQ.0.OR.kc2.EQ.0.OR.kc3.EQ.0) CALL luerrm(12,
26  &'(LU3ENT:) unknown flavour code')
27 
28 C...Find masses. Reset K, P and V vectors.
29  pm1=0.
30  IF(mstu(10).EQ.1) pm1=p(ipa,5)
31  IF(mstu(10).GE.2) pm1=ulmass(kf1)
32  pm2=0.
33  IF(mstu(10).EQ.1) pm2=p(ipa+1,5)
34  IF(mstu(10).GE.2) pm2=ulmass(kf2)
35  pm3=0.
36  IF(mstu(10).EQ.1) pm3=p(ipa+2,5)
37  IF(mstu(10).GE.2) pm3=ulmass(kf3)
38  DO 100 i=ipa,ipa+2
39  DO 100 j=1,5
40  k(i,j)=0
41  p(i,j)=0.
42  100 v(i,j)=0.
43 
44 C...Check flavours.
45  kq1=kchg(kc1,2)*isign(1,kf1)
46  kq2=kchg(kc2,2)*isign(1,kf2)
47  kq3=kchg(kc3,2)*isign(1,kf3)
48  IF(kq1.EQ.0.AND.kq2.EQ.0.AND.kq3.EQ.0) THEN
49  ELSEIF(kq1.NE.0.AND.kq2.EQ.2.AND.(kq1+kq3.EQ.0.OR.kq1+kq3.EQ.4))
50  &THEN
51  ELSE
52  CALL luerrm(2,'(LU3ENT:) unphysical flavour combination')
53  ENDIF
54  k(ipa,2)=kf1
55  k(ipa+1,2)=kf2
56  k(ipa+2,2)=kf3
57 
58 C...Store partons/particles in K vectors for normal case.
59  IF(ip.GE.0) THEN
60  k(ipa,1)=1
61  IF(kq1.NE.0.AND.(kq2.NE.0.OR.kq3.NE.0)) k(ipa,1)=2
62  k(ipa+1,1)=1
63  IF(kq2.NE.0.AND.kq3.NE.0) k(ipa+1,1)=2
64  k(ipa+2,1)=1
65 
66 C...Store partons in K vectors for parton shower evolution.
67  ELSE
68  IF(kq1.EQ.0.OR.kq2.EQ.0.OR.kq3.EQ.0) CALL luerrm(2,
69  & '(LU3ENT:) requested flavours can not develop parton shower')
70  k(ipa,1)=3
71  k(ipa+1,1)=3
72  k(ipa+2,1)=3
73  kcs=4
74  IF(kq1.EQ.-1) kcs=5
75  k(ipa,kcs)=mstu(5)*(ipa+1)
76  k(ipa,9-kcs)=mstu(5)*(ipa+2)
77  k(ipa+1,kcs)=mstu(5)*(ipa+2)
78  k(ipa+1,9-kcs)=mstu(5)*ipa
79  k(ipa+2,kcs)=mstu(5)*ipa
80  k(ipa+2,9-kcs)=mstu(5)*(ipa+1)
81  ENDIF
82 
83 C...Check kinematics.
84  mkerr=0
85  IF(0.5*x1*pecm.LE.pm1.OR.0.5*(2.-x1-x3)*pecm.LE.pm2.OR.
86  &0.5*x3*pecm.LE.pm3) mkerr=1
87  pa1=sqrt(max(0.,(0.5*x1*pecm)**2-pm1**2))
88  pa2=sqrt(max(0.,(0.5*(2.-x1-x3)*pecm)**2-pm2**2))
89  pa3=sqrt(max(0.,(0.5*x3*pecm)**2-pm3**2))
90  cthe2=(pa3**2-pa1**2-pa2**2)/(2.*pa1*pa2)
91  cthe3=(pa2**2-pa1**2-pa3**2)/(2.*pa1*pa3)
92  IF(abs(cthe2).GE.1.001.OR.abs(cthe3).GE.1.001) mkerr=1
93  cthe3=max(-1.,min(1.,cthe3))
94  IF(mkerr.NE.0) CALL luerrm(13,
95  &'(LU3ENT:) unphysical kinematical variable setup')
96 
97 C...Store partons/particles in P vectors.
98  p(ipa,3)=pa1
99  p(ipa,4)=sqrt(pa1**2+pm1**2)
100  p(ipa,5)=pm1
101  p(ipa+2,1)=pa3*sqrt(1.-cthe3**2)
102  p(ipa+2,3)=pa3*cthe3
103  p(ipa+2,4)=sqrt(pa3**2+pm3**2)
104  p(ipa+2,5)=pm3
105  p(ipa+1,1)=-p(ipa+2,1)
106  p(ipa+1,3)=-p(ipa,3)-p(ipa+2,3)
107  p(ipa+1,4)=sqrt(p(ipa+1,1)**2+p(ipa+1,3)**2+pm2**2)
108  p(ipa+1,5)=pm2
109 
110 C...Set N. Optionally fragment/decay.
111  n=ipa+2
112  IF(ip.EQ.0) CALL luexec
113 
114  RETURN
115  END