forgotten change
[u/mrichter/AliRoot.git] / TAmpt / AMPT / art1f.f
1 c....................art1f.f
2 **************************************
3 *
4 *                           PROGRAM ART1.0 
5 *
6 *        A relativistic transport (ART) model for heavy-ion collisions
7 *
8 *   sp/01/04/2002
9 *   calculates K+K- from phi decay, dimuons from phi decay
10 *   has finite baryon density & possibilites of varying Kaon 
11 *   in-medium mass in phiproduction-annhilation channel only.
12 *
13 *
14 * RELEASING DATE: JAN., 1997 
15 ***************************************
16
17 * Bao-An Li & Che Ming Ko
18 * Cyclotron Institute, Texas A&M University.
19 * Phone: (409) 845-1411
20 * e-mail: Bali@comp.tamu.edu & Ko@comp.tamu.edu 
21 * http://wwwcyc.tamu.edu/~bali
22 ***************************************
23 * Speical notice on the limitation of the code:
24
25 * (1) ART is a hadronic transport model
26
27 * (2) E_beam/A <= 15 GeV
28
29 * (3) The mass of the colliding system is limited by the dimensions of arrays
30 *    which can be extended purposely. Presently the dimensions are large enough
31 *     for running Au+Au at 15 GeV/A.
32 *
33 * (4) The production and absorption of antiparticles (e.g., ki-, anti-nucleons,
34 *     etc) are not fully included in this version of the model. They, however, 
35 *     have essentially no effect on the reaction dynamics and observables 
36 *     related to nucleons, pions and kaons (K+) at and below AGS energies.
37
38 * (5) Bose enhancement for mesons and Pauli blocking for fermions are 
39 *     turned off.
40
41 *********************************
42 *
43 * USEFUL REFERENCES ON PHYSICS AND NUMERICS OF NUCLEAR TRANSPORT MODELS:
44 *     G.F. BERTSCH AND DAS GUPTA, PHYS. REP. 160 (1988) 189.
45 *     B.A. LI AND W. BAUER, PHYS. REV. C44 (1991) 450.
46 *     B.A. LI, W. BAUER AND G.F. BERTSCH, PHYS. REV. C44 (1991) 2095.
47 *     P. DANIELEWICZ AND G.F. BERTSCH, NUCL. PHYS. A533 (1991) 712.
48
49 * MAIN REFERENCES ON THIS VERSION OF ART MODEL:
50 *     B.A. LI AND C.M. KO, PHYS. REV. C52 (1995) 2037; 
51 *                          NUCL. PHYS. A601 (1996) 457. 
52 *
53 **********************************
54 **********************************
55 *  VARIABLES IN INPUT-SECTION:                                               * 
56 *                                                                      *
57 *  1) TARGET-RELATED QUANTITIES                                        *
58 *       MASSTA, ZTA -  TARGET MASS IN AMU, TARGET CHARGE  (INTEGER)    *
59 *                                                                      *
60 *  2) PROJECTILE-RELATED QUANTITIES                                    *
61 *       MASSPR, ZPR -  PROJECTILE MASS IN AMU, PROJ. CHARGE(INTEGER)   *
62 *       ELAB     -  BEAM ENERGY IN [MEV/NUCLEON]               (REAL)  *
63 *       ZEROPT   -  DISPLACEMENT OF THE SYSTEM IN Z-DIREC. [FM](REAL)  *
64 *       B        -  IMPACT PARAMETER [FM]                      (REAL)  *
65 *                                                                      *
66 *  3) PROGRAM-CONTROL PARAMETERS                                       *
67 *       ISEED    -  SEED FOR RANDOM NUMBER GENERATOR        (INTEGER)  *
68 *       DT       -  TIME-STEP-SIZE [FM/C]                      (REAL)  *
69 *       NTMAX    -  TOTAL NUMBER OF TIMESTEPS               (INTEGER)  *
70 *       ICOLL    -  (= 1 -> MEAN FIELD ONLY,                           *
71 *                -   =-1 -> CACADE ONLY, ELSE FULL ART)     (INTEGER)  *
72 *       NUM      -  NUMBER OF TESTPARTICLES PER NUCLEON     (INTEGER)  *
73 *       INSYS    -  (=0 -> LAB-SYSTEM, ELSE C.M. SYSTEM)    (INTEGER)  *
74 *       IPOT     -  1 -> SIGMA=2; 2 -> SIGMA=4/3; 3 -> SIGMA=7/6       *
75 *                   IN MEAN FIELD POTENTIAL                 (INTEGER)  *
76 *       MODE     -  (=1 -> interpolation for pauli-blocking,           *
77 *                    =2 -> local lookup, other -> unblocked)(integer)  *
78 *       DX,DY,DZ -  widths of cell for paulat in coor. sp. [fm](real)  *
79 *       DPX,DPY,DPZ-widths of cell for paulat in mom. sp.[GeV/c](real) *
80 *       IAVOID   -  (=1 -> AVOID FIRST COLL. WITHIN SAME NUCL.         *
81 *                    =0 -> ALLOW THEM)                      (INTEGER)  *
82 *       IMOMEN   -  FLAG FOR CHOICE OF INITIAL MOMENTUM DISTRIBUTION   *
83 *                   (=1 -> WOODS-SAXON DENSITY AND LOCAL THOMAS-FERMI  *
84 *                    =2 -> NUCLEAR MATTER DEN. AND LOCAL THOMAS-FERMI  *
85 *                    =3 -> COHERENT BOOST IN Z-DIRECTION)   (INTEGER)  *
86 *  4) CONTROL-PRINTOUT OPTIONS                                         *
87 *       NFREQ    -  NUMBER OF TIMSTEPS AFTER WHICH PRINTOUT            *
88 *                   IS REQUIRED OR ON-LINE ANALYSIS IS PERFORMED       *
89 *       ICFLOW      =1 PERFORM ON-LINE FLOW ANALYSIS EVERY NFREQ STEPS *
90 *       ICRHO       =1 PRINT OUT THE BARYON,PION AND ENERGY MATRIX IN  *
91 *                      THE REACTION PLANE EVERY NFREQ TIME-STEPS       *
92 *  5)
93 *       CYCBOX   -  ne.0 => cyclic boundary conditions;boxsize CYCBOX  *
94 *
95 **********************************
96 *               Lables of particles used in this code                     *
97 **********************************
98 *         
99 *         LB(I) IS USED TO LABEL PARTICLE'S CHARGE STATE
100 *    
101 *         LB(I)   =
102 clin-11/07/00:
103 *                -30 K*-
104 clin-8/29/00
105 *                -13 anti-N*(+1)(1535),s_11
106 *                -12 anti-N*0(1535),s_11
107 *                 -11 anti-N*(+1)(1440),p_11
108 *                 -10 anti-N*0(1440), p_11
109 *                  -9 anti-DELTA+2
110 *                  -8 anti-DELTA+1
111 *                  -7 anti-DELTA0
112 *                  -6 anti-DELTA-1
113 clin-8/29/00-end
114
115 cbali2/7/99 
116 *                  -2 antineutron 
117 *                             -1       antiproton
118 cbali2/7/99 end 
119 *                   0 eta
120 *                        1 PROTON
121 *                   2 NUETRON
122 *                   3 PION-
123 *                   4 PION0
124 *                   5 PION+
125 *                   6 DELTA-1
126 *                   7 DELTA0
127 *                   8 DELTA+1
128 *                   9 DELTA+2
129 *                   10 N*0(1440), p_11
130 *                   11 N*(+1)(1440),p_11
131 *                  12 N*0(1535),s_11
132 *                  13 N*(+1)(1535),s_11
133 *                  14 LAMBDA
134 *                   15 sigma-, since we used isospin averaged xsection for
135 *                   16 sigma0  sigma associated K+ production, sigma0 and 
136 *                   17 sigma+  sigma+ are counted as sigma-
137 *                   21 kaon-
138 *                   23 KAON+
139 *                   24 kaon0
140 *                   25 rho-
141 *                         26 rho0
142 *                   27 rho+
143 *                   28 omega meson
144 *                   29 phi
145 clin-11/07/00:
146 *                  30 K*+
147 * sp01/03/01
148 *                 -14 LAMBDA(bar)
149 *                  -15 sigma-(bar)
150 *                  -16 sigma0(bar)
151 *                  -17 sigma+(bar)
152 *                   31 eta-prime
153 *                   40 cascade-
154 *                  -40 cascade-(bar)
155 *                   41 cascade0
156 *                  -41 cascade0(bar)
157 *                   45 Omega baryon
158 *                  -45 Omega baryon(bar)
159 * sp01/03/01 end
160 clin-5/2008:
161 *                   42 Deuteron (same in ampt.dat)
162 *                  -42 anti-Deuteron (same in ampt.dat)
163 c
164 *                   ++  ------- SEE BAO-AN LI'S NOTE BOOK
165 **********************************
166 cbz11/16/98
167 c      PROGRAM ART
168        SUBROUTINE ARTMN
169 cbz11/16/98end
170 **********************************
171 * PARAMETERS:                                                           *
172 *  MAXPAR     - MAXIMUM NUMBER OF PARTICLES      PROGRAM CAN HANDLE     *
173 *  MAXP       - MAXIMUM NUMBER OF CREATED MESONS PROGRAM CAN HANDLE     *
174 *  MAXR       - MAXIMUM NUMBER OF EVENTS AT EACH IMPACT PARAMETER       *
175 *  MAXX       - NUMBER OF MESHPOINTS IN X AND Y DIRECTION = 2 MAXX + 1  *
176 *  MAXZ       - NUMBER OF MESHPOINTS IN Z DIRECTION       = 2 MAXZ + 1  *
177 *  AMU        - 1 ATOMIC MASS UNIT "GEV/C**2"                           *
178 *  MX,MY,MZ   - MESH SIZES IN COORDINATE SPACE [FM] FOR PAULI LATTICE   *
179 *  MPX,MPY,MPZ- MESH SIZES IN MOMENTUM SPACE [GEV/C] FOR PAULI LATTICE  *
180 *---------------------------------------------------------------------- *
181 clin      PARAMETER     (maxpar=200000,MAXR=50,AMU= 0.9383,
182       PARAMETER     (MAXSTR=150001,MAXR=1,AMU= 0.9383,
183      1               AKA=0.498,etaM=0.5475)
184       PARAMETER     (MAXX   =   20,  MAXZ  =    24)
185       PARAMETER     (ISUM   =   1001,  IGAM  =    1100)
186       parameter     (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
187 clin      PARAMETER (MAXP = 14000)
188 *----------------------------------------------------------------------*
189       INTEGER   OUTPAR, zta,zpr
190       COMMON  /AA/      R(3,MAXSTR)
191 cc      SAVE /AA/
192       COMMON  /BB/      P(3,MAXSTR)
193 cc      SAVE /BB/
194       COMMON  /CC/      E(MAXSTR)
195 cc      SAVE /CC/
196       COMMON  /DD/      RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
197      &                     RHOP(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
198      &                     RHON(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
199 cc      SAVE /DD/
200       COMMON  /EE/      ID(MAXSTR),LB(MAXSTR)
201 cc      SAVE /EE/
202       COMMON  /HH/  PROPER(MAXSTR)
203 cc      SAVE /HH/
204       common  /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
205 cc      SAVE /ff/
206       common  /gg/      dx,dy,dz,dpx,dpy,dpz
207 cc      SAVE /gg/
208       COMMON  /INPUT/ NSTAR,NDIRCT,DIR
209 cc      SAVE /INPUT/
210       COMMON  /PP/      PRHO(-20:20,-24:24)
211       COMMON  /QQ/      PHRHO(-MAXZ:MAXZ,-24:24)
212       COMMON  /RR/      MASSR(0:MAXR)
213 cc      SAVE /RR/
214       common  /ss/      inout(20)
215 cc      SAVE /ss/
216       common  /zz/      zta,zpr
217 cc      SAVE /zz/
218       COMMON  /RUN/     NUM
219 cc      SAVE /RUN/
220 clin-4/2008:
221 c      COMMON  /KKK/     TKAON(7),EKAON(7,0:200)
222       COMMON  /KKK/     TKAON(7),EKAON(7,0:2000)
223 cc      SAVE /KKK/
224       COMMON  /KAON/    AK(3,50,36),SPECK(50,36,7),MF
225 cc      SAVE /KAON/
226       COMMON/TABLE/ xarray(0:1000),earray(0:1000)
227 cc      SAVE /TABLE/
228       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
229 cc      SAVE /input1/
230       COMMON  /DDpi/    piRHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
231 cc      SAVE /DDpi/
232       common  /tt/  PEL(-maxx:maxx,-maxx:maxx,-maxz:maxz)
233      &,rxy(-maxx:maxx,-maxx:maxx,-maxz:maxz)
234 cc      SAVE /tt/
235 clin-4/2008:
236 c      DIMENSION TEMP(3,MAXSTR),SKAON(7),SEKAON(7,0:200)
237       DIMENSION TEMP(3,MAXSTR),SKAON(7),SEKAON(7,0:2000)
238 cbz12/2/98
239       COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, 
240      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
241 cc      SAVE /INPUT2/
242       COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
243 cc      SAVE /INPUT3/
244 cbz12/2/98end
245 cbz11/16/98
246       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
247 cc      SAVE /ARPRNT/
248
249 c.....note in the below, since a common block in ART is called EE,
250 c.....the variable EE in /ARPRC/is changed to PEAR.
251 clin-9/29/03 changed name in order to distinguish from /prec2/
252 c        COMMON /ARPRC/ ITYPAR(MAXSTR),
253 c     &       GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
254 c     &       PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
255 c     &       XMAR(MAXSTR)
256 cc      SAVE /ARPRC/
257 clin-9/29/03-end
258       COMMON /ARERCP/PRO1(MAXSTR, MAXR)
259 cc      SAVE /ARERCP/
260       COMMON /ARERC1/MULTI1(MAXR)
261 cc      SAVE /ARERC1/
262       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
263      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
264      &     FT1(MAXSTR, MAXR),
265      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
266      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
267 cc      SAVE /ARPRC1/
268 c
269       DIMENSION NPI(MAXR)
270       DIMENSION RT(3, MAXSTR, MAXR), PT(3, MAXSTR, MAXR)
271      &     , ET(MAXSTR, MAXR), LT(MAXSTR, MAXR), PROT(MAXSTR, MAXR)
272
273       EXTERNAL IARFLV, INVFLV
274 cbz11/16/98end
275       common /lastt/itimeh,bimp 
276 cc      SAVE /lastt/
277       common/snn/efrm,npart1,npart2
278 cc      SAVE /snn/
279       COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
280 cc      SAVE /hbt/
281       common/resdcy/NSAV,iksdcy
282 cc      SAVE /resdcy/
283       COMMON/RNDF77/NSEED
284 cc      SAVE /RNDF77/
285       COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
286       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
287      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
288      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
289       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
290 clin-4/2008 zet() expanded to avoid out-of-bound errors:
291       real zet(-45:45)
292       SAVE   
293       data zet /
294      4     1.,0.,0.,0.,0.,
295      3     1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
296      2     -1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
297      1     0.,0.,0.,-1.,0.,1.,0.,-1.,0.,-1.,
298      s     0.,-2.,-1.,0.,1.,0.,0.,0.,0.,-1.,
299      e     0.,
300      s     1.,0.,-1.,0.,1.,-1.,0.,1.,2.,0.,
301      1     1.,0.,1.,0.,-1.,0.,1.,0.,0.,0.,
302      2     -1.,0.,1.,0.,-1.,0.,1.,0.,0.,1.,
303      3     0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,
304      4     0.,0.,0.,0.,-1./
305
306       nlast=0
307       do 1002 i=1,MAXSTR
308          ftsv(i)=0.
309          do 1101 irun=1,maxr
310             ftsvt(i,irun)=0.
311  1101    continue
312          lblast(i)=999
313          do 1001 j=1,4
314 clin-4/2008 bugs pointed out by Vander Molen & Westfall:
315 c            xlast(i,j)=0.
316 c            plast(i,j)=0.
317             xlast(j,i)=0.
318             plast(j,i)=0.
319  1001    continue
320  1002 continue
321
322 *-------------------------------------------------------------------*
323 * Input information about the reaction system and contral parameters* 
324 *-------------------------------------------------------------------*
325 *              input section starts here                           *
326 *-------------------------------------------------------------------*
327
328 cbz12/2/98
329 c.....input section is moved to subroutine ARTSET
330 cbz12/2/98end
331
332 *-----------------------------------------------------------------------*
333 *                   input section ends here                            *
334 *-----------------------------------------------------------------------*
335 * read in the table for gengrating the transverse momentum
336 * IN THE NN-->DDP PROCESS
337        call tablem
338 * several control parameters, keep them fixed in this code. 
339        ikaon=1
340        nstar=1
341        ndirct=0
342        dir=0.02
343        asy=0.032
344        ESBIN=0.04
345        MF=36
346 *----------------------------------------------------------------------*
347 c      CALL FRONT(12,MASSTA,MASSPR,ELAB)
348 *----------------------------------------------------------------------*
349       RADTA  = 1.124 * FLOAT(MASSTA)**(1./3.)
350       RADPR  = 1.124 * FLOAT(MASSPR)**(1./3.)
351       ZDIST  = RADTA + RADPR
352 c      if ( cycbox.ne.0 ) zdist=0
353       BMAX   = RADTA + RADPR
354       MASS   = MASSTA + MASSPR
355       NTOTAL = NUM * MASS
356 *
357       IF (NTOTAL .GT. MAXSTR) THEN
358         WRITE(12,'(//10X,''**** FATAL ERROR: TOO MANY TEST PART. ****'//
359      & ' '')')
360         STOP
361       END IF
362 *
363 *-----------------------------------------------------------------------
364 *       RELATIVISTIC KINEMATICS
365 *
366 *       1) LABSYSTEM
367 *
368       ETA    = FLOAT(MASSTA) * AMU
369       PZTA   = 0.0
370       BETATA = 0.0
371       GAMMTA = 1.0
372 *
373       EPR    = FLOAT(MASSPR) * (AMU + 0.001 * ELAB)
374       PZPR   = SQRT( EPR**2 - (AMU * FLOAT(MASSPR))**2 )
375       BETAPR = PZPR / EPR
376       GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
377 *
378 * BETAC AND GAMMAC OF THE C.M. OBSERVED IN THE LAB. FRAME
379         BETAC=(PZPR+PZTA)/(EPR+ETA)
380         GAMMC=1.0 / SQRT(1.-BETAC**2)
381 *
382 c      WRITE(12,'(/10x,''****    KINEMATICAL PARAMETERS    ****''/)')
383 c      WRITE(12,'(10x,''1) LAB-FRAME:        TARGET PROJECTILE'')')
384 c      WRITE(12,'(10x,''   ETOTAL "GEV" '',2F11.4)') ETA, EPR
385 c      WRITE(12,'(10x,''   P "GEV/C"    '',2F11.4)') PZTA, PZPR
386 c      WRITE(12,'(10x,''   BETA         '',2F11.4)') BETATA, BETAPR
387 c      WRITE(12,'(10x,''   GAMMA        '',2F11.4)') GAMMTA, GAMMPR
388       IF (INSYS .NE. 0) THEN
389 *
390 *       2) C.M. SYSTEM
391 *
392         S      = (EPR+ETA)**2 - PZPR**2
393         xx1=4.*alog(float(massta))
394         xx2=4.*alog(float(masspr))
395         xx1=exp(xx1)
396         xx2=exp(xx2)
397         PSQARE = (S**2 + (xx1+ xx2) * AMU**4
398      &             - 2.0 * S * AMU**2 * FLOAT(MASSTA**2 + MASSPR**2)
399      &             - 2.0 * FLOAT(MASSTA**2 * MASSPR**2) * AMU**4)
400      &           / (4.0 * S)
401 *
402         ETA    = SQRT ( PSQARE + (FLOAT(MASSTA) * AMU)**2 )
403         PZTA   = - SQRT(PSQARE)
404         BETATA = PZTA / ETA
405         GAMMTA = 1.0 / SQRT( 1.0 - BETATA**2 )
406 *
407         EPR    = SQRT ( PSQARE + (FLOAT(MASSPR) * AMU)**2 )
408         PZPR   = SQRT(PSQARE)
409         BETAPR = PZPR/ EPR
410         GAMMPR = 1.0 / SQRT( 1.0 - BETAPR**2 )
411 *
412 c        WRITE(12,'(10x,''2) C.M.-FRAME:  '')')
413 c        WRITE(12,'(10x,''   ETOTAL "GEV" '',2F11.4)') ETA, EPR
414 c        WRITE(12,'(10x,''   P "GEV/C"    '',2F11.4)') PZTA, PZPR
415 c        WRITE(12,'(10x,''   BETA         '',2F11.4)') BETATA, BETAPR
416 c        WRITE(12,'(10x,''   GAMMA        '',2F11.4)') GAMMTA, GAMMPR
417 c        WRITE(12,'(10x,''S "GEV**2"      '',F11.4)')  S
418 c        WRITE(12,'(10x,''PSQARE "GEV/C"2 '',E14.3)')  PSQARE
419 c        WRITE(12,'(/10x,''*** CALCULATION DONE IN CM-FRAME ***''/)')
420       ELSE
421 c        WRITE(12,'(/10x,''*** CALCULATION DONE IN LAB-FRAME ***''/)')
422       END IF
423 * MOMENTUM PER PARTICLE
424       PZTA = PZTA / FLOAT(MASSTA)
425       PZPR = PZPR / FLOAT(MASSPR)
426 * total initial energy in the N-N cms frame
427       ECMS0=ETA+EPR
428 *-----------------------------------------------------------------------
429 *
430 * Start loop over many runs of different impact parameters
431 * IF MANYB=1, RUN AT A FIXED IMPACT PARAMETER B0, OTHERWISE GENERATE 
432 * MINIMUM BIAS EVENTS WITHIN THE IMPACT PARAMETER RANGE OF B_MIN AND B_MAX
433        DO 50000 IMANY=1,MANYB
434 *------------------------------------------------------------------------
435 * Initialize the impact parameter B
436        if (manyb. gt.1) then
437 111       BX=1.0-2.0*RANART(NSEED)
438        BY=1.0-2.0*RANART(NSEED)
439        B2=BX*BX+BY*BY
440        IF(B2.GT.1.0) GO TO 111       
441        B=SQRT(B2)*(BM-BI)+BI
442        ELSE
443        B=B0
444        ENDIF
445 c      WRITE(12,'(///10X,''RUN NUMBER:'',I6)') IMANY       
446 c      WRITE(12,'(//10X,''IMPACT PARAMETER B FOR THIS RUN:'',
447 c     &             F9.3,'' FM''/10X,49(''*'')/)') B
448 *
449 *-----------------------------------------------------------------------
450 *       INITIALIZATION
451 *1 INITIALIZATION IN ISOSPIN SPACE FOR BOTH THE PROJECTILE AND TARGET
452       call coulin(masspr,massta,NUM)
453 *2 INITIALIZATION IN PHASE SPACE FOR THE TARGET
454       CALL INIT(1       ,MASSTA   ,NUM     ,RADTA,
455      &          B/2.    ,ZEROPT+ZDIST/2.   ,PZTA,
456      &          GAMMTA  ,ISEED    ,MASS    ,IMOMEN)
457 *3.1 INITIALIZATION IN PHASE SPACE FOR THE PROJECTILE
458       CALL INIT(1+MASSTA,MASS     ,NUM     ,RADPR,
459      &          -B/2.   ,ZEROPT-ZDIST/2.   ,PZPR,
460      &          GAMMPR  ,ISEED    ,MASS    ,IMOMEN)
461 *3.2 OUTPAR IS THE NO. OF ESCAPED PARTICLES
462       OUTPAR = 0
463 *3.3 INITIALIZATION FOR THE NO. OF PARTICLES IN EACH SAMPLE
464 *    THIS IS NEEDED DUE TO THE FACT THAT PIONS CAN BE PRODUCED OR ABSORBED
465       MASSR(0)=0
466       DO 1003 IR =1,NUM
467       MASSR(IR)=MASS
468  1003 CONTINUE
469 *3.4 INITIALIZation FOR THE KAON SPECTRUM
470 *      CALL KSPEC0(BETAC,GAMMC)
471 * calculate the local baryon density matrix
472       CALL DENS(IPOT,MASS,NUM,OUTPAR)
473 *
474 *-----------------------------------------------------------------------
475 *       CONTROL PRINTOUT OF INITIAL CONFIGURATION
476 *
477 *      WRITE(12,'(''**********  INITIAL CONFIGURATION  **********''/)')
478 *
479 c print out the INITIAL density matrix in the reaction plane
480 c       do ix=-10,10
481 c       do iz=-10,10
482 c       write(1053,992)ix,iz,rho(ix,0,iz)/0.168
483 c       end do
484 c       end do
485 *-----------------------------------------------------------------------
486 *       CALCULATE MOMENTA FOR T = 0.5 * DT 
487 *       (TO OBTAIN 2ND DEGREE ACCURACY!)
488 *       "Reference: J. AICHELIN ET AL., PHYS. REV. C31, 1730 (1985)"
489 *
490       IF (ICOLL .NE. -1) THEN
491         DO 700 I = 1,NTOTAL
492           IX = NINT( R(1,I) )
493           IY = NINT( R(2,I) )
494           IZ = NINT( R(3,I) )
495 clin-4/2008 check bounds:
496           IF(IX.GE.MAXX.OR.IY.GE.MAXX.OR.IZ.GE.MAXZ
497      1         .OR.IX.LE.-MAXX.OR.IY.LE.-MAXX.OR.IZ.LE.-MAXZ) goto 700
498           CALL GRADU(IPOT,IX,IY,IZ,GRADX,GRADY,GRADZ)
499           P(1,I) = P(1,I) - (0.5 * DT) * GRADX
500           P(2,I) = P(2,I) - (0.5 * DT) * GRADY
501           P(3,I) = P(3,I) - (0.5 * DT) * GRADZ
502   700   CONTINUE
503       END IF
504 *-----------------------------------------------------------------------
505 *-----------------------------------------------------------------------
506 *4 INITIALIZATION OF TIME-LOOP VARIABLES
507 *4.1 COLLISION NUMBER COUNTERS
508 clin 51      RCNNE  = 0
509         RCNNE  = 0
510        RDD  = 0
511        RPP  = 0
512        rppk = 0
513        RPN  = 0
514        rpd  = 0
515        RKN  = 0
516        RNNK = 0
517        RDDK = 0
518        RNDK = 0
519       RCNND  = 0
520       RCNDN  = 0
521       RCOLL  = 0
522       RBLOC  = 0
523       RDIRT  = 0
524       RDECAY = 0
525       RRES   = 0
526 *4.11 KAON PRODUCTION PROBABILITY COUNTER FOR PERTURBATIVE CALCULATIONS ONLY
527       DO 1005 KKK=1,5
528          SKAON(KKK)  = 0
529          DO 1004 IS=1,2000
530             SEKAON(KKK,IS)=0
531  1004    CONTINUE
532  1005 CONTINUE
533 *4.12 anti-proton and anti-kaon counters
534        pr0=0.
535        pr1=0.
536        ska0=0.
537        ska1=0.
538 *       ============== LOOP OVER ALL TIME STEPS ================       *
539 *                             STARTS HERE                              *
540 *       ========================================================       *
541 cbz11/16/98
542       IF (IAPAR2(1) .NE. 1) THEN
543          DO 1016 I = 1, MAXSTR
544             DO 1015 J = 1, 3
545                R(J, I) = 0.
546                P(J, I) = 0.
547  1015       CONTINUE
548             E(I) = 0.
549             LB(I) = 0
550 cbz3/25/00
551             ID(I)=0
552 c     sp 12/19/00
553            PROPER(I) = 1.
554  1016   CONTINUE
555          MASS = 0
556 cbz12/22/98
557 c         MASSR(1) = 0
558 c         NP = 0
559 c         NPI = 1
560          NP = 0
561          DO 1017 J = 1, NUM
562             MASSR(J) = 0
563             NPI(J) = 1
564  1017    CONTINUE
565          DO 1019 I = 1, MAXR
566             DO 1018 J = 1, MAXSTR
567                RT(1, J, I) = 0.
568                RT(2, J, I) = 0.
569                RT(3, J, I) = 0.
570                PT(1, J, I) = 0.
571                PT(2, J, I) = 0.
572                PT(3, J, I) = 0.
573                ET(J, I) = 0.
574                LT(J, I) = 0
575 c     sp 12/19/00
576                PROT(J, I) = 1.
577  1018       CONTINUE
578  1019    CONTINUE
579 cbz12/22/98end
580       END IF
581 cbz11/16/98end
582         
583       DO 10000 NT = 1,NTMAX
584
585 *TEMPORARY PARTICLE COUNTERS
586 *4.2 PION COUNTERS : LP1,LP2 AND LP3 ARE THE NO. OF P+,P0 AND P-
587       LP1=0
588       LP2=0
589       LP3=0
590 *4.3 DELTA COUNTERS : LD1,LD2,LD3 AND LD4 ARE THE NO. OF D++,D+,D0 AND D-
591       LD1=0
592       LD2=0
593       LD3=0
594       LD4=0
595 *4.4 N*(1440) COUNTERS : LN1 AND LN2 ARE THE NO. OF N*+ AND N*0
596       LN1=0
597       LN2=0
598 *4.5 N*(1535) counters
599       LN5=0
600 *4.6 ETA COUNTERS
601       LE=0
602 *4.7 KAON COUNTERS
603       LKAON=0
604
605 clin-11/09/00:
606 * KAON* COUNTERS
607       LKAONS=0
608
609 *-----------------------------------------------------------------------
610         IF (ICOLL .NE. 1) THEN
611 * STUDYING BINARY COLLISIONS AMONG PARTICLES DURING THIS TIME INTERVAL *
612 clin-10/25/02 get rid of argument usage mismatch in relcol(.nt.):
613            numnt=nt
614           CALL RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk,
615      &    LPN,lpd,LRHO,LOMEGA,LKN,LNNK,LDDK,LNDK,LCNND,
616      &    LCNDN,LDIRT,LDECAY,LRES,LDOU,LDDRHO,LNNRHO,
617      &    LNNOM,numnt,ntmax,sp,akaon,sk)
618 c     &    LNNOM,NT,ntmax,sp,akaon,sk)
619 clin-10/25/02-end
620 *-----------------------------------------------------------------------
621
622 c dilepton production from Dalitz decay
623 c of pi0 at final time
624 *      if(nt .eq. ntmax) call dalitz_pi(nt,ntmax)
625 *                                                                      *
626 **********************************
627 *                Lables of collision channels                             *
628 **********************************
629 *         LCOLL   - NUMBER OF COLLISIONS              (INTEGER,OUTPUT) *
630 *         LBLOC   - NUMBER OF PULI-BLOCKED COLLISIONS (INTEGER,OUTPUT) *
631 *         LCNNE   - NUMBER OF ELASTIC COLLISION       (INTEGER,OUTPUT) *
632 *         LCNND   - NUMBER OF N+N->N+DELTA REACTION   (INTEGER,OUTPUT) *
633 *         LCNDN   - NUMBER OF N+DELTA->N+N REACTION   (INTEGER,OUTPUT) *
634 *         LDD     - NUMBER OF RESONANCE+RESONANCE COLLISIONS
635 *         LPP     - NUMBER OF PION+PION elastic COLIISIONS
636 *         lppk    - number of pion(RHO,OMEGA)+pion(RHO,OMEGA)
637 *                 -->K+K- collisions
638 *         LPN     - NUMBER OF PION+N-->KAON+X
639 *         lpd     - number of pion+n-->delta+pion
640 *         lrho    - number of pion+n-->Delta+rho
641 *         lomega  - number of pion+n-->Delta+omega
642 *         LKN     - NUMBER OF KAON RESCATTERINGS
643 *         LNNK    - NUMBER OF bb-->kAON PROCESS
644 *         LDDK    - NUMBER OF DD-->KAON PROCESS
645 *         LNDK    - NUMBER OF ND-->KAON PROCESS
646 ***********************************
647 * TIME-INTEGRATED COLLISIONS NUMBERS OF VARIOUS PROCESSES
648           RCOLL = RCOLL + FLOAT(LCOLL)/num
649           RBLOC = RBLOC + FLOAT(LBLOC)/num
650           RCNNE = RCNNE + FLOAT(LCNNE)/num
651          RDD   = RDD   + FLOAT(LDD)/num
652          RPP   = RPP   + FLOAT(LPP)/NUM
653          rppk  =rppk   + float(lppk)/num
654          RPN   = RPN   + FLOAT(LPN)/NUM
655          rpd   =rpd    + float(lpd)/num
656          RKN   = RKN   + FLOAT(LKN)/NUM
657          RNNK  =RNNK   + FLOAT(LNNK)/NUM
658          RDDK  =RDDK   + FLOAT(LDDK)/NUM
659          RNDK  =RNDK   + FLOAT(LNDK)/NUM
660           RCNND = RCNND + FLOAT(LCNND)/num
661           RCNDN = RCNDN + FLOAT(LCNDN)/num
662           RDIRT = RDIRT + FLOAT(LDIRT)/num
663           RDECAY= RDECAY+ FLOAT(LDECAY)/num
664           RRES  = RRES  + FLOAT(LRES)/num
665 * AVERAGE RATES OF VARIOUS COLLISIONS IN THE CURRENT TIME STEP
666           ADIRT=LDIRT/DT/num
667           ACOLL=(LCOLL-LBLOC)/DT/num
668           ACNND=LCNND/DT/num
669           ACNDN=LCNDN/DT/num
670           ADECAY=LDECAY/DT/num
671           ARES=LRES/DT/num
672          ADOU=LDOU/DT/NUM
673          ADDRHO=LDDRHO/DT/NUM
674          ANNRHO=LNNRHO/DT/NUM
675          ANNOM=LNNOM/DT/NUM
676          ADD=LDD/DT/num
677          APP=LPP/DT/num
678          appk=lppk/dt/num
679           APN=LPN/DT/num
680          apd=lpd/dt/num
681          arh=lrho/dt/num
682          aom=lomega/dt/num
683          AKN=LKN/DT/num
684          ANNK=LNNK/DT/num
685          ADDK=LDDK/DT/num
686          ANDK=LNDK/DT/num
687 * PRINT OUT THE VARIOUS COLLISION RATES
688 * (1)N-N COLLISIONS 
689 c       WRITE(1010,9991)NT*DT,ACNND,ADOU,ADIRT,ADDRHO,ANNRHO+ANNOM
690 c9991       FORMAT(6(E10.3,2X))
691 * (2)PION-N COLLISIONS
692 c       WRITE(1011,'(5(E10.3,2X))')NT*DT,apd,ARH,AOM,APN
693 * (3)KAON PRODUCTION CHANNELS
694 c        WRITE(1012,9993)NT*DT,ANNK,ADDK,ANDK,APN,Appk
695 * (4)D(N*)+D(N*) COLLISION
696 c       WRITE(1013,'(4(E10.3,2X))')NT*DT,ADDK,ADD,ADD+ADDK
697 * (5)MESON+MESON
698 c       WRITE(1014,'(4(E10.3,2X))')NT*DT,APPK,APP,APP+APPK
699 * (6)DECAY AND RESONANCE
700 c       WRITE(1016,'(3(E10.3,2X))')NT*DT,ARES,ADECAY
701 * (7)N+D(N*)
702 c       WRITE(1017,'(4(E10.3,2X))')NT*DT,ACNDN,ANDK,ACNDN+ANDK
703 c9992    FORMAT(5(E10.3,2X))
704 c9993    FORMAT(6(E10.3,2X))
705 * PRINT OUT TIME-INTEGRATED COLLISION INFORMATION
706 cbz12/28/98
707 c        write(1018,'(5(e10.3,2x),/, 4(e10.3,2x))')
708 c     &           RCNNE,RCNND,RCNDN,RDIRT,rpd,
709 c     &           RDECAY,RRES,RDD,RPP
710 c        write(1018,'(6(e10.3,2x),/, 5(e10.3,2x))')
711 c     &           NT*DT,RCNNE,RCNND,RCNDN,RDIRT,rpd,
712 c     &           NT*DT,RDECAY,RRES,RDD,RPP
713 cbz12/18/98end
714 * PRINT OUT TIME-INTEGRATED KAON MULTIPLICITIES FROM DIFFERENT CHANNELS
715 c       WRITE(1019,'(7(E10.3,2X))')NT*DT,RNNK,RDDK,RNDK,RPN,Rppk,
716 c     &                           RNNK+RDDK+RNDK+RPN+Rppk
717 *                                                                      *
718
719         END IF
720 *
721 *       UPDATE BARYON DENSITY
722 *
723         CALL DENS(IPOT,MASS,NUM,OUTPAR)
724 *
725 *       UPDATE POSITIONS FOR ALL THE PARTICLES PRESENT AT THIS TIME
726 *
727        sumene=0
728         ISO=0
729         DO 201 MRUN=1,NUM
730         ISO=ISO+MASSR(MRUN-1)
731         DO 201 I0=1,MASSR(MRUN)
732         I =I0+ISO
733         ETOTAL = SQRT( E(I)**2 + P(1,I)**2 + P(2,I)**2 +P(3,I)**2 )
734        sumene=sumene+etotal
735 C for kaons, if there is a potential
736 C CALCULATE THE ENERGY OF THE KAON ACCORDING TO THE IMPULSE APPROXIMATION
737 C REFERENCE: B.A. LI AND C.M. KO, PHYS. REV. C 54 (1996) 3283. 
738          if(kpoten.ne.0.and.lb(i).eq.23)then
739              den=0.
740               IX = NINT( R(1,I) )
741               IY = NINT( R(2,I) )
742               IZ = NINT( R(3,I) )
743 clin-4/2008:
744 c       IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
745 c     & ABS(IZ) .LT. MAXZ) den=rho(ix,iy,iz)
746               IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
747      1             .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
748      2             den=rho(ix,iy,iz)
749 c         ecor=0.1973**2*0.255*kmul*4*3.14159*(1.+0.4396/0.938)
750 c         etotal=sqrt(etotal**2+ecor*den)
751 c** G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV, m^*=m
752 c     GeV^2 fm^3
753           akg = 0.1727
754 c     GeV fm^3
755           bkg = 0.333
756          rnsg = den
757          ecor = - akg*rnsg + (bkg*den)**2
758          etotal = sqrt(etotal**2 + ecor)
759          endif
760 c
761          if(kpoten.ne.0.and.lb(i).eq.21)then
762              den=0.
763               IX = NINT( R(1,I) )
764               IY = NINT( R(2,I) )
765               IZ = NINT( R(3,I) )
766 clin-4/2008:
767 c       IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
768 c     & ABS(IZ) .LT. MAXZ) den=rho(ix,iy,iz)
769               IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
770      1             .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ)
771      2             den=rho(ix,iy,iz)
772 c* for song potential no effect on position
773 c** G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV, m^*=m
774 c     GeV^2 fm^3
775           akg = 0.1727
776 c     GeV fm^3
777           bkg = 0.333
778          rnsg = den
779          ecor = - akg*rnsg + (bkg*den)**2
780          etotal = sqrt(etotal**2 + ecor)
781           endif
782 c
783 C UPDATE POSITIONS
784           R(1,I) = R(1,I) + DT*P(1,I)/ETOTAL
785           R(2,I) = R(2,I) + DT*P(2,I)/ETOTAL
786           R(3,I) = R(3,I) + DT*P(3,I)/ETOTAL
787 c use cyclic boundary conitions
788             if ( cycbox.ne.0 ) then
789               if ( r(1,i).gt. cycbox/2 ) r(1,i)=r(1,i)-cycbox
790               if ( r(1,i).le.-cycbox/2 ) r(1,i)=r(1,i)+cycbox
791               if ( r(2,i).gt. cycbox/2 ) r(2,i)=r(2,i)-cycbox
792               if ( r(2,i).le.-cycbox/2 ) r(2,i)=r(2,i)+cycbox
793               if ( r(3,i).gt. cycbox/2 ) r(3,i)=r(3,i)-cycbox
794               if ( r(3,i).le.-cycbox/2 ) r(3,i)=r(3,i)+cycbox
795             end if
796 * UPDATE THE DELTA, N* AND PION COUNTERS
797           LB1=LB(I)
798 * 1. FOR DELTA++
799         IF(LB1.EQ.9)LD1=LD1+1
800 * 2. FOR DELTA+
801         IF(LB1.EQ.8)LD2=LD2+1
802 * 3. FOR DELTA0
803         IF(LB1.EQ.7)LD3=LD3+1
804 * 4. FOR DELTA-
805         IF(LB1.EQ.6)LD4=LD4+1
806 * 5. FOR N*+(1440)
807         IF(LB1.EQ.11)LN1=LN1+1
808 * 6. FOR N*0(1440)
809         IF(LB1.EQ.10)LN2=LN2+1
810 * 6.1 FOR N*(1535)
811        IF((LB1.EQ.13).OR.(LB1.EQ.12))LN5=LN5+1
812 * 6.2 FOR ETA
813        IF(LB1.EQ.0)LE=LE+1
814 * 6.3 FOR KAONS
815        IF(LB1.EQ.23)LKAON=LKAON+1
816 clin-11/09/00: FOR KAON*
817        IF(LB1.EQ.30)LKAONS=LKAONS+1
818
819 * UPDATE PION COUNTER
820 * 7. FOR PION+
821         IF(LB1.EQ.5)LP1=LP1+1
822 * 8. FOR PION0
823         IF(LB1.EQ.4)LP2=LP2+1
824 * 9. FOR PION-
825         IF(LB1.EQ.3)LP3=LP3+1
826 201     CONTINUE
827         LP=LP1+LP2+LP3
828         LD=LD1+LD2+LD3+LD4
829         LN=LN1+LN2
830         ALP=FLOAT(LP)/FLOAT(NUM)
831         ALD=FLOAT(LD)/FLOAT(NUM)
832         ALN=FLOAT(LN)/FLOAT(NUM)
833        ALN5=FLOAT(LN5)/FLOAT(NUM)
834         ATOTAL=ALP+ALD+ALN+0.5*ALN5
835        ALE=FLOAT(LE)/FLOAT(NUM)
836        ALKAON=FLOAT(LKAON)/FLOAT(NUM)
837 * UPDATE MOMENTUM DUE TO COULOMB INTERACTION 
838         if (icou .eq. 1) then
839 *       with Coulomb interaction
840           iso=0
841           do 1026 irun = 1,num
842             iso=iso+massr(irun-1)
843             do 1021 il = 1,massr(irun)
844                temp(1,il) = 0.
845                temp(2,il) = 0.
846                temp(3,il) = 0.
847  1021       continue
848             do 1023 il = 1, massr(irun)
849               i=iso+il
850               if (zet(lb(i)).ne.0) then
851                 do 1022 jl = 1,il-1
852                   j=iso+jl
853                   if (zet(lb(j)).ne.0) then
854                     ddx=r(1,i)-r(1,j)
855                     ddy=r(2,i)-r(2,j)
856                     ddz=r(3,i)-r(3,j)
857                     rdiff = sqrt(ddx**2+ddy**2+ddz**2)
858                     if (rdiff .le. 1.) rdiff = 1.
859                     grp=zet(lb(i))*zet(lb(j))/rdiff**3
860                     ddx=ddx*grp
861                     ddy=ddy*grp
862                     ddz=ddz*grp
863                     temp(1,il)=temp(1,il)+ddx
864                     temp(2,il)=temp(2,il)+ddy
865                     temp(3,il)=temp(3,il)+ddz
866                     temp(1,jl)=temp(1,jl)-ddx
867                     temp(2,jl)=temp(2,jl)-ddy
868                     temp(3,jl)=temp(3,jl)-ddz
869                   end if
870  1022          continue
871               end if
872  1023      continue
873             do 1025 il = 1,massr(irun)
874               i= iso+il
875               if (zet(lb(i)).ne.0) then
876                 do 1024 idir = 1,3
877                   p(idir,i) = p(idir,i) + temp(idir,il)
878      &                                    * dt * 0.00144
879  1024          continue
880               end if
881  1025      continue
882  1026   continue
883         end if
884 *       In the following, we shall:  
885 *       (1) UPDATE MOMENTA DUE TO THE MEAN FIELD FOR BARYONS AND KAONS,
886 *       (2) calculate the thermalization, temperature in a sphere of 
887 *           radius 2.0 fm AROUND THE CM
888 *       (3) AND CALCULATE THE NUMBER OF PARTICLES IN THE HIGH DENSITY REGION 
889        spt=0
890        spz=0
891        ncen=0
892        ekin=0
893           NLOST = 0
894           MEAN=0
895          nquark=0
896          nbaryn=0
897 csp06/18/01
898            rads = 2.
899            zras = 0.1
900            denst = 0.
901            edenst = 0.
902 csp06/18/01 end
903           DO 6000 IRUN = 1,NUM
904           MEAN=MEAN+MASSR(IRUN-1)
905           DO 5800 J = 1,MASSR(irun)
906           I=J+MEAN
907 c
908 csp06/18/01
909            radut = sqrt(r(1,i)**2+r(2,i)**2)
910        if( radut .le. rads )then
911         if( abs(r(3,i)) .le. zras*nt*dt )then
912 c         vols = 3.14159*radut**2*abs(r(3,i))      ! cylinder pi*r^2*l
913 c     cylinder pi*r^2*l
914          vols = 3.14159*rads**2*zras
915          engs=sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)
916          gammas=1.
917          if(e(i).ne.0.)gammas=engs/e(i)
918 c     rho
919          denst = denst + 1./gammas/vols
920 c     energy density
921          edenst = edenst + engs/gammas/gammas/vols
922         endif
923        endif
924 csp06/18/01 end
925 c
926          drr=sqrt(r(1,i)**2+r(2,i)**2+r(3,i)**2)
927          if(drr.le.2.0)then
928          spt=spt+p(1,i)**2+p(2,i)**2
929          spz=spz+p(3,i)**2
930          ncen=ncen+1
931          ekin=ekin+sqrt(p(1,i)**2+p(2,i)**2+p(3,i)**2+e(i)**2)-e(i)
932          endif
933               IX = NINT( R(1,I) )
934               IY = NINT( R(2,I) )
935               IZ = NINT( R(3,I) )
936 C calculate the No. of particles in the high density region
937 clin-4/2008:
938 c              IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
939 c     & ABS(IZ) .LT. MAXZ) THEN
940               IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
941      1          .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
942        if(rho(ix,iy,iz)/0.168.gt.dencut)go to 5800
943        if((rho(ix,iy,iz)/0.168.gt.5.).and.(e(i).gt.0.9))
944      &  nbaryn=nbaryn+1
945        if(pel(ix,iy,iz).gt.2.0)nquark=nquark+1
946        endif
947 c*
948 c If there is a kaon potential, propogating kaons 
949         if(kpoten.ne.0.and.lb(i).eq.23)then
950         den=0.
951 clin-4/2008:
952 c       IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
953 c     & ABS(IZ) .LT. MAXZ)then
954         IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
955      1       .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
956            den=rho(ix,iy,iz)
957 c        ecor=0.1973**2*0.255*kmul*4*3.14159*(1.+0.4396/0.938)
958 c       etotal=sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2+ecor*den)
959 c** for G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV
960 c     !! GeV^2 fm^3
961             akg = 0.1727
962 c     !! GeV fm^3
963             bkg = 0.333
964           rnsg = den
965           ecor = - akg*rnsg + (bkg*den)**2
966           etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
967           ecor = - akg + 2.*bkg**2*den + 2.*bkg*etotal
968 c** G.Q. Li potential (END)           
969         CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
970         P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
971         P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
972         P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
973         endif
974          endif
975 c
976         if(kpoten.ne.0.and.lb(i).eq.21)then
977          den=0.
978 clin-4/2008:
979 c           IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
980 c     &        ABS(IZ) .LT. MAXZ)then
981          IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
982      1        .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
983                den=rho(ix,iy,iz)
984         CALL GRADUK(IX,IY,IZ,GRADXk,GRADYk,GRADZk)
985 c        P(1,I) = P(1,I) - DT * GRADXk*(-0.12/0.168)    !! song potential
986 c        P(2,I) = P(2,I) - DT * GRADYk*(-0.12/0.168)
987 c        P(3,I) = P(3,I) - DT * GRADZk*(-0.12/0.168)
988 c** for G.Q Li potential form with n_s = n_b and pot(n_0)=29 MeV
989 c    !! GeV^2 fm^3
990             akg = 0.1727
991 c     !! GeV fm^3
992             bkg = 0.333
993           rnsg = den
994           ecor = - akg*rnsg + (bkg*den)**2
995           etotal = sqrt(P(1,i)**2+p(2,I)**2+p(3,i)**2+e(i)**2 + ecor)
996           ecor = - akg + 2.*bkg**2*den - 2.*bkg*etotal
997         P(1,I) = P(1,I) - DT * GRADXk*ecor/(2.*etotal)
998         P(2,I) = P(2,I) - DT * GRADYk*ecor/(2.*etotal)
999         P(3,I) = P(3,I) - DT * GRADZk*ecor/(2.*etotal)
1000 c** G.Q. Li potential (END)           
1001         endif
1002          endif
1003 c
1004 c for other mesons, there is no potential
1005        if(j.gt.mass)go to 5800         
1006 c  with mean field interaction for baryons   (open endif below) !!sp05
1007 **      if( (iabs(lb(i)).eq.1.or.iabs(lb(i)).eq.2) .or.
1008 **    &     (iabs(lb(i)).ge.6.and.iabs(lb(i)).le.17) .or.
1009 **    &      iabs(lb(i)).eq.40.or.iabs(lb(i)).eq.41 )then  
1010         IF (ICOLL .NE. -1) THEN
1011 * check if the baryon has run off the lattice
1012 *             IX0=NINT(R(1,I)/DX)
1013 *             IY0=NINT(R(2,I)/DY)
1014 *             IZ0=NINT(R(3,I)/DZ)
1015 *             IPX0=NINT(P(1,I)/DPX)
1016 *             IPY0=NINT(P(2,I)/DPY)
1017 *             IPZ0=NINT(P(3,I)/DPZ)
1018 *      if ( (abs(ix0).gt.mx) .or. (abs(iy0).gt.my) .or. (abs(iz0).gt.mz)
1019 *     &  .or. (abs(ipx0).gt.mpx) .or. (abs(ipy0) 
1020 *     &  .or. (ipz0.lt.-mpz) .or. (ipz0.gt.mpzp)) NLOST=NLOST+1
1021 clin-4/2008:
1022 c              IF (ABS(IX) .LT. MAXX .AND. ABS(IY) .LT. MAXX .AND.
1023 c     &                                    ABS(IZ) .LT. MAXZ     ) THEN
1024            IF(IX.LT.MAXX.AND.IY.LT.MAXX.AND.IZ.LT.MAXZ
1025      1          .AND.IX.GT.-MAXX.AND.IY.GT.-MAXX.AND.IZ.GT.-MAXZ) THEN
1026                 CALL GRADU(IPOT,IX,IY,IZ,GRADX,GRADY,GRADZ)
1027               TZ=0.
1028               GRADXN=0
1029               GRADYN=0
1030               GRADZN=0
1031               GRADXP=0
1032               GRADYP=0
1033               GRADZP=0
1034              IF(ICOU.EQ.1)THEN
1035                 CALL GRADUP(IX,IY,IZ,GRADXP,GRADYP,GRADZP)
1036                 CALL GRADUN(IX,IY,IZ,GRADXN,GRADYN,GRADZN)
1037                IF(ZET(LB(I)).NE.0)TZ=-1
1038                IF(ZET(LB(I)).EQ.0)TZ= 1
1039              END IF
1040            if(iabs(lb(i)).ge.14.and.iabs(lb(i)).le.17)then
1041               facl = 2./3.
1042             elseif(iabs(lb(i)).eq.40.or.iabs(lb(i)).eq.41)then
1043               facl = 1./3.
1044             else
1045               facl = 1.
1046             endif
1047         P(1,I) = P(1,I) - facl*DT * (GRADX+asy*(GRADXN-GRADXP)*TZ)
1048         P(2,I) = P(2,I) - facl*DT * (GRADY+asy*(GRADYN-GRADYP)*TZ)
1049         P(3,I) = P(3,I) - facl*DT * (GRADZ+asy*(GRADZN-GRADZP)*TZ)
1050                 end if                                                       
1051               ENDIF
1052 **          endif          !!sp05     
1053  5800       CONTINUE
1054  6000       CONTINUE
1055 c print out the average no. of particles in regions where the local 
1056 c baryon density is higher than 5*rho0 
1057 c       write(1072,'(e10.3,2x,e10.3)')nt*dt,float(nbaryn)/float(num)
1058 C print out the average no. of particles in regions where the local 
1059 c energy density is higher than 2 GeV/fm^3. 
1060 c       write(1073,'(e10.3,2x,e10.3)')nt*dt,float(nquark)/float(num)
1061 c print out the no. of particles that have run off the lattice
1062 *          IF (NLOST .NE. 0 .AND. (NT/NFREQ)*NFREQ .EQ. NT) THEN
1063 *            WRITE(12,'(5X,''***'',I7,'' TESTPARTICLES LOST AFTER '',
1064 *     &                   ''TIME STEP NUMBER'',I4)') NLOST, NT
1065 *         END IF
1066 *
1067 *       update phase space density
1068 *        call platin(mode,mass,num,dx,dy,dz,dpx,dpy,dpz,fnorm)
1069 *
1070 *       CONTROL-PRINTOUT OF CONFIGURATION (IF REQUIRED)
1071 *
1072 *        if (inout(5) .eq. 2) CALL ENERGY(NT,IPOT,NUM,MASS,EMIN,EMAX)
1073 *
1074
1075 * print out central baryon density as a function of time
1076        CDEN=RHO(0,0,0)/0.168
1077 cc        WRITE(1002,990)FLOAT(NT)*DT,CDEN
1078 c        WRITE(1002,1990)FLOAT(NT)*DT,CDEN,denst/real(num)
1079 * print out the central energy density as a function of time
1080 cc        WRITE(1003,990)FLOAT(NT)*DT,PEL(0,0,0)
1081 c        WRITE(1003,1990)FLOAT(NT)*DT,PEL(0,0,0),edenst/real(num)
1082 * print out the no. of pion-like particles as a function of time 
1083 c        WRITE(1004,9999)FLOAT(NT)*DT,ALD,ALN,ALP,ALN5,
1084 c     &               ALD+ALN+ALP+0.5*ALN5
1085 * print out the no. of eta-like particles as a function of time
1086 c        WRITE(1005,991)FLOAT(NT)*DT,ALN5,ALE,ALE+0.5*ALN5
1087 c990       FORMAT(E10.3,2X,E10.3)
1088 c1990       FORMAT(E10.3,2X,E10.3,2X,E10.3)
1089 c991       FORMAT(E10.3,2X,E10.3,2X,E10.3,2X,E10.3)
1090 c9999    FORMAT(e10.3,2X,e10.3,2X,E10.3,2X,E10.3,2X,
1091 c     1  E10.3,2X,E10.3)
1092 C THE FOLLOWING OUTPUTS CAN BE TURNED ON/OFF by setting icflow and icrho=0  
1093 c print out the baryon and meson density matrix in the reaction plane
1094         IF ((NT/NFREQ)*NFREQ .EQ. NT ) THEN
1095        if(icflow.eq.1)call flow(nt)
1096 cbz11/18/98
1097 c       if(icrho.ne.1)go to 10000 
1098 c       if (icrho .eq. 1) then 
1099 cbz11/18/98end
1100 c       do ix=-10,10
1101 c       do iz=-10,10
1102 c       write(1053,992)ix,iz,rho(ix,0,iz)/0.168
1103 c       write(1054,992)ix,iz,pirho(ix,0,iz)/0.168
1104 c       write(1055,992)ix,iz,pel(ix,0,iz)
1105 c       end do
1106 c       end do
1107 cbz11/18/98
1108 c        end if
1109 cbz11/18/98end
1110 c992       format(i3,i3,e11.4)
1111        endif
1112 c print out the ENERGY density matrix in the reaction plane
1113 C CHECK LOCAL MOMENTUM EQUILIBRIUM IN EACH CELL, 
1114 C AND PERFORM ON-LINE FLOW ANALYSIS AT A FREQUENCY OF NFREQ
1115 c        IF ((NT/NFREQ)*NFREQ .EQ. NT ) THEN
1116 c       call flow(nt)
1117 c       call equ(ipot,mass,num,outpar)
1118 c       do ix=-10,10
1119 c       do iz=-10,10
1120 c       write(1055,992)ix,iz,pel(ix,0,iz)
1121 c       write(1056,992)ix,iz,rxy(ix,0,iz)
1122 c       end do
1123 c       end do
1124 c       endif
1125 C calculate the volume of high BARYON AND ENERGY density 
1126 C matter as a function of time
1127 c       vbrho=0.
1128 c       verho=0.
1129 c       do ix=-20,20
1130 c       do iy=-20,20
1131 c       do iz=-20,20
1132 c       if(rho(ix,iy,iz)/0.168.gt.5.)vbrho=vbrho+1.
1133 c       if(pel(ix,iy,iz).gt.2.)verho=verho+1.
1134 c       end do
1135 c       end do
1136 c       end do
1137 c       write(1081,993)dt*nt,vbrho
1138 c       write(1082,993)dt*nt,verho
1139 c993       format(e11.4,2x,e11.4)
1140 *-----------------------------------------------------------------------
1141 cbz11/16/98
1142 c.....for read-in initial conditions produce particles from read-in 
1143 c.....common block.
1144 c.....note that this part is only for cascade with number of test particles
1145 c.....NUM = 1.
1146       IF (IAPAR2(1) .NE. 1) THEN
1147          CT = NT * DT
1148 cbz12/22/98
1149 c         NP = MASSR(1)
1150 c         DO WHILE (FTAR(NPI) .GT. CT - DT .AND. FTAR(NPI) .LE. CT)
1151 c            NP = NP + 1
1152 c            R(1, NP) = GXAR(NPI) + PXAR(NPI) / PEAR(NPI) * (CT - FTAR(NPI))
1153 c            R(2, NP) = GYAR(NPI) + PYAR(NPI) / PEAR(NPI) * (CT - FTAR(NPI))
1154 c            R(3, NP) = GZAR(NPI) + PZAR(NPI) / PEAR(NPI) * (CT - FTAR(NPI))
1155 c            P(1, NP) = PXAR(NPI)
1156 c            P(2, NP) = PYAR(NPI)
1157 c            P(3, NP) = PZAR(NPI)
1158 c            E(NP) = XMAR(NPI)
1159 c            LB(NP) = IARFLV(ITYPAR(NPI))
1160 c            NPI = NPI + 1
1161 c         END DO
1162 c         MASSR(1) = NP
1163          IA = 0
1164          DO 1028 IRUN = 1, NUM
1165             DO 1027 IC = 1, MASSR(IRUN)
1166                IE = IA + IC
1167                RT(1, IC, IRUN) = R(1, IE)
1168                RT(2, IC, IRUN) = R(2, IE)
1169                RT(3, IC, IRUN) = R(3, IE)
1170                PT(1, IC, IRUN) = P(1, IE)
1171                PT(2, IC, IRUN) = P(2, IE)
1172                PT(3, IC, IRUN) = P(3, IE)
1173                ET(IC, IRUN) = E(IE)
1174                LT(IC, IRUN) = LB(IE)
1175 c         !! sp 12/19/00
1176                PROT(IC, IRUN) = PROPER(IE)
1177 clin-5/2008:
1178                dpertt(IC, IRUN)=dpertp(IE)
1179  1027       CONTINUE
1180             NP = MASSR(IRUN)
1181             NP1 = NPI(IRUN)
1182
1183 cbz10/05/99
1184 c            DO WHILE (FT1(NP1, IRUN) .GT. CT - DT .AND. 
1185 c     &           FT1(NP1, IRUN) .LE. CT)
1186 cbz10/06/99
1187 c            DO WHILE (NPI(IRUN).LE.MULTI1(IRUN).AND.
1188 cbz10/06/99 end
1189 clin-11/13/00 finally read in all unformed particles and do the decays in ART:
1190 c           DO WHILE (NP1.LE.MULTI1(IRUN).AND.
1191 c    &           FT1(NP1, IRUN) .GT. CT - DT .AND. 
1192 c    &           FT1(NP1, IRUN) .LE. CT)
1193 c
1194                ctlong = ct
1195              if(nt .eq. (ntmax-1))then
1196                ctlong = 1.E30
1197              elseif(nt .eq. ntmax)then
1198                go to 1111
1199              endif
1200             DO WHILE (NP1.LE.MULTI1(IRUN).AND.
1201      &           FT1(NP1, IRUN) .GT. (CT - DT) .AND. 
1202      &           FT1(NP1, IRUN) .LE. ctlong)
1203                NP = NP + 1
1204                UDT = (CT - FT1(NP1, IRUN)) / EE1(NP1, IRUN)
1205 clin-10/28/03 since all unformed hadrons at time ct are read in at nt=ntmax-1, 
1206 c     their positions should not be propagated to time ct:
1207                if(nt.eq.(ntmax-1)) then
1208                   ftsvt(NP,IRUN)=FT1(NP1, IRUN)
1209                   if(FT1(NP1, IRUN).gt.ct) UDT=0.
1210                endif
1211                RT(1, NP, IRUN) = GX1(NP1, IRUN) + 
1212      &              PX1(NP1, IRUN) * UDT
1213                RT(2, NP, IRUN) = GY1(NP1, IRUN) + 
1214      &              PY1(NP1, IRUN) * UDT
1215                RT(3, NP, IRUN) = GZ1(NP1, IRUN) + 
1216      &              PZ1(NP1, IRUN) * UDT
1217                PT(1, NP, IRUN) = PX1(NP1, IRUN)
1218                PT(2, NP, IRUN) = PY1(NP1, IRUN)
1219                PT(3, NP, IRUN) = PZ1(NP1, IRUN)
1220                ET(NP, IRUN) = XM1(NP1, IRUN)
1221                LT(NP, IRUN) = IARFLV(ITYP1(NP1, IRUN))
1222 clin-5/2008:
1223                dpertt(NP,IRUN)=dpp1(NP1,IRUN)
1224 clin-4/30/03 ctest off 
1225 c     record initial phi,K*,Lambda(1520) resonances formed during the timestep:
1226 c               if(LT(NP, IRUN).eq.29.or.iabs(LT(NP, IRUN)).eq.30)
1227 c     1              write(17,112) 'formed',LT(NP, IRUN),PX1(NP1, IRUN),
1228 c     2 PY1(NP1, IRUN),PZ1(NP1, IRUN),XM1(NP1, IRUN),nt
1229 c 112           format(a10,1x,I4,4(1x,f9.3),1x,I4)
1230 c
1231                NP1 = NP1 + 1
1232 c     !! sp 12/19/00
1233                PROT(NP, IRUN) = 1.
1234             END DO
1235 *
1236  1111      continue
1237             NPI(IRUN) = NP1
1238             IA = IA + MASSR(IRUN)
1239             MASSR(IRUN) = NP
1240  1028    CONTINUE
1241          IA = 0
1242          DO 1030 IRUN = 1, NUM
1243             IA = IA + MASSR(IRUN - 1)
1244             DO 1029 IC = 1, MASSR(IRUN)
1245                IE = IA + IC
1246                R(1, IE) = RT(1, IC, IRUN)
1247                R(2, IE) = RT(2, IC, IRUN)
1248                R(3, IE) = RT(3, IC, IRUN)
1249                P(1, IE) = PT(1, IC, IRUN)
1250                P(2, IE) = PT(2, IC, IRUN)
1251                P(3, IE) = PT(3, IC, IRUN)
1252                E(IE) = ET(IC, IRUN)
1253                LB(IE) = LT(IC, IRUN)
1254 c     !! sp 12/19/00
1255                PROPER(IE) = PROT(IC, IRUN)
1256                if(nt.eq.(ntmax-1)) ftsv(IE)=ftsvt(IC,IRUN)
1257 clin-5/2008:
1258                dpertp(IE)=dpertt(IC, IRUN)
1259  1029       CONTINUE
1260 clin-3/2009 Moved here to better take care of freezeout spacetime:
1261             call hbtout(MASSR(IRUN),nt,ntmax)
1262  1030    CONTINUE
1263 cbz12/22/98end
1264       END IF
1265 cbz11/16/98end
1266
1267 clin-5/2009 ctest off:
1268 c      call flowh(ct) 
1269
1270 10000       continue
1271
1272 *                                                                      *
1273 *       ==============  END OF TIME STEP LOOP   ================       *
1274
1275 ************************************
1276 *     WRITE OUT particle's MOMENTA ,and/OR COORDINATES ,
1277 *     label and/or their local baryon density in the final state
1278         iss=0
1279         do 1032 lrun=1,num
1280            iss=iss+massr(lrun-1)
1281            do 1031 l0=1,massr(lrun)
1282               ipart=iss+l0
1283  1031      continue
1284  1032   continue
1285
1286 cbz11/16/98
1287       IF (IAPAR2(1) .NE. 1) THEN
1288 cbz12/22/98
1289 c        NSH = MASSR(1) - NPI + 1
1290 c        IAINT2(1) = IAINT2(1) + NSH
1291 c.....to shift the unformed particles to the end of the common block
1292 c        IF (NSH .GT. 0) THEN
1293 c           IB = IAINT2(1)
1294 c           IE = MASSR(1) + 1
1295 c           II = -1
1296 c        ELSE IF (NSH .LT. 0) THEN
1297 c           IB = MASSR(1) + 1
1298 c           IE = IAINT2(1)
1299 c           II = 1
1300 c        END IF
1301 c        IF (NSH .NE. 0) THEN
1302 c           DO I = IB, IE, II
1303 c              J = I - NSH
1304 c              ITYPAR(I) = ITYPAR(J)
1305 c              GXAR(I) = GXAR(J)
1306 c              GYAR(I) = GYAR(J)
1307 c              GZAR(I) = GZAR(J)
1308 c              FTAR(I) = FTAR(J)
1309 c              PXAR(I) = PXAR(J)
1310 c              PYAR(I) = PYAR(J)
1311 c              PZAR(I) = PZAR(J)
1312 c              PEAR(I) = PEAR(J)
1313 c              XMAR(I) = XMAR(J)
1314 c           END DO
1315 c        END IF
1316
1317 c.....to copy ART particle info to COMMON /ARPRC/
1318 c        DO I = 1, MASSR(1)
1319 c           ITYPAR(I) = INVFLV(LB(I))
1320 c           GXAR(I) = R(1, I)
1321 c           GYAR(I) = R(2, I)
1322 c           GZAR(I) = R(3, I)
1323 c           FTAR(I) = CT
1324 c           PXAR(I) = P(1, I)
1325 c           PYAR(I) = P(2, I)
1326 c           PZAR(I) = P(3, I)
1327 c           XMAR(I) = E(I)
1328 c           PEAR(I) = SQRT(PXAR(I) ** 2 + PYAR(I) ** 2 + PZAR(I) ** 2
1329 c     &        + XMAR(I) ** 2)
1330 c        END DO
1331         IA = 0
1332         DO 1035 IRUN = 1, NUM
1333            IA = IA + MASSR(IRUN - 1)
1334            NP1 = NPI(IRUN)
1335            NSH = MASSR(IRUN) - NP1 + 1
1336            MULTI1(IRUN) = MULTI1(IRUN) + NSH
1337 c.....to shift the unformed particles to the end of the common block
1338            IF (NSH .GT. 0) THEN
1339               IB = MULTI1(IRUN)
1340               IE = MASSR(IRUN) + 1
1341               II = -1
1342            ELSE IF (NSH .LT. 0) THEN
1343               IB = MASSR(IRUN) + 1
1344               IE = MULTI1(IRUN)
1345               II = 1
1346            END IF
1347            IF (NSH .NE. 0) THEN
1348               DO 1033 I = IB, IE, II
1349                  J = I - NSH
1350                  ITYP1(I, IRUN) = ITYP1(J, IRUN)
1351                  GX1(I, IRUN) = GX1(J, IRUN)
1352                  GY1(I, IRUN) = GY1(J, IRUN)
1353                  GZ1(I, IRUN) = GZ1(J, IRUN)
1354                  FT1(I, IRUN) = FT1(J, IRUN)
1355                  PX1(I, IRUN) = PX1(J, IRUN)
1356                  PY1(I, IRUN) = PY1(J, IRUN)
1357                  PZ1(I, IRUN) = PZ1(J, IRUN)
1358                  EE1(I, IRUN) = EE1(J, IRUN)
1359                  XM1(I, IRUN) = XM1(J, IRUN)
1360 c     !! sp 12/19/00
1361                  PRO1(I, IRUN) = PRO1(J, IRUN)
1362 clin-5/2008:
1363                  dpp1(I,IRUN)=dpp1(J,IRUN)
1364  1033         CONTINUE
1365            END IF
1366            
1367 c.....to copy ART particle info to COMMON /ARPRC1/
1368            DO 1034 I = 1, MASSR(IRUN)
1369               IB = IA + I
1370               ITYP1(I, IRUN) = INVFLV(LB(IB))
1371               GX1(I, IRUN) = R(1, IB)
1372               GY1(I, IRUN) = R(2, IB)
1373               GZ1(I, IRUN) = R(3, IB)
1374 clin-10/28/03:
1375 c since all unformed hadrons at time ct are read in at nt=ntmax-1, 
1376 c their formation time ft1 should be kept to determine their freezeout(x,t):
1377 c              FT1(I, IRUN) = CT
1378               if(FT1(I, IRUN).lt.CT) FT1(I, IRUN) = CT
1379               PX1(I, IRUN) = P(1, IB)
1380               PY1(I, IRUN) = P(2, IB)
1381               PZ1(I, IRUN) = P(3, IB)
1382               XM1(I, IRUN) = E(IB)
1383               EE1(I, IRUN) = SQRT(PX1(I, IRUN) ** 2 + 
1384      &             PY1(I, IRUN) ** 2 +
1385      &             PZ1(I, IRUN) ** 2 + 
1386      &             XM1(I, IRUN) ** 2)
1387 c     !! sp 12/19/00
1388               PRO1(I, IRUN) = PROPER(IB)
1389  1034      CONTINUE
1390  1035   CONTINUE
1391 cbz12/22/98end
1392       END IF
1393 cbz11/16/98end
1394 c
1395 **********************************
1396 *                                                                      *
1397 *       ======= END OF MANY LOOPS OVER IMPACT PARAMETERS ==========    *
1398 *                                                               *
1399 **********************************
1400 50000   CONTINUE
1401 *
1402 *-----------------------------------------------------------------------
1403 *                       ==== ART COMPLETED ====
1404 *-----------------------------------------------------------------------
1405 cbz11/16/98
1406 c      STOP
1407       RETURN
1408 cbz11/16/98end
1409       END
1410 **********************************
1411       subroutine coulin(masspr,massta,NUM)
1412 *                                                                      *
1413 *     purpose:   initialization of array zet() and lb() for all runs  *
1414 *                lb(i) = 1   =>  proton                               *
1415 *                lb(i) = 2   =>  neutron                              *
1416 **********************************
1417         integer  zta,zpr
1418         PARAMETER (MAXSTR=150001)
1419         common  /EE/ ID(MAXSTR),LB(MAXSTR)
1420 cc      SAVE /EE/
1421         COMMON  /ZZ/ ZTA,ZPR
1422 cc      SAVE /zz/
1423       SAVE   
1424         MASS=MASSTA+MASSPR
1425         DO 500 IRUN=1,NUM
1426         do 100 i = 1+(IRUN-1)*MASS,zta+(IRUN-1)*MASS
1427         LB(i) = 1
1428   100   continue
1429         do 200 i = zta+1+(IRUN-1)*MASS,massta+(IRUN-1)*MASS
1430         LB(i) = 2
1431   200   continue
1432         do 300 i = massta+1+(IRUN-1)*MASS,massta+zpr+(IRUN-1)*MASS
1433         LB(i) = 1
1434   300   continue
1435         do 400 i = massta+zpr+1+(IRUN-1)*MASS,
1436      1  massta+masspr+(IRUN-1)*MASS
1437         LB(i) = 2
1438   400   continue
1439   500   CONTINUE
1440         return
1441         end
1442 **********************************
1443 *                                                                      *
1444       SUBROUTINE RELCOL(LCOLL,LBLOC,LCNNE,LDD,LPP,lppk,
1445      &LPN,lpd,lrho,lomega,LKN,LNNK,LDDK,LNDK,LCNND,LCNDN,
1446      &LDIRT,LDECAY,LRES,LDOU,LDDRHO,LNNRHO,LNNOM,
1447      &NT,ntmax,sp,akaon,sk)
1448 *                                                                      *
1449 *       PURPOSE:    CHECK CONDITIONS AND CALCULATE THE KINEMATICS      * 
1450 *                   FOR BINARY COLLISIONS AMONG PARTICLES              *
1451 *                                 - RELATIVISTIC FORMULA USED          *
1452 *                                                                      *
1453 *       REFERENCES: HAGEDORN, RELATIVISTIC KINEMATICS (1963)           *
1454 *                                                                      *
1455 *       VARIABLES:                                                     *
1456 *         MASSPR  - NUMBER OF NUCLEONS IN PROJECTILE   (INTEGER,INPUT) *
1457 *         MASSTA  - NUMBER OF NUCLEONS IN TARGET       (INTEGER,INPUT) *
1458 *         NUM     - NUMBER OF TESTPARTICLES PER NUCLEON(INTEGER,INPUT) *
1459 *         ISEED   - SEED FOR RANDOM NUMBER GENERATOR   (INTEGER,INPUT) *
1460 *         IAVOID  - (= 1 => AVOID FIRST CLLISIONS WITHIN THE SAME      *
1461 *                   NUCLEUS, ELSE ALL COLLISIONS)      (INTEGER,INPUT) *
1462 *         DELTAR  - MAXIMUM SPATIAL DISTANCE FOR WHICH A COLLISION     *
1463 *                   STILL CAN OCCUR                       (REAL,INPUT) *
1464 *         DT      - TIME STEP SIZE                        (REAL,INPUT) *
1465 *         LCOLL   - NUMBER OF COLLISIONS              (INTEGER,OUTPUT) *
1466 *         LBLOC   - NUMBER OF PULI-BLOCKED COLLISIONS (INTEGER,OUTPUT) *
1467 *         LCNNE   - NUMBER OF ELASTIC COLLISION       (INTEGER,OUTPUT) *
1468 *         LCNND   - NUMBER OF N+N->N+DELTA REACTION   (INTEGER,OUTPUT) *
1469 *         LCNDN   - NUMBER OF N+DELTA->N+N REACTION   (INTEGER,OUTPUT) *
1470 *         LDD     - NUMBER OF RESONANCE+RESONANCE COLLISIONS
1471 *         LPP     - NUMBER OF PION+PION elastic COLIISIONS
1472 *         lppk    - number of pion(RHO,OMEGA)+pion(RHO,OMEGA)
1473 *                   -->K+K- collisions
1474 *         LPN     - NUMBER OF PION+N-->KAON+X
1475 *         lpd     - number of pion+n-->delta+pion
1476 *         lrho    - number of pion+n-->Delta+rho
1477 *         lomega  - number of pion+n-->Delta+omega
1478 *         LKN     - NUMBER OF KAON RESCATTERINGS
1479 *         LNNK    - NUMBER OF bb-->kAON PROCESS
1480 *         LDDK    - NUMBER OF DD-->KAON PROCESS
1481 *         LNDK    - NUMBER OF ND-->KAON PROCESS
1482 *         LB(I) IS USED TO LABEL PARTICLE'S CHARGE STATE
1483 *         LB(I)   = 
1484 cbali2/7/99 
1485 *                 -45 Omega baryon(bar)
1486 *                 -41 cascade0(bar)
1487 *                 -40 cascade-(bar)
1488 clin-11/07/00:
1489 *                 -30 K*-
1490 *                 -17 sigma+(bar)
1491 *                 -16 sigma0(bar)
1492 *                 -15 sigma-(bar)
1493 *                 -14 LAMBDA(bar)
1494 clin-8/29/00
1495 *                 -13 anti-N*(+1)(1535),s_11
1496 *                 -12 anti-N*0(1535),s_11
1497 *                 -11 anti-N*(+1)(1440),p_11
1498 *                 -10 anti-N*0(1440), p_11
1499 *                  -9 anti-DELTA+2
1500 *                  -8 anti-DELTA+1
1501 *                  -7 anti-DELTA0
1502 *                  -6 anti-DELTA-1
1503 *
1504 *                  -2 antineutron 
1505 *                  -1 antiproton
1506 cbali2/7/99end 
1507 *                   0 eta
1508 *                   1 PROTON
1509 *                   2 NUETRON
1510 *                   3 PION-
1511 *                   4 PION0
1512 *                   5 PION+          
1513 *                   6 DELTA-1
1514 *                   7 DELTA0
1515 *                   8 DELTA+1
1516 *                   9 DELTA+2
1517 *                   10 N*0(1440), p_11
1518 *                   11 N*(+1)(1440),p_11
1519 *                  12 N*0(1535),s_11
1520 *                  13 N*(+1)(1535),s_11
1521 *                  14 LAMBDA
1522 *                   15 sigma-
1523 *                   16 sigma0
1524 *                   17 sigma+
1525 *                   21 kaon-
1526 clin-2/23/03        22 Kaon0Long (converted at the last timestep)
1527 *                   23 KAON+
1528 *                   24 Kaon0short (converted at the last timestep then decay)
1529 *                   25 rho-
1530 *                   26 rho0
1531 *                   27 rho+
1532 *                   28 omega meson
1533 *                   29 phi
1534 *                   30 K*+
1535 * sp01/03/01
1536 *                   31 eta-prime
1537 *                   40 cascade-
1538 *                   41 cascade0
1539 *                   45 Omega baryon
1540 * sp01/03/01 end
1541 *                   
1542 *                   ++  ------- SEE NOTE BOOK
1543 *         NSTAR=1 INCLUDING N* RESORANCE
1544 *         ELSE DELTA RESORANCE ONLY
1545 *         NDIRCT=1 INCLUDING DIRECT PROCESS,ELSE NOT
1546 *         DIR - PERCENTAGE OF DIRECT PION PRODUCTION PROCESS
1547 **********************************
1548       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
1549       parameter      (MX=4,MY=4,MZ=8,MPX=4,MPY=4,mpz=10,mpzp=10)
1550       PARAMETER      (AKA=0.498,ALA=1.1157,ASA=1.1974,aks=0.895)
1551       PARAMETER      (AA1=1.26,APHI=1.02,AP1=0.13496)
1552       parameter            (maxx=20,maxz=24)
1553       parameter            (rrkk=0.6,prkk=0.3,srhoks=5.,ESBIN=0.04)
1554       DIMENSION MASSRN(0:MAXR),RT(3,MAXSTR),PT(3,MAXSTR),ET(MAXSTR)
1555       DIMENSION LT(MAXSTR), PROT(MAXSTR)
1556       COMMON   /AA/  R(3,MAXSTR)
1557 cc      SAVE /AA/
1558       COMMON   /BB/  P(3,MAXSTR)
1559 cc      SAVE /BB/
1560       COMMON   /CC/  E(MAXSTR)
1561 cc      SAVE /CC/
1562       COMMON  /DD/      RHO(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
1563      &                     RHOP(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ),
1564      &                     RHON(-MAXX:MAXX,-MAXX:MAXX,-MAXZ:MAXZ)
1565 cc      SAVE /DD/
1566       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
1567 cc      SAVE /EE/
1568       COMMON   /HH/  PROPER(MAXSTR)
1569 cc      SAVE /HH/
1570       common /ff/f(-mx:mx,-my:my,-mz:mz,-mpx:mpx,-mpy:mpy,-mpz:mpzp)
1571 cc      SAVE /ff/
1572       common   /gg/  dx,dy,dz,dpx,dpy,dpz
1573 cc      SAVE /gg/
1574       COMMON   /INPUT/ NSTAR,NDIRCT,DIR
1575 cc      SAVE /INPUT/
1576       COMMON   /NN/NNN
1577 cc      SAVE /NN/
1578       COMMON   /RR/  MASSR(0:MAXR)
1579 cc      SAVE /RR/
1580       common   /ss/  inout(20)
1581 cc      SAVE /ss/
1582       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
1583 cc      SAVE /BG/
1584       COMMON   /RUN/NUM
1585 cc      SAVE /RUN/
1586       COMMON   /PA/RPION(3,MAXSTR,MAXR)
1587 cc      SAVE /PA/
1588       COMMON   /PB/PPION(3,MAXSTR,MAXR)
1589 cc      SAVE /PB/
1590       COMMON   /PC/EPION(MAXSTR,MAXR)
1591 cc      SAVE /PC/
1592       COMMON   /PD/LPION(MAXSTR,MAXR)
1593 cc      SAVE /PD/
1594       COMMON   /PE/PROPI(MAXSTR,MAXR)
1595 cc      SAVE /PE/
1596       COMMON   /KKK/TKAON(7),EKAON(7,0:2000)
1597 cc      SAVE /KKK/
1598       COMMON  /KAON/    AK(3,50,36),SPECK(50,36,7),MF
1599 cc      SAVE /KAON/
1600       COMMON/TABLE/ xarray(0:1000),earray(0:1000)
1601 cc      SAVE /TABLE/
1602       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1603 cc      SAVE /input1/
1604       common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
1605      1 px1n,py1n,pz1n,dp1n
1606 cc      SAVE /leadng/
1607       COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
1608 cc      SAVE /tdecay/
1609       common /lastt/itimeh,bimp 
1610 cc      SAVE /lastt/
1611 c
1612       COMMON/ppbmas/niso(15),nstate,ppbm(15,2),thresh(15),weight(15)
1613 cc      SAVE /ppbmas/
1614       common/ppb1/ene,factr2(6),fsum,ppinnb,s,wtot
1615 cc      SAVE /ppb1/
1616       common/ppmm/pprr,ppee,pppe,rpre,xopoe,rree
1617 cc      SAVE /ppmm/
1618       COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
1619 cc      SAVE /hbt/
1620       common/resdcy/NSAV,iksdcy
1621 cc      SAVE /resdcy/
1622       COMMON/RNDF77/NSEED
1623 cc      SAVE /RNDF77/
1624       COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
1625       dimension ftpisv(MAXSTR,MAXR),fttemp(MAXSTR)
1626       common /dpi/em2,lb2
1627       common/phidcy/iphidcy,pttrig,ntrig,maxmiss
1628 clin-5/2008:
1629       DIMENSION dptemp(MAXSTR)
1630       common /para8/ idpert,npertd,idxsec
1631       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1632      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1633      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1634 c
1635       real zet(-45:45)
1636       SAVE   
1637       data zet /
1638      4     1.,0.,0.,0.,0.,
1639      3     1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1640      2     -1.,0.,0.,0.,0.,0.,0.,0.,0.,0.,
1641      1     0.,0.,0.,-1.,0.,1.,0.,-1.,0.,-1.,
1642      s     0.,-2.,-1.,0.,1.,0.,0.,0.,0.,-1.,
1643      e     0.,
1644      s     1.,0.,-1.,0.,1.,-1.,0.,1.,2.,0.,
1645      1     1.,0.,1.,0.,-1.,0.,1.,0.,0.,0.,
1646      2     -1.,0.,1.,0.,-1.,0.,1.,0.,0.,1.,
1647      3     0.,0.,0.,0.,0.,0.,0.,0.,0.,-1.,
1648      4     0.,0.,0.,0.,-1./
1649
1650 clin-2/19/03 initialize n and nsav for resonance decay at each timestep
1651 c     in order to prevent integer overflow:
1652       call inidcy
1653
1654 c OFF skip ART collisions to reproduce HJ:      
1655 cc       if(nt.ne.ntmax) return
1656
1657 clin-11/07/00 rrkk is assumed to be 0.6mb(default) for mm->KKbar 
1658 c     with m=rho or omega, estimated from Ko's paper:
1659 c      rrkk=0.6
1660 c prkk: cross section of pi (rho or omega) -> K* Kbar (AND) K*bar K:
1661 c      prkk=0.3
1662 c     cross section in mb for (rho or omega) K* -> pi K:
1663 c      srhoks=5.
1664 clin-11/07/00-end
1665 c      ESBIN=0.04
1666       RESONA=5.
1667 *-----------------------------------------------------------------------
1668 *     INITIALIZATION OF COUNTING VARIABLES
1669       NODELT=0
1670       SUMSRT =0.
1671       LCOLL  = 0
1672       LBLOC  = 0
1673       LCNNE  = 0
1674       LDD  = 0
1675       LPP  = 0
1676       lpd  = 0
1677       lpdr=0
1678       lrho = 0
1679       lrhor=0
1680       lomega=0
1681       lomgar=0
1682       LPN  = 0
1683       LKN  = 0
1684       LNNK = 0
1685       LDDK = 0
1686       LNDK = 0
1687       lppk =0
1688       LCNND  = 0
1689       LCNDN  = 0
1690       LDIRT  = 0
1691       LDECAY = 0
1692       LRES   = 0
1693       Ldou   = 0
1694       LDDRHO = 0
1695       LNNRHO = 0
1696       LNNOM  = 0
1697       MSUM   = 0
1698       MASSRN(0)=0
1699 * COM: MSUM IS USED TO COUNT THE TOTAL NO. OF PARTICLES 
1700 *      IN PREVIOUS IRUN-1 RUNS
1701 * KAON COUNTERS
1702       DO 1002 IL=1,5
1703          TKAON(IL)=0
1704          DO 1001 IS=1,2000
1705             EKAON(IL,IS)=0
1706  1001    CONTINUE
1707  1002 CONTINUE
1708 c sp 12/19/00
1709       DO 1004 i =1,NUM
1710          DO 1003 j =1,MAXSTR
1711             PROPI(j,i) = 1.
1712  1003    CONTINUE
1713  1004 CONTINUE
1714       
1715       do 1102 i=1,maxstr
1716          fttemp(i)=0.
1717          do 1101 irun=1,maxr
1718             ftpisv(i,irun)=0.
1719  1101    continue
1720  1102 continue
1721
1722 c sp 12/19/00 end
1723       sp=0
1724 * antikaon counters
1725       akaon=0
1726       sk=0
1727 *-----------------------------------------------------------------------
1728 *     LOOP OVER ALL PARALLEL RUNS
1729 cbz11/17/98
1730 c      MASS=MASSPR+MASSTA
1731       MASS = 0
1732 cbz11/17/98end
1733       DO 1000 IRUN = 1,NUM
1734          NNN=0
1735          MSUM=MSUM+MASSR(IRUN-1)
1736 *     LOOP OVER ALL PSEUDOPARTICLES 1 IN THE SAME RUN
1737          J10=2
1738          IF(NT.EQ.NTMAX)J10=1
1739 c
1740 ctest off skips the check of energy conservation after each timestep:
1741 c         enetot=0.
1742 c         do ip=1,MASSR(IRUN)
1743 c            if(e(ip).ne.0.or.lb(ip).eq.10022) enetot=enetot
1744 c     1           +sqrt(p(1,ip)**2+p(2,ip)**2+p(3,ip)**2+e(ip)**2)
1745 c         enddo
1746 c         write(91,*) 'A:',nt,enetot,massr(irun),bimp 
1747
1748          DO 800 J1 = J10,MASSR(IRUN)
1749             I1  = J1 + MSUM
1750 * E(I)=0 are for pions having been absorbed or photons which do not enter here:
1751             IF(E(I1).EQ.0.)GO TO 800
1752
1753 c     To include anti-(Delta,N*1440 and N*1535):
1754 c          IF ((LB(I1) .LT. -13 .OR. LB(I1) .GT. 28)
1755 c     1         .and.iabs(LB(I1)) .ne. 30 ) GOTO 800
1756             IF (LB(I1) .LT. -45 .OR. LB(I1) .GT. 45) GOTO 800
1757             X1  = R(1,I1)
1758             Y1  = R(2,I1)
1759             Z1  = R(3,I1)
1760             PX1 = P(1,I1)
1761             PY1 = P(2,I1)
1762             PZ1 = P(3,I1)
1763             EM1 = E(I1)
1764             am1= em1
1765             E1  = SQRT( EM1**2 + PX1**2 + PY1**2 + PZ1**2 )
1766             ID1 = ID(I1)
1767             LB1 = LB(I1)
1768
1769 c     generate k0short and k0long from K+ and K- at the last timestep:
1770             if(nt.eq.ntmax.and.(lb1.eq.21.or.lb1.eq.23)) then
1771                pk0=RANART(NSEED)
1772                if(pk0.lt.0.25) then
1773                   LB(I1)=22
1774                elseif(pk0.lt.0.50) then
1775                   LB(I1)=24
1776                endif
1777                LB1=LB(I1)
1778             endif
1779             
1780 clin-8/07/02 these particles don't decay strongly, so skip decay routines:     
1781 c            IF( (lb1.ge.-2.and.lb1.le.5) .OR. lb1.eq.31 .OR.
1782 c     &           (iabs(lb1).ge.14.and.iabs(lb1).le.24) .OR.
1783 c     &           (iabs(lb1).ge.40.and.iabs(lb1).le.45) .or. 
1784 c     &           lb1.eq.31)GO TO 1 
1785 c     only decay K0short when iksdcy=1:
1786             if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1787      &           .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1788      &           .or.(iabs(lb1).ge.6.and.iabs(lb1).le.13)
1789      &           .or.(iksdcy.eq.1.and.lb1.eq.24)
1790      &           .or.iabs(lb1).eq.16) then
1791                continue
1792             else
1793                goto 1
1794             endif
1795 * IF I1 IS A RESONANCE, CHECK WHETHER IT DECAYS DURING THIS TIME STEP
1796          IF(lb1.ge.25.and.lb1.le.27) then
1797              wid=0.151
1798          ELSEIF(lb1.eq.28) then
1799              wid=0.00841
1800          ELSEIF(lb1.eq.29) then
1801              wid=0.00443
1802           ELSEIF(iabs(LB1).eq.30) then
1803              WID=0.051
1804          ELSEIF(lb1.eq.0) then
1805              wid=1.18e-6
1806 c     to give K0short ct0=2.676cm:
1807          ELSEIF(iksdcy.eq.1.and.lb1.eq.24) then
1808              wid=7.36e-15
1809 clin-4/29/03 add Sigma0 decay to Lambda, ct0=2.22E-11m:
1810          ELSEIF(iabs(lb1).eq.16) then
1811              wid=8.87e-6
1812 csp-07/25/01 test a1 resonance:
1813 cc          ELSEIF(LB1.EQ.32) then
1814 cc             WID=0.40
1815           ELSEIF(LB1.EQ.32) then
1816              call WIDA1(EM1,rhomp,WID,iseed)
1817           ELSEIF(iabs(LB1).ge.6.and.iabs(LB1).le.9) then
1818              WID=WIDTH(EM1)
1819           ELSEIF((iabs(LB1).EQ.10).OR.(iabs(LB1).EQ.11)) then
1820              WID=W1440(EM1)
1821           ELSEIF((iabs(LB1).EQ.12).OR.(iabs(LB1).EQ.13)) then
1822              WID=W1535(EM1)
1823           ENDIF
1824
1825 * if it is the last time step, FORCE all resonance to strong-decay
1826 * and go out of the loop
1827           if(nt.eq.ntmax)then
1828              pdecay=1.1
1829 clin-5b/2008 forbid phi decay at the end of hadronic cascade:
1830              if(iphidcy.eq.0.and.iabs(LB1).eq.29) pdecay=0.
1831           else
1832              T0=0.19733/WID
1833              GFACTR=E1/EM1
1834              T0=T0*GFACTR
1835              IF(T0.GT.0.)THEN
1836                 PDECAY=1.-EXP(-DT/T0)
1837              ELSE
1838                 PDECAY=0.
1839              ENDIF
1840           endif
1841           XDECAY=RANART(NSEED)
1842
1843 cc dilepton production from rho0, omega, phi decay 
1844 cc        if(lb1.eq.26 .or. lb1.eq.28 .or. lb1.eq.29)
1845 cc     &   call dec_ceres(nt,ntmax,irun,i1)
1846 cc
1847           IF(XDECAY.LT.PDECAY) THEN
1848 clin-10/25/02 get rid of argument usage mismatch in rhocay():
1849              idecay=irun
1850              tfnl=nt*dt
1851 clin-10/28/03 keep formation time of hadrons unformed at nt=ntmax-1:
1852              if(nt.eq.ntmax.and.ftsv(i1).gt.((ntmax-1)*dt)) 
1853      1            tfnl=ftsv(i1)
1854              xfnl=x1
1855              yfnl=y1
1856              zfnl=z1
1857 * use PYTHIA to perform decays of eta,rho,omega,phi,K*,(K0s) and Delta:
1858              if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1859      &           .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1860      &           .or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1861      &           .or.(iksdcy.eq.1.and.lb1.eq.24)
1862      &           .or.iabs(lb1).eq.16) then
1863 c     previous rho decay performed in rhodecay():
1864 c                nnn=nnn+1
1865 c                call rhodecay(idecay,i1,nnn,iseed)
1866 c
1867 ctest off record decays of phi,K*,Lambda(1520) resonances:
1868 c                if(lb1.eq.29.or.iabs(lb1).eq.30) 
1869 c     1               write(18,112) 'decay',lb1,px1,py1,pz1,am1,nt
1870                 call resdec(i1,nt,nnn,wid,idecay)
1871                 p(1,i1)=px1n
1872                 p(2,i1)=py1n
1873                 p(3,i1)=pz1n
1874 clin-5/2008:
1875                 dpertp(i1)=dp1n
1876 c     add decay time to freezeout positions & time at the last timestep:
1877                 if(nt.eq.ntmax) then
1878                    R(1,i1)=xfnl
1879                    R(2,i1)=yfnl
1880                    R(3,i1)=zfnl
1881                    tfdcy(i1)=tfnl
1882                 endif
1883 c
1884 * decay number for baryon resonance or L/S decay
1885                 if(iabs(lb1).ge.6.and.iabs(lb1).le.9) then
1886                    LDECAY=LDECAY+1
1887                 endif
1888
1889 * for a1 decay 
1890 c             elseif(lb1.eq.32)then
1891 c                NNN=NNN+1
1892 c                call a1decay(idecay,i1,nnn,iseed,rhomp)
1893
1894 * FOR N*(1440)
1895              elseif(iabs(LB1).EQ.10.OR.iabs(LB1).EQ.11) THEN
1896                 NNN=NNN+1
1897                 LDECAY=LDECAY+1
1898                 PNSTAR=1.
1899                 IF(E(I1).GT.1.22)PNSTAR=0.6
1900                 IF(RANART(NSEED).LE.PNSTAR)THEN
1901 * (1) DECAY TO SINGLE PION+NUCLEON
1902                    CALL DECAYA(idecay,I1,NNN,ISEED,wid,nt)
1903                 ELSE
1904 * (2) DECAY TO TWO PIONS + NUCLEON
1905                    CALL DECAY2(idecay,I1,NNN,ISEED,wid,nt)
1906                    NNN=NNN+1
1907                 ENDIF
1908 c for N*(1535) decay
1909              elseif(iabs(LB1).eq.12.or.iabs(LB1).eq.13) then
1910                 NNN=NNN+1
1911                 CALL DECAYA(idecay,I1,NNN,ISEED,wid,nt)
1912                 LDECAY=LDECAY+1
1913              endif
1914 c
1915 *COM: AT HIGH ENERGIES WE USE VERY SHORT TIME STEPS,
1916 *     IN ORDER TO TAKE INTO ACCOUNT THE FINITE FORMATIOM TIME, WE
1917 *     DO NOT ALLOW PARTICLES FROM THE DECAY OF RESONANCE TO INTERACT 
1918 *     WITH OTHERS IN THE SAME TIME STEP. CHANGE 9000 TO REVERSE THIS 
1919 *     ASSUMPTION. EFFECTS OF THIS ASSUMPTION CAN BE STUDIED BY CHANGING 
1920 *     THE STATEMENT OF 9000. See notebook for discussions on effects of
1921 *     changing statement 9000.
1922 c
1923 c     kaons from K* decay are converted to k0short (and k0long), 
1924 c     phi decay may produce rho, K0S or eta, N*(1535) decay may produce eta,
1925 c     and these decay daughters need to decay again if at the last timestep:
1926 c     (note: these daughters have been assigned to lb(i1) only, not to lpion)
1927 c             if(nt.eq.ntmax.and.(lb1.eq.29.or.iabs(lb1).eq.30
1928 c     1            .iabs(lb1).eq.12.or.iabs(lb1).eq.13)) then
1929              if(nt.eq.ntmax) then
1930                 if(lb(i1).eq.25.or.lb(i1).eq.26.or.lb(i1).eq.27) then
1931                    wid=0.151
1932                 elseif(lb(i1).eq.0) then
1933                    wid=1.18e-6
1934                 elseif(lb(i1).eq.24.and.iksdcy.eq.1) then
1935                    wid=7.36e-17
1936                 else
1937                    goto 9000
1938                 endif
1939                 LB1=LB(I1)
1940                 PX1=P(1,I1)
1941                 PY1=P(2,I1)
1942                 PZ1=P(3,I1)
1943                 EM1=E(I1)
1944                 E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1945                 call resdec(i1,nt,nnn,wid,idecay)
1946                 p(1,i1)=px1n
1947                 p(2,i1)=py1n
1948                 p(3,i1)=pz1n
1949                 R(1,i1)=xfnl
1950                 R(2,i1)=yfnl
1951                 R(3,i1)=zfnl
1952                 tfdcy(i1)=tfnl
1953 clin-5/2008:
1954                 dpertp(i1)=dp1n
1955              endif
1956
1957 * negelecting the Pauli blocking at high energies
1958  9000        go to 800
1959           ENDIF
1960 * LOOP OVER ALL PSEUDOPARTICLES 2 IN THE SAME RUN
1961 * SAVE ALL THE COORDINATES FOR POSSIBLE CHANGE IN THE FOLLOWING COLLISION
1962  1        if(nt.eq.ntmax)go to 800
1963           X1 = R(1,I1)
1964           Y1 = R(2,I1)
1965           Z1 = R(3,I1)
1966 c
1967            DO 600 J2 = 1,J1-1
1968             I2  = J2 + MSUM
1969 * IF I2 IS A MESON BEING ABSORBED, THEN GO OUT OF THE LOOP
1970             IF(E(I2).EQ.0.) GO TO 600
1971 clin-5/2008 in case the first particle is already destroyed:
1972             IF(E(I1).EQ.0.) GO TO 800
1973             IF (LB(I2) .LT. -45 .OR. LB(I2) .GT. 45) GOTO 600
1974 clin-7/26/03 improve speed
1975             X2=R(1,I2)
1976             Y2=R(2,I2)
1977             Z2=R(3,I2)
1978             dr0max=5.
1979 clin-9/2008 deuteron+nucleon elastic cross sections could reach ~2810mb:
1980             ilb1=iabs(LB(I1))
1981             ilb2=iabs(LB(I2))
1982             IF(ilb1.EQ.42.or.ilb2.EQ.42) THEN
1983                if((ILB1.GE.1.AND.ILB1.LE.2)
1984      1              .or.(ILB1.GE.6.AND.ILB1.LE.13)
1985      2              .or.(ILB2.GE.1.AND.ILB2.LE.2)
1986      3              .or.(ILB2.GE.6.AND.ILB2.LE.13)) then
1987                   if((lb(i1)*lb(i2)).gt.0) dr0max=10.
1988                endif
1989             ENDIF
1990 c
1991             if(((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2).GT.dr0max**2)
1992      1           GO TO 600
1993             IF (ID(I1)*ID(I2).EQ.IAVOID) GOTO 400
1994             ID1=ID(I1)
1995             ID2 = ID(I2)
1996 c
1997             ix1= nint(x1/dx)
1998             iy1= nint(y1/dy)
1999             iz1= nint(z1/dz)
2000             PX1=P(1,I1)
2001             PY1=P(2,I1)
2002             PZ1=P(3,I1)
2003             EM1=E(I1)
2004             AM1=EM1
2005             LB1=LB(I1)
2006             E1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
2007             IPX1=NINT(PX1/DPX)
2008             IPY1=NINT(PY1/DPY)
2009             IPZ1=NINT(PZ1/DPZ)         
2010             LB2 = LB(I2)
2011             PX2 = P(1,I2)
2012             PY2 = P(2,I2)
2013             PZ2 = P(3,I2)
2014             EM2=E(I2)
2015             AM2=EM2
2016             lb1i=lb(i1)
2017             lb2i=lb(i2)
2018             px1i=P(1,I1)
2019             py1i=P(2,I1)
2020             pz1i=P(3,I1)
2021             em1i=E(I1)
2022             px2i=P(1,I2)
2023             py2i=P(2,I2)
2024             pz2i=P(3,I2)
2025             em2i=E(I2)
2026 clin-2/26/03 ctest off check energy conservation after each binary search:
2027             eini=SQRT(E(I1)**2+P(1,I1)**2+P(2,I1)**2+P(3,I1)**2)
2028      1           +SQRT(E(I2)**2+P(1,I2)**2+P(2,I2)**2+P(3,I2)**2)
2029             pxini=P(1,I1)+P(1,I2)
2030             pyini=P(2,I1)+P(2,I2)
2031             pzini=P(3,I1)+P(3,I2)
2032             nnnini=nnn
2033 c
2034 clin-4/30/03 initialize value:
2035             iblock=0
2036 c
2037 * TO SAVE COMPUTING TIME we do the following
2038 * (1) make a ROUGH estimate to see whether particle i2 will collide with
2039 * particle I1, and (2) skip the particle pairs for which collisions are 
2040 * not modeled in the code.
2041 * FOR MESON-BARYON AND MESON-MESON COLLISIONS, we use a maximum 
2042 * interaction distance DELTR0=2.6
2043 * for ppbar production from meson (pi rho omega) interactions:
2044 c
2045             DELTR0=3.
2046         if( (iabs(lb1).ge.14.and.iabs(lb1).le.17) .or.
2047      &      (iabs(lb1).ge.30.and.iabs(lb1).le.45) ) DELTR0=5.0
2048         if( (iabs(lb2).ge.14.and.iabs(lb2).le.17) .or.
2049      &      (iabs(lb2).ge.30.and.iabs(lb2).le.45) ) DELTR0=5.0
2050
2051             if(lb1.eq.28.and.lb2.eq.28) DELTR0=4.84
2052 clin-10/08/00 to include pi pi -> rho rho:
2053             if((lb1.ge.3.and.lb1.le.5).and.(lb2.ge.3.and.lb2.le.5)) then
2054                E2=SQRT(EM2**2+PX2**2+PY2**2+PZ2**2)
2055          spipi=(e1+e2)**2-(px1+px2)**2-(py1+py2)**2-(pz1+pz2)**2
2056                if(spipi.ge.(4*0.77**2)) DELTR0=3.5
2057             endif
2058
2059 c khyperon
2060         IF (LB1.EQ.23 .AND. (LB2.GE.14.AND.LB2.LE.17)) GOTO 3699
2061         IF (LB2.EQ.23 .AND. (LB1.GE.14.AND.LB1.LE.17)) GOTO 3699
2062
2063 * K(K*) + Kbar(K*bar) scattering including 
2064 *     K(K*) + Kbar(K*bar) --> phi + pi(rho,omega) and pi pi(rho,omega)
2065        if(lb1.eq.21.and.lb2.eq.23)go to 3699
2066        if(lb2.eq.21.and.lb1.eq.23)go to 3699
2067        if(lb1.eq.30.and.lb2.eq.21)go to 3699
2068        if(lb2.eq.30.and.lb1.eq.21)go to 3699
2069        if(lb1.eq.-30.and.lb2.eq.23)go to 3699
2070        if(lb2.eq.-30.and.lb1.eq.23)go to 3699
2071        if(lb1.eq.-30.and.lb2.eq.30)go to 3699
2072        if(lb2.eq.-30.and.lb1.eq.30)go to 3699
2073 c
2074 clin-12/15/00
2075 c     kaon+rho(omega,eta) collisions:
2076       if(lb1.eq.21.or.lb1.eq.23) then
2077          if(lb2.eq.0.or.(lb2.ge.25.and.lb2.le.28)) then
2078             go to 3699
2079          endif
2080       elseif(lb2.eq.21.or.lb2.eq.23) then
2081          if(lb1.eq.0.or.(lb1.ge.25.and.lb1.le.28)) then
2082             goto 3699
2083          endif
2084       endif
2085
2086 clin-8/14/02 K* (pi, rho, omega, eta) collisions:
2087       if(iabs(lb1).eq.30 .and.
2088      1     (lb2.eq.0.or.(lb2.ge.25.and.lb2.le.28)
2089      2     .or.(lb2.ge.3.and.lb2.le.5))) then
2090          go to 3699
2091       elseif(iabs(lb2).eq.30 .and.
2092      1        (lb1.eq.0.or.(lb1.ge.25.and.lb1.le.28)
2093      2        .or.(lb1.ge.3.and.lb1.le.5))) then
2094          goto 3699
2095 clin-8/14/02-end
2096 c K*/K*-bar + baryon/antibaryon collisions:
2097         elseif( iabs(lb1).eq.30 .and.
2098      1     (iabs(lb2).eq.1.or.iabs(lb2).eq.2.or.
2099      2     (iabs(lb2).ge.6.and.iabs(lb2).le.13)) )then
2100               go to 3699
2101            endif
2102          if( iabs(lb2).eq.30 .and.
2103      1         (iabs(lb1).eq.1.or.iabs(lb1).eq.2.or.
2104      2         (iabs(lb1).ge.6.and.iabs(lb1).le.13)) )then
2105                 go to 3699
2106         endif                                                              
2107 * K^+ baryons and antibaryons:
2108 c** K+ + B-bar  --> La(Si)-bar + pi
2109 * K^- and antibaryons, note K^- and baryons are included in newka():
2110 * note that we fail to satisfy charge conjugation for these cross sections:
2111         if((lb1.eq.23.or.lb1.eq.21).and.
2112      1       (iabs(lb2).eq.1.or.iabs(lb2).eq.2.or.
2113      2       (iabs(lb2).ge.6.and.iabs(lb2).le.13))) then
2114            go to 3699
2115         elseif((lb2.eq.23.or.lb2.eq.21).and.
2116      1       (iabs(lb1).eq.1.or.iabs(lb1).eq.2.or.
2117      2       (iabs(lb1).ge.6.and.iabs(lb1).le.13))) then
2118            go to 3699
2119         endif
2120 *
2121 * For anti-nucleons annihilations:
2122 * Assumptions: 
2123 * (1) for collisions involving a p_bar or n_bar,
2124 * we allow only collisions between a p_bar and a baryon or a baryon 
2125 * resonance (as well as a n_bar and a baryon or a baryon resonance),
2126 * we skip all other reactions involving a p_bar or n_bar, 
2127 * such as collisions between p_bar (n_bar) and mesons, 
2128 * and collisions between two p_bar's (n_bar's). 
2129 * (2) we introduce a new parameter rppmax: the maximum interaction 
2130 * distance to make the quick collision check,rppmax=3.57 fm 
2131 * corresponding to a cutoff of annihilation xsection= 400mb which is
2132 * also used consistently in the actual annihilation xsection to be 
2133 * used in the following as given in the subroutine xppbar(srt)
2134         rppmax=3.57   
2135 * anti-baryon on baryons
2136         if((lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))
2137      1 .and.(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13))) then
2138             DELTR0 = RPPMAX
2139             GOTO 2699
2140        else if((lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6))
2141      1 .and.(lb1.eq.1.or.lb1.eq.2.or.(lb1.ge.6.and.lb1.le.13))) then
2142             DELTR0 = RPPMAX
2143             GOTO 2699
2144          END IF
2145
2146 c*  ((anti) lambda, cascade, omega  should not be rejected)
2147         if( (iabs(lb1).ge.14.and.iabs(lb1).le.17) .or.
2148      &      (iabs(lb2).ge.14.and.iabs(lb2).le.17) )go to 3699
2149 c
2150 clin-9/2008 maximum sigma~2810mb for deuteron+nucleon elastic collisions:
2151          IF (iabs(LB1).EQ.42.or.iabs(LB2).EQ.42) THEN
2152             ilb1=iabs(LB1)
2153             ilb2=iabs(LB2)
2154             if((ILB1.GE.1.AND.ILB1.LE.2)
2155      1           .or.(ILB1.GE.6.AND.ILB1.LE.13)
2156      2           .or.(ILB2.GE.1.AND.ILB2.LE.2)
2157      3           .or.(ILB2.GE.6.AND.ILB2.LE.13)) then
2158                if((lb1*lb2).gt.0) deltr0=9.5
2159             endif
2160          ENDIF
2161 c
2162         if( (iabs(lb1).ge.40.and.iabs(lb1).le.45) .or. 
2163      &      (iabs(lb2).ge.40.and.iabs(lb2).le.45) )go to 3699
2164 c
2165 c* phi channel --> elastic + inelastic scatt.  
2166          IF( (lb1.eq.29 .and.((lb2.ge.1.and.lb2.le.13).or.  
2167      &       (lb2.ge.21.and.lb2.le.28).or.iabs(lb2).eq.30)) .OR.
2168      &     (lb2.eq.29 .and.((lb1.ge.1.and.lb1.le.13).or.
2169      &       (lb1.ge.21.and.lb1.le.28).or.iabs(lb1).eq.30)) )THEN
2170              DELTR0=3.0
2171              go to 3699
2172         endif
2173 c
2174 c  La/Si, Cas, Om (bar)-meson elastic colln
2175 * pion vs. La & Ca (bar) coll. are treated in resp. subroutines
2176
2177 * SKIP all other K* RESCATTERINGS
2178         If(iabs(lb1).eq.30.or.iabs(lb2).eq.30) go to 400
2179 * SKIP KAON(+) RESCATTERINGS WITH particles other than pions and baryons 
2180          If(lb1.eq.23.and.(lb2.lt.1.or.lb2.gt.17))go to 400
2181          If(lb2.eq.23.and.(lb1.lt.1.or.lb1.gt.17))go to 400
2182 c
2183 c anti-baryon proccess: B-bar+M, N-bar+R-bar, N-bar+N-bar, R-bar+R-bar
2184 c  R = (D,N*)
2185          if( ((lb1.le.-1.and.lb1.ge.-13)
2186      &        .and.(lb2.eq.0.or.(lb2.ge.3.and.lb2.le.5)
2187      &            .or.(lb2.ge.25.and.lb2.le.28))) 
2188      &      .OR.((lb2.le.-1.and.lb2.ge.-13)
2189      &         .and.(lb1.eq.0.or.(lb1.ge.3.and.lb1.le.5)
2190      &              .or.(lb1.ge.25.and.lb1.le.28))) ) then
2191          elseIF( ((LB1.eq.-1.or.lb1.eq.-2).
2192      &             and.(LB2.LT.-5.and.lb2.ge.-13))
2193      &      .OR. ((LB2.eq.-1.or.lb2.eq.-2).
2194      &             and.(LB1.LT.-5.and.lb1.ge.-13)) )then
2195          elseIF((LB1.eq.-1.or.lb1.eq.-2)
2196      &     .AND.(LB2.eq.-1.or.lb2.eq.-2))then
2197          elseIF((LB1.LT.-5.and.lb1.ge.-13).AND.
2198      &          (LB2.LT.-5.and.lb2.ge.-13)) then
2199 c        elseif((lb1.lt.0).or.(lb2.lt.0)) then
2200 c         go to 400
2201        endif               
2202
2203  2699    CONTINUE
2204 * for baryon-baryon collisions
2205          IF (LB1 .EQ. 1 .OR. LB1 .EQ. 2 .OR. (LB1 .GE. 6 .AND.
2206      &        LB1 .LE. 17)) THEN
2207             IF (LB2 .EQ. 1 .OR. LB2 .EQ. 2 .OR. (LB2 .GE. 6 .AND.
2208      &           LB2 .LE. 17)) THEN
2209                DELTR0 = 2.
2210             END IF
2211          END IF
2212 c
2213  3699   RSQARE = (X1-X2)**2 + (Y1-Y2)**2 + (Z1-Z2)**2
2214         IF (RSQARE .GT. DELTR0**2) GO TO 400
2215 *NOW PARTICLES ARE CLOSE ENOUGH TO EACH OTHER !
2216 * KEEP ALL COORDINATES FOR POSSIBLE PHASE SPACE CHANGE
2217             ix2 = nint(x2/dx)
2218             iy2 = nint(y2/dy)
2219             iz2 = nint(z2/dz)
2220             ipx2 = nint(px2/dpx)
2221             ipy2 = nint(py2/dpy)
2222             ipz2 = nint(pz2/dpz)
2223 * FIND MOMENTA OF PARTICLES IN THE CMS OF THE TWO COLLIDING PARTICLES
2224 * AND THE CMS ENERGY SRT
2225           CALL CMS(I1,I2,PCX,PCY,PCZ,SRT)
2226 clin-7/26/03 improve speed
2227           drmax=dr0max
2228           call distc0(drmax,deltr0,DT,
2229      1         Ifirst,PCX,PCY,PCZ,
2230      2         x1,y1,z1,px1,py1,pz1,em1,x2,y2,z2,px2,py2,pz2,em2)
2231           if(Ifirst.eq.-1) goto 400
2232
2233          ISS=NINT(SRT/ESBIN)
2234 clin-4/2008 use last bin if ISS is out of EKAON's upper bound of 2000:
2235          if(ISS.gt.2000) ISS=2000
2236 *Sort collisions
2237 c
2238 clin-8/2008 Deuteron+Meson->B+B; 
2239 c     meson=(pi,rho,omega,eta), B=(n,p,Delta,N*1440,N*1535):
2240          IF (iabs(LB1).EQ.42.or.iabs(LB2).EQ.42) THEN
2241             ilb1=iabs(LB1)
2242             ilb2=iabs(LB2)
2243             if(LB1.eq.0.or.(LB1.GE.3.AND.LB1.LE.5)
2244      1           .or.(LB1.GE.25.AND.LB1.LE.28)
2245      2           .or.
2246      3           LB2.eq.0.or.(LB2.GE.3.AND.LB2.LE.5)
2247      4           .or.(LB2.GE.25.AND.LB2.LE.28)) then
2248                GOTO 505
2249 clin-9/2008 Deuteron+Baryon or antiDeuteron+antiBaryon elastic collisions:
2250             elseif(((ILB1.GE.1.AND.ILB1.LE.2)
2251      1              .or.(ILB1.GE.6.AND.ILB1.LE.13)
2252      2              .or.(ILB2.GE.1.AND.ILB2.LE.2)
2253      3              .or.(ILB2.GE.6.AND.ILB2.LE.13))
2254      4              .and.(lb1*lb2).gt.0) then
2255                GOTO 506
2256             else
2257                GOTO 400
2258             endif
2259          ENDIF
2260 c
2261 * K+ + (N,N*,D)-bar --> L/S-bar + pi
2262           if( ((lb1.eq.23.or.lb1.eq.30).and.
2263      &         (lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6))) 
2264      &         .OR.((lb2.eq.23.or.lb2.eq.30).and.
2265      &         (lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))) )
2266      &         then
2267              bmass=0.938
2268              if(srt.le.(bmass+aka)) then
2269                 pkaon=0.
2270              else
2271                 pkaon=sqrt(((srt**2-(aka**2+bmass**2))
2272      1               /2./bmass)**2-aka**2)
2273              endif
2274 clin-10/31/02 cross sections are isospin-averaged, same as those in newka
2275 c     for K- + (N,N*,D) --> L/S + pi:
2276              sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
2277              SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
2278              SIG = sigela + SIGSGM + AKPLAM(PKAON)
2279              if(sig.gt.1.e-7) then
2280 c     ! K+ + N-bar reactions
2281                 icase=3
2282                 brel=sigela/sig
2283                 brsgm=sigsgm/sig
2284                 brsig = sig
2285                 nchrg = 1
2286                 go to 3555
2287              endif
2288              go to 400
2289           endif
2290 c
2291 c
2292 c  meson + hyperon-bar -> K+ + N-bar
2293           if(((lb1.ge.-17.and.lb1.le.-14).and.(lb2.ge.3.and.lb2.le.5)) 
2294      &         .OR.((lb2.ge.-17.and.lb2.le.-14)
2295      &         .and.(lb1.ge.3.and.lb1.le.5)))then
2296              nchrg=-100
2297  
2298 C*       first classify the reactions due to total charge.
2299              if((lb1.eq.-15.and.(lb2.eq.5.or.lb2.eq.27)).OR.
2300      &            (lb2.eq.-15.and.(lb1.eq.5.or.lb1.eq.27))) then
2301                 nchrg=-2
2302 c     ! D-(bar)
2303                 bmass=1.232
2304                 go to 110
2305              endif
2306              if( (lb1.eq.-15.and.(lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
2307      &            lb2.eq.28)).OR.(lb2.eq.-15.and.(lb1.eq.0.or.
2308      &            lb1.eq.4.or.lb1.eq.26.or.lb1.eq.28)).OR.
2309      &   ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.5.or.lb2.eq.27)).OR.
2310      &   ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.5.or.lb1.eq.27)) )then
2311                 nchrg=-1
2312 c     ! n-bar
2313                 bmass=0.938
2314                 go to 110
2315              endif
2316              if(  (lb1.eq.-15.and.(lb2.eq.3.or.lb2.eq.25)).OR.
2317      &            (lb2.eq.-15.and.(lb1.eq.3.or.lb1.eq.25)).OR.
2318      &            (lb1.eq.-17.and.(lb2.eq.5.or.lb2.eq.27)).OR.
2319      &            (lb2.eq.-17.and.(lb1.eq.5.or.lb1.eq.27)).OR.
2320      &            ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.0.or.lb2.eq.4
2321      &            .or.lb2.eq.26.or.lb2.eq.28)).OR.
2322      &            ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.0.or.lb1.eq.4
2323      &            .or.lb1.eq.26.or.lb1.eq.28)) )then
2324                nchrg=0
2325 c     ! p-bar
2326                 bmass=0.938
2327                 go to 110
2328              endif
2329              if( (lb1.eq.-17.and.(lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
2330      &            lb2.eq.28)).OR.(lb2.eq.-17.and.(lb1.eq.0.or.
2331      &            lb1.eq.4.or.lb1.eq.26.or.lb1.eq.28)).OR.
2332      &  ((lb1.eq.-14.or.lb1.eq.-16).and.(lb2.eq.3.or.lb2.eq.25)).OR.
2333      &  ((lb2.eq.-14.or.lb2.eq.-16).and.(lb1.eq.3.or.lb1.eq.25)))then
2334                nchrg=1
2335 c     ! D++(bar)
2336                 bmass=1.232
2337              endif
2338 c
2339 c 110     if(nchrg.ne.-100.and.srt.ge.(aka+bmass))then !! for elastic
2340  110         sig = 0.
2341 c !! for elastic
2342          if(nchrg.ne.-100.and.srt.ge.(aka+bmass))then
2343 cc110        if(nchrg.eq.-100.or.srt.lt.(aka+bmass)) go to 400
2344 c             ! PI + La(Si)-bar => K+ + N-bar reactions
2345             icase=4
2346 cc       pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
2347             pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
2348 c ! lambda-bar + Pi
2349             if(lb1.eq.-14.or.lb2.eq.-14) then
2350                if(nchrg.ge.0) sigma0=akPlam(pkaon)
2351                if(nchrg.lt.0) sigma0=akNlam(pkaon)
2352 c                ! sigma-bar + pi
2353             else
2354 c !K-p or K-D++
2355                if(nchrg.ge.0) sigma0=akPsgm(pkaon)
2356 c !K-n or K-D-
2357                if(nchrg.lt.0) sigma0=akNsgm(pkaon)
2358                SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
2359             endif
2360             sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
2361      &           (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
2362 c ! K0barD++, K-D-
2363             if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
2364 C*     the factor 2 comes from spin of delta, which is 3/2
2365 C*     detailed balance. copy from Page 423 of N.P. A614 1997
2366             IF (LB1 .EQ. -14 .OR. LB2 .EQ. -14) THEN
2367                SIG = 4.0 / 3.0 * SIG
2368             ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN
2369                SIG = 8.0 / 9.0 * SIG
2370             ELSE
2371                SIG = 4.0 / 9.0 * SIG
2372             END IF
2373 cc        brel=0.
2374 cc        brsgm=0.
2375 cc        brsig = sig
2376 cc          if(sig.lt.1.e-7) go to 400
2377 *-
2378          endif
2379 c                ! PI + La(Si)-bar => elastic included
2380          icase=4
2381          sigela = 10.
2382          sig = sig + sigela
2383          brel= sigela/sig
2384          brsgm=0.
2385          brsig = sig
2386 *-
2387          go to 3555
2388       endif
2389       
2390 ** MULTISTRANGE PARTICLE (Cas,Omega -bar) PRODUCTION - (NON)PERTURBATIVE
2391
2392 * K-/K*0bar + La/Si --> cascade + pi/eta
2393       if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.ge.14.and.lb2.le.17)).OR.
2394      &  ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.ge.14.and.lb1.le.17)) )then
2395           kp = 0
2396           go to 3455
2397         endif
2398 c K+/K*0 + La/Si(bar) --> cascade-bar + pi/eta
2399       if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.le.-14.and.lb2.ge.-17)).OR.
2400      &  ((lb2.eq.23.or.lb2.eq.30).and.(lb1.le.-14.and.lb1.ge.-17)) )then
2401           kp = 1
2402           go to 3455
2403         endif
2404 * K-/K*0bar + cascade --> omega + pi
2405        if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.eq.40.or.lb2.eq.41)).OR.
2406      & ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.eq.40.or.lb1.eq.41)) )then
2407           kp = 0
2408           go to 3455
2409         endif
2410 * K+/K*0 + cascade-bar --> omega-bar + pi
2411        if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.eq.-40.or.lb2.eq.-41)).OR.
2412      &  ((lb2.eq.23.or.lb2.eq.30).and.(lb1.eq.-40.or.lb1.eq.-41)) )then
2413           kp = 1
2414           go to 3455
2415         endif
2416 * Omega + Omega --> Di-Omega + photon(eta)
2417 cc        if( lb1.eq.45.and.lb2.eq.45 ) go to 3455
2418
2419 c annhilation of cascade(bar), omega(bar)
2420          kp = 3
2421 * K- + L/S <-- cascade(bar) + pi/eta
2422        if( (((lb1.ge.3.and.lb1.le.5).or.lb1.eq.0) 
2423      &       .and.(iabs(lb2).eq.40.or.iabs(lb2).eq.41))
2424      & .OR. (((lb2.ge.3.and.lb2.le.5).or.lb2.eq.0) 
2425      &       .and.(iabs(lb1).eq.40.or.iabs(lb1).eq.41)) )go to 3455
2426 * K- + cascade(bar) <-- omega(bar) + pi
2427 *         if(  (lb1.eq.0.and.iabs(lb2).eq.45)
2428 *    &       .OR. (lb2.eq.0.and.iabs(lb1).eq.45) )go to 3455
2429         if( ((lb1.ge.3.and.lb1.le.5).and.iabs(lb2).eq.45)
2430      &  .OR.((lb2.ge.3.and.lb2.le.5).and.iabs(lb1).eq.45) )go to 3455
2431 c
2432
2433 ***  MULTISTRANGE PARTICLE PRODUCTION  (END)
2434
2435 c* K+ + La(Si) --> Meson + B
2436         IF (LB1.EQ.23 .AND. (LB2.GE.14.AND.LB2.LE.17)) GOTO 5699
2437         IF (LB2.EQ.23 .AND. (LB1.GE.14.AND.LB1.LE.17)) GOTO 5699
2438 c* K- + La(Si)-bar --> Meson + B-bar
2439        IF (LB1.EQ.21 .AND. (LB2.GE.-17.AND.LB2.LE.-14)) GOTO 5699
2440        IF (LB2.EQ.21 .AND. (LB1.GE.-17.AND.LB1.LE.-14)) GOTO 5699
2441
2442 c La/Si-bar + B --> pi + K+
2443        IF( (((LB1.eq.1.or.LB1.eq.2).or.(LB1.ge.6.and.LB1.le.13))
2444      &       .AND.(LB2.GE.-17.AND.LB2.LE.-14)) .OR.
2445      &     (((LB2.eq.1.or.LB2.eq.2).or.(LB2.ge.6.and.LB2.le.13))
2446      &      .AND.(LB1.GE.-17.AND.LB1.LE.-14)) )go to 5999
2447 c La/Si + B-bar --> pi + K-
2448        IF( (((LB1.eq.-1.or.LB1.eq.-2).or.(LB1.le.-6.and.LB1.ge.-13))
2449      &       .AND.(LB2.GE.14.AND.LB2.LE.17)) .OR.
2450      &     (((LB2.eq.-1.or.LB2.eq.-2).or.(LB2.le.-6.and.LB2.ge.-13))
2451      &       .AND.(LB1.GE.14.AND.LB1.LE.17)) )go to 5999 
2452 *
2453 *
2454 * K(K*) + Kbar(K*bar) --> phi + pi(rho,omega), M + M (M=pi,rho,omega,eta)
2455        if(lb1.eq.21.and.lb2.eq.23) go to 8699
2456        if(lb2.eq.21.and.lb1.eq.23) go to 8699
2457        if(lb1.eq.30.and.lb2.eq.21) go to 8699
2458        if(lb2.eq.30.and.lb1.eq.21) go to 8699
2459        if(lb1.eq.-30.and.lb2.eq.23) go to 8699
2460        if(lb2.eq.-30.and.lb1.eq.23) go to 8699
2461        if(lb1.eq.-30.and.lb2.eq.30) go to 8699
2462        if(lb2.eq.-30.and.lb1.eq.30) go to 8699
2463 c* (K,K*)-bar + rho(omega) --> phi +(K,K*)-bar, piK and elastic
2464        IF( ((lb1.eq.23.or.lb1.eq.21.or.iabs(lb1).eq.30) .and.
2465      &      (lb2.ge.25.and.lb2.le.28)) .OR.
2466      &     ((lb2.eq.23.or.lb2.eq.21.or.iabs(lb2).eq.30) .and.
2467      &      (lb1.ge.25.and.lb1.le.28)) ) go to 8799
2468 c
2469 c* K*(-bar) + pi --> phi + (K,K*)-bar
2470        IF( (iabs(lb1).eq.30.and.(lb2.ge.3.and.lb2.le.5)) .OR.
2471      &     (iabs(lb2).eq.30.and.(lb1.ge.3.and.lb1.le.5)) )go to 8799
2472 *
2473 c
2474 c* phi + N --> pi+N(D),  rho+N(D),  K+ +La
2475 c* phi + D --> pi+N(D),  rho+N(D)
2476        IF( (lb1.eq.29 .and.(lb2.eq.1.or.lb2.eq.2.or.
2477      &       (lb2.ge.6.and.lb2.le.9))) .OR.
2478      &     (lb2.eq.29 .and.(lb1.eq.1.or.lb1.eq.2.or.
2479      &       (lb1.ge.6.and.lb1.le.9))) )go to 7222
2480 c
2481 c* phi + (pi,rho,ome,K,K*-bar) --> K+K, K+K*, K*+K*, (pi,rho,omega)+(K,K*-bar)
2482        IF( (lb1.eq.29 .and.((lb2.ge.3.and.lb2.le.5).or.
2483      &      (lb2.ge.21.and.lb2.le.28).or.iabs(lb2).eq.30)) .OR.
2484      &     (lb2.eq.29 .and.((lb1.ge.3.and.lb1.le.5).or.
2485      &      (lb1.ge.21.and.lb1.le.28).or.iabs(lb1).eq.30)) )THEN
2486              go to 7444
2487       endif
2488 *
2489 c
2490 * La/Si, Cas, Om (bar)-(rho,omega,phi) elastic colln
2491 * pion vs. La, Ca, Omega-(bar) elastic coll. treated in resp. subroutines
2492       if( ((iabs(lb1).ge.14.and.iabs(lb1).le.17).or.iabs(lb1).ge.40)
2493      &    .and.((lb2.ge.25.and.lb2.le.29).or.lb2.eq.0) )go to 888
2494       if( ((iabs(lb2).ge.14.and.iabs(lb2).le.17).or.iabs(lb2).ge.40)
2495      &    .and.((lb1.ge.25.and.lb1.le.29).or.lb1.eq.0) )go to 888
2496 c
2497 c K+/K* (N,R)  OR   K-/K*- (N,R)-bar  elastic scatt
2498         if( ((lb1.eq.23.or.lb1.eq.30).and.(lb2.eq.1.or.lb2.eq.2.or.
2499      &         (lb2.ge.6.and.lb2.le.13))) .OR.
2500      &      ((lb2.eq.23.or.lb2.eq.30).and.(lb1.eq.1.or.lb1.eq.2.or.
2501      &         (lb1.ge.6.and.lb1.le.13))) ) go to 888
2502         if( ((lb1.eq.21.or.lb1.eq.-30).and.(lb2.eq.-1.or.lb2.eq.-2.or.
2503      &       (lb2.ge.-13.and.lb2.le.-6))) .OR. 
2504      &      ((lb2.eq.21.or.lb2.eq.-30).and.(lb1.eq.-1.or.lb1.eq.-2.or.
2505      &       (lb1.ge.-13.and.lb1.le.-6))) ) go to 888
2506 c
2507 * L/S-baryon elastic collision 
2508        If( ((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.6.and.lb2.le.13))
2509      & .OR.((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.6.and.lb1.le.13)) )
2510      &   go to 7799
2511        If(((lb1.le.-14.and.lb1.ge.-17).and.(lb2.le.-6.and.lb2.ge.-13))
2512      &.OR.((lb2.le.-14.and.lb2.ge.-17).and.(lb1.le.-6.and.lb1.ge.-13)))
2513      &   go to 7799
2514 c
2515 c skip other collns with perturbative particles or hyperon-bar
2516        if( iabs(lb1).ge.40 .or. iabs(lb2).ge.40
2517      &    .or. (lb1.le.-14.and.lb1.ge.-17) 
2518      &    .or. (lb2.le.-14.and.lb2.ge.-17) )go to 400
2519 c
2520 c
2521 * anti-baryon on baryon resonaces 
2522         if((lb1.eq.-1.or.lb1.eq.-2.or.(lb1.ge.-13.and.lb1.le.-6))
2523      1 .and.(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13))) then
2524             GOTO 2799
2525        else if((lb2.eq.-1.or.lb2.eq.-2.or.(lb2.ge.-13.and.lb2.le.-6))
2526      1 .and.(lb1.eq.1.or.lb1.eq.2.or.(lb1.ge.6.and.lb1.le.13))) then
2527             GOTO 2799
2528          END IF
2529 c
2530 clin-10/25/02 get rid of argument usage mismatch in newka():
2531          inewka=irun
2532 c        call newka(icase,irun,iseed,dt,nt,
2533 clin-5/01/03 set iblock value in art1f.f, necessary for resonance studies:
2534 c        call newka(icase,inewka,iseed,dt,nt,
2535 c     &                  ictrl,i1,i2,srt,pcx,pcy,pcz)
2536         call newka(icase,inewka,iseed,dt,nt,
2537      &                  ictrl,i1,i2,srt,pcx,pcy,pcz,iblock)
2538
2539 clin-10/25/02-end
2540         IF (ICTRL .EQ. 1) GOTO 400
2541 c
2542 * SEPARATE NUCLEON+NUCLEON( BARYON RESONANCE+ BARYON RESONANCE ELASTIC
2543 * COLLISION), BARYON RESONANCE+NUCLEON AND BARYON-PION
2544 * COLLISIONS INTO THREE PARTS TO CHECK IF THEY ARE GOING TO SCATTER,
2545 * WE only allow L/S to COLLIDE elastically with a nucleon and meson
2546        if((iabs(lb1).ge.14.and.iabs(lb1).le.17).
2547      &  or.(iabs(lb2).ge.14.and.iabs(lb2).le.17))go to 400
2548 * IF PION+PION COLLISIONS GO TO 777
2549 * if pion+eta, eta+eta to create kaons go to 777 
2550        IF((lb1.ge.3.and.lb1.le.5).and.(lb2.ge.3.and.lb2.le.5))GO TO 777
2551        if(lb1.eq.0.and.(lb2.ge.3.and.lb2.le.5)) go to 777
2552        if(lb2.eq.0.and.(lb1.ge.3.and.lb1.le.5)) go to 777
2553        if(lb1.eq.0.and.lb2.eq.0)go to 777
2554 * we assume that rho and omega behave the same way as pions in 
2555 * kaon production
2556 * (1) rho(omega)+rho(omega)
2557        if( (lb1.ge.25.and.lb1.le.28).and.
2558      &     (lb2.ge.25.and.lb2.le.28) )goto 777
2559 * (2) rho(omega)+pion
2560       If((lb1.ge.25.and.lb1.le.28).and.(lb2.ge.3.and.lb2.le.5))go to 777
2561       If((lb2.ge.25.and.lb2.le.28).and.(lb1.ge.3.and.lb1.le.5))go to 777
2562 * (3) rho(omega)+eta
2563        if((lb1.ge.25.and.lb1.le.28).and.lb2.eq.0)go to 777
2564        if((lb2.ge.25.and.lb2.le.28).and.lb1.eq.0)go to 777
2565 c
2566 * if kaon+pion collisions go to 889
2567        if((lb1.eq.23.or.lb1.eq.21).and.(lb2.ge.3.and.lb2.le.5))go to 889
2568        if((lb2.eq.23.or.lb2.eq.21).and.(lb1.ge.3.and.lb1.le.5))go to 889
2569 c
2570 clin-2/06/03 skip all other (K K* Kbar K*bar) channels:
2571 * SKIP all other K and K* RESCATTERINGS
2572         If(iabs(lb1).eq.30.or.iabs(lb2).eq.30) go to 400
2573         If(lb1.eq.21.or.lb2.eq.21) go to 400
2574         If(lb1.eq.23.or.lb2.eq.23) go to 400
2575 c
2576 * IF PION+baryon COLLISION GO TO 3
2577            IF( (LB1.ge.3.and.LB1.le.5) .and. 
2578      &         (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2579      &          (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 3
2580            IF( (LB2.ge.3.and.LB2.le.5) .and. 
2581      &         (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2582      &          (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 3
2583 c
2584 * IF rho(omega)+NUCLEON (baryon resonance) COLLISION GO TO 33
2585            IF( (LB1.ge.25.and.LB1.le.28) .and. 
2586      &         (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2587      &          (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 33
2588            IF( (LB2.ge.25.and.LB2.le.28) .and. 
2589      &         (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2590      &          (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 33
2591 c
2592 * IF ETA+NUCLEON (baryon resonance) COLLISIONS GO TO 547
2593            IF( LB1.eq.0 .and. 
2594      &         (iabs(LB2).eq.1.or.iabs(LB2).eq.2.or.
2595      &          (iabs(LB2).ge.6.and.iabs(LB2).le.13)) )GO TO 547
2596            IF( LB2.eq.0 .and. 
2597      &         (iabs(LB1).eq.1.or.iabs(LB1).eq.2.or.
2598      &          (iabs(LB1).ge.6.and.iabs(LB1).le.13)) )GO TO 547
2599 c
2600 * IF NUCLEON+BARYON RESONANCE COLLISION GO TO 44
2601             IF((LB1.eq.1.or.lb1.eq.2).
2602      &        AND.(LB2.GT.5.and.lb2.le.13))GOTO 44
2603             IF((LB2.eq.1.or.lb2.eq.2).
2604      &        AND.(LB1.GT.5.and.lb1.le.13))GOTO 44
2605             IF((LB1.eq.-1.or.lb1.eq.-2).
2606      &        AND.(LB2.LT.-5.and.lb2.ge.-13))GOTO 44
2607             IF((LB2.eq.-1.or.lb2.eq.-2).
2608      &        AND.(LB1.LT.-5.and.lb1.ge.-13))GOTO 44
2609 c
2610 * IF NUCLEON+NUCLEON COLLISION GO TO 4
2611        IF((LB1.eq.1.or.lb1.eq.2).AND.(LB2.eq.1.or.lb2.eq.2))GOTO 4
2612        IF((LB1.eq.-1.or.lb1.eq.-2).AND.(LB2.eq.-1.or.lb2.eq.-2))GOTO 4
2613 c
2614 * IF BARYON RESONANCE+BARYON RESONANCE COLLISION GO TO 444
2615             IF((LB1.GT.5.and.lb1.le.13).AND.
2616      &         (LB2.GT.5.and.lb2.le.13)) GOTO 444
2617             IF((LB1.LT.-5.and.lb1.ge.-13).AND.
2618      &         (LB2.LT.-5.and.lb2.ge.-13)) GOTO 444
2619 c
2620 * if L/S+L/S or L/s+nucleon go to 400
2621 * otherwise, develop a model for their collisions
2622        if((lb1.lt.3).and.(lb2.ge.14.and.lb2.le.17))goto 400
2623        if((lb2.lt.3).and.(lb1.ge.14.and.lb1.le.17))goto 400
2624        if((lb1.ge.14.and.lb1.le.17).and.
2625      &  (lb2.ge.14.and.lb2.le.17))goto 400
2626 c
2627 * otherwise, go out of the loop
2628               go to 400
2629 *
2630 *
2631 547           IF(LB1*LB2.EQ.0)THEN
2632 * (1) FOR ETA+NUCLEON SYSTEM, we allow both elastic collision, 
2633 *     i.e. N*(1535) formation and kaon production
2634 *     the total kaon production cross section is
2635 *     ASSUMED to be THE SAME AS PION+NUCLEON COLLISIONS
2636 * (2) for eta+baryon resonance we only allow kaon production
2637            ece=(em1+em2+0.02)**2
2638            xkaon0=0.
2639            if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
2640            IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
2641 cbz3/7/99 neutralk
2642             XKAON0 = 2.0 * XKAON0
2643 cbz3/7/99 neutralk end
2644
2645 * Here we negelect eta+n inelastic collisions other than the 
2646 * kaon production, therefore the total inelastic cross section
2647 * xkaon equals to the xkaon0 (kaon production cross section)
2648            xkaon=xkaon0
2649 * note here the xkaon is in unit of fm**2
2650             XETA=XN1535(I1,I2,0)
2651         If((iabs(LB(I1)).ge.6.and.iabs(LB(I1)).le.13).or.
2652      &     (iabs(LB(I2)).ge.6.and.iabs(LB(I2)).le.13)) xeta=0.      
2653             IF((XETA+xkaon).LE.1.e-06)GO TO 400
2654             DSE=SQRT((XETA+XKAON)/PI)
2655            DELTRE=DSE+0.1
2656         px1cm=pcx
2657         py1cm=pcy
2658         pz1cm=pcz
2659 * CHECK IF N*(1535) resonance CAN BE FORMED
2660          CALL DISTCE(I1,I2,DELTRE,DSE,DT,ECE,SRT,IC,
2661      1   PCX,PCY,PCZ)
2662          IF(IC.EQ.-1) GO TO 400
2663          ekaon(4,iss)=ekaon(4,iss)+1
2664         IF(XKAON0/(XKAON+XETA).GT.RANART(NSEED))then
2665 * kaon production, USE CREN TO CALCULATE THE MOMENTUM OF L/S K+
2666         CALL CREN(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2667 * kaon production
2668        IF(IBLOCK.EQ.7) then
2669           LPN=LPN+1
2670        elseIF(IBLOCK.EQ.-7) then
2671        endif
2672 c
2673        em1=e(i1)
2674        em2=e(i2)
2675        GO TO 440
2676        endif
2677 * N*(1535) FORMATION
2678         resona=1.
2679          GO TO 98
2680          ENDIF
2681 *IF PION+NUCLEON (baryon resonance) COLLISION THEN
2682 3           CONTINUE
2683            px1cm=pcx
2684            py1cm=pcy
2685            pz1cm=pcz
2686 * the total kaon production cross section for pion+baryon (resonance) is
2687 * assumed to be the same as in pion+nucleon
2688            xkaon0=0.
2689            if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
2690            IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
2691             XKAON0 = 2.0 * XKAON0
2692 c
2693 c sp11/21/01  phi production: pi +N(D) -> phi + N(D)
2694          Xphi = 0.
2695        if( ( ((lb1.ge.1.and.lb1.le.2).or.
2696      &        (lb1.ge.6.and.lb1.le.9))
2697      &   .OR.((lb2.ge.1.and.lb2.le.2).or.
2698      &        (lb2.ge.6.and.lb2.le.9)) )
2699      &       .AND. srt.gt.1.958)
2700      &        call pibphi(srt,lb1,lb2,em1,em2,Xphi,xphin)
2701 c !! in fm^2 above
2702
2703 * if a pion collide with a baryon resonance, 
2704 * we only allow kaon production AND the reabsorption 
2705 * processes: Delta+pion-->N+pion, N*+pion-->N+pion
2706 * Later put in pion+baryon resonance elastic
2707 * cross through forming higher resonances implicitly.
2708 c          If(em1.gt.1.or.em2.gt.1.)go to 31
2709          If((iabs(LB(I1)).ge.6.and.iabs(LB(I1)).le.13).or.
2710      &      (iabs(LB(I2)).ge.6.and.iabs(LB(I2)).le.13)) go to 31
2711 * For pion+nucleon collisions: 
2712 * using the experimental pion+nucleon inelastic cross section, we assume it
2713 * is exhausted by the Delta+pion, Delta+rho and Delta+omega production 
2714 * and kaon production. In the following we first check whether 
2715 * inelastic pion+n collision can happen or not, then determine in 
2716 * crpn whether it is through pion production or through kaon production
2717 * note that the xkaon0 is the kaon production cross section
2718 * Note in particular that: 
2719 * xkaon in the following is the total pion+nucleon inelastic cross section
2720 * note here the xkaon is in unit of fm**2, xnpi is also in unit of fm**2
2721 * FOR PION+NUCLEON SYSTEM, THE MINIMUM S IS 1.2056 the minimum srt for 
2722 * elastic scattering, and it is 1.60 for pion production, 1.63 for LAMBDA+kaon 
2723 * production and 1.7 FOR SIGMA+KAON
2724 * (EC = PION MASS+NUCLEON MASS+20MEV)**2
2725             EC=(em1+em2+0.02)**2
2726            xkaon=0.
2727            if(srt.gt.1.23)xkaon=(pionpp(srt)+PIPP1(SRT))/2.
2728 * pion+nucleon elastic cross section is divided into two parts:
2729 * (1) forming D(1232)+N*(1440) +N*(1535)
2730 * (2) cross sections forming higher resonances are calculated as
2731 *     the difference between the total elastic and (1), this part is 
2732 *     treated as direct process since we do not explicitLY include
2733 *     higher resonances.
2734 * the following is the resonance formation cross sections.
2735 *1. PION(+)+PROTON-->DELTA++,PION(-)+NEUTRON-->DELTA(-)
2736            IF( (LB1*LB2.EQ.5.OR.((LB1*LB2.EQ.6).AND.
2737      &         (LB1.EQ.3.OR.LB2.EQ.3)))
2738      &    .OR. (LB1*LB2.EQ.-3.OR.((LB1*LB2.EQ.-10).AND.
2739      &         (LB1.EQ.5.OR.LB2.EQ.5))) )then    
2740               XMAX=190.
2741               xmaxn=0
2742               xmaxn1=0
2743               xdirct=dirct1(srt)
2744                go to 678
2745            endif
2746 *2. PION(-)+PROTON-->DELTA0,PION(+)+NEUTRON-->DELTA+ 
2747 *   or N*(+)(1440) or N*(+)(1535)
2748 * note the factor 2/3 is from the isospin consideration and
2749 * the factor 0.6 or 0.5 is the branching ratio for the resonance to decay
2750 * into pion+nucleon
2751             IF( (LB1*LB2.EQ.3.OR.((LB1*LB2.EQ.10).AND.
2752      &          (LB1.EQ.5.OR.LB2.EQ.5)))
2753      &     .OR. (LB1*LB2.EQ.-5.OR.((LB1*LB2.EQ.-6).AND.
2754      &          (LB1.EQ.3.OR.LB2.EQ.3))) )then      
2755               XMAX=27.
2756               xmaxn=2./3.*25.*0.6
2757                xmaxn1=2./3.*40.*0.5
2758               xdirct=dirct2(srt)
2759                go to 678
2760               endif
2761 *3. PION0+PROTON-->DELTA+,PION0+NEUTRON-->DELTA0, or N*(0)(1440) or N*(0)(1535)
2762             IF((LB1.EQ.4.OR.LB2.EQ.4).AND.
2763      &         (iabs(LB1*LB2).EQ.4.OR.iabs(LB1*LB2).EQ.8))then
2764               XMAX=50.
2765               xmaxn=1./3.*25*0.6
2766               xmaxn1=1/3.*40.*0.5
2767               xdirct=dirct3(srt)
2768                 go to 678
2769               endif
2770 678           xnpin1=0
2771            xnpin=0
2772             XNPID=XNPI(I1,I2,1,XMAX)
2773            if(xmaxn1.ne.0)xnpin1=XNPI(i1,i2,2,XMAXN1)
2774             if(xmaxn.ne.0)XNPIN=XNPI(I1,I2,0,XMAXN)
2775 * the following 
2776            xres=xnpid+xnpin+xnpin1
2777            xnelas=xres+xdirct 
2778            icheck=1
2779            go to 34
2780 * For pion + baryon resonance the reabsorption 
2781 * cross section is calculated from the detailed balance
2782 * using reab(i1,i2,srt,ictrl), ictrl=1, 2 and 3
2783 * for pion, rho and omega + baryon resonance
2784 31           ec=(em1+em2+0.02)**2
2785            xreab=reab(i1,i2,srt,1)
2786
2787 clin-12/02/00 to satisfy detailed balance, forbid N* absorptions:
2788           if((iabs(lb1).ge.10.and.iabs(lb1).le.13)
2789      1         .or.(iabs(lb2).ge.10.and.iabs(lb2).le.13)) XREAB = 0.
2790
2791            xkaon=xkaon0+xreab
2792 * a constant of 10 mb IS USED FOR PION + N* RESONANCE, 
2793         IF((iabs(LB1).GT.9.AND.iabs(LB1).LE.13) .OR.
2794      &      (iabs(LB2).GT.9.AND.iabs(LB2).LE.13))THEN
2795            Xnelas=1.0
2796         ELSE
2797            XNELAS=DPION(EM1,EM2,LB1,LB2,SRT)
2798         ENDIF
2799            icheck=2
2800 34          IF((Xnelas+xkaon+Xphi).LE.0.000001)GO TO 400
2801             DS=SQRT((Xnelas+xkaon+Xphi)/PI)
2802 csp09/20/01
2803 c           totcr = xnelas+xkaon
2804 c           if(srt .gt. 3.5)totcr = max1(totcr,3.)
2805 c           DS=SQRT(totcr/PI)
2806 csp09/20/01 end
2807             
2808            deltar=ds+0.1
2809          CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
2810      1   PCX,PCY,PCZ)
2811          IF(IC.EQ.-1) GO TO 400
2812        ekaon(4,iss)=ekaon(4,iss)+1
2813 c***
2814 * check what kind of collision has happened
2815 * (1) pion+baryon resonance
2816 * if direct elastic process
2817         if(icheck.eq.2)then
2818 c  !!sp11/21/01
2819       if(xnelas/(xnelas+xkaon+Xphi).ge.RANART(NSEED))then
2820 c               call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2)
2821                call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2822               go to 440
2823               else
2824 * for inelastic process, go to 96 to check
2825 * kaon production and pion reabsorption : pion+D(N*)-->pion+N
2826                go to 96
2827                 endif
2828               endif
2829 *(2) pion+n
2830 * CHECK IF inELASTIC COLLISION IS POSSIBLE FOR PION+N COLLISIONS
2831 clin-8/17/00 typo corrected, many other occurences:
2832 c        IF(XKAON/(XKAON+Xnelas).GT.RANART(NSEED))GO TO 95
2833        IF((XKAON+Xphi)/(XKAON+Xphi+Xnelas).GT.RANART(NSEED))GO TO 95
2834
2835 * direct process
2836         if(xdirct/xnelas.ge.RANART(NSEED))then
2837 c               call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2)
2838                call Crdir(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,IBLOCK)
2839               go to 440
2840               endif
2841 * now resonance formation or direct process (higher resonances)
2842            IF( (LB1*LB2.EQ.5.OR.((LB1*LB2.EQ.6).AND.
2843      &         (LB1.EQ.3.OR.LB2.EQ.3)))
2844      &    .OR. (LB1*LB2.EQ.-3.OR.((LB1*LB2.EQ.-10).AND.
2845      &         (LB1.EQ.5.OR.LB2.EQ.5))) )then    
2846 c
2847 * ONLY DELTA RESONANCE IS POSSIBLE, go to 99
2848         GO TO 99
2849        else
2850 * NOW BOTH DELTA AND N* RESORANCE ARE POSSIBLE
2851 * DETERMINE THE RESORANT STATE BY USING THE MONTRE CARLO METHOD
2852             XX=(XNPIN+xnpin1)/xres
2853             IF(RANART(NSEED).LT.XX)THEN
2854 * N* RESONANCE IS SELECTED
2855 * decide N*(1440) or N*(1535) formation
2856         xx0=xnpin/(xnpin+xnpin1)
2857         if(RANART(NSEED).lt.xx0)then
2858          RESONA=0.
2859 * N*(1440) formation
2860          GO TO 97
2861         else
2862 * N*(1535) formation
2863         resona=1.
2864          GO TO 98
2865         endif
2866          ELSE
2867 * DELTA RESONANCE IS SELECTED
2868          GO TO 99
2869          ENDIF
2870          ENDIF
2871 97       CONTINUE
2872             IF(RESONA.EQ.0.)THEN
2873 *N*(1440) IS PRODUCED,WE DETERMINE THE CHARGE STATE OF THE PRODUCED N*
2874             I=I1
2875             IF(EM1.LT.0.6)I=I2
2876 * (0.1) n+pion(+)-->N*(+)
2877            IF( (LB1*LB2.EQ.10.AND.(LB1.EQ.5.OR.LB2.EQ.5))
2878      &      .OR.(LB1*LB2.EQ.-6.AND.(LB1.EQ.3.OR.LB2.EQ.3)) )THEN
2879             LB(I)=11
2880            go to 303
2881             ENDIF
2882 * (0.2) p+pion(0)-->N*(+)
2883 c            IF(LB(I1)*LB(I2).EQ.4.AND.(LB(I1).EQ.1.OR.LB(I2).EQ.1))THEN
2884             IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.
2885      &         (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN    
2886             LB(I)=11
2887            go to 303
2888             ENDIF
2889 * (0.3) n+pion(0)-->N*(0)
2890 c            IF(LB(I1)*LB(I2).EQ.8.AND.(LB(I1).EQ.2.OR.LB(I2).EQ.2))THEN
2891             IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.
2892      &        (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN    
2893             LB(I)=10
2894            go to 303
2895             ENDIF
2896 * (0.4) p+pion(-)-->N*(0)
2897 c            IF(LB(I1)*LB(I2).EQ.3)THEN
2898             IF( (LB(I1)*LB(I2).EQ.3)
2899      &      .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
2900             LB(I)=10
2901             ENDIF
2902 303         CALL DRESON(I1,I2)
2903             if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I)
2904             lres=lres+1
2905             GO TO 101
2906 *COM: GO TO 101 TO CHANGE THE PHASE SPACE DENSITY OF THE NUCLEON
2907             ENDIF
2908 98          IF(RESONA.EQ.1.)THEN
2909 *N*(1535) IS PRODUCED, WE DETERMINE THE CHARGE STATE OF THE PRODUCED N*
2910             I=I1
2911             IF(EM1.LT.0.6)I=I2
2912 * note: this condition applies to both eta and pion
2913 * (0.1) n+pion(+)-->N*(+)
2914 c            IF(LB1*LB2.EQ.10.AND.(LB1.EQ.2.OR.LB2.EQ.2))THEN
2915             IF( (LB1*LB2.EQ.10.AND.(LB1.EQ.5.OR.LB2.EQ.5))
2916      &      .OR.(LB1*LB2.EQ.-6.AND.(LB1.EQ.3.OR.LB2.EQ.3)) )THEN
2917             LB(I)=13
2918            go to 304
2919             ENDIF
2920 * (0.2) p+pion(0)-->N*(+)
2921 c            IF(LB(I1)*LB(I2).EQ.4.AND.(LB(I1).EQ.1.OR.LB(I2).EQ.1))THEN
2922             IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.
2923      &           (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN 
2924             LB(I)=13
2925            go to 304
2926             ENDIF
2927 * (0.3) n+pion(0)-->N*(0)
2928 c            IF(LB(I1)*LB(I2).EQ.8.AND.(LB(I1).EQ.2.OR.LB(I2).EQ.2))THEN
2929             IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.
2930      &           (LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN      
2931             LB(I)=12
2932            go to 304
2933             ENDIF
2934 * (0.4) p+pion(-)-->N*(0)
2935 c            IF(LB(I1)*LB(I2).EQ.3)THEN
2936             IF( (LB(I1)*LB(I2).EQ.3)
2937      &      .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
2938             LB(I)=12
2939            go to 304
2940            endif
2941 * (0.5) p+eta-->N*(+)(1535),n+eta-->N*(0)(1535)
2942            if(lb(i1)*lb(i2).eq.0)then
2943 c            if((lb(i1).eq.1).or.(lb(i2).eq.1))then
2944             if(iabs(lb(i1)).eq.1.or.iabs(lb(i2)).eq.1)then
2945            LB(I)=13
2946            go to 304
2947            ELSE
2948            LB(I)=12
2949            ENDIF
2950            endif
2951 304         CALL DRESON(I1,I2)
2952             if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I) 
2953             lres=lres+1
2954             GO TO 101
2955 *COM: GO TO 101 TO CHANGE THE PHASE SPACE DENSITY OF THE NUCLEON
2956             ENDIF
2957 *DELTA IS PRODUCED,IN THE FOLLOWING WE DETERMINE THE
2958 *CHARGE STATE OF THE PRODUCED DELTA
2959 99      LRES=LRES+1
2960         I=I1
2961         IF(EM1.LE.0.6)I=I2
2962 * (1) p+pion(+)-->DELTA(++)
2963 c        IF(LB(I1)*LB(I2).EQ.5)THEN
2964             IF( (LB(I1)*LB(I2).EQ.5)
2965      &      .OR.(LB(I1)*LB(I2).EQ.-3) )THEN
2966         LB(I)=9
2967        go to 305
2968         ENDIF
2969 * (2) p+pion(0)-->delta(+)
2970 c        IF(LB(I1)*LB(I2).EQ.4.AND.(LB(I1).EQ.1.OR.LB(I2).EQ.1))then
2971        IF(iabs(LB(I1)*LB(I2)).EQ.4.AND.(LB(I1).EQ.4.OR.LB(I2).EQ.4))then
2972         LB(I)=8
2973        go to 305
2974         ENDIF
2975 * (3) n+pion(+)-->delta(+)
2976 c        IF(LB(I1)*LB(I2).EQ.10.AND.(LB(I1).EQ.2.OR.LB(I2).EQ.2))THEN
2977        IF( (LB(I1)*LB(I2).EQ.10.AND.(LB(I1).EQ.5.OR.LB(I2).EQ.5))
2978      & .OR.(LB(I1)*LB(I2).EQ.-6.AND.(LB(I1).EQ.3.OR.LB(I2).EQ.3)) )THEN
2979         LB(I)=8
2980        go to 305
2981         ENDIF
2982 * (4) n+pion(0)-->delta(0)
2983 c        IF(LB(I1)*LB(I2).EQ.8.AND.(LB(I1).EQ.2.OR.LB(I2).EQ.2))THEN
2984        IF(iabs(LB(I1)*LB(I2)).EQ.8.AND.(LB(I1).EQ.4.OR.LB(I2).EQ.4))THEN
2985         LB(I)=7
2986        go to 305
2987         ENDIF
2988 * (5) p+pion(-)-->delta(0)
2989 c        IF(LB(I1)*LB(I2).EQ.3)THEN
2990             IF( (LB(I1)*LB(I2).EQ.3)
2991      &      .OR.(LB(I1)*LB(I2).EQ.-5) )THEN
2992         LB(I)=7
2993        go to 305
2994         ENDIF
2995 * (6) n+pion(-)-->delta(-)
2996 c        IF(LB(I1)*LB(I2).EQ.6.AND.(LB(I1).EQ.2.OR.LB(I2).EQ.2))THEN
2997        IF( (LB(I1)*LB(I2).EQ.6.AND.(LB(I1).EQ.3.OR.LB(I2).EQ.3))
2998      & .OR.(LB(I1)*LB(I2).EQ.-10.AND.(LB(I1).EQ.5.OR.LB(I2).EQ.5)) )THEN 
2999         LB(I)=6
3000         ENDIF
3001 305     CALL DRESON(I1,I2)
3002         if(LB1.lt.0.or.LB2.lt.0) LB(I)=-LB(I) 
3003        GO TO 101
3004
3005 csp-11/08/01 K*
3006 * FOR kaON+pion COLLISIONS, form K* (bar) or
3007 c La/Si-bar + N <-- pi + K+
3008 c La/Si + N-bar <-- pi + K-                                             
3009 c phi + K <-- pi + K                                             
3010 clin (rho,omega) + K* <-- pi + K
3011 889       CONTINUE
3012         PX1CM=PCX
3013         PY1CM=PCY
3014         PZ1CM=PCZ
3015         EC=(em1+em2+0.02)**2
3016 * the cross section is from C.M. Ko, PRC 23, 2760 (1981).
3017        spika=60./(1.+4.*(srt-0.895)**2/(0.05)**2)
3018 c
3019 cc       if(lb(i1).eq.23.or.lb(i2).eq.23)then   !! block  K- + pi->La + B-bar
3020
3021         call Crkpla(PX1CM,PY1CM,PZ1CM,EC,SRT,spika,
3022      &                  emm1,emm2,lbp1,lbp2,I1,I2,icase,srhoks)
3023 cc
3024 c* only K* or K*bar formation
3025 c       else 
3026 c      DSkn=SQRT(spika/PI/10.)
3027 c      dsknr=dskn+0.1
3028 c      CALL DISTCE(I1,I2,dsknr,DSkn,DT,EC,SRT,IC,
3029 c    1     PX1CM,PY1CM,PZ1CM)
3030 c        IF(IC.EQ.-1) GO TO 400
3031 c       icase = 1
3032 c      endif
3033 c
3034          if(icase .eq. 0) then
3035             iblock=0
3036             go to 400
3037          endif
3038
3039        if(icase .eq. 1)then
3040              call KSRESO(I1,I2)
3041 clin-4/30/03 give non-zero iblock for resonance selections:
3042              iblock = 171
3043 ctest off for resonance (phi, K*) studies:
3044 c             if(iabs(lb(i1)).eq.30) then
3045 c             write(17,112) 'ks',lb(i1),p(1,i1),p(2,i1),p(3,i1),e(i1),nt
3046 c             elseif(iabs(lb(i2)).eq.30) then
3047 c             write(17,112) 'ks',lb(i2),p(1,i2),p(2,i2),p(3,i2),e(i2),nt
3048 c             endif
3049 c
3050               lres=lres+1
3051               go to 101
3052        elseif(icase .eq. 2)then
3053              iblock = 71
3054 c
3055 * La/Si (bar) formation                                                   
3056
3057        elseif(iabs(icase).eq.5)then
3058              iblock = 88
3059
3060        else
3061 *
3062 * phi formation
3063              iblock = 222
3064        endif
3065              LB(I1) = lbp1
3066              LB(I2) = lbp2
3067              E(I1) = emm1
3068              E(I2) = emm2
3069              em1=e(i1)
3070              em2=e(i2)
3071              ntag = 0
3072              go to 440
3073 c             
3074 33       continue
3075        em1=e(i1)
3076        em2=e(i2)
3077 * (1) if rho or omega collide with a nucleon we allow both elastic 
3078 *     scattering and kaon production to happen if collision conditions 
3079 *     are satisfied.
3080 * (2) if rho or omega collide with a baryon resonance we allow
3081 *     kaon production, pion reabsorption: rho(omega)+D(N*)-->pion+N
3082 *     and NO elastic scattering to happen
3083            xelstc=0
3084             if((lb1.ge.25.and.lb1.le.28).and.
3085      &    (iabs(lb2).eq.1.or.iabs(lb2).eq.2))
3086      &      xelstc=ERHON(EM1,EM2,LB1,LB2,SRT)
3087             if((lb2.ge.25.and.lb2.le.28).and.
3088      &   (iabs(lb1).eq.1.or.iabs(lb1).eq.2))
3089      &      xelstc=ERHON(EM1,EM2,LB1,LB2,SRT)
3090             ec=(em1+em2+0.02)**2
3091 * the kaon production cross section is
3092            xkaon0=0
3093            if(srt.ge.1.63.AND.SRT.LE.1.7)xkaon0=pnlka(srt)
3094            IF(SRT.GT.1.7)XKAON0=PNLKA(SRT)+pnska(srt)
3095            if(xkaon0.lt.0)xkaon0=0
3096
3097 cbz3/7/99 neutralk
3098             XKAON0 = 2.0 * XKAON0
3099 cbz3/7/99 neutralk end
3100
3101 * the total inelastic cross section for rho(omega)+N is
3102            xkaon=xkaon0
3103            ichann=0
3104 * the total inelastic cross section for rho (omega)+D(N*) is 
3105 * xkaon=xkaon0+reab(**) 
3106
3107 c sp11/21/01  phi production: rho + N(D) -> phi + N(D)
3108          Xphi = 0.
3109        if( ( (((lb1.ge.1.and.lb1.le.2).or.
3110      &         (lb1.ge.6.and.lb1.le.9))
3111      &         .and.(lb2.ge.25.and.lb2.le.27))
3112      &   .OR.(((lb2.ge.1.and.lb2.le.2).or.
3113      &         (lb2.ge.6.and.lb2.le.9))
3114      &        .and.(lb1.ge.25.and.lb1.le.27)) ).AND. srt.gt.1.958)
3115      &    call pibphi(srt,lb1,lb2,em1,em2,Xphi,xphin)
3116 c !! in fm^2 above
3117 c
3118         if((iabs(lb1).ge.6.and.lb2.ge.25).or.
3119      &    (lb1.ge.25.and.iabs(lb2).ge.6))then
3120            ichann=1
3121            ictrl=2
3122            if(lb1.eq.28.or.lb2.eq.28)ictrl=3
3123             xreab=reab(i1,i2,srt,ictrl)
3124
3125 clin-12/02/00 to satisfy detailed balance, forbid N* absorptions:
3126             if((iabs(lb1).ge.10.and.iabs(lb1).le.13)
3127      1           .or.(iabs(lb2).ge.10.and.iabs(lb2).le.13)) XREAB = 0.
3128
3129         if(xreab.lt.0)xreab=1.E-06
3130             xkaon=xkaon0+xreab
3131           XELSTC=1.0
3132            endif
3133             DS=SQRT((XKAON+Xphi+xelstc)/PI)
3134 c
3135 csp09/20/01
3136 c           totcr = xelstc+xkaon
3137 c           if(srt .gt. 3.5)totcr = max1(totcr,3.)
3138 c           DS=SQRT(totcr/PI)
3139 csp09/20/01 end
3140 c
3141         DELTAR=DS+0.1
3142        px1cm=pcx
3143        py1cm=pcy
3144        pz1cm=pcz
3145 * CHECK IF the collision can happen
3146          CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
3147      1   PCX,PCY,PCZ)
3148          IF(IC.EQ.-1) GO TO 400
3149         ekaon(4,iss)=ekaon(4,iss)+1
3150 c*
3151 * NOW rho(omega)+N or D(N*) COLLISION IS POSSIBLE
3152 * (1) check elastic collision
3153        if(xelstc/(xelstc+xkaon+Xphi).gt.RANART(NSEED))then
3154 c       call crdir(px1CM,py1CM,pz1CM,srt,I1,i2)
3155        call crdir(px1CM,py1CM,pz1CM,srt,I1,i2,IBLOCK)
3156        go to 440
3157        endif
3158 * (2) check pion absorption or kaon production
3159         CALL CRRD(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3160      1  IBLOCK,xkaon0,xkaon,Xphi,xphin)
3161
3162 * kaon production
3163 csp05/16/01
3164        IF(IBLOCK.EQ.7) then
3165           LPN=LPN+1
3166        elseIF(IBLOCK.EQ.-7) then
3167        endif
3168 csp05/16/01 end
3169 * rho obsorption
3170        if(iblock.eq.81) lrhor=lrhor+1
3171 * omega obsorption
3172        if(iblock.eq.82) lomgar=lomgar+1
3173        em1=e(i1)
3174        em2=e(i2)
3175        GO TO 440
3176 * for pion+n now using the subroutine crpn to change 
3177 * the particle label and set the new momentum of L/S+K final state
3178 95       continue
3179 * NOW PION+N INELASTIC COLLISION IS POSSIBLE
3180 * check pion production or kaon production
3181         CALL CRPN(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3182      1  IBLOCK,xkaon0,xkaon,Xphi,xphin)
3183
3184 * kaon production
3185 csp05/16/01
3186        IF(IBLOCK.EQ.7) then
3187           LPN=LPN+1
3188        elseIF(IBLOCK.EQ.-7) then
3189        endif
3190 csp05/16/01 end
3191 * pion production
3192        if(iblock.eq.77) lpd=lpd+1
3193 * rho production
3194        if(iblock.eq.78) lrho=lrho+1
3195 * omega production
3196        if(iblock.eq.79) lomega=lomega+1
3197        em1=e(i1)
3198        em2=e(i2)
3199        GO TO 440
3200 * for pion+D(N*) now using the subroutine crpd to 
3201 * (1) check kaon production or pion reabsorption 
3202 * (2) change the particle label and set the new 
3203 *     momentum of L/S+K final state
3204 96       continue
3205         CALL CRPD(PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3206      1  IBLOCK,xkaon0,xkaon,Xphi,xphin)
3207
3208 * kaon production
3209 csp05/16/01
3210        IF(IBLOCK.EQ.7) then
3211           LPN=LPN+1
3212        elseIF(IBLOCK.EQ.-7) then
3213        endif
3214 csp05/16/01 end
3215 * pion obserption
3216        if(iblock.eq.80) lpdr=lpdr+1
3217        em1=e(i1)
3218        em2=e(i2)
3219        GO TO 440
3220 * CALCULATE KAON PRODUCTION PROBABILITY FROM PION + N COLLISIONS
3221 C        IF(SRT.GT.1.615)THEN
3222 C        CALL PKAON(SRT,XXp,PK)
3223 C        TKAON(7)=TKAON(7)+PK 
3224 C        EKAON(7,ISS)=EKAON(7,ISS)+1
3225 c        CALL KSPEC1(SRT,PK)
3226 C        call LK(3,srt,iseed,pk)
3227 C        ENDIF
3228 * negelecting the pauli blocking at high energies
3229
3230 101       continue
3231         IF(E(I2).EQ.0.)GO TO 600
3232         IF(E(I1).EQ.0.)GO TO 800
3233 * IF NUCLEON+BARYON RESONANCE COLLISIONS
3234 44      CONTINUE
3235 * CALCULATE THE TOTAL CROSS SECTION OF NUCLEON+ BARYON RESONANCE COLLISION
3236 * WE ASSUME THAT THE ELASTIC CROSS SECTION IS THE SAME AS NUCLEON+NUCLEON
3237 * COM: WE USE THE PARAMETERISATION BY CUGNON FOR LOW ENERGIES
3238 *      AND THE PARAMETERIZATIONS FROM CERN DATA BOOK FOR HIGHER 
3239 *      ENERGIES. THE CUTOFF FOR THE TOTAL CROSS SECTION IS 55 MB 
3240        cutoff=em1+em2+0.02
3241        IF(SRT.LE.CUTOFF)GO TO 400
3242         IF(SRT.GT.2.245)THEN
3243        SIGNN=PP2(SRT)
3244        ELSE
3245         SIGNN = 35.0 / (1. + (SRT - CUTOFF) * 100.0)  +  20.0
3246        ENDIF 
3247         call XND(pcx,pcy,pcz,srt,I1,I2,xinel,
3248      &               sigk,xsk1,xsk2,xsk3,xsk4,xsk5)
3249        sig=signn+xinel
3250 * For nucleon+baryon resonance collision, the minimum cms**2 energy is
3251         EC=(EM1+EM2+0.02)**2
3252 * CHECK THE DISTENCE BETWEEN THE TWO PARTICLES
3253         PX1CM=PCX
3254         PY1CM=PCY
3255         PZ1CM=PCZ
3256
3257 clin-6/2008 Deuteron production:
3258         ianti=0
3259         if(lb(i1).lt.0 .and. lb(i2).lt.0) ianti=1
3260         call sbbdm(srt,sdprod,ianti,lbm,xmm,pfinal)
3261         sig=sig+sdprod
3262 clin-6/2008 perturbative treatment of deuterons:
3263         ipdflag=0
3264         if(idpert.eq.1) then
3265            ipert1=1
3266            sigr0=sig
3267            dspert=sqrt(sigr0/pi/10.)
3268            dsrpert=dspert+0.1
3269            CALL DISTCE(I1,I2,dsrpert,dspert,DT,EC,SRT,IC,
3270      1          PX1CM,PY1CM,PZ1CM)
3271            IF(IC.EQ.-1) GO TO 363
3272            signn0=0.
3273            CALL CRND(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3274      &  IBLOCK,SIGNN0,SIGr0,sigk,xsk1,xsk2,xsk3,xsk4,xsk5,NT,ipert1)
3275 c     &  IBLOCK,SIGNN,SIG,sigk,xsk1,xsk2,xsk3,xsk4,xsk5)
3276            ipdflag=1
3277  363       continue
3278            ipert1=0
3279         endif
3280         if(idpert.eq.2) ipert1=1
3281 c
3282         DS=SQRT(SIG/(10.*PI))
3283         DELTAR=DS+0.1
3284         CALL DISTCE(I1,I2,DELTAR,DS,DT,EC,SRT,IC,
3285      1  PX1CM,PY1CM,PZ1CM)
3286 c        IF(IC.EQ.-1)GO TO 400
3287         IF(IC.EQ.-1) then
3288            if(ipdflag.eq.1) iblock=501
3289            GO TO 400
3290         endif
3291
3292         ekaon(3,iss)=ekaon(3,iss)+1
3293 * CALCULATE KAON PRODUCTION PROBABILITY FROM NUCLEON + BARYON RESONANCE 
3294 * COLLISIONS
3295         go to 361
3296
3297 * CHECK WHAT KIND OF COLLISION HAS HAPPENED
3298  361    continue 
3299         CALL CRND(IRUN,PX1CM,PY1CM,PZ1CM,SRT,I1,I2,
3300      &     IBLOCK,SIGNN,SIG,sigk,xsk1,xsk2,xsk3,xsk4,xsk5,NT,ipert1)
3301 c     &  IBLOCK,SIGNN,SIG,sigk,xsk1,xsk2,xsk3,xsk4,xsk5)
3302         IF(iblock.eq.0.and.ipdflag.eq.1) iblock=501
3303         IF(IBLOCK.EQ.11)THEN
3304            LNDK=LNDK+1
3305            GO TO 400
3306 c        elseIF(IBLOCK.EQ.-11) then
3307         elseIF(IBLOCK.EQ.-11.or.iblock.eq.501) then
3308            GO TO 400
3309         ENDIF
3310         if(iblock .eq. 222)then
3311 c    !! sp12/17/01 
3312            GO TO 400
3313         ENDIF
3314         em1=e(i1)
3315         em2=e(i2)
3316         GO TO 440
3317 * IF NUCLEON+NUCLEON OR BARYON RESONANCE+BARYON RESONANCE COLLISIONS
3318 4       CONTINUE
3319 * PREPARE THE EALSTIC CROSS SECTION FOR BARYON+BARYON COLLISIONS
3320 * COM: WE USE THE PARAMETERISATION BY CUGNON FOR SRT LEQ 2.0 GEV
3321 *      AND THE PARAMETERIZATIONS FROM CERN DATA BOOK FOR HIGHER 
3322 *      ENERGIES. THE CUTOFF FOR THE TOTAL CROSS SECTION IS 55 MB 
3323 *      WITH LOW-ENERGY-CUTOFF
3324         CUTOFF=em1+em2+0.14
3325 * AT HIGH ENERGIES THE ISOSPIN DEPENDENCE IS NEGLIGIBLE
3326 * THE TOTAL CROSS SECTION IS TAKEN AS THAT OF THE PP 
3327 * ABOVE E_KIN=800 MEV, WE USE THE ISOSPIN INDEPENDNET XSECTION
3328         IF(SRT.GT.2.245)THEN
3329            SIG=ppt(srt)
3330            SIGNN=SIG-PP1(SRT)
3331         ELSE
3332 * AT LOW ENERGIES THE ISOSPIN DEPENDENCE FOR NN COLLISION IS STRONG
3333            SIG=XPP(SRT)
3334            IF(ZET(LB(I1))*ZET(LB(I2)).LE.0)SIG=XNP(SRT)
3335            IF(ZET(LB(I1))*ZET(LB(I2)).GT.0)SIG=XPP(SRT)
3336            IF(ZET(LB(I1)).EQ.0.
3337      &          AND.ZET(LB(I2)).EQ.0)SIG=XPP(SRT)
3338            if((lb(i1).eq.-1.and.lb(i2).eq.-2) .or.
3339      &          (lb(i2).eq.-1.and.lb(i1).eq.-2))sig=xnp(srt)
3340 *     WITH LOW-ENERGY-CUTOFF
3341            IF (SRT .LT. 1.897) THEN
3342               SIGNN = SIG
3343            ELSE