diff options
Diffstat (limited to 'src/libcf/cf_screen_burst.c')
| -rw-r--r-- | src/libcf/cf_screen_burst.c | 758 |
1 files changed, 758 insertions, 0 deletions
diff --git a/src/libcf/cf_screen_burst.c b/src/libcf/cf_screen_burst.c new file mode 100644 index 0000000..f80374e --- /dev/null +++ b/src/libcf/cf_screen_burst.c @@ -0,0 +1,758 @@ +/**************************************************************************** + * Johns Hopkins University + * Center for Astrophysical Sciences + * FUSE + **************************************************************************** + * + * Synopsis: cf_screen_burst(infits, nevents, ptime, x, y, locflags, + * gti, nsec, ttime, sflag, aic_rate, bkgd) + * + * Description: Determines the times for bursts and sets the burst + * flag in the timeline status array so that the busts + * can be removed later. + * + * Arguments: fitsfile infile Pointer to Intermediate Data File + * long nevents Number of photons + * float ptime Time of event + * float *x, *y X and Y positions of the photon + * byte locflags flags on photon location + * (LOC_FLGS column in IDF) + * GTI gti Structure with Good Time Intervals + * long nsec Length of timeline array + * float ttime Timeline array + * byte sflag Status flags in timeline table + * float aic_rate AIC count rate + * short bkgd Count rate for the background + * sample region + * + * Algorithm used: + * + * 1. Generate an array (flgarr) with 1 second intervals that + * contains flags indicating whether that time interval is good + * or not. All elements of this array are originally set to negative + * numbers, indicating bad times. + * 2. Read in the good time intervals from the initial screening + * analysis and set all good time elements in flgarr to positive + * values + * 3. Determine the median count per binned element in the entire time + * sequence (ignoring times with no data) and set all elements with + * a count > 5 times this median value to negative values (indicating + * a burst at that time). + * 4. Do a median filter of the modified time sequence (ignoring the neg + * numbers) and identify all times where the count per bin is greater + * than a specified value or is more than a specified number of + * standard deviations from the median. These points are flagged with + * negative numbers to indicate bursts. The screening parameters used + * are specified in the screening parameter file. + * 5. Iterate time-sequence analysis until no more points are removed. + * 6. Transfer the burst times from the cnt array to flgarr + * 7. Go through all detected photons. If the time of arrival for that + * photon has been flagged as bad in FLGARR then ignore it, + * otherwise copy the information to the output file. + * + * + * Calibration files: + * BCHR_CAL: location of the background sample regions + * SCRN_CAL: burst screening params + * + * Parameter file keywords used to control the burst screening process + * + * MNCNT : minimum enhancement of the background (in counts) needed to + * flag a burst (currently = 5) + * STDRE : minimum number of standard deviations from the local mean + * for a burst to be detected (currently=5) + * NBIN : binning factor (in seconds) for the tiome sequence before + * analysis. (currently=15 seconds) + * NSMED : Interval (in seconds) used in determining the median filter + * (currently=600) + * SRCFRAC : Fraction of the source intensity required for a + * burst to be removed. SCRFRAC=0.01 requires that a burst + * exceed 1% of the integrated source intensity before being + * removed (default 0.001) + * + * Returns: 0 on success + * + * History: 10/29/02 v1.1 RDR Adapted from cf_ttag_burst + * 11/14/02 v1.2 peb Tidied code and made it compatible + * with cf_screen_burst + * 12/09/02 v1.3 RDR - Adjusted procedure for filling + * burst flags in the timeline + * - Corrected error in the procedure + * that tabulates bursts properties. + * - Used LOCATION flags to exclude + * geocoronal emission and photons + * off the active detector. + * 12/20/02 v1.4 rdr changed status flag arrays to + * unsigned char + * 01/24/03 v1.5 rdr corrected problem with high + * count rates + * 03/01/03 v1.6 wvd Correct use of pointer in FITS_read_key + * 05/20/03 v1.7 rdr Add cf_proc_check + * 05/22/03 v1.8 wvd Exit if EXPTIME !> 0. + * Change cf_error_init to stderr. + * Implement cf_verbose throughout. + * 04/06/03 v1.9 rdr Exclude times which have previously + * had a screening bit set + * 09/03/03 v1.10 wvd Delete extraneous print statements + * 09/10/03 v1.11 wvd Change background array to type short + * Set burst flag only when flgarr[i] + * < -1, not 0, to prevent bursts + * from echoing LIMB/SAA/HV flags. + * 09/15/03 v1.12 wvd Use ttime[i] to populate bkgd array. + * Old scheme didn't work. + * 10/06/03 v1.13 wvd Change locflgs to locflags throughout. + * Use bit-wise logic to test sflag array. + * 06/04/04 v1.14 wvd Clean up i/o. + * 02/02/05 v1.15 wvd Read AIC count rate from timeline table + * rather than file header keywords. + * 02/02/05 v1.16 wvd Generate bkgd array for the entire + * exposure, not just good times. + * 08/30/05 v1.17 wvd Ignore photons with arrival times + * greater than MAX_EXPTIME. + * 04/09/06 v1.18 wvd In subroutine median_filter, require + * that nwin >= 1. + * 06/12/06 v1.19 wvd Call cf_verbose rather than + * cf_if_warning. + * + ***************************************************************************/ + +#include <math.h> +#include <string.h> +#include <stdio.h> +#include <stdlib.h> +#include <ctype.h> +#include "calfuse.h" + +/* + * Function to determine the median value of an array. + * Arguments are the array, + * the number of points in the array (npt) and the maximum count to + * use in the histogram analysis (nmax). The median value is returned. + * + * The analysis uses a histogram approach to determine the median + * value. Negative values of the array have been flagged as bursts + * and should be ignored in determining the median value. If all of + * the points in the array have been flagged, returns zero. + */ + +static int +median(int *array, int npts, int nmax) +{ + int i, ngood=0, *hist, medval; + + /* + * Set up an array to contain the histogram of the data. + */ + hist = (int *) cf_calloc((nmax+1), sizeof(int)); + /* + * Fill the histogram with the data. + * Negative points are not counted. + */ + for (i=0; i< npts; i++) { + int ndx = array[i]; + if (ndx > nmax-1) + ndx = nmax-1; + if (ndx >= 0 ) { + ngood ++; + hist[ndx]++; + } + } + + /* + * Determine the median value. Set the data to 0 if there + * are fewer than 2 points. + */ + if (ngood > 2 ) { + int hnum= (int) ngood/2, atot=0; + i=-1; + while ((atot <= hnum) && (i < nmax) ) + atot += hist[++i]; + medval = (int) i; + } + else + medval = 0; + + free(hist) ; + + return medval; +} + +/* + * Function to construct a median filter for an array: + * array = array for which the median filter is to be constructed + * npts = number of points in the array + * nflt = number of points over which to determine the median value + * medflt = array containing the median filter. + */ + +static void +median_filter(int *array, int npts, int nflt, int *medflt) +{ + int cntmax=0, i, hwidth, nwin, *win, medval, nsam ; + /* + * Determine the maximum value of the array. + */ + for (i = 0; i < npts; i++) + if (array[i] > cntmax) + cntmax=array[i]; + + cf_verbose(3,"Constructing median filter - max cnt = %d",cntmax) ; + + /* + * Determine the properties of the window over which to calculate + * the median values. Make sure that the window is not larger + * than the width of the data. + */ + if (nflt > npts) + hwidth = (int) ( npts/2. ) - 1; + else + hwidth = (int) (0.5 + nflt/2. ); + nwin = (int) 2*hwidth+1; + if (nwin < 1) nwin = 1; + win = (int *) cf_calloc(nwin, sizeof(int)); + /* + * Construct the median filter. + */ + for (i = 0; i < npts; i++) { + int j, nmin, nmax; + /* + * Determine the limits to the sample. + */ + nmin = i - hwidth; + if (nmin < 0) + nmin=0; + nmax = nmin+nwin; + if (nmax > npts-1 ) + nmin = npts-nwin-1; + + /* Fill the sample array */ + nsam=0 ; + for (j=0; j<=2*hwidth; j++) { + int ndx = nmin+j; + if (ndx > npts-1) + ndx = npts-1; + win[j] = array[ndx]; + if (array[ndx] > 0) nsam++ ; + } + + /* Get the median value for the sample array */ + medval = median(win, nwin, cntmax); + + /* + cf_verbose(6,"at point %d, nmin=%d, nmax=%d, nsam=%d, medval=%d", + i, nmin, nmax, nsam, medval) ; + */ + + /* Update the median filter */ + medflt[i] = medval; + } + + free(win) ; + +} + +/* + * Function to take an array marked with the times and intensities of + * the bursts and identify individual bursts events and output these + * to an ASCII file for later analysis. + */ + +static void +tab_burst(int *bursts, int npts, int nbin, int mjd, long tst, char *outfile) +{ + FILE *output; + int i, bflg=-1, bdur=0, bcnts=0, bst=0; + /* + * Go through the bursts array and determine the times and + * properties of the bursts. + */ + output = fopen(outfile, "w"); + + cf_verbose(2, "BURST SUMMARY") ; + for (i=0; i<npts; i++) { + if (bursts[i] > 0 && bflg < 0) { + bflg=1; + bst= i*nbin; + } + if (bursts[i] > 0 && bflg > 0) { + bdur++; + bcnts = bcnts + bursts[i]; + } + if (bursts[i] < 0 && bflg > 0) { + cf_verbose(2, "start=%d, duration=%ld seconds, counts=%d", + bst, bdur*nbin, bcnts) ; + fprintf(output, " %d %ld %d %d \n", + mjd, bst+tst, bdur*nbin, bcnts); + bflg=-1; + bdur=0; + bcnts=0; + } + } + + /* Check to see whether a burst is in progress at the end of the obs */ + if (bflg > 0) { + cf_verbose(2, "start=%d, duration=%ld seconds, counts=%d", + bst, bdur*nbin, bcnts) ; + fprintf(output, " %d %ld %d %d\n", + mjd, bst+tst, bdur*nbin, bcnts); + } + +} + + +int +cf_screen_burst(fitsfile *infits, long nevents, float *ptime, float *x, + float *y, unsigned char *locflags, GTI *gti, long nsec, + float *ttime, unsigned char *sflag, float *aic_rate, + short *bkgd) +{ + char CF_PRGM_ID[] = "cf_screen_burst"; + char CF_VER_NUM[] = "1.19"; + + short ylim[8], k; + int anynull, errflg=0, intnull=0, cflg=-1, nrowbkg=0, ngood; + int npts, nmax, nflt, iflg, numflg, ndxs, ndxe, ndx, ndxb; + int *bcnts, *cnts, *gcnts, *medflt, *flgarr, *bursts, lim_col; + int nflg, niter, mjd , status=0, medval; + long frow=1, felem=1, tst , i, j; + long stdrej, mncnt, nsmed, mnburst, nbin, nflgarr; + float exptime, max_aic_rate, bkgsf, bkgsft, srcfrac, mxtime; + double expstrt; + char burst_tab[FLEN_CARD]; + char expid[FLEN_CARD], progid[FLEN_CARD], targid[FLEN_CARD]; + char scobsid[FLEN_CARD], scrnfile[FLEN_CARD], aperf[FLEN_CARD]; + char det[FLEN_CARD], bchrfile[FLEN_CARD], aper[5]; + fitsfile *scrnfits, *bchrfits; + + /* + * Initialize error checking + */ + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing"); + if ((errflg = cf_proc_check(infits, CF_PRGM_ID))) return errflg; + + /* + * Read header keywords. + */ + FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL, &status); + FITS_read_key(infits, TDOUBLE, "EXPSTART", &expstrt, NULL, &status); + FITS_read_key(infits, TSTRING, "PRGRM_ID", progid, NULL, &status); + FITS_read_key(infits, TSTRING, "TARG_ID", targid, NULL, &status); + FITS_read_key(infits, TSTRING, "SCOBS_ID", scobsid, NULL, &status); + FITS_read_key(infits, TSTRING, "EXP_ID", expid, NULL, &status); + FITS_read_key(infits, TSTRING, "DETECTOR", det, NULL, &status); + + /* + * Exit if EXPTIME is 0. + */ + if (exptime < 1) { + cf_verbose(1, "EXPTIME = %f. Exiting.", exptime); + return status; + } + + /* + * Determine the Julian date and start time from expstrt. + */ + mjd = (int) expstrt; + tst = (int) ((expstrt - mjd)*86400.); + det[1] = tolower(det[1]); + + /* Ignore events with arrival times greater than MAX_EXPTIME. */ + i = nevents - 1; + while (ptime[i] > MAX_EXPTIME && i > 0) i--; + nevents = i + 1; + i = nsec - 1; + while (ttime[i] > MAX_EXPTIME && i > 0) i--; + nsec = i + 1; + + /* + * Determine the maximum count rate for the exposure. + * If the count rate is too high, there is the possibility of + * losing counts and of dropouts - Check the count rate and set + * a flag if it is too high. + */ + max_aic_rate = -999.; + for (i = 0; i < nsec; i++) + if (aic_rate[i] > max_aic_rate) max_aic_rate = aic_rate[i]; + cf_verbose(2, "Max AIC count rate = %d", cf_nint(max_aic_rate)); + if(max_aic_rate > 4000) cflg=1; + + /* + * Read the screening parameter file. + */ + FITS_read_key(infits, TSTRING, "SCRN_CAL", scrnfile, NULL, &status); + cf_verbose(3, "Screening parameters from %s\n", scrnfile); + FITS_open_file(&scrnfits, cf_parm_file(scrnfile), READONLY, &status); + FITS_read_key(scrnfits, TLONG, "MNCNT", &mncnt, NULL, &status); + FITS_read_key(scrnfits, TLONG, "STDREJ", &stdrej, NULL, &status); + FITS_read_key(scrnfits, TLONG, "NBIN", &nbin, NULL, &status); + /* + * Set the binning factor to 1 for high count rate data. + */ + if (cflg > 0) { + nbin=1; + cf_verbose(1, "High count-rate observation: using 1 sec time bins") ; + } + FITS_read_key(scrnfits,TLONG,"NSMED",&nsmed, NULL, &status); + /* + * srcfrac=required fraction of the source intensity for a + * burst to be counted, i.e. the burst must exceed this fraction + * before it is deemed to be significant. + */ + FITS_read_key(scrnfits, TFLOAT, "SRCFRAC", &srcfrac, NULL, &status); + FITS_read_key(infits, TSTRING, "APERTURE", aperf, NULL, &status); + FITS_close_file(scrnfits, &status); + + cf_verbose(3, "Screening parameters used:"); + cf_verbose(3, "Minimum count rate = %ld", mncnt); + cf_verbose(3, "Detection criterion (num sigma) = %ld", stdrej); + cf_verbose(3, "Number of seconds to bin the time sequence = %ld", nbin); + cf_verbose(3, "Interval (seconds) used in median filter = %ld", nsmed); + cf_verbose(3, "Fraction of source intensity for burst = %f\n", srcfrac); + /* + * Select the background extraction regions - RFPT implies + * that there is no source, so the entire detector can be used. + */ + strncpy(aper, aperf, 5); + cf_verbose(3, "Aperture designation = %s", aper); + if (strcmp ( aper , "RFPT") == 0 ) { + cf_verbose(3, "APERTURE = RFPT. No source in aperture."); + ylim[0] = 0; + ylim[1] = 1; + ylim[2] = 2; + ylim[3] = 1020; + ylim[4] = 1021; + ylim[5] = 1022; + ylim[6] = 1022; + ylim[7] = 1022; + } + else { + /* + * open the bchr parameter file + */ + FITS_read_key(infits, TSTRING, "BCHR_CAL", bchrfile, NULL, &status); + cf_verbose(3, "BCHR parameter file: %s", bchrfile); + FITS_open_file(&bchrfits, cf_parm_file(bchrfile), READONLY, &status); + + /* Move to the first extension of the BCHR file, which contains + the limits to the background regions */ + FITS_movabs_hdu(bchrfits, 2, NULL, &status); + + /* Get the column number for the data to be extracted */ + FITS_get_colnum(bchrfits, CASEINSEN, aper, &lim_col, &status); + + /* Read the data from the table */ + FITS_read_col(bchrfits, TSHORT, lim_col, frow, felem, 8, &intnull, + ylim, &anynull, &status); + + /* Close the bchr file */ + FITS_close_file(bchrfits, &status); + } + cf_verbose(3, "Y limits of background region: %d-%d, %d-%d, %d-%d\n", + ylim[0], ylim[1], ylim[2], ylim[3], ylim[4], ylim[5]); + /* + * Now get the time sequence and a good time interval flag array. + */ + cf_verbose(3, "nevents = %ld", nevents); + mxtime=ptime[nevents-1]; + nflgarr = cf_nint(mxtime); + npts= (int) mxtime/nbin; + cf_verbose(3, "Number of points in the (binned) time sequence = %d", npts); + cf_verbose(3, "Number of points in the good time flag array = %ld", nflgarr); + /* + * Generate arrays to contain the time sequence in the + * background sample regions (cnts) as well as the entire + * detector (gcnts). Also generate an array to contain only + * the bursts (bursts). + * Array bcnts will contain the background count rate. + */ + bcnts = (int *) cf_calloc(npts, sizeof(int)); + cnts = (int *) cf_calloc(npts, sizeof(int)); + gcnts = (int *) cf_calloc(npts, sizeof(int)); + bursts = (int *) cf_calloc(npts, sizeof(int)); + /* Generate an array to contain the good time flags - this is an array + with 1s time intervals. Each interval is set to a negative number (-1) + if it contains a bad time interval or a positive number (1) if it + contains a good time interval. The contents of this array are used to + screen the photons and to determine the good time intervals */ + flgarr = (int*) cf_calloc(nflgarr, sizeof(int) ); + /* Generate an array with 1s time intervals to contain the day/night + information - this will be 1 for night and 0 for day */ + + /* Initially flag all of the points with a negative number. + The good time intervals will be determined later */ + for (i=0; i<npts; i++) { + cnts[i] = -10 ; + gcnts[i] = 0 ; + bursts[i] = -10 ; + } + for (i=0; i<nflgarr; i++) flgarr[i] = -1 ; + /* + * Set the counts array to zero (i.e. unflag the point) for all + * good time points. + */ + ngood = 0 ; + for (i=0; i < gti->ntimes; i++) { + ndxs= (int) gti->start[i]/nbin; + if (ndxs < 0 ) ndxs = 0 ; + ndxe = (int) (0.5 + gti->stop[i]/nbin); + if (ndxe > npts-1) ndxe = npts-1; + for (j=ndxs; j<= ndxe; j++) { + ngood++; + cnts[j]=0; + } + } + cf_verbose(2, "Number of good time points = %d\n", ngood); + + /* Set the good time flags array to a positive value for all good + time points */ + + for (i=0; i < gti->ntimes; i++) { + ndxs= (int) gti->start[i]; + if (ndxs < 0) ndxs = 0 ; + ndxe = (int) gti->stop[i]; + if (ndxe > nflgarr-1) ndxe=nflgarr-1; + for (j=ndxs; j<= ndxe; j++) + flgarr[j]=1; + } + + /* Check the timeline-table screening flag and mark all times that have + a screening bit set (except for day/night) */ + for (i=0; i<nsec; i++) { + ndx = (int) (0.5+ttime[i]) ; + if( ndx > nflgarr-1) ndx = nflgarr-1; + if( ndx < 0) ndx = 0 ; + ndxb = (int) (0.5 + (float) ndx/ (float) nbin) ; + if(ndxb > npts-1) ndxb = npts-1 ; + if (ndxb < 0) ndxb=0 ; + if (sflag[i] & ~TEMPORAL_DAY) { + flgarr[ndx] = -1 ; + cnts[ndxb] = -10 ; + } + } + + /* Generate the time sequence for the total counts from the detector. + Use only the active area and ignore geocoronal emission. */ + + for (i=0; i < nevents; i+=1) + if( (locflags[i] & LOCATION_SHLD) == 0 && + (locflags[i] & LOCATION_AIR) == 0 ) { + int ndx = (int) ptime[i]/nbin; + if (ndx > npts-1) + ndx=npts-1; + if (ndx < 0) + ndx=0; + gcnts[ndx] += 1; + } + + /* Determine the number of rows in the background sample region */ + for (k=0; k < 3; k++) + nrowbkg += ylim[2*k+1] - ylim[2*k]; + /* + * Determine the scaling factor needed to convert the measured + * counts in the background regions to an estimate of the + * background counts over the active aperture. + */ + bkgsf= (float) ylim[6]/nrowbkg; + bkgsft = (float) ylim[7]/nrowbkg; + + cf_verbose(2, "Rows in target aperture = %d", ylim[6]); + cf_verbose(2, "Rows on the detector = %d", ylim[7]); + cf_verbose(2, "Rows in background region = %d", nrowbkg); + cf_verbose(2, "Scale factor to aperture = %f ", bkgsf); + cf_verbose(2, "Scale factor to entire detector = %f", bkgsft); + + /* Generate the time sequence using only the counts from the three + selected background regions - ignore geocoronal lines and times + that have been flagged as bad (with negative count values) */ + for (k=0; k<3; k++) { + short ylow=ylim[2*k]; + short yhigh=ylim[2*k+1]; + for (i=0; i < nevents; i+=1) { + if( (locflags[i] & LOCATION_SHLD) == 0 && + (locflags[i] & LOCATION_AIR) == 0 && + (y[i] > ylow) && (y[i] < yhigh) ) { + int ndx = (int) ptime[i]/nbin; + if (ndx > npts-1) + ndx=npts-1; + if (ndx < 0) + ndx=0; + if (cnts[ndx] >= 0) cnts[ndx] += 1; + bcnts[ndx] ++; + } + } + } + + /* Fill background count rate array with cnts array, and divide by nbin + to give count rate */ +/* THIS SCHEME WORKS ONLY IF THERE ARE NO TIME BREAKS IN THE TIMELINE + * TABLE. SINCE TIME BREAKS ARE COMMON, MUST USE TTIME ARRAY TO MAP + * CNTS ARRAY TO BKGD ARRAY. WVD 09/15/03 + * for (i=0; i<npts; i++) { + * int ndxs= i*nbin; + * int ndxe= ndxs + nbin; + * if (ndxe > nsec-1) + * ndxe=nsec-1; + * for (j=ndxs;j<ndxe; j++) + * bkgd[j] = (short) ((float) cnts[i] / nbin + 0.5); + * } + */ + for (i=0; i<nsec; i++) { + ndx = cf_nint(ttime[i]/nbin) ; + if( ndx > nflgarr-1) + ndx = nflgarr-1; +/* if (cnts[ndx] > 0.) + bkgd[i] = (short) ((float) cnts[ndx] / nbin + 0.5); + else + bkgd[i] = 0; +*/ + bkgd[i] = bcnts[ndx] / nbin; + } + free (bcnts); + + /* Remove the background counts (scaled to the entire detector) from the + measured gross count rate */ + for (i=0; i < npts; i++) if (cnts[i] > 0) { + gcnts[i] = gcnts[i] - bkgsft*cnts[i]; + if (gcnts[i] < 0) + gcnts[i] = 0; + } + + /* Check for zeros in the count array for high count rate data - these + indicate dropouts. If one is found, flag it with a negative number so + that the program will ignore it in future analysis */ + numflg=0; + if (cflg > 0) { + for (i=0; i < npts; i++) { + if (cnts[i] <= 0.1 ) { + numflg += 1; + cnts[i]=-10; + } + } + } + cf_verbose(2, "Number of dropouts found = %d", numflg); + + /* Determine the median value of the array */ + nmax=2000; + medval = median(cnts, npts, nmax); + cf_verbose(2, "Median value of the entire original array = %d ", medval); + + /* Screen the time series for very large count rates (> 3 x median) */ + cf_verbose(2, "Searching for large bursts ") ; + for ( i=0; i<npts; i++) + if (cnts[i] > 4*medval) { + cf_verbose(2,"point rej: time=%d, cnts=%d, 4*medval=%d ", + i*nbin, cnts[i], 4*medval) ; + bursts[i] = cnts[i]-medval; + cnts[i] = -100; + } + + /* Do a median filter on the pre-filtered array */ + + /* Generate the array to contain the median filter */ + medflt = (int *) cf_calloc(npts, sizeof(int) ); + + /* + * Now specify the parameters of the filter - + * nsmed = number of seconds over which to determine median values + * nflt = number of points in the cnts array over which to + * determine the median value + */ + nflt = nsmed / nbin; + niter=0; + + cf_verbose(3,"Searching for small bursts - nflt = %d", nflt) ; + + do { + median_filter(cnts, npts, nflt, medflt); + + /* Throw out all points which are more than (stdrej) standard + deviations above the median. Also require that the counts exceed + a certain minimum value. */ + nflg=0; + mnburst=mncnt*nbin; + for (i=0; i<npts; i++) { + int dcnt=cnts[i]-medflt[i]; + if ( (dcnt > mnburst) && + (dcnt > stdrej*sqrt((double) cnts[i])) && + (dcnt*bkgsf > srcfrac*gcnts[i]) ) { + nflg++; + cf_verbose(2, "point rej: time=%ld, cnts=%d, medflt=%d, " + "dif=%d, gcnts= %d", + i*nbin, cnts[i], medflt[i], dcnt, gcnts[i]); + cnts[i] = -100; + bursts[i] = dcnt; + } + } + niter++; + cf_verbose(3,"After iteration %d, we found %d times affected by bursts", + niter, nflg) ; + /* Iterate the median filter if any additional points have + * been modified*/ + } while(nflg > 0 && niter < 10); + + /* Print out the results of the screening + cf_verbose(5, " "); + for (i=0; i<npts; i++) + cf_verbose(5, "At time %4d, cnts= %4d, medflt= %4d, dif= %4d, gcnts= %4d", + i*nbin, cnts[i],medflt[i], cnts[i]-medflt[i], gcnts[i]); + cf_verbose(3, " "); + */ + + /* Use the regions of the cnts array equal to -100 to specify + the burst times in the flgarr array */ + for (i=0; i<npts; i++ ) { + if (cnts[i] < -20) { + int ndxs = i*nbin, ndxe; + if (ndxs > nflgarr-nbin-2) + ndxs=nflgarr-nbin-2; + ndxe = ndxs + nbin - 1; + for (j=ndxs; j<=ndxe; j++) + flgarr[j] = -10; + } + } + + /* tabulate the burst times */ + sprintf(burst_tab, "%4s%2s%2s%3s%2s_bursts.dat", progid, targid, scobsid, expid, det); + tab_burst(bursts, npts, nbin, mjd, tst, burst_tab); + cf_verbose(3, "Burst properties written to file: %s", burst_tab); + + /* Check the last 2*nbin elements of flgarr. If any of these are set as + bad, then set all the end points as bad. The binned count time sequence + will not sample the end of the ovarall time sequence if there is not an + integral number of bins present. We also do not need the last bin if the + second-to-last is bad. */ + + i=nflgarr-2*nbin; + iflg=10; + do{ + if (flgarr[i] < -1) { /* 09/12/03 - wvd - Change from 0 */ + for (j=i; j<nflgarr; j++) + flgarr[j]=-10; + iflg = -10; + }; + i++; + } while (iflg > 0 && i < nflgarr); + + /* set the timeline status flag (sflag - TEMPORAL_BRST) for all + times containing bursts */ + for (i=0; i<nsec; i++) { + ndx = (int) (0.5+ttime[i]) ; + if( ndx > nflgarr-1) + ndx = nflgarr-1; + if (flgarr[ndx] < -1) /* 09/12/03 - wvd - Change from 0 */ + sflag[i] = sflag[i] | TEMPORAL_BRST; + } + /* deallocate storage space */ + free(cnts) ; + free(gcnts) ; + free(flgarr) ; + free(medflt) ; + free(bursts) ; + + /* Enter a time stamp into the log */ + cf_proc_update(infits, CF_PRGM_ID, "COMPLETE"); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished processing"); + + return (status); +} |