#! /usr/bin/env python -tt -u #-3 #"""Run basic GRB processing on Fermi/GBM data. # # This script runs a simple GRB analysis on a set of Fermi/GBM event # files (.fit) for a given trigger. The output is one FITS light # curve (.lc) file and one FITS PHA2 (.pha) file for each input # detector event file. A spectral fit is performed. # #""" import math import os import re import subprocess import sys import numpy import pyfits import xspec #### def get_times(time_string): """Split a string of times into background and source intervals. Usage: times = get_times(time_string) Input: time_string -- (string) b1-b2:s1-s2:b3-b4 where b1,b2 are the stop and start times of the first background interval; s1,s2 are the start and stop times of the source interval; and b3,b4 are the start and stop times of the second background interval. Output: times -- (list) A list containing the six times that define the three intervals. Returns None if time_string cannot be parsed into six floats. """ match = re.search(r'^(\S+)-(\S+):(\S+)-(\S+):(\S+)-(\S+)$', time_string, re.IGNORECASE) if match: print 'get_times: {0} is a time string'.format(time_string) times = [] for i in range(6): t = match.group(i+1) if is_number(t): times.append( float(t) ) else: print 'get_times: error in time string {0}, {1} is a bad time'.format(time_string, t) return None else: return None if ( times[1] <= times[0] ) or \ ( times[3] <= times[2] ) or \ ( times[5] <= times[4] ): print 'get_times: times not increasing in {0}'.format(time_string) return 'ERROR' return times #### def is_number(s): """Check if a string is a number.""" try: float(s) except ValueError: return False else: return True #### def get_detectors(detector_string, trigger_name): """Parse the detectors string. Usage: detector_files = get_detectors(detector_string, trigger_name) Input: detector_string -- (string) ALL - Use all the detectors that had a detection. TRIGGER - Use the detectors that triggered. nnnnnnnnnnnn - Use detectors specified by the user. trigger_name -- (string) Fermi/GBM trigger id. Output: files -- (list) List of detector files, or None if none were found If detector_string = TRIGGER this function reads the DET_MASK keyword from the glg_tcat_all file to identify the triggered detectors. """ detector_files = get_list_of_detector_files(trigger_name) if detector_string == 'ALL': files = detector_files return files elif detector_string == 'TRIGGER': # Find the tcat file. working_directory = os.getcwd() files = [] for file in os.listdir(working_directory): if re.match(r'^glg_tcat_all_' + trigger_name + '_v\d+.fit', file, re.IGNORECASE): tcat_file = file break # Read which detectors triggered fits = pyfits.open(tcat_file, 'readonly') detmask = fits[0].header['DET_MASK'] fits.close() files = get_specified_detectors(trigger_name, detector_files, detmask) return files elif re.match(r'\d{12}', detector_string): detmask = detector_string files = get_specified_detectors(trigger_name, detector_files, detmask) return files else: return None #### def get_specified_detectors(trigger_name, detector_files, detmask): """Create a list of detector files to use. Usage: files = get_specified_detectors(trigger_name, detector_files, detmask) Input: trigger_name -- (string) Fermi/GBM trigger id. detector_files -- (list) List of detector files to check. detmask -- (string) Which detectors to use. Output: files -- (list) List of detector files to use. """ # Build the list of detectors. list = [] match = re.search(r'(\d)(\d)(\d)(\d)(\d)(\d)(\d)(\d)(\d)(\d)(\d)(\d)', detmask, re.IGNORECASE) # Clumsy! for i in range(12): bit = int( match.group(i+1) ) if bit == 1: if ( i < 10 ): det = 'n' + str(i) elif i == 10: det = 'na' elif i == 11: det = 'nb' else: print 'warning: get_detectors: bad detector number {0}'.format(i) det = '' list.append(det) files = [] for file in detector_files: for det in list: if re.match(r'.+_' + det + '.+.fit', file): files.append(file) return files #### def get_list_of_detector_files(trigger_name): """Return Fermi/GBM detector files in the current directory. Usage: files = get_list_of_detector_files(trigger_name) Input: trigger_name -- (string) Fermi/GBM trigger id Output: files -- (list) List of detector files """ working_directory = os.getcwd() files = [] for detector_file in os.listdir(working_directory): if re.match(r'^glg_' + '\w+' + '_\S{2}_' + trigger_name + '_v\d+.fit', detector_file, re.IGNORECASE): files.append(detector_file) return files #### def run_gtbin_lightcurve(detector_file, lightcurve_file, tstart, tstop, dtime): """Call gtbin to generate a Fermi/GBM light curve. Usage: run_gtbin_lightcurve(detector_file, lightcurve_file, tstart, tstop, dtime) Input: detector_file -- (string) Fermi/GBM detector file name lightcurve_file -- (string) Output lightcurve file name tstart -- (float) Start time for lightcurve in MET s. tstop -- (float) Stop time for lightcurve in MET s. dtime -- (float) Bin width for lightcurve in s Output: Returns None if gtbin ran without generating an error. Returns an error string if there was a problem. """ command = ['gtbin', 'evfile=' + detector_file, 'scfile=NONE', 'outfile=' + lightcurve_file, 'algorithm=LC', 'tbinalg=LIN', 'tstart=' + str(tstart), 'tstop=' + str(tstop), 'dtime=' + str(dtime), 'chatter=2', 'clobber=YES'] try: subprocess.check_call(command) except: print 'warning, subprocess call failed for {0}'.format(command) return 'gtbin failed' else: print 'created {0} (dtime={1})'.format(lightcurve_file, dtime) return None #### def time_intervals(lightcurve_file, trigtime): """Determine background and source intervals. Usage: intervals = time_intervals(lightcurve_file, trigtime) Input: lightcurve_file -- (string) lightcurve file name trigtime -- (float) trigger time of the burst Output: intervals -- (list) The start and stop times, in MET, of the three intervals (first background, source, second background. Returns None if no intervals were found This function finds preliminary intervals using the input light curve then uses the two background intervals to fit and subtract the background. It then finds refined intervals using the background-subtracted light curve. """ # Determine the preliminary background and source intervals intervals = get_time_intervals(lightcurve_file) if intervals is None: print 'warning, no source found in {0}, skipping'.format(lightcurve_file) return None # Fit and subtract the background from the preliminary light curve flat_lightcurve_file = subtract_background_curvature(lightcurve_file, trigtime, intervals) if flat_lightcurve_file is None: print 'warning, unable to flatten the lightcurve for {0}'.format(lightcurve_file) return None # Determine the background and source intervals from the # background-subtracted light curve. intervals = get_time_intervals(flat_lightcurve_file) if intervals is None: print 'warning, no source found in {0}, skipping'.format(flat_lightcurve_file) return None return intervals #### def fit_background(time, counts, trigtime, intervals): """Fit the background for a light curve. Usage: background = fit_background(time, counts, trigtime, intervals) Input: time -- (list) time in MET seconds counts -- (list) counts corresponding to time trigtime -- (float) trigger time of burst intervals -- (list) the background and source time intervals Output: background -- (list) the fitted background corresonding to time This function fits a polynomial to the two background intervals. The source interval is not used in the fit. """ polynomial_order = 2 bkg1_start = intervals[0] bkg1_stop = intervals[1] src_start = intervals[2] src_stop = intervals[3] bkg2_start = intervals[4] bkg2_stop = intervals[5] # Subtract the trigger time. t = [] for i in range(len(time)): t.append( time[i] - trigtime ) # Create arrays with the source removed tt = [] cc = [] for i in range(len(time)): if ( time[i] >= bkg1_start ) and ( time[i] < bkg1_stop ): tt.append(t[i]) cc.append(counts[i]) if ( time[i] >= bkg2_start ) and ( time[i] < bkg2_stop ): tt.append(t[i]) cc.append(counts[i]) # Fit the background coeffs = numpy.polyfit(tt, cc, polynomial_order) # Create the background array background = [] for i in range( len(t) ): bkg = 0.0 for j in range( polynomial_order + 1 ): power = polynomial_order - j bkg += coeffs[j] * t[i]**power background.append(bkg) return background #### def subtract_background_curvature(lightcurve, trigtime, intervals): """Fit and subtract curvature from a Fermi/GBM lightcurve. Usage: new_lightcurve = subtract_background_curvature(lightcurve, trigtime, intervals) Input: lightcurve -- (string) lightcurve filename trigtime -- (float) burst trigger time intervals -- (list) source and background intervals Output new_lightcurve -- (string) background-subtracted lightcurve file name Read the light curve, then fit and subtract the background. The mean background is added back to the light curve to preserve Poisson statistics. The final light curve is written to a new light curve file. """ # Get the light curve lc_file = pyfits.open(lightcurve, 'readonly') phu = lc_file[0] rate_header = lc_file['RATE'].header rate_data = lc_file['RATE'].data gti = lc_file['GTI'] time = lc_file['RATE'].data.field('TIME') counts = lc_file['RATE'].data.field('COUNTS') error = lc_file['RATE'].data.field('ERROR') lc_file.close() # The last data point can be incomplete, so set it equal to the # second-to-last data point. n = len(time) - 1 counts[n] = counts[n-1] error[n] = error[n-1] # Fit the background background = fit_background(time, counts, trigtime, intervals) # Compute the mean background and standard deviation sum = 0.0 for i in range( len(background) ): sum += background[i] mean_bkg = sum / len(background) sum = 0.0 for i in range( len(background) ): term = (background[i] - mean_bkg)**2 sum += math.sqrt(term) n = len(background) var_bkg = sum / ( n - 1 ) sd_bkg = math.sqrt(var_bkg) # Subtract the fitted trend from each time bin bkg_subtracted_counts = [] for i in range( len(counts) ): src = counts[i] - background[i] + mean_bkg bkg_subtracted_counts.append(src) # Create a new HDU with trend-subtracted counts. # Bewarned and beware: get_coldefs() is depreciated in # pyfits 3. #cdefs = rate_data.columns # pyfits 3.0 cdefs = lc_file['RATE'].get_coldefs() # pyfits 2.1.1 #print cdefs.info # Bewarned and beware: there is a bug in pyfits 2.2.2 and earlier # that causes del_col to raise a KeyError exception when deleting # a column from a ColDefs object. See # . #cdefs.del_col('COUNTS') # pyfits 3.0 #cdefs.change_name('COUNTS', 'OCOUNTS') # pyfits 2.1.1 for i in range( len(counts) ): rate_data.field('COUNTS')[i] = bkg_subtracted_counts[i] #print cdefs.info #col = pyfits.Column(name='COUNTS', format='J', unit='Counts', # array=bkg_subtracted_counts) # Bewarned and beware: the behaviour of add_col changed in pyfits # 2.3. #cdefs.add_col(col) # pyfits 2.3 #new_rate_hdu = pyfits.new_table(cdefs, header=rate_header) #new_cdefs = cdefs.add_col(col) # before pyfits 2.3 #print cdefs.info #new_rate_hdu = pyfits.new_table(new_cdefs, header=rate_header) # Write the trend-subtracted light curve to a new FITS file. #lc_file['RATE'] = new_rate_hdu lc_file['RATE'].data = rate_data #lc_file.info() new_lightcurve = lightcurve.replace('.lc', '_flatbkg.lc') # Update keywords lc_file[0].header['FILENAME'] = new_lightcurve lc_file[0].header['CREATOR'] = 'do_gbm.py' lc_file['RATE'].header['CREATOR'] = 'do_gbm.py' if os.path.exists(new_lightcurve): os.remove(new_lightcurve) try: lc_file.writeto(new_lightcurve) except ValueError: #print 'warning, apparently meaningless python ValueError exception when writing to {0} is being ignored in {1}'.format(new_lightcurve, __name__) pass except: print 'subtract_background_curvature: error writing new lightcurve file {0}'.format(new_lightcurve) return None finally: return new_lightcurve #### def get_time_intervals(event_file): """Return background and source start and stop times intervals. Usage: (bkg1_start,bkg2_stop,src_start,src_stop bkg2_start,bkg2_stop) = get_time_intervals(event_file) Input: event_file -- (string) Fermi/GBM detector event file name Output: bkg1_start -- (float) background interval 1 start time (MET) bkg1_stop -- (float) background interval 1 stop time (MET) src_start -- (float) source interval start time (MET) src_stop -- (float) source interval stop time (MET) bkg2_start -- (float) background interval 2 start time (MET) bkg2_stop -- (float) background interval 2 stop time (MET) Call gtburstfit to determine Bayesian blocks for the lightcurve. The first background interval is the first Bayesian block. The last time interval is the last Bayesian block unless it is shorter than min_block_duration, then it is the last two Bayesian blocks. The source interval is the interval between the two background intervals. """ # This is the minimum allowed duration of the second background # region. min_block_duration = 1.0 # s # Just a guess! log_file = open('TT_gtburstfit.LOG', 'w') command = ['gtburstfit', 'evfile=' + event_file, 'fitguess=AUTO', 'amp=0.0', 'time0=0.0', 'tau1=0.0', 'tau2=0.0', 'bckgnd=0.0', 'evtable=RATE', 'fit=YES', 'plotres=NO', 'ncpprior=9.0', 'plot=NO', 'title=DEFAULT', 'chatter=2', 'clobber=YES', 'debug=NO', 'gui=NO'] try: subprocess.check_call(command, stdout=log_file) except: print 'warning, unable to fit the light curve for: {0}'.format(event_file) return None finally: log_file.close() log_file = open('TT_gtburstfit.LOG', 'r') blocks = [] i = 0 found_bayesian_blocks = False for line in log_file.readlines(): if not found_bayesian_blocks: if 'Bayesian Blocks' in line: found_bayesian_blocks = True else: m = re.match(r'\[(\S+),\s*(\S+)\]', line, re.IGNORECASE) if m is not None: blocks.append( (m.group(1), m.group(2)) ) i += 1 log_file.close() i -= 1 # The definitions of the background interval and source interval # are the same as those used by the Swift/BAT tool battblocks. bkg1_start = float( blocks[0][0] ) # First Bayesian block is the bkg1_stop = float( blocks[0][1] ) # first background interval src_start = float( blocks[0][1] ) # Source interval is from end of 1st src_stop = float( blocks[i][0] ) # Bayesian block to start of last. # If the last Bayesian block is less than some minimum duration # then use the second-to-last Bayesian block to determine the stop # time of the source. last_block_duration = float(blocks[i][1]) - float(blocks[i][0]) if ( last_block_duration < min_block_duration ): src_stop = float(blocks[i-1][0]) # Set the second background interval to be from the source stop # time to the end of the last Bayesian block. bkg2_start = src_stop bkg2_stop = float( blocks[i][1] ) return [bkg1_start, bkg1_stop, src_start, src_stop, bkg2_start, bkg2_stop] #### def compute_signal_to_noise(lightcurve_file, intervals, user_times, trigtime): """Compute the signal-to-noise ratio of a lightcurve. Usage: sn = compute_signal_to_noise(lightcurve_file, intervals user_times, trigtime) Input: lightcurve_file -- (string) name of the lightcurve file intervals -- (list) start and stop times of the source and background intervals user_times -- (boolean) True if the intervals were set by the user and False if the intervals were determined by the script trigtime -- (float) trigger time (MET s) Output: mean_sn -- (float) the mean sigmal-to-noise ratio per bin in the source interval intervals -- (list) improved source and background intervals This function reads the lightcurve file then fits and subtracts the background. The signal-to-noise ratio is computed in each bin in the source interval then a mean signal-to-noise ratio is computed. If user_times is False (the user did not specify the intervals) then the source interval is refined so that it ends when max_consecutive_negative_sn bins with negative signal-to- noise ratios are encountered. The second backgroud region is also adjusted. """ # This is the maximum number of consecutive bins that can have # negative signal-to-noise ratios before the source is assumed to # have ended. This parameter is ignored if the user specified # the source interval. max_consecutive_negative_sn = 3 bkg1_start = intervals[0] bkg1_stop = intervals[1] src_start = intervals[2] src_stop = intervals[3] bkg2_start = intervals[4] bkg2_stop = intervals[5] # Get the light curve f = pyfits.open(lightcurve_file, 'readonly') time = f['RATE'].data.field('TIME') counts = f['RATE'].data.field('COUNTS') error = f['RATE'].data.field('ERROR') f.close() # Compute the background counts background = fit_background(time, counts, trigtime, intervals) # Compute the signal-to-noise spectrum n = 0 sum = 0.0 src_end_time = src_stop # Skip the last time bin because it may be truncated. consecutive_negative_sn = 0 for i in range( len(time) - 1 ): if ( time[i] >= src_start and time[i] <= src_stop ): sn = ( counts[i] - background[i] ) / error[i] if sn < 0.0: consecutive_negative_sn += 1 if consecutive_negative_sn > max_consecutive_negative_sn: dt = src_stop - time[i] print 'warning, source appears to have stopped {0} s before the fitted end of the burst ({1} s)'.format(dt, src_stop) src_end_time = time[i] if not user_times: print 'adjusting the source stop time to {0} (MET)'.format(src_end_time) break # The source has ended, stop summing n += 1 sum += sn elif time[i] > src_stop: break mean_src_sn = sum / n if mean_src_sn < 0.0: print 'warning, source S/N < 0 (S/N = {0}'.format(mean_src_sn) print ' sum={0} n={1}'.format(sum, n) if not user_times: intervals[3] = src_end_time # Revised source end time intervals[4] = src_end_time # Revise 2nd bkgd interval start time return [mean_src_sn, intervals] #### def compute_duration(lightcurve_file, src_start, src_stop, bkg_start, bkg_stop, fraction): """Compute the T_fraction duration of the source. Usage: [t_fractiont_fraction, t_fraction_err, t_fraction_start, t_fraction_stop] = compute_duration(lightcurve_file, src_start, src_stop, bkg_start, bkg_stop, fraction) Input: lightcurve_file -- (string) Name of lightcurve file src_start -- (float) Start time of source (MET) src_stop -- (float) Stop time of source (MET) bkg_start -- (float) Start time of background (MET) bkg_stop -- (float) Stop time of background (MET) fraction -- (float) T_fraction value (0-100) Output: t_fraction -- (float) T_fraction duration (s) t_fraction_err -- (float) T_fraction 1-sigma error (s) t_fraction_start -- (float) Start time of T_fraction interval (MET s) t_fraction_stop -- (float) Stop time of T_fraction interval (MET s) Computes the T_fraction (e.g. T_90, T_50) duration of a source as well as its 1-sigma error and the start and stop times of the T_fraction duration. This function uses the method described in the battblocks documentation and the TOTVAR method to compute the error. """ # Get the light curve f = pyfits.open(lightcurve_file, 'readonly') time = f['RATE'].data.field('TIME') counts = f['RATE'].data.field('COUNTS') error = f['RATE'].data.field('ERROR') f.close() # Compute the mean background counts n = 0 sum = 0.0 for i in range( len(time) ): if ( time[i] >= bkg_start and time[i] <= bkg_stop ): n += 1 sum += counts[i] bkg = sum / n if n <= 0: print 'error: no counts in specified background interval' # Compute the total source counts total_src_counts = 0.0 for i in range( len(time) ): if ( time[i] >= src_start and time[i] <= src_stop ): total_src_counts += counts[i] - bkg # Compute the cumulative fraction light curve sum = 0.0 cumulative_lc = [] f = [] f_error = [] for i in range( len(time) ): if ( time[i] < src_start ): cumulative_lc.append(0.0) f.append(0.0) f_error.append(0.0) elif ( time[i] >= src_start and time[i] <= src_stop ): src_counts = counts[i] - bkg sum += src_counts cumulative_lc.append(sum) frac = cumulative_lc[i] / total_src_counts f.append(frac) e = error[i] / total_src_counts f_error.append(e) elif ( time[i] > src_stop ): cumulative_lc.append(total_src_counts) f.append(1.0) f_error.append(0.0) # Compute the T_fraction time and error offset = (100.0 - fraction ) / 2.0 / 100.0 xp = [offset, 1.0-offset] yp = numpy.interp(xp, f, time) t_fraction = yp[1] - yp[0] # Record the start and stop times of the T_fraction interval. t_fraction_start = yp[0] t_fraction_stop = yp[1] sum = 0 for i in range( len(time) ): sum += f_error[i]**2 t_fraction_err = math.sqrt(sum) return [t_fraction, t_fraction_err, t_fraction_start, t_fraction_stop] #### def select_on_energy(infile, emin, emax, outfile): """Extract a light curve in a user-specified energy band. Usage: select_on_energy(infile, emin, emax, outfile) Input: infile -- (string) name of input events file emin -- (float) minimum energy (keV) emax -- (float) maximum energy (keV) outfile -- (string) name of output events file Call fselect to select events with energies in the range [emin,emax]. """ # Get the EBOUNDS data from the events file. event_file = pyfits.open(infile, 'readonly') channel = event_file['EBOUNDS'].data.field('CHANNEL') e_min = event_file['EBOUNDS'].data.field('E_MIN') e_max = event_file['EBOUNDS'].data.field('E_MAX') event_file.close() # Find the channel bounderies corresponding to [emin,emax]. cmin = -1 cmax = -1 for i in range( len(channel) ): if ( e_min[i] <= emin ) and ( e_max[i] >= emin ): cmin = channel[i] if ( e_min[i] <= emax ) and ( e_max[i] >= emax ): cmax = channel[i] # Run fselect in_events = infile + '[EVENTS]' expression = 'PHA.GE.' + str(cmin) + ' .AND. PHA.LE.' + str(cmax) command = ['fselect', 'infile=' + in_events, 'outfile=' + outfile, 'expr=' + expression, 'histkw=YES', 'copyall=YES', 'keycopy=YES', 'clobber=YES'] try: subprocess.check_call(command) except: print 'warning, unable to extract {0} - {1} keV events from: {2}'.format(emin, emax, in_event) print 'expression used: {0}'.format(expression) return #### def create_time_bins_file(bkg_start, bkg_stop, src_start, src_stop, fits_file,): """Create a time bin file for gtbin. Usage: create_time_bin_file(background_start_time, background_stop_time, source_start_time, source_stop_time, output_file_name) Input: background_start_time -- (float) in MET s background_stop_time -- (float) in MET s source_start_time -- (float) in MET s source_stop_time -- (float) in MET s fits_file -- (string) name of output time bins file Call gtbindef to create a FITS file with the background and source intervals that will be used to extract spectra from the Fermi/GBM detector files. """ text_file = fits_file + '.txt' f = open(text_file, 'w') print >>f, bkg_start, bkg_stop print >>f, src_start, src_stop f.close() command = ['gtbindef', 'bintype=T', 'binfile=' + text_file, 'outfile=' + fits_file, 'chatter=2', 'clobber=YES', 'debug=NO', 'gui=NO'] try: subprocess.check_call(command) except: print 'warning, subprocess call failed for {0}'.format(fits_file) os.remove(text_file) return #### def find_response_file(pha_file): """Find a .rsp file corresponding a .pha file. Usage: rsp_file = find_response_file(pha_file) Input: pha_file -- (string) Name of .pha file Output: rsp_file -- (string Name of corresponding .rsp file This function finds all of the .rsp files in the current directory. Then it looks for a .rsp file with exactly the same name as the .pha except that tte is replaced with cspec. Example: glg_tte_n6_bn121205000_v01.rsp and glg_cspec_n6_bn121205000_v01.pha """ spectrum_root = pha_file[0:pha_file.rfind('.pha')-18] s_len = len(spectrum_root) s_id = pha_file[s_len:len(pha_file)] s_id = s_id[0:s_id.rfind('.pha')-4] response_root = 'glg_cspec_' r_len = len(response_root) rsp_file = None working_directory = os.getcwd() spectrum_files = [] response_files = [] for file in os.listdir(working_directory): if re.match(r'\S+\.rsp', file, re.IGNORECASE): response_files.append(file) response_file_found = False for r in response_files: if 'cspec' in r: r_id = r[r_len:len(r)] r_id = r_id[0:r_id.rfind('.rsp')-4] if r_id == s_id: response_file_found = True rsp_file = r return rsp_file #### def fit_model_spectrum(statmethod, xspec_save_file): """Fit an XSpec model to spectral data. Usage: fit_model_spectrum(statmethod, xspec_save_file) Input: statmethod -- (string) XSpec statmethod (e.g., chi, cstat) xspec_save_file -- (string) File to save XSpec commands to This function uses pyxspec to fit an XSpec model. It assumes that the data and the model have already been set up. This functions just does the actual fitting. """ # Fit a spectral model. xspec.Fit.renorm() xspec_save_file.write('renorm\n') xspec.Fit.query = 'yes' xspec_save_file.write('query yes\n') xspec.Fit.statMethod = statmethod s = 'statistic ' + statmethod + '\n' xspec_save_file.write(s) status=xspec.Fit.perform() print'status: ',status xspec_save_file.write('fit\n') return #### def fit_spectra(trigger_name, spectrum_list, response_list, fit_statistic, e_flu_min, e_flu_max, duration, plot_file): """Do an XSpec spectral fit to Fermi/GBM burst data. Usage: [model_data, [goodness_of_fit, dof], [fluence, fluence_err]] = fit)spectra(trigger_name, spectrum_list, response_list, fit_statistic e_flu_min, e_flu_max, duration, plot_file) Input: trigger_name -- (string) Fermi/GBM trigger id spectrum_list -- (list) list of .pha files to use response_list -- (list) of .rsp files to use fit_statistic -- (string) XSpec fit statistic to use e_flu_min -- (float) minimum energy for fluence calculation (keV) e_flu_max -- (float) maximum energy for fluence calculation (keV) duration -- (float) duration of the spectrum (s) plot_file -- (string) name of output plot file Output: model_data -- (list) XSpec model information [name, par1, par1_err, par2, par2_err, ...] goodnes_of_fit -- (float) goodness of fit value dof -- (float) degrees of freedom in fit fluence -- (float) fluence (erg cm^-2) fluence_err -- (float) error in fluence (erg cm^-2) This function does the spectral fitting using pyxspec. It currently tries to fit a Band model and a cutoff power law model and chooses the one with the best fit. """ # Set a log file for the output from XSpec log_file_name = trigger_name + '_xspec.log' log_file = xspec.Xset.openLog(log_file_name) # Use the log file, not the screen, for output xspec.Xset.chatter = 0 # Create a file to save XSpec commands to. xspec_save_file_name = trigger_name + '_xspec_save' + '.xcm' xspec_save_file = open(xspec_save_file_name, 'w') xspec_save_file.write('# XSpec command file\n') s = '# Created by {0}'.format(__file__) xspec_save_file.write(s) s = '# ' + trigger_name + '\n' xspec_save_file.write(s) xspec.AllData.clear() # Read in the data for i in range( len(spectrum_list) ): j = i + 1 xspec.AllData += spectrum_list[i]+'{2}' s = 'data 1:' + str(j) + ' ' + spectrum_list[i] + '{2}\n' xspec_save_file.write(s) xspec.AllData(i+1).response = response_list[i] s = 'response 1:' + str(j) + ' ' + response_list[i] + '\n' xspec_save_file.write(s) xspec.AllData(i+1).background = spectrum_list[i]+'{1}' s = 'backgrnd 1:' + str(j) + ' ' + spectrum_list[i] + '{1}\n' xspec_save_file.write(s) xspec.AllData.show() # Determine the Fermi/GBM detector type. dettype_list = [] print 'Detector Types:' for i in range( len(spectrum_list) ): spectrum_file = spectrum_list[i] fits = pyfits.open(spectrum_file, 'readonly') detnam = fits[0].header['DETNAM'] fits.close() dettype_list.append(detnam[0:detnam.find('_')]) print ' {0}: {1}'.format(spectrum_file, dettype_list[i]) # Set the good channels. xspec.Plot.xAxis = 'keV' # The ignore values are in keV when the X axis has units of energy. # Ignore values must be decimals cast to strings. for i in range( len(spectrum_list) ): if 'NAI' in dettype_list[i]: ign1 = str(0.0) ign2 = str(8.0) ign3 = str(900.0) ign4 = '**' elif 'BGO' in dettype_list[i]: ign1 = str(0.0) ign2 = str(200.0) ign3 = str(40000.0) ign4 = '**' else: print '{0}, error, unknown detector type ({1})'.format(__file__, dettype_list[i]) xspec.AllData(i+1).ignore(ign1 + '-' + ign2 + ',' + ign3 + '-' + ign4) j = i + 1 s = 'ignore ' + str(j) + ':' + ign1 + '-' + ign2 + \ ' ' + ign3 + '-' + ign4 + '\n' xspec_save_file.write(s) # Try fitting both the grbm and the cutoffpl models and see which # gives the best fit. print 'fit_spectra: fitting cflux*(grbm) model' xspec.AllModels.clear() # Load the GRB model. model = xspec.Model('cflux*(grbm)') xspec_save_file.write('model cflux*(grbm)\n') xspec_save_file.write(str(e_flu_min)+'\n') # Emin model.cflux.Emin = [e_flu_min] xspec_save_file.write(str(e_flu_max) + '\n') # Emax model.cflux.Emax = [e_flu_max] xspec_save_file.write('-6\n') # lg10Flux model.cflux.lg10Flux = [-6] xspec_save_file.write('-1\n') # alpha model.grbm.alpha = [-1] xspec_save_file.write('-2\n') # beta model.grbm.beta = [-2] xspec_save_file.write('300\n') # tem model.grbm.tem = [300] xspec_save_file.write('1,-1\n') # norm model.grbm.norm = [1] model.grbm.norm.frozen = True # Fit the GRB model. try: fit_model_spectrum(fit_statistic, xspec_save_file) xspec.Fit.error("3 4 5 6") xspec_save_file.write('error 3 4 5 6\n') fit_model_spectrum(fit_statistic, xspec_save_file) grbm_lg10Flux = model.cflux.lg10Flux.values[0] grbm_alpha = model.grbm.alpha.values[0] grbm_beta = model.grbm.beta.values[0] grbm_tem = model.grbm.tem.values[0] goodness_of_fit_grbm = xspec.Fit.statistic dof_grbm = xspec.Fit.dof xspec.Fit.error("3") xspec_save_file.write('error 3\n') grbm_lg10Flux_err = xspec.AllModels(1)(3).error xspec.Fit.error("4") xspec_save_file.write('error 4\n') err = xspec.AllModels(1)(4).error grbm_alpha_err = ( err[1] - err[0] ) /2.0 xspec.Fit.error("5") xspec_save_file.write('error 5\n') err = xspec.AllModels(1)(5).error grbm_beta_err = ( err[1] - err[0] ) / 2.0 xspec.Fit.error("6") xspec_save_file.write('error 6\n') err = xspec.AllModels(1)(6).error grbm_tem_err = ( err[1] - err[0] ) / 2.0 # Plot the results plot_file1 = 'grbm_' + plot_file xspec.Plot.device = plot_file1+'/ps' xspec.Plot('uf') except Exception, msg: print msg model.show() # Compute the flux #s = str(e_flu_min) + ' ' + str(e_flu_max) + ' err' #xspec.AllModels.calcFlux(s) #xspec_save_file.write(s + '\n') #flux_cutoffpl = xspec.AllData(1).flux #flux_grbm = [-1, -1, -1] flux_grbm = [10**grbm_lg10Flux, 10**grbm_lg10Flux_err[0] , 10**grbm_lg10Flux_err[1]] # print flux_grbm print 'fit_spectra: fitting cflux*(cutoffpl) model' xspec.AllModels.clear() # Load the cutoff power law model. model = xspec.Model('cflux*(cutoffpl)') xspec_save_file.write('model cflux*(cutoffpl)\n') xspec_save_file.write(str(e_flu_min)+'\n') # Emin model.cflux.Emin = [e_flu_min] xspec_save_file.write(str(e_flu_max) + '\n') # Emax model.cflux.Emax = [e_flu_max] xspec_save_file.write('-6\n') # lg10Flux model.cflux.lg10Flux = [-6] xspec_save_file.write('1\n') # PhoIndex model.cutoffpl.PhoIndex = [1] xspec_save_file.write('100\n') # HighECut model.cutoffpl.HighECut = [100] xspec_save_file.write('1,-1\n') # norm model.cutoffpl.norm = [1] model.cutoffpl.norm.frozen = True # Fit the model. try: fit_model_spectrum(fit_statistic, xspec_save_file) xspec.Fit.error("3 4 5") xspec_save_file.write('error 3 4 5\n') fit_model_spectrum(fit_statistic, xspec_save_file) cutoffpl_lg10Flux = model.cflux.lg10Flux.values[0] cutoffpl_alpha = model.cutoffpl.PhoIndex.values[0] cutoffpl_ec = model.cutoffpl.HighECut.values[0] cutoffpl_chi2 = xspec.Fit.statistic cutoff_pl_dof = xspec.Fit.dof goodness_of_fit_cutoffpl = xspec.Fit.statistic dof_cutoffpl = xspec.Fit.dof xspec.Fit.error("3") xspec_save_file.write('error 3\n') cutoffpl_lg10Flux_err = xspec.AllModels(1)(3).error xspec.Fit.error("4") xspec_save_file.write('error 4\n') err = xspec.AllModels(1)(4).error cutoffpl_alpha_err = ( err[1] - err[0] ) /2.0 xspec.Fit.error("5") xspec_save_file.write('error 5\n') err = xspec.AllModels(1)(5).error cutoffpl_ec_err = ( err[1] - err[0] ) / 2.0 #Plot the results plot_file1 = 'cutoffpl_'+plot_file xspec.Plot.device = plot_file1+'/ps' xspec.Plot('uf') except Exception, msg: print msg model.show() # Compute the flux #s = str(e_flu_min) + ' ' + str(e_flu_max) + ' err' #xspec.AllModels.calcFlux(s) #xspec_save_file.write(s + '\n') #flux_cutoffpl = xspec.AllData(1).flux #flux_cutoffpl = [-1, -1, -1] flux_cutoffpl = [10**cutoffpl_lg10Flux, 10**cutoffpl_lg10Flux_err[0] , 10**cutoffpl_lg10Flux_err[1]] # print flux_cutoffpl # Pick the best fit print '' reduced_fit_stat_grbm = goodness_of_fit_grbm / dof_grbm reduced_fit_stat_cpl = goodness_of_fit_cutoffpl / dof_cutoffpl if reduced_fit_stat_grbm < reduced_fit_stat_cpl: print '==> Band model gives best fit' best_model = 'grbm' print '' print 'The Band model spectral parameters are' print ' alpha = {0} +/- {1}'.format(grbm_alpha, grbm_alpha_err) print ' beta = {0} +/- {1}'.format(grbm_beta, grbm_beta_err) print ' E_b = {0} +/- {1} keV'.format(grbm_tem, grbm_tem_err) model_data = [best_model, grbm_alpha, grbm_alpha_err, grbm_beta, grbm_beta_err, grbm_tem, grbm_tem_err] print ' {0}/dof = {1} / {2}'.format(fit_statistic, goodness_of_fit_grbm, dof_grbm) flux = flux_grbm[0] flux_err = ( flux_grbm[2] - flux_grbm[1] ) / 2.0 print ' flux = {0} +/- {1} erg/cm^2/s'.format(flux, flux_err) fluence = flux * duration fluence_err = flux_err * duration print ' fluence = {0} +/- {1} erg/cm^2 ({2} s)'.format(fluence, fluence_err, duration) goodness_of_fit = goodness_of_fit_grbm dof = dof_grbm else: print '==> Cutoff Power Law gives best fit' best_model = 'cutoffpl' print '' print 'The cutoff power law spectral parameters are' print ' alpha = {0} +/- {1}'.format(cutoffpl_alpha, cutoffpl_alpha_err) print ' E_c = {0} +/- {1} keV'.format(cutoffpl_ec, cutoffpl_ec_err) model_data = [best_model, cutoffpl_alpha, cutoffpl_alpha_err, cutoffpl_ec, cutoffpl_ec_err] print ' {0}/dof = {1} / {2}'.format(fit_statistic, goodness_of_fit_cutoffpl, dof_cutoffpl) flux = flux_cutoffpl[0] flux_err = ( flux_cutoffpl[2] - flux_cutoffpl[1] ) / 2.0 print ' flux = {0} +/- {1} erg/cm^2/s'.format(flux, flux_err) fluence = flux * duration fluence_err = flux_err * duration print ' fluence = {0} +/- {1} erg/cm^2 ({2} s)'.format(fluence, fluence_err, duration) goodness_of_fit = goodness_of_fit_cutoffpl dof = cutoff_pl_dof print '' # Close the XSpec log file xspec.Xset.closeLog() xspec_save_file.close() print '' print 'XSpec commands are stored in {0}'.format(xspec_save_file_name) print '' return [model_data, [goodness_of_fit, dof], [fluence, fluence_err]] ##################################################################### #### def main(): """Run basic GRB processing on Fermi/GBM data. This script takes a Fermi/GBM burst trigger directory as input. It constructs light curves, identifies the burst, computes burst durations, and fits a spectrum to the burst. The analysis is fully automated and is intended to provide a quick look analysis of a Fermi/GBM trigger. Usage: do_gbm.py trigger_id timestring detectors timezero trigger_id -- The Fermi/GBM trigger id of this burst. timestring -- NONE tells the script to identify the source and background intervals automatically. b1-b2:s1-s2:b3-b4 where b1, b2, s1, s2, b3, b4 are the start and stop times of the first background interval (b1,b2), the source interval (s1,s2), and the second background interval (b3,b4) respectively, in MET seconds. detectors -- ALL tells the script to use all the detectors, TRIGGER tells the script to use only the detectors that were triggered by the burst, nnnnnnnnnnnn is a string of 0s and 1s where 1 indicates that this detector is to be used. The string starts with the n0 detector and ends with the nb detector. timezero -- TRIGTIME tells the script to set the trigger time to the value of the TRIGTIME keyword in the input detector file. x.y tells the script to use the time x.y (MET seconds) as the trigger time. """ version = '0.1.1 beta' date = '2011-12-16' print print '==> {0}'.format(__file__) print ' version {0} date {1}'.format(version, date) print '' # Energy range for T_90 calculation. e_t90_min = 50 # keV e_t90_max = 300 # keV # Energy range for fluence calculation. e_flu_min = 10 # keV e_flu_max = 1000 # keV # Get command line arguments if len(sys.argv) != 5: print 'usage: {0} trigger_name bkg1_start-bkg1_stop:src_start-src_stop:bkg2_start-bkg2_stop detectors timezero'.format(__file__) sys.exit() trigger_name = sys.argv[1] time_string = sys.argv[2] detector_string = sys.argv[3] timezero = sys.argv[4] times = get_times(time_string) if times == None: user_times = False elif times == 'ERROR': print 'error, bad user-supplied time intervals' sys.exit() else: user_times = True print 'using user-supplied time intervals' intervals = times # For now assume that the user-supplied times are in MET. files = get_detectors(detector_string, trigger_name) if len(files) < 1: print 'error, no GBM detector files found for trigger {0}'.format(trigger_name) sys.exit() if timezero == 'TRIGTIME': print 'timezero = TRIGTIME from detector files' else: print 'timezero = {0} (s; MET)'.format(timezero) # Set up some working lists sn_ratio = [] t100 = [] bin_width = [] print '' print '*****' print 'Construct light curves for {0}'.format(trigger_name) print '' # Run gtibin on each detector data file to create light curves and # PHA2 files. file_number = -1 for detector_file in files: print '' print 'working on {0}'.format(detector_file) file_number += 1 # Get time information from this detector file fitsfile = pyfits.open(detector_file, 'readonly') tstart = fitsfile['EVENTS'].header['TSTART'] tstop = fitsfile['EVENTS'].header['TSTOP'] trigtime = fitsfile['EVENTS'].header['TRIGTIME'] fitsfile.close() print ' TSTART = {0} s (MET)'.format(tstart) print ' TSTOP = {0} s (MET)'.format(tstop) print 'TRIGTIME = {0} s (MET)'.format(trigtime) if timezero != 'TRIGTIME': trigtime = float(timezero) print 'setting TRIGTIME = {0} s (MET)'.format(timezero) file_root = detector_file[0:detector_file.rfind('.fit')] i = detector_file.find('_') + 1 j = detector_file.find('_',i) data_type = file_root[i:j] # Time bin width # Valerie's tutorial says # tte: 0.128 s default # cspec: 1.024 s default after TRIGTIME # ctime: 0.064 s default if data_type == 'tte': default_dtime = 0.128 # seconds elif data_type == 'cspec': default_dtime = 1.024 # seconds elif data_type == 'ctime': default_dtime = 0.064 # seconds else: nbins = 1024 default_dtime = ( tstop - tstart ) / nbins # Arbitrary bin size print 'warning, unknown data type ({0})'.format(data_type) print 'warning, arbitrarily binning data into {1} bins'.format(nbins) dtime = default_dtime print 'light curve bin width: {0} s'.format(dtime) # Create the light curve file lightcurve_file = detector_file.replace('.fit', '.lc') status = run_gtbin_lightcurve(detector_file, lightcurve_file, tstart, tstop, dtime) if status is not None: print 'warning, failed to generate light curve for {0}'.format(detector_file) continue # Skip this detector # Determine the background and source intervals. if user_times: duration = intervals[3] - intervals[2] print 'using user-defined time intervals' else: # Determine the intervals intervals = time_intervals(lightcurve_file, trigtime) if intervals is None: print 'warning, unable to determine source and background intervals for {0}'.format(lightcurve_file) sn_ratio.append(-1) t100.append(-1) bin_width.append(-1) continue # Skip this detector else: src_start = intervals[2] src_stop = intervals[3] duration = src_stop - src_start print 'preliminary burst duration = {0} s'.format(duration) # Compute the signal to noise of the source part of the spectrum. list = compute_signal_to_noise(lightcurve_file, intervals, user_times, trigtime) mean_sn = list[0] intervals = list[1] print 'mean S/N per bin = {0}'.format(mean_sn) sn_ratio.append( mean_sn ) src_start = intervals[2] src_stop = intervals[3] duration = src_stop - src_start print 'refined burst duration = {0} s'.format(duration) t100.append( duration ) bin_width.append( dtime ) print '' print '*****' print 'Find the best light curves for {0}'.format(trigger_name) print '' for i in range( len(files) ): print '{0}: {1} S/N = {2} t = {3} s'.format(i, files[i], sn_ratio[i], t100[i]) print '' # Determine which light curve has the best signal-to-noise # ratio for the source. threshold_sn = 1.5 print 'S/N per bin threshold: {0}'.format(threshold_sn) n_good = 0 i_best = 0 max_sn = sn_ratio[0] for i in range( len(sn_ratio) ): if sn_ratio[i] > threshold_sn: n_good += 1 if sn_ratio[i] > max_sn: max_sn = sn_ratio[i] i_best = i if n_good > 0: best_event_file = files[i_best] print 'Best source S/N: {0} S/N={1}'.format(files[i_best], max_sn) else: print 'error, no light curves with S/N > {0}'.format(threshold_sn) sys.exit() # Create a new events file and corresponding light curve using the # events in the T90 energy range (keV) from the best events # file. suffix = '_' + str(e_t90_min) + '_' + str(e_t90_max) + '.evt' event_file = best_event_file.replace('.fit', suffix) select_on_energy(best_event_file, e_t90_min, e_t90_max, event_file) # Create the light curve file from events in the specified energy # range. lightcurve_file = event_file.replace('.evt', '.lc') status = run_gtbin_lightcurve(event_file, lightcurve_file, tstart, tstop, dtime) if status is not None: print 'warning, failed to generate light curve for {0}'.format(detector_file) sys.exit() # Determine the background and source intervals. if user_times: duration = intervals[3] - intervals[2] print 'using user-defined time intervals' else: # Determine the intervals intervals = time_intervals(lightcurve_file, trigtime) if intervals is None: print 'warning, unable to determine source and background intervals for {0}'.format(lightcurve_file) sys.exit() else: src_start = intervals[2] src_stop = intervals[3] duration = src_stop - src_start print 'burst duration = {0} s'.format(duration) # Set the background and source intervals using the detector with # the best light curve. fits_bin_file = event_file.replace('.evt', '_tbins.fits') b1_start = intervals[0] b1_stop = intervals[1] s_start = intervals[2] s_stop = intervals[3] create_time_bins_file(b1_start, b1_stop, s_start, s_stop, fits_bin_file) duration = s_stop - s_start print 'burst duration ({0}-{1} keV): {2} s'.format(e_t90_min, e_t90_max, duration) # Compute the T90 duration.`< t = compute_duration(lightcurve_file, s_start, s_stop, b1_start, b1_stop, 90) t90 = t[0] t90_err = t[1] t90_start = t[2] - trigtime t90_stop = t[3] - trigtime print 'burst T_90: {0} +/- {1} s from {2} to {3}'.format(t90, t90_err, t90_start, t90_stop) print '' print 'Time Intervals' print ' 1st background start: {0}'.format(intervals[0]) print ' 1st background stop: {0}'.format(intervals[1]) print ' source start: {0}'.format(intervals[2]) print ' source stop: {0}'.format(intervals[3]) print ' 2nd background start: {0}'.format(intervals[4]) print ' 2nd background stop: {0}'.format(intervals[5]) print '' print '' print '*****' print 'Extract the spectra for {0}'.format(trigger_name) print '' # Extract spectra for each light curve where a source was found. for detector_file in files: print '' print 'working on {0}'.format(detector_file) file_root = detector_file[0:detector_file.rfind('.fit')] # Create the PHA2 files pha_file = file_root + '.pha' command = ['gtbin', 'evfile='+detector_file, 'scfile=NONE', 'outfile='+pha_file, 'algorithm=PHA2', 'tbinalg=FILE', 'tbinfile='+fits_bin_file, 'chatter=2', 'clobber=YES'] try: subprocess.check_call(command) except: print 'warning, subprocess call failed for {0}'.format(command) print '' continue print 'created spectrum {0}'.format(pha_file) print '' print '' print '*****' print 'Find the best spectra for {0}'.format(trigger_name) print '' # Determine which light curves have a mean S/N per bin # for the source that is above some threshold # Create the PHA2 files spectrum_list = [] for i in range( len(sn_ratio) ): if ( sn_ratio[i] > threshold_sn ): pha_file = files[i].replace('.fit', '.pha') spectrum_list.append(pha_file) print '{0} has S/N = {1}'.format(pha_file, sn_ratio[i]) #if ( '_b0_' in files[i] ) or ( '_b1_' in files[i] ): # file_root = files[i][0:files[i].rfind('.fit')] # pha_file = file_root + '.pha' # spectrum_list.append(pha_file) # print '{0} has S/N = {1}'.format(pha_file, sn_ratio[i]) # Add the response matrix data to the PHA file. #update_pha_keywords() Is this the right approach # For now just create a list of responses manually response_list = [] for spectrum in spectrum_list: response = find_response_file(spectrum) response_list.append(response) print '' print '# Spectrum Response' for i in range( len(spectrum_list) ): print '{0}: {1} {2}'.format(i, spectrum_list[i], response_list[i]) # The right solution is probably to put the # response matrix and background file information # into the FITS header. # Do the XSpec analysis. print '' print '*****' print 'Perform the joint spectral fit for {0}'.format(trigger_name) print '' #fit_statistic = 'cstat' # Valerie's tutorial says to use cstat. fit_statistic = 'pgstat' # try pgstat for non-poisson errors. print 'The fitting statistic is {0}'.format(fit_statistic) spectra_plot_file = trigger_name+'_xspec.ps' spectrum_results = fit_spectra(trigger_name, spectrum_list, response_list, fit_statistic, e_flu_min, e_flu_max, duration, spectra_plot_file) model_data = spectrum_results[0] model_type = model_data[0] if model_type == 'grbm': model = 'Band Function' alpha = model_data[1] alpha_err = model_data[2] beta = model_data[3] beta_err = model_data[4] epeak = model_data[5] * (2.0 + alpha) epeak_err = math.sqrt( (2.0 + alpha)**2 * model_data[6]**2 + model_data[5]**2 * alpha_err**2 ) elif model_type == 'cutoffpl': model = 'Power Law with Exponential Cutoff' alpha = model_data[1] alpha_err = model_data[2] epeak = model_data[3] * (2.0 + alpha) epeak_err = math.sqrt( (2.0 + alpha)**2 * model_data[4]**2 + model_data[3]**2 * alpha_err**2 ) gof_data = spectrum_results[1] goodness_of_fit = gof_data[0] dof = gof_data[1] fluence_data = spectrum_results[2] fluence = fluence_data[0] fluence_err = fluence_data[1] # Tell the user what the results are print '' print '*****' print '*****' print 'Summary Information for {0}'.format(trigger_name) print '' fitsfile = pyfits.open(best_event_file, 'readonly') trigtime = fitsfile['EVENTS'].header['TRIGTIME'] object = fitsfile['EVENTS'].header['OBJECT'] fitsfile.close() print ' Source: {0}'.format(object) print ' Trigger Time: {0} (MET)'.format(trigtime) if timezero != 'TRIGTIME': print ' Used Timezero: {0} (MET)'.format(timezero) print ' Best Data File: {0}'.format(best_event_file) print 'Signal-to-Noise Per Bin: {0:0.2f}'.format(max_sn) print ' T_90: {0:0.3f} +/- {1:0.3f} s ({2}-{3} keV)'.format(t90, t90_err, e_t90_min, e_t90_max) print ' Spectrum: {0}'.format(model) if model_type == 'grbm': print ' alpha: {0:5.2f} +/- {1:5.2f}'.format(alpha, alpha_err) print ' beta: {0:5.2f} +/- {1:5.2f}'.format(beta, beta_err) print ' E_peak: {0:0.1f} +/- {1:0.1f} keV'.format(epeak, epeak_err) elif model_type == 'cutoffpl': print ' alpha: {0:5.2f} +/- {1:5.2f}'.format(alpha, alpha_err) print ' E_peak: {0:0.1f} +/- {1:0.1f} keV'.format(epeak, epeak_err) print ' fluence ({0}-{1} keV): {2:8.2e} +/- {3:8.2e} erg/cm^-2'.format(e_flu_min, e_flu_max, fluence, fluence_err) print '' print 'goodness of fit / dof ({0}) : {1} / {2}'.format(fit_statistic, goodness_of_fit, dof) print '' #### if __name__ == '__main__': main()