EPOS/epos167/eposurqmd.f

   1 c----------------------------------------------------------------------------------
   2 c   urqmdinit ... transfer EPOS -> UrQMD
   3 c   urqmdepos ... main uUrQMD routine
   4 c                    essentially copied from urqmdepos (urqmd.f)
   5 c   urqmdexit ... transfer UrQMD -> EPOS
   6 c   file14outx .. analysis, essentially copied from file14out (output.f)
   7 c----------------------------------------------------------------------------------
   8 c modifs of urqmd routines:
   9 c   upmerge.f:
  10 c      comment the line "write(6,*) '(Info) pdg2ityp: ..."
  11 c----------------------------------------------------------------------------------
  12
  13 c----------------------------------------------------------------------------------
  14       subroutine urqmdinit
  15 c----------------------------------------------------------------------------------
  16
  17       implicit none
  18
  19
  20
  21       !urqmd common block
  22
  23 c     debug and validity range
  24
  25       integer nmax, nspl
  26       real*8 hit_sphere
  27       parameter (nmax = 40000) ! maximum number of particles
  28       parameter (nspl = 500)  ! dimension of spline arrays
  29       parameter (hit_sphere = 8.d0)  ! hard collision cutoff: 251 mbarn
  30
  31       integer Ap, At, Zp, Zt, npart, nbar, nmes, ctag
  32       integer nsteps,ranseed,event,eos,dectag,uid_cnt
  33       integer NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
  34       real*8  time,  acttime, bdist, ebeam, bimp,bmin,ecm
  35 c 7 integer
  36
  37       common /sys/ npart, nbar, nmes, ctag,nsteps,uid_cnt,
  38      +             ranseed,event,Ap,At,Zp,Zt,eos,dectag,
  39      +             NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
  40       common /rsys/ time,acttime,bdist,bimp,bmin,ebeam,ecm
  41
  42
  43
  44       real*8
  45      +     gw, sgw, delr, fdel, dt,
  46      +     da, db,
  47      +     Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky, gamYuk, drPau, dpPau,
  48      +     dtimestep
  49 c 19 real*8
  50
  51       real*8 cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
  52       common /cuts/ cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
  53
  54       real*8 spx(nspl), spPauy(nspl), outPau(nspl),
  55      +                spCby(nspl),  outCb(nspl),
  56      +                spYuky(nspl), outYuk(nspl),
  57      +                spSkyy(nspl), outSky(nspl),
  58      +                spdwwy(nspl), outdww(nspl)
  59
  60       common /spdata/ spx, spPauy, outPau, spCby,  outCb,
  61      +                     spYuky, outYuk, spSkyy, outSky,
  62      +                     spdwwy, outdww
  63
  64       real*8
  65      +     r0(nmax), rx(nmax), ry(nmax), rz(nmax),
  66      +     p0(nmax), px(nmax), py(nmax), pz(nmax),
  67      +     airx(nmax), airy(nmax), airz(nmax),
  68      +     aipx(nmax), aipy(nmax), aipz(nmax),
  69      +     aorx(nmax,4), aory(nmax,4), aorz(nmax,4),
  70      +     aopx(nmax,4), aopy(nmax,4), aopz(nmax,4),
  71      +     fmass(nmax), rww(nmax),
  72      +     dectime(nmax), tform(nmax), xtotfac(nmax)
  73
  74
  75       integer spin(nmax),ncoll(nmax),charge(nmax),strid(nmax),
  76      +        ityp(nmax),lstcoll(nmax),iso3(nmax),origin(nmax),uid(nmax)
  77       common/isys/spin,ncoll,charge,ityp,lstcoll,iso3,origin,strid,
  78      +            uid
  79
  80       common /coor/ r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime
  81       common /frag/ tform, xtotfac
  82
  83
  84       common /aios/ airx, airy, airz, aipx, aipy, aipz,
  85      +              aorx, aory, aorz, aopx, aopy, aopz
  86
  87       common /pots/ Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky,
  88      +              gamYuk, drPau, dpPau, gw, sgw, delr, fdel,
  89      +              dt,da, db,dtimestep
  90
  91
  92 c spectator arrays
  93         integer smax
  94         parameter(smax=500)  ! maximum number of spectators
  95         real*8 r0s(smax), rxs(smax), rys(smax), rzs(smax),
  96      +         p0s(smax), pxs(smax), pys(smax), pzs(smax),
  97      +         sfmass(smax)
  98
  99
 100         integer sspin(smax), scharge(smax), sityp(smax), siso3(smax),
 101      +          suid(smax)
 102
 103         integer nspec
 104
 105         common /scoor/ r0s, rxs, rys, rzs, p0s, pxs ,pys, pzs, sfmass
 106
 107         common /sisys/ sspin, scharge, sityp, siso3, suid
 108
 109         common /ssys/ nspec
 110
 111
 112         real*8 p0td(2,nmax),pxtd(2,nmax),pytd(2,nmax),pztd(2,nmax),
 113      +         fmasstd(2,nmax)
 114         integer ityptd(2,nmax),iso3td(2,nmax)
 115         integer itypt(2),uidt(2),origint(2),iso3t(2)
 116
 117         common /rtdelay/p0td,pxtd,pytd,pztd,fmasstd
 118         common /itdelay/ityptd,iso3td
 119         common /svinfo/itypt,uidt,origint,iso3t
 120         real*8 ffermpx(nmax), ffermpy(nmax), ffermpz(nmax)
 121         real*8 peq1, peq2
 122         common /ffermi/ ffermpx, ffermpy, ffermpz
 123         common /peq/ peq1,peq2
 124
 125
 126         integer numcto,numctp,maxstables
 127         parameter(numcto=400) ! maximum number of options
 128         parameter(numctp=400) ! maximum number of parameters
 129         parameter(maxstables=20) ! maximum number of stable particles
 130         integer   CTOption(numcto)
 131         real*8    CTParam(numctp)
 132         common /options/CTOption,CTParam
 133
 134         real*8 frr0(nmax), frrx(nmax), frry(nmax), frrz(nmax),
 135      +     frp0(nmax), frpx(nmax), frpy(nmax), frpz(nmax)
 136
 137       common /frcoor/ frr0, frrx, frry, frrz, frp0, frpx, frpy, frpz
 138               integer maxbar,maxbra,minbar
 139       integer offmeson,maxmeson,pimeson,maxbrm,minnuc,mindel
 140       integer maxbrs1,maxbrs2
 141       integer numnuc,numdel,nucleon,maxnuc,maxdel
 142       integer minmes,maxmes
 143
 144
 145       parameter (minnuc=1) ! lowest baryon particle ID
 146       parameter (minmes=100) ! lowest meson particle ID
 147       parameter (maxmes=132) ! hightest meson particle ID
 148
 149 c number of resonances of a kind
 150       parameter (numnuc=16) ! number of nucleon resonances
 151       parameter (numdel=10) ! number of delta resonances
 152 c indices of minimal and maximal itype of a kind (redundant but nice)
 153       parameter (maxnuc=minnuc+numnuc-1) ! highest nucleon ID
 154       parameter (mindel=minnuc+maxnuc)   ! lowest delta ID
 155       parameter (maxdel=mindel+numdel-1) ! highest delta ID
 156
 157 c minres & maxres define the range of nonstable & nonstrange baryons
 158       integer minres,maxres
 159       parameter (minres=minnuc+1) ! lowest baryon resonance ID
 160       parameter (maxres=maxdel)   ! highest (nonstrange) baryon
 161                                   ! resonance ID
 162
 163 c strangenes.ne.0 baryon resonances
 164       integer minlam,minsig,mincas,minome
 165       integer numlam,numsig,numcas,numome
 166       integer maxlam,maxsig,maxcas,maxome
 167       parameter (numlam=13) ! number of lambda states
 168       parameter (numsig=9)  ! number of sigma states
 169       parameter (numcas=6)  ! number of cascade states
 170       parameter (numome=1)  ! number of omega states
 171       parameter (minlam=mindel+numdel)   ! ID of lowest lambda state
 172       parameter (maxlam=minlam+numlam-1) ! ID of highest lambda state
 173       parameter (minsig=minlam+numlam)   ! ID of lowest sigma state
 174       parameter (maxsig=minsig+numsig-1) ! ID of highest sigma state
 175       parameter (mincas=minsig+numsig)   ! ID of lowest cascade state
 176       parameter (maxcas=mincas+numcas-1) ! ID of highest cascade state
 177       parameter (minome=mincas+numcas)   ! ID of lowest omega state
 178       parameter (maxome=minome+numome-1) ! ID of highest omega state
 179
 180 c minbar & maxbar define the range of all baryons
 181       parameter (minbar=minnuc) ! ID of lowest baryon state
 182       parameter (maxbar=maxome) ! ID of highest baryon state
 183
 184       parameter (offmeson=minmes) ! offset between zero and lowest
 185                                   ! meson state
 186       parameter (maxmeson=maxmes) ! ID of highest meson state
 187 c... these variables are in principal obsolete and should be exchanged
 188 c were referenced
 189
 190 c... avoid hard coded itypes
 191       integer itrho,itome,iteta,itkaon,itphi,itetapr
 192       parameter (itkaon=106)   ! ID of kaon
 193       parameter (itrho=104)    ! ID of rho meson
 194       parameter (itome=103)    ! ID of omega meson
 195       parameter (iteta=102)    ! ID of eta
 196       parameter (itphi=109)    ! ID of phi
 197       parameter (itetapr=107)  ! ID of eta'
 198       parameter (pimeson=101)  ! ID of $\pi$
 199       parameter (nucleon=minnuc) ! ID of nucleon
 200
 201       integer itmin,itmax
 202       parameter (itmin=minnuc)  ! lowest defined ID
 203       parameter (itmax=maxmes)  ! highest defined ID
 204 c
 205       parameter (maxbra=11)  ! decay channels for $s=0$ baryon resonances
 206       parameter (maxbrm=25) ! decay channels for meson resonances
 207       parameter (maxbrs1=10)! decay channels for $s=1$ baryon resonances
 208       parameter (maxbrs2=3) ! decay channels for $s=2$ baryon resonances
 209
 210 c
 211        integer mlt2it(maxmes-minmes) ! meson IDs sorted by multipletts
 212
 213
 214       real*8 massoff,mresmin,mresmax
 215       parameter (massoff=1d-4)      ! offset for mass generation
 216       parameter (mresmin=1.0765d0)  ! minimum baryon resonance mass
 217       parameter (mresmax=5d0)       ! maximum baryon resonance mass
 218
 219       character*45 versiontag
 220       common /versioning/ versiontag
 221
 222       real*8 massres(minbar:maxbar),widres(minbar:maxbar)
 223       real*8 branmes(0:maxbrm,minmes+1:maxmes)
 224       real*8 branres(0:maxbra,minnuc+1:maxdel)
 225       real*8 branbs1(0:maxbrs1,minlam+1:maxsig)
 226       real*8 branbs2(0:maxbrs2,mincas+1:maxcas)
 227       integer Jres(minbar:maxbar)
 228       integer Jmes(minmes:maxmes)
 229       integer pares(minbar:maxbar),pames(minmes:maxmes)
 230       integer Isores(minbar:maxbar), Isomes(minmes:maxmes)
 231       integer brtype(4,0:maxbra),bmtype(4,0:maxbrm)
 232       integer bs1type(4,0:maxbrs1),bs2type(4,0:maxbrs2)
 233       real*8 massmes(minmes:maxmes)
 234       real*8 mmesmn(minmes:maxmes)
 235       real*8 widmes(minmes:maxmes)
 236       integer strres(minbar:maxbar),strmes(minmes:maxmes)
 237
 238       integer lbr(0:maxbra,minnuc+1:maxdel)
 239       integer lbs1(0:maxbrs1,minlam+1:maxsig)
 240       integer lbs2(0:maxbrs2,mincas+1:maxcas)
 241       integer lbm(0:maxbrm,minmes+1:maxmes)
 242
 243       common /resonances/ massres,widres,massmes,widmes,mmesmn,
 244      ,                    branres,branmes,branbs1,branbs2,
 245      ,                    bs1type,bs2type,lbs1,lbs2,lbm,
 246      ,                    jres,jmes,lbr,brtype,pares,pames,
 247      ,                    bmtype,
 248      ,                    Isores,Isomes,strres,strmes,mlt2it
 249
 250
 251
 252
 253
 254
 255       !#############################################################################
 256       !epos common blocks for particle list
 257
 258       integer mmry,mxptl
 259       parameter (mmry=1)   !memory saving factor
 260       parameter (mxptl=200000/mmry) !max nr of particles in epos ptl list
 261
 262       integer iorptl(mxptl),idptl(mxptl),istptl(mxptl),
 263      *  ifrptl(2,mxptl),jorptl(mxptl),ibptl(4,mxptl),ityptl(mxptl)
 264       real pptl(5,mxptl),tivptl(2,mxptl),xorptl(4,mxptl)
 265
 266       common/cptl/nptl,pptl,iorptl,idptl
 267      *,istptl,tivptl,ifrptl,jorptl
 268      *,xorptl,ibptl,ityptl
 269
 270       integer nptl
 271
 272       integer nevt,kolevt,koievt,npjevt,ntgevt,npnevt,nppevt,
 273      *        ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt,
 274      *        nglevt,minfra,maxfra
 275       real phievt,bimevt,pmxevt,egyevt,xbjevt,qsqevt,zppevt,zptevt
 276
 277       common/cevt/phievt,nevt,bimevt,kolevt,koievt,pmxevt,egyevt,npjevt
 278      *,ntgevt,npnevt,nppevt,ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt
 279      *,xbjevt,qsqevt,nglevt,zppevt,zptevt,minfra,maxfra
 280
 281       integer laproj,maproj,latarg,matarg,nptlini,nptlbd
 282       real core,fctrmx
 283       common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
 284       common/c4ptl/nptlbd
 285
 286       integer istore,istmax,irescl,ntrymx,nclean,iopdg
 287       real gaumx
 288       common/othe1/istore,istmax,gaumx,irescl,ntrymx,nclean,iopdg
 289       !#############################################################################
 290
 291       integer i,nn,iok,idtmp,ityptmp,iso3tmp,itmp
 292
 293
 294
 295       integer idtrafo
 296       external idtrafo
 297       integer fchg
 298       external fchg
 299       real*8 dectim
 300       external dectim
 301
 302       real*8 mintime,eb
 303       integer j,k,icount,npold
 304       integer strcount
 305       common /inewpartxx/ strcount
 306
 307       call uinit(0)
 308       call osc_header
 309       call osc99_header
 310
 311       npart = 0
 312       npold = 0
 313       nbar=0
 314       nmes=0
 315       uid_cnt=0
 316 c reset counters
 317 c     all collisions/decays
 318       ctag  = 0
 319 c     all decays
 320       dectag = 0
 321 c     number of prod. hard resonances
 322       NHardRes=0
 323 c     number of prod. soft resonances
 324       NSoftRes=0
 325 c     number of prod. resonances via decay
 326       NDecRes=0
 327 c     number of blocked collisions
 328       NBlColl=0
 329 c     number of elastic collisions
 330       NElColl=0
 331 c     number of strings
 332       strcount=1
 333 c
 334       eb=0D0
 335 c icount is the number of EXTRAordinary pro/tar combinations (i.e. pion ...)
 336       icount = 0
 337 c reset particle vectors
 338       do 20 j=1,nmax
 339         spin(j)  = 0
 340         ncoll(j) = 0
 341         lstcoll(j)=0
 342         r0(j) = 0.0
 343         rx(j)    = 0.0
 344         ry(j)    = 0.0
 345         rz(j)    = 0.0
 346         p0(j)    = 0.0
 347         px(j)    = 0.0
 348         py(j)    = 0.0
 349         pz(j)    = 0.0
 350         frr0(j) = 0.0
 351         frrx(j)    = 0.0
 352         frry(j)    = 0.0
 353         frrz(j)    = 0.0
 354         frp0(j)    = 0.0
 355         frpx(j)    = 0.0
 356         frpy(j)    = 0.0
 357         frpz(j)    = 0.0
 358         fmass(j) = 0.0
 359         charge(j)= 0
 360         iso3(j)  = 0
 361         ityp(j)  = 0
 362         dectime(j)= 0.0
 363         origin(j)=0
 364         tform(j)=0.0
 365         xtotfac(j)=1.0
 366         strid(j)=0
 367         uid(j)=0
 368          ffermpx(j) = 0.0
 369          ffermpy(j) = 0.0
 370          ffermpz(j) = 0.0
 371
 372          do 21 k=1,2
 373             p0td(k,j)=0.d0
 374             pxtd(k,j)=0.d0
 375             pytd(k,j)=0.d0
 376             pztd(k,j)=0.d0
 377             fmasstd(k,j)=0.d0
 378             ityptd(k,j)=0
 379            iso3td(k,j)=0
 380  21      continue
 381  20   continue
 382
 383
 384 c epos event info to urqmd event info
 385       bimp = bimevt
 386
 387 c convert epos list to urqmd list here. or in a seperate routine?
 388       npart=0  !number of output particles
 389       mintime = 1d2 !the minimum formation time
 390       nptlini=maproj+matarg+1
 391
 392 c fill in the baryons first
 393       nbar = 0
 394       do nn=nptlini,nptlbd
 395         iok=1
 396         if(istptl(nn).gt.istmax)iok=0
 397         if(iok.eq.1)then
 398            idtmp=idtrafo('nxs','pdg',idptl(nn))
 399            call pdg2id(ityptmp,iso3tmp,idtmp)
 400            if(abs(ityptmp).le.maxbar) then
 401                nbar=nbar+1
 402                r0(nbar)=xorptl(4,nn)
 403                rx(nbar)=xorptl(1,nn)
 404                ry(nbar)=xorptl(2,nn)
 405                rz(nbar)=xorptl(3,nn)
 406                p0(nbar)=pptl(4,nn)
 407                px(nbar)=pptl(1,nn)
 408                py(nbar)=pptl(2,nn)
 409                pz(nbar)=pptl(3,nn)
 410                fmass(nbar)=pptl(5,nn)
 411                ityp(nbar)=ityptmp
 412                iso3(nbar)=iso3tmp
 413                charge(nbar)=fchg(iso3(nbar),ityp(nbar))
 414                lstcoll(nbar)=0
 415                ncoll(nbar)=0
 416                origin(nbar)=0
 417                tform(nbar)=r0(nbar)
 418                dectime(nbar)=dectim(nbar,1)+tform(nbar)
 419                xtotfac(nbar)=0d0
 420                if(r0(nbar).lt.mintime) mintime = r0(nbar)
 421            endif
 422         endif
 423       enddo
 424
 425 c then fill in the mesons
 426       nmes = 0
 427       do nn=nptlini,nptlbd
 428         iok=1
 429         if(istptl(nn).gt.istmax)iok=0
 430         if(iok.eq.1)then
 431            idtmp=idtrafo('nxs','pdg',idptl(nn))
 432            call pdg2id(ityptmp,iso3tmp,idtmp)
 433            if(abs(ityptmp).ge.minmes) then
 434                nmes=nmes+1
 435                itmp=nbar+nmes
 436                r0(itmp)=xorptl(4,nn)
 437                rx(itmp)=xorptl(1,nn)
 438                ry(itmp)=xorptl(2,nn)
 439                rz(itmp)=xorptl(3,nn)
 440                p0(itmp)=pptl(4,nn)
 441                px(itmp)=pptl(1,nn)
 442                py(itmp)=pptl(2,nn)
 443                pz(itmp)=pptl(3,nn)
 444                fmass(itmp)=pptl(5,nn)
 445                ityp(itmp)=ityptmp
 446                iso3(itmp)=iso3tmp
 447                charge(itmp)=fchg(iso3(itmp),ityp(itmp))
 448                lstcoll(itmp)=0
 449                ncoll(itmp)=0
 450                origin(itmp)=0
 451                tform(itmp)=r0(itmp)
 452                dectime(itmp)=dectim(itmp,1)+tform(itmp)
 453                xtotfac(itmp)=0d0
 454                if(r0(itmp).lt.mintime) mintime = r0(itmp)
 455            endif
 456         endif
 457       enddo
 458
 459       npart = nbar + nmes
 460
 461
 462 c back to the same starting time
 463       do i = 1, npart
 464          !save freeze-out configuration, in case of no further
 465          !rescatterings
 466          frr0(i) = r0(i)
 467          frrx(i) = rx(i)
 468          frry(i) = ry(i)
 469          frrz(i) = rz(i)
 470          frp0(i) = p0(i)
 471          frpx(i) = px(i)
 472          frpy(i) = py(i)
 473          frpz(i) = pz(i)
 474          rx(i)=rx(i)-px(i)/p0(i)*(r0(i)-mintime)
 475          ry(i)=ry(i)-py(i)/p0(i)*(r0(i)-mintime)
 476          rz(i)=rz(i)-pz(i)/p0(i)*(r0(i)-mintime)
 477          r0(i)=mintime
 478       enddo
 479
 480       acttime=mintime
 481
 482 c      write(*,*)'DEBUG INFO (epos.f): ',mintime,npart,istmax,nbar,nmes
 483
 484       return
 485       end
 486
 487
 488 c----------------------------------------------------------------------------------
 489       subroutine urqmdexit
 490 c----------------------------------------------------------------------------------
 491 c  transfer UrQMD -> EPOS
 492 c----------------------------------------------------------------------------------
 493
 494       implicit none
 495
 496
 497 !urqmd common block
 498
 499 c     debug and validity range
 500
 501       integer nmax, nspl
 502       real*8 hit_sphere
 503       parameter (nmax = 40000) ! maximum number of particles
 504       parameter (nspl = 500)  ! dimension of spline arrays
 505       parameter (hit_sphere = 8.d0)  ! hard collision cutoff: 251 mbarn
 506
 507       integer Ap, At, Zp, Zt, npart, nbar, nmes, ctag
 508       integer nsteps,ranseed,event,eos,dectag,uid_cnt
 509       integer NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
 510       real*8  time,  acttime, bdist, ebeam, bimp,bmin,ecm
 511 c 7 integer
 512
 513       common /sys/ npart, nbar, nmes, ctag,nsteps,uid_cnt,
 514      +             ranseed,event,Ap,At,Zp,Zt,eos,dectag,
 515      +             NHardRes,NSoftRes,NDecRes,NElColl,NBlColl
 516       common /rsys/ time,acttime,bdist,bimp,bmin,ebeam,ecm
 517
 518
 519
 520       real*8
 521      +     gw, sgw, delr, fdel, dt,
 522      +     da, db,
 523      +     Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky, gamYuk, drPau, dpPau,
 524      +     dtimestep
 525 c 19 real*8
 526
 527       real*8 cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
 528       common /cuts/ cutmax, cutPau, cutCb, cutYuk, cutSky, cutdww
 529
 530       real*8 spx(nspl), spPauy(nspl), outPau(nspl),
 531      +                spCby(nspl),  outCb(nspl),
 532      +                spYuky(nspl), outYuk(nspl),
 533      +                spSkyy(nspl), outSky(nspl),
 534      +                spdwwy(nspl), outdww(nspl)
 535
 536       common /spdata/ spx, spPauy, outPau, spCby,  outCb,
 537      +                     spYuky, outYuk, spSkyy, outSky,
 538      +                     spdwwy, outdww
 539
 540       real*8
 541      +     r0(nmax), rx(nmax), ry(nmax), rz(nmax),
 542      +     p0(nmax), px(nmax), py(nmax), pz(nmax),
 543      +     airx(nmax), airy(nmax), airz(nmax),
 544      +     aipx(nmax), aipy(nmax), aipz(nmax),
 545      +     aorx(nmax,4), aory(nmax,4), aorz(nmax,4),
 546      +     aopx(nmax,4), aopy(nmax,4), aopz(nmax,4),
 547      +     fmass(nmax), rww(nmax),
 548      +     dectime(nmax), tform(nmax), xtotfac(nmax)
 549
 550
 551       integer spin(nmax),ncoll(nmax),charge(nmax),strid(nmax),
 552      +        ityp(nmax),lstcoll(nmax),iso3(nmax),origin(nmax),uid(nmax)
 553       common/isys/spin,ncoll,charge,ityp,lstcoll,iso3,origin,strid,
 554      +            uid
 555
 556       common /coor/ r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime
 557       common /frag/ tform, xtotfac
 558
 559
 560       common /aios/ airx, airy, airz, aipx, aipy, aipz,
 561      +              aorx, aory, aorz, aopx, aopy, aopz
 562
 563       common /pots/ Cb0, Yuk0, Pau0, Sky20, Sky30, gamSky,
 564      +              gamYuk, drPau, dpPau, gw, sgw, delr, fdel,
 565      +              dt,da, db,dtimestep
 566
 567
 568 c spectator arrays
 569         integer smax
 570         parameter(smax=500)  ! maximum number of spectators
 571         real*8 r0s(smax), rxs(smax), rys(smax), rzs(smax),
 572      +         p0s(smax), pxs(smax), pys(smax), pzs(smax),
 573      +         sfmass(smax)
 574
 575
 576         integer sspin(smax), scharge(smax), sityp(smax), siso3(smax),
 577      +          suid(smax)
 578
 579         integer nspec
 580
 581         common /scoor/ r0s, rxs, rys, rzs, p0s, pxs ,pys, pzs, sfmass
 582
 583         common /sisys/ sspin, scharge, sityp, siso3, suid
 584
 585         common /ssys/ nspec
 586
 587
 588         real*8 p0td(2,nmax),pxtd(2,nmax),pytd(2,nmax),pztd(2,nmax),
 589      +         fmasstd(2,nmax)
 590         integer ityptd(2,nmax),iso3td(2,nmax)
 591         integer itypt(2),uidt(2),origint(2),iso3t(2)
 592
 593         common /rtdelay/p0td,pxtd,pytd,pztd,fmasstd
 594         common /itdelay/ityptd,iso3td
 595         common /svinfo/itypt,uidt,origint,iso3t
 596         real*8 ffermpx(nmax), ffermpy(nmax), ffermpz(nmax)
 597         real*8 peq1, peq2
 598         common /ffermi/ ffermpx, ffermpy, ffermpz
 599         common /peq/ peq1,peq2
 600
 601
 602         integer numcto,numctp,maxstables
 603         parameter(numcto=400) ! maximum number of options
 604         parameter(numctp=400) ! maximum number of parameters
 605         parameter(maxstables=20) ! maximum number of stable particles
 606         integer   CTOption(numcto)
 607         real*8    CTParam(numctp)
 608         common /options/CTOption,CTParam
 609
 610         real*8 frr0(nmax), frrx(nmax), frry(nmax), frrz(nmax),
 611      +     frp0(nmax), frpx(nmax), frpy(nmax), frpz(nmax)
 612
 613       common /frcoor/ frr0, frrx, frry, frrz, frp0, frpx, frpy, frpz
 614               integer maxbar,maxbra,minbar
 615       integer offmeson,maxmeson,pimeson,maxbrm,minnuc,mindel
 616       integer maxbrs1,maxbrs2
 617       integer numnuc,numdel,nucleon,maxnuc,maxdel
 618       integer minmes,maxmes
 619
 620
 621       parameter (minnuc=1) ! lowest baryon particle ID
 622       parameter (minmes=100) ! lowest meson particle ID
 623       parameter (maxmes=132) ! hightest meson particle ID
 624
 625 c number of resonances of a kind
 626       parameter (numnuc=16) ! number of nucleon resonances
 627       parameter (numdel=10) ! number of delta resonances
 628 c indices of minimal and maximal itype of a kind (redundant but nice)
 629       parameter (maxnuc=minnuc+numnuc-1) ! highest nucleon ID
 630       parameter (mindel=minnuc+maxnuc)   ! lowest delta ID
 631       parameter (maxdel=mindel+numdel-1) ! highest delta ID
 632
 633 c minres & maxres define the range of nonstable & nonstrange baryons
 634       integer minres,maxres
 635       parameter (minres=minnuc+1) ! lowest baryon resonance ID
 636       parameter (maxres=maxdel)   ! highest (nonstrange) baryon
 637                                   ! resonance ID
 638
 639 c strangenes.ne.0 baryon resonances
 640       integer minlam,minsig,mincas,minome
 641       integer numlam,numsig,numcas,numome
 642       integer maxlam,maxsig,maxcas,maxome
 643       parameter (numlam=13) ! number of lambda states
 644       parameter (numsig=9)  ! number of sigma states
 645       parameter (numcas=6)  ! number of cascade states
 646       parameter (numome=1)  ! number of omega states
 647       parameter (minlam=mindel+numdel)   ! ID of lowest lambda state
 648       parameter (maxlam=minlam+numlam-1) ! ID of highest lambda state
 649       parameter (minsig=minlam+numlam)   ! ID of lowest sigma state
 650       parameter (maxsig=minsig+numsig-1) ! ID of highest sigma state
 651       parameter (mincas=minsig+numsig)   ! ID of lowest cascade state
 652       parameter (maxcas=mincas+numcas-1) ! ID of highest cascade state
 653       parameter (minome=mincas+numcas)   ! ID of lowest omega state
 654       parameter (maxome=minome+numome-1) ! ID of highest omega state
 655
 656 c minbar & maxbar define the range of all baryons
 657       parameter (minbar=minnuc) ! ID of lowest baryon state
 658       parameter (maxbar=maxome) ! ID of highest baryon state
 659
 660       parameter (offmeson=minmes) ! offset between zero and lowest
 661                                   ! meson state
 662       parameter (maxmeson=maxmes) ! ID of highest meson state
 663 c... these variables are in principal obsolete and should be exchanged
 664 c were referenced
 665
 666 c... avoid hard coded itypes
 667       integer itrho,itome,iteta,itkaon,itphi,itetapr
 668       parameter (itkaon=106)   ! ID of kaon
 669       parameter (itrho=104)    ! ID of rho meson
 670       parameter (itome=103)    ! ID of omega meson
 671       parameter (iteta=102)    ! ID of eta
 672       parameter (itphi=109)    ! ID of phi
 673       parameter (itetapr=107)  ! ID of eta'
 674       parameter (pimeson=101)  ! ID of $\pi$
 675       parameter (nucleon=minnuc) ! ID of nucleon
 676
 677       integer itmin,itmax
 678       parameter (itmin=minnuc)  ! lowest defined ID
 679       parameter (itmax=maxmes)  ! highest defined ID
 680 c
 681       parameter (maxbra=11)  ! decay channels for $s=0$ baryon resonances
 682       parameter (maxbrm=25) ! decay channels for meson resonances
 683       parameter (maxbrs1=10)! decay channels for $s=1$ baryon resonances
 684       parameter (maxbrs2=3) ! decay channels for $s=2$ baryon resonances
 685
 686 c
 687        integer mlt2it(maxmes-minmes) ! meson IDs sorted by multipletts
 688
 689
 690       real*8 massoff,mresmin,mresmax
 691       parameter (massoff=1d-4)      ! offset for mass generation
 692       parameter (mresmin=1.0765d0)  ! minimum baryon resonance mass
 693       parameter (mresmax=5d0)       ! maximum baryon resonance mass
 694
 695       character*45 versiontag
 696       common /versioning/ versiontag
 697
 698       real*8 massres(minbar:maxbar),widres(minbar:maxbar)
 699       real*8 branmes(0:maxbrm,minmes+1:maxmes)
 700       real*8 branres(0:maxbra,minnuc+1:maxdel)
 701       real*8 branbs1(0:maxbrs1,minlam+1:maxsig)
 702       real*8 branbs2(0:maxbrs2,mincas+1:maxcas)
 703       integer Jres(minbar:maxbar)
 704       integer Jmes(minmes:maxmes)
 705       integer pares(minbar:maxbar),pames(minmes:maxmes)
 706       integer Isores(minbar:maxbar), Isomes(minmes:maxmes)
 707       integer brtype(4,0:maxbra),bmtype(4,0:maxbrm)
 708       integer bs1type(4,0:maxbrs1),bs2type(4,0:maxbrs2)
 709       real*8 massmes(minmes:maxmes)
 710       real*8 mmesmn(minmes:maxmes)
 711       real*8 widmes(minmes:maxmes)
 712       integer strres(minbar:maxbar),strmes(minmes:maxmes)
 713
 714       integer lbr(0:maxbra,minnuc+1:maxdel)
 715       integer lbs1(0:maxbrs1,minlam+1:maxsig)
 716       integer lbs2(0:maxbrs2,mincas+1:maxcas)
 717       integer lbm(0:maxbrm,minmes+1:maxmes)
 718
 719       common /resonances/ massres,widres,massmes,widmes,mmesmn,
 720      ,                    branres,branmes,branbs1,branbs2,
 721      ,                    bs1type,bs2type,lbs1,lbs2,lbm,
 722      ,                    jres,jmes,lbr,brtype,pares,pames,
 723      ,                    bmtype,
 724      ,                    Isores,Isomes,strres,strmes,mlt2it
 725
 726
 727       !epos common blocks for particle list
 728       integer mmry,mxptl
 729       parameter (mmry=1)   !memory saving factor
 730       parameter (mxptl=200000/mmry) !max nr of particles in epos ptl list
 731
 732       integer iorptl(mxptl),idptl(mxptl),istptl(mxptl),
 733      *  ifrptl(2,mxptl),jorptl(mxptl),ibptl(4,mxptl),ityptl(mxptl)
 734       real pptl(5,mxptl),tivptl(2,mxptl),xorptl(4,mxptl)
 735
 736       common/cptl/nptl,pptl,iorptl,idptl
 737      *,istptl,tivptl,ifrptl,jorptl
 738      *,xorptl,ibptl,ityptl
 739
 740
 741       integer nptl
 742
 743
 744       integer nevt,kolevt,koievt,npjevt,ntgevt,npnevt,nppevt,
 745      *        ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt,
 746      *        nglevt,minfra,maxfra
 747       real phievt,bimevt,pmxevt,egyevt,xbjevt,qsqevt,zppevt,zptevt
 748
 749       common/cevt/phievt,nevt,bimevt,kolevt,koievt,pmxevt,egyevt,npjevt
 750      *,ntgevt,npnevt,nppevt,ntnevt,ntpevt,jpnevt,jppevt,jtnevt,jtpevt
 751      *,xbjevt,qsqevt,nglevt,zppevt,zptevt,minfra,maxfra
 752
 753       integer istore,istmax,irescl,ntrymx,nclean,iopdg
 754       real gaumx
 755       common/othe1/istore,istmax,gaumx,irescl,ntrymx,nclean,iopdg
 756
 757       integer nn,idpdgg,idepos
 758
 759
 760
 761       integer idtrafo
 762       external idtrafo
 763       integer fchg
 764       external fchg
 765       integer pdgid
 766       external pdgid
 767       real*8 dectim
 768       external dectim
 769
 770       nbar=0
 771       do nn=1,npart
 772 csp        write(6,*) 'ityp', ityp(nn), 'iso3', iso3(nn)
 773            idpdgg=pdgid(ityp(nn),iso3(nn))
 774            idepos=idtrafo('pdg','nxs',idpdgg)
 775
 776                nbar=nbar+1
 777
 778                xorptl(4,nbar)= r0(nn)
 779                xorptl(1,nbar)=rx(nn)
 780                xorptl(2,nbar)=ry(nn)
 781                xorptl(3,nbar)=rz(nn)
 782                pptl(4,nbar)= p0(nn)
 783                pptl(1,nbar)=px(nn)
 784                pptl(2,nbar)=py(nn)
 785                pptl(3,nbar)=pz(nn)
 786                pptl(5,nbar)=fmass(nn)
 787                idptl(nbar)=idepos
 788
 789
 790         nptl=nbar
 791       enddo
 792
 793       end
 794
 795 c----------------------------------------------------------------------------------
 796       subroutine  urqmdepos
 797 c----------------------------------------------------------------------------------
 798 c   main urqmd modul, essentially copied from urqmdepos (urqmd.f)
 799 c----------------------------------------------------------------------------------
 800
 801       implicit none
 802       include '../urqmd23/coms.f'
 803       include '../urqmd23/comres.f'
 804       include '../urqmd23/options.f'
 805       include '../urqmd23/colltab.f'
 806       include '../urqmd23/inputs.f'
 807       include '../urqmd23/newpart.f'
 808       include '../urqmd23/boxinc.f'
 809       integer i,j,k,steps,ii,ocharge,ncharge, mc, mp, noc, it1,it2
 810       real*8 sqrts,otime,xdummy,st
 811       logical isstable
 812       integer stidx
 813       real*8 Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau
 814       common /energies/ Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau
 815       integer cti1sav,cti2sav
 816
 817       mc=0
 818       mp=0
 819       noc=0
 820
 821       time = 0.0  !time is the system time at the BEGINNING of every timestep
 822
 823       !initialize random number generator
 824       !call auto-seed generator only for first event and if no seed was fixed
 825       if(.not.firstseed.and.(.not.fixedseed)) then
 826          ranseed=-(1*abs(ranseed))
 827          call sseed(ranseed)
 828       else
 829          firstseed=.false.
 830       endif
 831
 832       !old time if an old fort.14 is used
 833       if(CTOption(40).eq.1)time=acttime
 834       if(CTOption(40).eq.3)time=acttime
 835
 836       !write headers to file
 837       call output(13)
 838       !call output(14)
 839       call output(15)
 840       call output(16)
 841       !if(event.eq.1)call output(17)
 842       call osc99_event(-1)
 843
 844       !for CTOption(4)=1 : output of initialization configuration
 845       if(CTOption(4).eq.1)call file14out(0)
 846       !participant/spectator model:
 847       if(CTOption(28).ne.0) call rmspec(0.5d0*bimp,-(0.5d0*bimp))
 848
 849       otime = outsteps*dtimestep  !compute time of output
 850
 851       steps = 0  !reset time step counter
 852
 853   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 854   ! loop over all timesteps
 855   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 856
 857       do 20  steps=1,nsteps
 858
 859          if (eos.ne.0) then
 860            do j=1,npart
 861                r0_t(j) = r0(j)
 862                rx_t(j) = rx(j)
 863                ry_t(j) = ry(j)
 864                rz_t(j) = rz(j)
 865            enddo
 866          end if
 867
 868          !we are at the beginning of the timestep, set current time (acttime)
 869          acttime = time
 870
 871          if(CTOption(16).ne.0) goto 103  !option for MD without collision term
 872
 873          call colload  ! Load collision table with next collisions in current timestep
 874
 875          ! check for collisions in time-step, nct = # of collisions in table
 876          if (nct.gt.0) then
 877
 878  101        continue              !entry-point for collision loop in case
 879             k = 0                 !      of full colload after every coll
 880  100        continue              !normal entry-point for collision loop
 881             call getnext(k)       !get next collision
 882             if (k.eq.0) goto 102  !exit collision loop if no collisions are left
 883
 884             !propagate all particles to next collision time
 885             !store actual time in acttime, propagation time st=cttime(k)-acttime
 886             st=cttime(k)-acttime
 887             call cascstep(acttime,st)
 888             acttime = cttime(k)   !new actual time (for upcoming collision)
 889
 890             !perform collision
 891
 892             if(cti2(k).gt.0.and.
 893      .           abs(sqrts(cti1(k),cti2(k))-ctsqrts(k)).gt.1d-3)then
 894                write(6,*)' ***(E) wrong collision update (col) ***'
 895                write(6,*)cti1(k),cti2(k),
 896      .              ctsqrts(k),sqrts(cti1(k),cti2(k))
 897             else if(cti2(k).eq.0.and.
 898      .              abs(fmass(cti1(k))-ctsqrts(k)).gt.1d-3) then
 899                write(6,*)' *** main(W) wrong collision update (decay)'
 900                write(6,*)ctag,cti1(k),ityp(cti1(k)),dectime(cti1(k)),
 901      .              fmass(cti1(k)),ctsqrts(k)
 902             endif
 903
 904             ocharge=charge(cti1(k))
 905             if(cti2(k).gt.0) ocharge=ocharge+charge(cti2(k))
 906
 907             !store quantities in local variables for charge conservation check
 908             it1= ityp(cti1(k))
 909             if(cti2(k).gt.0)it2= ityp(cti2(k))
 910
 911             !increment "dirty" collision counter
 912             if(cti2(k).gt.0)then !scatter
 913                mc=mc+1
 914             endif
 915             !perform scattering/decay
 916             cti1sav = cti1(k)
 917             cti2sav = cti2(k)
 918             call scatter(cti1(k),cti2(k),ctsigtot(k),ctsqrts(k),
 919      &                   ctcolfluc(k))
 920
 921             !update collision table
 922
 923             !normal update mode
 924             if(CTOption(17).eq.0) then
 925                if(nexit.eq.0) then
 926                  !new collision partners for pauli-blocked states (nexit=0)
 927                  if (cti1(k).ne.cti1sav.or.cti2(k).ne.cti2sav) then
 928                    cti1(k) = cti1sav
 929                    cti2(k) = cti2sav
 930                  endif
 931
 932                  call collupd(cti1(k),1)
 933                  if(cti2(k).gt.0) call collupd(cti2(k),1)
 934                else
 935                  ncharge=0
 936                  !new collision partners for scattered/produced particles (nexit><0)
 937                  do i=1,nexit
 938                    !ncharge is used for charge conservation check
 939                    ncharge=ncharge+charge(inew(i))
 940                    call collupd(inew(i),1)
 941                  enddo
 942                endif
 943                !update collisions for partners of annihilated particles
 944                do ii=1,nsav
 945                   call collupd(ctsav(ii),1)
 946                enddo
 947                nsav=0
 948             else ! (CTOption(17).ne.0)
 949               !full collision load
 950               call colload
 951             endif
 952
 953             if (CTOption(17).eq.0) goto 100
 954             goto 101
 955
 956             !this is the point to jump to after all collisions in the timestep
 957             !have been taken care of
 958  102        continue
 959
 960          endif ! (nct.gt.0)
 961
 962
 963          !After all collisions in the timestep are done, propagate to end of
 964          !the timestep.
 965
 966          !point to jump to in case of MD without collision term
 967  103     continue
 968
 969          time = time+dtimestep  !increment timestep
 970
 971          !After all collisions in the timestep are done, propagate to end of
 972          !the timestep.
 973          call cascstep(acttime,time-acttime)
 974
 975          !in case of potential interaction, do MD propagation step
 976          if (eos.ne.0) then
 977             ! set initial conditions for MD propagation-step
 978             do j=1,npart
 979                r0(j) = r0_t(j)
 980                rx(j) = rx_t(j)
 981                ry(j) = ry_t(j)
 982                rz(j) = rz_t(j)
 983             enddo
 984             !now molecular dynamics trajectories
 985             call proprk(time,dtimestep)
 986          endif ! (eos.ne.0)
 987
 988          !perform output if desired
 989          if(mod(steps,outsteps).eq.0.and.steps.lt.nsteps)then
 990             if(CTOption(28).eq.2)call spectrans(otime)
 991             call file14outx(steps)
 992          endif ! output handling
 993
 994  20   continue ! time step loop
 995
 996       acttime=time
 997
 998   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 999   ! optional decay of all unstable
1000   !  particles before final output
1001   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1002
1003       !DANGER: pauli-blocked decays are not performed !!!
1004       if(CTOption(18).eq.0.and.CTOption(51).eq.0) then
1005          !no do-loop is used because npart changes in loop-structure
1006          i=0
1007          nct=0
1008          actcol=0
1009          CTOption(10)=1  !disable Pauli-Blocker for final decays
1010  40      continue  !decay loop structure starts here
1011          i=i+1
1012          if(dectime(i).lt.1.d30) then !if particle unstable
1013  41         continue
1014             isstable = .false.
1015             do stidx=1,nstable
1016                if (ityp(i).eq.stabvec(stidx)) then
1017                   !write (6,*) 'no decay of particle ',ityp(i)
1018                   isstable = .true.
1019                endif
1020             enddo
1021             if (.not.isstable) then
1022                call scatter(i,0,0.d0,fmass(i),xdummy) !perform decay
1023                !backtracing if decay-product is unstable itself
1024                if(dectime(i).lt.1.d30) goto 41
1025             endif
1026          endif
1027          !check next particle
1028          if(i.lt.npart) goto 40
1029       endif ! final decay
1030
1031   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1032   !     final output
1033   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1034
1035       if(CTOption(28).eq.2)call spectrans(otime)
1036
1037       call file13out(nsteps)
1038       !call file14out(nsteps)
1039       call file16out
1040       call osc_event
1041       call osc99_event(1)
1042       call osc99_eoe
1043
1044       mp=mp+npart
1045       if(ctag.eq.0)then
1046          write(*,*)'(W) No collision in event ',event
1047          noc=noc+1
1048       endif
1049
1050       end
1051
1052 c-------------------------------------------------------------------------------------
1053       subroutine file14outx(timestep)
1054 c-------------------------------------------------------------------------------------
1055       implicit none
1056       include '../urqmd23/comres.f'
1057       include '../urqmd23/coms.f'
1058       include '../urqmd23/options.f'
1059       include '../urqmd23/inputs.f'
1060       include '../urqmd23/newpart.f'
1061       include '../urqmd23/freezeout.f'
1062       include '../urqmd23/boxinc.f'
1063       integer i,ttime,timestep,itotcoll,iinelcoll
1064       real*8 sigmatot,t
1065       common /outco2/sigmatot
1066       include '../urqmd23/outcom.f'
1067       save
1068
1069       if(bf14)return
1070       ttime=int(timestep*dtimestep+0.01)
1071       itotcoll=ctag-dectag
1072       iinelcoll=itotcoll-NBlColl-NElColl
1073       print*,'(file14outx)',ttime,npart
1074      @ ,itotcoll,NElColl,iinelcoll,NBlColl,dectag,
1075      @     NHardRes,NSoftRes,NDecRes
1076        !----------------------------------------------------------------------------
1077        !  r0(i), rx(i), ry(i), rz(i)   ........................................ x4
1078        !  p0(i),px(i)+ffermpx(i),py(i)+ffermpy(i),pz(i)+ffermpz(i),fmass(i) ... p5
1079        !  ityp(i)  ..... particle id (listed in the urqmd user guide)
1080        !  iso3(i) ...... 2 times the isospin of a particle
1081        !  charge(i) .... charge of the particle
1082        !  lstcoll(i) ... index of last collision partner
1083        !  ncoll(i) ..... number of collisions
1084        !  origin(i) ....
1085        !  dectime(i) ...
1086        !  tform(i) ..... formation time
1087        !  xtotfac(i) ... cross section (zero if the particle is not yet formed)
1088        !  uid(i) ......
1089        !---------------------------------------------------------------------------
1090        do i=1,npart
1091          t=r0(i)
1092        enddo
1093
1094       end
1095