#---------------------------------------------------------------------------
.help dst_remove_particles May92 hrs
.ih
NAME
dst_remove_particles -- Remove particle hits, record events, find averages.
.endhelp
#---------------------------------------------------------------------------

procedure dst_remove_particles (data, limit, n_diodes, grp_name, time,
                                event_table, event_row, sum,
                                sumsq, n_adds, n_particles, n_event_diodes,
                                total_event_rate, maximum, mean, sigma)

double  data[n_diodes]          # I:  The current group data.
double  limit[n_diodes]         # I:  Particle hit limits.
int     n_diodes                # I:  Number of diodes.
char    grp_name[ARB]           # I:  Name of current group.
char    time[ARB]               # I:  Time of exposure.
pointer event_table             # I:  Events table descriptor.
int     event_row               # IO: Current row of events table.
double  sum[n_diodes]           # O:  Diodes not affected by particle hits.
double  sumsq[n_diodes]         # O:  Diodes squared.
int     n_adds[n_diodes]        # O:  Diode not affected by particle hit.
int     n_particles             # O:  Number of events found.
int     n_event_diodes          # O:  Total number of diodes in all events.
double  total_event_rate        # O:  Total count rate for all events.
double  maximum                 # O:  The maximum for the useful data.
double  mean                    # O:  Mean of current group w/o particle hits.
double  sigma                   # O:  Sigma of current group w/o particle hits.

# Declarations.
double  event_counts            # Total count rate in current event.
double  n_data                  # Number of good diodes.

int     diode                   # Current diode.
int     first_diode             # First diode of an event.
int     n_diodes_per_event      # Number of diodes in current event.

bool    event                   # True if we are in an event.

pointer sp                      # Stack pointer.
pointer tmp_string              # Temporary string.

# Function prototypes.
double  asumd()
real    asumi()

begin
        call smark (sp)
        call salloc (tmp_string, SZ_LINE, TY_CHAR)
        
        # Loop through the diodes looking for particle events.
        event = false
        maximum = 0.0d0
        n_event_diodes = 0
        total_event_rate = 0.0d0
        n_particles = 0
        do diode = 1, n_diodes {
            
            # Check if current diode is a hit.
            if ((!IS_INDEFD(limit[diode])) && (data[diode] > limit[diode])) {
                
                sum[diode] = 0.0d0
                n_adds[diode] = 0
                
                # Check for new event.
                if (!event) {
                    event = true
                    first_diode = diode
                    n_diodes_per_event = 1
                    event_counts = data[diode]
                    n_particles = n_particles + 1
                } else {
                    n_diodes_per_event = n_diodes_per_event + 1
                    event_counts = event_counts + data[diode]
                }
            }
            
            # Not an event.
            else {
                
                # Indicate useful diode.
                sum[diode] = data[diode]
                n_adds[diode] = 1
                maximum = max (maximum, data[diode])
                
                # If there had been an event- write it out to the
                # event table.
                if (event) {
                    call dst_event_write (event_table, grp_name, first_diode,
                                          n_diodes_per_event, event_counts,
                                          time, event_row)
                    
                    n_event_diodes = n_event_diodes + n_diodes_per_event
                    total_event_rate = total_event_rate + event_counts
                    event = false
                }
            }
        }
        
        # Now figure out group statistics and write them out.
        n_data = double (asumi (n_adds, n_diodes))
        if (n_data > 0) {
            mean = asumd (sum, n_diodes) / n_data
            call apowkd (sum, 2, sumsq, n_diodes)
            sigma = sqrt ((asumd (sumsq, n_diodes) / n_data) -
                          (mean * mean))
        } else {
            call sprintf (Memc[tmp_string], SZ_LINE,
                          "No good diodes in grp %s.\n")
            call pargstr (grp_name)
            call eprintf (Memc[tmp_string])
        }
        
        call sfree (sp)
end
#---------------------------------------------------------------------------
# End of dst_remove_particles
#---------------------------------------------------------------------------