Coverage for src / lstautorta / acceptanceCalculation.py: 0%

1import astropy.units as u 

2import numpy as np 

3from astropy.coordinates import SkyCoord 

4from gammapy.datasets import MapDataset 

5from gammapy.irf import Background2D 

6from gammapy.makers import MapDatasetMaker 

7from gammapy.maps import WcsGeom, WcsNDMap 

8from regions import CircleAnnulusSkyRegion 

9 

10 

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

12 """ 

13 Calculate a radial acceptance map 

14 

15 Parameters 

16 ---------- 

17 observations : gammapy.data.observations.Observations 

18 The collection of observations used to make the acceptance map 

19 energy_axis : gammapy.maps.geom.MapAxis 

20 The energy axis for the acceptance model 

21 offset_axis : gammapy.maps.geom.MapAxis 

22 The offset axis for the acceptance model 

23 oversample_map : int 

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

25 exclude_regions : list of 'regions.SkyRegion' 

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

27 

28 Returns 

29 ------- 

30 background : gammapy.irf.background.Background2D 

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

32 """ 

33 n_bins_map = offset_axis.nbin * oversample_map * 2 

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

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

36 

37 geom = WcsGeom.create( 

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

39 ) 

40 

41 count_map_background = WcsNDMap(geom=geom) 

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

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

44 livetime = 0.0 * u.s 

45 

46 for obs in observations: 

47 geom_obs = WcsGeom.create( 

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

49 ) 

50 count_map_obs = MapDataset.create(geom=geom_obs) 

51 exp_map_obs = MapDataset.create(geom=geom_obs) 

52 exp_map_obs_total = MapDataset.create(geom=geom_obs) 

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

54 

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

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

57 

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

59 

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

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

62 

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

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

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

66 

67 count_map_background.data += count_map_obs.counts.data 

68 exp_map_background.data += exp_map_obs.counts.data 

69 exp_map_background_total.data += exp_map_obs_total.counts.data 

70 livetime += obs.observation_live_time_duration 

71 

72 # import matplotlib.pyplot as plt 

73 # plt.figure() 

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

75 # plt.figure() 

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

77 # plt.figure() 

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

79 # plot_x = [] 

80 # plot_y = [] 

81 

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

83 for i in range(offset_axis.nbin): 

84 selection_region = CircleAnnulusSkyRegion( 

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

86 ) 

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

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

89 for j in range(energy_axis.nbin): 

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

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

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

93 ) 

94 

95 # plt.figure() 

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

97 

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

99 value /= ( 

100 2.0 

101 * np.pi 

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

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

104 ) 

105 value /= livetime 

106 data_background[j, i] = value 

107 

108 # 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)) 

109 

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

111 # 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) 

112 

113 # plt.figure() 

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

115 

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

117 

118 return background