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