hijing/src/lujmas.f

0001
0002 C*********************************************************************
0003
0004       SUBROUTINE LUJMAS(PMH,PML)
0005
0006 C...Purpose: to determine, approximately, the two jet masses that
0007 C...minimize the sum m_H|2 + m_L|2, a la Clavelli and Wyler.
0008       COMMON/LUJETS/N,K(9000,5),P(9000,5),V(9000,5)
0009       SAVE /LUJETS/
0010       COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0011       SAVE /LUDAT1/
0012       COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
0013       SAVE /LUDAT2/
0014       DIMENSION SM(3,3),SAX(3),PS(3,5)
0015
0016 C...Reset.
0017       NP=0
0018       DO 110 J1=1,3
0019       DO 100 J2=J1,3
0020   100 SM(J1,J2)=0.
0021       DO 110 J2=1,4
0022   110 PS(J1,J2)=0.
0023       PSS=0.
0024
0025 C...Take copy of particles that are to be considered in mass analysis.
0026       DO 150 I=1,N
0027       IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
0028       IF(MSTU(41).GE.2) THEN
0029         KC=LUCOMP(K(I,2))
0030         IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
0031      &  KC.EQ.18) GOTO 150
0032         IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE(K(I,2)).EQ.0)
0033      &  GOTO 150
0034       ENDIF
0035       IF(N+NP+1.GE.MSTU(4)-MSTU(32)-5) THEN
0036         CALL LUERRM(11,'(LUJMAS:) no more memory left in LUJETS')
0037         PMH=-2.
0038         PML=-2.
0039         RETURN
0040       ENDIF
0041       NP=NP+1
0042       DO 120 J=1,5
0043   120 P(N+NP,J)=P(I,J)
0044       IF(MSTU(42).EQ.0) P(N+NP,5)=0.
0045       IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
0046       P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
0047
0048 C...Fill information in sphericity tensor and total momentum vector.
0049       DO 130 J1=1,3
0050       DO 130 J2=J1,3
0051   130 SM(J1,J2)=SM(J1,J2)+P(I,J1)*P(I,J2)
0052       PSS=PSS+(P(I,1)**2+P(I,2)**2+P(I,3)**2)
0053       DO 140 J=1,4
0054   140 PS(3,J)=PS(3,J)+P(N+NP,J)
0055   150 CONTINUE
0056
0057 C...Very low multiplicities (0 or 1) not considered.
0058       IF(NP.LE.1) THEN
0059         CALL LUERRM(8,'(LUJMAS:) too few particles for analysis')
0060         PMH=-1.
0061         PML=-1.
0062         RETURN
0063       ENDIF
0064       PARU(61)=SQRT(MAX(0.,PS(3,4)**2-PS(3,1)**2-PS(3,2)**2-PS(3,3)**2))
0065
0066 C...Find largest eigenvalue to matrix (third degree equation).
0067       DO 160 J1=1,3
0068       DO 160 J2=J1,3
0069   160 SM(J1,J2)=SM(J1,J2)/PSS
0070       SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
0071      &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
0072       SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
0073      &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
0074       SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
0075       SMA=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
0076
0077 C...Find largest eigenvector by solving equation system.
0078       DO 170 J1=1,3
0079       SM(J1,J1)=SM(J1,J1)-SMA
0080       DO 170 J2=J1+1,3
0081   170 SM(J2,J1)=SM(J1,J2)
0082       SMAX=0.
0083       DO 180 J1=1,3
0084       DO 180 J2=1,3
0085       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 180
0086       JA=J1
0087       JB=J2
0088       SMAX=ABS(SM(J1,J2))
0089   180 CONTINUE
0090       SMAX=0.
0091       DO 190 J3=JA+1,JA+2
0092       J1=J3-3*((J3-1)/3)
0093       RL=SM(J1,JB)/SM(JA,JB)
0094       DO 190 J2=1,3
0095       SM(J1,J2)=SM(J1,J2)-RL*SM(JA,J2)
0096       IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 190
0097       JC=J1
0098       SMAX=ABS(SM(J1,J2))
0099   190 CONTINUE
0100       JB1=JB+1-3*(JB/3)
0101       JB2=JB+2-3*((JB+1)/3)
0102       SAX(JB1)=-SM(JC,JB2)
0103       SAX(JB2)=SM(JC,JB1)
0104       SAX(JB)=-(SM(JA,JB1)*SAX(JB1)+SM(JA,JB2)*SAX(JB2))/SM(JA,JB)
0105
0106 C...Divide particles into two initial clusters by hemisphere.
0107       DO 200 I=N+1,N+NP
0108       PSAX=P(I,1)*SAX(1)+P(I,2)*SAX(2)+P(I,3)*SAX(3)
0109       IS=1
0110       IF(PSAX.LT.0.) IS=2
0111       K(I,3)=IS
0112       DO 200 J=1,4
0113   200 PS(IS,J)=PS(IS,J)+P(I,J)
0114       PMS=(PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2)+
0115      &(PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2)
0116
0117 C...Reassign one particle at a time; find maximum decrease of m|2 sum.
0118   210 PMD=0.
0119       IM=0
0120       DO 220 J=1,4
0121   220 PS(3,J)=PS(1,J)-PS(2,J)
0122       DO 230 I=N+1,N+NP
0123       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)
0124       IF(K(I,3).EQ.1) PMDI=2.*(P(I,5)**2-PPS)
0125       IF(K(I,3).EQ.2) PMDI=2.*(P(I,5)**2+PPS)
0126       IF(PMDI.LT.PMD) THEN
0127         PMD=PMDI
0128         IM=I
0129       ENDIF
0130   230 CONTINUE
0131
0132 C...Loop back if significant reduction in sum of m|2.
0133       IF(PMD.LT.-PARU(48)*PMS) THEN
0134         PMS=PMS+PMD
0135         IS=K(IM,3)
0136         DO 240 J=1,4
0137         PS(IS,J)=PS(IS,J)-P(IM,J)
0138   240   PS(3-IS,J)=PS(3-IS,J)+P(IM,J)
0139         K(IM,3)=3-IS
0140         GOTO 210
0141       ENDIF
0142
0143 C...Final masses and output.
0144       MSTU(61)=N+1
0145       MSTU(62)=NP
0146       PS(1,5)=SQRT(MAX(0.,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2))
0147       PS(2,5)=SQRT(MAX(0.,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2))
0148       PMH=MAX(PS(1,5),PS(2,5))
0149       PML=MIN(PS(1,5),PS(2,5))
0150
0151       RETURN
0152       END