Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyx3jt.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyx3jt.f
1 
2 C*********************************************************************
3 
4 C...PYX3JT
5 C...Selects the kinematical variables of three-jet events.
6 
7  SUBROUTINE pyx3jt(NJET,CUT,KFL,ECM,X1,X2)
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 array.
17  dimension zhup(5,12)
18 
19 C...Coefficients of Zhu second order parametrization.
20  DATA ((zhup(ic1,ic2),ic2=1,12),ic1=1,5)/
21  &18.29d0, 89.56d0, 4.541d0, -52.09d0, -109.8d0, 24.90d0,
22  &11.63d0, 3.683d0, 17.50d0,0.002440d0, -1.362d0,-0.3537d0,
23  &11.42d0, 6.299d0, -22.55d0, -8.915d0, 59.25d0, -5.855d0,
24  &-32.85d0, -1.054d0, -16.90d0,0.006489d0,-0.8156d0,0.01095d0,
25  &7.847d0, -3.964d0, -35.83d0, 1.178d0, 29.39d0, 0.2806d0,
26  &47.82d0, -12.36d0, -56.72d0, 0.04054d0,-0.4365d0, 0.6062d0,
27  &5.441d0, -56.89d0, -50.27d0, 15.13d0, 114.3d0, -18.19d0,
28  &97.05d0, -1.890d0, -139.9d0, 0.08153d0,-0.4984d0, 0.9439d0,
29  &-17.65d0, 51.44d0, -58.32d0, 70.95d0, -255.7d0, -78.99d0,
30  &476.9d0, 29.65d0, -239.3d0, 0.4745d0, -1.174d0, 6.081d0/
31 
32 C...Dilogarithm of x for x<0.5 (x>0.5 obtained by analytic trick).
33  dilog(x)=x+x**2/4d0+x**3/9d0+x**4/16d0+x**5/25d0+x**6/36d0+
34  &x**7/49d0
35 
36 C...Event type. Mass effect factors and other common constants.
37  mstj(120)=2
38  mstj(121)=0
39  pmq=pymass(kfl)
40  qme=(2d0*pmq/ecm)**2
41  IF(mstj(109).NE.1) THEN
42  cutl=log(cut)
43  cutd=log(1d0/cut-2d0)
44  IF(mstj(109).EQ.0) THEN
45  cf=4d0/3d0
46  cn=3d0
47  tr=2d0
48  wtmx=min(20d0,37d0-6d0*cutd)
49  IF(mstj(110).EQ.2) wtmx=2d0*(7.5d0+80d0*cut)
50  ELSE
51  cf=1d0
52  cn=0d0
53  tr=12d0
54  wtmx=0d0
55  ENDIF
56 
57 C...Alpha_strong and effects of optimized Q^2 scale. Maximum weight.
58  als2pi=paru(118)/paru(2)
59  wtopt=0d0
60  IF(mstj(111).EQ.1) wtopt=(33d0-2d0*mstu(112))/6d0*
61  & log(parj(169))*als2pi
62  wtmax=max(0d0,1d0+wtopt+als2pi*wtmx)
63 
64 C...Choose three-jet events in allowed region.
65  100 njet=3
66  110 y13l=cutl+cutd*pyr(0)
67  y23l=cutl+cutd*pyr(0)
68  y13=exp(y13l)
69  y23=exp(y23l)
70  y12=1d0-y13-y23
71  IF(y12.LE.cut) goto 110
72  IF(y13**2+y23**2+2d0*y12.LE.2d0*pyr(0)) goto 110
73 
74 C...Second order corrections.
75  IF(mstj(101).EQ.2.AND.mstj(110).LE.1) THEN
76  y12l=log(y12)
77  y13m=log(1d0-y13)
78  y23m=log(1d0-y23)
79  y12m=log(1d0-y12)
80  IF(y13.LE.0.5d0) y13i=dilog(y13)
81  IF(y13.GE.0.5d0) y13i=1.644934d0-y13l*y13m-dilog(1d0-y13)
82  IF(y23.LE.0.5d0) y23i=dilog(y23)
83  IF(y23.GE.0.5d0) y23i=1.644934d0-y23l*y23m-dilog(1d0-y23)
84  IF(y12.LE.0.5d0) y12i=dilog(y12)
85  IF(y12.GE.0.5d0) y12i=1.644934d0-y12l*y12m-dilog(1d0-y12)
86  wt1=(y13**2+y23**2+2d0*y12)/(y13*y23)
87  wt2=cf*(-2d0*(cutl-y12l)**2-3d0*cutl-1d0+3.289868d0+
88  & 2d0*(2d0*cutl-y12l)*cut/y12)+
89  & cn*((cutl-y12l)**2-(cutl-y13l)**2-(cutl-y23l)**2-
90  & 11d0*cutl/6d0+67d0/18d0+1.644934d0-(2d0*cutl-y12l)*cut/y12+
91  & (2d0*cutl-y13l)*cut/y13+(2d0*cutl-y23l)*cut/y23)+
92  & tr*(2d0*cutl/3d0-10d0/9d0)+
93  & cf*(y12/(y12+y13)+y12/(y12+y23)+(y12+y23)/y13+(y12+y13)/y23+
94  & y13l*(4d0*y12**2+2d0*y12*y13+4d0*y12*y23+y13*y23)/
95  & (y12+y23)**2+y23l*(4d0*y12**2+2d0*y12*y23+4d0*y12*y13+
96  & y13*y23)/(y12+y13)**2)/wt1+
97  & cn*(y13l*y13/(y12+y23)+y23l*y23/(y12+y13))/wt1+(cn-2d0*cf)*
98  & ((y12**2+(y12+y13)**2)*(y12l*y23l-y12l*y12m-y23l*
99  & y23m+1.644934d0-y12i-y23i)/(y13*y23)+(y12**2+(y12+y23)**2)*
100  & (y12l*y13l-y12l*y12m-y13l*y13m+1.644934d0-y12i-y13i)/
101  & (y13*y23)+(y13**2+y23**2)/(y13*y23*(y13+y23))-
102  & 2d0*y12l*y12**2/(y13+y23)**2-4d0*y12l*y12/(y13+y23))/wt1-
103  & cn*(y13l*y23l-y13l*y13m-y23l*y23m+1.644934d0-y13i-y23i)
104  IF(1d0+wtopt+als2pi*wt2.LE.0d0) mstj(121)=1
105  IF(1d0+wtopt+als2pi*wt2.LE.wtmax*pyr(0)) goto 110
106  parj(156)=(wtopt+als2pi*wt2)/(1d0+wtopt+als2pi*wt2)
107 
108  ELSEIF(mstj(101).EQ.2.AND.mstj(110).EQ.2) THEN
109 C...Second order corrections; Zhu parametrization of ERT.
110  zx=(y23-y13)**2
111  zy=1d0-y12
112  iza=0
113  DO 120 iy=1,5
114  IF(abs(cut-0.01d0*iy).LT.0.0001d0) iza=iy
115  120 CONTINUE
116  IF(iza.NE.0) THEN
117  iz=iza
118  wt2=zhup(iz,1)+zhup(iz,2)*zx+zhup(iz,3)*zx**2+(zhup(iz,4)+
119  & zhup(iz,5)*zx)*zy+(zhup(iz,6)+zhup(iz,7)*zx)*zy**2+
120  & (zhup(iz,8)+zhup(iz,9)*zx)*zy**3+zhup(iz,10)/(zx-zy**2)+
121  & zhup(iz,11)/(1d0-zy)+zhup(iz,12)/zy
122  ELSE
123  iz=100d0*cut
124  wtl=zhup(iz,1)+zhup(iz,2)*zx+zhup(iz,3)*zx**2+(zhup(iz,4)+
125  & zhup(iz,5)*zx)*zy+(zhup(iz,6)+zhup(iz,7)*zx)*zy**2+
126  & (zhup(iz,8)+zhup(iz,9)*zx)*zy**3+zhup(iz,10)/(zx-zy**2)+
127  & zhup(iz,11)/(1d0-zy)+zhup(iz,12)/zy
128  iz=iz+1
129  wtu=zhup(iz,1)+zhup(iz,2)*zx+zhup(iz,3)*zx**2+(zhup(iz,4)+
130  & zhup(iz,5)*zx)*zy+(zhup(iz,6)+zhup(iz,7)*zx)*zy**2+
131  & (zhup(iz,8)+zhup(iz,9)*zx)*zy**3+zhup(iz,10)/(zx-zy**2)+
132  & zhup(iz,11)/(1d0-zy)+zhup(iz,12)/zy
133  wt2=wtl+(wtu-wtl)*(100d0*cut+1d0-iz)
134  ENDIF
135  IF(1d0+wtopt+2d0*als2pi*wt2.LE.0d0) mstj(121)=1
136  IF(1d0+wtopt+2d0*als2pi*wt2.LE.wtmax*pyr(0)) goto 110
137  parj(156)=(wtopt+2d0*als2pi*wt2)/(1d0+wtopt+2d0*als2pi*wt2)
138  ENDIF
139 
140 C...Impose mass cuts (gives two jets). For fixed jet number new try.
141  x1=1d0-y23
142  x2=1d0-y13
143  x3=1d0-y12
144  IF(4d0*y23*y13*y12/x3**2.LE.qme) njet=2
145  IF(mod(mstj(103),4).GE.2.AND.iabs(mstj(101)).LE.1.AND.qme*x3+
146  & 0.5d0*qme**2+(0.5d0*qme+0.25d0*qme**2)*((1d0-x2)/(1d0-x1)+
147  & (1d0-x1)/(1d0-x2)).GT.(x1**2+x2**2)*pyr(0)) njet=2
148  IF(mstj(101).EQ.-1.AND.njet.EQ.2) goto 100
149 
150 C...Scalar gluon model (first order only, no mass effects).
151  ELSE
152  130 njet=3
153  140 x3=sqrt(4d0*cut**2+pyr(0)*((1d0-cut)**2-4d0*cut**2))
154  IF(log((x3-cut)/cut).LE.pyr(0)*log((1d0-2d0*cut)/cut)) goto 140
155  yd=sign(2d0*cut*((x3-cut)/cut)**pyr(0)-x3,pyr(0)-0.5d0)
156  x1=1d0-0.5d0*(x3+yd)
157  x2=1d0-0.5d0*(x3-yd)
158  IF(4d0*(1d0-x1)*(1d0-x2)*(1d0-x3)/x3**2.LE.qme) njet=2
159  IF(mstj(102).GE.2) THEN
160  IF(x3**2-2d0*(1d0+x3)*(1d0-x1)*(1d0-x2)*parj(171).LT.
161  & x3**2*pyr(0)) njet=2
162  ENDIF
163  IF(mstj(101).EQ.-1.AND.njet.EQ.2) goto 130
164  ENDIF
165 
166  RETURN
167  END