From dc08f63ca2bc35cb1dbb4020971ac86f85e3a1bf Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Fri, 10 Jun 2016 17:54:15 +0800 Subject: Add ciao_calc_bkg.py: calculate the background surface brightness --- ciao_calc_bkg.py | 152 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 152 insertions(+) create mode 100755 ciao_calc_bkg.py (limited to 'ciao_calc_bkg.py') diff --git a/ciao_calc_bkg.py b/ciao_calc_bkg.py new file mode 100755 index 0000000..cf5fdb5 --- /dev/null +++ b/ciao_calc_bkg.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python3 +# +# Calculate the *background surface brightness (SB)* level from the +# *corrected background spectrum*. +# The calculated background SB value is used to provide constraint for +# surface brightness profile (SBP) fitting, and is also adopted to +# subtract the background contribution before carrying out the SBP +# deprojection. +# +# Weitian LI +# Created: 2016-06-10 +# Updated: 2016-06-10 +# + +import argparse +import subprocess +import tempfile +import json +from collections import OrderedDict + +import numpy as np +from astropy.io import fits + + +def energy2channel(energy): + """ + Convert energy (eV) to Chandra ACIS channel number. + + Reference: http://cxc.harvard.edu/ciao/dictionary/pi.html + """ + return int((energy / 14.6) + 1) + + +def parse_erange(erange): + """ + Parse the given erange string, and return the lower and upper energies. + """ + e_low, e_high = map(float, erange.split(":")) + return (e_low, e_high) + + +def calc_exp(expmap, regfile): + """ + Calculate the area of the background spectrum extraction region, + and the mean exposure (used as the MEAN_SRC_EXPOSURE) of that region. + """ + tf = tempfile.NamedTemporaryFile() + cmd_args = [ + "dmextract", + "infile=%s[bin sky=region(%s)]" % (expmap, regfile), + "outfile=%s" % tf.name, + "opt=generic", "clobber=yes" + ] + subprocess.run(["punlearn", "dmextract"]) + subprocess.run(args=cmd_args) + with fits.open(tf.name) as hist_fits: + total_exp = hist_fits["HISTOGRAM"].data["COUNTS"][0] + total_exp_err = hist_fits["HISTOGRAM"].data["ERR_COUNTS"][0] + area = hist_fits["HISTOGRAM"].data["AREA"][0] + mean_exp = hist_fits["HISTOGRAM"].data["SUR_BRI"][0] + mean_exp_err = hist_fits["HISTOGRAM"].data["SUR_BRI_ERR"][0] + tf.close() + return { + "total_exp": total_exp, + "total_exp_err": total_exp_err, + "area": area, + "mean_exp": mean_exp, + "mean_exp_err": mean_exp_err, + } + + +def calc_spec_counts(spec, erange): + """ + Calculate the spectrum total counts within the specified energy range. + Also extract the EXPOSURE and BACKSCAL information. + """ + specfits = fits.open(spec) + channel = specfits["SPECTRUM"].data["CHANNEL"] + counts = specfits["SPECTRUM"].data["COUNTS"] + channel_low, channel_high = map(energy2channel, erange) + chan_idx = np.where(np.logical_and(channel >= channel_low, + channel <= channel_high)) + total_counts = np.sum(counts[chan_idx]) + return { + "energy_low": erange[0], + "energy_high": erange[1], + "channel_low": channel_low, + "channel_high": channel_high, + "counts": int(total_counts), + "exposure": specfits["SPECTRUM"].header["EXPOSURE"], + "backscal": specfits["SPECTRUM"].header["BACKSCAL"], + } + + +def main(): + parser = argparse.ArgumentParser( + description="Calculate the background surface brightness") + parser.add_argument("-b", "--orig-bkg", dest="orig_bkg", required=True, + help="original/uncorrected background spectrum " + + "(e.g., blanksky_lbkg.pi), from which to get the " + + "original EXPOSURE and BACKSCAL, etc.") + parser.add_argument("-B", "--corr-bkg", dest="corr_bkg", required=True, + help="corrected background spectrum for the " + + "Galactic and cosmic background radiations " + + "(e.g., bkgcorr_blanksky_lbkg.pi)") + parser.add_argument("-r", "--bkg-region", dest="bkg_region", required=True, + help="region used to extract the background " + + "spectrum (e.g., lbkg.reg") + parser.add_argument("-e", "--expmap", dest="expmap", required=True, + help="exposure map w.r.t the spectrum") + parser.add_argument("-E", "--erange", dest="erange", required=True, + help="energy range used for the exposure map " + + "(e.g., 700:7000)") + parser.add_argument("-o", "--outfile", dest="outfile", + required=False, default="sb_bkg.json", + help="json file to save the background SB results") + args = parser.parse_args() + + e_low, e_high = parse_erange(args.erange) + orig_bkg_results = calc_spec_counts(args.orig_bkg, erange=(e_low, e_high)) + corr_bkg_results = calc_spec_counts(args.corr_bkg, erange=(e_low, e_high)) + exp_results = calc_exp(args.expmap, regfile=args.bkg_region) + + counts = corr_bkg_results["counts"] + exposure = corr_bkg_results["exposure"] + backscal = corr_bkg_results["backscal"] + orig_backscal = orig_bkg_results["backscal"] + mean_exp = exp_results["mean_exp"] + bkg_sb = counts / exposure / mean_exp / backscal * orig_backscal + + results = OrderedDict([ + ("energy_low", e_low), + ("energy_high", e_high), + ("channel_low", corr_bkg_results["channel_low"]), + ("channel_high", corr_bkg_results["channel_high"]), + ("counts", counts), + ("exposure", exposure), + ("backscal", backscal), + ("total_exp", exp_results["total_exp"]), + ("total_exp_err", exp_results["total_exp_err"]), + ("area", exp_results["area"]), + ("mean_exp", exp_results["mean_exp"]), + ("mean_exp_err", exp_results["mean_exp_err"]), + ("bkg_sb", bkg_sb), + ]) + results_json = json.dumps(results, indent=2) + print(results_json) + open(args.outfile, "w").write(results_json+"\n") + + +if __name__ == "__main__": + main() -- cgit v1.2.2