pyadsh.f

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

 
C*********************************************************************
 
C...PYADSH
C...Administers the generation of successive final-state showers
C...in external processes.
 
      SUBROUTINE pyadsh(NFIN)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement for maximum size of showers.
      parameter(maxnur=1000)
C...Commonblocks.
      common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pyctag/nct,mct(4000,2)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pypart/,/pyjets/,/pyctag/,/pydat1/,/pypars/,/pyint1/
C...Local array.
      dimension ibeg(100),ksav(100,5),psum(4),beta(3)
 
C...Set primary vertex.
      DO 100 j=1,5
        v(mint(83)+5,j)=0d0
        v(mint(83)+6,j)=0d0
        v(mint(84)+1,j)=0d0
        v(mint(84)+2,j)=0d0
  100 CONTINUE
 
C...Isolate systems of particles with the same mother.
      nsys=0
      ims=-1
      DO 140 i=mint(84)+3,nfin
        im=k(i,3)
        IF(im.GT.0.AND.im.LE.mint(84)) im=k(im,3)
        IF(im.NE.ims) THEN
          nsys=nsys+1
          ibeg(nsys)=i
          ims=im
        ENDIF
 
C...Set production vertices.
        IF(im.LE.mint(83)+6.OR.(im.GT.mint(84).AND.im.LE.mint(84)+2))
     &  THEN
          DO 110 j=1,4
            v(i,j)=0d0
  110     CONTINUE
        ELSE
          DO 120 j=1,4
            v(i,j)=v(im,j)+v(im,5)*p(im,j)/p(im,5)
  120     CONTINUE
        ENDIF
        IF(mstp(125).GE.1) THEN
          idoc=i-mstp(126)+4
          DO 130 j=1,5
            v(idoc,j)=v(i,j)
  130     CONTINUE
        ENDIF
  140 CONTINUE
 
C...End loop over systems. Return if no showers to be performed.
      ibeg(nsys+1)=nfin+1
      IF(mstp(71).LE.0) RETURN
 
C...Loop through systems of particles; check that sensible size.
      DO 270 isys=1,nsys
        nsiz=ibeg(isys+1)-ibeg(isys)
        IF(mint(35).LE.1) THEN
          IF(nsiz.EQ.1.AND.isys.EQ.1) THEN
            goto 270
          ELSEIF(nsiz.LE.1) THEN
            CALL pyerrm(2,'(PYADSH:) only one particle in system')
            goto 270
          ELSEIF(nsiz.GT.80) THEN
            CALL pyerrm(2,'(PYADSH:) more than 80 particles in system')
            goto 270
          ENDIF
        ENDIF
 
C...Save status codes and daughters of showering particles; reset them.
        DO 150 j=1,4
          psum(j)=0d0
  150   CONTINUE
        DO 170 ii=1,nsiz
          i=ibeg(isys)-1+ii
          ksav(ii,1)=k(i,1)
          IF(k(i,1).GT.10) THEN
            k(i,1)=1
            IF(ksav(ii,1).EQ.14) k(i,1)=3
          ENDIF
          IF(ksav(ii,1).LE.10) THEN
          ELSEIF(k(i,1).EQ.1) THEN
            ksav(ii,4)=k(i,4)
            ksav(ii,5)=k(i,5)
            k(i,4)=0
            k(i,5)=0
          ELSE
            ksav(ii,4)=mod(k(i,4),mstu(5))
            ksav(ii,5)=mod(k(i,5),mstu(5))
            k(i,4)=k(i,4)-ksav(ii,4)
            k(i,5)=k(i,5)-ksav(ii,5)
          ENDIF
          DO 160 j=1,4
            psum(j)=psum(j)+p(i,j)
  160     CONTINUE
  170   CONTINUE
 
C...Perform shower.
        qmax=sqrt(max(0d0,psum(4)**2-psum(1)**2-psum(2)**2-
     &  psum(3)**2))
        IF(isys.EQ.1) qmax=min(qmax,sqrt(parp(71))*vint(55))
        nsav=n
        IF(mint(35).LE.1) THEN
          IF(nsiz.EQ.2) THEN
            CALL pyshow(ibeg(isys),ibeg(isys)+1,qmax)
          ELSE
            CALL pyshow(ibeg(isys),-nsiz,qmax)
          ENDIF
 
C...For external processes, first call, also ISR partons radiate.
C...Can use existing PYPART list, removing partons that radiate later.
        ELSEIF(isys.EQ.1) THEN
          npartn=0
          DO 175 ii=1,npart
            IF(ipart(ii).LT.ibeg(2).OR.ipart(ii).GE.ibeg(nsys+1)) THEN
              npartn=npartn+1
              ipart(npartn)=ipart(ii)
              ptpart(npartn)=ptpart(ii)
            ENDIF
 175      CONTINUE
          npart=npartn
          CALL pyptfs(1,0.5d0*qmax,0d0,ptgen)
        ELSE
C...For subsequent calls use the systems excluded above.
          npart=nsiz
          npartd=0
          DO 180 ii=1,nsiz
            i=ibeg(isys)-1+ii
            ipart(ii)=i
            ptpart(ii)=0.5d0*qmax
  180     CONTINUE
          CALL pyptfs(2,0.5d0*qmax,0d0,ptgen)
        ENDIF
 
C...Look up showered copies of original showering particles.
        DO 260 ii=1,nsiz
          i=ibeg(isys)-1+ii
          imv=i
C...Particles without daughters need not be studied.
          IF(ksav(ii,1).LE.10) goto 260
          IF(n.EQ.nsav.OR.k(i,1).LE.10) THEN
          ELSEIF(k(i,1).EQ.11) THEN
  190       imv=mod(k(imv,4),mstu(5))
            IF(k(imv,1).EQ.11) goto 190
          ELSE
            kda1=mod(k(i,4),mstu(5))
            IF(kda1.GT.0) THEN
              IF(k(kda1,2).EQ.21) kda1=k(kda1,5)/mstu(5)
            ENDIF
            kda2=mod(k(i,5),mstu(5))
            IF(kda2.GT.0) THEN
              IF(k(kda2,2).EQ.21) kda2=k(kda2,4)/mstu(5)
            ENDIF
            DO 200 i3=i+1,n
              IF(k(i3,2).EQ.k(i,2).AND.(i3.EQ.kda1.OR.i3.EQ.kda2))
     &        THEN
                imv=i3
                kda1=mod(k(i3,4),mstu(5))
                IF(kda1.GT.0) THEN
                  IF(k(kda1,2).EQ.21) kda1=k(kda1,5)/mstu(5)
                ENDIF
                kda2=mod(k(i3,5),mstu(5))
                IF(kda2.GT.0) THEN
                  IF(k(kda2,2).EQ.21) kda2=k(kda2,4)/mstu(5)
                ENDIF
              ENDIF
  200       CONTINUE
          ENDIF
 
C...Restore daughter info of original partons to showered copies.
          IF(ksav(ii,1).GT.10) k(imv,1)=ksav(ii,1)
          IF(ksav(ii,1).LE.10) THEN
          ELSEIF(k(i,1).EQ.1) THEN
            k(imv,4)=ksav(ii,4)
            k(imv,5)=ksav(ii,5)
          ELSE
            k(imv,4)=k(imv,4)+ksav(ii,4)
            k(imv,5)=k(imv,5)+ksav(ii,5)
          ENDIF
 
C...Reset mother info of existing daughters to showered copies.
          DO 210 i3=ibeg(isys+1),nfin
            IF(k(i3,3).EQ.i) k(i3,3)=imv
            IF(k(i3,1).EQ.3.OR.k(i3,1).EQ.14) THEN
              IF(k(i3,4)/mstu(5).EQ.i) k(i3,4)=k(i3,4)+mstu(5)*(imv-i)
              IF(k(i3,5)/mstu(5).EQ.i) k(i3,5)=k(i3,5)+mstu(5)*(imv-i)
            ENDIF
  210     CONTINUE
 
C...Boost all original daughters to new frame of showered copy.
C...Also update their colour tags.
          IF(imv.NE.i) THEN
            DO 220 j=1,3
              beta(j)=(p(imv,j)-p(i,j))/(p(imv,4)+p(i,4))
  220       CONTINUE
            fac=2d0/(1d0+beta(1)**2+beta(2)**2+beta(3)**2)
            DO 230 j=1,3
              beta(j)=fac*beta(j)
  230       CONTINUE
            DO 250 i3=ibeg(isys+1),nfin
              imo=i3
  240         imo=k(imo,3)
              IF(mstp(128).LE.0) THEN
                IF(imo.GT.0.AND.imo.NE.i.AND.imo.NE.k(i,3)) goto 240
                IF(imo.EQ.i.OR.(k(i,3).LE.mint(84).AND.imo.EQ.k(i,3)))
     &          THEN
                  CALL pyrobo(i3,i3,0d0,0d0,beta(1),beta(2),beta(3))
                  IF(mct(i3,1).EQ.mct(i,1)) mct(i3,1)=mct(imv,1)
                  IF(mct(i3,2).EQ.mct(i,2)) mct(i3,2)=mct(imv,2)
                ENDIF
              ELSE
                IF(imo.EQ.imv) THEN
                  CALL pyrobo(i3,i3,0d0,0d0,beta(1),beta(2),beta(3))
                  IF(mct(i3,1).EQ.mct(i,1)) mct(i3,1)=mct(imv,1)
                  IF(mct(i3,2).EQ.mct(i,2)) mct(i3,2)=mct(imv,2)
                ELSEIF(imo.GT.0.AND.imo.NE.i.AND.imo.NE.k(i,3)) THEN
                  goto 240
                ENDIF
              ENDIF
  250       CONTINUE
          ENDIF
  260   CONTINUE
 
C...End of loop over showering systems
  270 CONTINUE
 
      RETURN
      END