PHOS/shaker/lujmas.f

   1 *CMZ :          17/07/98  15.44.34  by  Federico Carminati
   2 *-- Author :
   3 C*********************************************************************
   4
   5       SUBROUTINE LUJMAS(PMH,PML)
   6
   7 C...Purpose: to determine, approximately, the two jet masses that
   8 C...minimize the sum m_H^2 + m_L^2, a la Clavelli and Wyler.
   9 *KEEP,LUJETS.
  10       COMMON /LUJETS/ N,K(200000,5),P(200000,5),V(200000,5)
  11       SAVE /LUJETS/
  12 *KEEP,LUDAT1.
  13       COMMON /LUDAT1/ MSTU(200),PARU(200),MSTJ(200),PARJ(200)
  14       SAVE /LUDAT1/
  15 *KEEP,LUDAT2.
  16       COMMON /LUDAT2/ KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
  17       SAVE /LUDAT2/
  18 *KEND.
  19       DIMENSION SM(3,3),SAX(3),PS(3,5)
  20
  21 C...Reset.
  22       NP=0
  23       DO 110 J1=1,3
  24       DO 100 J2=J1,3
  25   100 SM(J1,J2)=0.
  26       DO 110 J2=1,4
  27   110 PS(J1,J2)=0.
  28       PSS=0.
  29
  30 C...Take copy of particles that are to be considered in mass analysis.
  31       DO 150 I=1,N
  32       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
  33       IF(MSTU(41).GE.2) THEN
  34         KC=LUCOMP(K(I,2))
  35         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
  36      &  KC.EQ.18) GOTO 150
  37         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
  38      &  GOTO 150
  39       ENDIF
  40       IF(N+NP+1.GE.MSTU(4)-MSTU(32)-5) THEN
  41         CALL LUERRM(11,'(LUJMAS:) no more memory left in LUJETS')
  42         PMH=-2.
  43         PML=-2.
  44         RETURN
  45       ENDIF
  46       NP=NP+1
  47       DO 120 J=1,5
  48   120 P(N+NP,J)=P(I,J)
  49       IF(MSTU(42).EQ.0) P(N+NP,5)=0.
  50       IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
  51       P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
  52
  53 C...Fill information in sphericity tensor and total momentum vector.
  54       DO 130 J1=1,3
  55       DO 130 J2=J1,3
  56   130 SM(J1,J2)=SM(J1,J2)+P(I,J1)*P(I,J2)
  57       PSS=PSS+(P(I,1)**2+P(I,2)**2+P(I,3)**2)
  58       DO 140 J=1,4
  59   140 PS(3,J)=PS(3,J)+P(N+NP,J)
  60   150 CONTINUE
  61
  62 C...Very low multiplicities (0 or 1) not considered.
  63       IF(NP.LE.1) THEN
  64         CALL LUERRM(8,'(LUJMAS:) too few particles for analysis')
  65         PMH=-1.
  66         PML=-1.
  67         RETURN
  68       ENDIF
  69       PARU(61)=SQRT(MAX(0.,PS(3,4)**2-PS(3,1)**2-PS(3,2)**2-PS(3,3)**2))
  70
  71 C...Find largest eigenvalue to matrix (third degree equation).
  72       DO 160 J1=1,3
  73       DO 160 J2=J1,3
  74   160 SM(J1,J2)=SM(J1,J2)/PSS
  75       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
  76      &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
  77       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
  78      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
  79       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
  80       SMA=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
  81
  82 C...Find largest eigenvector by solving equation system.
  83       DO 170 J1=1,3
  84       SM(J1,J1)=SM(J1,J1)-SMA
  85       DO 170 J2=J1+1,3
  86   170 SM(J2,J1)=SM(J1,J2)
  87       SMAX=0.
  88       DO 180 J1=1,3
  89       DO 180 J2=1,3
  90       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 180
  91       JA=J1
  92       JB=J2
  93       SMAX=ABS(SM(J1,J2))
  94   180 CONTINUE
  95       SMAX=0.
  96       DO 190 J3=JA+1,JA+2
  97       J1=J3-3*((J3-1)/3)
  98       RL=SM(J1,JB)/SM(JA,JB)
  99       DO 190 J2=1,3
 100       SM(J1,J2)=SM(J1,J2)-RL*SM(JA,J2)
 101       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 190
 102       JC=J1
 103       SMAX=ABS(SM(J1,J2))
 104   190 CONTINUE
 105       JB1=JB+1-3*(JB/3)
 106       JB2=JB+2-3*((JB+1)/3)
 107       SAX(JB1)=-SM(JC,JB2)
 108       SAX(JB2)=SM(JC,JB1)
 109       SAX(JB)=-(SM(JA,JB1)*SAX(JB1)+SM(JA,JB2)*SAX(JB2))/SM(JA,JB)
 110
 111 C...Divide particles into two initial clusters by hemisphere.
 112       DO 200 I=N+1,N+NP
 113       PSAX=P(I,1)*SAX(1)+P(I,2)*SAX(2)+P(I,3)*SAX(3)
 114       IS=1
 115       IF(PSAX.LT.0.) IS=2
 116       K(I,3)=IS
 117       DO 200 J=1,4
 118   200 PS(IS,J)=PS(IS,J)+P(I,J)
 119       PMS=MAX(1E-10,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2)+
 120      &MAX(1E-10,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2)
 121
 122 C...Reassign one particle at a time; find maximum decrease of m^2 sum.
 123   210 PMD=0.
 124       IM=0
 125       DO 220 J=1,4
 126   220 PS(3,J)=PS(1,J)-PS(2,J)
 127       DO 230 I=N+1,N+NP
 128       PPS=P(I,4)*PS(3,4)-P(I,1)*PS(3,1)-P(I,2)*PS(3,2)-P(I,3)*PS(3,3)
 129       IF(K(I,3).EQ.1) PMDI=2.*(P(I,5)**2-PPS)
 130       IF(K(I,3).EQ.2) PMDI=2.*(P(I,5)**2+PPS)
 131       IF(PMDI.LT.PMD) THEN
 132         PMD=PMDI
 133         IM=I
 134       ENDIF
 135   230 CONTINUE
 136
 137 C...Loop back if significant reduction in sum of m^2.
 138       IF(PMD.LT.-PARU(48)*PMS) THEN
 139         PMS=PMS+PMD
 140         IS=K(IM,3)
 141         DO 240 J=1,4
 142         PS(IS,J)=PS(IS,J)-P(IM,J)
 143   240   PS(3-IS,J)=PS(3-IS,J)+P(IM,J)
 144         K(IM,3)=3-IS
 145         GOTO 210
 146       ENDIF
 147
 148 C...Final masses and output.
 149       MSTU(61)=N+1
 150       MSTU(62)=NP
 151       PS(1,5)=SQRT(MAX(0.,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2))
 152       PS(2,5)=SQRT(MAX(0.,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2))
 153       PMH=MAX(PS(1,5),PS(2,5))
 154       PML=MIN(PS(1,5),PS(2,5))
 155
 156       RETURN
 157       END