/***************************************************************************** * Johns Hopkins University * Center For Astrophysical Sciences * FUSE ***************************************************************************** * * Synopsis: cf_extract_spectra options intermediate_file * * Description: This module computes a background model, reads the various * calibration files, and performs either a standard or optimal * extraction for the target spectrum in each of the LiF and * SiC channels. * * Note: some of the routines called by this module modify values * in the IDF file header. We want those changes to appear in * the final, extracted spectral files, but not in the archived * IDF, so we make a copy in memory of the IDF header and pass it * to the various routines. We delete it at the end. * * Arguments: input_file Intermediate data file * * Calibration files: * * Returns: 0 on successful completion * * Calls: cf_apply_filters, cf_scale_bkgd, cf_make_wave_array, * cf_rebin_and_flux_calibrate, cf_rebin_probability_array, * cf_standard_or_optimal_extraction, cf_optimal_extraction * cf_write_extracted_spectrum * * History: 02/13/03 peb 1.1 Begin work * 03/02/03 wvd 1.3 Add cf_standard_or_optimal_extraction() * 03/10/03 peb 1.4 Added verbose_level, unistd.h for * getopt portability, and remove debug * print statements * 04/08/03 wvd 1.6 Copy header of IDF into memory * so that subsequent routines can * modify it without changing IDF. * Delete this virtual header when done. * Read arrays from timeline table and * pass to cf_apply_filters. * 05/06/03 rdr 1.7 Read from column ERGCM2 rather than * ERGCM2S * 05/07/03 wvd 1.8 Change ergcm2 to ergcm2 throughout * 05/28/03 rdr 1.9 Read pothole mask * 06/09/03 rdr 1.10 Incorporate the tscreen flag * 06/11/03 wvd 1.11 Pass datatype to cf_read_col * 08/25/03 wvd 1.12 Change coltype from string to int in * cf_read_col. * 09/29/03 wvd 1.13 Move standard_or_optimal so that it * runs just once. Don't check return * values from subroutines. * 10/08/03 wvd 1.14 Change counts_out array to type long. * 10/16/03 wvd 1.15 If EXPTIME = 0, generate a pair of * null-valued spectra and exit. * 10/31/03 wvd 1.16 Change channel to unsigned char. * 12/12/03 wvd 1.18 Clean up i/o. * 03/16/04 wvd 1.19 Don't pass wave_out to optimal * extraction subroutine. * Change pycent from int to float. * 03/24/04 wvd 1.20 Eliminate got_no_data(); fold into * main routine. If either EXPTIME = 0 * or NEVENTS = 0, write null spectrum. * 04/26/04 wvd 1.21 Replace cf_rebin_and_flux_calibrate * background with cf_rebin_background. * cf_optimal_extraction now returns * flux_out and sigma_out in units of * counts; must flux-calibrate each. * 06/10/04 wvd 1.22 Don't pass time array from timeline * table to cf_apply_filters. * 06/11/04 peb 1.23 Add the -r option to override the * default rootname. * 01/28/05 wvd 1.24 If EXP_STAT < 0, generate a pair of * null-valued spectra and exit. * 03/09/05 wvd 1.25 Change cf_ttag_bkgd to cf_scale_bkgd, * pass weights array to it. * 04/12/05 wvd 1.27 Add -n flag, which forces extraction * of night-only spectrum. Its argument * is the name of a night-only BPM file. * 06/15/05 wvd 1.28 BUG FIX: program always read the * point-source probability arrays from * the WGTS_CAL file. Now uses value * of extended to determine which HDU * to read. * 11/23/05 wvd 1.29 Add flag to disable optimal * extraction from the command line. * 01/27/06 wvd 1.30 Add flag to force use of optimal * extraction from the command line. * 05/19/06 wvd 1.31 Don't discard spectra with non-zero * values of EXP_STAT. * 06/12/06 wvd 1.32 Change cf_if_warning to cf_verbose. * ****************************************************************************/ #include #include #include #include "calfuse.h" static char CF_PRGM_ID[] = "cf_extract_spectra"; static char CF_VER_NUM[] = "1.32"; int main(int argc, char *argv[]) { char *bpm_file=NULL, *outrootname = NULL; unsigned char *channel=NULL, *timeflags=NULL, *loc_flags=NULL; unsigned char *time_status=NULL; int extended, status=0, optc, j, tscreen=1; int binx, biny, bnx[2], bny[2], bymin[2], aper[2]; int night_only=FALSE, optimal=FALSE; int force_optimal=FALSE, override_optimal=FALSE; long nevents=0, ntimes=0, ngood=0, dtime=0, ntime=0, *good_index=NULL; float *weight=NULL, *x=NULL, *y=NULL, *lambda=NULL; float exptime, *bpmask=NULL ; float *bimage[2]={NULL, NULL}; fitsfile *memp, *header; char opts[] = "hn:ofr:sv:"; char usage[] = "Usage:\n" " cf_extract_spectra [-hsof] [-n bpm_filename] [-r rootname] [-v level]" " idf_file\n"; char option[] = "Options:\n" " -h: this help message\n" " -n: extract night-only spectrum using given bad-pixel map\n" " -o: disable optimal extraction\n" " -f: force optimal extraction\n" " -r: override the default rootname\n" " -s: do not perform screening on time flags\n" " -v: verbosity level (=1; 0 is silent)\n" ; verbose_level = 1; while ((optc = getopt(argc, argv, opts)) != -1) { switch(optc) { case 'h': printf("%s\n%s", usage, option); return 0; case 'n': bpm_file = optarg; night_only = TRUE; break; case 'o': override_optimal = TRUE; break; case 'f': force_optimal = TRUE; break; case 'r': outrootname = optarg; break; case 'v': verbose_level = atoi(optarg); break; case 's': tscreen = 0 ; break; } } cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); if (argc <= optind) { printf("%s", usage); return -1; } cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing"); /* Open input data file. */ FITS_open_file(&header, argv[optind], READONLY, &status); FITS_read_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status); FITS_movabs_hdu(header, 2, NULL, &status); nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight); nevents = cf_read_col(header, TFLOAT, "X", (void **) &x); nevents = cf_read_col(header, TFLOAT, "Y", (void **) &y); nevents = cf_read_col(header, TBYTE, "CHANNEL", (void **) &channel); nevents = cf_read_col(header, TBYTE, "TIMEFLGS", (void **) &timeflags); nevents = cf_read_col(header, TBYTE, "LOC_FLGS", (void **) &loc_flags); nevents = cf_read_col(header, TFLOAT, "LAMBDA", (void **) &lambda); FITS_movabs_hdu(header, 4, NULL, &status); ntimes = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &time_status); /* * Copy header of IDF into memory, close IDF, * and pass copy of header to subsequent routines. */ FITS_movabs_hdu(header, 1, NULL, &status); FITS_create_file(&memp, "mem://", &status); FITS_copy_hdu(header, memp, 0, &status); FITS_close_file(header, &status); header = memp; /* If night-only spectrum is requested, modify header accordingly. */ if (night_only) { cf_verbose(3, "Night-only spectrum requested from command line."); FITS_read_key(header, TFLOAT, "EXPNIGHT", &exptime, NULL, &status); FITS_update_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status); FITS_update_key(header, TSTRING, "DAYNIGHT", "NIGHT", NULL, &status); FITS_update_key(header, TSTRING, "BPM_CAL", bpm_file, NULL, &status); } extended = cf_source_aper(header, aper); /* Determine whether to attempt optimal extraction. */ if (override_optimal) { cf_verbose (1, "User set -o flag. No optimal extraction."); optimal=FALSE; } else if (force_optimal) { cf_verbose (1, "User set -f flag. Forcing use of optimal extraction."); optimal=TRUE; } else cf_standard_or_optimal_extraction(header, &optimal); /* Make list of all photons that satisfy requested screenings. */ cf_apply_filters(header, tscreen, nevents, timeflags, loc_flags, ntimes, time_status, &dtime, &ntime, &ngood, &good_index); /* Generate model backgrounds for LiF and SiC channels. */ cf_scale_bkgd(header, nevents, x, y, weight, channel, timeflags, loc_flags, ngood, good_index, dtime, ntime, &binx, &biny, &bnx[0], &bny[0], &bymin[0], &bimage[0], &bnx[1], &bny[1], &bymin[1], &bimage[1]); /* Run once for each of LiF and SiC target channels. */ for (j=0; j<2; j++) { int pny, valid_spectrum=TRUE; float pycent, wpc; long *counts_out=NULL, i, nout=0, nphotons=0; float *wave_out=NULL, *bkgd_out=NULL, *flux_out=NULL, *sigma_out=NULL; float *weights_out=NULL, *barray=NULL, *parray=NULL; short *bpix_out=NULL; unsigned char *channel_temp=NULL; /* Generate output wavelength array. */ cf_make_wave_array(header, aper[j], &nout, &wave_out); /* If no data, generate null-valued output arrays. */ for (i = 0; i < ngood; i++) { if (channel[good_index[i]] == aper[j]) { nphotons++; break; } } if (exptime < 1. || nphotons < 1) { bkgd_out = (float *) cf_calloc(nout, sizeof(float)); bpix_out = (short *) cf_calloc(nout, sizeof(short)); counts_out = (long *) cf_calloc(nout, sizeof(long)); flux_out = (float *) cf_calloc(nout, sizeof(float)); sigma_out = (float *) cf_calloc(nout, sizeof(float)); weights_out = (float *) cf_calloc(nout, sizeof(float)); valid_spectrum=FALSE; cf_verbose(1, "No data for aperture %d. Generating " "null-valued output spectrum.", aper[j]); } else { /* Extract target spectrum. */ /* Bin background model to match output wavelength array. */ cf_rebin_background(header, aper[j], nout, wave_out, binx, biny, bnx[j], bny[j], bimage[j], &barray); /* Use QUAL_CAL files to create a bad-pixel mask for this exposure. */ cf_make_mask(header, aper[j], nout, wave_out, bny[j], bymin[j], &bpmask); /* Bin the weights/probability array to match output wave array. */ cf_rebin_probability_array(header, extended, aper[j], nout, wave_out, &pny, &pycent, &parray); /* Use optimal extraction if possible, standard extraction otherwise. */ cf_optimal_extraction(header, optimal, aper[j], weight, y, channel, lambda, ngood, good_index, barray, bpmask, pny, pycent, parray, nout, wave_out, &flux_out, &sigma_out, &counts_out, &weights_out, &bkgd_out, &bpix_out); /* Optimal-extraction routine returns flux_out and sigma_out in units of counts. Must flux-calibrate both arrays. */ channel_temp = (unsigned char *) cf_malloc(sizeof(unsigned char) * nout); memset (channel_temp, aper[j], sizeof(unsigned char) * nout); FITS_read_key(header, TFLOAT, "WPC", &wpc, NULL, &status); FITS_update_key(header, TSTRING, "FLUX_COR", "PERFORM", NULL, &status) ; cf_convert_to_ergs(header, nout, flux_out, flux_out, channel_temp, wave_out); FITS_update_key(header, TSTRING, "FLUX_COR", "PERFORM", NULL, &status) ; cf_convert_to_ergs(header, nout, sigma_out, sigma_out, channel_temp, wave_out); free(channel_temp); /* Divide flux and errors by EXPTIME * WPC to get units of erg/cm2/s/A. */ for (i = 0; i < nout; i++) { flux_out[i] /= exptime * wpc; sigma_out[i] /= exptime * wpc; } } /* Write extracted spectrum to output file. */ cf_write_extracted_spectrum(header, aper[j], valid_spectrum, nout, wave_out, flux_out, sigma_out, counts_out, weights_out, bkgd_out, bpix_out, outrootname); free(counts_out); free(bkgd_out); free(weights_out); free(sigma_out); free(flux_out); free(wave_out); free(bpix_out) ; free(parray); free(barray); } cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); FITS_delete_file(header, &status); free(bimage[1]); free(bimage[0]); free(lambda); free(loc_flags); free(timeflags); free(channel); free(y); free(x); free(weight); free(time_status); cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing"); return 0; }