lujmas.f

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

    
C*********************************************************************  
    
      SUBROUTINE lujmas(PMH,PML)    
    
C...Purpose: to determine, approximately, the two jet masses that   
C...minimize the sum m_H|2 + m_L|2, a la Clavelli and Wyler.    
      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 sm(3,3),sax(3),ps(3,5)  
    
C...Reset.  
      np=0  
      DO 110 j1=1,3 
      DO 100 j2=j1,3    
  100 sm(j1,j2)=0.  
      DO 110 j2=1,4 
  110 ps(j1,j2)=0.  
      pss=0.    
    
C...Take copy of particles that are to be considered in mass analysis.  
      DO 150 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 150  
      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 150  
        IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)   
     &  goto 150    
      ENDIF 
      IF(n+np+1.GE.mstu(4)-mstu(32)-5) THEN 
        CALL luerrm(11,'(LUJMAS:) no more memory left in LUJETS')   
        pmh=-2. 
        pml=-2. 
        RETURN  
      ENDIF 
      np=np+1   
      DO 120 j=1,5  
  120 p(n+np,j)=p(i,j)  
      IF(mstu(42).EQ.0) p(n+np,5)=0.    
      IF(mstu(42).EQ.1.AND.k(i,2).NE.22) p(n+np,5)=pmas(101,1)  
      p(n+np,4)=sqrt(p(n+np,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)    
    
C...Fill information in sphericity tensor and total momentum vector.    
      DO 130 j1=1,3 
      DO 130 j2=j1,3    
  130 sm(j1,j2)=sm(j1,j2)+p(i,j1)*p(i,j2)   
      pss=pss+(p(i,1)**2+p(i,2)**2+p(i,3)**2)   
      DO 140 j=1,4  
  140 ps(3,j)=ps(3,j)+p(n+np,j) 
  150 CONTINUE  
    
C...Very low multiplicities (0 or 1) not considered.    
      IF(np.LE.1) THEN  
        CALL luerrm(8,'(LUJMAS:) too few particles for analysis')   
        pmh=-1. 
        pml=-1. 
        RETURN  
      ENDIF 
      paru(61)=sqrt(max(0.,ps(3,4)**2-ps(3,1)**2-ps(3,2)**2-ps(3,3)**2))    
    
C...Find largest eigenvalue to matrix (third degree equation).  
      DO 160 j1=1,3 
      DO 160 j2=j1,3    
  160 sm(j1,j2)=sm(j1,j2)/pss   
      sq=(sm(1,1)*sm(2,2)+sm(1,1)*sm(3,3)+sm(2,2)*sm(3,3)-sm(1,2)**2-   
     &sm(1,3)**2-sm(2,3)**2)/3.-1./9.   
      sr=-0.5*(sq+1./9.+sm(1,1)*sm(2,3)**2+sm(2,2)*sm(1,3)**2+sm(3,3)*  
     &sm(1,2)**2-sm(1,1)*sm(2,2)*sm(3,3))+sm(1,2)*sm(1,3)*sm(2,3)+1./27.    
      sp=cos(acos(max(min(sr/sqrt(-sq**3),1.),-1.))/3.) 
      sma=1./3.+sqrt(-sq)*max(2.*sp,sqrt(3.*(1.-sp**2))-sp) 
    
C...Find largest eigenvector by solving equation system.    
      DO 170 j1=1,3 
      sm(j1,j1)=sm(j1,j1)-sma   
      DO 170 j2=j1+1,3  
  170 sm(j2,j1)=sm(j1,j2)   
      smax=0.   
      DO 180 j1=1,3 
      DO 180 j2=1,3 
      IF(abs(sm(j1,j2)).LE.smax) goto 180   
      ja=j1 
      jb=j2 
      smax=abs(sm(j1,j2))   
  180 CONTINUE  
      smax=0.   
      DO 190 j3=ja+1,ja+2   
      j1=j3-3*((j3-1)/3)    
      rl=sm(j1,jb)/sm(ja,jb)    
      DO 190 j2=1,3 
      sm(j1,j2)=sm(j1,j2)-rl*sm(ja,j2)  
      IF(abs(sm(j1,j2)).LE.smax) goto 190   
      jc=j1 
      smax=abs(sm(j1,j2))   
  190 CONTINUE  
      jb1=jb+1-3*(jb/3) 
      jb2=jb+2-3*((jb+1)/3) 
      sax(jb1)=-sm(jc,jb2)  
      sax(jb2)=sm(jc,jb1)   
      sax(jb)=-(sm(ja,jb1)*sax(jb1)+sm(ja,jb2)*sax(jb2))/sm(ja,jb)  
    
C...Divide particles into two initial clusters by hemisphere.   
      DO 200 i=n+1,n+np 
      psax=p(i,1)*sax(1)+p(i,2)*sax(2)+p(i,3)*sax(3)    
      is=1  
      IF(psax.LT.0.) is=2   
      k(i,3)=is 
      DO 200 j=1,4  
  200 ps(is,j)=ps(is,j)+p(i,j)  
      pms=(ps(1,4)**2-ps(1,1)**2-ps(1,2)**2-ps(1,3)**2)+    
     &(ps(2,4)**2-ps(2,1)**2-ps(2,2)**2-ps(2,3)**2) 
    
C...Reassign one particle at a time; find maximum decrease of m|2 sum.  
  210 pmd=0.    
      im=0  
      DO 220 j=1,4  
  220 ps(3,j)=ps(1,j)-ps(2,j)   
      DO 230 i=n+1,n+np 
      pps=p(i,4)*ps(3,4)-p(i,1)*ps(3,1)-p(i,2)*ps(3,2)-p(i,3)*ps(3,3)   
      IF(k(i,3).EQ.1) pmdi=2.*(p(i,5)**2-pps)   
      IF(k(i,3).EQ.2) pmdi=2.*(p(i,5)**2+pps)   
      IF(pmdi.LT.pmd) THEN  
        pmd=pmdi    
        im=i    
      ENDIF 
  230 CONTINUE  
    
C...Loop back if significant reduction in sum of m|2.   
      IF(pmd.LT.-paru(48)*pms) THEN 
        pms=pms+pmd 
        is=k(im,3)  
        DO 240 j=1,4    
        ps(is,j)=ps(is,j)-p(im,j)   
  240   ps(3-is,j)=ps(3-is,j)+p(im,j)   
        k(im,3)=3-is    
        goto 210    
      ENDIF 
    
C...Final masses and output.    
      mstu(61)=n+1  
      mstu(62)=np   
      ps(1,5)=sqrt(max(0.,ps(1,4)**2-ps(1,1)**2-ps(1,2)**2-ps(1,3)**2)) 
      ps(2,5)=sqrt(max(0.,ps(2,4)**2-ps(2,1)**2-ps(2,2)**2-ps(2,3)**2)) 
      pmh=max(ps(1,5),ps(2,5))  
      pml=min(ps(1,5),ps(2,5))  
    
      RETURN    
      END