]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-fra-163.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update timestamp for new data points simulation
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-fra-163.f
CommitLineData
9ef1c2d9 1c-----------------------------------------------------------------------
2 subroutine gakfra(iret)
3c-----------------------------------------------------------------------
4
5 include 'epos.inc'
6 include 'epos.incems'
7 include 'epos.incsem'
8 parameter (eta=0.4,etap=0.1)
9 parameter (mxnbr=500,mxnba=5000)
10 common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
11 $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
12 & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
13c common /cpptr/ pptr(4,mxptl),ystr(mxptl)
14 double precision p1,p12,p2
15 dimension co(12),ind(0:mxnbr),p1(5),p2(5)
16 dimension ic2(2),ic1(2),ic(2)
17 common/pb/pb /cnsbp/nsbp /cn8ptl/n8ptl
18 common/czz/kky,krm,zzod,zzop,kdrm,zzoq,zzos,zzrm,zzipp,zzipt
19 &,zzid,zzip,zzzzr,itp
20 double precision psg
21 common/zpsg/psg(5)
22 logical leadstr,go
23
24 pf(i,j)=(pst(4,i)-pst(3,i))*(pst(4,j)+pst(3,j))
25 & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
26
27 nob=0
28 call utpri('gakfra',ish,ishini,4)
29 iret=0
30c.....search string in pptl-list...................................
31 nkmax=nptl
32 nk1=maproj+matarg+1
33 2 continue
34 if(nclean.gt.0.and.nptl.gt.mxptl/2)then
35 nnnpt=0
36 do iii=maproj+matarg+1,nptl
37 go=.true.
38 if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
39 if(go.and.mod(istptl(iii),10).ne.0)then
40 istptl(iii)=99
41 nnnpt=nnnpt+1
42 endif
43 enddo
44 if(nnnpt.gt.mxptl-nptl)then
45 nptl0=nptl
46 call utclea(maproj+matarg+1,nptl0)
47 nkmax=nptl
48 nk1=maproj+matarg+1
49c write(ifch,'(a)')' ---------after utclea------------'
50c call blist('&',nk1,nkmax)
51 endif
52 endif
53
54 do while(nk1.le.nkmax)
55 if(istptl(nk1).eq.20)then
56 nk2=nk1+1
57 do while(idptl(nk2).eq.9)
58 nk2=nk2+1
59 enddo
60 goto 3 !ok, string from nk1 to nk2
61 else
62 nk1=nk1+1
63 endif
64 enddo
65
66 431 continue
67
68
69 goto 9999 !no more string
70
71
72c.......................................................................
73c decay string nk1 - nk2
74c.......................................................................
75
76 3 nob=-1
77 if(ish.ge.3)then
78 write(ifch,'(/1x,25a1/5x,a/1x,25a1/)')
79 * ('-',k=1,25),'string decay',('-',k=1,25)
80 write(ifch,'(a)')' ---------string------------'
81 call blist('&',nk1,nk2)
82 endif
83
84 kky=0
85 if(nk2-nk1.gt.1.or.iappl.eq.6)kky=1
86 itp=ityptl(nk1)
87 krm=0
88 if(itp.ge.50.and.itp.le.59.or.itp.ge.40.and.itp.le.49)krm=1
89 kdrm=0
90 if(itp.eq.43.or.itp.eq.53)kdrm=1
91c special behaviour of remnant of reactions with lost Pomerons
92 if(itp.eq.47.or.itp.eq.57)then
93 kdrm=1
94 krm=0
95 endif
96 zzzz=1.
97 zzzzp=1.
98 zzzzt=1.
99 zzod=1.
100 zzos=1.
101 zzop=1.
102 zzoq=1.
103 zzrm=1.
104 zzipp=1.
105 zzipt=1.
106 zzid=1.
107 zzip=1.
108 if(isplit.eq.1)then
109 zzzz=zpaptl(nk1)+zpaptl(nk2)
110 zzzzp=zpaptl(nk1)
111 zzzzt=zpaptl(nk2)
112 zzid=min( 1+zidmax , 1+zidinc *zzzz ) !<-----
113 zzip=min( 1+zopmax , 1+zopinc *zzzz ) !<-----
114c zzrm=min( zrmmax , zrminc *zzzz ) !<-----
115 if(itp.ge.40.and.itp.le.49)then
116 zzzzr=zzzzp
117 elseif(itp.ge.50.and.itp.le.59)then
118 zzzzr=zzzzt
119 endif
120 zzod=min( 1+zodmax , 1+zodinc *zzzzr ) !<-----
121 zzos=min( 1+zosmax , 1+zosinc *zzzzr ) !<-----
122 zzop=min( 1+zopmax , 1+zopinc *zzzzr ) !<-----
123 zzoq=min( 1+2.*zopmax , 1+ 2.*zopinc *zzzzr ) !<-----
124 zzipp=min( 1+zipmax , 1+zipinc *zzzzp ) !<-----
125 zzipt=min( 1+zipmax , 1+zipinc *zzzzt ) !<-----
126 endif
127 do k=1,5
128 p1(k)=0d0
129 enddo
130 do i=nk1,nk2
131 do k=1,5
132 p2(k)=dble(pptl(k,i))
133 enddo
134 p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
135 do k=1,4
136 p1(k)=p1(k)+p2(k)
137 enddo
138 enddo
139 p1(5)=(p1(4)-p1(3))*(p1(4)+p1(3))-p1(2)**2-p1(1)**2
140 if(p1(5).gt.0.d0)then
141 p1(5)=sqrt(p1(5))
142 else
143 iorrem=iorptl(nk1)
144 write(*,*)'Precision problem in gakfra, p:',
145 & p1(5),p1(4),p1(3),p1(2),p1(1)
146 write(*,*)'origin : ',iorrem
147 p1(5)=0.d0
148c if(iorrem.ne.0.and.
149c & (ityptl(iorrem).eq.40.or.ityptl(iorrem).eq.10))
150c & p1(5)=dble(pptl(5,iorrem))
151 if(iorrem.ne.0)then
152 write(*,*)'string mass : ',ityptl(iorrem),pptl(5,iorrem)
153 p1(5)=dble(pptl(5,iorrem))
154 endif
155 endif
156 do k=1,5
157 psg(k)=p1(k)
158 enddo
159
160 if(iappl.ne.6)qsqptl(nk1)=p1(5)*p1(5) !hadronic string energy
161 !(already defined for lepton)
162
163c.....light like ...............................................
164
165 if(ish.ge.6) call blist("before light like&",nk1,nk2)
166 if(rangen().lt.0.5)then
167 i1=nk1
168 i2=nk1+1
169 else
170 i1=nk2-1
171 i2=nk2
172 endif
173 ii=1
174 do while(ii.ge.0)
175 do i=1,4
176 p2(i)=dble(pptl(i,i1))+dble(pptl(i,i2))
177 enddo
178c p2(5)=sqrt(p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2)
179 am1=pptl(5,i1)**2
180 am2=pptl(5,i2)**2
181 am12=max(0d0,p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**2-am1-am2)/2
182 if(am12**2.gt.am1*am2)then
183 ak1=.5*((am2+am12)/sqrt(am12**2-am1*am2))-.5
184 ak2=.5*((am1+am12)/sqrt(am12**2-am1*am2))-.5
185 else
186 ak1=0.
187 ak2=0.
188 endif
189 if(ish.ge.6)write(ifch,'(a,2i4,9f12.6)') 'overlaps'
190 $ ,i1,i2,ak1,ak2,sqrt(2.*am12)
191 & ,am1,am2
192 do j=1,4
193 a1=(1.+ak1)*pptl(j,i1)-ak2*pptl(j,i2)
194 a2=(1.+ak2)*pptl(j,i2)-ak1*pptl(j,i1)
195 pptl(j,i1)=a1
196 pptl(j,i2)=a2
197 enddo
198 pptl(5,i1)=0.
199 pptl(5,i2)=0.
200 if(nk2-nk1.eq.1) ii=-1
201 ii=ii-1
202 if (nk1.eq.i1)then
203 i1=nk2-1
204 i2=nk2
205 else
206 i1=nk1
207 i2=nk1+1
208 endif
209 enddo
210 if(ish.ge.6) call blist("after light like&",nk1,nk2)
211
212c.....on-shell and copy ...............................................
213
214 if(ish.ge.6) write (ifch,*) 'on-shell check'
215 do i=nk1,nk2
216 do k=1,5
217 p2(k)=dble(pptl(k,i))
218 enddo
219 p12=p2(4)
220 p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
221 if(abs(p12-p2(4))/max(.1d0,p12,p2(4)).ge.1e-1.and.ish.ge.2)then
222 write(ifch,*)'warning: on-shell problem:'
223 & ,i, idptl(i),(pptl(k,i),k=1,4),' (',sngl(p2(4)),') '
224 & ,pptl(5,i), istptl(i)
225 endif
226 if(ish.ge.6) write (ifch,*) p12 ,'->',p2(4)
227 call utlob2(1,p1(1),p1(2),p1(3),p1(4),p1(5)
228 $ ,p2(1),p2(2),p2(3),p2(4),60)
229 f=1.0
230 if(i.ne.nk1.and.i.ne.nk2)f=0.5
231 nmax=1
232 ff=1.
233 aemax=1000. ! here max band mass
234 if(f*p2(4).gt.aemax) then
235 nmax = int(f*p2(4)/aemax)+1
236 ff=1./real(nmax)
237c print *,"nmax",nmax
238 endif
239 nn=1
240 do while(nn.le.nmax)
241 f=.5
242 if(i.eq.nk1.and.nn.eq. 1 ) f=1.
243 if(i.eq.nk2.and.nn.eq.nmax) f=1.
244 nob=nob+1
245 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
246c if(nmax.ge.2) print *, nob,ff,f,i
247 do k=1,4
248 pst(k,nob)=ff*f*p2(k)
249 ipst(nob)=i
250 enddo
251 nn=nn+1
252 enddo
253 enddo !i
254
255 do i1=nob-1,1,-1 ! copy again gluons
256 nob=nob+1
257 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
258 do k=1,4
259 pst(k,nob)=pst(k,i1)
260 ipst(nob)=ipst(i1)
261 enddo
262 enddo
263 nob=nob+1 ! nob is number of kinks
264 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
265 if(ish.ge.6) then
266 write (ifch,*) 'bands:'
267 do i=0,nob-1
268 write (ifch,'(i4,4g20.13)') i,(pst(k,i),k=1,4)
269 enddo
270 endif
271
272c.....total momentum...............................................
273 do k=1,4
274 pst(k,nob)=0.
275 enddo
276 do i=0,nob-1
277 do k=1,4
278 pst(k,nob)=pst(k,nob)+pst(k,i)
279 enddo
280 enddo
281
282
283c.....pbreak ................................
284
285 x=sqrt(qsqptl(nk1))
286 if(krm.eq.1)then
287 if(kdrm.eq.0)then
288 pbi=pbreakr
289 else
290 pbi=pbreakd
291 endif
292 elseif(kky.eq.0)then
293 pbi=pbreak
294 elseif(kky.eq.1)then
295 pbi=pbreakk
296 else
297 stop'in gakfra: unknown kky'
298 endif
299 if(pbi.gt.0.d0)then
300 pb=pbi
301 else
302 if(x.lt.14)then
303c pb=0.38
304 pb=0.4
305 else
306 pb=0.14369+ 854.965/x**3 - 111.1144/x**2 + 7.163/x
307c !! {{14,0.4},{22,0.32},{34,0.28},{91.2,0.21}}
308c pb=0.143+ 900/x**3 - 77/x**2 + 4.5/x
309c pb=0.165+ 1020/x**3 - 80/x**2 + 3.8/x
310c pb=0.17+ 950/x**3 - 76/x**2 + 3.8/x
311c pb=0.165+ 600/x**3 - 50/x**2 + 3.8/x
312c pb=0.135+ 855./x**3 - 111./x**2 + 7./x
313c pb=0.05+ 1100./x**3 - 80./x**2 + 5./x
314 endif
315 pb=pb*abs(pbi)
316 endif
317c write(*,*)'pb',kky,krm,kdrm,x,pb,zzrm,zzzz
318
319
320c.....write total string on pptl - list.............................
321 nptl=nptl+1
322 if(nptl.gt.mxptl)call utstop('gakfra: mxptl too small ...&')
323 call idtr5(idptl(nk1),ic1)
324 call idtr5(idptl(nk2),ic2)
325 ic(1)=ic1(1)+ic2(1)
326 ic(2)=ic1(2)+ic2(2)
327 idptl(nptl) = 8*10**8+ ic(1)*100 + ic(2)/100 !nbr+1
328 istptl(nptl)=10
329 n8ptl=nptl
330 do i=1,5
331 pptl(i,nptl)=p1(i)
332 enddo
333 if(ish.ge.3)then
334 write(ifch,'(a)')' ---------total string------------'
335 call blist('&',nptl,nptl)
336 endif
337
338c....................................................................
339 ijb(1,0)=-1
340 ijb(2,0)=-1
341 ijb(1,1)=-1
342 ijb(2,1)=-1
343 iflb(0)=-idptl(nk1)
344 iflb(1)= idptl(nk2)
345 do i=1,4
346 ptb(i,0)=0.
347 ptb(i,1)=0.
348 enddo
349c calculate isospin to conserve isospin symmetry in popcorn
350c p -> pi+ <=> n -> pi-
351 if(abs(iflb(0)).gt.6)then
352 iq01=mod(abs(iflb(0))/100,10)
353 iq02=abs(iflb(0))/1000
354 iso=(isospin(iq01)+isospin(iq02))*sign(1,-iflb(0))
355 else
356 iso=isospin(abs(iflb(0)))*sign(1,-iflb(0))
357 endif
358 if(abs(iflb(1)).gt.6)then
359 iq11=mod(abs(iflb(1))/100,10)
360 iq12=abs(iflb(1))/1000
361 iso=iso+(isospin(iq11)+isospin(iq12))*sign(1,iflb(1))
362 else
363 iso=iso+isospin(abs(iflb(1)))*sign(1,iflb(1))
364 endif
365
366c...................light string....................................
367 amm=ammin(idptl(nk1),idptl(nk2))
368
369 if(sqrt(max(0.,pf(nob,nob))).lt.amm)then
370 id=idsp( idptl(nk1),idptl(nk2) )
371 if(ish.ge.1)then
372 write (ifch,*)
373 . '-------string too light to decay--------'
374 write (ifch,*) id,p1(5),amm
375 write (ifch,*) id, sqrt(max(0.,pf(nob,nob))) ,amm
376 endif
377 am=sqrt(max(0.,pf(nob,nob)))
378 call idress(id,am,idr,iadj)
379 if(ish.ge.1)write (ifch,*) id,am,idr,iadj
380 nptl=nptl+1
381 if(nptl.gt.mxptl)
382 & call utstop('gakfra (light): mxptl too small ...&')
383 do i=1,5
384 pptl(i,nptl)=p1(i)
385 enddo
386 do i=nk1,nk2
387 istptl(i)=21
388 ifrptl(1,i)=n8ptl
389 ifrptl(2,i)=0
390 enddo
391 istptl(n8ptl)=29 !code for fragmented string
392 ityptl(n8ptl)=ityptl(nk1)
393 itsptl(n8ptl)=itsptl(nk1)
394 rinptl(n8ptl)=-9999
395 iorptl(n8ptl)=nk1
396 jorptl(n8ptl)=nk2
397 ifrptl(1,n8ptl)=n8ptl+1
398 ifrptl(2,n8ptl)=nptl
399 ityptl(nptl)=ityptl(nk1)
400 itsptl(nptl)=itsptl(nk1)
401 rinptl(nptl)=-9999
402 istptl(nptl)=0
403 iorptl(nptl)=n8ptl
404 jorptl(nptl)=0
405 ifrptl(1,nptl)=0
406 ifrptl(2,nptl)=0
407 if(idr.ne.0)then
408 idptl(nptl)=idr
409 else
410 idptl(nptl)=id
411 endif
412 do j=1,4
413 xorptl(j,nptl)=xorptl(j,nk1)
414 enddo
415 tivptl(1,nptl)=xorptl(4,nptl)
416 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
417 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
418 if(abs(p1(5)-am).gt.0.01)then
419c write (*,*) 'light string --- particle off-shell!!!!'
420c write (*,*) idr,' mass:',p1(5),' should be:',am
421 endif
422 goto 98
423 endif
424
425c................search breakpoints...................................
426
427 n=0
428 nsbp=1
429 nbr=0
430 iok=0
431 !----------------------!
432 11 call gaksbp(0,1,iok)
433 !----------------------!
434 nbrtry=0
435 goto 10
436
437c..............delete breakpoint....................................
438 9 if(ish.ge.4)write(ifch,*) 'delete breakpoint',n,' -> '
439 &,nbr-1,' breakpoints left'
440 call gakdel(n) !hier loeschen
441
442c..............no more breakpoints -> redo..........................
443 if(nbr.eq.0)then
444 nsbp=nsbp+1
445 if(ish.ge.3)write (ifch,*)
446 & 'no breakpoints left ... try again (',nsbp,')'
447 if(ish.ge.3)write (ifch,*)' '
448 goto11
449 endif
450
451c................make index list of breakpoints to adjust...........
452 10 continue
453 nbrtry=nbrtry+1
454 nind=0
455 xlastdiq=0.
456 do i=1,nbr
457 if(ip(i).eq.0.or.ip(i+1).eq.0)then
458 nind=nind+1
459 ind(nind)=i
460 endif
461 enddo
462 if(nbrtry.gt.1000000)then
463 if(ish.ge.1)then
464 write(*,*)'Gakfra : string can not be fragmented, redo event !'
465 write(*,*)nk1,nk2,nbr,nind,pb
466 endif
467 iret=1
468 goto 9999
469 endif
470
471c.....no more breakpoint to adjust...............................
472 if(nind.eq.0) goto 100
473
474c................................................................
475 if(ish.ge.5)then
476 write(ifch,*) 'breakpoints:'
477 write(ifch,'(i3,i5)') 0,-iflb(0)
478 do i=1,nbr
479 write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
480 & i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
481 & ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
482 enddo
483 write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
484 endif
485
486c.....choose breakpoint, calculate masses, lambda..................
487 r=rangen()
488 nn=1+int(r*real(nind))
489c nn=1+(nind-1)*int(r+0.5)
490 n=ind(nn)
491 if(ish.ge.4)write(ifch,*) 'choose breakpoint',n
492 nl=n-1
493 nr=n+1
494 do while (nl.gt.0.and.ip(nl).eq.0)
495 nl=nl-1
496 enddo
497 do while (nr.lt.nbr+1.and.ip(nr+1).eq.0)
498 nr=nr+1
499 enddo
500 if(ish.ge.6)write (ifch,*) 'nl,n,nr:',nl,n,nr
501c print *,'------------------------------------',1
502 call gaksco(n-1,n,n+1,ijb(1,n),ijb(2,n),x1,x2,y1,y2)
503 if(x2.le.x1.or.y2.le.y1)goto 9
504cc x=(xyb(1,n)-x1)/(x2-x1)
505cc y=(xyb(2,n)-y1)/(y2-y1)
506 x=xyb(1,n)
507 y=xyb(2,n)
508 aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
509 amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
510 ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
511
512c.....determine id of both particles..............................
513 aml=sqrt(max(0.,aml2))
514 idl=idsp( -iflb(n-1) , iflb(n) )
515 amr=sqrt(max(0.,amr2))
516 idr=idsp( -iflb(n) , iflb(n+1) )
517
518c.....if mass to small (because of spacelike pt) reject..........
519c amin=0.0
520c if (aml2.le.amin.or.amr2.le.amin) goto 9
521
522c.....if no left or right ptl id -> reject.........
523 if(idl.eq.0.or.idr.eq.0)then
524 if(ish.ge.5)write(ifch,*)'no left or right ptl id'
525 goto 9
526 endif
527
528 iadjl=0
529 iadjr=0
530 idrl=0
531 idrr=0
532
533c.....if no adjusted mass on left side, check for resonance...........
534 if(ip(n) .eq.0) then
535 aml=sqrt(max(0.,aml2+0.*min(0.,amr2))) !!!????
536 amlxx=aml
537 call idress(idl,aml,idrl,iadjl)
538 r=rangen()
539 if(idrl.eq.110.and.r.lt.0.5)then
540 if (r.gt.eta+etap) goto 9 !rare numerical errors
541 idl=220
542 aml=.549
543 if(r.lt.etap)aml=.958
544 amlxx=aml
545 call idress(idl,aml,idrl,iadjl)
546 iadjl=1
547 endif
548 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
549 & ' l: ',idl,amlxx,aml,idrl,iadjl
550 else
551 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
552 & ' l: ',idl,aml,aml,ip(n),'ok'
553 endif
554
555c.....if no adjusted mass on right side, check for resonance...........
556 if(ip(n+1).eq.0) then
557 amr=sqrt(max(0.,amr2+0.*min(0.,aml2))) !!!????
558 amrxx=amr
559 call idress(idr,amr,idrr,iadjr)
560 r=rangen()
561 if(idrr.eq.110.and.r.lt.0.5)then
562 if (r.gt.eta+etap) goto 9 !rare numerical errors
563 idr=220
564 amr=.549
565 if(r.lt.etap)amr=.958
566 amrxx=amr
567 call idress(idr,amr,idrr,iadjr)
568 iadjr=1
569 endif
570 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
571 & ' r: ',idr,amrxx,amr,idrr,iadjr
572 else
573 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
574 & ' r: ',idr,amr,amr,ip(n+1),'ok'
575 endif
576
577c.....mass adjustments.........................................
578 iok=0
579 if(ip(n+1).ne.0)then !.........adjusted mass on right side
580 if(idrl.eq.0)then
581 call gaksbp(n-1,n,iok)
582 if(iok.eq.1)goto 9
583 goto 10
584 endif
585 if(iadjl.eq.1)then
586 if(ish.ge.5)write(ifch,*)'mass adjustment 1'
587 n2=n+1
588 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
589 endif
590 if(iok.eq.1)goto 9
591 ip(n)=idrl
592 elseif(ip(n).ne.0)then !.........adjusted mass on left side
593 if(idrr.eq.0)then
594 call gaksbp(n,n+1,iok)
595 if(iok.eq.1)goto 9
596 goto 10
597 endif
598 if(iadjr.eq.1)then
599 if(ish.ge.5)write(ifch,*)'mass adjustment 2'
600 n2=n+1
601 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
602 endif
603 if(iok.eq.1)goto 9
604 ip(n+1)=idrr
605 else !.........adjusted mass neither left nor right
606 if(idrr.eq.0.and.idrl.eq.0)then
607 call gaksbp(n,n+1,iok)
608 if(iok.eq.1)goto 9
609 call gaksbp(n-1,n,iok)
610 if(iok.eq.1)goto 9
611 goto 10
612 elseif(idrl.ne.0.and.idrr.eq.0)then
613 if(iadjl.eq.1) then
614 if(ish.ge.5)write(ifch,*)'mass adjustment 3'
615 call gakanp(n-1,n,nr,aml**2,0.,ala2,iok)
616 endif
617 elseif(idrl.eq.0.and.idrr.ne.0)then
618 if(iadjr.eq.1) then
619 if(ish.ge.5)write(ifch,*)'mass adjustment 4'
620 n2=n+1
621 call gakanp(nl,n,n2,0.,amr**2,ala2,iok)
622 endif
623 else !if(idrl.ne.0.and.idrr.ne.0)then
624 if(iadjl.eq.1.or.iadjr.eq.1) then
625 if(ish.ge.5)write(ifch,*)'mass adjustment 5'
626 n2=n+1
627 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
628 endif
629 endif
630 if(iok.eq.1)goto 9
631 ip(n)=idrl
632 ip(n+1)=idrr
633 endif
634 if(ish.ge.4)then
635 write(ifch,*) 'left/right particles:'
636 & ,ip(n),aml,' ',ip(n+1),amr
637 endif
638 goto 10
639
640c................................................................
641c end of string decay
642c................................................................
643
644c.....final list...............................................
645 100 if(ish.ge.4)then
646 write(ifch,*) ' ************ OK **************'
647 write(ifch,*) 'final breakpoints:'
648 write(ifch,'(i3,i5)') 0,-iflb(0)
649 do i=1,nbr
650 write(ifch,'(i3,2i5,1x,4e12.4,2x,2i3,2f6.3)')
651 & i,iflb(i),-iflb(i),(ptb(j,i),j=1,4)
652 & ,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
653 enddo
654 write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
655 endif
656
657c.....write particles in pptl-list................................
658 if(ish.ge.3)then
659 write(ifch,'(a)')' ---------produced particles---------'
660 endif
661c to avoid leading strange baryon for diffractive string, rotate string
662
663 zz=0.
664 if(kdrm.eq.1.and.krm.ne.0.and.nbr.gt.0)then
665 iq1=abs(iflb(1))
666 iq2=abs(iflb(2))
667 if(iq1.gt.3.or.iq2.gt.2)then
668 am=0.
669 elseif(iq1.eq.3)then
670 am=qmass(iq1)
671 elseif(iso.gt.0)then !for u or d, use proper mass if isospin > 0 (proton)
672 am=qmass(iq1)
673 elseif(iso.lt.0)then !us d mass for u and u mass for u if isospin < 0
674 am=qmass(3-iq1)
675 else
676 am=0.5*(qmass(1)+qmass(2))
677 endif
678 zz=zrminc*pf(nob,nob)**2.
679 zz=zz*(1.+zrmmax*log(float(matarg+maproj-1)))
680c if(itp/10.eq.4)then
681c zz=zz*(1.+zrmmax*log(max(1.,float(matarg)/float(maproj))))
682c elseif(itp/10.eq.5)then
683c zz=zz*(1.+zrmmax*log(max(1.,float(maproj)/float(matarg))))
684c endif
685 zz=am*zz
686c print *,iso,itp,pf(nob,nob),am,zz,iq1,iq2,exp(-min(100.,zz))
687 endif
688c if(kdrm.eq.1)print *,'reminv',itp,am,pf(nob,nob),zz,iq1,iq2
689c & ,reminv+exp(-min(100.,zz))
690 if(rangen().gt.reminv+exp(-min(100.,zz)))then
691 leadstr=.true.
692 else
693 leadstr=.false.
694 endif
695
696 nptlini=nptl
697 if(nptlini+nbr+1.gt.mxptl)
698 & call utstop('gakfra (end): mxptl too small ...&')
699 do i=1,nbr+1
700 nptl=nptl+1
701 if(i.lt.nbr+1)then !particle = left side of breakpoints
702 call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
703 taubrr=pst(4,nob+7)+x*pst(4,nob+8)+y*pst(4,nob+9)
704 x=xyb(1,i)
705 y=xyb(2,i)
706 aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
707 amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
708
709 ala2=co(9)+x*co(10)+y*co(11)+x*y*co(12)
710 aml=sign(sqrt(abs(aml2)),aml2)
711 amr=sign(sqrt(abs(amr2)),amr2)
712 if(aml.le.0.d0.or.amr.le.0.d0)then
713 if(ish.ge.4)write(ifch,*)
714 & 'Negative mass, fragmentation not OK -> continue ...'
715 n=i
716 nptl=nptlini
717 goto 9
718 endif
719 qsqptl(nptl)=ala2
720 pptl(5,nptl)=aml
721 do j=1,4
722 pptl(j,nptl)=pst(j,nob+1)-
723 & x*pst(j,nob+2)+y*pst(j,nob+3)
724c pptr(j,nptl)=ptb(j,i)
725 enddo
726 else !last particle = right side of last breakpoint
727 pptl(5,nptl)=amr
728 qsqptl(nptl)=0.
729 do j=1,4
730 pptl(j,nptl)=pst(j,nob+4)+
731 & x*pst(j,nob+5)-y*pst(j,nob+6)
732c pptr(j,nptl)=ptb(j,i) !should be zero
733 enddo
734 taubrr=0.
735 endif
736 idptl(nptl)=ip(i)
737
738 enddo
739
740 if(leadstr)then
741 iptmp=idptl(nptlini+2)
742 idptl(nptlini+2)=idptl(nptlini+1)
743 idptl(nptlini+1)=iptmp
744 pptltmp=pptl(5,nptlini+2)
745 pptl(5,nptlini+2)=pptl(5,nptlini+1)
746 pptl(5,nptlini+1)=pptltmp
747 pptl(4,nptlini+1)=pptl(5,nptlini+1)*pptl(5,nptlini+1)
748 pptl(4,nptlini+2)=pptl(5,nptlini+2)*pptl(5,nptlini+2)
749 do k=1,3
750 pptl(4,nptlini+1)=pptl(4,nptlini+1)
751 & +pptl(k,nptlini+1)*pptl(k,nptlini+1)
752 pptl(4,nptlini+2)=pptl(4,nptlini+2)
753 & +pptl(k,nptlini+2)*pptl(k,nptlini+2)
754 enddo
755 pptl(4,nptlini+1)=sqrt(pptl(4,nptlini+1))
756 pptl(4,nptlini+2)=sqrt(pptl(4,nptlini+2))
757 endif
758 nptl=nptlini
759 do i=1,nbr+1
760 nptl=nptl+1
761
762 call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
763 $ ,pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl))
764c call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
765c $ ,pptr(1,nptl),pptr(2,nptl),pptr(3,nptl),pptr(4,nptl))
766
767c........Origin..................................................
768 istptl(nptl)=0
769 iorptl(nptl)=n8ptl
770 jorptl(nptl)=0
771 ifrptl(1,nptl)=0
772 ifrptl(2,nptl)=0
773
774 r=rangen()
775 tauran=-taurea*alog(r)*pptl(4,nptl)/pptl(5,nptl)
776c if(jpsi.ge.1)tauran=0.
777c tauran=max(taubrl,taubrr) !take formation time from string-theory
778 do j=1,4
779 xorptl(j,nptl)=xorptl(j,nk1)+pptl(j,nptl)
780 & /pptl(4,nptl)*tauran
781 enddo
782 tivptl(1,nptl)=xorptl(4,nptl)
783 call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
784 tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
785 ityptl(nptl)=ityptl(nk1)
786 itsptl(nptl)=itsptl(nk1)
787 rinptl(nptl)=-9999
788 if(ish.ge.3)call blist('&',nptl,nptl)
789ctp060829 taubrl=taubrr
790 enddo
791 iorptl(n8ptl)=nk1
792 jorptl(n8ptl)=nk2
793 ityptl(n8ptl)=ityptl(nk1)
794 do i=nk1,nk2
795 istptl(i)=21
796 ifrptl(1,i)=n8ptl
797 ifrptl(2,i)=0
798 enddo
799 istptl(n8ptl)=29 !code for fragmented string
800 ifrptl(1,n8ptl)=n8ptl+1
801 ifrptl(2,n8ptl)=nptl
802 if(ish.ge.5)then
803 write(ifch,*)'string momentum sum:'
804 do kk=1,5
805 pptl(kk,nptl+1)=0
806 do ii=n8ptl+1,nptl
807 pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
808 enddo
809 enddo
810 call alist2('&',n8ptl,n8ptl,nptl+1,nptl+1)
811 endif
812
813
814
815c.....another string?.........................................
816 98 nk1=nk2+1
817 goto 2 !next string
818
819c.....end of fragmentation.....................................
820 9999 continue
821 call utprix('gakfra',ish,ishini,4)
822 return
823 end
824
825c---------------------------------------------------------------------
826 subroutine gaksbp(ibr1,ibr2,iok)
827c---------------------------------------------------------------------
828 ! search break points
829 !-----------------------------------------
830 ! nbr ... number of break points
831 !
832 !
833 !
834 !--------------------------------------------------------------
835 include 'epos.inc'
836
837 parameter (mxnbr=500,mxnba=5000)
838 common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
839 $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
840 & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
841 common/pb/pb /cnsbp/nsbp /cn8ptl/n8ptl
842 double precision ax,ay,az,ae,am,bx,by,bz,be
843 dimension co(12)
844 logical weiter
845 common/czz/kky,krm,zzod,zzop,kdrm,zzoq,zzos,zzrm,zzipp,zzipt
846 &,zzid,zzip,zzzzr,itp
847 double precision psg
848 common/zpsg/psg(5)
849 dimension pleft(5),pright(5)
850
851 pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
852 & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
853 mmod(i)=mod(mod(i,nob)+nob,nob)
854
855 call utpri('gaksbp',ish,ishini,5)
856
857
858 ib1=ibr1
859 ib2=ibr2
860 iside=1
861 if(rangen().lt.0.5)iside=2
862 if(ish.ge.6)write(ifch,*)'iside:',iside
863 Amxx=80./pb
864 if(ish.ge.4)write(ifch,*)
865 &'search brk between ib1=',ib1,' and ib2=',ib2
866 ntry=0
867 nbrold=nbr
868 24 Amin=0.
869 25 Amax=Amxx
870 26 ntry=ntry+1
871 A0=-log(exp(-pb*Amin)+rangen()*(exp(-pb*Amax)-exp(-pb*Amin)))/pb
872 if(ish.ge.6)write(ifch,*)'pb, Amin, Amax, A0:',pb, Amin, Amax, A0
873 A=A0
874 Amaxn=0.
875 if(iside.eq.2)goto 51
876c.....................................................................
877 if(ib1.eq.0)then !startwert der j-schleife
878 jj=nob-1
879 else
880 jj=ijb(2,ib1)
881 endif
882 do while (jj.gt.0)
883 if(jj.eq.ijb(2,ib1) )then !linker brk im Streifen jj?
884 y1=xyb(2,ib1)
885 else !nein
886 y1=0.
887 endif
888 if(jj.eq.ijb(2,ib2))then !rechter brk im Streifen jj?
889 y2=xyb(2,ib2)
890 else !nein
891 y2=1.
892 endif
893 if(y1.eq.y2) goto 9999
894 if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
895 if(ish.ge.6)write(ifch,*)'y1,y2,A',y1,y2,A
896 !startwert der i-schleife
897 if(ib1.eq.0)then !ohne linken brkpt
898 ii=mmod(jj+1)
899 if(jj.lt.nob/2)ii=mmod(nob-jj+1)
900 if(ib2.le.nbr.and.mmod(ii+nob/2)
901 & .gt.mmod(ijb(1,ib2)+nob/2))then
902 if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
903 goto12
904 endif
905 else
906 ii=ijb(1,ib1)
907 if(jj.lt.nob/2 .and.
908 $ mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
909 $ )ii=mmod(nob-jj+1)
910 endif
911 weiter=.true.
912 aa=0.
913 if(ii.eq.jj) then
914 if(ish.ge.6)write(ifch,*) 'Rand erreicht'
915 goto 15
916 endif
917 do while(weiter)
918 if(ii.eq.ijb(1,ib1))then !linker brk im Feld (ii,jj)
919 x2=xyb(1,ib1)
920 else
921 x2=1.
922 endif
923
924 if(ii.eq.ijb(1,ib2))then !rechter brk im Feld (ii,jj)
925 x1=xyb(1,ib2)
926 else
927 x1=0.
928 endif
929 if(x1.eq.x2) goto 9999
930 if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
931 f=1.0
932 if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
933 if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
934 $ ,ipst(ii).ne.ipst(jj),x1,x2,aa
935 & ,pf(ii,jj)
936 if(ii.eq.ijb(1,ib2))then
937 weiter=.false.
938 else
939 ii=mmod(ii+1)
940 if(ii.eq.jj.or.mmod(ii+jj).eq.0)weiter=.false.
941 endif
942 enddo
943 Amaxn=Amaxn+aa
944 if(aa.gt.A)goto 10
945 A=A-aa
946 if(jj.eq.ijb(2,ib2)) then
947 if(ish.ge.6)write(ifch,*) 'brk erreicht'
948 goto 15
949 endif
950 12 jj=mmod(jj-1)
951 enddo
952 goto 15
953
954 10 continue
955 yb=A/aa*(y2-y1)+y1
956 if(ish.ge.6)write(ifch,*)'jj,yb ok:',jj,yb
957 r=rangen()
958 ra=aa*r
959 if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
960 !nochmal die letzte ii-Schleife
961 if(ib1.eq.0)then !ohne linken brkpt
962 ii=mmod(jj+1)
963 if(jj.lt.nob/2)ii=mmod(nob-jj+1)
964 else
965 ii=ijb(1,ib1)
966 if(jj.lt.nob/2 .and.
967 $ mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
968 $ )ii=mmod(nob-jj+1)
969 endif
970 do while(ra.gt.0.)
971 if(ii.eq.ijb(1,ib1))then !linker brk im Feld (ii,jj)
972 x2=xyb(1,ib1)
973 else
974 x2=1.
975 endif
976
977 if(ii.eq.ijb(1,ib2))then !rechter brk im Feld (ii,jj)
978 x1=xyb(1,ib2)
979 else
980 x1=0.
981 endif
982 f=1.0
983 ab=0.
984 if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
985 if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
986 $ ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
987 if(ab.gt.ra)then
988 xb=ra/ab*(x2-x1)+x1
989 else
990 ii=mmod(ii+1)
991 endif
992 ra=ra-ab
993 enddo
994 if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
995 & ,' at position ',xb,yb
996 goto 95
997
998c......................................................................
999 !von rechts
1000 51 if(ib2.eq.nbr+1)then !startwert der i-schleife
1001 ii=nob/2-1
1002 else
1003 ii=ijb(1,ib2)
1004 endif
1005 do while (ii.ne.nob/2)
1006 if(ii.eq.ijb(1,ib1))then !linker brk im Streifen (ii)
1007 x2=xyb(1,ib1)
1008 else
1009 x2=1.
1010 endif
1011 if(ii.eq.ijb(1,ib2))then !rechter brk im Streifen (ii)
1012 x1=xyb(1,ib2)
1013 else
1014 x1=0.
1015 endif
1016 if(x1.eq.x2) goto 9999
1017 if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
1018 if(ish.ge.6)write(ifch,*)'x1,x2 A',x1,x2,A
1019
1020 if(ib2.eq.nbr+1)then
1021 jj=mmod(ii+1)
1022 if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1023 if(ib1.ge.1.and.jj.gt.ijb(2,ib1))then
1024 if(ish.ge.6)write(ifch,*) 'very special case',ii,jj
1025 goto 13
1026 endif
1027 else
1028 jj=ijb(2,ib2)
1029 if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1030 endif
1031 weiter=.true.
1032 aa=0.
1033 if(ii.eq.jj) then
1034 if(ish.ge.6)write(ifch,*) 'Rand erreicht'
1035 goto 15
1036 endif
1037 do while(weiter)
1038 if(jj.eq.ijb(2,ib1))then !linker brk im Feld (ii,jj)
1039 y1=xyb(2,ib1)
1040 else
1041 y1=0.
1042 endif
1043 if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
1044 y2=xyb(2,ib2)
1045 else
1046 y2=1.
1047 endif
1048 if(y1.eq.y2) goto 9999
1049 if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
1050 f=1.0
1051 if(ipst(ii).ne.ipst(jj))aa=aa+2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1052 if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,aa:',ii,jj
1053 $ ,ipst(ii).ne.ipst(jj),x1,x2,aa
1054 & ,pf(ii,jj)
1055 if(jj.eq.ijb(2,ib1))then
1056 weiter=.false.
1057 else
1058 jj=mmod(jj+1)
1059 if(jj.eq.ii.or.mmod(ii+jj).eq.0)weiter=.false.
1060 endif
1061 enddo
1062 Amaxn=Amaxn+aa
1063 if(aa.gt.A)goto 14
1064 A=A-aa
1065 if(ii.eq.ijb(1,ib1)) then
1066 if(ish.ge.6)write(ifch,*) 'brk erreicht'
1067 goto 15
1068 endif
1069 13 ii=mmod(ii-1)
1070 enddo
1071 goto 15
1072
1073
1074
1075 14 continue
1076 xb=A/aa*(x2-x1)+x1
1077 if(ish.ge.6)write(ifch,*)'ii,xb ok:',ii,xb
1078 r=rangen()
1079 ra=aa*r
1080 if(ish.ge.6)write(ifch,*)'r,ra,aa',r,ra,aa
1081 !nochmal die letzte jj-Schleife
1082 if(ib2.eq.nbr+1)then
1083 jj=mmod(ii+1)
1084 if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1085 else
1086 jj=ijb(2,ib2)
1087 if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1088 endif
1089
1090 do while(ra.gt.0.)
1091 if(jj.eq.ijb(2,ib1))then !linker brk im Feld (ii,jj)
1092 y1=xyb(2,ib1)
1093 else
1094 y1=0.
1095 endif
1096
1097 if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
1098 y2=xyb(2,ib2)
1099 else
1100 y2=1.
1101 endif
1102 f=1.0
1103 ab=0.
1104 if(ipst(ii).ne.ipst(jj)) ab=2.*(x2-x1)*(y2-y1)*f*pf(ii,jj)
1105 if(ish.ge.6)write(ifch,*)'ii,jj,x1,x2,ab,ra:',ii,jj
1106 $ ,ipst(ii).ne.ipst(jj),x1,x2,ab,ra
1107 if(ab.gt.ra)then
1108 yb=ra/ab*(y2-y1)+y1
1109 else
1110 jj=mmod(jj+1)
1111 endif
1112 ra=ra-ab
1113 enddo
1114 if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
1115 & ,' at position ',xb,yb
1116
1117 95 continue
1118
1119c.....breakpoint accepted......................
1120 nbr=nbr+1
1121 if(nbr.gt.mxnbr-2) stop 'gaksbp: increase parameter mxnbr'
1122 do i=nbr+1,ib1+1,-1
1123 do j=1,2
1124 ijb(j,i)=ijb(j,i-1)
1125 xyb(j,i)=xyb(j,i-1)
1126 enddo
1127 do k=1,4
1128 ptb(k,i)=ptb(k,i-1)
1129 enddo
1130 iflb(i)=iflb(i-1)
1131 ip(i)=ip(i-1)
1132 enddo
1133 ip(ib1+1)=0
1134 ip(ib1+2)=0
1135 ijb(1,ib1+1)=ii
1136 ijb(2,ib1+1)=jj
1137 xyb(1,ib1+1)=xb
1138 xyb(2,ib1+1)=yb
1139
1140 if(kdrm.eq.1)then
1141 f=2.
1142 if(ib1.eq.0)f=1.
1143 if(abs(iflb(ib1)).gt.3)then
1144 iq1=mod(abs(iflb(ib1))/1000,10)
1145 iq2=mod(abs(iflb(ib1))/100,10)
1146 amlast=qmass(0)+f*(qmass(iq1)+qmass(iq2))
1147 else
1148 amlast=f*qmass(abs(iflb(ib1)))
1149 endif
1150 if(itp.ge.50)then
1151 ampt=amtarg
1152 elseif(itp.ge.40)then
1153 ampt=amproj
1154 endif
1155 xlastdiq=0.
1156 do i=0,ib1-1
1157 f=2.
1158 if(i.eq.0)f=1.
1159 if(abs(iflb(i)).gt.3)then
1160 iq1=mod(abs(iflb(i))/1000,10)
1161 iq2=mod(abs(iflb(i))/100,10)
1162 xlastdiq=xlastdiq+qmass(0)+f*(qmass(iq1)+qmass(iq2))
1163 else
1164 xlastdiq=xlastdiq+f*qmass(abs(iflb(i)))
1165 endif
1166 enddo
1167 if(pf(nob,nob).le.5.)amlast=1.e-10 !if very low energy, strangeness should not be totally killed
1168c if(pf(nob,nob).le.5.)then
1169c xlastdiq=-(ampt+2.*qmass(3)) !if very low energy, strangeness should not be totally killed
1170c amlast=strcut**2
1171c endif
1172 if(ish.ge.5)write(ifch,*) 'diff cut:',pf(nob,nob),xlastdiq,amlast
1173 & ,float(ib1+1)-ampt-xlastdiq+diqcut*sqrt(pf(nob,nob))
1174 & ,strcut/amlast/ampt*pf(nob,nob)
1175 endif
1176
1177c ..........diquark...............................................
1178
1179
1180 if(krm.eq.1)then
1181 if(kdrm.eq.0)then
1182 pdiqu=pdiquar*zzod
1183 else
1184 pdiqu=pdiquad
1185 endif
1186 elseif(kky.eq.1)then
1187 pdiqu=pdiquak*zzid
1188 else
1189 pdiqu=pdiqua!*zzid
1190 endif
1191
1192 if(kdrm.eq.1)pdiqu=pdiqu*(1.-exp(-min(100.,
1193 & strcut/amlast/qmass(0)/ampt
1194 & *max(0.,float(ib1+1)-ampt-xlastdiq-2.*qmass(0)+diqcut*
1195 & sqrt(pf(nob,nob)))*pf(nob,nob))))
1196c if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.diqcut)pdiqu=0.
1197
1198 if(rangen().lt.pdiqu.and.iabs(iflb(ib1)).lt.6
1199 & .and.iabs(iflb(ib2)).lt.6)then
1200 jqu=2
1201 else
1202 jqu=1
1203 endif
1204c ..........flavor...............................................
1205 77 continue
1206 if(krm.eq.1)then
1207 if(kdrm.eq.0)then
1208 pud1= ( 1 - min(0.33,zzos*(1-2*pudr)) ) /2.
1209c pud2= pud1
1210 else
1211 pud1= pudd
1212c pud2= pud1
1213 endif
1214 elseif(kky.eq.1)then
1215 pud1=pudk
1216c pud2=pud1
1217 else
1218 pud1= pud
1219c pud2= pud1
1220 endif
1221 if(kdrm.eq.1)pud1=0.5*(1.-(1.-2.*pud1)
1222 & *(1.-exp(-min(100.,strcut/amlast/qmass(3)/ampt
1223 & *max(0.,float(ib1+1)-ampt-xlastdiq-2.*qmass(3)+diqcut
1224 & *sqrt(pf(nob,nob)))*pf(nob,nob)))))
1225c if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.strcut)pud1=0.5
1226
1227 pud2=(1.-(1.-2.*pud1)*puds)*0.5 !from e+e- data
1228
1229 if(iabs(iflb(ib1)).gt.1000)then
1230 ifl1=int(rangen()**0.625/pud1)+1
1231 else
1232 ifl1=int(rangen()**1.15/pud1)+1 !assymmetric u and d relevant because of assymmetry Kch / K0 in e+e-
1233 endif
1234
1235c amk0=2.*qmass(ifl1) !mass for mt distribution
1236 if(ifl1.eq.4)ifl1=3
1237 ifl2=0
1238 if(jqu.eq.2)then !diquark ------
1239 ifl2=int(rangen()**0.625/pud2)+1
1240c amk0=amk0+2.*qmass(ifl2)+qmass(0) !mass for mt distribution with bounding energy for diquark (qmass(0))
1241 if(ifl2.eq.4)ifl2=3
1242 ifl=-min(ifl1,ifl2)*1000-max(ifl1,ifl2)*100
1243 else !quark ------
1244 ifl=ifl1
1245 endif
1246 if(iflb(ib1).lt.0.and.iflb(ib1).ge.-6)ifl=-ifl
1247 if(iflb(ib1).gt.1000)ifl=-ifl
1248 iflb(ib1+1)=ifl
1249 if(ish.ge.5)write(ifch,*) 'flavor,pud1/2:',ifl,pud1,pud2,pdiqu
1250c..............................pt.......................................
1251 if(krm.ne.1.and.kdrm.ne.1)then
1252 icub=ib1+1
1253 !---------------------------------------------------------------------------
1254 ! ib1+1 is the current break point index
1255 ! (between 1 and nbr)
1256 ! ijb(1,ib1) and ijb(2,ib1) are band indices
1257 ! each index from 0 to nob-1 (nob= number of bands)
1258 ! 0 is left, then come the gluons, then antiquark, then agin gluons
1259 ! like q - g - g - g - ~q - g - g - g
1260 !--------------------------------------------------------------------------
1261 call gaksco(icub-1,icub,icub+1
1262 & ,ijb(1,icub),ijb(2,icub),x1,x2,y1,y2)
1263 x=xyb(1,icub)
1264 y=xyb(2,icub)
1265 aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
1266 amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
1267 aml=sign(sqrt(abs(aml2)),aml2)
1268 amr=sign(sqrt(abs(amr2)),amr2)
1269 !------segment left of current breakpoint icub -----
1270 pleft(5)=aml
1271 do j=1,4
1272 pleft(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
1273 enddo
1274 call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1275 $ ,pleft(1),pleft(2),pleft(3),pleft(4))
1276 !------segment right of current breakpoint icub-----
1277 pright(5)=amr
1278 do j=1,4
1279 pright(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
1280 enddo
1281 call utlob4(-1,psg(1),psg(2),psg(3),psg(4),psg(5)
1282 $ ,pright(1),pright(2),pright(3),pright(4))
1283 !-------------------------
1284 amt=pleft(5)**2+pleft(1)**2+pleft(2)**2
1285 if(amt.gt.0..and.pleft(4)+abs(pleft(3)).gt.0.d0)then
1286 amt=sqrt(amt)
1287 yleft=sign(1.,pleft(3))*alog((pleft(4)+abs(pleft(3)))/amt)
1288 else
1289 yleft=0. !
1290 endif
1291 amt=pright(5)**2+pright(1)**2+pright(2)**2
1292 if(amt.gt.0..and.pright(4)+abs(pright(3)).gt.0.d0)then
1293 amt=sqrt(amt)
1294 yright=sign(1.,pright(3))*alog((pright(4)+abs(pright(3)))/amt)
1295 else
1296 yright=0. !
1297 endif
1298 ybrk=(yleft+yright)/2.
1299 yhax=2 !0.5*yhaha
1300 zzipx=1
1301 if(ybrk.gt.yhax)then
1302 zzipx=zzipx+(zzipp-1)
1303 elseif(ybrk.gt.-yhax)then
1304 zzipx=zzipx+(zzipp-1)*(ybrk+yhax)/(2*yhax)
1305 endif
1306 if(ybrk.lt.-yhax)then
1307 zzipx=zzipx+(zzipt-1)
1308 elseif(ybrk.lt.yhax)then
1309 zzipx=zzipx+(zzipt-1)*(yhax-ybrk)/(2*yhax)
1310 endif
1311 endif
1312 delptfra=0
1313 if(ifl1.eq.3.and.ifl2.eq.0)delptfra=delptfra+ptfrasr
1314c if(ifl1.eq.3.or.ifl2.eq.3)delptfra=delptfra+ptfrasr
1315 fnsbp=1
1316 if(nsbp.gt.9)fnsbp=0
1317
1318c for a better pt distribution, we generate mt=sqrt(pt2+m02) (with m0=amk0) instead of pt (do not give good results !)
1319 if(krm.eq.1)then
1320 if(kdrm.eq.0)then
1321c ptadd=ranpt()*(ptfrair+delptfra)*(zzop-1)*fnsbp !*jqu
1322c ptadd=ranpt()*ptfrair*(zzop-1)*fnsbp !*jqu
1323 ptadd=(ptfrair+delptfra)*(zzop-1)*fnsbp !*jqu
1324 if(jqu.eq.1)then
1325c pt=ranpt()*(ptfrair +delptfra ) + ptadd
1326c pt=ranpt()*ptfrair
1327 pt=ranpt()*(ptfrair + ptadd )
1328 else
1329 pt=ranpt()*(ptfrair + ptadd )
1330c pt=ranpt()*ptfraqq
1331c pt=ranpt()*(ptfraqq + ptadd )
1332 endif
1333 else
1334c ptadd=0.
1335 pt=ranpt()*ptfradr
1336 if(jqu.eq.1)then
1337c pt=ranpt()*(ptfradr +delptfra )
1338 pt=ranpt()*ptfradr
1339 else
1340 pt=ranpt()*ptfraqq
1341c pt=ranpt()*(ptfradr )
1342 endif
1343 endif
1344 elseif(kky.eq.1)then
1345c ptadd=ranptk()*(ptfrak+delptfra)*(zzipx-1)*fnsbp !*jqu
1346 pttra=(ptfrak+delptfra)*(zzipx-1)*fnsbp !/jqu
1347 if(jqu.eq.1)then
1348 !pt=ranptk()*(ptfrak+delptfra + pttra )
1349 !pt=ranptk()*(ptfrak + pttra)
1350 pt=ranptk()*ptfrak + ranptk()*pttra
1351 else
1352 !pt=ranptk()*(ptfrakd + pttra)
1353 pt=ranptk()*ptfrakd + ranptk()*pttra
1354 endif
1355 else
1356c ptadd=ranpt()*(ptfra+delptfra)*(zzip-1)*fnsbp !*jqu
1357c pttra=ranpt()*ptfra*(zzip-1)*fnsbp !*jqu
1358 pttra=(ptfra+delptfra)*(zzip-1)*fnsbp !*jqu
1359 if(jqu.eq.1)then
1360c pt=ranpt()*(ptfra+delptfra ) + pttra
1361c pt=ranpt()*ptfra
1362 pt=ranpt()*(ptfra + pttra )
1363 else
1364c pt=ranpt()*(ptfra+delptfra) + pttra
1365c pt=ranpt()*ptfrakd
1366 pt=ranpt()*(ptfra + pttra )
1367 endif
1368 endif
1369c pt=sqrt(pt*pt+2.*pt*amk0)+ptadd !sample mt-m0 instead of pt ... doesn't work at RHIC
1370 beta=2.*pi*rangen()
1371 60 continue
1372 if(ish.ge.5)then
1373 write(ifch,*) 'pt:',pt
1374 endif
1375 ptb(1,ib1+1)=pt*cos(beta)
1376 ptb(2,ib1+1)=pt*sin(beta)
1377 ptb(3,ib1+1)=0.
1378 ptb(4,ib1+1)=0.
1379 if(ish.ge.8)then
1380 write(ifch,*) 'the bands'
1381 write(ifch,'(4g12.6)') (pst(i,ii),i=1,4)
1382 write(ifch,'(4g12.6)') (pst(i,jj),i=1,4)
1383 write(ifch,*) 'pt before rotation'
1384 write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1385 endif
1386 ax=dble(pst(1,ii))+dble(pst(1,jj))
1387 ay=dble(pst(2,ii))+dble(pst(2,jj))
1388 az=dble(pst(3,ii))+dble(pst(3,jj))
1389 ae=dble(pst(4,ii))+dble(pst(4,jj))
1390 am=sqrt(max(1d-10,(ae+az)*(ae-az)-ax**2-ay**2)) !?????????
1391 if(ish.ge.8)then
1392 write(ifch,*) 'boost vector ( region ) '
1393 write(ifch,'(5g12.6)') ax,ay,az,ae,am,pf(ii,jj)
1394 endif
1395 bx=pst(1,ii)
1396 by=pst(2,ii)
1397 bz=pst(3,ii)
1398 be=pst(4,ii)
1399 call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1400 if(ish.ge.8) then
1401 write (ifch,*) 'boost of b'
1402 write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1403 endif
1404 call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1405 if(ish.ge.8) then
1406 write (ifch,*) 'rot of pt'
1407 write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1408 endif
1409
1410 call utlob4(-1,ax,ay,az,ae,am
1411 & ,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1),ptb(4,ib1+1))
1412
1413 if(ish.ge.8) then
1414 write (ifch,*) 'backboost of pt'
1415 write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
1416 endif
1417 if(az.eq.0..and.ay.eq.0.)az=1e-8
1418 if(ish.ge.8)write(ifch,*)'rot vector:',ax,ay,az,ae,am
1419
1420c.....call utrota(-1,ax,ay,az,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1421
1422 if(ish.ge.6)then
1423 write(ifch,*) 'pt'
1424 write(ifch,'(4g12.6)') (ptb(i,ib1+1),i=1,4)
1425 write (ifch,*) ijb(1,ib1+1),ijb(2,ib1+1),xyb(1,ib1+1)
1426 $ ,xyb(2,ib1+1)
1427 endif
1428
1429 if(iside.eq.1)then
1430 ib1=ib1+1
1431 ib2=ib2+1
1432 endif
1433
1434
1435c Amin=0.
1436c if(Amax.lt.Amxx) goto 15
1437c goto 25
1438
1439 15 continue
1440 if(ish.ge.6)write(ifch,*) 'Amax:',Amax,Amaxn
1441 if (nbr.eq.nbrold) then
1442 Amax=Amaxn
1443 Amin=0.
1444 if(pb*Amax.le.0..or.ntry.ge.1000)then
1445 goto 9999
1446 endif
1447 goto 26
1448 endif
1449
1450 if(ish.ge.6)then
1451 write(ifch,*) 0,iflb(0)
1452 do i=1,nbr
1453 if(i.eq.ib2) write(ifch,*) '.................'
1454 write(ifch,'(i3,2x,2(i3),2(" ",g14.7),3x,i5,4(" ",g12.6))')
1455 & i,ijb(1,i),ijb(2,i),xyb(1,i),xyb(2,i)
1456 & ,iflb(i),(ptb(j,i),j=1,4)
1457 if(i.eq.ibr1) write(ifch,*) '.................'
1458 enddo
1459 write(ifch,*) nbr+1,iflb(nbr+1)
1460 endif
1461
1462 9999 if(nbr.eq.nbrold)iok=1
1463 call utprix('gaksbp',ish,ishini,5)
1464 return
1465 end
1466
1467c----------------------------------------------------------------------
1468 function ranptcut(xcut)
1469c----------------------------------------------------------------------
1470c .........exp(-x**2)
1471 12 x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1472C 12 x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1473
1474
1475c if(xcut.gt.0.)then
1476c if(rangen().lt.x/xcut)goto 12
1477c endif
1478
1479 ranptcut=x
1480
1481 return
1482
1483c .........exp(-x)
1484c 12 xmx=50
1485c r=2.
1486c do while (r.gt.1.)
1487c 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1488c if(x.eq.0.)goto11
1489c r=rangen() / ( exp(-x)*(1+x**2) )
1490c enddo
1491c x=x/2.
1492
1493 end
1494
1495cc----------------------------------------------------------------------
1496c function ranpticut(xcut)
1497cc----------------------------------------------------------------------
1498c
1499cc .........exp(-x)
1500c
1501c 12 xmx=50
1502c r=2.
1503c do while (r.gt.1.)
1504c 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1505c if(x.eq.0.)goto11
1506c r=rangen() / ( exp(-x)*(1+x**2) )
1507c enddo
1508c x=x/2.
1509c if(rangen().gt.xcut/x)goto 12
1510c
1511c ranpticut=x
1512c
1513c end
1514
1515c----------------------------------------------------------------------
1516 function ranpt()
1517c----------------------------------------------------------------------
1518
1519c .........exp(-x)
1520 xmx=50
1521 r=2.
1522 do while (r.gt.1.)
1523 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1524 if(x.eq.0.)goto11
1525 r=rangen() / ( exp(-x)*(1+x**2) )
1526 enddo
1527 ranpte=x/2.
1528
1529c -------------
1530 ranpt=ranpte
1531 return
1532c -------------
1533c .........exp(-x**2)
1534 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1535
1536c .........exp(-sqrt(x))
1537 xmx=500
1538 r=2.
1539 do while (r.gt.1.)
1540 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1541 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1542 enddo
1543 ranpts=x/20.
1544
1545
1546 end
1547c----------------------------------------------------------------------
1548 function ranptk()
1549c----------------------------------------------------------------------
1550
1551c .........exp(-x)
1552 xmx=50
1553 r=2.
1554 do while (r.gt.1.)
1555 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1556 if(x.eq.0.)goto11
1557 r=rangen() / ( exp(-x)*(1+x**2) )
1558 enddo
1559 ranpte=x/2.
1560
1561c -------------
1562 ranptk=ranpte
1563 return
1564c -------------
1565
1566c .........exp(-x**2)
1567 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1568
1569c .........exp(-sqrt(x))
1570 xmx=500
1571 r=2.
1572 do while (r.gt.1.)
1573 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1574 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1575 enddo
1576 ranpts=x/20.
1577
1578 end
1579
1580c----------------------------------------------------------------------
1581 function ranptd()
1582c----------------------------------------------------------------------
1583
1584c .........exp(-x**2)
1585 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1586
1587
1588c -------------
1589 ranptd=ranptg
1590 return
1591c -------------
1592
1593c .........exp(-x)
1594 xmx=50
1595 r=2.
1596 do while (r.gt.1.)
1597 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1598 if(x.eq.0.)goto11
1599 r=rangen() / ( exp(-x)*(1+x**2) )
1600 enddo
1601 ranpte=x/2.
1602
1603c .........exp(-sqrt(x))
1604 xmx=500
1605 r=2.
1606 do while (r.gt.1.)
1607 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1608 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1609 enddo
1610 ranpts=x/20.
1611
1612
1613
1614
1615
1616 end
1617c----------------------------------------------------------------------
1618 subroutine gakdel(ibr)
1619c----------------------------------------------------------------------
1620 parameter (mxnbr=500,mxnba=5000)
1621 common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1622 $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1623 & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1624 dimension co(12)
1625
1626 do i=ibr,nbr+1
1627 do j=1,2
1628 ijb(j,i)=ijb(j,i+1)
1629 xyb(j,i)=xyb(j,i+1)
1630 enddo
1631 do k=1,4
1632 ptb(k,i)=ptb(k,i+1)
1633 enddo
1634 iflb(i)=iflb(i+1)
1635 ip(i)=ip(i+1)
1636 enddo
1637 ip(ibr)=0
1638 nbr=nbr-1
1639 end
1640
1641c----------------------------------------------------------------------
1642 subroutine gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1643c----------------------------------------------------------------------
1644 include 'epos.inc'
1645
1646 parameter (mxnbr=500,mxnba=5000)
1647 common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1648 $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1649 & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1650 dimension co(12)
1651 logical weiter
1652
1653 pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1654 & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1655 mmod(i)=mod(mod(i,nob)+nob,nob)
1656
1657 call utpri('gaksco',ish,ishini,8)
1658
1659 if(ish.ge.8)then
1660 write(ifch,*) 'zwischen brk:',ibr1,ibr,ibr2,'(',nob,')',nbr
1661 write(ifch,*) 'region:',ib,jb
1662 endif
1663 if(ib.eq.ijb(1,ibr1))then
1664 x2=xyb(1,ibr1)
1665 else
1666 x2=1.
1667 endif
1668 if(ib.eq.ijb(1,ibr2))then
1669 x1=xyb(1,ibr2)
1670 else
1671 x1=0.
1672 endif
1673 if(jb.eq.ijb(2,ibr1))then
1674 y1=xyb(2,ibr1)
1675 else
1676 y1=0.
1677 endif
1678 if(jb.eq.ijb(2,ibr2))then
1679 y2=xyb(2,ibr2)
1680 else
1681 y2=1.
1682 endif
1683
1684c.....left side...................................................
1685 n=nob+1
1686 if(ish.ge.8)write(ifch,*)'x1,x2',x1,x2
1687 do i=1,4
1688cc pst(i,n)=(x2-x1)*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)
1689 pst(i,n)= x2*pst(i,ib)+ptb(i,ibr)-ptb(i,ibr1)-y1*pst(i,jb)
1690 enddo
1691 if(ish.ge.8)write(ifch,*) 'add a 1 ii',1.,ib
1692 ii=mmod(ib-1)
1693 weiter=.true.
1694 if(ib.eq.ijb(1,ibr1))weiter=.false.
1695 do while(ii.ne.jb.and.weiter) !linker Rand??
1696 f1=1.
1697 if(ii.eq.ijb(1,ibr1))f1=xyb(1,ibr1)
1698 do i=1,4
1699 pst(i,n)=pst(i,n)+f1*pst(i,ii)
1700 enddo
1701 if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1702 if(ii.eq.ijb(1,ibr1))weiter=.false.
1703 ii=mmod(ii-1)
1704 enddo
1705 jj=mmod(jb+1)
1706 weiter=.not.weiter
1707 if(jb.eq.ijb(2,ibr1))weiter=.false.
1708 do while(weiter)
1709 f1=1.
1710 if(jj.eq.ijb(2,ibr1))f1=1.-xyb(2,ibr1)
1711 do i=1,4
1712 pst(i,n)=pst(i,n)+f1*pst(i,jj)
1713 enddo
1714 if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1715 if(jj.eq.ijb(2,ibr1))weiter=.false.
1716 jj=mmod(jj+1)
1717 enddo
1718 do i=1,4
1719cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1720 pst(i,n+1)= pst(i,ib)
1721cc pst(i,n+2)=(y2-y1)*pst(i,jb)
1722 pst(i,n+2)= pst(i,jb)
1723 enddo
1724 co(1)= pf(n,n)
1725 co(2)=-2.*pf(n,n+1)
1726 co(3)= 2.*pf(n,n+2)
1727 co(4)=-2.*pf(n+1,n+2)
1728 if(ish.ge.8) then
1729 do i=n,n+2
1730 write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1731 enddo
1732 endif
1733 if(ish.ge.8)write(ifch,*) 'co left:',co(1),co(2),co(3),co(4)
1734
1735c.....right side...................................................
1736 n=nob+4
1737 if(ish.ge.8)write(ifch,*)'y1,y2',y1,y2
1738 do i=1,4
1739cc pst(i,n)=(y2-y1)*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)
1740 pst(i,n)= y2*pst(i,jb)-ptb(i,ibr)+ptb(i,ibr2)-x1*pst(i,ib)
1741 enddo
1742 if(ish.ge.8)write(ifch,*) 'add a 1 jj',1.,jb
1743 ii=mmod(ib+1)
1744 weiter=.true.
1745 if(ib.eq.ijb(1,ibr2))weiter=.false.
1746 do while(ii.ne.jb.and.weiter)
1747 f1=1.
1748 if(ii.eq.ijb(1,ibr2))f1=1.-xyb(1,ibr2)
1749 do i=1,4
1750 pst(i,n)=pst(i,n)+f1*pst(i,ii)
1751 enddo
1752 if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1753 if(ii.eq.ijb(1,ibr2))weiter=.false.
1754 ii=mmod(ii+1)
1755 enddo
1756 jj=mmod(jb-1)
1757 weiter=.not.weiter
1758 if(jb.eq.ijb(2,ibr2))weiter=.false.
1759 do while(weiter)
1760 f1=1.
1761 if(jj.eq.ijb(2,ibr2))f1=xyb(2,ibr2)
1762 do i=1,4
1763 pst(i,n)=pst(i,n)+f1*pst(i,jj)
1764 enddo
1765 if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1766 if(jj.eq.ijb(2,ibr2))weiter=.false.
1767 jj=mmod(jj-1)
1768 enddo
1769 do i=1,4
1770cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1771 pst(i,n+1)= pst(i,ib)
1772cc pst(i,n+2)=(y2-y1)*pst(i,jb)
1773 pst(i,n+2)= pst(i,jb)
1774 enddo
1775 co(5)=pf(n,n)
1776 co(6)= 2.*pf(n,n+1)
1777 co(7)=-2.*pf(n,n+2)
1778 co(8)=-2.*pf(n+1,n+2)
1779 if(ish.ge.8) then
1780 do i=n,n+2
1781 write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1782 enddo
1783 endif
1784 if(ish.ge.8)write(ifch,*) 'co right:',co(5),co(6),co(7),co(8)
1785
1786c.....lambda (absolute past).....................................
1787 n=nob+7
1788 do i=1,4
1789cc pst(i,n)= x1*pst(i,ib)+y1*pst(i,jb)
1790 pst(i,n)= 0.
1791 enddo
1792 ii=mmod(ib+1)
1793 do while (mmod(ii+jb).ne.0)
1794 if(ish.ge.8)write(ifch,*) 'add lambda',ii
1795 do i=1,4
1796 pst(i,n)=pst(i,n)+pst(i,ii)
1797 enddo
1798 ii=mmod(ii+1)
1799 enddo
1800 do i=1,4
1801cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1802cc pst(i,n+2)=(y2-y1)*pst(i,jb)
1803 pst(i,n+1)= pst(i,ib)
1804 pst(i,n+2)= pst(i,jb)
1805 enddo
1806 co(9)= pf(n,n)
1807 co(10)= 2.*pf(n,n+1)
1808 co(11)= 2.*pf(n,n+2)
1809 co(12)= 2.*pf(n+1,n+2)
1810 if(ish.ge.8)write(ifch,*)'co lambda:',co(9),co(10),co(11),co(12)
1811 call utprix('gaksco',ish,ishini,8)
1812 end
1813
1814c---------------------------------------------------------------------
1815 subroutine gakanp(ibr1,ibr,ibrr2,aml2,amr2,ala2,iok)
1816c---------------------------------------------------------------------
1817c mass adjustment of fragments
1818c ibr1-ibr ibr-ibrr2
1819c where ibr1,ibr,ibrr2 are break point indices
1820c aml2,amr2 are the reqired squared masses (if zero -> any mass)
1821c iok=0 (ok) or ok=1 (error)
1822c---------------------------------------------------------------------
1823 include 'epos.inc'
1824 parameter (mxnbr=500,mxnba=5000,mxnin=2000)
1825 common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
1826 $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
1827 & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
1828 double precision ax,ay,az,ae,am,bx,by,bz,be,A,B,C
1829 dimension co(12),am2(0:2),nin(2,0:mxnin)
1830 logical weiter
1831c pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1832c & -pst(2,i)*pst(2,j)-pst(1,i)*pst(1,j)
1833 mmod(i)=mod(mod(i,nob)+nob,nob)
1834 call utpri('gakanp',ish,ishini,6)
1835
1836 ibr2=ibrr2
1837 if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1838 iok=0
1839 ib=ijb(1,ibr)
1840 jb=ijb(2,ibr)
1841 ni=0
1842 10 do i=1,ni
1843 if((nin(1,i).eq.ib.and.nin(2,i).eq.jb)
1844 $ .or.(ipst(ib).eq.ipst(jb)))then
1845 iok=1
1846 if(ish.ge.4)then
1847 write(ifch,*) 'error ... endless loop'
1848 if(ib.eq.ipst(jb)) write(ifch,*) ' in zero mass region'
1849 endif
1850 goto 9999
1851 endif
1852 enddo
1853 ni=ni+1
1854 if(ni.gt.mxnin)stop'gakanp: increase parameter mxnin '
1855 nin(1,ni)=ib
1856 nin(2,ni)=jb
1857 if(ish.ge.6)write(ifch,*)
1858 if(ish.ge.6)write(ifch,*) 'ib,jb=',ib,jb
1859 if(ni.ge.2)then
1860 if(ish.ge.6)write(ifch,*)'rotate pt to new band'
1861 if(ish.ge.6)write(ifch,*)'from',nin(1,ni-1),nin(2,ni-1)
1862 if(ish.ge.6)write(ifch,*)' to',ib,jb
1863 if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1864 ax=pst(1,nin(1,ni-1))+pst(1,nin(2,ni-1))
1865 ay=pst(2,nin(1,ni-1))+pst(2,nin(2,ni-1))
1866 az=pst(3,nin(1,ni-1))+pst(3,nin(2,ni-1))
1867 ae=pst(4,nin(1,ni-1))+pst(4,nin(2,ni-1))
1868 am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1869 $ -dble(ay)**2-dble(az)**2))) !???????????????????????
1870 bx=pst(1,nin(1,ni-1))
1871 by=pst(2,nin(1,ni-1))
1872 bz=pst(3,nin(1,ni-1))
1873 be=pst(4,nin(1,ni-1))
1874 if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1875 call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1876 if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1877 call utlob4( 1,ax,ay,az,ae,am
1878 & ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
1879 if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1880 call utrot4( 1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
1881 if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1882 ax=pst(1,ib)+pst(1,jb)
1883 ay=pst(2,ib)+pst(2,jb)
1884 az=pst(3,ib)+pst(3,jb)
1885 ae=pst(4,ib)+pst(4,jb)
1886 am=sngl(sqrt(max(1d-8,dble(ae)**2-dble(ax)**2
1887 $ -dble(ay)**2-dble(az)**2))) !???????????????????????
1888 if(am.le.1.1e-4)then
1889 if(ish.ge.5)write(ifch,*)'error ... am<1.1e-4'
1890 iok=1
1891 goto 9999
1892 endif
1893 bx=pst(1,ib)
1894 by=pst(2,ib)
1895 bz=pst(3,ib)
1896 be=pst(4,ib)
1897 if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1898 if(ish.ge.6)write (ifch,*) 'ax,ay,az,ae',ax,ay,az,ae,am
1899 call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1900 if(ish.ge.6)write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1901 call utrot4(-1,bx,by,bz,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr))
1902 if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1903 call utlob4(-1,ax,ay,az,ae,am
1904 & ,ptb(1,ibr),ptb(2,ibr),ptb(3,ibr),ptb(4,ibr))
1905 if(ish.ge.6)write(ifch,'(4f12.8)') (ptb(i,ibr),i=1,4)
1906 endif
1907 call gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1908c if(ni.eq.1)print *,'------------------------------------',2
1909cc x=(xyb(1,ibr)-x1)/(x2-x1)
1910cc y=(xyb(2,ibr)-y1)/(y2-y1)
1911 x=xyb(1,ibr)
1912 y=xyb(2,ibr)
1913 am2(0)=aml2
1914 am2(1)=amr2
1915 am2(2)=ala2
1916 if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1917 if(amr2.le.0.)then
1918 l1=2
1919 l2=0
1920 elseif(aml2.le.0.)then
1921 l1=1
1922 l2=2
1923 elseif(ala2.le.0.)then
1924 l1=1
1925 l2=0
1926 else
1927 stop' not like this , please...'
1928 endif
1929 if(ish.ge.6.and.amr2.le.0)write(ifch,*) 'adjust: 1',l1,l2
1930 if(ish.ge.6.and.aml2.le.0)write(ifch,*) 'adjust: 2' ,l1,l2
1931 if(ish.ge.6.and.ala2.le.0)write(ifch,*) 'adjust: 3',l1,l2
1932 i=4*l1
1933 j=4*l2
1934 A = dble(co(i+4))*dble(co(j+3)) - dble(co(j+4))*dble(co(i+3))
1935 B = dble(co(i+4))*dble(co(j+1)) - dble(co(i+3))*dble(co(j+2))
1936 & + dble(co(i+2))*dble(co(j+3)) - dble(co(i+1))*dble(co(j+4))
1937 & - dble(am2(l2))*dble(co(i+4)) + dble(am2(l1))*dble(co(j+4))
1938 C = dble(co(i+2))*dble(co(j+1)) - dble(co(i+1))*dble(co(j+2))
1939 & + dble(am2(l1))*dble(co(j+2)) - dble(am2(l2))*dble(co(i+2))
1940 if (ish.ge.7) then
1941 write(ifch,*) 'ABC,q ',A,B,C,B**2-4.*A*C
1942 if(abs(A).gt.0.d0)then
1943 write(ifch,*) sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
1944 write(ifch,*) -sqrt(max(0d0,B**2-4d0*A*C))/2d0/A-B/A/2d0
1945 endif
1946 endif
1947 x=0.
1948 y=0.
1949 xx=0.
1950 yy=0.
1951 if(abs(A).gt.1.d-20.and.B*B-4.*A*C.ge.0.d0)then
1952 y=sngl(sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
1953 if(abs(co(i+2)+y*co(i+4)).gt.0.)
1954 & x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
1955 elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
1956 y=-sngl(C/B)
1957 if(abs(co(i+2)+y*co(i+4)).gt.0.)
1958 & x=(am2(l1)-co(i+1)-y*co(i+3))/(co(i+2)+y*co(i+4))
1959 else
1960 if(ish.ge.5)write(ifch,*)'error ... no solution of quadr equ'
1961 iok=1
1962 goto 9999
1963 endif
1964 if(abs(A).gt.1.d-20.and.B**2-4.*A*C.ge.0.d0)then
1965 yy=sngl(-sqrt(max(0.d0,B**2-4.*A*C))/2.d0/A-B/A/2.d0)
1966 if(abs(co(i+2)+yy*co(i+4)).gt.0.)
1967 & xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
1968 elseif(abs(A).le.1.d-20.and.abs(B).gt.0.d0)then
1969 yy=-sngl(C/B)
1970 if(abs(co(i+2)+yy*co(i+4)).gt.0.)
1971 & xx=(am2(l1)-co(i+1)-yy*co(i+3))/(co(i+2)+yy*co(i+4))
1972 else
1973 if(ish.ge.5)write(ifch,*)'error ... no solution (2) '
1974 iok=1
1975 goto 9999
1976 endif
1977 if(ish.ge.6)then
1978 write(ifch,*) x ,y ,(co(i+2)+ y*co(i+4)),' OK '
1979 write(ifch,*) xx,yy,(co(i+2)+yy*co(i+4)),' OK '
1980 endif
1981 weiter=.true.
1982 50 if(x.gt.x1.and.x.lt.x2.and.y.gt.y1.and.y.lt.y2)then
1983cc xyb(1,ibr)=x1+(x2-x1)*x
1984cc xyb(2,ibr)=y1+(y2-y1)*y
1985 xyb(1,ibr)=x
1986 xyb(2,ibr)=y
1987 ijb(1,ibr)=ib
1988 ijb(2,ibr)=jb
1989 e1=pst(4,nob+1)-x*pst(4,nob+2)+y*pst(4,nob+3)
1990 e2=pst(4,nob+4)+x*pst(4,nob+5)-y*pst(4,nob+6)
1991 if( e1.lt.0. .or. e2.lt.0. ) then
1992 if(ish.ge.5)write(ifch,*)'error ... e1<0 or e2<0'
1993 iok=1
1994 goto 9999
1995 endif
1996 !amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
1997 !amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
1998 if(ish.ge.6)then
1999 amal2=co(1)+co(2)*x+co(3)*y+co(4)*x*y
2000 amar2=co(5)+co(6)*x+co(7)*y+co(8)*x*y
2001 write(ifch,*) 'brkshift:',xyb(1,ibr),xyb(2,ibr),ib,jb
2002 write (ifch,*)'E:',e1
2003 write (ifch,*)'E:',e2
2004 write(ifch,'(2(a6,1g12.6))') 'aml:'
2005 & ,sqrt(max(0.,amal2)),'amr:',sqrt(max(0.,amar2))
2006 endif
2007 i=ibr1+1
2008 do while(i.le.ibr-1)
2009 if((mmod(ijb(1,i)+nob/2) .lt. mmod(ijb(1,ibr)+nob/2)
2010 & .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2011 & .and.(xyb(1,i).gt.xyb(1,ibr))))
2012 & .and.
2013 & (ijb(2,i) .gt. ijb(2,ibr)
2014 & .or.(ijb(2,i) .eq. ijb(2,ibr)
2015 & .and.xyb(2,i).lt.xyb(2,ibr)))) goto 150
2016 if(ish.ge.6) then
2017 write(ifch,*) 'away:'
2018 & ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2019 endif
2020 call gakdel(i)
2021 i=i-1
2022 ibr=ibr-1
2023 ibr2=ibr2-1
2024 150 i=i+1
2025 enddo
2026 i=ibr+1
2027 do while (i.le.ibr2-1)
2028 if((mmod(ijb(1,i)+nob/2) .gt. mmod(ijb(1,ibr)+nob/2)
2029 & .or.(mmod(ijb(1,i)+nob/2) .eq. mmod(ijb(1,ibr)+nob/2)
2030 & .and.(xyb(1,i).lt.xyb(1,ibr))))
2031 & .and.
2032 & (ijb(2,i) .lt. ijb(2,ibr)
2033 & .or.(ijb(2,i) .eq. ijb(2,ibr)
2034 & .and.xyb(2,i).gt.xyb(2,ibr)))) goto 160
2035 if(ish.ge.6) then
2036 write(ifch,*) 'away:'
2037 & ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
2038 endif
2039 call gakdel(i)
2040 ibr2=ibr2-1
2041 i=i-1
2042 160 i=i+1
2043 enddo
2044 goto 9999
2045 else
2046 if(x.gt.x2
2047 & .and.ib.ne.ijb(1,ibr1) !brk-begrenzung
2048 & .and.mmod(ib-1).ne.jb !linker oder rechter Rand
2049 & .and.mmod(ib-1+jb).ne.0)then !oben oder unten
2050 ib=mmod(ib-1)
2051 goto 10
2052 endif
2053 if(x.lt.x1
2054 & .and.ib.ne.ijb(1,ibr2) !brk-begrenzung
2055 & .and.mmod(ib+1).ne.jb !linker oder rechter Rand
2056 & .and.mmod(ib+1+jb).ne.0)then !oben oder unten
2057 ib=mmod(ib+1)
2058 goto 10
2059 endif
2060 if(y.gt.y2
2061 & .and.jb.ne.ijb(2,ibr2) !brk-begrenzung
2062 & .and.mmod(jb-1).ne.ib !linker oder rechter Rand
2063 & .and.mmod(jb-1+ib).ne.0)then !oben oder unten
2064 jb=mmod(jb-1)
2065 goto 10
2066 endif
2067 if(y.lt.y1
2068 & .and.jb.ne.ijb(2,ibr1) !brk-begrenzung
2069 & .and.mmod(jb+1).ne.ib !linker oder rechter Rand
2070 & .and.mmod(jb+1+ib).ne.0)then !oben oder unten
2071 jb=mmod(jb+1)
2072 goto 10
2073 endif
2074 if(weiter)then
2075 weiter=.false.
2076 x=xx
2077 y=yy
2078 goto 50
2079 endif
2080 if(ish.ge.5)write(ifch,*)'error ... x,y not in allowed range'
2081 iok=1
2082 endif
2083 9999 if(amr2.eq.0.) ibrr2=ibr2
2084 call utprix('gakanp',ish,ishini,6)
2085 end
2086
2087cc----------------------------------------------------------------------
2088c subroutine gakstr(ifl)
2089cc----------------------------------------------------------------------
2090cc
2091cc calculates string-fragments by taking off pt of breakup
2092cc do with ifl=1 undo with ifl=-1
2093cc
2094cc----------------------------------------------------------------------
2095c include 'epos.inc'
2096c common /cpptr/ pptr(4,mxptl),ystr(mxptl)
2097c
2098c do i=1,nptl
2099c if(istptl(i).eq.29)then
2100c nk1=ifrptl(1,i)
2101c nk2=ifrptl(2,i)
2102c do j=nk1,nk2
2103c if ((istptl(j).eq.0.or.istptl(j-1).eq.0).and.j.ne.nk1) then
2104c do k=1,4
2105c pptl(k,j)=pptl(k,j)+pptr(k,j-1)*ifl
2106c enddo
2107c !write(ifch,*)"left side back to ",j,(pptr(k,j-1),k=1,4)
2108c endif
2109c if ((istptl(j).eq.0.or.istptl(j+1).eq.0).and.j.ne.nk2) then
2110c do k=1,4
2111c pptl(k,j)=pptl(k,j)-pptr(k,j)*ifl
2112c enddo
2113c !write(ifch,*)"right side back to ",j,(-pptr(k,j),k=1,4)
2114c endif
2115c if(ifl.eq.-1.and.istptl(j).eq.0)then
2116c e=pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
2117c & +pptl(5,j)**2
2118c e=sqrt(e)
2119c !dif=abs(e-pptl(4,j))
2120c !if(dif.gt.0.01.and.dif/e.gt.0.1)print*,j,e,pptl(4,j)
2121c pptl(4,j)=e
2122c endif
2123c enddo
2124c endif
2125c if(istptl(i).eq.0)then
2126c if ( ifl.eq.1 ) then
2127c ystr(i)=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
2128c * /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2) )
2129c endif
2130c endif
2131c enddo
2132c end
2133
2134c----------------------------------------------------------------------
2135 subroutine gakli2(nn1,nn2)
2136c----------------------------------------------------------------------
2137
2138 include 'epos.inc'
2139 double precision pgampr,rgampr
2140 common/cgampr/pgampr(5),rgampr(4)
2141c double precision db1,db2,db3,db4,db5
2142 character label*8,idlabl*8
2143c db1=0d0
2144c db2=0d0
2145c db3=rgampr(4)
2146c db4=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
2147c db5=sqrt( db4**2-rgampr(4)**2)
2148 n1=nn1
2149 n2=nn2
2150 if (n1.eq.0) n1=1
2151 if (n2.eq.0) n2=nptl
2152 write (ifch,'(1a4,5a12,4a4,a10,2a4)')
2153 &'no.','px','py','pz','E','m','ior','jor','if1','if2'
2154 &,'id','ist','ity'
2155 do i=n1,n2
2156 if (idptl(i).lt.10000)then
2157 label=' '
2158 label=idlabl(idptl(i))
2159 else
2160 label=' '
2161 endif
2162 chrg=0.
2163 if(iabs(idptl(i)).le.9999)call idchrg(idptl(i),chrg)
2164 write (ifch,125) i,(pptl(j,i),j=1,5),iorptl(i),jorptl(i)
2165 & ,ifrptl(1,i),ifrptl(2,i),idptl(i)
2166 $ ,chrg !charge
2167 & ,istptl(i),ityptl(i),label
2168 enddo
2169 125 format (1i4,5g18.10,4i6,1i10
2170 $ ,1f5.2 !charge
2171 $ ,2i4,' ',A8
2172c $ ,7g12.4,i5
2173 $ )
2174 126 format (A4,7g12.4,i5)
2175 return
2176 end
2177
2178c----------------------------------------------------------------------
2179c subroutine gakli4
2180cc----------------------------------------------------------------------
2181c
2182c include 'epos.inc'
2183c parameter (mxnbr=500,mxnba=5000)
2184c common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2185c $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2186c & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2187c dimension co(12)
2188c do i=0,nob-1
2189c write (ifch,10) i,(pst(j,i),j=1,4)
2190c enddo
2191c 10 format(1i4,5g18.10)
2192c return
2193c end
2194c
2195cc----------------------------------------------------------------------
2196c subroutine gakli3
2197cc----------------------------------------------------------------------
2198c
2199c include 'epos.inc'
2200c parameter (mxnbr=500,mxnba=5000)
2201c common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2202c $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2203c & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2204c dimension co(12),p1(5),p2(5)
2205c
2206c write(ifch,*) 'particle list of string decay'
2207c do i=1,5
2208c p1(i)=0.
2209c enddo
2210c do i=1,nbr+1
2211c if(i.lt.nbr+1)then
2212c call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
2213c if(x2.gt.x1)then
2214ccc x=(xyb(1,i)-x1)/(x2-x1)
2215c x=xyb(1,i)
2216c else
2217c x=0.
2218c endif
2219c if(y2.gt.y1)then
2220ccc y=(xyb(2,i)-y1)/(y2-y1)
2221c y=xyb(2,i)
2222c else
2223c y=0.
2224c endif
2225c aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
2226c amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
2227c aml=sign(sqrt(abs(aml2)),aml2)
2228c amr=sign(sqrt(abs(amr2)),amr2)
2229c do j=1,4
2230c p2(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
2231c p1(j)=p1(j)+p2(j)
2232c enddo
2233c p2(5)=aml
2234c else
2235c do j=1,4
2236c p2(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
2237c p1(j)=p1(j)+p2(j)
2238c enddo
2239c p2(5)=amr
2240c endif
2241c write(ifch,'(2i4,i6,a,i5,i10,5g14.6)') i-1,i
2242c & ,-iflb(i-1),'==',iflb(i),ip(i)
2243c & ,(p2(j),j=1,5)
2244c enddo
2245c am2=p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2
2246c p1(5)=sign(sqrt(abs(am2)),am2)
2247c write(ifch,'(12x,a60)')
2248c & '------------------------------------------------------------'
2249c write(ifch,'(14x,5g14.6)') (p1(j),j=1,5)
2250c write(ifch,*)
2251c
2252c end
2253c
2254
2255c---------------------------------------------------------------------
2256 subroutine idress(id,am,idr,iadj)
2257c---------------------------------------------------------------------
2258 include 'epos.inc'
2259 call idres(id,am,idr,iadj,1)
2260 if(idr.eq.0)then
2261 return
2262 endif
2263 ids=max(mod(iabs(id)/100,10),mod(iabs(id)/10,10))
2264 if(iabs(idr).le.999) then
2265c write (*,*) ' ',id,idr
2266 if(ids.le.3)return
2267 if(ids.le.2)then
2268 idr=sign(iabs(id)+int(rangen()+0.5),id)
2269 elseif(ids.eq.3)then
2270 idr=sign(iabs(id)+int(rangen()+0.6),id)
2271 else
2272 idr=sign(iabs(id)+int(rangen()+0.75),id)
2273 endif
2274c write (*,*) '->',id,idr
2275 call idmass(idr,am)
2276 elseif(iabs(idr).le.9999)then
2277 if(ids.le.3)return
2278 if(mod(iabs(idr),10).gt.1)then
2279 if(iabs(id).ne.1111.and.iabs(id).ne.2221.and.iabs(id).ne.3331)
2280 $ then
2281 idr=sign(iabs(id)+1,id)
2282 call idmass(idr,am)
2283 else
2284 idr=id
2285 call idmass(idr,am)
2286 endif
2287 endif
2288 endif
2289
2290 end
2291
2292
2293c---------------------------------------------------------------------
2294 SUBROUTINE gaksphe(sphe,r,mstu41)
2295
2296C...Purpose: to perform sphericity tensor analysis to give sphericity,
2297C...aplanarity and the related event axes. stolen from jetset ;-)
2298 include 'epos.inc'
2299 dimension sphe(4,3)
2300 DIMENSION SM(3,3),SV(3,3)
2301
2302C...Calculate matrix to be diagonalized.
2303 NP=0
2304c MSTU41=2
2305 DO 110 J1=1,3
2306 DO 100 J2=J1,3
2307 SM(J1,J2)=0.
2308 100 CONTINUE
2309 110 CONTINUE
2310 PS=0.
2311 DO 140 I=1,nptl
2312 IF(istptl(i).ne.0) GOTO 140
2313 IF(MSTU41.GE.2) THEN
2314 ida=iabs(idptl(i))
2315 IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15) GOTO 140
2316 IF(MSTU41.GE.3) then
2317 call idchrg(idptl(i),chrg)
2318 if (abs(chrg).le.0.1) goto 140
2319 endif
2320 ENDIF
2321 NP=NP+1
2322 PA=SQRT(pptl(1,i)**2+pptl(2,I)**2+pptl(3,i)**2)
2323 PWT=1.
2324 IF(ABS(r-2.).GT.0.001) PWT=MAX(1E-10,PA)**(r-2.)
2325 DO 130 J1=1,3
2326 DO 120 J2=J1,3
2327 SM(J1,J2)=SM(J1,J2)+PWT*pptl(j1,i)*pptl(j2,i)
2328 120 CONTINUE
2329 130 CONTINUE
2330 PS=PS+PWT*PA**2
2331 140 CONTINUE
2332
2333C...Very low multiplicities (0 or 1) not considered.
2334 IF(NP.LE.1) THEN
2335 if(ish.ge.1)then
2336 CALL utmsg('sphe ')
2337 write(ifch,*) 'too few particles for analysis'
2338 call utmsgf
2339 endif
2340 sphe(4,1)=-1.
2341 RETURN
2342 ENDIF
2343 DO 160 J1=1,3
2344 DO 150 J2=J1,3
2345 SM(J1,J2)=SM(J1,J2)/PS
2346 150 CONTINUE
2347 160 CONTINUE
2348
2349C...Find eigenvalues to matrix (third degree equation).
2350 SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
2351 & SM(1,3)**2-SM(2,3)**2)/3.-1./9.
2352 SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
2353 & SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))
2354 & +SM(1,2)*SM(1,3)*SM(2,3)+1./27.
2355 SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
2356 sphe(4,1)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
2357 sphe(4,3)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
2358 sphe(4,2)=1.-sphe(4,1)-sphe(4,3)
2359 IF(sphe(4,2).LT.1E-5) THEN
2360 if(ish.ge.1)then
2361 CALL utmsg('gaksphe')
2362 write(ifch,*) 'all particles back-to-back'
2363 call utmsgf
2364 endif
2365 sphe(4,1)=-1.
2366 RETURN
2367 ENDIF
2368
2369C...Find first and last eigenvector by solving equation system.
2370 DO 240 I=1,3,2
2371 DO 180 J1=1,3
2372 SV(J1,J1)=SM(J1,J1)-sphe(4,I)
2373 DO 170 J2=J1+1,3
2374 SV(J1,J2)=SM(J1,J2)
2375 SV(J2,J1)=SM(J1,J2)
2376 170 CONTINUE
2377 180 CONTINUE
2378 SMAX=0.
2379 DO 200 J1=1,3
2380 DO 190 J2=1,3
2381 IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
2382 JA=J1
2383 JB=J2
2384 SMAX=ABS(SV(J1,J2))
2385 190 CONTINUE
2386 200 CONTINUE
2387 SMAX=0.
2388 DO 220 J3=JA+1,JA+2
2389 J1=J3-3*((J3-1)/3)
2390 RL=SV(J1,JB)/SV(JA,JB)
2391 DO 210 J2=1,3
2392 SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
2393 IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
2394 JC=J1
2395 SMAX=ABS(SV(J1,J2))
2396 210 CONTINUE
2397 220 CONTINUE
2398 JB1=JB+1-3*(JB/3)
2399 JB2=JB+2-3*((JB+1)/3)
2400 sphe(JB1,I)=-SV(JC,JB2)
2401 sphe(jb2,I)=SV(JC,JB1)
2402 sphe(jb,I)=-(SV(JA,JB1)*sphe(jb1,I)+SV(JA,JB2)
2403 & *sphe(jb2,I))/SV(JA,JB)
2404 PA=SQRT(sphe(1,I)**2+sphe(2,I)**2+sphe(3,I)**2)
2405 SGN=(-1.)**INT(rangen()+0.5)
2406 DO 230 J=1,3
2407 sphe(j,I)=SGN*sphe(j,I)/PA
2408 230 CONTINUE
2409 240 CONTINUE
2410
2411C...Middle axis orthogonal to other two. Fill other codes.
2412 SGN=(-1.)**INT(rangen()+0.5)
2413 sphe(1,2)=SGN*(sphe(2,1)*sphe(3,3)-sphe(3,1)*sphe(2,3))
2414 sphe(2,2)=SGN*(sphe(3,1)*sphe(1,3)-sphe(1,1)*sphe(3,3))
2415 sphe(3,2)=SGN*(sphe(1,1)*sphe(2,3)-sphe(2,1)*sphe(1,3))
2416
2417 do i=1,3
2418 do j=1,4
2419 pptl(j,nptl+i)=sphe(j,i)
2420 enddo
2421 enddo
2422
2423
2424C...Calculate sphericity and aplanarity. Select storing option.
2425ccc SPH=1.5*(sphe(4,2)+sphe(4,3))
2426ccc APL=1.5*sphe(4,3)
2427
2428 RETURN
2429 END
2430
2431C*********************************************************************
2432
2433 SUBROUTINE gakthru(thru,mstu41)
2434
2435C...Purpose: to perform thrust analysis to give thrust, oblateness
2436C...and the related event axes. stolen from jetset ;-)
2437 include 'epos.inc'
2438 DIMENSION TDI(3),TPR(3)
2439 dimension thru(4,3)
2440
2441C...Take copy of particles that are to be considered in thrust analysis.
2442 NP=0
2443 PS=0.
2444c MSTU41=2
2445 MSTU44=4
2446 MSTU45=2
2447 PARU42=1.0
2448 PARU48=0.0000001
2449 DO 100 I=1,nptl
2450 IF(istptl(i).ne.0)goto 100
2451 ida=iabs(idptl(i))
2452 IF(ida.EQ.0.OR.ida.EQ.11.OR.ida.EQ.13.OR.ida.EQ.15)GOTO 100
2453 IF(MSTU41.GE.3) then
2454 call idchrg(idptl(i),chrg)
2455 if (abs(chrg).le.0.1) goto 100
2456 endif
2457
2458 IF(nptl+NP.GE.mxptl) THEN
2459 CALL utstop('gakthru: no more memory left in cptl&')
2460 thru(4,1)=-1.
2461 RETURN
2462 ENDIF
2463 NP=NP+1
2464c K(nptl+NP,1)=23
2465 pptl(1,nptl+NP)=pptl(1,I)
2466 pptl(2,nptl+NP)=pptl(2,I)
2467 pptl(3,nptl+NP)=pptl(3,I)
2468 pptl(4,nptl+NP)=SQRT(pptl(1,I)**2+pptl(2,I)**2+pptl(3,I)**2)
2469 pptl(5,nptl+NP)=1.
2470 IF(ABS(PARU42-1.).GT.0.001)
2471 & pptl(5,nptl+NP)=pptl(4,nptl+NP)**(PARU42-1.)
2472 PS=PS+pptl(4,nptl+NP)*pptl(5,nptl+NP)
2473 100 CONTINUE
2474
2475C...Very low multiplicities (0 or 1) not considered.
2476 IF(NP.LE.1) THEN
2477 CALL utmsg('thru ')
2478 write(ifch,*) 'too few particles for analysis'
2479 call utmsgf
2480 thru(4,1)=-1
2481 RETURN
2482 ENDIF
2483
2484C...Loop over thrust and major. T axis along z direction in latter case.
2485 DO 320 ILD=1,2
2486 IF(ILD.EQ.2) THEN
2487c PHI=GAKANG(pptl(1,nptl+NP+1),pptl(2,nptl+NP+1))
2488c CALL lurot(nptl+1,nptl+NP+1,0.,-PHI)
2489c THE=GAKANG(pptl(3,nptl+NP+1),pptl(1,nptl+NP+1))
2490c CALL lurot(nptl+1,nptl+NP+1,-THE,0.)
2491 ax=pptl(1,nptl+NP+1)
2492 ay=pptl(2,nptl+NP+1)
2493 az=pptl(3,nptl+NP+1)
2494 do irot=nptl+1,nptl+NP+1
2495 call utrota(1,ax,ay,az
2496 & ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2497 enddo
2498 if(np.eq.2)then
2499 pptl(1,nptl+NP+2)=1.
2500 pptl(2,nptl+NP+2)=0.
2501 pptl(3,nptl+NP+2)=0.
2502 pptl(4,nptl+NP+2)=0.
2503 goto 325
2504 endif
2505 ENDIF
2506
2507C...Find and order particles with highest p (pT for major).
2508 DO 110 ILF=nptl+NP+4,nptl+NP+MSTU44+4
2509 pptl(4,ILF)=0.
2510 110 CONTINUE
2511 DO 160 I=nptl+1,nptl+NP
2512 IF(ILD.EQ.2) pptl(4,I)=SQRT(pptl(1,I)**2+pptl(2,I)**2)
2513 DO 130 ILF=nptl+NP+MSTU44+3,nptl+NP+4,-1
2514 IF(pptl(4,I).LE.pptl(4,ILF)) GOTO 140
2515 DO 120 J=1,5
2516 pptl(j,ILF+1)=pptl(j,ILF)
2517 120 CONTINUE
2518 130 CONTINUE
2519 ILF=nptl+NP+3
2520 140 DO 150 J=1,5
2521 pptl(j,ILF+1)=pptl(j,I)
2522 150 CONTINUE
2523 160 CONTINUE
2524
2525C...Find and order initial axes with highest thrust (major).
2526 DO 170 ILG=nptl+NP+MSTU44+5,nptl+NP+MSTU44+15
2527 pptl(4,ILG)=0.
2528 170 CONTINUE
2529 NC=2**(MIN(MSTU44,NP)-1)
2530 DO 250 ILC=1,NC
2531 DO 180 J=1,3
2532 TDI(J)=0.
2533 180 CONTINUE
2534 DO 200 ILF=1,MIN(MSTU44,NP)
2535 SGN=pptl(5,nptl+NP+ILF+3)
2536 IF(2**ILF*((ILC+2**(ILF-1)-1)/2**ILF).GE.ILC) SGN=-SGN
2537 DO 190 J=1,4-ILD
2538 TDI(J)=TDI(J)+SGN*pptl(j,nptl+NP+ILF+3)
2539 190 CONTINUE
2540 200 CONTINUE
2541 TDS=TDI(1)**2+TDI(2)**2+TDI(3)**2
2542 DO 220 ILG=nptl+NP+MSTU44+MIN(ILC,10)+4,nptl+NP+MSTU44+5,-1
2543 IF(TDS.LE.pptl(4,ILG)) GOTO 230
2544 DO 210 J=1,4
2545 pptl(j,ILG+1)=pptl(j,ILG)
2546 210 CONTINUE
2547 220 CONTINUE
2548 ILG=nptl+NP+MSTU44+4
2549 230 DO 240 J=1,3
2550 pptl(j,ILG+1)=TDI(J)
2551 240 CONTINUE
2552 pptl(4,ILG+1)=TDS
2553 250 CONTINUE
2554
2555C... Iterate direction of axis until stable maximum.
2556 pptl(4,nptl+NP+ILD)=0.
2557 ILG=0
2558 260 ILG=ILG+1
2559 THP=0.
2560 270 THPS=THP
2561 DO 280 J=1,3
2562 IF(THP.LE.1E-10) TDI(J)=pptl(j,nptl+NP+MSTU44+4+ILG)
2563 IF(THP.GT.1E-10) TDI(J)=TPR(J)
2564 TPR(J)=0.
2565 280 CONTINUE
2566 DO 300 I=nptl+1,nptl+NP
2567 SGN=SIGN(pptl(5,I),TDI(1)*pptl(1,I)
2568 & +TDI(2)*pptl(2,I)+TDI(3)*pptl(3,I))
2569 DO 290 J=1,4-ILD
2570 TPR(J)=TPR(J)+SGN*pptl(j,I)
2571 290 CONTINUE
2572 300 CONTINUE
2573 THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
2574 IF(THP.GE.THPS+PARU48) GOTO 270
2575
2576C... Save good axis. Try new initial axis until a number of tries agree.
2577 IF(THP.LT.pptl(4,nptl+NP+ILD)-PARU48.AND.ILG.LT.MIN(10,NC))
2578 & GOTO 260
2579 IF(THP.GT.pptl(4,nptl+NP+ILD)+PARU48)
2580 $ THEN
2581 IAGR=0
2582 SGN=(-1.)**INT(rangen()+0.5)
2583 DO 310 J=1,3
2584 pptl(j,nptl+NP+ILD)=SGN*TPR(J)/(PS*THP)
2585 310 CONTINUE
2586 pptl(4,nptl+NP+ILD)=THP
2587 pptl(5,nptl+NP+ILD)=0.
2588 ENDIF
2589 IAGR=IAGR+1
2590 IF(IAGR.LT.MSTU45.AND.ILG.LT.MIN(10,NC)) GOTO 260
2591 320 CONTINUE
2592
2593C... Find minor axis and value by orthogonality.
2594 325 SGN=(-1.)**INT(rangen()+0.5)
2595 pptl(1,nptl+NP+3)=-SGN*pptl(2,nptl+NP+2)
2596 pptl(2,nptl+NP+3)=SGN*pptl(1,nptl+NP+2)
2597 pptl(3,nptl+NP+3)=0.
2598 THP=0.
2599 DO 330 I=nptl+1,nptl+NP
2600 THP=THP+pptl(5,I)*ABS(pptl(1,nptl+NP+3)*pptl(1,I)
2601 & +pptl(2,nptl+NP+3)*pptl(2,I))
2602 330 CONTINUE
2603 pptl(4,nptl+NP+3)=THP/PS
2604 pptl(5,nptl+NP+3)=0.
2605
2606
2607C... Fill axis information. Rotate back to original coordinate system.
2608 do irot=nptl+NP+1,nptl+NP+3
2609 call utrota(-1,ax,ay,az
2610 & ,pptl(1,irot),pptl(2,irot),pptl(3,irot))
2611 enddo
2612
2613 do ild=1,3
2614 do j=1,4
2615 thru(j,ild)=pptl(j,nptl+NP+ild)
2616 enddo
2617 enddo
2618
2619C...Calculate thrust and oblateness. Select storing option.
2620ccc THR=thru(4,1)
2621ccc OBL=thru(4,2)-thru(4,3)
2622
2623 RETURN
2624 END
2625
2626 subroutine gakjet(ijadu)
2627 include 'epos.inc'
2628 common/njjjj/njets(5,2,mxbins)
2629c nmin=xpar1
2630c nmax=xpar2
2631 if(nrevt.eq.1)then
2632 do i=1,5
2633 do j=1,mxbins
2634 njets(i,ijadu,j)=0
2635 enddo
2636 enddo
2637 endif
2638 do i=1,nrbins
2639 ycut=xminim*(xmaxim/xminim)**((real(i)-0.5)/nrbins)
2640c if(iologe.ne.1)ycut=xminim+(xmaxim-xminim)/nrbins*(nrbins)
2641 nj=gaknjt(ycut,ijadu)
2642 if(nj.ge.1.and.nj.le.5)then
2643 njets(nj,ijadu,i)=njets(nj,ijadu,i)+1
2644 endif
2645 enddo
2646 end
2647
2648 subroutine gakjto
2649 include 'epos.inc'
2650 common/njjjj/njets(5,2,mxbins)
2651 n=xpar4
2652 ijadu=xpar3
2653 do i=1,nrbins
2654 ycut=xminim*(xmaxim/xminim)**((real(i)-0.5)/nrbins)
2655 write (ifhi,*) ycut,real(njets(n,ijadu,i))/nrevt
2656 $ ,sqrt(1.*njets(n,ijadu,i))/nrevt
2657 enddo
2658 end
2659C*********************************************************************
2660 function gaknjt(ycut,ijadu)
2661c
2662c ijadu 1 = JADE 2=DURHAM
2663c ycut - max. distance
2664c
2665 include 'epos.inc'
2666
2667 a2j(i,j)=2.*pptl(4,i)*pptl(4,j)*(1.-(pptl(1,i)*pptl(1,j)
2668 & +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2669 & /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2670 & *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2671
2672 a2d(i,j)=2.*min(pptl(4,i)**2,pptl(4,j)**2)
2673 & *(1.-(pptl(1,i)*pptl(1,j)
2674 & +pptl(2,i)*pptl(2,j)+pptl(3,i)*pptl(3,j))
2675 & /(sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2676 & *sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2)))/evis**2
2677
2678 bet(i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2679
2680 ska(i,j)=pptl(1,i)*pptl(1,j)+pptl(2,i)*pptl(2,j)
2681 & +pptl(3,i)*pptl(3,j)
2682
2683 a2c(i,j)= ((bet(i)*bet(j)-ska(i,j))
2684 & *2.*bet(i)*bet(j) / (0.00001+bet(i)+bet(j))**2 )
2685
2686 evis=0.
2687 nn=0
2688 do i=1,nptl
2689 if (istptl(i).eq.0) then
2690 nn=nn+1
2691 do j=1,5
2692 pptl(j,nptl+nn)=pptl(j,i)
2693 enddo
2694 evis=evis+pptl(4,i)
2695 jorptl(i)=nn
2696 endif
2697 enddo
2698 iflag=0
2699 do while (iflag.eq.0.and.nn.ge.2)
2700 a2min=ycut
2701 iflag=1
2702 do i=nptl+1,nptl+nn-1
2703 do j=i+1,nptl+nn
2704 if(ijadu.eq.1)then
2705 a2=a2j(i,j)
2706 elseif(ijadu.eq.2) then
2707 a2=a2d(i,j)
2708 else
2709 a2=a2c(i,j)
2710 endif
2711 if (a2.lt.a2min) then
2712 iflag=0
2713 i1=i
2714 j1=j
2715 a2min=a2
2716 endif
2717 enddo
2718 enddo
2719 if (iflag.eq.0) then
2720 do i=1,4
2721 pptl(i,i1)=pptl(i,i1)+pptl(i,j1)
2722 enddo
2723 do i=1,nptl
2724 if(istptl(i).eq.0.and.jorptl(i).eq.j1-nptl)
2725 & jorptl(i)=i1-nptl
2726 if(istptl(i).eq.0.and.jorptl(i)+nptl.gt.j1)
2727 & jorptl(i)=jorptl(i)-1
2728 enddo
2729 do i=j1,nptl+nn
2730 do j=1,5
2731 pptl(j,i)=pptl(j,i+1)
2732 enddo
2733 enddo
2734 nn=nn-1
2735 endif
2736 enddo
2737 do i=nptl+1,nptl+nn
2738 istptl(i)=-1
2739 jorptl(i)=i-nptl
2740 pptl(5,i)=sqrt(max(0.,(pptl(4,i)+pptl(3,i))
2741 & *(pptl(4,i)-pptl(3,i))-pptl(2,i)**2-pptl(1,i)**2))
2742 enddo
2743 do i=nptl+1,nptl+nn-1
2744 do j=i+1,nptl+nn
2745 if(pptl(4,jorptl(i)+nptl).lt.pptl(4,jorptl(j)+nptl))then
2746 k=jorptl(i)
2747 jorptl(i)=jorptl(j)
2748 jorptl(j)=k
2749 endif
2750 enddo
2751 enddo
2752 do i=nptl+1,nptl+nn
2753 idptl(nptl+jorptl(i))=9910+i-nptl
2754 enddo
2755 do i=1,nptl
2756 jorptl(i)=0 !jorptl(nptl+jorptl(i))
2757 enddo
2758c nptl=nptl+nn
2759
2760 gaknjt=nn
2761 return
2762 end
2763
2764 subroutine idtr5(id,ic)
2765 integer ic(2)
2766 ic(1)=0
2767 ic(2)=0
2768 ii=1
2769 if(id.lt.0)ii=2
2770 i1=1
2771 if(iabs(id).gt.999)i1=3
2772 do i=i1,int(log(abs(real(id)))/log(10.))+1
2773 j=mod(iabs(id)/10**(i-1),10)
2774 if(j.gt.0)then
2775 ic(ii)=ic(ii)+10**(6-j)
2776 endif
2777 enddo
2778 return
2779 end
2780
2781 function ammin(id1,id2)
2782 dimension ic1(2),ic2(2),jc2(6,2),jc1(6,2)
2783 call idtr5(id1,ic1)
2784 call idtr5(id2,ic2)
2785 call idcomk(ic1)
2786 call idcomk(ic2)
2787 call iddeco(ic1,jc1)
2788 call iddeco(ic2,jc2)
2789 ammin=utamnx(jc1,jc2)
2790 end
2791
2792
2793c function idtr(id1,id2)
2794c integer ic(2),id(2)
2795c id(1)=id1
2796c id(2)=id2
2797c do i=1,2
2798c ic(i)=0
2799c enddo
2800c do j=1,2
2801c ii=1
2802c if(id(j).lt.0)ii=2
2803c i1=1
2804c if(iabs(id(j)).gt.999)i1=3
2805c do i=i1,int(log(abs(real(id(j))))/log(10.))+1
2806c jj=mod(iabs(id(j))/10**(i-1),10)
2807c if(jj.gt.0)then
2808c ic(ii)=ic(ii)+10**(6-jj)
2809c endif
2810c enddo
2811c enddo
2812c idtr=idtra(ic,0,0,3)
2813c if(idtr.ne.idsp(id1,id2))then
2814c write (*,'(4i6)') idtr,idsp(id1,id2),id1,id2
2815c endif
2816c return
2817c end
2818
2819 function idsp(id1,id2)
2820 ia1=iabs(id1)
2821 ia2=iabs(id2)
2822 if(ia1.ge.1000.and.ia2.ge.1000)then
2823 idsp=0
2824 isign=0
2825 elseif(ia1.le.1000.and.ia2.le.1000)then
2826 idsp=min(ia1,ia2)*100+max(ia1,ia2)*10
2827 isign=1
2828 if(max(ia1,ia2).ne.-min(id1,id2)) isign = -1
2829 if(idsp.eq.220)idsp=110
2830 if(idsp.eq.330)idsp=220
2831 else
2832 isign=1
2833 if(id1.lt.0.and.id2.lt.0)isign=-1
2834 idb=min(ia1,ia2)
2835 if(idb.eq.5)then
2836 idsp=0
2837 return
2838 endif
2839 ida=max(ia1,ia2)
2840 ida1=ida/1000
2841 ida2=mod(ida/100,10)
2842 if(idb.le.ida1)then
2843 idsp=idb*1000+ida/10
2844 elseif(idb.le.ida2)then
2845 idsp=ida1*1000+idb*100+ida2*10
2846 else
2847 idsp=ida+idb*10
2848 endif
2849 if(ida1.eq.ida2.and.ida2.eq.idb)idsp=idsp+1
2850 endif
2851 idsp=idsp*isign
2852 return
2853 end
2854