]> git.uio.no Git - u/mrichter/AliRoot.git/blob - 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
1 c
2 c => Original copy from /star/emc/shester/cl/pams/emc/jet/jet_finer_ua1.F
3 c
4 c:>------------------------------------------------------------------
5 C:ROUTINE:      subroutine jet_finder_ua1
6 C:DESCRIPTION:  UA1 jet algorithm from LUND JETSET
7 C:RETURN VALUE: ierror=1 on error
8 c:>------------------------------------------------------------------
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
90 c*-sort cells by Et decending order, store the order in index
91       call sortzv(etc,index,ncell,-1,1,0)   
92 c*-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 
105 c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
106  999  continue
107 c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
108 cfca      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
114 c*-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
126 C         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
138 c*-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
143 c*-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
155 c*-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
168 c*-New weighted eta and phi including this cell
169                            eta=eta0+etas/ets
170                            phi=phi0+phis/ets                          
171 c*-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
176 c*-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
183 c*-Store this loop information
184                            etab=eta
185                            phib=phi
186                         endif
187                      endif
188                   endif
189                endif
190             enddo
191  159        continue 
192             
193 c*-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             
231 c*-reject cluster below minimum Ej_min
232 c* 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
249 c*-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       
274 c*-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
285 c*-if background level changes more than prec_bg, go next iteration!!!
286 c*-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
297 c*-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