gammalib / gammalib

Commit 24b8267034dc450975e00b217f89a114964f21e6

Authored by Jürgen Knödlseder
1 parent 38765419

Prepare for code integration

Inline Side-by-side
Showing 6 changed files with 103 additions and 20 deletions
ChangeLog
  View file @24b8267
1   -2023-02-10
  1 +2023-02-16
2 2  
3 3 * Version 2.1.0 released
4 4 ========================
5 5  
  6 + Add COMPTEL DRW support to GCOMObservation and GCOMDris (#4209)
6 7 Add GSkyMap comparison operators (#4209)
7 8 Add GIntegral::adaptive_gauss_kronrod() method (#4204)
8 9 Add optional unit parameter to GApplication::log_value() methods (#4202)
... ...
1 1 New Features and Important Changes in GammaLib 2.1.0
2 2  
3   -10 February 2023
  3 +16 February 2023
4 4  
5 5  
6 6 1. Introduction
... ... @@ -12,7 +12,7 @@ made since the last release of GammaLib.
12 12 2. Public interface modifications
13 13 ---------------------------------
14 14 The following classes have been added:
15   -- none
  15 +- GCOMDris
16 16  
17 17 The following classes have been removed:
18 18 - none
... ... @@ -31,7 +31,10 @@ The following methods have been added:
31 31 - GResponseVectorCache::load()
32 32 - GResponseVectorCache::save()
33 33 - GResponseVectorCache::read()
  34 +- GCOMObservation(GCOMDri&, GCOMDri&, GCOMDri&, GCOMDri&, GCOMDri&)
34 35 - GCOMObservation::npred()
  36 +- GCOMObservation::drw()
  37 +- GCOMObservation::drwname()
35 38 - GCOMObservation::cachename()
36 39 - GCOMResponse::load_cache()
37 40 - GCOMResponse::save_cache()
... ... @@ -44,6 +47,10 @@ The following methods have been renamed:
44 47  
45 48 The arguments for the following methods have been changed:
46 49 - GApplication::log_value() methods now have optional unit parameters
  50 +- GCOMObservation() filename constructor now takes 5 filenames on
  51 + input; was four filenames before
  52 +- GCOMObservation::load() method now takes 5 filenames on input; was
  53 + four filenames before
47 54  
48 55 The return value of the following methods has been changed:
49 56 - none
... ... @@ -162,6 +169,16 @@ None
162 169  
163 170 22. COMPTEL interface
164 171 ---------------------
  172 +Added GCOMDris container class. This class specifically provides the
  173 +GCOMDris::compute_drws() that allows computation of weighting cubes for
  174 +all energy bins (#4209).
  175 +
  176 +Added support for weighting cubes in GCOMObservation class. The filename
  177 +constructor takes five filenames, with the third filename being in the DRW
  178 +filename. A GCOMObservation::drw() method was added to retrieve the DRW,
  179 +GCOMObservation::drwname() methods were added to set and retrieve the DRW
  180 +filename (#4209).
  181 +
165 182 Added GCOMObservation::cachename() method to set and retrieve the name of a
166 183 response cache file. If set, the response cache file will automatically
167 184 written or loaded when handling XML files. In addition, GCOMResponse methods
... ...
README.md
  View file @24b8267
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.1.0.dev (10 February 2023)
  3 +* Version: 2.1.0.dev (16 February 2023)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/2.1.rst
  View file @24b8267
... ... @@ -10,7 +10,7 @@ GammaLib 2.1 is a major release that adds significant functionality.
10 10  
11 11 In particular, this release provides:
12 12  
13   -* ...
  13 +* Weighting cube DRW support for COMPTEL interface
14 14  
15 15  
16 16 Bug fixes
... ... @@ -24,7 +24,8 @@ Improvements
24 24 ------------
25 25  
26 26 * [`4209 <https://cta-redmine.irap.omp.eu/issues/4209>`_] -
27   - Add ``GSkyMap`` comparison operators (#4209)
  27 + Add ``GSkyMap`` comparison operators.
  28 + Add COMPTEL DRW support to ``GCOMObservation`` and add ``GCOMDris`` class.
28 29 * [`4204 <https://cta-redmine.irap.omp.eu/issues/4204>`_] -
29 30 Add ``GIntegral::adaptive_gauss_kronrod()`` method
30 31 * [`4202 <https://cta-redmine.irap.omp.eu/issues/4202>`_] -
... ...
inst/com/include/GCOMDris.hpp
  View file @24b8267
... ... @@ -90,14 +90,14 @@ protected:
90 90 void free_members(void);
91 91 void compute_drws_energy(const GCOMObservation& obs,
92 92 const GCOMEventList* events,
93   - const GCOMSelection& select = GCOMSelection(),
94   - const double& zetamin = 5.0,
95   - const double& timebin = 300.0);
  93 + const GCOMSelection& select,
  94 + const double& zetamin,
  95 + const double& timebin);
96 96 void compute_drws_phibar(const GCOMObservation& obs,
97 97 const GCOMEventList* events,
98   - const GCOMSelection& select = GCOMSelection(),
99   - const double& zetamin = 5.0,
100   - const double& timebin = 300.0);
  98 + const GCOMSelection& select,
  99 + const double& zetamin,
  100 + const double& timebin);
101 101  
102 102 // Protected data members
103 103 std::vector<GCOMDri> m_dris; //!< Data space instances
... ...
inst/com/src/GCOMDris.cpp
  View file @24b8267
... ... @@ -28,6 +28,8 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
  31 +#include <fstream>
  32 +#include <iostream>
31 33 #include "GException.hpp"
32 34 #include "GMath.hpp"
33 35 #include "GCOMTools.hpp"
... ... @@ -53,12 +55,12 @@
53 55 /* __ Macros _____________________________________________________________ */
54 56  
55 57 /* __ Coding definitions _________________________________________________ */
  58 +#define G_COMPUTE_DRWS_EHACORR //!< Correct Earth horizon cut
56 59  
57 60 /* __ Debug definitions __________________________________________________ */
58   -#define G_DEBUG_COMPUTE_DRWS
  61 +//#define G_DEBUG_COMPUTE_DRWS
59 62  
60 63 /* __ Constants __________________________________________________________ */
61   -const double superpacket_duration = 16.384; // Duration of superpacket (s)
62 64  
63 65  
64 66 /*==========================================================================
... ... @@ -876,7 +878,7 @@ void GCOMDris::compute_drws_energy(const GCOMObservation&amp; obs,
876 878 * weighted geometry cubes which were multiplied by the solid angle of
877 879 * the (Chi,Psi) pixels.
878 880 *
879   - * For this method the event rates depend on time and Phibar angle. For
  881 + * For this method the event rates depend on time and Phibar angle. For
880 882 * each superpacket, the event rates are linearly interpolated to the
881 883 * relevant superpacket time.
882 884 *
... ... @@ -889,6 +891,9 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
889 891 const double& zetamin,
890 892 const double& timebin)
891 893 {
  894 + // Set constants
  895 + const double ehacorrmax = 10.0; //!< Maximum EHA correction factor
  896 +
892 897 // Debug
893 898 #if defined(G_DEBUG_COMPUTE_DRWS)
894 899 std::cout << "GCOMDris::compute_drws_phibar" << std::endl;
... ... @@ -940,9 +945,11 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
940 945 // Apply Earth horizon cut. Note that we need to correct later
941 946 // for this cut to remove the additional time variation coming
942 947 // from the cut.
  948 + #if defined(G_COMPUTE_DRWS_EHACORR)
943 949 if ((event->eha() - event->phibar()) < zetamin) {
944 950 continue;
945 951 }
  952 + #endif
946 953  
947 954 // If time exceeds time bin then add rates
948 955 while (event->time().secs() > (time + timebin)) {
... ... @@ -991,6 +998,21 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
991 998  
992 999 } // endif: appended remaining rates
993 1000  
  1001 + // Debug: write debug file
  1002 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1003 + std::ofstream debugfile;
  1004 + debugfile.open("compute_drws_phibar_tbin_rates.csv");
  1005 + debugfile.precision(12);
  1006 + for (int i = 0; i < times.size(); ++i) {
  1007 + debugfile << times[i];
  1008 + for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
  1009 + debugfile << "," << rates[iphibar][i];
  1010 + }
  1011 + debugfile << std::endl;
  1012 + }
  1013 + debugfile.close();
  1014 + #endif
  1015 +
994 1016 // Debug
995 1017 #if defined(G_DEBUG_COMPUTE_DRWS)
996 1018 std::cout << "Compute DRWs" << std::endl;
... ... @@ -1009,6 +1031,16 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1009 1031 tim.reduce(select.pulsar().validity());
1010 1032 }
1011 1033  
  1034 + // Debug: open debug files file
  1035 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1036 + std::ofstream debugfile1;
  1037 + std::ofstream debugfile2;
  1038 + debugfile1.open("compute_drws_phibar_sp_rates.csv");
  1039 + debugfile2.open("compute_drws_phibar_sp_ehacorr.csv");
  1040 + debugfile1.precision(12);
  1041 + debugfile2.precision(12);
  1042 + #endif
  1043 +
1012 1044 // Loop over Orbit Aspect Data
1013 1045 for (int i_oad = 0; i_oad < obs.oads().size(); ++i_oad) {
1014 1046  
... ... @@ -1031,10 +1063,13 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1031 1063 double time = oad.tstart().secs() + 0.5 * (oad.tstop() - oad.tstart());
1032 1064  
1033 1065 // Compute Phibar-dependent rates by interpolation of precomputed
1034   - // rate vectors
  1066 + // rate vectors. Make sure that rates never become negative.
1035 1067 std::vector<double> rate_oad(dri->nphibar());
1036 1068 for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
1037 1069 rate_oad[iphibar] = times.interpolate(time, rates[iphibar]);
  1070 + if (rate_oad[iphibar] < 0.0) {
  1071 + rate_oad[iphibar] = 0.0;
  1072 + }
1038 1073 }
1039 1074  
1040 1075 // Prepare Earth horizon angle computation. The celestial system
... ... @@ -1074,6 +1109,12 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1074 1109  
1075 1110 } // endfor: computed solid-angle-corrected geometry factors
1076 1111  
  1112 + // Debug: write time in debug files
  1113 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1114 + debugfile1 << time;
  1115 + debugfile2 << time;
  1116 + #endif
  1117 +
1077 1118 // Correct rates for Earth horizon cut
1078 1119 for (int iphibar = 0; iphibar < dri->nphibar(); ++iphibar) {
1079 1120  
... ... @@ -1092,14 +1133,31 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1092 1133 }
1093 1134 }
1094 1135  
1095   - // If something passes the Earth horizon cut then correct
1096   - // the rate
1097   - if (rates_sum_cut != 0.0) {
1098   - rate_oad[iphibar] *= rates_sum_all / rates_sum_cut;
  1136 + // Compute Earth horizon cut correction, limited to a maximum of ehacorrmax
  1137 + double ehacorr = (rates_sum_cut != 0.0) ? rates_sum_all / rates_sum_cut : ehacorrmax;
  1138 + if (ehacorr > ehacorrmax) {
  1139 + ehacorr = ehacorrmax;
1099 1140 }
1100 1141  
  1142 + // Correct rate
  1143 + #if defined(G_COMPUTE_DRWS_EHACORR)
  1144 + rate_oad[iphibar] *= ehacorr;
  1145 + #endif
  1146 +
  1147 + // Debug: write rate and correction in debug files
  1148 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1149 + debugfile1 << "," << rate_oad[iphibar];
  1150 + debugfile2 << "," << ehacorr;
  1151 + #endif
  1152 +
1101 1153 } // endfor: looped over Phibar
1102 1154  
  1155 + // Debug: write linefeeds in debug files
  1156 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1157 + debugfile1 << std::endl;
  1158 + debugfile2 << std::endl;
  1159 + #endif
  1160 +
1103 1161 // Loop over all DRWs
1104 1162 for (int i = 0; i < size(); ++i) {
1105 1163  
... ... @@ -1132,6 +1190,12 @@ void GCOMDris::compute_drws_phibar(const GCOMObservation&amp; obs,
1132 1190  
1133 1191 } // endfor: looped over Orbit Aspect Data
1134 1192  
  1193 + // Debug: close debug files
  1194 + #if defined(G_DEBUG_COMPUTE_DRWS)
  1195 + debugfile1.close();
  1196 + debugfile2.close();
  1197 + #endif
  1198 +
1135 1199 // Return
1136 1200 return;
1137 1201 }
... ...