IAP GITLAB

Commit 3a26d12e authored by Tanguy Pierog's avatar Tanguy Pierog

reactivate old DPMJET 3.0-6 and secure interface for Phojet, Pythia and DPMJET

parent b1e0b614
......@@ -10,12 +10,13 @@ OPTION (__QGSJET01__ "Build with model" OFF)
OPTION (__GHEISHA__ "Build with model" OFF)
OPTION (__PYTHIA__ "Build with model" OFF)
OPTION (__HIJING__ "Build with model" OFF)
OPTION (__SIBYLL__ "Build with model" ON)
OPTION (__SIBYLL__ "Build with model" OFF)
OPTION (__PHOJET__ "Build with model" OFF)
OPTION (__DPMJET06__ "Build with model" OFF)
OPTION (__DPMJET17__ "Build with model" OFF)
OPTION (__DPMJET19__ "Build with model" ON)
OPTION (__DPMJET19__ "Build with model" OFF)
OPTION (__QGSJETII03__ "Build with model" OFF)
OPTION (__QGSJETII04__ "Build with model" ON)
OPTION (__QGSJETII04__ "Build with model" OFF)
######################################ONLY EDIT THIS######################################
if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
......@@ -156,6 +157,10 @@ FILE(READ ${PROJECT_SOURCE_DIR}/src/dpmjet/3.2019-1/hepevt-HEPMCTEMPLATE templat
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template20 "${template19}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/dpmjet/3.2019-1/DPMJET-19.1/include/pythia/inc/hepevt "${template20}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/dpmjet/3.0-6/pythia6115dpm3v1.f-HEPMCTEMPLATE template21)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template22 "${template21}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/dpmjet/3.0-6/pythia6115dpm3v1.f "${template22}")
......@@ -225,6 +230,15 @@ get_property(HELPER SOURCE src/models.F PROPERTY COMPILE_FLAGS)
set_property(SOURCE src/models.F PROPERTY COMPILE_FLAGS "${HELPER} -D __PHOJET__")
ENDIF (__PHOJET__)
IF (__DPMJET06__)
ADD_SUBDIRECTORY ("${PROJECT_SOURCE_DIR}/src/dpmjet/3.0-6")
LIST(APPEND TABS DPMJET-0.6/dpmjet.dat)
LIST(APPEND TABS DPMJET-0.6/phojet_fitpar.dat)
SET(STATIC_LIBS ${STATIC_LIBS} Dpmjet06)
get_property(HELPER SOURCE src/models.F PROPERTY COMPILE_FLAGS)
set_property(SOURCE src/models.F PROPERTY COMPILE_FLAGS "${HELPER} -D __DPMJET__ -D __DPMJET06__")
ENDIF (__DPMJET06__)
IF (__DPMJET17__)
ADD_SUBDIRECTORY ("${PROJECT_SOURCE_DIR}/src/dpmjet/3.2017-1")
LIST(APPEND TABS DPMJET-17.1/dpmjet.dat)
......
......@@ -15,7 +15,7 @@
#cmakedefine __SIBYLL__
#cmakedefine __QGSJETII03__
#cmakedefine __QGSJETII04__
#cmakedefine __DPMJET__
#cmakedefine __DPMJET06__
#cmakedefine __DPMJET17__
#cmakedefine __DPMJET19__
......@@ -56,6 +56,9 @@ bool CRMCinterface::init(int HEmodel)
#ifdef __QGSJETII03__
case 11: libname << "QgsjetII03"; break;
#endif
#ifdef __DPMJET06__
case 12: libname << "Dpmjet06"; break;
#endif
#ifdef __DPMJET17__
case 12: libname << "Dpmjet17"; break;
#endif
......
......@@ -94,7 +94,7 @@ CRMCoptions::ParseOptions(int argc, char** argv)
<< ", 3=Gheisha"
#endif
#ifdef __PYTHIA__
<< ", 4=Pythia_6.115"
<< ", 4=Pythia_6.4.28"
#endif
#ifdef __HIJING__
<< ", 5=Hijing_1.38"
......@@ -111,6 +111,9 @@ CRMCoptions::ParseOptions(int argc, char** argv)
#ifdef __QGSJETII03__
<< ", 11=QGSJETII-03"
#endif
#ifdef __DPMJET06__
<< ", 12=DPMJet 3.0-6"
#endif
#ifdef __DPMJET17__
<< ", 12=DPMJet-III_2017.1"
#endif
......@@ -392,6 +395,9 @@ CRMCoptions::DumpConfig() const
case 7: cout << " (QGSJETII-04) \n"; break;
case 8: cout << " (Phojet) \n"; break;
case 11: cout << " (QGSJETII-03) \n"; break;
#ifdef __DPMJET06__
case 12: cout << " (DPMJet 3.0-6) \n"; break;
#endif
#ifdef __DPMJET17__
case 12: cout << " (DPMJet-III 2017.1) \n"; break;
#endif
......@@ -479,11 +485,14 @@ CRMCoptions::GetOutputFileName() const
case 7: crmcFileName << "qgsjetII04"; break;
case 8: crmcFileName << "phojet"; break;
case 11: crmcFileName << "qgsjetII03"; break;
#ifdef __DPMJET06__
case 12: crmcFileName << "dpmjet3.0-6"; break;
#endif
#ifdef __DPMJET17__
case 12: crmcFileName << "dpmjetIII17"; break;
case 12: crmcFileName << "dpmjetIII.17"; break;
#endif
#ifdef __DPMJET19__
case 12: crmcFileName << "dpmjetIII19"; break;
case 12: crmcFileName << "dpmjetIII.19"; break;
#endif
default:
cerr << " crmcOut: error - unknown model " << fHEModel << endl;
......
......@@ -34,6 +34,8 @@ nodecay -19 !uncomment not to decay antialpha
MinDecayLength 1. !minimum c.Tau to define stable particles (cm)
fdpmjet path @CRMCROOT@/tabs/
fdpmjet dat @CRMCROOT@/tabs/dpmjet.dat
fdpmjet pho @CRMCROOT@/tabs/phojet_fitpar.dat
fqgsjet dat @CRMCROOT@/tabs/qgsjet.dat
fqgsjet ncs @CRMCROOT@/tabs/qgsjet.ncs
fqgsjetII03 dat @CRMCROOT@/tabs/qgsdat-II-03.lzma
......
......@@ -3,13 +3,13 @@ FILE(GLOB files *.f)
set_source_files_properties( dpmjet3.0-6.f dpmjet-crmc.f phojet1.12-35c4.f pythia6115dpm3v1.f ../../crmc-aaa.f ../../models.F PROPERTIES COMPILE_FLAGS "-D __DPMJET__")
IF (__CRMCSTATIC__)
add_library(Dpmjet STATIC dpmjet3.0-6.f dpmjet-crmc.f phojet1.12-35c4.f pythia6115dpm3v1.f ../../crmc-aaa.f ../../models.F)
add_library(Dpmjet06 STATIC dpmjet3.0-6.f dpmjet-crmc.f phojet1.12-35c4.f pythia6115dpm3v1.f ../../crmc-aaa.f ../../models.F)
ELSE (__CRMCSTATIC__)
add_library(Dpmjet SHARED dpmjet3.0-6.f dpmjet-crmc.f phojet1.12-35c4.f pythia6115dpm3v1.f ../../crmc-aaa.f ../../models.F)
target_link_libraries(Dpmjet CrmcBasic)
add_library(Dpmjet06 SHARED dpmjet3.0-6.f dpmjet-crmc.f phojet1.12-35c4.f pythia6115dpm3v1.f ../../crmc-aaa.f ../../models.F)
target_link_libraries(Dpmjet06 CrmcBasic)
ENDIF (__CRMCSTATIC__)
INSTALL (TARGETS Dpmjet
INSTALL (TARGETS Dpmjet06
RUNTIME DESTINATION bin
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib/static
......
......@@ -88,13 +88,14 @@ c particle without breit-wigner decay will be decayed in EPOS to get full histor
idpdg = KCHG(I,4)
c don't touch decays with unknown code
eposid = idtrafostatus('pdg','nxs',idpdg,istatus)
idepos = idtrafostatus('pdg','nxs',idpdg,istatus)
if(abs(idpdg).eq.10323)print *,idpdg,idepos,istatus
if(istatus.ne.0) then
go to 100
endif
c don't touch decays with unknown width
call idtaustatus(eposid,0,0,taugm,istatus)
call idtaustatus(idepos,0,0,taugm,istatus)
c print *,i,(taugm.eq.indexOutOfRange.or.taugm.eq.iijjOutOfRange)
c + ,(PMAS(i,2).LT.PARP(41))
if(istatus.ne.0) then
......@@ -230,7 +231,8 @@ c PLEASE FILL!!
GOTO 1001
ENDIF
do k=1,NHKK
do 100 k=1,NHKK
nptlhep(k)=0
c LIST is the code of final particle, P - its 4-momentum and mass.
ic=IDHKK(k) !! copies IDs
......@@ -240,8 +242,6 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
$ , ' momentum :',ISTHKK(k),(sngl(PHKK(i,k)),i=1,5)
$ ,JMOHKK(1,k),JMOHKK(2,k),JDAHKK(1,k),JDAHKK(2,k)
nptl=nptl+1 !! add 1 particle to stack
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...
......@@ -250,8 +250,23 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
else
ist=abs(ISTHKK(k)) !!0 means last generation other codes are e.g
endif
id = idtrafostatus('pdg','nxs',ic,istatus)
if(istatus.ne.0)then
if(ist.eq.0)then
print *,ic,id,istatus
stop "pdg id from DPMJet unknown !"
else
goto 100
endif
endif
id=idtrafo('pdg','nxs',ic) !! convert to epos ids
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
nptlhep(k)=nptl
if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
......@@ -351,7 +366,7 @@ C $ nptlhep(JMOHKK(1,k))
$ , ' momentum :',(pptl(i,nptl),i=1,5)
enddo !! end of loop over particles for copying to epos block
100 continue !! end of loop over particles for copying to epos block
c update daughter index
do i=1,nptl
......
......@@ -168,7 +168,7 @@ ctp modifications for CRMC line 16784 and in DT_INIT for NCASES=-1
* flags for activated histograms
COMMON /DTHIS3/ IHISPP(50),IHISXS(50),IXSTBL
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
* LEPTO
**LUND single / double precision
REAL CUT,PARL,TMPX,TMPY,TMPW2,TMPQ2,TMPU
......@@ -2116,13 +2116,13 @@ C IF (IDP.EQ.27) IDP = 6
ENDIF
* initialization of Glauber-formalism (moved to xAEVT, sr 26.3.96)
IF (NCOMPO.LE.0) THEN
CALL DT_SHMAKI(IP,IPZ,IT,ITZ,IDP,PPN,IGLAU)
ELSE
DO 493 I=1,NCOMPO
CALL DT_SHMAKI(IP,IPZ,IEMUMA(I),IEMUCH(I),IDP,PPN,0)
493 CONTINUE
ENDIF
c IF (NCOMPO.LE.0) THEN
c CALL DT_SHMAKI(IP,IPZ,IT,ITZ,IDP,PPN,IGLAU)
c ELSE
c DO 493 I=1,NCOMPO
c CALL DT_SHMAKI(IP,IPZ,IEMUMA(I),IEMUCH(I),IDP,PPN,0)
c 493 CONTINUE
c ENDIF
* pre-tabulation of elastic cross-sections
CALL DT_SIGTBL(JDUM,JDUM,DUM,DUM,-1)
......@@ -14997,7 +14997,7 @@ C DO 1 I=NPOINT(5),NEND
COMMON /DTEVT2/ IDRES(NMXHKK),IDXRES(NMXHKK),NOBAM(NMXHKK),
& IDBAM(NMXHKK),IDCH(NMXHKK),NPOINT(10),
& IHIST(2,NMXHKK)
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
PARAMETER (MAXLND=4000)
COMMON/PYJETS/N,NPAD,K(MAXLND,5),P(MAXLND,5),V(MAXLND,5)
......@@ -15586,7 +15586,7 @@ C ISU = 4
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
* flags for particle decays
COMMON /DTFRPA/ MSTUX(20),PARUX(20),MSTJX(20),PARJX(20),
& IMSTU(20),IPARU(20),IMSTJ(20),IPARJ(20),
......@@ -20296,8 +20296,8 @@ C & ZSTLIN (2,IA) = MAX ( SQRT (ZSTLIN(2,IA)-ZSTLIN(1,IA)**2),
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
CHARACTER*500 FNDPMJET
COMMON/DPMJETFNAME/ FNDPMJET
CHARACTER*500 fndpmjet,fndpmjetpho
COMMON/DPMJETFNAME/ fndpmjet,fndpmjetpho
PARAMETER ( CSNNRM = 2.0D-15 )
PARAMETER ( ZERZER = 0.D+00 )
......@@ -28339,56 +28339,56 @@ C & '-------------------------------'
*
*===pohist=============================================================*
*
SUBROUTINE PHO_PHIST(IMODE,WEIGHT)
IMPLICIT DOUBLE PRECISION (A-H,O-X,Z)
SAVE
COMMON /DTIONT/ LINP,LOUT,LDAT
PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
* Glauber formalism: cross sections
COMMON /DTGLXS/ ECMNN(NEB),Q2G(NQB),ECMNOW,Q2,
& XSTOT(NEB,NQB,NCOMPX),XSELA(NEB,NQB,NCOMPX),
& XSQEP(NEB,NQB,NCOMPX),XSQET(NEB,NQB,NCOMPX),
& XSQE2(NEB,NQB,NCOMPX),XSPRO(NEB,NQB,NCOMPX),
& XSDEL(NEB,NQB,NCOMPX),XSDQE(NEB,NQB,NCOMPX),
& XETOT(NEB,NQB,NCOMPX),XEELA(NEB,NQB,NCOMPX),
& XEQEP(NEB,NQB,NCOMPX),XEQET(NEB,NQB,NCOMPX),
& XEQE2(NEB,NQB,NCOMPX),XEPRO(NEB,NQB,NCOMPX),
& XEDEL(NEB,NQB,NCOMPX),XEDQE(NEB,NQB,NCOMPX),
& BSLOPE,NEBINI,NQBINI
ILAB = 0
IF (IMODE.EQ.10) THEN
IMODE = 1
ILAB = 1
ENDIF
IF (ABS(IMODE).LT.1000) THEN
* PHOJET-statistics
C CALL POHISX(IMODE,WEIGHT)
IF (IMODE.EQ.-1) THEN
MODE = 1
XSTOT(1,1,1) = WEIGHT
ENDIF
IF (IMODE.EQ. 1) MODE = 2
IF (IMODE.EQ.-2) MODE = 3
IF (MODE.EQ.2) CALL DT_SWPPHO(ILAB)
C IF (MODE.EQ.3) WRITE(LOUT,*)
C & ' Sigma = ',XSPRO(1,1,1),' mb used for normalization'
CALL DT_HISTOG(MODE)
CALL DT_USRHIS(MODE)
ELSE
* DTUNUC-statistics
MODE = IMODE/1000
C IF (MODE.EQ.3) WRITE(LOUT,*)
C & ' Sigma = ',XSPRO(1,1,1),' mb used for normalization'
CALL DT_HISTOG(MODE)
CALL DT_USRHIS(MODE)
ENDIF
RETURN
END
c SUBROUTINE PHO_PHIST(IMODE,WEIGHT)
c
c IMPLICIT DOUBLE PRECISION (A-H,O-X,Z)
c SAVE
c
c COMMON /DTIONT/ LINP,LOUT,LDAT
c PARAMETER (NCOMPX=20,NEB=8,NQB= 5,KSITEB=50)
c* Glauber formalism: cross sections
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
c
c ILAB = 0
c IF (IMODE.EQ.10) THEN
c IMODE = 1
c ILAB = 1
c ENDIF
c IF (ABS(IMODE).LT.1000) THEN
c* PHOJET-statistics
cC CALL POHISX(IMODE,WEIGHT)
c IF (IMODE.EQ.-1) THEN
c MODE = 1
c XSTOT(1,1,1) = WEIGHT
c ENDIF
c IF (IMODE.EQ. 1) MODE = 2
c IF (IMODE.EQ.-2) MODE = 3
c IF (MODE.EQ.2) CALL DT_SWPPHO(ILAB)
cC IF (MODE.EQ.3) WRITE(LOUT,*)
cC & ' Sigma = ',XSPRO(1,1,1),' mb used for normalization'
c CALL DT_HISTOG(MODE)
c CALL DT_USRHIS(MODE)
c ELSE
c* DTUNUC-statistics
c MODE = IMODE/1000
cC IF (MODE.EQ.3) WRITE(LOUT,*)
cC & ' Sigma = ',XSPRO(1,1,1),' mb used for normalization'
c CALL DT_HISTOG(MODE)
c CALL DT_USRHIS(MODE)
c ENDIF
c
c RETURN
c END
c
*$ CREATE DT_SWPPHO.FOR
*COPY DT_SWPPHO
*
......@@ -33002,15 +33002,15 @@ C SUBROUTINE POHISX(I,X)
*COPY PHO_LHIST
*
*===poluhi=============================================================*
*
SUBROUTINE PHO_LHIST(I,X)
**
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
RETURN
END
c*
c SUBROUTINE PHO_LHIST(I,X)
c**
c
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c SAVE
c
c RETURN
c END
*$ CREATE PDFSET.FOR
*COPY PDFSET
......@@ -509,7 +509,7 @@ C hard cross sections and MC selection weights
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER PYCOMP
......@@ -34813,7 +34813,7 @@ C input/output channels
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
INTEGER MDCY,MDME,KFDP
DOUBLE PRECISION BRAT
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYDAT3/MDCY(500,3),MDME(8000,2),BRAT(8000),KFDP(8000,5)
INTEGER PYCOMP
......@@ -223,13 +223,13 @@ c particle without breit-wigner decay will be decayed in EPOS to get full histor
idpdgi = KCHG(I,4)
c don't touch decays with unknown code
eposid = idtrafostatus('pdg','nxs',idpdgi,istatus)
idepos = idtrafostatus('pdg','nxs',idpdgi,istatus)
if(istatus.ne.0) then
go to 100
endif
c don't touch decays with unknown width
call idtaustatus(eposid,0,0,taugm,istatus)
call idtaustatus(idepos,0,0,taugm,istatus)
c print *,i,(taugm.eq.indexOutOfRange.or.taugm.eq.iijjOutOfRange)
c + ,(PMAS(i,2).LT.PARP(41))
if(istatus.ne.0) then
......@@ -415,7 +415,8 @@ c PLEASE FILL!!
anintine=anintine+1.
c esum=0.d0
do k=1,NHKK
do 100 k=1,NHKK
nptlhep(k)=0
c LIST is the code of final particle, P - its 4-momentum and mass.
ic=IDHKK(k) !! copies IDs
......@@ -425,17 +426,7 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
$ , ' momentum :',ISTHKK(k),(sngl(PHKK(i,k)),i=1,5)
$ ,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...
......@@ -460,9 +451,25 @@ c esum=esum+PHKK(4,k)
elseif(ic.eq.66666)then !cascade
id=80000000
else
id=idtrafo('pdg','nxs',ic) !! convert to epos id
id = idtrafostatus('pdg','nxs',ic,istatus)
if(istatus.ne.0)then
if(ist.eq.0)then
print *,ic,id,istatus
stop "pdg id from DPMJet unknown !"
else
goto 100
endif
endif
c if(ist.eq.0)esum=esum+PHKK(4,k)
endif
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
nptlhep(k)=nptl
if(ish.ge.7)write(ifch,'(a,i5,a,i12,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
......@@ -571,7 +578,7 @@ C $ nptlhep(JMOHKK(1,k))
$ , ' momentum :',(pptl(i,nptl),i=1,5)
enddo !! end of loop over particles for copying to epos block
100 continue !! end of loop over particles for copying to epos block
c write(ifch,*)'esum',esum,0.5*engy*(maproj+matarg)
c update daughter index
......
......@@ -223,13 +223,13 @@ c particle without breit-wigner decay will be decayed in EPOS to get full histor
idpdgi = KCHG(I,4)
c don't touch decays with unknown code
eposid = idtrafostatus('pdg','nxs',idpdgi,istatus)
idepos = idtrafostatus('pdg','nxs',idpdgi,istatus)
if(istatus.ne.0) then
go to 100
endif
c don't touch decays with unknown width
call idtaustatus(eposid,0,0,taugm,istatus)
call idtaustatus(idepos,0,0,taugm,istatus)
c print *,i,(taugm.eq.indexOutOfRange.or.taugm.eq.iijjOutOfRange)
c + ,(PMAS(i,2).LT.PARP(41))
if(istatus.ne.0) then
......@@ -415,7 +415,8 @@ c PLEASE FILL!!
anintine=anintine+1.
c esum=0.d0
do k=1,NHKK
do 100 k=1,NHKK
nptlhep(k)=0
c LIST is the code of final particle, P - its 4-momentum and mass.
ic=IDHKK(k) !! copies IDs
......@@ -425,17 +426,7 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
$ , ' momentum :',ISTHKK(k),(sngl(PHKK(i,k)),i=1,5)
$ ,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...
......@@ -460,9 +451,25 @@ c esum=esum+PHKK(4,k)
elseif(ic.eq.66666)then !cascade
id=80000000
else
id=idtrafo('pdg','nxs',ic) !! convert to epos id
id = idtrafostatus('pdg','nxs',ic,istatus)
if(istatus.ne.0)then
if(ist.eq.0)then
print *,ic,id,istatus
stop "pdg id from DPMJet unknown !"
else
goto 100
endif
endif
c if(ist.eq.0)esum=esum+PHKK(4,k)
endif
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
nptlhep(k)=nptl
if(ish.ge.7)write(ifch,'(a,i5,a,i12,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
......@@ -571,7 +578,7 @@ C $ nptlhep(JMOHKK(1,k))
$ , ' momentum :',(pptl(i,nptl),i=1,5)
enddo !! end of loop over particles for copying to epos block
100 continue !! end of loop over particles for copying to epos block
c write(ifch,*)'esum',esum,0.5*engy*(maproj+matarg)
c update daughter index
......
......@@ -12,8 +12,8 @@ c-----------------------------------------------------------------------
common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi/cifset/ifset
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
c &,fndpmjet,fndpmjetpho,fndpmpath !qgs-II????????
c common/dpmjetfname/ fndpmjet,fndpmjetpho,fndpmpath
&,fndpmjet,fndpmjetpho!,fndpmpath !qgs-II????????
common/dpmjetfname/ fndpmjet,fndpmjetpho!,fndpmpath
common/qgsfname/ fndat, fnncs, ifdat, ifncs
common/qgsIIfname/fnIIdat, fnIIncs, ifIIdat, ifIIncs !qgs-II????????
common/qgsII03fname/fnII03dat, fnII03ncs, ifII03dat, ifII03ncs !qgs-II????????
......@@ -3498,8 +3498,8 @@ c---------------------------------------------------------------------
common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi
common /ems12/bidiba,amhdibar,iodiba ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
c &,fndpmjet,fndpmjetpho,fndpmpath
c common/dpmjetfname/ fndpmjet,fndpmjetpho,fndpmpath
&,fndpmjet,fndpmjetpho!,fndpmpath
common/dpmjetfname/ fndpmjet,fndpmjetpho!,fndpmpath
cdh datadir for path to the data sets to be read in by dpmjet/phojet
COMMON /DATADIR/ DATADIR
CHARACTER*132 DATADIR
......@@ -3765,8 +3765,8 @@ c if(line(i:j).eq.'xEmsPx')call xEmsPxNo(2,0.,0.,0,0)
call utworn(line,j,ne)
if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-name?'
call utword(line,i,j,0)
c if(linex(ix:jx).eq.'dat')fndpmjet(1:j-i+2)=line(i:j)//' '
c if(linex(ix:jx).eq.'pho')fndpmjetpho(1:j-i+2)=line(i:j)//' '
if(linex(ix:jx).eq.'dat')fndpmjet(1:j-i+1)=line(i:j)
if(linex(ix:jx).eq.'pho')fndpmjetpho(1:j-i+1)=line(i:j)
if(linex(ix:jx).eq.'path')DATADIR(1:j-i+2)=line(i:j)//' '
elseif(line(i:j).eq.'fqgsjet')then !QGSJet
......
......@@ -856,6 +856,9 @@ C...Format for error printout.
RETURN
END
#endif
#if __PYTHIA__ || __DPMJET__
c--------------------------------------------------------------------
subroutine PHO_PHIST(idum,dum)
......
......@@ -213,6 +213,7 @@ c set projectile and target as non final particles
if(ish.ge.5)write(ifch,'(a,i5)')
$ ' number of particles from PHOJET :',NHEPP
do 500 k=1,NHEPP
nptlhep(k)=0
c LLIST is the code of final particle, P - its 4-momentum and mass.
ic=IDHEPP(k)
......@@ -221,20 +222,31 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
$ ' PHOJET particle ',k,' id :',ic,' before conversion'
$ , ' momentum :',ISTHEPP(k),(sngl(PPHEP(i,k)),i=1,5)
IF(ISTHEPP(k).GE.1.AND.ISTHEPP(k).LE.2)THEN
id = idtrafostatus('pdg','nxs',ic,istatus)
if(istatus.ne.0)then
if(ISTHEPP(k).eq.1)then
print *,ic,id,istatus
stop "pdg id from Phojet unknown !"
else
if(ish.ge.7)write(ifch,'(a)')
& 'Not final particle with unknown id: skipped ...'
goto 500
endif
endif
c IF(ISTHEPP(k).GE.1.AND.ISTHEPP(k).LE.2)THEN
IF(ISTHEPP(k).EQ.1.OR.k.LE.2)THEN
c final particle
nptl=nptl+1
nptlhep(k)=nptl
istptl(nptl)=ISTHEPP(k)-1
else