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