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