EPOS/epos167/epos-bas-168.f

   1
   2 c-----------------------------------------------------------------------
   3       subroutine aaset(iop)
   4 c-----------------------------------------------------------------------
   5 c     sets parameters and options, initializes counters ...
   6 c-----------------------------------------------------------------------
   7
   8       include 'epos.inc'
   9       include 'epos.incpar'
  10       include 'epos.incsem'
  11       common/record/maxrec(2),irecty(30,2)
  12       common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi/cifset/ifset
  13       common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
  14       character*80 fndat,fnncs,fnIIdat,fnIIncs                 !qgs-II????????
  15       common/qgsfname/  fndat, fnncs, ifdat, ifncs
  16       common/qgsIIfname/fnIIdat, fnIIncs, ifIIdat, ifIIncs     !qgs-II????????
  17       common/ghecsquel/anquasiel,iquasiel
  18       common/cbincond/nozero,ibmin,ibmax
  19       !gauss weights
  20       data (tgss(2,i),i=1,7)/ .3399810436,.8611363116    ,5*0.     /
  21       data (wgss(2,i),i=1,7)/ .6521451549,.3478548451    ,5*0.     /
  22       data (tgss(3,i),i=1,7)/ .2386192,.6612094,.9324700  ,4*0.    /
  23       data (wgss(3,i),i=1,7)/ .4679139,.3607616,.1713245  ,4*0.    /
  24       data (tgss(5,i),i=1,7)/ .1488743389,.4333953941,.6794095682
  25      *                       ,.8650633666,.9739065285    ,2*0.     /
  26       data (wgss(5,i),i=1,7)/ .2955242247,.2692667193,.2190863625
  27      *                       ,.1494513491,.0666713443    ,2*0.     /
  28       data (tgss(7,i),i=1,7)/ .9862838,.9284349,.8272013,.6872929
  29      *                       ,.5152486,.3191124,.1080549           /
  30       data (wgss(7,i),i=1,7)/ .03511946,.08015809,.1215186,.1572032
  31      *                       ,.1855384,.2051985,.2152639           /
  32
  33       if(iop.eq.1)write(ifmt,'(a)')'default settings ...'
  34       if(iop.eq.1)goto1001
  35       if(iop.eq.2)goto1002
  36
  37 c  version
  38
  39       iversn=int(1.67*100) !version number
  40       iverso=int(1.67*100) !last official version
  41
  42 c  application
  43
  44       iappl=1          !choice for application (0,1,2,3,4)
  45
  46 c  model
  47
  48       model=1
  49       iquasiel=1       !allow (1) or not (0) quasi-elastic event in model 3
  50
  51 c  file names and units
  52
  53       fnnx='path/epos '      !path epos name
  54       fnch='zzz.check '         !check-file name
  55       fnhi='zzz.histo '         !histo-file name
  56       fndt='zzz.data '          !data-file name
  57       fncp='zzz.copy '          !copy-file name
  58       fnii='zzz.initl '       !initl-file name
  59       fnid='zzz.inidi '       !inidi-file name
  60       fndr='zzz.inidr '       !inidr-file name
  61       fnie='zzz.iniev '       !iniev-file name
  62       fnrj='zzz.inirj '       !inirj-file name
  63       fncs='zzz.inics '       !inics-file name
  64       nfnnx=index(fnnx,' ')-1   !length of path epos name
  65       nfnch=index(fnch,' ')-1   !length of check-file name
  66       nfnhi=index(fnhi,' ')-1   !length of histo-file name
  67       nfndt=index(fndt,' ')-1   !length of data-file name
  68       nfncp=index(fncp,' ')-1   !length of copy-file name
  69       nfnii=index(fnii,' ')-1   !length of initl-file name
  70       nfnid=index(fnid,' ')-1   !length of inidi-file name
  71       nfndr=index(fndr,' ')-1   !length of inidr-file name
  72       nfnie=index(fnie,' ')-1   !length of iniev-file name
  73       nfnrj=index(fnrj,' ')-1   !length of inirj-file name
  74       nfncs=index(fncs,' ')-1   !length of inics-file name
  75
  76       if(iop.eq.0)then
  77       ifop=5     !optns-file unit
  78       ifmt=6     !std output unit
  79       ifcx=31    !check-file unit (for open)
  80       ifch=31    !check-file unit (for write)
  81       ifhi=35    !histo-file unit
  82       ifdt=51    !data-file unit
  83       ifcp=52    !copy-file unit
  84       endif
  85 c  initial seed
  86
  87       seedi=0d0   !.ne.0.: seed for random number generator: at start program
  88       seedj=0d0   !.ne.0.: seed for random number generator: for first event
  89
  90
  91  1001 continue ! ----> entry iop=1 <-----
  92
  93       call factoriel
  94
  95
  96 c  some basic things
  97
  98       nevent=1    !number of events
  99       nfull=-1          !number of full events (ico+hydro+f.o.)(if different from -1)
 100       nfreeze=1    !number of freeze out events per epos event
 101       ninicon=1    ! +-number of events per initial condition
 102                     ! if positive: keep same b, otherwise generate b each time
 103       engy=-1      !energy
 104       elab=-1      !energy lab
 105       ecms=-1      !energy centre of mass
 106       ekin=-1      !energy kinetic
 107       pnll=-1      !momentum
 108       egymin=6.    !minimum energy
 109       egymax=2.E+06 !maximum energy
 110       noebin=1     !number of energy bins
 111       engmin=0     !minimum energy
 112       engmax=0     !maximum energy
 113       iologe=-1     !0=linear bins   1=log bins (cms engy)
 114       iologl=-1     !0=linear bins   1=log bins (Kinetic engy)
 115
 116 c  printout options
 117
 118       if(iop.eq.0)iprmpt=-2   !prompt (>0) or not (<0)
 119       ish=0      !1,2,3,4 ...: more and more detailed messages
 120       irandm=0   !prints all random numbers (1)
 121       irewch=0   !rewinds check file before each event (1)
 122       if(iop.eq.0)iecho=1    !verify option for input reading
 123       modsho=1   !message all modsho events
 124       idensi=0   !must be 1 when subr xjden1 is used
 125       ishevt=0   !minimum event number for certain print options
 126       iwseed=1   !print out seed (1) or not
 127
 128 c  fragmentation and decay parameters
 129
 130       pud=0.434  !prob u d
 131       pude=0.434  !prob u d for ee
 132       pudk=0.434  !prob u d for kinky string
 133       pudr=0.434 !prob u d for remnant string
 134       pudd=0.44   !prob u d for break diffractive
 135       puds=0.14   !s suppression in diquark break
 136       strcut=0.00012   !cut factor for diffractive string fragmentation
 137       diqcut=0.12       !cut factor for diffractive string fragmentation
 138       pdiqua= 0.21   !qq-qqbar probability
 139       pdiquae=0.09  !qq-qqbar probability for ee
 140       pdiquak=0.09   !qq-qqbar probability for kinky string
 141       pdiquar=0.09   !qq-qqbar probability for remnant string
 142       pdiquad=0.09   !qq-qqbar probability for diffractive remnant string
 143       ptfra= 0.30  !string break pt
 144       ptfrae=0.30  !string break pt for ee
 145       ptfrak=0.30  !string break pt for kinky string in hadronic ia
 146       ptfrakd=0.30  !string break pt for kinky string diquark break
 147       ptfrair=0.30 !string break pt for inel remnant string
 148       ptfradr=0.40 !string break pt for diffr remnant string
 149       ptfrasr=0.25 !string break del pt for strange breakup of remn string
 150       ptfraqq=0.40 !string break pt for diquark break
 151       pbreak=-1.  !break-parameter
 152       pbreake=-1.  !break-parameter ee
 153       pbreakk=-1.  !break-parameter for kinky string
 154       pbreakr=-1. !break-parameter for remnant string
 155       pbreakd=0.27 !break-parameter for diff remnant string
 156
 157 c  fragmentation and decay options
 158
 159       ndecay=0   !ndecay suppresses certain resonance decays
 160                  !0000001: all resonance decays suppressed
 161                  !0000010: k_short/long (+-20) decays suppressed
 162                  !0000100: lambda (+-2130) decays suppressed
 163                  !0001000: sigma (+-1130,+-2230) decays suppressed
 164                  !0010000: cascade (+-2330,+-1330) decays suppressed
 165                  !0100000: omega (+-3331) decays suppressed
 166                  !1000000: pi0 (110) decays suppressed
 167                  !also several "1"s may be used with obvious meaning
 168       maxres=99999 !maximum resonance spin
 169       aouni=0.   !technical parameter for development
 170
 171 c  lepton-nucleon and e+e- options
 172
 173       iolept=1     !q2-x ditribution calculated (1) or taken from table (<0)
 174       ydmin=0      ! range of y
 175       ydmax=0
 176       qdmin=0      ! range of q**2
 177       qdmax=0
 178       themin=0     !minimum scattering angle
 179       themax=0     !maximum scattering angle
 180       elomin=0     !minimum energy of outgoing lepton
 181       elepti=0     !incoming lepton energy
 182       elepto=0     !outgoing lepton energy
 183       angmue=3.9645/360.*2*3.14159 !mue angle
 184       icinpu=0
 185       itflav=0     ! initial flavor for e+e-
 186       idisco=0     !deep inelastic contributions
 187                    !0=all, 1=direct-light, 2=direct-charm, 3=resolved
 188
 189 c  hadron-hadron options +++
 190
 191       isetcs=2    !option to obtain pomeron parameters
 192                   ! 0.....determine parameters but do not use Kfit
 193                   ! 1.....determine parameters and use Kfit
 194                   ! else..get from table
 195                   !         should be sufficiently detailed
 196                   !          say iclegy1=1,iclegy2=99
 197                   !         table is always done, more or less detailed!!!
 198                   !and option to use cross ection tables
 199                   ! 2....tabulation
 200                   ! else...not
 201       iclegy1=1   !energy class lower limit ( 1 minimal    1 realistic    )
 202       iclegy2=99   !energy class upper limit ( 1  option   99 use of table )
 203       isigma=1    !option for xsection printing (always calculated now)
 204                   !  0=no, 1=yes
 205
 206 c  hadron-hadron options
 207
 208       idprojin=1120 !projectile particle id
 209       idtargin=1120 !target particle id
 210       idproj=1120 !projectile particle id
 211       idtarg=1120 !target particle id
 212       iregge=0    !consider reggeons (1=yes 0=no)
 213       isopom=1    !consider soft pomerons (1=yes 0=no)
 214       ishpom=1    !consider semihard pomerons (1=yes 0=no)
 215       iscreen=1   !consider screening corrections (1=yes 0=no)
 216       isplit=1    !consider splitting corrections (1=yes 0=no)
 217       irzptn=0    !recalculate Zptn (1=yes 0=no)  ????????maybe obsolete??????????
 218       irmdrop=1    !consider droplet for remnants (1=yes 0=no)
 219       nprmax=10000 !maximum number of pomerons/reggeons
 220       iemspl=0
 221       intpol=3     !number of points for interpolation in psfli, psfaz
 222       ioems=2
 223       iomega=1    !option for omega calculation (if 2, no diffraction in G)
 224         !hadron excitation classes (used in psvin)
 225       icdp=2     !projectile hadron excitation
 226       icdt=2     !target hadron excitation
 227               !hadron classes: 1=pions,2=nucleons,3=kaons
 228       iclpro1=1   !projectile hadron class lower limit
 229       iclpro2=3   !projectile hadron class upper limit
 230       icltar1=2   !target hadron class lower limit (should not be change (see epos-sem))
 231       icltar2=2   !target hadron class upper limit (should not be change (see epos-sem))
 232       egylow=1.5  !lower limit of lowest energy bin
 233       delh=0.25   !effective overcriticity for the hard pom (techn param)
 234       factgam=1         !enhancement factor for gammas
 235
 236 c  hadron-hadron parameters +++
 237
 238       alpfomi=   8.    !normalization of function fom for Phi^n (z0=alpfomi)
 239       betfom=    9.d0   !slope of function fom for Phi^n
 240       gamfom=    5.d0   !Z slope of function fom for Phi^n
 241       betpom=    0.25    !gluon structure function parameter for the so
 242       glusea=    0.1   !sea quarks / (sea quarks + gluons) at x -> 0
 243       r2had(2)=  1.   !r2 proton
 244       r2hads(2)= 1.06     !diff corr factor proton
 245       slopom=    0.06    !alpha prime pomeron
 246       slopoms=   0.175    !alpha prime pomeron dif
 247       gamhad(2)= 1.1    !gamma proton           increase->  sig up, n up, softPom up
 248       gamhadsi(2)=-1.   !gamma soft proton (<0 = same as gamhad)
 249       gamtil=    0.08   !                       increase -> sig up, n up, hard Pom up
 250       alppar=    0.5  !alpha particip (not 1 !) increase -> sig up, n down, width y down
 251       alppom=    1.06 !alpha pomeron
 252       ptsend=    0.8    !string end pt
 253       ptsendi=   0.05   !string end pt for non excited state
 254       ptsemx=    2.     !cut factor for gaussian distribution
 255       ptdiff=    0.27   !pt for diffraction
 256       q2min=     4.0    !q**2 cutoff
 257       q2ini=     0.25   !resolution scale
 258       q2fin=     0.02   !q**2 cutoff timelike      decrease ->  high pt down
 259       amdrmax=   2.5    !maximum mass leading droplet
 260       facdif=    0.35   !factor for diffractive profile
 261       reminv=    0.33     !remnant inversion probability for strangeness
 262       zrminc=    0.033  !remnant inversion probability string mass evolution
 263       zrmmax=    10.    !remnant inversion probability system mass evolution
 264       edmaxi=    1.e12   !defines edmax in epos-sem
 265       epmaxi=    1.e12   !defines epmax in epos-sem
 266
 267 c  hadron-hadron parameters
 268
 269       qcdlam=.04           !lambda_qcd squared
 270       naflav=4          !number of active flavors
 271       factk=2.         !k-factor value
 272       alfe=1./137.
 273       pt2cut=0.         !p_t-cutoff for the born process
 274       rstrau(1)=1.   !pion !effective ratio of u sea over u sea basis
 275       rstrad(1)=1.         !effective ratio of d sea over u sea basis
 276       rstras(1)=0.8         !effective ratio of s sea over u sea basis
 277       rstrau(2)=1.   !nucl !effective ratio of u sea over u sea basis
 278       rstrad(2)=1.         !effective ratio of d sea over u sea basis
 279       rstras(2)=0.8         !effective ratio of s sea over u sea basis
 280       rstrau(3)=1.   !kaon !effective ratio of u sea over u sea basis
 281       rstrad(3)=1.         !effective ratio of d sea over u sea basis
 282       rstras(3)=0.8         !effective ratio of s sea over u sea basis
 283       rstrau(4)=1.   !j/psi!effective ratio of u sea over u sea basis
 284       rstrad(4)=1.         !effective ratio of d sea over u sea basis
 285       rstras(4)=0.8         !effective ratio of s sea over u sea basis
 286       rstrasi=0.0      !effective ratio of strange sea over u sea increase
 287       wgtval=0.          !weight for valq - antiq as soft string ends
 288       wgtsea=1.          !weight for seaq - antiq as soft string ends
 289       wgtdiq=0.          !weight for diq - antidiq as soft string ends
 290       wgtqqq(1)=0.17          !weight for seaq - valqq as soft string end for pion
 291       wgtqqq(2)=0.17          !weight for seaq - valqq as soft string end for nucleon
 292       wgtqqq(3)=0.17          !weight for seaq - valqq as soft string end for kaon
 293       wgtqqq(4)=0.17          !weight for seaq - valqq as soft string end for J/Psi
 294       exmass=0.02         !excitation mass for remnant
 295       r3pom=0.01      !triple pomeron coupling
 296         r3pomi=r3pom  !store
 297       r4pom=0.001     !4-pomeron coupling
 298         r4pomi=r4pom  !store
 299       gfactor=0.4     !exp(-gfactor*Z) in front of eps in om5jk (to choose soft or hard)
 300       gwidth=0.85   !diff relative b-width
 301 !     gamhad(1) defined in psaini
 302 !     gamhad(3) defined in psaini
 303       gamhadsi(1)=-1.    !gamma soft pion
 304       gamhadsi(3)=-1.    !gamma soft kaon
 305       gamhadsi(4)=-1.    !gamma soft charm
 306       r2had(1)=1.    !r2 pion
 307       r2had(3)=1.30   !r2 kaon
 308       r2had(4)=0.     !r2 charm
 309       r2hads(1)=0.85  !diff corr factor pion
 310       r2hads(3)=1.15  !diff corr factor kaon
 311       r2hads(4)=1.25  !diff corr factor kaon
 312       chad(1)=1.      !c pion
 313       chad(2)=1.      !c proton
 314       chad(3)=1.      !c kaon
 315       chad(4)=1.      !c charm
 316       wdiff(1)=0.36   !diffractive probability
 317       wdiff(2)=0.645  !diffractive probability
 318       wdiff(3)=0.145  !diffractive probability
 319       wdiff(4)=0.1    !diffractive probability
 320       alplea(1)=0.75  !alpha leading pion
 321       alplea(2)=1.25  !alpha leading proton
 322       alplea(3)=1.    !alpha leading kaon
 323       alplea(4)=1.    !alpha leading jpsi
 324       rexndf=1.       !relative value of rexndi compare to rexdif if >0
 325 c following parameters recalculated in xsigma ...
 326       rexdifi(1)=-0.88  !remnant excitation probability diffractive pion
 327       rexdifi(2)=-0.81  !remnant excitation probability diffractive proton
 328       rexdifi(3)=-0.81  !remnant excitation probability diffractive kaon
 329       rexdifi(4)=1.     !remnant excitation probability diffractive jpsi
 330       rexndii(1)=0.5   !remnant excitation probability nondiffractive pion
 331       rexndii(2)=0.81   !remnant excitation probability nondiffractive proton
 332       rexndii(3)=0.65   !remnant excitation probability nondiffractive kaon
 333       rexndii(4)=1.   !remnant excitation probability nondiffractive jpsi
 334 c ... up to here.
 335       xmaxremn=0.6      !maximum momentum fraction allowed to be given for remnant mass determination (very important for multiplicity and xf distri at low energy)
 336       xmaxdiff=0.18      !maximum momentum fraction allowed to be given for diffractif remnant mass determination (very important for multiplicity and xf distri at low energy)
 337       rexres=0.5      !remnant excitation probability in nucleus
 338       zbarfl=1.0      !remnant baryon flow parameter (if negative)
 339       ptfrab=0.30    !string pt for remnant baryon (in conn with zbarfl)
 340       alpdif=0.6666     !alpha mass diffractive for cross section and metropolis
 341       alpdi=0.6666      !alpha mass diffractive
 342       alpndi=1.05       !alpha mass nondiffractive
 343       alpsea=0.3      !alpha string end x for sea parton
 344       alpval=0.3      !alpha string end x for valence parton
 345       alpdiq=0.3      !alpha string end x for sea diquark
 346       ammsqq=0.28     !minimum mass string quark quark
 347       ammsqd=1.08    !minimum mass string quark diquark
 348       ammsdd=1.88    !minimum mass string diquark diquark
 349       delrex=.01     !mass to be added to the above masses
 350       cumpom=0       !cutoff mass for virtual pomerons
 351       zbcut=1.e-10    !minimum probability for bkmxndif calculation (for Z_p/t)
 352       alpdro(1)=1.25   !factor for minimum mass of leading droplet
 353       alpdro(2)=0.6   !pt of leading droplet
 354       alpdro(3)=alpndi  !alpha mass of leading droplet
 355
 356       iodiba=0.      ! if iodiba=1, study H-Dibaryon (not used (see ProRef in epos-ems ????)
 357       bidiba=0.030   ! epsilon of H-Dibaryon
 358
 359 c screening splitting +++
 360
 361       epscrw=0.18      !overall factor for Z (Zsame,Zother)-> pp xsect .... w_Z
 362       epscrv=0.0       !low energy value for Zother   ..................... w_M
 363       epscrp=0.75       !b width param for Z     -> pp xsection ............ w_B
 364       egyscr=10.       !screening minimum energy -> pp xsection ........... s_M
 365       epscrs=0.02       !screening power increase soft  -> pp xsctn ........ alp_S
 366       epscrh=1.75      !screening power increase hard  -> pp xsctn ........ alp_H
 367       epscrx=0.35       !screening power maxi           -> pp xsctn
 368       zidinc= 0.2       !inner diquark modif factor increase
 369       zidmax= 0.6       !inner diquark modif factor max
 370       zodinc= 0.      !outer diquark modif factor increase
 371       zodmax= 0.       !outer diquark modif factor max
 372       zosinc= 0.      !outer strange modif factor increase
 373       zosmax= 0.       !outer strange modif factor max
 374       zipinc= 0.4    !inner pt modif factor increase
 375       zipmax= 0.9      !inner modif factor max
 376       zopinc= 0.07    !outer pt and end pt modif factor increase
 377       zopmax= 1.      !outer pt and end pt modif factor max
 378       zoeinc=0.44         !cutoff timelike evol modif factor increase
 379       zoemax=4.      !cutoff timelike evol modif factor max
 380
 381 c Reggeon parameters  (not used)
 382
 383       alpreg=0.734   !alpha_reggeon
 384       sloreg=0.499   !slope_reggeon
 385       gamreg=16.46   !gamma_reggeon
 386       r2reg=0.613    !r^2_reggeon
 387
 388 c  masses
 389
 390       amhadr(1)=.14            !pion mass
 391       amhadr(2)=.939           !nucleon mass
 392       amhadr(3)=.496           !kaon mass
 393       amhadr(4)=1.868          !d-meson mass
 394       amhadr(6)=1.116          !lambda mass
 395       amhadr(5)=2.27           !lambda_c mass
 396       amhadr(7)=.548           !eta mass
 397       amhadr(8)=3.5            !J/psi mass
 398       qcmass=1.6               !c-quark mass
 399 c      amhdibar=2.200           !h-dibaryon mass       !not used any more
 400       qmass(0)=0.980           !diquark effective bounding energy (for pt distribtions)
 401       isospin(0)=0
 402       qmass(1)=0.001           !u quark effective mass (for pt distribtions)
 403       isospin(1)=1
 404       qmass(2)=0.003           !d quark effective mass (for pt distribtions)
 405       isospin(2)=-1
 406       qmass(3)=0.140           !s quark effective mass (for pt distribtions)
 407       isospin(3)=0
 408       qmass(4)=qcmass          !c quark effective mass (for pt distribtions)
 409       isospin(4)=0
 410       call idmass(5,qbmass)
 411       qmass(5)=qbmass          !b quark effective mass (for pt distribtions)
 412       isospin(5)=0
 413       call idmass(6,qtmass)
 414       qmass(6)=qtmass          !t quark effective mass (for pt distribtions)
 415       isospin(6)=0
 416
 417 c  nucleus-nucleus
 418
 419       iokoll=0      !fix # of collisions (1)
 420       laproj=0      !projectile charge number
 421       maproj=0      !projectile mass number
 422       latarg=0      !target charge number
 423       matarg=0      !target mass number
 424       core=0.34     !hard core distance(2*0.17)
 425 c      ncolmx=100000 !maximum number of collisions
 426       fctrmx=10     !parameter determining range for density distribution
 427       bmaxim=10000  !maximum impact parameter
 428       bminim=0.     !minimum impact parameter
 429       phimax=2*3.1415927 !maximum phi
 430       phimin=0      !minimum phi
 431       ionudi=0      !nuclear diffraction included (1) or not (0)
 432
 433 c rescattering parameters +++
 434
 435       iorsce=0      !color exchange turned on(1) or off(0)
 436       iorsdf=3      !droplet formation turned on(>0) or off(0)
 437       iorshh=0      !other hadron-hadron int. turned on(1) or off(0)
 438       iocluin=1     !include inwards moving corona segments into clusters (1) or not (0)
 439       ioquen=0      !jet quenching option (0=no)
 440       taumin=1.0    !min tau for rescattering
 441       taurea=0.001    !reaction time (=formation time)
 442       nsegsu=30     !number of segments per subcluster
 443       nsegce=4      !number of segments per cell
 444       ptclu=3  !obove this: string segments escape cluster
 445       amimfs=1.0    !below this: elastic
 446       amimel=0.050  !below this: nothing
 447       delamf=1      !above this: color exch  !cutoffs for kinetic energy
 448       deuamf=1      !above this: nothing     !mass - minimum mass of pair
 449       epscri(1)=0.22   !critical energy density
 450       epscri(3)=0.08   !critical energy density
 451
 452 c rescattering parameters
 453
 454       amsiac=.8     !interaction mass
 455       amprif=0      !print option
 456       delvol=1      !print option
 457       deleps=1      !print option
 458       deltau=0.2    !initial delta_tau for space-time evolution
 459       factau=1.05   !factor for delta_tau for space-time evolution
 460       numtau=80     !number of tau steps for space-time evolution
 461       dlzeta=0.5    !delta_zeta for longitudinal droplet binning
 462       etafac=1.0    !factor determining inner range
 463       facnuc=0.0    !factor for nuclear size to determine inner range
 464       hacore=0.8    !hadron compressibility factor
 465       cepara=0.03   !parameter for excitation for color exchange
 466       dscale=1.    !scale parameter for hadron-hadron
 467       iceopt=1      !real color exchange (1) or just excitation (0)
 468
 469 c coupling with urqmd
 470
 471       iurqmd=0    ! call eposurqmd (1) or not
 472
 473 c initial conditions for hydro
 474
 475       ispherio=0    !call spherio
 476       cutico=3      !cutoff parameter for smoothing kernel in epos-ico
 477                     !  (as small as possible with stable results)
 478       dssico=0.2    !s step size for string integration in epos-ico
 479                     !  (as big as possible with stable results)
 480       icocore=0     !consider core initial condition (1 or 2)
 481       icotabm=0     !make table for initial condition stuff
 482       icotabr=0     !read table for initial condition stuff
 483
 484 c  cluster decay
 485
 486       ioclude=1     !cluster decay option
 487       amuseg=5.0    !min mass for radial boost
 488       yradmx= 0.77   !max radial collective boost
 489       yradmi= 0.15   !max radial collective boost increase
 490       yradpp=1.2    !max radial collective boost pp basis
 491       yradpi=0.4    !max radial collective boost pp increase
 492       yradpx=1.2    !max radial collective boost pp maxi
 493       facecc=0.45   !eccentricity parameter
 494       ylongmx=-1.   !max long collective boost ( < 0 -> take from jintpo )
 495       rcoll=0.0     !radial collective flow param
 496       bag4rt=0.200  !bag constant ^1/4
 497       taunll=1.0    !decay time (comoving frame)
 498       vrad=0.3
 499       facts=0.3       !gamma-s factor
 500       factb=1       !gamma-s factor baryons
 501       factq=1       !gamma-qqbar
 502       inbxxx=0      !1: fast nboby via hnbxxx   ????? maybe obsolete???????
 503
 504 c  droplet decay initializations
 505
 506          asuhax(1)=1.134  !two lowest multiplets
 507          asuhax(2)=1.301  !two lowest multiplets
 508          asuhax(3)=1.461  !two lowest multiplets
 509          asuhax(4)=1.673  !two lowest multiplets
 510          asuhax(5)=0.7700 !two lowest multiplets   rho
 511          asuhax(6)=0.8920 !two lowest multiplets   K*
 512          asuhax(7)=1.2320 !two lowest multiplets
 513          asuhay(1)=0.940  !lowest multiplet
 514          asuhay(2)=1.200  !lowest multiplet
 515          asuhay(3)=1.322  !lowest multiplet
 516          asuhay(4)=1.673  !lowest multiplet
 517          asuhay(5)=0.1400 !lowest multiplet
 518          asuhay(6)=0.4977 !lowest multiplet
 519          asuhay(7)=1.2320 !lowest multiplet
 520
 521 c  droplet specification
 522
 523       keu=0     !u flavour of droplet
 524       ked=0     !d flavour of droplet
 525       kes=0     !s flavour of droplet
 526       kec=0     !c flavour of droplet
 527       keb=0     !b flavour of droplet
 528       ket=0     !t flavour of droplet
 529       tecm=10   !energy of droplet
 530       volu=70   !volume of droplet
 531
 532 c  metropolis and grand canonical
 533
 534       iospec=10 !option for particle species
 535       iocova=1  !covariant (1) or noncovariant (2) phase space integral
 536       iopair=2  !single pair (1) or double pair (2) method
 537       iozero=65 !relative weight of zeros (compared to hadrons)
 538                 ! (-1) nspecs
 539                 ! (-2) nspecs/sqrt(tecm/volu)
 540       ioflac=1  !test multipl distr without (1) or with (2) flavour conserv
 541                 !  (2 only good for nspecs=3,7)
 542       iostat=1  !use boltzmann (1) or quantum (0) statistics in hgc-routines
 543       ioinco=1  !call hnbmin for initial configuration (0)
 544                 !call hgcnbi for initial configuration to generate better
 545                 !initial configuration (1)
 546                 !call hgcnbi for initial configuration to generate optimal
 547                 !initial configuration (2)
 548       iograc=1  !call hgcaaa in case of ioinco=0 (1)
 549       epsgc=2.  !required accuracy in hgcaaa 10**(-epsgc)
 550       iocite=0  !perform counting at metropolis iterations (1) or not (else)
 551       ioceau=0  !perform counting for exp. autocorrel. time (1) or not (else)
 552       iociau=0  !perform counting for int. autocorrel. time (1) or not (else)
 553       ioinct=0  !test grand canonical metropolis sampling (1)
 554                 !to plot results call xhgccc, xhgcfl and xhgcam
 555       ioinfl=1  !conserve flavor in initial configuration in hgcnbi (1)
 556                 !do not conserve flavor (0)
 557                 !do not conserve flavor and energy (-1)
 558       iowidn=2  !width of total multiplicity distribution in hgcnbi
 559                 ! sigma_tot -> sigma_tot/iowidn
 560                 ! >0 unnormalized
 561                 ! <0 normalized
 562       ionlat=2  !determine nlattc ,old method (0)
 563                 !or determine nlattc in hgcnbi as:
 564                 ! (1) max(1.3*<N>,<N>+2*sigma,6)
 565                 ! (2) max(1.5*<N>,<N>+3*sigma,6)
 566                 ! (3) max(2.0*<N>,<N>+4*sigma,6)
 567       iomom=1   !number of momenta to be changed in hnbodz
 568       ioobsv=0  !observable for autocorrelation time calculation
 569                 !0: total multiplicity
 570                 !else: particle id for particle species
 571       iosngl=0  !event # for which counting at metropolis iterations is done
 572       iorejz=0  !reject pair exchange with only zeros (1) or not (else)
 573       iompar=4  !parameter for windowing algorithm
 574       iozinc=0  !if iozevt>0: modifies iozero for testing (in sr hgcnbi)
 575       iozevt=0  !if >0: modifies iozero for testing (in sr hgcnbi)
 576       nadd=0    !number of pi0s added to minimum initial configuration
 577       iterma=-6 !>0: maximum number of iterations
 578                 !<0: - number of iterations per corr time
 579       iterpr=10000 !iter-increment for printout
 580       iterpl=1  !iter-increment for plot
 581       iternc=50 !iterations not counted for analysis
 582       epsr=1e-4 !required accuracy in sr hnbraw
 583       keepr=1   !keep most random numbers rnoz in hnbodz (1)
 584                 !  or update all (0)
 585
 586 c  strangelets
 587
 588       iopenu=1      !option for minimum energy
 589                     !1: sum of hadron masses
 590                     !2: bag model curve with minimum at nonzero strangen
 591       themas=.51225 !parameter theta in berger/jaffe mass formula
 592
 593 c tests
 594
 595       iotst1=0     !test
 596       iotst2=0     !test
 597       iotst3=0     !test
 598       iotst4=0     !test
 599
 600 c  jpsi
 601
 602       jpsi=0     !jpsi to be produced (1) or not (0)
 603       jpsifi=0   !jpsi final state interaction (1) or not (0)
 604       sigj=0.2   !jpsi nucleon cross section [fm**2]
 605       taumx=20   !max time for jpsi evolution
 606       nsttau=100 !time steps for jpsi evolution
 607       ijphis=0   !fill jpsi histograms (1) or not (0)
 608
 609 c  analysis: intermittency, space-time, droplets, formation time
 610
 611       ymximi=2   !upper limit for rapidity interval for intermittency analysis
 612       imihis=0   !fill intermittency histograms (1) or not (0)
 613       isphis=0   !fill space-time histograms (1) or not (0)
 614       iologb=0   !0=linear bins   1=log bins
 615       ispall=1   !xspace: all ptls (1) or only interacting ptls (else)
 616       wtmini=-3  !tmin in xspace
 617       wtstep=1   !t-step in xspace
 618       iwcent=0   !only central point (1) or longitudinal distr (else) in xspace
 619       iclhis=0   !fill droplet histograms (1) or not (0)
 620       iwtime=0   !fill formation time histogram (1) or not (else)
 621       wtimet=100 !max time in wtime
 622       wtimei=0   !max mass in wtime
 623       wtimea=1000 !max mass in wtime
 624
 625
 626 c  storing
 627       maxrec(1)=7
 628       irecty(1,1)=1
 629       irecty(2,1)=2
 630       irecty(3,1)=3
 631       irecty(4,1)=4
 632       irecty(5,1)=5
 633       irecty(6,1)=6
 634       irecty(7,1)=7
 635       maxrec(2)=14
 636       irecty(1,2)=1
 637       irecty(2,2)=2
 638       irecty(3,2)=3
 639       irecty(4,2)=4
 640       irecty(5,2)=5
 641       irecty(6,2)=6
 642       irecty(7,2)=7
 643       irecty(8,2)=8
 644       irecty(9,2)=9
 645       irecty(10,2)=10
 646       irecty(11,2)=11
 647       irecty(12,2)=12
 648       irecty(13,2)=13
 649       irecty(14,2)=14
 650
 651
 652 c  other
 653
 654       gaumx=8    !range for gauss distribution
 655       nclean=1   !clean /cptl/ if nclean > 0
 656                  !(not for part with istptl<istmax if nclean=1 (do not change analysis)
 657                  ! for all part with ist.ne.0 if nclean > 1 (particle nb reduce at max))
 658       istore=0   !0: no storage to data-file
 659                  !1: epos     standard
 660                  !2: OSC1997A standart
 661                  !3: OSC1999A standart
 662                  !5: calls routine ustore (modifiable by the user)
 663       ioidch=1   !id choice for storage to data-file
 664       iframe=0   !frame specification production run
 665       jframe=0   !frame specification analysis
 666       kframe=0   !frame specification analysis (2nd frame)
 667                  ! 1:total
 668                  !11:nucleon-nucleon
 669                  !12:target
 670                  !21:gamma-nucleon
 671                  !22:lab
 672                  !32:sphericity
 673                  !33:thrust
 674       irescl=1   !momentum rescaling (1) or not (0)
 675       ifrade=1   !suppression of fragmentation and decay (0)
 676       idecay=1   !suppression of decay (0)
 677       jdecay=1   !suppression of cluster decay (0), concerns only ity=60 cluster
 678       ntrymx=10  !try-again parameter
 679       istmax=1   !analyse only istptl <= istmax
 680       irdmpr=0   !random sign for projectile if 1
 681
 682 c  constants
 683
 684       pi=3.1415927
 685       pii=1./(2*pi**2)
 686       hquer=0.197327
 687       prom=.94
 688       piom=.14
 689       ainfin=1e31
 690
 691 c air
 692
 693       airanxs(1)=14.007
 694       airznxs(1)=7.
 695       airwnxs(1)=0.781
 696       airanxs(2)=15.999
 697       airznxs(2)=8.
 698       airwnxs(2)=.21
 699       airanxs(3)=39.948
 700       airznxs(3)=18.
 701       airwnxs(3)=0.009
 702       airavanxs=airanxs(1)*airwnxs(1)+airanxs(2)*airwnxs(2)
 703      &         +airanxs(3)*airwnxs(3)
 704       airavznxs=airznxs(1)*airwnxs(1)+airznxs(2)*airwnxs(2)
 705      &         +airznxs(3)*airwnxs(3)
 706
 707 c  zero initializations
 708
 709       ixgeometry=0
 710       ixtau=0
 711       iEmsB=0
 712       iEmsBg=0
 713       iEmsPm=0
 714       iEmsPx=0
 715       iEmsSe=0
 716       iEmsDr=0
 717       iEmsRx=0
 718       iEmsI2=0
 719       iEmsI1=0
 720       iSpaceTime=0
 721       nemsi=0
 722       facmss=1.
 723       nstmax=0
 724       do 6 i=1,99
 725       prob(i)=0
 726       do 6 j=1,2
 727       icbac(i,j)=0
 728 6     icfor(i,j)=0
 729       imsg=0
 730       do j=1,mxjerr
 731        jerr(j)=0
 732       enddo
 733       ntevt=0
 734       nrevt=0
 735       naevt=0
 736       nrstr=0
 737       nrptl=0
 738       nptlu=0
 739       do itau=1,mxtau
 740       volsum(itau)=0
 741       vo2sum(itau)=0
 742       nclsum(itau)=0
 743       do ivol=1,mxvol
 744       do ieps=1,mxeps
 745       clust(itau,ivol,ieps)=0
 746       enddo
 747       enddo
 748       enddo
 749       iutotc=0
 750       iutote=0
 751       nopen=0
 752       nopenr=0
 753       knxopen=0
 754       kchopen=0
 755       khiopen=0
 756       kdtopen=0
 757       klgopen=0
 758       ifdat=0
 759       ifncs=0
 760       xpar1=0.
 761       xpar2=0.
 762       xpar3=0.
 763       xpar4=0.
 764       xpar5=0.
 765       xpar6=0.
 766       xpar7=0.
 767       xpar8=0.
 768       if(iop.eq.0)khisto=0
 769       nrclu=0
 770       nrnody=0
 771       do n=1,mxnody
 772       nody(n)=0
 773       enddo
 774       nrpri=0
 775       do n=1,mxpri
 776       subpri(n)='      '
 777       ishpri(n)=0
 778       enddo
 779       nctcor=0
 780       ncttim=0
 781       do n=1,matau
 782       tauv(n)=0
 783       enddo
 784       ncnt=0
 785       nrnucl(1)=0
 786       nrnucl(2)=0
 787       do i=1,mxnucl
 788       rnucl(i,1)=0
 789       rnucl(i,2)=0
 790       rnuclo(i,1)=0
 791       rnuclo(i,2)=0
 792       enddo
 793       inicnt=0
 794       accept=0.
 795       reject=0.
 796       do n=1,matau
 797       tauv(n)=0.
 798       enddo
 799       anquasiel=0.
 800       nglacc=0
 801       ifset=1
 802
 803  1002 continue ! ----> entry iop=2 <----
 804
 805
 806 c  analysis
 807
 808       xvaria='numptl'
 809       yvaria='ycmptl'
 810       normal=11
 811       xminim=-100
 812       xmaxim=100
 813       nrbins=100
 814       hisfac=1
 815       do nr=1,mxbins
 816       do l=1,5
 817       ar(nr,l)=0
 818       enddo
 819       enddo
 820       nozero=0
 821       ibmin=1
 822       ibmax=1e8
 823
 824       return
 825       end
 826
 827 c-----------------------------------------------------------------------
 828       subroutine ustore
 829 c-----------------------------------------------------------------------
 830 c     writes the results of a simulation into the file with unit ifdt
 831 c     contains a description of the stored variables.
 832 c     modifiable by the user
 833 c-----------------------------------------------------------------------
 834       include 'epos.inc'
 835
 836 c  count the number of particles to be stored (--> nptevt)
 837
 838       nptevt=0
 839       do i=1,nptl
 840       if(istptl(i).le.istmax)nptevt=nptevt+1
 841       enddo
 842
 843 c  store event variables:
 844
 845       write(ifdt,*)nrevt,nptevt,bimevt,phievt,kolevt
 846      *,pmxevt,egyevt,npjevt,ntgevt
 847
 848 c     nrevt.......... event number
 849 c     nptevt ........ number of (stored!) particles per event
 850 c     bimevt ........ absolute value of impact parameter
 851 c     phievt ........ angle of impact parameter
 852 c     kolevt ........ number of collisions
 853 c     pmxevt ........ reference momentum
 854 c     egyevt ........ pp cm energy (hadron) or string energy (lepton)
 855 c     npjevt ........ number of primary projectile participants
 856 c     ntgevt ........ number of primary target participants
 857 c     npnevt ........ number of primary projectile neutron spectators
 858 c     nppevt ........ number of primary projectile proton spectators
 859 c     ntnevt ........ number of primary target neutron spectators
 860 c     ntpevt ........ number of primary target proton spectators
 861 c     jpnevt ........ number of absolute projectile neutron spectators
 862 c     jppevt ........ number of absolute projectile proton spectators
 863 c     jtnevt ........ number of absolute target neutron spectators
 864 c     jtpevt ........ number of absolute target proton spectators
 865
 866       do i=1,nptl
 867
 868       if(istptl(i).le.istmax)then !store events with istptl < istmax
 869
 870 c  store particle variables:
 871
 872       write(ifdt,*)
 873      *i,idptl(i),pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i),pptl(5,i)
 874      *,iorptl(i),jorptl(i),istptl(i)
 875      *,xorptl(1,i),xorptl(2,i),xorptl(3,i),xorptl(4,i)
 876      *,tivptl(1,i),tivptl(2,i)
 877
 878 c     i ............. particle number
 879 c     idptl(i) ...... particle id
 880 c     pptl(1,i) ..... x-component of particle momentum
 881 c     pptl(2,i) ..... y-component of particle momentum
 882 c     pptl(3,i) ..... z-component of particle momentum
 883 c     pptl(4,i) ..... particle energy
 884 c     pptl(5,i) ..... particle mass
 885 c     iorptl(i) ..... particle number of father (if .le. 0 : no father)
 886 c     jorptl(i) ..... particle number of mother (if .le. 0 : no mother)
 887 c     istptl(i) ..... generation flag: last gen. (0) or not (1)
 888 c     ityptl(i) ..... particle type (string, remnant ...)
 889 c     xorptl(1,i) ... x-component of formation point
 890 c     xorptl(2,i) ... y-component of formation point
 891 c     xorptl(3,i) ... z-component of formation point
 892 c     xorptl(4,i) ... formation time
 893 c     tivptl(1,i) ... formation time (always in the pp-cms!)
 894 c     tivptl(2,i) ... destruction time (always in the pp-cms!)
 895
 896       endif
 897       enddo
 898
 899       return
 900       end
 901
 902
 903 c-----------------------------------------------------------------------
 904       subroutine bstora
 905 c-----------------------------------------------------------------------
 906 c     writes the results of a simulation into the file with unit ifdt
 907 c     contains a description of the stored variables.
 908 c-----------------------------------------------------------------------
 909       include 'epos.inc'
 910       common/record/maxrec(2),irecty(30,2)
 911       character code*8,version*8,frame*4
 912       code='epos   '
 913       write(version,'(f5.2,3x)')iversn/100.
 914       if(iframe.eq. 1)frame='ttcm'
 915       if(iframe.eq.11)frame='nncm'
 916       if(iframe.eq.12)frame='targ'
 917       if(iframe.eq.21)frame='gncm'
 918       if(iframe.eq.22)frame='lncm'
 919       ntest=1
 920       if (istore.eq.2) then     ! OSC1997A
 921         write (ifdt,'(a)') 'OSC1997A'
 922         write (ifdt,'(a)') 'final_id_p_x'
 923         write(ifdt,100) code,version
 924      *       ,maproj,laproj,matarg,latarg,frame,engy,ntest
 925  100    format(2(a8,'  '),'(',i3,',',i3,')+(',i3,',',i3,')',
 926      *       '  ',a4,'  ',e10.4,'  ',i8)
 927         maxrec(1)=4
 928         irecty(1,1)=1           !nevt
 929         irecty(2,1)=2
 930         irecty(3,1)=3
 931         irecty(4,1)=4
 932         maxrec(2)=11
 933         irecty(1,2)=1           !nr
 934         irecty(2,2)=2           !id
 935         irecty(3,2)=3           !px
 936         irecty(4,2)=4           !py
 937         irecty(5,2)=5           !pz
 938         irecty(6,2)=6           !E
 939         irecty(7,2)=7           !M
 940         irecty(8,2)=11          !x
 941         irecty(9,2)=12          !y
 942         irecty(10,2)=13         !z
 943         irecty(11,2)=14         !t
 944       elseif(istore.eq.3) then
 945  201    format('# ',a)
 946         write (ifdt,201) 'OSC1999A'
 947         if(istmax.eq.0) then
 948           write (ifdt,201) 'final_id_p_x'
 949         elseif(istmax.ge.2) then
 950           write (ifdt,201) 'full_event_history'
 951         endif
 952  202    format('# ',a8,' ',a8)
 953         write(ifdt,202) code,version !3rd line
 954  203    format('# (',i3,',',i3,')+(',i3,',',i3,')',
 955      *       '  ',a4,'  ',e10.4,'  ',i8)
 956         write(ifdt,203) maproj,laproj,matarg,latarg,frame,engy,ntest
 957         maxrec(1)=5
 958         irecty(1,1)=2           !nevt
 959         irecty(2,1)=0           !zero
 960         irecty(3,1)=1           !additional information
 961         irecty(4,1)=3           !additional information
 962         irecty(5,1)=4           !additional information
 963         maxrec(2)=12
 964         irecty(1,2)=1           !nr
 965         irecty(2,2)=2           !id
 966         irecty(3,2)=10          !ist
 967         irecty(4,2)=3           !px
 968         irecty(5,2)=4           !py
 969         irecty(6,2)=5           !pz
 970         irecty(7,2)=6           !E
 971         irecty(8,2)=7           !M
 972         irecty(9,2)=11          !x
 973         irecty(10,2)=12         !y
 974         irecty(11,2)=13         !z
 975         irecty(12,2)=14         !t
 976                                 ! nin nout [optional information]
 977                                 ! ipart id ist px py pz p0 mass x y z t [optional information]
 978
 979       endif
 980
 981       end
 982
 983 c-----------------------------------------------------------------------
 984       subroutine bstore
 985 c-----------------------------------------------------------------------
 986 c     writes the results of a simulation into the file with unit ifdt
 987 c     contains a description of the stored variables.
 988 c-----------------------------------------------------------------------
 989
 990       include 'epos.inc'
 991       common/record/maxrec(2),irecty(30,2)
 992       common/dimensi/k2(100)
 993
 994       nptevt=0
 995       do n=1,nptl
 996         iok=1               !idcode simple
 997         if(istptl(n).gt.istmax) then
 998           iok=0
 999         endif
1000       if (iok.eq.1) nptevt=nptevt+1
1001       enddo
1002  11   format (i6,' ',$)
1003  12   format (e12.6,' ',$)
1004       do i=1,maxrec(1)
1005         l=irecty(i,1)
1006         if(l.eq.0)write(ifdt,21) 0
1007         if(l.eq.1)write(ifdt,11)nrevt
1008         if(l.eq.2)write(ifdt,11)nptevt
1009         if(l.eq.3)write(ifdt,12)bimevt
1010         if(l.eq.4)write(ifdt,12)phievt
1011         if(l.eq.5)write(ifdt,11)kolevt
1012         if(l.eq.6)write(ifdt,12)pmxevt
1013         if(l.eq.7)write(ifdt,12)egyevt
1014         if(l.eq.8)write(ifdt,11)npjevt
1015         if(l.eq.9)write(ifdt,11)ntgevt
1016         if(l.eq.10)write(ifdt,11)npnevt
1017         if(l.eq.11)write(ifdt,11)nppevt
1018         if(l.eq.12)write(ifdt,11)ntnevt
1019         if(l.eq.13)write(ifdt,11)ntpevt
1020         if(l.eq.14)write(ifdt,11)jpnevt
1021         if(l.eq.15)write(ifdt,11)jppevt
1022         if(l.eq.16)write(ifdt,11)jtnevt
1023         if(l.eq.17)write(ifdt,11)jtpevt
1024         if(l.eq.20)write(ifdt,12)amproj
1025         if(l.eq.21)write(ifdt,12)amtarg
1026       enddo
1027       write (ifdt,*)            !RETURN
1028  21   format (i6,' ',$)
1029  22   format (e12.6,' ',$)
1030  23   format (i10,' ',$)
1031       do n=1,nptl
1032         iok=1                   !idcode simple
1033         if(istptl(n).gt.istmax) then
1034           iok=0
1035         endif
1036         if (iok.eq.1) then
1037           id=idptl(n)
1038           if(istore.eq.2.or.istore.eq.3.or.ioidch.eq.2)then
1039             id=idtrafo('nxs','pdg',idptl(n))
1040           endif
1041           do i=1,maxrec(2)
1042             l=irecty(i,2)
1043             if(l.eq.0)write(ifdt,21) 0
1044             if(l.eq.1)write(ifdt,21) n
1045             if(l.eq.2)write(ifdt,23) id
1046             if(l.eq.3.or.l.eq.17)write(ifdt,22) pptl(1,n)
1047             if(l.eq.4.or.l.eq.17)write(ifdt,22) pptl(2,n)
1048             if(l.eq.5.or.l.eq.17)write(ifdt,22) pptl(3,n)
1049             if(l.eq.6.or.l.eq.17)write(ifdt,22) pptl(4,n)
1050             if(l.eq.7.or.l.eq.17)write(ifdt,22) pptl(5,n)
1051             if(l.eq.8)write(ifdt,21) iorptl(n)
1052             if(l.eq.9)write(ifdt,21) jorptl(n)
1053             if(l.eq.10)write(ifdt,21) istptl(n)
1054             if(l.eq.11.or.l.eq.18)write(ifdt,22) xorptl(1,n)
1055             if(l.eq.12.or.l.eq.18)write(ifdt,22) xorptl(2,n)
1056             if(l.eq.13.or.l.eq.18)write(ifdt,22) xorptl(3,n)
1057             if(l.eq.14.or.l.eq.18)write(ifdt,22) xorptl(4,n)
1058             if(l.eq.19)write(ifdt,22) dezptl(n)
1059             if(l.eq.21)write(ifdt,21) ifrptl(1,n)
1060             if(l.eq.22)write(ifdt,21) ifrptl(2,n)
1061             if(l.eq.23)write(ifdt,21) ityptl(n)
1062             if(l.eq.15) then
1063               if(iorptl(n).gt.0)then
1064                 write(ifdt,23) idptl(iorptl(n))
1065               else
1066                 write(ifdt,23) 0
1067               endif
1068             endif
1069             if(l.eq.16) then
1070               if(jorptl(n).gt.0)then
1071                 write(ifdt,23) idptl(jorptl(n))
1072               else
1073                 write(ifdt,23) 0
1074               endif
1075             endif
1076           enddo
1077           write (ifdt,*)        !RETURN
1078         endif
1079       enddo
1080  9999 return
1081       end
1082
1083
1084 c-----------------------------------------------------------------------
1085       subroutine bread
1086 c-----------------------------------------------------------------------
1087 c     writes the results of a simulation into the file with unit ifdt
1088 c     contains a description of the stored variables.
1089 c-----------------------------------------------------------------------
1090
1091       include 'epos.inc'
1092       common/record/maxrec(2),irecty(30,2)
1093       character*255 line
1094       dimension inptl(mxptl)
1095       logical info
1096
1097       info=.false.
1098       do n=1,mxptl
1099         inptl(n)=0
1100       enddo
1101       read(ifdt,'(a255)',end=1)line
1102  1    k=1
1103       nptevt=0
1104       do i=1,maxrec(1)
1105         l=irecty(i,1)
1106         if(l.eq.0)read(line(k:),'(i6)')ldummy !0
1107         if(l.eq.0) k=k+7
1108         if(l.eq.1)read(line(k:),'(i6)')ldummy !nrevt
1109         if(l.eq.1) k=k+7
1110         if(l.eq.2)read(line(k:),'(i6)')nptevt
1111         if(l.eq.2) k=k+7
1112         if(l.eq.3)read(line(k:),'(e12.6)')bimevt
1113         if(l.eq.3) k=k+13
1114         if(l.eq.4)read(line(k:),'(e12.6)')phievt
1115         if(l.eq.4) k=k+13
1116         if(l.eq.5)read(line(k:),'(i6)')kolevt
1117         if(l.eq.5) k=k+7
1118         if(l.eq.6)read(line(k:),'(e12.6)')pmxevt
1119         if(l.eq.6) k=k+13
1120         if(l.eq.7)read(line(k:),'(e12.6)')egyevt
1121         if(l.eq.7) k=k+13
1122         if(l.eq.8)read(line(k:),'(i6)')npjevt
1123         if(l.eq.8) k=k+7
1124         if(l.eq.9)read(line(k:),'(i6)')ntgevt
1125         if(l.eq.9) k=k+7
1126         if(l.eq.10)read(line(k:),'(i6)')npnevt
1127         if(l.eq.10) k=k+7
1128         if(l.eq.11)read(line(k:),'(i6)')nppevt
1129         if(l.eq.11) k=k+7
1130         if(l.eq.12)read(line(k:),'(i6)')ntnevt
1131         if(l.eq.12) k=k+7
1132         if(l.eq.13)read(line(k:),'(i6)')ntpevt
1133         if(l.eq.13) k=k+7
1134         if(l.eq.14)read(line(k:),'(i6)')jpnevt
1135         if(l.eq.14) k=k+7
1136         if(l.eq.15)read(line(k:),'(i6)')jppevt
1137         if(l.eq.15) k=k+7
1138         if(l.eq.16)read(line(k:),'(i6)')jtnevt
1139         if(l.eq.16) k=k+7
1140         if(l.eq.17)read(line(k:),'(i6)')jtpevt
1141         if(l.eq.17) k=k+7
1142         if(l.eq.20)read(line(k:),'(e12.6)')amproj
1143         if(l.eq.20) k=k+13
1144         if(l.eq.21)read(line(k:),'(e12.6)')amtarg
1145         if(l.eq.21) k=k+13
1146       enddo
1147       if(nptevt.eq.0)then
1148         print *,'sorry, '
1149         print *,'there is no particle number in the event record  '
1150         stop
1151       endif
1152       do n=1,nptevt
1153         read(ifdt,'(a255)',end=2)line
1154  2      k=1
1155         do i=1,maxrec(2)
1156           l=irecty(i,2)
1157           if(l.eq.0)read(line(k:),'(i6)') ldummy
1158           if(l.eq.0) k=k+7
1159           if(l.eq.1)then
1160             read(line(k:),'(i6)') nidp
1161             if(nidp.gt.0.and.nidp.le.mxptl)then
1162               info=.true.
1163               inptl(nidp)=n
1164             endif
1165             k=k+7
1166           endif
1167           if(l.eq.2)read(line(k:),'(i10)') idptl(n)
1168           if(l.eq.2) k=k+11
1169           if(l.eq.3.or.l.eq.17)read(line(k:),'(e12.6)') pptl(1,n)
1170           if(l.eq.3) k=k+13
1171           if(l.eq.4.or.l.eq.17)read(line(k:),'(e12.6)') pptl(2,n)
1172           if(l.eq.4) k=k+13
1173           if(l.eq.5.or.l.eq.17)read(line(k:),'(e12.6)') pptl(3,n)
1174           if(l.eq.5) k=k+13
1175           if(l.eq.6.or.l.eq.17)read(line(k:),'(e12.6)') pptl(4,n)
1176           if(l.eq.6) k=k+13
1177           if(l.eq.7.or.l.eq.17)read(line(k:),'(e12.6)') pptl(5,n)
1178           if(l.eq.7) k=k+13
1179           if(l.eq.8)then
1180             read(line(k:),'(i6)') iorptl(n)
1181             k=k+7
1182             if(info.and.iorptl(n).gt.0)iorptl(n)=inptl(iorptl(n))
1183           endif
1184           if(l.eq.9)then
1185             read(line(k:),'(i6)') jorptl(n)
1186             k=k+7
1187             if(info.and.jorptl(n).gt.0)jorptl(n)=inptl(jorptl(n))
1188           endif
1189           if(l.eq.10)read(line(k:),'(i6)') istptl(n)
1190           if(l.eq.10) k=k+7
1191           if(l.eq.11.or.l.eq.18)read(line(k:),'(e12.6)')xorptl(1,n)
1192           if(l.eq.11) k=k+13
1193           if(l.eq.12.or.l.eq.18)read(line(k:),'(e12.6)')xorptl(2,n)
1194           if(l.eq.12) k=k+13
1195           if(l.eq.13.or.l.eq.18)read(line(k:),'(e12.6)')xorptl(3,n)
1196           if(l.eq.13) k=k+13
1197           if(l.eq.14.or.l.eq.18)read(line(k:),'(e12.6)')xorptl(4,n)
1198           if(l.eq.14) k=k+13
1199 c     if(i.eq.15)read(line(k:),'(i6)') idiptl(n)
1200           if(l.eq.15) k=k+7
1201 c     if(i.eq.16)read(line(k:),'(i6)') idjptl(n)
1202           if(l.eq.16) k=k+7
1203           if(l.eq.19)read(line(k:),'(e12.6)') dezptl(n)
1204           if(l.eq.19) k=k+13
1205 c          if(l.eq.21)read(line(k:),'(I6)') ifrptl(1,n)
1206           if(l.eq.21) k=k+7
1207 c          if(l.eq.22)read(line(k:),'(I6)') ifrptl(2,n)
1208           if(l.eq.22) k=k+7
1209           if(l.eq.23)read(line(k:),'(I6)') ityptl(n)
1210           if(l.eq.23) k=k+7
1211         enddo
1212       enddo
1213       nptl=nptevt
1214       nevt=1
1215       end
1216
1217 c-----------------------------------------------------------------------
1218       subroutine aafinal
1219 c-----------------------------------------------------------------------
1220 c  * calculates xorptl(j,i), representing formation points.
1221 c    (xorptl(j,i),j=1,4) is the 4-vector representing the space-time of
1222 c    creation of particle i.
1223 c-----------------------------------------------------------------------
1224       include 'epos.inc'
1225       do i=1,nptl
1226         if(idptl(i).ne.0.and.istptl(i).le.1)then
1227           if(    abs(tivptl(1,i)).le.ainfin
1228      .    .and.abs(xorptl(1,i)).le.ainfin
1229      .    .and.abs(xorptl(2,i)).le.ainfin
1230      .    .and.abs(xorptl(3,i)).le.ainfin
1231      .    .and.abs(xorptl(4,i)).le.ainfin
1232      .    .and.pptl(5,i).le.ainfin
1233      .    .and.pptl(4,i).gt.0.)then
1234             if(ish.ge.4)call alistc('afinal&',i,i)
1235             t=tivptl(1,i)
1236             xorptl(1,i)=xorptl(1,i)+pptl(1,i)/pptl(4,i)*(t-xorptl(4,i))
1237             xorptl(2,i)=xorptl(2,i)+pptl(2,i)/pptl(4,i)*(t-xorptl(4,i))
1238             xorptl(3,i)=xorptl(3,i)+pptl(3,i)/pptl(4,i)*(t-xorptl(4,i))
1239             xorptl(4,i)=t
1240           else
1241             if(ish.ge.1)then
1242               if(iorptl(i).gt.0)idior=idptl(iorptl(i))
1243               write(ifmt,'(a)')
1244      .        '*** warning (afinal see check file): '
1245               write(ifch,'(a,i6,i10,i10,i3,1x,7(e7.1,1x))')
1246      .        '*** warning (afinal): ',
1247      .        i,idptl(i),idior,ityptl(i),tivptl(1,i), pptl(4,i)
1248      .        ,pptl(5,i),xorptl(1,i),xorptl(2,i),xorptl(3,i),xorptl(4,i)
1249             endif
1250             tivptl(1,i)=2*ainfin
1251             xorptl(1,i)=2*ainfin
1252             xorptl(2,i)=2*ainfin
1253             xorptl(3,i)=2*ainfin
1254             xorptl(4,i)=2*ainfin
1255           endif
1256         endif
1257       enddo
1258       end
1259
1260 c-----------------------------------------------------------------------
1261       subroutine afinal
1262 c-----------------------------------------------------------------------
1263 c  does some final calculations, to be called before call aasto.
1264 c  * calculates nptlu, the maximum nptl for all events.
1265 c  * in case of mod(iframe,10) .ne. 1, these vectors are transformed
1266 c    (being originally in the "natural frame",
1267 c  * calculates numbers of spectators:
1268 c    npnevt (number of primary proj neutron spectators)
1269 c    nppevt (number of primary proj proton spectators)
1270 c    ntnevt (number of primary targ neutron spectators)
1271 c    ntpevt (number of primary targ proton spectators)
1272 c    jpnevt (number of absolute proj neutron spectators)
1273 c    jppevt (number of absolute proj proton spectators)
1274 c    jtnevt (number of absolute targ neutron spectators)
1275 c    jtpevt (number of absolute targ proton spectators)
1276 c-----------------------------------------------------------------------
1277
1278       include 'epos.inc'
1279       common/geom/rmproj,rmtarg,bmax,bkmx
1280
1281       double precision pp1,pp2,pp3,pp4,pp5
1282       call utpri('afinal',ish,ishini,4)
1283
1284       nptlu=max0(nptl,nptlu)
1285
1286       if(mod(iframe,10).ne.1)then
1287       if(iframe.eq.12)then   !targ
1288       pp1=0d0
1289       pp2=0d0
1290       pp3=dsinh(dble(yhaha))
1291       pp4=dcosh(dble(yhaha))
1292       pp5=1d0
1293       else
1294       stop'transformation not yet defined'
1295       endif
1296       endif
1297
1298       do i=1,nptl
1299         if(idptl(i).ne.0.and.istptl(i).le.1)then
1300           if(    abs(tivptl(1,i)).le.ainfin
1301      .    .and.abs(xorptl(1,i)).le.ainfin
1302      .    .and.abs(xorptl(2,i)).le.ainfin
1303      .    .and.abs(xorptl(3,i)).le.ainfin
1304      .    .and.abs(xorptl(4,i)).le.ainfin
1305      .    .and.pptl(5,i).le.ainfin
1306      .    .and.pptl(4,i).gt.0.)then
1307             if(mod(iframe,10).ne.1)then
1308               if(iframe.eq.12)then
1309                 call utlob4(-1,pp1,pp2,pp3,pp4,pp5
1310      .          ,xorptl(1,i),xorptl(2,i),xorptl(3,i),xorptl(4,i))
1311                 call utlob5(-yhaha
1312      .          ,pptl(1,i), pptl(2,i), pptl(3,i), pptl(4,i), pptl(5,i))
1313               else
1314                 stop'transformation not yet defined'
1315               endif
1316             endif
1317           else
1318             tivptl(1,i)=ainfin
1319             xorptl(1,i)=ainfin
1320             xorptl(2,i)=ainfin
1321             xorptl(3,i)=ainfin
1322             xorptl(4,i)=ainfin
1323           endif
1324         endif
1325       enddo
1326
1327       if(ish.ge.2)then
1328         if(model.eq.1)call alistf('EPOS&')
1329         if(model.eq.2)call alistf('QGSJET01&')
1330         if(model.eq.3)call alistf('GHEISHA&')
1331         if(model.eq.4)call alistf('PYTHIA&')
1332         if(model.eq.5)call alistf('HIJING&')
1333         if(model.eq.6)call alistf('SIBYLL 2.1&')
1334         if(model.eq.7)call alistf('QGSJET II&')
1335       endif
1336
1337 c      if(isto.eq.1)stop
1338 c$$$      call testconex(2)
1339
1340       npnevt=0
1341       nppevt=0
1342       ntnevt=0
1343       ntpevt=0
1344       jpnevt=0
1345       jppevt=0
1346       jtnevt=0
1347       jtpevt=0
1348       if(ish.ge.6)write(ifch,'(/31a1/a/31a1)')('-',l=1,31)
1349      *,'primary and absolute spectators',('-',l=1,31)
1350       if(ish.ge.6)write(ifch,'(/a//a/)')'projectile nucleons:'
1351      *,'     i    id   ior   ist'
1352       do i=1,maproj
1353       if(ish.ge.6)write(ifch,'(4i6)')i,idptl(i),iorptl(i),istptl(i)
1354       io=iorptl(i)
1355       id=idptl(i)
1356       is=istptl(i)
1357       if(io.eq.0.and.id.eq.1220)npnevt=npnevt+1
1358       if(io.eq.0.and.id.eq.1120)nppevt=nppevt+1
1359       if(io.eq.0.and.is.eq.0.and.id.eq.1220)jpnevt=jpnevt+1
1360       if(io.eq.0.and.is.eq.0.and.id.eq.1120)jppevt=jppevt+1
1361       enddo
1362       if(ish.ge.6)write(ifch,'(/a//a/)')'target nucleons:'
1363      *,'     i    id   ior   ist'
1364       do i=maproj+1,maproj+matarg
1365       if(ish.ge.6)write(ifch,'(4i6)')i,idptl(i),iorptl(i),istptl(i)
1366       io=iorptl(i)
1367       id=idptl(i)
1368       is=istptl(i)
1369       if(io.eq.0.and.id.eq.1220)ntnevt=ntnevt+1
1370       if(io.eq.0.and.id.eq.1120)ntpevt=ntpevt+1
1371       if(io.eq.0.and.is.eq.0.and.id.eq.1220)jtnevt=jtnevt+1
1372       if(io.eq.0.and.is.eq.0.and.id.eq.1120)jtpevt=jtpevt+1
1373       enddo
1374       if(ish.ge.6)then
1375       write(ifch,'(/a/)')'numbers of participants and spectators:'
1376       write(ifch,'(a,i4,a,i4)')'primary participants:   projectile:'
1377      *,npjevt,'   target:',ntgevt
1378       write(ifch,'(a,i4,a,i4)')'primary spectators:     projectile:'
1379      *,npnevt+nppevt,'   target:',ntnevt+ntpevt
1380       write(ifch,'(a,i4,a,i4)')
1381      *'primary spectator neutrons:   projectile:',npnevt
1382      *,'   target:',ntnevt
1383       write(ifch,'(a,i4,a,i4)')
1384      *'primary spectator protons:    projectile:',nppevt
1385      *,'   target:',ntpevt
1386       write(ifch,'(a,i4,a,i4)')'absolute spectators:    projectile:'
1387      *,jpnevt+jppevt,'   target:',jtnevt+jtpevt
1388       endif
1389
1390       if(isigma.eq.1.and.ntevt.gt.0)then
1391         b1=bminim
1392         b2=min(bmax,bmaxim)
1393         a=pi*(b2**2-b1**2)
1394         sigineex=anintine/float(ntevt)*a*10
1395         sigdifex=anintdiff/float(ntevt)*a*10
1396         sigsdex=anintsdif/float(ntevt)*a*10
1397       endif
1398
1399
1400       if(ish.eq.27)stop'change this????????????' !call hnbcor(2)
1401       if(imihis.eq.1)call wimi
1402       if(imihis.eq.1.and.nrevt.eq.nevent)call wimino
1403       if(isphis.eq.1)call xspace(1)
1404       if(iclhis.eq.1)call wclu
1405       if(iclhis.eq.1.and.nrevt.eq.nevent)call wclufi
1406       if(iwtime.eq.1)call wtime(1)
1407       if(iwtime.eq.1.and.nrevt.eq.nevent)call wtime(2)
1408
1409       if(ish.ge.8)call alistc('afinal&',1,nptl)
1410
1411       call utprix('afinal',ish,ishini,4)
1412       return
1413       end
1414
1415 c-----------------------------------------------------------------------
1416       subroutine bfinal
1417 c-----------------------------------------------------------------------
1418       include 'epos.inc'
1419       if(jerr(1).gt.4.and.jerr(1).gt.0.01*nevent)then
1420         write(ifch,'(3x,70a1)')('#',i=1,70)
1421         write(ifch,*)'  #   number of events:',nevent
1422         write(ifch,*)'  #   number of (flav > 9) warnings:',jerr(1)
1423         write(ifch,*)'  #        (OK when happens rarely)'
1424         write(ifch,'(3x,70a1)')('#',i=1,70)
1425       endif
1426       if(jerr(3).gt.4.and.jerr(3).gt.0.01*jerr(2))then
1427         write(ifch,'(3x,70a1)')('#',i=1,70)
1428         write(ifch,*)'  #   number of clusters:',jerr(2)
1429         write(ifch,*)'  #   number of neg m^2 clusters:',jerr(3)
1430         write(ifch,*)'  #          (OK when happens rarely)'
1431         write(ifch,'(3x,70a1)')('#',i=1,70)
1432       endif
1433       if(jerr(5).gt.4.and.jerr(5).gt.0.01*jerr(4))then
1434         write(ifch,'(3x,70a1)')('#',i=1,70)
1435         write(ifch,*)'  #   number of successful remnant cluster'
1436      &                       ,' decays:',jerr(4)
1437         write(ifch,*)'  #   number of unsuccessful remnant cluster'
1438      &                       ,' decays:',jerr(5)
1439         write(ifch,*)'  #        (OK when happens rarely)'
1440         write(ifch,'(3x,70a1)')('#',i=1,70)
1441       endif
1442       if(jerr(6).gt.4.and.jerr(6).gt.0.01*nevent)then
1443         write(ifch,'(3x,70a1)')('#',i=1,70)
1444         write(ifch,*)'  #   number of events:',nevent
1445         write(ifch,*)'  #   number of low mass remnant clusters:'
1446      &                      ,jerr(6)
1447         write(ifch,*)'  #        (OK when happens rarely)'
1448         write(ifch,'(3x,70a)')('#',i=1,70)
1449       endif
1450       end
1451
1452 c-----------------------------------------------------------------------
1453       subroutine ainit
1454 c-----------------------------------------------------------------------
1455       include 'epos.inc'
1456       include 'epos.incsem'
1457       include 'epos.incpar'
1458       common/cquama/quama
1459       parameter (nptj=129)
1460       common /cptj/xptj(nptj),qptj(nptj),wptj(nptj)
1461       common/geom/rmproj,rmtarg,bmax,bkmx
1462       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat,seedp
1463       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctain/mtain
1464       double precision rcproj,rctarg
1465       common/geom1/rcproj,rctarg
1466
1467
1468       external sptj
1469
1470       call utpri('ainit ',ish,ishini,4)
1471
1472       inicnt=inicnt+1
1473       ntevt=0
1474
1475       if(inicnt.eq.1)write(ifmt,'(a)')'initializations ...'
1476
1477       if(seedi.ne.0d0)call ranfst(seedi)
1478       seedc=seedi
1479       call aseedi
1480
1481       if(isphis.eq.1)iframe=11  !nncm
1482       if(icinpu.ge.1)elepti=engy
1483 ctp060829      if(iopenu.eq.2)call smassi(themas)
1484       if(iopenu.eq.2.and.ish.eq.19)stop'change this?????????' !call smassp
1485
1486       if(iappl.eq.5)then
1487       yhaha=0
1488       ypjtl=0
1489       endif
1490
1491       if(ispherio.ne.0)ndecay=1
1492       if(ispherio.ne.0)idecay=0
1493       if(ispherio.ne.0)jdecay=0
1494       if(ifrade.eq.0)irescl=0
1495       idtarg=idtargin
1496       idproj=idprojin
1497       do 111 iii=1,4
1498  111  rexdif(iii)=abs(rexdifi(iii))
1499       if(noebin.gt.1)then
1500         engy=-1
1501         ekin=-1
1502         if(iologe.eq.1)engy=
1503      *       engmin*(engmax/engmin)**((real(nrebin)-0.5)/noebin)
1504         if(iologe.eq.0.or.(iologe.lt.0.and.iologl.lt.0))engy=
1505      *       engmin+(engmax-engmin)/noebin*(nrebin-0.5)
1506         if(iologl.eq.1)ekin=
1507      *       engmin*(engmax/engmin)**((real(nrebin-0.5))/noebin)
1508         if(iologl.eq.0)ekin=
1509      *       engmin+(engmax-engmin)/noebin*(real(nrebin)-0.5)
1510         elab=-1
1511         ecms=-1
1512         pnll=-1
1513         if(jpsi.lt.0)then
1514   11      z=0.19*sqrt(-2*alog(rangen()))*cos(2*pi*rangen())
1515           engy=abs(z)*engmax
1516           if(engy.lt.egymin)goto11
1517         endif
1518       endif
1519
1520            if(iappl.le.3)then
1521
1522         if(idtarg.eq.0)then   !in case of Air target, initialize with Argon nucleus
1523           if(model.eq.6)then     !no Argon in Sibyll
1524             latarg=7
1525             matarg=14
1526           else
1527             latarg=20
1528             matarg=40
1529           endif
1530         endif
1531
1532         if((idproj.ne.1120.and.(laproj.ne.-1.or.maproj.ne.1))
1533      &    .or.maproj.le.0)
1534      &  call utstop('Invalid projectile setup !&')
1535 c        if((idtarg.ne.1120.and.(latarg.ne.-1.or.matarg.ne.1))
1536 c     &    .or.matarg.le.0)
1537 c     &  call utstop('Invalid target setup !&')
1538
1539       if(iabs(idtarg).ne.1120.and.iabs(idtarg).ne.1220.and.idtarg.ne.0)
1540      &  call utstop('Invalid target !&')
1541         if((((idtarg.eq.-1120.or.iabs(idtarg).eq.1220)
1542      &    .and.(latarg.ne.-1.or.matarg.ne.1))
1543      &    .and.(idtarg.ne.1120.or.latarg.lt.0))
1544      &    .or.matarg.le.0)
1545      &  call utstop('Invalid target setup !&')
1546
1547
1548       call idmass(idproj,amproj)
1549       call idmass(idtarg,amtarg)
1550       nre=0
1551       if(engy.ge.0.)nre=nre+1
1552       if(pnll.ge.0.)nre=nre+1
1553       if(elab.ge.0.)nre=nre+1
1554       if(ekin.ge.0.)nre=nre+1
1555       if(ecms.ge.0.)nre=nre+1
1556       if(nre.ne.1)stop'invalid energy definition'
1557       ifirstghe=0
1558  101  continue
1559          if(engy.gt.0.)then
1560       pnll=sqrt(max(0., ((engy**2-amproj**2-amtarg**2)/2/amtarg)**2
1561      &                   -amproj**2) )
1562       elab=sqrt(pnll**2+amproj**2)
1563       ekin=elab-amproj
1564       ecms=engy
1565          elseif(ecms.gt.0.)then
1566       engy=ecms
1567       pnll=sqrt(max(0., ((engy**2-amproj**2-amtarg**2)/2/amtarg)**2
1568      &                   -amproj**2) )
1569       elab=sqrt(pnll**2+amproj**2)
1570       ekin=elab-amproj
1571          elseif(elab.gt.0)then
1572       pnll=sqrt(max(0.,elab**2-amproj**2))
1573       engy=sqrt( 2*elab*amtarg+amtarg**2+amproj**2 )
1574       ecms=engy
1575       ekin=elab-amproj
1576          elseif(pnll.gt.0)then
1577       elab=sqrt(pnll**2+amproj**2)
1578       engy=sqrt( 2*sqrt(pnll**2+amproj**2)*amtarg+amtarg**2+amproj**2 )
1579       ecms=engy
1580       ekin=elab-amproj
1581          elseif(ekin.gt.0.)then
1582       elab=ekin+amproj
1583       pnll=sqrt(max(0.,elab**2-amproj**2))
1584       engy=sqrt( 2*elab*amtarg+amtarg**2+amproj**2 )
1585       ecms=engy
1586          endif
1587
1588          if(model.eq.3.and.ifirstghe.eq.0)then    !det, trit and alp
1589            if(maproj.eq.2.and.laproj.eq.1)idproj=17
1590            if(maproj.eq.3.and.laproj.eq.1)idproj=18
1591            if(maproj.eq.4.and.laproj.eq.2)idproj=19
1592            if(idproj.ge.17.and.idproj.le.19)then
1593              elab=elab*maproj
1594              call idmass(idproj,amproj)
1595              maproj=1
1596              laproj=-1
1597              ifirstghe=1
1598              engy=-1
1599              ecms=-1
1600              pnll=-1
1601              ekin=-1
1602              goto 101
1603            endif
1604          endif
1605
1606       if(pnll.le.0.001)call utstop('ainit: energy too low&')
1607       if(engy.gt.egymax)call utstop('ainit: energy too high&')
1608       s=engy**2
1609       pnullx=utpcm(engy,amproj,amtarg)
1610       yhaha=alog((sqrt(pnll**2+s)+pnll)/sqrt(s))
1611       ypjtl=alog((sqrt(pnll**2+amproj**2)+pnll)/amproj)
1612
1613          elseif(iappl.eq.7)then
1614
1615       call idmass(idproj,amproj)
1616          if(elab.gt.0)then
1617       pnll=sqrt(max(0.,elab**2-amproj**2))
1618       engy=amproj
1619       ecms=engy
1620       ekin=elab-amproj
1621          elseif(pnll.gt.0)then
1622       elab=sqrt(pnll**2+amproj**2)
1623       engy=amproj
1624       ecms=engy
1625       ekin=elab-amproj
1626          elseif(ekin.gt.0.)then
1627       elab=ekin+amproj
1628       pnll=sqrt(max(0.,elab**2-amproj**2))
1629       engy=amproj
1630       ecms=engy
1631          else
1632       engy=amproj
1633       ecms=amproj
1634       elab=0.
1635       pnll=0.
1636       ekin=0.
1637         endif
1638
1639       pnullx=0.
1640       ypjtl=alog((sqrt(pnll**2+amproj**2)+pnll)/amproj)
1641       yhaha=ypjtl
1642
1643          endif
1644
1645       detap=(ypjtl-yhaha)*etafac
1646       detat=-yhaha*etafac
1647       tpro=dcosh(detap)
1648       zpro=dsinh(detap)
1649       ttar=dcosh(detat)
1650       ztar=dsinh(detat)
1651
1652       egyevt=engy
1653       pmxevt=pnll
1654
1655       if(iappl.gt.8)stop'update following statement'
1656       if(iappl.ge.5.and.iappl.le.8)then
1657       s=12.**2
1658       endif
1659
1660       if(iappl.le.3)then
1661        if(maproj.gt.1)then
1662         rpj=1.19*maproj**(1./3.)-1.61*maproj**(-1./3.)
1663         rmproj=rpj+fctrmx*.54
1664         rcproj=dble(rpj/cosh(yhaha)*facnuc)
1665        else
1666         rmproj=0
1667         rcproj=dble(0.8/cosh(yhaha)*facnuc)
1668        endif
1669        if(matarg.gt.1)then
1670         rtg=1.19*matarg**(1./3.)-1.61*matarg**(-1./3.)
1671         rmtarg=rtg+fctrmx*.54
1672         rctarg=dble(rtg/cosh(yhaha)*facnuc)
1673        else
1674         rmtarg=0
1675         rctarg=dble(0.8/cosh(yhaha)*facnuc)
1676        endif
1677
1678       endif
1679
1680       call iclass(idproj,iclpro)
1681       call iclass(idtarg,icltar)
1682
1683          if(inicnt.eq.1)then
1684
1685           call ranfgt(seedp)                     !not to change the seed ...
1686
1687       call hdecin(.false.)
1688
1689       if(iappl.eq.1.or.iappl.ge.5)then
1690       c=6
1691       call utquaf(sptj,nptj,xptj,qptj,wptj,0.,.33*c,.66*c,c)
1692       endif
1693
1694       if(iappl.ne.2)then
1695         call hnbspd(iospec)
1696         ktnbod=0
1697 c        if(model.eq.1)call hnbxxxini
1698         call hnbpajini
1699       endif
1700
1701       if(model.eq.1)then
1702         if(iclegy2.gt.1)then
1703           egyfac=(egymax*1.0001/egylow)**(1./float(iclegy2-1))
1704         else
1705           egyfac=1.
1706         endif
1707       endif
1708
1709       if(model.eq.1)then
1710         call conini
1711       else
1712         iorsce=0
1713         iorsdf=0
1714         iorshh=0
1715         iorsdf=0
1716         irescl=0       !rescaling
1717       endif
1718
1719       if(iappl.ne.6.and.model.eq.1)call psaini
1720
1721           call ranfst(seedp)                     ! ... after this initialization
1722
1723          endif
1724
1725 c$$$      if(idproj.eq.1120)icp=2        !????????????? for what ?
1726 c$$$      if(idproj.eq.-1120)icp=-2
1727 c$$$      if(idproj.eq.120)icp=1
1728 c$$$      if(idproj.eq.-120)icp=-1
1729 c$$$      if(idproj.eq.130)icp=4
1730 c$$$      if(idproj.eq.-130)icp=-4
1731
1732
1733
1734       if(model.eq.1)then                   !only for epos
1735
1736
1737       if(iappl.le.3)then
1738         call emsini(engy,idproj,idtarg)
1739         call paramini(1)
1740         if(ish.ge.4)then
1741           do i=idxD0,idxD1
1742             write(ifch,'(9(a,f8.4))')
1743      *      'AlpD:',alpD(i,iclpro,icltar)
1744      * ,'    AlpDp:',alpDp(i,iclpro,icltar)
1745      * ,'    AlpDpp:',alpDpp(i,iclpro,icltar)
1746      * ,'    BetD:',betD(i,iclpro,icltar)
1747      * ,'    BetDp:',betDp(i,iclpro,icltar)
1748      * ,'    BetDpp:',betDpp(i,iclpro,icltar)
1749      * ,'    GamD:',gamD(i,iclpro,icltar)
1750      * ,'    DelD:',delD(i,iclpro,icltar)
1751      * ,'    AlpPar:',alppar
1752      * ,'    bkmxdif:',bmxdif(iclpro,icltar)
1753           enddo
1754         endif
1755       endif
1756
1757       if(iappl.le.3)then
1758         bkmxndif=conbmxndif()
1759         bkmx=conbmx()
1760         if(ish.ge.3)write(ifch,*)'bkmx,bkmxndif',bkmx,bkmxndif
1761
1762         if(maproj.gt.1.or.matarg.gt.1)then
1763         bmax=rmproj+rmtarg
1764         else
1765         bmax=bkmx
1766         endif
1767       endif
1768
1769       if(ixtau.eq.1)call xtauev(0)
1770
1771       if(model.eq.1)then
1772       if(iEmsB.eq.1)call xEmsB(0,0,0)
1773       if(iEmsBg.eq.1)call xEmsBg(0,0,0)
1774       if(iEmsPm.eq.1)call xEmsPm(0,0,0)
1775       if(iEmsPx.eq.1)call xEmsPx(0,0.,0.,0)
1776       if(iEmsPBx.eq.1)call xEmsP2(0,0,0,0.,0.,0.,0.,0.,0.)
1777 c      if(iEmsPx.eq.1)call xEmsPxNo(0,0.,0.,0,0)
1778       if(iEmsSe.eq.1)call xEmsSe(0,0.,0.,0,1)
1779       if(iEmsSe.eq.1)call xEmsSe(0,0.,0.,0,2)
1780       if(iEmsDr.eq.1)call xEmsDr(0,0.,0.,0)
1781       if(iEmsRx.eq.1)call xEmsRx(0,0,0.,0.)
1782       endif
1783
1784 c G function parameters
1785
1786       if(iappl.eq.1)then
1787         call GfunPar(1,1,0.,s,alp,bet,betp,epsp,epst,epss,gamvv)  !b=0
1788         epszero=epss
1789         do i=1,nclha
1790          alpff(i)=   engy**epszero*gamhad(i)
1791         enddo
1792         betff(1)=   -alppar+epsp
1793         betff(2)=   -alppar+epst
1794       endif
1795
1796       endif
1797
1798 c additional initialization procedures
1799
1800       if(model.ne.1)then
1801         call IniEvtModel
1802       elseif(iappl.le.3)then
1803 c Cross section calculation
1804         call xsigma
1805       endif
1806
1807
1808       if(idtarg.eq.0)idtarg=1120 !air = nucleus
1809
1810
1811
1812       if(seedj.ne.0.)call ranfst(seedj)
1813       seedc=seedj
1814
1815
1816 ccc      call MakeFpartonTable
1817
1818 c$$$      call testconex(1)
1819
1820       call utprix('ainit ',ish,ishini,4)
1821       return
1822       end
1823
1824 c---------------------------------------------------------------------
1825       subroutine aread
1826 c---------------------------------------------------------------------
1827 c  reads and interprets input commands
1828 c---------------------------------------------------------------------
1829
1830       include 'epos.inc'
1831       include 'epos.incpar'
1832       include 'epos.incsem'
1833
1834       double precision histoweight
1835       common/chiswei/histoweight
1836       common/cyield/yield/cifset/ifset/caverg/averg
1837       common/csigma/sigma
1838       double precision val,val1,val2
1839       character*160 line,linex,cline
1840       data nappl /0/
1841       common/record/maxrec(2),irecty(30,2)
1842       common/cfacmss/facmss /cr3pomi/r3pomi,r4pomi
1843       common /ems12/iodiba,bidiba  ! defaut iodiba=0. if iodiba=1, study H-Dibaryon
1844       character*80 fndat,fnncs,fnIIdat,fnIIncs                 !qgs-II????????
1845       common/qgsfname/  fndat, fnncs, ifdat, ifncs
1846       common/qgsIIfname/fnIIdat, fnIIncs, ifIIdat, ifIIncs     !qgs-II????????
1847       common/qgsnfname/ nfndat, nfnncs
1848       common/qgsIInfname/ nfnIIdat, nfnIIncs     !qgs-II????????
1849       common/ghecsquel/anquasiel,iquasiel
1850       common/cjjj/jjj,cline
1851       character cmodel*21
1852       parameter(mxdefine=40)
1853       character w1define*100,w2define*100
1854       common/cdefine/ndefine,l1define(mxdefine),l2define(mxdefine)
1855      &               ,w1define(mxdefine),w2define(mxdefine)
1856       common/cbincond/nozero,ibmin,ibmax
1857
1858           j=-1
1859       nhsto=0
1860       ndefine=0
1861
1862     1 call utword(line,i,j,1)
1863
1864           if(line(i:j).eq.'#define')then
1865
1866       call utword(line,i,j,ne)
1867       if(line(i:j).eq.'bim3')stop'****** bim3->bim03 ****** '
1868       if(line(i:j).eq.'bim5')stop'****** bim5->bim05 ****** '
1869       if(line(i:j).eq.'bim6')stop'****** bim6->bim06 ****** '
1870       if(ndefine+1.gt.mxdefine)stop'too many `define` statements.      '
1871       l1=j-i+1
1872       if(l1.gt.99)stop'`define` argument 1 too long.            '
1873       w1define(ndefine+1)(1:l1)=line(i:j)
1874       w1define(ndefine+1)(l1+1:l1+1)=' '
1875       call utword(line,i,j,ne)
1876       l2=j-i+1
1877       if(l2.gt.99)stop'`define` argument 2 too long.            '
1878       w2define(ndefine+1)(1:l2)=line(i:j)
1879       w2define(ndefine+1)(l2+1:l2+1)=' '
1880       ndefine=ndefine+1
1881       l1define(ndefine)=l1
1882       l2define(ndefine)=l2
1883
1884           elseif(line(i:j).eq.'goto')then
1885
1886       call utword(line,i,j,ne)
1887       ix=i
1888       jx=j
1889       linex=line
1890       call utword(line,i,j,ne)
1891       do while(line(i:j).ne.linex(ix:jx))
1892       call utword(line,i,j,ne)
1893       enddo
1894       goto1
1895
1896            elseif(line(i:j).eq.'application')then
1897
1898       call utworn(line,j,ne)
1899       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'application?'
1900       call utword(line,i,j,0)
1901       if(nopen.ne.-1)then       !only first read
1902       if(line(i:j).eq.'analysis')  iappl=0
1903       if(line(i:j).eq.'hadron')    iappl=1
1904       if(line(i:j).eq.'geometry')  iappl=2
1905       if(line(i:j).eq.'read')      iappl=3
1906       if(line(i:j).eq.'micro')     iappl=4
1907       if(line(i:j).eq.'kinky')     iappl=5
1908       if(line(i:j).eq.'ee')        iappl=6
1909       if(line(i:j).eq.'decay')     iappl=7
1910       if(line(i:j).eq.'lepton')    iappl=8
1911       if(line(i:j).eq.'ee')    then
1912         naflav=5                ! number of flavors
1913       endif
1914       nappl=nappl+1
1915       if(iappl.gt.0.and.nappl.gt.1)call aaset(1)
1916       if(iappl.eq.0.and.nappl.gt.1)call aaset(2)
1917       if(iappl.eq.0)jframe=iframe
1918       if(iappl.eq.0)kframe=iframe
1919       if(iappl.eq.1)iframe=0
1920       if(iappl.eq.2)iframe=0
1921       if(iappl.eq.4)iframe=1
1922       if(iappl.eq.5)iframe=1
1923       if(iappl.eq.6)iframe=1
1924       if(iappl.eq.7)iframe=0
1925       if(iappl.eq.8)iframe=21        !gncm
1926       endif
1927
1928            elseif(line(i:j).eq.'call')then
1929
1930       call utworn(line,j,ne)
1931       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'subroutine?'
1932       call utword(line,i,j,0)
1933
1934       if(nopen.eq.-1)then       !-----only second run
1935
1936       if(line(i:j).eq.'xConThickProj')call xConThick(1)
1937       if(line(i:j).eq.'xConThickTarg')call xConThick(2)
1938       if(line(i:j).eq.'xConNuclDensProj')call xConNuclDens(1)
1939       if(line(i:j).eq.'xConNuclDensTarg')call xConNuclDens(2)
1940       if(line(i:j).eq.'xConNuclDensProjTarg')call xConNuclDens(1)
1941       if(line(i:j).eq.'xConNuclDensProjTarg')call xConNuclDens(2)
1942       if(line(i:j).eq.'xFom')call xfom
1943       if(line(i:j).eq.'xGeometry')call xGeometry(2)
1944       if(line(i:j).eq.'xEpsilon')call xEpsilon(2)
1945       if(line(i:j).eq.'xRanPt')call xRanPt
1946       if(line(i:j).eq.'xParGam')call xParGam
1947       if(line(i:j).eq.'xParGampp')call xParGampp
1948       if(line(i:j).eq.'xParOmega1xy')call xParOmega1xy
1949 c$$$      if(line(i:j).eq.'xParOmega3xyz')call xParOmega3xyz
1950       if(line(i:j).eq.'xParPro')call xParPro
1951       if(line(i:j).eq.'xParPro1')call xParPro1
1952       if(line(i:j).eq.'xParPomInc')call xParPomInc
1953       if(line(i:j).eq.'xParPomIncX')call xParPomIncX
1954       if(line(i:j).eq.'xParPomIncP')call xParPomIncP
1955       if(line(i:j).eq.'xParPomIncM')call xParPomIncM
1956       if(line(i:j).eq.'xParPomIncXI')call xParPomIncXI
1957       if(line(i:j).eq.'xParPomIncPI')call xParPomIncPI
1958       if(line(i:j).eq.'xParPomIncMI')call xParPomIncMI
1959       if(line(i:j).eq.'xParPomIncJ')call xParPomIncJ
1960       if(line(i:j).eq.'xParOmega1')call xParOmega1
1961 c$$$      if(line(i:j).eq.'xParOmega3')call xParOmega3
1962 c$$$      if(line(i:j).eq.'xParOmega5')call xParOmega5
1963       if(line(i:j).eq.'xParOmegaN')call xParOmegaN
1964       if(line(i:j).eq.'xParGauss')call xParGauss
1965       if(line(i:j).eq.'xParSigma')call xParSigma
1966 c$$$      if(line(i:j).eq.'xParSigma2')call xParSigma2
1967 c$$$      if(line(i:j).eq.'xScrD')call xScrD
1968       if(line(i:j).eq.'xFitD1')call xFitD1
1969 c$$$      if(line(i:j).eq.'xExaD2')call xExaD2
1970       if(line(i:j).eq.'xbExaD')call xbExaD
1971 c$$$      if(line(i:j).eq.'xbExaD2')call xbExaD2
1972       if(line(i:j).eq.'xbnExaD')call xbnExaD
1973 c$$$      if(line(i:j).eq.'xbnExaD2')call xbnExaD2
1974       if(line(i:j).eq.'xFitD2')call xFitD2
1975       if(line(i:j).eq.'xbParD')call xbParD
1976 c$$$      if(line(i:j).eq.'xParD2')call xParD2
1977       if(line(i:j).eq.'xGexaJ')call xGexaJ
1978       if(line(i:j).eq.'xbnParD')call xbnParD
1979       if(line(i:j).eq.'xsParD')call xsParD
1980 c$$$      if(line(i:j).eq.'xmParD2')call xmParD2
1981       if(line(i:j).eq.'xyParD')call xyParD
1982 c$$$      if(line(i:j).eq.'xyParD2')call xyParD2
1983       if(line(i:j).eq.'xParPhi1')call xParPhi1
1984       if(line(i:j).eq.'xParPhi')call xParPhi
1985       if(line(i:j).eq.'xParH')call xParH
1986       if(line(i:j).eq.'xParHPhiInt')call xParHPhiInt
1987       if(line(i:j).eq.'xParZ')call xParZ
1988       if(line(i:j).eq.'xtauev')call xtauev(2)
1989       if(line(i:j).eq.'xspace')call xspace(2)
1990       if(line(i:j).eq.'gakjto'   )call gakjto
1991       if(line(i:j).eq.'psaevp')call psaevp
1992 c     if(line(i:j).eq.'pyarea')call pyarea
1993       if(line(i:j).eq.'xhgcam')call xhgcam(0.,-1)
1994       if(line(i:j).eq.'xhgcfl')call xhgcfl(0.,0.,0.,-1)
1995       if(line(i:j).eq.'xhgccc')call xhgccc(-1.)
1996       if(line(i:j).eq.'xjden1')call xjden1(2,0,0.,0.,0.,0.,0.)
1997       if(line(i:j).eq.'xjden2')call xjden2(2,0,0.,0.,0.,0.)
1998 c     if(line(i:j).eq.'xjdis' )call xjdis(2,0,0)
1999       if(line(i:j).eq.'xhgcen')call xhgcen
2000       if(line(i:j).eq.'xhgcmt')call xhgcmt
2001       if(line(i:j).eq.'xhgcmu')call xhgcmu
2002       if(line(i:j).eq.'xhgcmx')call xhgcmx
2003       if(line(i:j).eq.'xhgcpt')call xhgcpt
2004       if(line(i:j).eq.'xhgcra')call xhgcra
2005       if(line(i:j).eq.'xhnben')call xhnben
2006       if(line(i:j).eq.'xhnbit')call xhnbit
2007       if(line(i:j).eq.'xhnbmu')call xhnbmu
2008       if(line(i:j).eq.'xhnbmz')call xhnbmz
2009       if(line(i:j).eq.'xhnbte')call xhnbte(-1)
2010       if(line(i:j).eq.'xhnbti')call xhnbti(-1)
2011       if(model.eq.1)then
2012       if(line(i:j).eq.'xEmsB' )call xEmsB(2,0,0)
2013       if(line(i:j).eq.'xEmsBg')call xEmsBg(2,0,0)
2014       if(line(i:j).eq.'xEmsPm')call xEmsPm(2,0,0)
2015       if(line(i:j).eq.'xEmsPx')call xEmsPx(2,0.,0.,0)
2016 c      if(line(i:j).eq.'xEmsPx')call xEmsPxNo(2,0.,0.,0,0)
2017       if(line(i:j).eq.'xEmsSe')call xEmsSe(2,0.,0.,0,1)
2018       if(line(i:j).eq.'xEmsSe')call xEmsSe(2,0.,0.,0,2)
2019       if(line(i:j).eq.'xEmsDr')call xEmsDr(2,0.,0.,0)
2020       if(line(i:j).eq.'xEmsRx')call xEmsRx(2,0,0.,0.)
2021       if(line(i:j-4).eq.'xEmsP2')then
2022         read(line(j-1:j-1),*)val
2023         idh=nint(val)
2024         read(line(j:j),*)val
2025         jex=nint(val)
2026         if(line(j-3:j-2).eq.'PE')
2027      &       call xEmsP2(2,idh,jex,0.,0.,0.,0.,0.,0.)
2028         if(line(j-3:j-2).eq.'IB')
2029      &       call xEmsP2(3,idh,jex,0.,0.,0.,0.,0.,0.)
2030         if(line(j-3:j-2).eq.'OB')
2031      &       call xEmsP2(4,idh,jex,0.,0.,0.,0.,0.,0.)
2032       endif
2033       if(line(i:j).eq.'xConxyzProj')
2034      &stop'xConxyzProj->xConNuclDensProj'
2035       if(line(i:j).eq.'xConxyzTarg')
2036      &stop'xConxyzTarg->xConNuclDensTarg'
2037       if(line(i:j).eq.'xConxyzProjTarg')
2038      &stop'xConxyzProjTarg->xConNuclDensProjTarg'
2039
2040       endif
2041
2042       elseif(model.eq.1)then  !first run and epos
2043
2044       if(line(i:j).eq.'xGeometry')then
2045        call xGeometry(0)
2046        ixgeometry=1
2047       elseif(line(i:j).eq.'xEpsilon')then
2048        call xEpsilon(0)
2049       elseif(line(i:j).eq.'xtauev')then
2050        ixtau=1
2051       elseif(line(i:j).eq.'xEmsB')then
2052        iEmsB=1
2053       elseif(line(i:j).eq.'xEmsBg')then
2054        iEmsBg=1
2055       elseif(line(i:j).eq.'xEmsPm')then
2056        iEmsPm=1
2057       elseif(line(i:j).eq.'xEmsPx')then
2058        iEmsPx=1
2059       elseif(line(i:j-4).eq.'xEmsP2')then
2060        iEmsPBx=1
2061       elseif(line(i:j).eq.'xEmsSe')then
2062        iEmsSe=1
2063       elseif(line(i:j).eq.'xEmsDr')then
2064        iEmsDr=1
2065       elseif(line(i:j).eq.'xEmsRx')then
2066        iEmsRx=1
2067       elseif(line(i:j).eq.'xEmsI1')then
2068        iEmsI1=1
2069        if(iEmsI1+iEmsI2.eq.1)write(ifhi,'(a)')'newpage zone 3 4 1'
2070       elseif(line(i:j).eq.'xEmsI2')then
2071        iEmsI2=1
2072        if(iEmsI1+iEmsI2.eq.1)write(ifhi,'(a)')'newpage zone 3 4 1'
2073       elseif(line(i:j).eq.'xSpaceTime')then
2074        iSpaceTime=1
2075       elseif(line(i:j).eq.'xxSpaceTime')then
2076        stop'xxSpaceTime->xSpaceTime.              '
2077       endif
2078       endif
2079
2080            elseif(line(i:j).eq.'decayall')then
2081
2082       nrnody=0
2083
2084            elseif(line(i:j).eq.'echo')then
2085
2086       call utworn(line,j,ne)
2087       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'on or off?'
2088       call utword(line,i,j,0)
2089       if(line(i:j).eq.'on')iecho=1
2090       if(line(i:j).eq.'off')iecho=0
2091       if(line(i:j).ne.'on'.and.line(i:j).ne.'off')stop'invalid option'
2092       if(nopen.eq.-1)iecho=0
2093
2094            elseif(line(i:j).eq.'fqgsjet')then              !QGSJet
2095
2096       call utworn(line,j,ne)
2097       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-type file-name?'
2098       call utword(line,i,j,0)
2099       linex=line
2100       ix=i
2101       jx=j
2102       call utworn(line,j,ne)
2103       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-name?'
2104       call utword(line,i,j,0)
2105       if(linex(ix:jx).eq.'dat')fndat(1:j-i+1)=line(i:j)
2106       if(linex(ix:jx).eq.'dat')nfndat=j-i+1             !length of qgsdat01-file name
2107       if(linex(ix:jx).eq.'ncs')fnncs(1:j-i+1)=line(i:j)
2108       if(linex(ix:jx).eq.'ncs')nfnncs=j-i+1             !length of sectnu-file name
2109       if(nfndat.gt.1)ifdat=1
2110       if(nfnncs.gt.1)ifncs=2
2111
2112            elseif(line(i:j).eq.'fqgsjetII')then              !QGSJET-II     !qgs-II????????
2113
2114       call utworn(line,j,ne)
2115       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-type file-name?'
2116       call utword(line,i,j,0)
2117       linex=line
2118       ix=i
2119       jx=j
2120       call utworn(line,j,ne)
2121       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-name?'
2122       call utword(line,i,j,0)
2123       if(linex(ix:jx).eq.'dat')fnIIdat(1:j-i+1)=line(i:j)
2124       if(linex(ix:jx).eq.'dat')nfnIIdat=j-i+1             !length of qgsjet-II.dat name
2125       if(linex(ix:jx).eq.'ncs')fnIIncs(1:j-i+1)=line(i:j)
2126       if(linex(ix:jx).eq.'ncs')nfnIIncs=j-i+1             !length of qgsjet-II.ncs name
2127       if(nfnIIdat.gt.1)ifIIdat=1
2128       if(nfnIIncs.gt.1)ifIIncs=2
2129
2130            elseif(line(i:j).eq.'fname')then
2131
2132       call utworn(line,j,ne)
2133       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-type file-name?'
2134       call utword(line,i,j,0)
2135       linex=line
2136       ix=i
2137       jx=j
2138       call utworn(line,j,ne)
2139       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'file-name?'
2140       call utword(line,i,j,0)
2141       if(linex(ix:jx).eq.'pathnx')fnnx(1:j-i+1)=line(i:j)
2142       if(linex(ix:jx).eq.'check')fnch(1:j-i+1)=line(i:j)
2143       if(linex(ix:jx).eq.'histo')fnhi(1:j-i+1)=line(i:j)
2144       if(linex(ix:jx).eq.'data') fndt(1:j-i+1)=line(i:j)
2145       if(linex(ix:jx).eq.'copy') fncp(1:j-i+1)=line(i:j)
2146       if(linex(ix:jx).eq.'initl') fnii(1:j-i+1)=line(i:j)
2147       if(linex(ix:jx).eq.'initl+') fnii(nfnii+1:nfnii+j-i+1)=line(i:j)
2148       if(linex(ix:jx).eq.'inidi') fnid(1:j-i+1)=line(i:j)
2149       if(linex(ix:jx).eq.'inidi+') fnid(nfnid+1:nfnid+j-i+1)=line(i:j)
2150       if(linex(ix:jx).eq.'inidr') fndr(1:j-i+1)=line(i:j)
2151       if(linex(ix:jx).eq.'inidr+') fndr(nfndr+1:nfndr+j-i+1)=line(i:j)
2152       if(linex(ix:jx).eq.'iniev') fnie(1:j-i+1)=line(i:j)
2153       if(linex(ix:jx).eq.'iniev+') fnie(nfnie+1:nfnie+j-i+1)=line(i:j)
2154       if(linex(ix:jx).eq.'inirj') fnrj(1:j-i+1)=line(i:j)
2155       if(linex(ix:jx).eq.'inirj+') fnrj(nfnrj+1:nfnrj+j-i+1)=line(i:j)
2156       if(linex(ix:jx).eq.'inics') fncs(1:j-i+1)=line(i:j)
2157       if(linex(ix:jx).eq.'inics+') fncs(nfncs+1:nfncs+j-i+1)=line(i:j)
2158       if(linex(ix:jx).eq.'inihy') fnhy(1:j-i+1)=line(i:j)
2159       if(linex(ix:jx).eq.'pathnx')nfnnx=j-i+1
2160       if(linex(ix:jx).eq.'check')nfnch=j-i+1
2161       if(linex(ix:jx).eq.'histo')nfnhi=j-i+1
2162       if(linex(ix:jx).eq.'data') nfndt=j-i+1
2163       if(linex(ix:jx).eq.'copy') nfncp=j-i+1
2164       if(linex(ix:jx).eq.'initl') nfnii=j-i+1
2165       if(linex(ix:jx).eq.'initl+')nfnii=nfnii+j-i+1
2166       if(linex(ix:jx).eq.'inidi') nfnid=j-i+1
2167       if(linex(ix:jx).eq.'inidi+')nfnid=nfnid+j-i+1
2168       if(linex(ix:jx).eq.'inidr') nfndr=j-i+1
2169       if(linex(ix:jx).eq.'inidr+')nfndr=nfndr+j-i+1
2170       if(linex(ix:jx).eq.'iniev') nfnie=j-i+1
2171       if(linex(ix:jx).eq.'iniev+')nfnie=nfnie+j-i+1
2172       if(linex(ix:jx).eq.'inirj') nfnrj=j-i+1
2173       if(linex(ix:jx).eq.'inirj+')nfnrj=nfnrj+j-i+1
2174       if(linex(ix:jx).eq.'inics') nfncs=j-i+1
2175       if(linex(ix:jx).eq.'inihy') nfnhy=j-i+1
2176       if(linex(ix:jx).eq.'inics+')nfncs=nfncs+j-i+1
2177       if(linex(ix:jx).eq.'check'.and.fnch(1:nfnch).ne.'none') then
2178       open(unit=ifcx,file=fnch(1:nfnch),status='unknown')
2179       kchopen=1
2180       elseif(linex(ix:jx).eq.'pathnx'.and.fnnx(1:nfnnx).ne.'none')then
2181       knxopen=1
2182       elseif(linex(ix:jx).eq.'histo'.and.fnhi(1:nfnhi).ne.'none')then
2183       open(unit=ifhi,file=fnhi(1:nfnhi),status='unknown')
2184       khiopen=1
2185       elseif(linex(ix:jx).eq.'data'.and.fndt(1:nfndt).ne.'none')then
2186       open(unit=ifdt,file=fndt(1:nfndt),status='unknown')
2187       kdtopen=1
2188       write(6,'(a,a)')'Opened data file ',fndt(1:nfndt)
2189       elseif(linex(ix:jx).eq.'copy'.and.fncp(1:nfncp).ne.'none')then
2190       open(unit=ifcp,file=fncp(1:nfncp),status='unknown')
2191       kcpopen=1
2192       endif
2193
2194            elseif(line(i:j).eq.'frame')then
2195
2196       call utworn(line,j,ne)
2197       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'frame?'
2198       call utword(line,i,j,0)
2199         if(nopen.ne.-1)then       ! event definition only, not analysis
2200         if(line(i:j).eq.'nucleon-nucleon')then
2201       if(iappl.eq.0)jframe=11
2202       if(iappl.eq.1)iframe=11
2203       if(iappl.eq.3)iframe=11
2204       if(iappl.gt.3.or.iappl.eq.2)stop'invalid option nucleon-nucleon'
2205         elseif(line(i:j).eq.'target')then
2206       if(iappl.eq.0)jframe=12
2207       if(iappl.eq.1)iframe=12
2208       if(iappl.eq.3)iframe=12
2209       if(iappl.gt.3.or.iappl.eq.2)stop'invalid option target'
2210         elseif(line(i:j).eq.'gamma-nucleon')then
2211       if(iappl.eq.0)jframe=21
2212       if((iappl+1)/2.eq.4)iframe=21
2213       if(iappl.eq.3)iframe=21
2214       if(iappl.ne.0.and.(iappl+1)/2.ne.4.and.iappl.ne.3)
2215      *stop'invalid option gamma-nucleon'
2216         elseif(line(i:j).eq.'lab')then
2217       if(iappl.eq.0)jframe=22
2218       if((iappl+1)/2.eq.4)iframe=22
2219       if(iappl.eq.3)iframe=22
2220       if(iappl.ne.0.and.(iappl+1)/2.ne.4.and.iappl.ne.3)
2221      *stop'invalid option lab'
2222         elseif(line(i:j).eq.'sphericity')then
2223       if(iappl.eq.0)jframe=32
2224       if(iappl.ne.0)stop'invalid option sphericity'
2225         elseif(line(i:j).eq.'thrust')then
2226       if(iappl.eq.0)jframe=33
2227       if(iappl.ne.0)stop'invalid option thrust'
2228         elseif(line(i:j).eq.'breit')then
2229           if(iappl.ne.0)stop'invalid option breit'
2230         else
2231       stop'frame not recognized'
2232         endif
2233         endif
2234
2235            elseif(line(i:j).eq.'frame+')then
2236
2237       call utword(line,i,j,0)
2238
2239            elseif(line(i:j).eq.'binning')then
2240
2241       call utworn(line,j,ne)
2242       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'log/lin?'
2243       call utword(line,i,j,0)
2244       if(line(i:j).eq.'lin')iologb=0
2245       if(line(i:j).eq.'log')iologb=1
2246
2247            elseif(line(i:j).eq.'beginhisto')then
2248
2249       if(nopen.ne.-1)then
2250         jjj=j
2251         cline=line
2252         call xini
2253         j=jjj
2254         line=cline
2255       endif
2256
2257            elseif(line(i:j).eq.'endhisto')then
2258
2259       if(nopen.eq.-1)then
2260         nhsto=nhsto+1
2261         call xhis(nhsto)
2262       endif
2263
2264            elseif(line(i:j).eq.'noweak')then
2265
2266       !only used in xini
2267
2268            elseif(line(i:j).eq.'histogram')then
2269
2270       call utword(line,i,j,0)
2271       call utword(line,i,j,0)
2272       call utword(line,i,j,0)
2273       call utword(line,i,j,0)
2274       call utword(line,i,j,0)
2275       call utword(line,i,j,0)
2276
2277            elseif(line(i:j).eq.'plot')then
2278
2279       stop' "plot" not used any more'
2280
2281            elseif(line(i:j).eq.'root')then
2282
2283       stop' "root" not used any more'
2284
2285            elseif(line(i:j).eq.'idcode')then
2286
2287       call utword(line,i,j,0)
2288
2289            elseif(line(i:j).eq.'xpara')then
2290
2291       call utword(line,i,j,0)
2292       call utword(line,i,j,0)
2293
2294            elseif(line(i:j).eq.'histoweight')then
2295
2296       if(nopen.eq.-1)then
2297       write(ifhi,'(a,d22.14)')'histoweight ',histoweight
2298       endif
2299
2300            elseif(line(i:j).eq.'input')then
2301
2302       call utworn(line,j,ne)
2303       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'input file?'
2304       call utword(line,i,j,0)
2305       if(nopen.ge.0)then
2306        nopen=nopen+1
2307        if(nopen.gt.9)stop'too many nested input commands'
2308        open(unit=20+nopen,file=line(i:j),status='old')
2309        if(iprmpt.eq.1)iprmpt=-1
2310       endif
2311
2312            elseif(line(i:j).eq.'nodecays')then
2313
2314       call utword(line,i,j,0)
2315       do while (line(i:j).ne.'end')
2316        if(nrnody.ge.mxnody)then
2317         write(ifmt,'(a)')'too many nodecays; command ignored'
2318        else
2319         nrnody=nrnody+1
2320         read(line(i:j),*)val
2321         nody(nrnody)=nint(val)
2322        endif
2323       call utword(line,i,j,0)
2324       enddo
2325
2326            elseif(line(i:j).eq.'nodecay')then
2327
2328       call utword(line,i,j,0)
2329       if(nopen.ne.-1)then
2330       if(nrnody.ge.mxnody)then
2331       write(ifmt,'(a)')'too many nodecay commands; command ignored'
2332       j=160
2333       i=j+1
2334       goto1
2335       endif
2336       nrnody=nrnody+1
2337       read(line(i:j),*)val
2338       nody(nrnody)=nint(val)
2339       endif
2340
2341            elseif(line(i:j).eq.'dodecay')then
2342
2343       call utword(line,i,j,0)
2344       if(nopen.ne.-1)then
2345       read(line(i:j),*)val
2346       idx=nint(val)
2347       nrn=0
2348       imv=0
2349       do while(nrn.lt.nrnody)
2350        nrn=nrn+1
2351        if(idx.eq.nody(nrn))then
2352          nrnody=nrnody-1
2353          imv=1
2354        endif
2355        if(imv.eq.1.and.nrn.le.nrnody)nody(nrn)=nody(nrn+1)
2356       enddo
2357       endif
2358
2359            elseif(line(i:j).eq.'print')then
2360
2361       call utworn(line,j,ne)
2362       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'subroutine?'
2363       call utword(line,i,j,0)
2364        if(line(i:j).eq.'*')then
2365       nrpri=0
2366       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'level?'
2367       call utword(line,i,j,0)
2368       read(line(i:j),*)val
2369       ish=nint(val)
2370        else
2371       nrpri=nrpri+1
2372       subpri(nrpri)='                    '
2373       subpri(nrpri)(1:j-i+1)=line(i:j)
2374       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'level?'
2375       call utword(line,i,j,0)
2376       read(line(i:j),*)val
2377       ishpri(nrpri)=nint(val)
2378        endif
2379
2380            elseif(line(i:j).eq.'printcheck')then
2381
2382       call utworn(line,j,ne)
2383       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'screen or file?'
2384       call utword(line,i,j,0)
2385       if(line(i:j).eq.'screen')ifch=ifmt
2386       if(line(i:j).eq.'file')ifch=ifcx
2387       if(line(i:j).ne.'screen'.and.line(i:j).ne.'file')
2388      *write(ifmt,'(a)')'invalid option; command ignored'
2389
2390            elseif(line(i:j).eq.'prompt')then
2391
2392       call utworn(line,j,ne)
2393       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'on or off or auto?'
2394       call utword(line,i,j,0)
2395       if(line(i:j).eq.'on')iprmpt=2
2396       if(line(i:j).eq.'off')iprmpt=-2
2397       if(line(i:j).eq.'auto'.and.nopen.eq.0)iprmpt=1
2398       if(line(i:j).eq.'auto'.and.nopen.eq.1)iprmpt=-1
2399       if(line(i:j).ne.'on'.and.line(i:j).ne.'off'
2400      *.and.line(i:j).ne.'auto')stop'invalid option'
2401
2402            elseif(line(i:j).eq.'return')then
2403
2404       if(nopen.ne.-1)then
2405         close(ifop)
2406         nopen=nopen-1
2407         if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2408 ccc      close(20+nopen)
2409 ccc      nopen=nopen-1
2410 ccc      if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2411       endif
2412
2413            elseif(line(i:j).eq.'run')then
2414
2415            elseif(line(i:j).eq.'runprogram')then
2416
2417       if(kchopen.eq.0.and.fnch(1:nfnch).ne.'none')then
2418         open(unit=ifcx,file=fnch(1:nfnch),status='unknown')
2419         kchopen=1
2420       endif
2421       if(khiopen.eq.0.and.fnhi(1:nfnhi).ne.'none')then
2422         open(unit=ifhi,file=fnhi(1:nfnhi),status='unknown')
2423         khiopen=1
2424       endif
2425       if(kdtopen.eq.0.and.fndt(1:nfndt).ne.'none')then
2426         open(unit=ifdt,file=fndt(1:nfndt),status='unknown')
2427         kdtopen=1
2428         write(6,'(a,a)')'Opened data file ',fndt(1:nfndt)
2429       endif
2430       if(kcpopen.eq.0.and.fncp(1:nfncp).ne.'none')then
2431         open(unit=ifcp,file=fncp(1:nfncp),status='unknown')
2432         kcpopen=1
2433       endif
2434       return
2435
2436            elseif(line(i:j).eq.'if')then
2437
2438       call utword(line,i,j,0)
2439       ix=i
2440       jx=j
2441       linex=line
2442       call utword(line,i,j,0)
2443       read(line(i:j),*)val1
2444       call utword(line,i,j,0)
2445       read(line(i:j),*)val2
2446       ifset=1
2447       if(linex(ix:jx).eq.'engy')then
2448         call idmass(idproj,amproj)
2449         call idmass(idtarg,amtarg)
2450         if(engy.gt.0.)then
2451           xxengy=engy
2452         elseif(ecms.gt.0.)then
2453           xxengy=ecms
2454         elseif(elab.gt.0)then
2455           xxengy=sqrt( 2*elab*amtarg+amtarg**2+amproj**2 )
2456         elseif(pnll.gt.0)then
2457           xxengy=sqrt( 2*sqrt(pnll**2+amproj**2)*amtarg
2458      *                   +amtarg**2+amproj**2 )
2459         elseif(ekin.gt.0.)then
2460           xxelab=ekin+amproj
2461           xxengy=sqrt( 2*xxelab*amtarg+amtarg**2+amproj**2 )
2462         endif
2463         if(xxengy.lt.val1.or.xxengy.gt.val2)ifset=0
2464       endif
2465
2466            elseif(line(i:j).eq.'set'.and.ifset.eq.1)then
2467
2468       call utworn(line,j,ne)
2469       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'p/o-name p/o-value?'
2470       call utword(line,i,j,0)
2471       linex=line
2472       ix=i
2473       jx=j
2474       call utworn(line,j,ne)
2475       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'p/o-value?'
2476       call utword(line,i,j,0)
2477       read(line(i:j),*)val
2478 c       general
2479       if(linex(ix:jx).eq.'model') model=nint(val)
2480       if(linex(ix:jx).eq.'iquasiel') iquasiel=nint(val)
2481       if(linex(ix:jx).eq.'iversn')iversn=nint(val)
2482       if(linex(ix:jx).eq.'iappl' )iappl=nint(val)
2483       if(linex(ix:jx).eq.'nevent')nevent=nint(val)
2484       if(linex(ix:jx).eq.'nfull') nfull=nint(val)
2485       if(linex(ix:jx).eq.'nfreeze')nfreeze=nint(val)
2486       if(linex(ix:jx).eq.'ninicon')ninicon=nint(val)
2487       if(linex(ix:jx).eq.'egymin' )egymin=sngl(val)
2488       if(linex(ix:jx).eq.'egymax' )egymax=sngl(val)
2489 c       constants
2490       if(linex(ix:jx).eq.'ainfin')ainfin=sngl(val)
2491 c       printout options
2492       if(linex(ix:jx).eq.'iprmpt')iprmpt=nint(val)
2493       if(linex(ix:jx).eq.  'ish' )ish=nint(val)
2494       if(linex(ix:jx).eq.'ishsub')ishsub=nint(val)
2495       if(linex(ix:jx).eq.'irandm')irandm=nint(val)
2496       if(linex(ix:jx).eq.'irewch')irewch=nint(val)
2497       if(linex(ix:jx).eq.'iecho ')iecho =nint(val)
2498       if(linex(ix:jx).eq.'modsho')modsho=nint(val)
2499       if(linex(ix:jx).eq.'idensi')idensi=nint(val)
2500       if(linex(ix:jx).eq.'ishevt')ishevt=nint(val)
2501       if(linex(ix:jx).eq.'iwseed')iwseed=nint(val)
2502 c       fragmentation and decay
2503       if(linex(ix:jx).eq.'pud'   )pud   =sngl(val)
2504       if(linex(ix:jx).eq.'pude'  )pude  =sngl(val)
2505       if(linex(ix:jx).eq.'pudk'  )pudk  =sngl(val)
2506       if(linex(ix:jx).eq.'puds ' )puds   =sngl(val)
2507       if(linex(ix:jx).eq.'pudr ' )pudr   =sngl(val)
2508       if(linex(ix:jx).eq.'pudd ' )pudd   =sngl(val)
2509       if(linex(ix:jx).eq.'strcut' )strcut   =sngl(val)
2510       if(linex(ix:jx).eq.'diqcut' )diqcut   =sngl(val)
2511       if(linex(ix:jx).eq.'pdiqua')pdiqua=sngl(val)
2512       if(linex(ix:jx).eq.'pdiquae')pdiquae=sngl(val)
2513       if(linex(ix:jx).eq.'pdiquak')pdiquak=sngl(val)
2514       if(linex(ix:jx).eq.'pdiquar')pdiquar=sngl(val)
2515       if(linex(ix:jx).eq.'pdiquad')pdiquad=sngl(val)
2516       if(linex(ix:jx).eq.'ioptf ')ioptf =nint(val)
2517       if(linex(ix:jx).eq.'delrex')delrex=sngl(val)
2518       if(linex(ix:jx).eq.'ndecay')ndecay=nint(val)
2519       if(linex(ix:jx).eq.'maxres')maxres=nint(val)
2520       if(linex(ix:jx).eq.'pbreak')pbreak=sngl(val)
2521       if(linex(ix:jx).eq.'pbreake')pbreake=sngl(val)
2522       if(linex(ix:jx).eq.'pbreakk')pbreakk=sngl(val)
2523       if(linex(ix:jx).eq.'pbreakr')pbreakr=sngl(val)
2524       if(linex(ix:jx).eq.'pbreakd')pbreakd=sngl(val)
2525       if(linex(ix:jx).eq.'ptfra')ptfra=sngl(val)
2526       if(linex(ix:jx).eq.'ptfrae')ptfrae=sngl(val)
2527       if(linex(ix:jx).eq.'ptfrak')ptfrak=sngl(val)
2528       if(linex(ix:jx).eq.'ptfrakd')ptfrakd=sngl(val)
2529       if(linex(ix:jx).eq.'ptfrair')ptfrair=sngl(val)
2530       if(linex(ix:jx).eq.'ptfradr')ptfradr=sngl(val)
2531       if(linex(ix:jx).eq.'ptfrasr')ptfrasr=sngl(val)
2532       if(linex(ix:jx).eq.'ptfraqq')ptfraqq=sngl(val)
2533       if(linex(ix:jx).eq.'ptfrab')ptfrab=sngl(val)
2534       if(linex(ix:jx).eq.'aouni ')aouni=sngl(val)
2535 c       lepton-nucleon and e+e-
2536       if(linex(ix:jx).eq.'iolept')iolept=nint(val)
2537       if(linex(ix:jx).eq.'ydmin')ydmin=sngl(val)
2538       if(linex(ix:jx).eq.'ydmax')ydmax=sngl(val)
2539       if(linex(ix:jx).eq.'qdmin')qdmin=sngl(val)
2540       if(linex(ix:jx).eq.'qdmax')qdmax=sngl(val)
2541       if(linex(ix:jx).eq.'themin')themin=sngl(val)
2542       if(linex(ix:jx).eq.'themax')themax=sngl(val)
2543       if(linex(ix:jx).eq.'elomin')elomin=sngl(val)
2544       if(linex(ix:jx).eq.'engy' )engy=sngl(val)
2545       if(linex(ix:jx).eq.'elab' )elab=sngl(val)
2546       if(linex(ix:jx).eq.'ekin' )ekin=sngl(val)
2547       if(linex(ix:jx).eq.'ecms' )ecms=sngl(val)
2548       if(linex(ix:jx).eq.'elepti')elepti=sngl(val)
2549       if(linex(ix:jx).eq.'elepto')elepto=sngl(val)
2550       if(linex(ix:jx).eq.'angmue')angmue=sngl(val)
2551       if(linex(ix:jx).eq.'noebin')noebin=nint(val)
2552       if(linex(ix:jx).eq.'engmin')engmin=sngl(val)
2553       if(linex(ix:jx).eq.'engmax')engmax=sngl(val)
2554       if(linex(ix:jx).eq.'iologe')iologe=nint(val)
2555       if(linex(ix:jx).eq.'iologl')iologl=nint(val)
2556       if(linex(ix:jx).eq.'itflav')itflav=nint(val)
2557       if(linex(ix:jx).eq.'idisco')idisco=nint(val)
2558 c       hadron-hadron
2559       if(linex(ix:jx).eq. 'pnll' )pnll=sngl(val)
2560       if(linex(ix:jx).eq.'idproj')idprojin=nint(val)
2561       idproj=idprojin
2562       if(linex(ix:jx).eq.'idtarg')idtargin=nint(val)
2563       idtarg=idtargin
2564       if(idtarg.eq.0)idtarg=1120
2565       if(linex(ix:jx).eq.'ptq   ')ptq   =sngl(val)
2566       if(linex(ix:jx).eq.'rstrau(1)')rstrau(1)=sngl(val)
2567       if(linex(ix:jx).eq.'rstrad(1)')rstrad(1)=sngl(val)
2568       if(linex(ix:jx).eq.'rstras(1)')rstras(1)=sngl(val)
2569       if(linex(ix:jx).eq.'rstrau(2)')rstrau(2)=sngl(val)
2570       if(linex(ix:jx).eq.'rstrad(2)')rstrad(2)=sngl(val)
2571       if(linex(ix:jx).eq.'rstras(2)')rstras(2)=sngl(val)
2572       if(linex(ix:jx).eq.'rstrau(3)')rstrau(3)=sngl(val)
2573       if(linex(ix:jx).eq.'rstrad(3)')rstrad(3)=sngl(val)
2574       if(linex(ix:jx).eq.'rstras(3)')rstras(3)=sngl(val)
2575       if(linex(ix:jx).eq.'rstrau(4)')rstrau(4)=sngl(val)
2576       if(linex(ix:jx).eq.'rstrad(4)')rstrad(4)=sngl(val)
2577       if(linex(ix:jx).eq.'rstras(4)')rstras(4)=sngl(val)
2578       if(linex(ix:jx).eq.'rstrasi')rstrasi=sngl(val)
2579       if(linex(ix:jx).eq.'wgtval')wgtval=sngl(val)
2580       if(linex(ix:jx).eq.'wgtsea')wgtsea=sngl(val)
2581       if(linex(ix:jx).eq.'wgtdiq')wgtdiq=sngl(val)
2582       if(linex(ix:jx).eq.'wgtqqq(1)')wgtqqq(1)=sngl(val)
2583       if(linex(ix:jx).eq.'wgtqqq(2)')wgtqqq(2)=sngl(val)
2584       if(linex(ix:jx).eq.'wgtqqq(3)')wgtqqq(3)=sngl(val)
2585       if(linex(ix:jx).eq.'wgtqqq(4)')wgtqqq(4)=sngl(val)
2586       if(linex(ix:jx).eq.'wproj ')wproj =sngl(val)
2587       if(linex(ix:jx).eq.'wtarg ')wtarg =sngl(val)
2588       if(linex(ix:jx).eq.'wexcit')wexcit=sngl(val)
2589 c      if(linex(ix:jx).eq.'cutmsq')cutmsq=sngl(val)
2590       if(linex(ix:jx).eq.'cutmss')cutmss=sngl(val)
2591       if(linex(ix:jx).eq.'exmass')exmass=sngl(val)
2592       if(linex(ix:jx).eq.'iregge')iregge=nint(val)
2593       if(linex(ix:jx).eq.'isopom')isopom=nint(val)
2594       if(linex(ix:jx).eq.'ishpom')ishpom=nint(val)
2595       if(linex(ix:jx).eq.'iscreen')iscreen=nint(val)
2596       if(linex(ix:jx).eq.'isplit')isplit=nint(val)
2597       if(linex(ix:jx).eq.'irzptn')irzptn=nint(val)
2598       if(linex(ix:jx).eq.'irmdrop')irmdrop=nint(val)
2599       if(linex(ix:jx).eq.'nprmax')nprmax=nint(val)
2600       if(linex(ix:jx).eq.'iemspl')iemspl=nint(val)
2601       if(linex(ix:jx).eq.'intpol')intpol=nint(val)
2602       if(linex(ix:jx).eq.'isigma')isigma=nint(val)
2603       if(linex(ix:jx).eq.'iomega')iomega=nint(val)
2604       if(linex(ix:jx).eq.'isetcs')isetcs=nint(val)
2605       if(linex(ix:jx).eq.'iemsb' )iemsb= nint(val)
2606       if(linex(ix:jx).eq.'iemspm')iemspm=nint(val)
2607       if(linex(ix:jx).eq.'iemspx')iemspx=nint(val)
2608       if(linex(ix:jx).eq.'iemsse')iemsse=nint(val)
2609       if(linex(ix:jx).eq.'iemsdr')iemsdr=nint(val)
2610       if(linex(ix:jx).eq.'iemsrx')iemsrx=nint(val)
2611       if(linex(ix:jx).eq.'iemsi2')iemsi2=nint(val)
2612       if(linex(ix:jx).eq.'iemsi1')iemsi1=nint(val)
2613       if(linex(ix:jx).eq.'ioems' )ioems= nint(val)
2614       if(linex(ix:jx).eq.'ispacetime')ispacetime= nint(val)
2615 c       unified parameters
2616       if(linex(ix:jx).eq.'iclpro1')iclpro1=nint(val)
2617       if(linex(ix:jx).eq.'iclpro2')iclpro2=nint(val)
2618       if(linex(ix:jx).eq.'icltar1')icltar1=nint(val)
2619       if(linex(ix:jx).eq.'icltar2')icltar2=nint(val)
2620       if(linex(ix:jx).eq.'iclegy1')iclegy1=nint(val)
2621       if(linex(ix:jx).eq.'iclegy2')iclegy2=nint(val)
2622       if(linex(ix:jx).eq.'egylow')egylow=sngl(val)
2623       if(linex(ix:jx).eq.'alppom')alppom=sngl(val)
2624       if(linex(ix:jx).eq.'gamhad(1)')stop'gamhad(1) not allowed'
2625       if(linex(ix:jx).eq.'gamhad(2)')gamhad(2)=sngl(val)
2626       if(linex(ix:jx).eq.'gamhad(3)')stop'gamhad(3) not allowed'
2627       if(linex(ix:jx).eq.'gamhads(1)')gamhadsi(1)=sngl(val)
2628       if(linex(ix:jx).eq.'gamhads(2)')gamhadsi(2)=sngl(val)
2629       if(linex(ix:jx).eq.'gamhads(3)')gamhadsi(3)=sngl(val)
2630       if(linex(ix:jx).eq.'gamhads(4)')gamhadsi(4)=sngl(val)
2631       if(linex(ix:jx).eq.'gamtil')gamtil=sngl(val)
2632       if(linex(ix:jx).eq.'slopom')slopom=sngl(val)
2633       if(linex(ix:jx).eq.'slopoms')slopoms=sngl(val)
2634       if(linex(ix:jx).eq.'r2had(1)' )r2had(1)= sngl(val)
2635       if(linex(ix:jx).eq.'r2had(2)' )r2had(2)= sngl(val)
2636       if(linex(ix:jx).eq.'r2had(3)' )r2had(3)= sngl(val)
2637       if(linex(ix:jx).eq.'r2had(4)' )r2had(4)= sngl(val)
2638       if(linex(ix:jx).eq.'r2hads(1)' )r2hads(1)= sngl(val)
2639       if(linex(ix:jx).eq.'r2hads(2)' )r2hads(2)= sngl(val)
2640       if(linex(ix:jx).eq.'r2hads(3)' )r2hads(3)= sngl(val)
2641       if(linex(ix:jx).eq.'r2hads(4)' )r2hads(4)= sngl(val)
2642       if(linex(ix:jx).eq.'r3pom'   )r3pom= sngl(val)
2643       if(linex(ix:jx).eq.'r4pom'   )r4pom= sngl(val)
2644       if(linex(ix:jx).eq.'egyscr'  )egyscr= sngl(val)
2645       if(linex(ix:jx).eq.'epscrw'  )epscrw= sngl(val)
2646       if(linex(ix:jx).eq.'epscrx'  )epscrx= sngl(val)
2647       if(linex(ix:jx).eq.'epscrs'  )epscrs= sngl(val)
2648       if(linex(ix:jx).eq.'epscrh'  )epscrh= sngl(val)
2649       if(linex(ix:jx).eq.'epscrp'  )epscrp= sngl(val)
2650       if(linex(ix:jx).eq.'epscrv'  )epscrv= sngl(val)
2651       if(linex(ix:jx).eq.'gfactor' )gfactor= sngl(val)
2652       if(linex(ix:jx).eq.'gwidth'  )gwidth= sngl(val)
2653       if(linex(ix:jx).eq.'chad(1)' )chad(1)=  sngl(val)
2654       if(linex(ix:jx).eq.'chad(2)' )chad(2)=  sngl(val)
2655       if(linex(ix:jx).eq.'chad(3)' )chad(3)=  sngl(val)
2656       if(linex(ix:jx).eq.'chad(4)' )chad(4)=  sngl(val)
2657       if(linex(ix:jx).eq.'wdiff(1)')wdiff(1)=sngl(val)
2658       if(linex(ix:jx).eq.'wdiff(2)')wdiff(2)=sngl(val)
2659       if(linex(ix:jx).eq.'wdiff(3)')wdiff(3)=sngl(val)
2660       if(linex(ix:jx).eq.'wdiff(4)')wdiff(4)=sngl(val)
2661       if(linex(ix:jx).eq.'facdif')facdif=sngl(val)
2662       if(linex(ix:jx).eq.'rexndf')   rexndf=sngl(val)
2663       if(linex(ix:jx).eq.'rexndi(1)')rexndii(1)=sngl(val)
2664       if(linex(ix:jx).eq.'rexndi(2)')rexndii(2)=sngl(val)
2665       if(linex(ix:jx).eq.'rexndi(3)')rexndii(3)=sngl(val)
2666       if(linex(ix:jx).eq.'rexndi(4)')rexndii(4)=sngl(val)
2667       if(linex(ix:jx).eq.'rexdif(1)')rexdifi(1)=sngl(val)
2668       if(linex(ix:jx).eq.'rexdif(2)')rexdifi(2)=sngl(val)
2669       if(linex(ix:jx).eq.'rexdif(3)')rexdifi(3)=sngl(val)
2670       if(linex(ix:jx).eq.'rexdif(4)')rexdifi(4)=sngl(val)
2671       if(linex(ix:jx).eq.'rexres')rexres=sngl(val)
2672       if(linex(ix:jx).eq.'zbarfl')zbarfl=sngl(val)
2673       if(linex(ix:jx).eq.'alpreg')alpreg=sngl(val)
2674       if(linex(ix:jx).eq.'sloreg')sloreg=sngl(val)
2675       if(linex(ix:jx).eq.'gamreg')gamreg=sngl(val)
2676       if(linex(ix:jx).eq.'r2reg' )r2reg= sngl(val)
2677 c      if(linex(ix:jx).eq.'amhdibar')amhdibar= sngl(val)
2678       if(linex(ix:jx).eq.'ptdiff')ptdiff=sngl(val)
2679       if(linex(ix:jx).eq.'alppar')alppar=sngl(val)
2680       if(linex(ix:jx).eq.'alpsea')alpsea=sngl(val)
2681       if(linex(ix:jx).eq.'alpval')alpval=sngl(val)
2682       if(linex(ix:jx).eq.'alpdiq')alpdiq=sngl(val)
2683       if(linex(ix:jx).eq.'alplea(1)')alplea(1)=sngl(val)
2684       if(linex(ix:jx).eq.'alplea(2)')alplea(2)=sngl(val)
2685       if(linex(ix:jx).eq.'alplea(3)')alplea(3)=sngl(val)
2686       if(linex(ix:jx).eq.'alplea(4)')alplea(4)=sngl(val)
2687       if(linex(ix:jx).eq.'alpdif')alpdif=sngl(val)
2688       if(linex(ix:jx).eq.'alpdi')alpdi=sngl(val)
2689       if(linex(ix:jx).eq.'alpndi')alpndi=sngl(val)
2690       if(linex(ix:jx).eq.'ammsqq')ammsqq=sngl(val)
2691       if(linex(ix:jx).eq.'ammsqd')ammsqd=sngl(val)
2692       if(linex(ix:jx).eq.'ammsdd')ammsdd=sngl(val)
2693       if(linex(ix:jx).eq.'cumpom')cumpom=sngl(val)
2694       if(linex(ix:jx).eq.'reminv')reminv=sngl(val)
2695       if(linex(ix:jx).eq.'ptsend')ptsend=sngl(val)
2696       if(linex(ix:jx).eq.'ptsendi')ptsendi=sngl(val)
2697       if(linex(ix:jx).eq.'ptsemx')ptsemx=sngl(val)
2698       if(linex(ix:jx).eq.'zrminc')zrminc=sngl(val)
2699       if(linex(ix:jx).eq.'zrmmax')zrmmax=sngl(val)
2700       if(linex(ix:jx).eq.'edmaxi')edmaxi=sngl(val)
2701       if(linex(ix:jx).eq.'epmaxi')epmaxi=sngl(val)
2702       if(linex(ix:jx).eq.'zopinc')zopinc=sngl(val)
2703       if(linex(ix:jx).eq.'zopmax')zopmax=sngl(val)
2704       if(linex(ix:jx).eq.'zipinc')zipinc=sngl(val)
2705       if(linex(ix:jx).eq.'zipmax')zipmax=sngl(val)
2706       if(linex(ix:jx).eq.'zidinc')zidinc=sngl(val)
2707       if(linex(ix:jx).eq.'zidmax')zidmax=sngl(val)
2708       if(linex(ix:jx).eq.'zodinc')zodinc=sngl(val)
2709       if(linex(ix:jx).eq.'zodmax')zodmax=sngl(val)
2710       if(linex(ix:jx).eq.'zosinc')zosinc=sngl(val)
2711       if(linex(ix:jx).eq.'zosmax')zosmax=sngl(val)
2712       if(linex(ix:jx).eq.'zoeinc')zoeinc=sngl(val)
2713       if(linex(ix:jx).eq.'zoemax')zoemax=sngl(val)
2714       if(linex(ix:jx).eq.'zbcut') zbcut=sngl(val)
2715       if(linex(ix:jx).eq.'xmaxremn')xmaxremn=sngl(val)
2716       if(linex(ix:jx).eq.'xmaxdiff')xmaxdiff=sngl(val)
2717       if(linex(ix:jx).eq.'alpdro(1)')alpdro(1)=sngl(val)
2718       if(linex(ix:jx).eq.'alpdro(2)')alpdro(2)=sngl(val)
2719       if(linex(ix:jx).eq.'alpdro(3)')alpdro(3)=sngl(val)
2720       if(linex(ix:jx).eq.'amdrmax')amdrmax=sngl(val)
2721       if(linex(ix:jx).eq.'iodiba')iodiba=nint(val)
2722       if(linex(ix:jx).eq.'bidiba')bidiba=sngl(val)
2723
2724 c       hard pomeron parameters
2725       if(linex(ix:jx).eq.'q2min' )q2min=sngl(val)
2726       if(linex(ix:jx).eq.'q2ini' )q2ini=sngl(val)
2727       if(linex(ix:jx).eq.'q2fin' )q2fin=sngl(val)
2728       if(linex(ix:jx).eq.'betpom')betpom=sngl(val)
2729       if(linex(ix:jx).eq.'alpfom')alpfomi=sngl(val)
2730       if(linex(ix:jx).eq.'betfom')betfom=dble(val)
2731       if(linex(ix:jx).eq.'gamfom')gamfom=dble(val)
2732       if(linex(ix:jx).eq.'glusea')glusea=sngl(val)
2733       if(linex(ix:jx).eq.'factk' )factk=sngl(val)
2734       if(linex(ix:jx).eq.'naflav')naflav=nint(val)
2735       if(linex(ix:jx).eq.'pt2cut')pt2cut=sngl(val)
2736       if(linex(ix:jx).eq.'factgam')factgam=sngl(val)
2737       if(linex(ix:jx).eq.'delh')delh=sngl(val)
2738 c       nucleus-nucleus
2739       if(linex(ix:jx).eq.'iokoll')iokoll=nint(val)
2740       if(linex(ix:jx).eq.'laproj')laproj=nint(val)
2741       if(linex(ix:jx).eq.'maproj')maproj=nint(val)
2742       if(linex(ix:jx).eq.'latarg')latarg=nint(val)
2743       if(linex(ix:jx).eq.'matarg')matarg=nint(val)
2744       if(linex(ix:jx).eq.'core'  )core  =sngl(val)
2745 c      if(linex(ix:jx).eq.'ncolmx')ncolmx=nint(val)
2746       if(linex(ix:jx).eq.'fctrmx')fctrmx=sngl(val)
2747       if(linex(ix:jx).eq.'bmaxim')bmaxim=sngl(val)
2748       if(linex(ix:jx).eq.'bminim')bminim=sngl(val)
2749       if(linex(ix:jx).eq.'phimax')phimax=sngl(val)
2750       if(linex(ix:jx).eq.'phimin')phimin=sngl(val)
2751 c       rescattering parameters
2752       if(linex(ix:jx).eq.'iorsce')iorsce=nint(val)
2753       if(linex(ix:jx).eq.'iorsdf')iorsdf=nint(val)
2754       if(linex(ix:jx).eq.'iorshh')iorshh=nint(val)
2755       if(linex(ix:jx).eq.'iocluin')iocluin=nint(val)
2756       if(linex(ix:jx).eq.'ioquen')ioquen=nint(val)
2757       if(linex(ix:jx).eq.'hacore')hacore=sngl(val)
2758       if(linex(ix:jx).eq.'amimfs')amimfs=sngl(val)
2759       if(linex(ix:jx).eq.'amimel')amimel=sngl(val)
2760       if(linex(ix:jx).eq.'cepara')cepara=sngl(val)
2761       if(linex(ix:jx).eq.'dscale')dscale=sngl(val)
2762       if(linex(ix:jx).eq.'iceopt')iceopt=nint(val)
2763       if(linex(ix:jx).eq.'delamf')delamf=sngl(val)
2764       if(linex(ix:jx).eq.'deuamf')deuamf=sngl(val)
2765       if(linex(ix:jx).eq.'taurea')taurea=sngl(val)
2766       if(linex(ix:jx).eq.'nsegsu')nsegsu=nint(val)
2767       if(linex(ix:jx).eq.'nsegce')nsegce=nint(val)
2768       if(linex(ix:jx).eq.'ptclu') ptclu=sngl(val)
2769       if(linex(ix:jx).eq.'epscri(1)')epscri(1)=sngl(val)
2770       if(linex(ix:jx).eq.'epscri(3)')epscri(3)=sngl(val)
2771       if(linex(ix:jx).eq.'amsiac')amsiac=sngl(val)
2772       if(linex(ix:jx).eq.'amprif')amprif=sngl(val)
2773       if(linex(ix:jx).eq.'delvol')delvol=sngl(val)
2774       if(linex(ix:jx).eq.'deleps')deleps=sngl(val)
2775       if(linex(ix:jx).eq.'taumin')taumin=sngl(val)
2776       if(linex(ix:jx).eq.'deltau')deltau=sngl(val)
2777       if(linex(ix:jx).eq.'factau')factau=sngl(val)
2778       if(linex(ix:jx).eq.'numtau')numtau=nint(val)
2779       if(linex(ix:jx).eq.'dlzeta')dlzeta=sngl(val)
2780       if(linex(ix:jx).eq.'etafac')etafac=sngl(val)
2781       if(linex(ix:jx).eq.'facnuc')facnuc=sngl(val)
2782 c       urqmd
2783       if(linex(ix:jx).eq.'iurqmd')iurqmd=int(val)
2784 c       spherio
2785       if(linex(ix:jx).eq.'ispherio')ispherio=int(val)
2786 c       ico
2787       if(linex(ix:jx).eq.'cutico') cutico=sngl(val)
2788       if(linex(ix:jx).eq.'dssico') dssico=sngl(val)
2789       if(linex(ix:jx).eq.'icocore')icocore=int(val)
2790       if(linex(ix:jx).eq.'icotabm')icotabm=int(val)
2791       if(linex(ix:jx).eq.'icotabr')icotabr=int(val)
2792 c       droplet decay
2793       if(linex(ix:jx).eq.'dezzer')dezzer=sngl(val)
2794       if(linex(ix:jx).eq.'ioclude')ioclude=nint(val)
2795       if(linex(ix:jx).eq.'amuseg')amuseg=sngl(val)
2796       if(linex(ix:jx).eq.'yradmx')yradmx=sngl(val)
2797       if(linex(ix:jx).eq.'yradmi')yradmi=sngl(val)
2798       if(linex(ix:jx).eq.'yradpp')yradpp=sngl(val)
2799       if(linex(ix:jx).eq.'yradpi')yradpi=sngl(val)
2800       if(linex(ix:jx).eq.'yradpx')yradpx=sngl(val)
2801       if(linex(ix:jx).eq.'facecc')facecc=sngl(val)
2802       if(linex(ix:jx).eq.'rcoll' )rcoll= sngl(val)
2803       if(linex(ix:jx).eq.'ylongmx' )ylongmx= sngl(val)
2804       if(linex(ix:jx).eq.'bag4rt')bag4rt=sngl(val)
2805       if(linex(ix:jx).eq.'taunll')taunll=sngl(val)
2806 c       droplet specification
2807       if(linex(ix:jx).eq. 'keu'  )keu=nint(val)
2808       if(linex(ix:jx).eq. 'ked'  )ked=nint(val)
2809       if(linex(ix:jx).eq. 'kes'  )kes=nint(val)
2810       if(linex(ix:jx).eq. 'kec'  )kec=nint(val)
2811       if(linex(ix:jx).eq. 'keb'  )keb=nint(val)
2812       if(linex(ix:jx).eq. 'ket'  )ket=nint(val)
2813       if(linex(ix:jx).eq. 'tecm' )tecm=sngl(val)
2814       if(linex(ix:jx).eq. 'volu' )volu=sngl(val)
2815       if(linex(ix:jx).eq. 'vrad' )vrad=sngl(val)
2816       if(linex(ix:jx).eq. 'facts')facts=sngl(val)
2817       if(linex(ix:jx).eq. 'factb')factb=sngl(val)
2818       if(linex(ix:jx).eq. 'factq')factq=sngl(val)
2819       if(linex(ix:jx).eq.'inbxxx')inbxxx=sngl(val)
2820 c       metropolis
2821       if(linex(ix:jx).eq.'iospec')iospec=nint(val)
2822       if(linex(ix:jx).eq.'iocova')iocova=nint(val)
2823       if(linex(ix:jx).eq.'iopair')iopair=nint(val)
2824       if(linex(ix:jx).eq.'iozero')iozero=nint(val)
2825       if(linex(ix:jx).eq.'ioflac')ioflac=nint(val)
2826       if(linex(ix:jx).eq.'iostat')iostat=nint(val)
2827       if(linex(ix:jx).eq.'ioinco')ioinco=nint(val)
2828       if(linex(ix:jx).eq.'iograc')iograc=nint(val)
2829       if(linex(ix:jx).eq.'epsgc' )epsgc=sngl(val)
2830       if(linex(ix:jx).eq.'iocite')iocite=nint(val)
2831       if(linex(ix:jx).eq.'ioceau')ioceau=nint(val)
2832       if(linex(ix:jx).eq.'iociau')iociau=nint(val)
2833       if(linex(ix:jx).eq.'ioinct')ioinct=nint(val)
2834       if(linex(ix:jx).eq.'ioinfl')ioinfl=nint(val)
2835       if(linex(ix:jx).eq.'iowidn')iowidn=nint(val)
2836       if(linex(ix:jx).eq.'ionlat')ionlat=nint(val)
2837       if(linex(ix:jx).eq.'iomom')iomom=nint(val)
2838       if(linex(ix:jx).eq.'ioobsv')ioobsv=nint(val)
2839       if(linex(ix:jx).eq.'iosngl')iosngl=nint(val)
2840       if(linex(ix:jx).eq.'iorejz')iorejz=nint(val)
2841       if(linex(ix:jx).eq.'iompar')iompar=nint(val)
2842       if(linex(ix:jx).eq.'iozinc')iozinc=nint(val)
2843       if(linex(ix:jx).eq.'iozevt')iozevt=nint(val)
2844       if(linex(ix:jx).eq. 'nadd' )nadd=nint(val)
2845       if(linex(ix:jx).eq.'iterma')iterma=nint(val)
2846       if(linex(ix:jx).eq.'itermx')stop'STOP: set iterma, not itermx'
2847       if(linex(ix:jx).eq.'iterpr')iterpr=nint(val)
2848       if(linex(ix:jx).eq.'iterpl')iterpl=nint(val)
2849       if(linex(ix:jx).eq.'iternc')iternc=nint(val)
2850       if(linex(ix:jx).eq.'epsr'  )epsr=sngl(val)
2851       if(linex(ix:jx).eq.'keepr' )keepr=nint(val)
2852 c       strangelets
2853       if(linex(ix:jx).eq.'iopenu')iopenu=nint(val)
2854       if(linex(ix:jx).eq.'themas')themas=sngl(val)
2855 c       tests
2856       if(linex(ix:jx).eq.'iotst1')iotst1=nint(val)
2857       if(linex(ix:jx).eq.'iotst2')iotst2=nint(val)
2858       if(linex(ix:jx).eq.'iotst3')iotst3=nint(val)
2859       if(linex(ix:jx).eq.'iotst4')iotst4=nint(val)
2860 c       jpsi
2861       if(linex(ix:jx).eq.'jpsi  ')jpsi  =nint(val)
2862       if(linex(ix:jx).eq.'jpsifi')jpsifi=nint(val)
2863       if(linex(ix:jx).eq.'sigj  ')sigj  =sngl(val)
2864       if(linex(ix:jx).eq.'taumx ')taumx =sngl(val)
2865       if(linex(ix:jx).eq.'nsttau')nsttau=nint(val)
2866       if(linex(ix:jx).eq.'ijphis')ijphis=nint(val)
2867 c       analysis: intermittency, space-time, droplets, formation time
2868       if(linex(ix:jx).eq.'ymximi')ymximi=sngl(val)
2869       if(linex(ix:jx).eq.'imihis')imihis=nint(val)
2870       if(linex(ix:jx).eq.'isphis')isphis=nint(val)
2871       if(linex(ix:jx).eq.'iologb')iologb=nint(val)
2872       if(linex(ix:jx).eq.'ispall')ispall=nint(val)
2873       if(linex(ix:jx).eq.'wtmini')wtmini=sngl(val)
2874       if(linex(ix:jx).eq.'wtstep')wtstep=sngl(val)
2875       if(linex(ix:jx).eq.'iwcent')iwcent=nint(val)
2876       if(linex(ix:jx).eq.'iclhis')iclhis=nint(val)
2877       if(linex(ix:jx).eq.'iwtime')iwtime=nint(val)
2878       if(linex(ix:jx).eq.'wtimet')wtimet=sngl(val)
2879       if(linex(ix:jx).eq.'wtimei')wtimei=sngl(val)
2880       if(linex(ix:jx).eq.'wtimea')wtimea=sngl(val)
2881 c       other
2882       if(linex(ix:jx).eq.'gaumx ')gaumx =sngl(val)
2883       if(linex(ix:jx).eq.'nclean')nclean=nint(val)
2884       if(linex(ix:jx).eq.'istore')istore=nint(val)
2885       if(linex(ix:jx).eq.'ioidch')ioidch=nint(val)
2886       if(linex(ix:jx).eq.'iframe')iframe=nint(val)
2887       if(linex(ix:jx).eq.'jframe')jframe=nint(val)
2888       if(linex(ix:jx).eq.'labsys')stop'labsys no longer supported'
2889       if(linex(ix:jx).eq.'irescl')irescl=nint(val)
2890       if(linex(ix:jx).eq.'ifrade')ifrade=nint(val)
2891       if(linex(ix:jx).eq.'idecay')idecay=nint(val)
2892       if(linex(ix:jx).eq.'jdecay')jdecay=nint(val)
2893       if(linex(ix:jx).eq.'ntrymx')ntrymx=nint(val)
2894       if(linex(ix:jx).eq.'istmax')istmax=nint(val)
2895       if(linex(ix:jx).eq.'ionudi')ionudi=nint(val)
2896       if(linex(ix:jx).eq.'seedi')seedi =val
2897       if(linex(ix:jx).eq.'seedj')seedj =val
2898       if(linex(ix:jx).eq.'xvaria')xvaria='      '
2899       if(linex(ix:jx).eq.'xvaria')xvaria(1:j-i+1)=line(i:j)
2900       if(linex(ix:jx).eq.'yvaria')yvaria='      '
2901       if(linex(ix:jx).eq.'yvaria')yvaria(1:j-i+1)=line(i:j)
2902       if(linex(ix:jx).eq.'normal')normal=nint(val)
2903       if(linex(ix:jx).eq.'xminim')xminim=sngl(val)
2904       if(linex(ix:jx).eq.'xmaxim')xmaxim=sngl(val)
2905       if(linex(ix:jx).eq.'nrbins')nrbins=nint(val)
2906       if(linex(ix:jx).eq.'hisfac')hisfac=sngl(val)
2907       if(linex(ix:jx).eq.'xshift')xshift=sngl(val)
2908       if(linex(ix:jx).eq.'etacut')etacut=sngl(val)
2909       if(linex(ix:jx).eq.'xpar1' )xpar1=sngl(val)
2910       if(linex(ix:jx).eq.'xpar2' )xpar2=sngl(val)
2911       if(linex(ix:jx).eq.'xpar3' )xpar3=sngl(val)
2912       if(linex(ix:jx).eq.'xpar4' )xpar4=sngl(val)
2913       if(linex(ix:jx).eq.'xpar5' )xpar5=sngl(val)
2914       if(linex(ix:jx).eq.'xpar6' )xpar6=sngl(val)
2915       if(linex(ix:jx).eq.'xpar7' )xpar7=sngl(val)
2916       if(linex(ix:jx).eq.'xpar8' )xpar8=sngl(val)
2917       if(linex(ix:jx).eq.'irdmpr')irdmpr=nint(val)
2918 c       frame definitions
2919       if(linex(ix:jx).eq.'engy'.and.iframe.eq.0)iframe=11
2920       if(linex(ix:jx).eq.'ecms'.and.iframe.eq.0)iframe=11
2921       if(linex(ix:jx).eq.'elab'.and.iframe.eq.0)iframe=12
2922       if(linex(ix:jx).eq.'ekin'.and.iframe.eq.0)iframe=12
2923       if(linex(ix:jx).eq.'pnll'.and.iframe.eq.0)iframe=12
2924       if(linex(ix:jx).eq.'noebin'.and.dabs(val-1.d0).gt.1.d-10)iframe=11
2925
2926            elseif(line(i:j).eq.'set')then
2927
2928       call utword(line,i,j,0)
2929       call utword(line,i,j,0)
2930       ifset=1
2931
2932            elseif(line(i:j).eq.'stop')then  !same as return
2933
2934       if(nopen.ne.-1)then
2935       close(20+nopen)
2936       nopen=nopen-1
2937       if(nopen.eq.0.and.iprmpt.eq.-1)iprmpt=1
2938       endif
2939
2940            elseif(line(i:j).eq.'stopprogram')then
2941
2942       close(unit=ifcx)
2943       close(unit=ifhi)
2944       close(unit=ifdt)
2945       stop
2946
2947            elseif(line(i:j).eq.'EndEposInput')then
2948
2949       return
2950
2951            elseif(line(i:j).eq.'string')then
2952
2953       nstmax=nstmax+1
2954       ns=nstmax
2955       icinpu=0
2956       call utworn(line,j,ne)
2957       if(ne.eq.0.and.iprmpt.gt.0)
2958      *write(ifmt,'(a)')'string: prob icbac1 icbac2 icfor1 icfor2?'
2959       call utword(line,i,j,0)
2960       read(line(i:j),*)val
2961       prob(ns)=sngl(val)
2962       call utworn(line,j,ne)
2963       if(ne.eq.0.and.iprmpt.gt.0)
2964      *write(ifmt,'(a)')'string: icbac1 icbac2 icfor1 icfor2?'
2965       call utword(line,i,j,0)
2966       read(line(i:j),*)val
2967       icbac(ns,1)=nint(val)
2968       call utworn(line,j,ne)
2969       if(ne.eq.0.and.iprmpt.gt.0)
2970      *write(ifmt,'(a)')'string: icbac2 icfor1 icfor2?'
2971       call utword(line,i,j,0)
2972       read(line(i:j),*)val
2973       icbac(ns,2)=nint(val)
2974       call utworn(line,j,ne)
2975       if(ne.eq.0.and.iprmpt.gt.0)
2976      *write(ifmt,'(a)')'string: icfor1 icfor2?'
2977       call utword(line,i,j,0)
2978       read(line(i:j),*)val
2979       icfor(ns,1)=nint(val)
2980       call utworn(line,j,ne)
2981       if(ne.eq.0.and.iprmpt.gt.0)
2982      *write(ifmt,'(a)')'string: icfor2?'
2983       call utword(line,i,j,0)
2984       read(line(i:j),*)val
2985       icfor(ns,2)=nint(val)
2986
2987            elseif(line(i:j).eq.'kinks')then
2988
2989       nptl=0
2990 ctp290806      nrow=0
2991       nel=0
2992  10   continue
2993       call utword(line,i,j,0)
2994       if(line(i:j).eq.'endkinks')goto 12
2995       nel=nel+1
2996 ctp290806      nrow=1+(nel-1)/4
2997       nc=mod(nel-1,4)+1
2998       read(line(i:j),*)a
2999       if(nc.eq.1)nptl=nptl+1
3000       if(nc.eq.1)idptl(nptl)=nint(a)
3001       if(nc.eq.2)pptl(1,nptl)=a
3002       if(nc.eq.3)pptl(2,nptl)=a
3003       if(nc.eq.4)then
3004         pptl(3,nptl)=a
3005         istptl(nptl)=20
3006         pptl(4,nptl)=sqrt(pptl(3,nptl)**2+pptl(2,nptl)**2
3007      $       +pptl(1,nptl)**2)
3008       endif
3009       goto10
3010  12   continue
3011
3012            elseif(line(i:j).eq.'record')then
3013
3014       call utworn(line,j,ne)
3015 c      if(ne.eq.0.and.iprmpt.gt.0)
3016 c     *     write(ifmt,'(a)')'kinks: icbac1 icbac2 icfor1 icfor2?'
3017       call utword(line,i,j,0)
3018       if(line(i:j).eq.'event')ir=1
3019       if(line(i:j).eq.'particle')ir=2
3020       maxrec(ir)=0
3021  20   call utworn(line,j,ne)
3022 c      if(ne.eq.0.and.iprmpt.gt.0)
3023 c     *     write(6,'(a)')'<kinks-data (px-py-pz)>? (End=endkinks)'
3024       call utword(line,i,j,0)
3025       if(line(i:j).eq.'endrecord')then
3026          goto 22
3027       endif
3028       maxrec(ir)=maxrec(ir)+1
3029       irecty(maxrec(ir),ir)=0
3030       if(ir.eq.1)then
3031         if(line(i:j).eq.'0') irecty(maxrec(ir),ir)=0
3032         if(line(i:j).eq.'nevt') irecty(maxrec(ir),ir)=1
3033         if(line(i:j).eq.'nptl') irecty(maxrec(ir),ir)=2
3034         if(line(i:j).eq.'b') irecty(maxrec(ir),ir)=3
3035         if(line(i:j).eq.'phi') irecty(maxrec(ir),ir)=4
3036         if(line(i:j).eq.'ncol') irecty(maxrec(ir),ir)=5
3037         if(line(i:j).eq.'pmx') irecty(maxrec(ir),ir)=6
3038         if(line(i:j).eq.'egy') irecty(maxrec(ir),ir)=7
3039         if(line(i:j).eq.'npj') irecty(maxrec(ir),ir)=8
3040         if(line(i:j).eq.'ntg') irecty(maxrec(ir),ir)=9
3041         if(line(i:j).eq.'npn') irecty(maxrec(ir),ir)=10
3042         if(line(i:j).eq.'npp') irecty(maxrec(ir),ir)=11
3043         if(line(i:j).eq.'ntn') irecty(maxrec(ir),ir)=12
3044         if(line(i:j).eq.'ntp') irecty(maxrec(ir),ir)=13
3045         if(line(i:j).eq.'jpn') irecty(maxrec(ir),ir)=14
3046         if(line(i:j).eq.'jpp') irecty(maxrec(ir),ir)=15
3047         if(line(i:j).eq.'jtn') irecty(maxrec(ir),ir)=16
3048         if(line(i:j).eq.'jtp') irecty(maxrec(ir),ir)=17
3049         if(line(i:j).eq.'amproj') irecty(maxrec(ir),ir)=20
3050         if(line(i:j).eq.'amtarg') irecty(maxrec(ir),ir)=21
3051       else
3052         if(line(i:j).eq.'0') irecty(maxrec(ir),ir)=0
3053         if(line(i:j).eq.'i') irecty(maxrec(ir),ir)=1
3054         if(line(i:j).eq.'id') irecty(maxrec(ir),ir)=2
3055         if(line(i:j).eq.'p1') irecty(maxrec(ir),ir)=3
3056         if(line(i:j).eq.'p2') irecty(maxrec(ir),ir)=4
3057         if(line(i:j).eq.'p3') irecty(maxrec(ir),ir)=5
3058         if(line(i:j).eq.'p4') irecty(maxrec(ir),ir)=6
3059         if(line(i:j).eq.'p5') irecty(maxrec(ir),ir)=7
3060         if(line(i:j).eq.'fa') irecty(maxrec(ir),ir)=8
3061         if(line(i:j).eq.'mo') irecty(maxrec(ir),ir)=9
3062         if(line(i:j).eq.'st') irecty(maxrec(ir),ir)=10
3063         if(line(i:j).eq.'x1') irecty(maxrec(ir),ir)=11
3064         if(line(i:j).eq.'x2') irecty(maxrec(ir),ir)=12
3065         if(line(i:j).eq.'x3') irecty(maxrec(ir),ir)=13
3066         if(line(i:j).eq.'x4') irecty(maxrec(ir),ir)=14
3067         if(line(i:j).eq.'idfa') irecty(maxrec(ir),ir)=15
3068         if(line(i:j).eq.'idmo') irecty(maxrec(ir),ir)=16
3069         if(line(i:j).eq.'p') irecty(maxrec(ir),ir)=17
3070         if(line(i:j).eq.'x') irecty(maxrec(ir),ir)=18
3071         if(line(i:j).eq.'dez') irecty(maxrec(ir),ir)=19
3072         if(line(i:j).eq.'c1') irecty(maxrec(ir),ir)=21
3073         if(line(i:j).eq.'c2') irecty(maxrec(ir),ir)=22
3074         if(line(i:j).eq.'ty') irecty(maxrec(ir),ir)=23
3075       endif
3076       if(irecty(maxrec(ir),ir).eq.0)then
3077         write(*,*) 'unknown variable ',line(i:j)
3078         stop
3079       endif
3080       goto 20
3081  22   continue
3082
3083            elseif(line(i:j).eq.'switch')then
3084
3085       call utworn(line,j,ne)
3086       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'option on/off?'
3087       call utword(line,i,j,0)
3088       call utworn(line,j,ne)
3089       if(ne.eq.0.and.iprmpt.gt.0)write(ifmt,'(a)')'on/off?'
3090         if(line(i:j).eq.'droplet')then
3091       call utword(line,i,j,0)
3092       if(line(i:j).eq.'on' )iorsdf=1
3093       if(line(i:j).eq.'off')iorsdf=0
3094         elseif(line(i:j).eq.'cascade')then
3095       call utword(line,i,j,0)
3096 c      if(line(i:j).eq.'on' )iorsce=1
3097       if(line(i:j).eq.'on' )iorshh=1
3098 c      if(line(i:j).eq.'off')iorsce=0
3099       if(line(i:j).eq.'off')iorshh=0
3100         elseif(line(i:j).eq.'soft')then
3101       call utword(line,i,j,0)
3102       if(line(i:j).eq.'on' )isopom=1
3103       if(line(i:j).eq.'off')isopom=0
3104         elseif(line(i:j).eq.'hard')then
3105       call utword(line,i,j,0)
3106       if(line(i:j).eq.'on' )ishpom=1
3107       if(line(i:j).eq.'off')ishpom=0
3108         elseif(line(i:j).eq.'splitting')then
3109       call utword(line,i,j,0)
3110       if(line(i:j).eq.'on' )isplit=1
3111       if(line(i:j).eq.'on' )iscreen=1
3112       if(line(i:j).eq.'off')isplit=0
3113       if(line(i:j).eq.'off')iscreen=0
3114         elseif(line(i:j).eq.'fusion')then
3115       call utword(line,i,j,0)
3116       if(line(i:j).eq.'on' )iorsce=0
3117       if(line(i:j).eq.'on' )iorsdf=3
3118       if(line(i:j).eq.'on' )iorshh=0
3119       if(line(i:j).eq.'off')iorsce=0
3120       if(line(i:j).eq.'off')iorsdf=0
3121       if(line(i:j).eq.'off')iorshh=0
3122         elseif(line(i:j).eq.'urqmd')then
3123       call utword(line,i,j,0)
3124       if(line(i:j).eq.'on' ) iurqmd=1
3125       if(line(i:j).eq.'off') iurqmd=0
3126         elseif(line(i:j).eq.'spherio')then
3127       call utword(line,i,j,0)
3128       if(line(i:j).eq.'on' ) ispherio=1
3129       if(line(i:j).eq.'off') ispherio=0
3130         elseif(line(i:j).eq.'decay')then
3131       call utword(line,i,j,0)
3132       if(line(i:j).eq.'on' ) ndecay=0
3133       if(line(i:j).eq.'off') ndecay=1
3134       if(line(i:j).eq.'on' ) idecay=1
3135       if(line(i:j).eq.'off') idecay=0
3136         elseif(line(i:j).eq.'clusterdecay')then
3137       call utword(line,i,j,0)
3138       if(line(i:j).eq.'on' ) jdecay=1
3139       if(line(i:j).eq.'off') jdecay=0
3140         elseif(line(i:j).eq.'fragdecay')then
3141       call utword(line,i,j,0)
3142       if(line(i:j).eq.'on' ) ifrade=1
3143       if(line(i:j).eq.'off') ifrade=0
3144         elseif(line(i:j).eq.'icocore')then
3145           stop'switch icocore not supported any more.    '
3146         endif
3147
3148            elseif(line(i:j).eq.'idchoice')then
3149
3150       call utword(line,i,j,0)
3151       if(line(i:j).eq.'nxs')then
3152         ioidch=1
3153       elseif(line(i:j).eq.'pdf')then
3154         ioidch=2
3155       else
3156         stop'invalid idchoice.     '
3157       endif
3158
3159            elseif(line(i:j).eq.'make')then
3160
3161       call utword(line,i,j,0)
3162       if(line(i:j).eq.'icotable')icotabm=1
3163
3164            elseif(line(i:j).eq.'read')then
3165
3166       call utword(line,i,j,0)
3167       if(line(i:j).eq.'icotable')icotabr=1
3168
3169            elseif(line(i:j).eq.'output')then
3170
3171       call utword(line,i,j,0)
3172       if(line(i:j).eq.'epos' )     istore=1
3173       if(line(i:j).eq.'osc1997a' )  istore=2
3174       if(line(i:j).eq.'osc1999a' )  istore=3
3175       if(line(i:j).eq.'ustore' )    istore=5
3176
3177            elseif(line(i:j).eq.'model')then
3178
3179       call utword(line,i,j,0)
3180       if(line(i:j).eq.'epos')then
3181         model=1
3182       else
3183         nij=j-i+1
3184         if(nij.gt.20)stop'cmodel too small'
3185         cmodel(1:nij)=line(i:j)
3186         cmodel(nij+1:nij+1)=' '
3187         call NumberModel(cmodel,model)
3188       endif
3189       if(iappl.ne.1.and.iappl.ne.3.and.model.ne.1)
3190      &call utstop('Application not possible with this model&')
3191
3192            elseif(line(i:j).eq.'trigger')then
3193
3194       call utword(line,i,j,0)
3195       ntc=1
3196       if(line(i:j).eq.'or'.or.line(i:j).eq.'contr')then
3197         call utword(line,i,j,0)
3198         read(line(i:j),*)ztc
3199         ntc=nint(ztc)
3200         call utword(line,i,j,1)
3201       endif
3202       do n=1,ntc
3203         if(n.ne.1)call utword(line,i,j,0)
3204         call utword(line,i,j,0)
3205         call utword(line,i,j,0)
3206       enddo
3207
3208            elseif(line(i:j).eq.'noerrorbut')then
3209
3210       call utword(line,i,j,0)
3211
3212            elseif(line(i:j).eq.'b')then
3213
3214       nbarray=nbarray+1
3215       call utword(line,i,j,0)
3216       read(line(i:j),*)val
3217       barray(nbarray)=val
3218
3219            elseif(line(i:j).eq.'message')then
3220
3221       call utword(line,i,j,0)
3222       if(nopen.eq.-1)then      !only write in second read
3223       write(ifmt,'(a,$)')line(i:j)
3224       endif
3225
3226            elseif(line(i:j).eq.'endmessage')then
3227
3228       if(nopen.eq.-1)then      !only write in second read
3229       write(ifmt,'(a)')' '
3230       endif
3231
3232            elseif(line(i:j).eq.'write'.or.line(i:j).eq.'writex')then
3233
3234       if(line(i:j).eq.'write')ii=1
3235       if(line(i:j).eq.'writex')ii=2
3236       call utword(line,i,j,0)
3237       idol=0
3238       if(line(i:j).eq.'$')then
3239        idol=1
3240        call utword(line,i,j,0)
3241       endif
3242       !write: only write in second read; writex: only write in first read
3243       if(ii.eq.1.and.nopen.eq.-1.or.ii.eq.2.and.nopen.ne.-1)then
3244        if(j-6.ge.i)then
3245         do l=i,j-6
3246          if(line(l:l+6).eq.'$iversn')then
3247           if(mod(iversn,100).le.9)then
3248            write(line(l:l+6),'(i1,a,i1,3x)')
3249      *     int(iversn/100),'.0',mod(iversn,100)
3250           else
3251            write(line(l:l+6),'(i1,a,i2,3x)')
3252      *     int(iversn/100),'.',mod(iversn,100)
3253           endif
3254          endif
3255         enddo
3256        endif
3257        if(j-8.ge.i)then
3258         do l=i,j-8
3259          if(line(l:l+8).eq.'$xxxyield')then
3260           write(line(l:l+8),'(f8.6,1x)')yield
3261          endif
3262          if(line(l:l+7).eq.'$xxyield')then
3263           write(line(l:l+8),'(f7.5,1x)')yield
3264          endif
3265          if(line(l:l+6).eq.'$xyield')then
3266           write(line(l:l+8),'(f6.4,1x)')yield
3267          endif
3268         enddo
3269        endif
3270        if(j-5.ge.i)then
3271         do l=i,j-5
3272          if(line(l:l+5).eq.'$yield')then
3273           if(yield.lt.1.0)then
3274            write(line(l:l+5),'(f5.3,1x)')yield
3275           elseif(yield.lt.100.)then
3276            write(line(l:l+5),'(f5.2,1x)')yield
3277           elseif(yield.lt.1000.)then
3278            write(line(l:l+5),'(f5.1,1x)')yield
3279           elseif(yield.lt.10000.)then
3280            write(line(l:l+5),'(f6.1)')yield
3281           else
3282            write(line(l:l+5),'(i6)')nint(yield)
3283           endif
3284          endif
3285         enddo
3286         do l=i,j-5
3287          if(line(l:l+5).eq.'$averg')then
3288           if(averg.lt.1.0)then
3289            write(line(l:l+5),'(f5.3,1x)')averg
3290           elseif(averg.lt.100.)then
3291            write(line(l:l+5),'(f5.2,1x)')averg
3292           elseif(averg.lt.1000.)then
3293            write(line(l:l+5),'(f5.1,1x)')averg
3294           elseif(averg.lt.10000.)then
3295            write(line(l:l+5),'(f6.1)')averg
3296           else
3297            write(line(l:l+5),'(i6)')nint(averg)
3298           endif
3299          endif
3300         enddo
3301         do l=i,j-5
3302          if(line(l:l+5).eq.'$sigma')then
3303           if(sigma.lt.1.0)then
3304            write(line(l:l+5),'(f5.3,1x)')sigma
3305           elseif(sigma.lt.100.)then
3306            write(line(l:l+5),'(f5.2,1x)')sigma
3307           elseif(sigma.lt.1000.)then
3308            write(line(l:l+5),'(f5.1,1x)')sigma
3309           elseif(sigma.lt.10000.)then
3310            write(line(l:l+5),'(f6.1)')sigma
3311           else
3312            write(line(l:l+5),'(i6)')nint(sigma)
3313           endif
3314          endif
3315         enddo
3316        endif
3317        if(idol.eq.0)then
3318         write(ifhi,'(a)')line(i:j)
3319        else
3320         write(ifhi,'(a,a,$)')line(i:j),' '
3321        endif
3322       endif
3323
3324            elseif(line(i:j).eq.'nozero')then
3325
3326       nozero=1
3327
3328            elseif(line(i:j).eq.'ibmin')then
3329
3330       call utword(line,i,j,0)
3331       read(line(i:j),*)val
3332       ibmin=nint(val)
3333
3334            elseif(line(i:j).eq.'ibmax')then
3335
3336       call utword(line,i,j,0)
3337       read(line(i:j),*)val
3338       ibmax=nint(val)
3339
3340             elseif(line(i:j).eq.'writearray'
3341      $     .or.line(i:j).eq.'writehisto')then
3342
3343       if(nopen.eq.-1)then !second run
3344        ih=0
3345        if(line(i:j).eq.'writearray') ih=1
3346        call utword(line,i,j,0)
3347        if(line(i:j).eq.'s')then
3348         call utword(line,i,j,0)
3349         linex=line
3350         ix=i
3351         jx=j
3352         call utword(line,i,j,0)
3353         if(linex(ix:jx).eq.'inicon')stop'error 060307'
3354        else
3355         ioint=0
3356         iocontr=0
3357         if(line(i:j).eq.'int')then
3358          ioint=1
3359          call utword(line,i,j,0)
3360         endif
3361         if(line(i:j).eq.'contr')then
3362          iocontr=1
3363          call utword(line,i,j,0)
3364          read(line(i:j),*)val
3365          ncontr=nint(val)
3366          call utword(line,i,j,0)
3367         endif
3368         read(line(i:j),*)val
3369         nco=nint(val)
3370         if(ih.eq.1)write(ifhi,'(a,i3)')'array',nco
3371         if(ioint.eq.0)then
3372          sum=0
3373          averg=0
3374          do k=1,nrbins
3375           if(iocontr.eq.0.and.ionoerr.eq.0)then
3376             ar3=ar(k,3)
3377             ar4=ar(k,4)
3378           elseif(ionoerr.eq.1)then
3379             ar3=ar(k,3)
3380             ar4=ar(k,4)
3381           elseif(ionoerr.eq.2)then
3382             ar3=ar(k,3)
3383             ar4=ar(k,4)
3384             ar5=ar(k,5)
3385           else
3386             ar3=ary(k,ncontr)
3387             ar4=ardy(k,ncontr)
3388           endif
3389           iok=1
3390           if(k.lt.ibmin.or.k.gt.ibmax)iok=0
3391           sum=sum+ar3
3392           averg=averg+ar(k,1)*ar3
3393           if(nco.eq.2)then
3394             if(nozero.eq.1.and.ar3.eq.0.)iok=0
3395             if(iok.eq.1)write(ifhi,'(3e12.4)')ar(k,1),ar3
3396           elseif(nco.eq.3)then
3397             if(nozero.eq.1.and.ar3.eq.0..and.ar4.eq.0.)iok=0
3398             if(iok.eq.1)write(ifhi,'(3e12.4)')ar(k,1),ar3,ar4
3399           elseif(nco.eq.4)then
3400           if(nozero.eq.1.and.ar3.eq.0..and.ar4.eq.0..and.ar5.eq.0.)iok=0
3401             if(iok.eq.1)write(ifhi,'(4e12.4)')ar(k,1),ar3,ar4,ar5
3402           endif
3403          enddo
3404          if(sum.gt.0.)averg=averg/sum
3405         else
3406          sum=0.
3407          sum2=0.
3408          sum3=0.
3409          err2=0.
3410          do k=1,nrbins
3411           if(iocontr.eq.0.and.ionoerr.eq.0)then
3412             ar3=ar(k,3)
3413             ar4=ar(k,4)
3414           elseif(ionoerr.eq.1)then
3415             ar3=ar(k,3)
3416             ar4=ar(k,4)
3417           elseif(ionoerr.eq.2)then
3418             ar3=ar(k,3)
3419             ar4=ar(k,4)
3420             ar5=ar(k,5)
3421           else
3422             ar3=ary(k,ncontr)
3423             ar4=ardy(k,ncontr)
3424           endif
3425           sum=sum+ar3*(ar(2,1)-ar(1,1))
3426           if(nco.eq.2)write(ifhi,'(3e12.4)')ar(k,1),sum
3427           if(ionoerr.eq.0)then
3428             err2=err2+(ar4*(ar(2,1)-ar(1,1)))**2
3429             if(nco.eq.3)write(ifhi,'(3e12.4)')ar(k,1),sum,sqrt(err2)
3430           elseif(ionoerr.eq.1)then
3431             sum2=sum2+(ar4*(ar(2,1)-ar(1,1)))
3432             if(nco.eq.3)write(ifhi,'(3e12.4)')ar(k,1),sum,sum2
3433           elseif(ionoerr.eq.2)then
3434             sum2=sum2+(ar4*(ar(2,1)-ar(1,1)))
3435             sum3=sum3+(ar5*(ar(2,1)-ar(1,1)))
3436             if(nco.eq.3)write(ifhi,'(3e12.4)')ar(k,1),sum,sum2
3437             if(nco.eq.4)write(ifhi,'(3e12.4)')ar(k,1),sum,sum2,sum3
3438           endif
3439          enddo
3440         endif
3441         if(ih.eq.1)write(ifhi,'(a)')'  endarray'
3442        endif
3443       else !nopen .ge. 0 -- first run
3444         call utword(line,i,j,0)
3445         if(line(i:j).eq.'s')then
3446           call utword(line,i,j,0)
3447           call utword(line,i,j,0)
3448         elseif(line(i:j).eq.'int')then
3449           call utword(line,i,j,0)
3450         elseif(line(i:j).eq.'contr')then
3451           call utword(line,i,j,0)
3452           call utword(line,i,j,0)
3453         endif
3454       endif
3455       nozero=0
3456       ibmin=1
3457       ibmax=1e8
3458
3459            else
3460
3461       write(ifmt,'(a,a,a)')'command "',line(i:j),'" not found'
3462       j=160
3463       stop
3464
3465            endif
3466
3467       i=j+1
3468       goto 1
3469
3470       end
3471
3472 c-----------------------------------------------------------------------
3473       subroutine aseed(modus)
3474 c-----------------------------------------------------------------------
3475
3476       include 'epos.inc'
3477       double precision seedf
3478       call utpri('aseed ',ish,ishini,3)
3479
3480       call ranfgt(seedf)
3481       if(iwseed.eq.1)then
3482         if(nrevt.eq.0)then
3483           write(ifmt,'(a,d27.16)')'seedj:',seedf
3484         elseif(mod(nrevt,modsho).eq.0)then
3485           if(modus.eq.1)
3486      *   write(ifmt,'(a,i10,10x,a,d27.16)')'nrevt:',nrevt,'seedf:',seedf
3487           if(modus.eq.2)
3488      *   write(ifmt,'(a,d27.16)')'seed:',seedf
3489         endif
3490       endif
3491       seedc=seedf
3492
3493       call utprix('aseed ',ish,ishini,3)
3494       return
3495       end
3496
3497 c-----------------------------------------------------------------------
3498       subroutine aseedi
3499 c-----------------------------------------------------------------------
3500
3501       include 'epos.inc'
3502       double precision seedini
3503       call utpri('aseedi',ish,ishini,3)
3504
3505       call ranfgt(seedini)
3506
3507       if(ish.ge.1)write(ifmt,'(a,d26.15)')'seedi:',seedini
3508
3509       call utprix('aseedi',ish,ishini,3)
3510       return
3511       end
3512
3513 c$$$c-----------------------------------------------------------------------
3514 c$$$        subroutine aseed(modus)        !Flush ????
3515 c$$$c-----------------------------------------------------------------------
3516 c$$$
3517 c$$$      include 'epos.inc'
3518 c$$$      double precision seedf
3519 c$$$      call utpri('aseed',ish,ishini,4)
3520 c$$$
3521 c$$$      call ranfgt(seedf)
3522 c$$$      if(modus.eq.2)then
3523 c$$$        write(ifmt,'(a,d26.15)')'seed:',seedf
3524 c$$$      elseif(modus.eq.1)then
3525 c$$$        if(mod(nrevt,modsho).eq.0)then
3526 c$$$          write(ifmt,100)'nrevt:',nrevt,'seedf:',seedf
3527 c$$$          call flush(ifmt)
3528 c$$$        endif
3529 c$$$      endif
3530 c$$$      seedc=seedf
3531 c$$$
3532 c$$$  100 format(a,i10,10x,a,d26.15)
3533 c$$$      call utprix('aseed',ish,ishini,4)
3534 c$$$      return
3535 c$$$      end
3536 c$$$
3537 c-----------------------------------------------------------------------
3538       subroutine astati
3539 c-----------------------------------------------------------------------
3540
3541       include 'epos.inc'
3542       common/geom/rmproj,rmtarg,bmax,bkmx
3543       common/ghecsquel/anquasiel,iquasiel
3544
3545       call utpri('astati',ish,ishini,1)
3546       if(ish.ge.1.and.iappl.eq.1.)then
3547         if(abs(accept+reject).gt.1.e-5)write(ifch,'(a,f9.5)')
3548      *'EMS acc.rate:',accept/(accept+reject)
3549         if(antot.ne.0.)write(ifch,*)'initial soft,hard(%)'
3550      *                           ,ansf/antot*100.
3551      *                           ,ansh/antot*100.,' of' ,antot
3552         if(antotf.ne.0.)write(ifch,*)'final soft,hard(%)'
3553      *                           ,ansff/antotf*100.
3554      *                           ,anshf/antotf*100.,' of' ,antotf
3555         if(antotre.ne.0.)write(ifch,*)
3556      *                  'droplet,string(+d),reson(+d), (had)(%) '
3557         if(antotre.ne.0.)write(ifch,*)'     '
3558      *                           ,andropl/antotre*100.
3559      *                           ,anstrg0/antotre*100.
3560      *                           ,'(',anstrg1/antotre*100.,') '
3561      *                           ,anreso0/antotre*100.
3562      *                           ,'(',anreso1/antotre*100.,') '
3563         if(antotre.ne.0.)write(ifch,*)'     '
3564      *             ,' (',anghadr/antotre*100.,')',' of' ,antotre
3565        if(pp4ini.gt.0.)write(ifch,*)'Energy loss',(pp4ini-pp4max)/pp4ini
3566       write(ifch,*)'ine cross section:',sigineex
3567       write(ifch,*)'diffr cross section:',sigdifex
3568       write(ifch,*)'SD cross section:',sigsdex
3569 c      if(model.eq.3)write(ifch,*)'quasi-elastic cross section:'
3570 c     &,anquasiel/float(ntevt)*a*10
3571          endif
3572
3573 c$$$      call testconex(3)
3574
3575       if(iprmpt.le.0)goto1000
3576
3577       write(ifch,'(77a1)')('-',i=1,77)
3578       write(ifch,'(a)')'statistics'
3579       write(ifch,'(77a1)')('-',i=1,77)
3580       write(ifch,'(a,i6)')'nr of messages:',imsg
3581       write(ifch,'(a,i8)')'maximum nptl:',nptlu
3582       write(ifch,'(77a1)')('-',i=1,77)
3583
3584 1000  continue
3585       call utprix('astati',ish,ishini,1)
3586
3587       return
3588       end
3589
3590 c-----------------------------------------------------------------------
3591       subroutine atitle
3592 c-----------------------------------------------------------------------
3593
3594       include 'epos.inc'
3595
3596       write(ifmt,'(67a1/a1,8x,a,f5.2,5x,a,7x,a1/a,22x,a,10x,a1/67a1)')
3597      *('#',l=1,68),'EPOS',iversn/100.
3598      *,'K. WERNER, T. PIEROG, S. PORTEBOEUF.','#'
3599      *,'#','Contact: werner@subatech.in2p3.fr'
3600      *,('#',l=1,68)
3601       if(iversn.eq.iverso)return
3602       write(ifmt,'(a,8x,a,10x,a/a,5x,a,6x,a/67a1)')
3603      * '#','WARNING: This is a non-official beta version!!!','#'
3604      *,'#','Do not publish results without contacting the authors.','#'
3605      *,('#',l=1,67)
3606
3607       return
3608       end
3609
3610 c-----------------------------------------------------------------------
3611       subroutine avehep
3612 c-----------------------------------------------------------------------
3613
3614       include 'epos.inc'
3615
3616       call utpri('avehep',ish,ishini,4)
3617
3618
3619       call utprix('avehep',ish,ishini,4)
3620       end
3621
3622
3623 c-----------------------------------------------------------------------
3624       subroutine aepos(nin)
3625 c-----------------------------------------------------------------------
3626
3627       include 'epos.inc'
3628       include 'epos.incems'
3629       double precision eppass,etpass
3630       common/emnpass/eppass(mamx,4),etpass(mamx,4)
3631       integer iutime(5)
3632       call utpri('aepos',ish,ishini,4)
3633
3634       call timer(iutime)
3635       timeini=iutime(3)+float(iutime(4))/1000.
3636
3637       if(ish.ge.2)then
3638       call alist('start event&',0,0)
3639       write(ifch,*)'event number:',nrevt+1
3640       endif
3641
3642 c if random sign for projectile, set it here
3643       if(irdmpr.ne.0.and.maproj.eq.1)
3644      &   idproj=idprojin*(1-2*int(rangen()+0.5d0))
3645
3646 c for Air target, set the target nucleus
3647       if(idtargin.eq.0.and.model.ne.6)then
3648         call getairmol(latarg,matarg)
3649         if(ish.ge.2)write(ifch,*)'Air Target, select (Z,A) :'
3650      &                           ,latarg,matarg
3651       endif
3652
3653       if(iappl.eq.4)then
3654       call amicro
3655       if(ish.ge.2)call alist('list before int/decays&',1,nptl)
3656       nevt=1
3657       nbdky=nptl
3658       call bjinta(ier)
3659       if(ier.eq.1)stop'error in bjinta'
3660       if(ish.ge.2)call alist('list after int/decays&',1,nptl)
3661       goto1000
3662       endif
3663
3664       ntry=0
3665       if(nin.le.1)bimevt=-1
3666       ntevt0=ntevt
3667  3    ntevt=ntevt0
3668  2    ntevt=ntevt+1
3669       iret=0
3670       ntry=ntry+1
3671       if(iappl.eq.1.or.iappl.eq.2)naevt=naevt+1
3672  5    nevt=0
3673       if(nrevt.eq.0)nptly=nptl
3674       call cleanup
3675       nptl=0
3676
3677       minfra=mxptl
3678       maxfra=0
3679       if(iappl.eq.1.or.iappl.eq.2)then !---hadron---geometry---
3680       if(ntry.lt.10000.and.engy.ge.egymin)then    !if no inel scattering -> nothing !
3681         if(model.eq.1)then
3682           call emsaaa(iret)
3683         else
3684           call emsaaaModel(model,iret)
3685         endif
3686         if(iret.gt.0)goto2
3687       else
3688         ntevt=ntevt0+1
3689         iret=0
3690         nevt=1
3691         call conre     !define projectile and target (elastic scattering)
3692         call conwr     !when the MC is suppose to produce something but failed
3693       endif
3694       if(iappl.eq.2)then
3695         nevt=1
3696         goto 1000
3697       endif
3698
3699       elseif(iappl.eq.3) then
3700         call bread
3701         if(ish.ge.2)call alist('list after reading&',1,nptl)
3702         goto 1000
3703
3704       elseif(iappl.eq.5) then !---kinky---
3705          nptl=nptly
3706          nevt=1
3707          do i=1,nptl
3708            istptl(i)=20
3709          enddo
3710
3711       elseif(iappl.eq.6)then !---ee---
3712
3713         call timann
3714         nevt=1
3715
3716       elseif(iappl.eq.7)then !---decay---
3717
3718         call conwr
3719         nevt=1
3720
3721       elseif(iappl.eq.8)then !---lepton---
3722
3723         call psadis(iret)
3724         if(iret.gt.0)goto5
3725         nevt=1
3726
3727       endif
3728
3729       if(nevt.eq.0)stop'************ should not be ***************'
3730
3731       if(ish.ge.2)call alist('list before fragmentation&',1,nptl)
3732       nptlx=nptl+1
3733       if(iappl.ne.2.and.iappl.ne.7.and.nevt.eq.1.and.ifrade.ne.0)then
3734         call gakfra(iret)
3735         if(iret.gt.0)goto 3
3736         maxfra=nptl
3737         if(ish.ge.2.and.model.eq.1)
3738      &              call alist('list after fragmentation&',nptlx,nptl)
3739         if(irescl.eq.1)then
3740           call utghost(iret)
3741           if(iret.gt.0)goto 3
3742         endif
3743 c       nptlx=nptl+1
3744       endif
3745
3746       if(ispherio.eq.1.and.irescl.eq.1)then
3747         call utrsph(iret)
3748         if(iret.gt.0)goto 3
3749       endif
3750
3751
3752
3753       if(iappl.ne.2.and.nevt.eq.1)then
3754         nbdky=nptl
3755         call bjinta(ier)
3756         if(ier.eq.1)goto3
3757         if(iappl.eq.1.and.irescl.eq.1)then
3758           call utresc(iret)
3759           if(iret.gt.0)goto 3
3760         endif
3761         if(ish.ge.1)then
3762           !if(ish.ge.2.and.ifrade.ne.0)
3763           !&    call alist('list after int/decays&',1,nptl)
3764           if(iappl.eq.1)then
3765             numbar=0
3766             pp4=0.
3767             do j=1,nptl
3768               if(istptl(j).eq.0)then
3769              if(idptl(j).gt. 1000.and.idptl(j).lt. 10000)numbar=numbar+1
3770              if(idptl(j).lt.-1000.and.idptl(j).gt.-10000)numbar=numbar-1
3771              if((((idptl(j).eq.1120.or.idptl(j).eq.1220)
3772      *           .and.idproj.gt.1000).or.(iabs(idptl(j)).gt.100
3773      *           .and.idproj.lt.1000)).and.pptl(4,j)
3774      *           .gt.pp4.and.pptl(3,j).gt.0.)pp4=pptl(4,j)
3775               endif
3776             enddo
3777             pp4max=pp4max+pp4
3778             pp4ini=pp4ini+pptl(4,1)
3779             nvio=isign(matarg,idtarg)-numbar
3780             if(iabs(idproj).gt.1000)nvio=nvio+isign(maproj,idproj)
3781             if(ish.ge.2)write (ifch,*)'- Baryon number conservation : '
3782      &                  ,nvio,' -'
3783
3784           endif
3785         endif
3786       endif
3787
3788
3789       if((iappl.eq.1.or.iappl.eq.2).and.nevt.eq.0)then
3790         if(nin.le.1)bimevt=-1
3791         goto2
3792       endif
3793
3794       if(ifrade.ne.0.and.iappl.eq.2.and.
3795      $     idproj.eq.1120.and.idtarg.eq.1120)then
3796        numbar=0
3797        do j=1,nptl
3798         if(istptl(j).eq.0)then
3799          if(idptl(j).gt. 1000.and.idptl(j).lt. 10000)numbar=numbar+1
3800          if(idptl(j).lt.-1000.and.idptl(j).gt.-10000)numbar=numbar-1
3801         endif
3802        enddo
3803        nvio=maproj+matarg-numbar
3804        if(nvio.ne.0)then
3805         call alist('complete list&',1,nptl)
3806         write(6,'(//10x,a,i3//)')'ERROR: baryon number violation:',nvio
3807         write(6,'(10x,a//)')
3808      *        'a complete list has been printed into the check-file'
3809         stop
3810        endif
3811       endif
3812
3813
3814       ifirst=0
3815       if(nrevt+1.eq.1)ifirst=1
3816       if(jpsi.gt.0)then
3817         npjpsi=0
3818         do i=1,jpsi
3819           call jpsifo(npjpsi)
3820           call jpsian(ifirst)
3821         enddo
3822         if(ish.ge.1)call jtauan(0,0)
3823         if(nrevt+1.eq.nevent)call jpsihi
3824       endif
3825
3826       if(ixtau.eq.1)call xtauev(1)
3827
3828 1000  continue
3829       if(iabs(nin).eq.iabs(ninicon))nrevt=nrevt+1
3830       nglacc=nglacc+nglevt
3831       call timer(iutime)
3832       timefin=iutime(3)+float(iutime(4))/1000.
3833       call utprix('aepos',ish,ishini,4)
3834       return
3835       end
3836
3837 c-----------------------------------------------------------------------
3838       subroutine cleanup
3839 c-----------------------------------------------------------------------
3840       include 'epos.inc'
3841       do i=1,nptl
3842         do  k=1,5
3843           pptl(k,i)=0
3844         enddo
3845         iorptl(i)  =0
3846         jorptl(i)  =0
3847         idptl(i)   =0
3848         istptl(i)  =0
3849         tivptl(1,i)=0
3850         tivptl(2,i)=0
3851         ifrptl(1,i)=0
3852         ifrptl(2,i)=0
3853         ityptl(i)  =0
3854         iaaptl(i)  =0
3855         radptl(i)  =0
3856         dezptl(i)  =0
3857         itsptl(i)  =0
3858         rinptl(i)  =-9999
3859         do  k=1,4
3860           xorptl(k,i)=0
3861           ibptl(k,i) =0
3862         enddo
3863       enddo
3864       end
3865
3866 c-----------------------------------------------------------------------
3867       subroutine emsaaa(iret)
3868 c-----------------------------------------------------------------------
3869 c  basic EMS routine to determine Pomeron configuration
3870 c-----------------------------------------------------------------------
3871
3872       include 'epos.inc'
3873       common/col3/ncol,kolpt
3874
3875       call utpri('emsaaa',ish,ishini,4)
3876       if(ish.ge.3)call alist('Determine Pomeron Configuration&',0,0)
3877
3878       iret=0
3879
3880 1     nptl=0
3881       call conaa(iret)
3882       if(iret.gt.0)goto 1001   !no interaction
3883       if(iappl.eq.2.and.ixgeometry.eq.1)call xGeometry(1)
3884       if(iappl.eq.2)goto1000
3885       call conre
3886       call conwr
3887       call GfunParK(iret)
3888       if(iret.gt.0)goto 1    !redo
3889       if(ionudi.eq.0
3890      &  .and.(maproj.ne.1.or.matarg.ne.1).and.nglevt.eq.0)goto 1001
3891       call integom1(iret)
3892       if(iret.gt.0)goto 1    !redo
3893       call emsaa(iret)
3894       if(iret.gt.0)goto 1    !redo
3895       if(ncol.eq.0)goto 1001 !no interaction
3896
3897 1000  call utprix('emsaaa',ish,ishini,4)
3898       return
3899
3900 1001  iret=1
3901       goto 1000
3902
3903       end
3904
3905
3906 c----------------------------------------------------------------------
3907       subroutine alist(text,n1,n2)
3908 c----------------------------------------------------------------------
3909 c    ior  jor  i  ifr1  ifr2     id  ist  ity      pt  m  y
3910 c----------------------------------------------------------------------
3911 c       ist                                     ity
3912 c                                  light cluster ........ 19
3913 c   ptl ...  0 1                   soft pom ............. 20-23   25(reggeon)
3914 c   clu ... 10 11                  hard pom low mass .... 30
3915 c   ptn ... 20 21                  proj remnant ......... 40-49
3916 c   str ... 29                     targ remnant ......... 50-59
3917 c   pom ... 30 31 32(virtual)      cluster .............. 60
3918 c   rem ... 40 41                  direct photon ........ 71,72
3919 c----------------------------------------------------------------------
3920       include 'epos.inc'
3921       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
3922      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
3923       parameter(itext=40)
3924       character  text*40
3925       dimension pp(5)
3926       if(n1.gt.n2)return
3927       imax=itext+1
3928       do i=itext,1,-1
3929       if(text(i:i).eq.'&')imax=i
3930       enddo
3931       if(imax.gt.1)then
3932       write(ifch,'(/1x,89a1/1x,a,a,a,90a1)')
3933      *('#',k=1,89),'############  ',text(1:imax-1),'  '
3934      *,('#',k=1,74-imax)
3935       write(ifch,'(1x,89a1/)')('#',k=1,89)
3936       endif
3937       if(n1.eq.0.and.n2.eq.0)return
3938       if(imax.gt.1)then
3939       write(ifch,'(1x,a,a/1x,89a1)')
3940      *'   ior   jor     i  ifr1  ifr2       id ist ity',
3941      *'        pt         m         E         y'
3942      *,('-',k=1,89)
3943       endif
3944
3945       do j=1,5
3946         pp(j)=0.
3947       enddo
3948       nqu=0
3949       nqd=0
3950       nqs=0
3951       do i=n1,n2
3952         ptptl=pptl(1,i)**2.+pptl(2,i)**2.
3953         if(ptptl.le.0.d0)then
3954           ptptl=0
3955         else
3956           ptptl=sqrt(ptptl)
3957         endif
3958         amtptl=pptl(1,i)**2.+pptl(2,i)**2.+pptl(5,i)**2.
3959         if(amtptl.le.0.d0)then
3960           if(idptl(i).lt.10000)then
3961             call idmass(idptl(i),amtptl)
3962           endif
3963           amtptl=sqrt(amtptl*amtptl+pptl(1,i)**2.+pptl(2,i)**2.)
3964         else
3965           amtptl=sqrt(amtptl)
3966         endif
3967         rap=0.
3968         if(amtptl.gt.0.d0.and.pptl(4,i).gt.0.)
3969      &  rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amtptl)
3970          write(ifch,'(1x,i6,i6,i6,i6,i6,i10,2i3,2x,4(e9.3,1x),$)')
3971      &          iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
3972      &    ,idptl(i),istptl(i),ityptl(i),ptptl,pptl(5,i),pptl(4,i),rap
3973         write(ifch,*)' '
3974 c        if(istptl(i).ne.12)write(ifch,*)' '
3975 c        if(istptl(i).eq.12)write(ifch,'(1x,3(e9.3,1x))')
3976 c     &           sptl(i),sqrt(uptl(i)-xorptl(1,i)**2)
3977 c     &            ,sqrt(optl(i)-xorptl(2,i)**2)
3978         if(mod(istptl(i),10).eq.0.and.n1.eq.1.and.n2.eq.nptl)then
3979           do j=1,4
3980             pp(j)=pp(j)+pptl(j,i)
3981           enddo
3982           if(istptl(i).ne.40.and.istptl(i).ne.30)then
3983             call idqufl(i,idptl(i),ifl1,ifl2,ifl3)
3984             nqu=nqu+ifl1
3985             nqd=nqd+ifl2
3986             nqs=nqs+ifl3
3987           endif
3988         endif
3989       enddo
3990       end
3991
3992 c----------------------------------------------------------------------
3993       subroutine blist(text,n1,n2)
3994 c----------------------------------------------------------------------
3995       include 'epos.inc'
3996       parameter(itext=40)
3997       character  text*40
3998       dimension pp(5)
3999       if(n1.gt.n2)return
4000       imax=itext+1
4001       do i=itext,1,-1
4002       if(text(i:i).eq.'&')imax=i
4003       enddo
4004       if(imax.gt.1)then
4005       write(ifch,'(/1x,89a1/1x,a,a,a,90a1)')
4006      *('#',k=1,89),'#############  ',text(1:imax-1),'  '
4007      *,('#',k=1,74-imax)
4008       write(ifch,'(1x,89a1/)')('#',k=1,89)
4009       endif
4010       if(n1.eq.0.and.n2.eq.0)return
4011       if(imax.gt.1)then
4012       write(ifch,'(1x,a,a,a/1x,90a1)')
4013      *'   ior   jor     i  ifr1   ifr2      id ist ity',
4014      *'        pt      mass    energy','       rap'
4015      *,('-',k=1,90)
4016       endif
4017
4018       do j=1,5
4019         pp(j)=0.
4020       enddo
4021       nqu=0
4022       nqd=0
4023       nqs=0
4024       do i=n1,n2
4025         amtptl=pptl(1,i)**2.+pptl(2,i)**2.+pptl(5,i)**2.
4026         if(amtptl.le.0.d0)then
4027           if(idptl(i).lt.10000)then
4028             call idmass(idptl(i),amtptl)
4029           endif
4030           amtptl=sqrt(amtptl*amtptl+pptl(1,i)**2.+pptl(2,i)**2.)
4031         else
4032           amtptl=sqrt(amtptl)
4033         endif
4034         pt=pptl(1,i)**2.+pptl(2,i)**2.
4035         if(pt.gt.0.)pt=sqrt(pt)
4036         rap=0.
4037         if(amtptl.gt.0.d0.and.pptl(4,i).gt.0.)
4038      &  rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amtptl)
4039         write(ifch,125)iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
4040      &       ,idptl(i),istptl(i),ityptl(i)
4041      &       ,pt,pptl(5,i),pptl(4,i),rap
4042  125  format (1x,i6,i6,i6,i6,i6,i10,2i3,2x,5(e9.3,1x)
4043      *     ,f9.2,4x,5(e8.2,1x))
4044         if(mod(istptl(i),10).eq.0.and.n1.eq.1.and.n2.eq.nptl)then
4045           do j=1,4
4046             pp(j)=pp(j)+pptl(j,i)
4047           enddo
4048           if(istptl(i).ne.40.and.istptl(i).ne.30)then
4049             call idqufl(i,idptl(i),ifl1,ifl2,ifl3)
4050             nqu=nqu+ifl1
4051             nqd=nqd+ifl2
4052             nqs=nqs+ifl3
4053           endif
4054         endif
4055       enddo
4056       end
4057
4058 c----------------------------------------------------------------------
4059       subroutine alistf(text)
4060 c----------------------------------------------------------------------
4061       include 'epos.inc'
4062       parameter(itext=40)
4063       character  text*40
4064       dimension pp(5),erest(5),errp(4)
4065       n1=1
4066       n2=nptl
4067       imax=itext+1
4068       do i=itext,1,-1
4069       if(text(i:i).eq.'&')imax=i
4070       enddo
4071       if(imax.gt.1)then
4072       write(ifch,'(/1x,89a1/1x,a,a,a,90a1)')
4073      *('#',k=1,89),'#############  ',text(1:imax-1),'  '
4074      *,('#',k=1,74-imax)
4075       write(ifch,'(1x,75a1/)')('#',k=1,89)
4076       endif
4077       if(imax.gt.1)then
4078       write(ifch,'(1x,a,a,a/1x,124a1)')
4079      *'   ior   jor        i     ifr1   ifr2         id ist ity',
4080      *'            px         py         pz         p0       mass',
4081      *'       rap'
4082      *,('-',k=1,124)
4083       endif
4084
4085       do j=1,4
4086         pp(j)=0.
4087         errp(j)=0.
4088       enddo
4089       pp(5)=0.
4090       do i=n1,n2
4091         if(istptl(i).eq.0)then
4092         amtptl=pptl(1,i)**2.+pptl(2,i)**2.+pptl(5,i)**2.
4093         if(amtptl.le.0.d0)then
4094           if(idptl(i).lt.10000)then
4095             call idmass(idptl(i),amtptl)
4096           endif
4097           amtptl=sqrt(amtptl*amtptl+pptl(1,i)**2.+pptl(2,i)**2.)
4098         else
4099           amtptl=sqrt(amtptl)
4100         endif
4101         rap=0.
4102         if(amtptl.gt.0.d0.and.pptl(4,i).gt.0.)
4103      &  rap=sign(1.,pptl(3,i))*alog((pptl(4,i)+abs(pptl(3,i)))/amtptl)
4104         write(ifch,125)iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
4105      &       ,idptl(i),istptl(i),ityptl(i),(pptl(j,i),j=1,5),rap
4106 c     &,(xorptl(j,i),j=1,4)
4107         do j=1,4
4108           pp(j)=pp(j)+pptl(j,i)
4109         enddo
4110         endif
4111       enddo
4112  125  format (1x,i6,i6,3x,i6,3x,i6,i6,i12,2i4,4x,5(e10.4,1x)
4113      *     ,f9.2,4x,4(e8.2,1x))
4114  126  format (51x,5(e10.4,1x))
4115  128  format (51x,65('-'))
4116       pp(5)=(pp(4)-pp(3))*(pp(4)+pp(3))-pp(2)**2-pp(1)**2
4117       if(pp(5).gt.0.)then
4118         pp(5)=sqrt(pp(5))
4119       else
4120         pp(5)=0.
4121       endif
4122       write (ifch,128)
4123       write (ifch,126) (pp(i),i=1,5)
4124       erest(1)=0.
4125       erest(2)=0.
4126       erest(3)=maproj*pptl(3,1)+matarg*pptl(3,maproj+1)
4127       erest(4)=maproj*pptl(4,1)
4128      &        +matarg*pptl(4,maproj+1)
4129       erest(5)=amproj
4130       write (ifch,129)  (erest(j),j=1,5)
4131  129  format (50x,'(',5(e10.4,1x),')')
4132       do j=1,4
4133       if(abs(pp(j)).gt.0.d0)errp(j)=100.*(pp(j)-erest(j))/pp(j)
4134       enddo
4135       write (ifch,130)  (errp(j),j=1,4)
4136  130  format (50x,'(',3x,4(f7.2,4x),2x,'err(%))')
4137       end
4138
4139 c----------------------------------------------------------------------
4140       subroutine alist2(text,n1,n2,n3,n4)
4141 c----------------------------------------------------------------------
4142       include 'epos.inc'
4143       parameter(itext=40)
4144       character  text*40
4145       if(n1.gt.n2)return
4146       imax=itext+1
4147       do i=itext,1,-1
4148       if(text(i:i).eq.'&')imax=i
4149       enddo
4150       write(ifch,'(1x,a,a,a)')
4151      *'--------------- ',text(1:imax-1),' ---------------  '
4152       do i=n1,n2
4153       write(ifch,125)iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
4154      &,idptl(i),istptl(i),ityptl(i),(pptl(j,i),j=1,5)
4155 c     &,(xorptl(j,i),j=1,4)
4156       enddo
4157       write(ifch,'(1x,a)')'----->'
4158       do i=n3,n4
4159       write(ifch,125)iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
4160      &,idptl(i),istptl(i),ityptl(i),(pptl(j,i),j=1,5)
4161 c     &,(xorptl(j,i),j=1,4)
4162       enddo
4163  125  format (1x,i6,i6,3x,i6,3x,i6,i6,i12,2i4,4x,5(e8.2,1x))
4164 c     *,4x,4(e8.2,1x))
4165       end
4166
4167 c----------------------------------------------------------------------
4168       subroutine alistc(text,n1,n2)
4169 c----------------------------------------------------------------------
4170       include 'epos.inc'
4171       parameter(itext=40)
4172       character  text*40
4173       if(n1.gt.n2)return
4174       imax=itext+1
4175       do i=itext,1,-1
4176       if(text(i:i).eq.'&')imax=i
4177       enddo
4178       if(n1.ne.n2)write(ifch,'(1x,a,a,a)')
4179      *'--------------- ',text(1:imax-1),' ---------------  '
4180       do i=n1,n2
4181       write(ifch,130)iorptl(i),jorptl(i),i,ifrptl(1,i),ifrptl(2,i)
4182      &,idptl(i),istptl(i),ityptl(i),(pptl(j,i),j=1,5)
4183      &,(xorptl(j,i),j=1,4),tivptl(1,i),tivptl(2,i)
4184       enddo
4185  130  format (1x,i6,i6,3x,i6,3x,i6,i6,i12,2i4,4x,5(e8.2,1x)
4186      *,4x,6(e8.2,1x))
4187       end
4188
4189 c-----------------------------------------------------------------------
4190       subroutine sigmaint(g0,gz,sigdo)
4191 c-----------------------------------------------------------------------
4192 c hadron-hadron cross sections integration
4193 c-----------------------------------------------------------------------
4194       common /ar3/  x1(7),a1(7)
4195       include 'epos.inc'
4196       include 'epos.incsem'
4197       double precision PhiExact,vvv11,vvv12,vvv21,om1intbi!,PomNbri
4198      *,vvv22,ww01,ww02,ww11,ww12,ww21,ww22,gz(0:3)
4199      *,vvv11e,vvv12e,vvv21e,vvv22e,PhiExpo
4200
4201
4202 c       rs=r2had(iclpro)+r2had(icltar)+slopom*log(engy**2)
4203         rs=r2had(iclpro)+r2had(icltar)+max(slopom,slopoms)*log(engy**2)
4204      &     +gwidth*(r2had(iclpro)+r2had(icltar))
4205      &     +bmxdif(iclpro,icltar)/4./0.0389
4206         rpom=4.*.0389*rs
4207         e1=exp(-1.)
4208 c        cpt=chad(iclpro)*chad(icltar)
4209
4210         gz(0)=0.d0
4211         gz(1)=0.d0
4212         gz(2)=0.d0
4213         gz(3)=0.d0
4214
4215         sigdo=0.
4216         do i1=1,7
4217         do m=1,2
4218
4219           z=.5+x1(i1)*(m-1.5)
4220           zv1=exp(-z)
4221           zv2=(e1*z)
4222           b1=sqrt(-rpom*log(zv1))
4223           b2=sqrt(-rpom*log(zv2))
4224
4225           vvv11=max(0.d0,PhiExpo(.5,1.d0,1.d0,engy**2,b1))
4226           vvv12=max(0.d0,PhiExpo(.5,1.d0,1.d0,engy**2,b2))
4227           vvv21=max(0.d0,PhiExpo(1.,1.d0,1.d0,engy**2,b1))
4228           vvv22=max(0.d0,PhiExpo(1.,1.d0,1.d0,engy**2,b2))
4229 c          vvv11=sqrt(vvv21)
4230 c          vvv12=sqrt(vvv21)
4231           if(isetcs.eq.0)then
4232             vvv11e=max(0.d0,PhiExact(.5,1.d0,1.d0,engy**2,b1))
4233             vvv12e=max(0.d0,PhiExact(.5,1.d0,1.d0,engy**2,b2))
4234             vvv21e=max(0.d0,PhiExact(1.,1.d0,1.d0,engy**2,b1))
4235             vvv22e=max(0.d0,PhiExact(1.,1.d0,1.d0,engy**2,b2))
4236             vvv11=0.5d0*(vvv11+min(vvv11,vvv11e))  !to be as close as possible
4237             vvv12=0.5d0*(vvv12+min(vvv12,vvv12e))  !than the value with
4238             vvv21=0.5d0*(vvv21+min(vvv21,vvv21e))  !isetcs > 0
4239             vvv22=0.5d0*(vvv22+min(vvv22,vvv22e))
4240           endif
4241           ww11=1.d0-vvv11
4242           ww12=1.d0-vvv12
4243           ww21=1.d0-vvv21
4244           ww22=1.d0-vvv22
4245           ww01=vvv21-2d0*vvv11+1d0
4246           ww02=vvv22-2d0*vvv12+1d0
4247
4248           gz(0)=gz(0)+a1(i1)*(ww01+ww02/z)
4249           gz(1)=gz(1)+a1(i1)*(ww11+ww12/z)
4250           gz(2)=gz(2)+a1(i1)*(ww21+ww22/z)
4251           gz(3)=gz(3)+a1(i1)*(ww11*z+ww12/z*(1.-log(z)))
4252
4253           phi1=vvv21
4254           phi2=vvv22
4255           phi1x=real(min(50d0,exp(om1intbi(b1,2)/dble(r2hads(iclpro)
4256      &                                               +r2hads(icltar)))))
4257           phi2x=real(min(50d0,exp(om1intbi(b2,2)/dble(r2hads(iclpro)
4258      &                                               +r2hads(icltar)))))
4259           sigdo=sigdo+a1(i1)*(phi1*(phi1x-1.)+phi2*(phi2x-1.)/z)
4260
4261         enddo
4262         enddo
4263         g0=pi*rpom*10./2.                 !common factor (pi*rpom because of b->z, 10 to have mbarn and /2. because z is from -1 to 1 but we need 0 to 1.
4264
4265       return
4266       end
4267 c-----------------------------------------------------------------------
4268       subroutine xsigma
4269 c-----------------------------------------------------------------------
4270 c hadron-hadron and hadron-nucleus cross sections calculation
4271 c b - impact parameter squared (in case of hadron-nucleus interaction);
4272 c-----------------------------------------------------------------------
4273       include 'epos.inc'
4274       include 'epos.incsem'
4275       double precision gz(0:3),gzp(0:3),GZ0(2)
4276 c Model 2 Common
4277       COMMON /Q_AREA1/  IA(2),ICZ,ICP
4278       COMMON /Q_AREA6/  PIQGS,BM,AM
4279       COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
4280       COMMON /Q_AREA7/  RP1
4281       COMMON /Q_AREA16/ CC(5)
4282       double precision RP1,FP,RQ,CD,PIQGS,BM,AM,CC,GDP,GDT,GDD
4283
4284 c Model 3 Common
4285       COMMON/GSECTI/ AIEL(20),AIIN(20),AIFI(20),AICA(20),ALAM,K0FLAG
4286       INTEGER K0FLAG
4287       DOUBLE PRECISION AIEL,AIIN,AIFI,AICA,ALAM
4288       PARAMETER (KKK=3)
4289       COMMON /CGCOMP/ ACOMP,ZCOMP,WCOMP,KK
4290       REAL           ACOMP(KKK),ZCOMP(KKK),WCOMP(KKK)
4291
4292 C...Total cross sections in Pythia
4293 c      double precision SIGT
4294 c      COMMON/PYINT7/SIGT(0:6,0:6,0:5)
4295
4296 c Model 5 Common
4297       COMMON/HIPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50)
4298
4299 c theoretical cross sections
4300       sigcut=0.
4301       sigtot=0.
4302       sigela=0.
4303       sigine=0.
4304       sigtotold=0.
4305       sigtotf=0.
4306       sigelaold=0.
4307       sigelaf=0.
4308       sigineold=0.
4309       sigineaa=0.
4310       sigdif=0.
4311       sloela=0.
4312       sigsd=0.
4313 c simulated cross sections
4314       sigineex=0.          !calclated in ems if isigma=1
4315       sigdifex=0.
4316       sigsdex=0.
4317
4318
4319       call utpri('xsigma',ish,ishini,4)
4320
4321       if(model.eq.1)then                        !epos
4322
4323         if(icltar.ne.2)stop'xsigma: only icltar=2 allowed.'
4324
4325         call sigmaint(g0,gz,sigdifold)
4326
4327         sigelaold=g0*gz(0)               !elastic cross section
4328         rs=g0*0.4091    !=g0/pi/10.*2./4./.0389
4329         if(gz(1).gt.0d0)sloela=2.*rs*gz(3)/gz(1)
4330
4331         sigineold=g0*gz(2)               !inelastic pomerons cross-section
4332         sigtotold=2.*g0*gz(1)                  !tot cross-section
4333         sigdifold=sigdifold * g0 !xs in mb
4334         sigcut=sigineold-sigdifold             !cut cross section
4335         x=engy
4336         sigtotf=21*x**0.17+44.*x**(-0.8) !fit to data
4337         sigelaf=30.*(x-1)**(-3)+17*x**(-0.47)+0.3*alog(x)**2 !fit to data
4338         sigtotfp=sigtotf
4339         sigelafp=sigelaf
4340         if(iclpro.eq.1)then        !pi+p
4341           sigtotf=15.*x**0.15+55.*x**(-1.5)                     !fit to data
4342           sigelaf=20.*(x-1)**(-3)+6.*x**(-0.4)+0.2*alog(x)**2. !fit to data
4343         elseif(iclpro.eq.3)then    !K+p
4344           sigtotf=12.5*x**0.15+35.*x**(-1.5)                     !fit to data
4345           sigelaf=15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data
4346         elseif(iclpro.eq.4)then    !D+p
4347           sigtotf=0.!12.5*x**0.15+35.*x**(-1.5)                     !fit to data
4348           sigelaf=0.!15.*(x-1)**(-3)+3.*x**(-0.4)+0.2*alog(x)**2 !fit to data
4349         endif
4350         if(engy.lt.20.)then
4351           sigcoul=max(0.,sigtotf-sigtotold)
4352         else
4353           sigcoul=0.
4354         endif
4355
4356
4357         sigdif=sigdifold
4358         sigdelaf=max(0.,(sigelaf+sigineold-sigtotf))
4359 c        sigdelaf=max(0.,(sigelaf-sigelaold))
4360
4361         if(rexdifi(iclpro).lt.0..or.rexdifi(icltar).lt.0.)then
4362           sigdela=min(sigdelaf,sigdif)
4363
4364 c calculate rexdif for proton first (always needed)
4365           if(rexdifi(icltar).lt.0.)then
4366
4367             if(engy.lt.min(30.,-log(-rexdifi(icltar))/0.015))then !use fit of sigela to get rexdif
4368               if(iclpro+icltar.ne.4)then !pi or K - p, calculate sigdif for pp
4369                 iclprosave=iclpro
4370                 iclpro=2
4371                 call sigmaint(g0p,gzp,sigdifp)
4372                 siginep=g0p*gzp(2)
4373                 sigdifp=sigdifp * g0p
4374                 sigdelafp=max(0.,(sigelafp+siginep-sigtotfp))
4375                 sigdelap=min(sigdelafp,sigdifp)
4376                 iclpro=iclprosave
4377               else
4378                 sigdifp=sigdif
4379                 sigdelafp=sigdelaf
4380                 sigdelap=sigdela
4381               endif
4382               if(sigdifp.gt.0.)then
4383                 ratioSig=sigdelap/sigdifp
4384                 rexdif(icltar)=1.-sqrt(ratioSig)
4385                 if(rexdif(icltar).ge.exp(-0.015*engy))
4386      &               rexdif(icltar)=exp(-0.015*engy)
4387                 rexdif(icltar)=max(rexdif(icltar),abs(rexdifi(icltar)))
4388               else
4389                 rexdif(icltar)=1.
4390               endif
4391             else
4392 c         rexdif(icltar)=max(exp(-1.7/engy**0.3),abs(rexdifi(icltar)))  !strong reduction
4393 c        rexdif(icltar)=max(exp(-0.33/engy**0.066),abs(rexdifi(icltar))) !moderate one (constant sig NSD)
4394               rexdif(icltar)=abs(rexdifi(icltar))
4395             endif
4396           else
4397             rexdif(icltar)=rexdifi(icltar)
4398           endif
4399
4400          if(iclpro.ne.2)then  !pi or K rexdif knowing p rexdif
4401           if(rexdifi(iclpro).lt.0.)then
4402             if(engy.lt.min(8.,-log(-rexdifi(iclpro))/0.015)
4403      &         .and.sigdif.gt.0.)then !use fit of sigela to get rexdif
4404               ratioSig=sigdela/sigdif
4405               if(abs(1.-rexdif(icltar)).gt.1.e-6)then
4406                 rexdif(iclpro)=1.-ratioSig/(1.-rexdif(icltar))
4407               else
4408                 rexdif(iclpro)=abs(rexdifi(iclpro))
4409               endif
4410               if(rexdif(iclpro).ge.exp(-0.015*engy))
4411      &             rexdif(iclpro)=exp(-0.015*engy)
4412             elseif(sigdif.le.0.)then
4413               rexdif(iclpro)=1.
4414             else
4415 c          rexdif(iclpro)=max(exp(-1.7/engy**0.3),abs(rexdifi(iclpro))) !strong reduction
4416 c         rexdif(iclpro)=max(exp(-0.33/engy**0.066),abs(rexdifi(iclpro)))  !moderate one (constant sig NSD)
4417               rexdif(iclpro)=abs(rexdifi(iclpro))
4418             endif
4419            else
4420              rexdif(iclpro)=abs(rexdifi(iclpro))
4421            endif
4422          endif
4423          sigdela=(1.-rexdif(iclpro))
4424      &           *(1.-rexdif(icltar))  *sigdif
4425
4426         else
4427           rexdif(iclpro)=rexdifi(iclpro)
4428           rexdif(icltar)=rexdifi(icltar)
4429           sigdela=(1.-rexdif(iclpro))
4430      &           *(1.-rexdif(icltar))  *sigdif
4431         endif
4432
4433
4434         if(rexndf.gt.0..and.iclpro.eq.2)then
4435           rexndi(iclpro)=rexndf*rexdif(iclpro)
4436         else
4437           rexndi(iclpro)=rexndii(iclpro)
4438         endif
4439         if(rexndf.gt.0..and.icltar.eq.2)then
4440           rexndi(icltar)=rexndf*rexdif(icltar)
4441         else
4442           rexndi(icltar)=rexndii(icltar)
4443         endif
4444         if(ish.ge.2)write(ifch,*)'Xsigma : rexdif/ndi=',rexdif(iclpro)
4445      &                                                 ,rexdif(icltar)
4446      &                                                 ,rexndi(iclpro)
4447      &                                                 ,rexndi(icltar)
4448
4449         sigsd=( (1.-rexdif(icltar))*rexdif(iclpro)
4450      &         +(1.-rexdif(iclpro))*rexdif(icltar) )*sigdif
4451 c        if(engy.lt.10.)sigela=max(sigelaf,sigela)
4452         sigela=sigelaold+sigdela+sigcoul
4453         sigine=sigineold-sigdela
4454         sigtot=sigine+sigela
4455         sigineaa=eposcrse(ekin,maproj,matarg,idtarg)
4456
4457 c        write(ifmt,*)'Rexdif',rexdif(iclpro),rexdif(icltar)
4458
4459       elseif(model.eq.2)then
4460
4461         g0=real(PIQGS*RP1/CD(ICZ)*AM**2*10.D0)
4462         CALL m2XXFZ(0.D0,GZ0)
4463         gz(1)=GZ0(1)
4464         gz(2)=GZ0(2)
4465         gz(3)=0d0
4466         sigcut=g0*gz(2)/2.               !cut pomerons cross-section
4467         sigtot=g0*gz(1)                  !tot cross-section
4468         gz(0)=sigtot-sigcut
4469         sigela=gz(0)*CC(ICZ)*CC(2)       !elastic cross section
4470 c GDP - projectile diffraction cross section
4471         GDP=(1.D0-CC(ICZ))*CC(2)*gz(0)
4472 c GDT - target diffraction cross section
4473         GDT=(1.D0-CC(2))*CC(ICZ)*gz(0)
4474 c  GDD - double diffractive cross section
4475         GDD=(1.D0-CC(ICZ))*(1.D0-CC(2))*gz(0)
4476         sigsd=GDT+GDD
4477         sigdif=sigsd+GDP
4478         sigine=sigcut+sigdif
4479         if(gz(1).gt.0.)sloela=2.*rs*gz(3)/gz(1)
4480         sigdifold=sigtot-sigcut-sigela       !diffractive cross section
4481         sigineaa=crseModel(2,ekin,maproj,matarg,idtarg)
4482
4483       elseif(model.eq.3)then
4484
4485         call m3SIGMA(ekin,idproj,1120,1,1,sigi,sige)
4486         sigine=sigi
4487         sigela=sige
4488         sigcut=sigine
4489         sigtot=sigine+sigela
4490         sigdif=sigtot-sigcut-sigela       !diffractive cross section
4491         call m3SIGMA(ekin,idproj,idtarg,latarg,matarg,sigi,sige)
4492         sigineaa=sigi
4493
4494 c$$$      elseif(model.eq.4)then
4495 c$$$
4496 c$$$        sigsd=real(SIGT(0,0,2)+SIGT(0,0,3))
4497 c$$$        sigela=real(SIGT(0,0,1))
4498 c$$$        sigcut=real(SIGT(0,0,5))
4499 c$$$        sigtot=real(SIGT(0,0,0))
4500 c$$$        sigine=sigtot-sigela
4501 c$$$        sigdif=sigtot-sigcut-sigela       !diffractive cross section
4502
4503       elseif(model.eq.5)then
4504
4505         sigsd=0.
4506         sigdif=0.
4507         sigcut=0.
4508         sigtot=HINT1(13)
4509         sigine=HINT1(12)
4510         sigela=sigtot-sigine
4511         sigineaa=crseModel(5,ekin,maproj,matarg,idtarg)
4512
4513       elseif(model.eq.6)then                  !for Sibyll
4514
4515         call m6SIGMA(iclpro,engy,stot,sela,sine,sdifr,slela,Rho)
4516         sigtot=stot
4517         sigela=sela
4518         sigine=sine
4519         sigdif=sdifr
4520         sloela=slela
4521         sigcut=sigtot-sigdif-sigela     ! cut cross section
4522         sigsd=sigdif/2.
4523         sigineaa=crseModel(6,ekin,maproj,matarg,idtarg)
4524
4525       elseif(model.eq.7)then                  !for QGSJET-II
4526
4527         call m7SIGMA(stot,scut,sine,slela)
4528         sigtot=stot
4529         sigcut=scut
4530         sigine=sine
4531         sloela=slela
4532         sigela=sigtot-sigine     ! elastic cross section
4533         sigdifr=sigine-sigcut
4534         sigsd=sigdifr
4535         sigineaa=crseModel(7,ekin,maproj,matarg,idtarg)
4536
4537       endif
4538
4539              if(isigma.eq.1)then  !===============!
4540
4541       if(ish.ge.1)
4542      *write (ifmt,225)engy,ekin,sigtot,sigtotf,sigtotold
4543      *,sigine,sigtotf-sigelaf,sigineold
4544      *,sigela,sigelaf,sigelaold,sigcut,sloela,sigdif,sigsd
4545      *,sigineaa
4546       if(ish.ge.1)
4547      *write (ifch,225)engy,ekin,sigtot,sigtotf,sigtotold
4548      *,sigine,sigtotf-sigelaf,sigineold
4549      *,sigela,sigelaf,sigelaold,sigcut,sloela,sigdif,sigsd
4550      *,sigineaa
4551
4552       if(.not.(maproj.eq.1.and.matarg.eq.1))then  !(from tabulation) for pA/AA
4553
4554         sigtotaa=0.
4555         sigineaa=0.
4556         sigelaaa=0.
4557         if(model.eq.1)sigineaa=eposcrse(ekin,maproj,matarg,idtarg)
4558         if(model.eq.2)sigineaa=crseModel(2,ekin,maproj,matarg,idtarg)
4559         if(model.eq.3)then
4560           sigineaa=gheincs
4561           sigelaaa=0.
4562           do k=1,kk
4563             sigelaaa=sigelaaa+AIEL(k)
4564           enddo
4565           sigtotaa=sigineaa+sigelaaa
4566         endif
4567         if(model.eq.4)sigineaa=crseModel(4,ekin,maproj,matarg,idtarg)
4568         if(model.eq.5)sigineaa=crseModel(5,ekin,maproj,matarg,idtarg)
4569         if(model.eq.6)sigineaa=crseModel(6,ekin,maproj,matarg,idtarg)
4570         if(model.eq.7)sigineaa=crseModel(7,ekin,maproj,matarg,idtarg)
4571         if(ish.ge.1)
4572      &  write (ifmt,226)sigtotaa,sigineaa,sigelaaa
4573         if(ish.ge.4)write (ifch,226)sigtotaa,sigineaa,sigelaaa
4574
4575              endif  !================!
4576
4577
4578       endif
4579
4580  100  continue
4581
4582
4583 225   format(' hadron-proton cross sections for ',f10.2,' GeV',
4584      *'  (ekin:',f18.3,' GeV)'/
4585      *4x,'total cross section:           ',f8.2,3x,f8.2,3x,f8.2/
4586      *4x,'inelastic cross section:       ',f8.2,3x,f8.2,3x,f8.2/
4587      *4x,'elastic cross section:         ',f8.2,3x,f8.2,3x,f8.2/
4588      *4x,'cut cross section:             ',f8.2/
4589      *4x,'elastic slope parameter:       ',f8.2/
4590      *4x,'diffr. cross section:          ',f8.2,14x,f8.2/
4591      *4x,'inelastic (tab) cross section: ',f8.2)
4592  226  format(' hadron/nucleus-hadron/nucleus cross sections'/
4593      *4x,'total pA/AA cross section:     ',f8.2/
4594      *4x,'inelastic pA/AA cross section: ',f8.2/
4595      *4x,'elastic pA/AA cross section:   ',f8.2)
4596
4597       call utprix('xsigma',ish,ishini,4)
4598
4599       return
4600       end
4601