IAP GITLAB

Commit d92286ea authored by Tanguy Pierog's avatar Tanguy Pierog

correct link to QGSJETII-04 to get proper energy balance (problem with boost...

correct link to QGSJETII-04 to get proper energy balance (problem with boost and mass of particles different in QGSJETII)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4855 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 67586912
...@@ -13,8 +13,8 @@ OPTION (__HIJING__ "Build with model" OFF) ...@@ -13,8 +13,8 @@ OPTION (__HIJING__ "Build with model" OFF)
OPTION (__SIBYLL__ "Build with model" ON) OPTION (__SIBYLL__ "Build with model" ON)
OPTION (__PHOJET__ "Build with model" OFF) OPTION (__PHOJET__ "Build with model" OFF)
OPTION (__DPMJET__ "Build with model" OFF) OPTION (__DPMJET__ "Build with model" OFF)
OPTION (__QGSJETII03__ "Build with model" OFF) OPTION (__QGSJETII03__ "Build with model" ON)
OPTION (__QGSJETII04__ "Build with model" ON) OPTION (__QGSJETII04__ "Build with model" OFF)
######################################ONLY EDIT THIS###################################### ######################################ONLY EDIT THIS######################################
if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) if (CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
......
...@@ -61,7 +61,8 @@ c QGSJET-II Common ...@@ -61,7 +61,8 @@ c QGSJET-II Common
if(icp.eq.0)icp=1-2*int(rangen()+0.5) !pi0=pi+ or p- if(icp.eq.0)icp=1-2*int(rangen()+0.5) !pi0=pi+ or p-
if(abs(icp).gt.5) if(abs(icp).gt.5)
& call utstop('Projectile not allowed in QGSJET-II !&') & call utstop('Projectile not allowed in QGSJET-II !&')
e0=dble(elab) c Important to calculate e0 (kinetic energy) in double precision
e0=dble(elab)-dble(amproj)
call qgini(e0,icp,maproj,matarg) call qgini(e0,icp,maproj,matarg)
call qgini(e0,icp,maproj,matarg) !again to set bm properly call qgini(e0,icp,maproj,matarg) !again to set bm properly
bmax=BMAXQGS bmax=BMAXQGS
...@@ -95,6 +96,11 @@ c nsf - number of secondary fragments; ...@@ -95,6 +96,11 @@ c nsf - number of secondary fragments;
c iaf(i) - mass of the i-th fragment c iaf(i) - mass of the i-th fragment
common /qgarr13/ nsf,iaf(iapmax) common /qgarr13/ nsf,iaf(iapmax)
common /qgarr55/ nwt,nwp common /qgarr55/ nwt,nwp
c particle masses for proper boost
double precision am0,amn,amk,amc,amlamc,amlam,ameta,ammu
*,dmmin,wex,dmres,wdres
common /qgarr10/ am0,amn,amk,amc,amlamc,amlam,ameta,ammu
common /qgarr21/ dmmin(3),wex(3),dmres(3),wdres(3)
iret=0 iret=0
b1=bminim b1=bminim
...@@ -289,6 +295,8 @@ c restore target spectators ...@@ -289,6 +295,8 @@ c restore target spectators
if(NWT.lt.matarg)then if(NWT.lt.matarg)then
c number of participants c number of participants
ntgevt=NWT ntgevt=NWT
c if SD pro (no wounded target nucleons) one nucleon included in QII list
if(NWT.eq.0.and.nint(typevt).eq.4)ntgevt=1
do is=ntgevt+1,matarg !make the last target nucleon final do is=ntgevt+1,matarg !make the last target nucleon final
iorptl(maproj+is)=0 iorptl(maproj+is)=0
istptl(maproj+is)=0 istptl(maproj+is)=0
...@@ -314,13 +322,31 @@ c k0l ...@@ -314,13 +322,31 @@ c k0l
id=idtrafo('qgs','nxs',ic) id=idtrafo('qgs','nxs',ic)
if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)') if(ish.ge.7)write(ifch,'(a,i5,a,i5,a)')
$ ' epos particle ',nptl,' id :',id,' after conversion' $ ' epos particle ',nptl,' id :',id,' after conversion'
call idmass(id,am) call idmass(id,amepo)
am=amepo !just in case particle is not in the following list
ica=abs(ic)
pptl(1,nptl)=real(esp(3,is)) !P_x if(ica.le.1)then
pptl(2,nptl)=real(esp(4,is)) !P_y am=am0
pptl(3,nptl)=real(esp(2,is)) !P_z elseif(ica.eq.2.or.ica.eq.3)then
pptl(4,nptl)=real(esp(1,is)) !E am=amn
elseif(ica.eq.4.or.ica.eq.5)then
am=amk
elseif(ica.eq.6)then !lambda
am=amlam
elseif(ica.eq.7.or.ica.eq.8)then !D meson
am=amc
elseif(ica.eq.9)then !lambdac
am=amlamc
elseif(ica.eq.10)then !eta
am=ameta
elseif(ica.eq.19)then !rho0
am=dmmin(1)
endif
esum=esum+esp(1,is)
pptl(1,nptl)=sngl(esp(3,is)) !P_x
pptl(2,nptl)=sngl(esp(4,is)) !P_y
pptl(3,nptl)=sngl(esp(2,is)) !P_z
pptl(4,nptl)=sngl(esp(1,is)) !E
pptl(5,nptl)=am !mass pptl(5,nptl)=am !mass
istptl(nptl)=0 istptl(nptl)=0
ityptl(nptl)=0 ityptl(nptl)=0
...@@ -340,6 +366,11 @@ c boost in CMS frame ...@@ -340,6 +366,11 @@ c boost in CMS frame
call utlob5(yhaha, pptl(1,nptl), pptl(2,nptl) call utlob5(yhaha, pptl(1,nptl), pptl(2,nptl)
. , pptl(3,nptl), pptl(4,nptl), pptl(5,nptl)) . , pptl(3,nptl), pptl(4,nptl), pptl(5,nptl))
c give particle proper mass ...
pptl(5,nptl)=amepo
pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2
. +pptl(3,nptl)**2+pptl(5,nptl)**2)
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,i5,a,4(e10.4,1x),f6.3)')
$ ' particle from qgsjet ',nptl,' id :',idptl(nptl) $ ' particle from qgsjet ',nptl,' id :',idptl(nptl)
$ , ' momentum :',(pptl(k,nptl),k=1,5) $ , ' momentum :',(pptl(k,nptl),k=1,5)
......
...@@ -57,7 +57,7 @@ c QGSJET-II Common ...@@ -57,7 +57,7 @@ c QGSJET-II Common
call iclass(idtarg,icltar) call iclass(idtarg,icltar)
icp=idtrafo('nxs','qgs',idproj) icp=idtrafo('nxs','qgs',idproj)
if(icp.eq.0)icp=1-2*int(rangen()+0.5) !pi0=pi+ or p- if(icp.eq.0)icp=1-2*int(rangen()+0.5) !pi0=pi+ or p-
e0=dble(elab) e0=dble(elab)-dble(amproj)
call qgini(e0,icp,maproj,matarg) call qgini(e0,icp,maproj,matarg)
call qgini(e0,icp,maproj,matarg) !again to set bm properly call qgini(e0,icp,maproj,matarg) !again to set bm properly
bmax=BMAXQGS bmax=BMAXQGS
......
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