GeneratorInterface/AMPTInterface/src/amptsub.f - master - githubmirror/cmssw

GeneratorInterface/AMPTInterface/src/amptsub.f
5 042 строки153 KB
Kurt
clean-up 1 of compiler warnings
08 ноя 2015, 16:45
08 ноя 2015, 16:459eb58b7
c....................amptsub.f
c.....this file contains 4 sections:
c.....1. ART subroutines;
c.....2. ART functions;
c.....3. ART block data;
c.....4. subprocesses borrowed from other codes.
c.....5. the previous artana.f
c.....6. the previous zpcsub.f
c.....7. subroutine getnp
c.....Note that Parts1-4 are the previous artsub.f
c
c=======================================================================
c.....subroutine to set up ART parameters and analysis files
c.....before looping different events
      SUBROUTINE ARTSET
c
      PARAMETER (AMU= 0.9383,nxymax=10001)
      double precision dpcoal,drcoal,ecritl
      INTEGER ZTA, ZPR
      common  /gg/      dx,dy,dz,dpx,dpy,dpz
clin-10/03/03 
c     "SAVE   " (without argument) is used for most subroutines and functions,
c     this is important for the success when using "f77" to compile:
cc      SAVE /gg/
      common  /zz/      zta,zpr
cc      SAVE /zz/
      COMMON  /RUN/     NUM
cc      SAVE /RUN/
      common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc      SAVE /input1/
      COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, 
     &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
cc      SAVE /INPUT2/
      COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
cc      SAVE /INPUT3/
      common /imulst/ iperts
cc      SAVE /imulst/
      common /coal/dpcoal,drcoal,ecritl
      common/anim/nevent,isoft,isflag,izpc
      common /para7/ ioscar,nsmbbbar,nsmmeson
      common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd,
     1     psembd,tmaxembd,phidecomp
      common/xyembed/nxyjet,xyjet(nxymax,2)
      SAVE   
clin-10/03/03  ecritl: local energy density below which a parton 
c     will freeze out (in GeV/fm^3), for improvements on string melting, 
c     not used in this version of AMPT:
clin-4/2008
c      data ecritl/1.d0/
      ecritl=1.d0
c
c     combine ART initialization into ampt.ini:
c     (Note that the following values are relics from the old ART structure)
c.....input parameter file
c      OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN')
c      READ (13, *) MASSTA, ZTA
      MASSTA=1
      ZTA=1
c      write(12,*) massta, zta, ' massta, zta'
c      READ (13, *) MASSPR, ZPR
      MASSPR=1
      ZPR=1
c      write(12,*) masspr, zpr, ' masspr, zpr'
c      READ (13, *) PLAB, IPLAB
      PLAB=14.6 
      IPLAB=2
c      write(12,*) plab, iplab, ' plab, iplab'
      if(iplab.eq.2)then
         elab=sqrt(plab**2+amu**2)-amu
      else
         elab=plab
      endif
      elab=elab*1000.
c      READ (13, *) ZEROPT
      ZEROPT=0.
c      write(12,*) zeropt, ' zeropt'
clin-10/03/03 ISEED was used as a seed for random number inside ART, 
c     not used in AMPT:
      ISEED=700721
c     0/1: (Normal or Perturbative) multistrange partice production.
c     Perturbative option is disabled for now:
      iperts=0
c      READ (13, *) MANYB, B0, BI, BM
c     2/04/00 MANYB MUST BE SET TO 1 !
c     in order to skip impact parameter setting by ART, then B0 has no effect.
      MANYB=1
      B0=1
      BI=0
      BM=0
c      write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm'
c      READ (13, *) ISEED
c      write(12,*) iseed, ' iseed'
c      READ (13, *) DT
c      write(12,*) dt, ' dt'
c      READ (13, *) NTMAX
c      write(12,*) ntmax, ' ntmax'
c      READ (13, *) ICOLL
      ICOLL=-1
c      write(12,*) icoll, ' icoll'
c      READ (13, *) NUM
c     2/11/03 run events without test particles for now:
      NUM=1
c      write(12,*) num, ' num'
c      READ (13, *) INSYS
      INSYS=1
c      write(12,*) insys, ' insys'
c      READ (13, *) IPOT
      IPOT=3
c      write(12,*) ipot, ' ipot'
c      READ (13, *) MODE
      MODE=0
      IF(ICOLL.EQ.-1)IPOT=0
c      write(12,*) mode, ' mode'
c      READ (13, *) DX, DY, DZ
      DX=2.73
      DY=2.73
      DZ=2.73
c      write(12,*) dx,dy,dz,' dx,dy,dz'
c      READ (13, *) DPX, DPY, DPZ
      DPX=0.6
      DPY=0.6
      DPZ=0.6
c      write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz'
c      READ (13, *) IAVOID
      IAVOID=1
c      write(12,*) iavoid, ' iavoid'
c      READ (13, *) IMOMEN
      IMOMEN=1
c      write(12,*) imomen, ' imomen'
      if(icoll.eq.-1)imomen=3
c      READ (13, *) NFREQ
      NFREQ=10
c      write(12,*) nfreq, ' nfreq'
c      READ (13, *) ICFLOW
      ICFLOW=0
c      write(12,*) ICFLOW, ' ICFLOW'
c      READ (13, *) ICRHO
      ICRHO=0
c      write(12,*) ICRHO, ' ICRHO'
c      READ (13, *) ICOU
      ICOU=0
c      write(12,*)icou, ' icou'
* kaon potential control parameter
* KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH
c      READ (13, *) KPOTEN, KMUL
      KPOTEN=0
      KMUL=1
c      write(12,*)kpoten,kmul, ' kpoten, kmul'
* mean field control parameter FOR BARYONS
* no mean filed is used for baryons if their 
* local density is higher than dencut. 
c      READ (13, *) DENCUT
      DENCUT=15
c      write(12,*)dencut, ' dencut'
* test reactions in a box of side-length cycbox
* input cycbox
c      READ (13, *) CYCBOX
      CYCBOX=0
c      write(12,*) cycbox, ' cycbox'
c
clin-5b/2008
c      if(ioscar.eq.2) then
      if(ioscar.eq.2.or.ioscar.eq.3) then
c         OPEN (92,FILE='ana/parton-initial-afterPropagation.dat',
c     1        STATUS = 'UNKNOWN')
      endif
      if(ioscar.eq.3) then
clin-6/2009 write out full parton collision history:
c         OPEN (95,FILE='ana/parton-collisionsHistory.dat',
c     1        STATUS='UNKNOWN')
clin-6/2009 write out initial minijet information:
c         OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat',
c     1        STATUS='UNKNOWN')
clin-6/2009 write out parton info after coalescence:
         if(isoft.eq.4.or.isoft.eq.5) then
c            OPEN (85,FILE='ana/parton-after-coalescence.dat',
c     1           STATUS='UNKNOWN')
         endif
      endif
clin-6/2009 write out initial transverse positions of initial nucleons:
c      OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN')
c
clin-8/2009 In case that random positions are used to embed high-Pt jets:
      if(iembed.eq.3.or.iembed.eq.4) then
c         OPEN (97,FILE='embed-jet-xy.txt',STATUS = 'UNKNOWN')
c         read(97,*) nxyjet
c     Save positions in array to reuse when embedding more jet pairs 
c     than the number of entries in the position file:
         if(nevent.gt.nxyjet) then
            if(nxyjet.gt.nxymax) then
c              print *, 'Too many lines in embed-jet-xy.txt: 
c     1 increase value of the parameter nxymax'
               stop
            elseif(nxyjet.le.0) then
c               print *, 'Check number of entries in embed-jet-xy.txt' 
               stop
            endif
            do ixy=1,nxyjet
c               read(97,*) xyjet(ixy,1),xyjet(ixy,2)
            enddo
         endif
      endif

      RETURN
      END

c-----------------------------------------------------------------------

c.....subroutine to initialize cascade.

      SUBROUTINE ARINI

c.....before invoking ARINI:
c.....IAPAR2(1), IAINT2(1) must be set.
      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc      SAVE /ARPRNT/
      SAVE   

ctest off for resonance (phi, K*) studies:
c      OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN')

      IFLG = IAPAR2(1)
      GOTO (200, 200, 300) IFLG

c.....error choice of initialization
      PRINT *, 'IAPAR2(1) must be 1, 2, or 3'
      STOP

c.....to use default initial conditions generated by the cascade,
c.....or to read in initial conditions.
 200  RETURN

c.....to generate formation time and the position at formation time from 
c.....read-in initial conditions with an averaged formation proper time.
 300  CALL ARINI1
c.....ordering the particle label according to increasing order of 
c.....formation time.
      CALL ARTORD
      RETURN

      END

c-----------------------------------------------------------------------

c.....subroutine to generate formation time and position at formation time
c.....from read-in initial conditions with an averaged formation proper 
c.....time.

      SUBROUTINE ARINI1

c.....before invoking ARINI1:
c.....ARPAR1(1), IAINT2(1) must be set:
      PARAMETER (MAXSTR=150001)
      double precision  smearp,smearh

      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc      SAVE /ARPRNT/
      COMMON /ARPRC/ ITYPAR(MAXSTR),
     &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
     &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
     &     XMAR(MAXSTR)
cc      SAVE /ARPRC/
      COMMON /smearz/smearp,smearh
cc      SAVE /smearz/
      common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc      SAVE /input1/
      common/anim/nevent,isoft,isflag,izpc
cc      SAVE /anim/
      common /nzpc/nattzp
cc      SAVE /nzpc/
      COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc      SAVE /HPARNT/
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      common /para8/ idpert,npertd,idxsec
      SAVE   

clin-5/2008 for perturbatively-produced hadrons (currently only deuterons):
c      OPEN (91, FILE = 'ana/deuteron_processes.dat', 
c     1     STATUS = 'UNKNOWN')
      if(idpert.eq.1.or.idpert.eq.2) then
c         OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN')
      endif
c.....generate formation time and position at formation time.
      TAU0 = ARPAR1(1)
      NP = IAINT2(1)
clin-7/10/01     initial positions already given for hadrons 
c     formed from partons inside ZPC (from string melting):
      if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
         if(NP.le.nattzp) return
         do 1001 I = nattzp+1, NP
clin-9/2012 determine rapidity more generally 
c     to prevent overflow when Pt~=0 and E=|Pz|:
c            IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
c               PRINT *, ' IN ARINI1'
c               PRINT *, 'ABS(PZ) .GE. EE for particle ', I
c               PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), 
c     &              ' PY = ', PYAR(I)
c               PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
c               PRINT *, ' XM = ', XMAR(I)
c               RAP = 1000000.0
c               GOTO 50
c            END IF
cc            RAP=0.5*LOG((PEAR(I)+PZAR(I))/(PEAR(I)-PZAR(I)))
c            RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
c 50         CONTINUE
            if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then
               RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
            else
               PRINT *, ' IN ARINI1 mt=0'
               RAP = 1000000.0*sign(1.,PZAR(I))
            endif

            VX = PXAR(I) / PEAR(I)
            VY = PYAR(I) / PEAR(I)
            FTAR(I) = TAU0 * COSH(RAP)
            GXAR(I) = GXAR(I) + VX * FTAR(I)
            GYAR(I) = GYAR(I) + VY * FTAR(I)
            GZAR(I) = TAU0 * SINH(RAP)
clin-5/2009 No formation time for spectator projectile or target nucleons:
            if(PXAR(I).eq.0.and.PYAR(I).eq.0
     2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
clin-2/2013 for spectator target nucleons in LAB frame:
c     1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
              if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01)
     1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then
c
               TAUI=1.E-20
               FTAR(I)=TAUI*COSH(RAP)
               GZAR(I)=TAUI*SINH(RAP)
               endif
            endif
 1001    continue
clin-7/10/01-end
clin-3/2009 cleanup of program flow:
      else
         DO 1002 I = 1, NP
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
c               PRINT *, ' IN ARINI1'
c               PRINT *, 'ABS(PZ) .GE. EE for particle ', I
c               PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I), 
c     &              ' PY = ', PYAR(I)
c               PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
c               PRINT *, ' XM = ', XMAR(I)
c               RAP = 1000000.0
c               GOTO 100
cc               STOP
c            END IF
c 100        CONTINUE
c            RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
            if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then
               RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
            else
               PRINT *, ' IN ARINI1 mt=0'
               RAP = 1000000.0*sign(1.,PZAR(I))
            endif

            VX = PXAR(I) / PEAR(I)
            VY = PYAR(I) / PEAR(I)
c.....give initial formation time shift
            TAUI = FTAR(I) + TAU0
            FTAR(I) = TAUI * COSH(RAP)
            GXAR(I) = GXAR(I) + VX * TAU0 * COSH(RAP)
            GYAR(I) = GYAR(I) + VY * TAU0 * COSH(RAP)
c     4/25/03: hadron z-position upon formation determined the same way as x,y:
            GZAR(I) = TAUI * SINH(RAP)
c     the old prescription:
c            GZAR(I) = GZAR(I) + TAU0 * SINH(RAP)
            zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
            GZAR(I)=GZAR(I)+zsmear
cbz1/28/99end
c     10/05/01 no formation time for spectator projectile or target nucleons:
            if(PXAR(I).eq.0.and.PYAR(I).eq.0
     2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
clin-2/2013 for spectator target nucleons in LAB frame:
c     1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
              if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01)
     1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then
c
clin-5/2008:
c               TAUI=0.00001
               TAUI=1.E-20
               FTAR(I)=TAUI*COSH(RAP)
               GZAR(I)=TAUI*SINH(RAP)+zsmear
               endif
            endif
 1002    CONTINUE
clin-3/2009 cleanup of program flow:
      endif

clin-3/2009 Add initial hadrons before the hadron cascade starts:
      call addhad

      RETURN
      END

c-----------------------------------------------------------------------

c.....subroutine to order particle labels according to increasing 
c.....formation time

      SUBROUTINE ARTORD

c.....before invoking ARTORD:
c.....IAINT2(1) must be set:
      PARAMETER (MAXSTR=150001,MAXR=1)
      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc      SAVE /ARPRNT/
      COMMON /ARPRC/ ITYPAR(MAXSTR),
     &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
     &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
     &     XMAR(MAXSTR)
cc      SAVE /ARPRC/
clin-3/2009 Take care of particle weights when user inserts initial hadrons:
      COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
     1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
     2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
      DIMENSION dptemp(MAXSTR)
c
      DIMENSION ITYP0(MAXSTR), 
     &   GX0(MAXSTR), GY0(MAXSTR), GZ0(MAXSTR), FT0(MAXSTR),
     &   PX0(MAXSTR), PY0(MAXSTR), PZ0(MAXSTR), EE0(MAXSTR),
     &   XM0(MAXSTR)
      DIMENSION INDX(MAXSTR)
      EXTERNAL ARINDX
      SAVE   
c
      NPAR = 0
      NP = IAINT2(1)
      DO 1001 I = 1, NP
         ITYP0(I) = ITYPAR(I)
         GX0(I) = GXAR(I)
         GY0(I) = GYAR(I)
         GZ0(I) = GZAR(I)
         FT0(I) = FTAR(I)
         PX0(I) = PXAR(I)
         PY0(I) = PYAR(I)
         PZ0(I) = PZAR(I)
         EE0(I) = PEAR(I)
         XM0(I) = XMAR(I)
clin-3/2009:
         dptemp(I) = dpertp(I)
 1001 CONTINUE
      CALL ARINDX(MAXSTR, NP, FT0, INDX)
      DO 1002 I = 1, NP
cbz12/3/98
c         IF (ITYP0(INDX(I)) .EQ. 211) THEN
c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN
c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
c     &      ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR.
c     &      ITYP0(INDX(I)) .EQ. 111) THEN
c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
c     &      ITYP0(INDX(I)) .EQ. 2112) THEN
         NPAR = NPAR + 1
c         ITYPAR(I) = ITYP0(INDX(I))
c         GXAR(I) = GX0(INDX(I))
c         GYAR(I) = GY0(INDX(I))
c         GZAR(I) = GZ0(INDX(I))
c         FTAR(I) = FT0(INDX(I))
c         PXAR(I) = PX0(INDX(I))
c         PYAR(I) = PY0(INDX(I))
c         PZAR(I) = PZ0(INDX(I))
c         PEAR(I) = EE0(INDX(I))
c         XMAR(I) = XM0(INDX(I))
         ITYPAR(NPAR) = ITYP0(INDX(I))
         GXAR(NPAR) = GX0(INDX(I))
         GYAR(NPAR) = GY0(INDX(I))
         GZAR(NPAR) = GZ0(INDX(I))
         FTAR(NPAR) = FT0(INDX(I))
         PXAR(NPAR) = PX0(INDX(I))
         PYAR(NPAR) = PY0(INDX(I))
         PZAR(NPAR) = PZ0(INDX(I))
         PEAR(NPAR) = EE0(INDX(I))
         XMAR(NPAR) = XM0(INDX(I))
clin-3/2009:
         dpertp(NPAR)=dptemp(INDX(I))
c         END IF
cbz12/3/98end
 1002 CONTINUE
      IAINT2(1) = NPAR
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....subroutine to copy individually generated particle record into
c.....particle record for many test particle runs.

      SUBROUTINE ARINI2(K)

      PARAMETER (MAXSTR=150001,MAXR=1)
      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc      SAVE /ARPRNT/
      COMMON /ARPRC/ ITYPAR(MAXSTR),
     &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
     &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
     &     XMAR(MAXSTR)
cc      SAVE /ARPRC/
      COMMON /ARERC1/MULTI1(MAXR)
cc      SAVE /ARERC1/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cc      SAVE /ARPRC1/
      COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
cc      SAVE /tdecay/
      common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc      SAVE /input1/
      COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE, 
     &     IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
cc      SAVE /INPUT2/
      COMMON/RNDF77/NSEED
      COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
     1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
     2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
cc      SAVE /RNDF77/
      SAVE   

      MULTI1(K) = IAINT2(1)
      DO 1001 I = 1, MULTI1(K)
         ITYP1(I, K) = ITYPAR(I)
         GX1(I, K) = GXAR(I)
         GY1(I, K) = GYAR(I)
         GZ1(I, K) = GZAR(I)
         FT1(I, K) = FTAR(I)
         PX1(I, K) = PXAR(I)
         PY1(I, K) = PYAR(I)
         PZ1(I, K) = PZAR(I)
         EE1(I, K) = PEAR(I)
         XM1(I, K) = XMAR(I)
clin-3/2009 hadron weights are initialized in addhad():
clin-5/2008 all hadrons not perturbatively-produced have the weight of 1:
c         dpp1(I,K)=1.
         dpp1(I,K)=dpertp(I)
 1001 CONTINUE

c     initialize final time of each particle to ntmax*dt except for 
c     decay daughters, which have values given by tfdcy() and >(ntmax*dt):
      do 1002 ip=1,MAXSTR
         tfdcy(ip)=NTMAX*DT
         tft(ip)=NTMAX*DT
 1002 continue
c
      do 1004 irun=1,MAXR
         do 1003 ip=1,MAXSTR
            tfdpi(ip,irun)=NTMAX*DT
 1003    continue
 1004 continue

      RETURN
      END

c=======================================================================

c.....function to convert PDG flavor code into ART flavor code.

      FUNCTION IARFLV(IPDG)

      common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc      SAVE /input1/
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      SAVE   

c.....anti-Delta-
      IF (IPDG .EQ. -1114) THEN
         IARFLV = -6
         RETURN
      END IF

c.....anti-Delta0
      IF (IPDG .EQ. -2114) THEN
         IARFLV = -7
         RETURN
      END IF

c.....anti-Delta+
      IF (IPDG .EQ. -2214) THEN
         IARFLV = -8
         RETURN
      END IF

c.....anti-Delta++
      IF (IPDG .EQ. -2224) THEN
         IARFLV = -9
         RETURN
      END IF

cbzdbg2/23/99
c.....anti-proton
      IF (IPDG .EQ. -2212) THEN
         IARFLV = -1
         RETURN
      END IF

c.....anti-neutron
      IF (IPDG .EQ. -2112) THEN
         IARFLV = -2
         RETURN
      END IF
cbzdbg2/23/99end

c.....eta
      IF (IPDG .EQ. 221) THEN
         IARFLV = 0
         RETURN
      END IF

c.....proton
      IF (IPDG .EQ. 2212) THEN
         IARFLV = 1
         RETURN
      END IF

c.....neutron
      IF (IPDG .EQ. 2112) THEN
         IARFLV = 2
         RETURN
      END IF

c.....pi-
      IF (IPDG .EQ. -211) THEN
         IARFLV = 3
         RETURN
      END IF

c.....pi0
      IF (IPDG .EQ. 111) THEN
         IARFLV = 4
         RETURN
      END IF

c.....pi+
      IF (IPDG .EQ. 211) THEN
         IARFLV = 5
         RETURN
      END IF

c.....Delta-
      IF (IPDG .EQ. 1114) THEN
         IARFLV = 6
         RETURN
      END IF

c.....Delta0
      IF (IPDG .EQ. 2114) THEN
         IARFLV = 7
         RETURN
      END IF

c.....Delta+
      IF (IPDG .EQ. 2214) THEN
         IARFLV = 8
         RETURN
      END IF

c.....Delta++
      IF (IPDG .EQ. 2224) THEN
         IARFLV = 9
         RETURN
      END IF

c.....Lambda
      IF (IPDG .EQ. 3122) THEN
         IARFLV = 14
         RETURN
      END IF

c.....Lambda-bar
      IF (IPDG .EQ. -3122) THEN
         IARFLV = -14
         RETURN
      END IF

c.....Sigma-
      IF (IPDG .EQ. 3112) THEN
         IARFLV = 15
         RETURN
      END IF

c.....Sigma-bar
      IF (IPDG .EQ. -3112) THEN
         IARFLV = -15
         RETURN
      END IF 

c.....Sigma0
      IF (IPDG .EQ. 3212) THEN
         IARFLV = 16
         RETURN
      END IF

c.....Sigma0-bar
      IF (IPDG .EQ. -3212) THEN
         IARFLV = -16
         RETURN
      END IF 

c.....Sigma+
      IF (IPDG .EQ. 3222) THEN
         IARFLV = 17
         RETURN
      END IF

c.....Sigma+ -bar
      IF (IPDG .EQ. -3222) THEN
         IARFLV = -17
         RETURN
      END IF 

c.....K-
      IF (IPDG .EQ. -321) THEN
         IARFLV = 21
         RETURN
      END IF

c.....K+
      IF (IPDG .EQ. 321) THEN
         IARFLV = 23
         RETURN
      END IF

c.....temporary entry for K0
      IF (IPDG .EQ. 311) THEN
         IARFLV = 23
         RETURN
      END IF

c.....temporary entry for K0bar
      IF (IPDG .EQ. -311) THEN
         IARFLV = 21
         RETURN
      END IF

c.....temporary entry for K0S and K0L
      IF (IPDG .EQ. 310 .OR. IPDG .EQ. 130) THEN
         R = RANART(NSEED)
         IF (R .GT. 0.5) THEN
            IARFLV = 23
         ELSE
            IARFLV = 21
         END IF
         RETURN
      END IF

c.....rho-
      IF (IPDG .EQ. -213) THEN
         IARFLV = 25
         RETURN
      END IF

c.....rho0
      IF (IPDG .EQ. 113) THEN
         IARFLV = 26
         RETURN
      END IF

c.....rho+
      IF (IPDG .EQ. 213) THEN
         IARFLV = 27
         RETURN
      END IF

c.....omega
      IF (IPDG .EQ. 223) THEN
         IARFLV = 28
         RETURN
      END IF

c.....phi
      IF (IPDG .EQ. 333) THEN
         IARFLV = 29
         RETURN
      END IF

c.....K*+
      IF (IPDG .EQ. 323) THEN
         IARFLV = 30
         RETURN
      END IF
c.....K*-
      IF (IPDG .EQ. -323) THEN
         IARFLV = -30
         RETURN
      END IF
c.....temporary entry for K*0
      IF (IPDG .EQ. 313) THEN
         IARFLV = 30
         RETURN
      END IF
c.....temporary entry for K*0bar
      IF (IPDG .EQ. -313) THEN
         IARFLV = -30
         RETURN
      END IF

c...... eta-prime
      IF (IPDG .EQ. 331) THEN
         IARFLV = 31
         RETURN
      END IF
 
c...... a1
c     IF (IPDG .EQ. 777) THEN
c        IARFLV = 32
c        RETURN
c     END IF
                                
c... cascade-
      IF (IPDG .EQ. 3312) THEN
         IARFLV = 40
         RETURN
      END IF
 
c... cascade+ (bar)
      IF (IPDG .EQ. -3312) THEN
         IARFLV = -40
         RETURN
      END IF
 
c... cascade0
      IF (IPDG .EQ. 3322) THEN
         IARFLV = 41
         RETURN
      END IF
 
c... cascade0 -bar
      IF (IPDG .EQ. -3322) THEN
         IARFLV = -41
         RETURN
      END IF
 
c... Omega-
      IF (IPDG .EQ. 3334) THEN
         IARFLV = 45
         RETURN
      END IF 

c... Omega+ (bar)
      IF (IPDG .EQ. -3334) THEN
         IARFLV = -45
         RETURN
      END IF

c... Di-Omega
      IF (IPDG .EQ. 6666) THEN
         IARFLV = 44
         RETURN
      END IF
c sp06/05/01 end    

clin-3/2009 keep the same ID numbers in case there are initial deuterons:
      IF (IPDG .EQ. 42 .or. IPDG .EQ. -42) THEN
         IARFLV = IPDG
         RETURN
      END IF

c.....other
      IARFLV = IPDG + 10000

      RETURN
      END

c-----------------------------------------------------------------------

c.....function to convert ART flavor code into PDG flavor code.

      FUNCTION INVFLV(IART)

      common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
cc      SAVE /input1/
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      SAVE   

c.....anti-Delta-
      IF (IART .EQ. -6) THEN
         INVFLV = -1114
         RETURN
      END IF

c.....anti-Delta0
      IF (IART .EQ. -7) THEN
         INVFLV = -2114
         RETURN
      END IF

c.....anti-Delta+
      IF (IART .EQ. -8) THEN
         INVFLV = -2214
         RETURN
      END IF

c.....anti-Delta++
      IF (IART .EQ. -9) THEN
         INVFLV = -2224
         RETURN
      END IF

cbzdbg2/23/99
c.....anti-proton
      IF (IART .EQ. -1) THEN
         INVFLV = -2212
         RETURN
      END IF

c.....anti-neutron
      IF (IART .EQ. -2) THEN
         INVFLV = -2112
         RETURN
      END IF
cbzdbg2/23/99end

c.....eta
      IF (IART .EQ. 0) THEN
         INVFLV = 221
         RETURN
      END IF

c.....proton
      IF (IART .EQ. 1) THEN
         INVFLV = 2212
         RETURN
      END IF

c.....neutron
      IF (IART .EQ. 2) THEN
         INVFLV = 2112
         RETURN
      END IF

c.....pi-
      IF (IART .EQ. 3) THEN
         INVFLV = -211
         RETURN
      END IF

c.....pi0
      IF (IART .EQ. 4) THEN
         INVFLV = 111
         RETURN
      END IF

c.....pi+
      IF (IART .EQ. 5) THEN
         INVFLV = 211
         RETURN
      END IF

c.....Delta-
      IF (IART .EQ. 6) THEN
         INVFLV = 1114
         RETURN
      END IF

c.....Delta0
      IF (IART .EQ. 7) THEN
         INVFLV = 2114
         RETURN
      END IF

c.....Delta+
      IF (IART .EQ. 8) THEN
         INVFLV = 2214
         RETURN
      END IF

c.....Delta++
      IF (IART .EQ. 9) THEN
         INVFLV = 2224
         RETURN
      END IF

cc.....N*(1440), N*(1535) temporary entry
c      IF (IART .GE. 10 .AND. IART .LE.13) THEN
c         INVFLV = 0
c         RETURN
c      END IF

c.....Lambda
      IF (IART .EQ. 14) THEN
         INVFLV = 3122
         RETURN
      END IF
c.....Lambda-bar
      IF (IART .EQ. -14) THEN
         INVFLV = -3122
         RETURN
      END IF 

cbz3/12/99
c.....temporary entry for Sigma's
c      IF (IART .EQ. 15) THEN
c         R = RANART(NSEED)
c         IF (R .GT. 2. / 3.) THEN
c            INVFLV = 3112
c         ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN
c            INVFLV = 3212
c         ELSE
c            INVFLV = 3222
c         END IF
c         RETURN
c      END IF

c.....Sigma-
      IF (IART .EQ. 15) THEN
         INVFLV = 3112
         RETURN
      END IF

c.....Sigma- bar
      IF (IART .EQ. -15) THEN
         INVFLV = -3112
         RETURN
      END IF 

c.....Sigma0
      IF (IART .EQ. 16) THEN
         INVFLV = 3212
         RETURN
      END IF

c.....Sigma0 -bar
      IF (IART .EQ. -16) THEN
         INVFLV = -3212
         RETURN
      END IF

c.....Sigma+
      IF (IART .EQ. 17) THEN
         INVFLV = 3222
         RETURN
      END IF

c.....Sigma+ -bar
      IF (IART .EQ. -17) THEN
         INVFLV = -3222
         RETURN
      END IF 

clin-2/23/03 K0S and K0L are generated at the last timestep:
c.....temporary entry for K- and K0bar
      IF (IART .EQ. 21) THEN
c         R = RANART(NSEED)
c         IF (R .GT. 0.5) THEN
            INVFLV = -321
c         ELSE
c            INVFLV = -311
c            R = RANART(NSEED)
c            IF (R .GT. 0.5) THEN
c               INVFLV = 310
c            ELSE
c               INVFLV = 130
c            END IF
c         END IF
         RETURN
      END IF

c.....temporary entry for K+ and K0
      IF (IART .EQ. 23) THEN
c         R = RANART(NSEED)
c         IF (R .GT. 0.5) THEN
            INVFLV = 321
c         ELSE
c            INVFLV = 311
c            R = RANART(NSEED)
c            IF (R .GT. 0.5) THEN
c               INVFLV = 310
c            ELSE
c               INVFLV = 130
c            END IF
c         END IF
         RETURN
      END IF

c.....K0Long:
      IF (IART .EQ. 22) THEN
         INVFLV = 130
         RETURN
      ENDIF
c.....K0Short:
      IF (IART .EQ. 24) THEN
         INVFLV = 310
         RETURN
      ENDIF

c.....rho-
      IF (IART .EQ. 25) THEN
         INVFLV = -213
         RETURN
      END IF

c.....rho0
      IF (IART .EQ. 26) THEN
         INVFLV = 113
         RETURN
      END IF

c.....rho+
      IF (IART .EQ. 27) THEN
         INVFLV = 213
         RETURN
      END IF

c.....omega
      IF (IART .EQ. 28) THEN
         INVFLV = 223
         RETURN
      END IF

c.....phi
      IF (IART .EQ. 29) THEN
         INVFLV = 333
         RETURN
      END IF

c.....temporary entry for K*+ and K*0
      IF (IART .EQ. 30) THEN
         INVFLV = 323
         IF (RANART(NSEED).GT.0.5) INVFLV = 313
         RETURN
      END IF

c.....temporary entry for K*- and K*0bar
      IF (IART .EQ. -30) THEN
         INVFLV = -323
         IF (RANART(NSEED).GT.0.5) INVFLV = -313
         RETURN
      END IF

c... eta-prime (bar)
      IF (IART .EQ. 31) THEN
         INVFLV = 331
         RETURN
      END IF
 
c... a1
      IF (IART .EQ. 32) THEN
         INVFLV = 777
         RETURN
      END IF
 
c... cascade-
      IF (IART .EQ. 40) THEN
         INVFLV = 3312
         RETURN
      END IF                   

c... cascade+ (bar)
      IF (IART .EQ. -40) THEN
         INVFLV = -3312
         RETURN
      END IF
 
c... cascade0
      IF (IART .EQ. 41) THEN
         INVFLV = 3322
         RETURN
      END IF
 
c... cascade0 -bar
      IF (IART .EQ. -41) THEN
         INVFLV = -3322
         RETURN
      END IF
 
c... Omega-
      IF (IART .EQ. 45) THEN
         INVFLV = 3334
         RETURN
      END IF

c... Omega+ (bar)
      IF (IART .EQ. -45) THEN
         INVFLV = -3334
         RETURN
      END IF

c... Di-Omega
      IF (IART .EQ. 44) THEN
         INVFLV = 6666
         RETURN
      END IF
c sp 12/19/00 end           

clin-5/2008 deuteron ID numbers in ART and ampt.dat:
      IF (IART .EQ. 42) THEN
         INVFLV = 42
         RETURN
      ELSEIF (IART .EQ. -42) THEN         
         INVFLV = -42
         RETURN
      END IF
c
c.....other
      INVFLV = IART - 10000

      RETURN
      END

c=======================================================================

      BLOCK DATA ARDATA

      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
cc      SAVE /ARPRNT/
      SAVE   
      DATA ARPAR1/1.19, 99 * 0.0/
      DATA IAPAR2/3, 49 * 0/
      DATA ARINT1/100 * 0.0/
      DATA IAINT2/50 * 0/

      END

c=======================================================================

c.....Routine borrowed from ZPC.
c.....double precision  is modified to real*4.

cbz1/29/99
c      subroutine index1(n, m, arrin, indx)
      subroutine arindx(n, m, arrin, indx)
cbz1/29/99end
c     indexes the first m elements of ARRIN of length n, i.e., outputs INDX
c     such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m

c      implicit real*4 (a-h, o-z)

      dimension arrin(n), indx(n)
      SAVE   
      do 1001 j = 1, m
         indx(j) = j
 1001 continue
      l = m / 2 + 1
      ir = m
 10   continue
      if (l .gt. 1) then
         l = l - 1
         indxt = indx(l)
         q = arrin(indxt)
      else
         indxt = indx(ir)
         q = arrin(indxt)
         indx(ir) = indx(1)
         ir = ir - 1
         if (ir .eq. 1) then
            indx(1) = indxt
            return
         end if
      end if
      i = l
      j = l + l
 20   if (j .le. ir) then
         if (j .lt. ir) then
            if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
         end if
         if (q .lt. arrin(indx(j))) then
            indx(i) = indx(j)
            i = j
            j = j + j
         else
            j = ir + 1
         end if
      goto 20
      end if
      indx(i) = indxt
      goto 10

      end

c-----------------------------------------------------------------------

c.....extracted from G. Song's ART expasion including K- interactions
c.....file `NEWKAON.FOR'

c     5/01/03 send iblock value into art1f.f, necessary for resonance studies:
c        subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
c     &                                   srt,pcx,pcy,pcz)
        subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
     &                                   srt,pcx,pcy,pcz,iblock)
      PARAMETER      (MAXSTR=150001,MAXR=1)
      PARAMETER      (AKA=0.498)
      COMMON   /AA/  R(3,MAXSTR)
cc      SAVE /AA/
      COMMON   /BB/  P(3,MAXSTR)
cc      SAVE /BB/
      COMMON   /CC/  E(MAXSTR)
cc      SAVE /CC/
      COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
cc      SAVE /EE/
      COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
cc      SAVE /BG/
      COMMON   /NN/NNN
cc      SAVE /NN/
      COMMON   /RUN/NUM
cc      SAVE /RUN/
      COMMON   /PA/RPION(3,MAXSTR,MAXR)
cc      SAVE /PA/
      COMMON   /PB/PPION(3,MAXSTR,MAXR)
cc      SAVE /PB/
      COMMON   /PC/EPION(MAXSTR,MAXR)
cc      SAVE /PC/
      COMMON   /PD/LPION(MAXSTR,MAXR)
cc      SAVE /PD/
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      SAVE   
c
        logical lb1bn, lb2bn,lb1mn,lb2mn
cbz3/7/99 neutralk
c        logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0
        logical lb1bn1, lb2bn1, lb1bn0, lb2bn0
cbz3/7/99 neutralk end
        logical lb1mn0, lb2mn0, lb1mn1, lb2mn1
        logical lb1mn2, lb2mn2
        icase=-1
c        icase: flag for the type of reaction that is going to happen.
c        icase=-1,  no desired reaction, return to main program.
c              1,  NN,ND,DD
c              2,  PI+N, PI+D
c              3,  K(-) absorption.
        nchrg=-100
c        nchrg: Net charges of the two incoming particles.
        ictrl = 1
        lb1=lb(i1)
        lb2=lb(i2)
        em1=e(i1)
        em2=e(i2)
        lb1bn=lb1.eq.1.or.lb1.eq.2.or.(lb1.gt.5.and.lb1.le.13)
        lb2bn=lb2.eq.1.or.lb2.eq.2.or.(lb2.gt.5.and.lb2.le.13)
        lb1bn0=lb1.eq.2.or.lb1.eq.7.or.lb1.eq.10.or.lb1.eq.12
        lb2bn0=lb2.eq.2.or.lb2.eq.7.or.lb2.eq.10.or.lb2.eq.12
        lb1bn1=lb1.eq.1.or.lb1.eq.8.or.lb1.eq.11.or.lb1.eq.13
        lb2bn1=lb2.eq.1.or.lb2.eq.8.or.lb2.eq.11.or.lb2.eq.13
        lb1mn=em1.lt.0.2.or.lb1.eq.0.or.(lb1.ge.25.and.lb1.le.29)
        lb2mn=em2.lt.0.2.or.lb2.eq.0.or.(lb2.ge.25.and.lb2.le.29)
        lb1mn0=lb1.eq.0.or.lb1.eq.4.or.lb1.eq.26.or.
     &                        lb1.eq.28.or.lb1.eq.29
        lb2mn0=lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
     &                        lb2.eq.28.or.lb2.eq.29
        lb1mn1= lb1.eq.5.or.lb1.eq.27
        lb2mn1= lb2.eq.5.or.lb2.eq.27
        lb1mn2=lb1.eq.3.or.lb1.eq.25
        lb2mn2=lb2.eq.3.or.lb2.eq.25

c        1. consider N+N, N+Resonance, R + R reactions
        if(lb1bn.and.lb2bn) then
c     NN,ND,DD:
           icase=1
c     total cross section
           sig=40.
           if(lb1.eq.9.and.lb2.eq.9) then
                nchrg=4
           endif   
           if((lb1bn1.and.lb2.eq.9)
     &        .or.(lb2bn1.and.lb1.eq.9))then
                nchrg=3
           endif
           if((lb1bn0.and.lb2.eq.9)
     &        .or.(lb2bn0.and.lb1.eq.9)
     &        .or.(lb1bn1.and.lb2bn1)) then
                   nchrg=2
           endif
           if((lb1bn1.and.lb2bn0).or.(lb1.eq.6.and.lb2.eq.9)
     &        .or.(lb2bn1.and.lb1bn0)
     &        .or.(lb2.eq.6.and.lb1.eq.9))then
                   nchrg=1
           endif
           if((lb1bn0.and.lb2bn0).or.(lb1bn1.and.lb2.eq.6)
     &              .or.(lb2bn1.and.lb1.eq.6)) then
                   nchrg=0
           endif
           if((lb1bn0.and.lb2.eq.6)
     &        .or.(lb2bn0.and.lb1.eq.6))then
                nchrg=-1
           endif
           if(lb1.eq.6.and.lb2.eq.6) then
                nchrg=-2
           endif
c     brsig = x2kaon_no_isospin(srt)
           if(nchrg.ge.-1.and.nchrg.le.2) then
c     K,Kbar prduction x sect.
                   brsig = x2kaon(srt)
           else
                   brsig=0.0
c                if(nchrg.eq.-2.or.nchrg.eq.3) then
c                   brsig = x2kaon(srt+0.938-1.232)
c                else
c     nchrg=4
c                   brsig = x2kaon(srt+2.*(0.938-1.232))
c                endif
           endif

cbz3/7/99 neutralk
           BRSIG = 2.0 * BRSIG
cbz3/7/99 neutralk end

        endif

c        2. consider PI(meson:eta,omega,rho,phi) + N(N*,D)
        if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then
c     PN,PD
          icase=2
          sig=20.
          sigma0 = piNsg0(srt)
          brsig=0.0
          if((lb1bn1.and.lb2mn0)
     &       .or.(lb2bn1.and.lb1mn0).
     & or.(lb1bn0.and.lb2mn1).or.(lb2bn0.and.lb1mn1).
     & or.(lb1.eq.9.and.lb2mn2).or.(lb2.eq.9.and.lb1mn2))then
                nchrg=1
cbz3/2/99/song
c                if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0
c                if(lb1bn0.or.lb2bn0) brsig=0.5*sigma0
                if(lb1bn1.or.lb2bn1) brsig=0.5*sigma0
                if(lb1bn0.or.lb2bn0) brsig=2.0*sigma0
cbz3/2/99/song end
c                if(lb1.eq.9.or.lb2.eq.9) brsig=1.5*sigma0
          endif
          if( (lb1bn0.and.lb2mn0 )
     &       .or.(lb2bn0.and.lb1mn0)
     &  .or.(lb1bn1.and.lb2mn2).or.(lb2bn1.and.lb1mn2)
     &  .or.(lb1.eq.6.and.lb2mn1).or.(lb2.eq.6.and.lb1mn1)) then
                nchrg=0
                if(lb1bn1.or.lb2bn1) then
cbz3/2/99/song
c                  brsig=1.5*sigma0
                  brsig=3.0*sigma0
cbz3/2/99/song end
cbz3/11/99/song
c                  ratiok = 1./3.
                  ratiok = 2./3.
cbz3/11/99/song end

c                  ratiok: the ratio of channels: ->nK+k- vs. -> pK0K-
                endif
                if(lb1bn0.or.lb2bn0) then
                  brsig=2.5*sigma0
cbz3/2/99/song
c                  ratiok = 0.8
                  ratiok = 0.2
cbz3/2/99/song end
                endif
c                if(lb1.eq.6.or.lb2.eq.6) then
c     lb=6 : D-
c                  brsig=1.5*sigma0
c                  ratiok = 0.5
c                endif
          endif
          if( (lb1bn0.and.lb2mn2)
     &       .or.(lb2bn0.and.lb1mn2)
     & .or.(lb1.eq.6.and.lb2mn0).or.(lb2.eq.6.and.lb1mn0)) then
                nchrg=-1
                if(lb1bn0.or.lb2bn0) brsig=sigma0
c                if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0
          endif
c          if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then
c                nchrg=-2
c          endif
c          if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1)
c    &           .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
c                nchrg=2
c          endif

cbz3/11/99 neutralk
          if((lb1.eq.6.and.lb2mn2)
     &       .or.(lb2.eq.6.and.lb1mn2))then
                nchrg=-2
          endif
cbz3/11/99 neutralk
cbz3/8/99 neutralk
          if((lb1bn1.and.lb2mn1)
     &       .or.(lb2bn1.and.lb1mn1)
     & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
                nchrg=2
          endif
cbz3/8/99 neutralk end

cbz3/7/99 neutralk
          IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN
             BRSIG = 3.0 * SIGMA0
          END IF
cbz3/7/99 neutralk end

        endif

c        3. consider K- + N(N*,D) absorption.
c        if((lb1bn.and.lb2.eq.21).OR.(lb2bn.and.lb1.eq.21)) then
        if( (lb1bn.and.(lb2.eq.21.or.lb2.eq.-30)).OR.
     &     (lb2bn.and.(lb1.eq.21.or.lb1.eq.-30)) )then 
c          bmass=em1+em2-aka
          bmass=0.938
          if(srt.le.(bmass+aka)) then
cbz3/2/99
c                write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt
cbz3/2/99end
                pkaon=0.
          else
            pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
          endif
          sig=0.
          if(lb1.eq.1.or.lb2.eq.1.or.lb1.eq.8.or.lb2.eq.8.or.
     &    lb1.eq.11.or.lb2.eq.11.or.lb1.eq.13.or.lb2.eq.13) then
c          K- + (D+,N*+)p ->
              nchrg=0
              sigela=akPel(pkaon)
              sigsgm=3.*akPsgm(pkaon)
              sig=sigela+sigsgm+akPlam(pkaon)
          endif
          if(lb1.eq.2.or.lb2.eq.2.or.lb1.eq.7.or.lb2.eq.7.or.
     &    lb1.eq.10.or.lb2.eq.10.or.lb1.eq.12.or.lb2.eq.12) then
c          K- + (D0, N*0)n ->
              nchrg=-1
              sigela=akNel(pkaon)
              sigsgm=2.*akNsgm(pkaon)
              sig=sigela+sigsgm+akNlam(pkaon)
          endif
          if(lb1.eq.6.or.lb2.eq.6) then
c     K- + D-
              nchrg=-2
              sigela=akNel(pkaon)
              sigsgm=akNsgm(pkaon)
              sig=sigela+sigsgm
          endif
          if(lb1.eq.9.or.lb2.eq.9) then
c     K- + D++
              nchrg=1
              sigela=akPel(pkaon)
              sigsgm=2.*akPsgm(pkaon)
              sig=sigela+sigsgm+akPlam(pkaon)
          endif

cbz3/8/99 neutralk
          sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
          SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
          SIG = sigela + SIGSGM + AKPLAM(PKAON)
cbz3/8/99 neutralk end

          if(sig.gt.1.e-7) then
c     K(-) + N reactions
              icase=3
              brel=sigela/sig
              brsgm=sigsgm/sig
c              branch_lambda=akNlam(pkaon)/sig
              brsig = sig
          endif
        endif

c        4. meson + hyperon -> K- + N
c        if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR.
c     &     ((lb2.ge.14.and.lb2.le.17).and.lb1mn)) then
        if(((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.3.and.lb2.le.5)).OR.
     &     ((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.3.and.lb1.le.5)))then
c        first classify the reactions due to total charge.
           nchrg=-100
           if((lb1.eq.15.and.(lb2.eq.3.or.lb2.eq.25)).OR.
     &              (lb2.eq.15.and.(lb1.eq.3.or.lb1.eq.25))) then
                nchrg=-2
c     D-
                  bmass=1.232
           endif
           if((lb1.eq.15.and.lb2mn0).or.(lb2.eq.15.and.lb1mn0).OR.
     &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.3.or.lb2.eq.25)).OR.
     &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.3.or.lb1.eq.25)))then
                nchrg=-1
c     n
                 bmass=0.938
           endif
           if((lb1.eq.15.and.(lb2.eq.5.or.lb2.eq.27)).OR.
     &              (lb2.eq.15.and.(lb1.eq.5.or.lb1.eq.27)).or.
     &        (lb1.eq.17.and.(lb2.eq.3.or.lb2.eq.25)).OR.
     &              (lb2.eq.17.and.(lb1.eq.3.or.lb1.eq.25)).or.
     &       ((lb1.eq.14.or.lb1.eq.16).and.lb2mn0).OR.
     &       ((lb2.eq.14.or.lb2.eq.16).and.lb1mn0)) then
                nchrg=0
c     p
                 bmass=0.938
           endif
           if((lb1.eq.17.and.lb2mn0).or.(lb2.eq.17.and.lb1mn0).OR.
     &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.5.or.lb2.eq.27)).OR.
     &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.5.or.lb1.eq.27)))then
                nchrg=1
c     D++
                 bmass=1.232
           endif
           sig = 0.
           if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then
c     PI+sigma or PI + Lambda => Kbar + N reactions
             icase=4
c             pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
             pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
c     lambda + Pi
             if(lb1.eq.14.or.lb2.eq.14) then
                if(nchrg.ge.0) sigma0=akPlam(pkaon)
                if(nchrg.lt.0) sigma0=akNlam(pkaon)
c     sigma + pi
             else
c     K-p or K-D++
                if(nchrg.ge.0) sigma0=akPsgm(pkaon)
c     K-n or K-D-
                if(nchrg.lt.0) sigma0=akNsgm(pkaon)

cbz3/8/99 neutralk
                SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
cbz3/8/99 neutralk end

             endif
             sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
     &         (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
cbz3/8/99 neutralk
c     if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D-
c     K0barD++, K-D-
             if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig

cbz3/8/99 neutralk end

c             the factor 2 comes from spin of delta, which is 3/2
c             detailed balance. copy from Page 423 of N.P. A614 1997

cbz3/8/99 neutralk
             IF (LB1 .EQ. 14 .OR. LB2 .EQ. 14) THEN
                SIG = 4.0 / 3.0 * SIG
             ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN
                SIG = 8.0 / 9.0 * SIG
             ELSE
                SIG = 4.0 / 9.0 * SIG
             END IF
cbz3/8/99 neutralk end
             brsig = sig
             if(sig.lt.1.e-7) sig = 1.e-7
           endif
csp05/07/01
* comment icase=4 statement below if only inelastic
c     PI+L/Si => Kbar + N  OR ELASTIC SCATTERING
           icase=4
           brsig = sig
c     elastic xsecn of 10mb
           sigela = 10.
           sig = sig + sigela
           brel = sigela/sig
cc          brsig = sig
csp05/07/01 end   
        endif
c
c        if(em2.lt.0.2.and.em1.lt.0.2) then
c     PI + PI 
c             icase=5
c     assumed PI PI total x section.
c              sig=50.
c     Mk + Mkbar
c              s0=aka+aka
c              brsig = 0.
c              if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76
c              x section for PIPI->KKbar   PRC43 (1991) 1881
c        endif
        if(icase.eq.-1) then
           ictrl = -1
           return
        endif
        px1cm=pcx
        py1cm=pcy
        pz1cm=pcz
        ds=sqrt(sig/31.4)
        dsr=ds+0.1
        ec=(em1+em2+0.02)**2
c        ec=3.59709
c        if((e(i1).ge.1.).and.(e(i2).ge.1.)) ec = 4.75

        call distce(i1,i2,dsr,ds,dt,ec,srt,ic,px1cm,py1cm,pz1cm)
        if(ic.eq.-1) then
c     no anti-kaon production
           ictrl = -1
c           in=in+1
c           write(60,*)'--------------distance-----',in
           return
        endif

clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi, 
c     sgsum,sgsum1,sgsum3:
        ik=0
        ik0=0
        ik1=0
        ik2=0
        ik3=0
        il=0
        im=0
        im3=0
        im4=0
        in=0
        inpion=0
        ipipi=0
        sgsum=0.
        sgsum1=0.
        sgsum3=0.
        if(icase.eq.1) then
           ik=ik+1
           if(srt.gt.2.8639) then
                ik0=ik0+1
                if(em1.lt.1.0.and.em2.lt.1.0) then
                        ik1=ik1+1
                        sgsum1=sgsum1+brsig
c                        ratio_1=sgsum1/ik1/40.
                endif
                if(em1.gt.1.0.and.em2.gt.1.0) then
                        ik3=ik3+1
                        sgsum3=sgsum3+brsig
c                        ratio_3=sgsum3/ik3/40.
                endif
                if(em1.gt.1.0.and.em2.lt.1.0) ik2=ik2+1
                if(em1.lt.1.0.and.em2.gt.1.0) ik2=ik2+1
                sgsum=sgsum+brsig
c                ratio=sgsum/ik0/40.
           endif
        endif
        if(icase.eq.2) inpion=inpion+1
        if(icase.eq.5) ipipi=ipipi+1
c        write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio
c        write(62,*)'inpion,ipipi',inpion,ipipi
        if(RANART(NSEED).gt.(brsig/sig)) then
c     no anti-kaon production
           ictrl = -1
           return
        endif
        il=il+1
c        kaons could be created now.
        if(icase.eq.1) then
          in=in+1
c          write(60,*)'------in,s2kaon,sig=',in,brsig,sig,lb1,lb2
          call nnkaon(irun,iseed,
     &          ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
        endif
        if(icase.eq.2) then
          im=im+1
c          call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
c     &              pcx,pcy,pcz,nchrg,ratiok)
          call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
     &              pcx,pcy,pcz,nchrg,ratiok,iblock)
        endif
c
        if(icase.eq.3) then
          im3=im3+1
c          write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon
c          write(63,*)'sig,el,sigma',sig,brel,brsgm
c          write(63,*)'srt,pcx,pcy,pcz,em1,em2',srt,pcx,pcy,pcz,em1,em2
          call kaonN(brel,brsgm,irun,iseed,dt,nt,ictrl,
     &                i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
c         this subroutine format is diff. since three final states are possible
        endif
c

        if(icase.eq.4) then
          im4=im4+1
c          write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig
c          write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon

csp06/07/01
      if(RANART(NSEED).lt.brel) then
         ielstc = 1
      else
         ielstc = 0
      endif                  
c          call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
c     &                   pcx,pcy,pcz,nchrg)
          call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
     &                   pcx,pcy,pcz,nchrg,iblock)

csp06/07/01 end
        endif
c        if(icase.eq.5) then
c          im5=im5+1
c          write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig
c          call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz)
c        endif
cbz3/2/99
c        write(101,*)lb1,lb2,lb(i1),lb(i2)
c        write(101,*)em1,em2,e(i1),e(i2),srt
cbz3/2/99end

        return
        end

******************************************
* for pp-->pp + kaon + anti-kaon
c      real*4 function X2kaon(srt)
      real function X2kaon(srt)
      SAVE   
*  This function contains the experimental total pp->pp+K(+)K(-) Xsections    *
*  srt    = DSQRT(s) in GeV                                                   *
*  xsec   = production cross section in mb                                    *
*                                                                             *
******************************************
c     minimum c.m.s. energy to create 2 kaon. = 2*(mp+mk)        
        smin = 2.8639
        x2kaon=0.0000001
        if(srt.lt.smin)return
        sigma1 = 2.8
        sigma2 = 7.7
        sigma3 = 3.9
        x = srt**2/smin**2 + 0.0000001
        f1 = (1.+1./sqrt(x))*alog(x) - 4.*(1.-1./sqrt(x))
        f2 = 1. - (1./sqrt(x))*(1.+alog(sqrt(x)))
        f3 = ((x-1.)/x**2)**3.5
        x2kaon = (1.-1./x)**3*(sigma1*f1 + sigma2*f2) + sigma3*f3
        return
        END

        real function piNsg0(srt)
      SAVE   
* cross section in mb for PI- + P -> P + K0 + K-
c     Mn + 2* Mk
        srt0 = 0.938 + 2.*0.498
        if(srt.lt.srt0) then
           piNsg0 = 0.0
           return
        endif
        ratio = srt0**2/srt**2
        piNsg0=1.121*(1.-ratio)**1.86*ratio**2
        return
        end

        real function akNel(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c        the following data come from PRC 41 (1701)
c        sigma1: K(-) + neutron elastic
        if(pkaon.lt.0.5.or. pkaon.ge.4.0) sigma1=0.
        if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma1=20.*pkaon**2.74
        if(pkaon.ge.1.0.and.pkaon.lt.4.0) sigma1=20.*pkaon**(-1.8)
        akNel=sigma1
        return
        end

        real function akPel(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c        the following data come from PRC 41 (1701)
c        sigma2: K(-) + proton elastic
        if(pkaon.lt.0.25.or. pkaon.ge.4.0) sigma2=0.
        if(pkaon.ge.0.25.and.pkaon.lt.4.0) sigma2=13.*pkaon**(-0.9)
        akPel=sigma2
        return
        end

        real function akNsgm(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c        sigma2: x section for K- + n -> sigma0 + PI-
        if(pkaon.lt.0.5.or. pkaon.ge.6.0) sigma2=0.
        if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma2=1.2*pkaon**(-1.3)
        if(pkaon.ge.1.0.and.pkaon.lt.6.0) sigma2=1.2*pkaon**(-2.3)
        akNsgm=sigma2
        return
        end

        real function akPsgm(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c        sigma1: x section for K- + p -> sigma0 + PI0
        if(pkaon.lt.0.2.or. pkaon.ge.1.5) sigma1=0.
        if(pkaon.ge.0.2.and.pkaon.lt.1.5) sigma1=0.6*pkaon**(-1.8)
        akPsgm=sigma1
        return
        end

        real function akPlam(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c        sigma: x section for K- + p -> lambda + PI0
        p=pkaon
        if(pkaon.lt.0.2.or. pkaon.ge.10.0) sigma=0.
        if(pkaon.ge.0.2.and.pkaon.lt.0.9) sigma=50.*p**2-67.*p+24.
        if(pkaon.ge.0.9.and.pkaon.lt.10.0) sigma=3.0*pkaon**(-2.6)
        akPlam=sigma
        return
        end

        real function akNlam(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
        akNlam=akPlam(pkaon)
        return
        end

* GQ Li parametrization (without resonance)
        real function akNPsg(pkaon)
      SAVE   
*cross section in mb for K- + N reactions.
c       sigma1: x section for K- + p/n -> sigma0 + PI0
         if(pkaon.le.0.345)then
           sigma1=0.624*pkaon**(-1.83)
         else
           sigma1=0.7*pkaon**(-2.09)
         endif
        akNPsg=sigma1
        return
        end   

c-----------------------------------------------------------------------

c.....extracted from G. Song's ART expasion including K- interactions
c.....file `NEWNNK.FOR'

        subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock,
     &                                   srt,pcx,pcy,pcz,nchrg)
c        <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997
c     CANCELED also alpha=1 changed to alpha=3 to decrease the leadng effect.
      PARAMETER      (MAXSTR=150001,MAXR=1)
      PARAMETER      (AKA=0.498)
      COMMON   /AA/  R(3,MAXSTR)
cc      SAVE /AA/
      COMMON   /BB/  P(3,MAXSTR)
cc      SAVE /BB/
      COMMON   /CC/  E(MAXSTR)
cc      SAVE /CC/
      COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
cc      SAVE /EE/
      COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
cc      SAVE /BG/
      COMMON   /NN/NNN
cc      SAVE /NN/
      COMMON   /RUN/NUM
cc      SAVE /RUN/
      COMMON   /PA/RPION(3,MAXSTR,MAXR)
cc      SAVE /PA/
      COMMON   /PB/PPION(3,MAXSTR,MAXR)
cc      SAVE /PB/
      COMMON   /PC/EPION(MAXSTR,MAXR)
cc      SAVE /PC/
      COMMON   /PD/LPION(MAXSTR,MAXR)
cc      SAVE /PD/
      dimension px(4),py(4),pz(4)
      COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
     1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
     2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
      SAVE   
c      dm1=e(i1)
c      dm2=e(i2)
      dm3=0.938
      dm4=0.938
c     10/24/02 initialize n to 0:
      n=0

cbz3/11/99 neutralk
c        if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232
c        if(nchrg.eq.4) dm4=1.232
        if(nchrg.le.-1.or.nchrg.ge.3) dm3=1.232
        if(nchrg.eq.-2.or.nchrg.eq.4) dm4=1.232
cbz3/11/99 neutralk end
          iblock = 0 
        call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
        if(iflag.lt.0) then
c           write(60,*)'------------final state fail-------',n
c     no anti-kaon production
           ictrl = -1
           n=n+1
           return
        endif
        iblock = 12
* Rotate the momenta of particles in the cms of I1 & I2
* px(1), py(1), pz(1): momentum of I1
* px(2), py(2), pz(2): momentum of I2
* px(3), py(3), pz(3): momentum of anti-kaon
* px(4), py(4), pz(4): momentum of kaon


c     10/28/02 get rid of argument usage mismatch in rotate():
        pxrota=px(1)
        pyrota=py(1)
        pzrota=pz(1)
c        call rotate(pcx,pcy,pcz,px(1),py(1),pz(1))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        px(1)=pxrota
        py(1)=pyrota
        pz(1)=pzrota
c
        pxrota=px(2)
        pyrota=py(2)
        pzrota=pz(2)
c        call rotate(pcx,pcy,pcz,px(2),py(2),pz(2))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        px(2)=pxrota
        py(2)=pyrota
        pz(2)=pzrota
c
        pxrota=px(3)
        pyrota=py(3)
        pzrota=pz(3)
c        call rotate(pcx,pcy,pcz,px(3),py(3),pz(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        px(3)=pxrota
        py(3)=pyrota
        pz(3)=pzrota
c
        pxrota=px(4)
        pyrota=py(4)
        pzrota=pz(4)
c        call rotate(pcx,pcy,pcz,px(4),py(4),pz(4))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        px(4)=pxrota
        py(4)=pyrota
        pz(4)=pzrota

        nnn=nnn+2
c     K+
        lpion(nnn,irun)=23
        if(nchrg.eq.-1.or.nchrg.eq.-2) then
c        To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K-

cbz3/7/99 neutralk
c           lpion(nnn,irun)=24 ! K0
cbz3/7/99 neutralk end

        endif
c     aka: rest mass of K
        epion(nnn,irun)=aka
c     K-
        lpion(nnn-1,irun)=21
c     aka: rest mass of K
        epion(nnn-1,irun)=aka
* Find the momenta of particles in the final state in the nucleus_nucleus
* cms frame.   Lorentz transformation into lab frame.
        e1cm   = sqrt(dm3**2 + px(1)**2 + py(1)**2 + pz(1)**2)
        p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz
        transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
        pt1i1 = betax*transf + px(1)
        pt2i1 = betay*transf + py(1)
        pt3i1 = betaz*transf + pz(1)
        eti1  = dm3
c        lb1   = lb(i1)
        lb1   = 2
        if(nchrg.ge.-2.and.nchrg.le.1) lb1=2

cbz3/7/99 neutralk
        if (nchrg .eq. -2 .or. nchrg .eq. -1) then
           lb1 = 6
        end if
cbz3/7/99 neutralk end

cbz3/11/99 neutralk
c        if(nchrg.eq.2.or.nchrg.eq.3) lb1=1
c        if(nchrg.eq.4) lb1=9
        if(nchrg.eq.1.or.nchrg.eq.2) lb1=1
        if(nchrg.eq.3.or.nchrg.eq.4) lb1=9
cbz3/11/99 neutralk end

* For second nulceon, same
        e2cm   = sqrt(dm4**2 + px(2)**2 + py(2)**2 + pz(2)**2)
        p2beta = px(2)*betax + py(2)*betay + pz(2)*betaz
        transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
        pt1i2 = betax*transf + px(2)
        pt2i2 = betay*transf + py(2)
        pt3i2 = betaz*transf + pz(2)
        eti2  = dm4
c        lb2   = lb(i2)
        lb2   = 2

cbz3/11/99 neutralk
c        if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2
c        if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1
c        if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9
c        if(nchrg.eq.-2) lb2=6
        if(nchrg.ge.-1.or.nchrg.le.1) lb2=2
        if(nchrg.eq. 2.or.nchrg.eq.3) lb2=1
        if(nchrg.eq. 4) lb2=9
        if(nchrg.eq.-2) lb2=6
cbz3/11/99 neutralk end

c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
                p(1,i1)=pt1i1
                p(2,i1)=pt2i1
                p(3,i1)=pt3i1
                e(i1)=eti1
                lb(i1)=lb1
                p(1,i2)=pt1i2
                p(2,i2)=pt2i2
                p(3,i2)=pt3i2
                e(i2)=eti2
                lb(i2)=lb2

c                px1 = p(1,i1)
c                py1 = p(2,i1)
c                pz1 = p(3,i1)
c                em1 = e(i1)
c                id(i1) = 2
c                id(i2) = 2
c                id1 = id(i1)
c                iblock = 101  ! K(+)K(-) production
* Get anti-kaons' momenta and coordinates in nucleus-nucleus cms. frame.
        epcmk = sqrt(epion(nnn-1,irun)**2 + px(3)**2+py(3)**2+pz(3)**2)
        betak = px(3)*betax + py(3)*betay + pz(3)*betaz
        transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
        ppion(1,nnn-1,irun)=betax*transf + px(3)
        ppion(2,nnn-1,irun)=betay*transf + py(3)
        ppion(3,nnn-1,irun)=betaz*transf + pz(3)
        rpion(1,nnn-1,irun)=r(1,i1)
        rpion(2,nnn-1,irun)=r(2,i1)
        rpion(3,nnn-1,irun)=r(3,i1)
clin-5/2008:
        dppion(nnn-1,irun)=dpertp(i1)*dpertp(i2)
* Same thing for kaon **************************************
        epcmak = sqrt(epion(nnn,irun)**2 + px(4)**2 +py(4)**2+pz(4)**2)
        betaak = px(4)*betax + py(4)*betay + pz(4)*betaz
        transf= gamma*(gamma*betaak/(gamma+1.) + epcmak)
        ppion(1,nnn,irun)=betax*transf + px(4)
        ppion(2,nnn,irun)=betay*transf + py(4)
        ppion(3,nnn,irun)=betaz*transf + pz(4)
        rpion(1,nnn,irun)=r(1,i2)
        rpion(2,nnn,irun)=r(2,i2)
        rpion(3,nnn,irun)=r(3,i2)
clin-5/2008:
        dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
        return
        end

        subroutine lorntz(ilo,b,pi,pj)
c       It uses to perform Lorentz (or inverse Lorentz) transformation
        dimension pi(4),pj(4),b(3)
      SAVE   
c       dimension db(3)
        bb=b(1)*b(1)+b(2)*b(2)+b(3)*b(3)
        deno3=sqrt(1.-bb)
        if(deno3.eq.0.)deno3=1.e-10
        gam=1./deno3
        ga=gam*gam/(gam+1.)
        if(ilo.eq.1) goto 100
c       Lorentz transformation
        pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
        pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
c       drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3)
c       drdb=db(1)*b(1)+db(2)*b(2)+db(3)*b(3)
        do 1001 i=1,3
           pi(i)=pi(i)+b(i)*(ga*pib-gam*pi(4))
           pj(i)=pj(i)+b(i)*(ga*pjb-gam*pj(4))
c       drd(i)=drd(i)+b(i)*ga*drb
c       db(i)=db(i)+b(i)*ga*drdb
 1001   continue
        pi(4)=gam*(pi(4)-pib)
        pj(4)=gam*(pj(4)-pjb)
        return
100     continue
c       inverse Lorentz transformation
        pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
        pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
        do 1002 i=1,3
           pi(i)=pi(i)+b(i)*(ga*pib+gam*pi(4))
           pj(i)=pj(i)+b(i)*(ga*pjb+gam*pj(4))
 1002   continue
        pi(4)=gam*(pi(4)+pib)
        pj(4)=gam*(pj(4)+pjb)
        return
        end
        
        subroutine fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
*        function: decide final momentum for N,N,K(+),and K(-)        
        dimension px(4), py(4), pz(4), pe(4)
        COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
        SAVE   

        iflag=-1
c        iflag=-1: fail to find momenta
c             = 1: success
        pio=3.1415926
        aka=0.498
c        v=0.43
c        w=-0.84
c        b=3.78
c        c=0.47
c        d=3.60
c        fmax=1.056
c        gmax=1.+c

        icount=0
        ekmax=(srt-dm3-dm4)/2.
        if(ekmax.le.aka) return
        pkmax=sqrt(ekmax**2-aka**2)

        if(dm3.le.0.0.or.dm4.le.0.0) then
           write(1,*)'error: minus mass!!!'
           return
        endif

c        after we have the momenta for both nucleus, we sample the
c        transverse momentum for K-. 
c        dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on
c        page 72, fig 23i.
50        continue
        icount=icount+1
        if(icount.gt.10) return
        ptkmi2=-1./4.145*alog(RANART(NSEED))
        ptkm=sqrt(ptkmi2)
3        v1=RANART(NSEED)
        v2=RANART(NSEED)
        rsq=v1**2+v2**2
        if(rsq.ge.1.0.or.rsq.le.0.) goto 3
        fac=sqrt(-2.*alog(rsq)/rsq)
        guass=v1*fac
        if(guass.ge.5.) goto 3
        xstar=guass/5.
        pzkm=pkmax*xstar
        ekm=sqrt(aka**2+pzkm**2+ptkm**2)
        if(RANART(NSEED).gt.aka/ekm) goto 50
        bbb=RANART(NSEED)
        px(3)=ptkm*cos(2.*pio*bbb)
        py(3)=ptkm*sin(2.*pio*bbb)
        if(RANART(NSEED).gt.0.5) pzkm=-1.*pzkm
        pz(3)=pzkm
        pe(3)=ekm
150        ptkpl2=-1./3.68*alog(RANART(NSEED))
        ptkp=sqrt(ptkpl2)
13        v1=RANART(NSEED)
        v2=RANART(NSEED)
        rsq=v1**2+v2**2
        if(rsq.ge.1.0.or.rsq.le.0.) goto 13
        fac=sqrt(-2.*alog(rsq)/rsq)
        guass=v1*fac
        if(guass.ge.3.25) goto 13
        xstar=guass/3.25
        pzkp=pkmax*xstar
        ekp=sqrt(aka**2+pzkp**2+ptkp**2)
        if(RANART(NSEED).gt.aka/ekp) goto 150
        bbb=RANART(NSEED)
        px(4)=ptkp*cos(2.*pio*bbb)
        py(4)=ptkp*sin(2.*pio*bbb)
        if(RANART(NSEED).gt.0.5) pzkp=-1.*pzkp
        pz(4)=pzkp
        pe(4)=ekp

        resten=srt-pe(3)-pe(4)
        restpz=-pz(3)-pz(4)
c     resample
        if(resten.le.abs(restpz)) goto 50
        restms=sqrt(resten**2-restpz**2)
c     resample 
        if(restms.lt.(dm3+dm4)) goto 50
        ptp2=-1./2.76*alog(RANART(NSEED))
        ptp=sqrt(ptp2)
        bbb=RANART(NSEED)
        px(2)=ptp*cos(2.*pio*bbb)
        py(2)=ptp*sin(2.*pio*bbb)
        px(1)=-1.*(px(4)+px(3)+px(2))
        py(1)=-1.*(py(4)+py(3)+py(2))
c     transverse mass for K-
        rmt3=sqrt(dm3**2+px(1)**2+py(1)**2)
c     transverse mass for K+
        rmt4=sqrt(dm4**2+px(2)**2+py(2)**2)
        if(restms.lt.(rmt3+rmt4)) goto 50
c        else: sampling success!
        pzcms=sqrt((restms**2-(rmt3+rmt4)**2)*
     &             (restms**2-(rmt3-rmt4)**2))/2./restms
        if(RANART(NSEED).gt.0.5) then
           pz(1)=pzcms
           pz(2)=-pzcms
        else
           pz(1)=-pzcms
           pz(2)=pzcms
        endif
        beta=restpz/resten        
        gama=1./sqrt(1.-beta**2)
        pz(1)=pz(1)*gama + beta*gama*sqrt(rmt3**2+pz(1)**2)
        pz(2)=pz(2)*gama + beta*gama*sqrt(rmt4**2+pz(2)**2)
        pe(1)=sqrt(rmt3**2+pz(1)**2)
        pe(2)=sqrt(rmt4**2+pz(2)**2)

        iflag=1
        return
        end

c-----------------------------------------------------------------------

c.....extracted from G. Song's ART expasion including K- interactions
c.....file `NPIK.FOR'

****************************************
c        subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
c     &                  pcx,pcy,pcz,nchrg,ratiok)
        subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
     &                  pcx,pcy,pcz,nchrg,ratiok,iblock)
*
* Process: PI + N -> K(-) + ANYTHING
* 1.  PI- + P -> P + K0 + K-
* 2.  PI+ + N -> P + K+ + K- 
* 3.  PI0 + P -> P + K+ + K-
* 4.  PI0 + N -> P + K0 + K-
* 5.  PI0 + N -> N + K+ + K-
* 6.  PI- + P -> N + K+ + K-
* 7.  PI- + N -> N + K0 + K-
* NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
****************************************
      PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
      PARAMETER      (AKA=0.498)
      COMMON   /AA/  R(3,MAXSTR)
cc      SAVE /AA/
      COMMON   /BB/  P(3,MAXSTR)
cc      SAVE /BB/
      COMMON   /CC/  E(MAXSTR)
cc      SAVE /CC/
      COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
cc      SAVE /EE/
      COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
cc      SAVE /BG/
      COMMON   /NN/NNN
cc      SAVE /NN/
      COMMON   /RUN/NUM
cc      SAVE /RUN/
      COMMON   /PA/RPION(3,MAXSTR,MAXR)
cc      SAVE /PA/
      COMMON   /PB/PPION(3,MAXSTR,MAXR)
cc      SAVE /PB/
      COMMON   /PC/EPION(MAXSTR,MAXR)
cc      SAVE /PC/
      COMMON   /PD/LPION(MAXSTR,MAXR)
cc      SAVE /PD/
      dimension bb(3),p1(4),p2(4),p3(4)
cms      ,px(4),py(4),pz(4)
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
     1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
     2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
      SAVE   
        px1cm=pcx
        py1cm=pcy
        pz1cm=pcz
        ictrl = 1
        lb1=lb(i1)
        lb2=lb(i2)
        k1=i1
        k2=i2
c        k1 must be bayron. k2 be meson. If not, exchange.
        if(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13)) then
            k1=i2
            k2=i1
        endif
cbz3/8/99 neutralk
cbz10/12/99
c        LB(I1) = 1 + 2 * RANART(NSEED)
c        LB(I2) = 23
        LB(k1) = 1 + int(2*RANART(NSEED))
        LB(k2) = 23
c       pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2))
c     &           /2./srt
        pkmax=sqrt((srt**2-(aka+0.938+aka)**2)
     &           *(srt**2-(aka+0.938-aka)**2))/2./srt
        pk = RANART(NSEED)*pkmax
c-----------------------------------------------------
        css=1.-2.*RANART(NSEED)
        sss=sqrt(1.-css**2)
        fai=2*3.1415926*RANART(NSEED)
        p3(1)=pk*sss*cos(fai)
        p3(2)=pk*sss*sin(fai)
        p3(3)=pk*css
        eip = srt - sqrt(aka**2 + pk**2)
        rmnp=sqrt(eip**2-pk**2)
        do 1001 i= 1, 3
           bb(i) = -1.*p3(i)/eip
 1001   continue
c        bb: velocity of the other two particles as a whole.
        pznp=sqrt((rmnp**2-(aka+0.938)**2)
     c  *(rmnp**2-(0.938-aka)**2))/2./rmnp    
c-----------------------------------------------------
        css=1.-2.*RANART(NSEED)
        sss=sqrt(1.-css**2)
        fai=2*3.1415926*RANART(NSEED)
        p1(1)=pznp*sss*cos(fai)
        p1(2)=pznp*sss*sin(fai)
        p1(3)=pznp*css
        p1(4)=sqrt(0.938**2+pznp**2)
        p2(4)=sqrt(aka**2+pznp**2)
        do 1002 i=1,3
           p2(i)=-1.*p1(i)
 1002   continue
c        p1,p2: the momenta of the two particles in their cms
c        p1: momentum of N or P
c        p2: momentum of anti_kaon
c        p3: momentum of K0 or K+
        ilo=1
c        write(61,*)'--------p1,p2',p1,p2
c        write(61,*)'--------bb',bb
        call lorntz(ilo,bb,p1,p2)
c******* Checking *************
c        pxsum = p1(1)+p2(1)+p3(1)
c        pysum = p1(2)+p2(2)+p3(2)
c        pzsum = p1(3)+p2(3)+p3(3)
c        pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt
c        write(61,*)'---p1,pxsum',p1,pxsum
c        write(61,*)'---p2,pysum',p2,pysum
c        write(61,*)'---p3,pzsum',p3,pzsum
c        write(61,*)'---pesum',pesum
c***********************************

* Rotate the momenta of particles in the cms of I1 & I2
* px(1), py(1), pz(1): momentum of I1
* px(2), py(2), pz(2): momentum of I2
* px(3), py(3), pz(3): momentum of anti-kaon

c     10/28/02 get rid of argument usage mismatch in rotate():
        pxrota=p1(1)
        pyrota=p1(2)
        pzrota=p1(3)
c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p1(1)=pxrota
        p1(2)=pyrota
        p1(3)=pzrota
c
        pxrota=p2(1)
        pyrota=p2(2)
        pzrota=p2(3)
c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p2(1)=pxrota
        p2(2)=pyrota
        p2(3)=pzrota
c
        pxrota=p3(1)
        pyrota=p3(2)
        pzrota=p3(3)
c        call rotate(pcx,pcy,pcz,p3(1),p3(2),p3(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p3(1)=pxrota
        p3(2)=pyrota
        p3(3)=pzrota

        nnn=nnn+1
c     K(-)
        lpion(nnn,irun)=21
c     aka: rest mass of K
        epion(nnn,irun)=aka
* Find the momenta of particles in the final state in the nucleus_nucleus
* cms frame.   Lorentz transformation into lab frame.
        e1cm   = sqrt(0.938**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
        p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
        transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
        pt1i1 = betax*transf + p1(1)
        pt2i1 = betay*transf + p1(2)
        pt3i1 = betaz*transf + p1(3)
        eti1  = 0.938
        lb1   = lb(k1)
         
* For second nulceon, same
        e2cm   = sqrt(aka**2 + p3(1)**2 + p3(2)**2 + p3(3)**2)
        p2beta = p3(1)*betax + p3(2)*betay + p3(3)*betaz
        transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
        pt1i2 = betax*transf + p3(1)
        pt2i2 = betay*transf + p3(2)
        pt3i2 = betaz*transf + p3(3)
        eti2  = aka
        lb2   = lb(k2)

c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
*       k1 stand for nucleon, k2 stand for kaon. lpion stand for Kbar.
                p(1,k1)=pt1i1
                p(2,k1)=pt2i1
                p(3,k1)=pt3i1
                e(k1)=eti1
                lb(k1)=lb1
                p(1,k2)=pt1i2
                p(2,k2)=pt2i2
                p(3,k2)=pt3i2
                e(k2)=eti2
                lb(k2)=lb2

c                px1 = p(1,i1)
c                py1 = p(2,i1)
c                pz1 = p(3,i1)
c                em1 = e(i1)
c                id(i1) = 2
c                id(i2) = 2
c                id1 = id(i1)
c     K(+)K(-) production
                iblock = 101
* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
c  p2:  momentum of anti-kaon.
c        epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
        epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2+p2(2)**2+p2(3)**2)
        betak = p2(1)*betax + p2(2)*betay + p2(3)*betaz
        transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
        ppion(1,nnn,irun)=betax*transf + p2(1)
        ppion(2,nnn,irun)=betay*transf + p2(2)
        ppion(3,nnn,irun)=betaz*transf + p2(3)
clin-5/2008:
        dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
cbz3/2/99
c        write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun),
c     &                    ppion(3,nnn,irun), dt*nt, srt
cbz3/2/99end
c        write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun),
c     &                    ppion(3,nnn,irun), dt*nt, srt
        k=i2
        if(lb(i1).eq.1.or.lb(i1).eq.2) k=i1
        rpion(1,nnn,irun)=r(1,k)
        rpion(2,nnn,irun)=r(2,k)
        rpion(3,nnn,irun)=r(3,k)
        return
        end

c-----------------------------------------------------------------------

c.....extracted from G. Song's ART expasion including K- interactions
c.....file `PIHYPN.FOR'

******************************************
        subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,
     &     srt,pcx,pcy,pcz,nchrg,iblock)
*
* Process: PI + sigma(or Lambda) -> Kbar + N
* NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
******************************************

c NOTE: for PI + Hyperon: the produced kaons have mass 0.498
      PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
      PARAMETER      (AKA=0.498)
      COMMON   /AA/  R(3,MAXSTR)
cc      SAVE /AA/
      COMMON   /BB/  P(3,MAXSTR)
cc      SAVE /BB/
      COMMON   /CC/  E(MAXSTR)
cc      SAVE /CC/
      COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
cc      SAVE /EE/
      COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
cc      SAVE /BG/
      COMMON   /NN/NNN
cc      SAVE /NN/
      COMMON   /RUN/NUM
cc      SAVE /RUN/
      COMMON   /PA/RPION(3,MAXSTR,MAXR)
cc      SAVE /PA/
      COMMON   /PB/PPION(3,MAXSTR,MAXR)
cc      SAVE /PB/
      COMMON   /PC/EPION(MAXSTR,MAXR)
cc      SAVE /PC/
      COMMON   /PD/LPION(MAXSTR,MAXR)
cc      SAVE /PD/
      dimension p1(4),p2(4)
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      SAVE   
        px1cm=pcx
        py1cm=pcy
        pz1cm=pcz
        ictrl = 1
csp06/07/01
        if(ielstc .eq. 1) then
*    L/Si + meson -> L/Si + meson
             k1=i1
             k2=i2
           dm3=e(k1)
           dm4=e(k2)
           iblock = 10
        else
           iblock = 12
csp06/07/01 end  
c        PI + Sigma(or Lambda) -> Kbar + N
        k1=i1
        k2=i2
c        k1 must be bayron! So if I1 is PI, exchange k1 & k2.
        if(lb(i1).lt.14.or.lb(i1).gt.17) then
           k1=i2
           k2=i1
        endif
cbz3/8/99 neutralk
        LB(K1) = 1 + int(2*RANART(NSEED))
        if(nchrg.eq.-2) lb(k1)=6
c     if(nchrg.eq.-1) lb(k1)=2
c     if(nchrg.eq. 0) lb(k1)=1
c     if(nchrg.eq. 1) lb(k1)=9
        IF (NCHRG .EQ. 2) LB(K1) = 9
cbz3/8/99 neutralk end

c     K-
        lb(k2)=21
        dm3=0.938
        if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232
        dm4=aka
c        dm3,dm4: the mass of final state particles.
         endif
    
********Now, antikaon will be created.
c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
c        pkmax: the maximum momentum of anti-kaon
        pkmax=sqrt((srt**2-(dm3+dm4)**2)*(srt**2-(dm3-dm4)**2))
     &         /2./srt
        pk=pkmax
c-----------------------------------------------------
        css=1.-2.*RANART(NSEED)
        sss=sqrt(1.-css**2)
        fai=2*3.1415926*RANART(NSEED)
        p1(1)=pk*sss*cos(fai)
        p1(2)=pk*sss*sin(fai)
        p1(3)=pk*css
        do 1001 i=1,3
           p2(i)=-1.*p1(i)
 1001   continue
c        p1,p2: the momenta of the two particles in their cms
c        p1: momentum of kaon
c        p2: momentum of Kbar

* Rotate the momenta of particles in the cms of I1 & I2
clin-10/28/02 get rid of argument usage mismatch in rotate():
        pxrota=p1(1)
        pyrota=p1(2)
        pzrota=p1(3)
c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p1(1)=pxrota
        p1(2)=pyrota
        p1(3)=pzrota
c
        pxrota=p2(1)
        pyrota=p2(2)
        pzrota=p2(3)
c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p2(1)=pxrota
        p2(2)=pyrota
        p2(3)=pzrota
clin-10/28/02-end

* Find the momenta of particles in the final state in the nucleus_nucleus
* cms frame.   Lorentz transformation into lab frame.
        e1cm   = sqrt(dm3**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
        p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
        transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
        pt1i1 = betax*transf + p1(1)
        pt2i1 = betay*transf + p1(2)
        pt3i1 = betaz*transf + p1(3)
        eti1  = dm3
        lb1   = lb(k1)
         
* For second kaon, same
        e2cm   = sqrt(dm4**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
        p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
        transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
        pt1i2 = betax*transf + p2(1)
        pt2i2 = betay*transf + p2(2)
        pt3i2 = betaz*transf + p2(3)
cbz3/2/99
c        write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt
cbz3/2/99end
c        write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt
        eti2  = dm4
        lb2   = lb(k2)

c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
c        k1=i1
c        k2=i2
*       k1 stand for nucleon, k2 stand for kaon.
                p(1,k1)=pt1i1
                p(2,k1)=pt2i1
                p(3,k1)=pt3i1
                e(k1)=eti1
                lb(k1)=lb1
                p(1,k2)=pt1i2
                p(2,k2)=pt2i2
                p(3,k2)=pt3i2
                e(k2)=eti2
                lb(k2)=lb2

cc                iblock = 101  ! K(+)K(-) production
* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
        return
        end

c-----------------------------------------------------------------------

c.....extracted from G. Song's ART expasion including K- interactions
c.....file `KAONN.FOR'

****************************************
        subroutine kaonN(brel,brsgm,irun,iseed,dt,nt,
     &     ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
*
* Process: PI + sigma(or Lambda) <- Kbar + N
* NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
****************************************
      PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
      PARAMETER      (AKA=0.498,ALA=1.1157,ASA=1.1974)
      COMMON   /AA/  R(3,MAXSTR)
cc      SAVE /AA/
      COMMON   /BB/  P(3,MAXSTR)
cc      SAVE /BB/
      COMMON   /CC/  E(MAXSTR)
cc      SAVE /CC/
      COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
cc      SAVE /EE/
      COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
cc      SAVE /BG/
      COMMON   /NN/NNN
cc      SAVE /NN/
      COMMON   /RUN/NUM
cc      SAVE /RUN/
      COMMON   /PA/RPION(3,MAXSTR,MAXR)
cc      SAVE /PA/
      COMMON   /PB/PPION(3,MAXSTR,MAXR)
cc      SAVE /PB/
      COMMON   /PC/EPION(MAXSTR,MAXR)
cc      SAVE /PC/
      COMMON   /PD/LPION(MAXSTR,MAXR)
cc      SAVE /PD/
      dimension p1(4),p2(4)
      COMMON/RNDF77/NSEED
cc      SAVE /RNDF77/
      SAVE   
        px1cm=pcx
        py1cm=pcy
        pz1cm=pcz
        ictrl = 1
c        ratio: used for isospin decision.
        k1=i1
        k2=i2
c        k1 must be bayron! So if I1 is Kaon, exchange k1 & k2.
        if(e(i1).lt.0.5.and.e(i1).gt.0.01) then
           k1=i2
           k2=i1
        endif
*** note: for print out only *******************************
c     record kaon's mass
        eee=e(k2)
*** end **************
        rrr=RANART(NSEED)
        if(rrr.lt.brel) then
c       Kbar + N -> Kbar + N
           lb1=lb(k1)
           lb2=lb(k2)
           em1=e(k1)
           em2=e(k2)
           iblock = 10
        else 
           iblock = 12
        if(rrr.lt.(brel+brsgm)) then
c        nchrg: Net charges of the two incoming particles.
c           Kbar + N -> Sigma + PI
           em1=asa
           em2=0.138

cbz3/8/99 neutralk
           LB1 = 15 + int(3*RANART(NSEED))
           LB2 = 3 + int(3*RANART(NSEED))
        else
c           Kbar + N -> Lambda + PI
           em1=ala
           em2=0.138
c     LAmbda
           lb1=14
cbz3/8/99 neutralk
           LB2 = 3 + int(3*RANART(NSEED))
c           if(nchrg.eq.1)  lb2=5  ! K- + D++ -> Lambda + PI+
c           if(nchrg.eq.0)  lb2=4  ! K- + p(D+,N*+) -> Lambda + PI0
c          if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI-
cbz3/8/99 neutralk

        endif
        endif
        lb(k1)=lb1
        lb(k2)=lb2
    
********Now, antikaon will be created.
c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
c        pkmax: the maximum momentum of anti-kaon
c        write(63,*)'srt,em1,em2',srt,em1,em2
c        write(63,*)'-srt,em1,em2',srt,em1,em2
        pkmax=sqrt((srt**2-(em1+em2)**2)*(srt**2-(em1-em2)**2))
     &         /2./srt
        pk=pkmax
c-----------------------------------------------------
        css=1.-2.*RANART(NSEED)
        sss=sqrt(1.-css**2)
        fai=2*3.1415926*RANART(NSEED)
        p1(1)=pk*sss*cos(fai)
        p1(2)=pk*sss*sin(fai)
        p1(3)=pk*css
        do 1001 i=1,3
           p2(i)=-1.*p1(i)
 1001   continue
c        p1,p2: the momenta of the two particles in their cms
c        p1: momentum of kaon
c        p2: momentum of Kbar

* Rotate the momenta of particles in the cms of I1 & I2

clin-10/28/02 get rid of argument usage mismatch in rotate():
        pxrota=p1(1)
        pyrota=p1(2)
        pzrota=p1(3)
c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p1(1)=pxrota
        p1(2)=pyrota
        p1(3)=pzrota
c
        pxrota=p2(1)
        pyrota=p2(2)
        pzrota=p2(3)
c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
        call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
        p2(1)=pxrota
        p2(2)=pyrota
        p2(3)=pzrota
clin-10/28/02-end

* Find the momenta of particles in the final state in the nucleus_nucleus
* cms frame.   Lorentz transformation into lab frame.
        e1cm   = sqrt(em1**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
        p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
        transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
        pt1i1 = betax*transf + p1(1)
        pt2i1 = betay*transf + p1(2)
        pt3i1 = betaz*transf + p1(3)
        eti1  = em1
         
* For second kaon, same
        e2cm   = sqrt(em2**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
        p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
        transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
        pt1i2 = betax*transf + p2(1)
        pt2i2 = betay*transf + p2(2)
        pt3i2 = betaz*transf + p2(3)
        eti2  = em2

c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
c        k1=i1
c        k2=i2
*       k1 stand for bayron, k2 stand for meson.
                p(1,k1)=pt1i1
                p(2,k1)=pt2i1
                p(3,k1)=pt3i1
                e(k1)=eti1
                p(1,k2)=pt1i2
                p(2,k2)=pt2i2
                p(3,k2)=pt3i2
                e(k2)=eti2

cc                iblock = 101  ! K(+)K(-) production
* Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
        return
        end

c=======================================================================

clin Below is the previous artana.f:
c=======================================================================

c.....analysis subroutine before the hadronic space-time evolution

      SUBROUTINE ARTAN1
      PARAMETER (MAXSTR=150001, MAXR=1)
c.....y cut for mt spectrum
cbz3/17/99
c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
      PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
cbz3/17/99 end
c.....bin width for mt spectrum and y spectrum
c.... clin-9/26/03 no symmetrization in y (or eta) for ana dat files
c      PARAMETER (BMT = 0.05, BY = 0.2)
      PARAMETER (BMT = 0.05, BY = 0.4)
      COMMON /RUN/ NUM
cc      SAVE /RUN/
      COMMON /ARERC1/MULTI1(MAXR)
cc      SAVE /ARERC1/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cbz3/17/99
c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
cc      SAVE /ARPRC1/
      COMMON /ARANA1/
     &     dy1ntb(50), dy1ntp(50), DY1HM(50), 
     &     DY1KP(50), DY1KM(50), DY1K0S(50),
     &     DY1LA(50), DY1LB(50), DY1PHI(50),
     &     dm1pip(50), dm1pim(50), DMT1PR(50),
     &     DMT1PB(50), DMT1KP(50), dm1km(50),
     &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
     &     dy1msn(50), DY1PIP(50), DY1PIM(50), 
     &     DY1PI0(50), DY1PR(50), DY1PB(50)
     &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
cc      SAVE /ARANA1/
      SAVE   

cbz3/17/99 end
      DO 1002 J = 1, NUM
         DO 1001 I = 1, MULTI1(J)
            ITYP = ITYP1(I, J)
            PX = PX1(I, J)
            PY = PY1(I, J)
            PZ = PZ1(I, J)
            EE = EE1(I, J)
            XM = XM1(I, J)
c     2/24/03 leptons and photons:
            if(xm.lt.0.01) goto 200
            ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
clin-9/2012 determine pseudo-rapidity more generally:
c            eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
            if((PX**2+PY**2).gt.0.) then
               eta=asinh(PZ/sqrt(PX**2+PY**2))
            else
               eta = 1000000.0*sign(1.,PZ)
clin-2/2013 for spectator target nucleons in LAB frame,
c     note that finite precision of HBOOST 
c     would give spectator target nucleons a small but non-zero pz:
               if(abs(pz).le.1e-3) eta=0.
            endif

            XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
            DXMT = XMT - XM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. EE) THEN
c               PRINT *, 'IN ARTAN1'
c               PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
ccbzdbg2/16/99
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
ccbzdbg2/16/99
ccbzdbg2/15/99
c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
ccbzdbg2/16/99
c               PRINT *, ' XM = ', XM
ccbzdbg2/16/99end
c               GOTO 200
c            else
cc            Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
c               Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ +1e-5))
cc               STOP
ccbzdbg2/15/99end
c            END IF
            if(XMT.gt.0.) then
               Y=asinh(PZ/XMT)
            else
               PRINT *, ' IN ARTAN1 mt=0'
               Y = 1000000.0*sign(1.,PZ)
            endif

c.....rapidity cut for the rapidity distribution
            IF (ABS(Y) .GE. 10.0) GOTO 100
c.....clin-9/26/03 no symmetrization in y (or eta) for ana dat files
c            IY = 1 + int(ABS(Y) / BY)
c            Ieta = 1 + int(ABS(eta) / BY)
            IF (ABS(eta) .GE. 10.0) GOTO 100
            IY = 1 + int((Y+10.) / BY)
            Ieta = 1 + int((eta+10.) / BY)

            IF (ITYP .LT. -1000) THEN
               dy1ntb(IY) = dy1ntb(IY) - 1.0
            END IF
            IF (ITYP .GT. 1000) THEN
               dy1ntb(IY) = dy1ntb(IY) + 1.0
            END IF
            IF (ITYP .EQ. -2212) THEN
               dy1ntp(IY) = dy1ntp(IY) - 1.0
            END IF
            IF (ITYP .EQ. 2212) THEN
               dy1ntp(IY) = dy1ntp(IY) + 1.0
            END IF
c            IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR.
c     &         ITYP .EQ. -2212) THEN
            IF (ITYP .EQ. -2112) THEN
               DY1HM(IY) = DY1HM(IY) + 1.0
            END IF
c
            IF (LUCHGE(ITYP).ne.0) THEN
               DY1CH(IY) = DY1CH(IY) + 1.0
               DE1CH(Ieta) = DE1CH(Ieta) + 1.0
               IF (LUCHGE(ITYP).lt.0) THEN
                  DY1NEG(IY) = DY1NEG(IY) + 1.0
                  DE1NEG(Ieta) = DE1NEG(Ieta) + 1.0
               endif
            END IF

cbz3/17/99
            IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. 
     &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
               dy1msn(IY) = dy1msn(IY) + 1.0
            END IF
            IF (ITYP .EQ. 211) THEN
               DY1PIP(IY) = DY1PIP(IY) + 1.0
            END IF
            IF (ITYP .EQ. -211) THEN
               DY1PIM(IY) = DY1PIM(IY) + 1.0
            END IF
            IF (ITYP .EQ. 111) THEN
               DY1PI0(IY) = DY1PI0(IY) + 1.0
            END IF
            IF (ITYP .EQ. 2212) THEN
               DY1PR(IY) = DY1PR(IY) + 1.0
            END IF
            IF (ITYP .EQ. -2212) THEN
               DY1PB(IY) = DY1PB(IY) + 1.0
            END IF
cbz3/17/99 end
            IF (ITYP .EQ. 321) THEN
               DY1KP(IY) = DY1KP(IY) + 1.0
            END IF
            IF (ITYP .EQ. -321) THEN
               DY1KM(IY) = DY1KM(IY) + 1.0
            END IF
clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S:
c            IF (ITYP .EQ. 310) THEN
            IF (ITYP .EQ. 130) THEN
               DY1K0S(IY) = DY1K0S(IY) + 1.0
            END IF
            IF (ITYP .EQ. 3122) THEN
               DY1LA(IY) = DY1LA(IY) + 1.0
            END IF
            IF (ITYP .EQ. -3122) THEN
               DY1LB(IY) = DY1LB(IY) + 1.0
            END IF
            IF (ITYP .EQ. 333) THEN
               DY1PHI(IY) = DY1PHI(IY) + 1.0
            END IF

c.....insert rapidity cut for mt spectrum here
 100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
            IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
            IMT = 1 + int(DXMT / BMT)
            IF (ITYP .EQ. 211) THEN
               dm1pip(IMT) = dm1pip(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -211) THEN
               dm1pim(IMT) = dm1pim(IMT) + 
     &            1.0 / XMT
            END IF
            IF (ITYP .EQ. 2212) THEN
               DMT1PR(IMT) = DMT1PR(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -2212) THEN
               DMT1PB(IMT) = DMT1PB(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. 321) THEN
               DMT1KP(IMT) = DMT1KP(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -321) THEN
               dm1km(IMT) = dm1km(IMT) + 1.0 / XMT
            END IF
clin-4/24/03:
c            IF (ITYP .EQ. 310) THEN
            IF (ITYP .EQ. 130) THEN
               dm1k0s(IMT) = dm1k0s(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. 3122) THEN
               DMT1LA(IMT) = DMT1LA(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -3122) THEN
               DMT1LB(IMT) = DMT1LB(IMT) + 1.0 / XMT
            END IF

 200        CONTINUE
 1001    CONTINUE
 1002 CONTINUE

      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine after the hadronic space-time evolution

      SUBROUTINE ARTAN2

      PARAMETER (MAXSTR=150001, MAXR=1)
c.....y cut for mt spectrum
cbz3/17/99
c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
      PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
cbz3/17/99 end
c.....bin width for mt spectrum and y spectrum
c      PARAMETER (BMT = 0.05, BY = 0.2)
      PARAMETER (BMT = 0.05, BY = 0.4)
      COMMON /RUN/ NUM
cc      SAVE /RUN/
      COMMON /ARERC1/MULTI1(MAXR)
cc      SAVE /ARERC1/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cbz3/17/99
c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
cc      SAVE /ARPRC1/
      COMMON /ARANA2/
     &     dy2ntb(50), dy2ntp(50), DY2HM(50), 
     &     DY2KP(50), DY2KM(50), DY2K0S(50),
     &     DY2LA(50), DY2LB(50), DY2PHI(50),
     &     dm2pip(50), dm2pim(50), DMT2PR(50),
     &     DMT2PB(50), DMT2KP(50), dm2km(50),
     &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
     &     dy2msn(50), DY2PIP(50), DY2PIM(50), 
     &     DY2PI0(50), DY2PR(50), DY2PB(50)
     &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
cbz3/17/99 end
cc      SAVE /ARANA2/
      SAVE   

      DO 1002 J = 1, NUM
         DO 1001 I = 1, MULTI1(J)
            ITYP = ITYP1(I, J)
            PX = PX1(I, J)
            PY = PY1(I, J)
            PZ = PZ1(I, J)
            EE = EE1(I, J)
            XM = XM1(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
c     2/24/03 leptons and photons:
            if(xm.lt.0.01) goto 200
            DXMT = XMT - XM
            ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
clin-9/2012 determine pseudo-rapidity more generally:
c            eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
            if((PX**2+PY**2).gt.0.) then
               eta=asinh(PZ/sqrt(PX**2+PY**2))
            else
               eta = 1000000.0*sign(1.,PZ)
clin-2/2013 for spectator target nucleons in LAB frame:
               if(abs(pz).le.1e-3) eta=0.
            endif

clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. EE) THEN
c               PRINT *, 'IN ARTAN2'
c               PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
ccbzdbg2/16/99
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
ccbzdbg2/16/99
ccbzdbg2/15/99
c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
ccbzdbg2/16/99
c               PRINT *, ' XM = ', XM
ccbzdbg2/16/99end
c               GOTO 200
cc               STOP
ccbzdbg2/15/99end
c            END IF
c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
            if(XMT.gt.0.) then
               Y=asinh(PZ/XMT)
            else
               PRINT *, ' IN ARTAN2 mt=0'
               Y = 1000000.0*sign(1.,PZ)
            endif

c.....rapidity cut for the rapidity distribution
            IF (ABS(Y) .GE. 10.0) GOTO 100
c            IY = 1 + int(ABS(Y) / BY)
c            Ieta = 1 + int(ABS(eta) / BY)
            IF (ABS(eta) .GE. 10.0) GOTO 100
            IY = 1 + int((Y+10.) / BY)
            Ieta = 1 + int((eta+10.) / BY)

            IF (ITYP .LT. -1000) THEN
               dy2ntb(IY) = dy2ntb(IY) - 1.0
            END IF
            IF (ITYP .GT. 1000) THEN
               dy2ntb(IY) = dy2ntb(IY) + 1.0
            END IF
            IF (ITYP .EQ. -2212) THEN
               dy2ntp(IY) = dy2ntp(IY) - 1.0
            END IF
            IF (ITYP .EQ. 2212) THEN
               dy2ntp(IY) = dy2ntp(IY) + 1.0
            END IF
            IF (ITYP .EQ. -2112) THEN
               DY2HM(IY) = DY2HM(IY) + 1.0
            END IF

            IF (LUCHGE(ITYP).ne.0) THEN
               DY2CH(IY) = DY2CH(IY) + 1.0
               DE2CH(Ieta) = DE2CH(Ieta) + 1.0
               IF (LUCHGE(ITYP).lt.0) THEN
                  DY2NEG(IY) = DY2NEG(IY) + 1.0
                  DE2NEG(Ieta) = DE2NEG(Ieta) + 1.0
               endif
            END IF

cbz3/17/99
            IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR. 
     &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
               dy2msn(IY) = dy2msn(IY) + 1.0
            END IF
            IF (ITYP .EQ. 211) THEN
               DY2PIP(IY) = DY2PIP(IY) + 1.0
            END IF
            IF (ITYP .EQ. -211) THEN
               DY2PIM(IY) = DY2PIM(IY) + 1.0
            END IF
            IF (ITYP .EQ. 111) THEN
               DY2PI0(IY) = DY2PI0(IY) + 1.0
            END IF
            IF (ITYP .EQ. 2212) THEN
               DY2PR(IY) = DY2PR(IY) + 1.0
            END IF
            IF (ITYP .EQ. -2212) THEN
               DY2PB(IY) = DY2PB(IY) + 1.0
            END IF
cbz3/17/99 end
            IF (ITYP .EQ. 321) THEN
               DY2KP(IY) = DY2KP(IY) + 1.0
            END IF
            IF (ITYP .EQ. -321) THEN
               DY2KM(IY) = DY2KM(IY) + 1.0
            END IF
clin-4/24/03:
c            IF (ITYP .EQ. 310) THEN
            IF (ITYP .EQ. 130) THEN
               DY2K0S(IY) = DY2K0S(IY) + 1.0
            END IF
            IF (ITYP .EQ. 3122) THEN
               DY2LA(IY) = DY2LA(IY) + 1.0
            END IF
            IF (ITYP .EQ. -3122) THEN
               DY2LB(IY) = DY2LB(IY) + 1.0
            END IF
            IF (ITYP .EQ. 333) THEN
               DY2PHI(IY) = DY2PHI(IY) + 1.0
            END IF

c.....insert rapidity cut for mt spectrum here
 100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
            IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
            IMT = 1 + int(DXMT / BMT)
            IF (ITYP .EQ. 211) THEN
               dm2pip(IMT) = dm2pip(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -211) THEN
               dm2pim(IMT) = dm2pim(IMT) + 
     &            1.0 / XMT
            END IF
            IF (ITYP .EQ. 2212) THEN
               DMT2PR(IMT) = DMT2PR(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -2212) THEN
               DMT2PB(IMT) = DMT2PB(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. 321) THEN
               DMT2KP(IMT) = DMT2KP(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -321) THEN
               dm2km(IMT) = dm2km(IMT) + 1.0 / XMT
            END IF
clin-4/24/03:
c            IF (ITYP .EQ. 310) THEN
            IF (ITYP .EQ. 130) THEN
               dm2k0s(IMT) = dm2k0s(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. 3122) THEN
               DMT2LA(IMT) = DMT2LA(IMT) + 1.0 / XMT
            END IF
            IF (ITYP .EQ. -3122) THEN
               DMT2LB(IMT) = DMT2LB(IMT) + 1.0 / XMT
            END IF

 200        CONTINUE
 1001    CONTINUE
 1002 CONTINUE

      RETURN
      END

c-----------------------------------------------------------------------

c.....output analysis results at the end of the simulation

      SUBROUTINE ARTOUT(NEVNT)

      PARAMETER (MAXSTR=150001, MAXR=1)
c.....y cut for mt spectrum
cbz3/17/99
c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
      PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
cbz3/17/99 end
c.....bin width for mt spectrum and y spectrum
c      PARAMETER (BMT = 0.05, BY = 0.2)
      PARAMETER (BMT = 0.05, BY = 0.4)
      COMMON /RUN/ NUM
cc      SAVE /RUN/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cbz3/17/99
c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
cc      SAVE /ARPRC1/
      COMMON /ARANA1/
     &     dy1ntb(50), dy1ntp(50), DY1HM(50), 
     &     DY1KP(50), DY1KM(50), DY1K0S(50),
     &     DY1LA(50), DY1LB(50), DY1PHI(50),
     &     dm1pip(50), dm1pim(50), DMT1PR(50),
     &     DMT1PB(50), DMT1KP(50), dm1km(50),
     &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
     &     dy1msn(50), DY1PIP(50), DY1PIM(50), 
     &     DY1PI0(50), DY1PR(50), DY1PB(50)
     &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
cbz3/17/99 end
cc      SAVE /ARANA1/
cbz3/17/99
c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
      COMMON /ARANA2/
     &     dy2ntb(50), dy2ntp(50), DY2HM(50), 
     &     DY2KP(50), DY2KM(50), DY2K0S(50),
     &     DY2LA(50), DY2LB(50), DY2PHI(50),
     &     dm2pip(50), dm2pim(50), DMT2PR(50),
     &     DMT2PB(50), DMT2KP(50), dm2km(50),
     &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
     &     dy2msn(50), DY2PIP(50), DY2PIM(50), 
     &     DY2PI0(50), DY2PR(50), DY2PB(50)
     &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
cc      SAVE /ARANA2/
      SAVE   
cbz3/17/99 end
cms      OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN')
cms      OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN')
cms      OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN')
cms      OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN')
cms      OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN')
cms      OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN')
cms      OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN')
cms      OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN')
cms      OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN')
clin-4/24/03
c      OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN')
cms      OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN')
cms      OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN')
cms      OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN')
cms      OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN')
cbz3/17/99
cms      OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN')
cms      OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN')
cms      OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN')
cms      OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN')
cms      OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN')
cms      OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN')
cbz3/17/99 end

cms      OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN')
cms      OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN')
cms      OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN')
cms      OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN')
cms      OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN')
cms      OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN')
cms      OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN')
cms      OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN')
cms      OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN')
clin-9/26/03 no symmetrization in y (or eta) for ana dat files
c      SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0
      SCALE1 = 1. / REAL(NEVNT * NUM) / BY
      SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1)
c
      DO 1001 I = 1, 50
         ymid=-10.+BY * (I - 0.5)
cms         WRITE (30, 333) ymid, SCALE1 * dy1ntb(I)
cms         WRITE (31, 333) ymid, SCALE1 * dy1ntp(I)
cms         WRITE (32, 333) ymid, SCALE1 * DY1HM(I)
cms         WRITE (37, 333) ymid, SCALE1 * DY1KP(I)
cms         WRITE (38, 333) ymid, SCALE1 * DY1KM(I)
cms         WRITE (39, 333) ymid, SCALE1 * DY1K0S(I)
cms         WRITE (40, 333) ymid, SCALE1 * DY1LA(I)
cms         WRITE (41, 333) ymid, SCALE1 * DY1LB(I)
cms         WRITE (42, 333) ymid, SCALE1 * DY1PHI(I)
cms         WRITE (33, 333) ymid, SCALE1 * DY1NEG(I)
cms         WRITE (34, 333) ymid, SCALE1 * DY1CH(I)
cms         WRITE (35, 333) ymid, SCALE1 * DE1NEG(I)
cms         WRITE (36, 333) ymid, SCALE1 * DE1CH(I)
cms         WRITE (43, 333) ymid, SCALE1 * dy1msn(I)
cms         WRITE (44, 333) ymid, SCALE1 * DY1PIP(I)
cms         WRITE (45, 333) ymid, SCALE1 * DY1PIM(I)
cms         WRITE (46, 333) ymid, SCALE1 * DY1PI0(I)
cms         WRITE (47, 333) ymid, SCALE1 * DY1PR(I)
cms         WRITE (48, 333) ymid, SCALE1 * DY1PB(I)

         IF (dm1pip(I) .NE. 0.0) THEN
c            WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I)
         END IF
         IF (dm1pim(I) .NE. 0.0) THEN
c            WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * 
c     &         dm1pim(I)
         END IF
         IF (DMT1PR(I) .NE. 0.0) THEN
c            WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I)
         END IF
         IF (DMT1PB(I) .NE. 0.0) THEN
c            WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I)
         END IF
         IF (DMT1KP(I) .NE. 0.0) THEN
c            WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I)
         END IF
         IF (dm1km(I) .NE. 0.0) THEN
c            WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
c     &         dm1km(I)
         END IF
         IF (dm1k0s(I) .NE. 0.0) THEN
c            WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I)
         END IF
         IF (DMT1LA(I) .NE. 0.0) THEN
c            WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I)
         END IF
         IF (DMT1LB(I) .NE. 0.0) THEN
c            WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I)
         END IF
 1001 CONTINUE
c
      DO 1002 I = 30, 48
c         WRITE (I, *) 'after hadron evolution'
 1002    CONTINUE
      DO 1003 I = 50, 58
c         WRITE (I, *) 'after hadron evolution'
 1003 CONTINUE

      DO 1004 I = 1, 50
         ymid=-10.+BY * (I - 0.5)
cms         WRITE (30, 333) ymid, SCALE1 * dy2ntb(I)
cms         WRITE (31, 333) ymid, SCALE1 * dy2ntp(I)
cms         WRITE (32, 333) ymid, SCALE1 * DY2HM(I)
cms         WRITE (37, 333) ymid, SCALE1 * DY2KP(I)
cms         WRITE (38, 333) ymid, SCALE1 * DY2KM(I)
cms         WRITE (39, 333) ymid, SCALE1 * DY2K0S(I)
cms         WRITE (40, 333) ymid, SCALE1 * DY2LA(I)
cms         WRITE (41, 333) ymid, SCALE1 * DY2LB(I)
cms         WRITE (42, 333) ymid, SCALE1 * DY2PHI(I)
cms         WRITE (33, 333) ymid, SCALE1 * DY2NEG(I)
cms         WRITE (34, 333) ymid, SCALE1 * DY2CH(I)
cms         WRITE (35, 333) ymid, SCALE1 * DE2NEG(I)
cms         WRITE (36, 333) ymid, SCALE1 * DE2CH(I)
cms         WRITE (43, 333) ymid, SCALE1 * dy2msn(I)
cms         WRITE (44, 333) ymid, SCALE1 * DY2PIP(I)
cms         WRITE (45, 333) ymid, SCALE1 * DY2PIM(I)
cms         WRITE (46, 333) ymid, SCALE1 * DY2PI0(I)
cms         WRITE (47, 333) ymid, SCALE1 * DY2PR(I)
cms         WRITE (48, 333) ymid, SCALE1 * DY2PB(I)
c
         IF (dm2pip(I) .NE. 0.0) THEN
c            WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I)
         END IF
         IF (dm2pim(I) .NE. 0.0) THEN
c            WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 * 
c     &         dm2pim(I)
         END IF
         IF (DMT2PR(I) .NE. 0.0) THEN
c            WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I)
         END IF
         IF (DMT2PB(I) .NE. 0.0) THEN
c            WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I)
         END IF
         IF (DMT2KP(I) .NE. 0.0) THEN
c            WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I)
         END IF
         IF (dm2km(I) .NE. 0.0) THEN
c            WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 * 
c     &         dm2km(I)
         END IF
         IF (dm2k0s(I) .NE. 0.0) THEN
c            WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I)
         END IF
         IF (DMT2LA(I) .NE. 0.0) THEN
c            WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I)
         END IF
         IF (DMT2LB(I) .NE. 0.0) THEN
c            WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I)
         END IF
 1004 CONTINUE
c      format(2(f12.5,1x))

      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine in HIJING before parton cascade evolution
      SUBROUTINE HJANA1

      PARAMETER (YMAX = 1.0, YMIN = -1.0)
      PARAMETER (DMT = 0.05, DY = 0.2)
      PARAMETER (DR = 0.2)
      PARAMETER (MAXPTN=400001,MAXSTR=150001)
      DIMENSION dyp1(50), DMYP1(200), DEYP1(50)
      DIMENSION dyg1(50), DMYG1(200), DEYG1(50)
      DIMENSION SNYP1(50), SMYP1(200), SEYP1(50)
      DIMENSION SNYG1(50), SMYG1(200), SEYG1(50)
      DIMENSION dnrpj1(50), dnrtg1(50), dnrin1(50),
     &   dnrtt1(50)
      DIMENSION dyg1c(50), dmyg1c(50), deyg1c(50)
      DIMENSION snrpj1(50), snrtg1(50), snrin1(50),
     &   snrtt1(50)
      DIMENSION snyg1c(50), smyg1c(50), seyg1c(50)
      DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0

      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc      SAVE /HPARNT/
      COMMON/hjcrdn/YP(3,300),YT(3,300)
cc      SAVE /hjcrdn/
      COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
     &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
     &   PJPM(300,500),NTJ(300),KFTJ(300,500),
     &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
     &   PJTE(300,500),PJTM(300,500)
cc      SAVE /HJJET1/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
      COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
     &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
     &     XMASS0(MAXPTN), ITYP0(MAXPTN)
cc      SAVE /prec1/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 200
            DMYP1(I) = SMYP1(I)
            DMYG1(I) = SMYG1(I)
 1001    CONTINUE

         DO 1002 I = 1, 50
            dyp1(I) = SNYP1(I)
            DEYP1(I) = SEYP1(I)
            dyg1(I) = SNYG1(I)
            DEYG1(I) = SEYG1(I)
            dnrpj1(I) = snrpj1(I)
            dnrtg1(I) = snrtg1(I)
            dnrin1(I) = snrin1(I)
            dnrtt1(I) = snrtt1(I)
            dyg1c(I) = snyg1c(I)
            dmyg1c(I) = smyg1c(I)
            deyg1c(I) = seyg1c(I)
 1002    CONTINUE
         nsubp = nsubpS
         nsubg = nsubgS
         NISG = NISGS
      ELSE
         DO 1003 I = 1, 200
            SMYP1(I) = DMYP1(I)
            SMYG1(I) = DMYG1(I)
 1003    CONTINUE

         DO 1004 I = 1, 50
            SNYP1(I) = dyp1(I)
            SEYP1(I) = DEYP1(I)
            SNYG1(I) = dyg1(I)
            SEYG1(I) = DEYG1(I)
            snrpj1(I) = dnrpj1(I)
            snrtg1(I) = dnrtg1(I)
            snrin1(I) = dnrin1(I)
            snrtt1(I) = dnrtt1(I)
            snyg1c(I) = dyg1c(I)
            smyg1c(I) = dmyg1c(I)
            seyg1c(I) = deyg1c(I)
 1004    CONTINUE
         nsubpS = nsubp
         nsubgS = nsubg
         NISGS = NISG
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF
c.....analysis
      DO 1006 I = 1, IHNT2(1)
         DO 1005 J = 1, NPJ(I)
            ITYP = KFPJ(I, J)
            PX = PJPX(I, J)
            PY = PJPY(I, J)
            PZ = PJPZ(I, J)
            PE = PJPE(I, J)
            PM = PJPM(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA1, PROJ STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 200
c            END IF
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 100
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
            dyp1(IY) = dyp1(IY) + 1.0
            DEYP1(IY) = DEYP1(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg1(IY) = dyg1(IY) + 1.0
               DEYG1(IY) = DEYG1(IY) + XMT
            END IF
 100        CONTINUE
            IMT = 1 + int(DXMT / DMT)
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200
            IF (IMT .GT. 200) GOTO 200
            DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
            END IF
 200        CONTINUE
 1005    CONTINUE
 1006 CONTINUE

      DO 1008 I = 1, IHNT2(3)
         DO 1007 J = 1, NTJ(I)
            ITYP = KFTJ(I, J)
            PX = PJTX(I, J)
            PY = PJTY(I, J)
            PZ = PJTZ(I, J)
            PE = PJTE(I, J)
            PM = PJTM(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA1, TARG STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 400
c            END IF
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA1 mt=0'
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 300
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 300
            dyp1(IY) = dyp1(IY) + 1.0
            DEYP1(IY) = DEYP1(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg1(IY) = dyg1(IY) + 1.0
               DEYG1(IY) = DEYG1(IY) + XMT
            END IF
 300        CONTINUE
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 200) GOTO 400
            DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
            END IF
 400        CONTINUE
 1007    CONTINUE
 1008 CONTINUE

      DO 1010 I = 1, NSG
         DO 1009 J = 1, NJSG(I)
            ITYP = K2SG(I, J)
            PX = PXSG(I, J)
            PY = PYSG(I, J)
            PZ = PZSG(I, J)
            PE = PESG(I, J)
            PM = PMSG(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA1, INDP STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 600
c            END IF
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA1 mt=0'
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 500
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 500
            dyp1(IY) = dyp1(IY) + 1.0
            DEYP1(IY) = DEYP1(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg1(IY) = dyg1(IY) + 1.0
               DEYG1(IY) = DEYG1(IY) + XMT
            END IF
 500        CONTINUE
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 200) GOTO 600
            DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
            END IF
 600        CONTINUE
 1009    CONTINUE
 1010 CONTINUE

      DO 1011 I = 1, IHNT2(1)
         YR = SQRT(YP(1, I) ** 2 + YP(2, I) ** 2)
         IR = 1 + int(YR / DR)
clin-4/2008 protect against out-of-bound errors:
c         IF (IR .GT. 50) GOTO 601
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601
         dnrpj1(IR) = dnrpj1(IR) + 1.0
         dnrtt1(IR) = dnrtt1(IR) + 1.0
 601     CONTINUE
 1011 CONTINUE

      DO 1012 I = 1, IHNT2(3)
         YR = SQRT(YT(1, I) ** 2 + YT(2, I) ** 2)
         IR = 1 + int(YR / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 602
         dnrtg1(IR) = dnrtg1(IR) + 1.0
         dnrtt1(IR) = dnrtt1(IR) + 1.0
 602     CONTINUE
 1012 CONTINUE

      DO 1013 I = 1, NSG
         Y1 = 0.5 * (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))
         Y2 = 0.5 * (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))
         YR = SQRT(Y1 ** 2 + Y2 ** 2)
         IR = 1 + int(YR / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603
         dnrin1(IR) = dnrin1(IR) + 1.0
         dnrtt1(IR) = dnrtt1(IR) + 1.0
 603     CONTINUE
 1013 CONTINUE

      DO 1014 I = 1, MUL
         ITYP = ITYP0(I)
         PX = sngl(PX0(I))
         PY = sngl(PY0(I))
         PZ = sngl(PZ0(I))
         PE = sngl(E0(I))
         PM = sngl(XMASS0(I))
         XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
         DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
c         IF (ABS(PZ) .GE. PE) THEN
c            PRINT *, ' IN HJANA1, GLUON ', I
c            PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c            PRINT *, ' PZ = ', PZ, ' EE = ', PE
c            PRINT *, ' XM = ', PM
c            GOTO 800
c         END IF
c         RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
         if(XMT.gt.0.) then
            RAP=asinh(PZ/XMT)
         else
            PRINT *, ' IN HJANA1 mt=0'
            RAP = 1000000.0*sign(1.,PZ)
         endif

         IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c         IF (IY .GT. 50) GOTO 700
         IF (IY.lt.1 .or.IY .GT. 50) GOTO 700
         dyg1c(IY) = dyg1c(IY) + 1.0
         deyg1c(IY) = deyg1c(IY) + XMT
 700     CONTINUE
         IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800
         IMT = 1 + int(DXMT / DMT)
         IF (IMT .GT. 50) GOTO 800
         dmyg1c(IMT) = dmyg1c(IMT) + 1.0 / XMT
 800     CONTINUE
 1014 CONTINUE
c.....count number of particles
      DO 1016 I = 1, IHNT2(1)
         DO 1015 J = 1, NPJ(I)
            nsubp = nsubp + 1
            IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
 1015    CONTINUE
 1016 CONTINUE

      DO 1018 I = 1, IHNT2(3)
         DO 1017 J = 1, NTJ(I)
            nsubp = nsubp + 1
            IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
 1017    CONTINUE
 1018 CONTINUE

      DO 1020 I = 1, NSG
         DO 1019 J = 1, NJSG(I)
            nsubp = nsubp + 1
            IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
 1019    CONTINUE
 1020 CONTINUE
      NISG = NISG + NSG
      IF (IOUT .EQ. 1) THEN
cbzdbg2/16/99
c      PRINT *, ' in HJANA1 '
c      PRINT *, ' total number of partons = ', nsubp
c      PRINT *, ' total number of gluons = ', nsubg, MUL
c      PRINT *, ' number of projectile strings = ', IHNT2(1)
c      PRINT *, ' number of target strings = ', IHNT2(3)
c      PRINT *, ' number of independent strings = ', NSG
      PRINT *, ' in HJANA1 '
      PRINT *, ' total number of partons = ', nsubp / IW
      PRINT *, ' total number of gluons = ', nsubg / IW
c      PRINT *, ' number of projectile strings = ', IHNT2(1)
c      PRINT *, ' number of target strings = ', IHNT2(3)
      PRINT *, ' number of independent strings = ', NISG / IW
cbzdbg2/16/99end
      END IF
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine in ZPC after generation of additional initial
c.....phase space distributions.

      SUBROUTINE HJAN1A
      PARAMETER (MAXPTN=400001)
      PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2)
      DIMENSION dgxg1a(50), dgyg1a(50), dtg1a(50)
      DIMENSION sgxg1a(50), sgyg1a(50), stg1a(50)
      DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 50
            dgxg1a(I) = sgxg1a(I)
            dgyg1a(I) = sgyg1a(I)
            dtg1a(I) = stg1a(I)
 1001    CONTINUE
      ELSE
         DO 1002 I = 1, 50
            sgxg1a(I) = dgxg1a(I)
            sgyg1a(I) = dgyg1a(I)
            stg1a(I) = dtg1a(I)
 1002    CONTINUE
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF
c.....analysis
      DO 1003 I = 1, MUL
         IGX = 1 + int(sngl(ABS(GX5(I))) / DGX)
clin-4/2008 protect against out-of-bound errors:
c         IF (IGX .GT. 50) GOTO 100
         IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100
         dgxg1a(IGX) = dgxg1a(IGX) + 1.0
 100     CONTINUE
         IGY = 1 + int(sngl(ABS(GY5(I))) / DGY)
         IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200
         dgyg1a(IGY) = dgyg1a(IGY) + 1.0
 200     CONTINUE
         IT = 1 + int(sngl(SQRT(FT5(I) ** 2 - GZ5(I) ** 2)) / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 300
         dtg1a(IT) = dtg1a(IT) + 1.0
 300     CONTINUE
 1003 CONTINUE
      CALL HJAN1B
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine in HJAN1A

      SUBROUTINE HJAN1B
      PARAMETER (MAXPTN=400001,MAXSTR=150001)
      PARAMETER (DR = 0.2, DT = 0.2)
      DIMENSION DNRG1B(50), dtg1b(50)
      DIMENSION SNRG1B(50), stg1b(50)
      DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
cc      SAVE /ilist8/
      COMMON /SREC1/ NSP, NST, NSI
cc      SAVE /SREC1/
      COMMON/hjcrdn/YP(3,300),YT(3,300)
cc      SAVE /hjcrdn/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 50
            DNRG1B(I) = SNRG1B(I)
            dtg1b(I) = stg1b(I)
 1001    CONTINUE
      ELSE
         DO 1002 I = 1, 50
            SNRG1B(I) = DNRG1B(I)
            stg1b(I) = dtg1b(I)
 1002    CONTINUE
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF
c.....analysis
      DO 1003 I = 1, MUL
         J = LSTRG1(I)

         IF (J .LE. NSP) THEN
            K = J
            GX0 = YP(1, J)
            GY0 = YP(2, J)
         ELSE IF (J .LE. NSP + NST) THEN
            K = J - NSP
            GX0 = YT(1, K)
            GY0 = YT(2, K)
         ELSE
            K = J - NSP - NST
            GX0 = 0.5 * (YP(1, IASG(K, 1)) + YT(1, IASG(K, 2)))
            GY0 = 0.5 * (YP(2, IASG(K, 1)) + YT(2, IASG(K, 2)))
         END IF
         R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2)
         IR = 1 + int(R0 / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100
         DNRG1B(IR) = DNRG1B(IR) + 1.0
 100     CONTINUE
         TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2))
         IT = 1 + int(TAU7 / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 200
         dtg1b(IT) = dtg1b(IT) + 1.0
 200     CONTINUE
 1003 CONTINUE
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine in HIJING after parton cascade evolution
      SUBROUTINE HJANA2
c
      PARAMETER (YMAX = 1.0, YMIN = -1.0)
      PARAMETER (DMT = 0.05, DY = 0.2)
      PARAMETER (DR = 0.2, DT = 0.2)
      PARAMETER (MAXPTN=400001)
      PARAMETER (MAXSTR=150001)
      DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
     1     GXSGS,GYSGS,GZSGS,FTSGS
      DIMENSION dyp2(50), DMYP2(200), DEYP2(50)
      DIMENSION dyg2(50), DMYG2(200), DEYG2(50)
      DIMENSION SNYP2(50), SMYP2(200), SEYP2(50)
      DIMENSION SNYG2(50), SMYG2(200), SEYG2(50)
      DIMENSION dnrpj2(50), dnrtg2(50), dnrin2(50),
     &   dnrtt2(50)
      DIMENSION dtpj2(50), dttg2(50), dtin2(50),
     &   dttot2(50)
      DIMENSION dyg2c(50), dmyg2c(50), deyg2c(50)
      DIMENSION snrpj2(50), snrtg2(50), snrin2(50),
     &   snrtt2(50)
      DIMENSION stpj2(50), sttg2(50), stin2(50),
     &   sttot2(50)
      DIMENSION snyg2c(50), smyg2c(50), seyg2c(50)
      DOUBLE PRECISION  ATAUI, ZT1, ZT2, ZT3
      DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc      SAVE /HPARNT/
      COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
cc      SAVE /SREC2/
      COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
     &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
     &   PJPM(300,500),NTJ(300),KFTJ(300,500),
     &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
     &   PJTE(300,500),PJTM(300,500)
cc      SAVE /HJJET1/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      common/anim/nevent,isoft,isflag,izpc
cc      SAVE /anim/
      COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
     &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
     &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
     &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
cc      SAVE /SOFT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 200
            DMYP2(I) = SMYP2(I)
            DMYG2(I) = SMYG2(I)
 1001    CONTINUE

         DO 1002 I = 1, 50
            dyp2(I) = SNYP2(I)
            DEYP2(I) = SEYP2(I)
            dyg2(I) = SNYG2(I)
            DEYG2(I) = SEYG2(I)
            dnrpj2(I) = snrpj2(I)
            dnrtg2(I) = snrtg2(I)
            dnrin2(I) = snrin2(I)
            dnrtt2(I) = snrtt2(I)
            dtpj2(I) = stpj2(I)
            dttg2(I) = sttg2(I)
            dtin2(I) = stin2(I)
            dttot2(I) = sttot2(I)
            dyg2c(I) = snyg2c(I)
            dmyg2c(I) = smyg2c(I)
            deyg2c(I) = seyg2c(I)
 1002    CONTINUE
         nsubp = nsubpS
         nsubg = nsubgS
         NISG = NISGS
      ELSE
         DO 1003 I = 1, 200
            SMYP2(I) = DMYP2(I)
            SMYG2(I) = DMYG2(I)
 1003    CONTINUE

         DO 1004 I = 1, 50
            SNYP2(I) = dyp2(I)
            SEYP2(I) = DEYP2(I)
            SNYG2(I) = dyg2(I)
            SEYG2(I) = DEYG2(I)
            snrpj2(I) = dnrpj2(I)
            snrtg2(I) = dnrtg2(I)
            snrin2(I) = dnrin2(I)
            snrtt2(I) = dnrtt2(I)
            stpj2(I) = dtpj2(I)
            sttg2(I) = dttg2(I)
            stin2(I) = dtin2(I)
            sttot2(I) = dttot2(I)
            snyg2c(I) = dyg2c(I)
            smyg2c(I) = dmyg2c(I)
            seyg2c(I) = deyg2c(I)
 1004    CONTINUE
         nsubpS = nsubp
         nsubgS = nsubg
         NISGS = NISG
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF

clin-4/28/01:
      if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 510

c.....analysis
      DO 1006 I = 1, IHNT2(1)
         DO 1005 J = 1, NPJ(I)
            ITYP = KFPJ(I, J)
            PX = PJPX(I, J)
            PY = PJPY(I, J)
            PZ = PJPZ(I, J)
            PE = PJPE(I, J)
            PM = PJPM(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
ccbzdbg2/16/99
cc            IF (ABS(PZ) .GE. PE) GOTO 200
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 200
c            END IF
ccbzdbg2/16/99end
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA2 mt=0'
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 100
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
            dyp2(IY) = dyp2(IY) + 1.0
            DEYP2(IY) = DEYP2(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg2(IY) = dyg2(IY) + 1.0
               DEYG2(IY) = DEYG2(IY) + XMT
            END IF
 100        CONTINUE
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 200) GOTO 200
            DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
            END IF
 200        CONTINUE
 1005    CONTINUE
 1006 CONTINUE

      DO 1008 I = 1, IHNT2(3)
         DO 1007 J = 1, NTJ(I)
            ITYP = KFTJ(I, J)
            PX = PJTX(I, J)
            PY = PJTY(I, J)
            PZ = PJTZ(I, J)
            PE = PJTE(I, J)
            PM = PJTM(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
clin-9/2012 determine rapidity more generally:
ccbzdbg2/16/99
cc            IF (ABS(PZ) .GE. PE) GOTO 400
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 400
c            END IF
ccbzdbg2/16/99end
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA2 mt=0'
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 300
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 300
            dyp2(IY) = dyp2(IY) + 1.0
            DEYP2(IY) = DEYP2(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg2(IY) = dyg2(IY) + 1.0
               DEYG2(IY) = DEYG2(IY) + XMT
            END IF
 300        CONTINUE
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 200) GOTO 400
            DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
            END IF
 400        CONTINUE
 1007    CONTINUE
 1008 CONTINUE

clin-4/28/01:
 510  continue

      DO 1010 I = 1, NSG
clin-4/25/01 soft3:
c         DO J = 1, NJSG(I)
         NJ=NJSG(I)
         if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
         DO 1009 J = 1, NJ
clin-4/25/01-end

            ITYP = K2SG(I, J)
            PX = PXSG(I, J)
            PY = PYSG(I, J)
            PZ = PZSG(I, J)
            PE = PESG(I, J)
            PM = PMSG(I, J)
clin-4/25/01 soft3:
            if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
               ITYP = K2SGS(I, J)
               PX = sngl(PXSGS(I, J))
               PY = sngl(PYSGS(I, J))
               PZ = sngl(PZSGS(I, J))
               PE = sngl(PESGS(I, J))
               PM = sngl(PMSGS(I, J))
            endif
clin-4/25/01-end

clin-9/2012 determine rapidity more generally:
            XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
            DXMT = XMT - PM
ccbzdbg2/16/99
cc            IF (ABS(PZ) .GE. PE) GOTO 600
c            IF (ABS(PZ) .GE. PE) THEN
c               PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
c               PRINT *, ' XM = ', PM
c               GOTO 600
c            END IF
ccbzdbg2/16/99end
c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
            if(XMT.gt.0.) then
               RAP=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA2 mt=0'
               RAP = 1000000.0*sign(1.,PZ)
            endif

            IY = 1 + int(ABS(RAP) / DY)
clin-8/2014 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 500
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 500
            dyp2(IY) = dyp2(IY) + 1.0
            DEYP2(IY) = DEYP2(IY) + XMT
            IF (ITYP .EQ. 21) THEN
               dyg2(IY) = dyg2(IY) + 1.0
               DEYG2(IY) = DEYG2(IY) + XMT
            END IF
 500        CONTINUE
            IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 200) GOTO 600
            DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
            IF (ITYP .EQ. 21) THEN
               DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
            END IF
 600        CONTINUE
 1009    CONTINUE
 1010 CONTINUE

clin-4/28/01:
      if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 520

      DO 1011 I = 1, IHNT2(1)
         J = I
         YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
         IR = 1 + int(YR / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601
         dnrpj2(IR) = dnrpj2(IR) + 1.0
         dnrtt2(IR) = dnrtt2(IR) + 1.0
 601     CONTINUE
         IT = 1 + int(sngl(ATAUI(J)) / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 602
         dtpj2(IT) = dtpj2(IT) + 1.0
         dttot2(IT) = dttot2(IT) + 1.0
 602     CONTINUE
 1011 CONTINUE

      DO 1012 I = 1, IHNT2(3)
         J = I + IHNT2(1)
         YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
         IR = 1 + int(YR / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603
         dnrtg2(IR) = dnrtg2(IR) + 1.0
         dnrtt2(IR) = dnrtt2(IR) + 1.0
 603     CONTINUE
         IT = 1 + int(sngl(ATAUI(J)) / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 604
         dttg2(IT) = dttg2(IT) + 1.0
         dttot2(IT) = dttot2(IT) + 1.0
 604     CONTINUE
 1012 CONTINUE

clin-4/28/01:
 520  continue

      DO 1013 I = 1, NSG
         J = I + IHNT2(1) + IHNT2(3)
clin-4/28/01:
         if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) J = I

         YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
         IR = 1 + int(YR / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 605
         dnrin2(IR) = dnrin2(IR) + 1.0
         dnrtt2(IR) = dnrtt2(IR) + 1.0
 605     CONTINUE
         IT = 1 + int(sngl(ATAUI(J)) / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 606
         dtin2(IT) = dtin2(IT) + 1.0
         dttot2(IT) = dttot2(IT) + 1.0
 606     CONTINUE
 1013 CONTINUE

      DO 1014 I = 1, MUL
         ITYP = ITYP5(I)
         PX = sngl(PX5(I))
         PY = sngl(PY5(I))
         PZ = sngl(PZ5(I))
         PE = sngl(E5(I))
         PM = sngl(XMASS5(I))
clin-9/2012 determine rapidity more generally:
         XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
         DXMT = XMT - PM
ccbzdbg2/16/99
cc            IF (ABS(PZ) .GE. PE) GOTO 800
c         IF (ABS(PZ) .GE. PE) THEN
c            PRINT *, ' IN HJANA2, GLUON ', I
c            PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c            PRINT *, ' PZ = ', PZ, ' EE = ', PE
c            PRINT *, ' XM = ', PM
c            GOTO 800
c         END IF
ccbzdbg2/16/99end
c         RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
         if(XMT.gt.0.) then
            RAP=asinh(PZ/XMT)
         else
            PRINT *, ' IN HJANA2 mt=0'
            RAP = 1000000.0*sign(1.,PZ)
         endif

         IY = 1 + int(ABS(RAP) / DY)
clin-9/2012 prevent possible segmentation fault (due to IY<=0):
c         IF (IY .GT. 50) GOTO 700
         IF (IY.lt.1 .or.IY .GT. 50) GOTO 700
         dyg2c(IY) = dyg2c(IY) + 1.0
         deyg2c(IY) = deyg2c(IY) + XMT
 700     CONTINUE
         IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800
         IMT = 1 + int(DXMT / DMT)
         IF (IMT .GT. 50) GOTO 800
         dmyg2c(IMT) = dmyg2c(IMT) + 1.0 / XMT
 800     CONTINUE
 1014 CONTINUE

clin-4/25/01 soft3:
      if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 530

c.....count number of particles
      DO 1016 I = 1, IHNT2(1)
         DO 1015 J = 1, NPJ(I)
            nsubp = nsubp + 1
            IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
 1015    CONTINUE
 1016 CONTINUE

      DO 1018 I = 1, IHNT2(3)
         DO 1017 J = 1, NTJ(I)
            nsubp = nsubp + 1
            IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
 1017    CONTINUE
 1018 CONTINUE

clin-4/25/01 soft3:
 530  continue

      DO 1020 I = 1, NSG
clin-4/25/01 soft3:
c         DO J = 1, NJSG(I)
         NJ=NJSG(I)
         if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
         DO 1019 J = 1, NJ
clin-4/25/01-end

            nsubp = nsubp + 1

clin-4/25/01
c            IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
            if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
               IF(K2SGS(I, J) .EQ. 21) nsubg = nsubg + 1
            else
               IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
            endif
clin-4/25/01-end
 1019    CONTINUE
 1020 CONTINUE
cbzdbg2/16/99
      NISG = NISG + NSG

      IF (IOUT .EQ. 1) THEN
cbzdbg2/16/99end
cbzdbg2/16/99
c      PRINT *, ' in HJANA2 '
c      PRINT *, ' total number of partons = ', nsubp
c      PRINT *, ' total number of gluons = ', nsubg, MUL
c      PRINT *, ' number of projectile strings = ', IHNT2(1)
c      PRINT *, ' number of target strings = ', IHNT2(3)
c      PRINT *, ' number of independent strings = ', NSG
      PRINT *, ' in HJANA2 '
      PRINT *, ' total number of partons = ', nsubp / IW
      PRINT *, ' total number of gluons = ', nsubg / IW
c      PRINT *, ' number of projectile strings = ', IHNT2(1)
c      PRINT *, ' number of target strings = ', IHNT2(3)
      PRINT *, ' number of independent strings = ', NISG / IW
      END IF

      CALL HJAN2A
      CALL HJAN2B

      RETURN
      END

c-----------------------------------------------------------------------

c.....subroutine called by HJANA2
      SUBROUTINE HJAN2A

      PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2)
      PARAMETER (MAXPTN=400001,MAXSTR=150001)
      DIMENSION dgxp2a(50), dgyp2a(50), dtp2a(50)
      DIMENSION dgxg2a(50), dgyg2a(50), dtg2a(50)
      DIMENSION sgxp2a(50), sgyp2a(50), stp2a(50)
      DIMENSION sgxg2a(50), sgyg2a(50), stg2a(50)
      DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc      SAVE /HPARNT/
      COMMON/hjcrdn/YP(3,300),YT(3,300)
cc      SAVE /hjcrdn/
      COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
     &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
     &   PJPM(300,500),NTJ(300),KFTJ(300,500),
     &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
     &   PJTE(300,500),PJTM(300,500)
cc      SAVE /HJJET1/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 50
            dgxp2a(I) = sgxp2a(I)
            dgyp2a(I) = sgyp2a(I)
            dtp2a(I) = stp2a(I)
            dgxg2a(I) = sgxg2a(I)
            dgyg2a(I) = sgyg2a(I)
            dtg2a(I) = stg2a(I)
 1001    CONTINUE
      ELSE
         DO 1002 I = 1, 50
            sgxp2a(I) = dgxp2a(I)
            sgyp2a(I) = dgyp2a(I)
            stp2a(I) = dtp2a(I)
            sgxg2a(I) = dgxg2a(I)
            sgyg2a(I) = dgyg2a(I)
            stg2a(I) = dtg2a(I)
 1002    CONTINUE
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF
c.....analysis
      DO 1004 I = 1, IHNT2(1)
         DO 1003 J = 1, NPJ(I)
            IF (KFPJ(I, J) .NE. 21) THEN
               IGX = 1 + int(ABS(YP(1, I)) / DGX)
               IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100
               dgxp2a(IGX) = dgxp2a(IGX) + 1.0
 100           CONTINUE
               IGY = 1 + int(ABS(YP(2, I)) / DGY)
               IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200
               dgyp2a(IGY) = dgyp2a(IGY) + 1.0
 200           CONTINUE
               IT = 1
               dtp2a(IT) = dtp2a(IT) + 1.0
            END IF
 1003    CONTINUE
 1004 CONTINUE

      DO 1006 I = 1, IHNT2(3)
         DO 1005 J = 1, NTJ(I)
            IF (KFTJ(I, J) .NE. 21) THEN
               IGX = 1 + int(ABS(YT(1, I)) / DGX)
               IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 300
               dgxp2a(IGX) = dgxp2a(IGX) + 1.0
 300           CONTINUE
               IGY = 1 + int(ABS(YT(2, I)) / DGY)
               IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 400
               dgyp2a(IGY) = dgyp2a(IGY) + 1.0
 400           CONTINUE
               IT = 1
               dtp2a(IT) = dtp2a(IT) + 1.0
            END IF
 1005    CONTINUE
 1006 CONTINUE

      DO 1008 I = 1, NSG
         DO 1007 J = 1, NJSG(I)
            IF (K2SG(I, J) .NE. 21) THEN
               IGX = 1 + int(ABS(0.5 * 
     &            (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))) / DGX)
               IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 500
               dgxp2a(IGX) = dgxp2a(IGX) + 1.0
 500           CONTINUE
               IGY = 1 + int(ABS(0.5 * 
     &            (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))) / DGY)
               IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 600
               dgyp2a(IGY) = dgyp2a(IGY) + 1.0
 600           CONTINUE
               IT = 1
               dtp2a(IT) = dtp2a(IT) + 1.0               
            END IF
 1007    CONTINUE
 1008 CONTINUE

      DO 1009 I = 1, MUL
         IGX = 1 + int(ABS(sngl(GX5(I))) / DGX)
         IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 700
         dgxg2a(IGX) = dgxg2a(IGX) + 1.0
         dgxp2a(IGX) = dgxp2a(IGX) + 1.0
 700     CONTINUE
         IGY = 1 + int(ABS(sngl(GY5(I))) / DGY)
         IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 800
         dgyg2a(IGY) = dgyg2a(IGY) + 1.0
         dgyp2a(IGY) = dgyp2a(IGY) + 1.0
 800     CONTINUE
         IT = 1 + int(SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2)) / DT)
         IF (IT .GT. 50 .or. IT .LT. 1) GOTO 900
         dtg2a(IT) = dtg2a(IT) + 1.0
         dtp2a(IT) = dtp2a(IT) + 1.0
 900     CONTINUE
 1009 CONTINUE
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine in HJANA2

      SUBROUTINE HJAN2B

      PARAMETER (MAXPTN=400001)
      PARAMETER (MAXSTR=150001)
      PARAMETER (DR = 0.2, DT = 0.2)
      DIMENSION DNRG2B(50), dtg2b(-24:25)
      DIMENSION SNRG2B(50), stg2b(-24:25)
      DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
      DOUBLE PRECISION  ATAUI, ZT1, ZT2, ZT3
      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
cc      SAVE /ilist8/
      COMMON /SREC1/ NSP, NST, NSI
cc      SAVE /SREC1/
      COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
cc      SAVE /SREC2/
      COMMON/hjcrdn/YP(3,300),YT(3,300)
cc      SAVE /hjcrdn/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
      COMMON /AREVT/ IAEVT, IARUN, MISS
cc      SAVE /AREVT/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      SAVE   
      DATA IW/0/

      IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
         DO 1001 I = 1, 50
            DNRG2B(I) = SNRG2B(I)
            dtg2b(I - 25) = stg2b(I - 25)
 1001    CONTINUE
      ELSE
         DO 1002 I = 1, 50
            SNRG2B(I) = DNRG2B(I)
            stg2b(I - 25) = dtg2b(I - 25)
 1002    CONTINUE
         isevt = IAEVT
         isrun = IARUN
         IW = IW + 1
      END IF
c.....analysis
      DO 1003 I = 1, MUL
         J = LSTRG1(I)
         GX0 = sngl(ZT1(J))
         GY0 = sngl(ZT2(J))
         R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2)
         IR = 1 + int(R0 / DR)
         IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100
         DNRG2B(IR) = DNRG2B(IR) + 1.0
 100     CONTINUE
         TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2))
         DTAU=TAU7 - sngl(ATAUI(J))
         IT = 1 + int(DTAU / DT)
cbzdbg2/21/99
c         IF (ABS(IT) .GT. 25) GOTO 200
         IF (IT .GT. 25 .OR. IT .LT. -24) GOTO 200
cbzdbg2/21/99end
         dtg2b(IT) = dtg2b(IT) + 1.0
 200     CONTINUE
 1003 CONTINUE
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine before ARTMN
      SUBROUTINE HJANA3
c
      PARAMETER (MAXSTR=150001, MAXR=1)
c.....y cut for mt spectrum
      PARAMETER (YMIN = -1.0, YMAX = 1.0)
cbz11/7/99 end
c.....bin width for mt spectrum and y spectrum
      PARAMETER (DMT = 0.05, DY = 0.2)
      DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult,
     1     v2bi,s2bi,xbmult
      DIMENSION dndyh3(50), DMYH3(50), DEYH3(50)
      COMMON /RUN/ NUM
cc      SAVE /RUN/
      COMMON /ARERC1/MULTI1(MAXR)
cc      SAVE /ARERC1/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cc      SAVE /ARPRC1/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
cc      SAVE /iflow/
      SAVE   
      DATA IW/0/

      IW = IW + 1
      DO 1002 J = 1, NUM
         DO 1001 I = 1, MULTI1(J)
            ITYP = ITYP1(I, J)
            IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
            PX = PX1(I, J)
            PY = PY1(I, J)
            PZ = PZ1(I, J)
            EE = EE1(I, J)
            XM = XM1(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
            DXMT = XMT - XM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. EE) THEN
c               PRINT *, 'IN HJANA3'
c               PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR'
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
c               PRINT *, ' XM = ', XM
c               GOTO 200
c            END IF
c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
            if(XMT.gt.0.) then
               Y=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA3 mt=0'
               Y = 1000000.0*sign(1.,PZ)
            endif

c.....rapidity cut for the rapidity distribution
c            IY = 1 + int(ABS(Y) / DY)
clin-8/2014 no rapidity shift here:
c            IY = 1 + int((Y+10.) / DY)
            IY = 1 + int(Y/DY)
clin-9/2012 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 100
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
            dndyh3(IY) = dndyh3(IY) + 1.0
            DEYH3(IY) = DEYH3(IY) + XMT
 100        CONTINUE
c.....insert rapidity cut for mt spectrum here
            IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 50) GOTO 200
            DMYH3(IMT) = DMYH3(IMT) + 1.0 / XMT
 200        CONTINUE
 1001    CONTINUE
 1002 CONTINUE
c
      RETURN
      END

c-----------------------------------------------------------------------

c.....analysis subroutine after ARTMN
      SUBROUTINE HJANA4
      PARAMETER (MAXSTR=150001, MAXR=1)
c.....y cut for mt spectrum
cbz11/7/99
c      PARAMETER (YMIN = -0.5, YMAX = 0.5)
      PARAMETER (YMIN = -1.0, YMAX = 1.0)
cbz11/7/99 end
c.....bin width for mt spectrum and y spectrum
      PARAMETER (DMT = 0.05, DY = 0.2)

      DIMENSION dndyh4(50), DMYH4(50), DEYH4(50)
      COMMON /RUN/ NUM
cc      SAVE /RUN/
      COMMON /ARERC1/MULTI1(MAXR)
cc      SAVE /ARERC1/
      COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
     &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR), 
     &     FT1(MAXSTR, MAXR),
     &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
     &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
cc      SAVE /ARPRC1/
      COMMON /AROUT/ IOUT
cc      SAVE /AROUT/
      COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi
cc      SAVE /fflow/
      SAVE   
      DATA IW/0/

      IW = IW + 1
      DO 1002 J = 1, NUM
         DO 1001 I = 1, MULTI1(J)
            ITYP = ITYP1(I, J)
            IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
            PX = PX1(I, J)
            PY = PY1(I, J)
            PZ = PZ1(I, J)
            EE = EE1(I, J)
            XM = XM1(I, J)
            XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
            DXMT = XMT - XM
clin-9/2012 determine rapidity more generally:
c            IF (ABS(PZ) .GE. EE) THEN
c               PRINT *, 'IN HJANA4'
c               PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR'
c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
c               PRINT *, ' XM = ', XM
c               GOTO 200
c            END IF
c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
            if(XMT.gt.0.) then
               Y=asinh(PZ/XMT)
            else
               PRINT *, ' IN HJANA4 mt=0'
               Y = 1000000.0*sign(1.,PZ)
            endif

c.....rapidity cut for the rapidity distribution
c            IY = 1 + int(ABS(Y) / DY)
clin-8/2014 no rapidity shift here:
c            IY = 1 + int((Y+10.) / DY)
            IY = 1 + int(Y/DY)
clin-9/2012 prevent possible segmentation fault (due to IY<=0):
c            IF (IY .GT. 50) GOTO 100
            IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
            dndyh4(IY) = dndyh4(IY) + 1.0
            DEYH4(IY) = DEYH4(IY) + XMT
 100        CONTINUE
c.....insert rapidity cut for mt spectrum here
            IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
            IMT = 1 + int(DXMT / DMT)
            IF (IMT .GT. 50) GOTO 200
            DMYH4(IMT) = DMYH4(IMT) + 1.0 / XMT
 200        CONTINUE
 1001    CONTINUE
 1002 CONTINUE
c
      RETURN
      END

c=======================================================================

c.....subroutine to get average values for different strings

      SUBROUTINE zpstrg

      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      PARAMETER (MAXPTN=400001)
      PARAMETER (MAXSTR=150001)
c      REAL*4 YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB
      REAL YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB

      COMMON /PARA1/ MUL
cc      SAVE /PARA1/
      COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
     &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
     &   XMASS5(MAXPTN), ITYP5(MAXPTN)
cc      SAVE /prec2/
      COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
cc      SAVE /ilist8/
      COMMON /SREC1/ NSP, NST, NSI
cc      SAVE /SREC1/
      COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
cc      SAVE /SREC2/
      COMMON/hjcrdn/YP(3,300),YT(3,300)
cc      SAVE /hjcrdn/
      COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
     &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
     &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
cc      SAVE /HJJET2/
cbz6/28/99 flow1
      COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
cc      SAVE /HPARNT/
cbz6/28/99 flow1 end
      common/anim/nevent,isoft,isflag,izpc
cc      SAVE /anim/
      common/strg/np(maxstr)
cc      SAVE /strg/
clin-6/06/02 test local freezeout:
      common /frzprc/ 
     &     gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
     &     pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
     &     xmfrz(MAXPTN), 
     &     tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
cc      SAVE /frzprc/
      SAVE   

clin-6/06/02 test local freezeout for string melting,
c     use space-time values at local freezeout saved in /frzprc/:
      if(isoft.eq.5) then
         do 1001 I = 1, MUL
            ITYP5(i)=idfrz(i)
            GX5(i)=gxfrz(i)
            GY5(i)=gyfrz(i)
            GZ5(i)=gzfrz(i)
            FT5(i)=ftfrz(i)
            PX5(i)=pxfrz(i)
            PY5(i)=pyfrz(i)
            PZ5(i)=pzfrz(i)
            E5(i)=efrz(i)
            XMASS5(i)=xmfrz(i)
 1001    continue
      endif
clin-6/06/02-end

      DO 1002 I = 1, MAXSTR
         ATAUI(I) = 0d0
         ZT1(I) = 0d0
         ZT2(I) = 0d0
clin-4/25/03 add zt3(I) to track longitudinal positions of partons/strings:
         ZT3(I) = 0d0
         NP(I) = 0
 1002 CONTINUE
      DO 1003 I = 1, MUL
         ISTRG = LSTRG1(I)
         TAU7 = SQRT(FT5(I) ** 2 - GZ5(I) ** 2)
         ATAUI(ISTRG) = ATAUI(ISTRG) + TAU7
         ZT1(ISTRG) = ZT1(ISTRG) + GX5(I)
         ZT2(ISTRG) = ZT2(ISTRG) + GY5(I)
         ZT3(ISTRG) = ZT3(ISTRG) + GZ5(I)
         NP(ISTRG) = NP(ISTRG) + 1
 1003 CONTINUE
      
      NSTR = NSP + NST + NSI

clin-7/03/01 correct averaging on transverse coordinates, no shift needed:
      if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
         DO 1004 I = 1, NSTR
            IF (NP(I) .NE. 0) THEN
               ATAUI(I) = ATAUI(I) / NP(I)
               ZT1(I) = ZT1(I) / NP(I)
               ZT2(I) = ZT2(I) / NP(I)
               ZT3(I) = ZT3(I) / NP(I)
            ENDIF
 1004    CONTINUE
         return
      endif
clin-7/03/01-end

      DO 1005 I = 1, NSTR
         IF (NP(I) .NE. 0) THEN
            ATAUI(I) = ATAUI(I) / NP(I)
            ZT1(I) = ZT1(I) / NP(I)
            ZT2(I) = ZT2(I) / NP(I)
            ZT3(I) = ZT3(I) / NP(I)
         ELSE
            IF (I .LE. NSP) THEN
               J = I
               ZT1(I) = dble(YP(1, J))
               ZT2(I) = dble(YP(2, J))
               ZT3(I) = 0d0
            ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN
               J = I - NSP
               ZT1(I) = dble(YT(1, J))
               ZT2(I) = dble(YT(2, J))
               ZT3(I) = 0d0
            ELSE
               J = I - NSP - NST
               ZT1(I) = 0.5d0*
     1              dble((YP(1, IASG(J, 1)) + YT(1, IASG(J, 2))))
               ZT2(I) = 0.5d0 *
     1              dble((YP(2, IASG(J, 1)) + YT(2, IASG(J, 2))))
               ZT3(I) = 0d0
            END IF
         END IF
 1005 CONTINUE

cbz6/28/99 flow1
      BB = HINT1(19)
      DO 1006 I = 1, NSTR
         IF (NP(I).NE.0) THEN
            SHIFT=0d0
         ELSE
            SHIFT=0.5d0*dble(BB)
         END IF
         IF (I .LE. NSP) THEN
            ZT1(I) = ZT1(I) + SHIFT
         ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN
            ZT1(I) = ZT1(I) - SHIFT
         END IF
 1006 CONTINUE
cbz6/28/99 flow1 end
c
      RETURN
      END

clin-10/01/03 random number generator for f77:
c      function RANART(NSEED)
c      SAVE   
clin-4/2008 ran(nseed) is renamed to avoid conflict with system functions:
c      ran=rand()
c      ranart=rand(0)
c     one may also use the following random number generator in PYTHIA/JETSET:
c      ranart=rlu(0)
c      return
c      end

clin-3/2009
c     Initialize hadron weights; 
c     Can add initial hadrons before the hadron cascade starts (but after ZPC).
      subroutine addhad
      PARAMETER (MAXSTR=150001,MAXR=1,xmd=1.8756)
      double precision  smearp,smearh
      COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
      COMMON /ARPRC/ ITYPAR(MAXSTR),
     &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
     &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
     &     XMAR(MAXSTR)
      COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
     1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
     2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
      COMMON /smearz/smearp,smearh
      COMMON/RNDF77/NSEED
      common /para8/ idpert,npertd,idxsec
      SAVE   

c     All hadrons at the start of hadron cascade have the weight of 1
c     except those inserted by the user in this subroutine:
      np0=IAINT2(1)
      DO i=1,np0
         dpertp(I)=1.
      ENDDO
c     Specify number, species, weight, initial x,p,m for inserted hadrons here:
      nadd=0
      tau0=ARPAR1(1)
      DO 100 i=np0+1,np0+nadd
         ITYPAR(I)=42
clin-5/2012 fix type mismatch:
c         dpertp(I)=1d0/dble(nadd)
         dpertp(I)=1./float(nadd)
         GXAR(I)=5.*(1.-2.*RANART(NSEED))
         GYAR(I)=5.*(1.-2.*RANART(NSEED))
         GZAR(I)=2.*(1.-2.*RANART(NSEED))
         FTAR(I)=0.
         PXAR(I)=1.
         PYAR(I)=0.
         PZAR(I)=1.
         XMAR(I)=xmd
c
         PEAR(I)=sqrt(PXAR(I)**2+PYAR(I)**2+PZAR(I)**2+XMAR(I)**2)
clin-9/2012 determine rapidity more generally:
c         RAP=0.5*alog((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
         RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
c
         VX=PXAR(I)/PEAR(I)
         VY=PYAR(I)/PEAR(I)
c.....give initial formation time shift and boost according to rapidity:
         TAUI=FTAR(I)+TAU0
         FTAR(I)=TAUI*COSH(RAP)
         GXAR(I)=GXAR(I)+VX*TAU0*COSH(RAP)
         GYAR(I)=GYAR(I)+VY*TAU0*COSH(RAP)
c     Allow the intial z-position to be different from the Bjorken picture:
         GZAR(I)=TAUI*SINH(RAP)+GZAR(I)
c         GZAR(I)=TAUI*SINH(RAP)
         zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
         GZAR(I)=GZAR(I)+zsmear
 100  CONTINUE
      IAINT2(1)=IAINT2(1)+nadd
c
      if(nadd.ge.1.and.idpert.ne.1.and.idpert.ne.2) then
         write(16,*) 'IDPERT must be 1 or 2 to add initial hadrons,
     1 set NPERTD to 0 if you do not need perturbative deuterons'
         stop
      endif
      if(IAINT2(1).gt.MAXSTR) then
         write(16,*) 'Too many initial hadrons, array size is exceeded!'
         stop
      endif
c
      return
      end

clin-8/2014 define function asinh():
      FUNCTION asinh(x)
      SAVE
      if(x.gt.0) then
         ASINH=alog(x+sqrt(x**2+1.))
      else
c     a la suggestion de YP Liu:
         ASINH=-alog(-x+sqrt(x**2+1.))
      endif
      return
      end