Adding TAmpt (Constantin)
[u/mrichter/AliRoot.git] / TAmpt / AMPT / amptsub.f
1 c....................amptsub.f
2 c.....this file contains 4 sections:
3 c.....1. ART subroutines;
4 c.....2. ART functions;
5 c.....3. ART block data;
6 c.....4. subprocesses borrowed from other codes.
7 c.....5. the previous artana.f
8 c.....6. the previous zpcsub.f
9 c.....7. subroutine getnp
10 c.....Note that Parts1-4 are the previous artsub.f
11 c
12 c=======================================================================
13 c.....subroutine to set up ART parameters and analysis files
14 c.....before looping different events
15 cms
16 cms   dlw & gsfs 8/2009 commented out lots of output files
17 cms
18       SUBROUTINE ARTSET
19 c
20       PARAMETER (AMU= 0.9383)
21       double precision dpcoal,drcoal,ecritl
22       INTEGER ZTA, ZPR
23       common  /gg/      dx,dy,dz,dpx,dpy,dpz
24 clin-10/03/03 
25 c     "SAVE   " (without argument) is used for most subroutines and functions,
26 c     this is important for the success when using "f77" to compile:
27 cc      SAVE /gg/
28       common  /zz/      zta,zpr
29 cc      SAVE /zz/
30       COMMON  /RUN/     NUM
31 cc      SAVE /RUN/
32       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
33 cc      SAVE /input1/
34       COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, 
35      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
36 cc      SAVE /INPUT2/
37       COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
38 cc      SAVE /INPUT3/
39       common /imulst/ iperts
40 cc      SAVE /imulst/
41       common /coal/dpcoal,drcoal,ecritl
42       common/anim/nevent,isoft,isflag,izpc
43       common /para7/ ioscar,nsmbbbar,nsmmeson
44       SAVE   
45 clin-10/03/03  ecritl: local energy density below which a parton 
46 c     will freeze out (in GeV/fm^3), for improvements on string melting, 
47 c     not used in this version of AMPT:
48 clin-4/2008
49 c      data ecritl/1.d0/
50       ecritl=1.d0
51 c
52 c     combine ART initialization into ampt.ini:
53 c     (Note that the following values are relics from the old ART structure)
54 c.....input parameter file
55 c      OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN')
56 c      READ (13, *) MASSTA, ZTA
57       MASSTA=1
58       ZTA=1
59 c      write(12,*) massta, zta, ' massta, zta'
60 c      READ (13, *) MASSPR, ZPR
61       MASSPR=1
62       ZPR=1
63 c      write(12,*) masspr, zpr, ' masspr, zpr'
64 c      READ (13, *) PLAB, IPLAB
65       PLAB=14.6 
66       IPLAB=2
67 c      write(12,*) plab, iplab, ' plab, iplab'
68       if(iplab.eq.2)then
69          elab=sqrt(plab**2+amu**2)-amu
70       else
71          elab=plab
72       endif
73       elab=elab*1000.
74 c      READ (13, *) ZEROPT
75       ZEROPT=0.
76 c      write(12,*) zeropt, ' zeropt'
77 clin-10/03/03 ISEED was used as a seed for random number inside ART, 
78 c     not used in AMPT:
79       ISEED=700721
80 c     0/1: (Normal or Perturbative) multistrange partice production.
81 c     Perturbative option is disabled for now:
82       iperts=0
83 c      READ (13, *) MANYB, B0, BI, BM
84 c     2/04/00 MANYB MUST BE SET TO 1 !
85 c     in order to skip impact parameter setting by ART, then B0 has no effect.
86       MANYB=1
87       B0=1
88       BI=0
89       BM=0
90 c      write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm'
91 c      READ (13, *) ISEED
92 c      write(12,*) iseed, ' iseed'
93 c      READ (13, *) DT
94 c      write(12,*) dt, ' dt'
95 c      READ (13, *) NTMAX
96 c      write(12,*) ntmax, ' ntmax'
97 c      READ (13, *) ICOLL
98       ICOLL=-1
99 c      write(12,*) icoll, ' icoll'
100 c      READ (13, *) NUM
101 c     2/11/03 run events without test particles for now:
102       NUM=1
103 c      write(12,*) num, ' num'
104 c      READ (13, *) INSYS
105       INSYS=1
106 c      write(12,*) insys, ' insys'
107 c      READ (13, *) IPOT
108       IPOT=3
109 c      write(12,*) ipot, ' ipot'
110 c      READ (13, *) MODE
111       MODE=0
112       IF(ICOLL.EQ.-1)IPOT=0
113 c      write(12,*) mode, ' mode'
114 c      READ (13, *) DX, DY, DZ
115       DX=2.73
116       DY=2.73
117       DZ=2.73
118 c      write(12,*) dx,dy,dz,' dx,dy,dz'
119 c      READ (13, *) DPX, DPY, DPZ
120       DPX=0.6
121       DPY=0.6
122       DPZ=0.6
123 c      write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz'
124 c      READ (13, *) IAVOID
125       IAVOID=1
126 c      write(12,*) iavoid, ' iavoid'
127 c      READ (13, *) IMOMEN
128       IMOMEN=1
129 c      write(12,*) imomen, ' imomen'
130       if(icoll.eq.-1)imomen=3
131 c      READ (13, *) NFREQ
132       NFREQ=10
133 c      write(12,*) nfreq, ' nfreq'
134 c      READ (13, *) ICFLOW
135       ICFLOW=0
136 c      write(12,*) ICFLOW, ' ICFLOW'
137 c      READ (13, *) ICRHO
138       ICRHO=0
139 c      write(12,*) ICRHO, ' ICRHO'
140 c      READ (13, *) ICOU
141       ICOU=0
142 c      write(12,*)icou, ' icou'
143 * kaon potential control parameter
144 * KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH
145 c      READ (13, *) KPOTEN, KMUL
146       KPOTEN=0
147       KMUL=1
148 c      write(12,*)kpoten,kmul, ' kpoten, kmul'
149 * mean field control parameter FOR BARYONS
150 * no mean filed is used for baryons if their 
151 * local density is higher than dencut. 
152 c      READ (13, *) DENCUT
153       DENCUT=15
154 c      write(12,*)dencut, ' dencut'
155 * test reactions in a box of side-length cycbox
156 * input cycbox
157 c      READ (13, *) CYCBOX
158       CYCBOX=0
159 c      write(12,*) cycbox, ' cycbox'
160 c
161 clin-5b/2008
162 c      if(ioscar.eq.2) then
163       if(ioscar.eq.2.or.ioscar.eq.3) then
164 cms      OPEN (92,FILE='ana/parton-initial-afterPropagation.dat',
165 cms  1        STATUS = 'UNKNOWN')
166       endif
167       if(ioscar.eq.3) then
168 clin-6/2009 write out full parton collision history:
169 cms      OPEN (95,FILE='ana/parton-collisionsHistory.dat',
170 cms  1        STATUS='UNKNOWN')
171 clin-6/2009 write out initial minijet information:
172 cms      OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat',
173 cms  1        STATUS='UNKNOWN')
174 clin-6/2009 write out parton info after coalescence:
175          if(isoft.eq.4.or.isoft.eq.5) then
176 cms         OPEN (85,FILE='ana/parton-after-coalescence.dat',
177 cms  1           STATUS='UNKNOWN')
178          endif
179       endif
180 clin-6/2009 write out initial transverse positions of initial nucleons:
181 cms   OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN')
182
183       RETURN
184       END
185
186 c-----------------------------------------------------------------------
187
188 c.....subroutine to initialize cascade.
189
190       SUBROUTINE ARINI
191
192 c.....before invoking ARINI:
193 c.....IAPAR2(1), IAINT2(1) must be set.
194       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
195 cc      SAVE /ARPRNT/
196       SAVE   
197
198 ctest off for resonance (phi, K*) studies:
199 c      OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN')
200
201       IFLG = IAPAR2(1)
202       GOTO (200, 200, 300) IFLG
203
204 c.....error choice of initialization
205       PRINT *, 'IAPAR2(1) must be 1, 2, or 3'
206       STOP
207
208 c.....to use default initial conditions generated by the cascade,
209 c.....or to read in initial conditions.
210  200  RETURN
211
212 c.....to generate formation time and the position at formation time from 
213 c.....read-in initial conditions with an averaged formation proper time.
214  300  CALL ARINI1
215 c.....ordering the particle label according to increasing order of 
216 c.....formation time.
217       CALL ARTORD
218       RETURN
219
220       END
221
222 c-----------------------------------------------------------------------
223
224 c.....subroutine to generate formation time and position at formation time
225 c.....from read-in initial conditions with an averaged formation proper 
226 c.....time.
227
228       SUBROUTINE ARINI1
229
230 c.....before invoking ARINI1:
231 c.....ARPAR1(1), IAINT2(1) must be set:
232       PARAMETER (MAXSTR=150001)
233       double precision  smearp,smearh
234
235       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
236 cc      SAVE /ARPRNT/
237       COMMON /ARPRC/ ITYPAR(MAXSTR),
238      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
239      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
240      &     XMAR(MAXSTR)
241 cc      SAVE /ARPRC/
242       COMMON /smearz/smearp,smearh
243 cc      SAVE /smearz/
244       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
245 cc      SAVE /input1/
246       common/anim/nevent,isoft,isflag,izpc
247 cc      SAVE /anim/
248       common /nzpc/nattzp
249 cc      SAVE /nzpc/
250       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
251 cc      SAVE /HPARNT/
252       COMMON/RNDF77/NSEED
253 cc      SAVE /RNDF77/
254       common /para8/ idpert,npertd,idxsec
255       SAVE   
256 clin-5/2008 for perturbatively-produced hadrons (currently only deuterons):
257 cms   OPEN (91, FILE = 'ana/deuteron_processes.dat', 
258 cms  1     STATUS = 'UNKNOWN')
259       if(idpert.eq.1.or.idpert.eq.2) then
260 cms      OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN')
261       endif
262 c.....generate formation time and position at formation time.
263       TAU0 = ARPAR1(1)
264       NP = IAINT2(1)
265 clin-7/10/01     initial positions already given for hadrons 
266 c     formed from partons inside ZPC (from string melting):
267       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
268          if(NP.le.nattzp) return
269          do 1001 I = nattzp+1, NP
270             IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
271                PRINT *, ' IN ARINI1'
272                PRINT *, 'ABS(PZ) .GE. EE for particle ', I
273                PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), 
274      &              ' PY = ', PYAR(I)
275                PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
276                PRINT *, ' XM = ', XMAR(I)
277                RAP = 1000000.0
278                GOTO 50
279             END IF
280             RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I)))
281  50         CONTINUE
282             VX = PXAR(I) / PEAR(I)
283             VY = PYAR(I) / PEAR(I)
284             FTAR(I) = TAU0 * COSH(RAP)
285             GXAR(I) = GXAR(I) + VX * FTAR(I)
286             GYAR(I) = GYAR(I) + VY * FTAR(I)
287             GZAR(I) = TAU0 * SINH(RAP)
288 clin-5/2009 No formation time for spectator projectile or target nucleons:
289             if(PXAR(I).eq.0.and.PYAR(I).eq.0
290      1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
291      2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
292                TAUI=1.E-20
293                FTAR(I)=TAUI*COSH(RAP)
294                GZAR(I)=TAUI*SINH(RAP)
295             endif
296  1001    continue
297 clin-7/10/01-end
298 clin-3/2009 cleanup of program flow:
299       else
300          DO 1002 I = 1, NP
301             IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
302                PRINT *, ' IN ARINI1'
303                PRINT *, 'ABS(PZ) .GE. EE for particle ', I
304                PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), 
305      &              ' PY = ', PYAR(I)
306                PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
307                PRINT *, ' XM = ', XMAR(I)
308                RAP = 1000000.0
309                GOTO 100
310 c               STOP
311             END IF
312             RAP = 0.5 * LOG((PEAR(I) + PZAR(I)) / (PEAR(I) - PZAR(I)))
313  100        CONTINUE
314             VX = PXAR(I) / PEAR(I)
315             VY = PYAR(I) / PEAR(I)
316 c.....give initial formation time shift
317             TAUI = FTAR(I) + TAU0
318             FTAR(I) = TAUI * COSH(RAP)
319             GXAR(I) = GXAR(I) + VX * TAU0 * COSH(RAP)
320             GYAR(I) = GYAR(I) + VY * TAU0 * COSH(RAP)
321 c     4/25/03: hadron z-position upon formation determined the same way as x,y:
322             GZAR(I) = TAUI * SINH(RAP)
323 c     the old prescription:
324 c            GZAR(I) = GZAR(I) + TAU0 * SINH(RAP)
325             zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
326             GZAR(I)=GZAR(I)+zsmear
327 cbz1/28/99end
328 c     10/05/01 no formation time for spectator projectile or target nucleons:
329             if(PXAR(I).eq.0.and.PYAR(I).eq.0
330      1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
331      2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
332 clin-5/2008:
333 c               TAUI=0.00001
334                TAUI=1.E-20
335                FTAR(I)=TAUI*COSH(RAP)
336                GZAR(I)=TAUI*SINH(RAP)+zsmear
337             endif
338  1002    CONTINUE
339 clin-3/2009 cleanup of program flow:
340       endif
341
342 clin-3/2009 Add initial hadrons before the hadron cascade starts:
343       call addhad
344
345       RETURN
346       END
347
348 c-----------------------------------------------------------------------
349
350 c.....subroutine to order particle labels according to increasing 
351 c.....formation time
352
353       SUBROUTINE ARTORD
354
355 c.....before invoking ARTORD:
356 c.....IAINT2(1) must be set:
357       PARAMETER (MAXSTR=150001,MAXR=1)
358       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
359 cc      SAVE /ARPRNT/
360       COMMON /ARPRC/ ITYPAR(MAXSTR),
361      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
362      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
363      &     XMAR(MAXSTR)
364 cc      SAVE /ARPRC/
365 clin-3/2009 Take care of particle weights when user inserts initial hadrons:
366       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
367      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
368      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
369       DIMENSION dptemp(MAXSTR)
370 c
371       DIMENSION ITYP0(MAXSTR), 
372      &   GX0(MAXSTR), GY0(MAXSTR), GZ0(MAXSTR), FT0(MAXSTR),
373      &   PX0(MAXSTR), PY0(MAXSTR), PZ0(MAXSTR), EE0(MAXSTR),
374      &   XM0(MAXSTR)
375       DIMENSION INDX(MAXSTR)
376       EXTERNAL ARINDX
377       SAVE   
378 c
379       NPAR = 0
380       NP = IAINT2(1)
381       DO 1001 I = 1, NP
382          ITYP0(I) = ITYPAR(I)
383          GX0(I) = GXAR(I)
384          GY0(I) = GYAR(I)
385          GZ0(I) = GZAR(I)
386          FT0(I) = FTAR(I)
387          PX0(I) = PXAR(I)
388          PY0(I) = PYAR(I)
389          PZ0(I) = PZAR(I)
390          EE0(I) = PEAR(I)
391          XM0(I) = XMAR(I)
392 clin-3/2009:
393          dptemp(I) = dpertp(I)
394  1001 CONTINUE
395       CALL ARINDX(MAXSTR, NP, FT0, INDX)
396       DO 1002 I = 1, NP
397 cbz12/3/98
398 c         IF (ITYP0(INDX(I)) .EQ. 211) THEN
399 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN
400 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
401 c     &      ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR.
402 c     &      ITYP0(INDX(I)) .EQ. 111) THEN
403 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
404 c     &      ITYP0(INDX(I)) .EQ. 2112) THEN
405          NPAR = NPAR + 1
406 c         ITYPAR(I) = ITYP0(INDX(I))
407 c         GXAR(I) = GX0(INDX(I))
408 c         GYAR(I) = GY0(INDX(I))
409 c         GZAR(I) = GZ0(INDX(I))
410 c         FTAR(I) = FT0(INDX(I))
411 c         PXAR(I) = PX0(INDX(I))
412 c         PYAR(I) = PY0(INDX(I))
413 c         PZAR(I) = PZ0(INDX(I))
414 c         PEAR(I) = EE0(INDX(I))
415 c         XMAR(I) = XM0(INDX(I))
416          ITYPAR(NPAR) = ITYP0(INDX(I))
417          GXAR(NPAR) = GX0(INDX(I))
418          GYAR(NPAR) = GY0(INDX(I))
419          GZAR(NPAR) = GZ0(INDX(I))
420          FTAR(NPAR) = FT0(INDX(I))
421          PXAR(NPAR) = PX0(INDX(I))
422          PYAR(NPAR) = PY0(INDX(I))
423          PZAR(NPAR) = PZ0(INDX(I))
424          PEAR(NPAR) = EE0(INDX(I))
425          XMAR(NPAR) = XM0(INDX(I))
426 clin-3/2009:
427          dpertp(NPAR)=dptemp(INDX(I))
428 c         END IF
429 cbz12/3/98end
430  1002 CONTINUE
431       IAINT2(1) = NPAR
432 c
433       RETURN
434       END
435
436 c-----------------------------------------------------------------------
437
438 c.....subroutine to copy individually generated particle record into
439 c.....particle record for many test particle runs.
440
441       SUBROUTINE ARINI2(K)
442
443       PARAMETER (MAXSTR=150001,MAXR=1)
444       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
445 cc      SAVE /ARPRNT/
446       COMMON /ARPRC/ ITYPAR(MAXSTR),
447      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
448      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
449      &     XMAR(MAXSTR)
450 cc      SAVE /ARPRC/
451       COMMON /ARERC1/MULTI1(MAXR)
452 cc      SAVE /ARERC1/
453       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
454      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
455      &     FT1(MAXSTR, MAXR),
456      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
457      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
458 cc      SAVE /ARPRC1/
459       COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
460 cc      SAVE /tdecay/
461       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
462 cc      SAVE /input1/
463       COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, 
464      &     IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
465 cc      SAVE /INPUT2/
466       COMMON/RNDF77/NSEED
467       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
468      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
469      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
470 cc      SAVE /RNDF77/
471       SAVE   
472
473       MULTI1(K) = IAINT2(1)
474       DO 1001 I = 1, MULTI1(K)
475          ITYP1(I, K) = ITYPAR(I)
476          GX1(I, K) = GXAR(I)
477          GY1(I, K) = GYAR(I)
478          GZ1(I, K) = GZAR(I)
479          FT1(I, K) = FTAR(I)
480          PX1(I, K) = PXAR(I)
481          PY1(I, K) = PYAR(I)
482          PZ1(I, K) = PZAR(I)
483          EE1(I, K) = PEAR(I)
484          XM1(I, K) = XMAR(I)
485 clin-3/2009 hadron weights are initialized in addhad():
486 clin-5/2008 all hadrons not perturbatively-produced have the weight of 1:
487 c         dpp1(I,K)=1.
488          dpp1(I,K)=dpertp(I)
489  1001 CONTINUE
490
491 c     initialize final time of each particle to ntmax*dt except for 
492 c     decay daughters, which have values given by tfdcy() and >(ntmax*dt):
493       do 1002 ip=1,MAXSTR
494          tfdcy(ip)=NTMAX*DT
495          tft(ip)=NTMAX*DT
496  1002 continue
497 c
498       do 1004 irun=1,MAXR
499          do 1003 ip=1,MAXSTR
500             tfdpi(ip,irun)=NTMAX*DT
501  1003    continue
502  1004 continue
503
504       RETURN
505       END
506
507 c=======================================================================
508
509 c.....function to convert PDG flavor code into ART flavor code.
510
511       FUNCTION IARFLV(IPDG)
512
513       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
514 cc      SAVE /input1/
515       COMMON/RNDF77/NSEED
516 cc      SAVE /RNDF77/
517       SAVE   
518
519 c.....anti-Delta-
520       IF (IPDG .EQ. -1114) THEN
521          IARFLV = -6
522          RETURN
523       END IF
524
525 c.....anti-Delta0
526       IF (IPDG .EQ. -2114) THEN
527          IARFLV = -7
528          RETURN
529       END IF
530
531 c.....anti-Delta+
532       IF (IPDG .EQ. -2214) THEN
533          IARFLV = -8
534          RETURN
535       END IF
536
537 c.....anti-Delta++
538       IF (IPDG .EQ. -2224) THEN
539          IARFLV = -9
540          RETURN
541       END IF
542
543 cbzdbg2/23/99
544 c.....anti-proton
545       IF (IPDG .EQ. -2212) THEN
546          IARFLV = -1
547          RETURN
548       END IF
549
550 c.....anti-neutron
551       IF (IPDG .EQ. -2112) THEN
552          IARFLV = -2
553          RETURN
554       END IF
555 cbzdbg2/23/99end
556
557 c.....eta
558       IF (IPDG .EQ. 221) THEN
559          IARFLV = 0
560          RETURN
561       END IF
562
563 c.....proton
564       IF (IPDG .EQ. 2212) THEN
565          IARFLV = 1
566          RETURN
567       END IF
568
569 c.....neutron
570       IF (IPDG .EQ. 2112) THEN
571          IARFLV = 2
572          RETURN
573       END IF
574
575 c.....pi-
576       IF (IPDG .EQ. -211) THEN
577          IARFLV = 3
578          RETURN
579       END IF
580
581 c.....pi0
582       IF (IPDG .EQ. 111) THEN
583          IARFLV = 4
584          RETURN
585       END IF
586
587 c.....pi+
588       IF (IPDG .EQ. 211) THEN
589          IARFLV = 5
590          RETURN
591       END IF
592
593 c.....Delta-
594       IF (IPDG .EQ. 1114) THEN
595          IARFLV = 6
596          RETURN
597       END IF
598
599 c.....Delta0
600       IF (IPDG .EQ. 2114) THEN
601          IARFLV = 7
602          RETURN
603       END IF
604
605 c.....Delta+
606       IF (IPDG .EQ. 2214) THEN
607          IARFLV = 8
608          RETURN
609       END IF
610
611 c.....Delta++
612       IF (IPDG .EQ. 2224) THEN
613          IARFLV = 9
614          RETURN
615       END IF
616
617 c.....Lambda
618       IF (IPDG .EQ. 3122) THEN
619          IARFLV = 14
620          RETURN
621       END IF
622
623 c.....Lambda-bar
624       IF (IPDG .EQ. -3122) THEN
625          IARFLV = -14
626          RETURN
627       END IF
628
629 c.....Sigma-
630       IF (IPDG .EQ. 3112) THEN
631          IARFLV = 15
632          RETURN
633       END IF
634
635 c.....Sigma-bar
636       IF (IPDG .EQ. -3112) THEN
637          IARFLV = -15
638          RETURN
639       END IF 
640
641 c.....Sigma0
642       IF (IPDG .EQ. 3212) THEN
643          IARFLV = 16
644          RETURN
645       END IF
646
647 c.....Sigma0-bar
648       IF (IPDG .EQ. -3212) THEN
649          IARFLV = -16
650          RETURN
651       END IF 
652
653 c.....Sigma+
654       IF (IPDG .EQ. 3222) THEN
655          IARFLV = 17
656          RETURN
657       END IF
658
659 c.....Sigma+ -bar
660       IF (IPDG .EQ. -3222) THEN
661          IARFLV = -17
662          RETURN
663       END IF 
664
665 c.....K-
666       IF (IPDG .EQ. -321) THEN
667          IARFLV = 21
668          RETURN
669       END IF
670
671 c.....K+
672       IF (IPDG .EQ. 321) THEN
673          IARFLV = 23
674          RETURN
675       END IF
676
677 c.....temporary entry for K0
678       IF (IPDG .EQ. 311) THEN
679          IARFLV = 23
680          RETURN
681       END IF
682
683 c.....temporary entry for K0bar
684       IF (IPDG .EQ. -311) THEN
685          IARFLV = 21
686          RETURN
687       END IF
688
689 c.....temporary entry for K0S and K0L
690       IF (IPDG .EQ. 310 .OR. IPDG .EQ. 130) THEN
691          R = RANART(NSEED)
692          IF (R .GT. 0.5) THEN
693             IARFLV = 23
694          ELSE
695             IARFLV = 21
696          END IF
697          RETURN
698       END IF
699
700 c.....rho-
701       IF (IPDG .EQ. -213) THEN
702          IARFLV = 25
703          RETURN
704       END IF
705
706 c.....rho0
707       IF (IPDG .EQ. 113) THEN
708          IARFLV = 26
709          RETURN
710       END IF
711
712 c.....rho+
713       IF (IPDG .EQ. 213) THEN
714          IARFLV = 27
715          RETURN
716       END IF
717
718 c.....omega
719       IF (IPDG .EQ. 223) THEN
720          IARFLV = 28
721          RETURN
722       END IF
723
724 c.....phi
725       IF (IPDG .EQ. 333) THEN
726          IARFLV = 29
727          RETURN
728       END IF
729
730 c.....K*+
731       IF (IPDG .EQ. 323) THEN
732          IARFLV = 30
733          RETURN
734       END IF
735 c.....K*-
736       IF (IPDG .EQ. -323) THEN
737          IARFLV = -30
738          RETURN
739       END IF
740 c.....temporary entry for K*0
741       IF (IPDG .EQ. 313) THEN
742          IARFLV = 30
743          RETURN
744       END IF
745 c.....temporary entry for K*0bar
746       IF (IPDG .EQ. -313) THEN
747          IARFLV = -30
748          RETURN
749       END IF
750
751 c...... eta-prime
752       IF (IPDG .EQ. 331) THEN
753          IARFLV = 31
754          RETURN
755       END IF
756  
757 c...... a1
758 c     IF (IPDG .EQ. 777) THEN
759 c        IARFLV = 32
760 c        RETURN
761 c     END IF
762                                 
763 c... cascade-
764       IF (IPDG .EQ. 3312) THEN
765          IARFLV = 40
766          RETURN
767       END IF
768  
769 c... cascade+ (bar)
770       IF (IPDG .EQ. -3312) THEN
771          IARFLV = -40
772          RETURN
773       END IF
774  
775 c... cascade0
776       IF (IPDG .EQ. 3322) THEN
777          IARFLV = 41
778          RETURN
779       END IF
780  
781 c... cascade0 -bar
782       IF (IPDG .EQ. -3322) THEN
783          IARFLV = -41
784          RETURN
785       END IF
786  
787 c... Omega-
788       IF (IPDG .EQ. 3334) THEN
789          IARFLV = 45
790          RETURN
791       END IF 
792
793 c... Omega+ (bar)
794       IF (IPDG .EQ. -3334) THEN
795          IARFLV = -45
796          RETURN
797       END IF
798
799 c... Di-Omega
800       IF (IPDG .EQ. 6666) THEN
801          IARFLV = 44
802          RETURN
803       END IF
804 c sp06/05/01 end    
805
806 clin-3/2009 keep the same ID numbers in case there are initial deuterons:
807       IF (IPDG .EQ. 42 .or. IPDG .EQ. -42) THEN
808          IARFLV = IPDG
809          RETURN
810       END IF
811
812 c.....other
813       IARFLV = IPDG + 10000
814
815       RETURN
816       END
817
818 c-----------------------------------------------------------------------
819
820 c.....function to convert ART flavor code into PDG flavor code.
821
822       FUNCTION INVFLV(IART)
823
824       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
825 cc      SAVE /input1/
826       COMMON/RNDF77/NSEED
827 cc      SAVE /RNDF77/
828       SAVE   
829
830 c.....anti-Delta-
831       IF (IART .EQ. -6) THEN
832          INVFLV = -1114
833          RETURN
834       END IF
835
836 c.....anti-Delta0
837       IF (IART .EQ. -7) THEN
838          INVFLV = -2114
839          RETURN
840       END IF
841
842 c.....anti-Delta+
843       IF (IART .EQ. -8) THEN
844          INVFLV = -2214
845          RETURN
846       END IF
847
848 c.....anti-Delta++
849       IF (IART .EQ. -9) THEN
850          INVFLV = -2224
851          RETURN
852       END IF
853
854 cbzdbg2/23/99
855 c.....anti-proton
856       IF (IART .EQ. -1) THEN
857          INVFLV = -2212
858          RETURN
859       END IF
860
861 c.....anti-neutron
862       IF (IART .EQ. -2) THEN
863          INVFLV = -2112
864          RETURN
865       END IF
866 cbzdbg2/23/99end
867
868 c.....eta
869       IF (IART .EQ. 0) THEN
870          INVFLV = 221
871          RETURN
872       END IF
873
874 c.....proton
875       IF (IART .EQ. 1) THEN
876          INVFLV = 2212
877          RETURN
878       END IF
879
880 c.....neutron
881       IF (IART .EQ. 2) THEN
882          INVFLV = 2112
883          RETURN
884       END IF
885
886 c.....pi-
887       IF (IART .EQ. 3) THEN
888          INVFLV = -211
889          RETURN
890       END IF
891
892 c.....pi0
893       IF (IART .EQ. 4) THEN
894          INVFLV = 111
895          RETURN
896       END IF
897
898 c.....pi+
899       IF (IART .EQ. 5) THEN
900          INVFLV = 211
901          RETURN
902       END IF
903
904 c.....Delta-
905       IF (IART .EQ. 6) THEN
906          INVFLV = 1114
907          RETURN
908       END IF
909
910 c.....Delta0
911       IF (IART .EQ. 7) THEN
912          INVFLV = 2114
913          RETURN
914       END IF
915
916 c.....Delta+
917       IF (IART .EQ. 8) THEN
918          INVFLV = 2214
919          RETURN
920       END IF
921
922 c.....Delta++
923       IF (IART .EQ. 9) THEN
924          INVFLV = 2224
925          RETURN
926       END IF
927
928 cc.....N*(1440), N*(1535) temporary entry
929 c      IF (IART .GE. 10 .AND. IART .LE.13) THEN
930 c         INVFLV = 0
931 c         RETURN
932 c      END IF
933
934 c.....Lambda
935       IF (IART .EQ. 14) THEN
936          INVFLV = 3122
937          RETURN
938       END IF
939 c.....Lambda-bar
940       IF (IART .EQ. -14) THEN
941          INVFLV = -3122
942          RETURN
943       END IF 
944
945 cbz3/12/99
946 c.....temporary entry for Sigma's
947 c      IF (IART .EQ. 15) THEN
948 c         R = RANART(NSEED)
949 c         IF (R .GT. 2. / 3.) THEN
950 c            INVFLV = 3112
951 c         ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN
952 c            INVFLV = 3212
953 c         ELSE
954 c            INVFLV = 3222
955 c         END IF
956 c         RETURN
957 c      END IF
958
959 c.....Sigma-
960       IF (IART .EQ. 15) THEN
961          INVFLV = 3112
962          RETURN
963       END IF
964
965 c.....Sigma- bar
966       IF (IART .EQ. -15) THEN
967          INVFLV = -3112
968          RETURN
969       END IF 
970
971 c.....Sigma0
972       IF (IART .EQ. 16) THEN
973          INVFLV = 3212
974          RETURN
975       END IF
976
977 c.....Sigma0 -bar
978       IF (IART .EQ. -16) THEN
979          INVFLV = -3212
980          RETURN
981       END IF
982
983 c.....Sigma+
984       IF (IART .EQ. 17) THEN
985          INVFLV = 3222
986          RETURN
987       END IF
988
989 c.....Sigma+ -bar
990       IF (IART .EQ. -17) THEN
991          INVFLV = -3222
992          RETURN
993       END IF 
994
995 clin-2/23/03 K0S and K0L are generated at the last timestep:
996 c.....temporary entry for K- and K0bar
997       IF (IART .EQ. 21) THEN
998 c         R = RANART(NSEED)
999 c         IF (R .GT. 0.5) THEN
1000             INVFLV = -321
1001 c         ELSE
1002 c            INVFLV = -311
1003 c            R = RANART(NSEED)
1004 c            IF (R .GT. 0.5) THEN
1005 c               INVFLV = 310
1006 c            ELSE
1007 c               INVFLV = 130
1008 c            END IF
1009 c         END IF
1010          RETURN
1011       END IF
1012
1013 c.....temporary entry for K+ and K0
1014       IF (IART .EQ. 23) THEN
1015 c         R = RANART(NSEED)
1016 c         IF (R .GT. 0.5) THEN
1017             INVFLV = 321
1018 c         ELSE
1019 c            INVFLV = 311
1020 c            R = RANART(NSEED)
1021 c            IF (R .GT. 0.5) THEN
1022 c               INVFLV = 310
1023 c            ELSE
1024 c               INVFLV = 130
1025 c            END IF
1026 c         END IF
1027          RETURN
1028       END IF
1029
1030 c.....K0Long:
1031       IF (IART .EQ. 22) THEN
1032          INVFLV = 130
1033          RETURN
1034       ENDIF
1035 c.....K0Short:
1036       IF (IART .EQ. 24) THEN
1037          INVFLV = 310
1038          RETURN
1039       ENDIF
1040
1041 c.....rho-
1042       IF (IART .EQ. 25) THEN
1043          INVFLV = -213
1044          RETURN
1045       END IF
1046
1047 c.....rho0
1048       IF (IART .EQ. 26) THEN
1049          INVFLV = 113
1050          RETURN
1051       END IF
1052
1053 c.....rho+
1054       IF (IART .EQ. 27) THEN
1055          INVFLV = 213
1056          RETURN
1057       END IF
1058
1059 c.....omega
1060       IF (IART .EQ. 28) THEN
1061          INVFLV = 223
1062          RETURN
1063       END IF
1064
1065 c.....phi
1066       IF (IART .EQ. 29) THEN
1067          INVFLV = 333
1068          RETURN
1069       END IF
1070
1071 c.....temporary entry for K*+ and K*0
1072       IF (IART .EQ. 30) THEN
1073          INVFLV = 323
1074          IF (RANART(NSEED).GT.0.5) INVFLV = 313
1075          RETURN
1076       END IF
1077
1078 c.....temporary entry for K*- and K*0bar
1079       IF (IART .EQ. -30) THEN
1080          INVFLV = -323
1081          IF (RANART(NSEED).GT.0.5) INVFLV = -313
1082          RETURN
1083       END IF
1084
1085 c... eta-prime (bar)
1086       IF (IART .EQ. 31) THEN
1087          INVFLV = 331
1088          RETURN
1089       END IF
1090  
1091 c... a1
1092       IF (IART .EQ. 32) THEN
1093          INVFLV = 777
1094          RETURN
1095       END IF
1096  
1097 c... cascade-
1098       IF (IART .EQ. 40) THEN
1099          INVFLV = 3312
1100          RETURN
1101       END IF                   
1102
1103 c... cascade+ (bar)
1104       IF (IART .EQ. -40) THEN
1105          INVFLV = -3312
1106          RETURN
1107       END IF
1108  
1109 c... cascade0
1110       IF (IART .EQ. 41) THEN
1111          INVFLV = 3322
1112          RETURN
1113       END IF
1114  
1115 c... cascade0 -bar
1116       IF (IART .EQ. -41) THEN
1117          INVFLV = -3322
1118          RETURN
1119       END IF
1120  
1121 c... Omega-
1122       IF (IART .EQ. 45) THEN
1123          INVFLV = 3334
1124          RETURN
1125       END IF
1126
1127 c... Omega+ (bar)
1128       IF (IART .EQ. -45) THEN
1129          INVFLV = -3334
1130          RETURN
1131       END IF
1132
1133 c... Di-Omega
1134       IF (IART .EQ. 44) THEN
1135          INVFLV = 6666
1136          RETURN
1137       END IF
1138 c sp 12/19/00 end           
1139
1140 clin-5/2008 deuteron ID numbers in ART and ampt.dat:
1141       IF (IART .EQ. 42) THEN
1142          INVFLV = 1000000000+10020
1143          RETURN
1144       ELSEIF (IART .EQ. -42) THEN         
1145          INVFLV = -1000000000-10020
1146          RETURN
1147       END IF
1148 c
1149 c.....other
1150       INVFLV = IART - 10000
1151
1152       RETURN
1153       END
1154
1155 c=======================================================================
1156
1157       BLOCK DATA ARDATA
1158
1159       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
1160 cc      SAVE /ARPRNT/
1161       SAVE   
1162       DATA ARPAR1/1.19, 99 * 0.0/
1163       DATA IAPAR2/3, 49 * 0/
1164       DATA ARINT1/100 * 0.0/
1165       DATA IAINT2/50 * 0/
1166
1167       END
1168
1169 c=======================================================================
1170
1171 c.....Routine borrowed from ZPC.
1172 c.....double precision  is modified to real*4.
1173
1174 cbz1/29/99
1175 c      subroutine index1(n, m, arrin, indx)
1176       subroutine arindx(n, m, arrin, indx)
1177 cbz1/29/99end
1178 c     indexes the first m elements of ARRIN of length n, i.e., outputs INDX
1179 c     such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
1180
1181 c      implicit real*4 (a-h, o-z)
1182
1183       dimension arrin(n), indx(n)
1184       SAVE   
1185       do 1001 j = 1, m
1186          indx(j) = j
1187  1001 continue
1188       l = m / 2 + 1
1189       ir = m
1190  10   continue
1191       if (l .gt. 1) then
1192          l = l - 1
1193          indxt = indx(l)
1194          q = arrin(indxt)
1195       else
1196          indxt = indx(ir)
1197          q = arrin(indxt)
1198          indx(ir) = indx(1)
1199          ir = ir - 1
1200          if (ir .eq. 1) then
1201             indx(1) = indxt
1202             return
1203          end if
1204       end if
1205       i = l
1206       j = l + l
1207  20   if (j .le. ir) then
1208          if (j .lt. ir) then
1209             if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
1210          end if
1211          if (q .lt. arrin(indx(j))) then
1212             indx(i) = indx(j)
1213             i = j
1214             j = j + j
1215          else
1216             j = ir + 1
1217          end if
1218       goto 20
1219       end if
1220       indx(i) = indxt
1221       goto 10
1222
1223       end
1224
1225 c-----------------------------------------------------------------------
1226
1227 c.....extracted from G. Song's ART expasion including K- interactions
1228 c.....file `NEWKAON.FOR'
1229
1230 c     5/01/03 send iblock value into art1f.f, necessary for resonance studies:
1231 c        subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1232 c     &                                   srt,pcx,pcy,pcz)
1233         subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1234      &                                   srt,pcx,pcy,pcz,iblock)
1235       PARAMETER      (MAXSTR=150001,MAXR=1)
1236       PARAMETER      (AKA=0.498)
1237       COMMON   /AA/  R(3,MAXSTR)
1238 cc      SAVE /AA/
1239       COMMON   /BB/  P(3,MAXSTR)
1240 cc      SAVE /BB/
1241       COMMON   /CC/  E(MAXSTR)
1242 cc      SAVE /CC/
1243       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
1244 cc      SAVE /EE/
1245       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
1246 cc      SAVE /BG/
1247       COMMON   /NN/NNN
1248 cc      SAVE /NN/
1249       COMMON   /RUN/NUM
1250 cc      SAVE /RUN/
1251       COMMON   /PA/RPION(3,MAXSTR,MAXR)
1252 cc      SAVE /PA/
1253       COMMON   /PB/PPION(3,MAXSTR,MAXR)
1254 cc      SAVE /PB/
1255       COMMON   /PC/EPION(MAXSTR,MAXR)
1256 cc      SAVE /PC/
1257       COMMON   /PD/LPION(MAXSTR,MAXR)
1258 cc      SAVE /PD/
1259       COMMON/RNDF77/NSEED
1260 cc      SAVE /RNDF77/
1261       SAVE   
1262 c
1263         logical lb1bn, lb2bn,lb1mn,lb2mn
1264 cbz3/7/99 neutralk
1265 c        logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0
1266         logical lb1bn1, lb2bn1, lb1bn0, lb2bn0
1267 cbz3/7/99 neutralk end
1268         logical lb1mn0, lb2mn0, lb1mn1, lb2mn1
1269         logical lb1mn2, lb2mn2
1270         icase=-1
1271 c        icase: flag for the type of reaction that is going to happen.
1272 c        icase=-1,  no desired reaction, return to main program.
1273 c              1,  NN,ND,DD
1274 c              2,  PI+N, PI+D
1275 c              3,  K(-) absorption.
1276         nchrg=-100
1277 c        nchrg: Net charges of the two incoming particles.
1278         ictrl = 1
1279         lb1=lb(i1)
1280         lb2=lb(i2)
1281         em1=e(i1)
1282         em2=e(i2)
1283         lb1bn=lb1.eq.1.or.lb1.eq.2.or.(lb1.gt.5.and.lb1.le.13)
1284         lb2bn=lb2.eq.1.or.lb2.eq.2.or.(lb2.gt.5.and.lb2.le.13)
1285         lb1bn0=lb1.eq.2.or.lb1.eq.7.or.lb1.eq.10.or.lb1.eq.12
1286         lb2bn0=lb2.eq.2.or.lb2.eq.7.or.lb2.eq.10.or.lb2.eq.12
1287         lb1bn1=lb1.eq.1.or.lb1.eq.8.or.lb1.eq.11.or.lb1.eq.13
1288         lb2bn1=lb2.eq.1.or.lb2.eq.8.or.lb2.eq.11.or.lb2.eq.13
1289         lb1mn=em1.lt.0.2.or.lb1.eq.0.or.(lb1.ge.25.and.lb1.le.29)
1290         lb2mn=em2.lt.0.2.or.lb2.eq.0.or.(lb2.ge.25.and.lb2.le.29)
1291         lb1mn0=lb1.eq.0.or.lb1.eq.4.or.lb1.eq.26.or.
1292      &                        lb1.eq.28.or.lb1.eq.29
1293         lb2mn0=lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
1294      &                        lb2.eq.28.or.lb2.eq.29
1295         lb1mn1= lb1.eq.5.or.lb1.eq.27
1296         lb2mn1= lb2.eq.5.or.lb2.eq.27
1297         lb1mn2=lb1.eq.3.or.lb1.eq.25
1298         lb2mn2=lb2.eq.3.or.lb2.eq.25
1299
1300 c        1. consider N+N, N+Resonance, R + R reactions
1301         if(lb1bn.and.lb2bn) then
1302 c     NN,ND,DD:
1303            icase=1
1304 c     total cross section
1305            sig=40.
1306            if(lb1.eq.9.and.lb2.eq.9) then
1307                 nchrg=4
1308            endif   
1309            if((lb1bn1.and.lb2.eq.9)
1310      &        .or.(lb2bn1.and.lb1.eq.9))then
1311                 nchrg=3
1312            endif
1313            if((lb1bn0.and.lb2.eq.9)
1314      &        .or.(lb2bn0.and.lb1.eq.9)
1315      &        .or.(lb1bn1.and.lb2bn1)) then
1316                    nchrg=2
1317            endif
1318            if((lb1bn1.and.lb2bn0).or.(lb1.eq.6.and.lb2.eq.9)
1319      &        .or.(lb2bn1.and.lb1bn0)
1320      &        .or.(lb2.eq.6.and.lb1.eq.9))then
1321                    nchrg=1
1322            endif
1323            if((lb1bn0.and.lb2bn0).or.(lb1bn1.and.lb2.eq.6)
1324      &              .or.(lb2bn1.and.lb1.eq.6)) then
1325                    nchrg=0
1326            endif
1327            if((lb1bn0.and.lb2.eq.6)
1328      &        .or.(lb2bn0.and.lb1.eq.6))then
1329                 nchrg=-1
1330            endif
1331            if(lb1.eq.6.and.lb2.eq.6) then
1332                 nchrg=-2
1333            endif
1334 c     brsig = x2kaon_no_isospin(srt)
1335            if(nchrg.ge.-1.and.nchrg.le.2) then
1336 c     K,Kbar prduction x sect.
1337                    brsig = x2kaon(srt)
1338            else
1339                    brsig=0.0
1340 c                if(nchrg.eq.-2.or.nchrg.eq.3) then
1341 c                   brsig = x2kaon(srt+0.938-1.232)
1342 c                else
1343 c     nchrg=4
1344 c                   brsig = x2kaon(srt+2.*(0.938-1.232))
1345 c                endif
1346            endif
1347
1348 cbz3/7/99 neutralk
1349            BRSIG = 2.0 * BRSIG
1350 cbz3/7/99 neutralk end
1351
1352         endif
1353
1354 c        2. consider PI(meson:eta,omega,rho,phi) + N(N*,D)
1355         if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then
1356 c     PN,PD
1357           icase=2
1358           sig=20.
1359           sigma0 = piNsg0(srt)
1360           brsig=0.0
1361           if((lb1bn1.and.lb2mn0)
1362      &       .or.(lb2bn1.and.lb1mn0).
1363      & or.(lb1bn0.and.lb2mn1).or.(lb2bn0.and.lb1mn1).
1364      & or.(lb1.eq.9.and.lb2mn2).or.(lb2.eq.9.and.lb1mn2))then
1365                 nchrg=1
1366 cbz3/2/99/song
1367 c                if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0
1368 c                if(lb1bn0.or.lb2bn0) brsig=0.5*sigma0
1369                 if(lb1bn1.or.lb2bn1) brsig=0.5*sigma0
1370                 if(lb1bn0.or.lb2bn0) brsig=2.0*sigma0
1371 cbz3/2/99/song end
1372 c                if(lb1.eq.9.or.lb2.eq.9) brsig=1.5*sigma0
1373           endif
1374           if( (lb1bn0.and.lb2mn0 )
1375      &       .or.(lb2bn0.and.lb1mn0)
1376      &  .or.(lb1bn1.and.lb2mn2).or.(lb2bn1.and.lb1mn2)
1377      &  .or.(lb1.eq.6.and.lb2mn1).or.(lb2.eq.6.and.lb1mn1)) then
1378                 nchrg=0
1379                 if(lb1bn1.or.lb2bn1) then
1380 cbz3/2/99/song
1381 c                  brsig=1.5*sigma0
1382                   brsig=3.0*sigma0
1383 cbz3/2/99/song end
1384 cbz3/11/99/song
1385 c                  ratiok = 1./3.
1386                   ratiok = 2./3.
1387 cbz3/11/99/song end
1388
1389 c                  ratiok: the ratio of channels: ->nK+k- vs. -> pK0K-
1390                 endif
1391                 if(lb1bn0.or.lb2bn0) then
1392                   brsig=2.5*sigma0
1393 cbz3/2/99/song
1394 c                  ratiok = 0.8
1395                   ratiok = 0.2
1396 cbz3/2/99/song end
1397                 endif
1398 c                if(lb1.eq.6.or.lb2.eq.6) then
1399 c     lb=6 : D-
1400 c                  brsig=1.5*sigma0
1401 c                  ratiok = 0.5
1402 c                endif
1403           endif
1404           if( (lb1bn0.and.lb2mn2)
1405      &       .or.(lb2bn0.and.lb1mn2)
1406      & .or.(lb1.eq.6.and.lb2mn0).or.(lb2.eq.6.and.lb1mn0)) then
1407                 nchrg=-1
1408                 if(lb1bn0.or.lb2bn0) brsig=sigma0
1409 c                if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0
1410           endif
1411 c          if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then
1412 c                nchrg=-2
1413 c          endif
1414 c          if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1)
1415 c    &           .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1416 c                nchrg=2
1417 c          endif
1418
1419 cbz3/11/99 neutralk
1420           if((lb1.eq.6.and.lb2mn2)
1421      &       .or.(lb2.eq.6.and.lb1mn2))then
1422                 nchrg=-2
1423           endif
1424 cbz3/11/99 neutralk
1425 cbz3/8/99 neutralk
1426           if((lb1bn1.and.lb2mn1)
1427      &       .or.(lb2bn1.and.lb1mn1)
1428      & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1429                 nchrg=2
1430           endif
1431 cbz3/8/99 neutralk end
1432
1433 cbz3/7/99 neutralk
1434           IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN
1435              BRSIG = 3.0 * SIGMA0
1436           END IF
1437 cbz3/7/99 neutralk end
1438
1439         endif
1440
1441 c        3. consider K- + N(N*,D) absorption.
1442 c        if((lb1bn.and.lb2.eq.21).OR.(lb2bn.and.lb1.eq.21)) then
1443         if( (lb1bn.and.(lb2.eq.21.or.lb2.eq.-30)).OR.
1444      &     (lb2bn.and.(lb1.eq.21.or.lb1.eq.-30)) )then 
1445 c          bmass=em1+em2-aka
1446           bmass=0.938
1447           if(srt.le.(bmass+aka)) then
1448 cbz3/2/99
1449 c                write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt
1450 cbz3/2/99end
1451                 pkaon=0.
1452           else
1453             pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
1454           endif
1455           sig=0.
1456           if(lb1.eq.1.or.lb2.eq.1.or.lb1.eq.8.or.lb2.eq.8.or.
1457      &    lb1.eq.11.or.lb2.eq.11.or.lb1.eq.13.or.lb2.eq.13) then
1458 c          K- + (D+,N*+)p ->
1459               nchrg=0
1460               sigela=akPel(pkaon)
1461               sigsgm=3.*akPsgm(pkaon)
1462               sig=sigela+sigsgm+akPlam(pkaon)
1463           endif
1464           if(lb1.eq.2.or.lb2.eq.2.or.lb1.eq.7.or.lb2.eq.7.or.
1465      &    lb1.eq.10.or.lb2.eq.10.or.lb1.eq.12.or.lb2.eq.12) then
1466 c          K- + (D0, N*0)n ->
1467               nchrg=-1
1468               sigela=akNel(pkaon)
1469               sigsgm=2.*akNsgm(pkaon)
1470               sig=sigela+sigsgm+akNlam(pkaon)
1471           endif
1472           if(lb1.eq.6.or.lb2.eq.6) then
1473 c     K- + D-
1474               nchrg=-2
1475               sigela=akNel(pkaon)
1476               sigsgm=akNsgm(pkaon)
1477               sig=sigela+sigsgm
1478           endif
1479           if(lb1.eq.9.or.lb2.eq.9) then
1480 c     K- + D++
1481               nchrg=1
1482               sigela=akPel(pkaon)
1483               sigsgm=2.*akPsgm(pkaon)
1484               sig=sigela+sigsgm+akPlam(pkaon)
1485           endif
1486
1487 cbz3/8/99 neutralk
1488           sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
1489           SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1490           SIG = sigela + SIGSGM + AKPLAM(PKAON)
1491 cbz3/8/99 neutralk end
1492
1493           if(sig.gt.1.e-7) then
1494 c     K(-) + N reactions
1495               icase=3
1496               brel=sigela/sig
1497               brsgm=sigsgm/sig
1498 c              branch_lambda=akNlam(pkaon)/sig
1499               brsig = sig
1500           endif
1501         endif
1502
1503 c        4. meson + hyperon -> K- + N
1504 c        if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR.
1505 c     &     ((lb2.ge.14.and.lb2.le.17).and.lb1mn)) then
1506         if(((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.3.and.lb2.le.5)).OR.
1507      &     ((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.3.and.lb1.le.5)))then
1508 c        first classify the reactions due to total charge.
1509            nchrg=-100
1510            if((lb1.eq.15.and.(lb2.eq.3.or.lb2.eq.25)).OR.
1511      &              (lb2.eq.15.and.(lb1.eq.3.or.lb1.eq.25))) then
1512                 nchrg=-2
1513 c     D-
1514                   bmass=1.232
1515            endif
1516            if((lb1.eq.15.and.lb2mn0).or.(lb2.eq.15.and.lb1mn0).OR.
1517      &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.3.or.lb2.eq.25)).OR.
1518      &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.3.or.lb1.eq.25)))then
1519                 nchrg=-1
1520 c     n
1521                  bmass=0.938
1522            endif
1523            if((lb1.eq.15.and.(lb2.eq.5.or.lb2.eq.27)).OR.
1524      &              (lb2.eq.15.and.(lb1.eq.5.or.lb1.eq.27)).or.
1525      &        (lb1.eq.17.and.(lb2.eq.3.or.lb2.eq.25)).OR.
1526      &              (lb2.eq.17.and.(lb1.eq.3.or.lb1.eq.25)).or.
1527      &       ((lb1.eq.14.or.lb1.eq.16).and.lb2mn0).OR.
1528      &       ((lb2.eq.14.or.lb2.eq.16).and.lb1mn0)) then
1529                 nchrg=0
1530 c     p
1531                  bmass=0.938
1532            endif
1533            if((lb1.eq.17.and.lb2mn0).or.(lb2.eq.17.and.lb1mn0).OR.
1534      &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.5.or.lb2.eq.27)).OR.
1535      &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.5.or.lb1.eq.27)))then
1536                 nchrg=1
1537 c     D++
1538                  bmass=1.232
1539            endif
1540            sig = 0.
1541            if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then
1542 c     PI+sigma or PI + Lambda => Kbar + N reactions
1543              icase=4
1544 c             pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
1545              pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
1546 c     lambda + Pi
1547              if(lb1.eq.14.or.lb2.eq.14) then
1548                 if(nchrg.ge.0) sigma0=akPlam(pkaon)
1549                 if(nchrg.lt.0) sigma0=akNlam(pkaon)
1550 c     sigma + pi
1551              else
1552 c     K-p or K-D++
1553                 if(nchrg.ge.0) sigma0=akPsgm(pkaon)
1554 c     K-n or K-D-
1555                 if(nchrg.lt.0) sigma0=akNsgm(pkaon)
1556
1557 cbz3/8/99 neutralk
1558                 SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1559 cbz3/8/99 neutralk end
1560
1561              endif
1562              sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
1563      &         (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
1564 cbz3/8/99 neutralk
1565 c     if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D-
1566 c     K0barD++, K-D-
1567              if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
1568
1569 cbz3/8/99 neutralk end
1570
1571 c             the factor 2 comes from spin of delta, which is 3/2
1572 c             detailed balance. copy from Page 423 of N.P. A614 1997
1573
1574 cbz3/8/99 neutralk
1575              IF (LB1 .EQ. 14 .OR. LB2 .EQ. 14) THEN
1576                 SIG = 4.0 / 3.0 * SIG
1577              ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN
1578                 SIG = 8.0 / 9.0 * SIG
1579              ELSE
1580                 SIG = 4.0 / 9.0 * SIG
1581              END IF
1582 cbz3/8/99 neutralk end
1583              brsig = sig
1584              if(sig.lt.1.e-7) sig = 1.e-7
1585            endif
1586 csp05/07/01
1587 * comment icase=4 statement below if only inelastic
1588 c     PI+L/Si => Kbar + N  OR ELASTIC SCATTERING
1589            icase=4
1590            brsig = sig
1591 c     elastic xsecn of 10mb
1592            sigela = 10.
1593            sig = sig + sigela
1594            brel = sigela/sig
1595 cc          brsig = sig
1596 csp05/07/01 end   
1597         endif
1598 c
1599 c        if(em2.lt.0.2.and.em1.lt.0.2) then
1600 c     PI + PI 
1601 c             icase=5
1602 c     assumed PI PI total x section.
1603 c              sig=50.
1604 c     Mk + Mkbar
1605 c              s0=aka+aka
1606 c              brsig = 0.
1607 c              if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76
1608 c              x section for PIPI->KKbar   PRC43 (1991) 1881
1609 c        endif
1610         if(icase.eq.-1) then
1611            ictrl = -1
1612            return
1613         endif
1614         px1cm=pcx
1615         py1cm=pcy
1616         pz1cm=pcz
1617         ds=sqrt(sig/31.4)
1618         dsr=ds+0.1
1619         ec=(em1+em2+0.02)**2
1620 c        ec=3.59709
1621 c        if((e(i1).ge.1.).and.(e(i2).ge.1.)) ec = 4.75
1622
1623         call distce(i1,i2,dsr,ds,dt,ec,srt,ic,px1cm,py1cm,pz1cm)
1624         if(ic.eq.-1) then
1625 c     no anti-kaon production
1626            ictrl = -1
1627 c           in=in+1
1628 c           write(60,*)'--------------distance-----',in
1629            return
1630         endif
1631
1632 clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi, 
1633 c     sgsum,sgsum1,sgsum3:
1634         ik=0
1635         ik0=0
1636         ik1=0
1637         ik2=0
1638         ik3=0
1639         il=0
1640         im=0
1641         im3=0
1642         im4=0
1643         in=0
1644         inpion=0
1645         ipipi=0
1646         sgsum=0.
1647         sgsum1=0.
1648         sgsum3=0.
1649         if(icase.eq.1) then
1650            ik=ik+1
1651            if(srt.gt.2.8639) then
1652                 ik0=ik0+1
1653                 if(em1.lt.1.0.and.em2.lt.1.0) then
1654                         ik1=ik1+1
1655                         sgsum1=sgsum1+brsig
1656 c                        ratio_1=sgsum1/ik1/40.
1657                 endif
1658                 if(em1.gt.1.0.and.em2.gt.1.0) then
1659                         ik3=ik3+1
1660                         sgsum3=sgsum3+brsig
1661 c                        ratio_3=sgsum3/ik3/40.
1662                 endif
1663                 if(em1.gt.1.0.and.em2.lt.1.0) ik2=ik2+1
1664                 if(em1.lt.1.0.and.em2.gt.1.0) ik2=ik2+1
1665                 sgsum=sgsum+brsig
1666 c                ratio=sgsum/ik0/40.
1667            endif
1668         endif
1669         if(icase.eq.2) inpion=inpion+1
1670         if(icase.eq.5) ipipi=ipipi+1
1671 c        write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio
1672 c        write(62,*)'inpion,ipipi',inpion,ipipi
1673         if(RANART(NSEED).gt.(brsig/sig)) then
1674 c     no anti-kaon production
1675            ictrl = -1
1676            return
1677         endif
1678         il=il+1
1679 c        kaons could be created now.
1680         if(icase.eq.1) then
1681           in=in+1
1682 c          write(60,*)'------in,s2kaon,sig=',in,brsig,sig,lb1,lb2
1683           call nnkaon(irun,iseed,
1684      &          ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1685         endif
1686         if(icase.eq.2) then
1687           im=im+1
1688 c          call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1689 c     &              pcx,pcy,pcz,nchrg,ratiok)
1690           call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1691      &              pcx,pcy,pcz,nchrg,ratiok,iblock)
1692         endif
1693 c
1694         if(icase.eq.3) then
1695           im3=im3+1
1696 c          write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon
1697 c          write(63,*)'sig,el,sigma',sig,brel,brsgm
1698 c          write(63,*)'srt,pcx,pcy,pcz,em1,em2',srt,pcx,pcy,pcz,em1,em2
1699           call kaonN(brel,brsgm,irun,iseed,dt,nt,ictrl,
1700      &                i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1701 c         this subroutine format is diff. since three final states are possible
1702         endif
1703 c
1704
1705         if(icase.eq.4) then
1706           im4=im4+1
1707 c          write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig
1708 c          write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon
1709
1710 csp06/07/01
1711       if(RANART(NSEED).lt.brel) then
1712          ielstc = 1
1713       else
1714          ielstc = 0
1715       endif                  
1716 c          call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1717 c     &                   pcx,pcy,pcz,nchrg)
1718           call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1719      &                   pcx,pcy,pcz,nchrg,iblock)
1720
1721 csp06/07/01 end
1722         endif
1723 c        if(icase.eq.5) then
1724 c          im5=im5+1
1725 c          write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig
1726 c          call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz)
1727 c        endif
1728 cbz3/2/99
1729 c        write(101,*)lb1,lb2,lb(i1),lb(i2)
1730 c        write(101,*)em1,em2,e(i1),e(i2),srt
1731 cbz3/2/99end
1732
1733         return
1734         end
1735
1736 ******************************************
1737 * for pp-->pp + kaon + anti-kaon
1738 c      real*4 function X2kaon(srt)
1739       real function X2kaon(srt)
1740       SAVE   
1741 *  This function contains the experimental total pp->pp+K(+)K(-) Xsections    *
1742 *  srt    = DSQRT(s) in GeV                                                   *
1743 *  xsec   = production cross section in mb                                    *
1744 *                                                                             *
1745 ******************************************
1746 c     minimum c.m.s. energy to create 2 kaon. = 2*(mp+mk)        
1747         smin = 2.8639
1748         x2kaon=0.0000001
1749         if(srt.lt.smin)return
1750         sigma1 = 2.8
1751         sigma2 = 7.7
1752         sigma3 = 3.9
1753         x = srt**2/smin**2 + 0.0000001
1754         f1 = (1.+1./sqrt(x))*alog(x) - 4.*(1.-1./sqrt(x))
1755         f2 = 1. - (1./sqrt(x))*(1.+alog(sqrt(x)))
1756         f3 = ((x-1.)/x**2)**3.5
1757         x2kaon = (1.-1./x)**3*(sigma1*f1 + sigma2*f2) + sigma3*f3
1758         return
1759         END
1760
1761         real function piNsg0(srt)
1762       SAVE   
1763 * cross section in mb for PI- + P -> P + K0 + K-
1764 c     Mn + 2* Mk
1765         srt0 = 0.938 + 2.*0.498
1766         if(srt.lt.srt0) then
1767            piNsg0 = 0.0
1768            return
1769         endif
1770         ratio = srt0**2/srt**2
1771         piNsg0=1.121*(1.-ratio)**1.86*ratio**2
1772         return
1773         end
1774
1775         real function akNel(pkaon)
1776       SAVE   
1777 *cross section in mb for K- + N reactions.
1778 c        the following data come from PRC 41 (1701)
1779 c        sigma1: K(-) + neutron elastic
1780         if(pkaon.lt.0.5.or. pkaon.ge.4.0) sigma1=0.
1781         if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma1=20.*pkaon**2.74
1782         if(pkaon.ge.1.0.and.pkaon.lt.4.0) sigma1=20.*pkaon**(-1.8)
1783         akNel=sigma1
1784         return
1785         end
1786
1787         real function akPel(pkaon)
1788       SAVE   
1789 *cross section in mb for K- + N reactions.
1790 c        the following data come from PRC 41 (1701)
1791 c        sigma2: K(-) + proton elastic
1792         if(pkaon.lt.0.25.or. pkaon.ge.4.0) sigma2=0.
1793         if(pkaon.ge.0.25.and.pkaon.lt.4.0) sigma2=13.*pkaon**(-0.9)
1794         akPel=sigma2
1795         return
1796         end
1797
1798         real function akNsgm(pkaon)
1799       SAVE   
1800 *cross section in mb for K- + N reactions.
1801 c        sigma2: x section for K- + n -> sigma0 + PI-
1802         if(pkaon.lt.0.5.or. pkaon.ge.6.0) sigma2=0.
1803         if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma2=1.2*pkaon**(-1.3)
1804         if(pkaon.ge.1.0.and.pkaon.lt.6.0) sigma2=1.2*pkaon**(-2.3)
1805         akNsgm=sigma2
1806         return
1807         end
1808
1809         real function akPsgm(pkaon)
1810       SAVE   
1811 *cross section in mb for K- + N reactions.
1812 c        sigma1: x section for K- + p -> sigma0 + PI0
1813         if(pkaon.lt.0.2.or. pkaon.ge.1.5) sigma1=0.
1814         if(pkaon.ge.0.2.and.pkaon.lt.1.5) sigma1=0.6*pkaon**(-1.8)
1815         akPsgm=sigma1
1816         return
1817         end
1818
1819         real function akPlam(pkaon)
1820       SAVE   
1821 *cross section in mb for K- + N reactions.
1822 c        sigma: x section for K- + p -> lambda + PI0
1823         p=pkaon
1824         if(pkaon.lt.0.2.or. pkaon.ge.10.0) sigma=0.
1825         if(pkaon.ge.0.2.and.pkaon.lt.0.9) sigma=50.*p**2-67.*p+24.
1826         if(pkaon.ge.0.9.and.pkaon.lt.10.0) sigma=3.0*pkaon**(-2.6)
1827         akPlam=sigma
1828         return
1829         end
1830
1831         real function akNlam(pkaon)
1832       SAVE   
1833 *cross section in mb for K- + N reactions.
1834         akNlam=akPlam(pkaon)
1835         return
1836         end
1837
1838 * GQ Li parametrization (without resonance)
1839         real function akNPsg(pkaon)
1840       SAVE   
1841 *cross section in mb for K- + N reactions.
1842 c       sigma1: x section for K- + p/n -> sigma0 + PI0
1843          if(pkaon.le.0.345)then
1844            sigma1=0.624*pkaon**(-1.83)
1845          else
1846            sigma1=0.7*pkaon**(-2.09)
1847          endif
1848         akNPsg=sigma1
1849         return
1850         end   
1851
1852 c-----------------------------------------------------------------------
1853
1854 c.....extracted from G. Song's ART expasion including K- interactions
1855 c.....file `NEWNNK.FOR'
1856
1857         subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock,
1858      &                                   srt,pcx,pcy,pcz,nchrg)
1859 c        <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997
1860 c     CANCELED also alpha=1 changed to alpha=3 to decrease the leadng effect.
1861       PARAMETER      (MAXSTR=150001,MAXR=1)
1862       PARAMETER      (AKA=0.498)
1863       COMMON   /AA/  R(3,MAXSTR)
1864 cc      SAVE /AA/
1865       COMMON   /BB/  P(3,MAXSTR)
1866 cc      SAVE /BB/
1867       COMMON   /CC/  E(MAXSTR)
1868 cc      SAVE /CC/
1869       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
1870 cc      SAVE /EE/
1871       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
1872 cc      SAVE /BG/
1873       COMMON   /NN/NNN
1874 cc      SAVE /NN/
1875       COMMON   /RUN/NUM
1876 cc      SAVE /RUN/
1877       COMMON   /PA/RPION(3,MAXSTR,MAXR)
1878 cc      SAVE /PA/
1879       COMMON   /PB/PPION(3,MAXSTR,MAXR)
1880 cc      SAVE /PB/
1881       COMMON   /PC/EPION(MAXSTR,MAXR)
1882 cc      SAVE /PC/
1883       COMMON   /PD/LPION(MAXSTR,MAXR)
1884 cc      SAVE /PD/
1885       dimension px(4),py(4),pz(4)
1886       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1887      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1888      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1889       SAVE   
1890 c      dm1=e(i1)
1891 c      dm2=e(i2)
1892       dm3=0.938
1893       dm4=0.938
1894 c     10/24/02 initialize n to 0:
1895       n=0
1896
1897 cbz3/11/99 neutralk
1898 c        if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232
1899 c        if(nchrg.eq.4) dm4=1.232
1900         if(nchrg.le.-1.or.nchrg.ge.3) dm3=1.232
1901         if(nchrg.eq.-2.or.nchrg.eq.4) dm4=1.232
1902 cbz3/11/99 neutralk end
1903           iblock = 0 
1904         call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
1905         if(iflag.lt.0) then
1906 c           write(60,*)'------------final state fail-------',n
1907 c     no anti-kaon production
1908            ictrl = -1
1909            n=n+1
1910            return
1911         endif
1912         iblock = 12
1913 * Rotate the momenta of particles in the cms of I1 & I2
1914 * px(1), py(1), pz(1): momentum of I1
1915 * px(2), py(2), pz(2): momentum of I2
1916 * px(3), py(3), pz(3): momentum of anti-kaon
1917 * px(4), py(4), pz(4): momentum of kaon
1918
1919
1920 c     10/28/02 get rid of argument usage mismatch in rotate():
1921         pxrota=px(1)
1922         pyrota=py(1)
1923         pzrota=pz(1)
1924 c        call rotate(pcx,pcy,pcz,px(1),py(1),pz(1))
1925         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1926         px(1)=pxrota
1927         py(1)=pyrota
1928         pz(1)=pzrota
1929 c
1930         pxrota=px(2)
1931         pyrota=py(2)
1932         pzrota=pz(2)
1933 c        call rotate(pcx,pcy,pcz,px(2),py(2),pz(2))
1934         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1935         px(2)=pxrota
1936         py(2)=pyrota
1937         pz(2)=pzrota
1938 c
1939         pxrota=px(3)
1940         pyrota=py(3)
1941         pzrota=pz(3)
1942 c        call rotate(pcx,pcy,pcz,px(3),py(3),pz(3))
1943         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1944         px(3)=pxrota
1945         py(3)=pyrota
1946         pz(3)=pzrota
1947 c
1948         pxrota=px(4)
1949         pyrota=py(4)
1950         pzrota=pz(4)
1951 c        call rotate(pcx,pcy,pcz,px(4),py(4),pz(4))
1952         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1953         px(4)=pxrota
1954         py(4)=pyrota
1955         pz(4)=pzrota
1956
1957         nnn=nnn+2
1958 c     K+
1959         lpion(nnn,irun)=23
1960         if(nchrg.eq.-1.or.nchrg.eq.-2) then
1961 c        To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K-
1962
1963 cbz3/7/99 neutralk
1964 c           lpion(nnn,irun)=24 ! K0
1965 cbz3/7/99 neutralk end
1966
1967         endif
1968 c     aka: rest mass of K
1969         epion(nnn,irun)=aka
1970 c     K-
1971         lpion(nnn-1,irun)=21
1972 c     aka: rest mass of K
1973         epion(nnn-1,irun)=aka
1974 * Find the momenta of particles in the final state in the nucleus_nucleus
1975 * cms frame.   Lorentz transformation into lab frame.
1976         e1cm   = sqrt(dm3**2 + px(1)**2 + py(1)**2 + pz(1)**2)
1977         p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz
1978         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
1979         pt1i1 = betax*transf + px(1)
1980         pt2i1 = betay*transf + py(1)
1981         pt3i1 = betaz*transf + pz(1)
1982         eti1  = dm3
1983 c        lb1   = lb(i1)
1984         lb1   = 2
1985         if(nchrg.ge.-2.and.nchrg.le.1) lb1=2
1986
1987 cbz3/7/99 neutralk
1988         if (nchrg .eq. -2 .or. nchrg .eq. -1) then
1989            lb1 = 6
1990         end if
1991 cbz3/7/99 neutralk end
1992
1993 cbz3/11/99 neutralk
1994 c        if(nchrg.eq.2.or.nchrg.eq.3) lb1=1
1995 c        if(nchrg.eq.4) lb1=9
1996         if(nchrg.eq.1.or.nchrg.eq.2) lb1=1
1997         if(nchrg.eq.3.or.nchrg.eq.4) lb1=9
1998 cbz3/11/99 neutralk end
1999
2000 * For second nulceon, same
2001         e2cm   = sqrt(dm4**2 + px(2)**2 + py(2)**2 + pz(2)**2)
2002         p2beta = px(2)*betax + py(2)*betay + pz(2)*betaz
2003         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2004         pt1i2 = betax*transf + px(2)
2005         pt2i2 = betay*transf + py(2)
2006         pt3i2 = betaz*transf + pz(2)
2007         eti2  = dm4
2008 c        lb2   = lb(i2)
2009         lb2   = 2
2010
2011 cbz3/11/99 neutralk
2012 c        if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2
2013 c        if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1
2014 c        if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9
2015 c        if(nchrg.eq.-2) lb2=6
2016         if(nchrg.ge.-1.or.nchrg.le.1) lb2=2
2017         if(nchrg.eq. 2.or.nchrg.eq.3) lb2=1
2018         if(nchrg.eq. 4) lb2=9
2019         if(nchrg.eq.-2) lb2=6
2020 cbz3/11/99 neutralk end
2021
2022 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2023                 p(1,i1)=pt1i1
2024                 p(2,i1)=pt2i1
2025                 p(3,i1)=pt3i1
2026                 e(i1)=eti1
2027                 lb(i1)=lb1
2028                 p(1,i2)=pt1i2
2029                 p(2,i2)=pt2i2
2030                 p(3,i2)=pt3i2
2031                 e(i2)=eti2
2032                 lb(i2)=lb2
2033
2034 c                px1 = p(1,i1)
2035 c                py1 = p(2,i1)
2036 c                pz1 = p(3,i1)
2037 c                em1 = e(i1)
2038 c                id(i1) = 2
2039 c                id(i2) = 2
2040 c                id1 = id(i1)
2041 c                iblock = 101  ! K(+)K(-) production
2042 * Get anti-kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2043         epcmk = sqrt(epion(nnn-1,irun)**2 + px(3)**2+py(3)**2+pz(3)**2)
2044         betak = px(3)*betax + py(3)*betay + pz(3)*betaz
2045         transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2046         ppion(1,nnn-1,irun)=betax*transf + px(3)
2047         ppion(2,nnn-1,irun)=betay*transf + py(3)
2048         ppion(3,nnn-1,irun)=betaz*transf + pz(3)
2049         rpion(1,nnn-1,irun)=r(1,i1)
2050         rpion(2,nnn-1,irun)=r(2,i1)
2051         rpion(3,nnn-1,irun)=r(3,i1)
2052 clin-5/2008:
2053         dppion(nnn-1,irun)=dpertp(i1)*dpertp(i2)
2054 * Same thing for kaon **************************************
2055         epcmak = sqrt(epion(nnn,irun)**2 + px(4)**2 +py(4)**2+pz(4)**2)
2056         betaak = px(4)*betax + py(4)*betay + pz(4)*betaz
2057         transf= gamma*(gamma*betaak/(gamma+1.) + epcmak)
2058         ppion(1,nnn,irun)=betax*transf + px(4)
2059         ppion(2,nnn,irun)=betay*transf + py(4)
2060         ppion(3,nnn,irun)=betaz*transf + pz(4)
2061         rpion(1,nnn,irun)=r(1,i2)
2062         rpion(2,nnn,irun)=r(2,i2)
2063         rpion(3,nnn,irun)=r(3,i2)
2064 clin-5/2008:
2065         dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2066         return
2067         end
2068
2069         subroutine lorntz(ilo,b,pi,pj)
2070 c       It uses to perform Lorentz (or inverse Lorentz) transformation
2071         dimension pi(4),pj(4),b(3)
2072       SAVE   
2073 c       dimension db(3)
2074         bb=b(1)*b(1)+b(2)*b(2)+b(3)*b(3)
2075         deno3=sqrt(1.-bb)
2076         if(deno3.eq.0.)deno3=1.e-10
2077         gam=1./deno3
2078         ga=gam*gam/(gam+1.)
2079         if(ilo.eq.1) goto 100
2080 c       Lorentz transformation
2081         pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2082         pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2083 c       drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3)
2084 c       drdb=db(1)*b(1)+db(2)*b(2)+db(3)*b(3)
2085         do 1001 i=1,3
2086            pi(i)=pi(i)+b(i)*(ga*pib-gam*pi(4))
2087            pj(i)=pj(i)+b(i)*(ga*pjb-gam*pj(4))
2088 c       drd(i)=drd(i)+b(i)*ga*drb
2089 c       db(i)=db(i)+b(i)*ga*drdb
2090  1001   continue
2091         pi(4)=gam*(pi(4)-pib)
2092         pj(4)=gam*(pj(4)-pjb)
2093         return
2094 100     continue
2095 c       inverse Lorentz transformation
2096         pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2097         pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2098         do 1002 i=1,3
2099            pi(i)=pi(i)+b(i)*(ga*pib+gam*pi(4))
2100            pj(i)=pj(i)+b(i)*(ga*pjb+gam*pj(4))
2101  1002   continue
2102         pi(4)=gam*(pi(4)+pib)
2103         pj(4)=gam*(pj(4)+pjb)
2104         return
2105         end
2106         
2107         subroutine fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
2108 *        function: decide final momentum for N,N,K(+),and K(-)        
2109         dimension px(4), py(4), pz(4), pe(4)
2110         COMMON/RNDF77/NSEED
2111 cc      SAVE /RNDF77/
2112         SAVE   
2113
2114         iseed=iseed
2115         iflag=-1
2116 c        iflag=-1: fail to find momenta
2117 c             = 1: success
2118         pio=3.1415926
2119         aka=0.498
2120 c        v=0.43
2121 c        w=-0.84
2122 c        b=3.78
2123 c        c=0.47
2124 c        d=3.60
2125 c        fmax=1.056
2126 c        gmax=1.+c
2127
2128         icount=0
2129         ekmax=(srt-dm3-dm4)/2.
2130         if(ekmax.le.aka) return
2131         pkmax=sqrt(ekmax**2-aka**2)
2132
2133         if(dm3.le.0.0.or.dm4.le.0.0) then
2134            write(1,*)'error: minus mass!!!'
2135            return
2136         endif
2137
2138 c        after we have the momenta for both nucleus, we sample the
2139 c        transverse momentum for K-. 
2140 c        dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on
2141 c        page 72, fig 23i.
2142 50        continue
2143         icount=icount+1
2144         if(icount.gt.10) return
2145         ptkmi2=-1./4.145*alog(RANART(NSEED))
2146         ptkm=sqrt(ptkmi2)
2147 3        v1=RANART(NSEED)
2148         v2=RANART(NSEED)
2149         rsq=v1**2+v2**2
2150         if(rsq.ge.1.0.or.rsq.le.0.) goto 3
2151         fac=sqrt(-2.*alog(rsq)/rsq)
2152         guass=v1*fac
2153         if(guass.ge.5.) goto 3
2154         xstar=guass/5.
2155         pzkm=pkmax*xstar
2156         ekm=sqrt(aka**2+pzkm**2+ptkm**2)
2157         if(RANART(NSEED).gt.aka/ekm) goto 50
2158         bbb=RANART(NSEED)
2159         px(3)=ptkm*cos(2.*pio*bbb)
2160         py(3)=ptkm*sin(2.*pio*bbb)
2161         if(RANART(NSEED).gt.0.5) pzkm=-1.*pzkm
2162         pz(3)=pzkm
2163         pe(3)=ekm
2164 150        ptkpl2=-1./3.68*alog(RANART(NSEED))
2165         ptkp=sqrt(ptkpl2)
2166 13        v1=RANART(NSEED)
2167         v2=RANART(NSEED)
2168         rsq=v1**2+v2**2
2169         if(rsq.ge.1.0.or.rsq.le.0.) goto 13
2170         fac=sqrt(-2.*alog(rsq)/rsq)
2171         guass=v1*fac
2172         if(guass.ge.3.25) goto 13
2173         xstar=guass/3.25
2174         pzkp=pkmax*xstar
2175         ekp=sqrt(aka**2+pzkp**2+ptkp**2)
2176         if(RANART(NSEED).gt.aka/ekp) goto 150
2177         bbb=RANART(NSEED)
2178         px(4)=ptkp*cos(2.*pio*bbb)
2179         py(4)=ptkp*sin(2.*pio*bbb)
2180         if(RANART(NSEED).gt.0.5) pzkp=-1.*pzkp
2181         pz(4)=pzkp
2182         pe(4)=ekp
2183
2184         resten=srt-pe(3)-pe(4)
2185         restpz=-pz(3)-pz(4)
2186 c     resample
2187         if(resten.le.abs(restpz)) goto 50
2188         restms=sqrt(resten**2-restpz**2)
2189 c     resample 
2190         if(restms.lt.(dm3+dm4)) goto 50
2191         ptp2=-1./2.76*alog(RANART(NSEED))
2192         ptp=sqrt(ptp2)
2193         bbb=RANART(NSEED)
2194         px(2)=ptp*cos(2.*pio*bbb)
2195         py(2)=ptp*sin(2.*pio*bbb)
2196         px(1)=-1.*(px(4)+px(3)+px(2))
2197         py(1)=-1.*(py(4)+py(3)+py(2))
2198 c     transverse mass for K-
2199         rmt3=sqrt(dm3**2+px(1)**2+py(1)**2)
2200 c     transverse mass for K+
2201         rmt4=sqrt(dm4**2+px(2)**2+py(2)**2)
2202         if(restms.lt.(rmt3+rmt4)) goto 50
2203 c        else: sampling success!
2204         pzcms=sqrt((restms**2-(rmt3+rmt4)**2)*
2205      &             (restms**2-(rmt3-rmt4)**2))/2./restms
2206         if(RANART(NSEED).gt.0.5) then
2207            pz(1)=pzcms
2208            pz(2)=-pzcms
2209         else
2210            pz(1)=-pzcms
2211            pz(2)=pzcms
2212         endif
2213         beta=restpz/resten        
2214         gama=1./sqrt(1.-beta**2)
2215         pz(1)=pz(1)*gama + beta*gama*sqrt(rmt3**2+pz(1)**2)
2216         pz(2)=pz(2)*gama + beta*gama*sqrt(rmt4**2+pz(2)**2)
2217         pe(1)=sqrt(rmt3**2+pz(1)**2)
2218         pe(2)=sqrt(rmt4**2+pz(2)**2)
2219
2220         iflag=1
2221         return
2222         end
2223
2224 c-----------------------------------------------------------------------
2225
2226 c.....extracted from G. Song's ART expasion including K- interactions
2227 c.....file `NPIK.FOR'
2228
2229 ****************************************
2230 c        subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2231 c     &                  pcx,pcy,pcz,nchrg,ratiok)
2232         subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2233      &                  pcx,pcy,pcz,nchrg,ratiok,iblock)
2234 *
2235 * Process: PI + N -> K(-) + ANYTHING
2236 * 1.  PI- + P -> P + K0 + K-
2237 * 2.  PI+ + N -> P + K+ + K- 
2238 * 3.  PI0 + P -> P + K+ + K-
2239 * 4.  PI0 + N -> P + K0 + K-
2240 * 5.  PI0 + N -> N + K+ + K-
2241 * 6.  PI- + P -> N + K+ + K-
2242 * 7.  PI- + N -> N + K0 + K-
2243 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2244 ****************************************
2245       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2246       PARAMETER      (AKA=0.498)
2247       COMMON   /AA/  R(3,MAXSTR)
2248 cc      SAVE /AA/
2249       COMMON   /BB/  P(3,MAXSTR)
2250 cc      SAVE /BB/
2251       COMMON   /CC/  E(MAXSTR)
2252 cc      SAVE /CC/
2253       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2254 cc      SAVE /EE/
2255       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2256 cc      SAVE /BG/
2257       COMMON   /NN/NNN
2258 cc      SAVE /NN/
2259       COMMON   /RUN/NUM
2260 cc      SAVE /RUN/
2261       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2262 cc      SAVE /PA/
2263       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2264 cc      SAVE /PB/
2265       COMMON   /PC/EPION(MAXSTR,MAXR)
2266 cc      SAVE /PC/
2267       COMMON   /PD/LPION(MAXSTR,MAXR)
2268 cc      SAVE /PD/
2269       dimension bb(3),p1(4),p2(4),p3(4),px(4),py(4),pz(4)
2270       COMMON/RNDF77/NSEED
2271 cc      SAVE /RNDF77/
2272       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
2273      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
2274      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
2275       SAVE   
2276         iseed=iseed
2277         dt=dt
2278         nchrg=nchrg
2279         nt=nt
2280         ratiok=ratiok
2281         px(1)=px(1)
2282         py(1)=py(1)
2283         pz(1)=pz(1)
2284         px1cm=pcx
2285         py1cm=pcy
2286         pz1cm=pcz
2287         ictrl = 1
2288         lb1=lb(i1)
2289         lb2=lb(i2)
2290         k1=i1
2291         k2=i2
2292 c        k1 must be bayron. k2 be meson. If not, exchange.
2293         if(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13)) then
2294             k1=i2
2295             k2=i1
2296         endif
2297 cbz3/8/99 neutralk
2298 cbz10/12/99
2299 c        LB(I1) = 1 + 2 * RANART(NSEED)
2300 c        LB(I2) = 23
2301         LB(k1) = 1 + int(2*RANART(NSEED))
2302         LB(k2) = 23
2303 c       pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2))
2304 c     &           /2./srt
2305         pkmax=sqrt((srt**2-(aka+0.938+aka)**2)
2306      &           *(srt**2-(aka+0.938-aka)**2))/2./srt
2307         pk = RANART(NSEED)*pkmax
2308 c-----------------------------------------------------
2309         css=1.-2.*RANART(NSEED)
2310         sss=sqrt(1.-css**2)
2311         fai=2*3.1415926*RANART(NSEED)
2312         p3(1)=pk*sss*cos(fai)
2313         p3(2)=pk*sss*sin(fai)
2314         p3(3)=pk*css
2315         eip = srt - sqrt(aka**2 + pk**2)
2316         rmnp=sqrt(eip**2-pk**2)
2317         do 1001 i= 1, 3
2318            bb(i) = -1.*p3(i)/eip
2319  1001   continue
2320 c        bb: velocity of the other two particles as a whole.
2321         pznp=sqrt((rmnp**2-(aka+0.938)**2)
2322      c  *(rmnp**2-(0.938-aka)**2))/2./rmnp    
2323 c-----------------------------------------------------
2324         css=1.-2.*RANART(NSEED)
2325         sss=sqrt(1.-css**2)
2326         fai=2*3.1415926*RANART(NSEED)
2327         p1(1)=pznp*sss*cos(fai)
2328         p1(2)=pznp*sss*sin(fai)
2329         p1(3)=pznp*css
2330         p1(4)=sqrt(0.938**2+pznp**2)
2331         p2(4)=sqrt(aka**2+pznp**2)
2332         do 1002 i=1,3
2333            p2(i)=-1.*p1(i)
2334  1002   continue
2335 c        p1,p2: the momenta of the two particles in their cms
2336 c        p1: momentum of N or P
2337 c        p2: momentum of anti_kaon
2338 c        p3: momentum of K0 or K+
2339         ilo=1
2340 c        write(61,*)'--------p1,p2',p1,p2
2341 c        write(61,*)'--------bb',bb
2342         call lorntz(ilo,bb,p1,p2)
2343 c******* Checking *************
2344 c        pxsum = p1(1)+p2(1)+p3(1)
2345 c        pysum = p1(2)+p2(2)+p3(2)
2346 c        pzsum = p1(3)+p2(3)+p3(3)
2347 c        pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt
2348 c        write(61,*)'---p1,pxsum',p1,pxsum
2349 c        write(61,*)'---p2,pysum',p2,pysum
2350 c        write(61,*)'---p3,pzsum',p3,pzsum
2351 c        write(61,*)'---pesum',pesum
2352 c***********************************
2353
2354 * Rotate the momenta of particles in the cms of I1 & I2
2355 * px(1), py(1), pz(1): momentum of I1
2356 * px(2), py(2), pz(2): momentum of I2
2357 * px(3), py(3), pz(3): momentum of anti-kaon
2358
2359 c     10/28/02 get rid of argument usage mismatch in rotate():
2360         pxrota=p1(1)
2361         pyrota=p1(2)
2362         pzrota=p1(3)
2363 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2364         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2365         p1(1)=pxrota
2366         p1(2)=pyrota
2367         p1(3)=pzrota
2368 c
2369         pxrota=p2(1)
2370         pyrota=p2(2)
2371         pzrota=p2(3)
2372 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2373         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2374         p2(1)=pxrota
2375         p2(2)=pyrota
2376         p2(3)=pzrota
2377 c
2378         pxrota=p3(1)
2379         pyrota=p3(2)
2380         pzrota=p3(3)
2381 c        call rotate(pcx,pcy,pcz,p3(1),p3(2),p3(3))
2382         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2383         p3(1)=pxrota
2384         p3(2)=pyrota
2385         p3(3)=pzrota
2386
2387         nnn=nnn+1
2388 c     K(-)
2389         lpion(nnn,irun)=21
2390 c     aka: rest mass of K
2391         epion(nnn,irun)=aka
2392 * Find the momenta of particles in the final state in the nucleus_nucleus
2393 * cms frame.   Lorentz transformation into lab frame.
2394         e1cm   = sqrt(0.938**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2395         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2396         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2397         pt1i1 = betax*transf + p1(1)
2398         pt2i1 = betay*transf + p1(2)
2399         pt3i1 = betaz*transf + p1(3)
2400         eti1  = 0.938
2401         lb1   = lb(k1)
2402          
2403 * For second nulceon, same
2404         e2cm   = sqrt(aka**2 + p3(1)**2 + p3(2)**2 + p3(3)**2)
2405         p2beta = p3(1)*betax + p3(2)*betay + p3(3)*betaz
2406         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2407         pt1i2 = betax*transf + p3(1)
2408         pt2i2 = betay*transf + p3(2)
2409         pt3i2 = betaz*transf + p3(3)
2410         eti2  = aka
2411         lb2   = lb(k2)
2412
2413 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2414 *       k1 stand for nucleon, k2 stand for kaon. lpion stand for Kbar.
2415                 p(1,k1)=pt1i1
2416                 p(2,k1)=pt2i1
2417                 p(3,k1)=pt3i1
2418                 e(k1)=eti1
2419                 lb(k1)=lb1
2420                 p(1,k2)=pt1i2
2421                 p(2,k2)=pt2i2
2422                 p(3,k2)=pt3i2
2423                 e(k2)=eti2
2424                 lb(k2)=lb2
2425
2426 c                px1 = p(1,i1)
2427 c                py1 = p(2,i1)
2428 c                pz1 = p(3,i1)
2429 c                em1 = e(i1)
2430 c                id(i1) = 2
2431 c                id(i2) = 2
2432 c                id1 = id(i1)
2433 c     K(+)K(-) production
2434                 iblock = 101
2435 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2436 c  p2:  momentum of anti-kaon.
2437 c        epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2438         epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2+p2(2)**2+p2(3)**2)
2439         betak = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2440         transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2441         ppion(1,nnn,irun)=betax*transf + p2(1)
2442         ppion(2,nnn,irun)=betay*transf + p2(2)
2443         ppion(3,nnn,irun)=betaz*transf + p2(3)
2444 clin-5/2008:
2445         dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2446 cbz3/2/99
2447 c        write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun),
2448 c     &                    ppion(3,nnn,irun), dt*nt, srt
2449 cbz3/2/99end
2450 c        write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun),
2451 c     &                    ppion(3,nnn,irun), dt*nt, srt
2452         k=i2
2453         if(lb(i1).eq.1.or.lb(i1).eq.2) k=i1
2454         rpion(1,nnn,irun)=r(1,k)
2455         rpion(2,nnn,irun)=r(2,k)
2456         rpion(3,nnn,irun)=r(3,k)
2457         return
2458         end
2459
2460 c-----------------------------------------------------------------------
2461
2462 c.....extracted from G. Song's ART expasion including K- interactions
2463 c.....file `PIHYPN.FOR'
2464
2465 ******************************************
2466         subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,
2467      &     srt,pcx,pcy,pcz,nchrg,iblock)
2468 *
2469 * Process: PI + sigma(or Lambda) -> Kbar + N
2470 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2471 ******************************************
2472
2473 c NOTE: for PI + Hyperon: the produced kaons have mass 0.498
2474       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2475       PARAMETER      (AKA=0.498)
2476       COMMON   /AA/  R(3,MAXSTR)
2477 cc      SAVE /AA/
2478       COMMON   /BB/  P(3,MAXSTR)
2479 cc      SAVE /BB/
2480       COMMON   /CC/  E(MAXSTR)
2481 cc      SAVE /CC/
2482       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2483 cc      SAVE /EE/
2484       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2485 cc      SAVE /BG/
2486       COMMON   /NN/NNN
2487 cc      SAVE /NN/
2488       COMMON   /RUN/NUM
2489 cc      SAVE /RUN/
2490       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2491 cc      SAVE /PA/
2492       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2493 cc      SAVE /PB/
2494       COMMON   /PC/EPION(MAXSTR,MAXR)
2495 cc      SAVE /PC/
2496       COMMON   /PD/LPION(MAXSTR,MAXR)
2497 cc      SAVE /PD/
2498       dimension p1(4),p2(4)
2499       COMMON/RNDF77/NSEED
2500 cc      SAVE /RNDF77/
2501       SAVE   
2502         irun=irun
2503         iseed=iseed
2504         nt=nt
2505         dt=dt
2506         px1cm=pcx
2507         py1cm=pcy
2508         pz1cm=pcz
2509         ictrl = 1
2510 csp06/07/01
2511         if(ielstc .eq. 1) then
2512 *    L/Si + meson -> L/Si + meson
2513              k1=i1
2514              k2=i2
2515            dm3=e(k1)
2516            dm4=e(k2)
2517            iblock = 10
2518         else
2519            iblock = 12
2520 csp06/07/01 end  
2521 c        PI + Sigma(or Lambda) -> Kbar + N
2522         k1=i1
2523         k2=i2
2524 c        k1 must be bayron! So if I1 is PI, exchange k1 & k2.
2525         if(lb(i1).lt.14.or.lb(i1).gt.17) then
2526            k1=i2
2527            k2=i1
2528         endif
2529 cbz3/8/99 neutralk
2530         LB(K1) = 1 + int(2*RANART(NSEED))
2531         if(nchrg.eq.-2) lb(k1)=6
2532 c     if(nchrg.eq.-1) lb(k1)=2
2533 c     if(nchrg.eq. 0) lb(k1)=1
2534 c     if(nchrg.eq. 1) lb(k1)=9
2535         IF (NCHRG .EQ. 2) LB(K1) = 9
2536 cbz3/8/99 neutralk end
2537
2538 c     K-
2539         lb(k2)=21
2540         dm3=0.938
2541         if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232
2542         dm4=aka
2543 c        dm3,dm4: the mass of final state particles.
2544          endif
2545     
2546 ********Now, antikaon will be created.
2547 c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2548 c        pkmax: the maximum momentum of anti-kaon
2549         pkmax=sqrt((srt**2-(dm3+dm4)**2)*(srt**2-(dm3-dm4)**2))
2550      &         /2./srt
2551         pk=pkmax
2552 c-----------------------------------------------------
2553         css=1.-2.*RANART(NSEED)
2554         sss=sqrt(1.-css**2)
2555         fai=2*3.1415926*RANART(NSEED)
2556         p1(1)=pk*sss*cos(fai)
2557         p1(2)=pk*sss*sin(fai)
2558         p1(3)=pk*css
2559         do 1001 i=1,3
2560            p2(i)=-1.*p1(i)
2561  1001   continue
2562 c        p1,p2: the momenta of the two particles in their cms
2563 c        p1: momentum of kaon
2564 c        p2: momentum of Kbar
2565
2566 * Rotate the momenta of particles in the cms of I1 & I2
2567 clin-10/28/02 get rid of argument usage mismatch in rotate():
2568         pxrota=p1(1)
2569         pyrota=p1(2)
2570         pzrota=p1(3)
2571 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2572         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2573         p1(1)=pxrota
2574         p1(2)=pyrota
2575         p1(3)=pzrota
2576 c
2577         pxrota=p2(1)
2578         pyrota=p2(2)
2579         pzrota=p2(3)
2580 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2581         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2582         p2(1)=pxrota
2583         p2(2)=pyrota
2584         p2(3)=pzrota
2585 clin-10/28/02-end
2586
2587 * Find the momenta of particles in the final state in the nucleus_nucleus
2588 * cms frame.   Lorentz transformation into lab frame.
2589         e1cm   = sqrt(dm3**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2590         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2591         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2592         pt1i1 = betax*transf + p1(1)
2593         pt2i1 = betay*transf + p1(2)
2594         pt3i1 = betaz*transf + p1(3)
2595         eti1  = dm3
2596         lb1   = lb(k1)
2597          
2598 * For second kaon, same
2599         e2cm   = sqrt(dm4**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2600         p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2601         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2602         pt1i2 = betax*transf + p2(1)
2603         pt2i2 = betay*transf + p2(2)
2604         pt3i2 = betaz*transf + p2(3)
2605 cbz3/2/99
2606 c        write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt
2607 cbz3/2/99end
2608 c        write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt
2609         eti2  = dm4
2610         lb2   = lb(k2)
2611
2612 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2613 c        k1=i1
2614 c        k2=i2
2615 *       k1 stand for nucleon, k2 stand for kaon.
2616                 p(1,k1)=pt1i1
2617                 p(2,k1)=pt2i1
2618                 p(3,k1)=pt3i1
2619                 e(k1)=eti1
2620                 lb(k1)=lb1
2621                 p(1,k2)=pt1i2
2622                 p(2,k2)=pt2i2
2623                 p(3,k2)=pt3i2
2624                 e(k2)=eti2
2625                 lb(k2)=lb2
2626
2627 cc                iblock = 101  ! K(+)K(-) production
2628 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2629         return
2630         end
2631
2632 c-----------------------------------------------------------------------
2633
2634 c.....extracted from G. Song's ART expasion including K- interactions
2635 c.....file `KAONN.FOR'
2636
2637 ****************************************
2638         subroutine kaonN(brel,brsgm,irun,iseed,dt,nt,
2639      &     ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
2640 *
2641 * Process: PI + sigma(or Lambda) <- Kbar + N
2642 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2643 ****************************************
2644       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2645       PARAMETER      (AKA=0.498,ALA=1.1157,ASA=1.1974)
2646       COMMON   /AA/  R(3,MAXSTR)
2647 cc      SAVE /AA/
2648       COMMON   /BB/  P(3,MAXSTR)
2649 cc      SAVE /BB/
2650       COMMON   /CC/  E(MAXSTR)
2651 cc      SAVE /CC/
2652       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2653 cc      SAVE /EE/
2654       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2655 cc      SAVE /BG/
2656       COMMON   /NN/NNN
2657 cc      SAVE /NN/
2658       COMMON   /RUN/NUM
2659 cc      SAVE /RUN/
2660       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2661 cc      SAVE /PA/
2662       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2663 cc      SAVE /PB/
2664       COMMON   /PC/EPION(MAXSTR,MAXR)
2665 cc      SAVE /PC/
2666       COMMON   /PD/LPION(MAXSTR,MAXR)
2667 cc      SAVE /PD/
2668       dimension p1(4),p2(4)
2669       COMMON/RNDF77/NSEED
2670 cc      SAVE /RNDF77/
2671       SAVE   
2672         dt=dt
2673         irun=irun
2674         iseed=iseed
2675         nchrg=nchrg
2676         nt=nt
2677         px1cm=pcx
2678         py1cm=pcy
2679         pz1cm=pcz
2680         ictrl = 1
2681 c        ratio: used for isospin decision.
2682         k1=i1
2683         k2=i2
2684 c        k1 must be bayron! So if I1 is Kaon, exchange k1 & k2.
2685         if(e(i1).lt.0.5.and.e(i1).gt.0.01) then
2686            k1=i2
2687            k2=i1
2688         endif
2689 *** note: for print out only *******************************
2690 c     record kaon's mass
2691         eee=e(k2)
2692 *** end **************
2693         rrr=RANART(NSEED)
2694         if(rrr.lt.brel) then
2695 c       Kbar + N -> Kbar + N
2696            lb1=lb(k1)
2697            lb2=lb(k2)
2698            em1=e(k1)
2699            em2=e(k2)
2700            iblock = 10
2701         else 
2702            iblock = 12
2703         if(rrr.lt.(brel+brsgm)) then
2704 c        nchrg: Net charges of the two incoming particles.
2705 c           Kbar + N -> Sigma + PI
2706            em1=asa
2707            em2=0.138
2708
2709 cbz3/8/99 neutralk
2710            LB1 = 15 + int(3*RANART(NSEED))
2711            LB2 = 3 + int(3*RANART(NSEED))
2712         else
2713 c           Kbar + N -> Lambda + PI
2714            em1=ala
2715            em2=0.138
2716 c     LAmbda
2717            lb1=14
2718 cbz3/8/99 neutralk
2719            LB2 = 3 + int(3*RANART(NSEED))
2720 c           if(nchrg.eq.1)  lb2=5  ! K- + D++ -> Lambda + PI+
2721 c           if(nchrg.eq.0)  lb2=4  ! K- + p(D+,N*+) -> Lambda + PI0
2722 c          if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI-
2723 cbz3/8/99 neutralk
2724
2725         endif
2726         endif
2727         lb(k1)=lb1
2728         lb(k2)=lb2
2729     
2730 ********Now, antikaon will be created.
2731 c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2732 c        pkmax: the maximum momentum of anti-kaon
2733 c        write(63,*)'srt,em1,em2',srt,em1,em2
2734 c        write(63,*)'-srt,em1,em2',srt,em1,em2
2735         pkmax=sqrt((srt**2-(em1+em2)**2)*(srt**2-(em1-em2)**2))
2736      &         /2./srt
2737         pk=pkmax
2738 c-----------------------------------------------------
2739         css=1.-2.*RANART(NSEED)
2740         sss=sqrt(1.-css**2)
2741         fai=2*3.1415926*RANART(NSEED)
2742         p1(1)=pk*sss*cos(fai)
2743         p1(2)=pk*sss*sin(fai)
2744         p1(3)=pk*css
2745         do 1001 i=1,3
2746            p2(i)=-1.*p1(i)
2747  1001   continue
2748 c        p1,p2: the momenta of the two particles in their cms
2749 c        p1: momentum of kaon
2750 c        p2: momentum of Kbar
2751
2752 * Rotate the momenta of particles in the cms of I1 & I2
2753
2754 clin-10/28/02 get rid of argument usage mismatch in rotate():
2755         pxrota=p1(1)
2756         pyrota=p1(2)
2757         pzrota=p1(3)
2758 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2759         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2760         p1(1)=pxrota
2761         p1(2)=pyrota
2762         p1(3)=pzrota
2763 c
2764         pxrota=p2(1)
2765         pyrota=p2(2)
2766         pzrota=p2(3)
2767 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2768         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2769         p2(1)=pxrota
2770         p2(2)=pyrota
2771         p2(3)=pzrota
2772 clin-10/28/02-end
2773
2774 * Find the momenta of particles in the final state in the nucleus_nucleus
2775 * cms frame.   Lorentz transformation into lab frame.
2776         e1cm   = sqrt(em1**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2777         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2778         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2779         pt1i1 = betax*transf + p1(1)
2780         pt2i1 = betay*transf + p1(2)
2781         pt3i1 = betaz*transf + p1(3)
2782         eti1  = em1
2783          
2784 * For second kaon, same
2785         e2cm   = sqrt(em2**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2786         p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2787         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2788         pt1i2 = betax*transf + p2(1)
2789         pt2i2 = betay*transf + p2(2)
2790         pt3i2 = betaz*transf + p2(3)
2791         eti2  = em2
2792
2793 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2794 c        k1=i1
2795 c        k2=i2
2796 *       k1 stand for bayron, k2 stand for meson.
2797                 p(1,k1)=pt1i1
2798                 p(2,k1)=pt2i1
2799                 p(3,k1)=pt3i1
2800                 e(k1)=eti1
2801                 p(1,k2)=pt1i2
2802                 p(2,k2)=pt2i2
2803                 p(3,k2)=pt3i2
2804                 e(k2)=eti2
2805
2806 cc                iblock = 101  ! K(+)K(-) production
2807 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2808         return
2809         end
2810
2811 c=======================================================================
2812
2813 clin Below is the previous artana.f:
2814 c=======================================================================
2815
2816 c.....analysis subroutine before the hadronic space-time evolution
2817
2818       SUBROUTINE ARTAN1
2819       PARAMETER (MAXSTR=150001, MAXR=1)
2820 c.....y cut for mt spectrum
2821 cbz3/17/99
2822 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
2823       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
2824 cbz3/17/99 end
2825 c.....bin width for mt spectrum and y spectrum
2826 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2827 c      PARAMETER (BMT = 0.05, BY = 0.2)
2828       PARAMETER (BMT = 0.05, BY = 0.4)
2829       COMMON /RUN/ NUM
2830 cc      SAVE /RUN/
2831       COMMON /ARERC1/MULTI1(MAXR)
2832 cc      SAVE /ARERC1/
2833       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2834      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
2835      &     FT1(MAXSTR, MAXR),
2836      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2837      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2838 cbz3/17/99
2839 c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
2840 cc      SAVE /ARPRC1/
2841       COMMON /ARANA1/
2842      &     dy1ntb(50), dy1ntp(50), DY1HM(50), 
2843      &     DY1KP(50), DY1KM(50), DY1K0S(50),
2844      &     DY1LA(50), DY1LB(50), DY1PHI(50),
2845      &     dm1pip(50), dm1pim(50), DMT1PR(50),
2846      &     DMT1PB(50), DMT1KP(50), dm1km(50),
2847      &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
2848      &     dy1msn(50), DY1PIP(50), DY1PIM(50), 
2849      &     DY1PI0(50), DY1PR(50), DY1PB(50)
2850      &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
2851 cc      SAVE /ARANA1/
2852       SAVE   
2853
2854 cbz3/17/99 end
2855       DO 1002 J = 1, NUM
2856          DO 1001 I = 1, MULTI1(J)
2857             ITYP = ITYP1(I, J)
2858             PX = PX1(I, J)
2859             PY = PY1(I, J)
2860             PZ = PZ1(I, J)
2861             EE = EE1(I, J)
2862             XM = XM1(I, J)
2863 c     2/24/03 leptons and photons:
2864             if(xm.lt.0.01) goto 200
2865             ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
2866             eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
2867
2868             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
2869             IF (ABS(PZ) .GE. EE) THEN
2870                PRINT *, 'IN ARTAN1'
2871                PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
2872 cbzdbg2/16/99
2873                PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
2874 cbzdbg2/16/99
2875 cbzdbg2/15/99
2876                PRINT *, ' PZ = ', PZ, ' EE = ', EE
2877 cbzdbg2/16/99
2878                PRINT *, ' XM = ', XM
2879 cbzdbg2/16/99end
2880                GOTO 200
2881 c               STOP
2882 cbzdbg2/15/99end
2883             END IF
2884             DXMT = XMT - XM
2885             Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
2886 c.....rapidity cut for the rapidity distribution
2887             IF (ABS(Y) .GE. 10.0) GOTO 100
2888 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
2889 c            IY = 1 + int(ABS(Y) / BY)
2890 c            Ieta = 1 + int(ABS(eta) / BY)
2891             IF (ABS(eta) .GE. 10.0) GOTO 100
2892             IY = 1 + int((Y+10.) / BY)
2893             Ieta = 1 + int((eta+10.) / BY)
2894
2895             IF (ITYP .LT. -1000) THEN
2896                dy1ntb(IY) = dy1ntb(IY) - 1.0
2897             END IF
2898             IF (ITYP .GT. 1000) THEN
2899                dy1ntb(IY) = dy1ntb(IY) + 1.0
2900             END IF
2901             IF (ITYP .EQ. -2212) THEN
2902                dy1ntp(IY) = dy1ntp(IY) - 1.0
2903             END IF
2904             IF (ITYP .EQ. 2212) THEN
2905                dy1ntp(IY) = dy1ntp(IY) + 1.0
2906             END IF
2907 c            IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR.
2908 c     &         ITYP .EQ. -2212) THEN
2909             IF (ITYP .EQ. -2112) THEN
2910                DY1HM(IY) = DY1HM(IY) + 1.0
2911             END IF
2912 c
2913             IF (LUCHGE(ITYP).ne.0) THEN
2914                DY1CH(IY) = DY1CH(IY) + 1.0
2915                DE1CH(Ieta) = DE1CH(Ieta) + 1.0
2916                IF (LUCHGE(ITYP).lt.0) THEN
2917                   DY1NEG(IY) = DY1NEG(IY) + 1.0
2918                   DE1NEG(Ieta) = DE1NEG(Ieta) + 1.0
2919                endif
2920             END IF
2921
2922 cbz3/17/99
2923             IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. 
2924      &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
2925                dy1msn(IY) = dy1msn(IY) + 1.0
2926             END IF
2927             IF (ITYP .EQ. 211) THEN
2928                DY1PIP(IY) = DY1PIP(IY) + 1.0
2929             END IF
2930             IF (ITYP .EQ. -211) THEN
2931                DY1PIM(IY) = DY1PIM(IY) + 1.0
2932             END IF
2933             IF (ITYP .EQ. 111) THEN
2934                DY1PI0(IY) = DY1PI0(IY) + 1.0
2935             END IF
2936             IF (ITYP .EQ. 2212) THEN
2937                DY1PR(IY) = DY1PR(IY) + 1.0
2938             END IF
2939             IF (ITYP .EQ. -2212) THEN
2940                DY1PB(IY) = DY1PB(IY) + 1.0
2941             END IF
2942 cbz3/17/99 end
2943             IF (ITYP .EQ. 321) THEN
2944                DY1KP(IY) = DY1KP(IY) + 1.0
2945             END IF
2946             IF (ITYP .EQ. -321) THEN
2947                DY1KM(IY) = DY1KM(IY) + 1.0
2948             END IF
2949 clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S:
2950 c            IF (ITYP .EQ. 310) THEN
2951             IF (ITYP .EQ. 130) THEN
2952                DY1K0S(IY) = DY1K0S(IY) + 1.0
2953             END IF
2954             IF (ITYP .EQ. 3122) THEN
2955                DY1LA(IY) = DY1LA(IY) + 1.0
2956             END IF
2957             IF (ITYP .EQ. -3122) THEN
2958                DY1LB(IY) = DY1LB(IY) + 1.0
2959             END IF
2960             IF (ITYP .EQ. 333) THEN
2961                DY1PHI(IY) = DY1PHI(IY) + 1.0
2962             END IF
2963
2964 c.....insert rapidity cut for mt spectrum here
2965  100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
2966             IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
2967             IMT = 1 + int(DXMT / BMT)
2968             IF (ITYP .EQ. 211) THEN
2969                dm1pip(IMT) = dm1pip(IMT) + 1.0 / XMT
2970             END IF
2971             IF (ITYP .EQ. -211) THEN
2972                dm1pim(IMT) = dm1pim(IMT) + 
2973      &            1.0 / XMT
2974             END IF
2975             IF (ITYP .EQ. 2212) THEN
2976                DMT1PR(IMT) = DMT1PR(IMT) + 1.0 / XMT
2977             END IF
2978             IF (ITYP .EQ. -2212) THEN
2979                DMT1PB(IMT) = DMT1PB(IMT) + 1.0 / XMT
2980             END IF
2981             IF (ITYP .EQ. 321) THEN
2982                DMT1KP(IMT) = DMT1KP(IMT) + 1.0 / XMT
2983             END IF
2984             IF (ITYP .EQ. -321) THEN
2985                dm1km(IMT) = dm1km(IMT) + 1.0 / XMT
2986             END IF
2987 clin-4/24/03:
2988 c            IF (ITYP .EQ. 310) THEN
2989             IF (ITYP .EQ. 130) THEN
2990                dm1k0s(IMT) = dm1k0s(IMT) + 1.0 / XMT
2991             END IF
2992             IF (ITYP .EQ. 3122) THEN
2993                DMT1LA(IMT) = DMT1LA(IMT) + 1.0 / XMT
2994             END IF
2995             IF (ITYP .EQ. -3122) THEN
2996                DMT1LB(IMT) = DMT1LB(IMT) + 1.0 / XMT
2997             END IF
2998
2999  200        CONTINUE
3000  1001    CONTINUE
3001  1002 CONTINUE
3002
3003       RETURN
3004       END
3005
3006 c-----------------------------------------------------------------------
3007
3008 c.....analysis subroutine after the hadronic space-time evolution
3009
3010       SUBROUTINE ARTAN2
3011
3012       PARAMETER (MAXSTR=150001, MAXR=1)
3013 c.....y cut for mt spectrum
3014 cbz3/17/99
3015 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3016       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3017 cbz3/17/99 end
3018 c.....bin width for mt spectrum and y spectrum
3019 c      PARAMETER (BMT = 0.05, BY = 0.2)
3020       PARAMETER (BMT = 0.05, BY = 0.4)
3021       COMMON /RUN/ NUM
3022 cc      SAVE /RUN/
3023       COMMON /ARERC1/MULTI1(MAXR)
3024 cc      SAVE /ARERC1/
3025       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
3026      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
3027      &     FT1(MAXSTR, MAXR),
3028      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
3029      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
3030 cbz3/17/99
3031 c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
3032 cc      SAVE /ARPRC1/
3033       COMMON /ARANA2/
3034      &     dy2ntb(50), dy2ntp(50), DY2HM(50), 
3035      &     DY2KP(50), DY2KM(50), DY2K0S(50),
3036      &     DY2LA(50), DY2LB(50), DY2PHI(50),
3037      &     dm2pip(50), dm2pim(50), DMT2PR(50),
3038      &     DMT2PB(50), DMT2KP(50), dm2km(50),
3039      &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
3040      &     dy2msn(50), DY2PIP(50), DY2PIM(50), 
3041      &     DY2PI0(50), DY2PR(50), DY2PB(50)
3042      &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
3043 cbz3/17/99 end
3044 cc      SAVE /ARANA2/
3045       SAVE   
3046
3047       DO 1002 J = 1, NUM
3048          DO 1001 I = 1, MULTI1(J)
3049             ITYP = ITYP1(I, J)
3050             PX = PX1(I, J)
3051             PY = PY1(I, J)
3052             PZ = PZ1(I, J)
3053             EE = EE1(I, J)
3054             XM = XM1(I, J)
3055             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
3056 c     2/24/03 leptons and photons:
3057             if(xm.lt.0.01) goto 200
3058             ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
3059             eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
3060
3061             IF (ABS(PZ) .GE. EE) THEN
3062                PRINT *, 'IN ARTAN2'
3063                PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
3064 cbzdbg2/16/99
3065                PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3066 cbzdbg2/16/99
3067 cbzdbg2/15/99
3068                PRINT *, ' PZ = ', PZ, ' EE = ', EE
3069 cbzdbg2/16/99
3070                PRINT *, ' XM = ', XM
3071 cbzdbg2/16/99end
3072                GOTO 200
3073 c               STOP
3074 cbzdbg2/15/99end
3075             END IF
3076             DXMT = XMT - XM
3077             Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
3078 c.....rapidity cut for the rapidity distribution
3079             IF (ABS(Y) .GE. 10.0) GOTO 100
3080 c            IY = 1 + int(ABS(Y) / BY)
3081 c            Ieta = 1 + int(ABS(eta) / BY)
3082             IF (ABS(eta) .GE. 10.0) GOTO 100
3083             IY = 1 + int((Y+10.) / BY)
3084             Ieta = 1 + int((eta+10.) / BY)
3085
3086             IF (ITYP .LT. -1000) THEN
3087                dy2ntb(IY) = dy2ntb(IY) - 1.0
3088             END IF
3089             IF (ITYP .GT. 1000) THEN
3090                dy2ntb(IY) = dy2ntb(IY) + 1.0
3091             END IF
3092             IF (ITYP .EQ. -2212) THEN
3093                dy2ntp(IY) = dy2ntp(IY) - 1.0
3094             END IF
3095             IF (ITYP .EQ. 2212) THEN
3096                dy2ntp(IY) = dy2ntp(IY) + 1.0
3097             END IF
3098             IF (ITYP .EQ. -2112) THEN
3099                DY2HM(IY) = DY2HM(IY) + 1.0
3100             END IF
3101
3102             IF (LUCHGE(ITYP).ne.0) THEN
3103                DY2CH(IY) = DY2CH(IY) + 1.0
3104                DE2CH(Ieta) = DE2CH(Ieta) + 1.0
3105                IF (LUCHGE(ITYP).lt.0) THEN
3106                   DY2NEG(IY) = DY2NEG(IY) + 1.0
3107                   DE2NEG(Ieta) = DE2NEG(Ieta) + 1.0
3108                endif
3109             END IF
3110
3111 cbz3/17/99
3112             IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. 
3113      &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
3114                dy2msn(IY) = dy2msn(IY) + 1.0
3115             END IF
3116             IF (ITYP .EQ. 211) THEN
3117                DY2PIP(IY) = DY2PIP(IY) + 1.0
3118             END IF
3119             IF (ITYP .EQ. -211) THEN
3120                DY2PIM(IY) = DY2PIM(IY) + 1.0
3121             END IF
3122             IF (ITYP .EQ. 111) THEN
3123                DY2PI0(IY) = DY2PI0(IY) + 1.0
3124             END IF
3125             IF (ITYP .EQ. 2212) THEN
3126                DY2PR(IY) = DY2PR(IY) + 1.0
3127             END IF
3128             IF (ITYP .EQ. -2212) THEN
3129                DY2PB(IY) = DY2PB(IY) + 1.0
3130             END IF
3131 cbz3/17/99 end
3132             IF (ITYP .EQ. 321) THEN
3133                DY2KP(IY) = DY2KP(IY) + 1.0
3134             END IF
3135             IF (ITYP .EQ. -321) THEN
3136                DY2KM(IY) = DY2KM(IY) + 1.0
3137             END IF
3138 clin-4/24/03:
3139 c            IF (ITYP .EQ. 310) THEN
3140             IF (ITYP .EQ. 130) THEN
3141                DY2K0S(IY) = DY2K0S(IY) + 1.0
3142             END IF
3143             IF (ITYP .EQ. 3122) THEN
3144                DY2LA(IY) = DY2LA(IY) + 1.0
3145             END IF
3146             IF (ITYP .EQ. -3122) THEN
3147                DY2LB(IY) = DY2LB(IY) + 1.0
3148             END IF
3149             IF (ITYP .EQ. 333) THEN
3150                DY2PHI(IY) = DY2PHI(IY) + 1.0
3151             END IF
3152
3153 c.....insert rapidity cut for mt spectrum here
3154  100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
3155             IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
3156             IMT = 1 + int(DXMT / BMT)
3157             IF (ITYP .EQ. 211) THEN
3158                dm2pip(IMT) = dm2pip(IMT) + 1.0 / XMT
3159             END IF
3160             IF (ITYP .EQ. -211) THEN
3161                dm2pim(IMT) = dm2pim(IMT) + 
3162      &            1.0 / XMT
3163             END IF
3164             IF (ITYP .EQ. 2212) THEN
3165                DMT2PR(IMT) = DMT2PR(IMT) + 1.0 / XMT
3166             END IF
3167             IF (ITYP .EQ. -2212) THEN
3168                DMT2PB(IMT) = DMT2PB(IMT) + 1.0 / XMT
3169             END IF
3170             IF (ITYP .EQ. 321) THEN
3171                DMT2KP(IMT) = DMT2KP(IMT) + 1.0 / XMT
3172             END IF
3173             IF (ITYP .EQ. -321) THEN
3174                dm2km(IMT) = dm2km(IMT) + 1.0 / XMT
3175             END IF
3176 clin-4/24/03:
3177 c            IF (ITYP .EQ. 310) THEN
3178             IF (ITYP .EQ. 130) THEN
3179                dm2k0s(IMT) = dm2k0s(IMT) + 1.0 / XMT
3180             END IF
3181             IF (ITYP .EQ. 3122) THEN
3182                DMT2LA(IMT) = DMT2LA(IMT) + 1.0 / XMT
3183             END IF
3184             IF (ITYP .EQ. -3122) THEN
3185                DMT2LB(IMT) = DMT2LB(IMT) + 1.0 / XMT
3186             END IF
3187
3188  200        CONTINUE
3189  1001    CONTINUE
3190  1002 CONTINUE
3191
3192       RETURN
3193       END
3194
3195 c-----------------------------------------------------------------------
3196
3197 c.....output analysis results at the end of the simulation
3198
3199       SUBROUTINE ARTOUT(NEVNT)
3200
3201       PARAMETER (MAXSTR=150001, MAXR=1)
3202 c.....y cut for mt spectrum
3203 cbz3/17/99
3204 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3205       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3206 cbz3/17/99 end
3207 c.....bin width for mt spectrum and y spectrum
3208 c      PARAMETER (BMT = 0.05, BY = 0.2)
3209       PARAMETER (BMT = 0.05, BY = 0.4)
3210       COMMON /RUN/ NUM
3211 cc      SAVE /RUN/
3212       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
3213      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
3214      &     FT1(MAXSTR, MAXR),
3215      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
3216      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
3217 cbz3/17/99
3218 c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
3219 cc      SAVE /ARPRC1/
3220       COMMON /ARANA1/
3221      &     dy1ntb(50), dy1ntp(50), DY1HM(50), 
3222      &     DY1KP(50), DY1KM(50), DY1K0S(50),
3223      &     DY1LA(50), DY1LB(50), DY1PHI(50),
3224      &     dm1pip(50), dm1pim(50), DMT1PR(50),
3225      &     DMT1PB(50), DMT1KP(50), dm1km(50),
3226      &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
3227      &     dy1msn(50), DY1PIP(50), DY1PIM(50), 
3228      &     DY1PI0(50), DY1PR(50), DY1PB(50)
3229      &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
3230 cbz3/17/99 end
3231 cc      SAVE /ARANA1/
3232 cbz3/17/99
3233 c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
3234       COMMON /ARANA2/
3235      &     dy2ntb(50), dy2ntp(50), DY2HM(50), 
3236      &     DY2KP(50), DY2KM(50), DY2K0S(50),
3237      &     DY2LA(50), DY2LB(50), DY2PHI(50),
3238      &     dm2pip(50), dm2pim(50), DMT2PR(50),
3239      &     DMT2PB(50), DMT2KP(50), dm2km(50),
3240      &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
3241      &     dy2msn(50), DY2PIP(50), DY2PIM(50), 
3242      &     DY2PI0(50), DY2PR(50), DY2PB(50)
3243      &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
3244 cc      SAVE /ARANA2/
3245       SAVE   
3246 cbz3/17/99 end
3247 cms   OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN')
3248 cms   OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN')
3249 cms   OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN')
3250 cms   OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN')
3251 cms   OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN')
3252 cms   OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN')
3253 cms   OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN')
3254 cms   OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN')
3255 cms   OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN')
3256 clin-4/24/03
3257 c      OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN')
3258 cms   OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN')
3259 cms   OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN')
3260 cms   OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN')
3261 cms   OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN')
3262 cbz3/17/99
3263 cms   OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN')
3264 cms   OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN')
3265 cms   OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN')
3266 cms   OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN')
3267 cms   OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN')
3268 cms   OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN')
3269 cbz3/17/99 end
3270
3271 cms   OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN')
3272 cms   OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN')
3273 cms   OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN')
3274 cms   OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN')
3275 cms   OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN')
3276 cms   OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN')
3277 cms   OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN')
3278 cms   OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN')
3279 cms   OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN')
3280 clin-9/26/03 no symmetrization in y (or eta) for ana/*.dat:
3281 c      SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0
3282       SCALE1 = 1. / REAL(NEVNT * NUM) / BY
3283       SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1)
3284 c
3285       DO 1001 I = 1, 50
3286          ymid=-10.+BY * (I - 0.5)
3287 cms      WRITE (30, 333) ymid, SCALE1 * dy1ntb(I)
3288 cms      WRITE (31, 333) ymid, SCALE1 * dy1ntp(I)
3289 cms      WRITE (32, 333) ymid, SCALE1 * DY1HM(I)
3290 cms      WRITE (37, 333) ymid, SCALE1 * DY1KP(I)
3291 cms      WRITE (38, 333) ymid, SCALE1 * DY1KM(I)
3292 cms      WRITE (39, 333) ymid, SCALE1 * DY1K0S(I)
3293 cms      WRITE (40, 333) ymid, SCALE1 * DY1LA(I)
3294 cms      WRITE (41, 333) ymid, SCALE1 * DY1LB(I)
3295 cms      WRITE (42, 333) ymid, SCALE1 * DY1PHI(I)
3296 cms      WRITE (33, 333) ymid, SCALE1 * DY1NEG(I)
3297 cms      WRITE (34, 333) ymid, SCALE1 * DY1CH(I)
3298 cms      WRITE (35, 333) ymid, SCALE1 * DE1NEG(I)
3299 cms      WRITE (36, 333) ymid, SCALE1 * DE1CH(I)
3300 cms      WRITE (43, 333) ymid, SCALE1 * dy1msn(I)
3301 cms      WRITE (44, 333) ymid, SCALE1 * DY1PIP(I)
3302 cms      WRITE (45, 333) ymid, SCALE1 * DY1PIM(I)
3303 cms      WRITE (46, 333) ymid, SCALE1 * DY1PI0(I)
3304 cms      WRITE (47, 333) ymid, SCALE1 * DY1PR(I)
3305 cms      WRITE (48, 333) ymid, SCALE1 * DY1PB(I)
3306
3307          IF (dm1pip(I) .NE. 0.0) THEN
3308 cms         WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I)
3309          END IF
3310          IF (dm1pim(I) .NE. 0.0) THEN
3311 cms         WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * 
3312 cms  &         dm1pim(I)
3313          END IF
3314          IF (DMT1PR(I) .NE. 0.0) THEN
3315 cms         WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I)
3316          END IF
3317          IF (DMT1PB(I) .NE. 0.0) THEN
3318 cms         WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I)
3319          END IF
3320          IF (DMT1KP(I) .NE. 0.0) THEN
3321 cms         WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I)
3322          END IF
3323          IF (dm1km(I) .NE. 0.0) THEN
3324 cms         WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3325 cms  &         dm1km(I)
3326          END IF
3327          IF (dm1k0s(I) .NE. 0.0) THEN
3328 cms         WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I)
3329          END IF
3330          IF (DMT1LA(I) .NE. 0.0) THEN
3331 cms         WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I)
3332          END IF
3333          IF (DMT1LB(I) .NE. 0.0) THEN
3334 cms         WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I)
3335          END IF
3336  1001 CONTINUE
3337 c
3338       DO 1002 I = 30, 48
3339 cms      WRITE (I, *) 'after hadron evolution'
3340  1002    CONTINUE
3341       DO 1003 I = 50, 58
3342 cms      WRITE (I, *) 'after hadron evolution'
3343  1003 CONTINUE
3344
3345       DO 1004 I = 1, 50
3346          ymid=-10.+BY * (I - 0.5)
3347 cms      WRITE (30, 333) ymid, SCALE1 * dy2ntb(I)
3348 cms      WRITE (31, 333) ymid, SCALE1 * dy2ntp(I)
3349 cms      WRITE (32, 333) ymid, SCALE1 * DY2HM(I)
3350 cms      WRITE (37, 333) ymid, SCALE1 * DY2KP(I)
3351 cms      WRITE (38, 333) ymid, SCALE1 * DY2KM(I)
3352 cms      WRITE (39, 333) ymid, SCALE1 * DY2K0S(I)
3353 cms      WRITE (40, 333) ymid, SCALE1 * DY2LA(I)
3354 cms      WRITE (41, 333) ymid, SCALE1 * DY2LB(I)
3355 cms      WRITE (42, 333) ymid, SCALE1 * DY2PHI(I)
3356 cms      WRITE (33, 333) ymid, SCALE1 * DY2NEG(I)
3357 cms      WRITE (34, 333) ymid, SCALE1 * DY2CH(I)
3358 cms      WRITE (35, 333) ymid, SCALE1 * DE2NEG(I)
3359 cms      WRITE (36, 333) ymid, SCALE1 * DE2CH(I)
3360 cms      WRITE (43, 333) ymid, SCALE1 * dy2msn(I)
3361 cms      WRITE (44, 333) ymid, SCALE1 * DY2PIP(I)
3362 cms      WRITE (45, 333) ymid, SCALE1 * DY2PIM(I)
3363 cms      WRITE (46, 333) ymid, SCALE1 * DY2PI0(I)
3364 cms      WRITE (47, 333) ymid, SCALE1 * DY2PR(I)
3365 cms      WRITE (48, 333) ymid, SCALE1 * DY2PB(I)
3366 c
3367          IF (dm2pip(I) .NE. 0.0) THEN
3368 cms         WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I)
3369          END IF
3370          IF (dm2pim(I) .NE. 0.0) THEN
3371 cms         WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * 
3372 cms  &         dm2pim(I)
3373          END IF
3374          IF (DMT2PR(I) .NE. 0.0) THEN
3375 cms         WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I)
3376          END IF
3377          IF (DMT2PB(I) .NE. 0.0) THEN
3378 cms         WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I)
3379          END IF
3380          IF (DMT2KP(I) .NE. 0.0) THEN
3381 cms         WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I)
3382          END IF
3383          IF (dm2km(I) .NE. 0.0) THEN
3384 cms         WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * 
3385 cms  &         dm2km(I)
3386          END IF
3387          IF (dm2k0s(I) .NE. 0.0) THEN
3388 cms         WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I)
3389          END IF
3390          IF (DMT2LA(I) .NE. 0.0) THEN
3391 cms         WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I)
3392          END IF
3393          IF (DMT2LB(I) .NE. 0.0) THEN
3394 cms         WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I)
3395          END IF
3396  1004 CONTINUE
3397 cms 333  format(2(f12.5,1x))
3398
3399       RETURN
3400       END
3401
3402 c-----------------------------------------------------------------------
3403
3404 c.....analysis subroutine in HIJING before parton cascade evolution
3405       SUBROUTINE HJANA1
3406
3407       PARAMETER (YMAX = 1.0, YMIN = -1.0)
3408       PARAMETER (DMT = 0.05, DY = 0.2)
3409       PARAMETER (DR = 0.2)
3410       PARAMETER (MAXPTN=400001,MAXSTR=150001)
3411       DIMENSION dyp1(50), DMYP1(200), DEYP1(50)
3412       DIMENSION dyg1(50), DMYG1(200), DEYG1(50)
3413       DIMENSION SNYP1(50), SMYP1(200), SEYP1(50)
3414       DIMENSION SNYG1(50), SMYG1(200), SEYG1(50)
3415       DIMENSION dnrpj1(50), dnrtg1(50), dnrin1(50),
3416      &   dnrtt1(50)
3417       DIMENSION dyg1c(50), dmyg1c(50), deyg1c(50)
3418       DIMENSION snrpj1(50), snrtg1(50), snrin1(50),
3419      &   snrtt1(50)
3420       DIMENSION snyg1c(50), smyg1c(50), seyg1c(50)
3421       DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
3422
3423       COMMON /PARA1/ MUL
3424 cc      SAVE /PARA1/
3425       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3426 cc      SAVE /HPARNT/
3427       COMMON/hjcrdn/YP(3,300),YT(3,300)
3428 cc      SAVE /hjcrdn/
3429       COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3430      &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3431      &   PJPM(300,500),NTJ(300),KFTJ(300,500),
3432      &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3433      &   PJTE(300,500),PJTM(300,500)
3434 cc      SAVE /HJJET1/
3435       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3436      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3437      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3438 cc      SAVE /HJJET2/
3439       COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
3440      &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
3441      &     XMASS0(MAXPTN), ITYP0(MAXPTN)
3442 cc      SAVE /prec1/
3443       COMMON /AREVT/ IAEVT, IARUN, MISS
3444 cc      SAVE /AREVT/
3445       COMMON /AROUT/ IOUT
3446 cc      SAVE /AROUT/
3447       SAVE   
3448       DATA IW/0/
3449       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3450          DO 1001 I = 1, 200
3451             DMYP1(I) = SMYP1(I)
3452             DMYG1(I) = SMYG1(I)
3453  1001    CONTINUE
3454
3455          DO 1002 I = 1, 50
3456             dyp1(I) = SNYP1(I)
3457             DEYP1(I) = SEYP1(I)
3458             dyg1(I) = SNYG1(I)
3459             DEYG1(I) = SEYG1(I)
3460             dnrpj1(I) = snrpj1(I)
3461             dnrtg1(I) = snrtg1(I)
3462             dnrin1(I) = snrin1(I)
3463             dnrtt1(I) = snrtt1(I)
3464             dyg1c(I) = snyg1c(I)
3465             dmyg1c(I) = smyg1c(I)
3466             deyg1c(I) = seyg1c(I)
3467  1002    CONTINUE
3468          nsubp = nsubpS
3469          nsubg = nsubgS
3470          NISG = NISGS
3471       ELSE
3472          DO 1003 I = 1, 200
3473             SMYP1(I) = DMYP1(I)
3474             SMYG1(I) = DMYG1(I)
3475  1003    CONTINUE
3476
3477          DO 1004 I = 1, 50
3478             SNYP1(I) = dyp1(I)
3479             SEYP1(I) = DEYP1(I)
3480             SNYG1(I) = dyg1(I)
3481             SEYG1(I) = DEYG1(I)
3482             snrpj1(I) = dnrpj1(I)
3483             snrtg1(I) = dnrtg1(I)
3484             snrin1(I) = dnrin1(I)
3485             snrtt1(I) = dnrtt1(I)
3486             snyg1c(I) = dyg1c(I)
3487             smyg1c(I) = dmyg1c(I)
3488             seyg1c(I) = deyg1c(I)
3489  1004    CONTINUE
3490          nsubpS = nsubp
3491          nsubgS = nsubg
3492          NISGS = NISG
3493          isevt = IAEVT
3494          isrun = IARUN
3495          IW = IW + 1
3496       END IF
3497 c.....analysis
3498       DO 1006 I = 1, IHNT2(1)
3499          DO 1005 J = 1, NPJ(I)