]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EPOS/epos167/epos-bas-168.f
Merge branch 'master' of http://git.cern.ch/pub/AliRoot
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-bas-168.f
CommitLineData
9ef1c2d9 1
2c-----------------------------------------------------------------------
3 subroutine aaset(iop)
4c-----------------------------------------------------------------------
5c sets parameters and options, initializes counters ...
6c-----------------------------------------------------------------------
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
37c version
38
39 iversn=int(1.67*100) !version number
40 iverso=int(1.67*100) !last official version
41
42c application
43
44 iappl=1 !choice for application (0,1,2,3,4)
45
46c model
47
48 model=1
49 iquasiel=1 !allow (1) or not (0) quasi-elastic event in model 3
50
51c 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
85c 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
96c 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
116c 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
128c 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
157c 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
171c 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
189c 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
206c 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
236c 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
267c 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
325c 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
334c ... 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
359c 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
381c 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
388c 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
399c 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
417c 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)
425c 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
433c 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
452c 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
469c coupling with urqmd
470
471 iurqmd=0 ! call eposurqmd (1) or not
472
473c 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
484c 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
504c 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
521c 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
532c 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
586c 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
593c tests
594
595 iotst1=0 !test
596 iotst2=0 !test
597 iotst3=0 !test
598 iotst4=0 !test
599
600c 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
609c 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
626c 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
652c 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
682c 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
691c 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
707c 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
7286 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
806c 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
827c-----------------------------------------------------------------------
828 subroutine ustore
829c-----------------------------------------------------------------------
830c writes the results of a simulation into the file with unit ifdt
831c contains a description of the stored variables.
832c modifiable by the user
833c-----------------------------------------------------------------------
834 include 'epos.inc'
835
836c 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
843c store event variables:
844
845 write(ifdt,*)nrevt,nptevt,bimevt,phievt,kolevt
846 *,pmxevt,egyevt,npjevt,ntgevt
847
848c nrevt.......... event number
849c nptevt ........ number of (stored!) particles per event
850c bimevt ........ absolute value of impact parameter
851c phievt ........ angle of impact parameter
852c kolevt ........ number of collisions
853c pmxevt ........ reference momentum
854c egyevt ........ pp cm energy (hadron) or string energy (lepton)
855c npjevt ........ number of primary projectile participants
856c ntgevt ........ number of primary target participants
857c npnevt ........ number of primary projectile neutron spectators
858c nppevt ........ number of primary projectile proton spectators
859c ntnevt ........ number of primary target neutron spectators
860c ntpevt ........ number of primary target proton spectators
861c jpnevt ........ number of absolute projectile neutron spectators
862c jppevt ........ number of absolute projectile proton spectators
863c jtnevt ........ number of absolute target neutron spectators
864c 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
870c 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
878c i ............. particle number
879c idptl(i) ...... particle id
880c pptl(1,i) ..... x-component of particle momentum
881c pptl(2,i) ..... y-component of particle momentum
882c pptl(3,i) ..... z-component of particle momentum
883c pptl(4,i) ..... particle energy
884c pptl(5,i) ..... particle mass
885c iorptl(i) ..... particle number of father (if .le. 0 : no father)
886c jorptl(i) ..... particle number of mother (if .le. 0 : no mother)
887c istptl(i) ..... generation flag: last gen. (0) or not (1)
888c ityptl(i) ..... particle type (string, remnant ...)
889c xorptl(1,i) ... x-component of formation point
890c xorptl(2,i) ... y-component of formation point
891c xorptl(3,i) ... z-component of formation point
892c xorptl(4,i) ... formation time
893c tivptl(1,i) ... formation time (always in the pp-cms!)
894c tivptl(2,i) ... destruction time (always in the pp-cms!)
895
896 endif
897 enddo
898
899 return
900 end
901
902
903c-----------------------------------------------------------------------
904 subroutine bstora
905c-----------------------------------------------------------------------
906c writes the results of a simulation into the file with unit ifdt
907c contains a description of the stored variables.
908c-----------------------------------------------------------------------
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
983c-----------------------------------------------------------------------
984 subroutine bstore
985c-----------------------------------------------------------------------
986c writes the results of a simulation into the file with unit ifdt
987c contains a description of the stored variables.
988c-----------------------------------------------------------------------
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
1084c-----------------------------------------------------------------------
1085 subroutine bread
1086c-----------------------------------------------------------------------
1087c writes the results of a simulation into the file with unit ifdt
1088c contains a description of the stored variables.
1089c-----------------------------------------------------------------------
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
1199c if(i.eq.15)read(line(k:),'(i6)') idiptl(n)
1200 if(l.eq.15) k=k+7
1201c 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
1205c if(l.eq.21)read(line(k:),'(I6)') ifrptl(1,n)
1206 if(l.eq.21) k=k+7
1207c 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
1217c-----------------------------------------------------------------------
1218 subroutine aafinal
1219c-----------------------------------------------------------------------
1220c * calculates xorptl(j,i), representing formation points.
1221c (xorptl(j,i),j=1,4) is the 4-vector representing the space-time of
1222c creation of particle i.
1223c-----------------------------------------------------------------------
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
1260c-----------------------------------------------------------------------
1261 subroutine afinal
1262c-----------------------------------------------------------------------
1263c does some final calculations, to be called before call aasto.
1264c * calculates nptlu, the maximum nptl for all events.
1265c * in case of mod(iframe,10) .ne. 1, these vectors are transformed
1266c (being originally in the "natural frame",
1267c * calculates numbers of spectators:
1268c npnevt (number of primary proj neutron spectators)
1269c nppevt (number of primary proj proton spectators)
1270c ntnevt (number of primary targ neutron spectators)
1271c ntpevt (number of primary targ proton spectators)
1272c jpnevt (number of absolute proj neutron spectators)
1273c jppevt (number of absolute proj proton spectators)
1274c jtnevt (number of absolute targ neutron spectators)
1275c jtpevt (number of absolute targ proton spectators)
1276c-----------------------------------------------------------------------
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
1337c if(isto.eq.1)stop
1338c$$$ 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
1415c-----------------------------------------------------------------------
1416 subroutine bfinal
1417c-----------------------------------------------------------------------
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
1452c-----------------------------------------------------------------------
1453 subroutine ainit
1454c-----------------------------------------------------------------------
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
1483ctp060829 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 !&')
1535c if((idtarg.ne.1120.and.(latarg.ne.-1.or.matarg.ne.1))
1536c & .or.matarg.le.0)
1537c & 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
1697c 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
1725c$$$ if(idproj.eq.1120)icp=2 !????????????? for what ?
1726c$$$ if(idproj.eq.-1120)icp=-2
1727c$$$ if(idproj.eq.120)icp=1
1728c$$$ if(idproj.eq.-120)icp=-1
1729c$$$ if(idproj.eq.130)icp=4
1730c$$$ 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.)
1777c 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
1784c 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
1798c additional initialization procedures
1799
1800 if(model.ne.1)then
1801 call IniEvtModel
1802 elseif(iappl.le.3)then
1803c 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
1816ccc call MakeFpartonTable
1817
1818c$$$ call testconex(1)
1819
1820 call utprix('ainit ',ish,ishini,4)
1821 return
1822 end
1823
1824c---------------------------------------------------------------------
1825 subroutine aread
1826c---------------------------------------------------------------------
1827c reads and interprets input commands
1828c---------------------------------------------------------------------
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
1949c$$$ 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
1961c$$$ if(line(i:j).eq.'xParOmega3')call xParOmega3
1962c$$$ 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
1966c$$$ if(line(i:j).eq.'xParSigma2')call xParSigma2
1967c$$$ if(line(i:j).eq.'xScrD')call xScrD
1968 if(line(i:j).eq.'xFitD1')call xFitD1
1969c$$$ if(line(i:j).eq.'xExaD2')call xExaD2
1970 if(line(i:j).eq.'xbExaD')call xbExaD
1971c$$$ if(line(i:j).eq.'xbExaD2')call xbExaD2
1972 if(line(i:j).eq.'xbnExaD')call xbnExaD
1973c$$$ 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
1976c$$$ 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
1980c$$$ if(line(i:j).eq.'xmParD2')call xmParD2
1981 if(line(i:j).eq.'xyParD')call xyParD
1982c$$$ 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
1992c 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.)
1998c 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)
2016c 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
2408ccc close(20+nopen)
2409ccc nopen=nopen-1
2410ccc 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
2478c 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)
2489c constants
2490 if(linex(ix:jx).eq.'ainfin')ainfin=sngl(val)
2491c 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)
2502c 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)
2535c 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)
2558c 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)
2589c 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)
2615c 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)
2677c 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
2724c 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)
2738c 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)
2745c 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)
2751c 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)
2782c urqmd
2783 if(linex(ix:jx).eq.'iurqmd')iurqmd=int(val)
2784c spherio
2785 if(linex(ix:jx).eq.'ispherio')ispherio=int(val)
2786c 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)
2792c 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)
2806c 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)
2820c 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)
2852c strangelets
2853 if(linex(ix:jx).eq.'iopenu')iopenu=nint(val)
2854 if(linex(ix:jx).eq.'themas')themas=sngl(val)
2855c 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)
2860c 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)
2867c 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)
2881c 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)
2918c 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
2990ctp290806 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
2996ctp290806 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)
3015c if(ne.eq.0.and.iprmpt.gt.0)
3016c * 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)
3022c if(ne.eq.0.and.iprmpt.gt.0)
3023c * 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)
3096c if(line(i:j).eq.'on' )iorsce=1
3097 if(line(i:j).eq.'on' )iorshh=1
3098c 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
3472c-----------------------------------------------------------------------
3473 subroutine aseed(modus)
3474c-----------------------------------------------------------------------
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
3497c-----------------------------------------------------------------------
3498 subroutine aseedi
3499c-----------------------------------------------------------------------
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
3513c$$$c-----------------------------------------------------------------------
3514c$$$ subroutine aseed(modus) !Flush ????
3515c$$$c-----------------------------------------------------------------------
3516c$$$
3517c$$$ include 'epos.inc'
3518c$$$ double precision seedf
3519c$$$ call utpri('aseed',ish,ishini,4)
3520c$$$
3521c$$$ call ranfgt(seedf)
3522c$$$ if(modus.eq.2)then
3523c$$$ write(ifmt,'(a,d26.15)')'seed:',seedf
3524c$$$ elseif(modus.eq.1)then
3525c$$$ if(mod(nrevt,modsho).eq.0)then
3526c$$$ write(ifmt,100)'nrevt:',nrevt,'seedf:',seedf
3527c$$$ call flush(ifmt)
3528c$$$ endif
3529c$$$ endif
3530c$$$ seedc=seedf
3531c$$$
3532c$$$ 100 format(a,i10,10x,a,d26.15)
3533c$$$ call utprix('aseed',ish,ishini,4)
3534c$$$ return
3535c$$$ end
3536c$$$
3537c-----------------------------------------------------------------------
3538 subroutine astati
3539c-----------------------------------------------------------------------
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
3569c if(model.eq.3)write(ifch,*)'quasi-elastic cross section:'
3570c &,anquasiel/float(ntevt)*a*10
3571 endif
3572
3573c$$$ 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
35841000 continue
3585 call utprix('astati',ish,ishini,1)
3586
3587 return
3588 end
3589
3590c-----------------------------------------------------------------------
3591 subroutine atitle
3592c-----------------------------------------------------------------------
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
3610c-----------------------------------------------------------------------
3611 subroutine avehep
3612c-----------------------------------------------------------------------
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
3623c-----------------------------------------------------------------------
3624 subroutine aepos(nin)
3625c-----------------------------------------------------------------------
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
3642c 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
3646c 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
3743c 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
38281000 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
3837c-----------------------------------------------------------------------
3838 subroutine cleanup
3839c-----------------------------------------------------------------------
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
3866c-----------------------------------------------------------------------
3867 subroutine emsaaa(iret)
3868c-----------------------------------------------------------------------
3869c basic EMS routine to determine Pomeron configuration
3870c-----------------------------------------------------------------------
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
38801 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
38971000 call utprix('emsaaa',ish,ishini,4)
3898 return
3899
39001001 iret=1
3901 goto 1000
3902
3903 end
3904
3905
3906c----------------------------------------------------------------------
3907 subroutine alist(text,n1,n2)
3908c----------------------------------------------------------------------
3909c ior jor i ifr1 ifr2 id ist ity pt m y
3910c----------------------------------------------------------------------
3911c ist ity
3912c light cluster ........ 19
3913c ptl ... 0 1 soft pom ............. 20-23 25(reggeon)
3914c clu ... 10 11 hard pom low mass .... 30
3915c ptn ... 20 21 proj remnant ......... 40-49
3916c str ... 29 targ remnant ......... 50-59
3917c pom ... 30 31 32(virtual) cluster .............. 60
3918c rem ... 40 41 direct photon ........ 71,72
3919c----------------------------------------------------------------------
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,*)' '
3974c if(istptl(i).ne.12)write(ifch,*)' '
3975c if(istptl(i).eq.12)write(ifch,'(1x,3(e9.3,1x))')
3976c & sptl(i),sqrt(uptl(i)-xorptl(1,i)**2)
3977c & ,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
3992c----------------------------------------------------------------------
3993 subroutine blist(text,n1,n2)
3994c----------------------------------------------------------------------
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
4058c----------------------------------------------------------------------
4059 subroutine alistf(text)
4060c----------------------------------------------------------------------
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
4106c &,(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
4139c----------------------------------------------------------------------
4140 subroutine alist2(text,n1,n2,n3,n4)
4141c----------------------------------------------------------------------
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)
4155c &,(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)
4161c &,(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))
4164c *,4x,4(e8.2,1x))
4165 end
4166
4167c----------------------------------------------------------------------
4168 subroutine alistc(text,n1,n2)
4169c----------------------------------------------------------------------
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
4189c-----------------------------------------------------------------------
4190 subroutine sigmaint(g0,gz,sigdo)
4191c-----------------------------------------------------------------------
4192c hadron-hadron cross sections integration
4193c-----------------------------------------------------------------------
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
4202c 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.)
4208c 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))
4229c vvv11=sqrt(vvv21)
4230c 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
4267c-----------------------------------------------------------------------
4268 subroutine xsigma
4269c-----------------------------------------------------------------------
4270c hadron-hadron and hadron-nucleus cross sections calculation
4271c b - impact parameter squared (in case of hadron-nucleus interaction);
4272c-----------------------------------------------------------------------
4273 include 'epos.inc'
4274 include 'epos.incsem'
4275 double precision gz(0:3),gzp(0:3),GZ0(2)
4276c 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
4284c 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
4292C...Total cross sections in Pythia
4293c double precision SIGT
4294c COMMON/PYINT7/SIGT(0:6,0:6,0:5)
4295
4296c Model 5 Common
4297 COMMON/HIPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50)
4298
4299c 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.
4313c 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))
4359c sigdelaf=max(0.,(sigelaf-sigelaold))
4360
4361 if(rexdifi(iclpro).lt.0..or.rexdifi(icltar).lt.0.)then
4362 sigdela=min(sigdelaf,sigdif)
4363
4364c 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
4392c rexdif(icltar)=max(exp(-1.7/engy**0.3),abs(rexdifi(icltar))) !strong reduction
4393c 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
4415c rexdif(iclpro)=max(exp(-1.7/engy**0.3),abs(rexdifi(iclpro))) !strong reduction
4416c 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
4451c 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
4457c 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
4470c GDP - projectile diffraction cross section
4471 GDP=(1.D0-CC(ICZ))*CC(2)*gz(0)
4472c GDT - target diffraction cross section
4473 GDT=(1.D0-CC(2))*CC(ICZ)*gz(0)
4474c 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
4494c$$$ elseif(model.eq.4)then
4495c$$$
4496c$$$ sigsd=real(SIGT(0,0,2)+SIGT(0,0,3))
4497c$$$ sigela=real(SIGT(0,0,1))
4498c$$$ sigcut=real(SIGT(0,0,5))
4499c$$$ sigtot=real(SIGT(0,0,0))
4500c$$$ sigine=sigtot-sigela
4501c$$$ 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
4583225 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