From d71b071624b19cfb18e0c0a31796ad30f9895bc3 Mon Sep 17 00:00:00 2001 From: "dr.hualinxiao@gmail.com" Date: Tue, 13 Feb 2024 15:40:15 +0100 Subject: [PATCH] add compression error calculator --- setup.py | 4 +- stixdcpy/compression_error.py | 88 +++++++++++++++++++ stixdcpy/science.py | 159 ++++++++++++++++++++++++++++++++-- 3 files changed, 240 insertions(+), 11 deletions(-) create mode 100644 stixdcpy/compression_error.py diff --git a/setup.py b/setup.py index a968797..1acd77a 100755 --- a/setup.py +++ b/setup.py @@ -3,11 +3,11 @@ setup( name='stixdcpy', description='STIX data center APIs and data analysis tools', - version='2.2', + version='3.0', author='Hualin Xiao', author_email='hualin.xiao@fhnw.ch', long_description=open('README.md').read(), - install_requires=['numpy', 'requests', 'python-dateutil', + install_requires=['numpy', 'requests', 'python-dateutil', 'ipython', 'astropy', 'matplotlib','tqdm','pandas','joblib', 'roentgen', 'simplejson', 'sunpy','wget'], long_description_content_type='text/markdown', #packages=find_packages(where='stixdcpy'), diff --git a/stixdcpy/compression_error.py b/stixdcpy/compression_error.py new file mode 100644 index 0000000..61a6211 --- /dev/null +++ b/stixdcpy/compression_error.py @@ -0,0 +1,88 @@ +""" +This script provides functions for decompressing integer values and creating error lookup tables. + +It includes the following functions: + +- decompress(x, K, M): Decompresses a compressed integer value. +- make_lut(k, m): Creates a lookup table for error calculation based on k and m values. + +Additionally, it defines error lookup tables for specific combinations of k and m values. + +Author: Hualin Xiao(hualin.xiao@fhnw.ch) + +""" + +import numpy as np + +MAX_STORED_INTEGER = 1e8 + +def decompress(x, K, M): + """ + Decompresses a compressed integer value. + + Parameters: + x (int): The compressed integer value to decompress. + K (int): The number of bits reserved for the exponent. + M (int): The number of bits reserved for the mantissa. + + Returns: + tuple: A tuple containing the error and the decompressed value. + The error represents the uncertainty in the decompressed value. + The decompressed value is the original integer value. + """ + S = 0 + if S + K + M > 8 or S not in (0, 1) or K > 7 or M > 7: + return None, None + if K == 0 or M == 0: + return None, None + + sign = 1 + if S == 1: # signed + MSB = x & (1 << 7) + if MSB != 0: + sign = -1 + x = x & ((1 << 7) - 1) + + x0 = 1 << (M + 1) + if x < x0: + return None, x + mask1 = (1 << M) - 1 + mask2 = (1 << M) + mantissa1 = x & mask1 + exponent = (x >> M) - 1 + # number of shifted bits + mantissa2 = mask2 | mantissa1 # add 1 before mantissa + low = mantissa2 << exponent # minimal possible value + high = low | ((1 << exponent) - 1) # maximal possible value + mean = (low + high) >> 1 # mean value + error = np.sqrt((high - low) ** 2 / 12) + + if mean > MAX_STORED_INTEGER: + return error, float(mean) + + return error, sign * mean + +def make_lut(k, m): + """ + Creates a lookup table for error calculation based on k and m values. + + Parameters: + k (int): The number of bits reserved for the exponent. + m (int): The number of bits reserved for the mantissa. + + Returns: + dict: A dictionary mapping decompressed values to their respective errors. + """ + res = {} + for i in range(256): + err, val = decompress(i, k, m) + if err is not None: + res[val] = np.sqrt(err ** 2 + val) + return res + +error_luts = {'035': make_lut(3, 5), '044': make_lut(4, 4), '053': make_lut(5, 3)} + +if __name__ == '__main__': + from pprint import pprint + pprint(error_luts) + diff --git a/stixdcpy/science.py b/stixdcpy/science.py index 3ee020c..c0413d4 100644 --- a/stixdcpy/science.py +++ b/stixdcpy/science.py @@ -281,9 +281,15 @@ def make_spectra(self, pixel_counts=None): @property def pixel_counts_error(self): - return error_computation(self.hdul['data'].data['counts_err'], + try: + counts_err=self.hdul['data'].data['counts_err'] + except KeyError: + counts_err=self.hdul['data'].data['counts_comp_err'] + return error_computation(counts_err, self.pixel_counts) + + def correct_dead_time(self): """ dead time correction Returns: @@ -295,22 +301,30 @@ def correct_dead_time(self): live_ratio: np.array """ - def correct(triggers, counts_arr, time_bins): + def correct(triggers, counts_arr, counts_err_arr, time_bins): """ Live time correction Args triggers: ndarray triggers in the spectrogram counts_arr:ndarray counts in the spectrogram + counts_err_arr:ndarray + counts error in the spectrogram time_bins: ndarray time_bins in the spectrogram Returns live_time_ratio: ndarray live time ratio of detectors - count_rate: + corrected_rates: corrected count rate + count_rate: + count rate before dead time correction + time: + timestamp photons_in: rate of photons illuminating the detector group + corrected_counts: + dead-time corrected counts """ time_bins = time_bins[:, None] @@ -321,31 +335,154 @@ def correct(triggers, counts_arr, time_bins): time_bins = time_bins[:, :, None, None] count_rate = counts_arr / time_bins + count_rate_err = counts_err_arr / time_bins # print(counts_arr.shape) for det in range(32): trig_idx = inst.detector_id_to_trigger_index(det) nin = photons_in[:, trig_idx] live_ratio[:, det] = np.exp( -BETA * nin * ASIC_TAU) / (1 + nin * TRIG_TAU) + + + + #live_ratio=np.zeros_like(live_ratio)+1 + #disable live time correction + corrected_rate = count_rate / live_ratio[:, :, None, None] + + corrected_rate_err = count_rate_err / live_ratio[:, :, None, None] + + corrected_counts = corrected_rate * time_bins + corrected_counts_err = corrected_rate_err * time_bins + return { 'corrected_rates': corrected_rate, + 'corrected_rate_err': corrected_rate_err, 'count_rate': count_rate, 'photons_in': photons_in, + 'corrected_counts':corrected_counts, + 'corrected_counts_err':corrected_counts_err, + 'time': self.datetime, + 'time_bins': time_bins.flatten(), 'live_ratio': live_ratio } - self.corrected = correct(self.triggers, self.pixel_counts, + self.corrected = correct(self.triggers, self.pixel_counts, self.pixel_counts_error, self.timedel) # approximate the live ratio error like Ewan does - above = correct(self.triggers + self.trigger_error, self.pixel_counts, + above = correct(self.triggers + self.trigger_error, self.pixel_counts,self.pixel_counts_error, self.timedel) - below = correct(self.triggers - self.trigger_error, self.pixel_counts, + below = correct(self.triggers - self.trigger_error, self.pixel_counts,self.pixel_counts_error, self.timedel) self.corrected['live_error'] = np.abs(above['live_ratio'] - below['live_ratio']) / 2 return self.corrected + def get_sum_counts(self, start_utc=None, end_utc=None) : + """ + Calculate the total counts in different regions of a pixel data file within a specified time range. + + Parameters: + - start_utc (str or None, optional): The start time in UTC format. If not provided, the beginning of the observation is used. + - end_utc (str or None, optional): The end time in UTC format. If not provided, the end of the observation is used. + + Returns: + - dict: A dictionary containing total counts in different regions: + - 'top': Total counts in the top row(channels 1-4). + - 'bottom': Total counts in the bottom row(channels 5-8). + - 'small': Total counts in the small pixels. + - 'total': Total counts, sum of 'top' and 'bottom'. + Note: + - The function relies on external modules 'sci' and 'st' for pixel data manipulation and time conversions. + - The time range for counts calculation is determined by start_utc and end_utc. + If these are not provided, the entire observation duration is considered. + """ + cl1=self.correct_dead_time() + data_start=self.T0_unix + data_end=self.T0_unix+self.duration + pixel_counts=cl1['corrected_counts'] + pixel_counts_err=cl1['corrected_counts_err'] + + if start_utc is None and end_utc is None: + start_i_tbin, end_i_tbin=0, pixel_counts.shape[0]-1 + duration = self.duration + + else: + start_unix=sdt.utc2unix(start_utc) + end_unix=sdt.utc2unix(end_utc) + if start_unix > data_end or end_unix < data_start: + return None + + start_unix=max(data_start, start_unix) + end_unix=min(data_end, end_unix) + + + + #duration=end_unix-start_unix + + start_time = start_unix - self.T0_unix #relative start time + end_time = end_unix - self.T0_unix + + start_i_tbin=np.argmax( + self.time - 0.5 * self.timedel >= start_time) if ( + 0 <= start_time <= self.duration) else 0 + end_i_tbin=np.argmin( + self.time + 0.5 * self.timedel <= end_time) if ( + start_time <= end_time <= self.duration)+1 else len( + self.time) + + duration=np.sum(self.timedel[start_i_tbin:end_i_tbin]) + + + + sum_counts = {'top': np.sum(pixel_counts[start_i_tbin:end_i_tbin, :,0:4,:], axis=(0,1,2) ), + 'bottom': np.sum(pixel_counts[start_i_tbin:end_i_tbin, :,4:8,:],axis=(0,1,2) ), + 'small': np.sum(pixel_counts[start_i_tbin:end_i_tbin, :,8:,:],axis=(0,1,2) ), + 'duration':duration, + + 'top_err': np.sqrt(np.sum(pixel_counts_err[start_i_tbin:end_i_tbin, :,0:4,:]**2, axis=(0,1,2) )), + 'bottom_err': np.sqrt(np.sum(pixel_counts_err[start_i_tbin:end_i_tbin, :,4:8,:]**2,axis=(0,1,2) )), + 'small_err': np.sqrt(np.sum(pixel_counts_err[start_i_tbin:end_i_tbin, :,8:,:]**2,axis=(0,1,2) )) + } + + sum_counts['big'] = sum_counts['top']+sum_counts['bottom'] + sum_counts['big_err'] = np.sqrt(sum_counts['top_err']**2 + + sum_counts['bottom_err']**2 ) + + sum_counts['total'] = sum_counts['top']+sum_counts['bottom']+sum_counts['small'] + sum_counts['total_err'] = np.sqrt(sum_counts['top_err']**2 + + sum_counts['bottom_err']**2 + sum_counts['small_err']**2) + + return sum_counts + + + def get_mean_rate(self,start_utc =None, end_utc=None): + """ + Calculate the mean count rate in different regions of a pixel data file within a specified time range. + + Parameters: + - start_utc (str or None, optional): The start time in UTC format. If not provided, the beginning of the observation is used. + - end_utc (str or None, optional): The end time in UTC format. If not provided, the end of the observation is used. + + Returns: + - dict: A dictionary containing total counts in different regions: + - 'top': mean count rate in the top row(channels 1-4). + - 'bottom': mean count rate in the bottom row(channels 5-8). + - 'small': mean count rate in the small pixels. + - 'total': mean rate of all pixels 'top' and 'bottom'. + Note: + - The function relies on external modules 'sci' and 'st' for pixel data manipulation and time conversions. + - The time range for counts calculation is determined by start_utc and end_utc. If these are not provided, the entire observation duration is considered. + """ + + cnts=get_sum_counts(start_utc , end_utc) + result={} + for key,val in cnts.items(): + norm = cnts['duration'] if key !='duration' else 1 + result[key]=val/norm + return result + + def peek(self, plots=['spg', 'lc', 'spec', 'tbin', 'qllc'], @@ -450,7 +587,7 @@ def __init__(self, l1sig: PixelData, l1bkg: PixelData): """ self.l1sig = l1sig self.l1bkg = l1bkg - print(self.l1sig.energy_bin_mask) + #print(self.l1sig.energy_bin_mask) dmask = self.l1bkg.energy_bin_mask - self.l1sig.energy_bin_mask if np.any(dmask < 0): @@ -465,7 +602,7 @@ def __init__(self, l1sig: PixelData, l1bkg: PixelData): ]) # set counts beyond the signal energy range to 0 self.subtracted_counts = (self.l1sig.counts - self.pixel_bkg_counts) - print(self.l1sig.inversed_energy_bin_mask) + #print(self.l1sig.inversed_energy_bin_mask) self.subtracted_counts *= self.l1sig.energy_bin_mask # Dead time correction needs to be included in the future @@ -660,7 +797,7 @@ def correct(triggers, counts_arr, time_bins, num_detectors): except KeyError: num_detectors = self.hdul[1].data['detector_mask'].sum() - self.corrected = correct(self.triggers, self.counts, self.timedel, + self.corrected = correct(self.triggers, self.counts, self.timedel, num_detectors) # approximate the live ratio error like Ewan does @@ -1005,3 +1142,7 @@ def spec_fits_concatenate(fitsfile1, hdul.writeto(outfilename, overwrite=True) return outfilename + + + +