IAP GITLAB

Commit 89aa2496 authored by Colin Baus's avatar Colin Baus

dpmjet decays handled in EPOS (for history) only for those decays that have...

dpmjet decays handled in EPOS (for history) only for those decays that have known id, known width and large enough width

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4940 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent f8e2825c
......@@ -48,6 +48,9 @@ c fresh common block
COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5)
COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
integer istatus
real taugm
C general initialization
NCASES = -1 !skip reading steering cards
EPROJ = dble(elab)
......@@ -71,11 +74,34 @@ c set decay flag in Pythia for DPMJET (after DT_INIT otherwise default is used)
call IniDkyJetset
c particle without breit-wigner decay will be decayed in EPOS to get full history
MDCY(13,1)=1 !force muon decay in Jetset
do i=100,500
if(PMAS(i,2).LT.PARP(41))MDCY(i,1)=0
c print *,i,PMAS(i,1),PMAS(i,2),KCHG(I,4),PMAS(i,2).GE.PARP(41)
enddo
idpdg = KCHG(I,4)
c don't touch decays with unknown code
eposid = idtrafostatus('pdg','nxs',idpdg,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)
c print *,i,(taugm.eq.indexOutOfRange.or.taugm.eq.iijjOutOfRange)
c + ,(PMAS(i,2).LT.PARP(41))
if(istatus.ne.0) then
go to 100
endif
c don't touch decays with too small width
if(PMAS(i,2).LT.PARP(41)) then
go to 100
endif
if(ish.ge.2) write(ifch,*)
+ "CF Id:",i, PMAS(i,1),PMAS(i,2),KCHG(I,4),
+ "Width ok:", PMAS(i,2).GE.PARP(41),
+ "decay in epos"
100 enddo !loop over 100-500 particle ids
c initialize cross-sections by calling epos-bas.f function -> models.F -> dpmjet-crmc.f
call xsigma
......@@ -207,7 +233,7 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
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
ist=31 !convert to status for epos string
else
ist=abs(ISTHKK(k)) !!0 means last generation other codes are e.g
endif
......@@ -475,7 +501,7 @@ c xprod seems to be prod + quasi-ela. this is why we subtract it again
* 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,
......
......@@ -1576,14 +1576,21 @@ c endif
c return
c end
c
c-----------------------------------------------------------------------
subroutine idtau(id,p4,p5,taugm)
subroutine idtaustatus(id,p4,p5,taugm,istatus)
c returns lifetime(c*tau(fm))*gamma for id with energy p4, mass p5
c-----------------------------------------------------------------------
include 'epos.inc'
parameter (mxindx=1000,mxre=100,mxma=11,mxmx=6)
common/crema/indx(mxindx),rema(mxre,mxma),rewi(mxre,mxma)
*,idmx(mxma,mxmx),icre1(mxre,mxma),icre2(mxre,mxma)
parameter (indexOutOfRange=-999.)
parameter (iijjOutOfRange=-9999.)
integer istatus=0
if(iabs(id).eq.14)then
wi=.197/658.654e15
elseif(iabs(id).eq.16)then
......@@ -1598,8 +1605,8 @@ c-----------------------------------------------------------------------
elseif(iabs(id).lt.1e8)then
ix=iabs(id)/10
if(ix.lt.1.or.ix.gt.mxindx)then
write(ifch,*)'id:',id
call utstop('idtau: ix out of range.&')
istatus=indexOutOfRange
return
endif
ii=indx(ix)
jj=mod(iabs(id),10)+2
......@@ -1616,8 +1623,8 @@ c-----------------------------------------------------------------------
endif
75 continue
if(ii.lt.1.or.ii.gt.mxre.or.jj.lt.1.or.jj.gt.mxma)then
write(ifch,*)'id,ii,jj:',id,' ',ii,jj
call utstop('idtau: ii or jj out of range&')
istatus=iijjOutOfRange
return
endif
wi=rewi(ii,jj)
else
......@@ -1644,6 +1651,29 @@ c-c tauz=amax1(.2,tauz)
return
end
c-----------------------------------------------------------------------
subroutine idtau(id,p4,p5,taugm)
c wrapper function for idtaustatus that also checks wether IDs exist
c returns lifetime(c*tau(fm))*gamma for id with energy p4, mass p5
c fails when index out of range. idtaustatus will return negative
c-----------------------------------------------------------------------
include 'epos.inc'
parameter (indexOutOfRange=-999.)
parameter (iijjOutOfRange=-9999.)
integer istatus=0
call idtaustatus(id,p4,p5,taugm,istatus)
if(istatus.eq.indexOutOfRange) then
write(ifch,*)'id:',id
call utstop('idtau: ix out of range.&')
endif
if(istatus.eq.iijjOutOfRange) then
write(ifch,*)'id,ii,jj:',id,' ',ii,jj
call utstop('idtau: ii or jj out of range&')
endif
return
end
c-----------------------------------------------------------------------
subroutine idtr4(id,ic)
c transforms generalized paige_id -> werner_id (for < 4 flv)
......@@ -1896,10 +1926,25 @@ c-----------------------------------------------------------------------
endif
return
end
c------------------------------------------------------------------------------
integer function idtrafo(code1,code2,idi)
c------------------------------------------------------------------------------
c wrapper for idtrafo to catch if ID not found
character*3 code1,code2
integer istatus=0
idtrafo = idtrafostatus(code1,code2,idi,istatus)
if (istatus.ne.0) then
print *,'idtrafo: ',code1,' -> ', code2,idi,' not found.'
stop
endif
return
end
c------------------------------------------------------------------------------
integer function idtrafostatus(code1,code2,idi,istatus)
+ result(idtrafo)
c------------------------------------------------------------------------------
c.....tranforms id of code1 (=idi) into id of code2 (=idtrafo)
c.....supported codes:
c.....'nxs' = epos
......@@ -1910,6 +1955,8 @@ c.....'sib' = Sibyll
c.....'cor' = Corsika (GEANT)
c.....'flk' = Fluka
c.....returns status 0 if ID could be converted and 1 if it was not found in table
C --- ighenex(I)=EPOS CODE CORRESPONDING TO GHEISHA CODE I ---
common /ighnx/ ighenex(35)
......@@ -1971,6 +2018,7 @@ C IFCTABL CONVERTS FLUKA PARTICLES INTO CORSIKA PARTICLES
* 156, 157, 0, 0, 36*0/
c-------------------------------------------------------------------------------
integer istatus
character*3 code1,code2
parameter (ncode=5,nidt=341)
integer idt(ncode,nidt)
......@@ -2324,7 +2372,8 @@ c nxs|pdg|qgs|cor|sib
c print *,'idtrafo',' ',code1,' ',code2,idi
istatus=0
idtrafo=0
nidtmx=68
id1=idi
if(code1.eq.'nxs')then
......@@ -2473,10 +2522,9 @@ c print *,'idtrafo',' ',code1,' ',code2,idi
return
endif
print *, 'idtrafo: ',code1,' -> ', code2,id1,' not found. '
stop
istatus=1
c idtrafocx=0
c return
return
100 if(j.eq.4)then !corsika
if(idtrafo.eq.201)then
......
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