]> git.uio.no Git - u/mrichter/AliRoot.git/blame - JETAN/ua1_jet_finder.F
First commit of new jet reconstruction and analysis code to be used for the
[u/mrichter/AliRoot.git] / JETAN / ua1_jet_finder.F
CommitLineData
99e5fe42 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
7C:RETURN VALUE: ierror=1 on error
8c:>------------------------------------------------------------------
9 subroutine ua1_jet_finder(
10 + ncell, ncell_tot, etc, etac, phic,
11 + min_move, max_move, mode, prec_bg, ierror)
12 implicit none ! 5-oct-2001
13
14 real C_PI
15 real C_2PI
16 real etaCellSize
17 real phiCellSize
18 real arg
19 INTEGER NMAX,JMAX
20 parameter(NMAX=60000,JMAX=100) ! 10-oct-2201
21 integer ncell, ierror, mode, ncell_tot
22 real etc(NMAX),etac(NMAX),phic(NMAX)
23 real cone_rad, et_seed, ej_min, et_min
24 real min_move, max_move, prec_bg
25 integer flag(NMAX),index(NMAX)
26 integer n_iter,i,j,k,l,nc
27 real et_sum,et_ave,et_sum_old
28 real et,eta,phi,eta0,phi0,etab,phib,ets,etas,phis
29 real deta,dphi,r_int
30!
31 real occupationAll, occupationInJet ! 3-oct-2001 by PAI
32 integer maxTowerInJet ! 5-oct-2001 by PAI
33 save maxTowerInJet ! This is obligatory
34 integer idPerfomance
35 data idPerfomance /119/
36 save idPerfomance
37! print parameter - 27-sep-2001 by PAI
38 integer kpri
39 data kpri /0/
40 integer njet, ncellj
41 real etj, etaj, phij
42* Results
43 COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2),
44 + NCELLJ(100)
45* Cell Geometry
46 COMMON /UA1CELL/ etaCellSize, phiCellSize
47* Parameters
48 COMMON /UA1PARA/ 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
104 print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell
105c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
106 999 continue
107c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
108cfca call vzero(flag,ncell)
109 do i=1, ncell
110 flag(i)=0
111 if(etc(i)-et_ave .le. et_min) flag(i)=2
112 enddo
113 njet = 0
114c*-Initiator cell is the one with highest Et of not yet used ones
115 i=1
116 j=index(i)
117 if(i.eq.1. and. etc(j).lt.et_seed) then
118 if(mode.eq.0 .or. (mode.eq.1 .and. n_iter.eq.1)) then
119 call hf1(idPerfomance, 2., 1.)
120 print *,' no cell with Et higher than et_seed ->', etc(j)
121 return
122 endif
123 endif
124 do while(etc(j) .ge. et_seed)
125 if(flag(j) .eq. 0) then
126C write(6,*) j, etc(j), et_seed, etac(j), phic(j), flag(j)
127
128 et =etc(j)-et_ave
129 eta=etac(j)
130 phi=phic(j)
131 eta0=eta
132 phi0=phi
133 etab=eta
134 phib=phi
135 ets =0.0
136 etas=0.0
137 phis=0.0
138c*-weighted eta and phi.
139 do k = 1, ncell
140 l = index(k)
141 if(flag(l).eq.0) then
142 deta = etac(l)-eta
143c*-Is the cell is in the cone?
144 if(abs(deta).le.cone_rad)then
145 dphi=phic(l)-phi
146 do while(dphi .gt. C_PI)
147 dphi=dphi-C_2PI
148 enddo
149 do while(dphi .le. -C_PI)
150 dphi=dphi+C_2PI
151 enddo
152 if(abs(dphi).le.cone_rad) then
153 r_int=sqrt(deta**2+dphi**2)
154 if(r_int.le.cone_rad)then
155c*-calculate offset from initiate cell
156 deta=etac(l)-eta0
157 dphi=phic(l)-phi0
158 do while(dphi .gt. C_PI)
159 dphi=dphi-C_2PI
160 enddo
161 do while(dphi .lt. -C_PI)
162 dphi=dphi+C_2PI
163 enddo
164 et=etc(l)-et_ave
165 etas=etas+abs(et)*deta
166 phis=phis+abs(et)*dphi
167 ets =ets +et
168c*-New weighted eta and phi including this cell
169 eta=eta0+etas/ets
170 phi=phi0+phis/ets
171c*-If cone does not move much from previous cone, just go next step
172 r_int=sqrt((eta-etab)**2+(phi-phib)**2)
173 if(r_int .le. min_move) then
174 goto 159
175 endif
176c*-Cone should not move more than MAX_CONE_MOVE from initiator cell
177 r_int=sqrt((etas/ets)**2+(phis/ets)**2)
178 if(r_int .ge. max_move) then
179 eta=etab
180 phi=phib
181 goto 159
182 endif
183c*-Store this loop information
184 etab=eta
185 phib=phi
186 endif
187 endif
188 endif
189 endif
190 enddo
191 159 continue
192
193c*-sum up unused cells within required distance of given eta/phi
194 nc=0
195 ets=0.0
196 etas=0.0
197 phis=0.0
198 do k=1,ncell
199 l=index(k)
200 if(flag(l) .eq. 0) then
201 deta=etac(l)-eta
202 if(abs(deta).le.cone_rad)then
203 dphi=phic(l)-phi
204 do while(dphi .gt. C_PI)
205 dphi=dphi-C_2PI
206 enddo
207 do while(dphi .le. -C_PI)
208 dphi=dphi+C_2PI
209 enddo
210 if(abs(dphi).le.cone_rad) then
211 r_int=sqrt(deta**2+dphi**2)
212 if(r_int.le.cone_rad)then
213 flag(l)=-1
214 et =etc(l)-et_ave
215 ets =ets +et
216 etas=etas+et*deta
217 phis=phis+et*dphi
218 nc = nc + 1
219 endif
220 endif
221 endif
222 endif
223 enddo ! do k=1,ncell
224! 5-oct-2001 by PAI - remove 20-feb-2002 by PAI
225! 20-feb-2002 - it is work if you apply cut on eT before jet finder !!!
226! if(maxTowerInJet .gt. nc) then
227! ets = ets - et_ave*(maxTowerInJet - nc)
228! endif
229! 5-oct-2001 by PAI
230
231c*-reject cluster below minimum Ej_min
232c* protection (am)
233 etas=eta+etas/ets
234 arg = 0.
235 if (ets .ne. 0.) then
236 if (abs(etas/ets) .lt. 23.719) then
237 arg = ets * cosh(etas/ets)
238 else
239 arg = 1.e10
240 endif
241 endif
242
243 if(arg .lt. ej_min) then
244 do k=1,ncell
245 if(flag(k).le.0) flag(k)=0
246 enddo
247 if(njet.eq.0) call hf1(idPerfomance, 3., 1.)
248 else
249c*-eles, store flags and jet variables
250 do k=1,ncell
251 if(flag(k).eq.-1) flag(k)=1
252 enddo
253 phi=phi+phis/ets
254 do while(phi .ge. C_2PI)
255 phi=phi-C_2PI
256 enddo
257 do while(phi .lt. 0.0)
258 phi=phi+C_2PI
259 enddo
260 njet=njet+1
261 etj(njet) =ets
262 etaj(njet,1)=eta0
263 phij(njet,1)=phi0
264 etaj(njet,2)=etas
265 phij(njet,2)=phi
266 ncellj(njet)=nc
267 call hf1(112, float(nc)/float(maxTowerInJet), 1.) ! 8-oct-2001
268 endif
269 endif
270 i=i+1
271 j=index(i)
272 enddo
273
274c*-recalculate energy sum excluding used cells.
275 if(mode.eq.1)then
276 et_sum=0.0
277 nc=0 ! #cells in jets
278 do i=1,ncell
279 if(flag(i).ne.1) then ! 1 if cell in jet
280 et_sum=et_sum+etc(i)
281 else
282 nc=nc+1
283 endif
284 enddo
285c*-if background level changes more than prec_bg, go next iteration!!!
286c*-after 10 iteration, stop working and finish
287 if( (et_sum .gt. 0.)
288 + .and. (abs(et_sum-et_sum_old)/et_sum.gt.prec_bg
289 + .and. n_iter.le.10)
290 + .or. n_iter.eq.1) then ! minimum 2 iteration - 10-oct-2001 by pai
291 et_ave=et_sum/float(ncell_tot-nc)
292 n_iter=n_iter+1
293 et_sum_old = et_sum
294 print *,'End of iteration : et_ave ', et_ave, ' nc ', nc
295 + , ' Jet energy ', ets
296 goto 999
297c*-Watch out!!! Here is a big jump!!!
298 endif
299 occupationInJet = float(ncell) / float(ncell_tot)
300 call hf1(111, occupationInJet, 1.)
301 write(*,*) njet,' jet found in ',n_iter,
302 + ' iteration(s) : EtSum, EtAve =',
303 + et_sum,et_sum/float(ncell_tot-nc)
304 + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
305 + , ' occupationAll ', occupationAll
306 + , ' occupationInJet ', occupationInJet
307 endif
308
309 if(njet.gt.100)then
310 write(*,*)'UA1:Problem:More than 100 jets found!'
311 ierror = 1
312 elseif(njet.eq.0)then
313 write(*,*)'UA1:Done:No jet found!'
314 else
315 write(*,*)
316 +'UA1:Done:Found ',njet,' jet(s)'
317 end if
318 return
319 end
320
321
322