Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyklim.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyklim.f
1 
2 C***********************************************************************
3 
4 C...PYKLIM
5 C...Checks generated variables against pre-set kinematical limits;
6 C...also calculates limits on variables used in generation.
7 
8  SUBROUTINE pyklim(ILIM)
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...Commonblocks.
15  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
17  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
19  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
20  common/pypars/mstp(200),parp(200),msti(200),pari(200)
21  common/pyint1/mint(400),vint(400)
22  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
23  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/,/pysubs/,/pypars/,
24  &/pyint1/,/pyint2/
25 
26 C...Common kinematical expressions.
27  mint(51)=0
28  isub=mint(1)
29  istsb=iset(isub)
30  IF(isub.EQ.96) goto 100
31  sqm3=vint(63)
32  sqm4=vint(64)
33  IF(ilim.NE.0) THEN
34  IF(abs(sqm3).LT.1d-4.AND.abs(sqm4).LT.1d-4) THEN
35  ckin09=max(ckin(9),ckin(13))
36  ckin10=min(ckin(10),ckin(14))
37  ckin11=max(ckin(11),ckin(15))
38  ckin12=min(ckin(12),ckin(16))
39  ELSE
40  ckin09=max(ckin(9),min(0d0,ckin(13)))
41  ckin10=min(ckin(10),max(0d0,ckin(14)))
42  ckin11=max(ckin(11),min(0d0,ckin(15)))
43  ckin12=min(ckin(12),max(0d0,ckin(16)))
44  ENDIF
45  ENDIF
46  IF(ilim.NE.1) THEN
47  tau=vint(21)
48  rm3=sqm3/(tau*vint(2))
49  rm4=sqm4/(tau*vint(2))
50  be34=sqrt(max(1d-20,(1d0-rm3-rm4)**2-4d0*rm3*rm4))
51  ENDIF
52  pthmin=ckin(3)
53  IF(min(sqm3,sqm4).LT.ckin(6)**2.AND.istsb.NE.1.AND.istsb.NE.3)
54  &pthmin=max(ckin(3),ckin(5))
55 
56  IF(ilim.EQ.0) THEN
57 C...Check generated values of tau, y*, cos(theta-hat), and tau' against
58 C...pre-set kinematical limits.
59  yst=vint(22)
60  cth=vint(23)
61  taup=vint(26)
62  taue=tau
63  IF(istsb.GE.3.AND.istsb.LE.5) taue=taup
64  x1=sqrt(taue)*exp(yst)
65  x2=sqrt(taue)*exp(-yst)
66  xf=x1-x2
67  IF(mint(47).NE.1) THEN
68  IF(tau*vint(2).LT.ckin(1)**2) mint(51)=1
69  IF(ckin(2).GE.0d0.AND.tau*vint(2).GT.ckin(2)**2) mint(51)=1
70  IF(yst.LT.ckin(7).OR.yst.GT.ckin(8)) mint(51)=1
71  IF(xf.LT.ckin(25).OR.xf.GT.ckin(26)) mint(51)=1
72  ENDIF
73  IF(mint(45).NE.1) THEN
74  IF(x1.LT.ckin(21).OR.x1.GT.ckin(22)) mint(51)=1
75  ENDIF
76  IF(mint(46).NE.1) THEN
77  IF(x2.LT.ckin(23).OR.x2.GT.ckin(24)) mint(51)=1
78  ENDIF
79  IF(mint(45).EQ.2) THEN
80  IF(x1.GT.1d0-2d0*parp(111)/vint(1)) mint(51)=1
81  ENDIF
82  IF(mint(46).EQ.2) THEN
83  IF(x2.GT.1d0-2d0*parp(111)/vint(1)) mint(51)=1
84  ENDIF
85  IF(istsb.EQ.2.OR.istsb.EQ.4) THEN
86  pth=0.5d0*be34*sqrt(tau*vint(2)*max(0d0,1d0-cth**2))
87  expy3=max(1d-20,(1d0+rm3-rm4+be34*cth)/
88  & max(1d-20,(1d0+rm3-rm4-be34*cth)))
89  expy4=max(1d-20,(1d0-rm3+rm4-be34*cth)/
90  & max(1d-20,(1d0-rm3+rm4+be34*cth)))
91  y3=yst+0.5d0*log(expy3)
92  y4=yst+0.5d0*log(expy4)
93  ylarge=max(y3,y4)
94  ysmall=min(y3,y4)
95  etalar=20d0
96  etasma=-20d0
97  sth=sqrt(max(0d0,1d0-cth**2))
98  exsq3=sqrt(max(1d-20,((1d0+rm3-rm4)*cosh(yst)+be34*sinh(yst)*
99  & cth)**2-4d0*rm3))
100  exsq4=sqrt(max(1d-20,((1d0-rm3+rm4)*cosh(yst)-be34*sinh(yst)*
101  & cth)**2-4d0*rm4))
102  IF(sth.GE.1d-10) THEN
103  expet3=((1d0+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth+exsq3)/
104  & (be34*sth)
105  expet4=((1d0-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth+exsq4)/
106  & (be34*sth)
107  eta3=log(min(1d10,max(1d-10,expet3)))
108  eta4=log(min(1d10,max(1d-10,expet4)))
109  etalar=max(eta3,eta4)
110  etasma=min(eta3,eta4)
111  ENDIF
112  cts3=((1d0+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth)/exsq3
113  cts4=((1d0-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth)/exsq4
114  ctslar=min(1d0,max(-1d0,cts3,cts4))
115  ctssma=max(-1d0,min(1d0,cts3,cts4))
116  sh=tau*vint(2)
117  rpts=4d0*vint(71)**2/sh
118  be34l=sqrt(max(0d0,(1d0-rm3-rm4)**2-4d0*rm3*rm4-rpts))
119  rm34=max(1d-20,2d0*rm3*rm4)
120  IF(2d0*vint(71)**2/(vint(21)*vint(2)).LT.0.0001d0)
121  & rm34=max(rm34,2d0*vint(71)**2/(vint(21)*vint(2)))
122  rthm=(4d0*rm3*rm4+rpts)/(1d0-rm3-rm4+be34l)
123  tha=0.5d0*sh*max(rthm,1d0-rm3-rm4-be34*cth)
124  uha=0.5d0*sh*max(rthm,1d0-rm3-rm4+be34*cth)
125  IF(pth.LT.pthmin) mint(51)=1
126  IF(ckin(4).GE.0d0.AND.pth.GT.ckin(4)) mint(51)=1
127  IF(ylarge.LT.ckin(9).OR.ylarge.GT.ckin(10)) mint(51)=1
128  IF(ysmall.LT.ckin(11).OR.ysmall.GT.ckin(12)) mint(51)=1
129  IF(etalar.LT.ckin(13).OR.etalar.GT.ckin(14)) mint(51)=1
130  IF(etasma.LT.ckin(15).OR.etasma.GT.ckin(16)) mint(51)=1
131  IF(ctslar.LT.ckin(17).OR.ctslar.GT.ckin(18)) mint(51)=1
132  IF(ctssma.LT.ckin(19).OR.ctssma.GT.ckin(20)) mint(51)=1
133  IF(cth.LT.ckin(27).OR.cth.GT.ckin(28)) mint(51)=1
134  IF(tha.LT.ckin(35)) mint(51)=1
135  IF(ckin(36).GE.0d0.AND.tha.GT.ckin(36)) mint(51)=1
136  IF(uha.LT.ckin(37)) mint(51)=1
137  IF(ckin(38).GE.0d0.AND.uha.GT.ckin(38)) mint(51)=1
138  ENDIF
139  IF(istsb.GE.3.AND.istsb.LE.5) THEN
140  IF(taup*vint(2).LT.ckin(31)**2) mint(51)=1
141  IF(ckin(32).GE.0d0.AND.taup*vint(2).GT.ckin(32)**2) mint(51)=1
142  ENDIF
143 
144 C...Additional cuts on W2 (approximately) in DIS.
145  IF(isub.EQ.10.AND.mint(43).GE.2) THEN
146  xbj=x2
147  IF(iabs(mint(12)).LT.20) xbj=x1
148  q2bj=tha
149  w2bj=q2bj*(1d0-xbj)/xbj
150  IF(w2bj.LT.ckin(39)) mint(51)=1
151  IF(ckin(40).GT.0d0.AND.w2bj.GT.ckin(40)) mint(51)=1
152  ENDIF
153 
154  ELSEIF(ilim.EQ.1) THEN
155 C...Calculate limits on tau
156 C...0) due to definition
157  taumn0=0d0
158  taumx0=1d0
159 C...1) due to limits on subsystem mass
160  taumn1=ckin(1)**2/vint(2)
161  taumx1=1d0
162  IF(ckin(2).GE.0d0) taumx1=ckin(2)**2/vint(2)
163 C...2) due to limits on pT-hat (and non-overlapping rapidity intervals)
164  tm3=sqrt(sqm3+pthmin**2)
165  tm4=sqrt(sqm4+pthmin**2)
166  ydcosh=1d0
167  IF(ckin09.GT.ckin12) ydcosh=cosh(ckin09-ckin12)
168  taumn2=(tm3**2+2d0*tm3*tm4*ydcosh+tm4**2)/vint(2)
169  taumx2=1d0
170 C...3) due to limits on pT-hat and cos(theta-hat)
171  cth2mn=min(ckin(27)**2,ckin(28)**2)
172  cth2mx=max(ckin(27)**2,ckin(28)**2)
173  taumn3=0d0
174  IF(ckin(27)*ckin(28).GT.0d0) taumn3=
175  & (sqrt(sqm3+pthmin**2/(1d0-cth2mn))+
176  & sqrt(sqm4+pthmin**2/(1d0-cth2mn)))**2/vint(2)
177  taumx3=1d0
178  IF(ckin(4).GE.0d0.AND.cth2mx.LT.1d0) taumx3=
179  & (sqrt(sqm3+ckin(4)**2/(1d0-cth2mx))+
180  & sqrt(sqm4+ckin(4)**2/(1d0-cth2mx)))**2/vint(2)
181 C...4) due to limits on x1 and x2
182  taumn4=ckin(21)*ckin(23)
183  taumx4=ckin(22)*ckin(24)
184 C...5) due to limits on xF
185  taumn5=0d0
186  taumx5=max(1d0-ckin(25),1d0+ckin(26))
187 C...6) due to limits on that and uhat
188  taumn6=(sqm3+sqm4+ckin(35)+ckin(37))/vint(2)
189  taumx6=1d0
190  IF(ckin(36).GT.0d0.AND.ckin(38).GT.0d0) taumx6=
191  & (sqm3+sqm4+ckin(36)+ckin(38))/vint(2)
192 
193 C...Net effect of all separate limits.
194  vint(11)=max(taumn0,taumn1,taumn2,taumn3,taumn4,taumn5,taumn6)
195  vint(31)=min(taumx0,taumx1,taumx2,taumx3,taumx4,taumx5,taumx6)
196  IF(mint(47).EQ.1.AND.(istsb.EQ.1.OR.istsb.EQ.2)) THEN
197  vint(11)=1d0-1d-9
198  vint(31)=1d0+1d-9
199  ELSEIF(mint(47).EQ.5) THEN
200  vint(31)=min(vint(31),1d0-2d-10)
201  ELSEIF(mint(47).GE.6) THEN
202  vint(31)=min(vint(31),1d0-1d-10)
203  ENDIF
204  IF(vint(31).LE.vint(11)) mint(51)=1
205 
206  ELSEIF(ilim.EQ.2) THEN
207 C...Calculate limits on y*
208  taue=tau
209  IF(istsb.GE.3.AND.istsb.LE.5) taue=vint(26)
210  taurt=sqrt(taue)
211 C...0) due to kinematics
212  ystmn0=log(taurt)
213  ystmx0=-ystmn0
214 C...1) due to explicit limits
215  ystmn1=ckin(7)
216  ystmx1=ckin(8)
217 C...2) due to limits on x1
218  ystmn2=log(max(taue,ckin(21))/taurt)
219  ystmx2=log(max(taue,ckin(22))/taurt)
220 C...3) due to limits on x2
221  ystmn3=-log(max(taue,ckin(24))/taurt)
222  ystmx3=-log(max(taue,ckin(23))/taurt)
223 C...4) due to limits on xF
224  yepmn4=0.5d0*abs(ckin(25))/taurt
225  ystmn4=sign(log(max(1d-20,sqrt(1d0+yepmn4**2)+yepmn4)),ckin(25))
226  yepmx4=0.5d0*abs(ckin(26))/taurt
227  ystmx4=sign(log(max(1d-20,sqrt(1d0+yepmx4**2)+yepmx4)),ckin(26))
228 C...5) due to simultaneous limits on y-large and y-small
229  yepsmn=(rm3-rm4)*sinh(ckin09-ckin11)
230  yepsmx=(rm3-rm4)*sinh(ckin10-ckin12)
231  ydifmn=abs(log(max(1d-20,sqrt(1d0+yepsmn**2)-yepsmn)))
232  ydifmx=abs(log(max(1d-20,sqrt(1d0+yepsmx**2)-yepsmx)))
233  ystmn5=0.5d0*(ckin09+ckin11-ydifmn)
234  ystmx5=0.5d0*(ckin10+ckin12+ydifmx)
235 C...6) due to simultaneous limits on cos(theta-hat) and y-large or
236 C... y-small
237  cthlim=sqrt(max(0d0,1d0-4d0*pthmin**2/(be34**2*taue*vint(2))))
238  rzmn=be34*max(ckin(27),-cthlim)
239  rzmx=be34*min(ckin(28),cthlim)
240  yex3mx=(1d0+rm3-rm4+rzmx)/max(1d-10,1d0+rm3-rm4-rzmx)
241  yex4mx=(1d0+rm4-rm3-rzmn)/max(1d-10,1d0+rm4-rm3+rzmn)
242  yex3mn=max(1d-10,1d0+rm3-rm4+rzmn)/(1d0+rm3-rm4-rzmn)
243  yex4mn=max(1d-10,1d0+rm4-rm3-rzmx)/(1d0+rm4-rm3+rzmx)
244  ystmn6=ckin09-0.5d0*log(max(yex3mx,yex4mx))
245  ystmx6=ckin12-0.5d0*log(min(yex3mn,yex4mn))
246 
247 C...Net effect of all separate limits.
248  vint(12)=max(ystmn0,ystmn1,ystmn2,ystmn3,ystmn4,ystmn5,ystmn6)
249  vint(32)=min(ystmx0,ystmx1,ystmx2,ystmx3,ystmx4,ystmx5,ystmx6)
250  IF(mint(47).EQ.1) THEN
251  vint(12)=-1d-9
252  vint(32)=1d-9
253  ELSEIF(mint(47).EQ.2.OR.mint(47).EQ.6) THEN
254  vint(12)=(1d0-1d-9)*ystmx0
255  vint(32)=(1d0+1d-9)*ystmx0
256  ELSEIF(mint(47).EQ.3.OR.mint(47).EQ.7) THEN
257  vint(12)=-(1d0+1d-9)*ystmx0
258  vint(32)=-(1d0-1d-9)*ystmx0
259  ELSEIF(mint(47).EQ.5) THEN
260  ystee=log((1d0-1d-10)/taurt)
261  vint(12)=max(vint(12),-ystee)
262  vint(32)=min(vint(32),ystee)
263  ENDIF
264  IF(vint(32).LE.vint(12)) mint(51)=1
265 
266  ELSEIF(ilim.EQ.3) THEN
267 C...Calculate limits on cos(theta-hat)
268  yst=vint(22)
269 C...0) due to definition
270  ctnmn0=-1d0
271  ctnmx0=0d0
272  ctpmn0=0d0
273  ctpmx0=1d0
274 C...1) due to explicit limits
275  ctnmn1=min(0d0,ckin(27))
276  ctnmx1=min(0d0,ckin(28))
277  ctpmn1=max(0d0,ckin(27))
278  ctpmx1=max(0d0,ckin(28))
279 C...2) due to limits on pT-hat
280  ctnmn2=-sqrt(max(0d0,1d0-4d0*pthmin**2/(be34**2*tau*vint(2))))
281  ctpmx2=-ctnmn2
282  ctnmx2=0d0
283  ctpmn2=0d0
284  IF(ckin(4).GE.0d0) THEN
285  ctnmx2=-sqrt(max(0d0,1d0-4d0*ckin(4)**2/
286  & (be34**2*tau*vint(2))))
287  ctpmn2=-ctnmx2
288  ENDIF
289 C...3) due to limits on y-large and y-small
290  ctnmn3=min(0d0,max((1d0+rm3-rm4)/be34*tanh(ckin11-yst),
291  & -(1d0-rm3+rm4)/be34*tanh(ckin10-yst)))
292  ctnmx3=min(0d0,(1d0+rm3-rm4)/be34*tanh(ckin12-yst),
293  & -(1d0-rm3+rm4)/be34*tanh(ckin09-yst))
294  ctpmn3=max(0d0,(1d0+rm3-rm4)/be34*tanh(ckin09-yst),
295  & -(1d0-rm3+rm4)/be34*tanh(ckin12-yst))
296  ctpmx3=max(0d0,min((1d0+rm3-rm4)/be34*tanh(ckin10-yst),
297  & -(1d0-rm3+rm4)/be34*tanh(ckin11-yst)))
298 C...4) due to limits on that
299  ctnmn4=-1d0
300  ctnmx4=0d0
301  ctpmn4=0d0
302  ctpmx4=1d0
303  sh=tau*vint(2)
304  IF(ckin(35).GT.0d0) THEN
305  ctlim=(1d0-rm3-rm4-2d0*ckin(35)/sh)/be34
306  IF(ctlim.GT.0d0) THEN
307  ctpmx4=ctlim
308  ELSE
309  ctpmx4=0d0
310  ctnmx4=ctlim
311  ENDIF
312  ENDIF
313  IF(ckin(36).GT.0d0) THEN
314  ctlim=(1d0-rm3-rm4-2d0*ckin(36)/sh)/be34
315  IF(ctlim.LT.0d0) THEN
316  ctnmn4=ctlim
317  ELSE
318  ctnmn4=0d0
319  ctpmn4=ctlim
320  ENDIF
321  ENDIF
322 C...5) due to limits on uhat
323  ctnmn5=-1d0
324  ctnmx5=0d0
325  ctpmn5=0d0
326  ctpmx5=1d0
327  IF(ckin(37).GT.0d0) THEN
328  ctlim=(2d0*ckin(37)/sh-(1d0-rm3-rm4))/be34
329  IF(ctlim.LT.0d0) THEN
330  ctnmn5=ctlim
331  ELSE
332  ctnmn5=0d0
333  ctpmn5=ctlim
334  ENDIF
335  ENDIF
336  IF(ckin(38).GT.0d0) THEN
337  ctlim=(2d0*ckin(38)/sh-(1d0-rm3-rm4))/be34
338  IF(ctlim.GT.0d0) THEN
339  ctpmx5=ctlim
340  ELSE
341  ctpmx5=0d0
342  ctnmx5=ctlim
343  ENDIF
344  ENDIF
345 
346 C...Net effect of all separate limits.
347  vint(13)=max(ctnmn0,ctnmn1,ctnmn2,ctnmn3,ctnmn4,ctnmn5)
348  vint(33)=min(ctnmx0,ctnmx1,ctnmx2,ctnmx3,ctnmx4,ctnmx5)
349  vint(14)=max(ctpmn0,ctpmn1,ctpmn2,ctpmn3,ctpmn4,ctpmn5)
350  vint(34)=min(ctpmx0,ctpmx1,ctpmx2,ctpmx3,ctpmx4,ctpmx5)
351  IF(vint(33).LE.vint(13).AND.vint(34).LE.vint(14)) mint(51)=1
352 
353  IF(vint(14).GT.vint(34)) vint(34)=vint(14)
354  IF(vint(13).GT.vint(33)) vint(33)=vint(13)
355 
356  ELSEIF(ilim.EQ.4) THEN
357 C...Calculate limits on tau'
358 C...0) due to kinematics
359  tapmn0=tau
360  IF(istsb.EQ.5.AND.vint(201).GT.0d0) THEN
361  pqrat=(vint(201)+vint(206))/vint(1)
362  tapmn0=(sqrt(tau)+pqrat)**2
363  ENDIF
364  tapmx0=1d0
365 C...1) due to explicit limits
366  tapmn1=ckin(31)**2/vint(2)
367  tapmx1=1d0
368  IF(ckin(32).GE.0d0) tapmx1=ckin(32)**2/vint(2)
369 
370 C...Net effect of all separate limits.
371  vint(16)=max(tapmn0,tapmn1)
372  vint(36)=min(tapmx0,tapmx1)
373  IF(mint(47).EQ.1) THEN
374  vint(16)=1d0-1d-9
375  vint(36)=1d0+1d-9
376  ELSEIF(mint(47).EQ.5) THEN
377  vint(36)=min(vint(36),1d0-2d-10)
378  ELSEIF(mint(47).EQ.6.OR.mint(47).EQ.7) THEN
379  vint(36)=min(vint(36),1d0-1d-10)
380  ENDIF
381  IF(vint(36).LE.vint(16)) mint(51)=1
382 
383  ENDIF
384  RETURN
385 
386 C...Special case for low-pT and multiple interactions:
387 C...effective kinematical limits for tau, y*, cos(theta-hat).
388  100 IF(ilim.EQ.0) THEN
389  ELSEIF(ilim.EQ.1) THEN
390  IF(mstp(82).LE.1) THEN
391  vint(11)=4d0*(parp(81)*(vint(1)/parp(89))**parp(90))**2/
392  & vint(2)
393  ELSE
394  vint(11)=(parp(82)*(vint(1)/parp(89))**parp(90))**2/vint(2)
395  ENDIF
396  vint(31)=1d0
397  ELSEIF(ilim.EQ.2) THEN
398  vint(12)=0.5d0*log(vint(21))
399  vint(32)=-vint(12)
400  ELSEIF(ilim.EQ.3) THEN
401  IF(mstp(82).LE.1) THEN
402  st2eff=4d0*(parp(81)*(vint(1)/parp(89))**parp(90))**2/
403  & (vint(21)*vint(2))
404  ELSE
405  st2eff=0.01d0*(parp(82)*(vint(1)/parp(89))**parp(90))**2/
406  & (vint(21)*vint(2))
407  ENDIF
408  vint(13)=-sqrt(max(0d0,1d0-st2eff))
409  vint(33)=0d0
410  vint(14)=0d0
411  vint(34)=-vint(13)
412  ENDIF
413 
414  RETURN
415  END