[u/mrichter/AliRoot.git] / EMCAL / jet_finder_ua1.F

c
c => Original copy from /star/emc/shester/cl/pams/emc/jet/jet_finer_ua1.F
c
c:>------------------------------------------------------------------
C:ROUTINE:      subroutine jet_finder_ua1
C:DESCRIPTION:  UA1 jet algorithm from LUND JETSET called from EMC_erj
C:ARGUMENTS:    emc_energy_h, emc_energy,
C:ARGUMENTS:    et_tot, EMC_jetpar, ierror
C:RETURN VALUE: ierror=1 on error
c:>------------------------------------------------------------------
      subroutine jet_finder_ua1(
     +        ncell, ncell_tot, etc, etac, phic,
     +        min_move, max_move, mode, prec_bg, ierror)
      implicit none ! 5-oct-2001
c#include "math_constants.inc"
c#include "calorSize.inc"
      real C_PI
      real C_2PI 
      real etaCellSize 
      real phiCellSize 
      real arg
      INTEGER NMAX,JMAX
      parameter(NMAX=30000,JMAX=100) ! 10-oct-2201
      integer ncell, ierror, mode, ncell_tot
      real etc(NMAX),etac(NMAX),phic(NMAX)
      real cone_rad, et_seed, ej_min, et_min
      real min_move, max_move, prec_bg
      integer flag(NMAX),index(NMAX)
      integer n_iter,i,j,k,l,nc
      real et_sum,et_ave,et_sum_old
      real et,eta,phi,eta0,phi0,etab,phib,ets,etas,phis
      real deta,dphi,r_int
!
      real occupationAll, occupationInJet  ! 3-oct-2001 by PAI 
      integer maxTowerInJet                ! 5-oct-2001 by PAI
      save    maxTowerInJet                ! This is obligatory
      integer idPerfomance
      data    idPerfomance /119/
      save    idPerfomance
!  print parameter - 27-sep-2001 by PAI
      integer kpri
      data    kpri /0/
      integer njet, ncellj
      real    etj, etaj, phij
*    Results
      COMMON /EMCALJETS/ NJET, ETJ(100), ETAJ(100,2), PHIJ(100,2), 
     +     NCELLJ(100)
*    Cell Geometry
      COMMON /EMCALCELLGEO/ etaCellSize, phiCellSize
*    Parameters
      COMMON /EMCALJETPARAM/ cone_rad, et_seed, ej_min, et_min

      C_PI  = 3.1415926
      C_2PI = 2.*C_PI

      if(kpri.eq.0) then
        kpri = 1
! for add. correction of jet energy if ocupation in jet not 100% !!
! may differ from real number because grid 
        maxTowerInJet=
     +       int((C_PI*cone_rad*cone_rad)/(etaCellSize*phiCellSize)+0.5)
        print 1, ncell_tot
     +       ,cone_rad, et_seed, ej_min, et_min
     +       ,min_move, max_move, mode, prec_bg
     +       ,maxTowerInJet
 1    format(/
     +    '    == jet_finder_UA1 ==  '/
     +    ' ncell_tot                   ', i5/
     +    ' cone rad                    ', f5.2/
     +    ' et_seed                     ', f5.2,' GeV/C'/
     +    ' ej_min                      ', f5.2,' GeV/C'/
     +    ' et_min(tower after bg sub.) ', f5.2,' GeV/C'/
     +    ' min_cone_move               ', f5.3/
     +    ' max_cone_move               ', f5.3/
     +    ' Mode for BG subtraction     ', i5/
     +    ' prec_bg                     ', f5.4/
     +    ' -- PAI"s addition --'/
     +    ' maxTowerInJet               ', i5/
     +    ' -- '/
     +    )
        if(NMAX .lt. ncell_tot) then
          print 2, NMAX, ncell_tot
 2    format('<E> you must increase memory -> NMAX ',i6
     +      ,' ncell_tot ', i6)
        endif 
      endif
      call hf1(idPerfomance, 1., 1.)
      occupationAll = float(ncell) / float(ncell_tot)
      call hf1(111, occupationAll, 1.)
! print parameter - 27-sep-2001 by PAI
      ierror =0
      n_iter =0
c*-sort cells by Et decending order, store the order in index
      call sortzv(etc,index,ncell,-1,1,0)   
c*-sum up total energy of all cells
      if(mode .eq. 1) then
         n_iter=1
         et_sum=0.0
         do i=1, ncell
            et_sum=et_sum+etc(i)
         enddo
         et_sum_old=et_sum
         et_ave=et_sum/float(ncell_tot)
      else
         et_ave=0.0
      endif
      print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell 
c*-Watch out!!! For mode=1, it can jump back to here for next iteration!!!
 999  continue
c*-kill cells (flag=2) with Et below ET_MIN after background subtraction
cfca      call vzero(flag,ncell)
      do i=1, ncell
         flag(i)=0
         if(etc(i)-et_ave .le. et_min) flag(i)=2
      enddo
      njet = 0
c*-Initiator cell is the one with highest Et of not yet used ones
      i=1
      j=index(i)
      if(i.eq.1. and. etc(j).lt.et_seed) then
        if(mode.eq.0 .or. (mode.eq.1 .and. n_iter.eq.1)) then
          call hf1(idPerfomance, 2., 1.)
          print *,' no cell with Et higher than et_seed ->', etc(j)       
        endif
      endif 
      do while(etc(j) .ge. et_seed)
         if(flag(j) .eq. 0) then
            et =etc(j)-et_ave
            eta=etac(j)
            phi=phic(j)
            eta0=eta
            phi0=phi
            etab=eta
            phib=phi
            ets =0.0
            etas=0.0
            phis=0.0
c*-weighted eta and phi. 
            do k= 1, ncell
               l=index(k)
               if(flag(l).eq.0) then
                  deta=etac(l)-eta
c*-Is the cell is in the cone?
                  if(abs(deta).le.cone_rad)then
                     dphi=phic(l)-phi
                     do while(dphi .gt. C_PI)
                        dphi=dphi-C_2PI
                     enddo
                     do while(dphi .le. -C_PI)
                        dphi=dphi+C_2PI
                     enddo
                     if(abs(dphi).le.cone_rad) then
                        r_int=sqrt(deta**2+dphi**2)
                        if(r_int.le.cone_rad)then
c*-calculate offset from initiate cell
                           deta=etac(l)-eta0
                           dphi=phic(l)-phi0
                           do while(dphi .gt. C_PI)
                              dphi=dphi-C_2PI
                           enddo
                           do while(dphi .lt. -C_PI)
                              dphi=dphi+C_2PI
                           enddo
                           et=etc(l)-et_ave
                           etas=etas+abs(et)*deta
                           phis=phis+abs(et)*dphi
                           ets =ets +et
c*-New weighted eta and phi including this cell
                           eta=eta0+etas/ets
                           phi=phi0+phis/ets                          
c*-If cone does not move much from previous cone, just go next step
                           r_int=sqrt((eta-etab)**2+(phi-phib)**2)
                           if(r_int .le. min_move) then
                              goto 159
                           endif
c*-Cone should not move more than MAX_CONE_MOVE from initiator cell
                           r_int=sqrt((etas/ets)**2+(phis/ets)**2)              
                           if(r_int .ge. max_move) then
                              eta=etab
                              phi=phib
                              goto 159
                           endif
c*-Store this loop information
                           etab=eta
                           phib=phi
                        endif
                     endif
                  endif
               endif
            enddo
 159        continue 
            
c*-sum up unused cells within required distance of given eta/phi
            nc=0
            ets=0.0
            etas=0.0
            phis=0.0
            do k=1,ncell
               l=index(k)
               if(flag(l) .eq. 0) then
                  deta=etac(l)-eta
                  if(abs(deta).le.cone_rad)then
                     dphi=phic(l)-phi
                     do while(dphi .gt. C_PI)
                        dphi=dphi-C_2PI
                     enddo
                     do while(dphi .le. -C_PI)
                        dphi=dphi+C_2PI
                     enddo
                     if(abs(dphi).le.cone_rad) then
                        r_int=sqrt(deta**2+dphi**2)
                        if(r_int.le.cone_rad)then
                           flag(l)=-1
                           et  =etc(l)-et_ave
                           ets =ets +et
                           etas=etas+et*deta
                           phis=phis+et*dphi
                           nc  = nc + 1
                        endif
                     endif
                  endif
               endif
            enddo  ! do k=1,ncell
!  5-oct-2001 by PAI - remove 20-feb-2002 by PAI
! 20-feb-2002 - it is work if you apply cut on eT before jet finder !!!
!            if(maxTowerInJet .gt. nc) then
!              ets = ets - et_ave*(maxTowerInJet - nc) 
!            endif
! 5-oct-2001 by PAI
            
c*-reject cluster below minimum Ej_min
c* protection (am)
            etas=eta+etas/ets
            arg = 0.
            if (ets .ne. 0.) then
               if (abs(etas/ets) .lt. 23.719) then
                  arg = ets * cosh(etas/ets)
               else
                  arg = 1.e10
               endif
            endif
            
            if(arg .lt. ej_min) then
               do k=1,ncell
                  if(flag(k).le.0) flag(k)=0
               enddo
               if(njet.eq.0) call hf1(idPerfomance, 3., 1.)
            else
c*-eles, store flags and jet variables
               do k=1,ncell
                  if(flag(k).eq.-1) flag(k)=1
               enddo
               phi=phi+phis/ets
               do while(phi .ge. C_2PI)
                  phi=phi-C_2PI
               enddo
               do while(phi .lt. 0.0)
                  phi=phi+C_2PI
               enddo
               njet=njet+1
               etj(njet) =ets
               etaj(njet,1)=eta0
               phij(njet,1)=phi0
               etaj(njet,2)=etas
               phij(njet,2)=phi
               ncellj(njet)=nc
               call hf1(112, float(nc)/float(maxTowerInJet), 1.) ! 8-oct-2001
            endif 
         endif
         i=i+1
         j=index(i)        
      enddo
      
c*-recalculate energy sum excluding used cells.
      if(mode.eq.1)then
         et_sum=0.0
         nc=0       ! #cells in jets
         do i=1,ncell
            if(flag(i).ne.1) then  ! 1 if cell in jet
               et_sum=et_sum+etc(i)
            else
               nc=nc+1
            endif
         enddo
c*-if background level changes more than prec_bg, go next iteration!!!
c*-after 10 iteration, stop working and finish
         if( (abs(et_sum-et_sum_old)/et_sum.gt.prec_bg 
     +        .and. n_iter.le.10)
     +        .or. n_iter.eq.1) then ! minimum 2 iteration - 10-oct-2001 by pai
            et_ave=et_sum/float(ncell_tot-nc)
            n_iter=n_iter+1
            et_sum_old = et_sum
            print *,'End of iteration : et_ave ', et_ave, ' nc ', nc
     + , ' Jet energy ', ets     
            goto 999
c*-Watch out!!! Here is a big jump!!! 
         endif
         occupationInJet = float(ncell) / float(ncell_tot)
         call hf1(111, occupationInJet, 1.)
         write(*,*) njet,' jet found in ',n_iter,
     +     ' iteration(s) : EtSum, EtAve =',
     +     et_sum,et_sum/float(ncell_tot-nc)
     + , ' ncell_tot ', ncell_tot, ' #cell in jet ', nc
     + , ' occupationAll ', occupationAll
     + , ' occupationInJet ', occupationInJet
      endif
      
      if(njet.gt.100)then
         write(*,*)'UA1:Problem:More than 100 jets found!'
         ierror = 1
      elseif(njet.eq.0)then
         write(*,*)'UA1:Done:No jet found!'
      else
         write(*,*)
     +'UA1:Done:Found ',njet,' jet(s)'         
      end if
      return
      end
Commit	Line	Data
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
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
2398fd93	37	integer idPerfomance
2398fd93	38	data idPerfomance /119/
471f69dc	39	save idPerfomance
471f69dc	40	! print parameter - 27-sep-2001 by PAI
2398fd93	41	integer kpri
2398fd93	42	data kpri /0/
471f69dc	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_PIcone_radcone_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
ae70b416	107	print *,'Iter ', n_iter, ' et_ave ', et_ave, ' #cells ', ncell
471f69dc	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
b9d0a01d	111	cfca call vzero(flag,ncell)
471f69dc	112	do i=1, ncell
b9d0a01d	113	flag(i)=0
471f69dc	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(deta2+dphi2)
	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(deta2+dphi2)
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
ae70b416	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
471f69dc	229	! 5-oct-2001 by PAI
	230
	231	c*-reject cluster below minimum Ej_min
1a434afd	232	c* protection (am)
471f69dc	233	etas=eta+etas/ets
1a434afd	234	arg = 0.
a13dfde0	235	if (ets .ne. 0.) then
1a434afd	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
471f69dc	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
ae70b416	292	et_sum_old = et_sum
	293	print *,'End of iteration : et_ave ', et_ave, ' nc ', nc
	294	+ , ' Jet energy ', ets
471f69dc	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