IAP GITLAB

Commit 1a5d5e0c authored by Tanguy Pierog's avatar Tanguy Pierog

put all particles on-shell in double precision using the mass stored in...

put all particles on-shell in double precision using the mass stored in phep(5) (solve the problems seen at high energy)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@5263 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent edb0decc
CRMC v1.5.6 Last modifications 2015/10/02
CRMC v1.5.7 Last modifications 2016/03/23
**********************************************************
The Program "crmc"
......@@ -36,7 +36,7 @@ Reference : C. Baus, T. Pierog and R. Ulrich. To be published (2016)
* Post LHC :
EPOS LHC (-m0) : T. Pierog et al.
arXiv:1306.0121 [hep-ph]
Phys.Rev. C92 (2015) 034906
QGSJETII-04 (-m7) : S.Ostapchenko,
Phys.Rev. D83 (2011) 014018
......
......@@ -60,7 +60,7 @@ c Calculations of energy of the center-of-mass in the detector frame
e1=dsqrt(dble(m1)**2+pproj**2)
e2=dsqrt(dble(m2)**2+ptarg**2)
decms=dsqrt(((e1+e2)+(pproj+ptarg))*(e1+e2)-(pproj+ptarg)))
decms=dsqrt(((e1+e2)+(pproj+ptarg))*((e1+e2)-(pproj+ptarg)))
c Later a rapidity boost back into the detector system will be performed
doBoost = .true.
c ycm2det defines this rapidity
......@@ -72,6 +72,7 @@ c ycm2det defines this rapidity
ycm2det=0.5d0*dlog(((e1+e2)+(pproj+ptarg))/
+ ((e1+e2)-(pproj+ptarg)))
endif
if (pproj .le. 0d0) then
ycm2det=-ycm2det
endif
......@@ -138,10 +139,12 @@ c Output quantities
integer outstat(*)
double precision boostvec1,boostvec2,boostvec3,boostvec4,boostvec5
double precision mass,amt
double precision mass,ppp,xcount
double precision ycm2det
logical doBoost
common/boostvars/ycm2det,doBoost
data xcount / 0d0 /
save
integer i!,k
......@@ -175,35 +178,36 @@ c define vec to boost from cm. to cms frame
boostvec5=1d0
c write(*,*)nevhep,nhep,boostvec3,boostvec4,ycm2det
do i=1,nhep
c Put particle is on-shell (needed because of single/double precision pb)
ppp = sqrt( phep(1,i)**2 + phep(2,i)**2 + phep(3,i)**2)
mass = (phep(4,i)+ppp)*(phep(4,i)-ppp)
if(abs(mass-phep(5,i)**2)/max(1d2,ppp**2).gt.1d-4) !not to count precision problems
+ xcount=xcount+1d0
mass=phep(5,i)
phep(4,i)=sqrt(phep(3,i)**2+phep(2,i)**2+phep(1,i)**2
+ +mass**2) !force particles to be on-shell
c boost output to cms frame
if(doBoost.eqv..true..and.ycm2det.ne.0d0) then
if(model.eq.6) then ! sibyll needs calculated mass because itlob5 does not work with off-shell particles
amt = dsqrt( phep(1,i)**2 + phep(2,i)**2 + phep(3,i)**2)
mass = dsqrt((phep(4,i)+amt)*(phep(4,i)-amt))
else
mass = phep(5,i)
endif
if(doBoost.eqv..true..and.ycm2det.ne.0d0) then
call utlob2(-1,boostvec1,boostvec2,boostvec3
+ ,boostvec4,boostvec5
+ ,vhep(1,i),vhep(2,i),vhep(3,i),vhep(4,i),-99)
call utlob5dbl(-ycm2det
+ ,phep(1,i), phep(2,i), phep(3,i), phep(4,i), mass)
else
phep(4,i)=sqrt(phep(3,i)**2+phep(2,i)**2+phep(2,i)**2
+ +phep(5,i)**2)
endif
outpart(i)=idhep(i)
outpx(i)=phep(1,i)
outpy(i)=phep(2,i)
outpz(i)=phep(3,i)
oute(i)=phep(4,i)
outm(i)=phep(5,i)
outstat(i)=isthep(i)
endif
outpart(i)=idhep(i)
outpx(i)=phep(1,i)
outpy(i)=phep(2,i)
outpz(i)=phep(3,i)
oute(i)=phep(4,i)
outm(i)=phep(5,i)
outstat(i)=isthep(i)
c write(*,'(4x,i6,1x,4(e12.6,1x))')idhep(i),(vhep(k,i),k=1,4)
c write(*,'(i5,3x,i2,2x,2i5,2x,2i5)')i,isthep(i)
c * ,jmohep(1,i),jmohep(2,i),jdahep(1,i),jdahep(2,i)
c write(*,'(i10,1x,4(e12.6,1x))')idhep(i),(phep(k,i),k=1,4)
enddo
if(ievent.eq.nevent.and.xcount.gt.0d0)print *,
+ 'Warning : negative mass for ',xcount,' particles !'
c Write lhe file
if(iout.eq.1)call lhesave(ievent)
......
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