IAP GITLAB

Commit 72002831 authored by Tanguy Pierog's avatar Tanguy Pierog

correct nucleus definition (for input) to be compatible with PDG (and output...

correct nucleus definition (for input) to be compatible with PDG (and output code for nuclei) as : Z*10000+A*10
use realisitc nuclear fragmentation as default and get nuclear fragment as output particles (in EPOS, QGSJET01, QGSJETII and SIBYLL)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@3041 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 6c37d84a
......@@ -87,8 +87,8 @@ CRMCoptions::ParseOptions(int argc, char** argv)
("model,m", po::value<int>(), model_desc.str().c_str())
("projectile-momentum,p", po::value<double>(), "momentum/(GeV/c)")
("target-momentum,P", po::value<double>(), "momentum/(GeV/c)")
("projectile-id,i", po::value<int>(), "PDG or A*10000+Z*10")
("target-id,I", po::value<int>(), "PDG or A*10000+Z*10")
("projectile-id,i", po::value<int>(), "PDG or Z*10000+A*10")
("target-id,I", po::value<int>(), "PDG or Z*10000+A*10")
("config,c", po::value<string>(), "config file")
("out,f", po::value<string>(), "output file name (auto if none provided)")
;
......
......@@ -274,6 +274,8 @@ c open(ifcx,file=fnch(1:nfnch),status='unknown')
ecms=sngl(iecms) !center of mass energy in GeV/c2
c pnll=pproj !beam momentum GeV/c
infragm=2 !nuclear fragmentation (realistic)
c Projecticle definitions
if (abs(ipart) .eq. 1) then
c proton
......@@ -303,8 +305,8 @@ c pi-
c nuclei
elseif (ipart.gt.10000)then
idprojin=1120
laproj=mod(ipart,10000)/10 !proj Z
maproj=mod(ipart,10000000)/10000 !proj A
maproj=mod(ipart,10000)/10 !proj A
laproj=mod(ipart,10000000)/10000 !proj Z
c PDG
else
idprojin=idtrafo('pdg','nxs',ipart)
......@@ -342,8 +344,8 @@ c lead
c nuclei
elseif (itarg.gt.10000)then
idtargin=1120
latarg=mod(itarg,10000)/10 !targ Z
matarg=mod(itarg,10000000)/10000 !targ A
matarg=mod(itarg,10000)/10 !targ A
latarg=mod(itarg,10000000)/10000 !targ Z
c PDG
elseif (abs(itarg).eq.2112.or.abs(itarg).eq.2212)then
idtargin=idtrafo('pdg','nxs',itarg)
......
......@@ -49,6 +49,10 @@ c entry
ttauz=ttaus
c skip nuclei
if(idptl(i).gt.1e9)return
c small droplet decay
if(iabs(idptl(i)).gt.1e8)then
......
......@@ -1718,36 +1718,42 @@ c-----------------------------------------------------------------------
do 27 imx=1,mxmx
do 27 ima=2,mxma
if(iabs(id).eq.idmx(ima,imx))jj=ima
27 continue
if(id.gt.0)then
ic(1)=icre1(ii,jj)
ic(2)=icre2(ii,jj)
else
ic(2)=icre1(ii,jj)
ic(1)=icre2(ii,jj)
endif
if(ic(1).eq.100000.and.ic(2).eq.100000.and.rangen().lt.0.5)
$ then
ic(1)=010000
ic(2)=010000
endif
elseif(mod(id/10**8,10).eq.8)then
ic(1)=mod(id,10**8)/10000*100
ic(2)=mod(id,10**4)*100
27 continue
if(id.gt.0)then
ic(1)=icre1(ii,jj)
ic(2)=icre2(ii,jj)
else
write(ifch,*)'***** id: ',id
call utstop('idtr4: unrecognized id&')
ic(2)=icre1(ii,jj)
ic(1)=icre2(ii,jj)
endif
return
9998 continue
write(ifch,*)'id: ',id
call utstop('idtr4: indx=0.&')
if(ic(1).eq.100000.and.ic(2).eq.100000.and.rangen().lt.0.5)
$ then
ic(1)=010000
ic(2)=010000
endif
elseif(mod(id/10**8,10).eq.8)then
ic(1)=mod(id,10**8)/10000*100
ic(2)=mod(id,10**4)*100
elseif(id/10**9.eq.1.and.mod(id,10).eq.0)then !nuclei
nstr=mod(id,10**8)/10000000
npro=mod(id,10**7)/10000
nneu=mod(id,10**4)/10
ic(1)=(2*npro+nneu)*10**5+(2*nneu+npro)*10**4+nstr*10**3
ic(2)=0
else
write(ifch,*)'***** id: ',id
call utstop('idtr4: unrecognized id&')
endif
return
9999 continue
write(ifch,*)'id: ',id
call utstop('idtr4: ix out of range.&')
end
9998 continue
write(ifch,*)'id: ',id
call utstop('idtr4: indx=0.&')
9999 continue
write(ifch,*)'id: ',id
call utstop('idtr4: ix out of range.&')
end
c-----------------------------------------------------------------------
integer function idtra(ic,ier,ires,imix)
......@@ -2342,8 +2348,8 @@ c print *,'idtrafo',' ',code1,' ',code2,idi
idtrafo=1000000000+mod(id1,100)*10000+(id1/100)*10
return
elseif(i.eq.5.and.id1.gt.1004)then !nucleus from Sibyll
id1=(id1-1000)*100
idtrafo=1000000000+id1/2*10000+(id1/100)*10
id1=(id1-1000)
idtrafo=1000000000+id1/2*10000+id1*10
return
elseif(id1.eq.130.and.i.eq.2)then
idtrafo=-20
......@@ -2361,8 +2367,8 @@ c print *,'idtrafo',' ',code1,' ',code2,idi
idtrafo=1000000000+mod(id1,100)*10000+(id1/100)*10
return
elseif(i.eq.5.and.id1.gt.1004)then !nucleus from Sibyll
id1=(id1-1000)*100
idtrafo=1000000000+id1/2*10000+(id1/100)*10
id1=(id1-1000)
idtrafo=1000000000+id1/2*10000+id1*10
return
elseif(id1.eq.-20.and.i.eq.1)then
idtrafo=130
......
......@@ -5472,7 +5472,7 @@ c-----------------------------------------------------------------------
COMMON /Q_DEBUG/ DEBUG
SAVE
IF(DEBUG.GE.2)WRITE (MONIOU,201)NS-1,AW,G
IF(DEBUG.GE.2)WRITE (MONIOU,201)AW,NS-1,G
201 FORMAT(2X,'IXXSON - POISSON DITR.: AVERAGE AW=',E10.3,
* ' MAXIMAL VALUE NS=',I2,' RANDOM NUMBER G=',E10.3)
W=EXP(-AW)
......
......@@ -125,15 +125,111 @@ c IAF(i) - mass of the i-th fragment
call conqgs
c keep the projectile spectators as fragments
if(infragm.eq.2)then
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
if(ish.ge.7)write(ifch,'(a,i5,a,i5)')
$ ' Projecticle Fragment ',is,' Mass :',IAF(is)
nptl=nptl+1
istptl(nptl)=0
if(IAF(is).eq.1)then
id=idptl(is)
pptl(3,nptl)=pptl(3,is)
pptl(4,nptl)=pptl(4,is)
pptl(5,nptl)=pptl(5,is)
else
if(IAF(is).eq.2)then
id=17
elseif(IAF(is).eq.3)then
id=18
elseif(IAF(is).eq.4)then
id=19
else
inucl=IAF(is)
iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
id=1000000000+iprot*10000+inucl*10 !code for nuclei
endif
call idmass(id,am)
pptl(4,nptl)=dble(IAF(is))*pptl(4,is) !Etot
pptl(5,nptl)=am !mass
pz2tmp=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
if(pz2tmp.gt.0.d0)then
pptl(3,nptl)=sqrt(pz2tmp) !Pz
else
write(*,*)'Warning in emsqgs !'
write(*,*)'energy of fragment too small :',IAF(is),am
& ,pptl(4,nptl)
pptl(3,nptl)=pptl(4,nptl)
endif
endif
pptl(1,nptl)=0.d0 !P_x
pptl(2,nptl)=0.d0 !P_y
ityptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj+matarg
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=0.d0
xorptl(2,nptl)=0.d0
xorptl(3,nptl)=0.d0
xorptl(4,nptl)=0.d0
tivptl(1,nptl)=0.d0
tivptl(2,nptl)=0.d0
idptl(nptl)=id
if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
$ ' Fragment from qgsjet ',nptl,' id :',idptl(nptl)
$ , ' momentum :',(pptl(k,nptl),k=1,5)
enddo
endif
c make the projectile spectators as free nucleons
ns=0
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
ns=ns+IAF(is)
enddo
do is=1,ns !make the ns first projectile nucleon actives
istptl(is)=0
enddo
else
ns=0
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
ns=ns+IAF(is)
enddo
if(infragm.eq.1)then
c remaining nucleus is one fragment
nptl=nptl+1
istptl(nptl)=0
pptl(1,nptl)=0.d0
pptl(2,nptl)=0.d0
pptl(4,nptl)=0.d0
inucl=0
do is=1,ns
inucl=inucl+1
pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
enddo
iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
call idmass(idnucl,am)
pptl(5,nptl)=am !mass
ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
pptl(3,nptl)=sqrt(ptot)
ityptl(nptl)=0
istptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=xorptl(1,1)
xorptl(2,nptl)=xorptl(2,1)
xorptl(3,nptl)=xorptl(3,1)
xorptl(4,nptl)=xorptl(4,1)
tivptl(1,nptl)=tivptl(1,1)
tivptl(2,nptl)=tivptl(2,1)
idptl(nptl)=idnucl
else
do is=1,ns !make the ns first projectile nucleon actives
istptl(is)=0
enddo
endif
endif
endif
c restore target spectators
......@@ -447,7 +543,7 @@ c--------------------------------------------------------------------
c--------------------------------------------------------------------
c returns number ne of nonempty characters of line(j+1:160)
c Random number generator
c--------------------------------------------------------------------
double precision b10,drangen
include 'epos.inc'
......
......@@ -132,15 +132,112 @@ c iaf(i) - mass of the i-th fragment
call conqgsII
c keep the projectile spectators as fragments
if(infragm.eq.2)then
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
if(ish.ge.7)write(ifch,'(a,i5,a,i5)')
$ ' Projecticle Fragment ',is,' Mass :',IAF(is)
nptl=nptl+1
istptl(nptl)=0
if(IAF(is).eq.1)then
id=idptl(is)
pptl(3,nptl)=pptl(3,is)
pptl(4,nptl)=pptl(4,is)
pptl(5,nptl)=pptl(5,is)
else
if(IAF(is).eq.2)then
id=17
elseif(IAF(is).eq.3)then
id=18
elseif(IAF(is).eq.4)then
id=19
else
inucl=IAF(is)
iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
id=1000000000+iprot*10000+inucl*10 !code for nuclei
endif
call idmass(id,am)
pptl(4,nptl)=dble(IAF(is))*pptl(4,is) !Etot
pptl(5,nptl)=am !mass
pz2tmp=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
if(pz2tmp.gt.0.d0)then
pptl(3,nptl)=sqrt(pz2tmp) !Pz
else
write(*,*)'Warning in emsqgsII !'
write(*,*)'energy of fragment too small :',IAF(is),am
& ,pptl(4,nptl)
pptl(3,nptl)=pptl(4,nptl)
endif
endif
pptl(1,nptl)=0.d0 !P_x
pptl(2,nptl)=0.d0 !P_y
ityptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj+matarg
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=0.d0
xorptl(2,nptl)=0.d0
xorptl(3,nptl)=0.d0
xorptl(4,nptl)=0.d0
tivptl(1,nptl)=0.d0
tivptl(2,nptl)=0.d0
idptl(nptl)=id
if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
$ ' Fragment from qgsjetII ',nptl,' id :',idptl(nptl)
$ , ' momentum :',(pptl(k,nptl),k=1,5)
enddo
endif
c make the projectile spectators as free nucleons
ns=0
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
ns=ns+IAF(is)
enddo
do is=1,ns !make the ns first projectile nucleon actives
istptl(is)=0
enddo
else
ns=0
if(NSF.gt.0)then
do is=1,NSF !count the number of spectators
ns=ns+IAF(is)
enddo
if(infragm.eq.1)then
c remaining nucleus is one fragment
nptl=nptl+1
istptl(nptl)=0
pptl(1,nptl)=0.d0
pptl(2,nptl)=0.d0
pptl(4,nptl)=0.d0
inucl=0
do is=1,ns
inucl=inucl+1
pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
enddo
iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
call idmass(idnucl,am)
pptl(5,nptl)=am !mass
ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
pptl(3,nptl)=sqrt(ptot)
ityptl(nptl)=0
istptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=xorptl(1,1)
xorptl(2,nptl)=xorptl(2,1)
xorptl(3,nptl)=xorptl(3,1)
xorptl(4,nptl)=xorptl(4,1)
tivptl(1,nptl)=tivptl(1,1)
tivptl(2,nptl)=tivptl(2,1)
idptl(nptl)=idnucl
else
do is=1,ns !make the ns first projectile nucleon actives
istptl(is)=0
enddo
endif
endif
endif
c restore target spectators
......@@ -392,7 +489,7 @@ c--------------------------------------------------------------------
c--------------------------------------------------------------------
c returns number ne of nonempty characters of line(j+1:160)
c random number generator
c--------------------------------------------------------------------
double precision b10,drangen
include 'epos.inc'
......
......@@ -246,7 +246,7 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
CALL SIBNUC (IAP, itrg, engy)
if(ish.ge.5)write(ifch,'(a,i5)')
$ ' number of particles from Sibyll :',NPA
do k=1,NPA
do 100 k=1,NPA
c LLIST is the code of final particle, P - its 4-momentum and mass.
ic=LLA(k)
......@@ -256,60 +256,116 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
$ , ' momentum :',(PA(i,k),i=1,5)
if(ic.ge.1001) then !count spectators
ns=ns+ic-1000
else
nptl=nptl+1
if(nptl.gt.mxptl)call utstop('Sibyll: mxptl too small&')
id=idtrafo('sib','nxs',ic)
if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
nNuc=0
if(ic.ge.1001) then !count spectators
nNuc=ic-1000
if(infragm.le.1
& .or.dble(PA(4,k))/dble(nNuc).lt.egymin)then !nuclear interaction only above min energy, otherwise : fragmentation
ns=ns+nNuc
goto 100
elseif(ic.eq.1001)then
if(drangen(dummy).lt.0.45d0) then
ic = 13
else
ic = 14
endif
nNuc=0
else
ptm=sqrt(PA(1,k)*PA(1,k)+PA(2,k)*PA(2,k)+PA(5,k)*PA(5,k))
PA(4,k)=PA(4,k)*float(nNuc) !energy by nucleon
PA(3,k)=sign(sqrt((PA(4,k)+ptm)*(PA(4,k)-ptm)),PA(3,k))
endif
endif
nptl=nptl+1
if(nptl.gt.mxptl)call utstop('Sibyll: mxptl too small&')
id=idtrafo('sib','nxs',ic)
if(ish.ge.7)write(ifch,'(a,i5,a,i10,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion'
if(abs(id).gt.1000)nbar=nbar+sign(1,id)
pptl(1,nptl)=PA(1,k) !P_x
pptl(2,nptl)=PA(2,k) !P_y
pptl(3,nptl)=PA(3,k) !P_z
pptl(4,nptl)=PA(4,k) !E
pptl(5,nptl)=PA(5,k) !mass
istptl(nptl)=0
ityptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj+matarg
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=0.
xorptl(2,nptl)=0.
xorptl(3,nptl)=0.
xorptl(4,nptl)=0.
tivptl(1,nptl)=0.
tivptl(2,nptl)=0.
idptl(nptl)=id
nbar=nbar+nNuc
if(abs(id).gt.1000.and.nNuc.eq.0)nbar=nbar+sign(1,id)
pptl(1,nptl)=PA(1,k) !P_x
pptl(2,nptl)=PA(2,k) !P_y
pptl(3,nptl)=PA(3,k) !P_z
pptl(4,nptl)=PA(4,k) !E
pptl(5,nptl)=PA(5,k) !mass
istptl(nptl)=0
ityptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj+matarg
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=0.
xorptl(2,nptl)=0.
xorptl(3,nptl)=0.
xorptl(4,nptl)=0.
tivptl(1,nptl)=0.
tivptl(2,nptl)=0.
idptl(nptl)=id
if(ish.ge.5)write(ifch,'(a,i5,a,i5,a,4(e10.4,1x),f6.3)')
if(ish.ge.5)write(ifch,'(a,i5,a,i10,a,4(e10.4,1x),f6.3)')
$ ' particle from Sibyll ',nptl,' id :',idptl(nptl)
$ , ' momentum :',(pptl(i,nptl),i=1,5)
endif
100 continue
enddo
ntw=nbar-(maproj-ns)
nsf=0
if(ntw.lt.matarg)then
nts=matarg-ntw
do is=maproj+1,maproj+nts !make the nts first target nucleon actives (not wounded)
istptl(is)=0
enddo
else
ns=maproj+matarg-nbar
nsf=maproj+matarg-nbar
endif
if(ns.le.maproj)then
do is=1,ns !make the ns first projectile nucleon actives (not wounded)
istptl(is)=0
enddo
if(ish.ge.5)write(ifch,'(a,i3,a,i3,a,i2,a)')
$ ' target spectators :',matarg-ntw
$ ,' projectile spectators (ns) :',nsf,' (',ns,')'
if((infragm.le.1.or.nsf.gt.ns).and.nsf.le.maproj)then
if(infragm.eq.2)ns=nsf-ns
if(infragm.eq.1.and.ns.gt.0)then
c remaining nucleus is one fragment
nptl=nptl+1
istptl(nptl)=0
pptl(1,nptl)=0.d0
pptl(2,nptl)=0.d0
pptl(4,nptl)=0.d0
inucl=0
do is=1,ns
inucl=inucl+1
pptl(4,nptl)=pptl(4,nptl)+pptl(4,is)
enddo
iprot= int(dble(inucl) / 2.15d0 + 0.7d0)
idnucl=1000000000+iprot*10000+inucl*10 !code for nuclei
call idmass(idnucl,am)
pptl(5,nptl)=am !mass
ptot=(pptl(4,nptl)+am)*(pptl(4,nptl)-am)
pptl(3,nptl)=sqrt(ptot)
ityptl(nptl)=0
istptl(nptl)=0
iorptl(nptl)=1
jorptl(nptl)=maproj
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
xorptl(1,nptl)=xorptl(1,1)
xorptl(2,nptl)=xorptl(2,1)
xorptl(3,nptl)=xorptl(3,1)
xorptl(4,nptl)=xorptl(4,1)
tivptl(1,nptl)=tivptl(1,1)
tivptl(2,nptl)=tivptl(2,1)
idptl(nptl)=idnucl
else
do is=1,ns !make the nsf first projectile nucleon actives (not wounded)
istptl(is)=0
enddo
endif
endif
endif
1000 return
1001 iret=-1
......
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