Progress report

My work on the last 3 weeks has been mainly in the container class for the cube background models (X, Y, energy). The class is called CubeBackgroundModel and the code has recently been merged to the master branch of Gammapy.

The class has remodeled after a few code reviews from its first draft as in the post on Friday, June 19, 2015. For instance it can read/write 2 different kind of FITS formats:

• FITS binary tables: more convenient for storing and data analysis.
• FITS images: more convenient for the visualization, using for instance DS9.

For the records, FITS is a standard data format largely used in astronomy.

In addition, the plotting methods have been also simplified to allow a more customizable API for the user. Now only one plot is returned by the methods, and the user can easily combine the plots as desired with only a few lines of code using matplotlib.

A new function has been added to the repository as well for creating dummy background cube models called make_test_bg_cube_model. This function creates a background following a 2D symmetric Gaussian model for the spatial coordinates (X, Y) and a power-law in energy. The Gaussian width varies in energy from sigma/2 to sigma. An option is also available to mask 1/4th of the Gaussian images. This option will be useful in the future, when testing the still-to-come reprojection methods, necessary for applying the background model to the analysis data to subtract the background. Since the models are produced in the detector coordinate system (a.k.a. nominal system), the models need to be projected to sky coordinates (i.e. Galactic, or RA/Dec) in order to apply them to the data.

The work on the CubeBackgroundModel class has also triggered the development of other utility functions, for instance to create WCS coordinate objects for describing detector coordinates in FITS format or a converter of Astropy Table objects to FITS binary table ones.

Moreover, a test file with a dummy background cube produced with the make_test_bg_cube_model tool has been placed in the gammapy-extra repository here for testing the input/output (read/write) methods of the class.

This work has also triggered some discussions about some methods and classes in both the Astropy and Gammapy repositories. As a matter of fact, I am currently solving some of them, especially for the preparation of the release of the Gammapy 0.3 stable version in the coming weeks.

In parallel I am also currently working on a script that should become a command-line program to produce background models using the data of a given gamma-ray astronomy experiment. The script is still on a first draft version, but the idea is to have a program that:

1. looks for the data (all observations of a given experiment)
2. filters out the observations taken on known sources
3. divides the data into groups of similar observation conditions
4. creates the background models and stores them to file

In order to create the model, the following steps are necessary:

• stack events and bin then (fill a histogram)
• apply livetime correction
• apply bin volume correction
• smooth histogram (not yet implemented)

A first glimpse on such a background model is shown in the following animated image (please click on the animation for an enlarged view):

The movie shows a sequence of 4 images (X, Y), one for each energy bin slice of the cube. The image spans 10 deg on each direction, and the energy binning is defined between 0.01 TeV and 100 TeV, equidistant in logarithmic scale. The model is performed for a zenith angle range between 0 deg and 20 deg.

There is still much work to do in order to polish the script and move most of the functionality into Gammapy classes and functions, until the script is only a few high-level calls to the necessary methods in the correct order.

Mid-term summary

Mid-term has arrived and quite some work has been done for Gammapy, especially in the observation, dataset and background modules. At the same time I have learnt a lot about Gammapy, Astropy (especially tables, quantities, angles, times and fits files handling), and python (especially numpy and matplotlib.pyplot). But the most useful thing I'm learning is to produce good code via code reviews. The code review process is sometimes hard and frustrating, but very necessary in order to produce clear code that can be read and used by others.

The last week I have been working on a method to filter observations tables as the one presented in the figure on the first report. The method is intended to be used to select observations according to different criteria (for instance data quality, or within a certain region in the sky) that should be used for a particular analysis.

In the case of background modeling this is important to separate observations taken on or close to known sources or far from them. In addition, the observations can be grouped according to similar observation conditions. For instance observations taken under a similar zenith angle. This parameter is very important in gamma-ray observations.

The zenith angle of the telescopes is defined as the angle between the vertical (zenith) and the direction where the telescopes are pointing. The smaller the zenith angle is, the more vertical the telescopes are pointing, and the thinner is the atmosphere layer. This has large consequences in the amount and properties of the gamma-rays detected by the telescopes. Gamma-rays interact in the upper atmosphere and produce Cherenkov light, which is detected by the telescopes. The amount of light produced is directly proportional to the energy of the gamma-ray. In addition, the light is emitted in a narrow cone along the direction of the gamma-ray.

At lower zenith angles the Cherenkov light has to travel a smaller distance through the atmosphere, so there is less absorption. This means that lower energy gamma-rays can be detected.

At higher zenith angles the Cherenkov light of low-energy gamma-rays is totally absorbed, but the Cherenkov light cones of the high-energy ones are longer, and hence the section of ground covered is larger, so particles that fall further away from the telescopes can be detected, increasing the amount of detected high-energy gamma-rays.

The zenith angle is maybe the most important parameter, when grouping the observations in order to produce models of the background.

The method implemented can filter the observations according to this (and other) parameters. An example using a dummy observation table generated with the tool presented on the first report is presented here (please click on the picture for an enlarged view):

Please notice that instead of the mentioned zenith angle, altitude as the zenith's complementary angle (altitude_angle = 90 deg - zenith_angle) is used.
In this case, the first table was generated with random altitude angles between 45 deg and 90 deg (or 0 deg to 45 deg in zenith), while the second table is filtered to keep only zenith angles in the range of 20 deg to 30 deg (or 60 deg to 70 deg in altitude).

The tool can be used to apply selections in any variable present in the observation table. In addition, an 'inverted' flag has been programmed in order o be able to apply the filter to keep the values outside the selection range, instead of inside.

Recapitulating the progress done until now, the next steps will be to finish the tools that I am implementing now: the filter observations method described before and the background cube model class on the previous report. In both cases there is still some work to do: an inline application for filtering observations and more methods to create cube background models.

The big milestone is to have a working chain to produce cube background models from existing event lists within a couple of weeks.

Progress report

The last two weeks have been a bit tough. After finishing the observation table generator mentioned in the previous post, I started working on the background module of Gammapy.

I am currently working on a container class for 3D background models (a.k.a. cube background models). The three dimensions of the model are detector coordinates (X and Y) and energy. These kind of background models are largely in use by Fermi. The development of tools for creating such background models is an important milestone in my project and so implementing a class to handle them is crucial. For now, the class can read cube models from fits files and slice the 3D models to produce 2D plots of the background rate.

Two kinds of plots are produced:

1. Color maps of the background rate for each energy slice defined in the model.
   
2. 1D curves of the spectrum of the background rate for each detector bin (X, Y) defined in the model.
   

The figures attached have been performed using a sample background model fits file from the Gammalib repository. (Please click for an enlarged view).

More functionality should come soon into this class, for instance, methods to create the bg models using event lists and the smoothing of the models to attenuate the effects of the statistical nature of the event detection. As I mentioned, I had some trouble developing a more complex class in python and this task is taking more time than expected.I am working hard to keep on track.

First report

Wow! Two weeks are already gone?! Time goes fast, when working on an interesting project!

After a few iterations of the API with my mentors, I started working on the observation and datataset modules of Gammapy. These modules are essential for observation bookkeeping and data access. Since the background subtraction methods should act on the data, a good implementation of these modules is important for developing the toolbox mentioned in my proposal.

Before explaining what I am working on exactly, I will try to clarify how TeV gamma-ray observations are carried and organized.

A typical observation would be to point the telescope system (for instance H.E.S.S. or CTA) to a particular position in the sky (i.e. a source like the Crab nebula) for a certain amount of time (typically ~30 min.). Since most sources in gamma-ray need deep exposures (typically several hours to ~100 h) to obtain a significant emission, and hence actually "see something", an object is observed during multiple observations of ~30 min each.

Each of these observations delivers a large amount of data (raw data): telescope images of interesting events (plus also a lot of noise), atmosphere monitoring information, telescope calibration information, etc. Each experiment or observatory has a team of experts that filter the data and derive the properties of the primary particles (ideally gamma-rays, but also a lot of gamma-like background events) and produce the so called event list files: one file per observation.

These files contain ordered list of events with their properties (arrival direction, energy, etc.) and a header with the information common to all of them, for being in the same observation (observation time, telescope pointing direction, etc.). These event list files are the data, that is delivered to the observer/analyzer and are the input for the analysis tools.

In order to keep track of these observations and organize the event lists needed for a specific analysis, observation tables are used. These are but a list of observations IDs that should identify the corresponding event lists, together with interesting properties that help the analyzer in the selection of the observations needed for a specific analysis. These properties are for instance the observation time, the pointing position, the observation duration, the quality of the data obtained, etc.

My work on the last 2 weeks has been especially in the development of a tool that generates dummy observation tables that I will be able to use for testing some of the tools that I need to develop for the background module of Gammapy. In addition, I am working together with my mentors on defining the data format that Gammapy observation tables will support. A link to the pull request can be found here.

The work is not yet finished, but the tool generates already test observation tables with a few columns. Here is a screenshot with a 10 observation example table (please click on the picture for an enlarged view):

The algorithm of the tool generates random pairs of azimuth and altitude evenly distributed in spherical coordinates. Azimuth values are generated in the interval (0 deg, 360 deg), and altitude values in the interval (45 deg, 90 deg). It also generates random time values between the start of 2010 and the end of 2014. The az, alt pairs are then translated into celestial coordinates (RA, dec) for a specific observatory. More on coordinates and its transformations can be found in the Astropy docs here, or in a less technical way in Wikipedia here.

This tool could also be used in the future to generate example observation lists.

Items that need to be improved:

• Write the doc specifying the supported format for observation tables in Gammapy.
• Add additional columns to the dummy observation table generator.
• Add a header to the observation table with for instance the name of the observatory/experiment (i.e. H.E.S.S. or CTA).
• The times in the observation table generator should be written in a more convenient format.
• Optional: since TeV gammas-ray observations are carried at night, it would be nice if the randomly generated observation times could be restricted to night hours (i.e. from 22h to 4h).
• Optional: develop a tool to filter observation tables depending on a few parameters like:
• Time interval: define t_min, t_max
• Pointing position: define a circle/box in the sky.

GSoC bonding period

The bonding period of the Google Summer of Code is on going, and I have been in touch with the mentors of my project about once a week.

During this time (and also during the application process) I have been getting familiarized with the tools needed for the project: install an appropriate Conda environment, learn how to use git and Github, install the development versions of Astropy and Gammapy and make some pull requests for merging some code, add some documentations, and correct some (minor) bugs.

In addition, together with my mentors, I have started working on the API (application programming interface) for the toolbox for background modeling for Gammapy. The approach is to start writing down some typical use cases and from there write a high-level pseudo-code that should work with minimal modifications once the necessary tools are implemented.

One use case is for instance: having a list of off run event lists from a IACT system (such as H.E.S.S. or CTA), develop the necessary tools to divide the list of runs into specific bins according to the observation properties (for instance altitude and azimuth angles of the IACT system) and build bg templates for each bin.

Project description

Gamma-ray astronomy has experienced a fast development in the past 2 decades with both ground-based imaging atmospheric Cherenkov telescope (IACT) experiments like H.E.S.S., MAGIC and VERITAS (ref: ) and satellites like Fermi. In addition the next generation of IACT experiment CTA is in its prototype phase. 

The instruments developed for detecting gamma-rays accumulate many background events. The majority is rejected either by using an intelligent triggering system for the detectors, or in early stages of the data analysis. Unfortunately, there is still a dominant fraction of events passing the gamma-ray selection cuts, called the gamma-like background. In order to extract the gamma-ray signal, clever algorithms for modeling the gamma-like background are essential.

During the GSoC 2015 I will be working for the Gammapy project. Gammapy is an open source (BSD licensed) gamma-ray astronomy Python package. It is an in-development affiliated package of Astropy, a community effort to develop a single core package for Astronomy in Python. Gammapy builds on the core scientific Python stack to provide tools to simulate and analyze the gamma-ray sky for telescopes such as Fermi, H.E.S.S. and CTA.

Specifically I will implement the most successful background modeling methods largely in use by the gamma-ray community in the Astropy/Gammapy framework.
The background methods can be classified in two categories, according to the observation strategy:

1. Background models from OFF observations, where the background is modeled from observations far away from any known sources. IACT experiments require dedicated OFF-source observations for modeling the background.
2.  Background models from ON observations, where the background is modeled using observation within or close-by to the region of interest (a.k.a. ON region).

For more details about background modeling you can read my proposal or the paper by Berge 2007.
In a first step I will implement tools to create background model templates from observations with no or only a few gamma-ray sources in the field of view. In a second step, I will develop algorithms to estimate the background in observations containing gamma-ray sources to detect them and measure their spatial shape and energy spectrum, in some cases using the model templates from the first step.

GSoC application

Link to my application in Astropy wiki.