From 01102115eb1ea40bdc3377372f4ee94f372506f9 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Tue, 20 Aug 2024 15:39:18 +0200 Subject: [PATCH 01/11] Possible to use Poisson error in spectrum --- ogip/spec.py | 543 ++++++++++++++++++++++++--------------------------- 1 file changed, 251 insertions(+), 292 deletions(-) diff --git a/ogip/spec.py b/ogip/spec.py index a7528e2..f718815 100644 --- a/ogip/spec.py +++ b/ogip/spec.py @@ -1,25 +1,33 @@ import logging -import numpy as np # type: ignore -import astropy.io.fits as fits # type: ignore +import numpy as np # type: ignore +import astropy.io.fits as fits # type: ignore logger = logging.getLogger() + def log_bins(n_bins_per_decade, global_e_min, global_e_max): n_bin_edges_per_decade = n_bins_per_decade + 1 - new_e_bins = np.logspace( - np.log10(global_e_min), - np.log10(global_e_max), - int(np.log10(global_e_max / global_e_min) * n_bin_edges_per_decade), - ) + # if new_e_min is None: + # new_e_min = e_min[0] + # logger.info("choosing new e_min: %s", new_e_min) + + # if new_e_max is None: + # n_bins_in_range = np.ceil(np.log10(e_max[-1] / new_e_min) * n_bin_edges_per_decade) + # new_e_max = new_e_min * 10**(n_bins_in_range/n_bin_edges_per_decade) + # logger.info("bins in nrange: %s", n_bins_in_range) + # logger.info("choosing new e_max: %s", new_e_max) + + new_e_bins = np.logspace(np.log10(global_e_min), np.log10(global_e_max), int(np.log10(global_e_max/global_e_min)*n_bin_edges_per_decade)) logger.debug("new e bins: %s", new_e_bins) return new_e_bins - + class Spectrum: + @staticmethod def from_file_name(fn): raise NotImplementedError @@ -29,24 +37,19 @@ def to_long_string(self): ### - class Readable: @staticmethod def from_file_name(cls, fn): methods_tried = {} for n, m in cls.__dict__.items(): - if n.startswith("from_file_name_"): + if n.startswith('from_file_name_'): try: return m.__func__(fn) except Exception as e: - logger.debug( - "failed to read with %s.%s %s: %s", cls.__name__, n, m, e - ) + logger.debug("failed to read with %s.%s %s: %s", cls.__name__, n, m, e) methods_tried[n] = e - raise Exception( - f"failed to read with any method, tried: {methods_tried}") - + raise Exception(f"failed to read with any method, tried: {methods_tried}") class PHAI(Spectrum): _rate = None @@ -67,12 +70,12 @@ def from_file_name_osa(fn): for e in f: try: return PHAI.from_arrays( - exposure=e.header["EXPOSURE"], - rate=e.data["RATE"], - stat_err=e.data["STAT_ERR"], - sys_err=e.data["SYS_ERR"], - filename=fn, - ) + exposure=e.header['EXPOSURE'], + rate=e.data['RATE'], + stat_err=e.data['STAT_ERR'], + sys_err=e.data['SYS_ERR'], + filename=fn + ) except Exception as ex: logger.debug("failed to read from %s: %s", e, ex) @@ -80,25 +83,44 @@ def from_file_name_osa(fn): @staticmethod def from_file_name_normal(fn): + f = fits.open(fn) + if 'RATE' in f['SPECTRUM'].data.names: + rate = f['SPECTRUM'].data['RATE'] + elif 'COUNTS' in f['SPECTRUM'].data.names: + rate = f['SPECTRUM'].data['COUNTS'] / f['SPECTRUM'].header['EXPOSURE'] + else: + rate = None + + if 'STAT_ERR' in f['SPECTRUM'].data.names: + stat_err = f['SPECTRUM'].data['STAT_ERR'] + elif 'POISSERR' in f['SPECTRUM'].header and f['SPECTRUM'].header['POISSERR']: + if 'COUNTS' in f['SPECTRUM'].data.names: + stat_err = np.sqrt(f['SPECTRUM'].data['COUNTS']) / f['SPECTRUM'].header['EXPOSURE'] + elif 'RATE' in f['SPECTRUM'].data.names: + stat_err = np.sqrt(f['SPECTRUM'].data['RATE'] * f['SPECTRUM'].header['EXPOSURE'])/f['SPECTRUM'].header['EXPOSURE'] + else: + stat_err = None + else: + stat_err= None + + if 'SYS_ERR' in f['SPECTRUM'].data.names: + sys_err = f['SPECTRUM'].data['SYS_ERR'] + else: + sys_err = None + return PHAI.from_arrays( - exposure=f["SPECTRUM"].header["EXPOSURE"], - rate=f["SPECTRUM"].data["RATE"], - stat_err=f["SPECTRUM"].data["STAT_ERR"], - sys_err=f["SPECTRUM"].data["SYS_ERR"], - filename=fn, - ) + exposure=f['SPECTRUM'].header['EXPOSURE'], + rate=rate, + stat_err = stat_err, + sys_err=sys_err, + filename=fn + ) + + @staticmethod - def from_arrays( - exposure, - rate=None, - stat_err=None, - sys_err=None, - quality=None, - counts=None, - filename=None, - ): + def from_arrays(exposure, rate=None, stat_err=None, sys_err=None, quality=None, counts=None, filename=None): self = PHAI() self.filename = filename @@ -106,13 +128,11 @@ def from_arrays( if counts is not None: logger.warning("counts found: converting to rate") - rate = counts / exposure + rate = counts/exposure if stat_err is None: - stat_err = counts**0.5 - logger.warning( - "assuming poisson stat errors. Is it really what you want?" - ) + _stat_err = counts**0.5 + logger.warning("assuming poisson stat errors. Is it really what you want?") if rate is None: raise Exception("need rate or counts") @@ -139,59 +159,42 @@ def from_arrays( @property def spectrum_hdu(self): - return fits.BinTableHDU.from_columns( - [ - fits.Column(name="RATE", array=self._rate, format="1E"), - fits.Column(name="STAT_ERR", - array=self._stat_err, format="1E"), - fits.Column(name="SYS_ERR", array=self._sys_err, format="1E"), - fits.Column(name="QUALITY", array=self._quality, format="1I"), - ], - # https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html - header=fits.Header( - cards=dict( - # - the name (i.e. type) of the extension - EXTNAME="SPECTRUM", - # - the "telescope" (i.e. mission/satellite name). - TELESCOP="", - INSTRUME="", # - the instrument/detector. - # FILTER - the instrument filter in use (if any) - AREASCAL=1.0, - BACKSCAL=1.0, - # the integration time (in seconds) for the PHA data (assumed to be corrected for deadtime, data drop-outs etc. ) - EXPOSURE=self._exposure, - # - the name of the corresponding background file (if any) - BACKFILE="", - # CORRFILE - the name of the corresponding correction file (if any) - # CORRSCAL - the correction scaling factor. - # - the name of the corresponding (default) redistribution matrix file (RMF; see George et al. 1992a). - RESPFILE="", - # - the name of the corresponding (default) ancillary response file (ARF; see George et al. 1992a). - ANCRFILE="", - # - should contain the string "OGIP" to indicate that this is an OGIP style file. - HDUCLASS="OGIP", - # - should contain the string "SPECTRUM" to indicate this is a spectrum. - HDUCLAS1="SPECTRUM", - # - the version number of the format (this document describes version 1.2.1) - HDUVERS="1.2.1", - # POISSERR #- whether Poissonian errors are appropriate to the data (see below). # noqa: E800 - # CHANTYPE #- whether the channels used in the file have been corrected in anyway (see below). # noqa: E800 - DETCHANS=len( - self._rate - ), # - the total number of detector channels available. - ) - ), - ) + return fits.BinTableHDU.from_columns([ + fits.Column(name='RATE', array=self._rate, format='1E'), + fits.Column(name='STAT_ERR', array=self._stat_err, format='1E'), + fits.Column(name='SYS_ERR', array=self._sys_err, format='1E'), + fits.Column(name='QUALITY', array=self._quality, format='1I'), + ], + # https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html + header=fits.Header(cards=dict( + EXTNAME="SPECTRUM", # - the name (i.e. type) of the extension + TELESCOP="", # - the "telescope" (i.e. mission/satellite name). + INSTRUME="", #- the instrument/detector. + #FILTER - the instrument filter in use (if any) + AREASCAL=1., + BACKSCAL=1., + EXPOSURE=self._exposure, # the integration time (in seconds) for the PHA data (assumed to be corrected for deadtime, data drop-outs etc. ) + BACKFILE="", #- the name of the corresponding background file (if any) + #CORRFILE - the name of the corresponding correction file (if any) + #CORRSCAL - the correction scaling factor. + RESPFILE="", # - the name of the corresponding (default) redistribution matrix file (RMF; see George et al. 1992a). + ANCRFILE="", # - the name of the corresponding (default) ancillary response file (ARF; see George et al. 1992a). + HDUCLASS="OGIP", # - should contain the string "OGIP" to indicate that this is an OGIP style file. + HDUCLAS1="SPECTRUM", # - should contain the string "SPECTRUM" to indicate this is a spectrum. + HDUVERS="1.2.1", #- the version number of the format (this document describes version 1.2.1) + #POISSERR #- whether Poissonian errors are appropriate to the data (see below). + #CHANTYPE #- whether the channels used in the file have been corrected in anyway (see below). + DETCHANS=len(self._rate), #- the total number of detector channels available. + )), + ) def to_fits(self, fn: str): logger.info("store to fits %s %s", self, fn) - fits.HDUList( - [ + fits.HDUList([ fits.PrimaryHDU(), self.spectrum_hdu, - ] - ).writeto(fn, overwrite=True) + ]).writeto(fn, overwrite=True) @property def total_counts(self): @@ -206,208 +209,173 @@ class PHAII(Spectrum): class RMF: - _telescop = "not-a-telescope" - _instrume = "not-an-instrument" + _telescop="not-a-telescope" + _instrume="not-an-instrument" - _energ_lo = None # type: ignore - _energ_hi = None # type: ignore - _matrix = None # type: ignore - _e_min = None # type: ignore - _e_max = None # type: ignore + _energ_lo = None # type: ignore + _energ_hi = None # type: ignore + _matrix = None # type: ignore + _e_min = None # type: ignore + _e_max = None # type: ignore + def to_long_string(self): return f"{self.__class__.__name__}: {len(self._e_min)} x {len(self._energ_lo)} channels " @staticmethod def from_file_name(fn): return Readable.from_file_name(RMF, fn) - + @staticmethod def from_file_name_osaisgri(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f["ISGR-RMF.-RSP"].data["ENERG_LO"], - energ_hi=f["ISGR-RMF.-RSP"].data["ENERG_HI"], - matrix=np.stack(f["ISGR-RMF.-RSP"].data["MATRIX"]), - e_min=f["ISGR-EBDS-MOD"].data["E_MIN"], - e_max=f["ISGR-EBDS-MOD"].data["E_MAX"], - ) - + energ_lo=f['ISGR-RMF.-RSP'].data['ENERG_LO'], + energ_hi=f['ISGR-RMF.-RSP'].data['ENERG_HI'], + matrix=np.stack(f['ISGR-RMF.-RSP'].data['MATRIX']), + e_min=f['ISGR-EBDS-MOD'].data['E_MIN'], + e_max=f['ISGR-EBDS-MOD'].data['E_MAX'], + ) + @staticmethod def from_file_name_osaspi(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f["SPI.-RMF.-RSP"].data["ENERG_LO"], - energ_hi=f["SPI.-RMF.-RSP"].data["ENERG_HI"], - matrix=np.stack(f["SPI.-RMF.-RSP"].data["MATRIX"]), - e_min=f["SPI.-EBDS-SET"].data["E_MIN"], - e_max=f["SPI.-EBDS-SET"].data["E_MAX"], - ) - + energ_lo=f['SPI.-RMF.-RSP'].data['ENERG_LO'], + energ_hi=f['SPI.-RMF.-RSP'].data['ENERG_HI'], + matrix=np.stack(f['SPI.-RMF.-RSP'].data['MATRIX']), + e_min=f['SPI.-EBDS-SET'].data['E_MIN'], + e_max=f['SPI.-EBDS-SET'].data['E_MAX'], + ) + @staticmethod def from_file_name_normal(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f["MATRIX"].data["ENERG_LO"], - energ_hi=f["MATRIX"].data["ENERG_HI"], - matrix=np.vstack(f["MATRIX"].data["MATRIX"]), - e_min=f["EBOUNDS"].data["E_MIN"], - e_max=f["EBOUNDS"].data["E_MAX"], - ) + energ_lo=f['MATRIX'].data['ENERG_LO'], + energ_hi=f['MATRIX'].data['ENERG_HI'], + matrix=np.vstack(f['MATRIX'].data['MATRIX']), + e_min=f['EBOUNDS'].data['E_MIN'], + e_max=f['EBOUNDS'].data['E_MAX'], + ) @staticmethod def from_file_name_alt_normal(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f["SPECRESP MATRIX"].data["ENERG_LO"], - energ_hi=f["SPECRESP MATRIX"].data["ENERG_HI"], - matrix=np.vstack(f["SPECRESP MATRIX"].data["MATRIX"]), - e_min=f["EBOUNDS"].data["E_MIN"], - e_max=f["EBOUNDS"].data["E_MAX"], - ) + energ_lo=f['SPECRESP MATRIX'].data['ENERG_LO'], + energ_hi=f['SPECRESP MATRIX'].data['ENERG_HI'], + matrix=np.vstack(f['SPECRESP MATRIX'].data['MATRIX']), + e_min=f['EBOUNDS'].data['E_MIN'], + e_max=f['EBOUNDS'].data['E_MAX'], + ) @staticmethod def from_arrays(energ_lo, energ_hi, matrix, e_min, e_max): self = RMF() - if not (len(energ_lo) == len(energ_hi) == matrix.shape[0]): - raise Exception( - f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {matrix.shape[0]}!" - ) + if not (len(energ_lo) == len(energ_hi) == matrix.shape[0]): + raise Exception(f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {matrix.shape[0]}!") + if not (len(e_min) == len(e_max) == matrix.shape[1]): - raise Exception( - f"incompatible dimensions of channel energy, bounds {len(e_min)} {len(e_max)} but matrix {matrix.shape[1]}!" - ) + raise Exception(f"incompatible dimensions of channel energy, bounds {len(e_min)} {len(e_max)} but matrix {matrix.shape[1]}!") - self._energ_lo = energ_lo # type: ignore - self._energ_hi = energ_hi # type: ignore - self._matrix = matrix # type: ignore - self._e_min = e_min # type: ignore - self._e_max = e_max # type: ignore + self._energ_lo = energ_lo # type: ignore + self._energ_hi = energ_hi # type: ignore + self._matrix = matrix # type: ignore + self._e_min = e_min # type: ignore + self._e_max = e_max # type: ignore return self @property def matrix_hdu(self): - return fits.BinTableHDU.from_columns( - [ - fits.Column(name="ENERG_LO", - array=self._energ_lo, format="1E"), - fits.Column(name="ENERG_HI", - array=self._energ_hi, format="1E"), - fits.Column( - name="N_GRP", array=np.ones_like(self._energ_lo), format="1I" - ), - fits.Column( - name="F_CHAN", array=0 * np.ones_like(self._energ_lo), format="1I" - ), - fits.Column( - name="N_CHAN", - array=len(self._e_min) * np.ones_like(self._energ_lo), - format="1I", - ), - fits.Column(name="MATRIX", array=self._matrix, format="PE"), - # fits.Column(name='MATRIX', array=self._matrix, format=f'{len(self._e_min)}E'), # noqa: E800 - ], - header=fits.Header( - cards=dict( + return fits.BinTableHDU.from_columns([ + fits.Column(name='ENERG_LO', array=self._energ_lo, format='1E'), + fits.Column(name='ENERG_HI', array=self._energ_hi, format='1E'), + fits.Column(name='N_GRP', array=np.ones_like(self._energ_lo), format='1I'), + fits.Column(name='F_CHAN', array=0*np.ones_like(self._energ_lo), format='1I'), + fits.Column(name='N_CHAN', array=len(self._e_min)*np.ones_like(self._energ_lo), format='1I'), + fits.Column(name='MATRIX', array=self._matrix, format='PE'), + #fits.Column(name='MATRIX', array=self._matrix, format=f'{len(self._e_min)}E'), + ], + header=fits.Header(cards=dict( # https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc7.1.1 - EXTNAME="MATRIX", # / name of this binary table extension - # BITPIX = 8 / 8-bit bytes - # NAXIS = 2 / 2-dimensional binary table - # NAXIS1 = 34 / width of table in bytes - # NAXIS2 = 1180 / number of rows in table - # PCOUNT = 1031160 / Number of bytes acumulated in heap - # GCOUNT = 1 / one data group (required keyword) - # TFIELDS = 6 / number of fields in each row - # TTYPE1 = 'ENERG_LO' / label for field 1 - # TFORM1 = 'E ' / data format of field: 4-byte REAL - # TUNIT1 = 'keV ' / physical unit of field - # TTYPE2 = 'ENERG_HI' / label for field 2 - # TFORM2 = 'E ' / data format of field: 4-byte REAL - # TUNIT2 = 'keV ' / physical unit of field - # TTYPE3 = 'N_GRP ' / label for field 3 - # TFORM3 = 'I ' / data format of field: 2-byte INTEGER - # TTYPE4 = 'F_CHAN ' / label for field 4 - # TFORM4 = 'PI(2) ' / data format of field: variable length array - # TTYPE5 = 'N_CHAN ' / label for field 5 - # TFORM5 = 'PI(2) ' / data format of field: variable length array - # TTYPE6 = 'MATRIX ' / label for field 6 - # TFORM6 = 'PE(418) ' / data format of field: variable length array - TLMIN4=0, # / First legal channel number - TLMAX4=len(self._e_min), # / Highest legal channel number - TELESCOP=self._telescop, # / mission/satellite name - INSTRUME=self._instrume, # / instrument/detector - # FILTER = 'NONE ' / filter information - CHANTYPE="PI", # / Type of channels (PHA, PI etc) - DETCHANS=len( - self._e_min - ), # / Total number of detector PHA channels - LO_THRES=1.00e-07, # / Lower probability density threshold for matrix - HDUCLASS="OGIP", # /KeywordinformationforCaltoolsSoftware. - # /KeywordinformationforCaltoolsSoftware. - HDUCLAS1="RESPONSE", - # /KeywordinformationforCaltoolsSoftware. - HDUCLAS2="RSP_MATRIX", - HDUVERS="1.3.0", # /KeywordinformationforCaltoolsSoftware. - # /KeywordinformationforCaltoolsSoftware. - HDUCLAS3="DETECTOR", - ) - ), - ) - + EXTNAME = 'MATRIX', # / name of this binary table extension + #BITPIX = 8 / 8-bit bytes + #NAXIS = 2 / 2-dimensional binary table + #NAXIS1 = 34 / width of table in bytes + #NAXIS2 = 1180 / number of rows in table + #PCOUNT = 1031160 / Number of bytes acumulated in heap + #GCOUNT = 1 / one data group (required keyword) + #TFIELDS = 6 / number of fields in each row + #TTYPE1 = 'ENERG_LO' / label for field 1 + #TFORM1 = 'E ' / data format of field: 4-byte REAL + #TUNIT1 = 'keV ' / physical unit of field + #TTYPE2 = 'ENERG_HI' / label for field 2 + #TFORM2 = 'E ' / data format of field: 4-byte REAL + #TUNIT2 = 'keV ' / physical unit of field + #TTYPE3 = 'N_GRP ' / label for field 3 + #TFORM3 = 'I ' / data format of field: 2-byte INTEGER + #TTYPE4 = 'F_CHAN ' / label for field 4 + #TFORM4 = 'PI(2) ' / data format of field: variable length array + #TTYPE5 = 'N_CHAN ' / label for field 5 + #TFORM5 = 'PI(2) ' / data format of field: variable length array + #TTYPE6 = 'MATRIX ' / label for field 6 + #TFORM6 = 'PE(418) ' / data format of field: variable length array + TLMIN4=0, #/ First legal channel number + TLMAX4=len(self._e_min), #/ Highest legal channel number + TELESCOP=self._telescop, # / mission/satellite name + INSTRUME=self._instrume, # / instrument/detector + #FILTER = 'NONE ' / filter information + CHANTYPE='PI', # / Type of channels (PHA, PI etc) + DETCHANS=len(self._e_min), # / Total number of detector PHA channels + LO_THRES=1.00E-07, # / Lower probability density threshold for matrix + HDUCLASS='OGIP', #/KeywordinformationforCaltoolsSoftware. + HDUCLAS1='RESPONSE', #/KeywordinformationforCaltoolsSoftware. + HDUCLAS2='RSP_MATRIX', #/KeywordinformationforCaltoolsSoftware. + HDUVERS='1.3.0', #/KeywordinformationforCaltoolsSoftware. + HDUCLAS3='DETECTOR', #/KeywordinformationforCaltoolsSoftware. + )), + ) + @property def ebounds_hdu(self): - return fits.BinTableHDU.from_columns( - [ - fits.Column( - name="CHANNEL", array=np.arange(self._e_min.shape[0]), format="1I" - ), - fits.Column(name="E_MIN", array=self._e_min, format="1E"), - fits.Column(name="E_MAX", array=self._e_max, format="1E"), - ], - header=fits.Header( - cards=dict( - EXTNAME="EBOUNDS", # / name of this binary table extension - TLMIN1=0, # / First legal channel number - TLMAX1=len( - self._e_min - ), # 511 / Highest legal channel number - TELESCOP=self._telescop, # / mission/satellite name - INSTRUME=self._instrume, # / instrument/detector - # FILTER = 'NONE ' / filter information - CHANTYPE="PI", # / Type of channels (PHA, PI etc) - DETCHANS=len( - self._e_min - ), # / Total number of detector PHA channels - # SMOOTHED= 0 / 0 = raw, 1-12 = smooth, -1 = ep-lin, -2 = mean- - # / Keyword information for Caltools Software. - HDUCLASS="OGIP", - # / Keyword information for Caltools Software. - HDUCLAS1="RESPONSE", - # / Keyword information for Caltools Software. - HDUCLAS2="EBOUNDS", - # / Keyword information for Caltools Software. - HDUVERS="1.2.0", - ) - ), - ) - + return fits.BinTableHDU.from_columns([ + fits.Column(name='CHANNEL', array=np.arange(self._e_min.shape[0]), format='1I'), + fits.Column(name='E_MIN', array=self._e_min, format='1E'), + fits.Column(name='E_MAX', array=self._e_max, format='1E'), + ], + header=fits.Header(cards=dict( + EXTNAME='EBOUNDS', # / name of this binary table extension + TLMIN1=0, #/ First legal channel number + TLMAX1=len(self._e_min), # 511 / Highest legal channel number + TELESCOP=self._telescop, # / mission/satellite name + INSTRUME=self._instrume, # / instrument/detector + #FILTER = 'NONE ' / filter information + CHANTYPE='PI', # / Type of channels (PHA, PI etc) + DETCHANS=len(self._e_min), # / Total number of detector PHA channels + #SMOOTHED= 0 / 0 = raw, 1-12 = smooth, -1 = ep-lin, -2 = mean- + HDUCLASS='OGIP', # / Keyword information for Caltools Software. + HDUCLAS1='RESPONSE', # / Keyword information for Caltools Software. + HDUCLAS2='EBOUNDS', # / Keyword information for Caltools Software. + HDUVERS ='1.2.0', # / Keyword information for Caltools Software. + )), + ) + def to_fits(self, fn): logger.info("store to fits %s %s", self, fn) - fits.HDUList( - [ + fits.HDUList([ fits.PrimaryHDU(), self.ebounds_hdu, self.matrix_hdu, - ] - ).writeto(fn, overwrite=True) + ]).writeto(fn, overwrite=True) @property def d_e(self): @@ -419,105 +387,96 @@ def d_energ(self): @property def c_e(self): - return (self._e_max + self._e_min) / 2 + return (self._e_max + self._e_min)/2 @property def c_energ(self): - return (self._energ_hi + self._energ_lo) / 2 + return (self._energ_hi + self._energ_lo)/2 + def rebin(pha: PHAI, rmf: RMF, new_e_bins) -> (PHAI, RMF): - new_e_bins_assigned = [] + new_e_bins_assigned = [] if pha._rate.shape[0] != rmf._matrix[0].shape[0]: raise RuntimeError() for new_bin_req_e1, new_bin_req_e2 in zip(new_e_bins[:-1], new_e_bins[1:]): m = rmf._e_min >= new_bin_req_e1 - m &= rmf._e_min < new_bin_req_e2 # really e_min here - - new_e_bins_assigned.append( - dict( - mask=m, - matrix_row=rmf._matrix[:, m].sum(1), - e_min=np.array(rmf._e_min)[m].min(), - e_max=np.array(rmf._e_max)[m].max(), - rate=pha._rate[m].sum(), - stat_err=np.sum(pha._stat_err[m] ** 2) ** 0.5, - sys_err=np.sum(pha._sys_err[m] ** 2) ** 0.5, - ) - ) - - logger.info( - "for %s - %s new bin assigned: %s", - new_bin_req_e1, - new_bin_req_e2, - {k: v for k, - v in new_e_bins_assigned[-1].items() if len(str(v)) < 50}, - ) + m &= rmf._e_min < new_bin_req_e2 # really e_min here + + new_e_bins_assigned.append(dict( + mask=m, + matrix_row=rmf._matrix[:, m].sum(1), + e_min=np.array(rmf._e_min)[m].min(), + e_max=np.array(rmf._e_max)[m].max(), + rate=pha._rate[m].sum(), + stat_err=np.sum(pha._stat_err[m]**2)**0.5, + sys_err=np.sum(pha._sys_err[m]**2)**0.5, + )) + + logger.info("for %s - %s new bin assigned: %s", new_bin_req_e1, new_bin_req_e2, {k: v for k, v in new_e_bins_assigned[-1].items() if len(str(v)) < 50}) return ( PHAI.from_arrays( exposure=pha._exposure, - rate=np.array([r["rate"] for r in new_e_bins_assigned]), - stat_err=np.array([r["stat_err"] for r in new_e_bins_assigned]), - sys_err=np.array([r["stat_err"] for r in new_e_bins_assigned]), + rate=np.array([r['rate'] for r in new_e_bins_assigned]), + stat_err=np.array([r['stat_err'] for r in new_e_bins_assigned]), + sys_err=np.array([r['stat_err'] for r in new_e_bins_assigned]), ), - RMF.from_arrays( + RMF.from_arrays( energ_lo=rmf._energ_lo.copy(), energ_hi=rmf._energ_hi.copy(), - matrix=np.vstack( - [r["matrix_row"] for r in new_e_bins_assigned] - ).transpose(), - e_min=np.array([r["e_min"] for r in new_e_bins_assigned]), - e_max=np.array([r["e_max"] for r in new_e_bins_assigned]), + matrix=np.vstack([r['matrix_row'] for r in new_e_bins_assigned]).transpose(), + e_min=np.array([r['e_min'] for r in new_e_bins_assigned]), + e_max=np.array([r['e_max'] for r in new_e_bins_assigned]), ), ) + + + class ARF: _arf = None @staticmethod def from_file_name(fn): return Readable.from_file_name(ARF, fn) - + @staticmethod def from_file_name_normal(fn): f = fits.open(fn) return ARF.from_arrays( - energ_lo=f["SPECRESP"].data["ENERG_LO"], - energ_hi=f["SPECRESP"].data["ENERG_HI"], - arf=f["SPECRESP"].data["SPECRESP"], - ) + energ_lo=f['SPECRESP'].data['ENERG_LO'], + energ_hi=f['SPECRESP'].data['ENERG_HI'], + arf=f['SPECRESP'].data['SPECRESP'], + ) @staticmethod def from_file_name_osa(fn): f = fits.open(fn) return ARF.from_arrays( - energ_lo=f["ISGR-ARF.-RSP"].data["ENERG_LO"], - energ_hi=f["ISGR-ARF.-RSP"].data["ENERG_HI"], - arf=f["ISGR-ARF.-RSP"].data["SPECRESP"], - ) + energ_lo=f['ISGR-ARF.-RSP'].data['ENERG_LO'], + energ_hi=f['ISGR-ARF.-RSP'].data['ENERG_HI'], + arf=f['ISGR-ARF.-RSP'].data['SPECRESP'], + ) @staticmethod def from_arrays(energ_lo, energ_hi, arf): self = ARF() if not (len(energ_lo) == len(energ_hi) == len(arf)): - raise Exception( - f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {len(arf)}!" - ) - - self._energ_lo = energ_lo # type: ignore - self._energ_hi = energ_hi # type: ignore - self._arf = arf # type: ignore + raise Exception(f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {len(arf)}!") + + self._energ_lo = energ_lo # type: ignore + self._energ_hi = energ_hi # type: ignore + self._arf = arf # type: ignore return self def to_long_string(self): - return ( - f"{self.__class__.__name__}: {len(self._arf)} energies {np.max(self._arf)}" - ) + return f"{self.__class__.__name__}: {len(self._arf)} energies {np.max(self._arf)}" + From 6dde2dbe9a8289eee3e49447b0190fd3b9e62415 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Tue, 20 Aug 2024 18:09:48 +0200 Subject: [PATCH 02/11] merged from master --- ogip/spec.py | 524 +++++++++++++++++++++++++++++---------------------- 1 file changed, 295 insertions(+), 229 deletions(-) diff --git a/ogip/spec.py b/ogip/spec.py index f718815..8d33783 100644 --- a/ogip/spec.py +++ b/ogip/spec.py @@ -1,33 +1,25 @@ import logging -import numpy as np # type: ignore -import astropy.io.fits as fits # type: ignore +import numpy as np # type: ignore +import astropy.io.fits as fits # type: ignore logger = logging.getLogger() - def log_bins(n_bins_per_decade, global_e_min, global_e_max): n_bin_edges_per_decade = n_bins_per_decade + 1 - # if new_e_min is None: - # new_e_min = e_min[0] - # logger.info("choosing new e_min: %s", new_e_min) - - # if new_e_max is None: - # n_bins_in_range = np.ceil(np.log10(e_max[-1] / new_e_min) * n_bin_edges_per_decade) - # new_e_max = new_e_min * 10**(n_bins_in_range/n_bin_edges_per_decade) - # logger.info("bins in nrange: %s", n_bins_in_range) - # logger.info("choosing new e_max: %s", new_e_max) - - new_e_bins = np.logspace(np.log10(global_e_min), np.log10(global_e_max), int(np.log10(global_e_max/global_e_min)*n_bin_edges_per_decade)) + new_e_bins = np.logspace( + np.log10(global_e_min), + np.log10(global_e_max), + int(np.log10(global_e_max / global_e_min) * n_bin_edges_per_decade), + ) logger.debug("new e bins: %s", new_e_bins) return new_e_bins - -class Spectrum: +class Spectrum: @staticmethod def from_file_name(fn): raise NotImplementedError @@ -37,19 +29,24 @@ def to_long_string(self): ### + class Readable: @staticmethod def from_file_name(cls, fn): methods_tried = {} for n, m in cls.__dict__.items(): - if n.startswith('from_file_name_'): + if n.startswith("from_file_name_"): try: return m.__func__(fn) except Exception as e: - logger.debug("failed to read with %s.%s %s: %s", cls.__name__, n, m, e) + logger.debug( + "failed to read with %s.%s %s: %s", cls.__name__, n, m, e + ) methods_tried[n] = e - raise Exception(f"failed to read with any method, tried: {methods_tried}") + raise Exception( + f"failed to read with any method, tried: {methods_tried}") + class PHAI(Spectrum): _rate = None @@ -70,12 +67,12 @@ def from_file_name_osa(fn): for e in f: try: return PHAI.from_arrays( - exposure=e.header['EXPOSURE'], - rate=e.data['RATE'], - stat_err=e.data['STAT_ERR'], - sys_err=e.data['SYS_ERR'], - filename=fn - ) + exposure=e.header["EXPOSURE"], + rate=e.data["RATE"], + stat_err=e.data["STAT_ERR"], + sys_err=e.data["SYS_ERR"], + filename=fn, + ) except Exception as ex: logger.debug("failed to read from %s: %s", e, ex) @@ -83,7 +80,7 @@ def from_file_name_osa(fn): @staticmethod def from_file_name_normal(fn): - + f = fits.open(fn) if 'RATE' in f['SPECTRUM'].data.names: rate = f['SPECTRUM'].data['RATE'] @@ -98,11 +95,11 @@ def from_file_name_normal(fn): if 'COUNTS' in f['SPECTRUM'].data.names: stat_err = np.sqrt(f['SPECTRUM'].data['COUNTS']) / f['SPECTRUM'].header['EXPOSURE'] elif 'RATE' in f['SPECTRUM'].data.names: - stat_err = np.sqrt(f['SPECTRUM'].data['RATE'] * f['SPECTRUM'].header['EXPOSURE'])/f['SPECTRUM'].header['EXPOSURE'] + stat_err = np.sqrt(f['SPECTRUM'].data['RATE'] * f['SPECTRUM'].header['EXPOSURE']) / f['SPECTRUM'].header['EXPOSURE'] else: stat_err = None else: - stat_err= None + stat_err = None if 'SYS_ERR' in f['SPECTRUM'].data.names: sys_err = f['SPECTRUM'].data['SYS_ERR'] @@ -110,17 +107,23 @@ def from_file_name_normal(fn): sys_err = None return PHAI.from_arrays( - exposure=f['SPECTRUM'].header['EXPOSURE'], - rate=rate, - stat_err = stat_err, - sys_err=sys_err, - filename=fn - ) - - + exposure=f['SPECTRUM'].header['EXPOSURE'], + rate=rate, + stat_err=stat_err, + sys_err=sys_err, + filename=fn + ) @staticmethod - def from_arrays(exposure, rate=None, stat_err=None, sys_err=None, quality=None, counts=None, filename=None): + def from_arrays( + exposure, + rate=None, + stat_err=None, + sys_err=None, + quality=None, + counts=None, + filename=None, + ): self = PHAI() self.filename = filename @@ -128,11 +131,13 @@ def from_arrays(exposure, rate=None, stat_err=None, sys_err=None, quality=None, if counts is not None: logger.warning("counts found: converting to rate") - rate = counts/exposure + rate = counts / exposure if stat_err is None: - _stat_err = counts**0.5 - logger.warning("assuming poisson stat errors. Is it really what you want?") + stat_err = counts**0.5 + logger.warning( + "assuming poisson stat errors. Is it really what you want?" + ) if rate is None: raise Exception("need rate or counts") @@ -159,42 +164,59 @@ def from_arrays(exposure, rate=None, stat_err=None, sys_err=None, quality=None, @property def spectrum_hdu(self): - return fits.BinTableHDU.from_columns([ - fits.Column(name='RATE', array=self._rate, format='1E'), - fits.Column(name='STAT_ERR', array=self._stat_err, format='1E'), - fits.Column(name='SYS_ERR', array=self._sys_err, format='1E'), - fits.Column(name='QUALITY', array=self._quality, format='1I'), - ], - # https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html - header=fits.Header(cards=dict( - EXTNAME="SPECTRUM", # - the name (i.e. type) of the extension - TELESCOP="", # - the "telescope" (i.e. mission/satellite name). - INSTRUME="", #- the instrument/detector. - #FILTER - the instrument filter in use (if any) - AREASCAL=1., - BACKSCAL=1., - EXPOSURE=self._exposure, # the integration time (in seconds) for the PHA data (assumed to be corrected for deadtime, data drop-outs etc. ) - BACKFILE="", #- the name of the corresponding background file (if any) - #CORRFILE - the name of the corresponding correction file (if any) - #CORRSCAL - the correction scaling factor. - RESPFILE="", # - the name of the corresponding (default) redistribution matrix file (RMF; see George et al. 1992a). - ANCRFILE="", # - the name of the corresponding (default) ancillary response file (ARF; see George et al. 1992a). - HDUCLASS="OGIP", # - should contain the string "OGIP" to indicate that this is an OGIP style file. - HDUCLAS1="SPECTRUM", # - should contain the string "SPECTRUM" to indicate this is a spectrum. - HDUVERS="1.2.1", #- the version number of the format (this document describes version 1.2.1) - #POISSERR #- whether Poissonian errors are appropriate to the data (see below). - #CHANTYPE #- whether the channels used in the file have been corrected in anyway (see below). - DETCHANS=len(self._rate), #- the total number of detector channels available. - )), - ) + return fits.BinTableHDU.from_columns( + [ + fits.Column(name="RATE", array=self._rate, format="1E"), + fits.Column(name="STAT_ERR", + array=self._stat_err, format="1E"), + fits.Column(name="SYS_ERR", array=self._sys_err, format="1E"), + fits.Column(name="QUALITY", array=self._quality, format="1I"), + ], + # https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html + header=fits.Header( + cards=dict( + # - the name (i.e. type) of the extension + EXTNAME="SPECTRUM", + # - the "telescope" (i.e. mission/satellite name). + TELESCOP="", + INSTRUME="", # - the instrument/detector. + # FILTER - the instrument filter in use (if any) + AREASCAL=1.0, + BACKSCAL=1.0, + # the integration time (in seconds) for the PHA data (assumed to be corrected for deadtime, data drop-outs etc. ) + EXPOSURE=self._exposure, + # - the name of the corresponding background file (if any) + BACKFILE="", + # CORRFILE - the name of the corresponding correction file (if any) + # CORRSCAL - the correction scaling factor. + # - the name of the corresponding (default) redistribution matrix file (RMF; see George et al. 1992a). + RESPFILE="", + # - the name of the corresponding (default) ancillary response file (ARF; see George et al. 1992a). + ANCRFILE="", + # - should contain the string "OGIP" to indicate that this is an OGIP style file. + HDUCLASS="OGIP", + # - should contain the string "SPECTRUM" to indicate this is a spectrum. + HDUCLAS1="SPECTRUM", + # - the version number of the format (this document describes version 1.2.1) + HDUVERS="1.2.1", + # POISSERR #- whether Poissonian errors are appropriate to the data (see below). # noqa: E800 + # CHANTYPE #- whether the channels used in the file have been corrected in anyway (see below). # noqa: E800 + DETCHANS=len( + self._rate + ), # - the total number of detector channels available. + ) + ), + ) def to_fits(self, fn: str): logger.info("store to fits %s %s", self, fn) - fits.HDUList([ + fits.HDUList( + [ fits.PrimaryHDU(), self.spectrum_hdu, - ]).writeto(fn, overwrite=True) + ] + ).writeto(fn, overwrite=True) @property def total_counts(self): @@ -209,173 +231,208 @@ class PHAII(Spectrum): class RMF: - _telescop="not-a-telescope" - _instrume="not-an-instrument" + _telescop = "not-a-telescope" + _instrume = "not-an-instrument" - _energ_lo = None # type: ignore - _energ_hi = None # type: ignore - _matrix = None # type: ignore - _e_min = None # type: ignore - _e_max = None # type: ignore + _energ_lo = None # type: ignore + _energ_hi = None # type: ignore + _matrix = None # type: ignore + _e_min = None # type: ignore + _e_max = None # type: ignore - def to_long_string(self): return f"{self.__class__.__name__}: {len(self._e_min)} x {len(self._energ_lo)} channels " @staticmethod def from_file_name(fn): return Readable.from_file_name(RMF, fn) - + @staticmethod def from_file_name_osaisgri(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f['ISGR-RMF.-RSP'].data['ENERG_LO'], - energ_hi=f['ISGR-RMF.-RSP'].data['ENERG_HI'], - matrix=np.stack(f['ISGR-RMF.-RSP'].data['MATRIX']), - e_min=f['ISGR-EBDS-MOD'].data['E_MIN'], - e_max=f['ISGR-EBDS-MOD'].data['E_MAX'], - ) - + energ_lo=f["ISGR-RMF.-RSP"].data["ENERG_LO"], + energ_hi=f["ISGR-RMF.-RSP"].data["ENERG_HI"], + matrix=np.stack(f["ISGR-RMF.-RSP"].data["MATRIX"]), + e_min=f["ISGR-EBDS-MOD"].data["E_MIN"], + e_max=f["ISGR-EBDS-MOD"].data["E_MAX"], + ) + @staticmethod def from_file_name_osaspi(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f['SPI.-RMF.-RSP'].data['ENERG_LO'], - energ_hi=f['SPI.-RMF.-RSP'].data['ENERG_HI'], - matrix=np.stack(f['SPI.-RMF.-RSP'].data['MATRIX']), - e_min=f['SPI.-EBDS-SET'].data['E_MIN'], - e_max=f['SPI.-EBDS-SET'].data['E_MAX'], - ) - + energ_lo=f["SPI.-RMF.-RSP"].data["ENERG_LO"], + energ_hi=f["SPI.-RMF.-RSP"].data["ENERG_HI"], + matrix=np.stack(f["SPI.-RMF.-RSP"].data["MATRIX"]), + e_min=f["SPI.-EBDS-SET"].data["E_MIN"], + e_max=f["SPI.-EBDS-SET"].data["E_MAX"], + ) + @staticmethod def from_file_name_normal(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f['MATRIX'].data['ENERG_LO'], - energ_hi=f['MATRIX'].data['ENERG_HI'], - matrix=np.vstack(f['MATRIX'].data['MATRIX']), - e_min=f['EBOUNDS'].data['E_MIN'], - e_max=f['EBOUNDS'].data['E_MAX'], - ) + energ_lo=f["MATRIX"].data["ENERG_LO"], + energ_hi=f["MATRIX"].data["ENERG_HI"], + matrix=np.vstack(f["MATRIX"].data["MATRIX"]), + e_min=f["EBOUNDS"].data["E_MIN"], + e_max=f["EBOUNDS"].data["E_MAX"], + ) @staticmethod def from_file_name_alt_normal(fn): f = fits.open(fn) return RMF.from_arrays( - energ_lo=f['SPECRESP MATRIX'].data['ENERG_LO'], - energ_hi=f['SPECRESP MATRIX'].data['ENERG_HI'], - matrix=np.vstack(f['SPECRESP MATRIX'].data['MATRIX']), - e_min=f['EBOUNDS'].data['E_MIN'], - e_max=f['EBOUNDS'].data['E_MAX'], - ) + energ_lo=f["SPECRESP MATRIX"].data["ENERG_LO"], + energ_hi=f["SPECRESP MATRIX"].data["ENERG_HI"], + matrix=np.vstack(f["SPECRESP MATRIX"].data["MATRIX"]), + e_min=f["EBOUNDS"].data["E_MIN"], + e_max=f["EBOUNDS"].data["E_MAX"], + ) @staticmethod def from_arrays(energ_lo, energ_hi, matrix, e_min, e_max): self = RMF() - if not (len(energ_lo) == len(energ_hi) == matrix.shape[0]): - raise Exception(f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {matrix.shape[0]}!") - + raise Exception( + f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {matrix.shape[0]}!" + ) + if not (len(e_min) == len(e_max) == matrix.shape[1]): - raise Exception(f"incompatible dimensions of channel energy, bounds {len(e_min)} {len(e_max)} but matrix {matrix.shape[1]}!") + raise Exception( + f"incompatible dimensions of channel energy, bounds {len(e_min)} {len(e_max)} but matrix {matrix.shape[1]}!" + ) - self._energ_lo = energ_lo # type: ignore - self._energ_hi = energ_hi # type: ignore - self._matrix = matrix # type: ignore - self._e_min = e_min # type: ignore - self._e_max = e_max # type: ignore + self._energ_lo = energ_lo # type: ignore + self._energ_hi = energ_hi # type: ignore + self._matrix = matrix # type: ignore + self._e_min = e_min # type: ignore + self._e_max = e_max # type: ignore return self @property def matrix_hdu(self): - return fits.BinTableHDU.from_columns([ - fits.Column(name='ENERG_LO', array=self._energ_lo, format='1E'), - fits.Column(name='ENERG_HI', array=self._energ_hi, format='1E'), - fits.Column(name='N_GRP', array=np.ones_like(self._energ_lo), format='1I'), - fits.Column(name='F_CHAN', array=0*np.ones_like(self._energ_lo), format='1I'), - fits.Column(name='N_CHAN', array=len(self._e_min)*np.ones_like(self._energ_lo), format='1I'), - fits.Column(name='MATRIX', array=self._matrix, format='PE'), - #fits.Column(name='MATRIX', array=self._matrix, format=f'{len(self._e_min)}E'), - ], - header=fits.Header(cards=dict( + return fits.BinTableHDU.from_columns( + [ + fits.Column(name="ENERG_LO", + array=self._energ_lo, format="1E"), + fits.Column(name="ENERG_HI", + array=self._energ_hi, format="1E"), + fits.Column( + name="N_GRP", array=np.ones_like(self._energ_lo), format="1I" + ), + fits.Column( + name="F_CHAN", array=0 * np.ones_like(self._energ_lo), format="1I" + ), + fits.Column( + name="N_CHAN", + array=len(self._e_min) * np.ones_like(self._energ_lo), + format="1I", + ), + fits.Column(name="MATRIX", array=self._matrix, format="PE"), + # fits.Column(name='MATRIX', array=self._matrix, format=f'{len(self._e_min)}E'), # noqa: E800 + ], + header=fits.Header( + cards=dict( # https://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/docs/memos/cal_gen_92_002/cal_gen_92_002.html#tth_sEc7.1.1 - EXTNAME = 'MATRIX', # / name of this binary table extension - #BITPIX = 8 / 8-bit bytes - #NAXIS = 2 / 2-dimensional binary table - #NAXIS1 = 34 / width of table in bytes - #NAXIS2 = 1180 / number of rows in table - #PCOUNT = 1031160 / Number of bytes acumulated in heap - #GCOUNT = 1 / one data group (required keyword) - #TFIELDS = 6 / number of fields in each row - #TTYPE1 = 'ENERG_LO' / label for field 1 - #TFORM1 = 'E ' / data format of field: 4-byte REAL - #TUNIT1 = 'keV ' / physical unit of field - #TTYPE2 = 'ENERG_HI' / label for field 2 - #TFORM2 = 'E ' / data format of field: 4-byte REAL - #TUNIT2 = 'keV ' / physical unit of field - #TTYPE3 = 'N_GRP ' / label for field 3 - #TFORM3 = 'I ' / data format of field: 2-byte INTEGER - #TTYPE4 = 'F_CHAN ' / label for field 4 - #TFORM4 = 'PI(2) ' / data format of field: variable length array - #TTYPE5 = 'N_CHAN ' / label for field 5 - #TFORM5 = 'PI(2) ' / data format of field: variable length array - #TTYPE6 = 'MATRIX ' / label for field 6 - #TFORM6 = 'PE(418) ' / data format of field: variable length array - TLMIN4=0, #/ First legal channel number - TLMAX4=len(self._e_min), #/ Highest legal channel number - TELESCOP=self._telescop, # / mission/satellite name - INSTRUME=self._instrume, # / instrument/detector - #FILTER = 'NONE ' / filter information - CHANTYPE='PI', # / Type of channels (PHA, PI etc) - DETCHANS=len(self._e_min), # / Total number of detector PHA channels - LO_THRES=1.00E-07, # / Lower probability density threshold for matrix - HDUCLASS='OGIP', #/KeywordinformationforCaltoolsSoftware. - HDUCLAS1='RESPONSE', #/KeywordinformationforCaltoolsSoftware. - HDUCLAS2='RSP_MATRIX', #/KeywordinformationforCaltoolsSoftware. - HDUVERS='1.3.0', #/KeywordinformationforCaltoolsSoftware. - HDUCLAS3='DETECTOR', #/KeywordinformationforCaltoolsSoftware. - )), - ) - + EXTNAME="MATRIX", # / name of this binary table extension + # BITPIX = 8 / 8-bit bytes + # NAXIS = 2 / 2-dimensional binary table + # NAXIS1 = 34 / width of table in bytes + # NAXIS2 = 1180 / number of rows in table + # PCOUNT = 1031160 / Number of bytes acumulated in heap + # GCOUNT = 1 / one data group (required keyword) + # TFIELDS = 6 / number of fields in each row + # TTYPE1 = 'ENERG_LO' / label for field 1 + # TFORM1 = 'E ' / data format of field: 4-byte REAL + # TUNIT1 = 'keV ' / physical unit of field + # TTYPE2 = 'ENERG_HI' / label for field 2 + # TFORM2 = 'E ' / data format of field: 4-byte REAL + # TUNIT2 = 'keV ' / physical unit of field + # TTYPE3 = 'N_GRP ' / label for field 3 + # TFORM3 = 'I ' / data format of field: 2-byte INTEGER + # TTYPE4 = 'F_CHAN ' / label for field 4 + # TFORM4 = 'PI(2) ' / data format of field: variable length array + # TTYPE5 = 'N_CHAN ' / label for field 5 + # TFORM5 = 'PI(2) ' / data format of field: variable length array + # TTYPE6 = 'MATRIX ' / label for field 6 + # TFORM6 = 'PE(418) ' / data format of field: variable length array + TLMIN4=0, # / First legal channel number + TLMAX4=len(self._e_min), # / Highest legal channel number + TELESCOP=self._telescop, # / mission/satellite name + INSTRUME=self._instrume, # / instrument/detector + # FILTER = 'NONE ' / filter information + CHANTYPE="PI", # / Type of channels (PHA, PI etc) + DETCHANS=len( + self._e_min + ), # / Total number of detector PHA channels + LO_THRES=1.00e-07, # / Lower probability density threshold for matrix + HDUCLASS="OGIP", # /KeywordinformationforCaltoolsSoftware. + # /KeywordinformationforCaltoolsSoftware. + HDUCLAS1="RESPONSE", + # /KeywordinformationforCaltoolsSoftware. + HDUCLAS2="RSP_MATRIX", + HDUVERS="1.3.0", # /KeywordinformationforCaltoolsSoftware. + # /KeywordinformationforCaltoolsSoftware. + HDUCLAS3="DETECTOR", + ) + ), + ) + @property def ebounds_hdu(self): - return fits.BinTableHDU.from_columns([ - fits.Column(name='CHANNEL', array=np.arange(self._e_min.shape[0]), format='1I'), - fits.Column(name='E_MIN', array=self._e_min, format='1E'), - fits.Column(name='E_MAX', array=self._e_max, format='1E'), - ], - header=fits.Header(cards=dict( - EXTNAME='EBOUNDS', # / name of this binary table extension - TLMIN1=0, #/ First legal channel number - TLMAX1=len(self._e_min), # 511 / Highest legal channel number - TELESCOP=self._telescop, # / mission/satellite name - INSTRUME=self._instrume, # / instrument/detector - #FILTER = 'NONE ' / filter information - CHANTYPE='PI', # / Type of channels (PHA, PI etc) - DETCHANS=len(self._e_min), # / Total number of detector PHA channels - #SMOOTHED= 0 / 0 = raw, 1-12 = smooth, -1 = ep-lin, -2 = mean- - HDUCLASS='OGIP', # / Keyword information for Caltools Software. - HDUCLAS1='RESPONSE', # / Keyword information for Caltools Software. - HDUCLAS2='EBOUNDS', # / Keyword information for Caltools Software. - HDUVERS ='1.2.0', # / Keyword information for Caltools Software. - )), - ) - + return fits.BinTableHDU.from_columns( + [ + fits.Column( + name="CHANNEL", array=np.arange(self._e_min.shape[0]), format="1I" + ), + fits.Column(name="E_MIN", array=self._e_min, format="1E"), + fits.Column(name="E_MAX", array=self._e_max, format="1E"), + ], + header=fits.Header( + cards=dict( + EXTNAME="EBOUNDS", # / name of this binary table extension + TLMIN1=0, # / First legal channel number + TLMAX1=len( + self._e_min + ), # 511 / Highest legal channel number + TELESCOP=self._telescop, # / mission/satellite name + INSTRUME=self._instrume, # / instrument/detector + # FILTER = 'NONE ' / filter information + CHANTYPE="PI", # / Type of channels (PHA, PI etc) + DETCHANS=len( + self._e_min + ), # / Total number of detector PHA channels + # SMOOTHED= 0 / 0 = raw, 1-12 = smooth, -1 = ep-lin, -2 = mean- + # / Keyword information for Caltools Software. + HDUCLASS="OGIP", + # / Keyword information for Caltools Software. + HDUCLAS1="RESPONSE", + # / Keyword information for Caltools Software. + HDUCLAS2="EBOUNDS", + # / Keyword information for Caltools Software. + HDUVERS="1.2.0", + ) + ), + ) + def to_fits(self, fn): logger.info("store to fits %s %s", self, fn) - fits.HDUList([ + fits.HDUList( + [ fits.PrimaryHDU(), self.ebounds_hdu, self.matrix_hdu, - ]).writeto(fn, overwrite=True) + ] + ).writeto(fn, overwrite=True) @property def d_e(self): @@ -387,96 +444,105 @@ def d_energ(self): @property def c_e(self): - return (self._e_max + self._e_min)/2 + return (self._e_max + self._e_min) / 2 @property def c_energ(self): - return (self._energ_hi + self._energ_lo)/2 - + return (self._energ_hi + self._energ_lo) / 2 def rebin(pha: PHAI, rmf: RMF, new_e_bins) -> (PHAI, RMF): - new_e_bins_assigned = [] + new_e_bins_assigned = [] if pha._rate.shape[0] != rmf._matrix[0].shape[0]: raise RuntimeError() for new_bin_req_e1, new_bin_req_e2 in zip(new_e_bins[:-1], new_e_bins[1:]): m = rmf._e_min >= new_bin_req_e1 - m &= rmf._e_min < new_bin_req_e2 # really e_min here - - new_e_bins_assigned.append(dict( - mask=m, - matrix_row=rmf._matrix[:, m].sum(1), - e_min=np.array(rmf._e_min)[m].min(), - e_max=np.array(rmf._e_max)[m].max(), - rate=pha._rate[m].sum(), - stat_err=np.sum(pha._stat_err[m]**2)**0.5, - sys_err=np.sum(pha._sys_err[m]**2)**0.5, - )) - - logger.info("for %s - %s new bin assigned: %s", new_bin_req_e1, new_bin_req_e2, {k: v for k, v in new_e_bins_assigned[-1].items() if len(str(v)) < 50}) + m &= rmf._e_min < new_bin_req_e2 # really e_min here + + new_e_bins_assigned.append( + dict( + mask=m, + matrix_row=rmf._matrix[:, m].sum(1), + e_min=np.array(rmf._e_min)[m].min(), + e_max=np.array(rmf._e_max)[m].max(), + rate=pha._rate[m].sum(), + stat_err=np.sum(pha._stat_err[m] ** 2) ** 0.5, + sys_err=np.sum(pha._sys_err[m] ** 2) ** 0.5, + ) + ) + + logger.info( + "for %s - %s new bin assigned: %s", + new_bin_req_e1, + new_bin_req_e2, + {k: v for k, + v in new_e_bins_assigned[-1].items() if len(str(v)) < 50}, + ) return ( PHAI.from_arrays( exposure=pha._exposure, - rate=np.array([r['rate'] for r in new_e_bins_assigned]), - stat_err=np.array([r['stat_err'] for r in new_e_bins_assigned]), - sys_err=np.array([r['stat_err'] for r in new_e_bins_assigned]), + rate=np.array([r["rate"] for r in new_e_bins_assigned]), + stat_err=np.array([r["stat_err"] for r in new_e_bins_assigned]), + sys_err=np.array([r["stat_err"] for r in new_e_bins_assigned]), ), - RMF.from_arrays( + RMF.from_arrays( energ_lo=rmf._energ_lo.copy(), energ_hi=rmf._energ_hi.copy(), - matrix=np.vstack([r['matrix_row'] for r in new_e_bins_assigned]).transpose(), - e_min=np.array([r['e_min'] for r in new_e_bins_assigned]), - e_max=np.array([r['e_max'] for r in new_e_bins_assigned]), + matrix=np.vstack( + [r["matrix_row"] for r in new_e_bins_assigned] + ).transpose(), + e_min=np.array([r["e_min"] for r in new_e_bins_assigned]), + e_max=np.array([r["e_max"] for r in new_e_bins_assigned]), ), ) - - - class ARF: _arf = None @staticmethod def from_file_name(fn): return Readable.from_file_name(ARF, fn) - + @staticmethod def from_file_name_normal(fn): f = fits.open(fn) return ARF.from_arrays( - energ_lo=f['SPECRESP'].data['ENERG_LO'], - energ_hi=f['SPECRESP'].data['ENERG_HI'], - arf=f['SPECRESP'].data['SPECRESP'], - ) + energ_lo=f["SPECRESP"].data["ENERG_LO"], + energ_hi=f["SPECRESP"].data["ENERG_HI"], + arf=f["SPECRESP"].data["SPECRESP"], + ) @staticmethod def from_file_name_osa(fn): f = fits.open(fn) return ARF.from_arrays( - energ_lo=f['ISGR-ARF.-RSP'].data['ENERG_LO'], - energ_hi=f['ISGR-ARF.-RSP'].data['ENERG_HI'], - arf=f['ISGR-ARF.-RSP'].data['SPECRESP'], - ) + energ_lo=f["ISGR-ARF.-RSP"].data["ENERG_LO"], + energ_hi=f["ISGR-ARF.-RSP"].data["ENERG_HI"], + arf=f["ISGR-ARF.-RSP"].data["SPECRESP"], + ) @staticmethod def from_arrays(energ_lo, energ_hi, arf): self = ARF() if not (len(energ_lo) == len(energ_hi) == len(arf)): - raise Exception(f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {len(arf)}!") - - self._energ_lo = energ_lo # type: ignore - self._energ_hi = energ_hi # type: ignore - self._arf = arf # type: ignore + raise Exception( + f"incompatible dimensions of mc energy, bounds {len(energ_lo)} {len(energ_hi)} but matrix {len(arf)}!" + ) + + self._energ_lo = energ_lo # type: ignore + self._energ_hi = energ_hi # type: ignore + self._arf = arf # type: ignore return self def to_long_string(self): - return f"{self.__class__.__name__}: {len(self._arf)} energies {np.max(self._arf)}" - + return ( + f"{self.__class__.__name__}: {len(self._arf)} energies {np.max(self._arf)}" + ) From d79628aa3a920e37a4f423b94c0666561cc6a8c0 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:13:45 +0200 Subject: [PATCH 03/11] update gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 1bb2e45..08767e5 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,5 @@ build .vscode dist *.egg-info +.idea +.ipynb_checkpoints From 7c3de76d7e403b4e5a39c31e2d8ee34b70a8ee36 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:14:05 +0200 Subject: [PATCH 04/11] Add example in tools.ipynb --- tools.ipynb | 107 +++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 76 insertions(+), 31 deletions(-) diff --git a/tools.ipynb b/tools.ipynb index 5ddf4fe..093ffd0 100644 --- a/tools.ipynb +++ b/tools.ipynb @@ -3,6 +3,22 @@ { "cell_type": "code", "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "ename": "Exception", + "evalue": "failed to read with any method, tried: {'from_file_name_osa': FileNotFoundError(2, 'No such file or directory'), 'from_file_name_normal': FileNotFoundError(2, 'No such file or directory')}", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mException\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[1], line 6\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mogip\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtools\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m plot, transform_rmf, crab_ph_cm2_s_kev, convolve\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpylab\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[0;32m----> 6\u001b[0m crab_pha \u001b[38;5;241m=\u001b[39m \u001b[43mogip\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mspec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPHAI\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_file_name\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcrab_pha.fits\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 7\u001b[0m crab_rmf \u001b[38;5;241m=\u001b[39m ogip\u001b[38;5;241m.\u001b[39mspec\u001b[38;5;241m.\u001b[39mRMF\u001b[38;5;241m.\u001b[39mfrom_file_name(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcrab_rmf.fits\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 8\u001b[0m crab_arf \u001b[38;5;241m=\u001b[39m ogip\u001b[38;5;241m.\u001b[39mspec\u001b[38;5;241m.\u001b[39mARF\u001b[38;5;241m.\u001b[39mfrom_file_name(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcrab_arf.fits\u001b[39m\u001b[38;5;124m\"\u001b[39m) \n", + "File \u001b[0;32m~/Soft/ogip/ogip/spec.py:62\u001b[0m, in \u001b[0;36mPHAI.from_file_name\u001b[0;34m(fn)\u001b[0m\n\u001b[1;32m 60\u001b[0m \u001b[38;5;129m@staticmethod\u001b[39m\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfrom_file_name\u001b[39m(fn):\n\u001b[0;32m---> 62\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mReadable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_file_name\u001b[49m\u001b[43m(\u001b[49m\u001b[43mPHAI\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfn\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Soft/ogip/ogip/spec.py:47\u001b[0m, in \u001b[0;36mReadable.from_file_name\u001b[0;34m(cls, fn)\u001b[0m\n\u001b[1;32m 42\u001b[0m logger\u001b[38;5;241m.\u001b[39mdebug(\n\u001b[1;32m 43\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfailed to read with \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mcls\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, n, m, e\n\u001b[1;32m 44\u001b[0m )\n\u001b[1;32m 45\u001b[0m methods_tried[n] \u001b[38;5;241m=\u001b[39m e\n\u001b[0;32m---> 47\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m(\n\u001b[1;32m 48\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfailed to read with any method, tried: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmethods_tried\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[0;31mException\u001b[0m: failed to read with any method, tried: {'from_file_name_osa': FileNotFoundError(2, 'No such file or directory'), 'from_file_name_normal': FileNotFoundError(2, 'No such file or directory')}" + ] + } + ], "source": [ "import numpy as np\n", "import ogip\n", @@ -52,46 +68,75 @@ " assert np.nanstd(np.abs((s - t_s)/s)) < 0.05\n", "\n", "t_rmf.to_fits(\"transformed_rmf.fits\")\n" - ], - "outputs": [ - { - "output_type": "display_data", - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+3ElEQVR4nO3de1xUdf748deHu6AikZqlX0VX8QZyE3JNJfNWmaWm5T5q1Syz1rbcXc1Wy/JX20UftZVaaanb1paZaWrtalakdlHByHsohYJmmgRxldvn98cw48wwAwMDzAHez8eDR845Z875MMu+z2fe5/N5f5TWGiGEEM2fl6cbIIQQonFIwBdCiBZCAr4QQrQQEvCFEKKFkIAvhBAthAR8IYRoIXw83YDqXH755bpbt26eboYQLikoKCAoKMjTzRCClJSUX7TW7e23Gzrgd+vWjeTkZE83QwiXJCUlkZiY6OlmCIFS6qSj7ZLSEUKIFsKQAV8pdZNSamVubq6nmyKEEM2GIQO+1nqL1npmcHCwp5sihBDNhiEDvhBCiPonAV8IIVoIQ4/S0RWanevSOH/yN8rLqlb19PZRtO/alvixYQQE+XqghUII0XQYOuDn/FzI2XTnD27LyzRn03PZ/GJqteeRG4MQQjRiwFdK9QEeBC4HPtVav1LTe0qLz3P20KtuXzswuC/lZUMtNwa5AQghWiK3Ar5SajUwFjinte5vtX0M8CLgDbyutX5Ga30UmKWU8gLeBGoM+PWhtPxncrN/pjD3iGWb/Q0A5CYghGj+3O3hrwWWYQrgACilvIHlwEggC9inlNqstT6ilBoH3Af825WT+wa054r+s9xq4G+ZO22CvaMbADi+CViTG4IQoqlzK+BrrXcqpbrZbY4HTmitfwBQSr0L3Awc0VpvBjYrpT4C/uPonEqpmcBMgM7te5CXl4fyglaXQfv+Ch9/RdlFzflDmqJs0BXVt1G1iyaoXbTldcn5fZQUfI/10o5lFefIzf6ZgpzDTs/jFxROzq8DOZ6aZTqvXZuEyM/PJykpydPNEMKphsjhXwVkWr3OAhKUUonABMAf+NjZm7XWK5VSPwE3BQSr2BlPD3d84OiaG1JcUMrerT/ajvJpMxywPaf9twB7peU/U5Z3jtLCNMu2wOC+eAcNJTvF9Fq+AQippSOMrtEe2mqtk4AkF4/dAmyJi4u7x51rBgT5MvS2XpbXDm8AQNsuQ2nbZajT87iSFrJOCUnwF0IYUUME/NNAF6vXnSu3uUwpdRNw0+9+97v6bFeVGwA4vwlYs78h1HQDkOAvhDAiZZ3LrtMJTDn8reZROkopHyANuA5ToN8H/EFr7TxB7kRcXJz2dHlkV24I1jeA0vKfAfD17giYgn/bLkO5okdwlZuNaF4kpSOMQimVorWOs9/u7rDMd4BE4HKlVBawSGv9hlJqNrAN07DM1bUN9g3Vw68LV9JC1t8ArIN/WXk2hblHaNtlKGfTc3nvH/tszi29fyFEY3K7h9+QjNDDr0l13wDOHnqVsvJsfLwvc/he6f03L9LDF0bRID38hmKkHn5NrL8B2Af/wOC+Tkf/VNf7l56/EKIhSA+/gexcl1ZtHaDqev/S82+apIcvjEJ6+I0sfmxYtQ97nfX+pecvhGgo0sP3gOp6/9Lzb7qkhy+Mokn18Ju76nr/te35S69fCOEqQwb85pDSqY6jCWBg6vmD41m/5p6/fblo8ySvvVt/lF6/EKJaktIxkOqGeDqq92NO/bhSUVS+CTQ8SekIo5CUThNQ256/s16/I/JNQAghi5g3AfFjw7iiRzDePrZlmAOD+zqd1GXNnP8HOH/ytwZpoxDC+AzZw2/uOfzactTzry7fb8/6G0B5ma5S4sFM0j5CNG+GDPj1VR65OatpnL89V1I/kvYRonkzZMAXNXOW77e34bnkaks8mNVU6A3kG4AQTZ0E/Gaufde2lJfVnPpx5QGwfAMQommTh7bNnLMHvvZqegAsD36FaPoM2cOXh7b1x5XUjysPgM8eehUqKij58UcA3rrnB5v9viHBdIoNk3SPEAZmyIAvD20blysPgJWPD6UXz5FdsK3qzooKAkp6oNIvk3SPEAZmyIAvGldN3wI2PJdMUPtoCn896HB/ScEZii+mU/Ljj2T9WMHJrU/a7G87diwht02u1zYLIWpPAr6okenBbyytQ2Md7v/52OvosjIAKuweCxUfOwYgAV8IA5CAL2pUU8pH+figfHzwCwsDYA8PWPaVBPyIFxWcXJcm+X0hPKzRAr5S6hbgRqAt8IbWentjXVu4p7qUj3m8fknRWc6dWFtlf0VxMQG+YZxN7yH5fSE8zK2Ar5RaDYwFzmmt+1ttHwO8CHgDr2utn9FabwI2KaVCgKWABPxmwNtHERgS4XR/WXk2xZWje+zz+5LbF6JxudvDXwssA940b1BKeQPLgZFAFrBPKbVZa22e6rmwcr9oBuqa33eU26+uPDTITF8h3OVWwNda71RKdbPbHA+c0Fr/AKCUehe4WSl1FHgG+K/Wer871xXGUdf8vqPc/t6tP3Ji72dORwMFhkRQXhbL5hdTLdvkJiCE6xoih38VkGn1OgtIAB4ARgDBSqnfaa0dzuFXSs0EZgJ07NiRpKSkBmiiqFcdIbRj1c3p2yooKy+nvOQ8Z75/w2afKikhwKcbx1M7cDw1y7I975dUykvO4+3X3ub48uIsLhac5LdfUm22+7XuTc6vkZzOyqJTrGcnjufn58vfqzC0Rntoq7V+CXjJheNWKqV+Am5q06ZNrKwg1HRd2JuMKoly2GO/WHGO/JJz+P/8s832ipLz+Le6gg6/m2azPf9CSpXzlBSdpaLImzZXDca7XJGYWGWBn0YlK14Jo2uIgH8a6GL1unPlNpfJTNvmobr8fs732yi+mE5FcbHNdh/VjgD/7lWObx1a9TzWo4Ks6/xLmkcIxxoi4O8DeiqlwjAF+tuBP9TmBFJLp3moLr/fpsNAAnOqDtGsKCyA36Dkxx8JqThH/5K9lOBHml8UuV6X2zz4rSgupqzi1yrDQc25fhkGKoQtd4dlvgMkApcrpbKARVrrN5RSs4FtmIZlrtZaH67NeaWH3zw4G79f3Wic8l+z0Tm/Elxxjl4lqQD4UUL/kr02x30ZcAP++krQ5TbfEsrKs9FlZbQOjZWqnkLYcXeUzhQn2z8GPq7reaWH37y5unhLdU6uS+N0cgjlObk227Pz/msZBippHiFsGbK0gvTwRU3ix4axF6p+Szjk5XAhF/PiLeYhnXIDEC2RIQO+9PBFTRx9S3C2nGNp+c/kZv9ss726GwBQY7louWGIpkhpXfMC2J4SFxenk5OTPd0M0UTsXJfG2fTcKtt/y9xpE+xLy01DQX29L00eCAzua7MAjP177NkfD5CXl0e7kLZyIxAep5RK0VpXGadsyIBv1cO/5/jx455ujmgiairNYObKDcDRtpr2aa0JatePtl2GckWPYBkhJDymSQV8M+nhC3e5chNw1pt31Iuv7viy8mx8vC/jiv6zbLZL+kc0Ngn4QuD6t4C6+OngK5RX/FplMXjzjUN6/aKxNKmALykd0ViquwE46plXd/yFHz6jtDDNZpujXr/0+EVDa1IB30x6+KKp2LkujeOpWbRp08Zm+9lDr1qCvpn0+EVDcxbwDTksU4imJn5sGKezsvAuVzY9f/thomXl2RTmHqFtl6GcTc+ViWGiUUnAF6IeBAT50inWy6Zi5851acBQmwe/5h6/s4lhUv9HNCTPFhB3Qil1k1JqZW5u1THVQjQV8WPDuKJHMN4+yrItMLhvlYe65l4/YOn1b3gumZ3r0iguKG3UNovmTXL4QtSTmurhO5sYJnl+Ud+c5fAN2cMXojly1OOHqr1+6fGLhmL4HP5z78wkOa+yNG5Fuem/Xt4A/ObVjl+9Q22OH3rlKJaMvrcxmyiESxzV/3GW57dXXqY5m54rOX7hFkMH/KO/pPPvkovgD30uelsCPUBQRQFBFQW0Lbtg2uDlzVH/cv53No2Try6zbAMYXR7IjICQSyeOuBXipjfWryGEU84WibF/sGtO8UiNf+EOQwZ888SrkP9rRXRxKcMqLmPGfV/ZHpS8Bg6+b7PpjcJf2eZdaJOoOupfzlHy2FaUB0CAKmXYvlXMkIAvDMC+179zXRq/ZTofymmu8S/DOEVdGPuhbbdgnbzoGrd65HO3vcbOM9strwu9TDMh+1z0rpIisrl2m3jmTVlZp2uKlqk+FjF3NJNXHuqK2pKZtpVsngk4cdTfdCPoc/HSjUDSQqIm9RHw7W14Lplff/yiSo9fyjWI6kjArwX7m4LlBlBkeh2EaQ3VAq8gyzHyjUA0RMB3NJRTevyiJhLw3WCfFgopv0DbihxLSuhoK9N28zcCCf4tU0MEfEcpHvvyzNLjF/Y8HvCVUt2BBUCw1vpWV95jlIBfk+femUly7tcAZPlVEF5Swhr/yl6WpH5ajIYI+Pakxy9c0SATr5RSq5VS55RSh+y2j1FKfa+UOqGUmg+gtf5Baz3DnesZ1bwpK3lv1kHem3WQthXtOewXwOTS0/yxJJ039q3ydPNEM+JKuQb7iVsyaUuYuTsscy2wDHjTvEEp5Q0sB0YCWcA+pdRmrbXzBUKbkYguU9h5ZjsnfU0jgr4lj20roxweK6kfUVuOhnHWVKBNCrMJM7cCvtZ6p1Kqm93meOCE1voHAKXUu8DNQIsI+KZZvqaZvtWNCMr0KYUaRgsJURNHE7esSzI7Kscs+f2WqyEmXl0FZFq9zgISlFKhwFNAtFLqEa31047erJSaCcwE6NixI0lJSQ3QxMYR3+kPxHf6g8N9y79/EI1u0r+fsJWfn++Z/z07Qmjleuo/pVSg2kUT1C7a1KaMtygrz+ang69YDvcLCifn14GczsqiU6yU02pJGm2mrdb6AjDLheNWKqV+Am5q06ZNbEM/BPOUFWmKTJ9SVqQ9VGUCmKR6mqbGeGhbk+KBtqN6dLt+VUb0lBamEdpmOFyEC3ult9+SNETAPw10sXrduXKby7TWW4AtcXFx99Rnw4wkrk38pZSO1UxfSfUId1jn+CW/L+y5PSyzMoe/VWvdv/K1D5AGXIcp0O8D/qC1PlyLc7bYRcwnr4wi06eULmVVe1vS8zc2I/TwrdU0ht9+/L7k9puPBlnTVin1DpAIXK6UygIWaa3fUErNBrYB3sDq2gT7ls6m529Fev6itmoa0SO9/ZZHZto2Ec56/tLrNw6j9fDt2ff4pbfffDVID7+hWKV0PN0Uw3DU85dev6iN6vL70ttvGQwZ8FvCQ9vactSLN/f6JzuZ2GVNvgkIa/bj952N3ZcFV5oXQwZ84Rpn+X578k1A2Kupt29mXnDFTFI9TZshA76kdFzjao/dlW8AouVyNFvXfolFM0n1NG2GDPiS0ql/NqkfJyt9SdqnZbLu7W94LtkmvWNNUj1NnyEDvqhfVVI/DpZ0lLSPAGjftS3lZbaTtcycpXokzdN0GDLgS0qnfrnSa3flAbB8A2j+HKV3rDlK9Uiap+kwZMCXlE7jq+kBsHwDaBnsJ2uZOUv1SJqnaTFkwBeNr6aee03fAKT337w5S/VImqdpMWTAl5SO8VT3DeCofzlHS74mWW4GzVZ1qR5J8zQdhgz4ktIxnuoCtiz00vw5SvVImqfpMWTAF01LdTeD6lJB0vNv2lxN8wjjkIAvGpSzVFCWTylBuZ/DmhurviniVoib3gitE+6oaUSPmeT1jUMCvmhQznrwY15P5LDfBSaXnraZCBagixm2bxUzJOAbnqM0jzm4S17fmCTgC4+I6DKFnWe2c9Kuo6d1GpRmM0N6/k2St4+SvL6BGTLgyyid5m/J6HuBe6tsH/N6Isd8peffVMnwTWMz5JL1WustWuuZwcHBnm6KaGQRXaagVC9O+vbgpH8v049vD475+rLNu9DTzRM1iB8bxhU9gvH2US4dX16mOZuey96tPzZwywQYtIcvWi5nPX9no31kpI+xOBu+Cc7z+pLmaTyG7OELYS+uTbxpeceKckuqJ9On1On4f2Ec7bu2JTC4Lz7el9lsN+f1QYZvNhbp4YsmwdmKX8L44seGATdx/uQwm8Bu39t/7x/7JKffwBot4CulgoAVQAmQpLV+u7GuLZovSfMYX22Gb8rQzYblVkpHKbVaKXVOKXXIbvsYpdT3SqkTSqn5lZsnAO9rre8BxrlzXSHAKs0DllSPpHmaBvPwTes0j3WK52x6Lu/9Yx8bnktm57o0igtKPdXUZsXdHv5aYBnwpnmDUsobWA6MBLKAfUqpzUBn4GDlYeVuXlcISfM0YY6Gb5499Kr0+BuYWwFfa71TKdXNbnM8cEJr/QOAUupd4GZMwb8zkIo8LBYNKMOnnIQ1E222Db1yVOUIIGEEjsoy2E/Ysp6sZe7xS47fPQ2Rw78KyLR6nQUkAC8By5RSNwJbnL1ZKTUTmAnQsWNHkpKSGqCJorkaXuTDV/7FFJMGgFZenPStICPjBDkvvGlzbOv8HwDIb92dnzsO5acrRzs8Z6cz2+j4806bbY6Oz8/Pl7/X2ugIoR1N//wppQLVLpqgdtGW3fkZb1FWns1PB1+xbPMLCifn14GczsqiU6z0G2ur0R7aaq0LgBqnSWqtVyqlfgJuatOmTWxiYmKDt000H4mt/8Ssg+/bbJtcepoM33Jmtqoc6+3lzejyQGYUm/782xVn0u7iAcITn3Z80jVLoDgTMI8wUQ6PT0pKQv5e66Z4YGmVHr9u169Kj7+0MI3QNsPhIpz76tL7pefvmoYI+KeBLlavO1duc5nUwxd1Fje9Sr2drtte4+SZ7Zz0N70u9ErjKHmsDIg37S9NZ3Txr8ywes/cba+x88z2yv2nIfQyTvr2cHq8cI/9SJ6d69KAqjl+Z6xn7Equ3zmltXsTHipz+Fu11v0rX/sAacB1mAL9PuAPWuvDtTinuZbOPcePH3erfULYsw7mYLoBAPS56G2Z1HW0lWlfnyIIoIRir1aWgK91Gr1LS3nTrwecPWA68IpIcnJyaDfkbinwVg+KC6r2+M0Pde0ncMGlGbvePoqJ8+Iau7mGo5RK0VpX+SDc6uErpd4BEoHLlVJZwCKt9RtKqdnANsAbWF2bYC9EQ7Mv3+Boxa4+FyvTPjmVj6NGzrcEcpsCb6Ghpv2lpwloJQXe6oujHv9vmVWrcILtw12ZsVs9d0fpTHGy/WPgYzfOKykd0WhqO0mrutLOxWWFkuppAM5m60L1qR5hy5ClFaQ8sjAyZwXeJr02gEzfUia/GmHa4OUNyMzf+uBoti5cmrErXGPIcU1SHlk0Rb3pc2nmbyWZ+SuMxJA9fCGaouvCZ/H/7IZlWso6W/X6pcffsOx7/TJk8xJDBnxJ6Yjm4tIi7qbRP04Xb5flG93iqNa+mZRnuMSQAV8e2ormwr4n72jx9gBVKqN76sjZGrpmspauLUMGfCGaK0cjfGR0T905W0PXzH4t3ZbOkAFfUjqiuXI0wmfyyqgqBd+k2JtrHBVhE84ZMuBLSke0JKPLA/lCZ1OsTTN+M/zgYOYFHA39FLacDdcEGbLpiCEDvhAtyYyB9zDj4PuWMg3TQwKA0/B0ZUmqKyJND3XNReGmf+SZhoomz5ABX1I6okWxK/j20+uJZHtdYHJIAABBF9O44YsFTCoxTeRizY0yqqeOXO31N9ehnIYM+JLSES2Z5cFuZXXPUv0DBV5BTLpwwbThbOXCcRLwXVLdkE1nmutQTkMG/OqUlpaSlZVFcXGxp5siWrCAgAA6d+6Mr2/99/7sH+wmrJnIUZ9TJFxxqZzzhJyj3G6f8pEbgI2ahmw605yHcja5gJ+VlUWbNm3o1q0bSilPN0e0QFprLly4QFZWFmFhYQ1+vaFXjrIp53zCR/FBUFtuL5Yef3VqGrLpTHMeytnkAn5xcbEEe+FRSilCQ0M5f/58o1zPaY//sn4AdFVnZUEWB2TIZlVNLuADEuyFx3nyb9DS4/czvc7Q5WxDJm7Zq27IpjPNfSinIQO+jNIRwjn7Hr9M3BKuMmTAb86jdNauXUtycjLLli3zyPUff/xxWrduzd/+9je3jhHGMbo8kG0UclIXAFDoe4GdGRthzWYZs18P6rPX7+nhnoYM+E1dWVkZPj7y0YrGMSMghBmEWF6PKQ+irT4NZ3MujdmXSVu1UpehnK7w9HBPQy6AYnRvvvkmkZGRDBgwgDvvvBOAadOmMWvWLBISEpg3bx579+5l0KBBREdH8/vf/57vv//e8v7MzEwSExPp2bMnTzzxRI3XS0xMZM6cOcTFxdGnTx/27dvHhAkT6NmzJwsXLrQc9/zzz9O/f3/69+/PP//5T8v2p556il69enHNNdfYtCM9PZ0xY8YQGxvLkCFDOHbsWD18OsIjzh40zdQ9e4BOZafJ8qtgckg7pl9MY/0XCyz7WHMjJK/xdGsNyzyU09FC6e4yD/cEPDbcs0l3Q5/YcpgjZ+r3g+t7ZVsW3dTP6f7Dhw/z5JNP8tVXX3H55ZeTnZ1t2ZeVlcVXX32Ft7c3v/32G7t27cLHx4cdO3bw97//nQ0bNgCwd+9eDh06RGBgIAMHDuTGG28kLi6OG264gddff50rr7yyynX9/PxITk7mxRdf5OabbyYlJYXLLruMHj16MGfOHDIyMlizZg179uxBa01CQgLDhg2joqKCd999l9TUVMrKyoiJiSE2NhaAmTNn8uqrr9KzZ0/27NnD/fffz2effVavn6doBBG3mv5bWZrhGt2en0rgpGpLqe/PMmmrFuo6lNMVRhju2WgBXynVHVgABGutb22s69a3zz77jEmTJnH55ZcDcNlll3oCkyZNwtvbNP09NzeXqVOncvz4cZRSlJaWWo4bOXIkoaGhAEyYMIHdu3cTFxfHxx87X/d93LhxAERERNCvXz86deoEQPfu3cnMzGT37t2MHz+eoKAgy3l37dpFRUUF48ePJzAw0OY8+fn5fPXVV0yaNMlyjYsXL7r34QjPMJdmqFxUZcb0j5hR+e8Eulc/hNO8EIukeoDmP5TTpYCvlFoNjAXOaa37W20fA7wIeAOva62fcXYOrfUPwAyl1PvuNfmS6nrinmAOtgCPPvoo1157LRs3biQjI4NEq6Xv7If0uTLEz9/fNM/ey8vL8m/z67Kyslq3taKignbt2pGamlrr9wqDsg7alf8euu21yklbpge6GT4yhLM6dRnK6QqjDPd0NYe/FhhjvUEp5Q0sB64H+gJTlFJ9lVIRSqmtdj8d6rXVHjR8+HDWr1/PhcqvyNYpHWu5ublcddVVgGlkjrVPPvmE7OxsioqK2LRpE4MHD3a7XUOGDGHTpk0UFhZSUFDAxo0bGTJkCEOHDmXTpk0UFRWRl5fHli1bAGjbti1hYWGsX78eMM0e/e6779xuhzCWJaPvZc/0Dey593/sufd/dCvzNg3hfG0MY15P5I1ff5DcfgviUsDXWu8E7CNbPHBCa/2D1roEeBe4WWt9UGs91u7nXD2322P69evHggULGDZsGAMGDOAvf/mLw+PmzZvHI488QnR0dJUeeHx8PBMnTiQyMpKJEycSFxcHwA033MCZM2fq1K6YmBimTZtGfHw8CQkJ3H333URHRxMTE8Ntt93GgAEDuP766xk4cKDlPW+//TZvvPEGAwYMoF+/fnz44Yd1urZoOkaXB9KtxPTvbK8LfBFYYXpx9uClkTyi2VJau5anUkp1A7aaUzpKqVuBMVrruytf3wkkaK1nO3l/KPAUMBJT+udpJ8fNBGYCdOzYMfbdd9+12R8cHIxMyBJGcOLECXJzcy2v8/Pzad26tQdbVLOobxfQOv8H8lt3557AfE76VkDJlXRVZ7mu2JfeA5Z6uonNUvq2CvIz3gKgdbc7AOgxuuEGSV577bUpWus4++2N9tBWa30BmOXCcSuVUj8BN7Vp0ybWOvcNcPToUdq0adMwjRSiFgICAoiOjra8TkpKwv7v1XB+bAe/XaRd8Skm4ssHrYI4qRQnfSv4TJUxy5X2y4PeWjv31T4KKp/VmeNXYuLA6t7SINy5xZwGuli97ly5TQhhVBG3gp9pcMHt7frwXokfe/y8L+X210wkYc1E5m57zcMNFQ3BnYC/D+iplApTSvkBtwOb66NRWustWuuZwcHB9XE6IYRZ3HRT/fwrIk099D+nwPSPLuX2Swoo1KdsyjFbJK8x9e6tJ3FZP+w1vxaG5VLAV0q9A3wNhCulspRSM7TWZcBsYBtwFHhPa324PhqllLpJKbXSOj8qhGg4MwJCeO/XHPb4eROo/u/SDusgfvD9SxO3zORhb5PiUg5faz3FyfaPAeezheqoORdPE6IpKNSnTNU3S8pNG9ZMpGvpaUa368yMgMq6PdM/quzxH7zU8wdZc9fADFlLR3r4QnjO0CtH2fbyK2X4lLPNu9B2Y8StcEWE7Tbp9RuWIQN+c87hr127ltmzHY5cdUlGRgb9+/ev8Zj//Oc/db6GaOamf+R4hE1JAZw9wJIzm9lDsekn+7Dph2J6W5UHsYibful85mcD9jcAYRiGDPhNXV1KHdQnCfii1qxG7zhTrAL4zaud8xuGMDxDVss0+opXb775JkuXLkUpRWRkJP/+97+ZNm0aAQEBfPvttwwePJjbb7+dBx98kOLiYlq1asWaNWsIDw8HLpVHPn36NHfccQeLFi2q9nopKSncddddAIwaNcqyPSMjgzvvvJOCAlOdlGXLlvH73/+e+fPnc/ToUaKiopg6dSrjx493eJwQFnHTHdfMtxpzf3LNxEu5fSuyulbTYciA7/JD2//OrzpqwF1XRMD1TmvAeaQ88vTp01m2bBlDhw5l7ty5lu0dOnTgk08+ISAggOPHjzNlyhSSk5N55plnWLp0KVu3bgWgsLDQ4XFC1IZlLV0rl4Zw3nvpRlHd0EyZtOVRhgz4Ru7hN3Z55JycHHJychg61FSf+8477+S///0vAKWlpcyePZvU1FS8vb1JS0tz2GZXjxOiOvZr6QJVevvC2AwZ8F3u4VfTE/eEhiyP7MgLL7xAx44d+e6776ioqCAgIMCt44SoC/s0T9fS04wuD5QSzDXwRMlkQwZ8Ixs+fDjjx4/nL3/5C6GhoWRnZ9v08s1cKY/cqlUrNm3axOrVq51er127drRr147du3dzzTXX8Pbbb9tco3Pnznh5efGvf/2L8nLTmOk2bdqQl5dX43FC1KiG1IujNI/U3HeuodbKdZWM0qklT5RHXrNmDX/605+IiorCurrp/fffz7/+9S8GDBjAsWPHLN8wIiMj8fb2ZsCAAbzwwgtOjxPCXZZ6+1Y/3cq8qx5YXVmGFlCSoSHXyq0Nl8sjNyarHP49x48ft9l39OhR+vTp45mGCWHF/m+xSVTLrE49PVCdvDIKgPdmptqe++xBoDLeXBFpem09Zr8ZP8jduS6Ns+mNN5H0tgXxni2PXBtSWkEID6jHgGuuvGnmtCyDK2pzIzLoKCCjrJVryIAvhGi6RpcHso1CTvpe2nYprx/iuYZ5UEOtlevMbQscb5eAL4SoVzMCQkyBffoGyzZzmkd4lgR8IYSx1DYtk7zm0ixh64qdZlK508KQAd/IE6+EEHWT4VNOgl255QkFv3F7ceWkRHNZ5doy1+l3VLTNPBNfAj5g0GGZzblaphAt0ejywCrDNU/4KD4Iantpgztlla+IsK3YaangKZU7rRky4BtZTk4OK1asaJRrTZkyhcjISF544YVGuZ49d0s5u2LatGm8/77naqcnJibWWFfIlWNE9WYEhPCe71Xs8fM2/UzfgK/qzknfHnUrq9wCxu43BEOmdIzMHPDvv//+KvvKysrw8amfj/Ts2bPs27ePEydOuPye+rx+fSovL7fUGBItjH1Qrs1wSXNuXlbSqjfSw6+l+fPnk56eTlRUFHPnziUpKYkhQ4Ywbtw4+vbtC8Att9xCbGws/fr1Y+XKlZb3tm7dmgULFjBgwACuvvpqfv75ZwDWr19P//79GTBggKVI2qhRozh9+jRRUVHs2rWL1NRUrr76aiIjIxk/fjy//vorYOp9PvTQQ8TFxfHiiy+SmJjInDlziIuLo0+fPuzbt48JEybQs2dPFi5caGnLW2+9RXx8PFFRUdx7772Wcgtr1qyhV69exMfH8+WXX9b4eZSXl/O3v/2N/v37ExkZycsvvwxAt27dePjhh4mJiWH9+vWsWrWKgQMHMmDAACZOnEhh4aWVk3bs2EFcXBy9evWyVPisTuvWrZk7dy79+vVjxIgR7N27l8TERLp3787mzZsBKC4uZvr06URERBAdHc3nn38OQFFREbfffjt9+vRh/PjxFBUVWc67fft2Bg0aRExMDJMmTSI/P7/Gtgj3FOpTJJSYcvuTS0/zRvGvl3bar6ErK2m5zXjdwVp4du+zHMs+Vq/n7H1Zbx6Of9jp/meeeYZDhw6RmpoKmGZX7t+/n0OHDhEWFgbA6tWrueyyyygqKmLgwIFMnDiR0NBQCgoKuPrqq3nqqaeYN28eq1atYuHChSxevJht27Zx1VVXkZOTA8DmzZsZO3as5TrmYDps2DAee+wxnnjiCf75z38CUFJSYkk5bNmyBT8/P5KTk3nxxRe5+eabSUlJ4bLLLqNHjx7MmTOHc+fOsW7dOr788kt8fX25//77efvttxk5ciSLFi0iJSWF4OBgrr32WqKjo6v9vFauXElGRgapqan4+PjYlIsODQ1l//79AFy4cIF77jHNo1u4cCFvvPEGDzzwAGCq6793717S09O59tprOXHiBNnZ2dx9990OK4gWFBQwfPhwlixZwvjx41m4cCGffPIJR44cYerUqYwbN47ly5ejlOLgwYMcO3aMUaNGkZaWxiuvvEJgYCBHjx7lwIEDxMTEAPDLL7/w5JNPsmPHDoKCgnj22Wd5/vnneeyxx6r9/UXdXarDY1qnweFYfcnB16tGC/hKqVuAG4G2wBta6+3Vv6PpiI+PtwR7gJdeeomNGzcCpsVOjh8/TmhoKH5+fowdOxaA2NhYPvnkEwAGDx7MtGnTmDx5MhMmTKhy/tzcXHJychg2bBgAU6dOZdKkSZb9t912m83x48aNAyAiIoJ+/frRqVMnALp3705mZia7d+8mJSWFgQMHAqZeb4cOHdizZw+JiYm0b9/ect6aSinv2LGDWbNmWVJJ1oXkrNt16NAhFi5cSE5ODvn5+YwePdqyb/LkyXh5edGzZ0+6d+/OsWPHiIqKchjsAfz8/BgzZozld/T398fX15eIiAgyMjIA2L17t+WG0rt3b7p27UpaWho7d+7kz3/+M2C6iUZGRgLwzTffcOTIEQYPHgyYbqKDBg2q9ncX7rGUW65M+0wu9bYZydNVnTVV3QyoxWQtg82wNRqXAr5SajUwFjinte5vtX0M8CLgDbyutXZar1hrvQnYpJQKAZYCbgf86nrijcm6GFlSUhI7duzg66+/JjAwkMTERIqLiwHw9fW1lEL29va2FFV79dVX2bNnDx999BGxsbGkpKTU+foA/v7+AHh5eVn+bX5dVlaG1pqpU6fy9NNP27xv06ZNtbpubdo1bdo0Nm3axIABA1i7di1JSUmWfbUtF239OVr/jubfry601owcOZJ33nmnTu8X7rPM0K183eCzcw1ahqEhuZrDXwuMsd6glPIGlgPXA32BKUqpvkqpCKXUVrufDlZvXVj5vibJvvSwvdzcXEJCQggMDOTYsWN88803NZ4zPT2dhIQEFi9eTPv27cnMzLTZHxwcTEhICLt27QLg3//+t6W3XxfXXXcd77//PufOnQMgOzubkydPkpCQwBdffMGFCxcoLS1l/fr1lvds3LiRRx55pMq5Ro4cyWuvvWYJtNYpHWt5eXl06tSJ0tJSmxLPYHqGUVFRQXp6Oj/88INlKUh3DBkyxHKdtLQ0Tp06RXh4OEOHDrWs93vo0CEOHDA9ELz66qv58ssvLQ/JCwoKZKGYRmY/kqdb2aUev8Mcv6g1l3r4WuudSqludpvjgRNa6x8AlFLvAjdrrZ/G9G3AhjJ1yZ4B/qu13u/sWkqpmcBMgI4dO9r0BMEU/KoLuA3Nz8+P+Ph4+vbty8iRIxk9ejRlZWWWNg0ePJhly5YRHh5Oz549GThwIIWFhZb95v8WFRVRWlpKXl4ec+bMIT09Ha01w4YNo3v37pw6dYqKigrL8StWrOChhx6iqKiIbt26sWLFCvLy8igvL6egoMBynPXrwsJCm7aZ94WHh7NgwQJGjBhBRUUFvr6+LF26lPj4eObPn09CQgLBwcFERkZSUlJCXl4eR44cwd/fv8pnf9ttt3Ho0CH69++Pr68vU6dO5d5770VrTX5+vqX3vWDBAuLj4wkNDSUuLo78/Hzy8vIoLS2lU6dOxMXF8dtvv/H8889TWlrKqVOnmD17tmVZSHvmdly8eBFfX1+bduXl5XHnnXcyZ84c+vXrh4+PDytWrKCkpIQ77riD++67j/DwcMLDw4mKiqKgoICAgABWrFjB5MmTKSkpAUyL2HTq1KnKZ2xWXFxs8/eZn59f5e+1JYqqfA6VmpRk+bdZqoPPx/qY1vk/Yq6oeUuQPxuDgjlZWdE3w7ec/+kCelidN9Xu39W1pTb7miuXyyNXBvyt5pSOUupWYIzW+u7K13cCCVprhwO3lVJ/BqYC+4BUrXWNqwDExcVp+/HPUh7ZM+644w5eeOEFS35fNMPyyPXFOlXiyrBM8zERt9oOw7wi8tI2YHLpaaCy7LKja1R37trua+KUUp4tj6y1fgl4yZVjpbSC8bz11lueboJo7uKmm37sA7HVUExL2WW7Eg2ypKJr3An4p4EuVq87V24TQoh6N7o8kC90NsUqHXRlwC/1JtOnVJZUdJE7AX8f0FMpFYYp0N8O/KE+GiULoAgh7M0YeA8z7KtiXhHJ5NLT/ObVzvkbzx50XIbBUWVNs2Y6o9elUTpKqXeAr4FwpVSWUmqG1roMmA1sA44C72mtD9dHo5RSNymlVubmNt6SYEIIg4ubblUU7VKRtJO+PfjVO9TxeyJurf3krWY8o9fVUTpTnGz/GHA8O8YN0sMXQtQL83MBR5w9tG3GRdkMWUtHevhCCFH/DBnwjV4Pv3Xr1gCcOXOGW2+tw4INjSgpKclSzsGdY4RoMOY0jWhwhiye1lSGZV555ZUereUuhDAp1KdMwzUdGHrlqMq6PUJ6+G7IyMigf39TaaG1a9cyYcIExowZQ8+ePZk3b57luNqW3X388ceZOnUqQ4YMoWvXrnzwwQfMmzePiIgIxowZQ2mpaUm4Tz/9lOjoaCIiIrjrrru4ePEiAP/73//o3bs3MTExfPDBB5bzFhQUcNdddxEfH090dDQffvhhfX8kQjS6oVeOIlD9n8N9hfpUZUVOAQbt4bvq7D/+wcWj9Vse2b9Pb674+9/r9N7U1FS+/fZb/P39CQ8P54EHHqBVq1ZOy+4+9thjxMXFWapbWktPT+fzzz/nyJEjDBo0iA0bNvDcc88xfvx4PvroI8aMGcO0adP49NNP6dWrF3/84x955ZVXmDVrFvfccw+fffYZv/vd72wqVj711FMMHz6c1atXk5OTQ3x8PCNGjKjzZyWEEViqbjrgrNffUhky4DeVlI696667DvO3kr59+3Ly5ElycnKclt1dvHix03Ndf/31lpK/5eXlNuWAMzIy+P777wkLC6NXr16AqWTy8uXLSUxMJCwsjJ49ewKmkgjmRVi2b9/O5s2bWbp0KWCqBXPq1KkG+CSEEEZkyIDv6rDMuvbEG4p1KWJz+eO6lt21LvlrXw7YnRLAGzZsqFKN0rzylhCieTNkDr85aaiyu+Hh4WRkZFjOay6Z3Lt3bzIyMkhPTwewudGMHj2al19+GXPBvG+//dbtdgghmg4J+A2sffv2rF27lilTphAZGcmgQYM4dsz03OGxxx6zrMFaWwEBAaxZs4ZJkyYRERGBl5cXs2bNIiAggJUrV3LjjTcSExNDhw6XliJ49NFHKS0tJTIykn79+vHoo4/Wy+8ohGgaXC6P3Jiscvj3HD9+3GaflEcWRiHlkZ2obXnk6t7vzjGYHtoW6lOOR/GUmNbSxc92xbiupemm6pv3fVVzWw3KWXlkQ/bwm8qwTCGEsVU3ZNOZDJ9ytnkXNlCLPMuQD22FEKI+VDdk05nJK6MapC1GIAFfCFH/zCWJzSWIUbWvWikcc2OlLgn4Qoj6FeGgvtQVEY63i0ZlyIDfVCdeCSGwLUncjNeNbYrkoa0QQrQQhuzhG1lOTg7/+c9/uP/++xv8WlOmTOHw4cNMnz6dOXPmNPj17K1du5bk5GSWLVvm9JikpCSWLl3K1q1b2bx5M0eOHGH+/PmN2Mraefzxx2ndujV/+9vf3DpGNG+WxdLrmacrdxqyh29kOTk5rFixwuG+upY8cOTs2bPs27ePAwcOuBzs6/P6dTFu3DhDB3shXDG6PJBuZd71fl4jVO5ssj389/6xr8GvMfnvA6tsmz9/Punp6URFRTFy5EhuvPFGHn30UUJCQjh27BhpaWnccsstZGZmUlxczIMPPsjMmTMB08IpDz74IFu3bqVVq1Z8+OGHdOzYkfXr1/PEE0/g7e1NcHAwO3fuZNSoUZw+fZqoqChefvll2rRpw6xZsygsLKRHjx6sXr2akJAQEhMTiYqKYvfu3UyZMoUtW7YQHR3Nrl27KCgo4M033+Tpp5/m4MGD3HbbbTz55JMAvPXWW7z00kuUlJSQkJDAihUr8Pb2Zs2aNTz99NO0a9eOAQMG2NQHqon1N4Jp06bRtm1bkpOTOXv2LM8995xlsZglS5bw3nvvcfHiRcaPH88TTzxR7XkTExNd+p2ef/55Vq9eDcDdd9/NQw89BJiqhP7rX/+iQ4cOdOnShdjYWMBUkfRPf/oT58+fJzAwkFWrVtG7d2+Xf1/RPM0ICGEGITB9Q72e1wiVO6WHX0vPPPMMPXr0IDU1lSVLlgCwf/9+XnzxRUuNnNWrV5OSkkJycjIvvfQSFy5cAEx1dK6++mq+++47hg4dyqpVqwBT1cxt27bx3XffWUotbN682XKdIUOG8Mc//pFnn32WAwcOEBERYRMkS0pKSE5O5q9//SsAfn5+JCcnM2vWLG6++WaWL1/OoUOHWLt2LRcuXODo0aOsW7eOL7/8ktTUVLy9vXn77bf56aefWLRoEV9++SW7d+/myJEjbn1WP/30E7t372br1q2Wnv/27ds5fvw4e/fuJTU1lZSUFHbu3AnADTfcwJkzZxyeq6bfKSUlhTVr1rBnzx6++eYbVq1axbfffktKSgrvvvsuqampfPzxx+zbd6mjMHPmTF5++WVSUlJYunRpo6TphPCkRuvhK6X6AA8ClwOfaq1faaxrN7T4+HjCwsIsr1966SU2btwIQGZmJsePHyc0NBQ/Pz/LUoKxsbF88sknAAwePJhp06YxefJkJkyYUOX8ubm55OTkMGzYMMBUCnnSpEmW/dY17wFLff2IiAj69etHp06dAOjevTuZmZns3r2blJQUBg40fYMpKiqiQ4cO7Nmzh8TERNq3b285rzuF3m655Ra8vLzo27evpSLn9u3b2b59O9HR0QDk5+dz/Phxhg4dyscff+z0XK78TuPHjycoyDRNfsKECezatYuKigrGjx9PYGCgzXny8/P56quvbD5H8wIyQjRXLgV8pdRqYCxwTmvd32r7GOBFwBt4XWv9jLNzaK2PArOUUl7Am0CzCfjmIAOmh5g7duzg66+/JjAwkMTERIqLiwFsyhybyycDvPrqq+zZs4ePPvqI2NhYUlJS6nx9sC2tbJ2SMZdW1lozdepUnn76aZv3bdq0qVbXrYn1tc01m7TWPPLII9x7b+0eXNX0O9VWRUUF7dq1IzU1tdbvFaIuQsov0LYip2p9odoyT2arw3lcTemsBcZYb1BKeQPLgeuBvsAUpVRfpVSEUmqr3U+HyveMAz4CnHflDK5Nmzbk5eU53Z+bm0tISAiBgYEcO3aMb775psZzpqenk5CQwOLFi2nfvj2ZmZk2+4ODgwkJCWHXrl3ApVLIdXXdddfx/vvvc+7cOQCys7M5efIkCQkJfPHFF1y4cIHS0lLWr19vec/GjRt55JFH6nxNs9GjR7N69WrLMo+nT5+2tMMdQ4YMYdOmTRQWFlJQUMDGjRsZMmQIQ4cOZdOmTRQVFZGXl8eWLVsAaNu2LWFhYZbfUWvNd99953Y7hHCmbUUOAbrYo21wqYevtd6plOpmtzkeOKG1/gFAKfUucLPW+mlM3wYcnWczsFkp9RHwH0fHKKVmAjMBOnbsSFJSks3+4OBg8vLyuP6Bhn+45iiw+/n5ER8fT9++fRk5ciSjR4+mrKzMcuzgwYNZtmwZ4eHh9OzZk4EDB1JYWGjZb/5vUVERpaWl5OXlMWfOHNLT09FaM2zYMLp3786pU6eoqKiwHL9ixQoeeughioqK6NatGytWrCAvL4/y8nIKCgosx1m/LiwstGmbeV94eDgLFixgxIgRVFRU4Ovry9KlS4mPj2f+/PkkJCQQHBxMZGQkJSUl5OXlceTIEfz9/at8JtbXKC4uthxfWlpqCbLWn+egQYOYMGECCQkJgOnbyapVq2jVqhUTJ05k2bJllnSNmSu/U0xMDFOmTCEuzlQg8I9//CPmiXu33HILERERtG/fnqioKC5evEheXh6vvfYac+bMYfHixZSWljJx4kS6d+/OxYsX8fX1rfbGDqYVw6z/PvPz86v8vbZ0UTk5AKTW4nNx5T11OW99Xr8utNYU4U9S2Fy3zhOVswCA1GrP47hP7XJ55MqAv9Wc0lFK3QqM0VrfXfn6TiBBaz3byfsTgQmAP3BAa728pmvGxcXp5ORkm21SHtkz7rjjDl544QVLfl9IeWSX1GWmbT2WR66TBjq3uSjbezNT3TuRC+1zVh650R7aaq2TgCRXjpXSCsbz1ltveboJQgg3uTMs8zTQxep158ptQgghDMidgL8P6KmUClNK+QG3A3Vbr8+O1NIRQjQJa250f9RNI3Ip4Cul3gG+BsKVUllKqRla6zJgNrANOAq8p7U+XB+NUkrdpJRamZubWx+nE0IIgeujdKY42f4xDTDEUmu9BdgSFxd3T32fWwghWipDllaQHr4QQtQ/Qwb8ppjD79atG7/88kuV7Zs3b+aZZ0wTkM+fP09CQoKlEJizqpvVefzxx1m6dGm12x977DF27NhR63M3psTEROyH3NblGCGE6wwZ8JtTD9+6ZPCnn35KREQE3377LV26dKlTwHfF4sWLGTFiRIOcWwjRdBky4Bu5h19QUMCNN97IgAED6N+/P+vWrbPse/nll4mJiSEiIoJjx44BppLBs2fPJjU1lXnz5vHhhx8SFRXFww8/bCmzPHeuacbckiVLGDhwIJGRkSxatMhy3qeeeopevXpxzTXX8P3339fYxmnTpvH+++8Dpm8eixYtqtKugoIC7rrrLuLj44mOjubDDz+s8bytW7dm7ty59OvXjxEjRrB3714SExPp3r27pcpncXEx06dPJyIigujoaD7//HPANLP49ttvp0+fPowfP56ioiLLebdv386gQYOIiYlh0qRJlrILQoj61WTr4QN8vnYl507+UK/n7NC1O9dOm+l0///+9z+uvPJKPvrINMvN+lvI5Zdfzv79+1mxYgVLly7l9ddft+yLiopi8eLFlnrxGRkZHD582FK8y7pssNaacePGsXPnToKCgizlfcvKyoiJibHUc3eVo3Y99dRTDB8+nNWrV5OTk0N8fDwjRowgNzeXu+++22HlyoKCAoYPH86SJUsYP348Cxcu5JNPPuHIkSNMnTqVcePGsXz5cpRSHDx4kGPHjjFq1CjS0tJ45ZVXCAwM5OjRoxw4cICYmBgAfvnlF5588kl27NhBUFAQzz77LM8//zyPPfZYrX5HIUTNDBnwjTzTNiIigr/+9a88/PDDjB07liFDhlj2mUsbx8bG8sEHH9TqvM7KBufl5Tks71sbjtq1fft2Nm/ebMn7FxcXc+rUKfr06eO0TLGfnx9jxphq6EVERODv74+vry8RERFkZGQAsHv3bh544AEAevfuTdeuXUlLS2Pnzp38+c9/BiAyMpLIyEgAvvnmG44cOcLgwYMBU23/QYMG1fp3FELUzJAB39VhmdX1xBtKr1692L9/Px9//DELFy7kuuuus/RGzWV7rUsfu8pZ2eB//vOfbrfZUbu01mzYsIHw8HCXz2Nd3tm6THFdSxSb2zFy5EjeeeedOr1fCOE6Q+bwjezMmTMEBgZyxx13MHfuXPbv31+n89iXWXZWNthZeV93jR49mpdfftlSp/7bb7+tl/MOGTKEt99+G4C0tDROnTpFeHg4Q4cO5T//MRVIPXToEAcOmGp6X3311Xz55ZecOHECMKWN3Fl0RQjhnCF7+EZ28OBB5s6di5eXF76+vrzySt3WcQkNDWXw4MH079+f66+/niVLlnD06FFLOqN169a89dZbxMTEcNtttzFgwAA6dOhgWaXKXY8++igPPfQQkZGRVFRUEBYWxtatWzlz5ozTHL4r7r//fu677z4iIiLw8fFh7dq1+Pv7c9999zF9+nT69OlDnz59LM8h2rdvz9q1a5kyZYplxaknn3ySXr161cvvKYSRZPiUu7+2bUm56b91OI/L5ZEbk1UO/57jx4/b7JPyyMIopDyyC5p7eeRaHPvGK79nm3chJ317uNE4oKTA9F+/IKeH7L3rA8+WR64NKa0ghGhuZgSEMIMQmL7BvRNZbjLOz6PuUg63Sw5fCCFaiCYZ8I2YhhIti/wNiqaoyQX8gIAALly4IP+HEx6jtebChQsEBAR4uilC1Iohc/jVTbzq3LkzWVlZnD9/vvEbJkSlgIAAOnfu7OlmCFErhgz41T209fX1JSwszAOtEkKIpq3JpXSEEELUjQR8IYRoISTgCyFEC2HImbZmSqnzwMlGuFQwYLTVVhqrTfV9nfo4X13PUdv3uXq8q8ddDlRd9qxlkP8PGet8PbXWVRcU0Vq3+B9gpafb4Kk21fd16uN8dT1Hbd/n6vG1OC65sf4+jPYj/x8y1vmcnUNSOib1U4KyfjVWm+r7OvVxvrqeo7bvc/V4I/59GI0RPyP5/5AdQ6d0hGhKlFLJ2kHBKiGMQnr4QtSflZ5ugBDVkR6+EEK0ENLDF0KIFkICvhBCtBAS8IUQooUwZPE0IZoDpdQtwI1AW+ANrfV2z7ZItHTSwxeiFpRSq5VS55RSh+y2j1FKfa+UOqGUmg+gtd6ktb4HmAXc5on2CmFNAr4QtbMWGGO9QSnlDSwHrgf6AlOUUn2tDllYuV8Ij5KAL0QtaK13Atl2m+OBE1rrH7TWJcC7wM3K5Fngv1rr/Y3dViHsSQ5fCPddBWRavc4CEoAHgBFAsFLqd1rrVz3ROCHMJOAL0UC01i8BL3m6HUKYSUpHCPedBrpYve5cuU0IQ5GAL4T79gE9lVJhSik/4HZgs4fbJEQVEvCFqAWl1DvA10C4UipLKTVDa10GzAa2AUeB97TWhz3ZTiEckeJpQgjRQkgPXwghWggJ+EII0UJIwBdCiBZCAr4QQrQQEvCFEKKFkIAvhBAthAR8IYRoISTgCyFECyEBXwghWoj/D3aaO7Pwk0a2AAAAAElFTkSuQmCC", - "text/plain": [ - "
" - ] - }, - "metadata": { - "needs_background": "light" - } - } - ], - "metadata": {} + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## How Poissonian counts are read" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import ogip.spec\n", + "pha = ogip.spec.PHAI.from_file_name(\"tests/data/MOS1source_spectrum_150_rbn.pi\")\n", + "from astropy.io import fits as pf\n", + "import numpy as np\n", + "ff = pf.open(\"tests/data/MOS1source_spectrum_150_rbn.pi\")\n", + "counts = ff[1].data['COUNTS']\n", + "assert(np.abs(np.sum(counts - (pha._rate)*pha._exposure)) < 1e-12)\n", + "assert(np.sum(pha._rate*pha._exposure - (pha._stat_err*pha._exposure)**2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { - "orig_nbformat": 4, + "interpreter": { + "hash": "6c0c14c3bc85a7c6c1cb359b38aedade0f64d7248b8647c7b1bd974c3ac87f4c" + }, + "kernelspec": { + "display_name": "py10", + "language": "python", + "name": "py10" + }, "language_info": { - "name": "python", - "version": "3.8.11", - "mimetype": "text/x-python", "codemirror_mode": { "name": "ipython", "version": 3 }, - "pygments_lexer": "ipython3", + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", "nbconvert_exporter": "python", - "file_extension": ".py" - }, - "kernelspec": { - "name": "python3", - "display_name": "Python 3.8.11 64-bit ('3.8.11': pyenv)" - }, - "interpreter": { - "hash": "6c0c14c3bc85a7c6c1cb359b38aedade0f64d7248b8647c7b1bd974c3ac87f4c" + "pygments_lexer": "ipython3", + "version": "3.10.14" } }, "nbformat": 4, - "nbformat_minor": 2 -} \ No newline at end of file + "nbformat_minor": 4 +} From d004721c89a118e9dc56adc66f223c8c522a191a Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:14:23 +0200 Subject: [PATCH 05/11] Binary file for test --- tests/data/MOS1source_spectrum_150_rbn.pi | 50 +++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 tests/data/MOS1source_spectrum_150_rbn.pi diff --git a/tests/data/MOS1source_spectrum_150_rbn.pi b/tests/data/MOS1source_spectrum_150_rbn.pi new file mode 100644 index 0000000..5aba1e4 --- /dev/null +++ b/tests/data/MOS1source_spectrum_150_rbn.pi @@ -0,0 +1,50 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extensions XPROC0 = 'arfgen spectrumset=MOS1source_spectrum_150.fits rmfset=response.ds &'CONTINUE 'withrmfset=no arfset=deletearf.ds detmaptype=flat detmaparray=detma&'CONTINUE 'pfile.ds: detxoffset=1200 detyoffset=1200 withdetbounds=no detxbins&'CONTINUE '=1 detybins=1 withdetbins=yes psfenergy=2 filterdss=yes filteredset&'CONTINUE '=filteredpixellist.ds withfilteredset=no sourcecoords=eqpos sourcex&'CONTINUE '=0 sourcey=0 withsourcepos=no extendedsource=yes modeleffarea=no mo&'CONTINUE 'delquantumeff=no modelfiltertrans=no modelcontamination=yes modelee&'CONTINUE '=no modelootcorr=yes applyxcaladjustment=no eegridfactor=100 withba&'CONTINUE 'dpixcorr=yes badpixlocation=MOS1clean.fits psfmodel=ELLBETA badpixe&'CONTINUE 'lresolution=2 withbadpixres=no setbackscale=yes keeparfset=no useod&'CONTINUE 'fatt=no ignoreoutoffov=yes crossreg_spectrumset='''' crossregionarf&'CONTINUE '=no # (arfgen-1.90.4) [xmmsas_20141104_1833-14.0.0]' XPROC1 = 'arfgen spectrumset=MOS1source_spectrum_150.fits rmfset=response.ds &'CONTINUE 'withrmfset=no arfset=deletearf.ds detmaptype=flat detmaparray=detma&'CONTINUE 'pfile.ds: detxoffset=1200 detyoffset=1200 withdetbounds=no detxbins&'CONTINUE '=1 detybins=1 withdetbins=yes psfenergy=2 filterdss=yes filteredset&'CONTINUE '=filteredpixellist.ds withfilteredset=no sourcecoords=eqpos sourcex&'CONTINUE '=0 sourcey=0 withsourcepos=no extendedsource=yes modeleffarea=no mo&'CONTINUE 'delquantumeff=no modelfiltertrans=no modelcontamination=yes modelee&'CONTINUE '=no modelootcorr=yes applyxcaladjustment=no eegridfactor=100 withba&'CONTINUE 'dpixcorr=yes badpixlocation=MOS1clean.fits psfmodel=ELLBETA badpixe&'CONTINUE 'lresolution=2 withbadpixres=no setbackscale=yes keeparfset=no useod&'CONTINUE 'fatt=no ignoreoutoffov=yes crossreg_spectrumset='''' crossregionarf&'CONTINUE '=no # (arfgen-1.90.4) [xmmsas_20141104_1833-14.0.0]' XDAL0 = 'MOS1source_spectrum_150.fits 2015-04-01T11:41:32.000 Modify arfgen &'CONTINUE '(arfgen-1.90.4) [xmmsas_20141104_1833-14.0.0] High SAS_MEMORY_MODEL&'CONTINUE '= SAS_ROWS= SAS_ZERO_ROWS= SAS_COLUMN_WISE= ' CREATOR = 'evselect (evselect-3.62) [xmmsas_20141104_1833-14.0.0]' / name of codDATE = '2015-04-01T11:41:28.000' / creation date LONGSTRN= 'OGIP 1.0' DATAMODE= 'IMAGING ' / Instrument mode (IMAGING or TIMING) TELESCOP= 'XMM ' / XMM mission INSTRUME= 'EMOS1 ' / EPIC MOS Instrument OBS_ID = '0745240201' / Observation Identifier EXP_ID = '0745240201001' / Exposure Identifier DATE-OBS= '2014-05-29T23:19:08' / Start time of exposure DATE-END= '2014-05-30T08:23:24' / End time of exposure OBS_MODE= 'POINTING' / Observation mode (pointing or slew) REVOLUT = 2650 / Revolution number OBJECT = 'EXO 2030+375' / Name of observed object OBSERVER= 'Dr Carlo Ferrigno' / Name of observer RA_OBJ = 3.08063707500000E+02 / [deg] RA of target DEC_OBJ = 3.76375000000000E+01 / [deg] Dec of target RA_NOM = 3.08063707500000E+02 / [deg] RA of nominal boresight DEC_NOM = 3.76375000000000E+01 / [deg] Dec of nominal boresight EXPIDSTR= 'S001 ' / Exposure ID FILTER = 'Thin1 ' / Filter ID ATT_SRC = 'AHF ' / Source of attitude data (AHF|RAF|THF) ORB_RCNS= T / Reconstructed orbit data used? TFIT_RPD= F / Recalculated signal propagation delays? TFIT_DEG= 4 / Degree of OBT-MET fit polynomial TFIT_RMS= 7.21618977381168E-07 / RMS value of OBT-MET polynomial fit TFIT_PFR= 0.00000000000000E+00 / Fraction of disregarded TCS data points TFIT_IGH= T / Ignored GS handover in TC data? SUBMODE = 'PrimePartialW2' / Guest Observer mode EQUINOX = 2.00000000000000E+03 / Equinox for sky coordinate x/y axes RADECSYS= 'FK5 ' / World coord. system for this file REFXCTYP= 'RA---TAN' / WCS Coord. type: RA tangent plane projection REFXCRPX= 25921 / WCS axis reference pixel of projected image REFXCRVL= 3.08029916666667E+02 / [deg] WCS coord. at X axis ref. pixel REFXCDLT= -1.38888888888889E-05 / [deg/pix] WCS X increment at ref. pixel REFXLMIN= 1 / WCS minimum legal QPOE projected image X axis vREFXLMAX= 51840 / WCS maximum legal QPOE projected image X axis vREFXDMIN= 29567 / WCS minimum projected image X axis data value REFXDMAX= 47911 / WCS maximum projected image X axis data value REFXCUNI= 'deg ' / WCS Physical units of X axis REFYCTYP= 'DEC--TAN' / WCS Coord. type: DEC tangent plane projection REFYCRPX= 25921 / WCS axis reference pixel of projected image REFYCRVL= 3.76486111111111E+01 / [deg] WCS coord. at Y axis ref. pixel REFYCDLT= 1.38888888888889E-05 / [deg/pix] WCS Y increment at ref. pixel REFYLMIN= 1 / WCS minimum legal QPOE projected image Y axis vREFYLMAX= 51840 / WCS maximum legal QPOE projected image Y axis vREFYDMIN= 14128 / WCS minimum projected image Y axis data value REFYDMAX= 32485 / WCS maximum projected image Y axis data value REFYCUNI= 'deg ' / WCS Physical units of Y axis AVRG_PNT= 'MEDIAN ' / Meaning of PNT values (mean or median) RA_PNT = 3.08029916666667E+02 / [deg] Actual (mean) pointing RA of the optical DEC_PNT = 3.76486111111111E+01 / [deg] Actual (mean) pointing Dec of the opticalPA_PNT = 5.54859504699707E+01 / [deg] Actual (mean) measured position angle of CONTENT = 'EPIC MOS IMAGING MODE EVENT LIST' HISTORY Created by evselect (evselect-3.62) [xmmsas_20141104_1833-14.0.0] at 201HISTORY 5-04-01T11:41:28 HISTORY Modified by arfgen (arfgen-1.90.4) [xmmsas_20141104_1833-14.0.0] at 2015HISTORY -04-01T11:41:28 HISTORY Modified by arfgen (arfgen-1.90.4) [xmmsas_20141104_1833-14.0.0] at 2015HISTORY -04-01T11:41:32 END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 10 / width of table in bytes NAXIS2 = 2400 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 4 / number of fields in each row TTYPE1 = 'CHANNEL ' / Pulse Invarient (PI) Channel TFORM1 = 'I ' / data format of field: 2-byte INTEGER TTYPE2 = 'COUNTS ' / Counts per channel TFORM2 = 'J ' / data format of field: 4-byte INTEGER TUNIT2 = 'count ' / physical unit of field TTYPE3 = 'QUALITY ' / Quality flag of this channel (0=good) TFORM3 = 'I ' / data format of field: 2-byte INTEGER TTYPE4 = 'GROUPING' / Grouping flag for channel (0=undefined) TFORM4 = 'I ' / data format of field: 2-byte INTEGER EXTNAME = 'SPECTRUM' / The name of this table HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'SPECTRUM' / PHA dataset (OGIP memo OGIP-92-007) HDUVERS1= '1.1.0 ' / Obsolete - included for backwards compatibilityHDUVERS = '1.1.0 ' / Version of format (OGIP memo OGIP-92-007) HDUCLAS2= 'TOTAL ' / Gross PHA Spectrum (source + bkgd) HDUCLAS3= 'COUNT ' / PHA data stored as Counts (not count/s) TLMIN1 = 0 / Lowest legal channel number TLMAX1 = 2399 / Highest legal channel number TELESCOP= 'XMM ' / Telescope (mission) name INSTRUME= 'EMOS1 ' / Instrument name FILTER = 'Thin1 ' / Instrument filter in use EXPOSURE= 3.167416E+04 / Exposure time AREASCAL= 1.000000E+00 / Area scaling factor BACKFILE= 'MOS1back_spectrum.fits' / Background FITS file BACKSCAL= 4.684692E+06 / Background scale factor CORRFILE= 'NONE ' / Correlation FITS file CORRSCAL= 1.000000E+00 / Correlation scale factor RESPFILE= 'MOS1_150.rmf' / redistribution matrix ANCRFILE= 'MOS1_150.arf' / ancillary response PHAVERSN= '1992a ' / obsolete DETCHANS= 2400 / total number possible channels CHANTYPE= 'PI ' / channel type (PHA, PI etc) POISSERR= T / Poissonian errors to be assumed STAT_ERR= 0 / no statistical error specified SYS_ERR = 0 / no systematic error specified HISTORY FITS SPECTRUM extension written by UPDPHA 1.0.0 SLCTEXPR= '(FLAG & 0x766ba000) == 0 && (PATTERN<=12) && ((X,Y) IN ANNULUS(2404&'CONTINUE '0,25080,150.0,1750.0))' / Filtering expression used by evselect SPECDELT= 5 / Spectral channel size, ie the binning factor SPECPIX = 0 / The rebinned channel correspondsing to SPECVAL SPECVAL = 2.00000000000000E+00 / Original chan value at center of rebinned chan DSTYP1 = 'CCDNR ' / data subspace descriptor: name DSTYP2 = 'FLAG ' / data subspace descriptor: name DSTYP3 = 'FLAG ' / data subspace descriptor: name DSTYP4 = 'PI ' / data subspace descriptor: name DSUNI4 = 'CHAN ' / data subspace descriptor: units DSTYP5 = 'TIME ' / data subspace descriptor: name DSUNI5 = 's ' / data subspace descriptor: units DSTYP6 = 'FLAG ' / data subspace descriptor: name DSTYP7 = 'PATTERN ' / data subspace descriptor: name DSTYP8 = 'POS(X,Y)' / data subspace descriptor: name DSVAL1 = '1 ' / data subspace descriptor: value DSFORM2 = 'X ' / data subspace descriptor: data type DSVAL2 = 'b000x00xxx0x0x0x0x0xxxxxxxxxxxxx' / data subspace descriptor: value DSFORM3 = 'X ' / data subspace descriptor: data type DSVAL3 = 'b000x00xx00x0x000x0xxxxxxxxxxxxx' / data subspace descriptor: value DSVAL4 = '(150: ' / data subspace descriptor: value DSVAL5 = 'TABLE ' / data subspace descriptor: value DSREF5 = ':GTI00004' / data subspace descriptor: reference DSFORM6 = 'X ' / data subspace descriptor: data type DSVAL6 = 'b000x00xx00x0x000x0xxxxxxxxxxxxx' / data subspace descriptor: value DSVAL7 = ':12 ' / data subspace descriptor: value DSVAL8 = 'TABLE ' / data subspace descriptor: value DSREF8 = ':REG00108' / data subspace descriptor: reference 2DSVAL1 = '2 ' / data subspace descriptor: value 2DSREF5 = ':GTI00104' / data subspace descriptor: reference 3DSVAL1 = '4 ' / data subspace descriptor: value 3DSREF5 = ':GTI00204' / data subspace descriptor: reference 4DSVAL1 = '5 ' / data subspace descriptor: value 4DSREF5 = ':GTI00304' / data subspace descriptor: reference 5DSVAL1 = '7 ' / data subspace descriptor: value 5DSREF5 = ':GTI00404' / data subspace descriptor: reference MTYPE1 = 'POS ' / DM Keyword: Descriptor name MFORM1 = 'X,Y ' CREATOR = 'grppha 3.0.1' / s/w task which wrote this dataset DATE = '2016-04-01T06:43:42' / file creation date (YYYY-MM-DDThh:mm:ss UT) HISTORY This extension has been written by WT_SPEC Ver 1.7.0 HISTORY The original fits file was ../MOS1source_spectrum_150.fits END   +    !"# $ %& ' ( ) * +, -./ 01 2 3 45 678 9 :;<=>? @AB C DEF GH +I J KLMNOPQRSTUVWX YZ [ \] +^ +_`ab cd eџџf џџgџџhџџiџџjџџklџџm џџnџџoџџpqџџrџџsџџtџџuџџv +wџџxџџyџџzџџ{џџ| } џџ~ +џџ џџ€џџџџ‚ ƒ џџ„џџ… џџ†џџ‡ +ˆџџ‰ +џџŠ +џџ‹џџŒџџ Žџџ џџџџ‘џџ’ +џџ“” џџ•џџ– +џџ—џџ˜џџ™šџџ› џџœ +џџџџžџџŸ  џџЁџџЂ +џџЃџџЄџџЅІџџЇ џџЈџџЉџџЊџџЋ +Ќ џџ­џџЎџџЏџџА џџБ В +џџГџџДџџЕ џџЖџџЗ ИџџЙ џџКџџЛ џџМ џџН +џџОПџџРџџС џџТџџУџџФ ХџџЦџџЧ џџШџџЩџџЪЫ џџЬџџЭ џџЮџџЯџџабџџвџџгџџдџџе џџжзџџиџџйџџкџџлџџмџџноџџпџџр џџсџџт!џџуфџџхџџц&џџч)џџш&џџщџџъ$ыџџь+џџэ$џџю6џџяџџ№'џџё+ђ'џџѓ'џџє'џџѕ2џџі+џџї/џџј4љ8џџњ@џџћ-џџќ;џџ§'џџў4џџџ56џџ2џџCџџ<џџCџџ4џџ<6џџ1џџ Fџџ +6џџ Dџџ Eџџ JLџџ@џџMџџPџџTџџ>џџ?`џџGџџFџџTџџLџџSџџDVџџ[џџYџџRџџ Lџџ!Qџџ"U#mџџ$Zџџ%Rџџ&Uџџ'Vџџ(Vџџ)h*bџџ+cџџ,Xџџ-Xџџ.nџџ/oџџ0}1pџџ2dџџ3‚џџ4wџџ5^џџ6Šџџ7z8nџџ9kџџ:nџџ;lџџ<mџџ=џџ>xџџ?„@„џџAџџB~џџC‚џџD{џџEkџџF†џџGƒHŒџџI|џџJ‡џџKoџџL•џџM}џџN~џџOPzџџQŠџџRџџSџџTџџU‘џџV џџW‰X€џџY–џџZ”џџ[žџџ\Ѕџџ]Јџџ^›џџ_›`џџaЂџџb‘џџcŸџџdЋџџežџџf‹џџgУhЊџџiџџj€џџkŠџџlџџm™џџnŠџџoЄp“џџq…џџr™џџs†џџtџџu›џџvџџw‰xwџџyˆџџzџџ{†џџ|ˆџџ}€џџ~ џџt€šџџžџџ‚џџƒџџ„•џџ…™џџ†Ž‡Ѕџџˆ–џџ‰–џџŠ„џџ‹…џџŒŽџџ‡џџŽœџџ™џџ‘‡џџ’”џџ“~џџ”Љџџ•Џџџ––—›џџ˜’џџ™’џџšџџ›ЇџџœДџџЧџџžœџџŸЁ ЄџџЁЃџџЂІџџЃžџџЄЉџџЅŸџџІЕџџЇМџџЈРЉЄџџЊИџџЋЂџџЌЌџџ­­џџЎЋџџЏЄџџАЄБІџџВДџџГРџџДЕџџЕВџџЖЗџџЗЉџџИЦџџЙДКФџџЛЄџџМГџџНАџџОИџџПЅџџР­џџСЁџџТŒУ†џџФŽџџХ”џџЦЅџџЧˆџџШ•џџЩ”џџЪnџџЫœЬœџџЭ€џџЮŽџџЯџџа‘џџб‘џџвrџџгŒџџдzеŒџџжtџџз„џџи’џџй”џџк™џџлџџмˆн†џџо‰џџп„џџр†џџсšџџтŽџџу”џџф€џџх˜ц˜џџч–џџшџџщЅџџъ…џџы†џџьЋџџэџџюя‹џџ№ˆџџё‡џџђЅџџѓ…џџєŸџџѕ џџіЎџџї‰јЄџџљŸџџњ–џџћІџџќžџџ§žџџўŸџџџЅЈџџ›џџšџџ“џџЁџџџџЎџџŸџџЂџџ ’ +šџџ ˜џџ ”џџ ЁџџАџџџџАџџЁџџ•ВџџЄџџŸџџ›џџЉџџ–џџ›џџ•џџ­џџЛЄџџЉџџТџџ Еџџ!“џџ"žџџ#Љџџ$­џџ%­џџ&Џ'Иџџ(žџџ)Бџџ*™џџ+Џџџ,Мџџ-Єџџ.Іџџ/Ѕ0–џџ1Ќџџ2Ёџџ3Њџџ4жџџ5Јџџ6Ѕџџ7Йџџ8Ёџџ9­:Жџџ;Јџџ<Ѕџџ=Јџџ>Їџџ?Аџџ@ЄџџAЦџџBШџџCЦDЎџџEЦџџF­џџGџџHГџџIИџџJЯџџKВџџLЏMМџџNЈџџOМџџPЦџџQГџџRЏџџSšџџTЬџџUГџџVАWКџџXЫџџYМџџZНџџ[Фџџ\Ъџџ]Иџџ^Жџџ_Щ`вџџaБџџbНџџcЎџџdФџџeПџџfИџџgРџџhЙџџiЗjЕџџkеџџlОџџmИџџnДџџoЄџџpжџџqЮџџrзџџsТtЩџџuШџџvдџџwЮџџxЦџџyНџџzЪџџ{Юџџ|Цџџ}в~Пџџйџџ€ЧџџСџџ‚УџџƒЗџџ„иџџ…Юџџ†Цџџ‡РˆФџџ‰НџџŠШџџ‹ЬџџŒФџџжџџŽФџџЗџџЩџџ‘У’Пџџ“Ыџџ”оџџ•Ъџџ–Йџџ—иџџ˜зџџ™Хџџšоџџ›сџџœааџџžПџџŸЭџџ КџџЁГџџЂхџџЃЦџџЄРџџЅНџџІРЇКџџЈАџџЉяџџЊЛџџЋфџџЌфџџ­УџџЎаџџЏЧџџАФџџБшВмџџГЫџџДвџџЕПџџЖЖџџЗЫџџИжџџЙБџџКНџџЛУџџМдНТџџОбџџПЫџџРЧџџССџџТеџџУуџџФХџџХРџџЦСЧНџџШЬџџЩЗџџЪМџџЫлџџЬьџџЭЖџџЮЈџџЯЮџџажџџбгвдџџгЙџџдУџџеЯџџжСџџзПџџиЬџџймџџкАџџлТмЭџџннџџощџџпЛџџртџџсТџџтПџџуЯџџфЫџџх№џџцйчРџџшМџџщеџџъЫџџысџџьмџџэЮџџюхџџяЭџџ№кёФџџђюџџѓХџџєЕџџѕЯџџідџџїєџџјеџџљЛџџњЖџџћЩќфџџ§иџџўЬџџџуџџЮџџЩџџЗџџЛџџрџџлџџЩУџџРџџ Ьџџ +Тџџ Фџџ Юџџ гџџйџџШџџвџџТПџџфџџзџџПџџЎџџхџџмџџЭџџСџџУџџЧЗџџЕџџМџџ Эџџ!уџџ"Яџџ#Пџџ$Оџџ%Чџџ&Ыџџ'Ц(Шџџ)Пџџ*Єџџ+Сџџ,Щџџ-вџџ.аџџ/Бџџ0Оџџ1лџџ2гџџ3Ч4Нџџ5Эџџ6Нџџ7™џџ8Гџџ9Пџџ:Кџџ;Рџџ<лџџ=Уџџ>Д?Кџџ@нџџAУџџBЮџџCЪџџDФџџEЦџџFоџџGФџџHХџџIЎџџJЙKПџџLЏџџMЩџџNЕџџOЛџџPФџџQЊџџRЙџџSЌџџTКџџUСVГџџWИџџXЬџџYЛџџZЫџџ[еџџ\Кџџ]Юџџ^Тџџ_Рџџ`­aИџџbДџџcЗџџdвџџeЛџџfфџџgЄџџhСџџi‘џџjЧџџkЎџџlТmЫџџnТџџoЦџџpУџџqИџџrІџџsЮџџtРџџuИџџvЛџџwШxЗџџyЛџџzЎџџ{šџџ|Дџџ}Ќџџ~БџџЗџџ€УџџДџџ‚ЛџџƒМ„Оџџ…Вџџ†Іџџ‡ЕџџˆЛџџ‰ЏџџŠВџџ‹НџџŒ­џџМџџŽКџџЁšџџ‘Ћџџ’Еџџ“Нџџ”Сџџ•Џџџ–Ќџџ—Бџџ˜Ўџџ™ЃџџšКџџ›ЏœЈџџСџџž­џџŸЋџџ ГџџЁžџџЂАџџЃЫџџЄЊџџЅ›џџІІџџЇ­ЈЖџџЉБџџЊЁџџЋŒџџЌКџџ­ЊџџЎЎџџЏРџџАžџџБЄџџВ—Г­џџДЇџџЕБџџЖНџџЗЖџџИІџџЙЏџџК–џџЛЖџџ̘џџНЂџџОœџџПВР џџСŸџџТІџџУГџџФМџџХЇџџЦЃџџЧЉџџШЏџџЩЃџџЪЌџџЫ’ЬЈџџЭ”џџЮЈџџЯЊџџаџџбˆџџвЊџџгšџџдЊџџеАџџжџџзРи‹џџйЏџџкЇџџлœџџмЊџџн€џџоšџџпЅџџр‚џџсŒџџтЂџџу—џџфЅх™џџцŽџџчœџџшЉџџщ“џџъŠџџы˜џџь—џџэИџџюІџџяЃџџ№†ёŽџџђ џџѓ“џџє•џџѕœџџі‰џџїЏџџјˆџџљ‘џџњ‰џџћ›џџќ§“џџўџџџ‹џџ„џџ‹џџšџџŒџџšџџŒџџ•џџЂџџ џџ +џџ {џџ Žџџ џџ”џџ‡џџˆџџ€џџŽџџwџџŽџџ†џџrџџ|џџ…џџtџџ‘џџ џџ…џџ‰џџpџџ џџ!t"џџ#ˆџџ$Іџџ%}џџ&‚џџ'jџџ(џџ)†џџ*Ђџџ+ˆџџ,wџџ-pџџ.x/•џџ0wџџ1wџџ2|џџ3pџџ4ƒџџ5ƒџџ6…џџ7vџџ8†џџ9ƒџџ:hџџ;ƒ<Šџџ=oџџ>hџџ?|џџ@{џџAwџџBzџџCzџџD€џџE…џџFџџGH|џџI…џџJhџџKџџLgџџMvџџNџџOzџџPyџџQvџџRmџџSqџџT‡U€џџVpџџWiџџXyџџY{џџZ|џџ[rџџ\dџџ]џџ^~џџ_iџџ`…џџaubџџcqџџd}џџe~џџf}џџgnџџhpџџigџџjaџџkpџџlrџџmmџџncoxџџpVџџqzџџr_џџssџџtuџџu‚џџvpџџweџџxnџџyrџџzuџџ{T|eџџ}gџџ~rџџhџџ€sџџjџџ‚dџџƒoџџ„cџџ…oџџ†Zџџ‡iџџˆm‰nџџŠkџџ‹VџџŒlџџOџџŽjџџ_џџnџџ‘hџџ’mџџ“xџџ”`џџ•sџџ–\—_џџ˜]џџ™bџџšXџџ›]џџœ]џџ^џџž_џџŸlџџ gџџЁfџџЂ_џџЃTЄQџџЅ`џџІrџџЇgџџЈXџџЉVџџЊnџџЋPџџЌVџџ­SџџЎ_џџЏ]џџАQБVџџВcџџГRџџДvџџЕHџџЖcџџЗ`џџИhџџЙQџџКRџџЛ_џџМYџџНWОVџџПTџџРNџџСdџџТrџџУhџџФeџџХeџџЦ]џџЧXџџШXџџЩiџџЪYџџЫ]ЬYџџЭYџџЮXџџЯWџџаNџџб]џџв[џџгKџџдOџџеNџџжKџџзWџџи]йeџџкEџџлUџџмJџџнVџџоKџџпPџџрKџџсZџџтEџџуUџџфSџџхSџџцWч[џџшOџџщTџџъGџџыPџџьSџџэ[џџюGџџяeџџ№MџџёRџџђGџџѓLєMџџѕSџџіAџџїLџџј]џџљ@џџњ>џџћDџџќRџџ§\џџўUџџџOџџQџџLQџџZџџ\џџWџџNџџTџџOџџ Pџџ +Dџџ Lџџ dџџ IџџXџџFPџџOџџEџџ=џџQџџNџџDџџPџџ;џџEџџ;џџKџџQ?џџOџџ@џџ Fџџ!Dџџ"9џџ#9џџ$?џџ%Sџџ&Hџџ'=џџ(=џџ)Bџџ*@+Nџџ,?џџ-7џџ.Gџџ/7џџ0?џџ1.џџ2Cџџ3@џџ4Dџџ5Iџџ6.џџ76џџ8891џџ:9џџ;Eџџ<;џџ=<џџ>4џџ?7џџ@FџџA3џџB/џџC:џџD=џџE5џџF:G9џџH3џџI2џџJ8џџK@џџL.џџM9џџN>џџO>џџP4џџQ7џџR<џџS1џџT4UBџџV?џџW,џџX2џџY>џџZ4џџ[;џџ\6џџ]0џџ^5џџ_6џџ`0џџa5џџb:c(џџd0џџe=џџf*џџg2џџh*џџiHџџj-џџk6џџl1џџm.џџn9џџo3џџp/џџq,r"џџs8џџt8џџu1џџv0џџw+џџx0џџy'џџz3џџ{)џџ|5џџ}1џџ~3џџ,€3џџ1џџ‚)џџƒ1џџ„1џџ…'џџ†.џџ‡,џџˆ-џџ‰/џџŠ7џџ‹*џџŒ/џџџџŽ&џџ0џџ‘&џџ’"џџ“)џџ”&џџ•,џџ–.џџ—!џџ˜3џџ™(џџš6џџ›#џџœ0џџž/џџŸ5џџ џџЁ"џџЂџџЃ+џџЄ.џџЅ(џџІ.џџЇ џџЈ-џџЉ!џџЊ"џџЋ%Ќ%џџ­(џџЎ.џџЏ#џџА џџБ%џџВ)џџГџџД$џџЕ$џџЖ)џџЗ&џџИ(џџЙК%џџЛ$џџМ$џџНџџО#џџП џџРџџСџџТџџУџџФџџХ$џџЦџџЧ+џџШЩџџЪџџЫ"џџЬџџЭџџЮџџЯ#џџа#џџбџџв џџгџџд!џџеџџжзџџи"џџйџџк&џџлџџм'џџн'џџоџџпџџрџџс$џџтџџуџџфџџх!цџџчџџшџџщ"џџъџџыџџь"џџэџџюџџя"џџ№џџёџџђџџѓє џџѕџџіџџїџџјџџљ"џџњџџћџџќџџ§џџўџџџџџџџџџ!џџџџџџџџџџџџ џџ +џџ џџ џџ џџџџџџ џџџџџџџџџџџџџџџџџџџџџџџџџџџџџџ !џџ"џџ#џџ$џџ%џџ&џџ'џџ(џџ)џџ*џџ+џџ,џџ-џџ.џџ/џџ01џџ2џџ3џџ4џџ5 џџ6џџ7џџ8џџ9џџ:џџ;џџ<џџ=џџ>џџ?@џџAџџBџџCџџDџџEџџFџџGџџHџџIџџJџџKџџLџџMџџN O џџPџџQџџRџџSџџTџџUџџVџџWџџX џџYџџZ џџ[џџ\џџ]џџ^ _џџ`џџa џџb џџc џџdџџeџџfџџg џџhџџi џџjџџkџџlџџmnџџoџџpџџq џџr +џџsџџtџџuџџv џџw џџxџџyџџz џџ{џџ| } џџ~ џџџџ€џџџџ‚џџƒџџ„џџ…џџ†џџ‡џџˆ џџ‰џџŠџџ‹џџŒџџŽџџ џџ џџ‘ +џџ’џџ“џџ”џџ•џџ–џџ—џџ˜ џџ™ џџšџџ›œџџ џџžџџŸ џџ џџЁ +џџЂ +џџЃ џџЄ џџЅџџІџџЇџџЈџџЉ +џџЊЋџџЌџџ­џџЎ џџЏ џџА џџБ џџВ џџГџџД џџЕ џџЖџџЗџџИџџЙџџКЛ џџМ џџН џџО џџП џџРџџС џџТ џџУџџФџџХ џџЦџџЧ џџШ џџЩ Ъ џџЫџџЬ џџЭџџЮ џџЯ џџаџџбџџв +џџгџџдџџеџџжџџзџџиџџйк џџлџџм џџнџџоџџпџџрџџс џџтџџу џџф џџх +џџцџџч џџшџџщъџџыџџьџџэ +џџю џџяџџ№ џџёџџђ џџѓ џџєџџѕџџіџџїџџј џџљњ џџћ џџќџџ§џџўџџџџџ +џџ џџџџ +џџ џџ џџџџџџџџ  +џџ џџ џџ џџџџ +џџџџџџџџ џџџџ џџ џџџџџџџџ +џџџџџџџџ џџ џџ! +џџ"џџ#џџ$ џџ%џџ& џџ' +џџ(џџ)џџ* + џџ, џџ-џџ.џџ/џџ0џџ1 +џџ2џџ3џџ4џџ5 џџ6џџ7џџ8џџ9џџ:;џџ<џџ= џџ> +џџ?џџ@џџA џџBџџCџџD +џџEџџFџџGџџHџџI џџJKџџLџџMџџNџџOџџPџџQ џџR џџSџџTџџUџџVџџWџџXџџYџџZџџ[ \џџ] +џџ^ џџ_џџ`џџaџџb џџc џџdџџeџџfџџg џџh џџiџџjџџklџџmџџnџџoџџp +џџqџџrџџsџџtџџuџџvџџwџџx +џџyџџzџџ{|џџ}џџ~џџџџ€џџџџ‚џџƒџџ„џџ…џџ† џџ‡џџˆџџ‰џџŠџџ‹ŒџџџџŽ +џџџџ +џџ‘џџ’џџ“џџ”џџ•џџ–џџ—џџ˜џџ™џџš +џџ›œџџџџžџџŸ џџ џџЁџџЂџџЃџџЄџџЅџџІџџЇџџЈџџЉџџЊџџЋџџЌ­џџЎџџЏ џџАџџБџџВџџГџџДџџЕџџЖџџЗџџИџџЙџџКџџЛџџМНџџО џџП +џџР џџСџџТџџУџџФџџХ џџЦџџЧџџШџџЩџџЪџџЫџџЬџџЭЮџџЯџџабвгдежзийклмнопрстуфхцчшщъыьэюя№ёђѓєѕіїјљњћќ§ўџ  +    !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~€‚ƒ„…†‡ˆ‰Š‹ŒŽ‘’“”•–—˜™š›œžŸ ЁЂЃЄЅІЇЈЉЊЋЌ­ЎЏАБВГДЕЖЗИЙКЛМНОПРСТУФХЦЧШЩЪЫЬЭЮЯабвгдежзийклмнопрстуфхцчшщъыьэюя№ёђѓєѕіїјљњћќ§ўџ           +                       ! " # $ % & ' ( ) * + , - . / 0 1 2 3 4 5 6 7 8 9 : ; < = > ? @ A B C D E F G H I J K L M N O P Q R S T U V W X Y Z [ \ ] ^ _XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 16 / width of table in bytes NAXIS2 = 2 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 2 / number of fields in each row TTYPE1 = 'START ' / lower GTI boundary TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'STOP ' / upper GTI boundary TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field EXTNAME = 'GTI00004' / The name of this table CCDID = 1 / CCD ID HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'GTI ' / table contains Good Time Intervals HDUCLAS2= 'STANDARD' / standard Good Time Interval table ONTIME = 3.26511665188074E+04 / [s] sum of all Good Time Intervals TSTART = 5.17792813525046E+08 / [s] Lower bound of first GTI TSTOP = 5.17825468591598E+08 / [s] Uppler bound of last GTI TIMEUNIT= 's ' / All times in s unless specified otherwise TIMESYS = 'TT ' / XMM time will be TT (Terrestrial Time) MJDREF = 5.08140000000000E+04 / [d] 1998-01-01T00:00:00 (TT) expressed in MJD TIMEREF = 'LOCAL ' / Reference location of photon arrival times TASSIGN = 'SATELLITE' / Location of time assignment TIMEZERO= 0.00000000000000E+00 / [s] Clock correction (if not zero) CLOCKAPP= T / Clock correction applied END AОмф-†ihAОнV9@дAОнV='<ЋAОнcМ—rєXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 33 / width of table in bytes NAXIS2 = 1 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 5 / number of fields in each row TTYPE1 = 'SHAPE ' / shape of element TFORM1 = '16A ' / data format of field: ASCII Character TTYPE2 = 'X ' / x-coordinate vector TFORM2 = 'E ' / data format of field: 4-byte REAL TTYPE3 = 'Y ' / y-coordinate vector TFORM3 = 'E ' / data format of field: 4-byte REAL TTYPE4 = 'R ' / radius vector TFORM4 = '2E ' / data format of field: 4-byte REAL TTYPE5 = 'COMPONENT' / component number that shape belongs to TFORM5 = 'B ' / data format of field: BYTE EXTNAME = 'REG00108' / The name of this table HDUCLASS= 'ASC ' HDUCLAS1= 'REGION ' HDUCLAS2= 'STANDARD' MTYPE1 = 'pos ' MFORM1 = 'X,Y ' END ANNULUS FЛаFУ№CDкРXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 16 / width of table in bytes NAXIS2 = 1 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 2 / number of fields in each row TTYPE1 = 'START ' / lower GTI boundary TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'STOP ' / upper GTI boundary TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field EXTNAME = 'GTI00104' / The name of this table CCDID = 2 / CCD ID HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'GTI ' / table contains Good Time Intervals HDUCLAS2= 'STANDARD' / standard Good Time Interval table ONTIME = 3.26487833531499E+04 / [s] sum of all Good Time Intervals TSTART = 5.17792813534244E+08 / [s] Lower bound of first GTI TSTOP = 5.17825462317597E+08 / [s] Uppler bound of last GTI TIMEUNIT= 's ' / All times in s unless specified otherwise TIMESYS = 'TT ' / XMM time will be TT (Terrestrial Time) MJDREF = 5.08140000000000E+04 / [d] 1998-01-01T00:00:00 (TT) expressed in MJD TIMEREF = 'LOCAL ' / Reference location of photon arrival times TASSIGN = 'SATELLITE' / Location of time assignment TIMEZERO= 0.00000000000000E+00 / [s] Clock correction (if not zero) CLOCKAPP= T / Clock correction applied END AОмф-ˆФ2AОнcЖQNXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 16 / width of table in bytes NAXIS2 = 1 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 2 / number of fields in each row TTYPE1 = 'START ' / lower GTI boundary TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'STOP ' / upper GTI boundary TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field EXTNAME = 'GTI00204' / The name of this table CCDID = 4 / CCD ID HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'GTI ' / table contains Good Time Intervals HDUCLAS2= 'STANDARD' / standard Good Time Interval table ONTIME = 3.26487833932042E+04 / [s] sum of all Good Time Intervals TSTART = 5.17792813543401E+08 / [s] Lower bound of first GTI TSTOP = 5.17825462326795E+08 / [s] Uppler bound of last GTI TIMEUNIT= 's ' / All times in s unless specified otherwise TIMESYS = 'TT ' / XMM time will be TT (Terrestrial Time) MJDREF = 5.08140000000000E+04 / [d] 1998-01-01T00:00:00 (TT) expressed in MJD TIMEREF = 'LOCAL ' / Reference location of photon arrival times TASSIGN = 'SATELLITE' / Location of time assignment TIMEZERO= 0.00000000000000E+00 / [s] Clock correction (if not zero) CLOCKAPP= T / Clock correction applied END AОмф-‹\AОнcЖSЈбXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 16 / width of table in bytes NAXIS2 = 1 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 2 / number of fields in each row TTYPE1 = 'START ' / lower GTI boundary TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'STOP ' / upper GTI boundary TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field EXTNAME = 'GTI00304' / The name of this table CCDID = 5 / CCD ID HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'GTI ' / table contains Good Time Intervals HDUCLAS2= 'STANDARD' / standard Good Time Interval table ONTIME = 3.26487833531499E+04 / [s] sum of all Good Time Intervals TSTART = 5.17792813534244E+08 / [s] Lower bound of first GTI TSTOP = 5.17825462317597E+08 / [s] Uppler bound of last GTI TIMEUNIT= 's ' / All times in s unless specified otherwise TIMESYS = 'TT ' / XMM time will be TT (Terrestrial Time) MJDREF = 5.08140000000000E+04 / [d] 1998-01-01T00:00:00 (TT) expressed in MJD TIMEREF = 'LOCAL ' / Reference location of photon arrival times TASSIGN = 'SATELLITE' / Location of time assignment TIMEZERO= 0.00000000000000E+00 / [s] Clock correction (if not zero) CLOCKAPP= T / Clock correction applied END AОмф-ˆФ2AОнcЖQNXTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 16 / width of table in bytes NAXIS2 = 1 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 2 / number of fields in each row TTYPE1 = 'START ' / lower GTI boundary TFORM1 = 'D ' / data format of field: 8-byte DOUBLE TUNIT1 = 's ' / physical unit of field TTYPE2 = 'STOP ' / upper GTI boundary TFORM2 = 'D ' / data format of field: 8-byte DOUBLE TUNIT2 = 's ' / physical unit of field EXTNAME = 'GTI00404' / The name of this table CCDID = 7 / CCD ID HDUCLASS= 'OGIP ' / format conforms to OGIP standard HDUCLAS1= 'GTI ' / table contains Good Time Intervals HDUCLAS2= 'STANDARD' / standard Good Time Interval table ONTIME = 3.26487833532095E+04 / [s] sum of all Good Time Intervals TSTART = 5.17792813543401E+08 / [s] Lower bound of first GTI TSTOP = 5.17825462326755E+08 / [s] Uppler bound of last GTI TIMEUNIT= 's ' / All times in s unless specified otherwise TIMESYS = 'TT ' / XMM time will be TT (Terrestrial Time) MJDREF = 5.08140000000000E+04 / [d] 1998-01-01T00:00:00 (TT) expressed in MJD TIMEREF = 'LOCAL ' / Reference location of photon arrival times TASSIGN = 'SATELLITE' / Location of time assignment TIMEZERO= 0.00000000000000E+00 / [s] Clock correction (if not zero) CLOCKAPP= T / Clock correction applied END AОмф-‹\AОнcЖSІ2 \ No newline at end of file From 31d65660e8109c57a9b6e435a142c2787c9c29bb Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:14:58 +0200 Subject: [PATCH 06/11] improve logging and deal with Poissonian uncertainties --- ogip/spec.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/ogip/spec.py b/ogip/spec.py index 8d33783..c739aca 100644 --- a/ogip/spec.py +++ b/ogip/spec.py @@ -76,26 +76,28 @@ def from_file_name_osa(fn): except Exception as ex: logger.debug("failed to read from %s: %s", e, ex) - raise Exception("unable to read from any extension") + raise Exception("unable to read from any extension in OSA") @staticmethod def from_file_name_normal(fn): + logger.debug('Reading OGIP format') f = fits.open(fn) + + exposure = f['SPECTRUM'].header['EXPOSURE'] if 'RATE' in f['SPECTRUM'].data.names: rate = f['SPECTRUM'].data['RATE'] elif 'COUNTS' in f['SPECTRUM'].data.names: - rate = f['SPECTRUM'].data['COUNTS'] / f['SPECTRUM'].header['EXPOSURE'] + rate = f['SPECTRUM'].data['COUNTS'] / exposure else: rate = None if 'STAT_ERR' in f['SPECTRUM'].data.names: stat_err = f['SPECTRUM'].data['STAT_ERR'] - elif 'POISSERR' in f['SPECTRUM'].header and f['SPECTRUM'].header['POISSERR']: + elif f['SPECTRUM'].header.get('POISSERR', False): + logger.debug('Reading Poissonian uncertainties') if 'COUNTS' in f['SPECTRUM'].data.names: - stat_err = np.sqrt(f['SPECTRUM'].data['COUNTS']) / f['SPECTRUM'].header['EXPOSURE'] - elif 'RATE' in f['SPECTRUM'].data.names: - stat_err = np.sqrt(f['SPECTRUM'].data['RATE'] * f['SPECTRUM'].header['EXPOSURE']) / f['SPECTRUM'].header['EXPOSURE'] + stat_err = np.sqrt(f['SPECTRUM'].data['COUNTS']) / exposure else: stat_err = None else: @@ -106,12 +108,18 @@ def from_file_name_normal(fn): else: sys_err = None + if 'GROUPING' in f['SPECTRUM'].data.names: + grouping = f['SPECTRUM'].data['GROUPING'] + else: + grouping = None + return PHAI.from_arrays( - exposure=f['SPECTRUM'].header['EXPOSURE'], + exposure=exposure, rate=rate, stat_err=stat_err, sys_err=sys_err, - filename=fn + filename=fn, + grouping=grouping ) @staticmethod @@ -123,6 +131,7 @@ def from_arrays( quality=None, counts=None, filename=None, + grouping=None ): self = PHAI() @@ -159,6 +168,11 @@ def from_arrays( self._quality = quality else: self._quality = np.zeros_like(rate) + + if grouping is not None: + self._grouping = grouping + else: + self._grouping = np.zeros_like(rate) return self From ce397af9cca0e4fd1f5d833d091b3703b3456081 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:15:12 +0200 Subject: [PATCH 07/11] add test to read Poissonian uncertainties --- tests/test_spectra.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/test_spectra.py b/tests/test_spectra.py index 8af2aa8..bae1e16 100644 --- a/tests/test_spectra.py +++ b/tests/test_spectra.py @@ -171,3 +171,14 @@ def model_gen(p): # TODO: check that it all looks the same plt.savefig("unf_png.png") + +def test_read_poisson(): + import ogip.spec + pha = ogip.spec.PHAI.from_file_name("tests/data/MOS1source_spectrum_150_rbn.pi") + from astropy.io import fits as pf + import numpy as np + ff = pf.open("tests/data/MOS1source_spectrum_150_rbn.pi") + counts = ff[1].data['COUNTS'] + + assert(np.abs(np.sum(counts - (pha._rate)*pha._exposure)) < 1e-12) + assert(np.sum(pha._rate*pha._exposure - (pha._stat_err*pha._exposure)**2)) From 19e12da2c0bf0aa383bafb6de13f182ce259b9ce Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:16:05 +0200 Subject: [PATCH 08/11] missing save in notebook --- tools.ipynb | 23 +++-------------------- 1 file changed, 3 insertions(+), 20 deletions(-) diff --git a/tools.ipynb b/tools.ipynb index 093ffd0..e61ab1c 100644 --- a/tools.ipynb +++ b/tools.ipynb @@ -2,23 +2,9 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "ename": "Exception", - "evalue": "failed to read with any method, tried: {'from_file_name_osa': FileNotFoundError(2, 'No such file or directory'), 'from_file_name_normal': FileNotFoundError(2, 'No such file or directory')}", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mException\u001b[0m Traceback (most recent call last)", - "Cell \u001b[0;32mIn[1], line 6\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mogip\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mtools\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m plot, transform_rmf, crab_ph_cm2_s_kev, convolve\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mmatplotlib\u001b[39;00m\u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mpylab\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mplt\u001b[39;00m\n\u001b[0;32m----> 6\u001b[0m crab_pha \u001b[38;5;241m=\u001b[39m \u001b[43mogip\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mspec\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPHAI\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_file_name\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mcrab_pha.fits\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m 7\u001b[0m crab_rmf \u001b[38;5;241m=\u001b[39m ogip\u001b[38;5;241m.\u001b[39mspec\u001b[38;5;241m.\u001b[39mRMF\u001b[38;5;241m.\u001b[39mfrom_file_name(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcrab_rmf.fits\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m 8\u001b[0m crab_arf \u001b[38;5;241m=\u001b[39m ogip\u001b[38;5;241m.\u001b[39mspec\u001b[38;5;241m.\u001b[39mARF\u001b[38;5;241m.\u001b[39mfrom_file_name(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcrab_arf.fits\u001b[39m\u001b[38;5;124m\"\u001b[39m) \n", - "File \u001b[0;32m~/Soft/ogip/ogip/spec.py:62\u001b[0m, in \u001b[0;36mPHAI.from_file_name\u001b[0;34m(fn)\u001b[0m\n\u001b[1;32m 60\u001b[0m \u001b[38;5;129m@staticmethod\u001b[39m\n\u001b[1;32m 61\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mfrom_file_name\u001b[39m(fn):\n\u001b[0;32m---> 62\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mReadable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrom_file_name\u001b[49m\u001b[43m(\u001b[49m\u001b[43mPHAI\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfn\u001b[49m\u001b[43m)\u001b[49m\n", - "File \u001b[0;32m~/Soft/ogip/ogip/spec.py:47\u001b[0m, in \u001b[0;36mReadable.from_file_name\u001b[0;34m(cls, fn)\u001b[0m\n\u001b[1;32m 42\u001b[0m logger\u001b[38;5;241m.\u001b[39mdebug(\n\u001b[1;32m 43\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfailed to read with \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m.\u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m: \u001b[39m\u001b[38;5;132;01m%s\u001b[39;00m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mcls\u001b[39m\u001b[38;5;241m.\u001b[39m\u001b[38;5;18m__name__\u001b[39m, n, m, e\n\u001b[1;32m 44\u001b[0m )\n\u001b[1;32m 45\u001b[0m methods_tried[n] \u001b[38;5;241m=\u001b[39m e\n\u001b[0;32m---> 47\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mException\u001b[39;00m(\n\u001b[1;32m 48\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfailed to read with any method, tried: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mmethods_tried\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", - "\u001b[0;31mException\u001b[0m: failed to read with any method, tried: {'from_file_name_osa': FileNotFoundError(2, 'No such file or directory'), 'from_file_name_normal': FileNotFoundError(2, 'No such file or directory')}" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", "import ogip\n", @@ -116,13 +102,10 @@ } ], "metadata": { - "interpreter": { - "hash": "6c0c14c3bc85a7c6c1cb359b38aedade0f64d7248b8647c7b1bd974c3ac87f4c" - }, "kernelspec": { "display_name": "py10", "language": "python", - "name": "py10" + "name": "python3" }, "language_info": { "codemirror_mode": { From 5cdc06f3caa415d297036b50f7a292de96feacee Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:19:30 +0200 Subject: [PATCH 09/11] style change in test_spectra --- tests/test_spectra.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_spectra.py b/tests/test_spectra.py index bae1e16..79b4bc5 100644 --- a/tests/test_spectra.py +++ b/tests/test_spectra.py @@ -172,6 +172,7 @@ def model_gen(p): plt.savefig("unf_png.png") + def test_read_poisson(): import ogip.spec pha = ogip.spec.PHAI.from_file_name("tests/data/MOS1source_spectrum_150_rbn.pi") From d0c5ca98e910a655c1a8e163727bf640f3941973 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:28:29 +0200 Subject: [PATCH 10/11] more formatting --- tests/test_spectra.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_spectra.py b/tests/test_spectra.py index 79b4bc5..79c7a99 100644 --- a/tests/test_spectra.py +++ b/tests/test_spectra.py @@ -181,5 +181,5 @@ def test_read_poisson(): ff = pf.open("tests/data/MOS1source_spectrum_150_rbn.pi") counts = ff[1].data['COUNTS'] - assert(np.abs(np.sum(counts - (pha._rate)*pha._exposure)) < 1e-12) - assert(np.sum(pha._rate*pha._exposure - (pha._stat_err*pha._exposure)**2)) + assert (np.abs(np.sum(counts - (pha._rate) * pha._exposure)) < 1e-12) + assert (np.sum(pha._rate * pha._exposure - (pha._stat_err * pha._exposure)**2)) From 66adc7b51ab8172bbdc4b80a5822d2184e140229 Mon Sep 17 00:00:00 2001 From: Carlo Ferrigno Date: Thu, 5 Sep 2024 10:42:44 +0200 Subject: [PATCH 11/11] style again on spec.py --- ogip/spec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ogip/spec.py b/ogip/spec.py index c739aca..24fbe17 100644 --- a/ogip/spec.py +++ b/ogip/spec.py @@ -168,7 +168,7 @@ def from_arrays( self._quality = quality else: self._quality = np.zeros_like(rate) - + if grouping is not None: self._grouping = grouping else: