PYTHIA/jetset/lujmas.F

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