include include # Number of characters of root that should match for similar observations. define TEMPLATE_SIZE 6 # How many diodes that represent background/radiation diodes. define N_BACK_DIODES 4 define N_RAD_DIODES 2 #--------------------------------------------------------------------------- .help dst_darkstats May92 hrs .ih NAME dst_darkstats -- Compute various statistics for GHRS dark observations. .endhelp #--------------------------------------------------------------------------- procedure dst_darkstats() # Declarations double dnum_obs_sum # Sum of isolated diode per observation. double dnum_tot_sum # Sum of isolated diode for all observations. double dnum_value # Mean of the isolated diode. double maximum, mean, sigma # Statistics. double obs_arate # Anti-coincidence rate sum per obs. double obs_brate # Background diode rate sum per obs. double obs_event_rate # Total event count rate per observation. double obs_max # Maximum per observation. double obs_rrate # Radiation diode counts per obs. double grp_arate # Anticoincidence rate per group. double grp_brate # Background diode rate per group. double grp_event_rate # Total event count rate per group. double grp_rrate # Radiation counts per group. double tot_arate # Anticoincidence rate sum for all obs. double tot_brate # Background diode rate sum for all obs. double tot_event_rate # Event count rate per related observations. double tot_max # Maximum per total of observations. double tot_rrate # Radiation counts for all observations. real datamin, datamax # Mininum, maximum value of current group. int detector # Detector from which data came from. # Must always be the same. int event_row # Current row of events table. int group # Current group number. int grp_row # Current row of the group table. int i # Generic. int n_diodes # Number of diodes per observation. Must # always be the same. int n_event_diodes # Number of diodes in events per group. int n_files # Number of input files specified. int n_groups # Number of groups per observation. int n_grps # Total number of groups for current sums. int n_particles # Number of particle hits per group. int noisy_row # Table row for current noisy diode. int obs_event_diodes # Number of diodes in events per observation. int obs_row # Current row of the observation table. int tot_event_diodes # Number of diodes in events per related obs. int tot_particles # Total number of particle hits found. int tot_row # Current row of the total table. bool first # TRUE if first time through image list loop. pointer dark_table # Table descriptor for the dark averages. pointer event_table # Table descriptor for events. pointer extracted_image # Image descriptor for extracted data file. pointer file # Current file opened. pointer first_file # Name of the first file in the related group. pointer list # Observation list structure. pointer n_adds # Number of groups in final sums. pointer noisy_table # Table descriptor for noisy diodes. pointer obs_n_data # Number of values in observation sums. pointer obs_sum # Sum of diodes per observation. pointer obs_sumsq # Sum of square of diodes per observation. pointer obs_table # Table descriptor for observation stats. pointer particle_limit # Thresholds indicating a particle hit. pointer root # Root of the file name. pointer grp_data # Observational data. pointer grp_image # Observation image descriptor. pointer grp_table # Table descriptor for single group stats. pointer sp # Stack pointer. pointer special_image # Image descriptor for special diode file. pointer sum # Sum per diode. pointer sumsq # Sum squared per diode. pointer template # Template observation name. pointer time # Time of observation. pointer tmp_string # Temporary string. pointer tot_n_data # Number of points in total sums. pointer tot_sum # Sum of diodes for all observations. pointer tot_sumsq # Sum of square of diodes for all observations. pointer tot_table # Table descriptor for related obs. stats. include "darkstat_params.com" # Function prototypes. int fnroot(), imgeti(), imtgetim(), imtlen(), strlen(), strncmp() pointer imgl1d(), immap(), imtopen() begin call smark(sp) call salloc (file, SZ_LINE, TY_CHAR) call salloc (first_file, SZ_LINE, TY_CHAR) call salloc (root, SZ_LINE, TY_CHAR) call salloc (time, SZ_LINE, TY_CHAR) call salloc (tmp_string, SZ_LINE, TY_CHAR) # Open the list of observations to get stats on. list = imtopen (input) n_files = imtlen (list) if (n_files == 0) call error (1, "darkstats: No input files specified!") # Open the statistics tables. call dst_stat_open (grp_stats, grp_table, grp_row) call dst_stat_open (obs_stats, obs_table, obs_row) call dst_stat_open (total_stats, tot_table, tot_row) # Open the noisy diode table. call dst_noisy_open (noisy, noisy_table, noisy_row) # Open the event table. call dst_event_open (events, event_table, event_row) # Open the dark table. call dst_dark_open (dark, dark_table) # First pass through the data. Collect statistics. first = true n_grps = 0 while (imtgetim (list, Memc[file], SZ_LINE) != EOF) { # Open the image. iferr (grp_image = immap (Memc[file], READ_ONLY, 0)) { call sprintf (Memc[tmp_string], SZ_LINE, "could not open file %s") call pargstr (Memc[file]) call error (1, Memc[tmp_string]) } # Check dimensionality. if (IM_NDIM(grp_image) != 1) { call sprintf (Memc[tmp_string], SZ_LINE, "image %s has too many dimensions (%d)") call pargstr (Memc[file]) call pargi (IM_NDIM(grp_image)) call error (1, Memc[tmp_string]) } # Initialize memory and values if first time through. if (first) { n_diodes = IM_LEN(grp_image, 1) detector = imgeti (grp_image, "DETECTOR") n_groups = IM_CLSIZE(grp_image) call salloc (particle_limit, n_diodes, TY_DOUBLE) call salloc (sum, n_diodes, TY_DOUBLE) call salloc (sumsq, n_diodes, TY_DOUBLE) call amovkd (0.0d0, Memd[sum], n_diodes) call amovkd (0.0d0, Memd[sumsq], n_diodes) # Allocate arrays that will be used later. call salloc (n_adds, n_diodes, TY_INT) call salloc (obs_sum, n_diodes, TY_DOUBLE) call salloc (obs_sumsq, n_diodes, TY_DOUBLE) call salloc (obs_n_data, n_diodes, TY_INT) call salloc (tot_sum, n_diodes, TY_DOUBLE) call salloc (tot_sumsq, n_diodes, TY_DOUBLE) call salloc (tot_n_data, n_diodes, TY_INT) call salloc (particle_limit, n_diodes, TY_DOUBLE) first = false } # Check to make sure current image is consistent with previous. else { if (IM_LEN(grp_image, 1) != n_diodes) { call sprintf (Memc[tmp_string], SZ_LINE, "image %s has odd data length: %d != %d") call pargstr (Memc[file]) call pargi (IM_LEN(grp_image, 1)) call pargi (n_diodes) call error (1, Memc[tmp_string]) } if (imgeti (grp_image, "DETECTOR") != detector) { call sprintf (Memc[tmp_string], SZ_LINE, "image %s is the wrong detector: %d != %d") call pargstr (Memc[file]) call pargi (imgeti (grp_image, "DETECTOR")) call pargi (detector) } if (IM_CLSIZE(grp_image) != n_groups) { call sprintf (Memc[tmp_string], SZ_LINE, "image %s has wrong number of groups: %d != %d") call pargstr (Memc[file]) call pargi (IM_CLSIZE(grp_image)) call pargi (n_groups) } } # Loop through all the groups. do group = 1, n_groups { # If we are into the next group, reopen the next group. if (group > 1) iferr (call gf_opengr (grp_image, group, datamin, datamax, 0)) { call sprintf (Memc[tmp_string], SZ_LINE, "could not get group %d of file %s") call pargi (group) call pargstr (Memc[file]) } # Accumulate statistics. grp_data = imgl1d (grp_image) call aaddd (Memd[sum], Memd[grp_data], Memd[sum], n_diodes) call apowkd (Memd[grp_data], 2, Memd[grp_data], n_diodes) call aaddd (Memd[sumsq], Memd[grp_data], Memd[sumsq], n_diodes) n_grps = n_grps + 1 } # Close the image. call imunmap (grp_image) } # Compute particle hit threshold for each diode. call dst_particle (Memd[sum], Memd[sumsq], n_grps, pthresh, pnsig, Memd[particle_limit], n_diodes) # Next pass on data. This is where the real stats come from. # Boy, ain't this a mess! template = NULL Memc[first_file] = EOS call amovkd (0.0d0, Memd[tot_sum], n_diodes) call amovkd (0.0d0, Memd[tot_sumsq], n_diodes) call amovki (0, Memi[tot_n_data], n_diodes) tot_max = 0.0d0 dnum_tot_sum = 0.0d0 n_grps = 0 tot_particles = 0 tot_event_diodes = 0 tot_event_rate = 0.0d0 tot_arate = 0.0d0 tot_brate = 0.0d0 tot_rrate = 0.0d0 grp_image = NULL extracted_image = NULL special_image = NULL call imtrew (list) while (imtgetim (list, Memc[file], SZ_LINE) != EOF) { # Store the name of the file if its the first in this group. if (strlen (Memc[first_file]) == 0) call strcpy (Memc[file], Memc[first_file], SZ_LINE) # Retrieve just the root part of the name. i = fnroot (Memc[file], Memc[root], SZ_LINE) if (template == NULL) { call salloc (template, TEMPLATE_SIZE+1, TY_CHAR) call strcpy (Memc[root], Memc[template], TEMPLATE_SIZE) } # Loop through all the groups. call amovkd (0.0d0, Memd[obs_sum], n_diodes) call amovkd (0.0d0, Memd[obs_sumsq], n_diodes) call amovki (0, Memi[obs_n_data], n_diodes) obs_max = 0.0d0 dnum_obs_sum = 0.0d0 obs_event_diodes = 0 obs_event_rate = 0 obs_arate = 0.0d0 obs_brate = 0.0d0 obs_rrate = 0.0d0 do group = 1, n_groups { # Open the image and retrieve some of the extracted # data. call dst_open_data (Memc[file], group, grp_image, extracted_image, special_image) # Get end-of-exposure time. iferr (call imgstr (grp_image, "PKTTIME", Memc[time], SZ_LINE)) { call sprintf (Memc[tmp_string], SZ_LINE, "Could not get time for observation %s\n") call pargstr (Memc[file]) call eprintf (Memc[tmp_string]) Memc[time] = EOS } # Get data, determine grp statistics. call sprintf (Memc[tmp_string], SZ_LINE, "%s[%d]") call pargstr (Memc[root]) call pargi (group) grp_data = imgl1d (grp_image) call dst_remove_particles (Memd[grp_data], Memd[particle_limit], n_diodes, Memc[tmp_string], Memc[time], event_table, event_row, Memd[sum], Memd[sumsq], Memi[n_adds], n_particles, n_event_diodes, grp_event_rate, maximum, mean, sigma) # Get the special diodes. call dst_special (extracted_image, special_image, grp_arate, grp_brate, grp_rrate) # Determine the isolated diode value. if (IS_INDEFI(dnum)) dnum_value = INDEFD else dnum_value = Memd[grp_data+dnum-1] # Add sums to observation sums. obs_max = max (obs_max, maximum) call aaddd (Memd[sum], Memd[obs_sum], Memd[obs_sum], n_diodes) call aaddd (Memd[sumsq], Memd[obs_sumsq], Memd[obs_sumsq], n_diodes) call aaddi (Memi[n_adds], Memi[obs_n_data], Memi[obs_n_data], n_diodes) if (!IS_INDEFI(dnum)) dnum_obs_sum = dnum_obs_sum + Memd[grp_data+dnum-1] obs_event_diodes = obs_event_diodes + n_event_diodes obs_event_rate = obs_event_rate + grp_event_rate obs_arate = obs_arate + grp_arate obs_brate = obs_brate + grp_brate obs_rrate = obs_rrate + grp_rrate tot_particles = tot_particles + n_particles # Write out the information. grp_brate = grp_brate / N_BACK_DIODES grp_rrate = grp_rrate / N_RAD_DIODES call dst_stat_write (grp_table, Memc[tmp_string], mean, sigma, maximum, n_event_diodes, grp_event_rate, grp_arate, grp_brate, grp_rrate, dnum_value, grp_row) } # Now that a group has been dealt with, find statistics. call dst_stat (Memd[obs_sum], Memd[obs_sumsq], Memi[obs_n_data], n_diodes, mean, sigma) if (IS_INDEFI(dnum)) dnum_value = INDEFD else dnum_value = dnum_obs_sum / n_groups # Check if this is no longer part of the same observation, If it # isn't, then calculate the total statistics and reset. if (strncmp (Memc[template], Memc[root], TEMPLATE_SIZE) != 0) { if (IS_INDEFI(dnum)) dnum_value = INDEFD else { dnum_value = dnum_tot_sum / n_grps call printf ("Dark rate for diode %d = %6.6g\n") call pargi (dnum) call pargd (dnum_value) } call dst_stat (Memd[tot_sum], Memd[tot_sumsq], Memi[tot_n_data], n_diodes, mean, sigma) tot_arate = tot_arate tot_brate = tot_brate / N_BACK_DIODES / n_grps tot_rrate = tot_rrate / N_RAD_DIODES / n_grps call dst_stat_write (tot_table, Memc[template], mean, sigma, tot_max, tot_event_diodes, tot_event_rate, tot_arate, tot_brate, tot_rrate, dnum_value, tot_row) call dst_find_noisy (Memc[template], Memd[tot_sum], Memi[tot_n_data], n_diodes, dthresh, dnsig, verbose, noisy_table, noisy_row) call dst_dark_write (dark_table, Memc[template], Memd[tot_sum], Memi[tot_n_data], n_diodes) call strcpy (Memc[root], Memc[template], TEMPLATE_SIZE) call amovd (Memd[obs_sum], Memd[tot_sum], n_diodes) call amovd (Memd[obs_sumsq], Memd[tot_sumsq], n_diodes) call amovi (Memi[obs_n_data], Memi[tot_n_data], n_diodes) tot_max = obs_max n_grps = n_groups dnum_tot_sum = dnum_obs_sum tot_event_diodes = obs_event_diodes tot_event_rate = obs_event_rate tot_arate = obs_arate tot_brate = obs_brate tot_rrate = obs_rrate } # Else, add into the total. else { call aaddd (Memd[tot_sum], Memd[obs_sum], Memd[tot_sum], n_diodes) call aaddd (Memd[tot_sumsq], Memd[obs_sumsq], Memd[tot_sumsq], n_diodes) call aaddi (Memi[tot_n_data], Memi[obs_n_data], Memi[tot_n_data], n_diodes) tot_max = max (obs_max, tot_max) n_grps = n_grps + n_groups dnum_tot_sum = dnum_tot_sum + dnum_obs_sum tot_event_rate = tot_event_rate + obs_event_rate tot_event_diodes = tot_event_diodes + obs_event_diodes tot_arate = tot_arate + obs_arate tot_brate = tot_brate + obs_brate tot_rrate = tot_rrate + obs_rrate } # Write the observations statistics. obs_arate = obs_arate obs_brate = obs_brate / N_BACK_DIODES / n_groups obs_rrate = obs_rrate / N_RAD_DIODES / n_groups call dst_stat_write (obs_table, Memc[root], mean, sigma, obs_max, obs_event_diodes, obs_event_rate, obs_arate, obs_brate, obs_rrate, dnum_value, obs_row) } # Close all potential images. call imunmap (grp_image) if (extracted_image != NULL) call imunmap (extracted_image) if (special_image != NULL) call imunmap (special_image) # Make sure everything is written to the total statistics table. # If an isolated diode is to be examined, write out is current rate. if (IS_INDEFI(dnum)) dnum_value = INDEFD else { dnum_value = dnum_tot_sum / n_grps call printf ("Dark rate for diode %d = %6.6g\n") call pargi (dnum) call pargd (dnum_value) } call dst_stat (Memd[tot_sum], Memd[tot_sumsq], Memi[tot_n_data], n_diodes, mean, sigma) tot_arate = tot_arate tot_brate = tot_brate / N_BACK_DIODES / n_grps tot_rrate = tot_rrate / N_RAD_DIODES / n_grps call dst_stat_write (tot_table, Memc[template], mean, sigma, tot_max, tot_event_diodes, tot_event_rate, tot_arate, tot_brate, tot_rrate, dnum_value, tot_row) call dst_find_noisy (Memc[template], Memd[tot_sum], Memi[tot_n_data], n_diodes, dthresh, dnsig, verbose, noisy_table, noisy_row) call dst_dark_write (dark_table, Memc[template], Memd[tot_sum], Memi[tot_n_data], n_diodes) # Write out number of particle events. if (verbose) if (tot_particles > 0) { call printf ("%d particle events found.\n") call pargi (tot_particles) } else call printf ("No particle events found.\n") # That's all folks. call tbtclo (dark_table) call tbtclo (event_table) call tbtclo (noisy_table) call tbtclo (tot_table) call tbtclo (obs_table) call tbtclo (grp_table) call imtclose (list) call sfree (sp) end #--------------------------------------------------------------------------- # End of dst_darkstats #---------------------------------------------------------------------------