pyboei.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyboei.f

 
C*********************************************************************
 
C...PYBOEI
C...Modifies an event so as to approximately take into account
C...Bose-Einstein effects according to a simple phenomenological
C...parametrization.
 
      SUBROUTINE pyboei(NSAV)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pyint1/mint(400),vint(400)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pyint1/
C...Local arrays and data.
      dimension dps(4),kfbe(9),nbe(0:10),bei(100),bei3(100),
     &beiw(100),bei3w(100)
      DATA kfbe/211,-211,111,321,-321,130,310,221,331/
C...Statement function: squared invariant mass.
      sdip(i,j)=((p(i,4)+p(j,4))**2-(p(i,3)+p(j,3))**2-
     &(p(i,2)+p(j,2))**2-(p(i,1)+p(j,1))**2)
 
C...Boost event to overall CM frame. Calculate CM energy.
      IF((mstj(51).NE.1.AND.mstj(51).NE.2).OR.n-nsav.LE.1) RETURN
      DO 100 j=1,4
        dps(j)=0d0
  100 CONTINUE
      DO 120 i=1,n
        kfa=iabs(k(i,2))
        IF(k(i,1).LE.10.AND.((kfa.GT.10.AND.kfa.LE.20).OR.kfa.EQ.22)
     &  .AND.k(i,3).GT.0) THEN
          kfma=iabs(k(k(i,3),2))
          IF(kfma.GT.10.AND.kfma.LE.80) k(i,1)=-k(i,1)
        ENDIF
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120
        DO 110 j=1,4
          dps(j)=dps(j)+p(i,j)
  110   CONTINUE
  120 CONTINUE
      CALL pyrobo(0,0,0d0,0d0,-dps(1)/dps(4),-dps(2)/dps(4),
     &-dps(3)/dps(4))
      pecm=0d0
      DO 130 i=1,n
        IF(k(i,1).GE.1.AND.k(i,1).LE.10) pecm=pecm+p(i,4)
  130 CONTINUE
 
C...Check if we have separated strings
 
C...Reserve copy of particles by species at end of record.
      iwp=0
      iwn=0
      nbe(0)=n+mstu(3)
      nmax=nbe(0)
      smmin=pecm
      DO 190 ibe=1,min(10,mstj(52)+1)
        nbe(ibe)=nbe(ibe-1)
        DO 180 i=nsav+1,n
          IF(ibe.EQ.min(10,mstj(52)+1)) THEN
            DO 140 iibe=1,ibe-1
              IF(k(i,2).EQ.kfbe(iibe)) goto 180
  140       CONTINUE
          ELSE
            IF(k(i,2).NE.kfbe(ibe)) goto 180
          ENDIF
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 180
          IF(nbe(ibe).GE.mstu(4)-mstu(32)-5) THEN
            CALL pyerrm(11,'(PYBOEI:) no more memory left in PYJETS')
            RETURN
          ENDIF
          nbe(ibe)=nbe(ibe)+1
          nmax=nbe(ibe)
          k(nbe(ibe),1)=i
          k(nbe(ibe),2)=0
          k(nbe(ibe),3)=0
          k(nbe(ibe),4)=0
          k(nbe(ibe),5)=0
          p(nbe(ibe),1)=0.0d0
          p(nbe(ibe),2)=0.0d0
          p(nbe(ibe),3)=0.0d0
          p(nbe(ibe),4)=0.0d0
          p(nbe(ibe),5)=0.0d0
          smmin=min(smmin,p(i,5))
C...Check if particles comes from different W's or Z's
          IF((mstj(53).NE.0.OR.mstj(56).GT.0).AND.mint(32).EQ.0) THEN
            im=i
  150       IF(k(im,3).GT.0) THEN
              im=k(im,3)
              IF(abs(k(im,2)).NE.24.AND.k(im,2).NE.23) goto 150
              k(nbe(ibe),5)=im
              IF(iwp.EQ.0.AND.k(im,2).EQ.24) iwp=im
              IF(iwn.EQ.0.AND.k(im,2).EQ.-24) iwn=im
              IF(iwp.EQ.0.AND.k(im,2).EQ.23) iwp=im
              IF(iwn.EQ.0.AND.k(im,2).EQ.23.AND.im.NE.iwp) iwn=im
            ENDIF
          ENDIF
C...Check if particles comes from different strings.
          IF(parj(94).GT.0.0d0) THEN
            im=i
  160       IF(k(im,3).GT.0) THEN
              im=k(im,3)
              IF(k(im,2).NE.92.AND.k(im,2).NE.91) goto 160
              k(nbe(ibe),5)=im
            ENDIF
          ENDIF
          DO 170 j=1,3
            p(nbe(ibe),j)=0d0
            v(nbe(ibe),j)=0d0
  170     CONTINUE
          p(nbe(ibe),5)=-1.0d0
  180   CONTINUE
  190 CONTINUE
      IF(nbe(min(9,mstj(52)))-nbe(0).LE.1) goto 510
 
C...Calculate separation between W+ and W- or between two Z0's.
C...No separation if there has been re-connections.
      sigw=parj(93)
      IF(iwp.GT.0.AND.iwn.GT.0.AND.mstj(56).GT.0.AND.mint(32).EQ.0) THEN
        IF(k(iwp,2).EQ.23) THEN
          dmw=pmas(23,1)
          dgw=pmas(23,2)
        ELSE
          dmw=pmas(24,1)
          dgw=pmas(24,2)
        ENDIF
        dmp=p(iwp,5)
        dmn=p(iwn,5)
        taupd=dmp/sqrt((dmp**2-dmw**2)**2+(dgw*(dmp**2)/dmw)**2)
        taund=dmn/sqrt((dmn**2-dmw**2)**2+(dgw*(dmn**2)/dmw)**2)
        taup=-taupd*log(pyr(idum))
        taun=-taund*log(pyr(idum))
        dxp=taup*pyp(iwp,8)/dmp
        dxn=taun*pyp(iwn,8)/dmn
        dx=dxp+dxn
        sigw=1.0d0/(1.0d0/parj(93)+REAL(mstj(56))*dx)
        IF(parj(94).LT.0.0d0) sigw=1.0d0/(1.0d0/sigw-1.0d0/parj(94))
      ENDIF
 
C...Add separation between strings.
      IF(parj(94).GT.0.0d0) THEN
        sigw=1.0d0/(1.0d0/sigw+1.0d0/parj(94))
        iwp=-1
        iwn=-1
      ENDIF
 
      IF(mstj(57).EQ.1.AND.mstj(54).LT.0) THEN
        DO 220 ibe=1,min(9,mstj(52))
          DO 210 i1m=nbe(ibe-1)+1,nbe(ibe)
            q2min=pecm**2
            i1=k(i1m,1)
            DO 200 i2m=nbe(ibe-1)+1,nbe(ibe)
              IF(i2m.EQ.i1m) goto 200
              i2=k(i2m,1)
              q2=(p(i1,4)+p(i2,4))**2-(p(i1,1)+p(i2,1))**2-
     &        (p(i1,2)+p(i2,2))**2-(p(i1,3)+p(i2,3))**2-
     &        (p(i1,5)+p(i2,5))**2
              IF(q2.GT.0.0d0.AND.q2.LT.q2min) THEN
                q2min=q2
              ENDIF
  200       CONTINUE
            p(i1m,5)=q2min
  210     CONTINUE
  220   CONTINUE
      ENDIF
 
C...Tabulate integral for subsequent momentum shift.
      DO 400 ibe=1,min(9,mstj(52))
        IF(ibe.NE.1.AND.ibe.NE.4.AND.ibe.LE.7) goto 270
        IF(ibe.EQ.1.AND.max(nbe(1)-nbe(0),nbe(2)-nbe(1),nbe(3)-nbe(2))
     &  .LE.1) goto 270
        IF(ibe.EQ.4.AND.max(nbe(4)-nbe(3),nbe(5)-nbe(4),nbe(6)-nbe(5),
     &  nbe(7)-nbe(6)).LE.1) goto 270
        IF(ibe.GE.8.AND.nbe(ibe)-nbe(ibe-1).LE.1) goto 270
        IF(ibe.EQ.1) pmhq=2d0*pymass(211)
        IF(ibe.EQ.4) pmhq=2d0*pymass(321)
        IF(ibe.EQ.8) pmhq=2d0*pymass(221)
        IF(ibe.EQ.9) pmhq=2d0*pymass(331)
        qdel=0.1d0*min(pmhq,parj(93))
        qdel3=0.1d0*min(pmhq,parj(93)*3.0d0)
        qdelw=0.1d0*min(pmhq,sigw)
        qdel3w=0.1d0*min(pmhq,sigw*3.0d0)
        IF(mstj(51).EQ.1) THEN
          nbin=min(100,nint(9d0*parj(93)/qdel))
          nbin3=min(100,nint(27d0*parj(93)/qdel3))
          nbinw=min(100,nint(9d0*sigw/qdelw))
          nbin3w=min(100,nint(27d0*sigw/qdel3w))
          beex=exp(0.5d0*qdel/parj(93))
          beex3=exp(0.5d0*qdel3/(3.0d0*parj(93)))
          beexw=exp(0.5d0*qdelw/sigw)
          beex3w=exp(0.5d0*qdel3w/(3.0d0*sigw))
          bert=exp(-qdel/parj(93))
          bert3=exp(-qdel3/(3.0d0*parj(93)))
          bertw=exp(-qdelw/sigw)
          bert3w=exp(-qdel3w/(3.0d0*sigw))
        ELSE
          nbin=min(100,nint(3d0*parj(93)/qdel))
          nbin3=min(100,nint(9d0*parj(93)/qdel3))
          nbinw=min(100,nint(3d0*sigw/qdelw))
          nbin3w=min(100,nint(9d0*sigw/qdel3w))
        ENDIF
        DO 230 ibin=1,nbin
          qbin=qdel*(ibin-0.5d0)
          bei(ibin)=qdel*(qbin**2+qdel**2/12d0)/sqrt(qbin**2+pmhq**2)
          IF(mstj(51).EQ.1) THEN
            beex=beex*bert
            bei(ibin)=bei(ibin)*beex
          ELSE
            bei(ibin)=bei(ibin)*exp(-(qbin/parj(93))**2)
          ENDIF
          IF(ibin.GE.2) bei(ibin)=bei(ibin)+bei(ibin-1)
  230   CONTINUE
        DO 240 ibin=1,nbin3
          qbin=qdel3*(ibin-0.5d0)
          bei3(ibin)=qdel3*(qbin**2+qdel3**2/12d0)/sqrt(qbin**2+pmhq**2)
          IF(mstj(51).EQ.1) THEN
            beex3=beex3*bert3
            bei3(ibin)=bei3(ibin)*beex3
          ELSE
            bei3(ibin)=bei3(ibin)*exp(-(qbin/(3.0d0*parj(93)))**2)
          ENDIF
          IF(ibin.GE.2) bei3(ibin)=bei3(ibin)+bei3(ibin-1)
  240   CONTINUE
        DO 250 ibin=1,nbinw
          qbin=qdelw*(ibin-0.5d0)
          beiw(ibin)=qdelw*(qbin**2+qdelw**2/12d0)/sqrt(qbin**2+pmhq**2)
          IF(mstj(51).EQ.1) THEN
            beexw=beexw*bertw
            beiw(ibin)=beiw(ibin)*beexw
          ELSE
            beiw(ibin)=beiw(ibin)*exp(-(qbin/sigw)**2)
          ENDIF
          IF(ibin.GE.2) beiw(ibin)=beiw(ibin)+beiw(ibin-1)
  250   CONTINUE
        DO 260 ibin=1,nbin3w
          qbin=qdel3w*(ibin-0.5d0)
          bei3w(ibin)=qdel3w*(qbin**2+qdel3w**2/12d0)/
     &    sqrt(qbin**2+pmhq**2)
          IF(mstj(51).EQ.1) THEN
            beex3w=beex3w*bert3w
            bei3w(ibin)=bei3w(ibin)*beex3w
          ELSE
            bei3w(ibin)=bei3w(ibin)*exp(-(qbin/(3.0d0*sigw))**2)
          ENDIF
          IF(ibin.GE.2) bei3w(ibin)=bei3w(ibin)+bei3w(ibin-1)
  260   CONTINUE
 
C...Loop through particle pairs and find old relative momentum.
  270   DO 390 i1m=nbe(ibe-1)+1,nbe(ibe)-1
          i1=k(i1m,1)
          DO 380 i2m=i1m+1,nbe(ibe)
            IF(mstj(53).EQ.1.AND.k(i1m,5).NE.k(i2m,5)) goto 380
            IF(mstj(53).EQ.2.AND.k(i1m,5).EQ.k(i2m,5)) goto 380
            i2=k(i2m,1)
            q2old=(p(i1,4)+p(i2,4))**2-(p(i1,1)+p(i2,1))**2-(p(i1,2)+
     &      p(i2,2))**2-(p(i1,3)+p(i2,3))**2-(p(i1,5)+p(i2,5))**2
            IF(q2old.LE.0.0d0) goto 380
            qold=sqrt(q2old)
 
C...Calculate new relative momentum.
            qmov=0.0d0
            qmov3=0.0d0
            qmovw=0.0d0
            qmov3w=0.0d0
            IF(qold.LT.1d-3*qdel) THEN
              goto 280
            ELSEIF(qold.LE.qdel) THEN
              qmov=qold/3d0
            ELSEIF(qold.LT.(nbin-0.1d0)*qdel) THEN
              rbin=qold/qdel
              ibin=rbin
              rinp=(rbin**3-ibin**3)/(3*ibin*(ibin+1)+1)
              qmov=(bei(ibin)+rinp*(bei(ibin+1)-bei(ibin)))*
     &        sqrt(q2old+pmhq**2)/q2old
            ELSE
              qmov=bei(nbin)*sqrt(q2old+pmhq**2)/q2old
            ENDIF
  280       q2new=q2old*(qold/(qold+3d0*parj(92)*qmov))**(2d0/3d0)
            IF(qold.LT.1d-3*qdel3) THEN
              goto 290
            ELSEIF(qold.LE.qdel3) THEN
              qmov3=qold/3d0
            ELSEIF(qold.LT.(nbin3-0.1d0)*qdel3) THEN
              rbin3=qold/qdel3
              ibin3=rbin3
              rinp3=(rbin3**3-ibin3**3)/(3*ibin3*(ibin3+1)+1)
              qmov3=(bei3(ibin3)+rinp3*(bei3(ibin3+1)-bei3(ibin3)))*
     &        sqrt(q2old+pmhq**2)/q2old
            ELSE
              qmov3=bei3(nbin3)*sqrt(q2old+pmhq**2)/q2old
            ENDIF
  290       q2new3=q2old*(qold/(qold+3d0*parj(92)*qmov3))**(2d0/3d0)
            rscale=1.0d0
            IF(mstj(54).EQ.2)
     &      rscale=1.0d0-exp(-(qold/(2d0*parj(93)))**2)
            IF((iwp.NE.-1.AND.mstj(56).LE.0).OR.iwp.EQ.0.OR.iwn.EQ.0.OR.
     &      k(i1m,5).EQ.k(i2m,5)) goto 320
 
            IF(qold.LT.1d-3*qdelw) THEN
              goto 300
            ELSEIF(qold.LE.qdelw) THEN
              qmovw=qold/3d0
            ELSEIF(qold.LT.(nbinw-0.1d0)*qdelw) THEN
              rbinw=qold/qdelw
              ibinw=rbinw
              rinpw=(rbinw**3-ibinw**3)/(3*ibinw*(ibinw+1)+1)
              qmovw=(beiw(ibinw)+rinpw*(beiw(ibinw+1)-beiw(ibinw)))*
     &        sqrt(q2old+pmhq**2)/q2old
            ELSE
              qmovw=beiw(nbinw)*sqrt(q2old+pmhq**2)/q2old
            ENDIF
  300       q2new=q2old*(qold/(qold+3d0*parj(92)*qmovw))**(2d0/3d0)
            IF(qold.LT.1d-3*qdel3w) THEN
              goto 310
            ELSEIF(qold.LE.qdel3w) THEN
              qmov3w=qold/3d0
            ELSEIF(qold.LT.(nbin3w-0.1d0)*qdel3w) THEN
              rbin3w=qold/qdel3w
              ibin3w=rbin3w
              rinp3w=(rbin3w**3-ibin3w**3)/(3*ibin3w*(ibin3w+1)+1)
              qmov3w=(bei3w(ibin3w)+rinp3w*(bei3w(ibin3w+1)-
     &        bei3w(ibin3w)))*sqrt(q2old+pmhq**2)/q2old
            ELSE
              qmov3w=bei3w(nbin3w)*sqrt(q2old+pmhq**2)/q2old
            ENDIF
  310       q2new3=q2old*(qold/(qold+3d0*parj(92)*qmov3w))**(2d0/3d0)
            IF(mstj(54).EQ.2)
     &      rscale=1.0d0-exp(-(qold/(2d0*sigw))**2)
 
  320       CALL pybesq(i1,i2,nmax,q2old,q2new)
            DO 330 j=1,3
              p(i1m,j)=p(i1m,j)+p(nmax+1,j)
              p(i2m,j)=p(i2m,j)+p(nmax+2,j)
  330       CONTINUE
            IF(mstj(54).GE.1) THEN
              CALL pybesq(i1,i2,nmax,q2old,q2new3)
              DO 340 j=1,3
                v(i1m,j)=v(i1m,j)+p(nmax+1,j)*rscale
                v(i2m,j)=v(i2m,j)+p(nmax+2,j)*rscale
  340         CONTINUE
            ELSEIF(mstj(54).LE.-1) THEN
              edel=p(i1,4)+p(i2,4)-
     &        sqrt(max(q2new-q2old+(p(i1,4)+p(i2,4))**2,0.0d0))
              a2=(p(i1,1)-p(i2,1))**2+(p(i1,2)-p(i2,2))**2+
     &        (p(i1,3)-p(i2,3))**2
              wmax=-1.0d20
              mi3=0
              mi4=0
              s12=sdip(i1,i2)
              sm1=(p(i1,5)+smmin)**2
              DO 360 i3m=nbe(0)+1,nbe(min(10,mstj(52)+1))
                IF(i3m.EQ.i1m.OR.i3m.EQ.i2m) goto 360
                IF(mstj(53).EQ.1.AND.k(i3m,5).NE.k(i1m,5)) goto 360
                IF(mstj(53).EQ.-2.AND.k(i1m,5).EQ.k(i2m,5).AND.
     &          k(i3m,5).NE.k(i1m,5)) goto 360
                i3=k(i3m,1)
                IF(k(i3,2).EQ.k(i1,2)) goto 360
                s13=sdip(i1,i3)
                s23=sdip(i2,i3)
                sm3=(p(i3,5)+smmin)**2
                IF(mstj(54).EQ.-2) THEN
                  wi=(min(s12*sm3,s13*min(sm1,sm3),
     &            s23*min(sm1,sm3))*sm1)
                ELSE
                  wi=((p(i1,4)+p(i2,4)+p(i3,4))**2-
     &            (p(i1,3)+p(i2,3)+p(i3,3))**2-
     &            (p(i1,2)+p(i2,2)+p(i3,2))**2-
     &            (p(i1,1)+p(i2,1)+p(i3,1))**2)
                ENDIF
                IF(mstj(57).EQ.1.AND.p(i3m,5).GT.0) THEN
                  IF (wmax*wi.GE.(1.0d0-exp(-p(i3m,5)/(parj(93)**2))))
     &                 goto 360
                ELSE
                  IF(wmax*wi.GE.1.0) goto 360
                ENDIF
                DO 350 i4m=i3m+1,nbe(min(10,mstj(52)+1))
                  IF(i4m.EQ.i1m.OR.i4m.EQ.i2m) goto 350
                  IF(mstj(53).EQ.1.AND.k(i4m,5).NE.k(i1m,5)) goto 350
                  IF(mstj(53).EQ.-2.AND.k(i1m,5).EQ.k(i2m,5).AND.
     &            k(i4m,5).NE.k(i1m,5)) goto 350
                  i4=k(i4m,1)
                  IF(k(i3,2).EQ.k(i4,2).OR.k(i4,2).EQ.k(i1,2))
     &            goto 350
                  IF((p(i3,4)+p(i4,4)+edel)**2.LT.
     &            (p(i3,1)+p(i4,1))**2+(p(i3,2)+p(i4,2))**2+
     &            (p(i3,3)+p(i4,3))**2+(p(i3,5)+p(i4,5))**2)
     &            goto 350
                  IF(mstj(54).EQ.-2) THEN
                    s14=sdip(i1,i4)
                    s24=sdip(i2,i4)
                    s34=sdip(i3,i4)
                    w=s12*min(min(s23,s24),min(s13,s14))*s34
                    w=min(w,s13*min(min(s23,s34),s12)*s24)
                    w=min(w,s14*min(min(s24,s34),s12)*s23)
                    w=min(w,min(s23,s24)*s13*s14)
                    w=1.0d0/w
                  ELSE
C...weight=1-cos(theta)/mtot2
                    s1234=(p(i1,4)+p(i2,4)+p(i3,4)+p(i4,4))**2-
     &              (p(i1,3)+p(i2,3)+p(i3,3)+p(i4,3))**2-
     &              (p(i1,2)+p(i2,2)+p(i3,2)+p(i4,2))**2-
     &              (p(i1,1)+p(i2,1)+p(i3,1)+p(i4,1))**2
                    w=1.0d0/s1234
                    IF(w.LE.wmax) goto 350
                  ENDIF
                  IF(mstj(57).EQ.1.AND.p(i3m,5).GT.0)
     &            w=w*(1.0d0-exp(-p(i3m,5)/(parj(93)**2)))
                  IF(mstj(57).EQ.1.AND.p(i4m,5).GT.0)
     &            w=w*(1.0d0-exp(-p(i4m,5)/(parj(93)**2)))
                  IF(w.LE.wmax) goto 350
                  mi3=i3m
                  mi4=i4m
                  wmax=w
  350           CONTINUE
  360         CONTINUE
              IF(mi4.EQ.0) goto 380
              i3=k(mi3,1)
              i4=k(mi4,1)
              eold=p(i3,4)+p(i4,4)
              enew=eold+edel
              p2=(p(i3,1)+p(i4,1))**2+(p(i3,2)+p(i4,2))**2+
     &        (p(i3,3)+p(i4,3))**2
              q2newp=max(0.0d0,enew**2-p2-(p(i3,5)+p(i4,5))**2)
              q2oldp=max(0.0d0,eold**2-p2-(p(i3,5)+p(i4,5))**2)
              CALL pybesq(i3,i4,nmax,q2oldp,q2newp)
              DO 370 j=1,3
                v(mi3,j)=v(mi3,j)+p(nmax+1,j)
                v(mi4,j)=v(mi4,j)+p(nmax+2,j)
  370         CONTINUE
            ENDIF
  380     CONTINUE
  390   CONTINUE
  400 CONTINUE
 
C...Shift momenta and recalculate energies.
      esump=0.0d0
      esum=0.0d0
      prod=0.0d0
      DO 430 im=nbe(0)+1,nbe(min(10,mstj(52)+1))
        i=k(im,1)
        esump=esump+p(i,4)
        DO 410 j=1,3
          p(i,j)=p(i,j)+p(im,j)
  410   CONTINUE
        p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
        esum=esum+p(i,4)
        DO 420 j=1,3
          prod=prod+v(im,j)*p(i,j)/p(i,4)
  420   CONTINUE
  430 CONTINUE
 
      parj(96)=0.0d0
      IF(mstj(54).NE.0.AND.prod.NE.0.0d0) THEN
  440   alpha=(esump-esum)/prod
        parj(96)=parj(96)+alpha
        prod=0.0d0
        esum=0.0d0
        DO 470 im=nbe(0)+1,nbe(min(10,mstj(52)+1))
          i=k(im,1)
          DO 450 j=1,3
            p(i,j)=p(i,j)+alpha*v(im,j)
  450     CONTINUE
          p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
          esum=esum+p(i,4)
          DO 460 j=1,3
            prod=prod+v(im,j)*p(i,j)/p(i,4)
  460     CONTINUE
  470   CONTINUE
        IF(prod.NE.0.0d0.AND.abs(esump-esum)/pecm.GT.0.00001d0)
     &  goto 440
      ENDIF
 
C...Rescale all momenta for energy conservation.
      pes=0d0
      pqs=0d0
      DO 480 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 480
        pes=pes+p(i,4)
        pqs=pqs+p(i,5)**2/p(i,4)
  480 CONTINUE
      parj(95)=pes-pecm
      fac=(pecm-pqs)/(pes-pqs)
      DO 500 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 500
        DO 490 j=1,3
          p(i,j)=fac*p(i,j)
  490   CONTINUE
        p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
  500 CONTINUE
 
C...Boost back to correct reference frame.
  510 CALL pyrobo(0,0,0d0,0d0,dps(1)/dps(4),dps(2)/dps(4),dps(3)/dps(4))
      DO 520 i=1,n
        IF(k(i,1).LT.0) k(i,1)=-k(i,1)
  520 CONTINUE
 
      RETURN
      END