GammaLib

I converted the SDA, SDB and ICT files that are needed to compute the IAQs from the native COMPASS format into FITS and added these files to the calibration database.

Added GCOMD1Response and GCOMD2Response classes to deal with SDA and SDB datasets.

Note that geometrical smearing is not implemented since this neads a NAG routine. I’m not sure that this is very standard.

The continuum response has also been implemented, allowing to generate a response for any arbitrary law:

void set(const GModelSpectral& spectrum);

Here the list of things by class that are still missing:

• GCOMIaq::GCOMIaq: test input argument validity
• GCOMIaq::set(const GEnergy& energy): implement geometrical smearing
• GCOMIaq::set(const GModelSpectral& spectrum): implement geometrical smearing
• GCOMIaq::compute_iaq_bin: study impact of integration precision (needs full D2 response)
• GCOMIaq::klein_nishina_integral: check this formula

• GCOMInstChars::prob_D2inter: verify formula
• GCOMInstChars::prob_multihit: Implement method. Need to find out what exactly is computed.
• GCOMInstChars::atten_D1: Implement method
• GCOMInstChars::atten_D2: Implement method
• GCOMInstChars::atten_selfveto: Implement method
• GCOMInstChars::multi_scatter: Implement method

• GCOMD2Response::operator(): Add non photo peak response components
• GCOMD2Response::update_cache: Assure normalization

• Use GCOMIaq class

Advancing. Here the new list:

• GCOMIaq::GCOMIaq: test input argument validity
• GCOMIaq::set(const GEnergy& energy): implement geometrical smearing
• GCOMIaq::set(const GModelSpectral& spectrum): implement geometrical smearing
• GCOMIaq::compute_iaq_bin: study impact of integration precision (needs full D2 response)
• GCOMIaq::klein_nishina_integral: check this formula

• GCOMInstChars::multi_scatter: Implement method

• GCOMD2Response::operator(): Add Compton response components
• GCOMD2Response::update_cache: Assure normalization

• Use GCOMIaq class

A big step forward, I finally managed to implement the D2 response. Here a visualisation of the D2 response at different energies. This figure can be compared to Fig. 2.11 in Rob van Dijk’s PhD thesis. Note that the 20 MeV shape is a bit different since Rob has apparently used the input energy and not the photo peak position for the computation of the Klein-Nishina Compton tail.

The plot has be generated using the respd2.py script.

Note that so far energy thresholds are neither implement for the D1, nor the D2 response. However, response normalisation is now assured.

Got it working. Here now a couple of ratio plots of the IAQ I computed using the GCOMIaq class and the UNH simulated IAQ that is included in GammaLib. The plots are from left to right for the standard energy bins: 0.75-1 MeV, 1-3 MeV, 3-10 MeV, 10-30 MeV. From top to bottom the number of IAQ’s computed for the continuum response has double, and is 10, 20, 40, 80.

For all energy bands the computation stabilises from 40 bins on (for higher energies even earlier), so I will select for the moment 50 bins as default.

There is basically no impact of the integration precision in GCOMIaq::compute_iaq_bin on the results. Here the integral probability and computing time as function of precision, and below some ratio plots. Note that the integral probability of the UNH response is 0.044088.

I therefore keep the precision at 1e-4.

And here a nice comparison summary plot, generated using iaqcmp.py:

I now implemented the zenith-angle dependent location smearing into the code. Below a few comparison plots for zenith angles of 0°, 20°, 30°, 40° (top left, top right, bottom left, bottom right). A zenith angle of 30° seems to provide the best match with the UNH IRFs. I use this value as default.

I implemented the last missing correction which was the multi-scatter correction. Here now the final results:

And here the summary of the response integral for the four energy band, comparing UNH and GammaLib IAQs. The largest uncertainty exists for the lowest energy range.

Here the list of remaining things to be done:

• GCOMIaq::GCOMIaq: test input argument validity

• Use GCOMIaq class

I redid the analysis cookbook analysis with these new IAQs (see http://cta.irap.omp.eu/ctools/users/tutorials/howto/comptel/howto_comptel_fitting.html).

Below the ctlike results. The prefactor is about the same, the index a bit steeper:

Below the comparison of the spectra. On the left the spectrum derived with the UNH response files, on the right the GammaLib IAQs. The main difference in the spectral points is for the first energy bin that is higher (due to the -32% difference in IAQ normalisation as derived above).

I should double check the code to see whether everything is okay.

I cross-checked extensively the multi-scattering code and it seems to be correct. Here the essential functions that are used:

Here all the other probabilities that come into play in the computation of the IAQs:

It turned out that I misinterpreted the ICT data. The table show in fact the linear attenuation coefficient for the D1 & D2 modules, not the mass attenuation. The interaction probabilities are shown below.

The Phigeo and Phibar shapes of the IAQ are now better (look in fact very good except for 10-30 MeV), but overall the IAQ’s are too low. Here a table with the underestimations:

Removing the V1 and V2+V3 veto dome transmission did improve things. Below the control plots that compare the GammaLib generated IAQs to the simulated UNH IAQs:

And here the comparison of the integrated probabilities:

Below the Crab spectrum obtained with the GammaLib IAQs (red) compared to the spectrum obtained with the UNH IAQs (green). The black points is the Crab spectrum published by Kuiper.

A powerlaw fit gives the following results:

Code merged into devel.