IAP GITLAB

Commit 3eb5b9c3 authored by Colin Baus's avatar Colin Baus

decays by pythia. init to treat properly the pdgid. copy particles status<istmax parameter

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4137 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent bddc3963
......@@ -8,30 +8,108 @@ c-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
include 'epos.inc'
c change to center of mass system
c common block var iframe renamed to dframe since it is used by epos
c common block var iframe renamed to idframe since it is used by epos
LOGICAL LEMCCK,LHADRO,LSEADI,LEVAPO
COMMON /DTFLG1/ IFRAG(2),IRESCO,IMSHL,IRESRJ,IOULEV(6),
& LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,DFRAME,ITRSPT
& LEMCCK,LHADRO(0:9),LSEADI,LEVAPO,IDFRAME,ITRSPT
COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
C switch off decays
MSTJ(21)=0
DFRAME=2
C set some variables for decays handled by epos
egymin=10. !?
egymax=1e7 !?
irescl=0 !momentum rescale
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
egymin=10. !min energy for model
egymax=5e4 !max energy for model
irescl=0 !don't rescale/skip events with wrong energy
END
c-----------------------------------------------------------------------
subroutine IniEvtDpm
c-----------------------------------------------------------------------
c Setting energy, primaries,... for each event. Usefull for e.g. conex
c Setting energy, primaries,... for each event class. Useful for e.g. conex
c-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
include 'epos.inc'
......@@ -41,22 +119,32 @@ c-----------------------------------------------------------------------
COMMON /DTPART/ ANAME(210),AAM(210),GA(210),TAU(210),
& IICH(210),IIBAR(210),K1(210),K2(210)
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDP
c fresh common block
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
DOUBLE PRECISION EPROJ
COMMON /DPMEVTINI/ EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP
COMMON /DPMEVTINI/ EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
C general initialization
NCASES = -1
NCASES = -1 !skip reading steering cards
EPROJ = dble(elab)
NPMASS = maproj
NPCHAR = laproj
NTMASS = matarg
NTCHAR = latarg
IDP = 1
IF(laproj.eq.-1) THEN !not nucleus
IDPDG = idtrafo('nxs','pdg',idproj)
ELSE
c will be treated as nucleus in DT_INIT and NPMASS,... will be used
IDPDG = 0
ENDIF
IGLAU = 0
EPROJ = (ECMS**2-AAM(NPMASS)**2-AAM(NTMASS)**2)
& /(2.0D0*AAM(NTMASS))
CALL DT_INIT(NCASES,EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,IGLAU)
CALL DT_INIT(NCASES,EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
+,IGLAU)
c initialize cross-sections by calling epos-bas.f function -> models.F -> dpmjet-crmc.f
c initialize cross-sections by calling epos-bas.f function -> models.F -> dpmjet-crmc.f
call xsigma
END
......@@ -87,10 +175,10 @@ c-----------------------------------------------------------------------
* properties of interacting particles
COMMON /DTPRTA/ IT,ITZ,IP,IPZ,IJPROJ,IBPROJ,IJTARG,IBTARG
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDP
INTEGER NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
INTEGER IREJ, KKMAT
DOUBLE PRECISION EPROJ
COMMON /DPMEVTINI/ EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDP
COMMON /DPMEVTINI/ EPROJ,NPMASS,NPCHAR,NTMASS,NTCHAR,IDPDG
integer nptlhep(NMXHKK)
anintine=anintine+1.
......@@ -102,6 +190,11 @@ c-----------------------------------------------------------------------
KKMAT = -1
C call dpmjet event simulation
IF (IDPDG.EQ.0) THEN
IDP = 1
ELSE
IDP = IDT_ICIHAD(IDPDG)
ENDIF
CALL dt_kkinc(NPMASS,NPCHAR,NTMASS,NTCHAR,IDP,EPROJ, KKMAT,IREJ)
IF(IREJ.NE.0) THEN
WRITE(ifch,'(a)')
......@@ -109,6 +202,8 @@ C call dpmjet event simulation
GOTO 1001
ENDIF
call dt_evtout(1)
c nevt .......... error code. 1=valid event, 0=invalid event
c bimevt ........ absolute value of impact parameter
c phievt ........ angle of impact parameter
......@@ -166,8 +261,9 @@ c LIST is the code of final particle, P - its 4-momentum and mass.
$ ' DPMJET particle ',k,' id :',ic,' before conversion'
$ , ' momentum :',(sngl(PHKK(i,k)),i=1,5)
IF(ISTHKK(k).EQ.1)THEN !! if final particle
c final particle
imaxhepstatus = MIN(ISTMAX,2)
print*,imaxhepstatus
IF(ISTHKK(k).GE.1 .AND. ISTHKK(k).LE.imaxhepstatus)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...
......@@ -187,10 +283,6 @@ c final particle
pptl(4,nptl)=sngl(PHKK(4,k)) !E
pptl(5,nptl)=sngl(PHKK(5,k)) !mass
ityptl(nptl)=0
iorptl(nptl)=1
c if(JMOHKK(1,k).ne.0)iorptl(nptl)=nptlhep(JMOHKK(1,k))
jorptl(nptl)=maproj+matarg
c if(JMOHKK(2,k).ne.0)jorptl(nptl)=nptlhep(JMOHKK(2,k))
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=sngl(VHKK(1,k))*1e12
......@@ -201,6 +293,20 @@ c if(JMOHKK(2,k).ne.0)jorptl(nptl)=nptlhep(JMOHKK(2,k))
tivptl(2,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
c if(JMOHKK(1,k).ne.0.AND.nptlhep(JMOHKK(1,k)).NE.0) then !mother
c iorptl(nptl)=nptlhep(JMOHKK(1,k))
c endif
c if(JMOHKK(2,k).ne.0.AND.nptlhep(JMOHKK(2,k)).NE.0) then !father
c jorptl(nptl)=nptlhep(JMOHKK(2,k))
c endif
print*,istmax,k,JMOHKK(1,k),nptlhep(k),iorptl(nptl),
$ nptlhep(JMOHKK(1,k))
if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,i5,a,4(e10.4,1x),f6.3)')
$ ' particle from DPMJET ',nptl,' id :',idptl(nptl)
......
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