common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
*,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
- common /cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
+ common /cttaus/ tpro,zpro,ttar,ztar,ttaus,detap,detat
integer jc(nflav,2),ic(2)
call utpri('amicro',ish,ishini,4)
$ .lt. rangen() ) goto 15
kes=nint(x)
endif
- if(real((keu+ked+kes)/3).ne.real(keu+kes+ked)/3.) goto 10
+ if(real((keu+ked+kes)/3).ne.real(keu+kes+ked)/3.) goto 10
if(keu.ge.0)jc(1,1)=keu
if(ked.ge.0)jc(2,1)=ked
if(kes.ge.0)jc(3,1)=kes
idr=8*10**8+ic(1)*100+ic(2)/100
if(ish.ge.5)write(ifch,'(a,i9)')' id:',idr
dez=0.5e-4
- else
+ else
idr=idr
* +mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
* +mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
radptl(nptl)=(3*volu/pi/4.)**0.33333
istptl(nptl)=10
tivptl(2,1)=1.
-
+
nptlb=nptl
call hnbaaa(nptl,iret)
if(iret.eq.1)stop'STOP: sr amicro: hnbaaa return code = 1. '
istptl(n)=0
ifrptl(1,n)=0
ifrptl(2,n)=0
- xorptl(1,n)=x
- xorptl(2,n)=y
+ xorptl(1,n)=x
+ xorptl(2,n)=y
xorptl(3,n)=z
xorptl(4,n)=t
tivptl(1,n)=t
subroutine hgcaaa
c-----------------------------------------------------------------------
c hadronic resonance gas in grand canonical treatment
-c returns T, chemical potentials and hadronic yield
-c (hadron chemical potentials as combinations of quark chemical potentials)
+c returns T, chemical potentials and hadronic yield
+c (hadron chemical potentials as combinations of quark chemical potentials)
c
c input:
c iostat: 1: Boltzmann approximation, 0: quantum statistics /metr3/
c and employing the same algorithm as for massive
c-----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cbol/rmsbol(mspecs),ptlbol(mspecs),chebol(mspecs),tembol
gfac=gfac+gspecs(i)
endif
enddo
- if(iabs(ispecs(nspecs)).lt.10)gfac=gfac+16.
+ if(iabs(ispecs(nspecs)).lt.10)gfac=gfac+16.
tem=(tecm/volu*hquer**3*30./pi**2/gfac)**.25
else
do i=1,nspecs
if(ish.ge.5.and.iospec.ne.iug)then
write(ifch,*)'inversion in Boltzmann approx. :'
- elseif(ish.ge.5.and.iospec.eq.iug)then
+ elseif(ish.ge.5.and.iospec.eq.iug)then
write(ifch,*)'inversion for massless hadrons :'
endif
- if(ish.ge.5)then
+ if(ish.ge.5)then
if(nflavs.eq.1)write(ifch,'(3x,a,8x,a)')
*'T:','chemu:'
if(nflavs.eq.2)write(ifch,'(3x,a,8x,a,5x,a)')
if(i.eq.3.and.ish.ge.5)
*write(ifch,'(1x,a,1x,f9.6)')'chems:',chem(3)
enddo
-
+
c checking results
c ----------------
ish=isho
return
endif
-
+
c quantum statistics
c ------------------
if(ish.ge.5)write(ifch,*)'quantum statistics:'
- if(ish.ge.5.and.nflavs.eq.1)write(ifch,'(3x,a,8x,a)')
+ if(ish.ge.5.and.nflavs.eq.1)write(ifch,'(3x,a,8x,a)')
*'T:','chemu:'
- if(ish.ge.5.and.nflavs.eq.2)write(ifch,'(3x,a,8x,a,6x,a)')
+ if(ish.ge.5.and.nflavs.eq.2)write(ifch,'(3x,a,8x,a,6x,a)')
*'T:','chemu:','chemd:'
- if(ish.ge.5.and.nflavs.eq.3)write(ifch,'(3x,a,8x,a,6x,a,6x,a)')
+ if(ish.ge.5.and.nflavs.eq.3)write(ifch,'(3x,a,8x,a,6x,a,6x,a)')
*'T:','chemu:','chemd:','chems:'
k=1
c checking results
c ----------------
- if(ish.ge.5)call hgcchh(i)
+ if(ish.ge.5)call hgcchh(i)
c particle yield
c --------------
if(ish.ge.5)write(ifch,*)('-',i=1,30)
*,' exit sr hgcaaa ',('-',i=1,10)
ish=isho
- return
+ return
endif
if(k.gt.300)then
c-------------------------------------------------------------------
function hgcbk(n,x)
-c------------------------------------------------------------------
+c------------------------------------------------------------------
tox=2.0/x
bkm=hgcbk0(x)
bk=hgcbk1(x)
END
-c----------------------------------------------------------------
+c----------------------------------------------------------------
subroutine hgccbo(iba)
c----------------------------------------------------------------
c returns new chem(iafs) for boltzmann statistics
c-----------------------------------------------------------------------
common/cnsta/pi,pii,hquer,prom,piom,ainfin
common/drop6/tecm,volu
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
parameter(nflav=6)
11 continue
fd=0.0
call hgchac(0)
-
+
do i=1,nspecs
if(ifok(iafs,i).ne.0)then
if(k.gt.300)return
goto10
-
+
end
c plots iteration values for T and mu_i
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
parameter (nbin=500)
write(cu,'(i3)')keu
write(cd,'(i3)')ked
write(cs,'(i3)')kes
-
+
write(ifhi,'(a)') 'newpage zone 1 4 1 openhisto'
write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
enddo
write(ifhi,'(a)') ' endarray'
write(ifhi,'(a)') 'closehisto plot 0'
-
+
write(ifhi,'(a)') 'openhisto'
write(ifhi,'(a)') 'htyp lin xmod lin ymod lin'
write(ifhi,'(a)') 'text 0 0 "xaxis Iteration"'
write(ifhi,'(a)') 'closehisto plot 0'
endif
-
+
return
-
+
end
c-----------------------------------------------------------------------
c chem(iafs)
c-----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cflavs/nflavs,kef(nflav),chem(nflav)
common/cflac/ifok(nflav,mspecs),ifoa(nflav)
if(ifok(iafs,ians).ne.0)then
call hgchac(0)
- call hgclim(a,b)
+ call hgclim(a,b)
if(b.eq.0.0)then
hpd=0.0
else
c checks flavor conservation in particle yield
c------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
if(dkef.le.1.e-2)then
if(i.eq.1.and.ish.ge.5)write(ifch,*)'u conserved'
if(i.eq.2.and.ish.ge.5)write(ifch,*)'d conserved'
- if(i.eq.3.and.ish.ge.5)write(ifch,*)'s conserved'
+ if(i.eq.3.and.ish.ge.5)write(ifch,*)'s conserved'
else
if(i.eq.1.and.ish.ge.5)write(ifch,*)'u not conserved'
if(i.eq.2.and.ish.ge.5)write(ifch,*)'d not conserved'
endif
enddo
- return
+ return
end
c----------------------------------------------------------------
c checks results by numerical integration
c----------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
c checks results by numerical integration
c----------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
endif
hfd=ifok(i,ians)*hpd*gspecs(ians)/2./pi**2/hquer**3
if(ish.ge.9)write(ifch,*)'hfd:',hfd
- cfd=cfd+hfd
+ cfd=cfd+hfd
enddo
if(i.eq.1.and.ish.ge.5)write(ifch,5)'flavor density u :',cfd
if(i.eq.2.and.ish.ge.5)write(ifch,5)'flavor density d :',cfd
if(ish.ge.5)write(ifch,*)'results disagree'
endif
- return
+ return
end
c chem
c---------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
c endif
c hpd=hpd*gspecs(i)*tem**3/pi**2/hquer**3
c if(chemgc(i).eq.abs(chemgc(i)))then
-c hpd=gspecs(i)*(chemgc(i)*tem**2+chemgc(i)**3/pi**2)/6./hquer**3
+c hpd=gspecs(i)*(chemgc(i)*tem**2+chemgc(i)**3/pi**2)/6./hquer**3
c *-hpd
c endif
c else
c hpd=3.*gspecs(i)*tem**3*z3/4./pi**2/hquer**3
c endif
-
+
else
hpd=gspecs(i)*tem**3*z3/pi**2/hquer**3
c-----------------------------------------------------------------------
c integrand of energy density
c------------------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
c-----------------------------------------------------------------
c integrand of mean square variance of energy
c----------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
c-----------------------------------------------------------------
c integrand of hadron density
c-----------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
c-----------------------------------------------------------------------
c integrand of energy density
c------------------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
if(mod(igsp,2).ne.0)then
d=-1.0
- if(eex.lt.1.e-10)return
+ if(eex.lt.1.e-10)return
else
d=1.0
- endif
+ endif
hgcfhe=sq*x**2/(exp(eex)+d)
- return
+ return
end
c-----------------------------------------------------------------
c-----------------------------------------------------------------
c integrand of mean square variance of energy
c----------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
c-----------------------------------------------------------------
c integrand of hadron density
c-----------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
if(mod(igsp,2).ne.0)then
d=-1.0
- if(eex.lt.1.e-10)return
+ if(eex.lt.1.e-10)return
else
d=1.0
- endif
+ endif
hgcfhn=x**2/(exp(eex)+d)
c-----------------------------------------------------------------
c integrand of mean square variance of hadron yield
c----------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
sq=sqrt(x**2+aspecs(ians)**2)
if(tem.ne.0.0)eex=(sq-chemgc(ians))/tem
- if(eex.gt.60.)return
- if(eex.lt.(-60.))return
+ if(eex.gt.60.)return
+ if(eex.lt.(-60.))return
if(mod(igsp,2).ne.0)then
d=-1.0
- if(eex.lt.1.0e-10.and.eex.gt.(-1.0e-10))return
+ if(eex.lt.1.0e-10.and.eex.gt.(-1.0e-10))return
else
d=1.0
- endif
+ endif
hgcfhw=x**2/(exp(eex)+2.0*d+exp(-eex))
c chemical potentials
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
c and energy densities using quantum statistics
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
include 'epos.inc'
parameter(maxp=500)
common/chnbin/nump,ihadro(maxp)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cgctot/rmstot,ptltot
write(ifch,*)'impossible to generate hadron configuration'
call utstop('hgcnbi: tecm less than two pi0 masses&')
endif
-
+
kk=1
100 continue
enddo
if(ish.ge.7)write(ifch,*)
- *'sample hadron multiplicities and total mass:'
+ *'sample hadron multiplicities and total mass:'
kbar=keu+ked+kes
kpar=iabs(keu)+iabs(ked)+iabs(kes)
if(ish.ge.9)write(ifch,*)'pnlog:',pnlog
if(pnlog.lt.60)then
pn=exp(pnlog)
- else
+ else
pn=1.1
endif
102 continue
-
+
ndd=0
c if(nbb.lt.nb)then
c nba=nb-nbb
c if(nbar.gt.0)then
-c if(ish.ge.7)write(ifch,*)'add protons: nba:',nba
+c if(ish.ge.7)write(ifch,*)'add protons: nba:',nba
c nptlgc(19)=nptlgc(19)+nba
c n=n+nba
c amtot=amtot+aspecs(19)*nba
c elseif(nbar.lt.0)then
-c if(ish.ge.7)write(ifch,*)'add aprotons: nba:',nba
+c if(ish.ge.7)write(ifch,*)'add aprotons: nba:',nba
c nptlgc(20)=nptlgc(20)+nba
c n=n+nba
c amtot=amtot+aspecs(20)*nba
kk=kk+1
goto100
endif
-
+
iii=0
if(ish.ge.7)then
pna=0.0
endif
endif
-
+
if(ptlngc(ib).gt.5.0)then
pnbli=hgcpnl(ib,0)
pnblf=hgcpnl(ib,1)
else
pnb=1.1
endif
- else
+ else
pnb=ptlngc(ib)/(nptlgc(ib)+1)
endif
else
r=rangen()
endif
-
+
c else
c r=1.0
enddo
endif
-
+
if(k(1).ne.0.or.k(2).ne.0.or.k(3).ne.0)then
ll=ll+1
if(ll.le.llmax)then
chitot=chitot+chi**2
enddo
call xhgccc(chitot)
-
+
u=0
d=0
s=0
c returns random multiplicity from gaussian distribution for species i
c---------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cgctot/rmstot,ptltot
common/clatt/nlattc,npmax
endif
else
-
+
2 continue
kk=kk+1
p=0.0
function hgcpml(i1,n1,i2,n2)
c--------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/camgc/amgc,samgc,amtot
common/cgcnb/nptlgc(mspecs)
function hgcpnl(i,n)
c--------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cgcnb/nptlgc(mspecs)
if(ish.ge.9)write(ifch,*)'i:',i,' n:',n
c xpar8 strange chem.pot.
c--------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
call uttraq(hgcfbn,a,b,hden)
endif
hd=hden*gspecs(ians)/2./pi**2/hquer**3
-
+
if(ish.ge.7)write(ifch,*)'i:',ians,' n_u:',ifok(1,ians),' hd:',hd
- qd=qd+ifok(1,ians)*hd+ifok(2,ians)*hd
+ qd=qd+ifok(1,ians)*hd+ifok(2,ians)*hd
if(qd.gt.ymax)qd=ymax
c if(qd.gt.ymax)qd=0.0
if(qd.lt.-ymax)qd=-ymax
c if(qd.lt.-ymax)qd=0.0
-
+
if(b.eq.0.0)then
edi=0.0
if(ish.ge.7)write(ifch,*)'i:',ians,' mu:',chemgc(ians)
* ,' edi:',edi
-
+
ed=ed+edi
if(ed.gt.ymax)ed=ymax
c if(ed.gt.ymax)ed=0.0
if(ish.ge.5)write(ifch,*)' ed:',ed,' qd:',qd
edensi(i,ii)=ed
qdensi(i,ii)=qd
-
+
enddo
enddo
c xpar8 strange chem.pot.
c--------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cflavs/nflavs,kef(nflav),chem(nflav)
hn=hn*volu*gspecs(ians)/2./pi**2/hquer**3
hv=hv*volu*gspecs(ians)/2./pi**2/hquer**3
if(ish.ge.5)write(ifch,*)'hn:',hn,' hv:',hv
-
+
hn=max(hn,1.e-15)
qv=qv+hv
qe=qe+hn
if(qv.gt.ymax)qv=ymax
if(qe.gt.ymax)qe=ymax
-
+
if(b.eq.0.0)then
eei=0.0
evi=evi*volu*gspecs(ians)/2./pi**2/hquer**3
if(ish.ge.5)write(ifch,*)'eei:',eei,' evi:',evi
-
+
eei=max(eei,1.e-15)
ev=ev+evi
ee=ee+eei
wn(i)=qfl(i,ii)
v(i)=volu
endif
-
+
enddo
enddo
write(ifhi,'(2e13.5)')v(j),we(j)
enddo
write(ifhi,'(a)') ' endarray'
- write(ifhi,'(a)') 'closehisto plot 0'
+ write(ifhi,'(a)') 'closehisto plot 0'
write(ifhi,'(a)') 'openhisto'
write(ifhi,'(a,2e11.3)')'xrange',xpar4,xpar5
c------------------------------------------------------------------
subroutine hgcpyi(ist)
c------------------------------------------------------------------
-c returns particle yield
+c returns particle yield
c input:
c tem : temperature
c chemgc: chemical potentials
c ist=0 quantum statistics
c--------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cgctot/rmstot,ptltot
rmsngc(ians)=0.0
if(ist.eq.0)then
-
+
call uttraq(hgcfhw,a,b,var)
var=var*gspecs(ians)*volu/2./pi**2/hquer**3
vartot=vartot+var
if(var.ge.0.0)rmsngc(ians)=sqrt(var)
samgc=samgc+var*aspecs(ians)
-
+
else
if(ptlngc(ians).ge.0.0)rmsngc(ians)=sqrt(ptlngc(ians))
*'<N(',ispecs(ians),')> :',ptlngc(ians),'sigma :',rmsngc(ians)
enddo
-
+
if(vartot.ge.0.0)rmstot=sqrt(vartot)
if(samgc.ge.0.0)samgc=sqrt(samgc)
if(amgc.ge.tecm)samgc=sqrt(amgc)
if(ish.ge.5)write(ifch,'(1x,a,2x,f8.4,2x,a,2x,f10.4)')
*'<M_tot> :',amgc,'sigma :',samgc
- return
+ return
end
c------------------------------------------------------------------------
c tem
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
k=1
t1=0.0
t2=1.0
-
+
goto15
10 tem=t1+.5*(t2-t1)
endif
if(k.gt.300)return
-
+
k=k+1
goto10
end
c tem
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
15 continue
if(tem.le.1.e-6)return
eden=0.0
-
+
do ians=1,nspecs
call hgclim(a,b)
if(b.eq.0.0)then
*write(ifch,*)'failure in tex'
return
endif
-
+
k=k+1
goto10
end
c----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/ciakt/gen,iafs,ians,genm
k=1
-
+
t1=0.0
t2=1.0
10 tem=t1+.5*(t2-t1)
*write(ifch,*)'failure in tm0'
return
endif
-
+
k=k+1
goto10
end
integer ifl(nflav)
double precision p(5),c(5)
real u(3),pin(4),pout(4),poutx(4)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
!-------------------------------------------------------------
! 110, 120, -120, 130, -130, 230, -230, 220, 330
!, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
!, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330
!, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
- !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
+ !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
!-------------------------------------------------------------
integer ianti(mspecs)
data ianti/ 1, 3, 2, 5, 4, 7, 6, 8, 9
kes=jc(3,1)-jc(3,2)
n=ndrop(keu,ked,kes)
if(n.eq.0)stop'hnbxxx: n=0'
-
+
c...fill wspecs
e=pptl(5,ip)
if(n.lt.0)ii=ianti(i)
w1=wwspecs(abs(n),k-1,ii)
w2=wwspecs(abs(n),k,ii)
- wspecs(i)=w1+xi*(w2-w1)
+ wspecs(i)=w1+xi*(w2-w1)
enddo
-
+
if(ish.ge.5)then
write(ifch,*)'keu,ked,kes,n:',keu,ked,kes,n
write(ifch,'(9x, 9f6.3)')(wspecs(i),i=1,9)
w12=w12+wspecs(i)
enddo
w32=0 !35,36,41,42,53,54 excluded
- do i=37,40
+ do i=37,40
w32=w32+wspecs(i)
enddo
- do i=43,52
+ do i=43,52
w32=w32+wspecs(i)
enddo
sum=w12+w32
do nf=1,nflav
ifl(nf)=jc(nf,1)-jc(nf,2)
enddo
-
-c...print
-
+
+c...print
+
if(ish.ge.5)then
write(ifch,*)'ip=',ip,' id=',idptl(ip),' e=',sngl(c(5))
write(ifch,*)'jc=',jc
endif
c...generate number of hadrons
-
+
wfac=1.05 !mean increased by factor wfac
aa=wtot*wfac
if(aa.le.70.)then
nhad=max(2,nhad)
endif
if(ish.ge.5)write(ifch,*)'-----> ',nhad,' hadrons'
-
-c...generate first n-2 hadrons
+
+c...generate first n-2 hadrons
do n=1,nhad-2
1 sum=0
elseif(ispecs(i).lt.-1000)then
nbari=-1
else
- nbari=0
+ nbari=0
endif
if(nbari*nbarini.gt.0
* .and.nbari*(nbar-nbari).lt.0.and.rangen().gt.wbar(-nbari))then
goto1
elseif(nbari*(nbar-nbari).lt.0)then
if(ish.ge.5)write(ifch,*)'+++++',wbar(-nbari)
- endif
- nbar=nbar-nbari
+ endif
+ nbar=nbar-nbari
nptl=nptl+1
id=ispecs(i)
- idptl(nptl)=id
+ idptl(nptl)=id
call idquac(nptl,nq,ns,na,jc1)
do nf=1,nflav
ifl(nf)=ifl(nf)-jc1(nf,1)+jc1(nf,2)
enddo
if(ish.ge.5)
- * write(ifch,*)'nptl=',nptl,' id=',id,' ifl=',ifl
+ * write(ifch,*)'nptl=',nptl,' id=',id,' ifl=',ifl
enddo
do nf=1,nflav
if(ifl(nf).ge.0)then
else
jc(nf,1)=0
jc(nf,2)=-ifl(nf)
- endif
+ endif
enddo
if(ish.ge.5)then
write(ifch,*)'jc=',jc
write(ifch,*)'hadrons:',(idptl(n),n=nptlb+1,nptl)
* ,' --> nbar=',nbar
endif
-
+
c...last two hadrons
-
+
if(nbar.ne.0)then
do n=1,abs(nbar)
- ii=nbar/abs(nbar)
+ ii=nbar/abs(nbar)
i1=idraflz(jc,(3-ii)/2)
i2=idraflz(jc,(3-ii)/2)
i3=idraflz(jc,(3-ii)/2)
if(i1.eq.i2.and.i2.eq.i3)then
id=ii*(i1*1000+i2*100+i3*10+1)
- else
+ else
if(i2.lt.i1)then
- ix=i1
+ ix=i1
i1=i2
i2=ix
- endif
+ endif
if(i3.lt.i2)then
- ix=i2
+ ix=i2
i2=i3
i3=ix
- endif
+ endif
if(i2.lt.i1)then
- ix=i1
+ ix=i1
i1=i2
i2=ix
- endif
+ endif
ispin=0
if(rangen().lt.w32)ispin=1
id=ii*(i1*1000+i2*100+i3*10+ispin)
- endif
+ endif
nptl=nptl+1
- idptl(nptl)=id
+ idptl(nptl)=id
if(ish.ge.5)
- * write(ifch,*)'nptl=',nptl,' baryon=',id,' jc=',jc
+ * write(ifch,*)'nptl=',nptl,' baryon=',id,' jc=',jc
enddo
endif
endif
ii=1
if(i2.lt.i1)then
- ix=i1
+ ix=i1
i1=i2
i2=ix
ii=-1
if(rangen().lt.w1)ispin=1
id=ii*(i1*100+i2*10+ispin)
nptl=nptl+1
- idptl(nptl)=id
+ idptl(nptl)=id
if(ish.ge.5)write(ifch,*)'nptl=',nptl,' nqu=',nqu
- & ,' naq=',naq,' --> meson',id
+ & ,' naq=',naq,' --> meson',id
call idquacjc(jc,nqu,naq)
enddo
write(ifch,*)nmiss,' hadron(s) missing'
write(ifch,*)'hadrons:',(idptl(n),n=nptlb+1,nptl)
endif
-
-c nsechad=lkfoi(1,ifl(1),ifl(2),ifl(3),ifl(4))
-c if(nsechad.gt.0)then
+
+c nsechad=lkfoi(1,ifl(1),ifl(2),ifl(3),ifl(4))
+c if(nsechad.gt.0)then
c i2x=min(nsechad,1+rangen()*nsechad)
c i2=lkfoi(1+i2x,ifl(1),ifl(2),ifl(3),ifl(4))
c !print*,'secnd chosen hadron:',ispecs(i2),wzspecs(i2)
-
-c... generate momenta
+
+c... generate momenta
if(ish.ge.5)write(ifch,*)'hadron momenta:'
do n=nptlb+1,nptl
call idmass(id,am)
!prepare proposal function
! f_prop: f0(x)=0, x<am
- ! f1(x)=const=f2(b), am<x<b,
+ ! f1(x)=const=f2(b), am<x<b,
! f2(x)=x**2*exp(7-a*x), x>b
a=11
b=0.9
call utprix('hnbxxx',ish,ishini,4)
return
end
-
+
c----------------------------------------------------------------------
subroutine hnbxxxini
c----------------------------------------------------------------------
include 'epos.inc'
logical lcalc
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
!-------------------------------------------------------------
! 110, 120, -120, 130, -130, 230, -230, 220, 330
!, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
!, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330
!, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
- !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
+ !, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
!-------------------------------------------------------------
common/xxxspecs/wtot,wspecs(mspecs),zspecs(mspecs)
- integer ittspecs(mspecs)
+ integer ittspecs(mspecs)
parameter(mxdrop=35,mxe=10)
common/xxxspecsy/ndrop(-4:4,-4:4,-4:4)
common/xxxspecsx/ee(mxe),wwspecs(mxdrop,mxe,mspecs)
write(ifmt,'(3a)')'read from ',fndr(1:nfndr),' ...'
open(1,file=fndr(1:nfndr),status='old')
read(1,*)mxxdrop
- if(mxxdrop.ne.mxdrop)stop'hnbxxxini: wrong nr of droplets'
+ if(mxxdrop.ne.mxdrop)stop'hnbxxxini: wrong nr of droplets'
do n=1,mxdrop
read(1,*)ku,kd,ks
if(abs(ku).gt.4)stop'hnbxxxini: ku out of range'
else
stop'hnbxxxini: file not found. '
endif
- end
+ end
c----------------------------------------------------------------------
subroutine hnbaaa(ip,iret)
c----------------------------------------------------------------------
- include 'epos.inc'
- if(ioclude.eq.1)call hnbaaa156(ip,iret)
+ include 'epos.inc'
+ if(ioclude.eq.1)call hnbaaa156(ip,iret)
if(ioclude.eq.2)stop'ioclude.eq.2 no longer supported. '
if(ioclude.eq.3)call hnbaaanew(ip,iret)
end
-
+
c----------------------------------------------------------------------
subroutine hnbaaanew(ip,iret)
c----------------------------------------------------------------------
common/cen/ncentr /ctauhac/ntauhac(ncenthy,netahy)
common/cranphi/ranphi,ranecc,weiecc
character txt*40
- data icnthnb /0/ !vv2 /0./ nvv2 /0/ vv3 /0./
+ data icnthnb /0/ !vv2 /0./ nvv2 /0/ vv3 /0./
!save vv2,nvv2,vv3
save icnthnb
-
+
call utpri('hnbaaa',ish,ishini,4)
-
+
if(ish.ge.3)then
write(ifch,140)
140 format(/' ----------------------------------'/
write(ifch,*)'droplet:'
call alist('&',ip,ip)
endif
-
+
iret=0
do j=1,5
c(j)=pptl(j,ip)
enddo
-
+
call idquac(ip,nqi,nsi,nai,jc)
keu=jc(1,1)-jc(1,2)
ked=jc(2,1)-jc(2,2)
* 'id:',idptl(ip),' r:',radptl(ip),' m:',pptl(5,ip)
call utmsgf
endif
-
- !~~~~~~~~~read in freeze out surface properties from hydro~~~~~~~~~~~~
- if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
- icnthnb=icnthnb+1
+
+ !~~~~~~~~~read in freeze out surface properties from hydro~~~~~~~~~~~~
+ if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
+ icnthnb=icnthnb+1
if(icnthnb.eq.1)then
- !here we use epos.iniXXX (like epos.ini1fc)
+ !here we use epos.iniXXX (like epos.ini1fc)
!instead of the generic filename epos.inihy
open(3,file=fnnx(1:nfnnx)//fnhy(1:nfnhy),status='old',err=99)
read(3,*)txt
if(epscrix.ne.epscri(ioclude))then
print*,'hnbaaanew: epscrix.ne.epscri(',ioclude,'). '
stop
- endif
+ endif
read(3,*)ncenthyx,netahyx,ntauhyx,nphihyx,nradhyx
if(ncenthyx.ne.ncenthy)stop'hnbaaanew: ncenthyx.ne.ncenthy. '
if(netahyx.ne.netahy)stop'hnbaaanew: netahyx.ne.netahy. '
endif
endif
- !~~~~~~~~~define womi yomi romi~~~~~~~~~~~~
- if(iorsdf.eq.3.and.icnthnb.eq.1)then
+ !~~~~~~~~~define womi yomi romi~~~~~~~~~~~~
+ if(iorsdf.eq.3.and.icnthnb.eq.1)then
if(ioclude.eq.3)then
do ncent=1,ncenthy
do neta=1,netahy
elseif(ioclude.eq.2)then
stop'in hnbaaanew: ioclude=2 not supported any more.'
else
- stop'in hnbaaa: invalid ioclude. '
+ stop'in hnbaaa: invalid ioclude. '
endif
endif
!~~tau partition function paut(ncent,neta,ntau)
!~~phi partition function pauf(ncent,neta,ntau,nphi)
- if(iorsdf.eq.3.and.icnthnb.eq.1)then
+ if(iorsdf.eq.3.and.icnthnb.eq.1)then
do ncent=1,ncenthy
do neta=1,netahy
womax=0
ntauhac(ncent,neta)=0
- do ntau=1,ntauhoc(ncent)
- if(womi(ncent,neta,ntau,1).gt.womax )then
+ do ntau=1,ntauhoc(ncent)
+ if(womi(ncent,neta,ntau,1).gt.womax )then
womax=womi(ncent,neta,ntau,1)
ntauhac(ncent,neta)=ntau
endif
c1=cos(2*phihy(nphi-1))
c2=cos(2*phihy(nphi))
w1=womi(ncent,neta,ntau ,1)+c1*womi(ncent,neta,ntau ,2)
- . -womi(ncent,neta,ntau-1,1)-c1*womi(ncent,neta,ntau-1,2)
+ . -womi(ncent,neta,ntau-1,1)-c1*womi(ncent,neta,ntau-1,2)
w2=womi(ncent,neta,ntau ,1)+c2*womi(ncent,neta,ntau ,2)
- . -womi(ncent,neta,ntau-1,1)-c2*womi(ncent,neta,ntau-1,2)
+ . -womi(ncent,neta,ntau-1,1)-c2*womi(ncent,neta,ntau-1,2)
pauf(ncent,neta,ntau,nphi)
- . =pauf(ncent,neta,ntau,nphi-1)+0.5*(w1+w2)*dphi
+ . =pauf(ncent,neta,ntau,nphi-1)+0.5*(w1+w2)*dphi
enddo
w=pauf(ncent,neta,ntau,nphihy)
if(w.eq.0.)stop'hnbaaanew: w.eq.0. '
enddo
enddo
endif
-
+
!~~~~~~~~~determine ncentr~~~~~~~~~~~~
if(iorsdf.eq.3.and.ityptl(ip).eq.60)then !!!fusion on!!!
ncentr=0
if(db.lt.dbmin)then
dbmin=db
ncentr=ncent
- endif
+ endif
enddo
!print*,ncentr,bimevt,centhy(ncentr)
endif
- !~~~~~define masses~~~~~~~~~~~~~~~~
+ !~~~~~define masses~~~~~~~~~~~~~~~~
amin=utamnu(keu,ked,kes,kec,keb,ket,5)
aumin=amuseg
ipo=ip
tecmxx=tecm
!~~~~~~~~~determine netar~~~~~~~~~~~~
- if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
+ if(iorsdf.eq.3.and.ityptl(ip).eq.60)then
z=xorptl(3,ipo)
t=xorptl(4,ipo)
!print*,z,t,ityptl(ip),tecm
if(deta.lt.detamin)then
detamin=deta
netar=neta
- endif
+ endif
enddo
!print*,netar,zetaor,etahy(netar)
endif
fradflo=1.
-
- !~~~~~redefine energy in case of radial flow~~~~~~~~~~~~~~~~
+
+ !~~~~~redefine energy in case of radial flow~~~~~~~~~~~~~~~~
if(iappl.ne.4.and.iorsdf.eq.3.and.tecmor.gt.aumin
..and.ityptl(ip).eq.60)then
am=0.
fradflo=1.
if(en.ne.0.)fradflo=am/en
if(tecm*fradflo.lt.amin)fradflo=1.
- endif
+ endif
tecm=tecm*fradflo
-
- !~~~~~redefine energy in case of long coll flow
+
+ !~~~~~redefine energy in case of long coll flow
if(iappl.eq.4.or.iorsdf.ne.3
- &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
+ &.or.ityptl(ip).eq.40.or.ityptl(ip).eq.50)then !not for droplets from remnants
yco=0
else
if(ylongmx.lt.0.)then
else
yco=ylongmx
endif
- endif
+ endif
tecmx=tecm
if(yco.gt.0..and.tecmor.gt.aumin) then
- tecm=tecm/sinh(yco)*yco
- else
- yco=0.
+ tecm=tecm/sinh(yco)*yco
+ else
+ yco=0.
endif
!print*,'========= cluster energy: ',pptl(5,ip),tecmx,tecm
!~~~~~~~~~redefine volume~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- vocri=tecm/epscri(ioclude)
- volu=max(vocri,vocell)
+ vocri=tecm/epscri(ioclude)
+ volu=max(vocri,vocell)
!~~~~~~~~~decay~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
call hnbini(iret)
!if(iret.ne.0)write(ifch,*)'***** unsucessfull hnbini *****'
if(iret.ne.0)goto1000
if(ioinct.ge.1)goto1
-
+
do iter=1,itermx
naccit(iter)=0
call hnbmet
enddo
-
+
1 continue
if(ioceau.eq.1.and.iappl.eq.1)call xhnbte(ip)
do j=1,4
pe(j)=0.
enddo
- do i=1,np
+ do i=1,np
uu(1)= 0
uu(2)= 0
uu(3)= sinh(yrad(i))
do j=1,4
pa(j)=0.
enddo
- do i=1,np
+ do i=1,np
call utlob3(1,pe(1),pe(2),pe(3),pe(4),pe(5)
* , pcm(1,i), pcm(2,i), pcm(3,i), pcm(4,i))
do j=1,4
pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
* +amass(j)**2)
sum=sum+pcm(4,j)
- enddo
+ enddo
scal=esoll/sum
!write(6,*)'ipass,scal,e,esoll:'
- ! $ ,ipass,scal,sum,esoll
+ ! $ ,ipass,scal,sum,esoll
if(abs(scal-1.).le.errlim) goto301
enddo
- 301 continue
+ 301 continue
do j=1,4
pa(j)=0.
enddo
- do i=1,np
+ do i=1,np
do j=1,4
pa(j)=pa(j)+pcm(j,i)
enddo
bb=0
cc=0
dd=1
- if(ityptl(ip).eq.60)then
+ if(ityptl(ip).eq.60)then
ipo=iorptl(ip)
xx=uptl(ipo) ! <x**2>
yy=optl(ipo) ! <y**2>
endif
errlim=0.0001
tecm=tecmxx
- phinull=phievt+ranphi
+ phinull=phievt+ranphi
do n=1,np
!~~determine random tau from paut(ncentr,netar,ntau)
r=rangen()
f=f1/(f1-f2)
fx=f
tau= tau1*(1-f) + tau2*f
- taufop(n)=tau
+ taufop(n)=tau
!~~determine phifop~~~~~
r=rangen()
if(pauf(ncentr,netar,ntau-1,nphihy).gt.0.
endif
phifop(n)=phi
!~~determine yrad~~~~
- yr1=yomi(ncentr,netar,ntau1,1)
- . +yomi(ncentr,netar,ntau1,2)*cos(2*phi)
- yr2=yomi(ncentr,netar,ntau2,1)
- . +yomi(ncentr,netar,ntau2,2)*cos(2*phi)
- yr= yr1*(1-fx)+ fx*yr2
- yrad(n)=yr
+ yr1=yomi(ncentr,netar,ntau1,1)
+ . +yomi(ncentr,netar,ntau1,2)*cos(2*phi)
+ yr2=yomi(ncentr,netar,ntau2,1)
+ . +yomi(ncentr,netar,ntau2,2)*cos(2*phi)
+ yr= yr1*(1-fx)+ fx*yr2
+ yrad(n)=yr
!~~determine radfop~~~~
rad1=romi(ncentr,netar,ntau1,1)
. +romi(ncentr,netar,ntau1,2)*cos(2*phi)
radfop(n)=rad
enddo
energ=0.
- do n=1,np
- uu(1)=sinh(yrad(n))*cos(phifop(n)+phinull)
- uu(2)=sinh(yrad(n))*sin(phifop(n)+phinull)
- uu(3)=0d0
+ do n=1,np
+ uu(1)=sinh(yrad(n))*cos(phifop(n)+phinull)
+ uu(2)=sinh(yrad(n))*sin(phifop(n)+phinull)
+ uu(3)=0d0
uu(4)=sqrt(1+uu(1)**2+uu(2)**2)
!px=pcm(1,n)
!py=pcm(2,n)
pcm(4,j)=sqrt(pcm(1,j)**2+pcm(2,j)**2+pcm(3,j)**2
* +amass(j)**2)
sum=sum+pcm(4,j)
- enddo
+ enddo
scal=esoll/sum
!write(6,*)'ipass,scal,e,esoll:'
- ! $ ,ipass,scal,sum,esoll
+ ! $ ,ipass,scal,sum,esoll
if(abs(scal-1.).le.errlim) goto300
enddo
- 300 continue
+ 300 continue
!print*, scal
else
do n=1,np
!zeta=zetaor-0.5*delzet+delzet*rangen()
z=tau*sinh(zeta)
t=tau*cosh(zeta)
- xorptl(1,nptl)=r*cos(phifop(n)+phinull)
- xorptl(2,nptl)=r*sin(phifop(n)+phinull)
+ xorptl(1,nptl)=r*cos(phifop(n)+phinull)
+ xorptl(2,nptl)=r*sin(phifop(n)+phinull)
xorptl(3,nptl)=z
xorptl(4,nptl)=t
else
xorptl(2,nptl)=xorptl(2,ip)
xorptl(3,nptl)=xorptl(3,ip)
xorptl(4,nptl)=xorptl(4,ip)
- endif
+ endif
endif
- enddo
+ enddo
if(ish.ge.3)then
write(ifch,*)'decay products:'
call alist('&',nptlb+1,nptl)
call utprix('hnbaaa',ish,ishini,4)
return
- 99 print*,'hnbaaanew: error opening hydro table'
+ 99 print*,'hnbaaanew: error opening hydro table'
print*,' file=',fnnx(1:nfnnx)//fnhy(1:nfnhy)
print*,'maybe "fname inihy ..." forgotten in optns file ???'
print*,' this is necessary in case of "set ioclude 3"'
stop'070817'
-
+
end
c------------------------------------------------------------------------------
if(iSpaceTime.eq.1.and.ioclude.gt.1)then
call xCoreCorona(0,0)
do neta=1,5,2
- call xFoMass(neta)
- call xFoRadius(neta)
- call xFoRadRapidity(neta)
- call xFreezeOutTauX(neta)
+ call xFoMass(neta)
+ call xFoRadius(neta)
+ call xFoRadRapidity(neta)
+ call xFreezeOutTauX(neta)
enddo
call xFreezeOutTauEta
call xFreezeOutTZ
elseif(iSpaceTime.eq.1)then
call xCoreCorona(0,0)
!stop'bjinta: space-time plots require ioclude>1. '
- endif
+ endif
end
c------------------------------------------------------------------------------
npl=0
endif
endif
- endif
+ endif
endif
- endif
- endif
+ endif
+ endif
enddo
if(nplx.gt.0)write(ifhi,'(a)') ' endarray closehisto plot 0-'
!..........................................................................
npl=0
endif
endif
- endif
- endif
+ endif
+ endif
enddo
write(ifhi,'(a)') ' endarray closehisto plot 0'
!..........................................................................
end
-
+
c------------------------------------------------------------------------------
subroutine xFreezeOutTauEta
c------------------------------------------------------------------------------
write(ifhi,'(a,f5.1)') 'yrange 0 ',taumax
write(ifhi,'(a)') 'txt "xaxis [c] "'
write(ifhi,'(a)') 'txt "yaxis [t] (fm/c)"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.77 0.92 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.77 0.92 "'//cbim//'"'
write(ifhi,'(a)') 'array 2'
endif
eta=
npl=0
endif
endif
- endif
- endif
+ endif
+ endif
enddo
if(npl.ne.0)write(ifhi,'(a)') ' endarray closehisto plot 0'
if(npl.eq.0)stop'xFreezeOutTZ: no particles!!!!! '
write(ifhi,'(a)') 'yrange 0 25 '
write(ifhi,'(a)') 'txt "xaxis z (fm)"'
write(ifhi,'(a)') 'txt "yaxis t (fm/c)"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.22 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.22 "'//cbim//'"'
write(ifhi,'(a)') 'array 2'
endif
write(ifhi,'(2e11.3)') xorptl(3,n),xorptl(4,n)
nhis=nhis+1
npl=0
endif
- endif
- endif
+ endif
+ endif
enddo
if(npl.ne.0)write(ifhi,'(a)') ' endarray closehisto plot 0'
if(npl.eq.0)stop'xFreezeOutTZ: no particles!!!!! '
do ii=1,2
write(ifhi,'(a,3i1)')'openhisto htyp lin name w-',neta,nn,ii
if(ii.eq.1)then !----------------------
- write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
- write(ifhi,'(a)') 'txt "xaxis [t] (fm/c)"'
+ write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
+ write(ifhi,'(a)') 'txt "xaxis [t] (fm/c)"'
write(ifhi,'(a)') 'ymod lin yrange auto auto '
write(ifhi,'(a,f4.2,a)') 'text 0.1 0.9 " [c]=',etahy(neta),'"'
write(ifhi,'(a)')'txt "yaxis w?0! w?2! (GeVc/fm) "'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
endif !-------------------------------
write(ifhi,'(a)')'array 2'
deltau=tauhoc(ncentr,2)-tauhoc(ncentr,1)
write(ifhi,'(a)') 'endarray closehisto '
if(ii.ne.2)write(ifhi,'(a)') 'plot 0-'
if(ii.eq.2)write(ifhi,'(a)') 'plot 0'
- enddo
+ enddo
end
c------------------------------------------------------------------------------
do ii=1,2
write(ifhi,'(a,3i1)')'openhisto htyp lin name r-',neta,nn,ii
if(ii.eq.1)then !----------------------
- write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
+ write(ifhi,'(a,f4.1)')'xmod lin xrange 0. ',taumax
write(ifhi,'(a)')'txt "xaxis [t] (fm/c)"'
write(ifhi,'(a)') 'ymod lin yrange auto auto '
write(ifhi,'(a,f4.2,a)')'text 0.1 0.9 " [c]=',etahy(neta),'"'
write(ifhi,'(a)')'txt "yaxis r?0! r?2! (fm) "'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
endif !-------------------------------
write(ifhi,'(a)')'array 2'
do ntau=2,ntauhac(ncentr,neta)
if(ii.eq.2)write(ifhi,'(a)') 'plot 0'
enddo
end
-
+
c------------------------------------------------------------------------------
subroutine xFoRadRapidity(neta)
c------------------------------------------------------------------------------
write(ifhi,'(a)') 'ymod lin yrange auto auto '
write(ifhi,'(a,f4.2,a)')'text 0.1 0.9 " [c]=',etahy(neta),'"'
write(ifhi,'(a)')'txt "yaxis y?0! y?2! "'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.50 0.9 "'//cbim//' "'
endif !-------------------------------
write(ifhi,'(a)')'array 2'
do ntau=2,ntauhac(ncentr,neta)
& ,w1define(mxdefine),w2define(mxdefine)
data iperc /0,5,10,15,20,25,30,35,40,45,50,60,70,80,92,100/
-
+
do n=1,maxb
bim(n)=0
enddo
read(w2define(m)(1:l2define(m)),*)bim(14)
elseif(w1define(m)(1:l1define(m)).eq.'bim92')then
read(w2define(m)(1:l2define(m)),*)bim(15)
- endif
+ endif
enddo
bim(16)=100.
do n=2,maxb
print*,'******* ERROR in subroutine centrality: '
. ,' #define bim?? ??? missing. ******* '
print*,' n=',n
- stop'14082007'
+ stop'14082007'
endif
enddo
if(abs(ch).gt.0.1.and.abs(rap).le.1.)multy1=multy1+1
endif
endif
- enddo
+ enddo
ih1=jjj/100
ih2=mod(jjj,100)
if(0.5*multy1.lt.ih1.or.0.5*multy1.gt.ih2)return
write(ifhi,'(a)') 'text 0.05 0.90 "core"'
write(ifhi,'(a)') 'text 0.05 0.80 "corona"'
write(ifhi,'(a,f4.1,a)')'text 0.82 0.07 "[c]=0"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.77 0.92 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.77 0.92 "'//cbim//'"'
write(ifhi,'(a)') 'array 2'
ncore=0
do i=1,nptl
do j=-50,50
phi=j/50.*pi*0.55
write(ifhi,'(2e11.3)')r1*cos(phi)-b,r1*sin(phi)
- enddo
+ enddo
write(ifhi,'(a)') ' endarray'
write(ifhi,'(a)') 'closehisto'
- endif
+ endif
if(r2.ne.0.0)then
write(ifhi,'(a)') 'plot 0-'
write(ifhi,'(a)') 'openhisto name stc1 htyp lyu'
do j=-50,50
phi=j/50.*pi*0.55
write(ifhi,'(2e11.3)')-r1*cos(phi)+b,r1*sin(phi)
- enddo
+ enddo
write(ifhi,'(a)') ' endarray'
write(ifhi,'(a)') 'closehisto'
- endif
-
+ endif
+
write(ifhi,'(a)') 'plot 0'
!........................................................................................
m=(rapx+rapmax)/delrap+1
if(m.gt.myy)m=myy
yy(m)=yy(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! core segment energy per d[c]dxdy '
write(ifhi,'(a)')'! vs space-time rapidity rapx '
write(ifhi,'(a)') 'txt "title initial energy "'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "x=0"'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "y=0"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
write(ifhi,'(a)') 'txt "xaxis space-time rapidity [c] "'
write(ifhi,'(a)') 'txt "yaxis dE/d[c]dxdy "'
write(ifhi,'(a)') 'array 2'
endif
m=(rapx+rapmax)/delrap+1
if(m.ge.1.and.m.le.myy)yy(m)=yy(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! corona segment energy per d[c]dxdy '
write(ifhi,'(a)')'! vs space-time rapidity rapx '
write(ifhi,'(2e11.3)')-rapmax+(m-0.5)*delrap, yy(m)/4./delrap
enddo
write(ifhi,'(a)') ' endarray closehisto plot 0-'
- call xEiniEta(1)
+ call xEiniEta(1)
write(ifhi,'(a)') 'plot 0'
!........................................................................................
delrad=2*radmax/float(mrr)
m=(rinp+radmax)/delrad+1
if(m.gt.mrr)m=mrr
rr(m)=rr(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! core segment energy per d[c]dxdy vs x '
write(ifhi,'(a)')'! (same as histogram rinp eco... in optns) '
write(ifhi,'(a)') 'txt "title initial energy "'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "[c]=0"'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "y=0"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
write(ifhi,'(a)') 'txt "xaxis x (fm)"'
write(ifhi,'(a)') 'txt "yaxis dE/d[c]dxdy "'
write(ifhi,'(a)') 'array 2'
m=(rinp+radmax)/delrad+1
if(m.gt.mrr)m=mrr
rr(m)=rr(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! corona segment energy per d[c]dxdy vs x '
write(ifhi,'(a)')'! (same as histogram rinp eco... in optns) '
write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
enddo
write(ifhi,'(a)') ' endarray closehisto plot 0-'
- call xEiniX(1)
+ call xEiniX(1)
write(ifhi,'(a)') 'plot 0'
!........................................................................................
delrad=2*radmax/float(mrr)
m=(routp+radmax)/delrad+1
if(m.gt.mrr)m=mrr
rr(m)=rr(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! core segment energy per d[c]dxdy vs y '
write(ifhi,'(a)')'! (same as histogram routp eco... in optns) '
write(ifhi,'(a)') 'txt "title initial energy "'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.70 "[c]=0"'
write(ifhi,'(a,f4.1,a)')'text 0.05 0.60 "x=0"'
- write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
+ write(ifhi,'(a,f4.1,a,f4.1,a)')'text 0.75 0.90 "'//cbim//'"'
write(ifhi,'(a)') 'txt "xaxis y (fm)"'
write(ifhi,'(a)') 'txt "yaxis dE/d[c]dxdy "'
write(ifhi,'(a)') 'array 2'
m=(routp+radmax)/delrad+1
if(m.gt.mrr)m=mrr
rr(m)=rr(m)+eco
- endif
- endif
- enddo
+ endif
+ endif
+ enddo
write(ifhi,'(a)')'!---------------------------------------------'
write(ifhi,'(a)')'! corona segment energy per d[c]dxdy vs y '
write(ifhi,'(a)')'! (same as histogram routp eco... in optns) '
write(ifhi,'(2e11.3)')-radmax+(m-0.5)*delrad, rr(m)/4./delrad
enddo
write(ifhi,'(a)') ' endarray closehisto plot 0-'
- call xEiniY(1)
+ call xEiniY(1)
write(ifhi,'(a)') 'plot 0'
end
enddo
write(ifhi,'(a)') 'endarray closehisto '
return
-
+
entry xEiniY(ii)
write(ifhi,'(a)') '!----------------------------------------'
write(ifhi,'(a,i3)') '! hydro initial energy vs y '
c------------------------------------------------------------------------------
subroutine hnbcor(mode)
c------------------------------------------------------------------------------
-c determines(mode=1) and plots (mode=2) two particle correlations
+c determines(mode=1) and plots (mode=2) two particle correlations
c for the configurations /confg/
c------------------------------------------------------------------------------
include 'epos.inc'
if(abs(cs).gt.1.)then
cs=aint(cs)
- ang=acos(cs)
+ ang=acos(cs)
else
ang=acos(cs)
endif
nk=1
nw=bns
else
- nw=1+aint(ang/pi*bns)
- nk=1+aint((cs+1.)/2.*bns)
+ nw=1+aint(ang/pi*bns)
+ nk=1+aint((cs+1.)/2.*bns)
endif
nw=min(nw,bns)
nk=min(nk,bns)
flog=flog+alog(gg*am*volu/4/pi**3/hquer**3/(nlattc+1-i))
enddo
faclog=faclog+flog
-
+
return
end
end
cc----------------------------------------------------------------------
-c subroutine hnbids(jc,ids,iwts,i)
+c subroutine hnbids(jc,ids,iwts,i)
cc----------------------------------------------------------------------
cc returns i id-codes ids() corr to jc and their weights iwts()
cc----------------------------------------------------------------------
c i=i+1
c ids(i)=0
c iwts(i)=iozero
-c 1 continue
+c 1 continue
c
c do j=1,nspecs
c do n=1,nflav
c i=i+1
c ids(i)=ispecs(j)
c iwts(i)=1
-c 2 continue
+c 2 continue
c enddo
c
c return
subroutine hnbiiw(x,f,df)
c----------------------------------------------------------------------
c returns fctn value and first derivative at x of the
-c i-th integrated weight fctn minus random number
+c i-th integrated weight fctn minus random number
c for the asympotic phase space integral.
c input:
c x: x-value
include 'epos.inc'
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/crnoz/rnoz(maxp-1)
common/citer/iter,itermx
- common/cfact/faclog
+ common/cfact/faclog
common/chnbin/nump,ihadro(maxp)
common /clatt/nlattc,npmax
parameter(maxit=50000)
call hnbmin(keu,ked,kes,kec)
elseif(tecm.lt.2.0.and.nk.ne.0)then
call hnbmin(keu,ked,kes,kec)
- else
+ else
call hgcaaa
call hgcnbi(iret)
if(iret.eq.1)then
c *'n:',n,'weight:',wt,'wts/n:',wts/n,'error:',wts/n/sqrt(1.*n)
c enddo
c return
-c end
+c end
cc----------------------------------------------------------------------
subroutine hnbmet
c----------------------------------------------------------------------
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
common/crnoz/rnoz(maxp-1)
real rnozo(maxp-1)
- common/cfact/faclog
+ common/cfact/faclog
dimension amasso(maxp),idento(maxp),pcmo(5,maxp)
integer jc(nflav,2),jc1(nflav,2),jc2(nflav,2)
common/citer/iter,itermx
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (literm=500)
common/cmet/kspecs(mspecs),liter,lspecs(literm,mspecs)
enddo
call hnbzmu(-1)
endif
-
+
c remember old configuration
c --------------------------
wtlo=wtlog
*call utstop('hnbmet: invalid particle species&')
enddo
keepr=0
-c-c call hnbolo(1000) !instead of "call hnbody" for testing
+c-c call hnbolo(1000) !instead of "call hnbody" for testing
keepr=1
if(iocova.eq.1)call hnbody
if(iocova.eq.2)call hnbodz
if(wtlog-wtlo.lt.30.)then
q=exp(wtlog-wtlo)*xba/xab
r=rangen()
- if(r.le.q)iacc=1
+ if(r.le.q)iacc=1
if(ish.ge.7)write(ifch,*)'new weight / old weight:',q,' '
*,'random number:',r
else
if(ish.ge.7)write(ifch,*)'log new weight / old weight:'
*,wtlog-wtlo
endif
- if(iacc.eq.1)then
+ if(iacc.eq.1)then
if(ish.ge.7)write(ifch,*)'new configuration accepted'
nacc=nacc+1
naccit(iter)=1
c do j=1,nspecs
c lspecs(liter,j)=lspecs(liter-1,j)
c enddo
-c endif
+c endif
endif
endif
if(iter.le.iternc)return
nump=0
f1='(4i3,i7,i6)'
ke=iabs(keux+kedx+kesx+kecx)
-
+
if(keux+kedx+kesx+kecx.ge.0)then
keu=keux
ked=kedx
endif
if(wri)write(ifch,'(4i3)')keux,kedx,kesx,kecx
if(wri)write(ifch,'(4i3)')keu,ked,kes,kec
-
+
c get rid of anti-c and c (140, 240, -140, -240)
if(kec.ne.0)then
10 continue
goto11
endif
endif
-
+
c get rid of anti-s (130,230)
5 continue
if(kes.lt.0)then
goto5
endif
-c get rid of anti-d (120, -230)
+c get rid of anti-d (120, -230)
6 continue
if(ked.lt.0)then
ked=ked+1
endif
goto7
endif
-
+
if(keu+ked+kes+kec.ne.ke)call utstop('hnbmin: sum_kei /= ke&')
keq=keu+ked
-c get rid of s (3331, x330, xx30)
+c get rid of s (3331, x330, xx30)
i=4
2 i=i-1
3 continue
if(keu+ked.ne.keq)call utstop('hnbmin: keu+ked /= keq&')
-c get rid of d (2221, 1220, 1120)
+c get rid of d (2221, 1220, 1120)
i=4
12 i=i-1
13 continue
if(ked.ne.0)call utstop('hnbmin: ked .ne. 0&')
-c get rid of u (1111)
+c get rid of u (1111)
9 continue
if(keu.gt.0)then
keu=keu-3
c-------------------------------------------------------------
subroutine hnbody
-c-------------------------------------------------------------
+c-------------------------------------------------------------
c formerly subr genbod from genlib (cernlib).
c modified by K. Werner, march 94.
c subr to generate n-body event
c according to fermi lorentz-invariant phase space.
-c the phase space integral is the sum over the weights wt divided
+c the phase space integral is the sum over the weights wt divided
c by the number of events (sum wt / n).
c adapted from fowl (cern w505) sept. 1974 by f. james.
c events are generated in their own center-of-mass,
c
c input to and output from subr thru common block config.
c input:
-c np=number of outgoing particles
+c np=number of outgoing particles
c tecm=total energy in center-of-mass
c amass(i)=mass of ith outgoing particle
c output:
external hnbiiw
ctp060829 nas=5 !must be at least 3
wri=.false.
- if(ish.ge.7)wri=.true.
+ if(ish.ge.7)wri=.true.
if(wri)then
write(ifch,*)('-',i=1,10)
*,' entry sr hnbody ',('-',i=1,30)
if(ntm2) 9,5,4
4 continue
call flpsore(rno,ntm2)
-
+
c...calculate emm().......M_i
do 6 j=2,ntm1
6 emm(j)=rno(j-1)*tecmtm+sm(j)
-
+
c...calculate wtlog
5 continue
pcm(5,j)=sqrt(aa)
pcm(4,j)=sqrt(aa+ems(j))
call hnbrt2(c,s,cb,sb,pcm,j)
- enddo
+ enddo
endif
enddo
R=J
GO TO 20
50 continue
-
+
do i=1,n-1
if(a(i).gt.a(i+1))stop'FLPSORE: ERROR. '
- enddo
-
+ enddo
+
RETURN
END
c-------------------------------------------------------------
subroutine hnbodz
-c-------------------------------------------------------------
+c-------------------------------------------------------------
c subr to generate n-body event
c according to non-invariant phase space.
c the phase space integral is the sum over the weights exp(wtxlog)
c
c input to and output from subr is thru common block config.
c input:
-c np=number of outgoing particles
+c np=number of outgoing particles
c tecm=total energy in center-of-mass
c amass(i)=mass of ith outgoing particle
c output:
if(ish.ge.6)write(ifch,1200)np,tecm
if(ish.ge.6)write(ifch,*)'particle masses:'
if(ish.ge.6)write(ifch,'(1x,10f6.3)')(amass(n),n=1,np)
-
+
c initialization ktnbod=1
ktnbod=ktnbod + 1
if(ktnbod.gt.1) goto 1
ffqlog(n)=ffqlog(n-1)+alog(4*pi/(n-1))
enddo
1 continue
-c set wtxlog -infinity for np<2
+c set wtxlog -infinity for np<2
if(np.lt.2) goto 1001
c special treatment for np=2
if(np.eq.2)then
wtxlog=alog(tt)*(np-1) + ffqlog(np)
if(ish.ge.7)
*write(ifch,*)'wtxlog:',wtxlog,' (prefactor)'
-c fill rnoz with np-1 random numbers
+c fill rnoz with np-1 random numbers
if(keepr.eq.0)then
do 3 i= 1, np-1
3 rnoz(i)=rangen()
pcm(3,i)=0
pcm(4,i)=ti(i)+amass(i)
p52=ti(i)*(ti(i)+2*amass(i))
- if(p52.gt.0)then
+ if(p52.gt.0)then
pcm(5,i)=sqrt(p52)
else
pcm(5,i)=ti(i)*sqrt(1+2*amass(i)/ti(i))
endif
if(w.le.0.)goto1111
c complete specification of event (random rotations and then deformations)
- call hnbrot
+ call hnbrot
if(ish.ge.7)write(ifch,*)'momenta after rotations:'
call hnbrop(96,0)
call hnbrod
include 'epos.inc'
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
- a=0
- k=0
+ a=0
+ k=0
do j=1,loops
c-c if(mod(j,iterpr).eq.0)write(ifmt,*)' iteration:',iter,j
if(iocova.eq.1)call hnbody
if(iocova.eq.2)call hnbodz
if(ish.ge.8)write(ifch,*)'j:',j,' wtxlog:',wtxlog
- if(wtxlog.gt.-1e30)then
- k=k+1
+ if(wtxlog.gt.-1e30)then
+ k=k+1
if(k.eq.1)c=wtxlog
if(a.gt.0.)then
if(alog(a).lt.wtxlog-c-20)then
c=wtxlog
endif
endif
- a=a+exp(wtxlog-c)
- endif
+ a=a+exp(wtxlog-c)
+ endif
if(ish.ge.8)write(ifch,*)'k:',k,' c:',c
- enddo
+ enddo
a=a/loops
- wtxlog=alog(a)+c
+ wtxlog=alog(a)+c
return
end
c returns momentum p for twobody decay a --> b + c
c a, b, c are the three masses
c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-c this p is related to twobody phase space as R2 = pi * p /a
+c this p is related to twobody phase space as R2 = pi * p /a
c-----------------------------------------------------------------------
double precision aa,bb,cc,a2,b2,c2
aa=a
if(n2.lt.n1)then
n1r=n1
n1=n2
- n2=n1r
+ n2=n1r
endif
if(k.eq.2)then
if(n1.eq.k1.or.n1.eq.k2.or.n2.eq.k1.or.n2.eq.k2)goto1
endif
- if(ident(n1).ne.0.and.ident(n2).ne.0)mm=1 ! hadron-hadron
- if(ident(n1).ne.0.and.ident(n2).eq.0)mm=2 ! hadron-empty
- if(ident(n1).eq.0.and.ident(n2).ne.0)mm=2 ! empty-hadron
- if(ident(n1).eq.0.and.ident(n2).eq.0)mm=3 ! empty-empty
+ if(ident(n1).ne.0.and.ident(n2).ne.0)mm=1 ! hadron-hadron
+ if(ident(n1).ne.0.and.ident(n2).eq.0)mm=2 ! hadron-empty
+ if(ident(n1).eq.0.and.ident(n2).ne.0)mm=2 ! empty-hadron
+ if(ident(n1).eq.0.and.ident(n2).eq.0)mm=3 ! empty-empty
if(ish.ge.7)then
write(ifch,'(a,i2)')' mm:',mm
write(ifch,*)'to be replaced:',n1,ident(n1)
write(ifch,*)'to be replaced:',n2,ident(n2)
endif
-c flavour of n1+n2 --> jc
+c flavour of n1+n2 --> jc
c -----------------------
- if(mm.eq.1)then
+ if(mm.eq.1)then
call idtr4(ident(n1),ic1)
call iddeco(ic1,jc1)
call idtr4(ident(n2),ic2)
do j=1,2
jc(i,j)=0
enddo
- enddo
+ enddo
endif
if(k.eq.2)then
c----------------------------------------------------------------------
subroutine hnbpai(id1,id2,jc)
c----------------------------------------------------------------------
-c returns arbitrary hadron pair id1,id2, flavour written to jc
+c returns arbitrary hadron pair id1,id2, flavour written to jc
c----------------------------------------------------------------------
include 'epos.inc'
integer jc(nflav,2),jc1(nflav,2),ic1(2),jc2(nflav,2),ic2(2)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
c construct pair id1,id2
c construct possible pairs id1,id2
c --------------------------------
-
+
ipair=0
iwpair=0
idx=0
if(abs(jcmi(1)).gt.3)goto3
if(abs(jcmi(2)).gt.3)goto3
if(abs(jcmi(3)).gt.3)goto3
- if(abs(jcmi(4)).gt.3)goto3 !-charm
+ if(abs(jcmi(4)).gt.3)goto3 !-charm
if(jcmi(1).ne.0)goto111
if(jcmi(2).ne.0)goto111
if(jcmi(3).ne.0)goto111
- if(jcmi(4).ne.0)goto111 !-charm
+ if(jcmi(4).ne.0)goto111 !-charm
nids=nids+1
ids(nids)=0
iwts(nids)=iozero
111 continue
-
+
lkfok1=lkfok(1,jcmi(1),jcmi(2),jcmi(3),jcmi(4))
if(lkfok1.gt.0)then
nids=nids+1
subroutine hnbpajini
c----------------------------------------------------------------------
c initialize array to speed up hnbpaj calculation
-c store sum of weights iwpair of possible pairs in an array
+c store sum of weights iwpair of possible pairs in an array
c for any combinations of quarks
c----------------------------------------------------------------------
include 'epos.inc'
do iqs=0,2
do iqd=0,2
do iqu=0,2
-
+
idx=idx+1
idxpair(iqu,iqd,iqs,iaqu,iaqd,iaqs)=idx
ids(nids)=0
iwts(nids)=iozero
111 continue
-
+
lkfok1=lkfok(1,jcmi(1),jcmi(2),jcmi(3),jcmi(4))
if(lkfok1.gt.0)then
nids=nids+1
c input: dimension np and momenta p_i=pcm(5,i) via /confg/
c 1 < np <= npx : hagedorn method
c npx < np <= npy : integral method
-c npy < np : asymptotic method
+c npy < np : asymptotic method
c--------------------------------------------------------------------
include 'epos.inc'
parameter(maxp=500)
integer ii(maxp),isi(maxp)
double precision ppcm(maxp),ww,ppsum,ppmax
external hnbrax
- common/cepsr/nepsr
+ common/cepsr/nepsr
if(ish.ge.9)write(ifch,*)('-',i=1,10)
*,' entry sr hnbraw ',('-',i=1,30)
if(ish.ge.7)write(ifch,'(1x,a,e12.5,4x)')
*'sum p_i - 2*p_max not positive --> w:',w
goto1000
- endif
+ endif
c asymptotic method
c -----------------
if(ish.ge.9)write(ifch,*)'it:',it
b=b*5./3.
wio=win
- call uttrap(hnbrax,0.,b,win)
+ call uttrap(hnbrax,0.,b,win)
iok=0
if(abs(win-wio).le.epsr*abs((win+wio)/2))iok=1
- if(it.eq.itmax)iok=1
+ if(it.eq.itmax)iok=1
if(ish.ge.8.or.ish.ge.7.and.iok.eq.1)
*write(ifch,'(1x,2(a,e12.5,2x),a,i2,2x,a,i4)')
*'integral method: win:',win
ww=ww*pmax/ppcm(i)/2./i
enddo
ww=-ww/pmax**3/pi/2.*np*(np-1)*(np-2)
- whd=ww
+ whd=ww
if(ish.ge.7)write(ifch,'(1x,a,e12.5,4x,a)')
*'hagedorn method: whd:',whd,'double precision'
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
c integer identx(maxp)
- common /clatt/nlattc,npmax
+ common /clatt/nlattc,npmax
if(ish.ge.9)write(ifch,*)('-',i=1,10)
*,' entry sr hnbrmz ',('-',i=1,30)
if(np.eq.0)goto1000
endif
if(i.eq.np+1)goto1000
-
+
ident(i)=ident(np)
ident(np)=0
goto1
c input: pcm(5,i)
c output: pcm(1-3,i)
c----------------------------------------------------------------------
- common/cnsta/pi,pii,hquer,prom,piom,ainfin
+ common/cnsta/pi,pii,hquer,prom,piom,ainfin
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
real u(3)
cc-------------------------------------------------------------------
c subroutine hnbrt2old(c,s,c2,s2,pr,i)
cc-------------------------------------------------------------------
-cc formerly subr rotes2 from cernlib
+cc formerly subr rotes2 from cernlib
cc this subr now does two rotations (xy and xz)
cc-------------------------------------------------------------------
c parameter(maxp=500)
c-------------------------------------------------------------------
subroutine hnbrt2(c,s,c2,s2,pr,i)
c-------------------------------------------------------------------
-c formerly subr rotes2 from cernlib
+c formerly subr rotes2 from cernlib
c this subr now does two rotations (xy and xz)
c-------------------------------------------------------------------
parameter(maxp=500)
c defines particle species and masses and degeneracies.
c input:
c iopt=odd number: massless
-c iopt=even number: same as iopt-1, but massive
+c iopt=even number: same as iopt-1, but massive
c iopt= 1: pi0 (massless)
c iopt= 2: pi0
c iopt= 3: pi-,pi0,pi+ (massless)
c iopt= 8: 25 hadrons
c iopt= 9: 54 hadrons (massless)
c iopt=10: 54 hadrons
-c iopt=11: 3 quarks (massless)
-c iopt=12: 3 quarks
-c iopt=13: 54 hadrons + J/psi (massless)
-c iopt=14: 54 hadrons + J/psi
-c iopt=15: 54 hadrons + J/psi + H (massless)
-c iopt=16: 54 hadrons + J/psi + H
+c iopt=11: 3 quarks (massless)
+c iopt=12: 3 quarks
+c iopt=13: 54 hadrons + J/psi (massless)
+c iopt=14: 54 hadrons + J/psi
+c iopt=15: 54 hadrons + J/psi + H (massless)
+c iopt=16: 54 hadrons + J/psi + H
c output:
c nspecs: nr of species
c ispecs: id's
c aspecs: masses
c-----------------------------------------------------------------------
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nflav=6)
integer jc(nflav,2),ic(2)
data jspe01/ 110 /
data jspe03/ 110, 120, -120 /
data jspe05/ 110, 120, -120, 1120,-1120, 1220,-1220 /
- data jspe07/
+ data jspe07/
* 110, 120, -120, 130, -130, 230, -230, 220, 330
*, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
*, 1230,-1230, 2230,-2230, 1330,-1330, 2330,-2330 /
- data jspe09/
+ data jspe09/
* 110, 120, -120, 130, -130, 230, -230, 220, 330
*, 111, 121, -121, 131, -131, 231, -231, 221, 331
*, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
*, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331 /
data jspe11/
* 1, -1, 2, -2, 3, -3 /
- data jspe13/
+ data jspe13/
* 110, 120, -120, 130, -130, 230, -230, 220, 330
*, 111, 121, -121, 131, -131, 231, -231, 221, 331
*, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
*, 1111,-1111, 1121,-1121, 1221,-1221, 2221,-2221, 1131,-1131
*, 1231,-1231, 2231,-2231, 1331,-1331, 2331,-2331, 3331,-3331
*, 441 /
- data jspe15/
+ data jspe15/
* 110, 120, -120, 130, -130, 230, -230, 220, 330
*, 111, 121, -121, 131, -131, 231, -231, 221, 331
*, 1120,-1120, 1220,-1220, 1130,-1130, 2130,-2130
if(ioptx.eq.13)nspecs=nspe13
if(ioptx.eq.15)nspecs=nspe15
do i=1,nspecs
- if(ioptx.eq.1)ispecs(i)=jspe01(i)
- if(ioptx.eq.3)ispecs(i)=jspe03(i)
- if(ioptx.eq.5)ispecs(i)=jspe05(i)
- if(ioptx.eq.7)ispecs(i)=jspe07(i)
- if(ioptx.eq.9)ispecs(i)=jspe09(i)
+ if(ioptx.eq.1)ispecs(i)=jspe01(i)
+ if(ioptx.eq.3)ispecs(i)=jspe03(i)
+ if(ioptx.eq.5)ispecs(i)=jspe05(i)
+ if(ioptx.eq.7)ispecs(i)=jspe07(i)
+ if(ioptx.eq.9)ispecs(i)=jspe09(i)
if(ioptx.eq.11)ispecs(i)=jspe11(i)
if(ioptx.eq.13)ispecs(i)=jspe13(i)
if(ioptx.eq.15)ispecs(i)=jspe15(i)
c-------------------------------------------------------------
subroutine hnbspf(ku,kd,ks,kc,kb,kt,j,n,spelog)
-c-------------------------------------------------------------
+c-------------------------------------------------------------
c returns spelog = log of factor for consid. different species
c spelog is double precision
c option ioflac determines the method:
c ioflac=2: flavour conservation implemented straightforward
c (only for nspecs=3,7)
c ioflac=3: flavour conservation via generating fctn
-c further input:
+c further input:
c ku,...,kt (integer) : flavour
-c j (integer) : excluded species
-c n (integer) : multiplicity
-c-------------------------------------------------------------
+c j (integer) : excluded species
+c n (integer) : multiplicity
+c-------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cflac/ifok(nflav,mspecs),ifoa(nflav)
integer m(7),l(7),ifot(nflav)
n3=n-n1-n2
do 5 nf=1,nflav
if(ifoa(nf).eq.0.and.ifot(nf).eq.0)goto5
- if(n1*ifok(nf,1)+n2*ifok(nf,2)+n3*ifok(nf,3).ne.ifot(nf))goto2
+ if(n1*ifok(nf,1)+n2*ifok(nf,2)+n3*ifok(nf,3).ne.ifot(nf))goto2
5 continue
spe=spe+utgam2(1.d0+n)
&/utgam2(1.d0+n1)/utgam2(1.d0+n2)/utgam2(1.d0+n3)
do 15 nf=1,nflav
if(ifoa(nf).eq.0.and.ifot(nf).eq.0)goto15
if(n1*ifok(nf,1)+n2*ifok(nf,2)+n3*ifok(nf,3)+n4*ifok(nf,4)
- *+n5*ifok(nf,5)+n6*ifok(nf,6)+n7*ifok(nf,7).ne.ifot(nf))goto12
+ *+n5*ifok(nf,5)+n6*ifok(nf,6)+n7*ifok(nf,7).ne.ifot(nf))goto12
15 continue
spe=spe+1d0/faci(n)*faci(n1)*faci(n2)*faci(n3)*faci(n4)
&*faci(n5)*faci(n6)*faci(n7)
do 16 nf=1,nflav
if(ifoa(nf).eq.0.and.ifot(nf).eq.0)goto16
if(n1*ifok(nf,1)+n2*ifok(nf,2)+n3*ifok(nf,3)+n4*ifok(nf,4)
- *+n5*ifok(nf,5)+n6*ifok(nf,6)+n7*ifok(nf,7).ne.ifot(nf))goto13
+ *+n5*ifok(nf,5)+n6*ifok(nf,6)+n7*ifok(nf,7).ne.ifot(nf))goto13
16 continue
spe=spe+1d0/faci(n)*faci(n1)*faci(n2)*faci(n3)*faci(n4)
&*faci(n5)*faci(n6)*faci(n7)
c-------------------------------------------------------------
subroutine hnbspg(ku,kd,ks,kc,kb,kt,j,n,spelog)
-c-------------------------------------------------------------
+c-------------------------------------------------------------
include 'epos.inc'
double precision spelog,spalog
if(ioflac.ne.0)return
call hnbspf(ku,kd,ks,kc,kb,kt,j,n,spalog)
ioflac=3
call hnbspf(ku,kd,ks,kc,kb,kt,j,n,spelog)
- ioflac=0
+ ioflac=0
write(ifch,*)'ioflac=2/3:',spalog,spelog
return
end
*, 8*2.
*,10*4.
*, 8*2.
- *,10*4.
+ *,10*4.
*,1*3
*,1*3/
do i=1,nspec
fac=fac*(1+facts)
elseif(abs(id).lt.1000)then
fac=fac*(1-facts)
- endif
- endif
+ endif
+ endif
spideg=spideg*fac
goto1
endif
*,' entry sr hnbtst ',('-',i=1,30)
if(ish.ge.7)write(ifch,*)'configuration:'
if(ish.ge.7)write(ifch,*)(ident(i),i=1,np)
- if(ish.ge.7)write(ifch,*)'n_l:',nlattc,' n_0:',nlattc-np
+ if(ish.ge.7)write(ifch,*)'n_l:',nlattc,' n_0:',nlattc-np
c log of prod m_i*volu/4/pi**3/hquer**3 -> f5log
f5log=0
f4log=0
if(ish.ge.7)write(ifch,*)'log(f4):',f4log
-c log of 1/prod n_alpha! -> f3log
+c log of 1/prod n_alpha! -> f3log
dbllog=0
n1=1
nx=1
c log of f3 * f4 * f5
f35log=f5log+f4log+f3log
if(ish.ge.7)write(ifch,*)'log(f3*f4*f5):',f35log
-
+
c log of phase space integral --> psilog
if(iocova.eq.1)then
psilog=alog(2.*np*np*(np-1)/tecm**4/pi)
psulog=psilog
wtulog=w35log
-
+
if(ish.ge.7)write(ifch,*)('-',i=1,30)
*,' exit sr hnbtst ',('-',i=1,10)
ish=ish0
cc x --> x*10.**e with x.lt.10.**10.
cc----------------------------------------------------------------------
c if(x.eq.0.)then
-c e=0.
+c e=0.
c else
c e=int(alog10(abs(x)))/10*10
c x=x/10.**e
c----------------------------------------------------------------------
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nhise=100)
common/chise/hise(mspecs,nhise)
c----------------------------------------------------------------------
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nhismu=500)
common/chismu/hismu(mspecs,0:nhismu),hismus(nhismu)
c-----------------------------------------------------------------------
subroutine xhgcam(amt,iii)
c-----------------------------------------------------------------------
-c creates unnormalized histogram for total mass of grand
+c creates unnormalized histogram for total mass of grand
c canonically generated sample
c xpar1: nr. of bins
c xpar2: m_1 (lower boundary)
character cen*6,cvol*6
if(iii.eq.0)then
-
+
am(nrclu)=amt
return
-
+
elseif(iii.lt.0)then
nbin=nint(xpar3)
parameter(mxclu=10000)
common/cchi/chi2(mxclu)
character cnu*2,cinco*1,cen*6,cvol*6
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
if(chi.ge.0.0)then
c xpar5 specifies statistics to be used ,(0) same as iostat
c (1) boltzmann
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
c-----------------------------------------------------------------------
include 'epos.inc'
common/citer/iter,itermx
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cbol/rmsbol(mspecs),ptlbol(mspecs),chebol(mspecs),tembol
endif
daty(j)=daty(j)+dnde*gspecs(i)*volu/hquer**3/8./pi**3
10 continue
-
+
else
-
+
dnde=0.0
if(datx(j).ge.aspecs(id))then
x=100.
c-----------------------------------------------------------------------
subroutine xhgcfl(u,d,s,iii)
c-----------------------------------------------------------------------
-c creates unnormalized histogram for net flavor content of grand
+c creates unnormalized histogram for net flavor content of grand
c canonically generated sample
c xpar1: specifies width of plot, netflavor centered
c-----------------------------------------------------------------------
character cfl*3,cen*6,cvol*6
if(iii.eq.0)then
-
+
ku(nrclu)=u
kd(nrclu)=d
ks(nrclu)=s
return
-
+
elseif(iii.lt.0)then
kwid=nint(xpar1)
c xpar2 and xpar3 specify xrange of plot
c xpar4 specifies line type : dashed (0), dotted (1), full (2)
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
c-----------------------------------------------------------------------
include 'epos.inc'
common/citer/iter,itermx
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
character cen*6,cvo*6,cit*5,ctem*5
c xpar7 specifies statistics : same as iostat (0)
c boltzmann (1)
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
-c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
include 'epos.inc'
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cbol/rmsbol(mspecs),ptlbol(mspecs),chebol(mspecs),tembol
common/cgctot/rmstot,ptltot
character cyield*8,cen*6,cvo*6,cinco*1
-
+
idpa=nint(xpar1)
ixra=nint(xpar2)
ist=nint(xpar7)
if(ist.eq.0.and.iostat.eq.1)ist=1
-
+
id=0
jx=100
ymin=1./nevent/10.
do i=1,nspecs
if(ispecs(i).eq.idpa)id=i
enddo
-
+
if(ixra.eq.1)then
x1=anint(xpar3)
x2=anint(xpar4)
if(id.eq.0)then
x1=anint(ptltot-iwid*rmstot)
x2=anint(ptltot+iwid*rmstot)
- else
+ else
x1=anint(ptlngc(id)-iwid*rmsngc(id))
x2=anint(ptlngc(id)+iwid*rmsngc(id))
endif
c ------------------
x=100.
if(rmstot.ge.1.e-10)x=(datx(j)-ptltot)**2/rmstot**2/2.
-
+
if(x.ge.60)then
pn=0.0
else
daty(j)=pn
else
-
+
c one species (specified by id)
-c ------------------------------
+c ------------------------------
x=100.
if(rmsngc(id).ge.1.e-10.and.ist.eq.0)
*x=(datx(j)-ptlngc(id))**2/rmsngc(id)**2/2.
do j=1,jx
write(ifhi,'(2e12.4)')datx(j),daty(j)
- enddo
+ enddo
write(ifhi,'(a)') ' endarray'
write(ifhi,'(a)') 'closehisto'
- return
+ return
end
c xpar7 specifies statistics : same as iostat (0)
c boltzmann (1)
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
-c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
include 'epos.inc'
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cbol/rmsbol(mspecs),ptlbol(mspecs),chebol(mspecs),tembol
common/cgctot/rmstot,ptltot
character cyield*8,cen*6,cvo*6,cinco*1
-
+
idpa=nint(xpar1)
ixra=nint(xpar2)
ist=nint(xpar7)
if(ist.eq.0.and.iostat.eq.1)ist=1
-
+
id=0
ymin=1./nevent/10.
if(nevent.le.10)ymin=ymin/10.
do i=1,nspecs
if(ispecs(i).eq.idpa)id=i
enddo
-
+
if(ixra.eq.1)then
n1=nint(xpar3)
n2=nint(xpar4)
if(id.eq.0)then
n1=nint(ptltot-iwid*rmstot)
n2=nint(ptltot+iwid*rmstot)
- else
+ else
n1=nint(ptlngc(id)-iwid*rmsngc(id))
n2=nint(ptlngc(id)+iwid*rmsngc(id))
endif
daty(j)=1./jf*ptltot**(j-1)*exp(-ptltot)
else
-
+
c one species (specified by id)
-c ------------------------------
+c ------------------------------
if(ist.eq.0)pn=1./jf*ptlngc(id)**(j-1)*exp(-ptlngc(id))
if(ist.eq.1)pn=1./jf*ptlbol(id)**(j-1)*exp(-ptlbol(id))
do j=1,jx
write(ifhi,'(2e12.4)')datx(j),daty(j)
- enddo
+ enddo
write(ifhi,'(a)') ' endarray'
write(ifhi,'(a)') 'closehisto'
- return
+ return
end
c-----------------------------------------------------------------------
c xpar3 and xpar4 specify xrange of plot
c xpar5 specifies line type : dashed (0), dotted (1), full (2)
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
c-----------------------------------------------------------------------
include 'epos.inc'
common/citer/iter,itermx
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
character crap*5,cen*6,cvo*6,cit*5
c xpar2 and xpar3 specify xrange of plot
c xpar4 specifies line type : dashed (0), dotted (1), full (2)
c output:
-c histo-file
+c histo-file
c newpage, zone and plot commands not included !!!
c-----------------------------------------------------------------------
include 'epos.inc'
parameter (nbin=200)
real datx(nbin),daty(nbin)
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
common/cgctot/rmstot,ptltot
daty(j)=daty(j)+dndy
-10 continue
+10 continue
else
endif
enddo
-
+
write(cen,'(f6.1)')tecm
write(cvo,'(f6.1)')volu
if(id.eq.0)then
c xpar3: 1: de/d3p 2: ede/d3e
c-----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nhise=100)
common/chise/hise(mspecs,nhise)
write(ifhi,'(a,a)') 'text 0 0 "title id='//chid
* ,' N='//cyield//' T='//ctem//'"'
write(ifhi,'(a)') 'text 0 0 "xaxis energy (GeV)"'
- write(ifhi,'(a)') 'text 0 0 "yaxis '//ch//' dn/d3p (GeV-3)"'
+ write(ifhi,'(a)') 'text 0 0 "yaxis '//ch//' dn/d3p (GeV-3)"'
write(ifhi,'(a)') 'array 2'
do i=1,nhise
if(mode.eq.1)write(ifhi,'(2e12.4)')datx(i),daty(i)
c xpar2: 1:actual multiplicity 2:average multiplicity 3:grand canonical
c-----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (literm=500)
common/cmet/kspecs(mspecs),liter,lspecs(literm,mspecs)
c xpar3,4: xrange
c-----------------------------------------------------------------------
include 'epos.inc'
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nhismu=500)
common/chismu/hismu(mspecs,0:nhismu),hismus(nhismu)
* ,cvolu//'"'
write(ifhi,'(a)') 'text 0 0 "xaxis multiplicity n "'
write(ifhi,'(a)') 'text 0 0 "yaxis dN/dn"'
- write(ifhi,'(a)') 'text 0.30 0.25 "N?MC!='//cyield//'"'
+ write(ifhi,'(a)') 'text 0.30 0.25 "N?MC!='//cyield//'"'
write(ifhi,'(a)') 'array 2'
do i=1,jx
write(ifhi,'(2e12.4)') datx(i),daty(i)
write(ifhi,'(a)') 'text 0 0 "title id='//chid//'"'
write(ifhi,'(a)') 'text 0 0 "xaxis multiplicity n "'
write(ifhi,'(a)') 'text 0 0 "yaxis dN/dn"'
- write(ifhi,'(a)') 'text 0.30 0.25 "N?MC!='//cyield//'"'
+ write(ifhi,'(a)') 'text 0.30 0.25 "N?MC!='//cyield//'"'
write(ifhi,'(a)') 'array 2'
do i=1,jx
write(ifhi,'(2e12.4)') datx(i),daty(i)
c-----------------------------------------------------------------------
subroutine xhnbmz
c-----------------------------------------------------------------------
-c produces histogram of multiplicity distribution from droplet decay
+c produces histogram of multiplicity distribution from droplet decay
c or average multiplicity versus iterations
c for massless hadrons
c complete histogram: openhisto ... closehisto
c xpar3: upper limit multiplicity
c xpar4: lower limit total multiplicity (also necc for xpar1.ne.0)
c xpar5: upper limit " " (also necc for xpar1.ne.0)
-c xpar6: sets htyp: 1->lfu, 2->ldo, 3->lda, 4->ldd
+c xpar6: sets htyp: 1->lfu, 2->ldo, 3->lda, 4->ldd
c xpar7: 0: multiplicity distribution
c >0: av multiplicity vs iterations (itermx=xpar7)
c-----------------------------------------------------------------------
parameter(maxp=500)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
common/ctst/psulog,wtulog
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
parameter (nhismu=500)
common/cflac/ifok(nflav,mspecs),ifoa(nflav)
call hnbtst(0)
wtzlog=wtulog
if(ioflac.eq.0)call hnbspg(keu,ked,kes,kec,keb,ket,0,np,spelog)
- if(ioflac.ne.0)call hnbspf(keu,ked,kes,kec,keb,ket,0,np,spelog)
+ if(ioflac.ne.0)call hnbspf(keu,ked,kes,kec,keb,ket,0,np,spelog)
wtulog=wtulog+spelog
else
wtzlog=-1000
datyu(l)=exp(datyu(l))
else
datyu(l)=exp(-50.)
- endif
+ endif
yield=yield+i*datyu(l)
su=su+datyu(l)
enddo
if(idcode.eq.0.and.itmax.eq.0)then
write(ifhi,'(a,2e11.3)')'openhisto xrange',x1,x2
write(ifhi,'(a)') 'htyp '//htyp//' xmod lin ymod log'
- write(ifhi,'(a)') 'text 0.30 0.15 "N?ana!='//cyieur//'"'
+ write(ifhi,'(a)') 'text 0.30 0.15 "N?ana!='//cyieur//'"'
write(ifhi,'(a)') 'array 2'
do i=1,jx
write(ifhi,'(2e12.4)') datx(i),datyu(i)
if(itmax.eq.0)then
write(ifhi,'(a,2e11.3)')'openhisto xrange',x1,x2
write(ifhi,'(a)') 'htyp '//htyp//' xmod lin ymod log'
- write(ifhi,'(a)') 'text 0.30 0.15 "N?ana!='//cyieur//'"'
+ write(ifhi,'(a)') 'text 0.30 0.15 "N?ana!='//cyieur//'"'
write(ifhi,'(a)') 'array 2'
do i=1,jx
write(ifhi,'(2e12.4)') datx(i),datyu(i)
c regarding exponential autocorrelation time and acceptance rate
c
c input:
-c requires complete run with application hadron (iappl=1)
+c requires complete run with application hadron (iappl=1)
c or application metropolis (iappl=4)
c ioceau=1 necessary
c
c for iii=0 (only valid for iappl=4):
c data(nrevt): nrevt (event number) /cdat/
c datb(nrevt): taui (calculated corr time) /cdat/
-c datc(nrevt): accrat (acceptance rate) /cdat/
+c datc(nrevt): accrat (acceptance rate) /cdat/
c datd(nrevt): taue (parametrized corr time) /cdat/
c for iii>0 (only valid for iappl=1):
c nrclu=nrclu+1 /cnrclu/
c data(nrclu): nrclu (droplet number) /cdat/
c datb(nrclu): taui-taue (calc - param corr time) /cdat/
-c datc(nrclu): accrat (acceptance rate) /cdat/
+c datc(nrclu): accrat (acceptance rate) /cdat/
c datd(nrclu): avnp (average particle number) /cdat/
c for iii<0:
c writes complete histogram (openhisto ... closehisto) to histofile
common/citer/iter,itermx
common /clatt/nlattc,npmax
common/cgctot/rmstot,ptltot
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
parameter (nbin=500)
datc(nrclu)=accrat
datd(nrclu)=avnp
endif
-
-c -----------------------------------
+
+c -----------------------------------
elseif(iii.lt.0.and.iappl.eq.4)then
c -----------------------------------
endif
write(cvol,'(f7.3)')volu
write(clatt,'(i3)')nlattc
- write(cit,'(i5)')itermx
+ write(cit,'(i5)')itermx
if(ioobsv.eq.0)then
write(cobs,'(a)')'all'
else
write(cnc,'(i5)')iternc
if(iozevt.eq.0)write(czer,'(i5)')iozero
if(iozevt.gt.0)write(cdz,'(i5)')iozinc
-
+
x1=1
x2=nevent
endif
-c -----------------------------------
+c -----------------------------------
elseif(iii.lt.0.and.iappl.eq.1)then
c -----------------------------------
-
+
if(ioobsv.eq.0)then
write(cobs,'(a)')'all'
else
c for iii=0 (only valid for iappl=4):
c data(nrevt): nrevt (event number) /cdat/
c datb(nrevt): tau (calculated int corr time) /cdat/
-c datc(nrevt): stau (variance tau) /cdat/
+c datc(nrevt): stau (variance tau) /cdat/
c datd(nrevt): avnp (multiplicity) /cdat/
c date(nrevt): sobs (variance multiplicity) /cdat/
c datf(nrevt): (gc multiplicity) /cdat/
c for iii=0 and iosngl>0:
c writes complete set of histograms (newpage zone 1 3 1
c openhisto ... closehisto plot0 ... openhisto ... closehisto plot 0)
-c concerning acceptance rate, rejection rate, correlation function
+c concerning acceptance rate, rejection rate, correlation function
c for specific event, specified by value of iosngl (=nrevt+1)
c for iii<0:
c writes complete histogram (openhisto ... closehisto) to histofile
-c xpar1=1: (data,datb,datc)
-c xpar1=2: (data,datd,date,datf)
+c xpar1=1: (data,datb,datc)
+c xpar1=2: (data,datd,date,datf)
c------------------------------------------------------------------------
include 'epos.inc'
parameter(maxit=50000)
common/confg/np,amass(maxp),ident(maxp),pcm(5,maxp),wtxlog,wtlog
common /clatt/nlattc,npmax
common/cgctot/rmstot,ptltot
- parameter (mspecs=56)
+ parameter (mspecs=56)
common/cspecs/nspecs,ispecs(mspecs),aspecs(mspecs),gspecs(mspecs)
common/cgchg/rmsngc(mspecs),ptlngc(mspecs),chemgc(mspecs),tem
parameter (nbin=500)
c ----------------
if(iii.eq.0)then
c ----------------
-
+
c mean
c ----
xnptot=nptot
c -----
endif
-c -----
+c -----
- return
+ return
end
c-----------------------------------------------------------------------
subroutine utresc(iret)
c-----------------------------------------------------------------------
-c if irescl=1 rescaling is done, otherwise the purpose of going through
+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'
call utpri('utresc',ish,ishini,4)
errlim=0.001 !max(0.001,1./engy)
-
+
iret=0
scal=1.
nptlpt=iabs(maproj)+iabs(matarg)
call ranfgt(seedp) !not to change the seed ...
if(nptl.le.nptlpt) goto 9999
-
+
esoll=0.d0
p1(1)=0.d0
p1(2)=0.d0
if(ish.ge.2)write (ifch,*) 'Problem in utresc (0): redo event'
write (ifmt,*) 'Problem in utresc (0): redo event'
c call utstop('utresc&')
- goto 9999
+ goto 9999
endif
esoll=p1(5)
if(ish.ge.4) write (ifch,'(a,5g13.6)') 'boost-vector: ',p1
if(ish.ge.6)then
call alistf('list after boost&')
- endif
+ endif
if(ish.ge.5)write(ifch,'(a)')'--------rescale momenta----------'
pini(j)=pptl(3,j)
else
pptl3new=0.
- if( force .or.(
+ 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)
+ & 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(ish.ge.6)write(ifch,*)
$ 'ipass,scal,diff/pnullx,e,esoll,pz,pmax,ndif,f:'
$ ,ipass,scal,diff/pnullx,sum,esoll,sum3,pnullx,ndif,force
- if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).lt.errlim)
+ if(abs(scal-1.).le.errlim.and.abs(diff/pnullx).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/pnullx).lt.5e-2)then
+ elseif(abs(scal-1.).le.1e-2.and.abs(diff/pnullx).lt.5e-2)then
goto 300
elseif(force)then
if(ish.ge.2)
iret=1
goto 9999
endif
-
+
c trafo
c -----
300 continue
endif
endif
enddo
-
+
if(ish.ge.4)call alist('list after rescaling&',1,nptl)
-
+
9999 continue
if(ish.ge.2)then
scalmean=scalmean+scal
write(ifch,*)' average rescaling factor: ',scalmean
& /float(nrevt+1)
- endif
- call ranfst(seedp) ! ... after this subroutine
+ endif
+ call ranfst(seedp) ! ... after this subroutine
call utprix('utresc',ish,ishini,4)
end
nptlpt=iabs(maproj)+iabs(matarg)
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
+ do j=nptlpt+1,nptl
if(mod(istptl(j),10).eq.0.and.pptl(4,j).gt.0.d0)then
amass=pptl(5,j)
call idmass(idptl(j),amass)
$ .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
+ if(ish.ge.5)then
write(ifch,'(a,$)')'ghost:'
call alistc("&",j,j)
endif
ityptl(j)=100+ityptl(j)/10
- endif
+ endif
elseif(mod(istptl(j),10).eq.0)then
if(ish.ge.1)then
write(ifmt,*)'Lost particle (E=0)'
enddo
if(ish.ge.5)write(ifch,'(a)')'---------treat ghosts---------'
-
+
c treat ghosts
c ------------
ifirst=1
psum=0
esum=0.
ntry=ntry+1
- do j=nptlpt+1,nptl
+ do j=nptlpt+1,nptl
if(mod(istptl(j),10).eq.0)then
if(ityptl(j).ge.101.and.ityptl(j).le.105)then
nfif=nfif+1
- if(ifirst.eq.1)then
+ 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(irescl.ge.1) then
if(ifirst.gt.1)then
if(pptl(4,j).gt.0.)then
Einv=1./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
endif
enddo
if ( ish.ge.5 ) write (ifch,*) 'tot',nfif,efif,pfif,esum,psum
c Check Ghost list
- if(ish.ge.5)then
- do j=nptlpt+1,nptl
+ 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:'
enddo
endif
-
+
9999 continue
- call ranfst(seedp) ! ... after this subroutine
+ 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 if irescl=1 and ispherio=1 rescaling is done for particle used by
c spherio as initial condition.
c-----------------------------------------------------------------------
include 'epos.inc'
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
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
+ $ .or.istptl(i).eq.20.or.istptl(i).eq.21)then
do j=1,2
p1(j)=p1(j)+dble(pptl(j,i))
enddo
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
+ $ .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
sum=0.
sum3=0.
ndif=0
- do j=1,nptl
+ 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
+ $ .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)
call utmsgf
endif
-
+
c trafo
c -----
300 continue
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
+ $ .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
endif
enddo
-
+
9999 call utprix('utrsph',ish,ishini,4)
end
c dddlog=-1d50
c if(xxx.gt.0d0)dddlog=dlog(xxx)
c end
-c
+c
ccc-----------------------------------------------------------------------
c subroutine randfl(jc,iqa0,iflav,ic,isame)
cc-----------------------------------------------------------------------
c *,' exit sr randfl ',('-',i=1,10)
c return
c end
-c
-c
+c
+c
cc-----------------------------------------------------------------------
c subroutine ranhvy(x,eps)
cc-----------------------------------------------------------------------
c x=z
c return
c end
-c
+c
c-----------------------------------------------------------------------
function ransig()
c-----------------------------------------------------------------------
if(rangen().gt.0.5)ransig=-1
return
end
-
+
cc-----------------------------------------------------------------------
c function ranxq(n,x,q,xmin)
cc-----------------------------------------------------------------------
c qran=q(imin)+rangen()*(q(n)-q(imin))
c ranxq=utinvt(n,x,q,qran)
c4 continue
-c
+c
c if(ranxq.lt.xmin)then
c call utmsg('ranxq ')
c write(ifch,*)'***** ranxq=',ranxq,' < xmin=',xmin
c call utmsgf
c ranxq=xmin
c endif
-c
+c
c return
c end
-c
+c
cc ***** end r-routines
-cc ***** beg s-routines
+cc ***** beg s-routines
c
cc-----------------------------------------------------------------------
c function sbet(z,w)
c sbet=utgam1(z)*utgam1(w)/utgam1(z+w)
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c function smass(a,y,z)
cc-----------------------------------------------------------------------
c *+dy/a/2.*(y-ymin)**2
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine smassi(theta)
cc-----------------------------------------------------------------------
cc-----------------------------------------------------------------------
c common/cmass/thet,epsi,as,ac,dy,dz,ym,cz,zm,sigma,rzero
c thet=theta
-c
+c
c astr=.150
c pi=3.14159
c alp=1./137.
-c
+c
c co=cos(theta)
c si=sin(theta)
c bet=(1+co**3)/2.
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
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 *(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
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine smassp
cc-----------------------------------------------------------------------
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
c stop
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine smasst(kux,kdx,ksx,kcx,a,y,z)
cc-----------------------------------------------------------------------
c z=nz/3
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine smassu(ax,yx,zx,ku,kd,ks,kc)
cc-----------------------------------------------------------------------
c kc=0
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c function spoc(a,b,c,d,x)
cc-----------------------------------------------------------------------
utacos=acos(argum)
return
end
-
+
c----------------------------------------------------------------------
function utamnu(keux,kedx,kesx,kecx,kebx,ketx,modus)
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'
1 format(' flavours:',6i5 )
100 format(' flavours+mass:',6i5,f8.2 )
c write(ifch,1)keux,kedx,kesx,kecx,kebx,ketx c write(ifmt,*)'wrong mass in gethad ',damss
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
ket=-ketx
endif
-c write(ifch,*)keu,ked,kes,kec,keb,ket
+c write(ifch,*)keu,ked,kes,kec,keb,ket
-c removing top mesons to remove t quarks or antiquarks
+c removing top mesons to remove t quarks or antiquarks
if(ket.ne.0)then
12 continue
ii=sign(1,ket)
if(ket.ne.0)goto12
endif
-c removing bottom mesons to remove b quarks or antiquarks
+c removing bottom mesons to remove b quarks or antiquarks
if(keb.ne.0)then
11 continue
ii=sign(1,keb)
amnull=amnull+1.87 ! (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
+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)
kes=kes+1
goto5
endif
-
-c removing mesons to remove d antiquarks
+
+c removing mesons to remove d antiquarks
6 continue
if(ked.lt.0)then
if(keu.ge.kes)then
ked=ked+1
goto6
endif
-
-c removing mesons to remove u antiquarks
+
+c removing mesons to remove u antiquarks
7 continue
if(keu.lt.0)then
if(ked.ge.kes)then
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
-
+
+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
+
+c removing strange baryons
i=4
2 i=i-1
3 continue
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
+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)
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
-
+
+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
-
+
+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)
-
+
1000 utamnu=amnull
return
end
-
+
c-----------------------------------------------------------------------
function utamnx(jcp,jcm)
c-----------------------------------------------------------------------
end
-
+
cc-----------------------------------------------------------------------
c function utamny(jcp,jcm)
cc-----------------------------------------------------------------------
utamnz=utamnu(keu,ked,kes,kec,keb,ket,modus)
return
end
-
+
c-----------------------------------------------------------------------
subroutine utar(i1,i2,i3,x0,x1,x2,x3,xx)
c-----------------------------------------------------------------------
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 *,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 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 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 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 call utmsgf
c endif
c a=dsqrt( (pj1-pi1)**2 + (pj2-pi2)**2 + (pj3-pi3)**2 )
c if(a.eq.0.d0)return
c a3=sngl((pi3-pj3)/a)
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine utchm(arp,arm,ii)
cc-----------------------------------------------------------------------
c endif
c return
c end
-c
+c
c-----------------------------------------------------------------------
subroutine utclea(nptlii,nptl0)
c-----------------------------------------------------------------------
c starting from nptlii
-c overwrites istptl=99 particles in /cptl/, reduces so nptl
+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)
newptl(i)=i
c oldptl(i)=i
goto 1
-
+
2 i=i-1
j=i
3 i=i+1
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)
+ 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
+ nptl0=0
goto 1000
endif
c11 continue
c do 12 k=1,nptl
c12 iorptl(k)=ii(k)
-c
+c
c do 13 k=1,nptl
c jo=jorptl(k)
c if(jo.le.0)ii(k)=jo
c13 continue
c do 14 k=1,nptl
c14 jorptl(k)=ii(k)
-c
+c
c do 15 k=1,nptl
c if1=ifrptl(1,k)
c if(if1.le.0)ii(k)=if1
c15 continue
c do 16 k=1,nptl
c16 ifrptl(1,k)=ii(k)
-c
+c
c do 17 k=1,nptl
c if2=ifrptl(2,k)
c if(if2.le.0)ii(k)=if2
c17 continue
c do 18 k=1,nptl
c18 ifrptl(2,k)=ii(k)
-c
+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
*,istptl(n),ityptl(n)
35 continue
endif
-
+
if(ish.ge.2)write(ifch,*)'exiting subr utclea:',nptl
& ,minfra,maxfra
c---------------------------------------------------------------------
subroutine utfit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
c---------------------------------------------------------------------
-c linear fit to data
+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---------------------------------------------------------------------
+c---------------------------------------------------------------------
INTEGER mwt,ndata
REAL a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
-CU USES utgmq
+CU USES utgmq
INTEGER i
REAL sigdat,ss,st2,sx,sxoss,sy,t,wt,utgmq
sx=0.
endif
return
END
-
+
c-----------------------------------------------------------------------
function utgam1(x)
c-----------------------------------------------------------------------
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)
+ 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
1 gamser=sum*exp(-x+a*log(x)-gln)
return
END
-
+
c-------------------------------------------------------------------------
subroutine uticpl(ic,ifla,iqaq,iret)
-c-------------------------------------------------------------------------
+c-------------------------------------------------------------------------
c adds a quark (iqaq=1) or antiquark (iqaq=2) of flavour ifla
c to 2-id ic
-c-------------------------------------------------------------------------
+c-------------------------------------------------------------------------
include 'epos.inc'
integer jc(nflav,2),ic(2)
iret=0
cc-----------------------------------------------------------------------
c subroutine utindx(n,xar,x,i)
cc-----------------------------------------------------------------------
-cc input: dimension n
-cc array xar(n) with xar(i) > xar(i-1)
+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 i=lu
c return
c end
-c
+c
c-----------------------------------------------------------------------
function utinvt(n,x,q,y)
c-----------------------------------------------------------------------
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 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 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-----------------------------------------------------------------------
return
end
-
+
c-----------------------------------------------------------------------
subroutine utlobo(isig,p1,p2,p3,p4,p5,x1,x2,x3,x4)
c-----------------------------------------------------------------------
x4=z(4)
return
end
-
+
c-----------------------------------------------------------------------
subroutine utloc(ar,n,a,l)
c-----------------------------------------------------------------------
l=n
return
end
-
+
cc-----------------------------------------------------------------------
c subroutine utlow(cone)
cc-----------------------------------------------------------------------
c if(cone.eq.'Z')cone='z'
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine utlow3(cthree)
cc-----------------------------------------------------------------------
c1 call utlow(cthree(i:i))
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c subroutine utlow6(csix)
cc-----------------------------------------------------------------------
c1 call utlow(csix(i:i))
c return
c end
-c
+c
cc-----------------------------------------------------------------------
c function utmom(k,n,x,q)
cc-----------------------------------------------------------------------
c utmom=utmom/q(n)
c return
c end
-c
+c
c-----------------------------------------------------------------------
function utpcm(a,b,c)
c-----------------------------------------------------------------------
do nr=1,nrpri
if(subpri(nr)(1:6).eq.text)then
ish=ishpri(nr)
- endif
+ endif
enddo
endif
if(ish.ge.ishx)then
endif
return
end
-
+
c-----------------------------------------------------------------------
subroutine utprix(text,idmmy,ishini,ishx)
c-----------------------------------------------------------------------
ish=ishini
return
end
-
+
c-----------------------------------------------------------------------
subroutine utprj(text,idmmy,ishini,ishx)
c-----------------------------------------------------------------------
do nr=1,nrpri
if(subpri(nr)(1:idx).eq.text(1:idx))then
ish=ishpri(nr)
- endif
+ endif
enddo
endif
if(ish.ge.ishx)then
endif
return
end
-
+
c-----------------------------------------------------------------------
subroutine utprjx(text,idmmy,ishini,ishx)
c-----------------------------------------------------------------------
1 utquad=utquad+(f(i)+f(i+1))/2*(x(i+1)-x(i))
return
end
-
+
c-----------------------------------------------------------------------
subroutine utquaf(fu,n,x,q,wo,x0,x1,x2,x3)
c-----------------------------------------------------------------------
rinptl(i) =rinptl(j)
return
end
-
+
cc-----------------------------------------------------------------------
c subroutine utresm(icp1,icp2,icm1,icm2,amp,idpr,iadj,ireten)
cc-----------------------------------------------------------------------
cc funcd: subr returning fctn value and first derivative
cc x1,x2: x-interval
cc xacc: accuracy
-cc output:
+cc output:
cc utroot: root
cc--------------------------------------------------------------------
c include 'epos.inc'
z=zs
return
end
-
+
c-----------------------------------------------------------------------
subroutine utrot4(isig,ax,ay,az,x,y,z)
subroutine utrota(isig,ax,ay,az,x,y,z)
c-----------------------------------------------------------------------
c performs a rotation
-c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
if(az.ge.0.)then
rx=ax
ry=ay
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 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-----------------------------------------------------------------------
if(ish.eq.1)return
write(ifch,'(1x,74a1)')('*',j=1,72)
end
-
+
c-----------------------------------------------------------------
subroutine uttrap(func,a,b,s)
c-----------------------------------------------------------------
if (ds.lt.epsr*abs(olds)) return
olds=s
11 continue
-
+
c-c nepsr=nepsr+1
if(ish.ge.9)then
call utmsg('uttrap')
c finds the first word of the character string line(j+1:160).
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 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 input:
c line: character string (*160)
c i: integer between 1 and 160
c iqu: for iqu=1 a "?" is written to output before reading a line,
c otherwise (iqu/=1) nothing is typed
-c output:
+c output:
c i,j: left and right end of word (word=line(i:j))
c-------------------------------------------------------------------
include 'epos.inc'
character w1define*100,w2define*100
common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
& ,w1define(mxdefine),w2define(mxdefine)
-
+
if(j.ge.0)then
i=j
goto1
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
+ 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
+ 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
endif
enddo
enddo
return
end
-c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
subroutine getairmol(iz,ia)
-c-----------------------------------------------------------------------
+c-----------------------------------------------------------------------
include 'epos.inc'
- i=0
- r=rangen()
+ 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))
+ 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----------------------------------------------------------------------
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-----------------------------------------------------------------------
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff