]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-omg-160.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of http://git.cern.ch/pub/AliRoot
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-omg-160.f
CommitLineData
9ef1c2d9 1c----------------------------------------------------------------------
2c The two elementary fuctions of our approach, the profile fuction G
3c and the Phi exponent G~, are here referred to as Gpro and Gtilde.
4c Both functions can be written as
5c
6c G = \sum_type alp * xp**bet * xm**betp
7c
8c The parameters alp, bet, betp are calculated in GfunParK (k-mode,
9c b is takento the one of pair k) or GfunPar (b-mode: arbitrary b) as
10c
11c Gpro: bet = betD' + epsGp + gamD*b**2 + epsG -alppar
12c bet' = betD'' + epsGt + gamD*b**2 + epsG -alppar
13c
14c alp = alpD*f * s**(betD+gamD*b**2+epsG) * exp(-b**2/delD)
15c
16c Gtilde: bet~ = bet + 1
17c bet~' = bet' + 1
18c
19c alp~ = alp * gam(bet~) * gam(bet~')
20c * gam(1+alppro) * gam(1+alptar)
21c * gam(1+alppro+bet~) * gam(1+alptar+bet~')
22c * (1+epsGt') * (1+epsGt')
23c
24c The parameters depend implicitely on s.
25c
26c In the program we use om1 = Gpro
27c (they differ by a constant which is actually one)
28c and om5 = om1 * 0.5
29c All functions related to om1 are called om1... .
30c
31c The inclusive Pomeron distributions are
32c
33c PomInc(xp,xm) = Gpro(xp,xm) * (1-xp)**alppro * (1-xm)**alptar
34c
35c----------------------------------------------------------------------
36
37
38c----------------------------------------------------------------------
39 subroutine GfunParK(irea) !---MC---
40c----------------------------------------------------------------------
41c calculates parameters alp,bet,betp of the G functions (k-mode)
42c and the screening exponents epsilongp(k,i), epsilongt(k,i), epsilongs(k,i)
43c----------------------------------------------------------------------
44c Gpro parameters written to /comtilde/atilde(,)btildep(,),btildepp(,)
45c Gtilde parameters written to /cgtilde/atildg(,)btildgp(,),btildgpp(,)
46c two subscripts: first=type, second=collision k
47c Certain pieces of this routine are only done if irea is <= or = zero.
48c----------------------------------------------------------------------
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
118ctp060829 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)
132c epsilongp(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
133c & epscrs*min(epscrx,zkp*abs(betD(0,iclpro,icltar))))
134c epsilongp(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
135c & 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
146c 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)
158c epsilongt(k,0)=max(-betDp(0,iclpro,icltar)-0.95+alppar,
159c & epscrs*min(epscrx,zkt*abs(betD(0,iclpro,icltar))))
160c epsilongt(k,1)=max(-betDp(1,iclpro,icltar)-0.95+alppar,
161c & 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
172c 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
190c 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))
205c & epscrs*x)
206c & epscrs*x*abs(betDp(0,iclpro,icltar)))
207c & 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))))
210cc epsilongp(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar)))
211cc 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))
224c & epscrs*x)
225c & epscrs*x*abs(betDpp(0,iclpro,icltar)))
226c & 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))))
229cc epsilongt(k,0)=epscrs*min(epscrx,x*abs(betD(0,iclpro,icltar)))
230cc epsilongt(k,1)=epscrh*min(epscrx,x*abs(betD(1,iclpro,icltar)))
231cc 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
304c 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
456c----------------------------------------------------------------------
457 subroutine GfunPar(m,i,b,spp,alp,bet,betp,epsp,epst,epss,gamvv)
458c----------------------------------------------------------------------
459c calculates parameters alp,bet,betp of the G functions for pp (b-mode)
460c and the screening exponents epsp,epst,epss,gamvv
461c----------------------------------------------------------------------
462c m=1: profile function Gpro, i=0: soft i=2: diff
463c m=2: Gtilde, i=1: semi (no screening for diff)
464c b: impact param, spp: pp energy squared
465c----------------------------------------------------------------------
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
509ctp060829 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))
519c & ,epscrs*x)
520c & ,epscrs*x*abs(betDp(i,iclpro,icltar)))
521c & ,epscrs*min(epscrx,x*abs(betDp(i,iclpro,icltar))))
522c 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))
534c & ,epscrs*x)
535c & ,epscrs*x*abs(betDpp(i,iclpro,icltar)))
536c & ,epscrs*min(epscrx,x*abs(betDpp(i,iclpro,icltar))))
537c 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
544c 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.
591c 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
618c 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
667c----------------------------------------------------------------------
668 subroutine GfomPar(b,spp)
669c----------------------------------------------------------------------
670c calculates parameters of the fom functions for pp (b-mode)
671c----------------------------------------------------------------------
672c b: impact param, spp: pp energy squared
673c----------------------------------------------------------------------
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
707ctp060829 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)))
723c 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)))
726c 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
737c----------------------------------------------------------------------
738 function fscra(x)
739c----------------------------------------------------------------------
740 fscra=0
741 x2=x*x
742 if(x2.gt.1.)fscra=log(x2)
743 end
744
745c----------------------------------------------------------------------
746 function fscro(x)
747c----------------------------------------------------------------------
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
754c----------------------------------------------------------------------
755 subroutine recalcZPtn !???????????? not updated !!!
756c----------------------------------------------------------------------
757
758 include 'epos.inc'
759 include 'epos.incems'
760 stop'recalcZPtn not valid any more!!!!!!!'
761c if(koll.eq.1.and.maproj.eq.1.and.matarg.eq.1)then
762c npom=nprt(1)
763c k=1
764c ip=iproj(1)
765c it=itarg(1)
766c zparpro(k)=max(0,npom-1)*0.2
767c zpartar(k)=0
768c zpartar(k)=max(0,npom-1)*0.2
769c ztav=zpartar(k)
770c zpav=zparpro(k)
771c zppevt=zpav
772c zptevt=ztav
773c endif
774 end
775
776c----------------------------------------------------------------------
777 double precision function om1(xh,yp,b) !---test---
778c----------------------------------------------------------------------
779c om1 = G * C * gamd (C and gamd usually 1)
780c xh - fraction of the energy squared s for the pomeron;
781c b - impact parameter between the pomeron ends;
782c yp - rapidity for the pomeron;
783c----------------------------------------------------------------------
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
804c----------------------------------------------------------------------
805 double precision function om1intb(b) !---MC---
806c----------------------------------------------------------------------
807c om1 integrated over xp and xm for given b
808c Calculation by analytical integration
809c----------------------------------------------------------------------
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
848c----------------------------------------------------------------------
849 double precision function om1intbi(b,iqq) !---MC---
850c----------------------------------------------------------------------
851c om1 integrated over xp and xm for given b
852c Calculation by analytical integration of contribution iqq
853c----------------------------------------------------------------------
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
885c----------------------------------------------------------------------
886 double precision function om1intbc(b) !---MC---
887c----------------------------------------------------------------------
888c om1*F*F integrated over xp and xm for given b
889c Calculation by analytical integration
890c----------------------------------------------------------------------
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
928cc----------------------------------------------------------------------
929c double precision function om1intbci(b,iqq) !---MC---
930cc----------------------------------------------------------------------
931cc om1*F*F integrated over xp and xm for given b and given Pomeron type iqq
932cc Calculation by analytical integration
933cc----------------------------------------------------------------------
934c include 'epos.inc'
935c include 'epos.incems'
936c include 'epos.incsem'
937c double precision cint,gamom,deltap,deltam
938c double precision utgam2,Fp,Fm,eps
939c
940c spp=engy**2
941c om1intbci=0.d0
942c Fp=dble(ucfpro) !gamma(1+alplea)
943c Fm=dble(ucftar)
944c
945c i=iqq
946c call GfunPar(1,i,b,spp,alp,bet,betp,epsp,epst,epss,gamv)
947c gamom=dble(alp*gamv)*chad(iclpro)*chad(icltar)
948c deltap=dble(bet)
949c deltam=dble(betp)
950c cint=gamom*utgam2(deltap+1.d0)*utgam2(deltam+1.d0)
951c & /utgam2(2.d0+deltap+dble(alplea(iclpro)))
952c & /utgam2(2.d0+deltam+dble(alplea(icltar)))
953c
954c om1intbci=cint*Fp*Fm
955c
956c if(om1intbci.lt.-1.d-10)then
957c write(*,*) 'WARNING ! om1intbci in epos-omg is <0 !!!!!'
958c write(*,*) 'WARNING ! => om1intbci set to 0. !!!!!'
959c om1intbci=0.d0
960c endif
961c
962c return
963c end
964c
965c----------------------------------------------------------------------
966 double precision function om1intgck(k,xprem,xmrem) !---MC---
967c----------------------------------------------------------------------
968c om1*(xprem-xp)*(xmrem-xm) integrated over xp and xm for given k
969c Calculation by analytical integration
970c----------------------------------------------------------------------
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
996c----------------------------------------------------------------------
997 double precision function om1intgc(b) !---test---
998c----------------------------------------------------------------------
999c om1*(1-xp)*(1-xm) integrated over xp and xm for given b
1000c Calculation by analytical integration
1001c----------------------------------------------------------------------
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
1048c----------------------------------------------------------------------
1049 subroutine integom1(irea)
1050c----------------------------------------------------------------------
1051c om1 integrated over xp and xm for all b=bk(k) if irea=0
1052c result written to :
1053c om1int(k) = om1intb(bk(k)) = \int om1
1054c om1intc(k) = om1intgc(bk(k)... = \int om1*(1-xp)*(1-xm)
1055c----------------------------------------------------------------------
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
1081c----------------------------------------------------------------------
1082 double precision function om1xpk(xp,xpr1i,xmr1i,k) !---test---
1083c----------------------------------------------------------------------
1084c \int dxm om1*(1-xp)*(1-xm) (normalized)
1085c k - pair indice;
1086c----------------------------------------------------------------------
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
1140c----------------------------------------------------------------------
1141 double precision function om1xmk(xp,xm,xpr1i,xmr1i,k) !---test---
1142c----------------------------------------------------------------------
1143c om1(xp,xm)*(1-xp)*(1-xm) (normalized for fixed xp)
1144c k - pair indice;
1145c----------------------------------------------------------------------
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
1189c----------------------------------------------------------------------
1190 double precision function om1xprk(k,xpremi,xmremi,ir) !---MC---
1191c----------------------------------------------------------------------
1192c Random number generated from the function om1xpk. We solve the equation
1193c which give om1xprk by Newton-Raphson + secant method.
1194c k - pair indice;
1195c ir - 1 to get xp, -1 to get xm (be carrefull to inverse xpremi et xmremi
1196c when calling with ir=-1)
1197c----------------------------------------------------------------------
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
1342c----------------------------------------------------------------------
1343 double precision function om1xmrk(k,xp,xpremi,xmremi,ir) !---MC---
1344c----------------------------------------------------------------------
1345c Random number generated from the function om1xmk. We solve the equation
1346c which give om1xmrk by Newton-Raphson + secant method.
1347c k - pair indice;
1348c ir - 1 to get xm, -1 to get xp (be carrefull to inverse xpremi et xmremi
1349c when calling with ir=-1)
1350c----------------------------------------------------------------------
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
1486c----------------------------------------------------------------------
1487 double precision function om1xk(xh,k) !---test---
1488c----------------------------------------------------------------------
1489c \int dxp om1 (normalised)
1490c k - pair indice;
1491c----------------------------------------------------------------------
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
1548c----------------------------------------------------------------------
1549 double precision function om1yk(xh,yp,k) !---test---
1550c----------------------------------------------------------------------
1551c om1 normalized for fixed xp
1552c xh - fraction of the energy squared s for the pomeron;
1553c k - pair indice;
1554c----------------------------------------------------------------------
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
1601c----------------------------------------------------------------------
1602 double precision function om1xrk(k) !---test---
1603c----------------------------------------------------------------------
1604c Random number generated from the function om1xk. We solve the equation
1605c which give om1xrk by Newton-Raphson + secant method.
1606c k - pair indice;
1607c----------------------------------------------------------------------
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
1703c if(f1*f0.gt.eps)then
1704c call utmsg('om1xrk')
1705c write(ifch,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
1706c write(ifmt,*)'Poblem with x0, no root ! --> om1xrk=xminDf'
1707c write(ifmt,*)f0,f1,cint,r
1708c call utmsgf
1709c om1xrk=x0
1710c return
1711c 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
1720c 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
1802c----------------------------------------------------------------------
1803 double precision function om1yrk(xh,k) !---test---
1804c----------------------------------------------------------------------
1805c Random number generated from the function om1yk(xh). We solve the
1806c equation which give om1yrk by Newton-Raphson + secant method.
1807c xh - fraction of the energy squared s for the pomeron;
1808c k - pair indice;
1809c----------------------------------------------------------------------
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
1822c----------------------------------------------------------------------
1823 function ffom12aii(iq,je1,je2) !---test---
1824c----------------------------------------------------------------------
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
1848c----------------------------------------------------------------------
1849 function ffom12ai(xp,iq1,iq2,je1,je2) !---test---
1850c----------------------------------------------------------------------
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
1866c----------------------------------------------------------------------
1867 function ffom12a(xp,xm,iq1,iq2,je1,je2) !---test---
1868c----------------------------------------------------------------------
1869c
1870c 2*om52*F*F == PomInc
1871c
1872c xp - xplus
1873c xm - xminus
1874c iq=1 .... sea-sea
1875c iq1 - min iq iq=2 .... val-sea
1876c iq2 - max iq iq=3 .... sea-val
1877c iq=4 .... val-val
1878c je = emission type (projectile and target side)
1879c 0 ... no emissions
1880c 1 ... emissions
1881c else ... all
1882c
1883c already b-averaged (\int d2b /sigine*10)
1884c----------------------------------------------------------------------
1885 include 'epos.inc'
1886
1887 sy=engy*engy
1888 xh=xm*xp
1889ctp060829 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
1908c----------------------------------------------------------------------
1909 function ffom11a(xp,xm,iq1,iq2) !---test---
1910c----------------------------------------------------------------------
1911c
1912c int(db) om1ff /sigine*10
1913c
1914c xp - xplus iq=-1 ... fit
1915c xm - xminus iq=0 .... soft
1916c iq=1 .... gg
1917c iq1 - min iq iq=2 .... qg
1918c iq2 - max iq iq=3 .... gq
1919c iq=4 .... qq
1920c iq=5 .... diff
1921c----------------------------------------------------------------------
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
1938c----------------------------------------------------------------------
1939 function ffom11(xp,xm,b,iq1,iq2) !---test---
1940c----------------------------------------------------------------------
1941c
1942c 2*om5*F*F == PomInc
1943c
1944c xp - xplus iq=-1 ... fit
1945c xm - xminus iq=0 .... soft
1946c b - impact parameter iq=1 .... gg
1947c iq1 - min iq iq=2 .... qg
1948c iq2 - max iq iq=3 .... gq
1949c iq=4 .... qq
1950c iq=5 .... diff
1951c----------------------------------------------------------------------
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
1984c----------------------------------------------------------------------
1985 double precision function om51(xh,yp,b,iq1,iq2) !---test---
1986c----------------------------------------------------------------------
1987c xh - xplus*xminus iq=-1 ... fit (om1 * 0.5)
1988c yp - rapidity iq=0 .... soft
1989c b - impact param iq=1 .... gg
1990c iq1 - min iq iq=2 .... qg
1991c iq2 - max iq iq=3 .... gq
1992c iq=4 .... qq
1993c iq=5 .... diff
1994c----------------------------------------------------------------------
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
2028c----------------------------------------------------------------------
2029 double precision function om5s(sx,xh,yp,b,iq1,iq2) !---test---
2030c----------------------------------------------------------------------
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
2041c----------------------------------------------------------------------
2042 double precision function om5Jk(k,xh,yp,iqq) !---MC---
2043c----------------------------------------------------------------------
2044c partial om5
2045c xh - fraction of the energy squared s for the pomeron;
2046c b - impact parameter between the pomeron ends;
2047c yp - rapidity for the pomeron;
2048c iqq=0 - soft
2049c iqq=1 - gg
2050c iqq=2 - qg
2051c iqq=3 - gq
2052c iqq=4 - qq
2053c iqq=5 - diffractif
2054c----------------------------------------------------------------------
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
2094c----------------------------------------------------------------------
2095 double precision function omIgamint(b,iqq) !---test---
2096c----------------------------------------------------------------------
2097c - integrated chi~(b)FF/2 for cut I diagram (simple Pomeron)
2098c b - impact parameter between the pomeron ends;
2099c yp - rapidity for the pomeron;
2100c iqq=0 effective one
2101c iqq=1 soft
2102c iqq=2 gg
2103c----------------------------------------------------------------------
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
2141c-----------------------------------------------------------------------
2142 subroutine WomTy(w,xh,yp,k)
2143c-----------------------------------------------------------------------
2144c - w(ity) for group iqq of cut enhanced diagram giving
2145c the probability of the type of the same final state.
2146c k - pair indice;
2147c xh - fraction of the energy squared s for the pomeron;
2148c yp - rapidity for the pomeron;
2149c xpr,xmr impulsion fraction of remnant
2150c ity = 0 - soft
2151c ity = 1 - gg
2152c ity = 2 - qg
2153c ity = 3 - gq
2154c ity = 4 - qq
2155c-----------------------------------------------------------------------
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
2198c-----------------------------------------------------------------------
2199 double precision function Womegak(xp,xm,xprem,xmrem,k) !---MC---
2200c-----------------------------------------------------------------------
2201c - sum(omGam(xp,xm))*(1-xp)*(1-xm) for group of cut enhanced
2202c diagram giving the same final state (without nuclear effect).
2203c xp,xm - fraction of the loght cone momenta of the pomeron;
2204c k - pair index
2205c-----------------------------------------------------------------------
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
2226cc----------------------------------------------------------------------
2227c double precision function omNpcut(xp,xm,xprem,xmrem,bh,iqq) !---test---
2228cc----------------------------------------------------------------------
2229cc Sum of all cut diagrams
2230cc iqq=0 ideal G
2231cc iqq=1 exact G + diff
2232cc----------------------------------------------------------------------
2233c include "epos.inc"
2234c double precision om51,xh,yp,xprem,xmrem,xp,xm!,omYcutI
2235c
2236c omNpcut=0.d0
2237c xh=xp*xm
2238c if(abs(xh).gt.1.d-10)then
2239c yp=0.5d0*log(xp/xm)
2240c else
2241c yp=0.d0
2242c endif
2243c
2244c if(iqq.eq.0)omNpcut=om51(xh,yp,bh,-1,-1)
2245c if(iqq.eq.1)omNpcut=om51(xh,yp,bh,0,5)
2246c
2247c omNpcut=omNpcut*2.d0
2248c
2249c return
2250c end
2251c
2252c----------------------------------------------------------------------
2253 double precision function omGam(xp,xm,bh) !---test---
2254c-----------------------------------------------------------------------
2255c Cut diagram part for calculation of probability distribution
2256c xp,xm impulsion fraction of remnant
2257c bh - impact parameter between the pomeron ends;
2258c-----------------------------------------------------------------------
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
2280c----------------------------------------------------------------------
2281 double precision function omGamk(k,xp,xm) !---MC---
2282c-----------------------------------------------------------------------
2283c Cut diagram part for calculation of probability distribution (for omega)
2284c xp,xm impulsion fraction of remnant
2285c bh - impact parameter between the pomeron ends;
2286c-----------------------------------------------------------------------
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
2300c----------------------------------------------------------------------
2301 double precision function omGamint(bh) !---test---
2302c-----------------------------------------------------------------------
2303c Integrated cut diagram part for calculation of probability distribution
2304c bh - impact parameter between the pomeron ends;
2305c-----------------------------------------------------------------------
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
2318c----------------------------------------------------------------------
2319 block data dgdata
2320c----------------------------------------------------------------------
2321c constants for numerical integration (gaussian weights)
2322c----------------------------------------------------------------------
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
2354c----------------------------------------------------------------------
2355 double precision function PhiExact(fj,xp,xm,s,b) !---test---
2356c----------------------------------------------------------------------
2357c Exact expression of the Phi function for pp collision
2358c----------------------------------------------------------------------
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
2446c----------------------------------------------------------------------
2447 double precision function PhiExpoK(k,xp,xm) !---MC---
2448c----------------------------------------------------------------------
2449c Exponential expression of the Phi function for pp collision
2450c for given k
2451c----------------------------------------------------------------------
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
2480c----------------------------------------------------------------------
2481 double precision function PhiExpo(fj,xp,xm,s,b) !---MC---
2482c----------------------------------------------------------------------
2483c Exponential expression of the Phi function for pp collision
2484c for given b
2485c----------------------------------------------------------------------
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
2518c----------------------------------------------------------------------
2519 double precision function PhiUnit(xp,xm) !---test---
2520c----------------------------------------------------------------------
2521c Exponential expression of the Phi function for pp collision
2522c for given b
2523c----------------------------------------------------------------------
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
2542c 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
2556cc----------------------------------------------------------------------
2557c double precision function PhiUnit(xp,xm,s,b) !---inu---
2558cc----------------------------------------------------------------------
2559c include 'epos.inc'
2560c double precision xp,xm,PhiExpo,Znorm
2561c
2562c PhiUnit=PhiExpo(1.,xp,xm,s,b)
2563c & /Znorm(s,b)
2564c
2565c return
2566c end
2567c
2568
2569c----------------------------------------------------------------------
2570 double precision function Hrst(s,b,xp,xm) !test
2571c----------------------------------------------------------------------
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.))
2605c 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)
2626c 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
2693c----------------------------------------------------------------------
2694 double precision function HrstI(s,b,xp,xm) !test
2695c----------------------------------------------------------------------
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.))
2731c 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
2819cc----------------------------------------------------------------------
2820c double precision function HrstI(s,xp,xm) !---inu---
2821cc----------------------------------------------------------------------
2822c include 'epos.inc'
2823c include 'epos.incems'
2824c include 'epos.incsem'
2825c include 'epos.incpar'
2826c double precision al(idxD0:idxD1),betp(idxD0:idxD1)
2827c *,z,xJrstI!,ffacto
2828c double precision zp(idxD0:idxD1),Htmp,betpp(idxD0:idxD1)
2829c *,yp,ym,xp,xm
2830c dimension ipr(idxD0:idxD1),imax(idxD0:idxD1)
2831c
2832c if(idxD0.ne.0.or.idxD1.ne.2) stop "Problem in HrstI"
2833c HrstI=0.d0
2834c do i=idxD0,idxD1
2835c imax(i)=0
2836c ipr(i)=0
2837c zp(i)=1.d0
2838c al(i)=0.d0
2839c enddo
2840c
2841c if(xp.ge.0.d0.and.xm.ge.0.d0.and.xp.lt.1.d0.and.xm.lt.1.d0)then
2842c
2843c HrstI=0.d0
2844c
2845c imax0=idxD1
2846c if(iomega.eq.2)imax0=1
2847c
2848c do i=idxDmin,imax0
2849c imax(i)=max(5,int(log10(s)))
2850cc if(i.ge.2)imax(i)=imax(i)*2
2851c imax(i)=min(30,imax(i))
2852c enddo
2853c Htmp=0.d0
2854c do i=idxDmin,imax0
2855c betp(i)=betUni(i,1)+1.d0
2856c betpp(i)=betpUni(i,1)+1.d0
2857c al(i)=alpUni(i,1)*dble(chad(iclpro)*chad(icltar))
2858c enddo
2859c
2860c do ipr0=0,imax(0)
2861c ipr(0)=ipr0
2862c zp(0)=1.d0
2863c if (ipr(0).ne.0) zp(0)=al(0)**ipr(0)*facto(ipr(0))
2864c do ipr1=0,imax(1)
2865c ipr(1)=ipr1
2866c zp(1)=1.d0
2867c if (ipr(1).ne.0) zp(1)=al(1)**ipr(1)*facto(ipr(1))
2868c do ipr2=0,imax(2)
2869c ipr(2)=ipr2
2870c zp(2)=1.d0
2871c if (ipr(2).ne.0) zp(2)=al(2)**ipr(2)*facto(ipr(2))
2872c if (ipr(0)+ipr(1)+ipr(2).ne.0) then
2873c yp=0.d0
2874c ym=0.d0
2875c z=1.d0
2876c do i=idxDmin,imax0
2877c yp=yp+dble(ipr(i))*betp(i)
2878c ym=ym+dble(ipr(i))*betpp(i)
2879c z=z*zp(i)
2880c enddo
2881c z=z*xJrstI(xp,yp,betp,ipr)
2882c z=z*xJrstI(xm,ym,betpp,ipr)
2883c Htmp=Htmp+z
2884c endif
2885c enddo
2886c enddo
2887c enddo
2888c
2889c HrstI=Htmp
2890c
2891c endif
2892c
2893c return
2894c end
2895c
2896
2897cc----------------------------------------------------------------------
2898c double precision function ffacto(n) !---test---
2899cc----------------------------------------------------------------------
2900c
2901c ffacto=1.D0
2902c do i=1,n
2903c ffacto=ffacto*dble(i)
2904c enddo
2905c return
2906c end
2907c
2908
2909c----------------------------------------------------------------------
2910 double precision function xIrst(id,x,y,bet,ipr) !---test---
2911c----------------------------------------------------------------------
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
2936c----------------------------------------------------------------------
2937 double precision function xJrst(x,y,Gbeta,ipr) !---inu---
2938c----------------------------------------------------------------------
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
2959c 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
2973c----------------------------------------------------------------------
2974 double precision function xJrstI(x,y,Gbeta,ipr) !---inu---
2975c----------------------------------------------------------------------
2976c Function used for the integration of H*Phi. We do the changement of
2977c variable (1-x)=z**alpha. The power alpha can be change if necessary.
2978c----------------------------------------------------------------------
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
3013c----------------------------------------------------------------------
3014 double precision function HPhiInt(s,b) !---inu---
3015c----------------------------------------------------------------------
3016c Set integrated over xp and xm (x and y) H(x,y)*Phi(x,y) for a
3017c given b by gauss method
3018c----------------------------------------------------------------------
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
3026c double precision zp2,zm2,HrstI,eps
3027c 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
3068c do i=1,7
3069 x=xhm+dble((2*m-3)*x9(i))*xhm
3070c write(ifmt,*)'HPhiInt, xp int :',x
3071 do n=1,2
3072 do j=1,3
3073c 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)
3077c & *PhiExact(1.,x,y,s,b)
3078 enddo
3079 enddo
3080 enddo
3081 enddo
3082
3083 HPhiInt=w*xhm*yhm
3084
3085
3086c w=0.d0
3087c xhm=.5d0*eps
3088c yhm=.5d0*eps
3089c do m=1,2
3090c do i=1,7
3091c x=1d0-eps+xhm+dble((2*m-3)*x1(i))*xhm
3092c do n=1,2
3093c do j=1,7
3094c y=1d0-epsyhm+dble((2*n-3)*x1(j))*yhm
3095c zp2=1.d0-x**4
3096c zm2=1.d0-y**4
3097c w=w+dble(a1(i)*a1(j))*HrstI(s,x,y)
3098cc & *PhiUnit(zp2,zm2)
3099c & *PhiExact(1.,zp2,zm2,s,b)
3100c enddo
3101c enddo
3102c enddo
3103c enddo
3104c
3105c HPhiInt=HPhiInt+w*xhm*yhm
3106
3107 return
3108 end
3109
3110
3111
3112c----------------------------------------------------------------------
3113 subroutine Kfit(iiclegy)
3114c----------------------------------------------------------------------
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
3266c----------------------------------------------------------------------
3267 double precision function Znorm(s,b) !---inu---
3268c----------------------------------------------------------------------
3269 include 'epos.inc'
3270 common /kwrite/ xkapZ
3271 double precision HPhiInt,PhiUnit!,PhiExact
3272
3273c 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)
3282c write(ifch,*)'int',Znorm,' phi',PhiExact(1.,1.d0,1.d0,s,b)
3283 Znorm=Znorm
3284 & +PhiUnit(1.d0,1.d0)
3285c & +PhiExact(1.,1.d0,1.d0,s,b)
3286
3287 !write(ifmt,*)'Z=',Znorm,xkapZ,b
3288 return
3289 end
3290
3291
3292c------------------------------------------------------------
3293 double precision function gammag(iclrem,x) !---test---
3294c----------------------------------------------------------------------
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
3307cc----------------------------------------------------------------------
3308c double precision function PomNbri(iqq) !---xsigma---
3309cc----------------------------------------------------------------------
3310cc integral d2b om1intbci
3311cc iqq, Pomeron type
3312cc----------------------------------------------------------------------
3313c include 'epos.inc'
3314c double precision om1intbci
3315c common/geom/rmproj,rmtarg,bmax,bkmx
3316c common /ar3/ x1(7),a1(7)
3317c
3318c bmid=bkmx/2.
3319c PomNbri=0.d0
3320c do i=1,7
3321c do m=1,2
3322c bb=bmid*(1.+(2.*m-3)*x1(i))
3323c PomNbri=PomNbri+dble(bb*a1(i))*om1intbci(bb,iqq)
3324c enddo
3325c enddo
3326c
3327c PomNbri=PomNbri*dble(2.*pi*bmid)
3328c
3329c return
3330c end
3331c
3332c
3333c
3334
3335
3336c####################################################################################
3337c############# former chk #########################################################
3338c####################################################################################
3339
3340
3341
3342
3343
3344
3345
3346
3347cc----------------------------------------------------------------------
3348c double precision function PomIncII(b) !---check---
3349cc----------------------------------------------------------------------
3350cc integral_dx_dy om1*F_remn*F_remn for a given b !---check---
3351cc----------------------------------------------------------------------
3352c include 'epos.inc'
3353c include 'epos.incems'
3354c include 'epos.incsem'
3355c include 'epos.incpar'
3356c double precision cint,gamom(idxD0:idxD1),deltap(idxD0:idxD1)
3357c &,deltapp(idxD0:idxD1),utgam2
3358c
3359cc Calculation by analytical integration (perfect but it changes
3360cc if om1 change):
3361c
3362c s=engy**2
3363c imax=1
3364c if(iomega.eq.2)imax=1
3365c do i=idxDmin,imax
3366c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3367c gamom(i)=dble(alp*gamv)*chad(iclpro)*chad(icltar)
3368c deltap(i)=dble(bet)
3369c deltapp(i)=dble(betp)
3370c
3371cc Integration possible only if delta(i)>-1
3372c
3373c if(deltap(i).le.-1.d0.or.deltapp(i).le.-1.d0)
3374c & stop 'Error in epos-par-300 in PomIncII'
3375c enddo
3376c
3377c cint=0.d0
3378c do i=idxDmin,imax
3379c cint=cint+gamom(i)*utgam2(deltap(i)+1.d0)*utgam2(deltapp(i)+1.d0)
3380c & *dble(ucfpro*ucftar)
3381c & /utgam2(dble(alplea(iclpro))+deltap(i)+2.d0)
3382c & /utgam2(dble(alplea(icltar))+deltapp(i)+2.d0)
3383c enddo
3384c
3385c PomIncII=cint
3386c
3387c return
3388c end
3389c
3390
3391c----------------------------------------------------------------------
3392 double precision function PomIncXIExact(x) !---check---
3393c----------------------------------------------------------------------
3394c integral d2b PomIncXExact
3395c----------------------------------------------------------------------
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
3415c----------------------------------------------------------------------
3416 double precision function PomIncXIUnit(x) !---check---
3417c----------------------------------------------------------------------
3418c integral d2b PomIncXUnit
3419c----------------------------------------------------------------------
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
3439c----------------------------------------------------------------------
3440 double precision function PomIncPIExact(x) !---check---
3441c----------------------------------------------------------------------
3442c integral d2b PomIncPExact
3443c----------------------------------------------------------------------
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
3463c----------------------------------------------------------------------
3464 double precision function PomIncPIUnit(x) !---check---
3465c----------------------------------------------------------------------
3466c integral d2b PomIncPUnit
3467c----------------------------------------------------------------------
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
3487c----------------------------------------------------------------------
3488 double precision function PomIncMIExact(x) !---check---
3489c----------------------------------------------------------------------
3490c integral d2b PomIncMExact
3491c----------------------------------------------------------------------
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
3511c----------------------------------------------------------------------
3512 double precision function PomIncMIUnit(x) !---check---
3513c----------------------------------------------------------------------
3514c integral d2b PomIncMUnit
3515c----------------------------------------------------------------------
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
3535c----------------------------------------------------------------------
3536 double precision function PomIncMExact(xm,b) !---check---
3537c----------------------------------------------------------------------
3538c incluse Pomeron distribution \int dx+ { 2G F_remn F_remn }
3539c----------------------------------------------------------------------
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
3566c----------------------------------------------------------------------
3567 double precision function PomIncMUnit(xm,b) !---check---
3568c----------------------------------------------------------------------
3569c incluse Unitarized Pomeron distribution \int dx+
3570c----------------------------------------------------------------------
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
3580c 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
3587ctp060829 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
3607c----------------------------------------------------------------------
3608 double precision function PomIncPExact(xp,b) !---check---
3609c----------------------------------------------------------------------
3610c incluse Pomeron distribution \int dx- { 2G F_remn F_remn }
3611c----------------------------------------------------------------------
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
3637c----------------------------------------------------------------------
3638 double precision function PomIncPUnit(xp,b) !---check---
3639c----------------------------------------------------------------------
3640c incluse Unitarized Pomeron distribution \int dx-
3641c----------------------------------------------------------------------
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
3652c 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
3659ctp060829 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
3680c----------------------------------------------------------------------
3681 double precision function PomIncJExact(b) !---check---
3682c----------------------------------------------------------------------
3683c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn }
3684c----------------------------------------------------------------------
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
3704c----------------------------------------------------------------------
3705 double precision function PomIncJUnit(b) !---check---
3706c----------------------------------------------------------------------
3707c integral of Pomeron distribution \int dy dx { 2G F_remn F_remn }
3708c----------------------------------------------------------------------
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
3728cc----------------------------------------------------------------------
3729c double precision function PomIncXExact1(xh,b) !---check---
3730cc----------------------------------------------------------------------
3731cc incluse Pomeron distribution \int dy { 2G F_remn F_remn }
3732cc----------------------------------------------------------------------
3733c include 'epos.inc'
3734c include 'epos.incsem'
3735c double precision xh,Df,xp,xm,w,ymax
3736c common /ar3/ x1(7),a1(7)
3737c
3738c imax=1
3739c if(iomega.eq.2)imax=1
3740c s=engy**2
3741ccctp060829 sy=s*real(xh)
3742cc Calculation by numeric integration :
3743c w=0.d0
3744c ymax=-.5d0*log(xh)
3745c do m=1,2
3746c do j=1,7
3747c xp=sqrt(xh)*exp(dble((2.*m-3.)*x1(j))*ymax)
3748c xm=xh/xp
3749c Df=0.d0
3750c do i=idxDmin,imax
3751c call GfunPar(1,i,b,s,alp,bet,betp,epsp,epst,epss,gamv)
3752c Df=Df+alp*gamv*real(xp)**bet*real(xm)**betp
3753c * *(1.d0-xp)**dble(alplea(iclpro))
3754c * *(1.d0-xm)**dble(alplea(icltar))
3755c enddo
3756c w=w+dble(a1(j))*Df
3757c enddo
3758c enddo
3759c w=w*ymax*xh**dble(-alppar)
3760c
3761c
3762c PomIncXExact1=w
3763c
3764c return
3765c end
3766
3767c----------------------------------------------------------------------
3768 double precision function PomIncXExact(xh,b) !---check---
3769c----------------------------------------------------------------------
3770c incluse Pomeron distribution \int dy { 2G F_remn F_remn }
3771c (optimized integration but with no y dependance)
3772c----------------------------------------------------------------------
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
3832c----------------------------------------------------------------------
3833 double precision function PomIncXUnit(xh,b) !---check---
3834c----------------------------------------------------------------------
3835c incluse Unitarized Pomeron distribution \int dy
3836c----------------------------------------------------------------------
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
3846ctp060829 sy=s*real(xh)
3847c 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
3874c----------------------------------------------------------------------
3875 double precision function PoInU(xp,xm,s,b) !---check---
3876c----------------------------------------------------------------------
3877c Function : PhiU(1-xp,1-xm) + /int(H(z+ + x+,z- + x-)PhiU(z+,z-)dz+dz-)
3878c----------------------------------------------------------------------
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)
3908c & *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)
3938c & *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)
3947c & +PhiExact(1.,1.d0-xp,1.d0-xm,s,b)
3948
3949 return
3950 end
3951
3952c----------------------------------------------------------------------
3953 double precision function PomIncExact(xp,xm,b) !---check---
3954c----------------------------------------------------------------------
3955c inclusive Pomeron distribution { 2G F_remn F_remn }
3956c----------------------------------------------------------------------
3957 include "epos.inc"
3958 include "epos.incsem"
3959 double precision Df,xp,xm!,xh
3960
3961 Df=0.d0
3962 xh=xp*xm
3963ctp060829 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
3981c----------------------------------------------------------------------
3982 double precision function PomIncUnit(xp,xm,b) !---check---
3983c----------------------------------------------------------------------
3984c inclusive Pomeron distribution { Sum{int{G*Phi} }
3985c----------------------------------------------------------------------
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
4002cc----------------------------------------------------------------------
4003c double precision function PomIncUnitMC(xp,xm,b) !---check---
4004cc----------------------------------------------------------------------
4005cc inclusive Pomeron distribution { Sum{int{G*Phi} }
4006cc----------------------------------------------------------------------
4007c include "epos.inc"
4008c include "epos.incpar"
4009c include "epos.incsem"
4010c include "epos.incems"
4011c parameter(mmax=20)
4012c double precision Gtab,Phitab,xxpu(mmax),xxmu(mmax)
4013c double precision Zm,xp,xm,pGtab,Z,omNpcut,xprem,xmrem,
4014c * sxp,sxm,PhiExpo
4015c
4016c PomIncUnitMC=0.d0
4017c if(xp.lt.1.d0.and.xm.lt.1.d0)then
4018c m=10
4019c
4020c sy=engy*engy
4021c nmcint=2000
4022c nmax=nmcint
4023c do i=1,mmax
4024c xxpu(i)=0.d0
4025c xxmu(i)=0.d0
4026c enddo
4027c xprem=1.d0
4028c xmrem=1.d0
4029c sxp=xprem-xp
4030c sxm=xmrem-xm
4031c
4032c Gtab=omNpcut(xp,xm,sxp,sxm,b,0)
4033c Phitab=PhiExpo(1.,sxp,sxm,sy,b)
4034c Z=Gtab*Phitab
4035c Zm=0.d0
4036c
4037c do mtmp=2,m
4038c
4039c write(*,*)"GPhi",mtmp-1,Zm,Z
4040c Zm=0.d0
4041c n=0
4042c
4043c 10 continue
4044c n=n+1
4045c if(mod(n,1000000).eq.0)write(*,*)
4046c & "Calculation of PomIncUnit(",mtmp,")->",n
4047c xxpu(1)=xp
4048c xxmu(1)=xm
4049c sxp=xxpu(1)
4050c sxm=xxmu(1)
4051c pGtab=1.d0
4052c do i=2,mtmp
4053c rnau=rangen()*real(xprem-sxp)
4054c xxpu(i)=dble(rnau)
4055c sxp=sxp+xxpu(i)
4056c rnau=rangen()*real(xmrem-sxm)
4057c xxmu(i)=dble(rnau)
4058c sxm=sxm+xxmu(i)
4059c enddo
4060c if(sxp.lt.xprem.and.sxm.lt.xmrem)then
4061c do i=1,mtmp
4062c Gtab=omNpcut(xxpu(i),xxmu(i),xprem-sxp,xmrem-sxm,b,0)
4063c pGtab=pGtab*Gtab
4064c enddo
4065c Zm=Zm+pGtab*PhiExpo(1.,xprem-sxp,xmrem-sxm,sy,b)
4066c endif
4067c if(n.lt.nmax)goto 10
4068c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(mtmp)
4069c Z=Z+Zm
4070c enddo
4071c
4072c PomIncUnitMC=Z/dble(chad(iclpro)*chad(icltar))
4073c endif
4074c
4075c return
4076c end
4077c
4078c
4079cc----------------------------------------------------------------------
4080c double precision function PhiMCExact(xp,xm,b) !---check---
4081cc----------------------------------------------------------------------
4082cc virtual emissions { Sum{int{-GFF} }
4083cc----------------------------------------------------------------------
4084c include "epos.inc"
4085c include "epos.incpar"
4086c include "epos.incsem"
4087c include "epos.incems"
4088c parameter(mmax=20)
4089c double precision Gtab,xxpu(mmax),xxmu(mmax)
4090c double precision Zm,xp,xm,pGtab,Z,om51,sxp,sxm,xh,yp
4091cc * ,omNpuncut
4092c
4093c PhiMCExact=0.d0
4094c if(xp.le.1.d0.and.xm.le.1.d0)then
4095c m=6
4096c
4097c sy=engy*engy
4098c nmcint=50000
4099c nmax=nmcint
4100c do i=1,mmax
4101c xxpu(i)=0.d0
4102c xxmu(i)=0.d0
4103c enddo
4104c
4105c Z=xp**dble(alplea(iclpro))
4106c * *xm**dble(alplea(icltar))
4107c Zm=0.d0
4108c
4109c do mtmp=1,m
4110c
4111c write(*,*)"GPhi",mtmp-1,Zm,Z/xp**dble(alplea(iclpro))
4112c * /xm**dble(alplea(icltar))
4113c Zm=0.d0
4114c n=0
4115c
4116c 10 continue
4117c n=n+1
4118c if(mod(n,1000000).eq.0)write(*,*)
4119c & "Calculation of PhiExact(",mtmp,")->",n
4120c sxp=0.d0
4121c sxm=0.d0
4122c pGtab=1.d0
4123c do i=1,mtmp
4124c rnau=rangen()!*real(xp-sxp)
4125c xxpu(i)=dble(rnau)
4126c sxp=sxp+xxpu(i)
4127c rnau=rangen()!*real(xm-sxm)
4128c xxmu(i)=dble(rnau)
4129c sxm=sxm+xxmu(i)
4130c enddo
4131c if(sxp.lt.xp.and.sxm.lt.xm)then
4132c do i=1,mtmp
4133c xh=xxpu(i)*xxmu(i)
4134c if(abs(xh).gt.1.d-30)then
4135c yp=0.5d0*log(xxpu(i)/xxmu(i))
4136c else
4137c yp=0.d0
4138c endif
4139c Gtab=2*om51(xh,yp,b,0,4)
4140cc * +omNpuncut(sy*real(xh),xh,yp,b,1) !om1(xh,yp,b)
4141c pGtab=pGtab*(-Gtab)
4142c enddo
4143c Zm=Zm+pGtab*(xp-sxp)**dble(alplea(iclpro))
4144c * *(xm-sxm)**dble(alplea(icltar))
4145c endif
4146c if(n.lt.nmax)goto 10
4147c Zm=Zm/dble(nmax)*fctrl(m-mtmp)*facto(m)
4148c Z=Z+Zm
4149c enddo
4150c
4151c PhiMCExact=Z
4152c endif
4153c
4154c return
4155c end
4156
4157c----------------------------------------------------------------------
4158 double precision function Gammapp(sy,b,mtmp) !---check---
4159c----------------------------------------------------------------------
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
4218cc----------------------------------------------------------------------
4219c double precision function GammaMCnew(sy,b,mtmp) !---check---
4220cc----------------------------------------------------------------------
4221c include "epos.inc"
4222c include "epos.incpar"
4223c include "epos.incsem"
4224c include "epos.incems"
4225c parameter(mmax=20)
4226c common /psar7/ delx,alam3p,gam3p
4227c double precision Gtab,xxpu(mmax),xxmu(mmax)
4228c double precision Zm,xp,xm,pGtab,omGam,Zmtot,
4229c * sxp,sxm,PhiExpo!,om1,yp,xh
4230c
4231c GammaMCnew=0.d0
4232c Zmtot=0.d0
4233c xp=1.d0
4234c xm=1.d0
4235c nmcint=1000
4236c nmax=nmcint
4237c
4238c
4239c do i=1,mmax
4240c xxpu(i)=0.d0
4241c xxmu(i)=0.d0
4242c enddo
4243c
4244c Zm=0.d0
4245c
4246c n=0
4247c
4248c 10 continue
4249c n=n+1
4250c if(mod(n,1000000).eq.0)write(*,*)
4251c & "Calculation of GammaMCnew(",mtmp,")->",n
4252c sxp=0.d0
4253c sxm=0.d0
4254c pGtab=1.d0
4255c do i=1,mtmp
4256c rnau=rangen()
4257c xxpu(i)=dble(rnau)
4258c sxp=sxp+xxpu(i)
4259c rnau=rangen()
4260c xxmu(i)=dble(rnau)
4261c sxm=sxm+xxmu(i)
4262c enddo
4263c if(sxp.lt.xp.and.sxm.lt.xm)then
4264c i=0
4265c do k=1,mtmp
4266c i=i+1
4267c Gtab=omGam(xxpu(i),xxmu(i),b) !om1(xh,yp,b)
4268c pGtab=pGtab*Gtab
4269c enddo
4270c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4271c endif
4272c if(n.lt.nmax)goto 10
4273c
4274c Zmtot=Zmtot+Zm/dble(nmax)
4275c
4276c GammaMCnew=Zmtot
4277c
4278c return
4279c end
4280cc----------------------------------------------------------------------
4281c double precision function GammaMC(sy,b,mtmp) !---check---
4282cc----------------------------------------------------------------------
4283c include "epos.inc"
4284c include "epos.incpar"
4285c include "epos.incsem"
4286c include "epos.incems"
4287c parameter(mmax=20)
4288c double precision Gtab,xxpu(mmax),xxmu(mmax)
4289c double precision Zm,xp,xm,pGtab,om1,
4290c * sxp,sxm,xh,yp,PhiExpo!,omNpcut
4291c
4292c GammaMC=0.d0
4293c
4294c xp=1.d0
4295c xm=1.d0
4296c nmcint=50000
4297c nmax=nmcint
4298c do i=1,mmax
4299c xxpu(i)=0.d0
4300c xxmu(i)=0.d0
4301c enddo
4302c
4303c Zm=0.d0
4304c
4305c n=0
4306c
4307c 10 continue
4308c n=n+1
4309c if(mod(n,1000000).eq.0)write(*,*)
4310c & "Calculation of GammaMC(",mtmp,")->",n
4311c sxp=0.d0
4312c sxm=0.d0
4313c pGtab=1.d0
4314c do i=1,mtmp
4315c rnau=rangen()!*real(xp-sxp)
4316c xxpu(i)=dble(rnau)
4317c sxp=sxp+xxpu(i)
4318c rnau=rangen()!*real(xm-sxm)
4319c xxmu(i)=dble(rnau)
4320c sxm=sxm+xxmu(i)
4321c enddo
4322c if(sxp.lt.xp.and.sxm.lt.xm)then
4323c do i=1,mtmp
4324c xh=xxpu(i)*xxmu(i)
4325c if(abs(xh).gt.1.d-30)then
4326c yp=0.5d0*log(xxpu(i)/xxmu(i))
4327c else
4328c yp=0.d0
4329c endif
4330c Gtab=om1(xh,yp,b)!omNpcut(xxpu(i),xxmu(i),xp-sxp,xm-sxm,b,0) !om1(xh,yp,b)
4331c pGtab=pGtab*Gtab
4332c enddo
4333c Zm=Zm+pGtab*PhiExpo(1.,xp-sxp,xm-sxm,sy,b)
4334c endif
4335c if(n.lt.nmax)goto 10
4336c Zm=Zm/dble(nmax)*fctrl(n-mtmp)*facto(mtmp)
4337c
4338c GammaMC=Zm
4339c
4340c return
4341c end
4342c
4343c
4344c
4345
4346c----------------------------------------------------------------------
4347 double precision function GammaGauss(sy,b,mtmp) !---check---
4348c----------------------------------------------------------------------
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
4359c *,PhiExact
4360 common /ar3/ x1(7),a1(7)
4361c double precision dgx1,dga1
4362c 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
4427c 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
4452c write(*,*)'bornes1=',exp(-zpmax(1)),exp(-zpmin(1))
4453c &,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))
4499c &dble(a1(jm1))*PhiExact(1.,sxp,sxm,sy,b)
4500c &*exp(-zm(1))
4501 else
4502 Gp1=Gp1+pGtab*
4503 &dble(a1(jp1))*PhiExpo(1.,sxp,sxm,sy,b)
4504 &*exp(-zp(1))
4505c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)
4506c &*exp(-zp(1))
4507 endif
4508c 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
4535c write(*,*)'bornes2=',exp(-zpmax(2)),exp(-zpmin(2))
4536c &,exp(-zmmax(2)),exp(-zmmin(2)),xpmax(2),1.d0-exp(-zp(1))
4537c &,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))
4583c &dble(a1(jm2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4584c &*exp(-zm(2))
4585 else
4586 Gp2=Gp2+pGtab*
4587 &dble(a1(jp2))*PhiExpo(1.,sxp,sxm,sy,b)
4588 &*exp(-zp(2))
4589c &dble(a1(jp2))*PhiExact(1.,sxp,sxm,sy,b,mk)
4590c &*exp(-zp(2))
4591 endif
4592c 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
4610c write(*,*)'bornes3=',exp(-zpmax(3)),exp(-zpmin(3))
4611c &,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
4703c-----------------------------------------------------------------------
4704 double precision function omWi(sy,b) !---check---
4705c-----------------------------------------------------------------------
4706c cut enhanced diagram integrated over xp, xm, xpr,xmr
4707c (with ideal G)
4708c b - impact parameter between the pomeron ends;
4709c sy- total energy
4710c-----------------------------------------------------------------------
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
4720c *,PhiExact
4721 common /ar3/ x1(7),a1(7)
4722c double precision dgx1,dga1
4723c 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)
4732ctp060829 gamp=dble(gamb*b2/2.-alppar)
4733ctp060829 gamm=dble(gamb*b2/2.-alppar)
4734ctp060829 gampp=1.d0+2.d0*gamp
4735ctp060829 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
4785c 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
4816c 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)))
4821c 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)))
4834c 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
4858c &dga1(jm1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacm
4859c &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
4863c &dga1(jp1)*PhiExpo(1.,sxp,sxm,sy,b)*xjacp
4864c &dble(a1(jp1))*PhiExact(1.,sxp,sxm,sy,b)*xjacp
4865 endif
4866c 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
4872c 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
4886c-----------------------------------------------------------------------
4887 double precision function Womint(sy,bh) !---check---
4888c-----------------------------------------------------------------------
4889c - chi~(xp,xm)/2. for group of cut enhanced diagram giving
4890c the same final state integrated over xpr and xmr (with ideal G)
4891c bh - impact parameter between the pomeron ends;
4892c xh - fraction of the energy squared s for the pomeron;
4893c yp - rapidity for the pomeron;
4894c-----------------------------------------------------------------------
4895 include 'epos.inc'
4896 double precision omWi
4897
4898 Womint=omWi(sy,bh)
4899
4900 return
4901 end
4902
4903
4904
4905c-----------------------------------------------------------------------
4906 double precision function WomGamint(bh) !---check---
4907c-----------------------------------------------------------------------
4908c - chi~(xp,xm)/2. for group of integrated cut enhanced diagram giving
4909c the same final for proposal.
4910c bh - impact parameter between the pomeron ends;
4911c xh - fraction of the energy squared s for the pomeron;
4912c yp - rapidity for the pomeron;
4913c-----------------------------------------------------------------------
4914 include 'epos.inc'
4915 double precision omGamint
4916
4917 WomGamint=omGamint(bh)
4918
4919 return
4920 end
4921