Commit a5b019ec authored by Marcel Strzys's avatar Marcel Strzys
Browse files

- Added the function for generating the background based on Yoshiki's scripts...

- Added the function for generating the background based on Yoshiki's scripts (note that due to a bug in gammapy the current output bkg hdu cannot be read with gammapy)
- changed 'type' to 'tel_type' in the TelescopeDescription since 'type' gave me errors
parent 7686b08f
...@@ -490,25 +490,92 @@ class IRFGenerator: ...@@ -490,25 +490,92 @@ class IRFGenerator:
return aeff_hdu return aeff_hdu
def _generate_background_hdu(self): def _generate_bkg_hdu(self):
# columns = [ bkg_shower_data = self.bkg_shower_data.query(self.cuts)
# ]
# colDefs = pyfits.ColDefs(columns) # Compute elapsed observation time
# bkg_hdu = pyfits.BinTableHDU.from_columns(colDefs) elapsed_time = np.array([])
# bkg_hdu.name = 'BACKGROUND' obs_id_list = np.array(bkg_shower_data.index.levels[0])
for obs_item in obs_id_list:
obs_item_events = bkg_shower_data.loc[(obs_item, slice(None), slice(None))]
obs_event_mean_arr_time = obs_item_events.groupby(['obs_id', 'event_id'])['mjd'].mean()
time_diff = np.diff(obs_event_mean_arr_time)*u.day.to(u.s)
# excludes gaps of possible technical problems
time_diff = time_diff[np.where(time_diff < 3e-1)]
elapsed_time = np.append(elapsed_time, np.sum(time_diff))
elapsed_time = np.sum(elapsed_time)
# Computing camera off-center angle for background events
offcenter = angular_separation(bkg_shower_data['az_reco_mean'].values * u.rad,
bkg_shower_data['alt_reco_mean'].values * u.rad,
bkg_shower_data['tel_az'].values * u.rad,
bkg_shower_data['tel_alt'].values * u.rad)
offcenter = offcenter.to(u.deg)
bkg_shower_data = bkg_shower_data.loc[slice(None), ['energy_reco_mean']]
bkg_shower_data['offcenter'] = offcenter
# Binning in energy
energy_edges = scipy.logspace(scipy.log10(self.min_energy),
scipy.log10(self.max_energy),
self.n_energy_bins+1)
energ_lo = energy_edges[:-1]
energ_hi = energy_edges[1:]
# Binning in off-center distance
theta_edges = scipy.linspace(self.min_theta,
self.max_theta,
self.n_theta_bins+1)
theta_lo = theta_edges[:-1]
theta_hi = theta_edges[1:]
# bkg_hdu.header['HDUDOC'] = 'https://github.com/open-gamma-ray-astro/gamma-astro-data-formats' bkg_event_matrix, _, _ = scipy.histogram2d(bkg_shower_data['energy_reco_mean'].values,
# bkg_hdu.header['HDUVERS'] = '0.2' bkg_shower_data['offcenter'].values,
# bkg_hdu.header['HDUCLASS'] = 'GADF' bins=[energy_edges, theta_edges])
# bkg_hdu.header['HDUCLAS1'] = 'RESPONSE'
# bkg_hdu.header['HDUCLAS2'] = 'BKG' # Compute bin sizes for density
# bkg_hdu.header['HDUCLAS3'] = 'FULL-ENCLOSURE' theta_area = np.pi * np.diff(theta_edges**2)
# bkg_hdu.header['HDUCLAS4'] = 'BKG_2D' energy_width = np.diff(energy_edges)
bkg_matrix = bkg_event_matrix / elapsed_time / theta_area / energy_width.reshape((-1, 1))
# --------------------------
# --- Converting to FITS ---
col_energ_lo = pyfits.Column(name='ENERG_LO', unit='TeV', format=f'{energ_lo.size}E', array=[energ_lo])
col_energ_hi = pyfits.Column(name='ENERG_HI', unit='TeV', format=f'{energ_hi.size}E', array=[energ_hi])
col_theta_lo = pyfits.Column(name='THETA_LO', unit='deg', format=f'{theta_lo.size}E', array=[theta_lo])
col_theta_hi = pyfits.Column(name='THETA_HI', unit='deg', format=f'{theta_hi.size}E', array=[theta_hi])
with pyfits.open('irf_file.fits') as fin: col_bkg_matrix = pyfits.Column(name='BKG', unit='s^-1 MeV^-1 sr^-1', format=f"{bkg_matrix.size}E",
bkg_hdu = fin['BACKGROUND'].copy() array=[bkg_matrix.transpose()],
dim=str(bkg_matrix.shape))
columns = [
col_energ_lo,
col_energ_hi,
col_theta_lo,
col_theta_hi,
col_bkg_matrix
]
# Aeff HDU
colDefs = pyfits.ColDefs(columns)
bkg_hdu = pyfits.BinTableHDU.from_columns(colDefs)
bkg_hdu.name = 'BACKGROUND'
bkg_hdu.header['HDUDOC'] = 'https://github.com/open-gamma-ray-astro/gamma-astro-data-formats'
bkg_hdu.header['HDUVERS'] = '0.2'
bkg_hdu.header['HDUCLASS'] = 'GADF'
bkg_hdu.header['HDUCLAS1'] = 'RESPONSE'
bkg_hdu.header['HDUCLAS2'] = 'BKG'
bkg_hdu.header['HDUCLAS3'] = 'FULL-ENCLOSURE'
bkg_hdu.header['HDUCLAS4'] = 'BKG_2D'
# --------------------------
return bkg_hdu return bkg_hdu
def generate_irf(self, output_name): def generate_irf(self, output_name):
...@@ -574,7 +641,7 @@ magic_tel_positions = { ...@@ -574,7 +641,7 @@ magic_tel_positions = {
magic_optics = OpticsDescription.from_name('MAGIC') magic_optics = OpticsDescription.from_name('MAGIC')
magic_cam = CameraGeometry.from_name('MAGICCam') magic_cam = CameraGeometry.from_name('MAGICCam')
magic_tel_description = TelescopeDescription(name='MAGIC', magic_tel_description = TelescopeDescription(name='MAGIC',
type='MAGIC', tel_type='MAGIC',
optics=magic_optics, optics=magic_optics,
camera=magic_cam) camera=magic_cam)
magic_tel_descriptions = {1: magic_tel_description, magic_tel_descriptions = {1: magic_tel_description,
...@@ -582,7 +649,8 @@ magic_tel_descriptions = {1: magic_tel_description, ...@@ -582,7 +649,8 @@ magic_tel_descriptions = {1: magic_tel_description,
# ----------------- # -----------------
mc_file_name = '../../../MCs/MAGIC/ST.03.07/za05to35/Test_sample/3.Reco/reco_m1.h5' mc_file_name = '../../../MCs/MAGIC/ST.03.07/za05to35/Test_sample/3.Reco/reco_m1.h5'
irf_generator = IRFGenerator(mc_file_name) bkg_file_name = '/remote/ceph/group/magic/MAGIC-LST/Data/MAGIC/Off/Test_sample/ctapipe/reco/reco_m1_magic_clean_step_20170109.h5'
irf_generator = IRFGenerator(mc_file_name, bkg_file_name)
irf_generator.set_energy_binning(min_energy=0.1, max_energy=30, n_energy_bins=10) irf_generator.set_energy_binning(min_energy=0.1, max_energy=30, n_energy_bins=10)
irf_generator.set_theta_binning(min_theta=0.0, max_theta=1.5, n_theta_bins=5) irf_generator.set_theta_binning(min_theta=0.0, max_theta=1.5, n_theta_bins=5)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment