lufowo.f

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

    
C*********************************************************************  
    
      SUBROUTINE lufowo(H10,H20,H30,H40)    
    
C...Purpose: to calculate the first few Fox-Wolfram moments.    
      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...Copy momenta for particles and calculate H0.    
      np=0  
      h0=0. 
      hd=0. 
      DO 110 i=1,n  
      IF(k(i,1).LE.0.OR.k(i,1).GT.10) 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 
      IF(n+np.GE.mstu(4)-mstu(32)-5) THEN   
        CALL luerrm(11,'(LUFOWO:) no more memory left in LUJETS')   
        h10=-1. 
        h20=-1. 
        h30=-1. 
        h40=-1. 
        RETURN  
      ENDIF 
      np=np+1   
      DO 100 j=1,3  
  100 p(n+np,j)=p(i,j)  
      p(n+np,4)=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2) 
      h0=h0+p(n+np,4)   
      hd=hd+p(n+np,4)**2    
  110 CONTINUE  
      h0=h0**2  
    
C...Very low multiplicities (0 or 1) not considered.    
      IF(np.LE.1) THEN  
        CALL luerrm(8,'(LUFOWO:) too few particles for analysis')   
        h10=-1. 
        h20=-1. 
        h30=-1. 
        h40=-1. 
        RETURN  
      ENDIF 
    
C...Calculate H1 - H4.  
      h10=0.    
      h20=0.    
      h30=0.    
      h40=0.    
      DO 120 i1=n+1,n+np    
      DO 120 i2=i1+1,n+np   
      cthe=(p(i1,1)*p(i2,1)+p(i1,2)*p(i2,2)+p(i1,3)*p(i2,3))/   
     &(p(i1,4)*p(i2,4)) 
      h10=h10+p(i1,4)*p(i2,4)*cthe  
      h20=h20+p(i1,4)*p(i2,4)*(1.5*cthe**2-0.5) 
      h30=h30+p(i1,4)*p(i2,4)*(2.5*cthe**3-1.5*cthe)    
      h40=h40+p(i1,4)*p(i2,4)*(4.375*cthe**4-3.75*cthe**2+0.375)    
  120 CONTINUE  
    
C...Calculate H1/H0 - H4/H0. Output.    
      mstu(61)=n+1  
      mstu(62)=np   
      h10=(hd+2.*h10)/h0    
      h20=(hd+2.*h20)/h0    
      h30=(hd+2.*h30)/h0    
      h40=(hd+2.*h40)/h0    
    
      RETURN    
      END