/***************************************************************************** * 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 #include #include #include #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; }