Feature #2187

Implement COMPTEL DRI computation software

Added by Knödlseder Jürgen about 7 years ago. Updated about 7 years ago.

Status:ClosedStart date:08/27/2017
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.5.0
Duration:

Description

A software should be implement for COMPTEL that converts an event list into a DRI dataset. In COMPASS this job is done using the skydri software.

Here the layout of skydrs12:

skydir
|
+-- UITPBG (initialise user shell)
|
+-- SETSW (set global switches)
|
+-- RDTPT (read TPT)
|
+-- AUTOCD (generate automatically coordinates)
|
+-- PNDRD (print DRI definitions)
|
+-- SETX (generate DRX dataset)
|   |
|   +-- ZERO (zero array)
|   |
|   +-- SETXT (generate DRX for one response group, i.e. energy)
|   |   |
|   |   +-- GALCEL (coordinate conversion)
|   |   |
|   |   +-- P7OPEN (open TIM dataset)
|   |   +-- POADOS (open OAD dataset)
|   |   +-- POADRS (read OAD data)
|   |   +-- P7TIME (check against time in TIM dataset)
|   |   |
|   |   +-- COORD1 (coordinate conversion)
|   |   |
|   |   +-- EXPOLM (compute exposure array)
|   |   |   |
|   |   |   +-- COORD1 (coordinate conversion)
|   |   |   |
|   |   |   +-- SENSAR (compute exposure summed over all D1)
|   |   |
|   |   +-- P7CLOS (close TIM dataset)
|   |
|   +-- DRXOUT (output DRX array to FITS file)
|
+-- SETE (generate DRE dataset)
|   |
|   +-- EVSTUP (initialise event statistics)
|   |
|   +-- ZERO (zero array)
|   |
|   +-- SETET (generate DRE dataset for one response group, i.e. energy)
|   |   |
|   |   +-- GALCEL (coordinate conversion)
|   |   |
|   |   +-- P7OPEN (open TIM dataset)
|   |   +-- POADOS (open OAD dataset)
|   |   +-- PEVPOR (open EVP dataset)
|   |   +-- POADRS (read OAD data)
|   |   +-- P7TIME (check against time in TIM dataset)
|   |   |
|   |   +-- COORD1 (coordinate conversion)
|   |   |
|   |   +-- GEOANG (compute angle subtended by Earth)
|   |   |
|   |   +-- MODSTA (update module status table)
|   |   |   |
|   |   |   +-- PMODST (??? read module status from FPM ???)
|   |   |
|   |   +-- EVENLM (bin events)
|   |   |   |
|   |   |   +-- COORD1 (coordinate conversion)
|   |   |   |
|   |   |   +-- PEVPSR (read next event)
|   |   |   |
|   |   |   +-- PSSEVP (apply event selection) ??? where is this code ???
|   |   |   |
|   |   |   +-- COORD1 (coordinate conversion)
|   |   |   |
|   |   |   +-- EVSTUP (update event statistics)
|   |   |
|   |   +-- P7CLOS (close TIM dataset)
|   |   +-- PEVPCR (close EVP dataset)
|   |
|   +-- DREOUT (output DRE array to FITS file)
|
+-- SETG (generate DRE dataset)
|
+-- UITPND (terminate user shell)

dre-difference.png (53 KB) Knödlseder Jürgen, 09/06/2017 10:18 AM

dre-difference-plot.png (106 KB) Knödlseder Jürgen, 09/06/2017 11:17 AM

dre-cmp-noeha.png (108 KB) Knödlseder Jürgen, 09/06/2017 11:19 AM

vp1-crab-comparison.png (73 KB) Knödlseder Jürgen, 09/08/2017 01:27 PM

dre-drb-vp5-1.8MeV.png (101 KB) Knödlseder Jürgen, 09/08/2017 04:57 PM

crab-nodes-vp1.png (23.9 KB) Knödlseder Jürgen, 09/08/2017 05:02 PM

Dre-difference Dre-difference-plot Dre-cmp-noeha Vp1-crab-comparison Dre-drb-vp5-1.8mev Crab-nodes-vp1

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen about 7 years ago

  • Status changed from New to In Progress
  • Assigned To set to Knödlseder Jürgen
  • Target version set to 1.5.0
  • % Done changed from 0 to 50

The issue makes good progress. I have added GCOMEventAtom, @GCOMEventList and GCOMRoi classes to have EVP files, and I also added the GCOMTim class to handle TIM files. Now need the OAD file classes and things are in place. Since an EVP and TIM files typically covers 14 days, but OAD files cover only one day I decided to implement the following XML file structure for an unbinned analysis:

<?xml version="1.0" standalone="no"?>
<observation_list title="observation library">
  <observation name="Crab" id="000001" instrument="COM">
    <parameter name="EVP" file="m16992_evp.fits"/>
    <parameter name="TIM" file="m10695_tim.fits"/>
    <parameter name="OAD" file="m20039_oad.fits"/>
    <parameter name="OAD" file="m20041_oad.fits"/>
  </observation>
</observation_list>

This allows to append an arbitrary list of OAD files to a given observation.

#2 Updated by Knödlseder Jürgen about 7 years ago

I made good progress, a first event binning method is now in GCOMObservation. Running on the entire VP1, using the following XML file

<?xml version="1.0" standalone="no"?>
<observation_list title="observation library">
  <observation name="Crab" id="000001" instrument="COM">
    <parameter name="EVP" file="/project-data/comptel/old.2/phase01/vp0001_0/m16992_evp.fits"/>
    <parameter name="TIM" file="/project-data/comptel/old.2/phase01/vp0001_0/m10695_tim.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20037_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20039_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20041_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20045_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20048_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20050_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20054_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20058_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20064_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20066_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20068_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20071_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20073_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20078_oad.fits"/>
    <parameter name="OAD" file="/project-data/comptel/old.2/phase01/vp0001_0/m20081_oad.fits"/>
  </observation>
</observation_list>

I obtain the following result for the 1-3 MeV DRE:
GCOMObservation::compute_dre
============================
Used superpackets ............: 52873
Skipped superpackets .........: 26210
Used events ..................: 483723
Events outside superpacket ...: 135838
Energy too low ...............: 78054
Energy too high ..............: 294900
Phibar too low ...............: 0
Phibar too high ..............: 0
Outside D1 selection .........: 0
Outside D2 selection .........: 0
Outside TOF selection ........: 0
Outside PSD selection ........: 0
Outside rejection flag sel. ..: 0
Outside DRE cube .............: 0
No scatter angle .............: 0
Bad module combination .......: 0
Earth horizon angle too small : 31853

with the following event cube
=== GCOMEventCube ===
 Number of events ..........: 483723
 Number of elements ........: 140600
 Size (Chi x Psi x Phi) ....: 76 x 74 x 25
 Energy range ..............: 1 - 3 MeV
 Mean energy ...............: 1.73205080756888 MeV
 Energy bin width ..........: 2 MeV
 MJD interval ..............: 48392.9035427778 - 48406.7861383319 days
 Mean time .................: -587274271.960063 s (TT)
 Ontime ....................: 1199456.25587511 s

For reference, the file m50439_dre.fits gives the following result:
=== GCOMEventCube ===
 Number of events ..........: 316141
 Number of elements ........: 140600
 Size (Chi x Psi x Phi) ....: 76 x 74 x 25
 Energy range ..............: 1 - 3 MeV
 Mean energy ...............: 1.73205080756888 MeV
 Energy bin width ..........: 2 MeV
 MJD interval ..............: 48392.9035427778 - 48406.7861383319 days
 Mean time .................: -587274271.960063 s (TT)
 Ontime ....................: 1199456.25587511 s

Maybe the difference comes from the module status. Here the number of events per module after accumulating 20000 event for each time step:
([3418, 3240, 3106, 3502, 3830, 2905, 0], [0, 0, 1761, 1872, 0, 1925, 1880, 1755, 1866, 1855, 1912, 1697, 1641, 1837])
([3301, 3325, 2996, 3523, 3879, 2977, 0], [0, 0, 1799, 1833, 0, 1882, 1946, 1793, 1719, 1791, 1990, 1808, 1662, 1778])
([3417, 3378, 2917, 3435, 3881, 2973, 0], [0, 0, 1850, 1911, 0, 1797, 1820, 1853, 1799, 1711, 1995, 1813, 1671, 1781])
([3372, 3157, 2942, 3554, 4024, 2952, 0], [0, 0, 1902, 1819, 0, 1879, 1943, 1801, 1683, 1783, 2047, 1694, 1646, 1804])
([3356, 3147, 3000, 3556, 4027, 2915, 0], [0, 0, 1890, 1834, 0, 1858, 1972, 1726, 1815, 1837, 1928, 1694, 1623, 1824])
([3375, 3422, 3142, 3423, 3840, 2799, 0], [0, 0, 1852, 1899, 0, 1821, 1898, 1834, 1701, 1810, 1900, 1827, 1699, 1760])
([3386, 3389, 2966, 3440, 3881, 2939, 0], [0, 0, 1860, 1913, 0, 1801, 1934, 1862, 1728, 1737, 2007, 1783, 1639, 1737])
([3310, 3059, 2978, 3656, 4073, 2925, 0], [0, 0, 1980, 1762, 0, 1823, 1992, 1710, 1782, 1772, 2040, 1729, 1582, 1829])
([3508, 3258, 2968, 3449, 3788, 3030, 0], [0, 0, 1796, 1874, 0, 1807, 1777, 1840, 1809, 1891, 1932, 1756, 1703, 1816])
([3369, 3457, 3068, 3395, 3779, 2933, 0], [0, 0, 1765, 1913, 0, 1858, 1916, 1861, 1846, 1732, 1927, 1737, 1681, 1765])
([3368, 3448, 2888, 3470, 3909, 2918, 0], [0, 0, 1815, 1882, 0, 1785, 1831, 1927, 1773, 1857, 1915, 1809, 1717, 1690])
([3326, 3123, 2895, 3613, 4052, 2992, 0], [0, 0, 2011, 1755, 0, 1921, 2095, 1707, 1664, 1768, 2145, 1586, 1545, 1804])
([3416, 3315, 3021, 3465, 3899, 2885, 0], [0, 0, 1849, 1864, 0, 1909, 1869, 1808, 1788, 1786, 1974, 1725, 1601, 1828])
([3491, 3464, 3095, 3256, 3755, 2940, 0], [0, 0, 1766, 1970, 0, 1772, 1778, 2037, 1801, 1741, 1916, 1749, 1678, 1793])
([3529, 3327, 2957, 3358, 3897, 2933, 0], [0, 0, 1936, 1815, 0, 1836, 1966, 1748, 1759, 1802, 2019, 1761, 1605, 1754])
([3413, 3082, 2828, 3717, 4016, 2945, 0], [0, 0, 1977, 1762, 0, 1915, 2149, 1679, 1654, 1798, 2119, 1621, 1596, 1731])
([3452, 3303, 2988, 3563, 3822, 2873, 0], [0, 0, 1807, 1951, 0, 1908, 1891, 1793, 1756, 1816, 1927, 1768, 1587, 1797])
([3504, 3377, 3079, 3377, 3753, 2911, 0], [0, 0, 1936, 1907, 0, 1798, 1856, 1836, 1797, 1674, 1946, 1816, 1663, 1772])
([3458, 3480, 2942, 3419, 3856, 2846, 0], [0, 0, 1861, 1898, 0, 1776, 1914, 1909, 1794, 1673, 1937, 1769, 1654, 1816])
([3338, 3084, 2899, 3654, 4116, 2910, 0], [0, 0, 2044, 1800, 0, 1869, 2144, 1535, 1664, 1822, 2116, 1579, 1612, 1816])
([3428, 3315, 3033, 3506, 3803, 2916, 0], [0, 0, 1771, 1881, 0, 1840, 1951, 1824, 1807, 1818, 1987, 1734, 1642, 1746])
([3476, 3461, 3032, 3414, 3678, 2940, 0], [0, 0, 1803, 1890, 0, 1819, 1872, 1945, 1762, 1692, 1941, 1829, 1683, 1765])
([3461, 3331, 2959, 3401, 3933, 2916, 0], [0, 0, 1927, 1906, 0, 1812, 1940, 1776, 1735, 1710, 1936, 1814, 1637, 1808])
([3355, 3021, 2913, 3635, 4218, 2859, 0], [0, 0, 1988, 1801, 0, 1916, 2050, 1610, 1744, 1838, 2192, 1592, 1622, 1648])
([3434, 3297, 2926, 3363, 3972, 3009, 0], [676, 0, 1917, 2069, 0, 1941, 2011, 1876, 1850, 1913, 2062, 1921, 1765, 0])
([3423, 3315, 3029, 3422, 3877, 2935, 0], [1818, 0, 1758, 1936, 0, 1802, 1886, 1923, 1780, 1715, 1936, 1807, 1640, 0])
([3287, 3179, 2747, 3605, 4163, 3020, 0], [1781, 0, 2062, 1812, 0, 1771, 2039, 1744, 1644, 1732, 2107, 1729, 1580, 0])
([3362, 3021, 2878, 3586, 4053, 3101, 0], [1833, 0, 1929, 1738, 0, 1971, 2014, 1643, 1740, 1811, 2081, 1640, 1601, 0])
([3404, 3394, 2945, 3533, 3672, 3053, 0], [1883, 0, 1793, 1872, 0, 1720, 1831, 1927, 1786, 1801, 1814, 1776, 1798, 0])
([3415, 3330, 2792, 3499, 3864, 3101, 0], [1787, 0, 1956, 1893, 0, 1873, 1940, 1783, 1765, 1757, 1871, 1751, 1625, 0])
([3422, 3061, 2719, 3602, 4225, 2972, 0], [1774, 0, 1959, 1743, 0, 1902, 2117, 1710, 1727, 1725, 2097, 1669, 1578, 0])
([3436, 3193, 2864, 3656, 3910, 2942, 0], [1738, 0, 1924, 1799, 0, 1866, 1943, 1717, 1768, 1887, 2024, 1649, 1686, 0])
([3453, 3421, 3000, 3419, 3776, 2932, 0], [1841, 0, 1776, 1946, 0, 1745, 1802, 1905, 1774, 1730, 2028, 1793, 1661, 0])
([3445, 3422, 2946, 3399, 3852, 2937, 0], [1817, 0, 1777, 1856, 0, 1822, 1870, 1895, 1771, 1777, 1879, 1781, 1756, 0])
([3395, 3007, 2868, 3591, 4130, 3010, 0], [1781, 0, 1997, 1717, 0, 1863, 2053, 1669, 1685, 1844, 2155, 1629, 1608, 0])
([3433, 3161, 2907, 3549, 3925, 3026, 0], [1804, 0, 1788, 1911, 0, 1829, 1786, 1816, 1827, 1962, 1994, 1668, 1616, 0])
([3376, 3282, 2952, 3463, 3981, 2947, 0], [1809, 0, 1834, 1763, 0, 1811, 1841, 1914, 1747, 1759, 2032, 1792, 1699, 0])
([3375, 3358, 2899, 3437, 3979, 2953, 0], [1737, 0, 1824, 1861, 0, 1802, 1855, 1829, 1785, 1842, 2039, 1799, 1628, 0])
([3503, 3401, 2851, 3433, 3927, 2886, 0], [1802, 0, 1728, 1871, 0, 1835, 1768, 1948, 1806, 1809, 1971, 1748, 1715, 0])
([3406, 3287, 2967, 3398, 3951, 2992, 0], [1810, 0, 1838, 1897, 0, 1829, 1806, 1770, 1767, 1747, 1985, 1812, 1740, 0])
([3457, 3112, 2814, 3519, 4144, 2955, 0], [1801, 0, 1889, 1705, 0, 1845, 2033, 1754, 1687, 1755, 2124, 1783, 1625, 0])
([3422, 3091, 2903, 3570, 3992, 3023, 0], [1771, 0, 1836, 1768, 0, 1940, 1907, 1727, 1779, 1802, 2130, 1690, 1651, 0])
([3490, 3364, 2904, 3318, 4000, 2925, 0], [1785, 0, 1825, 1855, 0, 1840, 1906, 1806, 1764, 1750, 1981, 1818, 1671, 0])
([3505, 3382, 2953, 3440, 3837, 2884, 0], [1915, 0, 1746, 1832, 0, 1730, 1835, 1860, 1751, 1795, 1954, 1860, 1723, 0])
([3367, 3020, 2719, 3663, 4248, 2984, 0], [1789, 0, 2046, 1655, 0, 1934, 2083, 1596, 1620, 1889, 2175, 1582, 1632, 0])
([3415, 3258, 2963, 3481, 3930, 2954, 0], [1807, 0, 1795, 1884, 0, 1816, 1847, 1890, 1723, 1838, 1916, 1802, 1683, 0])
([3390, 3350, 2868, 3523, 3950, 2920, 0], [1786, 0, 1893, 1879, 0, 1790, 1808, 1859, 1763, 1806, 1906, 1824, 1687, 0])
([3344, 3262, 2921, 3514, 4011, 2949, 0], [1845, 0, 1986, 1718, 0, 1834, 1922, 1882, 1704, 1770, 2068, 1674, 1598, 0])
([3384, 3020, 2910, 3650, 4107, 2930, 0], [1833, 0, 1266, 1819, 0, 1991, 2134, 1789, 1723, 1910, 2184, 1691, 1661, 0])
([3324, 3466, 2959, 3385, 4011, 2856, 0], [2018, 0, 0, 2010, 0, 1983, 2052, 2034, 1969, 2005, 2060, 1979, 1891, 0])
([3105, 3106, 2736, 3272, 3860, 2708, 1214], [2006, 0, 0, 2069, 0, 1985, 2077, 1991, 1952, 1975, 2231, 1904, 1811, 0])
([2924, 2818, 2487, 3107, 3450, 2515, 2700], [1934, 0, 0, 2023, 0, 2104, 2126, 1954, 1908, 2011, 2196, 1907, 1838, 0])

For D1, module 7 was dead most of the time. For D2 there is quite some intermittence in modules 1, 2, 3, 5 and 14. Excluding all these modules gives 396509 events, closer to the 316141 in the reference DRE dataset, but still too large.

#3 Updated by Knödlseder Jürgen about 7 years ago

  • % Done changed from 50 to 60

I just received the D1 & D2 module status table:

                // 8392,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8393,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8394,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8395,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8396,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8397,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,1
                // 8398,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8399,1,1,1,1,1,1,0, 0,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8400,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8401,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8402,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8403,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8404,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8405,1,1,1,1,1,1,0, 1,0,1,1,0,1,1,1,1,1,1,1,1,0
                // 8406,1,1,1,1,1,1,1, 1,0,1,1,0,1,1,1,1,1,1,1,1,0

Implementing this table hard-wired for the moment leads to a rejection of 3734 events, this is a minor correction.

I forgot however that the TOF selection is more stringent when compared to the EVP data. The data comprise the window 110 - 150, the standard selection is 115 - 130. Here the summary of the standard selection criteria that are now hard-wired:
Parameter Selection Interval
E1 70 - 20000 keV
E2 650 - 30000 keV
Phibar 0 - 50 deg
TOF 115 - 130
PSD 0 - 110
Zeta > 5 deg

And here the number of events in the standard VP1 DRE, compared to the numbers I get:

Energy Standard GCOMObservation::compute_dre Difference
0.75 - 1 MeV 48806 40788 -16.4%
1 - 3 MeV 316141 290487 -8.1%
3 - 10 MeV 140395 123601 -12.0%
10 - 30 MeV 22253 21473 -3.5 %

Interestingly, there is an energy dependent difference.

#4 Updated by Knödlseder Jürgen about 7 years ago

Interestingly, the difference is statistical in the sense that for some pixels there are more events in the standard DRE, is other there are more events in the GammaLib DRE. See image below as an example.

#5 Updated by Knödlseder Jürgen about 7 years ago

Making a more detailed comparison it look like if the difference is Phibar dependent:

#6 Updated by Knödlseder Jürgen about 7 years ago

That’s what happens without Earth horizon cut. Now the difference also exists at high Phibar values, hence this does not seem to be the problem.

#7 Updated by Knödlseder Jürgen about 7 years ago

Got it. The event data shipped with HEASARC are apparently not compatible with the DRE datasets.

I now compute also the DRG and DRX in GammaLib and there things are consistent. For illustration, the picture below shows a TS map for DRE, DRG and DRX dataset produced by GammaLib (left) compared to a TS map for the datasets provided by HEASARC. The image is for viewing period 1, the source that is visible is the Crab. The energy band 1-3 MeV was used for the image.

#8 Updated by Knödlseder Jürgen about 7 years ago

Here the comparison of a DRE and a DRG (used as DRB) for the viewing period 5 and the 1.8 MeV line band. There was an issue with a sign reversion in the Earth horizon computation, but now things seem to be okay.

#9 Updated by Knödlseder Jürgen about 7 years ago

Here the Crab standard band spectrum based on the GammaLib DRE, DRG and DRX files. Red are the GammaLib points, green are the points obtained with the HEASARC event cubes, and black is the spectrum from Kuiper. All is for viewing period 1. Things look again reasonable.

#10 Updated by Knödlseder Jürgen about 7 years ago

  • Status changed from Pull request to Closed

Code merged into devel.

Also available in: Atom PDF