lucell.f

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

    
C*********************************************************************  
    
      SUBROUTINE lucell(NJET)   
    
C...Purpose: to provide a simple way of jet finding in an eta-phi-ET    
C...coordinate frame, as used for calorimeters at hadron colliders. 
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
    
C...Loop over all particles. Find cell that was hit by given particle.  
      nce2=2*mstu(51)*mstu(52)  
      ptlrat=1./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=lucomp(k(i,2))   
        IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.  
     &  kc.EQ.18) goto 110  
        IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(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.5*(eta/paru(51)+1.))))   
      phi=ulangl(p(i,1),p(i,2)) 
      iphi=max(1,min(mstu(52),1+int(mstu(52)*0.5*(phi/paru(1)+1.))))    
      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 luerrm(11,'(LUCELL:) no more memory left in LUJETS')   
        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(-2.*log(max(1e-10,rlu(0)))*pei)*  
     &  cos(paru(2)*rlu(0)) 
        IF(pef.LT.0..OR.pef.GT.paru(56)*pei) goto 120   
        p(ic,5)=pef 
  130   IF(mstu(53).EQ.2) p(ic,5)=pef*cosh(p(ic,1)) 
      ENDIF 
    
C...Find initiator cell: the one with highest pT of not yet used ones.  
      nj=nc 
  140 etmax=0.  
      DO 150 ic=n+1,nc  
      IF(k(ic,5).NE.2) goto 150 
      IF(p(ic,5).LE.etmax) goto 150 
      icmax=ic  
      eta=p(ic,1)   
      phi=p(ic,2)   
      etmax=p(ic,5) 
  150 CONTINUE  
      IF(etmax.LT.paru(52)) goto 210    
      IF(nj.GE.mstu(4)-mstu(32)-5) THEN 
        CALL luerrm(11,'(LUCELL:) no more memory left in LUJETS')   
        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)=0.    
      p(nj,4)=0.    
      p(nj,5)=0.    
    
C...Sum up unused cells within required distance of initiator.  
      DO 160 ic=n+1,nc  
      IF(k(ic,5).EQ.0) goto 160 
      IF(abs(p(ic,1)-eta).GT.paru(54)) goto 160 
      dphia=abs(p(ic,2)-phi)    
      IF(dphia.GT.paru(54).AND.dphia.LT.paru(2)-paru(54)) goto 160  
      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 160    
      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)   
  160 CONTINUE  
    
C...Reject cluster below minimum ET, else accept.   
      IF(p(nj,5).LT.paru(53)) THEN  
        nj=nj-1 
        DO 170 ic=n+1,nc    
  170   IF(k(ic,5).LT.0) k(ic,5)=-k(ic,5)   
      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 180 ic=n+1,nc    
  180   IF(k(ic,1).LT.0) k(ic,1)=0  
      ELSE  
        DO 190 j=1,4    
  190   p(nj,j)=0.  
        DO 200 ic=n+1,nc    
        IF(k(ic,5).GE.0) goto 200   
        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   
  200   CONTINUE    
      ENDIF 
      goto 140  
    
C...Arrange clusters in falling ET sequence.    
  210 DO 230 i=1,nj-nc  
      etmax=0.  
      DO 220 ij=nc+1,nj 
      IF(k(ij,5).EQ.0) goto 220 
      IF(p(ij,5).LT.etmax) goto 220 
      ijmax=ij  
      etmax=p(ij,5) 
  220 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 230 j=1,5  
      p(n+i,j)=p(ijmax,j)   
  230 v(n+i,j)=0.   
      njet=nj-nc    
    
C...Convert to massless or massive four-vectors.    
      IF(mstu(54).EQ.2) THEN    
        DO 240 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) 
  240   p(i,5)=0.   
      ELSEIF(mstu(54).GE.3) THEN    
        DO 250 i=n+1,n+njet 
  250   p(i,5)=sqrt(max(0.,p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2))    
      ENDIF 
    
C...Information about storage.  
      mstu(61)=n+1  
      mstu(62)=np   
      mstu(63)=nc-n 
      IF(mstu(43).LE.1) mstu(3)=njet    
      IF(mstu(43).GE.2) n=n+njet    
    
      RETURN    
      END