luedit.f

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

    
C*********************************************************************  
    
      SUBROUTINE luedit(MEDIT)  
    
C...Purpose: to perform global manipulations on the event record,   
C...in particular to exclude unstable or undetectable partons/particles.    
      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/ 
      dimension ns(2),pts(2),pls(2) 
    
C...Remove unwanted partons/particles.  
      IF((medit.GE.0.AND.medit.LE.3).OR.medit.EQ.5) THEN    
        imax=n  
        IF(mstu(2).GT.0) imax=mstu(2)   
        i1=max(1,mstu(1))-1 
        DO 110 i=max(1,mstu(1)),imax    
        IF(k(i,1).EQ.0.OR.k(i,1).GT.20) goto 110    
        IF(medit.EQ.1) THEN 
          IF(k(i,1).GT.10) goto 110 
        ELSEIF(medit.EQ.2) THEN 
          IF(k(i,1).GT.10) goto 110 
          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  
        ELSEIF(medit.EQ.3) THEN 
          IF(k(i,1).GT.10) goto 110 
          kc=lucomp(k(i,2)) 
          IF(kc.EQ.0) goto 110  
          IF(kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0) goto 110  
        ELSEIF(medit.EQ.5) THEN 
          IF(k(i,1).EQ.13.OR.k(i,1).EQ.14) goto 110 
          kc=lucomp(k(i,2)) 
          IF(kc.EQ.0) goto 110  
          IF(k(i,1).GE.11.AND.kchg(kc,2).EQ.0) goto 110 
        ENDIF   
    
C...Pack remaining partons/particles. Origin no longer known.   
        i1=i1+1 
        DO 100 j=1,5    
        k(i1,j)=k(i,j)  
        p(i1,j)=p(i,j)  
  100   v(i1,j)=v(i,j)  
        k(i1,3)=0   
  110   CONTINUE    
        n=i1    
    
C...Selective removal of class of entries. New position of retained.    
      ELSEIF(medit.GE.11.AND.medit.LE.15) THEN  
        i1=0    
        DO 120 i=1,n    
        k(i,3)=mod(k(i,3),mstu(5))  
        IF(medit.EQ.11.AND.k(i,1).LT.0) goto 120    
        IF(medit.EQ.12.AND.k(i,1).EQ.0) goto 120    
        IF(medit.EQ.13.AND.(k(i,1).EQ.11.OR.k(i,1).EQ.12.OR.    
     &  k(i,1).EQ.15).AND.k(i,2).NE.94) goto 120    
        IF(medit.EQ.14.AND.(k(i,1).EQ.13.OR.k(i,1).EQ.14.OR.    
     &  k(i,2).EQ.94)) goto 120 
        IF(medit.EQ.15.AND.k(i,1).GE.21) goto 120   
        i1=i1+1 
        k(i,3)=k(i,3)+mstu(5)*i1    
  120   CONTINUE    
    
C...Find new event history information and replace old. 
        DO 140 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GT.20.OR.k(i,3)/mstu(5).EQ.0) goto 140 
        id=i    
  130   im=mod(k(id,3),mstu(5)) 
        IF(medit.EQ.13.AND.im.GT.0.AND.im.LE.n) THEN    
          IF((k(im,1).EQ.11.OR.k(im,1).EQ.12.OR.k(im,1).EQ.15).AND. 
     &    k(im,2).NE.94) THEN   
            id=im   
            goto 130    
          ENDIF 
        ELSEIF(medit.EQ.14.AND.im.GT.0.AND.im.LE.n) THEN    
          IF(k(im,1).EQ.13.OR.k(im,1).EQ.14.OR.k(im,2).EQ.94) THEN  
            id=im   
            goto 130    
          ENDIF 
        ENDIF   
        k(i,3)=mstu(5)*(k(i,3)/mstu(5)) 
        IF(im.NE.0) k(i,3)=k(i,3)+k(im,3)/mstu(5)   
        IF(k(i,1).NE.3.AND.k(i,1).NE.13.AND.k(i,1).NE.14) THEN  
          IF(k(i,4).GT.0.AND.k(i,4).LE.mstu(4)) k(i,4)= 
     &    k(k(i,4),3)/mstu(5)   
          IF(k(i,5).GT.0.AND.k(i,5).LE.mstu(4)) k(i,5)= 
     &    k(k(i,5),3)/mstu(5)   
        ELSE    
          kcm=mod(k(i,4)/mstu(5),mstu(5))   
          IF(kcm.GT.0.AND.kcm.LE.mstu(4)) kcm=k(kcm,3)/mstu(5)  
          kcd=mod(k(i,4),mstu(5))   
          IF(kcd.GT.0.AND.kcd.LE.mstu(4)) kcd=k(kcd,3)/mstu(5)  
          k(i,4)=mstu(5)**2*(k(i,4)/mstu(5)**2)+mstu(5)*kcm+kcd 
          kcm=mod(k(i,5)/mstu(5),mstu(5))   
          IF(kcm.GT.0.AND.kcm.LE.mstu(4)) kcm=k(kcm,3)/mstu(5)  
          kcd=mod(k(i,5),mstu(5))   
          IF(kcd.GT.0.AND.kcd.LE.mstu(4)) kcd=k(kcd,3)/mstu(5)  
          k(i,5)=mstu(5)**2*(k(i,5)/mstu(5)**2)+mstu(5)*kcm+kcd 
        ENDIF   
  140   CONTINUE    
    
C...Pack remaining entries. 
        i1=0    
        DO 160 i=1,n    
        IF(k(i,3)/mstu(5).EQ.0) goto 160    
        i1=i1+1 
        DO 150 j=1,5    
        k(i1,j)=k(i,j)  
        p(i1,j)=p(i,j)  
  150   v(i1,j)=v(i,j)  
        k(i1,3)=mod(k(i1,3),mstu(5))    
  160   CONTINUE    
        n=i1    
    
C...Save top entries at bottom of LUJETS commonblock.   
      ELSEIF(medit.EQ.21) THEN  
        IF(2*n.GE.mstu(4)) THEN 
          CALL luerrm(11,'(LUEDIT:) no more memory left in LUJETS') 
          RETURN    
        ENDIF   
        DO 170 i=1,n    
        DO 170 j=1,5    
        k(mstu(4)-i,j)=k(i,j)   
        p(mstu(4)-i,j)=p(i,j)   
  170   v(mstu(4)-i,j)=v(i,j)   
        mstu(32)=n  
    
C...Restore bottom entries of commonblock LUJETS to top.    
      ELSEIF(medit.EQ.22) THEN  
        DO 180 i=1,mstu(32) 
        DO 180 j=1,5    
        k(i,j)=k(mstu(4)-i,j)   
        p(i,j)=p(mstu(4)-i,j)   
  180   v(i,j)=v(mstu(4)-i,j)   
        n=mstu(32)  
    
C...Mark primary entries at top of commonblock LUJETS as untreated. 
      ELSEIF(medit.EQ.23) THEN  
        i1=0    
        DO 190 i=1,n    
        kh=k(i,3)   
        IF(kh.GE.1) THEN    
          IF(k(kh,1).GT.20) kh=0    
        ENDIF   
        IF(kh.NE.0) goto 200    
        i1=i1+1 
  190   IF(k(i,1).GT.10.AND.k(i,1).LE.20) k(i,1)=k(i,1)-10  
  200   n=i1    
    
C...Place largest axis along z axis and second largest in xy plane. 
      ELSEIF(medit.EQ.31.OR.medit.EQ.32) THEN   
        CALL ludbrb(1,n+mstu(3),0.,-ulangl(p(mstu(61),1),   
     &  p(mstu(61),2)),0d0,0d0,0d0) 
        CALL ludbrb(1,n+mstu(3),-ulangl(p(mstu(61),3),  
     &  p(mstu(61),1)),0.,0d0,0d0,0d0)  
        CALL ludbrb(1,n+mstu(3),0.,-ulangl(p(mstu(61)+1,1), 
     &  p(mstu(61)+1,2)),0d0,0d0,0d0)   
        IF(medit.EQ.31) RETURN  
    
C...Rotate to put slim jet along +z axis.   
        DO 210 is=1,2   
        ns(is)=0    
        pts(is)=0.  
  210   pls(is)=0.  
        DO 220 i=1,n    
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 220    
        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 220    
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0) 
     &    goto 220  
        ENDIF   
        is=2.-sign(0.5,p(i,3))  
        ns(is)=ns(is)+1 
        pts(is)=pts(is)+sqrt(p(i,1)**2+p(i,2)**2)   
  220   CONTINUE    
        IF(ns(1)*pts(2)**2.LT.ns(2)*pts(1)**2)  
     &  CALL ludbrb(1,n+mstu(3),paru(1),0.,0d0,0d0,0d0) 
    
C...Rotate to put second largest jet into -z,+x quadrant.   
        DO 230 i=1,n    
        IF(p(i,3).GE.0.) goto 230   
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 230    
        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 230    
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0) 
     &    goto 230  
        ENDIF   
        is=2.-sign(0.5,p(i,1))  
        pls(is)=pls(is)-p(i,3)  
  230   CONTINUE    
        IF(pls(2).GT.pls(1)) CALL ludbrb(1,n+mstu(3),0.,paru(1),    
     &  0d0,0d0,0d0)    
      ENDIF 
    
      RETURN    
      END