JETAN/ua1_jet_finder.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
   7 C:RETURN VALUE: ierror=1 on error
   8 c:>------------------------------------------------------------------
   9       subroutine ua1_jet_finder(
  10      +        ncell, ncell_tot, etc, etac, phic,
  11      +        min_move, max_move, mode, prec_bg, ierror)
  12       implicit none ! 5-oct-2001
  13
  14       real C_PI
  15       real C_2PI
  16       real etaCellSize
  17       real phiCellSize
  18       real arg
  19       INTEGER NMAX,JMAX
  20       parameter(NMAX=60000,JMAX=100) ! 10-oct-2201
  21       integer ncell, ierror, mode, ncell_tot
  22       real etc(NMAX),etac(NMAX),phic(NMAX)
  23       real cone_rad, et_seed, ej_min, et_min
  24       real min_move, max_move, prec_bg
  25       integer flag(NMAX),index(NMAX)
  26       integer n_iter,i,j,k,l,nc
  27       real et_sum,et_ave,et_sum_old
  28       real et,eta,phi,eta0,phi0,etab,phib,ets,etas,phis
  29       real deta,dphi,r_int
  30 !
  31       real occupationAll, occupationInJet  ! 3-oct-2001 by PAI
  32       integer maxTowerInJet                ! 5-oct-2001 by PAI
  33       save    maxTowerInJet                ! This is obligatory
  34       integer idPerfomance
  35       data    idPerfomance /119/
  36       save    idPerfomance
  37 !  print parameter - 27-sep-2001 by PAI
  38       integer kpri
  39       data    kpri /0/
  40       integer njet, ncellj
  41       real    etj, etaj, phij, etavg
  42 *    Results
  43       COMMON /UA1JETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2),
  44      +     NCELLJ(100), ETAVG
  45 *    Cell Geometry
  46       COMMON /UA1CELL/ etaCellSize, phiCellSize
  47 *    Parameters
  48       COMMON /UA1PARA/ 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       occupationAll = float(ncell) / float(ncell_tot)
  85 ! print parameter - 27-sep-2001 by PAI
  86       ierror =0
  87       n_iter =0
  88 c*-sort cells by Et decending order, store the order in index
  89       call sortzv(etc,index,ncell,-1,1,0)
  90 c*-sum up total energy of all cells
  91       if(mode .eq. 1) then
  92          n_iter=1
  93          et_sum=0.0
  94          do i=1, ncell
  95             et_sum=et_sum+etc(i)
  96          enddo
  97          et_sum_old=et_sum
  98          et_ave=et_sum/float(ncell_tot)
  99       else
 100          et_ave=0.0
 101       endif
 102       print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell
 103 c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
 104  999  continue
 105 c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
 106 cfca      call vzero(flag,ncell)
 107       do i=1, ncell
 108          flag(i)=0
 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           print *,' no cell with Et higher than et_seed ->', etc(j)
 118           return
 119         endif
 120       endif
 121       do while(etc(j) .ge. et_seed)
 122          if(flag(j) .eq. 0) then
 123 C         write(6,*) j, etc(j), et_seed, etac(j), phic(j), flag(j)
 124
 125             et =etc(j)-et_ave
 126             eta=etac(j)
 127             phi=phic(j)
 128             eta0=eta
 129             phi0=phi
 130             etab=eta
 131             phib=phi
 132             ets =0.0
 133             etas=0.0
 134             phis=0.0
 135 c*-weighted eta and phi.
 136             do k = 1, ncell
 137                l = index(k)
 138                if(flag(l).eq.0) then
 139                   deta = etac(l)-eta
 140 c*-Is the cell is in the cone?
 141                   if(abs(deta).le.cone_rad)then
 142                      dphi=phic(l)-phi
 143                      do while(dphi .gt. C_PI)
 144                         dphi=dphi-C_2PI
 145                      enddo
 146                      do while(dphi .le. -C_PI)
 147                         dphi=dphi+C_2PI
 148                      enddo
 149                      if(abs(dphi).le.cone_rad) then
 150                         r_int=sqrt(deta**2+dphi**2)
 151                         if(r_int.le.cone_rad)then
 152 c*-calculate offset from initiate cell
 153                            deta=etac(l)-eta0
 154                            dphi=phic(l)-phi0
 155                            do while(dphi .gt. C_PI)
 156                               dphi=dphi-C_2PI
 157                            enddo
 158                            do while(dphi .lt. -C_PI)
 159                               dphi=dphi+C_2PI
 160                            enddo
 161                            et=etc(l)-et_ave
 162                            etas=etas+abs(et)*deta
 163                            phis=phis+abs(et)*dphi
 164                            ets =ets +et
 165 c*-New weighted eta and phi including this cell
 166                            eta=eta0+etas/ets
 167                            phi=phi0+phis/ets
 168 c*-If cone does not move much from previous cone, just go next step
 169                            r_int=sqrt((eta-etab)**2+(phi-phib)**2)
 170                            if(r_int .le. min_move) then
 171                               goto 159
 172                            endif
 173 c*-Cone should not move more than MAX_CONE_MOVE from initiator cell
 174                            r_int=sqrt((etas/ets)**2+(phis/ets)**2)
 175                            if(r_int .ge. max_move) then
 176                               eta=etab
 177                               phi=phib
 178                               goto 159
 179                            endif
 180 c*-Store this loop information
 181                            etab=eta
 182                            phib=phi
 183                         endif
 184                      endif
 185                   endif
 186                endif
 187             enddo
 188  159        continue
 189
 190 c*-sum up unused cells within required distance of given eta/phi
 191             nc=0
 192             ets=0.0
 193             etas=0.0
 194             phis=0.0
 195             do k=1,ncell
 196                l=index(k)
 197                if(flag(l) .eq. 0) then
 198                   deta=etac(l)-eta
 199                   if(abs(deta).le.cone_rad)then
 200                      dphi=phic(l)-phi
 201                      do while(dphi .gt. C_PI)
 202                         dphi=dphi-C_2PI
 203                      enddo
 204                      do while(dphi .le. -C_PI)
 205                         dphi=dphi+C_2PI
 206                      enddo
 207                      if(abs(dphi).le.cone_rad) then
 208                         r_int=sqrt(deta**2+dphi**2)
 209                         if(r_int.le.cone_rad)then
 210                            flag(l)=-1
 211                            et  =etc(l)-et_ave
 212                            ets =ets +et
 213                            etas=etas+et*deta
 214                            phis=phis+et*dphi
 215                            nc  = nc + 1
 216                         endif
 217                      endif
 218                   endif
 219                endif
 220             enddo  ! do k=1,ncell
 221 !  5-oct-2001 by PAI - remove 20-feb-2002 by PAI
 222 ! 20-feb-2002 - it is work if you apply cut on eT before jet finder !!!
 223 !            if(maxTowerInJet .gt. nc) then
 224 !              ets = ets - et_ave*(maxTowerInJet - nc)
 225 !            endif
 226 ! 5-oct-2001 by PAI
 227
 228 c*-reject cluster below minimum Ej_min
 229 c* protection (am)
 230
 231 c            arg = 0.
 232 c            if (ets .ne. 0.) then
 233 c               if (abs(etas/ets) .lt. 23.719) then
 234 c                  arg = ets * cosh(etas/ets)
 235 c               else
 236 c                  arg = 1.e10
 237 c               endif
 238 c            endif
 239
 240             if(ets .lt. ej_min) then
 241                do k=1,ncell
 242                   if(flag(k).le.0) flag(k)=0
 243                enddo
 244             else
 245 c*-eles, store flags and jet variables
 246                do k=1,ncell
 247                   if(flag(k).eq.-1) flag(k)=1
 248                enddo
 249                etas=eta+etas/ets
 250                phi=phi+phis/ets
 251                do while(phi .ge. C_2PI)
 252                   phi=phi-C_2PI
 253                enddo
 254                do while(phi .lt. 0.0)
 255                   phi=phi+C_2PI
 256                enddo
 257                njet=njet+1
 258                etj(njet) =ets
 259                etaj(njet,1)=eta0
 260                phij(njet,1)=phi0
 261                etaj(njet,2)=etas
 262                phij(njet,2)=phi
 263                ncellj(njet)=nc
 264                etavg = et_ave
 265             endif
 266          endif
 267          i=i+1
 268          j=index(i)
 269       enddo
 270
 271 c*-recalculate energy sum excluding used cells.
 272       if(mode.eq.1)then
 273          et_sum=0.0
 274          nc=0       ! #cells in jets
 275          do i=1,ncell
 276             if(flag(i).ne.1) then  ! 1 if cell in jet
 277                et_sum=et_sum+etc(i)
 278             else
 279                nc=nc+1
 280             endif
 281          enddo
 282 c*-if background level changes more than prec_bg, go next iteration!!!
 283 c*-after 10 iteration, stop working and finish
 284          if( (et_sum .gt. 0.)
 285      +        .and. (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          write(*,*) njet,' jet found in ',n_iter,
 298      +     ' iteration(s) : EtSum, EtAve =',
 299      +     et_sum,et_sum/float(ncell_tot-nc)
 300      + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
 301      + , ' occupationAll ', occupationAll
 302      + , ' occupationInJet ', occupationInJet
 303       endif
 304
 305       if(njet.gt.100)then
 306          write(*,*)'UA1:Problem:More than 100 jets found!'
 307          ierror = 1
 308       elseif(njet.eq.0)then
 309          write(*,*)'UA1:Done:No jet found!'
 310       else
 311          write(*,*)
 312      +'UA1:Done:Found ',njet,' jet(s)'
 313       end if
 314       return
 315       end
 316
 317
 318