IAP GITLAB

Commit 24d782a0 authored by Tanguy Pierog's avatar Tanguy Pierog

add urqmd 1.3 as in CONEX/CORSIKA. No test done and interface not updated (only EPOS part)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@7467 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 3c2dc888
......@@ -105,7 +105,7 @@ LIST(APPEND CRMC_SOURCES src/OutputPolicyHepMC.cc)
add_definitions(-DWITH_HEPMC)
MESSAGE("Build HEPMC Output Interface")
FIND_PACKAGE (Fastjet REQUIRED COMPONENTS fastjet)
FIND_PACKAGE (Fastjet COMPONENTS fastjet)
if (FASTJET_FOUND)
include_directories("${Fastjet_INCLUDE_DIRS}")
add_definitions(-DWITH_FASTJET)
......
......@@ -56,6 +56,12 @@
#endif
model=9
endif
if(cmodel(1:n).eq.'urqmd' )then
#ifndef __URQMD__
stop'please compile with requested model'
#endif
model=10
endif
if(cmodel(1:n).eq.'dpmjet' )then
#ifndef __DPMJET__
stop'please compile with requested model'
......@@ -126,6 +132,13 @@
stop'please compile with requested model'
#else
call IniFluka
#endif
endif
if(model.eq.10) then
#ifndef __URQMD__
stop'please compile with requested model'
#else
call IniUrQMD
#endif
endif
if(model.eq.12) then
......@@ -207,6 +220,13 @@
stop'please compile with requested model'
#else
call IniEvtFlu
#endif
endif
if(model.eq.10)then
#ifndef __URQMD__
stop'please compile with requested model'
#else
call IniEvtUrq
#endif
endif
if(model.eq.12)then
......@@ -284,6 +304,13 @@
stop'please compile with requested model'
#else
call emsflu(iret)
#endif
endif
if(model.eq.10) then
#ifndef __URQMD__
stop'please compile with requested model'
#else
call emsurq(iret)
#endif
endif
if(model.eq.12) then
......@@ -321,6 +348,7 @@
subroutine crseaaModel0(sigt,sigi,sigc,sige)
include 'epos.inc'
common/geom/rmproj,rmtarg,bmax,bkmx
double precision GTOT,GPROD,GABS,GDD,GQEL,GCOH
double precision e0
......@@ -424,6 +452,15 @@
sige=0.
sigc=0.
sigi=flucrse(ekin,maproj,matarg,idtargin)
#endif
elseif(model.eq.10)then
#ifndef __URQMD__
stop'please compile with requested model'
#else
sigi=0.
sige=0.
sigc=0.
sigt=urqcrse(ekin,idproj,maproj,matarg,idtarg,bmax)
#endif
elseif(model.eq.12)then
#ifndef __DPMJET__
......
Copyright
UrQMD source and documentation are provided freely for the purpose of
checking and reproducing published results of the authors.
The Open Standard Codes and Routines (OSCAR)-Group has established -
for good reasons - guidelines for reproducablity, usage and quality
control of simlulations codes for pA and AA collsions.
UrQMD is a complex model. In order to ensure that it is used correctly
that all results are reproducible and that the proper credits are
given we ask for your agreement to the following copyright and
safeguard mechanisms in the OSCAR spirit.
The UrQMD collaboration favors cooperation and joint projects with
outside researchers. We encourage experimental collaborations to
compare their results to UrQMD. We support you and/or cooperate on any
sensible project related to UrQMD
If you are interested in a project, please contact us.
Projects without the participation of the UrQMD-Collaboration are
accepted, if the project is not a current thesis topic of any
UrQMD-Collaboration member.
We expect that the code authors are informed about any changes and
modifications made to the code. Any changes to the official version
must be documented.
The code or any fragments of it shall not be given away to third
parties. Similarily, events generated with UrQMD shall not be given
to third parties without consent of the code authors.
10.01.2008 UrQMD 1.3 to be used with Corsika.
The stand-alone version of the UrQMD code is available at http://th.physik.uni-frankfurt.de/~urqmd/
This diff is collapsed.
c $Id: addpart.f,v 1.4 2000/01/12 16:02:32 bass Exp $
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
subroutine addpart(index)
c
c Revision : 1.0
c
cinput index : index for slot to create in particle arrays
c
c This subroutine creates an entry for a particle with index {\tt index} in all
c particle arrays.
c \danger{The {\tt nbar} and {\tt nmes} counters must be updated by the
c calling routine !}
c
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
include 'coms.f'
include 'newpart.f'
include 'freezeout.f'
integer ind,i,j,index
ind=index
c now shift vectors upwards
do 10 i=npart,ind,-1
r0(i+1)=r0(i)
rx(i+1)=rx(i)
ry(i+1)=ry(i)
rz(i+1)=rz(i)
p0(i+1)=p0(i)
px(i+1)=px(i)
py(i+1)=py(i)
pz(i+1)=pz(i)
fmass(i+1)=fmass(i)
ityp(i+1)=ityp(i)
iso3(i+1)=iso3(i)
ncoll(i+1)=ncoll(i)
lstcoll(i+1)=lstcoll(i)
charge(i+1)=charge(i)
spin(i+1)=spin(i)
dectime(i+1)=dectime(i)
tform(i+1)=tform(i)
xtotfac(i+1)=xtotfac(i)
origin(i+1)=origin(i)
strid(i+1)=strid(i)
uid(i+1)=uid(i)
frr0(i+1)=frr0(i)
frrx(i+1)=frrx(i)
frry(i+1)=frry(i)
frrz(i+1)=frrz(i)
frp0(i+1)=frp0(i)
frpx(i+1)=frpx(i)
frpy(i+1)=frpy(i)
frpz(i+1)=frpz(i)
ffermpx(i+1)=ffermpx(i)
ffermpy(i+1)=ffermpy(i)
ffermpz(i+1)=ffermpz(i)
r0_t(i+1)=r0_t(i)
rx_t(i+1)=rx_t(i)
ry_t(i+1)=ry_t(i)
rz_t(i+1)=rz_t(i)
do 11 j=1,2
p0td(j,i+1)=p0td(j,i)
pxtd(j,i+1)=pxtd(j,i)
pytd(j,i+1)=pytd(j,i)
pztd(j,i+1)=pztd(j,i)
fmasstd(j,i+1)=fmasstd(j,i)
ityptd(j,i+1)=ityptd(j,i)
iso3td(j,i+1)=iso3td(j,i)
11 continue
10 continue
c increment npart
npart=npart+1
c
if(npart.ge.nmax) then
write(6,*)'*** error in addpart:too much particles>',nmax
write(6,*)' -> increase nmax in coms.f '
call file14out(999)
stop
endif
c initialize new slot
r0(ind)=0.0
rx(ind)=0.0
ry(ind)=0.0
rz(ind)=0.0
p0(ind)=0.0
px(ind)=0.0
py(ind)=0.0
pz(ind)=0.0
fmass(ind)=0.0
ityp(ind)=0
iso3(ind)=0
ncoll(ind)=0
lstcoll(ind)=0
charge(ind)=0
spin(ind)=-1
dectime(ind)=0.d0
tform(ind)=0.0d0
xtotfac(ind)=1.0d0
origin(ind)=0
strid(ind)=0
uid(ind)=0
frr0(ind)=0.d0
frrx(ind)=0.d0
frry(ind)=0.d0
frrz(ind)=0.d0
frp0(ind)=0.d0
frpx(ind)=0.d0
frpy(ind)=0.d0
frpz(ind)=0.d0
ffermpx(ind)=0.d0
ffermpy(ind)=0.d0
ffermpz(ind)=0.d0
cpot
r0_t(ind)=0.d0
rx_t(ind)=0.d0
ry_t(ind)=0.d0
rz_t(ind)=0.d0
ctd
do 12 j=1,2
p0td(j,ind)=0.d0
pxtd(j,ind)=0.d0
pytd(j,ind)=0.d0
pztd(j,ind)=0.d0
fmasstd(j,ind)=0.d0
ityptd(j,ind)=0
iso3td(j,ind)=0
12 continue
c rescan collision table - all particle indices which have been
c shifted upwards must be modified in the collision tables, too.
call scantab(ind,1)
c the lstcoll array must also be shifted due to the new particle slot
do 15 i=1,npart
if(lstcoll(i).le.nmax) then
if(lstcoll(i).eq.ind) then
lstcoll(i)=0
elseif(lstcoll(i).gt.ind) then
lstcoll(i)=lstcoll(i) + 1
endif
endif
15 continue
return
end
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
c $Id: boxinc.f,v 1.3 1999/01/18 09:56:53 ernst Exp $
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Unit : all modules using the box
c Version : 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c max count of different species
integer bptmax
parameter(bptmax=20)
c counter
integer cbox
c Flags
integer boxflag
integer edensflag, mbflag
integer para,solid,mtest
c number of different species
integer mbox
c particle counters
integer bptpart(bptmax),bptityp(bptmax),bptiso3(bptmax)
real*8 bptpmax(bptmax)
real*8 edens
c edge length of a cube in fm
real*8 lbox
c half edge lengt of a cube
real*8 lboxhalbe
c double edge length of a cube
real*8 lboxd
c momenta
real*8 mbp0, mbpx, mbpy, mbpz
common /boxic/ cbox,boxflag, mbox, bptityp, bptiso3, bptpart
common /boxic/ edensflag, para, solid, mbflag, mtest
common /boxrc/ lbox, lboxhalbe, lboxd, bptpmax, edens
common /boxrc/ mbp0, mbpx, mbpy, mbpz
c $Id: boxprg.f,v 1.8 1999/01/18 09:56:53 ernst Exp $
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine bptinit(ibox)
c
c Unit : Initis all the particles setted by the bpt command
c Version : 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'coms.f'
include 'comres.f'
include 'boxinc.f'
include 'options.f'
c Var
c counter, spin
integer i,ibox,fchg,getspin
c randomnumbergenerator, particlemass, deacy times
real*8 ranf,massit,dectim
c momentum, angle distribution
real*8 P,cost,sint,phi
ecm=0.d0
ebeam=0.d0
c main program
c loop over all particles
do 42 i=npart+1,npart+bptpart(ibox)
c configuration space
r0(i)=0.d0
rx(i)=lboxhalbe*(1-2*ranf(0))
ry(i)=lboxhalbe*(1-2*ranf(0))
rz(i)=lboxhalbe*(1-2*ranf(0))
c isospin and ityp
iso3(i)=bptiso3(ibox)
ityp(i)=bptityp(ibox)
c set baryon and meson numbers
if(abs(ityp(i)).le.maxbar) then
nbar=nbar+1
else
nmes=nmes+1
endif
c charge
charge(i)=fchg(iso3(i),ityp(i))
c massarray
fmass(i)=massit(ityp(i))
c Spin
spin(i)=getspin(ityp(i),-1)
c decaytime
dectime(i)=dectim(i,1)
42 continue
c End of loop
if (edensflag.le.0) then
c homogenious momentum distribution, randomly distributed
c max momentum is a parameter
do 45 i=npart+1,npart+bptpart(ibox)
P=bptpmax(ibox)*ranf(0)**(1./3.)
cost = 1.-2.*ranf(0)
sint = sqrt(1.-cost**2)
phi = 2.*Pi*ranf(0)
px(i) = P*sint*cos(phi)
py(i) = P*sint*sin(phi)
pz(i) = P*cost
call setonshell(i)
45 continue
elseif (edensflag.ge.1) then
c energiedensity
c loop over all particles
do 60 i=npart+1,npart+bptpart(ibox)
P=bptpmax(ibox)/bptpart(ibox)*ranf(0)**(1./3.)
cost = 1.-2.*ranf(0)
sint = sqrt(1.-cost**2)
phi = 2.*Pi*ranf(0)
c different initialisations
if (para.eq.0) then
c Boxmode
if (i.eq.1) write(*,*) 'Boxmode'
px(i) = P*sint*cos(phi)
py(i) = P*sint*sin(phi)
pz(i) = P*cost
elseif(para.eq.1) then
c stream over stream
if (i.eq.1) write(*,*) 'streammode'
px(i) = 0.d0
py(i) = 0.d0
pz(i) = bptpmax(ibox)/bptpart(ibox)*(-1.d0)**i
elseif(para.eq.2) then
c slab on slab
if (i.eq.1) write(*,*) 'slabmode'
px(i)=0.d0
py(i)=0.d0
if (rz(i).gt.0) then
pz(i)=(-1.0d0)*bptpmax(ibox)/bptpart(ibox)
else
pz(i)=bptpmax(ibox)/bptpart(ibox)
endif
endif
60 continue
endif
c sum of particles
npart=npart+bptpart(ibox)
Write(*,*) 'Particles = ',npart
cc
return
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
function swapi(x,dx)
c
c Version: 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 swapi, x, dx
swapi = x
1 if (swapi.lt.-dx) then
swapi = swapi + 2.0d0*dx
goto 1
end if
2 if (swapi.gt.dx) then
swapi = swapi - 2.0d0*dx
goto 2
end if
end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
Function Energie(alpha,max)
c
c Unit : calculate the energy
c Version : 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'coms.f'
integer i
real*8 alpha,max
real*8 Energie, E
E=0
Do 42 i=1,npart
E=E+sqrt((alpha**2)*(px(i)*px(i)+py(i)*py(i)+pz(i)*pz(i))+
$ fmass(i)**2)
42 continue
Energie=E-max
Return
End
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
Function Regula(me)
c
c Unit : Searches for the zero of the function
c Version : 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'coms.f'
real*8 Regula,xn,xu,x0,me,Energie
integer i
real*8 E1,E2,E3
c init
xu=0.d0
x0=3.d0
c main
Write(*,*) 'Regula is running!'
i=0
10 Continue
i=i+1
E1=Energie(x0,me)
E2=Energie(xu,me)
xn=x0-(E1*(x0-xu))/(E1-E2)
E3=Energie(xn,me)
IF ((E2*E3).LE.0) then
x0=xn
else
xu=xn
EndIF
IF ((ABS(x0-xu).GE.1.D-12).and.(i.le.1000).and.(
& ((E3.ge.1.D-12).or.(-E3.ge.1.D-12)))) goto 10
Regula=xn
End
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
real*8 function wallcoll(ipart,wall)
c
c Unit : Collisions with an imaginary wall
c Version : 1.0
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'coms.f'
include 'boxinc.f'
include 'options.f'
c var
real*8 betax,betay,betaz
real*8 ty,tz
real*8 tn
integer wall,ipart
integer wally,wallz
c Mainprogram
c init the variables
wall=0
tn=0
c velocity
betax=px(ipart)/p0(ipart)
betay=py(ipart)/p0(ipart)
betaz=pz(ipart)/p0(ipart)
c check which wall is reached next and sort by impact time
c wall presents the wall and tn the time
if (betax.lt.0) then
wall=-4
tn=(-lboxhalbe-rx(ipart))/(-max(-betax,1.d-13))
else
wall=-1
tn=((lboxhalbe-rx(ipart))/max(betax,1.d-13))
endif
if (betay.lt.0) then
wally=-5
ty=(-lboxhalbe-ry(ipart))/(-max(-betay,1.d-13))
else
wally=-2
ty=((lboxhalbe-ry(ipart))/max(betay,1.d-13))
endif
if(ty.lt.tn) then
tn=ty
wall=wally
endif
if (betaz.lt.0) then
wallz=-6
tz=(-lboxhalbe-rz(ipart))/(-max(-betaz,1.d-13))
else
wallz=-3
tz=((lboxhalbe-rz(ipart))/max(betaz,1.d-13))
endif
if(tz.lt.tn) then
tn=tz
wall=wallz
endif
c sets the time for the earliest wall collision