lusphe.f

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

    
C*********************************************************************  
    
      SUBROUTINE lusphe(SPH,APL)    
    
C...Purpose: to perform sphericity tensor analysis to give sphericity,  
C...aplanarity and the related event axes.  
      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),sv(3,3) 
    
C...Calculate matrix to be diagonalized.    
      np=0  
      DO 100 j1=1,3 
      DO 100 j2=j1,3    
  100 sm(j1,j2)=0.  
      ps=0. 
      DO 120 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120  
      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 120  
        IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)   
     &  goto 120    
      ENDIF 
      np=np+1   
      pa=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)    
      pwt=1.    
      IF(abs(paru(41)-2.).GT.0.001) pwt=max(1e-10,pa)**(paru(41)-2.)    
      DO 110 j1=1,3 
      DO 110 j2=j1,3    
  110 sm(j1,j2)=sm(j1,j2)+pwt*p(i,j1)*p(i,j2)   
      ps=ps+pwt*pa**2   
  120 CONTINUE  
    
C...Very low multiplicities (0 or 1) not considered.    
      IF(np.LE.1) THEN  
        CALL luerrm(8,'(LUSPHE:) too few particles for analysis')   
        sph=-1. 
        apl=-1. 
        RETURN  
      ENDIF 
      DO 130 j1=1,3 
      DO 130 j2=j1,3    
  130 sm(j1,j2)=sm(j1,j2)/ps    
    
C...Find eigenvalues to matrix (third degree equation). 
      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.) 
      p(n+1,4)=1./3.+sqrt(-sq)*max(2.*sp,sqrt(3.*(1.-sp**2))-sp)    
      p(n+3,4)=1./3.+sqrt(-sq)*min(2.*sp,-sqrt(3.*(1.-sp**2))-sp)   
      p(n+2,4)=1.-p(n+1,4)-p(n+3,4) 
      IF(p(n+2,4).LT.1e-5) THEN 
        CALL luerrm(8,'(LUSPHE:) all particles back-to-back')   
        sph=-1. 
        apl=-1. 
        RETURN  
      ENDIF 
    
C...Find first and last eigenvector by solving equation system. 
      DO 170 i=1,3,2    
      DO 140 j1=1,3 
      sv(j1,j1)=sm(j1,j1)-p(n+i,4)  
      DO 140 j2=j1+1,3  
      sv(j1,j2)=sm(j1,j2)   
  140 sv(j2,j1)=sm(j1,j2)   
      smax=0.   
      DO 150 j1=1,3 
      DO 150 j2=1,3 
      IF(abs(sv(j1,j2)).LE.smax) goto 150   
      ja=j1 
      jb=j2 
      smax=abs(sv(j1,j2))   
  150 CONTINUE  
      smax=0.   
      DO 160 j3=ja+1,ja+2   
      j1=j3-3*((j3-1)/3)    
      rl=sv(j1,jb)/sv(ja,jb)    
      DO 160 j2=1,3 
      sv(j1,j2)=sv(j1,j2)-rl*sv(ja,j2)  
      IF(abs(sv(j1,j2)).LE.smax) goto 160   
      jc=j1 
      smax=abs(sv(j1,j2))   
  160 CONTINUE  
      jb1=jb+1-3*(jb/3) 
      jb2=jb+2-3*((jb+1)/3) 
      p(n+i,jb1)=-sv(jc,jb2)    
      p(n+i,jb2)=sv(jc,jb1) 
      p(n+i,jb)=-(sv(ja,jb1)*p(n+i,jb1)+sv(ja,jb2)*p(n+i,jb2))/ 
     &sv(ja,jb) 
      pa=sqrt(p(n+i,1)**2+p(n+i,2)**2+p(n+i,3)**2)  
      sgn=(-1.)**int(rlu(0)+0.5)    
      DO 170 j=1,3  
  170 p(n+i,j)=sgn*p(n+i,j)/pa  
    
C...Middle axis orthogonal to other two. Fill other codes.  
      sgn=(-1.)**int(rlu(0)+0.5)    
      p(n+2,1)=sgn*(p(n+1,2)*p(n+3,3)-p(n+1,3)*p(n+3,2))    
      p(n+2,2)=sgn*(p(n+1,3)*p(n+3,1)-p(n+1,1)*p(n+3,3))    
      p(n+2,3)=sgn*(p(n+1,1)*p(n+3,2)-p(n+1,2)*p(n+3,1))    
      DO 180 i=1,3  
      k(n+i,1)=31   
      k(n+i,2)=95   
      k(n+i,3)=i    
      k(n+i,4)=0    
      k(n+i,5)=0    
      p(n+i,5)=0.   
      DO 180 j=1,5  
  180 v(i,j)=0. 
    
C...Select storing option. Calculate sphericity and aplanarity. 
      mstu(61)=n+1  
      mstu(62)=np   
      IF(mstu(43).LE.1) mstu(3)=3   
      IF(mstu(43).GE.2) n=n+3   
      sph=1.5*(p(n+2,4)+p(n+3,4))   
      apl=1.5*p(n+3,4)  
    
      RETURN    
      END