]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-qsh-165.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
ALIROOT-5836 AliESDpid not respecting the AliVTrack interface (patch from Mihaela)
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-qsh-165.f
CommitLineData
9ef1c2d9 1c reshuffled from sem, sto, sha
2
3c contains DIS, and unused 3P stuff
4c --- ---------------
5
6
7
8
9c###########################################################################
10c###########################################################################
11c###########################################################################
12c###########################################################################
13c
14c DIS
15c
16c###########################################################################
17c###########################################################################
18c###########################################################################
19c###########################################################################
20
21
22
23
24
25
26
27
28
29
30c-----------------------------------------------------------------------
31 subroutine lepexp(rxbj,rqsq)
32c-----------------------------------------------------------------------
33c generates x_bjorken and q**2 according to an experimental
34c distribution ( given in array xq(nxbj,nqsq) ).
35c-----------------------------------------------------------------------
36 parameter (nxbj=10,nqsq=10)
37 parameter (xbjmin=0.,qsqmin=4.)
38 parameter (xbjwid=0.025, qsqwid=4.)
39 dimension xq(nxbj,nqsq),vxq(nxbj*nqsq)
40 equivalence (xq(1,1),vxq(1))
41
42 data (vxq(i),i=1,50)/
43 & 1304.02, 366.40, 19.84, 10.79, 6.42,
44 & 4.54, 4.15, 3.38, 2.03, 1.56,
45 & 241.63, 1637.26, 427.36, 164.51, 73.72,
46 & 43.07, 20.73, 12.78, 9.34, 5.83,
47 & 0.01, 724.66, 563.79, 275.08, 176.13,
48 & 106.44, 85.82, 54.52, 37.12, 28.65,
49 & 0.01, 202.40, 491.10, 245.13, 157.07,
50 & 104.43, 61.05, 49.42, 37.84, 26.79,
51 & 0.01, 3.77, 316.38, 226.92, 133.45,
52 & 90.30, 63.67, 48.42, 35.73, 28.04/
53 data (vxq(i),i=51,100)/
54 & 0.01, 0.01, 153.74, 213.09, 114.14,
55 & 76.26, 60.02, 43.15, 43.47, 25.60,
56 & 0.01, 0.01, 39.31, 185.74, 108.56,
57 & 88.40, 47.29, 39.35, 31.80, 22.91,
58 & 0.01, 0.01, 0.01, 104.61, 107.01,
59 & 66.24, 45.34, 37.45, 33.44, 23.78,
60 & 0.01, 0.01, 0.01, 56.58, 99.39,
61 & 67.78, 43.28, 35.98, 34.63, 18.31,
62 & 0.01, 0.01, 0.01, 13.56, 76.25,
63 & 64.30, 42.80, 28.56, 21.19, 20.75 /
64
65 data init/0/
66 init=init+1
67 if(init.eq.1) then
68 n=nxbj*nqsq
69 sum=0.
70 do 1 i=1,n
71 sum=sum+vxq(i)
721 continue
73 do 2 i=2,n
742 vxq(i)=vxq(i)+vxq(i-1)
75 do 3 i=1,n
763 vxq(i)=vxq(i)/sum
77 endif
78
79 n=nxbj*nqsq
80 r=rangen()
81 call utloc(vxq,n,r,iloc)
82 if(iloc.ge.n) iloc=iloc-1
83 i=mod(iloc,nxbj)+1
84 if(i.eq.0) i=nxbj
85 j=iloc/nxbj + 1
86 dxint=vxq(1)
87 if(iloc.gt.0) dxint=vxq(iloc+1)-vxq(iloc)
88 dxbj=xbjwid*abs(r-vxq(iloc+1))/dxint
89 dy=qsqwid*rangen()
90 rxbj=xbjmin+xbjwid*float(i-1)+dxbj
91 rqsq=qsqmin+qsqwid*float(j-1)+dy
92 return
93 end
94
95c-----------------------------------------------------------------------
96 subroutine fremny(wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0)
97c-----------------------------------------------------------------------
98c treats remnant from deep inelastic process;
99c-----------------------------------------------------------------------
100 include 'epos.inc'
101 include 'epos.incsem'
102 dimension coord(6),ic(2),ep(4),ey(3),ey0(3),ep3(4)
103 double precision ept(4),ept1(4)
104
105 call utpri('fremny',ish,ishini,5)
106 if(ish.ge.5)write(ifch,*)'writing remnant'
107 if(ish.ge.5)write(ifch,*)
108 *'wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0:'
109 if(ish.ge.5)write(ifch,*)
110 *wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0
111
112 if(ic3.eq.0.and.ic4.eq.0)then
113
114 ic(1)=ic1
115 ic(2)=ic2
116 nptl=nptl+1
117 ep3(3)=pnx
118 ep3(4)=pny
119 ep3(2)=(wp1-wm1)/2
120 ep3(1)=(wp1+wm1)/2
121 call pstrans(ep3,ey0,-1)
122 pptl(1,nptl)=ep3(3)
123 pptl(2,nptl)=ep3(4)
124 pptl(3,nptl)=ep3(2)
125 pptl(4,nptl)=ep3(1)
126 pptl(5,nptl)=sqrt(sm)
127 idptl(nptl)=idtra(ic,0,0,3)
128 iorptl(nptl)=1
129 istptl(nptl)=0
130 jorptl(nptl)=0
131 do i=1,4
132 xorptl(i,nptl)=coord(i)
133 enddo
134 tivptl(1,nptl)=coord(5)
135 tivptl(2,nptl)=coord(6)
136 ityptl(nptl)=40
137
138 if(ish.ge.7)then
139 write(ifch,*)'proj: nptl, mass**2,id',nptl,sm,idptl(nptl)
140 write(ifch,*)'ept',
141 *pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl)
142 endif
143
144 else
145
146 ic(1)=ic1
147 ic(2)=ic2
148 nptl=nptl+1
149 idptl(nptl)=idtra(ic,0,0,3)
150 istptl(nptl)=20
151 iorptl(nptl)=1
152 jorptl(nptl)=0
153 do i=1,4
154 xorptl(i,nptl)=coord(i)
155 enddo
156 tivptl(1,nptl)=coord(5)
157 tivptl(2,nptl)=coord(6)
158 ityptl(nptl)=40
159
160 ic(1)=ic3
161 ic(2)=ic4
162 idptl(nptl+1)=idtra(ic,0,0,3)
163 istptl(nptl+1)=20
164 iorptl(nptl+1)=1
165 jorptl(nptl+1)=0
166 do i=1,4
167 xorptl(i,nptl+1)=coord(i)
168 enddo
169 tivptl(1,nptl+1)=coord(5)
170 tivptl(2,nptl+1)=coord(6)
171 ityptl(nptl+1)=40
172
173 ep3(3)=pnx
174 ep3(4)=pny
175 ep3(2)=(wp1-wm1)/2
176 ep3(1)=(wp1+wm1)/2
177 call pstrans(ep3,ey0,-1) !boost to hadronic c.m.s.
178 ept(1)=ep3(3)
179 ept(2)=ep3(4)
180 ept(3)=ep3(2)
181 ept(4)=ep3(1)
182 do i=1,4
183 ept1(i)=ep3(i)
184 enddo
185
186 sww=sqrt(sm)
187 call psdeftr(sm,ept1,ey)
188 ep(1)=.5*sww
189 ep(2)=.5*sww
190 ep(3)=0.
191 ep(4)=0.
192 call pstrans(ep,ey,1)
193 pptl(1,nptl)=ep(3)
194 pptl(2,nptl)=ep(4)
195 pptl(3,nptl)=ep(2)
196 pptl(4,nptl)=ep(1)
197 do i=1,4
198 pptl(i,nptl+1)=ept(i)-pptl(i,nptl)
199 enddo
200 nptl=nptl+1
201 endif
202
203 1000 continue
204 if(ish.ge.5)write(ifch,*)'fremny: final nptl',nptl
205 call utprix('fremny',ish,ishini,5)
206 return
207 end
208
209c-----------------------------------------------------------------------
210 subroutine psadis(iret)
211c-----------------------------------------------------------------------
212c psadis - DIS interaction
213c-----------------------------------------------------------------------
214 double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,psutz
215 *,psuds
216 dimension ep3(4),ey(3),ey0(3),bx(6),
217 *qmin(2),iqc(2),nqc(2),ncc(2,2),gdv(2),gds(2),dfp(4)
218 parameter (mjstr=20000)
219 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
220 common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
221 double precision pgampr,rgampr
222 common/cgampr/pgampr(5),rgampr(4)
223 parameter (ntim=1000)
224 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
225 &,idaprt(2,ntim)
226 common/ciptl/iptl
227 include 'epos.inc'
228 include 'epos.incsem'
229
230 call utpri('psadis',ish,ishini,3)
231 if(ish.ge.3)write (ifch,*)'engy,elepti,iolept:'
232 if(ish.ge.3)write (ifch,*)engy,elepti,iolept
233 nptl=nptl+1
234 idptl(nptl)=1220
235 istptl(nptl)=1
236 nptlh=nptl
237 iptl=nptl
238 s00=1.
239
240 pptl(1,nptl)=0.
241 pptl(2,nptl)=0.
242 pptl(3,nptl)=-engy/2
243 pptl(4,nptl)=engy/2
244 pptl(5,nptl)=0
245
2461 continue
247 if(iolept.eq.1)then
248 wtot=engy**2
249 engypr=wtot/4./elepti
250 gdv01=psdh(ydmax*wtot,qdmin,iclpro,0)
251 gdv02=psdh(ydmax*wtot,qdmin,iclpro,1)
252 gds01=psdsh(ydmax*wtot,qdmin,iclpro,dqsh,0)
253 gds02=psdsh(ydmax*wtot,qdmin,iclpro,dqsh1,1)
254 gb0=(1.+(1.-ydmax)**2)*(gdv01+gds01)
255 * +2.*(1.-ydmax)*(gdv02+gds02)
256
2572 continue
258 qq=qdmin*(qdmax/qdmin)**rangen()
259 yd=ydmin*(ydmax/ydmin)**rangen()
260 wd=yd*wtot
261 if(ish.ge.4)write (ifch,*)'qq,wd,yd,ydmin,ydmax:'
262 if(ish.ge.4)write (ifch,*)qq,wd,yd,ydmin,ydmax
263 if(wd.lt.qq)goto 2
264
265 gdv(1)=psdh(wd,qq,iclpro,0)
266 gdv(2)=psdh(wd,qq,iclpro,1)
267 gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
268 gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
269 gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
270 gblong=2.*(1.-yd)*(gdv(2)+gds(2))
271c gblong=0. !???????
272 gb=(gbtr+gblong)/gb0*.7
273 if(ish.ge.4)then
274 if(gb.gt.1.)write(ifmt,*)'gb,qq,yd,wd',gb,qq,yd,wd
275 write (ifch,*)'gb,gdv,gds,gdv0,gds0,yd:'
276 write (ifch,*)gb,gdv,gds,gdv01,gds01,
277 * gdv02,gds02,yd
278 endif
279 if(rangen().gt.gb)goto 2
280
281 long=int(rangen()+gblong/(gbtr+gblong))
282 elepto=qq/elepti/4.+elepti*(1.-yd)
283 costhet=1.-qq/elepti/elepto/2.
284 theta=acos(costhet)
285 if(theta/pi*180..lt.themin)goto 2
286 if(theta/pi*180..gt.themax)goto 2
287 if(elepto.lt.elomin)goto 2
288 if(ish.ge.3)write (ifch,*)'theta,elepto,elepti,iclpro:'
289 if(ish.ge.3)write (ifch,*)theta/pi*180.,elepto,elepti,iclpro
290 xbjevt=qq/wd
291 qsqevt=qq
292
293 call pscs(bcos,bsin)
294 rgampr(1)=-elepto*sin(theta)*bcos
295 rgampr(2)=-elepto*sin(theta)*bsin
296 rgampr(3)=elepti-elepto*costhet
297 rgampr(4)=elepti-elepto
298
299 pgampr(1)=rgampr(1)
300 pgampr(2)=rgampr(2)
301 pgampr(3)=rgampr(3)-engypr
302 pgampr(4)=rgampr(4)+engypr
303 sm2=pgampr(4)*pgampr(4)
304 * -pgampr(1)*pgampr(1)-pgampr(2)*pgampr(2)-pgampr(3)*pgampr(3)
305 pgampr(5)=sqrt(sm2)
306 call utlob2(1,pgampr(1),pgampr(2),pgampr(3),pgampr(4),pgampr(5)
307 * ,rgampr(1),rgampr(2),rgampr(3),rgampr(4),40)
308 if(ish.ge.4)write (ifch,*)'rgampr:',rgampr
309
310 elseif(iolept.lt.0)then
31121 call lepexp(xbjevt,qsq)
312 qq=qsq
313 wd=qq/xbjevt
314 if(qq.lt.qdmin.or.qq.gt.qdmax)goto21
315
316 gdv(1)=psdh(wd,qq,iclpro,0)
317 gdv(2)=psdh(wd,qq,iclpro,1)
318 gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
319 gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
320 yd=wd/engy**2
321 gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
322 gblong=2.*(1.-yd)*(gdv(2)+gds(2))
323 gblong=0. !????????????
324 long=int(rangen()+gblong/(gbtr+gblong))
325 else
326 stop'wrong iolept'
327 endif
328 if(ish.ge.3)write (ifch,*)'qq,xbj,wd,gdv,gds,dqsh:'
329 if(ish.ge.3)write (ifch,*)qq,xbjevt,wd,gdv,gds,dqsh
330
331 egyevt=sqrt(wd-qq)
332 pmxevt=.5*egyevt
333
334 wp0=sqrt(qq) !breit frame
335 wm0=(wd-qq)/wp0
336 ey0(1)=egyevt/wp0 !boost to the hadronic c.m.s.
337 ey0(2)=1.
338 ey0(3)=1.
339 do i=1,6
340 bx(i)=0.
341 enddo
342
343 if(long.eq.0)then
344 sdmin=qq/(1.-sqrt(q2ini/qq))
345 sqmin=sdmin
346 else
347 sdmin=4.*max(q2min,qcmass**2)+qq !minimal mass for born
348 xmm=(5.*sdmin-qq)/4.
349 sqmin=1.1*(xmm+sqrt(xmm**2-qq*(sdmin-qq-4.*q2ini)))
350 * /2./(1.-4.*q2ini/(sdmin-qq))
351 endif
352 if(long.eq.1.and.wd.lt.1.001*sdmin)goto 1
353
354 proja=210000.
355 projb=0.
356 call fremnu(ammin,proja,projb,proja,projb,
357 *icp1,icp2,icp3,icp4)
358
359 nj=0
360
361 if((rangen().lt.gdv(long+1)/(gdv(long+1)+gds(long+1)).or.
362 *egyevt.lt.1.0001*(ammin+sqrt(sdmin-qq))).and.
363 *(long.eq.0.or.wd.gt.sqmin))then
364 if(long.eq.0)then
365 xd=qq/wd
366 tu=psdfh4(xd,q2min,0.,iclpro,1)/2.25
367 td=psdfh4(xd,q2min,0.,iclpro,2)/9.
368 gdv0=(tu+td)*4.*pi**2*alfe/qq
369 * *sngl(psuds(qq,1)/psuds(q2min,1))
370 if(ish.ge.4)write (ifch,*)'gdv0:',gdv0,sdmin
371
372 if(rangen().lt.gdv0/gdv(1).or.wd.le.1.0001*sdmin)then !?????
373 if(ish.ge.3)write (ifch,*)'no cascade,gdv0,gdv',gdv0,gdv
374 if(rangen().lt.tu/(tu+td))then
375 iq1=1
376 izh=3
377 else
378 iq1=2
379 izh=6
380 endif
381 jq=1
382 if(ish.ge.8)write (ifch,*)'before call timsh2: ',
383 * 'qq,egyevt,iq1',qq,egyevt,iq1
384 call timsh2(qq,0.,egyevt,iq1,-iq1)
385
386 nj=nj+1
387 iqj(nj)=izh
388 nqc(1)=nj
389 nqc(2)=0
390
391 ep3(1)=pprt(4,2)
392 ep3(2)=pprt(3,2)
393 ep3(3)=0.
394 ep3(4)=0.
395 call pstrans(ep3,ey0,1)
396 do i=1,4
397 eqj(i,nj)=ep3(i)
398 enddo
399 s0h=0.
400 c0h=1.
401 s0xh=0.
402 c0xh=1.
403 call psreti(nqc,jq,1,ey0,s0xh,c0xh,s0h,c0h)
404 goto 17
405 endif
406 endif
407
408 call psdint(wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,1,long)
409 if(ish.ge.3)write (ifch,*)'wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
410 if(ish.ge.3)write (ifch,*)wd,qq,sds0,sdn0,sdb0,sdt0
411 gb10=(sdn0+sdt0)*(1.-qq/wd)
412 xdmin=sqmin/wd
413
4143 continue
415 xd=(xdmin-qq/wd)/((xdmin-qq/wd)/(1.-qq/wd))
416 * **rangen()+qq/wd
417 call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,1,long)
418 if(ish.ge.3)write (ifch,*)'wdhard,qq,sds,sdn,sdb,sdt:'
419 if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt
420 tu=psdfh4(xd,q2min,0.,iclpro,1)
421 td=psdfh4(xd,q2min,0.,iclpro,2)
422 gb1=(sdn*(tu/2.25+td/9.)+sdt*(tu+td)/4.5)
423 * *(1.-qq/wd/xd)/gb10
424 if(gb1.gt.1..and.ish.ge.1)write(ifmt,*)'gb1,xd,wd,qq,sdt0,sdt',
425 * gb1,xd,wd,qq,sdt0,sdt
426 if(ish.ge.6)write (ifch,*)'gb1,xd,wd,qq,sdt0,sdt:'
427 if(ish.ge.6)write (ifch,*)gb1,xd,wd,qq,sdt0,sdt
428 if(rangen().gt.gb1)goto 3
429
430 gdres=(sdt-sds)/4.5
431 gdrga=sdr/4.5
432 gdsin=sds/4.5
433 dtu=tu*(sdn/2.25+sdt/4.5)
434 dtd=td*(sdn/9.+sdt/4.5)
435 if(rangen().lt.dtu/(dtu+dtd))then
436 iq1=1
437 izh=3
438 gdbor=sdb/2.25
439 gdnon=sdn/2.25
440 else
441 iq1=2
442 izh=6
443 gdbor=sdb/9.
444 gdnon=sdn/9.
445 endif
446
447 wpi=wp0
448 wmi=(xd*wd-qq)/wpi
449 iqc(2)=iq1
450 nj=nj+1
451 iqj(nj)=izh
452 eqj(1,nj)=.5*(wm0-wmi)
453 eqj(2,nj)=-eqj(1,nj)
454 eqj(3,nj)=0.
455 eqj(4,nj)=0.
456 ncc(1,2)=nj
457 ncc(2,2)=0
458 if(ish.ge.3)write (ifch,*)'wp0,wm0,wpi,wmi,iqc(2),eqj'
459 if(ish.ge.3)write (ifch,*)wp0,wm0,wpi,wmi,iqc(2),eqj(2,nj)
460
461 else
462 xdmin=sdmin/wd
463 xpmax=((egyevt-ammin)**2+qq)/wd
464 iq1=int(3.*rangen()+1.)*(2.*int(.5+rangen())-1.)
465
466 aks=rangen()
467 if(long.eq.0.and.aks.lt.dqsh/gds(1).and.
468 * egyevt.gt.ammin+sqrt(s00))then
469 if(ish.ge.3)write (ifch,*)'no cascade for q_s',
470 * aks,dqsh/gds(1)
471 xd=qq/wd
472 xpmin=xd+s00/wd
473 jcasc=0
474 if(iq1.gt.0)then
475 jq=1
476 else
477 jq=2
478 endif
479 else
480 jcasc=1
481 call psdint(xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,0,long)
482 call psdint(xpmax*wd,qq,sdsq0,sdnq0,sdbq0,sdtq0,sdrq0,1,long)
483 if(ish.ge.3)write (ifch,*)
484 * 'xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
485 if(ish.ge.3)write (ifch,*)
486 * xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0
487 gb10=sdt0*fzeroGluZZ(0.,iclpro)+(sdnq0+sdtq0)
488 * *fzeroSeaZZ(0.,iclpro)
489 gb10=gb10*15.
490
4914 xd=xdmin*(xpmax/xdmin)**rangen()
492 xpmin=xd
493 call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,0,long)
494 call psdint(xd*wd,qq,sdsq,sdnq,sdbq,sdtq,sdrq,1,long)
495 if(ish.ge.3)write (ifch,*)'xd*wd,qq,sds,sdn,sdb,sdt,sdr:'
496 if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt,sdr
497 wwg=sdt*fzeroGluZZ(xd,iclpro)
498 wwq=(sdnq+sdtq)*fzeroSeaZZ(xd,iclpro)
499 gb12=(wwq+wwg)/gb10*(xpmax/xd)**dels
500 if(gb12.gt.1..and.ish.ge.1)write(ifmt,*)
501 * 'gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq',
502 * gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq,
503 * wwq,wwg,(xpmax/xd)**dels,gb10
504 if(ish.ge.5)write (ifch,*)'gb12,xd,xpmax,wwq,wwg:'
505 if(ish.ge.5)write (ifch,*)gb12,xd,xpmax,wwq,wwg
506 if(rangen().gt.gb12)goto 4
507 endif
508
509 if(jcasc.ne.0)then
510 gb20=(1.-xd/xpmax)**betpom*sdt*(1.-glusea)+
511 * EsoftQZero(xd/xpmax)*(sdnq+sdtq)*glusea
512 else
513 gb20=EsoftQZero(xd/xpmax)
514 endif
515 if(1.+2.*(-alpqua)+dels.ge.0.)then
516 xpminl=(1.-xpmax)**(alplea(iclpro)+1.)
517 xpmaxl=(1.-xpmin)**(alplea(iclpro)+1.)
518
5195 xp=1.-(xpminl+(xpmaxl-xpminl)*rangen())**
520 * (1./(alplea(iclpro)+1.))
521 if(jcasc.ne.0)then
522 gb2=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
523 * EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*(xp/xpmax)**
524 * (1.+2.*(-alpqua)+dels)/gb20
525 else
526 gb2=EsoftQZero(xd/xp)*(xp/xpmax)**(1.+2.*(-alpqua)+dels)/gb20
527 endif
528 if(gb2.gt.1..and.ish.ge.1)then
529 write(ifmt,*)'gb2,xp:',gb2,xp
530c read (*,*)
531 endif
532 if(rangen().gt.gb2)goto 5
533 else
534 xpmaxl=xpmax**(2.+2.*(-alpqua)+dels)
535 xpminl=xpmin**(2.+2.*(-alpqua)+dels)
536
5376 xp=(xpminl+(xpmaxl-xpminl)*rangen())**
538 * (1./(2.+2.*(-alpqua)+dels))
539 if(jcasc.ne.0)then
540 gb21=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
541 * EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*
542 * ((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
543 else
544 gb21=EsoftQZero(xd/xp)*((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
545 endif
546 if(gb21.gt.1..and.ish.ge.1)then
547 write(ifmt,*)'gb21,xp:',gb21,xp
548c read (*,*)
549 endif
550 if(rangen().gt.gb21)goto 6
551 endif
552
553 wwh=xd*wd-qq
554 wwsh=xp*wd-qq
555 ammax=(egyevt-sqrt(wwsh))**2
55622 call fremnx(ammax,ammin,sm,icp3,icp4,iret)
557 if(iret.ne.0.and.ish.ge.1)write(ifmt,*)'iret.ne.0!'
558 * ,ammax,ammin**2
559 wmn=(1.-xp)*wd/wp0
560 wpn=sm/wmn
561 pnx=0.
562 pny=0.
563 wpp=wp0-wpn
564 wmp=wm0-wmn
565 if(ish.ge.5)write(ifch,*)'wp0,wm0,wpn,wmn,wpp,wmp:'
566 if(ish.ge.5)write(ifch,*)wp0,wm0,wpn,wmn,wpp,wmp
567
568 if(jcasc.eq.0.or.rangen().lt.wwq/(wwg+wwq).
569 * and.xd*wd.gt.sqmin.and.wwsh.gt.
570 * (sqrt(wwh)+sqrt(s00))**2)then
571 zgmin=xd/xp
572 zgmax=1./(1.+wp0/xd/wd/(wpp-wwh/wmp))
573 if(zgmin.gt.zgmax)goto 22
57423 zg=zgmin-rangen()*(zgmin-zgmax)
575 if(rangen().gt.zg**dels*((1.-xd/xp/zg)/ (1.-xd/xp))**betpom)
576 * goto 23
577 xg=xd/zg !w- share for the struck quark
578 wmq=wd/wp0*(xg-xd) !w- for its counterpart
579 wpq=s00/wmq !1. gev^2 / wmq
580 wmq=0.
581 wpp=wpp-wpq
582 wmp=wmp-wmq
583 sxx=wpp*wmp
584 if(ish.ge.5)write (ifch,*)'wpq,wmq,wpp,wmp,sxx:'
585 if(ish.ge.5)write (ifch,*)wpq,wmq,wpp,wmp,sxx
586
587 if(jcasc.eq.0)then
588 if(ish.ge.6)write (ifch,*)'before call timsh2: qq,sxx,iq1',
589 * qq,sxx,iq1
590 call timsh2(qq,0.,sqrt(sxx),iq1,-iq1)
591 ept(1)=.5*(wpp+wmp)
592 ept(2)=.5*(wpp-wmp)
593 ept(3)=0.
594 ept(4)=0.
595 call psdeftr(sxx,ept,ey)
596 ep3(1)=pprt(4,2)
597 ep3(2)=pprt(3,2)
598 ep3(3)=0.
599 ep3(4)=0.
600
601 call pstrans(ep3,ey,1)
602 wmp=ep3(1)-ep3(2)
603 goto 24
604 endif
605 else
606 iq1=0
607 sxx=wpp*wmp
608 endif
609
610 if(ish.ge.3)write (ifch,*)'wwh,wwsh,sxx,wpp,wmp:',
611 * wwh,wwsh,sxx,wpp,wmp
612
613 wpi=wpp
614 wmi=wwh/wpp
615 wmp=wmp-wmi
61624 call fremny(wpn,wmn,pnx,pny,sm,icp1,icp2,icp3,icp4,bx,ey0)
617
618 if((-alpqua).eq.-1.)stop'dis does not work for 1/x'
61925 aks=rangen()
620 z=.5*aks**(1./(1.+(-alpqua)))
621 if(z.lt.1.e-5.or.rangen().gt.(2.*(1.-z))**(-alpqua))goto 25
622 if(rangen().gt..5)z=1.-z
623 wm2=wmp*z
624 wm1=wmp-wm2
625
626 iqc(2)=iq1
627 nj=nj+1
628 iqj(nj)=-int(2.*rangen()+1.)
629 iqj(nj+1)=-iqj(nj)
630 eqj(1,nj)=.5*wm1
631 eqj(2,nj)=-eqj(1,nj)
632 eqj(3,nj)=0.
633 eqj(4,nj)=0.
634 eqj(1,nj+1)=.5*wm2
635 eqj(2,nj+1)=-eqj(1,nj+1)
636 eqj(3,nj+1)=0.
637 eqj(4,nj+1)=0.
638 nj=nj+1
639
640 if(iq1.eq.0)then
641 ncc(1,2)=nj-1
642 ncc(2,2)=nj
643 gdres=sdt-sds
644 gdrga=sdr
645 gdsin=sds
646 gdbor=sdb
647 gdnon=sdn
648 else
649 nj=nj+1
650 if(iabs(iq1).eq.3)then
651 iqj(nj)=-iq1*4/3
652 else
653 iqj(nj)=-iq1
654 endif
655 eqj(1,nj)=.5*(wpq+wmq)
656 eqj(2,nj)=.5*(wpq-wmq)
657 eqj(3,nj)=0.
658 eqj(4,nj)=0.
659 if(iq1.gt.0)then
660 ncj(1,nj)=nj-1
661 ncj(1,nj-1)=nj
662 ncj(2,nj)=0
663 ncj(2,nj-1)=0
664 else
665 ncj(1,nj)=nj-2
666 ncj(1,nj-2)=nj
667 ncj(2,nj)=0
668 ncj(2,nj-2)=0
669 endif
670
671 if(jcasc.eq.0)then
672 if(iq1.gt.0)then
673 nqc(1)=nj-2
674 nqc(2)=0
675 else
676 nqc(1)=nj-1
677 nqc(2)=0
678 endif
679 s0h=0.
680 c0h=1.
681 s0xh=0.
682 c0xh=1.
683 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
684 goto 17
685 else
686 gdres=(sdtq-sdsq)/4.5
687 gdrga=sdrq/4.5
688 gdsin=sdsq/4.5
689 gdbor=sdbq/4.5
690 gdnon=sdnq/4.5
691 if(iq1.gt.0)then
692 ncc(1,2)=nj-2
693 ncc(2,2)=0
694 else
695 ncc(1,2)=nj-1
696 ncc(2,2)=0
697 endif
698 endif
699 endif
700
701 if(ish.ge.3)write (ifch,*)'wpn,wmn,wpi,wmi,wm1,wm2,nj'
702 if(ish.ge.3)write (ifch,*)wpn,wmn,wpi,wmi,wm1,wm2,nj
703 endif
704
705 si=wpi*wmi+qq
706 qmin(2)=q2min !effective momentum cutoff below
707 s2min=max(4.*qq,16.*q2min) !mass cutoff for born scattering
708
709 if(rangen().gt.gdres/(gdres+gdsin+gdnon).or.
710 *si.lt.(s2min+qq))goto 12
711
712c---------------------------------------
713c hard pomeron (resolved photon)
714c---------------------------------------
715 if(ish.ge.3)write(ifmt,*)'resolved,gdrga,gdres',gdrga,gdres
716
717 jj=1
718 if(rangen().gt.gdrga/gdres.and.si.gt.1.1*s2min+qq)then
719 if(ish.ge.3)write(ifmt,*)'dir-res,si,qq',si,qq
720 pt=0.
721 pt2=0.
722 iqc(1)=0
723 ept(1)=.5*(wpi+wmi)
724 ept(2)=.5*(wpi-wmi)
725 ept(3)=0.
726 ept(4)=0.
727 wpt(1)=wpi !lc+ for the current jet emission
728 wpt(2)=si/wpi !lc- for the current jet emission
729
730 qqmin=max(q2min,s2min/(si/qq-1.))
731 qqmax=min(si/2.,si-s2min)
732 qmin(1)=qqmin
733 xmax=1.
734 xmin=(s2min+qq)/si
735 if(qqmin.ge.qqmax.or.xmin.ge.xmax)stop'min>max'
736 gb0=psjti(qmin(1),qq,si-qq,7,iqc(2),1)*psfap(1.d0,0,1)
737
738 ncc(1,1)=0
739 ncc(2,1)=0
740 jgamma=1
741 ntry=0
742 goto 9
743 else
744 if(ish.ge.3)write(ifmt,*)'res,si,qq',si,qq
745 qmin(1)=q2min !effective momentum cutoff above
746 si=si-qq
747 zmin=s2min/si
748 dft0=psdfh4(zmin,q2min,qq,0,0)*psjti(q2min,qq,si,0,iqc(2),1)
749 * +(psdfh4(zmin,q2min,qq,0,1)+psdfh4(zmin,q2min,qq,0,2)+
750 * psdfh4(zmin,q2min,qq,0,3))*psjti(q2min,qq,si,7,iqc(2),1)
751
7527 continue
753 z=zmin**rangen()
754 do i=1,4
755 dfp(i)=psdfh4(z,q2min,qq,0,i-1)
756 enddo
757 dfp(1)=dfp(1)*psjti(q2min,qq,z*si,0,iqc(2),1)
758 dfptot=dfp(1)
759 if(iqc(2).eq.0)then
760 sjq=psjti(q2min,qq,z*si,1,0,1)
761 do i=2,4
762 dfp(i)=dfp(i)*sjq
763 dfptot=dfptot+dfp(i)
764 enddo
765 else
766 sjqqp=psjti(q2min,qq,z*si,1,2,1)
767 do i=2,4
768 if(iabs(iqc(2)).eq.i-1)then
769 dfp(i)=dfp(i)*(psjti(q2min,qq,z*si,1,1,1)+
770 * psjti(q2min,qq,z*si,1,-1,1))/2.
771 else
772 dfp(i)=dfp(i)*sjqqp
773 endif
774 dfptot=dfptot+dfp(i)
775 enddo
776 endif
777
778 if(rangen().gt.dfptot/dft0)goto 7
779
780 wpq=wpi*(1.-z)
781 wpi=wpi*z
782 aks=dfptot*rangen()
783 if(aks.lt.dfp(1))then
784 iqc(1)=0
785 nj=nj+1
786 ncc(1,1)=nj
787 ncc(2,1)=nj+1
788
789 iqj(nj)=-int(2.*rangen()+1.)
790 iqj(nj+1)=-iqj(nj)
791 wpq1=wpq*rangen()
792 eqj(1,nj)=.5*wpq1
793 eqj(2,nj)=eqj(1,nj)
794 eqj(3,nj)=0.
795 eqj(4,nj)=0.
796 eqj(1,nj+1)=.5*(wpq-wpq1)
797 eqj(2,nj+1)=eqj(1,nj+1)
798 eqj(3,nj+1)=0.
799 eqj(4,nj+1)=0.
800 nj=nj+1
801
802 else
803 if(aks.lt.dfp(1)+dfp(2))then
804 iqc(1)=1
805 elseif(aks.lt.dfp(1)+dfp(2)+dfp(3))then
806 iqc(1)=2
807 else
808 iqc(1)=3
809 endif
810 iqc(1)=iqc(1)*(2*int(2.*rangen())-1)
811 nj=nj+1
812 ncc(1,1)=nj
813 ncc(2,1)=0
814
815 iqj(nj)=-iqc(1)
816 eqj(1,nj)=.5*wpq
817 eqj(2,nj)=eqj(1,nj)
818 eqj(3,nj)=0.
819 eqj(4,nj)=0.
820 endif
821
822 ept(1)=.5*(wpi+wmi)
823 ept(2)=.5*(wpi-wmi)
824 ept(3)=0.
825 ept(4)=0.
826 jgamma=0
827 ntry=0
828 endif
829
8308 continue
831
832c ladder rung
833c---------------------------------------
834 pt2=ept(3)**2+ept(4)**2
835 pt=sqrt(pt2)
836
837 wpt(1)=ept(1)+ept(2) !lc+ for the current jet emissi
838 wpt(2)=ept(1)-ept(2) !lc- for the current jet emissi
839
840 s2min=max(qmin(1),16.*qmin(2)) !mass cutoff for born
841 s2min=max(s2min,4.*qq)
842
843 if(jj.eq.1)then
844 wwmin=2.*s2min-2.*pt*sqrt(q2ini)
845 wwmin=(wwmin+sqrt(wwmin**2+4.*pt2*(s2min-q2ini)))
846 * /(1.-q2ini/s2min)/2.
847 sj=psjti(qmin(1),qq,si,iqc(1),iqc(2),1) !total jet
848 sj2=psjti1(q2min,qmin(1),qq,si,iqc(2),iqc(1),1)
849 if(ish.ge.3)write(ifch,*)'resolved - si,wwmin,s2min,sj,sj2:'
850 if(ish.ge.3)write(ifch,*)si,wwmin,s2min,sj,sj2
851 if(sj.eq.0.)stop'sj=0'
852 if(rangen().gt.sj2/sj.and.si.gt.1.1*wwmin)goto 26
853 jj=2
854 endif
855 sj=psjti1(qmin(2),qmin(1),qq,si,iqc(2),iqc(1),1)
856 sjb=psbint(qmin(1),qmin(2),qq,si,iqc(1),iqc(2),1) !born parton-parton
857 wwmin=17./16*s2min-2.*pt*sqrt(q2ini)
858 wwmin=(wwmin+sqrt(wwmin**2+pt2*(s2min/4.-4.*q2ini)))
859 */(1.-16.*q2ini/s2min)/2.
860 if(rangen().lt.sjb/sj.or.si.lt.1.1*wwmin)goto 10
861
86226 continue
863 wpt(jj)=wpt(jj)-pt2/wpt(3-jj)
864
865 if(jj.eq.1)then
866 discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(2.*si+pt2)
867 if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
868 * discr,si,pt,wwmin
869 discr=sqrt(discr)
870 qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(2.+pt2/si)
871 else
872 discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(17.*si+pt2)
873 if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
874 * discr,si,pt,wwmin
875 discr=sqrt(discr)
876 qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(17.+pt2/si)
877 endif
878 qqmin=2.*q2ini*si/(si+2.*pt*sqrt(q2ini)+discr)
879 if(jj.eq.1.and.s2min.gt.qqmin.or.
880 *jj.eq.2.and.s2min.gt.16.*qqmin)then
881 xmm=.5*(si-s2min+2.*pt*sqrt(q2ini))
882 discr=xmm**2-q2ini*(si+pt2)
883 if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
884 * discr,si,pt,wwmin
885 qqmin=q2ini*si/(xmm+sqrt(discr))
886 endif
887
888 xmin=1.-q2ini/qqmin
889 xmax=1.-q2ini/qqmax
890 if(ish.ge.6)write(ifch,*)'qqmin,qqmax,xmin,xmax',
891 *qqmin,qqmax,xmin,xmax
892 if(qqmin.lt.qmin(jj))then
893 qqmin=qmin(jj)
894 xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
895 * si*(si-s2min-qqmin*(1.+pt2/si))))/si)**2,
896 * (s2min+qqmin*(1.+pt2/si)-2.*pt*sqrt(qqmin))/si)
897 xmin=max(xmin,xmi)
898 if(xmin.le.0.)xmin=(s2min+qqmin*(1.+pt2/si))/si
899 if(ish.ge.6)write(ifch,*)'qqmin,qmin(jj),xmin,s2min',
900 * qqmin,qmin(jj),xmin,s2min
901 endif
902
903 qm0=qmin(jj)
904 xm0=1.-q2ini/qm0
905 if(xm0.gt.xmax.or.xm0.lt.xmin)then
906 xm0=.5*(xmax+xmin)
907 endif
908c s2max=xm0*si
909 s2max=xm0*si-qm0*(1.+pt2/si)+2.*pt*sqrt(q2ini) !new ladder mass squared
910 xx=xm0
911
912 if(jj.eq.1)then
913 sj0=psjti(qm0,qq,s2max,0,iqc(2),1)*psfap(xx,iqc(1),0)+
914 * psjti(qm0,qq,s2max,7,iqc(2),1)*psfap(xx,iqc(1),1)
915 gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(1)))*qm0*2.
916 else
917 sj0=psjti1(qm0,qmin(1),qq,s2max,0,iqc(1),1)*psfap(xx,iqc(2),0)
918 * +psjti1(qm0,qmin(1),qq,s2max,7,iqc(1),1)*psfap(xx,iqc(2),1)
919 gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*2.
920 endif
921 if(gb0.le.0.)then
922 write(ifmt,*)'gb0.le.0. si,qq,pt2:',si,qq,pt2
923 iret=1
924 goto 9999
925 endif
926 if(xm0.le..5)then
927 gb0=gb0*xm0**(1.-delh)
928 else
929 gb0=gb0*(1.-xm0)*2.**delh
930 endif
931
932 xmin2=max(.5,xmin)
933 xmin1=xmin**delh !xmin, xmax are put into powe
934 xmax1=min(xmax,.5)**delh !to simulate x value below
935 if(xmin.ge..5)then
936 djl=1.
937 elseif(xmax.lt..5)then
938 djl=0.
939 else
940 djl=1./(1.+((2.*xmin)**delh-1.)/delh/
941 * log(2.*(1.-xmax)))
942 endif
943
944 ntry=0
9459 continue
946 ntry=ntry+1
947 if(ntry.ge.10000)then
948 print *,"ntry.ge.10000"
949 iret=1
950 goto 9999
951 endif
952 if(jgamma.ne.1)then
953 if(rangen().gt.djl)then !lc momentum share in the cur
954 x=(xmin1+rangen()*(xmax1-xmin1))**(1./delh)
955 else
956 x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
957 endif
958 q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
959 qt2=q2*(1.-x)
960 if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
961 if(qt2.lt.q2ini)goto 9
962 else
963 x=xmin+rangen()*(xmax-xmin)
964 q2=qqmin*(qqmax/qqmin)**rangen()
965 qt2=(q2-x*qq)*(1.-x)
966 if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
967 if(qt2.lt.0.)goto 9
968 endif
969
970 qt=sqrt(qt2)
971 call pscs(bcos,bsin)
972c ep3 is now 4-vector for s-channel gluon produced in current ladder run
973 ep3(3)=qt*bcos
974 ep3(4)=qt*bsin
975 ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
976 if(jj.eq.1)then
977 s2min2=max(q2,s2min)
978 else
979 s2min2=max(s2min,16.*q2)
980 endif
981
982 if(jgamma.ne.1)then
983 s2=x*si-q2*(1.+pt2/si)-ptnew+pt2+qt2 !new ladder mass squared
984 if(s2.lt.s2min2)goto 9 !rejection in case of too low mass
985 xx=x
986
987 if(jj.eq.1)then
988 sj1=psjti(q2,qq,s2,0,iqc(2),1)
989 if(iqc(1).ne.0)then
990 sj2=psjti(q2,qq,s2,iqc(1),iqc(2),1)
991 elseif(iqc(2).eq.0)then
992 sj2=psjti(q2,qq,s2,1,0,1)
993 else
994 sj2=psjti(q2,qq,s2,1,1,1)/6.+
995 * psjti(q2,qq,s2,-1,1,1)/6.+
996 * psjti(q2,qq,s2,2,1,1)/1.5
997 endif
998 else
999 sj1=psjti1(q2,qmin(1),qq,s2,0,iqc(1),1)
1000 if(iqc(2).ne.0)then
1001 sj2=psjti1(q2,qmin(1),qq,s2,iqc(2),iqc(1),1)
1002 elseif(iqc(1).eq.0)then
1003 sj2=psjti1(q2,qmin(1),qq,s2,1,0,1)
1004 else
1005 sj2=psjti1(q2,qmin(1),qq,s2,1,1,1)/6.+
1006 * psjti1(q2,qmin(1),qq,s2,-1,1,1)/6.+
1007 * psjti1(q2,qmin(1),qq,s2,2,1,1)/1.5
1008 endif
1009 endif
1010c gb7 is the rejection function for x and q**2 simulation
1011 gb7=(sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))
1012 * /log(qt2/qcdlam)*sngl(psuds(q2,iqc(jj)))*q2/gb0
1013
1014 if(x.le..5)then
1015 gb7=gb7*x**(1.-delh)
1016 else
1017 gb7=gb7*(1.-x)*2.**delh
1018 endif
1019 else
1020 s2=x*si-q2 !new ladder mass squared
1021 if(s2.lt.s2min2)goto 9 !rejection in case of too low mass
1022
1023 sj1=0.
1024 xx=x
1025 if(iqc(2).eq.0)then
1026 sj2=psjti(q2,qq,s2,1,0,1)
1027 else
1028 sj2=psjti(q2,qq,s2,1,1,1)/naflav/2.+
1029 * psjti(q2,qq,s2,-1,1,1)/naflav/2.+
1030 * psjti(q2,qq,s2,2,1,1)*(1.-1./naflav)
1031 endif
1032 gb7=sj2*psfap(xx,0,1)/gb0 !????*(1.-x*qq/q2)
1033 endif
1034 if(gb7.gt.1..or.gb7.lt.0..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,gb0'
1035 *,gb7,q2,x,gb0
1036 if(rangen().gt.gb7)goto 9
1037
1038 if(ish.ge.6)write(ifch,*)'res: jj,iqc,ncc:',
1039 *jj,iqc(jj),ncc(1,jj),ncc(2,jj)
1040
1041 nqc(2)=0
1042 iqnew=iqc(jj)
1043 if(jgamma.ne.1)then
1044 if(rangen().lt.sj1/(sj1+sj2))then
1045 if(iqc(jj).eq.0)then
1046 jt=1
1047 jq=int(1.5+rangen())
1048 nqc(1)=ncc(jq,jj)
1049 else
1050 jt=2
1051 if(iqc(jj).gt.0)then
1052 jq=1
1053 else
1054 jq=2
1055 endif
1056 nqc(1)=0
1057 iqnew=0
1058 endif
1059 iq1=iqc(jj)
1060 else
1061 if(iqc(jj).ne.0)then
1062 iq1=0
1063 jt=3
1064 if(iqc(jj).gt.0)then
1065 jq=1
1066 else
1067 jq=2
1068 endif
1069 nqc(1)=ncc(1,jj)
1070 else
1071 jt=4
1072 jq=int(1.5+rangen())
1073 iq1=int(naflav*rangen()+1.)*(3-2*jq)
1074 nqc(1)=ncc(jq,jj)
1075 iqnew=-iq1
1076 endif
1077 endif
1078 else
1079 jt=5
1080 jq=int(1.5+rangen())
1081 iq1=int(naflav*rangen()+1.)*(3-2*jq)
1082 iqnew=-iq1
1083 nqc(1)=0
1084 endif
1085 eprt=max(1.d0*qt,
1086 *.5d0*((1.d0-x)*wpt(jj)+qt2/(1.d0-x)/wpt(jj)))
1087 pl=((1.d0-x)*wpt(jj)-eprt)*(3-2*jj)
1088 if(iq1.eq.0)then
1089 iq2ini=9
1090 else
1091 iq2ini=iq1
1092 endif
109327 call timsh1(q2,sngl(eprt),iq2ini)
1094 amprt=pprt(5,1)**2
1095 plprt=eprt**2-amprt-qt2
1096 if(plprt.lt.-1d-6)goto 27
1097 ep3(1)=eprt
1098 ep3(2)=dsqrt(max(0.d0,plprt))
1099 if(pl.lt.0.d0)ep3(2)=-ep3(2)
1100 ey(1)=1.
1101 ey(2)=1.
1102 ey(3)=1.
1103 do i=1,4
1104 ept1(i)=ept(i)-ep3(i)
1105 enddo
1106 call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1107 if(ish.ge.6)then
1108 write(ifch,*)'q2,amprt,qt2',q2,amprt,qt2
1109 write(ifch,*)'eprt,plprt',eprt,plprt
1110 write(ifch,*)'ep3',ep3
1111 write(ifch,*)'ept',ept
1112 write(ifch,*)'ept1',ept1
1113 endif
1114 s2new=psnorm(ept1)
1115
1116 if(s2new.gt.s2min2)then
1117 if(jj.eq.1)then
1118 gb=psjti(q2,qq,s2new,iqnew,iqc(2),1)
1119 else
1120 gb=psjti1(q2,qmin(1),qq,s2new,iqnew,iqc(1),1)
1121 endif
1122 if(iqnew.eq.0)then
1123 gb=gb/sj1
1124 else
1125 gb=gb/sj2
1126 endif
1127 if(ish.ge.1)then
1128 if(gb.gt.1.)write (ifch,*)'gb,s2new,s2,q2,iqnew',
1129 * gb,s2new,s2,q2,iqnew
1130 endif
1131 if(rangen().gt.gb)goto 9
1132 else
1133 goto 9
1134 endif
1135 jgamma=0
1136
1137 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1138
1139 if(jt.eq.1)then
1140 ncc(jq,jj)=nqc(2)
1141 elseif(jt.eq.2)then
1142 ncc(jq,jj)=ncc(1,jj)
1143 ncc(3-jq,jj)=nqc(1)
1144 elseif(jt.eq.3)then
1145 ncc(1,jj)=nqc(2)
1146 elseif(jt.eq.4)then
1147 ncc(1,jj)=ncc(3-jq,jj)
1148 elseif(jt.eq.5)then
1149 ncc(1,jj)=nqc(1)
1150 ncc(2,jj)=0
1151 endif
1152 iqc(jj)=iqnew
1153 if(ish.ge.6)write(ifch,*)'qt2,amprt,ncc:',
1154 *qt2,amprt,ncc(1,jj),ncc(2,jj)
1155
1156 do i=1,4
1157 ept(i)=ept1(i)
1158 enddo
1159c c.m. energy squared, minimal 4-momentum transfer square and gluon 4-v
1160c for the next ladder run
1161 qmin(jj)=q2
1162 si=s2new
1163 if(ish.ge.3)write (ifch,*)'res: new jet - iqj,ncj,ep3,ept',
1164 *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1165
1166 goto 8 !next simulation step will be considered
1167
116810 continue
1169 if(ish.ge.3)write(ifch,*)'res: iqc,si,ept:',iqc,si,ept
1170
1171c highest virtuality subprocess in the ladder
1172c---------------------------------------
1173 qqs=max(qmin(1)/4.,4.*qmin(2))
1174 qqs=max(qqs,qq)
1175 call psabor(si,qqs,iqc,ncc,ept,1,nptlh,bx)
1176 goto 17
1177
117812 continue
1179c---------------------------------------
1180c hard pomeron (direct photon)
1181c---------------------------------------
1182 ept(1)=.5*(wpi+wmi)
1183 ept(2)=.5*(wpi-wmi)
1184 ept(3)=0.
1185 ept(4)=0.
1186 if(ish.ge.3)write (ifch,*)'direct photon - ept,si,qq:',ept,si,qq
1187
118813 continue
1189
1190c ladder rung
1191c---------------------------------------
1192 pt2=ept(3)**2+ept(4)**2
1193 pt=sqrt(pt2)
1194 wpt(1)=ept(1)+ept(2)
1195 wpt(2)=si/wpt(1)
1196
1197 gdbor=psdbin(qmin(2),qq,si,iqc(2),long)
1198 gdtot=psdsin(qmin(2),qq,si,iqc(2),long)
1199 if(iqc(2).ne.0)then
1200 if(ish.ge.8)write (ifch,*)'qmin(2),qq,si',qmin(2),qq,si
1201 gdnon=psdnsi(qmin(2),qq,si,long)
1202 if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1203 gdbor=gdbor/2.25
1204 gdtot=gdnon/2.25+gdtot/4.5
1205 else
1206 gdbor=gdbor/9.
1207 gdtot=gdnon/9.+gdtot/4.5
1208 endif
1209 else
1210 gdnon=0.
1211 endif
1212
1213 if(long.ne.0.or.qmin(2).ge.qq)then
1214 s2min=qq+4.*max(qmin(2),qcmass**2)
1215 wwmin=(5.*s2min-qq)/4.-2.*pt*sqrt(q2ini)
1216 wwmin=(wwmin+sqrt(wwmin**2-(qq-pt2)*(s2min-qq-4.*q2ini)))
1217 * /2./(1.-4.*q2ini/(s2min-qq))
1218 else
1219 s2min=qq/(1.-sqrt(q2ini/qq))
1220 wwmin=s2min+qq-2.*pt*sqrt(q2ini)
1221 wwmin=(wwmin+sqrt(wwmin**2-4.*(qq-pt2)*(qq-q2ini)))
1222 * /2./(1.-q2ini/qq)
1223 endif
1224
1225 if(ish.ge.3)write(ifch,*)'si,s2min,wwmin,qmin(2),gdtot,gdbor:'
1226 if(ish.ge.3)write(ifch,*)si,s2min,wwmin,qmin(2),gdtot,gdbor
1227
1228 if((rangen().lt.gdbor/gdtot.or.si.lt.1.1*wwmin).and.
1229 *(long.eq.0.and.qmin(2).lt.qq.or.iqc(2).eq.0))goto 15
1230 if(si.lt.1.1*wwmin)stop'si<1.1*wwmin'
1231
1232 qqmax=0.
1233 qqmin=0.
1234
1235 xmm=si+2.*sqrt(q2ini)*pt-qq
1236 discr=xmm**2-4.*q2ini*(5.*si-qq+pt2)
1237 if(discr.lt.0.)goto 29
1238 discr=sqrt(discr)
1239 qqmax=(xmm+discr)/2./(5.-(qq-pt2)/si)
1240 qqmin=2.*q2ini*si/(xmm+discr)
1241
124229 continue
1243 if(4.*qqmin.lt.s2min-qq.or.long.eq.0.and.
1244 *qmin(2).lt.qq)then
1245 xmm=si-s2min+2.*sqrt(q2ini)*pt
1246 qqmin=2.*q2ini*si/(xmm+sqrt(xmm**2-4.*q2ini*(si-qq+pt2)))
1247 endif
1248 xmin=1.-q2ini/qqmin
1249
1250 if(qqmin.lt.qmin(2))then
1251 qqmin=qmin(2)
1252 xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
1253 * si*(si-s2min-qqmin*(1.-(qq-pt2)/si))))/si)**2,
1254 * (s2min+qqmin*(1.-(qq-pt2)/si)-2.*pt*sqrt(qqmin))/si)
1255 xmin=max(xmin,xmi)
1256 endif
1257 if(xmin.le.qq/si)xmin=1.001*qq/si
1258
1259 if(long.eq.0.and.qmin(2).lt.qq)qqmax=max(qqmax,qq)
1260 xmax=1.-q2ini/qqmax
1261
1262 if(ish.ge.6)write(ifch,*)'qqmax,qqmin,xmax,xmin:',
1263 *qqmax,qqmin,xmax,xmin
1264 if(qqmax.lt.qqmin)stop'qqmax<qqmin'
1265
1266 qm0=qqmin
1267 xm0=1.-q2ini/qm0
1268 s2max=si*xm0-qm0*(1.-qq/si)
1269
1270 sds=psdsin(qm0,qq,s2max,0,long)/4.5
1271 sdv=psdsin(qm0,qq,s2max,1,long)/4.5
1272
1273 sdn=psdnsi(qm0,qq,s2max,long)
1274 if(iqc(2).eq.0)then
1275 sdn=sdn/4.5
1276 elseif(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1277 sdn=sdn/2.25
1278 else
1279 sdn=sdn/9.
1280 endif
1281 sdv=sdv+sdn
1282 xx=xm0
1283
1284 sj0=sds*psfap(xx,iqc(2),0)+sdv*psfap(xx,iqc(2),1)
1285 gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*5.
1286 if(gb0.le.0.)then
1287 write(ifmt,*)'gb0.le.0. si,qq,pt2:',si,qq,pt2
1288 iret=1
1289 goto 9999
1290 endif
1291
1292 if(xm0.le..5)then
1293 gb0=gb0*(xm0-qq/si)/(1.-2.*qq/si)
1294 else
1295 gb0=gb0*(1.-xm0)
1296 endif
1297
1298 xmin2=max(.5,xmin)
1299 xmax1=min(xmax,.5)
1300 if(xmin.ge..5)then
1301 djl=1.
1302 elseif(xmax.lt..5)then
1303 djl=0.
1304 else
1305 djl=1./(1.-(1.-2.*qq/si)*log((.5-qq/si)/(xmin-qq/si))/
1306 * log(2.*(1.-xmax)))
1307 endif
1308
130914 continue
1310 if(rangen().gt.djl)then !lc momentum share in the cur
1311 x=(xmin-qq/si)*((xmax1-qq/si)/(xmin-qq/si))**rangen()+qq/si
1312 else
1313 x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
1314 endif
1315
1316 q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
1317
1318 qt2=q2*(1.-x)
1319 if(ish.ge.9)write(ifch,*)'q2,x,qt2,qq,qqmin,qqmax:',
1320 *q2,x,qt2,qq,qqmin,qqmax
1321 if(qt2.lt.q2ini)goto 14 !p_t check
1322
1323 if(long.ne.0.or.q2.ge.qq)then
1324 s2min2=max(4.*q2+qq,s2min)
1325 else
1326 s2min2=s2min
1327 endif
1328 qt=sqrt(qt2)
1329 call pscs(bcos,bsin)
1330c ep3 is now 4-vector for s-channel gluon produced in current ladder run
1331 ep3(3)=qt*bcos
1332 ep3(4)=qt*bsin
1333 ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
1334
1335 s2=x*si-ptnew+pt2-q2*(x-(qq-pt2)/si)
1336 if(s2.lt.s2min2)goto 14 !check of the kinematics
1337 sds=psdsin(q2,qq,s2,0,long)/4.5
1338 sdv0=psdsin(q2,qq,s2,1,long)/4.5
1339 if(ish.ge.8)write (ifch,*)'q2,qq,s2',q2,qq,s2
1340 sdn0=psdnsi(q2,qq,s2,long)
1341
1342 if(iqc(2).eq.0)then
1343 sdn=sdn0/4.5
1344 else
1345 if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1346 sdn=sdn0/2.25
1347 else
1348 sdn=sdn0/9.
1349 endif
1350 endif
1351 sdv=sdv0+sdn
1352
1353 xx=x
1354 sj1=sds*psfap(xx,iqc(2),0)
1355 sj2=sdv*psfap(xx,iqc(2),1)
1356
1357c gb7 is the rejection function for x and q**2 simulation.
1358 gb7=(sj1+sj2)/log(qt2/qcdlam)*sngl(psuds(q2,iqc(2)))/gb0*q2
1359 if(x.le..5)then
1360 gb7=gb7*(x-qq/si)/(1.-2.*qq/si)
1361 else
1362 gb7=gb7*(1.-x)
1363 endif
1364
1365 if(gb7.gt.1..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,qt2,iqc(2),'
1366 * ,'gb0,sj1,sj2',gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2
1367 if(ish.ge.3)write (ifch,*)'gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long',
1368 * gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long
1369 if(rangen().gt.gb7)goto 14
1370
1371
1372 if(ish.ge.6)write(ifch,*)'iqc,ncc:',iqc(2),ncc(1,2),ncc(2,2)
1373 iqcnew=iqc(2)
1374 nqc(2)=0 !emitted parton color connections
1375 if(rangen().lt.sj1/(sj1+sj2).or.(long.ne.0.or.q2.ge.qq).and.
1376 *s2.lt.1.5*s2min2)then
1377 if(iqc(2).eq.0)then
1378 jt=1
1379 jq=int(1.5+rangen())
1380 nqc(1)=ncc(jq,2)
1381 else
1382 jt=2
1383 if(iqc(2).gt.0)then
1384 jq=1
1385 else
1386 jq=2
1387 endif
1388 nqc(1)=0
1389 endif
1390 iq1=iqc(2)
1391 iqcnew=0
1392
1393 else
1394 if(iqc(2).ne.0)then
1395 jt=3
1396 iq1=0
1397 if(iqc(2).gt.0)then
1398 jq=1
1399 else
1400 jq=2
1401 endif
1402 nqc(1)=ncc(1,2)
1403
1404 else
1405 tu=sdn0/2.25+sdv0
1406 if(naflav.eq.4)tu=tu*2.
1407 td=sdn0/9.+sdv0
1408 if(rangen().lt.tu/(tu+2.*td))then
1409 if(naflav.eq.3)then
1410 iq1=1
1411 else
1412 iq1=1+3*int(.5+rangen())
1413 endif
1414 else
1415 iq1=int(2.5+rangen())
1416 endif
1417 jq=int(1.5+rangen())
1418 iq1=iq1*(3-2*jq)
1419 iqcnew=-iq1
1420 jt=4
1421 nqc(1)=ncc(jq,2)
1422 endif
1423 endif
1424
1425 eprt=max(1.d0*qt,
1426 *.5d0*((1.d0-x)*wpt(2)+qt2/(1.d0-x)/wpt(2)))
1427 pl=eprt-(1.d0-x)*wpt(2)
1428 if(iq1.eq.0)then
1429 iq2ini=9
1430 else
1431 iq2ini=iq1
1432 endif
143328 call timsh1(q2,sngl(eprt),iq2ini)
1434 amprt=pprt(5,1)**2
1435 plprt=eprt**2-amprt-qt2
1436 if(plprt.lt.-1d-6)goto 28
1437 ep3(1)=eprt
1438 ep3(2)=dsqrt(max(0.d0,plprt))
1439 if(pl.lt.0.d0)ep3(2)=-ep3(2)
1440 ey(1)=1.
1441 ey(2)=1.
1442 ey(3)=1.
1443 do i=1,4
1444 ept1(i)=ept(i)-ep3(i)
1445 enddo
1446 call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1447 call psrotat(ep3,s0xh,c0xh,s0h,c0h)
1448 s2new=psnorm(ept1)+qq
1449
1450 if((long.ne.0.or.q2.ge.qq).and.iqcnew.ne.0)then
1451 xmm=(5.*s2min2-qq)/4.-2.*sqrt(ptnew*q2ini)
1452 s2min2=1.1*(xmm+sqrt(xmm**2-(qq-ptnew)*
1453 * (s2min2-qq-4.*q2ini)))/2./(1.-4.*q2ini/(s2min2-qq))
1454 endif
1455 if(s2new.gt.s2min2)then
1456 sds1=psdsin(q2,qq,s2new,iqcnew,long)/4.5
1457 if(iqcnew.eq.0)then
1458 gb=sds1/sds
1459 else
1460 if(ish.ge.8)write (ifch,*)'q2,qq,s2new',q2,qq,s2new
1461 sdn1=psdnsi(q2,qq,s2new,long)
1462 if(iabs(iqcnew).eq.1.or.iabs(iqcnew).eq.4)then
1463 sdn1=sdn1/2.25
1464 sdv=sdv0+sdn0/2.25
1465 else
1466 sdn1=sdn1/9.
1467 sdv=sdv0+sdn0/9.
1468 endif
1469 gb=.9999*(sds1+sdn1)/sdv
1470 endif
1471 if(ish.ge.3.and.gb.gt.1..and.ish.ge.1)write(ifmt,*)'gbs2',gb
1472 if(rangen().gt.gb)goto 14
1473 else
1474 goto 14
1475 endif
1476
1477 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1478
1479 iqc(2)=iqcnew
1480 if(jt.eq.1)then !current parton color connections
1481 ncc(jq,2)=nqc(2)
1482 elseif(jt.eq.2)then
1483 ncc(jq,2)=ncc(1,2)
1484 ncc(3-jq,2)=nqc(1)
1485 elseif(jt.eq.3)then
1486 ncc(1,2)=nqc(2)
1487 elseif(jt.eq.4)then
1488 ncc(1,2)=ncc(3-jq,2)
1489 ncc(2,2)=0
1490 endif
1491
1492 do i=1,4
1493 ept(i)=ept1(i)
1494 enddo
1495 if(ish.ge.3)write (ifch,*)'new jet - iqj,ncj,ep3,ept',
1496 *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1497c c.m. energy squared, minimal 4-momentum transfer square and gluon 4-v
1498c for the next ladder run
1499 qmin(2)=q2
1500 si=s2new
1501 goto 13 !next simulation step will be considered
1502
150315 continue
1504 if(ish.ge.3)write(ifch,*)'iqc,si,qmin(2),nj:',
1505 *iqc(2),si,qmin(2),nj
1506c highest virtuality subprocess in the ladder
1507c---------------------------------------
1508 gb01=0.
1509 tmax=0.
1510 tmin=si
1511 if(iqc(2).eq.0.and.si.gt.qq+4.*max(qcmass**2,qmin(2)))then
1512 qminn=max(qcmass**2,qmin(2))
1513 tmin1=2.*qminn/(1.-qq/si)/(1.+sqrt(1.-4.*qminn/(si-qq)))
1514 tmin=tmin1
1515 tmax=si/2.
1516 fb01=psdbom(si,si/2.,si/2.,qq,long)
1517 if(long.eq.0)fb01=fb01*si/2.
1518 gb01=fb01/log(qminn/qcdlam)*sngl(psuds(qminn,iqc(2)))/si**2
1519 gb0=gb01
1520 endif
1521
1522 if(long.eq.0.and.qmin(2).lt.qq)then
1523 tmax=max(tmax,qq)
1524 tmin=max(qmin(2),
1525 * 2.*q2ini/(1.-qq/si)/(1.+sqrt(1.-4.*q2ini/(si-qq))))
1526 ze=qq/si+tmin/si*(1.-qq/si)
1527 xx=ze
1528 qt2=tmin*(1.-ze)
1529 if(qt2.lt..999*q2ini.and.ish.ge.1)write(ifmt,*)'bor-dir:qt20'
1530 * ,qt2
1531 gb0=gb01+psfap(xx,iqc(2),1)/log(qt2/qcdlam)
1532 * *sngl(psuds(tmin,iqc(2))/psuds(tmin,1)*psuds(qq,1))
1533 * /si*(1.-tmin*qq/si**2/ze)
1534 endif
1535 gb0=gb0*2.
1536
1537 call psdeftr(si-qq,ept,ey)
1538
1539 if(ish.ge.6)write(ifch,*)'tmin,tmax,qq,si-qq,gb0:'
1540 if(ish.ge.6)write(ifch,*)tmin,tmax,qq,si-qq,psnorm(ept),gb0
1541
1542c------------------------------------------------
154316 continue
1544 if(long.eq.0)then
1545 t=tmin*(tmax/tmin)**rangen()
1546 else
1547 t=tmin+(tmax-tmin)*rangen()
1548 endif
1549
1550 u=si-t
1551 ze=qq/si+t/si*(1.-qq/si)
1552 qt2=t*(1.-ze)
1553 if(t.le.qq.and.long.eq.0)then
1554 xx=ze
1555 gb=psfap(xx,iqc(2),1)/log(qt2/qcdlam)*sngl(psuds(t,iqc(2))
1556 * /psuds(t,1)*psuds(qq,1))/si*(1.-t*qq/si**2/ze)/gb0
1557 else
1558 gb=0.
1559 endif
1560
1561 gb1=0.
1562 if(iqc(2).eq.0..and.si.gt.qq+4.*max(qcmass**2,qmin(2)).
1563 *and.qt2.gt.qcmass**2.and.t.le.si/2..and.t.ge.tmin1)then
1564 fb1=psdbom(si,t,u,qq,long)
1565 if(long.eq.0)fb1=fb1*t
1566 gb1=fb1/log(qt2/qcdlam)*sngl(psuds(qt2,iqc(2)))/si**2/gb0
1567c gb1=0. !???????????????????????
1568 gb=gb+gb1
1569 endif
1570
1571
1572 if(ish.ge.6)write(ifch,*)'gb,t,iqc(2),si,qq,qmin(2),long:',
1573 *gb,t,iqc(2),si,qq,qmin(2),long
1574 if (ish.ge.1) then
1575 if(gb.gt.1.)write(*,*)'gb,gb1,gb0,gb01',
1576 * ',t,iqc(2),si,qq,qmin(2),long:',
1577 * gb,gb1,gb0,gb01,fb1,fb01,t,iqc(2),si,qq,qmin(2),long
1578 endif
1579
1580 if(rangen().gt.gb)goto 16
1581 if(ish.ge.3)write(ifch,*)'born:t,qt2:',t,qt2
1582
1583 nqc(2)=0
1584 if(iqc(2).eq.0)then
1585 jq=int(1.5+rangen())
1586 jq2=3-jq
1587 if(rangen().gt.gb1/gb)then
1588 iq1=(1+int(3.*rangen()))*(3-2*jq)
1589 else
1590 iq1=4*(3-2*jq)
1591 endif
1592 iq2=-iq1 !quark flavors
1593 nqc(1)=ncc(jq,2)
1594 else
1595 if(iqc(2).gt.0)then
1596 jq=1
1597 else
1598 jq=2
1599 endif
1600 jq2=jq
1601 iq1=0
1602 iq2=iqc(2)
1603 nqc(1)=ncc(1,2)
1604 endif
1605
1606 call pscs(bcos,bsin)
1607 z=sngl(psutz(dble(si-qq),dble(qt2),dble(qt2)))
1608 if(t.lt..5*si)z=1.-z
1609 wp3=z*sqrt(si-qq)
1610 wm3=qt2/wp3
1611 if(iabs(iq1).eq.4)qt2=qt2-qcmass**2
1612 qt=sqrt(qt2)
1613 ep3(1)=.5*(wp3+wm3)
1614 ep3(2)=.5*(wp3-wm3)
1615 ep3(3)=qt*bcos
1616 ep3(4)=qt*bsin
1617 call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1618 if(iq1.eq.0)then
1619 iq2ini1=9
1620 else
1621 iq2ini1=iq1
1622 endif
1623 if(iq2.eq.0)then
1624 iq2ini2=9
1625 else
1626 iq2ini2=iq2
1627 endif
1628 if(ish.ge.5)write (ifch,*)'jq,jt,iq2ini1,iq2ini2',
1629 *jq,jt,iq2ini1,iq2ini2
1630
1631 if(t.lt.qq.and.iabs(iq1).ne.4)then
1632 qq1=t*(1.-ze)
1633 qq2=qq
1634 else
1635 qq1=qt2
1636 qq2=qt2
1637 endif
1638 call timsh2(qq1,qq2,sqrt(si-qq),iq2ini1,iq2ini2)
1639 nfprt=1
1640 call psreti(nqc,jq,nfprt,ey,s0xh,c0xh,s0h,c0h)
1641
1642 if(iqc(2).eq.0)then
1643 nqc(1)=ncc(3-jq,2)
1644 nqc(2)=0
1645 else
1646 nqc(1)=nqc(2)
1647 nqc(2)=0
1648 endif
1649
1650 nfprt=2
1651 call psreti(nqc,jq2,nfprt,ey,s0xh,c0xh,s0h,c0h)
1652
165317 continue
1654 if(ish.ge.3)write (ifch,*)'nj',nj
1655 if(nj.gt.0)then
1656
1657 ityj(i)=30
1658 iorj(i)=nptlh
1659 do n=1,nj
1660 do i=1,4
1661 ep3(i)=eqj(i,n)
1662 enddo
1663 call pstrans(ep3,ey0,-1) !boost to the c.m. system
1664 do i=1,4
1665 eqj(i,n)=ep3(i)
1666 enddo
1667 do l=1,6
1668 bxj(l,n)=bx(l)
1669 enddo
1670 ityj(n)=0
1671 iorj(n)=1
1672 enddo
1673 endif
1674 call psjarr(jfl) !kinky strings formation
1675 if(ish.ge.3)write (ifch,*)'jfl',jfl
1676 if(jfl.eq.0)then
1677 iret=1
1678 else
1679 iret=0
1680 ep3(4)=egyevt
1681 ep3(2)=0.
1682 ep3(3)=0.
1683 ep3(1)=0.
1684 do i=2,nptl
1685 do l=1,4
1686 ep3(l)=ep3(l)-pptl(l,i)
1687 enddo
1688 enddo
1689 if(ish.ge.3)write(ifch,*)'energy-momentum balance:'
1690 if(ish.ge.3)write(ifch,*)ep3
1691 if(abs(ep3(4)).gt.3.e-2)write(*,*)'energy-momentum balance:',ep3
1692 endif
1693 9999 call utprix('psadis',ish,ishini,3)
1694 return
1695 end
1696
1697c-----------------------------------------------------------------------
1698 subroutine psaevc
1699c-----------------------------------------------------------------------
1700 include 'epos.inc'
1701c structure functions calculation
1702 logical lcalc
1703 double precision xx,xmax
1704 dimension evs(21,21,135)
1705 common /psar2/ edmax,epmax
1706 common /psar31/ evk0(21,21,54)
1707 common /psar32/ evk(21,21,135)
1708 include 'epos.incsem'
1709
1710 inquire(file=fnie(1:nfnie),exist=lcalc)
1711 if(lcalc)then
1712 if(inicnt.eq.1)then
1713 write(ifmt,'(3a)')'read from ',fnie(1:nfnie),' ...'
1714 open(1,file=fnie(1:nfnie),status='old')
1715 read (1,*)qcdlam0,q2min0,q2ini0,naflav0,epmax0
1716 if(qcdlam0.ne.qcdlam)write(ifmt,'(a)')'iniev: wrong qcdlam'
1717 if(q2min0 .ne.q2min )write(ifmt,'(a)')'iniev: wrong q2min'
1718 if(q2ini0 .ne.q2ini )write(ifmt,'(a)')'iniev: wrong q2ini'
1719 if(naflav0.ne.naflav)write(ifmt,'(a)')'iniev: wrong naflav'
1720 if(epmax0 .ne.epmax )write(ifmt,'(a)')'iniev: wrong epmax'
1721 if(qcdlam0.ne.qcdlam.or.q2min0.ne.q2min.or.q2ini0.ne.q2ini
1722 * .or.naflav0.ne.naflav.or.epmax0.ne.epmax)then
1723 write(6,'(//a//)')' iniev has to be reinitialized!!!'
1724 stop
1725 endif
1726 read (1,*)evk0,evk
1727 close(1)
1728 endif
1729 goto 101
1730 endif
1731
1732 write(ifmt,'(a)')'iniev does not exist -> calculate tables ...'
1733 xmax=1.d0-2.d0*q2ini/epmax
1734 do l=1,27
1735 if(l.le.12)then
1736 xx=.1d0*exp(l-13.d0)
1737 elseif(l.le.21)then
1738 xx=.1d0*(l-12.d0)
1739 else
1740 xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
1741 endif
1742
1743 qmin=max(1.d0*q2min,q2ini/(1.-xx))
1744 do i=1,21
1745 qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
1746 do j=1,21
1747 qj=qmin*(qq/qmin)**((j-1)/20.)
1748 if(l.eq.27.or.i.eq.1.or.j.eq.21)then
1749 evk0(i,j,l)=0.
1750 evk0(i,j,l+27)=0.
1751 do k=1,5
1752 evk(i,j,l+27*(k-1))=0.
1753 enddo
1754 else
1755 do k=1,2
1756 evk0(i,j,l+27*(k-1))=log(psev0(qj,qq,xx,k))
1757 enddo
1758 endif
1759 enddo
1760 enddo
1761 enddo
1762
1763 n=1
1764
17651 n=n+1
1766 write(ifmt,2)n
17672 format(5x,i2,'-th order contribution')
1768
1769 do l=1,26
1770 write(ifmt,*)'l',l
1771 if(l.le.12)then
1772 xx=.1d0*exp(l-13.d0)
1773 elseif(l.le.21)then
1774 xx=.1d0*(l-12.d0)
1775 else
1776 xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
1777 endif
1778
1779 qmin=max(1.d0*q2min,q2ini/(1.d0-xx))
1780 do i=2,21
1781 qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
1782 do j=1,20
1783 qj=qmin*(qq/qmin)**((j-1)/20.)
1784 do m=1,3
1785 do k=1,2
1786 if(m.ne.3)then
1787 ev=psev(qj,qq,xx,m,k,n)
1788 ev0=psevi0(qj,qq,xx,m,k)
1789 evs(i,j,l+27*(m-1)+54*(k-1))=log((ev+ev0)/psfap(xx,m-1,k-1)
1790 * /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
1791 elseif(k.ne.1)then
1792 evs(i,j,l+108)=log((psev(qj,qq,xx,m,k,n)+
1793 * psevi0(qj,qq,xx,2,2))/psfap(xx,2,2)
1794 * /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
1795 endif
1796 enddo
1797 enddo
1798 enddo
1799 enddo
1800 enddo
1801
1802 jec=0
1803 do i=2,21
1804 do j=1,20
1805 do l=1,26
1806 do k=1,5
1807 if(n.eq.2.or.evs(i,j,l+27*(k-1)).ne.0..and.
1808 * abs(1.-evk(i,j,l+27*(k-1))/evs(i,j,l+27*(k-1))).gt.1.e-2)then
1809 jec=1
1810 evk(i,j,l+27*(k-1))=evs(i,j,l+27*(k-1))
1811 endif
1812 enddo
1813 enddo
1814 enddo
1815 enddo
1816
1817 if(jec.ne.0)goto 1
1818
1819 write(ifmt,'(a)')'write to iniev ...'
1820 open(1,file=fnie(1:nfnie),status='unknown')
1821 write (1,*)qcdlam,q2min,q2ini,naflav,epmax
1822 write (1,*)evk0,evk
1823 close(1)
1824
1825101 continue
1826 return
1827 end
1828
1829c------------------------------------------------------------------------
1830 function psdbom(s,t,u,qq,long)
1831c-----------------------------------------------------------------------
1832c psdbom - integrand for DIS c-quark cross-sections (matrix element squared)
1833c s - total c.m. energy squared for the scattering (for n=2: s+qq),
1834c t - invariant variable for the scattering |(p1-p3)**2|
1835c u - invariant variable for the scattering |(p1-p4)**2|
1836c qq - photon virtuality
1837c long: 0 - contr. to (F2-F_L), 1 - contr. to F_L
1838c-----------------------------------------------------------------------
1839 include 'epos.incsem'
1840 if(long.eq.0)then !F2-F_L
1841 psdbom=(2.*(t/u+u/t)*(qq**2+(s-qq)**2)/s**2+
1842 * 4.*(qcmass*s/t/u)**2*(qq-2.*qcmass**2)+
1843 * 8.*qcmass**2/t/u*(s-2.*qq)) *2. !=4.5/2.25
1844 else !F_L_C
1845 psdbom=16.*qq*((s-qq)/s**2-qcmass**2/t/u) *2. !=4.5/2.25
1846 endif
1847 return
1848 end
1849
1850c------------------------------------------------------------------------
1851 function psdbin(q1,qq,s,m1,long)
1852c-----------------------------------------------------------------------
1853c psdbin - DIS born cross-section
1854c q1 - virtuality cutoff for current end of the ladder
1855c qq - photon virtuality
1856c s=2(pq) - s_true + qq
1857c s2min - mass cutoff for born scattering
1858c m1 - incoming parton type (0 - g, 1,2 - q)
1859c-----------------------------------------------------------------------
1860 double precision xx
1861 include 'epos.incsem'
1862 include 'epos.inc'
1863
1864 psdbin=0.
1865 q2mass=qcmass**2
1866 s2min=4.*max(q1,q2mass)+qq
1867 if(m1.eq.0.and.s.gt.s2min.and.(idisco.eq.0.or.idisco.eq.2))then
1868 tmax=s/2.
1869 qtq=4.*max(q2mass,q1)/(s-qq)
1870 if(qtq.lt.1.)then
1871 tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1872 else
1873 tmin=.5*s
1874 endif
1875 psdbin=psdbin+psdbor(q1,qq,s,long)*(1./tmin-1./tmax)
1876 endif
1877
1878 if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
1879 *.and.(idisco.eq.0.or.idisco.eq.1))then
1880 m=min(1,iabs(m1))+1
1881 xx=qq/s
1882 psdbin=psdbin+psevi0(q1,qq,xx,m,2)*4.*pi**2*alfe/s
1883 endif
1884 return
1885 end
1886
1887c------------------------------------------------------------------------
1888 function psdbor(q1,qq,s,long)
1889c-----------------------------------------------------------------------
1890c psdbor - DIS born cross-section
1891c q1 - virtuality cutoff for current end of the ladder
1892c qq - photon virtuality
1893c s=2(pq) - s_true + qq
1894c s2min - mass cutoff for born scattering
1895c-----------------------------------------------------------------------
1896 common /ar3/ x1(7),a1(7)
1897 include 'epos.inc'
1898 include 'epos.incsem'
1899 double precision psuds
1900
1901 psdbor=0.
1902 q2mass=qcmass**2
1903 qtq=4.*max(q2mass,q1)/(s-qq)
1904 j=0 !Gluon
1905
1906 tmax=s/2.
1907 if(qtq.lt.1.)then
1908 tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1909 else
1910 tmin=.5*s
1911 endif
1912 if(tmax.lt.tmin.and.ish.ge.1)write(ifmt,*)'s,q1,qq,tmin,tmax',
1913 *s,q1,qq,tmin,tmax
1914
1915 ft=0.
1916 do i=1,7
1917 do m=1,2
1918 t=2.*tmin/(1.+tmin/tmax+(2*m-3)*x1(i)*(1.-tmin/tmax))
1919 u=s-t
1920
1921 qt=t*u/s*(1.-qq/s)
1922 if(qt.lt..999*max(q2mass,q1).and.ish.ge.1)
1923 & write(ifmt,*)'psdbor:qt,q1',qt,q1
1924 fb=psdbom(s,t,u,qq,long)*t**2
1925 ft=ft+a1(i)*fb*pssalf(qt/qcdlam)*sngl(psuds(qt,j))
1926 enddo
1927 enddo
1928 psdbor=ft/s**2*pi**2*alfe/sngl(psuds(q1,j))
1929 return
1930 end
1931
1932c------------------------------------------------------------------------
1933 subroutine psdint(s,qq,sds,sdn,sdb,sdt,sdr,m1,long)
1934c-----------------------------------------------------------------------
1935c psdint - dis cross-sections interpolation - for minimal
1936c effective momentum cutoff in the ladder
1937c s - total c.m. energy squared for the ladder,
1938c qq - photon virtuality,
1939c sds - dis singlet cross-section,
1940c sdn - dis nonsinglet cross-section,
1941c sdb - dis born cross-section,
1942c sdt - dis singlet+resolved cross-section,
1943c m1 - parton type at current end of the ladder (0 - g, 1,2 - q)
1944c-----------------------------------------------------------------------
1945 double precision xx
1946 dimension wk(3),wj(3)
1947 common /psar2/ edmax,epmax
1948 common /psar27/ csds(21,26,4),csdt(21,26,2),csdr(21,26,2)
1949 include 'epos.incsem'
1950 include 'epos.inc'
1951
1952 sds=0.
1953 sdn=0.
1954 sdt=0.
1955 sdr=0.
1956 sdb=psdbin(q2min,qq,s,m1,long)
1957
1958 m=min(1,iabs(m1))+1
1959 qlj=log(qq/q2min)*2.+1.
1960 j=int(qlj)
1961 if(j.lt.1)j=1
1962 if(j.gt.19)j=19
1963 wj(2)=qlj-j
1964 wj(3)=wj(2)*(wj(2)-1.)*.5
1965 wj(1)=1.-wj(2)+wj(3)
1966 wj(2)=wj(2)-2.*wj(3)
1967
1968 s2min=4.*max(q2min,qcmass**2)+qq
1969 if(m1.ne.0)s2min=s2min/(1.-4.*q2ini/(s2min-qq))
1970 if(s.le.s2min.or.idisco.ne.0.and.idisco.ne.2)goto 1
1971
1972 qtq=4.*max(q2min,qcmass**2)/(s-qq)
1973 if(qtq.lt.1.)then
1974 tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
1975 else
1976 tmin=.5*s
1977 endif
1978 tmax=s/2.
1979
1980 sl=log(s/s2min)/log(edmax/s2min)*25.+1.
1981 k=int(sl)
1982 if(k.lt.1)k=1
1983 if(k.gt.24)k=24
1984 wk(2)=sl-k
1985 wk(3)=wk(2)*(wk(2)-1.)*.5
1986 wk(1)=1.-wk(2)+wk(3)
1987 wk(2)=wk(2)-2.*wk(3)
1988
1989 do k1=1,3
1990 k2=k+k1-1
1991 do j1=1,3
1992 sds=sds+csds(j+j1-1,k2,m+2*long)*wj(j1)*wk(k1)
1993 enddo
1994 enddo
1995 if(m.eq.1)then
1996 sds=exp(sds)*(1./tmin-1./tmax)
1997 else
1998 sds=max(sds,0.)
1999 endif
2000
20011 continue
2002 s2min=max(4.*qq,16.*q2min)+qq
2003 if(s.le.s2min.or.long.ne.0.or.idisco.ne.0.and.idisco.ne.3)then
2004 sdt=sds
2005 goto 2
2006 endif
2007
2008 sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2009 k=int(sl)
2010 if(k.lt.1)k=1
2011 if(k.gt.24)k=24
2012 wk(2)=sl-k
2013 wk(3)=wk(2)*(wk(2)-1.)*.5
2014 wk(1)=1.-wk(2)+wk(3)
2015 wk(2)=wk(2)-2.*wk(3)
2016
2017 do k1=1,3
2018 k2=k+k1-1
2019 do j1=1,3
2020 sdr=sdr+csdr(j+j1-1,k2,m)*wj(j1)*wk(k1)
2021 sdt=sdt+csdt(j+j1-1,k2,m)*wj(j1)*wk(k1)
2022 enddo
2023 enddo
2024
2025 sdr=max(sdr,0.)
2026 sdt=max(sds,sds+sdt)
2027 sdt=sdt+sdr
2028
20292 continue
2030 if(long.eq.0.and.q2min.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
2031 *.and.(idisco.eq.0.or.idisco.eq.1))then
2032 xx=qq/s
2033 dsi=psevi(q2min,qq,xx,m,2)*4.*pi**2*alfe/s
2034 if(m1.eq.0)then
2035 sds=sds+dsi
2036 sdt=sdt+dsi
2037 else
2038 dnsi=psevi(q2min,qq,xx,3,2)*4.*pi**2*alfe/s
2039 sdn=sdn+dnsi
2040 sds=sds+max(dsi-dnsi,0.)
2041 sdt=sdt+max(dsi-dnsi,0.)
2042 endif
2043 endif
2044
2045 if(m1.eq.0)then
2046 sds=max(sds,sdb)
2047 sdt=max(sdt,sdb)
2048 else
2049 sdn=max(sdn,sdb)
2050 endif
2051 return
2052 end
2053
2054c-----------------------------------------------------------------------
2055 function psdnsi(q1,qq,s,long)
2056c-----------------------------------------------------------------------
2057c psdnsi - DIS nonsinglet cross-section interpolation
2058c q1 - effective momentum cutoff for current end of the ladder,
2059c qq - photon virtuality,
2060c s - total c.m. energy squared for the ladder,
2061c-----------------------------------------------------------------------
2062 double precision xx
2063 include 'epos.incsem'
2064 include 'epos.inc'
2065
2066 psdnsi=0.
2067 if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2068 xx=qq/s
2069 psdnsi=psdnsi+max(0.,psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s)
2070 endif
2071 return
2072 end
2073
2074c-----------------------------------------------------------------------
2075 function psdrga(qq,s,s2min,j)
2076c-----------------------------------------------------------------------
2077c psdrga - DIS resolved cross-section (photon sf)
2078c qq - photon virtuality
2079c s - total c.m. energy squared for the process
2080c s2min - mass cutoff for born scattering
2081c j - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2082c-----------------------------------------------------------------------
2083 common /ar3/ x1(7),a1(7)
2084 include 'epos.incsem'
2085
2086 psdrga=0.
2087 if(s.le.s2min)return
2088
2089 xmin=s2min/s
2090 do i=1,7
2091 do m=1,2
2092 z=xmin**(.5+(m-1.5)*x1(i))
2093 tu=psdfh4(z,q2min,qq,0,1)
2094 td=psdfh4(z,q2min,qq,0,2)
2095 ts=psdfh4(z,q2min,qq,0,3)
2096 tg=psdfh4(z,q2min,qq,0,0)
2097 if(j.eq.0)then
2098 sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2099 * (tu+td+ts)*psjti(q2min,qq,z*s,1,j,1)
2100 else
2101 sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2102 * (tu+td)*(psjti(q2min,qq,z*s,1,1,1)/4.+
2103 * psjti(q2min,qq,z*s,-1,1,1)/4.+
2104 * psjti(q2min,qq,z*s,2,1,1)/2.)+
2105 * ts*psjti(q2min,qq,z*s,2,1,1)
2106 endif
2107 psdrga=psdrga+a1(i)*sj
2108 enddo
2109 enddo
2110 psdrga=-psdrga*log(xmin)*alfe/2. *4.5 !mean e^2 is taken out
2111 return
2112 end
2113
2114c-----------------------------------------------------------------------
2115 function psdres(qq,s,s2min,j)
2116c-----------------------------------------------------------------------
2117c psdres - DIS resolved photon cross-section
2118c qq - photon virtuality
2119c s - total w squared for the ladder (s+qq)
2120c s2min - mass cutoff for born scattering
2121c j - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2122c-----------------------------------------------------------------------
2123 double precision xx
2124 common /ar3/ x1(7),a1(7)
2125 include 'epos.inc'
2126 include 'epos.incsem'
2127
2128 psdres=0.
2129 if(s.le.s2min+qq)return
2130
2131 qmin=max(q2min,s2min/(s/qq-1.))
2132 qmax=min(s-s2min,s/2.)
2133
2134c numerical integration over transverse momentum squared;
2135c gaussian integration is used
2136 do i=1,7
2137 do m=1,2
2138 qi=2.*qmin/(1.+qmin/qmax+(2*m-3)*x1(i)*(1.-qmin/qmax))
2139
2140 zmax=min(1.,qi/qq)
2141 zmin=(max(qi,s2min)+qi)/s
2142
2143 fsj=0.
2144 if(zmax.gt.zmin)then
2145 do i1=1,7
2146 do m1=1,2
2147 z=.5*(zmax+zmin+(2*m1-3)*x1(i1)*(zmax-zmin))
2148 s2=z*s-qi
2149 xx=z
2150 if(j.eq.0)then
2151 sj=psfap(xx,0,1)*psjti(qi,qq,s2,1,j,1)
2152 else
2153 sj=psfap(xx,0,1)*(psjti(qi,qq,s2,1,1,1)/6.+
2154 * psjti(qi,qq,s2,-1,1,1)/6.+
2155 * psjti(qi,qq,s2,2,1,1)/1.5)
2156 endif
2157 fsj=fsj+a1(i1)*sj*qi !????????(qi-z*qq)
2158 enddo
2159 enddo
2160 fsj=fsj*(zmax-zmin)
2161 elseif(ish.ge.1)then
2162 write(ifmt,*)'psdres:zmax,zmin',zmax,zmin
2163 endif
2164 psdres=psdres+a1(i)*fsj
2165 enddo
2166 enddo
2167 psdres=psdres*(1./qmin-1./qmax)*alfe*.75/pi !alpha_s -> 6 alpha_e
2168 return
2169 end
2170
2171c------------------------------------------------------------------------
2172 function psds(q1,qq,s,j,long)
2173c-----------------------------------------------------------------------
2174c psds - DIS singlet cross-section
2175c q1 - virtuality cutoff for current end of the ladder
2176c qq - photon virtuality
2177c s=2(pq) - s_true + qq
2178c s2min - mass cutoff for born scattering
2179c-----------------------------------------------------------------------
2180 double precision xxe,xmax,xmin,xmax1,xmin1
2181 common /ar3/ x1(7),a1(7)
2182 include 'epos.inc'
2183 include 'epos.incsem'
2184
2185 psds=0.
2186 q2mass=qcmass**2
2187 s2min=4.*max(q1,q2mass)
2188 smin=(s2min+qq)/(1.-4.*q2ini/s2min)
2189 if(s.le.1.001*smin)return
2190
2191 xmax=.5d0*(1.d0+qq/s+dsqrt((1.d0-qq/s)**2-16.d0*q2ini/s))
2192 xmin=max(1.d0+qq/s-xmax,1.d0*(s2min+qq)/s)
2193 if(xmin.gt.xmax.and.ish.ge.1)write(ifmt,*)'xmin,xmax,q1,qq,s,smin'
2194 *,xmin,xmax,q1,qq,s,smin
2195
2196 fx1=0.
2197 fx2=0.
2198 if(xmax.gt..9d0)then
2199 xmin1=max(xmin,.9d0)
2200 do i=1,7
2201 do m=1,2
2202 xxe=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
2203 * (.5d0-x1(i)*(m-1.5))
2204 xx=xxe
2205
2206 sh=xx*s
2207 qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2208 qtq=4.*qtmin/(sh-qq)
2209
2210 tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2211 tmax=.5*sh
2212 if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2213 & ,tmin,tmax
2214
2215 ft=0.
2216 do i1=1,7
2217 do m1=1,2
2218 t=.5*(tmin+tmax+(2*m1-3)*x1(i1)*(tmin-tmax))
2219 u=sh-t
2220 qt=t*u/sh*(1.-qq/sh)
2221 if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2222 & ,qt,qtmin
2223
2224 fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2225 ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2226 enddo
2227 enddo
2228 ft=ft*(tmax-tmin)
2229 fx1=fx1+a1(i)*ft*(1.-xx)/sh**2
2230 enddo
2231 enddo
2232 fx1=fx1*log((1.d0-xmin1)/(1.d0-xmax))
2233 endif
2234
2235 if(xmin.lt..9d0)then
2236 xmax1=min(xmax,.9d0)
2237 do i=1,7
2238 do m=1,2
2239 xxe=xmin*(xmax1/xmin)**(.5-x1(i)*(m-1.5))
2240 xx=xxe
2241
2242 sh=xx*s
2243 qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2244 qtq=4.*qtmin/(sh-qq)
2245
2246 tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2247 tmax=.5*sh
2248 if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2249 * ,tmin,tmax
2250
2251 ft=0.
2252 do i1=1,7
2253 do m1=1,2
2254 t=(.5*(tmin+tmax+(2*m1-3)*x1(i1)*
2255 * (tmin-tmax)))
2256 u=sh-t
2257 qt=t*u/sh*(1.-qq/sh)
2258 if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2259 * ,qt,qtmin
2260
2261 fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2262 ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2263 enddo
2264 enddo
2265 ft=ft*(tmax-tmin)
2266 fx2=fx2+a1(i)*ft*xx/sh**2
2267 enddo
2268 enddo
2269 fx2=fx2*log(xmax1/xmin)
2270 endif
2271 psds=(fx1+fx2)*pi**2*alfe
2272 return
2273 end
2274
2275c-----------------------------------------------------------------------
2276 function psdsj(q1,xx,s,qt,t,u,qq,j,long)
2277c-----------------------------------------------------------------------
2278c psdsj - integrand for dis singlet cross-section
2279c q1 - virtuality cutoff for current end of the ladder
2280c xx - lc momentum ratio between initial (j) and final (l) partons
2281c s - c.m. energy squared for the born scattering,
2282c t - invariant variable for the born scattering |(p1-p3)**2|
2283c u - invariant variable for the born scattering |(p1-p4)**2|
2284c qq - photon virtuality
2285c j - initial parton at the end of the ladder (0 - g, 1,2 - q)
2286c-----------------------------------------------------------------------
2287 double precision xx
2288 include 'epos.incsem'
2289
2290 fb=psdbom(s,t,u,qq,long)
2291 psdsj=psevi(q1,qt,xx,min(1,iabs(j))+1,1)*fb
2292 return
2293 end
2294
2295c------------------------------------------------------------------------
2296 function psdsin(q1,qq,s,m1,long)
2297c-----------------------------------------------------------------------
2298c psdsin - DIS singlet cross-section interpolation
2299c q1 - effective momentum cutoff for current end of the ladder,
2300c qq - photon virtuality,
2301c s - total c.m. energy squared for the ladder,
2302c m1 - parton type at current end of the ladder (0 - g, 1,2 - q)
2303c-----------------------------------------------------------------------
2304 double precision xx
2305 dimension wi(3),wj(3),wk(3)
2306 common /psar2/ edmax,epmax
2307 common /psar25/ csdsi(21,21,104)
2308 include 'epos.incsem'
2309 include 'epos.inc'
2310
2311 psdsin=0.
2312 m=min(1,iabs(m1))+1
2313
2314 q2mass=qcmass**2
2315 s2min=4.*max(q2min,q2mass)+qq
2316 sdmin=4.*max(q1,q2mass)+qq
2317 if(m1.ne.0)then
2318 s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2319 sdmin=sdmin/(1.-4.*q2ini/(sdmin-qq))
2320 endif
2321c if(s.le.1.e8*sdmin)goto 2 !????????????????
2322 if(s.le.sdmin)goto 2
2323
2324 qmin=q2min
2325 qmax=(s-qq)/4.
2326 if(m1.ne.0)qmax=(s-qq+sqrt((s-qq)**2-16.*s*q2ini))/8.
2327 qtq=4.*max(q2mass,q1)/(s-qq)
2328 if(qtq.lt.1.)then
2329 tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2330 else
2331 tmin=.5*s
2332 endif
2333 tmax=s/2.
2334
2335 qlj=log(qq/q2min)*2.+1.
2336 j=int(qlj)
2337 if(j.lt.1)j=1
2338 if(j.gt.19)j=19
2339 wj(2)=qlj-j
2340 wj(3)=wj(2)*(wj(2)-1.)*.5
2341 wj(1)=1.-wj(2)+wj(3)
2342 wj(2)=wj(2)-2.*wj(3)
2343
2344 qli=log(q1/qmin)/log(qmax/qmin)*20.+1.
2345 i=int(qli)
2346 if(i.lt.1)i=1
2347 if(i.gt.19)i=19
2348 wi(2)=qli-i
2349 wi(3)=wi(2)*(wi(2)-1.)*.5
2350 wi(1)=1.-wi(2)+wi(3)
2351 wi(2)=wi(2)-2.*wi(3)
2352
2353 sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2354 k=int(sl)
2355 if(k.lt.1)k=1
2356 if(k.gt.24)k=24
2357 wk(2)=sl-k
2358 wk(3)=wk(2)*(wk(2)-1.)*.5
2359 wk(1)=1.-wk(2)+wk(3)
2360 wk(2)=wk(2)-2.*wk(3)
2361
2362 dsin1=0.
2363 do k1=1,3
2364 k2=k+k1-1+26*(m-1)+52*long
2365 do i1=1,3
2366 do j1=1,3
2367 dsin1=dsin1+csdsi(i+i1-1,j+j1-1,k2)*wi(i1)*wj(j1)*wk(k1)
2368 enddo
2369 enddo
2370 enddo
2371 if(m1.eq.0)then
2372 psdsin=psdsin+exp(dsin1)*(1./tmin-1./tmax)
2373 else
2374 psdsin=psdsin+max(0.,dsin1)
2375 endif
2376
23772 continue
2378 if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2379 xx=qq/s
2380 dsi=psevi(q1,qq,xx,m,2)*4.*pi**2*alfe/s
2381 if(m1.eq.0)then
2382 psdsin=psdsin+max(dsi,0.)
2383 else
2384 dnsi=psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s
2385 psdsin=psdsin+max(dsi-dnsi,0.)
2386 endif
2387 endif
2388 return
2389 end
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404c###########################################################################
2405c###########################################################################
2406c###########################################################################
2407c###########################################################################
2408c
2409c unused 3P
2410c
2411c###########################################################################
2412c###########################################################################
2413c###########################################################################
2414c###########################################################################
2415
2416
2417c------------------------------------------------------------------------
2418 function psvy(xpp0,xpr0,xpm0,xmr0,b,iqq)
2419c-----------------------------------------------------------------------
2420c psvy - 3p-contributions to the interaction eikonal
2421c xpp - lc+ for the pomeron,
2422c xpr - lc+ for the remnant,
2423c xpm - lc- for the pomeron,
2424c xpr - lc- for the remnant,
2425c b - impact parameter,
2426c iqq=1 - Y-proj-uncut
2427c iqq=2 - Y-proj-1-cut
2428c iqq=3 - Y-proj-2-cut
2429c iqq=4 - Y-proj-soft-cut
2430c iqq=5 - Y-proj-gss-cut
2431c iqq=6 - Y-proj-qss-cut
2432c iqq=7 - Y-proj-ssg-cut
2433c iqq=8 - Y-proj-ssq-cut
2434c iqq=9 - Y-proj-difr
2435c iqq=-1 - Y-targ-uncut
2436c iqq=-2 - Y-targ-1-cut
2437c iqq=-3 - Y-targ-2-cut
2438c iqq=-4 - Y-targ-soft-cut
2439c iqq=-5 - Y-targ-gss-cut
2440c iqq=-6 - Y-targ-qss-cut
2441c iqq=-7 - Y-targ-ssg-cut
2442c iqq=-8 - Y-targ-ssq-cut
2443c iqq=-9 - Y-targ-difr
2444c------------------------------------------------------------------------
2445
2446 psvy=0.
2447 return
2448 end
2449
2450c------------------------------------------------------------------------
2451 function psvx(xpp,xpr,xpm,xmr,b,iqq)
2452c-----------------------------------------------------------------------
2453c psvx - 4p-contributions to the interaction eikonal
2454c xpp - lc+ for the pomeron,
2455c xpr - lc+ for the remnant,
2456c xpm - lc- for the pomeron,
2457c xpr - lc- for the remnant,
2458c b - impact parameter,
2459c iqq=0 - X-uncut
2460c iqq=1 - X-1-cut
2461c iqq=2 - X-Y+cut
2462c iqq=-2 - X-Y-cut
2463c iqq=3 - X-1-cut-soft
2464c iqq=4 - X-1-cut-gss
2465c iqq=-4 - X-1-cut-gss
2466c iqq=5 - X-1-cut-qss
2467c iqq=-5 - X-1-cut-qss
2468c iqq=6 - X-difr+
2469c iqq=-6 - X-difr-
2470c------------------------------------------------------------------------
2471
2472 psvx=0.
2473 return
2474 end
2475
2476
2477
2478
2479
2480c
2481c
2482c
2483c
2484c
2485cc------------------------------------------------------------------------
2486c function psftig(x,xpomr,jj)
2487cc-----------------------------------------------------------------------
2488c
2489c psftig=0.
2490c end
2491c
2492c
2493cc------------------------------------------------------------------------
2494c function psftih(zz,ddd)
2495cc-----------------------------------------------------------------------
2496c
2497c psftih=0.
2498c return
2499c end
2500c
2501cc------------------------------------------------------------------------
2502c function psftij(zz,ddd)
2503cc------------------------------------------------------------------------
2504c
2505c psftij=0.
2506c return
2507c end
2508c
2509cc------------------------------------------------------------------------
2510c function psftik(xp,del1,rh1,rh2,rh3,rp,z)
2511cc-----------------------------------------------------------------------
2512c
2513c psftik=0.
2514c return
2515c end
2516c
2517cc------------------------------------------------------------------------
2518c function psftim(xp1,xp,del1,rh1,rh2,rh3,rp,z)
2519cc-----------------------------------------------------------------------
2520c
2521c psftim=0.
2522c return
2523c end
2524c
2525cc------------------------------------------------------------------------
2526c function psftil(xp,del1,rh1,rh2,rh3,rp,z)
2527cc------------------------------------------------------------------------
2528c
2529c psftil=0.
2530c return
2531c end
2532c
2533cc------------------------------------------------------------------------
2534c function psftigt(x)
2535cc-----------------------------------------------------------------------
2536c
2537c psftigt=0.
2538c return
2539c end
2540c
2541cc------------------------------------------------------------------------
2542c function psftist(x)
2543cc-----------------------------------------------------------------------
2544c psftist=0.
2545c return
2546c end
2547c
2548cc------------------------------------------------------------------------
2549c function psftig1(x,xpomr,jj)
2550cc-----------------------------------------------------------------------
2551c psftig1=0.
2552c return
2553c end
2554c
2555cc------------------------------------------------------------------------
2556c function psftis1(x,xpomr,jj)
2557cc-----------------------------------------------------------------------
2558c psftis1=0.
2559c return
2560c end
2561c
2562cc------------------------------------------------------------------------
2563c function psd3p1(xd,xpomr,qq,jj)
2564cc-----------------------------------------------------------------------
2565cc psd3p1 - df2difr/dx_pomr
2566cc xd - bjorken x,
2567cc xpomr - pomeron x,
2568cc qq - photon virtuality
2569cc jj=1 - 1st order
2570cc-----------------------------------------------------------------------
2571c
2572c psd3p1=0.
2573c
2574c return
2575c end
2576c
2577cc------------------------------------------------------------------------
2578c function psv3p(sy,xpp,xpm,zb)
2579cc-----------------------------------------------------------------------
2580cc psv3p - 3p-contributions to the interaction eikonal
2581cc sy - energy squared for the hard interaction,
2582cc xpp - lc+ for the sh pomeron,
2583cc xpm - lc- for the sh pomeron,
2584cc z - impact parameter factor, z=exp(-b**2/4*rp),
2585cc------------------------------------------------------------------------
2586c
2587c psv3p=0.
2588c return
2589c end
2590c
2591cc------------------------------------------------------------------------
2592c function psfro(xpomr,zb,iclp,icdpro)
2593cc-----------------------------------------------------------------------
2594cc psfro - generalized froissaron between proj. and 3p-vertex
2595cc xpp - lc+ for the proj. side,
2596cc xpomr - lc+ for the vertex,
2597cc zb - impact parameter factor, z=exp(-b**2/4*rp),
2598cc------------------------------------------------------------------------
2599c psfro=0.
2600c return
2601c end
2602c
2603cc------------------------------------------------------------------------
2604c function psv2(xpomr,zb,iclp,iqq)
2605cc-----------------------------------------------------------------------
2606cc psv2 - 2-pom contribution to the froissaron
2607cc xpomr - lc+ for the vertex,
2608cc zb - impact parameter factor, z=exp(-b**2/4*rp),
2609cc------------------------------------------------------------------------
2610c
2611c psv2=0.
2612c return
2613c end
2614c
2615cc------------------------------------------------------------------------
2616c function psvfro(xpp,xpr,xpomr,zb,iclp,icdpro,iqq)
2617cc-----------------------------------------------------------------------
2618cc psvfro - effective froissaron contributions
2619cc xpomr - lc+ for the vertex,
2620cc zb - impact parameter factor, z=exp(-b**2/4*rp),
2621cc iqq=0 - total uncut
2622cc iqq=1 - total 1-cut
2623cc iqq=2 - total 2-cut
2624cc iqq=3 - soft 1-cut
2625cc iqq=4 - gg 1-cut
2626cc iqq=5 - qg 1-cut
2627cc iqq=6 - difr
2628cc------------------------------------------------------------------------
2629c
2630c psvfro=0.
2631c return
2632c end
2633c
2634cc------------------------------------------------------------------------
2635c function psfroin(xpomr,z,iclp,icdpro)
2636cc-----------------------------------------------------------------------
2637cc psfroin - interpolation of effective froissaron corrections
2638cc xpomr - lc+ for the 3p-vertex,
2639cc z - impact parameter factor, z=exp(-b**2/4*rp),
2640cc-----------------------------------------------------------------------
2641c
2642c psfroin=0.
2643c return
2644c end
2645c
2646cc------------------------------------------------------------------------
2647c function psvnorm(b)
2648cc-----------------------------------------------------------------------
2649cc psvnorm - X-contribution normalization
2650cc b - impact parameter
2651cc-----------------------------------------------------------------------
2652c
2653c psvnorm=0.
2654c return
2655c end
2656c
2657cc------------------------------------------------------------------------
2658c function psvxb(dxpp,dxpr,dxpm,dxmr,dxpomr,bb1,bb2
2659c *,iclp,iclt,icdpro,icdtar,iqq)
2660cc-----------------------------------------------------------------------
2661cc psvxb - integrand for X-contributions
2662cc dxpomr - lc+ for the vertex,
2663cc bb1,bb2 - impact parameters to proj(targ),
2664cc iqq=0 - uncut
2665cc iqq=1 - 1-cut
2666cc iqq=2 - Y-cut
2667cc iqq=3 - soft 1-cut
2668cc iqq=4 - gg 1-cut
2669cc iqq=5 - qg 1-cut
2670cc iqq=6 - difr
2671cc------------------------------------------------------------------------
2672c double precision dxpp,dxpr,dxpm,dxmr,dxpomr
2673c
2674c psvxb=0.
2675c return
2676c end
2677c
2678cc------------------------------------------------------------------------
2679c function pscoef(zz,alp1,alp2,iclp,iqq)
2680cc-----------------------------------------------------------------------
2681cc pscoef - integrated vertexes
2682cc zz=xpomr/xpr(xpp)
2683cc iqq=0 - 2-cut
2684cc iqq=1 - 1-uncut
2685cc iqq=2 - 2-uncut
2686cc------------------------------------------------------------------------
2687c
2688c pscoef=0.
2689c return
2690c end
2691c
2692cc------------------------------------------------------------------------
2693c function pscoefi(z,i1,i2,iclp,iqq)
2694cc-----------------------------------------------------------------------
2695cc pscoefi - interpolation of integrated vertexes
2696cc z=xpomr/xpr(xpp)
2697cc iqq=0 - 2-cut
2698cc iqq=1 - 1-uncut
2699cc iqq=2 - 2-uncut
2700cc------------------------------------------------------------------------
2701c pscoefi=0.
2702c return
2703c end
2704c
2705c