pyfowo.f

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

 
C*********************************************************************
 
C...PYFOWO
C...Calculates the first few Fox-Wolfram moments.
 
      SUBROUTINE pyfowo(H10,H20,H30,H40)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      SAVE /pyjets/,/pydat1/,/pydat2/
 
C...Copy momenta for particles and calculate H0.
      np=0
      h0=0d0
      hd=0d0
      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=pycomp(k(i,2))
          IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
     &    kc.EQ.18.OR.k(i,2).EQ.ksusy1+22.OR.k(i,2).EQ.39.OR.
     &    k(i,2).EQ.ksusy1+39) goto 110
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0)
     &    goto 110
        ENDIF
        IF(n+np.GE.mstu(4)-mstu(32)-5) THEN
          CALL pyerrm(11,'(PYFOWO:) no more memory left in PYJETS')
          h10=-1d0
          h20=-1d0
          h30=-1d0
          h40=-1d0
          RETURN
        ENDIF
        np=np+1
        DO 100 j=1,3
          p(n+np,j)=p(i,j)
  100   CONTINUE
        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 pyerrm(8,'(PYFOWO:) too few particles for analysis')
        h10=-1d0
        h20=-1d0
        h30=-1d0
        h40=-1d0
        RETURN
      ENDIF
 
C...Calculate H1 - H4.
      h10=0d0
      h20=0d0
      h30=0d0
      h40=0d0
      DO 130 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.5d0*cthe**2-0.5d0)
          h30=h30+p(i1,4)*p(i2,4)*(2.5d0*cthe**3-1.5d0*cthe)
          h40=h40+p(i1,4)*p(i2,4)*(4.375d0*cthe**4-3.75d0*cthe**2+
     &    0.375d0)
  120   CONTINUE
  130 CONTINUE
 
C...Calculate H1/H0 - H4/H0. Output.
      mstu(61)=n+1
      mstu(62)=np
      h10=(hd+2d0*h10)/h0
      h20=(hd+2d0*h20)/h0
      h30=(hd+2d0*h30)/h0
      h40=(hd+2d0*h40)/h0
 
      RETURN
      END