diff --git a/acceptance_modelisation.egg-info/PKG-INFO b/acceptance_modelisation.egg-info/PKG-INFO new file mode 100644 index 0000000..d020c72 --- /dev/null +++ b/acceptance_modelisation.egg-info/PKG-INFO @@ -0,0 +1,273 @@ +Metadata-Version: 2.1 +Name: acceptance-modelisation +Version: 0.2.0 +Summary: Calculate 2D and 3d acceptance model for creating maps with gammapy (IACT analysis) +Author: Mathieu de Bony de Lavergne +Author-email: mathieu.debony@cea.fr +Description-Content-Type: text/markdown + +# Description + +This package create acceptance model to be used for IACT analysis with gammapy + +# Installation + + +```bash +git clone https://github.com/mdebony/acceptance_modelisation.git +cd acceptance_modelisation +python setup.py install +``` + +Dependencies : +- numpy +- scipy +- astropy +- gammapy 1.1 +- regions 0.7 + +# Example of use + +## Basic use + +You could first create the acceptance model + +```python +from gammapy.maps import MapAxis +from gammapy.data import DataStore +from regions import CircleSkyRegion +import astropy.units as u +import numpy as np +from astropy.coordinates import SkyCoord +from acceptance_modelisation import RadialAcceptanceMapCreator + +# The observations to use for creating the acceptance model +data_store = DataStore.from_dir("$GAMMAPY_DATA/hess-dl3-dr1") +obs_collection = data_store.get_observations([23523, 23526, 23559, 23592]) + +# The exclusion regions to apply during acceptance model calculation +exclude_regions = [CircleSkyRegion(center=SkyCoord.from_name('Crab'), + radius=0.2 * u.deg), ] + +# Define the binning of the model +e_min, e_max = 0.1 * u.TeV, 10. * u.TeV +size_fov = 2.5 * u.deg +offset_axis_acceptance = MapAxis.from_bounds(0. * u.deg, size_fov, nbin=6, name='offset') +energy_axis_acceptance = MapAxis.from_energy_bounds(e_min, e_max, nbin=6, name='energy') + +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10) +acceptance_model = acceptance_model_creator.create_acceptance_map(obs_collection) + +``` + +You can then check the acceptance model by plotting it using + +```python +acceptance_model.peek() +``` + +To use it with gammapy, you could first save it on a FITS file + +```python +hdu_acceptance = acceptance_model.to_table_hdu() +hdu_acceptance.writeto('acceptance.fits', overwrite=True) +``` + +It's then possible to load the acceptance model in the current gammapy DataStore with this code. +You would need then to recreate you gammapy Observations object in order than the acceptance model is taken into account for the analysis. + +```python +data_store.hdu_table.remove_rows(data_store.hdu_table['HDU_TYPE']=='bkg') +for obs_id in np.unique(data_store.hdu_table['OBS_ID']): + data_store.hdu_table.add_row({'OBS_ID': obs_id, + 'HDU_TYPE': 'bkg', + "HDU_CLASS": "bkg_2d", + "FILE_DIR": "", + "FILE_NAME": 'acceptance.fits', + "HDU_NAME": "BACKGROUND", + "SIZE": hdu_acceptance.size}) +data_store.hdu_table = data_store.hdu_table.copy() + +obs_collection = data_store.get_observations([23523, 23526, 23559, 23592]) + +data_store.hdu_table +``` + +## Telescope position + +**The observations should contain the telescope position in order to have the algorithm working.** +If the information is missing in the DL3, you could either add it or it possible to add it directly to the observation as shown in the example below. +```python +from astropy.coordinates import EarthLocation + +# Your telescope position in an EarthLocation object +loc = EarthLocation.of_site('Roque de los Muchachos') + +# Add telescope position to observations +for i in obs_collection: + obs_collection[i].obs_info['GEOLON'] = loc.lon.value + obs_collection[i].obs_info['GEOLAT'] = loc.lat.value + obs_collection[i].obs_info['GEOALT'] = loc.height.value + obs_collection[i]._location = loc +``` + +## Runwise norm of the model + +It's also possible to fit the normalisation of the model per run. For this use the method create_acceptance_map_per_observation . +In that case the output is a dictionary containing the acceptance model of each observations (with the observation Id as index). +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection) +``` + +## Zenith binned model + +It's also possible to create model binned per cos zenith. For this use the method create_acceptance_map_per_observation but with the option `zenith_binning` set at True. +The width of zenith bin could be control at the creation of the object with the parameter `initial_cos_zenith_binning`. +You can chose the `cos_zenith_binning_method`: `min_livetime` or `min_n_observation` to give a condition on minimum livetime or number of observations for each bin. The algorithm will then automatically rebin to larger bin in order to have in each bin at least `cos_zenith_binning_parameter_value` livetime (in seconds) or observation per bin. If you add `_per_wobble` to the method name, wobbles will be identified and the binning will require the condition to be fulfilled for each identified wobble. +In that case the output is a dictionary containing the acceptance model of each observations (with the observation ID as index). +Set `verbose` to True to get the binning result and a plot of the binned livetime. +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10, + cos_zenith_binning_method='min_livetime_per_wobble', + cos_zenith_binning_parameter_value=3600, + initial_cos_zenith_binning=0.01, + verbose=True) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection, + zenith_binning=True) +``` + +## Zenith interpolated model + +It's also possible to create model interpolated between the binned model per cos zenith. For this use the method create_acceptance_map_per_observation but with the option `zenith_interpolation` set at True. All the parameters controlling the cos zenith binning remain active. +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10, + min_run_per_cos_zenith_bin=3, + initial_cos_zenith_binning=0.01) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection, + zenith_binning=True, + zenith_interpolation=True) +``` + +## Using OFF runs for background model + +It is also possible to create a model from OFF runs and to apply in ON runs you want to analyse. The exclusions regions should cover potential sources both in the ON and OFF runs at the same time. The OFF runs don't need to be spatially connected. +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10, + min_run_per_cos_zenith_bin=3, + initial_cos_zenith_binning=0.01) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection, + off_observations=obs_collection_off, + zenith_binning=True, + zenith_interpolation=True) +``` + +## Store background model for later application + +It could be in some case useful to precompute a model and to apply it on data later. The example below cover the case where you want to use a model per run using zenith interpolation and OFF runs. + +However, it's possible to use this functionality without OFF runs or zenith interpolation. In the last case you just need to provide a model in a gammapy format. + +There are at this stage no define file format for storing the intermediate results, we suggest to store the BackgroundCollectionZenith object created directly using pickle. + +### Creating the model + +The example below creates the `BackgroundCollectionZenith` object containing the zenith binned model. +The `obs_collection` provided could either be your ON runs if you want to compute background directly from the data or OFF runs if you want to use other data for the background model. + +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10, + min_run_per_cos_zenith_bin=3, + initial_cos_zenith_binning=0.01) +base_model = acceptance_model_creator.create_model_cos_zenith_binned(obs_collection) +``` + +### Storing and loading a model with pickle + +The `BackgroundCollectionZenith` object containing the models could then be store using pickle. +```python +import pickle +with open('my_bkg_model.pck', mode='wb') as f: + pickle.dump(base_model, f) +``` + +It is then possible to retrieve it later using pickle. +```python +import pickle +with open('my_bkg_model.pck', mode='rb') as f: + base_model = pickle.load(f) +``` + +### Applying a model in memory + +When you have a precomputed model in memory, it is possible to apply it directly on a given set of runs by using the `base_mode` parameter. +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection, + base_model=base_model, + zenith_binning=True, + zenith_interpolation=True) +``` + +## Compute background model with a higher time resolution than the observation run + +If the background evolve quickly, like at high zenith angle, you could compute in the case of zenith binned or interpolated background the model at a smaller time scale than the observation run. The background model for the run will then correspond to the average of all the model computed for each part of the run. +It should improve accuracy of the model at the expanse of a larger compute time. +For this you need to set `use_mini_irf_computation = True` and you could control the time resolution used for computation with the parameter `mini_irf_time_resolution`. + +```python +acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions, + oversample_map=10, + min_run_per_cos_zenith_bin=3, + initial_cos_zenith_binning=0.01, + use_mini_irf_computation = True, + mini_irf_time_resolution = 1. * u.min) +acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection, + zenith_binning=True, + zenith_interpolation=True) +``` + +# Available model + +All models have an identical interface. You just need to change the class used to change the model created. + +There are two model currently available : + +- A 2D model with hypothesis of a radial symmetry of the background across the FoV. This is the class `RadialAcceptanceMapCreator`. + ````python + from acceptance_modelisation import RadialAcceptanceMapCreator + acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions) + acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection) + ```` +- A 3D model with a regular grid describing the FoV. This is the class `Grid3DAcceptanceMapCreator`. + ````python + from acceptance_modelisation import Grid3DAcceptanceMapCreator + acceptance_model_creator = Grid3DAcceptanceMapCreator(energy_axis_acceptance, + offset_axis_acceptance, + exclude_regions=exclude_regions) + acceptance_models = acceptance_model_creator.create_acceptance_map_per_observation(obs_collection) + ```` diff --git a/acceptance_modelisation.egg-info/SOURCES.txt b/acceptance_modelisation.egg-info/SOURCES.txt new file mode 100644 index 0000000..94c0efa --- /dev/null +++ b/acceptance_modelisation.egg-info/SOURCES.txt @@ -0,0 +1,16 @@ +.gitignore +README.md +setup.py +acceptance_modelisation/__init__.py +acceptance_modelisation/base_acceptance_map_creator.py +acceptance_modelisation/bkg_collection.py +acceptance_modelisation/exception.py +acceptance_modelisation/grid3d_acceptance_map_creator.py +acceptance_modelisation/modeling.py +acceptance_modelisation/radial_acceptance_map_creator.py +acceptance_modelisation/toolbox.py +acceptance_modelisation.egg-info/PKG-INFO +acceptance_modelisation.egg-info/SOURCES.txt +acceptance_modelisation.egg-info/dependency_links.txt +acceptance_modelisation.egg-info/requires.txt +acceptance_modelisation.egg-info/top_level.txt \ No newline at end of file diff --git a/acceptance_modelisation.egg-info/dependency_links.txt b/acceptance_modelisation.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/acceptance_modelisation.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/acceptance_modelisation.egg-info/requires.txt b/acceptance_modelisation.egg-info/requires.txt new file mode 100644 index 0000000..cccac70 --- /dev/null +++ b/acceptance_modelisation.egg-info/requires.txt @@ -0,0 +1,5 @@ +gammapy==1.2 +numpy +scipy +astropy>=4.0 +regions<0.8,>=0.7 diff --git a/acceptance_modelisation.egg-info/top_level.txt b/acceptance_modelisation.egg-info/top_level.txt new file mode 100644 index 0000000..f040c68 --- /dev/null +++ b/acceptance_modelisation.egg-info/top_level.txt @@ -0,0 +1 @@ +acceptance_modelisation diff --git a/acceptance_modelisation/__init__.py b/acceptance_modelisation/__init__.py index 2953564..557c3eb 100644 --- a/acceptance_modelisation/__init__.py +++ b/acceptance_modelisation/__init__.py @@ -1,8 +1,8 @@ -from .base_radial_acceptance_map_creator import BaseRadialAcceptanceMapCreator from .base_acceptance_map_creator import BaseAcceptanceMapCreator from .grid3d_acceptance_map_creator import Grid3DAcceptanceMapCreator from .radial_acceptance_map_creator import RadialAcceptanceMapCreator from .bkg_collection import BackgroundCollectionZenith import logging + logging.basicConfig() diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 0841d59..a95b6cf 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -5,7 +5,7 @@ import astropy.units as u import numpy as np -from astropy.coordinates import SkyCoord, AltAz, SkyOffsetFrame +from astropy.coordinates import SkyCoord, AltAz, SkyOffsetFrame, Angle from astropy.coordinates.erfa_astrom import erfa_astrom, ErfaAstromInterpolator from astropy.time import Time from gammapy.data import Observations, Observation @@ -20,35 +20,39 @@ from .bkg_collection import BackgroundCollectionZenith from .exception import BackgroundModelFormatException -from .toolbox import (compute_rotation_speed_fov, - get_unique_wobble_pointings, - get_time_mini_irf, - generate_irf_from_mini_irf, - compute_neighbour_condition_validation) +from .toolbox import ( + compute_rotation_speed_fov, + get_unique_wobble_pointings, + get_time_mini_irf, + generate_irf_from_mini_irf, + compute_neighbour_condition_validation, +) +from sklearn.decomposition import PCA logger = logging.getLogger(__name__) class BaseAcceptanceMapCreator(ABC): - - def __init__(self, - energy_axis: MapAxis, - max_offset: u.Quantity, - spatial_resolution: u.Quantity, - exclude_regions: Optional[List[SkyRegion]] = None, - cos_zenith_binning_method: str = 'min_livetime', - cos_zenith_binning_parameter_value: int = 3600, - initial_cos_zenith_binning: float = 0.01, - max_angular_separation_wobble: u.Quantity = 0.4 * u.deg, - zenith_binning_run_splitting: bool = False, - max_fraction_pixel_rotation_fov: float = 0.5, - time_resolution: u.Quantity = 0.1 * u.s, - use_mini_irf_computation: bool = False, - mini_irf_time_resolution: u.Quantity = 1. * u.min, - interpolation_type: str = 'linear', - activate_interpolation_cleaning: bool = False, - interpolation_cleaning_energy_relative_threshold: float = 1e-4, - interpolation_cleaning_spatial_relative_threshold: float = 1e-2) -> None: + def __init__( + self, + energy_axis: MapAxis, + max_offset: u.Quantity, + spatial_resolution: u.Quantity, + exclude_regions: Optional[List[SkyRegion]] = None, + cos_zenith_binning_method: str = "min_livetime", + cos_zenith_binning_parameter_value: int = 3600, + initial_cos_zenith_binning: float = 0.01, + max_angular_separation_wobble: u.Quantity = 0.4 * u.deg, + zenith_binning_run_splitting: bool = False, + max_fraction_pixel_rotation_fov: float = 0.5, + time_resolution: u.Quantity = 0.1 * u.s, + use_mini_irf_computation: bool = False, + mini_irf_time_resolution: u.Quantity = 1.0 * u.min, + interpolation_type: str = "linear", + activate_interpolation_cleaning: bool = False, + interpolation_cleaning_energy_relative_threshold: float = 1e-4, + interpolation_cleaning_spatial_relative_threshold: float = 1e-2, + ) -> None: """ Create the class for calculating radial acceptance model. @@ -103,14 +107,23 @@ def __init__(self, self.exclude_regions = exclude_regions # Calculate map parameter - self.n_bins_map = 2 * int(np.rint((self.max_offset / spatial_resolution).to(u.dimensionless_unscaled))) + self.n_bins_map = 2 * int( + np.rint((self.max_offset / spatial_resolution).to(u.dimensionless_unscaled)) + ) self.spatial_bin_size = self.max_offset / (self.n_bins_map / 2) - self.center_map = SkyCoord(ra=0. * u.deg, dec=0. * u.deg, frame='icrs') - self.geom = WcsGeom.create(skydir=self.center_map, npix=(self.n_bins_map, self.n_bins_map), - binsz=self.spatial_bin_size, frame="icrs", axes=[self.energy_axis]) + self.center_map = SkyCoord(ra=0.0 * u.deg, dec=0.0 * u.deg, frame="icrs") + self.geom = WcsGeom.create( + skydir=self.center_map, + npix=(self.n_bins_map, self.n_bins_map), + binsz=self.spatial_bin_size, + frame="icrs", + axes=[self.energy_axis], + ) logger.info( - 'Computation will be made with a bin size of {:.3f} arcmin'.format( - self.spatial_bin_size.to_value(u.arcmin))) + "Computation will be made with a bin size of {:.3f} arcmin".format( + self.spatial_bin_size.to_value(u.arcmin) + ) + ) # Store computation parameters for run splitting self.max_fraction_pixel_rotation_fov = max_fraction_pixel_rotation_fov @@ -129,8 +142,12 @@ def __init__(self, # Store cleaning parameters for models created from interpolation self.activate_interpolation_cleaning = activate_interpolation_cleaning - self.interpolation_cleaning_energy_relative_threshold = interpolation_cleaning_energy_relative_threshold - self.interpolation_cleaning_spatial_relative_threshold = interpolation_cleaning_spatial_relative_threshold + self.interpolation_cleaning_energy_relative_threshold = ( + interpolation_cleaning_energy_relative_threshold + ) + self.interpolation_cleaning_spatial_relative_threshold = ( + interpolation_cleaning_spatial_relative_threshold + ) self.max_cleaning_iteration = 50 # Store mini irf computation parameters @@ -138,7 +155,11 @@ def __init__(self, self.mini_irf_time_resolution = mini_irf_time_resolution @staticmethod - def _transform_obs_to_camera_frame(obs: Observation) -> Observation: + def _transform_obs_to_camera_frame( + obs: Observation, + rotate_to_obs: Optional[Observation] = None, + derotate: Optional[bool] = False, + ) -> Observation: """ Transform events and pointing of an obs from a sky frame to camera frame @@ -146,6 +167,12 @@ def _transform_obs_to_camera_frame(obs: Observation) -> Observation: ---------- obs : gammapy.data.observations.Observation The observation to transform + rotate_to_obs : gammapy.data.observations.Observation + All events of `obs` are rotate in camera coordinates to match + the azimuth angle of `rotate_to_obs`. + For stereoscopic telescopes, like MAGIC, the morphology changes with azimuth. + Therefore the statistics can be enhacned with this instead of creating a background + model for multiple azimuth bins. Returns ------- @@ -153,37 +180,149 @@ def _transform_obs_to_camera_frame(obs: Observation) -> Observation: The observation transformed for reference in camera frame """ + def magic_theory_angle(az: Angle): + return az - 34.23 * u.deg + + def gamma_pca_corr_magic(angle: Angle, az: Angle): + gamma_theory = magic_theory_angle(az) + if gamma_theory <= 90 * u.deg and gamma_theory >= 0 * u.deg: + return angle + if gamma_theory <= 270 * u.deg and gamma_theory > 90 * u.deg: + return angle + 180 * u.deg + if gamma_theory <= 360 * u.deg and gamma_theory > 270 * u.deg: + return angle + 360 * u.deg + # Transform to altaz frame - altaz_frame = AltAz(obstime=obs.events.time, - location=obs.observatory_earth_location) + altaz_frame = AltAz(obstime=obs.events.time, location=obs.meta.location) events_altaz = obs.events.radec.transform_to(altaz_frame) - pointing_altaz = obs.get_pointing_icrs(obs.events.time).transform_to(altaz_frame) + pointing_altaz = obs.get_pointing_icrs(obs.events.time).transform_to( + altaz_frame + ) # Rotation to transform to camera frame - camera_frame = SkyOffsetFrame(origin=AltAz(alt=pointing_altaz.alt, - az=pointing_altaz.az, - obstime=obs.events.time, - location=obs.observatory_earth_location), - rotation=[0., ] * len(obs.events.time) * u.deg) + camera_frame = SkyOffsetFrame( + origin=AltAz( + alt=pointing_altaz.alt, + az=pointing_altaz.az, + obstime=obs.events.time, + location=obs.meta.location, + ), + rotation=[ + 0.0, + ] + * len(obs.events.time) + * u.deg, + ) events_camera_frame = events_altaz.transform_to(camera_frame) + events_camera_frame = SkyCoord( + lon=-events_camera_frame.lon, + lat=events_camera_frame.lat, + frame=camera_frame, + ) + + # Rotation to one observation + if rotate_to_obs is not None: + altaz_frame_obs = AltAz( + obstime=rotate_to_obs.events.time, location=rotate_to_obs.meta.location + ) + events_altaz_obs = rotate_to_obs.events.radec.transform_to(altaz_frame_obs) + pointing_altaz_obs = rotate_to_obs.get_pointing_icrs( + rotate_to_obs.events.time + ).transform_to(altaz_frame_obs) + + # Rotation to transform to camera frame + camera_frame_obs = SkyOffsetFrame( + origin=AltAz( + alt=pointing_altaz_obs.alt, + az=pointing_altaz_obs.az, + obstime=rotate_to_obs.events.time, + location=rotate_to_obs.meta.location, + ), + rotation=[ + 0.0, + ] + * len(rotate_to_obs.events.time) + * u.deg, + ) + events_camera_frame_obs = events_altaz_obs.transform_to(camera_frame_obs) + events_camera_frame_obs = SkyCoord( + lon=-events_camera_frame_obs.lon, + lat=events_camera_frame_obs.lat, + frame=camera_frame_obs, + ) + + pca_obs = PCA(n_components=2).fit( + np.stack([events_camera_frame_obs.lon, events_camera_frame_obs.lat]).T + ) + component0_x, component0_y = pca_obs.components_[0] + derot_angle_obs = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg + # mid Az of rotate_to_obs + az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az + + derot_angle_obs = gamma_pca_corr_magic(derot_angle_obs, az_obs) + + ##### + frame_centers = pointing_altaz.transform_to(camera_frame) + + sep = frame_centers.separation(events_camera_frame) + pos_angle = frame_centers.position_angle(events_camera_frame) + + pca = PCA(n_components=2).fit( + np.stack([events_camera_frame.lon, events_camera_frame.lat]).T + ) + component0_x, component0_y = pca.components_[0] + derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg + + az = obs.get_pointing_altaz(obs.tmid).az + derot_angle = gamma_pca_corr_magic(derot_angle, az) + + theory_diff = az_obs - az + rot_angle = derot_angle_obs - derot_angle + if derotate: + rot_angle = -derot_angle + print(f"Rot_angle{rot_angle}") + print(f"Theory angle {theory_diff.deg}") + # if np.abs(theory_diff) > 90 * u.deg: + # rot_angle += 180 * u.deg + + # rot_angle = derot_angle_obs - derot_angle + # print(derot_angle_obs, derot_angle) + # rot_angle = angle_rot - derot_angle + # while rot_angle > 90 * u.deg: + # rot_angle = rot_angle - 180 * u.deg + # print(rot_angle) + + # pos angle North to East + events_camera_frame = frame_centers.directional_offset_by( + position_angle=pos_angle - rot_angle, separation=sep + ) # Formatting data for the output camera_frame_events = obs.events.copy() - camera_frame_events.table['RA'] = events_camera_frame.lon - camera_frame_events.table['DEC'] = events_camera_frame.lat + camera_frame_events.table["RA"] = events_camera_frame.lon + camera_frame_events.table["DEC"] = events_camera_frame.lat camera_frame_obs_info = copy.deepcopy(obs.obs_info) - camera_frame_obs_info['RA_PNT'] = 0. - camera_frame_obs_info['DEC_PNT'] = 0. - obs_camera_frame = Observation(obs_id=obs.obs_id, - obs_info=camera_frame_obs_info, - events=camera_frame_events, - gti=obs.gti, - aeff=obs.aeff) - obs_camera_frame._location = obs.observatory_earth_location - - return obs_camera_frame - - def _transform_exclusion_region_to_camera_frame(self, pointing_altaz: AltAz) -> List[SkyRegion]: + camera_frame_obs_info["RA_PNT"] = 0.0 + camera_frame_obs_info["DEC_PNT"] = 0.0 + obs_camera_frame = Observation( + obs_id=obs.obs_id, + obs_info=camera_frame_obs_info, + events=camera_frame_events, + gti=obs.gti, + aeff=obs.aeff, + ) + obs_camera_frame._location = obs.meta.location + if rotate_to_obs is not None: + return obs_camera_frame, rot_angle # , component0_x, component0_y + + return obs_camera_frame, None + + def _transform_exclusion_region_to_camera_frame( + self, + pointing_altaz: AltAz, + rotate_to_obs: Optional[Observation] = None, + pca_angle: Optional[Angle] = None, + ) -> List[SkyRegion]: """ Transform the list of exclusion regions in sky frame into a list in camera frame. @@ -191,6 +330,12 @@ def _transform_exclusion_region_to_camera_frame(self, pointing_altaz: AltAz) -> ---------- pointing_altaz : astropy.coordinates.AltAz The pointing position of the telescope. + rotate_to_obs : gammapy.data.observations.Observation + All events of `obs` are rotate in camera coordinates to match + the azimuth angle of `rotate_to_obs`. + For stereoscopic telescopes, like MAGIC, the morphology changes with azimuth. + Therefore the statistics can be enhacned with this instead of creating a background + model for multiple azimuth bins. Returns ------- @@ -203,42 +348,83 @@ def _transform_exclusion_region_to_camera_frame(self, pointing_altaz: AltAz) -> If the region type is not supported. """ - camera_frame = SkyOffsetFrame(origin=pointing_altaz, - rotation=[0., ] * u.deg) + camera_frame = SkyOffsetFrame( + origin=pointing_altaz, + rotation=[ + 0.0, + ] + * u.deg, + ) exclude_region_camera_frame = [] for region in self.exclude_regions: if isinstance(region, CircleSkyRegion): center_coordinate = region.center center_coordinate_altaz = center_coordinate.transform_to(pointing_altaz) - center_coordinate_camera_frame = center_coordinate_altaz.transform_to(camera_frame) - center_coordinate_camera_frame_arb = SkyCoord(ra=center_coordinate_camera_frame.lon[0], - dec=center_coordinate_camera_frame.lat[0]) - exclude_region_camera_frame.append(CircleSkyRegion(center=center_coordinate_camera_frame_arb, - radius=region.radius)) + center_coordinate_camera_frame = center_coordinate_altaz.transform_to( + camera_frame + ) + if rotate_to_obs is not None: + frame_center = pointing_altaz.transform_to(camera_frame) + # rot_angle = ( + # rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az + # - pca_angle + # ) + sep = frame_center.separation(center_coordinate_camera_frame) + pos_angle = frame_center.position_angle( + center_coordinate_camera_frame + ) + center_coordinate_camera_frame = frame_center.directional_offset_by( + position_angle=pos_angle + pca_angle, separation=sep + ) + center_coordinate_camera_frame_arb = SkyCoord( + ra=-center_coordinate_camera_frame.lon[0], + dec=center_coordinate_camera_frame.lat[0], + ) + exclude_region_camera_frame.append( + CircleSkyRegion( + center=center_coordinate_camera_frame_arb, radius=region.radius + ) + ) elif isinstance(region, EllipseSkyRegion): center_coordinate = region.center center_coordinate_altaz = center_coordinate.transform_to(pointing_altaz) - center_coordinate_camera_frame = center_coordinate_altaz.transform_to(camera_frame) - width_coordinate = center_coordinate.directional_offset_by(region.angle,region.width) + center_coordinate_camera_frame = center_coordinate_altaz.transform_to( + camera_frame + ) + width_coordinate = center_coordinate.directional_offset_by( + region.angle, region.width + ) width_coordinate_altaz = width_coordinate.transform_to(pointing_altaz) - width_coordinate_camera_frame = width_coordinate_altaz.transform_to(camera_frame) - angle_camera_frame = center_coordinate_camera_frame.position_angle(width_coordinate_camera_frame).to(u.deg)[0] - center_coordinate_camera_frame_arb = SkyCoord(ra=center_coordinate_camera_frame.lon[0], - dec=center_coordinate_camera_frame.lat[0]) - exclude_region_camera_frame.append(EllipseSkyRegion(center=center_coordinate_camera_frame_arb, - width=region.width, height=region.height, - angle=angle_camera_frame)) + width_coordinate_camera_frame = width_coordinate_altaz.transform_to( + camera_frame + ) + angle_camera_frame = center_coordinate_camera_frame.position_angle( + width_coordinate_camera_frame + ).to(u.deg)[0] + center_coordinate_camera_frame_arb = SkyCoord( + ra=center_coordinate_camera_frame.lon[0], + dec=center_coordinate_camera_frame.lat[0], + ) + exclude_region_camera_frame.append( + EllipseSkyRegion( + center=center_coordinate_camera_frame_arb, + width=region.width, + height=region.height, + angle=angle_camera_frame, + ) + ) else: - raise Exception(f'{type(region)} region type not supported') + raise Exception(f"{type(region)} region type not supported") return exclude_region_camera_frame - def _create_map(self, - obs: Observation, - geom: WcsGeom, - exclude_regions: List[SkyRegion], - add_bkg: bool = False - ) -> Tuple[MapDataset, WcsNDMap]: + def _create_map( + self, + obs: Observation, + geom: WcsGeom, + exclude_regions: List[SkyRegion], + add_bkg: bool = False, + ) -> Tuple[MapDataset, WcsNDMap]: """ Create a map and the associated exclusion mask based on the given geometry and exclusion region. @@ -265,21 +451,25 @@ def _create_map(self, if add_bkg: maker = MapDatasetMaker(selection=["counts", "background"]) - maker_safe_mask = SafeMaskMaker(methods=["offset-max"], offset_max=self.max_offset) + maker_safe_mask = SafeMaskMaker( + methods=["offset-max"], offset_max=self.max_offset + ) geom_image = geom.to_image() - exclusion_mask = ~geom_image.region_mask(exclude_regions) if len(exclude_regions) > 0 else ~Map.from_geom( - geom_image) + exclusion_mask = ( + ~geom_image.region_mask(exclude_regions) + if len(exclude_regions) > 0 + else ~Map.from_geom(geom_image) + ) map_obs = maker.run(MapDataset.create(geom=geom), obs) map_obs = maker_safe_mask.run(map_obs, obs) return map_obs, exclusion_mask - def _create_sky_map(self, - obs: Observation, - add_bkg: bool = False - ) -> Tuple[MapDataset, WcsNDMap]: + def _create_sky_map( + self, obs: Observation, add_bkg: bool = False + ) -> Tuple[MapDataset, WcsNDMap]: """ Create the sky map and the associated exclusion mask based on the observation and the exclusion regions. @@ -298,12 +488,16 @@ def _create_sky_map(self, The exclusion mask. """ - geom_obs = WcsGeom.create(skydir=obs.get_pointing_icrs(obs.tmid), - npix=(self.n_bins_map, self.n_bins_map), - binsz=self.spatial_bin_size, - frame="icrs", - axes=[self.energy_axis]) - map_obs, exclusion_mask = self._create_map(obs, geom_obs, self.exclude_regions, add_bkg=add_bkg) + geom_obs = WcsGeom.create( + skydir=obs.get_pointing_icrs(obs.tmid), + npix=(self.n_bins_map, self.n_bins_map), + binsz=self.spatial_bin_size, + frame="icrs", + axes=[self.energy_axis], + ) + map_obs, exclusion_mask = self._create_map( + obs, geom_obs, self.exclude_regions, add_bkg=add_bkg + ) return map_obs, exclusion_mask @@ -323,18 +517,44 @@ def _compute_time_intervals_based_on_fov_rotation(self, obs: Observation) -> Tim """ # Determine time interval for cutting the obs as function of the rotation of the Fov - n_bin = max(2, int(np.rint( - ((obs.tstop - obs.tstart) / self.time_resolution).to_value(u.dimensionless_unscaled)))) + n_bin = max( + 2, + int( + np.rint( + ((obs.tstop - obs.tstart) / self.time_resolution).to_value( + u.dimensionless_unscaled + ) + ) + ), + ) time_axis = np.linspace(obs.tstart, obs.tstop, num=n_bin) - rotation_speed_fov = compute_rotation_speed_fov(time_axis, obs.get_pointing_icrs(obs.tmid), - obs.observatory_earth_location) - rotation_fov = cumulative_trapezoid(x=time_axis.unix_tai, - y=rotation_speed_fov.to_value(u.rad / u.s), - initial=0.) * u.rad + rotation_speed_fov = compute_rotation_speed_fov( + time_axis, obs.get_pointing_icrs(obs.tmid), obs.meta.location + ) + rotation_fov = ( + cumulative_trapezoid( + x=time_axis.unix_tai, + y=rotation_speed_fov.to_value(u.rad / u.s), + initial=0.0, + ) + * u.rad + ) distance_rotation_fov = rotation_fov.to_value(u.rad) * np.pi * self.max_offset - node_obs = distance_rotation_fov // (self.spatial_bin_size * self.max_fraction_pixel_rotation_fov) + node_obs = distance_rotation_fov // ( + self.spatial_bin_size * self.max_fraction_pixel_rotation_fov + ) change_node = node_obs[2:] != node_obs[1:-1] - time_interval = Time([obs.tstart, ] + [time_axis[1:-1][change_node], ] + [obs.tstop, ]) + time_interval = Time( + [ + obs.tstart, + ] + + [ + time_axis[1:-1][change_node], + ] + + [ + obs.tstop, + ] + ) return time_interval @@ -381,9 +601,9 @@ def _create_base_computation_map(self, observations: Observation) -> Tuple: """ pass - def _normalised_model_per_run(self, - observations: Observations, - acceptance_map: dict[int, BackgroundIRF]) -> dict[int, BackgroundIRF]: + def _normalised_model_per_run( + self, observations: Observations, acceptance_map: dict[int, BackgroundIRF] + ) -> dict[int, BackgroundIRF]: """ Normalised the acceptance model associated to each run to the events associated with the run @@ -410,33 +630,49 @@ def _normalised_model_per_run(self, modified_observation.bkg = acceptance_map[id_observation] # Fit the background model - logger.info('Fit to model to run ' + str(id_observation)) - map_obs, exclusion_mask = self._create_sky_map(modified_observation, add_bkg=True) - maker_FoV_background = FoVBackgroundMaker(method='fit', exclusion_mask=exclusion_mask) + logger.info("Fit to model to run " + str(id_observation)) + map_obs, exclusion_mask = self._create_sky_map( + modified_observation, add_bkg=True + ) + maker_FoV_background = FoVBackgroundMaker( + method="fit", exclusion_mask=exclusion_mask + ) map_obs = maker_FoV_background.run(map_obs) # Extract the normalisation parameters = map_obs.models.to_parameters_table() - norm_background = parameters[parameters['name'] == 'norm']['value'][0] + norm_background = parameters[parameters["name"] == "norm"]["value"][0] - if norm_background < 0.: + if norm_background < 0.0: logger.error( - 'Invalid normalisation value for run ' + str(id_observation) + ' : ' + str( - norm_background) + ', normalisation set back to 1') - norm_background = 1. + "Invalid normalisation value for run " + + str(id_observation) + + " : " + + str(norm_background) + + ", normalisation set back to 1" + ) + norm_background = 1.0 elif norm_background > 1.5 or norm_background < 0.5: logger.warning( - 'High correction of the background normalisation for run ' + str(id_observation) + ' : ' + str( - norm_background)) + "High correction of the background normalisation for run " + + str(id_observation) + + " : " + + str(norm_background) + ) # Apply normalisation to the background model - normalised_acceptance_map[id_observation] = copy.deepcopy(acceptance_map[id_observation]) - normalised_acceptance_map[id_observation].data = normalised_acceptance_map[ - id_observation].data * norm_background + normalised_acceptance_map[id_observation] = copy.deepcopy( + acceptance_map[id_observation] + ) + normalised_acceptance_map[id_observation].data = ( + normalised_acceptance_map[id_observation].data * norm_background + ) return normalised_acceptance_map - def _compute_time_intervals_based_on_zenith_bin(self, obs: Observation, edge_zenith_bin: u.Quantity) -> Time: + def _compute_time_intervals_based_on_zenith_bin( + self, obs: Observation, edge_zenith_bin: u.Quantity + ) -> Time: """ Calculate time intervals based on an input zenith binning @@ -454,27 +690,49 @@ def _compute_time_intervals_based_on_zenith_bin(self, obs: Observation, edge_zen """ # Create the time axis - n_bin = max(2, int(np.rint( - ((obs.tstop - obs.tstart) / self.time_resolution).to_value(u.dimensionless_unscaled)))) + n_bin = max( + 2, + int( + np.rint( + ((obs.tstop - obs.tstart) / self.time_resolution).to_value( + u.dimensionless_unscaled + ) + ) + ), + ) time_axis = np.linspace(obs.tstart, obs.tstop, num=n_bin) # Compute the zenith for each evaluation time with erfa_astrom.set(ErfaAstromInterpolator(1000 * u.s)): altaz_coordinates = obs.get_pointing_altaz(time_axis) zenith_values = altaz_coordinates.zen - if np.any(zenith_values < np.min(edge_zenith_bin)) or np.any(zenith_values > np.max(edge_zenith_bin)): - logger.error('Run with zenith value outside of the considered range for zenith binning') + if np.any(zenith_values < np.min(edge_zenith_bin)) or np.any( + zenith_values > np.max(edge_zenith_bin) + ): + logger.error( + "Run with zenith value outside of the considered range for zenith binning" + ) # Split the time interval to transition between zenith bin id_bin = np.digitize(zenith_values, edge_zenith_bin) bin_transition = id_bin[2:] != id_bin[1:-1] - time_interval = Time([obs.tstart, ] + [time_axis[1:-1][bin_transition], ] + [obs.tstop, ]) + time_interval = Time( + [ + obs.tstart, + ] + + [ + time_axis[1:-1][bin_transition], + ] + + [ + obs.tstop, + ] + ) return time_interval - def create_model_cos_zenith_binned(self, - observations: Observations - ) -> BackgroundCollectionZenith: + def create_model_cos_zenith_binned( + self, observations: Observations + ) -> BackgroundCollectionZenith: """ Calculate a model for each cos zenith bin @@ -491,35 +749,63 @@ def create_model_cos_zenith_binned(self, """ # Determine binning method. Convention : per_wobble methods have negative values - methods = {'min_livetime': 1, 'min_livetime_per_wobble': -1, 'min_n_observation': 2, - 'min_n_observation_per_wobble': -2} + methods = { + "min_livetime": 1, + "min_livetime_per_wobble": -1, + "min_n_observation": 2, + "min_n_observation_per_wobble": -2, + } try: i_method = methods[self.cos_zenith_binning_method] except KeyError: - logger.error(f" KeyError : {self.cos_zenith_binning_method} not a valid zenith binning method.\nValid " - f"methods are {[*methods]}") + logger.error( + f" KeyError : {self.cos_zenith_binning_method} not a valid zenith binning method.\nValid " + f"methods are {[*methods]}" + ) raise per_wobble = i_method < 0 # Initial bins edges - cos_zenith_bin = np.sort(np.arange(1.0, 0. - self.initial_cos_zenith_binning, -self.initial_cos_zenith_binning)) + cos_zenith_bin = np.sort( + np.arange( + 1.0, + 0.0 - self.initial_cos_zenith_binning, + -self.initial_cos_zenith_binning, + ) + ) zenith_bin = np.rad2deg(np.arccos(cos_zenith_bin)) * u.deg # Cut observations if requested if self.zenith_binning_run_splitting: if abs(i_method) == 2: - logger.warning('Using a zenith binning requirement based on n_observation while using run splitting is not recommended and could lead to poor models. We recommend switching to a binning requirement based on livetime.') + logger.warning( + "Using a zenith binning requirement based on n_observation while using run splitting is not recommended and could lead to poor models. We recommend switching to a binning requirement based on livetime." + ) compute_observations = Observations() for obs in observations: - time_interval = self._compute_time_intervals_based_on_zenith_bin(obs, zenith_bin) + time_interval = self._compute_time_intervals_based_on_zenith_bin( + obs, zenith_bin + ) for i in range(len(time_interval) - 1): - compute_observations.append(obs.select_time(Time([time_interval[i], time_interval[i + 1]]))) + compute_observations.append( + obs.select_time(Time([time_interval[i], time_interval[i + 1]])) + ) else: compute_observations = observations # Determine initial bins values - cos_zenith_observations = np.array([np.cos(obs.get_pointing_altaz(obs.tmid).zen) for obs in compute_observations]) - livetime_observations = np.array([obs.observation_live_time_duration.to_value(u.s) for obs in compute_observations]) + cos_zenith_observations = np.array( + [ + np.cos(obs.get_pointing_altaz(obs.tmid).zen) + for obs in compute_observations + ] + ) + livetime_observations = np.array( + [ + obs.observation_live_time_duration.to_value(u.s) + for obs in compute_observations + ] + ) # Select the quantity used to count observations if i_method in [-1, 1]: @@ -530,19 +816,27 @@ def create_model_cos_zenith_binned(self, # Gather runs per separation angle or all together. Define the minimum multiplicity (-1) to create a zenith bin. if per_wobble: wobble_observations = np.array( - get_unique_wobble_pointings(compute_observations, self.max_angular_separation_wobble)) + get_unique_wobble_pointings( + compute_observations, self.max_angular_separation_wobble + ) + ) multiplicity_wob = 1 else: - wobble_observations = np.full(len(cos_zenith_observations), 'any', dtype=np.object_) + wobble_observations = np.full( + len(cos_zenith_observations), "any", dtype=np.object_ + ) multiplicity_wob = 0 cumsum_variable = {} for wobble in np.unique(wobble_observations): # Create an array of cumulative weight of the selected variable vs cos(zenith) - cumsum_variable[wobble] = np.cumsum(np.histogram(cos_zenith_observations[wobble_observations == wobble], - bins=cos_zenith_bin, - weights=cut_variable_weights[ - wobble_observations == wobble])[0]) + cumsum_variable[wobble] = np.cumsum( + np.histogram( + cos_zenith_observations[wobble_observations == wobble], + bins=cos_zenith_bin, + weights=cut_variable_weights[wobble_observations == wobble], + )[0] + ) # Initiate the list of index of selected zenith bin edges zenith_selected = [0] @@ -552,8 +846,12 @@ def create_model_cos_zenith_binned(self, while i < n: # For each wobble, find the index of the first zenith which fulfills the zd binning criteria if any # Then concatenate and sort the results for all wobbles - candidate_i_per_wobble = [np.nonzero(cum_cut_variable >= self.cos_zenith_binning_parameter_value)[0][:1] - for cum_cut_variable in cumsum_variable.values()] # nonzero is assumed sorted + candidate_i_per_wobble = [ + np.nonzero(cum_cut_variable >= self.cos_zenith_binning_parameter_value)[ + 0 + ][:1] + for cum_cut_variable in cumsum_variable.values() + ] # nonzero is assumed sorted candidate_i = np.sort(np.concatenate(candidate_i_per_wobble)) if len(candidate_i) > multiplicity_wob: @@ -579,11 +877,17 @@ def create_model_cos_zenith_binned(self, for i in range((len(cos_zenith_bin) - 1)): binned_observations.append(Observations()) for obs in compute_observations: - binned_observations[np.digitize(np.cos(obs.get_pointing_altaz(obs.tmid).zen), cos_zenith_bin) - 1].append( - obs) + binned_observations[ + np.digitize( + np.cos(obs.get_pointing_altaz(obs.tmid).zen), cos_zenith_bin + ) + - 1 + ].append(obs) # Compute the model for each bin - binned_model = [self.create_acceptance_map(binned_obs) for binned_obs in binned_observations] + binned_model = [ + self.create_acceptance_map(binned_obs) for binned_obs in binned_observations + ] # Determine the center of the bin (weighted as function of the livetime of each observation) bin_center = [] @@ -592,39 +896,56 @@ def create_model_cos_zenith_binned(self, livetime_per_obs = [] for obs in binned_observations[i]: weighted_cos_zenith_bin_per_obs.append( - obs.observation_live_time_duration * np.cos(obs.get_pointing_altaz(obs.tmid).zen)) + obs.observation_live_time_duration + * np.cos(obs.get_pointing_altaz(obs.tmid).zen) + ) livetime_per_obs.append(obs.observation_live_time_duration) - bin_center.append(np.sum([wcos.value for wcos in weighted_cos_zenith_bin_per_obs]) / np.sum( - [livet.value for livet in livetime_per_obs])) + bin_center.append( + np.sum([wcos.value for wcos in weighted_cos_zenith_bin_per_obs]) + / np.sum([livet.value for livet in livetime_per_obs]) + ) logger.info(f"cos zenith bin edges: {list(np.round(cos_zenith_bin, 2))}") logger.info(f"cos zenith bin centers: {list(np.round(bin_center, 2))}") - logger.info(f"observation per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin)[0])}") - logger.info(f"livetime per bin [s]: " + - f"{list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=livetime_observations)[0].astype(int))}") + logger.info( + f"observation per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin)[0])}" + ) + logger.info( + f"livetime per bin [s]: " + + f"{list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=livetime_observations)[0].astype(int))}" + ) if per_wobble: - wobble_observations_bool_arr = [(np.array(wobble_observations.tolist()) == wobble) for wobble in - np.unique(np.array(wobble_observations))] - livetime_observations_and_wobble = [np.array(livetime_observations) * wobble_bool for wobble_bool in - wobble_observations_bool_arr] + wobble_observations_bool_arr = [ + (np.array(wobble_observations.tolist()) == wobble) + for wobble in np.unique(np.array(wobble_observations)) + ] + livetime_observations_and_wobble = [ + np.array(livetime_observations) * wobble_bool + for wobble_bool in wobble_observations_bool_arr + ] for i, wobble in enumerate(np.unique(np.array(wobble_observations))): logger.info( - f"{wobble} observation per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=1 * wobble_observations_bool_arr[i])[0])}") + f"{wobble} observation per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=1 * wobble_observations_bool_arr[i])[0])}" + ) logger.info( - f"{wobble} livetime per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=livetime_observations_and_wobble[i])[0].astype(int))}") + f"{wobble} livetime per bin: {list(np.histogram(cos_zenith_observations, bins=cos_zenith_bin, weights=livetime_observations_and_wobble[i])[0].astype(int))}" + ) # Create the dict for output of the function collection_binned_model = BackgroundCollectionZenith() for i in range(len(binned_model)): - collection_binned_model[np.rad2deg(np.arccos(bin_center[i]))] = binned_model[i] + collection_binned_model[np.rad2deg(np.arccos(bin_center[i]))] = ( + binned_model[i] + ) return collection_binned_model - def create_acceptance_map_cos_zenith_binned(self, - observations: Observations, - off_observations: Observations = None, - base_model: BackgroundCollectionZenith = None - ) -> dict[int, BackgroundIRF]: + def create_acceptance_map_cos_zenith_binned( + self, + observations: Observations, + off_observations: Observations = None, + base_model: BackgroundCollectionZenith = None, + ) -> dict[int, BackgroundIRF]: """ Calculate an acceptance map per run using cos zenith binning @@ -648,14 +969,22 @@ def create_acceptance_map_cos_zenith_binned(self, if off_observations is None: off_observations = observations elif base_model is not None: - logger.warning('The off observations provided will be ignored as a base model has been provided.') + logger.warning( + "The off observations provided will be ignored as a base model has been provided." + ) # If needed produce the zenith binned model - if base_model is not None and not isinstance(base_model, BackgroundCollectionZenith): - error_message = 'The models should be provided as a BackgroundCollectionZenith object' + if base_model is not None and not isinstance( + base_model, BackgroundCollectionZenith + ): + error_message = ( + "The models should be provided as a BackgroundCollectionZenith object" + ) logger.error(error_message) raise BackgroundModelFormatException(error_message) - dict_binned_model = base_model or self.create_model_cos_zenith_binned(off_observations) + dict_binned_model = base_model or self.create_model_cos_zenith_binned( + off_observations + ) cos_zenith_model = [] key_model = [] @@ -671,53 +1000,72 @@ def create_acceptance_map_cos_zenith_binned(self, # Find the closest model for each observation and associate it to each observation acceptance_map = {} if len(cos_zenith_model) <= 1: - logger.warning('Only one zenith bin, zenith binning deactivated') + logger.warning("Only one zenith bin, zenith binning deactivated") for obs in observations: if self.use_mini_irf_computation: - evaluation_time, observation_time = get_time_mini_irf(obs, self.mini_irf_time_resolution) - - data_obs_all = np.zeros( - tuple([len(evaluation_time), ] + list(dict_binned_model[key_model[0]].data.shape))) * \ - dict_binned_model[key_model[0]].unit + evaluation_time, observation_time = get_time_mini_irf( + obs, self.mini_irf_time_resolution + ) + + data_obs_all = ( + np.zeros( + tuple( + [ + len(evaluation_time), + ] + + list(dict_binned_model[key_model[0]].data.shape) + ) + ) + * dict_binned_model[key_model[0]].unit + ) for i in range(len(evaluation_time)): - cos_zenith_observation = np.cos(obs.get_pointing_altaz(evaluation_time[i]).zen) - key_closest_model = key_model[(np.abs(cos_zenith_model - cos_zenith_observation)).argmin()] + cos_zenith_observation = np.cos( + obs.get_pointing_altaz(evaluation_time[i]).zen + ) + key_closest_model = key_model[ + (np.abs(cos_zenith_model - cos_zenith_observation)).argmin() + ] selected_model_bin = dict_binned_model[key_closest_model] data_obs_all[i] = selected_model_bin.data * selected_model_bin.unit data_obs = generate_irf_from_mini_irf(data_obs_all, observation_time) if type_model is Background2D: - acceptance_map[obs.obs_id] = Background2D(axes=axes_model, - data=data_obs) + acceptance_map[obs.obs_id] = Background2D( + axes=axes_model, data=data_obs + ) elif type_model is Background3D: - acceptance_map[obs.obs_id] = Background3D(axes=axes_model, - data=data_obs, - fov_alignment=FoVAlignment.ALTAZ) + acceptance_map[obs.obs_id] = Background3D( + axes=axes_model, data=data_obs, fov_alignment=FoVAlignment.ALTAZ + ) else: - raise Exception('Unknown background format') + raise Exception("Unknown background format") else: cos_zenith_observation = np.cos(obs.get_pointing_altaz(obs.tmid).zen) - key_closest_model = key_model[(np.abs(cos_zenith_model - cos_zenith_observation)).argmin()] + key_closest_model = key_model[ + (np.abs(cos_zenith_model - cos_zenith_observation)).argmin() + ] acceptance_map[obs.obs_id] = dict_binned_model[key_closest_model] return acceptance_map - def _create_interpolation_function(self, base_model: BackgroundCollectionZenith) -> interp1d: + def _create_interpolation_function( + self, base_model: BackgroundCollectionZenith + ) -> interp1d: """ - Create the function that will perform the interpolation + Create the function that will perform the interpolation - Parameters - ---------- - base_model : dict of gammapy.irf.background.BackgroundIRF - The binned base model - Each key of the dictionary should correspond to the zenith in degree of the model + Parameters + ---------- + base_model : dict of gammapy.irf.background.BackgroundIRF + The binned base model + Each key of the dictionary should correspond to the zenith in degree of the model - Returns - ------- - interp_func : scipy.interpolate.interp1d - The object that could be call directly for performing the interpolation + Returns + ------- + interp_func : scipy.interpolate.interp1d + The object that could be call directly for performing the interpolation """ # Reshape the base model @@ -728,20 +1076,36 @@ def _create_interpolation_function(self, base_model: BackgroundCollectionZenith) cos_zenith_model.append(np.cos(np.deg2rad(k))) cos_zenith_model = np.array(cos_zenith_model) - data_cube = np.zeros(tuple([len(binned_model), ] + list(binned_model[0].data.shape))) * binned_model[0].unit + data_cube = ( + np.zeros( + tuple( + [ + len(binned_model), + ] + + list(binned_model[0].data.shape) + ) + ) + * binned_model[0].unit + ) for i in range(len(binned_model)): data_cube[i] = binned_model[i].data * binned_model[i].unit - if self.interpolation_type == 'log': - interp_func = interp1d(x=cos_zenith_model, - y=np.log10(data_cube.to_value( - binned_model[0].unit) + self.threshold_value_log_interpolation), - axis=0, - fill_value='extrapolate') - elif self.interpolation_type == 'linear': - interp_func = interp1d(x=cos_zenith_model, - y=data_cube.to_value(binned_model[0].unit), - axis=0, - fill_value='extrapolate') + if self.interpolation_type == "log": + interp_func = interp1d( + x=cos_zenith_model, + y=np.log10( + data_cube.to_value(binned_model[0].unit) + + self.threshold_value_log_interpolation + ), + axis=0, + fill_value="extrapolate", + ) + elif self.interpolation_type == "linear": + interp_func = interp1d( + x=cos_zenith_model, + y=data_cube.to_value(binned_model[0].unit), + axis=0, + fill_value="extrapolate", + ) else: raise Exception("Unknown interpolation type") @@ -749,62 +1113,77 @@ def _create_interpolation_function(self, base_model: BackgroundCollectionZenith) def _background_cleaning(self, background_model): """ - Is cleaning the background model from suspicious values not compatible with neighbour pixels. + Is cleaning the background model from suspicious values not compatible with neighbour pixels. - Parameters - ---------- - background_model : numpy.array - The background model to be cleaned + Parameters + ---------- + background_model : numpy.array + The background model to be cleaned - Returns - ------- - background_model : numpy.array - The background model cleaned + Returns + ------- + background_model : numpy.array + The background model cleaned """ base_model = background_model.copy() final_model = background_model.copy() i = 0 - while (i < 1 or not np.allclose(base_model, final_model)) and (i < self.max_cleaning_iteration): + while (i < 1 or not np.allclose(base_model, final_model)) and ( + i < self.max_cleaning_iteration + ): base_model = final_model.copy() i += 1 - count_valid_neighbour_condition_energy = compute_neighbour_condition_validation(base_model, axis=0, - relative_threshold=self.interpolation_cleaning_energy_relative_threshold) - count_valid_neighbour_condition_spatial = compute_neighbour_condition_validation(base_model, axis=1, - relative_threshold=self.interpolation_cleaning_spatial_relative_threshold) + count_valid_neighbour_condition_energy = compute_neighbour_condition_validation( + base_model, + axis=0, + relative_threshold=self.interpolation_cleaning_energy_relative_threshold, + ) + count_valid_neighbour_condition_spatial = compute_neighbour_condition_validation( + base_model, + axis=1, + relative_threshold=self.interpolation_cleaning_spatial_relative_threshold, + ) if base_model.ndim == 3: - count_valid_neighbour_condition_spatial += compute_neighbour_condition_validation(base_model, axis=2, - relative_threshold=self.interpolation_cleaning_spatial_relative_threshold) + count_valid_neighbour_condition_spatial += compute_neighbour_condition_validation( + base_model, + axis=2, + relative_threshold=self.interpolation_cleaning_spatial_relative_threshold, + ) mask_energy = count_valid_neighbour_condition_energy > 0 - mask_spatial = count_valid_neighbour_condition_spatial > (1 if base_model.ndim == 3 else 0) + mask_spatial = count_valid_neighbour_condition_spatial > ( + 1 if base_model.ndim == 3 else 0 + ) mask_valid = np.logical_and(mask_energy, mask_spatial) - final_model[~mask_valid] = 0. + final_model[~mask_valid] = 0.0 return final_model - def _get_interpolated_background(self, interp_func: interp1d, zenith: u.Quantity) -> np.array: + def _get_interpolated_background( + self, interp_func: interp1d, zenith: u.Quantity + ) -> np.array: """ - Create the function that will perform the interpolation + Create the function that will perform the interpolation - Parameters - ---------- - interp_func : scipy.interpolate.interp1d - The interpolation function - zenith : u.Quantity - The zenith for which the interpolated background should be created + Parameters + ---------- + interp_func : scipy.interpolate.interp1d + The interpolation function + zenith : u.Quantity + The zenith for which the interpolated background should be created - Returns - ------- - interp_bkg : numpy.array - The object that could be call directly for performing the interpolation + Returns + ------- + interp_bkg : numpy.array + The object that could be call directly for performing the interpolation """ - if self.interpolation_type == 'log': - interp_bkg = (10. ** interp_func(np.cos(zenith))) - interp_bkg[interp_bkg < 100 * self.threshold_value_log_interpolation] = 0. - elif self.interpolation_type == 'linear': + if self.interpolation_type == "log": + interp_bkg = 10.0 ** interp_func(np.cos(zenith)) + interp_bkg[interp_bkg < 100 * self.threshold_value_log_interpolation] = 0.0 + elif self.interpolation_type == "linear": interp_bkg = interp_func(np.cos(zenith)) else: raise Exception("Unknown interpolation type") @@ -814,11 +1193,12 @@ def _get_interpolated_background(self, interp_func: interp1d, zenith: u.Quantity return interp_bkg - def create_acceptance_map_cos_zenith_interpolated(self, - observations: Observations, - off_observations: Observations = None, - base_model: BackgroundCollectionZenith = None - ) -> dict[int, BackgroundIRF]: + def create_acceptance_map_cos_zenith_interpolated( + self, + observations: Observations, + off_observations: Observations = None, + base_model: BackgroundCollectionZenith = None, + ) -> dict[int, BackgroundIRF]: """ Calculate an acceptance map per run using cos zenith binning and interpolation @@ -842,20 +1222,30 @@ def create_acceptance_map_cos_zenith_interpolated(self, if off_observations is None: off_observations = observations elif base_model is not None: - logger.warning('The off observations provided will be ignored as a base model has been provided.') + logger.warning( + "The off observations provided will be ignored as a base model has been provided." + ) # If needed produce the zenith binned model - if base_model is not None and not isinstance(base_model, BackgroundCollectionZenith): - error_message = 'The models should be provided as a BackgroundCollectionZenith object' + if base_model is not None and not isinstance( + base_model, BackgroundCollectionZenith + ): + error_message = ( + "The models should be provided as a BackgroundCollectionZenith object" + ) logger.error(error_message) raise BackgroundModelFormatException(error_message) - dict_binned_model = base_model or self.create_model_cos_zenith_binned(off_observations) + dict_binned_model = base_model or self.create_model_cos_zenith_binned( + off_observations + ) acceptance_map = {} if len(dict_binned_model) <= 1: - logger.warning('Only one zenith bin, zenith interpolation deactivated') + logger.warning("Only one zenith bin, zenith interpolation deactivated") for obs in observations: - acceptance_map[obs.obs_id] = dict_binned_model[dict_binned_model.zenith[0]] + acceptance_map[obs.obs_id] = dict_binned_model[ + dict_binned_model.zenith[0] + ] else: # Determine model properties type_model = type(dict_binned_model[dict_binned_model.zenith[0]]) @@ -867,39 +1257,57 @@ def create_acceptance_map_cos_zenith_interpolated(self, interp_func = self._create_interpolation_function(dict_binned_model) for obs in observations: if self.use_mini_irf_computation: - evaluation_time, observation_time = get_time_mini_irf(obs, self.mini_irf_time_resolution) - - data_obs_all = np.zeros(tuple([len(evaluation_time), ] + list(shape_model))) + evaluation_time, observation_time = get_time_mini_irf( + obs, self.mini_irf_time_resolution + ) + + data_obs_all = np.zeros( + tuple( + [ + len(evaluation_time), + ] + + list(shape_model) + ) + ) for i in range(len(evaluation_time)): - data_obs_bin = self._get_interpolated_background(interp_func, - obs.get_pointing_altaz(evaluation_time[i]).zen) + data_obs_bin = self._get_interpolated_background( + interp_func, obs.get_pointing_altaz(evaluation_time[i]).zen + ) data_obs_all[i, :, :] = data_obs_bin - data_obs = generate_irf_from_mini_irf(data_obs_all, observation_time) + data_obs = generate_irf_from_mini_irf( + data_obs_all, observation_time + ) else: - data_obs = self._get_interpolated_background(interp_func, obs.get_pointing_altaz(obs.tmid).zen) + data_obs = self._get_interpolated_background( + interp_func, obs.get_pointing_altaz(obs.tmid).zen + ) if type_model is Background2D: - acceptance_map[obs.obs_id] = Background2D(axes=axes_model, - data=data_obs * unit_model) + acceptance_map[obs.obs_id] = Background2D( + axes=axes_model, data=data_obs * unit_model + ) elif type_model is Background3D: - acceptance_map[obs.obs_id] = Background3D(axes=axes_model, - data=data_obs * unit_model, - fov_alignment=FoVAlignment.ALTAZ) + acceptance_map[obs.obs_id] = Background3D( + axes=axes_model, + data=data_obs * unit_model, + fov_alignment=FoVAlignment.ALTAZ, + ) else: - raise Exception('Unknown background format') + raise Exception("Unknown background format") return acceptance_map - def create_acceptance_map_per_observation(self, - observations: Observations, - zenith_binning: bool = False, - zenith_interpolation: bool = False, - runwise_normalisation: bool = True, - off_observations: Observations = None, - base_model: Union[BackgroundCollectionZenith, BackgroundIRF] = None, - ) -> dict[int, BackgroundIRF]: + def create_acceptance_map_per_observation( + self, + observations: Observations, + zenith_binning: bool = False, + zenith_interpolation: bool = False, + runwise_normalisation: bool = True, + off_observations: Observations = None, + base_model: Union[BackgroundCollectionZenith, BackgroundIRF] = None, + ) -> dict[int, BackgroundIRF]: """ Calculate an acceptance map with the norm adjusted for each run @@ -928,31 +1336,43 @@ def create_acceptance_map_per_observation(self, acceptance_map = {} if zenith_interpolation: - acceptance_map = self.create_acceptance_map_cos_zenith_interpolated(observations=observations, - off_observations=off_observations, - base_model=base_model) + acceptance_map = self.create_acceptance_map_cos_zenith_interpolated( + observations=observations, + off_observations=off_observations, + base_model=base_model, + ) elif zenith_binning: - acceptance_map = self.create_acceptance_map_cos_zenith_binned(observations=observations, - off_observations=off_observations, - base_model=base_model) + acceptance_map = self.create_acceptance_map_cos_zenith_binned( + observations=observations, + off_observations=off_observations, + base_model=base_model, + ) else: if off_observations is None: off_observations = observations elif base_model is not None: - logger.warning('The off observations provided will be ignored as a base model has been provided.') + logger.warning( + "The off observations provided will be ignored as a base model has been provided." + ) if base_model is not None: if not isinstance(base_model, BackgroundIRF): - error_message = 'The model provided should be a gammapy BackgroundIRF object' + error_message = ( + "The model provided should be a gammapy BackgroundIRF object" + ) logger.error(error_message) raise BackgroundModelFormatException(error_message) unique_base_acceptance_map = base_model else: - unique_base_acceptance_map = self.create_acceptance_map(off_observations) + unique_base_acceptance_map = self.create_acceptance_map( + off_observations + ) for obs in observations: acceptance_map[obs.obs_id] = unique_base_acceptance_map if runwise_normalisation: - acceptance_map = self._normalised_model_per_run(observations, acceptance_map) + acceptance_map = self._normalised_model_per_run( + observations, acceptance_map + ) return acceptance_map diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index ebffa74..136c569 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -5,12 +5,13 @@ from astropy.coordinates.erfa_astrom import erfa_astrom, ErfaAstromInterpolator import astropy.units as u import numpy as np -from gammapy.data import Observations +from gammapy.data import Observations, Observation from gammapy.datasets import MapDataset from gammapy.irf import FoVAlignment, Background3D from gammapy.maps import WcsNDMap, Map, MapAxis, RegionGeom from iminuit import Minuit from regions import SkyRegion +from scipy.ndimage import rotate from .base_acceptance_map_creator import BaseAcceptanceMapCreator from .modeling import FIT_FUNCTION, log_factorial, log_poisson @@ -19,29 +20,30 @@ class Grid3DAcceptanceMapCreator(BaseAcceptanceMapCreator): - - def __init__(self, - energy_axis: MapAxis, - offset_axis: MapAxis, - oversample_map: int = 10, - exclude_regions: Optional[List[SkyRegion]] = None, - cos_zenith_binning_method: str = 'min_livetime', - cos_zenith_binning_parameter_value: int = 3600, - initial_cos_zenith_binning: float = 0.01, - max_angular_separation_wobble: u.Quantity = 0.4 * u.deg, - zenith_binning_run_splitting: bool = False, - max_fraction_pixel_rotation_fov: float = 0.5, - time_resolution: u.Quantity = 0.1 * u.s, - use_mini_irf_computation: bool = False, - mini_irf_time_resolution: u.Quantity = 1. * u.min, - method='stack', - fit_fnc='gaussian2d', - fit_seeds=None, - fit_bounds=None, - interpolation_type: str = 'linear', - activate_interpolation_cleaning: bool = False, - interpolation_cleaning_energy_relative_threshold: float = 1e-4, - interpolation_cleaning_spatial_relative_threshold: float = 1e-2) -> None: + def __init__( + self, + energy_axis: MapAxis, + offset_axis: MapAxis, + oversample_map: int = 10, + exclude_regions: Optional[List[SkyRegion]] = None, + cos_zenith_binning_method: str = "min_livetime", + cos_zenith_binning_parameter_value: int = 3600, + initial_cos_zenith_binning: float = 0.01, + max_angular_separation_wobble: u.Quantity = 0.4 * u.deg, + zenith_binning_run_splitting: bool = False, + max_fraction_pixel_rotation_fov: float = 0.5, + time_resolution: u.Quantity = 0.1 * u.s, + use_mini_irf_computation: bool = False, + mini_irf_time_resolution: u.Quantity = 1.0 * u.min, + method="stack", + fit_fnc="gaussian2d", + fit_seeds=None, + fit_bounds=None, + interpolation_type: str = "linear", + activate_interpolation_cleaning: bool = False, + interpolation_cleaning_energy_relative_threshold: float = 1e-4, + interpolation_cleaning_spatial_relative_threshold: float = 1e-2, + ) -> None: """ Create the class for calculating 3D grid acceptance model @@ -102,12 +104,14 @@ def __init__(self, # Compute parameters for internal map self.offset_axis = offset_axis if not np.allclose(self.offset_axis.bin_width, self.offset_axis.bin_width[0]): - raise Exception('Support only regular linear bin for offset axis') - if not np.isclose(self.offset_axis.edges[0], 0. * u.deg): - raise Exception('Offset axis need to start at 0') + raise Exception("Support only regular linear bin for offset axis") + if not np.isclose(self.offset_axis.edges[0], 0.0 * u.deg): + raise Exception("Offset axis need to start at 0") self.oversample_map = oversample_map - spatial_resolution = np.min( - np.abs(self.offset_axis.edges[1:] - self.offset_axis.edges[:-1])) / self.oversample_map + spatial_resolution = ( + np.min(np.abs(self.offset_axis.edges[1:] - self.offset_axis.edges[:-1])) + / self.oversample_map + ) max_offset = np.max(self.offset_axis.edges) self.method = method @@ -116,23 +120,25 @@ def __init__(self, self.fit_bounds = fit_bounds # Initiate upper instance - super().__init__(energy_axis=energy_axis, - max_offset=max_offset, - spatial_resolution=spatial_resolution, - exclude_regions=exclude_regions, - cos_zenith_binning_method=cos_zenith_binning_method, - cos_zenith_binning_parameter_value=cos_zenith_binning_parameter_value, - initial_cos_zenith_binning=initial_cos_zenith_binning, - max_angular_separation_wobble=max_angular_separation_wobble, - zenith_binning_run_splitting=zenith_binning_run_splitting, - max_fraction_pixel_rotation_fov=max_fraction_pixel_rotation_fov, - time_resolution=time_resolution, - use_mini_irf_computation=use_mini_irf_computation, - mini_irf_time_resolution=mini_irf_time_resolution, - interpolation_type=interpolation_type, - activate_interpolation_cleaning=activate_interpolation_cleaning, - interpolation_cleaning_energy_relative_threshold=interpolation_cleaning_energy_relative_threshold, - interpolation_cleaning_spatial_relative_threshold=interpolation_cleaning_spatial_relative_threshold) + super().__init__( + energy_axis=energy_axis, + max_offset=max_offset, + spatial_resolution=spatial_resolution, + exclude_regions=exclude_regions, + cos_zenith_binning_method=cos_zenith_binning_method, + cos_zenith_binning_parameter_value=cos_zenith_binning_parameter_value, + initial_cos_zenith_binning=initial_cos_zenith_binning, + max_angular_separation_wobble=max_angular_separation_wobble, + zenith_binning_run_splitting=zenith_binning_run_splitting, + max_fraction_pixel_rotation_fov=max_fraction_pixel_rotation_fov, + time_resolution=time_resolution, + use_mini_irf_computation=use_mini_irf_computation, + mini_irf_time_resolution=mini_irf_time_resolution, + interpolation_type=interpolation_type, + activate_interpolation_cleaning=activate_interpolation_cleaning, + interpolation_cleaning_energy_relative_threshold=interpolation_cleaning_energy_relative_threshold, + interpolation_cleaning_spatial_relative_threshold=interpolation_cleaning_spatial_relative_threshold, + ) def fit_background(self, count_map, exp_map_total, exp_map): centers = self.offset_axis.center.to_value(u.deg) @@ -143,7 +149,9 @@ def fit_background(self, count_map, exp_map_total, exp_map): try: fnc = FIT_FUNCTION[self.fit_fnc] except KeyError: - logger.error(f"Invalid built-in fit_fnc. Use {FIT_FUNCTION.keys()} or a custom function.") + logger.error( + f"Invalid built-in fit_fnc. Use {FIT_FUNCTION.keys()} or a custom function." + ) raise raw_seeds = fnc.default_seeds.copy() bounds = fnc.default_bounds.copy() @@ -154,15 +162,19 @@ def fit_background(self, count_map, exp_map_total, exp_map): if self.fit_bounds is not None: bounds.update(self.fit_bounds) - mask = exp_map > 0 # Handles fully overlapping exclusion regions in the 'size' seed computation + mask = ( + exp_map > 0 + ) # Handles fully overlapping exclusion regions in the 'size' seed computation # Seeds the charge normalisation to the observed counts corrected for exclusion region reduction to exposure - raw_seeds['size'] = np.sum(count_map[mask] * exp_map_total[mask] / exp_map[mask]) / np.mean(mask) - bounds['size'] = (raw_seeds['size'] * 0.1, raw_seeds['size'] * 10) + raw_seeds["size"] = np.sum( + count_map[mask] * exp_map_total[mask] / exp_map[mask] + ) / np.mean(mask) + bounds["size"] = (raw_seeds["size"] * 0.1, raw_seeds["size"] * 10) # reorder seeds to fnc parameter order - param_fnc = list(fnc.__code__.co_varnames[:fnc.__code__.co_argcount]) - param_fnc.remove('x') - param_fnc.remove('y') + param_fnc = list(fnc.__code__.co_varnames[: fnc.__code__.co_argcount]) + param_fnc.remove("x") + param_fnc.remove("y") seeds = {key: raw_seeds[key] for key in param_fnc} x, y = np.meshgrid(centers, centers) @@ -171,12 +183,15 @@ def fit_background(self, count_map, exp_map_total, exp_map): def f(*args): return -np.sum( - log_poisson(count_map, fnc(x, y, *args) * exp_map / exp_map_total, log_factorial_count_map)) + log_poisson( + count_map, + fnc(x, y, *args) * exp_map / exp_map_total, + log_factorial_count_map, + ) + ) logger.info(f"seeds :\n{seeds}") - m = Minuit(f, - name=seeds.keys(), - *seeds.values()) + m = Minuit(f, name=seeds.keys(), *seeds.values()) for key, bound in bounds.items(): if bound is None: m.fixed[key] = True @@ -188,14 +203,25 @@ def f(*args): func = fnc(x, y, **m.values.to_dict()) * exp_map / exp_map_total func[exp_map == 0] = 1 rel_residuals = 100 * (count_map - func) / func - logger.info(f"Fit valid : {m.valid}\n" - f"Results ({fnc.__name__}) :\n{m.values.to_dict()}") - logger.info("Average relative residuals : %.1f %%," % (np.mean(rel_residuals)) + - "Std = %.2f %%" % (np.std(rel_residuals)) + "\n") + logger.info( + f"Fit valid : {m.valid}\n" + f"Results ({fnc.__name__}) :\n{m.values.to_dict()}" + ) + logger.info( + "Average relative residuals : %.1f %%," % (np.mean(rel_residuals)) + + "Std = %.2f %%" % (np.std(rel_residuals)) + + "\n" + ) return fnc(x, y, **m.values.to_dict()) - def create_acceptance_map(self, observations: Observations) -> Background3D: + def create_acceptance_map( + self, + observations: Observations, + rotate_all_to_obs: Optional[Observation] = None, + zd_correction: Optional[Background3D] = None, + derotate: Optional[bool] = False, + ) -> Background3D: """ Calculate a 3D grid acceptance map @@ -203,57 +229,99 @@ def create_acceptance_map(self, observations: Observations) -> Background3D: ---------- observations : gammapy.data.observations.Observations The collection of observations used to make the acceptance map - + rotate_all_to_obs : gammapy.data.observations.Observation + All events of `observations` are rotate in camera coordinates to match + the azimuth angle of `rotate_all_to_obs`. + For stereoscopic telescopes, like MAGIC, the morphology changes with azimuth. + Therefore the statistics can be enhacned with this instead of creating a background + model for multiple azimuth bins. Returns ------- acceptance_map : gammapy.irf.background.Background3D """ # Compute base data - count_map_background, exp_map_background, exp_map_background_total, livetime = self._create_base_computation_map( - observations) + count_map_background, exp_map_background, exp_map_background_total, livetime = ( + self._create_base_computation_map( + observations, rotate_all_to_obs, zd_correction, derotate + ) + ) # Downsample map to bkg model resolution - count_map_background_downsample = count_map_background.downsample(self.oversample_map, preserve_counts=True) - exp_map_background_downsample = exp_map_background.downsample(self.oversample_map, preserve_counts=True) - exp_map_background_total_downsample = exp_map_background_total.downsample(self.oversample_map, - preserve_counts=True) + count_map_background_downsample = count_map_background.downsample( + self.oversample_map, preserve_counts=True + ) + exp_map_background_downsample = exp_map_background.downsample( + self.oversample_map, preserve_counts=True + ) + exp_map_background_total_downsample = exp_map_background_total.downsample( + self.oversample_map, preserve_counts=True + ) # Create axis for bkg model edges = self.offset_axis.edges extended_edges = np.concatenate((-np.flip(edges), edges[1:]), axis=None) - extended_offset_axis_x = MapAxis.from_edges(extended_edges, name='fov_lon') - bin_width_x = np.repeat(extended_offset_axis_x.bin_width[:, np.newaxis], extended_offset_axis_x.nbin, axis=1) - extended_offset_axis_y = MapAxis.from_edges(extended_edges, name='fov_lat') - bin_width_y = np.repeat(extended_offset_axis_y.bin_width[np.newaxis, :], extended_offset_axis_y.nbin, axis=0) + extended_offset_axis_x = MapAxis.from_edges(extended_edges, name="fov_lon") + bin_width_x = np.repeat( + extended_offset_axis_x.bin_width[:, np.newaxis], + extended_offset_axis_x.nbin, + axis=1, + ) + extended_offset_axis_y = MapAxis.from_edges(extended_edges, name="fov_lat") + bin_width_y = np.repeat( + extended_offset_axis_y.bin_width[np.newaxis, :], + extended_offset_axis_y.nbin, + axis=0, + ) # Compute acceptance_map - if self.method == 'stack': - corrected_counts = count_map_background_downsample.data * (exp_map_background_total_downsample.data / - exp_map_background_downsample.data) - elif self.method == 'fit': + if self.method == "stack": + corrected_counts = count_map_background_downsample.data * ( + exp_map_background_total_downsample.data + / exp_map_background_downsample.data + ) + elif self.method == "fit": logger.info(f"Performing the background fit using {self.fit_fnc}.") corrected_counts = np.empty(count_map_background_downsample.data.shape) for e in range(count_map_background_downsample.data.shape[0]): - logger.info(f"Energy bin : [{self.energy_axis.edges[e]:.2f},{self.energy_axis.edges[e + 1]:.2f}]") - corrected_counts[e] = self.fit_background(count_map_background_downsample.data[e].astype(int), - exp_map_background_total_downsample.data[e], - exp_map_background_downsample.data[e], - ) + logger.info( + f"Energy bin : [{self.energy_axis.edges[e]:.2f},{self.energy_axis.edges[e + 1]:.2f}]" + ) + corrected_counts[e] = self.fit_background( + count_map_background_downsample.data[e].astype(int), + exp_map_background_total_downsample.data[e], + exp_map_background_downsample.data[e], + ) else: raise NotImplementedError(f"Requested method '{self.method}' is not valid.") - solid_angle = 4. * (np.sin(bin_width_x / 2.) * np.sin(bin_width_y / 2.)) * u.steradian - data_background = corrected_counts / solid_angle[np.newaxis, :, :] / self.energy_axis.bin_width[:, np.newaxis, - np.newaxis] / livetime - - acceptance_map = Background3D(axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], - data=data_background.to(u.Unit('s-1 MeV-1 sr-1')), - fov_alignment=FoVAlignment.ALTAZ) + solid_angle = ( + 4.0 * (np.sin(bin_width_x / 2.0) * np.sin(bin_width_y / 2.0)) * u.steradian + ) + data_background = ( + corrected_counts + / solid_angle[np.newaxis, :, :] + / self.energy_axis.bin_width[:, np.newaxis, np.newaxis] + / livetime + ) + # data_background = np.swapaxes(np.flip(data_background, 1), 1, 2) + data_background = np.swapaxes(np.flip(data_background, 2), 1, 2) + + acceptance_map = Background3D( + axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], + data=data_background.to(u.Unit("s-1 MeV-1 sr-1")), + fov_alignment=FoVAlignment.ALTAZ, + ) return acceptance_map - def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsNDMap, WcsNDMap, WcsNDMap, u.Quantity]: + def _create_base_computation_map( + self, + observations: Observations, + rotate_all_to_obs: Optional[Observation] = None, + zd_correction: Optional[Background3D] = None, + derotate: Optional[bool] = False, + ) -> Tuple[WcsNDMap, WcsNDMap, WcsNDMap, u.Quantity]: """ From a list of observations return a stacked finely binned counts and exposure map in camera frame to compute a model @@ -262,6 +330,12 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN ---------- observations : gammapy.data.observations.Observations The list of observations + rotate_all_to_obs : gammapy.data.observations.Observation + All events of `observations` are rotate in camera coordinates to match + the azimuth angle of `rotate_all_to_obs`. + For stereoscopic telescopes, like MAGIC, the morphology changes with azimuth. + Therefore the statistics can be enhacned with this instead of creating a background + model for multiple azimuth bins. Returns ------- @@ -277,7 +351,7 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN count_map_background = WcsNDMap(geom=self.geom) exp_map_background = WcsNDMap(geom=self.geom, unit=u.s) exp_map_background_total = WcsNDMap(geom=self.geom, unit=u.s) - livetime = 0. * u.s + livetime = 0.0 * u.s with erfa_astrom.set(ErfaAstromInterpolator(1000 * u.s)): for obs in observations: @@ -286,13 +360,22 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN mask = geom.contains(obs.events.radec) obs._events = obs.events.select_row_subset(~mask) # Create a count map in camera frame - camera_frame_obs = self._transform_obs_to_camera_frame(obs) - count_map_obs, _ = self._create_map(camera_frame_obs, self.geom, [], add_bkg=False) + camera_frame_obs, rot_angle = self._transform_obs_to_camera_frame( + obs, rotate_all_to_obs, derotate + ) + count_map_obs, _ = self._create_map( + camera_frame_obs, self.geom, [], add_bkg=False + ) + # Create exposure maps and fill them with the obs livetime - exp_map_obs = MapDataset.create(geom=count_map_obs.geoms['geom']) - exp_map_obs_total = MapDataset.create(geom=count_map_obs.geoms['geom']) - exp_map_obs.counts.data = camera_frame_obs.observation_live_time_duration.value - exp_map_obs_total.counts.data = camera_frame_obs.observation_live_time_duration.value + exp_map_obs = MapDataset.create(geom=count_map_obs.geoms["geom"]) + exp_map_obs_total = MapDataset.create(geom=count_map_obs.geoms["geom"]) + exp_map_obs.counts.data = ( + camera_frame_obs.observation_live_time_duration.value + ) + exp_map_obs_total.counts.data = ( + camera_frame_obs.observation_live_time_duration.value + ) # Evaluate the average exclusion mask in camera frame # by evaluating it on time intervals short compared to the field of view rotation @@ -302,15 +385,26 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN # Compute the exclusion region in camera frame for the average time dtime = time_interval[i + 1] - time_interval[i] time = time_interval[i] + dtime / 2 - average_alt_az_frame = AltAz(obstime=time, - location=obs.observatory_earth_location) - average_alt_az_pointing = obs.get_pointing_icrs(time).transform_to(average_alt_az_frame) - exclusion_region_camera_frame = self._transform_exclusion_region_to_camera_frame( - average_alt_az_pointing) + average_alt_az_frame = AltAz( + obstime=time, location=obs.meta.location + ) + average_alt_az_pointing = obs.get_pointing_icrs(time).transform_to( + average_alt_az_frame + ) + exclusion_region_camera_frame = ( + self._transform_exclusion_region_to_camera_frame( + average_alt_az_pointing, + rotate_all_to_obs, + rot_angle, + ) + ) geom_image = self.geom.to_image() - exclusion_mask_t = ~geom_image.region_mask(exclusion_region_camera_frame) if len( - exclusion_region_camera_frame) > 0 else ~Map.from_geom(geom_image) + exclusion_mask_t = ( + ~geom_image.region_mask(exclusion_region_camera_frame) + if len(exclusion_region_camera_frame) > 0 + else ~Map.from_geom(geom_image) + ) # Add the exclusion mask in camera frame weighted by the time interval duration exclusion_mask += exclusion_mask_t * (dtime).value # Normalise the exclusion mask by the full observation duration @@ -318,12 +412,94 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN # Correct the exposure map by the exclusion region for j in range(count_map_obs.counts.data.shape[0]): - exp_map_obs.counts.data[j, :, :] = exp_map_obs.counts.data[j, :, :] * exclusion_mask - + exp_map_obs.counts.data[j, :, :] = ( + exp_map_obs.counts.data[j, :, :] * exclusion_mask + ) + + # Remove Gradient with zd_correction + if zd_correction is not None: + zd_data = zd_correction.data + rot_zd_corr = rotate( + zd_data, + -rot_angle.to_value("deg"), + axes=[1, 2], + reshape=False, + order=3, + ) + # import matplotlib.pyplot as plt + # + # plt.imshow( + # rot_zd_corr[0], origin="lower", cmap="RdBu_r", vmin=-1, vmax=1 + # ) + # plt.colorbar() + # plt.show() + # plt.imshow( + # count_map_obs.counts.data[0], + # origin="lower", + # vmin=-1, + # vmax=1, + # ) + # plt.colorbar() + # plt.show() + flat = ( + count_map_obs.counts.data + - rot_zd_corr / 2 * count_map_obs.counts.data + ) + delta = np.sum(flat, 0) - rotate( + np.sum(flat, 0), 180, reshape=False + ) + # import matplotlib.pyplot as plt + # + # plt.imshow( + # np.sum(count_map_obs.counts.data, 0), + # origin="lower", + # cmap="Reds", + # ) + # plt.colorbar() + # plt.show() + # plt.imshow( + # rot_zd_corr[0], + # origin="lower", + # cmap="RdBu_r", + # vmin=-1, + # vmax=1, + # ) + # plt.colorbar() + # plt.show() + # + # plt.imshow( + # delta / np.sum(flat, 0), + # origin="lower", + # cmap="RdBu_r", + # vmin=-1, + # vmax=1, + # ) + # plt.colorbar() + # plt.show() + # plt.imshow( + # np.sum(flat, 0), + # origin="lower", + # cmap="RdBu_r", + # vmin=-1, + # vmax=1, + # ) + # plt.colorbar() + # plt.show() + + count_map_obs.counts.data = ( + count_map_obs.counts.data + - rot_zd_corr / 2 * count_map_obs.counts.data + + zd_data / 2 * count_map_obs.counts.data + ) # Stack counts and exposure maps and livetime of all observations count_map_background.data += count_map_obs.counts.data exp_map_background.data += exp_map_obs.counts.data exp_map_background_total.data += exp_map_obs_total.counts.data livetime += camera_frame_obs.observation_live_time_duration - return count_map_background, exp_map_background, exp_map_background_total, livetime + return ( + count_map_background, + exp_map_background, + exp_map_background_total, + livetime, + ) diff --git a/setup.py b/setup.py index d199c0b..ff1da02 100644 --- a/setup.py +++ b/setup.py @@ -3,23 +3,23 @@ # read the contents of your README file this_directory = path.abspath(path.dirname(__file__)) -with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f: +with open(path.join(this_directory, "README.md"), encoding="utf-8") as f: long_description = f.read() setup( - name='acceptance_modelisation', + name="acceptance_modelisation", packages=find_packages(), - version='0.2.0', - description='Calculate 2D and 3d acceptance model for creating maps with gammapy (IACT analysis)', + version="0.2.0", + description="Calculate 2D and 3d acceptance model for creating maps with gammapy (IACT analysis)", long_description=long_description, - long_description_content_type='text/markdown', + long_description_content_type="text/markdown", install_requires=[ - 'gammapy>=1.1,<1.2', - 'numpy', - 'scipy', - 'astropy>=4.0', - 'regions>=0.7,<0.8' + "gammapy==1.2", + "numpy", + "scipy", + "astropy>=4.0", + "regions>=0.7,<0.8", ], - author='Mathieu de Bony de Lavergne', - author_email='mathieu.debony@cea.fr' + author="Mathieu de Bony de Lavergne", + author_email="mathieu.debony@cea.fr", )