pyreco.f

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

 
C***********************************************************************
 
C...PYRECO
C...Handles the possibility of colour reconnection in W+W- events,
C...Based on the main scenarios of the Sjostrand and Khoze study:
C...I, II, II', intermediate and instantaneous; plus one model
C...along the lines of the Gustafson and Hakkinen: GH.
C...Note: also handles Z0 Z0 and W-W+ events, but notation below
C...is as if first resonance is W+ and second W-.
 
      SUBROUTINE pyreco(IW1,IW2,NSD1,NAFT1)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter value; number of points in MC integration.
      parameter(npt=100)
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 nbeg(2),nend(2),inp(50),inm(50),beww(3),xp(3),xm(3),
     &v1(3),v2(3),betp(50,4),dirp(50,3),betm(50,4),dirm(50,3),
     &xd(4),xb(4),iap(npt),iam(npt),wta(npt),v1p(3),v2p(3),v1m(3),
     &v2m(3),q(4,3),xpp(3),xmm(3),ipc(20),imc(20),tc(0:20),tpc(20),
     &tmc(20),ijoin(100)
 
C...Functions to give four-product and to do determinants.
      four(i,j)=p(i,4)*p(j,4)-p(i,1)*p(j,1)-p(i,2)*p(j,2)-p(i,3)*p(j,3)
      deter(i,j,l)=q(i,1)*q(j,2)*q(l,3)-q(i,1)*q(l,2)*q(j,3)+
     &q(j,1)*q(l,2)*q(i,3)-q(j,1)*q(i,2)*q(l,3)+
     &q(l,1)*q(i,2)*q(j,3)-q(l,1)*q(j,2)*q(i,3)
 
C...Only allow fraction of recoupling for GH, intermediate and
C...instantaneous.
      IF(mstp(115).EQ.5.OR.mstp(115).EQ.11.OR.mstp(115).EQ.12) THEN
        IF(pyr(0).GT.parp(120)) RETURN
      ENDIF
      isub=mint(1)
 
C...Common part for scenarios I, II, II', and GH.
      IF(mstp(115).EQ.1.OR.mstp(115).EQ.2.OR.mstp(115).EQ.3.OR.
     &mstp(115).EQ.5) THEN
 
C...Read out frequently-used parameters.
        pi=paru(1)
        hbar=paru(3)
        pmw=pmas(24,1)
        IF(isub.EQ.22) pmw=pmas(23,1)
        pgw=pmas(24,2)
        IF(isub.EQ.22) pgw=pmas(23,2)
        tfrag=parp(115)
        rhad=parp(116)
        fact=parp(117)
        blowr=parp(118)
        blowt=parp(119)
 
C...Find range of decay products of the W's.
C...Background: the W's are stored in IW1 and IW2.
C...Their direct decay products in NSD1+1 through NSD1+4.
C...Products after shower (if any) in NSD1+5 through NAFT1
C...for first W and in NAFT1+1 through N for the second.
        IF(naft1.GT.nsd1+4) THEN
          nbeg(1)=nsd1+5
          nend(1)=naft1
        ELSE
          nbeg(1)=nsd1+1
          nend(1)=nsd1+2
        ENDIF
        IF(n.GT.naft1) THEN
          nbeg(2)=naft1+1
          nend(2)=n
        ELSE
          nbeg(2)=nsd1+3
          nend(2)=nsd1+4
        ENDIF
 
C...Rearrange parton shower products along strings.
        nold=n
        CALL pyprep(nsd1+1)
        IF(mint(51).NE.0) RETURN
 
C...Find partons pointing back to W+ and W-; store them with quark
C...end of string first.
        nnp=0
        nnm=0
        isgp=0
        isgm=0
        DO 120 i=nold+1,n
          IF(k(i,1).NE.1.AND.k(i,1).NE.2) goto 120
          IF(iabs(k(i,2)).GE.22) goto 120
          IF(k(i,3).GE.nbeg(1).AND.k(i,3).LE.nend(1)) THEN
            IF(isgp.EQ.0) isgp=isign(1,k(i,2))
            nnp=nnp+1
            IF(isgp.EQ.1) THEN
              inp(nnp)=i
            ELSE
              DO 100 i1=nnp,2,-1
                inp(i1)=inp(i1-1)
  100         CONTINUE
              inp(1)=i
            ENDIF
            IF(k(i,1).EQ.1) isgp=0
          ELSEIF(k(i,3).GE.nbeg(2).AND.k(i,3).LE.nend(2)) THEN
            IF(isgm.EQ.0) isgm=isign(1,k(i,2))
            nnm=nnm+1
            IF(isgm.EQ.1) THEN
              inm(nnm)=i
            ELSE
              DO 110 i1=nnm,2,-1
                inm(i1)=inm(i1-1)
  110         CONTINUE
              inm(1)=i
            ENDIF
            IF(k(i,1).EQ.1) isgm=0
          ENDIF
  120   CONTINUE
 
C...Boost to W+W- rest frame (not strictly needed).
        DO 130 j=1,3
          beww(j)=(p(iw1,j)+p(iw2,j))/(p(iw1,4)+p(iw2,4))
  130   CONTINUE
        CALL pyrobo(iw1,iw1,0d0,0d0,-beww(1),-beww(2),-beww(3))
        CALL pyrobo(iw2,iw2,0d0,0d0,-beww(1),-beww(2),-beww(3))
        CALL pyrobo(nold+1,n,0d0,0d0,-beww(1),-beww(2),-beww(3))
 
C...Select decay vertices of W+ and W-.
        tp=hbar*(-log(pyr(0)))*p(iw1,4)/
     &  sqrt((p(iw1,5)**2-pmw**2)**2+(p(iw1,5)**2*pgw/pmw)**2)
        tm=hbar*(-log(pyr(0)))*p(iw2,4)/
     &  sqrt((p(iw2,5)**2-pmw**2)**2+(p(iw2,5)**2*pgw/pmw)**2)
        gtmax=max(tp,tm)
        DO 140 j=1,3
          xp(j)=tp*p(iw1,j)/p(iw1,4)
          xm(j)=tm*p(iw2,j)/p(iw2,4)
  140   CONTINUE
 
C...Begin scenario I specifics.
        IF(mstp(115).EQ.1) THEN
 
C...Reconstruct velocity and direction of W+ string pieces.
          DO 170 iip=1,nnp-1
            IF(k(inp(iip),2).LT.0) goto 170
            i1=inp(iip)
            i2=inp(iip+1)
            p1a=sqrt(p(i1,1)**2+p(i1,2)**2+p(i1,3)**2)
            p2a=sqrt(p(i2,1)**2+p(i2,2)**2+p(i2,3)**2)
            DO 150 j=1,3
              v1(j)=p(i1,j)/p1a
              v2(j)=p(i2,j)/p2a
              betp(iip,j)=0.5d0*(v1(j)+v2(j))
              dirp(iip,j)=v1(j)-v2(j)
  150       CONTINUE
            betp(iip,4)=1d0/sqrt(1d0-betp(iip,1)**2-betp(iip,2)**2-
     &      betp(iip,3)**2)
            dirl=sqrt(dirp(iip,1)**2+dirp(iip,2)**2+dirp(iip,3)**2)
            DO 160 j=1,3
              dirp(iip,j)=dirp(iip,j)/dirl
  160       CONTINUE
  170     CONTINUE
 
C...Reconstruct velocity and direction of W- string pieces.
          DO 200 iim=1,nnm-1
            IF(k(inm(iim),2).LT.0) goto 200
            i1=inm(iim)
            i2=inm(iim+1)
            p1a=sqrt(p(i1,1)**2+p(i1,2)**2+p(i1,3)**2)
            p2a=sqrt(p(i2,1)**2+p(i2,2)**2+p(i2,3)**2)
            DO 180 j=1,3
              v1(j)=p(i1,j)/p1a
              v2(j)=p(i2,j)/p2a
              betm(iim,j)=0.5d0*(v1(j)+v2(j))
              dirm(iim,j)=v1(j)-v2(j)
  180       CONTINUE
            betm(iim,4)=1d0/sqrt(1d0-betm(iim,1)**2-betm(iim,2)**2-
     &      betm(iim,3)**2)
            dirl=sqrt(dirm(iim,1)**2+dirm(iim,2)**2+dirm(iim,3)**2)
            DO 190 j=1,3
              dirm(iim,j)=dirm(iim,j)/dirl
  190       CONTINUE
  200     CONTINUE
 
C...Loop over number of space-time points.
          nacc=0
          sum=0d0
          DO 250 ipt=1,npt
 
C...Pick x,y,z,t Gaussian (width RHAD and TFRAG, respectively).
            r=sqrt(-log(pyr(0)))
            phi=2d0*pi*pyr(0)
            x=blowr*rhad*r*cos(phi)
            y=blowr*rhad*r*sin(phi)
            r=sqrt(-log(pyr(0)))
            phi=2d0*pi*pyr(0)
            z=blowr*rhad*r*cos(phi)
            t=gtmax+blowt*sqrt(0.5d0)*tfrag*r*abs(sin(phi))
 
C...Reject impossible points. Weight for sample distribution.
            IF(t**2-x**2-y**2-z**2.LT.0d0) goto 250
            wtsmp=exp(-(x**2+y**2+z**2)/(blowr*rhad)**2)*
     &      exp(-2d0*(t-gtmax)**2/(blowt*tfrag)**2)
 
C...Loop over W+ string pieces and find one with largest weight.
            imaxp=0
            wtmaxp=1d-10
            xd(1)=x-xp(1)
            xd(2)=y-xp(2)
            xd(3)=z-xp(3)
            xd(4)=t-tp
            DO 220 iip=1,nnp-1
              IF(k(inp(iip),2).LT.0) goto 220
              bed=betp(iip,1)*xd(1)+betp(iip,2)*xd(2)+betp(iip,3)*xd(3)
              bedg=betp(iip,4)*(betp(iip,4)*bed/(1d0+betp(iip,4))-xd(4))
              DO 210 j=1,3
                xb(j)=xd(j)+bedg*betp(iip,j)
  210         CONTINUE
              xb(4)=betp(iip,4)*(xd(4)-bed)
              sr2=xb(1)**2+xb(2)**2+xb(3)**2
              sz2=(dirp(iip,1)*xb(1)+dirp(iip,2)*xb(2)+
     &        dirp(iip,3)*xb(3))**2
              wtp=exp(-(sr2-sz2)/(2d0*rhad**2))*exp(-(xb(4)**2-sz2)/
     &        tfrag**2)
              IF(xb(4)-sqrt(sr2).LT.0d0) wtp=0d0
              IF(wtp.GT.wtmaxp) THEN
                imaxp=iip
                wtmaxp=wtp
              ENDIF
  220       CONTINUE
 
C...Loop over W- string pieces and find one with largest weight.
            imaxm=0
            wtmaxm=1d-10
            xd(1)=x-xm(1)
            xd(2)=y-xm(2)
            xd(3)=z-xm(3)
            xd(4)=t-tm
            DO 240 iim=1,nnm-1
              IF(k(inm(iim),2).LT.0) goto 240
              bed=betm(iim,1)*xd(1)+betm(iim,2)*xd(2)+betm(iim,3)*xd(3)
              bedg=betm(iim,4)*(betm(iim,4)*bed/(1d0+betm(iim,4))-xd(4))
              DO 230 j=1,3
                xb(j)=xd(j)+bedg*betm(iim,j)
  230         CONTINUE
              xb(4)=betm(iim,4)*(xd(4)-bed)
              sr2=xb(1)**2+xb(2)**2+xb(3)**2
              sz2=(dirm(iim,1)*xb(1)+dirm(iim,2)*xb(2)+
     &        dirm(iim,3)*xb(3))**2
              wtm=exp(-(sr2-sz2)/(2d0*rhad**2))*exp(-(xb(4)**2-sz2)/
     &        tfrag**2)
              IF(xb(4)-sqrt(sr2).LT.0d0) wtm=0d0
              IF(wtm.GT.wtmaxm) THEN
                imaxm=iim
                wtmaxm=wtm
              ENDIF
  240       CONTINUE
 
C...Result of integration.
            wt=0d0
            IF(imaxp.NE.0.AND.imaxm.NE.0) THEN
              wt=wtmaxp*wtmaxm/wtsmp
              sum=sum+wt
              nacc=nacc+1
              iap(nacc)=imaxp
              iam(nacc)=imaxm
              wta(nacc)=wt
            ENDIF
  250     CONTINUE
          res=blowr**3*blowt*sum/npt
 
C...Decide whether to reconnect and, if so, where.
          iacc=0
          prec=1d0-exp(-fact*res)
          IF(prec.GT.pyr(0)) THEN
            rsum=pyr(0)*sum
            DO 260 ia=1,nacc
              iacc=ia
              rsum=rsum-wta(ia)
              IF(rsum.LE.0d0) goto 270
  260       CONTINUE
  270       iip=iap(iacc)
            iim=iam(iacc)
          ENDIF
 
C...Begin scenario II and II' specifics.
        ELSEIF(mstp(115).EQ.2.OR.mstp(115).EQ.3) THEN
 
C...Loop through all string pieces, one from W+ and one from W-.
          ncross=0
          tc(0)=0d0
          DO 340 iip=1,nnp-1
            IF(k(inp(iip),2).LT.0) goto 340
            i1p=inp(iip)
            i2p=inp(iip+1)
            DO 330 iim=1,nnm-1
              IF(k(inm(iim),2).LT.0) goto 330
              i1m=inm(iim)
              i2m=inm(iim+1)
 
C...Find endpoint velocity vectors.
              DO 280 j=1,3
                v1p(j)=p(i1p,j)/p(i1p,4)
                v2p(j)=p(i2p,j)/p(i2p,4)
                v1m(j)=p(i1m,j)/p(i1m,4)
                v2m(j)=p(i2m,j)/p(i2m,4)
  280         CONTINUE
 
C...Define q matrix and find t.
              DO 290 j=1,3
                q(1,j)=v2p(j)-v1p(j)
                q(2,j)=-(v2m(j)-v1m(j))
                q(3,j)=xp(j)-xm(j)-tp*v1p(j)+tm*v1m(j)
                q(4,j)=v1p(j)-v1m(j)
  290         CONTINUE
              t=-deter(1,2,3)/deter(1,2,4)
 
C...Find alpha and beta; i.e. coordinates of crossing point.
              s11=q(1,1)*(t-tp)
              s12=q(2,1)*(t-tm)
              s13=q(3,1)+q(4,1)*t
              s21=q(1,2)*(t-tp)
              s22=q(2,2)*(t-tm)
              s23=q(3,2)+q(4,2)*t
              den=s11*s22-s12*s21
              alp=(s12*s23-s22*s13)/den
              bet=(s21*s13-s11*s23)/den
 
C...Check if solution acceptable.
              iansw=1
              IF(t.LT.gtmax) iansw=0
              IF(alp.LT.0d0.OR.alp.GT.1d0) iansw=0
              IF(bet.LT.0d0.OR.bet.GT.1d0) iansw=0
 
C...Find point of crossing and check that not inconsistent.
              DO 300 j=1,3
                xpp(j)=xp(j)+(v1p(j)+alp*(v2p(j)-v1p(j)))*(t-tp)
                xmm(j)=xm(j)+(v1m(j)+bet*(v2m(j)-v1m(j)))*(t-tm)
  300         CONTINUE
              d2pm=(xpp(1)-xmm(1))**2+(xpp(2)-xmm(2))**2+
     &        (xpp(3)-xmm(3))**2
              d2p=xpp(1)**2+xpp(2)**2+xpp(3)**2
              d2m=xmm(1)**2+xmm(2)**2+xmm(3)**2
              IF(d2pm.GT.1d-4*(d2p+d2m)) iansw=-1
 
C...Find string eigentimes at crossing.
              IF(iansw.EQ.1) THEN
                taup=sqrt(max(0d0,(t-tp)**2-(xpp(1)-xp(1))**2-
     &          (xpp(2)-xp(2))**2-(xpp(3)-xp(3))**2))
                taum=sqrt(max(0d0,(t-tm)**2-(xmm(1)-xm(1))**2-
     &          (xmm(2)-xm(2))**2-(xmm(3)-xm(3))**2))
              ELSE
                taup=0d0
                taum=0d0
              ENDIF
 
C...Order crossings by time. End loop over crossings.
              IF(iansw.EQ.1.AND.ncross.LT.20) THEN
                ncross=ncross+1
                DO 310 i1=ncross,1,-1
                  IF(t.GT.tc(i1-1).OR.i1.EQ.1) THEN
                    ipc(i1)=iip
                    imc(i1)=iim
                    tc(i1)=t
                    tpc(i1)=taup
                    tmc(i1)=taum
                    goto 320
                  ELSE
                    ipc(i1)=ipc(i1-1)
                    imc(i1)=imc(i1-1)
                    tc(i1)=tc(i1-1)
                    tpc(i1)=tpc(i1-1)
                    tmc(i1)=tmc(i1-1)
                  ENDIF
  310           CONTINUE
  320           CONTINUE
              ENDIF
  330       CONTINUE
  340     CONTINUE
 
C...Loop over crossings; find first (if any) acceptable one.
          iacc=0
          IF(ncross.GE.1) THEN
            DO 350 ic=1,ncross
              pnfrag=exp(-(tpc(ic)**2+tmc(ic)**2)/tfrag**2)
              IF(pnfrag.GT.pyr(0)) THEN
C...Scenario II: only compare with fragmentation time.
                IF(mstp(115).EQ.2) THEN
                  iacc=ic
                  iip=ipc(iacc)
                  iim=imc(iacc)
                  goto 360
C...Scenario II': also require that string length decreases.
                ELSE
                  iip=ipc(ic)
                  iim=imc(ic)
                  i1p=inp(iip)
                  i2p=inp(iip+1)
                  i1m=inm(iim)
                  i2m=inm(iim+1)
                  elold=four(i1p,i2p)*four(i1m,i2m)
                  elnew=four(i1p,i2m)*four(i1m,i2p)
                  IF(elnew.LT.elold) THEN
                    iacc=ic
                    iip=ipc(iacc)
                    iim=imc(iacc)
                    goto 360
                  ENDIF
                ENDIF
              ENDIF
  350       CONTINUE
  360       CONTINUE
          ENDIF
 
C...Begin scenario GH specifics.
        ELSEIF(mstp(115).EQ.5) THEN
 
C...Loop through all string pieces, one from W+ and one from W-.
          iacc=0
          elmin=1d0
          DO 380 iip=1,nnp-1
            IF(k(inp(iip),2).LT.0) goto 380
            i1p=inp(iip)
            i2p=inp(iip+1)
            DO 370 iim=1,nnm-1
              IF(k(inm(iim),2).LT.0) goto 370
              i1m=inm(iim)
              i2m=inm(iim+1)
 
C...Look for largest decrease of (exponent of) Lambda measure.
              elold=four(i1p,i2p)*four(i1m,i2m)
              elnew=four(i1p,i2m)*four(i1m,i2p)
              eldif=elnew/max(1d-10,elold)
              IF(eldif.LT.elmin) THEN
                iacc=iip+iim
                elmin=eldif
                ipc(1)=iip
                imc(1)=iim
              ENDIF
  370       CONTINUE
  380     CONTINUE
          iip=ipc(1)
          iim=imc(1)
        ENDIF
 
C...Common for scenarios I, II, II' and GH: reconnect strings.
        IF(iacc.NE.0) THEN
          mint(32)=1
          njoin=0
          DO 390 is=1,nnp+nnm
            njoin=njoin+1
            IF(is.LE.iip) THEN
              i=inp(is)
            ELSEIF(is.LE.iip+nnm-iim) THEN
              i=inm(is-iip+iim)
            ELSEIF(is.LE.iip+nnm) THEN
              i=inm(is-iip-nnm+iim)
            ELSE
              i=inp(is-nnm)
            ENDIF
            ijoin(njoin)=i
            IF(k(i,2).LT.0) THEN
              CALL pyjoin(njoin,ijoin)
              njoin=0
            ENDIF
  390     CONTINUE
 
C...Restore original event record if no reconnection.
        ELSE
          DO 400 i=nsd1+1,nold
            IF(k(i,1).EQ.13.OR.k(i,1).EQ.14) THEN
              k(i,4)=mod(k(i,4),mstu(5)**2)
              k(i,5)=mod(k(i,5),mstu(5)**2)
            ENDIF
  400     CONTINUE
          DO 410 i=nold+1,n
            k(k(i,3),1)=3
  410     CONTINUE
          n=nold
        ENDIF
 
C...Boost back system.
        CALL pyrobo(iw1,iw1,0d0,0d0,beww(1),beww(2),beww(3))
        CALL pyrobo(iw2,iw2,0d0,0d0,beww(1),beww(2),beww(3))
        IF(n.GT.nold) CALL pyrobo(nold+1,n,0d0,0d0,
     &  beww(1),beww(2),beww(3))
 
C...Common part for intermediate and instantaneous scenarios.
      ELSEIF(mstp(115).EQ.11.OR.mstp(115).EQ.12) THEN
        mint(32)=1
 
C...Remove old shower products and reset showering ones.
        n=nsd1+4
        DO 420 i=nsd1+1,nsd1+4
          k(i,1)=3
          k(i,4)=mod(k(i,4),mstu(5)**2)
          k(i,5)=mod(k(i,5),mstu(5)**2)
  420   CONTINUE
 
C...Identify quark-antiquark pairs.
        iq1=nsd1+1
        iq2=nsd1+2
        iq3=nsd1+3
        IF(k(iq1,2)*k(iq3,2).LT.0) iq3=nsd1+4
        iq4=2*nsd1+7-iq3
 
C...Reconnect strings.
        ijoin(1)=iq1
        ijoin(2)=iq4
        CALL pyjoin(2,ijoin)
        ijoin(1)=iq3
        ijoin(2)=iq2
        CALL pyjoin(2,ijoin)
 
C...Do new parton showers in intermediate scenario.
        IF(mstp(71).GE.1.AND.mstp(115).EQ.11) THEN
          mstj50=mstj(50)
          mstj(50)=0
          CALL pyshow(iq1,iq2,p(iw1,5))
          CALL pyshow(iq3,iq4,p(iw2,5))
          mstj(50)=mstj50
 
C...Do new parton showers in instantaneous scenario.
        ELSEIF(mstp(71).GE.1.AND.mstp(115).EQ.12) THEN
          ppm2=(p(iq1,4)+p(iq4,4))**2-(p(iq1,1)+p(iq4,1))**2-
     &    (p(iq1,2)+p(iq4,2))**2-(p(iq1,3)+p(iq4,3))**2
          ppm=sqrt(max(0d0,ppm2))
          CALL pyshow(iq1,iq4,ppm)
          ppm2=(p(iq3,4)+p(iq2,4))**2-(p(iq3,1)+p(iq2,1))**2-
     &    (p(iq3,2)+p(iq2,2))**2-(p(iq3,3)+p(iq2,3))**2
          ppm=sqrt(max(0d0,ppm2))
          CALL pyshow(iq3,iq2,ppm)
        ENDIF
      ENDIF
 
      RETURN
      END