IAP GITLAB

Commit 497eb2ed authored by Tanguy Pierog's avatar Tanguy Pierog

new strategy for mother/daughter for hepmc or lhe output :

* istmax=0 or NOT EPOS : 2 beam particles are created and all real particles are then saved (first particles which are coming from the primary vertex and then decay products)
* istmax=1 AND EPOS : 2 beam particles are created to be the mother of all projectile and target nucleons (in case of nuclei) which are the mother of virtual particles defining the particle production type. Stable and unstable particles are daughters of these virtual particles.
In both cases, the link is complete from beam to final particles for all models.
In case of ROOT output we can't use nuclei as beam particle but the mother/daughter link not being currently saved it doesn't matter.
Tested for QGSJETII and EPOS and OK.


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4185 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent 673a015c
......@@ -48,7 +48,7 @@ CRMC<OutputPolicy>::init()
fCfg.fTargetId,
fCfg.fHEModel,
fCfg.fProduceTables,
fCfg.fLHEout,
fCfg.fTypout,
fCfg.GetOutputFileName().c_str(),
fCfg.fParamFileName.c_str());
......@@ -78,7 +78,7 @@ CRMC<OutputPolicy>::run()
cout << " ==[crmc]==> Collision number " << iColl+1 << endl;
// loop over collisions
fInterface.crmc_generate(fCfg.fLHEout,iColl+1,
fInterface.crmc_generate(fCfg.fTypout,iColl+1,
gCRMC_data.fNParticles,
gCRMC_data.fImpactParameter,
gCRMC_data.fPartId[0],
......
......@@ -23,7 +23,7 @@ CRMCoptions::CRMCoptions(int argc, char** argv)
fTargetMomentum(-3500),
fParamFileName("crmc.param"),
fOutputFileName(""),
fLHEout(false),
fTypout(0),
fProduceTables(false),
fSeedProvided(false),
fTest(false)
......@@ -144,8 +144,11 @@ CRMCoptions::ParseOptions(int argc, char** argv)
}
}
if (fOutputMode == eLHE)
fLHEout = true;
if (fOutputMode == eLHE || fOutputMode == eLHEGZ)
fTypout = 1;
if (fOutputMode == eROOT)
fTypout = -1;
// parameter readout
if (opt.count("seed")) {
......@@ -215,7 +218,7 @@ const
{
if (pid < 10000) {
switch (fProjectileId) {
switch (pid) {
case 120 : return "pi"; break;
case -120 : return "antipi"; break;
case 1 : return "p"; break;
......
......@@ -44,12 +44,12 @@ class CRMCoptions {
int fProjectileId;
int fTargetId;
int fHEModel;
int fTypout;
double fProjectileMomentum;
double fTargetMomentum;
std::string fParamFileName;
std::string fOutputFileName;
bool fLHEout;
bool fProduceTables;
bool fSeedProvided;
bool fTest;
......
......@@ -2,7 +2,7 @@
c 15.01.2009 Simplified Main program and random number generator for epos
subroutine crmc_set_f(iEvent,iSeed,pproj,ptarg,
$ ipart,itarg,imodel,itab,ilheout,lheoutfile,param)
$ ipart,itarg,imodel,itab,itypout,lheoutfile,param)
***************************************************************
*
......@@ -16,7 +16,7 @@ c 15.01.2009 Simplified Main program and random number generator for epos
* itarg - target particle type
* imodel - HE model switch
* itab - force tables production or stop if missing
* ilheout - output type
* itypout - output type
* lheoutfile - output file name
* param - param file name
*
......@@ -25,10 +25,10 @@ c 15.01.2009 Simplified Main program and random number generator for epos
include "epos.inc"
c Input values
integer iSeed,ipart,itarg,imodel,itab,ilheout,iEvent,iout
integer iSeed,ipart,itarg,imodel,itab,itypout,iEvent,iout
double precision pproj, ptarg
character*1000 param,lheoutfile,output
common/lheoutput/iout,output
common/typoutput/iout,output
real m1,m2
double precision decms,e1,e2
......@@ -47,7 +47,7 @@ c Stop program if missing tables (after aaset)
if(itab.eq.1)producetables=.true.
c Set common for crmc_init
iout=ilheout
iout=itypout
output=lheoutfile
c Calculations of energy of the center-of-mass in the detector frame
......@@ -80,7 +80,7 @@ c Update some parameters value to run correctly
c The parameters can be changed optionnaly by reading a file
c (example.param) using the following subroutine call
call EposInput(param) !(it can be commented)
if ( model.ne.1 .and. model.ne.12 )istmax=0 !for most models daughter/mother link doesn't work
if ( model.ne.1 )istmax=0 !for most models virtual mothers are not defined
c if you put what is in input.optns in example.param, you can even run
c exactly the same way (coded parameters are overwritten). Don't forget
......@@ -97,7 +97,7 @@ c will not run.
implicit none
integer iout
character*1000 output
common/lheoutput/iout,output
common/typoutput/iout,output
c initialization for the given energy
call ainit
......@@ -150,7 +150,7 @@ c Fix final particles and some event parameters
call afinal
c Fill HEP common
if(model.ne.4.and.model.ne.5)call hepmcstore
if(model.ne.4.and.model.ne.5)call hepmcstore(iout) !use hepmcstore for non-HEP compatible events or non h-p models (Nucleus as projectile or target).
c optional Statistic information (only with debug level ish=1)
call astati
......@@ -375,7 +375,6 @@ c PDG
stop
endif
c Target definitions : for nucleons, idtarg does not exist
c Mass number matarg as well as charge, latarg, must be defined
......@@ -413,12 +412,8 @@ c PDG
endif
if ( model.eq.1 ) then !model variable has no eposLHC
istmax = 1 !final and mother particles (istmax=1 includes
!mother particles)
else
istmax=0
endif
istmax = 0 !final and mother particles (istmax=0 includes
!real mother particles)
end
......
......@@ -4,7 +4,7 @@
!switch fusion off !nuclear effects due to high density (QGP) in EPOS
!more realistic but slow (can be switched off)
!set istmax 1 !uncomment to get full daughter/mother particle chain with EPOS
!set istmax 1 !include virtual mother particles with EPOS to identify particle source
!set isigma 2 !uncomment to get correct inelastic cross-section for heavy ions with EPOS
......
......@@ -263,8 +263,7 @@ 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)
imaxhepstatus = MAX(1,MIN(ISTMAX+1,2))
IF(ISTHKK(k).GE.1 .AND. ISTHKK(k).LE.imaxhepstatus)THEN !! if final particle
IF(ISTHKK(k).GE.1 .AND. ISTHKK(k).LE.2)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...
......
......@@ -1078,10 +1078,11 @@ c count the number of particles to be stored (--> nptevt)
end
c-----------------------------------------------------------------------
subroutine hepmcstore
subroutine hepmcstore(iout)
c-----------------------------------------------------------------------
c writes the results of a simulation into the file with unit ifdt
c contains a description of the stored variables.
c iout < 0 : no nuclei as beam (not recognize by ROOT ???)
c-----------------------------------------------------------------------
include 'epos.inc'
double precision phep2,vhep2
......@@ -1091,6 +1092,8 @@ c-----------------------------------------------------------------------
&,jdahep2(2,nmxhep),phep2(5,nmxhep),vhep2(4,nmxhep)
integer iepo2hep(mxptl),ihep2epo(nmxhep),ihep2epo2(nmxhep)
double precision pprojin,ptargin
logical lrcore,lrcor0,lclean
c count the number of particles to be stored
......@@ -1126,17 +1129,33 @@ c jppevt ........ number of absolute projectile proton spectators
c jtnevt ........ number of absolute target neutron spectators
c jtpevt ........ number of absolute target proton spectators
c final particles (all mother/daughter + the one absorbed in core)
if(istmax.gt.0.and.ioclude.ne.0)then
istmaxhep=2
else
istmaxhep=min(istmax,1)
istmaxhep=1
endif
c fill first hep stack before ordering
nhep2=0
ipro=0
itar=0
pprojin=0d0
ptargin=0d0
do i=1,nptl
if(iout.ge.0)then
c initialize beam momenta
if(i.le.maproj)then
pprojin=pprojin+dble(pptl(3,i))
elseif(i.le.maproj+matarg)then
ptargin=ptargin+dble(pptl(3,i))
endif
endif
if(istptl(i).le.istmaxhep.or.i.le.maproj+matarg)then !store events with istptl < istmax
......@@ -1162,9 +1181,11 @@ c i ............. particle number
if(io.gt.0)then
if(istptl(io).le.1.and.i.gt.maproj+matarg.and..not.lclean)then!mother is normal particle (incl. spectators and fragments)
iadd=1
jm1hep=iepo2hep(io)
if(jorptl(i).gt.0)jm2hep=iepo2hep(jorptl(i))
elseif(istmaxhep.gt.0.and.(iorptl(io).gt.0.or.lclean))then
if(iepo2hep(io).gt.maproj+matarg)then
jm1hep=iepo2hep(io)
if(jorptl(i).gt.0)jm2hep=iepo2hep(jorptl(i))
endif
elseif(istmax.gt.0.and.(iorptl(io).gt.0.or.lclean))then
c create special father/mother to have the complete chain from beam to final particle
if(lclean)then
istptlio=99 !if cleaning defined: use all remnants
......@@ -1298,7 +1319,7 @@ c print *,'clean',i,iorptl(io),pptl(3,i),jm1io
goto 100
endif
if(nhep2+iadd.gt.nmxhep)then
if(nhep2+iadd.gt.nmxhep-2)then
print *,'Warning : produced number of particles is too high'
print *,' Particle list is truncated at, ', nmxhep
print *,' Skip event !'
......@@ -1343,7 +1364,7 @@ c pptl(5,i) ..... particle mass (GeV/c2)
phep2(5,nhep2)=dble(pptl(5,ii))
c istptl(i) ..... generation flag: last gen. (0) or not (1)
isthep2(nhep2)=min(2,istptl(ii)+1) !in hep:1=final, 2=decayed
if(i.le.maproj+matarg)isthep2(nhep2)=4 !beam particles
if(i.le.maproj+matarg.and.iout.lt.0)isthep2(nhep2)=4 !beam particles
c ityptl(i) ..... particle type (string, remnant ...)
c xorptl(1,i) ... x-component of formation point (fm)
vhep2(1,nhep2)=xorptl(1,ix)*1e-12 !conversion to mm
......@@ -1375,15 +1396,84 @@ c copy first list in final list to define daughters of beam particles
nhep=0
nhepio=0
if(istmaxhep.ne.0.and.model.eq.1)nhep=maproj+matarg
lrcor0=.true. !link spectator remnants to core only once
lrcore=.false.
c start with beam particles (except spectators producing fragments)
c define only 2 beam particles (projectile and target)
if(iout.ge.0)then
c projectile
c store initial target
if(maproj.gt.1)then
idprin=1000000000+maproj*10+laproj*10000
else
idprin=idprojin
endif
nhep=1
nhepio=nhepio+1
ii=1
id=idtrafo('nxs','pdg',idprin)
call idmass(idprin,amass)
idhep(nhep)=id
phep(1,nhep)=0d0
phep(2,nhep)=0d0
phep(3,nhep)=pprojin
phep(4,nhep)=sqrt(pprojin**2+dble(amass)**2)
phep(5,nhep)=dble(amass)
isthep(nhep)=4 !in hep:beam particle
vhep(1,nhep)=xorptl(1,ii)*1e-12 !conversion to mm
vhep(2,nhep)=xorptl(2,ii)*1e-12 !conversion to mm
vhep(3,nhep)=xorptl(3,ii)*1e-12 !conversion to mm
vhep(4,nhep)=xorptl(4,ii)*1E-12 !conversion to mm/c
jdahep(1,nhep)=3
jdahep(2,nhep)=maproj+matarg+2 !updated later if needed
jmohep(1,nhep)=-1
jmohep(2,nhep)=-1
c Target
if(matarg.gt.1)then
idtgin=1000000000+matarg*10+latarg*10000
else
idtgin=idtargin
endif
nhep=2
nhepio=nhepio+1
ii=maproj+1
id=idtrafo('nxs','pdg',idtgin)
call idmass(idtgin,amass)
idhep(nhep)=id
phep(1,nhep)=0d0
phep(2,nhep)=0d0
phep(3,nhep)=ptargin
phep(4,nhep)=sqrt(ptargin**2+dble(amass)**2)
phep(5,nhep)=dble(amass)
isthep(nhep)=4 !in hep:beam particle
vhep(1,nhep)=xorptl(1,ii)*1e-12 !conversion to mm
vhep(2,nhep)=xorptl(2,ii)*1e-12 !conversion to mm
vhep(3,nhep)=xorptl(3,ii)*1e-12 !conversion to mm
vhep(4,nhep)=xorptl(4,ii)*1E-12 !conversion to mm/c
jdahep(1,nhep)=3
jdahep(2,nhep)=maproj+matarg+2 !updated later if needed
jmohep(1,nhep)=-1
jmohep(2,nhep)=-1
endif
if(istmax.gt.0.and.model.eq.1)then
nhep=maproj+matarg !beam particles will be copied later
if(iout.ge.0)nhep=nhep+2 !add new beam particles
endif
nhep0=0
nskip=0
if(iout.ge.0)nskip=-2 !add 2 beam particles and may be skip nucleons
c reorder individual beam particles to make list of mothers for a given daughter
do j=1,maproj+matarg
if(istmaxhep.eq.0.or.model.ne.1)then
if(istmax.eq.0.or.model.ne.1)then
if(iout.lt.0)then
c when no daughter/mother informations, simply copy beam particles
nhep=nhep+1
......@@ -1407,9 +1497,16 @@ c when no daughter/mother informations, simply copy beam particles
isthep2(j)=-isthep2(j)
ihep2epo2(j)=-nhep
else
c skip individual beam particles in case of short list
nskip=nskip+1
endif
else
nhep0=nhep
nhep0=nhep !index of daughter list of current beam particle
nhepi0=nhepio+1
isthep(nhepi0)=0
nio=0
......@@ -1426,13 +1523,19 @@ c copy all daughters after the mother
phep(3,nhepio)=phep2(3,j)
phep(4,nhepio)=phep2(4,j)
phep(5,nhepio)=phep2(5,j)
isthep(nhepio)=4
if(iout.lt.0)then
isthep(nhepio)=4 !beam
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
else
isthep(nhepio)=3 !daughter of beam
jmohep(1,nhepio)=1
jmohep(2,nhepio)=2
endif
vhep(1,nhepio)=vhep2(1,j)
vhep(2,nhepio)=vhep2(2,j)
vhep(3,nhepio)=vhep2(3,j)
vhep(4,nhepio)=vhep2(4,j)
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
iepo2hep(j)=nhepio
isthep2(j)=-isthep2(j)
ihep2epo2(j)=-nhepio
......@@ -1457,13 +1560,19 @@ c copy all daughters after the mother
phep(3,nhepio)=phep2(3,i)
phep(4,nhepio)=phep2(4,i)
phep(5,nhepio)=phep2(5,i)
isthep(nhepio)=4
if(iout.lt.0)then
isthep(nhepio)=4 !beam
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
else
isthep(nhepio)=3 !daughter of beam
jmohep(1,nhepio)=1
jmohep(2,nhepio)=2
endif
vhep(1,nhepio)=vhep2(1,i)
vhep(2,nhepio)=vhep2(2,i)
vhep(3,nhepio)=vhep2(3,i)
vhep(4,nhepio)=vhep2(4,i)
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
isthep2(i)=-isthep2(i)
iepo2hep(i)=nhepio
ihep2epo2(i)=-nhepio
......@@ -1511,18 +1620,22 @@ c copy all daughters after the mother
enddo
iround=1
if(iout.ge.0.and.nhep0.eq.0)iround=0 !copy first all mothers and then daughters if not full list from EPOS (in that case nhep0.ne.0)
do iii=iround,1
c copy all other particles (secondary particles and spectators)
do k=maproj+matarg+1,nhep2
if(isthep2(k).gt.0)then
if(isthep2(k).gt.0.and.(iii.eq.1
& .or.(iii.eq.0.and.jmohep2(1,k).le.0)))then
c look for mother of fragments
nhep0=nhep+1
nhepi0=nhepio+1
if(jmohep2(1,k).gt.0.and.jmohep2(1,k).le.maproj+matarg)then
c copy all mothers before the daughter
c copy all mothers before the daughter to complete the beam remnant list
do j=1,maproj+matarg
if(isthep2(j).gt.0.and.jdahep2(1,j).eq.ihep2epo2(k))then
......@@ -1533,13 +1646,19 @@ c copy all mothers before the daughter
phep(3,nhepio)=phep2(3,j)
phep(4,nhepio)=phep2(4,j)
phep(5,nhepio)=phep2(5,j)
isthep(nhepio)=4
if(iout.lt.0)then
isthep(nhepio)=4 !beam
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
else
isthep(nhepio)=3 !daughter of beam
jmohep(1,nhepio)=1
jmohep(2,nhepio)=2
endif
vhep(1,nhepio)=vhep2(1,j)
vhep(2,nhepio)=vhep2(2,j)
vhep(3,nhepio)=vhep2(3,j)
vhep(4,nhepio)=vhep2(4,j)
jmohep(1,nhepio)=-1
jmohep(2,nhepio)=-1
jdahep(1,nhepio)=0
if(jdahep2(1,j).gt.0)jdahep(1,nhepio)=-jdahep2(1,j)
jdahep(2,nhepio)=0
......@@ -1571,36 +1690,46 @@ c copy all mothers before the daughter
if(jdahep2(2,k).gt.0)jdahep(2,nhep)=-ihep2epo2(jdahep2(2,k))
ihep2epo(nhep)=ihep2epo2(k)
if(ihep2epo(nhep).gt.0)iepo2hep(ihep2epo(nhep))=nhep
if(nhepio.lt.nhepi0)then
if(jmohep2(1,k).gt.0)then
if(ihep2epo2(jmohep2(1,k)).le.0)then
jmohep(1,nhep)=-ihep2epo2(jmohep2(1,k))
if(iii.eq.0)then !this particle is a first mother directly link to beam
jmohep(1,nhep)=1
jmohep(2,nhep)=2
jdahep(2,1)=nhep
jdahep(2,2)=nhep
else
if(nhepio.lt.nhepi0)then
if(jmohep2(1,k).gt.0)then
if(ihep2epo2(jmohep2(1,k)).le.0)then
jmohep(1,nhep)=-ihep2epo2(jmohep2(1,k))
else
jmohep(1,nhep)=iepo2hep(ihep2epo2(jmohep2(1,k)))
endif
else
jmohep(1,nhep)=iepo2hep(ihep2epo2(jmohep2(1,k)))
jmohep(1,nhep)=0
endif
else
jmohep(1,nhep)=0
endif
if(jmohep2(2,k).gt.0)then
if(ihep2epo2(jmohep2(2,k)).le.0)then
jmohep(2,nhep)=-ihep2epo2(jmohep2(2,k))
if(jmohep2(2,k).gt.0)then
if(ihep2epo2(jmohep2(2,k)).le.0)then
jmohep(2,nhep)=-ihep2epo2(jmohep2(2,k))
else
jmohep(2,nhep)=iepo2hep(ihep2epo2(jmohep2(2,k)))
endif
else
jmohep(2,nhep)=iepo2hep(ihep2epo2(jmohep2(2,k)))
jmohep(2,nhep)=0
endif
else
jmohep(2,nhep)=0
endif
else !for nuclear fragments
else !for nuclear fragments
jmohep(1,nhep)=nhepi0
jmohep(2,nhep)=nhepio
endif
endif
isthep2(k)=-isthep2(k)
isthep2(k)=-isthep2(k) !mark particle as copied to final list
endif
enddo
if(nhep.ne.nhep2.or.nhepio.ne.maproj+matarg)then
enddo !iround
if(nhep.ne.nhep2-nskip.or.nhepio.ne.maproj+matarg-nskip)then
print *,'Warning : number of particles changed after copy'
nrem1=0
do k=1,nhep2
......@@ -1615,16 +1744,18 @@ c copy all mothers before the daughter
print *,' ',nhep2,'->',nhep
nrem2=0
do k=1,nhep
if(isthep(k).eq.4)then
if((iout.ge.0.and.isthep(k).eq.3)
& .or.(iout.lt.0.and.isthep(k).eq.4))then
nrem2=nrem2+1
endif
c print *,' ',k,idhep(k),isthep(k),'from',ihep2epo(k)
print *,' ',k,idhep(k),isthep(k),'from',ihep2epo(k)
enddo
print *,' Particle list not consistent, skip event !'
print *,' ',nrem1,'->',nrem2
goto 10000
endif
c update daughter list with correct index
do j=1,nhep
......
......@@ -260,7 +260,7 @@ c LLIST is the code of final particle, P - its 4-momentum and mass.
$ ' PHOJET particle ',k,' id :',ic,' before conversion'
$ , ' momentum :',(sngl(PPHEP(i,k)),i=1,5)
IF(ISTHEPP(k).EQ.1)THEN
IF(ISTHEPP(k).GE.1.AND.ISTHEPP(k).LE.2)THEN
c final particle
nptl=nptl+1
nptlhep(k)=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