SDigits has become a copy of hits and threshold to associate primary particles has...
[u/mrichter/AliRoot.git] / EMCAL / jet_finder_ua1.F
CommitLineData
471f69dc 1c
2c => Original copy from /star/emc/shester/cl/pams/emc/jet/jet_finer_ua1.F
3c
4c:>------------------------------------------------------------------
5C:ROUTINE: subroutine jet_finder_ua1
6C:DESCRIPTION: UA1 jet algorithm from LUND JETSET called from EMC_erj
7C:ARGUMENTS: emc_energy_h, emc_energy,
8C:ARGUMENTS: et_tot, EMC_jetpar, ierror
9C:RETURN VALUE: ierror=1 on error
10c:>------------------------------------------------------------------
11 subroutine jet_finder_ua1(
12 + ncell, ncell_tot, etc, etac, phic,
13 + min_move, max_move, mode, prec_bg, ierror)
14 implicit none ! 5-oct-2001
15c#include "math_constants.inc"
16c#include "calorSize.inc"
17 real C_PI
18 real C_2PI
19 real etaCellSize
20 real phiCellSize
21 INTEGER NMAX,JMAX
22 parameter(NMAX=30000,JMAX=100) ! 10-oct-2201
23 integer ncell, ierror, mode, ncell_tot
24 real etc(NMAX),etac(NMAX),phic(NMAX)
25 real cone_rad, et_seed, ej_min, et_min
26 real min_move, max_move, prec_bg
27 integer flag(NMAX),index(NMAX)
28 integer n_iter,i,j,k,l,nc
29 real et_sum,et_ave,et_sum_old
30 real et,eta,phi,eta0,phi0,etab,phib,ets,etas,phis
31 real deta,dphi,r_int
32!
33 real occupationAll, occupationInJet ! 3-oct-2001 by PAI
34 integer maxTowerInJet ! 5-oct-2001 by PAI
35 save maxTowerInJet ! This is obligatory
36 integer idPerfomance /119/
37 save idPerfomance
38! print parameter - 27-sep-2001 by PAI
39 integer kpri/0/
40 integer njet, ncellj
41 real etj, etaj, phij
42* Results
43 COMMON /EMCALJETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2),
44 + NCELLJ(100)
45* Cell Geometry
46 COMMON /EMCALCELLGEO/ etaCellSize, phiCellSize
47* Parameters
48 COMMON /EMCALJETPARAM/ cone_rad, et_seed, ej_min, et_min
49
50 C_PI = 3.1415926
51 C_2PI = 2.*C_PI
52
53 if(kpri.eq.0) then
54 kpri = 1
55! for add. correction of jet energy if ocupation in jet not 100% !!
56! may differ from real number because grid
57 maxTowerInJet=
58 + int((C_PI*cone_rad*cone_rad)/(etaCellSize*phiCellSize)+0.5)
59 print 1, ncell_tot
60 + ,cone_rad, et_seed, ej_min, et_min
61 + ,min_move, max_move, mode, prec_bg
62 + ,maxTowerInJet
63 1 format(/
64 + ' == jet_finder_UA1 == '/
65 + ' ncell_tot ', i5/
66 + ' cone rad ', f5.2/
67 + ' et_seed ', f5.2,' GeV/C'/
68 + ' ej_min ', f5.2,' GeV/C'/
69 + ' et_min(tower after bg sub.) ', f5.2,' GeV/C'/
70 + ' min_cone_move ', f5.3/
71 + ' max_cone_move ', f5.3/
72 + ' Mode for BG subtraction ', i5/
73 + ' prec_bg ', f5.4/
74 + ' -- PAI"s addition --'/
75 + ' maxTowerInJet ', i5/
76 + ' -- '/
77 + )
78 if(NMAX .lt. ncell_tot) then
79 print 2, NMAX, ncell_tot
80 2 format('<E> you must increase memory -> NMAX ',i6
81 + ,' ncell_tot ', i6)
82 endif
83 endif
84 call hf1(idPerfomance, 1., 1.)
85 occupationAll = float(ncell) / float(ncell_tot)
86 call hf1(111, occupationAll, 1.)
87! print parameter - 27-sep-2001 by PAI
88 ierror =0
89 n_iter =0
90c*-sort cells by Et decending order, store the order in index
91 call sortzv(etc,index,ncell,-1,1,0)
92c*-sum up total energy of all cells
93 if(mode .eq. 1) then
94 n_iter=1
95 et_sum=0.0
96 do i=1, ncell
97 et_sum=et_sum+etc(i)
98 enddo
99 et_sum_old=et_sum
100 et_ave=et_sum/float(ncell_tot)
101 else
102 et_ave=0.0
103 endif
104c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
105 999 continue
106c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
107 call vzero(flag,ncell)
108 do i=1, ncell
109 if(etc(i)-et_ave .le. et_min) flag(i)=2
110 enddo
111 njet = 0
112c*-Initiator cell is the one with highest Et of not yet used ones
113 i=1
114 j=index(i)
115 if(i.eq.1. and. etc(j).lt.et_seed) then
116 if(mode.eq.0 .or. (mode.eq.1 .and. n_iter.eq.1)) then
117 call hf1(idPerfomance, 2., 1.)
118 print *,' no cell with Et higher than et_seed ->', etc(j)
119 endif
120 endif
121 do while(etc(j) .ge. et_seed)
122 if(flag(j) .eq. 0) then
123 et =etc(j)-et_ave
124 eta=etac(j)
125 phi=phic(j)
126 eta0=eta
127 phi0=phi
128 etab=eta
129 phib=phi
130 ets =0.0
131 etas=0.0
132 phis=0.0
133c*-weighted eta and phi.
134 do k= 1, ncell
135 l=index(k)
136 if(flag(l).eq.0) then
137 deta=etac(l)-eta
138c*-Is the cell is in the cone?
139 if(abs(deta).le.cone_rad)then
140 dphi=phic(l)-phi
141 do while(dphi .gt. C_PI)
142 dphi=dphi-C_2PI
143 enddo
144 do while(dphi .le. -C_PI)
145 dphi=dphi+C_2PI
146 enddo
147 if(abs(dphi).le.cone_rad) then
148 r_int=sqrt(deta**2+dphi**2)
149 if(r_int.le.cone_rad)then
150c*-calculate offset from initiate cell
151 deta=etac(l)-eta0
152 dphi=phic(l)-phi0
153 do while(dphi .gt. C_PI)
154 dphi=dphi-C_2PI
155 enddo
156 do while(dphi .lt. -C_PI)
157 dphi=dphi+C_2PI
158 enddo
159 et=etc(l)-et_ave
160 etas=etas+abs(et)*deta
161 phis=phis+abs(et)*dphi
162 ets =ets +et
163c*-New weighted eta and phi including this cell
164 eta=eta0+etas/ets
165 phi=phi0+phis/ets
166c*-If cone does not move much from previous cone, just go next step
167 r_int=sqrt((eta-etab)**2+(phi-phib)**2)
168 if(r_int .le. min_move) then
169 goto 159
170 endif
171c*-Cone should not move more than MAX_CONE_MOVE from initiator cell
172 r_int=sqrt((etas/ets)**2+(phis/ets)**2)
173 if(r_int .ge. max_move) then
174 eta=etab
175 phi=phib
176 goto 159
177 endif
178c*-Store this loop information
179 etab=eta
180 phib=phi
181 endif
182 endif
183 endif
184 endif
185 enddo
186 159 continue
187
188c*-sum up unused cells within required distance of given eta/phi
189 nc=0
190 ets=0.0
191 etas=0.0
192 phis=0.0
193 do k=1,ncell
194 l=index(k)
195 if(flag(l) .eq. 0) then
196 deta=etac(l)-eta
197 if(abs(deta).le.cone_rad)then
198 dphi=phic(l)-phi
199 do while(dphi .gt. C_PI)
200 dphi=dphi-C_2PI
201 enddo
202 do while(dphi .le. -C_PI)
203 dphi=dphi+C_2PI
204 enddo
205 if(abs(dphi).le.cone_rad) then
206 r_int=sqrt(deta**2+dphi**2)
207 if(r_int.le.cone_rad)then
208 flag(l)=-1
209 et =etc(l)-et_ave
210 ets =ets +et
211 etas=etas+et*deta
212 phis=phis+et*dphi
213 nc = nc + 1
214 endif
215 endif
216 endif
217 endif
218 enddo ! do k=1,ncell
219! 5-oct-2001 by PAI
220 if(maxTowerInJet .gt. nc) then
221 ets = ets - et_ave*(maxTowerInJet - nc)
222 endif
223! 5-oct-2001 by PAI
224
225c*-reject cluster below minimum Ej_min
226 etas=eta+etas/ets
227 if(ets*cosh(etas/ets).lt.ej_min) then
228 do k=1,ncell
229 if(flag(k).le.0) flag(k)=0
230 enddo
231 if(njet.eq.0) call hf1(idPerfomance, 3., 1.)
232 else
233c*-eles, store flags and jet variables
234 do k=1,ncell
235 if(flag(k).eq.-1) flag(k)=1
236 enddo
237 phi=phi+phis/ets
238 do while(phi .ge. C_2PI)
239 phi=phi-C_2PI
240 enddo
241 do while(phi .lt. 0.0)
242 phi=phi+C_2PI
243 enddo
244 njet=njet+1
245 etj(njet) =ets
246 etaj(njet,1)=eta0
247 phij(njet,1)=phi0
248 etaj(njet,2)=etas
249 phij(njet,2)=phi
250 ncellj(njet)=nc
251 call hf1(112, float(nc)/float(maxTowerInJet), 1.) ! 8-oct-2001
252 endif
253 endif
254 i=i+1
255 j=index(i)
256 enddo
257
258c*-recalculate energy sum excluding used cells.
259 if(mode.eq.1)then
260 et_sum=0.0
261 nc=0 ! #cells in jets
262 do i=1,ncell
263 if(flag(i).ne.1) then ! 1 if cell in jet
264 et_sum=et_sum+etc(i)
265 else
266 nc=nc+1
267 endif
268 enddo
269c*-if background level changes more than prec_bg, go next iteration!!!
270c*-after 10 iteration, stop working and finish
271 if( (abs(et_sum-et_sum_old)/et_sum.gt.prec_bg
272 + .and. n_iter.le.10)
273 + .or. n_iter.eq.1) then ! minimum 2 iteration - 10-oct-2001 by pai
274 et_ave=et_sum/float(ncell_tot-nc)
275 n_iter=n_iter+1
276 et_sum_old=et_sum
277 goto 999
278c*-Watch out!!! Here is a big jump!!!
279 endif
280 occupationInJet = float(ncell) / float(ncell_tot)
281 call hf1(111, occupationInJet, 1.)
282 write(*,*) njet,' jet found in ',n_iter,
283 + ' iteration(s) : EtSum, EtAve =',
284 + et_sum,et_sum/float(ncell_tot-nc)
285 + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
286 + , ' occupationAll ', occupationAll
287 + , ' occupationInJet ', occupationInJet
288 endif
289
290 if(njet.gt.100)then
291 write(*,*)'UA1:Problem:More than 100 jets found!'
292 ierror = 1
293 elseif(njet.eq.0)then
294 write(*,*)'UA1:Done:No jet found!'
295 else
296 write(*,*)
297 +'UA1:Done:Found ',njet,' jet(s)'
298 end if
299 return
300 end
301
302
303