IAP GITLAB

Skip to content

GitLab

C
CRMC
Commit bd00a54a authored by Tanguy Pierog's avatar Tanguy Pierog
Browse files
Options

update old DPMJET interface. Limit mass range to light nuclei to avoid problems wit DPMJET 3.0-6

parent 3a26d12e
Hide whitespace changes
Inline Side-by-side
Showing with 84 additions and 36 deletions
src/dpmjet/3.0-6/dpmjet-crmc.f
View file @ bd00a54a
......@@ -65,10 +65,22 @@ C general initialization
IF(laproj.eq.-1) THEN !not nucleus
IDPDG = idtrafo('nxs','pdg',idproj)
IDP=IDT_ICIHAD(IDPDG)
NPCHAR = IICH(IDP)
elseif(maproj.gt.20)then
call utstop("Only nucleus with A<20 as projectile !&")
ELSE
c will be treated as nucleus in DT_INIT and NPMASS,... will be used
IDPDG = 0
IDP = 0
ENDIF
if(abs(idtarg).ne.1120.and.idtarg.ne.0)then
call utstop("Only (anti)proton or nucleus as target !&")
elseif(matarg.gt.20)then
call utstop("Only nucleus with A<20 as target !&")
else
NTCHAR = sign(abs(latarg),idtarg)
endif
IGLAU = 0
CALL DT_INIT(NCASES,EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
......@@ -85,13 +97,12 @@ c set decay flag in Pythia for DPMJET (after DT_INIT otherwise default is used)
c particle without breit-wigner decay will be decayed in EPOS to get full history
do i=100,500
idpdg = KCHG(I,4)
idpdgi = KCHG(I,4)
c don't touch decays with unknown code
idepos = idtrafostatus('pdg','nxs',idpdg,istatus)
if(abs(idpdg).eq.10323)print *,idpdg,idepos,istatus
idepos = idtrafostatus('pdg','nxs',idpdgi,istatus)
if(istatus.ne.0) then
go to 100
go to 100
endif
c don't touch decays with unknown width
......@@ -166,13 +177,24 @@ c call conwr
IREJ = 0
KKMAT = -1
C call dpmjet event simulation
NDPMEVENT = NDPMEVENT+1 !needs to be updated for internal init. sets NEVHKK
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
NPMASS = maproj
IF(laproj.eq.-1) THEN !not nucleus
NPCHAR = IICH(IDP)
else
NPCHAR = laproj
endif
NTMASS = matarg
NTCHAR = sign(abs(latarg),idtarg)
CALL dt_kkinc(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPROJ, KKMAT,IREJ)
IF(IREJ.NE.0) THEN
WRITE(ifch,'(a)')
......@@ -247,35 +269,55 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
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
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 !"
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
goto 100
ist=19
endif
elseif(ic.eq.66666)then !cascade
id=80000000
else
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
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,i5,a)')
if(ish.ge.7)write(ifch,'(a,i5,a,i12,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
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
......@@ -297,18 +339,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)=0
iorptl(nptl)=0
jorptl(nptl)=0
istptl(nptl)=1
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
......@@ -320,18 +363,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)=0
iorptl(nptl)=0
jorptl(nptl)=0
istptl(nptl)=1
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
......@@ -344,11 +388,15 @@ c fix beam particles which are registered at rest in DPMJET
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
if(JMOHKK(2,k).ne.0.AND.nptlhep(JMOHKK(2,k)).NE.0) then !father
jorptl(nptl)=nptlhep(JMOHKK(2,k))
if(JMOHKK(1,k).ne.0)then
if(nptlhep(JMOHKK(1,k)).NE.0) then !mother
iorptl(nptl)=nptlhep(JMOHKK(1,k))
endif
endif
if(JMOHKK(2,k).ne.0)then
if(nptlhep(JMOHKK(2,k)).NE.0) then !father
jorptl(nptl)=nptlhep(JMOHKK(2,k))
endif
endif
if(jorptl(nptl).gt.0.and.iorptl(nptl).gt.jorptl(nptl))then !reorder mother/father to avoid problem with hepmc
io=iorptl(nptl)
......@@ -360,7 +408,7 @@ c fix beam particles which are registered at rest in DPMJET
C print*,istmax,k,JMOHKK(1,k),nptlhep(k),iorptl(nptl),
C $ nptlhep(JMOHKK(1,k))
if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,i5,a,4(e10.4,1x),f6.3)')
if(ish.ge.5)write(ifch,'(a,i5,a,i12,a,i5,a,4(e10.4,1x),f6.3)')
$ ' particle from DPMJET ',nptl,' id :',idptl(nptl)
$ ,' st :',istptl(nptl)
$ , ' momentum :',(pptl(i,nptl),i=1,5)
......
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

mailto ikp-admins [at] lists [dot] kit [dot] edu