New generator: EPOS
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-rsh-165.f
CommitLineData
9ef1c2d9 1c reshuffled from sem, sto, sha
2
3c contains psahot and related stuff
4c ------
5
6
7
8c-----------------------------------------------------------------------
9 subroutine psahot(kcol,ncolp,iret)
10c-----------------------------------------------------------------------
11c psahot - showering (semihard parton-parton interaction)
12c-----------------------------------------------------------------------
13 include 'epos.inc'
14 include 'epos.incsem'
15 include 'epos.incems'
16 double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,wplc
17 *,wp0,wm0,s,si,smin,xmin,xp1,wpi,wmi,xp,xm,wp1,wm1,wp2,wm2
18 *,wpq,wmq,wpi1,wmi1,pxh,pyh,pth,xmm,xp2,xg,zg,smax,xmax
19 *,xmax1,xmin1,zp0,psutz,xpmax0,xpmin0,gb0,tmax0,tmin0,zpm,gb
20 *,gbyj,tmax,tmin,t,gb7,x,s2,discr,qt,qt2,x1min,x1max,t1min,t1max
21 *,t1,xq1,qq,qqmax,qqmin,pt2,xmin2,ptnew,xpmin,xpmax,xm0,psuds
22c *,xprh,xmrh
23 dimension ep3(4),ey(3),bx(6)
24 *,qmin(2),iqc(2),nqc(2),ncc(2,2),amqt(4)
25 parameter (mjstr=20000)
26 common /psar7/ delx,alam3p,gam3p
27 common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
28 common /psar30/iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
29 common /testj/ ajeth(4),ajete(5),ajet0(7)
30 parameter (ntim=1000)
31 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
32 &,idaprt(2,ntim)
33 common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
34 integer icp(2),ict(2),icini(2)
35 integer jcp(nflav,2),jct(nflav,2),nemis(2)
36 common/cprtx/nprtjx,pprtx(5,2)/ciptl/iptl
37
38 call utpri('psahot',ish,ishini,3)
39
40 iret=0
41 alpq=-(alppar+1.)/2.
42 qqcut=q2min !????????????pt2cut
43
44
45 nptl1=nptl
46 iptl=nppr(ncolp,kcol)
47 ip=iproj(kcol)
48 it=itarg(kcol)
49 do i=1,2
50 icp(i)=icproj(i,ip)
51 ict(i)=ictarg(i,it)
52 enddo
53 idpomr=idhpr(ncolp,kcol)
54 bpomr=bhpr(ncolp,kcol)
55
56 q2finsave=q2fin
57 zzzz=zzremn(ip,1)+zzremn(it,2)
58 zz=min( 1+zoemax , 1+zoeinc *zzzz ) !<-----
59 q2fin=q2fin*zz
60
61 ajeth(idpomr+1)=ajeth(idpomr+1)+1.
62c write(*,*)ajeth
63 idfpomr=idfpr(ncolp,kcol)
64 if(ish.ge.3)write(ifch,*)'Start psahot (icp,ict):',ip,icp,it,ict
65 *,ncolp,kcol,iptl,idpomr,idfpomr,bpomr
66
67
68 if(idfpomr.eq.0)stop'idfpomr??????'
69 if(ish.ge.3)then
70 write(ifch,20)iptl
71 * ,sqrt(pptl(1,iptl)**2+pptl(2,iptl)**2),pptl(3,iptl)
72 * ,pptl(4,iptl),pptl(5,iptl)
73 * ,istptl(iptl),ityptl(iptl)
7420 format(1x,i4,3x,4(e11.5,1x),i2,1x,i3)
75 endif
76 istptl(iptl)=31
77
781 nj=0
79 nptl=nptl1
80
81 wp0=dsqrt(xpr(ncolp,kcol))*dexp(ypr(ncolp,kcol))*dble(engy) !???? new
82 wm0=dsqrt(xpr(ncolp,kcol))*dexp(-ypr(ncolp,kcol))*dble(engy) !double
83
84
85 amqt(1)=sqrt(sngl(xxp1pr(ncolp,kcol)**2+xyp1pr(ncolp,kcol)**2))
86 amqt(2)=sqrt(sngl(xxp2pr(ncolp,kcol)**2+xyp2pr(ncolp,kcol)**2))
87 amqt(3)=sqrt(sngl(xxm2pr(ncolp,kcol)**2+xym2pr(ncolp,kcol)**2))
88 amqt(4)=sqrt(sngl(xxm1pr(ncolp,kcol)**2+xym1pr(ncolp,kcol)**2))
89 amqpt=amqt(1)+amqt(2)+amqt(3)+amqt(4)
90
91 s2min=4.*q2min
92 if(sngl(wp0*wm0).le.(sqrt(s2min)+amqpt)**2)then
93 if(ish.ge.1)then
94 call utmsg('psahot: insufficient pomeron mass&')
95 write (ifch,*)'mass:',dsqrt(wp0*wm0),amqpt+sqrt(s2min)
96 call utmsgf
97 endif
98 iret=1
99 goto 16
100 endif
101
102 ih=iproj(kcol)
103 jh=itarg(kcol)
104c xprh=xpp(ih)
105c xmrh=xmt(jh)
106 rp=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
107 z=exp(-bpomr**2/(4.*.0389*rp))
108
109 if(z.eq.0)then
110 write(ifch,*)'psahot : z,ih,jh ! -> ',z,ih,jh
111 call gakli2(ih,ih)
112 call gakli2(jh,jh)
113 call gakli2(iptl,iptl)
114 stop
115 endif
116
117 do l=1,4
118 bx(l)=xorptl(l,iptl)
119 enddo
120 bx(5)=tivptl(1,iptl)
121 bx(6)=tivptl(2,iptl)
122 ity=ityptl(iptl)
123
124 if(idpomr.eq.0)then !gg-pomeron
125 iqq=0 !type of the hard interaction: 0 - gg, 1 - qg, 2 - gq, 3 - qq
126 pxh=0.d0 !p_x for sh pomeron
127 pyh=0.d0 !p_y for sh pomeron
128 elseif(idpomr.eq.1)then !qg-pomeron
129 iqq=1
130 pxh=xxp1pr(ncolp,kcol)
131 pyh=xyp1pr(ncolp,kcol)
132 amqpt=amqpt-amqt(1)
133 elseif(idpomr.eq.2)then !gq-pomeron
134 iqq=2
135 pxh=xxm2pr(ncolp,kcol)
136 pyh=xym2pr(ncolp,kcol)
137 amqpt=amqpt-amqt(3)
138 elseif(idpomr.eq.3)then !qq-pomeron
139 iqq=3
140 pxh=xxp1pr(ncolp,kcol)+xxm2pr(ncolp,kcol)
141 pyh=xyp1pr(ncolp,kcol)+xym2pr(ncolp,kcol)
142 amqpt=amqpt-amqt(1)-amqt(3)
143 else
144 stop'unknown pomeron'
145 endif
146 pth=pxh**2+pyh**2
147
148 nj0=nj
149 if(ish.ge.6)write(ifch,*)'iptl,nptl,wp0,wm0,z,iqq,bx:'
150 if(ish.ge.6)write(ifch,*) iptl,nptl,wp0,wm0,z,iqq,bx
151
152 s=wp0*wm0 !lc+*lc- for the semihard interaction
153 smin=dble(s2min)+pth !mass cutoff for the hard pomeron
154 xmin=smin/s
155 smax=(dsqrt(s)-dble(amqpt))**2 !max mass for the hard pomeron
156 xmax=smax/s
157 if(smax.le.smin)then
158 write (ifmt,*)'smax,smin',smax,smin
159 iret=1
160 goto 16
161 endif
162 xpmax0=psutz(s,smin,dble(amqpt**2)) !max x+ for the hard pomeron
163 xpmin0=xmin/xpmax0 !min x+ for the hard pomeron
164 xp1=wp0/dble(engy) !lc+ share for the semihard interaction
165 xp2=wm0/dble(engy) !lc- share for the semihard interaction
166
167
168c-----------------------------------------------------------------
169c determine LC momenta wpi,wmi for hard Pomeron
170c-----------------------------------------------------------------
171
172 if(ish.ge.4)write(ifch,*)
173 & 'determine LC momenta wpi,wmi for hard Pomeron'
174
175 !------------------------------------------
176 if(iqq.eq.3)then ! val-val
177 !------------------------------------------
178
179 if(ish.ge.4)write(ifch,*)'val-val'
180 xmin1=xmin**dble(delh+.4)
181 xmax1=xmax**dble(delh+.4)
182 zp0=dsqrt(xmin)
183 if(zp0.ge.1.d0)call utstop('zp0 in sem&')
184 !........ kinematical bounds
185 tmax0=dlog((1.d0+dsqrt(1.d0-zp0))/(1.d0-dsqrt(1.d0-zp0)))
186 tmin0=dlog((1.d0+dsqrt(1.d0-xpmax0))
187 * /(1.d0-dsqrt(1.d0-xpmax0)))
188 if(iclpro.ne.4)then
189 call psjti0(sngl(smax-pth),sqq,sqqb,1,1)
190 call psjti0(sngl(smax-pth),sqqp,sqqpb,1,2)
191 call psjti0(sngl(smax-pth),sqaq,sqaqb,-1,1)
192 else
193 call psjti0(sngl(smax-pth),sqqp,sqqpb,4,2)
194 sqq=0.
195 sqaq=0.
196 endif
197 uv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,1)
198 dv1=psdfh4(sngl(zp0*xp1),q2min,0.,iclpro,2)
199 uv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,1)
200 dv2=psdfh4(sngl(zp0*xp2),q2min,0.,icltar,2)
201 wwuu=uv1*uv2*sqq
202 if(iclpro.eq.2)then
203 wwdd=dv1*dv2*sqq
204 elseif(iclpro.eq.1)then
205 wwdd=dv1*dv2*sqaq
206 elseif(iclpro.eq.3)then
207 wwdd=dv1*dv2*sqqp
208 elseif(iclpro.eq.4)then
209 wwuu=uv1*uv2*sqqp
210 wwdd=0.
211 endif
212 wwud=uv1*dv2*sqqp
213 wwdu=dv1*uv2*sqqp
214 wudt=wwuu+wwdd+wwud+wwdu
215 gb0=dble(wudt)/xmax**dble(delh)/xmin**0.4*
216 * (1.d0-zp0*xp1)**dble(-1.-alpq-alplea(iclpro))*
217 * (1.d0-zp0*xp2)**dble(-1.-alpq-alplea(icltar))*(tmax0-tmin0)
218 * *(1.d0-zp0)**dble(.5+alpq)*(1.d0-zp0)**dble(alpq) *5.d0
219
2203 zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh+.4)) !zpm
221 if(iclpro.ne.4)then
222 call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
223 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
224 call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
225 else
226 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,4,2)
227 sqq=0.
228 sqaq=0.
229 endif
230 xpmax=psutz(s,zpm*s,dble(amqpt**2)) !max x+ for sh pomeron
231 tmax=dlog((1.d0+dsqrt(1.d0-dsqrt(zpm)))
232 * /(1.d0-dsqrt(1.d0-dsqrt(zpm))))
233 tmin=dlog((1.d0+dsqrt(1.d0-xpmax))/(1.d0-dsqrt(1.d0-xpmax)))
234
235 t=(tmin+dble(rangen())*(tmax-tmin))
236 xp=1.d0-((1.d0-dexp(-t))/(1.d0+dexp(-t)))**2 !x+_v
237 xm=zpm/xp !x-_v
238 if(xm.gt.xp.and.ish.ge.1)write(ifmt,*)'xm,xp',xm,xp
239 gb=(1.d0-xm)**alpq*(1.d0-xp)**(.5+alpq)*(tmax-tmin)
240 if(rangen().lt..5)then
241 xp=xm
242 xm=zpm/xp
243 endif
244
245 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
246 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
247 uv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
248 dv2=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
249 wwuu=uv1*uv2*sqq
250 if(iclpro.eq.2)then
251 wwdd=dv1*dv2*sqq
252 elseif(iclpro.eq.1)then
253 wwdd=dv1*dv2*sqaq
254 elseif(iclpro.eq.3)then
255 wwdd=dv1*dv2*sqqp
256 elseif(iclpro.eq.4)then
257 wwuu=uv1*uv2*sqqp
258 wwdd=0.
259 endif
260 wwud=uv1*dv2*sqqp
261 wwdu=dv1*uv2*sqqp
262 wudt=wwuu+wwdd+wwud+wwdu
263
264 gb=gb*dble(wudt)/zpm**dble(delh+0.4)
265 * *(1.d0-xp*xp1)**dble(-1.-alpq-alplea(iclpro))
266 * *(1.d0-xm*xp2)**dble(-1.-alpq-alplea(icltar))/gb0
267c if(ish.ge.4)then
268 if(gb.gt.1.d0.and.ish.ge.1)write (ifch,*)
269 * 'gb-qq,iclpro,zpm,xp,tmax,tmin,xpmax',
270 * gb,iclpro,zpm,xp,tmax,tmin,xpmax
271c endif
272 if(dble(rangen()).gt.gb)goto 3
273
274 aks=rangen()*wudt
275 if(aks.lt.wwuu)then
276 if(iclpro.ne.4)then
277 iq1=1
278 else
279 iq1=4*(2.*int(.5+rangen())-1.)
280 endif
281 iq2=1
282 elseif(aks.lt.wwuu+wwdd)then
283 if(iclpro.eq.2)then
284 iq1=2
285 elseif(iclpro.eq.1)then
286 iq1=-2
287 elseif(iclpro.eq.3)then
288 iq1=-3
289 else
290 stop'?????'
291 endif
292 iq2=2
293 elseif(aks.lt.wwuu+wwdd+wwud)then
294 if(iclpro.ne.4)then
295 iq1=1
296 else
297 iq1=4*(2.*int(.5+rangen())-1.)
298 endif
299 iq2=2
300 else
301 if(iclpro.eq.2)then
302 iq1=2
303 elseif(iclpro.eq.1)then
304 iq1=-2
305 elseif(iclpro.eq.3)then
306 iq1=-3
307 else
308 stop'?????'
309 endif
310 iq2=1
311 endif
312
313 wpi=xp*wp0 !lc+ for the semihard interaction
314 wmi=xm*wm0 !lc- for the semihard interaction
315 wp1=(wp0-wpi)
316 wm1=(wm0-wmi)
317 wp1=wp1*psutz(wp1*wm1,dble(amqt(2)**2),dble(amqt(4)**2))
318 wm1=wm1-amqt(2)**2/wp1
319
320 !-------------------------------------
321 else ! sea-sea val-sea sea-val
322 !-------------------------------------
323
324 if(ish.ge.4)write(ifch,*)'sea-sea val-sea sea-val'
325 xmin1=xmin**(delh-dels)
326 xmax1=xmax**(delh-dels)
327
328 if(iqq.eq.0)then !rejection function normalization
329 gb0=dlog(xpmax0/xpmin0)*(1.d0-dsqrt(xmin))**(2.*betpom) !y_soft =
330 else
331 tmax0=acos(dsqrt(xpmin0)) !kinematical limits for t=cos(x+-)**2
332 tmin0=acos(dsqrt(xpmax0))
333 if(iqq.eq.1)then
334 uv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,1)
335 * *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
336 dv1=psdfh4(sngl(xpmin0*xp1),q2min,0.,iclpro,2)
337 * *sngl(1.d0-xpmin0*xp1)**(-1.-alpq-alplea(iclpro))
338 else
339 uv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,1)
340 * *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
341 dv1=psdfh4(sngl(xpmin0*xp2),q2min,0.,icltar,2)
342 * *sngl(1.d0-xpmin0*xp2)**(-1.-alpq-alplea(icltar))
343 endif
344 gb0=(1.d0-xmin/xpmax0)**betpom*dble(uv1+dv1)
345 * *xpmin0**(-0.5+dels)
346 * *(1.d0-xpmin0)**(0.5+alpq)*(tmax0-tmin0)
347 if(ish.ge.6)write (ifch,*)
348 * 'gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2',
349 * gb0,tmax0,tmin0,xpmax0,xpmin0,xmin,xp1,xp2
350 endif
351 if(iclpro.ne.4.or.iqq.ne.1)then
352 call psjti0(sngl(smax-pth),sj,sjb,iqq,0) !inclusive (sj) and born
353 else
354 call psjti0(sngl(smax-pth),sj,sjb,4,0)
355 endif
356 gb0=gb0*dble(sj)/xmax**delh *1.5d0 !rejection function norm.
357 if(ish.ge.6)write (ifch,*)'gb0,sj,z',gb0,sj,z
358
359 if(gb0.le.0.)then
360 write (ifmt,*)'gb0<0, smax,pth',smax,pth
361 iret=1
362 goto 16
363 endif
364
365 ! sharing of light cone momenta between soft preevolution and
366 ! hard interaction: ( first energy-momentum is shared according to
367 ! f_hard(yj)~zpm**(delh-dels-1) and then rejected as
368 ! w_rej ~sigma_hard_tot(yj) / exp(delh*yj)
369
370 4 continue
371 zpm=(xmin1+dble(rangen())*(xmax1-xmin1))**dble(1./(delh-dels)) !zpm
372 if(iclpro.ne.4.or.iqq.ne.1)then
373 call psjti0(sngl(zpm*s-pth),sgg,sggb,0,0)!inclusive (sgg) and born
374 call psjti0(sngl(zpm*s-pth),sgq,sgqb,0,1)
375 call psjti0(sngl(zpm*s-pth),sqq,sqqb,1,1)
376 call psjti0(sngl(zpm*s-pth),sqaq,sqaqb,-1,1)
377 call psjti0(sngl(zpm*s-pth),sqqp,sqqpb,1,2)
378 sqq=(sqq+sqaq+2.*(naflav-1)*sqqp)/naflav/2.
379 else
380 call psjti0(sngl(zpm*s-pth),sgq,sgqb,4,0)
381 call psjti0(sngl(zpm*s-pth),sqq,sqqb,4,1)
382 endif
383 xpmax=psutz(s,zpm*s,dble(amqpt**2)) !max x+ for sh pomeron
384 xpmin=zpm/xpmax
385 if(ish.ge.8)write (ifch,*)'zpm,xpmax,xpmin',zpm,xpmax,xpmin
386
387 if(iqq.eq.0)then
388 xp=xpmin*(xpmax/xpmin)**rangen() !lc+ share for the hard interaction
389 xm=zpm/xp !lc- share for the hard interaction
390 else
391 tmax=acos(dsqrt(xpmin)) !kinematical limits for t=cos(x+-)**2
392 tmin=acos(dsqrt(xpmax))
393 t=tmin+dble(rangen())*(tmax-tmin)
394 xp=cos(t)**2
395 xm=zpm/xp
396 endif
397
398 if(ish.ge.8)write(ifch,*)'zpm,xp,xm,xpmax,xpmin:',
399 * zpm,xp,xm,xpmax,xpmin
400
401 if(iqq.eq.0)then ! --------- sea-sea -----------
402
403 dwwgg1=0.
404 dwwgq1=0.
405 dwwqg1=0.
406 dwwqq1=0.
407
408 dwwgg2=0.
409 dwwgq2=0.
410 dwwqg2=0.
411 dwwqq2=0.
412
413 wwgg=sngl((1.d0-xp)*(1.d0-xm))**betpom
414 wwgq=sngl(1.d0-xp)**betpom*EsoftQZero(sngl(xm))
415 wwqg=sngl(1.d0-xm)**betpom*EsoftQZero(sngl(xp))
416 wwqq=EsoftQZero(sngl(xp))*EsoftQZero(sngl(xm))
417
418 if(idfpomr.eq.1)then
419 rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(zpm))
420 wwgg=(wwgg*z**(rp/rh)/rh+dwwgg1+dwwgg2)
421 * *(r2had(iclpro)+r2had(icltar))/z
422 wwgq=(wwgq*z**(rp/rh)/rh+dwwgq1+dwwgq2)
423 * *(r2had(iclpro)+r2had(icltar))/z
424 wwqg=(wwqg*z**(rp/rh)/rh+dwwqg1+dwwqg2)
425 * *(r2had(iclpro)+r2had(icltar))/z
426 wwqq=(wwqq*z**(rp/rh)/rh+dwwqq1+dwwqq2)
427 * *(r2had(iclpro)+r2had(icltar))/z
428 elseif(idfpomr.eq.2)then
429 rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(zpm))
430 wwgg=wwgg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
431 wwqg=wwqg/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
432 wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
433 wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
434 elseif(idfpomr.eq.3)then
435 rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(zpm))
436 wwgg=wwgg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
437 wwqg=wwqg/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
438 wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
439 wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
440 else
441 stop'psahot-idfpomr????'
442 endif
443
444 wwgg=wwgg*sgg*(1.-glusea)**2
445 wwgq=wwgq*sgq*(1.-glusea)*glusea
446 wwqg=wwqg*sgq*(1.-glusea)*glusea
447 wwqq=wwqq*sqq*glusea**2
448 gbyj=dlog(xpmax/xpmin)*dble(wwgg+wwgq+wwqg+wwqq)
449 wpi=wp0*xp !lc+ for the hard interaction
450 wmi=wm0*xm !lc+-for the hard interaction
451 gbyj=gbyj/zpm**dble(delh)/gb0 !rejection fu
452 if(gbyj.gt.1.d0.and.ish.ge.2)write(ifmt,*)'gbyj',gbyj
453 if(dble(rangen()).gt.gbyj)goto 4 !rejection
454 wp1=wp0-wpi
455 wm1=wm0-wmi
456 call pslcsh(wp1,wm1,wp2,wm2,amqt,dble(amqpt))
457
458 else ! --------- val-sea sea-val -----------
459
460 dwwg=0.
461 dwwq=0.
462
463 wwgq=sngl(1.d0-xm)**betpom
464 wwqq=EsoftQZero(sngl(xm))
465 if(idfpomr.eq.1)then
466 rh=r2had(iclpro)+r2had(icltar)-slopom*sngl(dlog(xm))
467 wwgq=(wwgq*z**(rp/rh)/rh+dwwg)
468 * *(r2had(iclpro)+r2had(icltar))/z
469 wwqq=(wwqq*z**(rp/rh)/rh+dwwq)
470 * *(r2had(iclpro)+r2had(icltar))/z
471 else
472 if(iqq.eq.1)then
473 rh=r2had(iclpro)+alam3p/2.-slopom*sngl(dlog(xm))
474 wwgq=wwgq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
475 wwqq=wwqq/rh*(r2had(iclpro)+alam3p/2.)*z**(rp/rh-1.)
476 else
477 rh=r2had(icltar)+alam3p/2.-slopom*sngl(dlog(xm))
478 wwgq=wwgq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
479 wwqq=wwqq/rh*(r2had(icltar)+alam3p/2.)*z**(rp/rh-1.)
480 endif
481 endif
482
483 wwgq=wwgq*sgq*(1.-glusea)*sngl(xp)**(-0.5+dels)
484 * *sngl(1.d0-xp)**(0.5+alpq)
485 wwqq=wwqq*sqq*glusea*sngl(xp)**(-0.5+dels)
486 * *sngl(1.d0-xp)**(0.5+alpq)
487
488 if(iqq.eq.1)then !valence quark-gluon hard interaction
489 uv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,1)
490 dv1=psdfh4(sngl(xp*xp1),q2min,0.,iclpro,2)
491 wpi=wp0*xp !lc+ for the hard interaction
492 wmi=wm0*xm !lc+-for the hard interaction
493 aks=rangen()
494 if(aks.lt.uv1/(uv1+dv1))then
495 if(iclpro.ne.4)then
496 iq1=1
497 else
498 iq1=4*(2.*int(.5+rangen())-1.)
499 endif
500 else
501 iq1=2
502 endif
503 gbyj=dble((wwgq+wwqq)*(uv1+dv1))*(1.d0-xp*xp1)**
504 * (-1.-alpq-alplea(iclpro))
505 else !gluon-valence quark hard interaction
506 xm=xp
507 xp=zpm/xm
508 uv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,1)
509 dv1=psdfh4(sngl(xm*xp2),q2min,0.,icltar,2)
510 wpi=wp0*xp !lc+ for the hard interaction
511 wmi=wm0*xm !lc+-for the hard interaction
512 aks=rangen()
513 if(aks.lt.uv1/(uv1+dv1))then
514 iq2=1
515 else
516 iq2=2
517 endif
518 gbyj=dble(wwgq+wwqq)*dble(uv1+dv1)*
519 * (1.d0-xm*xp2)**(-1.-alpq-alplea(icltar))
520 endif
521
522 gbyj=gbyj*(tmax-tmin)/zpm**delh /gb0 /2.1d0 !rejection
523 if(ish.ge.6)write (ifch,*)
524 * 'gbyj,zpm,xp,tmax,tmin,xpmax,xpmin',
525 * gbyj,zpm,xp,tmax,tmin,xpmax,xpmin
526 if(dble(rangen()).gt.gbyj)goto 4
527
528 wp1=wp0-wpi
529 wm1=wm0-wmi
530 if(ish.ge.8)write (ifch,*)'q_sea mass check',wp1*wm1,amqpt
531 if(iqq.eq.1)then !valence quark-gluon hard interaction
532 amq1=amqt(3)**2
533 s24=(amqt(2)+amqt(4))**2
534 else
535 amq1=amqt(1)**2
536 s24=(amqt(2)+amqt(4))**2
537 xp=xm
538 xm=zpm/xp
539 endif
540 x1max=psutz(wp1*wm1,dble(amq1),dble(s24))
541 x1min=dble(amq1)/x1max/wp1/wm1
542 t1min=(1.d0/x1max-1.d0)
543 t1max=(1.d0/x1min-1.d0)
5445 t1=t1min*(t1max/t1min)**dble(rangen())
545 if(ish.ge.8)write (ifch,*)'t1,t1min,t1max',t1,t1min,t1max
546 xq1=1.d0/(1.d0+t1)
547 if(dble(rangen()).gt.(xq1*(1.d0-xq1))**(1.+(-alpqua)))goto 5
548 if(iqq.eq.1)then !valence quark-gluon hard interacti
549 wm2=wm1*(1.d0-xq1)
550 wm1=wm1*xq1
551 wp1=wp1-dble(amq1)/wm1
552 if(ish.ge.8)write (ifch,*)'q_sea+ mass check',
553 * wp1*wm2,s24
554 wp1=wp1*psutz(wp1*wm2,dble(amqt(2)**2),dble(amqt(4)**2))
555 wm2=wm2-dble(amqt(2))**2/wp1
556 else
557 wp2=wp1*(1.d0-xq1)
558 wp1=wp1*xq1
559 wm1=wm1-dble(amq1)/wp1
560 if(ish.ge.8)write (ifch,*)'q_sea- mass check',
561 * wp2*wm1,s24
562 wm1=wm1*psutz(wm1*wp2,dble(amqt(4)**2),dble(amqt(2)**2))
563 wp2=wp2-amqt(4)**2/wm1
564 endif
565
566 endif ! -------------------
567
568 !-------------------------------
569 endif
570 !-------------------------------
571
572c-------------------------------------------------------------------------
573c flavor and momenta of end partons of the hard Pomeron
574c-------------------------------------------------------------------------
575
576 if(ish.ge.4)write(ifch,*)
577 & 'flavor and momenta of end partons of the hard Pomeron'
578
579 6 continue
580 wpi1=wpi
581 wmi1=wmi
582 nj=nj0 !initialization for the number of jets
583 if(ish.ge.8)write (ifch,*)'5-ww,smin',wpi*wmi,smin
584
585 rrr=rangen()
586 jqq=0
587
588 call iddeco(icp,jcp)
589 call iddeco(ict,jct)
590
591 iret1=0
592 iret2=0
593
594 if(iqq.eq.1.or.iqq.eq.2)then
595 if(rrr.lt.wwqq/(wwgq+wwqq))jqq=1
596 elseif(iqq.eq.0)then
597 if(rrr.lt.wwqg/(wwgg+wwqg+wwgq+wwqq))then
598 jqq=1
599 elseif(rrr.lt.(wwqg+wwgq)/(wwgg+wwqg+wwgq+wwqq))then
600 jqq=2
601 elseif(rrr.lt.(wwqg+wwgq+wwqq)/(wwgg+wwqg+wwgq+wwqq))then
602 jqq=3
603 endif
604 endif
605
606 if((iqq-1)*(iqq-3).eq.0)then
607
608 if(iabs(idptl(ip)).eq.1220)iq1=3-iq1 !proj=neutron
609 iqc(1)=iq1 !proj=particle
610 if(idptl(ip).lt.-1000)iqc(1)=-iq1 !proj=antiparticle
611 nj=nj+1
612 ifl1=iabs(iq1)
613
614 if(iqc(1).gt.0)then
615 call idsufl(ifl1,1,jcp,iret1) !remove valence quark from remnant
616 if(iret1.ne.1)ifl=idrafl(iclpro,jcp,2,'s',iret2) !take sea antiquark
617 elseif(iqc(1).lt.0)then
618 call idsufl(ifl1,2,jcp,iret1)
619 if(iret1.ne.1)ifl=idrafl(iclpro,jcp,1,'s',iret2)
620 else
621 call utstop('No quark for hard Pomeron+ in psahot!&')
622 endif
623
624 if(iret1.eq.1.or.iret2.eq.1)then
625 ifl=ifl1
626 if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 1)'
627 if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 1)'
628 call iddeco(icp,jcp)
629 iret1=0
630 iret2=0
631 endif
632
633 if(ish.ge.5)write(ifch,*)'flavor vq+,sqb+',iq1,-ifl
634
635 if(ifl.eq.3)then
636 iqj(nj)=-isign(ifl,iqc(1))*4/3
637 elseif(ifl.eq.4)then
638 iqj(nj)=-isign(ifl,iqc(1))/4
639 else
640 iqj(nj)=-isign(ifl,iqc(1))
641 endif
642 eqj(1,nj)=.5*sngl(wp1+dble(amqt(2))**2/wp1)
643 eqj(2,nj)=wp1-eqj(1,nj)
644 eqj(3,nj)=xxp2pr(ncolp,kcol)
645 eqj(4,nj)=xyp2pr(ncolp,kcol)
646 if(ish.ge.8)write (ifch,*)'q_v+ mass check',(eqj(1,nj)-
647 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
648 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
649 ncc(1,1)=nj
650 ncc(2,1)=0
651
652 else
653
654 nj=nj+1
655 if(idfpomr.ne.3)then
656 iflq=idrafl(iclpro,jcp,1,'s',iret1) !Choose flavor of sea quark.
657 if(iret1.ne.1)iflqb=idrafl(iclpro,jcp,2,'s',iret2) !Choose antiquark.
658 else
659 iflq=jdrafl(iclpro,jcp,1,iret1) !Choose flavor of sea quark.
660 iflqb=iflq !antiquark=quark (vertex end)
661 endif
662
663 if(iret1.eq.1.or.iret2.eq.1)then
664 iflqb=iflq
665 if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 2)'
666 if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 2)'
667 call iddeco(icp,jcp)
668 iret1=0
669 iret2=0
670 endif
671
672 if(ish.ge.5)write(ifch,*)'flavor sq+,sqb+',iflq,-iflqb
673
674 if(iflqb.eq.3)then
675 iqj(nj)=-iflqb*4/3
676 elseif(iflqb.eq.4)then
677 iqj(nj)=-iflqb/4
678 else
679 iqj(nj)=-iflqb
680 endif
681 if(iflq.eq.3)then
682 iqj(nj+1)=iflq*4/3
683 elseif(iflq.eq.4)then
684 iqj(nj+1)=iflq/4
685 else
686 iqj(nj+1)=iflq
687 endif
688 eqj(1,nj)=.5*sngl(wp1+dble(amqt(1))**2/wp1)
689 eqj(2,nj)=wp1-eqj(1,nj)
690 eqj(3,nj)=xxp1pr(ncolp,kcol)
691 eqj(4,nj)=xyp1pr(ncolp,kcol)
692 if(ish.ge.8)write (ifch,*)'q_s1+ mass check',(eqj(1,nj)-
693 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
694 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
695 eqj(1,nj+1)=.5*sngl(wp2+dble(amqt(2))**2/wp2)
696 eqj(2,nj+1)=wp2-eqj(1,nj+1)
697 eqj(3,nj+1)=xxp2pr(ncolp,kcol)
698 eqj(4,nj+1)=xyp2pr(ncolp,kcol)
699 nj=nj+1
700 if(ish.ge.8)write (ifch,*)'q_s2+ mass check',(eqj(1,nj)-
701 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
702 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
703
704 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.2)then
705 iqc(1)=0
706 ncc(1,1)=nj-1
707 ncc(2,1)=nj
708 else
709 iqc(1)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
7107 zg=1.d0-dble(rangen())*(1.d0-xp)
711 if(ish.ge.8)write (ifch,*)'6-zg,xp',zg,xp
712 if(dble(rangen()).gt.zg**dels*((1.d0-xp/zg)/(1.d0-xp))
713 * **betpom)goto 7
714 xg=xp/zg
715 wpq=wp0*(xg-xp)
716 wmq=1.d0/wpq
717 wmi1=wmi1-wmq
718 if(wmi1*wpi1.le.smin)goto 6
719
720 nj=nj+1
721 if(iabs(iqc(1)).eq.3)then
722 iqj(nj)=-iqc(1)*4/3
723 else
724 iqj(nj)=-iqc(1)
725 endif
726 eqj(1,nj)=.5*wmq
727 eqj(2,nj)=-eqj(1,nj)
728 eqj(3,nj)=0.
729 eqj(4,nj)=0.
730 if(ish.ge.8)write (ifch,*)'q_s3+ mass check',eqj(1,nj)**2-
731 * eqj(2,nj)**2-eqj(3,nj)**2-eqj(4,nj)**2
732 if(iqc(1).gt.0)then
733 ncj(1,nj)=nj-1
734 ncj(1,nj-1)=nj
735 ncj(2,nj)=0
736 ncj(2,nj-1)=0
737 ncc(1,1)=nj-2
738 ncc(2,1)=0
739 else
740 ncj(1,nj)=nj-2
741 ncj(1,nj-2)=nj
742 ncj(2,nj)=0
743 ncj(2,nj-2)=0
744 ncc(1,1)=nj-1
745 ncc(2,1)=0
746 endif
747 endif
748 endif
749
750 if((iqq-2)*(iqq-3).eq.0)then
751 if(iabs(idptl(maproj+it)).eq.1220)iq2=3-iq2 !targ=neutron
752 iqc(2)=iq2 !tar=particle (can not be pion or kaon !)
753 if(idptl(maproj+it).lt.-1000)iqc(2)=-iq2 !targ=antiparticle
754 ifl2=iabs(iq2)
755
756 nj=nj+1
757 if(iqc(2).gt.0)then
758 call idsufl(ifl2,1,jct,iret1) !remove valence quark from remnant
759 if(iret1.ne.1)ifl=idrafl(icltar,jct,2,'s',iret2) !take sea antiquark
760 elseif(iqc(2).lt.0)then
761 call idsufl(ifl2,2,jct,iret1)
762 if(iret1.ne.1)ifl=idrafl(icltar,jct,1,'s',iret2)
763 else
764 call utstop('No quark for hard Pomeron- in psahot!&')
765 endif
766
767 if(iret1.eq.1.or.iret2.eq.1)then
768 ifl=ifl2
769 if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 3)'
770 if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 3)'
771 call iddeco(ict,jct)
772 iret1=0
773 iret2=0
774 endif
775
776 if(ish.ge.5)write(ifch,*)'flavor vq-,sqb-',iq2,-ifl
777
778 if(ifl.eq.3)then
779 iqj(nj)=-isign(ifl,iqc(2))*4/3
780 elseif(ifl.eq.4)then
781 iqj(nj)=-isign(ifl,iqc(2))/4
782 else
783 iqj(nj)=-isign(ifl,iqc(2))
784 endif
785
786 eqj(1,nj)=.5*sngl(wm1+dble(amqt(4))**2/wm1)
787 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
788 eqj(3,nj)=xxm1pr(ncolp,kcol)
789 eqj(4,nj)=xym1pr(ncolp,kcol)
790 if(ish.ge.8)write (ifch,*)'q_v- mass check',(eqj(1,nj)
791 * +eqj(2,nj))*(eqj(1,nj)-eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
792 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
793 ncc(1,2)=nj
794 ncc(2,2)=0
795 else
796 nj=nj+1
797 if(idfpomr.ne.2)then
798 iflq=idrafl(icltar,jct,1,'s',iret1) !Choose flavor of sea quark.
799 if(iret1.ne.1)iflqb=idrafl(icltar,jct,2,'s',iret2) !Choose antiquark.
800 else
801 iflq=jdrafl(icltar,jct,1,iret1) !Choose flavor of sea quark.
802 iflqb=iflq !antiquark=quark (vertex end)
803 endif
804
805 if(iret1.eq.1.or.iret2.eq.1)then
806 iflqb=iflq
807 if(ish.ge.3)write(ifmt,*)'Not enough place in rem (psahot 4)'
808 if(ish.ge.5)write(ifch,*)'Not enough place in rem (psahot 4)'
809 call iddeco(ict,jct)
810 iret1=0
811 iret2=0
812 endif
813
814 if(ish.ge.5)write(ifch,*)'flavor sq-,sqb-',iflq,-iflqb
815
816 if(iflqb.eq.3)then
817 iqj(nj)=-iflqb*4/3
818 elseif(iflqb.eq.4)then
819 iqj(nj)=-iflqb/4
820 else
821 iqj(nj)=-iflqb
822 endif
823 if(iflq.eq.3)then
824 iqj(nj+1)=iflq*4/3
825 elseif(iflq.eq.4)then
826 iqj(nj+1)=iflq/4
827 else
828 iqj(nj+1)=iflq
829 endif
830 eqj(1,nj)=.5*sngl(wm1+dble(amqt(3))**2/wm1)
831 eqj(2,nj)=eqj(1,nj)-sngl(wm1)
832 eqj(3,nj)=xxm2pr(ncolp,kcol)
833 eqj(4,nj)=xym2pr(ncolp,kcol)
834 if(ish.ge.8)write (ifch,*)'q_s1- mass check',(eqj(1,nj)-
835 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
836 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
837 eqj(1,nj+1)=.5*sngl(wm2+dble(amqt(4))**2/wm2)
838 eqj(2,nj+1)=eqj(1,nj+1)-wm2
839 eqj(3,nj+1)=xxm1pr(ncolp,kcol)
840 eqj(4,nj+1)=xym1pr(ncolp,kcol)
841 nj=nj+1
842 if(ish.ge.8)write (ifch,*)'q_s2- mass check',(eqj(1,nj)-
843 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
844 eqj(1,nj)=sqrt(eqj(2,nj)**2+eqj(3,nj)**2+eqj(4,nj)**2)
845
846 if(jqq.eq.0.or.iqq.eq.0.and.jqq.eq.1)then
847 iqc(2)=0
848 ncc(1,2)=nj-1
849 ncc(2,2)=nj
850 else
851 iqc(2)=int(3*rangen()+1.)*(2.*int(.5+rangen())-1.)
8528 zg=1.d0-dble(rangen())*(1.d0-xm)
853 if(ish.ge.8)write (ifch,*)'7-zg,xm',zg,xm
854 if(rangen().gt.zg**dels*((1.d0-xm/zg)/(1.d0-xm))**betpom)
855 * goto 8
856 xg=xm/zg
857 wmq=wm0*(xg-xm)
858 wpq=1.d0/wmq
859 wpi1=wpi1-wpq
860 if(wmi1*wpi1.le.smin)goto 6
861
862 nj=nj+1
863 if(iabs(iqc(2)).eq.3)then
864 iqj(nj)=-iqc(2)*4/3
865 else
866 iqj(nj)=-iqc(2)
867 endif
868 eqj(1,nj)=.5*sngl(wpq)
869 eqj(2,nj)=eqj(1,nj)
870 eqj(3,nj)=0.
871 eqj(4,nj)=0.
872 if(ish.ge.8)write (ifch,*)'q_s3- mass check',(eqj(1,nj)-
873 * eqj(2,nj))*(eqj(1,nj)+eqj(2,nj))-eqj(3,nj)**2-eqj(4,nj)**2
874 if(iqc(2).gt.0)then
875 ncj(1,nj)=nj-1
876 ncj(1,nj-1)=nj
877 ncj(2,nj)=0
878 ncj(2,nj-1)=0
879 ncc(1,2)=nj-2
880 ncc(2,2)=0
881 else
882 ncj(1,nj)=nj-2
883 ncj(1,nj-2)=nj
884 ncj(2,nj)=0
885 ncj(2,nj-2)=0
886 ncc(1,2)=nj-1
887 ncc(2,2)=0
888 endif
889 endif
890 endif
891 if(jqq.ne.0)then
892 if(iqq.ne.0.or.iqq.eq.0.and.jqq.eq.3)then
893 if(iclpro.ne.4.or.iqq.ne.1)then
894 call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,1,1)
895 call psjti0(sngl(wpi1*wmi1-pth),sqaq1,sqaqb1,-1,1)
896 call psjti0(sngl(wpi1*wmi1-pth),sqqp1,sqqpb1,1,2)
897 sqq1=(sqq1+sqaq1+2.*(naflav-1)*sqqp1)/naflav/2.
898 else
899 call psjti0(sngl(wpi1*wmi1-pth),sqq1,sqqb1,4,1)
900 endif
901 gbs=sqq1/sqq
902 else
903 call psjti0(sngl(wpi1*wmi1-pth),sgq1,sgqb1,0,1)
904 gbs=sgq1/sgq
905 endif
906 if(ish.ge.8)write (ifch,*)'gbs',gbs
907 if(rangen().gt.gbs)goto 6
908 endif
909
910
911c---------------------------------------------------------------
912c inner partons of the hard Pomeron
913c---------------------------------------------------------------
914
915 if(ish.ge.4)write(ifch,*)
916 & 'inner partons of the hard Pomeron'
917 nj00=nj
918 wpi=wpi1
919 wmi=wmi1
920 si=wpi*wmi-pxh**2-pyh**2 !total energy squared for the hard
921 if(ish.ge.7)write (ifch,*)'si,wpi,wmi',si,wpi,wmi
922
923 ept(1)=.5d0*(wpi+wmi)
924 ept(2)=.5d0*(wpi-wmi)
925 ept(3)=pxh
926 ept(4)=pyh
927 qmin(1)=q2min !effective momentum cutoff above current la
928 qmin(2)=q2min !effective momentum cutoff below current la
929 qminn=max(qmin(1),qmin(2)) !effective momentum cutoff for the bo
930c si=psnorm(ept) !4-momentum squared for the hard pomeron
931 jfirst=1
932 jj=int(1.5+rangen())
933 nemis(1)=0
934 nemis(2)=0
935
936 9 continue ! <<<<----------- ladder rung ---------------------------
937
938 if(ish.ge.4)write(ifch,*)'ladder rung'
939 pt2=ept(3)**2+ept(4)**2
940 pt=sqrt(sngl(pt2))
941 if(iabs(iqc(1)).ne.4)then
942 q2mass=0.
943 else
944 q2mass=qcmass**2
945 si=si-dble(q2mass)
946 endif
947 s2min=4.*qminn+q2mass !mass cutoff for born scattering
948 wwmin=5.*qminn+q2mass-2.*pt*sqrt(q2ini)
949 wwmin=(wwmin+sqrt(wwmin**2+4.*(q2mass+pt2)*(qminn-q2ini)))
950 */(1.-q2ini/qminn)/2.
951 if(ish.ge.5)write(ifch,*)'qminn,q2mass,pt,wwmin:'
952 if(ish.ge.5)write(ifch,*)qminn,q2mass,pt,wwmin
953
954 wpt(1)=ept(1)+ept(2) !lc+ for the current jet emissi
955 wpt(2)=ept(1)-ept(2) !lc- for the current jet emissi
956
957 if(iabs(iqc(1)).ne.4)then
958 if(jfirst.eq.1)then
959 sj=psjti(qmin(jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
960 sj2=psjti1(qmin(3-jj),qmin(jj),qqcut
961 * ,sngl(si),iqc(3-jj),iqc(jj),0)
962 if(ish.ge.5)write(ifch,*)'si,sj,sj2:',si,sj,sj2
963 if(rangen().gt.sj2/sj.and.si.gt.1.1d0*dble(wwmin))goto 112
964 jfirst=0
965 jj=3-jj
966 sj=sj2
967 goto 111
968 endif
969 sj=psjti1(qmin(jj),qmin(3-jj),qqcut,sngl(si),iqc(jj),iqc(3-jj),0)
970111 sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si),iqc(1),iqc(2),0) !born parton-parton
971 else
972 sjord=psjci(qmin(2),sngl(si),iqc(2))
973 sj=sjord
974 sjb=psbint(qmin(1),qmin(2),qqcut,sngl(si)+q2mass,iqc(1),iqc(2),0)
975 if(qmin(2).eq.q2min)then
976 wwmins=2.5*q2min*(1.+sqrt(1.+
977 * 4.*(pt2+q2mass)/q2min))
978 if(si.gt.dble(wwmins))call psjti0(sngl(si)+q2mass,sj,sjb
979 * ,iqc(1),iqc(2))
980 endif
981 endif
982 if(ish.ge.5)write(ifch,*)'si,pt2,wwmin,s2min,sj,sjb,iqc:'
983 if(ish.ge.5)write(ifch,*)si,pt2,wwmin,s2min,sj,sjb,iqc
984
985 if(si.lt.1.1d0*dble(wwmin))goto 12 !------>>>>>>>
986 if(rangen().lt.sjb/sj)goto 12 !------>>>>>>>
987
988 if(iabs(iqc(1)).eq.4)jj=min(2,int(sjord/sj+rangen())+1) !?????????
989
990112 continue
991
992 if(iabs(iqc(jj)).ne.4)then
993
994 discr=dble((sngl(si)+2.*pt*sqrt(q2ini)-q2mass)**2
995 * -4.*q2ini*(5.*sngl(si)+q2mass+sngl(pt2)))
996 if(discr.lt.0.d0.and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
997 * discr,si,pt,wwmin
998 discr=dsqrt(discr)
999 qqmax=(si+2.d0*dble(pt*sqrt(q2ini))-dble(q2mass)+discr)/2.d0
1000 * /(5.d0+(dble(q2mass)+pt2)/si)
1001 qqmin=qqmax-discr/(5.d0+(dble(q2mass)+pt2)/si)
1002 if(s2min.gt.4.d0*qqmin+dble(q2mass))then
1003 xmm=.5d0*(si-s2min+2.d0*dble(pt*sqrt(q2ini)))
1004 discr=xmm**2-si*dble(q2ini)*(1.d0+(dble(q2mass)+pt2)/si)
1005 if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
1006 * discr,si,pt,wwmin
1007 qqmin=(xmm-dsqrt(discr))/(1.d0+(dble(q2mass)+pt2)/si)
1008 endif
1009
1010 xmax=min(1.d0-dble(q2ini)/qqmax,.9999d0)
1011 if(qqmin.lt.dble(qmin(jj)))then
1012 qqmin=dble(qmin(jj))
1013 xmin=max(1.d0-((dble(pt)*dsqrt(qqmin)+dsqrt(pt2*qqmin+
1014 * si*(si-s2min-qqmin*(1.d0+(dble(q2mass)+pt2)/si))))/si)**2,
1015 * (s2min+qqmin*(1.d0+(dble(q2mass)+pt2)/si)-2.d0*dble(pt)
1016 * *dsqrt(qqmin))/si)
1017 if(xmin.le.0.d0)xmin=s2min/si
1018 else
1019 xmin=1.d0-dble(q2ini)/qqmin
1020 endif
1021
1022 qm0=qmin(jj)
1023 xm0=1.0-dble(q2ini/qm0)
1024 if(xm0.gt..999*xmax.or.xm0.lt.1.001*xmin)then
1025 xm0=.5d0*(xmax+xmin)
1026 endif
1027 s2max=sngl(xm0*si-qm0*(dble(q2mass)+pt2)/si)
1028 * +2.*pt*sqrt(q2ini)
1029 xx=xm0
1030
1031 if(iabs(iqc(1)).ne.4)then
1032 if(jfirst.eq.1)then
1033 sj0=psjti(qm0,qqcut,s2max,0,iqc(3-jj),0)*
1034 * psfap(xx,iqc(jj),0)+
1035 * psjti(qm0,qqcut,s2max,7,iqc(3-jj),0)
1036 * *psfap(xx,iqc(jj),1)
1037 else
1038 sj0=psjti1(qm0,qmin(3-jj),qqcut,s2max,0,iqc(3-jj),0)*
1039 * psfap(xx,iqc(jj),0)+
1040 * psjti1(qm0,qmin(3-jj),qqcut,s2max,7,iqc(3-jj),0)
1041 * *psfap(xx,iqc(jj),1)
1042 endif
1043 else
1044 sj0=psjci(qm0,s2max,0)*psfap(xx,iqc(jj),0)+
1045 * psjci(qm0,s2max,1)*psfap(xx,iqc(jj),1)
1046 endif
1047
1048 gb0=dble(sj0/log(q2ini/qcdlam)*qm0*5.)*psuds(qm0,iqc(jj))
1049 if(ish.ge.5)write(ifch,*)'gb0,qm0,xm0,s2max:',
1050 * gb0,qm0,xm0,s2max
1051c if(gb0.le.0.)stop'gb0<=0'
1052 if(gb0.le.0.d0)then
1053 write(ifmt,*)'gb0.le.0. si,pt2:',si,pt2
1054 iret=1
1055 goto 16
1056 endif
1057
1058 if(xm0.le..5d0)then
1059 gb0=gb0*xm0**(1.-delh)
1060 else
1061 gb0=gb0*(1.d0-xm0)*2.d0**delh
1062 endif
1063
1064 xmin2=max(.5d0,xmin)
1065 xmin1=xmin**delh
1066 xmax1=min(xmax,.5d0)**delh
1067 if(xmin.ge..5d0)then
1068 djl=1.
1069 elseif(xmax.le..5d0)then
1070 djl=0.
1071 else
1072 djl=1./(1.+((2.*sngl(xmin))**delh-1.)/delh/
1073 * log(2.*(1.-sngl(xmax))))
1074 endif
1075
1076 else
1077
1078 xmin=5.d0*dble(q2min)/si
1079 xmax=min(si/(si+5.0*(pt2+dble(q2mass))),.9999d0)
1080 qqmax=xmax*si/5.d0
1081 qqmin=dble(q2min)
1082
1083 qm0=sngl(qqmin)
1084 xm0=2.d0/(1.d0+sqrt(1.d0+4.d0*(pt2+dble(q2mass))/qm0))
1085 if(xm0.gt..999d0*xmax.or.xm0.lt.1.001d0*xmin)then
1086 xm0=.5d0*(xmax+xmin)
1087 endif
1088 s2max=sngl(xm0*si)
1089 xx=xm0
1090
1091 sj0=psjti(qm0,qmin(3-jj),s2max,0,iqc(3-jj),0)*
1092 * psfap(xx,iqc(jj),0)
1093 gb0=dble(sj0/log(qm0/qcdlam)*qm0 *5.)
1094 gb0=gb0*xm0**(1.-delh)
1095 if(gb0.le.0.d0)then
1096 if(ish.ge.2)write(ifch,*)'gb0.le.0. (charm) si,pt2:',si,pt2
1097 iret=1
1098 goto 16
1099 endif
1100 djl=0.
1101 xmin1=xmin**delh
1102 xmax1=xmax**delh
1103
1104 endif
1105
1106 if(ish.ge.5)write(ifch,*)'xmin,xmax,qqmin,qqmax:',
1107 *xmin,xmax,qqmin,qqmax
1108
1109 ntry=0
111010 continue
1111 ntry=ntry+1
1112 if(ntry.gt.5000000)then
1113 if(ish.ge.1)write(*,*)'Reject hard string (too many rejection)'
1114 &,kcol,ncolp,nptl,gb7
1115 iret=1
1116 goto 16
1117 endif
1118 if(rangen().gt.djl)then !lc momentum share in the cur
1119 x=(xmin1+dble(rangen())*(xmax1-xmin1))**(1./delh)
1120 else
1121 x=1.d0-(1.d0-xmin2)*((1.d0-xmax)/(1.d0-xmin2))**rangen()
1122 endif
1123 qq=qqmin/(1.d0+dble(rangen())*(qqmin/qqmax-1.d0))
1124 if(ish.ge.7)write(ifch,*)'x,qq:',x,qq,ntry
1125
1126 if(iabs(iqc(jj)).ne.4)then
1127
1128 qt2=qq*(1.d0-x)
1129 if(qt2.lt.dble(q2ini))then
1130 if(ish.ge.7)write(ifch,*)'qt2:',qt2
1131 goto 10
1132 endif
1133
1134 qmin2=max(qminn,sngl(qq))
1135 qt=dsqrt(qt2)
1136 call pscs(bcos,bsin)
1137 ep3(3)=sngl(qt)*bcos !new parton pt-s
1138 ep3(4)=sngl(qt)*bsin
1139 ptnew=(ept(3)-dble(ep3(3)))**2+(ept(4)-dble(ep3(4)))**2
1140 s2min2=4.*qmin2+q2mass
1141
1142 s2=x*(si-qq)-ptnew-qq*(dble(q2mass)+pt2)/si+pt2 !new ladder mass
1143 if(s2.lt.dble(s2min2))then
1144 if(ish.ge.7)write(ifch,*)'s2,s2min2:',s2,s2min2
1145 goto 10 !rejection in case of too low mass
1146 endif
1147
1148 xx=x
1149 if(iabs(iqc(1)).ne.4)then
1150 if(jfirst.eq.1)then
1151 sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1152 if(iqc(jj).ne.0)then !f1 - f2
1153 sj2=psjti(sngl(qq),qqcut,sngl(s2),iqc(jj),iqc(3-jj),0)
1154 elseif(iqc(3-jj).eq.0)then !f1 - g
1155 sj2=psjti(sngl(qq),qqcut,sngl(s2),1,0,0)
1156 else !g - f2
1157 sj2=psjti(sngl(qq),qqcut,sngl(s2),1,1,0)/naflav/2. !q - q
1158 * +psjti(sngl(qq),qqcut,sngl(s2),-1,1,0)/naflav/2. !q~ - q
1159 * +psjti(sngl(qq),qqcut,sngl(s2),1,2,0)*(1.-1./naflav) !q' - q
1160 endif
1161 else
1162 sj1=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),0,iqc(3-jj),0)
1163 if(iqc(jj).ne.0)then !f1 - f2
1164 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1165 * ,sngl(s2),iqc(jj),iqc(3-jj),0)
1166 elseif(iqc(3-jj).eq.0)then !f1 - g
1167 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut,sngl(s2),1,0,0)
1168 else !g - f2
1169 sj2=psjti1(sngl(qq),qmin(3-jj),qqcut
1170 * ,sngl(s2),1,1,0)/naflav/2. !q - q
1171 * +psjti1(sngl(qq),qmin(3-jj),qqcut
1172 * ,sngl(s2),-1,1,0)/naflav/2. !q~ - q
1173 * +psjti1(sngl(qq),qmin(3-jj),qqcut
1174 * ,sngl(s2),1,2,0)*(1.-1./naflav)!q' - q
1175 endif
1176 endif
1177 else
1178 sj1=psjci(sngl(qq),sngl(s2),0)
1179 sj2=psjci(sngl(qq),sngl(s2),1)
1180 endif
1181
1182 !...gb7 - rejection function for x and q**2 simulation
1183 gb7=dble((sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))/
1184 * log(sngl(qt2)/qcdlam))*psuds(sngl(qq),iqc(jj))*qq/gb0
1185
1186 if(x.le..5d0)then
1187 gb7=gb7*x**(1.-delh)
1188 else
1189 gb7=gb7*(1.d0-x)*2.d0**delh
1190 endif
1191 if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1192 * gb7,qq,x,qm0,xm0
1193 if(ish.ge.7)write(ifch,*)'gb7:',gb7
1194 if(dble(rangen()).gt.gb7)goto 10
1195 else
1196 qmin2=max(qminn,sngl(qq))
1197 s2min2=4.*qmin2
1198 s2=x*si-qq !new ladder mass
1199 if(s2.lt.dble(s2min2))goto 10 !rejection in case of too low mass
1200
1201 call pscs(bcos,bsin)
1202 xmm=x*(ept(3)*dble(bcos)+ept(4)*dble(bsin))
1203 discr=xmm**2+qq*(1.d0-x)-x**2*(pt2+dble(q2mass))
1204 if(discr.lt.0.d0)goto 10
1205 qt=xmm+dsqrt(discr)
1206 ep3(3)=sngl(ept(3)-qt*dble(bcos)) !new parton pt-s
1207 ep3(4)=sngl(ept(4)-qt*dble(bsin))
1208 qt2=dble(ep3(3)**2+ep3(4)**2)
1209
1210 xx=x
1211 sj1=psjti(sngl(qq),qqcut,sngl(s2),0,iqc(3-jj),0)
1212 sj2=0.
1213 !.... gb7 - rejection function for x and q**2 simulation
1214 gb7=dble(sj1*psfap(xx,iqc(jj),0)/
1215 * log(sngl(qq)/qcdlam))*qq/gb0
1216 gb7=gb7*x**(1.-delh)
1217 if(gb7.gt.1..and.ish.ge.2)write(ifmt,*)'gb7,qq,x,qm0,xm0',
1218 * gb7,qq,x,qm0,xm0
1219 if(ish.ge.7)write(ifch,*)'gb7:',gb7
1220 if(dble(rangen()).gt.gb7)goto 10
1221 endif
1222
1223 nqc(2)=0
1224 iqnew=iqc(jj)
1225 if(rangen().lt.sj1/(sj1+sj2))then
1226 if(iqc(jj).eq.0)then
1227 jt=1
1228 jq=int(1.5+rangen())
1229 nqc(1)=ncc(jq,jj)
1230 else
1231 jt=2
1232 if(iqc(jj).gt.0)then
1233 jq=1
1234 else
1235 jq=2
1236 endif
1237 nqc(1)=0
1238 iqnew=0
1239 endif
1240 iq1=iqc(jj)
1241
1242 else
1243 if(iqc(jj).ne.0)then
1244 iq1=0
1245 jt=3
1246 if(iqc(jj).gt.0)then
1247 jq=1
1248 else
1249 jq=2
1250 endif
1251 nqc(1)=ncc(1,jj)
1252 else
1253 jt=4
1254 jq=int(1.5+rangen())
1255 iq1=int(naflav*rangen()+1.)*(3-2*jq)
1256 iqnew=-iq1
1257 nqc(1)=ncc(jq,jj)
1258 endif
1259 endif
1260 if(iabs(iqc(jj)).ne.4)then
1261 q2part=0.
1262c write(*,*)'sem:',wpt(jj),pt2,q2mass,wpt(3-jj),jj
1263 if(dabs(wpt(3-jj)).gt.1.d-20)then
1264 wplc=wpt(jj)-(pt2+dble(q2mass))/wpt(3-jj)
1265 else
1266 if(ish.gt.1)write(ifmt,*)'Problem with wpt in sem',wpt(3-jj)
1267 wplc=wpt(jj)-(pt2+dble(q2mass))
1268 endif
1269 qp2max=max(0.,sngl(qt2))
1270 else
1271 q2part=q2mass
1272 wplc=wpt(jj)
1273 qp2max=max(0.,sngl(qq))
1274 endif
1275 eprt=max(dsqrt(qt2+dble(q2part)),.5d0*((1.d0-x)*wplc+
1276 *(qt2+dble(q2part))/(1.d0-x)/wplc))
1277 pl=((1.d0-x)*wplc-eprt)*dble(3-2*jj)
1278 if(iq1.eq.0)then
1279 iq2ini=9
1280 else
1281 iq2ini=iq1
1282 endif
1283 if(ish.ge.5)write(ifch,*)'qq,eprt,iq2ini',qq,eprt,iq2ini,
1284 *eprt**2-qt2
1285 ntest=0
128611 ntest=ntest+1
1287 call timsh1(qp2max,sngl(eprt),iq2ini)
1288 amprt=pprt(5,1)**2
1289 plprt=max(0.d0,eprt**2-qt2)-dble(amprt)
1290 if(plprt.lt.0.d0)then
1291 if(ntest.lt.10000)then
1292 goto 11
1293 else
1294 iret=1
1295 goto 16
1296 endif
1297 endif
1298 ep3(1)=sngl(eprt)
1299 ep3(2)=sngl(dsqrt(plprt))
1300 if(pl.lt.0.d0)ep3(2)=-ep3(2)
1301 ey(1)=1.
1302 ey(2)=1.
1303 ey(3)=1.
1304 do i=1,4
1305 ept1(i)=ept(i)-dble(ep3(i))
1306 enddo
1307 call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1308 s2new=psnorm(ept1)
1309
1310 if(iabs(iqc(jj)).ne.4.and.s2new-q2mass.gt.s2min2.or.
1311 *iabs(iqc(jj)).eq.4.and.s2new.gt.s2min2)then
1312 if(iabs(iqc(1)).ne.4.or.jj.eq.1)then
1313 if(jfirst.eq.1)then
1314 gb=dble(psjti(sngl(qq),qqcut,s2new,iqnew,iqc(3-jj),0))
1315 else
1316 gb=dble(psjti1(sngl(qq),qmin(3-jj),qqcut,s2new,iqnew,
1317 * iqc(3-jj),0))
1318 endif
1319 else
1320 gb=dble(psjci(sngl(qq),s2new-q2mass,iqnew))
1321 endif
1322 if(iqnew.eq.0)then
1323 gb=gb/dble(sj1)
1324 else
1325 gb=gb/dble(sj2)
1326 endif
1327 if(dble(rangen()).gt.gb)goto 10
1328 else
1329 goto 10
1330 endif
1331
1332 if(ish.ge.6)write(ifch,*)'jt,jj,jq,nqc:',jt,jj,jq,nqc
1333 nprtjx=0
1334 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1335
1336 if(jt.eq.1)then
1337 ncc(jq,jj)=nqc(2)
1338 elseif(jt.eq.2)then
1339 ncc(jq,jj)=ncc(1,jj)
1340 ncc(3-jq,jj)=nqc(1)
1341 elseif(jt.eq.3)then
1342 ncc(1,jj)=nqc(2)
1343 elseif(jt.eq.4)then
1344 ncc(1,jj)=ncc(3-jq,jj)
1345 endif
1346 iqc(jj)=iqnew
1347
1348 do i=1,4
1349 ept(i)=ept1(i)
1350 enddo
1351 ! c.m. energy squared, minimal 4-momentum transfer square and gluon 4-v
1352 ! for the next ladder run
1353 qmin(jj)=sngl(qq)
1354 qminn=qmin2
1355 si=dble(s2new)
1356 nemis(jj)=nemis(jj)+1
1357
1358 goto 9 ! ---------------next ladder rung ------------>>>>>>>>>>>
1359
1360
1361 12 continue !------------------- Born process------------------------
1362
1363 if(ish.ge.4)write(ifch,*)'Born process'
1364 if(ish.ge.6)write(ifch,*)'iqc,si:',iqc,si
1365 qq=dble(qminn)
1366
1367 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1368 xpprbor(ncolp,kcol)=(ept(1)+ept(2))/dble(engy)
1369 xmprbor(ncolp,kcol)=(ept(1)-ept(2))/dble(engy)
1370 nemispr(1,ncolp,kcol)=nemis(1)
1371 nemispr(2,ncolp,kcol)=nemis(2)
1372 ! write(*,'(a,2f8.3,i3,3x,2i3)')'------------'
1373 ! * ,xpprbor(ncolp,kcol),xmprbor(ncolp,kcol),nj-nj00,nemis
1374 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1375
1376 call psabor(sngl(si),sngl(qq),iqc,ncc,ept,0,iptl,bx)
1377
1378 !kkkkkkkkkkkkk out Born partons without timelike cascade
1379 if(nprtjx.eq.2)then
1380 do ii=1,2
1381 ptprboo(ii,ncolp,kcol)=sqrt(pprtx(1,ii)**2+pprtx(2,ii)**2)
1382 enddo
1383 else
1384 stop'psahot: should not happen!!!! '
1385 endif
1386 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1387 ! ptxxx=max( ptprboo(1,ncolp,kcol) , ptprboo(2,ncolp,kcol) )
1388 ! print*,sqrt(wm0*wp0),sqrt(wmi*wpi),ptxxx !++++++++++++++++++++++++++
1389 ! if(ptxxx.lt.10)goto1
1390 ! print*,' ++++++++++++++++++++++++++++++++++++ '
1391 !kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
1392 if(nj.ne.0.)then
1393 do i=1,nj
1394 do l=1,6
1395 bxj(l,i)=bx(l)
1396 enddo
1397 ityj(i)=ity
1398 iorj(i)=iptl
1399 enddo
1400 endif
1401
1402 call psjarr(jfl) !kinky strings formation
1403
1404 if(jfl.eq.0.and.ish.ge.4)write(ifch,*)
1405 *'jfl,nj,nptl',jfl,nj,nptl
1406 if(jfl.eq.0)goto 1
1407
1408c --- update remnant flavour ---
1409
1410 iret1=0
1411 call idenco(jcp,icp,iret1)
1412 if(iret1.eq.1)call utstop('Problem with proj rem in psahot!&')
1413 iret2=0
1414 call idenco(jct,ict,iret2)
1415 if(iret2.eq.1)call utstop('Problem with targ rem in psahot!&')
1416
1417 do i=1,2
1418 icproj(i,ip)=icp(i)
1419 ictarg(i,it)=ict(i)
1420 enddo
1421
1422 if(ish.ge.3)write(ifch,*)'End psahot (icp,ict):',icp,ict
1423
1424 if(iep(ip).ne.1)then
1425 call idtr4(idptl(ip),icini) !excited remnant ?
1426 if(ish.ge.5)write(ifch,*)'icini proj',icini
1427 & ,(icp(1)-icini(1)),(icp(2)-icini(2))
1428 if((icp(1)-icini(1))+(icp(2)-icini(2)).ne.0)iep(ip)=1
1429 endif
1430 if(iet(it).ne.1)then
1431 call idtr4(idptl(maproj+it),icini)
1432 if(ish.ge.5)write(ifch,*)'icini targ',icini
1433 & ,(ict(1)-icini(1)),(ict(2)-icini(2))
1434 if((ict(1)-icini(1))+(ict(2)-icini(2)).ne.0)iet(it)=1
1435 endif
1436 if(ish.ge.5)write(ifch,*)'iep,iet ',iep(ip),iet(it)
1437
1438c ------------------------------
1439
144016 continue
1441 call utprix('psahot',ish,ishini,3)
1442
1443 q2fin=q2finsave
1444 return
1445 end
1446
1447c------------------------------------------------------------------------
1448 subroutine psjarr(jfl)
1449c-----------------------------------------------------------------------
1450c
1451c final jets rearrangement according to their colour connection
1452c and write to /cptl/
1453c
1454c jfl - flag for the rejection (in case of jfl=0)
1455c-----------------------------------------------------------------------
1456c Input:
1457
1458 parameter (mjstr=20000)
1459 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1460
1461c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
1462c bxj(1:4,k) - coordinates for k-th parton formation point;
1463c iqj(k) - ID (qgs) for k-th parton;
1464c ncj(1:2,k) - colour connected partons indexes for k-th parton;
1465c nj - number of partons
1466c-----------------------------------------------------------------------
1467 dimension mark(mjstr)
1468 double precision ept(4)
1469 include 'epos.inc'
1470 include 'epos.incsem'
1471
1472 if(nj.eq.0)then
1473 jfl=1
1474 return
1475 endif
1476c do k=1,nj !???????????????
1477c eqj(1,k)=dsqrt(0d0+eqj(2,k)**2+eqj(3,k)**2+eqj(4,k)**2)
1478c enddo
1479 if(ish.ge.3)then
1480 write (ifch,*)'psjarr: nj',nj
1481 do k=1,nj
1482 if(iqj(k).ne.0)ncj(2,k)=0
1483 write(ifch,'(a,i4)')' parton',k
1484 write(ifch,'(i3,2x,4e10.3,2x,2i3)')iqj(k)
1485 * ,(eqj(j,k),j=1,4),(ncj(j,k),j=1,2)
1486 enddo
1487 endif
1488
1489 jfl=0
1490 do i=1,nj
1491 mark(i)=1
1492 enddo
1493 nptl0=nptl
1494 if(ish.ge.5)write (ifch,*)'nptl,nj',nptl,nj
1495
14961 continue
1497 do ij=1,nj
1498 if(mark(ij).ne.0.and.iqj(ij).ne.0)goto 2
1499 enddo
15002 continue
1501 jfirst=1
1502 if(iabs(iqj(ij)).le.2)then
1503 am1=amhadr(1)
1504 elseif(iabs(iqj(ij)).eq.4)then
1505 am1=amhadr(3)
1506 elseif(iabs(iqj(ij)).eq.40)then
1507 am1=qcmass
1508 else
1509 am1=amhadr(2)
1510 endif
1511
1512 do i=1,4
1513 ept(i)=0.
1514 enddo
1515
15163 continue
1517 call pspawr(ij)
1518 mark(ij)=0
1519
1520 do i=1,4
1521 ept(i)=ept(i)+eqj(i,ij)
1522 enddo
1523
1524 if(iqj(ij).ne.0)then
1525 if(jfirst.ne.1)then
1526 if(iabs(iqj(ij)).le.2)then
1527 am2=amhadr(1)
1528 elseif(iabs(iqj(ij)).eq.4)then
1529 am2=amhadr(3)
1530 elseif(iabs(iqj(ij)).eq.40)then
1531 am2=qcmass
1532 else
1533 am2=amhadr(2)
1534 endif
1535 amj=(am1+am2+stmass)**2
1536 sm=psnorm(ept)
1537 if(sm.lt.amj)then
1538 nptl=nptl0
1539 goto999
1540 endif
1541
1542 if(nptl-nptl0.lt.nj)then
1543 goto 1
1544 else
1545 if(ish.ge.3)then
1546 write (ifch,*)'psjarr: nptl',nptl
1547 do k=nptl0,nptl
1548 write(ifch,'(a,i4)')' particle',k
1549 write(ifch,'(i5,2x,5e10.3)')idptl(k)
1550 * ,(pptl(j,k),j=1,5)
1551 enddo
1552 endif
1553 jfl=1
1554 goto999
1555 endif
1556 else
1557 jfirst=0
1558 njpar=ij
1559 ij=ncj(1,ij)
1560 goto 3
1561 endif
1562
1563 else
1564 if(ncj(1,ij).eq.njpar)then
1565 njdau=ncj(2,ij)
1566 else
1567 njdau=ncj(1,ij)
1568 endif
1569 njpar=ij
1570 ij=njdau
1571 goto 3
1572 endif
1573 999 continue
1574 end
1575
1576c------------------------------------------------------------------------
1577 subroutine pspawr(kj)
1578c-----------------------------------------------------------------------
1579c pspawr - writing final parton kj into particle list
1580c------------------------------------------------------------------------
1581c Input:
1582
1583 parameter (mjstr=20000)
1584 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1585 common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr)
1586
1587c eqj(1:4,k) - 4-momentum (qgs) for k-th parton;
1588c bxj(1:4,k) - coordinates for k-th parton formation point;
1589c iqj(k) - ID (qgs) for k-th parton;
1590c ncj(1:2,k) - colour connected partons indexes for k-th parton;
1591c nj - number of partons
1592c-----------------------------------------------------------------------
1593 include 'epos.inc'
1594 include 'epos.incsem'
1595 common/ciptl/iptl
1596
1597 if(ish.ge.9)write (ifch,*)'nptl,kj (sto)',nptl,kj
1598 if(nptl.ge.mxptl)then
1599 write (ifmt,*)'nptl,kj',nptl,kj
1600 call alist('Error in pspawr: nptl out of bounds &',1,nptl)
1601 stop'nptl out of bounds'
1602 endif
1603
1604 if(ifrptl(1,iptl).eq.0)ifrptl(1,iptl)=nptl+1
1605 ifrptl(2,iptl)=nptl+1
1606
1607 nptl=nptl+1
1608 pptl(1,nptl)=eqj(3,kj)
1609 pptl(2,nptl)=eqj(4,kj)
1610 pptl(3,nptl)=eqj(2,kj)
1611 pptl(4,nptl)=eqj(1,kj)
1612 pptl(5,nptl)=0.
1613 idptl(nptl)=psidd(iqj(kj))
1614 iorptl(nptl)=iorj(kj)
1615 jorptl(nptl)=0
1616 istptl(nptl)=20
1617 do i=1,4
1618 xorptl(i,nptl)=bxj(i,kj)
1619 enddo
1620 tivptl(1,nptl)=bxj(5,kj)
1621 tivptl(2,nptl)=bxj(6,kj)
1622 ityptl(nptl)=ityj(kj)
1623 !kkkkkkkkkkkkkkkkkkkkkkk
1624 ! write(*,'(a,2i4,i6,f8.3)')'.... ',kj,nptl,idptl(nptl)
1625 !* ,sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2)
1626 return
1627 end
1628
1629c------------------------------------------------------------------------
1630 subroutine psabor(si,qq,iqc,ncc,ept,jdis,iptl,coordo)
1631c------------------------------------------------------------------------
1632c psabor - highest virtuality subprocess in the ladder
1633c si - c.m. energy squared for the process,
1634c qq - p_t-cutoff for the process due to the evolution,
1635c iqc(i), i=1,2 - incoming parton types(0-g,(+-)1-u(u~),(+-)2-d(d~)etc.),
1636c ncc(i,j), i,j=1,2 - incoming partons color connections,
1637c ept(4) - total 4-momentum for the system of the 2 partons
1638c jdis=0 - hadronic process in pp; 1 - resolved photon process
1639c------------------------------------------------------------------------
1640 double precision ept(4),psutz,psuds
1641 dimension ep3(4),ey(3),wsub(5),iqc(2),ncc(2,2),nqc(2),coordo(6)
1642 parameter (mjstr=20000)
1643 common /psar29/eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
1644 parameter (ntim=1000)
1645 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
1646 &,idaprt(2,ntim)
1647 common/cprtx/nprtjx,pprtx(5,2)
1648 include 'epos.incsem'
1649 include 'epos.inc'
1650 call utpri('psabor',ish,ishini,5)
1651
1652 if(iabs(iqc(1)).ne.4)then !gluon or light quark
1653 q2mass=0.
1654 else !c-quark
1655 q2mass=qcmass**2
1656 endif
1657 p1=si/(1.+q2mass/si)
1658 if(p1.gt.4.*qq)then !|t|-cutoff (p^2_t>qq)
1659 tmin=2.*qq/(1.+sqrt(1.-4.*qq/p1))
1660 else
1661 tmin=2.*qq
1662 endif
1663 tmax=.5*p1
1664
1665 fborn=0.
1666 qt2=tmin*(1.d0-tmin/p1)
1667 if(qt2.lt..999d0*qq.and.ish.ge.2)write(ifmt,*)'qt20,qq',qt2,qq
1668 if(iqc(1).ne.0.or.iqc(2).ne.0)then
1669 do l=1,5 !sum over different subprocesses
1670 wsub(l)=psbori(si,tmin,iqc(1),iqc(2),l) !matrix element
1671 if(l.le.3)then
1672 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1673 elseif(l.le.4)then
1674 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1675 else
1676 wsub(l)=wsub(l)*(alfe/2/pi)**2
1677 endif
1678 fborn=fborn+wsub(l)
1679 enddo
1680 fborn=tmin**2*fborn
1681 else
1682 do l=1,5
1683 wsub(l)=psbori(si,.5*si,iqc(1),iqc(2),l)
1684 if(l.le.3)then
1685 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1686 elseif(l.le.4)then
1687 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1688 else
1689 wsub(l)=wsub(l)*(alfe/2/pi)**2
1690 endif
1691 fborn=fborn+wsub(l)
1692 enddo
1693 fborn=.25*si**2*fborn
1694 endif
1695 if(jdis.eq.0)then
1696 scale1=qt2
1697 else
1698 scale1=4.*qt2
1699 endif
1700 gb0=dble(fborn)*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))
1701 if(ish.ge.7)write(ifch,*)'tmin,gb0:',tmin,gb0
1702
1703c------------------------------------------------
1704c 4-momentum transfer squared is simulated first as dq_t**2/q_t**4 from
1705c tmin to s/2
170614 q2=tmin/(1.-rangen()*(1.-tmin/tmax)) !q2=min(|t|,|u|)
1707 qt2=q2*(1.-q2/p1) !qt2=p_t^2 for the process
1708 if(qt2.lt.qq.and.ish.ge.2)write(ifmt,*)'qt2,qq',qt2,qq
1709 if(rangen().lt..5)then !|u|=q2, |t|=p1-q2
1710 jm=2 !first parton to be considered
1711 tq=p1-q2
1712 else !|t|=q2, |u|=p1-q2
1713 jm=1 !first parton to be considered
1714 tq=q2
1715 endif
1716
1717 fborn=0.
1718 do l=1,5 !sum over different subprocesses
1719 wsub(l)=psbori(si,tq,iqc(1),iqc(2),l)
1720 if(l.le.3)then
1721 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)**2
1722 elseif(l.le.4)then
1723 wsub(l)=wsub(l)*pssalf(qt2/qcdlam)*alfe/2/pi
1724 else
1725 wsub(l)=wsub(l)*(alfe/2/pi)**2
1726 endif
1727 fborn=fborn+wsub(l)
1728 enddo
1729 if(jdis.eq.0)then
1730 scale1=qt2
1731 else
1732 scale1=4.*qt2
1733 endif
1734 gb=dble(q2**2*fborn)
1735 &*psuds(scale1,iqc(1))*psuds(qt2,iqc(2))/gb0 !rejection function
1736 if(ish.ge.7)write(ifch,*)'q2,qt2,gb:',q2,qt2,gb
1737
1738 if(dble(rangen()).gt.gb)goto 14 !rejection
1739
1740c determination of the color configuration
1741 nqc(2)=0
1742 if(iqc(1).eq.0.and.iqc(2).eq.0)then !g+g
1743 jq=int(1.5+rangen()) !jq=1(2) - transfer of color (anticolor)
1744 nqc(1)=ncc(jq,jm)
1745
1746 if(rangen().lt.wsub(1)/fborn)then !gg->gg
1747 if(rangen().lt..5)then
1748 jt=1 !anticolor-color annihilation
1749 nqc(2)=0
1750 njc1=ncc(3-jq,jm)
1751 njc2=ncc(jq,3-jm)
1752 if(iqj(njc1).ne.0)then
1753 ncj(1,njc1)=njc2
1754 else
1755 ncj(jq,njc1)=njc2
1756 endif
1757 if(iqj(njc2).ne.0)then
1758 ncj(1,njc2)=njc1
1759 else
1760 ncj(3-jq,njc2)=njc1
1761 endif
1762 else
1763 jt=2 !produced gluons get color and
1764 nqc(2)=ncc(3-jq,3-jm) !anticolor from the 2 parents
1765 endif
1766 else !gg->qq~
1767 jt=9 !anticolor-color annihilation
1768 iqc(jm)=int(naflav*rangen()+1)*(3-2*jq) !(anti)quark flavor
1769 iqc(3-jm)=-iqc(jm)
1770 njc1=ncc(3-jq,jm)
1771 njc2=ncc(jq,3-jm)
1772 if(iqj(njc1).ne.0)then
1773 ncj(1,njc1)=njc2
1774 else
1775 ncj(jq,njc1)=njc2
1776 endif
1777 if(iqj(njc2).ne.0)then
1778 ncj(1,njc2)=njc1
1779 else
1780 ncj(3-jq,njc2)=njc1
1781 endif
1782 endif
1783
1784 elseif(iqc(1)*iqc(2).eq.0)then !q(q~)+g
1785 if(iqc(1)+iqc(2).gt.0)then
1786 jq=1 !q
1787 else
1788 jq=2 !q~
1789 endif
1790 if(rangen().lt.wsub(1)/fborn)then !q(q~)g->q(q~)g
1791 if(rangen().lt..5)then !anticolor-color annihilation
1792 if(iqc(jm).eq.0)then
1793 jt=3 !first parton=g
1794 nqc(1)=ncc(jq,jm)
1795 njc1=ncc(3-jq,jm)
1796 njc2=ncc(1,3-jm)
1797 if(iqj(njc1).ne.0)then
1798 ncj(1,njc1)=njc2
1799 else
1800 ncj(jq,njc1)=njc2
1801 endif
1802 if(iqj(njc2).ne.0)then
1803 ncj(1,njc2)=njc1
1804 else
1805 ncj(3-jq,njc2)=njc1
1806 endif
1807
1808 else
1809 jt=4 !first parton=q(q~)
1810 nqc(1)=0
1811 njc1=ncc(1,jm)
1812 njc2=ncc(3-jq,3-jm)
1813 if(iqj(njc1).ne.0)then
1814 ncj(1,njc1)=njc2
1815 else
1816 ncj(3-jq,njc1)=njc2
1817 endif
1818 if(iqj(njc2).ne.0)then
1819 ncj(1,njc2)=njc1
1820 else
1821 ncj(jq,njc2)=njc1
1822 endif
1823 endif
1824
1825 else !color transfer
1826 if(iqc(jm).eq.0)then
1827 jt=5 !first parton=g
1828 nqc(2)=ncc(3-jq,jm)
1829 nqc(1)=ncc(1,3-jm)
1830 else !first parton=q(q~)
1831 jt=6
1832 nqc(1)=ncc(jq,3-jm)
1833 endif
1834 endif
1835
1836 else !q(q~)g->q(q~)-gamma (+-color annihilation)
1837 if(iqc(jm).eq.0)then
1838 jt=11 !first parton=g
1839 nqc(1)=ncc(jq,jm)
1840 njc1=ncc(3-jq,jm)
1841 njc2=ncc(1,3-jm)
1842 if(iqj(njc1).ne.0)then
1843 ncj(1,njc1)=njc2
1844 else
1845 ncj(jq,njc1)=njc2
1846 endif
1847 if(iqj(njc2).ne.0)then
1848 ncj(1,njc2)=njc1
1849 else
1850 ncj(3-jq,njc2)=njc1
1851 endif
1852 iqc(jm)=iqc(3-jm)
1853 iqc(3-jm)=10 !make the second output is gamma.
1854
1855 else
1856 jt=12 !first parton=q(q~)
1857 nqc(1)=ncc(jq,3-jm) !here nqc(1) is gluon.
1858 njc1=ncc(1,jm)
1859 njc2=ncc(3-jq,3-jm)
1860 if(iqj(njc1).ne.0)then
1861 ncj(1,njc1)=njc2
1862 else
1863 ncj(3-jq,njc1)=njc2
1864 endif
1865 if(iqj(njc2).ne.0)then
1866 ncj(1,njc2)=njc1
1867 else
1868 ncj(jq,njc2)=njc1
1869 endif
1870 iqc(3-jm)=10
1871 endif
1872 endif
1873
1874 elseif(iqc(1)*iqc(2).gt.0)then
1875 jt=7 !qq->qq (q~q~->q~q~)
1876 if(iqc(1).gt.0)then
1877 jq=1
1878 else
1879 jq=2
1880 endif
1881 nqc(1)=ncc(1,3-jm)
1882
1883 else ! qq~ ->
1884 if(iqc(jm).gt.0)then
1885 jq=1
1886 else
1887 jq=2
1888 endif
1889 aks=rangen()
1890 if(aks.lt.(wsub(1)+wsub(2))/fborn)then
1891 jt=8 ! qq~->qq~ (anticolor-color annihilation)
1892 if(aks.gt.wsub(1)/fborn)then
1893 iqa=iabs(iqc(jm))
1894 iq=int((naflav-1)*rangen())+1
1895 if(iq.eq.iqa)iq=naflav
1896 iqc(jm)=iq*iqc(jm)/iqa
1897 iqc(3-jm)=-iqc(jm)
1898 endif
1899 nqc(1)=0
1900 njc1=ncc(1,jm)
1901 njc2=ncc(1,3-jm)
1902 if(iqj(njc1).ne.0)then
1903 ncj(1,njc1)=njc2
1904 else
1905 ncj(3-jq,njc1)=njc2
1906 endif
1907 if(iqj(njc2).ne.0)then
1908 ncj(1,njc2)=njc1
1909 else
1910 ncj(jq,njc2)=njc1
1911 endif
1912 elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3))/fborn)then
1913 jt=10 !color transfer qq~->gg
1914 iqc(1)=0
1915 iqc(2)=0
1916 nqc(1)=ncc(1,jm)
1917 nqc(2)=0
1918 elseif(aks.lt.(wsub(1)+wsub(2)+wsub(3)+wsub(4))/fborn)then
1919 jt=13 ! qq~->g+gamma
1920 nqc(1)=ncc(1,jm)
1921 nqc(2)=ncc(1,3-jm)
1922 iqc(jm)=0
1923 iqc(3-jm)=10
1924 else
1925 jt=14 ! qq~->gamma+gamma
1926 njc1=ncc(1,jm)
1927 njc2=ncc(1,3-jm)
1928 if(iqj(njc1).ne.0)then
1929 ncj(1,njc1)=njc2
1930 else
1931 ncj(jq,njc1)=njc2
1932 endif
1933 if(iqj(njc2).ne.0)then
1934 ncj(1,njc2)=njc1
1935 else
1936 ncj(3-jq,njc2)=njc1
1937 endif
1938 iqc(jm)=10
1939 iqc(3-jm)=10
1940 endif
1941 endif
1942
1943 if(jt.ne.8.and.jt.ne.9)then
1944 jq2=jq
1945 else
1946 jq2=3-jq
1947 endif
1948
1949 call psdeftr(si+q2mass,ept,ey) !lorentz boost to c.m. frame
1950
1951 qt=sqrt(qt2) !p_t
1952 call pscs(bcos,bsin) !cos and sin of the polar angle
1953 if(iabs(iqc(1)).ne.4)then
1954clight cone momentum share for the first parton
1955 z=sngl(psutz(dble(si),dble(qt2),dble(qt2)))
1956 if((jt.eq.11.and.jm.eq.1).or.(jt.eq.12.and.jm.eq.2)
1957 $ .or.(jt.eq.13.and.jm.eq.2))z=1-z
1958 wp3=z*sqrt(si)
1959 wm3=qt2/wp3
1960 elseif(jm.eq.1)then
1961 z=sngl(psutz(dble(si),dble(qt2+q2mass),dble(qt2)))
1962 wp3=z*sqrt(si)
1963 wm3=(qt2+q2mass)/wp3
1964 else
1965 z=sngl(psutz(dble(si),dble(qt2),dble(qt2+dble(q2mass))))
1966 wp3=z*sqrt(si)
1967 wm3=qt2/wp3
1968 endif
1969 ep3(1)=.5*(wp3+wm3) !parton 4-momentum
1970 ep3(2)=.5*(wp3-wm3)
1971 ep3(3)=qt*bcos
1972 ep3(4)=qt*bsin
1973 call psdefrot(ep3,s0xh,c0xh,s0h,c0h) !spacial rotation to z-axis
1974 if(iqc(jm).eq.0)then
1975 iq2ini1=9
1976 else
1977 iq2ini1=iqc(jm)
1978 endif
1979 if(iqc(3-jm).eq.0)then
1980 iq2ini2=9
1981 else
1982 iq2ini2=iqc(3-jm)
1983 endif
1984 if(jt.le.10)then
1985 qq1=qt2
1986 qq2=qt2
1987 elseif(jt.le.13)then
1988 qq1=qt2
1989 qq2=0
1990 else
1991 qq1=0
1992 qq2=0
1993 endif
1994
1995 call timsh2(qq1,qq2,sqrt(si),iq2ini1,iq2ini2) !final state cascade
1996
1997 if(jt.le.10)then !color connection for the 2nd parton
1998 if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
1999 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2000 if(jt.eq.1)then
2001 nqc(1)=nqc(2)
2002 nqc(2)=ncc(3-jq,3-jm)
2003 elseif(jt.eq.2)then
2004 nqc(2)=ncc(3-jq,jm)
2005 nqc(1)=ncc(jq,3-jm)
2006 elseif(jt.eq.3)then
2007 nqc(1)=nqc(2)
2008 elseif(jt.eq.4)then
2009 nqc(2)=nqc(1)
2010 nqc(1)=ncc(jq,3-jm)
2011 elseif(jt.eq.5)then
2012 nqc(1)=ncc(jq,jm)
2013 elseif(jt.eq.6)then
2014 nqc(2)=ncc(3-jq,3-jm)
2015 nqc(1)=ncc(1,jm)
2016 elseif(jt.eq.7)then
2017 nqc(1)=ncc(1,jm)
2018 elseif(jt.eq.9)then
2019 nqc(1)=ncc(3-jq,3-jm)
2020 elseif(jt.eq.10)then
2021 nqc(1)=nqc(2)
2022 nqc(2)=ncc(1,3-jm)
2023 endif
2024 if(ish.ge.6)write(ifch,*)'jt,jq2,nqc:',jt,jq2,nqc
2025 call psreti(nqc,jq2,2,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstr.
2026 elseif(jt.le.13)then
2027 if(ish.ge.6)write(ifch,*)'jt,jq,nqc:',jt,jq,nqc
2028 call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h) !color conn. reconstruction
2029 ep3(1)=pprt(4,2)
2030 ep3(2)=pprt(3,2)
2031 ep3(3)=pprt(1,2)
2032 ep3(4)=pprt(2,2)
2033 call psrotat(ep3,s0xh,c0xh,s0h,c0h) !special rotation for photon.
2034 call pstrans(ep3,ey,1)
2035 nptl=nptl+1
2036 pptl(1,nptl)=ep3(3)
2037 pptl(2,nptl)=ep3(4)
2038 pptl(3,nptl)=ep3(2)
2039 pptl(4,nptl)=ep3(1)
2040 pptl(5,nptl)=0
2041 idptl(nptl)=10
2042 iorptl(nptl)=iptl
2043 istptl(nptl)=0
2044 jorptl(nptl)=0
2045 do i=1,4
2046 xorptl(i,nptl)=coordo(i)
2047 enddo
2048 tivptl(1,nptl)=coordo(5)
2049 tivptl(2,nptl)=coordo(6)
2050 ityptl(nptl)=71
2051 ifrptl(1,nptl)=0
2052 ifrptl(2,nptl)=0
2053 else
2054 do j=1,2
2055 ep3(1)=pprt(4,j)
2056 ep3(2)=pprt(3,j)
2057 ep3(3)=pprt(1,j)
2058 ep3(4)=pprt(2,j)
2059 call psrotat(ep3,s0xh,c0xh,s0h,c0h) !special rotation for photon.
2060 call pstrans(ep3,ey,1)
2061 nptl=nptl+1
2062 pptl(1,nptl)=ep3(3)
2063 pptl(2,nptl)=ep3(4)
2064 pptl(3,nptl)=ep3(2)
2065 pptl(4,nptl)=ep3(1)
2066 pptl(5,nptl)=0
2067 idptl(nptl)=10
2068 iorptl(nptl)=iptl
2069 istptl(nptl)=0
2070 jorptl(nptl)=0
2071 do i=1,4
2072 xorptl(i,nptl)=coordo(i)
2073 enddo
2074 tivptl(1,nptl)=coordo(5)
2075 tivptl(2,nptl)=coordo(6)
2076 ityptl(nptl)=72
2077 ifrptl(1,nptl)=0
2078 ifrptl(2,nptl)=0
2079 enddo
2080 endif
2081 call utprix('psabor',ish,ishini,5)
2082 return
2083 end
2084
2085c------------------------------------------------------------------------
2086 subroutine psreti(nqc,jort,nfprt,ey,s0xh,c0xh,s0h,c0h)
2087c-----------------------------------------------------------------------
2088c jet reconstructuring procedure - 4-momenta for all final jets
2089c nqc(i) - colour connections for the jet
2090c jort - color orientation for gluons (=1 if +color goes first, =-1 otherwise)
2091
2092 parameter (ntim=1000)
2093 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
2094 &,idaprt(2,ntim)
2095c nprtj - number of partons in the jet (including virtual ones)
2096c pprt - 5-momenta for the partons
2097c idprt - parton id
2098c iorprt - parent parton position in the list
2099c idaprt - daughter partons positions in the list
2100
2101c output:
2102 parameter (mjstr=20000)
2103 common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),nj
2104c nj - number of final jets
2105c eqj(i,j) - 4-momentum for the final jet j
2106c iqj(j) - flavour for the final jet j
2107c ncj(m,j) - colour connections for the final jet j
2108c-----------------------------------------------------------------------
2109 dimension ep3(4),nqc(2),ncc(2,ntim),ey(3)
2110 include 'epos.inc'
2111 include 'epos.incsem'
2112 common/cprtx/nprtjx,pprtx(5,2)
2113
2114 if(ish.ge.6)then
2115 write (ifch,*)'nprtj',nprtj
2116 do i=1,nprtj
2117 write (ifch,*)'i,ic,np,ndd',i,idprt(i),iorprt(i),
2118 * idaprt(1,i),idaprt(2,i)
2119 enddo
2120 endif
2121
2122 ncc(1,nfprt)=nqc(1)
2123 if(idprt(nfprt).eq.9)ncc(2,nfprt)=nqc(2)
2124 iprt=nfprt
2125
2126 if(nprtjx.eq.2)then !out Born before timelike cascade
2127 ep3(1)=pprtx(4,iprt)
2128 ep3(2)=pprtx(3,iprt)
2129 ep3(3)=pprtx(1,iprt)
2130 ep3(4)=pprtx(2,iprt)
2131 call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2132 call pstrans(ep3,ey,1)
2133 pprtx(4,iprt)=ep3(1)
2134 pprtx(3,iprt)=ep3(2)
2135 pprtx(1,iprt)=ep3(3)
2136 pprtx(2,iprt)=ep3(4)
2137 endif
2138
21391 continue
2140
2141 idau1=idaprt(1,iprt)
2142 idau2=idaprt(2,iprt)
2143 icp=idprt(iprt)
2144
2145 if(ish.ge.6)then
2146 write (ifch,*)'1-iprt,icp,idau1,idau2',iprt,icp,idau1,idau2,
2147 * ncc(1,iprt)
2148 if(icp.eq.9)write (ifch,*)'ncc2',ncc(2,iprt)
2149 endif
2150
2151 if(idau1.ne.0.)then !virtual parton
2152 icd1=idprt(idau1)
2153
2154 if(icp.eq.9)then
2155 if(icd1.ne.9)then !g -> qq~
2156 ncc(1,idau1)=ncc(jort,iprt)
2157 ncc(1,idau2)=ncc(3-jort,iprt)
2158 else !g -> gg
2159 ncc(1,idau1)=ncc(1,iprt)
2160 ncc(2,idau1)=0
2161 ncc(2,idau2)=ncc(2,iprt)
2162 ncc(1,idau2)=0
2163 endif
2164 else !q -> qg
2165 ncc(1,idau1)=0
2166 if(icp*(3-2*jort).gt.0)then
2167 ncc(1,idau2)=ncc(1,iprt)
2168 ncc(2,idau2)=0
2169 else
2170 ncc(1,idau2)=0
2171 ncc(2,idau2)=ncc(1,iprt)
2172 endif
2173 endif
2174 iprt=idau1
2175 goto 1
2176 else
2177
2178 nj=nj+1
2179 ep3(1)=pprt(4,iprt)
2180 ep3(2)=pprt(3,iprt)
2181 ep3(3)=pprt(1,iprt)
2182 ep3(4)=pprt(2,iprt)
2183 call psrotat(ep3,s0xh,c0xh,s0h,c0h)
2184 call pstrans(ep3,ey,1)
2185 do i=1,4
2186 eqj(i,nj)=ep3(i)
2187 enddo
2188
2189 if(icp.eq.9)then
2190 iqj(nj)=0
2191 elseif(iabs(icp).lt.3)then
2192 iqj(nj)=icp
2193 elseif(iabs(icp).eq.3)then
2194 iqj(nj)=icp*4/3
2195 else
2196 iqj(nj)=icp*10
2197 endif
2198
2199 if(iqj(nj).ne.0)then
2200 njc=ncc(1,iprt)
2201 if(njc.ne.0)then
2202 ncj(1,nj)=njc
2203 iqc=iqj(njc)
2204 if(iqc.ne.0)then
2205 ncj(1,njc)=nj
2206 else
2207 if(iqj(nj).gt.0)then
2208 ncj(2,njc)=nj
2209 else
2210 ncj(1,njc)=nj
2211 endif
2212 endif
2213 else
2214 ncc(1,iprt)=nj
2215 endif
2216 else
2217
2218 do m=1,2
2219 if(jort.eq.1)then
2220 m1=m
2221 else
2222 m1=3-m
2223 endif
2224 njc=ncc(m1,iprt)
2225 if(njc.ne.0)then
2226 ncj(m,nj)=njc
2227 iqc=iqj(njc)
2228 if(iqc.ne.0)then
2229 ncj(1,njc)=nj
2230 else
2231 ncj(3-m,njc)=nj
2232 endif
2233 else
2234 ncc(m1,iprt)=nj
2235 endif
2236 enddo
2237 endif
2238 if(ish.ge.6)then
2239 write (ifch,*)'jet-nj,iprt,icp,iqj(nj),ncj',
2240 * nj,iprt,icp,iqj(nj),ncj(1,nj)
2241 if(icp.eq.9)write (ifch,*)'ncj2',ncj(2,nj)
2242 endif
2243 endif
2244
22452 continue
2246 if(iprt.ne.nfprt)then
2247 icp=idprt(iprt)
2248 ipar=iorprt(iprt)
2249 idau1=idaprt(1,ipar)
2250 idau2=idaprt(2,ipar)
2251 if(ish.ge.6)then
2252 write (ifch,*)'2-iprt,icp,idau1,idau2,ipar',
2253 * iprt,icp,idau1,idau2,ipar,ncc(1,iprt)
2254 if(icp.eq.9)write (ifch,*)ncc(2,iprt)
2255 endif
2256
2257 if(idau1.eq.iprt)then
2258 if(icp.eq.9)then !g -> gg
2259 ncc(1,ipar)=ncc(1,iprt)
2260 ncc(1,idau2)=ncc(2,iprt)
2261 else
2262 icpar=idprt(ipar)
2263 if(icpar.eq.9)then !g -> qq~
2264 ncc(jort,ipar)=ncc(1,iprt)
2265 else !q -> qg
2266 if(icp*(3-2*jort).gt.0)then
2267 ncc(2,idau2)=ncc(1,iprt)
2268 else
2269 ncc(1,idau2)=ncc(1,iprt)
2270 endif
2271 endif
2272 endif
2273 iprt=idau2
2274 goto 1
2275
2276 else
2277 if(icp.eq.9)then
2278 icpar=idprt(ipar)
2279 if(icpar.eq.9)then !g -> gg
2280 ncc(2,ipar)=ncc(2,iprt)
2281 ncc(2,idau1)=ncc(1,iprt)
2282 else !q -> qg
2283 if(icpar*(3-2*jort).gt.0)then
2284 ncc(1,ipar)=ncc(1,iprt)
2285 ncc(1,idau1)=ncc(2,iprt)
2286 else
2287 ncc(1,ipar)=ncc(2,iprt)
2288 ncc(1,idau1)=ncc(1,iprt)
2289 endif
2290 endif
2291 else
2292 ncc(3-jort,ipar)=ncc(1,iprt)
2293 endif
2294 iprt=ipar
2295 goto 2
2296 endif
2297 else
2298 if(ish.ge.6)write (ifch,*)'3-iprt,ncc',iprt,ncc(1,iprt)
2299 nqc(1)=ncc(1,nfprt)
2300 if(idprt(nfprt).eq.9)nqc(2)=ncc(2,nfprt)
2301 endif
2302 return
2303 end
2304