pyremn.f

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

 
C*********************************************************************
 
C...PYREMN
C...Adds on target remnants (one or two from each side) and
C...includes primordial kT for hadron beams.
 
      SUBROUTINE pyremn(IPU1,IPU2)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
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/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/
C...Local arrays.
      dimension kflch(2),kflsp(2),chi(2),pms(0:6),is(2),isn(2),robo(5),
     &psys(0:2,5),pmin(0:2),qold(4),qnew(4),dbe(3),psum(4)
 
C...Find event type and remaining energy.
      isub=mint(1)
      ns=n
      IF(mint(50).EQ.0.OR.mod(mstp(81),10).LE.0) THEN
        vint(143)=1d0-vint(141)
        vint(144)=1d0-vint(142)
      ENDIF
 
C...Define initial partons.
      ntry=0
  100 ntry=ntry+1
      DO 130 jt=1,2
        i=mint(83)+jt+2
        IF(jt.EQ.1) ipu=ipu1
        IF(jt.EQ.2) ipu=ipu2
        k(i,1)=21
        k(i,2)=k(ipu,2)
        k(i,3)=i-2
        pms(jt)=0d0
        vint(156+jt)=0d0
        vint(158+jt)=0d0
        IF(mint(47).EQ.1) THEN
          DO 110 j=1,5
            p(i,j)=p(i-2,j)
  110     CONTINUE
        ELSEIF(isub.EQ.95) THEN
          k(i,2)=21
        ELSE
          p(i,5)=p(ipu,5)
 
C...No primordial kT, or chosen according to truncated Gaussian or
C...exponential, or (for photon) predetermined or power law.
  120     IF(mint(40+jt).EQ.2.AND.mint(10+jt).NE.22) THEN
            IF(mstp(91).LE.0) THEN
              pt=0d0
            ELSEIF(mstp(91).EQ.1) THEN
              pt=parp(91)*sqrt(-log(pyr(0)))
            ELSE
              rpt1=pyr(0)
              rpt2=pyr(0)
              pt=-parp(92)*log(rpt1*rpt2)
            ENDIF
            IF(pt.GT.parp(93)) goto 120
          ELSEIF(mint(106+jt).EQ.3) THEN
            pta=sqrt(vint(282+jt))
            ptb=0d0
            IF(mstp(66).EQ.5.AND.mstp(93).EQ.1) THEN
              ptb=parp(99)*sqrt(-log(pyr(0)))
            ELSEIF(mstp(66).EQ.5.AND.mstp(93).EQ.2) THEN
              rpt1=pyr(0)
              rpt2=pyr(0)
              ptb=-parp(99)*log(rpt1*rpt2)
            ENDIF
            IF(ptb.GT.parp(100)) goto 120
            pt=sqrt(pta**2+ptb**2+2d0*pta*ptb*cos(paru(2)*pyr(0)))
            pt=pt*0.8d0**mint(57)
            IF(ntry.GT.10) pt=pt*0.8d0**(ntry-10)
          ELSEIF(iabs(mint(14+jt)).LE.8.OR.mint(14+jt).EQ.21) THEN
            IF(mstp(93).LE.0) THEN
              pt=0d0
            ELSEIF(mstp(93).EQ.1) THEN
              pt=parp(99)*sqrt(-log(pyr(0)))
            ELSEIF(mstp(93).EQ.2) THEN
              rpt1=pyr(0)
              rpt2=pyr(0)
              pt=-parp(99)*log(rpt1*rpt2)
            ELSEIF(mstp(93).EQ.3) THEN
              ha=parp(99)**2
              hb=parp(100)**2
              pt=sqrt(max(0d0,ha*(ha+hb)/(ha+hb-pyr(0)*hb)-ha))
            ELSE
              ha=parp(99)**2
              hb=parp(100)**2
              IF(mstp(93).EQ.5) hb=min(vint(48),parp(100)**2)
              pt=sqrt(max(0d0,ha*((ha+hb)/ha)**pyr(0)-ha))
            ENDIF
            IF(pt.GT.parp(100)) goto 120
          ELSE
            pt=0d0
          ENDIF
          vint(156+jt)=pt
          phi=paru(2)*pyr(0)
          p(i,1)=pt*cos(phi)
          p(i,2)=pt*sin(phi)
          pms(jt)=p(i,5)**2+p(i,1)**2+p(i,2)**2
        ENDIF
  130 CONTINUE
      IF(mint(47).EQ.1) RETURN
 
C...Kinematics construction for initial partons.
      i1=mint(83)+3
      i2=mint(83)+4
      IF(isub.EQ.95) THEN
        shs=0d0
        shr=0d0
      ELSE
        shs=vint(141)*vint(142)*vint(2)+(p(i1,1)+p(i2,1))**2+
     &  (p(i1,2)+p(i2,2))**2
        shr=sqrt(max(0d0,shs))
        IF((shs-pms(1)-pms(2))**2-4d0*pms(1)*pms(2).LE.0d0) goto 100
        p(i1,4)=0.5d0*(shr+(pms(1)-pms(2))/shr)
        p(i1,3)=sqrt(max(0d0,p(i1,4)**2-pms(1)))
        p(i2,4)=shr-p(i1,4)
        p(i2,3)=-p(i1,3)
 
C...Transform partons to overall CM-frame.
        robo(3)=(p(i1,1)+p(i2,1))/shr
        robo(4)=(p(i1,2)+p(i2,2))/shr
        CALL pyrobo(i1,i2,0d0,0d0,-robo(3),-robo(4),0d0)
        robo(2)=pyangl(p(i1,1),p(i1,2))
        CALL pyrobo(i1,i2,0d0,-robo(2),0d0,0d0,0d0)
        robo(1)=pyangl(p(i1,3),p(i1,1))
        CALL pyrobo(i1,i2,-robo(1),0d0,0d0,0d0,0d0)
        CALL pyrobo(i2+1,mint(52),0d0,-robo(2),0d0,0d0,0d0)
        CALL pyrobo(i1,mint(52),robo(1),robo(2),robo(3),robo(4),0d0)
        robo(5)=(vint(141)-vint(142))/(vint(141)+vint(142))
        CALL pyrobo(i1,mint(52),0d0,0d0,0d0,0d0,robo(5))
      ENDIF
 
C...Optionally fix up x and Q2 definitions for leptoproduction.
      idisxq=0
      IF((mint(43).EQ.2.OR.mint(43).EQ.3).AND.((isub.EQ.10.AND.
     &mstp(23).GE.1).OR.(isub.EQ.83.AND.mstp(23).GE.2))) idisxq=1
      IF(idisxq.EQ.1) THEN
 
C...Find where incoming and outgoing leptons/partons are sitting.
        lesd=1
        IF(mint(42).EQ.1) lesd=2
        lpin=mint(83)+3-lesd
        lein=mint(84)+lesd
        lqin=mint(84)+3-lesd
        leout=mint(84)+2+lesd
        lqout=mint(84)+5-lesd
        IF(k(lein,3).GT.lein) lein=k(lein,3)
        IF(k(lqin,3).GT.lqin) lqin=k(lqin,3)
        lscms=0
        DO 140 i=mint(84)+5,n
          IF(k(i,2).EQ.94) THEN
            lscms=i
            leout=i+lesd
            lqout=i+3-lesd
          ENDIF
  140   CONTINUE
        lqbg=ipu1
        IF(lesd.EQ.1) lqbg=ipu2
 
C...Calculate actual and wanted momentum transfer.
        xnom=vint(43-lesd)
        q2nom=-vint(45)
        hpk=2d0*(p(lpin,4)*p(lein,4)-p(lpin,1)*p(lein,1)-
     &  p(lpin,2)*p(lein,2)-p(lpin,3)*p(lein,3))*
     &  (p(mint(83)+lesd,4)*vint(40+lesd)/p(lein,4))
        hpt2=max(0d0,q2nom*(1d0-q2nom/(xnom*hpk)))
        fac=sqrt(hpt2/(p(leout,1)**2+p(leout,2)**2))
        p(n+1,1)=fac*p(leout,1)
        p(n+1,2)=fac*p(leout,2)
        p(n+1,3)=0.25d0*((hpk-q2nom/xnom)/p(lpin,4)-
     &  q2nom/(p(mint(83)+lesd,4)*vint(40+lesd)))*(-1)**(lesd+1)
        p(n+1,4)=sqrt(p(leout,5)**2+p(n+1,1)**2+p(n+1,2)**2+
     &  p(n+1,3)**2)
        DO 150 j=1,4
          qold(j)=p(lein,j)-p(leout,j)
          qnew(j)=p(lein,j)-p(n+1,j)
  150   CONTINUE
 
C...Boost outgoing electron and daughters.
        IF(lscms.EQ.0) THEN
          DO 160 j=1,4
            p(leout,j)=p(n+1,j)
  160     CONTINUE
        ELSE
          DO 170 j=1,3
            p(n+2,j)=(p(n+1,j)-p(leout,j))/(p(n+1,4)+p(leout,4))
  170     CONTINUE
          pinv=2d0/(1d0+p(n+2,1)**2+p(n+2,2)**2+p(n+2,3)**2)
          DO 180 j=1,3
            dbe(j)=pinv*p(n+2,j)
  180     CONTINUE
          DO 200 i=lscms+1,n
            iorig=i
  190       iorig=k(iorig,3)
            IF(iorig.GT.leout) goto 190
            IF(i.EQ.leout.OR.iorig.EQ.leout)
     &      CALL pyrobo(i,i,0d0,0d0,dbe(1),dbe(2),dbe(3))
  200     CONTINUE
        ENDIF
 
C...Copy shower initiator and all outgoing partons.
        ncop=n+1
        k(ncop,3)=lqbg
        DO 210 j=1,5
          p(ncop,j)=p(lqbg,j)
  210   CONTINUE
        DO 240 i=mint(84)+1,n
          icop=0
          IF(k(i,1).GT.10) goto 240
          IF(i.EQ.lqbg.OR.i.EQ.lqout) THEN
            icop=i
          ELSE
            iorig=i
  220       iorig=k(iorig,3)
            IF(iorig.EQ.lqbg.OR.iorig.EQ.lqout) THEN
              icop=iorig
            ELSEIF(iorig.GT.mint(84).AND.iorig.LE.n) THEN
              goto 220
            ENDIF
          ENDIF
          IF(icop.NE.0) THEN
            ncop=ncop+1
            k(ncop,3)=i
            DO 230 j=1,5
              p(ncop,j)=p(i,j)
  230       CONTINUE
          ENDIF
  240   CONTINUE
 
C...Calculate relative rescaling factors.
        slc=3-2*lesd
        plcsum=0d0
        DO 250 i=n+2,ncop
          plcsum=plcsum+(p(i,4)+slc*p(i,3))
  250   CONTINUE
        DO 260 i=n+2,ncop
          v(i,1)=(p(i,4)+slc*p(i,3))/plcsum
  260   CONTINUE
 
C...Transfer extra three-momentum of current.
        DO 280 i=n+2,ncop
          DO 270 j=1,3
            p(i,j)=p(i,j)+v(i,1)*(qnew(j)-qold(j))
  270     CONTINUE
          p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
  280   CONTINUE
 
C...Iterate change of initiator momentum to get energy right.
        iter=0
  290   iter=iter+1
        peex=-p(n+1,4)-qnew(4)
        pemv=-p(n+1,3)/p(n+1,4)
        DO 300 i=n+2,ncop
          peex=peex+p(i,4)
          pemv=pemv+v(i,1)*p(i,3)/p(i,4)
  300   CONTINUE
        IF(abs(pemv).LT.1d-10) THEN
          mint(51)=1
          mint(57)=mint(57)+1
          RETURN
        ENDIF
        pzch=-peex/pemv
        p(n+1,3)=p(n+1,3)+pzch
        p(n+1,4)=sqrt(p(n+1,5)**2+p(n+1,1)**2+p(n+1,2)**2+p(n+1,3)**2)
        DO 310 i=n+2,ncop
          p(i,3)=p(i,3)+v(i,1)*pzch
          p(i,4)=sqrt(p(i,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
  310   CONTINUE
        IF(iter.LT.10.AND.abs(peex).GT.1d-6*p(n+1,4)) goto 290
 
C...Modify momenta in event record.
        hbe=2d0*(p(n+1,4)+p(lqbg,4))*(p(n+1,3)-p(lqbg,3))/
     &  ((p(n+1,4)+p(lqbg,4))**2+(p(n+1,3)-p(lqbg,3))**2)
        IF(abs(hbe).GE.1d0) THEN
          mint(51)=1
          mint(57)=mint(57)+1
          RETURN
        ENDIF
        i=mint(83)+5-lesd
        CALL pyrobo(i,i,0d0,0d0,0d0,0d0,hbe)
        DO 330 i=n+1,ncop
          icop=k(i,3)
          DO 320 j=1,4
            p(icop,j)=p(i,j)
  320     CONTINUE
  330   CONTINUE
      ENDIF
 
C...Check minimum invariant mass of remnant system(s).
      psys(0,4)=p(i1,4)+p(i2,4)+0.5d0*vint(1)*(vint(151)+vint(152))
      psys(0,3)=p(i1,3)+p(i2,3)+0.5d0*vint(1)*(vint(151)-vint(152))
      pms(0)=max(0d0,psys(0,4)**2-psys(0,3)**2)
      pmin(0)=sqrt(pms(0))
      DO 340 jt=1,2
        psys(jt,4)=0.5d0*vint(1)*vint(142+jt)
        psys(jt,3)=psys(jt,4)*(-1)**(jt-1)
        pmin(jt)=0d0
        IF(mint(44+jt).EQ.1) goto 340
        mint(105)=mint(102+jt)
        mint(109)=mint(106+jt)
        CALL pyspli(mint(10+jt),mint(12+jt),kflch(jt),kflsp(jt))
        IF(mint(51).NE.0) THEN
          mint(57)=mint(57)+1
          RETURN
        ENDIF
        IF(kflch(jt).NE.0) pmin(jt)=pmin(jt)+pymass(kflch(jt))
        IF(kflsp(jt).NE.0) pmin(jt)=pmin(jt)+pymass(kflsp(jt))
        IF(kflch(jt)*kflsp(jt).NE.0) pmin(jt)=pmin(jt)+0.5d0*parp(111)
        pmin(jt)=sqrt(pmin(jt)**2+p(mint(83)+jt+2,1)**2+
     &  p(mint(83)+jt+2,2)**2)
  340 CONTINUE
      IF(pmin(0)+pmin(1)+pmin(2).GT.vint(1).OR.(mint(45).GE.2.AND.
     &pmin(1).GT.psys(1,4)).OR.(mint(46).GE.2.AND.pmin(2).GT.
     &psys(2,4))) THEN
        mint(51)=1
        mint(57)=mint(57)+1
        RETURN
      ENDIF
 
C...Loop over two remnants; skip if none there.
      i=ns
      DO 410 jt=1,2
        isn(jt)=0
        IF(mint(44+jt).EQ.1) goto 410
        IF(jt.EQ.1) ipu=ipu1
        IF(jt.EQ.2) ipu=ipu2
 
C...Store first remnant parton.
        i=i+1
        is(jt)=i
        isn(jt)=1
        DO 350 j=1,5
          k(i,j)=0
          p(i,j)=0d0
          v(i,j)=0d0
  350   CONTINUE
        k(i,1)=1
        k(i,2)=kflsp(jt)
        k(i,3)=mint(83)+jt
        p(i,5)=pymass(k(i,2))
 
C...First parton colour connections and kinematics.
        kcol=kchg(pycomp(kflsp(jt)),2)
        IF(kcol.EQ.2) THEN
          k(i,1)=3
          k(i,4)=mstu(5)*ipu+ipu
          k(i,5)=mstu(5)*ipu+ipu
          k(ipu,4)=mod(k(ipu,4),mstu(5))+mstu(5)*i
          k(ipu,5)=mod(k(ipu,5),mstu(5))+mstu(5)*i
        ELSEIF(kcol.NE.0) THEN
          k(i,1)=3
          kfls=(3-kcol*isign(1,kflsp(jt)))/2
          k(i,kfls+3)=ipu
          k(ipu,6-kfls)=mod(k(ipu,6-kfls),mstu(5))+mstu(5)*i
        ENDIF
        IF(kflch(jt).EQ.0) THEN
          p(i,1)=-p(mint(83)+jt+2,1)
          p(i,2)=-p(mint(83)+jt+2,2)
          pms(jt)=p(i,5)**2+p(i,1)**2+p(i,2)**2
          psys(jt,3)=sqrt(max(0d0,psys(jt,4)**2-pms(jt)))*(-1)**(jt-1)
          p(i,3)=psys(jt,3)
          p(i,4)=psys(jt,4)
 
C...When extra remnant parton or hadron: store extra remnant.
        ELSE
          i=i+1
          isn(jt)=2
          DO 360 j=1,5
            k(i,j)=0
            p(i,j)=0d0
            v(i,j)=0d0
  360     CONTINUE
          k(i,1)=1
          k(i,2)=kflch(jt)
          k(i,3)=mint(83)+jt
          p(i,5)=pymass(k(i,2))
 
C...Find parton colour connections of extra remnant.
          kcol=kchg(pycomp(kflch(jt)),2)
          IF(kcol.EQ.2) THEN
            k(i,1)=3
            k(i,4)=mstu(5)*ipu+ipu
            k(i,5)=mstu(5)*ipu+ipu
            k(ipu,4)=mod(k(ipu,4),mstu(5))+mstu(5)*i
            k(ipu,5)=mod(k(ipu,5),mstu(5))+mstu(5)*i
          ELSEIF(kcol.NE.0) THEN
            k(i,1)=3
            kfls=(3-kcol*isign(1,kflch(jt)))/2
            k(i,kfls+3)=ipu
            k(ipu,6-kfls)=mod(k(ipu,6-kfls),mstu(5))+mstu(5)*i
          ENDIF
 
C...Relative transverse momentum when two remnants.
          loop=0
  370     loop=loop+1
          CALL pyptdi(1,p(i-1,1),p(i-1,2))
          IF(iabs(mint(10+jt)).LT.20) THEN
            p(i-1,1)=0d0
            p(i-1,2)=0d0
          ELSE
            p(i-1,1)=p(i-1,1)-0.5d0*p(mint(83)+jt+2,1)
            p(i-1,2)=p(i-1,2)-0.5d0*p(mint(83)+jt+2,2)
          ENDIF
          pms(jt+2)=p(i-1,5)**2+p(i-1,1)**2+p(i-1,2)**2
          p(i,1)=-p(mint(83)+jt+2,1)-p(i-1,1)
          p(i,2)=-p(mint(83)+jt+2,2)-p(i-1,2)
          pms(jt+4)=p(i,5)**2+p(i,1)**2+p(i,2)**2
 
C...Meson or baryon; photon as meson. For splitup below.
          imb=1
          IF(mod(mint(10+jt)/1000,10).NE.0) imb=2
 
C***Relative distribution for electron into two electrons. Temporary!
          IF(iabs(mint(10+jt)).LT.20.AND.mint(14+jt).EQ.-mint(10+jt))
     &    THEN
            chi(jt)=pyr(0)
 
C...Relative distribution of electron energy into electron plus parton.
          ELSEIF(iabs(mint(10+jt)).LT.20) THEN
            xhrd=vint(140+jt)
            xe=vint(154+jt)
            chi(jt)=(xe-xhrd)/(1d0-xhrd)
 
C...Relative distribution of energy for particle into two jets.
          ELSEIF(iabs(kflch(jt)).LE.10.OR.kflch(jt).EQ.21) THEN
            chik=parp(92+2*imb)
            IF(mstp(92).LE.1) THEN
              IF(imb.EQ.1) chi(jt)=pyr(0)
              IF(imb.EQ.2) chi(jt)=1d0-sqrt(pyr(0))
            ELSEIF(mstp(92).EQ.2) THEN
              chi(jt)=1d0-pyr(0)**(1d0/(1d0+chik))
            ELSEIF(mstp(92).EQ.3) THEN
              cut=2d0*0.3d0/vint(1)
  380         chi(jt)=pyr(0)**2
              IF((chi(jt)**2/(chi(jt)**2+cut**2))**0.25d0*
     &        (1d0-chi(jt))**chik.LT.pyr(0)) goto 380
            ELSEIF(mstp(92).EQ.4) THEN
              cut=2d0*0.3d0/vint(1)
              cutr=(1d0+sqrt(1d0+cut**2))/cut
  390         chir=cut*cutr**pyr(0)
              chi(jt)=(chir**2-cut**2)/(2d0*chir)
              IF((1d0-chi(jt))**chik.LT.pyr(0)) goto 390
            ELSE
              cut=2d0*0.3d0/vint(1)
              cuta=cut**(1d0-parp(98))
              cutb=(1d0+cut)**(1d0-parp(98))
  400         chi(jt)=(cuta+pyr(0)*(cutb-cuta))**(1d0/(1d0-parp(98)))
              IF(((chi(jt)+cut)**2/(2d0*(chi(jt)**2+cut**2)))**
     &        (0.5d0*parp(98))*(1d0-chi(jt))**chik.LT.pyr(0)) goto 400
            ENDIF
 
C...Relative distribution of energy for particle into jet plus particle.
          ELSE
            IF(mstp(94).LE.1) THEN
              IF(imb.EQ.1) chi(jt)=pyr(0)
              IF(imb.EQ.2) chi(jt)=1d0-sqrt(pyr(0))
              IF(mod(kflch(jt)/1000,10).NE.0) chi(jt)=1d0-chi(jt)
            ELSEIF(mstp(94).EQ.2) THEN
              chi(jt)=1d0-pyr(0)**(1d0/(1d0+parp(93+2*imb)))
              IF(mod(kflch(jt)/1000,10).NE.0) chi(jt)=1d0-chi(jt)
            ELSEIF(mstp(94).EQ.3) THEN
              CALL pyzdis(1,0,pms(jt+4),zz)
              chi(jt)=zz
            ELSE
              CALL pyzdis(1000,0,pms(jt+4),zz)
              chi(jt)=zz
            ENDIF
          ENDIF
 
C...Construct total transverse mass; reject if too large.
          chi(jt)=max(1d-8,min(1d0-1d-8,chi(jt)))
          pms(jt)=pms(jt+4)/chi(jt)+pms(jt+2)/(1d0-chi(jt))
          IF(pms(jt).GT.psys(jt,4)**2) THEN
            IF(loop.LT.100) THEN
              goto 370
            ELSE
              mint(51)=1
              mint(57)=mint(57)+1
              RETURN
            ENDIF
          ENDIF
          psys(jt,3)=sqrt(max(0d0,psys(jt,4)**2-pms(jt)))*(-1)**(jt-1)
          vint(158+jt)=chi(jt)
 
C...Subdivide longitudinal momentum according to value selected above.
          pw1=chi(jt)*(psys(jt,4)+abs(psys(jt,3)))
          p(is(jt)+1,4)=0.5d0*(pw1+pms(jt+4)/pw1)
          p(is(jt)+1,3)=0.5d0*(pw1-pms(jt+4)/pw1)*(-1)**(jt-1)
          p(is(jt),4)=psys(jt,4)-p(is(jt)+1,4)
          p(is(jt),3)=psys(jt,3)-p(is(jt)+1,3)
        ENDIF
  410 CONTINUE
      n=i
 
C...Check if longitudinal boosts needed - if so pick two systems.
      pdev=abs(psys(0,4)+psys(1,4)+psys(2,4)-vint(1))+
     &abs(psys(0,3)+psys(1,3)+psys(2,3))
      IF(pdev.LE.1d-6*vint(1)) RETURN
      IF(isn(1).EQ.0) THEN
        ir=0
        il=2
      ELSEIF(isn(2).EQ.0) THEN
        ir=1
        il=0
      ELSEIF(vint(143).GT.0.2d0.AND.vint(144).GT.0.2d0) THEN
        ir=1
        il=2
      ELSEIF(vint(143).GT.0.2d0) THEN
        ir=1
        il=0
      ELSEIF(vint(144).GT.0.2d0) THEN
        ir=0
        il=2
      ELSEIF(pms(1)/psys(1,4)**2.GT.pms(2)/psys(2,4)**2) THEN
        ir=1
        il=0
      ELSE
        ir=0
        il=2
      ENDIF
      ig=3-ir-il
 
C...E+-pL wanted for system to be modified.
      IF((ig.EQ.1.AND.isn(1).EQ.0).OR.(ig.EQ.2.AND.isn(2).EQ.0)) THEN
        ppb=vint(1)
        pnb=vint(1)
      ELSE
        ppb=vint(1)-(psys(ig,4)+psys(ig,3))
        pnb=vint(1)-(psys(ig,4)-psys(ig,3))
      ENDIF
 
C...To keep x and Q2 in leptoproduction: do not count scattered lepton.
      IF(idisxq.EQ.1.AND.ig.NE.0) THEN
        ppb=ppb-(psys(0,4)+psys(0,3))
        pnb=pnb-(psys(0,4)-psys(0,3))
        DO 420 j=1,4
          psys(0,j)=0d0
  420   CONTINUE
        DO 450 i=mint(84)+1,ns
          IF(k(i,1).GT.10) goto 450
          incl=0
          iorig=i
  430     IF(iorig.EQ.lqout.OR.iorig.EQ.lpin+2) incl=1
          iorig=k(iorig,3)
          IF(iorig.GT.lpin) goto 430
          IF(incl.EQ.0) goto 450
          DO 440 j=1,4
            psys(0,j)=psys(0,j)+p(i,j)
  440     CONTINUE
  450   CONTINUE
        pms(0)=max(0d0,psys(0,4)**2-psys(0,3)**2)
        ppb=ppb+(psys(0,4)+psys(0,3))
        pnb=pnb+(psys(0,4)-psys(0,3))
      ENDIF
 
C...Construct longitudinal boosts.
      dpmtb=ppb*pnb
      dpmtr=pms(ir)
      dpmtl=pms(il)
      dsqlam=sqrt(max(0d0,(dpmtb-dpmtr-dpmtl)**2-4d0*dpmtr*dpmtl))
      IF(dsqlam.LE.1d-6*dpmtb) THEN
        mint(51)=1
        mint(57)=mint(57)+1
        RETURN
      ENDIF
      dsqsgn=sign(1d0,psys(ir,3)*psys(il,4)-psys(il,3)*psys(ir,4))
      drkr=(dpmtb+dpmtr-dpmtl+dsqlam*dsqsgn)/
     &(2d0*(psys(ir,4)+psys(ir,3))*pnb)
      drkl=(dpmtb+dpmtl-dpmtr+dsqlam*dsqsgn)/
     &(2d0*(psys(il,4)-psys(il,3))*ppb)
      dber=(drkr**2-1d0)/(drkr**2+1d0)
      dbel=-(drkl**2-1d0)/(drkl**2+1d0)
 
C...Perform longitudinal boosts.
      IF(ir.EQ.1.AND.isn(1).EQ.1.AND.dber.LE.-0.99999999d0) THEN
        p(is(1),3)=0d0
        p(is(1),4)=sqrt(p(is(1),5)**2+p(is(1),1)**2+p(is(1),2)**2)
      ELSEIF(ir.EQ.1) THEN
        CALL pyrobo(is(1),is(1)+isn(1)-1,0d0,0d0,0d0,0d0,dber)
      ELSEIF(idisxq.EQ.1) THEN
        DO 470 i=i1,ns
          incl=0
          iorig=i
  460     IF(iorig.EQ.lqout.OR.iorig.EQ.lpin+2) incl=1
          iorig=k(iorig,3)
          IF(iorig.GT.lpin) goto 460
          IF(incl.EQ.1) CALL pyrobo(i,i,0d0,0d0,0d0,0d0,dber)
  470   CONTINUE
      ELSE
        CALL pyrobo(i1,ns,0d0,0d0,0d0,0d0,dber)
      ENDIF
      IF(il.EQ.2.AND.isn(2).EQ.1.AND.dbel.GE.0.99999999d0) THEN
        p(is(2),3)=0d0
        p(is(2),4)=sqrt(p(is(2),5)**2+p(is(2),1)**2+p(is(2),2)**2)
      ELSEIF(il.EQ.2) THEN
        CALL pyrobo(is(2),is(2)+isn(2)-1,0d0,0d0,0d0,0d0,dbel)
      ELSEIF(idisxq.EQ.1) THEN
        DO 490 i=i1,ns
          incl=0
          iorig=i
  480     IF(iorig.EQ.lqout.OR.iorig.EQ.lpin+2) incl=1
          iorig=k(iorig,3)
          IF(iorig.GT.lpin) goto 480
          IF(incl.EQ.1) CALL pyrobo(i,i,0d0,0d0,0d0,0d0,dbel)
  490   CONTINUE
      ELSE
        CALL pyrobo(i1,ns,0d0,0d0,0d0,0d0,dbel)
      ENDIF
 
C...Final check that energy-momentum conservation worked.
      pesum=0d0
      pzsum=0d0
      DO 500 i=mint(84)+1,n
        IF(k(i,1).GT.10) goto 500
        pesum=pesum+p(i,4)
        pzsum=pzsum+p(i,3)
  500 CONTINUE
      pdev=abs(pesum-vint(1))+abs(pzsum)
      IF(pdev.GT.1d-4*vint(1)) THEN
        mint(51)=1
        mint(57)=mint(57)+1
        RETURN
      ENDIF
 
C...Calculate rotation and boost from overall CM frame to
C...hadronic CM frame in leptoproduction.
      mint(91)=0
      IF(mint(82).EQ.1.AND.(mint(43).EQ.2.OR.mint(43).EQ.3)) THEN
        mint(91)=1
        lesd=1
        IF(mint(42).EQ.1) lesd=2
        lpin=mint(83)+3-lesd
 
C...Sum upp momenta of everything not lepton or photon to define boost.
        DO 510 j=1,4
          psum(j)=0d0
  510   CONTINUE
        DO 530 i=1,n
          IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 530
          IF(iabs(k(i,2)).GE.11.AND.iabs(k(i,2)).LE.20) goto 530
          IF(k(i,2).EQ.22) goto 530
          DO 520 j=1,4
            psum(j)=psum(j)+p(i,j)
  520     CONTINUE
  530   CONTINUE
        vint(223)=-psum(1)/psum(4)
        vint(224)=-psum(2)/psum(4)
        vint(225)=-psum(3)/psum(4)
 
C...Boost incoming hadron to hadronic CM frame to determine rotations.
        k(n+1,1)=1
        DO 540 j=1,5
          p(n+1,j)=p(lpin,j)
          v(n+1,j)=v(lpin,j)
  540   CONTINUE
        CALL pyrobo(n+1,n+1,0d0,0d0,vint(223),vint(224),vint(225))
        vint(222)=-pyangl(p(n+1,1),p(n+1,2))
        CALL pyrobo(n+1,n+1,0d0,vint(222),0d0,0d0,0d0)
        IF(lesd.EQ.2) THEN
          vint(221)=-pyangl(p(n+1,3),p(n+1,1))
        ELSE
          vint(221)=pyangl(-p(n+1,3),p(n+1,1))
        ENDIF
      ENDIF
 
      RETURN
      END