]>
Commit | Line | Data |
---|---|---|
1 | c---------------------------------------------------------------------- | |
2 | c The two elementary fuctions of our approach, the profile fuction G | |
3 | c and the Phi exponent G~, are here referred to as Gpro and Gtilde. | |
4 | c Both functions can be written as | |
5 | c | |
6 | c G = \sum_type alp * xp**bet * xm**betp | |
7 | c | |
8 | c The parameters alp, bet, betp are calculated in GfunParK (k-mode, | |
9 | c b is takento the one of pair k) or GfunPar (b-mode: arbitrary b) as | |
10 | c | |
11 | c Gpro: bet = betD' + epsGp + gamD*b**2 + epsG -alppar | |
12 | c bet' = betD'' + epsGt + gamD*b**2 + epsG -alppar | |
13 | c | |
14 | c alp = alpD*f * s**(betD+gamD*b**2+epsG) * exp(-b**2/delD) | |
15 | c | |
16 | c Gtilde: bet~ = bet + 1 | |
17 | c bet~' = bet' + 1 | |
18 | c | |
19 | c alp~ = alp * gam(bet~) * gam(bet~') | |
20 | c * gam(1+alppro) * gam(1+alptar) | |
21 | c * gam(1+alppro+bet~) * gam(1+alptar+bet~') | |
22 | c * (1+epsGt') * (1+epsGt') | |
23 | c | |
24 | c The parameters depend implicitely on s. | |
25 | c | |
26 | c In the program we use om1 = Gpro | |
27 | c (they differ by a constant which is actually one) | |
28 | c and om5 = om1 * 0.5 | |
29 | c All functions related to om1 are called om1... . | |
30 | c | |
31 | c The inclusive Pomeron distributions are | |
32 | c | |
33 | c PomInc(xp,xm) = Gpro(xp,xm) * (1-xp)**alppro * (1-xm)**alptar | |
34 | c | |
35 | c---------------------------------------------------------------------- | |
36 | ||
37 | ||
38 | c---------------------------------------------------------------------- | |
39 | subroutine GfunParK(irea) !---MC--- | |
40 | c---------------------------------------------------------------------- | |
41 | c calculates parameters alp,bet,betp of the G functions (k-mode) | |
42 | c and the screening exponents epsilongp(k,i), epsilongt(k,i), epsilongs(k,i) | |
43 | c---------------------------------------------------------------------- | |
44 | c Gpro parameters written to /comtilde/atilde(,)btildep(,),btildepp(,) | |
45 | c Gtilde parameters written to /cgtilde/atildg(,)btildgp(,),btildgpp(,) | |
46 | c two subscripts: first=type, second=collision k | |
47 | c Certain pieces of this routine are only done if irea is <= or = zero. | |
48 | c---------------------------------------------------------------------- | |
49 | include 'epos.inc' | |
50 | include 'epos.incsem' | |
51 | include 'epos.incems' | |
52 | include 'epos.incpar' | |
53 | double precision atildg,btildgp,btildgpp | |
54 | common/cgtilde/atildg(idxD0:idxD1,kollmx) | |
55 | *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx) | |
56 | double precision utgam2,coefgdp,coefgdt | |
57 | parameter(nbkbin=40) | |
58 | common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin | |
59 | common /cgtilnu/ cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro | |
60 | &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam | |
61 | double precision cfbetpnp,cfbetppnp,cfbetpnm,cfbetppnm,cfalpro | |
62 | &,cfaltar,cfbpap,cfbpam,cfbppap,cfbppam,gamv,eps | |
63 | parameter (eps=1.d-25) | |
64 | parameter(nxeps=20,nyeps=32) | |
65 | common/cxeps1/w(0:nxeps,nyeps),y1(nyeps),y2(nyeps) | |
66 | common/cxeps2/db,b1,b2 | |
67 | common/geom/rmproj,rmtarg,bmax,bkmx | |
68 | common/nucl3/phi,bimp | |
69 | ||
70 | b1=0.03 | |
71 | b2=bkmx*1.2 | |
72 | db=(b2-b1)/nyeps | |
73 | ||
74 | call utpri('GfunParK',ish,ishini,5) | |
75 | ||
76 | cfbetpnp=0.d0 | |
77 | cfbetppnp=0.d0 | |
78 | cfbetpnm=0.d0 | |
79 | cfbetppnm=0.d0 | |
80 | cfalpro=dble(ucfpro) | |
81 | cfaltar=dble(ucftar) | |
82 | cfbpap=1.d0 | |
83 | cfbppap=1.d0 | |
84 | cfbpam=1.d0 | |
85 | cfbppam=1.d0 | |
86 | alpfom=0. | |
87 | sy=engy*engy | |
88 | ||
89 | do k=1,koll | |
90 | do i=ntymi,ntymx | |
91 | atilde(i,k)=0.d0 | |
92 | btildep(i,k)=0.d0 | |
93 | btildepp(i,k)=0.d0 | |
94 | enddo | |
95 | do i=idxD0,idxD1 | |
96 | atildg(i,k)=0.d0 | |
97 | btildgp(i,k)=0.d0 | |
98 | btildgpp(i,k)=0.d0 | |
99 | enddo | |
100 | enddo | |
101 | ||
102 | ||
103 | ! calculate collision number according to Glauber ----------------------- | |
104 | ||
105 | bglaub2=max(1.,sigine/10./pi) !10= fm^2 -> mb | |
106 | bglaub=sqrt(bglaub2) | |
107 | nglevt=0 | |
108 | do ko=1,koll | |
109 | r=bk(ko) | |
110 | if(r.le.bglaub)nglevt=nglevt+1 | |
111 | enddo | |
112 | ||
113 | ! Z_parton_projectile (zparpro) and Z_parton_target (zpartar)----------- | |
114 | ||
115 | if(iscreen.eq.1.or.isplit.eq.1)then | |
116 | ||
117 | bcut=0. !zbcut | |
118 | ctp060829 bzero=bglaub | |
119 | b2x=epscrp*epscrp*bglaub2 | |
120 | alpfom=alpfomi*epscrw*fscra(engy/egyscr) | |
121 | ztav=0. | |
122 | zpav=0. | |
123 | ||
124 | do k=1,koll | |
125 | ip=iproj(k) | |
126 | it=itarg(k) | |
127 | !....... targ partons seen by proj | |
128 | if(lproj(ip).ge.1)then | |
129 | absb=max(1.e-9,bk(k)) | |
130 | b2=absb*absb | |
131 | zkp=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr) | |
132 | c epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar, | |
133 | c & epscrs*min(epscrx,zkp*abs(betD(0,iclpro,icltar)))) | |
134 | c epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar, | |
135 | c & epscrh*min(epscrx,zkp*abs(betD(1,iclpro,icltar)))) | |
136 | if(lproj3(ip).ge.1)then | |
137 | zpartar(k)=zkp | |
138 | do lt=1,lproj3(ip) | |
139 | kp=kproj3(ip,lt) | |
140 | if(kp.ne.k)then | |
141 | absb=max(1.e-9,abs(bk(kp))-bcut) | |
142 | b2=absb*absb | |
143 | zkp=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr) | |
144 | endif | |
145 | zpartar(k)=zpartar(k)+zkp | |
146 | c alpfom=max(alpfom,dble(zpartar(k))) | |
147 | enddo | |
148 | endif | |
149 | else | |
150 | zpartar(k)=0. | |
151 | endif | |
152 | ztav=ztav+zpartar(k) | |
153 | !...........proj partons seen by targ | |
154 | if(ltarg(it).ge.1)then | |
155 | absb=max(1.e-9,bk(k)) | |
156 | b2=absb*absb | |
157 | zkt=epscrw*exp(-b2/2./b2x)*fscra(engy/egyscr) | |
158 | c epsilongt(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar, | |
159 | c & epscrs*min(epscrx,zkt*abs(betD(0,iclpro,icltar)))) | |
160 | c epsilongt(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar, | |
161 | c & epscrh*min(epscrx,zkt*abs(betD(1,iclpro,icltar)))) | |
162 | if(ltarg3(it).ge.1)then | |
163 | zparpro(k)=zkt | |
164 | do lp=1,ltarg3(it) | |
165 | kt=ktarg3(it,lp) | |
166 | if(kt.ne.k)then | |
167 | absb=max(1.e-9,abs(bk(kt))-bcut) | |
168 | b2=absb*absb | |
169 | zkt=epscrw*exp(-b2/2./b2x)*fscro(engy/egyscr) | |
170 | endif | |
171 | zparpro(k)=zparpro(k)+zkt | |
172 | c alpfom=max(alpfom,dble(zparpro(k))) | |
173 | enddo | |
174 | endif | |
175 | else | |
176 | zparpro(k)=0. | |
177 | endif | |
178 | zpav=zpav+zparpro(k) | |
179 | enddo | |
180 | ||
181 | else ! no screening | |
182 | ||
183 | do k=1,koll | |
184 | zparpro(k)=0. | |
185 | zpartar(k)=0. | |
186 | enddo | |
187 | ||
188 | endif | |
189 | ||
190 | c calculation of epsilongp epsilongt | |
191 | ||
192 | if(iscreen.eq.1)then | |
193 | ||
194 | !ip=0 | |
195 | do k=1,koll | |
196 | !ipp=ip | |
197 | epsG=0. | |
198 | epsilongs(k,0)=epsG | |
199 | epsilongs(k,1)=epsG | |
200 | ip=iproj(k) !...........projectile | |
201 | if(lproj(ip).gt.0)then | |
202 | x=zparpro(k) | |
203 | epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar, | |
204 | & min(epscrx,epscrs*x)) | |
205 | c & epscrs*x) | |
206 | c & epscrs*x*abs(betDp(0,iclpro,icltar))) | |
207 | c & epscrs*min(epscrx,x*abs(betDp(0,iclpro,icltar)))) | |
208 | epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar, | |
209 | & epscrh*min(epscrx,x*abs(betDp(1,iclpro,icltar)))) | |
210 | cc epsilongp(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar))) | |
211 | cc epsilongp(k,1)=epscrh*min(epscrx,x*abs(betD(1,iclpro,icltar))) | |
212 | gammaV(k)=1.d0 | |
213 | else | |
214 | epsilongp(k,0)=0. | |
215 | epsilongp(k,1)=0. | |
216 | gammaV(k)=1.d0 | |
217 | endif | |
218 | ||
219 | it=itarg(k) !...........target | |
220 | if(ltarg(it).gt.0)then | |
221 | x=zpartar(k) | |
222 | epsilongt(k,0)=max(-betDpp(0,iclpro,icltar)-0.95+alppar, | |
223 | & min(epscrx,epscrs*x)) | |
224 | c & epscrs*x) | |
225 | c & epscrs*x*abs(betDpp(0,iclpro,icltar))) | |
226 | c & epscrs*min(epscrx,x*abs(betDpp(0,iclpro,icltar)))) | |
227 | epsilongt(k,1)=max(-betDpp(1,iclpro,icltar)-0.95+alppar, | |
228 | & epscrh*min(epscrx,x*abs(betDpp(1,iclpro,icltar)))) | |
229 | cc epsilongt(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar))) | |
230 | cc epsilongt(k,1)=epscrh*min(epscrx,x*abs(betD(1,iclpro,icltar))) | |
231 | cc gammaV(k)=gammaV(k) | |
232 | else | |
233 | epsilongt(k,0)=0. | |
234 | epsilongt(k,1)=0. | |
235 | gammaV(k)=gammaV(k) | |
236 | endif | |
237 | ||
238 | enddo | |
239 | ||
240 | else ! no screening | |
241 | ||
242 | do k=1,koll | |
243 | epsilongs(k,0)=0. | |
244 | epsilongs(k,1)=0. | |
245 | epsilongp(k,0)=0. | |
246 | epsilongp(k,1)=0. | |
247 | epsilongt(k,0)=0. | |
248 | epsilongt(k,1)=0. | |
249 | gammaV(k)=1.d0 | |
250 | enddo | |
251 | ||
252 | endif | |
253 | ||
254 | ||
255 | !..............alpha beta betap for Gtilde (->PhiExpo)....................... | |
256 | ||
257 | imax=idxD1 | |
258 | if(iomega.eq.2)imax=1 | |
259 | ||
260 | do k=1,koll | |
261 | ||
262 | b=bk(k) | |
263 | b2=bk(k)*bk(k) | |
264 | ip=iproj(k) | |
265 | it=itarg(k) | |
266 | ||
267 | if(b.lt.(nbkbin-1)*bkbin)then | |
268 | ibk=int(bk(k)/bkbin)+1 | |
269 | if(isetcs.gt.1.and.iclegy.lt.iclegy2)then | |
270 | egy0=egylow*egyfac**real(iclegy-1) | |
271 | xkappa1=xkappafit(iclegy,iclpro,icltar,ibk) | |
272 | * +(bk(k)-bkbin*real(ibk-1))/bkbin | |
273 | * *(xkappafit(iclegy,iclpro,icltar,ibk+1) | |
274 | * -xkappafit(iclegy,iclpro,icltar,ibk)) | |
275 | xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk) | |
276 | * +(bk(k)-bkbin*real(ibk-1))/bkbin | |
277 | * *(xkappafit(iclegy+1,iclpro,icltar,ibk+1) | |
278 | * -xkappafit(iclegy+1,iclpro,icltar,ibk)) | |
279 | xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac) | |
280 | * *(log(engy)-log(egy0)) | |
281 | else | |
282 | xkappa=xkappafit(iclegy,iclpro,icltar,ibk) | |
283 | * +(bk(k)-bkbin*real(ibk-1))/bkbin | |
284 | * *(xkappafit(iclegy,iclpro,icltar,ibk+1) | |
285 | * -xkappafit(iclegy,iclpro,icltar,ibk)) | |
286 | endif | |
287 | else | |
288 | xkappa=1. | |
289 | endif | |
290 | gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub) | |
291 | gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub) | |
292 | ||
293 | do i=idxDmin,imax | |
294 | gamV=gammaV(k) | |
295 | if(i.lt.2)then | |
296 | epsG=epsilongs(k,i) | |
297 | else | |
298 | epsG=0. | |
299 | endif | |
300 | gamb=gamD(i,iclpro,icltar)*b2+epsG | |
301 | atildg(i,k)=dble(alpD(i,iclpro,icltar)) | |
302 | * *cfalpro*cfaltar | |
303 | * *gamv | |
304 | c if(i.eq.0) atildg(i,k)=atildg(i,k) | |
305 | atildg(i,k)=atildg(i,k) | |
306 | * *dble(xkappa*xkappa) | |
307 | if(i.lt.2)then | |
308 | atildg(i,k)=atildg(i,k)*dble( | |
309 | * chad(iclpro)*chad(icltar) | |
310 | * *exp(-b2/delD(i,iclpro,icltar)) | |
311 | * *sy**(betD(i,iclpro,icltar)+gamb)) | |
312 | * *gfactorp *gfactort | |
313 | epsGp=epsilongp(k,i) | |
314 | epsGt=epsilongt(k,i) | |
315 | btildgp(i,k)=dble(betDp(i,iclpro,icltar) | |
316 | * +epsGp | |
317 | * +gamb-alppar)+1.d0 | |
318 | btildgpp(i,k)=dble(betDpp(i,iclpro,icltar) | |
319 | * +epsGt | |
320 | * +gamb-alppar)+1.d0 | |
321 | else | |
322 | absb=abs(bk(k))-bmxdif(iclpro,icltar) | |
323 | b2a=absb*absb | |
324 | atildg(i,k)=atildg(i,k)*dble( | |
325 | * sy**betD(i,iclpro,icltar) | |
326 | * *exp(-b2a/delD(i,iclpro,icltar))) | |
327 | btildgp(i,k)=dble(betDp(i,iclpro,icltar)-alppar)+1.d0 | |
328 | btildgpp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar)+1.d0 | |
329 | endif | |
330 | coefgdp=utgam2(1.d0+dble(alplea(iclpro))+btildgp(i,k)) | |
331 | coefgdt=utgam2(1.d0+dble(alplea(icltar))+btildgpp(i,k)) | |
332 | atildg(i,k)=atildg(i,k) | |
333 | * *utgam2(btildgp(i,k))*utgam2(btildgpp(i,k)) | |
334 | * /coefgdp/coefgdt | |
335 | !...........prepare plot in xepsilon | |
336 | if(irea.eq.0.and.i.lt.2)then | |
337 | kk=max(1,int((bk(k)-b1)/db)+1) | |
338 | if(i.eq.0)w(0,kk)=w(0,kk)+1 | |
339 | if(i.eq.0)w(1,kk)=w(1,kk)+epsGp | |
340 | if(i.eq.0)w(2,kk)=w(2,kk)+epsGt | |
341 | if(i.eq.0)w(3+i,kk)=w(3+i,kk)+0 | |
342 | if(i.eq.1)w(3+i,kk)=w(3+i,kk)+0 | |
343 | !...5-8 soft ... 9-12 semi | |
344 | w(5+4*i,kk)=w(5+4*i,kk) | |
345 | * +betDp(i,iclpro,icltar) !prj eff | |
346 | * +epsGp+gamb | |
347 | w(6+4*i,kk)=w(6+4*i,kk) | |
348 | * +betDpp(i,iclpro,icltar) !tgt eff | |
349 | * +epsGt+gamb | |
350 | w(7+4*i,kk)=w(7+4*i,kk) | |
351 | * +betDp(i,iclpro,icltar) !prj unscr | |
352 | * +gamb | |
353 | w(8+4*i,kk)=w(8+4*i,kk) | |
354 | * +betDpp(i,iclpro,icltar) !tgt unscr | |
355 | * +gamb | |
356 | if(i.eq.0)w(13,kk)=w(13,kk)+zparpro(k) | |
357 | if(i.eq.0)w(14,kk)=w(14,kk)+zpartar(k) | |
358 | endif | |
359 | !................ | |
360 | enddo | |
361 | enddo | |
362 | ||
363 | !........................................................................... | |
364 | ||
365 | zppevt=zpav/koll | |
366 | zptevt=ztav/koll | |
367 | ktot=int(bimp)+1 | |
368 | if(ktot.le.nyeps)then | |
369 | w(15,ktot)=w(15,ktot)+zpav/koll | |
370 | w(16,ktot)=w(16,ktot)+ztav/koll | |
371 | w(17,ktot)=w(17,ktot)+1 | |
372 | endif | |
373 | n=1+float(nglevt)/(0.1*maproj*matarg)*(nyeps-1) | |
374 | if(nglevt.ge.1.and.n.ge.1.and.n.le.nyeps)then | |
375 | w(18,n)=w(18,n)+zpav/koll | |
376 | w(19,n)=w(19,n)+ztav/koll | |
377 | w(20,n)=w(20,n)+1 | |
378 | endif | |
379 | ||
380 | ||
381 | !........alpha beta betap for Gpro........................................... | |
382 | ||
383 | if(irea.le.0)then | |
384 | do k=1,koll | |
385 | ip=iproj(k) | |
386 | it=itarg(k) | |
387 | ||
388 | b2=bk(k)*bk(k) | |
389 | imax=ntymx | |
390 | if(iomega.eq.2)imax=1 | |
391 | do i=ntymin,imax | |
392 | ||
393 | if(i.lt.2)then | |
394 | epsG=epsilongs(k,i) | |
395 | else | |
396 | epsG=0. | |
397 | endif | |
398 | gamb=gamD(i,iclpro,icltar)*b2+epsG | |
399 | ||
400 | atilde(i,k)=dble(alpD(i,iclpro,icltar)) | |
401 | if(i.lt.2)then | |
402 | atilde(i,k)=atilde(i,k)*dble( | |
403 | * exp(-b2/delD(i,iclpro,icltar)) | |
404 | * *sy**(betD(i,iclpro,icltar) | |
405 | * +gamb) | |
406 | * *chad(iclpro)*chad(icltar)) | |
407 | epsGp=epsilongp(k,i) | |
408 | epsGt=epsilongt(k,i) | |
409 | btildep(i,k)=dble(betDp(i,iclpro,icltar) | |
410 | * +epsGp | |
411 | * +gamb-alppar) | |
412 | btildepp(i,k)=dble(betDpp(i,iclpro,icltar) | |
413 | * +epsGt | |
414 | * +gamb-alppar) | |
415 | else | |
416 | absb=abs(bk(k))-bmxdif(iclpro,icltar) | |
417 | b2a=absb*absb | |
418 | atilde(i,k)=atilde(i,k)*dble( | |
419 | * sy**betD(i,iclpro,icltar) | |
420 | * *exp(-b2a/delD(i,iclpro,icltar))) | |
421 | ||
422 | btildep(i,k)=dble(betDp(i,iclpro,icltar)-alppar) | |
423 | btildepp(i,k)=dble(betDpp(i,iclpro,icltar)-alppar) | |
424 | endif | |
425 | ||
426 | if(btildep(i,k)+1.d0.lt.-eps.or.btildepp(i,k)+1.d0.lt.-eps) | |
427 | & call utstop('Error in epos-omg in GfunPark&') | |
428 | enddo | |
429 | ||
430 | enddo | |
431 | endif | |
432 | ||
433 | !........................................................................... | |
434 | ||
435 | if(ish.ge.8)then | |
436 | do k=1,koll | |
437 | ip=iproj(k) | |
438 | it=itarg(k) | |
439 | write(ifch,*)' k,b,ip,it,',k,bk(k),ip,it | |
440 | write(ifch,*)' gammaV,epsilonGP1/2,epsilonGP1/2' | |
441 | * ,gammaV(k),epsilongp(k,0),epsilongp(k,1) | |
442 | * ,epsilongt(k,0),epsilongt(k,1) | |
443 | write(ifch,*)'*******************************************' | |
444 | write(ifch,*)" atilde,btildep,btildepp " | |
445 | do i=ntymin,ntymx | |
446 | write(ifch,*)i,atilde(i,k),btildep(i,k),btildepp(i,k) | |
447 | enddo | |
448 | enddo | |
449 | endif | |
450 | ||
451 | call utprix('GfunParK',ish,ishini,5) | |
452 | ||
453 | return | |
454 | end | |
455 | ||
456 | c---------------------------------------------------------------------- | |
457 | subroutine GfunPar(m,i,b,spp,alp,bet,betp,epsp,epst,epss,gamvv) | |
458 | c---------------------------------------------------------------------- | |
459 | c calculates parameters alp,bet,betp of the G functions for pp (b-mode) | |
460 | c and the screening exponents epsp,epst,epss,gamvv | |
461 | c---------------------------------------------------------------------- | |
462 | c m=1: profile function Gpro, i=0: soft i=2: diff | |
463 | c m=2: Gtilde, i=1: semi (no screening for diff) | |
464 | c b: impact param, spp: pp energy squared | |
465 | c---------------------------------------------------------------------- | |
466 | ||
467 | include 'epos.inc' | |
468 | include 'epos.incsem' | |
469 | include 'epos.incpar' | |
470 | include 'epos.incems' | |
471 | parameter(nbkbin=40) | |
472 | common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin | |
473 | common /kwrite/ xkapZ | |
474 | parameter(eps=1.e-20) | |
475 | double precision utgam2 | |
476 | ||
477 | call utprj('GfunPar ',ish,ishini,6) | |
478 | ||
479 | !???????????????????????????????????????????????????estimation | |
480 | x=sqrt(spp) | |
481 | if(sigine.eq.0.)then | |
482 | if(iclpro+icltar.eq.3)then !pi+p | |
483 | sigtotf=15.*x**0.15+55.*x**(-1.5) !fit to data | |
484 | sigelaf=20.*(x-1)**(-3)+6.*x**(-0.4)+0.2*alog(x)**2. !fit to data | |
485 | elseif(iclpro+icltar.eq.5)then !K+p | |
486 | sigtotf=12.5*x**0.15+35.*x**(-1.5) !fit to data | |
487 | sigelaf=15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data | |
488 | else | |
489 | sigtotf=21*x**0.17+44.*x**(-0.8) !fit to data | |
490 | sigelaf=30.*(x-1)**(-3)+17*x**(-0.47)+0.3*alog(x)**2 !fit to data | |
491 | endif | |
492 | siginex=sigtotf-sigelaf | |
493 | else | |
494 | siginex=sigine | |
495 | endif | |
496 | bglaub2=max(1.,siginex/10./pi) | |
497 | bglaub=sqrt(bglaub2) | |
498 | if(ish.ge.6)write(ifch,*)'Gf in:',m,i,x,b,sigine,siginex,bglaub | |
499 | & ,iclpro,icltar | |
500 | !!!!!print*,'++++++++++++++++++',x,siginex,bglaub | |
501 | b2=b*b | |
502 | ee=sqrt(spp) | |
503 | cfalpro=ucfpro | |
504 | cfaltar=ucftar | |
505 | gamb=gamD(i,iclpro,icltar)*b2 | |
506 | ||
507 | if(iscreen.ne.0)then | |
508 | bcut=0.!zbcut | |
509 | ctp060829 bzero=bglaub | |
510 | absb=max(1.e-9,abs(b)-bcut) | |
511 | b2a=absb*absb | |
512 | b2x=epscrp*epscrp*bglaub2 | |
513 | zzp=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr) | |
514 | zzt=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr) | |
515 | x=zzp | |
516 | if(i.eq.0)then | |
517 | epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar | |
518 | & ,min(epscrx,epscrs*x)) | |
519 | c & ,epscrs*x) | |
520 | c & ,epscrs*x*abs(betDp(i,iclpro,icltar))) | |
521 | c & ,epscrs*min(epscrx,x*abs(betDp(i,iclpro,icltar)))) | |
522 | c epsG=epsG+epsGp | |
523 | elseif(i.eq.1)then | |
524 | epsGp=max(-betDp(i,iclpro,icltar)-0.95+alppar | |
525 | & ,epscrh*min(epscrx,x*abs(betDp(i,iclpro,icltar)))) | |
526 | else | |
527 | epsGp=0. | |
528 | endif | |
529 | gamV=1. | |
530 | x=zzt | |
531 | if(i.eq.0)then | |
532 | epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar | |
533 | & ,min(epscrx,epscrs*x)) | |
534 | c & ,epscrs*x) | |
535 | c & ,epscrs*x*abs(betDpp(i,iclpro,icltar))) | |
536 | c & ,epscrs*min(epscrx,x*abs(betDpp(i,iclpro,icltar)))) | |
537 | c epsG=epsG+epsGt | |
538 | elseif(i.eq.1)then | |
539 | epsGt=max(-betDpp(i,iclpro,icltar)-0.95+alppar | |
540 | & ,epscrh*min(epscrx,x*abs(betDpp(i,iclpro,icltar)))) | |
541 | else | |
542 | epsGt=0. | |
543 | endif | |
544 | c gamV=gamV | |
545 | epsG=0. | |
546 | else | |
547 | zzp=0. | |
548 | zzt=0. | |
549 | epsGp=0. | |
550 | epsGt=0. | |
551 | epsG=0. | |
552 | gamV=1. | |
553 | endif | |
554 | ||
555 | gfactorp=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub) | |
556 | gfactort=1.!+(gfactor-1)*exp(-5*b/gwidth/bglaub) | |
557 | ||
558 | rho=betD(i,iclpro,icltar) | |
559 | * +gamb+epsG | |
560 | ||
561 | ||
562 | if(m.eq.1)then | |
563 | ||
564 | alp=alpD(i,iclpro,icltar) | |
565 | if(i.lt.2)then | |
566 | alp=alp | |
567 | * *exp(rho*log(spp)-b2/delD(i,iclpro,icltar)) | |
568 | bet=betDp(i,iclpro,icltar) | |
569 | * +epsGp | |
570 | * +gamb-alppar | |
571 | betp=betDpp(i,iclpro,icltar) | |
572 | * +epsGt | |
573 | * +gamb-alppar | |
574 | else | |
575 | absb=abs(b)-bmxdif(iclpro,icltar) | |
576 | b2a=absb*absb | |
577 | alp=alp | |
578 | * *exp(betD(i,iclpro,icltar)*log(spp) | |
579 | * -b2a/delD(i,iclpro,icltar)) | |
580 | bet=betDp(i,iclpro,icltar)-alppar | |
581 | betp=betDpp(i,iclpro,icltar)-alppar | |
582 | endif | |
583 | ||
584 | if((bet+1.0).lt.-eps.or.(betp+1.0).lt.-eps)then | |
585 | write(*,*)'m,i,b,spp,alp,bet,betp',m,i,b,spp,alp,bet,betp | |
586 | call utstop('Error : beta < -1 in Gfunpar in epos-omg&') | |
587 | endif | |
588 | ||
589 | elseif(m.eq.2)then | |
590 | xkappa=1. | |
591 | c if(i.eq.0.and.b.lt.(nbkbin-1)*bkbin)then | |
592 | if(b.lt.(nbkbin-1)*bkbin)then | |
593 | ibk=int(b/bkbin)+1 | |
594 | if(isetcs.gt.1.and.iclegy.lt.iclegy2)then | |
595 | egy0=egylow*egyfac**real(iclegy-1) | |
596 | xkappa1=xkappafit(iclegy,iclpro,icltar,ibk) | |
597 | * +(b-bkbin*real(ibk-1))/bkbin | |
598 | * *(xkappafit(iclegy,iclpro,icltar,ibk+1) | |
599 | * -xkappafit(iclegy,iclpro,icltar,ibk)) | |
600 | xkappa2=xkappafit(iclegy+1,iclpro,icltar,ibk) | |
601 | * +(b-bkbin*real(ibk-1))/bkbin | |
602 | * *(xkappafit(iclegy+1,iclpro,icltar,ibk+1) | |
603 | * -xkappafit(iclegy+1,iclpro,icltar,ibk)) | |
604 | xkappa=xkappa1+(xkappa2-xkappa1)/log(egyfac) | |
605 | * *(log(ee)-log(egy0)) | |
606 | else | |
607 | xkappa=xkappafit(iclegy,iclpro,icltar,ibk) | |
608 | * +(b-bkbin*real(ibk-1))/bkbin | |
609 | * *(xkappafit(iclegy,iclpro,icltar,ibk+1) | |
610 | * -xkappafit(iclegy,iclpro,icltar,ibk)) | |
611 | endif | |
612 | xkapZ=xkappa | |
613 | endif | |
614 | ||
615 | alp=alpD(i,iclpro,icltar) | |
616 | * *cfalpro*cfaltar | |
617 | * *gamV | |
618 | c if(i.eq.0)alp=alp | |
619 | alp=alp | |
620 | * *xkappa*xkappa | |
621 | ||
622 | if(i.lt.2)then | |
623 | alp=alp | |
624 | * *exp(rho*log(spp)-b2/delD(i,iclpro,icltar)) | |
625 | * *chad(iclpro)*chad(icltar) | |
626 | * *gfactorp *gfactort | |
627 | bet=betDp(i,iclpro,icltar) | |
628 | * +epsGp | |
629 | * +gamb-alppar+1. | |
630 | betp=betDpp(i,iclpro,icltar) | |
631 | * +epsGt | |
632 | * +gamb-alppar+1. | |
633 | else | |
634 | absb=abs(b)-bmxdif(iclpro,icltar) | |
635 | b2a=absb*absb | |
636 | alp=alp | |
637 | * *exp(betD(i,iclpro,icltar)*log(spp) | |
638 | * -b2a/delD(i,iclpro,icltar)) | |
639 | bet=betDp(i,iclpro,icltar)-alppar+1. | |
640 | betp=betDpp(i,iclpro,icltar)-alppar+1. | |
641 | endif | |
642 | coefgdp=utgam2(1.d0+dble(alplea(iclpro))+bet) | |
643 | coefgdt=utgam2(1.d0+dble(alplea(icltar))+betp) | |
644 | alp=alp*utgam2(dble(bet))*utgam2(dble(betp))/coefgdp/coefgdt | |
645 | else | |
646 | ||
647 | stop'GproPar: wrong m value. ' | |
648 | ||
649 | endif | |
650 | ||
651 | epsp=epsGp | |
652 | epst=epsGt | |
653 | epss=epsG | |
654 | gamvv=gamV | |
655 | ||
656 | alpUni(i,m)=alp | |
657 | betUni(i,m)=bet | |
658 | betpUni(i,m)=betp | |
659 | ||
660 | ||
661 | if(ish.ge.6)write(ifch,*)' GfunPar :',alp,bet,betp,epsp,epst | |
662 | & ,epss,gamvv | |
663 | ||
664 | call utprjx('GfunPar ',ish,ishini,6) | |
665 | end | |
666 | ||
667 | c---------------------------------------------------------------------- | |
668 | subroutine GfomPar(b,spp) | |
669 | c---------------------------------------------------------------------- | |
670 | c calculates parameters of the fom functions for pp (b-mode) | |
671 | c---------------------------------------------------------------------- | |
672 | c b: impact param, spp: pp energy squared | |
673 | c---------------------------------------------------------------------- | |
674 | ||
675 | include 'epos.inc' | |
676 | include 'epos.incsem' | |
677 | include 'epos.incpar' | |
678 | include 'epos.incems' | |
679 | parameter(eps=1.e-20) | |
680 | ||
681 | call utprj('GfomPar ',ish,ishini,6) | |
682 | ||
683 | !???????????????????????????????????????????????????estimation | |
684 | x=sqrt(spp) | |
685 | if(sigine.eq.0.)then | |
686 | if(iclpro+icltar.eq.3)then !pi+p | |
687 | sigtotf=15.*x**0.15+55.*x**(-1.5) !fit to data | |
688 | sigelaf=20.*(x-1)**(-3)+6.*x**(-0.4)+0.2*alog(x)**2. !fit to data | |
689 | elseif(iclpro+icltar.eq.5)then !K+p | |
690 | sigtotf=12.5*x**0.15+35.*x**(-1.5) !fit to data | |
691 | sigelaf=15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data | |
692 | else | |
693 | sigtotf=21*x**0.17+44.*x**(-0.8) !fit to data | |
694 | sigelaf=30.*(x-1)**(-3)+17*x**(-0.47)+0.3*alog(x)**2 !fit to data | |
695 | endif | |
696 | siginex=sigtotf-sigelaf | |
697 | else | |
698 | siginex=sigine | |
699 | endif | |
700 | bglaub2=max(1.,siginex/10./pi) | |
701 | bglaub=sqrt(bglaub2) | |
702 | b2=b*b | |
703 | ee=sqrt(spp) | |
704 | ||
705 | if(iscreen.ne.0)then | |
706 | bcut=0.!zbcut | |
707 | ctp060829 bzero=bglaub | |
708 | absb=max(1.e-9,abs(b)-bcut) | |
709 | b2a=absb*absb | |
710 | b2x=epscrp*epscrp*bglaub2 | |
711 | zzp=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr) | |
712 | zzt=(epscrw*exp(-b2a/2./b2x))*fscra(ee/egyscr) | |
713 | else | |
714 | zzp=0. | |
715 | zzt=0. | |
716 | endif | |
717 | ||
718 | ||
719 | z0=alpfomi!*epscrw*fscra(ee/egyscr) | |
720 | if(z0.gt.0.)then | |
721 | z1=zzp | |
722 | zpUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar))) | |
723 | c zpUni=dble(4.*z0*(z1/z0)**1.5) | |
724 | z1=zzt | |
725 | ztUni=dble(z1**gamfom/z0)*exp(-dble(b*b/delD(1,iclpro,icltar))) | |
726 | c ztUni=dble(4.*z0*(z1/z0)**1.5) | |
727 | else | |
728 | zpUni=0d0 | |
729 | ztUni=0d0 | |
730 | endif | |
731 | ||
732 | if(ish.ge.6)write(ifch,*)' GfomPar :',zpUni,ztUni | |
733 | ||
734 | call utprjx('GfomPar ',ish,ishini,6) | |
735 | end | |
736 | ||
737 | c---------------------------------------------------------------------- | |
738 | function fscra(x) | |
739 | c---------------------------------------------------------------------- | |
740 | fscra=0 | |
741 | x2=x*x | |
742 | if(x2.gt.1.)fscra=log(x2) | |
743 | end | |
744 | ||
745 | c---------------------------------------------------------------------- | |
746 | function fscro(x) | |
747 | c---------------------------------------------------------------------- | |
748 | include 'epos.incpar' | |
749 | fscro=epscrv | |
750 | x2=x*x | |
751 | if(x2.gt.1.)fscro=sqrt(log(x2)**2+epscrv**2) | |
752 | end | |
753 | ||
754 | c---------------------------------------------------------------------- | |
755 | subroutine recalcZPtn !???????????? not updated !!! | |
756 | c---------------------------------------------------------------------- | |
757 | ||
758 | include 'epos.inc' | |
759 | include 'epos.incems' | |
760 | stop'recalcZPtn not valid any more!!!!!!!' | |
761 | c if(koll.eq.1.and.maproj.eq.1.and.matarg.eq.1)then | |
762 | c npom=nprt(1) | |
763 | c k=1 | |
764 | c ip=iproj(1) | |
765 | c it=itarg(1) | |
766 | c zparpro(k)=max(0,npom-1)*0.2 | |
767 | c zpartar(k)=0 | |
768 | c zpartar(k)=max(0,npom-1)*0.2 | |
769 | c ztav=zpartar(k) | |
770 | c zpav=zparpro(k) | |
771 | c zppevt=zpav | |
772 | c zptevt=ztav | |
773 | c endif | |
774 | end | |
775 | ||
776 | c---------------------------------------------------------------------- | |
777 | double precision function om1(xh,yp,b) !---test--- | |
778 | c---------------------------------------------------------------------- | |
779 | c om1 = G * C * gamd (C and gamd usually 1) | |
780 | c xh - fraction of the energy squared s for the pomeron; | |
781 | c b - impact parameter between the pomeron ends; | |
782 | c yp - rapidity for the pomeron; | |
783 | c---------------------------------------------------------------------- | |
784 | include 'epos.inc' | |
785 | include 'epos.incsem' | |
786 | include 'epos.incpar' | |
787 | ||
788 | double precision Gf,xp,xm,xh,yp | |
789 | ||
790 | Gf=0.d0 | |
791 | xp=sqrt(xh)*exp(yp) | |
792 | xm=xh/xp | |
793 | spp=engy**2 | |
794 | imax=idxD1 | |
795 | if(iomega.eq.2)imax=1 | |
796 | do i=idxDmin,imax | |
797 | call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
798 | Gf=Gf+alp*real(xp)**bet*real(xm)**betp | |
799 | enddo | |
800 | om1=Gf | |
801 | * * dble(chad(iclpro)*chad(icltar)) | |
802 | end | |
803 | ||
804 | c---------------------------------------------------------------------- | |
805 | double precision function om1intb(b) !---MC--- | |
806 | c---------------------------------------------------------------------- | |
807 | c om1 integrated over xp and xm for given b | |
808 | c Calculation by analytical integration | |
809 | c---------------------------------------------------------------------- | |
810 | include 'epos.inc' | |
811 | include 'epos.incsem' | |
812 | include 'epos.incpar' | |
813 | parameter(eps=1.e-20) | |
814 | ||
815 | spp=engy*engy | |
816 | imax=idxD1 | |
817 | if(iomega.eq.2)imax=1 | |
818 | cint=0 | |
819 | do i=idxDmin,imax | |
820 | call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
821 | cint2=gamv*alp | |
822 | if((bet+1.0).gt.eps)then | |
823 | cint2=cint2/(bet+1.0) | |
824 | else | |
825 | cint2=-cint2*log(xminDf) | |
826 | endif | |
827 | if((betp+1.0).gt.eps)then | |
828 | cint2=cint2/(betp+1.0) | |
829 | else | |
830 | cint2=-cint2*log(xminDf) | |
831 | endif | |
832 | cint=cint+cint2 | |
833 | enddo | |
834 | ||
835 | if(cint.lt.-eps)then | |
836 | write(*,*) 'WARNING ! om1intb in epos-omg is <0 !!!!!' | |
837 | write(*,*) 'WARNING ! => om1intb set to 1e-3 !!!!!' | |
838 | write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!' | |
839 | cint=1.e-3 | |
840 | endif | |
841 | ||
842 | om1intb=cint | |
843 | * *chad(iclpro)*chad(icltar) | |
844 | ||
845 | return | |
846 | end | |
847 | ||
848 | c---------------------------------------------------------------------- | |
849 | double precision function om1intbi(b,iqq) !---MC--- | |
850 | c---------------------------------------------------------------------- | |
851 | c om1 integrated over xp and xm for given b | |
852 | c Calculation by analytical integration of contribution iqq | |
853 | c---------------------------------------------------------------------- | |
854 | include 'epos.inc' | |
855 | include 'epos.incsem' | |
856 | include 'epos.incpar' | |
857 | parameter(eps=1.e-20) | |
858 | ||
859 | spp=engy*engy | |
860 | call GfunPar(1,iqq,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
861 | cint=gamv*alp | |
862 | if((bet+1.0).gt.eps)then | |
863 | cint=cint/(bet+1.0) | |
864 | else | |
865 | cint=-cint*log(xminDf) | |
866 | endif | |
867 | if((betp+1.0).gt.eps)then | |
868 | cint=cint/(betp+1.0) | |
869 | else | |
870 | cint=-cint*log(xminDf) | |
871 | endif | |
872 | if(cint.lt.-eps)then | |
873 | write(*,*) 'WARNING ! om1intbi in epos-omg is <0 !!!!!' | |
874 | write(*,*) 'WARNING ! => om1intbi set to 1e-3 !!!!!' | |
875 | write(*,*) 'WARNING ! => bmax=3.5 fm !!!!!' | |
876 | cint=1.e-3 | |
877 | endif | |
878 | ||
879 | om1intbi=cint | |
880 | * *chad(iclpro)*chad(icltar) | |
881 | ||
882 | return | |
883 | end | |
884 | ||
885 | c---------------------------------------------------------------------- | |
886 | double precision function om1intbc(b) !---MC--- | |
887 | c---------------------------------------------------------------------- | |
888 | c om1*F*F integrated over xp and xm for given b | |
889 | c Calculation by analytical integration | |
890 | c---------------------------------------------------------------------- | |
891 | include 'epos.inc' | |
892 | include 'epos.incems' | |
893 | include 'epos.incsem' | |
894 | double precision cint,gamom,deltap,deltam | |
895 | double precision utgam2,Fp,Fm | |
896 | ||
897 | spp=engy**2 | |
898 | om1intbc=0.d0 | |
899 | Fp=dble(ucfpro) !gamma(1+alplea) | |
900 | Fm=dble(ucftar) | |
901 | ||
902 | imax=idxD1 | |
903 | if(iomega.eq.2)imax=1 | |
904 | ||
905 | cint=0.d0 | |
906 | ||
907 | do i=idxDmin,imax | |
908 | call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
909 | gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar) | |
910 | deltap=dble(bet) | |
911 | deltam=dble(betp) | |
912 | cint=cint+gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0) | |
913 | & /utgam2(2.d0+deltap+dble(alplea(iclpro))) | |
914 | & /utgam2(2.d0+deltam+dble(alplea(icltar))) | |
915 | enddo | |
916 | ||
917 | om1intbc=cint*Fp*Fm | |
918 | ||
919 | if(om1intbc.lt.-1.d-10)then | |
920 | write(*,*) 'WARNING ! om1intbc in epos-omg is <0 !!!!!' | |
921 | write(*,*) 'WARNING ! => om1intbc set to 0. !!!!!' | |
922 | om1intbc=0.d0 | |
923 | endif | |
924 | ||
925 | return | |
926 | end | |
927 | ||
928 | cc---------------------------------------------------------------------- | |
929 | c double precision function om1intbci(b,iqq) !---MC--- | |
930 | cc---------------------------------------------------------------------- | |
931 | cc om1*F*F integrated over xp and xm for given b and given Pomeron type iqq | |
932 | cc Calculation by analytical integration | |
933 | cc---------------------------------------------------------------------- | |
934 | c include 'epos.inc' | |
935 | c include 'epos.incems' | |
936 | c include 'epos.incsem' | |
937 | c double precision cint,gamom,deltap,deltam | |
938 | c double precision utgam2,Fp,Fm,eps | |
939 | c | |
940 | c spp=engy**2 | |
941 | c om1intbci=0.d0 | |
942 | c Fp=dble(ucfpro) !gamma(1+alplea) | |
943 | c Fm=dble(ucftar) | |
944 | c | |
945 | c i=iqq | |
946 | c call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
947 | c gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar) | |
948 | c deltap=dble(bet) | |
949 | c deltam=dble(betp) | |
950 | c cint=gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0) | |
951 | c & /utgam2(2.d0+deltap+dble(alplea(iclpro))) | |
952 | c & /utgam2(2.d0+deltam+dble(alplea(icltar))) | |
953 | c | |
954 | c om1intbci=cint*Fp*Fm | |
955 | c | |
956 | c if(om1intbci.lt.-1.d-10)then | |
957 | c write(*,*) 'WARNING ! om1intbci in epos-omg is <0 !!!!!' | |
958 | c write(*,*) 'WARNING ! => om1intbci set to 0. !!!!!' | |
959 | c om1intbci=0.d0 | |
960 | c endif | |
961 | c | |
962 | c return | |
963 | c end | |
964 | c | |
965 | c---------------------------------------------------------------------- | |
966 | double precision function om1intgck(k,xprem,xmrem) !---MC--- | |
967 | c---------------------------------------------------------------------- | |
968 | c om1*(xprem-xp)*(xmrem-xm) integrated over xp and xm for given k | |
969 | c Calculation by analytical integration | |
970 | c---------------------------------------------------------------------- | |
971 | include 'epos.inc' | |
972 | include 'epos.incems' | |
973 | include 'epos.incsem' | |
974 | double precision cint,gamom,deltap,deltam,xprem,xmrem | |
975 | ||
976 | om1intgck=0.d0 | |
977 | ||
978 | imax=idxD1 | |
979 | if(iomega.eq.2)imax=1 | |
980 | ||
981 | cint=0.d0 | |
982 | do i=idxDmin,imax | |
983 | gamom=dble(atilde(i,k)) | |
984 | deltap=dble(btildep(i,k)) | |
985 | deltam=dble(btildepp(i,k)) | |
986 | cint=cint+gamom/(deltap+1.d0)/(deltam+1.d0) | |
987 | & /(2.d0+deltap) /(2.d0+deltam) | |
988 | & *xprem**(deltap+2.d0) | |
989 | & *xmrem**(deltam+2.d0) | |
990 | enddo | |
991 | om1intgck=cint | |
992 | ||
993 | return | |
994 | end | |
995 | ||
996 | c---------------------------------------------------------------------- | |
997 | double precision function om1intgc(b) !---test--- | |
998 | c---------------------------------------------------------------------- | |
999 | c om1*(1-xp)*(1-xm) integrated over xp and xm for given b | |
1000 | c Calculation by analytical integration | |
1001 | c---------------------------------------------------------------------- | |
1002 | include 'epos.inc' | |
1003 | include 'epos.incpar' | |
1004 | include 'epos.incsem' | |
1005 | double precision cint,gamom,deltap,deltam,eps | |
1006 | parameter(eps=1.d-20) | |
1007 | ||
1008 | spp=engy**2 | |
1009 | om1intgc=0.d0 | |
1010 | ||
1011 | imax=idxD1 | |
1012 | if(iomega.eq.2)imax=1 | |
1013 | ||
1014 | ||
1015 | ||
1016 | cint=0.d0 | |
1017 | ||
1018 | do i=idxDmin,imax | |
1019 | call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv) | |
1020 | gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar) | |
1021 | deltap=dble(bet) | |
1022 | deltam=dble(betp) | |
1023 | if((deltap+1.d0).gt.eps)then | |
1024 | gamom=gamom/(deltap+1.0) | |
1025 | else | |
1026 | gamom=-gamom*log(xminDf) | |
1027 | endif | |
1028 | if((deltam+1.d0).gt.eps)then | |
1029 | gamom=gamom/(deltam+1.0) | |
1030 | else | |
1031 | gamom=-gamom*log(xminDf) | |
1032 | endif | |
1033 | cint=cint+gamom /(2.d0+deltap) /(2.d0+deltam) | |
1034 | enddo | |
1035 | om1intgc=cint | |
1036 | ||
1037 | if(om1intgc.lt.eps)then | |
1038 | write(*,*) b,deltap,deltam,gamom | |
1039 | write(*,*) 'WARNING ! om1intgc in epos-omg is <0 !!!!!' | |
1040 | write(*,*) 'WARNING ! => om1intgc set to 0. !!!!!' | |
1041 | om1intgc=0.d0 | |
1042 | endif | |
1043 | ||
1044 | return | |
1045 | end | |
1046 | ||
1047 | ||
1048 | c---------------------------------------------------------------------- | |
1049 | subroutine integom1(irea) | |
1050 | c---------------------------------------------------------------------- | |
1051 | c om1 integrated over xp and xm for all b=bk(k) if irea=0 | |
1052 | c result written to : | |
1053 | c om1int(k) = om1intb(bk(k)) = \int om1 | |
1054 | c om1intc(k) = om1intgc(bk(k)... = \int om1*(1-xp)*(1-xm) | |
1055 | c---------------------------------------------------------------------- | |
1056 | include 'epos.inc' | |
1057 | include 'epos.incems' | |
1058 | include 'epos.incsem' | |
1059 | double precision om1intb,om1intgc | |
1060 | ||
1061 | call utpri('intom1',ish,ishini,5) | |
1062 | ||
1063 | if(irea.eq.0)then | |
1064 | ||
1065 | ||
1066 | do k=1,koll | |
1067 | ||
1068 | om1int(k)=om1intb(bk(k)) | |
1069 | om1intc(k)=om1intgc(bk(k)) | |
1070 | ||
1071 | enddo | |
1072 | ||
1073 | endif | |
1074 | call utprix('intom1',ish,ishini,5) | |
1075 | ||
1076 | ||
1077 | return | |
1078 | end | |
1079 | ||
1080 | ||
1081 | c---------------------------------------------------------------------- | |
1082 | double precision function om1xpk(xp,xpr1i,xmr1i,k) !---test--- | |
1083 | c---------------------------------------------------------------------- | |
1084 | c \int dxm om1*(1-xp)*(1-xm) (normalized) | |
1085 | c k - pair indice; | |
1086 | c---------------------------------------------------------------------- | |
1087 | include 'epos.inc' | |
1088 | include 'epos.incpar' | |
1089 | include 'epos.incems' | |
1090 | double precision xp,gamomx(ntymi:ntymx),cint,gamom | |
1091 | double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps | |
1092 | & ,xpr1,xmr1,xpr1i,xmr1i | |
1093 | parameter(eps=1.d-20) | |
1094 | ||
1095 | om1xpk=0.d0 | |
1096 | if(xp.ge.xpr1i)return | |
1097 | ||
1098 | xpr1=1.d0 | |
1099 | xmr1=1.d0 | |
1100 | imin=ntymin | |
1101 | imax=ntymx | |
1102 | if(iomega.eq.2)imax=1 | |
1103 | ||
1104 | do i=imin,imax | |
1105 | deltap(i)=btildep(i,k) | |
1106 | deltam(i)=btildepp(i,k) | |
1107 | gamomx(i)=atilde(i,k)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i)) | |
1108 | if((deltam(i)+1.d0).gt.eps)then | |
1109 | gamomx(i)=gamomx(i)/(deltam(i)+1.0) | |
1110 | else | |
1111 | gamomx(i)=-gamomx(i)*log(xminDf) | |
1112 | endif | |
1113 | enddo | |
1114 | ||
1115 | ||
1116 | cint=0.d0 | |
1117 | do i=imin,imax | |
1118 | gamom=gamomx(i)*xpr1**(deltap(i)+2.d0)/(2.d0+deltap(i)) | |
1119 | if((deltap(i)+1.d0).gt.eps)then | |
1120 | gamom=gamom/(deltap(i)+1.d0) | |
1121 | else | |
1122 | gamom=-gamom*log(xminDf) | |
1123 | endif | |
1124 | cint=cint+gamom | |
1125 | enddo | |
1126 | ||
1127 | ||
1128 | do i=imin,imax | |
1129 | om1xpk=om1xpk+gamomx(i)*xp**deltap(i) | |
1130 | & *(xpr1-xp) | |
1131 | ||
1132 | enddo | |
1133 | ||
1134 | om1xpk=om1xpk/cint | |
1135 | ||
1136 | return | |
1137 | end | |
1138 | ||
1139 | ||
1140 | c---------------------------------------------------------------------- | |
1141 | double precision function om1xmk(xp,xm,xpr1i,xmr1i,k) !---test--- | |
1142 | c---------------------------------------------------------------------- | |
1143 | c om1(xp,xm)*(1-xp)*(1-xm) (normalized for fixed xp) | |
1144 | c k - pair indice; | |
1145 | c---------------------------------------------------------------------- | |
1146 | include 'epos.inc' | |
1147 | include 'epos.incpar' | |
1148 | include 'epos.incems' | |
1149 | double precision xp,xm,gamomx(ntymi:ntymx),cint,gamom | |
1150 | double precision deltam(ntymi:ntymx),eps,xpr1,xmr1,xpr1i,xmr1i | |
1151 | parameter(eps=1.d-20) | |
1152 | ||
1153 | om1xmk=0.d0 | |
1154 | if(xp.ge.xpr1i)return | |
1155 | if(xm.ge.xmr1i)return | |
1156 | xpr1=1.d0 | |
1157 | xmr1=1.d0 | |
1158 | ||
1159 | imin=ntymin | |
1160 | imax=ntymx | |
1161 | if(iomega.eq.2)imax=1 | |
1162 | ||
1163 | do i=imin,imax | |
1164 | gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xpr1-xp) | |
1165 | deltam(i)=btildepp(i,k) | |
1166 | enddo | |
1167 | ||
1168 | cint=0.d0 | |
1169 | do i=imin,imax | |
1170 | gamom=gamomx(i)*xmr1**(deltam(i)+2.d0)/(2.d0+deltam(i)) | |
1171 | if((deltam(i)+1.d0).gt.eps)then | |
1172 | gamom=gamom/(deltam(i)+1.0) | |
1173 | else | |
1174 | gamom=-gamom*log(xminDf) | |
1175 | endif | |
1176 | cint=cint+gamom | |
1177 | enddo | |
1178 | ||
1179 | ||
1180 | do i=imin,imax | |
1181 | om1xmk=om1xmk+gamomx(i)*xm**deltam(i)*(xmr1-xm) | |
1182 | enddo | |
1183 | ||
1184 | om1xmk=om1xmk/cint | |
1185 | ||
1186 | return | |
1187 | end | |
1188 | ||
1189 | c---------------------------------------------------------------------- | |
1190 | double precision function om1xprk(k,xpremi,xmremi,ir) !---MC--- | |
1191 | c---------------------------------------------------------------------- | |
1192 | c Random number generated from the function om1xpk. We solve the equation | |
1193 | c which give om1xprk by Newton-Raphson + secant method. | |
1194 | c k - pair indice; | |
1195 | c ir - 1 to get xp, -1 to get xm (be carrefull to inverse xpremi et xmremi | |
1196 | c when calling with ir=-1) | |
1197 | c---------------------------------------------------------------------- | |
1198 | include 'epos.inc' | |
1199 | include 'epos.incpar' | |
1200 | include 'epos.incems' | |
1201 | double precision x,x0,x1,gamomx(ntymi:ntymx),eps,f0t,f1t,f00 | |
1202 | double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,drangen,xmrem | |
1203 | double precision deltap(ntymi:ntymx),deltam(ntymi:ntymx),xprem | |
1204 | parameter (eps=1.d-20) | |
1205 | double precision xpremi,xmremi,lxmin | |
1206 | ||
1207 | ||
1208 | om1xprk=0.d0 | |
1209 | x0=log(xminDf) | |
1210 | x1=log(xpremi) | |
1211 | xprem=1.d0 | |
1212 | xmrem=1.d0 | |
1213 | imin=ntymin | |
1214 | imax=ntymx | |
1215 | lxmin=1.d0 | |
1216 | if(iomega.eq.2)imax=1 | |
1217 | ||
1218 | do i=imin,imax | |
1219 | if(ir.gt.0)then | |
1220 | deltap(i)=btildep(i,k) | |
1221 | deltam(i)=btildepp(i,k) | |
1222 | else | |
1223 | deltap(i)=btildepp(i,k) | |
1224 | deltam(i)=btildep(i,k) | |
1225 | endif | |
1226 | gamomx(i)=atilde(i,k)*xmrem**(deltam(i)+2.d0)/(2.d0+deltam(i)) | |
1227 | if((deltam(i)+1.d0).gt.eps)then | |
1228 | gamomx(i)=gamomx(i)/(deltam(i)+1.0) | |
1229 | else | |
1230 | lxmin=log(xminDf) | |
1231 | gamomx(i)=-gamomx(i)*lxmin | |
1232 | endif | |
1233 | enddo | |
1234 | ||
1235 | f0=0.d0 | |
1236 | f1=0.d0 | |
1237 | do i=imin,imax | |
1238 | if((deltap(i)+1.d0).gt.eps)then | |
1239 | f0=f0+gamomx(i) | |
1240 | & *(xprem*exp(x0)**(1.d0+deltap(i))/(1.d0+deltap(i)) | |
1241 | & -exp(x0)**(2.d0+deltap(i))/(2.d0+deltap(i))) | |
1242 | f1=f1+gamomx(i) | |
1243 | & *(xprem*exp(x1)**(1.d0+deltap(i))/(1.d0+deltap(i)) | |
1244 | & -exp(x1)**(2.d0+deltap(i))/(2.d0+deltap(i))) | |
1245 | else | |
1246 | lxmin=log(xminDf) | |
1247 | f0=f0+gamomx(i)*(xprem*(x0-lxmin)-exp(x0)+xminDf) | |
1248 | f1=f1+gamomx(i)*(xprem*(x1-lxmin)-exp(x1)+xminDf) | |
1249 | endif | |
1250 | enddo | |
1251 | f00=f0 | |
1252 | cint=f1-f00 | |
1253 | f0=-(f0-f00)/cint | |
1254 | f1=-(f1-f00)/cint | |
1255 | ntry=0 | |
1256 | 11 ntry=ntry+1 | |
1257 | r=drangen(dble(ntry)) | |
1258 | f0t=f0+r | |
1259 | f1t=f1+r | |
1260 | if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11 | |
1261 | if(f1t*f0t.ge.eps)then | |
1262 | do i=imin,imax | |
1263 | write(ifmt,*)i,gamomx(i),deltap(i),deltam(i) | |
1264 | enddo | |
1265 | write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k | |
1266 | call utstop('om1xprk (2)&') | |
1267 | endif | |
1268 | f0=f0t | |
1269 | f1=f1t | |
1270 | if(abs(f0).le.eps) then | |
1271 | om1xprk=exp(x0) | |
1272 | return | |
1273 | endif | |
1274 | if(abs(f1).le.eps) then | |
1275 | om1xprk=exp(x1) | |
1276 | return | |
1277 | endif | |
1278 | x=0.5d0*(x1+x0) | |
1279 | deltx=abs(x1-x0) | |
1280 | ||
1281 | ||
1282 | ntry=0 | |
1283 | ||
1284 | 111 continue | |
1285 | if(ntry.le.1000)then | |
1286 | fx=0.d0 | |
1287 | fpx=0.d0 | |
1288 | do i=imin,imax | |
1289 | if((deltap(i)+1.d0).gt.eps)then | |
1290 | fx=fx+gamomx(i) | |
1291 | & *(xprem*exp(x)**(1.d0+deltap(i))/(1.d0+deltap(i)) | |
1292 | & -exp(x)**(2.d0+deltap(i))/(2.d0+deltap(i))) | |
1293 | fpx=fpx+gamomx(i)*exp(x)**deltap(i)*(xprem-exp(x)) | |
1294 | else | |
1295 | fx=fx+gamomx(i)*(xprem*(x-lxmin)-exp(x)+xminDf) | |
1296 | fpx=fpx+gamomx(i)*(xprem/exp(x)-1.d0) | |
1297 | endif | |
1298 | enddo | |
1299 | fx=-(fx-f00)/cint+r | |
1300 | fpx=fpx/cint | |
1301 | xt=x-fx/fpx | |
1302 | ||
1303 | if (f0*fx.lt.0.D0) then | |
1304 | f1=fx | |
1305 | x1=x | |
1306 | else | |
1307 | f0=fx | |
1308 | x0=x | |
1309 | endif | |
1310 | if ((xt.lt.x0.or.xt.gt.x1).and.abs(f1-f0).gt.eps) then | |
1311 | xt=x1-f1*(x1-x0)/(f1-f0) | |
1312 | endif | |
1313 | ||
1314 | else | |
1315 | ||
1316 | write(ifmt,*)'Warning in om1xprk, to much try !' | |
1317 | ||
1318 | endif | |
1319 | ||
1320 | ||
1321 | if(abs(x-xt).gt.deltx*0.5d0) then | |
1322 | xt=(x1+x0)*0.5D0 | |
1323 | endif | |
1324 | deltx=abs(x-xt) | |
1325 | if(abs(x).gt.eps)then | |
1326 | prec=abs(deltx/x) | |
1327 | else | |
1328 | call utstop('Problem in om1xprk&') | |
1329 | endif | |
1330 | ||
1331 | if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then | |
1332 | x=xt | |
1333 | ntry=ntry+1 | |
1334 | goto 111 | |
1335 | endif | |
1336 | ||
1337 | om1xprk=exp(x) | |
1338 | ||
1339 | return | |
1340 | end | |
1341 | ||
1342 | c---------------------------------------------------------------------- | |
1343 | double precision function om1xmrk(k,xp,xpremi,xmremi,ir) !---MC--- | |
1344 | c---------------------------------------------------------------------- | |
1345 | c Random number generated from the function om1xmk. We solve the equation | |
1346 | c which give om1xmrk by Newton-Raphson + secant method. | |
1347 | c k - pair indice; | |
1348 | c ir - 1 to get xm, -1 to get xp (be carrefull to inverse xpremi et xmremi | |
1349 | c when calling with ir=-1) | |
1350 | c---------------------------------------------------------------------- | |
1351 | include 'epos.inc' | |
1352 | include 'epos.incpar' | |
1353 | include 'epos.incems' | |
1354 | double precision x,x0,x1,gamomx(ntymi:ntymx),eps,xp,f0t,f1t | |
1355 | double precision xt,fx,fpx,r,f1,f0,cint,deltx,prec,f00 | |
1356 | double precision deltam(ntymi:ntymx),drangen,xprem,xmrem | |
1357 | double precision xpremi,xmremi | |
1358 | parameter (eps=1.d-20) | |
1359 | ||
1360 | om1xmrk=0.d0 | |
1361 | if(xp.ge.xpremi)return | |
1362 | ||
1363 | xprem=1.d0 | |
1364 | xmrem=1.d0 | |
1365 | x0=log(xminDf) | |
1366 | x1=log(xmremi) | |
1367 | imin=ntymin | |
1368 | imax=ntymx | |
1369 | if(iomega.eq.2)imax=1 | |
1370 | ||
1371 | do i=imin,imax | |
1372 | if(ir.gt.0)then | |
1373 | gamomx(i)=atilde(i,k)*xp**btildep(i,k)*(xprem-xp) | |
1374 | deltam(i)=btildepp(i,k) | |
1375 | else | |
1376 | gamomx(i)=atilde(i,k)*xp**btildepp(i,k)*(xprem-xp) | |
1377 | deltam(i)=btildep(i,k) | |
1378 | endif | |
1379 | enddo | |
1380 | ||
1381 | ||
1382 | ||
1383 | f0=0.d0 | |
1384 | f1=0.d0 | |
1385 | do i=imin,imax | |
1386 | if(abs(deltam(i)+1.d0).gt.eps)then | |
1387 | f0=f0+gamomx(i) | |
1388 | & *(xmrem*exp(x0)**(1.d0+deltam(i))/(1.d0+deltam(i)) | |
1389 | & -exp(x0)**(2.d0+deltam(i))/(2.d0+deltam(i))) | |
1390 | f1=f1+gamomx(i) | |
1391 | & *(xmrem*exp(x1)**(1.d0+deltam(i))/(1.d0+deltam(i)) | |
1392 | & -exp(x1)**(2.d0+deltam(i))/(2.d0+deltam(i))) | |
1393 | else | |
1394 | lxmin=log(xminDf) | |
1395 | f0=f0+gamomx(i)*(xmrem*(x0-lxmin)-exp(x0)+xminDf) | |
1396 | f1=f1+gamomx(i)*(xmrem*(x1-lxmin)-exp(x1)+xminDf) | |
1397 | endif | |
1398 | enddo | |
1399 | f00=f0 | |
1400 | cint=f1-f00 | |
1401 | f0=-(f0-f00)/cint | |
1402 | f1=-(f1-f00)/cint | |
1403 | ntry=0 | |
1404 | 11 ntry=ntry+1 | |
1405 | r=drangen(dble(ntry)) | |
1406 | f0t=f0+r | |
1407 | f1t=f1+r | |
1408 | if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11 | |
1409 | if(f1t*f0t.ge.eps)then | |
1410 | write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry | |
1411 | call utstop('Error(2) in epos-omg in om1xmrk&') | |
1412 | endif | |
1413 | f0=f0t | |
1414 | f1=f1t | |
1415 | if(abs(f0).lt.eps) then | |
1416 | om1xmrk=exp(x0) | |
1417 | return | |
1418 | endif | |
1419 | if(abs(f1).lt.eps) then | |
1420 | om1xmrk=exp(x1) | |
1421 | return | |
1422 | endif | |
1423 | x=0.5d0*(x1+x0) | |
1424 | deltx=abs(x1-x0) | |
1425 | ||
1426 | ||
1427 | ntry=0 | |
1428 | ||
1429 | 111 continue | |
1430 | if(ntry.le.1000)then | |
1431 | fx=0.d0 | |
1432 | fpx=0.d0 | |
1433 | do i=imin,imax | |
1434 | if(abs(deltam(i)+1.d0).gt.eps)then | |
1435 | fx=fx+gamomx(i) | |
1436 | & *(xmrem*exp(x)**(1.d0+deltam(i))/(1.d0+deltam(i)) | |
1437 | & -exp(x)**(2.d0+deltam(i))/(2.d0+deltam(i))) | |
1438 | fpx=fpx+gamomx(i)*exp(x)**deltam(i)*(xmrem-exp(x)) | |
1439 | else | |
1440 | fx=fx+gamomx(i)*(xmrem*(x-lxmin)-exp(x)+xminDf) | |
1441 | fpx=fpx+gamomx(i)*(xmrem/exp(x)-1.d0) | |
1442 | endif | |
1443 | enddo | |
1444 | fx=-(fx-f00)/cint+r | |
1445 | fpx=fpx/cint | |
1446 | xt=x-fx/fpx | |
1447 | ||
1448 | if (f0*fx.lt.-eps) then | |
1449 | f1=fx | |
1450 | x1=x | |
1451 | else | |
1452 | f0=fx | |
1453 | x0=x | |
1454 | endif | |
1455 | if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then | |
1456 | xt=x1-f1*(x1-x0)/(f1-f0) | |
1457 | endif | |
1458 | ||
1459 | else | |
1460 | ||
1461 | write(ifmt,*)'Warning in om1xmrk, to much try !' | |
1462 | ||
1463 | endif | |
1464 | ||
1465 | if(abs(x-xt).gt.deltx*0.5d0) then | |
1466 | xt=(x1+x0)*0.5D0 | |
1467 | endif | |
1468 | deltx=abs(x-xt) | |
1469 | if(abs(x).gt.eps)then | |
1470 | prec=abs(deltx/x) | |
1471 | else | |
1472 | call utstop('Problem in om1xmrk&') | |
1473 | endif | |
1474 | ||
1475 | if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000) then | |
1476 | x=xt | |
1477 | ntry=ntry+1 | |
1478 | goto 111 | |
1479 | endif | |
1480 | ||
1481 | om1xmrk=exp(x) | |
1482 | ||
1483 | return | |
1484 | end | |
1485 | ||
1486 | c---------------------------------------------------------------------- | |
1487 | double precision function om1xk(xh,k) !---test--- | |
1488 | c---------------------------------------------------------------------- | |
1489 | c \int dxp om1 (normalised) | |
1490 | c k - pair indice; | |
1491 | c---------------------------------------------------------------------- | |
1492 | include 'epos.inc' | |
1493 | include 'epos.incpar' | |
1494 | include 'epos.incems' | |
1495 | ||
1496 | double precision xh,gamomx(ntymi:ntymx),cint,alpp(ntymi:ntymx) | |
1497 | &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx),eps | |
1498 | &,gamom | |
1499 | parameter(eps=1.d-20) | |
1500 | ||
1501 | ||
1502 | om1xk=0.d0 | |
1503 | ||
1504 | imin=ntymin | |
1505 | imax=ntymx | |
1506 | if(iomega.eq.2)imax=1 | |
1507 | ||
1508 | do i=imin,imax | |
1509 | gamomx(i)=atilde(i,k) | |
1510 | deltap(i)=btildep(i,k) | |
1511 | deltam(i)=btildepp(i,k) | |
1512 | ||
1513 | delta(i)=(deltap(i)+deltam(i))*0.5d0 | |
1514 | alpp(i)=deltap(i)-deltam(i) | |
1515 | ||
1516 | enddo | |
1517 | ||
1518 | cint=0.d0 | |
1519 | do i=imin,imax | |
1520 | gamom=gamomx(i) | |
1521 | if((deltap(i)+1.d0).gt.eps)then | |
1522 | gamom=gamom/(deltap(i)+1.d0) | |
1523 | else | |
1524 | gamom=-gamom*log(xminDf) | |
1525 | endif | |
1526 | if((deltam(i)+1.d0).gt.eps)then | |
1527 | gamom=gamom/(deltam(i)+1.d0) | |
1528 | else | |
1529 | gamom=-gamom*log(xminDf) | |
1530 | endif | |
1531 | cint=cint+gamom | |
1532 | enddo | |
1533 | ||
1534 | ||
1535 | do i=imin,imax | |
1536 | if(abs(alpp(i)).gt.eps)then | |
1537 | om1xk=om1xk+gamomx(i)/alpp(i)*xh**deltam(i)*(1.d0-xh**alpp(i)) | |
1538 | else | |
1539 | om1xk=om1xk-gamomx(i)*xh**delta(i)*log(xh) | |
1540 | endif | |
1541 | enddo | |
1542 | ||
1543 | om1xk=om1xk/cint | |
1544 | ||
1545 | return | |
1546 | end | |
1547 | ||
1548 | c---------------------------------------------------------------------- | |
1549 | double precision function om1yk(xh,yp,k) !---test--- | |
1550 | c---------------------------------------------------------------------- | |
1551 | c om1 normalized for fixed xp | |
1552 | c xh - fraction of the energy squared s for the pomeron; | |
1553 | c k - pair indice; | |
1554 | c---------------------------------------------------------------------- | |
1555 | include 'epos.inc' | |
1556 | include 'epos.incems' | |
1557 | double precision xh,yp,gamomy(ntymi:ntymx),alpp(ntymi:ntymx),cint | |
1558 | double precision deltap,deltam,eps | |
1559 | parameter(eps=1.d-20) | |
1560 | ||
1561 | om1yk=0.d0 | |
1562 | ||
1563 | imin=ntymin | |
1564 | imax=ntymx | |
1565 | if(iomega.eq.2)imax=1 | |
1566 | ||
1567 | do i=imin,imax | |
1568 | gamomy(i)=atilde(i,k) | |
1569 | deltap=btildep(i,k) | |
1570 | deltam=btildepp(i,k) | |
1571 | ||
1572 | alpp(i)=deltap-deltam | |
1573 | gamomy(i)=gamomy(i)*xh**((deltap+deltam)*0.5d0) | |
1574 | ||
1575 | enddo | |
1576 | ||
1577 | cint=0.d0 | |
1578 | do i=imin,imax | |
1579 | if(abs(alpp(i)).gt.eps)then | |
1580 | cint=cint-gamomy(i)/alpp(i)*xh**(alpp(i)*0.5d0) | |
1581 | & *(1.d0-xh**(-alpp(i))) | |
1582 | else | |
1583 | cint=cint-gamomy(i)*log(xh) | |
1584 | endif | |
1585 | enddo | |
1586 | ||
1587 | do i=imin,imax | |
1588 | if(abs(alpp(i)).gt.eps)then | |
1589 | om1yk=om1yk+gamomy(i)*exp(alpp(i)*yp) | |
1590 | else | |
1591 | om1yk=om1yk+gamomy(i) | |
1592 | endif | |
1593 | enddo | |
1594 | ||
1595 | om1yk=om1yk/cint | |
1596 | ||
1597 | ||
1598 | return | |
1599 | end | |
1600 | ||
1601 | c---------------------------------------------------------------------- | |
1602 | double precision function om1xrk(k) !---test--- | |
1603 | c---------------------------------------------------------------------- | |
1604 | c Random number generated from the function om1xk. We solve the equation | |
1605 | c which give om1xrk by Newton-Raphson + secant method. | |
1606 | c k - pair indice; | |
1607 | c---------------------------------------------------------------------- | |
1608 | include 'epos.inc' | |
1609 | include 'epos.incems' | |
1610 | include 'epos.incpar' | |
1611 | ||
1612 | double precision x,x0,x1,gamomx(ntymi:ntymx),eps,prec,drangen | |
1613 | double precision xt,fx,fpx,r,f1,f0,cint,deltx,alpp(ntymi:ntymx) | |
1614 | &,delta(ntymi:ntymx),deltap(ntymi:ntymx),deltam(ntymi:ntymx) | |
1615 | parameter (eps=1.d-20) | |
1616 | ||
1617 | ||
1618 | om1xrk=0.d0 | |
1619 | ||
1620 | imin=ntymin | |
1621 | imax=ntymx | |
1622 | if(iomega.eq.2)imax=1 | |
1623 | ||
1624 | do i=imin,imax | |
1625 | gamomx(i)=atilde(i,k) | |
1626 | deltap(i)=btildep(i,k) | |
1627 | deltam(i)=btildepp(i,k) | |
1628 | ||
1629 | delta(i)=(deltap(i)+deltam(i))*0.5d0 | |
1630 | alpp(i)=deltap(i)-deltam(i) | |
1631 | ||
1632 | enddo | |
1633 | ||
1634 | cint=0.d0 | |
1635 | do i=imin,imax | |
1636 | gamom=gamomx(i) | |
1637 | if((deltap(i)+1.d0).gt.eps)then | |
1638 | gamom=gamom/(deltap(i)+1.d0) | |
1639 | else | |
1640 | gamom=-gamom*log(xminDf) | |
1641 | endif | |
1642 | if((deltam(i)+1.d0).gt.eps)then | |
1643 | gamom=gamom/(deltam(i)+1.d0) | |
1644 | else | |
1645 | gamom=-gamom*log(xminDf) | |
1646 | endif | |
1647 | cint=cint+gamom | |
1648 | enddo | |
1649 | ||
1650 | x0=eps | |
1651 | x1=1.d0 | |
1652 | f0=0.d0 | |
1653 | f1=0.d0 | |
1654 | do i=imin,imax | |
1655 | ||
1656 | if(abs(alpp(i)).lt.eps)then | |
1657 | if(delta(i)+1.d0.gt.eps)then | |
1658 | f0=f0-gamomx(i)/(delta(i)+1.d0)*x0**(delta(i)+1.d0) | |
1659 | & *(log(x0)-1.d0/(delta(i)+1.d0)) | |
1660 | f1=f1-gamomx(i)/(delta(i)+1.d0)*x1**(delta(i)+1.d0) | |
1661 | & *(log(x1)-1.d0/(delta(i)+1.d0)) | |
1662 | else | |
1663 | f0=f0-0.5d0*gamomx(i)*(log(x0)**2-log(xminDf)**2) | |
1664 | f1=f1-0.5d0*gamomx(i)*(log(x1)**2-log(xminDf)**2) | |
1665 | endif | |
1666 | else | |
1667 | if(abs(deltap(i)+1.d0).gt.eps | |
1668 | & .and.abs(deltam(i)+1.d0).gt.eps)then | |
1669 | f0=f0+gamomx(i)/alpp(i)*(x0**(deltam(i)+1.d0)/(deltam(i)+1.d0) | |
1670 | & -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1671 | f1=f1+gamomx(i)/alpp(i)*(x1**(deltam(i)+1.d0)/(deltam(i)+1.d0) | |
1672 | & -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1673 | elseif(abs(deltap(i)+1.d0).gt.eps)then | |
1674 | f0=f0+gamomx(i)/alpp(i)*(log(x0/xminDf) | |
1675 | & -x0**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1676 | f1=f1+gamomx(i)/alpp(i)*(log(x1/xminDf) | |
1677 | & -x1**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1678 | elseif(abs(deltam(i)+1.d0).gt.eps)then | |
1679 | f0=f0-gamomx(i)/alpp(i)*(log(x0/xminDf) | |
1680 | & -x0**(deltam(i)+1.d0)/(deltam(i)+1.d0)) | |
1681 | f1=f1-gamomx(i)/alpp(i)*(log(x1/xminDf) | |
1682 | & -x1**(deltam(i)+1.d0)/(deltam(i)+1.d0)) | |
1683 | endif | |
1684 | endif | |
1685 | enddo | |
1686 | f0=-f0/cint | |
1687 | f1=-f1/cint | |
1688 | ntry=0 | |
1689 | 11 ntry=ntry+1 | |
1690 | r=drangen(dble(ntry)) | |
1691 | f0t=f0+r | |
1692 | f1t=f1+r | |
1693 | if(f1t*f0t.ge.eps.and.ntry.lt.100)goto 11 | |
1694 | if(f1t*f0t.ge.eps)then | |
1695 | do i=imin,imax | |
1696 | write(ifmt,*)i,gamomx(i),deltap(i),deltam(i),alpp(i),delta(i) | |
1697 | enddo | |
1698 | write(ifmt,*)x0,f0,f0t,x1,f1,f1t,r,cint,ntry,bk(k),k | |
1699 | call utstop('om1xrk (1)&') | |
1700 | endif | |
1701 | f0=f0t | |
1702 | f1=f1t | |
1703 | c if(f1*f0.gt.eps)then | |
1704 | c call utmsg('om1xrk') | |
1705 | c write(ifch,*)'Poblem with x0, no root ! --> om1xrk=xminDf' | |
1706 | c write(ifmt,*)'Poblem with x0, no root ! --> om1xrk=xminDf' | |
1707 | c write(ifmt,*)f0,f1,cint,r | |
1708 | c call utmsgf | |
1709 | c om1xrk=x0 | |
1710 | c return | |
1711 | c endif | |
1712 | if(abs(f0).lt.eps) then | |
1713 | om1xrk=x0 | |
1714 | return | |
1715 | endif | |
1716 | if(abs(f1).lt.eps) then | |
1717 | om1xrk=x1 | |
1718 | return | |
1719 | endif | |
1720 | c x=(x1+x0)*0.5D0 | |
1721 | x=sqrt(x1*x0) | |
1722 | deltx=abs(x1-x0) | |
1723 | ||
1724 | ntry=0 | |
1725 | fx=0.d0 | |
1726 | fpx=0.d0 | |
1727 | xt=x | |
1728 | 111 continue | |
1729 | ||
1730 | if(ntry.le.1000)then | |
1731 | fx=0.d0 | |
1732 | fpx=0.d0 | |
1733 | do i=imin,imax | |
1734 | if(abs(alpp(i)).lt.eps)then | |
1735 | if(delta(i)+1.d0.gt.eps)then | |
1736 | fx=fx-gamomx(i)/(delta(i)+1.d0)*x**(delta(i)+1.d0) | |
1737 | & *(log(x)-1.d0/(delta(i)+1.d0)) | |
1738 | fpx=fpx-gamomx(i)*x**delta(i)*log(x) | |
1739 | else | |
1740 | fx=fx-0.5d0*gamomx(i)*(log(x)**2-log(xminDf)**2) | |
1741 | fpx=fpx-gamomx(i)*log(x)/x | |
1742 | endif | |
1743 | else | |
1744 | if(abs(deltap(i)+1.d0).gt.eps | |
1745 | & .and.abs(deltam(i)+1.d0).gt.eps)then | |
1746 | fx=fx+gamomx(i)/alpp(i)*(x**(deltam(i)+1.d0)/(deltam(i)+1.d0) | |
1747 | & -x**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1748 | fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i)) | |
1749 | elseif(abs(deltap(i)+1.d0).gt.eps)then | |
1750 | fx=fx+gamomx(i)/alpp(i)*(log(x/xminDf) | |
1751 | & -x**(deltap(i)+1.d0)/(deltap(i)+1.d0)) | |
1752 | fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i)) | |
1753 | elseif(abs(deltam(i)+1.d0).gt.eps)then | |
1754 | fx=fx-gamomx(i)/alpp(i)*(log(x/xminDf) | |
1755 | & -x**(deltam(i)+1.d0)/(deltam(i)+1.d0)) | |
1756 | fpx=fpx+gamomx(i)/alpp(i)*x**deltam(i)*(1.d0-x**alpp(i)) | |
1757 | endif | |
1758 | endif | |
1759 | enddo | |
1760 | fx=-fx/cint+r | |
1761 | fpx=fpx/cint | |
1762 | xt=x-fx/fpx | |
1763 | ||
1764 | if (f0*fx.lt.-eps) then | |
1765 | f1=fx | |
1766 | x1=x | |
1767 | else | |
1768 | f0=fx | |
1769 | x0=x | |
1770 | endif | |
1771 | if ((xt.lt.x0-eps.or.xt.gt.x1+eps).and.abs(f1-f0).gt.eps) then | |
1772 | xt=x1-f1*(x1-x0)/(f1-f0) | |
1773 | endif | |
1774 | ||
1775 | else | |
1776 | ||
1777 | write(ifmt,*)'Warning in om1xrk, to much try !' | |
1778 | ||
1779 | endif | |
1780 | ||
1781 | if(abs(x-xt).gt.deltx*0.5d0) then | |
1782 | xt=sqrt(x1*x0) | |
1783 | endif | |
1784 | deltx=abs(x-xt) | |
1785 | if(abs(x).gt.eps)then | |
1786 | prec=deltx/x | |
1787 | else | |
1788 | call utstop('Problem in om1xrk&') | |
1789 | endif | |
1790 | ||
1791 | if (prec.gt.1.d-3.and.abs(f1-f0).gt.eps.and.ntry.le.1000)then | |
1792 | x=xt | |
1793 | ntry=ntry+1 | |
1794 | goto 111 | |
1795 | endif | |
1796 | ||
1797 | om1xrk=x | |
1798 | ||
1799 | return | |
1800 | end | |
1801 | ||
1802 | c---------------------------------------------------------------------- | |
1803 | double precision function om1yrk(xh,k) !---test--- | |
1804 | c---------------------------------------------------------------------- | |
1805 | c Random number generated from the function om1yk(xh). We solve the | |
1806 | c equation which give om1yrk by Newton-Raphson + secant method. | |
1807 | c xh - fraction of the energy squared s for the pomeron; | |
1808 | c k - pair indice; | |
1809 | c---------------------------------------------------------------------- | |
1810 | include 'epos.inc' | |
1811 | include 'epos.incems' | |
1812 | ||
1813 | double precision xh,r!,y0,y1,y,gamomy(ntymi:ntymx),eps,ymin,prec,yt | |
1814 | ||
1815 | r=dble(rangen()) | |
1816 | ||
1817 | om1yrk=(0.5d0-r)*log(xh) | |
1818 | return | |
1819 | ||
1820 | end | |
1821 | ||
1822 | c---------------------------------------------------------------------- | |
1823 | function ffom12aii(iq,je1,je2) !---test--- | |
1824 | c---------------------------------------------------------------------- | |
1825 | include 'epos.inc' | |
1826 | ig=5 | |
1827 | xmin=0.01/engy | |
1828 | xmax=1 | |
1829 | r2=0 | |
1830 | do i2=1,ig | |
1831 | do m2=1,2 | |
1832 | xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5)) | |
1833 | r1=0 | |
1834 | do i1=1,ig | |
1835 | do m1=1,2 | |
1836 | xp=xmin+(xmax-xmin)*(.5+tgss(ig,i1)*(m1-1.5)) | |
1837 | f=ffom12a(xp,xm,iq,iq,je1,je2) | |
1838 | r1=r1+wgss(ig,i1)*f | |
1839 | enddo | |
1840 | enddo | |
1841 | f=r1*0.5*(xmax-xmin) | |
1842 | r2=r2+wgss(ig,i2)*f | |
1843 | enddo | |
1844 | enddo | |
1845 | ffom12aii=r2*0.5*(xmax-xmin) | |
1846 | end | |
1847 | ||
1848 | c---------------------------------------------------------------------- | |
1849 | function ffom12ai(xp,iq1,iq2,je1,je2) !---test--- | |
1850 | c---------------------------------------------------------------------- | |
1851 | include 'epos.inc' | |
1852 | ig=5 | |
1853 | xmin=0.01/engy | |
1854 | xmax=1 | |
1855 | r2=0 | |
1856 | do i2=1,ig | |
1857 | do m2=1,2 | |
1858 | xm=xmin+(xmax-xmin)*(.5+tgss(ig,i2)*(m2-1.5)) | |
1859 | f=ffom12a(xp,xm,iq1,iq2,je1,je2) | |
1860 | r2=r2+wgss(ig,i2)*f | |
1861 | enddo | |
1862 | enddo | |
1863 | ffom12ai=r2*0.5*(xmax-xmin) | |
1864 | end | |
1865 | ||
1866 | c---------------------------------------------------------------------- | |
1867 | function ffom12a(xp,xm,iq1,iq2,je1,je2) !---test--- | |
1868 | c---------------------------------------------------------------------- | |
1869 | c | |
1870 | c 2*om52*F*F == PomInc | |
1871 | c | |
1872 | c xp - xplus | |
1873 | c xm - xminus | |
1874 | c iq=1 .... sea-sea | |
1875 | c iq1 - min iq iq=2 .... val-sea | |
1876 | c iq2 - max iq iq=3 .... sea-val | |
1877 | c iq=4 .... val-val | |
1878 | c je = emission type (projectile and target side) | |
1879 | c 0 ... no emissions | |
1880 | c 1 ... emissions | |
1881 | c else ... all | |
1882 | c | |
1883 | c already b-averaged (\int d2b /sigine*10) | |
1884 | c---------------------------------------------------------------------- | |
1885 | include 'epos.inc' | |
1886 | ||
1887 | sy=engy*engy | |
1888 | xh=xm*xp | |
1889 | ctp060829 yp=0.5*log(xp/xm) | |
1890 | ffom12a=0 | |
1891 | do i=iq1,iq2 | |
1892 | if(i.eq.1)then | |
1893 | ffom12a=ffom12a+2*om52pi(sy*xh,1.,1.,0,je1,je2) | |
1894 | elseif(i.eq.2)then | |
1895 | ffom12a=ffom12a+2*om52pi(sy*xh,xp,1.,1,je1,je2) | |
1896 | elseif(i.eq.3)then | |
1897 | ffom12a=ffom12a+2*om52pi(sy*xh,xm,1.,2,je1,je2) | |
1898 | elseif(i.eq.4)then | |
1899 | ffom12a=ffom12a+2*om52pi(sy*xh,xp,xm,3,je1,je2) | |
1900 | endif | |
1901 | enddo | |
1902 | ffom12a=ffom12a | |
1903 | * *alpff(iclpro)*xp**betff(1)*(1-xp)**alplea(iclpro) | |
1904 | * *alpff(icltar)*xm**betff(2)*(1-xm)**alplea(icltar) | |
1905 | ||
1906 | end | |
1907 | ||
1908 | c---------------------------------------------------------------------- | |
1909 | function ffom11a(xp,xm,iq1,iq2) !---test--- | |
1910 | c---------------------------------------------------------------------- | |
1911 | c | |
1912 | c int(db) om1ff /sigine*10 | |
1913 | c | |
1914 | c xp - xplus iq=-1 ... fit | |
1915 | c xm - xminus iq=0 .... soft | |
1916 | c iq=1 .... gg | |
1917 | c iq1 - min iq iq=2 .... qg | |
1918 | c iq2 - max iq iq=3 .... gq | |
1919 | c iq=4 .... qq | |
1920 | c iq=5 .... diff | |
1921 | c---------------------------------------------------------------------- | |
1922 | include 'epos.inc' | |
1923 | common/geom/rmproj,rmtarg,bmax,bkmx | |
1924 | ig=5 | |
1925 | bmid=bkmx/2. | |
1926 | r=0.d0 | |
1927 | do i=1,ig | |
1928 | do m=1,2 | |
1929 | bb=bmid*(1.+(2.*m-3)*tgss(ig,i)) | |
1930 | f=ffom11(xp,xm,bb,iq1,iq2) | |
1931 | r=r+bb*wgss(ig,i)*f | |
1932 | enddo | |
1933 | enddo | |
1934 | ffom11a=r*2.*pi*bmid /sigine*10 | |
1935 | return | |
1936 | end | |
1937 | ||
1938 | c---------------------------------------------------------------------- | |
1939 | function ffom11(xp,xm,b,iq1,iq2) !---test--- | |
1940 | c---------------------------------------------------------------------- | |
1941 | c | |
1942 | c 2*om5*F*F == PomInc | |
1943 | c | |
1944 | c xp - xplus iq=-1 ... fit | |
1945 | c xm - xminus iq=0 .... soft | |
1946 | c b - impact parameter iq=1 .... gg | |
1947 | c iq1 - min iq iq=2 .... qg | |
1948 | c iq2 - max iq iq=3 .... gq | |
1949 | c iq=4 .... qq | |
1950 | c iq=5 .... diff | |
1951 | c---------------------------------------------------------------------- | |
1952 | include 'epos.inc' | |
1953 | double precision om51 | |
1954 | ||
1955 | if(xm.ge.0.)then | |
1956 | ||
1957 | xh=xm*xp | |
1958 | yp=0.5*log(xp/xm) | |
1959 | ffom11=2*om51(dble(xh),dble(yp),b,iq1,iq2) | |
1960 | * *(1-xm)**alplea(icltar)*(1-xp)**alplea(iclpro) | |
1961 | ||
1962 | else !xm integration | |
1963 | ||
1964 | ig=5 | |
1965 | xmin=0.01/engy | |
1966 | xmax=1 | |
1967 | r=0 | |
1968 | do i=1,ig | |
1969 | do m=1,2 | |
1970 | xmm=xmin*(xmax/xmin)**(.5+tgss(ig,i)*(m-1.5)) | |
1971 | xh=xmm*xp | |
1972 | yp=0.5*log(xp/xmm) | |
1973 | f=2*om51(dble(xh),dble(yp),b,iq1,iq2) | |
1974 | * *(1-xmm)**alplea(icltar)*(1-xp)**alplea(iclpro) | |
1975 | r=r+wgss(ig,i)*f*xmm | |
1976 | enddo | |
1977 | enddo | |
1978 | ffom11=r*0.5*log(xmax/xmin) | |
1979 | ||
1980 | endif | |
1981 | ||
1982 | end | |
1983 | ||
1984 | c---------------------------------------------------------------------- | |
1985 | double precision function om51(xh,yp,b,iq1,iq2) !---test--- | |
1986 | c---------------------------------------------------------------------- | |
1987 | c xh - xplus*xminus iq=-1 ... fit (om1 * 0.5) | |
1988 | c yp - rapidity iq=0 .... soft | |
1989 | c b - impact param iq=1 .... gg | |
1990 | c iq1 - min iq iq=2 .... qg | |
1991 | c iq2 - max iq iq=3 .... gq | |
1992 | c iq=4 .... qq | |
1993 | c iq=5 .... diff | |
1994 | c---------------------------------------------------------------------- | |
1995 | include 'epos.inc' | |
1996 | include 'epos.incsem' | |
1997 | include 'epos.incpar' | |
1998 | double precision xp,xm,xh,yp,om51p,om1 | |
1999 | om51=0.d0 | |
2000 | if(xh.le.0.d0.or.xh.gt.1.d0)return | |
2001 | ||
2002 | sy=engy*engy | |
2003 | xp=sqrt(xh)*exp(yp) | |
2004 | xm=xh/xp | |
2005 | ||
2006 | if(iq1.eq.-1.and.iq2.eq.-1)then | |
2007 | om51=0.5d0*om1(xh,yp,b) | |
2008 | elseif(iq1.ge.0)then | |
2009 | om51=0.d0 | |
2010 | do i=iq1,iq2 | |
2011 | if(i.ne.5)then | |
2012 | i1=min(i,1) | |
2013 | call GfunPar(1,i1,b,sy,alp,bet,betp,epsp,epst,epss,gamv) | |
2014 | om51=om51+om51p(sy*real(xh),xh,yp,b,i) | |
2015 | * *xp**dble(epsp)*xm**dble(epst) | |
2016 | * *dble(sy)**dble(epss) | |
2017 | else | |
2018 | call GfunPar(1,2,b,sy,alp,bet,betp,epsp,epst,epss,gamv) | |
2019 | om51=om51+0.5d0*alp*real(xp)**bet*real(xm)**betp | |
2020 | endif | |
2021 | enddo | |
2022 | else | |
2023 | stop'om5: choice of iq1 and iq2 is nonsense. ' | |
2024 | endif | |
2025 | ||
2026 | end | |
2027 | ||
2028 | c---------------------------------------------------------------------- | |
2029 | double precision function om5s(sx,xh,yp,b,iq1,iq2) !---test--- | |
2030 | c---------------------------------------------------------------------- | |
2031 | include 'epos.inc' | |
2032 | double precision om51 | |
2033 | double precision xh,yp | |
2034 | ss=sx/xh | |
2035 | engyx=engy | |
2036 | engy=sqrt(ss) | |
2037 | om5s=om51(xh,yp,b,iq1,iq2) | |
2038 | engy=engyx | |
2039 | end | |
2040 | ||
2041 | c---------------------------------------------------------------------- | |
2042 | double precision function om5Jk(k,xh,yp,iqq) !---MC--- | |
2043 | c---------------------------------------------------------------------- | |
2044 | c partial om5 | |
2045 | c xh - fraction of the energy squared s for the pomeron; | |
2046 | c b - impact parameter between the pomeron ends; | |
2047 | c yp - rapidity for the pomeron; | |
2048 | c iqq=0 - soft | |
2049 | c iqq=1 - gg | |
2050 | c iqq=2 - qg | |
2051 | c iqq=3 - gq | |
2052 | c iqq=4 - qq | |
2053 | c iqq=5 - diffractif | |
2054 | c---------------------------------------------------------------------- | |
2055 | include 'epos.inc' | |
2056 | include 'epos.incems' | |
2057 | include 'epos.incsem' | |
2058 | ||
2059 | double precision xh,yp,om51p | |
2060 | double precision plc,s | |
2061 | common/cems5/plc,s | |
2062 | ||
2063 | sy=real(s*xh) | |
2064 | b=bk(k) | |
2065 | ||
2066 | if(iqq.ne.5)then | |
2067 | om5Jk=om51p(sy,xh,yp,b,iqq) | |
2068 | if(iscreen.ne.0)then | |
2069 | xp=sqrt(xh)*exp(yp) | |
2070 | xm=xh/xp | |
2071 | if(iqq.eq.0)then | |
2072 | epsG=epsilongs(k,0) | |
2073 | epsGp=epsilongp(k,0) | |
2074 | epsGt=epsilongt(k,0) | |
2075 | else | |
2076 | epsG=epsilongs(k,1) | |
2077 | epsGp=epsilongp(k,1) | |
2078 | epsGt=epsilongt(k,1) | |
2079 | endif | |
2080 | epsG=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsG | |
2081 | epsGp=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsGp | |
2082 | epsGt=exp(-gfactor*(zparpro(k)*zpartar(k)))*epsGt | |
2083 | om5Jk=om5Jk*xp**dble(epsGp)*xm**dble(epsGt)*s**dble(epsG) | |
2084 | endif | |
2085 | else | |
2086 | xp=sqrt(xh)*exp(yp) | |
2087 | xm=xh/xp | |
2088 | om5Jk=0.5d0*atilde(2,k)*xp**btildep(2,k)*xm**btildepp(2,k) | |
2089 | endif | |
2090 | ||
2091 | return | |
2092 | end | |
2093 | ||
2094 | c---------------------------------------------------------------------- | |
2095 | double precision function omIgamint(b,iqq) !---test--- | |
2096 | c---------------------------------------------------------------------- | |
2097 | c - integrated chi~(b)FF/2 for cut I diagram (simple Pomeron) | |
2098 | c b - impact parameter between the pomeron ends; | |
2099 | c yp - rapidity for the pomeron; | |
2100 | c iqq=0 effective one | |
2101 | c iqq=1 soft | |
2102 | c iqq=2 gg | |
2103 | c---------------------------------------------------------------------- | |
2104 | include 'epos.inc' | |
2105 | include 'epos.incems' | |
2106 | include 'epos.incsem' | |
2107 | include 'epos.incpar' | |
2108 | ||
2109 | double precision Df | |
2110 | ||
2111 | Df=0.d0 | |
2112 | sy=engy*engy | |
2113 | omIgamint=0.d0 | |
2114 | imax=idxD1 | |
2115 | if(iomega.eq.2)imax=1 | |
2116 | ||
2117 | if(iqq.eq.0)then | |
2118 | coefp=1.+alplea(iclpro) | |
2119 | coeft=1.+alplea(icltar) | |
2120 | ||
2121 | do i=idxDmin,imax | |
2122 | call GfunPar(1,i,b,sy,alpx,betx,betpx,epsp,epst,epss,gamv) | |
2123 | betp=1.+betx | |
2124 | betpp=1.+betpx | |
2125 | Df=alpx*dble(utgam1(betp)*utgam1(betpp)*ucfpro | |
2126 | * *ucftar/utgam1(betp+coefp)/utgam1(betpp+coeft)) | |
2127 | omIgamint=omIgamint+Df | |
2128 | enddo | |
2129 | else | |
2130 | call utstop('Wrong iqq in omIgamint&') | |
2131 | endif | |
2132 | ||
2133 | omIgamint=omIgamint | |
2134 | * *dble(chad(iclpro)*chad(icltar)) | |
2135 | ||
2136 | omIgamint=omIgamint*0.5d0 | |
2137 | ||
2138 | return | |
2139 | end | |
2140 | ||
2141 | c----------------------------------------------------------------------- | |
2142 | subroutine WomTy(w,xh,yp,k) | |
2143 | c----------------------------------------------------------------------- | |
2144 | c - w(ity) for group iqq of cut enhanced diagram giving | |
2145 | c the probability of the type of the same final state. | |
2146 | c k - pair indice; | |
2147 | c xh - fraction of the energy squared s for the pomeron; | |
2148 | c yp - rapidity for the pomeron; | |
2149 | c xpr,xmr impulsion fraction of remnant | |
2150 | c ity = 0 - soft | |
2151 | c ity = 1 - gg | |
2152 | c ity = 2 - qg | |
2153 | c ity = 3 - gq | |
2154 | c ity = 4 - qq | |
2155 | c----------------------------------------------------------------------- | |
2156 | include 'epos.inc' | |
2157 | include 'epos.incems' | |
2158 | include 'epos.incsem' | |
2159 | doubleprecision xh,yp,om5Jk,w(0:7) | |
2160 | ||
2161 | do i=0,7 | |
2162 | w(i)=0.d0 | |
2163 | enddo | |
2164 | ||
2165 | do i=0,5 | |
2166 | w(i)=om5Jk(k,xh,yp,i) | |
2167 | enddo | |
2168 | ||
2169 | if(iotst1.gt.0)then !???????????????????????????????????????????????? | |
2170 | corfac=float(iotst1) | |
2171 | ww=w(0)+w(1)+w(2)+w(3)+w(4)+w(5) | |
2172 | whard=0. | |
2173 | do i=1,4 | |
2174 | w(i)=w(i)*corfac | |
2175 | whard=whard+w(i) | |
2176 | enddo | |
2177 | if(whard.gt.ww)then | |
2178 | do i=1,4 | |
2179 | w(i)=w(i)/whard*ww | |
2180 | enddo | |
2181 | whard=ww | |
2182 | endif | |
2183 | w05=w(0)+w(5) | |
2184 | if(whard.lt.ww)then | |
2185 | w(0)=w(0)/w05*(ww-whard) | |
2186 | w(5)=w(5)/w05*(ww-whard) | |
2187 | else | |
2188 | w(0)=0 | |
2189 | w(5)=0 | |
2190 | endif | |
2191 | !write(*,'(2f11.4)')(ww-w05)/ww,(ww-w(0)-w(5))/ww | |
2192 | endif !????????????????????????????????????????????????????????????????? | |
2193 | ||
2194 | return | |
2195 | end | |
2196 | ||
2197 | ||
2198 | c----------------------------------------------------------------------- | |
2199 | double precision function Womegak(xp,xm,xprem,xmrem,k) !---MC--- | |
2200 | c----------------------------------------------------------------------- | |
2201 | c - sum(omGam(xp,xm))*(1-xp)*(1-xm) for group of cut enhanced | |
2202 | c diagram giving the same final state (without nuclear effect). | |
2203 | c xp,xm - fraction of the loght cone momenta of the pomeron; | |
2204 | c k - pair index | |
2205 | c----------------------------------------------------------------------- | |
2206 | include 'epos.inc' | |
2207 | include 'epos.incems' | |
2208 | double precision xp,xm,xprem,xmrem | |
2209 | ||
2210 | Womegak=0.d0 | |
2211 | ||
2212 | imax=ntymx | |
2213 | if(iomega.eq.2)imax=1 | |
2214 | ||
2215 | do i=ntymin,imax | |
2216 | Womegak=Womegak+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k) | |
2217 | enddo | |
2218 | ||
2219 | ||
2220 | Womegak=Womegak*(xprem-xp)*(xmrem-xm) | |
2221 | ||
2222 | return | |
2223 | end | |
2224 | ||
2225 | ||
2226 | cc---------------------------------------------------------------------- | |
2227 | c double precision function omNpcut(xp,xm,xprem,xmrem,bh,iqq) !---test--- | |
2228 | cc---------------------------------------------------------------------- | |
2229 | cc Sum of all cut diagrams | |
2230 | cc iqq=0 ideal G | |
2231 | cc iqq=1 exact G + diff | |
2232 | cc---------------------------------------------------------------------- | |
2233 | c include "epos.inc" | |
2234 | c double precision om51,xh,yp,xprem,xmrem,xp,xm!,omYcutI | |
2235 | c | |
2236 | c omNpcut=0.d0 | |
2237 | c xh=xp*xm | |
2238 | c if(abs(xh).gt.1.d-10)then | |
2239 | c yp=0.5d0*log(xp/xm) | |
2240 | c else | |
2241 | c yp=0.d0 | |
2242 | c endif | |
2243 | c | |
2244 | c if(iqq.eq.0)omNpcut=om51(xh,yp,bh,-1,-1) | |
2245 | c if(iqq.eq.1)omNpcut=om51(xh,yp,bh,0,5) | |
2246 | c | |
2247 | c omNpcut=omNpcut*2.d0 | |
2248 | c | |
2249 | c return | |
2250 | c end | |
2251 | c | |
2252 | c---------------------------------------------------------------------- | |
2253 | double precision function omGam(xp,xm,bh) !---test--- | |
2254 | c----------------------------------------------------------------------- | |
2255 | c Cut diagram part for calculation of probability distribution | |
2256 | c xp,xm impulsion fraction of remnant | |
2257 | c bh - impact parameter between the pomeron ends; | |
2258 | c----------------------------------------------------------------------- | |
2259 | include "epos.inc" | |
2260 | include "epos.incems" | |
2261 | double precision om51,xp,xm,xh,yp,eps!,omYgam | |
2262 | parameter (eps=1.d-20) | |
2263 | ||
2264 | omGam=0.d0 | |
2265 | if(xp.lt.eps.or.xm.lt.eps)return | |
2266 | xh=xp*xm | |
2267 | if(abs(xh).gt.1.d-10)then | |
2268 | yp=0.5d0*log(xp/xm) | |
2269 | else | |
2270 | yp=0.d0 | |
2271 | endif | |
2272 | ||
2273 | omGam=om51(xh,yp,bh,-1,-1) | |
2274 | ||
2275 | omGam=2.d0*omGam | |
2276 | ||
2277 | return | |
2278 | end | |
2279 | ||
2280 | c---------------------------------------------------------------------- | |
2281 | double precision function omGamk(k,xp,xm) !---MC--- | |
2282 | c----------------------------------------------------------------------- | |
2283 | c Cut diagram part for calculation of probability distribution (for omega) | |
2284 | c xp,xm impulsion fraction of remnant | |
2285 | c bh - impact parameter between the pomeron ends; | |
2286 | c----------------------------------------------------------------------- | |
2287 | include "epos.inc" | |
2288 | include "epos.incems" | |
2289 | double precision xp,xm | |
2290 | omGamk=0.d0 | |
2291 | imax=ntymx | |
2292 | if(iomega.eq.2)imax=1 | |
2293 | do i=ntymin,imax | |
2294 | omGamk=omGamk+atilde(i,k)*xp**btildep(i,k)*xm**btildepp(i,k) | |
2295 | enddo | |
2296 | ||
2297 | return | |
2298 | end | |
2299 | ||
2300 | c---------------------------------------------------------------------- | |
2301 | double precision function omGamint(bh) !---test--- | |
2302 | c----------------------------------------------------------------------- | |
2303 | c Integrated cut diagram part for calculation of probability distribution | |
2304 | c bh - impact parameter between the pomeron ends; | |
2305 | c----------------------------------------------------------------------- | |
2306 | include "epos.inc" | |
2307 | double precision omIgamint!,omYgamint | |
2308 | ||
2309 | omGamint=2.d0*omIgamint(bh,0) | |
2310 | ||
2311 | return | |
2312 | end | |
2313 | ||
2314 | ||
2315 | ||
2316 | ||
2317 | ||
2318 | c---------------------------------------------------------------------- | |
2319 | block data dgdata | |
2320 | c---------------------------------------------------------------------- | |
2321 | c constants for numerical integration (gaussian weights) | |
2322 | c---------------------------------------------------------------------- | |
2323 | double precision dgx1,dga1 | |
2324 | common /dga20/ dgx1(10),dga1(10) | |
2325 | ||
2326 | ||
2327 | data dgx1/ | |
2328 | & .765265211334973D-01, | |
2329 | & .227785851141645D+00, | |
2330 | & .373706088715420D+00, | |
2331 | & .510867001950827D+00, | |
2332 | & .636053680726515D+00, | |
2333 | & .746331906460151D+00, | |
2334 | & .839116971822219D+00, | |
2335 | & .912234428251326D+00, | |
2336 | & .963971927277914D+00, | |
2337 | & .993128599185095D+00/ | |
2338 | data dga1/ | |
2339 | & .152753387130726D+00, | |
2340 | & .149172986472604D+00, | |
2341 | & .142096109318382D+00, | |
2342 | & .131688638449177D+00, | |
2343 | & .118194531961518D+00, | |
2344 | & .101930119817233D+00, | |
2345 | & .832767415767047D-01, | |
2346 | & .626720483341090D-01, | |
2347 | & .406014298003871D-01, | |
2348 | & .176140071391506D-01/ | |
2349 | ||
2350 | end | |
2351 | ||
2352 | ||
2353 | ||
2354 | c---------------------------------------------------------------------- | |
2355 | double precision function PhiExact(fj,xp,xm,s,b) !---test--- | |
2356 | c---------------------------------------------------------------------- | |
2357 | c Exact expression of the Phi function for pp collision | |
2358 | c---------------------------------------------------------------------- | |
2359 | include 'epos.inc' | |
2360 | include 'epos.incsem' | |
2361 | include 'epos.incems' | |
2362 | double precision al(idxD0:idxD1),betp(idxD0:idxD1) | |
2363 | *,z,xIrst!,ffacto | |
2364 | double precision zp(idxD0:idxD1),Phitmp,betpp(idxD0:idxD1) | |
2365 | *,yp,ym,xm,xp | |
2366 | double precision eps | |
2367 | parameter(eps=1.d-20) | |
2368 | dimension ipr(idxD0:idxD1),imax(idxD0:idxD1) | |
2369 | ||
2370 | if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in PhiExact" | |
2371 | PhiExact=0.d0 | |
2372 | ||
2373 | if(xp.gt.eps.and.xm.gt.eps.and.xp.le.1.d0+eps | |
2374 | & .and.xm.le.1.d0+eps)then | |
2375 | ||
2376 | ||
2377 | do i=idxD0,idxD1 | |
2378 | imax(i)=0 | |
2379 | ipr(i)=0 | |
2380 | zp(i)=1.d0 | |
2381 | al(i)=0.d0 | |
2382 | betp(i)=0.d0 | |
2383 | betpp(i)=0.d0 | |
2384 | enddo | |
2385 | ||
2386 | imax0=idxD1 | |
2387 | if(iomega.eq.2)imax0=1 | |
2388 | ||
2389 | do i=idxDmin,imax0 | |
2390 | imax(i)=10+max(5,int(log10(s))) | |
2391 | if(b.ge.1.)imax(i)=4+max(3,int(log10(sqrt(s)))) | |
2392 | imax(i)=min(30,imax(i)) | |
2393 | enddo | |
2394 | Phitmp=0.d0 | |
2395 | do i=idxDmin,imax0 | |
2396 | call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
2397 | betp(i)=dble(betx)+1.d0 | |
2398 | betpp(i)=dble(betpx)+1.d0 | |
2399 | al(i)=dble(alpx*gamv) | |
2400 | * *dble(chad(iclpro)*chad(icltar)) | |
2401 | enddo | |
2402 | ||
2403 | do ipr0=0,imax(0) | |
2404 | ipr(0)=ipr0 | |
2405 | zp(0)=1.d0 | |
2406 | if (ipr(0).ne.0) zp(0)=(-dble(fj)*al(0))**ipr(0)*facto(ipr(0)) | |
2407 | do ipr1=0,imax(1) | |
2408 | ipr(1)=ipr1 | |
2409 | zp(1)=1.d0 | |
2410 | if (ipr(1).ne.0) zp(1)=(-dble(fj)*al(1))**ipr(1)*facto(ipr(1)) | |
2411 | do ipr2=0,imax(2) | |
2412 | ipr(2)=ipr2 | |
2413 | zp(2)=1.d0 | |
2414 | if (ipr(2).ne.0) zp(2)=(-dble(fj)*al(2))**ipr(2)*facto(ipr(2)) | |
2415 | yp=0.d0 | |
2416 | ym=0.d0 | |
2417 | z=1.d0 | |
2418 | isum=0 | |
2419 | do i=idxDmin,imax0 | |
2420 | yp=yp+dble(ipr(i))*betp(i) | |
2421 | ym=ym+dble(ipr(i))*betpp(i) | |
2422 | isum=isum+ipr(i) | |
2423 | z=z*zp(i) | |
2424 | enddo | |
2425 | ||
2426 | z=z*xIrst(1,xp,yp,betp,ipr) | |
2427 | z=z*xIrst(2,xm,ym,betpp,ipr) | |
2428 | ||
2429 | Phitmp=Phitmp+z | |
2430 | ||
2431 | enddo | |
2432 | enddo | |
2433 | enddo | |
2434 | ||
2435 | ||
2436 | ||
2437 | endif | |
2438 | ||
2439 | PhiExact=Phitmp | |
2440 | ||
2441 | ||
2442 | return | |
2443 | end | |
2444 | ||
2445 | ||
2446 | c---------------------------------------------------------------------- | |
2447 | double precision function PhiExpoK(k,xp,xm) !---MC--- | |
2448 | c---------------------------------------------------------------------- | |
2449 | c Exponential expression of the Phi function for pp collision | |
2450 | c for given k | |
2451 | c---------------------------------------------------------------------- | |
2452 | include 'epos.inc' | |
2453 | include 'epos.incsem' | |
2454 | include 'epos.incems' | |
2455 | ||
2456 | double precision xp,xm,Phitmp,Gt1 | |
2457 | double precision atildg,btildgp,btildgpp | |
2458 | common/cgtilde/atildg(idxD0:idxD1,kollmx) | |
2459 | *,btildgp(idxD0:idxD1,kollmx),btildgpp(idxD0:idxD1,kollmx) | |
2460 | ||
2461 | ||
2462 | Phitmp=0.d0 | |
2463 | ||
2464 | imax=idxD1 | |
2465 | if(iomega.eq.2)imax=1 | |
2466 | ||
2467 | Phitmp=0.d0 | |
2468 | Gt1=0.d0 | |
2469 | do i=idxDmin,imax | |
2470 | Gt1=Gt1+atildg(i,k)*xp**btildgp(i,k)*xm**btildgpp(i,k) | |
2471 | enddo | |
2472 | ||
2473 | Phitmp=exp(-Gt1) | |
2474 | ||
2475 | PhiExpoK=Phitmp | |
2476 | ||
2477 | return | |
2478 | end | |
2479 | ||
2480 | c---------------------------------------------------------------------- | |
2481 | double precision function PhiExpo(fj,xp,xm,s,b) !---MC--- | |
2482 | c---------------------------------------------------------------------- | |
2483 | c Exponential expression of the Phi function for pp collision | |
2484 | c for given b | |
2485 | c---------------------------------------------------------------------- | |
2486 | include 'epos.inc' | |
2487 | include 'epos.incsem' | |
2488 | include 'epos.incems' | |
2489 | include 'epos.incpar' | |
2490 | ||
2491 | parameter(nbkbin=40) | |
2492 | common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin | |
2493 | double precision AlTi | |
2494 | double precision BeTip,BeTipp | |
2495 | double precision xp,xm,Phitmp,Gt1 | |
2496 | ||
2497 | imax=idxD1 | |
2498 | if(iomega.eq.2)imax=1 | |
2499 | ||
2500 | Gt1=0.d0 | |
2501 | do i=idxDmin,imax | |
2502 | call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
2503 | BeTip =dble(betx) | |
2504 | BeTipp=dble(betpx) | |
2505 | AlTi =dble(alpx) | |
2506 | Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp*dble(fj*xkappa**(2.5*fj-2.5)) | |
2507 | enddo | |
2508 | ||
2509 | Phitmp=exp(-Gt1) | |
2510 | ||
2511 | PhiExpo=Phitmp | |
2512 | & *xp**dble(alplea(iclpro)) | |
2513 | & *xm**dble(alplea(icltar)) | |
2514 | ||
2515 | return | |
2516 | end | |
2517 | ||
2518 | c---------------------------------------------------------------------- | |
2519 | double precision function PhiUnit(xp,xm) !---test--- | |
2520 | c---------------------------------------------------------------------- | |
2521 | c Exponential expression of the Phi function for pp collision | |
2522 | c for given b | |
2523 | c---------------------------------------------------------------------- | |
2524 | include 'epos.inc' | |
2525 | include 'epos.incsem' | |
2526 | include 'epos.incems' | |
2527 | include 'epos.incpar' | |
2528 | ||
2529 | double precision AlTi | |
2530 | double precision BeTip,BeTipp | |
2531 | double precision xp,xm,Phitmp,Gt1 | |
2532 | ||
2533 | imax=idxD1 | |
2534 | if(iomega.eq.2)imax=1 | |
2535 | ||
2536 | Gt1=0.d0 | |
2537 | do i=idxDmin,imax | |
2538 | BeTip =betUni(i,2) | |
2539 | BeTipp=betpUni(i,2) | |
2540 | AlTi =alpUni(i,2) | |
2541 | Gt1=Gt1+AlTi*xp**BeTip*xm**BeTipp | |
2542 | c write(ifch,*)'Phiunit',i,xp,xm,Gt1,AlTi,BeTip,BeTipp | |
2543 | enddo | |
2544 | ||
2545 | Phitmp=exp(-Gt1) | |
2546 | ||
2547 | PhiUnit=Phitmp | |
2548 | & *xp**dble(alplea(iclpro)) | |
2549 | & *xm**dble(alplea(icltar)) | |
2550 | ||
2551 | ||
2552 | return | |
2553 | end | |
2554 | ||
2555 | ||
2556 | cc---------------------------------------------------------------------- | |
2557 | c double precision function PhiUnit(xp,xm,s,b) !---inu--- | |
2558 | cc---------------------------------------------------------------------- | |
2559 | c include 'epos.inc' | |
2560 | c double precision xp,xm,PhiExpo,Znorm | |
2561 | c | |
2562 | c PhiUnit=PhiExpo(1.,xp,xm,s,b) | |
2563 | c & /Znorm(s,b) | |
2564 | c | |
2565 | c return | |
2566 | c end | |
2567 | c | |
2568 | ||
2569 | c---------------------------------------------------------------------- | |
2570 | double precision function Hrst(s,b,xp,xm) !test | |
2571 | c---------------------------------------------------------------------- | |
2572 | include 'epos.inc' | |
2573 | include 'epos.incems' | |
2574 | include 'epos.incsem' | |
2575 | include 'epos.incpar' | |
2576 | parameter(idxD2=8) | |
2577 | double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni | |
2578 | common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2), | |
2579 | & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2), | |
2580 | & HalpUni(idxD0:idxD2) | |
2581 | double precision al(idxD0:idxD2),betp(idxD0:idxD2) | |
2582 | *,z,xJrst!,ffacto | |
2583 | double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2) | |
2584 | *,yp,ym,xp,xm | |
2585 | dimension ipr(idxD0:idxD2),imax(idxD0:idxD2) | |
2586 | ||
2587 | if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in Hrst" | |
2588 | Hrst=0.d0 | |
2589 | do i=idxD0,idxD2 | |
2590 | imax(i)=0 | |
2591 | ipr(i)=0 | |
2592 | zp(i)=1.d0 | |
2593 | al(i)=0.d0 | |
2594 | enddo | |
2595 | ||
2596 | if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then | |
2597 | ||
2598 | imax0=idxD1 | |
2599 | if(iomega.eq.2)imax0=1 | |
2600 | imax1=idxD2 | |
2601 | if(iomega.eq.2)imax1=imax1-1 | |
2602 | ||
2603 | do i=idxDmin,imax1 | |
2604 | imax(i)=max(2,int(log10(100.*s)/3.)) | |
2605 | c if(i.ge.2)imax(i)=imax(i)*2 | |
2606 | if(b.ge.1.5)imax(i)=2 !max(2,imax(i)/2) | |
2607 | imax(i)=min(30,imax(i)) | |
2608 | if(i.gt.imax0)then | |
2609 | if((zpUni*ztUni.lt.1.d-6) | |
2610 | & .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then | |
2611 | imax(i)=0 | |
2612 | else | |
2613 | imax(i)=1 !imax(i)/3 | |
2614 | endif | |
2615 | endif | |
2616 | enddo | |
2617 | ||
2618 | Htmp=0.d0 | |
2619 | do i=idxDmin,imax1 | |
2620 | betp(i)=HbetUni(i) | |
2621 | betpp(i)=HbetpUni(i) | |
2622 | al(i)=HalpUni(i) | |
2623 | enddo | |
2624 | ||
2625 | do ipr0=0,imax(0) | |
2626 | c write(ifmt,*)'Hrst ipr0,xp,xm :',ipr0,xp,xm | |
2627 | ipr(0)=ipr0 | |
2628 | zp(0)=1.d0 | |
2629 | if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0)) | |
2630 | do ipr1=0,imax(1) | |
2631 | ipr(1)=ipr1 | |
2632 | zp(1)=1.d0 | |
2633 | if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1)) | |
2634 | do ipr2=0,imax(2) | |
2635 | ipr(2)=ipr2 | |
2636 | zp(2)=1.d0 | |
2637 | if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2)) | |
2638 | do ipr3=0,imax(3) | |
2639 | ipr(3)=ipr3 | |
2640 | zp(3)=1.d0 | |
2641 | if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3)) | |
2642 | do ipr4=0,imax(4) | |
2643 | ipr(4)=ipr4 | |
2644 | zp(4)=1.d0 | |
2645 | if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4)) | |
2646 | do ipr5=0,imax(5) | |
2647 | ipr(5)=ipr5 | |
2648 | zp(5)=1.d0 | |
2649 | if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5)) | |
2650 | do ipr6=0,imax(6) | |
2651 | ipr(6)=ipr6 | |
2652 | zp(6)=1.d0 | |
2653 | if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6)) | |
2654 | do ipr7=0,imax(7) | |
2655 | ipr(7)=ipr7 | |
2656 | zp(7)=1.d0 | |
2657 | if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7)) | |
2658 | do ipr8=0,imax(8) | |
2659 | ipr(8)=ipr8 | |
2660 | zp(8)=1.d0 | |
2661 | if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8)) | |
2662 | if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5) | |
2663 | & +ipr(6)+ipr(7)+ipr(8).ne.0) then | |
2664 | yp=0.d0 | |
2665 | ym=0.d0 | |
2666 | z=1.d0 | |
2667 | do i=idxDmin,imax1 | |
2668 | yp=yp+dble(ipr(i))*betp(i) | |
2669 | ym=ym+dble(ipr(i))*betpp(i) | |
2670 | z=z*zp(i) | |
2671 | enddo | |
2672 | z=z*xJrst(xp,yp,GbetUni,ipr) | |
2673 | z=z*xJrst(xm,ym,GbetpUni,ipr) | |
2674 | Htmp=Htmp+z | |
2675 | endif | |
2676 | enddo | |
2677 | enddo | |
2678 | enddo | |
2679 | enddo | |
2680 | enddo | |
2681 | enddo | |
2682 | enddo | |
2683 | enddo | |
2684 | enddo | |
2685 | ||
2686 | endif | |
2687 | ||
2688 | Hrst=Htmp | |
2689 | ||
2690 | return | |
2691 | end | |
2692 | ||
2693 | c---------------------------------------------------------------------- | |
2694 | double precision function HrstI(s,b,xp,xm) !test | |
2695 | c---------------------------------------------------------------------- | |
2696 | include 'epos.inc' | |
2697 | include 'epos.incems' | |
2698 | include 'epos.incsem' | |
2699 | include 'epos.incpar' | |
2700 | parameter(idxD2=8) | |
2701 | double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni | |
2702 | common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2), | |
2703 | & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2), | |
2704 | & HalpUni(idxD0:idxD2) | |
2705 | double precision al(idxD0:idxD2),betp(idxD0:idxD2) | |
2706 | *,z,xJrstI!,ffacto | |
2707 | double precision zp(idxD0:idxD2),Htmp,betpp(idxD0:idxD2) | |
2708 | *,yp,ym,xp,xm | |
2709 | dimension ipr(idxD0:idxD2),imax(idxD0:idxD2) | |
2710 | ||
2711 | if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI" | |
2712 | HrstI=0.d0 | |
2713 | do i=idxD0,idxD2 | |
2714 | imax(i)=0 | |
2715 | ipr(i)=0 | |
2716 | zp(i)=1.d0 | |
2717 | al(i)=0.d0 | |
2718 | enddo | |
2719 | ||
2720 | ||
2721 | if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.le.1.d0)then | |
2722 | ||
2723 | ||
2724 | imax0=idxD1 | |
2725 | if(iomega.eq.2)imax0=1 | |
2726 | imax1=idxD2 | |
2727 | if(iomega.eq.2)imax1=imax1-1 | |
2728 | ||
2729 | do i=idxDmin,imax1 | |
2730 | imax(i)=max(3,int(log10(s)/2.)) | |
2731 | c if(i.ge.2)imax(i)=imax(i)*2 | |
2732 | if(b.ge.1.5)imax(i)=max(2,imax(i)/2) | |
2733 | imax(i)=min(30,imax(i)) | |
2734 | if(i.gt.imax0)then | |
2735 | if((zpUni*ztUni.lt.1.d-6) | |
2736 | & .or.(xp.lt.0.1d0.and.xm.lt.0.1d0))then | |
2737 | imax(i)=0 | |
2738 | else | |
2739 | imax(i)=1 !imax(i)/3 | |
2740 | endif | |
2741 | endif | |
2742 | enddo | |
2743 | ||
2744 | Htmp=0.d0 | |
2745 | do i=idxDmin,imax1 | |
2746 | betp(i)=HbetUni(i) | |
2747 | betpp(i)=HbetpUni(i) | |
2748 | al(i)=HalpUni(i) | |
2749 | enddo | |
2750 | do ipr0=0,imax(0) | |
2751 | ipr(0)=ipr0 | |
2752 | zp(0)=1.d0 | |
2753 | if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0)) | |
2754 | do ipr1=0,imax(1) | |
2755 | ipr(1)=ipr1 | |
2756 | zp(1)=1.d0 | |
2757 | if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1)) | |
2758 | do ipr2=0,imax(2) | |
2759 | ipr(2)=ipr2 | |
2760 | zp(2)=1.d0 | |
2761 | if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2)) | |
2762 | do ipr3=0,imax(3) | |
2763 | ipr(3)=ipr3 | |
2764 | zp(3)=1.d0 | |
2765 | if (ipr(3).ne.0) zp(3)=al(3)**ipr(3)*facto(ipr(3)) | |
2766 | do ipr4=0,imax(4) | |
2767 | ipr(4)=ipr4 | |
2768 | zp(4)=1.d0 | |
2769 | if (ipr(4).ne.0) zp(4)=al(4)**ipr(4)*facto(ipr(4)) | |
2770 | do ipr5=0,imax(5) | |
2771 | ipr(5)=ipr5 | |
2772 | zp(5)=1.d0 | |
2773 | if (ipr(5).ne.0) zp(5)=al(5)**ipr(5)*facto(ipr(5)) | |
2774 | do ipr6=0,imax(6) | |
2775 | ipr(6)=ipr6 | |
2776 | zp(6)=1.d0 | |
2777 | if (ipr(6).ne.0) zp(6)=al(6)**ipr(6)*facto(ipr(6)) | |
2778 | do ipr7=0,imax(7) | |
2779 | ipr(7)=ipr7 | |
2780 | zp(7)=1.d0 | |
2781 | if (ipr(7).ne.0) zp(7)=al(7)**ipr(7)*facto(ipr(7)) | |
2782 | do ipr8=0,imax(8) | |
2783 | ipr(8)=ipr8 | |
2784 | zp(8)=1.d0 | |
2785 | if (ipr(8).ne.0) zp(8)=al(8)**ipr(8)*facto(ipr(8)) | |
2786 | if (ipr(0)+ipr(1)+ipr(2)+ipr(3)+ipr(4)+ipr(5) | |
2787 | & +ipr(6)+ipr(7)+ipr(8).ne.0) then | |
2788 | yp=0.d0 | |
2789 | ym=0.d0 | |
2790 | z=1.d0 | |
2791 | do i=idxDmin,imax1 | |
2792 | yp=yp+dble(ipr(i))*betp(i) | |
2793 | ym=ym+dble(ipr(i))*betpp(i) | |
2794 | z=z*zp(i) | |
2795 | enddo | |
2796 | z=z*xJrstI(xp,yp,GbetUni,ipr) | |
2797 | z=z*xJrstI(xm,ym,GbetpUni,ipr) | |
2798 | Htmp=Htmp+z | |
2799 | endif | |
2800 | enddo | |
2801 | enddo | |
2802 | enddo | |
2803 | enddo | |
2804 | enddo | |
2805 | enddo | |
2806 | enddo | |
2807 | enddo | |
2808 | enddo | |
2809 | ||
2810 | endif | |
2811 | ||
2812 | HrstI=Htmp | |
2813 | ||
2814 | return | |
2815 | end | |
2816 | ||
2817 | ||
2818 | ||
2819 | cc---------------------------------------------------------------------- | |
2820 | c double precision function HrstI(s,xp,xm) !---inu--- | |
2821 | cc---------------------------------------------------------------------- | |
2822 | c include 'epos.inc' | |
2823 | c include 'epos.incems' | |
2824 | c include 'epos.incsem' | |
2825 | c include 'epos.incpar' | |
2826 | c double precision al(idxD0:idxD1),betp(idxD0:idxD1) | |
2827 | c *,z,xJrstI!,ffacto | |
2828 | c double precision zp(idxD0:idxD1),Htmp,betpp(idxD0:idxD1) | |
2829 | c *,yp,ym,xp,xm | |
2830 | c dimension ipr(idxD0:idxD1),imax(idxD0:idxD1) | |
2831 | c | |
2832 | c if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI" | |
2833 | c HrstI=0.d0 | |
2834 | c do i=idxD0,idxD1 | |
2835 | c imax(i)=0 | |
2836 | c ipr(i)=0 | |
2837 | c zp(i)=1.d0 | |
2838 | c al(i)=0.d0 | |
2839 | c enddo | |
2840 | c | |
2841 | c if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.lt.1.d0)then | |
2842 | c | |
2843 | c HrstI=0.d0 | |
2844 | c | |
2845 | c imax0=idxD1 | |
2846 | c if(iomega.eq.2)imax0=1 | |
2847 | c | |
2848 | c do i=idxDmin,imax0 | |
2849 | c imax(i)=max(5,int(log10(s))) | |
2850 | cc if(i.ge.2)imax(i)=imax(i)*2 | |
2851 | c imax(i)=min(30,imax(i)) | |
2852 | c enddo | |
2853 | c Htmp=0.d0 | |
2854 | c do i=idxDmin,imax0 | |
2855 | c betp(i)=betUni(i,1)+1.d0 | |
2856 | c betpp(i)=betpUni(i,1)+1.d0 | |
2857 | c al(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar)) | |
2858 | c enddo | |
2859 | c | |
2860 | c do ipr0=0,imax(0) | |
2861 | c ipr(0)=ipr0 | |
2862 | c zp(0)=1.d0 | |
2863 | c if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0)) | |
2864 | c do ipr1=0,imax(1) | |
2865 | c ipr(1)=ipr1 | |
2866 | c zp(1)=1.d0 | |
2867 | c if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1)) | |
2868 | c do ipr2=0,imax(2) | |
2869 | c ipr(2)=ipr2 | |
2870 | c zp(2)=1.d0 | |
2871 | c if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2)) | |
2872 | c if (ipr(0)+ipr(1)+ipr(2).ne.0) then | |
2873 | c yp=0.d0 | |
2874 | c ym=0.d0 | |
2875 | c z=1.d0 | |
2876 | c do i=idxDmin,imax0 | |
2877 | c yp=yp+dble(ipr(i))*betp(i) | |
2878 | c ym=ym+dble(ipr(i))*betpp(i) | |
2879 | c z=z*zp(i) | |
2880 | c enddo | |
2881 | c z=z*xJrstI(xp,yp,betp,ipr) | |
2882 | c z=z*xJrstI(xm,ym,betpp,ipr) | |
2883 | c Htmp=Htmp+z | |
2884 | c endif | |
2885 | c enddo | |
2886 | c enddo | |
2887 | c enddo | |
2888 | c | |
2889 | c HrstI=Htmp | |
2890 | c | |
2891 | c endif | |
2892 | c | |
2893 | c return | |
2894 | c end | |
2895 | c | |
2896 | ||
2897 | cc---------------------------------------------------------------------- | |
2898 | c double precision function ffacto(n) !---test--- | |
2899 | cc---------------------------------------------------------------------- | |
2900 | c | |
2901 | c ffacto=1.D0 | |
2902 | c do i=1,n | |
2903 | c ffacto=ffacto*dble(i) | |
2904 | c enddo | |
2905 | c return | |
2906 | c end | |
2907 | c | |
2908 | ||
2909 | c---------------------------------------------------------------------- | |
2910 | double precision function xIrst(id,x,y,bet,ipr) !---test--- | |
2911 | c---------------------------------------------------------------------- | |
2912 | include 'epos.inc' | |
2913 | double precision y,gammag,utgam2,x,bet(idxD0:idxD1) | |
2914 | dimension ipr(idxD0:idxD1) | |
2915 | ||
2916 | if(id.eq.1)iclrem=iclpro | |
2917 | if(id.eq.2)iclrem=icltar | |
2918 | imax=idxD1 | |
2919 | if(iomega.eq.2)imax=1 | |
2920 | if(y.le.160.)then | |
2921 | xIrst=gammag(iclrem,y)*x**dble(alplea(iclrem)) | |
2922 | else | |
2923 | xIrst=0 | |
2924 | endif | |
2925 | if(xIrst.gt.0.d0)then | |
2926 | do i=idxDmin,imax | |
2927 | if(ipr(i).ne.0.and.bet(i).gt.1.d-10) | |
2928 | & xIrst=xIrst*utgam2(bet(i))**dble(ipr(i)) | |
2929 | enddo | |
2930 | if (abs(y).gt.1.d-10) xIrst=xIrst*x**y | |
2931 | endif | |
2932 | return | |
2933 | end | |
2934 | ||
2935 | ||
2936 | c---------------------------------------------------------------------- | |
2937 | double precision function xJrst(x,y,Gbeta,ipr) !---inu--- | |
2938 | c---------------------------------------------------------------------- | |
2939 | include 'epos.inc' | |
2940 | parameter(idxD2=8) | |
2941 | double precision y,utgam2,x,Gbeta(idxD0:idxD2),eps,gam | |
2942 | dimension ipr(idxD0:idxD2) | |
2943 | ||
2944 | eps=1.d-10 | |
2945 | ||
2946 | imax=idxD2 | |
2947 | if(iomega.eq.2)imax=imax-1 | |
2948 | ||
2949 | gam=utgam2(y) | |
2950 | ||
2951 | if(gam.lt.1.d99)then | |
2952 | ||
2953 | if ((x-1.d0).gt.eps.or.(y-1.d0).gt.eps) then | |
2954 | xJrst=(1.d0-x)**(y-1.d0)/gam | |
2955 | do i=idxDmin,imax | |
2956 | if (ipr(i).ne.0) xJrst=xJrst*Gbeta(i)**dble(ipr(i)) | |
2957 | enddo | |
2958 | else | |
2959 | c write (*,*) 'Warning in xJrst, infinite value !' | |
2960 | xJrst=(1.d0-x+eps)**(y-1.d0)/gam | |
2961 | do i=idxDmin,imax | |
2962 | if (ipr(i).ne.0) xJrst=xJrst*Gbeta(i)**dble(ipr(i)) | |
2963 | enddo | |
2964 | endif | |
2965 | else | |
2966 | xJrst=0.d0 | |
2967 | endif | |
2968 | ||
2969 | return | |
2970 | end | |
2971 | ||
2972 | ||
2973 | c---------------------------------------------------------------------- | |
2974 | double precision function xJrstI(x,y,Gbeta,ipr) !---inu--- | |
2975 | c---------------------------------------------------------------------- | |
2976 | c Function used for the integration of H*Phi. We do the changement of | |
2977 | c variable (1-x)=z**alpha. The power alpha can be change if necessary. | |
2978 | c---------------------------------------------------------------------- | |
2979 | include 'epos.inc' | |
2980 | parameter(idxD2=8) | |
2981 | double precision y,utgam2,x,Gbeta(idxD0:idxD2),alpha,w,gam | |
2982 | dimension ipr(idxD0:idxD2) | |
2983 | ||
2984 | alpha=4.d0 | |
2985 | w=alpha*(y-1.d0)+alpha-1.d0 | |
2986 | imax=idxD2 | |
2987 | if(iomega.eq.2)imax=imax-1 | |
2988 | ||
2989 | gam=utgam2(y) | |
2990 | ||
2991 | if(gam.lt.1.d99)then | |
2992 | ||
2993 | if (w.ge.0)then | |
2994 | ||
2995 | xJrstI=alpha*x**w/gam | |
2996 | do i=idxDmin,imax | |
2997 | if (ipr(i).ne.0) xJrstI=xJrstI*Gbeta(i)**dble(ipr(i)) | |
2998 | enddo | |
2999 | ||
3000 | else | |
3001 | write(*,*) 'x,y,bet,ipr,w',x,y,Gbeta,ipr,w | |
3002 | stop 'Error in xJrstI in epos-omg, integration not possible' | |
3003 | endif | |
3004 | ||
3005 | else | |
3006 | xJrstI=0.d0 | |
3007 | endif | |
3008 | ||
3009 | ||
3010 | return | |
3011 | end | |
3012 | ||
3013 | c---------------------------------------------------------------------- | |
3014 | double precision function HPhiInt(s,b) !---inu--- | |
3015 | c---------------------------------------------------------------------- | |
3016 | c Set integrated over xp and xm (x and y) H(x,y)*Phi(x,y) for a | |
3017 | c given b by gauss method | |
3018 | c---------------------------------------------------------------------- | |
3019 | include 'epos.inc' | |
3020 | parameter(idxD2=8) | |
3021 | double precision GbetUni,GbetpUni,HbetUni,HbetpUni,HalpUni | |
3022 | common/DGamUni/GbetUni( idxD0:idxD2),HbetUni( idxD0:idxD2), | |
3023 | & GbetpUni(idxD0:idxD2),HbetpUni(idxD0:idxD2), | |
3024 | & HalpUni(idxD0:idxD2) | |
3025 | double precision xhm,x,y,yhm,w,Hrst,utgam2,PhiUnit!,PhiExact | |
3026 | c double precision zp2,zm2,HrstI,eps | |
3027 | c common /ar3/ x1(7),a1(7) | |
3028 | common /ar9/ x9(3),a9(3) | |
3029 | ||
3030 | eps=0d0 !1.d-5 | |
3031 | ||
3032 | imax0=idxD1 | |
3033 | imax1=idxD2 | |
3034 | if(iomega.eq.2)then | |
3035 | imax0=1 | |
3036 | imax1=imax1-1 | |
3037 | endif | |
3038 | do i=idxDmin,imax0 | |
3039 | HbetUni(i)=betUni(i,1)+1.d0 | |
3040 | HbetpUni(i)=betpUni(i,1)+1.d0 | |
3041 | GbetUni(i)=utgam2(HbetUni(i)) | |
3042 | GbetpUni(i)=utgam2(HbetpUni(i)) | |
3043 | HalpUni(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar)) | |
3044 | enddo | |
3045 | do i=0,1 | |
3046 | HbetUni(imax0+1+i)=betUni(i,1)+1.d0+betfom | |
3047 | HbetUni(imax0+3+i)=betUni(i,1)+1.d0 | |
3048 | HbetUni(imax0+5+i)=betUni(i,1)+1.d0+betfom | |
3049 | HbetpUni(imax0+1+i)=betpUni(i,1)+1.d0 | |
3050 | HbetpUni(imax0+3+i)=betpUni(i,1)+1.d0+betfom | |
3051 | HbetpUni(imax0+5+i)=betpUni(i,1)+1.d0+betfom | |
3052 | GbetUni(imax0+1+i)=utgam2(HbetUni(imax0+1+i)) | |
3053 | GbetUni(imax0+3+i)=utgam2(HbetUni(imax0+3+i)) | |
3054 | GbetUni(imax0+5+i)=utgam2(HbetUni(imax0+5+i)) | |
3055 | GbetpUni(imax0+1+i)=utgam2(HbetpUni(imax0+1+i)) | |
3056 | GbetpUni(imax0+3+i)=utgam2(HbetpUni(imax0+3+i)) | |
3057 | GbetpUni(imax0+5+i)=utgam2(HbetpUni(imax0+5+i)) | |
3058 | HalpUni(imax0+1+i)=ztUni*alpUni(i,1) | |
3059 | HalpUni(imax0+3+i)=zpUni*alpUni(i,1) | |
3060 | HalpUni(imax0+5+i)=zpUni*ztUni*alpUni(i,1) | |
3061 | enddo | |
3062 | ||
3063 | w=0.d0 | |
3064 | xhm=.5d0*(1d0-eps) | |
3065 | yhm=.5d0*(1d0-eps) | |
3066 | do m=1,2 | |
3067 | do i=1,3 | |
3068 | c do i=1,7 | |
3069 | x=xhm+dble((2*m-3)*x9(i))*xhm | |
3070 | c write(ifmt,*)'HPhiInt, xp int :',x | |
3071 | do n=1,2 | |
3072 | do j=1,3 | |
3073 | c do j=1,7 | |
3074 | y=yhm+dble((2*n-3)*x9(j))*yhm | |
3075 | w=w+dble(a9(i)*a9(j))*Hrst(s,b,x,y) | |
3076 | & *PhiUnit(x,y) | |
3077 | c & *PhiExact(1.,x,y,s,b) | |
3078 | enddo | |
3079 | enddo | |
3080 | enddo | |
3081 | enddo | |
3082 | ||
3083 | HPhiInt=w*xhm*yhm | |
3084 | ||
3085 | ||
3086 | c w=0.d0 | |
3087 | c xhm=.5d0*eps | |
3088 | c yhm=.5d0*eps | |
3089 | c do m=1,2 | |
3090 | c do i=1,7 | |
3091 | c x=1d0-eps+xhm+dble((2*m-3)*x1(i))*xhm | |
3092 | c do n=1,2 | |
3093 | c do j=1,7 | |
3094 | c y=1d0-epsyhm+dble((2*n-3)*x1(j))*yhm | |
3095 | c zp2=1.d0-x**4 | |
3096 | c zm2=1.d0-y**4 | |
3097 | c w=w+dble(a1(i)*a1(j))*HrstI(s,x,y) | |
3098 | cc & *PhiUnit(zp2,zm2) | |
3099 | c & *PhiExact(1.,zp2,zm2,s,b) | |
3100 | c enddo | |
3101 | c enddo | |
3102 | c enddo | |
3103 | c enddo | |
3104 | c | |
3105 | c HPhiInt=HPhiInt+w*xhm*yhm | |
3106 | ||
3107 | return | |
3108 | end | |
3109 | ||
3110 | ||
3111 | ||
3112 | c---------------------------------------------------------------------- | |
3113 | subroutine Kfit(iiclegy) | |
3114 | c---------------------------------------------------------------------- | |
3115 | include "epos.inc" | |
3116 | include "epos.incsem" | |
3117 | double precision Znorm | |
3118 | parameter(nbkbin=40) | |
3119 | common /kfitd/ xkappafit(nclegy,nclha,nclha,nbkbin),xkappa,bkbin | |
3120 | parameter (nmax=30) | |
3121 | logical lnoch | |
3122 | ||
3123 | if(iiclegy.eq.-1.or.iiclegy.gt.iclegy2)then | |
3124 | do iiiegy=1,nclegy | |
3125 | do iiipro=1,nclha | |
3126 | do iiitar=1,nclha | |
3127 | do iiibk=1,nbkbin | |
3128 | xkappafit(iiiegy,iiipro,iiitar,iiibk)=1. | |
3129 | enddo | |
3130 | enddo | |
3131 | enddo | |
3132 | enddo | |
3133 | ||
3134 | else | |
3135 | ||
3136 | if(isetcs.le.1)then | |
3137 | s=engy*engy | |
3138 | eps=0.05 | |
3139 | else | |
3140 | s=(egylow*egyfac**(iiclegy-1))**2. | |
3141 | eps=0.001 | |
3142 | endif | |
3143 | ||
3144 | write(ifmt,*)"Fit xkappa ..." | |
3145 | if(ish.ge.2)then | |
3146 | write(ifmt,*)"Kfit s,bkbin,iclegy,ipro,itar" | |
3147 | * ,s,bkbin,iiclegy,iclpro,icltar | |
3148 | endif | |
3149 | ||
3150 | ||
3151 | b=0. | |
3152 | if(isetcs.le.1.or.iiclegy.eq.iclegy2)then | |
3153 | xkf=1. | |
3154 | else | |
3155 | xkf=xkappafit(iiclegy+1,iclpro,icltar,1) | |
3156 | endif | |
3157 | delta=0. | |
3158 | ||
3159 | ||
3160 | do 5 ib=1,nbkbin-1 | |
3161 | b=real(ib-1)*bkbin | |
3162 | xkappafit(iiclegy,iclpro,icltar,ib)=1. | |
3163 | if(b.gt.3.+0.05*log(s).or.s.le.20.*q2min)then | |
3164 | xkf=1. | |
3165 | goto 5 | |
3166 | endif | |
3167 | if(ib.gt.1.and.ish.ge.3)write(ifch,*)" End",delta,xkf | |
3168 | delta=1.-real(Znorm(s,b)) | |
3169 | if(delta.le.0d0)then | |
3170 | if(xkf.ne.1.)then | |
3171 | xkappafit(iiclegy,iclpro,icltar,ib)=xkf | |
3172 | delta=1.-real(Znorm(s,b)) | |
3173 | endif | |
3174 | else!if(xkf.ne.1.)then | |
3175 | goto 5 | |
3176 | endif | |
3177 | if(abs(delta).lt.eps)then | |
3178 | if(delta.lt.0d0)then | |
3179 | xkfs=xkf-delta | |
3180 | deltas=delta | |
3181 | endif | |
3182 | xkf=1. | |
3183 | goto 5 | |
3184 | elseif(ib.le.nbkbin-1)then | |
3185 | ||
3186 | if(delta.gt.0.d0)then | |
3187 | xkf0=1. | |
3188 | xkf1=xkf | |
3189 | delta0=delta | |
3190 | xkf2=xkf-delta0 | |
3191 | xkappafit(iiclegy,iclpro,icltar,ib)=xkf2 | |
3192 | delta1=1.-real(Znorm(s,b)) | |
3193 | if(delta1.lt.0.d0)then | |
3194 | xkf0=xkf2 | |
3195 | xkf1=xkf | |
3196 | delta=delta1 | |
3197 | xkf=xkf0 | |
3198 | else | |
3199 | xkf1=max(delta0,xkf2) | |
3200 | xkf0=0. | |
3201 | xkf=xkf1 | |
3202 | endif | |
3203 | else | |
3204 | xkf0=xkf | |
3205 | xkf1=1.-delta | |
3206 | xkf2=xkf | |
3207 | delta1=delta | |
3208 | endif | |
3209 | ||
3210 | if(ib.eq.1)then | |
3211 | deltas=delta | |
3212 | xkfs=max(0.00001,1.-delta) | |
3213 | endif | |
3214 | ||
3215 | if(delta.le.deltas)xkf=xkfs | |
3216 | if(ish.ge.3)write(ifch,*)" Start",ib,b,delta,xkf,xkf0,xkf1 | |
3217 | if(xkf.eq.xkf2)delta=delta1 | |
3218 | ||
3219 | n=0 | |
3220 | delta0=delta | |
3221 | lnoch=.true. | |
3222 | 10 continue | |
3223 | n=n+1 | |
3224 | if(n.le.nmax.and.xkf1.ne.xkf0)then | |
3225 | if(abs(xkf-xkf2).gt.1e-6.or.abs(delta).gt.abs(deltas))then | |
3226 | xkappafit(iiclegy,iclpro,icltar,ib)=xkf | |
3227 | delta=1.-real(Znorm(s,b)) | |
3228 | endif | |
3229 | if(ish.ge.5)write(ifch,*)" step",ib,n,delta,xkf,delta0 | |
3230 | if(delta*delta0.ge.0.)then | |
3231 | if(lnoch.and.abs(delta).gt.abs(delta0))goto 5 | |
3232 | else | |
3233 | lnoch=.false. | |
3234 | endif | |
3235 | if(abs(delta).gt.eps)then | |
3236 | if(delta.gt.0.)then | |
3237 | xkf1=xkf | |
3238 | xkf=(xkf1+xkf0)*0.5 | |
3239 | delta0=delta | |
3240 | else | |
3241 | xkf0=xkf | |
3242 | xkf=(xkf1+xkf0)*0.5 | |
3243 | delta0=delta | |
3244 | endif | |
3245 | goto 10 | |
3246 | endif | |
3247 | else | |
3248 | if(ish.ge.2) | |
3249 | * write(ifmt,*)"Warning in Kfit, nmax reached : xkappafit=1." | |
3250 | xkappafit(iiclegy,iclpro,icltar,ib)=xkf | |
3251 | endif | |
3252 | endif | |
3253 | ||
3254 | 5 continue | |
3255 | ||
3256 | if(ish.ge.3)write(ifch,*)" End",delta,xkf | |
3257 | 100 if(xkf.gt.1.+eps)write(ifmt,*) | |
3258 | * "Warning in Kfit, xkappafit not yet 1" | |
3259 | xkappafit(iiclegy,iclpro,icltar,nbkbin)=1. | |
3260 | ||
3261 | endif | |
3262 | ||
3263 | return | |
3264 | end | |
3265 | ||
3266 | c---------------------------------------------------------------------- | |
3267 | double precision function Znorm(s,b) !---inu--- | |
3268 | c---------------------------------------------------------------------- | |
3269 | include 'epos.inc' | |
3270 | common /kwrite/ xkapZ | |
3271 | double precision HPhiInt,PhiUnit!,PhiExact | |
3272 | ||
3273 | c write(ifmt,*)'Z calculation for (s,b) :',s,b | |
3274 | imax=idxD1 | |
3275 | if(iomega.eq.2)imax=1 | |
3276 | do i=idxDmin,imax | |
3277 | call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
3278 | call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
3279 | enddo | |
3280 | call GfomPar(b,s) | |
3281 | Znorm=HPhiInt(s,b) | |
3282 | c write(ifch,*)'int',Znorm,' phi',PhiExact(1.,1.d0,1.d0,s,b) | |
3283 | Znorm=Znorm | |
3284 | & +PhiUnit(1.d0,1.d0) | |
3285 | c & +PhiExact(1.,1.d0,1.d0,s,b) | |
3286 | ||
3287 | !write(ifmt,*)'Z=',Znorm,xkapZ,b | |
3288 | return | |
3289 | end | |
3290 | ||
3291 | ||
3292 | c------------------------------------------------------------ | |
3293 | double precision function gammag(iclrem,x) !---test--- | |
3294 | c---------------------------------------------------------------------- | |
3295 | include 'epos.inc' | |
3296 | include 'epos.incsem' | |
3297 | double precision x,utgam2 | |
3298 | ||
3299 | gammag=utgam2(dble(alplea(iclrem))+1.D0) | |
3300 | & /utgam2(dble(alplea(iclrem))+1.D0+x) | |
3301 | ||
3302 | return | |
3303 | end | |
3304 | ||
3305 | ||
3306 | ||
3307 | cc---------------------------------------------------------------------- | |
3308 | c double precision function PomNbri(iqq) !---xsigma--- | |
3309 | cc---------------------------------------------------------------------- | |
3310 | cc integral d2b om1intbci | |
3311 | cc iqq, Pomeron type | |
3312 | cc---------------------------------------------------------------------- | |
3313 | c include 'epos.inc' | |
3314 | c double precision om1intbci | |
3315 | c common/geom/rmproj,rmtarg,bmax,bkmx | |
3316 | c common /ar3/ x1(7),a1(7) | |
3317 | c | |
3318 | c bmid=bkmx/2. | |
3319 | c PomNbri=0.d0 | |
3320 | c do i=1,7 | |
3321 | c do m=1,2 | |
3322 | c bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3323 | c PomNbri=PomNbri+dble(bb*a1(i))*om1intbci(bb,iqq) | |
3324 | c enddo | |
3325 | c enddo | |
3326 | c | |
3327 | c PomNbri=PomNbri*dble(2.*pi*bmid) | |
3328 | c | |
3329 | c return | |
3330 | c end | |
3331 | c | |
3332 | c | |
3333 | c | |
3334 | ||
3335 | ||
3336 | c#################################################################################### | |
3337 | c############# former chk ######################################################### | |
3338 | c#################################################################################### | |
3339 | ||
3340 | ||
3341 | ||
3342 | ||
3343 | ||
3344 | ||
3345 | ||
3346 | ||
3347 | cc---------------------------------------------------------------------- | |
3348 | c double precision function PomIncII(b) !---check--- | |
3349 | cc---------------------------------------------------------------------- | |
3350 | cc integral_dx_dy om1*F_remn*F_remn for a given b !---check--- | |
3351 | cc---------------------------------------------------------------------- | |
3352 | c include 'epos.inc' | |
3353 | c include 'epos.incems' | |
3354 | c include 'epos.incsem' | |
3355 | c include 'epos.incpar' | |
3356 | c double precision cint,gamom(idxD0:idxD1),deltap(idxD0:idxD1) | |
3357 | c &,deltapp(idxD0:idxD1),utgam2 | |
3358 | c | |
3359 | cc Calculation by analytical integration (perfect but it changes | |
3360 | cc if om1 change): | |
3361 | c | |
3362 | c s=engy**2 | |
3363 | c imax=1 | |
3364 | c if(iomega.eq.2)imax=1 | |
3365 | c do i=idxDmin,imax | |
3366 | c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3367 | c gamom(i)=dble(alp*gamv)*chad(iclpro)*chad(icltar) | |
3368 | c deltap(i)=dble(bet) | |
3369 | c deltapp(i)=dble(betp) | |
3370 | c | |
3371 | cc Integration possible only if delta(i)>-1 | |
3372 | c | |
3373 | c if(deltap(i).le.-1.d0.or.deltapp(i).le.-1.d0) | |
3374 | c & stop 'Error in epos-par-300 in PomIncII' | |
3375 | c enddo | |
3376 | c | |
3377 | c cint=0.d0 | |
3378 | c do i=idxDmin,imax | |
3379 | c cint=cint+gamom(i)*utgam2(deltap(i)+1.d0)*utgam2(deltapp(i)+1.d0) | |
3380 | c & *dble(ucfpro*ucftar) | |
3381 | c & /utgam2(dble(alplea(iclpro))+deltap(i)+2.d0) | |
3382 | c & /utgam2(dble(alplea(icltar))+deltapp(i)+2.d0) | |
3383 | c enddo | |
3384 | c | |
3385 | c PomIncII=cint | |
3386 | c | |
3387 | c return | |
3388 | c end | |
3389 | c | |
3390 | ||
3391 | c---------------------------------------------------------------------- | |
3392 | double precision function PomIncXIExact(x) !---check--- | |
3393 | c---------------------------------------------------------------------- | |
3394 | c integral d2b PomIncXExact | |
3395 | c---------------------------------------------------------------------- | |
3396 | include 'epos.inc' | |
3397 | double precision x,PomIncXExact | |
3398 | common /ar3/ x1(7),a1(7) | |
3399 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3400 | ||
3401 | bmid=bkmx/2. | |
3402 | PomIncXIExact=0.d0 | |
3403 | do i=1,7 | |
3404 | do m=1,2 | |
3405 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3406 | PomIncXIExact=PomIncXIExact+dble(bb*a1(i))*PomIncXExact(x,bb) | |
3407 | enddo | |
3408 | enddo | |
3409 | ||
3410 | PomIncXIExact=PomIncXIExact*dble(2.*pi*bmid) | |
3411 | ||
3412 | return | |
3413 | end | |
3414 | ||
3415 | c---------------------------------------------------------------------- | |
3416 | double precision function PomIncXIUnit(x) !---check--- | |
3417 | c---------------------------------------------------------------------- | |
3418 | c integral d2b PomIncXUnit | |
3419 | c---------------------------------------------------------------------- | |
3420 | include 'epos.inc' | |
3421 | double precision x,PomIncXUnit | |
3422 | common /ar3/ x1(7),a1(7) | |
3423 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3424 | ||
3425 | bmid=bkmx/2. | |
3426 | PomIncXIUnit=0.d0 | |
3427 | do i=1,7 | |
3428 | do m=1,2 | |
3429 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3430 | PomIncXIUnit=PomIncXIUnit+dble(bb*a1(i))*PomIncXUnit(x,bb) | |
3431 | enddo | |
3432 | enddo | |
3433 | ||
3434 | PomIncXIUnit=PomIncXIUnit*dble(2.*pi*bmid) | |
3435 | ||
3436 | return | |
3437 | end | |
3438 | ||
3439 | c---------------------------------------------------------------------- | |
3440 | double precision function PomIncPIExact(x) !---check--- | |
3441 | c---------------------------------------------------------------------- | |
3442 | c integral d2b PomIncPExact | |
3443 | c---------------------------------------------------------------------- | |
3444 | include 'epos.inc' | |
3445 | double precision x,PomIncPExact | |
3446 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3447 | common /ar3/ x1(7),a1(7) | |
3448 | ||
3449 | bmid=bkmx/2. | |
3450 | PomIncPIExact=0.d0 | |
3451 | do i=1,7 | |
3452 | do m=1,2 | |
3453 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3454 | PomIncPIExact=PomIncPIExact+dble(bb*a1(i))*PomIncPExact(x,bb) | |
3455 | enddo | |
3456 | enddo | |
3457 | ||
3458 | PomIncPIExact=PomIncPIExact*dble(2.*pi*bmid) | |
3459 | ||
3460 | return | |
3461 | end | |
3462 | ||
3463 | c---------------------------------------------------------------------- | |
3464 | double precision function PomIncPIUnit(x) !---check--- | |
3465 | c---------------------------------------------------------------------- | |
3466 | c integral d2b PomIncPUnit | |
3467 | c---------------------------------------------------------------------- | |
3468 | include 'epos.inc' | |
3469 | double precision x,PomIncPUnit | |
3470 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3471 | common /ar3/ x1(7),a1(7) | |
3472 | ||
3473 | bmid=bkmx/2. | |
3474 | PomIncPIUnit=0.d0 | |
3475 | do i=1,7 | |
3476 | do m=1,2 | |
3477 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3478 | PomIncPIUnit=PomIncPIUnit+dble(bb*a1(i))*PomIncPUnit(x,bb) | |
3479 | enddo | |
3480 | enddo | |
3481 | ||
3482 | PomIncPIUnit=PomIncPIUnit*dble(2.*pi*bmid) | |
3483 | ||
3484 | return | |
3485 | end | |
3486 | ||
3487 | c---------------------------------------------------------------------- | |
3488 | double precision function PomIncMIExact(x) !---check--- | |
3489 | c---------------------------------------------------------------------- | |
3490 | c integral d2b PomIncMExact | |
3491 | c---------------------------------------------------------------------- | |
3492 | include 'epos.inc' | |
3493 | double precision x,PomIncMExact | |
3494 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3495 | common /ar3/ x1(7),a1(7) | |
3496 | ||
3497 | bmid=bkmx/2. | |
3498 | PomIncMIExact=0.d0 | |
3499 | do i=1,7 | |
3500 | do m=1,2 | |
3501 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3502 | PomIncMIExact=PomIncMIExact+dble(bb*a1(i))*PomIncMExact(x,bb) | |
3503 | enddo | |
3504 | enddo | |
3505 | ||
3506 | PomIncMIExact=PomIncMIExact*dble(2.*pi*bmid) | |
3507 | ||
3508 | return | |
3509 | end | |
3510 | ||
3511 | c---------------------------------------------------------------------- | |
3512 | double precision function PomIncMIUnit(x) !---check--- | |
3513 | c---------------------------------------------------------------------- | |
3514 | c integral d2b PomIncMUnit | |
3515 | c---------------------------------------------------------------------- | |
3516 | include 'epos.inc' | |
3517 | double precision x,PomIncMUnit | |
3518 | common/geom/rmproj,rmtarg,bmax,bkmx | |
3519 | common /ar3/ x1(7),a1(7) | |
3520 | ||
3521 | bmid=bkmx/2. | |
3522 | PomIncMIUnit=0.d0 | |
3523 | do i=1,7 | |
3524 | do m=1,2 | |
3525 | bb=bmid*(1.+(2.*m-3)*x1(i)) | |
3526 | PomIncMIUnit=PomIncMIUnit+dble(bb*a1(i))*PomIncMUnit(x,bb) | |
3527 | enddo | |
3528 | enddo | |
3529 | ||
3530 | PomIncMIUnit=PomIncMIUnit*dble(2.*pi*bmid) | |
3531 | ||
3532 | return | |
3533 | end | |
3534 | ||
3535 | c---------------------------------------------------------------------- | |
3536 | double precision function PomIncMExact(xm,b) !---check--- | |
3537 | c---------------------------------------------------------------------- | |
3538 | c incluse Pomeron distribution \int dx+ { 2G F_remn F_remn } | |
3539 | c---------------------------------------------------------------------- | |
3540 | include 'epos.inc' | |
3541 | include 'epos.incsem' | |
3542 | include 'epos.incems' | |
3543 | double precision AlTiP,BeTip,al,bep,bepp,xpInt,utgam2,xm | |
3544 | ||
3545 | s=engy**2 | |
3546 | PomIncMExact=0.d0 | |
3547 | imax=1 | |
3548 | if(iomega.eq.2)imax=1 | |
3549 | do i=idxDmin,imax | |
3550 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3551 | bep =dble(bet) | |
3552 | bepp=dble(betp) | |
3553 | al =dble(alp*gamv) | |
3554 | ||
3555 | BeTip=bep+1.d0 | |
3556 | xpInt=utgam2(BeTip)*dble(ucfpro) | |
3557 | * /utgam2(1.d0+dble(alplea(iclpro))+BeTip) | |
3558 | AlTiP=al*xpInt | |
3559 | PomIncMExact=PomIncMExact+AlTiP*xm**bepp | |
3560 | * *(1.d0-xm)**dble(alplea(icltar)) | |
3561 | enddo | |
3562 | ||
3563 | return | |
3564 | end | |
3565 | ||
3566 | c---------------------------------------------------------------------- | |
3567 | double precision function PomIncMUnit(xm,b) !---check--- | |
3568 | c---------------------------------------------------------------------- | |
3569 | c incluse Unitarized Pomeron distribution \int dx+ | |
3570 | c---------------------------------------------------------------------- | |
3571 | include 'epos.inc' | |
3572 | include 'epos.incsem' | |
3573 | include 'epos.incems' | |
3574 | double precision xh,Df,xp,xm,G2,w,xpm | |
3575 | double precision PoInU!,Znorm | |
3576 | common /ar3/ x1(7),a1(7) | |
3577 | ||
3578 | s=engy**2 | |
3579 | ||
3580 | c Calculation by numeric integration : | |
3581 | w=0.d0 | |
3582 | xpm=.5d0 | |
3583 | do m=1,2 | |
3584 | do j=1,7 | |
3585 | xp=xpm*(1.d0+dble((2.*m-3.)*x1(j))) | |
3586 | xh=xp*xm | |
3587 | ctp060829 sy=s*sngl(xh) | |
3588 | Df=0.d0 | |
3589 | do i=idxDmin,idxD1 | |
3590 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3591 | Df=Df+dble(alp)*xp**dble(bet)*xm**dble(betp) | |
3592 | enddo | |
3593 | call GfomPar(b,s) | |
3594 | G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom) | |
3595 | w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2 | |
3596 | enddo | |
3597 | enddo | |
3598 | w=w*xpm | |
3599 | ||
3600 | ||
3601 | PomIncMUnit=w!/Znorm(s,b) | |
3602 | ||
3603 | return | |
3604 | end | |
3605 | ||
3606 | ||
3607 | c---------------------------------------------------------------------- | |
3608 | double precision function PomIncPExact(xp,b) !---check--- | |
3609 | c---------------------------------------------------------------------- | |
3610 | c incluse Pomeron distribution \int dx- { 2G F_remn F_remn } | |
3611 | c---------------------------------------------------------------------- | |
3612 | include 'epos.inc' | |
3613 | include 'epos.incsem' | |
3614 | include 'epos.incems' | |
3615 | double precision AlTiP,BeTipp,al,bep,bepp,xmInt,utgam2,xp | |
3616 | ||
3617 | s=engy**2 | |
3618 | PomIncPExact=0.d0 | |
3619 | imax=1 | |
3620 | if(iomega.eq.2)imax=1 | |
3621 | do i=idxDmin,imax | |
3622 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3623 | bep=dble(bet) | |
3624 | bepp=dble(betp) | |
3625 | al=dble(alp*gamv) | |
3626 | BeTipp=bepp+1.d0 | |
3627 | xmInt=utgam2(BeTipp)*dble(ucftar) | |
3628 | * /utgam2(1.d0+dble(alplea(icltar))+BeTipp) | |
3629 | AlTiP=al*xmInt | |
3630 | PomIncPExact=PomIncPExact+AlTiP*xp**bep | |
3631 | * *(1.d0-xp)**dble(alplea(iclpro)) | |
3632 | enddo | |
3633 | ||
3634 | return | |
3635 | end | |
3636 | ||
3637 | c---------------------------------------------------------------------- | |
3638 | double precision function PomIncPUnit(xp,b) !---check--- | |
3639 | c---------------------------------------------------------------------- | |
3640 | c incluse Unitarized Pomeron distribution \int dx- | |
3641 | c---------------------------------------------------------------------- | |
3642 | include 'epos.inc' | |
3643 | include 'epos.incsem' | |
3644 | double precision xh,Df,xp,xm,G2,w,xmm | |
3645 | double precision PoInU!,Znorm | |
3646 | common /ar3/ x1(7),a1(7) | |
3647 | ||
3648 | s=engy**2 | |
3649 | imax=1 | |
3650 | if(iomega.eq.2)imax=1 | |
3651 | ||
3652 | c Calculation by numeric integration : | |
3653 | w=0.d0 | |
3654 | xmm=.5d0 | |
3655 | do m=1,2 | |
3656 | do j=1,7 | |
3657 | xm=xmm*(1.d0+dble((2.*m-3.)*x1(j))) | |
3658 | xh=xp*xm | |
3659 | ctp060829 sy=s*real(xh) | |
3660 | Df=0.d0 | |
3661 | do i=idxDmin,imax | |
3662 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3663 | Df=Df+alp*real(xp)**bet*real(xm)**betp | |
3664 | enddo | |
3665 | call GfomPar(b,s) | |
3666 | G2=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom) | |
3667 | w=w+dble(a1(j))*PoInU(xp,xm,s,b)*G2 | |
3668 | enddo | |
3669 | enddo | |
3670 | w=w*xmm | |
3671 | ||
3672 | ||
3673 | PomIncPUnit=w!/Znorm(s,b) | |
3674 | ||
3675 | return | |
3676 | end | |
3677 | ||
3678 | ||
3679 | ||
3680 | c---------------------------------------------------------------------- | |
3681 | double precision function PomIncJExact(b) !---check--- | |
3682 | c---------------------------------------------------------------------- | |
3683 | c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn } | |
3684 | c---------------------------------------------------------------------- | |
3685 | include 'epos.inc' | |
3686 | double precision allea,PomIncXExact,xh | |
3687 | common /ar3/ x1(7),a1(7) | |
3688 | ||
3689 | allea=2.d0+dble(alplea(iclpro)+alplea(icltar)) | |
3690 | PomIncJExact=0.d0 | |
3691 | do i=1,7 | |
3692 | do m=1,2 | |
3693 | xh=1.d0-(.5d0+dble(x1(i)*(real(m)-1.5)))**(1.d0/allea) | |
3694 | PomIncJExact=PomIncJExact+dble(a1(i)) | |
3695 | & *PomIncXExact(xh,b)/(1.d0-xh)**(allea-1.d0) | |
3696 | enddo | |
3697 | enddo | |
3698 | PomIncJExact=PomIncJExact/allea/2.d0 | |
3699 | ||
3700 | return | |
3701 | end | |
3702 | ||
3703 | ||
3704 | c---------------------------------------------------------------------- | |
3705 | double precision function PomIncJUnit(b) !---check--- | |
3706 | c---------------------------------------------------------------------- | |
3707 | c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn } | |
3708 | c---------------------------------------------------------------------- | |
3709 | include 'epos.inc' | |
3710 | double precision PomIncXUnit,xh,xhm | |
3711 | common /ar3/ x1(7),a1(7) | |
3712 | ||
3713 | PomIncJUnit=0.d0 | |
3714 | xhm=.5d0 | |
3715 | do i=1,7 | |
3716 | do m=1,2 | |
3717 | xh=xhm*(1.d0+dble(x1(i)*(2.*real(m)-3.))) | |
3718 | PomIncJUnit=PomIncJUnit+dble(a1(i)) | |
3719 | & *PomIncXUnit(xh,b) | |
3720 | enddo | |
3721 | enddo | |
3722 | PomIncJUnit=PomIncJUnit*xhm | |
3723 | ||
3724 | return | |
3725 | end | |
3726 | ||
3727 | ||
3728 | cc---------------------------------------------------------------------- | |
3729 | c double precision function PomIncXExact1(xh,b) !---check--- | |
3730 | cc---------------------------------------------------------------------- | |
3731 | cc incluse Pomeron distribution \int dy { 2G F_remn F_remn } | |
3732 | cc---------------------------------------------------------------------- | |
3733 | c include 'epos.inc' | |
3734 | c include 'epos.incsem' | |
3735 | c double precision xh,Df,xp,xm,w,ymax | |
3736 | c common /ar3/ x1(7),a1(7) | |
3737 | c | |
3738 | c imax=1 | |
3739 | c if(iomega.eq.2)imax=1 | |
3740 | c s=engy**2 | |
3741 | ccctp060829 sy=s*real(xh) | |
3742 | cc Calculation by numeric integration : | |
3743 | c w=0.d0 | |
3744 | c ymax=-.5d0*log(xh) | |
3745 | c do m=1,2 | |
3746 | c do j=1,7 | |
3747 | c xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax) | |
3748 | c xm=xh/xp | |
3749 | c Df=0.d0 | |
3750 | c do i=idxDmin,imax | |
3751 | c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3752 | c Df=Df+alp*gamv*real(xp)**bet*real(xm)**betp | |
3753 | c * *(1.d0-xp)**dble(alplea(iclpro)) | |
3754 | c * *(1.d0-xm)**dble(alplea(icltar)) | |
3755 | c enddo | |
3756 | c w=w+dble(a1(j))*Df | |
3757 | c enddo | |
3758 | c enddo | |
3759 | c w=w*ymax*xh**dble(-alppar) | |
3760 | c | |
3761 | c | |
3762 | c PomIncXExact1=w | |
3763 | c | |
3764 | c return | |
3765 | c end | |
3766 | ||
3767 | c---------------------------------------------------------------------- | |
3768 | double precision function PomIncXExact(xh,b) !---check--- | |
3769 | c---------------------------------------------------------------------- | |
3770 | c incluse Pomeron distribution \int dy { 2G F_remn F_remn } | |
3771 | c (optimized integration but with no y dependance) | |
3772 | c---------------------------------------------------------------------- | |
3773 | include 'epos.inc' | |
3774 | include 'epos.incsem' | |
3775 | double precision AlTiP,bep,bepp,factor,factor1 | |
3776 | double precision xpmin,xh,xp,xm,ymax,y | |
3777 | common /ar3/ x1(7),a1(7) | |
3778 | ||
3779 | imax=1 | |
3780 | if(iomega.eq.2)imax=1 | |
3781 | s=engy**2 | |
3782 | PomIncXExact=0.d0 | |
3783 | do i=idxDmin,imax | |
3784 | call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
3785 | bep =betx | |
3786 | bepp =betpx | |
3787 | AlTiP=alpx*gamv | |
3788 | PomIncXExact=PomIncXExact+AlTiP*xh**((bep+bepp)/2.d0) | |
3789 | enddo | |
3790 | ||
3791 | factor=0.d0 | |
3792 | allea=min(alplea(iclpro),alplea(icltar))+1. | |
3793 | xpmin=max(sqrt(xh),exp(-1.d0)) | |
3794 | do i=1,7 | |
3795 | do m=1,2 | |
3796 | xp=1.d0-(1.d0-xpmin)*(.5d0+dble(x1(i)*(real(m)-1.5))) | |
3797 | * **(1.d0/dble(allea)) | |
3798 | xm=xh/xp | |
3799 | factor=factor+dble(a1(i)) | |
3800 | * *((1.d0-xp)**dble(alplea(iclpro)-allea+1.) | |
3801 | * *(1.d0-xm)**dble(alplea(icltar)) | |
3802 | * +(1.d0-xp)**dble(alplea(icltar)-allea+1.) | |
3803 | * *(1.d0-xm)**dble(alplea(iclpro)))/xp | |
3804 | enddo | |
3805 | enddo | |
3806 | factor=factor*(1.d0-xpmin)**dble(allea)/dble(allea) | |
3807 | ||
3808 | ||
3809 | if(xpmin.gt.1.00001d0*sqrt(xh))then | |
3810 | ymax=-log(xh)-2.d0 | |
3811 | factor1=0.d0 | |
3812 | do i=1,7 | |
3813 | do m=1,2 | |
3814 | y=ymax*dble(x1(i)*(2*m-3)) | |
3815 | xp=sqrt(xh*exp(y)) | |
3816 | xm=xh/xp | |
3817 | factor1=factor1+dble(a1(i))*(1.d0-xp)**dble(alplea(iclpro)) | |
3818 | * *(1.d0-xm)**dble(alplea(icltar)) | |
3819 | enddo | |
3820 | enddo | |
3821 | factor=factor+factor1*ymax | |
3822 | endif | |
3823 | ||
3824 | factor=factor/2.d0 | |
3825 | ||
3826 | PomIncXExact=PomIncXExact*factor | |
3827 | ||
3828 | return | |
3829 | end | |
3830 | ||
3831 | ||
3832 | c---------------------------------------------------------------------- | |
3833 | double precision function PomIncXUnit(xh,b) !---check--- | |
3834 | c---------------------------------------------------------------------- | |
3835 | c incluse Unitarized Pomeron distribution \int dy | |
3836 | c---------------------------------------------------------------------- | |
3837 | include 'epos.inc' | |
3838 | include 'epos.incsem' | |
3839 | double precision xh,Df,xp,xm,w | |
3840 | double precision PoInU,ymax!,Znorm | |
3841 | common /ar3/ x1(7),a1(7) | |
3842 | ||
3843 | imax=1 | |
3844 | if(iomega.eq.2)imax=1 | |
3845 | s=engy**2 | |
3846 | ctp060829 sy=s*real(xh) | |
3847 | c Calculation by numeric integration : | |
3848 | w=0.d0 | |
3849 | ymax=-.5d0*log(xh) | |
3850 | do m=1,2 | |
3851 | do j=1,7 | |
3852 | xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax) | |
3853 | xm=xh/xp | |
3854 | Df=0.d0 | |
3855 | do i=idxDmin,imax | |
3856 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3857 | Df=Df+alp*real(xp)**bet*real(xm)**betp | |
3858 | enddo | |
3859 | call GfomPar(b,s) | |
3860 | Df=Df*(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom) | |
3861 | w=w+dble(a1(j))*Df*PoInU(xp,xm,s,b) | |
3862 | enddo | |
3863 | enddo | |
3864 | w=w*ymax | |
3865 | ||
3866 | ||
3867 | PomIncXUnit=w!/Znorm(s,b) | |
3868 | ||
3869 | return | |
3870 | end | |
3871 | ||
3872 | ||
3873 | ||
3874 | c---------------------------------------------------------------------- | |
3875 | double precision function PoInU(xp,xm,s,b) !---check--- | |
3876 | c---------------------------------------------------------------------- | |
3877 | c Function : PhiU(1-xp,1-xm) + /int(H(z+ + x+,z- + x-)PhiU(z+,z-)dz+dz-) | |
3878 | c---------------------------------------------------------------------- | |
3879 | include 'epos.inc' | |
3880 | double precision xp,xm,zp,zm,zp2,zm2,zpmin,zmmin,deltzp,deltzm | |
3881 | double precision zpm,zmm,w,HrstI,PhiUnit,Hrst,eps!,PhiExact | |
3882 | common /ar3/ x1(7),a1(7) | |
3883 | ||
3884 | eps=1.d-5 | |
3885 | ||
3886 | imax=idxD1 | |
3887 | if(iomega.eq.2)imax=1 | |
3888 | do i=idxDmin,imax | |
3889 | call GfunPar(1,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
3890 | call GfunPar(2,i,b,s,alpx,betx,betpx,epsp,epst,epss,gamv) | |
3891 | enddo | |
3892 | call GfomPar(b,s) | |
3893 | ||
3894 | if (1.d0-xp-eps.gt.0.d0.and.1.d0-xm-eps.gt.0.d0) then | |
3895 | w=0.d0 | |
3896 | zpmin=1.d0-xp-eps | |
3897 | zmmin=1.d0-xm-eps | |
3898 | zpm=.5d0*zpmin | |
3899 | zmm=.5d0*zmmin | |
3900 | do m=1,2 | |
3901 | do i=1,7 | |
3902 | zp=zpm+dble((2*m-3)*x1(i))*zpm | |
3903 | do n=1,2 | |
3904 | do j=1,7 | |
3905 | zm=zmm+dble((2*n-3)*x1(j))*zmm | |
3906 | w=w+dble(a1(i)*a1(j))*Hrst(s,b,zp+xp,zm+xm) | |
3907 | & *PhiUnit(zp,zm) | |
3908 | c & *PhiExact(1.,zp,zm,s,b) | |
3909 | enddo | |
3910 | enddo | |
3911 | enddo | |
3912 | enddo | |
3913 | PoInU=w*zpm*zmm | |
3914 | deltzp=eps | |
3915 | deltzm=eps | |
3916 | else | |
3917 | PoInU=0.d0 | |
3918 | zpmin=0.d0 | |
3919 | zmmin=0.d0 | |
3920 | deltzp=1.d0-xp | |
3921 | deltzm=1.d0-xm | |
3922 | endif | |
3923 | ||
3924 | w=0.d0 | |
3925 | if(abs(deltzp).gt.1.d-10.and.abs(deltzm).gt.1.d-10)then | |
3926 | zpm=.5d0*deltzp | |
3927 | zmm=.5d0*deltzm | |
3928 | do m=1,2 | |
3929 | do i=1,7 | |
3930 | zp=zpmin+zpm+dble((2*m-3)*x1(i))*zpm | |
3931 | do n=1,2 | |
3932 | do j=1,7 | |
3933 | zm=zmmin+zmm+dble((2*n-3)*x1(j))*zmm | |
3934 | zp2=1.d0-xp-zp**2 | |
3935 | zm2=1.d0-xm-zm**2 | |
3936 | w=w+dble(a1(i)*a1(j))*HrstI(s,b,zp,zm) | |
3937 | & *PhiUnit(zp2,zm2) | |
3938 | c & *PhiExact(1.,zp2,zm2,s,b) | |
3939 | enddo | |
3940 | enddo | |
3941 | enddo | |
3942 | enddo | |
3943 | endif | |
3944 | ||
3945 | PoInU=PoInU+w*zpm*zmm | |
3946 | & +PhiUnit(1.d0-xp,1.d0-xm) | |
3947 | c & +PhiExact(1.,1.d0-xp,1.d0-xm,s,b) | |
3948 | ||
3949 | return | |
3950 | end | |
3951 | ||
3952 | c---------------------------------------------------------------------- | |
3953 | double precision function PomIncExact(xp,xm,b) !---check--- | |
3954 | c---------------------------------------------------------------------- | |
3955 | c inclusive Pomeron distribution { 2G F_remn F_remn } | |
3956 | c---------------------------------------------------------------------- | |
3957 | include "epos.inc" | |
3958 | include "epos.incsem" | |
3959 | double precision Df,xp,xm!,xh | |
3960 | ||
3961 | Df=0.d0 | |
3962 | xh=xp*xm | |
3963 | ctp060829 sy=engy**2*real(xh) | |
3964 | s=engy**2 | |
3965 | imax=1 | |
3966 | if(iomega.eq.2)imax=1 | |
3967 | do i=idxDmin,imax | |
3968 | call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv) | |
3969 | Df=Df+alp*gamv*real(xp)**bet*real(xm)**betp | |
3970 | enddo | |
3971 | Df=dble(chad(iclpro)*chad(icltar)) | |
3972 | * *Df | |
3973 | ||
3974 | PomIncExact=Df | |
3975 | * *(1.d0-xp)**dble(alplea(iclpro)) | |
3976 | * *(1.d0-xm)**dble(alplea(icltar)) | |
3977 | ||
3978 | return | |
3979 | end | |
3980 | ||
3981 | c---------------------------------------------------------------------- | |
3982 | double precision function PomIncUnit(xp,xm,b) !---check--- | |
3983 | c---------------------------------------------------------------------- | |
3984 | c inclusive Pomeron distribution { Sum{int{G*Phi} } | |
3985 | c---------------------------------------------------------------------- | |
3986 | include "epos.inc" | |
3987 | include "epos.incpar" | |
3988 | include "epos.incsem" | |
3989 | double precision PoInU,xp,xm,om1,xh,yp | |
3990 | ||
3991 | ||
3992 | xh=xp*xm | |
3993 | yp=0.d0 | |
3994 | if(xm.ne.0.d0)yp=0.5d0*log(xp/xm) | |
3995 | PomIncUnit=om1(xh,yp,b)*PoInU(xp,xm,engy*engy,b) | |
3996 | & *(1.d0+ztUni*xp**betfom)*(1.d0+zpUni*xm**betfom) | |
3997 | ||
3998 | return | |
3999 | end | |
4000 | ||
4001 | ||
4002 | cc---------------------------------------------------------------------- | |
4003 | c double precision function PomIncUnitMC(xp,xm,b) !---check--- | |
4004 | cc---------------------------------------------------------------------- | |
4005 | cc inclusive Pomeron distribution { Sum{int{G*Phi} } | |
4006 | cc---------------------------------------------------------------------- | |
4007 | c include "epos.inc" | |
4008 | c include "epos.incpar" | |
4009 | c include "epos.incsem" | |
4010 | c include "epos.incems" | |
4011 | c parameter(mmax=20) | |
4012 | c double precision Gtab,Phitab,xxpu(mmax),xxmu(mmax) | |
4013 | c double precision Zm,xp,xm,pGtab,Z,omNpcut,xprem,xmrem, | |
4014 | c * sxp,sxm,PhiExpo | |
4015 | c | |
4016 | c PomIncUnitMC=0.d0 | |
4017 | c if(xp.lt.1.d0.and.xm.lt.1.d0)then | |
4018 | c m=10 | |
4019 | c | |
4020 | c sy=engy*engy | |
4021 | c nmcint=2000 | |
4022 | c nmax=nmcint | |
4023 | c do i=1,mmax | |
4024 | c xxpu(i)=0.d0 | |
4025 | c xxmu(i)=0.d0 | |
4026 | c enddo | |
4027 | c xprem=1.d0 | |
4028 | c xmrem=1.d0 | |
4029 | c sxp=xprem-xp | |
4030 | c sxm=xmrem-xm | |
4031 | c | |
4032 | c Gtab=omNpcut(xp,xm,sxp,sxm,b,0) | |
4033 | c Phitab=PhiExpo(1.,sxp,sxm,sy,b) | |
4034 | c Z=Gtab*Phitab | |
4035 | c Zm=0.d0 | |
4036 | c | |
4037 | c do mtmp=2,m | |
4038 | c | |
4039 | c write(*,*)"GPhi",mtmp-1,Zm,Z | |
4040 | c Zm=0.d0 | |
4041 | c n=0 | |
4042 | c | |
4043 | c 10 continue | |
4044 | c n=n+1 | |
4045 | c if(mod(n,1000000).eq.0)write(*,*) | |
4046 | c & "Calculation of PomIncUnit(",mtmp,")->",n | |
4047 | c xxpu(1)=xp | |
4048 | c xxmu(1)=xm | |
4049 | c sxp=xxpu(1) | |
4050 | c sxm=xxmu(1) | |
4051 | c pGtab=1.d0 | |
4052 | c do i=2,mtmp | |
4053 | c rnau=rangen()*real(xprem-sxp) | |
4054 | c xxpu(i)=dble(rnau) | |
4055 | c sxp=sxp+xxpu(i) | |
4056 | c rnau=rangen()*real(xmrem-sxm) | |
4057 | c xxmu(i)=dble(rnau) | |
4058 | c sxm=sxm+xxmu(i) | |
4059 | c enddo | |
4060 | c if(sxp.lt.xprem.and.sxm.lt.xmrem)then | |
4061 | c do i=1,mtmp | |
4062 | c Gtab=omNpcut(xxpu(i),xxmu(i),xprem-sxp,xmrem-sxm,b,0) | |
4063 | c pGtab=pGtab*Gtab | |
4064 | c enddo | |
4065 | c Zm=Zm+pGtab*PhiExpo(1.,xprem-sxp,xmrem-sxm,sy,b) | |
4066 | c endif | |
4067 | c if(n.lt.nmax)goto 10 | |
4068 | c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(mtmp) | |
4069 | c Z=Z+Zm | |
4070 | c enddo | |
4071 | c | |
4072 | c PomIncUnitMC=Z/dble(chad(iclpro)*chad(icltar)) | |
4073 | c endif | |
4074 | c | |
4075 | c return | |
4076 | c end | |
4077 | c | |
4078 | c | |
4079 | cc---------------------------------------------------------------------- | |
4080 | c double precision function PhiMCExact(xp,xm,b) !---check--- | |
4081 | cc---------------------------------------------------------------------- | |
4082 | cc virtual emissions { Sum{int{-GFF} } | |
4083 | cc---------------------------------------------------------------------- | |
4084 | c include "epos.inc" | |
4085 | c include "epos.incpar" | |
4086 | c include "epos.incsem" | |
4087 | c include "epos.incems" | |
4088 | c parameter(mmax=20) | |
4089 | c double precision Gtab,xxpu(mmax),xxmu(mmax) | |
4090 | c double precision Zm,xp,xm,pGtab,Z,om51,sxp,sxm,xh,yp | |
4091 | cc * ,omNpuncut | |
4092 | c | |
4093 | c PhiMCExact=0.d0 | |
4094 | c if(xp.le.1.d0.and.xm.le.1.d0)then | |
4095 | c m=6 | |
4096 | c | |
4097 | c sy=engy*engy | |
4098 | c nmcint=50000 | |
4099 | c nmax=nmcint | |
4100 | c do i=1,mmax | |
4101 | c xxpu(i)=0.d0 | |
4102 | c xxmu(i)=0.d0 | |
4103 | c enddo | |
4104 | c | |
4105 | c Z=xp**dble(alplea(iclpro)) | |
4106 | c * *xm**dble(alplea(icltar)) | |
4107 | c Zm=0.d0 | |
4108 | c | |
4109 | c do mtmp=1,m | |
4110 | c | |
4111 | c write(*,*)"GPhi",mtmp-1,Zm,Z/xp**dble(alplea(iclpro)) | |
4112 | c * /xm**dble(alplea(icltar)) | |
4113 | c Zm=0.d0 | |
4114 | c n=0 | |
4115 | c | |
4116 | c 10 continue | |
4117 | c n=n+1 | |
4118 | c if(mod(n,1000000).eq.0)write(*,*) | |
4119 | c & "Calculation of PhiExact(",mtmp,")->",n | |
4120 | c sxp=0.d0 | |
4121 | c sxm=0.d0 | |
4122 | c pGtab=1.d0 | |
4123 | c do i=1,mtmp | |
4124 | c rnau=rangen()!*real(xp-sxp) | |
4125 | c xxpu(i)=dble(rnau) | |
4126 | c sxp=sxp+xxpu(i) | |
4127 | c rnau=rangen()!*real(xm-sxm) | |
4128 | c xxmu(i)=dble(rnau) | |
4129 | c sxm=sxm+xxmu(i) | |
4130 | c enddo | |
4131 | c if(sxp.lt.xp.and.sxm.lt.xm)then | |
4132 | c do i=1,mtmp | |
4133 | c xh=xxpu(i)*xxmu(i) | |
4134 | c if(abs(xh).gt.1.d-30)then | |
4135 | c yp=0.5d0*log(xxpu(i)/xxmu(i)) | |
4136 | c else | |
4137 | c yp=0.d0 | |
4138 | c endif | |
4139 | c Gtab=2*om51(xh,yp,b,0,4) | |
4140 | cc * +omNpuncut(sy*real(xh),xh,yp,b,1) !om1(xh,yp,b) | |
4141 | c pGtab=pGtab*(-Gtab) | |
4142 | c enddo | |
4143 | c Zm=Zm+pGtab*(xp-sxp)**dble(alplea(iclpro)) | |
4144 | c * *(xm-sxm)**dble(alplea(icltar)) | |
4145 | c endif | |
4146 | c if(n.lt.nmax)goto 10 | |
4147 | c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(m) | |
4148 | c Z=Z+Zm | |
4149 | c enddo | |
4150 | c | |
4151 | c PhiMCExact=Z | |
4152 | c endif | |
4153 | c | |
4154 | c return | |
4155 | c end | |
4156 | ||
4157 | c---------------------------------------------------------------------- | |
4158 | double precision function Gammapp(sy,b,mtmp) !---check--- | |
4159 | c---------------------------------------------------------------------- | |
4160 | include "epos.inc" | |
4161 | include "epos.incpar" | |
4162 | include "epos.incsem" | |
4163 | include "epos.incems" | |
4164 | parameter(mmax=20) | |
4165 | double precision Gtab,xxpu(mmax),xxmu(mmax),PhiExpo | |
4166 | double precision Zm,xp,xm,pGtab,om1,sxp,sxm,xh,yp | |
4167 | ||
4168 | Gammapp=0.d0 | |
4169 | ||
4170 | xp=1.d0 | |
4171 | xm=1.d0 | |
4172 | nmcint=20000 | |
4173 | nmax=nmcint | |
4174 | do i=1,mmax | |
4175 | xxpu(i)=0.d0 | |
4176 | xxmu(i)=0.d0 | |
4177 | enddo | |
4178 | Zm=0.d0 | |
4179 | ||
4180 | n=0 | |
4181 | ||
4182 | 10 continue | |
4183 | n=n+1 | |
4184 | if(mod(n,10000).eq.0)write(*,*) | |
4185 | & "Calculation of Gammapp(",mtmp,")->",n | |
4186 | sxp=0.d0 | |
4187 | sxm=0.d0 | |
4188 | pGtab=1.d0 | |
4189 | do i=1,mtmp | |
4190 | rnau=rangen()!*real(xp-sxp) | |
4191 | xxpu(i)=dble(rnau) | |
4192 | sxp=sxp+xxpu(i) | |
4193 | rnau=rangen()!*real(xm-sxm) | |
4194 | xxmu(i)=dble(rnau) | |
4195 | sxm=sxm+xxmu(i) | |
4196 | enddo | |
4197 | if(sxp.lt.xp.and.sxm.lt.xm)then | |
4198 | do i=1,mtmp | |
4199 | xh=xxpu(i)*xxmu(i) | |
4200 | if(abs(xh).gt.1.d-30)then | |
4201 | yp=0.5d0*log(xxpu(i)/xxmu(i)) | |
4202 | else | |
4203 | yp=0.d0 | |
4204 | endif | |
4205 | Gtab=om1(xh,yp,b) | |
4206 | pGtab=pGtab*Gtab | |
4207 | enddo | |
4208 | Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b) | |
4209 | endif | |
4210 | if(n.lt.nmax)goto 10 | |
4211 | Zm=Zm/dble(nmax)!**2.d0*(xp*xm)**dble(mtmp) | |
4212 | ||
4213 | Gammapp=Zm | |
4214 | ||
4215 | return | |
4216 | end | |
4217 | ||
4218 | cc---------------------------------------------------------------------- | |
4219 | c double precision function GammaMCnew(sy,b,mtmp) !---check--- | |
4220 | cc---------------------------------------------------------------------- | |
4221 | c include "epos.inc" | |
4222 | c include "epos.incpar" | |
4223 | c include "epos.incsem" | |
4224 | c include "epos.incems" | |
4225 | c parameter(mmax=20) | |
4226 | c common /psar7/ delx,alam3p,gam3p | |
4227 | c double precision Gtab,xxpu(mmax),xxmu(mmax) | |
4228 | c double precision Zm,xp,xm,pGtab,omGam,Zmtot, | |
4229 | c * sxp,sxm,PhiExpo!,om1,yp,xh | |
4230 | c | |
4231 | c GammaMCnew=0.d0 | |
4232 | c Zmtot=0.d0 | |
4233 | c xp=1.d0 | |
4234 | c xm=1.d0 | |
4235 | c nmcint=1000 | |
4236 | c nmax=nmcint | |
4237 | c | |
4238 | c | |
4239 | c do i=1,mmax | |
4240 | c xxpu(i)=0.d0 | |
4241 | c xxmu(i)=0.d0 | |
4242 | c enddo | |
4243 | c | |
4244 | c Zm=0.d0 | |
4245 | c | |
4246 | c n=0 | |
4247 | c | |
4248 | c 10 continue | |
4249 | c n=n+1 | |
4250 | c if(mod(n,1000000).eq.0)write(*,*) | |
4251 | c & "Calculation of GammaMCnew(",mtmp,")->",n | |
4252 | c sxp=0.d0 | |
4253 | c sxm=0.d0 | |
4254 | c pGtab=1.d0 | |
4255 | c do i=1,mtmp | |
4256 | c rnau=rangen() | |
4257 | c xxpu(i)=dble(rnau) | |
4258 | c sxp=sxp+xxpu(i) | |
4259 | c rnau=rangen() | |
4260 | c xxmu(i)=dble(rnau) | |
4261 | c sxm=sxm+xxmu(i) | |
4262 | c enddo | |
4263 | c if(sxp.lt.xp.and.sxm.lt.xm)then | |
4264 | c i=0 | |
4265 | c do k=1,mtmp | |
4266 | c i=i+1 | |
4267 | c Gtab=omGam(xxpu(i),xxmu(i),b) !om1(xh,yp,b) | |
4268 | c pGtab=pGtab*Gtab | |
4269 | c enddo | |
4270 | c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b) | |
4271 | c endif | |
4272 | c if(n.lt.nmax)goto 10 | |
4273 | c | |
4274 | c Zmtot=Zmtot+Zm/dble(nmax) | |
4275 | c | |
4276 | c GammaMCnew=Zmtot | |
4277 | c | |
4278 | c return | |
4279 | c end | |
4280 | cc---------------------------------------------------------------------- | |
4281 | c double precision function GammaMC(sy,b,mtmp) !---check--- | |
4282 | cc---------------------------------------------------------------------- | |
4283 | c include "epos.inc" | |
4284 | c include "epos.incpar" | |
4285 | c include "epos.incsem" | |
4286 | c include "epos.incems" | |
4287 | c parameter(mmax=20) | |
4288 | c double precision Gtab,xxpu(mmax),xxmu(mmax) | |
4289 | c double precision Zm,xp,xm,pGtab,om1, | |
4290 | c * sxp,sxm,xh,yp,PhiExpo!,omNpcut | |
4291 | c | |
4292 | c GammaMC=0.d0 | |
4293 | c | |
4294 | c xp=1.d0 | |
4295 | c xm=1.d0 | |
4296 | c nmcint=50000 | |
4297 | c nmax=nmcint | |
4298 | c do i=1,mmax | |
4299 | c xxpu(i)=0.d0 | |
4300 | c xxmu(i)=0.d0 | |
4301 | c enddo | |
4302 | c | |
4303 | c Zm=0.d0 | |
4304 | c | |
4305 | c n=0 | |
4306 | c | |
4307 | c 10 continue | |
4308 | c n=n+1 | |
4309 | c if(mod(n,1000000).eq.0)write(*,*) | |
4310 | c & "Calculation of GammaMC(",mtmp,")->",n | |
4311 | c sxp=0.d0 | |
4312 | c sxm=0.d0 | |
4313 | c pGtab=1.d0 | |
4314 | c do i=1,mtmp | |
4315 | c rnau=rangen()!*real(xp-sxp) | |
4316 | c xxpu(i)=dble(rnau) | |
4317 | c sxp=sxp+xxpu(i) | |
4318 | c rnau=rangen()!*real(xm-sxm) | |
4319 | c xxmu(i)=dble(rnau) | |
4320 | c sxm=sxm+xxmu(i) | |
4321 | c enddo | |
4322 | c if(sxp.lt.xp.and.sxm.lt.xm)then | |
4323 | c do i=1,mtmp | |
4324 | c xh=xxpu(i)*xxmu(i) | |
4325 | c if(abs(xh).gt.1.d-30)then | |
4326 | c yp=0.5d0*log(xxpu(i)/xxmu(i)) | |
4327 | c else | |
4328 | c yp=0.d0 | |
4329 | c endif | |
4330 | c Gtab=om1(xh,yp,b)!omNpcut(xxpu(i),xxmu(i),xp-sxp,xm-sxm,b,0) !om1(xh,yp,b) | |
4331 | c pGtab=pGtab*Gtab | |
4332 | c enddo | |
4333 | c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b) | |
4334 | c endif | |
4335 | c if(n.lt.nmax)goto 10 | |
4336 | c Zm=Zm/dble(nmax)*fctrl(n-mtmp)*facto(mtmp) | |
4337 | c | |
4338 | c GammaMC=Zm | |
4339 | c | |
4340 | c return | |
4341 | c end | |
4342 | c | |
4343 | c | |
4344 | c | |
4345 | ||
4346 | c---------------------------------------------------------------------- | |
4347 | double precision function GammaGauss(sy,b,mtmp) !---check--- | |
4348 | c---------------------------------------------------------------------- | |
4349 | include "epos.inc" | |
4350 | include "epos.incpar" | |
4351 | include "epos.incsem" | |
4352 | include "epos.incems" | |
4353 | parameter(mmax=3) | |
4354 | common /psar7/ delx,alam3p,gam3p | |
4355 | double precision xpmin,xmmin,Gst,zm(mmax),zp(mmax) | |
4356 | *,xpmax,xmmax,zpmin(mmax),zmmin(mmax),zpmax(mmax) | |
4357 | double precision xp,xm,pGtab,omGam,dzp(mmax),Gp1,Gm1,xmin,eps | |
4358 | *,sxp,sxm,PhiExpo,zmmax(mmax),dzm(mmax),Gp2,Gm2,Gp3,Gm3,G0 | |
4359 | c *,PhiExact | |
4360 | common /ar3/ x1(7),a1(7) | |
4361 | c double precision dgx1,dga1 | |
4362 | c common /dga20/ dgx1(10),dga1(10) | |
4363 | ||
4364 | GammaGauss=0.d0 | |
4365 | xp=1.d0 | |
4366 | xm=1.d0 | |
4367 | xmin=1.d-13 | |
4368 | eps=1.d-15 | |
4369 | ||
4370 | if(mtmp.eq.0)then | |
4371 | nmax1=0 | |
4372 | jmax1=0 | |
4373 | nmax2=0 | |
4374 | jmax2=0 | |
4375 | nmax3=0 | |
4376 | jmax3=0 | |
4377 | elseif(mtmp.eq.1)then | |
4378 | nmax1=2 | |
4379 | jmax1=7 | |
4380 | nmax2=0 | |
4381 | jmax2=0 | |
4382 | nmax3=0 | |
4383 | jmax3=0 | |
4384 | elseif(mtmp.eq.2)then | |
4385 | nmax1=2 | |
4386 | jmax1=7 | |
4387 | nmax2=2 | |
4388 | jmax2=7 | |
4389 | nmax3=0 | |
4390 | jmax3=0 | |
4391 | elseif(mtmp.eq.3)then | |
4392 | nmax1=2 | |
4393 | jmax1=7 | |
4394 | nmax2=2 | |
4395 | jmax2=7 | |
4396 | nmax3=2 | |
4397 | jmax3=7 | |
4398 | else | |
4399 | write(*,*)"m not between 0 and 3, return ..." | |
4400 | return | |
4401 | endif | |
4402 | ||
4403 | xpmin=xmin | |
4404 | xmmin=xmin | |
4405 | xpmax=1.d0 | |
4406 | xmmax=1.d0 | |
4407 | do i=1,mmax | |
4408 | zp(i)=0.d0 | |
4409 | zm(i)=0.d0 | |
4410 | dzp(i)=0.d0 | |
4411 | dzm(i)=0.d0 | |
4412 | zmmin(i)=0.d0 | |
4413 | zpmax(i)=0.d0 | |
4414 | zpmin(i)=0.d0 | |
4415 | zmmax(i)=0.d0 | |
4416 | enddo | |
4417 | ||
4418 | G0=1.d0 | |
4419 | ||
4420 | if(mtmp.eq.0)then | |
4421 | sxp=xp | |
4422 | sxm=xm | |
4423 | G0=PhiExpo(1.,sxp,sxm,sy,b) | |
4424 | endif | |
4425 | ||
4426 | ||
4427 | c write(*,*)'x+/-',xmmin,xmmax,xpmin,xpmax | |
4428 | ||
4429 | ||
4430 | dzm(1)=0.d0 | |
4431 | if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.1)then | |
4432 | zmmax(1)=-log(xmmin) | |
4433 | zmmin(1)=-log(xmmax) | |
4434 | if(abs(xmmin-xmin).lt.eps)then | |
4435 | zmmin(1)=-log(min(xmmax,1.d0-xmmin-xmmin)) | |
4436 | zmmax(1)=-log(max(xmmin,1.d0-xmmax-xmmax)) | |
4437 | endif | |
4438 | dzm(1)=(zmmax(1)-zmmin(1))/2.d0 | |
4439 | endif | |
4440 | ||
4441 | dzp(1)=0.d0 | |
4442 | if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.1)then | |
4443 | zpmax(1)=-log(xpmin) | |
4444 | zpmin(1)=-log(xpmax) | |
4445 | if(abs(xpmin-xmin).lt.eps)then | |
4446 | zpmin(1)=-log(min(xpmax,1.d0-xpmin-xpmin)) | |
4447 | zpmax(1)=-log(max(xpmin,1.d0-xpmax-xpmax)) | |
4448 | endif | |
4449 | dzp(1)=(zpmax(1)-zpmin(1))/2.d0 | |
4450 | endif | |
4451 | ||
4452 | c write(*,*)'bornes1=',exp(-zpmax(1)),exp(-zpmin(1)) | |
4453 | c &,exp(-zmmax(1)),exp(-zmmin(1)) | |
4454 | ||
4455 | Gp1=0.d0 | |
4456 | do np1=1,nmax1 | |
4457 | do jp1=1,jmax1 | |
4458 | zp(1)=zpmin(1)+dzp(1)*(1.d0+dble(2.*np1-3.)*dble(x1(jp1))) | |
4459 | Gm1=0.d0 | |
4460 | if(dzm(1).eq.0.d0)then | |
4461 | nmax1=1 | |
4462 | jmax1=1 | |
4463 | endif | |
4464 | do nm1=1,nmax1 | |
4465 | do jm1=1,jmax1 | |
4466 | if(dzm(1).ne.0.d0)then | |
4467 | zm(1)=zmmin(1)+dzm(1)*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1))) | |
4468 | else | |
4469 | zm(1)=zp(1) | |
4470 | endif | |
4471 | ||
4472 | if(mtmp.eq.1)then | |
4473 | sxp=xp | |
4474 | sxm=xm | |
4475 | do i=1,mtmp | |
4476 | sxp=sxp-exp(-zp(i)) | |
4477 | sxm=sxm-exp(-zm(i)) | |
4478 | enddo | |
4479 | pGtab=1.d0 | |
4480 | k=0 | |
4481 | do l=1,mtmp | |
4482 | k=k+1 | |
4483 | if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then | |
4484 | Gst=omGam(exp(-zp(k)),exp(-zm(k)),b) | |
4485 | pGtab=pGtab*Gst | |
4486 | if(Gst.eq.0.d0) | |
4487 | & write(*,*)'j1=',k,exp(-zp(k)),exp(-zm(k)) | |
4488 | & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1 | |
4489 | else | |
4490 | pGtab=0.d0 | |
4491 | write(*,*)'error1 ?',dzp(k),dzm(k) | |
4492 | endif | |
4493 | enddo | |
4494 | if(sxp.gt.0.d0.and.sxm.gt.0.d0)then | |
4495 | if(dzm(1).ne.0.d0)then | |
4496 | Gm1=Gm1+pGtab* | |
4497 | &dble(a1(jm1))*PhiExpo(1.,sxp,sxm,sy,b) | |
4498 | &*exp(-zm(1)) | |
4499 | c &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b) | |
4500 | c &*exp(-zm(1)) | |
4501 | else | |
4502 | Gp1=Gp1+pGtab* | |
4503 | &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b) | |
4504 | &*exp(-zp(1)) | |
4505 | c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b) | |
4506 | c &*exp(-zp(1)) | |
4507 | endif | |
4508 | c write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm | |
4509 | endif | |
4510 | endif | |
4511 | ||
4512 | dzp(2)=0.d0 | |
4513 | if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.2)then | |
4514 | zpmin(2)=-log(min(min(xpmax,1.d0-exp(-zp(1))), | |
4515 | & 1.d0-xpmin-exp(-zp(1)))) | |
4516 | zpmax(2)=-log(max(xpmin,1.d0-xpmax-exp(-zp(1)))) | |
4517 | if(abs(xpmax+xpmax+xpmax-3.d0*dble(1./delx)).lt.eps)then | |
4518 | zpmin(2)=-log(xpmax) | |
4519 | zpmax(2)=-log(xpmin) | |
4520 | endif | |
4521 | dzp(2)=(zpmax(2)-zpmin(2))/2.d0 | |
4522 | endif | |
4523 | ||
4524 | dzm(2)=0.d0 | |
4525 | if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.2)then | |
4526 | zmmin(2)=-log(min(min(xmmax,1.d0-exp(-zm(1))), | |
4527 | & 1.d0-xmmin-exp(-zm(1)))) | |
4528 | zmmax(2)=-log(max(xmmin,1.d0-xmmax-exp(-zm(1)))) | |
4529 | if(abs(xmmax+xmmax+xmmax-3.d0*dble(1./delx)).lt.eps)then | |
4530 | zmmin(2)=-log(xmmax) | |
4531 | zmmax(2)=-log(xmmin) | |
4532 | endif | |
4533 | dzm(2)=(zmmax(2)-zmmin(2))/2.d0 | |
4534 | endif | |
4535 | c write(*,*)'bornes2=',exp(-zpmax(2)),exp(-zpmin(2)) | |
4536 | c &,exp(-zmmax(2)),exp(-zmmin(2)),xpmax(2),1.d0-exp(-zp(1)) | |
4537 | c &,1.d0-xpmin(3)-exp(-zp(1)),xpmin(2),1.d0-xpmax(3)-exp(-zp(1)) | |
4538 | ||
4539 | Gp2=0.d0 | |
4540 | do np2=1,nmax2 | |
4541 | do jp2=1,jmax2 | |
4542 | zp(2)=zpmin(2)+dzp(2)*(1.d0+dble(2.*np2-3.)*dble(x1(jp2))) | |
4543 | Gm2=0.d0 | |
4544 | if(dzm(2).eq.0.d0)then | |
4545 | nmax2=1 | |
4546 | jmax2=1 | |
4547 | endif | |
4548 | do nm2=1,nmax2 | |
4549 | do jm2=1,jmax2 | |
4550 | if(dzm(2).ne.0.d0)then | |
4551 | zm(2)=zmmin(2)+dzm(2)*(1.d0+dble(2.*nm2-3.)*dble(x1(jm2))) | |
4552 | else | |
4553 | zm(2)=zp(2) | |
4554 | endif | |
4555 | ||
4556 | if(mtmp.eq.2)then | |
4557 | sxp=xp | |
4558 | sxm=xm | |
4559 | do i=1,mtmp | |
4560 | sxp=sxp-exp(-zp(i)) | |
4561 | sxm=sxm-exp(-zm(i)) | |
4562 | enddo | |
4563 | pGtab=1.d0 | |
4564 | k=0 | |
4565 | do l=1,mtmp | |
4566 | k=k+1 | |
4567 | if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then | |
4568 | Gst=omGam(exp(-zp(k)),exp(-zm(k)),b) | |
4569 | pGtab=pGtab*Gst | |
4570 | if(Gst.eq.0.d0) | |
4571 | & write(*,*)'j2=',k,exp(-zp(k)),exp(-zm(k)) | |
4572 | & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2 | |
4573 | else | |
4574 | pGtab=0.d0 | |
4575 | write(*,*)'error2 ?',dzp(k),dzm(k) | |
4576 | endif | |
4577 | enddo | |
4578 | if(sxp.gt.0.d0.and.sxm.gt.0.d0)then | |
4579 | if(dzm(2).ne.0.d0)then | |
4580 | Gm2=Gm2+pGtab* | |
4581 | &dble(a1(jm2))*PhiExpo(1.,sxp,sxm,sy,b) | |
4582 | &*exp(-zm(2)) | |
4583 | c &dble(a1(jm2))*PhiExact(1.,sxp,sxm,sy,b,mk) | |
4584 | c &*exp(-zm(2)) | |
4585 | else | |
4586 | Gp2=Gp2+pGtab* | |
4587 | &dble(a1(jp2))*PhiExpo(1.,sxp,sxm,sy,b) | |
4588 | &*exp(-zp(2)) | |
4589 | c &dble(a1(jp2))*PhiExact(1.,sxp,sxm,sy,b,mk) | |
4590 | c &*exp(-zp(2)) | |
4591 | endif | |
4592 | c write(*,*)'m=2',mtmp,Gm2,Gp2,pGtab,sxp,sxm | |
4593 | endif | |
4594 | endif | |
4595 | ||
4596 | dzp(3)=0.d0 | |
4597 | if(abs(xpmin-xpmax).ge.eps.and.mtmp.ge.3)then | |
4598 | zpmin(3)=-log(min(xpmax,1.d0-exp(-zp(1))-exp(-zp(2)))) | |
4599 | zpmax(3)=-log(xpmin) | |
4600 | dzp(3)=(zpmax(3)-zpmin(3))/2.d0 | |
4601 | endif | |
4602 | ||
4603 | dzm(3)=0.d0 | |
4604 | if(abs(xmmin-xmmax).ge.eps.and.mtmp.ge.3)then | |
4605 | zmmin(3)=-log(min(xmmax,1.d0-exp(-zm(1))-exp(-zm(2)))) | |
4606 | zmmax(3)=-log(xmmin) | |
4607 | dzm(3)=(zmmax(3)-zmmin(3))/2.d0 | |
4608 | endif | |
4609 | ||
4610 | c write(*,*)'bornes3=',exp(-zpmax(3)),exp(-zpmin(3)) | |
4611 | c &,exp(-zmmax(3)),exp(-zmmin(3)) | |
4612 | ||
4613 | Gp3=0.d0 | |
4614 | do np3=1,nmax3 | |
4615 | do jp3=1,jmax3 | |
4616 | zp(3)=zpmin(3)+dzp(3)*(1.d0+dble(2.*np3-3.)*dble(x1(jp3))) | |
4617 | Gm3=0.d0 | |
4618 | if(dzm(3).eq.0.d0)then | |
4619 | nmax3=1 | |
4620 | jmax3=1 | |
4621 | endif | |
4622 | do nm3=1,nmax3 | |
4623 | do jm3=1,jmax3 | |
4624 | if(dzm(3).ne.0.d0)then | |
4625 | zm(3)=zmmin(3)+dzm(3)*(1.d0+dble(2.*nm3-3.)*dble(x1(jm3))) | |
4626 | else | |
4627 | zm(3)=zp(3) | |
4628 | endif | |
4629 | ||
4630 | sxp=xp | |
4631 | sxm=xm | |
4632 | do i=1,mtmp | |
4633 | sxp=sxp-exp(-zp(i)) | |
4634 | sxm=sxm-exp(-zm(i)) | |
4635 | enddo | |
4636 | pGtab=1.d0 | |
4637 | k=0 | |
4638 | do l=1,mtmp | |
4639 | k=k+1 | |
4640 | if(dzp(k).ge.0.d0.and.dzm(k).ge.0.d0)then | |
4641 | Gst=omGam(exp(-zp(k)),exp(-zm(k)),b) | |
4642 | pGtab=pGtab*Gst | |
4643 | if(Gst.eq.0.d0) | |
4644 | & write(*,*)'j3=',k,exp(-zp(k)),exp(-zm(k)) | |
4645 | & ,exp(-zpmin(k)),exp(-zpmax(k)),dzp(k),dzm(k),jp1,jp2,jp3 | |
4646 | else | |
4647 | pGtab=0.d0 | |
4648 | write(*,*)'error3 ?',k,dzp(k),dzm(k) | |
4649 | endif | |
4650 | enddo | |
4651 | if(sxp.gt.0.d0.and.sxm.gt.0.d0)then | |
4652 | if(dzm(3).ne.0.d0)then | |
4653 | Gm3=Gm3+pGtab | |
4654 | &*dble(a1(jm3))*PhiExpo(1.,sxp,sxm,sy,b) | |
4655 | &*exp(-zm(3)) | |
4656 | else | |
4657 | Gp3=Gp3+pGtab | |
4658 | &*dble(a1(jp3))*PhiExpo(1.,sxp,sxm,sy,b) | |
4659 | &*exp(-zp(3)) | |
4660 | endif | |
4661 | endif | |
4662 | enddo | |
4663 | enddo | |
4664 | if(dzm(3).ne.0.d0)Gp3=Gp3+Gm3*dble(a1(jp3))*exp(-zp(3))*dzm(3) | |
4665 | nmax3=2 | |
4666 | jmax3=7 | |
4667 | enddo | |
4668 | enddo | |
4669 | if(mtmp.gt.2.and.dzm(2).ne.0.d0)then | |
4670 | Gm2=Gm2+Gp3*dble(a1(jm2))*exp(-zm(2))*dzp(3) | |
4671 | elseif(mtmp.gt.2)then | |
4672 | Gp2=Gp2+Gp3*dble(a1(jp2))*exp(-zp(2))*dzp(3) | |
4673 | endif | |
4674 | enddo | |
4675 | enddo | |
4676 | if(dzm(2).ne.0.d0)Gp2=Gp2+Gm2*dble(a1(jp2))*exp(-zp(2))*dzm(2) | |
4677 | nmax2=2 | |
4678 | jmax2=7 | |
4679 | enddo | |
4680 | enddo | |
4681 | if(mtmp.gt.1.and.dzm(1).ne.0.d0)then | |
4682 | Gm1=Gm1+Gp2*dble(a1(jm1))*exp(-zm(1))*dzp(2) | |
4683 | elseif(mtmp.gt.1)then | |
4684 | Gp1=Gp1+Gp2*dble(a1(jp1))*exp(-zp(1))*dzp(2) | |
4685 | endif | |
4686 | enddo | |
4687 | enddo | |
4688 | if(dzm(1).ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*exp(-zp(1))*dzm(1) | |
4689 | nmax1=2 | |
4690 | jmax1=7 | |
4691 | enddo | |
4692 | enddo | |
4693 | ||
4694 | if(mtmp.gt.0)G0=Gp1*dzp(1) | |
4695 | write(*,*)"int:",G0 | |
4696 | ||
4697 | GammaGauss=GammaGauss+G0 | |
4698 | ||
4699 | return | |
4700 | ||
4701 | end | |
4702 | ||
4703 | c----------------------------------------------------------------------- | |
4704 | double precision function omWi(sy,b) !---check--- | |
4705 | c----------------------------------------------------------------------- | |
4706 | c cut enhanced diagram integrated over xp, xm, xpr,xmr | |
4707 | c (with ideal G) | |
4708 | c b - impact parameter between the pomeron ends; | |
4709 | c sy- total energy | |
4710 | c----------------------------------------------------------------------- | |
4711 | include "epos.inc" | |
4712 | include "epos.incpar" | |
4713 | include "epos.incsem" | |
4714 | include "epos.incems" | |
4715 | common /psar7/ delx,alam3p,gam3p | |
4716 | double precision xpmin,xmmin,zp,zm,alpp,alpm,xjacp,xjacm | |
4717 | *,xpmax,xmmax,zpmin,zmmin,zpmax,chg | |
4718 | double precision xp,xm,pGtab,omGam,dzp,Gp1,Gm1,xmin,eps | |
4719 | *,sxp,sxm,PhiExpo,zmmax,dzm!,gamp,gamm,gampp,gammp | |
4720 | c *,PhiExact | |
4721 | common /ar3/ x1(7),a1(7) | |
4722 | c double precision dgx1,dga1 | |
4723 | c common /dga20/ dgx1(10),dga1(10) | |
4724 | ||
4725 | omWi=0.d0 | |
4726 | ||
4727 | xmin=1.d-30 | |
4728 | eps=1.d-15 | |
4729 | chg=1.d0/dble(delx) | |
4730 | b2=b*b | |
4731 | gamb=gamD(1,iclpro,icltar) | |
4732 | ctp060829 gamp=dble(gamb*b2/2.-alppar) | |
4733 | ctp060829 gamm=dble(gamb*b2/2.-alppar) | |
4734 | ctp060829 gampp=1.d0+2.d0*gamp | |
4735 | ctp060829 gammp=1.d0+2.d0*gamm | |
4736 | ||
4737 | nmax=2 | |
4738 | jmax=7 | |
4739 | ||
4740 | xpmin=xmin | |
4741 | xmmin=xmin | |
4742 | xpmax=1.d0 | |
4743 | xmmax=1.d0 | |
4744 | zpmin=0.d0 | |
4745 | zmmin=0.d0 | |
4746 | zpmax=0.d0 | |
4747 | zmmax=0.d0 | |
4748 | zp=0.d0 | |
4749 | zm=0.d0 | |
4750 | dzp=0.d0 | |
4751 | dzm=0.d0 | |
4752 | ||
4753 | ||
4754 | ||
4755 | do intp=1,2 | |
4756 | do intm=1,2 | |
4757 | ||
4758 | if(intp.eq.1)then | |
4759 | ||
4760 | xpmin=xmin | |
4761 | xpmax=chg | |
4762 | alpp=(1.d0+2.d0*dble(gamb*b2/2.)) | |
4763 | ||
4764 | else | |
4765 | ||
4766 | xpmin=chg | |
4767 | xpmax=1.d0 | |
4768 | alpp=1.d0!(1.d0+2.d0*dble(gamb*b2/2.)) | |
4769 | ||
4770 | endif | |
4771 | ||
4772 | if(intm.eq.1)then | |
4773 | ||
4774 | xmmin=xmin | |
4775 | xmmax=chg | |
4776 | alpm=(1.d0+2.d0*dble(gamb*b2/2.)) | |
4777 | ||
4778 | else | |
4779 | ||
4780 | xmmin=chg | |
4781 | xmmax=1.d0 | |
4782 | alpm=1.d0!(1.d0+2.d0*dble(gamb*b2/2.)) | |
4783 | ||
4784 | endif | |
4785 | c write(*,*)'x+/-',intp,intm,xmmin,xmmax,xpmin,xpmax,alpp,alpm | |
4786 | ||
4787 | ||
4788 | dzm=0.d0 | |
4789 | if(abs(xmmin-xmmax).ge.eps)then | |
4790 | if(alpm.eq.0.d0)then | |
4791 | zmmax=-log(xmmin) | |
4792 | zmmin=-log(xmmax) | |
4793 | else | |
4794 | zmmin=xmmin**alpm | |
4795 | zmmax=xmmax**alpm | |
4796 | endif | |
4797 | dzm=(zmmax-zmmin)/2.d0 | |
4798 | endif | |
4799 | ||
4800 | dzp=0.d0 | |
4801 | if(abs(xpmin-xpmax).ge.eps)then | |
4802 | if(alpp.eq.0.d0)then | |
4803 | zpmax=-log(xpmin) | |
4804 | zpmin=-log(xpmax) | |
4805 | else | |
4806 | zpmin=xpmin**alpp | |
4807 | zpmax=xpmax**alpp | |
4808 | endif | |
4809 | dzp=(zpmax-zpmin)/2.d0 | |
4810 | endif | |
4811 | ||
4812 | ||
4813 | Gp1=0.d0 | |
4814 | ||
4815 | if(abs(dzp).gt.eps.and.abs(dzm).gt.eps)then | |
4816 | c write(*,*)'Ca passe ...' | |
4817 | ||
4818 | do np1=1,nmax | |
4819 | do jp1=1,jmax | |
4820 | zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dble(x1(jp1))) | |
4821 | c zp=zpmin+dzp*(1.d0+dble(2.*np1-3.)*dgx1(jp1)) | |
4822 | if(alpp.eq.0.d0)then | |
4823 | xp=exp(-zp) | |
4824 | xjacp=xp | |
4825 | else | |
4826 | xp=zp**(1.d0/alpp) | |
4827 | xjacp=zp**(1.d0/alpp-1.d0)/alpp | |
4828 | endif | |
4829 | ||
4830 | Gm1=0.d0 | |
4831 | do nm1=1,nmax | |
4832 | do jm1=1,jmax | |
4833 | zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dble(x1(jm1))) | |
4834 | c zm=zmmin+dzm*(1.d0+dble(2.*nm1-3.)*dgx1(jm1)) | |
4835 | if(alpm.eq.0.d0)then | |
4836 | xm=exp(-zm) | |
4837 | xjacm=xm | |
4838 | else | |
4839 | xm=zm**(1.d0/alpm) | |
4840 | xjacm=zm**(1.d0/alpm-1.d0)/alpm | |
4841 | endif | |
4842 | ||
4843 | sxp=1.d0-xp | |
4844 | sxm=1.d0-xm | |
4845 | pGtab=1.d0 | |
4846 | if(dzp.ge.0.d0.and.dzm.ge.0.d0)then | |
4847 | pGtab=omGam(xp,xm,b) | |
4848 | if(pGtab.eq.0.d0) | |
4849 | & write(*,*)'j1=',xp,xm,xmmin,xmmax,dzp,dzm,jp1 | |
4850 | else | |
4851 | pGtab=0.d0 | |
4852 | write(*,*)'error ?',dzp,dzm | |
4853 | endif | |
4854 | if(sxp.gt.0.d0.and.sxm.gt.0.d0)then | |
4855 | if(dzm.ne.0.d0)then | |
4856 | Gm1=Gm1+pGtab* | |
4857 | &dble(a1(jm1))*PhiExpo(1.,sxp,sxm,sy,b)*xjacm | |
4858 | c &dga1(jm1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacm | |
4859 | c &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b)*xjacm | |
4860 | else | |
4861 | Gp1=Gp1+pGtab* | |
4862 | &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b)*xjacp | |
4863 | c &dga1(jp1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacp | |
4864 | c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)*xjacp | |
4865 | endif | |
4866 | c write(*,*)'m=1',mtmp,Gm1,Gp1,pGtab,sxp,sxm | |
4867 | endif | |
4868 | ||
4869 | enddo | |
4870 | enddo | |
4871 | if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dble(a1(jp1))*dzm*xjacp | |
4872 | c if(dzm.ne.0.d0)Gp1=Gp1+Gm1*dga1(jp1)*dzm*xjacp | |
4873 | enddo | |
4874 | enddo | |
4875 | endif | |
4876 | ||
4877 | omWi=omWi+Gp1*dzp | |
4878 | ||
4879 | enddo | |
4880 | enddo | |
4881 | ||
4882 | ||
4883 | return | |
4884 | end | |
4885 | ||
4886 | c----------------------------------------------------------------------- | |
4887 | double precision function Womint(sy,bh) !---check--- | |
4888 | c----------------------------------------------------------------------- | |
4889 | c - chi~(xp,xm)/2. for group of cut enhanced diagram giving | |
4890 | c the same final state integrated over xpr and xmr (with ideal G) | |
4891 | c bh - impact parameter between the pomeron ends; | |
4892 | c xh - fraction of the energy squared s for the pomeron; | |
4893 | c yp - rapidity for the pomeron; | |
4894 | c----------------------------------------------------------------------- | |
4895 | include 'epos.inc' | |
4896 | double precision omWi | |
4897 | ||
4898 | Womint=omWi(sy,bh) | |
4899 | ||
4900 | return | |
4901 | end | |
4902 | ||
4903 | ||
4904 | ||
4905 | c----------------------------------------------------------------------- | |
4906 | double precision function WomGamint(bh) !---check--- | |
4907 | c----------------------------------------------------------------------- | |
4908 | c - chi~(xp,xm)/2. for group of integrated cut enhanced diagram giving | |
4909 | c the same final for proposal. | |
4910 | c bh - impact parameter between the pomeron ends; | |
4911 | c xh - fraction of the energy squared s for the pomeron; | |
4912 | c yp - rapidity for the pomeron; | |
4913 | c----------------------------------------------------------------------- | |
4914 | include 'epos.inc' | |
4915 | double precision omGamint | |
4916 | ||
4917 | WomGamint=omGamint(bh) | |
4918 | ||
4919 | return | |
4920 | end | |
4921 |