ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhiklim.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiklim.f
1 
2 C***********************************************************************
3 
4  SUBROUTINE pyhiklim(ILIM)
5 
6 C...Checks generated variables against pre-set kinematical limits;
7 C...also calculates limits on variables used in generation.
8  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
9  SAVE /ludat1/
10  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
11  SAVE /ludat2/
12  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
13  SAVE /ludat3/
14  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
15  SAVE /pyhipars/
16  common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200)
17  SAVE /pyhisubs/
18  common/pyhiint1/mint(400),vint(400)
19  SAVE /pyhiint1/
20  common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2)
21  SAVE /pyhiint2/
22 
23 C...Common kinematical expressions.
24  isub=mint(1)
25  IF(isub.EQ.96) goto 110
26  sqm3=vint(63)
27  sqm4=vint(64)
28  IF(ilim.NE.1) THEN
29  tau=vint(21)
30  rm3=sqm3/(tau*vint(2))
31  rm4=sqm4/(tau*vint(2))
32  be34=sqrt((1.-rm3-rm4)**2-4.*rm3*rm4)
33  ENDIF
34  pthmin=ckin(3)
35  IF(min(sqm3,sqm4).LT.ckin(6)**2) pthmin=max(ckin(3),ckin(5))
36 
37  IF(ilim.EQ.0) THEN
38 C...Check generated values of tau, y*, cos(theta-hat), and tau' against
39 C...pre-set kinematical limits.
40  yst=vint(22)
41  cth=vint(23)
42  taup=vint(26)
43  IF(iset(isub).LE.2) THEN
44  x1=sqrt(tau)*exp(yst)
45  x2=sqrt(tau)*exp(-yst)
46  ELSE
47  x1=sqrt(taup)*exp(yst)
48  x2=sqrt(taup)*exp(-yst)
49  ENDIF
50  xf=x1-x2
51  IF(tau*vint(2).LT.ckin(1)**2) mint(51)=1
52  IF(ckin(2).GE.0..AND.tau*vint(2).GT.ckin(2)**2) mint(51)=1
53  IF(x1.LT.ckin(21).OR.x1.GT.ckin(22)) mint(51)=1
54  IF(x2.LT.ckin(23).OR.x2.GT.ckin(24)) mint(51)=1
55  IF(xf.LT.ckin(25).OR.xf.GT.ckin(26)) mint(51)=1
56  IF(yst.LT.ckin(7).OR.yst.GT.ckin(8)) mint(51)=1
57  IF(iset(isub).EQ.2.OR.iset(isub).EQ.4) THEN
58  pth=0.5*be34*sqrt(tau*vint(2)*(1.-cth**2))
59  y3=yst+0.5*log((1.+rm3-rm4+be34*cth)/(1.+rm3-rm4-be34*cth))
60  y4=yst+0.5*log((1.+rm4-rm3-be34*cth)/(1.+rm4-rm3+be34*cth))
61  ylarge=max(y3,y4)
62  ysmall=min(y3,y4)
63  etalar=10.
64  etasma=-10.
65  sth=sqrt(1.-cth**2)
66  IF(sth.LT.1.e-6) goto 100
67  expet3=((1.+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth+
68  & sqrt(((1.+rm3-rm4)*cosh(yst)+be34*sinh(yst)*cth)**2-4.*rm3))/
69  & (be34*sth)
70  expet4=((1.-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth+
71  & sqrt(((1.-rm3+rm4)*cosh(yst)-be34*sinh(yst)*cth)**2-4.*rm4))/
72  & (be34*sth)
73  eta3=log(min(1.e10,max(1.e-10,expet3)))
74  eta4=log(min(1.e10,max(1.e-10,expet4)))
75  etalar=max(eta3,eta4)
76  etasma=min(eta3,eta4)
77  100 cts3=((1.+rm3-rm4)*sinh(yst)+be34*cosh(yst)*cth)/
78  & sqrt(((1.+rm3-rm4)*cosh(yst)+be34*sinh(yst)*cth)**2-4.*rm3)
79  cts4=((1.-rm3+rm4)*sinh(yst)-be34*cosh(yst)*cth)/
80  & sqrt(((1.-rm3+rm4)*cosh(yst)-be34*sinh(yst)*cth)**2-4.*rm4)
81  ctslar=max(cts3,cts4)
82  ctssma=min(cts3,cts4)
83  IF(pth.LT.pthmin) mint(51)=1
84  IF(ckin(4).GE.0..AND.pth.GT.ckin(4)) mint(51)=1
85  IF(ylarge.LT.ckin(9).OR.ylarge.GT.ckin(10)) mint(51)=1
86  IF(ysmall.LT.ckin(11).OR.ysmall.GT.ckin(12)) mint(51)=1
87  IF(etalar.LT.ckin(13).OR.etalar.GT.ckin(14)) mint(51)=1
88  IF(etasma.LT.ckin(15).OR.etasma.GT.ckin(16)) mint(51)=1
89  IF(ctslar.LT.ckin(17).OR.ctslar.GT.ckin(18)) mint(51)=1
90  IF(ctssma.LT.ckin(19).OR.ctssma.GT.ckin(20)) mint(51)=1
91  IF(cth.LT.ckin(27).OR.cth.GT.ckin(28)) mint(51)=1
92  ENDIF
93  IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) THEN
94  IF(taup*vint(2).LT.ckin(31)**2) mint(51)=1
95  IF(ckin(32).GE.0..AND.taup*vint(2).GT.ckin(32)**2) mint(51)=1
96  ENDIF
97 
98  ELSEIF(ilim.EQ.1) THEN
99 C...Calculate limits on tau
100 C...0) due to definition
101  taumn0=0.
102  taumx0=1.
103 C...1) due to limits on subsystem mass
104  taumn1=ckin(1)**2/vint(2)
105  taumx1=1.
106  IF(ckin(2).GE.0.) taumx1=ckin(2)**2/vint(2)
107 C...2) due to limits on pT-hat (and non-overlapping rapidity intervals)
108  tm3=sqrt(sqm3+pthmin**2)
109  tm4=sqrt(sqm4+pthmin**2)
110  ydcosh=1.
111  IF(ckin(9).GT.ckin(12)) ydcosh=cosh(ckin(9)-ckin(12))
112  taumn2=(tm3**2+2.*tm3*tm4*ydcosh+tm4**2)/vint(2)
113  taumx2=1.
114 C...3) due to limits on pT-hat and cos(theta-hat)
115  cth2mn=min(ckin(27)**2,ckin(28)**2)
116  cth2mx=max(ckin(27)**2,ckin(28)**2)
117  taumn3=0.
118  IF(ckin(27)*ckin(28).GT.0.) taumn3=
119  & (sqrt(sqm3+pthmin**2/(1.-cth2mn))+
120  & sqrt(sqm4+pthmin**2/(1.-cth2mn)))**2/vint(2)
121  taumx3=1.
122  IF(ckin(4).GE.0..AND.cth2mx.LT.1.) taumx3=
123  & (sqrt(sqm3+ckin(4)**2/(1.-cth2mx))+
124  & sqrt(sqm4+ckin(4)**2/(1.-cth2mx)))**2/vint(2)
125 C...4) due to limits on x1 and x2
126  taumn4=ckin(21)*ckin(23)
127  taumx4=ckin(22)*ckin(24)
128 C...5) due to limits on xF
129  taumn5=0.
130  taumx5=max(1.-ckin(25),1.+ckin(26))
131  vint(11)=max(taumn0,taumn1,taumn2,taumn3,taumn4,taumn5)
132  vint(31)=min(taumx0,taumx1,taumx2,taumx3,taumx4,taumx5)
133  IF(mint(43).EQ.1.AND.(iset(isub).EQ.1.OR.iset(isub).EQ.2)) THEN
134  vint(11)=0.99999
135  vint(31)=1.00001
136  ENDIF
137  IF(vint(31).LE.vint(11)) mint(51)=1
138 
139  ELSEIF(ilim.EQ.2) THEN
140 C...Calculate limits on y*
141  IF(iset(isub).EQ.3.OR.iset(isub).EQ.4) tau=vint(26)
142  taurt=sqrt(tau)
143 C...0) due to kinematics
144  ystmn0=log(taurt)
145  ystmx0=-ystmn0
146 C...1) due to explicit limits
147  ystmn1=ckin(7)
148  ystmx1=ckin(8)
149 C...2) due to limits on x1
150  ystmn2=log(max(tau,ckin(21))/taurt)
151  ystmx2=log(max(tau,ckin(22))/taurt)
152 C...3) due to limits on x2
153  ystmn3=-log(max(tau,ckin(24))/taurt)
154  ystmx3=-log(max(tau,ckin(23))/taurt)
155 C...4) due to limits on xF
156  yepmn4=0.5*abs(ckin(25))/taurt
157  ystmn4=sign(log(sqrt(1.+yepmn4**2)+yepmn4),ckin(25))
158  yepmx4=0.5*abs(ckin(26))/taurt
159  ystmx4=sign(log(sqrt(1.+yepmx4**2)+yepmx4),ckin(26))
160 C...5) due to simultaneous limits on y-large and y-small
161  yepsmn=(rm3-rm4)*sinh(ckin(9)-ckin(11))
162  yepsmx=(rm3-rm4)*sinh(ckin(10)-ckin(12))
163  ydifmn=abs(log(sqrt(1.+yepsmn**2)-yepsmn))
164  ydifmx=abs(log(sqrt(1.+yepsmx**2)-yepsmx))
165  ystmn5=0.5*(ckin(9)+ckin(11)-ydifmn)
166  ystmx5=0.5*(ckin(10)+ckin(12)+ydifmx)
167 C...6) due to simultaneous limits on cos(theta-hat) and y-large or
168 C... y-small
169  cthlim=sqrt(1.-4.*pthmin**2/(be34*tau*vint(2)))
170  rzmn=be34*max(ckin(27),-cthlim)
171  rzmx=be34*min(ckin(28),cthlim)
172  yex3mx=(1.+rm3-rm4+rzmx)/max(1e-10,1.+rm3-rm4-rzmx)
173  yex4mx=(1.+rm4-rm3-rzmn)/max(1e-10,1.+rm4-rm3+rzmn)
174  yex3mn=max(1e-10,1.+rm3-rm4+rzmn)/(1.+rm3-rm4-rzmn)
175  yex4mn=max(1e-10,1.+rm4-rm3-rzmx)/(1.+rm4-rm3+rzmx)
176  ystmn6=ckin(9)-0.5*log(max(yex3mx,yex4mx))
177  ystmx6=ckin(12)-0.5*log(min(yex3mn,yex4mn))
178  vint(12)=max(ystmn0,ystmn1,ystmn2,ystmn3,ystmn4,ystmn5,ystmn6)
179  vint(32)=min(ystmx0,ystmx1,ystmx2,ystmx3,ystmx4,ystmx5,ystmx6)
180  IF(mint(43).EQ.1) THEN
181  vint(12)=-0.00001
182  vint(32)=0.00001
183  ELSEIF(mint(43).EQ.2) THEN
184  vint(12)=0.99999*ystmx0
185  vint(32)=1.00001*ystmx0
186  ELSEIF(mint(43).EQ.3) THEN
187  vint(12)=-1.00001*ystmx0
188  vint(32)=-0.99999*ystmx0
189  ENDIF
190  IF(vint(32).LE.vint(12)) mint(51)=1
191 
192  ELSEIF(ilim.EQ.3) THEN
193 C...Calculate limits on cos(theta-hat)
194  yst=vint(22)
195 C...0) due to definition
196  ctnmn0=-1.
197  ctnmx0=0.
198  ctpmn0=0.
199  ctpmx0=1.
200 C...1) due to explicit limits
201  ctnmn1=min(0.,ckin(27))
202  ctnmx1=min(0.,ckin(28))
203  ctpmn1=max(0.,ckin(27))
204  ctpmx1=max(0.,ckin(28))
205 C...2) due to limits on pT-hat
206  ctnmn2=-sqrt(1.-4.*pthmin**2/(be34**2*tau*vint(2)))
207  ctpmx2=-ctnmn2
208  ctnmx2=0.
209  ctpmn2=0.
210  IF(ckin(4).GE.0.) THEN
211  ctnmx2=-sqrt(max(0.,1.-4.*ckin(4)**2/(be34**2*tau*vint(2))))
212  ctpmn2=-ctnmx2
213  ENDIF
214 C...3) due to limits on y-large and y-small
215  ctnmn3=min(0.,max((1.+rm3-rm4)/be34*tanh(ckin(11)-yst),
216  & -(1.-rm3+rm4)/be34*tanh(ckin(10)-yst)))
217  ctnmx3=min(0.,(1.+rm3-rm4)/be34*tanh(ckin(12)-yst),
218  & -(1.-rm3+rm4)/be34*tanh(ckin(9)-yst))
219  ctpmn3=max(0.,(1.+rm3-rm4)/be34*tanh(ckin(9)-yst),
220  & -(1.-rm3+rm4)/be34*tanh(ckin(12)-yst))
221  ctpmx3=max(0.,min((1.+rm3-rm4)/be34*tanh(ckin(10)-yst),
222  & -(1.-rm3+rm4)/be34*tanh(ckin(11)-yst)))
223  vint(13)=max(ctnmn0,ctnmn1,ctnmn2,ctnmn3)
224  vint(33)=min(ctnmx0,ctnmx1,ctnmx2,ctnmx3)
225  vint(14)=max(ctpmn0,ctpmn1,ctpmn2,ctpmn3)
226  vint(34)=min(ctpmx0,ctpmx1,ctpmx2,ctpmx3)
227  IF(vint(33).LE.vint(13).AND.vint(34).LE.vint(14)) mint(51)=1
228 
229  ELSEIF(ilim.EQ.4) THEN
230 C...Calculate limits on tau'
231 C...0) due to kinematics
232  tapmn0=tau
233  tapmx0=1.
234 C...1) due to explicit limits
235  tapmn1=ckin(31)**2/vint(2)
236  tapmx1=1.
237  IF(ckin(32).GE.0.) tapmx1=ckin(32)**2/vint(2)
238  vint(16)=max(tapmn0,tapmn1)
239  vint(36)=min(tapmx0,tapmx1)
240  IF(mint(43).EQ.1) THEN
241  vint(16)=0.99999
242  vint(36)=1.00001
243  ENDIF
244  IF(vint(36).LE.vint(16)) mint(51)=1
245 
246  ENDIF
247  RETURN
248 
249 C...Special case for low-pT and multiple interactions:
250 C...effective kinematical limits for tau, y*, cos(theta-hat).
251  110 IF(ilim.EQ.0) THEN
252  ELSEIF(ilim.EQ.1) THEN
253  IF(mstp(82).LE.1) vint(11)=4.*parp(81)**2/vint(2)
254  IF(mstp(82).GE.2) vint(11)=parp(82)**2/vint(2)
255  vint(31)=1.
256  ELSEIF(ilim.EQ.2) THEN
257  vint(12)=0.5*log(vint(21))
258  vint(32)=-vint(12)
259  ELSEIF(ilim.EQ.3) THEN
260  IF(mstp(82).LE.1) st2eff=4.*parp(81)**2/(vint(21)*vint(2))
261  IF(mstp(82).GE.2) st2eff=0.01*parp(82)**2/(vint(21)*vint(2))
262  vint(13)=-sqrt(max(0.,1.-st2eff))
263  vint(33)=0.
264  vint(14)=0.
265  vint(34)=-vint(13)
266  ENDIF
267 
268  RETURN
269  END