aboutsummaryrefslogtreecommitdiffstats
path: root/acispy
diff options
context:
space:
mode:
Diffstat (limited to 'acispy')
-rw-r--r--acispy/__init__.py0
-rw-r--r--acispy/acis.py91
-rw-r--r--acispy/ds9.py30
-rw-r--r--acispy/pfiles.py38
-rw-r--r--acispy/region.py251
-rw-r--r--acispy/spectrum.py44
6 files changed, 454 insertions, 0 deletions
diff --git a/acispy/__init__.py b/acispy/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/acispy/__init__.py
diff --git a/acispy/acis.py b/acispy/acis.py
new file mode 100644
index 0000000..91d9fcf
--- /dev/null
+++ b/acispy/acis.py
@@ -0,0 +1,91 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Chandra ACIS utilities
+"""
+
+import math
+import subprocess
+import re
+
+
+class ACIS:
+ """
+ Chandra ACIS detector properties and utilities.
+
+ References
+ ----------
+ [1] CIAO - dictionary - ACIS (advanced camera for imaging and spectroscopy)
+ http://cxc.harvard.edu/ciao/dictionary/acis.html
+ [2] CIAO - Dictionary - PI (pulse invariant)
+ http://cxc.harvard.edu/ciao/dictionary/pi.html
+ """
+ # Pixel size
+ pixel2arcsec = 0.492 # [arcsec]
+ # Number of channels
+ nchannel = 1024
+ # Channel energy width
+ echannel = 14.6 # [eV]
+
+ @classmethod
+ def energy2channel(self, energy):
+ """
+ Convert energy [eV] to channel number.
+ """
+ return math.floor(energy/self.echannel + 1)
+
+ @classmethod
+ def get_type(self, filepath):
+ """
+ Determine the Chandra ACIS type (``I`` or ``S``) according the
+ active ACIS chips.
+
+ Parameters
+ ----------
+ filepath : str
+ Path to the input FITS file
+
+ Returns
+ -------
+ acis_type : str
+ ``I`` if ACIS-I, ``S`` if ACIS-S;
+ otherwise, ``ValueError`` raised.
+ """
+ subprocess.check_call(["punlearn", "dmkeypar"])
+ detnam = subprocess.check_output([
+ "dmkeypar", "infile=%s" % filepath, "keyword=DETNAM", "echo=yes"
+ ]).decode("utf-8").strip()
+ if re.match(r"^ACIS-0123", detnam):
+ return "I"
+ elif re.match(r"^ACIS-[0-6]*7", detnam):
+ return "S"
+ else:
+ raise ValueError("unknown chip combination: %s" % detnam)
+
+ @classmethod
+ def get_chips_str(self, filepath, sep=":"):
+ """
+ Return a string of the chips of interest according to the
+ active ACIS type.
+
+ Parameters
+ ----------
+ filepath : str
+ Path to the input FITS file
+ sep : str, optional
+ Separator to join the chip ranges, e.g., 0:3, 0-3
+
+ Returns
+ -------
+ chips : str
+ ``0:3`` if ACIS-I, ``7`` if ACIS-S;
+ otherwise, ``ValueError`` raised.
+ """
+ acis_type = self.get_type(filepath)
+ if acis_type == "I":
+ return sep.join(["0", "3"])
+ elif acis_type == "S":
+ return "7"
+ else:
+ raise ValueError("unknown ACIS type")
diff --git a/acispy/ds9.py b/acispy/ds9.py
new file mode 100644
index 0000000..e8c8e88
--- /dev/null
+++ b/acispy/ds9.py
@@ -0,0 +1,30 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Wrapper function to view FITS files using DS9.
+"""
+
+import subprocess
+
+
+def ds9_view(filename, regfile=None, regformat="ciao", regsystem="physical",
+ cmap="he", binfactor=2, scale="linear", smooth=None):
+ """
+ Wrapper function to view FITS files using DS9.
+ """
+ cmd = [
+ "ds9", filename,
+ "-cmap", cmap,
+ "-bin", "factor", str(binfactor),
+ "-scale", scale,
+ ]
+ if regfile:
+ cmd += [
+ "-regions", "format", regformat,
+ "-regions", "system", regsystem,
+ "-regions", regfile,
+ ]
+ if smooth:
+ cmd += ["-smooth", "yes", "-smooth", "radius", str(smooth)]
+ subprocess.check_call(cmd)
diff --git a/acispy/pfiles.py b/acispy/pfiles.py
new file mode 100644
index 0000000..8a80a14
--- /dev/null
+++ b/acispy/pfiles.py
@@ -0,0 +1,38 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+#
+# Weitian LI
+# 2017-02-06
+
+"""
+Prepare the CIAO parameter files and setup the PFILES environment
+variable to keep the pfiles locally, in order to avoid the conflicts
+between multiple instance of the same CIAO tools.
+"""
+
+import os
+import subprocess
+import shutil
+
+
+def setup_pfiles(tools):
+ """
+ Copy the parameter files of the specified tools to the current
+ working directory, and setup the ``PFILES`` environment variable.
+
+ Parameters
+ ----------
+ tools : list[str]
+ Name list of the tools to be set up
+ """
+ for tool in tools:
+ pfile = subprocess.check_output([
+ "paccess", tool
+ ]).decode("utf-8").strip()
+ subprocess.check_call(["punlearn", tool])
+ try:
+ shutil.copy(pfile, ".")
+ except shutil.SameFileError:
+ pass
+ # Setup the ``PFILES`` environment variable
+ os.environ["PFILES"] = "./:" + os.environ["PFILES"]
diff --git a/acispy/region.py b/acispy/region.py
new file mode 100644
index 0000000..09e42fa
--- /dev/null
+++ b/acispy/region.py
@@ -0,0 +1,251 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Region utilities.
+
+NOTE:
+Only support the commonly used region shapes in CIAO format.
+While more complex shapes and DS9 format are not considered
+at the moment.
+"""
+
+import os
+import re
+
+
+class Point:
+ """
+ Point region: point(xc,yc)
+ """
+ shape = "point"
+
+ def __init__(self, xc=None, yc=None):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f)" % (self.shape, self.xc, self.yc)
+
+
+class Circle:
+ """
+ Circle region: circle(xc,yc,radius)
+ """
+ shape = "circle"
+
+ def __init__(self, xc=None, yc=None, radius=None):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+ self.radius = radius
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+ self.radius = float(g[3])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f,%.4f)" % (self.shape, self.xc, self.yc,
+ self.radius)
+
+
+class Annulus:
+ """
+ Annulus region:
+ annulus(xc,yc,radius,radius2)
+ + xc, yc : center
+ + radius, radius2 : inner, outer radius
+ """
+ shape = "annulus"
+
+ def __init__(self, xc=None, yc=None, radius=None, radius2=None):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+ self.radius = radius
+ self.radius2 = radius2
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+ self.radius = float(g[3])
+ self.radius2 = float(g[4])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f,%.4f,%.4f)" % (self.shape, self.xc, self.yc,
+ self.radius, self.radius2)
+
+
+class Pie:
+ """
+ Pie region:
+ pie(xc,yc,radius,radius2,angle1,angle2)
+ + xc, yc : center
+ + radius, radius2 : inner, outer radius
+ + angle1, angle2: start, end angle [0, 360)
+ """
+ shape = "pie"
+
+ def __init__(self, xc=None, yc=None, radius=None, radius2=None,
+ angle1=0.0, angle2=360.0):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+ self.radius = radius
+ self.radius2 = radius2
+ self.angle1 = angle1
+ self.angle2 = angle2
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+ self.radius = float(g[3])
+ self.radius2 = float(g[4])
+ self.angle1 = float(g[5])
+ self.angle2 = float(g[6])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f,%.4f,%.4f,%.4f,%.4f)" % (
+ self.shape, self.xc, self.yc, self.radius, self.radius2,
+ self.angle1, self.angle2)
+
+
+class Ellipse:
+ """
+ Ellipse region:
+ ellipse(xc,yc,radius,radius2,rotation)
+ + xc, yc : center
+ + radius, radius2 : semi-major / semi-minor axis
+ + rotation: rotation angle [0, 360)
+ """
+ shape = "ellipse"
+
+ def __init__(self, xc=None, yc=None, radius=None, radius2=None,
+ rotation=0.0):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+ self.radius = radius
+ self.radius2 = radius2
+ self.rotation = rotation
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+ self.radius = float(g[3])
+ self.radius2 = float(g[4])
+ self.rotation = float(g[5])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f,%.4f,%.4f,%.4f)" % (
+ self.shape, self.xc, self.yc, self.radius,
+ self.radius2, self.rotation)
+
+
+class Box:
+ """
+ Box region: box(xc,yc,width,height,rotation)
+ """
+ shape = "box"
+
+ def __init__(self, xc=None, yc=None, width=None, height=None,
+ rotation=0.0):
+ if isinstance(xc, str):
+ self.parse(regstr=xc)
+ else:
+ self.xc = xc
+ self.yc = yc
+ self.width = width
+ self.height = height
+ self.rotation = rotation
+
+ def parse(self, regstr):
+ g = re.split(r"[(),]", regstr)
+ self.xc = float(g[1])
+ self.yc = float(g[2])
+ self.width = float(g[3])
+ self.height = float(g[4])
+ self.rotation = float(g[5])
+
+ def __str__(self):
+ return "%s(%.4f,%.4f,%.4f,%.4f,%.4f)" % (
+ self.shape, self.xc, self.yc, self.width,
+ self.height, self.rotation)
+
+
+class Regions:
+ """
+ Manipulate the region files (CIAO format), as well as parse regions
+ from strings.
+ """
+ REGION_SHAPES = {
+ "point": Point,
+ "circle": Circle,
+ "annulus": Annulus,
+ "pie": Pie,
+ "ellipse": Ellipse,
+ "box": Box,
+ }
+
+ def __init__(self, regfile=None):
+ if regfile:
+ self.regions = self.load(regfile)
+ else:
+ self.regions = []
+
+ def load(self, infile):
+ regstr = []
+ for line in open(infile):
+ if not (re.match(r"^\s*#.*$", line) or re.match(r"^\s*$", line)):
+ regstr.append(line.strip())
+ self.regions = self.parse(regstr)
+
+ def save(self, outfile, clobber=False):
+ if (not clobber) and os.path.exists(outfile):
+ raise OSError("output file already exists: %s" % outfile)
+ regstr = [str(reg) for reg in self.regions]
+ open(outfile, "w").write("\n".join(regstr) + "\n")
+
+ @classmethod
+ def parse(self, regstr):
+ """
+ Parse the given (list of) region string(s), and return the
+ corresponding parsed region objects.
+ """
+ if isinstance(regstr, list):
+ return [self.parse_single(reg) for reg in regstr]
+ else:
+ return self.parse_single(regstr)
+
+ @classmethod
+ def parse_single(self, regstr):
+ """
+ Parse the given single region string to its corresponding
+ region object.
+ """
+ shape = regstr.strip().split('(')[0].lower()
+ if shape in self.REGION_SHAPES:
+ return self.REGION_SHAPES[shape](regstr)
+ else:
+ raise ValueError("unknown region shape: %s" % shape)
diff --git a/acispy/spectrum.py b/acispy/spectrum.py
new file mode 100644
index 0000000..ad38b0d
--- /dev/null
+++ b/acispy/spectrum.py
@@ -0,0 +1,44 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Chandra ACIS spectrum.
+"""
+
+
+from astropy.io import fits
+
+from acis import ACIS
+
+
+class Spectrum:
+ """
+ Chandra ACIS spectrum
+ """
+ def __init__(self, filepath):
+ self.filepath = filepath
+ self.fitsobj = fits.open(filepath)
+ ext_spec = self.fitsobj["SPECTRUM"]
+ self.header = ext_spec.header
+ # spectral data
+ self.channel = ext_spec.data.columns["CHANNEL"].array
+ self.counts = ext_spec.data.columns["COUNTS"].array
+ # spectral keywords
+ self.EXPOSURE = self.header.get("EXPOSURE")
+ self.BACKSCAL = self.header.get("BACKSCAL")
+
+ def calc_pb_flux(self, elow=9500, ehigh=12000):
+ """
+ Calculate the particle background flux:
+ flux = counts / exposure / area
+
+ Parameters
+ ----------
+ elow, ehigh : float, optional
+ Lower and upper energy limit for the particle background.
+ """
+ chlow = ACIS.energy2channel(elow)
+ chhigh = ACIS.energy2channel(ehigh)
+ counts = self.counts[(chlow-1):chhigh].sum()
+ flux = counts / self.EXPOSURE / self.BACKSCAL
+ return flux
Copyright © 2017-2018 Aaron LI. All Rights Reserved.