c----------------------------------------------------------------------- subroutine utresc(iret) c----------------------------------------------------------------------- c if irescl=1 rescaling is done, otherwise the purpose of going through c this routine is to not change the seed in case of no rescaling c----------------------------------------------------------------------- include 'epos.inc' include 'epos.incems' double precision p1,esoll,ppp,seedp,psoll,pt1soll,pt2soll dimension p1(5),p2(4),p0(5,mamx+mamx),pini(mxptl) logical force,nolead(mxptl),lspec(mxptl),lim data scalmean/0./scalevt/0./ save scalmean,scalevt call utpri('utresc',ish,ishini,4) errlim=0.005 !max(0.001,1./engy) if(iLHC.eq.1)errlim=max(0.00005,0.5/engy) iret=0 nptlpt=iabs(maproj)+iabs(matarg) call ranfgt(seedp) !not to change the seed ... if(nptl.le.nptlpt) goto 9999 if(ish.ge.8)then call alistf('list before boost&') endif esoll=0.d0 psoll=0.d0 p1(1)=0.d0 p1(2)=0.d0 p1(3)=0.d0 p1(4)=0.d0 p2(3)=0.d0 p2(4)=0.d0 ipmax=4 imin=nptlpt+1 if(iappl.eq.1)then imin=1 ipmax=2 if(iLHC.eq.1)ipmax=0 endif c store projectile and target in p0 and sum pz an E in p1(3) and p1(4) do i=1,nptlpt nolead(i)=.false. do j=1,5 p0(j,i)=dble(pptl(j,i)) enddo c calculate total energy of primaries c if(mod(istptl(i),10).eq.1)then do j=ipmax+1,4 p1(j)=p1(j)+dble(pptl(j,i)) enddo c endif c calculate total energy of secondaries if(mod(istptl(i),10).eq.0)then do j=ipmax+1,4 p2(j)=p2(j)+dble(pptl(j,i)) enddo endif enddo c fix secondaries counted in the system do i=nptlpt+1,nptl if(mod(istptl(i),10).eq.0)then c check maximum energy if(iLHC.eq.1.and.pptl(4,i).gt.engy*0.51)then if(ish.ge.1)write(ifch,*)'Energy of particle too high !' & ,i,ityptl(i),pptl(4,i) c call utstop('utresc&') if(ityptl(i).eq.48.or.ityptl(i).eq.58 !diffractive resonance & .or.ityptl(i).eq.47.or.ityptl(i).eq.57)then !active spectators pptl(4,i)=0.5*engy amt=sqrt(pptl(1,i)**2+pptl(2,i)**2+pptl(5,i)**2) pptl(3,i)=(pptl(4,i)+amt)*(pptl(4,i)-amt) if(pptl(3,i).gt.0.)then pptl(3,i)=sqrt(pptl(3,i)) else iret=1 endif else iret=1 endif if(iret.eq.1)then if(ish.ge.1)write(ifch,*)'Warning in utresc: redo event...' c call utstop('Energy of particle too high !&') goto 9999 endif endif c fix pt (p1(1) and p2(1)) from secondaries do j=1,ipmax p1(j)=p1(j)+dble(pptl(j,i)) enddo c calculate total energy of secondaries do j=ipmax+1,4 p2(j)=p2(j)+dble(pptl(j,i)) enddo lspec(i)=.false. if(((ityptl(i).eq.45.or.ityptl(i).eq.48).and.maproj.ge.100) & .or.((ityptl(i).eq.55.or.ityptl(i).eq.58).and.matarg.ge.100)) & lspec(i)=.true. if((abs(pptl(3,i)/pnullx).le.0.9 & .and.abs(pptl(3,i)).gt.pptl(5,i)).or.lspec(i))then nolead(i)=.true. c write(ifch,*)'nolead',i else nolead(i)=.false. c write(ifch,*)'lead',i endif endif enddo psoll=max(dble(pnullx),abs(p1(3))) c check if energy is conserved before boost if(iappl.eq.1)then diff=abs(p2(3)-p1(3)) scal=p2(4)/p1(4) if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim & .and.(iLHC.eq.0.or. & (abs(p2(1)).lt.errlim.and.abs(p2(2)).lt.errlim)))then if(ish.ge.4) & write (ifch,'(2x,a,2g14.6)') 'Energy OK: ',scal,abs(diff/psoll) goto 9999 else diff=0. scal=1. endif endif c calculate boost vector to have sum of pt=0 and sum of pz=0 ppp=(p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)*p1(2)-p1(1)*p1(1) if(ppp.gt.0.d0)then p1(5)=sqrt(ppp) else iret=1 write(ifch,*)'p1',p1(1),p1(2),p1(3),p1(4),ppp if(ish.ge.2)write (ifch,*) 'Problem in utresc (1): redo event' write (ifmt,*) 'Problem in utresc (1): redo event' c call utstop('utresc&') goto 9999 endif esoll=p1(5) if(ish.ge.4) write (ifch,'(a,5g14.6)') 'boost-vector: ',p1 c trafo c ----- pmax=0.d0 npart=0 pt1soll=0.d0 pt2soll=0.d0 do i=imin,nptl if(mod(istptl(i),10).le.1)then call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) if(mod(istptl(i),10).eq.0.and.i.gt.nptlpt)then npart=npart+1 pt1soll=pt1soll+dble(pptl(1,i)) pt2soll=pt2soll+dble(pptl(2,i)) endif endif if(i.gt.nptlpt.and.nolead(i))pmax=max(pmax,abs(pptl(3,i))) enddo if(ish.ge.6)then call alistf('list after boost&') endif if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------' if(iLHC.eq.1)then if(ish.ge.6)write(ifch,*)'ptsoll:',pt1soll,pt2soll,npart pt1soll=pt1soll/dble(npart) pt2soll=pt2soll/dble(npart) do i=nptlpt+1,nptl if(mod(istptl(i),10).eq.0)then pptl(1,i)=pptl(1,i)-sngl(pt1soll) pptl(2,i)=pptl(2,i)-sngl(pt2soll) pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2 & +pptl(3,i)**2+pptl(5,i)**2) endif enddo endif if(ish.ge.6)write(ifch,*)'esoll,psoll,pmax:',esoll,psoll,pmax c rescale momenta in rest frame c ----------------------------- scal=1. diff0=0. c ndif0=1 ferr=0.05 force=.false. npart=nptl-imin+1 do ipass=1,300 sum=0. sum3=0. difft=diff0 ndif=0 nfirst=int(rangen()*float(npart)) !start list randomly do i=0,npart-1 j=imin+i+nfirst if(j.gt.nptl)j=imin+i+nfirst-npart if(mod(istptl(j),10).eq.0)then c modified particles if(nolead(j))then c if(j.gt.nptlpt)then c if(abs(pptl(3,j))/pnullx.lt.0.9)then !not spectator or diffraction if(scal.eq.1..and.abs(diff0).lt.1.e-6)then ndif=ndif+1 pini(j)=pptl(3,j) else pptl3new=0. if( force .or.( & ityptl(j)/10.eq.4.or.ityptl(j)/10.eq.5 & ))then !try just remnant first ndif=ndif+1 diff=sign(min(0.3*abs(pini(j)), & rangen()*abs(difft)),difft) pptl3new=scal*(pptl(3,j)-diff) c write(ifch,*)'par',j,pptl3new,pptl(3,j),diff,difft c & ,ndif,pmax,scal if(abs(pptl3new).lt.pmax)then c particle should not be too fast or too modified if(abs(pptl3new-pini(j)).lt.ferr*abs(pini(j)) & .or.(lspec(j).and.abs(pptl3new).lt.abs(0.8*pini(j))))then c write(ifch,*)'used' difft=difft-diff pptl(3,j)=scal*(pptl(3,j)-diff) pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2 * +pptl(3,j)**2+pptl(5,j)**2) endif endif endif endif endif c sum over all particles sum=sum+pptl(4,j) sum3=sum3+pptl(3,j) endif enddo diff=sum3 scal=sngl(esoll)/(sum-diff) if(ish.ge.6)write(ifch,*) $ 'ipass,scal,diff/psoll,e,pz,ndif,f:' $ ,ipass,scal,diff/psoll,sum,sum3,ndif,force if(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim) $ goto 300 if(ndif.gt.0.and.(force.or.ipass.lt.150))then c ndif0=ndif diff0=diff elseif(abs(scal-1.).le.1e-2.and.abs(diff/psoll).lt.5e-2 & .and.iLHC.eq.0)then goto 300 elseif(force)then if(ish.ge.2) $ write(ifmt,*)'Warning in utresc: no more allowed particle' goto 302 else force=.true. ferr=0.1 diff=0. endif enddo 302 if(iLHC.eq.1)then lim=.not.(abs(scal-1.).le.errlim.and.abs(diff/psoll).lt.errlim) else lim=abs(scal)+abs(diff/psoll).gt.2.5 endif if(ish.ge.1)then call utmsg('utrescl') write(ifch,*)'***** scal=',scal,diff/psoll,lim call utmsgf endif if(lim)then if(ish.ge.1)then write(ifmt,*)'Warning in utresc !' write(ifch,'(a,i10,d25.15)')'redo EPOS event ...' & ,nint(seedj),seedc endif iret=1 goto 9999 endif c trafo c ----- 300 continue do i=1,nptl if(i.le.nptlpt)then do j=1,5 pptl(j,i)=p0(j,i) enddo else if(mod(istptl(i),10).le.1)then call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) endif endif enddo if(ish.ge.4)call alist('list after rescaling&',1,nptl) 9999 continue if(ish.ge.2)then scalevt=scalevt+1. scalmean=scalmean+scal write(ifch,*)' average rescaling factor: ',scalmean & /scalevt endif call ranfst(seedp) ! ... after this subroutine call utprix('utresc',ish,ishini,4) end c----------------------------------------------------------------------- subroutine utghost(iret) c----------------------------------------------------------------------- c if irescl=1 make particle on-shell if not c----------------------------------------------------------------------- include 'epos.inc' include 'epos.incems' double precision seedp call utpri('ughost',ish,ishini,4) iret=0 nptlpt=iabs(maproj)+iabs(matarg) if(iappl.eq.6.or.iappl.eq.8)nptlpt=3 ! ee or DIS call ranfgt(seedp) !not to change the seed ... if(nptl.le.nptlpt) goto 9999 if(ish.ge.5)write(ifch,'(a)')'---------mark ghosts---------' c mark ghosts c ----------- do j=nptlpt+1,nptl if(istptl(j).le.1.and.pptl(4,j).gt.0.d0)then if(iLHC.eq.0.or.mod(abs(idptl(j)),10).le.1)then !for LHC tune don't fix mass of resonnances (to keep width) amass=pptl(5,j) call idmass(idptl(j),amass) if(abs(idptl(j)).gt.100.and. & abs(pptl(5,j)-amass).gt.0.01*amass)then if(ish.ge.5)write(ifch,*)'wrong particle mass',j,idptl(j) & ,pptl(5,j),amass amass=pptl(5,j) call idres(idptl(j),amass,idr,iadj) if(idr.ne.0)then pptl(5,j)=amass idptl(j)=idr else call idmass(idptl(j),amass) pptl(5,j)=amass endif call idtau(idptl(j),pptl(4,j),pptl(5,j),taugm) tivptl(2,j)=tivptl(1,j)+taugm*(-alog(rangen())) else pptl(5,j)=amass endif endif if(abs((pptl(4,j)+pptl(3,j))*(pptl(4,j)-pptl(3,j)) $ -pptl(2,j)**2-pptl(1,j)**2-pptl(5,j)**2).gt.0.3 $ .and.abs(1.-abs(pptl(3,j))/pptl(4,j)).gt.0.01)then !print*,'ghost',ityptl(j),idptl(j) if(ish.ge.1)write(ifmt,*)'ghost:',j,idptl(j),ityptl(j) if(ish.ge.5)then write(ifch,'(a,$)')'ghost:' call alistc("&",j,j) endif ityptl(j)=100+ityptl(j)/10 elseif(irescl.ge.1)then c ensure that all particles are really on-shell pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2 * +pptl(3,j)**2+pptl(5,j)**2) endif elseif(mod(istptl(j),10).eq.0)then c if not droplet with fusion if(istptl(j).ne.10.or.iorsdf.ne.3)then if(ish.ge.1)then write(ifmt,*)'Lost particle (E=0)' write(ifch,*)'Lost particle (E=0) :' call alistc("utghost&",j,j) endif istptl(j)=istptl(j)+2 endif endif enddo if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------' c treat ghosts c ------------ ifirst=1 scal=1. pfif=0. efif=0. ntry=0 132 nfif=0 psum=0 esum=0. ntry=ntry+1 do j=nptlpt+1,nptl if(mod(istptl(j),10).eq.0)then if(ityptl(j).gt.100)then nfif=nfif+1 if(ifirst.eq.1)then pfif=pfif+pptl(3,j) if(pptl(4,j).gt.0.)efif=efif+pptl(4,j) endif if(irescl.ge.1) then if(ifirst.gt.1)then if(pptl(4,j).gt.0.)then Einv=1./pptl(4,j) amt=1.-(pptl(5,j)*Einv)**2+(pptl(1,j)*Einv)**2 $ +(pptl(2,j)*Einv)**2 else amt=-1. endif if(amt.gt.0.)then pptl(3,j)=sign(pptl(4,j),pptl(3,j))*sqrt(amt) else y=(rangen()+rangen()+rangen()+rangen()-2.)/2.*yhaha y=sign(abs(y),pptl(3,j)) pptl(3,j) $ =sqrt(pptl(5,j)**2+pptl(1,j)**2 $ +pptl(2,j)**2)*sinh(y) pptl(4,j) $ =sqrt(pptl(5,j)**2+pptl(1,j)**2 $ +pptl(2,j)**2)*cosh(y) efif=efif+pptl(4,j) endif ifirst=0 else c do k=1,3 do k=3,3 pptl(k,j)=pptl(k,j)*scal enddo pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2 * +pptl(5,j)**2) endif endif psum=psum+pptl(3,j) esum=esum+pptl(4,j) if(ish.ge.5) $ write (ifch,*) 'nrevt,psum,esum,pfif,efif,nfif,scal' $ ,nrevt,psum,esum,pfif,efif,nfif,scal endif endif enddo if ( ish.ge.5 ) write (ifch,*) 'tot',nfif,efif,pfif,esum,psum if(nfif.gt.5.or.(esum.gt.0.05*engy.and.nfif.ne.1))then if(ifirst.eq.0)then do j=nptlpt+1,nptl if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then if((psum-pfif)*(1.-scal).ge.0) & pptl(3,j)=pptl(3,j)-(psum-pfif)/nfif endif enddo else ifirst=2 goto 132 endif scal=efif/esum if ( ish.ge.5 ) write (ifch,*) 'scal',scal if ( abs(scal-1.) .gt. 0.05 ) then if(ntry.le.1000)then goto 132 else iret=1 if(ish.ge.2)write (ifch,*) 'Problem in utghost : redo event' if(ish.ge.1)write (ifmt,*) 'Problem in utghost : redo event' goto 9999 endif endif else do j=nptlpt+1,nptl if ( ityptl(j).ge.101 .and. ityptl(j).le.105 )then pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2 * +pptl(5,j)**2) endif enddo endif if(ish.ge.5)write(ifch,'(a)')'---------Check Ghost list---------' c Check Ghost list if(ish.ge.5)then do j=nptlpt+1,nptl if(mod(istptl(j),10).eq.0)then if(ityptl(j).le.105.and.ityptl(j).ge.101)then write(ifch,'(a,$)')'ghost:' call alistc("&",j,j) endif endif enddo endif 9999 continue call ranfst(seedp) ! ... after this subroutine call utprix('ughost',ish,ishini,4) end c----------------------------------------------------------------------- subroutine utrsph(iret) c----------------------------------------------------------------------- c if irescl=1 and ispherio=1 rescaling is done for particle used by c spherio as initial condition. c----------------------------------------------------------------------- include 'epos.inc' include 'epos.incems' double precision p1,esoll dimension p1(5),p0(5,mamx+mamx) call utpri('utrsph',ish,ishini,4) errlim=0.0001 iret=0 nptlpt=iabs(maproj)+iabs(matarg) if(nptl.le.nptlpt) goto 9999 esoll=0.d0 p1(1)=0.d0 p1(2)=0.d0 p1(3)=0.d0 p1(4)=0.d0 do i=nptlpt+1,nptl if((istptl(i).le.11 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41)) $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then do j=1,2 p1(j)=p1(j)+dble(pptl(j,i)) enddo endif enddo do i=1,nptlpt do j=1,5 p0(j,i)=pptl(j,i) enddo do j=3,4 p1(j)=p1(j)+dble(pptl(j,i)) enddo enddo p1(5)=dsqrt((p1(4)+p1(3))*(p1(4)-p1(3))-p1(2)**2.d0-p1(1)**2.d0) esoll=p1(5) if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector',p1 c trafo c ----- do i=1,nptl if((istptl(i).le.11 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41)) $ .or.istptl(i).eq.20.or.istptl(i).eq.21 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then call utlob4(1,p1(1),p1(2),p1(3),p1(4),p1(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) endif enddo if(ish.ge.5)write(ifch,'(a)')'------------------' c rescale momenta in rest frame c ----------------------------- scal=1. diff=0. do ipass=1,1000 sum=0. sum3=0. ndif=0 do j=1,nptl if((istptl(j).le.11 $ .and.(iorptl(j).ge.1.and.istptl(iorptl(j)).eq.41)) $ .or.istptl(j).eq.20.or.istptl(j).eq.21 $ .or.(istptl(j).eq.0.and.j.le.nptlpt))then if(j.gt.nptlpt)then ndif=ndif+1 pptl(3,j)=scal*(pptl(3,j)-diff) pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2+pptl(3,j)**2 * +pptl(5,j)**2) endif sum=sum+pptl(4,j) sum3=sum3+pptl(3,j) endif enddo diff=sum3/real(ndif) scal=real(esoll)/sum if(ish.ge.6)write(ifch,*)'ipass,scal,diff,e,esoll,pz,ndif:' $ ,ipass,scal,diff,sum,esoll,sum3,ndif if(abs(scal-1.).le.errlim.and.abs(diff).lt.10.*errlim) goto300 enddo if(ish.ge.1)then call utmsg('hresph') write(ifch,*)'***** scal=',scal,diff call utmsgf endif c trafo c ----- 300 continue c do i=nptlpt+1,nptl do i=1,nptl if((istptl(i).le.11 $ .and.(iorptl(i).ge.1.and.istptl(iorptl(i)).eq.41)) $ .or.istptl(i).eq.20.or.istptl(i).eq.21 $ .or.(istptl(i).eq.0.and.i.le.nptlpt))then call utlob4(-1,p1(1),p1(2),p1(3),p1(4),p1(5) $ ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i)) endif if(i.le.nptlpt)then do j=1,5 pptl(j,i)=p0(j,i) enddo endif enddo 9999 call utprix('utrsph',ish,ishini,4) end cc----------------------------------------------------------------------- c double precision function dddlog(xxx) cc----------------------------------------------------------------------- c double precision xxx c dddlog=-1d50 c if(xxx.gt.0d0)dddlog=dlog(xxx) c end c ccc----------------------------------------------------------------------- c subroutine randfl(jc,iqa0,iflav,ic,isame) cc----------------------------------------------------------------------- cc returns random flavour ic(2) (iqa0=1:quark,2:antiquark,11:diquark) cc----------------------------------------------------------------------- c include 'epos.inc' c real probab(nflav),probsu(nflav+1) c integer jc(nflav,2),jc0(nflav,2),ic(2) c if(ish.ge.6)then c write(ifch,*)('-',i=1,10) c *,' entry sr randfl ',('-',i=1,30) c write(ifch,*)'iqa0:',iqa0 c write(ifch,*)'jc:' c write(ifch,*)jc c endif c iflav=0 c ic(1)=0 c ic(2)=0 c do 10 n=1,nflav c do 10 i=1,2 c10 jc0(n,i)=0 c iqa1=iqa0*10 c9999 iqa1=iqa1/10 c if(iqa1.eq.0)goto9998 c iqa=mod(iqa1,10) c su=0 c do 20 i=1,nflav c probab(i)=jc(i,iqa)-jc0(i,iqa) c if(isame.eq.1)probab(i)=probab(i)*(jc(i,3-iqa)-jc0(i,3-iqa)) c20 su=su+probab(i) c if(su.lt..5)then c iflav=0 c ic(1)=0 c ic(2)=0 c goto9998 c endif c probsu(1)=0. c do 30 i=1,nflav c probsu(i+1)=probsu(i)+probab(i)/su c if(probsu(i+1)-probsu(i).lt.1e-5)probsu(i+1)=probsu(i) c30 continue c r=rangen()*probsu(nflav+1) c do 50 i=1,nflav c if(probsu(i).le.r.and.r.lt.probsu(i+1))iflav=i c50 continue c jc0(iflav,iqa)=jc0(iflav,iqa)+1 c if(isame.eq.1)jc0(iflav,3-iqa)=jc0(iflav,3-iqa)+1 c call idenco(jc0,ic,ireten) c if(ireten.eq.1)call utstop('randfl: idenco ret code = 1&') c if(ish.ge.6)then c write(ifch,*)'probab:' c write(ifch,*)probab c write(ifch,*)'probsu:' c write(ifch,*)probsu c write(ifch,*)'ran#:',r,' flav:',iflav c endif c goto9999 c9998 continue c if(ish.ge.6)write(ifch,*)('-',i=1,30) c *,' exit sr randfl ',('-',i=1,10) c return c end c c cc----------------------------------------------------------------------- c subroutine ranhvy(x,eps) cc----------------------------------------------------------------------- cc generates x for heavy particle fragmentation according to cc the peterson form cc d(x)=1/(x*(1-1/x-eps/(1-x))**2) cc =d0(x)*d1(x)*d2(x) cc d0(x)=(1-x)**2/((1-x)**2+eps)**2 cc d1(x)=x cc d2(x)=(((1-x)**2+eps)/((1-x)**2+eps*x))**2 cc using x=1-y**pow cc generates flat in x if eps>1. cc----------------------------------------------------------------------- c data aln4/1.3863/ c if(eps.lt.1.) then c pow=alog((3.+eps)/eps)/aln4 c ymx=(eps*(3.*pow-1.)/(pow+1.))**(.5/pow) c zmx=1-ymx**pow c d0mx=(1-zmx)**2/((1.-zmx)**2+eps)**2*pow*ymx**(pow-1.) c d2mx=2./(2.-sqrt(eps)) c else c pow=1. c zmx=0. c d0mx=(1.-zmx)**2/((1.-zmx)**2+eps)**2 c d2mx=1.+eps c endif cc cc generate z according to (1-z)**2/((1-z)**2+eps*z)**2 c1 continue c y=rangen() c z=1.-y**pow cc c d0z=(1.-z)**2/((1.-z)**2+eps)**2*pow*y**(pow-1.) c if(d0z.lt.rangen()*d0mx) goto1 cc cc check remaining factors c d1=z c d2=(((1.-z)**2+eps)/((1.-z)**2+eps*z))**2 c if(d1*d2.lt.rangen()*d2mx) goto1 cc cc good x c x=z c return c end c *-- Author : D. HECK IK FZK KARLSRUHE 27/04/1994 C======================================================================= SUBROUTINE EPOVAPOR( MAPRO,INEW,JFIN,ITYP,PFRX,PFRY,PFRZ ) C----------------------------------------------------------------------- C (E)VAPOR(ATION OF NUCLEONS AND ALPHA PARTICLES FROM FRAGMENT) C C TREATES THE REMAINING UNFRAGMENTED NUCLEUS C EVAPORATION FOLLOWING CAMPI APPROXIMATION. C SEE: X. CAMPI AND J. HUEFNER, PHYS.REV. C24 (1981) 2199 C AND J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990) C THIS SUBROUTINE IS CALLED FROM SDPM, DPMJST, NSTORE, AND VSTORE. C ARGUMENTS INPUT: C MAPRO = NUMBER OF NUCLEONS OF PROJECTILE C INEW = PARTICLE TYPE OF SPECTATOR FRAGMENT C ARGUMENTS OUTPUT: C JFIN = NUMBER OF FRAGMENTS C ITYP(1:JFIN) = NATURE (PARTICLE CODE) OF FRAGMENTS C PFRX(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN X-DIRECTION C PFRY(1:JFIN) = TRANSVERSE MOMENTUM OF FRAGMENTS IN Y-DIRECTION C PFRZ(1:JFIN) = LONGITUDINAL MOMENTUM OF FRAGMENTS IN Z-DIRECTION C C FROM CORSIKA AND ADAPTED BY T. PIEROG TO INCLUDE LONG MOMENTUM AND C MORE REALISTIC FRAGMENTS C----------------------------------------------------------------------- IMPLICIT NONE include 'epos.inc' common/eporansto2/irndmseq integer irndmseq DOUBLE PRECISION PFR(mamxx),PFRX(mamxx),PFRY(mamxx),PFRZ(mamxx) * ,RD(2*mamxx),SPFRY,SPFRZ,drangen DOUBLE PRECISION AFIN,AGLH,APRF,BGLH,EEX,PHIFR,RANNORM,SPFRX INTEGER ITYP(mamxx),IARM,INEW,ITYPRM,INRM,IS,IZRM,JC,lseq * ,JFIN,K,L,LS,MAPRO,MF,NFIN,NINTA,NNUC,NPRF,NNSTEP SAVE EXTERNAL RANNORM C----------------------------------------------------------------------- IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : MAPRO,INEW=',MAPRO,INEW JC = 0 JFIN = 0 lseq = irndmseq ITYPRM = INEW NPRF = INEW/100 NINTA = MAPRO - NPRF IF ( NINTA .EQ. 0 ) THEN C NO NUCLEON HAS INTERACTED JFIN = 1 PFRX(1) = 0.D0 PFRY(1) = 0.D0 PFRZ(1) = 0.D0 ITYP(1) = ITYPRM IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,NINTA=',JFIN,NINTA GOTO 50 ENDIF C EXCITATION ENERGY EEX OF PREFRAGMENT C SEE: J.J. GAIMARD, THESE UNIVERSITE PARIS 7, (1990), CHPT. 4.2 EEX = 0.D0 CALL RMMARD( RD,2*NINTA,lseq ) DO L = 1, NINTA IF ( RD(NINTA+L) .LT. RD(L) ) RD(L) = 1.D0 - RD(L) EEX = EEX + RD(L) ENDDO C DEPTH OF WOODS-SAXON POTENTIAL TO FERMI SURFACE IS 0.040 GEV IF(ish.ge.7)WRITE(ifch,*)'EPOVAPOR : EEX=',SNGL(EEX*0.04D0), & ' GEV' C EVAPORATION: EACH EVAPORATION STEP NEEDS ABOUT 0.020 GEV, THEREFORE C NNSTEP IS EEX * 0.04/0.02 = EEX * 2. NNSTEP = INT( EEX*2.D0 ) IF ( NNSTEP .LE. 0 ) THEN C EXCITATION ENERGY TOO SMALL, NO EVAPORATION JFIN = 1 PFRX(1) = 0.D0 PFRY(1) = 0.D0 PFRZ(1) = 0.D0 ITYP(1) = ITYPRM IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : JFIN,EEX=',JFIN,SNGL(EEX) GOTO 50 ENDIF C AFIN IS ATOMIC NUMBER OF FINAL NUCLEUS APRF = DBLE(NPRF) AFIN = APRF - 1.6D0 * DBLE(NNSTEP) NFIN = MAX( 0, INT( AFIN+0.5D0 ) ) C CORRESPONDS TO DEFINITION; FRAGMENTATION-EVAPORATION C CONVOLUTION EMU07 /MODEL ABRASION EVAPORATION (JNC FZK APRIL 94) C NNUC IS NUMBER OF EVAPORATING NUCLEONS NNUC = NPRF - NFIN IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NFIN,NNUC=',NFIN,NNUC IF ( NNUC .LE. 0 ) THEN C NO EVAPORATION JFIN = 1 PFRX(1) = 0.D0 PFRY(1) = 0.D0 PFRZ(1) = 0.D0 ITYP(1) = ITYPRM GOTO 50 ELSEIF ( NNUC .GE. 4 ) THEN C EVAPORATION WITH FORMATION OF ALPHA PARTICLES POSSIBLE C IARM, IZRM, INRM ARE NUMBER OF NUCLEONS, PROTONS, NEUTRONS OF C REMAINDER DO LS = 1, NNSTEP IARM = ITYPRM/100 IF ( IARM .LE. 0 ) GOTO 100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 CALL RMMARD( RD,2,lseq ) IF ( RD(1) .LT. 0.2D0 .AND. IZRM .GE. 2 * .AND. INRM .GE. 2 ) THEN ITYP(JC) = -402 !alpha NNUC = NNUC - 4 ITYPRM = ITYPRM - 402 ELSE IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. INRM ) THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSEIF(RD(2)*(IZRM+INRM).LT.IZRM.AND.IZRM.GE.INRM)THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSE ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ENDIF NNUC = NNUC - 1 ENDIF IF ( NNUC .LE. 0 ) GOTO 50 ENDDO ENDIF IF ( NNUC .LT. 4 ) THEN C EVAPORATION WITHOUT FORMATION OF ALPHA PARTICLES CALL RMMARD( RD,NNUC,lseq ) DO IS = 1, NNUC IARM = ITYPRM/100 IF ( IARM .LE. 0 ) GOTO 100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 IF ( IZRM .EQ. 1 .AND. INRM .GT. IZRM ) THEN ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ELSEIF ( INRM .EQ. 1 .AND. IZRM .GT. IZRM ) THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSEIF ( RD(IS)*IARM .LT. IZRM .AND. IZRM .GE. INRM ) THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSE ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ENDIF ENDDO ENDIF 50 CONTINUE IARM = ITYPRM/100 IF ( IARM .LE. 0 ) GOTO 100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 C CLEAN FRAGMENT TO HAVE ALWAYS AT LEAST TWO TIMES MORE N THAN P CALL RMMARD( RD,2,lseq ) DO WHILE ( INRM .GT. INT( 1.15D0 * DBLE(IZRM) + RD(1) ) ) ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 IARM = ITYPRM/100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 ENDDO C CLEAN FRAGMENT NOT TO HAVE TOO MANY P DO WHILE ( IZRM .GE. NINT( (DBLE(INRM)+1.D0+RD(2)) / 1.15D0 ) ) ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 IARM = ITYPRM/100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 ENDDO IF ( IARM .EQ. 8 ) THEN !EXCLUDED DO WHILE ( IZRM .GE. 2 .AND. INRM .GE. 2 ) ITYP(JC) = -402 ITYPRM = ITYPRM - 402 IARM = ITYPRM/100 IZRM = MOD(ITYPRM,100) INRM = IARM - IZRM JC = JC + 1 ENDDO IF ( ITYPRM .GT. 0 ) THEN IF ( ITYPRM/100 .EQ. 5 ) THEN !EXCLUDED BUT SHOULD NOT HAPPEN HERE IF ( IZRM .GE. INRM ) THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSE ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ENDIF JC = JC + 1 ENDIF ITYP(JC) = -ITYPRM ELSE JC = JC - 1 ENDIF ELSEIF ( IARM .EQ. 5 ) THEN !EXCLUDED IF ( IZRM .GE. INRM ) THEN ITYP(JC) = 1120 ITYPRM = ITYPRM - 101 ELSE ITYP(JC) = 1220 ITYPRM = ITYPRM - 100 ENDIF JC = JC + 1 ITYP(JC) = -ITYPRM ELSEIF ( ITYPRM .GT. 200 ) THEN ITYP(JC) = -ITYPRM ELSEIF ( ITYPRM .EQ. 101 ) THEN ITYP(JC) = 1120 ELSEIF ( ITYPRM .EQ. 100 ) THEN ITYP(JC) = 1220 ELSE JC = JC - 1 IF ( ITYPRM .NE. 0 ) WRITE(*,*) * 'EPOVAPOR : ILLEGAL PARTICLE ITYPRM =',ITYPRM ENDIF 100 CONTINUE IF ( JC .GT. JFIN ) THEN JFIN = JC IF(ish.ge.7)WRITE(ifch,*) * 'EPOVAPOR : NO ITYP PFR PFL' IF ( infragm .EQ. 2 ) THEN C EVAPORATION WITH PT AFTER PARAMETRIZED JACEE DATA DO MF = 1, JFIN IF(ITYP(MF).LT.0)THEN IARM=-ITYP(MF)/100 ELSE IARM=1 ENDIF PFR(MF) = RANNORM(0.088D0,0.044D0) PFRZ(MF)= (2*int(drangen(PFR(MF))+0.5d0)-1) & *RANNORM(0.300D0/DBLE(IARM),0.100D0/SQRT(DBLE(IARM))) !Fermi motion about 300 MeV IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF)) & ,SNGL(PFRZ(MF)) ENDDO ELSEIF ( infragm .EQ. 3 ) THEN C EVAPORATION WITH PT AFTER GOLDHABER''S MODEL (PHYS.LETT.53B(1974)306) DO MF = 1, JFIN K = MAX( 1, -ITYP(MF)/100 ) BGLH = K * (MAPRO - K) / DBLE(MAPRO-1) C THE VALUE 0.103 [GEV] IS SIGMA(0)=P(FERMI)/SQRT(5.) * AGLH = 0.103D0 * SQRT( BGLH ) C THE VALUE 0.090 [GEV] IS EXPERIMENTALLY DETERMINED SIGMA(0) AGLH = 0.090D0 * SQRT( BGLH ) PFR(MF) = RANNORM(0.D0,AGLH) PFRZ(MF)= RANNORM(0.000D0,0.500D0) !from pAg at 100 GeV IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF)) & ,SNGL(PFRZ(MF)) ENDDO ELSE C EVAPORATION WITHOUT TRANSVERSE MOMENTUM DO MF = 1, JFIN PFR(MF) = 0.D0 PFRZ(MF)= 0.D0 IF(ish.ge.7)WRITE(ifch,*) MF,ITYP(MF),SNGL(PFR(MF)) & ,SNGL(PFRZ(MF)) ENDDO ENDIF C CALCULATE RESIDUAL TRANSVERSE MOMENTUM SPFRX = 0.D0 SPFRY = 0.D0 SPFRZ = 0.D0 CALL RMMARD( RD,JFIN,lseq ) DO MF = 1, JFIN PHIFR = PI * RD(MF) PFRX(MF) = PFR(MF) * COS( PHIFR ) PFRY(MF) = PFR(MF) * SIN( PHIFR ) SPFRY = SPFRY + PFRY(MF) SPFRX = SPFRX + PFRX(MF) SPFRZ = SPFRZ + PFRZ(MF) ENDDO C CORRECT ALL TRANSVERSE MOMENTA FOR MOMENTUM CONSERVATION SPFRX = SPFRX / JFIN SPFRY = SPFRY / JFIN SPFRZ = SPFRZ / JFIN DO MF = 1, JFIN PFRX(MF) = PFRX(MF) - SPFRX PFRY(MF) = PFRY(MF) - SPFRY PFRZ(MF) = PFRZ(MF) - SPFRZ ENDDO ENDIF IF(ish.ge.7)WRITE(ifch,*) 'EPOVAPOR : NINTA,JFIN=',NINTA,JFIN RETURN END *-- Author : The CORSIKA development group 21/04/1994 C======================================================================= DOUBLE PRECISION FUNCTION RANNORM( A,B ) C----------------------------------------------------------------------- C RAN(DOM NUMBER) NOR(MALLY DISTRIBUTED) C C GENERATES NORMAL DISTRIBUTED RANDOM NUMBER C DELIVERS 2 UNCORRELATED RANDOM NUMBERS, C THEREFORE RANDOM CALLS ARE ONLY NECESSARY EVERY SECOND TIME. c but to be used with CONEX/CORSIKA we should always use 2 new number c to be able to reproduce the shower C REFERENCE : NUMERICAL RECIPES, W.H. PRESS ET AL., C CAMBRIDGE UNIVERSITY PRESS, 1992 ISBN 0 521 43064 X C ARGUMENTS: C A = MEAN VALUE C B = STANDARD DEVIATION C C FROM CORSIKA C----------------------------------------------------------------------- IMPLICIT NONE double precision facrdm,u1rdm,u2rdm,drangen c logical knordm c data knordm/.true./ DOUBLE PRECISION A,B,RR SAVE facrdm,u1rdm,u2rdm!,knordm C----------------------------------------------------------------------- c IF ( KNORdm ) THEN 1 CONTINUE U1rdm = 2.D0*drangen(a) - 1.D0 U2rdm = 2.D0*drangen(b) - 1.D0 RR = U1rdm**2 + U2rdm**2 IF ( RR .GE. 1.D0 .OR. RR .EQ. 0.D0 ) GOTO 1 FACrdm = SQRT( (-2.D0) * LOG(RR) / RR ) RANNORM = FACrdm * U1rdm * B + A c KNORdm = .FALSE. c ELSE c RANNORM = FACrdm * U2rdm * B + A c KNORdm = .TRUE. c ENDIF RETURN END c----------------------------------------------------------------------- function ransig() c----------------------------------------------------------------------- c returns randomly +1 or -1 c----------------------------------------------------------------------- ransig=1 if(rangen().gt.0.5)ransig=-1 return end cc----------------------------------------------------------------------- c function ranxq(n,x,q,xmin) cc----------------------------------------------------------------------- cc returns random number according to x(i) q(i) with x>=xmin cc----------------------------------------------------------------------- c include 'epos.inc' c real x(n),q(n) c imin=1 c if(xmin.eq.0.)goto3 c i1=1 c i2=n c1 i=i1+(i2-i1)/2 c if(x(i).lt.xmin)then c i1=i c elseif(x(i).gt.xmin)then c i2=i c else c imin=i c goto3 c endif c if(i2-i1.gt.1)goto1 c imin=i2 c3 continue c if(q(imin).gt.q(n)*.9999)then c ranxq=xmin c goto4 c endif c qran=q(imin)+rangen()*(q(n)-q(imin)) c ranxq=utinvt(n,x,q,qran) c4 continue c c if(ranxq.lt.xmin)then c call utmsg('ranxq ') c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin c write(ifch,*)'q(imin) q q(n):',q(imin),qran,q(n) c write(ifch,*)'x(imin) x x(n):',x(imin),ranxq,x(n) c call utmsgf c ranxq=xmin c endif c c return c end c cc ***** end r-routines cc ***** beg s-routines c cc----------------------------------------------------------------------- c function sbet(z,w) cc----------------------------------------------------------------------- c sbet=utgam1(z)*utgam1(w)/utgam1(z+w) c return c end c cc----------------------------------------------------------------------- c function smass(a,y,z) cc----------------------------------------------------------------------- cc returns droplet mass (in gev) (per droplet, not (!) per nucleon) cc according to berger/jaffe mass formula, prc35(1987)213 eq.2.31, cc see also c. dover, BNL-46322, intersections-meeting, tucson, 91. cc a: massnr, y: hypercharge, z: charge, cc----------------------------------------------------------------------- c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero c ymin=ym*a c zmin=cz/(dz/a+zm/a**(1./3.)) c smass=epsi*a+as*a**(2./3.)+(ac/a**(1./3.)+dz/a/2.)*(z-zmin)**2 c *+dy/a/2.*(y-ymin)**2 c return c end c cc----------------------------------------------------------------------- c subroutine smassi(theta) cc----------------------------------------------------------------------- cc initialization for smass. cc calculates parameters for berger/jaffe mass formula cc (prc35(1987)213 eq.2.31, see also c. dover, BNL-46322). cc theta: parameter that determines all parameters in mass formula. cc----------------------------------------------------------------------- c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero c thet=theta c c astr=.150 c pi=3.14159 c alp=1./137. c c co=cos(theta) c si=sin(theta) c bet=(1+co**3)/2. c rzero=si/astr/( 2./3./pi*(1+co**3) )**(1./3.) cctp060829 cs=astr/si c cz=-astr/si*( ( .5*(1+co**3) )**(1./3.)-1 ) c sigma=6./8./pi*(astr/si)**3*(co**2/6.-si**2*(1-si)/3.- c *1./3./pi*(pi/2.-theta-sin(2*theta)+si**3*alog((1+co)/si))) c c epsi=astr*((.5*(1+co**3))**(1./3.)+2)/si c as=4*pi*sigma*rzero**2 c ac=3./5.*alp/rzero c dz=astr/si*bet**(1./3.)*co**2* c *(co**4*(1+bet**(2./3.))+(1+bet)**2)/ c *( (2*co**2+bet**(1./3.))*(co**4*(1+bet**(2./3.))+(1+bet)**2)- c *(co**4+bet**(1./3.)*(1+bet))*((2*bet**(2./3.)-1)*co**2+1+bet) ) c dy=astr/6.*(1+co**3)**3/si* c *( 1+(1+co)/(4*(1+co**3))**(2./3.) )/ c *(co**6+co+co*(.5*(1+co**3))**(4./3.)) c zm=6*alp/(5*rzero) c ym=(1-co**3)/(1+co**3) c c return c end c cc----------------------------------------------------------------------- c subroutine smassp cc----------------------------------------------------------------------- cc prints smass. cc----------------------------------------------------------------------- c include 'epos.inc' c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero c real eng(14),ymi(14),zmi(14) c pi=3.14159 c write(ifch,*)'parameters of mass formula:' c write(ifch,*)'theta=',thet,' epsi=',epsi c write(ifch,*)'as=',as,' ac=',ac c write(ifch,*)'dy=',dy,' dz=',dz c write(ifch,*)'ym=',ym c write(ifch,*)'cz dz zm=',cz,dz,zm c write(ifch,*)'sigma**1/3=',sigma**(1./3.),' rzero=',rzero c write(ifch,*)'mass:' c write(ifch,5000)(j,j=1,14) c5000 format(5x,'a:',14i5) c do 4 j=1,14 c a=j c ymi(j)=ym*a c4 zmi(j)=cz/(dz/a+zm/a**(1./3.)) c write(ifch,5002)(ymi(j),j=1,14) c5002 format(1x,'ymin: ',14f5.2) c write(ifch,5003)(zmi(j),j=1,14) c5003 format(1x,'zmin: ',14f5.2) c do 2 i=1,15 c ns=11-i c do 3 j=1,14 c a=j c y=a-ns c z=0. c3 eng(j)=smass(a,y,z)/a c write(ifch,5001)ns,(eng(j),j=1,14) c5001 format(1x,'s=',i2,2x,14f5.2) c2 continue c write(ifch,*)'mass-mass(free):' c write(ifch,5000)(j,j=1,14) c do 5 i=1,15 c ns=11-i c do 6 j=1,14 c a=j c y=a-ns c z=0. c call smassu(a,y,z,ku,kd,ks,kc) c6 eng(j)=(smass(a,y,z)-utamnu(ku,kd,ks,kc,0,0,3))/a c write(ifch,5001)ns,(eng(j),j=1,14) c5 continue c c stop c end c cc----------------------------------------------------------------------- c subroutine smasst(kux,kdx,ksx,kcx,a,y,z) cc----------------------------------------------------------------------- cc input: kux,kdx,ksx,kcx = net quark numbers (for u,d,s,c quarks). cc output: massnr a, hypercharge y and charge z. cc----------------------------------------------------------------------- c sg=1 c if(kux+kdx+ksx+kcx.lt.0.)sg=-1 c ku=sg*kux c kd=sg*kdx c ks=sg*ksx c kc=sg*kcx c k=ku+kd+ks+kc c if(mod(k,3).ne.0)stop'noninteger baryon number' c a=k/3 c y=a-ks c nz=2*ku-kd-ks+2*kc c if(mod(nz,3).ne.0)stop'noninteger charge' c z=nz/3 c return c end c cc----------------------------------------------------------------------- c subroutine smassu(ax,yx,zx,ku,kd,ks,kc) cc----------------------------------------------------------------------- cc input: massnr ax, hypercharge yx and charge zx. cc output: ku,kd,ks,kc = net quark numbers (for u,d,s,c quarks). cc----------------------------------------------------------------------- c sg=1 c if(ax.lt.0.)sg=-1 c a=sg*ax c y=sg*yx c z=sg*zx c ku=nint(a+z) c kd=nint(a-z+y) c ks=nint(a-y) c kc=0 c return c end c cc----------------------------------------------------------------------- c function spoc(a,b,c,d,x) cc----------------------------------------------------------------------- cc power fctn with cutoff cc----------------------------------------------------------------------- c spoc=0 c if(a.eq.0..and.b.eq.0.)return c spoc =a+b*x**c c spoc0=a+b*d**c c spoc=amin1(spoc,spoc0) c spoc=amax1(0.,spoc) c return c end c c----------------------------------------------------------------------- function utacos(x) c----------------------------------------------------------------------- c returns acos(x) for -1 <= x <= 1 , acos(+-1) else c----------------------------------------------------------------------- include 'epos.inc' argum=x if(x.lt.-1.)then if(ish.ge.1)then call utmsg('utacos') write(ifch,*)'***** argum = ',argum,' set -1' call utmsgf endif argum=-1. elseif(x.gt.1.)then if(ish.ge.1)then call utmsg('utacos') write(ifch,*)'***** argum = ',argum,' set 1' call utmsgf endif argum=1. endif utacos=acos(argum) return end c---------------------------------------------------------------------- function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus) c---------------------------------------------------------------------- c returns min mass of droplet with given u,d,s,c content c keux: net u quark number c kedx: net d quark number c kesx: net s quark number c kecx: net c quark number c kebx: net b quark number c ketx: net t quark number c modus: 4=two lowest multiplets; 5=lowest multiplet c---------------------------------------------------------------------- common/files/ifop,ifmt,ifch,ifcx,ifhi,ifdt,ifcp,ifdr common/csjcga/amnull,asuha(7) common/drop4/asuhax(7),asuhay(7) if(modus.lt.4.or.modus.gt.5)stop'UTAMNU: not supported' c 1 format(' flavours:',6i5 ) c 100 format(' flavours+mass:',6i5,f8.2 ) c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx amnull=0. do i=1,7 if(modus.eq.4)asuha(i)=asuhax(i) !two lowest multiplets if(modus.eq.5)asuha(i)=asuhay(i) !lowest multiplet enddo ke=iabs(keux+kedx+kesx+kecx+kebx+ketx) if(keux+kedx+kesx+kecx+kebx+ketx.ge.0)then keu=keux ked=kedx kes=kesx kec=kecx keb=kebx ket=ketx else keu=-keux ked=-kedx kes=-kesx kec=-kecx keb=-kebx ket=-ketx endif c write(ifch,*)keu,ked,kes,kec,keb,ket c removing top mesons to remove t quarks or antiquarks if(ket.ne.0)then 12 continue ii=sign(1,ket) ket=ket-ii if(ii*keu.le.ii*ked)then keu=keu+ii else ked=ked+ii endif amnull=amnull+200. ! ??????? if(ket.ne.0)goto12 endif c removing bottom mesons to remove b quarks or antiquarks if(keb.ne.0)then 11 continue ii=sign(1,keb) keb=keb-ii if(ii*keu.le.ii*ked)then keu=keu+ii else ked=ked+ii endif amnull=amnull+6. !5.28 ! (more than B-meson) if(keb.ne.0)goto11 endif c removing charm mesons to remove c quarks or antiquarks if(kec.ne.0)then 10 continue ii=sign(1,kec) kec=kec-ii if(keu*ii.le.ked*ii)then keu=keu+ii else ked=ked+ii endif amnull=amnull+2.2 !1.87 ! (more than D-meson) if(kec.ne.0)goto10 endif c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull c removing mesons to remove s antiquarks 5 continue if(kes.lt.0)then amnull=amnull+asuha(6) if(keu.ge.ked)then keu=keu-1 else ked=ked-1 endif kes=kes+1 goto5 endif c removing mesons to remove d antiquarks 6 continue if(ked.lt.0)then if(keu.ge.kes)then amnull=amnull+asuha(5) keu=keu-1 else amnull=amnull+asuha(6) kes=kes-1 endif ked=ked+1 goto6 endif c removing mesons to remove u antiquarks 7 continue if(keu.lt.0)then if(ked.ge.kes)then amnull=amnull+asuha(5) ked=ked-1 else amnull=amnull+asuha(6) kes=kes-1 endif keu=keu+1 goto7 endif c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull c print*,keu,ked,kes,kec,keb,ket,amnull if(keu+ked+kes+kec+keb+ket.ne.ke) *call utstop('utamnu: sum_kei /= ke&') keq=keu+ked keqx=keq amnux=0 c removing strange baryons i=4 2 i=i-1 3 continue if((4-i)*kes.gt.(i-1)*keq)then amnux=amnux+asuha(1+i) kes=kes-i keq=keq-3+i if(kes.lt.0)call utstop('utamnu: negative kes&') if(keq.lt.0)call utstop('utamnu: negative keq&') goto3 endif if(i.gt.1)goto2 if(keqx.gt.keq)then do 8 k=1,keqx-keq if(keu.ge.ked)then keu=keu-1 else ked=ked-1 endif 8 continue endif if(keu+ked.ne.keq)call utstop('utamnu: keu+ked /= keq&') c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux c print*,keu,ked,kes,kec,keb,ket,amnull+amnux c removing nonstrange baryons 9 continue if(keu.gt.2*ked)then amnux=amnux+asuha(7) keu=keu-3 if(keu.lt.0)call utstop('utamnu: negative keu&') goto9 endif if(ked.gt.2*keu)then amnux=amnux+asuha(7) ked=ked-3 if(ked.lt.0)call utstop('utamnu: negative ked&') goto9 endif keq=keu+ked c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux c print*,keu,ked,kes,kec,keb,ket,amnull+amnux if(mod(keq,3).ne.0)call utstop('utamnu: mod(keq,3) /= 0&') amnux=amnux+asuha(1)*keq/3 c write(ifch,100)keu,ked,kes,kec,keb,ket,amnull+amnux c print*,keu,ked,kes,kec,keb,ket,amnull+amnux amnull=amnull+amnux if(amnull.eq.0)amnull=asuha(5) utamnu=amnull return end c----------------------------------------------------------------------- function utamnx(jcp,jcm) c----------------------------------------------------------------------- c returns minimum mass for the decay of jcp---jcm (by calling utamnu). c----------------------------------------------------------------------- parameter (nflav=6) integer jcp(nflav,2),jcm(nflav,2) do i=1,nflav do j=1,2 if(jcp(i,j).ne.0)goto1 enddo enddo keu=jcm(1,1)-jcm(1,2) ked=jcm(2,1)-jcm(2,2) kes=jcm(3,1)-jcm(3,2) kec=jcm(4,1)-jcm(4,2) keb=jcm(5,1)-jcm(5,2) ket=jcm(6,1)-jcm(6,2) utamnx=utamnu(keu,ked,kes,kec,keb,ket,5) return 1 continue do i=1,nflav do j=1,2 if(jcm(i,j).ne.0)goto2 enddo enddo keu=jcp(1,1)-jcp(1,2) ked=jcp(2,1)-jcp(2,2) kes=jcp(3,1)-jcp(3,2) kec=jcp(4,1)-jcp(4,2) keb=jcp(5,1)-jcp(5,2) ket=jcp(6,1)-jcp(6,2) utamnx=utamnu(keu,ked,kes,kec,keb,ket,5) return 2 continue keu=jcp(1,1)-jcp(1,2) ked=jcp(2,1)-jcp(2,2) kes=jcp(3,1)-jcp(3,2) kec=jcp(4,1)-jcp(4,2) keb=jcp(5,1)-jcp(5,2) ket=jcp(6,1)-jcp(6,2) ke=keu+ked+kes+kec+keb+ket if(mod(ke+1,3).eq.0)then keu=keu+1 amms1=utamnu(keu,ked,kes,kec,keb,ket,5) keu=keu-1 ked=ked+1 amms2=utamnu(keu,ked,kes,kec,keb,ket,5) elseif(mod(ke-1,3).eq.0)then keu=keu-1 amms1=utamnu(keu,ked,kes,kec,keb,ket,5) keu=keu+1 ked=ked-1 amms2=utamnu(keu,ked,kes,kec,keb,ket,5) else amms1=0 amms2=0 amms3=0 amms4=0 call utstop('utamnx: no singlet possible (1)&') endif keu=jcm(1,1)-jcm(1,2) ked=jcm(2,1)-jcm(2,2) kes=jcm(3,1)-jcm(3,2) kec=jcm(4,1)-jcm(4,2) keb=jcm(5,1)-jcm(5,2) ket=jcm(6,1)-jcm(6,2) ke=keu+ked+kes+kec+keb+ket if(mod(ke+1,3).eq.0)then keu=keu+1 amms3=utamnu(keu,ked,kes,kec,keb,ket,5) keu=keu-1 ked=ked+1 amms4=utamnu(keu,ked,kes,kec,keb,ket,5) elseif(mod(ke-1,3).eq.0)then keu=keu-1 amms3=utamnu(keu,ked,kes,kec,keb,ket,5) keu=keu+1 ked=ked-1 amms4=utamnu(keu,ked,kes,kec,keb,ket,5) else call utstop('utamnx: no singlet possible (2)&') endif utamnx=min(amms1+amms3,amms2+amms4) c print *,amms1,amms3,amms2,amms4,jcp,jcm return end cc----------------------------------------------------------------------- c function utamny(jcp,jcm) cc----------------------------------------------------------------------- cc returns minimum mass of jcp+jcm (by calling utamnu). cc----------------------------------------------------------------------- c parameter (nflav=6) c integer jcp(nflav,2),jcm(nflav,2),jc(nflav,2) c do 7 nf=1,nflav c jc(nf,1)=jcp(nf,1)+jcm(nf,1) c7 jc(nf,2)=jcp(nf,2)+jcm(nf,2) c keu=jc(1,1)-jc(1,2) c ked=jc(2,1)-jc(2,2) c kes=jc(3,1)-jc(3,2) c kec=jc(4,1)-jc(4,2) c keb=jc(5,1)-jc(5,2) c ket=jc(6,1)-jc(6,2) c utamny=utamnu(keu,ked,kes,kec,keb,ket,5) c return c end c c----------------------------------------------------------------------- function utamnz(jc,modus) c----------------------------------------------------------------------- c returns minimum mass of jc (by calling utamnu). c----------------------------------------------------------------------- parameter (nflav=6) integer jc(nflav,2) keu=jc(1,1)-jc(1,2) ked=jc(2,1)-jc(2,2) kes=jc(3,1)-jc(3,2) kec=jc(4,1)-jc(4,2) keb=jc(5,1)-jc(5,2) ket=jc(6,1)-jc(6,2) utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus) return end c----------------------------------------------------------------------- subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx) c----------------------------------------------------------------------- c returns the array xx with xx(1)=x0 <= xx(i) <= xx(i3)=x3 c----------------------------------------------------------------------- real xx(i3) do 1 i=1,i1-1 1 xx(i)=x0+(i-1.)/(i1-1.)*(x1-x0) do 2 i=i1,i2-1 2 xx(i)=x1+(i-i1*1.)/(i2-i1*1.)*(x2-x1) do 3 i=i2,i3 3 xx(i)=x2+(i-i2*1.)/(i3-i2*1.)*(x3-x2) return end cc--------------------------------------------------------------------- c subroutine utaxis(i,j,a1,a2,a3) cc----------------------------------------------------------------------- cc calculates the axis defined by the ptls i,j in the i,j cm system cc--------------------------------------------------------------------- c include 'epos.inc' c double precision pi1,pi2,pi3,pi4,pj1,pj2,pj3,pj4,p1,p2,p3,p4,p5 c *,err,a c a1=0 c a2=0 c a3=1 c pi1=dble(pptl(1,i)) c pi2=dble(pptl(2,i)) c pi3=dble(pptl(3,i)) c pi4=dble(pptl(4,i)) c pj1=dble(pptl(1,j)) c pj2=dble(pptl(2,j)) c pj3=dble(pptl(3,j)) c pj4=dble(pptl(4,j)) c p1=pi1+pj1 c p2=pi2+pj2 c p3=pi3+pj3 c p4=pi4+pj4 c p5=dsqrt(p4**2-p3**2-p2**2-p1**2) c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50) c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51) c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2 c if(err.gt.1d-3)then c call utmsg('utaxis') c write(ifch,*)'***** err=',err c write(ifch,*)'pi:',pi1,pi2,pi3,pi4 c write(ifch,*)'pj:',pj1,pj2,pj3,pj4 c call utmsgf c endif c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 ) c if(a.eq.0.d0)return c a1=sngl((pi1-pj1)/a) c a2=sngl((pi2-pj2)/a) c a3=sngl((pi3-pj3)/a) c return c end c cc--------------------------------------------------------------------- c subroutine uthalf(i,j,zz1,zz2,iret) cc----------------------------------------------------------------------- cc give equal energy (E_i+E_j)/2 to particle i+j in their cm system cc--------------------------------------------------------------------- c include 'epos.inc' c double precision pi1,pi2,pi3,pi4,pi5,pj1,pj2,pj3,pj4,pj5 c *,p1,p2,p3,p4,p5,err,pt,pti,sinp,cosp,pmax,phi,drangen!,rrr c iret=0 c pi1=dble(pptl(1,i)) c pi2=dble(pptl(2,i)) c pi3=dble(pptl(3,i)) c pi5=dble(pptl(5,i)) c pi4=sqrt(pi1**2+pi2**2+pi3**2+pi5**2) c pj1=dble(pptl(1,j)) c pj2=dble(pptl(2,j)) c pj3=dble(pptl(3,j)) c pj5=dble(pptl(5,j)) c pj4=sqrt(pj1**2+pj2**2+pj3**2+pj5**2) c if(ish.ge.6)then c write(ifch,*)'uthalf for ',i,' and ',j c write(ifch,*)'in ',idptl(i),pi1,pi2,pi3,pi4,pi5 c write(ifch,*)'<> ',idptl(j),pj1,pj2,pj3,pj4,pj5 c endif c p1=pi1+pj1 c p2=pi2+pj2 c p3=pi3+pj3 c p4=pi4+pj4 c p5=(p4-p3)*(p4+p3)-p2**2-p1**2 c if(p5.lt.0d0.or.(pi3.lt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.4) c & .or.(pi3.gt.0.99d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then cc if(p5.lt.0d0)then c if(ish.ge.7)write(ifch,*)'Inversion not possible (1)',p5,pi3,pj3 c iret=1 c return c else c p5=sqrt(p5) c endif c call utlob2(1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,50) c call utlob2(1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,51) c err=(pi1+pj1)**2+(pi2+pj2)**2+(pi3+pj3)**2 c if(err.gt.1d-3)then c call utmsg('uthalf') c write(ifch,*)'***** err=',err c write(ifch,*)'pi:',pi1,pi2,pi3,pi4 c write(ifch,*)'pj:',pj1,pj2,pj3,pj4 c call utmsgf c endif c if(ish.ge.8)then c write(ifch,*)'pi 1:',pi1,pi2,pi3,pi4,pi5 c & ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2) c write(ifch,*)'pj 1:',pj1,pj2,pj3,pj4,pj5 c & ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2) c endif c c phi=drangen(p5)*2d0*dble(pi) c pti=sqrt(pi1*pi1+pi2*pi2) c cc sinp=abs(asin(pi2/pti)) c c cosp=cos(phi) c sinp=sin(phi) cc cosp=pi1/pti cc sinp=pi2/pti c pini=abs(pi3) c pmax=pini c c ntry=0 c 10 ntry=ntry+1 c pi1=-pj1 c pi2=-pj2 c pi3=-pj3 c c cc rrr=dble(ranptcut(4.))*max(0.1d0,dble(zz))+0.01d0 cc rrr=dble(ranptcut(zz*max(1.,float(ntry/10)**3))) !flip if pt too large cc rrr=dble(min(1.,ranptcut(zz/float(max(1,ntry/10))))) !return if pt too large c c rrr=rangen() c phi=dble(pi*zz1*(1.-rrr*zz2/float(max(1,ntry/10)))) c call utroa2(phi,cosp,sinp,0d0,pi1,pi2,pi3) !rotation around an axis perpendicular to p3 c pt=pi1*pi1+pi2*pi2 c ccc pi3=pini*(1d0-(pmax/pini+1d0)*min(1d0,rrr)) ccc & *sign(1.d0,pi3) cc pi3=pini*(1d0-(pmax/pini+1d0)*(1d0-min(1d0,rrr))) cc & *sign(1.d0,pi3) cccc pi3=(1d0-exp(-drangen(pti)))*pi3 cccc pi3=min(1d0,(drangen()+exp(-min(100d0,dble(zz)))))*pi3 cccc pi3=(dble(reminv)+exp(-min(100d0,dble(zz))))*pi3 cc pt=((pi4+pi3)*(pi4-pi3)-pi5*pi5) c if(ish.ge.8)write(ifch,*)'ut',ntry,zz,rrr,-pi3/pj3,pt,pti*pti c &,idptl(i),idptl(j) c if((pt.lt.0d0.or.pt.gt.max(2.*pti*pti,1.d0)) cc if(pt.lt.0d0 c & .and.ntry.lt.1000)then c goto 10 c elseif(ntry.ge.1000)then c if(ish.ge.7)write(ifch,*)'Inversion not possible (2)',pt,pti*pti c iret=1 c return !pion distribution fall at xf=1 cc pi3=pj3 !if flip all particle with large pt cc pt=pti*pti !then very very hard pi+ spectra (flat at xf=1 like in pp400 !) c !but to hard for K, rho, eta, etc .. c endif cc print *,'ut',ntry,phi/(pi-sinp),zz,rrr,-pi3/pj3,pt/pti/pti cc pt=sqrt(pt) cc pi1=cosp*pt cc pi2=sinp*pt c c pj1=-pi1 c pj2=-pi2 c pj3=-pi3 cc pi4=sqrt(pi3*pi3+pt*pt+pi5*pi5) cc pj4=sqrt(pj3*pj3+pt*pt+pj5*pj5) c pi4=sqrt(pi3*pi3+pt+pi5*pi5) c pj4=sqrt(pj3*pj3+pt+pj5*pj5) c if(ish.ge.8)then c write(ifch,*)'pi 2:',pi1,pi2,pi3,pi4,pi5 c & ,sqrt(pi1**2+pi2**2+pi3**2+pi5**2) c write(ifch,*)'pj 2:',pj1,pj2,pj3,pj4,pj5 c & ,sqrt(pj1**2+pj2**2+pj3**2+pj5**2) c endif c call utlob2(-1,p1,p2,p3,p4,p5,pi1,pi2,pi3,pi4,-50) c call utlob2(-1,p1,p2,p3,p4,p5,pj1,pj2,pj3,pj4,-51) c if(pi3/dble(pptl(3,i)).gt.1.00001d0 c &.or.(pi3.gt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.4) c &.or.(pi3.lt.1.001d0*pj3.and.mod(ityptl(j)/10,10).eq.5))then c if(ish.ge.7)write(ifch,*)'Inversion not possible (3)',pi3,pj3 c iret=1 c return c endif c id=idptl(i) c pptl(1,i)=sngl(pi1) c pptl(2,i)=sngl(pi2) c pptl(3,i)=sngl(pi3) c pptl(4,i)=sngl(pi4) c pptl(5,i)=sngl(pi5) cc pptl(5,i)=sngl(pj5) cc idptl(i)=idptl(j) c pptl(1,j)=sngl(pj1) c pptl(2,j)=sngl(pj2) c pptl(3,j)=sngl(pj3) c pptl(4,j)=sngl(pj4) c pptl(5,j)=sngl(pj5) cc pptl(5,j)=sngl(pi5) cc idptl(j)=id c if(ish.ge.6)then c write(ifch,*)'out ',idptl(i),pi1,pi2,pi3,pi4 c write(ifch,*)' <> ',idptl(j),pj1,pj2,pj3,pj4 c endif c c return c end c cc----------------------------------------------------------------------- c subroutine utchm(arp,arm,ii) cc----------------------------------------------------------------------- cc checks whether arp**2=0 and arm**2=0. cc----------------------------------------------------------------------- c include 'epos.inc' c double precision arp(4),arm(4),difp,difm c difp=arp(4)**2-arp(1)**2-arp(2)**2-arp(3)**2 c difm=arm(4)**2-arm(1)**2-arm(2)**2-arm(3)**2 c if(dabs(difp).gt.1e-3*arp(4)**2 c *.or.dabs(difm).gt.1e-3*arm(4)**2)then c call utmsg('utchm ') c write(ifch,*)'***** mass non zero - ',ii c write(ifch,*)'jet-mass**2`s: ',difp,difm c write(ifch,*)'energy**2`s: ',arp(4)**2,arm(4)**2 c write(ifch,*)(sngl(arp(i)),i=1,4) c write(ifch,*)(sngl(arm(i)),i=1,4) c call utmsgf c endif c return c end c c----------------------------------------------------------------------- subroutine utclea(nptlii,nptl0) c----------------------------------------------------------------------- c starting from nptlii c overwrites istptl=99 particles in /cptl/, reduces so nptl c and update minfra and maxfra c----------------------------------------------------------------------- include 'epos.inc' integer newptl(mxptl)!,oldptl(mxptl),ii(mxptl) ish0=ish if(ishsub/100.eq.18)ish=mod(ishsub,100) call utpri('utclea',ish,ishini,2) nptli=max(maproj+matarg+1,nptlii) minfra0=minfra maxfra0=maxfra minfra1=maxfra maxfra1=minfra if(ish.ge.2)write(ifch,*)'entering subr utclea:',nptl & ,minfra,maxfra if(ish.ge.7)then write(ifch,*)('-',l=1,68) write(ifch,*)'sr utclea. initial.' write(ifch,*)('-',l=1,68) do 34 n=nptli,nptl write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n) *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n) *,istptl(n),ityptl(n) 34 continue 116 format(1x,i6,i6,4x,i6,4x,i6,i6,i12,3x,3(e8.2,1x),i3,i3) endif c ish=ish0 c ish0=ish c if(ishsub/100.eq.18)ish=mod(ishsub,100) i=nptli-1 1 i=i+1 if(i.gt.nptl)goto 1000 if(istptl(i).eq.99)goto 2 newptl(i)=i c oldptl(i)=i goto 1 2 i=i-1 j=i 3 i=i+1 4 j=j+1 if(j.gt.nptl)goto 5 newptl(j)=0 if(istptl(j).eq.99)goto 4 newptl(j)=i c oldptl(i)=j c write(ifch,*)'move',j,' to ',i c write(ifch,*)idptl(i),ityptl(i),idptl(j),ityptl(j),minfra,maxfra call utrepl(i,j) if(j.ge.minfra0.and.j.le.maxfra0)then minfra1=min(minfra1,i) maxfra1=max(maxfra1,i) endif goto 3 5 nptl=i-1 if(nptl.eq.0)then nptl0=0 goto 1000 endif 20 n0=newptl(nptl0) if(n0.gt.0)then nptl0=n0 else nptl0=nptl0-1 if(nptl0.gt.0)goto 20 endif c do 11 k=1,nptl c io=iorptl(k) c if(io.le.0)ii(k)=io c if(io.gt.0)ii(k)=newptl(io) c11 continue c do 12 k=1,nptl c12 iorptl(k)=ii(k) c c do 13 k=1,nptl c jo=jorptl(k) c if(jo.le.0)ii(k)=jo c if(jo.gt.0)ii(k)=newptl(jo) c13 continue c do 14 k=1,nptl c14 jorptl(k)=ii(k) c c do 15 k=1,nptl c if1=ifrptl(1,k) c if(if1.le.0)ii(k)=if1 c if(if1.gt.0)ii(k)=newptl(if1) c15 continue c do 16 k=1,nptl c16 ifrptl(1,k)=ii(k) c c do 17 k=1,nptl c if2=ifrptl(2,k) c if(if2.le.0)ii(k)=if2 c if(if2.gt.0)ii(k)=newptl(if2) c17 continue c do 18 k=1,nptl c18 ifrptl(2,k)=ii(k) c c do 19 k=1,nptl c if(ifrptl(1,k).eq.0.and.ifrptl(2,k).gt.0)ifrptl(1,k)=ifrptl(2,k) c if(ifrptl(2,k).eq.0.and.ifrptl(1,k).gt.0)ifrptl(2,k)=ifrptl(1,k) c19 continue 1000 continue if(minfra1.lt.minfra0)minfra=minfra1 if(maxfra1.ge.minfra1)maxfra=maxfra1 if(ish.ge.2)then write(ifch,*)'exiting subr utclea:' do 35 n=1,nptl write(ifch,116)iorptl(n),jorptl(n),n,ifrptl(1,n),ifrptl(2,n) *,idptl(n),sqrt(pptl(1,n)**2+pptl(2,n)**2),pptl(3,n),pptl(5,n) *,istptl(n),ityptl(n) 35 continue endif if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl & ,minfra,maxfra call utprix('utclea',ish,ishini,2) ish=ish0 return end c--------------------------------------------------------------------- subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q) c--------------------------------------------------------------------- c linear fit to data c input: c ndata: nr of data points c x(),y(),sig(): data c mwt: unweighted (0) or weighted (else) data points c output: c a,b: parameters of linear fit a+b*x c--------------------------------------------------------------------- INTEGER mwt,ndata REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata) CU USES utgmq INTEGER i REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq sx=0. sy=0. st2=0. b=0. if(mwt.ne.0) then ss=0. do 11 i=1,ndata wt=1./(sig(i)**2) ss=ss+wt sx=sx+x(i)*wt sy=sy+y(i)*wt 11 continue else do 12 i=1,ndata sx=sx+x(i) sy=sy+y(i) 12 continue ss=float(ndata) endif sxoss=sx/ss if(mwt.ne.0) then do 13 i=1,ndata t=(x(i)-sxoss)/sig(i) st2=st2+t*t b=b+t*y(i)/sig(i) 13 continue else do 14 i=1,ndata t=x(i)-sxoss st2=st2+t*t b=b+t*y(i) 14 continue endif b=b/st2 a=(sy-sx*b)/ss siga=sqrt((1.+sx*sx/(ss*st2))/ss) sigb=sqrt(1./st2) chi2=0. if(mwt.eq.0) then do 15 i=1,ndata chi2=chi2+(y(i)-a-b*x(i))**2 15 continue q=1. sigdat=sqrt(chi2/(ndata-2)) siga=siga*sigdat sigb=sigb*sigdat else do 16 i=1,ndata chi2=chi2+((y(i)-a-b*x(i))/sig(i))**2 16 continue q=utgmq(0.5*(ndata-2),0.5*chi2) endif return END c----------------------------------------------------------------------- function utgam1(x) c----------------------------------------------------------------------- c gamma fctn tabulated c single precision c----------------------------------------------------------------------- double precision utgamtab,utgam,al,dl common/gamtab/utgamtab(10000) if(x.gt.0.01.and.x.lt.99.99)then al=100.d0*dble(x) k1=int(al) k2=k1+1 dl =al-dble(k1) utgam1=real(utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl)) elseif(x.eq.0.)then utgam1=0. else utgam1=real(utgam(dble(x))) endif end c----------------------------------------------------------------------- double precision function utgam2(x) c----------------------------------------------------------------------- c gamma fctn tabulated c double precision c----------------------------------------------------------------------- double precision utgamtab,x,al,dl,utgam common/gamtab/utgamtab(10000) if(x.gt.0.01d0.and.x.le.99.99d0)then al=100.d0*x k1=int(al) k2=k1+1 dl =al-dble(k1) utgam2=utgamtab(k2)*dl+utgamtab(k1)*(1.d0-dl) elseif(x.eq.0.d0)then utgam2=0.d0 else utgam2=utgam(x) endif end c----------------------------------------------------------------------- double precision function utgam(x) c----------------------------------------------------------------------- c gamma fctn c double precision c----------------------------------------------------------------------- include 'epos.inc' double precision c(13),x,z,f data c 1/ 0.00053 96989 58808, 0.00261 93072 82746, 0.02044 96308 23590, 2 0.07309 48364 14370, 0.27964 36915 78538, 0.55338 76923 85769, 3 0.99999 99999 99998,-0.00083 27247 08684, 0.00469 86580 79622, 4 0.02252 38347 47260,-0.17044 79328 74746,-0.05681 03350 86194, 5 1.13060 33572 86556/ utgam=0d0 z=x if(x .gt. 170.d0) goto6 if(x .gt. 0.0d0) goto1 if(x .eq. int(x)) goto5 z=1.0d0-z 1 f=1.0d0/z if(z .le. 1.0d0) goto4 f=1.0d0 2 continue if(z .lt. 2.0d0) goto3 z=z-1.0d0 f=f*z goto2 3 z=z-1.0d0 4 utgam= 1 f*((((((c(1)*z+c(2))*z+c(3))*z+c(4))*z+c(5))*z+c(6))*z+c(7))/ 2 ((((((c(8)*z+c(9))*z+c(10))*z+c(11))*z+c(12))*z+c(13))*z+1.0d0) if(x .gt. 0.0d0) return utgam=3.141592653589793d0/(sin(3.141592653589793d0*x)*utgam) return 5 write(ifch,10)sngl(x) 10 format(1x,'argument of gamma fctn = ',e20.5) call utstop('utgam : negative integer argument&') 6 write(ifch,11)sngl(x) 11 format(1x,'argument of gamma fctn = ',e20.5) call utstop('utgam : argument too large&') end c--------------------------------------------------------------------- subroutine utgcf(gammcf,a,x,gln) c--------------------------------------------------------------------- INTEGER ITMAX REAL a,gammcf,gln,x,EPS,FPMIN PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30) CU USES utgmln INTEGER i REAL an,b,c,d,del,h,utgmln gln=utgmln(a) b=x+1.-a c=1./FPMIN d=1./b h=d do 11 i=1,ITMAX an=-i*(i-a) b=b+2. d=an*d+b if(abs(d).lt.FPMIN)d=FPMIN c=b+an/c if(abs(c).lt.FPMIN)c=FPMIN d=1./d del=d*c h=h*del if(abs(del-1.).lt.EPS)goto 1 11 continue call utstop("a too large, ITMAX too small in utgcf&") 1 gammcf=exp(-x+a*log(x)-gln)*h return END c--------------------------------------------------------------------- function utgmln(xx) c--------------------------------------------------------------------- REAL utgmln,xx INTEGER j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, *-.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue utgmln=tmp+log(stp*ser/x) return END c--------------------------------------------------------------------- function utgmq(a,x) c--------------------------------------------------------------------- REAL a,utgmq,x CU USES utgcf,utgser REAL gammcf,gamser,gln if(x.lt.0..or.a.le.0.) call utstop("bad arguments in utgmq&") if(x.lt.a+1.)then call utgser(gamser,a,x,gln) utgmq=1.-gamser else call utgcf(gammcf,a,x,gln) utgmq=gammcf endif return END c--------------------------------------------------------------------- subroutine utgser(gamser,a,x,gln) c--------------------------------------------------------------------- INTEGER ITMAX REAL a,gamser,gln,x,EPS PARAMETER (ITMAX=100,EPS=3.e-7) CU USES utgmln INTEGER n REAL ap,del,sum,utgmln gln=utgmln(a) if(x.le.0.)then if(x.lt.0.)call utstop("x < 0 in utgser&") gamser=0. return endif ap=a sum=1./a del=sum do 11 n=1,ITMAX ap=ap+1. del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*EPS)goto 1 11 continue call utstop("a too large, ITMAX too small in utgser&") 1 gamser=sum*exp(-x+a*log(x)-gln) return END c------------------------------------------------------------------------- subroutine uticpl(ic,ifla,iqaq,iret) c------------------------------------------------------------------------- c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla c to 2-id ic c------------------------------------------------------------------------- include 'epos.inc' integer jc(nflav,2),ic(2) iret=0 if(ifla.eq.0)return call iddeco(ic,jc) if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc jqaq=3-iqaq if(jc(ifla,jqaq).gt.0)then jc(ifla,jqaq)=jc(ifla,jqaq)-1 else jc(ifla,iqaq)=jc(ifla,iqaq)+1 endif call idcomj(jc) call idenco(jc,ic,ireten) if(ish.ge.8)write(ifch,'(2i8,12i3)')ic,jc if(ireten.eq.1)iret=1 if(ic(1).eq.0.and.ic(2).eq.0.and.ireten.eq.0)then ic(1)=100000 ic(2)=100000 endif return end cc----------------------------------------------------------------------- c subroutine utindx(n,xar,x,i) cc----------------------------------------------------------------------- cc input: dimension n cc array xar(n) with xar(i) > xar(i-1) cc some number x between xar(1) and xar(n) cc output: the index i such that x is between xar(i) and xar(i+1) cc----------------------------------------------------------------------- c include 'epos.inc' c real xar(n) c if(x.lt.xar(1))then c if(ish.ge.5)then c call utmsg('utindx') c write(ifch,*)'***** x=',x,' < xar(1)=',xar(1) c call utmsgf c endif c i=1 c return c elseif(x.gt.xar(n))then c if(ish.ge.5)then c call utmsg('utindx') c write(ifch,*)'***** x=',x,' > xar(n)=',xar(n) c call utmsgf c endif c i=n c return c endif c lu=1 c lo=n c1 lz=(lo+lu)/2 c if((xar(lu).le.x).and.(x.le.xar(lz)))then c lo=lz c elseif((xar(lz).lt.x).and.(x.le.xar(lo)))then c lu=lz c else c call utstop('utindx: no interval found&') c endif c if((lo-lu).ge.2) goto1 c if(lo.le.lu)call utstop('utinvt: lo.le.lu&') c i=lu c return c end c c----------------------------------------------------------------------- function utinvt(n,x,q,y) c----------------------------------------------------------------------- c returns x with y=q(x) c----------------------------------------------------------------------- include 'epos.inc' real x(n),q(n) if(q(n).eq.0.)call utstop('utinvt: q(n)=0&') if(y.lt.0.)then if(ish.ge.1)then call utmsg('utinvt') write(ifch,*)'***** y=',y,' < 0' call utmsgf endif y=0. elseif(y.gt.q(n))then if(ish.ge.1)then call utmsg('utinvt') write(ifch,*)'***** y=',y,' > ',q(n) call utmsgf endif y=q(n) endif lu=1 lo=n 1 lz=(lo+lu)/2 if((q(lu).le.y).and.(y.le.q(lz)))then lo=lz elseif((q(lz).lt.y).and.(y.le.q(lo)))then lu=lz else write(ifch,*)'q(1),y,q(n):',q(1),y,q(n) write(ifch,*)'lu,lz,lo:',lu,lz,lo write(ifch,*)'q(lu),q(lz),q(lo):',q(lu),q(lz),q(lo) call utstop('utinvt: no interval found&') endif if((lo-lu).ge.2) goto1 if(lo.le.lu)call utstop('utinvt: lo.le.lu&') utinvt=x(lu)+(y-q(lu))*(x(lo)-x(lu))/(q(lo)-q(lu)) return end c----------------------------------------------------------------------- subroutine utlob2(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4,idi) c----------------------------------------------------------------------- c performs a lorentz boost, double prec. c isig=+1 is to boost the four vector x1,x2,x3,x4 such as to obtain it c in the frame specified by the 5-vector p1...p5 (5-vector=4-vector+mass). c isig=-1: the other way round, that means, c if the 4-vector x1...x4 is given in some frame characterized by c p1...p5 with respect to to some lab-frame, utlob2 returns the 4-vector c x1...x4 in the lab frame. c idi is a call identifyer (integer) to identify the call in case of problem c----------------------------------------------------------------------- include 'epos.inc' double precision beta(4),z(4),p1,p2,p3,p4,p5,pp,bp,x1,x2,x3,x4 *,xx0,x10,x20,x30,x40,x4x,x0123 if(ish.ge.2)then if(ish.ge.9)then write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1 write(ifch,301)p1,p2,p3,p4,p5,(p4-p3)*(p4+p3)-p2*p2-p1*p1 101 format(' utlob2: x = ',5e13.5) 301 format(' p = ',6e13.5) endif pp=(p4-p3)*(p4+p3)-p2*p2-p1*p1 if(dabs(pp-p5*p5).gt.1e-3*p4*p4.and.dabs(pp-p5*p5).gt.1e-3)then call utmsg('utlob2') write(ifch,*)'***** p**2 .ne. p5**2' write(ifch,*)'call identifyer:',idi write(ifch,*)'p**2,p5**2: ',pp,p5*p5 write(ifch,*)'p: ',p1,p2,p3,p4,p5 call utmsgf endif x10=x1 x20=x2 x30=x3 x40=x4 endif xx0=(x4-x3)*(x4+x3)-x2*x2-x1*x1 if(p5.le.0.)then call utmsg('utlob2') write(ifch,*)'***** p5 negative.' write(ifch,*)'call identifyer:',idi write(ifch,*)'p(5): ',p1,p2,p3,p4,p5 write(ifmt,*)'call identifyer:',idi write(ifmt,*)'p(5): ',p1,p2,p3,p4,p5 call utmsgf call utstop('utlob2: p5 negative.&') endif z(1)=x1 z(2)=x2 z(3)=x3 z(4)=x4 beta(1)=-p1/p5 beta(2)=-p2/p5 beta(3)=-p3/p5 beta(4)= p4/p5 bp=0. do 220 k=1,3 220 bp=bp+z(k)*isig*beta(k) do 230 k=1,3 230 z(k)=z(k)+isig*beta(k)*z(4) *+isig*beta(k)*bp/(beta(4)+1.) z(4)=beta(4)*z(4)+bp x1=z(1) x2=z(2) x3=z(3) x4=z(4) if(ish.ge.9) *write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1 x4x=x4 x0123=xx0+x1*x1+x2*x2+x3*x3 if(x0123.gt.0.)then x4=sign( dsqrt(x0123) , x4x ) else x4=0 endif if(ish.ge.9)then write(ifch,101)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1 endif if(ish.ge.2)then if(ish.ge.9)write(ifch,*)'check x**2_ini -- x**2_fin' if(dabs(x4-x4x).gt.1d-2*dabs(x4).and.dabs(x4-x4x).gt.1d-2)then call utmsg('utlob2') write(ifch,*)'***** x**2_ini .ne. x**2_fin.' write(ifch,*)'call identifyer:',idi write(ifch,*)'x1 x2 x3 x4 x**2 (initial/final/corrected):' 102 format(5e13.5) write(ifch,102)x10,x20,x30,x40,(x40-x30)*(x40+x30)-x20*x20-x10*x10 write(ifch,102)x1,x2,x3,x4x,(x4x-x3)*(x4x+x3)-x2*x2-x1*x1 write(ifch,102)x1,x2,x3,x4,(x4-x3)*(x4+x3)-x2*x2-x1*x1 call utmsgf endif endif if(ish.ge.9)write(ifch,*)'return from utlob2' return end c----------------------------------------------------------------------- subroutine utlob3(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4) c----------------------------------------------------------------------- c performs a lorentz boost, double prec. c but arguments are single precision c----------------------------------------------------------------------- double precision xx1,xx2,xx3,xx4 xx1=dble(x1) xx2=dble(x2) xx3=dble(x3) xx4=dble(x4) call utlob2(isig *,dble(p1),dble(p2),dble(p3),dble(p4),dble(p5) *,xx1,xx2,xx3,xx4,52) x1=sngl(xx1) x2=sngl(xx2) x3=sngl(xx3) x4=sngl(xx4) return end c----------------------------------------------------------------------- subroutine utlob5(yboost,x1,x2,x3,x4,x5) c----------------------------------------------------------------------- amt=sqrt(x5**2+x1**2+x2**2) y=sign(1.,x3)*alog((x4+abs(x3))/amt) y=y-yboost x4=amt*cosh(y) x3=amt*sinh(y) return end c----------------------------------------------------------------------- subroutine utlob4(isig,pp1,pp2,pp3,pp4,pp5,x1,x2,x3,x4) c----------------------------------------------------------------------- c performs a lorentz boost, double prec. c but arguments are partly single precision c----------------------------------------------------------------------- double precision xx1,xx2,xx3,xx4,pp1,pp2,pp3,pp4,pp5 xx1=dble(x1) xx2=dble(x2) xx3=dble(x3) xx4=dble(x4) call utlob2(isig,pp1,pp2,pp3,pp4,pp5,xx1,xx2,xx3,xx4,53) x1=sngl(xx1) x2=sngl(xx2) x3=sngl(xx3) x4=sngl(xx4) return end c----------------------------------------------------------------------- subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4) c----------------------------------------------------------------------- c performs a lorentz boost c----------------------------------------------------------------------- include 'epos.inc' real beta(4),z(4) if(p5.le.0.)then call utmsg('utlobo') write(ifch,*)'***** mass <= 0.' write(ifch,*)'p(5): ',p1,p2,p3,p4,p5 call utmsgf call utstop('utlobo: mass <= 0.&') endif z(1)=x1 z(2)=x2 z(3)=x3 z(4)=x4 beta(1)=-p1/p5 beta(2)=-p2/p5 beta(3)=-p3/p5 beta(4)= p4/p5 bp=0. do 220 k=1,3 220 bp=bp+z(k)*isig*beta(k) do 230 k=1,3 230 z(k)=z(k)+isig*beta(k)*z(4) *+isig*beta(k)*bp/(beta(4)+1.) z(4)=beta(4)*z(4)+bp x1=z(1) x2=z(2) x3=z(3) x4=z(4) return end c----------------------------------------------------------------------- subroutine utloc(ar,n,a,l) c----------------------------------------------------------------------- real ar(n) do 1 i=1,n l=i-1 if(a.lt.ar(i))return 1 continue l=n return end cc----------------------------------------------------------------------- c subroutine utlow(cone) cc----------------------------------------------------------------------- c character*1 cone c if(cone.eq.'A')cone='a' c if(cone.eq.'B')cone='b' c if(cone.eq.'C')cone='c' c if(cone.eq.'D')cone='d' c if(cone.eq.'E')cone='e' c if(cone.eq.'F')cone='f' c if(cone.eq.'G')cone='g' c if(cone.eq.'H')cone='h' c if(cone.eq.'I')cone='i' c if(cone.eq.'J')cone='j' c if(cone.eq.'K')cone='k' c if(cone.eq.'L')cone='l' c if(cone.eq.'M')cone='m' c if(cone.eq.'N')cone='n' c if(cone.eq.'O')cone='o' c if(cone.eq.'P')cone='p' c if(cone.eq.'Q')cone='q' c if(cone.eq.'R')cone='r' c if(cone.eq.'S')cone='s' c if(cone.eq.'T')cone='t' c if(cone.eq.'U')cone='u' c if(cone.eq.'V')cone='v' c if(cone.eq.'W')cone='w' c if(cone.eq.'X')cone='x' c if(cone.eq.'Y')cone='y' c if(cone.eq.'Z')cone='z' c return c end c cc----------------------------------------------------------------------- c subroutine utlow3(cthree) cc----------------------------------------------------------------------- c character cthree*3 c do 1 i=1,3 c1 call utlow(cthree(i:i)) c return c end c cc----------------------------------------------------------------------- c subroutine utlow6(csix) cc----------------------------------------------------------------------- c character csix*6 c do 1 i=1,6 c1 call utlow(csix(i:i)) c return c end c cc----------------------------------------------------------------------- c function utmom(k,n,x,q) cc----------------------------------------------------------------------- cc calculates kth moment for f(x) with q(i)=int[0,x(i)]f(z)dz cc----------------------------------------------------------------------- c real x(n),q(n) c if(n.lt.2)call utstop('utmom : dimension too small&') c utmom=0 c do 1 i=2,n c1 utmom=utmom+((x(i)+x(i-1))/2)**k*(q(i)-q(i-1)) c utmom=utmom/q(n) c return c end c c----------------------------------------------------------------------- function utpcm(a,b,c) c----------------------------------------------------------------------- c calculates cm momentum for a-->b+c c----------------------------------------------------------------------- val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c) if(val.lt.0..and.val.gt.-1e-4)then utpcm=0 return endif utpcm=sqrt(val)/(2.*a) return end c----------------------------------------------------------------------- double precision function utpcmd(a,b,c,iret) c----------------------------------------------------------------------- c calculates cm momentum for a-->b+c c----------------------------------------------------------------------- double precision a,b,c,val iret=0 val=(a*a-b*b-c*c)*(a*a-b*b-c*c)-(2.*b*c)*(2.*b*c) utpcmd=0d0 if(val.lt.0d0.and.val.gt.-1d-4)then return elseif(val.lt.0d0)then iret=1 return endif utpcmd=sqrt(val)/(2.d0*a) return end c----------------------------------------------------------------------- subroutine utpri(text,ishi,ishini,ishx) c----------------------------------------------------------------------- include 'epos.inc' character*(*) text c double precision seedx !!! ishini=ishi if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return if(nrpri.gt.0)then do nr=1,nrpri if(subpri(nr)(1:6).eq.text)then ishi=ishpri(nr) endif enddo endif if(ish.ge.ishx)then write(ifch,'(1x,43a)') * ('-',i=1,10),' entry ',text,' ',('-',i=1,30) c call ranfgt(seedx) !!! c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!! endif return end c----------------------------------------------------------------------- subroutine utprix(text,ishi,ishini,ishx) c----------------------------------------------------------------------- include 'epos.inc' character*(*) text if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return if(ish.ge.ishx)write(ifch,'(1x,44a)') *('-',i=1,30),' exit ',text,' ',('-',i=1,11) ishi=ishini return end c----------------------------------------------------------------------- subroutine utprj(text,ishi,ishini,ishx) c----------------------------------------------------------------------- include 'epos.inc' character*(*) text c double precision seedx !!! idx=index(text,' ')-1 ishini=ishi if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return if(nrpri.gt.0)then do nr=1,nrpri if(subpri(nr)(1:idx).eq.text(1:idx))then ishi=ishpri(nr) endif enddo endif if(ish.ge.ishx)then write(ifch,'(1x,43a)') * ('-',i=1,10),' entry ',text(1:idx),' ',('-',i=1,30) c call ranfgt(seedx) !!! c if(ish.ge.ishx)write(ifch,*)'seed:',seedx !!! endif return end c----------------------------------------------------------------------- subroutine utprjx(text,ishi,ishini,ishx) c----------------------------------------------------------------------- include 'epos.inc' character*(*) text idx=index(text,' ')-1 if(ishevt.ne.0.and.nrevt+1.ne.ishevt)return if(ish.ge.ishx)write(ifch,'(1x,44a)') *('-',i=1,30),' exit ',text(1:idx),' ',('-',i=1,11) ishi=ishini return end c----------------------------------------------------------------------- function utquad(m,x,f,k) c----------------------------------------------------------------------- c performs an integration according to simpson c----------------------------------------------------------------------- real x(m),f(m) utquad=0 do 1 i=1,k-1 1 utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i)) return end c----------------------------------------------------------------------- subroutine utquaf(fu,n,x,q,x0,x1,x2,x3) c----------------------------------------------------------------------- c returns q(i) = integral [x(1)->x(i)] fu(x) dx c----------------------------------------------------------------------- include 'epos.inc' real x(n),q(n) parameter (m=10) real xa(m),fa(m) external fu if(x1.lt.x0.or.x2.lt.x1.or.x3.lt.x2)then if(ish.ge.1)then call utmsg('utquaf') write(ifch,*)' xi=',x0,x1,x2,x3 call utmsgf endif endif call utar(n/3,n*2/3,n,x0,x1,x2,x3,x) q(1)=0 do 2 i=2,n do 3 k=1,m z=x(i-1)+(k-1.)/(m-1.)*(x(i)-x(i-1)) xa(k)=z 3 fa(k)=fu(z) q(i)=q(i-1)+utquad(m,xa,fa,m) 2 continue return end c----------------------------------------------------------------------- subroutine utrepl(i,j) c----------------------------------------------------------------------- c i is replaced by j in /cptl/ c----------------------------------------------------------------------- include 'epos.inc' do 1 k=1,5 1 pptl(k,i) =pptl(k,j) iorptl(i) = 0 !iorptl(j) idptl(i) =idptl(j) istptl(i) =istptl(j) do 2 k=1,2 2 tivptl(k,i)=tivptl(k,j) do 3 k=1,2 3 ifrptl(k,i)= 0 !ifrptl(k,j) jorptl(i) = 0 !jorptl(j) do 4 k=1,4 4 xorptl(k,i)=xorptl(k,j) do 5 k=1,4 5 ibptl(k,i) =ibptl(k,j) ityptl(i) =ityptl(j) iaaptl(i) =iaaptl(j) radptl(i) =radptl(j) desptl(i) =desptl(j) dezptl(i) =dezptl(j) qsqptl(i) =qsqptl(j) zpaptl(1,i)=zpaptl(1,j) zpaptl(2,i)=zpaptl(2,j) itsptl(i) =itsptl(j) rinptl(i) =rinptl(j) return end c----------------------------------------------------------------------- subroutine utrepla(i,j) c----------------------------------------------------------------------- c i is replaced by j in /cptl/ c----------------------------------------------------------------------- include 'epos.inc' do 1 k=1,5 1 pptl(k,i) =pptl(k,j) iorptl(i) = iorptl(j) idptl(i) =idptl(j) istptl(i) =istptl(j) do 2 k=1,2 2 tivptl(k,i)=tivptl(k,j) do 3 k=1,2 3 ifrptl(k,i)= ifrptl(k,j) jorptl(i) = jorptl(j) do 4 k=1,4 4 xorptl(k,i)=xorptl(k,j) do 5 k=1,4 5 ibptl(k,i) =ibptl(k,j) ityptl(i) =ityptl(j) iaaptl(i) =iaaptl(j) radptl(i) =radptl(j) desptl(i) =desptl(j) dezptl(i) =dezptl(j) qsqptl(i) =qsqptl(j) zpaptl(1,i)=zpaptl(1,j) zpaptl(2,i)=zpaptl(2,j) itsptl(i) =itsptl(j) rinptl(i) =rinptl(j) return end cc----------------------------------------------------------------------- c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten) cc----------------------------------------------------------------------- c parameter (nflav=6) c integer icm(2),icp(2),jcm(nflav,2),jcp(nflav,2) c icm(1)=icm1 c icm(2)=icm2 c icp(1)=icp1 c icp(2)=icp2 c CALL IDDECO(ICM,JCM) c CALL IDDECO(ICP,JCP) c do 37 nf=1,nflav c do 37 k=1,2 c37 jcP(nf,k)=jcp(nf,k)+jcm(nf,k) c CALL IDENCO(JCP,ICP,IRETEN) c IDP=IDTRA(ICP,0,0,3) c call idres(idp,amp,idpr,iadj) c return c end c cc----------------------------------------------------------------------- c subroutine utroa1(phi,a1,a2,a3,x1,x2,x3) cc----------------------------------------------------------------------- cc rotates x by angle phi around axis a. cc normalization of a is irrelevant. cc----------------------------------------------------------------------- c double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi c dphi=phi c xx(1)=x1 c xx(2)=x2 c xx(3)=x3 c aa(1)=a1 c aa(2)=a2 c aa(3)=a3 c aaa=0 c xxx=0 c do i=1,3 c aaa=aaa+aa(i)**2 c xxx=xxx+xx(i)**2 c enddo c if(xxx.eq.0d0)return c if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&') c aaa=dsqrt(aaa) c xxx=dsqrt(xxx) cc e3 = a / !a! c do i=1,3 c e3(i)=aa(i)/aaa c enddo cc x_parallel c xp=0 c do i=1,3 c xp=xp+xx(i)*e3(i) c enddo cc x_transverse c if(xxx**2-xp**2.le.0.)return c xt=dsqrt(xxx**2-xp**2) cc e1 = vector x_transverse / absolute value x_transverse c do i=1,3 c e1(i)=(xx(i)-e3(i)*xp)/xt c enddo cc e2 orthogonal e3,e1 c call utvec2(e3,e1,e2) cc rotate x c do i=1,3 c xx(i)=xp*e3(i)+xt*dcos(dphi)*e1(i)+xt*dsin(dphi)*e2(i) c enddo cc back to single precision c x1=xx(1) c x2=xx(2) c x3=xx(3) c return c end c c----------------------------------------------------------------------- subroutine utroa1(phi,a1,a2,a3,x1,x2,x3) c----------------------------------------------------------------------- c rotates x by angle phi around axis a (argument single precision) c normalization of a is irrelevant. c----------------------------------------------------------------------- double precision aa(3),xx(3),dphi dphi=phi xx(1)=x1 xx(2)=x2 xx(3)=x3 aa(1)=a1 aa(2)=a2 aa(3)=a3 call utroa2(dphi,aa(1),aa(2),aa(3),xx(1),xx(2),xx(3)) c back to single precision x1=sngl(xx(1)) x2=sngl(xx(2)) x3=sngl(xx(3)) return end c----------------------------------------------------------------------- subroutine utroa2(phi,a1,a2,a3,x1,x2,x3) c----------------------------------------------------------------------- c rotates x by angle phi around axis a. c normalization of a is irrelevant. c double precision phi,a1,a2,a3,x1,x2,x3 c----------------------------------------------------------------------- double precision phi,a1,a2,a3,x1,x2,x3 double precision aaa,aa(3),xxx,xx(3),e1(3),e2(3),e3(3),xp,xt,dphi dphi=phi xx(1)=x1 xx(2)=x2 xx(3)=x3 aa(1)=a1 aa(2)=a2 aa(3)=a3 aaa=0d0 xxx=0d0 do i=1,3 aaa=aaa+aa(i)**2 xxx=xxx+xx(i)**2 enddo if(xxx.eq.0d0)return if(aaa.eq.0d0)call utstop('utroa1: zero rotation axis&') aaa=1.0/dsqrt(aaa) c e3 = a / !a! do i=1,3 e3(i)=aa(i)*aaa enddo c x_parallel xp=0 do i=1,3 xp=xp+xx(i)*e3(i) enddo c x_transverse if(xxx-xp**2.le.0.)return xt=dsqrt(xxx-xp**2) c e1 = vector x_transverse / absolute value x_transverse do i=1,3 e1(i)=(xx(i)-e3(i)*xp)/xt enddo c e2 orthogonal e3,e1 call utvec2(e3,e1,e2) c rotate x do i=1,3 xx(i)=xp*e3(i)+xt*cos(dphi)*e1(i)+xt*sin(dphi)*e2(i) enddo xxx=0d0 do i=1,3 xxx=xxx+xx(i)**2 enddo c back to single precision x1=xx(1) x2=xx(2) x3=xx(3) return end cc-------------------------------------------------------------------- c function utroot(funcd,x1,x2,xacc) cc-------------------------------------------------------------------- cc combination of newton-raphson and bisection method for root finding cc input: cc funcd: subr returning fctn value and first derivative cc x1,x2: x-interval cc xacc: accuracy cc output: cc utroot: root cc-------------------------------------------------------------------- c include 'epos.inc' c INTEGER MAXIT c REAL utroot,x1,x2,xacc c EXTERNAL funcd c PARAMETER (MAXIT=100) c INTEGER j c REAL df,dx,dxold,f,fh,fl,temp,xh,xl c call funcd(x1,fl,df) c call funcd(x2,fh,df) c if((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.)) c *call utstop('utroot: root must be bracketed&') c if(fl.eq.0.)then c utroot=x1 c return c else if(fh.eq.0.)then c utroot=x2 c return c else if(fl.lt.0.)then c xl=x1 c xh=x2 c else c xh=x1 c xl=x2 c endif c utroot=.5*(x1+x2) c dxold=abs(x2-x1) c dx=dxold c call funcd(utroot,f,df) c do 11 j=1,MAXIT c if(((utroot-xh)*df-f)*((utroot-xl)*df-f).ge.0..or. abs(2.* c *f).gt.abs(dxold*df) ) then c dxold=dx c dx=0.5*(xh-xl) c utroot=xl+dx c if(xl.eq.utroot)return c else c dxold=dx c dx=f/df c temp=utroot c utroot=utroot-dx c if(temp.eq.utroot)return c endif c if(abs(dx).lt.xacc) return c call funcd(utroot,f,df) c if(f.lt.0.) then c xl=utroot c else c xh=utroot c endif c11 continue c call utmsg('utroot') c write(ifch,*)'***** exceeding maximum iterations' c write(ifch,*)'dx:',dx c call utmsgf c return c END c c----------------------------------------------------------------------- subroutine utrot2(isig,ax,ay,az,x,y,z) c----------------------------------------------------------------------- c performs a rotation, double prec. c----------------------------------------------------------------------- include 'epos.inc' double precision ax,ay,az,x,y,z,rx,ry,rz *,alp,bet,cosa,sina,cosb,sinb,xs,ys,zs if(ax**2.eq.0.and.ay**2.eq.0.and.az**2.eq.0.)then write(ifch,*)'ax**2,ay**2,az**2:',ax**2,ay**2,az**2 write(ifch,*)'ax,ay,az:',ax,ay,az call utstop('utrot2: zero vector.&') endif if(az.ge.0.)then rx=ax ry=ay rz=az else rx=-ax ry=-ay rz=-az endif if(rz**2+ry**2.ne.0.)then alp=dabs(dacos(rz/dsqrt(rz**2+ry**2)))*sign(1.,sngl(ry)) bet= *dabs(dacos(dsqrt(rz**2+ry**2)/dsqrt(rz**2+ry**2+rx**2)))* *sign(1.,sngl(rx)) else alp=3.1415927d0/2d0 bet=3.1415927d0/2d0 endif cosa=dcos(alp) sina=dsin(alp) cosb=dcos(bet) sinb=dsin(bet) if(isig.ge.0)then xs=x*cosb-y*sina*sinb-z*cosa*sinb ys= y*cosa -z*sina zs=x*sinb+y*sina*cosb+z*cosa*cosb else !if(isig.lt.0)then xs= x*cosb +z*sinb ys=-x*sinb*sina+y*cosa+z*cosb*sina zs=-x*sinb*cosa-y*sina+z*cosb*cosa endif x=xs y=ys z=zs return end c----------------------------------------------------------------------- subroutine utrot4(isig,ax,ay,az,x,y,z) c----------------------------------------------------------------------- c performs a rotation, double prec. c arguments partly single c----------------------------------------------------------------------- double precision ax,ay,az,xx,yy,zz xx=dble(x) yy=dble(y) zz=dble(z) call utrot2(isig,ax,ay,az,xx,yy,zz) x=sngl(xx) y=sngl(yy) z=sngl(zz) return end c----------------------------------------------------------------------- subroutine utrota(isig,ax,ay,az,x,y,z) c----------------------------------------------------------------------- c performs a rotation c----------------------------------------------------------------------- if(az.ge.0.)then rx=ax ry=ay rz=az else rx=-ax ry=-ay rz=-az endif if(rz.eq.0..and.ry.eq.0.)then alp=0. stop else alp=abs(utacos(rz/sqrt(rz**2+ry**2)))*sign(1.,ry) endif bet= *abs(utacos(sqrt(rz**2+ry**2)/sqrt(rz**2+ry**2+rx**2)))*sign(1.,rx) cosa=cos(alp) sina=sin(alp) cosb=cos(bet) sinb=sin(bet) if(isig.ge.0)then xs=x*cosb-y*sina*sinb-z*cosa*sinb ys= y*cosa -z*sina zs=x*sinb+y*sina*cosb+z*cosa*cosb else !if(isig.lt.0)then xs= x*cosb +z*sinb ys=-x*sinb*sina+y*cosa+z*cosb*sina zs=-x*sinb*cosa-y*sina+z*cosb*cosa endif x=xs y=ys z=zs return end c----------------------------------------------------------------------- subroutine utstop(text) c----------------------------------------------------------------------- c returns error message and stops execution. c text is an optonal text to appear in the error message. c text is a character string of length 40; c for shorter text, it has to be terminated by &; c example: call utstop('error in subr xyz&') c----------------------------------------------------------------------- include 'epos.inc' c parameter(itext=40) character text*(*) ,txt*6 imax=index(text,'&') do 1 j=1,2 if(j.eq.1)then ifi=ifch else !if(j.eq.2) ifi=ifmt endif if(imax.gt.1)then write(ifi,101)('*',k=1,72),text(1:imax-1) *,nrevt+1,nint(seedj),seedc,('*',k=1,72) else write(ifi,101)('*',k=1,72),' ' *,nrevt+1,nint(seedj),seedc,('*',k=1,72) endif 101 format( *1x,72a1 */1x,'***** stop in ',a */1x,'***** current event number: ',i12 */1x,'***** initial seed for current run:',i10 */1x,'***** initial seed for current event:',d25.15 */1x,72a1) 1 continue c c=0. c b=a/c stop entry utmsg(txt) imsg=imsg+1 write(ifch,'(1x,74a1)')('*',j=1,72) write(ifch,100)txt,nrevt+1,nint(seedj),seedc 100 format(1x,'***** msg from ',a6,'. es:',i7,2x,i9,2x,d23.17) return entry utmsgf if(ish.eq.1)return write(ifch,'(1x,74a1)')('*',j=1,72) end c----------------------------------------------------------------- subroutine uttrap(func,a,b,s) c----------------------------------------------------------------- c trapezoidal method for integration. c input: fctn func and limits a,b c output: value s of the integral c----------------------------------------------------------------- include 'epos.inc' INTEGER JMAX REAL a,b,func,s EXTERNAL func PARAMETER (JMAX=10) CU USES uttras INTEGER j REAL olds olds=-1.e30 do 11 j=1,JMAX if(ish.ge.9)write(ifch,*)'sr uttrap: j:',j call uttras(func,a,b,s,j) ds=abs(s-olds) if (ds.lt.epsr*abs(olds)) return olds=s 11 continue c-c nepsr=nepsr+1 if(ish.ge.9)then call utmsg('uttrap') write(ifch,*) *'***** requested accuracy could not be achieved' write(ifch,*)'achieved accuracy: ',ds/abs(olds) write(ifch,*)'requested accuracy:',epsr call utmsgf endif END c----------------------------------------------------------------- subroutine uttraq(func,a,b,s) c----------------------------------------------------------------- c trapezoidal method for integration. c input: function func and limits a,b c output: value s of the integral c----------------------------------------------------------------- REAL a,b,func,s EXTERNAL func PARAMETER (eps=1.e-6) CU USES uttras INTEGER j REAL olds olds=-1.e30 j=1 10 call uttras(func,a,b,s,j) ds=abs(s-olds) if (ds.le.eps*abs(olds)) return olds=s if(j.ge.15)return j=j+1 goto10 END c----------------------------------------------------------------- subroutine uttras(func,a,b,s,n) c----------------------------------------------------------------- c performs one iteration of the trapezoidal method for integration c----------------------------------------------------------------- INTEGER n REAL a,b,s,func EXTERNAL func INTEGER it,j REAL del,sum,tnm,x if (n.eq.1) then s=0.5*(b-a)*(func(a)+func(b)) else it=2**(n-2) tnm=it del=(b-a)/tnm x=a+0.5*del sum=0. do 11 j=1,it sum=sum+func(x) x=x+del 11 continue s=0.5*(s+(b-a)*sum/tnm) endif return END cc----------------------------------------------------------------------- c subroutine utvec1(a,b,c) cc----------------------------------------------------------------------- cc returns vector product c = a x b . cc----------------------------------------------------------------------- c real a(3),b(3),c(3) c c(1)=a(2)*b(3)-a(3)*b(2) c c(2)=a(3)*b(1)-a(1)*b(3) c c(3)=a(1)*b(2)-a(2)*b(1) c return c end c c----------------------------------------------------------------------- subroutine utvec2(a,b,c) c----------------------------------------------------------------------- c returns vector product c = a x b . c a,b,c double precision. c----------------------------------------------------------------------- double precision a(3),b(3),c(3) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end c------------------------------------------------------------------- subroutine utword(line,i,j,iqu) c------------------------------------------------------------------- c finds the first word of the character string line(j+1:1000). c the word is line(i:j) (with new i and j). c if j<0 or if no word found --> new line read. c a text between quotes "..." is considered one word; c stronger: a text between double quotes ""..."" is consid one word c stronger: a text between "{ and }" is considered one word c------------------------------------------------------------------- c input: c line: character string (*1000) c i: integer between 1 and 1000 c iqu: for iqu=1 a "?" is written to output before reading a line, c otherwise (iqu/=1) nothing is typed c output: c i,j: left and right end of word (word=line(i:j)) c------------------------------------------------------------------- include 'epos.inc' parameter(mempty=2) character*1 empty(mempty),mk character line*1000 character*2 mrk data empty/' ',','/ parameter(mxdefine=40) character w1define*100,w2define*100 common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine) & ,w1define(mxdefine),w2define(mxdefine) j0=0 if(j.ge.0)then i=j goto 1 endif 5 continue if(iqu.eq.1.and.iprmpt.gt.0)write(ifmt,'(a)')'?' if(nopen.eq.0)then ifopx=ifop elseif(nopen.gt.0)then ifopx=20+nopen else !if(nopen.lt.0) ifopx=ifcp endif read(ifopx,'(a1000)',end=9999)line if(iecho.eq.1.or.(nopen.ge.0.and.kcpopen.eq.1))then kmax=2 do k=3,1000 if(line(k:k).ne.' ')kmax=k enddo else kmax=0 endif if(nopen.ge.0.and.kcpopen.eq.1) & write(ifcp,'(a)')line(1:kmax) if(iecho.eq.1) & write(ifmt,'(a)')line(1:kmax) i=0 1 i=i+1 if(i.gt.1000)goto 5 if(line(i:i).eq.'!')goto 5 do ne=1,mempty if(line(i:i).eq.empty(ne))goto 1 enddo nbla=1 mrk=' ' mk=' ' if(line(i:i).eq.'~')mk='~' if(line(i:i+1).eq.'"{')mrk='}"' if(line(i:i+1).eq.'""')mrk='""' if(mrk.ne.' ')goto 10 if(line(i:i).eq.'"')mk='"' if(mk.ne.' ')goto 8 j=i-1 6 j=j+1 if(j.gt.1000)goto 7 if(line(j:j).eq.'!')goto 7 do ne=1,mempty if(line(j:j).eq.empty(ne))goto 7 enddo goto 6 8 continue if(i.ge.1000-1)stop'utword: make line shorter!!! ' i=i+1 j=i if(line(j:j).eq.mk)stop'utword: empty string!!! ' 9 j=j+1 if(j.gt.1000)then !reach the end of the line j=j-nbla+2 goto 7 endif if(line(j:j).eq.' ')then nbla=nbla+1 else nbla=2 endif if(line(j:j).eq.mk)then line(i-1:i-1)=' ' line(j:j)=' ' goto 7 endif goto 9 10 continue if(i.ge.1000-3)stop'utword: make line shorter!!!! ' i=i+2 j=i if(line(j:j+1).eq.mrk)stop'utword: empty string!!!! ' 11 j=j+1 if(j.gt.1000-1)then !reach the end of the line j=j-nbla+2 goto 7 endif if(line(j:j+1).eq.mrk)then line(i-2:i-1)=' ' line(j:j+1)=' ' goto 7 endif if(line(j:j).eq.' ')then nbla=nbla+1 else nbla=2 endif goto 11 7 j=j-1 !--------#define--------------- if(ndefine.gt.0)then do ndf=1,ndefine l1=l1define(ndf) l2=l2define(ndf) do i0=i,j+1-l1 if(line(i0:i0-1+l1).eq.w1define(ndf)(1:l1))then if(l2.eq.l1)then line(i0:i0-1+l1)=w2define(ndf)(1:l2) elseif(l2.lt.l1)then line(i0:i0+l2-1)=w2define(ndf)(1:l2) do k=i0+l2,i0-1+l1 line(k:k)=' ' enddo elseif(l2.gt.l1)then do k=i0+l1,i0+l2-1 if(line(k:k).ne.' ') & stop'utword: no space for `define` replacement. ' enddo line(i0:i0+l2-1)=w2define(ndf)(1:l2) j=i0+l2-1 endif endif enddo enddo do k=i,j if(line(k:k).ne.' ')j0=j enddo j=j0 endif !-------- return 9999 close(ifopx) nopen=nopen-1 if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1 goto 5 end c-------------------------------------------------------------------- subroutine utworn(line,j,ne) c-------------------------------------------------------------------- c returns number ne of nonempty characters of line(j+1:1000) c-------------------------------------------------------------------- character line*1000 ne=0 do l=j+1,1000 if(line(l:l).ne.' ')ne=ne+1 enddo return end c----------------------------------------------------------------------- subroutine getairmol(iz,ia) c----------------------------------------------------------------------- include 'epos.inc' i=0 r=rangen() do while(r.gt.0.) ! choose air-molecule i=i+1 r=r-airwnxs(i) enddo iz = nint(airznxs(i)) ia = nint(airanxs(i)) end c---------------------------------------------------------------------- subroutine factoriel c---------------------------------------------------------------------- c tabulation of fctrl(n)=n!, facto(n)=1/n! and utgamtab(x) for x=0 to 50 c---------------------------------------------------------------------- include 'epos.incems' double precision utgamtab,utgam,x common/gamtab/utgamtab(10000) nfctrl=100 fctrl(0)=1.D0 facto(0)=1.D0 do i=1,min(npommx,nfctrl) fctrl(i)=fctrl(i-1)*dble(i) facto(i)=1.d0/fctrl(i) enddo do k=1,10000 x=dble(k)/100.d0 utgamtab(k)=utgam(x) enddo return end c----------------------------------------------------------------------- subroutine fremnu(amin,ca,cb,ca0,cb0,ic1,ic2,ic3,ic4) c----------------------------------------------------------------------- common/hadr2/iomodl,idproj,idtarg,wexcit real pnll,ptq common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg parameter (nflav=6) integer ic(2),jc(nflav,2) ic(1)=ca ic(2)=cb call iddeco(ic,jc) keu=jc(1,1)-jc(1,2) ked=jc(2,1)-jc(2,2) kes=jc(3,1)-jc(3,2) kec=jc(4,1)-jc(4,2) keb=jc(5,1)-jc(5,2) ket=jc(6,1)-jc(6,2) amin=utamnu(keu,ked,kes,kec,keb,ket,5) !???4=2mults, 5=1mult if(ca-ca0.eq.0..and.cb-cb0.eq.0..and.rangen().gt.wexcit)then ic3=0 ic4=0 ic1=ca ic2=cb else amin=amin+exmass n=0 do i=1,4 do j=1,2 n=n+jc(i,j) enddo enddo k=1+rangen()*n do i=1,4 do j=1,2 k=k-jc(i,j) if(k.le.0)goto 1 enddo enddo 1 if(j.eq.1)then ic3=10**(6-i) ic4=0 else ic3=0 ic4=10**(6-i) endif ic1=int(ca)-ic3 ic2=int(cb)-ic4 endif return end c----------------------------------------------------------------------- function fremnux(jc) c----------------------------------------------------------------------- real pnll,ptq common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg parameter (nflav=6) integer jc(nflav,2)!,ic(2) c ic(1)=210000 c ic(2)=0 c call iddeco(ic,jc) keu=jc(1,1)-jc(1,2) ked=jc(2,1)-jc(2,2) kes=jc(3,1)-jc(3,2) kec=jc(4,1)-jc(4,2) keb=jc(5,1)-jc(5,2) ket=jc(6,1)-jc(6,2) fremnux=utamnu(keu,ked,kes,kec,keb,ket,4) !+exmass !???4=2mults, 5=1mult return end c----------------------------------------------------------------------- function fremnux2(jc) c----------------------------------------------------------------------- real pnll,ptq common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg parameter (nflav=6) integer jc(nflav,2)!,ic(2) c ic(1)=210000 c ic(2)=0 c call iddeco(ic,jc) keu=jc(1,1)-jc(1,2) ked=jc(2,1)-jc(2,2) kes=jc(3,1)-jc(3,2) kec=jc(4,1)-jc(4,2) keb=jc(5,1)-jc(5,2) ket=jc(6,1)-jc(6,2) fremnux2=utamnu(keu,ked,kes,kec,keb,ket,5) !+exmass !???4=2mults, 5=1mult return end c----------------------------------------------------------------------- function fremnux3(jci) c----------------------------------------------------------------------- c minimum mass from ic counting all quarks c----------------------------------------------------------------------- include 'epos.inc' integer jc(nflav,2),jci(nflav,2)!,ic(2) c ic(1)=210000 c ic(2)=0 c print *,'start',ic fremnux3=0. do j=1,2 do i=1,nflav jc(i,j)=jci(i,j) enddo enddo c call iddeco(ic,jc) call idquacjc(jc,nqua,naqu) do ii=1,2 if(ii.eq.1)then nqu=nqua else nqu=naqu endif if(nqu.ge.3)then do while(jc(3,ii).ne.0.and.nqu.ge.3) !count baryons with s quark jc(3,ii)=jc(3,ii)-1 if(jc(3,ii).gt.0)then jc(3,ii)=jc(3,ii)-1 if(jc(3,ii).gt.0)then jc(3,ii)=jc(3,ii)-1 fremnux3=fremnux3+asuhax(4) elseif(jc(2,ii).gt.0)then jc(2,ii)=jc(2,ii)-1 fremnux3=fremnux3+asuhax(4) elseif(jc(1,ii).gt.0)then jc(1,ii)=jc(1,ii)-1 fremnux3=fremnux3+asuhax(4) endif elseif(jc(2,ii).gt.0)then jc(2,ii)=jc(2,ii)-1 if(jc(1,ii).gt.0)then jc(1,ii)=jc(1,ii)-1 fremnux3=fremnux3+asuhax(3) elseif(jc(2,ii).gt.0)then jc(2,ii)=jc(2,ii)-1 fremnux3=fremnux3+asuhax(3) endif elseif(jc(1,ii).gt.0)then jc(1,ii)=jc(1,ii)-2 fremnux3=fremnux3+asuhay(3) endif nqu=nqu-3 enddo do while(jc(2,ii).ne.0.and.nqu.ge.3) !count baryons with d quark jc(2,ii)=jc(2,ii)-1 if(jc(1,ii).gt.0)then jc(1,ii)=jc(1,ii)-1 if(jc(2,ii).gt.0)then jc(2,ii)=jc(2,ii)-1 fremnux3=fremnux3+asuhay(2) elseif(jc(1,ii).gt.0)then jc(1,ii)=jc(1,ii)-1 fremnux3=fremnux3+asuhay(2) endif elseif(jc(2,ii).gt.0)then jc(2,ii)=jc(2,ii)-2 fremnux3=fremnux3+asuhay(3) endif nqu=nqu-3 enddo do while(jc(1,ii).ne.0.and.nqu.ge.3) !count baryons with s quark jc(1,ii)=jc(1,ii)-3 fremnux3=fremnux3+asuhay(3) nqu=nqu-3 enddo if(ii.eq.1)then nqua=nqu else naqu=nqu endif endif c print *,ii,nqua,naqu,jc,fremnux3 enddo if(nqua+naqu.ne.0)then do while(jc(3,1).ne.0) !count mesons with s quark jc(3,1)=jc(3,1)-1 if(jc(3,2).gt.0)then jc(3,2)=jc(3,2)-1 fremnux3=fremnux3+asuhax(6) elseif(jc(2,2).gt.0)then jc(2,2)=jc(2,2)-1 fremnux3=fremnux3+asuhay(6) elseif(jc(1,2).gt.0)then jc(1,2)=jc(1,2)-1 fremnux3=fremnux3+asuhay(6) endif enddo do while(jc(2,1).ne.0) !count mesons with d quark jc(2,1)=jc(2,1)-1 if(jc(2,2).gt.0)then jc(2,2)=jc(2,2)-1 fremnux3=fremnux3+asuhay(5) elseif(jc(1,2).gt.0)then jc(1,2)=jc(1,2)-1 fremnux3=fremnux3+asuhay(5) endif enddo do while(jc(1,1).ne.0) !count mesons with s quark jc(1,1)=jc(1,1)-1 if(jc(1,2).gt.0)then jc(1,2)=jc(1,2)-1 fremnux3=fremnux3+asuhay(5) endif enddo endif c fremnux3=fremnux3+0.5 c print *,'stop',nqua,naqu,fremnux3 return end c----------------------------------------------------------------------- subroutine fremnx(ammax,amin,sm,ic3,ic4,iret) c----------------------------------------------------------------------- common/psar9/ alpr include 'epos.inc' iret=0 if(ic3.eq.0.and.ic4.eq.0)then if(ammax.lt.amin**2)then iret=1 return endif sm=amin**2 else c ammax1=min(ammax,(engy/4.)**2) ammax1=ammax if(ammax1.lt.amin**2)then iret=1 return endif if(alpr.eq.-1.)then sm=amin**2*(ammax1/amin**2)**rangen() else sm=amin**2*(1.+((ammax1/amin**2)**(1.+alpr)-1.) * *rangen())**(1./(1.+alpr)) endif endif return end SUBROUTINE gaulag(x,w,n,alf) INTEGER n,MAXIT REAL alf,w(n),x(n) DOUBLE PRECISION EPS PARAMETER (EPS=3.D-14,MAXIT=10) CU USES gammln INTEGER i,its,j REAL ai,gammln DOUBLE PRECISION p1,p2,p3,pp,z,z1 z=0. do 13 i=1,n if(i.eq.1)then z=(1.+alf)*(3.+.92*alf)/(1.+2.4*n+1.8*alf) else if(i.eq.2)then z=z+(15.+6.25*alf)/(1.+.9*alf+2.5*n) else ai=i-2 z=z+((1.+2.55*ai)/(1.9*ai)+1.26*ai*alf/(1.+3.5*ai))* *(z-x(i-2))/(1.+.3*alf) endif do 12 its=1,MAXIT p1=1.d0 p2=0.d0 do 11 j=1,n p3=p2 p2=p1 p1=((2*j-1+alf-z)*p2-(j-1+alf)*p3)/j 11 continue pp=(n*p1-(n+alf)*p2)/z z1=z z=z1-p1/pp if(abs(z-z1).le.EPS)goto 1 12 continue call utstop("too many iterations in gaulag") 1 x(i)=z w(i)=-exp(gammln(alf+n)-gammln(float(n)))/(pp*n*p2) 13 continue return END C (C) Copr. 1986-92 Numerical Recipes Software 4+1$!]. FUNCTION gammln(xx) REAL gammln,xx INTEGER j DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) SAVE cof,stp DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, *24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, *-.5395239384953d-5,2.5066282746310005d0/ x=xx y=x tmp=x+5.5d0 tmp=(x+0.5d0)*log(tmp)-tmp ser=1.000000000190015d0 do 11 j=1,6 y=y+1.d0 ser=ser+cof(j)/y 11 continue gammln=tmp+log(stp*ser/x) return END C (C) Copr. 1986-92 Numerical Recipes Software 4+1$!]. function polar(x,y) pi=3.1415927 if(abs(x).gt.1.e-6)then phi=atan(y/x) if(x.lt.0.)phi=pi+phi if(phi.lt.0)phi=2*pi+phi else phi=0.5*pi if(y.lt.0)phi=phi+pi endif polar=phi end subroutine getJKNcentr end