1 c-----------------------------------------------------------------------
2 subroutine gakfra(iret)
3 c-----------------------------------------------------------------------
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
13 c 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
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)
28 call utpri('gakfra',ish,ishini,4)
30 c.....search string in pptl-list...................................
34 if(nclean.gt.0.and.nptl.gt.mxptl/2)then
36 do iii=maproj+matarg+1,nptl
38 if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
39 if(go.and.mod(istptl(iii),10).ne.0)then
44 if(nnnpt.gt.mxptl-nptl)then
46 call utclea(maproj+matarg+1,nptl0)
49 c write(ifch,'(a)')' ---------after utclea------------'
50 c call blist('&',nk1,nkmax)
54 do while(nk1.le.nkmax)
55 if(istptl(nk1).eq.20)then
57 do while(idptl(nk2).eq.9)
60 goto 3 !ok, string from nk1 to nk2
69 goto 9999 !no more string
72 c.......................................................................
73 c decay string nk1 - nk2
74 c.......................................................................
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)
85 if(nk2-nk1.gt.1.or.iappl.eq.6)kky=1
88 if(itp.ge.50.and.itp.le.59.or.itp.ge.40.and.itp.le.49)krm=1
90 if(itp.eq.43.or.itp.eq.53)kdrm=1
91 c special behaviour of remnant of reactions with lost Pomerons
92 if(itp.eq.47.or.itp.eq.57)then
109 zzzz=zpaptl(nk1)+zpaptl(nk2)
112 zzid=min( 1+zidmax , 1+zidinc *zzzz ) !<-----
113 zzip=min( 1+zopmax , 1+zopinc *zzzz ) !<-----
114 c zzrm=min( zrmmax , zrminc *zzzz ) !<-----
115 if(itp.ge.40.and.itp.le.49)then
117 elseif(itp.ge.50.and.itp.le.59)then
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 ) !<-----
132 p2(k)=dble(pptl(k,i))
134 p2(4)=sqrt(p2(3)**2+p2(2)**2+p2(1)**2+p2(5)**2)
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
144 write(*,*)'Precision problem in gakfra, p:',
145 & p1(5),p1(4),p1(3),p1(2),p1(1)
146 write(*,*)'origin : ',iorrem
148 c if(iorrem.ne.0.and.
149 c & (ityptl(iorrem).eq.40.or.ityptl(iorrem).eq.10))
150 c & p1(5)=dble(pptl(5,iorrem))
152 write(*,*)'string mass : ',ityptl(iorrem),pptl(5,iorrem)
153 p1(5)=dble(pptl(5,iorrem))
160 if(iappl.ne.6)qsqptl(nk1)=p1(5)*p1(5) !hadronic string energy
161 !(already defined for lepton)
163 c.....light like ...............................................
165 if(ish.ge.6) call blist("before light like&",nk1,nk2)
166 if(rangen().lt.0.5)then
176 p2(i)=dble(pptl(i,i1))+dble(pptl(i,i2))
178 c p2(5)=sqrt(p2(4)**2-p2(3)**2-p2(2)**2-p2(1)**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
189 if(ish.ge.6)write(ifch,'(a,2i4,9f12.6)') 'overlaps'
190 $ ,i1,i2,ak1,ak2,sqrt(2.*am12)
193 a1=(1.+ak1)*pptl(j,i1)-ak2*pptl(j,i2)
194 a2=(1.+ak2)*pptl(j,i2)-ak1*pptl(j,i1)
200 if(nk2-nk1.eq.1) ii=-1
210 if(ish.ge.6) call blist("after light like&",nk1,nk2)
212 c.....on-shell and copy ...............................................
214 if(ish.ge.6) write (ifch,*) 'on-shell check'
217 p2(k)=dble(pptl(k,i))
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)
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)
230 if(i.ne.nk1.and.i.ne.nk2)f=0.5
233 aemax=1000. ! here max band mass
234 if(f*p2(4).gt.aemax) then
235 nmax = int(f*p2(4)/aemax)+1
237 c print *,"nmax",nmax
242 if(i.eq.nk1.and.nn.eq. 1 ) f=1.
243 if(i.eq.nk2.and.nn.eq.nmax) f=1.
245 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
246 c if(nmax.ge.2) print *, nob,ff,f,i
248 pst(k,nob)=ff*f*p2(k)
255 do i1=nob-1,1,-1 ! copy again gluons
257 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
263 nob=nob+1 ! nob is number of kinks
264 if(nob.gt.mxnba-7)stop'gakfra: increase parameter mxnba '
266 write (ifch,*) 'bands:'
268 write (ifch,'(i4,4g20.13)') i,(pst(k,i),k=1,4)
272 c.....total momentum...............................................
278 pst(k,nob)=pst(k,nob)+pst(k,i)
283 c.....pbreak ................................
297 stop'in gakfra: unknown kky'
306 pb=0.14369+ 854.965/x**3 - 111.1144/x**2 + 7.163/x
307 c !! {{14,0.4},{22,0.32},{34,0.28},{91.2,0.21}}
308 c pb=0.143+ 900/x**3 - 77/x**2 + 4.5/x
309 c pb=0.165+ 1020/x**3 - 80/x**2 + 3.8/x
310 c pb=0.17+ 950/x**3 - 76/x**2 + 3.8/x
311 c pb=0.165+ 600/x**3 - 50/x**2 + 3.8/x
312 c pb=0.135+ 855./x**3 - 111./x**2 + 7./x
313 c pb=0.05+ 1100./x**3 - 80./x**2 + 5./x
317 c write(*,*)'pb',kky,krm,kdrm,x,pb,zzrm,zzzz
320 c.....write total string on pptl - list.............................
322 if(nptl.gt.mxptl)call utstop('gakfra: mxptl too small ...&')
323 call idtr5(idptl(nk1),ic1)
324 call idtr5(idptl(nk2),ic2)
327 idptl(nptl) = 8*10**8+ ic(1)*100 + ic(2)/100 !nbr+1
334 write(ifch,'(a)')' ---------total string------------'
335 call blist('&',nptl,nptl)
338 c....................................................................
349 c calculate isospin to conserve isospin symmetry in popcorn
350 c 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))
356 iso=isospin(abs(iflb(0)))*sign(1,-iflb(0))
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))
363 iso=iso+isospin(abs(iflb(1)))*sign(1,iflb(1))
366 c...................light string....................................
367 amm=ammin(idptl(nk1),idptl(nk2))
369 if(sqrt(max(0.,pf(nob,nob))).lt.amm)then
370 id=idsp( idptl(nk1),idptl(nk2) )
373 . '-------string too light to decay--------'
374 write (ifch,*) id,p1(5),amm
375 write (ifch,*) id, sqrt(max(0.,pf(nob,nob))) ,amm
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
382 & call utstop('gakfra (light): mxptl too small ...&')
391 istptl(n8ptl)=29 !code for fragmented string
392 ityptl(n8ptl)=ityptl(nk1)
393 itsptl(n8ptl)=itsptl(nk1)
397 ifrptl(1,n8ptl)=n8ptl+1
399 ityptl(nptl)=ityptl(nk1)
400 itsptl(nptl)=itsptl(nk1)
413 xorptl(j,nptl)=xorptl(j,nk1)
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
419 c write (*,*) 'light string --- particle off-shell!!!!'
420 c write (*,*) idr,' mass:',p1(5),' should be:',am
425 c................search breakpoints...................................
431 !----------------------!
432 11 call gaksbp(0,1,iok)
433 !----------------------!
437 c..............delete breakpoint....................................
438 9 if(ish.ge.4)write(ifch,*) 'delete breakpoint',n,' -> '
439 &,nbr-1,' breakpoints left'
440 call gakdel(n) !hier loeschen
442 c..............no more breakpoints -> redo..........................
445 if(ish.ge.3)write (ifch,*)
446 & 'no breakpoints left ... try again (',nsbp,')'
447 if(ish.ge.3)write (ifch,*)' '
451 c................make index list of breakpoints to adjust...........
457 if(ip(i).eq.0.or.ip(i+1).eq.0)then
462 if(nbrtry.gt.1000000)then
464 write(*,*)'Gakfra : string can not be fragmented, redo event !'
465 write(*,*)nk1,nk2,nbr,nind,pb
471 c.....no more breakpoint to adjust...............................
472 if(nind.eq.0) goto 100
474 c................................................................
476 write(ifch,*) 'breakpoints:'
477 write(ifch,'(i3,i5)') 0,-iflb(0)
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)
483 write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
486 c.....choose breakpoint, calculate masses, lambda..................
488 nn=1+int(r*real(nind))
489 c nn=1+(nind-1)*int(r+0.5)
491 if(ish.ge.4)write(ifch,*) 'choose breakpoint',n
494 do while (nl.gt.0.and.ip(nl).eq.0)
497 do while (nr.lt.nbr+1.and.ip(nr+1).eq.0)
500 if(ish.ge.6)write (ifch,*) 'nl,n,nr:',nl,n,nr
501 c 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
504 cc x=(xyb(1,n)-x1)/(x2-x1)
505 cc y=(xyb(2,n)-y1)/(y2-y1)
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)
512 c.....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) )
518 c.....if mass to small (because of spacelike pt) reject..........
520 c if (aml2.le.amin.or.amr2.le.amin) goto 9
522 c.....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'
533 c.....if no adjusted mass on left side, check for resonance...........
535 aml=sqrt(max(0.,aml2+0.*min(0.,amr2))) !!!????
537 call idress(idl,aml,idrl,iadjl)
539 if(idrl.eq.110.and.r.lt.0.5)then
540 if (r.gt.eta+etap) goto 9 !rare numerical errors
543 if(r.lt.etap)aml=.958
545 call idress(idl,aml,idrl,iadjl)
548 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
549 & ' l: ',idl,amlxx,aml,idrl,iadjl
551 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
552 & ' l: ',idl,aml,aml,ip(n),'ok'
555 c.....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))) !!!????
559 call idress(idr,amr,idrr,iadjr)
561 if(idrr.eq.110.and.r.lt.0.5)then
562 if (r.gt.eta+etap) goto 9 !rare numerical errors
565 if(r.lt.etap)amr=.958
567 call idress(idr,amr,idrr,iadjr)
570 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,1i10,i5)')
571 & ' r: ',idr,amrxx,amr,idrr,iadjr
573 if(ish.ge.5)write(ifch,'(a,i5,2f12.6,i10,a5)')
574 & ' r: ',idr,amr,amr,ip(n+1),'ok'
577 c.....mass adjustments.........................................
579 if(ip(n+1).ne.0)then !.........adjusted mass on right side
581 call gaksbp(n-1,n,iok)
586 if(ish.ge.5)write(ifch,*)'mass adjustment 1'
588 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
592 elseif(ip(n).ne.0)then !.........adjusted mass on left side
594 call gaksbp(n,n+1,iok)
599 if(ish.ge.5)write(ifch,*)'mass adjustment 2'
601 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
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)
609 call gaksbp(n-1,n,iok)
612 elseif(idrl.ne.0.and.idrr.eq.0)then
614 if(ish.ge.5)write(ifch,*)'mass adjustment 3'
615 call gakanp(n-1,n,nr,aml**2,0.,ala2,iok)
617 elseif(idrl.eq.0.and.idrr.ne.0)then
619 if(ish.ge.5)write(ifch,*)'mass adjustment 4'
621 call gakanp(nl,n,n2,0.,amr**2,ala2,iok)
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'
627 call gakanp(n-1,n,n2,aml**2,amr**2,0.,iok)
635 write(ifch,*) 'left/right particles:'
636 & ,ip(n),aml,' ',ip(n+1),amr
640 c................................................................
641 c end of string decay
642 c................................................................
644 c.....final list...............................................
646 write(ifch,*) ' ************ OK **************'
647 write(ifch,*) 'final breakpoints:'
648 write(ifch,'(i3,i5)') 0,-iflb(0)
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)
654 write(ifch,'(i3,i5)') nbr+1,iflb(nbr+1)
657 c.....write particles in pptl-list................................
659 write(ifch,'(a)')' ---------produced particles---------'
661 c to avoid leading strange baryon for diffractive string, rotate string
664 if(kdrm.eq.1.and.krm.ne.0.and.nbr.gt.0)then
667 if(iq1.gt.3.or.iq2.gt.2)then
671 elseif(iso.gt.0)then !for u or d, use proper mass if isospin > 0 (proton)
673 elseif(iso.lt.0)then !us d mass for u and u mass for u if isospin < 0
676 am=0.5*(qmass(1)+qmass(2))
678 zz=zrminc*pf(nob,nob)**2.
679 zz=zz*(1.+zrmmax*log(float(matarg+maproj-1)))
680 c if(itp/10.eq.4)then
681 c zz=zz*(1.+zrmmax*log(max(1.,float(matarg)/float(maproj))))
682 c elseif(itp/10.eq.5)then
683 c zz=zz*(1.+zrmmax*log(max(1.,float(maproj)/float(matarg))))
686 c print *,iso,itp,pf(nob,nob),am,zz,iq1,iq2,exp(-min(100.,zz))
688 c if(kdrm.eq.1)print *,'reminv',itp,am,pf(nob,nob),zz,iq1,iq2
689 c & ,reminv+exp(-min(100.,zz))
690 if(rangen().gt.reminv+exp(-min(100.,zz)))then
697 if(nptlini+nbr+1.gt.mxptl)
698 & call utstop('gakfra (end): mxptl too small ...&')
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)
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)
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 ...'
722 pptl(j,nptl)=pst(j,nob+1)-
723 & x*pst(j,nob+2)+y*pst(j,nob+3)
724 c pptr(j,nptl)=ptb(j,i)
726 else !last particle = right side of last breakpoint
730 pptl(j,nptl)=pst(j,nob+4)+
731 & x*pst(j,nob+5)-y*pst(j,nob+6)
732 c pptr(j,nptl)=ptb(j,i) !should be zero
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)
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)
755 pptl(4,nptlini+1)=sqrt(pptl(4,nptlini+1))
756 pptl(4,nptlini+2)=sqrt(pptl(4,nptlini+2))
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))
764 c call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5)
765 c $ ,pptr(1,nptl),pptr(2,nptl),pptr(3,nptl),pptr(4,nptl))
767 c........Origin..................................................
775 tauran=-taurea*alog(r)*pptl(4,nptl)/pptl(5,nptl)
776 c if(jpsi.ge.1)tauran=0.
777 c tauran=max(taubrl,taubrr) !take formation time from string-theory
779 xorptl(j,nptl)=xorptl(j,nk1)+pptl(j,nptl)
780 & /pptl(4,nptl)*tauran
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)
788 if(ish.ge.3)call blist('&',nptl,nptl)
789 ctp060829 taubrl=taubrr
793 ityptl(n8ptl)=ityptl(nk1)
799 istptl(n8ptl)=29 !code for fragmented string
800 ifrptl(1,n8ptl)=n8ptl+1
803 write(ifch,*)'string momentum sum:'
807 pptl(kk,nptl+1)=pptl(kk,nptl+1)+pptl(kk,ii)
810 call alist2('&',n8ptl,n8ptl,nptl+1,nptl+1)
815 c.....another string?.........................................
819 c.....end of fragmentation.....................................
821 call utprix('gakfra',ish,ishini,4)
825 c---------------------------------------------------------------------
826 subroutine gaksbp(ibr1,ibr2,iok)
827 c---------------------------------------------------------------------
828 ! search break points
829 !-----------------------------------------
830 ! nbr ... number of break points
834 !--------------------------------------------------------------
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
845 common/czz/kky,krm,zzod,zzop,kdrm,zzoq,zzos,zzrm,zzipp,zzipt
846 &,zzid,zzip,zzzzr,itp
849 dimension pleft(5),pright(5)
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)
855 call utpri('gaksbp',ish,ishini,5)
861 if(rangen().lt.0.5)iside=2
862 if(ish.ge.6)write(ifch,*)'iside:',iside
864 if(ish.ge.4)write(ifch,*)
865 &'search brk between ib1=',ib1,' and ib2=',ib2
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
875 if(iside.eq.2)goto 51
876 c.....................................................................
877 if(ib1.eq.0)then !startwert der j-schleife
883 if(jj.eq.ijb(2,ib1) )then !linker brk im Streifen jj?
888 if(jj.eq.ijb(2,ib2))then !rechter brk im Streifen jj?
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
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
908 $ mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
914 if(ish.ge.6)write(ifch,*) 'Rand erreicht'
918 if(ii.eq.ijb(1,ib1))then !linker brk im Feld (ii,jj)
924 if(ii.eq.ijb(1,ib2))then !rechter brk im Feld (ii,jj)
929 if(x1.eq.x2) goto 9999
930 if(abs(x1-x2)/abs(x1+x2).le.1e-7) goto 9999
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
936 if(ii.eq.ijb(1,ib2))then
940 if(ii.eq.jj.or.mmod(ii+jj).eq.0)weiter=.false.
946 if(jj.eq.ijb(2,ib2)) then
947 if(ish.ge.6)write(ifch,*) 'brk erreicht'
956 if(ish.ge.6)write(ifch,*)'jj,yb ok:',jj,yb
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
963 if(jj.lt.nob/2)ii=mmod(nob-jj+1)
967 $ mmod(nob-jj+1+nob/2).gt. mmod(ii+nob/2)
971 if(ii.eq.ijb(1,ib1))then !linker brk im Feld (ii,jj)
977 if(ii.eq.ijb(1,ib2))then !rechter brk im Feld (ii,jj)
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
994 if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
995 & ,' at position ',xb,yb
998 c......................................................................
1000 51 if(ib2.eq.nbr+1)then !startwert der i-schleife
1005 do while (ii.ne.nob/2)
1006 if(ii.eq.ijb(1,ib1))then !linker brk im Streifen (ii)
1011 if(ii.eq.ijb(1,ib2))then !rechter brk im Streifen (ii)
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
1020 if(ib2.eq.nbr+1)then
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
1029 if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1034 if(ish.ge.6)write(ifch,*) 'Rand erreicht'
1038 if(jj.eq.ijb(2,ib1))then !linker brk im Feld (ii,jj)
1043 if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
1048 if(y1.eq.y2) goto 9999
1049 if(abs(y1-y2)/abs(y1+y2).le.1e-7) goto 9999
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
1055 if(jj.eq.ijb(2,ib1))then
1059 if(jj.eq.ii.or.mmod(ii+jj).eq.0)weiter=.false.
1065 if(ii.eq.ijb(1,ib1)) then
1066 if(ish.ge.6)write(ifch,*) 'brk erreicht'
1077 if(ish.ge.6)write(ifch,*)'ii,xb ok:',ii,xb
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
1084 if(ii.gt.nob/2)jj=mmod(nob-ii+1)
1087 if(ii.gt.nob/2 .and. mmod(nob-ii+1).gt.jj)jj=mmod(nob-ii+1)
1091 if(jj.eq.ijb(2,ib1))then !linker brk im Feld (ii,jj)
1097 if(jj.eq.ijb(2,ib2))then !rechter brk im Feld (ii,jj)
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
1114 if(ish.ge.5)write(ifch,*) 'breakpoint in field ',ii,jj
1115 & ,' at position ',xb,yb
1119 c.....breakpoint accepted......................
1121 if(nbr.gt.mxnbr-2) stop 'gaksbp: increase parameter mxnbr'
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))
1148 amlast=f*qmass(abs(iflb(ib1)))
1152 elseif(itp.ge.40)then
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))
1164 xlastdiq=xlastdiq+f*qmass(abs(iflb(i)))
1167 if(pf(nob,nob).le.5.)amlast=1.e-10 !if very low energy, strangeness should not be totally killed
1168 c if(pf(nob,nob).le.5.)then
1169 c xlastdiq=-(ampt+2.*qmass(3)) !if very low energy, strangeness should not be totally killed
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)
1177 c ..........diquark...............................................
1186 elseif(kky.eq.1)then
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))))
1196 c if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.diqcut)pdiqu=0.
1198 if(rangen().lt.pdiqu.and.iabs(iflb(ib1)).lt.6
1199 & .and.iabs(iflb(ib2)).lt.6)then
1204 c ..........flavor...............................................
1208 pud1= ( 1 - min(0.33,zzos*(1-2*pudr)) ) /2.
1214 elseif(kky.eq.1)then
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)))))
1225 c if(kdrm.eq.1.and.sqrt(pf(nob,nob)).lt.strcut)pud1=0.5
1227 pud2=(1.-(1.-2.*pud1)*puds)*0.5 !from e+e- data
1229 if(iabs(iflb(ib1)).gt.1000)then
1230 ifl1=int(rangen()**0.625/pud1)+1
1232 ifl1=int(rangen()**1.15/pud1)+1 !assymmetric u and d relevant because of assymmetry Kch / K0 in e+e-
1235 c amk0=2.*qmass(ifl1) !mass for mt distribution
1238 if(jqu.eq.2)then !diquark ------
1239 ifl2=int(rangen()**0.625/pud2)+1
1240 c amk0=amk0+2.*qmass(ifl2)+qmass(0) !mass for mt distribution with bounding energy for diquark (qmass(0))
1242 ifl=-min(ifl1,ifl2)*1000-max(ifl1,ifl2)*100
1246 if(iflb(ib1).lt.0.and.iflb(ib1).ge.-6)ifl=-ifl
1247 if(iflb(ib1).gt.1000)ifl=-ifl
1249 if(ish.ge.5)write(ifch,*) 'flavor,pud1/2:',ifl,pud1,pud2,pdiqu
1250 c..............................pt.......................................
1251 if(krm.ne.1.and.kdrm.ne.1)then
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)
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 -----
1272 pleft(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
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-----
1279 pright(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
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
1287 yleft=sign(1.,pleft(3))*alog((pleft(4)+abs(pleft(3)))/amt)
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
1294 yright=sign(1.,pright(3))*alog((pright(4)+abs(pright(3)))/amt)
1298 ybrk=(yleft+yright)/2.
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)
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)
1313 if(ifl1.eq.3.and.ifl2.eq.0)delptfra=delptfra+ptfrasr
1314 c if(ifl1.eq.3.or.ifl2.eq.3)delptfra=delptfra+ptfrasr
1316 if(nsbp.gt.9)fnsbp=0
1318 c for a better pt distribution, we generate mt=sqrt(pt2+m02) (with m0=amk0) instead of pt (do not give good results !)
1321 c ptadd=ranpt()*(ptfrair+delptfra)*(zzop-1)*fnsbp !*jqu
1322 c ptadd=ranpt()*ptfrair*(zzop-1)*fnsbp !*jqu
1323 ptadd=(ptfrair+delptfra)*(zzop-1)*fnsbp !*jqu
1325 c pt=ranpt()*(ptfrair +delptfra ) + ptadd
1326 c pt=ranpt()*ptfrair
1327 pt=ranpt()*(ptfrair + ptadd )
1329 pt=ranpt()*(ptfrair + ptadd )
1330 c pt=ranpt()*ptfraqq
1331 c pt=ranpt()*(ptfraqq + ptadd )
1337 c pt=ranpt()*(ptfradr +delptfra )
1341 c pt=ranpt()*(ptfradr )
1344 elseif(kky.eq.1)then
1345 c ptadd=ranptk()*(ptfrak+delptfra)*(zzipx-1)*fnsbp !*jqu
1346 pttra=(ptfrak+delptfra)*(zzipx-1)*fnsbp !/jqu
1348 !pt=ranptk()*(ptfrak+delptfra + pttra )
1349 !pt=ranptk()*(ptfrak + pttra)
1350 pt=ranptk()*ptfrak + ranptk()*pttra
1352 !pt=ranptk()*(ptfrakd + pttra)
1353 pt=ranptk()*ptfrakd + ranptk()*pttra
1356 c ptadd=ranpt()*(ptfra+delptfra)*(zzip-1)*fnsbp !*jqu
1357 c pttra=ranpt()*ptfra*(zzip-1)*fnsbp !*jqu
1358 pttra=(ptfra+delptfra)*(zzip-1)*fnsbp !*jqu
1360 c pt=ranpt()*(ptfra+delptfra ) + pttra
1362 pt=ranpt()*(ptfra + pttra )
1364 c pt=ranpt()*(ptfra+delptfra) + pttra
1365 c pt=ranpt()*ptfrakd
1366 pt=ranpt()*(ptfra + pttra )
1369 c pt=sqrt(pt*pt+2.*pt*amk0)+ptadd !sample mt-m0 instead of pt ... doesn't work at RHIC
1373 write(ifch,*) 'pt:',pt
1375 ptb(1,ib1+1)=pt*cos(beta)
1376 ptb(2,ib1+1)=pt*sin(beta)
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)
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)) !?????????
1392 write(ifch,*) 'boost vector ( region ) '
1393 write(ifch,'(5g12.6)') ax,ay,az,ae,am,pf(ii,jj)
1399 call utlob2(1,ax,ay,az,ae,am,bx,by,bz,be,60)
1401 write (ifch,*) 'boost of b'
1402 write (ifch,*) 'bx,by,bz,be',bx,by,bz,be
1404 call utrot4(-1,bx,by,bz,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
1406 write (ifch,*) 'rot of pt'
1407 write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
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))
1414 write (ifch,*) 'backboost of pt'
1415 write(ifch,'(4f12.8)') (ptb(i,ib1+1),i=1,4)
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
1420 c.....call utrota(-1,ax,ay,az,ptb(1,ib1+1),ptb(2,ib1+1),ptb(3,ib1+1))
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)
1436 c if(Amax.lt.Amxx) goto 15
1440 if(ish.ge.6)write(ifch,*) 'Amax:',Amax,Amaxn
1441 if (nbr.eq.nbrold) then
1444 if(pb*Amax.le.0..or.ntry.ge.1000)then
1451 write(ifch,*) 0,iflb(0)
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,*) '.................'
1459 write(ifch,*) nbr+1,iflb(nbr+1)
1462 9999 if(nbr.eq.nbrold)iok=1
1463 call utprix('gaksbp',ish,ishini,5)
1467 c----------------------------------------------------------------------
1468 function ranptcut(xcut)
1469 c----------------------------------------------------------------------
1470 c .........exp(-x**2)
1471 12 x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1472 C 12 x=sqrt(-(log(rangen()))/(3.1415927/4.)/xcut) !gauss
1475 c if(xcut.gt.0.)then
1476 c if(rangen().lt.x/xcut)goto 12
1486 c do while (r.gt.1.)
1487 c 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1489 c r=rangen() / ( exp(-x)*(1+x**2) )
1495 cc----------------------------------------------------------------------
1496 c function ranpticut(xcut)
1497 cc----------------------------------------------------------------------
1503 c do while (r.gt.1.)
1504 c 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1506 c r=rangen() / ( exp(-x)*(1+x**2) )
1509 c if(rangen().gt.xcut/x)goto 12
1515 c----------------------------------------------------------------------
1517 c----------------------------------------------------------------------
1523 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1525 r=rangen() / ( exp(-x)*(1+x**2) )
1533 c .........exp(-x**2)
1534 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1536 c .........exp(-sqrt(x))
1540 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1541 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1547 c----------------------------------------------------------------------
1549 c----------------------------------------------------------------------
1555 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1557 r=rangen() / ( exp(-x)*(1+x**2) )
1566 c .........exp(-x**2)
1567 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1569 c .........exp(-sqrt(x))
1573 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1574 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1580 c----------------------------------------------------------------------
1582 c----------------------------------------------------------------------
1584 c .........exp(-x**2)
1585 ranptg=sqrt(-log(rangen())/(3.1415927/4.)) !gauss
1597 11 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1599 r=rangen() / ( exp(-x)*(1+x**2) )
1603 c .........exp(-sqrt(x))
1607 x=sqrt(exp(rangen()*log(1+xmx**2))-1)
1608 r=rangen() / ( exp(-sqrt(x))*(1+x**2)/5. )
1617 c----------------------------------------------------------------------
1618 subroutine gakdel(ibr)
1619 c----------------------------------------------------------------------
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
1641 c----------------------------------------------------------------------
1642 subroutine gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1643 c----------------------------------------------------------------------
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
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)
1657 call utpri('gaksco',ish,ishini,8)
1660 write(ifch,*) 'zwischen brk:',ibr1,ibr,ibr2,'(',nob,')',nbr
1661 write(ifch,*) 'region:',ib,jb
1663 if(ib.eq.ijb(1,ibr1))then
1668 if(ib.eq.ijb(1,ibr2))then
1673 if(jb.eq.ijb(2,ibr1))then
1678 if(jb.eq.ijb(2,ibr2))then
1684 c.....left side...................................................
1686 if(ish.ge.8)write(ifch,*)'x1,x2',x1,x2
1688 cc 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)
1691 if(ish.ge.8)write(ifch,*) 'add a 1 ii',1.,ib
1694 if(ib.eq.ijb(1,ibr1))weiter=.false.
1695 do while(ii.ne.jb.and.weiter) !linker Rand??
1697 if(ii.eq.ijb(1,ibr1))f1=xyb(1,ibr1)
1699 pst(i,n)=pst(i,n)+f1*pst(i,ii)
1701 if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1702 if(ii.eq.ijb(1,ibr1))weiter=.false.
1707 if(jb.eq.ijb(2,ibr1))weiter=.false.
1710 if(jj.eq.ijb(2,ibr1))f1=1.-xyb(2,ibr1)
1712 pst(i,n)=pst(i,n)+f1*pst(i,jj)
1714 if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1715 if(jj.eq.ijb(2,ibr1))weiter=.false.
1719 cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1720 pst(i,n+1)= pst(i,ib)
1721 cc pst(i,n+2)=(y2-y1)*pst(i,jb)
1722 pst(i,n+2)= pst(i,jb)
1727 co(4)=-2.*pf(n+1,n+2)
1730 write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1733 if(ish.ge.8)write(ifch,*) 'co left:',co(1),co(2),co(3),co(4)
1735 c.....right side...................................................
1737 if(ish.ge.8)write(ifch,*)'y1,y2',y1,y2
1739 cc 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)
1742 if(ish.ge.8)write(ifch,*) 'add a 1 jj',1.,jb
1745 if(ib.eq.ijb(1,ibr2))weiter=.false.
1746 do while(ii.ne.jb.and.weiter)
1748 if(ii.eq.ijb(1,ibr2))f1=1.-xyb(1,ibr2)
1750 pst(i,n)=pst(i,n)+f1*pst(i,ii)
1752 if(ish.ge.8)write(ifch,*) 'add a f1 ii',f1,ii
1753 if(ii.eq.ijb(1,ibr2))weiter=.false.
1758 if(jb.eq.ijb(2,ibr2))weiter=.false.
1761 if(jj.eq.ijb(2,ibr2))f1=xyb(2,ibr2)
1763 pst(i,n)=pst(i,n)+f1*pst(i,jj)
1765 if(ish.ge.8)write(ifch,*) 'add b f1 ii',f1,jj
1766 if(jj.eq.ijb(2,ibr2))weiter=.false.
1770 cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1771 pst(i,n+1)= pst(i,ib)
1772 cc pst(i,n+2)=(y2-y1)*pst(i,jb)
1773 pst(i,n+2)= pst(i,jb)
1778 co(8)=-2.*pf(n+1,n+2)
1781 write (ifch,'(4g12.5)') (pst(j,i),j=1,4)
1784 if(ish.ge.8)write(ifch,*) 'co right:',co(5),co(6),co(7),co(8)
1786 c.....lambda (absolute past).....................................
1789 cc pst(i,n)= x1*pst(i,ib)+y1*pst(i,jb)
1793 do while (mmod(ii+jb).ne.0)
1794 if(ish.ge.8)write(ifch,*) 'add lambda',ii
1796 pst(i,n)=pst(i,n)+pst(i,ii)
1801 cc pst(i,n+1)=(x2-x1)*pst(i,ib)
1802 cc 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)
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)
1814 c---------------------------------------------------------------------
1815 subroutine gakanp(ibr1,ibr,ibrr2,aml2,amr2,ala2,iok)
1816 c---------------------------------------------------------------------
1817 c mass adjustment of fragments
1818 c ibr1-ibr ibr-ibrr2
1819 c where ibr1,ibr,ibrr2 are break point indices
1820 c aml2,amr2 are the reqired squared masses (if zero -> any mass)
1821 c iok=0 (ok) or ok=1 (error)
1822 c---------------------------------------------------------------------
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)
1831 c pf(i,j)=pst(4,i)*pst(4,j)-pst(3,i)*pst(3,j)
1832 c & -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)
1837 if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1843 if((nin(1,i).eq.ib.and.nin(2,i).eq.jb)
1844 $ .or.(ipst(ib).eq.ipst(jb)))then
1847 write(ifch,*) 'error ... endless loop'
1848 if(ib.eq.ipst(jb)) write(ifch,*) ' in zero mass region'
1854 if(ni.gt.mxnin)stop'gakanp: increase parameter mxnin '
1857 if(ish.ge.6)write(ifch,*)
1858 if(ish.ge.6)write(ifch,*) 'ib,jb=',ib,jb
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'
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)
1907 call gaksco(ibr1,ibr,ibr2,ib,jb,x1,x2,y1,y2)
1908 c if(ni.eq.1)print *,'------------------------------------',2
1909 cc x=(xyb(1,ibr)-x1)/(x2-x1)
1910 cc y=(xyb(2,ibr)-y1)/(y2-y1)
1916 if(ish.ge.6)write(ifch,*) ibr1,ibr,ibr2,aml2,amr2,ala2,iok
1920 elseif(aml2.le.0.)then
1923 elseif(ala2.le.0.)then
1927 stop' not like this , please...'
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
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))
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
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
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))
1960 if(ish.ge.5)write(ifch,*)'error ... no solution of quadr equ'
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
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))
1973 if(ish.ge.5)write(ifch,*)'error ... no solution (2) '
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 '
1982 50 if(x.gt.x1.and.x.lt.x2.and.y.gt.y1.and.y.lt.y2)then
1983 cc xyb(1,ibr)=x1+(x2-x1)*x
1984 cc xyb(2,ibr)=y1+(y2-y1)*y
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'
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
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))
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))))
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
2017 write(ifch,*) 'away:'
2018 & ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
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))))
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
2036 write(ifch,*) 'away:'
2037 & ,i,xyb(1,i),xyb(2,i),ijb(1,i),ijb(2,i)
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
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
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
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
2080 if(ish.ge.5)write(ifch,*)'error ... x,y not in allowed range'
2083 9999 if(amr2.eq.0.) ibrr2=ibr2
2084 call utprix('gakanp',ish,ishini,6)
2087 cc----------------------------------------------------------------------
2088 c subroutine gakstr(ifl)
2089 cc----------------------------------------------------------------------
2091 cc calculates string-fragments by taking off pt of breakup
2092 cc do with ifl=1 undo with ifl=-1
2094 cc----------------------------------------------------------------------
2095 c include 'epos.inc'
2096 c common /cpptr/ pptr(4,mxptl),ystr(mxptl)
2099 c if(istptl(i).eq.29)then
2103 c if ((istptl(j).eq.0.or.istptl(j-1).eq.0).and.j.ne.nk1) then
2105 c pptl(k,j)=pptl(k,j)+pptr(k,j-1)*ifl
2107 c !write(ifch,*)"left side back to ",j,(pptr(k,j-1),k=1,4)
2109 c if ((istptl(j).eq.0.or.istptl(j+1).eq.0).and.j.ne.nk2) then
2111 c pptl(k,j)=pptl(k,j)-pptr(k,j)*ifl
2113 c !write(ifch,*)"right side back to ",j,(-pptr(k,j),k=1,4)
2115 c if(ifl.eq.-1.and.istptl(j).eq.0)then
2116 c e=pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2
2119 c !dif=abs(e-pptl(4,j))
2120 c !if(dif.gt.0.01.and.dif/e.gt.0.1)print*,j,e,pptl(4,j)
2125 c if(istptl(i).eq.0)then
2126 c if ( ifl.eq.1 ) then
2127 c ystr(i)=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))
2128 c * /sqrt(pptl(5,i)**2+pptl(1,i)**2+pptl(2,i)**2) )
2134 c----------------------------------------------------------------------
2135 subroutine gakli2(nn1,nn2)
2136 c----------------------------------------------------------------------
2139 double precision pgampr,rgampr
2140 common/cgampr/pgampr(5),rgampr(4)
2141 c double precision db1,db2,db3,db4,db5
2142 character label*8,idlabl*8
2146 c db4=sqrt(rgampr(1)**2+rgampr(2)**2+rgampr(3)**2)
2147 c db5=sqrt( db4**2-rgampr(4)**2)
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'
2156 if (idptl(i).lt.10000)then
2158 label=idlabl(idptl(i))
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)
2167 & ,istptl(i),ityptl(i),label
2169 125 format (1i4,5g18.10,4i6,1i10
2174 126 format (A4,7g12.4,i5)
2178 c----------------------------------------------------------------------
2180 cc----------------------------------------------------------------------
2182 c include 'epos.inc'
2183 c parameter (mxnbr=500,mxnba=5000)
2184 c common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2185 c $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2186 c & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2189 c write (ifch,10) i,(pst(j,i),j=1,4)
2191 c 10 format(1i4,5g18.10)
2195 cc----------------------------------------------------------------------
2197 cc----------------------------------------------------------------------
2199 c include 'epos.inc'
2200 c parameter (mxnbr=500,mxnba=5000)
2201 c common /gag/nob,pst(4,0:mxnba),ipst(0:mxnba)
2202 c $ ,nbr,ijb(2,0:mxnbr),xyb(2,0:mxnbr)
2203 c & ,ptb(4,0:mxnbr),iflb(0:mxnbr),ip(0:mxnbr),co,x,y
2204 c dimension co(12),p1(5),p2(5)
2206 c write(ifch,*) 'particle list of string decay'
2211 c if(i.lt.nbr+1)then
2212 c call gaksco(i-1,i,i+1,ijb(1,i),ijb(2,i),x1,x2,y1,y2)
2214 ccc x=(xyb(1,i)-x1)/(x2-x1)
2220 ccc y=(xyb(2,i)-y1)/(y2-y1)
2225 c aml2=co(1)+x*co(2)+y*co(3)+x*y*co(4)
2226 c amr2=co(5)+x*co(6)+y*co(7)+x*y*co(8)
2227 c aml=sign(sqrt(abs(aml2)),aml2)
2228 c amr=sign(sqrt(abs(amr2)),amr2)
2230 c p2(j)=pst(j,nob+1)-x*pst(j,nob+2)+y*pst(j,nob+3)
2236 c p2(j)=pst(j,nob+4)+x*pst(j,nob+5)-y*pst(j,nob+6)
2241 c write(ifch,'(2i4,i6,a,i5,i10,5g14.6)') i-1,i
2242 c & ,-iflb(i-1),'==',iflb(i),ip(i)
2245 c am2=p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2
2246 c p1(5)=sign(sqrt(abs(am2)),am2)
2247 c write(ifch,'(12x,a60)')
2248 c & '------------------------------------------------------------'
2249 c write(ifch,'(14x,5g14.6)') (p1(j),j=1,5)
2255 c---------------------------------------------------------------------
2256 subroutine idress(id,am,idr,iadj)
2257 c---------------------------------------------------------------------
2259 call idres(id,am,idr,iadj,1)
2263 ids=max(mod(iabs(id)/100,10),mod(iabs(id)/10,10))
2264 if(iabs(idr).le.999) then
2265 c write (*,*) ' ',id,idr
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)
2272 idr=sign(iabs(id)+int(rangen()+0.75),id)
2274 c write (*,*) '->',id,idr
2276 elseif(iabs(idr).le.9999)then
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)
2281 idr=sign(iabs(id)+1,id)
2293 c---------------------------------------------------------------------
2294 SUBROUTINE gaksphe(sphe,r,mstu41)
2296 C...Purpose: to perform sphericity tensor analysis to give sphericity,
2297 C...aplanarity and the related event axes. stolen from jetset ;-)
2300 DIMENSION SM(3,3),SV(3,3)
2302 C...Calculate matrix to be diagonalized.
2312 IF(istptl(i).ne.0) GOTO 140
2313 IF(MSTU41.GE.2) THEN
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
2322 PA=SQRT(pptl(1,i)**2+pptl(2,I)**2+pptl(3,i)**2)
2324 IF(ABS(r-2.).GT.0.001) PWT=MAX(1E-10,PA)**(r-2.)
2327 SM(J1,J2)=SM(J1,J2)+PWT*pptl(j1,i)*pptl(j2,i)
2333 C...Very low multiplicities (0 or 1) not considered.
2337 write(ifch,*) 'too few particles for analysis'
2345 SM(J1,J2)=SM(J1,J2)/PS
2349 C...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
2361 CALL utmsg('gaksphe')
2362 write(ifch,*) 'all particles back-to-back'
2369 C...Find first and last eigenvector by solving equation system.
2372 SV(J1,J1)=SM(J1,J1)-sphe(4,I)
2381 IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 190
2390 RL=SV(J1,JB)/SV(JA,JB)
2392 SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
2393 IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 210
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)
2407 sphe(j,I)=SGN*sphe(j,I)/PA
2411 C...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))
2419 pptl(j,nptl+i)=sphe(j,i)
2424 C...Calculate sphericity and aplanarity. Select storing option.
2425 ccc SPH=1.5*(sphe(4,2)+sphe(4,3))
2426 ccc APL=1.5*sphe(4,3)
2431 C*********************************************************************
2433 SUBROUTINE gakthru(thru,mstu41)
2435 C...Purpose: to perform thrust analysis to give thrust, oblateness
2436 C...and the related event axes. stolen from jetset ;-)
2438 DIMENSION TDI(3),TPR(3)
2441 C...Take copy of particles that are to be considered in thrust analysis.
2450 IF(istptl(i).ne.0)goto 100
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
2458 IF(nptl+NP.GE.mxptl) THEN
2459 CALL utstop('gakthru: no more memory left in cptl&')
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)
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)
2475 C...Very low multiplicities (0 or 1) not considered.
2478 write(ifch,*) 'too few particles for analysis'
2484 C...Loop over thrust and major. T axis along z direction in latter case.
2487 c PHI=GAKANG(pptl(1,nptl+NP+1),pptl(2,nptl+NP+1))
2488 c CALL lurot(nptl+1,nptl+NP+1,0.,-PHI)
2489 c THE=GAKANG(pptl(3,nptl+NP+1),pptl(1,nptl+NP+1))
2490 c 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))
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.
2507 C...Find and order particles with highest p (pT for major).
2508 DO 110 ILF=nptl+NP+4,nptl+NP+MSTU44+4
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
2516 pptl(j,ILF+1)=pptl(j,ILF)
2521 pptl(j,ILF+1)=pptl(j,I)
2525 C...Find and order initial axes with highest thrust (major).
2526 DO 170 ILG=nptl+NP+MSTU44+5,nptl+NP+MSTU44+15
2529 NC=2**(MIN(MSTU44,NP)-1)
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
2538 TDI(J)=TDI(J)+SGN*pptl(j,nptl+NP+ILF+3)
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
2545 pptl(j,ILG+1)=pptl(j,ILG)
2548 ILG=nptl+NP+MSTU44+4
2550 pptl(j,ILG+1)=TDI(J)
2555 C... Iterate direction of axis until stable maximum.
2556 pptl(4,nptl+NP+ILD)=0.
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)
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))
2570 TPR(J)=TPR(J)+SGN*pptl(j,I)
2573 THP=SQRT(TPR(1)**2+TPR(2)**2+TPR(3)**2)/PS
2574 IF(THP.GE.THPS+PARU48) GOTO 270
2576 C... 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))
2579 IF(THP.GT.pptl(4,nptl+NP+ILD)+PARU48)
2582 SGN=(-1.)**INT(rangen()+0.5)
2584 pptl(j,nptl+NP+ILD)=SGN*TPR(J)/(PS*THP)
2586 pptl(4,nptl+NP+ILD)=THP
2587 pptl(5,nptl+NP+ILD)=0.
2590 IF(IAGR.LT.MSTU45.AND.ILG.LT.MIN(10,NC)) GOTO 260
2593 C... 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.
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))
2603 pptl(4,nptl+NP+3)=THP/PS
2604 pptl(5,nptl+NP+3)=0.
2607 C... 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))
2615 thru(j,ild)=pptl(j,nptl+NP+ild)
2619 C...Calculate thrust and oblateness. Select storing option.
2621 ccc OBL=thru(4,2)-thru(4,3)
2626 subroutine gakjet(ijadu)
2628 common/njjjj/njets(5,2,mxbins)
2639 ycut=xminim*(xmaxim/xminim)**((real(i)-0.5)/nrbins)
2640 c 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
2650 common/njjjj/njets(5,2,mxbins)
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
2659 C*********************************************************************
2660 function gaknjt(ycut,ijadu)
2662 c ijadu 1 = JADE 2=DURHAM
2663 c ycut - max. distance
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
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
2678 bet(i)=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(3,i)**2)
2680 ska(i,j)=pptl(1,i)*pptl(1,j)+pptl(2,i)*pptl(2,j)
2681 & +pptl(3,i)*pptl(3,j)
2683 a2c(i,j)= ((bet(i)*bet(j)-ska(i,j))
2684 & *2.*bet(i)*bet(j) / (0.00001+bet(i)+bet(j))**2 )
2689 if (istptl(i).eq.0) then
2692 pptl(j,nptl+nn)=pptl(j,i)
2699 do while (iflag.eq.0.and.nn.ge.2)
2702 do i=nptl+1,nptl+nn-1
2706 elseif(ijadu.eq.2) then
2711 if (a2.lt.a2min) then
2719 if (iflag.eq.0) then
2721 pptl(i,i1)=pptl(i,i1)+pptl(i,j1)
2724 if(istptl(i).eq.0.and.jorptl(i).eq.j1-nptl)
2726 if(istptl(i).eq.0.and.jorptl(i)+nptl.gt.j1)
2727 & jorptl(i)=jorptl(i)-1
2731 pptl(j,i)=pptl(j,i+1)
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))
2743 do i=nptl+1,nptl+nn-1
2745 if(pptl(4,jorptl(i)+nptl).lt.pptl(4,jorptl(j)+nptl))then
2753 idptl(nptl+jorptl(i))=9910+i-nptl
2756 jorptl(i)=0 !jorptl(nptl+jorptl(i))
2764 subroutine idtr5(id,ic)
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)
2775 ic(ii)=ic(ii)+10**(6-j)
2781 function ammin(id1,id2)
2782 dimension ic1(2),ic2(2),jc2(6,2),jc1(6,2)
2787 call iddeco(ic1,jc1)
2788 call iddeco(ic2,jc2)
2789 ammin=utamnx(jc1,jc2)
2793 c function idtr(id1,id2)
2794 c integer ic(2),id(2)
2802 c if(id(j).lt.0)ii=2
2804 c if(iabs(id(j)).gt.999)i1=3
2805 c do i=i1,int(log(abs(real(id(j))))/log(10.))+1
2806 c jj=mod(iabs(id(j))/10**(i-1),10)
2808 c ic(ii)=ic(ii)+10**(6-jj)
2812 c idtr=idtra(ic,0,0,3)
2813 c if(idtr.ne.idsp(id1,id2))then
2814 c write (*,'(4i6)') idtr,idsp(id1,id2),id1,id2
2819 function idsp(id1,id2)
2822 if(ia1.ge.1000.and.ia2.ge.1000)then
2825 elseif(ia1.le.1000.and.ia2.le.1000)then
2826 idsp=min(ia1,ia2)*100+max(ia1,ia2)*10
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
2833 if(id1.lt.0.and.id2.lt.0)isign=-1
2841 ida2=mod(ida/100,10)
2843 idsp=idb*1000+ida/10
2844 elseif(idb.le.ida2)then
2845 idsp=ida1*1000+idb*100+ida2*10
2849 if(ida1.eq.ida2.and.ida2.eq.idb)idsp=idsp+1