From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/libcf/cf_scale_bkgd.c | 956 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 956 insertions(+) create mode 100644 src/libcf/cf_scale_bkgd.c (limited to 'src/libcf/cf_scale_bkgd.c') diff --git a/src/libcf/cf_scale_bkgd.c b/src/libcf/cf_scale_bkgd.c new file mode 100644 index 0000000..556c0c7 --- /dev/null +++ b/src/libcf/cf_scale_bkgd.c @@ -0,0 +1,956 @@ +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Synopsis: cf_scale_bkgd(infits, nevents, x_in, y_in, w_in, channel, + * timeflags, locflags, ngood, good_index, dtime, ntime, + * &binx, &biny, &nx_lif, &ny_lif, + * &ymin_lif, img_lif, &nx_sic, &ny_sic, &ymin_sic, &img_sic) + * + * Arguments: *fitsfile infits Input IDF file + * long nevents Number of photon events recorded + * float x_in, y_in Coordinates of each photon + * float w_in Weighting factor for each photon + * unsigned char channel Aperture ID for each photon + * unsigned char timeflags Photon status flags regarding time + * unsigned char locflags Photon status flags regarding location + * long ngood Number of good photons + * long good_index List of indices for the good photons + * long dtime Length of the daytime exposure + * long ntime Length of the night time exposure + * int binx, biny Binning factors of the background + * array along the X and Y axes + * int nx_lif, nx_sic X Size of background array + * int ny_lif, ny_sic Y Size of background array + * int ymin_lif, ymin_sic Y value corresponding to the + * first row of the bkgd image + * float *img_lif, *img_sic 2D background image (binned) + * + * Description: + * (a) Determine the total counts in a background sample region for data + * taken during the night and during the day. Use the data in the + * IDF event list, after screening out the geocoronal lines. + * (b) Read estimated intrinsic count rates from the background characterization + * file and estimate the intrinsic counts in the background region from + * the exposure times and number of pixels in the region. + * (c) Subtract the estimated intrinsic counts from the measured total to + * get an estimate of the scattered light contribution for both day + * and night portions of the observation + * (d) Use the scattered light counts to scale daytime and night background + * images obtained from the background calibration files + * (e) Combine the two background images to obtain the final background model + * + * A previous assumption that the intrinsic background is constant was + * found to be incorrect. The program now iterates the assumed intrinsic + * background count rate (assumed constant in day and night) and checks the + * calculated model against the observations by comparing the intensity + * variations along the x direction. In this comparison, we mask out regions + * that are affected by geocoronal contamination. + * + * Caveats: The "intrinsic" component of the background is neither intrinsic + * nor constant in day and night. It is better called the "spatially uniform" + * component and should be computed independently for day and night portions + * of the exposure. wvd 01/24/2006 + * + * History: 01/11/03 v1.1 rdr Started work - adapted program + * from cf_make_ttag_bkgd + * 03/01/03 v1.3 wvd Correct use of pointer in FITS_read_key + * 03/11/03 v1.4 wvd Update documentation regarding use of + * unsigned char. + * 03/28/03 v1.5 rdr corrected error in freeing storage + * 03/31/03 v1.6 rdr corrected error in assigning subarrays + * for lif and sic apertures + * 05/08/03 v1.7 rdr read limits to the background region + * from the header of the IDF. Also + * changed print statements to cf_verbose. + * 05/20/03 v1.8 rdr - Make a model using only exposure time + * if no background sample regions have + * been specified (generally for HIST data) + * - Look into parm file to determine +* whether to generate a background model + * 05/22/03 v1.9 rdr Correct problems occuring when exptime=0 + * 06/04/03 v1.10 rdr Change method for determining limits of + * output arrays + * 09/08/03 v1.12 wvd Use cf_read_col in get_limits() + * 09/29/03 v1.13 wvd Fix bug in population of 2-D bkgd + * array; clean up i/o. + * 10/07/03 v1.14 wvd Reject airglow photons using a + * bit-wise comparison of locflags. + * 10/31/03 v1.15 wvd Increase aperture pad from 10 to 15. + * Change channel to unsigned char. + * 10/31/03 v1.16 wvd Use cf_read_col in get_geocoronal_mask() + * 03/22/04 v1.17 wvd Change get_limits so that sizes of + * bkgd and probability arrays are equal. + * 04/08/04 v1.18 bjg Replace NXMAX by nx + * Remove unused variables + * Change format in printf to match + * arg type. + * 04/27/04 v1.19 wvd Replace RUN_MKBK with RUN_BKGD + * throughout. + * 07/01/04 v1.20 wvd Fix bug: nxb was undefined for + * RUN_BKGD = NO. + * 03/08/05 v1.21 wvd Convert to cf_scale_bkgd. + * Treat data and models the same way. + * Modify get_xvar_img. Add + * make_bad_pixel_mask and sum_good_pixels. + * 06/15/05 v1.22 wvd BUG FIX: get_limits always read the + * probability arrays for point-source + * targets from the WGTS_CAL file. Now + * it distinguishes between point-source + * and extended targets. + * 01/20/06 v1.23 wvd BUG FIX: Apply uniform X limits to + * data and background models. + * BUG FIX: Allow output bkgd images + * to fall off the detector if target + * centroid requires it. + * 01/24/06 v1.24 wvd BUG FIX: Test y coordinate of photon + * before adding it to the data image. + * 02/01/06 v1.25 wvd BUG FIX: Make npix a long in + * sum_good_pixels. Correct error in + * extraction of sub-arrays. + * + ****************************************************************************/ + +#include +#include +#include +#include +#include "calfuse.h" + + +/********************************************************************* + * + * GET_XVAR_IMG + * + * Subroutine to compute the intensity variation along the X direction + * of a background image. Pixels flagged as bad are excluded. + * bkgimg = image from which to extract the data + * maskimg = mask for data image (1 = good pixel) + * nxb = size of bkgimg in x + * nxsam = number of elements in the xvars array + * xvars = array of x variations + * + ***********************************************************************/ + + +int +get_xvar_img(float *bkgimg, float *maskimg, int nxb, int nxsam, float *xvars) +{ + int nbin; + long i, j; + + nbin = nxb/nxsam ; + for (i = 0; i < NYMAX; i++) + for (j = 0; j < nxb; j++) + xvars[j/nbin] += bkgimg[i*nxb+j] * maskimg[i*nxb+j]; + + return 0 ; +} + +/********************************************************************* + * + * MAKE_MODEL + * + * Subroutine to combine a day and night calibration image to generate the + * final background model: + * npix = number of pixels in the full background image + * nsam = number of pixels in the background regions of the detector + * nxb = size of the x dimension of the background images + * intcr = assumed intrinsic count rate + * expnight, expday = exposure times during the night and day + * data_sum_night, data_sum_day = sum of weights within the background + * sample regions for night and day exposures + * model_sum_night, model_sum_day = number of counts in the unscaled + * night and day background calibration images + * - background regions only + * bkgimgn = image containing the night calibration image. + * bkgimgd = image containing the daytime calibration image + * bkgimg = final background image + * maskimg = mask defining good background regions + * + *****************************************************************************/ + +int +make_model(long npix, long nsam, int nxb, float intcr, float expnight, + float expday, double data_sum_night, double data_sum_day, + double model_sum_night, double model_sum_day, + float *bkgimgn, float *bkgimgd, float *bkgimg, float *maskimg) +{ + int bin_size ; + long i; + float int_counts; + float sfd, sfn, intnight, intday, scatnight, scatday ; + + /* bin_size = binning factor of the background image */ + bin_size = NXMAX / nxb ; + + cf_verbose(3, "\n\tConstructing the background model\n") ; + cf_verbose(3, "total counts in the data: night= %.1f, day = %.1f", + data_sum_night, data_sum_day); + cf_verbose(3, "intrinsic count rate = %.2f x 10^-7 ",intcr ) ; + cf_verbose(3, "exptime night = %.1f, exptime day = %.1f ", + expnight, expday ); + + /* Sum the intrinsic background counts */ + intnight = nsam * bin_size * intcr * expnight * 1.e-7 ; + intday = nsam * bin_size * intcr * expday * 1.e-7 ; + cf_verbose(3, "total intrinsic counts: night= %.1f , day= %.1f ", + intnight, intday) ; + + /* Sum of scattered light in data (measured - intrinsic) */ + scatnight = data_sum_night - intnight ; + if (scatnight <= 0.) scatnight = 0. ; + scatday = data_sum_day - intday ; + if (scatday <= 0.) scatday = 0. ; + cf_verbose(3, "total scattered light count: night = %.1f, day = %.1f", + scatnight, scatday); + cf_verbose(3, "total counts in unscaled background model: " + "night = %.2f, day = %.2f", model_sum_night, model_sum_day); + + /* Calculate the scale factor for the night scattered light image */ + sfn = 0. ; + if (model_sum_night > 0) sfn = scatnight / model_sum_night ; + + /* Calculate the scale factor for the day scattered light image */ + sfd = 0. ; + if (model_sum_day > 0) sfd = scatday / model_sum_day ; + cf_verbose(3, "scale factors for background images: " + "night = %.1f, day = %.1f ", sfn, sfd) ; + + /* Calculate the intrinsic counts per binned pixel */ + int_counts = bin_size * intcr * (expnight + expday) * 1.e-7 ; + cf_verbose(3, "intrinsic counts per binned pixel = %.6f ", int_counts ) ; + + /* Use the scattered light scaling factors and the tabulated intrinsic + count rate to model the final background image */ + + for (i = 0; i < npix; i++) + bkgimg[i] = (sfn * bkgimgn[i]) + (sfd * bkgimgd[i]) + int_counts; + + return 0 ; +} + + +/********************************************************************* + * + * MAKE_MODEL_EXPTIME : + * + * Subroutine to combine a day and night calibration image based on + * exposure times only to generate the final background model: + * npix = total number of pixels in the image + * nxb = length of the x axis in the background images + * intcr = assumed intrinsic count rate + * expnight, expday = exposure times during the night and day + * bkging = image containing the night calibration image. + * bkgimgd = image containing the daytime calibration image + * bkgimg = final background model + * + *********************************************************************/ + + +int +make_model_exptime(long npix, int nxb, float intcr, float expnight, + float expday, float *bkgimgn, float *bkgimgd, float *bkgimg) +{ + + int i, bin_size ; + float intnight, intday ; + + /* Compute the binning factor in x */ + bin_size = NXMAX/nxb ; + + /* Calculate the expected intrinsic counts in the background */ + intnight = bin_size * intcr * expnight * 1.e-7 ; + intday = bin_size * intcr * expday * 1.e-7 ; + + /* Combine the day and night images based on exposure times */ + for (i=0 ; i 0) { + maskimg[ndx] = 0. ; + n_air++ ; + } + } + } + + free(gxmin); + free(gxmax); + free(gymin); + free(gymax); + + cf_verbose(3, "Number of bkgd pixels affected by geocoronal lines = %ld", + n_air) ; + + /* Count up and return the number of good pixels. */ + n_good = 0; + for (i = 0; i < npix; i++) if(maskimg[i] > 0) n_good++; + cf_verbose(3, "Number of bkgd pixels flagged as good = %ld", n_good) ; + + *nsam = n_good; + return 0 ; +} + + +/****************************************************************************** + * + * GET_LIMITS + * + * Set limits of background subimage to match those of probability array. + * + * + ******************************************************************************/ + +int get_limits(fitsfile *infits, int extended, int aperture, + int *ymin, int *ymax) +{ + + char wgtsfile[FLEN_VALUE], ycentname[FLEN_KEYWORD]; + int hdu, nywt, status=0; + float ycentv, wgt_cent; + fitsfile *wgtsfits; + + cf_verbose(3, "Determining Y limits for aperture %d background array.", + aperture); + + /* Read the measured centroid of the target spectrum. */ + sprintf(ycentname,"YCENT%1d", aperture); + FITS_movabs_hdu(infits, 1, NULL, &status); + FITS_read_key(infits, TFLOAT, ycentname, &ycentv, NULL, &status); + + /* Read the centroid and size of the probability array. */ + FITS_read_key(infits, TSTRING, "WGTS_CAL", wgtsfile, NULL, &status); + FITS_open_file(&wgtsfits, cf_cal_file(wgtsfile), READONLY, &status); + hdu = extended * 8 + aperture + 1; + cf_verbose(3, "WGTS_CAL: extended = %d, aperture = %d, hdu = %d", + extended, aperture, hdu); + FITS_movabs_hdu(wgtsfits, hdu, NULL, &status); + FITS_read_key(wgtsfits, TINT, "NAXIS2", &nywt, NULL, &status); + FITS_read_key(wgtsfits, TFLOAT, "WGT_CENT", &wgt_cent, NULL, &status); + FITS_close_file(wgtsfits, &status); + + cf_verbose(3, "Spectrum centroid=%6.2f, weights centroid=%5.2f, " + "weights ny=%4d", ycentv, wgt_cent, nywt) ; + + /* Note that we allow these values to exceed the limits of the + * detector. This is required if the bkgd arrays are to cover + * the same area as the 2-D probability and data arrays. */ + + *ymin = cf_nint (ycentv - wgt_cent); + *ymax = *ymin + nywt - 1; + + return status ; +} + + +/********************************************************************* + * + * SUM_GOOD_PIXELS + * + * Subroutine multiplies image and mask arrays, then sums the result. + * image = image from which to extract the data + * mask = mask for data image (1 = good pixels) + * npix = size of image and mask arrays + * + ***********************************************************************/ + + +int +sum_good_pixels(float *image, float *mask, long npix, double *image_sum) +{ + + long i; + + *image_sum = 0; + + for (i = 0; i < npix; i++) + *image_sum += image[i] * mask[i]; + + return 0 ; +} + + +/****************************************************************************** + * + * MAIN PROGRAM BEGINS HERE + * + ******************************************************************************/ + + +int cf_scale_bkgd(fitsfile *infits, long nevents, float *x_in, float *y_in, + float *w_in, unsigned char *channel, unsigned char *timeflags, + unsigned char *locflags, long ngood, long *good_index, long dtime, + long ntime, int *binx, int *biny, int *nx_lif, int *ny_lif, + int *ymin_lif, float **img_lif, int *nx_sic, int *ny_sic, int *ymin_sic, + float **img_sic) +{ + + char CF_PRGM_ID[] = "cf_scale_bkgd"; + char CF_VER_NUM[] = "1.25"; + + fitsfile *bkgdfits, *scrnfits, *bchrfits, *parmfits ; + char buffer[FLEN_CARD], bkgdfile[FLEN_CARD], scrnfile[FLEN_CARD]; + char bchrfile[FLEN_CARD], keyname[FLEN_CARD], parmfile[FLEN_CARD] ; + char run_bkgd[FLEN_VALUE]; + char comment[FLEN_CARD], datestr[FLEN_CARD]; + int nx, nxb, ny, status, hdutype, frow, felem, niter; + int src_ap[2], nullval, anynull, bkgd_num; + int ymin, ymax, phalow, nbinx, nxsam, bkgdtype ; + int bin_size, extended, keyval, timeref; + double data_sum_day, data_sum_night, data_sum; + double model_sum_night, model_sum_day, model_sum ; + float *xvarg, *bkgarr, *xvars0, *xvarm0 ; + float *bkgimg, *bkgimgd, *bkgimgn, *xvars, *xvarm; + float floatnull, exptime, expnight, expday, xvars_ave ; + float intcrarray[9], intcr, scattot, rms, rms0, dcr ; + float intcnt, scat_int_ratio, zero, mf, numtot, xvars_tot, xvarm_tot ; + float *dataimg, *maskimg; + short *ylim ; + long fpixel, npix, i, j, k, nsam; + long nbgood, npts, ndx, ndx1, ndx2, bkgd_ndx; + + /* Initialize error checking */ + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Started execution."); + + /* Initialize CFITSIO variables */ + anynull = 0 ; + hdutype = 0 ; + frow = 1 ; + felem = 1 ; + floatnull = 0. ; + fpixel = 1 ; + status = 0 ; + + /* Initialize other variables */ + data_sum_day=0.; + data_sum_night=0.; + zero = 0.; + + /* Open parameter file and check whether user wants a background model. */ + FITS_movabs_hdu(infits, 1, &hdutype, &status) ; + FITS_read_key(infits, TSTRING, "PARM_CAL", parmfile, NULL, &status) ; + cf_verbose(3, "parameter file is %s ", parmfile) ; + FITS_open_file(&parmfits,cf_parm_file(parmfile), READONLY, &status) ; + FITS_read_key(parmfits, TSTRING, "RUN_BKGD", run_bkgd, NULL, &status) ; + FITS_close_file(parmfits, &status) ; + cf_verbose(3, "RUN_BKGD = %s", run_bkgd) ; + + /* If the RUN_BKGD = NO, then fill background array with zeros. */ + if (!strncmp(run_bkgd, "N", 1) || !strncmp(run_bkgd, "n", 1)) { + cf_verbose(1, "RUN_BKGD = NO. No background subtraction.") ; + FITS_write_comment(infits, " ", &status); + FITS_write_comment(infits, "RUN_BKGD = NO. No background subtraction.", + &status); + fits_get_system_time(datestr, &timeref, &status); + sprintf(comment, "CalFUSE v%s %.10s", CALFUSE_VERSION, datestr); + FITS_write_comment(infits, comment, &status); + FITS_write_comment(infits, " ", &status); + + bin_size = 1 ; + nx=nxb=NXMAX ; + ny=NYMAX ; + npix=nx * ny ; + bkgimg = (float *) cf_calloc(npix, sizeof(float)) ; + goto fin ; + } + +/*********************************************************************** + * + * Set up the parameters for the analysis + * + ***********************************************************************/ + + /* Set day and night exposure times */ + expday = (float) dtime ; + expnight = (float) ntime ; + exptime = expday + expnight ; + + /* Get name of background characterization file */ + FITS_read_key(infits, TSTRING, "BCHR_CAL",bchrfile, NULL, &status) ; + cf_verbose(3, "BCHR parameter file is %s ",bchrfile) ; + + /* Get name of background calibration file */ + FITS_read_key(infits, TSTRING, "BKGD_CAL", bkgdfile, buffer, &status); + cf_verbose(3,"Background calibration file = %s ",bkgdfile ) ; + + /* Get name of screening file */ + FITS_read_key(infits, TSTRING, "SCRN_CAL",scrnfile, NULL, &status) ; + cf_verbose(4, "screening file is %s ",scrnfile) ; + + /* Read pha screening limit from the input file header */ + FITS_read_key(infits,TINT,"PHALOW",&phalow, NULL, &status) ; + if (phalow > 8) phalow=8 ; + cf_verbose(2, "pha screening limit = %d", phalow) ; + + /* Determine the type of background model to make by checking the + BKGDTYPE flag in the screening file. If BKGDTYPE > 0, do a full + model. If BKGDTYPE < 0, calculate a model based on exposure time. */ + FITS_open_file(&scrnfits,cf_parm_file(scrnfile), READONLY, &status) ; + FITS_read_key(scrnfits, TINT, "BKGDTYPE",&bkgdtype, NULL, &status) ; + FITS_close_file(scrnfits, &status) ; + + /* Read the intrinsic count rate from the BCHR file */ + FITS_open_file(&bchrfits,cf_parm_file(bchrfile), READONLY, &status) ; + FITS_movabs_hdu(bchrfits, 3, &hdutype, &status) ; + FITS_read_img(bchrfits, TFLOAT, fpixel, 9, &nullval, intcrarray, + &anynull, &status) ; + FITS_close_file(bchrfits, &status) ; + intcr=intcrarray[phalow] ; + cf_verbose(2,"Initial guess for the intrinsic count rate = " + "%.1f x 10^-7 c/s/pixel ", intcr) ; + + /* Read the limits of the background sample regions from the header + of the input file */ + FITS_read_key(infits, TINT, "BKGD_NUM", &bkgd_num, NULL, &status) ; + cf_verbose(3, "Number of background sample regions = %d ", bkgd_num) ; + + ylim = (short *) cf_calloc(bkgd_num*2, sizeof(short) ) ; + for (i=0; i 1) for (i=0; i= NYMAX) continue; + bkgd_ndx = j*nxb + cf_nint(x_in[ndx])/bin_size; + if (maskimg[bkgd_ndx] > 0) { + dataimg[bkgd_ndx] += w_in[ndx]; + if (timeflags[ndx] & TEMPORAL_DAY) data_sum_day += w_in[ndx] ; + else data_sum_night += w_in[ndx]; + nbgood++ ; + } + } + + cf_verbose(3, "Number of photons events in the background region = %ld", + nbgood); + cf_verbose(3, "Sum of daytime weights = %.1lf ", data_sum_day) ; + cf_verbose(3, "Sum of nighttime weights = %.1lf ", data_sum_night) ; + + /* If the background scattered light contribution is too large, the + background may be contaminated, perhaps by a source in one of the + other apertures. This is common in crowded fields or with nebulocity. + In this case, we calculate the background model based on the + exposure times, as done for the histogram data */ + + /* Check the level of the scattered-light background */ + intcnt = intcr * nsam * bin_size * exptime * 1.e-7 ; + scattot = data_sum_day + data_sum_night - intcnt ; + scat_int_ratio = 0. ; + if (intcnt > 0) scat_int_ratio = scattot / intcnt ; + cf_verbose(3, "int cnt = %.1f, scat cnt = %.1f , scat/int = %f ", + intcnt, scattot, scat_int_ratio) ; + + if (scat_int_ratio > 10) { + cf_verbose(1, "The scattered light background is too large."); + bkgdtype = -1; + } + else { + cf_verbose(2, + "Scattered light appears reasonable. Proceeding with analysis.") ; + } + + } /* Analysis of the data ends here. */ + + /* If the background model is to be calculated using only the exposure + time, then do it now. */ + if (bkgdtype < 0) { + cf_verbose(1, "Constructing background based only on exposure time.") ; + make_model_exptime( npix, nxb, intcr, expnight, expday, + bkgimgn, bkgimgd, bkgimg) ; + goto fin ; + } + +/**************************************************************************** + * + * In this section, we attempt to determine the intrinsic count rate. + * The scattered light shows little structure in the Y dimension, so + * we project all of the arrays onto the X axis. Only good pixels + * in the user-specified background regions are included in the sum. + * + ****************************************************************************/ + + /* Make a scattered light reference image assuming an intrinsic count + rate of zero. */ + make_model_exptime(npix, nxb, zero, expnight, expday, + bkgimgn, bkgimgd, bkgimg) ; + + /* We bin the data by another factor of 8 in X */ + /* nbinx = number of detector pixels in one bin */ + /* nxsam = number of samples in the array */ + nbinx=128 ; + nxsam = NXMAX/nbinx ; + cf_verbose(3, "X variation array information: " + "nbinx = %d, nxsam= %d ", nbinx, nxsam) ; + + xvars = (float *) cf_calloc(nxsam,sizeof(float)) ; + xvarm = (float *) cf_calloc(nxsam,sizeof(float)) ; + xvars0 = (float *) cf_calloc(nxsam,sizeof(float)) ; + xvarm0 = (float *) cf_calloc(nxsam,sizeof(float)) ; + xvarg = (float *) cf_calloc(nxsam,sizeof(float)) ; + + /* Project each image onto the X axis */ + get_xvar_img(bkgimg, maskimg, nxb, nxsam, xvarm0) ; + get_xvar_img(dataimg, maskimg, nxb, nxsam, xvars0) ; + get_xvar_img(maskimg, maskimg, nxb, nxsam, xvarg) ; + + /* Calculate the average counts per "good" bin in the data array. */ + xvars_ave = 0; + for (i = j = 0; i < nxsam; i++) + if (xvarg[i] > 0) { + xvars_ave += xvars0[i]; + j++; + } + xvars_ave /= j; + cf_verbose (3, "xvars_ave = %.1f (should be > 100)", xvars_ave); + + /* If the average number of counts in the xvars0 array (which holds + the observed intensity variation along the x direction) is small, + then we can't determine an intrinsic count rate from the data and + we use the tabulated rate. If the count is large enough, then we + can do a more detailed fit of the model to the observations. */ + + if (xvars_ave > 100) { /* Begin fit to intrinsic background */ + + cf_verbose(1, "Calculating a full background model.") ; + cf_verbose(3, "\n\tInterpolating to derive intrinsic count rate\n") ; + + /* Sum the counts in the data and model arrays. */ + sum_good_pixels(dataimg, maskimg, npix, &data_sum); + sum_good_pixels(bkgimg, maskimg, npix, &model_sum); + + cf_verbose(3, "Total counts in initial scattered-light model = %7.0lf", + model_sum); + cf_verbose(3, "Intrinsic counts summed over good regions = %7.0f", + (intcnt = intcr * nsam * bin_size * exptime * 1.e-7)); + cf_verbose(3, "Sum of above = %7.0lf", + model_sum + intcnt); + cf_verbose(3, "Total counts in good regions of data array = %7.0lf", + data_sum); + + /* Sum the counts in the data and model x variation arrays */ + xvarm_tot = 0.; + xvars_tot = 0.; + for (i = j = 0; i < nxsam; i++) if (xvarg[i] > 0) { + xvarm_tot += xvarm0[i]; + xvars_tot += xvars0[i]; + } + cf_verbose(3, "Total counts in scattered-light xvar array = %7.0f", + xvarm_tot); + cf_verbose(3, "Total counts in data x variation array = %7.0f", + xvars_tot); + + /* Convert the x variation arrays from total counts to units of + 1.e-7 c/s/pix, the units of the intrinsic count rate. */ + mf = bin_size * exptime * 1.e-7 ; + for (i=0 ; i 0) { + xvars0[i] /= xvarg[i] * mf ; + xvarm0[i] /= xvarg[i] * mf ; + } + + /* Sum the counts in the rescaled xvar arrays */ + xvarm_tot = 0.; + xvars_tot = 0.; + for (i = j = 0; i < nxsam; i++) if (xvarg[i] > 0) { + xvars_tot += xvars0[i]; + xvarm_tot += xvarm0[i]; + j++; + } + cf_verbose(3, "Mean xvar count rates in units of 10^-7 c/s/pixel:"); + cf_verbose(3, " Initial model = %.1f", xvarm_tot/j); + cf_verbose(3, " Intrinsic counts = %.1f", intcr); + cf_verbose(3, " Data array = %.1f", xvars_tot/j); + + /* Determine the value of the intrinsic background by comparing the X + variations of the data (after removing a guess for the intrinsic + background) with the model scattered light variations, making sure + that the final average count rate is the same for both. Iterate + until the intrinsic background changes by less than 0.1 units. */ + + dcr = 1. ; + niter = 0 ; + rms0 = 1.e10 ; + while ( fabs(dcr) > 0.1) { + + /* Subtract intrinsic counts from the data array and sum it. */ + xvars_tot = 0. ; + for (i=0; i 0) { + xvars[i] = xvars0[i] - intcr ; + xvars_tot += xvars[i] ; + } + + /* Scale the model array to match the data. Sum it. */ + numtot = 0. ; + for (i=0 ; i 0) { + xvarm[i] = xvarm0[i] * (xvars_tot / xvarm_tot) ; + numtot += xvarm[i] ; + } + + /* Sum the absolute difference between the data and model arrays. */ + rms = 0. ; + for (i=0 ; i 0) + rms += (fabs(xvars[i]-xvarm[i])) ; + cf_verbose(3, "intr cr= %5.2f, data-model = %6.2f, " + "dat cnts = %5.0f, mod cnts=%5.0f", intcr,rms,xvars_tot,numtot) ; + + niter++ ; + if (niter > 20) break ; + + /* if the current rsm value is greater than the preceeding one, + reverse direction and reduce the step size by 50% */ + if ( rms > rms0 ) dcr = -0.5*dcr ; + intcr += dcr ; + + /* initialize the next iteration */ + rms0 = rms ; + } + } /* End of fit to intrinsic background */ + + else { + cf_verbose(1, "Background count too small for a detailed fit."); + cf_verbose(1, "Assuming an intrinsic count rate of %.1f cts/s", intcr ) ; + } + + /* Sum the good pixels in the unscaled day and night background images. */ + sum_good_pixels(bkgimgn, maskimg, npix, &model_sum_night); + sum_good_pixels(bkgimgd, maskimg, npix, &model_sum_day); + + /* Make the final background model */ + make_model(npix, nsam, nxb, intcr, expnight, expday, + data_sum_night, data_sum_day, model_sum_night, model_sum_day, + bkgimgn, bkgimgd, bkgimg, maskimg) ; + + /* Sum the counts in the final background image. */ + sum_good_pixels(bkgimg, maskimg, npix, &model_sum); + + cf_verbose(3, "Total counts in final background model = %.0lf", model_sum); + cf_verbose(3, " measured counts = %.0lf", + data_sum_day+data_sum_night); + + free(bkgimgn) ; + free(bkgimgd) ; + free(dataimg) ; + free(maskimg) ; + free(xvars) ; + free(xvars0) ; + free(xvarm) ; + free(xvarm0) ; + free(xvarg) ; + + fin: + + /************************************************************************** + * + * Extract the subarrays covering the active LiF and SiC apertures + * + ***************************************************************************/ + + extended = cf_source_aper(infits, src_ap) ; + cf_verbose(3, "\n\tSelecting subarrays\n") ; + cf_verbose(3, "source apertures: lif=%d, sic=%d ",src_ap[0], src_ap[1]) ; + + /* Populate some output values */ + *binx=bin_size ; + *biny=1 ; + + for (k=0; k<2 ; k++) { /* Loop over LiF and SiC channels */ + + /* Read the y limits for the output array */ + get_limits(infits, extended, src_ap[k], &ymin, &ymax) ; + + /* Set up array to contain the image */ + npts = nxb * (ymax-ymin+1) ; + bkgarr = (float *) cf_calloc(npts, sizeof(float)) ; + + /* Fill the array */ + ndx1 = -1 ; + model_sum = 0 ; + for (i=ymin; i<=ymax ; i++) { + if (i < 0 || i >= ny) + ndx1 += nxb ; + else { + for (j=0; j= npix) break ; + ndx1++; + if (ndx1 >= npts) break ; + bkgarr[ndx1]=bkgimg[ndx2] ; + model_sum += bkgarr[ndx1] ; + } + } + } + + /* Populate the output arrays and values */ + if (k == 0) { + *nx_lif = nxb ; + *ny_lif = (ymax - ymin + 1) ; + *ymin_lif = ymin ; + *img_lif = bkgarr ; + } + else { + *nx_sic = nxb ; + *ny_sic = (ymax - ymin + 1) ; + *ymin_sic = ymin ; + *img_sic = bkgarr ; + } + + cf_verbose(2, "for aperture %d, nx=%d, ny=%d, ymin=%d, model_sum=%.0f", + src_ap[k], nxb, (ymax-ymin+1), ymin, model_sum) ; + + } /* End loop over LiF and SiC channels */ + + /* release memory for the background image */ + free(bkgimg) ; + + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution."); + return 0; +} -- cgit