Advanced Computing Platform for Theoretical Physics

Skip to content

GitLab

Commit c927d9f7 authored by Tanguy Pierog's avatar Tanguy Pierog
Browse files

correct typo and comments and cross-section behavior in EPOS. 

Final distributed version 1.0


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/tags/crmc.v1.0@3725 c7a5e08c-de06-0410-9364-b41cf42a0b17
parents 81ba7591 5dc814b9
Hide whitespace changes
Inline Side-by-side
README
View file @ c927d9f7
CRMC v1.0 Last modifications 2013/03/26
CRMC v1.0 Last modifications 2013/04/11
**********************************************************
The Program "crmc"
......@@ -30,7 +30,7 @@ Reference : C. Baus, T. Pierog and R. Ulrich. To be published (2013)
(please ask to colin.baus@kit.edu when needed)
-- Hadronic interacton models :
-- Hadronic interaction models :
* Post LHC :
......@@ -141,11 +141,13 @@ bin/crmc -T -m6
The details of the run can be controlled by the file "crmc.param".
In this file the 2 first commented options are valid for EPOS 1.99 and EPOS LHC only.
By defaults EPOS models use a simplified treatment of QGP in events where the energy density is high enough (including in pp). If you
In this file the 3 first commented options are valid for EPOS 1.99 and EPOS LHC only.
- By defaults EPOS models use a simplified treatment of QGP in events where the energy density is high enough (including in pp). If you
uncomment the line "!switch fusion off" this will be disable and running time will be shorter (but data description will be worth since
the pysics model will be incomplete). In that mode EPOS is comparable to PYTHIA model (no final state interaction).
By default only the final particles are recorded in the output file like for other cosmic ray models.
the physics model will be incomplete). In that mode EPOS is comparable to PYTHIA model (no final state interaction).
- By default only the final particles are recorded in the output file like for other cosmic ray models.
Uncommenting "!set istmax 1" allows the user to have the full chain of mother/daughter from the beam to the final particles with EPOS models.
The outfile is at least 2 times larger but includes the decayed particles and some special intermediate particles between the beam and
the real particles which allow the user to know where the particles were generated. The ids of such particles are the following :
......@@ -156,12 +158,16 @@ the real particles which allow the user to know where the particles were generat
93 : remnant. Mother of particles coming from the beam remnants.
Each primary particle has a mother with id 9x. For technical reasons all particles with the same id and same momentum are in fact the same
initial object (cluster, string or remnant) which was splitted in different local pieces to keep the correct production vertexes of the final
initial object (cluster, string or remnant) which was split in different local pieces to keep the correct production vertexes of the final
particles. If a set of particles with the same id (91 or 92) and the same momentum has 2 different mothers (one from the projectile and one
from the target) it means that the string was produced by the interaction of this pair of nucleons or for a cluster that it was the closest
pair of nucleon which participate at the formation of the cluster (a cluster is formed by string pieces of different pairs but hepmc format
allows only one mother in that case).
- By default the cross-section is calculated by a numerical method with is valid only for h-p or h-A
(h being pion, kaon or nucleon) but not AB (nucleus-nucleus). If you uncomment "!set isigma 2" the inelastic
cross-section will be always correct but it takes several minutes to compute it (see Note on simulating Heavy Ion events).
**********************************************************
Note on simulating heavy ion (HI) events
******I****************************************************
......@@ -173,8 +179,14 @@ state interactions).
EPOS 1.99 and LHC include a simplified treatment of final state interaction but doesn't include full hydro
simulations as already possible in EPOS 2. It gives a good overall description of HI data but it should not be
used to interpret QGP related observables (flow, jet quenching, etc ...) since the model is oversimplified in
this distribution.
For detailed simulations of HI collisions within EPOS 2, please contact K. Werner (werner@subatech.in2p3.fr).
this distribution.
* For detailed simulations of HI collisions within EPOS 2, please contact K. Werner (werner@subatech.in2p3.fr).*
Concerning the cross-section calculation in EPOS, to avoid unnecessary long time calculation needed
for HI interactions, the default type of cross-section calculation is fast and good for hA but not reliable for
AB collisions (OK for p-Pb, pi-C but not for C-C or Pb-Pb for instance). So if you really need to have
a correct cross-section for PbPb scattering for example, you should uncomment the line "set isigma 2" in the
crmc.param file but it will take several minutes to compute the cross-section each time you start CRMC.
The current IO library of HepMC (2.06) has a limitation of the number of produced particles.
In particular, it will return an error message if the number of particles is larger than 10000.
......
src/crmc-aaa.f
View file @ c927d9f7
......@@ -223,7 +223,7 @@ c-----------------------------------------------------------------------
lhct=".lhc"
iadd=4
else
isigma=2 !use pseudo-MC cross section for nuclear xs (slow)
isigma=1 !use pseudo-MC cross section for nuclear xs (slow)
ionudi=3
lhct=""
iadd=0
......
src/crmc.param.in
View file @ c927d9f7
......@@ -6,6 +6,8 @@
!set istmax 1 !uncomment to get full daughter/mother particle chain with EPOS
!set isigma 2 !uncomment to get correct inelastic cross-section for heavy ions with EPOS
!!Set up particle Decays
!switch decay off !no decay at all
......
src/epos/epos-bas.f
View file @ c927d9f7
......@@ -1603,7 +1603,7 @@ c copy all mothers before the daughter
nrem1=nrem1+1
endif
if(isthep2(k).gt.0)
& print *,' ',k,idhep2(k),jmohep2(1,k),isthep2(k) jm1hep=iepo2hep(io)
& print *,' ',k,idhep2(k),jmohep2(1,k),isthep2(k)
& ,'from',ihep2epo2(k)
enddo
......@@ -6406,6 +6406,9 @@ c (from tabulation) for pA/AA
sigelaaa=0.
if(model.eq.1)then
if(isigma.ne.2)then
if(maproj.gt.5.and.matarg.gt.5)write(ifmt,*)
& 'Cross-section may be wrong, use isigma=2 instead ! ',
& '(if you care about it ...)'
c eposcrse depends of ionudi while eposinecrse corresponds to ionudi=1 always
sigineaa=eposinecrse(ekin,maproj,matarg,idtarg)
sigelaaa=eposelacrse(ekin,maproj,matarg,idtarg)
......@@ -6445,6 +6448,8 @@ c here cut is absorbtion xs : cut + 95 % of excited diff.
& 'Cross-section can not be calculated with ionudi=0'
endif
else
write(ifmt,*)
& 'Compute Cross-section (can take a while...)'
call crseaaEpos(sigtotaa,sigineaa,sigcutaa,sigelaaa)
endif
elseif(isigma.eq.2.and.matarg.gt.1)then
......
src/qgsjetII-04/qgsjet-II-04.f
View file @ c927d9f7
......@@ -52,8 +52,9 @@ C (mainly to TOTEM data on total/elastic pp cross sections); C
C remnant treatment for pion-hadron/nucleus collisions improved C
C C
C 18.03.2013 - small corrections for A-p cross-section calculation C
C 09.04.2013 - diffractive flag correction (no physics change) C
C C
C last modification: 18.03.2013 C
C last modification: 09.04.2013 C
C Version qgsjet-II-04 (for CONEX) C
C C
C small corrections to adapt to CORSIKA : 25.07.2012 by T.Pierog C
......@@ -326,7 +327,7 @@ c-------------------------------------------------
* /,' ','| Publication to be cited when using this program: |',
* /,' ','| S.Ostapchenko, PRD 83 (2011) 014018 |',
* /,' ','| |',
* /,' ','| last modification: 18.03.2013 |',
* /,' ','| last modification: 09.04.2013 |',
* /,' ','| |',
* /,' ','| Any modification has to be approved by the author|',
* /,' ','====================================================',
......@@ -465,7 +466,7 @@ c reading cross sections from the file
if(debug.ge.0)write (moniou,*)'done'
goto 10
elseif(.not.producetables)then
write(moniou,*) "Missing QGSDAT-II-04 file !"
write(moniou,*) "Missing QGSDAT-II-04 file !"
write(moniou,*) "Please correct the defined path ",
&"or force production ..."
stop
......@@ -1972,7 +1973,7 @@ c nuclear cross sections
close(2)
 
elseif(.not.producetables)then
write(moniou,*) "Missing sectnu-II-04 file !"
write(moniou,*) "Missing sectnu-II-04 file !"
write(moniou,*) "Please correct the defined path ",
&"or force production ..."
stop
......@@ -5959,121 +5960,179 @@ c diffraction (hadron-nucleus & nucleus-nucleus)
endif
endif
8 continue
c check diffractive cross sections (hadron-proton only)
c check diffractive cross sections !so060413-beg
9 jdiff=0 !non-diffractive
if(ia(1).eq.1.and.ia(2).eq.1.and.(nwp.ne.0.or.nwt.ne.0)
*.and.nqs(1).eq.0)then
if(lqa(1).eq.0.and.lqb(1).eq.0)then
if(nbpom.eq.0.or.npomin(1).eq.0)then
if(iwp(1).ge.2.and.iwt(1).lt.2)then
nqst=0
nint=0
if(nbpom.ne.0)then
do i=1,nbpom
nqst=nqst+nqs(i)
nint=nint+npomin(i)
enddo
endif
if((nwp.ne.0.or.nwt.ne.0).and.nqst.eq.0)then !not elastic nor ND
lqat=0
do ip=1,ia(1)
lqat=lqat+lqa(ip)
enddo
lqbt=0
do it=1,ia(2)
lqbt=lqbt+lqb(it)
enddo
iwpt=0
do ip=1,ia(1)
if(iwp(ip).eq.1)then
iwpt=1
goto 10
elseif(iwp(ip).ge.2)then
iwpt=2
endif
enddo
10 continue
iwtt=0
do it=1,ia(2)
if(iwt(it).eq.1)then
iwtt=1
goto 11
elseif(iwt(it).ge.2)then
iwtt=2
endif
enddo
11 continue
if(lqat.eq.0.and.lqbt.eq.0)then
if(nbpom.eq.0.or.nint.eq.0)then
if(iwpt.eq.2.and.iwtt.ne.2)then
jdiff=6 !SD(LM)-proj
elseif(iwp(1).lt.2.and.iwt(1).ge.2)then
elseif(iwpt.ne.2.and.iwtt.eq.2)then
jdiff=7 !SD(LM)-targ
elseif(iwp(1).ge.2.and.iwt(1).ge.2)then
elseif(iwpt.eq.2.and.iwtt.eq.2)then
jdiff=8 !DD(LM)
else
goto 12
endif
else
if(iwp(1).lt.2.and.iwt(1).lt.2)then
else
goto 14
endif
else
if(iwpt.ne.2.and.iwtt.ne.2)then
jdiff=9 !CD(DPE)
else
else
jdiff=10 !CD+LMD
endif
endif
elseif(lqa(1).gt.0.and.lqb(1).eq.0.and.iwt(1).lt.2)then
endif
endif
elseif(lqat.gt.0.and.lqbt.eq.0.and.iwtt.ne.2)then
jdiff=1 !SD(HM)-proj
elseif(lqa(1).eq.0.and.lqb(1).gt.0.and.iwp(1).lt.2)then
elseif(lqat.eq.0.and.lqbt.gt.0.and.iwpt.ne.2)then
jdiff=2 !SD(HM)-targ
elseif(lqa(1).gt.0.and.lqb(1).eq.0.and.iwt(1).ge.2)then
elseif(lqat.gt.0.and.lqbt.eq.0.and.iwtt.eq.2)then
jdiff=3 !DD(LHM)-proj
elseif(lqa(1).eq.0.and.lqb(1).gt.0.and.iwp(1).ge.2)then
elseif(lqat.eq.0.and.lqbt.gt.0.and.iwpt.eq.2)then
jdiff=4 !DD(LHM)-targ
 
elseif(lqa(1).gt.0.and.lqb(1).gt.0)then
if(npompr(1).eq.0)stop'problem with npompr!!!'
xrapmax(1)=1.d0
do i=1,npompr(1)
xrapmax(1)=min(xrapmax(1),1.d0/xpompr(i,1)/scm)
enddo
if(npomtg(1).eq.0)stop'problem with npomtg!!!'
xrapmin(1)=1.d0/scm
do i=1,npomtg(1)
xrapmin(1)=max(xrapmin(1),xpomtg(i,1))
enddo
if(xrapmin(1).gt..999d0*xrapmax(1))goto 12
nraps=1
irap=1
11 if(nraps.gt.90)stop'nraps>90'
if(npomin(1).gt.0)then
do i=1,npomin(1)
if(xpomin(i,1).lt..999d0*xrapmax(irap)
* .and.xpopin(i,1).gt.1.001d0*xrapmin(irap))then
if(xpomin(i,1).lt.1.001d0*xrapmin(irap)
* .and.xpopin(i,1).gt..999d0*xrapmax(irap))then
nraps=nraps-1
if(nraps.eq.0)goto 12
irap=irap-1
goto 11
elseif(xpopin(i,1).gt..999d0*xrapmax(irap))then
xrapmax(irap)=xpomin(i,1)
if(xrapmin(irap).gt..999d0*xrapmax(irap))then
nraps=nraps-1
if(nraps.eq.0)goto 12
irap=irap-1
goto 11
endif
elseif(xpomin(i,1).lt.1.001d0*xrapmin(irap))then
xrapmin(irap)=xpopin(i,1)
if(xrapmin(irap).gt..999d0*xrapmax(irap))then
nraps=nraps-1
if(nraps.eq.0)goto 12
irap=irap-1
goto 11
endif
else
xrapmin(irap+1)=xrapmin(irap)
xrapmin(irap)=xpopin(i,1)
xrapmax(irap+1)=xpomin(i,1)
if(xrapmin(irap).lt..999d0*xrapmax(irap)
* .and.xrapmin(irap+1).lt..999d0*xrapmax(irap+1))then
irap=irap+1
nraps=nraps+1
goto 11
elseif(xrapmin(irap).lt..999d0*xrapmax(irap))then
goto 11
elseif(xrapmin(irap+1).lt..999d0*xrapmax(irap+1))then
xrapmin(irap)=xrapmin(irap+1)
xrapmax(irap)=xrapmax(irap+1)
goto 11
else
nraps=nraps-1
if(nraps.eq.0)goto 12
irap=irap-1
goto 11
endif
endif
endif
enddo !end of npin-loop
endif
elseif(lqat.gt.0.and.lqbt.gt.0)then
if(nbpom.eq.0)stop'problem with nbpom!!!'
xrapmax(1)=1.d0
xrapmin(1)=1.d0/scm
do ibpom=1,nbpom
if(npompr(ibpom).gt.0)then
do i=1,npompr(ibpom)
ip=ilpr(i,ibpom)
lpom=lnpr(i,ibpom)
xrapmax(1)=min(xrapmax(1),1.d0/xpompr(lpom,ip)/scm)
enddo
endif
if(npomtg(ibpom).gt.0)then
do i=1,npomtg(ibpom)
it=iltg(i,ibpom)
lpom=lntg(i,ibpom)
xrapmin(1)=max(xrapmin(1),xpomtg(lpom,it))
enddo
endif
enddo
if(xrapmin(1).gt..999d0*xrapmax(1))goto 14
nraps=1
12 if(nraps.gt.90)stop'nraps>90'
do ibpom=1,nbpom
if(npomin(ibpom).gt.0)then
do i=1,npomin(ibpom)
if(nraps.eq.1)then
if(1.d0/scm/xpomin(i,ibpom).lt..999d0*xrapmax(1)
* .and.xpopin(i,ibpom).gt.1.001d0*xrapmin(1))then !rap-gaps changed
if(1.d0/scm/xpomin(i,ibpom).lt.1.001d0*xrapmin(1)
* .and.xpopin(i,ibpom).gt..999d0*xrapmax(1))then !no rap-gap (filled)
goto 14
elseif(xpopin(i,ibpom).gt..999d0*xrapmax(1))then
xrapmax(1)=1.d0/scm/xpomin(i,ibpom)
elseif(1.d0/scm/xpomin(i,ibpom).lt.1.001d0*xrapmin(1))then
xrapmin(1)=xpopin(i,ibpom)
else
xrapmin(2)=xrapmin(1)
xrapmin(1)=xpopin(i,ibpom)
xrapmax(2)=1.d0/scm/xpomin(i,ibpom)
nraps=2
goto 12
endif
endif
else
if(1.d0/scm/xpomin(i,ibpom).lt..999d0*xrapmax(1)
* .and.xpopin(i,ibpom).gt.1.001d0*xrapmin(nraps))then !rap-gaps changed
if(1.d0/scm/xpomin(i,ibpom).lt.1.001d0*xrapmin(nraps)
* .and.xpopin(i,ibpom).gt..999d0*xrapmax(1))then !no rap-gaps (filled)
goto 14
else
do irap=1,nraps
if(xpopin(i,ibpom).gt..999d0*xrapmax(irap).and.1.d0/scm
* /xpomin(i,ibpom).lt.1.001d0*xrapmin(irap))then !gap filled
if(irap.lt.nraps)then
do j=irap,nraps-1
xrapmax(j)=xrapmax(j+1)
xrapmin(j)=xrapmin(j+1)
enddo
endif
nraps=nraps-1
goto 12
elseif(xpopin(i,ibpom).gt..999d0*xrapmax(irap))then
xrapmax(irap)=min(1.d0/scm/xpomin(i,ibpom)
* ,xrapmax(irap))
elseif(1.d0/scm/xpomin(i,ibpom)
* .lt.1.001d0*xrapmin(irap))then
xrapmin(irap)=max(xpopin(i,ibpom),xrapmin(irap))
elseif(1.d0/scm/xpomin(i,ibpom).gt.xrapmin(irap)
* .and.xpopin(i,ibpom).lt.xrapmax(irap))then
xrapmin(irap)=max(xpopin(i,ibpom),xrapmin(irap))
if(irap.lt.nraps)then
do j=1,nraps-irap
xrapmax(nraps-j+2)=xrapmax(nraps-j+1)
xrapmin(nraps-j+2)=xrapmin(nraps-j+1)
enddo
endif
xrapmin(irap+1)=xrapmin(irap)
xrapmin(irap)=xpopin(i,ibpom)
xrapmax(irap+1)=1.d0/scm/xpomin(i,ibpom)
nraps=nraps+1
goto 12
endif
enddo !end of irap-loop
endif
endif
endif
enddo !end of npin-loop
endif
enddo !end of ibpom-loop
jdiff=5 !DD(HM)
endif
endif !end of diffr. check
12 bdiff=b
14 bdiff=b
 
ctp define collision type
typevt=0 !no interaction
if(ia(1).eq.1.and.ia(2).eq.1.and.(nwp.gt.0.or.nwt.gt.0))then !only for h-h
if(nwp.gt.0.or.nwt.gt.0)then !so060413-end
if(jdiff.eq.0)then !ND (no rap-gaps)
typevt=1
elseif(jdiff.eq.8.or.jdiff.eq.10.or.
* (jdiff.gt.2.and.jdiff.lt.6))then !DD + (CD+LMD)
typevt=2
typevt=2
elseif(jdiff.eq.1.or.jdiff.eq.6)then !SD pro
typevt=4
typevt=4
elseif(jdiff.eq.2.or.jdiff.eq.7)then !SD tar
typevt=-4
typevt=-4
elseif(jdiff.eq.9)then !CD
typevt=3
else
......@@ -16577,7 +16636,7 @@ c-----------------------------------------------------------------------
 
if(debug.ge.3)write (moniou,201)ic,ep0,nsh
nsh=nsh+1
nstprev = nsh
 
if(nsh.gt.nptmax)stop'increase nptmax!!!'
......@@ -16585,9 +16644,9 @@ c-----------------------------------------------------------------------
do i=1,4
ep(i)=ep0(i)
enddo
c call qgtran(ep,ey0,1)
if(iab.eq.7.or.iab.eq.8)then !delta++(-)
call qgdec2(ep,ep1,ep2,dmmin(2)**2,amn**2,am0**2)
ich(nsh)=ic-5*ic/iab
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Provided by Institute of Theoretical Physics, Chinese Academy of Sciences. Please contact support@itp.ac.cn for more information.