Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lujmas.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14  dimension sm(3,3),sax(3),ps(3,5)
15 
16 C...Reset.
17  np=0
18  DO 110 j1=1,3
19  DO 100 j2=j1,3
20  100 sm(j1,j2)=0.
21  DO 110 j2=1,4
22  110 ps(j1,j2)=0.
23  pss=0.
24 
25 C...Take copy of particles that are to be considered in mass analysis.
26  DO 150 i=1,n
27  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 150
28  IF(mstu(41).GE.2) THEN
29  kc=lucomp(k(i,2))
30  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
31  & kc.EQ.18) goto 150
32  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
33  & goto 150
34  ENDIF
35  IF(n+np+1.GE.mstu(4)-mstu(32)-5) THEN
36  CALL luerrm(11,'(LUJMAS:) no more memory left in LUJETS')
37  pmh=-2.
38  pml=-2.
39  RETURN
40  ENDIF
41  np=np+1
42  DO 120 j=1,5
43  120 p(n+np,j)=p(i,j)
44  IF(mstu(42).EQ.0) p(n+np,5)=0.
45  IF(mstu(42).EQ.1.AND.k(i,2).NE.22) p(n+np,5)=pmas(101,1)
46  p(n+np,4)=sqrt(p(n+np,5)**2+p(i,1)**2+p(i,2)**2+p(i,3)**2)
47 
48 C...Fill information in sphericity tensor and total momentum vector.
49  DO 130 j1=1,3
50  DO 130 j2=j1,3
51  130 sm(j1,j2)=sm(j1,j2)+p(i,j1)*p(i,j2)
52  pss=pss+(p(i,1)**2+p(i,2)**2+p(i,3)**2)
53  DO 140 j=1,4
54  140 ps(3,j)=ps(3,j)+p(n+np,j)
55  150 CONTINUE
56 
57 C...Very low multiplicities (0 or 1) not considered.
58  IF(np.LE.1) THEN
59  CALL luerrm(8,'(LUJMAS:) too few particles for analysis')
60  pmh=-1.
61  pml=-1.
62  RETURN
63  ENDIF
64  paru(61)=sqrt(max(0.,ps(3,4)**2-ps(3,1)**2-ps(3,2)**2-ps(3,3)**2))
65 
66 C...Find largest eigenvalue to matrix (third degree equation).
67  DO 160 j1=1,3
68  DO 160 j2=j1,3
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)
76 
77 C...Find largest eigenvector by solving equation system.
78  DO 170 j1=1,3
79  sm(j1,j1)=sm(j1,j1)-sma
80  DO 170 j2=j1+1,3
81  170 sm(j2,j1)=sm(j1,j2)
82  smax=0.
83  DO 180 j1=1,3
84  DO 180 j2=1,3
85  IF(abs(sm(j1,j2)).LE.smax) goto 180
86  ja=j1
87  jb=j2
88  smax=abs(sm(j1,j2))
89  180 CONTINUE
90  smax=0.
91  DO 190 j3=ja+1,ja+2
92  j1=j3-3*((j3-1)/3)
93  rl=sm(j1,jb)/sm(ja,jb)
94  DO 190 j2=1,3
95  sm(j1,j2)=sm(j1,j2)-rl*sm(ja,j2)
96  IF(abs(sm(j1,j2)).LE.smax) goto 190
97  jc=j1
98  smax=abs(sm(j1,j2))
99  190 CONTINUE
100  jb1=jb+1-3*(jb/3)
101  jb2=jb+2-3*((jb+1)/3)
102  sax(jb1)=-sm(jc,jb2)
103  sax(jb2)=sm(jc,jb1)
104  sax(jb)=-(sm(ja,jb1)*sax(jb1)+sm(ja,jb2)*sax(jb2))/sm(ja,jb)
105 
106 C...Divide particles into two initial clusters by hemisphere.
107  DO 200 i=n+1,n+np
108  psax=p(i,1)*sax(1)+p(i,2)*sax(2)+p(i,3)*sax(3)
109  is=1
110  IF(psax.LT.0.) is=2
111  k(i,3)=is
112  DO 200 j=1,4
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)
116 
117 C...Reassign one particle at a time; find maximum decrease of m|2 sum.
118  210 pmd=0.
119  im=0
120  DO 220 j=1,4
121  220 ps(3,j)=ps(1,j)-ps(2,j)
122  DO 230 i=n+1,n+np
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)
126  IF(pmdi.LT.pmd) THEN
127  pmd=pmdi
128  im=i
129  ENDIF
130  230 CONTINUE
131 
132 C...Loop back if significant reduction in sum of m|2.
133  IF(pmd.LT.-paru(48)*pms) THEN
134  pms=pms+pmd
135  is=k(im,3)
136  DO 240 j=1,4
137  ps(is,j)=ps(is,j)-p(im,j)
138  240 ps(3-is,j)=ps(3-is,j)+p(im,j)
139  k(im,3)=3-is
140  goto 210
141  ENDIF
142 
143 C...Final masses and output.
144  mstu(61)=n+1
145  mstu(62)=np
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))
150 
151  RETURN
152  END