pysphe.f

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

 
C*********************************************************************
 
C...PYSPHE
C...Performs sphericity tensor analysis to give sphericity,
C...aplanarity and the related event axes.
 
      SUBROUTINE pysphe(SPH,APL)
 
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...Local arrays.
      dimension sm(3,3),sv(3,3)
 
C...Calculate matrix to be diagonalized.
      np=0
      DO 110 j1=1,3
        DO 100 j2=j1,3
          sm(j1,j2)=0d0
  100   CONTINUE
  110 CONTINUE
      ps=0d0
      DO 140 i=1,n
        IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 140
        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 140
          IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.pychge(k(i,2)).EQ.0)
     &    goto 140
        ENDIF
        np=np+1
        pa=sqrt(p(i,1)**2+p(i,2)**2+p(i,3)**2)
        pwt=1d0
        IF(abs(paru(41)-2d0).GT.0.001d0) pwt=
     &  max(1d-10,pa)**(paru(41)-2d0)
        DO 130 j1=1,3
          DO 120 j2=j1,3
            sm(j1,j2)=sm(j1,j2)+pwt*p(i,j1)*p(i,j2)
  120     CONTINUE
  130   CONTINUE
        ps=ps+pwt*pa**2
  140 CONTINUE
 
C...Very low multiplicities (0 or 1) not considered.
      IF(np.LE.1) THEN
        CALL pyerrm(8,'(PYSPHE:) too few particles for analysis')
        sph=-1d0
        apl=-1d0
        RETURN
      ENDIF
      DO 160 j1=1,3
        DO 150 j2=j1,3
          sm(j1,j2)=sm(j1,j2)/ps
  150   CONTINUE
  160 CONTINUE
 
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)/3d0-1d0/9d0
      sr=-0.5d0*(sq+1d0/9d0+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)+1d0/27d0
      sp=cos(acos(max(min(sr/sqrt(-sq**3),1d0),-1d0))/3d0)
      p(n+1,4)=1d0/3d0+sqrt(-sq)*max(2d0*sp,sqrt(3d0*(1d0-sp**2))-sp)
      p(n+3,4)=1d0/3d0+sqrt(-sq)*min(2d0*sp,-sqrt(3d0*(1d0-sp**2))-sp)
      p(n+2,4)=1d0-p(n+1,4)-p(n+3,4)
      IF(p(n+2,4).LT.1d-5) THEN
        CALL pyerrm(8,'(PYSPHE:) all particles back-to-back')
        sph=-1d0
        apl=-1d0
        RETURN
      ENDIF
 
C...Find first and last eigenvector by solving equation system.
      DO 240 i=1,3,2
        DO 180 j1=1,3
          sv(j1,j1)=sm(j1,j1)-p(n+i,4)
          DO 170 j2=j1+1,3
            sv(j1,j2)=sm(j1,j2)
            sv(j2,j1)=sm(j1,j2)
  170     CONTINUE
  180   CONTINUE
        smax=0d0
        DO 200 j1=1,3
          DO 190 j2=1,3
            IF(abs(sv(j1,j2)).LE.smax) goto 190
            ja=j1
            jb=j2
            smax=abs(sv(j1,j2))
  190     CONTINUE
  200   CONTINUE
        smax=0d0
        DO 220 j3=ja+1,ja+2
          j1=j3-3*((j3-1)/3)
          rl=sv(j1,jb)/sv(ja,jb)
          DO 210 j2=1,3
            sv(j1,j2)=sv(j1,j2)-rl*sv(ja,j2)
            IF(abs(sv(j1,j2)).LE.smax) goto 210
            jc=j1
            smax=abs(sv(j1,j2))
  210     CONTINUE
  220   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=(-1d0)**int(pyr(0)+0.5d0)
        DO 230 j=1,3
          p(n+i,j)=sgn*p(n+i,j)/pa
  230   CONTINUE
  240 CONTINUE
 
C...Middle axis orthogonal to other two. Fill other codes.
      sgn=(-1d0)**int(pyr(0)+0.5d0)
      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 260 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)=0d0
        DO 250 j=1,5
          v(i,j)=0d0
  250   CONTINUE
  260 CONTINUE
 
C...Calculate sphericity and aplanarity. Select storing option.
      sph=1.5d0*(p(n+2,4)+p(n+3,4))
      apl=1.5d0*p(n+3,4)
      mstu(61)=n+1
      mstu(62)=np
      IF(mstu(43).LE.1) mstu(3)=3
      IF(mstu(43).GE.2) n=n+3
 
      RETURN
      END