# Define memory management for CRATE array.
define CRATE    Memd[crate+$1-1]

#---------------------------------------------------------------------------
.help dst_find_noisy May92 hrs
.ih
NAME
dst_find_noisy -- Compute summary of all observations and find noisy diodes.
.endhelp
#---------------------------------------------------------------------------

procedure dst_find_noisy( root, sum, n_data, n_diodes, dthresh,
                          dnsig, verbose, noisy_table, noisy_row )

char    root[ARB]               # I:  Rootname of observations.
double  sum[n_diodes]           # I:  Sum of all observations per diode.
int     n_data[n_diodes]        # I:  Number of observations summed per diode.
int     n_diodes                # I:  Number of diodes.
double  dthresh                 # I:  Threshold for noisy diodes.
double  dnsig                   # I:  Number of sigma for noisy diode.
bool    verbose                 # I:  Write summary to standard output.
pointer noisy_table             # I:  Table of noisy diodes.
int     noisy_row               # IO: Row of noisy table to write to.

# Declarations.
double  mean                    # Mean
double  new_sum, new_sumsq      # New sums.
double  sigma                   # Sigma
double  thresh                  # Threshold level to determine noisy diodes.

int     diode                   # Current diode.
int     n_new                   # Number of new values.

pointer sp                      # Stack pointer.
pointer crate                   # Average rate for each diode.

begin
        call smark (sp)
        call salloc (crate, n_diodes, TY_DOUBLE)
        
        # Find noisy diodes.  First find average for each diode, then average
        # all the averages.
        n_new = 0
        new_sum = 0.0d0
        new_sumsq = 0.0d0
        do diode = 1, n_diodes {
            if (n_data[diode] == 0) {
                call eprintf ("No dark count can be computed for diode %d.\n")
                call pargi (diode)
                CRATE(diode) = INDEFD
            } else{
                CRATE(diode) = sum[diode] / n_data[diode]
                n_new = n_new + 1
                new_sum = new_sum + CRATE(diode)
                new_sumsq = new_sumsq + (CRATE(diode) * CRATE(diode))
            }
        }
        mean = new_sum / n_new
        sigma = sqrt ((new_sumsq / n_new) - (mean * mean))
        
        if (verbose) {
            call printf ("For observation group %s:\n")
            call pargstr (root)
            call printf ("Average dark count (all diodes) = %6.6g, Standard deviation = %6.6g\n")
            call pargd (mean)
            call pargd (sigma)
        }
        
        if (IS_INDEFD (dthresh)) {
            if (IS_INDEFD (dnsig)) {
                if (verbose)
                    call printf ("No noisy diodes found.\n")
                return
            } else
                thresh = mean + sigma * dnsig
        } else
            thresh = dthresh
        
        # Now found the noisy ones.
        n_new = 0
        new_sum = 0.0d0
        new_sumsq = 0.0d0
        do diode = 1, n_diodes {
            if (CRATE(diode) > thresh)
                call dst_noisy_write (noisy_table, root, diode, CRATE(diode),
                                      mean, sigma, noisy_row)
            else {
                n_new = n_new + 1
                new_sum = new_sum + CRATE(diode)
                new_sumsq = new_sumsq + (CRATE(diode) * CRATE(diode))
            }
        }
        
        if (n_new == 0)
            call eprintf ("All diodes are noisy, dthresh and dnsig may be too low!\n")
        else if (verbose) {
            if (n_new == n_diodes)
                call printf ("No noisy diodes.\n")
            else {
                call printf ("There were %d noisy diodes found.\n")
                call pargi (n_diodes - n_new)
                mean = new_sum / n_new
                sigma = sqrt ((new_sumsq / n_new) - (mean * mean))
                call printf ("Average dark rate (good diodes) = %6.6g, standard deviation = %6.6g\n")
                call pargd (mean)
                call pargd (sigma)
            }
        }
        
        call sfree(sp)
        
end
#---------------------------------------------------------------------------
# End of dst_find_noisy
#---------------------------------------------------------------------------