]> git.uio.no Git - u/mrichter/AliRoot.git/blob - HIJING/hipyset1_35/lux4jt_hijing.F
Protection against division by 0 added.
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / lux4jt_hijing.F
1 * $Id$
2     
3 C*********************************************************************  
4     
5       SUBROUTINE LUX4JT_HIJING(NJET,CUT,KFL,ECM,KFLN,X1,X2,X4,X12,X14) 
6     
7 C...Purpose: to select the kinematical variables of four-jet events.    
8 #include "ludat1_hijing.inc"
9       DIMENSION WTA(4),WTB(4),WTC(4),WTD(4),WTE(4)  
10     
11 C...Common constants. Colour factors for QCD and Abelian gluon theory.  
12       PMQ=ULMASS_HIJING(KFL)   
13       QME=(2.*PMQ/ECM)**2   
14       CT=LOG(1./CUT-5.) 
15       IF(MSTJ(109).EQ.0) THEN   
16         CF=4./3.    
17         CN=3.   
18         TR=2.   
19       ELSE  
20         CF=1.   
21         CN=0.   
22         TR=12.  
23       ENDIF 
24     
25 C...Choice of process (qqbargg or qqbarqqbar).  
26   100 NJET=4    
27       IT=1  
28       IF(PARJ(155).GT.RLU_HIJING(0)) IT=2  
29       IF(MSTJ(101).LE.-3) IT=-MSTJ(101)-2   
30       IF(IT.EQ.1) WTMX=0.7/CUT**2   
31       IF(IT.EQ.1.AND.MSTJ(109).EQ.2) WTMX=0.6/CUT**2    
32       IF(IT.EQ.2) WTMX=0.1125*CF*TR/CUT**2  
33       ID=1  
34     
35 C...Sample the five kinematical variables (for qqgg preweighted in y34).    
36   110 Y134=3.*CUT+(1.-6.*CUT)*RLU_HIJING(0)    
37       Y234=3.*CUT+(1.-6.*CUT)*RLU_HIJING(0)    
38       IF(IT.EQ.1) Y34=(1.-5.*CUT)*EXP(-CT*RLU_HIJING(0))   
39       IF(IT.EQ.2) Y34=CUT+(1.-6.*CUT)*RLU_HIJING(0)    
40       IF(Y34.LE.Y134+Y234-1..OR.Y34.GE.Y134*Y234) GOTO 110  
41       VT=RLU_HIJING(0) 
42       CP=COS(PARU(1)*RLU_HIJING(0))    
43       Y14=(Y134-Y34)*VT 
44       Y13=Y134-Y14-Y34  
45       VB=Y34*(1.-Y134-Y234+Y34)/((Y134-Y34)*(Y234-Y34)) 
46       Y24=0.5*(Y234-Y34)*(1.-4.*SQRT(MAX(0.,VT*(1.-VT)*VB*(1.-VB)))*    
47      &CP-(1.-2.*VT)*(1.-2.*VB)) 
48       Y23=Y234-Y34-Y24  
49       Y12=1.-Y134-Y23-Y24   
50       IF(MIN(Y12,Y13,Y14,Y23,Y24).LE.CUT) GOTO 110  
51       Y123=Y12+Y13+Y23  
52       Y124=Y12+Y14+Y24  
53     
54 C...Calculate matrix elements for qqgg or qqqq process. 
55       IC=0  
56       WTTOT=0.  
57   120 IC=IC+1   
58       IF(IT.EQ.1) THEN  
59         WTA(IC)=(Y12*Y34**2-Y13*Y24*Y34+Y14*Y23*Y34+3.*Y12*Y23*Y34+ 
60      &  3.*Y12*Y14*Y34+4.*Y12**2*Y34-Y13*Y23*Y24+2.*Y12*Y23*Y24-    
61      &  Y13*Y14*Y24-2.*Y12*Y13*Y24+2.*Y12**2*Y24+Y14*Y23**2+2.*Y12* 
62      &  Y23**2+Y14**2*Y23+4.*Y12*Y14*Y23+4.*Y12**2*Y23+2.*Y12*Y14**2+   
63      &  2.*Y12*Y13*Y14+4.*Y12**2*Y14+2.*Y12**2*Y13+2.*Y12**3)/(2.*Y13*  
64      &  Y134*Y234*Y24)+(Y24*Y34+Y12*Y34+Y13*Y24-Y14*Y23+Y12*Y13)/(Y13*  
65      &  Y134**2)+2.*Y23*(1.-Y13)/(Y13*Y134*Y24)+Y34/(2.*Y13*Y24)    
66         WTB(IC)=(Y12*Y24*Y34+Y12*Y14*Y34-Y13*Y24**2+Y13*Y14*Y24+2.*Y12* 
67      &  Y14*Y24)/(Y13*Y134*Y23*Y14)+Y12*(1.+Y34)*Y124/(Y134*Y234*Y14*   
68      &  Y24)-(2.*Y13*Y24+Y14**2+Y13*Y23+2.*Y12*Y13)/(Y13*Y134*Y14)+ 
69      &  Y12*Y123*Y124/(2.*Y13*Y14*Y23*Y24)  
70         WTC(IC)=-(5.*Y12*Y34**2+2.*Y12*Y24*Y34+2.*Y12*Y23*Y34+2.*Y12*   
71      &  Y14*Y34+2.*Y12*Y13*Y34+4.*Y12**2*Y34-Y13*Y24**2+Y14*Y23*Y24+    
72      &  Y13*Y23*Y24+Y13*Y14*Y24-Y12*Y14*Y24-Y13**2*Y24-3.*Y12*Y13*Y24-  
73      &  Y14*Y23**2-Y14**2*Y23+Y13*Y14*Y23-3.*Y12*Y14*Y23-Y12*Y13*Y23)/  
74      &  (4.*Y134*Y234*Y34**2)+(3.*Y12*Y34**2-3.*Y13*Y24*Y34+3.*Y12*Y24* 
75      &  Y34+3.*Y14*Y23*Y34-Y13*Y24**2-Y12*Y23*Y34+6.*Y12*Y14*Y34+2.*Y12*    
76      &  Y13*Y34-2.*Y12**2*Y34+Y14*Y23*Y24-3.*Y13*Y23*Y24-2.*Y13*Y14*    
77      &  Y24+4.*Y12*Y14*Y24+2.*Y12*Y13*Y24+3.*Y14*Y23**2+2.*Y14**2*Y23+  
78      &  2.*Y14**2*Y12+2.*Y12**2*Y14+6.*Y12*Y14*Y23-2.*Y12*Y13**2-   
79      &  2.*Y12**2*Y13)/(4.*Y13*Y134*Y234*Y34)   
80         WTC(IC)=WTC(IC)+(2.*Y12*Y34**2-2.*Y13*Y24*Y34+Y12*Y24*Y34+  
81      &  4.*Y13*Y23*Y34+4.*Y12*Y14*Y34+2.*Y12*Y13*Y34+2.*Y12**2*Y34- 
82      &  Y13*Y24**2+3.*Y14*Y23*Y24+4.*Y13*Y23*Y24-2.*Y13*Y14*Y24+    
83      &  4.*Y12*Y14*Y24+2.*Y12*Y13*Y24+2.*Y14*Y23**2+4.*Y13*Y23**2+  
84      &  2.*Y13*Y14*Y23+2.*Y12*Y14*Y23+4.*Y12*Y13*Y23+2.*Y12*Y14**2+4.*  
85      &  Y12**2*Y13+4.*Y12*Y13*Y14+2.*Y12**2*Y14)/(4.*Y13*Y134*Y24*Y34)- 
86      &  (Y12*Y34**2-2.*Y14*Y24*Y34-2.*Y13*Y24*Y34-Y14*Y23*Y34+Y13*Y23*  
87      &  Y34+Y12*Y14*Y34+2.*Y12*Y13*Y34-2.*Y14**2*Y24-4.*Y13*Y14*Y24-    
88      &  4.*Y13**2*Y24-Y14**2*Y23-Y13**2*Y23+Y12*Y13*Y14-Y12*Y13**2)/    
89      &  (2.*Y13*Y34*Y134**2)+(Y12*Y34**2-4.*Y14*Y24*Y34-2.*Y13*Y24*Y34- 
90      &  2.*Y14*Y23*Y34-4.*Y13*Y23*Y34-4.*Y12*Y14*Y34-4.*Y12*Y13*Y34-    
91      &  2.*Y13*Y14*Y24+2.*Y13**2*Y24+2.*Y14**2*Y23-2.*Y13*Y14*Y23-  
92      &  Y12*Y14**2-6.*Y12*Y13*Y14-Y12*Y13**2)/(4.*Y34**2*Y134**2)   
93         WTTOT=WTTOT+Y34*CF*(CF*WTA(IC)+(CF-0.5*CN)*WTB(IC)+CN*WTC(IC))/ 
94      &  8.  
95       ELSE  
96         WTD(IC)=(Y13*Y23*Y34+Y12*Y23*Y34-Y12**2*Y34+Y13*Y23*Y24+2.*Y12* 
97      &  Y23*Y24-Y14*Y23**2+Y12*Y13*Y24+Y12*Y14*Y23+Y12*Y13*Y14)/(Y13**2*    
98      &  Y123**2)-(Y12*Y34**2-Y13*Y24*Y34+Y12*Y24*Y34-Y14*Y23*Y34-Y12*   
99      &  Y23*Y34-Y13*Y24**2+Y14*Y23*Y24-Y13*Y23*Y24-Y13**2*Y24+Y14*  
100      &  Y23**2)/(Y13**2*Y123*Y134)+(Y13*Y14*Y12+Y34*Y14*Y12-Y34**2*Y12+ 
101      &  Y13*Y14*Y24+2.*Y34*Y14*Y24-Y23*Y14**2+Y34*Y13*Y24+Y34*Y23*Y14+  
102      &  Y34*Y13*Y23)/(Y13**2*Y134**2)-(Y34*Y12**2-Y13*Y24*Y12+Y34*Y24*  
103      &  Y12-Y23*Y14*Y12-Y34*Y14*Y12-Y13*Y24**2+Y23*Y14*Y24-Y13*Y14*Y24- 
104      &  Y13**2*Y24+Y23*Y14**2)/(Y13**2*Y134*Y123)   
105         WTE(IC)=(Y12*Y34*(Y23-Y24+Y14+Y13)+Y13*Y24**2-Y14*Y23*Y24+Y13*  
106      &  Y23*Y24+Y13*Y14*Y24+Y13**2*Y24-Y14*Y23*(Y14+Y23+Y13))/(Y13*Y23* 
107      &  Y123*Y134)-Y12*(Y12*Y34-Y23*Y24-Y13*Y24-Y14*Y23-Y14*Y13)/(Y13*  
108      &  Y23*Y123**2)-(Y14+Y13)*(Y24+Y23)*Y34/(Y13*Y23*Y134*Y234)+   
109      &  (Y12*Y34*(Y14-Y24+Y23+Y13)+Y13*Y24**2-Y23*Y14*Y24+Y13*Y14*Y24+  
110      &  Y13*Y23*Y24+Y13**2*Y24-Y23*Y14*(Y14+Y23+Y13))/(Y13*Y14*Y134*    
111      &  Y123)-Y34*(Y34*Y12-Y14*Y24-Y13*Y24-Y23*Y14-Y23*Y13)/(Y13*Y14*   
112      &  Y134**2)-(Y23+Y13)*(Y24+Y14)*Y12/(Y13*Y14*Y123*Y124)    
113         WTTOT=WTTOT+CF*(TR*WTD(IC)+(CF-0.5*CN)*WTE(IC))/16. 
114       ENDIF 
115     
116 C...Permutations of momenta in matrix element. Weighting.   
117   130 IF(IC.EQ.1.OR.IC.EQ.3.OR.ID.EQ.2.OR.ID.EQ.3) THEN 
118         YSAV=Y13    
119         Y13=Y14 
120         Y14=YSAV    
121         YSAV=Y23    
122         Y23=Y24 
123         Y24=YSAV    
124         YSAV=Y123   
125         Y123=Y124   
126         Y124=YSAV   
127       ENDIF 
128       IF(IC.EQ.2.OR.IC.EQ.4.OR.ID.EQ.3.OR.ID.EQ.4) THEN 
129         YSAV=Y13    
130         Y13=Y23 
131         Y23=YSAV    
132         YSAV=Y14    
133         Y14=Y24 
134         Y24=YSAV    
135         YSAV=Y134   
136         Y134=Y234   
137         Y234=YSAV   
138       ENDIF 
139       IF(IC.LE.3) GOTO 120  
140       IF(ID.EQ.1.AND.WTTOT.LT.RLU_HIJING(0)*WTMX) GOTO 110 
141       IC=5  
142     
143 C...qqgg events: string configuration and event type.   
144       IF(IT.EQ.1) THEN  
145         IF(MSTJ(109).EQ.0.AND.ID.EQ.1) THEN 
146           PARJ(156)=Y34*(2.*(WTA(1)+WTA(2)+WTA(3)+WTA(4))+4.*(WTC(1)+   
147      &    WTC(2)+WTC(3)+WTC(4)))/(9.*WTTOT) 
148           IF(WTA(2)+WTA(4)+2.*(WTC(2)+WTC(4)).GT.RLU_HIJING(0)*(WTA(1)
149      $         +WTA(2)+WTA(3)+WTA(4)+2.*(WTC(1)+WTC(2)+WTC(3)+WTC(4))))
150      $         ID=2 
151           IF(ID.EQ.2) GOTO 130  
152         ELSEIF(MSTJ(109).EQ.2.AND.ID.EQ.1) THEN 
153           PARJ(156)=Y34*(WTA(1)+WTA(2)+WTA(3)+WTA(4))/(8.*WTTOT)    
154           IF(WTA(2)+WTA(4).GT.RLU_HIJING(0)*(WTA(1)+WTA(2)+WTA(3)+WTA(4)
155      $         )) ID=2    
156           IF(ID.EQ.2) GOTO 130  
157         ENDIF   
158         MSTJ(120)=3 
159         IF(MSTJ(109).EQ.0.AND.0.5*Y34*(WTC(1)+WTC(2)+WTC(3)+WTC(4)).GT. 
160      &  RLU_HIJING(0)*WTTOT) MSTJ(120)=4   
161         KFLN=21 
162     
163 C...Mass cuts. Kinematical variables out.   
164         IF(Y12.LE.CUT+QME) NJET=2   
165         IF(NJET.EQ.2) GOTO 150  
166         Q12=0.5*(1.-SQRT(1.-QME/Y12))   
167         X1=1.-(1.-Q12)*Y234-Q12*Y134    
168         X4=1.-(1.-Q12)*Y134-Q12*Y234    
169         X2=1.-Y124  
170         X12=(1.-Q12)*Y13+Q12*Y23    
171         X14=Y12-0.5*QME 
172         IF(Y134*Y234/((1.-X1)*(1.-X4)).LE.RLU_HIJING(0)) NJET=2    
173     
174 C...qqbarqqbar events: string configuration, choose new flavour.    
175       ELSE  
176         IF(ID.EQ.1) THEN    
177           WTR=RLU_HIJING(0)*(WTD(1)+WTD(2)+WTD(3)+WTD(4))  
178           IF(WTR.LT.WTD(2)+WTD(3)+WTD(4)) ID=2  
179           IF(WTR.LT.WTD(3)+WTD(4)) ID=3 
180           IF(WTR.LT.WTD(4)) ID=4    
181           IF(ID.GE.2) GOTO 130  
182         ENDIF   
183         MSTJ(120)=5 
184         PARJ(156)=CF*TR*(WTD(1)+WTD(2)+WTD(3)+WTD(4))/(16.*WTTOT)   
185   140   KFLN=1+INT(5.*RLU_HIJING(0))   
186         IF(KFLN.NE.KFL.AND.0.2*PARJ(156).LE.RLU_HIJING(0)) GOTO 140    
187         IF(KFLN.EQ.KFL.AND.1.-0.8*PARJ(156).LE.RLU_HIJING(0)) GOTO 140 
188         IF(KFLN.GT.MSTJ(104)) NJET=2    
189         PMQN=ULMASS_HIJING(KFLN)   
190         QMEN=(2.*PMQN/ECM)**2   
191     
192 C...Mass cuts. Kinematical variables out.   
193         IF(Y24.LE.CUT+QME.OR.Y13.LE.1.1*QMEN) NJET=2    
194         IF(NJET.EQ.2) GOTO 150  
195         Q24=0.5*(1.-SQRT(1.-QME/Y24))   
196         Q13=0.5*(1.-SQRT(1.-QMEN/Y13))  
197         X1=1.-(1.-Q24)*Y123-Q24*Y134    
198         X4=1.-(1.-Q24)*Y134-Q24*Y123    
199         X2=1.-(1.-Q13)*Y234-Q13*Y124    
200         X12=(1.-Q24)*((1.-Q13)*Y14+Q13*Y34)+Q24*((1.-Q13)*Y12+Q13*Y23)  
201         X14=Y24-0.5*QME 
202         X34=(1.-Q24)*((1.-Q13)*Y23+Q13*Y12)+Q24*((1.-Q13)*Y34+Q13*Y14)  
203         IF(PMQ**2+PMQN**2+MIN(X12,X34)*ECM**2.LE.   
204      &  (PARJ(127)+PMQ+PMQN)**2) NJET=2 
205         IF(Y123*Y134/((1.-X1)*(1.-X4)).LE.RLU_HIJING(0)) NJET=2    
206       ENDIF 
207   150 IF(MSTJ(101).LE.-2.AND.NJET.EQ.2) GOTO 100    
208     
209       RETURN    
210       END