Numerical stability
[u/mrichter/AliRoot.git] / HIJING / hipyset1_35 / pyklim_hijing.F
1 * $Id$
2     
3 C***********************************************************************    
4     
5       SUBROUTINE PYKLIM_HIJING(ILIM)   
6     
7 C...Checks generated variables against pre-set kinematical limits;  
8 C...also calculates limits on variables used in generation. 
9 #include "ludat1_hijing.inc"
10 #include "ludat2_hijing.inc"
11 #include "ludat3_hijing.inc"
12 #include "pypars_hijing.inc"
13 #include "pysubs_hijing.inc"
14 #include "pyint1_hijing.inc"
15 #include "pyint2_hijing.inc"
16     
17 C...Common kinematical expressions. 
18       ISUB=MINT(1)  
19       IF(ISUB.EQ.96) GOTO 110   
20       SQM3=VINT(63) 
21       SQM4=VINT(64) 
22       IF(ILIM.NE.1) THEN    
23         TAU=VINT(21)    
24         RM3=SQM3/(TAU*VINT(2))  
25         RM4=SQM4/(TAU*VINT(2))  
26         BE34=SQRT((1.-RM3-RM4)**2-4.*RM3*RM4)   
27       ENDIF 
28       PTHMIN=CKIN(3)    
29       IF(MIN(SQM3,SQM4).LT.CKIN(6)**2) PTHMIN=MAX(CKIN(3),CKIN(5))  
30     
31       IF(ILIM.EQ.0) THEN    
32 C...Check generated values of tau, y*, cos(theta-hat), and tau' against 
33 C...pre-set kinematical limits. 
34         YST=VINT(22)    
35         CTH=VINT(23)    
36         TAUP=VINT(26)   
37         IF(ISET(ISUB).LE.2) THEN    
38           X1=SQRT(TAU)*EXP(YST) 
39           X2=SQRT(TAU)*EXP(-YST)    
40         ELSE    
41           X1=SQRT(TAUP)*EXP(YST)    
42           X2=SQRT(TAUP)*EXP(-YST)   
43         ENDIF   
44         XF=X1-X2    
45         IF(TAU*VINT(2).LT.CKIN(1)**2) MINT(51)=1    
46         IF(CKIN(2).GE.0..AND.TAU*VINT(2).GT.CKIN(2)**2) MINT(51)=1  
47         IF(X1.LT.CKIN(21).OR.X1.GT.CKIN(22)) MINT(51)=1 
48         IF(X2.LT.CKIN(23).OR.X2.GT.CKIN(24)) MINT(51)=1 
49         IF(XF.LT.CKIN(25).OR.XF.GT.CKIN(26)) MINT(51)=1 
50         IF(YST.LT.CKIN(7).OR.YST.GT.CKIN(8)) MINT(51)=1 
51         IF(ISET(ISUB).EQ.2.OR.ISET(ISUB).EQ.4) THEN 
52           PTH=0.5*BE34*SQRT(TAU*VINT(2)*(1.-CTH**2))    
53           Y3=YST+0.5*LOG((1.+RM3-RM4+BE34*CTH)/(1.+RM3-RM4-BE34*CTH))   
54           Y4=YST+0.5*LOG((1.+RM4-RM3-BE34*CTH)/(1.+RM4-RM3+BE34*CTH))   
55           YLARGE=MAX(Y3,Y4) 
56           YSMALL=MIN(Y3,Y4) 
57           ETALAR=10.    
58           ETASMA=-10.   
59           STH=SQRT(1.-CTH**2)   
60           IF(STH.LT.1.E-6) GOTO 100 
61           EXPET3=((1.+RM3-RM4)*SINH(YST)+BE34*COSH(YST)*CTH+    
62      &    SQRT(((1.+RM3-RM4)*COSH(YST)+BE34*SINH(YST)*CTH)**2-4.*RM3))/ 
63      &    (BE34*STH)    
64           EXPET4=((1.-RM3+RM4)*SINH(YST)-BE34*COSH(YST)*CTH+    
65      &    SQRT(((1.-RM3+RM4)*COSH(YST)-BE34*SINH(YST)*CTH)**2-4.*RM4))/ 
66      &    (BE34*STH)    
67           ETA3=LOG(MIN(1.E10,MAX(1.E-10,EXPET3)))   
68           ETA4=LOG(MIN(1.E10,MAX(1.E-10,EXPET4)))   
69           ETALAR=MAX(ETA3,ETA4) 
70           ETASMA=MIN(ETA3,ETA4) 
71   100     CTS3=((1.+RM3-RM4)*SINH(YST)+BE34*COSH(YST)*CTH)/ 
72      &    SQRT(((1.+RM3-RM4)*COSH(YST)+BE34*SINH(YST)*CTH)**2-4.*RM3)   
73           CTS4=((1.-RM3+RM4)*SINH(YST)-BE34*COSH(YST)*CTH)/ 
74      &    SQRT(((1.-RM3+RM4)*COSH(YST)-BE34*SINH(YST)*CTH)**2-4.*RM4)   
75           CTSLAR=MAX(CTS3,CTS4) 
76           CTSSMA=MIN(CTS3,CTS4) 
77           IF(PTH.LT.PTHMIN) MINT(51)=1  
78           IF(CKIN(4).GE.0..AND.PTH.GT.CKIN(4)) MINT(51)=1   
79           IF(YLARGE.LT.CKIN(9).OR.YLARGE.GT.CKIN(10)) MINT(51)=1    
80           IF(YSMALL.LT.CKIN(11).OR.YSMALL.GT.CKIN(12)) MINT(51)=1   
81           IF(ETALAR.LT.CKIN(13).OR.ETALAR.GT.CKIN(14)) MINT(51)=1   
82           IF(ETASMA.LT.CKIN(15).OR.ETASMA.GT.CKIN(16)) MINT(51)=1   
83           IF(CTSLAR.LT.CKIN(17).OR.CTSLAR.GT.CKIN(18)) MINT(51)=1   
84           IF(CTSSMA.LT.CKIN(19).OR.CTSSMA.GT.CKIN(20)) MINT(51)=1   
85           IF(CTH.LT.CKIN(27).OR.CTH.GT.CKIN(28)) MINT(51)=1 
86         ENDIF   
87         IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) THEN 
88           IF(TAUP*VINT(2).LT.CKIN(31)**2) MINT(51)=1    
89           IF(CKIN(32).GE.0..AND.TAUP*VINT(2).GT.CKIN(32)**2) MINT(51)=1 
90         ENDIF   
91     
92       ELSEIF(ILIM.EQ.1) THEN    
93 C...Calculate limits on tau 
94 C...0) due to definition    
95         TAUMN0=0.   
96         TAUMX0=1.   
97 C...1) due to limits on subsystem mass  
98         TAUMN1=CKIN(1)**2/VINT(2)   
99         TAUMX1=1.   
100         IF(CKIN(2).GE.0.) TAUMX1=CKIN(2)**2/VINT(2) 
101 C...2) due to limits on pT-hat (and non-overlapping rapidity intervals) 
102         TM3=SQRT(SQM3+PTHMIN**2)    
103         TM4=SQRT(SQM4+PTHMIN**2)    
104         YDCOSH=1.   
105         IF(CKIN(9).GT.CKIN(12)) YDCOSH=COSH(CKIN(9)-CKIN(12))   
106         TAUMN2=(TM3**2+2.*TM3*TM4*YDCOSH+TM4**2)/VINT(2)    
107         TAUMX2=1.   
108 C...3) due to limits on pT-hat and cos(theta-hat)   
109         CTH2MN=MIN(CKIN(27)**2,CKIN(28)**2) 
110         CTH2MX=MAX(CKIN(27)**2,CKIN(28)**2) 
111         TAUMN3=0.   
112         IF(CKIN(27)*CKIN(28).GT.0.) TAUMN3= 
113      &  (SQRT(SQM3+PTHMIN**2/(1.-CTH2MN))+  
114      &  SQRT(SQM4+PTHMIN**2/(1.-CTH2MN)))**2/VINT(2)    
115         TAUMX3=1.   
116         IF(CKIN(4).GE.0..AND.CTH2MX.LT.1.) TAUMX3=  
117      &  (SQRT(SQM3+CKIN(4)**2/(1.-CTH2MX))+ 
118      &  SQRT(SQM4+CKIN(4)**2/(1.-CTH2MX)))**2/VINT(2)   
119 C...4) due to limits on x1 and x2   
120         TAUMN4=CKIN(21)*CKIN(23)    
121         TAUMX4=CKIN(22)*CKIN(24)    
122 C...5) due to limits on xF  
123         TAUMN5=0.   
124         TAUMX5=MAX(1.-CKIN(25),1.+CKIN(26)) 
125         VINT(11)=MAX(TAUMN0,TAUMN1,TAUMN2,TAUMN3,TAUMN4,TAUMN5) 
126         VINT(31)=MIN(TAUMX0,TAUMX1,TAUMX2,TAUMX3,TAUMX4,TAUMX5) 
127         IF(MINT(43).EQ.1.AND.(ISET(ISUB).EQ.1.OR.ISET(ISUB).EQ.2)) THEN 
128           VINT(11)=0.99999  
129           VINT(31)=1.00001  
130         ENDIF   
131         IF(VINT(31).LE.VINT(11)) MINT(51)=1 
132     
133       ELSEIF(ILIM.EQ.2) THEN    
134 C...Calculate limits on y*  
135         IF(ISET(ISUB).EQ.3.OR.ISET(ISUB).EQ.4) TAU=VINT(26) 
136         TAURT=SQRT(TAU) 
137 C...0) due to kinematics    
138         YSTMN0=LOG(TAURT)   
139         YSTMX0=-YSTMN0  
140 C...1) due to explicit limits   
141         YSTMN1=CKIN(7)  
142         YSTMX1=CKIN(8)  
143 C...2) due to limits on x1  
144         YSTMN2=LOG(MAX(TAU,CKIN(21))/TAURT) 
145         YSTMX2=LOG(MAX(TAU,CKIN(22))/TAURT) 
146 C...3) due to limits on x2  
147         YSTMN3=-LOG(MAX(TAU,CKIN(24))/TAURT)    
148         YSTMX3=-LOG(MAX(TAU,CKIN(23))/TAURT)    
149 C...4) due to limits on xF  
150         YEPMN4=0.5*ABS(CKIN(25))/TAURT  
151         YSTMN4=SIGN(LOG(SQRT(1.+YEPMN4**2)+YEPMN4),CKIN(25))    
152         YEPMX4=0.5*ABS(CKIN(26))/TAURT  
153         YSTMX4=SIGN(LOG(SQRT(1.+YEPMX4**2)+YEPMX4),CKIN(26))    
154 C...5) due to simultaneous limits on y-large and y-small    
155         YEPSMN=(RM3-RM4)*SINH(CKIN(9)-CKIN(11)) 
156         YEPSMX=(RM3-RM4)*SINH(CKIN(10)-CKIN(12))    
157         YDIFMN=ABS(LOG(SQRT(1.+YEPSMN**2)-YEPSMN))  
158         YDIFMX=ABS(LOG(SQRT(1.+YEPSMX**2)-YEPSMX))  
159         YSTMN5=0.5*(CKIN(9)+CKIN(11)-YDIFMN)    
160         YSTMX5=0.5*(CKIN(10)+CKIN(12)+YDIFMX)   
161 C...6) due to simultaneous limits on cos(theta-hat) and y-large or  
162 C...   y-small  
163         CTHLIM=1.-4.*PTHMIN**2/(BE34*TAU*VINT(2))
164         if (CTHLIM .lt. 0) cthlim = 0.
165         CTHLIM=SQRT(cthlim) 
166         RZMN=BE34*MAX(CKIN(27),-CTHLIM) 
167         RZMX=BE34*MIN(CKIN(28),CTHLIM)  
168         YEX3MX=(1.+RM3-RM4+RZMX)/MAX(1E-10,1.+RM3-RM4-RZMX) 
169         YEX4MX=(1.+RM4-RM3-RZMN)/MAX(1E-10,1.+RM4-RM3+RZMN) 
170         YEX3MN=MAX(1E-10,1.+RM3-RM4+RZMN)/(1.+RM3-RM4-RZMN) 
171         YEX4MN=MAX(1E-10,1.+RM4-RM3-RZMX)/(1.+RM4-RM3+RZMX) 
172         YSTMN6=CKIN(9)-0.5*LOG(MAX(YEX3MX,YEX4MX))  
173         YSTMX6=CKIN(12)-0.5*LOG(MIN(YEX3MN,YEX4MN)) 
174         VINT(12)=MAX(YSTMN0,YSTMN1,YSTMN2,YSTMN3,YSTMN4,YSTMN5,YSTMN6)  
175         VINT(32)=MIN(YSTMX0,YSTMX1,YSTMX2,YSTMX3,YSTMX4,YSTMX5,YSTMX6)  
176         IF(MINT(43).EQ.1) THEN  
177           VINT(12)=-0.00001 
178           VINT(32)=0.00001  
179         ELSEIF(MINT(43).EQ.2) THEN  
180           VINT(12)=0.99999*YSTMX0   
181           VINT(32)=1.00001*YSTMX0   
182         ELSEIF(MINT(43).EQ.3) THEN  
183           VINT(12)=-1.00001*YSTMX0  
184           VINT(32)=-0.99999*YSTMX0  
185         ENDIF   
186         IF(VINT(32).LE.VINT(12)) MINT(51)=1 
187     
188       ELSEIF(ILIM.EQ.3) THEN    
189 C...Calculate limits on cos(theta-hat)  
190         YST=VINT(22)    
191 C...0) due to definition    
192         CTNMN0=-1.  
193         CTNMX0=0.   
194         CTPMN0=0.   
195         CTPMX0=1.   
196 C...1) due to explicit limits   
197         CTNMN1=MIN(0.,CKIN(27)) 
198         CTNMX1=MIN(0.,CKIN(28)) 
199         CTPMN1=MAX(0.,CKIN(27)) 
200         CTPMX1=MAX(0.,CKIN(28)) 
201 C...2) due to limits on pT-hat  
202         CTNMN2=1.-4.*PTHMIN**2/(BE34**2*TAU*VINT(2))
203         if (ctnmn2 .le. 0) ctnmn2 = 0
204         CTNMN2=-SQRT(ctnmn2) 
205         CTPMX2=-CTNMN2  
206         CTNMX2=0.   
207         CTPMN2=0.   
208         IF(CKIN(4).GE.0.) THEN  
209           CTNMX2=-SQRT(MAX(0.,1.-4.*CKIN(4)**2/(BE34**2*TAU*VINT(2))))  
210           CTPMN2=-CTNMX2    
211         ENDIF   
212 C...3) due to limits on y-large and y-small 
213         CTNMN3=MIN(0.,MAX((1.+RM3-RM4)/BE34*TANH(CKIN(11)-YST), 
214      &  -(1.-RM3+RM4)/BE34*TANH(CKIN(10)-YST))) 
215         CTNMX3=MIN(0.,(1.+RM3-RM4)/BE34*TANH(CKIN(12)-YST), 
216      &  -(1.-RM3+RM4)/BE34*TANH(CKIN(9)-YST))   
217         CTPMN3=MAX(0.,(1.+RM3-RM4)/BE34*TANH(CKIN(9)-YST),  
218      &  -(1.-RM3+RM4)/BE34*TANH(CKIN(12)-YST))  
219         CTPMX3=MAX(0.,MIN((1.+RM3-RM4)/BE34*TANH(CKIN(10)-YST), 
220      &  -(1.-RM3+RM4)/BE34*TANH(CKIN(11)-YST))) 
221         VINT(13)=MAX(CTNMN0,CTNMN1,CTNMN2,CTNMN3)   
222         VINT(33)=MIN(CTNMX0,CTNMX1,CTNMX2,CTNMX3)   
223         VINT(14)=MAX(CTPMN0,CTPMN1,CTPMN2,CTPMN3)   
224         VINT(34)=MIN(CTPMX0,CTPMX1,CTPMX2,CTPMX3)   
225         IF(VINT(33).LE.VINT(13).AND.VINT(34).LE.VINT(14)) MINT(51)=1    
226     
227       ELSEIF(ILIM.EQ.4) THEN    
228 C...Calculate limits on tau'    
229 C...0) due to kinematics    
230         TAPMN0=TAU  
231         TAPMX0=1.   
232 C...1) due to explicit limits   
233         TAPMN1=CKIN(31)**2/VINT(2)  
234         TAPMX1=1.   
235         IF(CKIN(32).GE.0.) TAPMX1=CKIN(32)**2/VINT(2)   
236         VINT(16)=MAX(TAPMN0,TAPMN1) 
237         VINT(36)=MIN(TAPMX0,TAPMX1) 
238         IF(MINT(43).EQ.1) THEN  
239           VINT(16)=0.99999  
240           VINT(36)=1.00001  
241         ENDIF   
242         IF(VINT(36).LE.VINT(16)) MINT(51)=1 
243     
244       ENDIF 
245       RETURN    
246     
247 C...Special case for low-pT and multiple interactions:  
248 C...effective kinematical limits for tau, y*, cos(theta-hat).   
249   110 IF(ILIM.EQ.0) THEN    
250       ELSEIF(ILIM.EQ.1) THEN    
251         IF(MSTP(82).LE.1) VINT(11)=4.*PARP(81)**2/VINT(2)   
252         IF(MSTP(82).GE.2) VINT(11)=PARP(82)**2/VINT(2)  
253         VINT(31)=1. 
254       ELSEIF(ILIM.EQ.2) THEN    
255         VINT(12)=0.5*LOG(VINT(21))  
256         VINT(32)=-VINT(12)  
257       ELSEIF(ILIM.EQ.3) THEN    
258         IF(MSTP(82).LE.1) ST2EFF=4.*PARP(81)**2/(VINT(21)*VINT(2))  
259         IF(MSTP(82).GE.2) ST2EFF=0.01*PARP(82)**2/(VINT(21)*VINT(2))    
260         VINT(13)=-SQRT(MAX(0.,1.-ST2EFF))   
261         VINT(33)=0. 
262         VINT(14)=0. 
263         VINT(34)=-VINT(13)  
264       ENDIF 
265     
266       RETURN    
267       END