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      &       .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1618      &       .or.iabs(lb1).eq.16) then
1619            kf=INVFLV(lb1)
1620         else
1621            return
1622         endif
1623 c
1624         IP=1
1625 c     label as undecayed and the only particle in the record:
1626         N=1
1627         K(IP,1)=1
1628         K(IP,3)=0
1629         K(IP,4)=0
1630         K(IP,5)=0
1631 c
1632         K(IP,2)=kf
1633         P(IP,1)=px1
1634         P(IP,2)=py1
1635         P(IP,3)=pz1
1636         em1a=em1
1637 c     eta or omega in ART may be below or too close to (pi+pi-pi0) mass,
1638 c     causing LUDECY error,thus increase their mass ADDM above this thresh,
1639 c     noting that rho (m=0.281) too close to 2pi thrshold fails to decay:
1640         if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then
1641            em1=2*APICH+API0+ADDM
1642 c     rho
1643         elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then
1644            em1=2*APICH+ADDM
1645 c     K*
1646         elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then
1647            em1=APICH+AK0+ADDM
1648 c     Delta created in ART may be below (n+pich) mass, causing LUDECY error:
1649         elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9
1650      1          .and.em1.lt.(APICH+AN+ADDM)) then
1651            em1=APICH+AN+ADDM
1652         endif
1653 c        if(em1.ge.(em1a+0.01)) write (6,*)
1654 c     1       'Mass increase in resdec():',nt,em1-em1a,lb1
1655         e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1656         P(IP,4)=e1
1657         P(IP,5)=em1
1658 clin-5/2008:
1659         dpdecp=dpertp(i1)
1660         call ludecy(IP)
1661 c     add decay time to daughter's formation time at the last timestep:
1662         if(nt.eq.ntmax) then
1663            tau0=hbarc/wid
1664            taudcy=tau0*(-1.)*alog(1.-RANART(NSEED))
1665            ndaut=n-nsav
1666            if(ndaut.le.1) then
1667               write(10,*) 'note: ndaut(<1)=',ndaut
1668               write(89,*) 'note: ndaut(<1)=',ndaut
1669               call lulist(2)
1670               stop
1671             endif
1672 c     lorentz boost:
1673            taudcy=taudcy*e1/em1
1674            tfnl=tfnl+taudcy
1675            xfnl=xfnl+px1/e1*taudcy
1676            yfnl=yfnl+py1/e1*taudcy
1677            zfnl=zfnl+pz1/e1*taudcy
1678 c     at the last timestep, assign rho, K0S or eta (decay daughter)
1679 c     to lb(i1) only (not to lpion) in order to decay them again:
1680            if(n.ge.(nsav+2)) then
1681               do 1001 idau=nsav+2,n
1682                  kdaut=K(idau,2)
1683                  if(kdaut.eq.221.or.kdaut.eq.113
1684      1                .or.kdaut.eq.213.or.kdaut.eq.-213
1685      2                .or.kdaut.eq.310) then
1686 c     switch idau and i1(nsav+1):
1687                     ksave=kdaut
1688                     pxsave=p(idau,1)
1689                     pysave=p(idau,2)
1690                     pzsave=p(idau,3)
1691                     esave=p(idau,4)
1692                     xmsave=p(idau,5)
1693                     K(idau,2)=K(nsav+1,2)
1694                     p(idau,1)=p(nsav+1,1)
1695                     p(idau,2)=p(nsav+1,2)
1696                     p(idau,3)=p(nsav+1,3)
1697                     p(idau,4)=p(nsav+1,4)
1698                     p(idau,5)=p(nsav+1,5)
1699                     K(nsav+1,2)=ksave
1700                     p(nsav+1,1)=pxsave
1701                     p(nsav+1,2)=pysave
1702                     p(nsav+1,3)=pzsave
1703                     p(nsav+1,4)=esave
1704                     p(nsav+1,5)=xmsave
1705 c     note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce
1706 c     eta, but only one daughter may be rho, K0s or eta:
1707                     goto 111
1708                  endif
1709  1001         continue
1710            endif
1711  111       continue
1712 c
1713            enet=0.
1714            do 1002 idau=nsav+1,n
1715               enet=enet+p(idau,4)
1716  1002      continue
1717 c           if(abs(enet-e1).gt.0.02)
1718 c     1          write(93,*) 'resdec(): nt=',nt,enet-e1,lb1
1719         endif
1720
1721 cyy 200    format(a20,3(1x,i6))
1722 cyy 210    format(i6,5(1x,f8.3))
1723 cyy 220    format(a2,i5,5(1x,f8.3))
1724
1725         do 1003 idau=nsav+1,n
1726            kdaut=K(idau,2)
1727            lbdaut=IARFLV(kdaut)
1728 c     K0S and K0L are named K+/K- during hadron cascade, and only
1729 c     at the last timestep they keep their real LB # before output;
1730 c     K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep:
1731            if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310
1732      1          .or.iabs(kdaut).eq.311)) then
1733               if(kdaut.eq.130) then
1734                  lbdaut=22
1735               elseif(kdaut.eq.310) then
1736                  lbdaut=24
1737               elseif(iabs(kdaut).eq.311) then
1738                  if(RANART(NSEED).lt.0.5) then
1739                     lbdaut=22
1740                  else
1741                     lbdaut=24
1742                  endif
1743               endif
1744            endif
1745 c
1746            if(idau.eq.(nsav+1)) then
1747               LB(i1)=lbdaut
1748               E(i1)=p(idau,5)
1749               px1n=p(idau,1)
1750               py1n=p(idau,2)
1751               pz1n=p(idau,3)
1752 clin-5/2008:
1753               dp1n=dpdecp
1754            else
1755               nnn=nnn+1
1756               LPION(NNN,IRUN)=lbdaut
1757               EPION(NNN,IRUN)=p(idau,5)
1758               PPION(1,NNN,IRUN)=p(idau,1)
1759               PPION(2,NNN,IRUN)=p(idau,2)
1760               PPION(3,NNN,IRUN)=p(idau,3)
1761               RPION(1,NNN,IRUN)=xfnl
1762               RPION(2,NNN,IRUN)=yfnl
1763               RPION(3,NNN,IRUN)=zfnl
1764               tfdpi(NNN,IRUN)=tfnl
1765 clin-5/2008:
1766               dppion(NNN,IRUN)=dpdecp
1767            endif
1768  1003   continue
1769 cyy 230    format(a2,i5,5(1x,e8.2))
1770         return
1771         end
1772
1773 c=======================================================================
1774         subroutine inidcy
1775
1776         COMMON/LUJETSA/N,K(9000,5),P(9000,5),V(9000,5)
1777 cc      SAVE /LUJETSA/
1778         common/resdcy/NSAV,iksdcy
1779 cc      SAVE /resdcy/
1780         SAVE
1781         N=1
1782         NSAV=N
1783         return
1784         end
1785
1786 c=======================================================================
1787 clin-6/06/02 local parton freezeout motivated from critical density:
1788         subroutine local(t)
1789 c
1790         implicit double precision  (a-h, o-z)
1791         PARAMETER (MAXPTN=400001)
1792         PARAMETER (r0=1d0)
1793         COMMON /para1/ mul
1794 cc      SAVE /para1/
1795         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1796      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1797      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
1798 cc      SAVE /prec2/
1799         common /frzprc/
1800      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1801      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1802      &       xmfrz(MAXPTN),
1803      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1804 cc      SAVE /frzprc/
1805         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1806 cc      SAVE /prec4/
1807         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1808 cc      SAVE /prec5/
1809         common /coal/dpcoal,drcoal,ecritl
1810 cc      SAVE /coal/
1811         SAVE
1812 c
1813       do 1001 it=1,301
1814          if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then
1815             if(it.eq.itlast) then
1816                return
1817             else
1818                itlast=it
1819                goto 50
1820             endif
1821          endif
1822  1001 continue
1823       write(1,*) 'local time out of range in LOCAL, stop',t,it
1824       stop
1825  50   continue
1826 c
1827       do 200 ip=1,mul
1828 c     skip partons which have frozen out:
1829          if(ifrz(ip).eq.1) goto 200
1830          if(it.eq.301) then
1831 c     freezeout all the left partons beyond the time of 3000 fm:
1832             etcrit=1d6
1833             goto 150
1834          else
1835 c     freezeout when transverse energy density < etcrit:
1836             etcrit=(ecritl*2d0/3d0)
1837          endif
1838 c     skip partons which have not yet formed:
1839          if(t.lt.FT5(ip)) goto 200
1840          rap0=rap(ip)
1841          eta0=eta(ip)
1842          x0=GX5(ip)+vx(ip)*(t-FT5(ip))
1843          y0=GY5(ip)+vy(ip)*(t-FT5(ip))
1844          detdy=0d0
1845          do 100 itest=1,mul
1846 c     skip self and partons which have not yet formed:
1847             if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100
1848             ettest=eta(itest)
1849             xtest=GX5(itest)+vx(itest)*(t-FT5(itest))
1850             ytest=GY5(itest)+vy(itest)*(t-FT5(itest))
1851             drt=sqrt((xtest-x0)**2+(ytest-y0)**2)
1852 c     count partons within drt<1 and -1<(eta-eta0)<1:
1853             if(dabs(ettest-eta0).le.1d0.and.drt.le.r0)
1854      1           detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2
1855      2           +XMASS5(itest)**2)*0.5d0
1856  100     continue
1857          detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0))
1858 c     when density is below critical density for phase transition, freeze out:
1859  150     if(detdy.le.etcrit) then
1860             ifrz(ip)=1
1861             idfrz(ip)=ITYP5(ip)
1862             pxfrz(ip)=PX5(ip)
1863             pyfrz(ip)=PY5(ip)
1864             pzfrz(ip)=PZ5(ip)
1865             efrz(ip)=E5(ip)
1866             xmfrz(ip)=XMASS5(ip)
1867             if(t.gt.FT5(ip)) then
1868                gxfrz(ip)=x0
1869                gyfrz(ip)=y0
1870                gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip))
1871                ftfrz(ip)=t
1872             else
1873 c     if this freezeout time < formation time, use formation time & positions.
1874 c     This ensures the recovery of default hadron when e_crit=infty:
1875                gxfrz(ip)=GX5(ip)
1876                gyfrz(ip)=GY5(ip)
1877                gzfrz(ip)=GZ5(ip)
1878                ftfrz(ip)=FT5(ip)
1879             endif
1880          endif
1881  200  continue
1882 c
1883         return
1884         end
1885
1886 c=======================================================================
1887 clin-6/06/02 initialization for local parton freezeout
1888         subroutine inifrz
1889 c
1890         implicit double precision  (a-h, o-z)
1891         PARAMETER (MAXPTN=400001)
1892         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1893 cc      SAVE /ilist5/
1894         common /frzprc/
1895      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1896      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1897      &       xmfrz(MAXPTN),
1898      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1899 cc      SAVE /frzprc/
1900         SAVE
1901 c
1902 c     for freezeout time 0-10fm, use interval of 0.1fm;
1903 c     for 10-100fm, use interval of 1fm;
1904 c     for 100-1000fm, use interval of 10fm;
1905 c     for 1000-3000fm, use interval of 100fm:
1906         step1=0.1d0
1907         step2=1d0
1908         step3=10d0
1909         step4=100d0
1910 c
1911         do 1001 it=1,101
1912            tfrz(it)=0d0+dble(it-1)*step1
1913  1001 continue
1914         do 1002 it=102,191
1915            tfrz(it)=10d0+dble(it-101)*step2
1916  1002   continue
1917         do 1003 it=192,281
1918            tfrz(it)=100d0+dble(it-191)*step3
1919  1003   continue
1920         do 1004 it=282,301
1921            tfrz(it)=1000d0+dble(it-281)*step4
1922  1004   continue
1923         tfrz(302)=tlarge
1924 c
1925         return
1926         end
1927
1928 clin-5/2009 ctest on v2 analysis
1929 c=======================================================================
1930 c     idd=0,1,2,3 specifies different subroutines for partonic flow analysis.
1931       subroutine flowp(idd)
1932 c
1933         implicit double precision  (a-h, o-z)
1934         real dt
1935         parameter (MAXPTN=400001)
1936 csp
1937         parameter (bmt=0.05d0)
1938         dimension nlfile(3),nsfile(3),nmfile(3)
1939 c
1940         dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3)
1941         dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3)
1942         common /ilist1/
1943      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1944      &     ictype, icsta(MAXPTN),
1945      &     nic(MAXPTN), icels(MAXPTN)
1946 cc      SAVE /ilist1/
1947         COMMON /para1/ mul
1948 cc      SAVE /para1/
1949         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1950      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1951      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
1952 cc      SAVE /prec2/
1953         COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3),
1954      1       s2p(30,3),v2p2(30,3),nevt(30)
1955 cc      SAVE /pflow/
1956         COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3),
1957      1                 xncoll(30),s2pf(30,3),v2pf2(30,3)
1958 cc      SAVE /pflowf/
1959         COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3),
1960      1       s2pfrz(30,3),v2p2fz(30,3),tscatt(31),
1961      2       nevtfz(30),iscatt(30)
1962 cc      SAVE /pfrz/
1963         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
1964      1 v2h2(30,3),s2h(30,3)
1965 cc      SAVE /hflow/
1966         COMMON /AREVT/ IAEVT, IARUN, MISS
1967 cc      SAVE /AREVT/
1968         common/anim/nevent,isoft,isflag,izpc
1969 cc      SAVE /anim/
1970         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1971 cc      SAVE /input1/
1972         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1973      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1974 cc      SAVE /INPUT2/
1975 cc      SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm
1976 cc      SAVE nfile,itanim,nlfile,nsfile,nmfile
1977         SAVE
1978 csp
1979         dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3),
1980      &       etlfrz(30,3),etsfrz(30,3),
1981      &       xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3),
1982      &       xnlfrz(30,3),xnsfrz(30,3),
1983      &       v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3),
1984      &       s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3),
1985      &       DMYil(50,3),DMYfl(50,3),
1986      &       DMYis(50,3),DMYfs(50,3)
1987         data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
1988      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
1989      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
1990 c     idd=0: initialization for flow analysis, called by artdri.f:
1991         if(idd.eq.0) then
1992            nfile(1)=60
1993            nfile(2)=64
1994            nfile(3)=20
1995 cms        OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN')
1996 cms        OPEN (nfile(1)+1,
1997 cms  1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN')
1998 cms        OPEN (nfile(1)+2,
1999 cms  1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN')
2000 cms        OPEN (nfile(1)+3,
2001 cms  1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN')
2002 cms        OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN')
2003 cms        OPEN (nfile(2)+1,
2004 cms  1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN')
2005 cms        OPEN (nfile(2)+2,
2006 cms  1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN')
2007 cms        OPEN (nfile(2)+3,
2008 cms  1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN')
2009 cms        OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN')
2010 cms        OPEN (nfile(3)+1,
2011 cms  1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN')
2012 cms        OPEN (nfile(3)+2,
2013 cms  1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN')
2014 cms        OPEN (nfile(3)+3,
2015 cms  1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN')
2016 cms        OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN')
2017 cms        write(49, *) '    ievt,  v2p,  v2p_y2,   v2p_y1'
2018 c
2019 cms        OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN')
2020 cms        OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN')
2021 cms        OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN')
2022 cms        OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN')
2023 cms        write(88, *) '    ievt,  v2h,  v2h_y2,   v2h_y1'
2024 csp07/05
2025            nlfile(1)=70
2026            nlfile(2)=72
2027            nlfile(3)=74
2028 cms        OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN')
2029 cms        OPEN (nlfile(1)+1,
2030 cms  1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN')
2031 cms        OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN')
2032 cms        OPEN (nlfile(2)+1,
2033 cms  1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN')
2034 cms        OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN')
2035 cms        OPEN (nlfile(3)+1,
2036 cms  1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN')
2037            nsfile(1)=76
2038            nsfile(2)=78
2039            nsfile(3)=80
2040 cms        OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN')
2041 cms        OPEN (nsfile(1)+1,
2042 cms  1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN')
2043 cms        OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN')
2044 cms        OPEN (nsfile(2)+1,
2045 cms  1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN')
2046 cms        OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN')
2047 cms        OPEN (nsfile(3)+1,
2048 cms  1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN')
2049            nmfile(1)=82
2050            nmfile(2)=83
2051            nmfile(3)=84
2052 cms        OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN')
2053 cms        OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN')
2054 cms        OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN')
2055 clin-11/27/00 for animation:
2056            if(nevent.eq.1) then
2057 cms        OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN')
2058 cms        write(91,*) ntmax, dt
2059 cms        OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN')
2060 cms        OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN')
2061            endif
2062 c
2063            itimep=0
2064            itanim=0
2065            iaevtp=0
2066 csp
2067            do 1002 ii=1,50
2068               do 1001 iy=1,3
2069                  DMYil(ii,iy) = 0d0
2070                  DMYfl(ii,iy) = 0d0
2071                  DMYis(ii,iy) = 0d0
2072                  DMYfs(ii,iy) = 0d0
2073  1001         continue
2074  1002      continue
2075 c
2076            do 1003 ii=1,31
2077               tscatt(ii)=0d0
2078  1003      continue
2079            do 1005 ii=1,30
2080               nevt(ii)=0
2081               xncoll(ii)=0d0
2082               nevtfz(ii)=0
2083               iscatt(ii)=0
2084               do 1004 iy=1,3
2085                  v2p(ii,iy)=0d0
2086                  v2p2(ii,iy)=0d0
2087                  s2p(ii,iy)=0d0
2088                  etp(ii,iy)=0d0
2089                  xnpart(ii,iy)=0d0
2090                  v2pf(ii,iy)=0d0
2091                  v2pf2(ii,iy)=0d0
2092                  s2pf(ii,iy)=0d0
2093                  etpf(ii,iy)=0d0
2094                  xnpf(ii,iy)=0d0
2095                  v2pfrz(ii,iy)=0d0
2096                  v2p2fz(ii,iy)=0d0
2097                  s2pfrz(ii,iy)=0d0
2098                  etpfrz(ii,iy)=0d0
2099                  xnpfrz(ii,iy)=0d0
2100 csp07/05
2101                  etpl(ii,iy)=0d0
2102                  etps(ii,iy)=0d0
2103                  etplf(ii,iy)=0d0
2104                  etpsf(ii,iy)=0d0
2105                  etlfrz(ii,iy)=0d0
2106                  etsfrz(ii,iy)=0d0
2107               xnpl(ii,iy)=0d0
2108               xnps(ii,iy)=0d0
2109               xnplf(ii,iy)=0d0
2110               xnpsf(ii,iy)=0d0
2111               xnlfrz(ii,iy)=0d0
2112               xnsfrz(ii,iy)=0d0
2113               v2pl(ii,iy)=0d0
2114               v2ps(ii,iy)=0d0
2115               v2plf(ii,iy)=0d0
2116               v2psf(ii,iy)=0d0
2117               s2pl(ii,iy)=0d0
2118               s2ps(ii,iy)=0d0
2119               s2plf(ii,iy)=0d0
2120               s2psf(ii,iy)=0d0
2121  1004      continue
2122  1005   continue
2123            do 1006 iy=1,3
2124               v2pevt(iy)=0d0
2125               v2pavg(iy)=0d0
2126               varv2p(iy)=0d0
2127               v2pp(iy)=0.d0
2128               xnpp(iy)=0d0
2129               v2psum(iy)=0.d0
2130               v2p2sm(iy)=0.d0
2131  1006      continue
2132 c     idd=1: calculate parton elliptic flow, called by zpc.f:
2133         else if(idd.eq.1) then
2134            t2time = FT5(iscat)
2135            do 1008 ianp = 1, 30
2136               if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then
2137 c     write flow info only once at each fixed time:
2138                  xncoll(ianp)=xncoll(ianp)+1d0
2139 c     to prevent an earlier t2time comes later in the same event
2140 c     and mess up nevt:
2141                  if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101
2142                  nevt(ianp)=nevt(ianp)+1
2143                  tscatt(ianp+1)=t2time
2144                  iscatt(ianp)=1
2145                  nevtfz(ianp)=nevtfz(ianp)+1
2146                  do 100 I=1,mul
2147             ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8))
2148                     pt2=PX5(I)**2+PY5(I)**2
2149 ctest off: initial (pt,y) and (x,y) distribution:
2150 c                    idtime=1
2151 c                    if(ianp.eq.idtime) then
2152 c                       iityp=iabs(ITYP5(I))
2153 c                       if(iityp.eq.1.or.iityp.eq.2) then
2154 c                          write(651,*) dsqrt(pt2),ypartn
2155 c                          write(654,*) GX5(I),GY5(I)
2156 c                       elseif(iityp.eq.1103.or.iityp.eq.2101
2157 c     1 .or.iityp.eq.2103.or.iityp.eq.2203.
2158 c     2 .or.iityp.eq.3101.or.iityp.eq.3103.
2159 c     3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303)
2160 c     4 then
2161 c                          write(652,*) dsqrt(pt2),ypartn
2162 c                          write(655,*) GX5(I),GY5(I)
2163 c                       elseif(iityp.eq.21) then
2164 c                          write(653,*) dsqrt(pt2),ypartn
2165 c                          write(656,*) GX5(I),GY5(I)
2166 c                       endif
2167 c                    endif
2168 ctest-end
2169 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2170 c                    gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2171 c                    gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2172 c                    gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2173 c                    rt_now=dsqrt(gx_now**2+gy_now**2)
2174 c                    zmax=0.1d0*t2time
2175 c                    volume=3.1416d0*(2d0**2)*(2*zmax)
2176 c                    if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax)
2177 c     1                   goto 100
2178 ctest-end
2179                     iloop=1
2180                     if(dabs(ypartn).le.1d0) then
2181                        iloop=2
2182                        if(dabs(ypartn).le.0.5d0) then
2183                           iloop=3
2184                        endif
2185                     endif
2186                     do 50 iy=1,iloop
2187                        if(pt2.gt.0.) then
2188                           v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2189                           if(abs(v2prtn).gt.1.)
2190      1 write(nfile(iy),*) 'v2prtn>1',v2prtn
2191                           v2p(ianp,iy)=v2p(ianp,iy)+v2prtn
2192                           v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2
2193                        endif
2194                        xperp2=GX5(I)**2+GY5(I)**2
2195                        if(xperp2.gt.0.)
2196      1        s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2197                        xnpart(ianp,iy)=xnpart(ianp,iy)+1d0
2198                        etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2199 ctest off density:
2200 c                       etp(ianp,iy)=etp(ianp,iy)
2201 c     1                  +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2202 clin-2/22/00 to write out parton info only for formed ones:
2203                        if(FT5(I).le.t2time) then
2204                           v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn
2205                           v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2
2206                           if(xperp2.gt.0.)
2207      1        s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2208                           xnpf(ianp,iy)=xnpf(ianp,iy)+1d0
2209                   etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2210 ctest off density:
2211 c                  etpf(ianp,iy)=etpf(ianp,iy)
2212 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2213                        endif
2214  50                    continue
2215  100                 continue
2216                  itimep=ianp
2217                  iaevtp=iaevt
2218 clin-3/30/00 ebe v2 variables:
2219                  if(ianp.eq.30) then
2220                     do 1007 iy=1,3
2221                        npartn=IDINT(xnpart(ianp,iy)-xnpp(iy))
2222                        if(npartn.ne.0) then
2223                           v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn
2224                           v2psum(iy)=v2psum(iy)+v2pevt(iy)
2225                           v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2
2226                           v2pp(iy)=v2p(ianp,iy)
2227                           xnpp(iy)=xnpart(ianp,iy)
2228                        endif
2229  1007               continue
2230                     write(49, 160) iaevt,v2pevt
2231                  endif
2232                  goto 101
2233               endif
2234  1008      continue
2235 clin-11/28/00 for animation:
2236  101           if(nevent.eq.1) then
2237               do 110 nt = 1, ntmax
2238                  time1=dble(nt*dt)
2239                  time2=dble((nt+1)*dt)
2240                  if (t2time.lt.time2.and.t2time.ge.time1) then
2241                     if(nt.le.itanim) return
2242 c                    write(92,*) t2time
2243                     iform=0
2244                     do 1009 I=1,mul
2245 c     write out parton info only for formed ones:
2246                        if(FT5(I).le.t2time) then
2247                           iform=iform+1
2248                        endif
2249  1009               continue
2250 c                    write(92,*) iform
2251                     do 120 I=1,mul
2252                        if(FT5(I).le.t2time) then
2253 clin-11/29/00-ctest off calculate parton coordinates after propagation:
2254 c                          gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2255 c                          gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2256 c                          gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2257 c          write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2258 c          write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i),
2259 c     1    PX5(i),PY5(i),PZ5(i),E5(i)
2260 ctest-end
2261                        endif
2262  120                    continue
2263                     itanim=nt
2264                  endif
2265  110              continue
2266            endif
2267 c
2268  140           format(i10,4(2x,f7.2))
2269  160           format(i10,3(2x,f9.5))
2270 c 180           format(i6,8(1x,f7.2))
2271 clin-5/17/01 calculate v2 for active partons (which have not frozen out):
2272 c     idd=3, called at end of zpc.f:
2273         else if(idd.eq.3) then
2274            do 1010 ianp=1,30
2275               if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp)
2276  1010      continue
2277            do 350 I=1,mul
2278               ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8)
2279      1 /(E5(i)-PZ5(i)+1.d-8))
2280               pt2=PX5(I)**2+PY5(I)**2
2281               iloop=1
2282               if(dabs(ypartn).le.1d0) then
2283                  iloop=2
2284                  if(dabs(ypartn).le.0.5d0) then
2285                     iloop=3
2286                  endif
2287               endif
2288 c
2289               do 325 ianp=1,30
2290                  if(iscatt(ianp).ne.0) then
2291                     if(FT5(I).lt.tscatt(ianp+1)
2292      1 .and.FT5(I).ge.tscatt(ianp)) then
2293                        do 1011 iy=1,iloop
2294                           if(pt2.gt.0.) then
2295                              v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2296                              v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn
2297                      v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2
2298                           endif
2299                           xperp2=GX5(I)**2+GY5(I)**2
2300                           if(xperp2.gt.0.) s2pfrz(ianp,iy)=
2301      1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2302         etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2303                           xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0
2304 ctest off density:
2305 c                    etpfrz(ianp,iy)=etpfrz(ianp,iy)
2306 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2307 csp07/05
2308             if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then
2309               etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2310               xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0
2311             elseif(ITYP5(I).eq.3)then
2312               etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2313               xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0
2314             endif
2315 csp07/05 end
2316  1011    continue
2317 c     parton freezeout info taken, proceed to next parton:
2318                        goto 350
2319                     endif
2320                  endif
2321  325          continue
2322  350       continue
2323 c     idd=2: calculate average partonic elliptic flow, called from artdri.f,
2324         else if(idd.eq.2) then
2325            do 1012 i=1,3
2326 cms           write(nfile(i),*) '   tsp,   v2p,     v2p2, '//
2327 cms  1 '   s2p,  etp,   xmult,    nevt,  xnctot'
2328 cms           write ((nfile(i)+1),*) '  tsp,   v2pf,   v2pf2, '//
2329 cms  1 '   s2pf, etpf,  xnform,  xcrate'
2330 cms           write ((nfile(i)+2),*) '  tsp,   v2pa,   v2pa2, '//
2331 cms  1 '   s2pa, etpa,  xmult_ap,  xnform,   nevt'
2332 cms           write ((nfile(i)+3),*) '  tsph,  v2ph,   v2ph2, '//
2333 cms  1 '   s2ph, etph,  xmult_(ap/2+h),xmult_ap/2,nevt'
2334 csp
2335 cms        write(nlfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2336 cms        write(nsfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2337 cms        write(nlfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2338 cms        write(nsfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2339 c
2340  1012   continue
2341 clin-3/30/00 ensemble average & variance of v2 (over particles in all events):
2342            do 150 ii=1, 30
2343               if(nevt(ii).eq.0) goto 150
2344               do 1014 iy=1,3
2345                  if(xnpart(ii,iy).gt.1) then
2346                     v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy)
2347                     v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy)
2348      1                    -v2p(ii,iy)**2)/(xnpart(ii,iy)-1))
2349                     s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy)
2350 c xmult and etp are multiplicity and et for an averaged event:
2351                     xmult=dble(xnpart(ii,iy)/dble(nevt(ii)))
2352                     etp(ii,iy)=etp(ii,iy)/dble(nevt(ii))
2353 csp
2354                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2355                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2356 c
2357                     xnctot=0d0
2358                     do 1013 inum=1,ii
2359                        xnctot=xnctot+xncoll(inum)
2360  1013               continue
2361                     if(nevt(1).ne.0) xnctot=xnctot/nevt(1)
2362 cms                 write (nfile(iy),200) tsp(ii),v2p(ii,iy),
2363 cms  1      v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot
2364                  endif
2365                  if(nevt(ii).ne.0)
2366      1                xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii)
2367 c
2368                  if(xnpf(ii,iy).gt.1) then
2369                     v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy)
2370                     v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy)
2371      1                    -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1))
2372                     s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy)
2373                     xnform=dble(xnpf(ii,iy)/dble(nevt(ii)))
2374                     etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii))
2375 csp
2376                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2377                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2378 c
2379 cms                 write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy),
2380 cms  1      v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate
2381                  endif
2382 csp
2383                  if(xnpl(ii,iy).gt.1) then
2384                     v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy)
2385                     s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy)
2386                     xmult=dble(xnpl(ii,iy)/dble(nevt(ii)))
2387                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2388 cms                 write (nlfile(iy),201) tsp(ii),v2pl(ii,iy),
2389 cms  1        s2pl(ii,iy),etpl(ii,iy),xmult
2390                  endif
2391                  if(xnps(ii,iy).gt.1) then
2392                     v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy)
2393                     s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy)
2394                     xmult=dble(xnps(ii,iy)/dble(nevt(ii)))
2395                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2396 cms                 write (nsfile(iy),201) tsp(ii),v2ps(ii,iy),
2397 cms  1        s2ps(ii,iy),etps(ii,iy),xmult
2398                  endif
2399                  if(xnplf(ii,iy).gt.1) then
2400                     v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy)
2401                     s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy)
2402                     xmult=dble(xnplf(ii,iy)/dble(nevt(ii)))
2403                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2404 cms                 write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy),
2405 cms  1        s2plf(ii,iy),etplf(ii,iy),xmult
2406                  endif
2407                  if(xnpsf(ii,iy).gt.1) then
2408                     v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy)
2409                     s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy)
2410                     xmult=dble(xnpsf(ii,iy)/dble(nevt(ii)))
2411                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2412 cms                 write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy),
2413 cms  1        s2psf(ii,iy),etpsf(ii,iy),xmult
2414                  endif
2415 csp-end
2416  1014         continue
2417  150           continue
2418 csp07/05 initial & final mt distrb
2419                scalei=0d0
2420                scalef=0d0
2421                if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT
2422                if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT
2423          do 1016 iy=2,3
2424            yra = 1d0
2425            if(iy .eq. 2)yra = 2d0
2426          do 1015 i=1,50
2427 cms        WRITE(nmfile(iy),251) BMT*dble(I - 0.5),
2428 cms  &     SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra,
2429 cms  &     SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra
2430  1015   continue
2431  1016 continue
2432 csp07/05 end
2433 clin-3/30/00 event-by-event average & variance of v2:
2434            if(nevt(30).ge.1) then
2435               do 1017 iy=1,3
2436                  v2pavg(iy)=v2psum(iy)/nevt(30)
2437                  v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2
2438                  if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0)
2439  1017 continue
2440               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg
2441               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p
2442            endif
2443 clin-11/28/00 for animation:
2444             if(nevent.eq.1) then
2445               do 1018 I=1,mul
2446                  if(FT5(I).le.t2time) then
2447          write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2448                  endif
2449  1018         continue
2450 clin-11/29/00 signal the end of animation file:
2451 c              write(91,*) -10.
2452 c              write(91,*) 0
2453 c              write(92,*) -10.
2454 c              write(92,*) 0
2455 c              close (91)
2456 c              close (92)
2457               close (93)
2458            endif
2459 clin-5/18/01 calculate v2 for active partons:
2460            do 450 ianp=1,30
2461               do 400 iy=1,3
2462                  v2pact=0d0
2463                  v2p2ac=0d0
2464                  s2pact=0d0
2465                  etpact=0d0
2466                  xnacti=0d0
2467                  if(xnpf(ianp,iy).gt.1) then
2468 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2469                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2470                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2471      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2472                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2473                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2474                     xnpact=xnpf(ianp,iy)
2475 c
2476                     do 1019 kanp=1,ianp
2477                        v2pact=v2pact-v2pfrz(kanp,iy)
2478                        v2p2ac=v2p2ac-v2p2fz(kanp,iy)
2479                        s2pact=s2pact-s2pfrz(kanp,iy)
2480                        etpact=etpact-etpfrz(kanp,iy)
2481                        xnpact=xnpact-xnpfrz(kanp,iy)
2482  1019               continue
2483 c     save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2484                     v2ph=v2pact
2485                     v2ph2=v2p2ac
2486                     s2ph=s2pact
2487                     etph=etpact
2488                     xnp2=xnpact/2d0
2489 c
2490                     if(xnpact.gt.1.and.nevt(ianp).ne.0) then
2491                        v2pact=v2pact/xnpact
2492                        v2p2ac=dsqrt((v2p2ac/xnpact
2493      1                    -v2pact**2)/(xnpact-1))
2494                        s2pact=s2pact/xnpact
2495                        xnacti=dble(xnpact/dble(nevt(ianp)))
2496                        etpact=etpact/dble(nevt(ianp))
2497 cms                    write (nfile(iy)+2, 250) tsp(ianp),v2pact,
2498 cms  1 v2p2ac,s2pact,etpact,xnacti,
2499 cms  2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp)
2500                     endif
2501                  endif
2502 c     To calculate combined v2 for active partons plus formed hadrons,
2503 c     add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons:
2504 c     scale the hadron part in case nevt(ianp) != nevent:
2505                  shadr=dble(nevt(ianp))/dble(nevent)
2506                  ianh=ianp
2507                  v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2508                  v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1)
2509      1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr
2510                  s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2511                  etph=etph+eth(ianh,iy)*dble(nevent)*shadr
2512                  xnph=xnpact+xnhadr(ianh,iy)*shadr
2513                  xnp2h=xnp2+xnhadr(ianh,iy)*shadr
2514                  if(xnph.gt.1) then
2515                     v2ph=v2ph/xnph
2516                     v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1))
2517                     s2ph=s2ph/xnph
2518                     etph=etph/dble(nevt(ianp))
2519                     xnp2=xnp2/dble(nevt(ianp))
2520                     xnp2h=xnp2h/dble(nevent)
2521 cms                 if(tsp(ianp).le.(ntmax*dt))
2522 cms  1                    write (nfile(iy)+3, 250) tsp(ianp),v2ph,
2523 cms  2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp)
2524                  endif
2525 c
2526  400              continue
2527  450       continue
2528            do 550 ianp=1,30
2529               do 500 iy=1,3
2530                  v2pact=0d0
2531                  v2p2ac=0d0
2532                  s2pact=0d0
2533                  etpact=0d0
2534                  xnacti=0d0
2535 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2536                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2537                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2538      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2539                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2540                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2541                     xnpact=xnpf(ianp,iy)
2542  500              continue
2543  550           continue
2544 cms        close (620)
2545 cms        close (630)
2546            do 1021 nf=1,3
2547               do 1020 ifile=0,3
2548 cms              close(nfile(nf)+ifile)
2549  1020        continue
2550  1021     continue
2551            do 1022 nf=1,3
2552 cms           close(740+nf)
2553  1022      continue
2554         endif
2555 cyy 200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2)
2556 cyy 210        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2))
2557  240        format(a30,3(2x,f9.5))
2558 cyy 250        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6)
2559 csp
2560 cyy 201        format(2x,f5.2,4(2x,f9.2))
2561 cyy 251        format(5e15.5)
2562 c
2563         return
2564         end
2565
2566 c=======================================================================
2567 c     Calculate flow from formed hadrons, called by art1e.f:
2568 c     Note: numbers in art not in double precision!
2569         subroutine flowh(ct)
2570         PARAMETER (MAXSTR=150001, MAXR=1)
2571         dimension tsh(31)
2572         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2573         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3)
2574         DOUBLE PRECISION  pt2, v2hadr
2575         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2576      1 v2h2(30,3),s2h(30,3)
2577 cc      SAVE /hflow/
2578         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2579 cc      SAVE /ebe/
2580         common /lastt/itimeh,bimp
2581 cc      SAVE /lastt/
2582         COMMON /RUN/ NUM
2583 cc      SAVE /RUN/
2584         COMMON  /AA/      R(3,MAXSTR)
2585 cc      SAVE /AA/
2586         COMMON  /BB/      P(3,MAXSTR)
2587 cc      SAVE /BB/
2588         COMMON  /CC/      E(MAXSTR)
2589 cc      SAVE /CC/
2590         COMMON  /EE/      ID(MAXSTR),LB(MAXSTR)
2591 cc      SAVE /EE/
2592         COMMON  /RR/      MASSR(0:MAXR)
2593 cc      SAVE /RR/
2594         common/anim/nevent,isoft,isflag,izpc
2595 cc      SAVE /anim/
2596         COMMON /AREVT/ IAEVT, IARUN, MISS
2597 cc      SAVE /AREVT/
2598         SAVE
2599 c
2600         do 1001 ii = 1, 31
2601            tsh(ii)=float(ii-1)
2602  1001   continue
2603 c
2604         do 1004 ianh = 1, 30
2605            if ((ct+0.0001).lt.tsh(ianh+1)
2606      1 .and.(ct+0.0001).ge.tsh(ianh)) then
2607               if(ianh.eq.itimeh) goto 101
2608               IA = 0
2609               DO 1002 J = 1, NUM
2610                  mult=MASSR(J)
2611                  IA = IA + MASSR(J - 1)
2612                  DO 100 IC = 1, mult
2613                     I = IA + IC
2614 c     5/04/01 exclude leptons and photons:
2615                     if(iabs(LB(I)-10000).lt.100) goto 100
2616                     px=p(1,i)
2617                     py=p(2,i)
2618                     pt2=dble(px)**2+dble(py)**2
2619 c     2/18/00 Note: e(i) gives the mass in ART:
2620                     ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2)
2621                     RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i)))
2622 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2623 c                rt_now=sqrt(r(1,i)**2+r(2,i)**2)
2624 c                gz_now=r(3,i)
2625 c                zmax=0.1*ct
2626 c                volume=3.1416*(2.**2)*2*zmax
2627 c                if(rt_now.gt.2.or.abs(gz_now).gt.zmax)
2628 c     1               goto 100
2629                     iloop=1
2630                     if(abs(rap).le.1) then
2631                        iloop=2
2632                        if(abs(rap).le.0.5) then
2633                           iloop=3
2634                        endif
2635                     endif
2636                     do 50 iy=1,iloop
2637                        if(pt2.gt.0d0) then
2638                           v2hadr=(dble(px)**2-dble(py)**2)/pt2
2639                           v2h(ianh,iy)=v2h(ianh,iy)+v2hadr
2640                           v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2
2641                           if(dabs(v2hadr).gt.1d0)
2642      1 write(1,*) 'v2hadr>1',v2hadr,px,py
2643                        endif
2644                        xperp2=r(1,I)**2+r(2,I)**2
2645                        if(xperp2.gt.0.)
2646      1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2)
2647                eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2)))
2648 ctest off density:
2649 c               eth(ianh,iy)=eth(ianh,iy)
2650 c     1                  +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume
2651                        xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0
2652  50                    continue
2653  100                 continue
2654  1002         CONTINUE
2655               itimeh=ianh
2656 clin-5/04/01 ebe v2 variables:
2657               if(ianh.eq.30) then
2658                  do 1003 iy=1,3
2659                     nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy))
2660                     if(nhadrn.ne.0) then
2661                v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn)
2662                        v2hsum(iy)=v2hsum(iy)+v2hevt(iy)
2663                        v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2
2664                        v2hp(iy)=v2h(ianh,iy)
2665                        xnhadp(iy)=xnhadr(ianh,iy)
2666                     endif
2667  1003            continue
2668                  write(88, 160) iaevt,v2hevt
2669               endif
2670               goto 101
2671            endif
2672  1004   continue
2673  160        format(i10,3(2x,f9.5))
2674 clin-11/27/00 for animation:
2675  101        if(nevent.eq.1) then
2676            IA = 0
2677            do 1005 J = 1, NUM
2678               mult=MASSR(J)
2679               IA = IA + MASSR(J - 1)
2680 c              write(91,*) ct
2681 c              write(91,*) mult
2682               DO 150 IC = 1, mult
2683                  I = IA + IC
2684 c                 write(91,210) LB(i),r(1,i),r(2,i),r(3,i),
2685 c     1              p(1,i),p(2,i),p(3,i),e(i)
2686  150              continue
2687  1005      continue
2688            return
2689         endif
2690 c 210        format(i6,7(1x,f8.3))
2691         return
2692         end
2693
2694 c=======================================================================
2695         subroutine flowh0(NEVNT,idd)
2696 c
2697         dimension tsh(31)
2698         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2699         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,
2700      1 v2havg(3),varv2h(3)
2701         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2702      1 v2h2(30,3),s2h(30,3)
2703 cc      SAVE /hflow/
2704         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2705 cc      SAVE /ebe/
2706         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2707 cc      SAVE /input1/
2708         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2709      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2710 cc      SAVE /INPUT2/
2711         common /lastt/itimeh,bimp
2712 cc      SAVE /lastt/
2713         SAVE
2714
2715 c     idd=0: initialization for flow analysis, called by artdri.f::
2716         if(idd.eq.0) then
2717            itimeh=0
2718 c
2719            do 1001 ii = 1, 31
2720               tsh(ii)=float(ii-1)
2721  1001      continue
2722 c
2723            do 1003 ii=1,30
2724               do 1002 iy=1,3
2725                  v2h(ii,iy)=0d0
2726                  xnhadr(ii,iy)=0d0
2727                  eth(ii,iy)=0d0
2728                  v2h2(ii,iy)=0d0
2729                  s2h(ii,iy)=0d0
2730  1002         continue
2731  1003      continue
2732            do 1004 iy=1,3
2733               v2hp(iy)=0d0
2734               xnhadp(iy)=0d0
2735               v2hsum(iy)=0d0
2736               v2h2sm(iy)=0d0
2737             if(iy.eq.1) then
2738                nunit=59
2739             elseif(iy.eq.2) then
2740                nunit=68
2741             else
2742                nunit=69
2743             endif
2744               write(nunit,*) '   tsh,   v2h,     v2h2,     s2h, '//
2745      1 ' eth,   xmulth'
2746  1004      continue
2747 c     idd=2: calculate average hadronic elliptic flow, called by artdri.f:
2748         else if(idd.eq.2) then
2749            do 100 ii=1, 30
2750               do 1005 iy=1,3
2751                  if(xnhadr(ii,iy).eq.0) then
2752                     xmulth=0.
2753                  elseif(xnhadr(ii,iy).gt.1) then
2754                     v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy)
2755                     eth(ii,iy)=eth(ii,iy)/dble(NEVNT)
2756                     v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy)
2757      1                    -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1))
2758                     s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy)
2759                     xmulth=sngl(xnhadr(ii,iy)/NEVNT)
2760                  endif
2761              if(iy.eq.1) then
2762                 nunit=59
2763              elseif(iy.eq.2) then
2764                 nunit=68
2765              else
2766                 nunit=69
2767              endif
2768                  if(tsh(ii).le.(ntmax*dt))
2769      1                    write (nunit,200) tsh(ii),v2h(ii,iy),
2770      2      v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth
2771  1005         continue
2772  100           continue
2773 c     event-by-event average & variance of v2h:
2774            do 1006 iy=1,3
2775               v2havg(iy)=v2hsum(iy)/dble(NEVNT)
2776       varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2)
2777  1006 continue
2778            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg
2779            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h
2780         endif
2781  200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2))
2782  240        format(a30,3(2x,f9.5))
2783         return
2784         end
2785
2786 c=======================================================================
2787 c     2/23/00 flow from all initial hadrons just before entering ARTMN:
2788         subroutine iniflw(NEVNT,idd)
2789         PARAMETER (MAXSTR=150001, MAXR=1)
2790         DOUBLE PRECISION  v2i,eti,xmulti,v2mi,s2mi,xmmult,
2791      1       v2bi,s2bi,xbmult
2792         COMMON /RUN/ NUM
2793 cc      SAVE /RUN/
2794         COMMON /ARERC1/MULTI1(MAXR)
2795 cc      SAVE /ARERC1/
2796         COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2797      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2798      &     FT1(MAXSTR, MAXR),
2799      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2800      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2801 cc      SAVE /ARPRC1/
2802         COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
2803 cc      SAVE /iflow/
2804         SAVE
2805 c
2806         if(idd.eq.0) then
2807            v2i=0d0
2808            eti=0d0
2809            xmulti=0d0
2810            v2mi=0d0
2811            s2mi=0d0
2812            xmmult=0d0
2813            v2bi=0d0
2814            s2bi=0d0
2815            xbmult=0d0
2816         else if(idd.eq.1) then
2817            do 1002 J = 1, NUM
2818               do 1001 I = 1, MULTI1(J)
2819                  ITYP = ITYP1(I, J)
2820 c     all hadrons:
2821                  IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100
2822                  xmulti=xmulti+1.d0
2823                  PX = PX1(I, J)
2824                  PY = PY1(I, J)
2825                  XM = XM1(I, J)
2826                  pt2=px**2+py**2
2827                  xh=gx1(I,J)
2828                  yh=gy1(I,J)
2829                  xt2=xh**2+yh**2
2830                  if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2)
2831                  eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2))
2832 c     baryons only:
2833                  IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then
2834                     xbmult=xbmult+1.d0
2835                     if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2)
2836                     if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2)
2837 c     mesons only:
2838                  else
2839                     xmmult=xmmult+1.d0
2840                     if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2)
2841                     if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2)
2842                  endif
2843  100                 continue
2844  1001         continue
2845  1002      continue
2846         else if(idd.eq.2) then
2847            if(xmulti.ne.0) v2i=v2i/xmulti
2848            eti=eti/dble(NEVNT)
2849            xmulti=xmulti/dble(NEVNT)
2850            if(xmmult.ne.0) then
2851               v2mi=v2mi/xmmult
2852               s2mi=s2mi/xmmult
2853            endif
2854            xmmult=xmmult/dble(NEVNT)
2855            if(xbmult.ne.0) then
2856               v2bi=v2bi/xbmult
2857               s2bi=s2bi/xbmult
2858            endif
2859            xbmult=xbmult/dble(NEVNT)
2860         endif
2861 c
2862         return
2863         end
2864
2865 c=======================================================================
2866 c     2/25/00 dN/dt analysis for production (before ZPCMN)
2867 c     and freezeout (right after ZPCMN) for all partons.
2868         subroutine frztm(NEVNT,idd)
2869 c
2870         implicit double precision  (a-h, o-z)
2871         PARAMETER (MAXPTN=400001)
2872         dimension tsf(31)
2873         COMMON /PARA1/ MUL
2874 cc      SAVE /PARA1/
2875         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
2876      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
2877      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
2878 cc      SAVE /prec1/
2879         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
2880      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
2881      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
2882 cc      SAVE /prec2/
2883         COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
2884      & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
2885 cc      SAVE /frzout/
2886         SAVE
2887         data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2888      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2889      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
2890 c
2891         if(idd.eq.0) then
2892            do 1001 ii=1,30
2893               xnprod(ii)=0d0
2894               etprod(ii)=0d0
2895               xnfrz(ii)=0d0
2896               etfrz(ii)=0d0
2897               dnprod(ii)=0d0
2898               detpro(ii)=0d0
2899               dnfrz(ii)=0d0
2900               detfrz(ii)=0d0
2901  1001      continue
2902 cms        OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN')
2903 cms        OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN')
2904         else if(idd.eq.1) then
2905            DO 100 ip = 1, MUL
2906               do 1002 ii=1,30
2907                  eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2)
2908                  eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2)
2909 c     total number and Et produced by time tsf(ii):
2910                  if (ft0(ip).lt.tsf(ii+1)) then
2911                     xnprod(ii)=xnprod(ii)+1d0
2912                     etprod(ii)=etprod(ii)+eth0
2913 c     number and Et produced from time tsf(ii) to tsf(ii+1):
2914                     if (ft0(ip).ge.tsf(ii)) then
2915                        dnprod(ii)=dnprod(ii)+1d0
2916                        detpro(ii)=detpro(ii)+eth0
2917                     endif
2918                  endif
2919 c     total number and Et freezed out by time tsf(ii):
2920                  if (FT5(ip).lt.tsf(ii+1)) then
2921                     xnfrz(ii)=xnfrz(ii)+1d0
2922                     etfrz(ii)=etfrz(ii)+eth2
2923 c     number and Et freezed out from time tsf(ii) to tsf(ii+1):
2924                     if (FT5(ip).ge.tsf(ii)) then
2925                        dnfrz(ii)=dnfrz(ii)+1d0
2926                        detfrz(ii)=detfrz(ii)+eth2
2927                     endif
2928                  endif
2929  1002         continue
2930  100           continue
2931         else if(idd.eq.2) then
2932 cms        write (86,*) '       t,       np,       dnp/dt,      etp '//
2933 cms  1 ' detp/dt'
2934 cms        write (87,*) '       t,       nf,       dnf/dt,      etf '//
2935 cms  1 ' detf/dt'
2936            do 1003 ii=1,30
2937               xnp=xnprod(ii)/dble(NEVNT)
2938               xnf=xnfrz(ii)/dble(NEVNT)
2939               etp=etprod(ii)/dble(NEVNT)
2940               etf=etfrz(ii)/dble(NEVNT)
2941               dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2942               dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2943               detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2944               detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
2945 cms           write (86, 200)
2946 cms  1        tsf(ii+1),xnp,dxnp,etp,detp
2947 cms           write (87, 200)
2948 cms  1        tsf(ii+1),xnf,dxnf,etf,detf
2949  1003      continue
2950         endif
2951 cyy 200    format(2x,f9.2,4(2x,f10.2))
2952 c
2953         return
2954         end
2955
2956 c=======================================================================
2957 clin-6/2009 write out initial minijet information
2958 c     before propagating to its formation time:
2959         subroutine minijet_out(BB)
2960         PARAMETER (MAXSTR=150001)
2961         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
2962         COMMON/hjcrdn/YP(3,300),YT(3,300)
2963         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
2964      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
2965      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
2966      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
2967      &                PJTE(300,500),PJTM(300,500)
2968         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
2969      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
2970      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
2971         COMMON /AREVT/ IAEVT, IARUN, MISS
2972         common /para7/ ioscar,nsmbbbar,nsmmeson
2973         common/phidcy/iphidcy,pttrig,ntrig,maxmiss
2974         SAVE
2975         ntrig=0
2976         do I = 1, IHNT2(1)
2977            do J = 1, NPJ(I)
2978               pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2)
2979               if(pt.ge.pttrig) ntrig=ntrig+1
2980            enddo
2981         enddo
2982         do I = 1, IHNT2(3)
2983            do J = 1, NTJ(I)
2984               pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2)
2985               if(pt.ge.pttrig) ntrig=ntrig+1
2986            enddo
2987         enddo
2988         do I = 1, NSG
2989            do J = 1, NJSG(I)
2990               pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2)
2991               if(pt.ge.pttrig) ntrig=ntrig+1
2992            enddo
2993         enddo
2994 c     Require at least 1 initial minijet parton above the trigger Pt value:
2995         if(ntrig.eq.0) return
2996
2997 c.....transfer data from HIJING to ZPC
2998         if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
2999         DO 1008 I = 1, IHNT2(1)
3000            DO 1007 J = 1, NPJ(I)
3001               ityp=KFPJ(I,J)
3002 c     write out not only gluons:
3003 c              if(ityp.ne.21) goto 1007
3004               gx=YP(1,I)+0.5*BB
3005               gy=YP(2,I)
3006               gz=0d0
3007               ft=0d0
3008               px=PJPX(I,J)
3009               py=PJPY(I,J)
3010               pz=PJPZ(I,J)
3011               xmass=PJPM(I,J)
3012               if(ioscar.eq.3) then
3013                  if(amax1(abs(gx),abs(gy),
3014      1                abs(gz),abs(ft)).lt.9999) then
3015                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3016                  else
3017                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3018                  endif
3019               endif
3020  1007      CONTINUE
3021  1008   CONTINUE
3022         DO 1010 I = 1, IHNT2(3)
3023            DO 1009 J = 1, NTJ(I)
3024               ityp=KFTJ(I,J)
3025 c              if(ityp.ne.21) goto 1009
3026               gx=YT(1,I)-0.5*BB
3027               gy=YT(2,I)
3028               gz=0d0
3029               ft=0d0
3030               px=PJTX(I,J)
3031               py=PJTY(I,J)
3032               pz=PJTZ(I,J)
3033               xmass=PJTM(I,J)
3034               if(ioscar.eq.3) then
3035                  if(amax1(abs(gx),abs(gy),
3036      1                abs(gz),abs(ft)).lt.9999) then
3037                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3038                  else
3039                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3040                  endif
3041               endif
3042  1009      CONTINUE
3043  1010   CONTINUE
3044         DO 1012 I = 1, NSG
3045            DO 1011 J = 1, NJSG(I)
3046               ityp=K2SG(I,J)
3047 c              if(ityp.ne.21) goto 1011
3048               gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2)))
3049               gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2)))
3050               gz=0d0
3051               ft=0d0
3052               px=PXSG(I,J)
3053               py=PYSG(I,J)
3054               pz=PZSG(I,J)
3055               xmass=PMSG(I,J)
3056               if(ioscar.eq.3) then
3057                  if(amax1(abs(gx),abs(gy),
3058      1                abs(gz),abs(ft)).lt.9999) then
3059                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3060                  else
3061                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3062                  endif
3063               endif
3064  1011      CONTINUE
3065  1012   CONTINUE
3066  200  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2)
3067  201  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2)
3068 c
3069         return
3070         end
3071
3072 c=======================================================================
3073 clin-6/2009 embed back-to-back high-Pt quark/antiquark pair
3074 c     via embedding back-to-back high-Pt pion pair then melting the pion pair
3075 c     by generating the internal quark and antiquark momentum parallel to
3076 c      the pion momentum (in order to produce a high-Pt and a low Pt parton):
3077       subroutine embedHighPt
3078       PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140)
3079       common/embed/iembed,pxqembd,pyqembd,xembd,yembd
3080       COMMON/RNDF77/NSEED
3081       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
3082       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
3083       COMMON /ARPRC/ ITYPAR(MAXSTR),
3084      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
3085      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
3086      &     XMAR(MAXSTR)
3087 cwei      DOUBLE PRECISION PATT
3088       SAVE
3089 c
3090       if(iembed.ne.1) return
3091       ptq=sqrt(pxqembd**2+pyqembd**2)
3092       if(ptq.lt.(pichmass/2.)) then
3093          print *, 'Embedded quark transverse momentum is too small'
3094          stop
3095       endif
3096 c     Randomly embed u/ubar or d/dbar at high Pt:
3097       idqembd=1+int(2*RANART(NSEED))
3098 c     Flavor content for the charged pion that contains the leading quark:
3099       if(idqembd.eq.1) then
3100          idqsoft=-2
3101          idpi=-211
3102       elseif(idqembd.eq.2) then
3103          idqsoft=-1
3104          idpi=211
3105       else
3106          print *, 'Wrong quark flavor embedded'
3107          stop
3108       endif
3109 c     Caculate transverse momentum of the parent charged pion:
3110       xmq=ulmass(idqembd)
3111       xmqsoft=ulmass(idqsoft)
3112       ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq
3113      1     -sqrt((xmq**2+ptq**2)*(pichmass**4
3114      2     -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2)))
3115      3     /(2.*xmq**2)
3116       pxpi=ptpi*pxqembd/ptq
3117       pypi=ptpi*pyqembd/ptq
3118 c     Embedded quark/antiquark are assumed to have pz=0:
3119       pzpi=0.
3120 c     Insert the two parent charged pions,
3121 c     ipion=1 for the pion containing the leading quark,
3122 c     ipion=2 for the pion containing the leading antiquark of the same flavor:
3123       do ipion=1,2
3124          if(ipion.eq.2) then
3125             idpi=-idpi
3126             pxpi=-pxpi
3127             pypi=-pypi
3128             pzpi=-pzpi
3129          endif
3130          NATT=NATT+1
3131          KATT(NATT,1)=idpi
3132          KATT(NATT,2)=40
3133          KATT(NATT,3)=0
3134          PATT(NATT,1)=pxpi
3135          PATT(NATT,2)=pypi
3136          PATT(NATT,3)=pzpi
3137          PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2)
3138          EATT=EATT+PATT(NATT,4)
3139          GXAR(NATT)=xembd
3140          GYAR(NATT)=yembd
3141          GZAR(NATT)=0.
3142          FTAR(NATT)=0.
3143          ITYPAR(NATT)=KATT(NATT,1)
3144          PXAR(NATT)=PATT(NATT,1)
3145          PYAR(NATT)=PATT(NATT,2)
3146          PZAR(NATT)=PATT(NATT,3)
3147          PEAR(NATT)=PATT(NATT,4)
3148          XMAR(NATT)=pichmass
3149       enddo
3150 c
3151       return
3152       end