3 C*********************************************************************
5 SUBROUTINE LUJMAS_HIJING(PMH,PML)
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 #include "lujets_hijing.inc"
10 #include "ludat1_hijing.inc"
11 #include "ludat2_hijing.inc"
12 DIMENSION SM(3,3),SAX(3),PS(3,5)
23 C...Take copy of particles that are to be considered in mass analysis.
25 IF(K(I,1).LE.0.OR.K(I,1).GT.10) GOTO 150
26 IF(MSTU(41).GE.2) THEN
27 KC=LUCOMP_HIJING(K(I,2))
28 IF(KC.EQ.0.OR.KC.EQ.12.OR.KC.EQ.14.OR.KC.EQ.16.OR.
30 IF(MSTU(41).GE.3.AND.KCHG(KC,2).EQ.0.AND.LUCHGE_HIJING(K(I,2))
33 IF(N+NP+1.GE.MSTU(4)-MSTU(32)-5) THEN
35 $ ,'(LUJMAS_HIJING:) no more memory left in LUJETS_HIJING')
43 IF(MSTU(42).EQ.0) P(N+NP,5)=0.
44 IF(MSTU(42).EQ.1.AND.K(I,2).NE.22) P(N+NP,5)=PMAS(101,1)
45 P(N+NP,4)=SQRT(P(N+NP,5)**2+P(I,1)**2+P(I,2)**2+P(I,3)**2)
47 C...Fill information in sphericity tensor and total momentum vector.
50 130 SM(J1,J2)=SM(J1,J2)+P(I,J1)*P(I,J2)
51 PSS=PSS+(P(I,1)**2+P(I,2)**2+P(I,3)**2)
53 140 PS(3,J)=PS(3,J)+P(N+NP,J)
56 C...Very low multiplicities (0 or 1) not considered.
59 $ ,'(LUJMAS_HIJING:) too few particles for analysis')
64 PARU(61)=SQRT(MAX(0.,PS(3,4)**2-PS(3,1)**2-PS(3,2)**2-PS(3,3)**2))
66 C...Find largest eigenvalue to matrix (third degree equation).
69 160 SM(J1,J2)=SM(J1,J2)/PSS
70 SQ=(SM(1,1)*SM(2,2)+SM(1,1)*SM(3,3)+SM(2,2)*SM(3,3)-SM(1,2)**2-
71 &SM(1,3)**2-SM(2,3)**2)/3.-1./9.
72 SR=-0.5*(SQ+1./9.+SM(1,1)*SM(2,3)**2+SM(2,2)*SM(1,3)**2+SM(3,3)*
73 &SM(1,2)**2-SM(1,1)*SM(2,2)*SM(3,3))+SM(1,2)*SM(1,3)*SM(2,3)+1./27.
74 SP=COS(ACOS(MAX(MIN(SR/SQRT(-SQ**3),1.),-1.))/3.)
75 SMA=1./3.+SQRT(-SQ)*MAX(2.*SP,SQRT(3.*(1.-SP**2))-SP)
77 C...Find largest eigenvector by solving equation system.
79 SM(J1,J1)=SM(J1,J1)-SMA
81 170 SM(J2,J1)=SM(J1,J2)
85 IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 180
93 RL=SM(J1,JB)/SM(JA,JB)
95 SM(J1,J2)=SM(J1,J2)-RL*SM(JA,J2)
96 IF(ABS(SM(J1,J2)).LE.SMAX) GOTO 190
101 JB2=JB+2-3*((JB+1)/3)
104 SAX(JB)=-(SM(JA,JB1)*SAX(JB1)+SM(JA,JB2)*SAX(JB2))/SM(JA,JB)
106 C...Divide particles into two initial clusters by hemisphere.
108 PSAX=P(I,1)*SAX(1)+P(I,2)*SAX(2)+P(I,3)*SAX(3)
113 200 PS(IS,J)=PS(IS,J)+P(I,J)
114 PMS=(PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2)+
115 &(PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2)
117 C...Reassign one particle at a time; find maximum decrease of m|2 sum.
121 220 PS(3,J)=PS(1,J)-PS(2,J)
123 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)
124 IF(K(I,3).EQ.1) PMDI=2.*(P(I,5)**2-PPS)
125 IF(K(I,3).EQ.2) PMDI=2.*(P(I,5)**2+PPS)
132 C...Loop back if significant reduction in sum of m|2.
133 IF(PMD.LT.-PARU(48)*PMS) THEN
137 PS(IS,J)=PS(IS,J)-P(IM,J)
138 240 PS(3-IS,J)=PS(3-IS,J)+P(IM,J)
143 C...Final masses and output.
146 PS(1,5)=SQRT(MAX(0.,PS(1,4)**2-PS(1,1)**2-PS(1,2)**2-PS(1,3)**2))
147 PS(2,5)=SQRT(MAX(0.,PS(2,4)**2-PS(2,1)**2-PS(2,2)**2-PS(2,3)**2))
148 PMH=MAX(PS(1,5),PS(2,5))
149 PML=MIN(PS(1,5),PS(2,5))