AMPTInterface/src/linana.f

0001 c.................... linana.f
0002 c=======================================================================
0003 c     10/26/01 update freezeout positions in case of interactions:
0004 clin-3/2009 Note: freezeout spacetime values cannot be trusted for K0S & K0L
0005 c     as K0S/K0L are converted from K+/K- by hand at the end of hadron cascade.
0006       subroutine hbtout(nnew,nt,ntmax)
0007 c
0008       PARAMETER  (MAXSTR=150001,MAXR=1)
0009 clin-5/2008 give tolerance to regular particles (perturbative probability 1):
0010       PARAMETER  (oneminus=0.99999,oneplus=1.00001)
0011       dimension lastkp(MAXSTR), newkp(MAXSTR),xnew(3)
0012       common /para7/ ioscar,nsmbbbar,nsmmeson
0013 cc      SAVE /para7/
0014       COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
0015 cc      SAVE /hbt/
0016       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0017 cc      SAVE /input1/
0018       COMMON   /AA/  R(3,MAXSTR)
0019 cc      SAVE /AA/
0020       COMMON   /BB/  P(3,MAXSTR)
0021 cc      SAVE /BB/
0022       COMMON   /CC/  E(MAXSTR)
0023 cc      SAVE /CC/
0024       COMMON   /EE/  ID(MAXSTR),LB(MAXSTR)
0025 cc      SAVE /EE/
0026       common /lastt/itimeh,bimp
0027 cc      SAVE /lastt/
0028       COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
0029 cc      SAVE /tdecay/
0030       COMMON /AREVT/ IAEVT, IARUN, MISS
0031 cc      SAVE /AREVT/
0032       common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
0033 cc      SAVE /snn/
0034       COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
0035 cc      SAVE /HJGLBR/
0036       COMMON/FTMAX/ftsv(MAXSTR),ftsvt(MAXSTR, MAXR)
0037       COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
0038      1     dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
0039      2     dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
0040 clin-12/14/03:
0041       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0042       EXTERNAL IARFLV, INVFLV
0043       common /para8/ idpert,npertd,idxsec
0044 clin-2/2012:
0045       common /phiHJ/iphirp,phiRP
0046       SAVE
0047 c
0048       do 1001 i=1,max0(nlast,nnew)
0049          lastkp(i)=0
0050  1001 continue
0051       do 1002 i=1,nnew
0052          newkp(i)=0
0053  1002 continue
0054 c     for each of the particles, search the freezeout record (common /hbt/)
0055 c     to find & keep those which do not have interactions during this timestep:
0056       do 100 ip=1,nnew
0057          do 1004 iplast=1,nlast
0058             if(p(1,ip).eq.plast(1,iplast).and.
0059      1           p(2,ip).eq.plast(2,iplast).and.
0060      2           p(3,ip).eq.plast(3,iplast).and.
0061      3           e(ip).eq.plast(4,iplast).and.
0062      4           lb(ip).eq.lblast(iplast).and.
0063      5      dpertp(ip).eq.dplast(iplast).and.lastkp(iplast).eq.0) then
0064 clin-5/2008 modified below to the above in case we have perturbative particles:
0065 c     5           lastkp(iplast).eq.0) then
0066                deltat=nt*dt-xlast(4,iplast)
0067                ene=sqrt(plast(1,iplast)**2+plast(2,iplast)**2
0068      1              +plast(3,iplast)**2+plast(4,iplast)**2)
0069 c     xnew gives the coordinate if a particle free-streams to current time:
0070                do 1003 ii=1,3
0071                   xnew(ii)=xlast(ii,iplast)+plast(ii,iplast)/ene*deltat
0072  1003          continue
0073                   dr=sqrt((r(1,ip)-xnew(1))**2+(r(2,ip)-xnew(2))**2
0074      1              +(r(3,ip)-xnew(3))**2)
0075 c     find particles with dp=0 and dr<0.01, considered to be those
0076 c     without any interactions during this timestep,
0077 c     thus keep their last positions and time:
0078                if(dr.le.0.01) then
0079                   lastkp(iplast)=1
0080                   newkp(ip)=1
0081 c                  if(lb(ip).eq.41) then
0082 c                write(95,*) 'nt,ip,px,x=',nt,ip,p(1,ip),r(1,ip),ftsv(ip)
0083 c                write(95,*) 'xnew=',xnew(1),xnew(2),xnew(3),xlast(4,ip)
0084 c                  endif
0085 clin-5/2009 Take care of formation time of particles read in at nt=ntmax-1:
0086                   if(nt.eq.ntmax.and.ftsv(ip).gt.((ntmax-1)*dt))
0087      1                 xlast(4,iplast)=ftsv(ip)
0088                   goto 100
0089                endif
0090             endif
0091  1004    continue
0092  100  continue
0093 c     for current particles with interactions, fill their current info in
0094 c     the freezeout record (if that record entry needs not to be kept):
0095       do 150 ip=1,nnew
0096          if(newkp(ip).eq.0) then
0097             do 1005 iplast=1,nnew
0098                if(lastkp(iplast).eq.0) then
0099 ctest off: write collision info
0100 c                  if(lb(ip).eq.41) then
0101 c                     write(95,*) 'nt,lb(ip)=',nt,lb(ip)
0102 c                  write(95,*) '  last p=',plast(1,iplast),
0103 c     1 plast(2,iplast),plast(3,iplast),plast(4,iplast)
0104 c                  write(95,*) '  after p=',p(1,ip),p(2,ip),p(3,ip),e(ip)
0105 c                  write(95,*) 'after x=',r(1,ip),r(2,ip),r(3,ip),ftsv(ip)
0106 c                  endif
0107 c
0108                   xlast(1,iplast)=r(1,ip)
0109                   xlast(2,iplast)=r(2,ip)
0110                   xlast(3,iplast)=r(3,ip)
0111                   xlast(4,iplast)=nt*dt
0112 c
0113                   if(nt.eq.ntmax) then
0114 c     freezeout time for decay daughters at the last timestep
0115 c     needs to include the decay time of the parent:
0116                      if(tfdcy(ip).gt.(ntmax*dt+0.001)) then
0117                         xlast(4,iplast)=tfdcy(ip)
0118 c     freezeout time for particles unformed at the next-to-last timestep
0119 c     needs to be their formation time instead of (ntmax*dt):
0120                      elseif(ftsv(ip).gt.((ntmax-1)*dt)) then
0121                         xlast(4,iplast)=ftsv(ip)
0122                      endif
0123                   endif
0124                   plast(1,iplast)=p(1,ip)
0125                   plast(2,iplast)=p(2,ip)
0126                   plast(3,iplast)=p(3,ip)
0127                   plast(4,iplast)=e(ip)
0128                   lblast(iplast)=lb(ip)
0129                   lastkp(iplast)=1
0130 clin-5/2008:
0131                   dplast(iplast)=dpertp(ip)
0132                   goto 150
0133                endif
0134  1005       continue
0135          endif
0136  150  continue
0137 c     if the current particle list is shorter than the freezeout record,
0138 c     condense the last-collision record by filling new record from 1 to nnew,
0139 c     and label these entries as keep:
0140       if(nnew.lt.nlast) then
0141          do 170 iplast=1,nlast
0142             if(lastkp(iplast).eq.0) then
0143                do 1006 ip2=iplast+1,nlast
0144                   if(lastkp(ip2).eq.1) then
0145                      xlast(1,iplast)=xlast(1,ip2)
0146                      xlast(2,iplast)=xlast(2,ip2)
0147                      xlast(3,iplast)=xlast(3,ip2)
0148                      xlast(4,iplast)=xlast(4,ip2)
0149                      plast(1,iplast)=plast(1,ip2)
0150                      plast(2,iplast)=plast(2,ip2)
0151                      plast(3,iplast)=plast(3,ip2)
0152                      plast(4,iplast)=plast(4,ip2)
0153                      lblast(iplast)=lblast(ip2)
0154                      lastkp(iplast)=1
0155 clin-5/2008:
0156                      dplast(iplast)=dplast(ip2)
0157                      goto 170
0158                   endif
0159  1006          continue
0160             endif
0161  170     continue
0162       endif
0163       nlast=nnew
0164 ctest off look inside each NT timestep (for debugging purpose):
0165 c      do ip=1,nlast
0166 c         write(99,*) ' p ',nt,ip,lblast(ip),plast(1,ip),
0167 c     1        plast(2,ip),plast(3,ip),plast(4,ip),dplast(ip)
0168 c         write(99,*) '  x ',nt,ip,lblast(ip),xlast(1,ip),
0169 c     1        xlast(2,ip),xlast(3,ip),xlast(4,ip),dplast(ip)
0170 c      enddo
0171 c
0172       if(nt.eq.ntmax) then
0173 clin-5/2008 find final number of perturbative particles (deuterons only):
0174          ndpert=0
0175          do ip=1,nlast
0176             if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
0177             else
0178                ndpert=ndpert+1
0179             endif
0180          enddo
0181 c
0182 c         write(16,190) IAEVT,IARUN,nlast,bimp,npart1,npart2,
0183 c     1 NELP,NINP,NELT,NINTHJ
0184 clin-2/2012:
0185 c         write(16,190) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
0186 c     1 NELP,NINP,NELT,NINTHJ
0187          write(16,191) IAEVT,IARUN,nlast-ndpert,bimp,npart1,npart2,
0188      1 NELP,NINP,NELT,NINTHJ,phiRP
0189 clin-5/2008 write out perturbatively-produced particles (deuterons only):
0190          if(idpert.eq.1.or.idpert.eq.2)
0191      1        write(90,190) IAEVT,IARUN,ndpert,bimp,npart1,npart2,
0192      2        NELP,NINP,NELT,NINTHJ
0193          do 1007 ip=1,nlast
0194 clin-12/14/03   No formation time for spectator projectile or target nucleons,
0195 c     see ARINI1 in 'amptsub.f':
0196 clin-3/2009 To be consistent with new particles produced in hadron cascade
0197 c     that are limited by the time-resolution (DT) of the hadron cascade,
0198 c     freezeout time of spectator projectile or target nucleons is written as
0199 c     DT as they are read at the 1st timestep and then propagated to time DT:
0200 c
0201 clin-9/2011 determine spectator nucleons consistently
0202 c            if(plast(1,ip).eq.0.and.plast(2,ip).eq.0
0203 c     1           .and.(sqrt(plast(3,ip)**2+plast(4,ip)**2)*2/HINT1(1))
0204 c     2           .gt.0.99.and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
0205             if(abs(plast(1,ip)).le.epsiPt.and.abs(plast(2,ip)).le.epsiPt
0206      1           .and.(plast(3,ip).gt.amax1(0.,PZPROJ-epsiPz)
0207      2                .or.plast(3,ip).lt.(-PZTARG+epsiPz))
0208      3           .and.(lblast(ip).eq.1.or.lblast(ip).eq.2)) then
0209 clin-5/2008 perturbatively-produced particles (currently only deuterons)
0210 c     are written to ana/ampt_pert.dat (without the column for the mass);
0211 c     ana/ampt.dat has regularly-produced particles (including deuterons);
0212 c     these two sets of deuteron data are close to each other(but not the same
0213 c     because of the bias from triggering the perturbative production);
0214 c     ONLY use one data set for analysis to avoid double-counting:
0215                if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
0216                   write(16,200) INVFLV(lblast(ip)), plast(1,ip),
0217      1                 plast(2,ip),plast(3,ip),plast(4,ip),
0218      2                 xlast(1,ip),xlast(2,ip),xlast(3,ip),
0219      3                 xlast(4,ip)
0220 clin-12/14/03-end
0221                else
0222                   if(idpert.eq.1.or.idpert.eq.2) then
0223                      write(90,250) INVFLV(lblast(ip)), plast(1,ip),
0224      1                 plast(2,ip),plast(3,ip),
0225      2                 xlast(1,ip),xlast(2,ip),xlast(3,ip),
0226      3                 xlast(4,ip)
0227                   else
0228                      write(99,*) 'Unexpected perturbative particles'
0229                   endif
0230                endif
0231             elseif(amax1(abs(xlast(1,ip)),abs(xlast(2,ip)),
0232      1              abs(xlast(3,ip)),abs(xlast(4,ip))).lt.9999) then
0233                if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
0234                  write(16,200) INVFLV(lblast(ip)), plast(1,ip),
0235      1           plast(2,ip),plast(3,ip),plast(4,ip),
0236      2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
0237                else
0238                   if(idpert.eq.1.or.idpert.eq.2) then
0239             write(90,250) INVFLV(lblast(ip)),plast(1,ip),
0240      1           plast(2,ip),plast(3,ip),
0241      2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
0242      3           dplast(ip)
0243                   else
0244                      write(99,*) 'Unexpected perturbative particles'
0245                   endif
0246                endif
0247             else
0248 c     change format for large numbers:
0249                if(dplast(ip).gt.oneminus.and.dplast(ip).lt.oneplus) then
0250                    write(16,201) INVFLV(lblast(ip)), plast(1,ip),
0251      1           plast(2,ip),plast(3,ip),plast(4,ip),
0252      2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip)
0253                else
0254                   if(idpert.eq.1.or.idpert.eq.2) then
0255                      write(90,251) INVFLV(lblast(ip)), plast(1,ip),
0256      1           plast(2,ip),plast(3,ip),
0257      2           xlast(1,ip),xlast(2,ip),xlast(3,ip),xlast(4,ip),
0258      3           dplast(ip)
0259                   else
0260                      write(99,*) 'Unexpected perturbative particles'
0261                   endif
0262                endif
0263             endif
0264  1007    continue
0265          if(ioscar.eq.1) call hoscar
0266       endif
0267  190  format(3(i7),f10.4,5x,6(i4))
0268  191  format(3(i7),f10.4,5x,6(i4),5x,f7.4)
0269 clin-3/2009 improve the output accuracy of Pz
0270  200  format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,f8.2))
0271  201  format(I6,2(1x,f8.3),1x,f11.4,1x,f6.3,4(1x,e8.2))
0272  250  format(I5,2(1x,f8.3),1x,f10.3,2(1x,f7.1),1x,f8.2,1x,f7.2,1x,e10.4)
0273  251  format(I5,2(1x,f8.3),1x,f10.3,4(1x,e8.2),1x,e10.4)
0274 c
0275         return
0276         end
0277
0278 c=======================================================================
0279         SUBROUTINE decomp(px0,py0,pz0,xm0,i,itq1)
0280 c
0281         IMPLICIT DOUBLE PRECISION(D)
0282         DOUBLE PRECISION  enenew, pxnew, pynew, pznew
0283         DOUBLE PRECISION  de0, beta2, gam
0284         common /lor/ enenew, pxnew, pynew, pznew
0285 cc      SAVE /lor/
0286         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0287 cc      SAVE /HPARNT/
0288         common /decom/ptwo(2,5)
0289 cc      SAVE /decom/
0290         COMMON/RNDF77/NSEED
0291 cc      SAVE /RNDF77/
0292         COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
0293         common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd,
0294      1       psembd,tmaxembd,phidecomp
0295         SAVE
0296 c
0297         dcth=dble(RANART(NSEED))*2.d0-1.d0
0298         dPHI=dble(RANART(NSEED)*HIPR1(40))*2.d0
0299 clin-6/2009 Added if embedding a high-Pt quark pair after string melting:
0300         if(iembed.ge.1.and.iembed.le.4) then
0301 c     Decompose the parent high-Pt pion to q and qbar with an internal momentum
0302 c     parallel to the pion direction so that one parton has ~the same hight Pt
0303 c     and the other parton has a very soft Pt:
0304 c     Note: htop() decomposes a meson to q as it(1) followed by qbar as it(2):
0305            if(i.eq.(natt-2*nsembd).or.i.eq.(natt-2*nsembd-1)) then
0306               dcth=0.d0
0307               dphi=dble(phidecomp)
0308            endif
0309         endif
0310 c
0311         ds=dble(xm0)**2
0312         dpcm=dsqrt((ds-dble(ptwo(1,5)+ptwo(2,5))**2)
0313      1 *(ds-dble(ptwo(1,5)-ptwo(2,5))**2)/ds/4d0)
0314         dpz=dpcm*dcth
0315         dpx=dpcm*dsqrt(1.d0-dcth**2)*dcos(dphi)
0316         dpy=dpcm*dsqrt(1.d0-dcth**2)*dsin(dphi)
0317         de1=dsqrt(dble(ptwo(1,5))**2+dpcm**2)
0318         de2=dsqrt(dble(ptwo(2,5))**2+dpcm**2)
0319 c
0320       de0=dsqrt(dble(px0)**2+dble(py0)**2+dble(pz0)**2+dble(xm0)**2)
0321         dbex=dble(px0)/de0
0322         dbey=dble(py0)/de0
0323         dbez=dble(pz0)/de0
0324 c     boost the reference frame up by beta (pznew=gam(pz+beta e)):
0325       beta2 = dbex ** 2 + dbey ** 2 + dbez ** 2
0326       gam = 1.d0 / dsqrt(1.d0 - beta2)
0327       if(beta2.ge.0.9999999999999d0) then
0328          write(6,*) '1',dbex,dbey,dbez,beta2,gam
0329       endif
0330 c
0331       call lorenz(de1,dpx,dpy,dpz,-dbex,-dbey,-dbez)
0332         ptwo(1,1)=sngl(pxnew)
0333         ptwo(1,2)=sngl(pynew)
0334         ptwo(1,3)=sngl(pznew)
0335         ptwo(1,4)=sngl(enenew)
0336       call lorenz(de2,-dpx,-dpy,-dpz,-dbex,-dbey,-dbez)
0337         ptwo(2,1)=sngl(pxnew)
0338         ptwo(2,2)=sngl(pynew)
0339         ptwo(2,3)=sngl(pznew)
0340         ptwo(2,4)=sngl(enenew)
0341 c
0342       RETURN
0343       END
0344
0345 c=======================================================================
0346       SUBROUTINE HTOP
0347 c
0348       PARAMETER (MAXSTR=150001)
0349       PARAMETER (MAXPTN=400001)
0350       PARAMETER (MAXIDL=4001)
0351       DOUBLE PRECISION  GX0, GY0, GZ0, FT0, PX0, PY0, PZ0, E0, XMASS0
0352       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
0353      1     GXSGS,GYSGS,GZSGS,FTSGS
0354       dimension it(4)
0355       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
0356 cc      SAVE /HMAIN2/
0357       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
0358 cc      SAVE /HMAIN1/
0359       COMMON /PARA1/ MUL
0360 cc      SAVE /PARA1/
0361       COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0362      &     PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0363      &     XMASS0(MAXPTN), ITYP0(MAXPTN)
0364 cc      SAVE /prec1/
0365       COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
0366 cc      SAVE /ilist7/
0367       COMMON /ARPRC/ ITYPAR(MAXSTR),
0368      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
0369      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
0370      &     XMAR(MAXSTR)
0371 cc      SAVE /ARPRC/
0372       common /decom/ptwo(2,5)
0373 cc      SAVE /decom/
0374       COMMON/RNDF77/NSEED
0375 cc      SAVE /RNDF77/
0376       COMMON /NOPREC/ NNOZPC, ITYPN(MAXIDL),
0377      &     GXN(MAXIDL), GYN(MAXIDL), GZN(MAXIDL), FTN(MAXIDL),
0378      &     PXN(MAXIDL), PYN(MAXIDL), PZN(MAXIDL), EEN(MAXIDL),
0379      &     XMN(MAXIDL)
0380 cc      SAVE /NOPREC/
0381       COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
0382 cc      SAVE /HPARNT/
0383 c     7/20/01: use double precision
0384 c     otherwise sometimes beta>1 and gamma diverge in lorenz():
0385       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
0386      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
0387      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
0388      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
0389 cc      SAVE /SOFT/
0390       common/anim/nevent,isoft,isflag,izpc
0391 cc      SAVE /anim/
0392       DOUBLE PRECISION  vxp0,vyp0,vzp0
0393       common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
0394 cc      SAVE /precpa/
0395       common /para7/ ioscar,nsmbbbar,nsmmeson
0396       COMMON /AREVT/ IAEVT, IARUN, MISS
0397       common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
0398       SAVE
0399 c
0400         npar=0
0401         nnozpc=0
0402 clin-5b/2008 calculate the number of hadrons to be converted to q/qbar:
0403         if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
0404      1       then
0405            nsmbbbar=0
0406            nsmmeson=0
0407            do i=1,natt
0408               id=ITYPAR(i)
0409               idabs=iabs(id)
0410               i2=MOD(idabs/10,10)
0411 clin-9/2011 determine spectator nucleons consistently
0412 c              if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
0413 c     1             .ge.(HINT1(1)/2*0.99).and.
0414 c     2             .and.(id.eq.2112.or.id.eq.2212)) then
0415               if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt
0416      1             .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz)
0417      2                .or.PZAR(i).lt.(-PZTARG+epsiPz))
0418      3             .and.(id.eq.2112.or.id.eq.2212)) then
0419 c     spectator proj or targ nucleons without interactions, do not enter ZPC:
0420               elseif(idabs.gt.1000.and.i2.ne.0) then
0421 c     baryons to be converted to q/qbar:
0422                  nsmbbbar=nsmbbbar+1
0423               elseif((idabs.gt.100.and.idabs.lt.1000)
0424      1                .or.idabs.gt.10000) then
0425 c     mesons to be converted to q/qbar:
0426                  nsmmeson=nsmmeson+1
0427               endif
0428            enddo
0429
0430 clin-6/2009:
0431            if(ioscar.eq.2.or.ioscar.eq.3) then
0432               write(92,*) iaevt,miss,3*nsmbbbar+2*nsmmeson,
0433      1             nsmbbbar,nsmmeson,natt,natt-nsmbbbar-nsmmeson
0434            endif
0435 c           write(92,*) iaevt, 3*nsmbbbar+2*nsmmeson
0436 c           write(92,*) ' event#, total # of initial partons after string
0437 c     1 melting'
0438 c           write(92,*) 'String melting converts ',nsmbbbar, ' baryons &'
0439 c     1, nsmmeson, ' mesons'
0440 c           write(92,*) 'Total # of initial particles= ',natt
0441 c           write(92,*) 'Total # of initial particles (gamma,e,muon,...)
0442 c     1 not entering ZPC= ',natt-nsmbbbar-nsmmeson
0443         endif
0444 clin-5b/2008-over
0445         do 100 i=1,natt
0446            id=ITYPAR(i)
0447            idabs=iabs(id)
0448            i4=MOD(idabs/1000,10)
0449            i3=MOD(idabs/100,10)
0450            i2=MOD(idabs/10,10)
0451            i1=MOD(idabs,10)
0452            rnum=RANART(NSEED)
0453            ftime=0.197*PEAR(i)/(PXAR(i)**2+PYAR(i)**2+XMAR(i)**2)
0454            inozpc=0
0455            it(1)=0
0456            it(2)=0
0457            it(3)=0
0458            it(4)=0
0459 c
0460 clin-9/2011 determine spectator nucleons consistently
0461 c           if(PXAR(i).eq.0.and.PYAR(i).eq.0.and.PEAR(i)
0462 c     1 .ge.(HINT1(1)/2*0.99).and.((id.eq.2112).or.(id.eq.2212))) then
0463               if(abs(PXAR(i)).le.epsiPt.and.abs(PYAR(i)).le.epsiPt
0464      1             .and.(PZAR(i).gt.amax1(0.,PZPROJ-epsiPz)
0465      2                .or.PZAR(i).lt.(-PZTARG+epsiPz))
0466      3             .and.(id.eq.2112.or.id.eq.2212)) then
0467 c     spectator proj or targ nucleons without interactions, do not enter ZPC:
0468               inozpc=1
0469            elseif(idabs.gt.1000.and.i2.ne.0) then
0470 c     baryons:
0471               if(((i4.eq.1.or.i4.eq.2).and.i4.eq.i3)
0472      1 .or.(i4.eq.3.and.i3.eq.3)) then
0473                  if(i1.eq.2) then
0474                     if(rnum.le.(1./2.)) then
0475                        it(1)=i4
0476                        it(2)=i3*1000+i2*100+1
0477                     elseif(rnum.le.(2./3.)) then
0478                        it(1)=i4
0479                        it(2)=i3*1000+i2*100+3
0480                     else
0481                        it(1)=i2
0482                        it(2)=i4*1000+i3*100+3
0483                     endif
0484                  elseif(i1.eq.4) then
0485                     if(rnum.le.(2./3.)) then
0486                        it(1)=i4
0487                        it(2)=i3*1000+i2*100+3
0488                     else
0489                        it(1)=i2
0490                        it(2)=i4*1000+i3*100+3
0491                     endif
0492                  endif
0493               elseif(i4.eq.1.or.i4.eq.2) then
0494                  if(i1.eq.2) then
0495                     if(rnum.le.(1./2.)) then
0496                        it(1)=i2
0497                        it(2)=i4*1000+i3*100+1
0498                     elseif(rnum.le.(2./3.)) then
0499                        it(1)=i2
0500                        it(2)=i4*1000+i3*100+3
0501                     else
0502                        it(1)=i4
0503                        it(2)=i3*1000+i2*100+3
0504                     endif
0505                  elseif(i1.eq.4) then
0506                     if(rnum.le.(2./3.)) then
0507                        it(1)=i2
0508                        it(2)=i4*1000+i3*100+3
0509                     else
0510                        it(1)=i4
0511                        it(2)=i3*1000+i2*100+3
0512                     endif
0513                  endif
0514               elseif(i4.ge.3) then
0515                  it(1)=i4
0516                  if(i3.lt.i2) then
0517                     it(2)=i2*1000+i3*100+1
0518                  else
0519                     it(2)=i3*1000+i2*100+3
0520                  endif
0521               endif
0522 c       antibaryons:
0523               if(id.lt.0) then
0524                  it(1)=-it(1)
0525                  it(2)=-it(2)
0526               endif
0527 c     isoft=4or5 decompose diquark flavor it(2) to two quarks it(3)&(4):
0528               if(isoft.eq.4.or.isoft.eq.5) then
0529                  it(3)=MOD(it(2)/1000,10)
0530                  it(4)=MOD(it(2)/100,10)
0531               endif
0532
0533            elseif((idabs.gt.100.and.idabs.lt.1000)
0534      1 .or.idabs.gt.10000) then
0535 c     mesons:
0536               if(i3.eq.i2) then
0537                  if(i3.eq.1.or.i3.eq.2) then
0538                     if(rnum.le.0.5) then
0539                        it(1)=1
0540                        it(2)=-1
0541                     else
0542                        it(1)=2
0543                        it(2)=-2
0544                     endif
0545                  else
0546                     it(1)=i3
0547                     it(2)=-i3
0548                  endif
0549               else
0550                  if((isign(1,id)*(-1)**i3).eq.1) then
0551                     it(1)=i3
0552                     it(2)=-i2
0553                  else
0554                     it(1)=i2
0555                     it(2)=-i3
0556                  endif
0557               endif
0558            else
0559 c     save other particles (leptons and photons) outside of ZPC:
0560               inozpc=1
0561            endif
0562 c
0563            if(inozpc.eq.1) then
0564               NJSGS(i)=0
0565               nnozpc=nnozpc+1
0566               itypn(nnozpc)=ITYPAR(i)
0567               pxn(nnozpc)=PXAR(i)
0568               pyn(nnozpc)=PYAR(i)
0569               pzn(nnozpc)=PZAR(i)
0570               een(nnozpc)=PEAR(i)
0571               xmn(nnozpc)=XMAR(i)
0572               gxn(nnozpc)=GXAR(i)
0573               gyn(nnozpc)=GYAR(i)
0574               gzn(nnozpc)=GZAR(i)
0575               ftn(nnozpc)=FTAR(i)
0576            else
0577               NJSGS(i)=2
0578               ptwo(1,5)=ulmass(it(1))
0579               ptwo(2,5)=ulmass(it(2))
0580               call decomp(patt(i,1),patt(i,2),patt(i,3),XMAR(i),i,it(1))
0581               ipamax=2
0582               if((isoft.eq.4.or.isoft.eq.5)
0583      1 .and.iabs(it(2)).gt.1000) ipamax=1
0584               do 1001 ipar=1,ipamax
0585                  npar=npar+1
0586                  ityp0(npar)=it(ipar)
0587                  px0(npar)=dble(ptwo(ipar,1))
0588                  py0(npar)=dble(ptwo(ipar,2))
0589                  pz0(npar)=dble(ptwo(ipar,3))
0590                  e0(npar)=dble(ptwo(ipar,4))
0591                  xmass0(npar)=dble(ptwo(ipar,5))
0592                  gx0(npar)=dble(GXAR(i))
0593                  gy0(npar)=dble(GYAR(i))
0594                  gz0(npar)=dble(GZAR(i))
0595                  ft0(npar)=dble(ftime)
0596                  lstrg0(npar)=i
0597                  lpart0(npar)=ipar
0598                  vxp0(npar)=dble(patt(i,1)/patt(i,4))
0599                  vyp0(npar)=dble(patt(i,2)/patt(i,4))
0600                  vzp0(npar)=dble(patt(i,3)/patt(i,4))
0601  1001     continue
0602 cms 200      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
0603 cms 201      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
0604 c
0605               if((isoft.eq.4.or.isoft.eq.5)
0606      1 .and.iabs(it(2)).gt.1000) then
0607                  NJSGS(i)=3
0608                  xmdq=ptwo(2,5)
0609                  ptwo(1,5)=ulmass(it(3))
0610                  ptwo(2,5)=ulmass(it(4))
0611 c     8/19/02 avoid actual argument in common blocks of DECOMP:
0612 c                 call decomp(ptwo(2,1),ptwo(2,2),ptwo(2,3),xmdq)
0613              ptwox=ptwo(2,1)
0614              ptwoy=ptwo(2,2)
0615              ptwoz=ptwo(2,3)
0616              call decomp(ptwox,ptwoy,ptwoz,xmdq,i,it(1))
0617 c
0618                  do 1002 ipar=1,2
0619                     npar=npar+1
0620                     ityp0(npar)=it(ipar+2)
0621                     px0(npar)=dble(ptwo(ipar,1))
0622                     py0(npar)=dble(ptwo(ipar,2))
0623                     pz0(npar)=dble(ptwo(ipar,3))
0624                     e0(npar)=dble(ptwo(ipar,4))
0625                     xmass0(npar)=dble(ptwo(ipar,5))
0626                     gx0(npar)=dble(GXAR(i))
0627                     gy0(npar)=dble(GYAR(i))
0628                     gz0(npar)=dble(GZAR(i))
0629                     ft0(npar)=dble(ftime)
0630                     lstrg0(npar)=i
0631                     lpart0(npar)=ipar+1
0632                     vxp0(npar)=dble(patt(i,1)/patt(i,4))
0633                     vyp0(npar)=dble(patt(i,2)/patt(i,4))
0634                     vzp0(npar)=dble(patt(i,3)/patt(i,4))
0635  1002        continue
0636               endif
0637 c
0638            endif
0639  100        continue
0640       MUL=NPAR
0641 c
0642 clin-5b/2008:
0643       if((isoft.eq.4.or.isoft.eq.5).and.(ioscar.eq.2.or.ioscar.eq.3))
0644      1     then
0645          if((natt-nsmbbbar-nsmmeson).ne.nnozpc)
0646      1        write(92,*) 'Problem with the total # of initial particles
0647      2 (gamma,e,muon,...) not entering ZPC'
0648          if((3*nsmbbbar+2*nsmmeson).ne.npar)
0649      1        write(92,*) 'Problem with the total # of initial partons
0650      2 after string melting'
0651       endif
0652 c
0653       RETURN
0654       END
0655
0656 c=======================================================================
0657       SUBROUTINE PTOH
0658 c
0659       PARAMETER (MAXSTR=150001)
0660       DOUBLE PRECISION  gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
0661       DOUBLE PRECISION  gxp0,gyp0,gzp0,ft0fom,drlocl
0662       DOUBLE PRECISION  enenew, pxnew, pynew, pznew, beta2, gam
0663       DOUBLE PRECISION  ftavg0,gxavg0,gyavg0,gzavg0,bex,bey,bez
0664       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
0665      1     GXSGS,GYSGS,GZSGS,FTSGS
0666       DOUBLE PRECISION  xmdiag,px1,py1,pz1,e1,px2,py2,pz2,e2,
0667      1     px3,py3,pz3,e3,xmpair,etot
0668 clin-9/2012: improve precision for argument in sqrt():
0669       DOUBLE PRECISION  p1,p2,p3
0670       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
0671      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
0672 cc      SAVE /loclco/
0673       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
0674 cc      SAVE /HMAIN1/
0675       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
0676 cc      SAVE /HMAIN2/
0677       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
0678      &     K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
0679      &     PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
0680 cc      SAVE /HJJET2/
0681       COMMON /ARPRNT/ ARPAR1(100), IAPAR2(50), ARINT1(100), IAINT2(50)
0682 cc      SAVE /ARPRNT/
0683       COMMON /ARPRC/ ITYPAR(MAXSTR),
0684      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
0685      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
0686      &     XMAR(MAXSTR)
0687 cc      SAVE /ARPRC/
0688       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
0689      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
0690      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
0691      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
0692 cc      SAVE /SOFT/
0693       COMMON/RNDF77/NSEED
0694 cc      SAVE /RNDF77/
0695       common/anim/nevent,isoft,isflag,izpc
0696 cc      SAVE /anim/
0697       common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
0698 cc      SAVE /prtn23/
0699       common /nzpc/nattzp
0700 cc      SAVE /nzpc/
0701       common /lor/ enenew, pxnew, pynew, pznew
0702 cc      SAVE /lor/
0703       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0704 cc      SAVE /LUDAT1/
0705 clin 4/19/2006
0706       common /lastt/itimeh,bimp
0707       COMMON/HJGLBR/NELT,NINTHJ,NELP,NINP
0708       COMMON /AREVT/ IAEVT, IARUN, MISS
0709       common /para7/ ioscar,nsmbbbar,nsmmeson
0710 clin-5/2011
0711       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
0712 c
0713       dimension xmdiag(MAXSTR),indx(MAXSTR),ndiag(MAXSTR)
0714       SAVE
0715 c
0716       call coales
0717 c     obtain particle mass here without broadening by Breit-Wigner width:
0718       mstj24=MSTJ(24)
0719       MSTJ(24)=0
0720         nuudd=0
0721         npich=0
0722         nrhoch=0
0723       ppi0=1.
0724       prho0=0.
0725 c     determine hadron flavor except (pi0,rho0,eta,omega):
0726       DO 1001 ISG = 1, NSG
0727            if(NJSGS(ISG).ne.0) then
0728               NATT=NATT+1
0729               K1=K2SGS(ISG,1)
0730               k1abs=iabs(k1)
0731               PX1=PXSGS(ISG,1)
0732               PY1=PYSGS(ISG,1)
0733               PZ1=PZSGS(ISG,1)
0734               K2=K2SGS(ISG,2)
0735               k2abs=iabs(k2)
0736               PX2=PXSGS(ISG,2)
0737               PY2=PYSGS(ISG,2)
0738               PZ2=PZSGS(ISG,2)
0739 c     5/02/01 try lowest spin states as first choices,
0740 c     i.e. octet baryons and pseudoscalar mesons (ibs=2*baryonspin+1):
0741               e1=PESGS(ISG,1)
0742               e2=PESGS(ISG,2)
0743               xmpair=dsqrt((e1+e2)**2-(px1+px2)**2-(py1+py2)**2
0744      1 -(pz1+pz2)**2)
0745               ibs=2
0746               imspin=0
0747               if(k1.eq.-k2.and.iabs(k1).le.2.
0748      1           and.NJSGS(ISG).eq.2) then
0749                nuudd=nuudd+1
0750                xmdiag(nuudd)=xmpair
0751                ndiag(nuudd)=natt
0752             endif
0753               K3=0
0754               if((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3) then
0755                K3=K2SGS(ISG,3)
0756                k3abs=iabs(k3)
0757                PX3=PXSGS(ISG,3)
0758                PY3=PYSGS(ISG,3)
0759                PZ3=PZSGS(ISG,3)
0760                e3=PESGS(ISG,3)
0761                xmpair=dsqrt((e1+e2+e3)**2-(px1+px2+px3)**2
0762      1              -(py1+py2+py3)**2-(pz1+pz2+pz3)**2)
0763               endif
0764 c*****     isoft=3 baryon decomposition is different:
0765               if(isoft.eq.3.and.
0766      1           (k1abs.gt.1000.or.k2abs.gt.1000)) then
0767                if(k1abs.gt.1000) then
0768                   kdq=k1abs
0769                   kk=k2abs
0770                else
0771                   kdq=k2abs
0772                   kk=k1abs
0773                endif
0774                ki=MOD(kdq/1000,10)
0775                kj=MOD(kdq/100,10)
0776                if(MOD(kdq,10).eq.1) then
0777                   idqspn=0
0778                else
0779                   idqspn=1
0780                endif
0781 c
0782                if(kk.gt.ki) then
0783                   ktemp=kk
0784                   kk=kj
0785                   kj=ki
0786                   ki=ktemp
0787                elseif(kk.gt.kj) then
0788                   ktemp=kk
0789                   kk=kj
0790                   kj=ktemp
0791                endif
0792 c
0793                if(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
0794                   if(idqspn.eq.0) then
0795                      kf=1000*ki+100*kk+10*kj+ibs
0796                   else
0797                      kf=1000*ki+100*kj+10*kk+ibs
0798                   endif
0799                elseif(ki.eq.kj.and.ki.eq.kk) then
0800 c     can only be decuplet baryons:
0801                   kf=1000*ki+100*kj+10*kk+4
0802                else
0803                   kf=1000*ki+100*kj+10*kk+ibs
0804                endif
0805 c     form a decuplet baryon if the q+diquark mass is closer to its mass
0806 c     (and if the diquark has spin 1):
0807 cc     for now only include Delta, which is present in ART:
0808 cc                 if(idqspn.eq.1.and.MOD(kf,10).eq.2) then
0809                if(kf.eq.2112.or.kf.eq.2212) then
0810                   if(abs(sngl(xmpair)-ULMASS(kf)).gt.
0811      1                 abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
0812                endif
0813                if(k1.lt.0) kf=-kf
0814 clin-6/22/01 isoft=4or5 baryons:
0815               elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
0816      1              then
0817                if(k1abs.gt.k2abs) then
0818                   ki=k1abs
0819                   kk=k2abs
0820                else
0821                   ki=k2abs
0822                   kk=k1abs
0823                endif
0824                if(k3abs.gt.ki) then
0825                   kj=ki
0826                   ki=k3abs
0827                elseif(k3abs.lt.kk) then
0828                   kj=kk
0829                   kk=k3abs
0830                else
0831                   kj=k3abs
0832                endif
0833 c
0834                if(ki.eq.kj.and.ki.eq.kk) then
0835 c     can only be decuplet baryons (Delta-,++, Omega):
0836                   ibs=4
0837                   kf=1000*ki+100*kj+10*kk+ibs
0838                elseif(ki.ne.kj.and.ki.ne.kk.and.kj.ne.kk) then
0839 c     form Lambda or Sigma according to 3-quark mass,
0840 c     for now neglect decuplet (Sigma*0 etc) which is absent in ART:
0841                   ibs=2
0842                   kf1=1000*ki+100*kj+10*kk+ibs
0843                   kf2=1000*ki+100*kk+10*kj+ibs
0844                   kf=kf1
0845                   if(abs(sngl(xmpair)-ULMASS(kf1)).gt.
0846      1                 abs(sngl(xmpair)-ULMASS(kf2))) kf=kf2
0847                else
0848                   ibs=2
0849                   kf=1000*ki+100*kj+10*kk+ibs
0850 cc     for now only include Delta0,+ as decuplets, which are present in ART:
0851                   if(kf.eq.2112.or.kf.eq.2212) then
0852                      if(abs(sngl(xmpair)-ULMASS(kf)).gt.
0853      1                    abs(sngl(xmpair)-ULMASS(kf+2))) kf=kf+2
0854                   endif
0855                endif
0856                if(k1.lt.0) kf=-kf
0857 c*****     mesons:
0858               else
0859                if(k1abs.eq.k2abs) then
0860                   if(k1abs.le.2) then
0861 c     treat diagonal mesons later in the subroutine:
0862                      kf=0
0863                   elseif(k1abs.le.3) then
0864 c     do not form eta', only form phi from s-sbar, since no eta' in ART:
0865                      kf=333
0866                   else
0867                      kf=100*k1abs+10*k1abs+2*imspin+1
0868                   endif
0869                else
0870                   if(k1abs.gt.k2abs) then
0871                      kmax=k1abs
0872                      kmin=k2abs
0873                   elseif(k1abs.lt.k2abs) then
0874                      kmax=k2abs
0875                      kmin=k1abs
0876                   endif
0877                   kf=(100*kmax+10*kmin+2*imspin+1)
0878      1                 *isign(1,k1+k2)*(-1)**kmax
0879 c     form a vector meson if the q+qbar mass is closer to its mass:
0880                   if(MOD(iabs(kf),10).eq.1) then
0881                      if(abs(sngl(xmpair)-ULMASS(iabs(kf))).gt.
0882      1                    abs(sngl(xmpair)-ULMASS(iabs(kf)+2)))
0883      2                    kf=(iabs(kf)+2)*isign(1,kf)
0884                   endif
0885                endif
0886               endif
0887               ITYPAR(NATT)=kf
0888               KATT(NATT,1)=kf
0889             if(iabs(kf).eq.211) then
0890                npich=npich+1
0891             elseif(iabs(kf).eq.213) then
0892                nrhoch=nrhoch+1
0893             endif
0894            endif
0895 clin-7/2011-check charm hadron flavors:
0896 c           if(k1abs.eq.4.or.k2abs.eq.4) then
0897 c              if(k3.eq.0) then
0898 c                 write(99,*) iaevt,k1,k2,kf,xmpair,
0899 c     1                ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg
0900 c              else
0901 c                 write(99,*) iaevt,k1,k2,k3,kf,xmpair,
0902 c     1                ULMASS(iabs(kf)),ULMASS(iabs(kf)+2),isg
0903 c              endif
0904 c           endif
0905 clin-7/2011-end
0906  1001   CONTINUE
0907 c     assume Npi0=(Npi+ + Npi-)/2, Nrho0=(Nrho+ + Nrho-)/2 on the average:
0908         if(nuudd.ne.0) then
0909          ppi0=float(npich/2)/float(nuudd)
0910          prho0=float(nrhoch/2)/float(nuudd)
0911       endif
0912 c     determine diagonal mesons (pi0,rho0,eta and omega) from uubar/ddbar:
0913       npi0=0
0914       DO 1002 ISG = 1, NSG
0915          if(K2SGS(ISG,1).eq.-K2SGS(ISG,2)
0916      1        .and.iabs(K2SGS(ISG,1)).le.2.and.NJSGS(ISG).eq.2) then
0917             if(RANART(NSEED).le.ppi0) npi0=npi0+1
0918          endif
0919  1002 CONTINUE
0920 c
0921       if(nuudd.gt.1) then
0922          call index1(MAXSTR,nuudd,xmdiag,indx)
0923       else
0924          indx(1)=1
0925       end if
0926 c
0927       DO 1003 ix=1,nuudd
0928          iuudd=indx(ix)
0929          inatt=ndiag(iuudd)
0930          if(ix.le.npi0) then
0931             kf=111
0932          elseif(RANART(NSEED).le.(prho0/(1-ppi0+0.00001))) then
0933             kf=113
0934          else
0935 c     at T=150MeV, thermal weights for eta and omega(spin1) are about the same:
0936             if(RANART(NSEED).le.0.5) then
0937                kf=221
0938             else
0939                kf=223
0940             endif
0941          endif
0942          ITYPAR(inatt)=kf
0943          KATT(inatt,1)=kf
0944  1003 CONTINUE
0945 c  determine hadron formation time, position and momentum:
0946       inatt=0
0947 clin-6/2009 write out parton info after coalescence:
0948       if(ioscar.eq.3) then
0949          WRITE (85, 395) IAEVT, 3*nsmbbbar+2*nsmmeson,nsmbbbar,nsmmeson,
0950      1     bimp, NELP,NINP,NELT,NINTHJ,MISS
0951       endif
0952  395  format(4I8,f10.4,5I5)
0953 c
0954       DO 1006 ISG = 1, NSG
0955            if(NJSGS(ISG).ne.0) then
0956             inatt=inatt+1
0957               K1=K2SGS(ISG,1)
0958               k1abs=iabs(k1)
0959               PX1=PXSGS(ISG,1)
0960               PY1=PYSGS(ISG,1)
0961               PZ1=PZSGS(ISG,1)
0962               K2=K2SGS(ISG,2)
0963               k2abs=iabs(k2)
0964               PX2=PXSGS(ISG,2)
0965               PY2=PYSGS(ISG,2)
0966               PZ2=PZSGS(ISG,2)
0967               e1=PESGS(ISG,1)
0968               e2=PESGS(ISG,2)
0969 c
0970               if(NJSGS(ISG).eq.2) then
0971                PXAR(inatt)=sngl(px1+px2)
0972                PYAR(inatt)=sngl(py1+py2)
0973                PZAR(inatt)=sngl(pz1+pz2)
0974                PATT(inatt,1)=PXAR(inatt)
0975                PATT(inatt,2)=PYAR(inatt)
0976                PATT(inatt,3)=PZAR(inatt)
0977                etot=e1+e2
0978 clin-9/2012: improve precision for argument in sqrt():
0979                p1=px1+px2
0980                p2=py1+py2
0981                p3=pz1+pz2
0982 c
0983               elseif((isoft.eq.4.or.isoft.eq.5).and.NJSGS(ISG).eq.3)
0984      1              then
0985                PX3=PXSGS(ISG,3)
0986                PY3=PYSGS(ISG,3)
0987                PZ3=PZSGS(ISG,3)
0988                e3=PESGS(ISG,3)
0989                PXAR(inatt)=sngl(px1+px2+px3)
0990                PYAR(inatt)=sngl(py1+py2+py3)
0991                PZAR(inatt)=sngl(pz1+pz2+pz3)
0992                PATT(inatt,1)=PXAR(inatt)
0993                PATT(inatt,2)=PYAR(inatt)
0994                PATT(inatt,3)=PZAR(inatt)
0995                etot=e1+e2+e3
0996 clin-9/2012: improve precision for argument in sqrt():
0997                p1=px1+px2+px3
0998                p2=py1+py2+py3
0999                p3=pz1+pz2+pz3
1000 c
1001               endif
1002               XMAR(inatt)=ULMASS(ITYPAR(inatt))
1003 clin-5/2011-add finite width to resonances (rho,omega,eta,K*,phi,Delta) after formation:
1004               kf=KATT(inatt,1)
1005               if(kf.eq.113.or.abs(kf).eq.213.or.kf.eq.221.or.kf.eq.223
1006      1             .or.abs(kf).eq.313.or.abs(kf).eq.323.or.kf.eq.333
1007      2             .or.abs(kf).eq.1114.or.abs(kf).eq.2114
1008      3             .or.abs(kf).eq.2214.or.abs(kf).eq.2224) then
1009                  XMAR(inatt)=resmass(kf)
1010               endif
1011 c
1012               PEAR(inatt)=sqrt(PXAR(inatt)**2+PYAR(inatt)**2
1013      1           +PZAR(inatt)**2+XMAR(inatt)**2)
1014               PATT(inatt,4)=PEAR(inatt)
1015               EATT=EATT+PEAR(inatt)
1016             ipartn=NJSGS(ISG)
1017             DO 1004 i=1,ipartn
1018                ftp(i)=ftsgs(isg,i)
1019                gxp(i)=gxsgs(isg,i)
1020                gyp(i)=gysgs(isg,i)
1021                gzp(i)=gzsgs(isg,i)
1022                pxp(i)=pxsgs(isg,i)
1023                pyp(i)=pysgs(isg,i)
1024                pzp(i)=pzsgs(isg,i)
1025                pmp(i)=pmsgs(isg,i)
1026                pep(i)=pesgs(isg,i)
1027  1004       CONTINUE
1028             call locldr(ipartn,drlocl)
1029 c
1030             tau0=ARPAR1(1)
1031             ftavg0=ft0fom+dble(tau0)
1032             gxavg0=0d0
1033             gyavg0=0d0
1034             gzavg0=0d0
1035             DO 1005 i=1,ipartn
1036                gxavg0=gxavg0+gxp0(i)/ipartn
1037                gyavg0=gyavg0+gyp0(i)/ipartn
1038                gzavg0=gzavg0+gzp0(i)/ipartn
1039  1005       CONTINUE
1040 clin-9/2012: improve precision for argument in sqrt():
1041 c            bex=dble(PXAR(inatt))/etot
1042 c            bey=dble(PYAR(inatt))/etot
1043 c            bez=dble(PZAR(inatt))/etot
1044             bex=p1/etot
1045             bey=p2/etot
1046             bez=p3/etot
1047 c
1048             beta2 = bex ** 2 + bey ** 2 + bez ** 2
1049             gam = 1.d0 / dsqrt(1.d0 - beta2)
1050             if(beta2.ge.0.9999999999999d0) then
1051                write(6,*) '2',bex,bey,bez,beta2,gam
1052             endif
1053 c
1054             call lorenz(ftavg0,gxavg0,gyavg0,gzavg0,-bex,-bey,-bez)
1055               GXAR(inatt)=sngl(pxnew)
1056               GYAR(inatt)=sngl(pynew)
1057               GZAR(inatt)=sngl(pznew)
1058               FTAR(inatt)=sngl(enenew)
1059 clin 4/19/2006 write out parton info after coalescence:
1060               if(ioscar.eq.3) then
1061                  WRITE (85, 313) K2SGS(ISG,1),px1,py1,pz1,PMSGS(ISG,1),
1062      1                inatt,katt(inatt,1),xmar(inatt)
1063                  WRITE (85, 312) K2SGS(ISG,2),px2,py2,pz2,PMSGS(ISG,2),
1064      1                inatt,katt(inatt,1)
1065                  if(NJSGS(ISG).eq.3) WRITE (85, 312) K2SGS(ISG,3),
1066      1                px3,py3,pz3,PMSGS(ISG,3),inatt,katt(inatt,1)
1067               endif
1068  312       FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6)
1069 clin-5/02/2011
1070  313          FORMAT(I6,4(1X,F10.3),1X,I6,1X,I6,1X,F10.3)
1071 c
1072            endif
1073  1006   CONTINUE
1074 c     number of hadrons formed from partons inside ZPC:
1075       nattzp=natt
1076       MSTJ(24)=mstj24
1077 c
1078       RETURN
1079       END
1080
1081 c=======================================================================
1082 clin-5/2011-add finite width to resonances (rho,omega,eta,K*,phi,Delta) after formation:
1083       FUNCTION resmass(kf)
1084
1085       PARAMETER  (arho=0.775,aomega=0.783,aeta=0.548,aks=0.894,
1086      1     aphi=1.019,adelta=1.232)
1087       PARAMETER  (wrho=0.149,womega=0.00849,weta=1.30E-6,wks=0.0498,
1088      1     wphi=0.00426,wdelta=0.118)
1089       common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1090       COMMON/RNDF77/NSEED
1091       SAVE
1092
1093       if(kf.eq.113.or.abs(kf).eq.213) then
1094          amass=arho
1095          wid=wrho
1096       elseif(kf.eq.221) then
1097          amass=aeta
1098          wid=weta
1099       elseif(kf.eq.223) then
1100          amass=aomega
1101          wid=womega
1102       elseif(abs(kf).eq.313.or.abs(kf).eq.323) then
1103          amass=aks
1104          wid=wks
1105       elseif(kf.eq.333) then
1106          amass=aphi
1107          wid=wphi
1108       elseif(abs(kf).eq.1114.or.abs(kf).eq.2114
1109      1        .or.abs(kf).eq.2214.or.abs(kf).eq.2224) then
1110          amass=adelta
1111          wid=wdelta
1112       endif
1113       dmin=amass-2*wid
1114       dmax=amass+2*wid
1115 c     Delta mass needs to be big enough to decay to N+pi:
1116       if(amass.eq.adelta) dmin=1.078
1117 c
1118       FM=1.
1119       NTRY1=0
1120  10   DM = RANART(NSEED) * (DMAX-DMIN) + DMIN
1121       NTRY1=NTRY1+1
1122       fmass=(amass*wid)**2/((DM**2-amass**2)**2+(amass*wid)**2)
1123 check      write (99,*) ntry1,kf,amass,wid,fmass,DM
1124       IF((RANART(NSEED) .GT. FMASS/FM).AND. (NTRY1.LE.10)) GOTO 10
1125 c
1126       resmass=DM
1127
1128       RETURN
1129       END
1130
1131 c=======================================================================
1132       SUBROUTINE coales
1133
1134       PARAMETER (MAXSTR=150001)
1135       IMPLICIT DOUBLE PRECISION(D)
1136       DOUBLE PRECISION  gxp,gyp,gzp,ftp,pxp,pyp,pzp,pep,pmp
1137       DIMENSION IOVER(MAXSTR),dp1(2:3),dr1(2:3)
1138       DOUBLE PRECISION  PXSGS,PYSGS,PZSGS,PESGS,PMSGS,
1139      1     GXSGS,GYSGS,GZSGS,FTSGS
1140       double precision  dpcoal,drcoal,ecritl
1141       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1142      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1143      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1144      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1145 cc      SAVE /SOFT/
1146       common /coal/dpcoal,drcoal,ecritl
1147 cc      SAVE /coal/
1148       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1149      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1150 cc      SAVE /loclco/
1151       COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
1152      &     K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
1153      &     PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
1154 cc      SAVE /HJJET2/
1155       SAVE
1156 c
1157       do 1001 ISG=1, NSG
1158          IOVER(ISG)=0
1159  1001 continue
1160 C1     meson q coalesce with all available qbar:
1161       do 150 ISG=1,NSG
1162          if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 150
1163 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1164          if(K2SGS(ISG,1).lt.0) then
1165             write(6,*) 'Antiquark appears in quark loop; stop'
1166             stop
1167          endif
1168 c
1169          do 1002 j=1,2
1170             ftp(j)=ftsgs(isg,j)
1171             gxp(j)=gxsgs(isg,j)
1172             gyp(j)=gysgs(isg,j)
1173             gzp(j)=gzsgs(isg,j)
1174             pxp(j)=pxsgs(isg,j)
1175             pyp(j)=pysgs(isg,j)
1176             pzp(j)=pzsgs(isg,j)
1177             pmp(j)=pmsgs(isg,j)
1178             pep(j)=pesgs(isg,j)
1179  1002    continue
1180          call locldr(2,drlocl)
1181          dr0=drlocl
1182 c     dp0^2 defined as (p1+p2)^2-(m1+m2)^2:
1183          dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1184      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1185 c
1186          do 120 JSG=1,NSG
1187 c     skip default or unavailable antiquarks:
1188             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 120
1189             if(NJSGS(JSG).eq.2) then
1190                ipmin=2
1191                ipmax=2
1192             elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).lt.0) then
1193                ipmin=1
1194                ipmax=3
1195             else
1196                goto 120
1197             endif
1198             do 100 ip=ipmin,ipmax
1199                dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1200      1              -pxp(1)*pxsgs(jsg,ip)
1201      2              -pyp(1)*pysgs(jsg,ip)
1202      3              -pzp(1)*pzsgs(jsg,ip)
1203      4              -pmp(1)*pmsgs(jsg,ip)))
1204 c     skip if outside of momentum radius:
1205                if(dplocl.gt.dpcoal) goto 120
1206                ftp(2)=ftsgs(jsg,ip)
1207                gxp(2)=gxsgs(jsg,ip)
1208                gyp(2)=gysgs(jsg,ip)
1209                gzp(2)=gzsgs(jsg,ip)
1210                pxp(2)=pxsgs(jsg,ip)
1211                pyp(2)=pysgs(jsg,ip)
1212                pzp(2)=pzsgs(jsg,ip)
1213                pmp(2)=pmsgs(jsg,ip)
1214                pep(2)=pesgs(jsg,ip)
1215                call locldr(2,drlocl)
1216 c     skip if outside of spatial radius:
1217                if(drlocl.gt.drcoal) goto 120
1218 c     q_isg coalesces with qbar_jsg:
1219                if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1220      1              .or.(drlocl.lt.dr0)) then
1221                   dp0=dplocl
1222                   dr0=drlocl
1223                   call exchge(isg,2,jsg,ip)
1224                endif
1225  100        continue
1226  120     continue
1227          if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1228  150  continue
1229 c
1230 C2     meson qbar coalesce with all available q:
1231       do 250 ISG=1,NSG
1232          if(NJSGS(ISG).ne.2.or.IOVER(ISG).eq.1) goto 250
1233 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1234          do 1003 j=1,2
1235             ftp(j)=ftsgs(isg,j)
1236             gxp(j)=gxsgs(isg,j)
1237             gyp(j)=gysgs(isg,j)
1238             gzp(j)=gzsgs(isg,j)
1239             pxp(j)=pxsgs(isg,j)
1240             pyp(j)=pysgs(isg,j)
1241             pzp(j)=pzsgs(isg,j)
1242             pmp(j)=pmsgs(isg,j)
1243             pep(j)=pesgs(isg,j)
1244  1003    continue
1245          call locldr(2,drlocl)
1246          dr0=drlocl
1247          dp0=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1248      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1249 c
1250          do 220 JSG=1,NSG
1251             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 220
1252             if(NJSGS(JSG).eq.2) then
1253                ipmin=1
1254                ipmax=1
1255             elseif(NJSGS(JSG).eq.3.and.K2SGS(JSG,1).gt.0) then
1256                ipmin=1
1257                ipmax=3
1258             else
1259                goto 220
1260             endif
1261             do 200 ip=ipmin,ipmax
1262                dplocl=dsqrt(2*(pep(2)*pesgs(jsg,ip)
1263      1              -pxp(2)*pxsgs(jsg,ip)
1264      2              -pyp(2)*pysgs(jsg,ip)
1265      3              -pzp(2)*pzsgs(jsg,ip)
1266      4              -pmp(2)*pmsgs(jsg,ip)))
1267 c     skip if outside of momentum radius:
1268                if(dplocl.gt.dpcoal) goto 220
1269                ftp(1)=ftsgs(jsg,ip)
1270                gxp(1)=gxsgs(jsg,ip)
1271                gyp(1)=gysgs(jsg,ip)
1272                gzp(1)=gzsgs(jsg,ip)
1273                pxp(1)=pxsgs(jsg,ip)
1274                pyp(1)=pysgs(jsg,ip)
1275                pzp(1)=pzsgs(jsg,ip)
1276                pmp(1)=pmsgs(jsg,ip)
1277                pep(1)=pesgs(jsg,ip)
1278                call locldr(2,drlocl)
1279 c     skip if outside of spatial radius:
1280                if(drlocl.gt.drcoal) goto 220
1281 c     qbar_isg coalesces with q_jsg:
1282                if((dp0.gt.dpcoal.or.dr0.gt.drcoal)
1283      1              .or.(drlocl.lt.dr0)) then
1284                   dp0=dplocl
1285                   dr0=drlocl
1286                   call exchge(isg,1,jsg,ip)
1287                endif
1288  200        continue
1289  220     continue
1290          if(dp0.le.dpcoal.and.dr0.le.drcoal) IOVER(ISG)=1
1291  250  continue
1292 c
1293 C3     baryon q (antibaryon qbar) coalesce with all available q (qbar):
1294       do 350 ISG=1,NSG
1295          if(NJSGS(ISG).ne.3.or.IOVER(ISG).eq.1) goto 350
1296          ibaryn=K2SGS(ISG,1)
1297 C     DETERMINE CURRENT RELATIVE DISTANCE AND MOMENTUM:
1298          do 1004 j=1,2
1299             ftp(j)=ftsgs(isg,j)
1300             gxp(j)=gxsgs(isg,j)
1301             gyp(j)=gysgs(isg,j)
1302             gzp(j)=gzsgs(isg,j)
1303             pxp(j)=pxsgs(isg,j)
1304             pyp(j)=pysgs(isg,j)
1305             pzp(j)=pzsgs(isg,j)
1306             pmp(j)=pmsgs(isg,j)
1307             pep(j)=pesgs(isg,j)
1308  1004    continue
1309          call locldr(2,drlocl)
1310          dr1(2)=drlocl
1311          dp1(2)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1312      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1313 c
1314          ftp(2)=ftsgs(isg,3)
1315          gxp(2)=gxsgs(isg,3)
1316          gyp(2)=gysgs(isg,3)
1317          gzp(2)=gzsgs(isg,3)
1318          pxp(2)=pxsgs(isg,3)
1319          pyp(2)=pysgs(isg,3)
1320          pzp(2)=pzsgs(isg,3)
1321          pmp(2)=pmsgs(isg,3)
1322          pep(2)=pesgs(isg,3)
1323          call locldr(2,drlocl)
1324          dr1(3)=drlocl
1325          dp1(3)=dsqrt(2*(pep(1)*pep(2)-pxp(1)*pxp(2)
1326      &        -pyp(1)*pyp(2)-pzp(1)*pzp(2)-pmp(1)*pmp(2)))
1327 c
1328          do 320 JSG=1,NSG
1329             if(JSG.eq.ISG.or.IOVER(JSG).eq.1) goto 320
1330             if(NJSGS(JSG).eq.2) then
1331                if(ibaryn.gt.0) then
1332                   ipmin=1
1333                else
1334                   ipmin=2
1335                endif
1336                ipmax=ipmin
1337             elseif(NJSGS(JSG).eq.3.and.
1338      1              (ibaryn*K2SGS(JSG,1)).gt.0) then
1339                ipmin=1
1340                ipmax=3
1341             else
1342                goto 320
1343             endif
1344             do 300 ip=ipmin,ipmax
1345                dplocl=dsqrt(2*(pep(1)*pesgs(jsg,ip)
1346      1              -pxp(1)*pxsgs(jsg,ip)
1347      2              -pyp(1)*pysgs(jsg,ip)
1348      3              -pzp(1)*pzsgs(jsg,ip)
1349      4              -pmp(1)*pmsgs(jsg,ip)))
1350 c     skip if outside of momentum radius:
1351                if(dplocl.gt.dpcoal) goto 320
1352                ftp(2)=ftsgs(jsg,ip)
1353                gxp(2)=gxsgs(jsg,ip)
1354                gyp(2)=gysgs(jsg,ip)
1355                gzp(2)=gzsgs(jsg,ip)
1356                pxp(2)=pxsgs(jsg,ip)
1357                pyp(2)=pysgs(jsg,ip)
1358                pzp(2)=pzsgs(jsg,ip)
1359                pmp(2)=pmsgs(jsg,ip)
1360                pep(2)=pesgs(jsg,ip)
1361                call locldr(2,drlocl)
1362 c     skip if outside of spatial radius:
1363                if(drlocl.gt.drcoal) goto 320
1364 c     q_isg may coalesce with q_jsg for a baryon:
1365                ipi=0
1366                if(dp1(2).gt.dpcoal.or.dr1(2).gt.drcoal) then
1367                   ipi=2
1368                   if((dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal)
1369      1                 .and.dr1(3).gt.dr1(2)) ipi=3
1370                elseif(dp1(3).gt.dpcoal.or.dr1(3).gt.drcoal) then
1371                   ipi=3
1372                elseif(dr1(2).lt.dr1(3)) then
1373                   if(drlocl.lt.dr1(3)) ipi=3
1374                elseif(dr1(3).le.dr1(2)) then
1375                   if(drlocl.lt.dr1(2)) ipi=2
1376                endif
1377                if(ipi.ne.0) then
1378                   dp1(ipi)=dplocl
1379                   dr1(ipi)=drlocl
1380                   call exchge(isg,ipi,jsg,ip)
1381                endif
1382  300        continue
1383  320     continue
1384          if(dp1(2).le.dpcoal.and.dr1(2).le.drcoal
1385      1        .and.dp1(3).le.dpcoal.and.dr1(3).le.drcoal)
1386      2        IOVER(ISG)=1
1387  350  continue
1388 c
1389       RETURN
1390       END
1391
1392 c=======================================================================
1393       SUBROUTINE exchge(isg,ipi,jsg,ipj)
1394 c
1395       implicit double precision  (a-h, o-z)
1396       PARAMETER (MAXSTR=150001)
1397       COMMON/SOFT/PXSGS(MAXSTR,3),PYSGS(MAXSTR,3),PZSGS(MAXSTR,3),
1398      &     PESGS(MAXSTR,3),PMSGS(MAXSTR,3),GXSGS(MAXSTR,3),
1399      &     GYSGS(MAXSTR,3),GZSGS(MAXSTR,3),FTSGS(MAXSTR,3),
1400      &     K1SGS(MAXSTR,3),K2SGS(MAXSTR,3),NJSGS(MAXSTR)
1401 cc      SAVE /SOFT/
1402       SAVE
1403 c
1404       k1=K1SGS(isg,ipi)
1405       k2=K2SGS(isg,ipi)
1406       px=PXSGS(isg,ipi)
1407       py=PYSGS(isg,ipi)
1408       pz=PZSGS(isg,ipi)
1409       pe=PESGS(isg,ipi)
1410       pm=PMSGS(isg,ipi)
1411       gx=GXSGS(isg,ipi)
1412       gy=GYSGS(isg,ipi)
1413       gz=GZSGS(isg,ipi)
1414       ft=FTSGS(isg,ipi)
1415       K1SGS(isg,ipi)=K1SGS(jsg,ipj)
1416       K2SGS(isg,ipi)=K2SGS(jsg,ipj)
1417       PXSGS(isg,ipi)=PXSGS(jsg,ipj)
1418       PYSGS(isg,ipi)=PYSGS(jsg,ipj)
1419       PZSGS(isg,ipi)=PZSGS(jsg,ipj)
1420       PESGS(isg,ipi)=PESGS(jsg,ipj)
1421       PMSGS(isg,ipi)=PMSGS(jsg,ipj)
1422       GXSGS(isg,ipi)=GXSGS(jsg,ipj)
1423       GYSGS(isg,ipi)=GYSGS(jsg,ipj)
1424       GZSGS(isg,ipi)=GZSGS(jsg,ipj)
1425       FTSGS(isg,ipi)=FTSGS(jsg,ipj)
1426       K1SGS(jsg,ipj)=k1
1427       K2SGS(jsg,ipj)=k2
1428       PXSGS(jsg,ipj)=px
1429       PYSGS(jsg,ipj)=py
1430       PZSGS(jsg,ipj)=pz
1431       PESGS(jsg,ipj)=pe
1432       PMSGS(jsg,ipj)=pm
1433       GXSGS(jsg,ipj)=gx
1434       GYSGS(jsg,ipj)=gy
1435       GZSGS(jsg,ipj)=gz
1436       FTSGS(jsg,ipj)=ft
1437 c
1438       RETURN
1439       END
1440
1441 c=======================================================================
1442       SUBROUTINE locldr(icall,drlocl)
1443 c
1444       implicit double precision (a-h, o-z)
1445       dimension ftp0(3),pxp0(3),pyp0(3),pzp0(3),pep0(3)
1446       common /loclco/gxp(3),gyp(3),gzp(3),ftp(3),
1447      1     pxp(3),pyp(3),pzp(3),pep(3),pmp(3)
1448 cc      SAVE /loclco/
1449       common /prtn23/ gxp0(3),gyp0(3),gzp0(3),ft0fom
1450 cc      SAVE /prtn23/
1451       common /lor/ enenew, pxnew, pynew, pznew
1452 cc      SAVE /lor/
1453       SAVE
1454 c     for 2-body kinematics:
1455       if(icall.eq.2) then
1456          etot=pep(1)+pep(2)
1457          bex=(pxp(1)+pxp(2))/etot
1458          bey=(pyp(1)+pyp(2))/etot
1459          bez=(pzp(1)+pzp(2))/etot
1460 c     boost the reference frame down by beta to get to the pair rest frame:
1461          do 1001 j=1,2
1462             beta2 = bex ** 2 + bey ** 2 + bez ** 2
1463             gam = 1.d0 / dsqrt(1.d0 - beta2)
1464             if(beta2.ge.0.9999999999999d0) then
1465                write(6,*) '4',pxp(1),pxp(2),pyp(1),pyp(2),
1466      1              pzp(1),pzp(2),pep(1),pep(2),pmp(1),pmp(2),
1467      2          dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2+pmp(1)**2)/pep(1),
1468      3          dsqrt(pxp(1)**2+pyp(1)**2+pzp(1)**2)/pep(1)
1469                write(6,*) '4a',pxp(1)+pxp(2),pyp(1)+pyp(2),
1470      1              pzp(1)+pzp(2),etot
1471                write(6,*) '4b',bex,bey,bez,beta2,gam
1472             endif
1473 c
1474             call lorenz(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
1475             gxp0(j)=pxnew
1476             gyp0(j)=pynew
1477             gzp0(j)=pznew
1478             ftp0(j)=enenew
1479             call lorenz(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
1480             pxp0(j)=pxnew
1481             pyp0(j)=pynew
1482             pzp0(j)=pznew
1483             pep0(j)=enenew
1484  1001    continue
1485 c
1486          if(ftp0(1).ge.ftp0(2)) then
1487             ilate=1
1488             iearly=2
1489          else
1490             ilate=2
1491             iearly=1
1492          endif
1493          ft0fom=ftp0(ilate)
1494 c
1495          dt0=ftp0(ilate)-ftp0(iearly)
1496          gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1497          gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1498          gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1499          drlocl=dsqrt((gxp0(ilate)-gxp0(iearly))**2
1500      1        +(gyp0(ilate)-gyp0(iearly))**2
1501      2        +(gzp0(ilate)-gzp0(iearly))**2)
1502 c     for 3-body kinematics, used for baryons formation:
1503       elseif(icall.eq.3) then
1504          etot=pep(1)+pep(2)+pep(3)
1505          bex=(pxp(1)+pxp(2)+pxp(3))/etot
1506          bey=(pyp(1)+pyp(2)+pyp(3))/etot
1507          bez=(pzp(1)+pzp(2)+pzp(3))/etot
1508          beta2 = bex ** 2 + bey ** 2 + bez ** 2
1509          gam = 1.d0 / dsqrt(1.d0 - beta2)
1510          if(beta2.ge.0.9999999999999d0) then
1511             write(6,*) '5',bex,bey,bez,beta2,gam
1512          endif
1513 c     boost the reference frame down by beta to get to the 3-parton rest frame:
1514          do 1002 j=1,3
1515             call lorenz(ftp(j),gxp(j),gyp(j),gzp(j),bex,bey,bez)
1516             gxp0(j)=pxnew
1517             gyp0(j)=pynew
1518             gzp0(j)=pznew
1519             ftp0(j)=enenew
1520             call lorenz(pep(j),pxp(j),pyp(j),pzp(j),bex,bey,bez)
1521             pxp0(j)=pxnew
1522             pyp0(j)=pynew
1523             pzp0(j)=pznew
1524             pep0(j)=enenew
1525  1002    continue
1526 c
1527          if(ftp0(1).gt.ftp0(2)) then
1528             ilate=1
1529             if(ftp0(3).gt.ftp0(1)) ilate=3
1530          else
1531             ilate=2
1532             if(ftp0(3).ge.ftp0(2)) ilate=3
1533          endif
1534          ft0fom=ftp0(ilate)
1535 c
1536          if(ilate.eq.1) then
1537             imin=2
1538             imax=3
1539             istep=1
1540          elseif(ilate.eq.2) then
1541             imin=1
1542             imax=3
1543             istep=2
1544          elseif(ilate.eq.3) then
1545             imin=1
1546             imax=2
1547             istep=1
1548          endif
1549 c
1550          do 1003 iearly=imin,imax,istep
1551             dt0=ftp0(ilate)-ftp0(iearly)
1552             gxp0(iearly)=gxp0(iearly)+pxp0(iearly)/pep0(iearly)*dt0
1553             gyp0(iearly)=gyp0(iearly)+pyp0(iearly)/pep0(iearly)*dt0
1554             gzp0(iearly)=gzp0(iearly)+pzp0(iearly)/pep0(iearly)*dt0
1555  1003    continue
1556       endif
1557 c
1558       RETURN
1559       END
1560
1561 c=======================================================================
1562         subroutine hoscar
1563 c
1564         parameter (MAXSTR=150001,AMN=0.939457,AMP=0.93828)
1565         character*8 code, reffra, FRAME
1566         character*25 amptvn
1567         common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
1568 cc      SAVE /snn/
1569         common /lastt/itimeh,bimp
1570 cc      SAVE /lastt/
1571         COMMON/hbt/lblast(MAXSTR),xlast(4,MAXSTR),plast(4,MAXSTR),nlast
1572 cc      SAVE /hbt/
1573         common/oscar1/iap,izp,iat,izt
1574 cc      SAVE /oscar1/
1575         common/oscar2/FRAME,amptvn
1576 cc      SAVE /oscar2/
1577         SAVE
1578         data nff/0/
1579 c
1580 c       file header
1581         if(nff.eq.0) then
1582            write (19, 101) 'OSCAR1997A'
1583            write (19, 111) 'final_id_p_x'
1584            code = 'AMPT'
1585            if(FRAME.eq.'CMS') then
1586               reffra = 'nncm'
1587               xmp=(amp*izp+amn*(iap-izp))/iap
1588               xmt=(amp*izt+amn*(iat-izt))/iat
1589               ebeam=(efrm**2-xmp**2-xmt**2)/2./xmt
1590            elseif(FRAME.eq.'LAB') then
1591               reffra = 'lab'
1592               ebeam=efrm
1593            else
1594               reffra = 'unknown'
1595               ebeam=0.
1596            endif
1597            ntestp = 1
1598            write (19, 102) code, amptvn, iap, izp, iat, izt,
1599      &        reffra, ebeam, ntestp
1600            nff = 1
1601            ievent = 1
1602            phi = 0.
1603            if(FRAME.eq.'CMS') write(19,112) efrm
1604         endif
1605 c       comment
1606 c       event header
1607         write (19, 103) ievent, nlast, bimp, phi
1608 c       particles
1609         do 99 i = 1, nlast
1610            ene=sqrt(plast(1,i)**2+plast(2,i)**2+plast(3,i)**2
1611      1          +plast(4,i)**2)
1612            write (19, 104) i, INVFLV(lblast(i)), plast(1,i),
1613      1          plast(2,i),plast(3,i),ene,plast(4,i),
1614      2          xlast(1,i),xlast(2,i),xlast(3,i),xlast(4,i)
1615  99     continue
1616         ievent = ievent + 1
1617  101        format (a10)
1618  111        format (a12)
1619  102        format (a4,1x,a20,1x,'(', i3, ',', i3, ')+(', i3, ',',
1620      &           i3, ')', 2x, a4, 2x, e10.4, 2x, i8)
1621  103        format (i10, 2x, i10, 2x, f8.3, 2x, f8.3)
1622  104        format (i10, 2x, i10, 2x, 9(e12.6, 2x))
1623  112        format ('# Center-of-mass energy/nucleon-pair is',
1624      & f12.3,'GeV')
1625 c
1626         return
1627         end
1628
1629 c=======================================================================
1630         subroutine getnp
1631
1632         PARAMETER (MAXSTR=150001)
1633         COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
1634 cc      SAVE /HMAIN1/
1635         COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
1636 cc      SAVE /HMAIN2/
1637         COMMON /HPARNT/HIPR1(100), IHPR2(50), HINT1(100), IHNT2(50)
1638 cc      SAVE /HPARNT/
1639         common/snn/efrm,npart1,npart2,epsiPz,epsiPt,PZPROJ,PZTARG
1640 cc      SAVE /snn/
1641         SAVE
1642
1643         if(NATT.eq.0) then
1644            npart1=0
1645            npart2=0
1646            return
1647         endif
1648 c
1649         PZPROJ=SQRT(HINT1(6)**2-HINT1(8)**2)
1650         PZTARG=SQRT(HINT1(7)**2-HINT1(9)**2)
1651         epsiPz=0.01
1652 clin-9/2011-add Pt tolerance in determining spectator nucleons
1653 c     (affect string melting runs when LAB frame is used):
1654         epsiPt=1e-6
1655 c
1656         nspec1=0
1657         nspec2=0
1658         DO 1000 I = 1, NATT
1659 clin-9/2011 determine spectator nucleons consistently
1660 c           if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212)
1661 c     1          .and.PATT(I, 1).eq.0.and.PATT(I, 2).eq.0) then
1662            if((KATT(I,1).eq.2112.or.KATT(I,1).eq.2212)
1663      1          .and.abs(PATT(I, 1)).le.epsiPt
1664      2          .and.abs(PATT(I, 2)).le.epsiPt) then
1665               if(PATT(I, 3).gt.amax1(0.,PZPROJ-epsiPz)) then
1666                  nspec1=nspec1+1
1667               elseif(PATT(I, 3).lt.(-PZTARG+epsiPz)) then
1668                  nspec2=nspec2+1
1669               endif
1670            endif
1671  1000   CONTINUE
1672         npart1=IHNT2(1)-nspec1
1673         npart2=IHNT2(3)-nspec2
1674
1675         return
1676         end
1677
1678 c=======================================================================
1679 c     2/18/03 use PYTHIA to decay eta,rho,omega,k*,phi and Delta
1680 c     4/2012 added pi0 decay flag:
1681 c       ipion=0: resonance or pi0 in lb(i1); >0: pi0 in lpion(ipion).
1682         subroutine resdec(i1,nt,nnn,wid,idecay,ipion)
1683
1684         PARAMETER (hbarc=0.19733)
1685         PARAMETER (AK0=0.498,APICH=0.140,API0=0.135,AN=0.940,ADDM=0.02)
1686         PARAMETER (MAXSTR=150001, MAXR=1)
1687         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
1688      &       IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
1689 cc      SAVE /INPUT2/
1690         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
1691 cc      SAVE /LUJETS/
1692         COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
1693 cc      SAVE /LUDAT1/
1694         COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
1695 cc      SAVE /LUDAT2/
1696         COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
1697 cc      SAVE /LUDAT3/
1698         COMMON /CC/ E(MAXSTR)
1699 cc      SAVE /CC/
1700         COMMON /EE/ ID(MAXSTR),LB(MAXSTR)
1701 cc      SAVE /EE/
1702         COMMON   /PA/RPION(3,MAXSTR,MAXR)
1703 cc      SAVE /PA/
1704         COMMON   /PB/PPION(3,MAXSTR,MAXR)
1705 cc      SAVE /PB/
1706         COMMON   /PC/EPION(MAXSTR,MAXR)
1707 cc      SAVE /PC/
1708         COMMON   /PD/LPION(MAXSTR,MAXR)
1709 cc      SAVE /PD/
1710         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
1711 cc      SAVE /input1/
1712         common/resdcy/NSAV,iksdcy
1713 cc      SAVE /resdcy/
1714         common/leadng/lb1,px1,py1,pz1,em1,e1,xfnl,yfnl,zfnl,tfnl,
1715      1       px1n,py1n,pz1n,dp1n
1716 cc      SAVE /leadng/
1717         EXTERNAL IARFLV, INVFLV
1718         COMMON/tdecay/tfdcy(MAXSTR),tfdpi(MAXSTR,MAXR),tft(MAXSTR)
1719 cc      SAVE /tdecay/
1720         COMMON/RNDF77/NSEED
1721         COMMON /dpert/dpertt(MAXSTR,MAXR),dpertp(MAXSTR),dplast(MAXSTR),
1722      1       dpdcy(MAXSTR),dpdpi(MAXSTR,MAXR),dpt(MAXSTR, MAXR),
1723      2       dpp1(MAXSTR,MAXR),dppion(MAXSTR,MAXR)
1724 cc      SAVE /RNDF77/
1725         common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy
1726         SAVE
1727         irun=idecay
1728 clin-4/2012 for option of pi0 decay:
1729         if(nt.eq.ntmax.and.ipi0dcy.eq.1
1730      &       .and.((lb1.eq.4.and.ipion.eq.0).or.ipion.ge.1)) then
1731            kf=111
1732 c        if(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1733         elseif(lb1.eq.0.or.lb1.eq.25.or.lb1.eq.26.or.lb1.eq.27
1734      &       .or.lb1.eq.28.or.lb1.eq.29.or.iabs(lb1).eq.30
1735      &       .or.lb1.eq.24.or.(iabs(lb1).ge.6.and.iabs(lb1).le.9)
1736      &       .or.iabs(lb1).eq.16) then
1737            kf=INVFLV(lb1)
1738         else
1739            return
1740         endif
1741 c
1742         IP=1
1743 c     label as undecayed and the only particle in the record:
1744         N=1
1745         K(IP,1)=1
1746         K(IP,3)=0
1747         K(IP,4)=0
1748         K(IP,5)=0
1749 c
1750         K(IP,2)=kf
1751 clin-4/2012 for option of pi0 decay:
1752         if(ipion.eq.0) then
1753 c
1754         P(IP,1)=px1
1755         P(IP,2)=py1
1756         P(IP,3)=pz1
1757 c        em1a=em1
1758 c     eta or omega in ART may be below or too close to (pi+pi-pi0) mass,
1759 c     causing LUDECY error,thus increase their mass ADDM above this thresh,
1760 c     noting that rho (m=0.281) too close to 2pi thrshold fails to decay:
1761         if((lb1.eq.0.or.lb1.eq.28).and.em1.lt.(2*APICH+API0+ADDM)) then
1762            em1=2*APICH+API0+ADDM
1763 c     rho
1764         elseif(lb1.ge.25.and.lb1.le.27.and.em1.lt.(2*APICH+ADDM)) then
1765            em1=2*APICH+ADDM
1766 c     K*
1767         elseif(iabs(lb1).eq.30.and.em1.lt.(APICH+AK0+ADDM)) then
1768            em1=APICH+AK0+ADDM
1769 c     Delta created in ART may be below (n+pich) mass, causing LUDECY error:
1770         elseif(iabs(lb1).ge.6.and.iabs(lb1).le.9
1771      1          .and.em1.lt.(APICH+AN+ADDM)) then
1772            em1=APICH+AN+ADDM
1773         endif
1774 c        if(em1.ge.(em1a+0.01)) write (6,*)
1775 c     1       'Mass increase in resdec():',nt,em1-em1a,lb1
1776         e1=SQRT(EM1**2+PX1**2+PY1**2+PZ1**2)
1777         P(IP,4)=e1
1778         P(IP,5)=em1
1779 clin-5/2008:
1780         dpdecp=dpertp(i1)
1781 clin-4/2012 for option of pi0 decay:
1782         elseif(nt.eq.ntmax.and.ipi0dcy.eq.1.and.ipion.ge.1) then
1783            P(IP,1)=PPION(1,ipion,IRUN)
1784            P(IP,2)=PPION(2,ipion,IRUN)
1785            P(IP,3)=PPION(3,ipion,IRUN)
1786            P(IP,5)=EPION(ipion,IRUN)
1787            P(IP,4)=SQRT(P(IP,5)**2+P(IP,1)**2+P(IP,2)**2+P(IP,3)**2)
1788            dpdecp=dppion(ipion,IRUN)
1789 ctest off
1790 c           write(99,*) P(IP,4), P(IP,5), dpdecp, ipion, wid
1791         else
1792            print *, 'stopped in resdec() a'
1793            stop
1794         endif
1795 c
1796         call ludecy(IP)
1797 c     add decay time to daughter's formation time at the last timestep:
1798         if(nt.eq.ntmax) then
1799            tau0=hbarc/wid
1800            taudcy=tau0*(-1.)*alog(1.-RANART(NSEED))
1801            ndaut=n-nsav
1802            if(ndaut.le.1) then
1803               write(10,*) 'note: ndaut(<1)=',ndaut
1804               call lulist(2)
1805               stop
1806             endif
1807 c     lorentz boost:
1808 clin-4/2012 for option of pi0 decay:
1809             if(ipion.eq.0) then
1810                taudcy=taudcy*e1/em1
1811                tfnl=tfnl+taudcy
1812                xfnl=xfnl+px1/e1*taudcy
1813                yfnl=yfnl+py1/e1*taudcy
1814                zfnl=zfnl+pz1/e1*taudcy
1815             elseif(ipion.ge.1) then
1816                taudcy=taudcy*P(IP,4)/P(IP,5)
1817                tfnl=tfdpi(ipion,IRUN)+taudcy
1818                xfnl=RPION(1,ipion,IRUN)+P(IP,1)/P(IP,4)*taudcy
1819                yfnl=RPION(2,ipion,IRUN)+P(IP,2)/P(IP,4)*taudcy
1820                zfnl=RPION(3,ipion,IRUN)+P(IP,3)/P(IP,4)*taudcy
1821             else
1822                print *, 'stopped in resdec() b',ipion,wid,P(ip,4)
1823                stop
1824             endif
1825 c     at the last timestep, assign rho, K0S or eta (decay daughter)
1826 c     to lb(i1) only (not to lpion) in order to decay them again:
1827 clin-4/2012 for option of pi0 decay:
1828 c           if(n.ge.(nsav+2)) then
1829            if(n.ge.(nsav+2).and.ipion.eq.0) then
1830               do 1001 idau=nsav+2,n
1831                  kdaut=K(idau,2)
1832                  if(kdaut.eq.221.or.kdaut.eq.113
1833      1                .or.kdaut.eq.213.or.kdaut.eq.-213
1834      2                .or.kdaut.eq.310) then
1835 c     switch idau and i1(nsav+1):
1836                     ksave=kdaut
1837                     pxsave=p(idau,1)
1838                     pysave=p(idau,2)
1839                     pzsave=p(idau,3)
1840                     esave=p(idau,4)
1841                     xmsave=p(idau,5)
1842                     K(idau,2)=K(nsav+1,2)
1843                     p(idau,1)=p(nsav+1,1)
1844                     p(idau,2)=p(nsav+1,2)
1845                     p(idau,3)=p(nsav+1,3)
1846                     p(idau,4)=p(nsav+1,4)
1847                     p(idau,5)=p(nsav+1,5)
1848                     K(nsav+1,2)=ksave
1849                     p(nsav+1,1)=pxsave
1850                     p(nsav+1,2)=pysave
1851                     p(nsav+1,3)=pzsave
1852                     p(nsav+1,4)=esave
1853                     p(nsav+1,5)=xmsave
1854 c     note: phi decay may produce rho, K0s or eta, N*(1535) decay may produce
1855 c     eta, but only one daughter may be rho, K0s or eta:
1856                     goto 111
1857                  endif
1858  1001         continue
1859            endif
1860  111       continue
1861 c
1862            enet=0.
1863            do 1002 idau=nsav+1,n
1864               enet=enet+p(idau,4)
1865  1002      continue
1866 c           if(abs(enet-e1).gt.0.02)
1867 c     1          write(93,*) 'resdec(): nt=',nt,enet-e1,lb1
1868         endif
1869
1870         do 1003 idau=nsav+1,n
1871            kdaut=K(idau,2)
1872            lbdaut=IARFLV(kdaut)
1873 c     K0S and K0L are named K+/K- during hadron cascade, and only
1874 c     at the last timestep they keep their real LB # before output;
1875 c     K0/K0bar (from K* decay) converted to K0S and K0L at the last timestep:
1876            if(nt.eq.ntmax.and.(kdaut.eq.130.or.kdaut.eq.310
1877      1          .or.iabs(kdaut).eq.311)) then
1878               if(kdaut.eq.130) then
1879                  lbdaut=22
1880               elseif(kdaut.eq.310) then
1881                  lbdaut=24
1882               elseif(iabs(kdaut).eq.311) then
1883                  if(RANART(NSEED).lt.0.5) then
1884                     lbdaut=22
1885                  else
1886                     lbdaut=24
1887                  endif
1888               endif
1889            endif
1890 c
1891            if(idau.eq.(nsav+1)) then
1892 clin-4/2012 for option of pi0 decay:
1893               if(ipion.eq.0) then
1894                  LB(i1)=lbdaut
1895                  E(i1)=p(idau,5)
1896                  px1n=p(idau,1)
1897                  py1n=p(idau,2)
1898                  pz1n=p(idau,3)
1899 clin-5/2008:
1900                  dp1n=dpdecp
1901               elseif(ipion.ge.1) then
1902                  LPION(ipion,IRUN)=lbdaut
1903                  EPION(ipion,IRUN)=p(idau,5)
1904                  PPION(1,ipion,IRUN)=p(idau,1)
1905                  PPION(2,ipion,IRUN)=p(idau,2)
1906                  PPION(3,ipion,IRUN)=p(idau,3)
1907                  RPION(1,ipion,IRUN)=xfnl
1908                  RPION(2,ipion,IRUN)=yfnl
1909                  RPION(3,ipion,IRUN)=zfnl
1910                  tfdpi(ipion,IRUN)=tfnl
1911                  dppion(ipion,IRUN)=dpdecp
1912               endif
1913 c
1914            else
1915               nnn=nnn+1
1916               LPION(NNN,IRUN)=lbdaut
1917               EPION(NNN,IRUN)=p(idau,5)
1918               PPION(1,NNN,IRUN)=p(idau,1)
1919               PPION(2,NNN,IRUN)=p(idau,2)
1920               PPION(3,NNN,IRUN)=p(idau,3)
1921               RPION(1,NNN,IRUN)=xfnl
1922               RPION(2,NNN,IRUN)=yfnl
1923               RPION(3,NNN,IRUN)=zfnl
1924               tfdpi(NNN,IRUN)=tfnl
1925 clin-5/2008:
1926               dppion(NNN,IRUN)=dpdecp
1927            endif
1928  1003   continue
1929         return
1930         end
1931
1932 c=======================================================================
1933         subroutine inidcy
1934
1935         COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
1936 cc      SAVE /LUJETS/
1937         common/resdcy/NSAV,iksdcy
1938 cc      SAVE /resdcy/
1939         SAVE
1940         N=1
1941         NSAV=N
1942         return
1943         end
1944
1945 c=======================================================================
1946 clin-6/06/02 local parton freezeout motivated from critical density:
1947         subroutine local(t)
1948 c
1949         implicit double precision  (a-h, o-z)
1950         PARAMETER (MAXPTN=400001)
1951         PARAMETER (r0=1d0)
1952         COMMON /para1/ mul
1953 cc      SAVE /para1/
1954         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
1955      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
1956      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
1957 cc      SAVE /prec2/
1958         common /frzprc/
1959      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1960      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1961      &       xmfrz(MAXPTN),
1962      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1963 cc      SAVE /frzprc/
1964         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1965 cc      SAVE /prec4/
1966         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1967 cc      SAVE /prec5/
1968         common /coal/dpcoal,drcoal,ecritl
1969 cc      SAVE /coal/
1970         SAVE
1971 c
1972       do 1001 it=1,301
1973          if(t.ge.tfrz(it).and.t.lt.tfrz(it+1)) then
1974             if(it.eq.itlast) then
1975                return
1976             else
1977                itlast=it
1978                goto 50
1979             endif
1980          endif
1981  1001 continue
1982       write(1,*) 'local time out of range in LOCAL, stop',t,it
1983       stop
1984  50   continue
1985 c
1986       do 200 ip=1,mul
1987 c     skip partons which have frozen out:
1988          if(ifrz(ip).eq.1) goto 200
1989          if(it.eq.301) then
1990 c     freezeout all the left partons beyond the time of 3000 fm:
1991             etcrit=1d6
1992             goto 150
1993          else
1994 c     freezeout when transverse energy density < etcrit:
1995             etcrit=(ecritl*2d0/3d0)
1996          endif
1997 c     skip partons which have not yet formed:
1998          if(t.lt.FT5(ip)) goto 200
1999          rap0=rap(ip)
2000          eta0=eta(ip)
2001          x0=GX5(ip)+vx(ip)*(t-FT5(ip))
2002          y0=GY5(ip)+vy(ip)*(t-FT5(ip))
2003          detdy=0d0
2004          do 100 itest=1,mul
2005 c     skip self and partons which have not yet formed:
2006             if(itest.eq.ip.or.t.lt.FT5(itest)) goto 100
2007             ettest=eta(itest)
2008             xtest=GX5(itest)+vx(itest)*(t-FT5(itest))
2009             ytest=GY5(itest)+vy(itest)*(t-FT5(itest))
2010             drt=sqrt((xtest-x0)**2+(ytest-y0)**2)
2011 c     count partons within drt<1 and -1<(eta-eta0)<1:
2012             if(dabs(ettest-eta0).le.1d0.and.drt.le.r0)
2013      1           detdy=detdy+dsqrt(PX5(itest)**2+PY5(itest)**2
2014      2           +XMASS5(itest)**2)*0.5d0
2015  100     continue
2016          detdy=detdy*(dcosh(eta0)**2)/(t*3.1416d0*r0**2*dcosh(rap0))
2017 c     when density is below critical density for phase transition, freeze out:
2018  150     if(detdy.le.etcrit) then
2019             ifrz(ip)=1
2020             idfrz(ip)=ITYP5(ip)
2021             pxfrz(ip)=PX5(ip)
2022             pyfrz(ip)=PY5(ip)
2023             pzfrz(ip)=PZ5(ip)
2024             efrz(ip)=E5(ip)
2025             xmfrz(ip)=XMASS5(ip)
2026             if(t.gt.FT5(ip)) then
2027                gxfrz(ip)=x0
2028                gyfrz(ip)=y0
2029                gzfrz(ip)=GZ5(ip)+vz(ip)*(t-FT5(ip))
2030                ftfrz(ip)=t
2031             else
2032 c     if this freezeout time < formation time, use formation time & positions.
2033 c     This ensures the recovery of default hadron when e_crit=infty:
2034                gxfrz(ip)=GX5(ip)
2035                gyfrz(ip)=GY5(ip)
2036                gzfrz(ip)=GZ5(ip)
2037                ftfrz(ip)=FT5(ip)
2038             endif
2039          endif
2040  200  continue
2041 c
2042         return
2043         end
2044
2045 c=======================================================================
2046 clin-6/06/02 initialization for local parton freezeout
2047         subroutine inifrz
2048 c
2049         implicit double precision  (a-h, o-z)
2050         PARAMETER (MAXPTN=400001)
2051         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2052 cc      SAVE /ilist5/
2053         common /frzprc/
2054      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
2055      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
2056      &       xmfrz(MAXPTN),
2057      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
2058 cc      SAVE /frzprc/
2059         SAVE
2060 c
2061 c     for freezeout time 0-10fm, use interval of 0.1fm;
2062 c     for 10-100fm, use interval of 1fm;
2063 c     for 100-1000fm, use interval of 10fm;
2064 c     for 1000-3000fm, use interval of 100fm:
2065         step1=0.1d0
2066         step2=1d0
2067         step3=10d0
2068         step4=100d0
2069 c
2070         do 1001 it=1,101
2071            tfrz(it)=0d0+dble(it-1)*step1
2072  1001 continue
2073         do 1002 it=102,191
2074            tfrz(it)=10d0+dble(it-101)*step2
2075  1002   continue
2076         do 1003 it=192,281
2077            tfrz(it)=100d0+dble(it-191)*step3
2078  1003   continue
2079         do 1004 it=282,301
2080            tfrz(it)=1000d0+dble(it-281)*step4
2081  1004   continue
2082         tfrz(302)=tlarge
2083 c
2084         return
2085         end
2086
2087 clin-5/2009 v2 analysis
2088 c=======================================================================
2089 c     idd=0,1,2,3 specifies different subroutines for partonic flow analysis.
2090       subroutine flowp(idd)
2091 c
2092         implicit double precision  (a-h, o-z)
2093         real dt
2094         parameter (MAXPTN=400001)
2095 csp
2096         parameter (bmt=0.05d0)
2097         dimension nlfile(3),nsfile(3),nmfile(3)
2098 c
2099         dimension v2pp(3),xnpp(3),v2psum(3),v2p2sm(3),nfile(3)
2100         dimension tsp(31),v2pevt(3),v2pavg(3),varv2p(3)
2101         common /ilist1/
2102      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2103      &     ictype, icsta(MAXPTN),
2104      &     nic(MAXPTN), icels(MAXPTN)
2105 cc      SAVE /ilist1/
2106         COMMON /para1/ mul
2107 cc      SAVE /para1/
2108         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
2109      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
2110      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
2111 cc      SAVE /prec2/
2112         COMMON /pflow/ v2p(30,3),xnpart(30,3),etp(30,3),
2113      1       s2p(30,3),v2p2(30,3),nevt(30)
2114 cc      SAVE /pflow/
2115         COMMON /pflowf/ v2pf(30,3),xnpf(30,3),etpf(30,3),
2116      1                 xncoll(30),s2pf(30,3),v2pf2(30,3)
2117 cc      SAVE /pflowf/
2118         COMMON /pfrz/ v2pfrz(30,3),xnpfrz(30,3),etpfrz(30,3),
2119      1       s2pfrz(30,3),v2p2fz(30,3),tscatt(31),
2120      2       nevtfz(30),iscatt(30)
2121 cc      SAVE /pfrz/
2122         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2123      1 v2h2(30,3),s2h(30,3)
2124 cc      SAVE /hflow/
2125         COMMON /AREVT/ IAEVT, IARUN, MISS
2126 cc      SAVE /AREVT/
2127         common/anim/nevent,isoft,isflag,izpc
2128 cc      SAVE /anim/
2129         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2130 cc      SAVE /input1/
2131         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2132      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2133 cc      SAVE /INPUT2/
2134 cc      SAVE itimep,iaevtp,v2pp,xnpp,v2psum,v2p2sm
2135 cc      SAVE nfile,itanim,nlfile,nsfile,nmfile
2136         SAVE
2137 csp
2138         dimension etpl(30,3),etps(30,3),etplf(30,3),etpsf(30,3),
2139      &       etlfrz(30,3),etsfrz(30,3),
2140      &       xnpl(30,3),xnps(30,3),xnplf(30,3),xnpsf(30,3),
2141      &       xnlfrz(30,3),xnsfrz(30,3),
2142      &       v2pl(30,3),v2ps(30,3),v2plf(30,3),v2psf(30,3),
2143      &       s2pl(30,3),s2ps(30,3),s2plf(30,3),s2psf(30,3),
2144      &       DMYil(50,3),DMYfl(50,3),
2145      &       DMYis(50,3),DMYfs(50,3)
2146         data tsp/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
2147      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
2148      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
2149 c     idd=0: initialization for flow analysis, called by artdri.f:
2150         if(idd.eq.0) then
2151            nfile(1)=60
2152            nfile(2)=64
2153            nfile(3)=20
2154            OPEN (nfile(1),FILE='ana1/v2p.dat', STATUS = 'UNKNOWN')
2155            OPEN (nfile(1)+1,
2156      1 FILE = 'ana1/v2p-formed.dat', STATUS = 'UNKNOWN')
2157            OPEN (nfile(1)+2,
2158      1 FILE = 'ana1/v2p-active.dat', STATUS = 'UNKNOWN')
2159            OPEN (nfile(1)+3,
2160      1 FILE = 'ana1/v2ph.dat', STATUS = 'UNKNOWN')
2161            OPEN (nfile(2),FILE='ana1/v2p-y2.dat', STATUS = 'UNKNOWN')
2162            OPEN (nfile(2)+1,
2163      1 FILE = 'ana1/v2p-formed2.dat', STATUS = 'UNKNOWN')
2164            OPEN (nfile(2)+2,
2165      1 FILE = 'ana1/v2p-active2.dat', STATUS = 'UNKNOWN')
2166            OPEN (nfile(2)+3,
2167      1 FILE = 'ana1/v2ph-y2.dat', STATUS = 'UNKNOWN')
2168            OPEN (nfile(3),FILE='ana1/v2p-y1.dat', STATUS = 'UNKNOWN')
2169            OPEN (nfile(3)+1,
2170      1 FILE = 'ana1/v2p-formed1.dat', STATUS = 'UNKNOWN')
2171            OPEN (nfile(3)+2,
2172      1 FILE = 'ana1/v2p-active1.dat', STATUS = 'UNKNOWN')
2173            OPEN (nfile(3)+3,
2174      1 FILE = 'ana1/v2ph-y1.dat', STATUS = 'UNKNOWN')
2175            OPEN (49, FILE = 'ana1/v2p-ebe.dat', STATUS = 'UNKNOWN')
2176            write(49, *) '    ievt,  v2p,  v2p_y2,   v2p_y1'
2177 c
2178            OPEN (59, FILE = 'ana1/v2h.dat', STATUS = 'UNKNOWN')
2179            OPEN (68, FILE = 'ana1/v2h-y2.dat', STATUS = 'UNKNOWN')
2180            OPEN (69, FILE = 'ana1/v2h-y1.dat', STATUS = 'UNKNOWN')
2181            OPEN (88, FILE = 'ana1/v2h-ebe.dat', STATUS = 'UNKNOWN')
2182            write(88, *) '    ievt,  v2h,  v2h_y2,   v2h_y1'
2183 csp07/05
2184            nlfile(1)=70
2185            nlfile(2)=72
2186            nlfile(3)=74
2187            OPEN (nlfile(1),FILE='ana1/mtl.dat', STATUS = 'UNKNOWN')
2188            OPEN (nlfile(1)+1,
2189      1 FILE = 'ana1/mtl-formed.dat', STATUS = 'UNKNOWN')
2190            OPEN (nlfile(2),FILE='ana1/mtl-y2.dat', STATUS = 'UNKNOWN')
2191            OPEN (nlfile(2)+1,
2192      1 FILE = 'ana1/mtl-formed2.dat', STATUS = 'UNKNOWN')
2193            OPEN (nlfile(3),FILE='ana1/mtl-y1.dat', STATUS = 'UNKNOWN')
2194            OPEN (nlfile(3)+1,
2195      1 FILE = 'ana1/mtl-formed1.dat', STATUS = 'UNKNOWN')
2196            nsfile(1)=76
2197            nsfile(2)=78
2198            nsfile(3)=80
2199            OPEN (nsfile(1),FILE='ana1/mts.dat', STATUS = 'UNKNOWN')
2200            OPEN (nsfile(1)+1,
2201      1 FILE = 'ana1/mts-formed.dat', STATUS = 'UNKNOWN')
2202            OPEN (nsfile(2),FILE='ana1/mts-y2.dat', STATUS = 'UNKNOWN')
2203            OPEN (nsfile(2)+1,
2204      1 FILE = 'ana1/mts-formed2.dat', STATUS = 'UNKNOWN')
2205            OPEN (nsfile(3),FILE='ana1/mts-y1.dat', STATUS = 'UNKNOWN')
2206            OPEN (nsfile(3)+1,
2207      1 FILE = 'ana1/mts-formed1.dat', STATUS = 'UNKNOWN')
2208            nmfile(1)=82
2209            nmfile(2)=83
2210            nmfile(3)=84
2211            OPEN (nmfile(1),FILE='ana1/Nmt.dat', STATUS = 'UNKNOWN')
2212            OPEN (nmfile(2),FILE='ana1/Nmt-y2.dat', STATUS = 'UNKNOWN')
2213            OPEN (nmfile(3),FILE='ana1/Nmt-y1.dat', STATUS = 'UNKNOWN')
2214 clin-11/27/00 for animation:
2215            if(nevent.eq.1) then
2216 c           OPEN (91, FILE = 'ana1/h-animate.dat', STATUS = 'UNKNOWN')
2217 c           write(91,*) ntmax, dt
2218 c           OPEN (92, FILE = 'ana1/p-animate.dat', STATUS = 'UNKNOWN')
2219            OPEN (93, FILE = 'ana1/p-finalft.dat', STATUS = 'UNKNOWN')
2220            endif
2221 c
2222            itimep=0
2223            itanim=0
2224            iaevtp=0
2225 csp
2226            do 1002 ii=1,50
2227               do 1001 iy=1,3
2228                  DMYil(ii,iy) = 0d0
2229                  DMYfl(ii,iy) = 0d0
2230                  DMYis(ii,iy) = 0d0
2231                  DMYfs(ii,iy) = 0d0
2232  1001         continue
2233  1002      continue
2234 c
2235            do 1003 ii=1,31
2236               tscatt(ii)=0d0
2237  1003      continue
2238            do 1005 ii=1,30
2239               nevt(ii)=0
2240               xncoll(ii)=0d0
2241               nevtfz(ii)=0
2242               iscatt(ii)=0
2243               do 1004 iy=1,3
2244                  v2p(ii,iy)=0d0
2245                  v2p2(ii,iy)=0d0
2246                  s2p(ii,iy)=0d0
2247                  etp(ii,iy)=0d0
2248                  xnpart(ii,iy)=0d0
2249                  v2pf(ii,iy)=0d0
2250                  v2pf2(ii,iy)=0d0
2251                  s2pf(ii,iy)=0d0
2252                  etpf(ii,iy)=0d0
2253                  xnpf(ii,iy)=0d0
2254                  v2pfrz(ii,iy)=0d0
2255                  v2p2fz(ii,iy)=0d0
2256                  s2pfrz(ii,iy)=0d0
2257                  etpfrz(ii,iy)=0d0
2258                  xnpfrz(ii,iy)=0d0
2259 csp07/05
2260                  etpl(ii,iy)=0d0
2261                  etps(ii,iy)=0d0
2262                  etplf(ii,iy)=0d0
2263                  etpsf(ii,iy)=0d0
2264                  etlfrz(ii,iy)=0d0
2265                  etsfrz(ii,iy)=0d0
2266               xnpl(ii,iy)=0d0
2267               xnps(ii,iy)=0d0
2268               xnplf(ii,iy)=0d0
2269               xnpsf(ii,iy)=0d0
2270               xnlfrz(ii,iy)=0d0
2271               xnsfrz(ii,iy)=0d0
2272               v2pl(ii,iy)=0d0
2273               v2ps(ii,iy)=0d0
2274               v2plf(ii,iy)=0d0
2275               v2psf(ii,iy)=0d0
2276               s2pl(ii,iy)=0d0
2277               s2ps(ii,iy)=0d0
2278               s2plf(ii,iy)=0d0
2279               s2psf(ii,iy)=0d0
2280  1004      continue
2281  1005   continue
2282            do 1006 iy=1,3
2283               v2pevt(iy)=0d0
2284               v2pavg(iy)=0d0
2285               varv2p(iy)=0d0
2286               v2pp(iy)=0.d0
2287               xnpp(iy)=0d0
2288               v2psum(iy)=0.d0
2289               v2p2sm(iy)=0.d0
2290  1006      continue
2291 c     idd=1: calculate parton elliptic flow, called by zpc.f:
2292         else if(idd.eq.1) then
2293            t2time = FT5(iscat)
2294            do 1008 ianp = 1, 30
2295               if (t2time.lt.tsp(ianp+1).and.t2time.ge.tsp(ianp)) then
2296 c     write flow info only once at each fixed time:
2297                  xncoll(ianp)=xncoll(ianp)+1d0
2298 c     to prevent an earlier t2time comes later in the same event
2299 c     and mess up nevt:
2300                  if(ianp.le.itimep.and.iaevt.eq.iaevtp) goto 101
2301                  nevt(ianp)=nevt(ianp)+1
2302                  tscatt(ianp+1)=t2time
2303                  iscatt(ianp)=1
2304                  nevtfz(ianp)=nevtfz(ianp)+1
2305                  do 100 I=1,mul
2306             ypartn=0.5d0*dlog((E5(i)+PZ5(i))/(E5(i)-PZ5(i)+1.d-8))
2307                     pt2=PX5(I)**2+PY5(I)**2
2308 ctest off: initial (pt,y) and (x,y) distribution:
2309 c                    idtime=1
2310 c                    if(ianp.eq.idtime) then
2311 c                       iityp=iabs(ITYP5(I))
2312 c                       if(iityp.eq.1.or.iityp.eq.2) then
2313 c                          write(651,*) dsqrt(pt2),ypartn
2314 c                          write(654,*) GX5(I),GY5(I)
2315 c                       elseif(iityp.eq.1103.or.iityp.eq.2101
2316 c     1 .or.iityp.eq.2103.or.iityp.eq.2203.
2317 c     2 .or.iityp.eq.3101.or.iityp.eq.3103.
2318 c     3 .or.iityp.eq.3201.or.iityp.eq.3203.or.iityp.eq.3303)
2319 c     4 then
2320 c                          write(652,*) dsqrt(pt2),ypartn
2321 c                          write(655,*) GX5(I),GY5(I)
2322 c                       elseif(iityp.eq.21) then
2323 c                          write(653,*) dsqrt(pt2),ypartn
2324 c                          write(656,*) GX5(I),GY5(I)
2325 c                       endif
2326 c                    endif
2327 ctest-end
2328 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2329 c                    gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2330 c                    gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2331 c                    gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2332 c                    rt_now=dsqrt(gx_now**2+gy_now**2)
2333 c                    zmax=0.1d0*t2time
2334 c                    volume=3.1416d0*(2d0**2)*(2*zmax)
2335 c                    if(rt_now.gt.2d0.or.dabs(gz_now).gt.zmax)
2336 c     1                   goto 100
2337 ctest-end
2338                     iloop=1
2339                     if(dabs(ypartn).le.1d0) then
2340                        iloop=2
2341                        if(dabs(ypartn).le.0.5d0) then
2342                           iloop=3
2343                        endif
2344                     endif
2345                     do 50 iy=1,iloop
2346 clin-5/2012:
2347 c                       if(pt2.gt.0.) then
2348                        if(pt2.gt.0d0) then
2349                           v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2350 clin-5/2012:
2351 c                          if(abs(v2prtn).gt.1.)
2352                           if(dabs(v2prtn).gt.1d0)
2353      1 write(nfile(iy),*) 'v2prtn>1',v2prtn
2354                           v2p(ianp,iy)=v2p(ianp,iy)+v2prtn
2355                           v2p2(ianp,iy)=v2p2(ianp,iy)+v2prtn**2
2356                        endif
2357                        xperp2=GX5(I)**2+GY5(I)**2
2358 clin-5/2012:
2359 c                       if(xperp2.gt.0.)
2360                        if(xperp2.gt.0d0)
2361      1        s2p(ianp,iy)=s2p(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2362                        xnpart(ianp,iy)=xnpart(ianp,iy)+1d0
2363                        etp(ianp,iy)=etp(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2364 ctest off density:
2365 c                       etp(ianp,iy)=etp(ianp,iy)
2366 c     1                  +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2367 clin-2/22/00 to write out parton info only for formed ones:
2368                        if(FT5(I).le.t2time) then
2369                           v2pf(ianp,iy)=v2pf(ianp,iy)+v2prtn
2370                           v2pf2(ianp,iy)=v2pf2(ianp,iy)+v2prtn**2
2371 clin-5/2012:
2372 c                          if(xperp2.gt.0.)
2373                           if(xperp2.gt.0d0)
2374      1        s2pf(ianp,iy)=s2pf(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2375                           xnpf(ianp,iy)=xnpf(ianp,iy)+1d0
2376                   etpf(ianp,iy)=etpf(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2377 ctest off density:
2378 c                  etpf(ianp,iy)=etpf(ianp,iy)
2379 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2380                        endif
2381  50                    continue
2382  100                 continue
2383                  itimep=ianp
2384                  iaevtp=iaevt
2385 clin-3/30/00 ebe v2 variables:
2386                  if(ianp.eq.30) then
2387                     do 1007 iy=1,3
2388                        npartn=IDINT(xnpart(ianp,iy)-xnpp(iy))
2389                        if(npartn.ne.0) then
2390                           v2pevt(iy)=(v2p(ianp,iy)-v2pp(iy))/npartn
2391                           v2psum(iy)=v2psum(iy)+v2pevt(iy)
2392                           v2p2sm(iy)=v2p2sm(iy)+v2pevt(iy)**2
2393                           v2pp(iy)=v2p(ianp,iy)
2394                           xnpp(iy)=xnpart(ianp,iy)
2395                        endif
2396  1007               continue
2397                     write(49, 160) iaevt,v2pevt
2398                  endif
2399                  goto 101
2400               endif
2401  1008      continue
2402 clin-11/28/00 for animation:
2403  101           if(nevent.eq.1) then
2404               do 110 nt = 1, ntmax
2405                  time1=dble(nt*dt)
2406                  time2=dble((nt+1)*dt)
2407                  if (t2time.lt.time2.and.t2time.ge.time1) then
2408                     if(nt.le.itanim) return
2409 c                    write(92,*) t2time
2410                     iform=0
2411                     do 1009 I=1,mul
2412 c     write out parton info only for formed ones:
2413                        if(FT5(I).le.t2time) then
2414                           iform=iform+1
2415                        endif
2416  1009               continue
2417 c                    write(92,*) iform
2418                     do 120 I=1,mul
2419                        if(FT5(I).le.t2time) then
2420 clin-11/29/00-ctest off calculate parton coordinates after propagation:
2421 c                          gx_now=GX5(i)+(t2time-FT5(i))*PX5(i)/E5(i)
2422 c                          gy_now=GY5(i)+(t2time-FT5(i))*PY5(i)/E5(i)
2423 c                          gz_now=GZ5(i)+(t2time-FT5(i))*PZ5(i)/E5(i)
2424 c          write(92,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2425 c          write(92,180) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i),
2426 c     1    PX5(i),PY5(i),PZ5(i),E5(i)
2427 ctest-end
2428                        endif
2429  120                    continue
2430                     itanim=nt
2431                  endif
2432  110              continue
2433            endif
2434 c
2435  140           format(i10,4(2x,f7.2))
2436  160           format(i10,3(2x,f9.5))
2437 c 180           format(i6,8(1x,f7.2))
2438 clin-5/17/01 calculate v2 for active partons (which have not frozen out):
2439 c     idd=3, called at end of zpc.f:
2440         else if(idd.eq.3) then
2441            do 1010 ianp=1,30
2442               if(iscatt(ianp).eq.0) tscatt(ianp+1)=tscatt(ianp)
2443  1010      continue
2444            do 350 I=1,mul
2445               ypartn=0.5d0*dlog((E5(i)+PZ5(i)+1.d-8)
2446      1 /(E5(i)-PZ5(i)+1.d-8))
2447               pt2=PX5(I)**2+PY5(I)**2
2448               iloop=1
2449               if(dabs(ypartn).le.1d0) then
2450                  iloop=2
2451                  if(dabs(ypartn).le.0.5d0) then
2452                     iloop=3
2453                  endif
2454               endif
2455 c
2456               do 325 ianp=1,30
2457                  if(iscatt(ianp).ne.0) then
2458                     if(FT5(I).lt.tscatt(ianp+1)
2459      1 .and.FT5(I).ge.tscatt(ianp)) then
2460                        do 1011 iy=1,iloop
2461 clin-5/2012:
2462 c                          if(pt2.gt.0.) then
2463                           if(pt2.gt.0d0) then
2464                              v2prtn=(PX5(I)**2-PY5(I)**2)/pt2
2465                              v2pfrz(ianp,iy)=v2pfrz(ianp,iy)+v2prtn
2466                      v2p2fz(ianp,iy)=v2p2fz(ianp,iy)+v2prtn**2
2467                           endif
2468                           xperp2=GX5(I)**2+GY5(I)**2
2469 clin-5/2012:
2470 c                          if(xperp2.gt.0.) s2pfrz(ianp,iy)=
2471                           if(xperp2.gt.0d0) s2pfrz(ianp,iy)=
2472      1 s2pfrz(ianp,iy)+(GX5(I)**2-GY5(I)**2)/xperp2
2473         etpfrz(ianp,iy)=etpfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2474                           xnpfrz(ianp,iy)=xnpfrz(ianp,iy)+1d0
2475 ctest off density:
2476 c                    etpfrz(ianp,iy)=etpfrz(ianp,iy)
2477 c     1                   +dsqrt(pt2+XMASS5(I)**2+PZ5(i)**2)/volume
2478 csp07/05
2479             if(ITYP5(I).eq.1.or.ITYP5(I).eq.2)then
2480               etlfrz(ianp,iy)=etlfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2481               xnlfrz(ianp,iy)=xnlfrz(ianp,iy)+1d0
2482             elseif(ITYP5(I).eq.3)then
2483               etsfrz(ianp,iy)=etsfrz(ianp,iy)+dsqrt(pt2+XMASS5(I)**2)
2484               xnsfrz(ianp,iy)=xnsfrz(ianp,iy)+1d0
2485             endif
2486 csp07/05 end
2487  1011    continue
2488 c     parton freezeout info taken, proceed to next parton:
2489                        goto 350
2490                     endif
2491                  endif
2492  325          continue
2493  350       continue
2494 c     idd=2: calculate average partonic elliptic flow, called from artdri.f,
2495         else if(idd.eq.2) then
2496            do 1012 i=1,3
2497               write(nfile(i),*) '   tsp,   v2p,     v2p2, '//
2498      1 '   s2p,  etp,   xmult,    nevt,  xnctot'
2499               write ((nfile(i)+1),*) '  tsp,   v2pf,   v2pf2, '//
2500      1 '   s2pf, etpf,  xnform,  xcrate'
2501               write ((nfile(i)+2),*) '  tsp,   v2pa,   v2pa2, '//
2502      1 '   s2pa, etpa,  xmult_ap,  xnform,   nevt'
2503               write ((nfile(i)+3),*) '  tsph,  v2ph,   v2ph2, '//
2504      1 '   s2ph, etph,  xmult_(ap/2+h),xmult_ap/2,nevt'
2505 csp
2506            write(nlfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2507            write(nsfile(i),*) '   tsp,    v2,     s2,    etp,    xmul'
2508            write(nlfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2509            write(nsfile(i)+1,*) '   tsp,    v2,     s2,    etp,    xmul'
2510 c
2511  1012   continue
2512 clin-3/30/00 ensemble average & variance of v2 (over particles in all events):
2513            do 150 ii=1, 30
2514               if(nevt(ii).eq.0) goto 150
2515               do 1014 iy=1,3
2516 clin-5/2012:
2517 c                 if(xnpart(ii,iy).gt.1) then
2518                  if(xnpart(ii,iy).gt.1d0) then
2519                     v2p(ii,iy)=v2p(ii,iy)/xnpart(ii,iy)
2520                     v2p2(ii,iy)=dsqrt((v2p2(ii,iy)/xnpart(ii,iy)
2521      1                    -v2p(ii,iy)**2)/(xnpart(ii,iy)-1))
2522                     s2p(ii,iy)=s2p(ii,iy)/xnpart(ii,iy)
2523 c xmult and etp are multiplicity and et for an averaged event:
2524                     xmult=dble(xnpart(ii,iy)/dble(nevt(ii)))
2525                     etp(ii,iy)=etp(ii,iy)/dble(nevt(ii))
2526 csp
2527                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2528                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2529 c
2530                     xnctot=0d0
2531                     do 1013 inum=1,ii
2532                        xnctot=xnctot+xncoll(inum)
2533  1013               continue
2534                     if(nevt(1).ne.0) xnctot=xnctot/nevt(1)
2535                     write (nfile(iy),200) tsp(ii),v2p(ii,iy),
2536      1      v2p2(ii,iy),s2p(ii,iy),etp(ii,iy),xmult,nevt(ii),xnctot
2537                  endif
2538                  if(nevt(ii).ne.0)
2539      1                xcrate=xncoll(ii)/(tsp(ii+1)-tsp(ii))/nevt(ii)
2540 c
2541 clin-5/2012:
2542 c                 if(xnpf(ii,iy).gt.1) then
2543                  if(xnpf(ii,iy).gt.1d0) then
2544                     v2pf(ii,iy)=v2pf(ii,iy)/xnpf(ii,iy)
2545                     v2pf2(ii,iy)=dsqrt((v2pf2(ii,iy)/xnpf(ii,iy)
2546      1                    -v2pf(ii,iy)**2)/(xnpf(ii,iy)-1))
2547                     s2pf(ii,iy)=s2pf(ii,iy)/xnpf(ii,iy)
2548                     xnform=dble(xnpf(ii,iy)/dble(nevt(ii)))
2549                     etpf(ii,iy)=etpf(ii,iy)/dble(nevt(ii))
2550 csp
2551                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2552                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2553 c
2554                     write (nfile(iy)+1, 210) tsp(ii),v2pf(ii,iy),
2555      1      v2pf2(ii,iy),s2pf(ii,iy),etpf(ii,iy),xnform,xcrate
2556                  endif
2557 csp
2558 clin-5/2012:
2559 c                 if(xnpl(ii,iy).gt.1) then
2560                  if(xnpl(ii,iy).gt.1d0) then
2561                     v2pl(ii,iy)=v2pl(ii,iy)/xnpl(ii,iy)
2562                     s2pl(ii,iy)=s2pl(ii,iy)/xnpl(ii,iy)
2563                     xmult=dble(xnpl(ii,iy)/dble(nevt(ii)))
2564                     etpl(ii,iy)=etpl(ii,iy)/dble(nevt(ii))
2565                     write (nlfile(iy),201) tsp(ii),v2pl(ii,iy),
2566      1        s2pl(ii,iy),etpl(ii,iy),xmult
2567                  endif
2568 clin-5/2012:
2569 c                 if(xnps(ii,iy).gt.1) then
2570                  if(xnps(ii,iy).gt.1d0) then
2571                     v2ps(ii,iy)=v2ps(ii,iy)/xnps(ii,iy)
2572                     s2ps(ii,iy)=s2ps(ii,iy)/xnps(ii,iy)
2573                     xmult=dble(xnps(ii,iy)/dble(nevt(ii)))
2574                     etps(ii,iy)=etps(ii,iy)/dble(nevt(ii))
2575                     write (nsfile(iy),201) tsp(ii),v2ps(ii,iy),
2576      1        s2ps(ii,iy),etps(ii,iy),xmult
2577                  endif
2578 clin-5/2012:
2579 c                 if(xnplf(ii,iy).gt.1) then
2580                  if(xnplf(ii,iy).gt.1d0) then
2581                     v2plf(ii,iy)=v2plf(ii,iy)/xnplf(ii,iy)
2582                     s2plf(ii,iy)=s2plf(ii,iy)/xnplf(ii,iy)
2583                     xmult=dble(xnplf(ii,iy)/dble(nevt(ii)))
2584                     etplf(ii,iy)=etplf(ii,iy)/dble(nevt(ii))
2585                     write (nlfile(iy)+1,201) tsp(ii),v2plf(ii,iy),
2586      1        s2plf(ii,iy),etplf(ii,iy),xmult
2587                  endif
2588 clin-5/2012:
2589 c                 if(xnpsf(ii,iy).gt.1) then
2590                  if(xnpsf(ii,iy).gt.1d0) then
2591                     v2psf(ii,iy)=v2psf(ii,iy)/xnpsf(ii,iy)
2592                     s2psf(ii,iy)=s2psf(ii,iy)/xnpsf(ii,iy)
2593                     xmult=dble(xnpsf(ii,iy)/dble(nevt(ii)))
2594                     etpsf(ii,iy)=etpsf(ii,iy)/dble(nevt(ii))
2595                     write (nsfile(iy)+1,201) tsp(ii),v2psf(ii,iy),
2596      1        s2psf(ii,iy),etpsf(ii,iy),xmult
2597                  endif
2598 csp-end
2599  1014         continue
2600  150           continue
2601 csp07/05 initial & final mt distrb
2602                scalei=0d0
2603                scalef=0d0
2604                if(nevt(1).ne.0) SCALEi = 1d0 / dble(nevt(1)) / BMT
2605                if(nevt(30).ne.0) SCALEf = 1d0 / dble(nevt(30)) / BMT
2606          do 1016 iy=2,3
2607            yra = 1d0
2608            if(iy .eq. 2)yra = 2d0
2609          do 1015 i=1,50
2610            WRITE(nmfile(iy),251) BMT*dble(I - 0.5),
2611      &     SCALEi*DMYil(I,iy)/yra, SCALEf*DMYfl(I,iy)/yra,
2612      &     SCALEi*DMYis(I,iy)/yra, SCALEf*DMYfs(I,iy)/yra
2613  1015   continue
2614  1016 continue
2615 csp07/05 end
2616 clin-3/30/00 event-by-event average & variance of v2:
2617            if(nevt(30).ge.1) then
2618               do 1017 iy=1,3
2619                  v2pavg(iy)=v2psum(iy)/nevt(30)
2620                  v2var0=v2p2sm(iy)/nevt(30)-v2pavg(iy)**2
2621 clin-5/2012:
2622 c                 if(v2var0.gt.0) varv2p(iy)=dsqrt(v2var0)
2623                  if(v2var0.gt.0d0) varv2p(iy)=dsqrt(v2var0)
2624  1017 continue
2625               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): avg=', v2pavg
2626               write(49, 240) 'EBE v2p,v2p(y2),v2p(y1): var=', varv2p
2627            endif
2628 clin-11/28/00 for animation:
2629             if(nevent.eq.1) then
2630               do 1018 I=1,mul
2631                  if(FT5(I).le.t2time) then
2632          write(93,140) ITYP5(i),GX5(i),GY5(i),GZ5(i),FT5(i)
2633                  endif
2634  1018         continue
2635 clin-11/29/00 signal the end of animation file:
2636 c              write(91,*) -10.
2637 c              write(91,*) 0
2638 c              write(92,*) -10.
2639 c              write(92,*) 0
2640 c              close (91)
2641 c              close (92)
2642               close (93)
2643            endif
2644 clin-5/18/01 calculate v2 for active partons:
2645            do 450 ianp=1,30
2646               do 400 iy=1,3
2647                  v2pact=0d0
2648                  v2p2ac=0d0
2649                  s2pact=0d0
2650                  etpact=0d0
2651                  xnacti=0d0
2652 clin-5/2012:
2653 c                 if(xnpf(ianp,iy).gt.1) then
2654                  if(xnpf(ianp,iy).gt.1d0) then
2655 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2656                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2657                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2658      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2659                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2660                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2661                     xnpact=xnpf(ianp,iy)
2662 c
2663                     do 1019 kanp=1,ianp
2664                        v2pact=v2pact-v2pfrz(kanp,iy)
2665                        v2p2ac=v2p2ac-v2p2fz(kanp,iy)
2666                        s2pact=s2pact-s2pfrz(kanp,iy)
2667                        etpact=etpact-etpfrz(kanp,iy)
2668                        xnpact=xnpact-xnpfrz(kanp,iy)
2669  1019               continue
2670 c     save the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2671                     v2ph=v2pact
2672                     v2ph2=v2p2ac
2673                     s2ph=s2pact
2674                     etph=etpact
2675                     xnp2=xnpact/2d0
2676 c
2677 clin-5/2012:
2678 c                    if(xnpact.gt.1.and.nevt(ianp).ne.0) then
2679                     if(xnpact.gt.1d0.and.nevt(ianp).ne.0) then
2680                        v2pact=v2pact/xnpact
2681                        v2p2ac=dsqrt((v2p2ac/xnpact
2682      1                    -v2pact**2)/(xnpact-1))
2683                        s2pact=s2pact/xnpact
2684                        xnacti=dble(xnpact/dble(nevt(ianp)))
2685                        etpact=etpact/dble(nevt(ianp))
2686                        write (nfile(iy)+2, 250) tsp(ianp),v2pact,
2687      1 v2p2ac,s2pact,etpact,xnacti,
2688      2 xnpf(ianp,iy)/dble(nevt(ianp)),nevt(ianp)
2689                     endif
2690                  endif
2691 c     To calculate combined v2 for active partons plus formed hadrons,
2692 c     add the sum of v2h, v2h2, s2h, eth, and xnh for formed hadrons:
2693 c     scale the hadron part in case nevt(ianp) != nevent:
2694                  shadr=dble(nevt(ianp))/dble(nevent)
2695                  ianh=ianp
2696                  v2ph=v2ph+v2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2697                  v2ph2=v2ph2+(v2h2(ianh,iy)**2*(xnhadr(ianh,iy)-1)
2698      1 +v2h(ianh,iy)**2)*xnhadr(ianh,iy)*shadr
2699                  s2ph=s2ph+s2h(ianh,iy)*xnhadr(ianh,iy)*shadr
2700                  etph=etph+eth(ianh,iy)*dble(nevent)*shadr
2701                  xnph=xnpact+xnhadr(ianh,iy)*shadr
2702                  xnp2h=xnp2+xnhadr(ianh,iy)*shadr
2703 clin-5/2012:
2704 c                 if(xnph.gt.1) then
2705                  if(xnph.gt.1d0) then
2706                     v2ph=v2ph/xnph
2707                     v2ph2=dsqrt((v2ph2/xnph-v2ph**2)/(xnph-1))
2708                     s2ph=s2ph/xnph
2709                     etph=etph/dble(nevt(ianp))
2710                     xnp2=xnp2/dble(nevt(ianp))
2711                     xnp2h=xnp2h/dble(nevent)
2712                     if(tsp(ianp).le.(ntmax*dt))
2713      1                    write (nfile(iy)+3, 250) tsp(ianp),v2ph,
2714      2 v2ph2,s2ph,etph,xnp2h,xnp2,nevt(ianp)
2715                  endif
2716 c
2717  400              continue
2718  450       continue
2719            do 550 ianp=1,30
2720               do 500 iy=1,3
2721                  v2pact=0d0
2722                  v2p2ac=0d0
2723                  s2pact=0d0
2724                  etpact=0d0
2725                  xnacti=0d0
2726 c     reconstruct the sum of v2p, v2p2, s2p, etp, and xnp for formed partons:
2727                     v2pact=v2pf(ianp,iy)*xnpf(ianp,iy)
2728                     v2p2ac=(v2pf2(ianp,iy)**2*(xnpf(ianp,iy)-1)
2729      1 +v2pf(ianp,iy)**2)*xnpf(ianp,iy)
2730                     s2pact=s2pf(ianp,iy)*xnpf(ianp,iy)
2731                     etpact=etpf(ianp,iy)*dble(nevt(ianp))
2732                     xnpact=xnpf(ianp,iy)
2733  500              continue
2734  550           continue
2735            close (620)
2736            close (630)
2737            do 1021 nf=1,3
2738               do 1020 ifile=0,3
2739                  close(nfile(nf)+ifile)
2740  1020        continue
2741  1021     continue
2742            do 1022 nf=1,3
2743               close(740+nf)
2744  1022      continue
2745         endif
2746  200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2),i6,2x,f9.2)
2747  210        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2))
2748  240        format(a30,3(2x,f9.5))
2749  250        format(2x,f5.2,3(2x,f7.4),3(2x,f9.2),i6)
2750 csp
2751  201        format(2x,f5.2,4(2x,f9.2))
2752  251        format(5e15.5)
2753 c
2754         return
2755         end
2756
2757 c=======================================================================
2758 c     Calculate flow from formed hadrons, called by art1e.f:
2759 c     Note: numbers in art not in double precision!
2760         subroutine flowh(ct)
2761         PARAMETER (MAXSTR=150001, MAXR=1)
2762         dimension tsh(31)
2763         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2764         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,v2hevt(3)
2765         DOUBLE PRECISION  pt2, v2hadr
2766         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2767      1 v2h2(30,3),s2h(30,3)
2768 cc      SAVE /hflow/
2769         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2770 cc      SAVE /ebe/
2771         common /lastt/itimeh,bimp
2772 cc      SAVE /lastt/
2773         COMMON /RUN/ NUM
2774 cc      SAVE /RUN/
2775         COMMON  /AA/      R(3,MAXSTR)
2776 cc      SAVE /AA/
2777         COMMON  /BB/      P(3,MAXSTR)
2778 cc      SAVE /BB/
2779         COMMON  /CC/      E(MAXSTR)
2780 cc      SAVE /CC/
2781         COMMON  /EE/      ID(MAXSTR),LB(MAXSTR)
2782 cc      SAVE /EE/
2783         COMMON  /RR/      MASSR(0:MAXR)
2784 cc      SAVE /RR/
2785         common/anim/nevent,isoft,isflag,izpc
2786 cc      SAVE /anim/
2787         COMMON /AREVT/ IAEVT, IARUN, MISS
2788 cc      SAVE /AREVT/
2789         SAVE
2790 c
2791         do 1001 ii = 1, 31
2792            tsh(ii)=float(ii-1)
2793  1001   continue
2794 c
2795         do 1004 ianh = 1, 30
2796            if ((ct+0.0001).lt.tsh(ianh+1)
2797      1 .and.(ct+0.0001).ge.tsh(ianh)) then
2798               if(ianh.eq.itimeh) goto 101
2799               IA = 0
2800               DO 1002 J = 1, NUM
2801                  mult=MASSR(J)
2802                  IA = IA + MASSR(J - 1)
2803                  DO 100 IC = 1, mult
2804                     I = IA + IC
2805 c     5/04/01 exclude leptons and photons:
2806                     if(iabs(LB(I)-10000).lt.100) goto 100
2807                     px=p(1,i)
2808                     py=p(2,i)
2809                     pt2=dble(px)**2+dble(py)**2
2810 c     2/18/00 Note: e(i) gives the mass in ART:
2811                     ene=sqrt(e(i)**2+sngl(pt2)+p(3,i)**2)
2812                     RAP=0.5*alog((ene+p(3,i))/(ene-p(3,i)))
2813 ctest off density with 2fm radius and z:(-0.1*t,0.1*t):
2814 c                rt_now=sqrt(r(1,i)**2+r(2,i)**2)
2815 c                gz_now=r(3,i)
2816 c                zmax=0.1*ct
2817 c                volume=3.1416*(2.**2)*2*zmax
2818 c                if(rt_now.gt.2.or.abs(gz_now).gt.zmax)
2819 c     1               goto 100
2820                     iloop=1
2821                     if(abs(rap).le.1) then
2822                        iloop=2
2823                        if(abs(rap).le.0.5) then
2824                           iloop=3
2825                        endif
2826                     endif
2827                     do 50 iy=1,iloop
2828                        if(pt2.gt.0d0) then
2829                           v2hadr=(dble(px)**2-dble(py)**2)/pt2
2830                           v2h(ianh,iy)=v2h(ianh,iy)+v2hadr
2831                           v2h2(ianh,iy)=v2h2(ianh,iy)+v2hadr**2
2832                           if(dabs(v2hadr).gt.1d0)
2833      1 write(1,*) 'v2hadr>1',v2hadr,px,py
2834                        endif
2835                        xperp2=r(1,I)**2+r(2,I)**2
2836                        if(xperp2.gt.0.)
2837      1 s2h(ianh,iy)=s2h(ianh,iy)+dble((r(1,I)**2-r(2,I)**2)/xperp2)
2838                eth(ianh,iy)=eth(ianh,iy)+dble(SQRT(e(i)**2+sngl(pt2)))
2839 ctest off density:
2840 c               eth(ianh,iy)=eth(ianh,iy)
2841 c     1                  +dble(SQRT(e(i)**2+sngl(pt2)+p(3,i)**2))/volume
2842                        xnhadr(ianh,iy)=xnhadr(ianh,iy)+1d0
2843  50                    continue
2844  100                 continue
2845  1002         CONTINUE
2846               itimeh=ianh
2847 clin-5/04/01 ebe v2 variables:
2848               if(ianh.eq.30) then
2849                  do 1003 iy=1,3
2850                     nhadrn=IDINT(xnhadr(ianh,iy)-xnhadp(iy))
2851                     if(nhadrn.ne.0) then
2852                v2hevt(iy)=(v2h(ianh,iy)-v2hp(iy))/dble(nhadrn)
2853                        v2hsum(iy)=v2hsum(iy)+v2hevt(iy)
2854                        v2h2sm(iy)=v2h2sm(iy)+v2hevt(iy)**2
2855                        v2hp(iy)=v2h(ianh,iy)
2856                        xnhadp(iy)=xnhadr(ianh,iy)
2857                     endif
2858  1003            continue
2859                  write(88, 160) iaevt,v2hevt
2860               endif
2861               goto 101
2862            endif
2863  1004   continue
2864  160        format(i10,3(2x,f9.5))
2865 clin-11/27/00 for animation:
2866  101        if(nevent.eq.1) then
2867            IA = 0
2868            do 1005 J = 1, NUM
2869               mult=MASSR(J)
2870               IA = IA + MASSR(J - 1)
2871 c              write(91,*) ct
2872 c              write(91,*) mult
2873               DO 150 IC = 1, mult
2874                  I = IA + IC
2875 c                 write(91,210) LB(i),r(1,i),r(2,i),r(3,i),
2876 c     1              p(1,i),p(2,i),p(3,i),e(i)
2877  150              continue
2878  1005      continue
2879            return
2880         endif
2881 c 210        format(i6,7(1x,f8.3))
2882         return
2883         end
2884
2885 c=======================================================================
2886         subroutine flowh0(NEVNT,idd)
2887 c
2888         dimension tsh(31)
2889         DOUBLE PRECISION  v2h,xnhadr,eth,v2h2,s2h
2890         DOUBLE PRECISION  v2hp,xnhadp,v2hsum,v2h2sm,
2891      1 v2havg(3),varv2h(3)
2892         COMMON /hflow/ v2h(30,3),xnhadr(30,3),eth(30,3),
2893      1 v2h2(30,3),s2h(30,3)
2894 cc      SAVE /hflow/
2895         common/ebe/v2hp(3),xnhadp(3),v2hsum(3),v2h2sm(3)
2896 cc      SAVE /ebe/
2897         common/input1/ MASSPR,MASSTA,ISEED,IAVOID,DT
2898 cc      SAVE /input1/
2899         COMMON /INPUT2/ ILAB, MANYB, NTMAX, ICOLL, INSYS, IPOT, MODE,
2900      &   IMOMEN, NFREQ, ICFLOW, ICRHO, ICOU, KPOTEN, KMUL
2901 cc      SAVE /INPUT2/
2902         common /lastt/itimeh,bimp
2903 cc      SAVE /lastt/
2904         SAVE
2905
2906 c     idd=0: initialization for flow analysis, called by artdri.f::
2907         if(idd.eq.0) then
2908            itimeh=0
2909 c
2910            do 1001 ii = 1, 31
2911               tsh(ii)=float(ii-1)
2912  1001      continue
2913 c
2914            do 1003 ii=1,30
2915               do 1002 iy=1,3
2916                  v2h(ii,iy)=0d0
2917                  xnhadr(ii,iy)=0d0
2918                  eth(ii,iy)=0d0
2919                  v2h2(ii,iy)=0d0
2920                  s2h(ii,iy)=0d0
2921  1002         continue
2922  1003      continue
2923            do 1004 iy=1,3
2924               v2hp(iy)=0d0
2925               xnhadp(iy)=0d0
2926               v2hsum(iy)=0d0
2927               v2h2sm(iy)=0d0
2928             if(iy.eq.1) then
2929                nunit=59
2930             elseif(iy.eq.2) then
2931                nunit=68
2932             else
2933                nunit=69
2934             endif
2935               write(nunit,*) '   tsh,   v2h,     v2h2,     s2h, '//
2936      1 ' eth,   xmulth'
2937  1004      continue
2938 c     idd=2: calculate average hadronic elliptic flow, called by artdri.f:
2939         else if(idd.eq.2) then
2940            do 100 ii=1, 30
2941               do 1005 iy=1,3
2942                  if(xnhadr(ii,iy).eq.0) then
2943                     xmulth=0.
2944                  elseif(xnhadr(ii,iy).gt.1) then
2945                     v2h(ii,iy)=v2h(ii,iy)/xnhadr(ii,iy)
2946                     eth(ii,iy)=eth(ii,iy)/dble(NEVNT)
2947                     v2h2(ii,iy)=dsqrt((v2h2(ii,iy)/xnhadr(ii,iy)
2948      1                    -v2h(ii,iy)**2)/(xnhadr(ii,iy)-1))
2949                     s2h(ii,iy)=s2h(ii,iy)/xnhadr(ii,iy)
2950                     xmulth=sngl(xnhadr(ii,iy)/NEVNT)
2951                  endif
2952              if(iy.eq.1) then
2953                 nunit=59
2954              elseif(iy.eq.2) then
2955                 nunit=68
2956              else
2957                 nunit=69
2958              endif
2959                  if(tsh(ii).le.(ntmax*dt))
2960      1                    write (nunit,200) tsh(ii),v2h(ii,iy),
2961      2      v2h2(ii,iy),s2h(ii,iy),eth(ii,iy),xmulth
2962  1005         continue
2963  100           continue
2964 c     event-by-event average & variance of v2h:
2965            do 1006 iy=1,3
2966               v2havg(iy)=v2hsum(iy)/dble(NEVNT)
2967       varv2h(iy)=dsqrt(v2h2sm(iy)/dble(NEVNT)-v2havg(iy)**2)
2968  1006 continue
2969            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): avg=', v2havg
2970            write(88, 240) 'EBE v2h,v2h(y2),v2h(y1): var=', varv2h
2971         endif
2972  200        format(2x,f5.2,3(2x,f7.4),2(2x,f9.2))
2973  240        format(a30,3(2x,f9.5))
2974         return
2975         end
2976
2977 c=======================================================================
2978 c     2/23/00 flow from all initial hadrons just before entering ARTMN:
2979         subroutine iniflw(NEVNT,idd)
2980         PARAMETER (MAXSTR=150001, MAXR=1)
2981         DOUBLE PRECISION  v2i,eti,xmulti,v2mi,s2mi,xmmult,
2982      1       v2bi,s2bi,xbmult
2983         COMMON /RUN/ NUM
2984 cc      SAVE /RUN/
2985         COMMON /ARERC1/MULTI1(MAXR)
2986 cc      SAVE /ARERC1/
2987         COMMON /ARPRC1/ITYP1(MAXSTR, MAXR),
2988      &     GX1(MAXSTR, MAXR), GY1(MAXSTR, MAXR), GZ1(MAXSTR, MAXR),
2989      &     FT1(MAXSTR, MAXR),
2990      &     PX1(MAXSTR, MAXR), PY1(MAXSTR, MAXR), PZ1(MAXSTR, MAXR),
2991      &     EE1(MAXSTR, MAXR), XM1(MAXSTR, MAXR)
2992 cc      SAVE /ARPRC1/
2993         COMMON/iflow/v2i,eti,xmulti,v2mi,s2mi,xmmult,v2bi,s2bi,xbmult
2994 cc      SAVE /iflow/
2995         SAVE
2996 c
2997         if(idd.eq.0) then
2998            v2i=0d0
2999            eti=0d0
3000            xmulti=0d0
3001            v2mi=0d0
3002            s2mi=0d0
3003            xmmult=0d0
3004            v2bi=0d0
3005            s2bi=0d0
3006            xbmult=0d0
3007         else if(idd.eq.1) then
3008            do 1002 J = 1, NUM
3009               do 1001 I = 1, MULTI1(J)
3010                  ITYP = ITYP1(I, J)
3011 c     all hadrons:
3012                  IF (ITYP .GT. -100 .AND. ITYP .LT. 100) GOTO 100
3013                  xmulti=xmulti+1.d0
3014                  PX = PX1(I, J)
3015                  PY = PY1(I, J)
3016                  XM = XM1(I, J)
3017                  pt2=px**2+py**2
3018                  xh=gx1(I,J)
3019                  yh=gy1(I,J)
3020                  xt2=xh**2+yh**2
3021                  if(pt2.gt.0) v2i=v2i+dble((px**2-py**2)/pt2)
3022                  eti=eti+dble(SQRT(PX ** 2 + PY ** 2 + XM ** 2))
3023 c     baryons only:
3024                  IF (ITYP .LT. -1000 .or. ITYP .GT. 1000) then
3025                     xbmult=xbmult+1.d0
3026                     if(pt2.gt.0) v2bi=v2bi+dble((px**2-py**2)/pt2)
3027                     if(xt2.gt.0) s2bi=s2bi+dble((xh**2-yh**2)/xt2)
3028 c     mesons only:
3029                  else
3030                     xmmult=xmmult+1.d0
3031                     if(pt2.gt.0) v2mi=v2mi+dble((px**2-py**2)/pt2)
3032                     if(xt2.gt.0) s2mi=s2mi+dble((xh**2-yh**2)/xt2)
3033                  endif
3034  100                 continue
3035  1001         continue
3036  1002      continue
3037         else if(idd.eq.2) then
3038            if(xmulti.ne.0) v2i=v2i/xmulti
3039            eti=eti/dble(NEVNT)
3040            xmulti=xmulti/dble(NEVNT)
3041            if(xmmult.ne.0) then
3042               v2mi=v2mi/xmmult
3043               s2mi=s2mi/xmmult
3044            endif
3045            xmmult=xmmult/dble(NEVNT)
3046            if(xbmult.ne.0) then
3047               v2bi=v2bi/xbmult
3048               s2bi=s2bi/xbmult
3049            endif
3050            xbmult=xbmult/dble(NEVNT)
3051         endif
3052 c
3053         return
3054         end
3055
3056 c=======================================================================
3057 c     2/25/00 dN/dt analysis for production (before ZPCMN)
3058 c     and freezeout (right after ZPCMN) for all partons.
3059         subroutine frztm(NEVNT,idd)
3060 c
3061         implicit double precision  (a-h, o-z)
3062         PARAMETER (MAXPTN=400001)
3063         dimension tsf(31)
3064         COMMON /PARA1/ MUL
3065 cc      SAVE /PARA1/
3066         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
3067      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
3068      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
3069 cc      SAVE /prec1/
3070         COMMON /prec2/GX5(MAXPTN),GY5(MAXPTN),GZ5(MAXPTN),FT5(MAXPTN),
3071      &       PX5(MAXPTN), PY5(MAXPTN), PZ5(MAXPTN), E5(MAXPTN),
3072      &       XMASS5(MAXPTN), ITYP5(MAXPTN)
3073 cc      SAVE /prec2/
3074         COMMON /frzout/ xnprod(30),etprod(30),xnfrz(30),etfrz(30),
3075      & dnprod(30),detpro(30),dnfrz(30),detfrz(30)
3076 cc      SAVE /frzout/
3077         SAVE
3078         data tsf/0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
3079      &       1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
3080      &       2  , 3,   4,   5,   6,   7,   8,   9,   10,  20,  30/
3081 c
3082         if(idd.eq.0) then
3083            do 1001 ii=1,30
3084               xnprod(ii)=0d0
3085               etprod(ii)=0d0
3086               xnfrz(ii)=0d0
3087               etfrz(ii)=0d0
3088               dnprod(ii)=0d0
3089               detpro(ii)=0d0
3090               dnfrz(ii)=0d0
3091               detfrz(ii)=0d0
3092  1001      continue
3093            OPEN (86, FILE = 'ana1/production.dat', STATUS = 'UNKNOWN')
3094            OPEN (87, FILE = 'ana1/freezeout.dat', STATUS = 'UNKNOWN')
3095         else if(idd.eq.1) then
3096            DO 100 ip = 1, MUL
3097               do 1002 ii=1,30
3098                  eth0=dSQRT(PX0(ip)**2+PY0(ip)**2+XMASS0(ip)**2)
3099                  eth2=dSQRT(PX5(ip)**2+PY5(ip)**2+XMASS5(ip)**2)
3100 c     total number and Et produced by time tsf(ii):
3101                  if (ft0(ip).lt.tsf(ii+1)) then
3102                     xnprod(ii)=xnprod(ii)+1d0
3103                     etprod(ii)=etprod(ii)+eth0
3104 c     number and Et produced from time tsf(ii) to tsf(ii+1):
3105                     if (ft0(ip).ge.tsf(ii)) then
3106                        dnprod(ii)=dnprod(ii)+1d0
3107                        detpro(ii)=detpro(ii)+eth0
3108                     endif
3109                  endif
3110 c     total number and Et freezed out by time tsf(ii):
3111                  if (FT5(ip).lt.tsf(ii+1)) then
3112                     xnfrz(ii)=xnfrz(ii)+1d0
3113                     etfrz(ii)=etfrz(ii)+eth2
3114 c     number and Et freezed out from time tsf(ii) to tsf(ii+1):
3115                     if (FT5(ip).ge.tsf(ii)) then
3116                        dnfrz(ii)=dnfrz(ii)+1d0
3117                        detfrz(ii)=detfrz(ii)+eth2
3118                     endif
3119                  endif
3120  1002         continue
3121  100           continue
3122         else if(idd.eq.2) then
3123            write (86,*) '       t,       np,       dnp/dt,      etp '//
3124      1 ' detp/dt'
3125            write (87,*) '       t,       nf,       dnf/dt,      etf '//
3126      1 ' detf/dt'
3127            do 1003 ii=1,30
3128               xnp=xnprod(ii)/dble(NEVNT)
3129               xnf=xnfrz(ii)/dble(NEVNT)
3130               etp=etprod(ii)/dble(NEVNT)
3131               etf=etfrz(ii)/dble(NEVNT)
3132               dxnp=dnprod(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
3133               dxnf=dnfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
3134               detp=detpro(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
3135               detf=detfrz(ii)/dble(NEVNT)/(tsf(ii+1)-tsf(ii))
3136               write (86, 200)
3137      1        tsf(ii+1),xnp,dxnp,etp,detp
3138               write (87, 200)
3139      1        tsf(ii+1),xnf,dxnf,etf,detf
3140  1003      continue
3141         endif
3142  200    format(2x,f9.2,4(2x,f10.2))
3143 c
3144         return
3145         end
3146
3147 c=======================================================================
3148 clin-6/2009 write out initial minijet information
3149 c     before propagating to its formation time:
3150 clin-2/2012:
3151 c        subroutine minijet_out(BB)
3152         subroutine minijet_out(BB,phiRP)
3153         PARAMETER (MAXSTR=150001)
3154         COMMON/HPARNT/HIPR1(100),IHPR2(50),HINT1(100),IHNT2(50)
3155         COMMON/hjcrdn/YP(3,300),YT(3,300)
3156         COMMON/HJJET1/NPJ(300),KFPJ(300,500),PJPX(300,500),
3157      &                PJPY(300,500),PJPZ(300,500),PJPE(300,500),
3158      &                PJPM(300,500),NTJ(300),KFTJ(300,500),
3159      &                PJTX(300,500),PJTY(300,500),PJTZ(300,500),
3160      &                PJTE(300,500),PJTM(300,500)
3161         COMMON/HJJET2/NSG,NJSG(MAXSTR),IASG(MAXSTR,3),K1SG(MAXSTR,100),
3162      &       K2SG(MAXSTR,100),PXSG(MAXSTR,100),PYSG(MAXSTR,100),
3163      &       PZSG(MAXSTR,100),PESG(MAXSTR,100),PMSG(MAXSTR,100)
3164         COMMON /AREVT/ IAEVT, IARUN, MISS
3165         common /para7/ ioscar,nsmbbbar,nsmmeson
3166         common/phidcy/iphidcy,pttrig,ntrig,maxmiss,ipi0dcy
3167         SAVE
3168         ntrig=0
3169         do I = 1, IHNT2(1)
3170            do J = 1, NPJ(I)
3171               pt=sqrt(PJPX(I,J)**2+PJPY(I,J)**2)
3172               if(pt.ge.pttrig) ntrig=ntrig+1
3173            enddo
3174         enddo
3175         do I = 1, IHNT2(3)
3176            do J = 1, NTJ(I)
3177               pt=sqrt(PJTX(I,J)**2+PJTY(I,J)**2)
3178               if(pt.ge.pttrig) ntrig=ntrig+1
3179            enddo
3180         enddo
3181         do I = 1, NSG
3182            do J = 1, NJSG(I)
3183               pt=sqrt(PXSG(I,J)**2+PYSG(I,J)**2)
3184               if(pt.ge.pttrig) ntrig=ntrig+1
3185            enddo
3186         enddo
3187 c     Require at least 1 initial minijet parton above the trigger Pt value:
3188         if(ntrig.eq.0) return
3189
3190 c.....transfer data from HIJING to ZPC
3191         if(ioscar.eq.3) write(96,*) IAEVT,MISS,IHNT2(1),IHNT2(3)
3192         DO 1008 I = 1, IHNT2(1)
3193            DO 1007 J = 1, NPJ(I)
3194               ityp=KFPJ(I,J)
3195 c     write out not only gluons:
3196 c              if(ityp.ne.21) goto 1007
3197 clin-2/2012:
3198 c              gx=YP(1,I)+0.5*BB
3199 c              gy=YP(2,I)
3200               gx=YP(1,I)+0.5*BB*cos(phiRP)
3201               gy=YP(2,I)+0.5*BB*sin(phiRP)
3202               gz=0.
3203               ft=0.
3204               px=PJPX(I,J)
3205               py=PJPY(I,J)
3206               pz=PJPZ(I,J)
3207               xmass=PJPM(I,J)
3208               if(ioscar.eq.3) then
3209                  if(amax1(abs(gx),abs(gy),
3210      1                abs(gz),abs(ft)).lt.9999) then
3211                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3212                  else
3213                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,1
3214                  endif
3215               endif
3216  1007      CONTINUE
3217  1008   CONTINUE
3218         DO 1010 I = 1, IHNT2(3)
3219            DO 1009 J = 1, NTJ(I)
3220               ityp=KFTJ(I,J)
3221 c              if(ityp.ne.21) goto 1009
3222 clin-2/2012:
3223 c              gx=YT(1,I)-0.5*BB
3224 c              gy=YT(2,I)
3225               gx=YT(1,I)-0.5*BB*cos(phiRP)
3226               gy=YT(2,I)-0.5*BB*sin(phiRP)
3227               gz=0.
3228               ft=0.
3229               px=PJTX(I,J)
3230               py=PJTY(I,J)
3231               pz=PJTZ(I,J)
3232               xmass=PJTM(I,J)
3233               if(ioscar.eq.3) then
3234                  if(amax1(abs(gx),abs(gy),
3235      1                abs(gz),abs(ft)).lt.9999) then
3236                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3237                  else
3238                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,2
3239                  endif
3240               endif
3241  1009      CONTINUE
3242  1010   CONTINUE
3243         DO 1012 I = 1, NSG
3244            DO 1011 J = 1, NJSG(I)
3245               ityp=K2SG(I,J)
3246 c              if(ityp.ne.21) goto 1011
3247               gx=0.5*(YP(1,IASG(I,1))+YT(1,IASG(I,2)))
3248               gy=0.5*(YP(2,IASG(I,1))+YT(2,IASG(I,2)))
3249               gz=0.
3250               ft=0.
3251               px=PXSG(I,J)
3252               py=PYSG(I,J)
3253               pz=PZSG(I,J)
3254               xmass=PMSG(I,J)
3255               if(ioscar.eq.3) then
3256                  if(amax1(abs(gx),abs(gy),
3257      1                abs(gz),abs(ft)).lt.9999) then
3258                     write(96,200) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3259                  else
3260                     write(96,201) ityp,px,py,pz,xmass,gx,gy,gz,ft,3
3261                  endif
3262               endif
3263  1011      CONTINUE
3264  1012   CONTINUE
3265  200  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,f8.2),2(2x,f2.0),2x,I2)
3266  201  format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,2(1x,e8.2),2(2x,f2.0),2x,I2)
3267 c
3268         return
3269         end
3270
3271 c=======================================================================
3272 clin-6/2009 embed back-to-back high-Pt quark/antiquark pair
3273 c     via embedding back-to-back high-Pt pion pair then melting the pion pair
3274 c     by generating the internal quark and antiquark momentum parallel to
3275 c      the pion momentum (in order to produce a high-Pt and a low Pt parton):
3276       subroutine embedHighPt
3277       PARAMETER (MAXSTR=150001,MAXR=1,pichmass=0.140,pi0mass=0.135,
3278      1     pi=3.1415926,nxymax=10001)
3279       common/embed/iembed,nsembd,pxqembd,pyqembd,xembd,yembd,
3280      1     psembd,tmaxembd,phidecomp
3281       COMMON/RNDF77/NSEED
3282       COMMON/HMAIN1/EATT,JATT,NATT,NT,NP,N0,N01,N10,N11
3283       COMMON/HMAIN2/KATT(MAXSTR,4),PATT(MAXSTR,4)
3284       COMMON /ARPRC/ ITYPAR(MAXSTR),
3285      &     GXAR(MAXSTR), GYAR(MAXSTR), GZAR(MAXSTR), FTAR(MAXSTR),
3286      &     PXAR(MAXSTR), PYAR(MAXSTR), PZAR(MAXSTR), PEAR(MAXSTR),
3287      &     XMAR(MAXSTR)
3288       common/anim/nevent,isoft,isflag,izpc
3289       COMMON /AREVT/ IAEVT, IARUN, MISS
3290       common/xyembed/nxyjet,xyjet(nxymax,2)
3291       SAVE
3292 c
3293       if(iembed.eq.1.or.iembed.eq.2) then
3294          xjet=xembd
3295          yjet=yembd
3296       elseif(iembed.eq.3.or.iembed.eq.4) then
3297          if(nevent.le.nxyjet) then
3298             read(97,*) xjet,yjet
3299          else
3300             ixy=mod(IAEVT,nxyjet)
3301             if(ixy.eq.0) ixy=nxyjet
3302             xjet=xyjet(ixy,1)
3303             yjet=xyjet(ixy,2)
3304          endif
3305       else
3306          return
3307       endif
3308 c
3309       ptq=sqrt(pxqembd**2+pyqembd**2)
3310       if(ptq.lt.(pichmass/2.)) then
3311          print *, 'Embedded quark transverse momentum is too small'
3312          stop
3313       endif
3314 c     Randomly embed u/ubar or d/dbar at high Pt:
3315       idqembd=1+int(2*RANART(NSEED))
3316 c     Flavor content for the charged pion that contains the leading quark:
3317       if(idqembd.eq.1) then
3318          idqsoft=-2
3319          idpi1=-211
3320       elseif(idqembd.eq.2) then
3321          idqsoft=-1
3322          idpi1=211
3323       else
3324          print *, 'Wrong quark flavor embedded'
3325          stop
3326       endif
3327 c     Caculate transverse momentum of the parent charged pion:
3328       xmq=ulmass(idqembd)
3329       xmqsoft=ulmass(idqsoft)
3330       ptpi=((pichmass**2+xmq**2-xmqsoft**2)*ptq
3331      1     -sqrt((xmq**2+ptq**2)*(pichmass**4
3332      2     -2.*pichmass**2*(xmq**2+xmqsoft**2)+(xmq**2-xmqsoft**2)**2)))
3333      3     /(2.*xmq**2)
3334       if(iembed.eq.1.or.iembed.eq.3) then
3335          pxpi1=ptpi*pxqembd/ptq
3336          pypi1=ptpi*pyqembd/ptq
3337          phidecomp=acos(pxqembd/ptq)
3338          if(pyqembd.lt.0) phidecomp=2.*pi-phidecomp
3339       else
3340          phidecomp=2.*pi*RANART(NSEED)
3341          pxpi1=ptpi*cos(phidecomp)
3342          pypi1=ptpi*sin(phidecomp)
3343       endif
3344 c     Embedded quark/antiquark are assumed to have pz=0:
3345       pzpi1=0.
3346 c     Insert the two parent charged pions,
3347 c     ipion=1 for the pion containing the leading quark,
3348 c     ipion=2 for the pion containing the leading antiquark of the same flavor:
3349       do ipion=1,2
3350          if(ipion.eq.1) then
3351             idpi=idpi1
3352             pxpi=pxpi1
3353             pypi=pypi1
3354             pzpi=pzpi1
3355          elseif(ipion.eq.2) then
3356             idpi=-idpi1
3357             pxpi=-pxpi1
3358             pypi=-pypi1
3359             pzpi=-pzpi1
3360          endif
3361          NATT=NATT+1
3362          KATT(NATT,1)=idpi
3363          KATT(NATT,2)=40
3364          KATT(NATT,3)=0
3365          PATT(NATT,1)=pxpi
3366          PATT(NATT,2)=pypi
3367          PATT(NATT,3)=pzpi
3368          PATT(NATT,4)=sqrt(pxpi**2+pypi**2+pzpi**2+pichmass**2)
3369          EATT=EATT+PATT(NATT,4)
3370          GXAR(NATT)=xjet
3371          GYAR(NATT)=yjet
3372          GZAR(NATT)=0.
3373          FTAR(NATT)=0.
3374          ITYPAR(NATT)=KATT(NATT,1)
3375          PXAR(NATT)=PATT(NATT,1)
3376          PYAR(NATT)=PATT(NATT,2)
3377          PZAR(NATT)=PATT(NATT,3)
3378          PEAR(NATT)=PATT(NATT,4)
3379          XMAR(NATT)=pichmass
3380       enddo
3381 c
3382 clin-8/2009
3383 c     Randomly embed a number of soft pions around each high-Pt quark in pair:
3384       if(nsembd.gt.0) then
3385          do ipion=1,2
3386             do ispion=1,nsembd
3387                idsart=3+int(3*RANART(NSEED))
3388                if(idsart.eq.3) then
3389                   pimass=pichmass
3390                   idpis=-211
3391                elseif(idsart.eq.4) then
3392                   pimass=pi0mass
3393                   idpis=111
3394                else
3395                   pimass=pichmass
3396                   idpis=211
3397                endif
3398                NATT=NATT+1
3399                KATT(NATT,1)=idpis
3400                KATT(NATT,2)=40
3401                KATT(NATT,3)=0
3402 c     theta: relative angle between soft pion & associated high-Pt q or qbar,
3403 c     generate theta and phi uniformly:
3404 c     Note: it is not generated uniformly in solid angle because that gives
3405 c     a valley at theta=0, unlike the jet-like correlation (a peak at theta=0).
3406                theta=tmaxembd*RANART(NSEED)
3407                phi=2.*pi*RANART(NSEED)
3408                pxspi=psembd*sin(theta)*cos(phi)
3409                pyspi=psembd*sin(theta)*sin(phi)
3410                pzspi=psembd*cos(theta)
3411                if(ipion.eq.1) then
3412                   call rotate(pxpi1,pypi1,pzpi1,pxspi,pyspi,pzspi)
3413                else
3414                   call rotate(-pxpi1,-pypi1,-pzpi1,pxspi,pyspi,pzspi)
3415                endif
3416 ctest off
3417 c               write(99,*) "2  ",pxspi,pyspi,pzspi
3418                PATT(NATT,1)=pxspi
3419                PATT(NATT,2)=pyspi
3420                PATT(NATT,3)=pzspi
3421                PATT(NATT,4)=sqrt(psembd**2+pimass**2)
3422                EATT=EATT+PATT(NATT,4)
3423                GXAR(NATT)=xjet
3424                GYAR(NATT)=yjet
3425                GZAR(NATT)=0.
3426                FTAR(NATT)=0.
3427                ITYPAR(NATT)=KATT(NATT,1)
3428                PXAR(NATT)=PATT(NATT,1)
3429                PYAR(NATT)=PATT(NATT,2)
3430                PZAR(NATT)=PATT(NATT,3)
3431                PEAR(NATT)=PATT(NATT,4)
3432                XMAR(NATT)=pimass
3433             enddo
3434          enddo
3435       endif
3436 clin-8/2009-end
3437 c
3438       return
3439       end