diff options
| -rwxr-xr-x | astro/spectrum/crosstalk_deprojection.py | 498 |
1 files changed, 256 insertions, 242 deletions
diff --git a/astro/spectrum/crosstalk_deprojection.py b/astro/spectrum/crosstalk_deprojection.py index 617b8ff..0c04d77 100755 --- a/astro/spectrum/crosstalk_deprojection.py +++ b/astro/spectrum/crosstalk_deprojection.py @@ -207,7 +207,6 @@ class ARF: # {{{ https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html """ filename = None - fitsobj = None # only consider the "SPECTRUM" extension header = None energ_lo = None @@ -225,12 +224,12 @@ class ARF: # {{{ def __init__(self, filename): self.filename = filename - self.fitsobj = fits.open(filename) - ext_specresp = self.fitsobj["SPECRESP"] - self.header = ext_specresp.header - self.energ_lo = ext_specresp.data["ENERG_LO"] - self.energ_hi = ext_specresp.data["ENERG_HI"] - self.specresp = ext_specresp.data["SPECRESP"] + with fits.open(filename) as fitsobj: + ext_specresp = fitsobj["SPECRESP"] + self.header = ext_specresp.header + self.energ_lo = ext_specresp.data["ENERG_LO"] + self.energ_hi = ext_specresp.data["ENERG_HI"] + self.specresp = ext_specresp.data["SPECRESP"] def get_data(self, groupped=False, group_squeeze=False, copy=True): if groupped: @@ -277,13 +276,15 @@ class ARF: # {{{ """ if not hasattr(self, "f_interp") or self.f_interp is None: energy = self.get_energy() - arf = self.get_data(copy=False) + arf = self.get_data(copy=False) if verbose: print("INFO: interpolating ARF '%s' (may take a while) ..." % self.filename, file=sys.stderr) - f_interp = scipy.interpolate.interp1d(energy, arf, kind="cubic", - bounds_error=False, fill_value=0.0, assume_sorted=True) + f_interp = scipy.interpolate.interp1d( + energy, arf, kind="cubic", bounds_error=False, + fill_value=0.0, assume_sorted=True) self.f_interp = f_interp + if x is not None: return self.f_interp(x) else: @@ -304,7 +305,7 @@ class ARF: # {{{ return if verbose: print("INFO: Grouping spectrum '%s' ..." % self.filename, - file=sys.stderr) + file=sys.stderr) self.energy_channel = energy_channel self.grouping = grouping # interpolate the ARF w.r.t the spectral channel energies @@ -349,7 +350,6 @@ class RMF: # {{{ https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html """ filename = None - fitsobj = None ## extension "MATRIX" hdr_matrix = None energ_lo = None @@ -369,23 +369,23 @@ class RMF: # {{{ rmfimg = None def __init__(self, filename): - self.filename = filename - self.fitsobj = fits.open(filename) - ## "MATRIX" extension - ext_matrix = self.fitsobj["MATRIX"] - self.hdr_matrix = ext_matrix.header - self.energ_lo = ext_matrix.data["ENERG_LO"] - self.energ_hi = ext_matrix.data["ENERG_HI"] - self.n_grp = ext_matrix.data["N_GRP"] - self.f_chan = ext_matrix.data["F_CHAN"] - self.n_chan = ext_matrix.data["N_CHAN"] - self.matrix = ext_matrix.data["MATRIX"] - ## "EBOUNDS" extension - ext_ebounds = self.fitsobj["EBOUNDS"] - self.hdr_ebounds = ext_ebounds.header - self.channel = ext_ebounds.data["CHANNEL"] - self.e_min = ext_ebounds.data["E_MIN"] - self.e_max = ext_ebounds.data["E_MAX"] + self.filename = filename + with fits.open(filename) as fitsobj: + # "MATRIX" extension + ext_matrix = fitsobj["MATRIX"] + self.hdr_matrix = ext_matrix.header + self.energ_lo = ext_matrix.data["ENERG_LO"] + self.energ_hi = ext_matrix.data["ENERG_HI"] + self.n_grp = ext_matrix.data["N_GRP"] + self.f_chan = ext_matrix.data["F_CHAN"] + self.n_chan = ext_matrix.data["N_CHAN"] + self.matrix = ext_matrix.data["MATRIX"] + # "EBOUNDS" extension + ext_ebounds = fitsobj["EBOUNDS"] + self.hdr_ebounds = ext_ebounds.header + self.channel = ext_ebounds.data["CHANNEL"] + self.e_min = ext_ebounds.data["E_MIN"] + self.e_max = ext_ebounds.data["E_MAX"] def get_energy(self, mean="geometric"): """ @@ -428,7 +428,9 @@ class RMF: # {{{ rmfimg = np.zeros(shape=(n_energy, n_channel), dtype=rmf_dtype) for i in np.arange(n_energy)[self.n_grp > 0]: rmfimg[i, :] = _make_rmfimg_row(n_channel, rmf_dtype, - self.f_chan[i], self.n_chan[i], self.matrix[i]) + self.f_chan[i], + self.n_chan[i], + self.matrix[i]) self.rmfimg = rmfimg return self.rmfimg @@ -446,13 +448,11 @@ class Spectrum: # {{{ """ Class that deals with the X-ray spectrum file (usually *.pi). """ - filename = None - # FITS object return by `fits.open()' - fitsobj = None + filename = None # header of "SPECTRUM" extension - header = None + header = None # "SPECTRUM" extension data - channel = None + channel = None # name of the spectrum data column (i.e., type, "COUNTS" or "RATE") spec_type = None # unit of the spectrum data ("count" for "COUNTS", "count/s" for "RATE") @@ -460,50 +460,49 @@ class Spectrum: # {{{ # spectrum data spec_data = None # estimated spectral errors for each channel/group - spec_err = None + spec_err = None # statistical errors for each channel/group - stat_err = None + stat_err = None # grouping and quality - grouping = None - quality = None + grouping = None + quality = None # whether the spectral data being groupped - groupped = False + groupped = False # several important keywords - EXPOSURE = None - BACKSCAL = None - AREASCAL = None - RESPFILE = None - ANCRFILE = None - BACKFILE = None + EXPOSURE = None + BACKSCAL = None + AREASCAL = None + RESPFILE = None + ANCRFILE = None + BACKFILE = None # numpy dtype and FITS format code of the spectrum data - spec_dtype = None + spec_dtype = None spec_fits_format = None - # output filename for writing the spectrum if no filename provided - outfile = None - def __init__(self, filename, outfile=None): + def __init__(self, filename): self.filename = filename - self.fitsobj = fits.open(filename) - ext_spec = self.fitsobj["SPECTRUM"] - self.header = ext_spec.header.copy(strip=True) - colnames = ext_spec.columns.names - if "COUNTS" in colnames: - self.spec_type = "COUNTS" - elif "RATE" in colnames: - self.spec_type = "RATE" - else: - raise ValueError("Invalid spectrum file") - self.channel = ext_spec.data.columns["CHANNEL"].array - col_spec_data = ext_spec.data.columns[self.spec_type] - self.spec_data = col_spec_data.array.copy() - self.spec_unit = col_spec_data.unit - self.spec_dtype = col_spec_data.dtype - self.spec_fits_format = col_spec_data.format - # grouping and quality - if "GROUPING" in colnames: - self.grouping = ext_spec.data.columns["GROUPING"].array - if "QUALITY" in colnames: - self.quality = ext_spec.data.columns["QUALITY"].array + with fits.open(filename) as fitsobj: + ext_spec = fitsobj["SPECTRUM"] + self.header = ext_spec.header.copy(strip=True) + colnames = ext_spec.columns.names + if "COUNTS" in colnames: + self.spec_type = "COUNTS" + elif "RATE" in colnames: + self.spec_type = "RATE" + else: + raise ValueError("Invalid spectrum file") + self.channel = ext_spec.data.columns["CHANNEL"].array + col_spec_data = ext_spec.data.columns[self.spec_type] + self.spec_data = col_spec_data.array.copy() + self.spec_unit = col_spec_data.unit + self.spec_dtype = col_spec_data.dtype + self.spec_fits_format = col_spec_data.format + # grouping and quality + if "GROUPING" in colnames: + self.grouping = ext_spec.data.columns["GROUPING"].array + if "QUALITY" in colnames: + self.quality = ext_spec.data.columns["QUALITY"].array + # keywords self.EXPOSURE = self.header.get("EXPOSURE") self.BACKSCAL = self.header.get("BACKSCAL") @@ -511,8 +510,6 @@ class Spectrum: # {{{ self.RESPFILE = self.header.get("RESPFILE") self.ANCRFILE = self.header.get("ANCRFILE") self.BACKFILE = self.header.get("BACKFILE") - # output filename - self.outfile = outfile def get_data(self, group_squeeze=False, copy=True): """ @@ -558,7 +555,7 @@ class Spectrum: # {{{ group, which are estimated by utilizing the Monte Carlo techniques. """ self.stat_err = np.zeros(self.spec_data.shape, - dtype=self.spec_data.dtype) + dtype=self.spec_data.dtype) if group_squeeze and self.groupped: assert sum(self.grouping == 1) == len(stat_err) self.stat_err[self.grouping == 1] = stat_err @@ -627,13 +624,13 @@ class Spectrum: # {{{ if self.groupped: idx = self.grouping == 1 self.spec_data[idx] = np.random.normal(self.spec_data[idx], - self.spec_err[idx]) + self.spec_err[idx]) else: self.spec_data = np.random.normal(self.spec_data, self.spec_err) return self def fix_header_keywords(self, - reset_kw=["ANCRFILE", "RESPFILE", "BACKFILE"]): + reset_kw=["ANCRFILE", "RESPFILE", "BACKFILE"]): """ Reset the keywords to "NONE" to avoid confusion or mistakes, and also add mandatory spectral keywords if missing. @@ -644,36 +641,37 @@ class Spectrum: # {{{ https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/ogip_92_007_summary.html """ default_keywords = { - ## Mandatory keywords - #"EXTNAME" : "SPECTRUM", - "TELESCOP" : "NONE", - "INSTRUME" : "NONE", - "FILTER" : "NONE", - #"EXPOSURE" : <integration_time (s)>, - "BACKFILE" : "NONE", - "CORRFILE" : "NONE", - "CORRSCAL" : 1.0, - "RESPFILE" : "NONE", - "ANCRFILE" : "NONE", - "HDUCLASS" : "OGIP", - "HDUCLAS1" : "SPECTRUM", - "HDUVERS" : "1.2.1", - "POISSERR" : True, - #"CHANTYPE" : "PI", - #"DETCHANS" : <total_number_of_detector_channels>, - ## Optional keywords for further information - "BACKSCAL" : 1.0, - "AREASCAL" : 1.0, - # Type of spectral data: - # (1) "TOTAL": gross spectrum (source+bkg); - # (2) "NET": background-subtracted spectrum - # (3) "BKG" background spectrum - #"HDUCLAS2" : "NET", - # Details of the type of data: - # (1) "COUNT": data stored as counts - # (2) "RATE": data stored as counts/s - "HDUCLAS3" : { "COUNTS":"COUNT", - "RATE":"RATE" }.get(self.spec_type), + # Mandatory keywords + # "EXTNAME": "SPECTRUM", + "TELESCOP": "NONE", + "INSTRUME": "NONE", + "FILTER": "NONE", + # "EXPOSURE" : <integration_time (s)>, + "BACKFILE": "NONE", + "CORRFILE": "NONE", + "CORRSCAL": 1.0, + "RESPFILE": "NONE", + "ANCRFILE": "NONE", + "HDUCLASS": "OGIP", + "HDUCLAS1": "SPECTRUM", + "HDUVERS": "1.2.1", + "POISSERR": True, + # "CHANTYPE": "PI", + # "DETCHANS": <total_number_of_detector_channels>, + + # Optional keywords for further information + "BACKSCAL": 1.0, + "AREASCAL": 1.0, + # Type of spectral data: + # (1) "TOTAL": gross spectrum (source+bkg); + # (2) "NET": background-subtracted spectrum + # (3) "BKG" background spectrum + # "HDUCLAS2": "NET", + # Details of the type of data: + # (1) "COUNT": data stored as counts + # (2) "RATE": data stored as counts/s + "HDUCLAS3": {"COUNTS": "COUNT", + "RATE": "RATE"}.get(self.spec_type), } # add mandatory keywords if missing for kw, value in default_keywords.items(): @@ -683,13 +681,13 @@ class Spectrum: # {{{ for kw in reset_kw: self.header[kw] = default_keywords.get(kw) - def write(self, filename=None, clobber=False): + def write(self, outfile, clobber=False): """ Create a new "SPECTRUM" table/extension and replace the original one, then write to output file. """ - if filename is None: - filename = self.outfile + # Open the original input spectrum as the base + fitsobj = fits.open(self.filename) columns = [ fits.Column(name="CHANNEL", format="I", array=self.channel), fits.Column(name=self.spec_type, format=self.spec_fits_format, @@ -698,7 +696,7 @@ class Spectrum: # {{{ if self.grouping is not None: columns.append(fits.Column(name="GROUPING", format="I", array=self.grouping)) - if self.quality is not None: + if self.quality is not None: columns.append(fits.Column(name="QUALITY", format="I", array=self.quality)) if self.stat_err is not None: @@ -708,8 +706,9 @@ class Spectrum: # {{{ ext_spec_cols = fits.ColDefs(columns) ext_spec = fits.BinTableHDU.from_columns(ext_spec_cols, header=self.header) - self.fitsobj["SPECTRUM"] = ext_spec - self.fitsobj.writeto(filename, clobber=clobber, checksum=True) + # Replace the original spectrum data + fitsobj["SPECTRUM"] = ext_spec + fitsobj.writeto(outfile, clobber=clobber, checksum=True) # class Spectrum }}} @@ -738,17 +737,16 @@ class SpectrumSet(Spectrum): # {{{ # numpy dtype and FITS format code to which the spectrum data be # converted if the data is "COUNTS" - #_spec_dtype = np.float32 - #_spec_fits_format = "E" - _spec_dtype = np.float64 + _spec_dtype = np.float64 _spec_fits_format = "D" def __init__(self, filename, outfile=None, arf=None, rmf=None, bkg=None): - super().__init__(filename, outfile) + super().__init__(filename=filename) + self.outfile = outfile # convert spectrum data type if necessary if self.spec_data.dtype != self._spec_dtype: - self.spec_data = self.spec_data.astype(self._spec_dtype) - self.spec_dtype = self._spec_dtype + self.spec_data = self.spec_data.astype(self._spec_dtype) + self.spec_dtype = self._spec_dtype self.spec_fits_format = self._spec_fits_format if arf is not None: if isinstance(arf, ARF): @@ -771,6 +769,11 @@ class SpectrumSet(Spectrum): # {{{ self.bkg.spec_dtype = self._spec_dtype self.bkg.spec_fits_format = self._spec_fits_format + def write(self, outfile=None, clobber=False): + if outfile is None: + outfile = self.outfile + super().write(outfile=outfile, clobber=clobber) + def get_energy(self, mean="geometric"): """ Get the energy values of each channel if RMF present. @@ -824,7 +827,7 @@ class SpectrumSet(Spectrum): # {{{ while True: try: angle_begin = float(self.header["XFLT%04d" % num]) - angle_end = float(self.header["XFLT%04d" % (num+1)]) + angle_end = float(self.header["XFLT%04d" % (num+1)]) num += 2 except KeyError: break @@ -936,7 +939,7 @@ class SpectrumSet(Spectrum): # {{{ return spec_data_subbkg def subtract(self, spectrumset, cross_arf, groupped=False, - group_squeeze=False, add_history=False, verbose=False): + group_squeeze=False, add_history=False, verbose=False): """ Subtract the photons that originate from the surrounding regions but were scattered into this spectrum due to the finite PSF. @@ -952,30 +955,32 @@ class SpectrumSet(Spectrum): # {{{ spec1_new = spec1 - spec2 * (cross_arf_2_to_1 / arf2) 2. The ARF are interpolated to match the energies of spetral channels. """ - operation = " SUBTRACT: %s - (%s/%s) * %s" % (self.filename, - cross_arf.filename, spectrumset.arf.filename, - spectrumset.filename) + operation = " SUBTRACT: %s - (%s/%s) * %s" % ( + self.filename, cross_arf.filename, spectrumset.arf.filename, + spectrumset.filename) if verbose: print(operation, file=sys.stderr) energy = self.get_energy() if groupped: spectrumset.arf.apply_grouping(energy_channel=energy, - grouping=self.grouping, verbose=verbose) + grouping=self.grouping, + verbose=verbose) cross_arf.apply_grouping(energy_channel=energy, - grouping=self.grouping, verbose=verbose) - arfresp_spec = spectrumset.arf.get_data(groupped=True, - group_squeeze=group_squeeze) - arfresp_cross = cross_arf.get_data(groupped=True, - group_squeeze=group_squeeze) + grouping=self.grouping, + verbose=verbose) + arfresp_spec = spectrumset.arf.get_data( + groupped=True, group_squeeze=group_squeeze) + arfresp_cross = cross_arf.get_data( + groupped=True, group_squeeze=group_squeeze) else: - arfresp_spec = spectrumset.arf.interpolate(x=energy, - verbose=verbose) + arfresp_spec = spectrumset.arf.interpolate( + x=energy, verbose=verbose) arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose) with np.errstate(divide="ignore", invalid="ignore"): arf_ratio = arfresp_cross / arfresp_spec # fix nan/inf values due to division by zero - arf_ratio[ ~ np.isfinite(arf_ratio) ] = 0.0 - spec_data = (self.get_data(group_squeeze=group_squeeze) - + arf_ratio[~np.isfinite(arf_ratio)] = 0.0 + spec_data = (self.get_data(group_squeeze=group_squeeze) - spectrumset.get_data(group_squeeze=group_squeeze)*arf_ratio) self.set_data(spec_data, group_squeeze=group_squeeze) # record history @@ -983,7 +988,7 @@ class SpectrumSet(Spectrum): # {{{ self.header.add_history(operation) def compensate(self, cross_arf, groupped=False, group_squeeze=False, - add_history=False, verbose=False): + add_history=False, verbose=False): """ Compensate the photons that originate from this regions but were scattered into the surrounding regions due to the finite PSF. @@ -991,25 +996,27 @@ class SpectrumSet(Spectrum): # {{{ formula: spec1_new = spec1 + spec1 * (cross_arf_1_to_2 / arf1) """ - operation = " COMPENSATE: %s + (%s/%s) * %s" % (self.filename, - cross_arf.filename, self.arf.filename, self.filename) + operation = " COMPENSATE: %s + (%s/%s) * %s" % ( + self.filename, cross_arf.filename, self.arf.filename, + self.filename) if verbose: print(operation, file=sys.stderr) energy = self.get_energy() if groupped: cross_arf.apply_grouping(energy_channel=energy, - grouping=self.grouping, verbose=verbose) - arfresp_this = self.arf.get_data(groupped=True, - group_squeeze=group_squeeze) - arfresp_cross = cross_arf.get_data(groupped=True, - group_squeeze=group_squeeze) + grouping=self.grouping, + verbose=verbose) + arfresp_this = self.arf.get_data( + groupped=True, group_squeeze=group_squeeze) + arfresp_cross = cross_arf.get_data( + groupped=True, group_squeeze=group_squeeze) else: - arfresp_this = self.arf.interpolate(x=energy, verbose=verbose) + arfresp_this = self.arf.interpolate(x=energy, verbose=verbose) arfresp_cross = cross_arf.interpolate(x=energy, verbose=verbose) with np.errstate(divide="ignore", invalid="ignore"): arf_ratio = arfresp_cross / arfresp_this # fix nan/inf values due to division by zero - arf_ratio[ ~ np.isfinite(arf_ratio) ] = 0.0 + arf_ratio[~np.isfinite(arf_ratio)] = 0.0 spec_data = (self.get_data(group_squeeze=group_squeeze) + self.get_data(group_squeeze=group_squeeze) * arf_ratio) self.set_data(spec_data, group_squeeze=group_squeeze) @@ -1034,7 +1041,7 @@ class SpectrumSet(Spectrum): # {{{ if verbose: if i == 1: print("*** Fixing negative channels: iter %d..." % i, - end="", file=sys.stderr) + end="", file=sys.stderr) else: print("%d..." % i, end="", file=sys.stderr) for ch in neg_channels: @@ -1095,8 +1102,6 @@ class Crosstalk: # {{{ """ # `SpectrumSet' object for the spectrum to be corrected spectrumset = None - # NOTE/XXX: do NOT use list (e.g., []) here, otherwise, all the - # instances will share these list properties. # `SpectrumSet' and `ARF' objects corresponding to the spectra from # which the photons were scattered into this spectrum. cross_in_specset = None @@ -1111,7 +1116,7 @@ class Crosstalk: # {{{ groupped = False def __init__(self, config, arf_dict={}, rmf_dict={}, - grouping=None, quality=None): + grouping=None, quality=None): """ Arguments: * config: a section of the whole config file (`ConfigObj' object) @@ -1120,17 +1125,19 @@ class Crosstalk: # {{{ self.cross_in_arf = [] self.cross_out_arf = [] # this spectrum to be corrected - self.spectrumset = SpectrumSet(filename=config["spec"], - outfile=config["outfile"], - arf=arf_dict.get(config["arf"], config["arf"]), - rmf=rmf_dict.get(config.get("rmf"), config.get("rmf")), - bkg=config.get("bkg")) + self.spectrumset = SpectrumSet( + filename=config["spec"], + outfile=config["outfile"], + arf=arf_dict.get(config["arf"], config["arf"]), + rmf=rmf_dict.get(config.get("rmf"), config.get("rmf")), + bkg=config.get("bkg")) # spectra and cross arf from which photons were scattered in for reg_in in config["cross_in"].values(): - specset = SpectrumSet(filename=reg_in["spec"], - arf=arf_dict.get(reg_in["arf"], reg_in["arf"]), - rmf=rmf_dict.get(reg_in.get("rmf"), reg_in.get("rmf")), - bkg=reg_in.get("bkg")) + specset = SpectrumSet( + filename=reg_in["spec"], + arf=arf_dict.get(reg_in["arf"], reg_in["arf"]), + rmf=rmf_dict.get(reg_in.get("rmf"), reg_in.get("rmf")), + bkg=reg_in.get("bkg")) self.cross_in_specset.append(specset) self.cross_in_arf.append(arf_dict.get(reg_in["cross_arf"], ARF(reg_in["cross_arf"]))) @@ -1145,11 +1152,11 @@ class Crosstalk: # {{{ def apply_grouping(self, verbose=False): self.spectrumset.apply_grouping(grouping=self.grouping, - quality=self.quality, verbose=verbose) + quality=self.quality, verbose=verbose) # also group the related surrounding spectra for specset in self.cross_in_specset: specset.apply_grouping(grouping=self.grouping, - quality=self.quality, verbose=verbose) + quality=self.quality, verbose=verbose) self.groupped = True def estimate_errors(self, gehrels=True, verbose=False): @@ -1161,7 +1168,7 @@ class Crosstalk: # {{{ specset.estimate_errors(gehrels=gehrels) def do_correction(self, subtract_bkg=True, fix_negative=False, - group_squeeze=True, add_history=False, verbose=False): + group_squeeze=True, add_history=False, verbose=False): """ Perform the crosstalk correction. The background contribution for each spectrum is subtracted first if `subtract_bkg' is True. @@ -1177,33 +1184,36 @@ class Crosstalk: # {{{ if verbose: print("INFO: subtract background ...", file=sys.stderr) self.spectrumset.subtract_bkg(inplace=True, - add_history=add_history, verbose=verbose) + add_history=add_history, + verbose=verbose) # also apply background subtraction to the surrounding spectra for specset in self.cross_in_specset: specset.subtract_bkg(inplace=True, - add_history=add_history, verbose=verbose) + add_history=add_history, + verbose=verbose) # subtractions if verbose: print("INFO: apply subtractions ...", file=sys.stderr) for specset, cross_arf in zip(self.cross_in_specset, - self.cross_in_arf): - self.spectrumset.subtract(spectrumset=specset, - cross_arf=cross_arf, groupped=self.groupped, - group_squeeze=group_squeeze, add_history=add_history, - verbose=verbose) + self.cross_in_arf): + self.spectrumset.subtract( + spectrumset=specset, cross_arf=cross_arf, + groupped=self.groupped, group_squeeze=group_squeeze, + add_history=add_history, verbose=verbose) # compensations if verbose: print("INFO: apply compensations ...", file=sys.stderr) for cross_arf in self.cross_out_arf: - self.spectrumset.compensate(cross_arf=cross_arf, - groupped=self.groupped, group_squeeze=group_squeeze, - add_history=add_history, verbose=verbose) + self.spectrumset.compensate( + cross_arf=cross_arf, groupped=self.groupped, + group_squeeze=group_squeeze, add_history=add_history, + verbose=verbose) # fix negative values in channels if fix_negative: if verbose: print("INFO: fix negative channel values ...", file=sys.stderr) self.spectrumset.fix_negative(add_history=add_history, - verbose=verbose) + verbose=verbose) if add_history: self.spectrumset.header.add_history("END Crosstalk Correction") @@ -1255,10 +1265,6 @@ class Deprojection: # {{{ https://www-xray.ast.cam.ac.uk/papers/dsdeproj/ Sanders & Fabian 2007, MNRAS, 381, 1381 """ - spectra = None - grouping = None - quality = None - def __init__(self, spectra, grouping=None, quality=None, verbose=False): """ Arguments: @@ -1273,20 +1279,20 @@ class Deprojection: # {{{ raise ValueError("Not a 'SpectrumSet' object") spec.read_xflt() self.spectra.append(spec) - self.spectra = spectra + self.spectra = spectra self.grouping = grouping - self.quality = quality + self.quality = quality # sort spectra by `radius_outer' self.spectra.sort(key=lambda x: x.radius_outer) # set the inner radii - radii_inner = [0.0] + [ x.radius_outer for x in self.spectra[:-1] ] + radii_inner = [0.0] + [x.radius_outer for x in self.spectra[:-1]] for spec, rin in zip(self.spectra, radii_inner): spec.set_radius_inner(rin) if verbose: print("Deprojection: loaded spectrum: radius: (%s, %s)" % (spec.radius_inner, spec.radius_outer), file=sys.stderr) # check EXPOSURE validity (all spectra must have the same exposures) - exposures = [ spec.EXPOSURE for spec in self.spectra ] + exposures = [spec.EXPOSURE for spec in self.spectra] assert np.allclose(exposures[:-1], exposures[1:]) def subtract_bkg(self, verbose=True): @@ -1299,7 +1305,7 @@ class Deprojection: # {{{ def apply_grouping(self, verbose=False): for spec in self.spectra: spec.apply_grouping(grouping=self.grouping, quality=self.quality, - verbose=verbose) + verbose=verbose) def estimate_errors(self, gehrels=True): for spec in self.spectra: @@ -1313,7 +1319,7 @@ class Deprojection: # {{{ spec.scale() def do_deprojection(self, group_squeeze=True, - add_history=False, verbose=False): + add_history=False, verbose=False): # # TODO/XXX: How to apply ARF correction here??? # @@ -1326,7 +1332,7 @@ class Deprojection: # {{{ for shellnum in reversed(range(num_spec)): if verbose: print("DEPROJECTION: deprojecting shell %d ..." % shellnum, - file=sys.stderr) + file=sys.stderr) spec = self.spectra[shellnum] # calculate projected spectrum of outlying shells proj_spec = np.zeros(spec_shape, spec_dtype) @@ -1357,8 +1363,8 @@ class Deprojection: # {{{ Extract the spectral data of each spectrum after deprojection performed. """ - return [ spec.get_data(group_squeeze=group_squeeze, copy=copy) - for spec in self.spectra ] + return [spec.get_data(group_squeeze=group_squeeze, copy=copy) + for spec in self.spectra] def set_spec_data(self, spec_data, group_squeeze=True): """ @@ -1397,7 +1403,7 @@ class Deprojection: # {{{ Write the deprojected spectra to output file. """ if filenames == []: - filenames = [ spec.outfile for spec in self.spectra ] + filenames = [spec.outfile for spec in self.spectra] for spec, outfile in zip(self.spectra, filenames): spec.write(filename=outfile, clobber=clobber) @@ -1445,10 +1451,10 @@ def calc_median_errors(results): # sort by the Monte Carlo simulation axis results.sort(0) mc_times = results.shape[0] - medians = results[ int(mc_times * 0.5) ] - lowerpcs = results[ int(mc_times * 0.1585) ] - upperpcs = results[ int(mc_times * 0.8415) ] - errors = np.sqrt(0.5 * ((medians-lowerpcs)**2 + (upperpcs-medians)**2)) + medians = results[int(mc_times * 0.5)] + lowerpcs = results[int(mc_times * 0.1585)] + upperpcs = results[int(mc_times * 0.8415)] + errors = np.sqrt(0.5 * ((medians-lowerpcs)**2 + (upperpcs-medians)**2)) return (medians, errors) @@ -1465,7 +1471,7 @@ def set_argument(name, default, cmdargs, config): # main routine {{{ def main(config, subtract_bkg, fix_negative, mc_times, - verbose=False, clobber=False): + verbose=False, clobber=False): # collect ARFs and RMFs into dictionaries (avoid interpolation every time) arf_files = set() rmf_files = set() @@ -1480,9 +1486,9 @@ def main(config, subtract_bkg, fix_negative, mc_times, for arf in config_reg["cross_out"].as_list("cross_arf"): arf_files.add(arf) arf_files = arf_files - set([None]) - arf_dict = { arf: ARF(arf) for arf in arf_files } + arf_dict = {arf: ARF(arf) for arf in arf_files} rmf_files = rmf_files - set([None]) - rmf_dict = { rmf: RMF(rmf) for rmf in rmf_files } + rmf_dict = {rmf: RMF(rmf) for rmf in rmf_files} if verbose: print("INFO: arf_files:", arf_files, file=sys.stderr) print("INFO: rmf_files:", rmf_files, file=sys.stderr) @@ -1490,7 +1496,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, # get the GROUPING and QUALITY data grouping_fits = fits.open(config["grouping"]) grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array - quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array + quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array # squeeze the groupped spectral data, etc. group_squeeze = True @@ -1504,8 +1510,8 @@ def main(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: processing '%s' ..." % region, file=sys.stderr) crosstalk = Crosstalk(config.get(region), - arf_dict=arf_dict, rmf_dict=rmf_dict, - grouping=grouping, quality=quality) + arf_dict=arf_dict, rmf_dict=rmf_dict, + grouping=grouping, quality=quality) crosstalk.apply_grouping(verbose=verbose) crosstalk.estimate_errors(verbose=verbose) # keep a (almost) clean copy of the crosstalk object @@ -1513,24 +1519,26 @@ def main(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: doing crosstalk correction ...", file=sys.stderr) crosstalk.do_correction(subtract_bkg=subtract_bkg, - fix_negative=fix_negative, group_squeeze=group_squeeze, - add_history=True, verbose=verbose) + fix_negative=fix_negative, + group_squeeze=group_squeeze, + add_history=True, verbose=verbose) cc_spectra.append(crosstalk.get_spectrum(copy=True)) # load back the crosstalk-corrected spectra for deprojection if verbose: print("INFO: preparing spectra for deprojection ...", file=sys.stderr) deprojection = Deprojection(spectra=cc_spectra, grouping=grouping, - quality=quality, verbose=verbose) + quality=quality, verbose=verbose) if verbose: print("INFO: scaling spectra according the region angular size...", - file=sys.stderr) + file=sys.stderr) deprojection.scale() if verbose: print("INFO: doing deprojection ...", file=sys.stderr) deprojection.do_deprojection(add_history=True, verbose=verbose) - deproj_results = [ deprojection.get_spec_data( - group_squeeze=group_squeeze, copy=True) ] + deproj_results = [] + deproj_results.append(deprojection.get_spec_data( + group_squeeze=group_squeeze, copy=True)) # Monte Carlo for spectral group error estimation print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % @@ -1544,12 +1552,14 @@ def main(config, subtract_bkg, fix_negative, mc_times, # copy and randomize crosstalk_copy = crosstalk.copy().randomize() crosstalk_copy.do_correction(subtract_bkg=subtract_bkg, - fix_negative=fix_negative, group_squeeze=group_squeeze, - add_history=False, verbose=False) + fix_negative=fix_negative, + group_squeeze=group_squeeze, + add_history=False, verbose=False) cc_spectra_copy.append(crosstalk_copy.get_spectrum(copy=True)) # deproject spectra deprojection_copy = Deprojection(spectra=cc_spectra_copy, - grouping=grouping, quality=quality, verbose=False) + grouping=grouping, + quality=quality, verbose=False) deprojection_copy.scale() deprojection_copy.do_deprojection(add_history=False, verbose=False) deproj_results.append(deprojection_copy.get_spec_data( @@ -1558,7 +1568,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: Calculating the median and errors for each spectrum ...", - file=sys.stderr) + file=sys.stderr) medians, errors = calc_median_errors(deproj_results) deprojection.set_spec_data(medians, group_squeeze=group_squeeze) deprojection.add_stat_err(errors, group_squeeze=group_squeeze) @@ -1583,9 +1593,9 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): arf_files.add(config_reg.get("arf")) rmf_files.add(config_reg.get("rmf")) arf_files = arf_files - set([None]) - arf_dict = { arf: ARF(arf) for arf in arf_files } + arf_dict = {arf: ARF(arf) for arf in arf_files} rmf_files = rmf_files - set([None]) - rmf_dict = { rmf: RMF(rmf) for rmf in rmf_files } + rmf_dict = {rmf: RMF(rmf) for rmf in rmf_files} if verbose: print("INFO: arf_files:", arf_files, file=sys.stderr) print("INFO: rmf_files:", rmf_files, file=sys.stderr) @@ -1593,7 +1603,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): # get the GROUPING and QUALITY data grouping_fits = fits.open(config["grouping"]) grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array - quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array + quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array # squeeze the groupped spectral data, etc. group_squeeze = True @@ -1603,24 +1613,24 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): proj_spectra = [] for region in config.sections: config_reg = config[region] - specset = SpectrumSet(filename=config_reg["spec"], - outfile=config_reg["outfile"], - arf=arf_dict.get(config_reg["arf"], config_reg["arf"]), - rmf=rmf_dict.get(config_reg["rmf"], config_reg["rmf"]), - bkg=config_reg["bkg"]) + specset = SpectrumSet( + filename=config_reg["spec"], outfile=config_reg["outfile"], + arf=arf_dict.get(config_reg["arf"], config_reg["arf"]), + rmf=rmf_dict.get(config_reg["rmf"], config_reg["rmf"]), + bkg=config_reg["bkg"]) proj_spectra.append(specset) deprojection = Deprojection(spectra=proj_spectra, grouping=grouping, - quality=quality, verbose=verbose) + quality=quality, verbose=verbose) deprojection.apply_grouping(verbose=verbose) deprojection.estimate_errors() if verbose: print("INFO: scaling spectra according the region angular size ...", - file=sys.stderr) + file=sys.stderr) deprojection.scale() # keep a (almost) clean copy of the input projected spectra - proj_spectra_cleancopy = [ spec.copy() for spec in proj_spectra ] + proj_spectra_cleancopy = [spec.copy() for spec in proj_spectra] if verbose: print("INFO: subtract the background ...", file=sys.stderr) @@ -1628,8 +1638,9 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): if verbose: print("INFO: doing deprojection ...", file=sys.stderr) deprojection.do_deprojection(add_history=True, verbose=verbose) - deproj_results = [ deprojection.get_spec_data( - group_squeeze=group_squeeze, copy=True) ] + deproj_results = [] + deproj_results.append(deprojection.get_spec_data( + group_squeeze=group_squeeze, copy=True)) # Monte Carlo for spectral group error estimation print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % @@ -1638,11 +1649,12 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): if i % 100 == 0: print("%d..." % i, end="", flush=True, file=sys.stderr) # copy and randomize the input projected spectra - proj_spectra_copy = [ spec.copy().randomize() - for spec in proj_spectra_cleancopy ] + proj_spectra_copy = [spec.copy().randomize() + for spec in proj_spectra_cleancopy] # deproject spectra deprojection_copy = Deprojection(spectra=proj_spectra_copy, - grouping=grouping, quality=quality, verbose=False) + grouping=grouping, + quality=quality, verbose=False) deprojection_copy.subtract_bkg(verbose=False) deprojection_copy.do_deprojection(add_history=False, verbose=False) deproj_results.append(deprojection_copy.get_spec_data( @@ -1651,7 +1663,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): if verbose: print("INFO: Calculating the median and errors for each spectrum ...", - file=sys.stderr) + file=sys.stderr) medians, errors = calc_median_errors(deproj_results) deprojection.set_spec_data(medians, group_squeeze=group_squeeze) deprojection.add_stat_err(errors, group_squeeze=group_squeeze) @@ -1665,7 +1677,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): # main_crosstalk routine {{{ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, - verbose=False, clobber=False): + verbose=False, clobber=False): """ Only perform the crosstalk correction. """ @@ -1683,9 +1695,9 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, for arf in config_reg["cross_out"].as_list("cross_arf"): arf_files.add(arf) arf_files = arf_files - set([None]) - arf_dict = { arf: ARF(arf) for arf in arf_files } + arf_dict = {arf: ARF(arf) for arf in arf_files} rmf_files = rmf_files - set([None]) - rmf_dict = { rmf: RMF(rmf) for rmf in rmf_files } + rmf_dict = {rmf: RMF(rmf) for rmf in rmf_files} if verbose: print("INFO: arf_files:", arf_files, file=sys.stderr) print("INFO: rmf_files:", rmf_files, file=sys.stderr) @@ -1694,11 +1706,11 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, if "grouping" in config.keys(): grouping_fits = fits.open(config["grouping"]) grouping = grouping_fits["SPECTRUM"].data.columns["GROUPING"].array - quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array + quality = grouping_fits["SPECTRUM"].data.columns["QUALITY"].array group_squeeze = True else: grouping = None - quality = None + quality = None group_squeeze = False # crosstalk objects (BEFORE background subtraction) @@ -1711,8 +1723,8 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: processing '%s' ..." % region, file=sys.stderr) crosstalk = Crosstalk(config.get(region), - arf_dict=arf_dict, rmf_dict=rmf_dict, - grouping=grouping, quality=quality) + arf_dict=arf_dict, rmf_dict=rmf_dict, + grouping=grouping, quality=quality) if grouping is not None: crosstalk.apply_grouping(verbose=verbose) crosstalk.estimate_errors(verbose=verbose) @@ -1721,15 +1733,16 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: doing crosstalk correction ...", file=sys.stderr) crosstalk.do_correction(subtract_bkg=subtract_bkg, - fix_negative=fix_negative, group_squeeze=group_squeeze, - add_history=True, verbose=verbose) + fix_negative=fix_negative, + group_squeeze=group_squeeze, + add_history=True, verbose=verbose) crosstalk.fix_header() cc_spectra.append(crosstalk.get_spectrum(copy=True)) # spectral data of the crosstalk-corrected spectra cc_results = [] - cc_results.append([ spec.get_data(group_squeeze=group_squeeze, copy=True) - for spec in cc_spectra ]) + cc_results = ([spec.get_data(group_squeeze=group_squeeze, copy=True) + for spec in cc_spectra]) # Monte Carlo for spectral group error estimation print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % @@ -1743,17 +1756,18 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, # copy and randomize crosstalk_copy = crosstalk.copy().randomize() crosstalk_copy.do_correction(subtract_bkg=subtract_bkg, - fix_negative=fix_negative, group_squeeze=group_squeeze, - add_history=False, verbose=False) + fix_negative=fix_negative, + group_squeeze=group_squeeze, + add_history=False, verbose=False) cc_spectra_copy.append(crosstalk_copy.get_spectrum(copy=True)) - cc_results.append([ spec.get_data(group_squeeze=group_squeeze, - copy=True) - for spec in cc_spectra_copy ]) + cc_results.append([spec.get_data(group_squeeze=group_squeeze, + copy=True) + for spec in cc_spectra_copy]) print("DONE!", flush=True, file=sys.stderr) if verbose: print("INFO: Calculating the median and errors for each spectrum ...", - file=sys.stderr) + file=sys.stderr) medians, errors = calc_median_errors(cc_results) if verbose: print("INFO: Writing the crosstalk-corrected spectra " + @@ -1807,16 +1821,16 @@ if __name__ == "__main__": if mode.lower() == "both": print("MODE: CROSSTALK + DEPROJECTION", file=sys.stderr) main(config, subtract_bkg=subtract_bkg, fix_negative=fix_negative, - mc_times=mc_times, verbose=verbose, clobber=clobber) + mc_times=mc_times, verbose=verbose, clobber=clobber) elif mode.lower() == "deprojection": print("MODE: DEPROJECTION", file=sys.stderr) main_deprojection(config, mc_times=mc_times, - verbose=verbose, clobber=clobber) + verbose=verbose, clobber=clobber) elif mode.lower() == "crosstalk": print("MODE: CROSSTALK", file=sys.stderr) main_crosstalk(config, subtract_bkg=subtract_bkg, - fix_negative=fix_negative, mc_times=mc_times, - verbose=verbose, clobber=clobber) + fix_negative=fix_negative, mc_times=mc_times, + verbose=verbose, clobber=clobber) else: raise ValueError("Invalid operation mode: %s" % mode) print(WARNING) |