Gammapy Overview Dortmund 2021

An overview of the package and analysis methods implemented in Gammapy, that I gave at the technical meeting of the Astroparticle Group Dortmund in January 2021.

Transcript

1. GAMMAPY Analysis methods & Package Axel Donath, Astroparticle Group Dortmund, Technical Meeting, Jan. 18th 2021 1
2. WHO AM I? • My name is Axel Donath, I am a PostDoc at MPIK Heidelberg • During my PhD I worked on the H.E.S.S. Galactic Plane Survey catalog. Out of this work the Gammapy package evolved. Now I am one of the lead-developers of Gammapy • Apologies if there is a bias to H.E.S.S jargon and data… • My Github profile: https://github.com/adonath • General interest in Python based astronomical software: co-organised Python in Gamma-Ray astronomy (e.g. PyGamma 2019), some contributions to the Astropy package • Thanks a lot for the invitation to present Gammapy to your group! 2
3. ANALYSIS METHODS Gamma-Ray Data Analysis in a Nutshell… 3
4. IACT OBSERVATION 4 “How CTA Detects Cherenkov Light”, taken from www.cta-observatory.org Basic definition
5. IACT OBSERVATION 5 “How CTA Detects Cherenkov Light”, taken from www.cta-observatory.org • Instrument is pointed for a given period of time to a position one the sky and measures events. • This is called an “observation” or “run”. • The instrument characteristics (“response”) is assumed to be stable during the whole observation run • For H.E.S.S. it’s typically ~30 min. • This is the fundamental “chunk” of data, from which the data is combined (“stacked”) or split, depending on the analysis case. • For data bookkeeping each observation is assigned a unique identifier Basic definition
6. IACT OBSERVATION 6 “How CTA Detects Cherenkov Light”, taken from www.cta-observatory.org • Instrument is pointed for a given period of time to a position one the sky and measures events. • This is called an “observation” or “run”. • The instrument characteristics (“response”) is assumed to be stable during the whole observation run • For H.E.S.S. it’s typically ~30 min. • This is the fundamental “chunk” of data, from which the data is combined (“stacked”) or split, depending on the analysis case. • For data bookkeeping each observation is assigned a unique identifier Basic definition
7. 7 Field of View (“FoV") + Pointing position Source Offset “How CTA Detects Cherenkov Light”, taken from www.cta-observatory.org “FoV coordinates” and conventions IACT OBSERVATION
8. 8 “How CTA Detects Cherenkov Light”, taken from www.cta-observatory.org Field of View (“FoV") + Source FoV - Lon FoV - Lat IACT OBSERVATION “FoV coordinates” and conventions Aligned with horizontal (“Alt-Az”) system,
 centered on pointing position Pointing position
9. DATA LEVELS 9 Reconstruction & Event Lists Telescope I Camera Telescope II Camera • Goal: estimate properties of the event: energy, arrival time, arrival direction, possibly reconstruction quality (“event class”)
10. DATA LEVELS 10 Reconstruction & Event Lists Telescope I Camera Telescope II Camera • Goal: estimate properties of the event: energy, arrival time, arrival direction, possibly reconstruction quality (“event class”) • Separation between Gamma-like and Hadron-like events • Only “Gamma-like” events are kept in the list
11. DATA LEVELS 11 Reconstruction & Event Lists Telescope I Camera Telescope II Camera • Goal: estimate properties of the event: energy, arrival time, arrival direction, possibly reconstruction quality (“event class”) • Separation between Gamma-like and Hadron-like events • Only “Gamma-like” events are kept in the list • However it is still >90% background, i.e. hadronic (or electron) events “mis-classified” as Gamma events (requires additional background estimation, comes later… ) • Reconstruction of events using various methods: template fitting, “Hillas Parameters”, experts are in this group…
12. DATA LEVELS 12 Reconstruction & Event Lists Telescope I Camera Telescope II Camera • Goal: estimate properties of the event: energy, arrival time, arrival direction, possibly reconstruction quality (“event class”) • Separation between Gamma-like and Hadron-like events • Only “Gamma-like” events are kept in the list • However it is still >90% background, i.e. hadronic (or electron) events “mis-classified” as Gamma events (requires additional background estimation, comes later… ) • Reconstruction of events using various methods: template fitting, “Hillas Parameters”, experts are in this group…
13. IRFS 13 • Effective area: gives number of expected (true) gamma-like events per area • Varies with offset from pointing position and true energy • Estimated from “Monte Carlo” simulations (again experts are in the group…) • Required to measure flux / brightness of Gamma-Rays sources Effective Area DL3 DR1
14. IRFS 14 Point Spread Function • Point spread function (PSF): angular resolution of the instrument, precision to reconstruct the arrival direction of an event • Estimated from “Monte Carlo” simulations and by binning events into offset and true energy • Typically stored as “radial profile” and varies with offset from pointing position and true energy • There are parametric models of the shape as well e.g. triple Gaussian or King profile, obtained by fitting the shape to the event distributions • Required to measure extension of Galactics sources & precise flux of point sources DL3 DR1
15. IRFS 15 Energy Dispersion • Energy dispersion: accuracy and precision to reconstruct the energy of an event • Varies with offset from pointing position and true energy • Often “migration” I used: ratio defined as ETrue / E • Required to measure precise spectra of source, especially at low energies (for IACTS) DL3 DR1
16. IRFS • Background model (sometime called “acceptance”…) derived from many “OFF” observations • Instrument is pointed to an “empty” (without known gamma-ray sources..) region in the sky • Done many times observations are grouped with similar observation conditions, mostly number of telescopes and bands of zenith angle • Taken from: L.Mohrmann et al. https://ui.adsabs.harvard.edu/abs/2019arXiv191008088M 16 Residual Hadronic Background FoV-Lon (deg) FoV-Lat (deg) FoV-Lon (deg) FoV-Lat (deg) DL3 DR1
17. DL3 DATA FORMATS • Events and IRFs are stored “per observation” in a dedicated FITS format “GADF” • Openly developed and documented data formats for Gamma-Ray astronomy • Initiative after the PyGamma15 conference • Based on the “Fermi-LAT” formats used for data distribution • https://gamma-astro-data-formats.readthedocs.io • https://github.com/open-gamma-ray-astro/gamma-astro-data-formats 17 A standardised FITS data format
18. DATA REDUCTION 18 Observation and / or time selection Bin events into ND sky maps • Selection of observation “ids” and / or time ranges Observation / Observations
19. DATA REDUCTION 19 Observation and / or time selection Bin events into ND sky maps • Selection of observation “ids” and / or time ranges Observation / Observations = Class names in Gammapy
20. DATA REDUCTION 20 Energy Lon Lat Observation and / or time selection Bin events into ND sky maps Bin selection: WCS & Energy • Selection of observation “ids” and / or time ranges • Most general case: selection of spatial & energy binning WcsGeom / MapAxis Observation / Observations
21. DATA REDUCTION 21 Energy Lon Lat Counts (NObs ) Energy Lon Lat Observation and / or time selection Bin selection: WCS & Energy Bin events into ND sky maps • Selection of observation “ids” and / or time ranges • Most general case: selection of spatial & energy binning WcsNDMap WcsGeom / MapAxis Observation / Observations Skymap / “Cube”
22. DATA REDUCTION 22 Energy Lon Lat Observation and / or time selection Bin selection: WCS & 1 Energy Bin Lon Lat Counts (NObs ) Bin events into ND sky maps • Image analysis is handled as a “cube” with one energy bin WcsNDMap WcsGeom / MapAxis Observation / Observations Skymap
23. DATA REDUCTION 23 Observation and / or time selection Bin selection: Region & Energy icrs;circle(0, 0, 0.1) Energy Counts Energy Lon Lat Bin events into ND sky maps Counts (NObs ) • Spectral analysis is handled as a “cube” with one spatial bin RegionNDMap RegionGeom / MapAxis Observation / Observations Spectrum
24. DATA REDUCTION 24 True Energy Lon Lat Exposure
 (eff. area x lifetime ) Reproject IRFs to ND sky maps True Energy Lon Lat Pointing position + Offset • Interpolate effective area on the chosen spatial and true energy binning • Multiply with the (dead time corrected) observation time, to get an exposure map EffectiveAreaTable2D WcsNDMap / RegionNDMap Bin selection: WCS & Energy WcsGeom / MapAxis
25. DATA REDUCTION 25 Energy Dispersion Matrix True Energy Energy Lon Lat 4th Dimension: Energy True Energy Energy dispersion Pointing position + Offset • Energy dispersion varies with offset: if multiple sources are fitted during the analysis a “map” of RMF matrices is required • Typically computed on a WCS with coarser bins (~0.2 deg) • During fitting / model evaluation the corresponding matrix for a given model component is extracted from the map and applied using matrix multiplication EDispKernelMap EDispKernel EnergyDispersion2D Bin selection: Coarse WCS & True Energy / Energy
26. DATA REDUCTION 26 Energy True Lon Lat 4th Dimension: Rad Rad PSF PSF Model Point spread function Pointing position + Offset • PSF varies with offset: if multiple sources are fitted during the analysis a “map” of PSFs is required • Typically computed on a WCS with coarser bins (~0.2 deg) • During fitting / model evaluation the corresponding PSF the current position of a model component is extracted from the map • From the radial profile a 3D “kernel” is computed and applied using convolution True Energy Lon Lat PSF Kernel PSFMap PSFKernel PSF3D Bin selection: Coarse WCS, True Energy & Rad
27. DATA REDUCTION 27 Background Energy Lon Lat Pointing position FoV-Lon (deg) FoV-Lat (deg) FoV-Lon FoV-Lat Bkg Template Energy Lon Lat • Background rate (events / s / sr / TeV) is reprojected on the sky and multiplied with solid angle, observation time and integrated in energy • For the later analysis the “template” is combined with a parametric correction such as norm and spectral ”tilt” WcsNDMap / RegionNDMap Bin selection: WCS & Energy WcsGeom / MapAxis Background3D / Background2D
28. “CLASSICAL” BKG • “Reflected regions”: used for spectral analysis, find “off regions” reflected from the pointing position, ignore exclusion regions • “Ring background”: mostly for image analysis, derive “off counts” per pixel from ring convolved counts, ignore exclusion regions. There is a variant called “adaptive ring”, while enlarges the ring to achieve better statistics. • See e.g. Berge et.al. https://ui.adsabs.harvard.edu/abs/2007A%26A...466.1219B/abstract 28 Measure “Off” counts from data ReflectedRegionsBackgroundMaker RingBackgroundMaker
29. DATA MODEL 29 Compute predicted counts NPred,Src = EDISPSrc (PSFSrc (ℰSrc ⋅ fSrc )) True Energy Lon Lat Source Model True Energy Lon Lat PSF Kernel True Energy Lon Lat Exposure
 (eff. area x lifetime ) Energy Dispersion Matrix True Energy Energy Sometimes called “foward folding” as well… fSrc = fSpectral (E) ⋅ fSpatial (E, l, b) ⋅ fTemporal (t)
30. CASH LIKELIHOOD 30 𝒞 = 2∑ i Ni Pred − Ni Obs ⋅ log Ni Pred “Cash” statistics: summed over all “bins” NPred = NBkg + ∑ Src NPred,Src “Fermi” style analysis Energy Lon Lat Counts NObs = Bkg Template Energy Lon Lat NBkg = ⋅ n0 ⋅ ( E E0 ) −Γ FoVBackgroundModel • Predicted counts (as defined on the previous slide) are computed per model component (“source”) and summed • A “global” background model with “correction parameters” (like Fermi-LAT analysis) is added • Summed per bin fit statistics= - 2 log(L) computed and optimised
31. ON-OFF LIKELIHOOD 31 “Classical” IACT analysis 𝒲 = 2∑ i (Ni Sig + (1 + 1/α)Ni Bkg − Ni On log(Ni Sig + Ni Bkg ) − Ni Off log(Ni Bkg /α)) “WStat” statistics: summed over all “bins” Background is “marginalised” NSig = ∑ i NPred,Src ∂𝒲 ∂NBkg = 0 Energy Lon Lat Counts NOn = Energy Lon Lat Counts Off NOff = !
32. JOINT LIKELIHOOD 32 Combining datasets … ℒ(x1 , θ) * ℒ(x2 , θ) ℒ(xn−1 , θ) ℒ(xn , θ) * ℒ(xn−2 , θ) * Model True Energy Lon Lat Source Model Exposure Map
 EDisp Kernel Map Counts Map PSF Map
 Background Map MapDataset / MapDatasetOnOff SpectrumDataset / SpectrumDatasetOnOff OBS=1 OBS=2
33. Parameter value Best fit value MORE STATISTICS 33 Parameter value Best fit value UL Parameter value Best fit value Neg. Error Pos. Error Asymmetric errors Delta TS = 1 { Upper limits Likelihood profile Delta TS = 4 • “Significance” based on model comparison: and (valid if the difference in degrees of freedom between the models is 1, Wilk’s theorem) • Errors and upper limits based on likelihood profiles, with or without “marginalisation” i.e. re- optimisation of other parameters, which are not of interested (often background norms, norms of other sources in the FoV,…) Significance, errors and upper limits Δ𝖳𝖲 = 𝒞0 − 𝒞1 σ = Δ𝖳𝖲
34. FLUX POINTS 34 Grouping energy bins Energy (TeV) Counts • Method based on the Fermi-LAT catalogs, no “unfolding” but rescaling a reference model (usually the best fit model) in a selected energy range FluxPointsEstimator
35. FLUX POINTS 35 Grouping energy bins Energy (TeV) Counts 10 TeV 0.3 TeV 1 TeV 3TeV • Method based on the Fermi-LAT catalogs, no “unfolding” but rescaling a reference model (usually the best fit model) in a selected energy range FluxPointsEstimator
36. FLUX POINTS 36 Grouping energy bins Energy (TeV) Counts 10 TeV 0.3 TeV 1 TeV 3TeV • Method based on the Fermi-LAT catalogs, no “unfolding” but rescaling a reference model (usually the best fit model) in a selected energy range FluxPointsEstimator
37. FLUX POINTS 37 Grouping energy bins Energy (TeV) Counts 10 TeV 0.3 TeV 3TeV 1 TeV • Method based on the Fermi-LAT catalogs, no “unfolding” but rescaling a reference model (usually the best fit model) in a selected energy range FluxPointsEstimator
38. FLUX POINTS 38 Grouping energy bins • Method based on the Fermi-LAT catalogs, no “unfolding” but rescaling a reference model (usually the best fit model) in a selected energy range • The fitting is done using a “norm” parameter , which scales the flux with respect to the reference model. Reference values for dnde (“differential flux”), flux (“integral flux”) and eflux (“energy flux”) are stored along with norm, wich allows to convert the flux points to a different SED type later Energy (TeV) Counts 10 TeV 0.3 TeV 3TeV 1 TeV FluxPointsEstimator
39. PACKAGE A Python Package Gamma-Ray Astronomy… 39
40. DEPENDENCIES 40 Coordinates, Quantities, Tables,
 FITS I/O, etc. Interpolation, minimisation,
 FFT convolution, etc. ND-data structures
 and computations Required dependencies Core dependencies: scientific Python stack
41. DEPENDENCIES 41 Coordinates, Quantities, Tables,
 FITS I/O, etc. Interpolation, minimisation,
 FFT convolution, etc. Command line tools pyyaml YAML I/O click Required dependencies ND-data structures
 and computations pydantic Configuration Required dependencies: make developer live a bit simpler
42. DEPENDENCIES 42 Coordinates, Quantities, Tables,
 FITS I/O, etc. Interpolation, minimisation,
 FFT convolution, etc. Plotting, visualisation Iminuit Command line tools healpy Healpix maps pydantic Configuration pyyaml YAML I/O click Optimisation, sampling Optional dependencies Required dependencies ND-data structures
 and computations Tutorial notebooks Optional dependencies: bring in useful functionality
43. GAMMAPY SUB-PACKAGE STRUCTURE 43 Mostly ordered by data level and transition between data levels
44. Coveralls used for monitoring
 of code test coverage DEVELOPMENT AND CI SETUP 44 Hosted and openly developed on Github:
 https://github.com/gammapy/gammapy Travis-CI and Azure Pipelines used
 for continuous integration Docs are build and deployed manually: https://docs.gammapy.org/

45. Coveralls used for monitoring
 of code test coverage DEVELOPMENT AND CI SETUP 45 Hosted and openly developed on Github:
 https://github.com/gammapy/gammapy Travis-CI and Azure Pipelines used
 for continuous integration Black code formatting used
 for a consistent code format. Sphinx used to build the
 documentation Pytest used for testing Docs are build and deployed manually: https://docs.gammapy.org/

46. RELATED PROJECTS 46 https://github.com/gammapy/gamma-sky gamma-cat “An open data collection
 and source catalog for VHE gamma-ray astronomy” gammasky.net “A portal to the gamma-ray sky” PyGamma 15 / 19 “Python and open data for
 Gamma-Ray Astronomy
 Workshop”, March 18th - 22nd Gamma-astro-data-formats “A place to propose and share data format descriptions for gamma-ray astronomy.” FITS https://github.com/open-gamma-ray-astro https://indico.cern.ch/event/783425/ https://github.com/gammapy/gamma-cat
47. INSTALLATION & SETUP 47 $ curl -O https://gammapy.org/download/install/gammapy-0.18.2-environment.yml $ conda env create -f gammapy-0.18.2-environment.yml $ conda activate gammapy-0.18.2 1. 2. 3. Install Gammapy: https://docs.gammapy.org/0.18.2/install/index.html
48. INSTALLATION & SETUP 48 Install Gammapy: Download tutorials: $ gammapy download tutorials $ cd gammapy-tutorials $ export GAMMAPY_DATA=$PWD/datasets 1. 2. 3. https://docs.gammapy.org/0.18.2/install/index.html $ curl -O https://gammapy.org/download/install/gammapy-0.18.2-environment.yml $ conda env create -f gammapy-0.18.2-environment.yml $ conda activate gammapy-0.18.2 1. 2. 3.
49. TUTORIALS • Best starting point are the tutorials notebooks: https://docs.gammapy.org/0.18.2/tutorials/ index.html • HESS data access: https://docs.gammapy.org/0.18.2/notebooks/cta.html • Spectral analysis: https://docs.gammapy.org/0.18.2/tutorials/spectrum_analysis.html • Lightcurve computation: https://docs.gammapy.org/0.18.2/tutorials/light_curve_flare.html • And many, many others… • Tutorials can be executed online as well (using Binder): 49 Click here!
50. DEMO You can execute the notebook along with me or just watch….
51. CONTACT AND SUPPORT • Slack: gammapy.slack.com (quick questions, immediate help,
 the @tu-dortmund.de domain can join without invitation…) • Github issues: https://github.com/gammapy/gammapy/issues
 (feature requests & bug reports) • Gammapy mailing list (not used much…): • gammapy@googlegroups.com • Archive: https://groups.google.com/forum/#!forum/gammapy • Gammapy dev meetings (technical meetings in case you’d like to contribute): • Every Friday 10am on Zoom (direct link in the Slack #dev channel) • Friday is Gammapy co-working day • Regular Gammapy CTA and non CTA user calls, will be announced on the CTA mailing list and Gammapy Slack 51