Protection against floating point exception in cosh.
[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
1a434afd 21 real arg
471f69dc 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
91c*-sort cells by Et decending order, store the order in index
92 call sortzv(etc,index,ncell,-1,1,0)
93c*-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
ae70b416 105 print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell
471f69dc 106c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
107 999 continue
108c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
b9d0a01d 109cfca call vzero(flag,ncell)
471f69dc 110 do i=1, ncell
b9d0a01d 111 flag(i)=0
471f69dc 112 if(etc(i)-et_ave .le. et_min) flag(i)=2
113 enddo
114 njet = 0
115c*-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
136c*-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
141c*-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
153c*-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
166c*-New weighted eta and phi including this cell
167 eta=eta0+etas/ets
168 phi=phi0+phis/ets
169c*-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
174c*-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
181c*-Store this loop information
182 etab=eta
183 phib=phi
184 endif
185 endif
186 endif
187 endif
188 enddo
189 159 continue
190
191c*-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
ae70b416 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
471f69dc 227! 5-oct-2001 by PAI
228
229c*-reject cluster below minimum Ej_min
1a434afd 230c* protection (am)
471f69dc 231 etas=eta+etas/ets
1a434afd 232 arg = 0.
233 if (arg .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
471f69dc 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
247c*-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
272c*-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
283c*-if background level changes more than prec_bg, go next iteration!!!
284c*-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
ae70b416 290 et_sum_old = et_sum
291 print *,'End of iteration : et_ave ', et_ave, ' nc ', nc
292 + , ' Jet energy ', ets
471f69dc 293 goto 999
294c*-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