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