Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luxjet.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luxjet.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luxjet(ECM,NJET,CUT)
5 
6 C...Purpose: to select number of jets in matrix element approach.
7  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
8  SAVE /ludat1/
9  dimension zhut(5)
10 
11 C...Relative three-jet rate in Zhu second order parametrization.
12  DATA zhut/3.0922, 6.2291, 7.4782, 7.8440, 8.2560/
13 
14 C...Trivial result for two-jets only, including parton shower.
15  IF(mstj(101).EQ.0.OR.mstj(101).EQ.5) THEN
16  cut=0.
17 
18 C...QCD and Abelian vector gluon theory: Q^2 for jet rate and R.
19  ELSEIF(mstj(109).EQ.0.OR.mstj(109).EQ.2) THEN
20  cf=4./3.
21  IF(mstj(109).EQ.2) cf=1.
22  IF(mstj(111).EQ.0) THEN
23  q2=ecm**2
24  q2r=ecm**2
25  ELSEIF(mstu(111).EQ.0) THEN
26  parj(169)=min(1.,parj(129))
27  q2=parj(169)*ecm**2
28  parj(168)=min(1.,max(parj(128),exp(-12.*paru(1)/
29  & ((33.-2.*mstu(112))*paru(111)))))
30  q2r=parj(168)*ecm**2
31  ELSE
32  parj(169)=min(1.,max(parj(129),(2.*paru(112)/ecm)**2))
33  q2=parj(169)*ecm**2
34  parj(168)=min(1.,max(parj(128),paru(112)/ecm,
35  & (2.*paru(112)/ecm)**2))
36  q2r=parj(168)*ecm**2
37  ENDIF
38 
39 C...alpha_strong for R and R itself.
40  alspi=(3./4.)*cf*ulalps(q2r)/paru(1)
41  IF(iabs(mstj(101)).EQ.1) THEN
42  rqcd=1.+alspi
43  ELSEIF(mstj(109).EQ.0) THEN
44  rqcd=1.+alspi+(1.986-0.115*mstu(118))*alspi**2
45  IF(mstj(111).EQ.1) rqcd=max(1.,rqcd+(33.-2.*mstu(112))/12.*
46  & log(parj(168))*alspi**2)
47  ELSE
48  rqcd=1.+alspi-(3./32.+0.519*mstu(118))*(4.*alspi/3.)**2
49  ENDIF
50 
51 C...alpha_strong for jet rate. Initial value for y cut.
52  alspi=(3./4.)*cf*ulalps(q2)/paru(1)
53  cut=max(0.001,parj(125),(parj(126)/ecm)**2)
54  IF(iabs(mstj(101)).LE.1.OR.(mstj(109).EQ.0.AND.mstj(111).EQ.0))
55  & cut=max(cut,exp(-sqrt(0.75/alspi))/2.)
56  IF(mstj(110).EQ.2) cut=max(0.01,min(0.05,cut))
57 
58 C...Parametrization of first order three-jet cross-section.
59  100 IF(mstj(101).EQ.0.OR.cut.GE.0.25) THEN
60  parj(152)=0.
61  ELSE
62  parj(152)=(2.*alspi/3.)*((3.-6.*cut+2.*log(cut))*
63  & log(cut/(1.-2.*cut))+(2.5+1.5*cut-6.571)*(1.-3.*cut)+
64  & 5.833*(1.-3.*cut)**2-3.894*(1.-3.*cut)**3+
65  & 1.342*(1.-3.*cut)**4)/rqcd
66  IF(mstj(109).EQ.2.AND.(mstj(101).EQ.2.OR.mstj(101).LE.-2))
67  & parj(152)=0.
68  ENDIF
69 
70 C...Parametrization of second order three-jet cross-section.
71  IF(iabs(mstj(101)).LE.1.OR.mstj(101).EQ.3.OR.mstj(109).EQ.2.OR.
72  & cut.GE.0.25) THEN
73  parj(153)=0.
74  ELSEIF(mstj(110).LE.1) THEN
75  ct=log(1./cut-2.)
76  parj(153)=alspi**2*ct**2*(2.419+0.5989*ct+0.6782*ct**2-
77  & 0.2661*ct**3+0.01159*ct**4)/rqcd
78 
79 C...Interpolation in second/first order ratio for Zhu parametrization.
80  ELSEIF(mstj(110).EQ.2) THEN
81  iza=0
82  DO 110 iy=1,5
83  110 IF(abs(cut-0.01*iy).LT.0.0001) iza=iy
84  IF(iza.NE.0) THEN
85  zhurat=zhut(iza)
86  ELSE
87  iz=100.*cut
88  zhurat=zhut(iz)+(100.*cut-iz)*(zhut(iz+1)-zhut(iz))
89  ENDIF
90  parj(153)=alspi*parj(152)*zhurat
91  ENDIF
92 
93 C...Shift in second order three-jet cross-section with optimized Q^2.
94  IF(mstj(111).EQ.1.AND.iabs(mstj(101)).GE.2.AND.mstj(101).NE.3.
95  & and.cut.LT.0.25) parj(153)=parj(153)+(33.-2.*mstu(112))/12.*
96  & log(parj(169))*alspi*parj(152)
97 
98 C...Parametrization of second order four-jet cross-section.
99  IF(iabs(mstj(101)).LE.1.OR.cut.GE.0.125) THEN
100  parj(154)=0.
101  ELSE
102  ct=log(1./cut-5.)
103  IF(cut.LE.0.018) THEN
104  xqqgg=6.349-4.330*ct+0.8304*ct**2
105  IF(mstj(109).EQ.2) xqqgg=(4./3.)**2*(3.035-2.091*ct+
106  & 0.4059*ct**2)
107  xqqqq=1.25*(-0.1080+0.01486*ct+0.009364*ct**2)
108  IF(mstj(109).EQ.2) xqqqq=8.*xqqqq
109  ELSE
110  xqqgg=-0.09773+0.2959*ct-0.2764*ct**2+0.08832*ct**3
111  IF(mstj(109).EQ.2) xqqgg=(4./3.)**2*(-0.04079+0.1340*ct-
112  & 0.1326*ct**2+0.04365*ct**3)
113  xqqqq=1.25*(0.003661-0.004888*ct-0.001081*ct**2+0.002093*
114  & ct**3)
115  IF(mstj(109).EQ.2) xqqqq=8.*xqqqq
116  ENDIF
117  parj(154)=alspi**2*ct**2*(xqqgg+xqqqq)/rqcd
118  parj(155)=xqqqq/(xqqgg+xqqqq)
119  ENDIF
120 
121 C...If negative three-jet rate, change y' optimization parameter.
122  IF(mstj(111).EQ.1.AND.parj(152)+parj(153).LT.0..AND.
123  & parj(169).LT.0.99) THEN
124  parj(169)=min(1.,1.2*parj(169))
125  q2=parj(169)*ecm**2
126  alspi=(3./4.)*cf*ulalps(q2)/paru(1)
127  goto 100
128  ENDIF
129 
130 C...If too high cross-section, use harder cuts, or fail.
131  IF(parj(152)+parj(153)+parj(154).GE.1) THEN
132  IF(mstj(110).EQ.2.AND.cut.GT.0.0499.AND.mstj(111).EQ.1.AND.
133  & parj(169).LT.0.99) THEN
134  parj(169)=min(1.,1.2*parj(169))
135  q2=parj(169)*ecm**2
136  alspi=(3./4.)*cf*ulalps(q2)/paru(1)
137  goto 100
138  ELSEIF(mstj(110).EQ.2.AND.cut.GT.0.0499) THEN
139  CALL luerrm(26,
140  & '(LUXJET:) no allowed y cut value for Zhu parametrization')
141  ENDIF
142  cut=0.26*(4.*cut)**(parj(152)+parj(153)+parj(154))**(-1./3.)
143  IF(mstj(110).EQ.2) cut=max(0.01,min(0.05,cut))
144  goto 100
145  ENDIF
146 
147 C...Scalar gluon (first order only).
148  ELSE
149  alspi=ulalps(ecm**2)/paru(1)
150  cut=max(0.001,parj(125),(parj(126)/ecm)**2,exp(-3./alspi))
151  parj(152)=0.
152  IF(cut.LT.0.25) parj(152)=(alspi/3.)*((1.-2.*cut)*
153  & log((1.-2.*cut)/cut)+0.5*(9.*cut**2-1.))
154  parj(153)=0.
155  parj(154)=0.
156  ENDIF
157 
158 C...Select number of jets.
159  parj(150)=cut
160  IF(mstj(101).EQ.0.OR.mstj(101).EQ.5) THEN
161  njet=2
162  ELSEIF(mstj(101).LE.0) THEN
163  njet=min(4,2-mstj(101))
164  ELSE
165  rnj=rlu(0)
166  njet=2
167  IF(parj(152)+parj(153)+parj(154).GT.rnj) njet=3
168  IF(parj(154).GT.rnj) njet=4
169  ENDIF
170 
171  RETURN
172  END