Initial commit

1 files changed, 142 insertions, 0 deletions

diff --git a/src/libcf/cf_apply_dead_time.c b/src/libcf/cf_apply_dead_time.c
new file mode 100644
index 0000000..a714a28
--- /dev/null
+++ b/src/libcf/cf_apply_dead_time.c
@@ -0,0 +1,142 @@
+/*****************************************************************************
+ *              Johns Hopkins University
+ *              Center For Astrophysical Sciences
+ *              FUSE
+ *****************************************************************************
+ *
+ * Synopsis:    cf_apply_dead_time(fitsfile *header, long nevents, float *time,
+ *                      float *weight, long nseconds, float *timeline, 
+ *			float *ids_dtc, float *elec_dtc)
+ *
+ * Description: Applies dead-time correction to all photons.  For HIST data,
+ *		uses mean of dead-time arrays.  Writes mean dead-time
+ *		values to file header.  Issue warning if DET_DEAD > 1.5.
+ *
+ * Arguments:   fitsfile  *header       Pointer to FITS file containing the
+ *                                      header of the intermediate data file
+ *              long      nevents       The number of events
+ *              float     *time         An array of event times
+ *              float     *weight       An array of event weights
+ *              long      nseconds      The number of timeline values
+ *              float     *timeline     An array of timeline times
+ *              float     *aic_rate     Active Image Counter rate
+ *              float     *ids_dtc      Output of cf_ids_dead_time
+ *              float     *elec_dtc     Output of cf_elec_dead_time
+ *
+ * Calls:
+ *
+ * Return:      0  on success
+ *
+ * History:     08/01/03   1.1    wvd   Begin work
+ *		08/04/03   1.2    wvd	Convert aic_rate to type short.
+ *		11/25/03   1.3    wvd	For HIST data, require that
+ *					mean_tot_dtc be >= 1.
+ *		11/26/03   1.4    wvd	Change aic_rate to type float.
+ *		02/09/04   1.5    wvd	Re-write with new logic.
+ *              04/09/04   1.6    bjg   Include stdlib.h
+ *		11/07/05   1.7    wvd	Add a test to catch bad DTC values.
+ *		03/10/06   1.8    wvd	Clean up i/o.
+ *		12/29/06   1.9    wvd	Check DET_DEAD rather than TOT_DEAD.
+ *		04/07/07   1.10   wvd	Clean up compiler errors.
+ *
+ ****************************************************************************/
+#include <stdlib.h>
+#include <string.h>
+#include "calfuse.h"
+
+int
+cf_apply_dead_time(fitsfile *header, long nevents, float *time, float *weight,
+	long nseconds, float *timeline, float *ids_dtc, float *elec_dtc)
+{
+    char CF_PRGM_ID[] = "cf_apply_dead_time";
+    char CF_VER_NUM[] = "1.10";
+
+    char comment[FLEN_CARD], datestr[FLEN_CARD], instmode[FLEN_VALUE];
+    char hkexists[FLEN_VALUE];
+    double ctimeb, ctimee, eng_time;
+    int  errflg=0, status=0, timeref;
+    int  det_flag=FALSE, ids_flag=FALSE;
+    long j, k;
+    float *tot_dtc, det_dead, ids_dead, mean_tot_dtc=0.;
+
+    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 header keywords.  Set flag if DET_DEAD > 1.5 */
+    FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
+    FITS_read_key(header, TSTRING, "HKEXISTS", hkexists, NULL, &status);
+    FITS_read_key(header, TFLOAT, "DET_DEAD", &det_dead, NULL, &status);
+    if (det_dead > 1.5) det_flag = TRUE;
+
+    /*
+     *  If there is no housekeeping file, AND the engineering snapshot times
+     *  are bad, AND the mean value of any DTC array is > 1.5, set the 
+     *  combined DTC array to unity and issue a warning.
+     */
+    if(!strncasecmp(hkexists, "N", 1)) {
+	FITS_read_key(header, TDOUBLE, "CTIME_B", &ctimeb, NULL, &status);
+        FITS_read_key(header, TDOUBLE, "CTIME_E", &ctimee, NULL, &status);
+        if ((eng_time = (ctimee-ctimeb) * 86400.) < 1.) {
+            FITS_read_key(header, TFLOAT, "IDS_DEAD", &ids_dead, NULL, &status);
+	    if (ids_dead > 1.5) ids_flag = TRUE;
+	}
+    }
+
+    /*
+     *  Compute total dead-time correction and its mean.
+     */
+    tot_dtc = (float *) cf_malloc(sizeof(float) * nseconds);
+    if (det_flag || ids_flag) {
+	for (k = 0; k < nseconds; k++) tot_dtc[k] = 1.0;
+	mean_tot_dtc = 1.0;
+    }
+    else {
+	for (k = 0; k < nseconds; k++) {
+	    tot_dtc[k] = elec_dtc[k] * ids_dtc[k];
+            mean_tot_dtc += tot_dtc[k];
+	}
+        mean_tot_dtc /= nseconds;
+    }
+
+    /* In TTAG mode, scale weights by tot_dtc as a function of time. */
+    if (!strncmp(instmode, "TTAG", 4)) for (j=k=0; j<nevents; j++) {
+        while(timeline[k+1]-FRAME_TOLERANCE < time[j] && k+1 < nseconds)
+            k++;
+        weight[j] *= tot_dtc[k];
+    }
+
+    /* In HIST mode, scale all weights by mean dead-time correction. */
+    else if (mean_tot_dtc > 1.)
+	for (j = 0; j < nevents; j++) 
+            weight[j] *= mean_tot_dtc;
+    else mean_tot_dtc = 1.;
+
+    cf_verbose(2, "Mean total dead-time correction: %f", mean_tot_dtc);
+    FITS_update_key(header, TFLOAT, "TOT_DEAD", &mean_tot_dtc, NULL, &status);
+
+    if (det_flag || ids_flag) {
+        FITS_write_comment(header, "  ", &status);
+	if (ids_flag) {
+	    sprintf(comment,
+		"IDS_DEAD > 1.5.  Setting dead-time correction to unity.");
+            cf_if_warning(comment);
+            FITS_write_comment(header, comment, &status);
+	}
+	if (det_flag) {
+	    sprintf(comment,
+		"DET_DEAD > 1.5.  Setting dead-time correction to unity.");
+            cf_if_warning(comment);
+            FITS_write_comment(header, comment, &status);
+	}
+        fits_get_system_time(datestr, &timeref, &status);
+        sprintf(comment, "CalFUSE v%s   %.10s", CALFUSE_VERSION, datestr);
+        FITS_write_comment(header, comment, &status);
+        FITS_write_comment(header, "  ", &status);
+    }
+
+    free (tot_dtc);
+    cf_proc_update(header, CF_PRGM_ID, "COMPLETE");
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+    return 0;
+}