TAmpt/AMPT/linana.f

   1 c.................... linana.f
   2 c=======================================================================
   3 c     10/26/01 update freezeout positions in case of interactions:
   4 clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L
   5 c     as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade.
   6 cms
   7 cms   dlw & gsfs Commented out output file writing
   8 cms
   9       subroutine hbtout(nnew,nt,ntmax)
  10 c
  11       PARAMETER  (MAXSTR=150001,MAXR=1)
  12 clin-5/2008 give tolerance to regular particles (perturbative probability 1):
  13       PARAMETER  (oneminus=0.99999,oneplus=1.00001)
  14       dimension lastkp(MAXSTR), newkp(MAXSTR),xnew(3)
  15       common /para7/ ioscar,nsmbbbar,nsmmeson
  16 cc      SAVE /para7/
  17       COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
  18 cc      SAVE /hbt/
  19       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
  20 cc      SAVE /input1/
  21       COMMON   /AA/  R(3,MAXSTR)
  22 cc      SAVE /AA/
  23       COMMON   /BB/  P(3,MAXSTR)
  24 cc      SAVE /BB/
  25       COMMON   /CC/  E(MAXSTR)
  26 cc      SAVE /CC/
  27       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
  28 cc      SAVE /EE/
  29       common /lastt/itimeh,bimp
  30 cc      SAVE /lastt/
  31       COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
  32 cc      SAVE /tdecay/
  33       COMMON /AREVT/ IAEVT, IARUN, MISS
  34 cc      SAVE /AREVT/
  35       common/snn/efrm,npart1,npart2
  36 cc      SAVE /snn/
  37       COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
  38 cc      SAVE /HJGLBR/
  39       COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
  40       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
  41      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
  42      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
  43 clin-12/14/03:
  44       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
  45       EXTERNAL IARFLV, INVFLV
  46       common /para8/ idpert,npertd,idxsec
  47       SAVE
  48 c
  49       do 1001 i=1,max0(nlast,nnew)
  50          lastkp(i)=0
  51  1001 continue
  52       do 1002 i=1,nnew
  53          newkp(i)=0
  54  1002 continue
  55 c     for each of the particles, search the freezeout record (common /hbt/)
  56 c     to find & keep those which do not have interactions during this timestep:
  57       do 100 ip=1,nnew
  58          do 1004 iplast=1,nlast
  59             if(p(1,ip).eq.plast(1,iplast).and.
  60      1           p(2,ip).eq.plast(2,iplast).and.
  61      2           p(3,ip).eq.plast(3,iplast).and.
  62      3           e(ip).eq.plast(4,iplast).and.
  63      4           lb(ip).eq.lblast(iplast).and.
  64      5      dpertp(ip).eq.dplast(iplast).and.lastkp(iplast).eq.0) then
  65 clin-5/2008 modified below to the above in case we have perturbative particles:
  66 c     5           lastkp(iplast).eq.0) then
  67                deltat=nt*dt-xlast(4,iplast)
  68                ene=sqrt(plast(1,iplast)**2+plast(2,iplast)**2
  69      1              +plast(3,iplast)**2+plast(4,iplast)**2)
  70 c     xnew gives the coordinate if a particle free-streams to current time:
  71                do 1003 ii=1,3
  72                   xnew(ii)=xlast(ii,iplast)+plast(ii,iplast)/ene*deltat
  73  1003          continue
  74                   dr=sqrt((r(1,ip)-xnew(1))**2+(r(2,ip)-xnew(2))**2
  75      1              +(r(3,ip)-xnew(3))**2)
  76 c     find particles with dp=0 and dr<0.01, considered to be those
  77 c     without any interactions during this timestep,
  78 c     thus keep their last positions and time:
  79                if(dr.le.0.01) then
  80                   lastkp(iplast)=1
  81                   newkp(ip)=1
  82 c                  if(lb(ip).eq.41) then
  83 c                write(95,*) 'nt,ip,px,x=',nt,ip,p(1,ip),r(1,ip),ftsv(ip)
  84 c                write(95,*) 'xnew=',xnew(1),xnew(2),xnew(3),xlast(4,ip)
  85 c                  endif
  86 clin-5/2009 Take care of formation time of particles read in at nt=ntmax-1:
  87                   if(nt.eq.ntmax.and.ftsv(ip).gt.((ntmax-1)*dt))
  88      1                 xlast(4,iplast)=ftsv(ip)
  89                   goto 100
  90                endif
  91             endif
  92  1004    continue
  93  100  continue
  94 c     for current particles with interactions, fill their current info in
  95 c     the freezeout record (if that record entry needs not to be kept):
  96       do 150 ip=1,nnew
  97          if(newkp(ip).eq.0) then
  98             do 1005 iplast=1,nnew
  99                if(lastkp(iplast).eq.0) then
 100 ctest off: write collision info
 101 c                  if(lb(ip).eq.41) then
 102 c                     write(95,*) 'nt,lb(ip)=',nt,lb(ip)
 103 c                  write(95,*) '  last p=',plast(1,iplast),
 104 c     1 plast(2,iplast),plast(3,iplast),plast(4,iplast)
 105 c                  write(95,*) '  after p=',p(1,ip),p(2,ip),p(3,ip),e(ip)
 106 c                  write(95,*) 'after x=',r(1,ip),r(2,ip),r(3,ip),ftsv(ip)
 107 c                  endif
 108 c
 109                   xlast(1,iplast)=r(1,ip)
 110                   xlast(2,iplast)=r(2,ip)
 111                   xlast(3,iplast)=r(3,ip)
 112                   xlast(4,iplast)=nt*dt
 113 c
 114                   if(nt.eq.ntmax) then
 115 c     freezeout time for decay daughters at the last timestep
 116 c     needs to include the decay time of the parent:
 117                      if(tfdcy(ip).gt.(ntmax*dt+0.001)) then
 118                         xlast(4,iplast)=tfdcy(ip)
 119 c     freezeout time for particles unformed at the next-to-last timestep
 120 c     needs to be their formation time instead of (ntmax*dt):
 121                      elseif(ftsv(ip).gt.((ntmax-1)*dt)) then
 122                         xlast(4,iplast)=ftsv(ip)
 123                      endif
 124                   endif
 125                   plast(1,iplast)=p(1,ip)
 126                   plast(2,iplast)=p(2,ip)
 127                   plast(3,iplast)=p(3,ip)
 128                   plast(4,iplast)=e(ip)
 129                   lblast(iplast)=lb(ip)
 130                   lastkp(iplast)=1
 131 clin-5/2008:
 132                   dplast(iplast)=dpertp(ip)
 133                   goto 150
 134                endif
 135  1005       continue
 136          endif
 137  150  continue
 138 c     if the current particle list is shorter than the freezeout record,
 139 c     condense the last-collision record by filling new record from 1 to nnew,
 140 c     and label these entries as keep:
 141       if(nnew.lt.nlast) then
 142          do 170 iplast=1,nlast
 143             if(lastkp(iplast).eq.0) then
 144                do 1006 ip2=iplast+1,nlast
 145                   if(lastkp(ip2).eq.1) then
 146                      xlast(1,iplast)=xlast(1,ip2)
 147                      xlast(2,iplast)=xlast(2,ip2)
 148                      xlast(3,iplast)=xlast(3,ip2)
 149                      xlast(4,iplast)=xlast(4,ip2)
 150                      plast(1,iplast)=plast(1,ip2)
 151                      plast(2,iplast)=plast(2,ip2)
 152                      plast(3,iplast)=plast(3,ip2)
 153                      plast(4,iplast)=plast(4,ip2)
 154                      lblast(iplast)=lblast(ip2)
 155                      lastkp(iplast)=1
 156 clin-5/2008:
 157                      dplast(iplast)=dplast(ip2)
 158                      goto 170
 159                   endif
 160  1006          continue
 161             endif
 162  170     continue
 163       endif
 164       nlast=nnew
 165 ctest off look inside each NT timestep (for debugging purpose):
 166 c      do ip=1,nlast
 167 c         if(nt.eq.5000) then
 168 c            write(95,*) ' p ',nt,ip,INVFLV(lblast(ip)),plast(1,ip),
 169 c     1           plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip)
 170 c            write(95,*) '  x ',nt,ip,INVFLV(lblast(ip)),xlast(1,ip),
 171 c     1           xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip)
 172 c         endif
 173 c      enddo
 174 c
 175       if(nt.eq.ntmax) then
 176 clin-5/2008 find final number of perturbative particles (deuterons only):
 177          ndpert=0
 178          do ip=1,nlast
 179             if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
 180             else
 181                ndpert=ndpert+1
 182             endif
 183          enddo
 184 c
 185 c         write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2,
 186 c     1 NELP,NINP,NELT,NINTHJ
 187 c         write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
 188 c     1 NELP,NINP,NELT,NINTHJ
 189 clin-5/2008 write out perturbatively-produced particles (deuterons only):
 190 c        if(idpert.eq.1.or.idpert.eq.2)
 191 c     1        write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2,
 192 c     2        NELP,NINP,NELT,NINTHJ
 193 c         do 1007 ip=1,nlast
 194 clin-12/14/03   No formation time for spectator projectile or target nucleons,
 195 c     see ARINI1 in 'amptsub.f':
 196 clin-3/2009 To be consistent with new particles produced in hadron cascade
 197 c     that are limited by the time-resolution (DT) of the hadron cascade,
 198 c     freezeout time of spectator projectile or target nucleons is written as
 199 c     DT as they are read at the 1st timestep and then propagated to time DT:
 200 c            if(plast(1,ip).eq.0.and.plast(2,ip).eq.0
 201 c     1           .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1))
 202 c     2           .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
 203 clin-5/2008 perturbatively-produced particles (currently only deuterons)
 204 c     are written to ana/ampt_pert.dat (without the column for the mass);
 205 c     ana/ampt.dat has regularly-produced particles (including deuterons);
 206 c     these two sets of deuteron data are close to each other(but not the same
 207 c     because of the bias from triggering the perturbative production);
 208 c     ONLY use one data set for analysis to avoid double-counting:
 209 c               if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
 210 c                  write(16,200) INVFLV(lblast(ip)), plast(1,ip),
 211 c     1                 plast(2,ip),plast(3,ip),plast(4,ip),
 212 c     2                 xlast(1,ip),xlast(2,ip),xlast(3,ip),
 213 c     3                 xlast(4,ip)
 214 clin-12/14/03-end
 215 c               else
 216 c                  if(idpert.eq.1.or.idpert.eq.2) then
 217 c                     write(90,250) INVFLV(lblast(ip)), plast(1,ip),
 218 c     1                 plast(2,ip),plast(3,ip),
 219 c     2                 xlast(1,ip),xlast(2,ip),xlast(3,ip),
 220 c     3                 xlast(4,ip)
 221 c                  else
 222 c                     write(99,*) 'Unexpected perturbative particles'
 223 c                  endif
 224 c               endif
 225 c            elseif(amax1(abs(xlast(1,ip)),abs(xlast(2,ip)),
 226 c     1              abs(xlast(3,ip)),abs(xlast(4,ip))).lt.9999) then
 227 c               if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
 228 c            write(16,200) INVFLV(lblast(ip)), plast(1,ip),
 229 c     1           plast(2,ip),plast(3,ip),plast(4,ip),
 230 c     2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
 231 c               else
 232 c                  if(idpert.eq.1.or.idpert.eq.2) then
 233 c            write(90,250) INVFLV(lblast(ip)),plast(1,ip),
 234 c     1           plast(2,ip),plast(3,ip),
 235 c     2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
 236 c     3           dplast(ip)
 237 c                  else
 238 c                     write(99,*) 'Unexpected perturbative particles'
 239 c                  endif
 240 c              endif
 241 c            else
 242 c     change format for large numbers:
 243 c               if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
 244 c            write(16,201) INVFLV(lblast(ip)), plast(1,ip),
 245 c     1           plast(2,ip),plast(3,ip),plast(4,ip),
 246 c     2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
 247 c               else
 248 c                  if(idpert.eq.1.or.idpert.eq.2) then
 249 c                     write(90,251) INVFLV(lblast(ip)), plast(1,ip),
 250 c     1           plast(2,ip),plast(3,ip),
 251 c     2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
 252 c     3           dplast(ip)
 253 c                  else
 254 c                     write(99,*) 'Unexpected perturbative particles'
 255 c                  endif
 256 c               endif
 257 c            endif
 258 c 1007    continue
 259          if(ioscar.eq.1) call hoscar
 260       endif
 261  190  format(3(i7),f10.4,5x,6(i4))
 262 clin-3/2009 improve the output accuracy of Pz
 263  200  format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,f8.2))
 264  201  format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,e8.2))
 265  250  format(I5,2(1x,f8.3),1x,f10.3,2(1x,f7.1),1x,f8.2,1x,f7.2,1x,e10.4)
 266  251  format(I5,2(1x,f8.3),1x,f10.3,4(1x,e8.2),1x,e10.4)
 267 c
 268         return
 269         end
 270
 271 c=======================================================================
 272         SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1)
 273 c
 274         IMPLICIT DOUBLE PRECISION(D)
 275         DOUBLE PRECISION  enenew, pxnew, pynew, pznew
 276         DOUBLE PRECISION  de0, beta2, gam
 277         common /lora/ enenew, pxnew, pynew, pznew
 278 cc      SAVE /lora/
 279         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
 280 cc      SAVE /HPARNT/
 281         common /decom/ptwo(2,5)
 282 cc      SAVE /decom/
 283         COMMON/RNDF77/NSEED
 284 cc      SAVE /RNDF77/
 285         COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
 286         common/embed/iembed,pxqembd,pyqembd,xembd,yembd
 287         SAVE
 288 c
 289         itq1=itq1
 290         dcth=dble(RANART(NSEED))*2.d0-1.d0
 291         dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0
 292 clin-6/2009 Added if embedding a high-Pt quark pair after string melting:
 293         if(iembed.eq.1) then
 294 c     Decompose the parent high-Pt pion to q and qbar with an internal momentum
 295 c     parallel to the pion direction so that one parton has ~the same hight Pt
 296 c     and the other parton has a very soft Pt:
 297 c     Note: htop() decomposes a meson to q as it(1) followed by qbar as it(2):
 298            if(i.eq.(natt-1).or.i.eq.natt) then
 299               dcth=0.d0
 300               dphi=dble(acos(pxqembd/sqrt(pxqembd**2+pyqembd**2)))
 301               if(pyqembd.lt.0) dphi=HIPR1(40)*2.d0-dphi
 302            endif
 303         endif
 304 c
 305         ds=dble(xm0)**2
 306         dpcm=dsqrt((ds-dble(ptwo(1,5)+ptwo(2,5))**2)
 307      1 *(ds-dble(ptwo(1,5)-ptwo(2,5))**2)/ds/4d0)
 308         dpz=dpcm*dcth
 309         dpx=dpcm*dsqrt(1.d0-dcth**2)*dcos(dphi)
 310         dpy=dpcm*dsqrt(1.d0-dcth**2)*dsin(dphi)
 311         de1=dsqrt(dble(ptwo(1,5))**2+dpcm**2)
 312         de2=dsqrt(dble(ptwo(2,5))**2+dpcm**2)
 313 c
 314       de0=dsqrt(dble(px0)**2+dble(py0)**2+dble(pz0)**2+dble(xm0)**2)
 315         dbex=dble(px0)/de0
 316         dbey=dble(py0)/de0
 317         dbez=dble(pz0)/de0
 318 c     boost the reference frame up by beta (pznew=gam(pz+beta e)):
 319       beta2 = dbex ** 2 + dbey ** 2 + dbez ** 2
 320       gam = 1.d0 / dsqrt(1.d0 - beta2)
 321       if(beta2.ge.0.9999999999999d0) then
 322          write(6,*) '1',dbex,dbey,dbez,beta2,gam
 323       endif
 324 c
 325       call lorenza(de1,dpx,dpy,dpz,-dbex,-dbey,-dbez)
 326         ptwo(1,1)=sngl(pxnew)
 327         ptwo(1,2)=sngl(pynew)
 328         ptwo(1,3)=sngl(pznew)
 329         ptwo(1,4)=sngl(enenew)
 330       call lorenza(de2,-dpx,-dpy,-dpz,-dbex,-dbey,-dbez)
 331         ptwo(2,1)=sngl(pxnew)
 332         ptwo(2,2)=sngl(pynew)
 333         ptwo(2,3)=sngl(pznew)
 334         ptwo(2,4)=sngl(enenew)
 335 c
 336       RETURN
 337       END
 338
 339 c=======================================================================
 340       SUBROUTINE HTOP
 341 c
 342       PARAMETER (MAXSTR=150001)
 343       PARAMETER (MAXPTN=400001)
 344       PARAMETER (MAXIDL=4001)
 345       DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
 346       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
 347      1     GXSGS,GYSGS,GZSGS,FTSGS
 348       dimension it(4)
 349       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
 350 cwei      DOUBLE PRECISION PATT
 351 cc      SAVE /HMAIN2/
 352       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
 353 cc      SAVE /HMAIN1/
 354       COMMON /PARA1/ MUL
 355 cc      SAVE /PARA1/
 356       COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 357      &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 358      &     XMASS0(MAXPTN), ITYP0(MAXPTN)
 359 cc      SAVE /prec1/
 360       COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
 361 cc      SAVE /ilist7/
 362       COMMON /ARPRC/ ITYPAR(MAXSTR),
 363      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
 364      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
 365      &     XMAR(MAXSTR)
 366 cc      SAVE /ARPRC/
 367       common /decom/ptwo(2,5)
 368 cc      SAVE /decom/
 369       COMMON/RNDF77/NSEED
 370 cc      SAVE /RNDF77/
 371       COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
 372      &     GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
 373      &     PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
 374      &     XMN(MAXIDL)
 375 cc      SAVE /NOPREC/
 376       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
 377 cc      SAVE /HPARNT/
 378 c     7/20/01: use double precision
 379 c     otherwise sometimes beta>1 and gamma diverge in lorenza():
 380       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
 381      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
 382      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
 383      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
 384 cc      SAVE /SOFT/
 385       common/anim/nevent,isoft,isflag,izpc
 386 cc      SAVE /anim/
 387       DOUBLE PRECISION  vxp0,vyp0,vzp0
 388       common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
 389 cc      SAVE /precpa/
 390       common /para7/ ioscar,nsmbbbar,nsmmeson
 391       COMMON /AREVT/ IAEVT, IARUN, MISS
 392       SAVE
 393 c
 394         npar=0
 395         nnozpc=0
 396 clin-5b/2008 calculate the number of hadrons to be converted to q/qbar:
 397         if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
 398      1       then
 399            nsmbbbar=0
 400            nsmmeson=0
 401            do i=1,natt
 402               id=ITYPAR(i)
 403               idabs=iabs(id)
 404               i2=MOD(idabs/10,10)
 405               if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
 406      1             .ge.(HINT1(1)/2*0.99).and.
 407      2             ((id.eq.2112).or.(id.eq.2212))) then
 408 c     proj or targ nucleons without interactions, do not enter ZPC:
 409               elseif(idabs.gt.1000.and.i2.ne.0) then
 410 c     baryons to be converted to q/qbar:
 411                  nsmbbbar=nsmbbbar+1
 412               elseif((idabs.gt.100.and.idabs.lt.1000)
 413      1                .or.idabs.gt.10000) then
 414 c     mesons to be converted to q/qbar:
 415                  nsmmeson=nsmmeson+1
 416               endif
 417            enddo
 418
 419 clin-6/2009:
 420            if(ioscar.eq.2.or.ioscar.eq.3) then
 421               write(92,*) iaevt,miss,3*nsmbbbar+2*nsmmeson,
 422      1             nsmbbbar,nsmmeson,natt,natt-nsmbbbar-nsmmeson
 423            endif
 424 c           write(92,*) iaevt, 3*nsmbbbar+2*nsmmeson
 425 c           write(92,*) ' event#, total # of initial partons after string
 426 c     1 melting'
 427 c           write(92,*) 'String melting converts ',nsmbbbar, ' baryons &'
 428 c     1, nsmmeson, ' mesons'
 429 c           write(92,*) 'Total # of initial particles= ',natt
 430 c           write(92,*) 'Total # of initial particles (gamma,e,muon,...)
 431 c     1 not entering ZPC= ',natt-nsmbbbar-nsmmeson
 432         endif
 433 clin-5b/2008-over
 434         do 100 i=1,natt
 435            id=ITYPAR(i)
 436            idabs=iabs(id)
 437            i4=MOD(idabs/1000,10)
 438            i3=MOD(idabs/100,10)
 439            i2=MOD(idabs/10,10)
 440            i1=MOD(idabs,10)
 441            rnum=RANART(NSEED)
 442            ftime=0.197*PEAR(i)/(PXAR(i)**2+PYAR(i)**2+XMAR(i)**2)
 443            inozpc=0
 444            it(1)=0
 445            it(2)=0
 446            it(3)=0
 447            it(4)=0
 448 c
 449            if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
 450      1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then
 451 c     proj or targ nucleons without interactions, do not enter ZPC:
 452               inozpc=1
 453            elseif(idabs.gt.1000.and.i2.ne.0) then
 454 c     baryons:
 455               if(((i4.eq.1.or.i4.eq.2).and.i4.eq.i3)
 456      1 .or.(i4.eq.3.and.i3.eq.3)) then
 457                  if(i1.eq.2) then
 458                     if(rnum.le.(1./2.)) then
 459                        it(1)=i4
 460                        it(2)=i3*1000+i2*100+1
 461                     elseif(rnum.le.(2./3.)) then
 462                        it(1)=i4
 463                        it(2)=i3*1000+i2*100+3
 464                     else
 465                        it(1)=i2
 466                        it(2)=i4*1000+i3*100+3
 467                     endif
 468                  elseif(i1.eq.4) then
 469                     if(rnum.le.(2./3.)) then
 470                        it(1)=i4
 471                        it(2)=i3*1000+i2*100+3
 472                     else
 473                        it(1)=i2
 474                        it(2)=i4*1000+i3*100+3
 475                     endif
 476                  endif
 477               elseif(i4.eq.1.or.i4.eq.2) then
 478                  if(i1.eq.2) then
 479                     if(rnum.le.(1./2.)) then
 480                        it(1)=i2
 481                        it(2)=i4*1000+i3*100+1
 482                     elseif(rnum.le.(2./3.)) then
 483                        it(1)=i2
 484                        it(2)=i4*1000+i3*100+3
 485                     else
 486                        it(1)=i4
 487                        it(2)=i3*1000+i2*100+3
 488                     endif
 489                  elseif(i1.eq.4) then
 490                     if(rnum.le.(2./3.)) then
 491                        it(1)=i2
 492                        it(2)=i4*1000+i3*100+3
 493                     else
 494                        it(1)=i4
 495                        it(2)=i3*1000+i2*100+3
 496                     endif
 497                  endif
 498               elseif(i4.ge.3) then
 499                  it(1)=i4
 500                  if(i3.lt.i2) then
 501                     it(2)=i2*1000+i3*100+1
 502                  else
 503                     it(2)=i3*1000+i2*100+3
 504                  endif
 505               endif
 506 c       antibaryons:
 507               if(id.lt.0) then
 508                  it(1)=-it(1)
 509                  it(2)=-it(2)
 510               endif
 511 c     isoft=4or5 decompose diquark flavor it(2) to two quarks it(3)&(4):
 512               if(isoft.eq.4.or.isoft.eq.5) then
 513                  it(3)=MOD(it(2)/1000,10)
 514                  it(4)=MOD(it(2)/100,10)
 515               endif
 516
 517            elseif((idabs.gt.100.and.idabs.lt.1000)
 518      1 .or.idabs.gt.10000) then
 519 c     mesons:
 520               if(i3.eq.i2) then
 521                  if(i3.eq.1.or.i3.eq.2) then
 522                     if(rnum.le.0.5) then
 523                        it(1)=1
 524                        it(2)=-1
 525                     else
 526                        it(1)=2
 527                        it(2)=-2
 528                     endif
 529                  else
 530                     it(1)=i3
 531                     it(2)=-i3
 532                  endif
 533               else
 534                  if((isign(1,id)*(-1)**i3).eq.1) then
 535                     it(1)=i3
 536                     it(2)=-i2
 537                  else
 538                     it(1)=i2
 539                     it(2)=-i3
 540                  endif
 541               endif
 542            else
 543 c     save other particles (leptons and photons) outside of ZPC:
 544               inozpc=1
 545            endif
 546 c
 547            if(inozpc.eq.1) then
 548               NJSGS(i)=0
 549               nnozpc=nnozpc+1
 550               itypn(nnozpc)=ITYPAR(i)
 551               pxn(nnozpc)=PXAR(i)
 552               pyn(nnozpc)=PYAR(i)
 553               pzn(nnozpc)=PZAR(i)
 554               een(nnozpc)=PEAR(i)
 555               xmn(nnozpc)=XMAR(i)
 556               gxn(nnozpc)=GXAR(i)
 557               gyn(nnozpc)=GYAR(i)
 558               gzn(nnozpc)=GZAR(i)
 559               ftn(nnozpc)=FTAR(i)
 560            else
 561               NJSGS(i)=2
 562               ptwo(1,5)=ulmass(it(1))
 563               ptwo(2,5)=ulmass(it(2))
 564               call decomp(patt(i,1),patt(i,2),patt(i,3),XMAR(i),i,it(1))
 565               ipamax=2
 566               if((isoft.eq.4.or.isoft.eq.5)
 567      1 .and.iabs(it(2)).gt.1000) ipamax=1
 568               do 1001 ipar=1,ipamax
 569                  npar=npar+1
 570                  ityp0(npar)=it(ipar)
 571                  px0(npar)=dble(ptwo(ipar,1))
 572                  py0(npar)=dble(ptwo(ipar,2))
 573                  pz0(npar)=dble(ptwo(ipar,3))
 574                  e0(npar)=dble(ptwo(ipar,4))
 575                  xmass0(npar)=dble(ptwo(ipar,5))
 576                  gx0(npar)=dble(GXAR(i))
 577                  gy0(npar)=dble(GYAR(i))
 578                  gz0(npar)=dble(GZAR(i))
 579                  ft0(npar)=dble(ftime)
 580                  lstrg0(npar)=i
 581                  lpart0(npar)=ipar
 582                  vxp0(npar)=dble(patt(i,1)/patt(i,4))
 583                  vyp0(npar)=dble(patt(i,2)/patt(i,4))
 584                  vzp0(npar)=dble(patt(i,3)/patt(i,4))
 585  1001     continue
 586 cyy 200      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
 587 cyy 201      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
 588 c
 589               if((isoft.eq.4.or.isoft.eq.5)
 590      1 .and.iabs(it(2)).gt.1000) then
 591                  NJSGS(i)=3
 592                  xmdq=ptwo(2,5)
 593                  ptwo(1,5)=ulmass(it(3))
 594                  ptwo(2,5)=ulmass(it(4))
 595 c     8/19/02 avoid actual argument in common blocks of DECOMP:
 596 c                 call decomp(ptwo(2,1),ptwo(2,2),ptwo(2,3),xmdq)
 597              ptwox=ptwo(2,1)
 598              ptwoy=ptwo(2,2)
 599              ptwoz=ptwo(2,3)
 600              call decomp(ptwox,ptwoy,ptwoz,xmdq,i,it(1))
 601 c
 602                  do 1002 ipar=1,2
 603                     npar=npar+1
 604                     ityp0(npar)=it(ipar+2)
 605                     px0(npar)=dble(ptwo(ipar,1))
 606                     py0(npar)=dble(ptwo(ipar,2))
 607                     pz0(npar)=dble(ptwo(ipar,3))
 608                     e0(npar)=dble(ptwo(ipar,4))
 609                     xmass0(npar)=dble(ptwo(ipar,5))
 610                     gx0(npar)=dble(GXAR(i))
 611                     gy0(npar)=dble(GYAR(i))
 612                     gz0(npar)=dble(GZAR(i))
 613                     ft0(npar)=dble(ftime)
 614                     lstrg0(npar)=i
 615                     lpart0(npar)=ipar+1
 616                     vxp0(npar)=dble(patt(i,1)/patt(i,4))
 617                     vyp0(npar)=dble(patt(i,2)/patt(i,4))
 618                     vzp0(npar)=dble(patt(i,3)/patt(i,4))
 619  1002        continue
 620               endif
 621 c
 622            endif
 623  100        continue
 624       MUL=NPAR
 625 c
 626 clin-5b/2008:
 627       if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
 628      1     then
 629          if((natt-nsmbbbar-nsmmeson).ne.nnozpc)
 630      1        write(92,*) 'Problem with the total # of initial particles
 631      2 (gamma,e,muon,...) not entering ZPC'
 632          if((3*nsmbbbar+2*nsmmeson).ne.npar)
 633      1        write(92,*) 'Problem with the total # of initial partons
 634      2 after string melting'
 635       endif
 636 c
 637       RETURN
 638       END
 639
 640 c=======================================================================
 641       SUBROUTINE PTOH
 642 c
 643       PARAMETER (MAXSTR=150001)
 644       DOUBLE PRECISION  gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
 645       DOUBLE PRECISION  gxp0,gyp0,gzp0,ft0fom,drlocl
 646       DOUBLE PRECISION  enenew, pxnew, pynew, pznew, beta2, gam
 647       DOUBLE PRECISION  ftavg0,gxavg0,gyavg0,gzavg0,bex,bey,bez
 648       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
 649      1     GXSGS,GYSGS,GZSGS,FTSGS
 650       DOUBLE PRECISION  xmdiag,px1,py1,pz1,e1,px2,py2,pz2,e2,
 651      1     px3,py3,pz3,e3,xmpair,etot
 652       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
 653      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
 654 cc      SAVE /loclco/
 655       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
 656 cc      SAVE /HMAIN1/
 657       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
 658 cc      SAVE /HMAIN2/
 659       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
 660      &     K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
 661      &     PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
 662 cc      SAVE /HJJET2/
 663       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
 664 cc      SAVE /ARPRNT/
 665       COMMON /ARPRC/ ITYPAR(MAXSTR),
 666      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
 667      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
 668      &     XMAR(MAXSTR)
 669 cc      SAVE /ARPRC/
 670       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
 671      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
 672      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
 673      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
 674 cc      SAVE /SOFT/
 675       COMMON/RNDF77/NSEED
 676 cc      SAVE /RNDF77/
 677       common/anim/nevent,isoft,isflag,izpc
 678 cc      SAVE /anim/
 679       common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
 680 cc      SAVE /prtn23/
 681       common /nzpc/nattzp
 682 cc      SAVE /nzpc/
 683       common /lora/ enenew, pxnew, pynew, pznew
 684 cc      SAVE /lora/
 685       COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
 686 cc      SAVE /LUDAT1A/
 687 clin 4/19/2006
 688       common /lastt/itimeh,bimp
 689       COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
 690       COMMON /AREVT/ IAEVT, IARUN, MISS
 691       common /para7/ ioscar,nsmbbbar,nsmmeson
 692 c
 693       dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR)
 694 cwei      DOUBLE PRECISION PATT
 695       SAVE
 696 c
 697       call coales
 698 c     obtain particle mass here without broadening by Breit-Wigner width:
 699       mstj24=MSTJ(24)
 700       MSTJ(24)=0
 701         nuudd=0
 702         npich=0
 703         nrhoch=0
 704       ppi0=1.
 705       prho0=0.
 706 c     determine hadron flavor except (pi0,rho0,eta,omega):
 707       DO 1001 ISG = 1, NSG
 708            if(NJSGS(ISG).ne.0) then
 709               NATT=NATT+1
 710               K1=K2SGS(ISG,1)
 711               k1abs=iabs(k1)
 712               PX1=PXSGS(ISG,1)
 713               PY1=PYSGS(ISG,1)
 714               PZ1=PZSGS(ISG,1)
 715               K2=K2SGS(ISG,2)
 716               k2abs=iabs(k2)
 717               PX2=PXSGS(ISG,2)
 718               PY2=PYSGS(ISG,2)
 719               PZ2=PZSGS(ISG,2)
 720 c     5/02/01 try lowest spin states as first choices,
 721 c     i.e. octet baryons and pseudoscalar mesons (ibs=2*baryonspin+1):
 722               e1=PESGS(ISG,1)
 723               e2=PESGS(ISG,2)
 724               xmpair=dsqrt((e1+e2)**2-(px1+px2)**2-(py1+py2)**2
 725      1 -(pz1+pz2)**2)
 726               ibs=2
 727               imspin=0
 728               if(k1.eq.-k2.and.iabs(k1).le.2.
 729      1           and.NJSGS(ISG).eq.2) then
 730                nuudd=nuudd+1
 731                xmdiag(nuudd)=xmpair
 732                ndiag(nuudd)=natt
 733             endif
 734               K3=0
 735               if((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3) then
 736                K3=K2SGS(ISG,3)
 737                k3abs=iabs(k3)
 738                PX3=PXSGS(ISG,3)
 739                PY3=PYSGS(ISG,3)
 740                PZ3=PZSGS(ISG,3)
 741                e3=PESGS(ISG,3)
 742                xmpair=dsqrt((e1+e2+e3)**2-(px1+px2+px3)**2
 743      1              -(py1+py2+py3)**2-(pz1+pz2+pz3)**2)
 744               endif
 745 c*****     isoft=3 baryon decomposition is different:
 746               if(isoft.eq.3.and.
 747      1           (k1abs.gt.1000.or.k2abs.gt.1000)) then
 748                if(k1abs.gt.1000) then
 749                   kdq=k1abs
 750                   kk=k2abs
 751                else
 752                   kdq=k2abs
 753                   kk=k1abs
 754                endif
 755                ki=MOD(kdq/1000,10)
 756                kj=MOD(kdq/100,10)
 757                if(MOD(kdq,10).eq.1) then
 758                   idqspn=0
 759                else
 760                   idqspn=1
 761                endif
 762 c
 763                if(kk.gt.ki) then
 764                   ktemp=kk
 765                   kk=kj
 766                   kj=ki
 767                   ki=ktemp
 768                elseif(kk.gt.kj) then
 769                   ktemp=kk
 770                   kk=kj
 771                   kj=ktemp
 772                endif
 773 c
 774                if(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
 775                   if(idqspn.eq.0) then
 776                      kf=1000*ki+100*kk+10*kj+ibs
 777                   else
 778                      kf=1000*ki+100*kj+10*kk+ibs
 779                   endif
 780                elseif(ki.eq.kj.and.ki.eq.kk) then
 781 c     can only be decuplet baryons:
 782                   kf=1000*ki+100*kj+10*kk+4
 783                else
 784                   kf=1000*ki+100*kj+10*kk+ibs
 785                endif
 786 c     form a decuplet baryon if the q+diquark mass is closer to its mass
 787 c     (and if the diquark has spin 1):
 788 cc     for now only include Delta, which is present in ART:
 789 cc                 if(idqspn.eq.1.and.MOD(kf,10).eq.2) then
 790                if(kf.eq.2112.or.kf.eq.2212) then
 791                   if(abs(sngl(xmpair)-ULMASS(kf)).gt.
 792      1                 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
 793                endif
 794                if(k1.lt.0) kf=-kf
 795 clin-6/22/01 isoft=4or5 baryons:
 796               elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
 797      1              then
 798                if(k1abs.gt.k2abs) then
 799                   ki=k1abs
 800                   kk=k2abs
 801                else
 802                   ki=k2abs
 803                   kk=k1abs
 804                endif
 805                if(k3abs.gt.ki) then
 806                   kj=ki
 807                   ki=k3abs
 808                elseif(k3abs.lt.kk) then
 809                   kj=kk
 810                   kk=k3abs
 811                else
 812                   kj=k3abs
 813                endif
 814 c
 815                if(ki.eq.kj.and.ki.eq.kk) then
 816 c     can only be decuplet baryons (Delta-,++, Omega):
 817                   ibs=4
 818                   kf=1000*ki+100*kj+10*kk+ibs
 819                elseif(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
 820 c     form Lambda or Sigma according to 3-quark mass,
 821 c     for now neglect decuplet (Sigma*0 etc) which is absent in ART:
 822                   ibs=2
 823                   kf1=1000*ki+100*kj+10*kk+ibs
 824                   kf2=1000*ki+100*kk+10*kj+ibs
 825                   kf=kf1
 826                   if(abs(sngl(xmpair)-ULMASS(kf1)).gt.
 827      1                 abs(sngl(xmpair)-ULMASS(kf2))) kf=kf2
 828                else
 829                   ibs=2
 830                   kf=1000*ki+100*kj+10*kk+ibs
 831 cc     for now only include Delta0,+ as decuplets, which are present in ART:
 832                   if(kf.eq.2112.or.kf.eq.2212) then
 833                      if(abs(sngl(xmpair)-ULMASS(kf)).gt.
 834      1                    abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
 835                   endif
 836                endif
 837                if(k1.lt.0) kf=-kf
 838 c*****     mesons:
 839               else
 840                if(k1abs.eq.k2abs) then
 841                   if(k1abs.le.2) then
 842 c     treat diagonal mesons later in the subroutine:
 843                      kf=0
 844                   elseif(k1abs.le.3) then
 845 c     do not form eta', only form phi from s-sbar, since no eta' in ART:
 846                      kf=333
 847                   else
 848                      kf=100*k1abs+10*k1abs+2*imspin+1
 849                   endif
 850                else
 851                   if(k1abs.gt.k2abs) then
 852                      kmax=k1abs
 853                      kmin=k2abs
 854                   elseif(k1abs.lt.k2abs) then
 855                      kmax=k2abs
 856                      kmin=k1abs
 857                   endif
 858                   kf=(100*kmax+10*kmin+2*imspin+1)
 859      1                 *isign(1,k1+k2)*(-1)**kmax
 860 c     form a vector meson if the q+qbar mass is closer to its mass:
 861                   if(MOD(iabs(kf),10).eq.1) then
 862                      if(abs(sngl(xmpair)-ULMASS(iabs(kf))).gt.
 863      1                    abs(sngl(xmpair)-ULMASS(iabs(kf)+2)))
 864      2                    kf=(iabs(kf)+2)*isign(1,kf)
 865                   endif
 866                endif
 867               endif
 868               ITYPAR(NATT)=kf
 869               KATT(NATT,1)=kf
 870             if(iabs(kf).eq.211) then
 871                npich=npich+1
 872             elseif(iabs(kf).eq.213) then
 873                nrhoch=nrhoch+1
 874             endif
 875            endif
 876  1001   CONTINUE
 877 c     assume Npi0=(Npi+ + Npi-)/2, Nrho0=(Nrho+ + Nrho-)/2 on the average:
 878         if(nuudd.ne.0) then
 879          ppi0=float(npich/2)/float(nuudd)
 880          prho0=float(nrhoch/2)/float(nuudd)
 881       endif
 882 c     determine diagonal mesons (pi0,rho0,eta and omega) from uubar/ddbar:
 883       npi0=0
 884       DO 1002 ISG = 1, NSG
 885          if(K2SGS(ISG,1).eq.-K2SGS(ISG,2)
 886      1        .and.iabs(K2SGS(ISG,1)).le.2.and.NJSGS(ISG).eq.2) then
 887             if(RANART(NSEED).le.ppi0) npi0=npi0+1
 888          endif
 889  1002 CONTINUE
 890 c
 891       if(nuudd.gt.1) then
 892          call index1(MAXSTR,nuudd,xmdiag,indx)
 893       else
 894          indx(1)=1
 895       end if
 896 c
 897       DO 1003 ix=1,nuudd
 898          iuudd=indx(ix)
 899          inatt=ndiag(iuudd)
 900          if(ix.le.npi0) then
 901             kf=111
 902          elseif(RANART(NSEED).le.(prho0/(1-ppi0+0.00001))) then
 903             kf=113
 904          else
 905 c     at T=150MeV, thermal weights for eta and omega(spin1) are about the same:
 906             if(RANART(NSEED).le.0.5) then
 907                kf=221
 908             else
 909                kf=223
 910             endif
 911          endif
 912          ITYPAR(inatt)=kf
 913          KATT(inatt,1)=kf
 914  1003 CONTINUE
 915 c  determine hadron formation time, position and momentum:
 916       inatt=0
 917 clin-6/2009 write out parton info after coalescence:
 918       if(ioscar.eq.3) then
 919          WRITE (85, 395) IAEVT, 3*nsmbbbar+2*nsmmeson,nsmbbbar,nsmmeson,
 920      1     bimp, NELP,NINP,NELT,NINTHJ,MISS
 921       endif
 922  395  format(4I8,f10.4,5I5)
 923 c
 924       DO 1006 ISG = 1, NSG
 925            if(NJSGS(ISG).ne.0) then
 926             inatt=inatt+1
 927               K1=K2SGS(ISG,1)
 928               k1abs=iabs(k1)
 929               PX1=PXSGS(ISG,1)
 930               PY1=PYSGS(ISG,1)
 931               PZ1=PZSGS(ISG,1)
 932               K2=K2SGS(ISG,2)
 933               k2abs=iabs(k2)
 934               PX2=PXSGS(ISG,2)
 935               PY2=PYSGS(ISG,2)
 936               PZ2=PZSGS(ISG,2)
 937               e1=PESGS(ISG,1)
 938               e2=PESGS(ISG,2)
 939 c
 940               if(NJSGS(ISG).eq.2) then
 941                PXAR(inatt)=sngl(px1+px2)
 942                PYAR(inatt)=sngl(py1+py2)
 943                PZAR(inatt)=sngl(pz1+pz2)
 944                PATT(inatt,1)=PXAR(inatt)
 945                PATT(inatt,2)=PYAR(inatt)
 946                PATT(inatt,3)=PZAR(inatt)
 947                etot=e1+e2
 948               elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
 949      1              then
 950                PX3=PXSGS(ISG,3)
 951                PY3=PYSGS(ISG,3)
 952                PZ3=PZSGS(ISG,3)
 953                e3=PESGS(ISG,3)
 954                PXAR(inatt)=sngl(px1+px2+px3)
 955                PYAR(inatt)=sngl(py1+py2+py3)
 956                PZAR(inatt)=sngl(pz1+pz2+pz3)
 957                PATT(inatt,1)=PXAR(inatt)
 958                PATT(inatt,2)=PYAR(inatt)
 959                PATT(inatt,3)=PZAR(inatt)
 960                etot=e1+e2+e3
 961               endif
 962               XMAR(inatt)=ULMASS(ITYPAR(inatt))
 963               PEAR(inatt)=sqrt(PXAR(inatt)**2+PYAR(inatt)**2
 964      1           +PZAR(inatt)**2+XMAR(inatt)**2)
 965               PATT(inatt,4)=PEAR(inatt)
 966               EATT=EATT+PEAR(inatt)
 967             ipartn=NJSGS(ISG)
 968             DO 1004 i=1,ipartn
 969                ftp(i)=ftsgs(isg,i)
 970                gxp(i)=gxsgs(isg,i)
 971                gyp(i)=gysgs(isg,i)
 972                gzp(i)=gzsgs(isg,i)
 973                pxp(i)=pxsgs(isg,i)
 974                pyp(i)=pysgs(isg,i)
 975                pzp(i)=pzsgs(isg,i)
 976                pmp(i)=pmsgs(isg,i)
 977                pep(i)=pesgs(isg,i)
 978  1004       CONTINUE
 979             call locldr(ipartn,drlocl)
 980 c
 981             tau0=ARPAR1(1)
 982             ftavg0=ft0fom+dble(tau0)
 983             gxavg0=0d0
 984             gyavg0=0d0
 985             gzavg0=0d0
 986             DO 1005 i=1,ipartn
 987                gxavg0=gxavg0+gxp0(i)/ipartn
 988                gyavg0=gyavg0+gyp0(i)/ipartn
 989                gzavg0=gzavg0+gzp0(i)/ipartn
 990  1005       CONTINUE
 991             bex=dble(PXAR(inatt))/etot
 992             bey=dble(PYAR(inatt))/etot
 993             bez=dble(PZAR(inatt))/etot
 994             beta2 = bex ** 2 + bey ** 2 + bez ** 2
 995             gam = 1.d0 / dsqrt(1.d0 - beta2)
 996             if(beta2.ge.0.9999999999999d0) then
 997                write(6,*) '2',bex,bey,bez,beta2,gam
 998             endif
 999 c
1000             call lorenza(ftavg0,gxavg0,gyavg0,gzavg0,-bex,-bey,-bez)
1001               GXAR(inatt)=sngl(pxnew)
1002               GYAR(inatt)=sngl(pynew)
1003               GZAR(inatt)=sngl(pznew)
1004               FTAR(inatt)=sngl(enenew)
1005 clin 4/19/2006 write out parton info after coalescence:
1006               if(ioscar.eq.3) then
1007                  WRITE (85, 312) K2SGS(ISG,1),px1,py1,pz1,PMSGS(ISG,1),
1008      1                inatt,katt(inatt,1)
1009                  WRITE (85, 312) K2SGS(ISG,2),px2,py2,pz2,PMSGS(ISG,2),
1010      1                inatt,katt(inatt,1)
1011                  if(NJSGS(ISG).eq.3) WRITE (85, 312) K2SGS(ISG,3),
1012      1                px3,py3,pz3,PMSGS(ISG,3),inatt,katt(inatt,1)
1013               endif
1014  312       FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6)
1015 c
1016            endif
1017  1006   CONTINUE
1018 c     number of hadrons formed from partons inside ZPC:
1019       nattzp=natt
1020       MSTJ(24)=mstj24
1021 c
1022       RETURN
1023       END
1024
1025 c=======================================================================
1026       SUBROUTINE coales
1027
1028       PARAMETER (MAXSTR=150001)
1029       IMPLICIT DOUBLE PRECISION(D)
1030       DOUBLE PRECISION  gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
1031       DIMENSION IOVER(MAXSTR),dp1(2:3),dr1(2:3)
1032       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
1033      1     GXSGS,GYSGS,GZSGS,FTSGS
1034       double precision  dpcoal,drcoal,ecritl
1035       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1036      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1037      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1038      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1039 cc      SAVE /SOFT/
1040       common /coal/dpcoal,drcoal,ecritl
1041 cc      SAVE /coal/
1042       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1043      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1044 cc      SAVE /loclco/
1045       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
1046      &     K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
1047      &     PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
1048 cc      SAVE /HJJET2/
1049       SAVE
1050 c
1051       do 1001 ISG=1, NSG
1052          IOVER(ISG)=0
1053  1001 continue
1054 C1     meson q coalesce with all available qbar:
1055       do 150 ISG=1,NSG
1056          if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 150
1057 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1058          if(K2SGS(ISG,1).lt.0) then
1059             write(6,*) 'Antiquark appears in quark loop; stop'
1060             stop
1061          endif
1062 c
1063          do 1002 j=1,2
1064             ftp(j)=ftsgs(isg,j)
1065             gxp(j)=gxsgs(isg,j)
1066             gyp(j)=gysgs(isg,j)
1067             gzp(j)=gzsgs(isg,j)
1068             pxp(j)=pxsgs(isg,j)
1069             pyp(j)=pysgs(isg,j)
1070             pzp(j)=pzsgs(isg,j)
1071             pmp(j)=pmsgs(isg,j)
1072             pep(j)=pesgs(isg,j)
1073  1002    continue
1074          call locldr(2,drlocl)
1075          dr0=drlocl
1076 c     dp0^2 defined as (p1+p2)^2-(m1+m2)^2:
1077          dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1078      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1079 c
1080          do 120 JSG=1,NSG
1081 c     skip default or unavailable antiquarks:
1082             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 120
1083             if(NJSGS(JSG).eq.2) then
1084                ipmin=2
1085                ipmax=2
1086             elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).lt.0) then
1087                ipmin=1
1088                ipmax=3
1089             else
1090                goto 120
1091             endif
1092             do 100 ip=ipmin,ipmax
1093                dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1094      1              -pxp(1)*pxsgs(jsg,ip)
1095      2              -pyp(1)*pysgs(jsg,ip)
1096      3              -pzp(1)*pzsgs(jsg,ip)
1097      4              -pmp(1)*pmsgs(jsg,ip)))
1098 c     skip if outside of momentum radius:
1099                if(dplocl.gt.dpcoal) goto 120
1100                ftp(2)=ftsgs(jsg,ip)
1101                gxp(2)=gxsgs(jsg,ip)
1102                gyp(2)=gysgs(jsg,ip)
1103                gzp(2)=gzsgs(jsg,ip)
1104                pxp(2)=pxsgs(jsg,ip)
1105                pyp(2)=pysgs(jsg,ip)
1106                pzp(2)=pzsgs(jsg,ip)
1107                pmp(2)=pmsgs(jsg,ip)
1108                pep(2)=pesgs(jsg,ip)
1109                call locldr(2,drlocl)
1110 c     skip if outside of spatial radius:
1111                if(drlocl.gt.drcoal) goto 120
1112 c     q_isg coalesces with qbar_jsg:
1113                if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1114      1              .or.(drlocl.lt.dr0)) then
1115                   dp0=dplocl
1116                   dr0=drlocl
1117                   call exchge(isg,2,jsg,ip)
1118                endif
1119  100        continue
1120  120     continue
1121          if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1122  150  continue
1123 c
1124 C2     meson qbar coalesce with all available q:
1125       do 250 ISG=1,NSG
1126          if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 250
1127 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1128          do 1003 j=1,2
1129             ftp(j)=ftsgs(isg,j)
1130             gxp(j)=gxsgs(isg,j)
1131             gyp(j)=gysgs(isg,j)
1132             gzp(j)=gzsgs(isg,j)
1133             pxp(j)=pxsgs(isg,j)
1134             pyp(j)=pysgs(isg,j)
1135             pzp(j)=pzsgs(isg,j)
1136             pmp(j)=pmsgs(isg,j)
1137             pep(j)=pesgs(isg,j)
1138  1003    continue
1139          call locldr(2,drlocl)
1140          dr0=drlocl
1141          dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1142      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1143 c
1144          do 220 JSG=1,NSG
1145             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 220
1146             if(NJSGS(JSG).eq.2) then
1147                ipmin=1
1148                ipmax=1
1149             elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).gt.0) then
1150                ipmin=1
1151                ipmax=3
1152             else
1153                goto 220
1154             endif
1155             do 200 ip=ipmin,ipmax
1156                dplocl=dsqrt(2*(pep(2)*pesgs(jsg,ip)
1157      1              -pxp(2)*pxsgs(jsg,ip)
1158      2              -pyp(2)*pysgs(jsg,ip)
1159      3              -pzp(2)*pzsgs(jsg,ip)
1160      4              -pmp(2)*pmsgs(jsg,ip)))
1161 c     skip if outside of momentum radius:
1162                if(dplocl.gt.dpcoal) goto 220
1163                ftp(1)=ftsgs(jsg,ip)
1164                gxp(1)=gxsgs(jsg,ip)
1165                gyp(1)=gysgs(jsg,ip)
1166                gzp(1)=gzsgs(jsg,ip)
1167                pxp(1)=pxsgs(jsg,ip)
1168                pyp(1)=pysgs(jsg,ip)
1169                pzp(1)=pzsgs(jsg,ip)
1170                pmp(1)=pmsgs(jsg,ip)
1171                pep(1)=pesgs(jsg,ip)
1172                call locldr(2,drlocl)
1173 c     skip if outside of spatial radius:
1174                if(drlocl.gt.drcoal) goto 220
1175 c     qbar_isg coalesces with q_jsg:
1176                if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1177      1              .or.(drlocl.lt.dr0)) then
1178                   dp0=dplocl
1179                   dr0=drlocl
1180                   call exchge(isg,1,jsg,ip)
1181                endif
1182  200        continue
1183  220     continue
1184          if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1185  250  continue
1186 c
1187 C3     baryon q (antibaryon qbar) coalesce with all available q (qbar):
1188       do 350 ISG=1,NSG
1189          if(NJSGS(ISG).ne.3.or.IOVER(ISG).eq.1) goto 350
1190          ibaryn=K2SGS(ISG,1)
1191 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1192          do 1004 j=1,2
1193             ftp(j)=ftsgs(isg,j)
1194             gxp(j)=gxsgs(isg,j)
1195             gyp(j)=gysgs(isg,j)
1196             gzp(j)=gzsgs(isg,j)
1197             pxp(j)=pxsgs(isg,j)
1198             pyp(j)=pysgs(isg,j)
1199             pzp(j)=pzsgs(isg,j)
1200             pmp(j)=pmsgs(isg,j)
1201             pep(j)=pesgs(isg,j)
1202  1004    continue
1203          call locldr(2,drlocl)
1204          dr1(2)=drlocl
1205          dp1(2)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1206      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1207 c
1208          ftp(2)=ftsgs(isg,3)
1209          gxp(2)=gxsgs(isg,3)
1210          gyp(2)=gysgs(isg,3)
1211          gzp(2)=gzsgs(isg,3)
1212          pxp(2)=pxsgs(isg,3)
1213          pyp(2)=pysgs(isg,3)
1214          pzp(2)=pzsgs(isg,3)
1215          pmp(2)=pmsgs(isg,3)
1216          pep(2)=pesgs(isg,3)
1217          call locldr(2,drlocl)
1218          dr1(3)=drlocl
1219          dp1(3)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1220      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1221 c
1222          do 320 JSG=1,NSG
1223             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 320
1224             if(NJSGS(JSG).eq.2) then
1225                if(ibaryn.gt.0) then
1226                   ipmin=1
1227                else
1228                   ipmin=2
1229                endif
1230                ipmax=ipmin
1231             elseif(NJSGS(JSG).eq.3.and.
1232      1              (ibaryn*K2SGS(JSG,1)).gt.0) then
1233                ipmin=1
1234                ipmax=3
1235             else
1236                goto 320
1237             endif
1238             do 300 ip=ipmin,ipmax
1239                dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1240      1              -pxp(1)*pxsgs(jsg,ip)
1241      2              -pyp(1)*pysgs(jsg,ip)
1242      3              -pzp(1)*pzsgs(jsg,ip)
1243      4              -pmp(1)*pmsgs(jsg,ip)))
1244 c     skip if outside of momentum radius:
1245                if(dplocl.gt.dpcoal) goto 320
1246                ftp(2)=ftsgs(jsg,ip)
1247                gxp(2)=gxsgs(jsg,ip)
1248                gyp(2)=gysgs(jsg,ip)
1249                gzp(2)=gzsgs(jsg,ip)
1250                pxp(2)=pxsgs(jsg,ip)
1251                pyp(2)=pysgs(jsg,ip)
1252                pzp(2)=pzsgs(jsg,ip)
1253                pmp(2)=pmsgs(jsg,ip)
1254                pep(2)=pesgs(jsg,ip)
1255                call locldr(2,drlocl)
1256 c     skip if outside of spatial radius:
1257                if(drlocl.gt.drcoal) goto 320
1258 c     q_isg may coalesce with q_jsg for a baryon:
1259                ipi=0
1260                if(dp1(2).gt.dpcoal.or.dr1(2).gt.drcoal) then
1261                   ipi=2
1262                   if((dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal)
1263      1                 .and.dr1(3).gt.dr1(2)) ipi=3
1264                elseif(dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal) then
1265                   ipi=3
1266                elseif(dr1(2).lt.dr1(3)) then
1267                   if(drlocl.lt.dr1(3)) ipi=3
1268                elseif(dr1(3).le.dr1(2)) then
1269                   if(drlocl.lt.dr1(2)) ipi=2
1270                endif
1271                if(ipi.ne.0) then
1272                   dp1(ipi)=dplocl
1273                   dr1(ipi)=drlocl
1274                   call exchge(isg,ipi,jsg,ip)
1275                endif
1276  300        continue
1277  320     continue
1278          if(dp1(2).le.dpcoal.and.dr1(2).le.drcoal
1279      1        .and.dp1(3).le.dpcoal.and.dr1(3).le.drcoal)
1280      2        IOVER(ISG)=1
1281  350  continue
1282 c
1283       RETURN
1284       END
1285
1286 c=======================================================================
1287       SUBROUTINE exchge(isg,ipi,jsg,ipj)
1288 c
1289       implicit double precision  (a-h, o-z)
1290       PARAMETER (MAXSTR=150001)
1291       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1292      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1293      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1294      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1295 cc      SAVE /SOFT/
1296       SAVE
1297 c
1298       k1=K1SGS(isg,ipi)
1299       k2=K2SGS(isg,ipi)
1300       px=PXSGS(isg,ipi)
1301       py=PYSGS(isg,ipi)
1302       pz=PZSGS(isg,ipi)
1303       pe=PESGS(isg,ipi)
1304       pm=PMSGS(isg,ipi)
1305       gx=GXSGS(isg,ipi)
1306       gy=GYSGS(isg,ipi)
1307       gz=GZSGS(isg,ipi)
1308       ft=FTSGS(isg,ipi)
1309       K1SGS(isg,ipi)=K1SGS(jsg,ipj)
1310       K2SGS(isg,ipi)=K2SGS(jsg,ipj)
1311       PXSGS(isg,ipi)=PXSGS(jsg,ipj)
1312       PYSGS(isg,ipi)=PYSGS(jsg,ipj)
1313       PZSGS(isg,ipi)=PZSGS(jsg,ipj)
1314       PESGS(isg,ipi)=PESGS(jsg,ipj)
1315       PMSGS(isg,ipi)=PMSGS(jsg,ipj)
1316       GXSGS(isg,ipi)=GXSGS(jsg,ipj)
1317       GYSGS(isg,ipi)=GYSGS(jsg,ipj)
1318       GZSGS(isg,ipi)=GZSGS(jsg,ipj)
1319       FTSGS(isg,ipi)=FTSGS(jsg,ipj)
1320       K1SGS(jsg,ipj)=k1
1321       K2SGS(jsg,ipj)=k2
1322       PXSGS(jsg,ipj)=px
1323       PYSGS(jsg,ipj)=py
1324       PZSGS(jsg,ipj)=pz
1325       PESGS(jsg,ipj)=pe
1326       PMSGS(jsg,ipj)=pm
1327       GXSGS(jsg,ipj)=gx
1328       GYSGS(jsg,ipj)=gy
1329       GZSGS(jsg,ipj)=gz
1330       FTSGS(jsg,ipj)=ft
1331 c
1332       RETURN
1333       END
1334
1335 c=======================================================================
1336       SUBROUTINE locldr(icall,drlocl)
1337 c
1338       implicit double precision (a-h, o-z)
1339       dimension ftp0(3),pxp0(3),pyp0(3),pzp0(3),pep0(3)
1340       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1341      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1342 cc      SAVE /loclco/
1343       common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
1344 cc      SAVE /prtn23/
1345       common /lora/ enenew, pxnew, pynew, pznew
1346 cc      SAVE /lora/
1347       SAVE
1348 c     for 2-body kinematics:
1349       if(icall.eq.2) then
1350          etot=pep(1)+pep(2)
1351          bex=(pxp(1)+pxp(2))/etot
1352          bey=(pyp(1)+pyp(2))/etot
1353          bez=(pzp(1)+pzp(2))/etot
1354 c     boost the reference frame down by beta to get to the pair rest frame:
1355          do 1001 j=1,2
1356             beta2 = bex ** 2 + bey ** 2 + bez ** 2
1357             gam = 1.d0 / dsqrt(1.d0 - beta2)
1358             if(beta2.ge.0.9999999999999d0) then
1359                write(6,*) '4',pxp(1),pxp(2),pyp(1),pyp(2),
1360      1              pzp(1),pzp(2),pep(1),pep(2),pmp(1),pmp(2),
1361      2          dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2+pmp(1)**2)/pep(1),
1362      3          dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2)/pep(1)
1363                write(6,*) '4a',pxp(1)+pxp(2),pyp(1)+pyp(2),
1364      1              pzp(1)+pzp(2),etot
1365                write(6,*) '4b',bex,bey,bez,beta2,gam
1366             endif
1367 c
1368             call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
1369             gxp0(j)=pxnew
1370             gyp0(j)=pynew
1371             gzp0(j)=pznew
1372             ftp0(j)=enenew
1373             call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
1374             pxp0(j)=pxnew
1375             pyp0(j)=pynew
1376             pzp0(j)=pznew
1377             pep0(j)=enenew
1378  1001    continue
1379 c
1380          if(ftp0(1).ge.ftp0(2)) then
1381             ilate=1
1382             iearly=2
1383          else
1384             ilate=2
1385             iearly=1
1386          endif
1387          ft0fom=ftp0(ilate)
1388 c
1389          dt0=ftp0(ilate)-ftp0(iearly)
1390          gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1391          gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1392          gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1393          drlocl=dsqrt((gxp0(ilate)-gxp0(iearly))**2
1394      1        +(gyp0(ilate)-gyp0(iearly))**2
1395      2        +(gzp0(ilate)-gzp0(iearly))**2)
1396 c     for 3-body kinematics, used for baryons formation:
1397       elseif(icall.eq.3) then
1398          etot=pep(1)+pep(2)+pep(3)
1399          bex=(pxp(1)+pxp(2)+pxp(3))/etot
1400          bey=(pyp(1)+pyp(2)+pyp(3))/etot
1401          bez=(pzp(1)+pzp(2)+pzp(3))/etot
1402          beta2 = bex ** 2 + bey ** 2 + bez ** 2
1403          gam = 1.d0 / dsqrt(1.d0 - beta2)
1404          if(beta2.ge.0.9999999999999d0) then
1405             write(6,*) '5',bex,bey,bez,beta2,gam
1406          endif
1407 c     boost the reference frame down by beta to get to the 3-parton rest frame:
1408          do 1002 j=1,3
1409             call lorenza(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
1410             gxp0(j)=pxnew
1411             gyp0(j)=pynew
1412             gzp0(j)=pznew
1413             ftp0(j)=enenew
1414             call lorenza(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
1415             pxp0(j)=pxnew
1416             pyp0(j)=pynew
1417             pzp0(j)=pznew
1418             pep0(j)=enenew
1419  1002    continue
1420 c
1421          if(ftp0(1).gt.ftp0(2)) then
1422             ilate=1
1423             if(ftp0(3).gt.ftp0(1)) ilate=3
1424          else
1425             ilate=2
1426             if(ftp0(3).ge.ftp0(2)) ilate=3
1427          endif
1428          ft0fom=ftp0(ilate)
1429 c
1430          if(ilate.eq.1) then
1431             imin=2
1432             imax=3
1433             istep=1
1434          elseif(ilate.eq.2) then
1435             imin=1
1436             imax=3
1437             istep=2
1438          elseif(ilate.eq.3) then
1439             imin=1
1440             imax=2
1441             istep=1
1442          endif
1443 c
1444          do 1003 iearly=imin,imax,istep
1445             dt0=ftp0(ilate)-ftp0(iearly)
1446             gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1447             gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1448             gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1449  1003    continue
1450       endif
1451 c
1452       RETURN
1453       END
1454
1455 c=======================================================================
1456         subroutine hoscar
1457 c
1458         parameter (MAXSTR=150001,AMN=0.939457,AMP=0.93828)
1459         character*8 code, reffra, FRAME
1460         character*25 amptvn
1461         common/snn/efrm,npart1,npart2
1462 cc      SAVE /snn/
1463         common /lastt/itimeh,bimp
1464 cc      SAVE /lastt/
1465         COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
1466 cc      SAVE /hbt/
1467         common/oscar1/iap,izp,iat,izt
1468 cc      SAVE /oscar1/
1469         common/oscar2/FRAME,amptvn
1470 cc      SAVE /oscar2/
1471         SAVE
1472         data nff/0/
1473 c
1474 c       file header
1475         if(nff.eq.0) then
1476            write (19, 101) 'OSCAR1997A'
1477            write (19, 111) 'final_id_p_x'
1478            code = 'AMPT'
1479            if(FRAME.eq.'CMS') then
1480               reffra = 'nncm'
1481               xmp=(amp*izp+amn*(iap-izp))/iap
1482               xmt=(amp*izt+amn*(iat-izt))/iat
1483               ebeam=(efrm**2-xmp**2-xmt**2)/2./xmt
1484            elseif(FRAME.eq.'LAB') then
1485               reffra = 'lab'
1486               ebeam=efrm
1487            else
1488               reffra = 'unknown'
1489               ebeam=0.
1490            endif
1491            ntestp = 1
1492            write (19, 102) code, amptvn, iap, izp, iat, izt,
1493      &        reffra, ebeam, ntestp
1494            nff = 1
1495            ievent = 1
1496            phi = 0.
1497            if(FRAME.eq.'CMS') write(19,112) efrm
1498         endif
1499 c       comment
1500 c       event header
1501         write (19, 103) ievent, nlast, bimp, phi
1502 c       particles
1503         do 99 i = 1, nlast
1504            ene=sqrt(plast(1,i)**2+plast(2,i)**2+plast(3,i)**2
1505      1          +plast(4,i)**2)
1506            write (19, 104) i, INVFLV(lblast(i)), plast(1,i),
1507      1          plast(2,i),plast(3,i),ene,plast(4,i),
1508      2          xlast(1,i),xlast(2,i),xlast(3,i),xlast(4,i)
1509  99     continue
1510         ievent = ievent + 1
1511  101        format (a10)
1512  111        format (a12)
1513  102        format (a4,1x,a20,1x,'(', i3, ',', i3, ')+(', i3, ',',
1514      &           i3, ')', 2x, a4, 2x, e10.4, 2x, i8)
1515  103        format (i10, 2x, i10, 2x, f8.3, 2x, f8.3)
1516  104        format (i10, 2x, i10, 2x, 9(e12.6, 2x))
1517  112        format ('# Center-of-mass energy/nucleon-pair is',
1518      & f12.3,'GeV')
1519 c
1520         return
1521         end
1522
1523 c=======================================================================
1524         subroutine getnp
1525
1526         PARAMETER (MAXSTR=150001)
1527         COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
1528 cc      SAVE /HMAIN1/
1529         COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
1530 cc      SAVE /HMAIN2/
1531         COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50)
1532 cc      SAVE /HPARNT/
1533         common/snn/efrm,npart1,npart2
1534 cc      SAVE /snn/
1535 cwei        DOUBLE PRECISION PATT
1536         COMMON/EXTRACENT/BB
1537         SAVE
1538
1539         BB = HINT1(19)
1540         if(NATT.eq.0) then
1541            npart1=0
1542            npart2=0
1543            return
1544         endif
1545 c
1546         PZPROJ=SQRT(HINT1(6)**2-HINT1(8)**2)
1547         PZTARG=SQRT(HINT1(7)**2-HINT1(9)**2)
1548         epsiln=0.01
1549 c
1550         nspec1=0
1551         nspec2=0
1552         DO 1000 I = 1, NATT
1553            if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212)
1554      1          .and.PATT(I, 1).eq.0.and.PATT(I, 2).eq.0) then
1555               if(PATT(I, 3).gt.amax1(0.,PZPROJ-epsiln)) then
1556                  nspec1=nspec1+1
1557               elseif(PATT(I, 3).lt.(-PZTARG+epsiln)) then
1558                  nspec2=nspec2+1
1559               endif
1560            endif
1561  1000    CONTINUE
1562         npart1=IHNT2(1)-nspec1
1563         npart2=IHNT2(3)-nspec2
1564
1565         return
1566         end
1567
1568 c=======================================================================
1569 c     2/18/03 use PYTHIA to decay eta,rho,omega,k*,phi and Delta
1570         subroutine resdec(i1,nt,nnn,wid,idecay)
1571
1572         PARAMETER (hbarc=0.19733)
1573         PARAMETER (AK0=0.498,APICH=0.140,API0=0.135,AN=0.940,ADDM=0.02)
1574         PARAMETER (MAXSTR=150001, MAXR=1)
1575         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1576      &       IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1577 cc      SAVE /INPUT2/
1578         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1579 cc      SAVE /LUJETSA/
1580         COMMON/LUDAT1A/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1581 cc      SAVE /LUDAT1A/
1582         COMMON/LUDAT2A/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
1583 cc      SAVE /LUDAT2A/
1584         COMMON/LUDAT3A/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
1585 cc      SAVE /LUDAT3A/
1586         COMMON /CC/ E(MAXSTR)
1587 cc      SAVE /CC/
1588         COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1589 cc      SAVE /EE/
1590         COMMON   /PA/RPION(3,MAXSTR,MAXR)
1591 cc      SAVE /PA/
1592         COMMON   /PB/PPION(3,MAXSTR,MAXR)
1593 cc      SAVE /PB/
1594         COMMON   /PC/EPION(MAXSTR,MAXR)
1595 cc      SAVE /PC/
1596         COMMON   /PD/LPION(MAXSTR,MAXR)
1597 cc      SAVE /PD/
1598         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1599 cc      SAVE /input1/
1600         common/resdcy/NSAV,iksdcy
1601 cc      SAVE /resdcy/
1602         common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
1603      1       px1n,py1n,pz1n,dp1n
1604 cc      SAVE /leadng/
1605         EXTERNAL IARFLV, INVFLV
1606         COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
1607 cc      SAVE /tdecay/
1608         COMMON/RNDF77/NSEED
1609         COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1610      1       dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1611      2       dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1612 cc      SAVE /RNDF77/
1613         SAVE
1614         irun=idecay
1615         if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1616      &       .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1617 c     &       .or.lb1.eq.28.or.lb1.eq.29
1618      &       .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1619      &       .or.iabs(lb1).eq.16) then
1620            kf=INVFLV(lb1)
1621         else
1622            return
1623         endif
1624 c
1625         IP=1
1626 c     label as undecayed and the only particle in the record:
1627         N=1
1628         K(IP,1)=1
1629         K(IP,3)=0
1630         K(IP,4)=0
1631         K(IP,5)=0
1632 c
1633         K(IP,2)=kf
1634         P(IP,1)=px1
1635         P(IP,2)=py1
1636         P(IP,3)=pz1
1637         em1a=em1
1638 c     eta or omega in ART may be below or too close to (pi+pi-pi0) mass,
1639 c     causing LUDECY error,thus increase their mass ADDM above this thresh,
1640 c     noting that rho (m=0.281) too close to 2pi thrshold fails to decay:
1641         if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then
1642            em1=2*APICH+API0+ADDM
1643 c     rho
1644         elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then
1645            em1=2*APICH+ADDM
1646 c     K*
1647         elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then
1648            em1=APICH+AK0+ADDM
1649 c     Delta created in ART may be below (n+pich) mass, causing LUDECY error:
1650         elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9
1651      1          .and.em1.lt.(APICH+AN+ADDM)) then
1652            em1=APICH+AN+ADDM
1653         endif
1654 c        if(em1.ge.(em1a+0.01)) write (6,*)
1655 c     1       'Mass increase in resdec():',nt,em1-em1a,lb1
1656         e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1657         P(IP,4)=e1
1658         P(IP,5)=em1
1659 clin-5/2008:
1660         dpdecp=dpertp(i1)
1661         call ludecy(IP)
1662 c     add decay time to daughter's formation time at the last timestep:
1663         if(nt.eq.ntmax) then
1664            tau0=hbarc/wid
1665            taudcy=tau0*(-1.)*alog(1.-RANART(NSEED))
1666            ndaut=n-nsav
1667            if(ndaut.le.1) then
1668               write(10,*) 'note: ndaut(<1)=',ndaut
1669               write(89,*) 'note: ndaut(<1)=',ndaut
1670               call lulist(2)
1671               stop
1672             endif
1673 c     lorentz boost:
1674            taudcy=taudcy*e1/em1
1675            tfnl=tfnl+taudcy
1676            xfnl=xfnl+px1/e1*taudcy
1677            yfnl=yfnl+py1/e1*taudcy
1678            zfnl=zfnl+pz1/e1*taudcy
1679 c     at the last timestep, assign rho, K0S or eta (decay daughter)
1680 c     to lb(i1) only (not to lpion) in order to decay them again:
1681            if(n.ge.(nsav+2)) then
1682               do 1001 idau=nsav+2,n
1683                  kdaut=K(idau,2)
1684                  if(kdaut.eq.221.or.kdaut.eq.113
1685      1                .or.kdaut.eq.213.or.kdaut.eq.-213
1686      2                .or.kdaut.eq.310) then
1687 c     switch idau and i1(nsav+1):
1688                     ksave=kdaut
1689                     pxsave=p(idau,1)
1690                     pysave=p(idau,2)
1691                     pzsave=p(idau,3)
1692                     esave=p(idau,4)
1693                     xmsave=p(idau,5)
1694                     K(idau,2)=K(nsav+1,2)
1695                     p(idau,1)=p(nsav+1,1)
1696                     p(idau,2)=p(nsav+1,2)
1697                     p(idau,3)=p(nsav+1,3)
1698                     p(idau,4)=p(nsav+1,4)
1699                     p(idau,5)=p(nsav+1,5)
1700                     K(nsav+1,2)=ksave
1701                     p(nsav+1,1)=pxsave
1702                     p(nsav+1,2)=pysave
1703                     p(nsav+1,3)=pzsave
1704                     p(nsav+1,4)=esave
1705                     p(nsav+1,5)=xmsave
1706 c     note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce
1707 c     eta, but only one daughter may be rho, K0s or eta:
1708                     goto 111
1709                  endif
1710  1001         continue
1711            endif
1712  111       continue
1713 c
1714            enet=0.
1715            do 1002 idau=nsav+1,n
1716               enet=enet+p(idau,4)
1717  1002      continue
1718 c           if(abs(enet-e1).gt.0.02)
1719 c     1          write(93,*) 'resdec(): nt=',nt,enet-e1,lb1
1720         endif
1721
1722 cyy 200    format(a20,3(1x,i6))
1723 cyy 210    format(i6,5(1x,f8.3))
1724 cyy 220    format(a2,i5,5(1x,f8.3))
1725
1726         do 1003 idau=nsav+1,n
1727            kdaut=K(idau,2)
1728            lbdaut=IARFLV(kdaut)
1729 c     K0S and K0L are named K+/K- during hadron cascade, and only
1730 c     at the last timestep they keep their real LB # before output;
1731 c     K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep:
1732            if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310
1733      1          .or.iabs(kdaut).eq.311)) then
1734               if(kdaut.eq.130) then
1735                  lbdaut=22
1736               elseif(kdaut.eq.310) then
1737                  lbdaut=24
1738               elseif(iabs(kdaut).eq.311) then
1739                  if(RANART(NSEED).lt.0.5) then
1740                     lbdaut=22
1741                  else
1742                     lbdaut=24
1743                  endif
1744               endif
1745            endif
1746 c
1747            if(idau.eq.(nsav+1)) then
1748               LB(i1)=lbdaut
1749               E(i1)=p(idau,5)
1750               px1n=p(idau,1)
1751               py1n=p(idau,2)
1752               pz1n=p(idau,3)
1753 clin-5/2008:
1754               dp1n=dpdecp
1755            else
1756               nnn=nnn+1
1757               LPION(NNN,IRUN)=lbdaut
1758               EPION(NNN,IRUN)=p(idau,5)
1759               PPION(1,NNN,IRUN)=p(idau,1)
1760               PPION(2,NNN,IRUN)=p(idau,2)
1761               PPION(3,NNN,IRUN)=p(idau,3)
1762               RPION(1,NNN,IRUN)=xfnl
1763               RPION(2,NNN,IRUN)=yfnl
1764               RPION(3,NNN,IRUN)=zfnl
1765               tfdpi(NNN,IRUN)=tfnl
1766 clin-5/2008:
1767               dppion(NNN,IRUN)=dpdecp
1768            endif
1769  1003   continue
1770 cyy 230    format(a2,i5,5(1x,e8.2))
1771         return
1772         end
1773
1774 c=======================================================================
1775         subroutine inidcy
1776
1777         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1778 cc      SAVE /LUJETSA/
1779         common/resdcy/NSAV,iksdcy
1780 cc      SAVE /resdcy/
1781         SAVE
1782         N=1
1783         NSAV=N
1784         return
1785         end
1786
1787 c=======================================================================
1788 clin-6/06/02 local parton freezeout motivated from critical density:
1789         subroutine local(t)
1790 c
1791         implicit double precision  (a-h, o-z)
1792         PARAMETER (MAXPTN=400001)
1793         PARAMETER (r0=1d0)
1794         COMMON /para1/ mul
1795 cc      SAVE /para1/
1796         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1797      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1798      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
1799 cc      SAVE /prec2/
1800         common /frzprc/
1801      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1802      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1803      &       xmfrz(MAXPTN),
1804      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1805 cc      SAVE /frzprc/
1806         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1807 cc      SAVE /prec4/
1808         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1809 cc      SAVE /prec5/
1810         common /coal/dpcoal,drcoal,ecritl
1811 cc      SAVE /coal/
1812         SAVE
1813 c
1814       do 1001 it=1,301
1815          if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then
1816             if(it.eq.itlast) then
1817                return
1818             else
1819                itlast=it
1820                goto 50
1821             endif
1822          endif
1823  1001 continue
1824       write(1,*) 'local time out of range in LOCAL, stop',t,it
1825       stop
1826  50   continue
1827 c
1828       do 200 ip=1,mul
1829 c     skip partons which have frozen out:
1830          if(ifrz(ip).eq.1) goto 200
1831          if(it.eq.301) then
1832 c     freezeout all the left partons beyond the time of 3000 fm:
1833             etcrit=1d6
1834             goto 150
1835          else
1836 c     freezeout when transverse energy density < etcrit:
1837             etcrit=(ecritl*2d0/3d0)
1838          endif
1839 c     skip partons which have not yet formed:
1840          if(t.lt.FT5(ip)) goto 200
1841          rap0=rap(ip)
1842          eta0=eta(ip)
1843          x0=GX5(ip)+vx(ip)*(t-FT5(ip))
1844          y0=GY5(ip)+vy(ip)*(t-FT5(ip))
1845          detdy=0d0
1846          do 100 itest=1,mul
1847 c     skip self and partons which have not yet formed:
1848             if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100
1849             ettest=eta(itest)
1850             xtest=GX5(itest)+vx(itest)*(t-FT5(itest))
1851             ytest=GY5(itest)+vy(itest)*(t-FT5(itest))
1852             drt=sqrt((xtest-x0)**2+(ytest-y0)**2)
1853 c     count partons within drt<1 and -1<(eta-eta0)<1:
1854             if(dabs(ettest-eta0).le.1d0.and.drt.le.r0)
1855      1           detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2
1856      2           +XMASS5(itest)**2)*0.5d0
1857  100     continue
1858          detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0))
1859 c     when density is below critical density for phase transition, freeze out:
1860  150     if(detdy.le.etcrit) then
1861             ifrz(ip)=1
1862             idfrz(ip)=ITYP5(ip)
1863             pxfrz(ip)=PX5(ip)
1864             pyfrz(ip)=PY5(ip)
1865             pzfrz(ip)=PZ5(ip)
1866             efrz(ip)=E5(ip)
1867             xmfrz(ip)=XMASS5(ip)
1868             if(t.gt.FT5(ip)) then
1869                gxfrz(ip)=x0
1870                gyfrz(ip)=y0
1871                gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip))
1872                ftfrz(ip)=t
1873             else
1874 c     if this freezeout time < formation time, use formation time & positions.
1875 c     This ensures the recovery of default hadron when e_crit=infty:
1876                gxfrz(ip)=GX5(ip)
1877                gyfrz(ip)=GY5(ip)
1878                gzfrz(ip)=GZ5(ip)
1879                ftfrz(ip)=FT5(ip)
1880             endif
1881          endif
1882  200  continue
1883 c
1884         return
1885         end
1886
1887 c=======================================================================
1888 clin-6/06/02 initialization for local parton freezeout
1889         subroutine inifrz
1890 c
1891         implicit double precision  (a-h, o-z)
1892         PARAMETER (MAXPTN=400001)
1893         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1894 cc      SAVE /ilist5/
1895         common /frzprc/
1896      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1897      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1898      &       xmfrz(MAXPTN),
1899      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1900 cc      SAVE /frzprc/
1901         SAVE
1902 c
1903 c     for freezeout time 0-10fm, use interval of 0.1fm;
1904 c     for 10-100fm, use interval of 1fm;
1905 c     for 100-1000fm, use interval of 10fm;
1906 c     for 1000-3000fm, use interval of 100fm:
1907         step1=0.1d0
1908         step2=1d0
1909         step3=10d0
1910         step4=100d0
1911 c
1912         do 1001 it=1,101
1913            tfrz(it)=0d0+dble(it-1)*step1
1914  1001 continue
1915         do 1002 it=102,191
1916            tfrz(it)=10d0+dble(it-101)*step2
1917  1002   continue
1918         do 1003 it=192,281
1919            tfrz(it)=100d0+dble(it-191)*step3
1920  1003   continue
1921         do 1004 it=282,301
1922            tfrz(it)=1000d0+dble(it-281)*step4
1923  1004   continue
1924         tfrz(302)=tlarge
1925 c
1926         return
1927         end
1928
1929 clin-5/2009 ctest on v2 analysis
1930 c=======================================================================
1931 c     idd=0,1,2,3 specifies different subroutines for partonic flow analysis.
1932       subroutine flowp(idd)
1933 c
1934         implicit double precision  (a-h, o-z)
1935         real dt
1936         parameter (MAXPTN=400001)
1937 csp
1938         parameter (bmt=0.05d0)
1939         dimension nlfile(3),nsfile(3),nmfile(3)
1940 c
1941         dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3)
1942         dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3)
1943         common /ilist1/
1944      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1945      &     ictype, icsta(MAXPTN),
1946      &     nic(MAXPTN), icels(MAXPTN)
1947 cc      SAVE /ilist1/
1948         COMMON /para1/ mul
1949 cc      SAVE /para1/
1950         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1951      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1952      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
1953 cc      SAVE /prec2/
1954         COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3),
1955      1       s2p(30,3),v2p2(30,3),nevt(30)
1956 cc      SAVE /pflow/
1957         COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3),
1958      1                 xncoll(30),s2pf(30,3),v2pf2(30,3)
1959 cc      SAVE /pflowf/
1960         COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3),
1961      1       s2pfrz(30,3),v2p2fz(30,3),tscatt(31),
1962      2       nevtfz(30),iscatt(30)
1963 cc      SAVE /pfrz/
1964         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
1965      1 v2h2(30,3),s2h(30,3)
1966 cc      SAVE /hflow/
1967         COMMON /AREVT/ IAEVT, IARUN, MISS
1968 cc      SAVE /AREVT/
1969         common/anim/nevent,isoft,isflag,izpc
1970 cc      SAVE /anim/
1971         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1972 cc      SAVE /input1/
1973         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1974      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1975 cc      SAVE /INPUT2/
1976 cc      SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm
1977 cc      SAVE nfile,itanim,nlfile,nsfile,nmfile
1978         SAVE
1979 csp
1980         dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3),
1981      &       etlfrz(30,3),etsfrz(30,3),
1982      &       xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3),
1983      &       xnlfrz(30,3),xnsfrz(30,3),
1984      &       v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3),
1985      &       s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3),
1986      &       DMYil(50,3),DMYfl(50,3),
1987      &       DMYis(50,3),DMYfs(50,3)
1988         data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1989      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
1990      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
1991 c     idd=0: initialization for flow analysis, called by artdri.f:
1992         if(idd.eq.0) then
1993            nfile(1)=60
1994            nfile(2)=64
1995            nfile(3)=20
1996 cms        OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN')
1997 cms        OPEN (nfile(1)+1,
1998 cms  1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN')
1999 cms        OPEN (nfile(1)+2,
2000 cms  1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN')
2001 cms        OPEN (nfile(1)+3,
2002 cms  1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN')
2003 cms        OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN')
2004 cms        OPEN (nfile(2)+1,
2005 cms  1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN')
2006 cms        OPEN (nfile(2)+2,
2007 cms  1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN')
2008 cms        OPEN (nfile(2)+3,
2009 cms  1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN')
2010 cms        OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN')
2011 cms        OPEN (nfile(3)+1,
2012 cms  1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN')
2013 cms        OPEN (nfile(3)+2,
2014 cms  1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN')
2015 cms        OPEN (nfile(3)+3,
2016 cms  1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN')
2017 cms        OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN')
2018 cms        write(49, *) '    ievt,  v2p,  v2p_y2,   v2p_y1'
2019 c
2020 cms        OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN')
2021 cms        OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN')
2022 cms        OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN')
2023 cms        OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN')
2024 cms        write(88, *) '    ievt,  v2h,  v2h_y2,   v2h_y1'
2025 csp07/05
2026            nlfile(1)=70
2027            nlfile(2)=72
2028            nlfile(3)=74
2029 cms        OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN')
2030 cms        OPEN (nlfile(1)+1,
2031 cms  1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN')
2032 cms        OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN')
2033 cms        OPEN (nlfile(2)+1,
2034 cms  1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN')
2035 cms        OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN')
2036 cms        OPEN (nlfile(3)+1,
2037 cms  1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN')
2038            nsfile(1)=76
2039            nsfile(2)=78
2040            nsfile(3)=80
2041 cms        OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN')
2042 cms        OPEN (nsfile(1)+1,
2043 cms  1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN')
2044 cms        OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN')
2045 cms        OPEN (nsfile(2)+1,
2046 cms  1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN')
2047 cms        OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN')
2048 cms        OPEN (nsfile(3)+1,
2049 cms  1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN')
2050            nmfile(1)=82
2051            nmfile(2)=83
2052            nmfile(3)=84
2053 cms        OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN')
2054 cms        OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN')
2055 cms        OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN')
2056 clin-11/27/00 for animation:
2057            if(nevent.eq.1) then
2058 cms        OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN')
2059 cms        write(91,*) ntmax, dt
2060 cms        OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN')
2061 cms        OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN')
2062            endif
2063 c
2064            itimep=0
2065            itanim=0
2066            iaevtp=0
2067 csp
2068            do 1002 ii=1,50
2069               do 1001 iy=1,3
2070                  DMYil(ii,iy) = 0d0
2071                  DMYfl(ii,iy) = 0d0
2072                  DMYis(ii,iy) = 0d0
2073                  DMYfs(ii,iy) = 0d0
2074  1001         continue
2075  1002      continue
2076 c
2077            do 1003 ii=1,31
2078               tscatt(ii)=0d0
2079  1003      continue
2080            do 1005 ii=1,30
2081               nevt(ii)=0
2082               xncoll(ii)=0d0
2083               nevtfz(ii)=0
2084               iscatt(ii)=0
2085               do 1004 iy=1,3
2086                  v2p(ii,iy)=0d0
2087                  v2p2(ii,iy)=0d0
2088                  s2p(ii,iy)=0d0
2089                  etp(ii,iy)=0d0
2090                  xnpart(ii,iy)=0d0
2091                  v2pf(ii,iy)=0d0
2092                  v2pf2(ii,iy)=0d0
2093                  s2pf(ii,iy)=0d0
2094                  etpf(ii,iy)=0d0
2095                  xnpf(ii,iy)=0d0
2096                  v2pfrz(ii,iy)=0d0
2097                  v2p2fz(ii,iy)=0d0
2098                  s2pfrz(ii,iy)=0d0
2099                  etpfrz(ii,iy)=0d0
2100                  xnpfrz(ii,iy)=0d0
2101 csp07/05
2102                  etpl(ii,iy)=0d0
2103                  etps(ii,iy)=0d0
2104                  etplf(ii,iy)=0d0
2105                  etpsf(ii,iy)=0d0
2106                  etlfrz(ii,iy)=0d0
2107                  etsfrz(ii,iy)=0d0
2108               xnpl(ii,iy)=0d0
2109               xnps(ii,iy)=0d0
2110               xnplf(ii,iy)=0d0
2111               xnpsf(ii,iy)=0d0
2112               xnlfrz(ii,iy)=0d0
2113               xnsfrz(ii,iy)=0d0
2114               v2pl(ii,iy)=0d0
2115               v2ps(ii,iy)=0d0
2116               v2plf(ii,iy)=0d0
2117               v2psf(ii,iy)=0d0
2118               s2pl(ii,iy)=0d0
2119               s2ps(ii,iy)=0d0
2120               s2plf(ii,iy)=0d0
2121               s2psf(ii,iy)=0d0
2122  1004      continue
2123  1005   continue
2124            do 1006 iy=1,3
2125               v2pevt(iy)=0d0
2126               v2pavg(iy)=0d0
2127               varv2p(iy)=0d0
2128               v2pp(iy)=0.d0
2129               xnpp(iy)=0d0
2130               v2psum(iy)=0.d0
2131               v2p2sm(iy)=0.d0
2132  1006      continue
2133 c     idd=1: calculate parton elliptic flow, called by zpc.f:
2134         else if(idd.eq.1) then
2135            t2time = FT5(iscat)
2136            do 1008 ianp = 1, 30
2137               if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then
2138 c     write flow info only once at each fixed time:
2139                  xncoll(ianp)=xncoll(ianp)+1d0
2140 c     to prevent an earlier t2time comes later in the same event
2141 c     and mess up nevt:
2142                  if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101
2143                  nevt(ianp)=nevt(ianp)+1
2144                  tscatt(ianp+1)=t2time
2145                  iscatt(ianp)=1
2146                  nevtfz(ianp)=nevtfz(ianp)+1
2147                  do 100 I=1,mul
2148             ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8))
2149                     pt2=PX5(I)**2+PY5(I)**2
2150 ctest off: initial (pt,y) and (x,y) distribution:
2151 c                    idtime=1
2152 c                    if(ianp.eq.idtime) then
2153 c                       iityp=iabs(ITYP5(I))
2154 c                       if(iityp.eq.1.or.iityp.eq.2) then
2155 c                          write(651,*) dsqrt(pt2),ypartn
2156 c                          write(654,*) GX5(I),GY5(I)
2157 c                       elseif(iityp.eq.1103.or.iityp.eq.2101
2158 c     1 .or.iityp.eq.2103.or.iityp.eq.2203.
2159 c     2 .or.iityp.eq.3101.or.iityp.eq.3103.
2160 c     3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303)
2161 c     4 then
2162 c                          write(652,*) dsqrt(pt2),ypartn
2163 c                          write(655,*) GX5(I),GY5(I)
2164 c                       elseif(iityp.eq.21) then
2165 c                          write(653,*) dsqrt(pt2),ypartn
2166 c                          write(656,*) GX5(I),GY5(I)
2167 c                       endif
2168 c                    endif
2169 ctest-end
2170 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2171 c                    gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2172 c                    gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2173 c                    gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2174 c                    rt_now=dsqrt(gx_now**2+gy_now**2)
2175 c                    zmax=0.1d0*t2time
2176 c                    volume=3.1416d0*(2d0**2)*(2*zmax)
2177 c                    if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax)
2178 c     1                   goto 100
2179 ctest-end
2180                     iloop=1
2181                     if(dabs(ypartn).le.1d0) then
2182                        iloop=2
2183                        if(dabs(ypartn).le.0.5d0) then
2184                           iloop=3
2185                        endif
2186                     endif
2187                     do 50 iy=1,iloop
2188                        if(pt2.gt.0.) then
2189                           v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2190                           if(abs(v2prtn).gt.1.)
2191      1 write(nfile(iy),*) 'v2prtn>1',v2prtn
2192                           v2p(ianp,iy)=v2p(ianp,iy)+v2prtn
2193                           v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2
2194                        endif
2195                        xperp2=GX5(I)**2+GY5(I)**2
2196                        if(xperp2.gt.0.)
2197      1        s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2198                        xnpart(ianp,iy)=xnpart(ianp,iy)+1d0
2199                        etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2200 ctest off density:
2201 c                       etp(ianp,iy)=etp(ianp,iy)
2202 c     1                  +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2203 clin-2/22/00 to write out parton info only for formed ones:
2204                        if(FT5(I).le.t2time) then
2205                           v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn
2206                           v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2
2207                           if(xperp2.gt.0.)
2208      1        s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2209                           xnpf(ianp,iy)=xnpf(ianp,iy)+1d0
2210                   etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2211 ctest off density:
2212 c                  etpf(ianp,iy)=etpf(ianp,iy)
2213 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2214                        endif
2215  50                    continue
2216  100                 continue
2217                  itimep=ianp
2218                  iaevtp=iaevt
2219 clin-3/30/00 ebe v2 variables:
2220                  if(ianp.eq.30) then
2221                     do 1007 iy=1,3
2222                        npartn=IDINT(xnpart(ianp,iy)-xnpp(iy))
2223                        if(npartn.ne.0) then
2224                           v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn
2225                           v2psum(iy)=v2psum(iy)+v2pevt(iy)
2226                           v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2
2227                           v2pp(iy)=v2p(ianp,iy)
2228                           xnpp(iy)=xnpart(ianp,iy)
2229                        endif
2230  1007               continue
2231                     write(49, 160) iaevt,v2pevt
2232                  endif
2233                  goto 101
2234               endif
2235  1008      continue
2236 clin-11/28/00 for animation:
2237  101           if(nevent.eq.1) then
2238               do 110 nt = 1, ntmax
2239                  time1=dble(nt*dt)
2240                  time2=dble((nt+1)*dt)
2241                  if (t2time.lt.time2.and.t2time.ge.time1) then
2242                     if(nt.le.itanim) return
2243 c                    write(92,*) t2time
2244                     iform=0
2245                     do 1009 I=1,mul
2246 c     write out parton info only for formed ones:
2247                        if(FT5(I).le.t2time) then
2248                           iform=iform+1
2249                        endif
2250  1009               continue
2251 c                    write(92,*) iform
2252                     do 120 I=1,mul
2253                        if(FT5(I).le.t2time) then
2254 clin-11/29/00-ctest off calculate parton coordinates after propagation:
2255 c                          gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2256 c                          gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2257 c                          gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2258 c          write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2259 c          write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i),
2260 c     1    PX5(i),PY5(i),PZ5(i),E5(i)
2261 ctest-end
2262                        endif
2263  120                    continue
2264                     itanim=nt
2265                  endif
2266  110              continue
2267            endif
2268 c
2269  140           format(i10,4(2x,f7.2))
2270  160           format(i10,3(2x,f9.5))
2271 c 180           format(i6,8(1x,f7.2))
2272 clin-5/17/01 calculate v2 for active partons (which have not frozen out):
2273 c     idd=3, called at end of zpc.f:
2274         else if(idd.eq.3) then
2275            do 1010 ianp=1,30
2276               if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp)
2277  1010      continue
2278            do 350 I=1,mul
2279               ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8)
2280      1 /(E5(i)-PZ5(i)+1.d-8))
2281               pt2=PX5(I)**2+PY5(I)**2
2282               iloop=1
2283               if(dabs(ypartn).le.1d0) then
2284                  iloop=2
2285                  if(dabs(ypartn).le.0.5d0) then
2286                     iloop=3
2287                  endif
2288               endif
2289 c
2290               do 325 ianp=1,30
2291                  if(iscatt(ianp).ne.0) then
2292                     if(FT5(I).lt.tscatt(ianp+1)
2293      1 .and.FT5(I).ge.tscatt(ianp)) then
2294                        do 1011 iy=1,iloop
2295                           if(pt2.gt.0.) then
2296                              v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2297                              v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn
2298                      v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2
2299                           endif
2300                           xperp2=GX5(I)**2+GY5(I)**2
2301                           if(xperp2.gt.0.) s2pfrz(ianp,iy)=
2302      1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2303         etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2304                           xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0
2305 ctest off density:
2306 c                    etpfrz(ianp,iy)=etpfrz(ianp,iy)
2307 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2308 csp07/05
2309             if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then
2310               etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2311               xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0
2312             elseif(ITYP5(I).eq.3)then
2313               etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2314               xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0
2315             endif
2316 csp07/05 end
2317  1011    continue
2318 c     parton freezeout info taken, proceed to next parton:
2319                        goto 350
2320                     endif
2321                  endif
2322  325          continue
2323  350       continue
2324 c     idd=2: calculate average partonic elliptic flow, called from artdri.f,
2325         else if(idd.eq.2) then
2326            do 1012 i=1,3
2327 cms           write(nfile(i),*) '   tsp,   v2p,     v2p2, '//
2328 cms  1 '   s2p,  etp,   xmult,    nevt,  xnctot'
2329 cms           write ((nfile(i)+1),*) '  tsp,   v2pf,   v2pf2, '//
2330 cms  1 '   s2pf, etpf,  xnform,  xcrate'
2331 cms           write ((nfile(i)+2),*) '  tsp,   v2pa,   v2pa2, '//
2332 cms  1 '   s2pa, etpa,  xmult_ap,  xnform,   nevt'
2333 cms           write ((nfile(i)+3),*) '  tsph,  v2ph,   v2ph2, '//
2334 cms  1 '   s2ph, etph,  xmult_(ap/2+h),xmult_ap/2,nevt'
2335 csp
2336 cms        write(nlfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2337 cms        write(nsfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2338 cms        write(nlfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2339 cms        write(nsfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2340 c
2341  1012   continue
2342 clin-3/30/00 ensemble average & variance of v2 (over particles in all events):
2343            do 150 ii=1, 30
2344               if(nevt(ii).eq.0) goto 150
2345               do 1014 iy=1,3
2346                  if(xnpart(ii,iy).gt.1) then
2347                     v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy)
2348                     v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy)
2349      1                    -v2p(ii,iy)**2)/(xnpart(ii,iy)-1))
2350                     s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy)
2351 c xmult and etp are multiplicity and et for an averaged event:
2352                     xmult=dble(xnpart(ii,iy)/dble(nevt(ii)))
2353                     etp(ii,iy)=etp(ii,iy)/dble(nevt(ii))
2354 csp
2355                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2356                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2357 c
2358                     xnctot=0d0
2359                     do 1013 inum=1,ii
2360                        xnctot=xnctot+xncoll(inum)
2361  1013               continue
2362                     if(nevt(1).ne.0) xnctot=xnctot/nevt(1)
2363 cms                 write (nfile(iy),200) tsp(ii),v2p(ii,iy),
2364 cms  1      v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot
2365                  endif
2366                  if(nevt(ii).ne.0)
2367      1                xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii)
2368 c
2369                  if(xnpf(ii,iy).gt.1) then
2370                     v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy)
2371                     v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy)
2372      1                    -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1))
2373                     s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy)
2374                     xnform=dble(xnpf(ii,iy)/dble(nevt(ii)))
2375                     etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii))
2376 csp
2377                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2378                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2379 c
2380 cms                 write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy),
2381 cms  1      v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate
2382                  endif
2383 csp
2384                  if(xnpl(ii,iy).gt.1) then
2385                     v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy)
2386                     s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy)
2387                     xmult=dble(xnpl(ii,iy)/dble(nevt(ii)))
2388                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2389 cms                 write (nlfile(iy),201) tsp(ii),v2pl(ii,iy),
2390 cms  1        s2pl(ii,iy),etpl(ii,iy),xmult
2391                  endif
2392                  if(xnps(ii,iy).gt.1) then
2393                     v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy)
2394                     s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy)
2395                     xmult=dble(xnps(ii,iy)/dble(nevt(ii)))
2396                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2397 cms                 write (nsfile(iy),201) tsp(ii),v2ps(ii,iy),
2398 cms  1        s2ps(ii,iy),etps(ii,iy),xmult
2399                  endif
2400                  if(xnplf(ii,iy).gt.1) then
2401                     v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy)
2402                     s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy)
2403                     xmult=dble(xnplf(ii,iy)/dble(nevt(ii)))
2404                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2405 cms                 write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy),
2406 cms  1        s2plf(ii,iy),etplf(ii,iy),xmult
2407                  endif
2408                  if(xnpsf(ii,iy).gt.1) then
2409                     v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy)
2410                     s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy)
2411                     xmult=dble(xnpsf(ii,iy)/dble(nevt(ii)))
2412                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2413 cms                 write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy),
2414 cms  1        s2psf(ii,iy),etpsf(ii,iy),xmult
2415                  endif
2416 csp-end
2417  1014         continue
2418  150           continue
2419 csp07/05 initial & final mt distrb
2420                scalei=0d0
2421                scalef=0d0
2422                if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT
2423                if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT
2424          do 1016 iy=2,3
2425            yra = 1d0
2426            if(iy .eq. 2)yra = 2d0
2427          do 1015 i=1,50
2428 cms        WRITE(nmfile(iy),251) BMT*dble(I - 0.5),
2429 cms  &     SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra,
2430 cms  &     SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra
2431  1015   continue
2432  1016 continue
2433 csp07/05 end
2434 clin-3/30/00 event-by-event average & variance of v2:
2435            if(nevt(30).ge.1) then
2436               do 1017 iy=1,3
2437                  v2pavg(iy)=v2psum(iy)/nevt(30)
2438                  v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2
2439                  if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0)
2440  1017 continue
2441               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg
2442               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p
2443            endif
2444 clin-11/28/00 for animation:
2445             if(nevent.eq.1) then
2446               do 1018 I=1,mul
2447                  if(FT5(I).le.t2time) then
2448          write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2449                  endif
2450  1018         continue
2451 clin-11/29/00 signal the end of animation file:
2452 c              write(91,*) -10.
2453 c              write(91,*) 0
2454 c              write(92,*) -10.
2455 c              write(92,*) 0
2456 c              close (91)
2457 c              close (92)
2458               close (93)
2459            endif
2460 clin-5/18/01 calculate v2 for active partons:
2461            do 450 ianp=1,30
2462               do 400 iy=1,3
2463                  v2pact=0d0
2464                  v2p2ac=0d0
2465                  s2pact=0d0
2466                  etpact=0d0
2467                  xnacti=0d0
2468                  if(xnpf(ianp,iy).gt.1) then
2469 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2470                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2471                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2472      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2473                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2474                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2475                     xnpact=xnpf(ianp,iy)
2476 c
2477                     do 1019 kanp=1,ianp
2478                        v2pact=v2pact-v2pfrz(kanp,iy)
2479                        v2p2ac=v2p2ac-v2p2fz(kanp,iy)
2480                        s2pact=s2pact-s2pfrz(kanp,iy)
2481                        etpact=etpact-etpfrz(kanp,iy)
2482                        xnpact=xnpact-xnpfrz(kanp,iy)
2483  1019               continue
2484 c     save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2485                     v2ph=v2pact
2486                     v2ph2=v2p2ac
2487                     s2ph=s2pact
2488                     etph=etpact
2489                     xnp2=xnpact/2d0
2490 c
2491                     if(xnpact.gt.1.and.nevt(ianp).ne.0) then
2492                        v2pact=v2pact/xnpact
2493                        v2p2ac=dsqrt((v2p2ac/xnpact
2494      1                    -v2pact**2)/(xnpact-1))
2495                        s2pact=s2pact/xnpact
2496                        xnacti=dble(xnpact/dble(nevt(ianp)))
2497                        etpact=etpact/dble(nevt(ianp))
2498 cms                    write (nfile(iy)+2, 250) tsp(ianp),v2pact,
2499 cms  1 v2p2ac,s2pact,etpact,xnacti,
2500 cms  2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp)
2501                     endif
2502                  endif
2503 c     To calculate combined v2 for active partons plus formed hadrons,
2504 c     add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons:
2505 c     scale the hadron part in case nevt(ianp) != nevent:
2506                  shadr=dble(nevt(ianp))/dble(nevent)
2507                  ianh=ianp
2508                  v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2509                  v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1)
2510      1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr
2511                  s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2512                  etph=etph+eth(ianh,iy)*dble(nevent)*shadr
2513                  xnph=xnpact+xnhadr(ianh,iy)*shadr
2514                  xnp2h=xnp2+xnhadr(ianh,iy)*shadr
2515                  if(xnph.gt.1) then
2516                     v2ph=v2ph/xnph
2517                     v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1))
2518                     s2ph=s2ph/xnph
2519                     etph=etph/dble(nevt(ianp))
2520                     xnp2=xnp2/dble(nevt(ianp))
2521                     xnp2h=xnp2h/dble(nevent)
2522 cms                 if(tsp(ianp).le.(ntmax*dt))
2523 cms  1                    write (nfile(iy)+3, 250) tsp(ianp),v2ph,
2524 cms  2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp)
2525                  endif
2526 c
2527  400              continue
2528  450       continue
2529            do 550 ianp=1,30
2530               do 500 iy=1,3
2531                  v2pact=0d0
2532                  v2p2ac=0d0
2533                  s2pact=0d0
2534                  etpact=0d0
2535                  xnacti=0d0
2536 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2537                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2538                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2539      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2540                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2541                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2542                     xnpact=xnpf(ianp,iy)
2543  500              continue
2544  550           continue
2545 cms        close (620)
2546 cms        close (630)
2547            do 1021 nf=1,3
2548               do 1020 ifile=0,3
2549 cms              close(nfile(nf)+ifile)
2550  1020        continue
2551  1021     continue
2552            do 1022 nf=1,3
2553 cms           close(740+nf)
2554  1022      continue
2555         endif
2556 cyy 200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2)
2557 cyy 210        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2))
2558  240        format(a30,3(2x,f9.5))
2559 cyy 250        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6)
2560 csp
2561 cyy 201        format(2x,f5.2,4(2x,f9.2))
2562 cyy 251        format(5e15.5)
2563 c
2564         return
2565         end
2566
2567 c=======================================================================
2568 c     Calculate flow from formed hadrons, called by art1e.f:
2569 c     Note: numbers in art not in double precision!
2570         subroutine flowh(ct)
2571         PARAMETER (MAXSTR=150001, MAXR=1)
2572         dimension tsh(31)
2573         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2574         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3)
2575         DOUBLE PRECISION  pt2, v2hadr
2576         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2577      1 v2h2(30,3),s2h(30,3)
2578 cc      SAVE /hflow/
2579         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2580 cc      SAVE /ebe/
2581         common /lastt/itimeh,bimp
2582 cc      SAVE /lastt/
2583         COMMON /RUN/ NUM
2584 cc      SAVE /RUN/
2585         COMMON  /AA/      R(3,MAXSTR)
2586 cc      SAVE /AA/
2587         COMMON  /BB/      P(3,MAXSTR)
2588 cc      SAVE /BB/
2589         COMMON  /CC/      E(MAXSTR)
2590 cc      SAVE /CC/
2591         COMMON  /EE/      ID(MAXSTR),LB(MAXSTR)
2592 cc      SAVE /EE/
2593         COMMON  /RR/      MASSR(0:MAXR)
2594 cc      SAVE /RR/
2595         common/anim/nevent,isoft,isflag,izpc
2596 cc      SAVE /anim/
2597         COMMON /AREVT/ IAEVT, IARUN, MISS
2598 cc      SAVE /AREVT/
2599         SAVE
2600 c
2601         do 1001 ii = 1, 31
2602            tsh(ii)=float(ii-1)
2603  1001   continue
2604 c
2605         do 1004 ianh = 1, 30
2606            if ((ct+0.0001).lt.tsh(ianh+1)
2607      1 .and.(ct+0.0001).ge.tsh(ianh)) then
2608               if(ianh.eq.itimeh) goto 101
2609               IA = 0
2610               DO 1002 J = 1, NUM
2611                  mult=MASSR(J)
2612                  IA = IA + MASSR(J - 1)
2613                  DO 100 IC = 1, mult
2614                     I = IA + IC
2615 c     5/04/01 exclude leptons and photons:
2616                     if(iabs(LB(I)-10000).lt.100) goto 100
2617                     px=p(1,i)
2618                     py=p(2,i)
2619                     pt2=dble(px)**2+dble(py)**2
2620 c     2/18/00 Note: e(i) gives the mass in ART:
2621                     ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2)
2622                     RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i)))
2623 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2624 c                rt_now=sqrt(r(1,i)**2+r(2,i)**2)
2625 c                gz_now=r(3,i)
2626 c                zmax=0.1*ct
2627 c                volume=3.1416*(2.**2)*2*zmax
2628 c                if(rt_now.gt.2.or.abs(gz_now).gt.zmax)
2629 c     1               goto 100
2630                     iloop=1
2631                     if(abs(rap).le.1) then
2632                        iloop=2
2633                        if(abs(rap).le.0.5) then
2634                           iloop=3
2635                        endif
2636                     endif
2637                     do 50 iy=1,iloop
2638                        if(pt2.gt.0d0) then
2639                           v2hadr=(dble(px)**2-dble(py)**2)/pt2
2640                           v2h(ianh,iy)=v2h(ianh,iy)+v2hadr
2641                           v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2
2642                           if(dabs(v2hadr).gt.1d0)
2643      1 write(1,*) 'v2hadr>1',v2hadr,px,py
2644                        endif
2645                        xperp2=r(1,I)**2+r(2,I)**2
2646                        if(xperp2.gt.0.)
2647      1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2)
2648                eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2)))
2649 ctest off density:
2650 c               eth(ianh,iy)=eth(ianh,iy)
2651 c     1                  +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume
2652                        xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0
2653  50                    continue
2654  100                 continue
2655  1002         CONTINUE
2656               itimeh=ianh
2657 clin-5/04/01 ebe v2 variables:
2658               if(ianh.eq.30) then
2659                  do 1003 iy=1,3
2660                     nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy))
2661                     if(nhadrn.ne.0) then
2662                v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn)
2663                        v2hsum(iy)=v2hsum(iy)+v2hevt(iy)
2664                        v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2
2665                        v2hp(iy)=v2h(ianh,iy)
2666                        xnhadp(iy)=xnhadr(ianh,iy)
2667                     endif
2668  1003            continue
2669                  write(88, 160) iaevt,v2hevt
2670               endif
2671               goto 101
2672            endif
2673  1004   continue
2674  160        format(i10,3(2x,f9.5))
2675 clin-11/27/00 for animation:
2676  101        if(nevent.eq.1) then
2677            IA = 0
2678            do 1005 J = 1, NUM
2679               mult=MASSR(J)
2680               IA = IA + MASSR(J - 1)
2681 c              write(91,*) ct
2682 c              write(91,*) mult
2683               DO 150 IC = 1, mult
2684                  I = IA + IC
2685 c                 write(91,210) LB(i),r(1,i),r(2,i),r(3,i),
2686 c     1              p(1,i),p(2,i),p(3,i),e(i)
2687  150              continue
2688  1005      continue
2689            return
2690         endif
2691 c 210        format(i6,7(1x,f8.3))
2692         return
2693         end
2694
2695 c=======================================================================
2696         subroutine flowh0(NEVNT,idd)
2697 c
2698         dimension tsh(31)
2699         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2700         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,
2701      1 v2havg(3),varv2h(3)
2702         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2703      1 v2h2(30,3),s2h(30,3)
2704 cc      SAVE /hflow/
2705         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2706 cc      SAVE /ebe/
2707         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2708 cc      SAVE /input1/
2709         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2710      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2711 cc      SAVE /INPUT2/
2712         common /lastt/itimeh,bimp
2713 cc      SAVE /lastt/
2714         SAVE
2715
2716 c     idd=0: initialization for flow analysis, called by artdri.f::
2717         if(idd.eq.0) then
2718            itimeh=0
2719 c
2720            do 1001 ii = 1, 31
2721               tsh(ii)=float(ii-1)
2722  1001      continue
2723 c
2724            do 1003 ii=1,30
2725               do 1002 iy=1,3
2726                  v2h(ii,iy)=0d0
2727                  xnhadr(ii,iy)=0d0
2728                  eth(ii,iy)=0d0
2729                  v2h2(ii,iy)=0d0
2730                  s2h(ii,iy)=0d0
2731  1002         continue
2732  1003      continue
2733            do 1004 iy=1,3
2734               v2hp(iy)=0d0
2735               xnhadp(iy)=0d0
2736               v2hsum(iy)=0d0
2737               v2h2sm(iy)=0d0
2738             if(iy.eq.1) then
2739                nunit=59
2740             elseif(iy.eq.2) then
2741                nunit=68
2742             else
2743                nunit=69
2744             endif
2745               write(nunit,*) '   tsh,   v2h,     v2h2,     s2h, '//
2746      1 ' eth,   xmulth'
2747  1004      continue
2748 c     idd=2: calculate average hadronic elliptic flow, called by artdri.f:
2749         else if(idd.eq.2) then
2750            do 100 ii=1, 30
2751               do 1005 iy=1,3
2752                  if(xnhadr(ii,iy).eq.0) then
2753                     xmulth=0.
2754                  elseif(xnhadr(ii,iy).gt.1) then
2755                     v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy)
2756                     eth(ii,iy)=eth(ii,iy)/dble(NEVNT)
2757                     v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy)
2758      1                    -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1))
2759                     s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy)
2760                     xmulth=sngl(xnhadr(ii,iy)/NEVNT)
2761                  endif
2762              if(iy.eq.1) then
2763                 nunit=59
2764              elseif(iy.eq.2) then
2765                 nunit=68
2766              else
2767                 nunit=69
2768              endif
2769                  if(tsh(ii).le.(ntmax*dt))
2770      1                    write (nunit,200) tsh(ii),v2h(ii,iy),
2771      2      v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth
2772  1005         continue
2773  100           continue
2774 c     event-by-event average & variance of v2h:
2775            do 1006 iy=1,3
2776               v2havg(iy)=v2hsum(iy)/dble(NEVNT)
2777       varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2)
2778  1006 continue
2779            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg
2780            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h
2781         endif
2782  200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2))
2783  240        format(a30,3(2x,f9.5))
2784         return
2785         end
2786
2787 c=======================================================================
2788 c     2/23/00 flow from all initial hadrons just before entering ARTMN:
2789         subroutine iniflw(NEVNT,idd)
2790         PARAMETER (MAXSTR=150001, MAXR=1)
2791         DOUBLE PRECISION  v2i,eti,xmulti,v2mi,s2mi,xmmult,
2792      1       v2bi,s2bi,xbmult
2793         COMMON /RUN/ NUM
2794 cc      SAVE /RUN/
2795         COMMON /ARERC1/MULTI1(MAXR)
2796 cc      SAVE /ARERC1/
2797         COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2798      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2799      &     FT1(MAXSTR, MAXR),
2800      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2801      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2802 cc      SAVE /ARPRC1/
2803         COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
2804 cc      SAVE /iflow/
2805         SAVE
2806 c
2807         if(idd.eq.0) then
2808            v2i=0d0
2809            eti=0d0
2810            xmulti=0d0
2811            v2mi=0d0
2812            s2mi=0d0
2813            xmmult=0d0
2814            v2bi=0d0
2815            s2bi=0d0
2816            xbmult=0d0
2817         else if(idd.eq.1) then
2818            do 1002 J = 1, NUM
2819               do 1001 I = 1, MULTI1(J)
2820                  ITYP = ITYP1(I, J)
2821 c     all hadrons:
2822                  IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100
2823                  xmulti=xmulti+1.d0
2824                  PX = PX1(I, J)
2825                  PY = PY1(I, J)
2826                  XM = XM1(I, J)
2827                  pt2=px**2+py**2
2828                  xh=gx1(I,J)
2829                  yh=gy1(I,J)
2830                  xt2=xh**2+yh**2
2831                  if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2)
2832                  eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2))
2833 c     baryons only:
2834                  IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then
2835                     xbmult=xbmult+1.d0
2836                     if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2)
2837                     if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2)
2838 c     mesons only:
2839                  else
2840                     xmmult=xmmult+1.d0
2841                     if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2)
2842                     if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2)
2843                  endif
2844  100                 continue
2845  1001         continue
2846  1002      continue
2847         else if(idd.eq.2) then
2848            if(xmulti.ne.0) v2i=v2i/xmulti
2849            eti=eti/dble(NEVNT)
2850            xmulti=xmulti/dble(NEVNT)
2851            if(xmmult.ne.0) then
2852               v2mi=v2mi/xmmult
2853               s2mi=s2mi/xmmult
2854            endif
2855            xmmult=xmmult/dble(NEVNT)
2856            if(xbmult.ne.0) then
2857               v2bi=v2bi/xbmult
2858               s2bi=s2bi/xbmult
2859            endif
2860            xbmult=xbmult/dble(NEVNT)
2861         endif
2862 c
2863         return
2864         end
2865
2866 c=======================================================================
2867 c     2/25/00 dN/dt analysis for production (before ZPCMN)
2868 c     and freezeout (right after ZPCMN) for all partons.
2869         subroutine frztm(NEVNT,idd)
2870 c
2871         implicit double precision  (a-h, o-z)
2872         PARAMETER (MAXPTN=400001)
2873         dimension tsf(31)
2874         COMMON /PARA1/ MUL
2875 cc      SAVE /PARA1/
2876         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
2877      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
2878      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
2879 cc      SAVE /prec1/
2880         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
2881      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
2882      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
2883 cc      SAVE /prec2/
2884         COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
2885      & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
2886 cc      SAVE /frzout/
2887         SAVE
2888         data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2889      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2890      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
2891 c
2892         if(idd.eq.0) then
2893            do 1001 ii=1,30
2894               xnprod(ii)=0d0
2895               etprod(ii)=0d0
2896               xnfrz(ii)=0d0
2897               etfrz(ii)=0d0
2898               dnprod(ii)=0d0
2899               detpro(ii)=0d0
2900               dnfrz(ii)=0d0
2901               detfrz(ii)=0d0
2902  1001      continue
2903 cms        OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN')
2904 cms        OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN')
2905         else if(idd.eq.1) then
2906            DO 100 ip = 1, MUL
2907               do 1002 ii=1,30
2908                  eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2)
2909                  eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2)
2910 c     total number and Et produced by time tsf(ii):
2911                  if (ft0(ip).lt.tsf(ii+1)) then
2912                     xnprod(ii)=xnprod(ii)+1d0
2913                     etprod(ii)=etprod(ii)+eth0
2914 c     number and Et produced from time tsf(ii) to tsf(ii+1):
2915                     if (ft0(ip).ge.tsf(ii)) then
2916                        dnprod(ii)=dnprod(ii)+1d0
2917                        detpro(ii)=detpro(ii)+eth0
2918                     endif
2919                  endif
2920 c     total number and Et freezed out by time tsf(ii):
2921                  if (FT5(ip).lt.tsf(ii+1)) then
2922                     xnfrz(ii)=xnfrz(ii)+1d0
2923                     etfrz(ii)=etfrz(ii)+eth2
2924 c     number and Et freezed out from time tsf(ii) to tsf(ii+1):
2925                     if (FT5(ip).ge.tsf(ii)) then
2926                        dnfrz(ii)=dnfrz(ii)+1d0
2927                        detfrz(ii)=detfrz(ii)+eth2
2928                     endif
2929                  endif
2930  1002         continue
2931  100           continue
2932         else if(idd.eq.2) then
2933 cms        write (86,*) '       t,       np,       dnp/dt,      etp '//
2934 cms  1 ' detp/dt'
2935 cms        write (87,*) '       t,       nf,       dnf/dt,      etf '//
2936 cms  1 ' detf/dt'
2937            do 1003 ii=1,30
2938               xnp=xnprod(ii)/dble(NEVNT)
2939               xnf=xnfrz(ii)/dble(NEVNT)
2940               etp=etprod(ii)/dble(NEVNT)
2941               etf=etfrz(ii)/dble(NEVNT)
2942               dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2943               dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2944               detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2945               detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2946 cms           write (86, 200)
2947 cms  1        tsf(ii+1),xnp,dxnp,etp,detp
2948 cms           write (87, 200)
2949 cms  1        tsf(ii+1),xnf,dxnf,etf,detf
2950  1003      continue
2951         endif
2952 cyy 200    format(2x,f9.2,4(2x,f10.2))
2953 c
2954         return
2955         end
2956
2957 c=======================================================================
2958 clin-6/2009 write out initial minijet information
2959 c     before propagating to its formation time:
2960         subroutine minijet_out(BB)
2961         PARAMETER (MAXSTR=150001)
2962         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2963         COMMON/hjcrdn/YP(3,300),YT(3,300)
2964         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2965      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2966      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
2967      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2968      &                PJTE(300,500),PJTM(300,500)
2969         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2970      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2971      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2972         COMMON /AREVT/ IAEVT, IARUN, MISS
2973         common /para7/ ioscar,nsmbbbar,nsmmeson
2974         common/phidcy/iphidcy,pttrig,ntrig,maxmiss
2975         SAVE
2976         ntrig=0
2977         do I = 1, IHNT2(1)
2978            do J = 1, NPJ(I)
2979               pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2)
2980               if(pt.ge.pttrig) ntrig=ntrig+1
2981            enddo
2982         enddo
2983         do I = 1, IHNT2(3)
2984            do J = 1, NTJ(I)
2985               pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2)
2986               if(pt.ge.pttrig) ntrig=ntrig+1
2987            enddo
2988         enddo
2989         do I = 1, NSG
2990            do J = 1, NJSG(I)
2991               pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2)
2992               if(pt.ge.pttrig) ntrig=ntrig+1
2993            enddo
2994         enddo
2995 c     Require at least 1 initial minijet parton above the trigger Pt value:
2996         if(ntrig.eq.0) return
2997
2998 c.....transfer data from HIJING to ZPC
2999         if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
3000         DO 1008 I = 1, IHNT2(1)
3001            DO 1007 J = 1, NPJ(I)
3002               ityp=KFPJ(I,J)
3003 c     write out not only gluons:
3004 c              if(ityp.ne.21) goto 1007
3005               gx=YP(1,I)+0.5*BB
3006               gy=YP(2,I)
3007               gz=0d0
3008               ft=0d0
3009               px=PJPX(I,J)
3010               py=PJPY(I,J)
3011               pz=PJPZ(I,J)
3012               xmass=PJPM(I,J)
3013               if(ioscar.eq.3) then
3014                  if(amax1(abs(gx),abs(gy),
3015      1                abs(gz),abs(ft)).lt.9999) then
3016                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3017                  else
3018                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3019                  endif
3020               endif
3021  1007      CONTINUE
3022  1008   CONTINUE
3023         DO 1010 I = 1, IHNT2(3)
3024            DO 1009 J = 1, NTJ(I)
3025               ityp=KFTJ(I,J)
3026 c              if(ityp.ne.21) goto 1009
3027               gx=YT(1,I)-0.5*BB
3028               gy=YT(2,I)
3029               gz=0d0
3030               ft=0d0
3031               px=PJTX(I,J)
3032               py=PJTY(I,J)
3033               pz=PJTZ(I,J)
3034               xmass=PJTM(I,J)
3035               if(ioscar.eq.3) then
3036                  if(amax1(abs(gx),abs(gy),
3037      1                abs(gz),abs(ft)).lt.9999) then
3038                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3039                  else
3040                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3041                  endif
3042               endif
3043  1009      CONTINUE
3044  1010   CONTINUE
3045         DO 1012 I = 1, NSG
3046            DO 1011 J = 1, NJSG(I)
3047               ityp=K2SG(I,J)
3048 c              if(ityp.ne.21) goto 1011
3049               gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2)))
3050               gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2)))
3051               gz=0d0
3052               ft=0d0
3053               px=PXSG(I,J)
3054               py=PYSG(I,J)
3055               pz=PZSG(I,J)
3056               xmass=PMSG(I,J)
3057               if(ioscar.eq.3) then
3058                  if(amax1(abs(gx),abs(gy),
3059      1                abs(gz),abs(ft)).lt.9999) then
3060                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3061                  else
3062                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3063                  endif
3064               endif
3065  1011      CONTINUE
3066  1012   CONTINUE
3067  200  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2)
3068  201  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2)
3069 c
3070         return
3071         end
3072
3073 c=======================================================================
3074 clin-6/2009 embed back-to-back high-Pt quark/antiquark pair
3075 c     via embedding back-to-back high-Pt pion pair then melting the pion pair
3076 c     by generating the internal quark and antiquark momentum parallel to
3077 c      the pion momentum (in order to produce a high-Pt and a low Pt parton):
3078       subroutine embedHighPt
3079       PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140)
3080       common/embed/iembed,pxqembd,pyqembd,xembd,yembd
3081       COMMON/RNDF77/NSEED
3082       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
3083       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
3084       COMMON /ARPRC/ ITYPAR(MAXSTR),
3085      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
3086      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
3087      &     XMAR(MAXSTR)
3088 cwei      DOUBLE PRECISION PATT
3089       SAVE
3090 c
3091       if(iembed.ne.1) return
3092       ptq=sqrt(pxqembd**2+pyqembd**2)
3093       if(ptq.lt.(pichmass/2.)) then
3094          print *, 'Embedded quark transverse momentum is too small'
3095          stop
3096       endif
3097 c     Randomly embed u/ubar or d/dbar at high Pt:
3098       idqembd=1+int(2*RANART(NSEED))
3099 c     Flavor content for the charged pion that contains the leading quark:
3100       if(idqembd.eq.1) then
3101          idqsoft=-2
3102          idpi=-211
3103       elseif(idqembd.eq.2) then
3104          idqsoft=-1
3105          idpi=211
3106       else
3107          print *, 'Wrong quark flavor embedded'
3108          stop
3109       endif
3110 c     Caculate transverse momentum of the parent charged pion:
3111       xmq=ulmass(idqembd)
3112       xmqsoft=ulmass(idqsoft)
3113       ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq
3114      1     -sqrt((xmq**2+ptq**2)*(pichmass**4
3115      2     -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2)))
3116      3     /(2.*xmq**2)
3117       pxpi=ptpi*pxqembd/ptq
3118       pypi=ptpi*pyqembd/ptq
3119 c     Embedded quark/antiquark are assumed to have pz=0:
3120       pzpi=0.
3121 c     Insert the two parent charged pions,
3122 c     ipion=1 for the pion containing the leading quark,
3123 c     ipion=2 for the pion containing the leading antiquark of the same flavor:
3124       do ipion=1,2
3125          if(ipion.eq.2) then
3126             idpi=-idpi
3127             pxpi=-pxpi
3128             pypi=-pypi
3129             pzpi=-pzpi
3130          endif
3131          NATT=NATT+1
3132          KATT(NATT,1)=idpi
3133          KATT(NATT,2)=40
3134          KATT(NATT,3)=0
3135          PATT(NATT,1)=pxpi
3136          PATT(NATT,2)=pypi
3137          PATT(NATT,3)=pzpi
3138          PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2)
3139          EATT=EATT+PATT(NATT,4)
3140          GXAR(NATT)=xembd
3141          GYAR(NATT)=yembd
3142          GZAR(NATT)=0.
3143          FTAR(NATT)=0.
3144          ITYPAR(NATT)=KATT(NATT,1)
3145          PXAR(NATT)=PATT(NATT,1)
3146          PYAR(NATT)=PATT(NATT,2)
3147          PZAR(NATT)=PATT(NATT,3)
3148          PEAR(NATT)=PATT(NATT,4)
3149          XMAR(NATT)=pichmass
3150       enddo
3151 c
3152       return
3153       end