Skip to content

Commit

Permalink
style
Browse files Browse the repository at this point in the history
Signed-off-by: Quentin Remy <quentin.remy@mpi-hd.mpg.de>
  • Loading branch information
QRemy committed Apr 13, 2023
1 parent c8c255d commit 975a5f7
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 37 deletions.
74 changes: 41 additions & 33 deletions examples/tutorials/data/hawc.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,16 +38,16 @@
"""

import numpy as np
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
from gammapy.data import DataStore, HDUIndexTable, ObservationTable
from gammapy.datasets import MapDataset
from gammapy.estimators import ExcessMapEstimator
from gammapy.makers import MapDatasetMaker, SafeMaskMaker
from gammapy.data import DataStore, HDUIndexTable, ObservationTable
from gammapy.maps import Map, MapAxis, WcsGeom
from gammapy.estimators import ExcessMapEstimator

######################################################################
# Check setup
Expand All @@ -60,7 +60,7 @@
######################################################################
# Chose which estimator we will use

energy_estimator = 'NN'
energy_estimator = "NN"


######################################################################
Expand All @@ -85,21 +85,21 @@
# Load the tables

data_path = "$GAMMAPY_DATA/hawc/crab_events_pass4/"
hdu_filename = f'hdu-index-table-{energy_estimator}-Crab.fits.gz'
obs_filename = f'obs-index-table-{energy_estimator}-Crab.fits.gz'
hdu_filename = f"hdu-index-table-{energy_estimator}-Crab.fits.gz"
obs_filename = f"obs-index-table-{energy_estimator}-Crab.fits.gz"

# there is only one observation table
obs_table = ObservationTable.read(data_path+obs_filename)
obs_table = ObservationTable.read(data_path + obs_filename)

# we read the hdu index table of fHit bin number 6
fHit = 6
hdu_table = HDUIndexTable.read(data_path+hdu_filename, hdu=fHit)
hdu_table = HDUIndexTable.read(data_path + hdu_filename, hdu=fHit)


######################################################################
# From the tables, we can create a Datastore

data_store = DataStore( hdu_table=hdu_table, obs_table=obs_table)
data_store = DataStore(hdu_table=hdu_table, obs_table=obs_table)

data_store.info()

Expand Down Expand Up @@ -148,35 +148,40 @@
# First we define the geometry and axes

# create the energy reco axis
energy_axis = MapAxis.from_edges([1.00,1.78,3.16,5.62,10.0,17.8,31.6,56.2,100,177,316] * u.TeV,
name="energy",
interp="log")
energy_axis = MapAxis.from_edges(
[1.00, 1.78, 3.16, 5.62, 10.0, 17.8, 31.6, 56.2, 100, 177, 316] * u.TeV,
name="energy",
interp="log",
)


# and energy true axis
energy_axis_true = MapAxis.from_energy_bounds(1e-3, 1e4, nbin=140, unit="TeV", name="energy_true")
energy_axis_true = MapAxis.from_energy_bounds(
1e-3, 1e4, nbin=140, unit="TeV", name="energy_true"
)


# create a geometry around the Crab location
geom = WcsGeom.create(skydir=SkyCoord(ra=83.63,dec=22.01, unit='deg', frame="icrs"),
width=6*u.deg,
axes=[energy_axis],
binsz=0.05)
geom = WcsGeom.create(
skydir=SkyCoord(ra=83.63, dec=22.01, unit="deg", frame="icrs"),
width=6 * u.deg,
axes=[energy_axis],
binsz=0.05,
)


######################################################################
# Define the makers we will use

maker = MapDatasetMaker(selection = ["counts", "background", "exposure", "edisp", "psf"])
safemask_maker = SafeMaskMaker(methods=['aeff-max'], aeff_percent=10)
maker = MapDatasetMaker(selection=["counts", "background", "exposure", "edisp", "psf"])
safemask_maker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)


######################################################################
# Create empty Mapdataset
dataset_empty = MapDataset.create(geom,
energy_axis_true= energy_axis_true,
name ='fHit ' + str(fHit),
reco_psf=True )
dataset_empty = MapDataset.create(
geom, energy_axis_true=energy_axis_true, name="fHit " + str(fHit), reco_psf=True
)
# run the map dataset maker
dataset = maker.run(dataset_empty, obs)

Expand All @@ -196,14 +201,12 @@
# with the event list. For convenience, the HAWC data release already
# included this quantity as a map

transit_map = Map.read(data_path + 'irfs/TransitsMap_Crab.fits.gz')
transit_map = Map.read(data_path + "irfs/TransitsMap_Crab.fits.gz")
transit_number = transit_map.get_by_coord(geom.center_skydir)

# Correct the quantities
dataset.background.data*=transit_number
dataset.exposure.data*=transit_number


dataset.background.data *= transit_number
dataset.exposure.data *= transit_number


######################################################################
Expand All @@ -218,22 +221,27 @@
dataset.peek()



######################################################################
# And make significance maps to check that the Crab is visible

excess_estimator = ExcessMapEstimator(correlation_radius='0.2 deg', selection_optional = [], energy_edges=energy_axis.edges)
excess_estimator = ExcessMapEstimator(
correlation_radius="0.2 deg", selection_optional=[], energy_edges=energy_axis.edges
)
excess = excess_estimator.run(dataset)

(dataset.mask*excess['sqrt_ts']).plot_grid(add_cbar=True, cmap='coolwarm', vmin=-5, vmax=5)
(dataset.mask * excess["sqrt_ts"]).plot_grid(
add_cbar=True, cmap="coolwarm", vmin=-5, vmax=5
)
plt.show()

######################################################################
# Combining all energies
excess_estimator_integrated = ExcessMapEstimator(correlation_radius='0.2 deg', selection_optional = [])
excess_estimator_integrated = ExcessMapEstimator(
correlation_radius="0.2 deg", selection_optional=[]
)
excess_integrated = excess_estimator_integrated.run(dataset)

excess_integrated['sqrt_ts'].plot(add_cbar=True)
excess_integrated["sqrt_ts"].plot(add_cbar=True)
plt.show()


Expand Down
4 changes: 2 additions & 2 deletions gammapy/datasets/map.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ def create_map_dataset_geoms(
migra_axis=None,
rad_axis=None,
binsz_irf=None,
reco_psf=False
reco_psf=False,
):
"""Create map geometries for a `MapDataset`
Expand Down Expand Up @@ -653,7 +653,7 @@ def create(
rad_axis=rad_axis,
migra_axis=migra_axis,
binsz_irf=binsz_irf,
reco_psf=reco_psf
reco_psf=reco_psf,
)

kwargs.update(geoms)
Expand Down
5 changes: 3 additions & 2 deletions gammapy/datasets/tests/test_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
EnergyDependentMultiGaussPSF,
EnergyDispersion2D,
PSFMap,
RecoPSFMap
RecoPSFMap,
)
from gammapy.makers.utils import make_map_exposure_true_energy, make_psf_map
from gammapy.maps import HpxGeom, Map, MapAxis, RegionGeom, WcsGeom, WcsNDMap
Expand Down Expand Up @@ -1875,6 +1875,7 @@ def test_peek(images):
with mpl_plot_check():
dataset.peek()


def test_create_psf_reco(geom):
dat = MapDataset.create(geom, reco_psf=True)
assert isinstance(dat.psf, RecoPSFMap)
assert isinstance(dat.psf, RecoPSFMap)

0 comments on commit 975a5f7

Please sign in to comment.