IAP GITLAB

Commit befa178d authored by Colin Baus's avatar Colin Baus

enable AA cross section calculation for dpmjet (requires isigma=2 option)

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4396 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent f258dbab
......@@ -10,11 +10,11 @@ OPTION (__QGSJET01__ "Build with model" OFF)
OPTION (__GHEISHA__ "Build with model" OFF)
OPTION (__PYTHIA__ "Build with model" OFF)
OPTION (__HIJING__ "Build with model" OFF) #doesn't work with gfortran !
OPTION (__SIBYLL__ "Build with model" OFF)
OPTION (__SIBYLL__ "Build with model" ON)
OPTION (__PHOJET__ "Build with model" OFF)
OPTION (__DPMJET__ "Build with model" ON)
OPTION (__DPMJET__ "Build with model" OFF)
OPTION (__QGSJETII03__ "Build with model" OFF)
OPTION (__QGSJETII04__ "Build with model" OFF)
OPTION (__QGSJETII04__ "Build with model" ON)
######################################ONLY EDIT THIS######################################
if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
......@@ -29,7 +29,7 @@ MESSAGE(STATUS "build ${CMAKE_BINARY_DIR}")
## The version number
SET (CRMC_VERSION_MAJOR 1)
SET (CRMC_VERSION_MINOR 5b)
SET (CRMC_VERSION_MINOR 4)
##Compiler options
# hotfix for lxbatch not sure why he selects wrong compiler
......
......@@ -199,6 +199,7 @@ for HI interactions, the default type of cross-section calculation is fast and g
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 number of loops is hard-coded (20000 check models.F) and does not depend on "-n" option.
The HepMC (2.06) IO library has a default limitation of the number of produced particles of 10000.
Thus, events with more than that number of secondaries are either truncated or skipped. For heavy ion
......
......@@ -2,11 +2,11 @@
!! a line starting with "!" is not read by the program
!switch fusion off !nuclear effects due to high density (QGP) in EPOS
!more realistic but slow (can be switched off)
!more realistic but slow (can be switched off)
!set istmax 1 !include virtual mother particles with EPOS to identify particle source
!set istmax 1 !include virtual mother particles with EPOS to identify particle source
!set isigma 2 !uncomment to get correct inelastic cross-section for heavy ions with EPOS
!set isigma 2 !uncomment to get correct inelastic cross-section for heavy ions with EPOS, QGSJET and DPMJET
!!Set up particle Decays
!switch decay off !no decay at all
......
......@@ -444,17 +444,38 @@ C CALL DT_EVTOUT(4)
c------------------------------------------------------------------------------
subroutine GetDPMJETSigma(stot,sine,sela)
c------------------------------------------------------------------------------
c return inelastic cross section of DPMJET
c dummy until returned from phojet
stot = 0.
sela = 0.
sine = 0.
end
c------------------------------------------------------------------------------
subroutine GetDPMJETSigmaAA(niter,stotaa,sineaa,selaaa)
c------------------------------------------------------------------------------
c return AA cross sections of DPMJET
c but also fill other cross section variables
c some dpmjet functions must be called for this
c this is experimental (no evaporation)
c in dpmjet the control card xs-table starts to calculate the table and then quits
c it calls subroutine dt_xstabl with the 6 parameters for binning and energy and binning and q
c in this subroutine the actual table COMMON/DTPART/ is filled
c this happens with dt_xsglau. emulsion is off? at least ncompo is 0 => no loops over table
c instead always (1,1,1) contians the value
c xprod seems to be prod + quasi-ela. this is why we subtract it again
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
include 'epos.inc'
real stot,sine,sela
integer niter !input
real stot,sine,sela,stotaa,sineaa,selaaa !subroutine return values
PARAMETER (TINY10=1.0D-10,TINY2=1.0D-2,ZERO=0.0D0,DLARGE=1.0D10,
& OHALF=0.5D0,ONE=1.0D0,TWO=2.0D0)
* properties of interacting particles
COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
* Glauber formalism: cross sections
PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
......@@ -468,17 +489,26 @@ c this is experimental (no evaporation)
& XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
& BSLOPE,NEBINI,NQBINI
* Glauber formalism: flags and parameters for statistics
LOGICAL LPROD
CHARACTER*8 CGLB
COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
JSTATB = niter !loops for glauber calculation !limited by ksiteb = 50
LPROD = .FALSE. !if set to false it enables the calculation of tot, el... otherwise only prod
XI = 0D0
Q2 = 0D0
CALL DT_XSGLAU(IP,IT,IJPROJ,XI,Q2,dble(ECMS),1,1,-1)
c #######from dt_xstabl##########
CALL DT_XSGLAU(IP,IT,IJPROJ,ZERO,ZERO,dble(ECMS),1,1,-1)
c #######from dt_xstabl##########
XTOT = XSTOT(1,1,1)
XELA = XSELA(1,1,1)
XINE = XSPRO(1,1,1)
XPRO = XTOT - XELA - XSQEP(1,1,1) - XSQET(1,1,1) - XSQE2(1,1,1)
stot = sngl(XTOT)
sine = sngl(XTOT-XELA)
sela = sngl(XELA)
stotaa = sngl(XTOT)
sineaa = sngl(XPRO)
selaaa = sngl(XELA)
end
......
......@@ -6594,12 +6594,7 @@ c GDD - double diffractive cross section
sigineaa=urqincs
elseif(model.eq.12)then !for DPMJet
call dpmjetSIGMA(stot,sine,sela)
sigtot=stot
sigine=sine
sigela=sela
call dpmjetSIGMA(sigtot,sigine,sigela)
endif
if(isigma.ge.1)then !===============!
......@@ -6619,8 +6614,6 @@ c (from tabulation) for pA/AA
if((isigma.eq.2.and.noebin.ge.0)
& .or..not.(maproj.eq.1.and.matarg.eq.1))then
sigtotaa=0.
sigelaaa=0.
if(model.eq.1)then
if(isigma.ne.2)then
if(maproj.gt.5.and.matarg.gt.5)then
......@@ -6669,18 +6662,18 @@ c excited diffractive part
sigelaaa=sigelaaa+sigqela
c here cut is absorbtion xs : cut + 95 % of excited diff.
sigcutaa=sigcutaa+0.95*difpart*sdpart
elseif(ionudi.eq.0)then
elseif(ionudi.eq.0)then
write(ifmt,*)
& 'Cross-section can not be calculated with ionudi=0'
endif
else
write(ifmt,*)
& 'Compute EPOS Cross-section (can take a while...)'
& 'Computing EPOS cross-section (can take a while...)'
call crseaaEpos(sigtotaa,sigineaa,sigcutaa,sigelaaa)
endif
elseif(isigma.eq.2.and.matarg.gt.1)then
write(ifmt,*)
& 'Compute model Cross-section (can take a while...)'
& 'Computing model cross-section (can take a while...)'
call crseaaModel(sigtotaa,sigineaa,sigcutaa,sigelaaa)
endif
if(ish.ge.1.and.noebin.ge.0)
......
......@@ -330,7 +330,7 @@
stop'please compile with requested model'
print *, GTOT,GPROD,GABS,GDD,GQEL,GCOH,sigt,sigi,sigc,sige,e0 !get rid of unused warning
#else
NITER=5000
NITER=20000
if(idtarg.eq.0)then
e0=dble(elab)
icp=idtrafo('nxs','qgs',idproj)
......@@ -378,7 +378,7 @@
#ifndef __QGSJETII__
stop'please compile with requested model'
#else
NITER=5000
NITER=20000
if(idtarg.eq.0)then
e0=dble(elab)
icp=idtrafo('nxs','qgs',idproj)
......@@ -398,6 +398,14 @@
sige=0.
sigc=0.
sigi=flucrse(ekin,maproj,matarg,idtargin)
#endif
elseif(model.eq.12)then
#ifndef __DPMJET__
stop'please compile with requested model'
#else
NITER=20000
call GetDPMJETSigmaAA(niter,sigt,sigi,sige)
sigc=0
#endif
else
sigt=0.
......@@ -522,16 +530,13 @@ C cross sections
#endif
end
subroutine dpmjetSIGMA(stot,sine,sela)
subroutine dpmjetSIGMA(stot,sine,sela,stotaa,sineaa,selaaa)
#ifndef __DPMJET__
stop'please compile with requested model'
print *, stot,scut,sine,sela,slela,ssd !get rid of unused warning
print *, stot,sine,sela,stotaa,sineaa,selaaa !get rid of unused warning
#else
call GetDPMJETSigma(tot,ine,ela)
stot=tot
sine=ine
sela=ela
call GetDPMJETSigma(stot,sine,sela)
#endif
end
......
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