Skip to content

Commit

Permalink
Automatically detect whether the response as PA info or not. Streamli…
Browse files Browse the repository at this point in the history 
…ne the header for all polarization and sparse combinations
  • Loading branch information
@israelmcmc
israelmcmc committed Nov 26, 2024
1 parent 715df97 commit 1fd9ed9
Showing 1 changed file with 66 additions and 190 deletions.
256 changes: 66 additions & 190 deletions cosipy/response/FullDetectorResponse.py
View file
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

0 comments on commit 1fd9ed9

Please sign in to comment.