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
  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