AMPTInterface/src/amptsub.f

0001 c....................amptsub.f
0002 c.....this file contains 4 sections:
0003 c.....1. ART subroutines;
0004 c.....2. ART functions;
0005 c.....3. ART block data;
0006 c.....4. subprocesses borrowed from other codes.
0007 c.....5. the previous artana.f
0008 c.....6. the previous zpcsub.f
0009 c.....7. subroutine getnp
0010 c.....Note that Parts1-4 are the previous artsub.f
0011 c
0012 c=======================================================================
0013 c.....subroutine to set up ART parameters and analysis files
0014 c.....before looping different events
0015       SUBROUTINE ARTSET
0016 c
0017       PARAMETER (AMU= 0.9383,nxymax=10001)
0018       double precision dpcoal,drcoal,ecritl
0019       INTEGER ZTA, ZPR
0020       common  /gg/      dx,dy,dz,dpx,dpy,dpz
0021 clin-10/03/03
0022 c     "SAVE   " (without argument) is used for most subroutines and functions,
0023 c     this is important for the success when using "f77" to compile:
0024 cc      SAVE /gg/
0025       common  /zz/      zta,zpr
0026 cc      SAVE /zz/
0027       COMMON  /RUN/     NUM
0028 cc      SAVE /RUN/
0029       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0030 cc      SAVE /input1/
0031       COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
0032      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
0033 cc      SAVE /INPUT2/
0034       COMMON /INPUT3/ PLAB, ELAB, ZEROPT, B0, BI, BM, DENCUT, CYCBOX
0035 cc      SAVE /INPUT3/
0036       common /imulst/ iperts
0037 cc      SAVE /imulst/
0038       common /coal/dpcoal,drcoal,ecritl
0039       common/anim/nevent,isoft,isflag,izpc
0040       common /para7/ ioscar,nsmbbbar,nsmmeson
0041       common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd,
0042      1     psembd,tmaxembd,phidecomp
0043       common/xyembed/nxyjet,xyjet(nxymax,2)
0044       SAVE
0045 clin-10/03/03  ecritl: local energy density below which a parton
0046 c     will freeze out (in GeV/fm^3), for improvements on string melting,
0047 c     not used in this version of AMPT:
0048 clin-4/2008
0049 c      data ecritl/1.d0/
0050       ecritl=1.d0
0051 c
0052 c     combine ART initialization into ampt.ini:
0053 c     (Note that the following values are relics from the old ART structure)
0054 c.....input parameter file
0055 c      OPEN(13, FILE = 'art1.ini', STATUS = 'UNKNOWN')
0056 c      READ (13, *) MASSTA, ZTA
0057       MASSTA=1
0058       ZTA=1
0059 c      write(12,*) massta, zta, ' massta, zta'
0060 c      READ (13, *) MASSPR, ZPR
0061       MASSPR=1
0062       ZPR=1
0063 c      write(12,*) masspr, zpr, ' masspr, zpr'
0064 c      READ (13, *) PLAB, IPLAB
0065       PLAB=14.6
0066       IPLAB=2
0067 c      write(12,*) plab, iplab, ' plab, iplab'
0068       if(iplab.eq.2)then
0069          elab=sqrt(plab**2+amu**2)-amu
0070       else
0071          elab=plab
0072       endif
0073       elab=elab*1000.
0074 c      READ (13, *) ZEROPT
0075       ZEROPT=0.
0076 c      write(12,*) zeropt, ' zeropt'
0077 clin-10/03/03 ISEED was used as a seed for random number inside ART,
0078 c     not used in AMPT:
0079       ISEED=700721
0080 c     0/1: (Normal or Perturbative) multistrange partice production.
0081 c     Perturbative option is disabled for now:
0082       iperts=0
0083 c      READ (13, *) MANYB, B0, BI, BM
0084 c     2/04/00 MANYB MUST BE SET TO 1 !
0085 c     in order to skip impact parameter setting by ART, then B0 has no effect.
0086       MANYB=1
0087       B0=1
0088       BI=0
0089       BM=0
0090 c      write(12,*) manyb, b0, bi, bm, ' manyb, b0, bi, bm'
0091 c      READ (13, *) ISEED
0092 c      write(12,*) iseed, ' iseed'
0093 c      READ (13, *) DT
0094 c      write(12,*) dt, ' dt'
0095 c      READ (13, *) NTMAX
0096 c      write(12,*) ntmax, ' ntmax'
0097 c      READ (13, *) ICOLL
0098       ICOLL=-1
0099 c      write(12,*) icoll, ' icoll'
0100 c      READ (13, *) NUM
0101 c     2/11/03 run events without test particles for now:
0102       NUM=1
0103 c      write(12,*) num, ' num'
0104 c      READ (13, *) INSYS
0105       INSYS=1
0106 c      write(12,*) insys, ' insys'
0107 c      READ (13, *) IPOT
0108       IPOT=3
0109 c      write(12,*) ipot, ' ipot'
0110 c      READ (13, *) MODE
0111       MODE=0
0112       IF(ICOLL.EQ.-1)IPOT=0
0113 c      write(12,*) mode, ' mode'
0114 c      READ (13, *) DX, DY, DZ
0115       DX=2.73
0116       DY=2.73
0117       DZ=2.73
0118 c      write(12,*) dx,dy,dz,' dx,dy,dz'
0119 c      READ (13, *) DPX, DPY, DPZ
0120       DPX=0.6
0121       DPY=0.6
0122       DPZ=0.6
0123 c      write(12,*) dpx,dpy,dpz,' dpx,dpy,dpz'
0124 c      READ (13, *) IAVOID
0125       IAVOID=1
0126 c      write(12,*) iavoid, ' iavoid'
0127 c      READ (13, *) IMOMEN
0128       IMOMEN=1
0129 c      write(12,*) imomen, ' imomen'
0130       if(icoll.eq.-1)imomen=3
0131 c      READ (13, *) NFREQ
0132       NFREQ=10
0133 c      write(12,*) nfreq, ' nfreq'
0134 c      READ (13, *) ICFLOW
0135       ICFLOW=0
0136 c      write(12,*) ICFLOW, ' ICFLOW'
0137 c      READ (13, *) ICRHO
0138       ICRHO=0
0139 c      write(12,*) ICRHO, ' ICRHO'
0140 c      READ (13, *) ICOU
0141       ICOU=0
0142 c      write(12,*)icou, ' icou'
0143 * kaon potential control parameter
0144 * KMUL IS A MULTIPLIER TO THE STANDARD K-N SCATTERING LENGTH
0145 c      READ (13, *) KPOTEN, KMUL
0146       KPOTEN=0
0147       KMUL=1
0148 c      write(12,*)kpoten,kmul, ' kpoten, kmul'
0149 * mean field control parameter FOR BARYONS
0150 * no mean filed is used for baryons if their
0151 * local density is higher than dencut.
0152 c      READ (13, *) DENCUT
0153       DENCUT=15
0154 c      write(12,*)dencut, ' dencut'
0155 * test reactions in a box of side-length cycbox
0156 * input cycbox
0157 c      READ (13, *) CYCBOX
0158       CYCBOX=0
0159 c      write(12,*) cycbox, ' cycbox'
0160 c
0161 clin-5b/2008
0162 c      if(ioscar.eq.2) then
0163       if(ioscar.eq.2.or.ioscar.eq.3) then
0164 c         OPEN (92,FILE='ana/parton-initial-afterPropagation.dat',
0165 c     1        STATUS = 'UNKNOWN')
0166       endif
0167       if(ioscar.eq.3) then
0168 clin-6/2009 write out full parton collision history:
0169 c         OPEN (95,FILE='ana/parton-collisionsHistory.dat',
0170 c     1        STATUS='UNKNOWN')
0171 clin-6/2009 write out initial minijet information:
0172 c         OPEN (96,FILE='ana/minijet-initial-beforePropagation.dat',
0173 c     1        STATUS='UNKNOWN')
0174 clin-6/2009 write out parton info after coalescence:
0175          if(isoft.eq.4.or.isoft.eq.5) then
0176 c            OPEN (85,FILE='ana/parton-after-coalescence.dat',
0177 c     1           STATUS='UNKNOWN')
0178          endif
0179       endif
0180 clin-6/2009 write out initial transverse positions of initial nucleons:
0181 c      OPEN (94,FILE='ana/npart-xy.dat',STATUS='UNKNOWN')
0182 c
0183 clin-8/2009 In case that random positions are used to embed high-Pt jets:
0184       if(iembed.eq.3.or.iembed.eq.4) then
0185 c         OPEN (97,FILE='embed-jet-xy.txt',STATUS = 'UNKNOWN')
0186 c         read(97,*) nxyjet
0187 c     Save positions in array to reuse when embedding more jet pairs
0188 c     than the number of entries in the position file:
0189          if(nevent.gt.nxyjet) then
0190             if(nxyjet.gt.nxymax) then
0191 c              print *, 'Too many lines in embed-jet-xy.txt:
0192 c     1 increase value of the parameter nxymax'
0193                stop
0194             elseif(nxyjet.le.0) then
0195 c               print *, 'Check number of entries in embed-jet-xy.txt'
0196                stop
0197             endif
0198             do ixy=1,nxyjet
0199 c               read(97,*) xyjet(ixy,1),xyjet(ixy,2)
0200             enddo
0201          endif
0202       endif
0203
0204       RETURN
0205       END
0206
0207 c-----------------------------------------------------------------------
0208
0209 c.....subroutine to initialize cascade.
0210
0211       SUBROUTINE ARINI
0212
0213 c.....before invoking ARINI:
0214 c.....IAPAR2(1), IAINT2(1) must be set.
0215       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0216 cc      SAVE /ARPRNT/
0217       SAVE
0218
0219 ctest off for resonance (phi, K*) studies:
0220 c      OPEN (89, FILE = 'ana/decay_rec.dat', STATUS = 'UNKNOWN')
0221
0222       IFLG = IAPAR2(1)
0223       GOTO (200, 200, 300) IFLG
0224
0225 c.....error choice of initialization
0226       PRINT *, 'IAPAR2(1) must be 1, 2, or 3'
0227       STOP
0228
0229 c.....to use default initial conditions generated by the cascade,
0230 c.....or to read in initial conditions.
0231  200  RETURN
0232
0233 c.....to generate formation time and the position at formation time from
0234 c.....read-in initial conditions with an averaged formation proper time.
0235  300  CALL ARINI1
0236 c.....ordering the particle label according to increasing order of
0237 c.....formation time.
0238       CALL ARTORD
0239       RETURN
0240
0241       END
0242
0243 c-----------------------------------------------------------------------
0244
0245 c.....subroutine to generate formation time and position at formation time
0246 c.....from read-in initial conditions with an averaged formation proper
0247 c.....time.
0248
0249       SUBROUTINE ARINI1
0250
0251 c.....before invoking ARINI1:
0252 c.....ARPAR1(1), IAINT2(1) must be set:
0253       PARAMETER (MAXSTR=150001)
0254       double precision  smearp,smearh
0255
0256       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0257 cc      SAVE /ARPRNT/
0258       COMMON /ARPRC/ ITYPAR(MAXSTR),
0259      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
0260      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
0261      &     XMAR(MAXSTR)
0262 cc      SAVE /ARPRC/
0263       COMMON /smearz/smearp,smearh
0264 cc      SAVE /smearz/
0265       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0266 cc      SAVE /input1/
0267       common/anim/nevent,isoft,isflag,izpc
0268 cc      SAVE /anim/
0269       common /nzpc/nattzp
0270 cc      SAVE /nzpc/
0271       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0272 cc      SAVE /HPARNT/
0273       COMMON/RNDF77/NSEED
0274 cc      SAVE /RNDF77/
0275       common /para8/ idpert,npertd,idxsec
0276       SAVE
0277
0278 clin-5/2008 for perturbatively-produced hadrons (currently only deuterons):
0279 c      OPEN (91, FILE = 'ana/deuteron_processes.dat',
0280 c     1     STATUS = 'UNKNOWN')
0281       if(idpert.eq.1.or.idpert.eq.2) then
0282 c         OPEN (90, FILE = 'ana/ampt_pert.dat', STATUS = 'UNKNOWN')
0283       endif
0284 c.....generate formation time and position at formation time.
0285       TAU0 = ARPAR1(1)
0286       NP = IAINT2(1)
0287 clin-7/10/01     initial positions already given for hadrons
0288 c     formed from partons inside ZPC (from string melting):
0289       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
0290          if(NP.le.nattzp) return
0291          do 1001 I = nattzp+1, NP
0292 clin-9/2012 determine rapidity more generally
0293 c     to prevent overflow when Pt~=0 and E=|Pz|:
0294 c            IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
0295 c               PRINT *, ' IN ARINI1'
0296 c               PRINT *, 'ABS(PZ) .GE. EE for particle ', I
0297 c               PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I),
0298 c     &              ' PY = ', PYAR(I)
0299 c               PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
0300 c               PRINT *, ' XM = ', XMAR(I)
0301 c               RAP = 1000000.0
0302 c               GOTO 50
0303 c            END IF
0304 cc            RAP=0.5*LOG((PEAR(I)+PZAR(I))/(PEAR(I)-PZAR(I)))
0305 c            RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
0306 c 50         CONTINUE
0307             if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then
0308                RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
0309             else
0310                PRINT *, ' IN ARINI1 mt=0'
0311                RAP = 1000000.0*sign(1.,PZAR(I))
0312             endif
0313
0314             VX = PXAR(I) / PEAR(I)
0315             VY = PYAR(I) / PEAR(I)
0316             FTAR(I) = TAU0 * COSH(RAP)
0317             GXAR(I) = GXAR(I) + VX * FTAR(I)
0318             GYAR(I) = GYAR(I) + VY * FTAR(I)
0319             GZAR(I) = TAU0 * SINH(RAP)
0320 clin-5/2009 No formation time for spectator projectile or target nucleons:
0321             if(PXAR(I).eq.0.and.PYAR(I).eq.0
0322      2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
0323 clin-2/2013 for spectator target nucleons in LAB frame:
0324 c     1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
0325               if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01)
0326      1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then
0327 c
0328                TAUI=1.E-20
0329                FTAR(I)=TAUI*COSH(RAP)
0330                GZAR(I)=TAUI*SINH(RAP)
0331                endif
0332             endif
0333  1001    continue
0334 clin-7/10/01-end
0335 clin-3/2009 cleanup of program flow:
0336       else
0337          DO 1002 I = 1, NP
0338 clin-9/2012 determine rapidity more generally:
0339 c            IF (ABS(PZAR(I)) .GE. PEAR(I)) THEN
0340 c               PRINT *, ' IN ARINI1'
0341 c               PRINT *, 'ABS(PZ) .GE. EE for particle ', I
0342 c               PRINT *, ' FLAV = ', ITYPAR(I), ' PX = ', PXAR(I),
0343 c     &              ' PY = ', PYAR(I)
0344 c               PRINT *, ' PZ = ', PZAR(I), ' EE = ', PEAR(I)
0345 c               PRINT *, ' XM = ', XMAR(I)
0346 c               RAP = 1000000.0
0347 c               GOTO 100
0348 cc               STOP
0349 c            END IF
0350 c 100        CONTINUE
0351 c            RAP=0.5*LOG((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
0352             if((XMAR(I)**2+PXAR(I)**2+PYAR(I)**2).gt.0.) then
0353                RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
0354             else
0355                PRINT *, ' IN ARINI1 mt=0'
0356                RAP = 1000000.0*sign(1.,PZAR(I))
0357             endif
0358
0359             VX = PXAR(I) / PEAR(I)
0360             VY = PYAR(I) / PEAR(I)
0361 c.....give initial formation time shift
0362             TAUI = FTAR(I) + TAU0
0363             FTAR(I) = TAUI * COSH(RAP)
0364             GXAR(I) = GXAR(I) + VX * TAU0 * COSH(RAP)
0365             GYAR(I) = GYAR(I) + VY * TAU0 * COSH(RAP)
0366 c     4/25/03: hadron z-position upon formation determined the same way as x,y:
0367             GZAR(I) = TAUI * SINH(RAP)
0368 c     the old prescription:
0369 c            GZAR(I) = GZAR(I) + TAU0 * SINH(RAP)
0370             zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
0371             GZAR(I)=GZAR(I)+zsmear
0372 cbz1/28/99end
0373 c     10/05/01 no formation time for spectator projectile or target nucleons:
0374             if(PXAR(I).eq.0.and.PYAR(I).eq.0
0375      2           .and.(ITYPAR(I).eq.2112.or.ITYPAR(I).eq.2212)) then
0376 clin-2/2013 for spectator target nucleons in LAB frame:
0377 c     1           .and.(PEAR(I)*2/HINT1(1)).gt.0.99
0378               if((PEAR(I)/HINT1(6).gt.0.99.and.PEAR(I)/HINT1(6).lt.1.01)
0379      1 .or.(PEAR(I)/HINT1(7).gt.0.99.and.PEAR(I)/HINT1(7).lt.1.01)) then
0380 c
0381 clin-5/2008:
0382 c               TAUI=0.00001
0383                TAUI=1.E-20
0384                FTAR(I)=TAUI*COSH(RAP)
0385                GZAR(I)=TAUI*SINH(RAP)+zsmear
0386                endif
0387             endif
0388  1002    CONTINUE
0389 clin-3/2009 cleanup of program flow:
0390       endif
0391
0392 clin-3/2009 Add initial hadrons before the hadron cascade starts:
0393       call addhad
0394
0395       RETURN
0396       END
0397
0398 c-----------------------------------------------------------------------
0399
0400 c.....subroutine to order particle labels according to increasing
0401 c.....formation time
0402
0403       SUBROUTINE ARTORD
0404
0405 c.....before invoking ARTORD:
0406 c.....IAINT2(1) must be set:
0407       PARAMETER (MAXSTR=150001,MAXR=1)
0408       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0409 cc      SAVE /ARPRNT/
0410       COMMON /ARPRC/ ITYPAR(MAXSTR),
0411      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
0412      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
0413      &     XMAR(MAXSTR)
0414 cc      SAVE /ARPRC/
0415 clin-3/2009 Take care of particle weights when user inserts initial hadrons:
0416       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
0417      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
0418      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
0419       DIMENSION dptemp(MAXSTR)
0420 c
0421       DIMENSION ITYP0(MAXSTR),
0422      &   GX0(MAXSTR), GY0(MAXSTR), GZ0(MAXSTR), FT0(MAXSTR),
0423      &   PX0(MAXSTR), PY0(MAXSTR), PZ0(MAXSTR), EE0(MAXSTR),
0424      &   XM0(MAXSTR)
0425       DIMENSION INDX(MAXSTR)
0426       EXTERNAL ARINDX
0427       SAVE
0428 c
0429       NPAR = 0
0430       NP = IAINT2(1)
0431       DO 1001 I = 1, NP
0432          ITYP0(I) = ITYPAR(I)
0433          GX0(I) = GXAR(I)
0434          GY0(I) = GYAR(I)
0435          GZ0(I) = GZAR(I)
0436          FT0(I) = FTAR(I)
0437          PX0(I) = PXAR(I)
0438          PY0(I) = PYAR(I)
0439          PZ0(I) = PZAR(I)
0440          EE0(I) = PEAR(I)
0441          XM0(I) = XMAR(I)
0442 clin-3/2009:
0443          dptemp(I) = dpertp(I)
0444  1001 CONTINUE
0445       CALL ARINDX(MAXSTR, NP, FT0, INDX)
0446       DO 1002 I = 1, NP
0447 cbz12/3/98
0448 c         IF (ITYP0(INDX(I)) .EQ. 211) THEN
0449 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 321) THEN
0450 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
0451 c     &      ITYP0(INDX(I)) .EQ. 2112 .OR. ITYP0(INDX(I)) .EQ. -211 .OR.
0452 c     &      ITYP0(INDX(I)) .EQ. 111) THEN
0453 c         IF (ITYP0(INDX(I)) .EQ. 211 .OR. ITYP0(INDX(I)) .EQ. 2212 .OR.
0454 c     &      ITYP0(INDX(I)) .EQ. 2112) THEN
0455          NPAR = NPAR + 1
0456 c         ITYPAR(I) = ITYP0(INDX(I))
0457 c         GXAR(I) = GX0(INDX(I))
0458 c         GYAR(I) = GY0(INDX(I))
0459 c         GZAR(I) = GZ0(INDX(I))
0460 c         FTAR(I) = FT0(INDX(I))
0461 c         PXAR(I) = PX0(INDX(I))
0462 c         PYAR(I) = PY0(INDX(I))
0463 c         PZAR(I) = PZ0(INDX(I))
0464 c         PEAR(I) = EE0(INDX(I))
0465 c         XMAR(I) = XM0(INDX(I))
0466          ITYPAR(NPAR) = ITYP0(INDX(I))
0467          GXAR(NPAR) = GX0(INDX(I))
0468          GYAR(NPAR) = GY0(INDX(I))
0469          GZAR(NPAR) = GZ0(INDX(I))
0470          FTAR(NPAR) = FT0(INDX(I))
0471          PXAR(NPAR) = PX0(INDX(I))
0472          PYAR(NPAR) = PY0(INDX(I))
0473          PZAR(NPAR) = PZ0(INDX(I))
0474          PEAR(NPAR) = EE0(INDX(I))
0475          XMAR(NPAR) = XM0(INDX(I))
0476 clin-3/2009:
0477          dpertp(NPAR)=dptemp(INDX(I))
0478 c         END IF
0479 cbz12/3/98end
0480  1002 CONTINUE
0481       IAINT2(1) = NPAR
0482 c
0483       RETURN
0484       END
0485
0486 c-----------------------------------------------------------------------
0487
0488 c.....subroutine to copy individually generated particle record into
0489 c.....particle record for many test particle runs.
0490
0491       SUBROUTINE ARINI2(K)
0492
0493       PARAMETER (MAXSTR=150001,MAXR=1)
0494       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0495 cc      SAVE /ARPRNT/
0496       COMMON /ARPRC/ ITYPAR(MAXSTR),
0497      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
0498      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
0499      &     XMAR(MAXSTR)
0500 cc      SAVE /ARPRC/
0501       COMMON /ARERC1/MULTI1(MAXR)
0502 cc      SAVE /ARERC1/
0503       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
0504      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
0505      &     FT1(MAXSTR, MAXR),
0506      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
0507      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
0508 cc      SAVE /ARPRC1/
0509       COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
0510 cc      SAVE /tdecay/
0511       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0512 cc      SAVE /input1/
0513       COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
0514      &     IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
0515 cc      SAVE /INPUT2/
0516       COMMON/RNDF77/NSEED
0517       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
0518      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
0519      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
0520 cc      SAVE /RNDF77/
0521       SAVE
0522
0523       MULTI1(K) = IAINT2(1)
0524       DO 1001 I = 1, MULTI1(K)
0525          ITYP1(I, K) = ITYPAR(I)
0526          GX1(I, K) = GXAR(I)
0527          GY1(I, K) = GYAR(I)
0528          GZ1(I, K) = GZAR(I)
0529          FT1(I, K) = FTAR(I)
0530          PX1(I, K) = PXAR(I)
0531          PY1(I, K) = PYAR(I)
0532          PZ1(I, K) = PZAR(I)
0533          EE1(I, K) = PEAR(I)
0534          XM1(I, K) = XMAR(I)
0535 clin-3/2009 hadron weights are initialized in addhad():
0536 clin-5/2008 all hadrons not perturbatively-produced have the weight of 1:
0537 c         dpp1(I,K)=1.
0538          dpp1(I,K)=dpertp(I)
0539  1001 CONTINUE
0540
0541 c     initialize final time of each particle to ntmax*dt except for
0542 c     decay daughters, which have values given by tfdcy() and >(ntmax*dt):
0543       do 1002 ip=1,MAXSTR
0544          tfdcy(ip)=NTMAX*DT
0545          tft(ip)=NTMAX*DT
0546  1002 continue
0547 c
0548       do 1004 irun=1,MAXR
0549          do 1003 ip=1,MAXSTR
0550             tfdpi(ip,irun)=NTMAX*DT
0551  1003    continue
0552  1004 continue
0553
0554       RETURN
0555       END
0556
0557 c=======================================================================
0558
0559 c.....function to convert PDG flavor code into ART flavor code.
0560
0561       FUNCTION IARFLV(IPDG)
0562
0563       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0564 cc      SAVE /input1/
0565       COMMON/RNDF77/NSEED
0566 cc      SAVE /RNDF77/
0567       SAVE
0568
0569 c.....anti-Delta-
0570       IF (IPDG .EQ. -1114) THEN
0571          IARFLV = -6
0572          RETURN
0573       END IF
0574
0575 c.....anti-Delta0
0576       IF (IPDG .EQ. -2114) THEN
0577          IARFLV = -7
0578          RETURN
0579       END IF
0580
0581 c.....anti-Delta+
0582       IF (IPDG .EQ. -2214) THEN
0583          IARFLV = -8
0584          RETURN
0585       END IF
0586
0587 c.....anti-Delta++
0588       IF (IPDG .EQ. -2224) THEN
0589          IARFLV = -9
0590          RETURN
0591       END IF
0592
0593 cbzdbg2/23/99
0594 c.....anti-proton
0595       IF (IPDG .EQ. -2212) THEN
0596          IARFLV = -1
0597          RETURN
0598       END IF
0599
0600 c.....anti-neutron
0601       IF (IPDG .EQ. -2112) THEN
0602          IARFLV = -2
0603          RETURN
0604       END IF
0605 cbzdbg2/23/99end
0606
0607 c.....eta
0608       IF (IPDG .EQ. 221) THEN
0609          IARFLV = 0
0610          RETURN
0611       END IF
0612
0613 c.....proton
0614       IF (IPDG .EQ. 2212) THEN
0615          IARFLV = 1
0616          RETURN
0617       END IF
0618
0619 c.....neutron
0620       IF (IPDG .EQ. 2112) THEN
0621          IARFLV = 2
0622          RETURN
0623       END IF
0624
0625 c.....pi-
0626       IF (IPDG .EQ. -211) THEN
0627          IARFLV = 3
0628          RETURN
0629       END IF
0630
0631 c.....pi0
0632       IF (IPDG .EQ. 111) THEN
0633          IARFLV = 4
0634          RETURN
0635       END IF
0636
0637 c.....pi+
0638       IF (IPDG .EQ. 211) THEN
0639          IARFLV = 5
0640          RETURN
0641       END IF
0642
0643 c.....Delta-
0644       IF (IPDG .EQ. 1114) THEN
0645          IARFLV = 6
0646          RETURN
0647       END IF
0648
0649 c.....Delta0
0650       IF (IPDG .EQ. 2114) THEN
0651          IARFLV = 7
0652          RETURN
0653       END IF
0654
0655 c.....Delta+
0656       IF (IPDG .EQ. 2214) THEN
0657          IARFLV = 8
0658          RETURN
0659       END IF
0660
0661 c.....Delta++
0662       IF (IPDG .EQ. 2224) THEN
0663          IARFLV = 9
0664          RETURN
0665       END IF
0666
0667 c.....Lambda
0668       IF (IPDG .EQ. 3122) THEN
0669          IARFLV = 14
0670          RETURN
0671       END IF
0672
0673 c.....Lambda-bar
0674       IF (IPDG .EQ. -3122) THEN
0675          IARFLV = -14
0676          RETURN
0677       END IF
0678
0679 c.....Sigma-
0680       IF (IPDG .EQ. 3112) THEN
0681          IARFLV = 15
0682          RETURN
0683       END IF
0684
0685 c.....Sigma-bar
0686       IF (IPDG .EQ. -3112) THEN
0687          IARFLV = -15
0688          RETURN
0689       END IF
0690
0691 c.....Sigma0
0692       IF (IPDG .EQ. 3212) THEN
0693          IARFLV = 16
0694          RETURN
0695       END IF
0696
0697 c.....Sigma0-bar
0698       IF (IPDG .EQ. -3212) THEN
0699          IARFLV = -16
0700          RETURN
0701       END IF
0702
0703 c.....Sigma+
0704       IF (IPDG .EQ. 3222) THEN
0705          IARFLV = 17
0706          RETURN
0707       END IF
0708
0709 c.....Sigma+ -bar
0710       IF (IPDG .EQ. -3222) THEN
0711          IARFLV = -17
0712          RETURN
0713       END IF
0714
0715 c.....K-
0716       IF (IPDG .EQ. -321) THEN
0717          IARFLV = 21
0718          RETURN
0719       END IF
0720
0721 c.....K+
0722       IF (IPDG .EQ. 321) THEN
0723          IARFLV = 23
0724          RETURN
0725       END IF
0726
0727 c.....temporary entry for K0
0728       IF (IPDG .EQ. 311) THEN
0729          IARFLV = 23
0730          RETURN
0731       END IF
0732
0733 c.....temporary entry for K0bar
0734       IF (IPDG .EQ. -311) THEN
0735          IARFLV = 21
0736          RETURN
0737       END IF
0738
0739 c.....temporary entry for K0S and K0L
0740       IF (IPDG .EQ. 310 .OR. IPDG .EQ. 130) THEN
0741          R = RANART(NSEED)
0742          IF (R .GT. 0.5) THEN
0743             IARFLV = 23
0744          ELSE
0745             IARFLV = 21
0746          END IF
0747          RETURN
0748       END IF
0749
0750 c.....rho-
0751       IF (IPDG .EQ. -213) THEN
0752          IARFLV = 25
0753          RETURN
0754       END IF
0755
0756 c.....rho0
0757       IF (IPDG .EQ. 113) THEN
0758          IARFLV = 26
0759          RETURN
0760       END IF
0761
0762 c.....rho+
0763       IF (IPDG .EQ. 213) THEN
0764          IARFLV = 27
0765          RETURN
0766       END IF
0767
0768 c.....omega
0769       IF (IPDG .EQ. 223) THEN
0770          IARFLV = 28
0771          RETURN
0772       END IF
0773
0774 c.....phi
0775       IF (IPDG .EQ. 333) THEN
0776          IARFLV = 29
0777          RETURN
0778       END IF
0779
0780 c.....K*+
0781       IF (IPDG .EQ. 323) THEN
0782          IARFLV = 30
0783          RETURN
0784       END IF
0785 c.....K*-
0786       IF (IPDG .EQ. -323) THEN
0787          IARFLV = -30
0788          RETURN
0789       END IF
0790 c.....temporary entry for K*0
0791       IF (IPDG .EQ. 313) THEN
0792          IARFLV = 30
0793          RETURN
0794       END IF
0795 c.....temporary entry for K*0bar
0796       IF (IPDG .EQ. -313) THEN
0797          IARFLV = -30
0798          RETURN
0799       END IF
0800
0801 c...... eta-prime
0802       IF (IPDG .EQ. 331) THEN
0803          IARFLV = 31
0804          RETURN
0805       END IF
0806
0807 c...... a1
0808 c     IF (IPDG .EQ. 777) THEN
0809 c        IARFLV = 32
0810 c        RETURN
0811 c     END IF
0812
0813 c... cascade-
0814       IF (IPDG .EQ. 3312) THEN
0815          IARFLV = 40
0816          RETURN
0817       END IF
0818
0819 c... cascade+ (bar)
0820       IF (IPDG .EQ. -3312) THEN
0821          IARFLV = -40
0822          RETURN
0823       END IF
0824
0825 c... cascade0
0826       IF (IPDG .EQ. 3322) THEN
0827          IARFLV = 41
0828          RETURN
0829       END IF
0830
0831 c... cascade0 -bar
0832       IF (IPDG .EQ. -3322) THEN
0833          IARFLV = -41
0834          RETURN
0835       END IF
0836
0837 c... Omega-
0838       IF (IPDG .EQ. 3334) THEN
0839          IARFLV = 45
0840          RETURN
0841       END IF
0842
0843 c... Omega+ (bar)
0844       IF (IPDG .EQ. -3334) THEN
0845          IARFLV = -45
0846          RETURN
0847       END IF
0848
0849 c... Di-Omega
0850       IF (IPDG .EQ. 6666) THEN
0851          IARFLV = 44
0852          RETURN
0853       END IF
0854 c sp06/05/01 end
0855
0856 clin-3/2009 keep the same ID numbers in case there are initial deuterons:
0857       IF (IPDG .EQ. 42 .or. IPDG .EQ. -42) THEN
0858          IARFLV = IPDG
0859          RETURN
0860       END IF
0861
0862 c.....other
0863       IARFLV = IPDG + 10000
0864
0865       RETURN
0866       END
0867
0868 c-----------------------------------------------------------------------
0869
0870 c.....function to convert ART flavor code into PDG flavor code.
0871
0872       FUNCTION INVFLV(IART)
0873
0874       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0875 cc      SAVE /input1/
0876       COMMON/RNDF77/NSEED
0877 cc      SAVE /RNDF77/
0878       SAVE
0879
0880 c.....anti-Delta-
0881       IF (IART .EQ. -6) THEN
0882          INVFLV = -1114
0883          RETURN
0884       END IF
0885
0886 c.....anti-Delta0
0887       IF (IART .EQ. -7) THEN
0888          INVFLV = -2114
0889          RETURN
0890       END IF
0891
0892 c.....anti-Delta+
0893       IF (IART .EQ. -8) THEN
0894          INVFLV = -2214
0895          RETURN
0896       END IF
0897
0898 c.....anti-Delta++
0899       IF (IART .EQ. -9) THEN
0900          INVFLV = -2224
0901          RETURN
0902       END IF
0903
0904 cbzdbg2/23/99
0905 c.....anti-proton
0906       IF (IART .EQ. -1) THEN
0907          INVFLV = -2212
0908          RETURN
0909       END IF
0910
0911 c.....anti-neutron
0912       IF (IART .EQ. -2) THEN
0913          INVFLV = -2112
0914          RETURN
0915       END IF
0916 cbzdbg2/23/99end
0917
0918 c.....eta
0919       IF (IART .EQ. 0) THEN
0920          INVFLV = 221
0921          RETURN
0922       END IF
0923
0924 c.....proton
0925       IF (IART .EQ. 1) THEN
0926          INVFLV = 2212
0927          RETURN
0928       END IF
0929
0930 c.....neutron
0931       IF (IART .EQ. 2) THEN
0932          INVFLV = 2112
0933          RETURN
0934       END IF
0935
0936 c.....pi-
0937       IF (IART .EQ. 3) THEN
0938          INVFLV = -211
0939          RETURN
0940       END IF
0941
0942 c.....pi0
0943       IF (IART .EQ. 4) THEN
0944          INVFLV = 111
0945          RETURN
0946       END IF
0947
0948 c.....pi+
0949       IF (IART .EQ. 5) THEN
0950          INVFLV = 211
0951          RETURN
0952       END IF
0953
0954 c.....Delta-
0955       IF (IART .EQ. 6) THEN
0956          INVFLV = 1114
0957          RETURN
0958       END IF
0959
0960 c.....Delta0
0961       IF (IART .EQ. 7) THEN
0962          INVFLV = 2114
0963          RETURN
0964       END IF
0965
0966 c.....Delta+
0967       IF (IART .EQ. 8) THEN
0968          INVFLV = 2214
0969          RETURN
0970       END IF
0971
0972 c.....Delta++
0973       IF (IART .EQ. 9) THEN
0974          INVFLV = 2224
0975          RETURN
0976       END IF
0977
0978 cc.....N*(1440), N*(1535) temporary entry
0979 c      IF (IART .GE. 10 .AND. IART .LE.13) THEN
0980 c         INVFLV = 0
0981 c         RETURN
0982 c      END IF
0983
0984 c.....Lambda
0985       IF (IART .EQ. 14) THEN
0986          INVFLV = 3122
0987          RETURN
0988       END IF
0989 c.....Lambda-bar
0990       IF (IART .EQ. -14) THEN
0991          INVFLV = -3122
0992          RETURN
0993       END IF
0994
0995 cbz3/12/99
0996 c.....temporary entry for Sigma's
0997 c      IF (IART .EQ. 15) THEN
0998 c         R = RANART(NSEED)
0999 c         IF (R .GT. 2. / 3.) THEN
1000 c            INVFLV = 3112
1001 c         ELSE IF (R .GT. 1./ 3. .AND. R .LE. 2. / 3.) THEN
1002 c            INVFLV = 3212
1003 c         ELSE
1004 c            INVFLV = 3222
1005 c         END IF
1006 c         RETURN
1007 c      END IF
1008
1009 c.....Sigma-
1010       IF (IART .EQ. 15) THEN
1011          INVFLV = 3112
1012          RETURN
1013       END IF
1014
1015 c.....Sigma- bar
1016       IF (IART .EQ. -15) THEN
1017          INVFLV = -3112
1018          RETURN
1019       END IF
1020
1021 c.....Sigma0
1022       IF (IART .EQ. 16) THEN
1023          INVFLV = 3212
1024          RETURN
1025       END IF
1026
1027 c.....Sigma0 -bar
1028       IF (IART .EQ. -16) THEN
1029          INVFLV = -3212
1030          RETURN
1031       END IF
1032
1033 c.....Sigma+
1034       IF (IART .EQ. 17) THEN
1035          INVFLV = 3222
1036          RETURN
1037       END IF
1038
1039 c.....Sigma+ -bar
1040       IF (IART .EQ. -17) THEN
1041          INVFLV = -3222
1042          RETURN
1043       END IF
1044
1045 clin-2/23/03 K0S and K0L are generated at the last timestep:
1046 c.....temporary entry for K- and K0bar
1047       IF (IART .EQ. 21) THEN
1048 c         R = RANART(NSEED)
1049 c         IF (R .GT. 0.5) THEN
1050             INVFLV = -321
1051 c         ELSE
1052 c            INVFLV = -311
1053 c            R = RANART(NSEED)
1054 c            IF (R .GT. 0.5) THEN
1055 c               INVFLV = 310
1056 c            ELSE
1057 c               INVFLV = 130
1058 c            END IF
1059 c         END IF
1060          RETURN
1061       END IF
1062
1063 c.....temporary entry for K+ and K0
1064       IF (IART .EQ. 23) THEN
1065 c         R = RANART(NSEED)
1066 c         IF (R .GT. 0.5) THEN
1067             INVFLV = 321
1068 c         ELSE
1069 c            INVFLV = 311
1070 c            R = RANART(NSEED)
1071 c            IF (R .GT. 0.5) THEN
1072 c               INVFLV = 310
1073 c            ELSE
1074 c               INVFLV = 130
1075 c            END IF
1076 c         END IF
1077          RETURN
1078       END IF
1079
1080 c.....K0Long:
1081       IF (IART .EQ. 22) THEN
1082          INVFLV = 130
1083          RETURN
1084       ENDIF
1085 c.....K0Short:
1086       IF (IART .EQ. 24) THEN
1087          INVFLV = 310
1088          RETURN
1089       ENDIF
1090
1091 c.....rho-
1092       IF (IART .EQ. 25) THEN
1093          INVFLV = -213
1094          RETURN
1095       END IF
1096
1097 c.....rho0
1098       IF (IART .EQ. 26) THEN
1099          INVFLV = 113
1100          RETURN
1101       END IF
1102
1103 c.....rho+
1104       IF (IART .EQ. 27) THEN
1105          INVFLV = 213
1106          RETURN
1107       END IF
1108
1109 c.....omega
1110       IF (IART .EQ. 28) THEN
1111          INVFLV = 223
1112          RETURN
1113       END IF
1114
1115 c.....phi
1116       IF (IART .EQ. 29) THEN
1117          INVFLV = 333
1118          RETURN
1119       END IF
1120
1121 c.....temporary entry for K*+ and K*0
1122       IF (IART .EQ. 30) THEN
1123          INVFLV = 323
1124          IF (RANART(NSEED).GT.0.5) INVFLV = 313
1125          RETURN
1126       END IF
1127
1128 c.....temporary entry for K*- and K*0bar
1129       IF (IART .EQ. -30) THEN
1130          INVFLV = -323
1131          IF (RANART(NSEED).GT.0.5) INVFLV = -313
1132          RETURN
1133       END IF
1134
1135 c... eta-prime (bar)
1136       IF (IART .EQ. 31) THEN
1137          INVFLV = 331
1138          RETURN
1139       END IF
1140
1141 c... a1
1142       IF (IART .EQ. 32) THEN
1143          INVFLV = 777
1144          RETURN
1145       END IF
1146
1147 c... cascade-
1148       IF (IART .EQ. 40) THEN
1149          INVFLV = 3312
1150          RETURN
1151       END IF
1152
1153 c... cascade+ (bar)
1154       IF (IART .EQ. -40) THEN
1155          INVFLV = -3312
1156          RETURN
1157       END IF
1158
1159 c... cascade0
1160       IF (IART .EQ. 41) THEN
1161          INVFLV = 3322
1162          RETURN
1163       END IF
1164
1165 c... cascade0 -bar
1166       IF (IART .EQ. -41) THEN
1167          INVFLV = -3322
1168          RETURN
1169       END IF
1170
1171 c... Omega-
1172       IF (IART .EQ. 45) THEN
1173          INVFLV = 3334
1174          RETURN
1175       END IF
1176
1177 c... Omega+ (bar)
1178       IF (IART .EQ. -45) THEN
1179          INVFLV = -3334
1180          RETURN
1181       END IF
1182
1183 c... Di-Omega
1184       IF (IART .EQ. 44) THEN
1185          INVFLV = 6666
1186          RETURN
1187       END IF
1188 c sp 12/19/00 end
1189
1190 clin-5/2008 deuteron ID numbers in ART and ampt.dat:
1191       IF (IART .EQ. 42) THEN
1192          INVFLV = 42
1193          RETURN
1194       ELSEIF (IART .EQ. -42) THEN
1195          INVFLV = -42
1196          RETURN
1197       END IF
1198 c
1199 c.....other
1200       INVFLV = IART - 10000
1201
1202       RETURN
1203       END
1204
1205 c=======================================================================
1206
1207       BLOCK DATA ARDATA
1208
1209       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
1210 cc      SAVE /ARPRNT/
1211       SAVE
1212       DATA ARPAR1/1.19, 99 * 0.0/
1213       DATA IAPAR2/3, 49 * 0/
1214       DATA ARINT1/100 * 0.0/
1215       DATA IAINT2/50 * 0/
1216
1217       END
1218
1219 c=======================================================================
1220
1221 c.....Routine borrowed from ZPC.
1222 c.....double precision  is modified to real*4.
1223
1224 cbz1/29/99
1225 c      subroutine index1(n, m, arrin, indx)
1226       subroutine arindx(n, m, arrin, indx)
1227 cbz1/29/99end
1228 c     indexes the first m elements of ARRIN of length n, i.e., outputs INDX
1229 c     such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
1230
1231 c      implicit real*4 (a-h, o-z)
1232
1233       dimension arrin(n), indx(n)
1234       SAVE
1235       do 1001 j = 1, m
1236          indx(j) = j
1237  1001 continue
1238       l = m / 2 + 1
1239       ir = m
1240  10   continue
1241       if (l .gt. 1) then
1242          l = l - 1
1243          indxt = indx(l)
1244          q = arrin(indxt)
1245       else
1246          indxt = indx(ir)
1247          q = arrin(indxt)
1248          indx(ir) = indx(1)
1249          ir = ir - 1
1250          if (ir .eq. 1) then
1251             indx(1) = indxt
1252             return
1253          end if
1254       end if
1255       i = l
1256       j = l + l
1257  20   if (j .le. ir) then
1258          if (j .lt. ir) then
1259             if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
1260          end if
1261          if (q .lt. arrin(indx(j))) then
1262             indx(i) = indx(j)
1263             i = j
1264             j = j + j
1265          else
1266             j = ir + 1
1267          end if
1268       goto 20
1269       end if
1270       indx(i) = indxt
1271       goto 10
1272
1273       end
1274
1275 c-----------------------------------------------------------------------
1276
1277 c.....extracted from G. Song's ART expasion including K- interactions
1278 c.....file `NEWKAON.FOR'
1279
1280 c     5/01/03 send iblock value into art1f.f, necessary for resonance studies:
1281 c        subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1282 c     &                                   srt,pcx,pcy,pcz)
1283         subroutine newka(icase,irun,iseed,dt,nt,ictrl,i1,i2,
1284      &                                   srt,pcx,pcy,pcz,iblock)
1285       PARAMETER      (MAXSTR=150001,MAXR=1)
1286       PARAMETER      (AKA=0.498)
1287       COMMON   /AA/  R(3,MAXSTR)
1288 cc      SAVE /AA/
1289       COMMON   /BB/  P(3,MAXSTR)
1290 cc      SAVE /BB/
1291       COMMON   /CC/  E(MAXSTR)
1292 cc      SAVE /CC/
1293       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
1294 cc      SAVE /EE/
1295       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
1296 cc      SAVE /BG/
1297       COMMON   /NN/NNN
1298 cc      SAVE /NN/
1299       COMMON   /RUN/NUM
1300 cc      SAVE /RUN/
1301       COMMON   /PA/RPION(3,MAXSTR,MAXR)
1302 cc      SAVE /PA/
1303       COMMON   /PB/PPION(3,MAXSTR,MAXR)
1304 cc      SAVE /PB/
1305       COMMON   /PC/EPION(MAXSTR,MAXR)
1306 cc      SAVE /PC/
1307       COMMON   /PD/LPION(MAXSTR,MAXR)
1308 cc      SAVE /PD/
1309       COMMON/RNDF77/NSEED
1310 cc      SAVE /RNDF77/
1311       SAVE
1312 c
1313         logical lb1bn, lb2bn,lb1mn,lb2mn
1314 cbz3/7/99 neutralk
1315 c        logical lb1bn1, lb2bayon1, lb1bn0, lb2bn0
1316         logical lb1bn1, lb2bn1, lb1bn0, lb2bn0
1317 cbz3/7/99 neutralk end
1318         logical lb1mn0, lb2mn0, lb1mn1, lb2mn1
1319         logical lb1mn2, lb2mn2
1320         icase=-1
1321 c        icase: flag for the type of reaction that is going to happen.
1322 c        icase=-1,  no desired reaction, return to main program.
1323 c              1,  NN,ND,DD
1324 c              2,  PI+N, PI+D
1325 c              3,  K(-) absorption.
1326         nchrg=-100
1327 c        nchrg: Net charges of the two incoming particles.
1328         ictrl = 1
1329         lb1=lb(i1)
1330         lb2=lb(i2)
1331         em1=e(i1)
1332         em2=e(i2)
1333         lb1bn=lb1.eq.1.or.lb1.eq.2.or.(lb1.gt.5.and.lb1.le.13)
1334         lb2bn=lb2.eq.1.or.lb2.eq.2.or.(lb2.gt.5.and.lb2.le.13)
1335         lb1bn0=lb1.eq.2.or.lb1.eq.7.or.lb1.eq.10.or.lb1.eq.12
1336         lb2bn0=lb2.eq.2.or.lb2.eq.7.or.lb2.eq.10.or.lb2.eq.12
1337         lb1bn1=lb1.eq.1.or.lb1.eq.8.or.lb1.eq.11.or.lb1.eq.13
1338         lb2bn1=lb2.eq.1.or.lb2.eq.8.or.lb2.eq.11.or.lb2.eq.13
1339         lb1mn=em1.lt.0.2.or.lb1.eq.0.or.(lb1.ge.25.and.lb1.le.29)
1340         lb2mn=em2.lt.0.2.or.lb2.eq.0.or.(lb2.ge.25.and.lb2.le.29)
1341         lb1mn0=lb1.eq.0.or.lb1.eq.4.or.lb1.eq.26.or.
1342      &                        lb1.eq.28.or.lb1.eq.29
1343         lb2mn0=lb2.eq.0.or.lb2.eq.4.or.lb2.eq.26.or.
1344      &                        lb2.eq.28.or.lb2.eq.29
1345         lb1mn1= lb1.eq.5.or.lb1.eq.27
1346         lb2mn1= lb2.eq.5.or.lb2.eq.27
1347         lb1mn2=lb1.eq.3.or.lb1.eq.25
1348         lb2mn2=lb2.eq.3.or.lb2.eq.25
1349
1350 c        1. consider N+N, N+Resonance, R + R reactions
1351         if(lb1bn.and.lb2bn) then
1352 c     NN,ND,DD:
1353            icase=1
1354 c     total cross section
1355            sig=40.
1356            if(lb1.eq.9.and.lb2.eq.9) then
1357                 nchrg=4
1358            endif
1359            if((lb1bn1.and.lb2.eq.9)
1360      &        .or.(lb2bn1.and.lb1.eq.9))then
1361                 nchrg=3
1362            endif
1363            if((lb1bn0.and.lb2.eq.9)
1364      &        .or.(lb2bn0.and.lb1.eq.9)
1365      &        .or.(lb1bn1.and.lb2bn1)) then
1366                    nchrg=2
1367            endif
1368            if((lb1bn1.and.lb2bn0).or.(lb1.eq.6.and.lb2.eq.9)
1369      &        .or.(lb2bn1.and.lb1bn0)
1370      &        .or.(lb2.eq.6.and.lb1.eq.9))then
1371                    nchrg=1
1372            endif
1373            if((lb1bn0.and.lb2bn0).or.(lb1bn1.and.lb2.eq.6)
1374      &              .or.(lb2bn1.and.lb1.eq.6)) then
1375                    nchrg=0
1376            endif
1377            if((lb1bn0.and.lb2.eq.6)
1378      &        .or.(lb2bn0.and.lb1.eq.6))then
1379                 nchrg=-1
1380            endif
1381            if(lb1.eq.6.and.lb2.eq.6) then
1382                 nchrg=-2
1383            endif
1384 c     brsig = x2kaon_no_isospin(srt)
1385            if(nchrg.ge.-1.and.nchrg.le.2) then
1386 c     K,Kbar prduction x sect.
1387                    brsig = x2kaon(srt)
1388            else
1389                    brsig=0.0
1390 c                if(nchrg.eq.-2.or.nchrg.eq.3) then
1391 c                   brsig = x2kaon(srt+0.938-1.232)
1392 c                else
1393 c     nchrg=4
1394 c                   brsig = x2kaon(srt+2.*(0.938-1.232))
1395 c                endif
1396            endif
1397
1398 cbz3/7/99 neutralk
1399            BRSIG = 2.0 * BRSIG
1400 cbz3/7/99 neutralk end
1401
1402         endif
1403
1404 c        2. consider PI(meson:eta,omega,rho,phi) + N(N*,D)
1405         if((lb1bn.and.lb2mn).OR.(lb2bn.and.lb1mn)) then
1406 c     PN,PD
1407           icase=2
1408           sig=20.
1409           sigma0 = piNsg0(srt)
1410           brsig=0.0
1411           if((lb1bn1.and.lb2mn0)
1412      &       .or.(lb2bn1.and.lb1mn0).
1413      & or.(lb1bn0.and.lb2mn1).or.(lb2bn0.and.lb1mn1).
1414      & or.(lb1.eq.9.and.lb2mn2).or.(lb2.eq.9.and.lb1mn2))then
1415                 nchrg=1
1416 cbz3/2/99/song
1417 c                if(lb1bn1.or.lb2bn1) brsig=2.0*sigma0
1418 c                if(lb1bn0.or.lb2bn0) brsig=0.5*sigma0
1419                 if(lb1bn1.or.lb2bn1) brsig=0.5*sigma0
1420                 if(lb1bn0.or.lb2bn0) brsig=2.0*sigma0
1421 cbz3/2/99/song end
1422 c                if(lb1.eq.9.or.lb2.eq.9) brsig=1.5*sigma0
1423           endif
1424           if( (lb1bn0.and.lb2mn0 )
1425      &       .or.(lb2bn0.and.lb1mn0)
1426      &  .or.(lb1bn1.and.lb2mn2).or.(lb2bn1.and.lb1mn2)
1427      &  .or.(lb1.eq.6.and.lb2mn1).or.(lb2.eq.6.and.lb1mn1)) then
1428                 nchrg=0
1429                 if(lb1bn1.or.lb2bn1) then
1430 cbz3/2/99/song
1431 c                  brsig=1.5*sigma0
1432                   brsig=3.0*sigma0
1433 cbz3/2/99/song end
1434 cbz3/11/99/song
1435 c                  ratiok = 1./3.
1436                   ratiok = 2./3.
1437 cbz3/11/99/song end
1438
1439 c                  ratiok: the ratio of channels: ->nK+k- vs. -> pK0K-
1440                 endif
1441                 if(lb1bn0.or.lb2bn0) then
1442                   brsig=2.5*sigma0
1443 cbz3/2/99/song
1444 c                  ratiok = 0.8
1445                   ratiok = 0.2
1446 cbz3/2/99/song end
1447                 endif
1448 c                if(lb1.eq.6.or.lb2.eq.6) then
1449 c     lb=6 : D-
1450 c                  brsig=1.5*sigma0
1451 c                  ratiok = 0.5
1452 c                endif
1453           endif
1454           if( (lb1bn0.and.lb2mn2)
1455      &       .or.(lb2bn0.and.lb1mn2)
1456      & .or.(lb1.eq.6.and.lb2mn0).or.(lb2.eq.6.and.lb1mn0)) then
1457                 nchrg=-1
1458                 if(lb1bn0.or.lb2bn0) brsig=sigma0
1459 c                if(lb1.eq.6.or.lb2.eq.6) brsig=sigma0
1460           endif
1461 c          if((lb1.eq.6.and.lb2mn2).or.(lb2.eq.6.and.lb1mn2))then
1462 c                nchrg=-2
1463 c          endif
1464 c          if((lb1bn1.and.lb2mn1).or.(lb2bn1.and.lb1mn1)
1465 c    &           .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1466 c                nchrg=2
1467 c          endif
1468
1469 cbz3/11/99 neutralk
1470           if((lb1.eq.6.and.lb2mn2)
1471      &       .or.(lb2.eq.6.and.lb1mn2))then
1472                 nchrg=-2
1473           endif
1474 cbz3/11/99 neutralk
1475 cbz3/8/99 neutralk
1476           if((lb1bn1.and.lb2mn1)
1477      &       .or.(lb2bn1.and.lb1mn1)
1478      & .or.(lb1.eq.9.and.lb2mn0).or.(lb2.eq.9.and.lb1mn0)) then
1479                 nchrg=2
1480           endif
1481 cbz3/8/99 neutralk end
1482
1483 cbz3/7/99 neutralk
1484           IF (NCHRG .GE. -2 .AND. NCHRG .LE. 2) THEN
1485              BRSIG = 3.0 * SIGMA0
1486           END IF
1487 cbz3/7/99 neutralk end
1488
1489         endif
1490
1491 c        3. consider K- + N(N*,D) absorption.
1492 c        if((lb1bn.and.lb2.eq.21).OR.(lb2bn.and.lb1.eq.21)) then
1493         if( (lb1bn.and.(lb2.eq.21.or.lb2.eq.-30)).OR.
1494      &     (lb2bn.and.(lb1.eq.21.or.lb1.eq.-30)) )then
1495 c          bmass=em1+em2-aka
1496           bmass=0.938
1497           if(srt.le.(bmass+aka)) then
1498 cbz3/2/99
1499 c                write(100,*)'--lb1,lb2,em1,em2,srt',lb1,lb2,em1,em2,srt
1500 cbz3/2/99end
1501                 pkaon=0.
1502           else
1503             pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
1504           endif
1505           sig=0.
1506           if(lb1.eq.1.or.lb2.eq.1.or.lb1.eq.8.or.lb2.eq.8.or.
1507      &    lb1.eq.11.or.lb2.eq.11.or.lb1.eq.13.or.lb2.eq.13) then
1508 c          K- + (D+,N*+)p ->
1509               nchrg=0
1510               sigela=akPel(pkaon)
1511               sigsgm=3.*akPsgm(pkaon)
1512               sig=sigela+sigsgm+akPlam(pkaon)
1513           endif
1514           if(lb1.eq.2.or.lb2.eq.2.or.lb1.eq.7.or.lb2.eq.7.or.
1515      &    lb1.eq.10.or.lb2.eq.10.or.lb1.eq.12.or.lb2.eq.12) then
1516 c          K- + (D0, N*0)n ->
1517               nchrg=-1
1518               sigela=akNel(pkaon)
1519               sigsgm=2.*akNsgm(pkaon)
1520               sig=sigela+sigsgm+akNlam(pkaon)
1521           endif
1522           if(lb1.eq.6.or.lb2.eq.6) then
1523 c     K- + D-
1524               nchrg=-2
1525               sigela=akNel(pkaon)
1526               sigsgm=akNsgm(pkaon)
1527               sig=sigela+sigsgm
1528           endif
1529           if(lb1.eq.9.or.lb2.eq.9) then
1530 c     K- + D++
1531               nchrg=1
1532               sigela=akPel(pkaon)
1533               sigsgm=2.*akPsgm(pkaon)
1534               sig=sigela+sigsgm+akPlam(pkaon)
1535           endif
1536
1537 cbz3/8/99 neutralk
1538           sigela = 0.5 * (AKPEL(PKAON) + AKNEL(PKAON))
1539           SIGSGM = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1540           SIG = sigela + SIGSGM + AKPLAM(PKAON)
1541 cbz3/8/99 neutralk end
1542
1543           if(sig.gt.1.e-7) then
1544 c     K(-) + N reactions
1545               icase=3
1546               brel=sigela/sig
1547               brsgm=sigsgm/sig
1548 c              branch_lambda=akNlam(pkaon)/sig
1549               brsig = sig
1550           endif
1551         endif
1552
1553 c        4. meson + hyperon -> K- + N
1554 c        if(((lb1.ge.14.and.lb1.le.17).and.lb2mn).OR.
1555 c     &     ((lb2.ge.14.and.lb2.le.17).and.lb1mn)) then
1556         if(((lb1.ge.14.and.lb1.le.17).and.(lb2.ge.3.and.lb2.le.5)).OR.
1557      &     ((lb2.ge.14.and.lb2.le.17).and.(lb1.ge.3.and.lb1.le.5)))then
1558 c        first classify the reactions due to total charge.
1559            nchrg=-100
1560            if((lb1.eq.15.and.(lb2.eq.3.or.lb2.eq.25)).OR.
1561      &              (lb2.eq.15.and.(lb1.eq.3.or.lb1.eq.25))) then
1562                 nchrg=-2
1563 c     D-
1564                   bmass=1.232
1565            endif
1566            if((lb1.eq.15.and.lb2mn0).or.(lb2.eq.15.and.lb1mn0).OR.
1567      &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.3.or.lb2.eq.25)).OR.
1568      &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.3.or.lb1.eq.25)))then
1569                 nchrg=-1
1570 c     n
1571                  bmass=0.938
1572            endif
1573            if((lb1.eq.15.and.(lb2.eq.5.or.lb2.eq.27)).OR.
1574      &              (lb2.eq.15.and.(lb1.eq.5.or.lb1.eq.27)).or.
1575      &        (lb1.eq.17.and.(lb2.eq.3.or.lb2.eq.25)).OR.
1576      &              (lb2.eq.17.and.(lb1.eq.3.or.lb1.eq.25)).or.
1577      &       ((lb1.eq.14.or.lb1.eq.16).and.lb2mn0).OR.
1578      &       ((lb2.eq.14.or.lb2.eq.16).and.lb1mn0)) then
1579                 nchrg=0
1580 c     p
1581                  bmass=0.938
1582            endif
1583            if((lb1.eq.17.and.lb2mn0).or.(lb2.eq.17.and.lb1mn0).OR.
1584      &       ((lb1.eq.14.or.lb1.eq.16).and.(lb2.eq.5.or.lb2.eq.27)).OR.
1585      &       ((lb2.eq.14.or.lb2.eq.16).and.(lb1.eq.5.or.lb1.eq.27)))then
1586                 nchrg=1
1587 c     D++
1588                  bmass=1.232
1589            endif
1590            sig = 0.
1591            if(nchrg.ne.-100.and.srt.gt.(aka+bmass)) then
1592 c     PI+sigma or PI + Lambda => Kbar + N reactions
1593              icase=4
1594 c             pkaon=sqrt(((srt**2-(aka**2+bmass**2))/2./bmass)**2-aka**2)
1595              pkaon=sqrt(((srt**2-(aka**2+0.938**2))/2./0.938)**2-aka**2)
1596 c     lambda + Pi
1597              if(lb1.eq.14.or.lb2.eq.14) then
1598                 if(nchrg.ge.0) sigma0=akPlam(pkaon)
1599                 if(nchrg.lt.0) sigma0=akNlam(pkaon)
1600 c     sigma + pi
1601              else
1602 c     K-p or K-D++
1603                 if(nchrg.ge.0) sigma0=akPsgm(pkaon)
1604 c     K-n or K-D-
1605                 if(nchrg.lt.0) sigma0=akNsgm(pkaon)
1606
1607 cbz3/8/99 neutralk
1608                 SIGMA0 = 1.5 * AKPSGM(PKAON) + AKNSGM(PKAON)
1609 cbz3/8/99 neutralk end
1610
1611              endif
1612              sig=(srt**2-(aka+bmass)**2)*(srt**2-(aka-bmass)**2)/
1613      &         (srt**2-(em1+em2)**2)/(srt**2-(em1-em2)**2)*sigma0
1614 cbz3/8/99 neutralk
1615 c     if(nchrg.eq.-2.or.nchrg.eq.1) sig=2.*sig K-D++, K-D-
1616 c     K0barD++, K-D-
1617              if(nchrg.eq.-2.or.nchrg.eq.2) sig=2.*sig
1618
1619 cbz3/8/99 neutralk end
1620
1621 c             the factor 2 comes from spin of delta, which is 3/2
1622 c             detailed balance. copy from Page 423 of N.P. A614 1997
1623
1624 cbz3/8/99 neutralk
1625              IF (LB1 .EQ. 14 .OR. LB2 .EQ. 14) THEN
1626                 SIG = 4.0 / 3.0 * SIG
1627              ELSE IF (NCHRG .EQ. -2 .OR. NCHRG .EQ. 2) THEN
1628                 SIG = 8.0 / 9.0 * SIG
1629              ELSE
1630                 SIG = 4.0 / 9.0 * SIG
1631              END IF
1632 cbz3/8/99 neutralk end
1633              brsig = sig
1634              if(sig.lt.1.e-7) sig = 1.e-7
1635            endif
1636 csp05/07/01
1637 * comment icase=4 statement below if only inelastic
1638 c     PI+L/Si => Kbar + N  OR ELASTIC SCATTERING
1639            icase=4
1640            brsig = sig
1641 c     elastic xsecn of 10mb
1642            sigela = 10.
1643            sig = sig + sigela
1644            brel = sigela/sig
1645 cc          brsig = sig
1646 csp05/07/01 end
1647         endif
1648 c
1649 c        if(em2.lt.0.2.and.em1.lt.0.2) then
1650 c     PI + PI
1651 c             icase=5
1652 c     assumed PI PI total x section.
1653 c              sig=50.
1654 c     Mk + Mkbar
1655 c              s0=aka+aka
1656 c              brsig = 0.
1657 c              if(srt.gt.s0) brsig = 2.7*(1.-s0**2/srt**2)**0.76
1658 c              x section for PIPI->KKbar   PRC43 (1991) 1881
1659 c        endif
1660         if(icase.eq.-1) then
1661            ictrl = -1
1662            return
1663         endif
1664         px1cm=pcx
1665         py1cm=pcy
1666         pz1cm=pcz
1667         ds=sqrt(sig/31.4)
1668         dsr=ds+0.1
1669         ec=(em1+em2+0.02)**2
1670 c        ec=3.59709
1671 c        if((e(i1).ge.1.).and.(e(i2).ge.1.)) ec = 4.75
1672
1673         call distce(i1,i2,dsr,ds,dt,ec,srt,ic,px1cm,py1cm,pz1cm)
1674         if(ic.eq.-1) then
1675 c     no anti-kaon production
1676            ictrl = -1
1677 c           in=in+1
1678 c           write(60,*)'--------------distance-----',in
1679            return
1680         endif
1681
1682 clin-10/24/02 set to 0: ik,ik0-3,il,im,im3-4,in,inpion,ipipi,
1683 c     sgsum,sgsum1,sgsum3:
1684         ik=0
1685         ik0=0
1686         ik1=0
1687         ik2=0
1688         ik3=0
1689         il=0
1690         im=0
1691         im3=0
1692         im4=0
1693         in=0
1694         inpion=0
1695         ipipi=0
1696         sgsum=0.
1697         sgsum1=0.
1698         sgsum3=0.
1699         if(icase.eq.1) then
1700            ik=ik+1
1701            if(srt.gt.2.8639) then
1702                 ik0=ik0+1
1703                 if(em1.lt.1.0.and.em2.lt.1.0) then
1704                         ik1=ik1+1
1705                         sgsum1=sgsum1+brsig
1706 c                        ratio_1=sgsum1/ik1/40.
1707                 endif
1708                 if(em1.gt.1.0.and.em2.gt.1.0) then
1709                         ik3=ik3+1
1710                         sgsum3=sgsum3+brsig
1711 c                        ratio_3=sgsum3/ik3/40.
1712                 endif
1713                 if(em1.gt.1.0.and.em2.lt.1.0) ik2=ik2+1
1714                 if(em1.lt.1.0.and.em2.gt.1.0) ik2=ik2+1
1715                 sgsum=sgsum+brsig
1716 c                ratio=sgsum/ik0/40.
1717            endif
1718         endif
1719         if(icase.eq.2) inpion=inpion+1
1720         if(icase.eq.5) ipipi=ipipi+1
1721 c        write(62,*)'ik1,ik2,ik3',ik1,ik2,ik3,ratio_1,ratio_3,ratio
1722 c        write(62,*)'inpion,ipipi',inpion,ipipi
1723         if(RANART(NSEED).gt.(brsig/sig)) then
1724 c     no anti-kaon production
1725            ictrl = -1
1726            return
1727         endif
1728         il=il+1
1729 c        kaons could be created now.
1730         if(icase.eq.1) then
1731           in=in+1
1732 c          write(60,*)'------in,s2kaon,sig=',in,brsig,sig,lb1,lb2
1733           call nnkaon(irun,iseed,
1734      &          ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1735         endif
1736         if(icase.eq.2) then
1737           im=im+1
1738 c          call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1739 c     &              pcx,pcy,pcz,nchrg,ratiok)
1740           call npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
1741      &              pcx,pcy,pcz,nchrg,ratiok,iblock)
1742         endif
1743 c
1744         if(icase.eq.3) then
1745           im3=im3+1
1746 c          write(63,*)'im3,lb1,lb2,pkaon',im3,lb1,lb2,pkaon
1747 c          write(63,*)'sig,el,sigma',sig,brel,brsgm
1748 c          write(63,*)'srt,pcx,pcy,pcz,em1,em2',srt,pcx,pcy,pcz,em1,em2
1749           call kaonN(brel,brsgm,irun,iseed,dt,nt,ictrl,
1750      &                i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
1751 c         this subroutine format is diff. since three final states are possible
1752         endif
1753 c
1754
1755         if(icase.eq.4) then
1756           im4=im4+1
1757 c          write(64,*)'im4,sigma0,branch,sig=',im4,sigma0,brsig,sig
1758 c          write(64,*)'lb1,lb2,em1,em2,pkaon=',lb1,lb2,em1,em2,pkaon
1759
1760 csp06/07/01
1761       if(RANART(NSEED).lt.brel) then
1762          ielstc = 1
1763       else
1764          ielstc = 0
1765       endif
1766 c          call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1767 c     &                   pcx,pcy,pcz,nchrg)
1768           call Pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,srt,
1769      &                   pcx,pcy,pcz,nchrg,iblock)
1770
1771 csp06/07/01 end
1772         endif
1773 c        if(icase.eq.5) then
1774 c          im5=im5+1
1775 c          write(65,*)'---im5,s2kaon,sig=',im5,brsig,sig
1776 c          call pipikaon(irun,iseed,dt,nt,ictrl,i1,i2,srt,pcx,pcy,pcz)
1777 c        endif
1778 cbz3/2/99
1779 c        write(101,*)lb1,lb2,lb(i1),lb(i2)
1780 c        write(101,*)em1,em2,e(i1),e(i2),srt
1781 cbz3/2/99end
1782
1783         return
1784         end
1785
1786 ******************************************
1787 * for pp-->pp + kaon + anti-kaon
1788 c      real*4 function X2kaon(srt)
1789       real function X2kaon(srt)
1790       SAVE
1791 *  This function contains the experimental total pp->pp+K(+)K(-) Xsections    *
1792 *  srt    = DSQRT(s) in GeV                                                   *
1793 *  xsec   = production cross section in mb                                    *
1794 *                                                                             *
1795 ******************************************
1796 c     minimum c.m.s. energy to create 2 kaon. = 2*(mp+mk)
1797         smin = 2.8639
1798         x2kaon=0.0000001
1799         if(srt.lt.smin)return
1800         sigma1 = 2.8
1801         sigma2 = 7.7
1802         sigma3 = 3.9
1803         x = srt**2/smin**2 + 0.0000001
1804         f1 = (1.+1./sqrt(x))*alog(x) - 4.*(1.-1./sqrt(x))
1805         f2 = 1. - (1./sqrt(x))*(1.+alog(sqrt(x)))
1806         f3 = ((x-1.)/x**2)**3.5
1807         x2kaon = (1.-1./x)**3*(sigma1*f1 + sigma2*f2) + sigma3*f3
1808         return
1809         END
1810
1811         real function piNsg0(srt)
1812       SAVE
1813 * cross section in mb for PI- + P -> P + K0 + K-
1814 c     Mn + 2* Mk
1815         srt0 = 0.938 + 2.*0.498
1816         if(srt.lt.srt0) then
1817            piNsg0 = 0.0
1818            return
1819         endif
1820         ratio = srt0**2/srt**2
1821         piNsg0=1.121*(1.-ratio)**1.86*ratio**2
1822         return
1823         end
1824
1825         real function akNel(pkaon)
1826       SAVE
1827 *cross section in mb for K- + N reactions.
1828 c        the following data come from PRC 41 (1701)
1829 c        sigma1: K(-) + neutron elastic
1830         if(pkaon.lt.0.5.or. pkaon.ge.4.0) sigma1=0.
1831         if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma1=20.*pkaon**2.74
1832         if(pkaon.ge.1.0.and.pkaon.lt.4.0) sigma1=20.*pkaon**(-1.8)
1833         akNel=sigma1
1834         return
1835         end
1836
1837         real function akPel(pkaon)
1838       SAVE
1839 *cross section in mb for K- + N reactions.
1840 c        the following data come from PRC 41 (1701)
1841 c        sigma2: K(-) + proton elastic
1842         if(pkaon.lt.0.25.or. pkaon.ge.4.0) sigma2=0.
1843         if(pkaon.ge.0.25.and.pkaon.lt.4.0) sigma2=13.*pkaon**(-0.9)
1844         akPel=sigma2
1845         return
1846         end
1847
1848         real function akNsgm(pkaon)
1849       SAVE
1850 *cross section in mb for K- + N reactions.
1851 c        sigma2: x section for K- + n -> sigma0 + PI-
1852         if(pkaon.lt.0.5.or. pkaon.ge.6.0) sigma2=0.
1853         if(pkaon.ge.0.5.and.pkaon.lt.1.0) sigma2=1.2*pkaon**(-1.3)
1854         if(pkaon.ge.1.0.and.pkaon.lt.6.0) sigma2=1.2*pkaon**(-2.3)
1855         akNsgm=sigma2
1856         return
1857         end
1858
1859         real function akPsgm(pkaon)
1860       SAVE
1861 *cross section in mb for K- + N reactions.
1862 c        sigma1: x section for K- + p -> sigma0 + PI0
1863         if(pkaon.lt.0.2.or. pkaon.ge.1.5) sigma1=0.
1864         if(pkaon.ge.0.2.and.pkaon.lt.1.5) sigma1=0.6*pkaon**(-1.8)
1865         akPsgm=sigma1
1866         return
1867         end
1868
1869         real function akPlam(pkaon)
1870       SAVE
1871 *cross section in mb for K- + N reactions.
1872 c        sigma: x section for K- + p -> lambda + PI0
1873         p=pkaon
1874         if(pkaon.lt.0.2.or. pkaon.ge.10.0) sigma=0.
1875         if(pkaon.ge.0.2.and.pkaon.lt.0.9) sigma=50.*p**2-67.*p+24.
1876         if(pkaon.ge.0.9.and.pkaon.lt.10.0) sigma=3.0*pkaon**(-2.6)
1877         akPlam=sigma
1878         return
1879         end
1880
1881         real function akNlam(pkaon)
1882       SAVE
1883 *cross section in mb for K- + N reactions.
1884         akNlam=akPlam(pkaon)
1885         return
1886         end
1887
1888 * GQ Li parametrization (without resonance)
1889         real function akNPsg(pkaon)
1890       SAVE
1891 *cross section in mb for K- + N reactions.
1892 c       sigma1: x section for K- + p/n -> sigma0 + PI0
1893          if(pkaon.le.0.345)then
1894            sigma1=0.624*pkaon**(-1.83)
1895          else
1896            sigma1=0.7*pkaon**(-2.09)
1897          endif
1898         akNPsg=sigma1
1899         return
1900         end
1901
1902 c-----------------------------------------------------------------------
1903
1904 c.....extracted from G. Song's ART expasion including K- interactions
1905 c.....file `NEWNNK.FOR'
1906
1907         subroutine nnkaon(irun,iseed,ictrl,i1,i2,iblock,
1908      &                                   srt,pcx,pcy,pcz,nchrg)
1909 c        <pt>=0.27+0.037*log(srt) was changed to 0.632 + ... on Aug. 14, 1997
1910 c     CANCELED also alpha=1 changed to alpha=3 to decrease the leadng effect.
1911       PARAMETER      (MAXSTR=150001,MAXR=1)
1912       PARAMETER      (AKA=0.498)
1913       COMMON   /AA/  R(3,MAXSTR)
1914 cc      SAVE /AA/
1915       COMMON   /BB/  P(3,MAXSTR)
1916 cc      SAVE /BB/
1917       COMMON   /CC/  E(MAXSTR)
1918 cc      SAVE /CC/
1919       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
1920 cc      SAVE /EE/
1921       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
1922 cc      SAVE /BG/
1923       COMMON   /NN/NNN
1924 cc      SAVE /NN/
1925       COMMON   /RUN/NUM
1926 cc      SAVE /RUN/
1927       COMMON   /PA/RPION(3,MAXSTR,MAXR)
1928 cc      SAVE /PA/
1929       COMMON   /PB/PPION(3,MAXSTR,MAXR)
1930 cc      SAVE /PB/
1931       COMMON   /PC/EPION(MAXSTR,MAXR)
1932 cc      SAVE /PC/
1933       COMMON   /PD/LPION(MAXSTR,MAXR)
1934 cc      SAVE /PD/
1935       dimension px(4),py(4),pz(4)
1936       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1937      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1938      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1939       SAVE
1940 c      dm1=e(i1)
1941 c      dm2=e(i2)
1942       dm3=0.938
1943       dm4=0.938
1944 c     10/24/02 initialize n to 0:
1945       n=0
1946
1947 cbz3/11/99 neutralk
1948 c        if(nchrg.eq.-2.or.nchrg.ge.3) dm3=1.232
1949 c        if(nchrg.eq.4) dm4=1.232
1950         if(nchrg.le.-1.or.nchrg.ge.3) dm3=1.232
1951         if(nchrg.eq.-2.or.nchrg.eq.4) dm4=1.232
1952 cbz3/11/99 neutralk end
1953           iblock = 0
1954         call fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
1955         if(iflag.lt.0) then
1956 c           write(60,*)'------------final state fail-------',n
1957 c     no anti-kaon production
1958            ictrl = -1
1959            n=n+1
1960            return
1961         endif
1962         iblock = 12
1963 * Rotate the momenta of particles in the cms of I1 & I2
1964 * px(1), py(1), pz(1): momentum of I1
1965 * px(2), py(2), pz(2): momentum of I2
1966 * px(3), py(3), pz(3): momentum of anti-kaon
1967 * px(4), py(4), pz(4): momentum of kaon
1968
1969
1970 c     10/28/02 get rid of argument usage mismatch in rotate():
1971         pxrota=px(1)
1972         pyrota=py(1)
1973         pzrota=pz(1)
1974 c        call rotate(pcx,pcy,pcz,px(1),py(1),pz(1))
1975         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1976         px(1)=pxrota
1977         py(1)=pyrota
1978         pz(1)=pzrota
1979 c
1980         pxrota=px(2)
1981         pyrota=py(2)
1982         pzrota=pz(2)
1983 c        call rotate(pcx,pcy,pcz,px(2),py(2),pz(2))
1984         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1985         px(2)=pxrota
1986         py(2)=pyrota
1987         pz(2)=pzrota
1988 c
1989         pxrota=px(3)
1990         pyrota=py(3)
1991         pzrota=pz(3)
1992 c        call rotate(pcx,pcy,pcz,px(3),py(3),pz(3))
1993         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
1994         px(3)=pxrota
1995         py(3)=pyrota
1996         pz(3)=pzrota
1997 c
1998         pxrota=px(4)
1999         pyrota=py(4)
2000         pzrota=pz(4)
2001 c        call rotate(pcx,pcy,pcz,px(4),py(4),pz(4))
2002         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2003         px(4)=pxrota
2004         py(4)=pyrota
2005         pz(4)=pzrota
2006
2007         nnn=nnn+2
2008 c     K+
2009         lpion(nnn,irun)=23
2010         if(nchrg.eq.-1.or.nchrg.eq.-2) then
2011 c        To keep charge conservation. D-n->nnK0K-, D-D- -> nD-K0K-
2012
2013 cbz3/7/99 neutralk
2014 c           lpion(nnn,irun)=24 ! K0
2015 cbz3/7/99 neutralk end
2016
2017         endif
2018 c     aka: rest mass of K
2019         epion(nnn,irun)=aka
2020 c     K-
2021         lpion(nnn-1,irun)=21
2022 c     aka: rest mass of K
2023         epion(nnn-1,irun)=aka
2024 * Find the momenta of particles in the final state in the nucleus_nucleus
2025 * cms frame.   Lorentz transformation into lab frame.
2026         e1cm   = sqrt(dm3**2 + px(1)**2 + py(1)**2 + pz(1)**2)
2027         p1beta = px(1)*betax + py(1)*betay + pz(1)*betaz
2028         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2029         pt1i1 = betax*transf + px(1)
2030         pt2i1 = betay*transf + py(1)
2031         pt3i1 = betaz*transf + pz(1)
2032         eti1  = dm3
2033 c        lb1   = lb(i1)
2034         lb1   = 2
2035         if(nchrg.ge.-2.and.nchrg.le.1) lb1=2
2036
2037 cbz3/7/99 neutralk
2038         if (nchrg .eq. -2 .or. nchrg .eq. -1) then
2039            lb1 = 6
2040         end if
2041 cbz3/7/99 neutralk end
2042
2043 cbz3/11/99 neutralk
2044 c        if(nchrg.eq.2.or.nchrg.eq.3) lb1=1
2045 c        if(nchrg.eq.4) lb1=9
2046         if(nchrg.eq.1.or.nchrg.eq.2) lb1=1
2047         if(nchrg.eq.3.or.nchrg.eq.4) lb1=9
2048 cbz3/11/99 neutralk end
2049
2050 * For second nulceon, same
2051         e2cm   = sqrt(dm4**2 + px(2)**2 + py(2)**2 + pz(2)**2)
2052         p2beta = px(2)*betax + py(2)*betay + pz(2)*betaz
2053         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2054         pt1i2 = betax*transf + px(2)
2055         pt2i2 = betay*transf + py(2)
2056         pt3i2 = betaz*transf + pz(2)
2057         eti2  = dm4
2058 c        lb2   = lb(i2)
2059         lb2   = 2
2060
2061 cbz3/11/99 neutralk
2062 c        if(nchrg.eq.-1.or.nchrg.eq.0) lb2=2
2063 c        if(nchrg.eq. 2.or.nchrg.eq.1) lb2=1
2064 c        if(nchrg.eq. 4.or.nchrg.eq.3) lb2=9
2065 c        if(nchrg.eq.-2) lb2=6
2066         if(nchrg.ge.-1.or.nchrg.le.1) lb2=2
2067         if(nchrg.eq. 2.or.nchrg.eq.3) lb2=1
2068         if(nchrg.eq. 4) lb2=9
2069         if(nchrg.eq.-2) lb2=6
2070 cbz3/11/99 neutralk end
2071
2072 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2073                 p(1,i1)=pt1i1
2074                 p(2,i1)=pt2i1
2075                 p(3,i1)=pt3i1
2076                 e(i1)=eti1
2077                 lb(i1)=lb1
2078                 p(1,i2)=pt1i2
2079                 p(2,i2)=pt2i2
2080                 p(3,i2)=pt3i2
2081                 e(i2)=eti2
2082                 lb(i2)=lb2
2083
2084 c                px1 = p(1,i1)
2085 c                py1 = p(2,i1)
2086 c                pz1 = p(3,i1)
2087 c                em1 = e(i1)
2088 c                id(i1) = 2
2089 c                id(i2) = 2
2090 c                id1 = id(i1)
2091 c                iblock = 101  ! K(+)K(-) production
2092 * Get anti-kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2093         epcmk = sqrt(epion(nnn-1,irun)**2 + px(3)**2+py(3)**2+pz(3)**2)
2094         betak = px(3)*betax + py(3)*betay + pz(3)*betaz
2095         transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2096         ppion(1,nnn-1,irun)=betax*transf + px(3)
2097         ppion(2,nnn-1,irun)=betay*transf + py(3)
2098         ppion(3,nnn-1,irun)=betaz*transf + pz(3)
2099         rpion(1,nnn-1,irun)=r(1,i1)
2100         rpion(2,nnn-1,irun)=r(2,i1)
2101         rpion(3,nnn-1,irun)=r(3,i1)
2102 clin-5/2008:
2103         dppion(nnn-1,irun)=dpertp(i1)*dpertp(i2)
2104 * Same thing for kaon **************************************
2105         epcmak = sqrt(epion(nnn,irun)**2 + px(4)**2 +py(4)**2+pz(4)**2)
2106         betaak = px(4)*betax + py(4)*betay + pz(4)*betaz
2107         transf= gamma*(gamma*betaak/(gamma+1.) + epcmak)
2108         ppion(1,nnn,irun)=betax*transf + px(4)
2109         ppion(2,nnn,irun)=betay*transf + py(4)
2110         ppion(3,nnn,irun)=betaz*transf + pz(4)
2111         rpion(1,nnn,irun)=r(1,i2)
2112         rpion(2,nnn,irun)=r(2,i2)
2113         rpion(3,nnn,irun)=r(3,i2)
2114 clin-5/2008:
2115         dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2116         return
2117         end
2118
2119         subroutine lorntz(ilo,b,pi,pj)
2120 c       It uses to perform Lorentz (or inverse Lorentz) transformation
2121         dimension pi(4),pj(4),b(3)
2122       SAVE
2123 c       dimension db(3)
2124         bb=b(1)*b(1)+b(2)*b(2)+b(3)*b(3)
2125         deno3=sqrt(1.-bb)
2126         if(deno3.eq.0.)deno3=1.e-10
2127         gam=1./deno3
2128         ga=gam*gam/(gam+1.)
2129         if(ilo.eq.1) goto 100
2130 c       Lorentz transformation
2131         pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2132         pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2133 c       drb=drd(1)*b(1)+drd(2)*b(2)+drd(3)*b(3)
2134 c       drdb=db(1)*b(1)+db(2)*b(2)+db(3)*b(3)
2135         do 1001 i=1,3
2136            pi(i)=pi(i)+b(i)*(ga*pib-gam*pi(4))
2137            pj(i)=pj(i)+b(i)*(ga*pjb-gam*pj(4))
2138 c       drd(i)=drd(i)+b(i)*ga*drb
2139 c       db(i)=db(i)+b(i)*ga*drdb
2140  1001   continue
2141         pi(4)=gam*(pi(4)-pib)
2142         pj(4)=gam*(pj(4)-pjb)
2143         return
2144 100     continue
2145 c       inverse Lorentz transformation
2146         pib=pi(1)*b(1)+pi(2)*b(2)+pi(3)*b(3)
2147         pjb=pj(1)*b(1)+pj(2)*b(2)+pj(3)*b(3)
2148         do 1002 i=1,3
2149            pi(i)=pi(i)+b(i)*(ga*pib+gam*pi(4))
2150            pj(i)=pj(i)+b(i)*(ga*pjb+gam*pj(4))
2151  1002   continue
2152         pi(4)=gam*(pi(4)+pib)
2153         pj(4)=gam*(pj(4)+pjb)
2154         return
2155         end
2156
2157         subroutine fstate(iseed,srt,dm3,dm4,px,py,pz,iflag)
2158 *        function: decide final momentum for N,N,K(+),and K(-)
2159         dimension px(4), py(4), pz(4), pe(4)
2160         COMMON/RNDF77/NSEED
2161 cc      SAVE /RNDF77/
2162         SAVE
2163
2164         iflag=-1
2165 c        iflag=-1: fail to find momenta
2166 c             = 1: success
2167         pio=3.1415926
2168         aka=0.498
2169 c        v=0.43
2170 c        w=-0.84
2171 c        b=3.78
2172 c        c=0.47
2173 c        d=3.60
2174 c        fmax=1.056
2175 c        gmax=1.+c
2176
2177         icount=0
2178         ekmax=(srt-dm3-dm4)/2.
2179         if(ekmax.le.aka) return
2180         pkmax=sqrt(ekmax**2-aka**2)
2181
2182         if(dm3.le.0.0.or.dm4.le.0.0) then
2183            write(1,*)'error: minus mass!!!'
2184            return
2185         endif
2186
2187 c        after we have the momenta for both nucleus, we sample the
2188 c        transverse momentum for K-.
2189 c        dsigma/dpt**2 = exp(-4.145*pt**2) obtained by fitting data on
2190 c        page 72, fig 23i.
2191 50        continue
2192         icount=icount+1
2193         if(icount.gt.10) return
2194         ptkmi2=-1./4.145*alog(RANART(NSEED))
2195         ptkm=sqrt(ptkmi2)
2196 3        v1=RANART(NSEED)
2197         v2=RANART(NSEED)
2198         rsq=v1**2+v2**2
2199         if(rsq.ge.1.0.or.rsq.le.0.) goto 3
2200         fac=sqrt(-2.*alog(rsq)/rsq)
2201         guass=v1*fac
2202         if(guass.ge.5.) goto 3
2203         xstar=guass/5.
2204         pzkm=pkmax*xstar
2205         ekm=sqrt(aka**2+pzkm**2+ptkm**2)
2206         if(RANART(NSEED).gt.aka/ekm) goto 50
2207         bbb=RANART(NSEED)
2208         px(3)=ptkm*cos(2.*pio*bbb)
2209         py(3)=ptkm*sin(2.*pio*bbb)
2210         if(RANART(NSEED).gt.0.5) pzkm=-1.*pzkm
2211         pz(3)=pzkm
2212         pe(3)=ekm
2213 150        ptkpl2=-1./3.68*alog(RANART(NSEED))
2214         ptkp=sqrt(ptkpl2)
2215 13        v1=RANART(NSEED)
2216         v2=RANART(NSEED)
2217         rsq=v1**2+v2**2
2218         if(rsq.ge.1.0.or.rsq.le.0.) goto 13
2219         fac=sqrt(-2.*alog(rsq)/rsq)
2220         guass=v1*fac
2221         if(guass.ge.3.25) goto 13
2222         xstar=guass/3.25
2223         pzkp=pkmax*xstar
2224         ekp=sqrt(aka**2+pzkp**2+ptkp**2)
2225         if(RANART(NSEED).gt.aka/ekp) goto 150
2226         bbb=RANART(NSEED)
2227         px(4)=ptkp*cos(2.*pio*bbb)
2228         py(4)=ptkp*sin(2.*pio*bbb)
2229         if(RANART(NSEED).gt.0.5) pzkp=-1.*pzkp
2230         pz(4)=pzkp
2231         pe(4)=ekp
2232
2233         resten=srt-pe(3)-pe(4)
2234         restpz=-pz(3)-pz(4)
2235 c     resample
2236         if(resten.le.abs(restpz)) goto 50
2237         restms=sqrt(resten**2-restpz**2)
2238 c     resample
2239         if(restms.lt.(dm3+dm4)) goto 50
2240         ptp2=-1./2.76*alog(RANART(NSEED))
2241         ptp=sqrt(ptp2)
2242         bbb=RANART(NSEED)
2243         px(2)=ptp*cos(2.*pio*bbb)
2244         py(2)=ptp*sin(2.*pio*bbb)
2245         px(1)=-1.*(px(4)+px(3)+px(2))
2246         py(1)=-1.*(py(4)+py(3)+py(2))
2247 c     transverse mass for K-
2248         rmt3=sqrt(dm3**2+px(1)**2+py(1)**2)
2249 c     transverse mass for K+
2250         rmt4=sqrt(dm4**2+px(2)**2+py(2)**2)
2251         if(restms.lt.(rmt3+rmt4)) goto 50
2252 c        else: sampling success!
2253         pzcms=sqrt((restms**2-(rmt3+rmt4)**2)*
2254      &             (restms**2-(rmt3-rmt4)**2))/2./restms
2255         if(RANART(NSEED).gt.0.5) then
2256            pz(1)=pzcms
2257            pz(2)=-pzcms
2258         else
2259            pz(1)=-pzcms
2260            pz(2)=pzcms
2261         endif
2262         beta=restpz/resten
2263         gama=1./sqrt(1.-beta**2)
2264         pz(1)=pz(1)*gama + beta*gama*sqrt(rmt3**2+pz(1)**2)
2265         pz(2)=pz(2)*gama + beta*gama*sqrt(rmt4**2+pz(2)**2)
2266         pe(1)=sqrt(rmt3**2+pz(1)**2)
2267         pe(2)=sqrt(rmt4**2+pz(2)**2)
2268
2269         iflag=1
2270         return
2271         end
2272
2273 c-----------------------------------------------------------------------
2274
2275 c.....extracted from G. Song's ART expasion including K- interactions
2276 c.....file `NPIK.FOR'
2277
2278 ****************************************
2279 c        subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2280 c     &                  pcx,pcy,pcz,nchrg,ratiok)
2281         subroutine npik(irun,iseed,dt,nt,ictrl,i1,i2,srt,
2282      &                  pcx,pcy,pcz,nchrg,ratiok,iblock)
2283 *
2284 * Process: PI + N -> K(-) + ANYTHING
2285 * 1.  PI- + P -> P + K0 + K-
2286 * 2.  PI+ + N -> P + K+ + K-
2287 * 3.  PI0 + P -> P + K+ + K-
2288 * 4.  PI0 + N -> P + K0 + K-
2289 * 5.  PI0 + N -> N + K+ + K-
2290 * 6.  PI- + P -> N + K+ + K-
2291 * 7.  PI- + N -> N + K0 + K-
2292 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2293 ****************************************
2294       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2295       PARAMETER      (AKA=0.498)
2296       COMMON   /AA/  R(3,MAXSTR)
2297 cc      SAVE /AA/
2298       COMMON   /BB/  P(3,MAXSTR)
2299 cc      SAVE /BB/
2300       COMMON   /CC/  E(MAXSTR)
2301 cc      SAVE /CC/
2302       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2303 cc      SAVE /EE/
2304       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2305 cc      SAVE /BG/
2306       COMMON   /NN/NNN
2307 cc      SAVE /NN/
2308       COMMON   /RUN/NUM
2309 cc      SAVE /RUN/
2310       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2311 cc      SAVE /PA/
2312       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2313 cc      SAVE /PB/
2314       COMMON   /PC/EPION(MAXSTR,MAXR)
2315 cc      SAVE /PC/
2316       COMMON   /PD/LPION(MAXSTR,MAXR)
2317 cc      SAVE /PD/
2318       dimension bb(3),p1(4),p2(4),p3(4)
2319 cms      ,px(4),py(4),pz(4)
2320       COMMON/RNDF77/NSEED
2321 cc      SAVE /RNDF77/
2322       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
2323      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
2324      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
2325       SAVE
2326         px1cm=pcx
2327         py1cm=pcy
2328         pz1cm=pcz
2329         ictrl = 1
2330         lb1=lb(i1)
2331         lb2=lb(i2)
2332         k1=i1
2333         k2=i2
2334 c        k1 must be bayron. k2 be meson. If not, exchange.
2335         if(lb2.eq.1.or.lb2.eq.2.or.(lb2.ge.6.and.lb2.le.13)) then
2336             k1=i2
2337             k2=i1
2338         endif
2339 cbz3/8/99 neutralk
2340 cbz10/12/99
2341 c        LB(I1) = 1 + 2 * RANART(NSEED)
2342 c        LB(I2) = 23
2343         LB(k1) = 1 + int(2*RANART(NSEED))
2344         LB(k2) = 23
2345 c       pkmax=sqrt((srt**2-(aka+0.938+aka)**2)*(srt**2-(aka+0.938-aka)**2))
2346 c     &           /2./srt
2347         pkmax=sqrt((srt**2-(aka+0.938+aka)**2)
2348      &           *(srt**2-(aka+0.938-aka)**2))/2./srt
2349         pk = RANART(NSEED)*pkmax
2350 c-----------------------------------------------------
2351         css=1.-2.*RANART(NSEED)
2352         sss=sqrt(1.-css**2)
2353         fai=2*3.1415926*RANART(NSEED)
2354         p3(1)=pk*sss*cos(fai)
2355         p3(2)=pk*sss*sin(fai)
2356         p3(3)=pk*css
2357         eip = srt - sqrt(aka**2 + pk**2)
2358         rmnp=sqrt(eip**2-pk**2)
2359         do 1001 i= 1, 3
2360            bb(i) = -1.*p3(i)/eip
2361  1001   continue
2362 c        bb: velocity of the other two particles as a whole.
2363         pznp=sqrt((rmnp**2-(aka+0.938)**2)
2364      c  *(rmnp**2-(0.938-aka)**2))/2./rmnp
2365 c-----------------------------------------------------
2366         css=1.-2.*RANART(NSEED)
2367         sss=sqrt(1.-css**2)
2368         fai=2*3.1415926*RANART(NSEED)
2369         p1(1)=pznp*sss*cos(fai)
2370         p1(2)=pznp*sss*sin(fai)
2371         p1(3)=pznp*css
2372         p1(4)=sqrt(0.938**2+pznp**2)
2373         p2(4)=sqrt(aka**2+pznp**2)
2374         do 1002 i=1,3
2375            p2(i)=-1.*p1(i)
2376  1002   continue
2377 c        p1,p2: the momenta of the two particles in their cms
2378 c        p1: momentum of N or P
2379 c        p2: momentum of anti_kaon
2380 c        p3: momentum of K0 or K+
2381         ilo=1
2382 c        write(61,*)'--------p1,p2',p1,p2
2383 c        write(61,*)'--------bb',bb
2384         call lorntz(ilo,bb,p1,p2)
2385 c******* Checking *************
2386 c        pxsum = p1(1)+p2(1)+p3(1)
2387 c        pysum = p1(2)+p2(2)+p3(2)
2388 c        pzsum = p1(3)+p2(3)+p3(3)
2389 c        pesum = p1(4)+p2(4)+sqrt(p3(1)**2+p3(2)**2+p3(3)**2+aka**2)-srt
2390 c        write(61,*)'---p1,pxsum',p1,pxsum
2391 c        write(61,*)'---p2,pysum',p2,pysum
2392 c        write(61,*)'---p3,pzsum',p3,pzsum
2393 c        write(61,*)'---pesum',pesum
2394 c***********************************
2395
2396 * Rotate the momenta of particles in the cms of I1 & I2
2397 * px(1), py(1), pz(1): momentum of I1
2398 * px(2), py(2), pz(2): momentum of I2
2399 * px(3), py(3), pz(3): momentum of anti-kaon
2400
2401 c     10/28/02 get rid of argument usage mismatch in rotate():
2402         pxrota=p1(1)
2403         pyrota=p1(2)
2404         pzrota=p1(3)
2405 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2406         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2407         p1(1)=pxrota
2408         p1(2)=pyrota
2409         p1(3)=pzrota
2410 c
2411         pxrota=p2(1)
2412         pyrota=p2(2)
2413         pzrota=p2(3)
2414 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2415         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2416         p2(1)=pxrota
2417         p2(2)=pyrota
2418         p2(3)=pzrota
2419 c
2420         pxrota=p3(1)
2421         pyrota=p3(2)
2422         pzrota=p3(3)
2423 c        call rotate(pcx,pcy,pcz,p3(1),p3(2),p3(3))
2424         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2425         p3(1)=pxrota
2426         p3(2)=pyrota
2427         p3(3)=pzrota
2428
2429         nnn=nnn+1
2430 c     K(-)
2431         lpion(nnn,irun)=21
2432 c     aka: rest mass of K
2433         epion(nnn,irun)=aka
2434 * Find the momenta of particles in the final state in the nucleus_nucleus
2435 * cms frame.   Lorentz transformation into lab frame.
2436         e1cm   = sqrt(0.938**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2437         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2438         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2439         pt1i1 = betax*transf + p1(1)
2440         pt2i1 = betay*transf + p1(2)
2441         pt3i1 = betaz*transf + p1(3)
2442         eti1  = 0.938
2443         lb1   = lb(k1)
2444
2445 * For second nulceon, same
2446         e2cm   = sqrt(aka**2 + p3(1)**2 + p3(2)**2 + p3(3)**2)
2447         p2beta = p3(1)*betax + p3(2)*betay + p3(3)*betaz
2448         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2449         pt1i2 = betax*transf + p3(1)
2450         pt2i2 = betay*transf + p3(2)
2451         pt3i2 = betaz*transf + p3(3)
2452         eti2  = aka
2453         lb2   = lb(k2)
2454
2455 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2456 *       k1 stand for nucleon, k2 stand for kaon. lpion stand for Kbar.
2457                 p(1,k1)=pt1i1
2458                 p(2,k1)=pt2i1
2459                 p(3,k1)=pt3i1
2460                 e(k1)=eti1
2461                 lb(k1)=lb1
2462                 p(1,k2)=pt1i2
2463                 p(2,k2)=pt2i2
2464                 p(3,k2)=pt3i2
2465                 e(k2)=eti2
2466                 lb(k2)=lb2
2467
2468 c                px1 = p(1,i1)
2469 c                py1 = p(2,i1)
2470 c                pz1 = p(3,i1)
2471 c                em1 = e(i1)
2472 c                id(i1) = 2
2473 c                id(i2) = 2
2474 c                id1 = id(i1)
2475 c     K(+)K(-) production
2476                 iblock = 101
2477 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2478 c  p2:  momentum of anti-kaon.
2479 c        epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2480         epcmk = sqrt(epion(nnn,irun)**2 + p2(1)**2+p2(2)**2+p2(3)**2)
2481         betak = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2482         transf= gamma*(gamma*betak/(gamma+1.) + epcmk)
2483         ppion(1,nnn,irun)=betax*transf + p2(1)
2484         ppion(2,nnn,irun)=betay*transf + p2(2)
2485         ppion(3,nnn,irun)=betaz*transf + p2(3)
2486 clin-5/2008:
2487         dppion(nnn,irun)=dpertp(i1)*dpertp(i2)
2488 cbz3/2/99
2489 c        write(400,*)'2 ', ppion(1,nnn,irun), ppion(2,nnn,irun),
2490 c     &                    ppion(3,nnn,irun), dt*nt, srt
2491 cbz3/2/99end
2492 c        write(420,*)ppion(1,nnn,irun), ppion(2,nnn,irun),
2493 c     &                    ppion(3,nnn,irun), dt*nt, srt
2494         k=i2
2495         if(lb(i1).eq.1.or.lb(i1).eq.2) k=i1
2496         rpion(1,nnn,irun)=r(1,k)
2497         rpion(2,nnn,irun)=r(2,k)
2498         rpion(3,nnn,irun)=r(3,k)
2499         return
2500         end
2501
2502 c-----------------------------------------------------------------------
2503
2504 c.....extracted from G. Song's ART expasion including K- interactions
2505 c.....file `PIHYPN.FOR'
2506
2507 ******************************************
2508         subroutine pihypn(ielstc,irun,iseed,dt,nt,ictrl,i1,i2,
2509      &     srt,pcx,pcy,pcz,nchrg,iblock)
2510 *
2511 * Process: PI + sigma(or Lambda) -> Kbar + N
2512 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2513 ******************************************
2514
2515 c NOTE: for PI + Hyperon: the produced kaons have mass 0.498
2516       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2517       PARAMETER      (AKA=0.498)
2518       COMMON   /AA/  R(3,MAXSTR)
2519 cc      SAVE /AA/
2520       COMMON   /BB/  P(3,MAXSTR)
2521 cc      SAVE /BB/
2522       COMMON   /CC/  E(MAXSTR)
2523 cc      SAVE /CC/
2524       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2525 cc      SAVE /EE/
2526       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2527 cc      SAVE /BG/
2528       COMMON   /NN/NNN
2529 cc      SAVE /NN/
2530       COMMON   /RUN/NUM
2531 cc      SAVE /RUN/
2532       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2533 cc      SAVE /PA/
2534       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2535 cc      SAVE /PB/
2536       COMMON   /PC/EPION(MAXSTR,MAXR)
2537 cc      SAVE /PC/
2538       COMMON   /PD/LPION(MAXSTR,MAXR)
2539 cc      SAVE /PD/
2540       dimension p1(4),p2(4)
2541       COMMON/RNDF77/NSEED
2542 cc      SAVE /RNDF77/
2543       SAVE
2544         px1cm=pcx
2545         py1cm=pcy
2546         pz1cm=pcz
2547         ictrl = 1
2548 csp06/07/01
2549         if(ielstc .eq. 1) then
2550 *    L/Si + meson -> L/Si + meson
2551              k1=i1
2552              k2=i2
2553            dm3=e(k1)
2554            dm4=e(k2)
2555            iblock = 10
2556         else
2557            iblock = 12
2558 csp06/07/01 end
2559 c        PI + Sigma(or Lambda) -> Kbar + N
2560         k1=i1
2561         k2=i2
2562 c        k1 must be bayron! So if I1 is PI, exchange k1 & k2.
2563         if(lb(i1).lt.14.or.lb(i1).gt.17) then
2564            k1=i2
2565            k2=i1
2566         endif
2567 cbz3/8/99 neutralk
2568         LB(K1) = 1 + int(2*RANART(NSEED))
2569         if(nchrg.eq.-2) lb(k1)=6
2570 c     if(nchrg.eq.-1) lb(k1)=2
2571 c     if(nchrg.eq. 0) lb(k1)=1
2572 c     if(nchrg.eq. 1) lb(k1)=9
2573         IF (NCHRG .EQ. 2) LB(K1) = 9
2574 cbz3/8/99 neutralk end
2575
2576 c     K-
2577         lb(k2)=21
2578         dm3=0.938
2579         if(nchrg.eq.-2.or.nchrg.eq.1) dm3=1.232
2580         dm4=aka
2581 c        dm3,dm4: the mass of final state particles.
2582          endif
2583
2584 ********Now, antikaon will be created.
2585 c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2586 c        pkmax: the maximum momentum of anti-kaon
2587         pkmax=sqrt((srt**2-(dm3+dm4)**2)*(srt**2-(dm3-dm4)**2))
2588      &         /2./srt
2589         pk=pkmax
2590 c-----------------------------------------------------
2591         css=1.-2.*RANART(NSEED)
2592         sss=sqrt(1.-css**2)
2593         fai=2*3.1415926*RANART(NSEED)
2594         p1(1)=pk*sss*cos(fai)
2595         p1(2)=pk*sss*sin(fai)
2596         p1(3)=pk*css
2597         do 1001 i=1,3
2598            p2(i)=-1.*p1(i)
2599  1001   continue
2600 c        p1,p2: the momenta of the two particles in their cms
2601 c        p1: momentum of kaon
2602 c        p2: momentum of Kbar
2603
2604 * Rotate the momenta of particles in the cms of I1 & I2
2605 clin-10/28/02 get rid of argument usage mismatch in rotate():
2606         pxrota=p1(1)
2607         pyrota=p1(2)
2608         pzrota=p1(3)
2609 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2610         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2611         p1(1)=pxrota
2612         p1(2)=pyrota
2613         p1(3)=pzrota
2614 c
2615         pxrota=p2(1)
2616         pyrota=p2(2)
2617         pzrota=p2(3)
2618 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2619         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2620         p2(1)=pxrota
2621         p2(2)=pyrota
2622         p2(3)=pzrota
2623 clin-10/28/02-end
2624
2625 * Find the momenta of particles in the final state in the nucleus_nucleus
2626 * cms frame.   Lorentz transformation into lab frame.
2627         e1cm   = sqrt(dm3**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2628         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2629         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2630         pt1i1 = betax*transf + p1(1)
2631         pt2i1 = betay*transf + p1(2)
2632         pt3i1 = betaz*transf + p1(3)
2633         eti1  = dm3
2634         lb1   = lb(k1)
2635
2636 * For second kaon, same
2637         e2cm   = sqrt(dm4**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2638         p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2639         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2640         pt1i2 = betax*transf + p2(1)
2641         pt2i2 = betay*transf + p2(2)
2642         pt3i2 = betaz*transf + p2(3)
2643 cbz3/2/99
2644 c        write(400,*)'3 ', pt1i2, pt2i2, pt3i2, dt*nt, srt
2645 cbz3/2/99end
2646 c        write(430,*)pt1i2, pt2i2, pt3i2, dt*nt, srt
2647         eti2  = dm4
2648         lb2   = lb(k2)
2649
2650 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2651 c        k1=i1
2652 c        k2=i2
2653 *       k1 stand for nucleon, k2 stand for kaon.
2654                 p(1,k1)=pt1i1
2655                 p(2,k1)=pt2i1
2656                 p(3,k1)=pt3i1
2657                 e(k1)=eti1
2658                 lb(k1)=lb1
2659                 p(1,k2)=pt1i2
2660                 p(2,k2)=pt2i2
2661                 p(3,k2)=pt3i2
2662                 e(k2)=eti2
2663                 lb(k2)=lb2
2664
2665 cc                iblock = 101  ! K(+)K(-) production
2666 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2667         return
2668         end
2669
2670 c-----------------------------------------------------------------------
2671
2672 c.....extracted from G. Song's ART expasion including K- interactions
2673 c.....file `KAONN.FOR'
2674
2675 ****************************************
2676         subroutine kaonN(brel,brsgm,irun,iseed,dt,nt,
2677      &     ictrl,i1,i2,iblock,srt,pcx,pcy,pcz,nchrg)
2678 *
2679 * Process: PI + sigma(or Lambda) <- Kbar + N
2680 * NOTE: the mass of K is assumed to be same as K0. ie. 0.498 NOT 0.494
2681 ****************************************
2682       PARAMETER      (MAXSTR=150001,MAXR=1,PI=3.1415926)
2683       PARAMETER      (AKA=0.498,ALA=1.1157,ASA=1.1974)
2684       COMMON   /AA/  R(3,MAXSTR)
2685 cc      SAVE /AA/
2686       COMMON   /BB/  P(3,MAXSTR)
2687 cc      SAVE /BB/
2688       COMMON   /CC/  E(MAXSTR)
2689 cc      SAVE /CC/
2690       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
2691 cc      SAVE /EE/
2692       COMMON   /BG/BETAX,BETAY,BETAZ,GAMMA
2693 cc      SAVE /BG/
2694       COMMON   /NN/NNN
2695 cc      SAVE /NN/
2696       COMMON   /RUN/NUM
2697 cc      SAVE /RUN/
2698       COMMON   /PA/RPION(3,MAXSTR,MAXR)
2699 cc      SAVE /PA/
2700       COMMON   /PB/PPION(3,MAXSTR,MAXR)
2701 cc      SAVE /PB/
2702       COMMON   /PC/EPION(MAXSTR,MAXR)
2703 cc      SAVE /PC/
2704       COMMON   /PD/LPION(MAXSTR,MAXR)
2705 cc      SAVE /PD/
2706       dimension p1(4),p2(4)
2707       COMMON/RNDF77/NSEED
2708 cc      SAVE /RNDF77/
2709       SAVE
2710         px1cm=pcx
2711         py1cm=pcy
2712         pz1cm=pcz
2713         ictrl = 1
2714 c        ratio: used for isospin decision.
2715         k1=i1
2716         k2=i2
2717 c        k1 must be bayron! So if I1 is Kaon, exchange k1 & k2.
2718         if(e(i1).lt.0.5.and.e(i1).gt.0.01) then
2719            k1=i2
2720            k2=i1
2721         endif
2722 *** note: for print out only *******************************
2723 c     record kaon's mass
2724         eee=e(k2)
2725 *** end **************
2726         rrr=RANART(NSEED)
2727         if(rrr.lt.brel) then
2728 c       Kbar + N -> Kbar + N
2729            lb1=lb(k1)
2730            lb2=lb(k2)
2731            em1=e(k1)
2732            em2=e(k2)
2733            iblock = 10
2734         else
2735            iblock = 12
2736         if(rrr.lt.(brel+brsgm)) then
2737 c        nchrg: Net charges of the two incoming particles.
2738 c           Kbar + N -> Sigma + PI
2739            em1=asa
2740            em2=0.138
2741
2742 cbz3/8/99 neutralk
2743            LB1 = 15 + int(3*RANART(NSEED))
2744            LB2 = 3 + int(3*RANART(NSEED))
2745         else
2746 c           Kbar + N -> Lambda + PI
2747            em1=ala
2748            em2=0.138
2749 c     LAmbda
2750            lb1=14
2751 cbz3/8/99 neutralk
2752            LB2 = 3 + int(3*RANART(NSEED))
2753 c           if(nchrg.eq.1)  lb2=5  ! K- + D++ -> Lambda + PI+
2754 c           if(nchrg.eq.0)  lb2=4  ! K- + p(D+,N*+) -> Lambda + PI0
2755 c          if(nchrg.eq.-1) lb2=3 ! K- + n(D,N*) -> Lambda + PI-
2756 cbz3/8/99 neutralk
2757
2758         endif
2759         endif
2760         lb(k1)=lb1
2761         lb(k2)=lb2
2762
2763 ********Now, antikaon will be created.
2764 c        call antikaon_fstate(iseed,srt,dm1,dm2,dm3,dm4,px,py,pz,icou1)
2765 c        pkmax: the maximum momentum of anti-kaon
2766 c        write(63,*)'srt,em1,em2',srt,em1,em2
2767 c        write(63,*)'-srt,em1,em2',srt,em1,em2
2768         pkmax=sqrt((srt**2-(em1+em2)**2)*(srt**2-(em1-em2)**2))
2769      &         /2./srt
2770         pk=pkmax
2771 c-----------------------------------------------------
2772         css=1.-2.*RANART(NSEED)
2773         sss=sqrt(1.-css**2)
2774         fai=2*3.1415926*RANART(NSEED)
2775         p1(1)=pk*sss*cos(fai)
2776         p1(2)=pk*sss*sin(fai)
2777         p1(3)=pk*css
2778         do 1001 i=1,3
2779            p2(i)=-1.*p1(i)
2780  1001   continue
2781 c        p1,p2: the momenta of the two particles in their cms
2782 c        p1: momentum of kaon
2783 c        p2: momentum of Kbar
2784
2785 * Rotate the momenta of particles in the cms of I1 & I2
2786
2787 clin-10/28/02 get rid of argument usage mismatch in rotate():
2788         pxrota=p1(1)
2789         pyrota=p1(2)
2790         pzrota=p1(3)
2791 c        call rotate(pcx,pcy,pcz,p1(1),p1(2),p1(3))
2792         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2793         p1(1)=pxrota
2794         p1(2)=pyrota
2795         p1(3)=pzrota
2796 c
2797         pxrota=p2(1)
2798         pyrota=p2(2)
2799         pzrota=p2(3)
2800 c        call rotate(pcx,pcy,pcz,p2(1),p2(2),p2(3))
2801         call rotate(pcx,pcy,pcz,pxrota,pyrota,pzrota)
2802         p2(1)=pxrota
2803         p2(2)=pyrota
2804         p2(3)=pzrota
2805 clin-10/28/02-end
2806
2807 * Find the momenta of particles in the final state in the nucleus_nucleus
2808 * cms frame.   Lorentz transformation into lab frame.
2809         e1cm   = sqrt(em1**2 + p1(1)**2 + p1(2)**2 + p1(3)**2)
2810         p1beta = p1(1)*betax + p1(2)*betay + p1(3)*betaz
2811         transf = gamma * ( gamma*p1beta / (gamma+1) + e1cm)
2812         pt1i1 = betax*transf + p1(1)
2813         pt2i1 = betay*transf + p1(2)
2814         pt3i1 = betaz*transf + p1(3)
2815         eti1  = em1
2816
2817 * For second kaon, same
2818         e2cm   = sqrt(em2**2 + p2(1)**2 + p2(2)**2 + p2(3)**2)
2819         p2beta = p2(1)*betax + p2(2)*betay + p2(3)*betaz
2820         transf = gamma * ( gamma*p2beta / (gamma+1) + e2cm)
2821         pt1i2 = betax*transf + p2(1)
2822         pt2i2 = betay*transf + p2(2)
2823         pt3i2 = betaz*transf + p2(3)
2824         eti2  = em2
2825
2826 c        if((pt1i1*px1+pt2i1*py1+pt3i1*pz1).gt.0.)then
2827 c        k1=i1
2828 c        k2=i2
2829 *       k1 stand for bayron, k2 stand for meson.
2830                 p(1,k1)=pt1i1
2831                 p(2,k1)=pt2i1
2832                 p(3,k1)=pt3i1
2833                 e(k1)=eti1
2834                 p(1,k2)=pt1i2
2835                 p(2,k2)=pt2i2
2836                 p(3,k2)=pt3i2
2837                 e(k2)=eti2
2838
2839 cc                iblock = 101  ! K(+)K(-) production
2840 * Get Kaons' momenta and coordinates in nucleus-nucleus cms. frame.
2841         return
2842         end
2843
2844 c=======================================================================
2845
2846 clin Below is the previous artana.f:
2847 c=======================================================================
2848
2849 c.....analysis subroutine before the hadronic space-time evolution
2850
2851       SUBROUTINE ARTAN1
2852       PARAMETER (MAXSTR=150001, MAXR=1)
2853 c.....y cut for mt spectrum
2854 cbz3/17/99
2855 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
2856       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
2857 cbz3/17/99 end
2858 c.....bin width for mt spectrum and y spectrum
2859 c.... clin-9/26/03 no symmetrization in y (or eta) for ana dat files
2860 c      PARAMETER (BMT = 0.05, BY = 0.2)
2861       PARAMETER (BMT = 0.05, BY = 0.4)
2862       COMMON /RUN/ NUM
2863 cc      SAVE /RUN/
2864       COMMON /ARERC1/MULTI1(MAXR)
2865 cc      SAVE /ARERC1/
2866       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2867      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2868      &     FT1(MAXSTR, MAXR),
2869      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2870      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2871 cbz3/17/99
2872 c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
2873 cc      SAVE /ARPRC1/
2874       COMMON /ARANA1/
2875      &     dy1ntb(50), dy1ntp(50), DY1HM(50),
2876      &     DY1KP(50), DY1KM(50), DY1K0S(50),
2877      &     DY1LA(50), DY1LB(50), DY1PHI(50),
2878      &     dm1pip(50), dm1pim(50), DMT1PR(50),
2879      &     DMT1PB(50), DMT1KP(50), dm1km(50),
2880      &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
2881      &     dy1msn(50), DY1PIP(50), DY1PIM(50),
2882      &     DY1PI0(50), DY1PR(50), DY1PB(50)
2883      &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
2884 cc      SAVE /ARANA1/
2885       SAVE
2886
2887 cbz3/17/99 end
2888       DO 1002 J = 1, NUM
2889          DO 1001 I = 1, MULTI1(J)
2890             ITYP = ITYP1(I, J)
2891             PX = PX1(I, J)
2892             PY = PY1(I, J)
2893             PZ = PZ1(I, J)
2894             EE = EE1(I, J)
2895             XM = XM1(I, J)
2896 c     2/24/03 leptons and photons:
2897             if(xm.lt.0.01) goto 200
2898             ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
2899 clin-9/2012 determine pseudo-rapidity more generally:
2900 c            eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
2901             if((PX**2+PY**2).gt.0.) then
2902                eta=asinh(PZ/sqrt(PX**2+PY**2))
2903             else
2904                eta = 1000000.0*sign(1.,PZ)
2905 clin-2/2013 for spectator target nucleons in LAB frame,
2906 c     note that finite precision of HBOOST
2907 c     would give spectator target nucleons a small but non-zero pz:
2908                if(abs(pz).le.1e-3) eta=0.
2909             endif
2910
2911             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
2912             DXMT = XMT - XM
2913 clin-9/2012 determine rapidity more generally:
2914 c            IF (ABS(PZ) .GE. EE) THEN
2915 c               PRINT *, 'IN ARTAN1'
2916 c               PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
2917 ccbzdbg2/16/99
2918 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
2919 ccbzdbg2/16/99
2920 ccbzdbg2/15/99
2921 c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
2922 ccbzdbg2/16/99
2923 c               PRINT *, ' XM = ', XM
2924 ccbzdbg2/16/99end
2925 c               GOTO 200
2926 c            else
2927 cc            Y = 0.5 * LOG((EE + PZ) / (EE - PZ))
2928 c               Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ +1e-5))
2929 cc               STOP
2930 ccbzdbg2/15/99end
2931 c            END IF
2932             if(XMT.gt.0.) then
2933                Y=asinh(PZ/XMT)
2934             else
2935                PRINT *, ' IN ARTAN1 mt=0'
2936                Y = 1000000.0*sign(1.,PZ)
2937             endif
2938
2939 c.....rapidity cut for the rapidity distribution
2940             IF (ABS(Y) .GE. 10.0) GOTO 100
2941 c.....clin-9/26/03 no symmetrization in y (or eta) for ana dat files
2942 c            IY = 1 + int(ABS(Y) / BY)
2943 c            Ieta = 1 + int(ABS(eta) / BY)
2944             IF (ABS(eta) .GE. 10.0) GOTO 100
2945             IY = 1 + int((Y+10.) / BY)
2946             Ieta = 1 + int((eta+10.) / BY)
2947
2948             IF (ITYP .LT. -1000) THEN
2949                dy1ntb(IY) = dy1ntb(IY) - 1.0
2950             END IF
2951             IF (ITYP .GT. 1000) THEN
2952                dy1ntb(IY) = dy1ntb(IY) + 1.0
2953             END IF
2954             IF (ITYP .EQ. -2212) THEN
2955                dy1ntp(IY) = dy1ntp(IY) - 1.0
2956             END IF
2957             IF (ITYP .EQ. 2212) THEN
2958                dy1ntp(IY) = dy1ntp(IY) + 1.0
2959             END IF
2960 c            IF (ITYP .EQ. -211 .OR. ITYP .EQ. -321 .OR.
2961 c     &         ITYP .EQ. -2212) THEN
2962             IF (ITYP .EQ. -2112) THEN
2963                DY1HM(IY) = DY1HM(IY) + 1.0
2964             END IF
2965 c
2966             IF (LUCHGE(ITYP).ne.0) THEN
2967                DY1CH(IY) = DY1CH(IY) + 1.0
2968                DE1CH(Ieta) = DE1CH(Ieta) + 1.0
2969                IF (LUCHGE(ITYP).lt.0) THEN
2970                   DY1NEG(IY) = DY1NEG(IY) + 1.0
2971                   DE1NEG(Ieta) = DE1NEG(Ieta) + 1.0
2972                endif
2973             END IF
2974
2975 cbz3/17/99
2976             IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR.
2977      &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
2978                dy1msn(IY) = dy1msn(IY) + 1.0
2979             END IF
2980             IF (ITYP .EQ. 211) THEN
2981                DY1PIP(IY) = DY1PIP(IY) + 1.0
2982             END IF
2983             IF (ITYP .EQ. -211) THEN
2984                DY1PIM(IY) = DY1PIM(IY) + 1.0
2985             END IF
2986             IF (ITYP .EQ. 111) THEN
2987                DY1PI0(IY) = DY1PI0(IY) + 1.0
2988             END IF
2989             IF (ITYP .EQ. 2212) THEN
2990                DY1PR(IY) = DY1PR(IY) + 1.0
2991             END IF
2992             IF (ITYP .EQ. -2212) THEN
2993                DY1PB(IY) = DY1PB(IY) + 1.0
2994             END IF
2995 cbz3/17/99 end
2996             IF (ITYP .EQ. 321) THEN
2997                DY1KP(IY) = DY1KP(IY) + 1.0
2998             END IF
2999             IF (ITYP .EQ. -321) THEN
3000                DY1KM(IY) = DY1KM(IY) + 1.0
3001             END IF
3002 clin-4/24/03 evaluate K0L instead of K0S, since sometimes we may decay K0S:
3003 c            IF (ITYP .EQ. 310) THEN
3004             IF (ITYP .EQ. 130) THEN
3005                DY1K0S(IY) = DY1K0S(IY) + 1.0
3006             END IF
3007             IF (ITYP .EQ. 3122) THEN
3008                DY1LA(IY) = DY1LA(IY) + 1.0
3009             END IF
3010             IF (ITYP .EQ. -3122) THEN
3011                DY1LB(IY) = DY1LB(IY) + 1.0
3012             END IF
3013             IF (ITYP .EQ. 333) THEN
3014                DY1PHI(IY) = DY1PHI(IY) + 1.0
3015             END IF
3016
3017 c.....insert rapidity cut for mt spectrum here
3018  100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
3019             IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
3020             IMT = 1 + int(DXMT / BMT)
3021             IF (ITYP .EQ. 211) THEN
3022                dm1pip(IMT) = dm1pip(IMT) + 1.0 / XMT
3023             END IF
3024             IF (ITYP .EQ. -211) THEN
3025                dm1pim(IMT) = dm1pim(IMT) +
3026      &            1.0 / XMT
3027             END IF
3028             IF (ITYP .EQ. 2212) THEN
3029                DMT1PR(IMT) = DMT1PR(IMT) + 1.0 / XMT
3030             END IF
3031             IF (ITYP .EQ. -2212) THEN
3032                DMT1PB(IMT) = DMT1PB(IMT) + 1.0 / XMT
3033             END IF
3034             IF (ITYP .EQ. 321) THEN
3035                DMT1KP(IMT) = DMT1KP(IMT) + 1.0 / XMT
3036             END IF
3037             IF (ITYP .EQ. -321) THEN
3038                dm1km(IMT) = dm1km(IMT) + 1.0 / XMT
3039             END IF
3040 clin-4/24/03:
3041 c            IF (ITYP .EQ. 310) THEN
3042             IF (ITYP .EQ. 130) THEN
3043                dm1k0s(IMT) = dm1k0s(IMT) + 1.0 / XMT
3044             END IF
3045             IF (ITYP .EQ. 3122) THEN
3046                DMT1LA(IMT) = DMT1LA(IMT) + 1.0 / XMT
3047             END IF
3048             IF (ITYP .EQ. -3122) THEN
3049                DMT1LB(IMT) = DMT1LB(IMT) + 1.0 / XMT
3050             END IF
3051
3052  200        CONTINUE
3053  1001    CONTINUE
3054  1002 CONTINUE
3055
3056       RETURN
3057       END
3058
3059 c-----------------------------------------------------------------------
3060
3061 c.....analysis subroutine after the hadronic space-time evolution
3062
3063       SUBROUTINE ARTAN2
3064
3065       PARAMETER (MAXSTR=150001, MAXR=1)
3066 c.....y cut for mt spectrum
3067 cbz3/17/99
3068 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3069       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3070 cbz3/17/99 end
3071 c.....bin width for mt spectrum and y spectrum
3072 c      PARAMETER (BMT = 0.05, BY = 0.2)
3073       PARAMETER (BMT = 0.05, BY = 0.4)
3074       COMMON /RUN/ NUM
3075 cc      SAVE /RUN/
3076       COMMON /ARERC1/MULTI1(MAXR)
3077 cc      SAVE /ARERC1/
3078       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
3079      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
3080      &     FT1(MAXSTR, MAXR),
3081      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
3082      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
3083 cbz3/17/99
3084 c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
3085 cc      SAVE /ARPRC1/
3086       COMMON /ARANA2/
3087      &     dy2ntb(50), dy2ntp(50), DY2HM(50),
3088      &     DY2KP(50), DY2KM(50), DY2K0S(50),
3089      &     DY2LA(50), DY2LB(50), DY2PHI(50),
3090      &     dm2pip(50), dm2pim(50), DMT2PR(50),
3091      &     DMT2PB(50), DMT2KP(50), dm2km(50),
3092      &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
3093      &     dy2msn(50), DY2PIP(50), DY2PIM(50),
3094      &     DY2PI0(50), DY2PR(50), DY2PB(50)
3095      &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
3096 cbz3/17/99 end
3097 cc      SAVE /ARANA2/
3098       SAVE
3099
3100       DO 1002 J = 1, NUM
3101          DO 1001 I = 1, MULTI1(J)
3102             ITYP = ITYP1(I, J)
3103             PX = PX1(I, J)
3104             PY = PY1(I, J)
3105             PZ = PZ1(I, J)
3106             EE = EE1(I, J)
3107             XM = XM1(I, J)
3108             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
3109 c     2/24/03 leptons and photons:
3110             if(xm.lt.0.01) goto 200
3111             DXMT = XMT - XM
3112             ptot = sqrt(PX ** 2 + PY ** 2 + pz ** 2)
3113 clin-9/2012 determine pseudo-rapidity more generally:
3114 c            eta = 0.5*alog((Ptot+pz+1e-5)/(ptot-pz+1e-5))
3115             if((PX**2+PY**2).gt.0.) then
3116                eta=asinh(PZ/sqrt(PX**2+PY**2))
3117             else
3118                eta = 1000000.0*sign(1.,PZ)
3119 clin-2/2013 for spectator target nucleons in LAB frame:
3120                if(abs(pz).le.1e-3) eta=0.
3121             endif
3122
3123 clin-9/2012 determine rapidity more generally:
3124 c            IF (ABS(PZ) .GE. EE) THEN
3125 c               PRINT *, 'IN ARTAN2'
3126 c               PRINT *, 'PARTICLE ', I, ' RUN ', J, 'PREC ERR'
3127 ccbzdbg2/16/99
3128 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3129 ccbzdbg2/16/99
3130 ccbzdbg2/15/99
3131 c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
3132 ccbzdbg2/16/99
3133 c               PRINT *, ' XM = ', XM
3134 ccbzdbg2/16/99end
3135 c               GOTO 200
3136 cc               STOP
3137 ccbzdbg2/15/99end
3138 c            END IF
3139 c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
3140             if(XMT.gt.0.) then
3141                Y=asinh(PZ/XMT)
3142             else
3143                PRINT *, ' IN ARTAN2 mt=0'
3144                Y = 1000000.0*sign(1.,PZ)
3145             endif
3146
3147 c.....rapidity cut for the rapidity distribution
3148             IF (ABS(Y) .GE. 10.0) GOTO 100
3149 c            IY = 1 + int(ABS(Y) / BY)
3150 c            Ieta = 1 + int(ABS(eta) / BY)
3151             IF (ABS(eta) .GE. 10.0) GOTO 100
3152             IY = 1 + int((Y+10.) / BY)
3153             Ieta = 1 + int((eta+10.) / BY)
3154
3155             IF (ITYP .LT. -1000) THEN
3156                dy2ntb(IY) = dy2ntb(IY) - 1.0
3157             END IF
3158             IF (ITYP .GT. 1000) THEN
3159                dy2ntb(IY) = dy2ntb(IY) + 1.0
3160             END IF
3161             IF (ITYP .EQ. -2212) THEN
3162                dy2ntp(IY) = dy2ntp(IY) - 1.0
3163             END IF
3164             IF (ITYP .EQ. 2212) THEN
3165                dy2ntp(IY) = dy2ntp(IY) + 1.0
3166             END IF
3167             IF (ITYP .EQ. -2112) THEN
3168                DY2HM(IY) = DY2HM(IY) + 1.0
3169             END IF
3170
3171             IF (LUCHGE(ITYP).ne.0) THEN
3172                DY2CH(IY) = DY2CH(IY) + 1.0
3173                DE2CH(Ieta) = DE2CH(Ieta) + 1.0
3174                IF (LUCHGE(ITYP).lt.0) THEN
3175                   DY2NEG(IY) = DY2NEG(IY) + 1.0
3176                   DE2NEG(Ieta) = DE2NEG(Ieta) + 1.0
3177                endif
3178             END IF
3179
3180 cbz3/17/99
3181             IF ((ITYP .GE. 100 .AND. ITYP .LT. 1000) .OR.
3182      &         (ITYP .GT. -1000 .AND. ITYP .LE. -100)) THEN
3183                dy2msn(IY) = dy2msn(IY) + 1.0
3184             END IF
3185             IF (ITYP .EQ. 211) THEN
3186                DY2PIP(IY) = DY2PIP(IY) + 1.0
3187             END IF
3188             IF (ITYP .EQ. -211) THEN
3189                DY2PIM(IY) = DY2PIM(IY) + 1.0
3190             END IF
3191             IF (ITYP .EQ. 111) THEN
3192                DY2PI0(IY) = DY2PI0(IY) + 1.0
3193             END IF
3194             IF (ITYP .EQ. 2212) THEN
3195                DY2PR(IY) = DY2PR(IY) + 1.0
3196             END IF
3197             IF (ITYP .EQ. -2212) THEN
3198                DY2PB(IY) = DY2PB(IY) + 1.0
3199             END IF
3200 cbz3/17/99 end
3201             IF (ITYP .EQ. 321) THEN
3202                DY2KP(IY) = DY2KP(IY) + 1.0
3203             END IF
3204             IF (ITYP .EQ. -321) THEN
3205                DY2KM(IY) = DY2KM(IY) + 1.0
3206             END IF
3207 clin-4/24/03:
3208 c            IF (ITYP .EQ. 310) THEN
3209             IF (ITYP .EQ. 130) THEN
3210                DY2K0S(IY) = DY2K0S(IY) + 1.0
3211             END IF
3212             IF (ITYP .EQ. 3122) THEN
3213                DY2LA(IY) = DY2LA(IY) + 1.0
3214             END IF
3215             IF (ITYP .EQ. -3122) THEN
3216                DY2LB(IY) = DY2LB(IY) + 1.0
3217             END IF
3218             IF (ITYP .EQ. 333) THEN
3219                DY2PHI(IY) = DY2PHI(IY) + 1.0
3220             END IF
3221
3222 c.....insert rapidity cut for mt spectrum here
3223  100        IF (Y .LT. YMT1 .OR. Y .GT. YMT2) GOTO 200
3224             IF (DXMT .GE. 50.0 * BMT .OR. DXMT .EQ. 0) GOTO 200
3225             IMT = 1 + int(DXMT / BMT)
3226             IF (ITYP .EQ. 211) THEN
3227                dm2pip(IMT) = dm2pip(IMT) + 1.0 / XMT
3228             END IF
3229             IF (ITYP .EQ. -211) THEN
3230                dm2pim(IMT) = dm2pim(IMT) +
3231      &            1.0 / XMT
3232             END IF
3233             IF (ITYP .EQ. 2212) THEN
3234                DMT2PR(IMT) = DMT2PR(IMT) + 1.0 / XMT
3235             END IF
3236             IF (ITYP .EQ. -2212) THEN
3237                DMT2PB(IMT) = DMT2PB(IMT) + 1.0 / XMT
3238             END IF
3239             IF (ITYP .EQ. 321) THEN
3240                DMT2KP(IMT) = DMT2KP(IMT) + 1.0 / XMT
3241             END IF
3242             IF (ITYP .EQ. -321) THEN
3243                dm2km(IMT) = dm2km(IMT) + 1.0 / XMT
3244             END IF
3245 clin-4/24/03:
3246 c            IF (ITYP .EQ. 310) THEN
3247             IF (ITYP .EQ. 130) THEN
3248                dm2k0s(IMT) = dm2k0s(IMT) + 1.0 / XMT
3249             END IF
3250             IF (ITYP .EQ. 3122) THEN
3251                DMT2LA(IMT) = DMT2LA(IMT) + 1.0 / XMT
3252             END IF
3253             IF (ITYP .EQ. -3122) THEN
3254                DMT2LB(IMT) = DMT2LB(IMT) + 1.0 / XMT
3255             END IF
3256
3257  200        CONTINUE
3258  1001    CONTINUE
3259  1002 CONTINUE
3260
3261       RETURN
3262       END
3263
3264 c-----------------------------------------------------------------------
3265
3266 c.....output analysis results at the end of the simulation
3267
3268       SUBROUTINE ARTOUT(NEVNT)
3269
3270       PARAMETER (MAXSTR=150001, MAXR=1)
3271 c.....y cut for mt spectrum
3272 cbz3/17/99
3273 c      PARAMETER (YMT1 = -0.4, YMT2 = 0.4)
3274       PARAMETER (YMT1 = -1.0, YMT2 = 1.0)
3275 cbz3/17/99 end
3276 c.....bin width for mt spectrum and y spectrum
3277 c      PARAMETER (BMT = 0.05, BY = 0.2)
3278       PARAMETER (BMT = 0.05, BY = 0.4)
3279       COMMON /RUN/ NUM
3280 cc      SAVE /RUN/
3281       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
3282      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
3283      &     FT1(MAXSTR, MAXR),
3284      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
3285      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
3286 cbz3/17/99
3287 c     &     dm1k0s(50), DMT1LA(50), DMT1LB(50)
3288 cc      SAVE /ARPRC1/
3289       COMMON /ARANA1/
3290      &     dy1ntb(50), dy1ntp(50), DY1HM(50),
3291      &     DY1KP(50), DY1KM(50), DY1K0S(50),
3292      &     DY1LA(50), DY1LB(50), DY1PHI(50),
3293      &     dm1pip(50), dm1pim(50), DMT1PR(50),
3294      &     DMT1PB(50), DMT1KP(50), dm1km(50),
3295      &     dm1k0s(50), DMT1LA(50), DMT1LB(50),
3296      &     dy1msn(50), DY1PIP(50), DY1PIM(50),
3297      &     DY1PI0(50), DY1PR(50), DY1PB(50)
3298      &     ,DY1NEG(50), DY1CH(50), DE1NEG(50), DE1CH(50)
3299 cbz3/17/99 end
3300 cc      SAVE /ARANA1/
3301 cbz3/17/99
3302 c     &     dm2k0s(50), DMT2LA(50), DMT2LB(50)
3303       COMMON /ARANA2/
3304      &     dy2ntb(50), dy2ntp(50), DY2HM(50),
3305      &     DY2KP(50), DY2KM(50), DY2K0S(50),
3306      &     DY2LA(50), DY2LB(50), DY2PHI(50),
3307      &     dm2pip(50), dm2pim(50), DMT2PR(50),
3308      &     DMT2PB(50), DMT2KP(50), dm2km(50),
3309      &     dm2k0s(50), DMT2LA(50), DMT2LB(50),
3310      &     dy2msn(50), DY2PIP(50), DY2PIM(50),
3311      &     DY2PI0(50), DY2PR(50), DY2PB(50)
3312      &     ,DY2NEG(50), DY2CH(50), DE2NEG(50), DE2CH(50)
3313 cc      SAVE /ARANA2/
3314       SAVE
3315 cbz3/17/99 end
3316 cms      OPEN (30, FILE = 'ana/dndy_netb.dat', STATUS = 'UNKNOWN')
3317 cms      OPEN (31, FILE = 'ana/dndy_netp.dat', STATUS = 'UNKNOWN')
3318 cms      OPEN (32, FILE = 'ana/dndy_nb.dat', STATUS = 'UNKNOWN')
3319 cms      OPEN (33, FILE = 'ana/dndy_neg.dat', STATUS = 'UNKNOWN')
3320 cms      OPEN (34, FILE = 'ana/dndy_ch.dat', STATUS = 'UNKNOWN')
3321 cms      OPEN (35, FILE = 'ana/dnde_neg.dat', STATUS = 'UNKNOWN')
3322 cms      OPEN (36, FILE = 'ana/dnde_ch.dat', STATUS = 'UNKNOWN')
3323 cms      OPEN (37, FILE = 'ana/dndy_kp.dat', STATUS = 'UNKNOWN')
3324 cms      OPEN (38, FILE = 'ana/dndy_km.dat', STATUS = 'UNKNOWN')
3325 clin-4/24/03
3326 c      OPEN (39, FILE = 'ana/dndy_k0s.dat', STATUS = 'UNKNOWN')
3327 cms      OPEN (39, FILE = 'ana/dndy_k0l.dat', STATUS = 'UNKNOWN')
3328 cms      OPEN (40, FILE = 'ana/dndy_la.dat', STATUS = 'UNKNOWN')
3329 cms      OPEN (41, FILE = 'ana/dndy_lb.dat', STATUS = 'UNKNOWN')
3330 cms      OPEN (42, FILE = 'ana/dndy_phi.dat', STATUS = 'UNKNOWN')
3331 cbz3/17/99
3332 cms      OPEN (43, FILE = 'ana/dndy_meson.dat', STATUS = 'UNKNOWN')
3333 cms      OPEN (44, FILE = 'ana/dndy_pip.dat', STATUS = 'UNKNOWN')
3334 cms      OPEN (45, FILE = 'ana/dndy_pim.dat', STATUS = 'UNKNOWN')
3335 cms      OPEN (46, FILE = 'ana/dndy_pi0.dat', STATUS = 'UNKNOWN')
3336 cms      OPEN (47, FILE = 'ana/dndy_pr.dat', STATUS = 'UNKNOWN')
3337 cms      OPEN (48, FILE = 'ana/dndy_pb.dat', STATUS = 'UNKNOWN')
3338 cbz3/17/99 end
3339
3340 cms      OPEN (50, FILE = 'ana/dndmtdy_pip.dat', STATUS = 'UNKNOWN')
3341 cms      OPEN (51, FILE = 'ana/dndmtdy_0_1_pim.dat', STATUS = 'UNKNOWN')
3342 cms      OPEN (52, FILE = 'ana/dndmtdy_pr.dat', STATUS = 'UNKNOWN')
3343 cms      OPEN (53, FILE = 'ana/dndmtdy_pb.dat', STATUS = 'UNKNOWN')
3344 cms      OPEN (54, FILE = 'ana/dndmtdy_kp.dat', STATUS = 'UNKNOWN')
3345 cms      OPEN (55, FILE = 'ana/dndmtdy_0_5_km.dat', STATUS = 'UNKNOWN')
3346 cms      OPEN (56, FILE = 'ana/dndmtdy_k0s.dat', STATUS = 'UNKNOWN')
3347 cms      OPEN (57, FILE = 'ana/dndmtdy_la.dat', STATUS = 'UNKNOWN')
3348 cms      OPEN (58, FILE = 'ana/dndmtdy_lb.dat', STATUS = 'UNKNOWN')
3349 clin-9/26/03 no symmetrization in y (or eta) for ana dat files
3350 c      SCALE1 = 1. / REAL(NEVNT * NUM) / BY / 2.0
3351       SCALE1 = 1. / REAL(NEVNT * NUM) / BY
3352       SCALE2 = 1. / REAL(NEVNT * NUM) / BMT / (YMT2 - YMT1)
3353 c
3354       DO 1001 I = 1, 50
3355          ymid=-10.+BY * (I - 0.5)
3356 cms         WRITE (30, 333) ymid, SCALE1 * dy1ntb(I)
3357 cms         WRITE (31, 333) ymid, SCALE1 * dy1ntp(I)
3358 cms         WRITE (32, 333) ymid, SCALE1 * DY1HM(I)
3359 cms         WRITE (37, 333) ymid, SCALE1 * DY1KP(I)
3360 cms         WRITE (38, 333) ymid, SCALE1 * DY1KM(I)
3361 cms         WRITE (39, 333) ymid, SCALE1 * DY1K0S(I)
3362 cms         WRITE (40, 333) ymid, SCALE1 * DY1LA(I)
3363 cms         WRITE (41, 333) ymid, SCALE1 * DY1LB(I)
3364 cms         WRITE (42, 333) ymid, SCALE1 * DY1PHI(I)
3365 cms         WRITE (33, 333) ymid, SCALE1 * DY1NEG(I)
3366 cms         WRITE (34, 333) ymid, SCALE1 * DY1CH(I)
3367 cms         WRITE (35, 333) ymid, SCALE1 * DE1NEG(I)
3368 cms         WRITE (36, 333) ymid, SCALE1 * DE1CH(I)
3369 cms         WRITE (43, 333) ymid, SCALE1 * dy1msn(I)
3370 cms         WRITE (44, 333) ymid, SCALE1 * DY1PIP(I)
3371 cms         WRITE (45, 333) ymid, SCALE1 * DY1PIM(I)
3372 cms         WRITE (46, 333) ymid, SCALE1 * DY1PI0(I)
3373 cms         WRITE (47, 333) ymid, SCALE1 * DY1PR(I)
3374 cms         WRITE (48, 333) ymid, SCALE1 * DY1PB(I)
3375
3376          IF (dm1pip(I) .NE. 0.0) THEN
3377 c            WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm1pip(I)
3378          END IF
3379          IF (dm1pim(I) .NE. 0.0) THEN
3380 c            WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3381 c     &         dm1pim(I)
3382          END IF
3383          IF (DMT1PR(I) .NE. 0.0) THEN
3384 c            WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT1PR(I)
3385          END IF
3386          IF (DMT1PB(I) .NE. 0.0) THEN
3387 c            WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT1PB(I)
3388          END IF
3389          IF (DMT1KP(I) .NE. 0.0) THEN
3390 c            WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT1KP(I)
3391          END IF
3392          IF (dm1km(I) .NE. 0.0) THEN
3393 c            WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3394 c     &         dm1km(I)
3395          END IF
3396          IF (dm1k0s(I) .NE. 0.0) THEN
3397 c            WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm1k0s(I)
3398          END IF
3399          IF (DMT1LA(I) .NE. 0.0) THEN
3400 c            WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT1LA(I)
3401          END IF
3402          IF (DMT1LB(I) .NE. 0.0) THEN
3403 c            WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT1LB(I)
3404          END IF
3405  1001 CONTINUE
3406 c
3407       DO 1002 I = 30, 48
3408 c         WRITE (I, *) 'after hadron evolution'
3409  1002    CONTINUE
3410       DO 1003 I = 50, 58
3411 c         WRITE (I, *) 'after hadron evolution'
3412  1003 CONTINUE
3413
3414       DO 1004 I = 1, 50
3415          ymid=-10.+BY * (I - 0.5)
3416 cms         WRITE (30, 333) ymid, SCALE1 * dy2ntb(I)
3417 cms         WRITE (31, 333) ymid, SCALE1 * dy2ntp(I)
3418 cms         WRITE (32, 333) ymid, SCALE1 * DY2HM(I)
3419 cms         WRITE (37, 333) ymid, SCALE1 * DY2KP(I)
3420 cms         WRITE (38, 333) ymid, SCALE1 * DY2KM(I)
3421 cms         WRITE (39, 333) ymid, SCALE1 * DY2K0S(I)
3422 cms         WRITE (40, 333) ymid, SCALE1 * DY2LA(I)
3423 cms         WRITE (41, 333) ymid, SCALE1 * DY2LB(I)
3424 cms         WRITE (42, 333) ymid, SCALE1 * DY2PHI(I)
3425 cms         WRITE (33, 333) ymid, SCALE1 * DY2NEG(I)
3426 cms         WRITE (34, 333) ymid, SCALE1 * DY2CH(I)
3427 cms         WRITE (35, 333) ymid, SCALE1 * DE2NEG(I)
3428 cms         WRITE (36, 333) ymid, SCALE1 * DE2CH(I)
3429 cms         WRITE (43, 333) ymid, SCALE1 * dy2msn(I)
3430 cms         WRITE (44, 333) ymid, SCALE1 * DY2PIP(I)
3431 cms         WRITE (45, 333) ymid, SCALE1 * DY2PIM(I)
3432 cms         WRITE (46, 333) ymid, SCALE1 * DY2PI0(I)
3433 cms         WRITE (47, 333) ymid, SCALE1 * DY2PR(I)
3434 cms         WRITE (48, 333) ymid, SCALE1 * DY2PB(I)
3435 c
3436          IF (dm2pip(I) .NE. 0.0) THEN
3437 c            WRITE (50, 333) BMT * (I - 0.5), SCALE2 * dm2pip(I)
3438          END IF
3439          IF (dm2pim(I) .NE. 0.0) THEN
3440 c            WRITE (51, 333) BMT * (I - 0.5), SCALE2 * 0.1 *
3441 c     &         dm2pim(I)
3442          END IF
3443          IF (DMT2PR(I) .NE. 0.0) THEN
3444 c            WRITE (52, 333) BMT * (I - 0.5), SCALE2 * DMT2PR(I)
3445          END IF
3446          IF (DMT2PB(I) .NE. 0.0) THEN
3447 c            WRITE (53, 333) BMT * (I - 0.5), SCALE2 * DMT2PB(I)
3448          END IF
3449          IF (DMT2KP(I) .NE. 0.0) THEN
3450 c            WRITE (54, 333) BMT * (I - 0.5), SCALE2 * DMT2KP(I)
3451          END IF
3452          IF (dm2km(I) .NE. 0.0) THEN
3453 c            WRITE (55, 333) BMT * (I - 0.5), SCALE2 * 0.5 *
3454 c     &         dm2km(I)
3455          END IF
3456          IF (dm2k0s(I) .NE. 0.0) THEN
3457 c            WRITE (56, 333) BMT * (I - 0.5), SCALE2 * dm2k0s(I)
3458          END IF
3459          IF (DMT2LA(I) .NE. 0.0) THEN
3460 c            WRITE (57, 333) BMT * (I - 0.5), SCALE2 * DMT2LA(I)
3461          END IF
3462          IF (DMT2LB(I) .NE. 0.0) THEN
3463 c            WRITE (58, 333) BMT * (I - 0.5), SCALE2 * DMT2LB(I)
3464          END IF
3465  1004 CONTINUE
3466 c      format(2(f12.5,1x))
3467
3468       RETURN
3469       END
3470
3471 c-----------------------------------------------------------------------
3472
3473 c.....analysis subroutine in HIJING before parton cascade evolution
3474       SUBROUTINE HJANA1
3475
3476       PARAMETER (YMAX = 1.0, YMIN = -1.0)
3477       PARAMETER (DMT = 0.05, DY = 0.2)
3478       PARAMETER (DR = 0.2)
3479       PARAMETER (MAXPTN=400001,MAXSTR=150001)
3480       DIMENSION dyp1(50), DMYP1(200), DEYP1(50)
3481       DIMENSION dyg1(50), DMYG1(200), DEYG1(50)
3482       DIMENSION SNYP1(50), SMYP1(200), SEYP1(50)
3483       DIMENSION SNYG1(50), SMYG1(200), SEYG1(50)
3484       DIMENSION dnrpj1(50), dnrtg1(50), dnrin1(50),
3485      &   dnrtt1(50)
3486       DIMENSION dyg1c(50), dmyg1c(50), deyg1c(50)
3487       DIMENSION snrpj1(50), snrtg1(50), snrin1(50),
3488      &   snrtt1(50)
3489       DIMENSION snyg1c(50), smyg1c(50), seyg1c(50)
3490       DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
3491
3492       COMMON /PARA1/ MUL
3493 cc      SAVE /PARA1/
3494       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3495 cc      SAVE /HPARNT/
3496       COMMON/hjcrdn/YP(3,300),YT(3,300)
3497 cc      SAVE /hjcrdn/
3498       COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3499      &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3500      &   PJPM(300,500),NTJ(300),KFTJ(300,500),
3501      &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3502      &   PJTE(300,500),PJTM(300,500)
3503 cc      SAVE /HJJET1/
3504       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3505      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3506      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3507 cc      SAVE /HJJET2/
3508       COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
3509      &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
3510      &     XMASS0(MAXPTN), ITYP0(MAXPTN)
3511 cc      SAVE /prec1/
3512       COMMON /AREVT/ IAEVT, IARUN, MISS
3513 cc      SAVE /AREVT/
3514       COMMON /AROUT/ IOUT
3515 cc      SAVE /AROUT/
3516       SAVE
3517       DATA IW/0/
3518
3519       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3520          DO 1001 I = 1, 200
3521             DMYP1(I) = SMYP1(I)
3522             DMYG1(I) = SMYG1(I)
3523  1001    CONTINUE
3524
3525          DO 1002 I = 1, 50
3526             dyp1(I) = SNYP1(I)
3527             DEYP1(I) = SEYP1(I)
3528             dyg1(I) = SNYG1(I)
3529             DEYG1(I) = SEYG1(I)
3530             dnrpj1(I) = snrpj1(I)
3531             dnrtg1(I) = snrtg1(I)
3532             dnrin1(I) = snrin1(I)
3533             dnrtt1(I) = snrtt1(I)
3534             dyg1c(I) = snyg1c(I)
3535             dmyg1c(I) = smyg1c(I)
3536             deyg1c(I) = seyg1c(I)
3537  1002    CONTINUE
3538          nsubp = nsubpS
3539          nsubg = nsubgS
3540          NISG = NISGS
3541       ELSE
3542          DO 1003 I = 1, 200
3543             SMYP1(I) = DMYP1(I)
3544             SMYG1(I) = DMYG1(I)
3545  1003    CONTINUE
3546
3547          DO 1004 I = 1, 50
3548             SNYP1(I) = dyp1(I)
3549             SEYP1(I) = DEYP1(I)
3550             SNYG1(I) = dyg1(I)
3551             SEYG1(I) = DEYG1(I)
3552             snrpj1(I) = dnrpj1(I)
3553             snrtg1(I) = dnrtg1(I)
3554             snrin1(I) = dnrin1(I)
3555             snrtt1(I) = dnrtt1(I)
3556             snyg1c(I) = dyg1c(I)
3557             smyg1c(I) = dmyg1c(I)
3558             seyg1c(I) = deyg1c(I)
3559  1004    CONTINUE
3560          nsubpS = nsubp
3561          nsubgS = nsubg
3562          NISGS = NISG
3563          isevt = IAEVT
3564          isrun = IARUN
3565          IW = IW + 1
3566       END IF
3567 c.....analysis
3568       DO 1006 I = 1, IHNT2(1)
3569          DO 1005 J = 1, NPJ(I)
3570             ITYP = KFPJ(I, J)
3571             PX = PJPX(I, J)
3572             PY = PJPY(I, J)
3573             PZ = PJPZ(I, J)
3574             PE = PJPE(I, J)
3575             PM = PJPM(I, J)
3576             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
3577             DXMT = XMT - PM
3578 clin-9/2012 determine rapidity more generally:
3579 c            IF (ABS(PZ) .GE. PE) THEN
3580 c               PRINT *, ' IN HJANA1, PROJ STR ', I, ' PART ', J
3581 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3582 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
3583 c               PRINT *, ' XM = ', PM
3584 c               GOTO 200
3585 c            END IF
3586 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
3587             if(XMT.gt.0.) then
3588                RAP=asinh(PZ/XMT)
3589             else
3590                RAP = 1000000.0*sign(1.,PZ)
3591             endif
3592
3593             IY = 1 + int(ABS(RAP) / DY)
3594 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
3595 c            IF (IY .GT. 50) GOTO 100
3596             IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
3597             dyp1(IY) = dyp1(IY) + 1.0
3598             DEYP1(IY) = DEYP1(IY) + XMT
3599             IF (ITYP .EQ. 21) THEN
3600                dyg1(IY) = dyg1(IY) + 1.0
3601                DEYG1(IY) = DEYG1(IY) + XMT
3602             END IF
3603  100        CONTINUE
3604             IMT = 1 + int(DXMT / DMT)
3605             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200
3606             IF (IMT .GT. 200) GOTO 200
3607             DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
3608             IF (ITYP .EQ. 21) THEN
3609                DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
3610             END IF
3611  200        CONTINUE
3612  1005    CONTINUE
3613  1006 CONTINUE
3614
3615       DO 1008 I = 1, IHNT2(3)
3616          DO 1007 J = 1, NTJ(I)
3617             ITYP = KFTJ(I, J)
3618             PX = PJTX(I, J)
3619             PY = PJTY(I, J)
3620             PZ = PJTZ(I, J)
3621             PE = PJTE(I, J)
3622             PM = PJTM(I, J)
3623             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
3624             DXMT = XMT - PM
3625 clin-9/2012 determine rapidity more generally:
3626 c            IF (ABS(PZ) .GE. PE) THEN
3627 c               PRINT *, ' IN HJANA1, TARG STR ', I, ' PART ', J
3628 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3629 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
3630 c               PRINT *, ' XM = ', PM
3631 c               GOTO 400
3632 c            END IF
3633 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
3634             if(XMT.gt.0.) then
3635                RAP=asinh(PZ/XMT)
3636             else
3637                PRINT *, ' IN HJANA1 mt=0'
3638                RAP = 1000000.0*sign(1.,PZ)
3639             endif
3640
3641             IY = 1 + int(ABS(RAP) / DY)
3642 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
3643 c            IF (IY .GT. 50) GOTO 300
3644             IF (IY.lt.1 .or.IY .GT. 50) GOTO 300
3645             dyp1(IY) = dyp1(IY) + 1.0
3646             DEYP1(IY) = DEYP1(IY) + XMT
3647             IF (ITYP .EQ. 21) THEN
3648                dyg1(IY) = dyg1(IY) + 1.0
3649                DEYG1(IY) = DEYG1(IY) + XMT
3650             END IF
3651  300        CONTINUE
3652             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400
3653             IMT = 1 + int(DXMT / DMT)
3654             IF (IMT .GT. 200) GOTO 400
3655             DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
3656             IF (ITYP .EQ. 21) THEN
3657                DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
3658             END IF
3659  400        CONTINUE
3660  1007    CONTINUE
3661  1008 CONTINUE
3662
3663       DO 1010 I = 1, NSG
3664          DO 1009 J = 1, NJSG(I)
3665             ITYP = K2SG(I, J)
3666             PX = PXSG(I, J)
3667             PY = PYSG(I, J)
3668             PZ = PZSG(I, J)
3669             PE = PESG(I, J)
3670             PM = PMSG(I, J)
3671             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
3672             DXMT = XMT - PM
3673 clin-9/2012 determine rapidity more generally:
3674 c            IF (ABS(PZ) .GE. PE) THEN
3675 c               PRINT *, ' IN HJANA1, INDP STR ', I, ' PART ', J
3676 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3677 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
3678 c               PRINT *, ' XM = ', PM
3679 c               GOTO 600
3680 c            END IF
3681 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
3682             if(XMT.gt.0.) then
3683                RAP=asinh(PZ/XMT)
3684             else
3685                PRINT *, ' IN HJANA1 mt=0'
3686                RAP = 1000000.0*sign(1.,PZ)
3687             endif
3688
3689             IY = 1 + int(ABS(RAP) / DY)
3690 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
3691 c            IF (IY .GT. 50) GOTO 500
3692             IF (IY.lt.1 .or.IY .GT. 50) GOTO 500
3693             dyp1(IY) = dyp1(IY) + 1.0
3694             DEYP1(IY) = DEYP1(IY) + XMT
3695             IF (ITYP .EQ. 21) THEN
3696                dyg1(IY) = dyg1(IY) + 1.0
3697                DEYG1(IY) = DEYG1(IY) + XMT
3698             END IF
3699  500        CONTINUE
3700             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600
3701             IMT = 1 + int(DXMT / DMT)
3702             IF (IMT .GT. 200) GOTO 600
3703             DMYP1(IMT) = DMYP1(IMT) + 1.0 / XMT
3704             IF (ITYP .EQ. 21) THEN
3705                DMYG1(IMT) = DMYG1(IMT) + 1.0 / XMT
3706             END IF
3707  600        CONTINUE
3708  1009    CONTINUE
3709  1010 CONTINUE
3710
3711       DO 1011 I = 1, IHNT2(1)
3712          YR = SQRT(YP(1, I) ** 2 + YP(2, I) ** 2)
3713          IR = 1 + int(YR / DR)
3714 clin-4/2008 protect against out-of-bound errors:
3715 c         IF (IR .GT. 50) GOTO 601
3716          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601
3717          dnrpj1(IR) = dnrpj1(IR) + 1.0
3718          dnrtt1(IR) = dnrtt1(IR) + 1.0
3719  601     CONTINUE
3720  1011 CONTINUE
3721
3722       DO 1012 I = 1, IHNT2(3)
3723          YR = SQRT(YT(1, I) ** 2 + YT(2, I) ** 2)
3724          IR = 1 + int(YR / DR)
3725          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 602
3726          dnrtg1(IR) = dnrtg1(IR) + 1.0
3727          dnrtt1(IR) = dnrtt1(IR) + 1.0
3728  602     CONTINUE
3729  1012 CONTINUE
3730
3731       DO 1013 I = 1, NSG
3732          Y1 = 0.5 * (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))
3733          Y2 = 0.5 * (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))
3734          YR = SQRT(Y1 ** 2 + Y2 ** 2)
3735          IR = 1 + int(YR / DR)
3736          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603
3737          dnrin1(IR) = dnrin1(IR) + 1.0
3738          dnrtt1(IR) = dnrtt1(IR) + 1.0
3739  603     CONTINUE
3740  1013 CONTINUE
3741
3742       DO 1014 I = 1, MUL
3743          ITYP = ITYP0(I)
3744          PX = sngl(PX0(I))
3745          PY = sngl(PY0(I))
3746          PZ = sngl(PZ0(I))
3747          PE = sngl(E0(I))
3748          PM = sngl(XMASS0(I))
3749          XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
3750          DXMT = XMT - PM
3751 clin-9/2012 determine rapidity more generally:
3752 c         IF (ABS(PZ) .GE. PE) THEN
3753 c            PRINT *, ' IN HJANA1, GLUON ', I
3754 c            PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
3755 c            PRINT *, ' PZ = ', PZ, ' EE = ', PE
3756 c            PRINT *, ' XM = ', PM
3757 c            GOTO 800
3758 c         END IF
3759 c         RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
3760          if(XMT.gt.0.) then
3761             RAP=asinh(PZ/XMT)
3762          else
3763             PRINT *, ' IN HJANA1 mt=0'
3764             RAP = 1000000.0*sign(1.,PZ)
3765          endif
3766
3767          IY = 1 + int(ABS(RAP) / DY)
3768 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
3769 c         IF (IY .GT. 50) GOTO 700
3770          IF (IY.lt.1 .or.IY .GT. 50) GOTO 700
3771          dyg1c(IY) = dyg1c(IY) + 1.0
3772          deyg1c(IY) = deyg1c(IY) + XMT
3773  700     CONTINUE
3774          IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800
3775          IMT = 1 + int(DXMT / DMT)
3776          IF (IMT .GT. 50) GOTO 800
3777          dmyg1c(IMT) = dmyg1c(IMT) + 1.0 / XMT
3778  800     CONTINUE
3779  1014 CONTINUE
3780 c.....count number of particles
3781       DO 1016 I = 1, IHNT2(1)
3782          DO 1015 J = 1, NPJ(I)
3783             nsubp = nsubp + 1
3784             IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
3785  1015    CONTINUE
3786  1016 CONTINUE
3787
3788       DO 1018 I = 1, IHNT2(3)
3789          DO 1017 J = 1, NTJ(I)
3790             nsubp = nsubp + 1
3791             IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
3792  1017    CONTINUE
3793  1018 CONTINUE
3794
3795       DO 1020 I = 1, NSG
3796          DO 1019 J = 1, NJSG(I)
3797             nsubp = nsubp + 1
3798             IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
3799  1019    CONTINUE
3800  1020 CONTINUE
3801       NISG = NISG + NSG
3802       IF (IOUT .EQ. 1) THEN
3803 cbzdbg2/16/99
3804 c      PRINT *, ' in HJANA1 '
3805 c      PRINT *, ' total number of partons = ', nsubp
3806 c      PRINT *, ' total number of gluons = ', nsubg, MUL
3807 c      PRINT *, ' number of projectile strings = ', IHNT2(1)
3808 c      PRINT *, ' number of target strings = ', IHNT2(3)
3809 c      PRINT *, ' number of independent strings = ', NSG
3810       PRINT *, ' in HJANA1 '
3811       PRINT *, ' total number of partons = ', nsubp / IW
3812       PRINT *, ' total number of gluons = ', nsubg / IW
3813 c      PRINT *, ' number of projectile strings = ', IHNT2(1)
3814 c      PRINT *, ' number of target strings = ', IHNT2(3)
3815       PRINT *, ' number of independent strings = ', NISG / IW
3816 cbzdbg2/16/99end
3817       END IF
3818 c
3819       RETURN
3820       END
3821
3822 c-----------------------------------------------------------------------
3823
3824 c.....analysis subroutine in ZPC after generation of additional initial
3825 c.....phase space distributions.
3826
3827       SUBROUTINE HJAN1A
3828       PARAMETER (MAXPTN=400001)
3829       PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2)
3830       DIMENSION dgxg1a(50), dgyg1a(50), dtg1a(50)
3831       DIMENSION sgxg1a(50), sgyg1a(50), stg1a(50)
3832       DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
3833       COMMON /PARA1/ MUL
3834 cc      SAVE /PARA1/
3835       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
3836      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
3837      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
3838 cc      SAVE /prec2/
3839       COMMON /AREVT/ IAEVT, IARUN, MISS
3840 cc      SAVE /AREVT/
3841       COMMON /AROUT/ IOUT
3842 cc      SAVE /AROUT/
3843       SAVE
3844       DATA IW/0/
3845
3846       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3847          DO 1001 I = 1, 50
3848             dgxg1a(I) = sgxg1a(I)
3849             dgyg1a(I) = sgyg1a(I)
3850             dtg1a(I) = stg1a(I)
3851  1001    CONTINUE
3852       ELSE
3853          DO 1002 I = 1, 50
3854             sgxg1a(I) = dgxg1a(I)
3855             sgyg1a(I) = dgyg1a(I)
3856             stg1a(I) = dtg1a(I)
3857  1002    CONTINUE
3858          isevt = IAEVT
3859          isrun = IARUN
3860          IW = IW + 1
3861       END IF
3862 c.....analysis
3863       DO 1003 I = 1, MUL
3864          IGX = 1 + int(sngl(ABS(GX5(I))) / DGX)
3865 clin-4/2008 protect against out-of-bound errors:
3866 c         IF (IGX .GT. 50) GOTO 100
3867          IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100
3868          dgxg1a(IGX) = dgxg1a(IGX) + 1.0
3869  100     CONTINUE
3870          IGY = 1 + int(sngl(ABS(GY5(I))) / DGY)
3871          IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200
3872          dgyg1a(IGY) = dgyg1a(IGY) + 1.0
3873  200     CONTINUE
3874          IT = 1 + int(sngl(SQRT(FT5(I) ** 2 - GZ5(I) ** 2)) / DT)
3875          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 300
3876          dtg1a(IT) = dtg1a(IT) + 1.0
3877  300     CONTINUE
3878  1003 CONTINUE
3879       CALL HJAN1B
3880 c
3881       RETURN
3882       END
3883
3884 c-----------------------------------------------------------------------
3885
3886 c.....analysis subroutine in HJAN1A
3887
3888       SUBROUTINE HJAN1B
3889       PARAMETER (MAXPTN=400001,MAXSTR=150001)
3890       PARAMETER (DR = 0.2, DT = 0.2)
3891       DIMENSION DNRG1B(50), dtg1b(50)
3892       DIMENSION SNRG1B(50), stg1b(50)
3893       DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
3894       COMMON /PARA1/ MUL
3895 cc      SAVE /PARA1/
3896       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
3897      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
3898      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
3899 cc      SAVE /prec2/
3900       COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
3901 cc      SAVE /ilist8/
3902       COMMON /SREC1/ NSP, NST, NSI
3903 cc      SAVE /SREC1/
3904       COMMON/hjcrdn/YP(3,300),YT(3,300)
3905 cc      SAVE /hjcrdn/
3906       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3907      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3908      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3909 cc      SAVE /HJJET2/
3910       COMMON /AREVT/ IAEVT, IARUN, MISS
3911 cc      SAVE /AREVT/
3912       COMMON /AROUT/ IOUT
3913 cc      SAVE /AROUT/
3914       SAVE
3915       DATA IW/0/
3916
3917       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
3918          DO 1001 I = 1, 50
3919             DNRG1B(I) = SNRG1B(I)
3920             dtg1b(I) = stg1b(I)
3921  1001    CONTINUE
3922       ELSE
3923          DO 1002 I = 1, 50
3924             SNRG1B(I) = DNRG1B(I)
3925             stg1b(I) = dtg1b(I)
3926  1002    CONTINUE
3927          isevt = IAEVT
3928          isrun = IARUN
3929          IW = IW + 1
3930       END IF
3931 c.....analysis
3932       DO 1003 I = 1, MUL
3933          J = LSTRG1(I)
3934
3935          IF (J .LE. NSP) THEN
3936             K = J
3937             GX0 = YP(1, J)
3938             GY0 = YP(2, J)
3939          ELSE IF (J .LE. NSP + NST) THEN
3940             K = J - NSP
3941             GX0 = YT(1, K)
3942             GY0 = YT(2, K)
3943          ELSE
3944             K = J - NSP - NST
3945             GX0 = 0.5 * (YP(1, IASG(K, 1)) + YT(1, IASG(K, 2)))
3946             GY0 = 0.5 * (YP(2, IASG(K, 1)) + YT(2, IASG(K, 2)))
3947          END IF
3948          R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2)
3949          IR = 1 + int(R0 / DR)
3950          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100
3951          DNRG1B(IR) = DNRG1B(IR) + 1.0
3952  100     CONTINUE
3953          TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2))
3954          IT = 1 + int(TAU7 / DT)
3955          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 200
3956          dtg1b(IT) = dtg1b(IT) + 1.0
3957  200     CONTINUE
3958  1003 CONTINUE
3959 c
3960       RETURN
3961       END
3962
3963 c-----------------------------------------------------------------------
3964
3965 c.....analysis subroutine in HIJING after parton cascade evolution
3966       SUBROUTINE HJANA2
3967 c
3968       PARAMETER (YMAX = 1.0, YMIN = -1.0)
3969       PARAMETER (DMT = 0.05, DY = 0.2)
3970       PARAMETER (DR = 0.2, DT = 0.2)
3971       PARAMETER (MAXPTN=400001)
3972       PARAMETER (MAXSTR=150001)
3973       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
3974      1     GXSGS,GYSGS,GZSGS,FTSGS
3975       DIMENSION dyp2(50), DMYP2(200), DEYP2(50)
3976       DIMENSION dyg2(50), DMYG2(200), DEYG2(50)
3977       DIMENSION SNYP2(50), SMYP2(200), SEYP2(50)
3978       DIMENSION SNYG2(50), SMYG2(200), SEYG2(50)
3979       DIMENSION dnrpj2(50), dnrtg2(50), dnrin2(50),
3980      &   dnrtt2(50)
3981       DIMENSION dtpj2(50), dttg2(50), dtin2(50),
3982      &   dttot2(50)
3983       DIMENSION dyg2c(50), dmyg2c(50), deyg2c(50)
3984       DIMENSION snrpj2(50), snrtg2(50), snrin2(50),
3985      &   snrtt2(50)
3986       DIMENSION stpj2(50), sttg2(50), stin2(50),
3987      &   sttot2(50)
3988       DIMENSION snyg2c(50), smyg2c(50), seyg2c(50)
3989       DOUBLE PRECISION  ATAUI, ZT1, ZT2, ZT3
3990       DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
3991       COMMON /PARA1/ MUL
3992 cc      SAVE /PARA1/
3993       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3994 cc      SAVE /HPARNT/
3995       COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
3996 cc      SAVE /SREC2/
3997       COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3998      &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3999      &   PJPM(300,500),NTJ(300),KFTJ(300,500),
4000      &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
4001      &   PJTE(300,500),PJTM(300,500)
4002 cc      SAVE /HJJET1/
4003       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
4004      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
4005      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
4006 cc      SAVE /HJJET2/
4007       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
4008      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
4009      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
4010 cc      SAVE /prec2/
4011       COMMON /AREVT/ IAEVT, IARUN, MISS
4012 cc      SAVE /AREVT/
4013       COMMON /AROUT/ IOUT
4014 cc      SAVE /AROUT/
4015       common/anim/nevent,isoft,isflag,izpc
4016 cc      SAVE /anim/
4017       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
4018      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
4019      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
4020      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
4021 cc      SAVE /SOFT/
4022       SAVE
4023       DATA IW/0/
4024
4025       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
4026          DO 1001 I = 1, 200
4027             DMYP2(I) = SMYP2(I)
4028             DMYG2(I) = SMYG2(I)
4029  1001    CONTINUE
4030
4031          DO 1002 I = 1, 50
4032             dyp2(I) = SNYP2(I)
4033             DEYP2(I) = SEYP2(I)
4034             dyg2(I) = SNYG2(I)
4035             DEYG2(I) = SEYG2(I)
4036             dnrpj2(I) = snrpj2(I)
4037             dnrtg2(I) = snrtg2(I)
4038             dnrin2(I) = snrin2(I)
4039             dnrtt2(I) = snrtt2(I)
4040             dtpj2(I) = stpj2(I)
4041             dttg2(I) = sttg2(I)
4042             dtin2(I) = stin2(I)
4043             dttot2(I) = sttot2(I)
4044             dyg2c(I) = snyg2c(I)
4045             dmyg2c(I) = smyg2c(I)
4046             deyg2c(I) = seyg2c(I)
4047  1002    CONTINUE
4048          nsubp = nsubpS
4049          nsubg = nsubgS
4050          NISG = NISGS
4051       ELSE
4052          DO 1003 I = 1, 200
4053             SMYP2(I) = DMYP2(I)
4054             SMYG2(I) = DMYG2(I)
4055  1003    CONTINUE
4056
4057          DO 1004 I = 1, 50
4058             SNYP2(I) = dyp2(I)
4059             SEYP2(I) = DEYP2(I)
4060             SNYG2(I) = dyg2(I)
4061             SEYG2(I) = DEYG2(I)
4062             snrpj2(I) = dnrpj2(I)
4063             snrtg2(I) = dnrtg2(I)
4064             snrin2(I) = dnrin2(I)
4065             snrtt2(I) = dnrtt2(I)
4066             stpj2(I) = dtpj2(I)
4067             sttg2(I) = dttg2(I)
4068             stin2(I) = dtin2(I)
4069             sttot2(I) = dttot2(I)
4070             snyg2c(I) = dyg2c(I)
4071             smyg2c(I) = dmyg2c(I)
4072             seyg2c(I) = deyg2c(I)
4073  1004    CONTINUE
4074          nsubpS = nsubp
4075          nsubgS = nsubg
4076          NISGS = NISG
4077          isevt = IAEVT
4078          isrun = IARUN
4079          IW = IW + 1
4080       END IF
4081
4082 clin-4/28/01:
4083       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 510
4084
4085 c.....analysis
4086       DO 1006 I = 1, IHNT2(1)
4087          DO 1005 J = 1, NPJ(I)
4088             ITYP = KFPJ(I, J)
4089             PX = PJPX(I, J)
4090             PY = PJPY(I, J)
4091             PZ = PJPZ(I, J)
4092             PE = PJPE(I, J)
4093             PM = PJPM(I, J)
4094             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
4095             DXMT = XMT - PM
4096 clin-9/2012 determine rapidity more generally:
4097 ccbzdbg2/16/99
4098 cc            IF (ABS(PZ) .GE. PE) GOTO 200
4099 c            IF (ABS(PZ) .GE. PE) THEN
4100 c               PRINT *, ' IN HJANA2, PROJ STR ', I, ' PART ', J
4101 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4102 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
4103 c               PRINT *, ' XM = ', PM
4104 c               GOTO 200
4105 c            END IF
4106 ccbzdbg2/16/99end
4107 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
4108             if(XMT.gt.0.) then
4109                RAP=asinh(PZ/XMT)
4110             else
4111                PRINT *, ' IN HJANA2 mt=0'
4112                RAP = 1000000.0*sign(1.,PZ)
4113             endif
4114
4115             IY = 1 + int(ABS(RAP) / DY)
4116 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
4117 c            IF (IY .GT. 50) GOTO 100
4118             IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
4119             dyp2(IY) = dyp2(IY) + 1.0
4120             DEYP2(IY) = DEYP2(IY) + XMT
4121             IF (ITYP .EQ. 21) THEN
4122                dyg2(IY) = dyg2(IY) + 1.0
4123                DEYG2(IY) = DEYG2(IY) + XMT
4124             END IF
4125  100        CONTINUE
4126             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 200
4127             IMT = 1 + int(DXMT / DMT)
4128             IF (IMT .GT. 200) GOTO 200
4129             DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
4130             IF (ITYP .EQ. 21) THEN
4131                DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
4132             END IF
4133  200        CONTINUE
4134  1005    CONTINUE
4135  1006 CONTINUE
4136
4137       DO 1008 I = 1, IHNT2(3)
4138          DO 1007 J = 1, NTJ(I)
4139             ITYP = KFTJ(I, J)
4140             PX = PJTX(I, J)
4141             PY = PJTY(I, J)
4142             PZ = PJTZ(I, J)
4143             PE = PJTE(I, J)
4144             PM = PJTM(I, J)
4145             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
4146             DXMT = XMT - PM
4147 clin-9/2012 determine rapidity more generally:
4148 ccbzdbg2/16/99
4149 cc            IF (ABS(PZ) .GE. PE) GOTO 400
4150 c            IF (ABS(PZ) .GE. PE) THEN
4151 c               PRINT *, ' IN HJANA2, TARG STR ', I, ' PART ', J
4152 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4153 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
4154 c               PRINT *, ' XM = ', PM
4155 c               GOTO 400
4156 c            END IF
4157 ccbzdbg2/16/99end
4158 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
4159             if(XMT.gt.0.) then
4160                RAP=asinh(PZ/XMT)
4161             else
4162                PRINT *, ' IN HJANA2 mt=0'
4163                RAP = 1000000.0*sign(1.,PZ)
4164             endif
4165
4166             IY = 1 + int(ABS(RAP) / DY)
4167 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
4168 c            IF (IY .GT. 50) GOTO 300
4169             IF (IY.lt.1 .or.IY .GT. 50) GOTO 300
4170             dyp2(IY) = dyp2(IY) + 1.0
4171             DEYP2(IY) = DEYP2(IY) + XMT
4172             IF (ITYP .EQ. 21) THEN
4173                dyg2(IY) = dyg2(IY) + 1.0
4174                DEYG2(IY) = DEYG2(IY) + XMT
4175             END IF
4176  300        CONTINUE
4177             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 400
4178             IMT = 1 + int(DXMT / DMT)
4179             IF (IMT .GT. 200) GOTO 400
4180             DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
4181             IF (ITYP .EQ. 21) THEN
4182                DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
4183             END IF
4184  400        CONTINUE
4185  1007    CONTINUE
4186  1008 CONTINUE
4187
4188 clin-4/28/01:
4189  510  continue
4190
4191       DO 1010 I = 1, NSG
4192 clin-4/25/01 soft3:
4193 c         DO J = 1, NJSG(I)
4194          NJ=NJSG(I)
4195          if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
4196          DO 1009 J = 1, NJ
4197 clin-4/25/01-end
4198
4199             ITYP = K2SG(I, J)
4200             PX = PXSG(I, J)
4201             PY = PYSG(I, J)
4202             PZ = PZSG(I, J)
4203             PE = PESG(I, J)
4204             PM = PMSG(I, J)
4205 clin-4/25/01 soft3:
4206             if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
4207                ITYP = K2SGS(I, J)
4208                PX = sngl(PXSGS(I, J))
4209                PY = sngl(PYSGS(I, J))
4210                PZ = sngl(PZSGS(I, J))
4211                PE = sngl(PESGS(I, J))
4212                PM = sngl(PMSGS(I, J))
4213             endif
4214 clin-4/25/01-end
4215
4216 clin-9/2012 determine rapidity more generally:
4217             XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
4218             DXMT = XMT - PM
4219 ccbzdbg2/16/99
4220 cc            IF (ABS(PZ) .GE. PE) GOTO 600
4221 c            IF (ABS(PZ) .GE. PE) THEN
4222 c               PRINT *, ' IN HJANA2, INDP STR ', I, ' PART ', J
4223 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4224 c               PRINT *, ' PZ = ', PZ, ' EE = ', PE
4225 c               PRINT *, ' XM = ', PM
4226 c               GOTO 600
4227 c            END IF
4228 ccbzdbg2/16/99end
4229 c            RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
4230             if(XMT.gt.0.) then
4231                RAP=asinh(PZ/XMT)
4232             else
4233                PRINT *, ' IN HJANA2 mt=0'
4234                RAP = 1000000.0*sign(1.,PZ)
4235             endif
4236
4237             IY = 1 + int(ABS(RAP) / DY)
4238 clin-8/2014 prevent possible segmentation fault (due to IY<=0):
4239 c            IF (IY .GT. 50) GOTO 500
4240             IF (IY.lt.1 .or.IY .GT. 50) GOTO 500
4241             dyp2(IY) = dyp2(IY) + 1.0
4242             DEYP2(IY) = DEYP2(IY) + XMT
4243             IF (ITYP .EQ. 21) THEN
4244                dyg2(IY) = dyg2(IY) + 1.0
4245                DEYG2(IY) = DEYG2(IY) + XMT
4246             END IF
4247  500        CONTINUE
4248             IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 600
4249             IMT = 1 + int(DXMT / DMT)
4250             IF (IMT .GT. 200) GOTO 600
4251             DMYP2(IMT) = DMYP2(IMT) + 1.0 / XMT
4252             IF (ITYP .EQ. 21) THEN
4253                DMYG2(IMT) = DMYG2(IMT) + 1.0 / XMT
4254             END IF
4255  600        CONTINUE
4256  1009    CONTINUE
4257  1010 CONTINUE
4258
4259 clin-4/28/01:
4260       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 520
4261
4262       DO 1011 I = 1, IHNT2(1)
4263          J = I
4264          YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
4265          IR = 1 + int(YR / DR)
4266          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 601
4267          dnrpj2(IR) = dnrpj2(IR) + 1.0
4268          dnrtt2(IR) = dnrtt2(IR) + 1.0
4269  601     CONTINUE
4270          IT = 1 + int(sngl(ATAUI(J)) / DT)
4271          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 602
4272          dtpj2(IT) = dtpj2(IT) + 1.0
4273          dttot2(IT) = dttot2(IT) + 1.0
4274  602     CONTINUE
4275  1011 CONTINUE
4276
4277       DO 1012 I = 1, IHNT2(3)
4278          J = I + IHNT2(1)
4279          YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
4280          IR = 1 + int(YR / DR)
4281          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 603
4282          dnrtg2(IR) = dnrtg2(IR) + 1.0
4283          dnrtt2(IR) = dnrtt2(IR) + 1.0
4284  603     CONTINUE
4285          IT = 1 + int(sngl(ATAUI(J)) / DT)
4286          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 604
4287          dttg2(IT) = dttg2(IT) + 1.0
4288          dttot2(IT) = dttot2(IT) + 1.0
4289  604     CONTINUE
4290  1012 CONTINUE
4291
4292 clin-4/28/01:
4293  520  continue
4294
4295       DO 1013 I = 1, NSG
4296          J = I + IHNT2(1) + IHNT2(3)
4297 clin-4/28/01:
4298          if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) J = I
4299
4300          YR = SQRT(sngl(ZT1(J) ** 2 + ZT2(J) ** 2))
4301          IR = 1 + int(YR / DR)
4302          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 605
4303          dnrin2(IR) = dnrin2(IR) + 1.0
4304          dnrtt2(IR) = dnrtt2(IR) + 1.0
4305  605     CONTINUE
4306          IT = 1 + int(sngl(ATAUI(J)) / DT)
4307          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 606
4308          dtin2(IT) = dtin2(IT) + 1.0
4309          dttot2(IT) = dttot2(IT) + 1.0
4310  606     CONTINUE
4311  1013 CONTINUE
4312
4313       DO 1014 I = 1, MUL
4314          ITYP = ITYP5(I)
4315          PX = sngl(PX5(I))
4316          PY = sngl(PY5(I))
4317          PZ = sngl(PZ5(I))
4318          PE = sngl(E5(I))
4319          PM = sngl(XMASS5(I))
4320 clin-9/2012 determine rapidity more generally:
4321          XMT = SQRT(PX ** 2 + PY ** 2 + PM ** 2)
4322          DXMT = XMT - PM
4323 ccbzdbg2/16/99
4324 cc            IF (ABS(PZ) .GE. PE) GOTO 800
4325 c         IF (ABS(PZ) .GE. PE) THEN
4326 c            PRINT *, ' IN HJANA2, GLUON ', I
4327 c            PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4328 c            PRINT *, ' PZ = ', PZ, ' EE = ', PE
4329 c            PRINT *, ' XM = ', PM
4330 c            GOTO 800
4331 c         END IF
4332 ccbzdbg2/16/99end
4333 c         RAP = 0.5 * LOG((PE + PZ +1e-5) / (PE - PZ + 1e-5))
4334          if(XMT.gt.0.) then
4335             RAP=asinh(PZ/XMT)
4336          else
4337             PRINT *, ' IN HJANA2 mt=0'
4338             RAP = 1000000.0*sign(1.,PZ)
4339          endif
4340
4341          IY = 1 + int(ABS(RAP) / DY)
4342 clin-9/2012 prevent possible segmentation fault (due to IY<=0):
4343 c         IF (IY .GT. 50) GOTO 700
4344          IF (IY.lt.1 .or.IY .GT. 50) GOTO 700
4345          dyg2c(IY) = dyg2c(IY) + 1.0
4346          deyg2c(IY) = deyg2c(IY) + XMT
4347  700     CONTINUE
4348          IF (RAP .GT. YMAX .OR. RAP .LE. YMIN) GOTO 800
4349          IMT = 1 + int(DXMT / DMT)
4350          IF (IMT .GT. 50) GOTO 800
4351          dmyg2c(IMT) = dmyg2c(IMT) + 1.0 / XMT
4352  800     CONTINUE
4353  1014 CONTINUE
4354
4355 clin-4/25/01 soft3:
4356       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) goto 530
4357
4358 c.....count number of particles
4359       DO 1016 I = 1, IHNT2(1)
4360          DO 1015 J = 1, NPJ(I)
4361             nsubp = nsubp + 1
4362             IF (KFPJ(I, J) .EQ. 21) nsubg = nsubg + 1
4363  1015    CONTINUE
4364  1016 CONTINUE
4365
4366       DO 1018 I = 1, IHNT2(3)
4367          DO 1017 J = 1, NTJ(I)
4368             nsubp = nsubp + 1
4369             IF (KFTJ(I, J) .EQ. 21) nsubg = nsubg + 1
4370  1017    CONTINUE
4371  1018 CONTINUE
4372
4373 clin-4/25/01 soft3:
4374  530  continue
4375
4376       DO 1020 I = 1, NSG
4377 clin-4/25/01 soft3:
4378 c         DO J = 1, NJSG(I)
4379          NJ=NJSG(I)
4380          if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) NJ=NJSGS(I)
4381          DO 1019 J = 1, NJ
4382 clin-4/25/01-end
4383
4384             nsubp = nsubp + 1
4385
4386 clin-4/25/01
4387 c            IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
4388             if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
4389                IF(K2SGS(I, J) .EQ. 21) nsubg = nsubg + 1
4390             else
4391                IF (K2SG(I, J) .EQ. 21) nsubg = nsubg + 1
4392             endif
4393 clin-4/25/01-end
4394  1019    CONTINUE
4395  1020 CONTINUE
4396 cbzdbg2/16/99
4397       NISG = NISG + NSG
4398
4399       IF (IOUT .EQ. 1) THEN
4400 cbzdbg2/16/99end
4401 cbzdbg2/16/99
4402 c      PRINT *, ' in HJANA2 '
4403 c      PRINT *, ' total number of partons = ', nsubp
4404 c      PRINT *, ' total number of gluons = ', nsubg, MUL
4405 c      PRINT *, ' number of projectile strings = ', IHNT2(1)
4406 c      PRINT *, ' number of target strings = ', IHNT2(3)
4407 c      PRINT *, ' number of independent strings = ', NSG
4408       PRINT *, ' in HJANA2 '
4409       PRINT *, ' total number of partons = ', nsubp / IW
4410       PRINT *, ' total number of gluons = ', nsubg / IW
4411 c      PRINT *, ' number of projectile strings = ', IHNT2(1)
4412 c      PRINT *, ' number of target strings = ', IHNT2(3)
4413       PRINT *, ' number of independent strings = ', NISG / IW
4414       END IF
4415
4416       CALL HJAN2A
4417       CALL HJAN2B
4418
4419       RETURN
4420       END
4421
4422 c-----------------------------------------------------------------------
4423
4424 c.....subroutine called by HJANA2
4425       SUBROUTINE HJAN2A
4426
4427       PARAMETER (DGX = 0.2, DGY = 0.2, DT = 0.2)
4428       PARAMETER (MAXPTN=400001,MAXSTR=150001)
4429       DIMENSION dgxp2a(50), dgyp2a(50), dtp2a(50)
4430       DIMENSION dgxg2a(50), dgyg2a(50), dtg2a(50)
4431       DIMENSION sgxp2a(50), sgyp2a(50), stp2a(50)
4432       DIMENSION sgxg2a(50), sgyg2a(50), stg2a(50)
4433       DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
4434       COMMON /PARA1/ MUL
4435 cc      SAVE /PARA1/
4436       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
4437      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
4438      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
4439 cc      SAVE /prec2/
4440       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4441 cc      SAVE /HPARNT/
4442       COMMON/hjcrdn/YP(3,300),YT(3,300)
4443 cc      SAVE /hjcrdn/
4444       COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
4445      &   PJPY(300,500),PJPZ(300,500),PJPE(300,500),
4446      &   PJPM(300,500),NTJ(300),KFTJ(300,500),
4447      &   PJTX(300,500),PJTY(300,500),PJTZ(300,500),
4448      &   PJTE(300,500),PJTM(300,500)
4449 cc      SAVE /HJJET1/
4450       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
4451      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
4452      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
4453 cc      SAVE /HJJET2/
4454       COMMON /AREVT/ IAEVT, IARUN, MISS
4455 cc      SAVE /AREVT/
4456       COMMON /AROUT/ IOUT
4457 cc      SAVE /AROUT/
4458       SAVE
4459       DATA IW/0/
4460
4461       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
4462          DO 1001 I = 1, 50
4463             dgxp2a(I) = sgxp2a(I)
4464             dgyp2a(I) = sgyp2a(I)
4465             dtp2a(I) = stp2a(I)
4466             dgxg2a(I) = sgxg2a(I)
4467             dgyg2a(I) = sgyg2a(I)
4468             dtg2a(I) = stg2a(I)
4469  1001    CONTINUE
4470       ELSE
4471          DO 1002 I = 1, 50
4472             sgxp2a(I) = dgxp2a(I)
4473             sgyp2a(I) = dgyp2a(I)
4474             stp2a(I) = dtp2a(I)
4475             sgxg2a(I) = dgxg2a(I)
4476             sgyg2a(I) = dgyg2a(I)
4477             stg2a(I) = dtg2a(I)
4478  1002    CONTINUE
4479          isevt = IAEVT
4480          isrun = IARUN
4481          IW = IW + 1
4482       END IF
4483 c.....analysis
4484       DO 1004 I = 1, IHNT2(1)
4485          DO 1003 J = 1, NPJ(I)
4486             IF (KFPJ(I, J) .NE. 21) THEN
4487                IGX = 1 + int(ABS(YP(1, I)) / DGX)
4488                IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 100
4489                dgxp2a(IGX) = dgxp2a(IGX) + 1.0
4490  100           CONTINUE
4491                IGY = 1 + int(ABS(YP(2, I)) / DGY)
4492                IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 200
4493                dgyp2a(IGY) = dgyp2a(IGY) + 1.0
4494  200           CONTINUE
4495                IT = 1
4496                dtp2a(IT) = dtp2a(IT) + 1.0
4497             END IF
4498  1003    CONTINUE
4499  1004 CONTINUE
4500
4501       DO 1006 I = 1, IHNT2(3)
4502          DO 1005 J = 1, NTJ(I)
4503             IF (KFTJ(I, J) .NE. 21) THEN
4504                IGX = 1 + int(ABS(YT(1, I)) / DGX)
4505                IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 300
4506                dgxp2a(IGX) = dgxp2a(IGX) + 1.0
4507  300           CONTINUE
4508                IGY = 1 + int(ABS(YT(2, I)) / DGY)
4509                IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 400
4510                dgyp2a(IGY) = dgyp2a(IGY) + 1.0
4511  400           CONTINUE
4512                IT = 1
4513                dtp2a(IT) = dtp2a(IT) + 1.0
4514             END IF
4515  1005    CONTINUE
4516  1006 CONTINUE
4517
4518       DO 1008 I = 1, NSG
4519          DO 1007 J = 1, NJSG(I)
4520             IF (K2SG(I, J) .NE. 21) THEN
4521                IGX = 1 + int(ABS(0.5 *
4522      &            (YP(1, IASG(I, 1)) + YT(1, IASG(I, 2)))) / DGX)
4523                IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 500
4524                dgxp2a(IGX) = dgxp2a(IGX) + 1.0
4525  500           CONTINUE
4526                IGY = 1 + int(ABS(0.5 *
4527      &            (YP(2, IASG(I, 1)) + YT(2, IASG(I, 2)))) / DGY)
4528                IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 600
4529                dgyp2a(IGY) = dgyp2a(IGY) + 1.0
4530  600           CONTINUE
4531                IT = 1
4532                dtp2a(IT) = dtp2a(IT) + 1.0
4533             END IF
4534  1007    CONTINUE
4535  1008 CONTINUE
4536
4537       DO 1009 I = 1, MUL
4538          IGX = 1 + int(ABS(sngl(GX5(I))) / DGX)
4539          IF (IGX .GT. 50 .or. IGX .LT. 1) GOTO 700
4540          dgxg2a(IGX) = dgxg2a(IGX) + 1.0
4541          dgxp2a(IGX) = dgxp2a(IGX) + 1.0
4542  700     CONTINUE
4543          IGY = 1 + int(ABS(sngl(GY5(I))) / DGY)
4544          IF (IGY .GT. 50 .or. IGY .LT. 1) GOTO 800
4545          dgyg2a(IGY) = dgyg2a(IGY) + 1.0
4546          dgyp2a(IGY) = dgyp2a(IGY) + 1.0
4547  800     CONTINUE
4548          IT = 1 + int(SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2)) / DT)
4549          IF (IT .GT. 50 .or. IT .LT. 1) GOTO 900
4550          dtg2a(IT) = dtg2a(IT) + 1.0
4551          dtp2a(IT) = dtp2a(IT) + 1.0
4552  900     CONTINUE
4553  1009 CONTINUE
4554 c
4555       RETURN
4556       END
4557
4558 c-----------------------------------------------------------------------
4559
4560 c.....analysis subroutine in HJANA2
4561
4562       SUBROUTINE HJAN2B
4563
4564       PARAMETER (MAXPTN=400001)
4565       PARAMETER (MAXSTR=150001)
4566       PARAMETER (DR = 0.2, DT = 0.2)
4567       DIMENSION DNRG2B(50), dtg2b(-24:25)
4568       DIMENSION SNRG2B(50), stg2b(-24:25)
4569       DOUBLE PRECISION  GX5, GY5, GZ5, FT5, PX5, PY5, PZ5, E5, XMASS5
4570       DOUBLE PRECISION  ATAUI, ZT1, ZT2, ZT3
4571       COMMON /PARA1/ MUL
4572 cc      SAVE /PARA1/
4573       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
4574      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
4575      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
4576 cc      SAVE /prec2/
4577       COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4578 cc      SAVE /ilist8/
4579       COMMON /SREC1/ NSP, NST, NSI
4580 cc      SAVE /SREC1/
4581       COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4582 cc      SAVE /SREC2/
4583       COMMON/hjcrdn/YP(3,300),YT(3,300)
4584 cc      SAVE /hjcrdn/
4585       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
4586      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
4587      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
4588 cc      SAVE /HJJET2/
4589       COMMON /AREVT/ IAEVT, IARUN, MISS
4590 cc      SAVE /AREVT/
4591       COMMON /AROUT/ IOUT
4592 cc      SAVE /AROUT/
4593       SAVE
4594       DATA IW/0/
4595
4596       IF (isevt .EQ. IAEVT .AND. isrun .EQ. IARUN) THEN
4597          DO 1001 I = 1, 50
4598             DNRG2B(I) = SNRG2B(I)
4599             dtg2b(I - 25) = stg2b(I - 25)
4600  1001    CONTINUE
4601       ELSE
4602          DO 1002 I = 1, 50
4603             SNRG2B(I) = DNRG2B(I)
4604             stg2b(I - 25) = dtg2b(I - 25)
4605  1002    CONTINUE
4606          isevt = IAEVT
4607          isrun = IARUN
4608          IW = IW + 1
4609       END IF
4610 c.....analysis
4611       DO 1003 I = 1, MUL
4612          J = LSTRG1(I)
4613          GX0 = sngl(ZT1(J))
4614          GY0 = sngl(ZT2(J))
4615          R0 = SQRT((sngl(GX5(I)) - GX0)**2 + (sngl(GY5(I)) - GY0)**2)
4616          IR = 1 + int(R0 / DR)
4617          IF (IR .GT. 50 .or. IR .LT. 1) GOTO 100
4618          DNRG2B(IR) = DNRG2B(IR) + 1.0
4619  100     CONTINUE
4620          TAU7 = SQRT(sngl(FT5(I) ** 2 - GZ5(I) ** 2))
4621          DTAU=TAU7 - sngl(ATAUI(J))
4622          IT = 1 + int(DTAU / DT)
4623 cbzdbg2/21/99
4624 c         IF (ABS(IT) .GT. 25) GOTO 200
4625          IF (IT .GT. 25 .OR. IT .LT. -24) GOTO 200
4626 cbzdbg2/21/99end
4627          dtg2b(IT) = dtg2b(IT) + 1.0
4628  200     CONTINUE
4629  1003 CONTINUE
4630 c
4631       RETURN
4632       END
4633
4634 c-----------------------------------------------------------------------
4635
4636 c.....analysis subroutine before ARTMN
4637       SUBROUTINE HJANA3
4638 c
4639       PARAMETER (MAXSTR=150001, MAXR=1)
4640 c.....y cut for mt spectrum
4641       PARAMETER (YMIN = -1.0, YMAX = 1.0)
4642 cbz11/7/99 end
4643 c.....bin width for mt spectrum and y spectrum
4644       PARAMETER (DMT = 0.05, DY = 0.2)
4645       DOUBLE PRECISION v2i,eti,xmulti,v2mi,s2mi,xmmult,
4646      1     v2bi,s2bi,xbmult
4647       DIMENSION dndyh3(50), DMYH3(50), DEYH3(50)
4648       COMMON /RUN/ NUM
4649 cc      SAVE /RUN/
4650       COMMON /ARERC1/MULTI1(MAXR)
4651 cc      SAVE /ARERC1/
4652       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
4653      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
4654      &     FT1(MAXSTR, MAXR),
4655      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
4656      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
4657 cc      SAVE /ARPRC1/
4658       COMMON /AROUT/ IOUT
4659 cc      SAVE /AROUT/
4660       COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
4661 cc      SAVE /iflow/
4662       SAVE
4663       DATA IW/0/
4664
4665       IW = IW + 1
4666       DO 1002 J = 1, NUM
4667          DO 1001 I = 1, MULTI1(J)
4668             ITYP = ITYP1(I, J)
4669             IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
4670             PX = PX1(I, J)
4671             PY = PY1(I, J)
4672             PZ = PZ1(I, J)
4673             EE = EE1(I, J)
4674             XM = XM1(I, J)
4675             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
4676             DXMT = XMT - XM
4677 clin-9/2012 determine rapidity more generally:
4678 c            IF (ABS(PZ) .GE. EE) THEN
4679 c               PRINT *, 'IN HJANA3'
4680 c               PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR'
4681 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4682 c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
4683 c               PRINT *, ' XM = ', XM
4684 c               GOTO 200
4685 c            END IF
4686 c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
4687             if(XMT.gt.0.) then
4688                Y=asinh(PZ/XMT)
4689             else
4690                PRINT *, ' IN HJANA3 mt=0'
4691                Y = 1000000.0*sign(1.,PZ)
4692             endif
4693
4694 c.....rapidity cut for the rapidity distribution
4695 c            IY = 1 + int(ABS(Y) / DY)
4696 clin-8/2014 no rapidity shift here:
4697 c            IY = 1 + int((Y+10.) / DY)
4698             IY = 1 + int(Y/DY)
4699 clin-9/2012 prevent possible segmentation fault (due to IY<=0):
4700 c            IF (IY .GT. 50) GOTO 100
4701             IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
4702             dndyh3(IY) = dndyh3(IY) + 1.0
4703             DEYH3(IY) = DEYH3(IY) + XMT
4704  100        CONTINUE
4705 c.....insert rapidity cut for mt spectrum here
4706             IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
4707             IMT = 1 + int(DXMT / DMT)
4708             IF (IMT .GT. 50) GOTO 200
4709             DMYH3(IMT) = DMYH3(IMT) + 1.0 / XMT
4710  200        CONTINUE
4711  1001    CONTINUE
4712  1002 CONTINUE
4713 c
4714       RETURN
4715       END
4716
4717 c-----------------------------------------------------------------------
4718
4719 c.....analysis subroutine after ARTMN
4720       SUBROUTINE HJANA4
4721       PARAMETER (MAXSTR=150001, MAXR=1)
4722 c.....y cut for mt spectrum
4723 cbz11/7/99
4724 c      PARAMETER (YMIN = -0.5, YMAX = 0.5)
4725       PARAMETER (YMIN = -1.0, YMAX = 1.0)
4726 cbz11/7/99 end
4727 c.....bin width for mt spectrum and y spectrum
4728       PARAMETER (DMT = 0.05, DY = 0.2)
4729
4730       DIMENSION dndyh4(50), DMYH4(50), DEYH4(50)
4731       COMMON /RUN/ NUM
4732 cc      SAVE /RUN/
4733       COMMON /ARERC1/MULTI1(MAXR)
4734 cc      SAVE /ARERC1/
4735       COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
4736      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
4737      &     FT1(MAXSTR, MAXR),
4738      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
4739      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
4740 cc      SAVE /ARPRC1/
4741       COMMON /AROUT/ IOUT
4742 cc      SAVE /AROUT/
4743       COMMON /fflow/ v2f,etf,xmultf,v2fpi,xmulpi
4744 cc      SAVE /fflow/
4745       SAVE
4746       DATA IW/0/
4747
4748       IW = IW + 1
4749       DO 1002 J = 1, NUM
4750          DO 1001 I = 1, MULTI1(J)
4751             ITYP = ITYP1(I, J)
4752             IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 200
4753             PX = PX1(I, J)
4754             PY = PY1(I, J)
4755             PZ = PZ1(I, J)
4756             EE = EE1(I, J)
4757             XM = XM1(I, J)
4758             XMT = SQRT(PX ** 2 + PY ** 2 + XM ** 2)
4759             DXMT = XMT - XM
4760 clin-9/2012 determine rapidity more generally:
4761 c            IF (ABS(PZ) .GE. EE) THEN
4762 c               PRINT *, 'IN HJANA4'
4763 c               PRINT *, ' PARTICLE ', I, ' RUN ', J, 'PREC ERR'
4764 c               PRINT *, ' FLAV = ', ITYP, ' PX = ', PX, ' PY = ', PY
4765 c               PRINT *, ' PZ = ', PZ, ' EE = ', EE
4766 c               PRINT *, ' XM = ', XM
4767 c               GOTO 200
4768 c            END IF
4769 c            Y = 0.5 * LOG((EE + PZ +1e-5) / (EE - PZ + 1e-5))
4770             if(XMT.gt.0.) then
4771                Y=asinh(PZ/XMT)
4772             else
4773                PRINT *, ' IN HJANA4 mt=0'
4774                Y = 1000000.0*sign(1.,PZ)
4775             endif
4776
4777 c.....rapidity cut for the rapidity distribution
4778 c            IY = 1 + int(ABS(Y) / DY)
4779 clin-8/2014 no rapidity shift here:
4780 c            IY = 1 + int((Y+10.) / DY)
4781             IY = 1 + int(Y/DY)
4782 clin-9/2012 prevent possible segmentation fault (due to IY<=0):
4783 c            IF (IY .GT. 50) GOTO 100
4784             IF (IY.lt.1 .or.IY .GT. 50) GOTO 100
4785             dndyh4(IY) = dndyh4(IY) + 1.0
4786             DEYH4(IY) = DEYH4(IY) + XMT
4787  100        CONTINUE
4788 c.....insert rapidity cut for mt spectrum here
4789             IF (Y. LT. YMIN .OR. Y .GE. YMAX) GOTO 200
4790             IMT = 1 + int(DXMT / DMT)
4791             IF (IMT .GT. 50) GOTO 200
4792             DMYH4(IMT) = DMYH4(IMT) + 1.0 / XMT
4793  200        CONTINUE
4794  1001    CONTINUE
4795  1002 CONTINUE
4796 c
4797       RETURN
4798       END
4799
4800 c=======================================================================
4801
4802 c.....subroutine to get average values for different strings
4803
4804       SUBROUTINE zpstrg
4805
4806       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
4807       PARAMETER (MAXPTN=400001)
4808       PARAMETER (MAXSTR=150001)
4809 c      REAL*4 YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB
4810       REAL YP, YT, PXSG, PYSG, PZSG, PESG, PMSG, HIPR1, HINT1, BB
4811
4812       COMMON /PARA1/ MUL
4813 cc      SAVE /PARA1/
4814       COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
4815      &   PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
4816      &   XMASS5(MAXPTN), ITYP5(MAXPTN)
4817 cc      SAVE /prec2/
4818       COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
4819 cc      SAVE /ilist8/
4820       COMMON /SREC1/ NSP, NST, NSI
4821 cc      SAVE /SREC1/
4822       COMMON /SREC2/ATAUI(MAXSTR),ZT1(MAXSTR),ZT2(MAXSTR),ZT3(MAXSTR)
4823 cc      SAVE /SREC2/
4824       COMMON/hjcrdn/YP(3,300),YT(3,300)
4825 cc      SAVE /hjcrdn/
4826       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
4827      &   K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
4828      &   PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
4829 cc      SAVE /HJJET2/
4830 cbz6/28/99 flow1
4831       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
4832 cc      SAVE /HPARNT/
4833 cbz6/28/99 flow1 end
4834       common/anim/nevent,isoft,isflag,izpc
4835 cc      SAVE /anim/
4836       common/strg/np(maxstr)
4837 cc      SAVE /strg/
4838 clin-6/06/02 test local freezeout:
4839       common /frzprc/
4840      &     gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
4841      &     pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
4842      &     xmfrz(MAXPTN),
4843      &     tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
4844 cc      SAVE /frzprc/
4845       SAVE
4846
4847 clin-6/06/02 test local freezeout for string melting,
4848 c     use space-time values at local freezeout saved in /frzprc/:
4849       if(isoft.eq.5) then
4850          do 1001 I = 1, MUL
4851             ITYP5(i)=idfrz(i)
4852             GX5(i)=gxfrz(i)
4853             GY5(i)=gyfrz(i)
4854             GZ5(i)=gzfrz(i)
4855             FT5(i)=ftfrz(i)
4856             PX5(i)=pxfrz(i)
4857             PY5(i)=pyfrz(i)
4858             PZ5(i)=pzfrz(i)
4859             E5(i)=efrz(i)
4860             XMASS5(i)=xmfrz(i)
4861  1001    continue
4862       endif
4863 clin-6/06/02-end
4864
4865       DO 1002 I = 1, MAXSTR
4866          ATAUI(I) = 0d0
4867          ZT1(I) = 0d0
4868          ZT2(I) = 0d0
4869 clin-4/25/03 add zt3(I) to track longitudinal positions of partons/strings:
4870          ZT3(I) = 0d0
4871          NP(I) = 0
4872  1002 CONTINUE
4873       DO 1003 I = 1, MUL
4874          ISTRG = LSTRG1(I)
4875          TAU7 = SQRT(FT5(I) ** 2 - GZ5(I) ** 2)
4876          ATAUI(ISTRG) = ATAUI(ISTRG) + TAU7
4877          ZT1(ISTRG) = ZT1(ISTRG) + GX5(I)
4878          ZT2(ISTRG) = ZT2(ISTRG) + GY5(I)
4879          ZT3(ISTRG) = ZT3(ISTRG) + GZ5(I)
4880          NP(ISTRG) = NP(ISTRG) + 1
4881  1003 CONTINUE
4882
4883       NSTR = NSP + NST + NSI
4884
4885 clin-7/03/01 correct averaging on transverse coordinates, no shift needed:
4886       if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
4887          DO 1004 I = 1, NSTR
4888             IF (NP(I) .NE. 0) THEN
4889                ATAUI(I) = ATAUI(I) / NP(I)
4890                ZT1(I) = ZT1(I) / NP(I)
4891                ZT2(I) = ZT2(I) / NP(I)
4892                ZT3(I) = ZT3(I) / NP(I)
4893             ENDIF
4894  1004    CONTINUE
4895          return
4896       endif
4897 clin-7/03/01-end
4898
4899       DO 1005 I = 1, NSTR
4900          IF (NP(I) .NE. 0) THEN
4901             ATAUI(I) = ATAUI(I) / NP(I)
4902             ZT1(I) = ZT1(I) / NP(I)
4903             ZT2(I) = ZT2(I) / NP(I)
4904             ZT3(I) = ZT3(I) / NP(I)
4905          ELSE
4906             IF (I .LE. NSP) THEN
4907                J = I
4908                ZT1(I) = dble(YP(1, J))
4909                ZT2(I) = dble(YP(2, J))
4910                ZT3(I) = 0d0
4911             ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN
4912                J = I - NSP
4913                ZT1(I) = dble(YT(1, J))
4914                ZT2(I) = dble(YT(2, J))
4915                ZT3(I) = 0d0
4916             ELSE
4917                J = I - NSP - NST
4918                ZT1(I) = 0.5d0*
4919      1              dble((YP(1, IASG(J, 1)) + YT(1, IASG(J, 2))))
4920                ZT2(I) = 0.5d0 *
4921      1              dble((YP(2, IASG(J, 1)) + YT(2, IASG(J, 2))))
4922                ZT3(I) = 0d0
4923             END IF
4924          END IF
4925  1005 CONTINUE
4926
4927 cbz6/28/99 flow1
4928       BB = HINT1(19)
4929       DO 1006 I = 1, NSTR
4930          IF (NP(I).NE.0) THEN
4931             SHIFT=0d0
4932          ELSE
4933             SHIFT=0.5d0*dble(BB)
4934          END IF
4935          IF (I .LE. NSP) THEN
4936             ZT1(I) = ZT1(I) + SHIFT
4937          ELSE IF (I .GT. NSP .AND. I .LE. NSP + NST) THEN
4938             ZT1(I) = ZT1(I) - SHIFT
4939          END IF
4940  1006 CONTINUE
4941 cbz6/28/99 flow1 end
4942 c
4943       RETURN
4944       END
4945
4946 clin-10/01/03 random number generator for f77:
4947 c      function RANART(NSEED)
4948 c      SAVE
4949 clin-4/2008 ran(nseed) is renamed to avoid conflict with system functions:
4950 c      ran=rand()
4951 c      ranart=rand(0)
4952 c     one may also use the following random number generator in PYTHIA/JETSET:
4953 c      ranart=rlu(0)
4954 c      return
4955 c      end
4956
4957 clin-3/2009
4958 c     Initialize hadron weights;
4959 c     Can add initial hadrons before the hadron cascade starts (but after ZPC).
4960       subroutine addhad
4961       PARAMETER (MAXSTR=150001,MAXR=1,xmd=1.8756)
4962       double precision  smearp,smearh
4963       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
4964       COMMON /ARPRC/ ITYPAR(MAXSTR),
4965      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
4966      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
4967      &     XMAR(MAXSTR)
4968       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
4969      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
4970      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
4971       COMMON /smearz/smearp,smearh
4972       COMMON/RNDF77/NSEED
4973       common /para8/ idpert,npertd,idxsec
4974       SAVE
4975
4976 c     All hadrons at the start of hadron cascade have the weight of 1
4977 c     except those inserted by the user in this subroutine:
4978       np0=IAINT2(1)
4979       DO i=1,np0
4980          dpertp(I)=1.
4981       ENDDO
4982 c     Specify number, species, weight, initial x,p,m for inserted hadrons here:
4983       nadd=0
4984       tau0=ARPAR1(1)
4985       DO 100 i=np0+1,np0+nadd
4986          ITYPAR(I)=42
4987 clin-5/2012 fix type mismatch:
4988 c         dpertp(I)=1d0/dble(nadd)
4989          dpertp(I)=1./float(nadd)
4990          GXAR(I)=5.*(1.-2.*RANART(NSEED))
4991          GYAR(I)=5.*(1.-2.*RANART(NSEED))
4992          GZAR(I)=2.*(1.-2.*RANART(NSEED))
4993          FTAR(I)=0.
4994          PXAR(I)=1.
4995          PYAR(I)=0.
4996          PZAR(I)=1.
4997          XMAR(I)=xmd
4998 c
4999          PEAR(I)=sqrt(PXAR(I)**2+PYAR(I)**2+PZAR(I)**2+XMAR(I)**2)
5000 clin-9/2012 determine rapidity more generally:
5001 c         RAP=0.5*alog((PEAR(I)+PZAR(I)+1e-5)/(PEAR(I)-PZAR(I)+1e-5))
5002          RAP=asinh(PZAR(I)/sqrt(XMAR(I)**2+PXAR(I)**2+PYAR(I)**2))
5003 c
5004          VX=PXAR(I)/PEAR(I)
5005          VY=PYAR(I)/PEAR(I)
5006 c.....give initial formation time shift and boost according to rapidity:
5007          TAUI=FTAR(I)+TAU0
5008          FTAR(I)=TAUI*COSH(RAP)
5009          GXAR(I)=GXAR(I)+VX*TAU0*COSH(RAP)
5010          GYAR(I)=GYAR(I)+VY*TAU0*COSH(RAP)
5011 c     Allow the intial z-position to be different from the Bjorken picture:
5012          GZAR(I)=TAUI*SINH(RAP)+GZAR(I)
5013 c         GZAR(I)=TAUI*SINH(RAP)
5014          zsmear=sngl(smearh)*(2.*RANART(NSEED)-1.)
5015          GZAR(I)=GZAR(I)+zsmear
5016  100  CONTINUE
5017       IAINT2(1)=IAINT2(1)+nadd
5018 c
5019       if(nadd.ge.1.and.idpert.ne.1.and.idpert.ne.2) then
5020          write(16,*) 'IDPERT must be 1 or 2 to add initial hadrons,
5021      1 set NPERTD to 0 if you do not need perturbative deuterons'
5022          stop
5023       endif
5024       if(IAINT2(1).gt.MAXSTR) then
5025          write(16,*) 'Too many initial hadrons, array size is exceeded!'
5026          stop
5027       endif
5028 c
5029       return
5030       end
5031
5032 clin-8/2014 define function asinh():
5033       FUNCTION asinh(x)
5034       SAVE
5035       if(x.gt.0) then
5036          ASINH=alog(x+sqrt(x**2+1.))
5037       else
5038 c     a la suggestion de YP Liu:
5039          ASINH=-alog(-x+sqrt(x**2+1.))
5040       endif
5041       return
5042       end