pyfscr.f

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

  
C*********************************************************************
 
C...PYFSCR
C...Performs colour annealing.
C...MSTP(95) : CR Type
C...         = 1  : old cut-and-paste reconnections, handled in PYMIHK
C...         = 2  : Type I(no gg loops); hadron-hadron only
C...         = 3  : Type I(no gg loops); all beams
C...         = 4  : Type II(gg loops)  ; hadron-hadron only
C...         = 5  : Type II(gg loops)  ; all beams
C...         = 6  : Type S             ; hadron-hadron only
C...         = 7  : Type S             ; all beams
C...Types I and II are described in Sandhoff+Skands, in hep-ph/0604120.
C...Type S is driven by starting only from free triplets, not octets.
C...A string piece remains unchanged with probability
C...    PKEEP = (1-PARP(78))**N
C...This scaling corresponds to each string piece having to go through
C...N other ones, each with probability PARP(78) for reconnection, where
C...N is here chosen simply as the number of multiple interactions,
C...for a rough scaling with the general level of activity.
 
      SUBROUTINE pyfscr(IP)
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      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/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pyint1/mint(400),vint(400)
C...The common block of colour tags.
      common/pyctag/nct,mct(4000,2)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pyint1/,/pyctag/,
     &/pypars/
C...MCN: Temporary storage of new colour tags
      DOUBLE PRECISION mcn(4000,2)
 
C...Function to give four-product.
      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)
 
C...Check valid range of MSTP(95), local copy
      IF (mstp(95).LE.1.OR.mstp(95).GE.8) RETURN
      mstp95=mod(mstp(95),10)
C...Set whether CR allowed inside resonance systems or not
C...(not implemented yet)
C      MRESCR=1
C      IF (MSTP(95).GE.10) MRESCR=0
 
C...Check whether colour tags already defined
      IF (mint(33).EQ.0) THEN
C...Erase any existing colour tags for this event
        DO 100 i=1,n
          mct(i,1)=0
          mct(i,2)=0
  100   CONTINUE
C...Create colour tags for this event
        DO 120 i=1,n
          IF (k(i,1).EQ.3) THEN
            DO 110 kcs=4,5
              kcsin=kcs
              IF (mct(i,kcsin-3).EQ.0) THEN
                CALL pycttr(i,kcsin,i)
              ENDIF
  110       CONTINUE
          ENDIF
  120 CONTINUE
C...Instruct PYPREP to use colour tags
        mint(33)=1
      ENDIF
 
C...For MSTP(95) even, only apply to hadron-hadron
      IF (mod(mstp(95),2).EQ.0) THEN
         ka1=iabs(mint(11))
         ka2=iabs(mint(12))
         IF (ka1.LT.100.OR.ka2.LT.100) goto 9999
      ENDIF
 
C...Initialize new tag array (but do not delete old yet)
      lct=nct
      DO 130 i=max(1,ip),n
         mcn(i,1)=0
         mcn(i,2)=0
  130 CONTINUE
 
C...For each final-state dipole, check whether string should be
C...preserved.
      DO 150 ict=1,nct
        ic=0
        ia=0
        DO 140 i=max(1,ip),n
          IF (k(i,1).EQ.3.AND.mct(i,1).EQ.ict) ic=i
          IF (k(i,1).EQ.3.AND.mct(i,2).EQ.ict) ia=i
  140   CONTINUE
        IF (ic.NE.0.AND.ia.NE.0) THEN
C...Chiefly consider large strings.
          pkeep=(1d0-parp(78))**mint(31)
          IF (pyr(0).LE.pkeep) THEN
            lct=lct+1
            mcn(ic,1)=lct
            mcn(ia,2)=lct
          ENDIF
        ENDIF
  150 CONTINUE
 
C...Loop over event record, starting from IP
C...(Ignore junctions for now.)
      nloop=0
  160 nloop=nloop+1
      mcimax=0
      mcjmax=0
      rlmax=0d0
      ilmax=0
      jlmax=0
      DO 230 i=max(1,ip),n
         IF (k(i,1).NE.3) goto 230
C...Check colour charge
         mci=kchg(pycomp(k(i,2)),2)*isign(1,k(i,2))
         IF (mci.EQ.0) goto 230
C...For Seattle algorithm, only start from partons with one dangling
C...colour tag
         IF (mstp(95).EQ.6.OR.mstp(95).EQ.7) THEN
           IF (mci.EQ.2.AND.mcn(i,1).EQ.0.AND.mcn(i,2).EQ.0) goto 230
         ENDIF
C...  Find optimal partner
         jlopt=0
         mcjopt=0
         mbropt=0
         mggopt=0
         rlopt=1d19
C...Loop over I colour/anticolour, check whether already connected
  170    DO 220 icl=1,2
            IF (mcn(i,icl).NE.0) goto 220
            IF (icl.EQ.1.AND.mci.EQ.-1) goto 220
            IF (icl.EQ.2.AND.mci.EQ.1) goto 220
C...Check whether this is a dangling colour tag (ie to junction!)
            ifound=0
            DO 180 j=max(1,ip),n
               IF (k(j,1).EQ.3.AND.mct(j,3-icl).EQ.mct(i,icl)) ifound=1
  180       CONTINUE
            IF (ifound.EQ.0) goto 220
            DO 210 j=max(1,ip),n
               IF (k(j,1).NE.3.OR.i.EQ.j) goto 210
C...Do not make direct connections between partons in same Beam Remnant
               mbrstr=0
               IF (k(i,3).LE.2.AND.k(j,3).LE.2.AND.k(i,3).EQ.k(j,3))
     &              mbrstr=1
C...Check colour charge
               mcj=kchg(pycomp(k(j,2)),2)*isign(1,k(j,2))
               IF (mcj.EQ.0.OR.(mcj.EQ.mci.AND.mci.NE.2)) goto 210
C...Check for gluon loops
               mggstr=0
               IF (mcj.EQ.2.AND.mci.EQ.2) THEN
                 icla=3-icl
                 IF (mcn(i,icla).EQ.mcn(j,icl).AND.mstp(95).LE.3.AND.
     &                mcn(i,icla).NE.0) mggstr=1
               ENDIF
C...Loop over J colour/anticolour, check whether already connected
               DO 200 jcl=1,2
                  IF (mcn(j,jcl).NE.0) goto 200
                  IF (jcl.EQ.icl) goto 200
                  IF (jcl.EQ.1.AND.mcj.EQ.-1) goto 200
                  IF (jcl.EQ.2.AND.mcj.EQ.1) goto 200
C...Check whether this is a dangling colour tag (ie to junction!)
                  ifound=0
                  DO 190 j2=max(1,ip),n
                     IF (k(j2,1).EQ.3.AND.mct(j2,3-jcl).EQ.mct(j,jcl))
     &                    ifound=1
  190             CONTINUE
                  IF (ifound.EQ.0) goto 200
C...Save connection with smallest lambda measure
C...If best so far was a BR string and this is not, also save.
C...If best so far was a gg string and this is not, also save.
                  rl=four(i,j)
                  IF (rl.LT.rlopt.OR.(rl.EQ.rlopt.AND.pyr(0).LE.0.5d0)
     &                 .OR.(mbropt.EQ.1.AND.mbrstr.EQ.0)
     &                 .OR.(mggopt.EQ.1.AND.mggstr.EQ.0)) THEN
                     rlopt=rl
                     jlopt=j
                     icopt=icl
                     jcopt=jcl
                     mcjopt=mcj
                     mbropt=mbrstr
                     mggopt=mggstr
                  ENDIF
  200          CONTINUE
  210       CONTINUE
  220    CONTINUE
         IF (jlopt.NE.0) THEN
C...Save pair with largest RLOPT so far
            IF (rlopt.GE.rlmax) THEN
               rlmax=rlopt
               ilmax=i
               jlmax=jlopt
               icmax=icopt
               jcmax=jcopt
               mcjmax=mcjopt
               mcimax=mci
            ENDIF
         ENDIF
  230 CONTINUE
C...Save and iterate
      IF (ilmax.GT.0) THEN
         lct=lct+1
         mcn(ilmax,icmax)=lct
         mcn(jlmax,jcmax)=lct
         IF (nloop.LE.2*(n-ip)) THEN
            goto 160
         ELSE
            CALL pyerrm(31,' PYFSCR: infinite loop in color annealing')
            CALL pystop(11)
         ENDIF
      ELSE
C...Save and exit. First check for leftover gluon(s)
         DO 260 i=max(1,ip),n
C...Check colour charge
            mci=kchg(pycomp(k(i,2)),2)*isign(1,k(i,2))
            IF (k(i,1).NE.3.OR.mci.NE.2) goto 260
            IF(mcn(i,1).EQ.0.AND.mcn(i,2).EQ.0) THEN
C...Decide where to put left-over gluon (minimal insertion)
               ilmax=0
               rlmax=1d19
               DO 250 kct=nct+1,lct
                  DO 240 it=max(1,ip),n
                     IF (it.EQ.i.OR.k(it,1).NE.3) goto 240
                     IF (mcn(it,1).EQ.kct) ic=it
                     IF (mcn(it,2).EQ.kct) ia=it
  240             CONTINUE
                  rl=four(ic,i)*four(ia,i)
                  IF (rl.LT.rlmax) THEN
                     rlmax=rl
                     icmax=ic
                     iamax=ia
                  ENDIF
  250          CONTINUE
               lct=lct+1
               mcn(i,1)=mcn(icmax,1)
               mcn(i,2)=lct
               mcn(icmax,1)=lct
            ENDIF
  260    CONTINUE
         DO 270 i=max(1,ip),n
C...Do not erase parton shower colour history
            IF (k(i,1).NE.3) goto 270
C...Check colour charge
            mci=kchg(pycomp(k(i,2)),2)*isign(1,k(i,2))
            IF (mci.EQ.0) goto 270
            IF (mcn(i,1).NE.0) mct(i,1)=mcn(i,1)
            IF (mcn(i,2).NE.0) mct(i,2)=mcn(i,2)
  270    CONTINUE
      ENDIF
 
 9999 RETURN
      END