HIJING/hipyset1_35/lusphe_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUSPHE_HIJING(SPH,APL)
   6
   7 C...Purpose: to perform sphericity tensor analysis to give sphericity,
   8 C...aplanarity and the related event axes.
   9 #include "lujets_hijing.inc"
  10 #include "ludat1_hijing.inc"
  11 #include "ludat2_hijing.inc"
  12       DIMENSION SM(3,3),SV(3,3)
  13
  14 C...Calculate matrix to be diagonalized.
  15       NP=0
  16       DO 100 J1=1,3
  17       DO 100 J2=J1,3
  18   100 SM(J1,J2)=0.
  19       PS=0.
  20       DO 120 I=1,N
  21       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 120
  22       IF(MSTU(41).GE.2) THEN
  23         KC=LUCOMP_HIJING(K(I,2))
  24         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  25      &  KC.EQ.18) GOTO 120
  26         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
  27      $       .EQ.0)GOTO 120
  28       ENDIF
  29       NP=NP+1
  30       PA=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  31       PWT=1.
  32       IF(ABS(PARU(41)-2.).GT.0.001) PWT=MAX(1E-10,PA)**(PARU(41)-2.)
  33       DO 110 J1=1,3
  34       DO 110 J2=J1,3
  35   110 SM(J1,J2)=SM(J1,J2)+PWT*P(I,J1)*P(I,J2)
  36       PS=PS+PWT*PA**2
  37   120 CONTINUE
  38
  39 C...Very low multiplicities (0 or 1) not considered.
  40       IF(NP.LE.1) THEN
  41          CALL LUERRM_HIJING(8
  42      $        ,'(LUSPHE_HIJING:) too few particles for analysis')
  43         SPH=-1.
  44         APL=-1.
  45         RETURN
  46       ENDIF
  47       DO 130 J1=1,3
  48       DO 130 J2=J1,3
  49   130 SM(J1,J2)=SM(J1,J2)/PS
  50
  51 C...Find eigenvalues to matrix (third degree equation).
  52       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
  53      &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
  54       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
  55      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
  56       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
  57       P(N+1,4)=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
  58       P(N+3,4)=1./3.+SQRT(-SQ)*MIN(2.*SP,-SQRT(3.*(1.-SP**2))-SP)
  59       P(N+2,4)=1.-P(N+1,4)-P(N+3,4)
  60       IF(P(N+2,4).LT.1E-5) THEN
  61          CALL LUERRM_HIJING(8
  62      $        ,'(LUSPHE_HIJING:) all particles back-to-back')
  63         SPH=-1.
  64         APL=-1.
  65         RETURN
  66       ENDIF
  67
  68 C...Find first and last eigenvector by solving equation system.
  69       DO 170 I=1,3,2
  70       DO 140 J1=1,3
  71       SV(J1,J1)=SM(J1,J1)-P(N+I,4)
  72       DO 140 J2=J1+1,3
  73       SV(J1,J2)=SM(J1,J2)
  74   140 SV(J2,J1)=SM(J1,J2)
  75       SMAX=0.
  76       DO 150 J1=1,3
  77       DO 150 J2=1,3
  78       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 150
  79       JA=J1
  80       JB=J2
  81       SMAX=ABS(SV(J1,J2))
  82   150 CONTINUE
  83       SMAX=0.
  84       DO 160 J3=JA+1,JA+2
  85       J1=J3-3*((J3-1)/3)
  86       RL=SV(J1,JB)/SV(JA,JB)
  87       DO 160 J2=1,3
  88       SV(J1,J2)=SV(J1,J2)-RL*SV(JA,J2)
  89       IF(ABS(SV(J1,J2)).LE.SMAX) GOTO 160
  90       JC=J1
  91       SMAX=ABS(SV(J1,J2))
  92   160 CONTINUE
  93       JB1=JB+1-3*(JB/3)
  94       JB2=JB+2-3*((JB+1)/3)
  95       P(N+I,JB1)=-SV(JC,JB2)
  96       P(N+I,JB2)=SV(JC,JB1)
  97       P(N+I,JB)=-(SV(JA,JB1)*P(N+I,JB1)+SV(JA,JB2)*P(N+I,JB2))/
  98      &SV(JA,JB)
  99       PA=SQRT(P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2)
 100       SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
 101       DO 170 J=1,3
 102   170 P(N+I,J)=SGN*P(N+I,J)/PA
 103
 104 C...Middle axis orthogonal to other two. Fill other codes.
 105       SGN=(-1.)**INT(RLU_HIJING(0)+0.5)
 106       P(N+2,1)=SGN*(P(N+1,2)*P(N+3,3)-P(N+1,3)*P(N+3,2))
 107       P(N+2,2)=SGN*(P(N+1,3)*P(N+3,1)-P(N+1,1)*P(N+3,3))
 108       P(N+2,3)=SGN*(P(N+1,1)*P(N+3,2)-P(N+1,2)*P(N+3,1))
 109       DO 180 I=1,3
 110       K(N+I,1)=31
 111       K(N+I,2)=95
 112       K(N+I,3)=I
 113       K(N+I,4)=0
 114       K(N+I,5)=0
 115       P(N+I,5)=0.
 116       DO 180 J=1,5
 117   180 V(I,J)=0.
 118
 119 C...Select storing option. Calculate sphericity and aplanarity.
 120       MSTU(61)=N+1
 121       MSTU(62)=NP
 122       IF(MSTU(43).LE.1) MSTU(3)=3
 123       IF(MSTU(43).GE.2) N=N+3
 124       SPH=1.5*(P(N+2,4)+P(N+3,4))
 125       APL=1.5*P(N+3,4)
 126
 127       RETURN
 128       END