]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PHOS/shaker/luxjet.f
Syntax problems on HP-UX corrected
[u/mrichter/AliRoot.git] / PHOS / shaker / luxjet.f
1 *CMZ :          17/07/98  15.44.35  by  Federico Carminati
2 *-- Author :
3 C*********************************************************************
4
5       SUBROUTINE LUXJET(ECM,NJET,CUT)
6
7 C...Purpose: to select number of jets in matrix element approach.
8 *KEEP,LUDAT1.
9       COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
10       SAVE /LUDAT1/
11 *KEND.
12       DIMENSION ZHUT(5)
13
14 C...Relative three-jet rate in Zhu second order parametrization.
15       DATA ZHUT/3.0922, 6.2291, 7.4782, 7.8440, 8.2560/
16
17 C...Trivial result for two-jets only, including parton shower.
18       IF(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.5) THEN
19         CUT=0.
20
21 C...QCD and Abelian vector gluon theory: Q^2 for jet rate and R.
22       ELSEIF(MSTJ(109).EQ.0.OR.MSTJ(109).EQ.2) THEN
23         CF=4./3.
24         IF(MSTJ(109).EQ.2) CF=1.
25         IF(MSTJ(111).EQ.0) THEN
26           Q2=ECM**2
27           Q2R=ECM**2
28         ELSEIF(MSTU(111).EQ.0) THEN
29           PARJ(169)=MIN(1.,PARJ(129))
30           Q2=PARJ(169)*ECM**2
31           PARJ(168)=MIN(1.,MAX(PARJ(128),EXP(-12.*PARU(1)/
32      &    ((33.-2.*MSTU(112))*PARU(111)))))
33           Q2R=PARJ(168)*ECM**2
34         ELSE
35           PARJ(169)=MIN(1.,MAX(PARJ(129),(2.*PARU(112)/ECM)**2))
36           Q2=PARJ(169)*ECM**2
37           PARJ(168)=MIN(1.,MAX(PARJ(128),PARU(112)/ECM,
38      &    (2.*PARU(112)/ECM)**2))
39           Q2R=PARJ(168)*ECM**2
40         ENDIF
41
42 C...alpha_strong for R and R itself.
43         ALSPI=(3./4.)*CF*ULALPS(Q2R)/PARU(1)
44         IF(IABS(MSTJ(101)).EQ.1) THEN
45           RQCD=1.+ALSPI
46         ELSEIF(MSTJ(109).EQ.0) THEN
47           RQCD=1.+ALSPI+(1.986-0.115*MSTU(118))*ALSPI**2
48           IF(MSTJ(111).EQ.1) RQCD=MAX(1.,RQCD+(33.-2.*MSTU(112))/12.*
49      &    LOG(PARJ(168))*ALSPI**2)
50         ELSE
51           RQCD=1.+ALSPI-(3./32.+0.519*MSTU(118))*(4.*ALSPI/3.)**2
52         ENDIF
53
54 C...alpha_strong for jet rate. Initial value for y cut.
55         ALSPI=(3./4.)*CF*ULALPS(Q2)/PARU(1)
56         CUT=MAX(0.001,PARJ(125),(PARJ(126)/ECM)**2)
57         IF(IABS(MSTJ(101)).LE.1.OR.(MSTJ(109).EQ.0.AND.MSTJ(111).EQ.0))
58      &  CUT=MAX(CUT,EXP(-SQRT(0.75/ALSPI))/2.)
59         IF(MSTJ(110).EQ.2) CUT=MAX(0.01,MIN(0.05,CUT))
60
61 C...Parametrization of first order three-jet cross-section.
62   100   IF(MSTJ(101).EQ.0.OR.CUT.GE.0.25) THEN
63           PARJ(152)=0.
64         ELSE
65           PARJ(152)=(2.*ALSPI/3.)*((3.-6.*CUT+2.*LOG(CUT))*
66      &    LOG(CUT/(1.-2.*CUT))+(2.5+1.5*CUT-6.571)*(1.-3.*CUT)+
67      &    5.833*(1.-3.*CUT)**2-3.894*(1.-3.*CUT)**3+
68      &    1.342*(1.-3.*CUT)**4)/RQCD
69           IF(MSTJ(109).EQ.2.AND.(MSTJ(101).EQ.2.OR.MSTJ(101).LE.-2))
70      &    PARJ(152)=0.
71         ENDIF
72
73 C...Parametrization of second order three-jet cross-section.
74         IF(IABS(MSTJ(101)).LE.1.OR.MSTJ(101).EQ.3.OR.MSTJ(109).EQ.2.OR.
75      &  CUT.GE.0.25) THEN
76           PARJ(153)=0.
77         ELSEIF(MSTJ(110).LE.1) THEN
78           CT=LOG(1./CUT-2.)
79           PARJ(153)=ALSPI**2*CT**2*(2.419+0.5989*CT+0.6782*CT**2-
80      &    0.2661*CT**3+0.01159*CT**4)/RQCD
81
82 C...Interpolation in second/first order ratio for Zhu parametrization.
83         ELSEIF(MSTJ(110).EQ.2) THEN
84           IZA=0
85           DO 110 IY=1,5
86   110     IF(ABS(CUT-0.01*IY).LT.0.0001) IZA=IY
87           IF(IZA.NE.0) THEN
88             ZHURAT=ZHUT(IZA)
89           ELSE
90             IZ=100.*CUT
91             ZHURAT=ZHUT(IZ)+(100.*CUT-IZ)*(ZHUT(IZ+1)-ZHUT(IZ))
92           ENDIF
93           PARJ(153)=ALSPI*PARJ(152)*ZHURAT
94         ENDIF
95
96 C...Shift in second order three-jet cross-section with optimized Q^2.
97         IF(MSTJ(111).EQ.1.AND.IABS(MSTJ(101)).GE.2.AND.MSTJ(101).NE.3.
98      &  AND.CUT.LT.0.25) PARJ(153)=PARJ(153)+(33.-2.*MSTU(112))/12.*
99      &  LOG(PARJ(169))*ALSPI*PARJ(152)
100
101 C...Parametrization of second order four-jet cross-section.
102         IF(IABS(MSTJ(101)).LE.1.OR.CUT.GE.0.125) THEN
103           PARJ(154)=0.
104         ELSE
105           CT=LOG(1./CUT-5.)
106           IF(CUT.LE.0.018) THEN
107             XQQGG=6.349-4.330*CT+0.8304*CT**2
108             IF(MSTJ(109).EQ.2) XQQGG=(4./3.)**2*(3.035-2.091*CT+
109      &      0.4059*CT**2)
110             XQQQQ=1.25*(-0.1080+0.01486*CT+0.009364*CT**2)
111             IF(MSTJ(109).EQ.2) XQQQQ=8.*XQQQQ
112           ELSE
113             XQQGG=-0.09773+0.2959*CT-0.2764*CT**2+0.08832*CT**3
114             IF(MSTJ(109).EQ.2) XQQGG=(4./3.)**2*(-0.04079+0.1340*CT-
115      &      0.1326*CT**2+0.04365*CT**3)
116             XQQQQ=1.25*(0.003661-0.004888*CT-0.001081*CT**2+0.002093*
117      &      CT**3)
118             IF(MSTJ(109).EQ.2) XQQQQ=8.*XQQQQ
119           ENDIF
120           PARJ(154)=ALSPI**2*CT**2*(XQQGG+XQQQQ)/RQCD
121           PARJ(155)=XQQQQ/(XQQGG+XQQQQ)
122         ENDIF
123
124 C...If negative three-jet rate, change y' optimization parameter.
125         IF(MSTJ(111).EQ.1.AND.PARJ(152)+PARJ(153).LT.0..AND.
126      &  PARJ(169).LT.0.99) THEN
127           PARJ(169)=MIN(1.,1.2*PARJ(169))
128           Q2=PARJ(169)*ECM**2
129           ALSPI=(3./4.)*CF*ULALPS(Q2)/PARU(1)
130           GOTO 100
131         ENDIF
132
133 C...If too high cross-section, use harder cuts, or fail.
134         IF(PARJ(152)+PARJ(153)+PARJ(154).GE.1) THEN
135           IF(MSTJ(110).EQ.2.AND.CUT.GT.0.0499.AND.MSTJ(111).EQ.1.AND.
136      &    PARJ(169).LT.0.99) THEN
137             PARJ(169)=MIN(1.,1.2*PARJ(169))
138             Q2=PARJ(169)*ECM**2
139             ALSPI=(3./4.)*CF*ULALPS(Q2)/PARU(1)
140             GOTO 100
141           ELSEIF(MSTJ(110).EQ.2.AND.CUT.GT.0.0499) THEN
142             CALL LUERRM(26,
143      &      '(LUXJET:) no allowed y cut value for Zhu parametrization')
144           ENDIF
145           CUT=0.26*(4.*CUT)**(PARJ(152)+PARJ(153)+PARJ(154))**(-1./3.)
146           IF(MSTJ(110).EQ.2) CUT=MAX(0.01,MIN(0.05,CUT))
147           GOTO 100
148         ENDIF
149
150 C...Scalar gluon (first order only).
151       ELSE
152         ALSPI=ULALPS(ECM**2)/PARU(1)
153         CUT=MAX(0.001,PARJ(125),(PARJ(126)/ECM)**2,EXP(-3./ALSPI))
154         PARJ(152)=0.
155         IF(CUT.LT.0.25) PARJ(152)=(ALSPI/3.)*((1.-2.*CUT)*
156      &  LOG((1.-2.*CUT)/CUT)+0.5*(9.*CUT**2-1.))
157         PARJ(153)=0.
158         PARJ(154)=0.
159       ENDIF
160
161 C...Select number of jets.
162       PARJ(150)=CUT
163       IF(MSTJ(101).EQ.0.OR.MSTJ(101).EQ.5) THEN
164         NJET=2
165       ELSEIF(MSTJ(101).LE.0) THEN
166         NJET=MIN(4,2-MSTJ(101))
167       ELSE
168         RNJ=RLU(0)
169         NJET=2
170         IF(PARJ(152)+PARJ(153)+PARJ(154).GT.RNJ) NJET=3
171         IF(PARJ(154).GT.RNJ) NJET=4
172       ENDIF
173
174       RETURN
175       END