luthru.f

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

    
C*********************************************************************  
    
      SUBROUTINE luthru(THR,OBL)    
    
C...Purpose: to perform thrust analysis to give thrust, oblateness  
C...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 tdi(3),tpr(3)   
    
C...Take copy of particles that are to be considered in thrust analysis.    
      np=0  
      ps=0. 
      DO 100 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 100  
      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 100  
        IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)   
     &  goto 100    
      ENDIF 
      IF(n+np+mstu(44)+15.GE.mstu(4)-mstu(32)-5) THEN   
        CALL luerrm(11,'(LUTHRU:) no more memory left in LUJETS')   
        thr=-2. 
        obl=-2. 
        RETURN  
      ENDIF 
      np=np+1   
      k(n+np,1)=23  
      p(n+np,1)=p(i,1)  
      p(n+np,2)=p(i,2)  
      p(n+np,3)=p(i,3)  
      p(n+np,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2) 
      p(n+np,5)=1.  
      IF(abs(paru(42)-1.).GT.0.001) p(n+np,5)=p(n+np,4)**(paru(42)-1.)  
      ps=ps+p(n+np,4)*p(n+np,5) 
  100 CONTINUE  
    
C...Very low multiplicities (0 or 1) not considered.    
      IF(np.LE.1) THEN  
        CALL luerrm(8,'(LUTHRU:) too few particles for analysis')   
        thr=-1. 
        obl=-1. 
        RETURN  
      ENDIF 
    
C...Loop over thrust and major. T axis along z direction in latter case.    
      DO 280 ild=1,2    
      IF(ild.EQ.2) THEN 
        k(n+np+1,1)=31  
        phi=ulangl(p(n+np+1,1),p(n+np+1,2)) 
        CALL ludbrb(n+1,n+np+1,0.,-phi,0d0,0d0,0d0) 
        the=ulangl(p(n+np+1,3),p(n+np+1,1)) 
        CALL ludbrb(n+1,n+np+1,-the,0.,0d0,0d0,0d0) 
      ENDIF 
    
C...Find and order particles with highest p (pT for major). 
      DO 110 ilf=n+np+4,n+np+mstu(44)+4 
  110 p(ilf,4)=0.   
      DO 150 i=n+1,n+np 
      IF(ild.EQ.2) p(i,4)=sqrt(p(i,1)**2+p(i,2)**2) 
      DO 120 ilf=n+np+mstu(44)+3,n+np+4,-1  
      IF(p(i,4).LE.p(ilf,4)) goto 130   
      DO 120 j=1,5  
  120 p(ilf+1,j)=p(ilf,j)   
      ilf=n+np+3    
  130 DO 140 j=1,5  
  140 p(ilf+1,j)=p(i,j) 
  150 CONTINUE  
    
C...Find and order initial axes with highest thrust (major).    
      DO 160 ilg=n+np+mstu(44)+5,n+np+mstu(44)+15   
  160 p(ilg,4)=0.   
      nc=2**(min(mstu(44),np)-1)    
      DO 220 ilc=1,nc   
      DO 170 j=1,3  
  170 tdi(j)=0. 
      DO 180 ilf=1,min(mstu(44),np) 
      sgn=p(n+np+ilf+3,5)   
      IF(2**ilf*((ilc+2**(ilf-1)-1)/2**ilf).GE.ilc) sgn=-sgn    
      DO 180 j=1,4-ild  
  180 tdi(j)=tdi(j)+sgn*p(n+np+ilf+3,j) 
      tds=tdi(1)**2+tdi(2)**2+tdi(3)**2 
      DO 190 ilg=n+np+mstu(44)+min(ilc,10)+4,n+np+mstu(44)+5,-1 
      IF(tds.LE.p(ilg,4)) goto 200  
      DO 190 j=1,4  
  190 p(ilg+1,j)=p(ilg,j)   
      ilg=n+np+mstu(44)+4   
  200 DO 210 j=1,3  
  210 p(ilg+1,j)=tdi(j) 
      p(ilg+1,4)=tds    
  220 CONTINUE  
    
C...Iterate direction of axis until stable maximum. 
      p(n+np+ild,4)=0.  
      ilg=0 
  230 ilg=ilg+1 
      thp=0.    
  240 thps=thp  
      DO 250 j=1,3  
      IF(thp.LE.1e-10) tdi(j)=p(n+np+mstu(44)+4+ilg,j)  
      IF(thp.GT.1e-10) tdi(j)=tpr(j)    
  250 tpr(j)=0. 
      DO 260 i=n+1,n+np 
      sgn=sign(p(i,5),tdi(1)*p(i,1)+tdi(2)*p(i,2)+tdi(3)*p(i,3))    
      DO 260 j=1,4-ild  
  260 tpr(j)=tpr(j)+sgn*p(i,j)  
      thp=sqrt(tpr(1)**2+tpr(2)**2+tpr(3)**2)/ps    
      IF(thp.GE.thps+paru(48)) goto 240 
    
C...Save good axis. Try new initial axis until a number of tries agree. 
      IF(thp.LT.p(n+np+ild,4)-paru(48).AND.ilg.LT.min(10,nc)) goto 230  
      IF(thp.GT.p(n+np+ild,4)+paru(48)) THEN    
        iagr=0  
        sgn=(-1.)**int(rlu(0)+0.5)  
        DO 270 j=1,3    
  270   p(n+np+ild,j)=sgn*tpr(j)/(ps*thp)   
        p(n+np+ild,4)=thp   
        p(n+np+ild,5)=0.    
      ENDIF 
      iagr=iagr+1   
  280 IF(iagr.LT.mstu(45).AND.ilg.LT.min(10,nc)) goto 230   
    
C...Find minor axis and value by orthogonality. 
      sgn=(-1.)**int(rlu(0)+0.5)    
      p(n+np+3,1)=-sgn*p(n+np+2,2)  
      p(n+np+3,2)=sgn*p(n+np+2,1)   
      p(n+np+3,3)=0.    
      thp=0.    
      DO 290 i=n+1,n+np 
  290 thp=thp+p(i,5)*abs(p(n+np+3,1)*p(i,1)+p(n+np+3,2)*p(i,2)) 
      p(n+np+3,4)=thp/ps    
      p(n+np+3,5)=0.    
    
C...Fill axis information. Rotate back to original coordinate system.   
      DO 300 ild=1,3    
      k(n+ild,1)=31 
      k(n+ild,2)=96 
      k(n+ild,3)=ild    
      k(n+ild,4)=0  
      k(n+ild,5)=0  
      DO 300 j=1,5  
      p(n+ild,j)=p(n+np+ild,j)  
  300 v(n+ild,j)=0. 
      CALL ludbrb(n+1,n+3,the,phi,0d0,0d0,0d0)  
    
C...Select storing option. Calculate thurst and oblateness. 
      mstu(61)=n+1  
      mstu(62)=np   
      IF(mstu(43).LE.1) mstu(3)=3   
      IF(mstu(43).GE.2) n=n+3   
      thr=p(n+1,4)  
      obl=p(n+2,4)-p(n+3,4) 
    
      RETURN    
      END