New functions and constructors added and some other fixes.
[u/mrichter/AliRoot.git] / PYTHIA / jetset / lujmas.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LUJMAS(PMH,PML)
5
6C...Purpose: to determine, approximately, the two jet masses that
7C...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
14C...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
26C...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
50C...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
62C...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
71C...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
84C...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
117C...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
130C...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
146C...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
158C...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