extend arrays to run central HI collisions

add sexaquark production


git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@7339 c7a5e08c-de06-0410-9364-b41cf42a0b17

CMakeLists.txt

@@ -12,7 +12,7 @@ OPTION (__PYTHIA__     "Build with model" OFF)
 OPTION (__HIJING__     "Build with model" OFF)
 OPTION (__SIBYLL__     "Build with model" ON)
 OPTION (__PHOJET__     "Build with model" OFF)
-OPTION (__DPMJET__     "Build with model" OFF)
+OPTION (__DPMJET__     "Build with model" ON)
 OPTION (__QGSJETII03__ "Build with model" OFF)
 OPTION (__QGSJETII04__ "Build with model" OFF)
 ######################################ONLY EDIT THIS######################################

README

-CRMC v1.6.0 Last modifications 2016/03/29
+CRMC v1.7.1 Last modifications 2019/10/25
 
 **********************************************************
                      The Program "crmc"
@@ -10,6 +10,8 @@ an easy-to-use command line interface. The output can be stored in
 different formats, i.e. in a root TTree or HepMC event file
 (http://lcgapp.cern.ch/project/simu/HepMC/).
 
+Version with extended arrays for heavy ion collisions.
+
 **********************************************************
                          Authors
 **********************************************************
@@ -222,7 +224,7 @@ The number of loops is hard-coded (20000 check models.F) and does not depend on
 
 The HepMC (2.06) IO library has a default limitation of the number of produced particles of 10000.
 Thus, events with more than that number of secondaries are either truncated or skipped. For heavy ion
-collisions this easily can become a problem in crmc. Please consider changing thie limitation to a
+collisions this easily can become a problem in crmc. Please consider changing this limitation to a
 higher value. We suggest 200000, which works for all models at LHC energies so far. We are happy for
 feedback on this issue if problems are encountered.
 The CRMC code automatically checks for the

src/epos.inc-HEPMC-TEMPLATE

@@ -238,8 +238,8 @@ c------------------------------------------------------------------------
       common/othe2/ifrade,iframe,idecay,jdecay,iremn
       integer       jframe,kframe
       common/othe3/jframe,kframe
-      integer      iselect
-      common/othe4/iselect
+      integer      iselect,iwoundg
+      common/othe4/iselect,iwoundg
       real         taumx,sigj
       integer      jpsi,jpsifi,nsttau,ijphis,ijtauan
       common/jpsif/jpsi,jpsifi,taumx,nsttau,sigj,ijphis,ijtauan

src/epos.incems

 
-      parameter (mmrye=4,mmryk=1)   !memory saving factor
+      parameter (mmrye=1,mmryk=1)   !memory saving factor
 
       parameter (kollmx=45000/mmrye,npommx=150)
       parameter (mamx=210/mmryk,ngrmx=1,ntymi=0,ntymx=2)

src/epos/epos-bas.f

@@ -10,7 +10,7 @@ c-----------------------------------------------------------------------
       include 'epos.incsem'
       common/record/maxrec(2),irecty(30,2)
       common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi/cifset/ifset
-      common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
+      common /ems12/bidiba,amhdibar,iodiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
       character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
 c     &,fndpmjet,fndpmjetpho,fndpmpath                  !qgs-II????????
 c      common/dpmjetfname/  fndpmjet,fndpmjetpho,fndpmpath
@@ -407,8 +407,8 @@ c ... up to here.
       alpdro(2)=0.3   !pt of leading droplet
       alpdro(3)=1.6  !alpha mass of leading droplet (iept=3)
 
-      iodiba=0.      ! if iodiba=1, study H-Dibaryon (not used (see ProRef in epos-ems ????)
-      bidiba=0.030   ! epsilon of H-Dibaryon
+      iodiba=0      ! if iodiba=1, study H-Dibaryon
+      bidiba=0.20   ! epsilon of H-Dibaryon
 
 c screening splitting +++
 
@@ -457,7 +457,7 @@ c  masses
       amhadr(7)=.548           !eta mass
       amhadr(8)=3.5            !J/psi mass
 c      qcmass=1.6               !c-quark mass  (in idmass = 1.2 )
-c      amhdibar=2.200           !h-dibaryon mass       !not used any more
+      amhdibar=1.80           !h-dibaryon mass       !not used any more
       qmass(0)=pmqq           !diquark effective bounding energy (for pt distribtions)
       isospin(0)=0
       call idmass(1,qumass)
@@ -767,6 +767,7 @@ c  other
       irdmpr=0   !random sign for projectile if 1
       ilprtg=1   !consider leading particle in projectile (1)
                  !or target  (-1) side
+      iwoundg=0  !save extended number of wounded nucleon in HEPEVT common
 
 c  constants
 
@@ -2996,6 +2997,7 @@ c-----------------------------------------------------------------------
       double precision rcproj,rctarg
       common/geom1/rcproj,rctarg
       common/photrans/phoele(4),ebeam
+      common /ems12/bidiba,amhdibar,iodiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
 
 
       external sptj
@@ -3343,6 +3345,7 @@ c          call ranfgt(seedp)   !not to change the seed ... not needed with 2 se
       endif
 
       if(iappl.ne.2)then
+        if(iodiba.gt.0)iospec=18
         call hnbspd(iospec)
         ktnbod=0
 c        if(model.eq.1)call hnbxxxini
@@ -3493,7 +3496,7 @@ c---------------------------------------------------------------------
       data nappl /0/
       common/record/maxrec(2),irecty(30,2)
       common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi
-      common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
+      common /ems12/bidiba,amhdibar,iodiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
       character*500 fndat,fnncs,fnIIdat,fnIIncs,fnII03dat,fnII03ncs
 c     &,fndpmjet,fndpmjetpho,fndpmpath
 c      common/dpmjetfname/  fndpmjet,fndpmjetpho,fndpmpath
@@ -4420,7 +4423,7 @@ c       unified parameters
       if(linex(ix:jx).eq.'sloreg')sloreg=sngl(val)
       if(linex(ix:jx).eq.'gamreg')gamreg=sngl(val)
       if(linex(ix:jx).eq.'r2reg' )r2reg= sngl(val)
-c      if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
+      if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
       if(linex(ix:jx).eq.'ptdiff')ptdiff=sngl(val)
       if(linex(ix:jx).eq.'alppar')alppar=sngl(val)
       if(linex(ix:jx).eq.'alpsea')alpsea=sngl(val)
@@ -4908,6 +4911,10 @@ c      if(line(i:j).eq.'off')iorsce=0
       call utword(line,i,j,0)
       if(line(i:j).eq.'on' ) ifrade=1
       if(line(i:j).eq.'off') ifrade=0
+        elseif(line(i:j).eq.'WoundedGlauber')then
+      call utword(line,i,j,0)
+      if(line(i:j).eq.'on' )iwoundg=1
+      if(line(i:j).eq.'off')iwoundg=0
         elseif(line(i:j).eq.'icocore')then
           stop'switch icocore not supported any more.    '
         endif

src/epos/epos-dro.f

@@ -51,11 +51,12 @@ c----------------------------------------------------------------------
       ncntmi=ncntmi+1
 
       if(ncntmi.eq.1)then
+        irescl=0
         call ranfgt(seedp)      !not to change the seed ...
         if(hydt.ne.'---')then
           call HydroTable2(0)
           call DefineParticles
-        else
+        elseif(ioclude.ne.3)then
           call ManiParticles(93,0)
         endif
         call ranfst(seedp)      ! ... after this initialization
@@ -67,6 +68,21 @@ c----------------------------------------------------------------------
 c      dez=0.5e-4
       call InitializeHyperbola !does not affect results but necessary
 
+      if(ioclude.eq.3)then
+        
+        nptl=1
+        idptl(nptl)=8*10**8
+        istptl(nptl)=1
+        ityptl(nptl)=11
+        pptl(1,nptl)=0.
+        pptl(2,nptl)=0.
+        pptl(3,nptl)=0.
+        pptl(4,nptl)=engy
+        pptl(5,nptl)=engy
+        call hnbaaa(1,iret)
+
+      else
+
 c  50  continue
       call GraCan
       energ=0
@@ -103,6 +119,8 @@ ccc      if(abs(energ-tecm).gt.0.1) goto50   !uncomment for energy conservation
         rinptl(n)=-9999
       enddo
 
+      endif
+
       call utprix('amicro',ish,ishini,4)
       return
       end

src/epos/epos-ems.f

@@ -5839,7 +5839,6 @@ c-----------------------------------------------------------------------
       logical gdrop, ghadr,gproj
       double precision ept(5),ep(4),aa(5),am2t,piq1,piq2,piq3
       common/emsptl/nppr(npommx,kollmx),npproj(mamx),nptarg(mamx)
-      common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
       character c*1,c1*1,c2*1
 
       call utpri('ProReF',ish,ishini,3)
@@ -6036,7 +6035,6 @@ c define masses
       xmdrmin=dble(fremnux(jcf)+amdrmin)**2
       xmdrmax=dble(fremnux(jcf)+amdrmax)**2
 
-
       if(ish.ge.4)write(ifch,*)'remnant masses:',am2t,amasini,amasmin
      &                ,xmdrmin,zz,iept
 
@@ -7088,7 +7086,7 @@ c-----------------------------------------------------------------
         double precision ep(5),a(5),re(5),p1(5)
         integer jc(nflav,2),jcv(nflav,2),jcini(nflav,2),jcvini(nflav,2)
      &    ,ifh(3),ic(2)
-        common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
+c        common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
         double precision ptm,qcm,u(3),utpcmd,ptt,phi,sxini,sxini0,strmas
      &                  ,ampt2dro,ampt2str,p5sq,amasex,drangen,xmaxrm
         logical gproj

src/epos/epos-ids.f

@@ -25,7 +25,7 @@ c     computes charge of particle with ident code id
 c     ichrg must be dimensioned nqlep+12
 c-----------------------------------------------------------------------
       dimension ichrg(53),ifl(3)
-      data ichrg/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,1,1,2,2*0
+      data ichrg/0,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,3,3,3,2*0
      *,2,-1,-1,2,-1,2,-1,2,0,0,0,-3,0,-3,0,-3,0,-3,3,0
      *,3,0,0,0,3,3,3,6,6,6,0/
       idabs=iabs(id)
@@ -315,7 +315,7 @@ c          quarks, leptons, etc
       ifl3=0
       jspin=0
       index=idabs
-      if(idabs.lt.20) return
+      if(idabs.lt.20.or.idabs.eq.30) return
 c          define index=20 for ks, index=21 for kl
       index=idabs+1
       if(id.eq.20) index=20
@@ -546,6 +546,7 @@ c-----------------------------------------------------------------------
 c     returns the mass of the particle with ident code id.
 c     (deuteron, triton and alpha mass come from Gheisha ???)
 c-----------------------------------------------------------------------
+      common /ems12/bidiba,amhdibar,iodiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
       dimension ammes0(15),ammes1(15),ambar0(30),ambar1(30)
       dimension amlep(52)
       parameter ( nqlep=41,nmes=2)
@@ -600,10 +601,10 @@ c     3/2+ baryon mass table
 c     entry
       id=idi
       amass=0.
-ctp060829      if(iabs(id).eq.30)then
-ctp060829        amass=amhdibar
-ctp060829        return
-ctp060829      endif
+      if(iabs(id).eq.30)then
+        amass=amhdibar
+        return
+      endif
       if(idi.eq.0)then
         id=1120                 !for air target
       elseif(abs(idi).ge.1000000000)then
@@ -1741,6 +1742,10 @@ c-----------------------------------------------------------------------
          ic(1)=222000
          ic(2)=000000
          return
+      elseif(id.eq.-30)then
+         ic(1)=000000
+         ic(2)=222000
+         return
       endif
       if(iabs(id).lt.1e8)then
         ix=iabs(id)/10
@@ -2430,9 +2435,12 @@ c      print *,'idtrafo',' ',code1,' ',code2,idi
       elseif(code2.eq.'pdg')then
         j=2
         ji=j
-        if(i.eq.1.and.id1.gt.1000000000)then !nucleus from NEXUS
+        if(i.eq.1.and.abs(id1).gt.1000000000)then !nucleus from NEXUS
           idtrafo=id1
           return
+        elseif(i.eq.1.and.abs(id1).eq.30)then !HDibaryon
+          idtrafo=sign(1020000020,id1)
+          return
         elseif(i.eq.4.and.id1.gt.402)then               !nucleus from Corsika
           idtrafo=1000000000+mod(id1,100)*10000+(id1/100)*10
           return

src/epos/epos-int.f

@@ -1909,7 +1909,7 @@ c----------------------------------------------------------------------
      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
       parameter(m1grid=10,kgrid=3,kegrid=3)
       parameter(m3xgrid=21*kgrid*kegrid)
-      parameter(mxcl=4000,mxcli=1000)
+      parameter(mxcl=10000,mxcli=1000)
       integer idropgrid(m1grid,m1grid,m3xgrid)
      &       ,jdropgrid(m1grid,m1grid,m3xgrid)
      &       ,jclu(m3xgrid),nsegp4(mxcl)
@@ -2114,8 +2114,14 @@ c...Start cluster formation
 
       nsegsuj=max(nsegce,nint(nsegsuj/1.08**ntry))
       ntry=ntry+1
-      if(ntry.gt.90)        !do not put more than 100 or change limit for p4mean
-     &call utstop('jintpo: cluster formation impossible ! &')
+c      if(ntry.gt.90)        !do not put more than 100 or change limit for p4mean
+c     &call utstop('jintpo: cluster formation impossible ! &')
+      if(ntry.gt.90)then        !do not put more than 100 or change limit for p4mean
+        write(*,*)'jintpo: cluster formation impossible ! redo ...'
+        write(*,*)'jintpo: (OK if rare)'
+        iret=1
+        goto 1000
+      endif
       nptl=nptla
 
       do 1 k=1,m3grid

src/epos/epos-uti.f

@@ -3675,6 +3675,9 @@ c-----------------------------------------------------------------------
 c-----------------------------------------------------------------------
       function fremnux(jc)
 c-----------------------------------------------------------------------
+      real bidiba,amhdibar
+      integer iodiba
+      common /ems12/bidiba,amhdibar,iodiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
       real pnll,ptq
       common/hadr1/pnll,ptq,exmass,cutmss,wproj,wtarg
       parameter (nflav=6)
@@ -3682,6 +3685,15 @@ c-----------------------------------------------------------------------
 c      ic(1)=210000
 c      ic(2)=0
 c      call iddeco(ic,jc)
+      if(iodiba.eq.1)then
+        if((jc(1,1).eq.2.and.jc(2,1).eq.2.and.jc(3,1).eq.2
+     .          .and.jc(1,2)+jc(2,2)+jc(3,2).eq.0)
+     . .or.(jc(1,2).eq.2.and.jc(2,2).eq.2.and.jc(3,2).eq.2
+     .       .and.jc(1,1)+jc(2,1)+jc(3,1).eq.0))then
+          fremnux=amhdibar-bidiba
+          return
+        endif
+      endif
       keu=jc(1,1)-jc(1,2)
       ked=jc(2,1)-jc(2,2)
       kes=jc(3,1)-jc(3,2)