Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Suggested changes for PR240 #1

Merged
merged 2 commits into from
Dec 16, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
256 changes: 66 additions & 190 deletions cosipy/response/FullDetectorResponse.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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.")
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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":
Expand Down Expand Up @@ -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}"
Expand Down Expand Up @@ -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 :
Expand Down
4 changes: 2 additions & 2 deletions cosipy/test_data/polarization_ori.ori
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
Type OrientationsGalactic
OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895
OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895
OG 0.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0
OG 15.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0
EN