Initial commit

1 files changed, 164 insertions, 0 deletions

diff --git a/src/libcf/cf_fifo_dead_time.c b/src/libcf/cf_fifo_dead_time.c
new file mode 100644
index 0000000..67ac065
--- /dev/null
+++ b/src/libcf/cf_fifo_dead_time.c
@@ -0,0 +1,164 @@
+/*****************************************************************************
+ *              Johns Hopkins University
+ *              Center For Astrophysical Sciences
+ *              FUSE
+ *****************************************************************************
+ *
+ * Synopsis:    cf_fifo_dead_time(fitsfile *header, long nevents,
+ *		float *ptime, long nseconds, float *ttime, float *aic_rate,
+ *		float *ids_dtc)
+ *
+ * Description: Searches high-count-rate observations for data drop-outs
+ *		due to FIFO overflows.  If found, calculates time-dependent
+ *		scale factor to correct for lost photon events.
+ *
+ * Arguments:   fitsfile  *header       Pointer to FITS file containing the
+ *                                      header of the intermediate data file
+ *              long      nevents       Number of photon events in the file
+ *              float     *ptime        Detection time for each photon
+ *              long      nseconds      The number of timeline values
+ *              float     *ttime        Tabulated times in the timeline
+ *              float     *aic_rate     Active Image Counter array
+ *              float     *ids_dtc      IDS dead-time correction array
+ *					(modified)
+ *
+ * Calls:
+ *
+ * Return:      0  on success
+ *
+ * History of cf_screen_fifo_overflow:
+ *
+ *		02/10/04   1.1    wvd   Begin work
+ *              02/17/05   1.2    wvd   Place parentheses around assignment
+ *                                      used as truth value.
+ *              11/09/05   1.3    wvd   In the current scheme, the rate_lif
+ *					and rate_sic arrays contain values
+ *					from the housekeeping file and do
+ *					not reflect data loss due to FIFO
+ *					overflows.  Now we calculate the
+ *					count-rate array directly from the
+ *					photon time array.
+ *              02/23/06   1.4    wvd   Change length of data_rate array
+ *					from nevents to nseconds.
+ *
+ * History:	12/29/06   1.1    wvd	Derived from cf_screen_fifo_overflow.
+ *              01/15/07   1.2    wvd   When count rate is zero, correct one
+ *					second before to two seconds after.
+ *              02/08/07   1.3    wvd   If TTPERIOD = 0, set it to 1.
+ *              04/07/07   1.4    wvd   Clean up compiler warnings.
+ *
+ ****************************************************************************/
+
+#include <stdlib.h>
+#include "calfuse.h"
+
+int
+compute_count_rate(long nevents, float *ptime, long nseconds,
+    float *ttime, long **data_rate)
+{
+    float delta_max = 1.0 + FRAME_TOLERANCE;
+    long j,k;
+
+    *data_rate = (long *) cf_calloc(nseconds, sizeof(long));
+    for (j=k=0; j<nevents; j++) {
+        while(ttime[k+1]-FRAME_TOLERANCE < ptime[j] && k+1 < nseconds) k++;
+        if ((ptime[j] - ttime[k]) < delta_max) (*data_rate)[k]++;
+    }
+    return 0;
+}
+
+int
+cf_fifo_dead_time(fitsfile *header, long nevents, float *ptime,
+    long nseconds, float *ttime, float *aic_rate, float *ids_dtc)
+{
+    char CF_PRGM_ID[] = "cf_fifo_dead_time";
+    char CF_VER_NUM[] = "1.4";
+
+    char elecfile[FLEN_VALUE];
+    int  errflg=0, status=0;
+    int  fifo_drain_rate;
+    int  *flag_array, found_drop_out = FALSE;
+    long *data_rate, k;
+    float mean_ids_rate, ttperiod;
+    float mean_ids_dtc;
+    fitsfile *elecfits;
+
+    /* Exit if data were obtained in HIST mode. */
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+    if ((errflg = cf_proc_check(header, CF_PRGM_ID))) return errflg;
+
+    /* Read IDS drain rate from calibration file. */
+    FITS_read_key(header, TSTRING, "ELEC_CAL", elecfile, NULL, &status);
+    FITS_open_file(&elecfits, cf_cal_file(elecfile), READONLY, &status);
+    FITS_read_key(elecfits, TINT, "FIFO_DRN", &fifo_drain_rate, NULL, &status);
+    FITS_close_file(elecfits, &status);
+    cf_verbose(4, "fifo_drain_rate = %d", fifo_drain_rate);
+
+    /* Read interval between time stamps. */
+    FITS_read_key(header, TFLOAT, "TTPERIOD", &ttperiod, NULL, &status);
+    if (ttperiod < 1E-6) ttperiod = 1.0;
+
+    /* Compute mean IDS count rate. */
+    mean_ids_rate = 0.;
+    for (k=0; k<nseconds; k++) 
+        mean_ids_rate += aic_rate[k];
+    mean_ids_rate /= nseconds;
+    mean_ids_rate += 1./ttperiod;
+    cf_verbose(4, "mean_ids_rate = %f", mean_ids_rate);
+
+    /* If mean IDS count rate > 3200, look for data drop-outs. */
+    if (mean_ids_rate > 3200 && nseconds > 4) {
+	compute_count_rate(nevents, ptime, nseconds, ttime, &data_rate);
+	flag_array = (int *) cf_calloc(nseconds, sizeof(int));
+
+        /* If count rate goes to zero, set flag_array.  Where possible,
+	   set flags from one second before to two seconds after. */
+	for (k=0; k<1; k++) if (data_rate[k] == 0) 
+	    flag_array[k] = flag_array[k+1] = flag_array[k+2] = TRUE;
+
+	for ( ; k<nseconds-2; k++) if (data_rate[k] == 0) 
+	    flag_array[k-1] = flag_array[k] = flag_array[k+1] = 
+		flag_array[k+2] = TRUE;
+
+	for ( ; k<nseconds-1; k++) if (data_rate[k] == 0) 
+	    flag_array[k-1] = flag_array[k] = flag_array[k+1] = TRUE;
+
+	for ( ; k<nseconds; k++) if (data_rate[k] == 0) 
+	    flag_array[k-1] = flag_array[k] = TRUE;
+
+	for (k=0; k<nseconds; k++) if (data_rate[k] == 0) {
+	   found_drop_out = TRUE;
+	   break;
+	}
+
+	/* If there are any drop-outs, modify IDS dead-time correction. */
+	if (found_drop_out)
+	     for ( ; k < nseconds; k++) if (flag_array[k])
+		ids_dtc[k] = (aic_rate[k] - 1./ttperiod) / fifo_drain_rate;
+
+	free(data_rate);
+	free(flag_array);
+    }
+
+    /*
+     *  Calculate mean dead-time correction.  Write to file header.
+     */
+    if (found_drop_out) {
+	mean_ids_dtc = 0;
+	for (k = 0; k < nseconds; k++) mean_ids_dtc += ids_dtc[k];
+	mean_ids_dtc /= nseconds;
+	cf_verbose(1, "IDS FIFO overflow detected.");
+        cf_verbose(2, "Modified mean IDS dead-time correction: %f",
+	    mean_ids_dtc);
+        FITS_update_key(header, TFLOAT, "IDS_DEAD", &mean_ids_dtc, NULL,
+	    &status);
+    }
+
+    /*
+     *  Clean up and go home.
+     */
+    cf_proc_update(header, CF_PRGM_ID, "COMPLETE");
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+    return 0;
+}