Initial commit

1 files changed, 396 insertions, 0 deletions

diff --git a/src/libcf/cf_find_spectra.c b/src/libcf/cf_find_spectra.c
new file mode 100644
index 0000000..412a1e4
--- /dev/null
+++ b/src/libcf/cf_find_spectra.c
@@ -0,0 +1,396 @@
+/*****************************************************************************
+ *              Johns Hopkins University
+ *              Center For Astrophysical Sciences
+ *              FUSE 
+ *****************************************************************************
+ *
+ * Synopsis:    cf_find_spectra(fitsfile *header, long nevents,
+ *		float *weight, float, *xfarf, float *yfarf, char *chan,
+ *		unsigned char *timeflags, unsigned char *locflags,
+ *		int airglow_centroid)
+ *
+ * Description: Determines the spectral centroid for each aperture.
+ *		First, we compute the airglow centroid for the LWRS
+ *		aperture and determine its offset relative to the
+ *		tabulated value.  For the MDRS and HIRS apertures, we
+ *		apply this offset to the tabulated centroid to get the
+ *		airglow centroid.  For each aperture, we search for a
+ *		target spectrum in a narrow range around the airglow
+ *		centroid.  If we find it, we use it; otherwise, we use
+ *		the airglow centroid.  If airglow_centroid = TRUE, use
+ *		airglow centroid for the target spectrum.
+ *
+ *
+ * Arguments:   fitsfile  *header       Pointer to the location of the FITS
+ *                                      header of the Intermediate Data File
+ *              long      nevents       Number of photons in the file
+ *              float     *weight       scale factor for each photon event
+ *              float     *xfarf        X position of each photon event
+ *              float     *yfarf        Y position of each photon event
+ *          unsigned char *chan         Pointer to the array containing
+ *                                      channel information for each event
+ *          unsigned char *timeflags    Pointer to the array containing
+ *					time flags for each event
+ *          unsigned char *locflags     Pointer to the array containing
+ *					location flags for each event
+ *		int  airglow_centroid	If TRUE, use airglow lines to compute
+ *					centroid of target spectrum.
+ *
+ * Calls:       
+ *
+ * Return:      0 on success
+ *
+ * History:     04/18/03    1.1  wvd    Initial work
+ *              05/20/03    1.2  rdr    Added call to cf_proc_check 
+ *              07/29/03    1.3  rdr    Deal with case where no source is
+ *                                      present (i.e. Aperture = RFPT)
+ *              08/21/03    1.6  wvd    Change channel to unsigned char.
+ *              09/08/03    1.7  wvd    Fix typo in verbose statement.
+ *              09/18/03    1.8  wvd    Stop iteration if
+ *					cf_calculate_y_centroid returns -1.
+ *              10/22/03    1.9  wvd    Rewrite program.
+ *              11/12/03    1.10 wvd    In HIST mode, use default centroid
+ *					for non-target channels.
+ *              04/09/04    1.11 bjg    Include stdlib.h and stdio.h
+ *                                      Change formats in printf to 
+ *                                      match arg types.
+ *              04/19/04    1.12 wvd    Sum background only over detector
+ *					active area.  Scale intrinsic count
+ *					rate by EXPTIME and detector width.
+ *              05/03/04    1.13 wvd    Move call to cf_identify_channel
+ *					from this routine to cf_remove_motions.
+ *              06/02/04    1.14 wvd    Don't need PHA array; can use locflags
+ *					instead.  Use time flags to exclude
+ *					bursts, etc.
+ *              04/21/05    1.15 wvd    Rewrite program using airglow lines
+ *					(if available) to constrain search
+ *					for target centroid.  Don't bother
+ *					reading model background files.
+ *              06/03/05    1.16 wvd    Fix bug in tabulation of exected Y
+ *					centroids. Remove unused variables.
+ *              02/01/06    1.17 wvd    Tighten region over which centroids
+ *					are computed.
+ *					If target centroid is too far from
+ *					airglow centroid, reject it.
+ *					If ycent is too far from default value,
+ *					use default, not airglow value.
+ *					Change ysum and norm to doubles.
+ *              08/30/06    1.18 wvd    Added airglow_centroid argument.
+ *		03/14/08    1.19 wvd	In HIST mode, there are no spectra in
+ *					other apertures to cause confusion,
+ *					so we need not require that the 
+ *					target centroid be within 30 pixels
+ *					of the airglow value.
+ *		11/25/08    1.20 wvd	Don't compare HIST centroids to
+ *					either airglow or tabulated centroid.
+ *
+ ****************************************************************************/
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include "calfuse.h"
+
+#define YLENGTH		NYMAX/16
+
+int
+cf_find_spectra(fitsfile *header, long nevents,
+	float *weight, float *xfarf, float *yfarf, unsigned char *channel,
+	unsigned char *timeflags, unsigned char *locflags, int airglow_centroid)
+{
+    char aper[FLEN_VALUE], instmode[FLEN_VALUE], quality[FLEN_VALUE];
+    char key[FLEN_KEYWORD], filename[FLEN_VALUE], detector[FLEN_VALUE];
+    double norm, ysum;
+    float bkgd, peak, ycent, yshift=0;
+    float yairg[NYMAX], ydist[NYMAX], ysigma[NYMAX];
+    float ybin[YLENGTH], sbin[YLENGTH], ycent_tab[8];
+    float dy, sigma, xmin, xmax, ycent_air;
+    int  active_ap[2], airglow=FALSE, target_ap, status=0;
+    int  bkgd_num, bkgd_min[8], bkgd_max[8];
+    int  chan_num, errflg = 0, lwrs;
+    int  chan_order[] = { 3, 2, 1, 7, 6, 5 };
+    int  center, high, low, npeak, hdu;
+    int  spex_sic, spex_lif, emax, emax_sic, emax_lif, extended;
+    long i, j, k;
+    fitsfile  *chidfits, *parmfits;
+
+    char CF_PRGM_ID[] = "cf_find_spectra";
+    char CF_VER_NUM[] = "1.20";
+
+    /* Initialize error checking */
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+    /* Enter a time stamp into the log */
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+    /* Check whether routine is appropriate for this data file. */
+    if ( (errflg = cf_proc_check(header, CF_PRGM_ID)) ) return errflg;
+
+    /* Initialize Y-distribution arrays. */
+    for (i = 0; i < NYMAX; i++)
+	yairg[i] = ydist[i] = ysigma[i] = 0.;
+    for (i = 0; i < YLENGTH; i++)
+	ybin[i] = sbin[i] = 0.;
+
+    /* Read header keywords from IDF. */
+    FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
+    FITS_read_key(header, TSTRING, "DETECTOR", detector, NULL, &status);
+    FITS_read_key(header, TSTRING, "APERTURE", aper, NULL, &status);
+    extended = cf_source_aper(header, active_ap);
+
+    FITS_read_key(header, TINT, "BKGD_NUM", &bkgd_num, NULL, &status);
+    for (i = 0; i < bkgd_num; i++) {
+	sprintf(key, "BKG_MIN%ld", i);
+        FITS_read_key(header, TINT, key, bkgd_min+i, NULL, &status);
+	sprintf(key, "BKG_MAX%ld", i);
+        FITS_read_key(header, TINT, key, bkgd_max+i, NULL, &status);
+    }
+
+    /* Read expected channel centroids from CHID_CAL file. */
+    FITS_read_key(header, TSTRING, "CHID_CAL", filename, NULL, &status);
+    FITS_open_file(&chidfits, cf_cal_file(filename), READONLY, &status);
+    for (i = 1; i < 8; i++) {
+	if (i==4) continue;
+        if ((i == active_ap[0] || i == active_ap[1]) && !extended) hdu=i+1;
+	else hdu=i+9;
+	FITS_movabs_hdu(chidfits, hdu, NULL, &status);
+	FITS_read_key(chidfits, TFLOAT, "CENTROID", ycent_tab+i, NULL, &status);
+    }
+    FITS_read_key(chidfits, TFLOAT, "XMIN", &xmin, NULL, &status);
+    FITS_read_key(chidfits, TFLOAT, "XMAX", &xmax, NULL, &status);
+    FITS_close_file(chidfits, &status);
+
+    /* 
+     * Generate photon and airglow arrays as a function of YFARF.
+     * Ignore events near the detector edge.
+     */
+    xmin += 250;
+    xmax -= 250;
+    for (i = 0; i < nevents; i++) {
+	if (!(timeflags[i] & ~TEMPORAL_DAY) &&
+	    (xfarf[i] > xmin) && (xfarf[i] < xmax)) {
+	    if (!locflags[i])
+	        ydist[cf_nint(yfarf[i])] += weight[i];
+	    else if (locflags[i] == LOCATION_AIR)
+	        yairg[cf_nint(yfarf[i])] += weight[i];
+	}
+    }
+    /* Set error array equal to photon counts (= variance). */
+    for (i = 0; i < NYMAX; i++)
+	ysigma[i] = ydist[i];
+
+    /* If exposure was taken in TTAG mode, estimate BKGD level */
+    bkgd = 0.;
+    if (!strcmp(instmode, "TTAG")) {
+	k = 0;
+	for (i = 0; i < bkgd_num; i++) {
+	    for (j = bkgd_min[i]; j <= bkgd_max[i]; j++) {
+		bkgd += ydist[j];
+		k++;
+	    }
+	}
+	bkgd /= k;
+	cf_verbose(3, "Mean background level: %.1f", bkgd);
+    }
+
+    /* Subtract mean background and bin data by 16 pixels. */
+    for (i = 0; i < NYMAX; i++) {
+	j = i / 16;
+	ybin[j] += ydist[i] - bkgd;
+	sbin[j] += ysigma[i];
+    }
+    for (i = 0; i < YLENGTH; i++) sbin[i] = sqrt(sbin[i]);
+
+    /* Read extraction parameters from PARM_CAL file. */
+    FITS_read_key(header, TSTRING, "PARM_CAL", filename, NULL, &status);
+    FITS_open_file(&parmfits, cf_parm_file(filename), READONLY, &status);
+    FITS_read_key(parmfits,TINT,"SPEX_SIC",&spex_sic,NULL,&status);
+    FITS_read_key(parmfits,TINT,"SPEX_LIF",&spex_lif,NULL,&status);
+    FITS_read_key(parmfits,TINT,"EMAX_SIC",&emax_sic,NULL,&status);
+    FITS_read_key(parmfits,TINT,"EMAX_LIF",&emax_lif,NULL,&status);
+    FITS_close_file(parmfits,&status);
+
+    /* Determine centroid for each channel. */
+    for (k = 0; k < 6; k++) {
+	chan_num = chan_order[k];
+
+	/* Is this the LWRS aperture? */
+	if (chan_num == 3 || chan_num == 7)
+	    lwrs = TRUE;
+	else 
+	    lwrs = FALSE;
+
+	/* Is it the target aperture? */
+	if (chan_num == active_ap[0] || chan_num == active_ap[1])
+	    target_ap = TRUE;
+	else
+	    target_ap = FALSE;
+
+        /* If we're in histogram mode, and this is not the target aperture,
+	   move to next aperture. */
+        if (!strncmp(instmode, "HIST", 4) && !target_ap) {
+	    if (lwrs) yshift = 0;
+	    continue;
+	}
+
+	/* First, we compute a centroid from the airglow lines.
+	   For the LWRS aperture, search within 70 pixels of
+	   the expected centroid.  Compute the offset between the
+	   measured and tabulated centroids.  For the MDRS and HIRS
+	   apertures, apply the offset computed from the LWRS aperture
+	   to the tabulated centroid.
+	   Disregard regions near the bottom and top of the detector. */
+
+	if (lwrs) {
+	    center = cf_nint(ycent_tab[chan_num]);
+	    dy = 70;
+            if ((low  = center-dy) < 16) low = 16;
+            if ((high = center+dy) > NYMAX-50) high = NYMAX-50;
+	    norm = ysum = 0.;
+	    for (i = low; i <= high; i++) {
+		ysum += i * yairg[i];
+		norm += yairg[i];
+	    }
+
+	    /* If the airglow lines contain at least 33 events, use their
+	       centroid.  If not, use the tabulated centroid. */
+	    if (norm > 33.) {
+		ycent = ysum / norm;
+		yshift = ycent - ycent_tab[chan_num];
+		airglow = TRUE;
+	    }
+	    else {
+		ycent = ycent_tab[chan_num];
+		yshift = 0;
+		airglow = FALSE;
+	    }
+	}
+ 		/* For MDRS and HIRS apertures */
+	else
+	    ycent = yshift + ycent_tab[chan_num];
+
+	/* Remember airglow centroid */
+	ycent_air = ycent;
+
+	/* Now we find the centroid for the target spectrum.  If we are
+	   in a target aperture and user has specified YCENT, use it. */
+
+	if (target_ap && ((chan_num<5 && spex_lif>=0 && spex_lif<1024) ||
+	    (chan_num>4 && spex_sic>=0 && spex_sic<1024))) {
+            if (chan_num<5) ycent=spex_lif;
+            else ycent=spex_sic;
+            sprintf(key, "YCENT%1d", chan_num);
+            FITS_update_key(header, TFLOAT, key, &ycent, NULL, &status);
+	    sprintf(quality, "HIGH");
+            sprintf(key, "YQUAL%1d", chan_num);
+            FITS_update_key(header, TSTRING, key, quality, NULL, &status);
+            cf_verbose(1, "Channel %d: User-specified Y centroid = %.2f",
+                chan_num, ycent);
+	    continue;
+        }
+
+	/* If we must find the target spectrum ourselves, use the binned
+	   data array to improve the S/N.  Search within 2 binned pixels
+	   of ycent, however it was computed. */
+
+	center = cf_nint(ycent / 16.);
+	dy = 2;
+	if ((low = center - dy) < 3) low = 3;
+	if ((high = center + dy) > YLENGTH-4) high = YLENGTH-4;
+
+	npeak = 0;
+	peak = -1.0E5;
+	for (i = low; i <= high; i++) {
+	    if (ybin[i] > peak) {
+		peak = ybin[i];
+		npeak = i;
+	    }
+	}
+
+	/* If peak is significant, compute Y centroid for events falling
+	   within 40 pixels of peak.  If this number is within 30 pixels
+	   of airglow centroid, use it.  If not, use airglow centroid.
+	   Note that we require a higher significance for the SiC LWRS
+	   channel on side 1, because it sits on an elevated background. */
+
+	if (!strncmp(detector, "1", 1) && chan_num == 7) sigma = 9.;
+	else sigma = 5.;
+
+	sprintf(quality, "UNKNOWN");
+
+	if (peak > sigma * sbin[npeak] && !airglow_centroid) {
+	    center = npeak * 16 + 8;
+	    dy = 40;
+	    if ((low = center - dy) < 32) low = 32;
+	    if ((high = center + dy) > NYMAX-50) high = NYMAX-50;
+	    norm = ysum = 0.;
+	    for (i = low; i <= high; i++) {
+		ysum += i * ydist[i];
+		norm += ydist[i];
+	    }
+	    ycent = ysum / norm;
+
+	    /* Don't test offset for HIST data. (wvd, 03/14/2008) */
+	    if (fabs(ycent - ycent_air) < 30 || !strncmp(instmode, "HIST", 4)) {
+	    	sprintf(quality, "HIGH");
+	    	cf_verbose(3, "Channel %d: default centroid: %.2f",
+			chan_num, ycent_tab[chan_num]);
+	    	if (airglow) cf_verbose(3, "Channel %d: airglow centroid: %.2f",
+			chan_num, ycent_air);
+	    	cf_verbose(3,"Channel %d: counts in peak = %.0f, limit = %.0f",
+			chan_num, peak, sigma * sbin[npeak]);
+	    	cf_verbose(2, "Channel %d: using target centroid = %.2lf",
+			chan_num, ycent);
+	    }
+	}
+	/* If peak is not significant or ycent differs too much from
+	   the airglow centroid, use the airglow centroid. */
+	if (!strncmp(quality, "UNKNOWN", 1)) {
+	    if (airglow) {
+		sprintf(quality, "MEDIUM");
+	        cf_verbose(3, "Channel %d: default centroid: %.2f",
+		    chan_num, ycent_tab[chan_num]);
+		cf_verbose(3,"Channel %d: peak = %.0f, limit = %.0f",
+		    chan_num, peak, sigma * sbin[npeak]);
+		cf_verbose(2, "Channel %d: using airglow centroid = %.2lf",
+		    chan_num, ycent);
+	    }
+	    else {
+	        sprintf(quality, "LOW");
+		cf_verbose(3,"Channel %d: peak = %f, limit = %f",
+		    chan_num, peak, sigma * sbin[npeak]);
+		cf_verbose(2, "Channel %d: using default centroid = %.2lf",
+		    chan_num, ycent);
+	    }
+	}
+
+	/* If centroid lies too far from the default value, use default. */
+        if (chan_num<5) emax=emax_lif; else emax=emax_sic;
+	    /* Don't test offset for HIST data. (wvd, 11/25/2008) */
+        if (fabs(ycent - ycent_tab[chan_num]) > emax && !strncmp(instmode, "TTAG", 4)) {
+	    if (target_ap) {
+        	cf_verbose(1, "Computed centroid (%.1lf) is too far from "
+		    "expected value.", ycent);
+	        cf_verbose(1, "Channel %d: using default centroid of %.1f",
+		    chan_num, ycent_tab[chan_num]);
+	    }
+	    else {
+	        cf_verbose(2, "Channel %d: using default centroid of %.1f",
+		    chan_num, ycent_tab[chan_num]);
+	    }
+            ycent = ycent_tab[chan_num];
+	    sprintf(quality, "LOW");
+        }
+
+	/* Write results to file header */
+        sprintf(key, "YCENT%1d", chan_num);
+        FITS_update_key(header, TFLOAT, key, &ycent, NULL, &status);
+        sprintf(key, "YQUAL%1d", chan_num);
+        FITS_update_key(header, TSTRING, key, quality, NULL, &status);
+    }
+
+    cf_proc_update(header, CF_PRGM_ID, "COMPLETE");
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done Processing");
+
+    return status;
+}