From 1fd9ed99807189b6a32d8da003f96c51f9f88664 Mon Sep 17 00:00:00 2001 From: Israel Martinez Date: Tue, 26 Nov 2024 18:10:30 -0500 Subject: [PATCH] Automatically detect whether the response as PA info or not. Streamline the header for all polarization and sparse combinations --- cosipy/response/FullDetectorResponse.py | 256 ++++++------------------ 1 file changed, 66 insertions(+), 190 deletions(-) diff --git a/cosipy/response/FullDetectorResponse.py b/cosipy/response/FullDetectorResponse.py index db295903..0e2bd41b 100644 --- a/cosipy/response/FullDetectorResponse.py +++ b/cosipy/response/FullDetectorResponse.py @@ -74,9 +74,6 @@ def open(cls, filename,Spectrumfile=None,norm="Linear" ,single_pixel = False,alp emin,emax : float emin/emax used in the simulation source file. - - polarization : bool, - True if response includes polarization angle. """ filename = Path(filename) @@ -85,7 +82,7 @@ def open(cls, filename,Spectrumfile=None,norm="Linear" ,single_pixel = False,alp if filename.suffix == ".h5": return cls._open_h5(filename) elif "".join(filename.suffixes[-2:]) == ".rsp.gz": - return cls._open_rsp(filename,Spectrumfile,norm ,single_pixel,alpha,emin,emax, polarization) + return cls._open_rsp(filename,Spectrumfile,norm ,single_pixel,alpha,emin,emax) else: raise ValueError( "Unsupported file format. Only .h5 and .rsp.gz extensions are supported.") @@ -145,7 +142,7 @@ def _open_h5(cls, filename): return new @classmethod - def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = False,alpha=0,emin=90,emax=10000, polarization=False): + def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = False,alpha=0,emin=90,emax=10000): """ Open a detector response rsp file. @@ -171,16 +168,11 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal emin,emax : float emin/emax used in the simulation source file. - - polarization : bool, - True if response includes polarization angle. """ - if polarization == True: - labels = ("Ei", "NuLambda", "Pol", "Em", "Phi", "PsiChi", "SigmaTau", "Dist") - else: - labels = ("Ei", "NuLambda", "Em", "Phi", "PsiChi", "SigmaTau", "Dist") - + + + axes_names = [] axes_edges = [] axes_types = [] sparse = None @@ -223,6 +215,9 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal if line[1] == "false" : sparse = False + elif key == 'AN': + axes_names += [" ".join(line[1:])] + elif key == 'AD': if axes_types[-1] == "FISBEL": @@ -255,6 +250,16 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal elif key == "StartStream": nbins = int(line[1]) break + + # Check axes names and relabel + if np.array_equal(axes_names, ['"Initial energy [keV]"', '"#nu [deg]" "#lambda [deg]"', '"Polarization Angle [deg]"', '"Measured energy [keV]"', '"#phi [deg]"', '"#psi [deg]" "#chi [deg]"', '"#sigma [deg]" "#tau [deg]"', '"Distance [cm]"']): + has_polarization = True + labels = ("Ei", "NuLambda", "Pol", "Em", "Phi", "PsiChi", "SigmaTau", "Dist") + elif np.array_equal(axes_names, ['"Initial energy [keV]"', '"#nu [deg]" "#lambda [deg]"', '"Measured energy [keV]"', '"#phi [deg]"', '"#psi [deg]" "#chi [deg]"', '"#sigma [deg]" "#tau [deg]"', '"Distance [cm]"']): + has_polarization = False + labels = ("Ei", "NuLambda", "Pol", "Em", "Phi", "PsiChi", "SigmaTau", "Dist") + else: + raise InputError("Unknown response format") #check if the type of spectrum is known assert norm=="powerlaw" or norm=="Mono" or norm=="Linear" or norm=="Gaussian",f"unknown normalisation ! {norm}" @@ -522,196 +527,67 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal # Header drm.attrs['UNIT'] = 'cm2' - #sparse - if sparse : - drm.attrs['SPARSE'] = True - - # Axes - axes = drm.create_group('AXES', track_order=True) - - if polarization == True: - for axis in dr.axes[['NuLambda', 'Ei', 'Pol', 'Em', 'Phi', 'PsiChi', 'SigmaTau', 'Dist']]: - - axis_dataset = axes.create_dataset(axis.label, - data=axis.edges) - - - if axis.label in ['NuLambda', 'PsiChi','SigmaTau']: - - # HEALPix - axis_dataset.attrs['TYPE'] = 'healpix' - - axis_dataset.attrs['NSIDE'] = nside - - axis_dataset.attrs['SCHEME'] = 'ring' - - else: - - # 1D - axis_dataset.attrs['TYPE'] = axis.axis_scale - - if axis.label in ['Ei', 'Em']: - axis_dataset.attrs['UNIT'] = 'keV' - axis_dataset.attrs['TYPE'] = 'log' - elif axis.label in ['Phi', 'Pol']: - axis_dataset.attrs['UNIT'] = 'deg' - axis_dataset.attrs['TYPE'] = 'linear' - elif axis.label in ['Dist']: - axis_dataset.attrs['UNIT'] = 'cm' - axis_dataset.attrs['TYPE'] = 'linear' - else: - raise ValueError("Shouldn't happend") - - axis_description = {'Ei': "Initial simulated energy", - 'NuLambda': "Location of the simulated source in the spacecraft coordinates", - 'Pol': "Polarization angle", - 'Em': "Measured energy", - 'Phi': "Compton angle", - 'PsiChi': "Location in the Compton Data Space", - 'SigmaTau': "Electron recoil angle", - 'Dist': "Distance from first interaction" - } - - axis_dataset.attrs['DESCRIPTION'] = axis_description[axis.label] - else: - for axis in dr.axes[['NuLambda', 'Ei', 'Em', 'Phi', 'PsiChi','SigmaTau','Dist']]: - - axis_dataset = axes.create_dataset(axis.label, - data=axis.edges) - - - if axis.label in ['NuLambda', 'PsiChi','SigmaTau']: - - # HEALPix - axis_dataset.attrs['TYPE'] = 'healpix' - - axis_dataset.attrs['NSIDE'] = nside - - axis_dataset.attrs['SCHEME'] = 'ring' + axis_description = {'Ei': "Initial simulated energy", + 'NuLambda': "Location of the simulated source in the spacecraft coordinates", + 'Pol': "Polarization angle", + 'Em': "Measured energy", + 'Phi': "Compton angle", + 'PsiChi': "Location in the Compton Data Space", + 'SigmaTau': "Electron recoil angle", + 'Dist': "Distance from first interaction" + } + + #keep the same dimension order of the data + axes_to_write = ['NuLambda', 'Ei'] + + if has_polarization: + axes_to_write += ['Pol'] - else: + axes_to_write += ['Em', 'Phi', 'PsiChi'] - # 1D - axis_dataset.attrs['TYPE'] = axis.axis_scale - - if axis.label in ['Ei', 'Em']: - axis_dataset.attrs['UNIT'] = 'keV' - axis_dataset.attrs['TYPE'] = 'log' - elif axis.label in ['Phi']: - axis_dataset.attrs['UNIT'] = 'deg' - axis_dataset.attrs['TYPE'] = 'linear' - elif axis.label in ['Dist']: - axis_dataset.attrs['UNIT'] = 'cm' - axis_dataset.attrs['TYPE'] = 'linear' - else: - raise ValueError("Shouldn't happend") - - axis_description = {'Ei': "Initial simulated energy", - 'NuLambda': "Location of the simulated source in the spacecraft coordinates", - 'Em': "Measured energy", - 'Phi': "Compton angle", - 'PsiChi': "Location in the Compton Data Space", - 'SigmaTau': "Electron recoil angle", - 'Dist': "Distance from first interaction" - } - - axis_dataset.attrs['DESCRIPTION'] = axis_description[axis.label] - - #non sparse - else : - drm.attrs['SPARSE'] = False + if sparse: + drm.attrs['SPARSE'] = True + + # singletos. Save space in dense + axes_to_write += ['SigmaTau', 'Dist'] + else: + drm.attrs['SPARSE'] = False + + axes = drm.create_group('AXES', track_order=True) - # Axes - axes = drm.create_group('AXES', track_order=True) + for axis in dr.axes[axes_to_write]: - #keep the same dimension order of the data - if polarization == True: - for axis in dr.axes[['NuLambda','Ei', 'Pol', 'Em', 'Phi', 'PsiChi']]:#'SigmaTau','Dist']]: + axis_dataset = axes.create_dataset(axis.label, + data=axis.edges) - axis_dataset = axes.create_dataset(axis.label, - data=axis.edges) - - if axis.label in ['NuLambda', 'PsiChi']:#,'SigmaTau']: + if axis.label in ['NuLambda', 'PsiChi','SigmaTau']: - # HEALPix - axis_dataset.attrs['TYPE'] = 'healpix' + # HEALPix + axis_dataset.attrs['TYPE'] = 'healpix' - axis_dataset.attrs['NSIDE'] = nside + axis_dataset.attrs['NSIDE'] = nside - axis_dataset.attrs['SCHEME'] = 'ring' - - else: + axis_dataset.attrs['SCHEME'] = 'ring' - # 1D - axis_dataset.attrs['TYPE'] = axis.axis_scale - - if axis.label in ['Ei', 'Em']: - axis_dataset.attrs['UNIT'] = 'keV' - axis_dataset.attrs['TYPE'] = 'log' - elif axis.label in ['Phi', 'Pol']: - axis_dataset.attrs['UNIT'] = 'deg' - axis_dataset.attrs['TYPE'] = 'linear' - #elif axis.label in ['Dist']: - # axis_dataset.attrs['UNIT'] = 'cm' - # axis_dataset.attrs['TYPE'] = 'linear' - else: - raise ValueError("Shouldn't happend") - - axis_description = {'Ei': "Initial simulated energy", - 'NuLambda': "Location of the simulated source in the spacecraft coordinates", - 'Pol': "Polarization angle", - 'Em': "Measured energy", - 'Phi': "Compton angle", - 'PsiChi': "Location in the Compton Data Space", - #'SigmaTau': "Electron recoil angle", - #'Dist': "Distance from first interaction" - } - - axis_dataset.attrs['DESCRIPTION'] = axis_description[axis.label] else: - for axis in dr.axes[['NuLambda','Ei', 'Em', 'Phi', 'PsiChi']]:#'SigmaTau','Dist']]: - - axis_dataset = axes.create_dataset(axis.label, - data=axis.edges) - - - if axis.label in ['NuLambda', 'PsiChi']:#,'SigmaTau']: - # HEALPix - axis_dataset.attrs['TYPE'] = 'healpix' - - axis_dataset.attrs['NSIDE'] = nside - - axis_dataset.attrs['SCHEME'] = 'ring' - - else: + # 1D + axis_dataset.attrs['TYPE'] = axis.axis_scale + + if axis.label in ['Ei', 'Em']: + axis_dataset.attrs['UNIT'] = 'keV' + axis_dataset.attrs['TYPE'] = 'log' + elif axis.label in ['Phi', 'Pol']: + axis_dataset.attrs['UNIT'] = 'deg' + axis_dataset.attrs['TYPE'] = 'linear' + elif axis.label in ['Dist']: + axis_dataset.attrs['UNIT'] = 'cm' + axis_dataset.attrs['TYPE'] = 'linear' + else: + raise ValueError("Shouldn't happend") - # 1D - axis_dataset.attrs['TYPE'] = axis.axis_scale - - if axis.label in ['Ei', 'Em']: - axis_dataset.attrs['UNIT'] = 'keV' - axis_dataset.attrs['TYPE'] = 'log' - elif axis.label in ['Phi']: - axis_dataset.attrs['UNIT'] = 'deg' - axis_dataset.attrs['TYPE'] = 'linear' - #elif axis.label in ['Dist']: - # axis_dataset.attrs['UNIT'] = 'cm' - # axis_dataset.attrs['TYPE'] = 'linear' - else: - raise ValueError("Shouldn't happend") - - axis_description = {'Ei': "Initial simulated energy", - 'NuLambda': "Location of the simulated source in the spacecraft coordinates", - 'Em': "Measured energy", - 'Phi': "Compton angle", - 'PsiChi': "Location in the Compton Data Space", - #'SigmaTau': "Electron recoil angle", - #'Dist': "Distance from first interaction" - } - - axis_dataset.attrs['DESCRIPTION'] = axis_description[axis.label] + axis_dataset.attrs['DESCRIPTION'] = axis_description[axis.label] #sparse matrice if sparse :