IAP GITLAB

Commit 27994e23 authored by Tanguy Pierog's avatar Tanguy Pierog

fix problem with hepevt common block size for various models

fix problems with all hadronic interaction models to have proper consistent final output
fix problem in DPMJET to get nuclear spectator (free by default but can form nuclear fragment if infragm is set to 2)
fix decay definition to use CRMC settings for all PYTHIA based model (in Phojet decay is done in EPOS otherwise no mother information available)
Note : Hijing may work with g77 but doesn't run with gfortran (endless loop)
Hopefully now all models can be used properly.


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4190 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 4a15e47d
......@@ -4,17 +4,17 @@ PROJECT (crmc)
######################################ONLY EDIT THIS######################################
# Enable/Disable models to be built
OPTION (__CRMCSTATIC__ "Build with static library" OFF)
OPTION (__CRMCSTATIC__ "Build with static library" OFF) #if ON should not combined DPMJET/PHOJET/PYTHIA because they use different version of pythia (for dynamic library no problem)
OPTION (__QGSJET01__ "Build with model" ON)
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 (__PHOJET__ "Build with model" OFF)
OPTION (__DPMJET__ "Build with model" OFF) #experimental
OPTION (__QGSJETII03__ "Build with model" ON)
OPTION (__QGSJETII04__ "Build with model" ON)
OPTION (__PYTHIA__ "Build with model" ON)
OPTION (__HIJING__ "Build with model" OFF) #doesn't work with gfortran !
OPTION (__SIBYLL__ "Build with model" OFF)
OPTION (__PHOJET__ "Build with model" ON)
OPTION (__DPMJET__ "Build with model" ON) #experimental
OPTION (__QGSJETII03__ "Build with model" OFF)
OPTION (__QGSJETII04__ "Build with model" OFF)
######################################ONLY EDIT THIS######################################
if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
......@@ -86,6 +86,28 @@ FILE(READ ${PROJECT_SOURCE_DIR}/src/CRMCinterface.h-HEPMC-TEMPLATE template3)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template4 "${template3}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/CRMCinterface.h "${template4}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/pythia/pythia6115.f-HEPMC-TEMPLATE template5)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template6 "${template5}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/pythia/pythia6115.f "${template6}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/pythia/pythia6215.f-HEPMC-TEMPLATE template7)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template8 "${template7}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/pythia/pythia6215.f "${template8}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/pythia/pythia-6.4.28.f-HEPMC-TEMPLATE template9)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template10 "${template9}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/pythia/pythia-6.4.28.f "${template10}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/hijing/hipyset1.35.f-HEPMC-TEMPLATE template11)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template12 "${template11}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/hijing/hipyset1.35.f "${template12}")
FILE(READ ${PROJECT_SOURCE_DIR}/src/dpmjet/3.0-6/pythia6115dpm3v1.f-HEPMCTEMPLATE template13)
STRING( REGEX REPLACE "HEPEVT_SIZE_REPLACE" "${Replace_String}" template14 "${template13}")
FILE(WRITE ${PROJECT_SOURCE_DIR}/src/dpmjet/3.0-6/pythia6115dpm3v1.f "${template14}")
## configure a header file to pass some of the CMake settings to the source code
CONFIGURE_FILE (
......
......@@ -150,7 +150,7 @@ c Fix final particles and some event parameters
call afinal
c Fill HEP common
if(model.ne.4.and.model.ne.5)call hepmcstore(iout) !use hepmcstore for non-HEP compatible events or non h-p models (Nucleus as projectile or target).
call hepmcstore(iout) !use hepmcstore for all models to be sure to get same vertex structure
c optional Statistic information (only with debug level ish=1)
call astati
......@@ -321,7 +321,7 @@ c nody(nrnody)=idtrafo('pdg','nxs',3122) !lambda using pdg code
c Debug
ish=0 !debug level
ifch=6 !debug output (screen)
c ifch=6 !debug output (screen)
c ifch=31 !debug output (file)
c fnch="epos.debug"
c nfnch=index(fnch,' ')-1
......@@ -334,6 +334,7 @@ c open(ifcx,file=fnch(1:nfnch),status='unknown')
c pnll=pproj !beam momentum GeV/c
infragm=2 !nuclear fragmentation (realistic)
if(model.eq.12)infragm=0 !no EPOS fragmentation by default but can be switch on !(evaporation in DPMJET switch off ???)
c Projecticle definitions
if (abs(ipart) .eq. 1) then
......
......@@ -13,6 +13,8 @@
nodecay 14 !uncomment not to decay mu- (PDG id = 13)
nodecay -14 !uncomment not to decay mu+ (PDG id = -13)
nodecay 1120 !uncomment not to decay proton (PDG id = 2212) (for pythia)
nodecay -1120 !uncomment not to decay aproton (PDG id = -2212) (for pythia)
nodecay 1220 !uncomment not to decay neutron (PDG id = 2112)
nodecay -1220 !uncomment not to decay aneutron (PDG id = -2112)
nodecay 120 !uncomment not to decay pi+ (PDG id = 211)
......@@ -46,6 +48,8 @@ fname inirj @CRMCROOT@/tabs/epos.inirj
fname inics @CRMCROOT@/tabs/epos.inics
fname inihy @CRMCROOT@/tabs/epos.inihy
set pytune 350 !possibility to change PYTHIA tune (for PYTHIA only !)
!!ImpactParameter
!set bminim 0 !works with epos
!set bmaxim 4
......
......@@ -13,90 +13,6 @@ c common block var iframe renamed to idframe since it is used by epos
COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
& LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IDFRAME,ITRSPT
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
character*160 NameDecay
c Integer function to translate normal id from pythia to compact one.
integer PYCOMP
if(idecay.eq.1.and.ctaumin.le.0.)then
if(mod(ndecay/10,10).eq.1)then !Kshort, Klong
ic1=PYCOMP(130)
ic2=PYCOMP(310)
write(NameDecay,'(a,i3,a,i3,a)')
& 'MDCY(',ic1,',1)=0;MDCY(',ic2,',1)=0'
call PYGIVE(NameDecay(1:27))
endif
if(mod(ndecay/100,10).eq.1)then !Lambda
ic=PYCOMP(3122)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/1000,10).eq.1)then !sigma+
ic=PYCOMP(3222)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/1000,10).eq.1)then !sigma-
ic=PYCOMP(3112)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/10000,10).eq.1)then !Xi+/-
ic=PYCOMP(3312)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/10000,10).eq.1)then !Xi0
ic=PYCOMP(3322)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/100000 ,10).eq.1)then !omega
ic=PYCOMP(3334)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(mod(ndecay/1000000,10).eq.1)then !pi0
ic=PYCOMP(111)
write(NameDecay,'(a,i3,a)')
& 'MDCY(',ic,',1)=0'
call PYGIVE(NameDecay(1:13))
endif
if(nrnody.gt.0)then !all other particle
do nod=1,nrnody
id=idtrafo('nxs','pdg',nody(nod))
ic=PYCOMP(id)
write(NameDecay,'(a,i3,a)')'MDCY(',ic,',1)=0)!'
call PYGIVE(NameDecay(1:13))
enddo
endif
elseif(ctaumin.gt.0.)then
MSTJ(21) = 1 !switch on decay according to life time
MSTJ(22) = 2
PARJ(71) = ctaumin*10. !cm->mm
else
MSTJ(21) = 0 !switch off all decays
endif
C switch on decays from pythia
MSTJ(21) = 1
MSTU(22) = 2 ! lifetime smaller than parj(72)
PARJ(71) = ctaumin*10. ! in mm
IDFRAME=2 !dpmjet iframe variable. nucleon-nucleon frame
......@@ -106,6 +22,7 @@ C switch on decays from pythia
irescl=0 !don't rescale/skip events with wrong energy
END
c-----------------------------------------------------------------------
subroutine IniEvtDpm
c-----------------------------------------------------------------------
......@@ -143,6 +60,8 @@ c will be treated as nucleus in DT_INIT and NPMASS,... will be used
CALL DT_INIT(NCASES,EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
+,IGLAU)
c set decay flag in Pythia for DPMJET (after DT_INIT otherwise default is used)
call IniDkyJetset
c initialize cross-sections by calling epos-bas.f function -> models.F -> dpmjet-crmc.f
call xsigma
......@@ -185,8 +104,8 @@ c-----------------------------------------------------------------------
anintine=anintine+1.
call conre !! adds beam particles that are known from intial call of crmc
call conwr
c call conre !! adds beam particles that are known from intial call of crmc
c call conwr
IREJ = 0
KKMAT = -1
......@@ -239,8 +158,8 @@ c PLEASE FILL!!
ncol=1
nevt=1
kolevt=ncol
npjevt=maproj
ntgevt=matarg
npjevt=0
ntgevt=0
pmxevt=pnll
egyevt=ecms
bimevt=BIMPAC
......@@ -259,18 +178,21 @@ c PLEASE FILL!!
c LIST is the code of final particle, P - its 4-momentum and mass.
ic=IDHKK(k) !! copies IDs
if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,4(e10.4,1x))')
if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,i3,1x,5(e10.4,1x),4i4)')
$ ' DPMJET particle ',k,' id :',ic,' before conversion'
$ , ' momentum :',(sngl(PHKK(i,k)),i=1,5)
IF(ISTHKK(k).GE.1 .AND. ISTHKK(k).LE.2)THEN !! if final particle
nptl=nptl+1 !! add 1 particle to stack
nptlhep(k)=nptl
istptl(nptl)=ISTHKK(k)-1 !!0 means last generation other codes are e.g. for pomerons, remnants...
$ , ' 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...
elseif(ic.eq.99999)then
ist=31
else
nptlhep(k)=0
cycle
endif !! all non final particles skipped
ist=abs(ISTHKK(k)) !!0 means last generation other codes are e.g
endif
id=idtrafo('pdg','nxs',ic) !! convert to epos ids
if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
......@@ -282,28 +204,82 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
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
ityptl(nptl)=0
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
ifrptl(1,nptl)=-JDAHKK(1,k) !updated later
ifrptl(2,nptl)=-JDAHKK(2,k) !updated later
xorptl(1,nptl)=sngl(VHKK(1,k))*1e12
xorptl(2,nptl)=sngl(VHKK(2,k))*1e12
xorptl(3,nptl)=sngl(VHKK(3,k))*1e12
xorptl(4,nptl)=sngl(VHKK(4,k))*1e12
tivptl(1,nptl)=0.
tivptl(2,nptl)=0.
ityptl(nptl)=0
idptl(nptl)=id
c treatment of mothers and daughters
iorptl(nptl)=1 !by default set to first and last beam particle
jorptl(nptl)=maproj+matarg
if(ist.ge.11.and.ist.le.18)then
c fix beam particles which are registered at rest in DPMJET
if(mod(ist,2).eq.0)then !target remnant (ist=12,14,16,18)
call utlob5(yhaha,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
ntgevt=ntgevt+1
elseif(ist.eq.14)then !spectator nucleon
istptl(nptl)=0
iorptl(nptl)=0
jorptl(nptl)=0
elseif(ist.eq.18)then !wounded nucleon
istptl(nptl)=41
iorptl(nptl)=0
jorptl(nptl)=0
else !ist=16 grey final particles
istptl(nptl)=0
iorptl(nptl)=0 !defined later
jorptl(nptl)=0 !defined later
endif
else !projectile remnant (ist=11,13,15,17)
call utlob5(-yhaha,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
npjevt=npjevt+1
elseif(ist.eq.13)then !spectator nucleon
istptl(nptl)=0
iorptl(nptl)=0
jorptl(nptl)=0
elseif(ist.eq.17)then !wounded nucleon
istptl(nptl)=41
iorptl(nptl)=0
jorptl(nptl)=0
else !ist=15 grey final particles
istptl(nptl)=0
iorptl(nptl)=0 !defined later
jorptl(nptl)=0 !defined later
endif
endif
else
istptl(nptl)=ist
iorptl(nptl)=0
jorptl(nptl)=0
endif
if(JMOHKK(1,k).ne.0.AND.nptlhep(JMOHKK(1,k)).NE.0) then !mother
iorptl(nptl)=nptlhep(JMOHKK(1,k))
endif
endif
if(JMOHKK(2,k).ne.0.AND.nptlhep(JMOHKK(2,k)).NE.0) then !father
jorptl(nptl)=nptlhep(JMOHKK(2,k))
endif
if(jorptl(nptl).gt.0.and.iorptl(nptl).gt.jorptl(nptl))then !reorder mother/father to avoid problem with hepmc
io=iorptl(nptl)
iorptl(nptl)=jorptl(nptl)
jorptl(nptl)=io
if(io.ne.0)istptl(io)=istptl(io)+2
endif
C print*,istmax,k,JMOHKK(1,k),nptlhep(k),iorptl(nptl),
C $ nptlhep(JMOHKK(1,k))
......@@ -316,6 +292,16 @@ C $ nptlhep(JMOHKK(1,k))
enddo !! end of loop over particles for copying to epos block
c update daughter index
do i=1,nptl
if(ifrptl(1,i).lt.0)then
ifrptl(1,i)=nptlhep(-ifrptl(1,i))
endif
if(ifrptl(2,i).lt.0)then
ifrptl(2,i)=nptlhep(-ifrptl(2,i))
endif
enddo
iret=0
1000 return
......
ctp modifications for CRMC line 16784 and in DT_INIT for NCASES=-1
*$ CREATE DT_INIT.FOR
*COPY DT_INIT
*
......@@ -245,6 +246,7 @@ C DIMENSION XPARA(5)
* bypass reading input cards (e.g. for use with Fluka)
* in this case Epn is expected to carry the beam momentum
IF (NCASES.EQ.-1) THEN
CALL DT_INITJS(0) !call it here first time to define correctly decay from outside
PPN = EPNSAV
EPN = ZERO
CMENER = ZERO
......@@ -16779,7 +16781,7 @@ C ENDIF
& 2,IDUM,IDUM)
* modify entry for interacting nucleons
IF (ISTHKK(IS).EQ.12+ICAS) ISTHKK(IS)=16+ICAS
IF (ISTHKK(IS).EQ.14+ICAS) ISTHKK(IS)=2
IF (ISTHKK(IS).EQ.14+ICAS) ISTHKK(IS)=16+ICAS !2 !ctp20131212 for crmc ????
IF (I.GE.2) THEN
JDAHKK(1,IS) = JDAHKK(1,IDXSPE(1))
JDAHKK(2,IS) = JDAHKK(2,IDXSPE(1))
......@@ -27313,16 +27315,16 @@ C WRITE(6,1000)
*
*===pyr================================================================*
*
DOUBLE PRECISION FUNCTION PYR(IDUMMY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DUMMY = DBLE(IDUMMY)
PYR = DT_RNDM(DUMMY)
RETURN
END
c DOUBLE PRECISION FUNCTION PYR(IDUMMY)
c
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c SAVE
c
c DUMMY = DBLE(IDUMMY)
c PYR = DT_RNDM(DUMMY)
c
c RETURN
c END
*$ CREATE DT_TITLE.FOR
*COPY DT_TITLE
......@@ -33018,26 +33020,26 @@ C
C dummy subroutines, remove to link PDFLIB
C
C**********************************************************************
SUBROUTINE PDFSET(PARAM,VALUE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION PARAM(20),VALUE(20)
CHARACTER*20 PARAM
END
*$ CREATE STRUCTM.FOR
*COPY STRUCTM
*
SUBROUTINE STRUCTM(XI,SCALE2,UV,DV,US,DS,SS,CS,BS,TS,GL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
END
*$ CREATE STRUCTP.FOR
*COPY STRUCTP
*
SUBROUTINE STRUCTP(XI,SCALE2,P2,IP2,UV,DV,US,DS,SS,CS,BS,TS,GL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
END
c SUBROUTINE PDFSET(PARAM,VALUE)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c DIMENSION PARAM(20),VALUE(20)
c CHARACTER*20 PARAM
c END
c
c*$ CREATE STRUCTM.FOR
c*COPY STRUCTM
c*
c SUBROUTINE STRUCTM(XI,SCALE2,UV,DV,US,DS,SS,CS,BS,TS,GL)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c END
c
c*$ CREATE STRUCTP.FOR
c*COPY STRUCTP
c*
c SUBROUTINE STRUCTP(XI,SCALE2,P2,IP2,UV,DV,US,DS,SS,CS,BS,TS,GL)
c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
c END
c
*$ CREATE DT_DIQBRK.FOR
*COPY DT_DIQBRK
*
......@@ -1941,7 +1941,7 @@ C...Commonblocks.
COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
C...HEPEVT commonblock.
PARAMETER (NMXHEP=9990)
PARAMETER (NMXHEP=HEPEVT_SIZE_REPLACE)
COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
&JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
DOUBLE PRECISION PHEP,VHEP
......@@ -530,5 +530,5 @@ c------------------------------------------------------------------------
* betfom,alpfom,alpfomi,gamfom
integer idlead,ilprtg
common/crvar/idlead,ilprtg
integer iLHC
common/LHCtune/iLHC
integer iLHC,ipytune
common/LHCtune/iLHC,ipytune
......@@ -131,7 +131,7 @@ c some basic things
engmax=0 !maximum energy
iologe=-1 !0=linear bins 1=log bins (cms engy)
iologl=-1 !0=linear bins 1=log bins (Kinetic engy)
infragm=0 !nuclear fragmentation
infragm=2 !nuclear fragmentation
c printout options
......@@ -492,6 +492,7 @@ c ncolmx=100000 !maximum number of collisions
bminim=0. !minimum impact parameter
phimax=2*3.1415927 !maximum phi
phimin=0 !minimum phi
ipytune=350 !Pythia Perugia Tune (2011) for PYTHIA run
iLHC=0 !LHC tune on(1)/off(0)
ionudi=3 !nuclear diffraction included (>0) or not (0)
! = 0 for RHIC nuclear data based on glauber
......@@ -1144,6 +1145,10 @@ c fill first hep stack before ordering
itar=0
pprojin=0d0
ptargin=0d0
imolim=0
if(istmax.eq.0
& .or.(model.ne.1.and.model.ne.5.and.mod(model,4).ne.0)) !model will full list of particles
&imolim=maproj+matarg
do i=1,nptl
......@@ -1181,7 +1186,7 @@ c i ............. particle number
if(io.gt.0)then
if(istptl(io).le.1.and.i.gt.maproj+matarg.and..not.lclean)then!mother is normal particle (incl. spectators and fragments)
iadd=1
if(iepo2hep(io).gt.maproj+matarg)then
if(iepo2hep(io).gt.imolim)then
jm1hep=iepo2hep(io)
if(jorptl(i).gt.0)jm2hep=iepo2hep(jorptl(i))
endif
......@@ -1307,11 +1312,15 @@ c print *,'clean',i,iorptl(io),pptl(3,i),jm1io
else
jm1hep=-1
jm2hep=-1
jd1hep=ifrptl(1,i)
jd2hep=ifrptl(2,i)
if(istptl(i).eq.2)then !beam remnant used in core
jm2hep=1
else
jd1hep=ifrptl(1,i)
jd2hep=ifrptl(2,i)
endif
endif
if(istptl(i).gt.1)goto 100 !skip non final particles (allowed before to define correctly some spectators going to core)
if(istptl(i).gt.1.and.i.gt.maproj+matarg)goto 100 !skip non final particles (allowed before to define correctly some spectators going to core)
idpdg=idtrafo('nxs','pdg',idptl(i))
if(idpdg.eq.99)then
......@@ -1364,7 +1373,7 @@ c pptl(5,i) ..... particle mass (GeV/c2)
phep2(5,nhep2)=dble(pptl(5,ii))
c istptl(i) ..... generation flag: last gen. (0) or not (1)
isthep2(nhep2)=min(2,istptl(ii)+1) !in hep:1=final, 2=decayed
if(i.le.maproj+matarg.and.iout.lt.0)isthep2(nhep2)=4 !beam particles
if(i.le.maproj+matarg)isthep2(nhep2)=4 !beam particles
c ityptl(i) ..... particle type (string, remnant ...)
c xorptl(1,i) ... x-component of formation point (fm)
vhep2(1,nhep2)=xorptl(1,ix)*1e-12 !conversion to mm
......@@ -1385,6 +1394,10 @@ c iorptl(i) ..... particle number of father (if .le. 0 : no father)
c jorptl(i) ..... particle number of mother (if .le. 0 : no mother)
jmohep2(2,nhep2)=jm2
c write(ifch,130)jmohep2(1,nhep2),jmohep2(2,nhep2),nhep2
c &,jdahep2(1,nhep2),jdahep2(2,nhep2),idhep2(nhep2)
c &,isthep2(nhep2),(phep2(k,nhep2),k=1,5),(vhep2(k,nhep2),k=1,4)
enddo
100 continue
......@@ -1458,12 +1471,11 @@ c Target
jmohep(2,nhep)=-1
endif
if(istmax.gt.0.and.model.eq.1)then
if(imolim.eq.0)then
nhep=maproj+matarg !beam particles will be copied later
if(iout.ge.0)nhep=nhep+2 !add new beam particles
endif
nhep0=0
nskip=0
if(iout.ge.0)nskip=-2 !add 2 beam particles and may be skip nucleons
......@@ -1471,7 +1483,7 @@ c Target
c reorder individual beam particles to make list of mothers for a given daughter
do j=1,maproj+matarg
if(istmax.eq.0.or.model.ne.1)then
if(imolim.ne.0)then
if(iout.lt.0)then
......@@ -1634,11 +1646,12 @@ c copy all other particles (secondary particles and spectators)
c look for mother of fragments
nhepi0=nhepio+1
if(jmohep2(1,k).gt.0.and.jmohep2(1,k).le.maproj+matarg)then
if(nhep0.ne.0.and.jmohep2(1,k).gt.0
& .and.jmohep2(1,k).le.maproj+matarg)then
c copy all mothers before the daughter to complete the beam remnant list
do j=1,maproj+matarg
if(isthep2(j).gt.0.and.jdahep2(1,j).eq.ihep2epo2(k))then
if(jdahep2(1,j).eq.ihep2epo2(k).and.isthep2(j).gt.0)then
nhepio=nhepio+1
idhep(nhepio)=idhep2(j)
phep(1,nhepio)=phep2(1,j)
......@@ -1729,19 +1742,20 @@ c copy all mothers before the daughter to complete the beam remnant list
enddo !iround
if(nhep.ne.nhep2-nskip.or.nhepio.ne.maproj+matarg-nskip)then
if(nhep+nskip.ne.nhep2.or.nhepio+nskip.ne.maproj+matarg)then
print *,'Warning : number of particles changed after copy'
& ,nhepio+nskip,maproj+matarg,nskip
nrem1=0
do k=1,nhep2
if(isthep2(k).eq.-4)then
nrem1=nrem1+1
endif
if(isthep2(k).eq.-4)then
nrem1=nrem1+1
endif
if(isthep2(k).gt.0)
& print *,' ',k,idhep2(k),jmohep2(1,k),isthep2(k)
& ,'from',ihep2epo2(k)
enddo
print *,' ',nhep2,'->',nhep
print *,' ',nhep2-nskip,'->',nhep
nrem2=0
do k=1,nhep
if((iout.ge.0.and.isthep(k).eq.3)
......@@ -1752,6 +1766,7 @@ c copy all mothers before the daughter to complete the beam remnant list
enddo
print *,' Particle list not consistent, skip event !'
print *,' ',nrem1,'->',nrem2
c stop
goto 10000
endif
......@@ -1760,6 +1775,7 @@ c update daughter list with correct index
do j=1,nhep