subroutine hec_crl_correl ( job )
c ----------------------------------------------------------------------
c-
c-    toy program to calculate
c-
c-      for each adc channel and for each time sampling:
c-      - average signal and r.m.s. (pedestal and noise)
c-
c-      for each adc channel:
c-      - correlation coefficients between different time samplings
c-
c-      for each time sampling:
c-      - correlation coefficients between different adc channels
c-
c-    input:  job - steering flag:
c-                  <0 - initialization
c-                  =0 - analysis of an event
c-                  >0 - final anlysis after all events
c-
c ----------------------------------------------------------------------
c
      implicit none
c
      integer job
c
c ****  include necessary parameters and arrays (in this order!!!)
c
      include 'hec_par.inc'
      include 'hec_adc.inc'
      include 'hec_geo.inc'
      include 'hec_crl.inc'
c
c ****  arrays to calculate average, r.m.s. and correlations
c
      double precision dpadc(crl_maxadc)
      double precision dptim(crl_maxsamp)
c
c ****  miscellaneous
c
      integer it, it2, ia, ia2
      double precision dpsum, dp1,dp2
c
      integer nchan, nsamp
c
c ----------------------------------------------------------------------
c
c ****  initialization
c
      if (job.lt.0) then
c
        nevtcrl=0
c
        do ia=1,crl_maxadc
          do it=1,crl_maxsamp
            dps1(it,ia)=0.d0
            dps2(it,ia)=0.d0
c
            do it2=1,crl_maxsamp
              dpsc1(it,it2,ia)=0.d0
            enddo
c
            do ia2=1,crl_maxadc
              dpsc2(ia,ia2,it)=0.d0
            enddo
          enddo
        enddo
c
c ****  analysis of the current event
c
      elseif (job.eq.0) then
c
c set number of channels/samplings to look at
c
        nchan = min(i_adc_used,crl_maxadc)
        nsamp = min(i_tim_used,crl_maxsamp)
c
        nevtcrl=nevtcrl+1
c
c ****  store data to calculate average and r.m.s.
c
        do ia=1,nchan !  loop over adc channels
          do it=1,nsamp
            dptim(it)=dble(i_adc_ic_t(ia,it))
            dps1(it,ia)=dps1(it,ia) + dptim(it)
            dps2(it,ia)=dps2(it,ia) + dptim(it)*dptim(it)
          enddo
c
c ****  store data to calculate correlations between time samplings
c
          do it=1,nsamp
            do it2=it,nsamp
              dpsc1(it,it2,ia)=dpsc1(it,it2,ia) + dptim(it)*dptim(it2)
            enddo
          enddo
        enddo !  loop over adc channels
c
c ****  store data to calculate correlations between adc channels
c
        do it=1,nsamp !  loop over time samplings
          do ia=1,nchan
            dpadc(ia)=dble(i_adc_ic_t(ia,it))
          enddo
          do ia=1,nchan
            do ia2=ia,nchan
              dpsc2(ia,ia2,it)=dpsc2(ia,ia2,it) + dpadc(ia)*dpadc(ia2)
            enddo
          enddo
        enddo !  loop over time samplings
c
c ****  final analysis: calculate average, r.m.s. and time correlation
c
      else
c
c set number of channels/samplings to look at
c
        nchan = min(i_adc_used,crl_maxadc)
        nsamp = min(i_tim_used,crl_maxsamp)
c
        if (nevtcrl.gt.0) then
          dpsum=dble(nevtcrl)
c
c ****  calculate average and r.m.s.
c
          do ia=1,nchan !  loop over adc channels
            do it=1,nsamp
              dp1 = dps1(it,ia) / dpsum
              dp2 = dps2(it,ia) / dpsum - dp1*dp1
              dp2=dsqrt(dp2)
c
              dpave(it,ia)=dp1
              dprms(it,ia)=dp2
            enddo
c
c ****  calculate correlations between time samplings
c
            do it=1,nsamp
              do it2=it,nsamp
                dp1 = dpsc1(it,it2,ia) / dpsum
                dp1 = dp1 - dpave(it,ia)*dpave(it2,ia)
                dp1 = dp1 /(dprms(it,ia)*dprms(it2,ia))
                dpcortim(it,it2,ia)=dp1
                if (it2.ne.it) then
                  dpcortim(it2,it,ia)=dp1
                endif
              enddo
            enddo
          enddo                 !  loop over adc channels
c
c ****  calculate correlations between adc channels
c
          do it=1,nsamp         !  loop over time samplings
            do ia=1,nchan
              do ia2=ia,nchan
                dp1 = dpsc2(ia,ia2,it) / dpsum
                dp1 = dp1 - dpave(it,ia)*dpave(it,ia2)
                dp1 = dp1 /(dprms(it,ia)*dprms(it,ia2))
                dpcoradc(ia,ia2,it)=dp1
                if (ia2.ne.ia) then
                  dpcoradc(ia2,ia,it)=dp1
                endif
              enddo
            enddo
          enddo                 !  loop over time samplings
c
          write (6,*) ' *** correl: ',nevtcrl,
     &    ' events for analysis *** '
c
        else
          write (6,*) ' *** correl:  no events for analysis *** '
          do ia=1,nchan
            do it=1,nsamp
              dpave(it,ia)=0.d0
              dprms(it,ia)=0.d0
              do it2=1,nsamp
                dpcortim(it,it2,ia)=0.d0
              enddo
              do ia2=1,nchan
                dpcoradc(ia,ia2,it)=0.d0
              enddo
            enddo
          enddo
c
        endif
      endif
c
      return
      end
c ======================================================================