Return to h_clusters_cal.f CVS log

Up to [HallC] / Analyzer / HTRACKING

Code for new Shower counter tubes

*=======================================================================
      subroutine h_clusters_cal(abort,errmsg)
*=======================================================================
*-
*-      Purpose: Finds clusters in the calorimeter and computes
*-               for each cluster it's size(number of hit blocks),
*-               position, energy deposition in the calorimeter
*-               columns and the total energy deposition.
*-               The energy depositions are not corrected yet for
*-               impact point coordinate dependence.
*-               A cluster is defined as a set of adjacent hit blocks
*-               which share a common edge or a corner. Any two hits
*-               from different clusters are separated by at least one 
*-               block which has not fired.
*-
*-      Input Banks: HMS_SPARSIFIED_CAL, HMS_DECODED_CAL
*-
*-      Output Bank: HMS_CLUSTERS_CAL
*-
*-      Created: 15 Mar 1994      Tsolak A. Amatuni
*-      Modified 25 Mar 1994      DFG
*-                                Change name of print routine
*                10 Apr 1994      DFG Protect for Et=0 division
* $Log: h_clusters_cal.f,v $
* Revision 1.5  1999/02/03 21:13:23  saw
* Code for new Shower counter tubes
*
* Revision 1.4  1995/05/22 19:39:07  cdaq
* (SAW) Split gen_data_data_structures into gen, hms, sos, and coin parts"
*
* Revision 1.3  1994/08/02  20:00:12  cdaq
* (JRA) Catch some out of bounds problems
*
* Revision 1.2  1994/04/13  05:31:37  cdaq
* *** empty log message ***
*
* Revision 1.1  1994/04/12  21:30:02  cdaq
* Initial revision
*
*-
*-----------------------------------------------------------------------
*
*
      implicit none
      save
*
      logical abort
      character*(*) errmsg
      character*14 here
      parameter (here='H_CLUSTERS_CAL')
*
      integer*4 nc                     !Cluster number
      integer*4 ihit,jhit,khit,nh      !Internal loop counters.
      integer*4 hits_tagged            !Current number of tagged hits.
      integer*4 irow,icol,jrow,jcol,col!Row and column indecies
      integer*4 d_row,d_col            !Distance between rows(columns)
      logical tagged
*
      include 'hms_data_structures.cmn'
      include 'hms_calorimeter.cmn'
*
*
      hnclusters_cal=0
      if(hcal_num_hits.le.0) go to 100   !Return
*
      do ihit=1,hmax_cal_blocks
        hcluster_hit(ihit)=0
      enddo
*
      nc                 = 1
      hcluster_hit(1)    =-1
      hcluster_size(1)   = 1
      hits_tagged        = 1
      tagged             =.true.
*
*     Find the clusters.
*
*
*-----Loop untill all the hits are tagged
      do while(hits_tagged.le.hcal_num_hits)
*
*-------Loop untill there are no more hits
*-------in the current cluster to be tagged
        do while(tagged)
          tagged=.false.
*
*---------Loop over all the hits
          do ihit=1,hcal_num_hits
*     
*-----------and find a hit("seed") which belongs to the
*-----------current cluster, but it's neighbors are not tagged
            if(hcluster_hit(ihit).lt.0) then
              irow  =hcal_rows(ihit)
              icol  =hcal_cols(ihit)
*
*-------------Loop over all the hits
              do jhit=1,hcal_num_hits
*
*---------------and find hits which are not tagged yet
                if(hcluster_hit(jhit).eq.0) then
                  jrow =hcal_rows(jhit)
                  jcol =hcal_cols(jhit)
                  d_row=iabs(jrow-irow)
                  d_col=iabs(jcol-icol)
*
*-----------------Are these hits a neighbor to "seed"?
                  if(d_row.le.1.and.d_col.le.1) then
*
*-------------------Assign them to the same current cluster
                    hcluster_hit(jhit)=hcluster_hit(ihit)
                    hcluster_size(nc) =hcluster_size(nc)+1
                    hits_tagged       =hits_tagged+1
                    tagged            =.true.
*
                  endif                 !End ... if neighbor of "seed"
*
                endif                   !End ... if not scanned yet
*
              enddo                     !End loop over all hits
*
*-------------All the neighbors of "seed" were scanned
              hcluster_hit(ihit)=-hcluster_hit(ihit)
*
            endif                       !End ... if "seed"
*
          enddo                         !End loop over all hits
*
        enddo                           !All the hits of the current cluster were tagged
*
*-------Initialize to start the search for the next cluster
        nc            =nc+1
        hits_tagged   =hits_tagged+1
        tagged        =.true.
*     
*-------Find a hit which is not tagged
        khit=1
        do while(hcluster_hit(khit).ne.0 .AND. KHIT.LT.HMAX_CAL_BLOCKS)
          khit=khit+1
        enddo
*
*-------This will be the new "seed"
	IF (NC.GT.HNCLUSTERS_MAX) NC=HNCLUSTERS_MAX !AVOID OUT/BOUNDS.
        hcluster_hit(khit)=-nc
        hcluster_size(nc) = 1
* 
      enddo                             !End. Now all the hits are assigned to some cluster
*
*-----Number of clusters found
      hnclusters_cal=nc-1
*
*     For each cluster found, compute the center of gravity in X
*     projection, the energy deposited in succesive calorimeter columns
*     and the total energy deposition
*     
      do nc=1,hnclusters_max
        hcluster_e1_pos(nc)=0.      
        hcluster_e1_neg(nc)=0.
        hcluster_e2_pos(nc)=0.      
        hcluster_e2_neg(nc)=0.     
*
        hcluster_e1(nc)=0.
        hcluster_e2(nc)=0.
        hcluster_e3(nc)=0.
        hcluster_e4(nc)=0.
        hcluster_et(nc)=0.
        hcluster_xc(nc)=0.
      enddo
*
*
      do nh=1,hcal_num_hits
        nc = MAX(1,hcluster_hit(nh) )   !THIS DOES NOT HELP THE ANALYSIS,
                                        !BUT IT AVOIDS SUBSCRIPT OUT/BOUNDS ERRS.
        col=hcal_cols(nh)
*     
        hcluster_xc(nc)=hcluster_xc(nc)+hblock_xc(nh)*hblock_de(nh)
*
        if(col.eq.1) then
          if(hcal_num_neg_columns.ge.1) then
            hcluster_e1_pos(nc)=hcluster_e1_pos(nc)+hblock_de_pos(nh)
            hcluster_e1_neg(nc)=hcluster_e1_neg(nc)+hblock_de_neg(nh)
            hcluster_e1(nc)=hcluster_e1_pos(nc)+hcluster_e1_neg(nc)
          else
            hcluster_e1(nc)=hcluster_e1(nc)+hblock_de(nh)
          endif
        else if (col.eq.2) then
          if(hcal_num_neg_columns.ge.2) then
            hcluster_e2_pos(nc)=hcluster_e2_pos(nc)+hblock_de_pos(nh)
            hcluster_e2_neg(nc)=hcluster_e2_neg(nc)+hblock_de_neg(nh)
            hcluster_e2(nc)=hcluster_e2_pos(nc)+hcluster_e2_neg(nc)
          else
            hcluster_e2(nc)=hcluster_e2(nc)+hblock_de(nh)
          endif
        else if(col.eq.3) then
          hcluster_e3(nc)=hcluster_e3(nc)+hblock_de(nh)
        else if(col.eq.4) then
          hcluster_e4(nc)=hcluster_e4(nc)+hblock_de(nh)
        endif
        hcluster_et(nc)=hcluster_et(nc)+hblock_de(nh)  ! Is hblock_de de_pos+de_neg?
*
      enddo
*
      do nc=1,hnclusters_cal
*     make sure hcluster_et .ne. zero so no divide by zero        
        if(hcluster_et(nc).gt.0.) then
          hcluster_xc(nc)=hcluster_xc(nc)/hcluster_et(nc)
        else
          hcluster_xc(nc)= -1.0         ! Set fraction negative for bad et
        endif
      enddo
*
 100  continue
      if(hdbg_clusters_cal.gt.0) call h_prt_cal_clusters
*
      return
      end