IAP GITLAB

Commit 2051e085 authored by Tanguy Pierog's avatar Tanguy Pierog

correct hepmc output to have beam particles first

correct number of participants and spectators for all models
add number of hard collisions for EPOS
correct warning due to text size in utstop a alist(blist...)
correct flow transition in EPOS LHC (better pA collisions)


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@3326 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 250a4222
......@@ -54,6 +54,7 @@ CRMCdata() :
int npjevt;
int ntgevt;
int kolevt;
int kohevt;
int npnevt;
int ntnevt;
int nppevt;
......@@ -114,6 +115,7 @@ extern "C"
float zptevt; // ........ average Z-parton-targ
int minfra; // ........
int maxfra; // ........
int kohevt; // ........ number of inelastic hard collisions
} cevt_; //epos.inc
}
......
......@@ -98,6 +98,9 @@ OutputPolicyHepMC::FillEvent(const CRMCoptions& cfg, const int nEvent)
throw std::runtime_error("!!!Could not read next event");
}
// if(!HepMC::HEPEVT_Wrapper::check_hepevt_consistency(std::cout))
// HepMC::HEPEVT_Wrapper::print_hepevt(std::cout);
if (cfg.fTest){
// Test mode : compute directly some observables
if(!HepMC::HEPEVT_Wrapper::check_hepevt_consistency(std::cout))
......@@ -160,7 +163,7 @@ OutputPolicyHepMC::FillEvent(const CRMCoptions& cfg, const int nEvent)
// float sigma_inel_NN // nucleon-nucleon inelastic
// (including diffractive) cross-section
HepMC::HeavyIon ion(-1, //gCRMC_data_.koievt, // FIXME
HepMC::HeavyIon ion(gCRMC_data_.kohevt,
gCRMC_data.npjevt,
gCRMC_data.ntgevt,
gCRMC_data.kolevt,
......
......@@ -32,12 +32,12 @@ c---------------------------------------------------------------------------
real phievt,bimevt,pmxevt,egyevt
*,xbjevt,qsqevt,zppevt,zptevt
integer nevt,kolevt,koievt,npjevt
integer nevt,kolevt,koievt,kohevt,npjevt
*,ntgevt,npnevt,nppevt,ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt
*,nglevt,minfra,maxfra
common/cevt/phievt,nevt,bimevt,kolevt,koievt,pmxevt,egyevt,npjevt
*,ntgevt,npnevt,nppevt,ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt
*,xbjevt,qsqevt,nglevt,zppevt,zptevt,minfra,maxfra
*,xbjevt,qsqevt,nglevt,zppevt,zptevt,minfra,maxfra,kohevt
real rglevt,sglevt,eglevt,fglevt,typevt
integer ng1evt,ng2evt,ikoevt
common/c2evt/ng1evt,ng2evt,rglevt,sglevt,eglevt,fglevt,ikoevt
......@@ -48,6 +48,7 @@ c bimevt ........ absolute value of impact parameter
c phievt ........ angle of impact parameter
c kolevt ........ number of collisions
c koievt ........ number of inelastic collisions
c kohevt ........ number of hard collisions
c pmxevt ........ reference momentum
c egyevt ........ pp cm energy (hadron) or string energy (lepton)
c npjevt ........ number of primary projectile participants
......
This diff is collapsed.
......@@ -829,7 +829,7 @@ c print *,'proj'
zproj(i)=XA(i,3)
istptl(nptl)=1
iorptl(nptl)=-1
elseif(model.eq.7)then !QGSJetII
elseif(model.eq.7.or.model.eq.11)then !QGSJetII
xproj(i)=xan(i,1)
yproj(i)=xan(i,2)
zproj(i)=xan(i,3)
......@@ -861,7 +861,7 @@ c print *,'targ'
ztarg(i)=XB(i,3)
istptl(nptl)=1
iorptl(nptl)=-1
elseif(model.eq.7)then !QGSJetII
elseif(model.eq.7.or.model.eq.11)then !QGSJetII
xtarg(i)=xbn(i,1)
ytarg(i)=xbn(i,2)
ztarg(i)=xbn(i,3)
......
......@@ -11,7 +11,7 @@ c-----------------------------------------------------------------------
double precision omega,omlog,oma,omb,wab,wba,wmatrix,wzero,nbar
*,wzerox,rrr,eps,xprem,xmrem,om1intgck
parameter(eps=1.d-30)
common/col3/ncol,kolpt
common/col3/ncol,kolpt,ncolh
c logical modu
common/cems5/plc,s
double precision s,px,py,pomass,plc!,PhiExpo
......@@ -713,11 +713,14 @@ c --- Write ---
c --- Treat hard Pomeron
ncolh=0
do k=1,koll
ncolhp=0
do n=1,nprmx(k)
if(idpr(n,k).eq.3)then
if(ishpom.eq.1)then
call psahot(k,n,iret)
if(iret.eq.0)ncolhp=ncolhp+1
if(iret.eq.1)then
if(nbkpr(n,k).ne.0)then
nn=nbkpr(n,k)
......@@ -773,6 +776,7 @@ c --- Treat hard Pomeron
endif
endif
enddo
if(ncolhp.gt.0)ncolh=ncolh+1 !count hard binary collisions
enddo
if(iremn.ge.2)then
......@@ -2631,7 +2635,7 @@ c-------------------------------------------------------------------------
include 'epos.inc'
include 'epos.incems'
include 'epos.incsem'
common/col3/ncol,kolpt
common/col3/ncol,kolpt,ncolh
common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
double precision plc,s
common/cems5/plc,s
......@@ -4095,7 +4099,7 @@ c------------------------------------------------------------------------
integer icp(2),ict(2),ic(2),icp1(2),icp2(2),icm1(2),icm2(2)
integer jcp(nflav,2),jct(nflav,2),jcpv(nflav,2),jctv(nflav,2)
integer jcp1(nflav,2),jcp2(nflav,2),jcm1(nflav,2),jcm2(nflav,2)
common/col3/ncol,kolpt /cfacmss/facmss /cts/its
common/col3/ncol,kolpt,ncolh /cfacmss/facmss /cts/its
c entry
......@@ -8374,7 +8378,7 @@ c-----------------------------------------------------------------------
include 'epos.inc'
include 'epos.incems'
common/nucl3/phi,bimp
common/col3/ncol,kolpt
common/col3/ncol,kolpt,ncolh
integer kolpz(mamx),koltz(mamx)
call utpri('emszz ',ish,ishini,6)
......@@ -8400,15 +8404,17 @@ c determine ncol
if(ish.ge.7)write(ifch,*)'k,itpr,ncol,ncolx',k,itpr(k),ncol,ncolx
if(itpr(k).eq.0)goto 8
if(abs(itpr(k)).eq.1)ncoli=ncoli+1
ncol=ncol+1
i=iproj(k)
j=itarg(k)
istptl(i)=1
iorptl(i)=-1
tivptl(2,i)=coord(4,k)
istptl(maproj+j)=1
iorptl(maproj+j)=-1
tivptl(2,maproj+j)=coord(4,k)
ncol=ncol+1
if(itpr(k).ne.3)then !empty pair, remnant not modified
i=iproj(k)
j=itarg(k)
istptl(i)=1
iorptl(i)=-1
tivptl(2,i)=coord(4,k)
istptl(maproj+j)=1
iorptl(maproj+j)=-1
tivptl(2,maproj+j)=coord(4,k)
endif
8 continue
if(ncolx.ne.ncol)write(6,*)'ncolx,ncol:', ncolx,ncol
if(ncolx.ne.ncol)call utstop('********ncolx.ne.ncol********&')
......@@ -8450,6 +8456,7 @@ c ------------
phievt=phi
kolevt=ncol
koievt=ncoli
kohevt=ncolh
npjevt=npj
ntgevt=ntg
pmxevt=pnll
......
......@@ -51,7 +51,7 @@ c write(ifch,*)'droplet uds=',keu,ked,kes,' E=',pptl(5,ip)
!~~~~~redefine energy in case of radial flow
amin=utamnu(keu,ked,kes,kec,keb,ket,4) !utamnu(...,4) and not utamnu(...,5)
aumin=amuseg+yrmaxi !otherwise droplet from remnant decay
aumin=amuseg+yrmaxi !for rad and long flow
ipo=ip !could be too light after flow
if(ityptl(ip).eq.60)ipo=iorptl(ip)
tecmor=pptl(5,ipo)
......@@ -109,7 +109,7 @@ c write(ifch,*)'droplet uds=',keu,ked,kes,' E=',pptl(5,ip)
endif
!~~~~~redefine energy in case of long coll flow
! MANDATORY if RAD FLOW USED BECAUSE IT SMOOTH THE ETA DISTRIBUTION
! because of the grid structure for the cluster, flucutations in eta
! because of the grid structure for the cluster, fluctuations in eta
! appears if there is no smearing with long flow !
tecmx=tecm
c if(yco.gt.0..and.tecmor.gt.aumin) then
......@@ -131,8 +131,9 @@ c if(yco.gt.0..and.tecmor.gt.aumin) then
else
yco=0.
endif
c print*,'===== cluster energy: ',pptl(5,ip),tecmx,tecm,amin
c & ,delzet,yrmax,yco,ityptl(ip)
if(ish.ge.4)write(ifch,*)'===== cluster energy: '
& ,pptl(5,ip),tecmx,tecm,amin,aumin
& ,delzet,yrmax,yco,ityptl(ip)
!~~~~~~~~~define volume~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
......
This diff is collapsed.
......@@ -56,6 +56,7 @@ c *,xprh,xmrh
q2finsave=q2fin
zzzz=zparpro(kcol)+zpartar(kcol)
nnn=ltarg3(it)+lproj3(ip)
zz=1.+zoeinc*zzzz**2 !<-----
q2fin=q2fin*zz
......@@ -178,7 +179,7 @@ c preevolution part (momentum of soft string ends)
c print *,xmin,zopinc*fegypp*exp(-bpomr**2/b2xscr),xmax
c *,smin,max(xmin,min(0.99*xmax,zzzz*zzsoft))*ss
xmin=max(xmin,min((sqrt(s)-dble(amqpt+ammsdd))**2/s
& ,dble(max(2.*fegypp*exp(-bpomr**2/b2xscr),zzzz)*zopinc)))
& ,dble(max(nnn*fegypp*exp(-bpomr**2/b2xscr),zzzz)*zopinc)))
c xmin=max(xmin,1d0-exp(-dble(zzsoft*zzzz**2))) !???????????????
smin=xmin*s
endif
......
......@@ -3200,19 +3200,21 @@ c example: call utstop('error in subr xyz&')
c-----------------------------------------------------------------------
include 'epos.inc'
parameter(itext=40)
character text*40 ,txt*6
imax=itext+1
do i=itext,1,-1
if(text(i:i).eq.'&')imax=i
enddo
character text*(*) ,txt*6
imax=index(text,'&')
do 1 j=1,2
if(j.eq.1)then
ifi=ifch
else !if(j.eq.2)
ifi=ifmt
endif
if(imax.gt.1)then
write(ifi,101)('*',k=1,72),text(1:imax-1)
*,nrevt+1,nint(seedj),seedc,('*',k=1,72)
else
write(ifi,101)('*',k=1,72),' '
*,nrevt+1,nint(seedj),seedc,('*',k=1,72)
endif
101 format(
*1x,72a1
*/1x,'***** stop in ',a
......
SUBROUTINE FASTJETPPgenkt(P,NPART,R,PALG,F77JETS,NJETS)
DOUBLE PRECISION P(4,*), R, PALG, F, F77JETS(4,*)
INTEGER NPART, NJETS
DOUBLE PRECISION P(4,*), R, PALG, F77JETS(4,*)
DOUBLE PRECISION dP, dR, dPALG
INTEGER NPART, NJETS,Nd
F77JETS(1,1)=0d0
NJETS=0
dP=P(1,1)
Nd=NPART
dR=R
dPALG=PALG
c write(*,*)"FastJet called with :",F77JETS(1,1),NJETS,R,PALG
c stop"But FastJet not installed !"
......
......@@ -166,8 +166,8 @@ C
egymin=0.01
egymax=1e7
egymin=1.076 !min cms energy using pion as proj (min ekin=0.01)
egymax=1e4
irescl=0
......@@ -197,7 +197,7 @@ c-----------------------------------------------------------------------
fluincs=fflucrse()
if(ekin.lt.egymin)fluincs=0. !below egymin, no interaction
if(engy.lt.egymin)fluincs=0. !below egymin, no interaction
c if(iclpro.eq.3.and.ekin.lt.0.1d0)fluincs=0. !elastic event for kaons below 100MeV (problem) !
c initialize cross-sections
call xsigma
......@@ -275,7 +275,6 @@ c *,xbjevt,qsqevt,nglevt,zppevt,zptevt,minfra,maxfra
double precision P0(5),P1,EKIN1,PPROJ,TXX,TYY,TZZ,WEE,POO,PRES(4)
iret=0
ncol=0
np=0
b1=bminim
b2=min(bmax,bmaxim)
......@@ -325,17 +324,22 @@ c write(ifch,*)'target used',ibtar,ichtar,mmmat,ibres
if(NP.eq.0.and.NPHEAV.eq.0)goto 1001 !no interaction
ncol=1
nevt=1
kolevt=ncol
npjevt=maproj
ntgevt=matarg
pmxevt=pnll
egyevt=engy
bimevt=0.
bimp=0.
phievt=0.
phi=0.
kolevt=-1 !information not defined
koievt=-1
kohevt=-1
npjevt=maproj
ntgevt=matarg
npns=0 !counter to be defined later
npps=0
ntns=0
ntps=0
anintine=anintine+1.
......@@ -361,7 +365,6 @@ C NOW TREAT THE REMAINING TARGET NUCLEUS
goto 1001
endif
ekin1=EKRES/dble(IBRES)
PRES(4)=sqrt(PXRES*PXRES+PYRES*PYRES+PZRES*PZRES)
if(PRES(4).gt.0d0)then
PRES(3)=PZRES/PRES(4) !Pz
......@@ -376,13 +379,37 @@ C NOW TREAT THE REMAINING TARGET NUCLEUS
* ,' ... skip event.'
goto 1001
endif
ntps=ICRES
ntns=IBRES-ICRES
if(infragm.gt.0)then
ekin1=EKRES
c final particle is a nucleus
nptl=nptl+1
istptl(nptl)=0
if(nptl.gt.mxptl)call utstop('Fluka: mxptl too small&')
id=1000000000+ICRES*10000+IBRES*10 !code for nuclei
P0(5)=AAM(1)*ICRES+AAM(8)*(IBRES-ICRES) !Mass
P0(4)=ekin1+P0(5) !Energy
P1=sqrt((P0(4)+P0(5))*(P0(4)-P0(5))) !Momentum
P0(3)=PRES(3)*P1 !Pz
P0(2)=PRES(2)*P1 !Py
P0(1)=PRES(1)*P1 !Px
call flk2epo(id,P0,1)
nptl=maproj+matarg+1 !register spectators
do k=1,IBRES
else
ekin1=EKRES/dble(IBRES)
do k=1,IBRES
c final particle
nptl=nptl-1 !go backwards
c final particles are free nucleons
nptl=nptl+1
istptl(nptl)=0
if(nptl.gt.mxptl)call utstop('Fluka: mxptl too small&')
......@@ -402,11 +429,11 @@ c final particle
call flk2epo(id,P0,1)
enddo
ntgevt=matarg-IBRES !update number of wounded nucleons
endif
enddo
nptl=maproj+matarg
endif
endif
if(ish.ge.5)write(ifch,'(a,i5)')
$ ' number of particles from Fluka :',NP
......@@ -416,7 +443,7 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
ic=KPART(k)
icm=IPTOKP(ic) !convert to internal FLUKA ID to use AAM
if(ish.ge.7)write(ifck,'(a,i5,a,i5,2a,5(e10.4,1x))')
if(ish.ge.7)write(ifch,'(a,i5,a,i5,2a,5(e10.4,1x))')
$ ' FLUKA particle ',k,' id :',ic,' before conversion'
$ , ' Ekin, cos and mass :',TKI(k),CXR(k),CYR(k),CZR(k),AAM(icm)
......@@ -479,6 +506,15 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
PRES(2)=CYHEAV(J) !Py
PRES(1)=CXHEAV(J) !Px
c count spectators
if(PRES(3).gt.sqrt(pnll))then !projectile
npns=npns+ia-iz
npps=npps+iz
else !target
ntns=ntns+ia-iz
ntps=ntps+iz
endif
if(ia.eq.2.and.iz.eq.1)then !deuterium
id=17
P0(5)=AAM(-3) !Mass
......@@ -526,6 +562,50 @@ c final particle
enddo
c number of participants
if(laproj.gt.1)then
npjevt=maproj-npps-npns
npppar=laproj-npps
npnpar=npjevt-npppar
c set participant projectile and target as non spectators
do i=1,maproj
if(idptl(i).eq.1120)then
if(npppar.gt.0)then
npppar=npppar-1
else !restore spectators
iorptl(i)=0
endif
endif
if(idptl(i).eq.1220)then
if(npnpar.gt.0)then
npnpar=npnpar-1
else !restore spectators
iorptl(i)=0
endif
endif
enddo
endif
if(latarg.gt.1)then
ntgevt=matarg-ntps-ntns
ntppar=latarg-ntps
ntnpar=ntgevt-ntppar
do j=maproj+1,maproj+matarg
if(idptl(j).eq.1120)then
if(ntppar.gt.0)then
ntppar=ntppar-1
else !restore spectators
iorptl(j)=0
endif
endif
if(idptl(j).eq.1220)then
if(ntnpar.gt.0)then
ntnpar=ntnpar-1
else !restore spectators
iorptl(j)=0
endif
endif
enddo
endif
1000 if(ish.ge.5)write(ifch,'(a,i4,i3)')
$ ' End of Fluka event :',NP,iret
......@@ -560,8 +640,8 @@ c------------------------------------------------------------------------------
iorptl(nptl)=1
jorptl(nptl)=maproj+matarg
else
iorptl(nptl)=-1
jorptl(nptl)=0
iorptl(nptl)=maproj+1
jorptl(nptl)=maproj+matarg
endif
ifrptl(1,nptl)=0
ifrptl(2,nptl)=0
......
......@@ -127,7 +127,9 @@ c-----------------------------------------------------------------------
nptl=0
ncol=1
nevt=1
kolevt=ncol
kolevt=-1
koievt=-1
kohevt=-1
npjevt=maproj
ntgevt=matarg
pmxevt=pnll
......
......@@ -51,7 +51,10 @@
model=8
endif
if(cmodel(1:n).eq.'fluka' )then
stop'wrong choice!!!!' ! model=9
#ifndef __FLUKA__
stop'please compile with requested model'
#endif
model=9
endif
end
......@@ -112,7 +115,13 @@
call IniPHOJET
#endif
endif
c if(model.eq.9)call IniFluka
if(model.eq.9) then
#ifndef __FLUKA__
stop'please compile with requested model'
#else
call IniFluka
#endif
endif
end
subroutine IniEvtModel
......@@ -180,7 +189,13 @@ c if(model.eq.9)call IniFluka
call IniEvtPho
#endif
endif
c if(model.eq.9)call IniEvtFlu
if(model.eq.9)then
##ifndef __FLUKA__
stop'please compile with requested model'
#else
call IniEvtFlu
#endif
endif
end
subroutine emsaaaModel(model,id,iret)
......@@ -244,7 +259,13 @@ c if(model.eq.9)call IniEvtFlu
call emspho(iret)
#endif
endif
c if(model.eq.9)call emsflu(iret)
if(model.eq.9) then
#ifndef __FLUKA__
stop'please compile with requested model'
#else
call emsflu(iret)
#endif
endif
end
......@@ -341,11 +362,15 @@ c if(model.eq.9)call emsflu(iret)
sigi=sngl(GPROD)
sigc=sngl(GABS)
sige=sigt-sigi
c elseif(model.eq.9)then
c sigt=0.
c sige=0.
c sigc=0.
c sigi=flucrse(ekin,maproj,matarg,idtarg)
#endif
elseif(model.eq.9)then
#ifndef __FLUKA__
stop'please compile with requested model'
#else
sigt=0.
sige=0.
sigc=0.
sigi=flucrse(ekin,maproj,matarg,idtargin)
#endif
else
sigt=0.
......@@ -445,26 +470,29 @@ C cross sections
end
subroutine m9SIGMA(stot,sine,sela)
c include "epos.inc"
c double precision PLA,EKIN1,SHPTOT,SHPINE,ZZ,AA,Sel,Zl
stot=0.
sine=0.
sela=0.
c if(iclpro.eq.1)then
c IP=13
c elseif(iclpro.eq.2)then
c IP=1
c else
c IP=15
c endif
c PLA=dble(pnll)
c EKIN1=dble(ekin)
c ZZ=1
c AA=1
c stot=SHPTOT(IP,PLA)
c sine=SHPINE(IP,PLA)
c call SIGELH(IP,ZZ,AA,EKIN1,PLA,Sel,Zl)
c sela=sngl(Sel)
include "epos.inc"
double precision PLA,EKIN1,SHPTOT,SHPINE,ZZ,AA,Sel,Zl
#ifndef __FLUKA__
stop'please compile with requested model'
print *, stot,sine,sela
* ,PLA,EKIN1,SHPTOT,SHPINE,ZZ,AA,Sel,Zl !get rid of unused warning
#else
if(iclpro.eq.1)then
IP=13
elseif(iclpro.eq.2)then
IP=1
else
IP=15
endif
PLA=dble(pnll)
EKIN1=dble(ekin)
ZZ=1
AA=1
stot=SHPTOT(IP,PLA)
sine=SHPINE(IP,PLA)
call SIGELH(IP,ZZ,AA,EKIN1,PLA,Sel,Zl)
sela=sngl(Sel)
#endif
end
subroutine decaymod(ip,iret)
......
......@@ -80,7 +80,6 @@ c-----------------------------------------------------------------------
double precision XA(210,3),XB(210,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
common/nucl3/phi,bimp
common/col3/ncol,kolpt
common /q_area12/ nsp
common /q_area14/ esp(4,95000),ich(95000)
double precision esp
......@@ -89,7 +88,6 @@ c IAF(i) - mass of the i-th fragment
COMMON /Q_AREA13/ NSF,IAF(210)
COMMON /Q_AREA99/ NWT
ncol=0
iret=0
b1=bminim
b2=min(bmax,bmaxim)
......@@ -102,9 +100,10 @@ c IAF(i) - mass of the i-th fragment
nsp=0
call psconf
ncol=1
nevt=1
kolevt=ncol
kolevt=-1
koievt=-1
kohevt=-1
npjevt=maproj