/**************************************************************************** * 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 #include #include #include #include #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 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; intimes; 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 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 nsec-1) * ndxe=nsec-1; * for (j=ndxs;j 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 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 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 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 0 && i < nflgarr); /* set the timeline status flag (sflag - TEMPORAL_BRST) for all times containing bursts */ for (i=0; i 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); }