C   11/09/81 807241304  MEMBER NAME  AXSHOW   (S)           FORTRAN
      SUBROUTINE AXSHOW
      IMPLICIT INTEGER*2 (H)
C---
C---     SUBROUTINE TO DISPLAY AXES AND RELATED RESULTS FROM SPHERICITY
C---     AND THRUST ROUTINES. 2,3 AND 4 JET ANALYSIS WITH THE PROGRAM
C---     MCGJET DUE TO M.C. (MONTE CARLO) GODDARD
C---           J.OLSSON 11.09.81      LAST CHANGE   23.06.82
C---       ERROR IN PROJECTIONS OF AXIS DISPLAY CORRECTED
C---       ACCORDING TO K.H. HELLENBRAND             5.11.85
C---    AN ERROR IN THE PRECEEDING CORRECTION CORRECTED
C---       ACCORDING TO K.H. HELLENBRAND            19.11.85
C---     DRAWING OF SPHERICITY AXIS ALSO IN VRX AND VC VIEWS
C---                    J. HAGEMANN                 06.03.86
C---     ALLOW ANALYSIS FOR DIFFERENT PATR-BANKS
C---                    J. HAGEMANN                 13.05.86
C---     DRAW SPHERICITY AXIS THROUGH RUN VERTEX
C---                    J. HAGEMANN                 16.05.86
C-
C-   INPUT FOR SPHERC IS STORED IN COMMON CWORK, ARRAY P
C-  OUTPUT FOR SPHERC IS STORED IN COMMON CWORK
C-
C-   INPUT FOR MCGJET IS STORED IN COMMON SENSE
C-   OUTPUT FOR MCGJET IS STORED IN COMMONS MARKET,COLD,CRUST
C-   NOW IN /CWORK/ FOR SPACE ECONOMY
C-
C-   MCGJET REQUIRES IN THE CASE OF 3-4 JET ANALYSIS ALSO THE INPUT
C-   OF A VECTOR DEFINING THE NORMAL TO THE EVENTPLANE. FOR THIS
C-   PURPOSE, SPHERC IS CALLED AND THE VECTOR CORRESPONDING TO SMALLEST
C-   EIGENVALUE IS USED.
C-
      DIMENSION HCOR(4,100),AXEL(3),AXD(3)
#include "cdata.for"
#include "cgeo1.for"
#include "cgraph.for"
      LOGICAL FL18,FL22,FL24
      COMMON/CPROJ/XMINR,XMAXR,YMINR,YMAXR,IPRJC,FL18,FL22,FL24
      LOGICAL DSPDTM
      COMMON /CGRAP2/ BCMD,DSPDTM(30)
      COMMON /CHEADR/ HEAD(108)
C
      LOGICAL FLVCDO
      COMMON / CGVCDO / FLVCDO(20)
      COMMON /CJTRIG/ PI,TWOPI,PIHALF
      COMMON /CWORK1/ NMOM,P(4,100),T(6),R(9),AXIS(9),SPHITY(3)
      EQUIVALENCE(SPHITY(3),PTHR),(HEAD(18),HRUN)
      COMMON /CWORK/ PAR(3,4),PP(4,100),INUM,DUM(124),
     $ AXTS(3),THR,
     $ IJ(4),PTH(3,100,4),IPJ(100,4)
      COMMON / CVX    / NNPATR,NDUMM,NNJETC,NNVTXC
C     COMMON /SENSE/ PP(4,100),INUM
      DIMENSION IPP(4,100)
      EQUIVALENCE (PP(1,1),IPP(1,1))
C     COMMON /MARKET/ PAR(3,4)
C     COMMON /COLD/ IJ(4),PTH(3,100,4),IPJ(100,4)
C     COMMON /CRUST/ AXTS(3),THR
      DATA ZDEEP /5800./,ICALL/0/
C
      IF(ICALL.NE.0) GO TO 3510
      ZETMX = ZLGPL + ZDEEP
      ICALL = 1
3510  IF( LASTVW .LE. 11  .OR.  LASTVW .EQ. 14  .OR.
     &    LASTVW .EQ. 17  .OR.  LASTVW .EQ. 20 ) GO TO 3511
      CALL TRMOUT(80,' NOT AVAILABLE IN THIS VIEW ...^')
      GO TO 2235
3511  NJ=ACMD
      NJ1=BCMD
      IF(NJ.NE.0.OR.FL18) GO TO 100
      CALL TRMOUT(80,'CODE 1   : DISPLAY SPHERICITY AXIS^')
      CALL TRMOUT(80,'CODE 1 -1: DISPLAY SPHERICITY AXIS,PRINT Q VALUES^
     $')
      CALL TRMOUT(80,'CODE 2   : 2-JET ANALYSIS; AXIS DISPLAY^')
      CALL TRMOUT(80,'CODE 3   : 3-JET ANALYSIS; AXIS DISPLAY^')
      CALL TRMOUT(80,'CODE 4   : 4-JET ANALYSIS; AXIS DISPLAY^')
      CALL TRMOUT(80,'CODE K L : K-JET ANALYSIS; DISPLAY OF TRACKS BELON
     $GING TO JET L^')
      CALL TRMOUT(80,'CODE K -1: K-JET ANALYSIS; DISPLAY OF THRUST AXIS^
     $')
      CALL TRMOUT(80,'ENTER OPTION:^')
      CALL FYRINT(NJ,NJ1,IDM3,IDM4)
100   IF(NJ.GE.1.AND.NJ.LE.4) GO TO 200
      GO TO 2531
200   IF(NJ1.GE.-1.AND.NJ1.LE.NJ) GO TO 201
      GO TO 2531
201   CONTINUE
C
C                   ---------- ***   SET UP INPUT VECTORS *** --------
C
      EBM = EBEAM(HRUN)*.001
C
      IF( NNPATR .LT. 1 .OR. NNPATR .GT. 99 ) GO TO 2220
         CALL CLOC( IPO, 'PATR', NNPATR )
         IF( IPO .GT. 0 ) GO TO 2224
            WRITE(6,2219) NNPATR
2219        FORMAT(' PATR BANK NR',I3,' (GIVEN BY COMMAND PATR) NOT EXIS
     &TING')
2220  IPO=IDATA(IBLN('PATR'))
      IF(IPO.NE.0) GO TO 2224
      CALL TRMOUT(80,'PATR BANK DOES NOT EXIST^')
      GO TO 2235
2224  IPLGCL = IDATA(IBLN('LGCL'))
      IF(IPLGCL.NE.0) GO TO 2223
      CALL TRMOUT(80,'LGCL BANK DOES NOT EXIST^')
      GO TO 2235
2223  IPALGN=IDATA(IBLN('ALGN'))
      IF(IDATA(IPLGCL+21).NE.2) CALL LGCDIR(IPO,IPALGN,IPLGCL)
      NWPCL = IDATA(IPLGCL+25)
      IPOG = IDATA(IPLGCL+3) + IPLGCL - 1 - NWPCL
      NCLST = IDATA(IPLGCL+7)
      LO = IDATA(IPO+1)
      NTR = IDATA(IPO+2)
      NMOM = 0
      INUM = 0
      IF(NTR.EQ.0) GO TO 523
      LTR = IDATA(IPO+3)
      IPZ = IPO + LO - LTR
      ITP = 1
      DO 1004  JTR = 1,NTR
      IPZ = IPZ + LTR
      CALL TRACK(IPZ,ITP,RMIN,NHT,FI,THETA,PTOT,PT,PL,PX,PY,PZ)
      IF(PTOT.LT.EBM) GO TO 3910
      WRITE(6,3911) JTR,PTOT,EBM
3911  FORMAT(' TRACK NR ',I2,' SCALED FROM ',E12.4,' TO ',E12.4)
      PTOTR = PTOT
      PTOT = EBM
      FACTOR = EBM/PTOTR
      PT = PT*FACTOR
      PL = PL*FACTOR
      PX = PX*FACTOR
      PY = PY*FACTOR
      PZ = PZ*FACTOR
3910  NMOM = NMOM + 1
      INUM = INUM + 1
      HCOR(1,NMOM) = ITP
      HCOR(2,NMOM) = IPZ
      HCOR(3,NMOM) = JTR
      HCOR(4,NMOM) = 0
      P(1,NMOM) = PX
      P(2,NMOM) = PY
      P(3,NMOM) = PZ
      P(4,NMOM) = PTOT
      PP(1,NMOM) = PX
      PP(2,NMOM) = PY
      PP(3,NMOM) = PZ
      IPP(4,NMOM) = NMOM
1004  CONTINUE
C  PHOTON PART
523   IF(NCLST.EQ.0) GO TO 1015
      IPZ = IPOG
      IGAM = 0
      ITP = 0
      DO 522  IGM = 1,NCLST
      IPZ = IPZ + NWPCL
      CALL TRACK(IPZ,ITP,RMIN,NHT,FI,THETA,PTOT,PT,PL,PX,PY,PZ)
      IF(IDATA(IPZ+8).EQ.0) GO TO 8800
      IF(PTOT.LT.1.0) GO TO 522
C SUBTRACK STANDARD 400 MEV FROM TRACK ASSOCIATED CLUSTERS ABOVE 1 GEV
      PTOT1 = PTOT - .400
      FACTOR = PTOT1/PTOT
      PX = PX*FACTOR
      PY = PY*FACTOR
      PZ = PZ*FACTOR
      PTOT = PTOT1
      IGAM = -1
      GO TO 8801
8800  IGAM = IGAM + 1
8801  CALL MOMGAM(IPZ,PX,PY,PZ,PT,PTOT,FI,THE)
      NMOM = NMOM + 1
      INUM = INUM + 1
      HCOR(1,NMOM) = ITP
      HCOR(2,NMOM) = IPZ
      HCOR(3,NMOM) = IGAM
      HCOR(4,NMOM) = 0
      P(1,NMOM) = PX
      P(2,NMOM) = PY
      P(3,NMOM) = PZ
      P(4,NMOM) = PTOT
      PP(1,NMOM) = PX
      PP(2,NMOM) = PY
      PP(3,NMOM) = PZ
      IPP(4,NMOM) = NMOM
522   CONTINUE
1015  IF(NMOM.GT.0) GO TO 1017
      CALL TRMOUT(80,'NO TRACKS, NO PHOTONS !!!^')
      GO TO 2235
1017  CALL SPHERC
      IF(NJ.GT.1) GO TO 2000
C
C--------------------------- *** SPHERICITY DISPLAY *** -------------
C
      SUM = SPHITY(1)+SPHITY(2)+SPHITY(3)
      Q1 = 1. - 2.*SPHITY(3)/SUM
      Q2 = 1. - 2.*SPHITY(2)/SUM
      Q3 = 1. - 2.*SPHITY(1)/SUM
      Q4 = (Q3-Q2)/SQRT(3.)
      Q5 = Q2 - Q1
C  WRITE OUTPUT OF SPHERC
      IF(.NOT.FL18) WRITE(6,7892) AXIS(1),AXIS(2),AXIS(3),SPHITY(1)
7892  FORMAT(' SPHERICITY:  AXES ',3F7.3,'   VALUE  ',F7.4)
      IF(NJ1.EQ.-1.AND..NOT.FL18) WRITE(6,7895) Q1,Q2,Q3
7895  FORMAT(' Q1-3 :  ',3E12.4)
      IF(NJ1.EQ.-1.AND..NOT.FL18) WRITE(6,7897) Q4,Q5
7897  FORMAT(' (Q3-Q2)/ROOT3 : ',E12.4,'    Q2-Q1 : ',E12.4)
C--
      DO 7575  I = 1,3
7575  AXD(I) = AXIS(I)
      IRETUR = 1
      GO TO 7000
C
C ------------ ***  JET AXIS ANALYSIS  WITH MCGJET *** ----------
C
2000  CONTINUE
      DO 2001  I = 1,3
2001  AXEL(I) = AXIS(6+I)
      CALL MCGJET(NJ,AXEL)
      IF(NJ1.NE.0) GO TO 3000
C
C ------------ ***  JET AXIS DISPLAY  *** ----------
C
      DO 2002  NJET = 1,NJ
      DO 7576  I = 1,3
7576  AXD(I) = PAR(I,NJET)
      IRETUR = 2
      GO TO 7000
2003  CONTINUE
2002  CONTINUE
      GO TO 2235
C
C ------------ ***  THRUST AXIS DISPLAY  *** ----------
C
3000  IF(NJ1.NE.-1) GO TO 4000
      DO 7577  I = 1,3
7577  AXD(I) = AXTS(I)
      IRETUR = 3
      GO TO 7000
3001  IF(.NOT.FL18) WRITE(6,3002) AXTS(1),AXTS(2),AXTS(3),THR
3002  FORMAT(' THRUST:  AXES ',3F7.3,'   VALUE  ',F7.4)
      GO TO 2235
C
C ------------ ***  DISPLAY TRACKS OF SINGLE JETS *** ---------
C
4000  RTS = RLG
      IF(LASTVW.EQ.1.OR.LASTVW.EQ.4.OR.LASTVW.EQ.8) RTS = RTOF
      VX = 0.
      VY = 0.
      VZ = 0.
      IJNRT = IJ(NJ1)
      DO 4001  IJK = 1,IJNRT
      IMOM = IPJ(IJK,NJ1)
      IF(HCOR(1,IMOM).EQ.0) GO TO 4010
C           CHARGED TRACKS
      IPO = HCOR(2,IMOM)
      ITR = HCOR(3,IMOM)
      ITYPLI = IDATA(IPO+29)
      INDEX = LASTVW
      IF(LASTVW.GT.3.AND.LASTVW.NE.14) GO TO 4002
C** DISPLAY TRACKS, RFI VIEWS
C                     PARABOLA FITS
      IF(IDATA(IPO+18).NE.1.AND.ITYPLI.EQ.1)
     $ CALL PARDIS(25,ADATA(IPO+19),ADATA(IPO+20),ADATA(IPO+21),
     $ ADATA(IPO+22),ADATA(IPO+5),ADATA(IPO+6),ADATA(IPO+12),
     $ ADATA(IPO+13),ADATA(IPO+31),ADATA(IPO+30),0.,0.)
C                     CIRCLE FITS
      IF(IDATA(IPO+18).EQ.1.AND.(ITYPLI.EQ.1.OR.ITYPLI.EQ.2))
     $ CALL CIRDIS(25,SIGN(ADATA(IPO+19),ADATA(IPO+25)),
     $ ADATA(IPO+20),ADATA(IPO+21),
     $ ADATA(IPO+5),ADATA(IPO+6),ADATA(IPO+12),ADATA(IPO+13),
     $ ADATA(IPO+31),ADATA(IPO+30),ITYPLI,0.,0.)
C** WRITE TRACK NUMBER
      CALL TRNUMB(ITR,0,ADATA(IPO+12),ADATA(IPO+13),ADATA(IPO+14))
      GO TO 4001
C** DISPLAY TRACKS, RZ VIEWS
4002  IF(ITYPLI.EQ.1) CALL LINDIS(IPO,INDEX,0.,0.,0.)
      GO TO 4001
4010  CONTINUE
C           NEUTRAL TRACKS
      IPO = HCOR(2,IMOM)
      IGAM = HCOR(3,IMOM)
C SKIP TRACK ASSOCIATED CLUSTERS
      IF(IGAM.EQ.-1) GO TO 4001
      JPART = IDATA(IPO+1)
      IF(JPART.NE.0) GO TO 4011
C****************************** BARREL PHOTONS *****
      PHI = ADATA(IPO+4)
      THE = ATAN2(ADATA(IPO+5)-VZ,RLG)
      THE = PIHALF - THE
      XP = RTS*COS(PHI)
      YP = RTS*SIN(PHI)
      IF(ABS(PIHALF-THE).LT..001) ZP = VZ
      IF(ABS(PIHALF-THE).GE..001) ZP = RTS/TAN(THE) + VZ
      VXP = VX
      VYP = VY
      VZP = VZ
      IF(LASTVW.NE.14) GO TO 4012
      FP = (ZDEEP - ZP)/(ZDEEP + ZLGPL)
      FP1 = (ZDEEP - VZP)/(ZDEEP + ZLGPL)
      XP = FP*XP
      YP = FP*YP
      VXP = FP1*VXP
      VYP = FP1*VYP
      GO TO 4012
C****************************** ENDCAP PHOTONS *****
4011  XP = ADATA(IPO+4)
      YP = ADATA(IPO+5)
      RR = SQRT(XP*XP + YP*YP)
      ZZ = ZENDPL
      IF(JPART.LT.0) ZZ = ZENDMI
      ZZ = ZZ - VZ
      THE = ATAN2(RR,ZZ)
      PHI = ATAN2(YP,XP)
      IF(PHI.LT.0.) PHI = PHI + TWOPI
      ZP = ZENDPL
      IF(JPART.LT.0) ZP = ZENDMI
      VXP = VX
      VYP = VY
      VZP = VZ
      IF(LASTVW.NE.14) GO TO 4012
      FP = (ZDEEP - ZP)/(ZDEEP + ZLGPL)
      FP1 = (ZDEEP - VZP)/(ZDEEP + ZLGPL)
      XP = FP*XP
      YP = FP*YP
      VXP = FP1*VXP
      VYP = FP1*VYP
4012  CONTINUE
C * * * * * * * * *  DISPLAY PART
      IF(LASTVW.GT.3.AND.LASTVW.NE.14) GO TO 4013
      XX = XP
      YY = YP
      XX1 = VXP
      YY1 = VYP
      CALL TRNUMB(IGAM,0,XX,YY,DUM1)
      XX = - XX
      XX1 = - XX1
      GO TO 4014
4013  CONTINUE
      XX = ZP
      XX1 = VZP
      IF(DSPDTL(9)) GO TO 4015
      YY = XP
      YY1 = VXP
      IF(LASTVW.GT.7) YY = YP
      IF(LASTVW.GT.7) YY1 = VYP
      GO TO 4016
4015  YY = SQRT(XP*XP+YP*YP)
      YY1 = SQRT(VXP*VXP+VYP*VYP)
      IF(YP.LT.0..AND.LASTVW.GT.7) YY = - YY
      IF(XP.LT.0..AND.LASTVW.LT.8) YY = - YY
      IF(VYP.LT.0..AND.LASTVW.GT.7) YY1 = - YY1
      IF(VXP.LT.0..AND.LASTVW.LT.8) YY1 = - YY1
4016  CALL TRNUMB(IGAM,0,-XX,YY,DUM1)
4014  CALL MOVEA(XX1,YY1)
      CALL DASHA(XX,YY,14)
4001  CONTINUE
      GO TO 2235
C---------------------------------------------------------------------
C
C    PROCEDURE FOR AXIS DISPLAY
C
7000  CONTINUE
      SINTHE = SQRT(1.-AXD(3)*AXD(3))
      RTF1 = RTOF - 50.
      IF( LASTVW .EQ. 17 ) RTF1 = 5.5
      IF( LASTVW .EQ. 20 ) RTF1 = RITNK + 110.0
      XX1 = RTF1*AXD(1)
      YY1 = RTF1*AXD(2)
      ZZ1 = RTF1*AXD(3)
      RR1 = 0.
      IF(SINTHE.GT..00004) RR1 = SQRT(XX1**2 + YY1**2)
      IF( LASTVW .GT. 3   .AND.  LASTVW .NE. 14  .AND.
     &    LASTVW .NE. 17  .AND.  LASTVW .NE. 20 ) GO TO 7001
C RFI VIEWS
      XRV = 0.0
      YRV = 0.0
      IF( .NOT. FLVCDO(12) ) GO TO 6997
         CALL VTXCRV( XRV, YRV, DXR, DYR )
6997  CONTINUE
      DLX = XX1
      DLY = YY1
      XX1 = XRV + DLX
      YY1 = YRV + DLY
      XX2 = XRV - DLX
      YY2 = YRV - DLY
      IF( SINTHE .LE. 0.00004 ) GO TO 6998
         DLX = DLX/SINTHE
         DLY = DLY/SINTHE
         XX1 = XRV + DLX
         XX2 = XRV - DLX
         YY1 = YRV + DLY
         YY2 = YRV - DLY
6998  CONTINUE
      IF( LASTVW .NE. 14 ) GO TO 6999
C RFI CYLINDER VIEW
         XX1 = XX1*(ZDEEP - ZZ1)/ZETMX
         YY1 = YY1*(ZDEEP - ZZ1)/ZETMX
         XX2 = XX2*(ZDEEP - ZZ1)/ZETMX
         YY2 = YY2*(ZDEEP - ZZ1)/ZETMX
6999  CALL MOVEA( -XRV,  YRV )
      WARRO = 20.0
      IF( LASTVW .EQ. 17 ) WARRO = 0.2
      IF( LASTVW .EQ. 20 ) WARRO = 10.0
      CALL DASHA( -XX1, YY1, 78 )
      CALL MOVEA( -XRV,  YRV )
      CALL DASHA( -XX1, YY1, 18 )
      PHIAR = ATAN2(DLY,-DLX)
      CALL ARROW( WARRO, -XX1, YY1,  PHIAR )
      IF(IRETUR.NE.1.AND.IRETUR.NE.3) GO TO 2237
      CALL MOVEA( -XRV,  YRV )
      CALL DASHA( -XX2, YY2, 78 )
      CALL MOVEA( -XRV,  YRV )
      CALL DASHA( -XX2, YY2, 18 )
      PHIAR = ATAN2(-DLY,DLX)
      CALL ARROW( WARRO, -XX2, YY2,  PHIAR )
      GO TO 2237
C ZX AND ZY VIEWS
7001  IF(LASTVW.GT.11) GO TO 2235
      IF(DSPDTL(9)) GO TO 7002
C TRUE PROJECTIONS
      XXX1 = ZZ1
      YYY1 = XX1
      IF(LASTVW.GT.7) YYY1 = YY1
      GO TO 7003
C ROTATED PROJECTIONS
7002  XXX1 = ZZ1
      YYY1 = RR1
      IF(LASTVW.GT.7.AND.YY1.LT.0.) YYY1 = -RR1
      IF(LASTVW.LT.8.AND.XX1.LT.0.) YYY1 = -RR1
7003  CALL MOVEA(0.,0.)
      CALL DASHA(XXX1,YYY1,78)
      CALL MOVEA(0.,0.)
      CALL DASHA(XXX1,YYY1,18)
      CALL ARROW(20.,XXX1,YYY1,-1000.0)
      IF(IRETUR.NE.1.AND.IRETUR.NE.3) GO TO 2237
      CALL MOVEA(0.,0.)
      CALL DASHA(-XXX1,-YYY1,78)
      CALL MOVEA(0.,0.)
      CALL DASHA(-XXX1,-YYY1,18)
      CALL ARROW(20.,-XXX1,-YYY1,-1000.0)
2237  CONTINUE
      GO TO (2235,2003,3001),IRETUR
2531  CALL TRMOUT(80,' WRONG INPUT...^')
2235  RETURN
      END
      SUBROUTINE ARROW(R,ADX,ADY,PHIAR)
      COMMON /CJTRIG/ PI,TWOPI
C     DRAW TRIANGLE
      DEFI = TWOPI/3.
      OFFSET = PHIAR
      IF( PHIAR .EQ. -1000.0 ) OFFSET = ATAN2(ADY,ADX)
      IF(OFFSET.LT.0.) OFFSET = OFFSET + TWOPI
      FISX = OFFSET
      CALL MOVEA(ADX+R*COS(OFFSET),ADY+R*SIN(OFFSET))
      DO 27 IIK = 1,3
      FISX = FISX + DEFI
 27   CALL DRAWA(ADX+R*COS(FISX),ADY+R*SIN(FISX))
      RETURN
      END
      SUBROUTINE SPHERC
C *---------------------------------------------------------
C *  VERSION OF SPHRCY, MODIFIED FOR GRAPHICS USE 21.07.79
C *                                  LAST CHANGE  21.07.79  (J.OLSSON)
C *  VERSION OF 23/02/79
C *  COMPUTE SPHERICITY FOR A SET OF MOMENTA
C *  P IS ARRAY OF MOMENTA P(1..4,1..NMOM), WHERE P(1..3)=PX,PY,PZ
C *  AND P(4)=PTOT
C *  SPHITY IS ARRAY OF NORMALIZED SPHERICITIES
C *  AXIS CONTAINS THE THREE EIGENVECTORS
C *  AXIS(1..3)  EIGENVECTOR WITH SMALLEST EIGENVALUE
C *  AXIS(4..6)  EIGENVECTOR
C *  AXIS(7..9)  EIGENVECTOR WITH BIGGEST EIGENVALUE
C *---------------------------------------------------------
C
      COMMON /CWORK1/ NMOM,P(4,100),T(6),R(9),AXIS(9),SPHITY(3)
C
C  LOAD MOMENTUM TENSOR
C
       DO 10 I=1,6
   10  T(I)=0.
       DO 20 J=1,NMOM
       T(1) = T(1) + P(2,J)*P(2,J) + P(3,J)*P(3,J)
       T(2) = T(2) - P(1,J)*P(2,J)
       T(3) = T(3) + P(1,J)*P(1,J) + P(3,J)*P(3,J)
       T(4) = T(4) - P(1,J)*P(3,J)
       T(5) = T(5) - P(2,J)*P(3,J)
   20  T(6) = T(6) + P(1,J)*P(1,J) + P(2,J)*P(2,J)
CCCCCCCCC
C   DIAGONALIZE MOMENTUM TENSOR T
C
       CALL EIGEN( T, R, 3, 0 )
C
C      LOAD SPHERICITY AND JET AXIS
C
       SUM = T(1) + T(3) + T(6)
       SPHITY(1) = 3.* T(6) / SUM
       SPHITY(2) = 3.* T(3) / SUM
       SPHITY(3) = 3.* T(1) / SUM
       DO 30 I=1,3
   30  AXIS(I) = R(I+6)
       DO 40 I=4,6
   40  AXIS(I)= R(I)
       DO 50 I=7,9
   50  AXIS(I)= R(I-6)
      RETURN
      END
