crosstalk_deprojection.py: v0.6.0, report spectral data changes; many fixes

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":