pycell.f

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

 
C*********************************************************************
 
C...PYCELL
C...Provides a simple way of jet finding in eta-phi-ET coordinates,
C...as used for calorimeters at hadron colliders.
 
      SUBROUTINE pycell(NJET)
 
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)
      SAVE /pyjets/,/pydat1/,/pydat2/
 
C...Loop over all particles. Find cell that was hit by given particle.
      ptlrat=1d0/sinh(paru(51))**2
      np=0
      nc=n
      DO 110 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 110
        IF(p(i,1)**2+p(i,2)**2.LE.ptlrat*p(i,3)**2) goto 110
        IF(mstu(41).GE.2) THEN
          kc=pycomp(k(i,2))
          IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
     &    kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
     &    k(i,2).EQ.ksusy1+39) goto 110
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0)
     &    goto 110
        ENDIF
        np=np+1
        pt=sqrt(p(i,1)**2+p(i,2)**2)
        eta=sign(log((sqrt(pt**2+p(i,3)**2)+abs(p(i,3)))/pt),p(i,3))
        ieta=max(1,min(mstu(51),1+int(mstu(51)*0.5d0*
     &  (eta/paru(51)+1d0))))
        phi=pyangl(p(i,1),p(i,2))
        iphi=max(1,min(mstu(52),1+int(mstu(52)*0.5d0*
     &  (phi/paru(1)+1d0))))
        ietph=mstu(52)*ieta+iphi
 
C...Add to cell already hit, or book new cell.
        DO 100 ic=n+1,nc
          IF(ietph.EQ.k(ic,3)) THEN
            k(ic,4)=k(ic,4)+1
            p(ic,5)=p(ic,5)+pt
            goto 110
          ENDIF
  100   CONTINUE
        IF(nc.GE.mstu(4)-mstu(32)-5) THEN
          CALL pyerrm(11,'(PYCELL:) no more memory left in PYJETS')
          njet=-2
          RETURN
        ENDIF
        nc=nc+1
        k(nc,3)=ietph
        k(nc,4)=1
        k(nc,5)=2
        p(nc,1)=(paru(51)/mstu(51))*(2*ieta-1-mstu(51))
        p(nc,2)=(paru(1)/mstu(52))*(2*iphi-1-mstu(52))
        p(nc,5)=pt
  110 CONTINUE
 
C...Smear true bin content by calorimeter resolution.
      IF(mstu(53).GE.1) THEN
        DO 130 ic=n+1,nc
          pei=p(ic,5)
          IF(mstu(53).EQ.2) pei=p(ic,5)*cosh(p(ic,1))
  120     pef=pei+paru(55)*sqrt(-2d0*log(max(1d-10,pyr(0)))*pei)*
     &    cos(paru(2)*pyr(0))
          IF(pef.LT.0d0.OR.pef.GT.paru(56)*pei) goto 120
          p(ic,5)=pef
          IF(mstu(53).EQ.2) p(ic,5)=pef/cosh(p(ic,1))
  130   CONTINUE
      ENDIF
 
C...Remove cells below threshold.
      IF(paru(58).GT.0d0) THEN
        ncc=nc
        nc=n
        DO 140 ic=n+1,ncc
          IF(p(ic,5).GT.paru(58)) THEN
            nc=nc+1
            k(nc,3)=k(ic,3)
            k(nc,4)=k(ic,4)
            k(nc,5)=k(ic,5)
            p(nc,1)=p(ic,1)
            p(nc,2)=p(ic,2)
            p(nc,5)=p(ic,5)
          ENDIF
  140   CONTINUE
      ENDIF
 
C...Find initiator cell: the one with highest pT of not yet used ones.
      nj=nc
  150 etmax=0d0
      DO 160 ic=n+1,nc
        IF(k(ic,5).NE.2) goto 160
        IF(p(ic,5).LE.etmax) goto 160
        icmax=ic
        eta=p(ic,1)
        phi=p(ic,2)
        etmax=p(ic,5)
  160 CONTINUE
      IF(etmax.LT.paru(52)) goto 220
      IF(nj.GE.mstu(4)-mstu(32)-5) THEN
        CALL pyerrm(11,'(PYCELL:) no more memory left in PYJETS')
        njet=-2
        RETURN
      ENDIF
      k(icmax,5)=1
      nj=nj+1
      k(nj,4)=0
      k(nj,5)=1
      p(nj,1)=eta
      p(nj,2)=phi
      p(nj,3)=0d0
      p(nj,4)=0d0
      p(nj,5)=0d0
 
C...Sum up unused cells within required distance of initiator.
      DO 170 ic=n+1,nc
        IF(k(ic,5).EQ.0) goto 170
        IF(abs(p(ic,1)-eta).GT.paru(54)) goto 170
        dphia=abs(p(ic,2)-phi)
        IF(dphia.GT.paru(54).AND.dphia.LT.paru(2)-paru(54)) goto 170
        phic=p(ic,2)
        IF(dphia.GT.paru(1)) phic=phic+sign(paru(2),phi)
        IF((p(ic,1)-eta)**2+(phic-phi)**2.GT.paru(54)**2) goto 170
        k(ic,5)=-k(ic,5)
        k(nj,4)=k(nj,4)+k(ic,4)
        p(nj,3)=p(nj,3)+p(ic,5)*p(ic,1)
        p(nj,4)=p(nj,4)+p(ic,5)*phic
        p(nj,5)=p(nj,5)+p(ic,5)
  170 CONTINUE
 
C...Reject cluster below minimum ET, else accept.
      IF(p(nj,5).LT.paru(53)) THEN
        nj=nj-1
        DO 180 ic=n+1,nc
          IF(k(ic,5).LT.0) k(ic,5)=-k(ic,5)
  180   CONTINUE
      ELSEIF(mstu(54).LE.2) THEN
        p(nj,3)=p(nj,3)/p(nj,5)
        p(nj,4)=p(nj,4)/p(nj,5)
        IF(abs(p(nj,4)).GT.paru(1)) p(nj,4)=p(nj,4)-sign(paru(2),
     &  p(nj,4))
        DO 190 ic=n+1,nc
          IF(k(ic,5).LT.0) k(ic,5)=0
  190   CONTINUE
      ELSE
        DO 200 j=1,4
          p(nj,j)=0d0
  200   CONTINUE
        DO 210 ic=n+1,nc
          IF(k(ic,5).GE.0) goto 210
          p(nj,1)=p(nj,1)+p(ic,5)*cos(p(ic,2))
          p(nj,2)=p(nj,2)+p(ic,5)*sin(p(ic,2))
          p(nj,3)=p(nj,3)+p(ic,5)*sinh(p(ic,1))
          p(nj,4)=p(nj,4)+p(ic,5)*cosh(p(ic,1))
          k(ic,5)=0
  210   CONTINUE
      ENDIF
      goto 150
 
C...Arrange clusters in falling ET sequence.
  220 DO 250 i=1,nj-nc
        etmax=0d0
        DO 230 ij=nc+1,nj
          IF(k(ij,5).EQ.0) goto 230
          IF(p(ij,5).LT.etmax) goto 230
          ijmax=ij
          etmax=p(ij,5)
  230   CONTINUE
        k(ijmax,5)=0
        k(n+i,1)=31
        k(n+i,2)=98
        k(n+i,3)=i
        k(n+i,4)=k(ijmax,4)
        k(n+i,5)=0
        DO 240 j=1,5
          p(n+i,j)=p(ijmax,j)
          v(n+i,j)=0d0
  240   CONTINUE
  250 CONTINUE
      njet=nj-nc
 
C...Convert to massless or massive four-vectors.
      IF(mstu(54).EQ.2) THEN
        DO 260 i=n+1,n+njet
          eta=p(i,3)
          p(i,1)=p(i,5)*cos(p(i,4))
          p(i,2)=p(i,5)*sin(p(i,4))
          p(i,3)=p(i,5)*sinh(eta)
          p(i,4)=p(i,5)*cosh(eta)
          p(i,5)=0d0
  260   CONTINUE
      ELSEIF(mstu(54).GE.3) THEN
        DO 270 i=n+1,n+njet
          p(i,5)=sqrt(max(0d0,p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2))
  270   CONTINUE
      ENDIF
 
C...Information about storage.
      mstu(61)=n+1
      mstu(62)=np
      mstu(63)=nc-n
      IF(mstu(43).LE.1) mstu(3)=max(0,njet)
      IF(mstu(43).GE.2) n=n+max(0,njet)
 
      RETURN
      END