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