IAP GITLAB

Commit 5b9dc6c3 authored by Tanguy Pierog's avatar Tanguy Pierog

update new DPMJETIII.2017-1 to a properly running version


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@6207 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 06c0246b
......@@ -223,8 +223,7 @@ ENDIF (__PHOJET__)
IF (__DPMJET__)
LIST(APPEND TABS dpmjpar.dat)
LIST(APPEND TABS dpmCT14LL.pds)
LIST(APPEND TABS glaubint.glb)
LIST(APPEND TABS glaubtar.glb)
LIST(APPEND TABS conextar.glb)
ADD_SUBDIRECTORY ("${PROJECT_SOURCE_DIR}/src/dpmjet/3.2017-1")
SET(STATIC_LIBS ${STATIC_LIBS} Dpmjet)
get_property(HELPER SOURCE src/models.F PROPERTY COMPILE_FLAGS)
......
This source diff could not be displayed because it is too large. You can view the blob instead.
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -18782,7 +18782,7 @@ cdh
* new patch for pre-initialized variable projectile/target/energy runs
IF (IOGLB.EQ.100) THEN
if (NIDX.eq.-2) then
if (NIDX.eq.-10) then
c write(0,*) ' -- dt_glaube -- skip call into dt_glbset()'
else
......@@ -18929,6 +18929,7 @@ c DATA ECMXPI / 100000.0D0 /
*--------------------------------------------------------------------------
* general initializations
*
LPRI = 20
LOUT = 6
* steps in projectile mass number for initialization
......@@ -19220,6 +19221,7 @@ cdh datadir for path to the data sets to be read in by dpmjet/phojet
*
* read data from file
*
IF (MODE.EQ.0) THEN
IF (LREAD) RETURN
......@@ -19403,7 +19405,7 @@ c modification for use with corsika using path to data file in DATADIR
*
ELSE
WRITE(*,*) 'DT_GLBSET called for ',IDPROJ,NA,NB,ELAB,MODE
c WRITE(*,*) 'DT_GLBSET called for ',IDPROJ,NA,NB,ELAB,MODE
*
* check for type of projectile and set index-offset to entry in
* Glauber data array correspondingly
......@@ -21549,14 +21551,14 @@ C WHAT(2) = 0
C CODEWD = 'START '
C GOTO 900
C ENDIF
IF (NCASES.EQ.-1) THEN
IF (NCASES.LE.-1) THEN !variable energy with air (TP20170630)
IP = NPMASS
IPZ = NPCHAR
IT = NTMASS
ITZ = NTCHAR
PPN = EPNSAV
VARELO = 10.D0
VAREHI = EPN*1.1D0
VAREHI = PPN*1.1D0
EPN = ZERO
CMENER = ZERO
LEINP = .TRUE.
......@@ -21565,7 +21567,27 @@ C ENDIF
WHAT(2) = 0
CODEWD = 'START '
LEVPRT = .TRUE.
IF(NCASES.EQ.-2)THEN
IOGLB = 0 ! don't use glauber tables
ELSE
IOGLB = 100 ! use glauber tables
ENDIF
GOTO 900
ELSEIF (NCASES.EQ.-100) THEN !make glauber table
if(ifirst.ne.1)stop
ifirst=2
IP = NPMASS
IPZ = NPCHAR
IT = NTMASS
ITZ = NTCHAR
PPN = EPNSAV
WHAT(1) = 10.D0
WHAT(2) = PPN*1.1D0
WHAT(3) = 7d0
WHAT(4) = 56d0
WHAT(5) = 7d0
CODEWD = 'GLAUB-INI'
goto 565
ENDIF
* read control card from input-unit LINP
READ(LINP,'(A78)',END=9999) CLINE
......@@ -26696,7 +26718,7 @@ Cf2py intent(out) IREJ
ELSE IF ((IT.EQ.1).AND.(ITZ.EQ.0)) THEN
IJTARG = 8
ELSE IF ((IT.EQ.1).AND.(ITZ.EQ.-1)) THEN
IJTARG = -2
IJTARG = 2
END IF
IBTARG = IIBAR(IJTARG)
......@@ -33931,7 +33953,7 @@ C STOP
1003 FORMAT(9X,F6.1,9X,F6.2,8X,F8.3,11X,F8.3)
1 CONTINUE
IVEOUT = 1
ELSE
ELSE
* hadron/photon/nucleus-nucleus
IF ((ABS(VAREHI).GT.ZERO).AND.
& (ABS(VAREHI).GT.ABS(VARELO))) THEN
......@@ -34170,7 +34192,6 @@ CDECK ID>, DT_SIGEMU
LOGICAL LPHOIN
COMMON /DTMODL/ CMODEL(4),ELOJET,MCGENE,LPHOIN
IF (MCGENE.NE.4) THEN
IF (LPRI.GT.4)
& WRITE(LOUT,'(A)') ' DT_SIGEMU: Combined cross sections'
......@@ -38722,7 +38743,6 @@ C DIBETA = SDIF1/STOT
SDQE = ZERO
SDQE2 = ZERO
FACN = ONE/DBLE(NSTATB)
IPNT = 0
RPNT = ZERO
......@@ -38759,6 +38779,7 @@ C CALL GSET(XAMLO,XAMHI,NPOINT,ABSZX,WEIGHT)
CALL DT_POILIK(NB,NTARG,ECMNN(IE),Q2,IPNT,RPNT,1)
ENDIF
* read pre-initialized profile-function from file
IF (IOGLB.EQ.1) THEN
READ(LDAT,'(5I10,E15.5)') KJPROJ,IA,IB,ISTATB,ISITEB,DUM
......@@ -33,37 +33,81 @@ c emulsion treatment
PARAMETER (NCOMPX=100,NEB=8,NQB= 5,KSITEB=50)
COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
& NCOMPO,IEMUL
LPRI=1+ish
c nucleon-nucleon event-generator
PARAMETER ( MXFFBK = 6 )
PARAMETER ( MXZFBK = 10 )
PARAMETER ( MXNFBK = 12 )
PARAMETER ( MXAFBK = 16 )
PARAMETER ( MXPSST = 700 )
PARAMETER ( MXPSFB = 43000 )
PARAMETER ( MXASST = 25)
PARAMETER ( NXAFBK = MXAFBK + 1 )
PARAMETER ( NXZFBK = INT(MXZFBK + MXFFBK / 3
& + MXASST - NXAFBK) )
PARAMETER ( NXNFBK = INT(MXNFBK + MXFFBK / 3
& + MXASST - NXAFBK) )
LOGICAL LFRMBK, LNCMSS
COMMON /FRBKCM/ AMUFBK, EEXFBK(MXPSST), AMFRBK(MXPSST),
& WEIFBK(MXPSST) ,GAMFBK(MXPSST), EXFRBK(MXPSFB),
& SDMFBK(MXPSFB), COUFBK(MXPSFB), CENFBK(MXPSFB),
& EXMXFB, R0FRBK, R0CFBK, C1CFBK, C2CFBK, FRBKLS,
& IFRBKN(MXPSST), IFRBKZ(MXPSST), IFBKSP(MXPSST),
& IFBKPR(MXPSST), IFBKST(MXPSST), IFBKLV(MXPSST),
& IPSIND(0:NXNFBK,0:NXZFBK,2), JPSIND(0:MXASST),
& IFBIND(0:NXNFBK,0:NXZFBK,2), JFBIND(0:NXAFBK),
& IFBCHA(9,MXPSFB), IPOSST,IPOSFB,IFBSTF, IFBPSF,
& IFBPSI, IFBFRB, IFBCHN, IFBNC1, IFBNC2, NBUFBK,
& LFRMBK, LNCMSS
c properties of interacting particles
LOGICAL LDIFFR,LINCTV,LEVPRT,LHEAVY,LDEEXG,LGDHPR,LPREEX,
& LHLFIX,LPRFIX,LPARWV,LPOWER,LSNGCH,LSCHDF,LHADRI,
& LNUCRI,LPEANU,LEVBME,LPHDRC,LATMSS,LISMRS,LCHDCY,
& LCHDCR,LMLCCR,LRVKIN,LVP2XX,LV2XNW,LNWV2X,LEVFIN
PARAMETER (NALLWP = 64)
COMMON /PAREVT/ DPOWER, FSPRD0, FSHPFN, RN1GSC, RN2GSC, RNSWTC,
& LDIFFR(NALLWP), LPOWER, LINCTV, LEVPRT, LHEAVY,
& LDEEXG, LGDHPR, LPREEX, LHLFIX, LPRFIX, LPARWV,
& LSNGCH, LSCHDF, LHADRI, LNUCRI, LPEANU, LEVBME,
& LPHDRC, LATMSS, LISMRS, LCHDCY, LCHDCR, LMLCCR,
& LRVKIN, LVP2XX, LV2XNW, LNWV2X, LEVFIN
LPRI = 1+ish
LOUT = ifch
call utpri('inidpm',ish,ishini,6)
write(ifmt,'(a,i6)')'initialize DPMjet ...'
LOUT = ifch
IDFRAME = 2 !dpmjet iframe variable. nucleon-nucleon frame
NDPMEVENT = 0 !event number needed in dpmjet
IFUSION = 1 ! ENABLE FUSION IN NUCLEUS-NUCLEUS COLLISIONS
if(idtargin.eq.0)then !Air
ITRSPT = 1 ! use glauber tables
CGLB = 'glaubtar'
C GLAUBER DATA SET CONTAINS TARGET EMULSION WITH 14N, 16O, 40AR
NCOMPO = 3
NB = 40
IEMUMA(1) = 14 ! NITROGEWN ARGET
IEMUCH(1) = 7
EMUFRA(1) = airwnxs(1)
IEMUMA(2) = 16 ! OXYGEN TARGET
IEMUCH(2) = 8
EMUFRA(2) = airwnxs(2)
IEMUMA(3) = 40 ! ARGON TARGET
IEMUCH(3) = 18
EMUFRA(3) = airwnxs(3)
else
ITRSPT = 0 ! don't use glauber tables
endif
egymin = 10. !min energy for model
egymax = 2.e6 !max energy for model
ICASCA = 1 ! NO STATISTICAL INFO NEEDED IN SHOWER SIMULATION
C EVAPORATION MODULES NOT AVAILABLE WITH THIS VERSION
LEVPRT = .FALSE.
LDEEXG = .FALSE.
LHEAVY = .FALSE.
LFRMBK = .FALSE.
CGLB = 'conextar'
c DO NOT USE EMULSION : CREATE SOME PROBLEM WITH ENERGY CONSERVATION (BUG IN INTERFACE ?)
C BUT TABULATION CAN BE USED ANY WAY ....
c if(idtargin.eq.0)then !Air
c CGLB = 'conextar'
cC GLAUBER DATA SET CONTAINS TARGET EMULSION WITH 14N, 16O, 40AR
c IEMUL = 1
c NCOMPO = 3
c NB = 40
c IEMUMA(1) = 14 ! NITROGEWN ARGET
c IEMUCH(1) = 7
c EMUFRA(1) = airwnxs(1)
c IEMUMA(2) = 16 ! OXYGEN TARGET
c IEMUCH(2) = 8
c EMUFRA(2) = airwnxs(2)
c IEMUMA(3) = 40 ! ARGON TARGET
c IEMUCH(3) = 18
c EMUFRA(3) = airwnxs(3)
c endif
egymin = 4. !min energy for model
egymax = 2.e6 !max energy for model
BIMIN = bminim !impact parameter range
BIMAX = bmaxim
......@@ -94,6 +138,11 @@ c-----------------------------------------------------------------------
COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
& IICH(210),IIBAR(210),K1(210),K2(210)
C model switches and parameters
CHARACTER*8 MDLNA
INTEGER ISWMDL,IPAMDL
DOUBLE PRECISION PARMDL
COMMON /POMDLS/ MDLNA(50),ISWMDL(50),PARMDL(400),IPAMDL(400)
c fresh common block
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
DOUBLE PRECISION EPROJ,EINID
......@@ -126,8 +175,7 @@ c will be treated as nucleus in DT_INIT and NPMASS,... will be used
IDPDG = 0
IDP = 0
ENDIF
if(idtarg.ne.1120)then
if(idtarg.ne.1120.and.idtarg.ne.0)then
call utstop("Only proton or nucleus as target !&")
else
NTCHAR = abs(NTCHAR)
......@@ -140,7 +188,6 @@ c will be treated as nucleus in DT_INIT and NPMASS,... will be used
C general initialization
NCASES = -1 !skip reading steering cards
if(noebin.gt.1)then !if more than one energy bin, initialize for
if(iologl.ge.0)then
EINID = dble(engmax)+1d0 !highest kin energy
......@@ -150,13 +197,17 @@ C general initialization
else
EINID = EPROJ
endif
NCASES = -2 !skip reading steering cards and fixed energy
if(idtargin.eq.0)NCASES = -1
c NCASES = -100 !make cross-section table (using iron as projectile and idtarg=0 for air)
IGLAU = 0
c write(ifch,*)'dt_init',NCASES,EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR
c +,IDP,IGLAU
CALL DT_INIT(NCASES,EINID,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
c CALL DT_INIT(NCASES,EINID,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
CALL DT_DTUINI(NCASES,EINID,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IEMU)
IPAMDL(179)=1000000000 !prevent Phojet to stop because of inconsistencies
* pre-initialization of profile function for nuclear collisions*
c ECM=dble(ecms)
......@@ -186,26 +237,35 @@ c + ,(PMAS(i,2).LT.PARP(41))
c don't touch decays with too small width
if(PMAS(i,2).LT.PARP(41)) then
go to 100
go to 101
endif
c because epos var ifrade is not 0 it will automatically use epos to decay when switched off in pythia
if(ish.ge.2) write(ifch,*)
+ "CF Id:",i,PMAS(i,1),PMAS(i,2),
+ "PDG ID:",KCHG(I,4),
+ "(decay in epos)"
+ "(decay in EPOS)"
MDCY(i,1)=0
100 enddo !loop over 100-500 particle ids
c initialize cross-sections by calling epos-bas.f function -> models.F -> dpmjet_epos.f
goto 101
c if unkown particle in EPOS, decay it in DPMJET
100 if(MDCY(i,1).ne.1.and.ish.ge.2) write(ifch,*)
+ "CF Id:",i,PMAS(i,1),PMAS(i,2),
+ "PDG ID:",KCHG(I,4),
+ "(decay in DPMJET)"
MDCY(i,1)=1
c 101 write(ifch,*)
c + "CF Id:",i,PMAS(i,1),PMAS(i,2),
c + "PDG ID:",KCHG(I,4),MDCY(i,1),istatus
c enddo
101 enddo !loop over 100-500 particle ids
c initialize cross-sections by calling epos-bas.f function -> models.f -> dpmjet-crmc.f
endif
call xsigma
dpmincs=sigineaa
c istmax=0 !uncomment only if you want to be sure to have only final particles
if(ish.ge.2)write(ifch,*)
& 'DPMJET used with (Elab,idproj,idtarg,xs) ',EPROJ,idproj,idtarg
& ,dpmincs,elab,ecms,engy,ekin
......@@ -246,6 +306,11 @@ c-----------------------------------------------------------------------
c event number
COMMON /DTEVNO/ NDPMEVENT,ICASCA
c emulsion treatment
PARAMETER (NCOMPX=100,NEB=8,NQB= 5,KSITEB=50)
COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
& NCOMPO,IEMUL
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
INTEGER IREJ, KKMAT
DOUBLE PRECISION EPROJ
......@@ -261,21 +326,30 @@ c event number
if(a.gt.0..and.rangen().gt.dpmincs/10./a)goto 1001 !no interaction
if(ish.ge.3)call alist('Determine DPMJET Production&',0,0)
c call conre !! adds beam particles that are known from intial call of crmc
c call conwr
call conre !! adds beam particles that are known from intial call of crmc
call conwr
IREJ = 0
KKMAT = -1
if (IEMUL.gt.0)then
CALL DT_GETEMU(NTMASS,NTCHAR,KKMAT,0)
KKMAT = -KKMAT !tables already in memory
endif
C call dpmjet event simulation
NDPMEVENT = NDPMEVENT+1 !needs to be updated for internal init. sets NEVHKK
c projecticle and target can be changed without reinitialization of the event
IF (IDPDG.EQ.0) THEN
IDP = 1
ELSE
if(irdmpr.ne.0.and.laproj.eq.-1)
. IDPDG = idtrafo('nxs','pdg',idproj)
IDP = IDT_ICIHAD(IDPDG)
ENDIF
c write(ifch,*)'dt_kkinit',NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPROJ
c &, KKMAT
NPMASS = maproj
NPCHAR = laproj
NTMASS = matarg
NTCHAR = latarg
CALL dt_kkinc(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPROJ, KKMAT,IREJ)
IF(IREJ.NE.0) THEN
WRITE(ifch,'(a)')
......@@ -334,6 +408,7 @@ c PLEASE FILL!!
GOTO 1001
ENDIF
anintine=anintine+1.
c esum=0.d0
do k=1,NHKK
......@@ -346,25 +421,42 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
$ ,JMOHKK(1,k),JMOHKK(2,k),JDAHKK(1,k),JDAHKK(2,k)
nptl=nptl+1 !! add 1 particle to stack
IF(nptl.GE.mxptl) THEN
if(ish.ge.2)WRITE(ifch,'(a)')
& 'mxptl too small for event. Skipping'
GOTO 1001
ENDIF
anintine=anintine+1.
nptlhep(k)=nptl
IF(ISTHKK(k).GE.1 .AND. ISTHKK(k).LE.2 .AND. ic.NE.99999 )THEN! if final particle
ist=ISTHKK(k)-1 !!0 means last generation other codes are e.g. for pomerons, remnants...
elseif(ic.eq.99999)then
ist=31 !convert to status for epos string
elseif(ic.eq.66666)then
ist=11 !convert to status for epos cluster
else
ist=abs(ISTHKK(k)) !!0 means last generation other codes are e.g
endif
am=PHKK(5,k)
if(ic.eq.80000)then !nucleus
id=1000000000+IDXRES(k)*10000+IDRES(k)*10
if(ist.eq.1000)then
ist=0
c esum=esum+PHKK(4,k)
if(abs(am-IDXRES(k)*0.938d0).gt.1d0)am=IDRES(k)*0.938d0
else
ist=19
endif
elseif(ic.eq.66666)then !cascade
id=80000000
else
id=idtrafo('pdg','nxs',ic) !! convert to epos id
c if(ist.eq.0)esum=esum+PHKK(4,k)
endif
if(ish.ge.7)write(ifch,'(a,i5,a,i12,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
......@@ -373,8 +465,11 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
pptl(1,nptl)=sngl(PHKK(1,k)) !P_x
pptl(2,nptl)=sngl(PHKK(2,k)) !P_y
pptl(3,nptl)=sngl(PHKK(3,k)) !P_z
pptl(4,nptl)=sngl(PHKK(4,k)) !E
pptl(5,nptl)=sngl(PHKK(5,k)) !mass
c pptl(4,nptl)=sngl(PHKK(4,k)) !E
pptl(5,nptl)=sngl(am) !mass
c put particle on-shell in single precision
pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2
. +pptl(3,nptl)**2+pptl(5,nptl)**2)
ifrptl(1,nptl)=-JDAHKK(1,k) !updated later
ifrptl(2,nptl)=-JDAHKK(2,k) !updated later
xorptl(1,nptl)=sngl(VHKK(1,k))*1e12
......@@ -396,18 +491,19 @@ c fix beam particles which are registered at rest in DPMJET
call utlob5(ytrg,pptl(1,nptl)
. ,pptl(2,nptl),pptl(3,nptl),pptl(4,nptl),pptl(5,nptl)) !boost in cms frame
if(ist.eq.12)then !inelastic nucleon
istptl(nptl)=1
iorptl(nptl)=-1
jorptl(nptl)=0
istptl(nptl)=41
iorptl(nptl)=maproj+1
jorptl(nptl)=maproj+matarg
idptl(nptl)=sign(abs(idptl(nptl))*100+99,idptl(nptl))
ntgevt=ntgevt+1
elseif(ist.eq.14)then !spectator nucleon
istptl(nptl)=1
iorptl(nptl)=0
jorptl(nptl)=0
iorptl(nptl)=maproj+1
jorptl(nptl)=maproj+matarg
elseif(ist.eq.18)then !wounded nucleon
istptl(nptl)=41
iorptl(nptl)=0
jorptl(nptl)=0
iorptl(nptl)=maproj+1
jorptl(nptl)=maproj+matarg
else !ist=16 grey final particles
istptl(nptl)=0
iorptl(nptl)=0 !defined later
......@@ -419,18 +515,19 @@ c fix beam particles which are registered at rest in DPMJET
call utlob5(-yprj,pptl(1,nptl)
. ,pptl(2,nptl),pptl(3,nptl),pptl(4,nptl),pptl(5,nptl)) !boost in cms frame
if(ist.eq.11)then !inelastic nucleon
istptl(nptl)=1
iorptl(nptl)=-1
jorptl(nptl)=0
istptl(nptl)=41
iorptl(nptl)=1
jorptl(nptl)=maproj
idptl(nptl)=sign(abs(idptl(nptl))*100+99,idptl(nptl))
npjevt=npjevt+1
elseif(ist.eq.13)then !spectator nucleon
istptl(nptl)=1
iorptl(nptl)=0
jorptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj
elseif(ist.eq.17)then !wounded nucleon
istptl(nptl)=41
iorptl(nptl)=0
jorptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj
else !ist=15 grey final particles
istptl(nptl)=0
iorptl(nptl)=0 !defined later
......@@ -470,6 +567,7 @@ C $ nptlhep(JMOHKK(1,k))
enddo !! end of loop over particles for copying to epos block
c write(ifch,*)'esum',esum,0.5*engy*(maproj+matarg)
c update daughter index
do i=1,nptl
......@@ -521,18 +619,6 @@ c & 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
c PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
c COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
c & XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
c & XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
c & XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
c & XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
c & XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
c & XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
c & XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
c & XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
c & BSLOPE,NEBINI,NQBINI
* Glauber formalism: flags and parameters for statistics
c LOGICAL LPROD
......@@ -559,7 +645,7 @@ c #######from dt_xstabl##########
end
c------------------------------------------------------------------------------
subroutine GetDPMJETSigmaAA(niter,stotaa,sineaa,selaaa)
subroutine GetDPMJETSigmaAA(niter,stotaa,sineaa,scutaa,selaaa)
c------------------------------------------------------------------------------
c return AA cross sections of DPMJET
c but also fill other cross section variables
......@@ -575,14 +661,11 @@ 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'
integer niter !input
real stotaa,sineaa,selaaa !subroutine return values
real stotaa,sineaa,scutaa,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=100,NEB=8,NQB= 5,KSITEB=50)
COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
......@@ -600,21 +683,98 @@ c xprod seems to be prod + quasi-ela. this is why we subtract it again
LOGICAL LPROD
CHARACTER*8 CGLB
COMMON /DTGLGP/ JSTATB,JBINSB,CGLB,IOGLB,LPROD
JSTATB = niter !loops for glauber calculation !limited by ksiteb = 50
c emulsion treatment
COMMON /DTCOMP/ EMUFRA(NCOMPX),IEMUMA(NCOMPX),IEMUCH(NCOMPX),
& NCOMPO,IEMUL
idum=niter
c IF(NCOMPO.EQ.0)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
c #######from dt_xstabl##########
CALL DT_XSGLAU(IP,IT,IJPROJ,ZERO,ZERO,dble(ECMS),1,1,-1)
c #######from dt_xstabl##########
c CALL DT_GLBSET(IJPROJ,maproj,matarg,dble(elab),2)
c call DT_SIGEMU !to print cross-section in dpmjet
IF(laproj.eq.-1) THEN !not nucleus
IDPDG = idtrafo('nxs','pdg',idproj)
IDP=IDT_ICIHAD(IDPDG)
ELSE
IDP = 1
ENDIF
c EMUFRA(1) = airwnxs(1)
c IF (NCOMPO.GT.0) THEN
c
c IE = 1
c IQ = 1
c XTOT = 0D0
c XELA = 0D0
c XINE = 0D0
c SIGQEP = 0D0
c SIGQET = 0D0
c SIGQE2 = 0D0
c SIGDEL = 0D0
c SIGDQE = 0D0
c IC = 1
c DO 3 ICI=1,NCOMPO
c CALL DT_GLBSET(IDP,maproj,IEMUMA(ICI),dble(elab),1)
c XTOT = XTOT+EMUFRA(ICI)*XSTOT(IE,IQ,IC)
c XELA = XELA+EMUFRA(ICI)*XSELA(IE,IQ,IC)
c XINE = XINE+EMUFRA(ICI)*XSPRO(IE,IQ,IC)
c SIGQEP = SIGQEP+EMUFRA(ICI)*XSQEP(IE,IQ,IC)
c SIGQET = SIGQET+EMUFRA(ICI)*XSQET(IE,IQ,IC)
c SIGQE2 = SIGQE2+EMUFRA(ICI)*XSQE2(IE,IQ,IC)
cc SIGDEL = SIGDEL+EMUFRA(IC)*XSDEL(IE,IQ,IC)
cc SIGDQE = SIGDQE+EMUFRA(IC)*XSDQE(IE,IQ,IC)
c
c 3 CONTINUE
IF (idtargin.eq.0) THEN
IE = 1
IQ = 1
XTOT = 0D0
XELA = 0D0
XINE = 0D0
SIGQEP = 0D0
SIGQET = 0D0
SIGQE2 = 0D0
SIGDEL = 0D0
SIGDQE = 0D0
IC = 1
DO 3 ICI=1,3
CALL DT_GLBSET(IDP,maproj,nint(airanxs(ici)),dble(elab),1)
c CALL DT_XSGLAU(maproj,nint(airanxs(ici)),IDP,ZERO,ZERO,
c . dble(ECMS),1,1,1)
XTOT = XTOT+airwnxs(ici)*XSTOT(IE,IQ,IC)
XELA = XELA+airwnxs(ici)*XSELA(IE,IQ,IC)
XINE = XINE+airwnxs(ici)*XSPRO(IE,IQ,IC)
SIGQEP = SIGQEP+airwnxs(ici)*XSQEP(IE,IQ,IC)
SIGQET = SIGQET+airwnxs(ici)*XSQET(IE,IQ,IC)
SIGQE2 = SIGQE2+airwnxs(ici)*XSQE2(IE,IQ,IC)
c SIGDEL = SIGDEL+airwnxs(ici)*XSDEL(IE,IQ,IC)
c SIGDQE = SIGDQE+airwnxs(ici)*XSDQE(IE,IQ,IC)
3 CONTINUE
else
CALL DT_XSGLAU(maproj,matarg,IDP,ZERO,ZERO,dble(ECMS),1,1,1)
XTOT = XSTOT(1,1,1)
XELA = XSELA(1,1,1)
XINE = XSPRO(1,1,1)
SIGQEP = XSQEP(1,1,1)
SIGQET = XSQET(1,1,1)
SIGQE2 = XSQE2(1,1,1)
endif
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)
XPRO = XTOT - XELA - SIGQEP - SIGQET - SIGQE2
stotaa = sngl(XTOT)
sineaa = sngl(XPRO)
scutaa = sngl(XPRO)
sineaa = sngl(XINE)
selaaa = sngl(XELA)
end
......