subroutine hec_his_ana1 
c
      implicit none
c
      include'hec_par.inc'      ! shared parameters
      include'hec_datacard.inc' ! datacard values
      include'hec_evth.inc'     ! event header common
      include'hec_geo.inc'      ! geometry correspondence tables
      include'hec_his_sys.inc'  ! histogram package
      include'hec_hvpar.inc'    !
      include'hec_imp.inc'      ! impact point common
      include'hec_ped.inc'      ! pedestal common
      include'hec_ped_sys.inc'  ! pedestal system common
      include'hec_runh.inc'     ! run header common
      include'hec_slow_pc.inc'  !
      include'hec_tdc.inc'      ! tdc common
c
c local variables
c
      character*4 ct
      character*2 cpart
      integer ic, i_t, nid, nid2, id, ibin
      integer i_z, i_mod, i_eta, i_phi
      real signal_sum(i_z_max,i_mod_max)
      real avetdc(ntdc_bin), errtdc(ntdc_bin)
      real error(i_adc_dim)
      real weight, average
      integer i,j
      real calv1(nchhv_cal_max),pmv1(nchhv_pm_max)
      real mwpcv1(nchhv_mwpc_max)
      real calv2(nchhv_cal_max),pmv2(nchhv_pm_max)
      real mwpcv2(nchhv_mwpc_max)
c
c fill adc pedestal and adc pedestal rms histograms
c
      nid = 0
      nid2 = 0
      call vzero(signal_sum,i_z_max*i_mod_max)
      call vzero(average_max,i_adc_dim)
      call vzero(errmax,i_adc_dim)
c
      do ic = 1, i_adc_used
c
c histo type 1 -- ped/noise/scaler
c
        if (histo_switch(1) .eq. 1) then
          call hf1(id_ped, adc_ped_run(ic), 1.)
          call hf1(id_ped_ic, real(ic), adc_ped_run(ic))
          call hf1(id_rms, adc_rms_run(ic), 1.)
c
c fill error array for rms vs. ic plot
c
          if (i_adc0_1(ic) .ge. 1) then
            error(ic) = adc_rms_run(ic)/sqrt(2.*real(i_adc0_1(ic)))
          else
            error(ic) = 0.
          endif
        endif
c
c skip disconnected channels
c
        if(ic_connected(ic))then
c
c histo type 2 -- geometric energy distribution
c
        if (histo_switch(2) .eq. 1) then
c
c fill arrays for histo of max vs. channel
c
          if (adc_rms_run(ic).lt.25. .and. maxmom(ic,1).ge.1.d0) then
            average_max(ic) = real( maxmom(ic,2)/maxmom(ic,1) )
            errmax(ic) = real( maxmom(ic,3)/maxmom(ic,1) - 
     +        dble(average_max(ic)**2) )
            if (errmax(ic) .gt. 0.) then
              errmax(ic) = real(dsqrt(dble(errmax(ic))/maxmom(ic,1)))
            else
              errmax(ic) = 0.
            endif
          else
            print *,'ic of bad cell is ',ic
            average_max(ic) = 0.
            errmax(ic) = 0.
          endif
c
c fill histo with lego map of weighted "hits" in each layer
c
          i_eta = ictoph(ic,1)
          i_phi = ictoph(ic,2)
          i_z = ictoph(ic,3)
          if (i_eta*i_phi*i_z .ne. 0) then
            if (maxmom(ic,1) .ge. 1.d0) then
              weight = real( maxmom(ic,2) / maxmom(ic,1) )
            else
              weight = 0.
            endif
            if (i_z .ge. 1 .and. i_z .le. i_z_used) then
              id = id_map - 1 + i_z
c              if (i_eta .le. 6.) then
              call hfill(id, real(i_phi), real(i_eta), weight)
c              endif
            endif
          endif
c
c sum to get average content of each depth for each module
c (includes bad cell cut)
c
          i_mod = ictodb(ic,1)
          if (i_z .ne. 0 .and. i_mod .ne. 0) then
            signal_sum(i_z,i_mod) = signal_sum(i_z,i_mod) + 
     +        average_max(ic)
          endif
c
        endif
c
c fill time profile histograms
c
c
c histo type 3 -- no tdc or calib corrections 
c
        if (histo_switch(3) .eq. 1) then
          do i_t = 1, i_tim_used
c
            if (adcmom(ic,i_t,1) .ge. 1.d0) then
              average = real( adcmom(ic,i_t,2) / adcmom(ic,i_t,1) )
            else
              average = 0.
            endif
c
            id = id_tprof_ic + ic
            call hff1(id, nid, real(i_t), average)
          enddo
        endif
c
c histo type 4 -- tdc corrected 
c
        if (histo_switch(4) .eq. 1) then
c
c errtdc is actually rms/rootn (error on the mean)
c
          do ibin = 1, ntdc_bin
            if (tdcmom(ic,ibin,1) .ne. 0.d0) then
              avetdc(ibin) = real( tdcmom(ic,ibin,2)/tdcmom(ic,ibin,1) )
              errtdc(ibin) = real( tdcmom(ic,ibin,3)/tdcmom(ic,ibin,1) 
     +          - dble(avetdc(ibin))**2 )
              if (errtdc(ibin) .gt. 0.) then
                errtdc(ibin) =sqrt(dble(errtdc(ibin))/tdcmom(ic,ibin,1))
              else
                errtdc(ibin) = 0.
              endif
            else
              avetdc(ibin) = 0.
              errtdc(ibin) = 0.
            endif
          enddo
          id = id_timns_ic + ic
          call hpak(id,avetdc)    
          call hpake(id,errtdc)
c
        endif
c
        endif
c
c end loop over channels
c
      enddo
c
c histo type 1 -- ped/noise 
c
      if (histo_switch(1) .eq. 1) then
c
c fill histogram of rms vs. channel number
c error is (bin content)/sqrt(2n)
c
        call hpak(id_rms_ic,adc_rms_run)    
        call hpake(id_rms_ic,error)
        call hpak(id_pedevt+1, adc_ped_time_ic1)
        call hpak(id_pedevt+2, adc_ped_time_ic2)
        call hpak(id_pedevt+3, adc_ped_time_ic3)
        call hpak(id_pedevt+4, adc_ped_time_ic4)
      endif
c
c histo type 2 -- geometric energy distribution 
c
      if (histo_switch(2) .eq. 1) then
c
c fill histogram of maximum adc vs. channel number
c
        call hpak(id_max_ic, average_max)    
        call hpake(id_max_ic, errmax)
c
c fill depth profile histos
c
        do i_mod = 1, i_mod_used
          do i_z = 1, i_z_used
            id = id_adc_z + i_z_used + i_mod
            call hfill(id,i_z-0.5,0.,signal_sum(i_z,i_mod))
          enddo
        enddo
c
      endif
c
c slow control histos      
c 
      if(irunpd.gt.7 .and. histo_switch(7).eq.1)then
         do i = 0, nmodhv_cal_max - 1
            do j=1,nchhv_cal_max
               calv1(j)=hv_cal_hv(i,j)
               calv2(j)=hv_cal_i(i,j)
            enddo
            call hpak(6600+i,calv1)
            call hpak(6610+i,calv2)
         enddo
         do i = 0, nmodhv_pm_max - 1
            do j=1,nchhv_pm_max
               pmv1(j)=hv_pm_hv(i,j)
               pmv2(j)=hv_pm_i(i,j)
            enddo
            call hpak(6620+i,pmv1)
            call hpak(6630+i,pmv2)
         enddo
         do i = 0, nmodhv_mwpc_max - 1
            do j=1,nchhv_mwpc_max
               mwpcv1(j)=hv_mwpc_hv(i,j)
               mwpcv2(j)=hv_mwpc_i(i,j)
            enddo
            call hpak(6640+i,mwpcv1)
            call hpak(6650+i,mwpcv2)
         enddo
         call hpak(6660,cryo_temp,0.,1.)
         call hfill(6661,cryo_press,0.,1.)
      endif
c
c energy cluster for impact point when appropriate: a mere copy
c of the corresponding cluster
c
      if(i_part.eq.1)then
         cpart='e-'
      else if(i_part.eq.2)then
         cpart='mu'
      else if(i_part.eq.3)then
         cpart='pi'
      else if(i_part.eq.0)then
         cpart='xx'
      endif
      j=0
      if(impact_point_index.ne.0)then
         do i=1,his_clus_num
            ct=his_clus_title(his_clus_index(i))
            if(impact_point_name(impact_point_index).eq.ct(1:1))then
               if(cpart.eq.ct(3:4))then
                  call hcopy(id_clus+i,9091,' ')
                  call hcopy(id_clus_cal+i,9092,' ')
                  call hcopy(id_clus_cal_hv+i,9093,' ')
               else if(cpart.eq.'xx')then
                  j=j+1
                  call hcopy(id_clus+i,9090+j,' ')
               endif
            endif
         enddo
      endif
c
c Also copy the Z clusters histograms (adc, calib, hv+calib) 
c for the actual particle type to ids 19091,19092,19092.
c In case of noise run they will contain adc distributions
c for the 3 kinds of particle. These copies are used to ease
c the online display/fitting and writing to files
c
      j=0
      do i=1,his_clus_num
         ct=his_clus_title(his_clus_index(i))
         if(ct(1:1).eq.'Z')then
            if(cpart.eq.ct(3:4))then
               call hcopy(id_clus+i,19091,' ')
               call hcopy(id_clus_cal+i,19092,' ')
               call hcopy(id_clus_cal_hv+i,19093,' ')
            else if(cpart.eq.'xx')then
               j=j+1
               call hcopy(id_clus+i,19090+j,' ')
            endif
         endif
      enddo
c
      return
      end