Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lux4jt.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lux4jt.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lux4jt(NJET,CUT,KFL,ECM,KFLN,X1,X2,X4,X12,X14)
5 
6 C...Purpose: to select the kinematical variables of four-jet events.
7  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
8  SAVE /ludat1/
9  dimension wta(4),wtb(4),wtc(4),wtd(4),wte(4)
10 
11 C...Common constants. Colour factors for QCD and Abelian gluon theory.
12  pmq=ulmass(kfl)
13  qme=(2.*pmq/ecm)**2
14  ct=log(1./cut-5.)
15  IF(mstj(109).EQ.0) THEN
16  cf=4./3.
17  cn=3.
18  tr=2.
19  ELSE
20  cf=1.
21  cn=0.
22  tr=12.
23  ENDIF
24 
25 C...Choice of process (qqbargg or qqbarqqbar).
26  100 njet=4
27  it=1
28  IF(parj(155).GT.rlu(0)) it=2
29  IF(mstj(101).LE.-3) it=-mstj(101)-2
30  IF(it.EQ.1) wtmx=0.7/cut**2
31  IF(it.EQ.1.AND.mstj(109).EQ.2) wtmx=0.6/cut**2
32  IF(it.EQ.2) wtmx=0.1125*cf*tr/cut**2
33  id=1
34 
35 C...Sample the five kinematical variables (for qqgg preweighted in y34).
36  110 y134=3.*cut+(1.-6.*cut)*rlu(0)
37  y234=3.*cut+(1.-6.*cut)*rlu(0)
38  IF(it.EQ.1) y34=(1.-5.*cut)*exp(-ct*rlu(0))
39  IF(it.EQ.2) y34=cut+(1.-6.*cut)*rlu(0)
40  IF(y34.LE.y134+y234-1..OR.y34.GE.y134*y234) goto 110
41  vt=rlu(0)
42  cp=cos(paru(1)*rlu(0))
43  y14=(y134-y34)*vt
44  y13=y134-y14-y34
45  vb=y34*(1.-y134-y234+y34)/((y134-y34)*(y234-y34))
46  y24=0.5*(y234-y34)*(1.-4.*sqrt(max(0.,vt*(1.-vt)*vb*(1.-vb)))*
47  &cp-(1.-2.*vt)*(1.-2.*vb))
48  y23=y234-y34-y24
49  y12=1.-y134-y23-y24
50  IF(min(y12,y13,y14,y23,y24).LE.cut) goto 110
51  y123=y12+y13+y23
52  y124=y12+y14+y24
53 
54 C...Calculate matrix elements for qqgg or qqqq process.
55  ic=0
56  wttot=0.
57  120 ic=ic+1
58  IF(it.EQ.1) THEN
59  wta(ic)=(y12*y34**2-y13*y24*y34+y14*y23*y34+3.*y12*y23*y34+
60  & 3.*y12*y14*y34+4.*y12**2*y34-y13*y23*y24+2.*y12*y23*y24-
61  & y13*y14*y24-2.*y12*y13*y24+2.*y12**2*y24+y14*y23**2+2.*y12*
62  & y23**2+y14**2*y23+4.*y12*y14*y23+4.*y12**2*y23+2.*y12*y14**2+
63  & 2.*y12*y13*y14+4.*y12**2*y14+2.*y12**2*y13+2.*y12**3)/(2.*y13*
64  & y134*y234*y24)+(y24*y34+y12*y34+y13*y24-y14*y23+y12*y13)/(y13*
65  & y134**2)+2.*y23*(1.-y13)/(y13*y134*y24)+y34/(2.*y13*y24)
66  wtb(ic)=(y12*y24*y34+y12*y14*y34-y13*y24**2+y13*y14*y24+2.*y12*
67  & y14*y24)/(y13*y134*y23*y14)+y12*(1.+y34)*y124/(y134*y234*y14*
68  & y24)-(2.*y13*y24+y14**2+y13*y23+2.*y12*y13)/(y13*y134*y14)+
69  & y12*y123*y124/(2.*y13*y14*y23*y24)
70  wtc(ic)=-(5.*y12*y34**2+2.*y12*y24*y34+2.*y12*y23*y34+2.*y12*
71  & y14*y34+2.*y12*y13*y34+4.*y12**2*y34-y13*y24**2+y14*y23*y24+
72  & y13*y23*y24+y13*y14*y24-y12*y14*y24-y13**2*y24-3.*y12*y13*y24-
73  & y14*y23**2-y14**2*y23+y13*y14*y23-3.*y12*y14*y23-y12*y13*y23)/
74  & (4.*y134*y234*y34**2)+(3.*y12*y34**2-3.*y13*y24*y34+3.*y12*y24*
75  & y34+3.*y14*y23*y34-y13*y24**2-y12*y23*y34+6.*y12*y14*y34+2.*y12*
76  & y13*y34-2.*y12**2*y34+y14*y23*y24-3.*y13*y23*y24-2.*y13*y14*
77  & y24+4.*y12*y14*y24+2.*y12*y13*y24+3.*y14*y23**2+2.*y14**2*y23+
78  & 2.*y14**2*y12+2.*y12**2*y14+6.*y12*y14*y23-2.*y12*y13**2-
79  & 2.*y12**2*y13)/(4.*y13*y134*y234*y34)
80  wtc(ic)=wtc(ic)+(2.*y12*y34**2-2.*y13*y24*y34+y12*y24*y34+
81  & 4.*y13*y23*y34+4.*y12*y14*y34+2.*y12*y13*y34+2.*y12**2*y34-
82  & y13*y24**2+3.*y14*y23*y24+4.*y13*y23*y24-2.*y13*y14*y24+
83  & 4.*y12*y14*y24+2.*y12*y13*y24+2.*y14*y23**2+4.*y13*y23**2+
84  & 2.*y13*y14*y23+2.*y12*y14*y23+4.*y12*y13*y23+2.*y12*y14**2+4.*
85  & y12**2*y13+4.*y12*y13*y14+2.*y12**2*y14)/(4.*y13*y134*y24*y34)-
86  & (y12*y34**2-2.*y14*y24*y34-2.*y13*y24*y34-y14*y23*y34+y13*y23*
87  & y34+y12*y14*y34+2.*y12*y13*y34-2.*y14**2*y24-4.*y13*y14*y24-
88  & 4.*y13**2*y24-y14**2*y23-y13**2*y23+y12*y13*y14-y12*y13**2)/
89  & (2.*y13*y34*y134**2)+(y12*y34**2-4.*y14*y24*y34-2.*y13*y24*y34-
90  & 2.*y14*y23*y34-4.*y13*y23*y34-4.*y12*y14*y34-4.*y12*y13*y34-
91  & 2.*y13*y14*y24+2.*y13**2*y24+2.*y14**2*y23-2.*y13*y14*y23-
92  & y12*y14**2-6.*y12*y13*y14-y12*y13**2)/(4.*y34**2*y134**2)
93  wttot=wttot+y34*cf*(cf*wta(ic)+(cf-0.5*cn)*wtb(ic)+cn*wtc(ic))/
94  & 8.
95  ELSE
96  wtd(ic)=(y13*y23*y34+y12*y23*y34-y12**2*y34+y13*y23*y24+2.*y12*
97  & y23*y24-y14*y23**2+y12*y13*y24+y12*y14*y23+y12*y13*y14)/(y13**2*
98  & y123**2)-(y12*y34**2-y13*y24*y34+y12*y24*y34-y14*y23*y34-y12*
99  & y23*y34-y13*y24**2+y14*y23*y24-y13*y23*y24-y13**2*y24+y14*
100  & y23**2)/(y13**2*y123*y134)+(y13*y14*y12+y34*y14*y12-y34**2*y12+
101  & y13*y14*y24+2.*y34*y14*y24-y23*y14**2+y34*y13*y24+y34*y23*y14+
102  & y34*y13*y23)/(y13**2*y134**2)-(y34*y12**2-y13*y24*y12+y34*y24*
103  & y12-y23*y14*y12-y34*y14*y12-y13*y24**2+y23*y14*y24-y13*y14*y24-
104  & y13**2*y24+y23*y14**2)/(y13**2*y134*y123)
105  wte(ic)=(y12*y34*(y23-y24+y14+y13)+y13*y24**2-y14*y23*y24+y13*
106  & y23*y24+y13*y14*y24+y13**2*y24-y14*y23*(y14+y23+y13))/(y13*y23*
107  & y123*y134)-y12*(y12*y34-y23*y24-y13*y24-y14*y23-y14*y13)/(y13*
108  & y23*y123**2)-(y14+y13)*(y24+y23)*y34/(y13*y23*y134*y234)+
109  & (y12*y34*(y14-y24+y23+y13)+y13*y24**2-y23*y14*y24+y13*y14*y24+
110  & y13*y23*y24+y13**2*y24-y23*y14*(y14+y23+y13))/(y13*y14*y134*
111  & y123)-y34*(y34*y12-y14*y24-y13*y24-y23*y14-y23*y13)/(y13*y14*
112  & y134**2)-(y23+y13)*(y24+y14)*y12/(y13*y14*y123*y124)
113  wttot=wttot+cf*(tr*wtd(ic)+(cf-0.5*cn)*wte(ic))/16.
114  ENDIF
115 
116 C...Permutations of momenta in matrix element. Weighting.
117  130 IF(ic.EQ.1.OR.ic.EQ.3.OR.id.EQ.2.OR.id.EQ.3) THEN
118  ysav=y13
119  y13=y14
120  y14=ysav
121  ysav=y23
122  y23=y24
123  y24=ysav
124  ysav=y123
125  y123=y124
126  y124=ysav
127  ENDIF
128  IF(ic.EQ.2.OR.ic.EQ.4.OR.id.EQ.3.OR.id.EQ.4) THEN
129  ysav=y13
130  y13=y23
131  y23=ysav
132  ysav=y14
133  y14=y24
134  y24=ysav
135  ysav=y134
136  y134=y234
137  y234=ysav
138  ENDIF
139  IF(ic.LE.3) goto 120
140  IF(id.EQ.1.AND.wttot.LT.rlu(0)*wtmx) goto 110
141  ic=5
142 
143 C...qqgg events: string configuration and event type.
144  IF(it.EQ.1) THEN
145  IF(mstj(109).EQ.0.AND.id.EQ.1) THEN
146  parj(156)=y34*(2.*(wta(1)+wta(2)+wta(3)+wta(4))+4.*(wtc(1)+
147  & wtc(2)+wtc(3)+wtc(4)))/(9.*wttot)
148  IF(wta(2)+wta(4)+2.*(wtc(2)+wtc(4)).GT.rlu(0)*(wta(1)+wta(2)+
149  & wta(3)+wta(4)+2.*(wtc(1)+wtc(2)+wtc(3)+wtc(4)))) id=2
150  IF(id.EQ.2) goto 130
151  ELSEIF(mstj(109).EQ.2.AND.id.EQ.1) THEN
152  parj(156)=y34*(wta(1)+wta(2)+wta(3)+wta(4))/(8.*wttot)
153  IF(wta(2)+wta(4).GT.rlu(0)*(wta(1)+wta(2)+wta(3)+wta(4))) id=2
154  IF(id.EQ.2) goto 130
155  ENDIF
156  mstj(120)=3
157  IF(mstj(109).EQ.0.AND.0.5*y34*(wtc(1)+wtc(2)+wtc(3)+wtc(4)).GT.
158  & rlu(0)*wttot) mstj(120)=4
159  kfln=21
160 
161 C...Mass cuts. Kinematical variables out.
162  IF(y12.LE.cut+qme) njet=2
163  IF(njet.EQ.2) goto 150
164  q12=0.5*(1.-sqrt(1.-qme/y12))
165  x1=1.-(1.-q12)*y234-q12*y134
166  x4=1.-(1.-q12)*y134-q12*y234
167  x2=1.-y124
168  x12=(1.-q12)*y13+q12*y23
169  x14=y12-0.5*qme
170  IF(y134*y234/((1.-x1)*(1.-x4)).LE.rlu(0)) njet=2
171 
172 C...qqbarqqbar events: string configuration, choose new flavour.
173  ELSE
174  IF(id.EQ.1) THEN
175  wtr=rlu(0)*(wtd(1)+wtd(2)+wtd(3)+wtd(4))
176  IF(wtr.LT.wtd(2)+wtd(3)+wtd(4)) id=2
177  IF(wtr.LT.wtd(3)+wtd(4)) id=3
178  IF(wtr.LT.wtd(4)) id=4
179  IF(id.GE.2) goto 130
180  ENDIF
181  mstj(120)=5
182  parj(156)=cf*tr*(wtd(1)+wtd(2)+wtd(3)+wtd(4))/(16.*wttot)
183  140 kfln=1+int(5.*rlu(0))
184  IF(kfln.NE.kfl.AND.0.2*parj(156).LE.rlu(0)) goto 140
185  IF(kfln.EQ.kfl.AND.1.-0.8*parj(156).LE.rlu(0)) goto 140
186  IF(kfln.GT.mstj(104)) njet=2
187  pmqn=ulmass(kfln)
188  qmen=(2.*pmqn/ecm)**2
189 
190 C...Mass cuts. Kinematical variables out.
191  IF(y24.LE.cut+qme.OR.y13.LE.1.1*qmen) njet=2
192  IF(njet.EQ.2) goto 150
193  q24=0.5*(1.-sqrt(1.-qme/y24))
194  q13=0.5*(1.-sqrt(1.-qmen/y13))
195  x1=1.-(1.-q24)*y123-q24*y134
196  x4=1.-(1.-q24)*y134-q24*y123
197  x2=1.-(1.-q13)*y234-q13*y124
198  x12=(1.-q24)*((1.-q13)*y14+q13*y34)+q24*((1.-q13)*y12+q13*y23)
199  x14=y24-0.5*qme
200  x34=(1.-q24)*((1.-q13)*y23+q13*y12)+q24*((1.-q13)*y34+q13*y14)
201  IF(pmq**2+pmqn**2+min(x12,x34)*ecm**2.LE.
202  & (parj(127)+pmq+pmqn)**2) njet=2
203  IF(y123*y134/((1.-x1)*(1.-x4)).LE.rlu(0)) njet=2
204  ENDIF
205  150 IF(mstj(101).LE.-2.AND.njet.EQ.2) goto 100
206 
207  RETURN
208  END