#!/usr/bin/env python3 # -*- coding: utf-8 -*- # # Calculate the surface brightness concentration (i.e., C_{SB}), which # is an index/indicator of the cool core, and may be defined as: # (1) brightness(<=40kpc) / brightness(<=400kpc) # (2) brightness(<=0.048R500) / brightness(<=0.45R500) # # References: # [1] Santos, J. S., Rosati, P., Tozi, P., et al. 2008, A&A, 483, 35 # # Aaron LI # Created: 2016-04-28 # Updated: 2016-05-16 # # Change log: # 2016-05-16: # * Add some background subtraction support # * Use `subprocess.run` instead of `subprocess.call` # * PEP8 fixes # 2016-04-29: # * Fix order of executing ds9 and print message # * Fix C_SB calculation # * Add reference # * Fix "cuspiness" to "concentration" # * Add "name" and "obsid" to results # 2016-04-28: # * Add "csb_type" to results # # TODO/XXX: # whether the background should be subtracted for C_SB calculation?? # import sys import os import glob import json import argparse import subprocess from collections import OrderedDict from astropy.io import fits from make_r500_regions import get_r500, get_center def make_csb_region(regfile, center, r1, r2): """ Make the regions for C_SB and save. """ regions = [ "pie(%.2f,%.2f,0,%.2f,0,360)" % (center[0], center[1], r1), "pie(%.2f,%.2f,0,%.2f,0,360)" % (center[0], center[1], r2), ] open(regfile, "w").write("\n".join(regions) + "\n") def calc_csb(evt, expmap, regfile, r1, r2, bkg=None): """ Calculate the C_SB If `bkg` is provided, then background subtraction is considered for `dmextract` to calculate the *net* counts and surface brightness. """ csbfile = os.path.splitext(regfile)[0] + ".fits" cmd_args = [ "dmextract", "infile=%s[bin sky=@%s]" % (evt, regfile), "outfile=%s" % csbfile, "exp=%s" % expmap, "opt=generic", "clobber=yes" ] if bkg is not None: # consider background subtraction subprocess.run(["punlearn", "dmkeypar"]) ret = subprocess.run(args=["dmkeypar", evt, "EXPOSURE", "echo=yes"], check=True, stdout=subprocess.PIPE) exposure_evt = float(ret.stdout.decode("utf-8")) ret = subprocess.run(args=["dmkeypar", bkg, "EXPOSURE", "echo=yes"], check=True, stdout=subprocess.PIPE) exposure_bkg = float(ret.stdout.decode("utf-8")) bkg_norm = exposure_evt / exposure_bkg cmd_args += [ "bkg=%s[energy=700:7000][bin sky=@%s]" % (bkg, regfile), "bkgnorm=%s" % bkg_norm, "bkgexp=)exp" ] subprocess.run(["punlearn", "dmextract"]) subprocess.run(args=cmd_args) # read calculate C_SB data from output FITS with fits.open(csbfile) as csb_fits: csb_s_val = csb_fits["HISTOGRAM"].data["SUR_BRI"] csb_s_err = csb_fits["HISTOGRAM"].data["SUR_BRI_ERR"] # if bkg is not None: # bkg_csb_s_val = csb_fits["HISTOGRAM"].data["BG_SUR_BRI"] # bkg_csb_s_err = csb_fits["HISTOGRAM"].data["BG_SUR_BRI_ERR"] # calculate C_SB and error area_ratio = (r2 / r1) ** 2 csb = csb_s_val[0] / csb_s_val[1] / area_ratio csb_err = csb * ((csb_s_err[0] / csb_s_val[0])**2 + (csb_s_err[1] / csb_s_val[1])**2) ** 0.5 results = OrderedDict([ ("csb_s1", csb_s_val[0]), ("csb_s1_err", csb_s_err[0]), ("csb_s2", csb_s_val[1]), ("csb_s2_err", csb_s_err[1]), ("csb", csb), ("csb_err", csb_err), ]) return results def main(): parser = argparse.ArgumentParser( description="Calculate the surface brightness concentration") # exclusive argument group for C_SB definition grp_csb = parser.add_mutually_exclusive_group(required=True) grp_csb.add_argument("-K", "--kpc", dest="kpc", action="store_true", help="C_SB = brightness(<=0.048R500) " + "/ brightness(<=0.45R500)") grp_csb.add_argument("-R", "--r500", dest="r500", action="store_true", help="C_SB = brightness(<=40kpc) " + "/ brightness(<=400kpc)") # parser.add_argument("-A", "--no-ask", dest="no_ask", required=False, action="store_true", help="do NOT check region and ask") parser.add_argument("-j", "--json", dest="json", required=False, help="the *_INFO.json file " + "(default: find ../*_INFO.json)") parser.add_argument("-r", "--region", dest="region", required=False, default="sbprofile.reg", help="region from which to extract the center " + "coordinate (default: sbprofile.reg)") parser.add_argument("-i", "--infile", dest="infile", required=True, help="input energy-filtered EVT2 used to " + "calculate the C_SB") parser.add_argument("-e", "--expmap", dest="expmap", required=True, help="exposure map of the input image") parser.add_argument("-b", "--bkg", dest="bkg", default=None, help="background corresponding to the input evt2") parser.add_argument("-o", "--outfile", dest="outfile", required=True, help="output json file to store the C_SB results") # args = parser.parse_args() # default "*_INFO.json" info_json = glob.glob("../*_INFO.json")[0] if args.json: info_json = args.json json_str = open(info_json).read().rstrip().rstrip(",") info = json.loads(json_str) name = info["Source Name"] obsid = int(info["Obs. ID"]) r500 = get_r500(info) r500_kpc = r500["r500_kpc"] r500_pix = r500["r500_pix"] kpc_per_pix = r500["kpc_per_pix"] print("R500: %.2f (kpc), %.2f (pixel)" % (r500_kpc, r500_pix)) # get center coordinate xc, yc = get_center(args.region) if args.r500: csb_type = "r500" r1 = 0.048 * r500_pix r2 = 0.450 * r500_pix elif args.kpc: csb_type = "kpc" r1 = 40.0 / kpc_per_pix r2 = 400.0 / kpc_per_pix else: raise ValueError("Unknown C_SB definition") # make regions for C_SB regfile = os.path.splitext(args.outfile)[0] + ".reg" make_csb_region(regfile, center=(xc, yc), r1=r1, r2=r2) # check region with DS9 if not args.no_ask: print("Check the C_SB regions; overwrite the region file " + "'%s' if modified" % regfile, flush=True, file=sys.stderr) cmd = "ds9 %s -cmap he " % args.infile + \ "-regions format ciao -regions %s" % regfile subprocess.call(cmd, shell=True) ans = input("C_SB regions exceed CCD (No/yes/modified)? ") if ans == "" or ans[0] in "nN": csb_region_note = "OK" elif ans[0] in "yY": csb_region_note = "EXCESS" elif ans[0] in "mM": csb_region_note = "MODIFIED" else: csb_region_note = "???" else: csb_region_note = None # calculate the C_SB csb = calc_csb(args.infile, expmap=args.expmap, regfile=regfile, r1=r1, r2=r2, bkg=args.bkg) csb_data = OrderedDict([ ("name", name), ("obsid", obsid), ("csb_type", csb_type), ("csb_r1", r1), ("csb_r2", r2), ]) csb_data.update(csb) csb_data["csb_region"] = csb_region_note csb_data_json = json.dumps(csb_data, indent=2) print(csb_data_json) open(args.outfile, "w").write(csb_data_json + "\n") if __name__ == "__main__": main()