Feature #2186
Implement COMPTEL IAQ computation software
Status: | Closed | Start date: | 08/27/2017 | |
---|---|---|---|---|
Priority: | Normal | Due date: | ||
Assigned To: | Knödlseder Jürgen | % Done: | 100% | |
Category: | - | |||
Target version: | 1.5.0 | |||
Duration: |
Description
The software that computes Instrument Response Functions (IAQ datasets) for COMPTEL data should be implemented in GammaLib to assure that instrument response functions can be generate for arbitrary energy intervals. This is also needed to assure the legacy of COMPTEL data as the IAQs are not stored in the HEASRC archive.
The code can be inspired by the COMPASS coderespit
. This code reads the following input files which would need to be integrated into GammaLib:
SDA
file containing D1 detector response informationSDB
file containing D2 detector response information
Here a call tree of the respit
code:
respsi | +-- UITPBG | +-- CPRINT | +-- RESPSR (get task parameters) | +-- RESIDI (get telescope constants) | +-- SPCIAQ (calculate IAQ) | | | +-- UTPRIV (read integer value from TPT) | | | +-- UTPRXV (read character value from TPT) | | | +-- DEFRAN | | | | | +-- RESPD1 (get D1 parameters from SDA database) | | | | | | | +-- PSDAOR (open SDA dataset) | | | +-- PSDARR (read SDA dataset) | | | +-- PSDACR (close SDA dataset) | | | | | +-- RESPD2 (get D2 parameters from SDB database) | | | | | | | +-- PSDBOR (open SDB dataset) | | | +-- PSDBRR (read SDB dataset) | | | +-- PSDBCR (close SDB dataset) | | | +-- UTPRIV (read integer value from TPT) | | | +-- GENWEI | | | +-- IAQWEI | | | | | +-- UTPRLV (read logical value from TPT) | | | | | +-- RESPD1 (get D1 parameters from SDA database) | | | | | | | +-- PSDAOR (open SDA dataset) | | | +-- PSDARR (read SDA dataset) | | | +-- PSDACR (close SDA dataset) | | | | | +-- RESPD2 (get D2 parameters from SDB database) | | | | | | | +-- PSDBOR (open SDB dataset) | | | +-- PSDBRR (read SDB dataset) | | | +-- PSDBCR (close SDB dataset) | | | | | +-- INITIL (initialise D2 response) | | | | | +-- RESPSC (compute response for one Phibar layer) | | | | | | | +-- RESINT (response integration) | | | | | +-- RESPSG (get Klein Nishina value for Phigeo) | | | | | | | +-- RESSLV (substitutes values) | | | | | +-- RESWD1 (compute total D1 interaction coefficient) | | | | | +-- RESAD1 (compute the attenuation coefficient above D1) | | | | | +-- RESSVE (compute the selfveto attenuation coefficient) | | | | | +-- RESMHT (compute the multi-hit probability) | | | | | +-- RESAD2 (compute the attenuation coefficient between D1 and D2) | | | | | +-- RESWD2 (compute the D2 interaction coefficient) | | | | | +-- MULTIS (compute the multi-scatter transmission inside D1 module) | | | +-- ITOFCO (Generates FAQ from IAQ) | | | +-- POLINT (Polynomial interpolation) | | | +-- UTPRRV (read real value from TPT) | | | +-- PFAQOW (open FAQ FITS file) | +-- PFAQWR (write FAQ FITS file) | +-- P3GDSI (...) | +-- PFAQCW (close FAQ FITS file) | | | +-- UTPRLV (read logical value from TPT) | | | +-- FNDSCL (...) | | | +-- PIAQOP (open IAQ FITS file) | +-- PIAQWR (write IAQ FITS file) | +-- PIAQCW (close IAQ FITS file) | +-- UITPND
Recurrence
No recurrence.
History
#1 Updated by Knödlseder Jürgen over 7 years ago
- Assigned To set to Knödlseder Jürgen
- % Done changed from 0 to 10
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.
#2 Updated by Knödlseder Jürgen over 7 years ago
- % Done changed from 10 to 20
Added GCOMD1Response
and GCOMD2Response
classes to deal with SDA
and SDB
datasets.
#3 Updated by Knödlseder Jürgen over 7 years ago
Added GCOMInstChar
class to deal with the ICT
dataset.
Added GCOMIaq
class to implement the COMPASS task respsi
. So far the folding with the D1 and D2 response is working. Normalizing the D2 response (since only the photo-peak is implemented so far), the following probability sums (over Phibar) are obtained for each Phigeo value (should be 1):
=== GCOMIaq === Number of Phigeo bins .....: 55 Number of Phibar bins .....: 25 Maximum Phigeo value ......: 55 deg Maximum Phibar value ......: 50 deg Phigeo bin size ...........: 1 deg Phibar bin size ...........: 2 deg Energy range ..............: 1 MeV - 3 MeV phigeo=0.5 etrue1=0.000298014 etrue2=1.9997 Prob=0 phigeo=1.5 etrue1=0.0026788 etrue2=1.99732 Prob=0 phigeo=2.5 etrue1=0.00742268 etrue2=1.99258 Prob=0 phigeo=3.5 etrue1=0.0144946 etrue2=1.98551 Prob=0 phigeo=4.5 etrue1=0.0238428 etrue2=1.97616 Prob=0 phigeo=5.5 etrue1=0.0353998 etrue2=1.9646 Prob=0 phigeo=6.5 etrue1=0.0490834 etrue2=1.95092 Prob=0.0208068 phigeo=7.5 etrue1=0.0647983 etrue2=1.9352 Prob=0.330892 phigeo=8.5 etrue1=0.0824376 etrue2=1.91756 Prob=0.819893 phigeo=9.5 etrue1=0.101885 etrue2=1.89812 Prob=0.981478 phigeo=10.5 etrue1=0.123015 etrue2=1.87698 Prob=0.9973 phigeo=11.5 etrue1=0.145698 etrue2=1.8543 Prob=0.9973 phigeo=12.5 etrue1=0.169797 etrue2=1.8302 Prob=0.9973 phigeo=13.5 etrue1=0.195176 etrue2=1.80482 Prob=0.9973 phigeo=14.5 etrue1=0.221696 etrue2=1.7783 Prob=0.997287 phigeo=15.5 etrue1=0.249218 etrue2=1.75078 Prob=0.997241 phigeo=16.5 etrue1=0.277608 etrue2=1.72239 Prob=0.997169 phigeo=17.5 etrue1=0.306731 etrue2=1.69327 Prob=0.997077 phigeo=18.5 etrue1=0.336461 etrue2=1.66354 Prob=0.99697 phigeo=19.5 etrue1=0.366674 etrue2=1.63333 Prob=0.997022 phigeo=20.5 etrue1=0.397255 etrue2=1.60275 Prob=0.997077 phigeo=21.5 etrue1=0.428092 etrue2=1.57191 Prob=0.997123 phigeo=22.5 etrue1=0.459083 etrue2=1.54092 Prob=0.997159 phigeo=23.5 etrue1=0.490132 etrue2=1.50987 Prob=0.997213 phigeo=24.5 etrue1=0.521151 etrue2=1.47885 Prob=0.997263 phigeo=25.5 etrue1=0.552059 etrue2=1.44794 Prob=0.99729 phigeo=26.5 etrue1=0.582783 etrue2=1.41722 Prob=0.9973 phigeo=27.5 etrue1=0.613257 etrue2=1.38674 Prob=0.997294 phigeo=28.5 etrue1=0.643422 etrue2=1.35658 Prob=0.997275 phigeo=29.5 etrue1=0.673225 etrue2=1.32677 Prob=0.997286 phigeo=30.5 etrue1=0.702621 etrue2=1.29738 Prob=0.997292 phigeo=31.5 etrue1=0.73157 etrue2=1.26843 Prob=0.997296 phigeo=32.5 etrue1=0.760036 etrue2=1.23996 Prob=0.997298 phigeo=33.5 etrue1=0.787992 etrue2=1.21201 Prob=0.9973 phigeo=34.5 etrue1=0.815414 etrue2=1.18459 Prob=0.9973 phigeo=35.5 etrue1=0.84228 etrue2=1.15772 Prob=0.9973 phigeo=36.5 etrue1=0.868577 etrue2=1.13142 Prob=0.9973 phigeo=37.5 etrue1=0.894291 etrue2=1.10571 Prob=0.997298 phigeo=38.5 etrue1=0.919416 etrue2=1.08058 Prob=0.997296 phigeo=39.5 etrue1=0.943944 etrue2=1.05606 Prob=0.997295 phigeo=40.5 etrue1=0.967874 etrue2=1.03213 Prob=0.997292 phigeo=41.5 etrue1=0.991205 etrue2=1.0088 Prob=0.99729 phigeo=42.5 etrue1=1.01394 etrue2=0.986062 Prob=0.997286 phigeo=43.5 etrue1=1.03608 etrue2=0.963921 Prob=0.997266 phigeo=44.5 etrue1=1.05763 etrue2=0.94237 Prob=0.997032 phigeo=45.5 etrue1=1.0786 etrue2=0.921399 Prob=0.994997 phigeo=46.5 etrue1=1.099 etrue2=0.901002 Prob=0.98327 phigeo=47.5 etrue1=1.11883 etrue2=0.881171 Prob=0.938375 phigeo=48.5 etrue1=1.13811 etrue2=0.861894 Prob=0.823054 phigeo=49.5 etrue1=1.15684 etrue2=0.843161 Prob=0.622092 phigeo=50.5 etrue1=1.17504 etrue2=0.824962 Prob=0.381694 phigeo=51.5 etrue1=1.19272 etrue2=0.807285 Prob=0.181873 phigeo=52.5 etrue1=1.20988 etrue2=0.790117 Prob=0.0655009 phigeo=53.5 etrue1=1.22655 etrue2=0.773448 Prob=0.0175808 phigeo=54.5 etrue1=1.24274 etrue2=0.757264 Prob=0.00349692
The difference from one for small phigeo values is due to the energy threshold, that at large phigeo values is due to the maximum phibar value.
#4 Updated by Knödlseder Jürgen over 7 years ago
- % Done changed from 20 to 30
Note that geometrical smearing is not implemented since this neads a NAG routine. I’m not sure that this is very standard.
#5 Updated by Knödlseder Jürgen over 7 years ago
- % Done changed from 30 to 50
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::GCOMIaq
: test input argument validityGCOMIaq::set(const GEnergy& energy)
: implement geometrical smearingGCOMIaq::set(const GModelSpectral& spectrum)
: implement geometrical smearingGCOMIaq::compute_iaq_bin
: study impact of integration precision (needs full D2 response)GCOMIaq::klein_nishina_integral
: check this formula
GCOMInstChars
:
GCOMInstChars::prob_D2inter
: verify formulaGCOMInstChars::prob_multihit
: Implement method. Need to find out what exactly is computed.GCOMInstChars::atten_D1
: Implement methodGCOMInstChars::atten_D2
: Implement methodGCOMInstChars::atten_selfveto
: Implement methodGCOMInstChars::multi_scatter
: Implement method
GCOMD2Response
:
GCOMD2Response::operator()
: Add non photo peak response componentsGCOMD2Response::update_cache
: Assure normalization
GCOMResponse
:
- Use
GCOMIaq
class
#6 Updated by Knödlseder Jürgen over 7 years ago
- % Done changed from 50 to 60
Advancing. Here the new list:
GCOMIaq
:
GCOMIaq::GCOMIaq
: test input argument validityGCOMIaq::set(const GEnergy& energy)
: implement geometrical smearingGCOMIaq::set(const GModelSpectral& spectrum)
: implement geometrical smearingGCOMIaq::compute_iaq_bin
: study impact of integration precision (needs full D2 response)GCOMIaq::klein_nishina_integral
: check this formula
GCOMInstChars
:
GCOMInstChars::multi_scatter
: Implement method
GCOMD2Response
:
GCOMD2Response::operator()
: Add Compton response componentsGCOMD2Response::update_cache
: Assure normalization
GCOMResponse
:
- Use
GCOMIaq
class
#7 Updated by Knödlseder Jürgen over 7 years ago
C* There is a call to RESD2 which has the following subroutines: * C* 1.COMPTT______________________ * C* | | | * C* RSD2NO RESFUN RESD2W * C* | | | * C* KLNFCT KLNFCT D01ALF * C* | * C* RESFUN * C* | * C* KLNFCT * C* 2.GAUSS *
#8 Updated by Knödlseder Jürgen over 7 years ago
- File D2_response.png added
- Status changed from New to In Progress
- Target version set to 1.5.0
- % Done changed from 60 to 70
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.
#9 Updated by Knödlseder Jürgen over 7 years ago
- File ratio-075-1MeV-10bins.png added
- File ratio-1-3MeV-10bins.png added
- File ratio-3-10MeV-10bins.png added
- File ratio-10-30MeV-10bins.png added
- File ratio-1-3MeV-20bins.png added
- File ratio-075-1MeV-20bins.png added
- File ratio-3-10MeV-20bins.png added
- File ratio-10-30MeV-20bins.png added
- File ratio-075-1MeV-40bins.png added
- File ratio-1-3MeV-40bins.png added
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.
#10 Updated by Knödlseder Jürgen over 7 years ago
- File ratio-3-10MeV-40bins.png added
- File ratio-10-30MeV-40bins.png added
- File ratio-075-1MeV-80bins.png added
- File ratio-1-3MeV-80bins.png added
- File ratio-3-10MeV-80bins.png added
- File ratio-10-30MeV-80bins.png added
#11 Updated by Knödlseder Jürgen over 7 years ago
Here the new list as it is now:
GCOMIaq
:
GCOMIaq::GCOMIaq
: test input argument validityGCOMIaq::set(const GEnergy& energy)
: implement geometrical smearingGCOMIaq::set(const GModelSpectral& spectrum)
: implement geometrical smearingGCOMIaq::compute_iaq_bin
: study impact of integration precision (needs full D2 response)GCOMIaq::GCOMIaq
: move energy boundaries toset()
method
GCOMInstChars
:
GCOMInstChars::multi_scatter
: Implement method
GCOMResponse
:
- Use
GCOMIaq
class
#12 Updated by Knödlseder Jürgen over 7 years ago
- File ratio-1-3MeV-1e-3.png added
- File ratio-1-3MeV-1e-4.png added
- File ratio-1-3MeV-1e-5.png added
- File ratio-1-3MeV-1e-6.png added
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.
Precision | Integral response | Computing time |
1e-3 | 0.04911986 | 27.124 s |
1e-4 | 0.04911092 | 25.411 s |
1e-5 | 0.04910129 | 29.122 s |
1e-6 | 0.0490938 | 48.814 s |
I therefore keep the precision at 1e-4.
#13 Updated by Knödlseder Jürgen over 7 years ago
- File comparison-1-3MeV.png added
And here a nice comparison summary plot, generated using iaqcmp.py
:
#14 Updated by Knödlseder Jürgen over 7 years ago
- File zenith-0.png added
- File zenith-20.png added
- File zenith-30.png added
- File zenith-40.png added
- % Done changed from 70 to 80
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.
#15 Updated by Knödlseder Jürgen over 7 years ago
- File final-075-1MeV.png added
- File final-1-3MeV.png added
- File final-3-10MeV.png added
- File final-10-30MeV.png added
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.
Energy range | UNH | GammaLib | Difference |
0.75 - 1 MeV | 0.024072 | 0.016430 | -31.7% |
1 - 3 MeV | 0.044088 | 0.038275 | -13.2% |
3 - 10 MeV | 0.033936 | 0.037473 | +10.4% |
10 - 30 MeV | 0.017635 | 0.018845 | +6.9% |
Here the list of remaining things to be done:
GCOMIaq
:
GCOMIaq::GCOMIaq
: test input argument validity
GCOMResponse
:
- Use GCOMIaq class
#16 Updated by Knödlseder Jürgen over 7 years ago
- File crab_iaq.png added
- File howto_comptel_spectrum.png added
- % Done changed from 80 to 90
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:
Parameter | UNH | GammaLib |
Optimized function value | -153218.344 | -153207.016 |
Test Statistic | 915.651 | 892.995 |
Prefactor | 0.002375 +/- 0.000107 | 0.002349 +/- 0.000115 |
Index | 2.472 +/- 0.050 | 2.654 +/- 0.055 |
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.
#17 Updated by Knödlseder Jürgen over 7 years ago
- File multi_scatter.png added
- File ne231a_mfpath.png added
- File kn_cross_section.png added
I cross-checked extensively the multi-scattering code and it seems to be correct. Here the essential functions that are used:
#18 Updated by Knödlseder Jürgen over 7 years ago
- File prob_D1inter.png added
- File prob_D2inter.png added
- File trans_D1.png added
- File trans_D2.png added
- File trans_V1.png added
- File transV23.png added
- File psd_correction.png added
- File psd_correction_vs_E1.png added
- File prob_no_multihit.png added
- File prob_no_selfveto.png added
Here all the other probabilities that come into play in the computation of the IAQs:
#19 Updated by Knödlseder Jürgen over 7 years ago
- File prob_D1inter-corrected.png added
- File prob_D2inter-corrected.png added
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:
Energies | Underestimation |
0.75-1 MeV | -35.6% |
1-3 MeV | -24.1% |
3-10 MeV | -14.3% |
10-30 MeV | -18.9% |
#20 Updated by Knödlseder Jürgen over 7 years ago
- File check-0.75-1MeV.png added
- File check-1-3MeV.png added
- File check-3-10MeV.png added
- File check-10-30MeV.png added
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:
Energy range | UNH | GammaLib | GammaLib-UNH |
0.75-1 MeV | 0.024072 | 0.024069 | -0.0% |
1-3 MeV | 0.044088 | 0.048042 | 9.0% |
3-10 MeV | 0.033936 | 0.036550 | 7.7% |
10-30 MeV | 0.017635 | 0.016684 | -5.4% |
#21 Updated by Knödlseder Jürgen over 7 years ago
- File crab-spectrum.png added
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:
Parameter | UNH | GammaLib |
Optimized function value | -153218.344 | -153250.177 |
Test Statistic | 915.651 | 979.317 |
Prefactor | 0.002375 +/- 0.000107 | 0.001859 +/- 0.000084 |
Index | 2.472 +/- 0.050 | 2.393 +/- 0.048 |
#22 Updated by Knödlseder Jürgen over 7 years ago
- Status changed from In Progress to Closed
- % Done changed from 90 to 100
Code merged into devel
.