Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lucell.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lucell.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE lucell(NJET)
5 
6 C...Purpose: to provide a simple way of jet finding in an eta-phi-ET
7 C...coordinate frame, as used for calorimeters at hadron colliders.
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 
15 C...Loop over all particles. Find cell that was hit by given particle.
16  nce2=2*mstu(51)*mstu(52)
17  ptlrat=1./sinh(paru(51))**2
18  np=0
19  nc=n
20  DO 110 i=1,n
21  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 110
22  IF(p(i,1)**2+p(i,2)**2.LE.ptlrat*p(i,3)**2) goto 110
23  IF(mstu(41).GE.2) THEN
24  kc=lucomp(k(i,2))
25  IF(kc.EQ.0.OR.kc.EQ.12.OR.kc.EQ.14.OR.kc.EQ.16.OR.
26  & kc.EQ.18) goto 110
27  IF(mstu(41).GE.3.AND.kchg(kc,2).EQ.0.AND.luchge(k(i,2)).EQ.0)
28  & goto 110
29  ENDIF
30  np=np+1
31  pt=sqrt(p(i,1)**2+p(i,2)**2)
32  eta=sign(log((sqrt(pt**2+p(i,3)**2)+abs(p(i,3)))/pt),p(i,3))
33  ieta=max(1,min(mstu(51),1+int(mstu(51)*0.5*(eta/paru(51)+1.))))
34  phi=ulangl(p(i,1),p(i,2))
35  iphi=max(1,min(mstu(52),1+int(mstu(52)*0.5*(phi/paru(1)+1.))))
36  ietph=mstu(52)*ieta+iphi
37 
38 C...Add to cell already hit, or book new cell.
39  DO 100 ic=n+1,nc
40  IF(ietph.EQ.k(ic,3)) THEN
41  k(ic,4)=k(ic,4)+1
42  p(ic,5)=p(ic,5)+pt
43  goto 110
44  ENDIF
45  100 CONTINUE
46  IF(nc.GE.mstu(4)-mstu(32)-5) THEN
47  CALL luerrm(11,'(LUCELL:) no more memory left in LUJETS')
48  njet=-2
49  RETURN
50  ENDIF
51  nc=nc+1
52  k(nc,3)=ietph
53  k(nc,4)=1
54  k(nc,5)=2
55  p(nc,1)=(paru(51)/mstu(51))*(2*ieta-1-mstu(51))
56  p(nc,2)=(paru(1)/mstu(52))*(2*iphi-1-mstu(52))
57  p(nc,5)=pt
58  110 CONTINUE
59 
60 C...Smear true bin content by calorimeter resolution.
61  IF(mstu(53).GE.1) THEN
62  DO 130 ic=n+1,nc
63  pei=p(ic,5)
64  IF(mstu(53).EQ.2) pei=p(ic,5)/cosh(p(ic,1))
65  120 pef=pei+paru(55)*sqrt(-2.*log(max(1e-10,rlu(0)))*pei)*
66  & cos(paru(2)*rlu(0))
67  IF(pef.LT.0..OR.pef.GT.paru(56)*pei) goto 120
68  p(ic,5)=pef
69  130 IF(mstu(53).EQ.2) p(ic,5)=pef*cosh(p(ic,1))
70  ENDIF
71 
72 C...Find initiator cell: the one with highest pT of not yet used ones.
73  nj=nc
74  140 etmax=0.
75  DO 150 ic=n+1,nc
76  IF(k(ic,5).NE.2) goto 150
77  IF(p(ic,5).LE.etmax) goto 150
78  icmax=ic
79  eta=p(ic,1)
80  phi=p(ic,2)
81  etmax=p(ic,5)
82  150 CONTINUE
83  IF(etmax.LT.paru(52)) goto 210
84  IF(nj.GE.mstu(4)-mstu(32)-5) THEN
85  CALL luerrm(11,'(LUCELL:) no more memory left in LUJETS')
86  njet=-2
87  RETURN
88  ENDIF
89  k(icmax,5)=1
90  nj=nj+1
91  k(nj,4)=0
92  k(nj,5)=1
93  p(nj,1)=eta
94  p(nj,2)=phi
95  p(nj,3)=0.
96  p(nj,4)=0.
97  p(nj,5)=0.
98 
99 C...Sum up unused cells within required distance of initiator.
100  DO 160 ic=n+1,nc
101  IF(k(ic,5).EQ.0) goto 160
102  IF(abs(p(ic,1)-eta).GT.paru(54)) goto 160
103  dphia=abs(p(ic,2)-phi)
104  IF(dphia.GT.paru(54).AND.dphia.LT.paru(2)-paru(54)) goto 160
105  phic=p(ic,2)
106  IF(dphia.GT.paru(1)) phic=phic+sign(paru(2),phi)
107  IF((p(ic,1)-eta)**2+(phic-phi)**2.GT.paru(54)**2) goto 160
108  k(ic,5)=-k(ic,5)
109  k(nj,4)=k(nj,4)+k(ic,4)
110  p(nj,3)=p(nj,3)+p(ic,5)*p(ic,1)
111  p(nj,4)=p(nj,4)+p(ic,5)*phic
112  p(nj,5)=p(nj,5)+p(ic,5)
113  160 CONTINUE
114 
115 C...Reject cluster below minimum ET, else accept.
116  IF(p(nj,5).LT.paru(53)) THEN
117  nj=nj-1
118  DO 170 ic=n+1,nc
119  170 IF(k(ic,5).LT.0) k(ic,5)=-k(ic,5)
120  ELSEIF(mstu(54).LE.2) THEN
121  p(nj,3)=p(nj,3)/p(nj,5)
122  p(nj,4)=p(nj,4)/p(nj,5)
123  IF(abs(p(nj,4)).GT.paru(1)) p(nj,4)=p(nj,4)-sign(paru(2),
124  & p(nj,4))
125  DO 180 ic=n+1,nc
126  180 IF(k(ic,1).LT.0) k(ic,1)=0
127  ELSE
128  DO 190 j=1,4
129  190 p(nj,j)=0.
130  DO 200 ic=n+1,nc
131  IF(k(ic,5).GE.0) goto 200
132  p(nj,1)=p(nj,1)+p(ic,5)*cos(p(ic,2))
133  p(nj,2)=p(nj,2)+p(ic,5)*sin(p(ic,2))
134  p(nj,3)=p(nj,3)+p(ic,5)*sinh(p(ic,1))
135  p(nj,4)=p(nj,4)+p(ic,5)*cosh(p(ic,1))
136  k(ic,5)=0
137  200 CONTINUE
138  ENDIF
139  goto 140
140 
141 C...Arrange clusters in falling ET sequence.
142  210 DO 230 i=1,nj-nc
143  etmax=0.
144  DO 220 ij=nc+1,nj
145  IF(k(ij,5).EQ.0) goto 220
146  IF(p(ij,5).LT.etmax) goto 220
147  ijmax=ij
148  etmax=p(ij,5)
149  220 CONTINUE
150  k(ijmax,5)=0
151  k(n+i,1)=31
152  k(n+i,2)=98
153  k(n+i,3)=i
154  k(n+i,4)=k(ijmax,4)
155  k(n+i,5)=0
156  DO 230 j=1,5
157  p(n+i,j)=p(ijmax,j)
158  230 v(n+i,j)=0.
159  njet=nj-nc
160 
161 C...Convert to massless or massive four-vectors.
162  IF(mstu(54).EQ.2) THEN
163  DO 240 i=n+1,n+njet
164  eta=p(i,3)
165  p(i,1)=p(i,5)*cos(p(i,4))
166  p(i,2)=p(i,5)*sin(p(i,4))
167  p(i,3)=p(i,5)*sinh(eta)
168  p(i,4)=p(i,5)*cosh(eta)
169  240 p(i,5)=0.
170  ELSEIF(mstu(54).GE.3) THEN
171  DO 250 i=n+1,n+njet
172  250 p(i,5)=sqrt(max(0.,p(i,4)**2-p(i,1)**2-p(i,2)**2-p(i,3)**2))
173  ENDIF
174 
175 C...Information about storage.
176  mstu(61)=n+1
177  mstu(62)=np
178  mstu(63)=nc-n
179  IF(mstu(43).LE.1) mstu(3)=njet
180  IF(mstu(43).GE.2) n=n+njet
181 
182  RETURN
183  END