diff options
| -rwxr-xr-x | astro/spectrum/crosstalk_deprojection.py | 361 |
1 files changed, 195 insertions, 166 deletions
diff --git a/astro/spectrum/crosstalk_deprojection.py b/astro/spectrum/crosstalk_deprojection.py index a7aa30a..8dd2a03 100755 --- a/astro/spectrum/crosstalk_deprojection.py +++ b/astro/spectrum/crosstalk_deprojection.py @@ -5,10 +5,12 @@ # # Change log: # 2017-11-12: -# * Various cleanups +# * Report spectrum data changes +# * Add 'regid' to identify spectrum/arf/rmf regions # * Add 'info_emin' and 'info_emax' configuration options -# * Remove command line arguments # * Add global configurations variable 'CONFIGS' +# * Remove command line arguments +# * Various cleanups # 2016-06-07: # * Explain the errors/uncertainties calculation approach # 2016-04-20: @@ -154,7 +156,7 @@ from astropy.io import fits from configobj import ConfigObj -__version__ = "0.5.6" +__version__ = "0.6.0" __date__ = "2017-11-12" WARNING = """ @@ -184,7 +186,7 @@ def group_data(data, grouping): else: # other channels of a group data_grp[i-1] += data_grp[i] - data_grp[i] = 0 + data_grp[i] = 0 assert np.isclose(sum(data_grp), sum(data)) return data_grp @@ -220,30 +222,31 @@ class ARF: # {{{ [2] Definition of RMF and ARF file formats https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html """ - filename = None - # only consider the "SPECTRUM" extension - header = None - energ_lo = None - energ_hi = None - specresp = None # function of the interpolated ARF f_interp = None # energies of the spectral channels energy_channel = None # spectral channel grouping specification - grouping = None - groupped = False + grouping = None + groupped = False # groupped ARF channels with respect to the grouping - specresp_grp = None + specresp_grp = None - def __init__(self, filename): + def __init__(self, filename, regid=None): self.filename = filename + self.regid = regid 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"] + self.energ_lo = ext_specresp.data["ENERG_LO"] # [keV] + self.energ_hi = ext_specresp.data["ENERG_HI"] # [keV] + self.specresp = ext_specresp.data["SPECRESP"] # [cm^2] + + def __str__(self): + if self.regid: + return self.regid + else: + return self.filename def get_data(self, groupped=False, group_squeeze=False, copy=True): if groupped: @@ -257,26 +260,26 @@ class ARF: # {{{ else: return specresp - def get_energy(self, mean="geometric"): + @property + def energy(self): """ Return the mean energy values of the ARF. + The geometric mean is adopted, i.e., e = sqrt(e_min*e_max) + Unit: [keV] + """ + return np.sqrt(self.energ_lo * self.energ_hi) - Arguments: - * mean: type of the mean energy: - + "geometric": geometric mean, i.e., e = sqrt(e_min*e_max) - + "arithmetic": arithmetic mean, i.e., e = 0.5*(e_min+e_max) - """ - if mean == "geometric": - energy = np.sqrt(self.energ_lo * self.energ_hi) - elif mean == "arithmetic": - energy = 0.5 * (self.energ_lo + self.energ_hi) + @property + def energy_groupped(self): + if self.groupped: + energy = self.energy + return energy[self.grouping == 1] else: - raise ValueError("Invalid mean type: %s" % mean) - return energy + return None def interpolate(self, x=None, verbose=False): """ - Cubic interpolate the ARF curve using `scipy.interpolate' + Interpolate the ARF curve using `scipy.interpolate' If the requested point is outside of the data range, the fill value of *zero* is returned. @@ -289,13 +292,12 @@ class ARF: # {{{ otherwise, the interpolated data are returned. """ if not hasattr(self, "f_interp") or self.f_interp is None: - energy = self.get_energy() arf = self.get_data(copy=False) if verbose: - print("INFO: interpolating ARF '%s' (may take a while) ..." % + 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, + self.energy, arf, kind="quadratic", bounds_error=False, fill_value=0.0, assume_sorted=True) self.f_interp = f_interp @@ -318,7 +320,7 @@ class ARF: # {{{ if self.groupped: return if verbose: - print("INFO: Grouping spectrum '%s' ..." % self.filename, + print("INFO: Grouping ARF '%s' ..." % self.filename, file=sys.stderr) self.energy_channel = energy_channel self.grouping = grouping @@ -363,33 +365,19 @@ class RMF: # {{{ [2] Definition of RMF and ARF file formats https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html """ - filename = None - ## extension "MATRIX" - hdr_matrix = None - energ_lo = None - energ_hi = None - n_grp = None - f_chan = None - n_chan = None - # raw squeezed RMF matrix data - matrix = None - ## extension "EBOUNDS" - hdr_ebounds = None - channel = None - e_min = None - e_max = None - ## converted 2D RMF matrix/image from the squeezed binary table - # size: len(energ_lo) x len(channel) - rmfimg = None - - def __init__(self, filename): + # converted 2D RMF matrix/image from the squeezed binary table + # size: len(energ_lo) x len(channel) + rmfimg = None + + def __init__(self, filename, regid=None): self.filename = filename + self.regid = regid 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.energ_lo = ext_matrix.data["ENERG_LO"] # [keV] + self.energ_hi = ext_matrix.data["ENERG_HI"] # [keV] self.n_grp = ext_matrix.data["N_GRP"] self.f_chan = ext_matrix.data["F_CHAN"] self.n_chan = ext_matrix.data["N_CHAN"] @@ -398,25 +386,23 @@ class RMF: # {{{ 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"] + self.e_min = ext_ebounds.data["E_MIN"] # [keV] + self.e_max = ext_ebounds.data["E_MAX"] # [keV] + + def __str__(self): + if self.regid: + return self.regid + else: + return self.filename - def get_energy(self, mean="geometric"): + @property + def energy(self): """ Return the mean energy values of the RMF "EBOUNDS". - - Arguments: - * mean: type of the mean energy: - + "geometric": geometric mean, i.e., e = sqrt(e_min*e_max) - + "arithmetic": arithmetic mean, i.e., e = 0.5*(e_min+e_max) - """ - if mean == "geometric": - energy = np.sqrt(self.e_min * self.e_max) - elif mean == "arithmetic": - energy = 0.5 * (self.e_min + self.e_max) - else: - raise ValueError("Invalid mean type: %s" % mean) - return energy + The geometric mean is adopted, i.e., e = sqrt(e_min*e_max) + Unit: [keV] + """ + return np.sqrt(self.e_min * self.e_max) def get_rmfimg(self): """ @@ -436,7 +422,7 @@ class RMF: # {{{ # if self.rmfimg is None: # Make the 2D RMF matrix/image - n_energy = len(self.energ_lo) + n_energy = len(self.energ_lo) n_channel = len(self.channel) rmf_dtype = self.matrix[0].dtype rmfimg = np.zeros(shape=(n_energy, n_channel), dtype=rmf_dtype) @@ -454,7 +440,10 @@ class RMF: # {{{ header = self.hdr_matrix.copy(strip=True) header.extend(self.hdr_ebounds.copy(strip=True)) outfits = fits.PrimaryHDU(data=rmfimg, header=header) - outfits.writeto(outfile, checksum=True, clobber=clobber) + try: + outfits.writeto(outfile, checksum=True, overwrite=clobber) + except TypeError: + outfits.writeto(outfile, checksum=True, clobber=clobber) # class RMF }}} @@ -729,7 +718,10 @@ class Spectrum: # {{{ header=self.header) # Replace the original spectrum data fitsobj["SPECTRUM"] = ext_spec - fitsobj.writeto(outfile, clobber=clobber, checksum=True) + try: + fitsobj.writeto(outfile, overwrite=clobber, checksum=True) + except TypeError: + fitsobj.writeto(outfile, clobber=clobber, checksum=True) # class Spectrum }}} @@ -796,7 +788,8 @@ class SpectrumSet(Spectrum): # {{{ outfile = self.outfile super().write(outfile=outfile, clobber=clobber) - def get_energy(self, mean="geometric"): + @property + def energy(self): """ Get the energy values of each channel if RMF present. @@ -810,7 +803,15 @@ class SpectrumSet(Spectrum): # {{{ if self.rmf is None: return None else: - return self.rmf.get_energy(mean=mean) + return self.rmf.energy + + @property + def energy_groupped(self): + energy = self.energy + if self.groupped and energy is not None: + return energy[self.grouping == 1] + else: + return None def get_arf(self, groupped=True, copy=True): """ @@ -827,7 +828,7 @@ class SpectrumSet(Spectrum): # {{{ else: return self.arf.get_data(groupped=groupped, copy=copy) - def read_xflt(self): + def read_xflt(self, verbose=False): """ Read the XFLT000# keywords from the header, check the validity (e.g., "XFLT0001" should equals "XFLT0002", "XFLT0003" should equals 0). @@ -857,8 +858,10 @@ class SpectrumSet(Spectrum): # {{{ # if NO additional XFLT000# keys exist, assume "annulus" region if self.angle < eps: self.angle = 360.0 + if verbose: + print("[%s] region angle: %.2f [deg]" % (str(self), self.angle)) - def scale(self): + def scale(self, verbose=False): """ Scale the spectral data (and spectral group errors if present) of the source spectrum (and background spectra if present) according @@ -870,14 +873,17 @@ class SpectrumSet(Spectrum): # {{{ The spectral data and errors (i.e., `self.spec_data', and `self.spec_err') is MODIFIED! """ - self.spec_data *= (360.0 / self.angle) + factor = 360.0 / self.angle + if verbose: + print("[%s] region scale factor: %.2f" % (str(self), factor)) + self.spec_data *= factor if self.spec_err is not None: - self.spec_err *= (360.0 / self.angle) + self.spec_err *= factor # also scale the background spectrum if present if self.bkg: - self.bkg.spec_data *= (360.0 / self.angle) + self.bkg.spec_data *= factor if self.bkg.spec_err is not None: - self.bkg.spec_err *= (360.0 / self.angle) + self.bkg.spec_err *= factor def apply_grouping(self, grouping=None, quality=None, verbose=False): """ @@ -894,7 +900,7 @@ class SpectrumSet(Spectrum): # {{{ """ super().apply_grouping(grouping=grouping, quality=quality) # also group the ARF accordingly - self.arf.apply_grouping(energy_channel=self.get_energy(), + self.arf.apply_grouping(energy_channel=self.energy, grouping=self.grouping, verbose=verbose) # group the background spectrum if present if self.bkg: @@ -944,12 +950,16 @@ class SpectrumSet(Spectrum): # {{{ (self.filename, ratio, self.bkg.filename)) if verbose: print(operation, file=sys.stderr) - datasum1 = np.sum(self.spec_data) + energy = self.energy + eidx = ((energy >= CONFIGS["info_emin"]) & + (energy <= CONFIGS["info_emax"])) + datasum1 = np.sum(self.spec_data[eidx]) spec_data_subbkg = self.spec_data - ratio * self.bkg.get_data() - datasum2 = np.sum(spec_data_subbkg) + datasum2 = np.sum(spec_data_subbkg[eidx]) datacp = 100 * (datasum2-datasum1) / datasum1 - print("[%s] bkg subtraction: %g -> %g (%.2f%%)" % - (str(self), datasum1, datasum2, datacp)) + if verbose: + print("[%s] bkg subtraction: %g -> %g (%.2f%%)" % + (str(self), datasum1, datasum2, datacp)) if inplace: self.spec_data = spec_data_subbkg self.spec_bkg_subtracted = True @@ -987,7 +997,7 @@ class SpectrumSet(Spectrum): # {{{ spectrumset.filename) if verbose: print(operation, file=sys.stderr) - energy = self.get_energy() + energy = self.energy if groupped: spectrumset.arf.apply_grouping(energy_channel=energy, grouping=self.grouping, @@ -1007,9 +1017,19 @@ class SpectrumSet(Spectrum): # {{{ 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) - - spectrumset.get_data(group_squeeze=group_squeeze)*arf_ratio) - self.set_data(spec_data, group_squeeze=group_squeeze) + spec_data1 = self.get_data(group_squeeze=group_squeeze) + spec_data2 = spectrumset.get_data(group_squeeze=group_squeeze) + spec_data_new = spec_data1 - spec_data2 * arf_ratio + self.set_data(spec_data_new, group_squeeze=group_squeeze) + energy_groupped = self.energy_groupped + eidx = ((energy_groupped >= CONFIGS["info_emin"]) & + (energy_groupped <= CONFIGS["info_emax"])) + datasum1 = np.sum(spec_data1[eidx]) + datasum2 = np.sum(spec_data_new[eidx]) + datacp = 100 * (datasum2-datasum1) / datasum1 + if verbose: + print("[%s] crosstalk/subtract -> [%s]: %g -> %g (%.2f%%)" % + (str(self), str(spectrumset), datasum1, datasum2, datacp)) # record history if add_history: self.header.add_history(operation) @@ -1028,7 +1048,7 @@ class SpectrumSet(Spectrum): # {{{ self.filename) if verbose: print(operation, file=sys.stderr) - energy = self.get_energy() + energy = self.energy if groupped: cross_arf.apply_grouping(energy_channel=energy, grouping=self.grouping, @@ -1044,9 +1064,18 @@ class SpectrumSet(Spectrum): # {{{ arf_ratio = arfresp_cross / arfresp_this # 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) + - self.get_data(group_squeeze=group_squeeze) * arf_ratio) - self.set_data(spec_data, group_squeeze=group_squeeze) + spec_data = self.get_data(group_squeeze=group_squeeze) + spec_data_new = spec_data * (1 + arf_ratio) + self.set_data(spec_data_new, group_squeeze=group_squeeze) + energy_groupped = self.energy_groupped + eidx = ((energy_groupped >= CONFIGS["info_emin"]) & + (energy_groupped <= CONFIGS["info_emax"])) + datasum1 = np.sum(spec_data[eidx]) + datasum2 = np.sum(spec_data_new[eidx]) + datacp = 100 * (datasum2-datasum1) / datasum1 + if verbose: + print("[%s] crosstalk/compensate <- [%s]: %g -> %g (%.2f%%)" % + (str(self), str(cross_arf), datasum1, datasum2, datacp)) # record history if add_history: self.header.add_history(operation) @@ -1132,7 +1161,7 @@ class Crosstalk: # {{{ # `SpectrumSet' and `ARF' objects corresponding to the spectra from # which the photons were scattered into this spectrum. cross_in_specset = None - cross_in_arf = None + cross_in_arf = None # `ARF' objects corresponding to the regions to which the photons of # this spectrum were scattered into. cross_out_arf = None @@ -1142,32 +1171,33 @@ class Crosstalk: # {{{ # whether the spectrum is groupped groupped = False - def __init__(self, config, arf_dict={}, rmf_dict={}, + def __init__(self, config, regid=None, arf_dict={}, rmf_dict={}, grouping=None, quality=None): """ Arguments: * config: a section of the whole config file (`ConfigObj' object) """ + self.regid = regid self.cross_in_specset = [] - self.cross_in_arf = [] - self.cross_out_arf = [] + self.cross_in_arf = [] + self.cross_out_arf = [] # this spectrum to be corrected self.spectrumset = SpectrumSet( - filename=config["spec"], + filename=config["spec"], regid=regid, outfile=config["outfile"], - arf=arf_dict.get(config["arf"], config["arf"]), - rmf=rmf_dict.get(config.get("rmf"), config.get("rmf")), + arf=arf_dict[config["arf"]], + rmf=rmf_dict.get(config.get("rmf"), None), bkg=config.get("bkg")) # spectra and cross arf from which photons were scattered in - for reg_in in config["cross_in"].values(): + for k, v in config["cross_in"].items(): 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")) + filename=v["spec"], regid=k, + arf=arf_dict[config["arf"]], + rmf=rmf_dict.get(config.get("rmf"), None), + bkg=v.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"]))) + self.cross_in_arf.append( + arf_dict.get(v["cross_arf"], ARF(v["cross_arf"], regid=k))) # regions into which the photons of this spectrum were scattered into if "cross_out" in config.sections: cross_arf = config["cross_out"].as_list("cross_arf") @@ -1304,7 +1334,7 @@ class Deprojection: # {{{ for spec in spectra: if not isinstance(spec, SpectrumSet): raise ValueError("Not a 'SpectrumSet' object") - spec.read_xflt() + spec.read_xflt(verbose=verbose) self.spectra.append(spec) self.spectra = spectra self.grouping = grouping @@ -1338,12 +1368,12 @@ class Deprojection: # {{{ for spec in self.spectra: spec.estimate_errors(gehrels=gehrels) - def scale(self): + def scale(self, verbose=False): """ Scale the spectral data according to the region angular size. """ for spec in self.spectra: - spec.scale() + spec.scale(verbose=verbose) def do_deprojection(self, group_squeeze=True, add_history=False, verbose=False): @@ -1425,14 +1455,14 @@ class Deprojection: # {{{ spec.fix_header_keywords( reset_kw=["RESPFILE", "ANCRFILE", "BACKFILE"]) - def write(self, filenames=[], clobber=False): + def write(self, outfiles=None, clobber=False): """ Write the deprojected spectra to output file. """ - if filenames == []: - filenames = [spec.outfile for spec in self.spectra] - for spec, outfile in zip(self.spectra, filenames): - spec.write(filename=outfile, clobber=clobber) + if outfiles is None: + outfiles = [spec.outfile for spec in self.spectra] + for spec, outfile in zip(self.spectra, outfiles): + spec.write(outfile=outfile, clobber=clobber) @staticmethod def projected_volume(r1, r2, R1, R2): @@ -1488,27 +1518,26 @@ def calc_median_errors(results): # main routine {{{ def main(config, subtract_bkg, fix_negative, mc_times, - verbose=False, clobber=False): + verbose=False, clobber=False, **kwargs): # collect ARFs and RMFs into dictionaries (avoid interpolation every time) arf_files = set() rmf_files = set() for region in config.sections: config_reg = config[region] - arf_files.add(config_reg.get("arf")) - rmf_files.add(config_reg.get("rmf")) - for reg_in in config_reg["cross_in"].values(): - arf_files.add(reg_in.get("arf")) - arf_files.add(reg_in.get("cross_arf")) + arf_files.add((config_reg["arf"], region)) + rmf_files.add((config_reg["rmf"], region)) + for k, v in config_reg["cross_in"].items(): + arf_files.add((v["arf"], k)) + arf_files.add((v["cross_arf"], k)) if "cross_out" in config_reg.sections: 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} - rmf_files = rmf_files - set([None]) - 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) + arf_files.add((arf, None)) + print("INFO: arf_files:", arf_files, file=sys.stderr) + print("INFO: rmf_files:", rmf_files, file=sys.stderr) + arf_dict = {arf: ARF(arf, regid=reg) + for arf, reg in arf_files if arf is not None} + rmf_dict = {rmf: RMF(rmf, regid=reg) + for rmf, reg in rmf_files if rmf is not None} # get the GROUPING and QUALITY data grouping_fits = fits.open(config["grouping"]) @@ -1526,7 +1555,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, for region in config.sections: if verbose: print("INFO: processing '%s' ..." % region, file=sys.stderr) - crosstalk = Crosstalk(config.get(region), + crosstalk = Crosstalk(config.get(region), regid=region, arf_dict=arf_dict, rmf_dict=rmf_dict, grouping=grouping, quality=quality) crosstalk.apply_grouping(verbose=verbose) @@ -1549,7 +1578,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, if verbose: print("INFO: scaling spectra according the region angular size...", file=sys.stderr) - deprojection.scale() + deprojection.scale(verbose=verbose) if verbose: print("INFO: doing deprojection ...", file=sys.stderr) deprojection.do_deprojection(add_history=True, verbose=verbose) @@ -1561,7 +1590,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % mc_times, file=sys.stderr) for i in range(mc_times): - if i % 100 == 0: + if i % 50 == 0: print("%d..." % i, end="", flush=True, file=sys.stderr) # correct crosstalk effects cc_spectra_copy = [] @@ -1577,7 +1606,7 @@ def main(config, subtract_bkg, fix_negative, mc_times, deprojection_copy = Deprojection(spectra=cc_spectra_copy, grouping=grouping, quality=quality, verbose=False) - deprojection_copy.scale() + deprojection_copy.scale(verbose=False) deprojection_copy.do_deprojection(add_history=False, verbose=False) deproj_results.append(deprojection_copy.get_spec_data( group_squeeze=group_squeeze, copy=True)) @@ -1598,7 +1627,8 @@ def main(config, subtract_bkg, fix_negative, mc_times, # main_deprojection routine {{{ -def main_deprojection(config, mc_times, verbose=False, clobber=False): +def main_deprojection(config, mc_times, verbose=False, clobber=False, + **kwargs): """ Only perform the spectral deprojection. """ @@ -1607,15 +1637,14 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): rmf_files = set() for region in config.sections: config_reg = config[region] - 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} - rmf_files = rmf_files - set([None]) - 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) + arf_files.add((config_reg["arf"], region)) + rmf_files.add((config_reg["rmf"], region)) + print("INFO: arf_files:", arf_files, file=sys.stderr) + print("INFO: rmf_files:", rmf_files, file=sys.stderr) + arf_dict = {arf: ARF(arf, regid=reg) + for arf, reg in arf_files if arf is not None} + rmf_dict = {rmf: RMF(rmf, regid=reg) + for rmf, reg in rmf_files if rmf is not None} # get the GROUPING and QUALITY data grouping_fits = fits.open(config["grouping"]) @@ -1631,7 +1660,8 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): for region in config.sections: config_reg = config[region] specset = SpectrumSet( - filename=config_reg["spec"], outfile=config_reg["outfile"], + filename=config_reg["spec"], regid=region, + 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"]) @@ -1644,7 +1674,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): if verbose: print("INFO: scaling spectra according the region angular size ...", file=sys.stderr) - deprojection.scale() + deprojection.scale(verbose=verbose) # keep a (almost) clean copy of the input projected spectra proj_spectra_cleancopy = [spec.copy() for spec in proj_spectra] @@ -1663,7 +1693,7 @@ def main_deprojection(config, mc_times, verbose=False, clobber=False): print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % mc_times, file=sys.stderr) for i in range(mc_times): - if i % 100 == 0: + if i % 50 == 0: print("%d..." % i, end="", flush=True, file=sys.stderr) # copy and randomize the input projected spectra proj_spectra_copy = [spec.copy().randomize() @@ -1694,7 +1724,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, **kwargs): """ Only perform the crosstalk correction. """ @@ -1703,21 +1733,20 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, rmf_files = set() for region in config.sections: config_reg = config[region] - arf_files.add(config_reg.get("arf")) - rmf_files.add(config_reg.get("rmf")) - for reg_in in config_reg["cross_in"].values(): - arf_files.add(reg_in.get("arf")) - arf_files.add(reg_in.get("cross_arf")) + arf_files.add((config_reg["arf"], region)) + rmf_files.add((config_reg["rmf"], region)) + for k, v in config_reg["cross_in"].items(): + arf_files.add((v["arf"], k)) + arf_files.add((v["cross_arf"], k)) if "cross_out" in config_reg.sections: 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} - rmf_files = rmf_files - set([None]) - 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) + arf_files.add((arf, None)) + print("INFO: arf_files:", arf_files, file=sys.stderr) + print("INFO: rmf_files:", rmf_files, file=sys.stderr) + arf_dict = {arf: ARF(arf, regid=reg) + for arf, reg in arf_files if arf is not None} + rmf_dict = {rmf: RMF(rmf, regid=reg) + for rmf, reg in rmf_files if rmf is not None} # get the GROUPING and QUALITY data if "grouping" in config.keys(): @@ -1739,7 +1768,7 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, for region in config.sections: if verbose: print("INFO: processing '%s' ..." % region, file=sys.stderr) - crosstalk = Crosstalk(config.get(region), + crosstalk = Crosstalk(config.get(region), regid=region, arf_dict=arf_dict, rmf_dict=rmf_dict, grouping=grouping, quality=quality) if grouping is not None: @@ -1765,7 +1794,7 @@ def main_crosstalk(config, subtract_bkg, fix_negative, mc_times, print("INFO: Monte Carlo to estimate spectral errors (%d times) ..." % mc_times, file=sys.stderr) for i in range(mc_times): - if i % 100 == 0: + if i % 50 == 0: print("%d..." % i, end="", flush=True, file=sys.stderr) # correct crosstalk effects cc_spectra_copy = [] @@ -1847,7 +1876,7 @@ if __name__ == "__main__": info_emax = 7.0 CONFIGS.update(info_emin=info_emin, info_emax=info_emax) - for k, v in CONFIGS: + for k, v in CONFIGS.items(): print("* %s: %s" % (k, v)) if mode.lower() == "both": |