Coverage for lst_auto_rta/acceptanceCalculation.py: 0%

1import astropy.units as u 

2import numpy as np 

3from astropy.coordinates import SkyCoord 

4from gammapy.data import DataStore 

5from gammapy.datasets import MapDataset 

6from gammapy.irf import Background2D 

7from gammapy.makers import MapDatasetMaker 

8from gammapy.maps import WcsGeom, WcsNDMap 

9from regions import CircleAnnulusSkyRegion 

10 

11 

12def create_radial_acceptance_map(observations, energy_axis, offset_axis, oversample_map=10, exclude_regions=[]): 

13 """ 

14 Calculate a radial acceptance map 

15 

16 Parameters 

17 ---------- 

18 observations : gammapy.data.observations.Observations 

19 The collection of observations used to make the acceptance map 

20 energy_axis : gammapy.maps.geom.MapAxis 

21 The energy axis for the acceptance model 

22 offset_axis : gammapy.maps.geom.MapAxis 

23 The offset axis for the acceptance model 

24 oversample_map : int 

25 oversample in number of pixel of the spatial axis used for the calculation 

26 exclude_regions : list of 'regions.SkyRegion' 

27 Region with sources, will be excluded of the calculation of the acceptance map 

28 

29 Returns 

30 ------- 

31 background : gammapy.irf.background.Background2D 

32 A bakground model that could be used as an acceptance model 

33 """ 

34 

35 n_bins_map = offset_axis.nbin * oversample_map * 2 

36 binsz = offset_axis.edges[-1] / (n_bins_map / 2) 

37 center_map = SkyCoord(ra=0.0 * u.deg, dec=0.0 * u.deg, frame="icrs") 

38 

39 geom = WcsGeom.create( 

40 skydir=center_map, npix=(n_bins_map, n_bins_map), binsz=binsz, frame="icrs", axes=[energy_axis] 

41 ) 

42 

43 count_map_background = WcsNDMap(geom=geom) 

44 exp_map_background = WcsNDMap(geom=geom, unit=u.s) 

45 exp_map_background_total = WcsNDMap(geom=geom, unit=u.s) 

46 livetime = 0.0 * u.s 

47 

48 for obs in observations: 

49 geom_obs = WcsGeom.create( 

50 skydir=obs.pointing_radec, npix=(n_bins_map, n_bins_map), binsz=binsz, frame="icrs", axes=[energy_axis] 

51 ) 

52 count_map_obs = MapDataset.create(geom=geom_obs) 

53 exp_map_obs = MapDataset.create(geom=geom_obs) 

54 exp_map_obs_total = MapDataset.create(geom=geom_obs) 

55 maker = MapDatasetMaker(selection=["counts"]) 

56 

57 exclusion_mask = geom_obs.to_image().region_mask(exclude_regions, inside=False) 

58 exclusion_map = WcsNDMap(geom_obs.to_image(), exclusion_mask) 

59 

60 count_map_obs = maker.run(MapDataset.create(geom=geom_obs), obs) 

61 

62 exp_map_obs.counts.data = obs.observation_live_time_duration.value 

63 exp_map_obs_total.counts.data = obs.observation_live_time_duration.value 

64 

65 for i in range(count_map_obs.counts.data.shape[0]): 

66 count_map_obs.counts.data[i, :, :] = count_map_obs.counts.data[i, :, :] * exclusion_map 

67 exp_map_obs.counts.data[i, :, :] = exp_map_obs.counts.data[i, :, :] * exclusion_map 

68 

69 count_map_background.data += count_map_obs.counts.data 

70 exp_map_background.data += exp_map_obs.counts.data 

71 exp_map_background_total.data += exp_map_obs_total.counts.data 

72 livetime += obs.observation_live_time_duration 

73 

74 # import matplotlib.pyplot as plt 

75 # plt.figure() 

76 # count_map_background.plot_grid(figsize=(10, 10))#, add_cbar=True) 

77 # plt.figure() 

78 # exp_map_background.plot_grid(figsize=(10, 10), add_cbar=True) 

79 # plt.figure() 

80 # exp_map_background_total.plot_grid(figsize=(10, 10), add_cbar=True) 

81 # plot_x = [] 

82 # plot_y = [] 

83 

84 data_background = np.zeros((energy_axis.nbin, offset_axis.nbin)) * u.Unit("s-1 MeV-1 sr-1") 

85 for i in range(offset_axis.nbin): 

86 selection_region = CircleAnnulusSkyRegion( 

87 center=center_map, inner_radius=offset_axis.edges[i], outer_radius=offset_axis.edges[i + 1] 

88 ) 

89 selection_mask = geom.to_image().region_mask([selection_region], inside=True) 

90 selection_map = WcsNDMap(geom.to_image(), selection_mask) 

91 for j in range(energy_axis.nbin): 

92 value = u.dimensionless_unscaled * np.sum(count_map_background.data[j, :, :] * selection_map) 

93 value *= np.sum(exp_map_background_total.data[j, :, :] * selection_map) / np.sum( 

94 exp_map_background.data[j, :, :] * selection_map 

95 ) 

96 

97 # plt.figure() 

98 # plt.imshow(count_map_background.data[j, :, :]*selection_map) 

99 

100 value /= energy_axis.edges[j + 1] - energy_axis.edges[j] 

101 value /= ( 

102 2.0 

103 * np.pi 

104 * (offset_axis.edges[i + 1].to("radian") - offset_axis.edges[i].to("radian")) 

105 * offset_axis.center[i].to("radian") 

106 ) 

107 value /= livetime 

108 data_background[j, i] = value 

109 

110 # print(np.sum(count_map_background.data[j, :, :]*selection_map), np.sum(count_map_background.data[j, :, :]*selection_map)*np.sum(exp_map_background_total.data[j, :, :]*selection_map)/np.sum(exp_map_background.data[j, :, :]*selection_map), (energy_axis.edges[j+1]-energy_axis.edges[j]), 2. * np.pi * (offset_axis.edges[j+1].to('radian')-offset_axis.edges[j].to('radian')) * offset_axis.center[j].to('radian'), livetime, value, np.sum(exp_map_background_total.data[j, :, :]*selection_map)/np.sum(exp_map_background.data[j, :, :]*selection_map), np.sum(exp_map_background_total.data[j, :, :]*selection_map), np.sum(exp_map_background.data[j, :, :]*selection_map)) 

111 

112 # plot_x.append(np.sum(exp_map_background.data[j, :, :]*selection_map)) 

113 # plot_y.append((2. * np.pi * (offset_axis.edges[i+1].to('radian')-offset_axis.edges[i].to('radian')) * offset_axis.center[i].to('radian')).value) 

114 

115 # plt.figure() 

116 # plt.plot(plot_x, plot_y, '+') 

117 

118 background = Background2D(energy_axis=energy_axis, offset_axis=offset_axis, data=data_background) 

119 

120 return background