]>
Commit | Line | Data |
---|---|---|
45a58699 | 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 | |
38 | data idPerfomance /119/ | |
39 | save idPerfomance | |
40 | ! print parameter - 27-sep-2001 by PAI | |
41 | integer kpri | |
42 | data kpri /0/ | |
43 | integer njet, ncellj | |
44 | real etj, etaj, phij | |
45 | * Results | |
46 | COMMON /EMCALJETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2), | |
47 | + NCELLJ(100) | |
48 | * Cell Geometry | |
49 | COMMON /EMCALCELLGEO/ etaCellSize, phiCellSize | |
50 | * Parameters | |
51 | COMMON /EMCALJETPARAM/ cone_rad, et_seed, ej_min, et_min | |
52 | ||
53 | C_PI = 3.1415926 | |
54 | C_2PI = 2.*C_PI | |
55 | ||
56 | if(kpri.eq.0) then | |
57 | kpri = 1 | |
58 | ! for add. correction of jet energy if ocupation in jet not 100% !! | |
59 | ! may differ from real number because grid | |
60 | maxTowerInJet= | |
61 | + int((C_PI*cone_rad*cone_rad)/(etaCellSize*phiCellSize)+0.5) | |
62 | print 1, ncell_tot | |
63 | + ,cone_rad, et_seed, ej_min, et_min | |
64 | + ,min_move, max_move, mode, prec_bg | |
65 | + ,maxTowerInJet | |
66 | 1 format(/ | |
67 | + ' == jet_finder_UA1 == '/ | |
68 | + ' ncell_tot ', i5/ | |
69 | + ' cone rad ', f5.2/ | |
70 | + ' et_seed ', f5.2,' GeV/C'/ | |
71 | + ' ej_min ', f5.2,' GeV/C'/ | |
72 | + ' et_min(tower after bg sub.) ', f5.2,' GeV/C'/ | |
73 | + ' min_cone_move ', f5.3/ | |
74 | + ' max_cone_move ', f5.3/ | |
75 | + ' Mode for BG subtraction ', i5/ | |
76 | + ' prec_bg ', f5.4/ | |
77 | + ' -- PAI"s addition --'/ | |
78 | + ' maxTowerInJet ', i5/ | |
79 | + ' -- '/ | |
80 | + ) | |
81 | if(NMAX .lt. ncell_tot) then | |
82 | print 2, NMAX, ncell_tot | |
83 | 2 format('<E> you must increase memory -> NMAX ',i6 | |
84 | + ,' ncell_tot ', i6) | |
85 | endif | |
86 | endif | |
87 | call hf1(idPerfomance, 1., 1.) | |
88 | occupationAll = float(ncell) / float(ncell_tot) | |
89 | call hf1(111, occupationAll, 1.) | |
90 | ! print parameter - 27-sep-2001 by PAI | |
91 | ierror =0 | |
92 | n_iter =0 | |
93 | c*-sort cells by Et decending order, store the order in index | |
94 | call sortzv(etc,index,ncell,-1,1,0) | |
95 | c*-sum up total energy of all cells | |
96 | if(mode .eq. 1) then | |
97 | n_iter=1 | |
98 | et_sum=0.0 | |
99 | do i=1, ncell | |
100 | et_sum=et_sum+etc(i) | |
101 | enddo | |
102 | et_sum_old=et_sum | |
103 | et_ave=et_sum/float(ncell_tot) | |
104 | else | |
105 | et_ave=0.0 | |
106 | endif | |
107 | print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell | |
108 | c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!! | |
109 | 999 continue | |
110 | c*-kill cells (flag=2) with Et below ET_MIN after background subtraction | |
111 | cfca call vzero(flag,ncell) | |
112 | do i=1, ncell | |
113 | flag(i)=0 | |
114 | if(etc(i)-et_ave .le. et_min) flag(i)=2 | |
115 | enddo | |
116 | njet = 0 | |
117 | c*-Initiator cell is the one with highest Et of not yet used ones | |
118 | i=1 | |
119 | j=index(i) | |
120 | if(i.eq.1. and. etc(j).lt.et_seed) then | |
121 | if(mode.eq.0 .or. (mode.eq.1 .and. n_iter.eq.1)) then | |
122 | call hf1(idPerfomance, 2., 1.) | |
123 | print *,' no cell with Et higher than et_seed ->', etc(j) | |
124 | endif | |
125 | endif | |
126 | do while(etc(j) .ge. et_seed) | |
127 | if(flag(j) .eq. 0) then | |
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( (abs(et_sum-et_sum_old)/et_sum.gt.prec_bg | |
288 | + .and. n_iter.le.10) | |
289 | + .or. n_iter.eq.1) then ! minimum 2 iteration - 10-oct-2001 by pai | |
290 | et_ave=et_sum/float(ncell_tot-nc) | |
291 | n_iter=n_iter+1 | |
292 | et_sum_old = et_sum | |
293 | print *,'End of iteration : et_ave ', et_ave, ' nc ', nc | |
294 | + , ' Jet energy ', ets | |
295 | goto 999 | |
296 | c*-Watch out!!! Here is a big jump!!! | |
297 | endif | |
298 | occupationInJet = float(ncell) / float(ncell_tot) | |
299 | call hf1(111, occupationInJet, 1.) | |
300 | write(*,*) njet,' jet found in ',n_iter, | |
301 | + ' iteration(s) : EtSum, EtAve =', | |
302 | + et_sum,et_sum/float(ncell_tot-nc) | |
303 | + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc | |
304 | + , ' occupationAll ', occupationAll | |
305 | + , ' occupationInJet ', occupationInJet | |
306 | endif | |
307 | ||
308 | if(njet.gt.100)then | |
309 | write(*,*)'UA1:Problem:More than 100 jets found!' | |
310 | ierror = 1 | |
311 | elseif(njet.eq.0)then | |
312 | write(*,*)'UA1:Done:No jet found!' | |
313 | else | |
314 | write(*,*) | |
315 | +'UA1:Done:Found ',njet,' jet(s)' | |
316 | end if | |
317 | return | |
318 | end | |
319 | ||
320 | ||
321 |