Protection against floating point exception in cosh.

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

diff --git a/EMCAL/jet_finder_ua1.F b/EMCAL/jet_finder_ua1.F
index 8f293cf42327eb43609c38a96bbd1ee5a43c8324..4cdc1cdeda9ad8a096bf2b2d9aee2765a3d6c655 100644 (file)
--- a/EMCAL/jet_finder_ua1.F
+++ b/EMCAL/jet_finder_ua1.F
@@ -18,6 +18,7 @@ c#include "calorSize.inc"
       real C_2PI 
       real etaCellSize 
       real phiCellSize 
       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
       INTEGER NMAX,JMAX
       parameter(NMAX=30000,JMAX=100) ! 10-oct-2201
       integer ncell, ierror, mode, ncell_tot


@@ -101,11 +102,13 @@ c*-sum up total energy of all cells
       else
          et_ave=0.0
       endif
       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
 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

-      call vzero(flag,ncell)
+cfca      call vzero(flag,ncell)

       do i=1, ncell
       do i=1, ncell

+         flag(i)=0

          if(etc(i)-et_ave .le. et_min) flag(i)=2
       enddo
       njet = 0
          if(etc(i)-et_ave .le. et_min) flag(i)=2
       enddo
       njet = 0


@@ -216,15 +219,26 @@ c*-sum up unused cells within required distance of given eta/phi
                   endif
                endif
             enddo  ! do k=1,ncell
                   endif
                endif
             enddo  ! do k=1,ncell

-! 5-oct-2001 by PAI
-            if(maxTowerInJet .gt. nc) then
-              ets = ets - et_ave*(maxTowerInJet - nc) 
-            endif
+!  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
 ! 5-oct-2001 by PAI
             
 c*-reject cluster below minimum Ej_min

+c* protection (am)

             etas=eta+etas/ets
             etas=eta+etas/ets

-            if(ets*cosh(etas/ets).lt.ej_min) then
+            arg = 0.
+            if (arg .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
                do k=1,ncell
                   if(flag(k).le.0) flag(k)=0
                enddo


@@ -273,7 +287,9 @@ c*-after 10 iteration, stop working and finish
      +        .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
      +        .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
+            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
             goto 999
 c*-Watch out!!! Here is a big jump!!! 
          endif