From bc51bf7a53e32b2df6710586f7f2aaa24fa0d237 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 12 Nov 2024 20:00:22 +0100 Subject: [PATCH 01/28] rot events, more statistics --- acceptance_modelisation/__init__.py | 2 +- .../base_acceptance_map_creator.py | 893 ++++++++++++------ .../grid3d_acceptance_map_creator.py | 302 +++--- setup.py | 24 +- 4 files changed, 823 insertions(+), 398 deletions(-) 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..472d33d 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -20,35 +20,38 @@ 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, +) 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 +106,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 +141,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 +154,9 @@ 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 + ) -> Observation: """ Transform events and pointing of an obs from a sky frame to camera frame @@ -146,6 +164,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 ------- @@ -154,36 +178,66 @@ def _transform_obs_to_camera_frame(obs: Observation) -> Observation: """ # 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) + # Rotation to one observation + if rotate_to_obs is not None: + frame_centers = pointing_altaz.transform_to(camera_frame) + + sep = frame_centers.separation(events_camera_frame) + pos_angle = frame_centers.position_angle(events_camera_frame) + + # mid Az of rotate_to_obs + az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az + rot_angle = az_obs - pointing_altaz.az + + 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 + 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 return obs_camera_frame - def _transform_exclusion_region_to_camera_frame(self, pointing_altaz: AltAz) -> List[SkyRegion]: + def _transform_exclusion_region_to_camera_frame( + self, + pointing_altaz: AltAz, + rotate_to_obs: Optional[Observation] = None, + ) -> List[SkyRegion]: """ Transform the list of exclusion regions in sky frame into a list in camera frame. @@ -191,6 +245,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 +263,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 + - pointing_altaz.az + ) + 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 + rot_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 +366,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 +403,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 +432,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 +516,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 +545,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 +605,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 +664,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 +731,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 +761,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 +792,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 +811,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 +884,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 +915,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 +991,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 +1028,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 +1108,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 +1137,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 +1172,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 +1251,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..b5fb4ca 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -5,7 +5,7 @@ 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 @@ -19,29 +19,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 +103,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 +119,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 +148,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 +161,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 +182,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 +202,23 @@ 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, + ) -> Background3D: """ Calculate a 3D grid acceptance map @@ -203,57 +226,93 @@ 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) + ) # 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 + ) + + 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, + ) -> 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 +321,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 +342,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 +351,21 @@ 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 = self._transform_obs_to_camera_frame( + obs, rotate_all_to_obs + ) + 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 +375,25 @@ 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, + ) + ) 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,7 +401,9 @@ 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 + ) # Stack counts and exposure maps and livetime of all observations count_map_background.data += count_map_obs.counts.data @@ -326,4 +411,9 @@ def _create_base_computation_map(self, observations: Observations) -> Tuple[WcsN 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", ) From 1ef7369a7c515a1acef93aa1f69e04662abc1b59 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 00:02:56 +0100 Subject: [PATCH 02/28] rot dir flip --- acceptance_modelisation/base_acceptance_map_creator.py | 4 ++-- acceptance_modelisation/grid3d_acceptance_map_creator.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 472d33d..bafba75 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -209,10 +209,10 @@ def _transform_obs_to_camera_frame( # mid Az of rotate_to_obs az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az - rot_angle = az_obs - pointing_altaz.az + rot_angle = az_obs + pointing_altaz.az events_camera_frame = frame_centers.directional_offset_by( - position_angle=pos_angle + rot_angle, separation=sep + position_angle=pos_angle - rot_angle, separation=sep ) # Formatting data for the output diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index b5fb4ca..99d372a 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -299,6 +299,7 @@ def create_acceptance_map( / self.energy_axis.bin_width[:, np.newaxis, np.newaxis] / livetime ) + # data_background = np.swapaxes(data_background, 1, 2) acceptance_map = Background3D( axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], From be42b5f0bdc0ce1dbaee72655c576b25ff6ab128 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 00:18:50 +0100 Subject: [PATCH 03/28] rot dir flip --- acceptance_modelisation/base_acceptance_map_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index bafba75..1b6ba63 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -289,7 +289,7 @@ def _transform_exclusion_region_to_camera_frame( center_coordinate_camera_frame ) center_coordinate_camera_frame = frame_center.directional_offset_by( - position_angle=pos_angle + rot_angle, separation=sep + position_angle=pos_angle - rot_angle, separation=sep ) center_coordinate_camera_frame_arb = SkyCoord( ra=center_coordinate_camera_frame.lon[0], From 311e02c9a24898b8a9d9d3b87fd645227b35ce07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 00:37:05 +0100 Subject: [PATCH 04/28] rot dir flip --- acceptance_modelisation/base_acceptance_map_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 1b6ba63..36c511f 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -209,7 +209,7 @@ def _transform_obs_to_camera_frame( # mid Az of rotate_to_obs az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az - rot_angle = az_obs + pointing_altaz.az + rot_angle = az_obs - pointing_altaz.az events_camera_frame = frame_centers.directional_offset_by( position_angle=pos_angle - rot_angle, separation=sep From 7cc451605e028cba6f017532d4f2bd7a990e1d9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 10:42:42 +0100 Subject: [PATCH 05/28] pca --- acceptance_modelisation/base_acceptance_map_creator.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 36c511f..7a27f6b 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -211,6 +211,13 @@ def _transform_obs_to_camera_frame( az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az rot_angle = az_obs - pointing_altaz.az + pca = PCA(n_components=2).fit( + np.stack([events_camera_frame.lon, events_camera_frame.lat]).T + ) + component0_x, component0_y = pca.components_[0] + angle = np.rad2deg(np.arctan(component0_y / component0_x)) + rot_angle = az_obs - angle * u.deg + events_camera_frame = frame_centers.directional_offset_by( position_angle=pos_angle - rot_angle, separation=sep ) From de30c2de69acad183a5f2e74ddaef822c15fad96 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 10:55:44 +0100 Subject: [PATCH 06/28] pca --- acceptance_modelisation/base_acceptance_map_creator.py | 1 + 1 file changed, 1 insertion(+) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 7a27f6b..5a5b35c 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -27,6 +27,7 @@ generate_irf_from_mini_irf, compute_neighbour_condition_validation, ) +from sklearn.decomposition import PCA logger = logging.getLogger(__name__) From b4fc209150e6bd08a3c8382caa5ec68848314d8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 12:18:29 +0100 Subject: [PATCH 07/28] fix pca angle exclusion --- acceptance_modelisation/base_acceptance_map_creator.py | 7 ++++--- acceptance_modelisation/grid3d_acceptance_map_creator.py | 3 ++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 5a5b35c..829d6be 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 @@ -239,12 +239,13 @@ def _transform_obs_to_camera_frame( ) obs_camera_frame._location = obs.meta.location - return obs_camera_frame + return obs_camera_frame, angle * u.deg 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. @@ -290,7 +291,7 @@ def _transform_exclusion_region_to_camera_frame( frame_center = pointing_altaz.transform_to(camera_frame) rot_angle = ( rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az - - pointing_altaz.az + - pca_angle ) sep = frame_center.separation(center_coordinate_camera_frame) pos_angle = frame_center.position_angle( diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 99d372a..61f8a94 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -352,7 +352,7 @@ def _create_base_computation_map( 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( + camera_frame_obs, pca_angle = self._transform_obs_to_camera_frame( obs, rotate_all_to_obs ) count_map_obs, _ = self._create_map( @@ -386,6 +386,7 @@ def _create_base_computation_map( self._transform_exclusion_region_to_camera_frame( average_alt_az_pointing, rotate_all_to_obs, + pca_angle, ) ) geom_image = self.geom.to_image() From 387976a0d1e300683cb3537ac86e944618495ffe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 14:23:12 +0100 Subject: [PATCH 08/28] fix rot angel --- .../base_acceptance_map_creator.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 829d6be..69df95a 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -216,11 +216,12 @@ def _transform_obs_to_camera_frame( np.stack([events_camera_frame.lon, events_camera_frame.lat]).T ) component0_x, component0_y = pca.components_[0] - angle = np.rad2deg(np.arctan(component0_y / component0_x)) - rot_angle = az_obs - angle * u.deg + derot_angle = -np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg + theory_angle = az_obs - 34.23 * u.deg + rot_angle = theory_angle - derot_angle events_camera_frame = frame_centers.directional_offset_by( - position_angle=pos_angle - rot_angle, separation=sep + position_angle=pos_angle + rot_angle, separation=sep ) # Formatting data for the output @@ -238,8 +239,10 @@ def _transform_obs_to_camera_frame( aeff=obs.aeff, ) obs_camera_frame._location = obs.meta.location + if rotate_to_obs is not None: + return obs_camera_frame, rot_angle - return obs_camera_frame, angle * u.deg + return obs_camera_frame, None def _transform_exclusion_region_to_camera_frame( self, From cd8d86714754b0daffb29920c26cbe26be0975f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 14:26:59 +0100 Subject: [PATCH 09/28] fix region --- acceptance_modelisation/base_acceptance_map_creator.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 69df95a..1ad6a0a 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -292,16 +292,16 @@ def _transform_exclusion_region_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 - ) + # 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 - rot_angle, separation=sep + position_angle=pos_angle + pca_angle, separation=sep ) center_coordinate_camera_frame_arb = SkyCoord( ra=center_coordinate_camera_frame.lon[0], From 73704a7a24e74a4ea658d787ce59af15b4cb2579 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 14:38:31 +0100 Subject: [PATCH 10/28] fix region 180 --- acceptance_modelisation/base_acceptance_map_creator.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 1ad6a0a..0882134 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -219,6 +219,8 @@ def _transform_obs_to_camera_frame( derot_angle = -np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg theory_angle = az_obs - 34.23 * u.deg rot_angle = theory_angle - derot_angle + while rot_angle > 90 * u.deg: + rot_angle = rot_angle - 180 * u.deg events_camera_frame = frame_centers.directional_offset_by( position_angle=pos_angle + rot_angle, separation=sep From ebb36467d9b8a497b21a60b553e072944695c366 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 16:30:17 +0100 Subject: [PATCH 11/28] pca --- .../base_acceptance_map_creator.py | 46 ++++++++++++++++--- 1 file changed, 39 insertions(+), 7 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 0882134..dae1448 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -203,6 +203,36 @@ def _transform_obs_to_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) + 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 + + ##### frame_centers = pointing_altaz.transform_to(camera_frame) sep = frame_centers.separation(events_camera_frame) @@ -210,20 +240,22 @@ def _transform_obs_to_camera_frame( # mid Az of rotate_to_obs az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az - rot_angle = az_obs - pointing_altaz.az 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 + derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg theory_angle = az_obs - 34.23 * u.deg - rot_angle = theory_angle - derot_angle - while rot_angle > 90 * u.deg: - rot_angle = 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) events_camera_frame = frame_centers.directional_offset_by( - position_angle=pos_angle + rot_angle, separation=sep + position_angle=pos_angle - rot_angle, separation=sep ) # Formatting data for the output @@ -303,7 +335,7 @@ def _transform_exclusion_region_to_camera_frame( center_coordinate_camera_frame ) center_coordinate_camera_frame = frame_center.directional_offset_by( - position_angle=pos_angle + pca_angle, separation=sep + position_angle=pos_angle - pca_angle, separation=sep ) center_coordinate_camera_frame_arb = SkyCoord( ra=center_coordinate_camera_frame.lon[0], From 5fc2d8873ae96fe47e447db8dc96c315cd91d920 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 14 Nov 2024 18:22:35 +0100 Subject: [PATCH 12/28] flip --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 61f8a94..817db96 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -299,7 +299,7 @@ def create_acceptance_map( / self.energy_axis.bin_width[:, np.newaxis, np.newaxis] / livetime ) - # data_background = np.swapaxes(data_background, 1, 2) + data_background = np.swapaxes(data_background, 1, 2) acceptance_map = Background3D( axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], From 041d44b3e7c2fc78602e845c681167c35745b53b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 19 Nov 2024 00:51:23 +0100 Subject: [PATCH 13/28] rot correction --- acceptance_modelisation.egg-info/PKG-INFO | 273 ++++++++++++++++++ acceptance_modelisation.egg-info/SOURCES.txt | 16 + .../dependency_links.txt | 1 + acceptance_modelisation.egg-info/requires.txt | 5 + .../top_level.txt | 1 + .../base_acceptance_map_creator.py | 11 +- 6 files changed, 305 insertions(+), 2 deletions(-) create mode 100644 acceptance_modelisation.egg-info/PKG-INFO create mode 100644 acceptance_modelisation.egg-info/SOURCES.txt create mode 100644 acceptance_modelisation.egg-info/dependency_links.txt create mode 100644 acceptance_modelisation.egg-info/requires.txt create mode 100644 acceptance_modelisation.egg-info/top_level.txt 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/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index dae1448..186cc4a 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -246,8 +246,15 @@ def _transform_obs_to_camera_frame( ) component0_x, component0_y = pca.components_[0] derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg - theory_angle = az_obs - 34.23 * u.deg - rot_angle = derot_angle_obs - derot_angle + az = obs.get_pointing_altaz(obs.tmid).az + theory_diff = az_obs - az + rot_angle = -derot_angle_obs + derot_angle + if theory_diff > 90 * u.deg: + rot_angle += 180 * u.deg + if 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: From bde8735938e29ac1c251cc5f4562d52b1c3eedca Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 19 Nov 2024 09:12:14 +0100 Subject: [PATCH 14/28] rot 180 --- acceptance_modelisation/base_acceptance_map_creator.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 186cc4a..185421b 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -248,11 +248,12 @@ def _transform_obs_to_camera_frame( derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg az = obs.get_pointing_altaz(obs.tmid).az theory_diff = az_obs - az - rot_angle = -derot_angle_obs + derot_angle - if theory_diff > 90 * u.deg: + rot_angle = derot_angle_obs - derot_angle + print(rot_angle) + if np.abs(theory_diff) > 90 * u.deg: rot_angle += 180 * u.deg - if theory_diff < -90 * u.deg: - rot_angle -= 180 * u.deg + print(theory_diff.deg) + print(rot_angle) # rot_angle = derot_angle_obs - derot_angle # print(derot_angle_obs, derot_angle) From f47da2f0fdf4b06a0453c114a19f2749577acc97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Wed, 20 Nov 2024 18:55:23 +0100 Subject: [PATCH 15/28] rem 180 --- acceptance_modelisation/base_acceptance_map_creator.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 185421b..07ee96c 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -250,10 +250,11 @@ def _transform_obs_to_camera_frame( theory_diff = az_obs - az rot_angle = derot_angle_obs - derot_angle print(rot_angle) - if np.abs(theory_diff) > 90 * u.deg: - rot_angle += 180 * u.deg + # if np.abs(theory_diff) > 90 * u.deg: + # rot_angle += 180 * u.deg print(theory_diff.deg) print(rot_angle) + print(derot_angle) # rot_angle = derot_angle_obs - derot_angle # print(derot_angle_obs, derot_angle) @@ -262,6 +263,7 @@ def _transform_obs_to_camera_frame( # 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 ) From 5dbd97db4e35efeec98b16e213ea23d4553a4c72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 01:07:49 +0100 Subject: [PATCH 16/28] zd corr --- .../grid3d_acceptance_map_creator.py | 33 +++++++++++++++++-- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 817db96..21c57ca 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -218,6 +218,7 @@ def create_acceptance_map( self, observations: Observations, rotate_all_to_obs: Optional[Observation] = None, + zd_correction: Optional[Background3D] = None, ) -> Background3D: """ Calculate a 3D grid acceptance map @@ -239,7 +240,9 @@ def create_acceptance_map( # Compute base data count_map_background, exp_map_background, exp_map_background_total, livetime = ( - self._create_base_computation_map(observations, rotate_all_to_obs) + self._create_base_computation_map( + observations, rotate_all_to_obs, zd_correction + ) ) # Downsample map to bkg model resolution @@ -313,6 +316,7 @@ def _create_base_computation_map( self, observations: Observations, rotate_all_to_obs: Optional[Observation] = None, + zd_correction: Optional[Background3D] = None, ) -> 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 @@ -352,12 +356,35 @@ def _create_base_computation_map( mask = geom.contains(obs.events.radec) obs._events = obs.events.select_row_subset(~mask) # Create a count map in camera frame - camera_frame_obs, pca_angle = self._transform_obs_to_camera_frame( + camera_frame_obs, rot_angle = self._transform_obs_to_camera_frame( obs, rotate_all_to_obs ) count_map_obs, _ = self._create_map( camera_frame_obs, self.geom, [], add_bkg=False ) + + # Remove Gradient with zd_correction + if zd_correction is not None: + from scipy.ndimage import rotate + + count_map_obs.counts.data = ( + count_map_obs.counts.data + - rotate( + zd_correction.data / 2, + rot_angle.to_value("deg"), + axes=[2, 1], + reshape=False, + order=1, + ) + * count_map_obs.counts.data + + zd_correction.data / 2 * count_map_obs.counts.data + ) + # import matplotlib.pyplot as plt + # + # plt.imshow(count_map_obs.counts.data[0]) + # plt.colorbar() + # plt.show() + # 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"]) @@ -386,7 +413,7 @@ def _create_base_computation_map( self._transform_exclusion_region_to_camera_frame( average_alt_az_pointing, rotate_all_to_obs, - pca_angle, + rot_angle, ) ) geom_image = self.geom.to_image() From ebf675f8f8df31b0eb5d290853cf5434231c916b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 11:55:37 +0100 Subject: [PATCH 17/28] zd corr --- .../grid3d_acceptance_map_creator.py | 23 ++++++++----------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 21c57ca..9420876 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -11,6 +11,7 @@ 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 @@ -365,25 +366,19 @@ def _create_base_computation_map( # Remove Gradient with zd_correction if zd_correction is not None: - from scipy.ndimage import rotate + rot_zd_corr = rotate( + zd_correction.data, + rot_angle.to_value("deg"), + axes=[2, 1], + reshape=False, + order=1, + ) count_map_obs.counts.data = ( count_map_obs.counts.data - - rotate( - zd_correction.data / 2, - rot_angle.to_value("deg"), - axes=[2, 1], - reshape=False, - order=1, - ) - * count_map_obs.counts.data + - rot_zd_corr * count_map_obs.counts.data + zd_correction.data / 2 * count_map_obs.counts.data ) - # import matplotlib.pyplot as plt - # - # plt.imshow(count_map_obs.counts.data[0]) - # plt.colorbar() - # plt.show() # Create exposure maps and fill them with the obs livetime exp_map_obs = MapDataset.create(geom=count_map_obs.geoms["geom"]) From 34186c3d071c653e07a068ca5914ae4a59171c19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 15:30:45 +0100 Subject: [PATCH 18/28] zd corr factor 2 --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 9420876..3175475 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -376,7 +376,7 @@ def _create_base_computation_map( count_map_obs.counts.data = ( count_map_obs.counts.data - - rot_zd_corr * count_map_obs.counts.data + - rot_zd_corr / 2 * count_map_obs.counts.data + zd_correction.data / 2 * count_map_obs.counts.data ) From 969c15da2b54d01e1669109034cb8f9b746af338 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 16:52:36 +0100 Subject: [PATCH 19/28] flip rot --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 3175475..24ac427 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -368,7 +368,7 @@ def _create_base_computation_map( if zd_correction is not None: rot_zd_corr = rotate( zd_correction.data, - rot_angle.to_value("deg"), + -rot_angle.to_value("deg"), axes=[2, 1], reshape=False, order=1, From 3dadbccef71dc25e5d04ea873eb1389896f5de94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 18:10:19 +0100 Subject: [PATCH 20/28] order=3 --- .../grid3d_acceptance_map_creator.py | 20 +++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 24ac427..0b5de24 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -368,11 +368,27 @@ def _create_base_computation_map( if zd_correction is not None: rot_zd_corr = rotate( zd_correction.data, - -rot_angle.to_value("deg"), + rot_angle.to_value("deg"), axes=[2, 1], reshape=False, - order=1, + 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( + # zd_correction.data[0], + # origin="lower", + # cmap="RdBu_r", + # vmin=-1, + # vmax=1, + # ) + # plt.colorbar() + # plt.show() count_map_obs.counts.data = ( count_map_obs.counts.data From 3e24993b1fdcc4311e14583cd458821f34a29ac7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 19:08:23 +0100 Subject: [PATCH 21/28] zd corr * exp --- .../grid3d_acceptance_map_creator.py | 66 ++++++++++--------- 1 file changed, 34 insertions(+), 32 deletions(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 0b5de24..ef978ae 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -364,38 +364,6 @@ def _create_base_computation_map( camera_frame_obs, self.geom, [], add_bkg=False ) - # Remove Gradient with zd_correction - if zd_correction is not None: - rot_zd_corr = rotate( - zd_correction.data, - rot_angle.to_value("deg"), - axes=[2, 1], - 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( - # zd_correction.data[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_correction.data / 2 * count_map_obs.counts.data - ) - # 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"]) @@ -445,6 +413,40 @@ def _create_base_computation_map( exp_map_obs.counts.data[j, :, :] * exclusion_mask ) + # Remove Gradient with zd_correction + if zd_correction is not None: + zd_data = ( + zd_correction.data * exp_map_obs.data / exp_map_obs_total.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( + # zd_correction.data[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 From 2b89b71b28353f33e397b25bb39cace9405c4afd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 19:10:38 +0100 Subject: [PATCH 22/28] counts.data --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index ef978ae..b683b43 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -416,7 +416,9 @@ def _create_base_computation_map( # Remove Gradient with zd_correction if zd_correction is not None: zd_data = ( - zd_correction.data * exp_map_obs.data / exp_map_obs_total.data + zd_correction.data + * exp_map_obs.counts.data + / exp_map_obs_total.counts.data ) rot_zd_corr = rotate( zd_data, From 591b1da2aac18fba77c6e6572d723152a5280569 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Thu, 21 Nov 2024 19:36:41 +0100 Subject: [PATCH 23/28] no exp --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index b683b43..dbf261d 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -415,11 +415,7 @@ def _create_base_computation_map( # Remove Gradient with zd_correction if zd_correction is not None: - zd_data = ( - zd_correction.data - * exp_map_obs.counts.data - / exp_map_obs_total.counts.data - ) + zd_data = zd_correction.data rot_zd_corr = rotate( zd_data, rot_angle.to_value("deg"), From a83629660c8db7fcb6f7ce606de34d9ffc20b95c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 26 Nov 2024 12:02:18 +0100 Subject: [PATCH 24/28] try solving rot definition --- .../base_acceptance_map_creator.py | 29 ++++++++-- .../grid3d_acceptance_map_creator.py | 53 +++++++++++++++++-- 2 files changed, 73 insertions(+), 9 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 07ee96c..0b10f8d 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -156,7 +156,9 @@ def __init__( @staticmethod def _transform_obs_to_camera_frame( - obs: Observation, rotate_to_obs: Optional[Observation] = None + 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 @@ -178,6 +180,18 @@ def _transform_obs_to_camera_frame( 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 + 270 * u.deg + # Transform to altaz frame altaz_frame = AltAz(obstime=obs.events.time, location=obs.meta.location) events_altaz = obs.events.radec.transform_to(altaz_frame) @@ -231,6 +245,10 @@ def _transform_obs_to_camera_frame( ) 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) @@ -238,17 +256,18 @@ def _transform_obs_to_camera_frame( sep = frame_centers.separation(events_camera_frame) pos_angle = frame_centers.position_angle(events_camera_frame) - # mid Az of rotate_to_obs - az_obs = rotate_to_obs.get_pointing_altaz(rotate_to_obs.tmid).az - 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 + derot_angle = np.rad2deg(np.arctan2(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(rot_angle) # if np.abs(theory_diff) > 90 * u.deg: # rot_angle += 180 * u.deg diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index dbf261d..3bb250c 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -220,6 +220,7 @@ def create_acceptance_map( 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 @@ -242,7 +243,7 @@ def create_acceptance_map( # Compute base data count_map_background, exp_map_background, exp_map_background_total, livetime = ( self._create_base_computation_map( - observations, rotate_all_to_obs, zd_correction + observations, rotate_all_to_obs, zd_correction, derotate ) ) @@ -318,6 +319,7 @@ def _create_base_computation_map( 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 @@ -358,7 +360,7 @@ def _create_base_computation_map( obs._events = obs.events.select_row_subset(~mask) # Create a count map in camera frame camera_frame_obs, rot_angle = self._transform_obs_to_camera_frame( - obs, rotate_all_to_obs + obs, rotate_all_to_obs, derotate ) count_map_obs, _ = self._create_map( camera_frame_obs, self.geom, [], add_bkg=False @@ -418,7 +420,7 @@ def _create_base_computation_map( zd_data = zd_correction.data rot_zd_corr = rotate( zd_data, - rot_angle.to_value("deg"), + -rot_angle.to_value("deg"), axes=[1, 2], reshape=False, order=3, @@ -431,7 +433,50 @@ def _create_base_computation_map( # plt.colorbar() # plt.show() # plt.imshow( - # zd_correction.data[0], + # 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, From f836457df03fbbc0a07679c5d9820eefa1141606 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 26 Nov 2024 12:21:49 +0100 Subject: [PATCH 25/28] try solving rot definition 2 --- .../base_acceptance_map_creator.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 0b10f8d..5269a46 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -260,7 +260,7 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): np.stack([events_camera_frame.lon, events_camera_frame.lat]).T ) component0_x, component0_y = pca.components_[0] - derot_angle = np.rad2deg(np.arctan2(component0_y, component0_x)) * u.deg + 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) @@ -268,12 +268,13 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): rot_angle = derot_angle_obs - derot_angle if derotate: rot_angle = -derot_angle - print(rot_angle) + print( + f"derot_angle_obs, derot_angle {derot_angle_obs.to_value(u.deg), derot_angle.to_value(u.deg)}" + ) + 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 - print(theory_diff.deg) - print(rot_angle) - print(derot_angle) # rot_angle = derot_angle_obs - derot_angle # print(derot_angle_obs, derot_angle) From 9df47a0c3e020513faa45c7b9a5184fb5f35fef9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 26 Nov 2024 14:49:39 +0100 Subject: [PATCH 26/28] - angle --- .../base_acceptance_map_creator.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index 5269a46..b08baaf 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -185,12 +185,14 @@ def magic_theory_angle(az: Angle): 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: + print(f"Gamma PCA: {angle}@{az}") + print(f"Gamma theory: {gamma_theory}") + 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 + 270 * u.deg + 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 - 270 * u.deg # Transform to altaz frame altaz_frame = AltAz(obstime=obs.events.time, location=obs.meta.location) @@ -261,6 +263,8 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): ) component0_x, component0_y = pca.components_[0] derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg + + print(component0_x, component0_y) az = obs.get_pointing_altaz(obs.tmid).az derot_angle = gamma_pca_corr_magic(derot_angle, az) From 4839c9e90c226125ccb31ad1f4f7fbe3f855eee4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Tue, 26 Nov 2024 23:05:45 +0100 Subject: [PATCH 27/28] reverse lon like gammapy --- .../base_acceptance_map_creator.py | 33 +++++++++++-------- .../grid3d_acceptance_map_creator.py | 2 +- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/acceptance_modelisation/base_acceptance_map_creator.py b/acceptance_modelisation/base_acceptance_map_creator.py index b08baaf..a95b6cf 100644 --- a/acceptance_modelisation/base_acceptance_map_creator.py +++ b/acceptance_modelisation/base_acceptance_map_creator.py @@ -185,14 +185,12 @@ def magic_theory_angle(az: Angle): def gamma_pca_corr_magic(angle: Angle, az: Angle): gamma_theory = magic_theory_angle(az) - print(f"Gamma PCA: {angle}@{az}") - print(f"Gamma theory: {gamma_theory}") - if gamma_theory < 90 * u.deg and gamma_theory >= 0 * u.deg: + 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 - 270 * u.deg + 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.meta.location) @@ -216,6 +214,11 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): * 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: @@ -242,6 +245,12 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): * 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 ) @@ -264,7 +273,6 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): component0_x, component0_y = pca.components_[0] derot_angle = np.rad2deg(np.arctan(component0_y / component0_x)) * u.deg - print(component0_x, component0_y) az = obs.get_pointing_altaz(obs.tmid).az derot_angle = gamma_pca_corr_magic(derot_angle, az) @@ -272,9 +280,6 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): rot_angle = derot_angle_obs - derot_angle if derotate: rot_angle = -derot_angle - print( - f"derot_angle_obs, derot_angle {derot_angle_obs.to_value(u.deg), derot_angle.to_value(u.deg)}" - ) print(f"Rot_angle{rot_angle}") print(f"Theory angle {theory_diff.deg}") # if np.abs(theory_diff) > 90 * u.deg: @@ -308,7 +313,7 @@ def gamma_pca_corr_magic(angle: Angle, az: Angle): ) obs_camera_frame._location = obs.meta.location if rotate_to_obs is not None: - return obs_camera_frame, rot_angle + return obs_camera_frame, rot_angle # , component0_x, component0_y return obs_camera_frame, None @@ -369,10 +374,10 @@ def _transform_exclusion_region_to_camera_frame( center_coordinate_camera_frame ) center_coordinate_camera_frame = frame_center.directional_offset_by( - position_angle=pos_angle - pca_angle, separation=sep + position_angle=pos_angle + pca_angle, separation=sep ) center_coordinate_camera_frame_arb = SkyCoord( - ra=center_coordinate_camera_frame.lon[0], + ra=-center_coordinate_camera_frame.lon[0], dec=center_coordinate_camera_frame.lat[0], ) exclude_region_camera_frame.append( diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 3bb250c..75b43d5 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -304,7 +304,7 @@ def create_acceptance_map( / self.energy_axis.bin_width[:, np.newaxis, np.newaxis] / livetime ) - data_background = np.swapaxes(data_background, 1, 2) + # data_background = np.swapaxes(data_background, 1, 2) acceptance_map = Background3D( axes=[self.energy_axis, extended_offset_axis_x, extended_offset_axis_y], From f7b04a51d5e1f2d3b185f0e8d051d1af02522bbe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20Fr=C3=B6se?= Date: Wed, 27 Nov 2024 10:10:38 +0100 Subject: [PATCH 28/28] flip --- acceptance_modelisation/grid3d_acceptance_map_creator.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/acceptance_modelisation/grid3d_acceptance_map_creator.py b/acceptance_modelisation/grid3d_acceptance_map_creator.py index 75b43d5..136c569 100644 --- a/acceptance_modelisation/grid3d_acceptance_map_creator.py +++ b/acceptance_modelisation/grid3d_acceptance_map_creator.py @@ -304,7 +304,8 @@ def create_acceptance_map( / self.energy_axis.bin_width[:, np.newaxis, np.newaxis] / livetime ) - # data_background = np.swapaxes(data_background, 1, 2) + # 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],