gammalib / gammalib

Commit 8f0a5f82ea69cbc5b679596ffc881ca756a27b71

Authored by Jürgen Knödlseder
1 parent 9261f320

Add DRW handling and implement BGDLIXF method based on DRW instead of DRG

Inline Side-by-side
Showing 3 changed files with 286 additions and 8 deletions
inst/com/include/GCOMObservation.hpp
  View file @8f0a5f8
... ... @@ -70,6 +70,10 @@ public:
70 70 explicit GCOMObservation(const GXmlElement& xml);
71 71 GCOMObservation(const GCOMDri& dre,
72 72 const GCOMDri& drb,
  73 + const GCOMDri& drg,
  74 + const GCOMDri& drx);
  75 + GCOMObservation(const GCOMDri& dre,
  76 + const GCOMDri& drb,
73 77 const GCOMDri& drw,
74 78 const GCOMDri& drg,
75 79 const GCOMDri& drx);
... ... @@ -184,6 +188,11 @@ protected:
184 188 const int& navgr = 3,
185 189 const int& nincl = 15,
186 190 const int& nexcl = 0);
  191 + void compute_drb_bgdlixf(const GCOMDri& drm,
  192 + const int& nrunav = 3,
  193 + const int& navgr = 3,
  194 + const int& nincl = 15,
  195 + const int& nexcl = 0);
187 196 GSkyMap get_weighted_drg_map(void) const;
188 197 void get_bgdlixa_phibar_indices(const int& iphibar,
189 198 const int& nincl,
... ...
inst/com/pyext/GCOMObservation.i
  View file @8f0a5f8
... ... @@ -41,6 +41,10 @@ public:
41 41 explicit GCOMObservation(const GXmlElement& xml);
42 42 GCOMObservation(const GCOMDri& dre,
43 43 const GCOMDri& drb,
  44 + const GCOMDri& drg,
  45 + const GCOMDri& drx);
  46 + GCOMObservation(const GCOMDri& dre,
  47 + const GCOMDri& drb,
44 48 const GCOMDri& drw,
45 49 const GCOMDri& drg,
46 50 const GCOMDri& drx);
... ...
inst/com/src/GCOMObservation.cpp
  View file @8f0a5f8
... ... @@ -65,6 +65,8 @@ const GObservationRegistry g_obs_com_registry(&g_obs_com_seed);
65 65 "GCOMDri&, int&, int&, int&, int&)"
66 66 #define G_COMPUTE_DRB_BGDLIXE "GCOMObservation::compute_drb_bgdlixe("\
67 67 "GCOMDri&, int&, int&, int&, int&)"
  68 +#define G_COMPUTE_DRB_BGDLIXF "GCOMObservation::compute_drb_bgdlixf("\
  69 + "GCOMDri&, int&, int&, int&, int&)"
68 70  
69 71 /* __ Macros _____________________________________________________________ */
70 72  
... ... @@ -120,6 +122,50 @@ GCOMObservation::GCOMObservation(const GXmlElement& xml) : GObservation()
120 122 *
121 123 * @param[in] dre Event cube.
122 124 * @param[in] drb Background cube.
  125 + * @param[in] drg Geometry cube.
  126 + * @param[in] drx Exposure map.
  127 + *
  128 + * Creates COMPTEL observation from DRI instances.
  129 + *
  130 + * The method fixes the deadtime correction factor deadc to 0.965.
  131 + ***************************************************************************/
  132 +GCOMObservation::GCOMObservation(const GCOMDri& dre,
  133 + const GCOMDri& drb,
  134 + const GCOMDri& drg,
  135 + const GCOMDri& drx) : GObservation()
  136 +{
  137 + // Initialise members
  138 + init_members();
  139 +
  140 + // Store DRIs
  141 + m_events = new GCOMEventCube(dre);
  142 + m_drb = drb;
  143 + m_drg = drg;
  144 + m_drx = drx;
  145 +
  146 + // Set attributes
  147 + m_obs_id = 0;
  148 + m_ontime = m_events->gti().ontime();
  149 + m_deadc = 0.965;
  150 + m_livetime = m_deadc * m_ontime;
  151 + m_ewidth = m_events->emax().MeV() - m_events->emin().MeV();
  152 + m_name = "unknown";
  153 + m_drename = "";
  154 + m_drbname = "";
  155 + m_drwname = "";
  156 + m_drgname = "";
  157 + m_drxname = "";
  158 +
  159 + // Return
  160 + return;
  161 +}
  162 +
  163 +
  164 +/***********************************************************************//**
  165 + * @brief Binned observation DRI constructor
  166 + *
  167 + * @param[in] dre Event cube.
  168 + * @param[in] drb Background cube.
123 169 * @param[in] drw Weighting cube.
124 170 * @param[in] drg Geometry cube.
125 171 * @param[in] drx Exposure map.
... ... @@ -901,15 +947,20 @@ GCOMDri GCOMObservation::drm(const GModels& models) const
901 947 /***********************************************************************//**
902 948 * @brief Compute DRB cube
903 949 *
904   - * @param[in] method Background method (PHINOR, BGDLIXA or BGDLIXE).
  950 + * @param[in] method Background method (PHINOR, BGDLIXA, BGDLIXE or BGDLIXF).
905 951 * @param[in] drm DRM cube.
906   - * @param[in] nrunav BGDLIXA: number of bins used for running average.
907   - * @param[in] navgr BGDLIXA: number of bins used for averaging.
908   - * @param[in] nincl BGDLIXA: number of Phibar layers to include.
909   - * @param[in] nexcl BGDLIXA: number of Phibar layers to exclude.
910   - *
911   - * Computes a COMPTEL DRB cube using either the PHINOR or BGDLIXA method.
912   - * See the protected methods compute_drb_phinor() and compute_drb_bgdlixa()
  952 + * @param[in] nrunav BGDLIXn: number of bins used for running average.
  953 + * @param[in] navgr BGDLIXn: number of bins used for averaging.
  954 + * @param[in] nincl BGDLIXn: number of Phibar layers to include.
  955 + * @param[in] nexcl BGDLIXn: number of Phibar layers to exclude.
  956 + *
  957 + * Computes a COMPTEL DRB cube using either the PHINOR or BGDLIX method.
  958 + * There are three variants of the BGFLIX method (BGDLIXA, BGDLIXE and
  959 + * BGDLIXF). See the protected methods
  960 + * compute_drb_phinor(),
  961 + * compute_drb_bgdlixa(),
  962 + * compute_drb_bgdlixe(), and
  963 + * compute_drb_bgdlixf()
913 964 * for more information.
914 965 ***************************************************************************/
915 966 void GCOMObservation::compute_drb(const std::string& method,
... ... @@ -929,6 +980,9 @@ void GCOMObservation::compute_drb(const std::string& method,
929 980 else if (method == "BGDLIXE") {
930 981 compute_drb_bgdlixe(drm, nrunav, navgr, nincl, nexcl);
931 982 }
  983 + else if (method == "BGDLIXF") {
  984 + compute_drb_bgdlixf(drm, nrunav, navgr, nincl, nexcl);
  985 + }
932 986 else {
933 987 std::string msg = "Unknown background method \""+method+"\". "
934 988 "Specify either \"PHINOR\", \"BGDLIXA\" or "
... ... @@ -2080,6 +2134,217 @@ void GCOMObservation::compute_drb_bgdlixe(const GCOMDri& drm,
2080 2134  
2081 2135  
2082 2136 /***********************************************************************//**
  2137 + * @brief Compute DRB cube using BGDLIXF method
  2138 + *
  2139 + * @param[in] drm DRM cube.
  2140 + * @param[in] nrunav Number of bins used for running average.
  2141 + * @param[in] navgr Number of bins used for averaging.
  2142 + * @param[in] nincl Number of Phibar layers to include.
  2143 + * @param[in] nexcl Number of Phibar layers to exclude.
  2144 + *
  2145 + * @exception GException::invalid_value
  2146 + * Observation does not contain an event cube
  2147 + * @exception GException::invalid_argument
  2148 + * DRM cube is incompatible with DRE
  2149 + *
  2150 + * Computes a DRB cube using the BGDLIXF method. This method differs from
  2151 + * the BGDLIXE method in the DRW instead of the DRG is used for background
  2152 + * model computation.
  2153 + ***************************************************************************/
  2154 +void GCOMObservation::compute_drb_bgdlixf(const GCOMDri& drm,
  2155 + const int& nrunav,
  2156 + const int& navgr,
  2157 + const int& nincl,
  2158 + const int& nexcl)
  2159 +{
  2160 + // Extract COMPTEL event cube
  2161 + const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(this->events());
  2162 + if (cube == NULL) {
  2163 + std::string msg = "Observation \""+this->name()+"\" ("+this->id()+") "
  2164 + "does not contain a COMPTEL event cube. Please "
  2165 + "specify a COMPTEL observation containing and event "
  2166 + "cube.";
  2167 + throw GException::invalid_value(G_COMPUTE_DRB_BGDLIXF, msg);
  2168 + }
  2169 +
  2170 + // Check DRM
  2171 + if (!check_dri(drm)) {
  2172 + std::string msg = "Specified DRM cube is incompatible with DRE. Please "
  2173 + "specify a DRM with a data-space definition that is "
  2174 + "identical to that of the DRE.";
  2175 + throw GException::invalid_argument(G_COMPUTE_DRB_BGDLIXE, msg);
  2176 + }
  2177 +
  2178 + // Initialise DRB and scratch DRI by cloning DRW
  2179 + m_drb = m_drw;
  2180 + GCOMDri dri = m_drw;
  2181 +
  2182 + // Get references to relevant sky maps
  2183 + const GSkyMap& map_dre = cube->dre().map();
  2184 + const GSkyMap& map_drm = drm.map();
  2185 + const GSkyMap& map_drw = m_drw.map();
  2186 + GSkyMap& map_drb = const_cast<GSkyMap&>(m_drb.map());
  2187 + GSkyMap& map_dri = const_cast<GSkyMap&>(dri.map());
  2188 +
  2189 + // Get data space dimensions
  2190 + int nchi = m_drw.nchi();
  2191 + int npsi = m_drw.npsi();
  2192 + int nphibar = m_drw.nphibar();
  2193 + int npix = nchi * npsi;
  2194 +
  2195 + // Precompute half running average lengths
  2196 + int navgr2 = int(navgr/2);
  2197 +
  2198 + // Initialise DRB and scratch DRI
  2199 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  2200 + for (int ipix = 0; ipix < npix; ++ipix) {
  2201 + map_drb(ipix, iphibar) = 0.0;
  2202 + map_dri(ipix, iphibar) = 0.0;
  2203 + }
  2204 + }
  2205 +
  2206 + // Phibar normalise DRW
  2207 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  2208 + double sum_dre = 0.0;
  2209 + double sum_drm = 0.0;
  2210 + double sum_drw = 0.0;
  2211 + for (int ipix = 0; ipix < npix; ++ipix) {
  2212 + sum_dre += map_dre(ipix, iphibar);
  2213 + sum_drm += map_drm(ipix, iphibar);
  2214 + sum_drw += map_drw(ipix, iphibar);
  2215 + }
  2216 + if (sum_drw > 0.0) {
  2217 + double norm = (sum_dre - sum_drm) / sum_drw;
  2218 + for (int ipix = 0; ipix < npix; ++ipix) {
  2219 + map_dri(ipix, iphibar) = map_drw(ipix, iphibar) * norm;
  2220 + }
  2221 + }
  2222 + }
  2223 +
  2224 + // Fit Phibar normalised DRW to DRE-DRM
  2225 + for (int ichi = 0; ichi < nchi; ++ichi) {
  2226 +
  2227 + // Compute running average window in Chi
  2228 + int kchi_min = ichi - navgr2;
  2229 + int kchi_max = ichi + navgr2;
  2230 + if (kchi_min < 0) {
  2231 + kchi_min = 0;
  2232 + }
  2233 + if (kchi_max >= nchi) {
  2234 + kchi_max = nchi - 1;
  2235 + }
  2236 +
  2237 + // Loop over Psi pixels
  2238 + for (int ipsi = 0; ipsi < npsi; ++ipsi) {
  2239 +
  2240 + // Compute running average window in Psi
  2241 + int kpsi_min = ipsi - navgr2;
  2242 + int kpsi_max = ipsi + navgr2;
  2243 + if (kpsi_min < 0) {
  2244 + kpsi_min = 0;
  2245 + }
  2246 + if (kpsi_max >= npsi) {
  2247 + kpsi_max = npsi - 1;
  2248 + }
  2249 +
  2250 + // Loop over Phibar layers
  2251 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  2252 +
  2253 + // Get Phibar index range for sums
  2254 + int ksel1 = 0;
  2255 + int kex1 = 0;
  2256 + int kex2 = 0;
  2257 + int ksel2 = 0;
  2258 + get_bgdlixa_phibar_indices(iphibar, nincl, nexcl,
  2259 + &ksel1, &kex1, &kex2, &ksel2);
  2260 +
  2261 + // Initialise sums
  2262 + double sum_dre = 0.0;
  2263 + double sum_drm = 0.0;
  2264 + double sum_dri = 0.0;
  2265 +
  2266 + // Compute running average
  2267 + for (int kchi = kchi_min; kchi <= kchi_max; ++kchi) {
  2268 + for (int kpsi = kpsi_min; kpsi <= kpsi_max; ++kpsi) {
  2269 +
  2270 + // Get index
  2271 + int kpix = kchi + kpsi * nchi;
  2272 +
  2273 + // Take sum over [ksel1, kex1]
  2274 + if (kex1 >= 0) {
  2275 + for (int kphibar = ksel1; kphibar <= kex1; ++kphibar) {
  2276 + sum_dre += map_dre(kpix, kphibar);
  2277 + sum_drm += map_drm(kpix, kphibar);
  2278 + sum_dri += map_dri(kpix, kphibar);
  2279 + }
  2280 + }
  2281 +
  2282 + // Take sum over [kxe2, ksel2]
  2283 + if (kex2 < nphibar) {
  2284 + for (int kphibar = kex2; kphibar <= ksel2; ++kphibar) {
  2285 + sum_dre += map_dre(kpix, kphibar);
  2286 + sum_drm += map_drm(kpix, kphibar);
  2287 + sum_dri += map_dri(kpix, kphibar);
  2288 + }
  2289 + }
  2290 +
  2291 + } // endfor: looped over kpsi
  2292 + } // endfor: looped over kchi
  2293 +
  2294 + // Renormalize DRI locally to DRE-DRM
  2295 + int ipix = ichi + ipsi * nchi;
  2296 + if (sum_dri > 0.0) {
  2297 + double norm = (sum_dre - sum_drm) / sum_dri;
  2298 + map_drb(ipix, iphibar) = map_dri(ipix, iphibar) * norm;
  2299 + }
  2300 + else {
  2301 + map_drb(ipix, iphibar) = 0.0;
  2302 + }
  2303 +
  2304 + } // endfor: looped over Phibar layers
  2305 +
  2306 + } // endfor: looped over Psi pixels
  2307 +
  2308 + } // endfor: looped over Chi pixels
  2309 +
  2310 + // Phibar normalise DRB
  2311 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  2312 + double sum_dre = 0.0;
  2313 + double sum_drm = 0.0;
  2314 + double sum_drb = 0.0;
  2315 + for (int ipix = 0; ipix < npix; ++ipix) {
  2316 + sum_dre += map_dre(ipix, iphibar);
  2317 + sum_drm += map_drm(ipix, iphibar);
  2318 + sum_drb += map_drb(ipix, iphibar);
  2319 + }
  2320 + if (sum_drb > 0.0) {
  2321 + double norm = (sum_dre - sum_drm) / sum_drb;
  2322 + for (int ipix = 0; ipix < npix; ++ipix) {
  2323 + map_drb(ipix, iphibar) *= norm;
  2324 + }
  2325 + }
  2326 + else {
  2327 + for (int ipix = 0; ipix < npix; ++ipix) {
  2328 + map_drb(ipix, iphibar) = 0.0;
  2329 + }
  2330 + }
  2331 + }
  2332 +
  2333 + // Make sure that DRB is non-negative
  2334 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  2335 + for (int ipix = 0; ipix < npix; ++ipix) {
  2336 + if (map_drb(ipix, iphibar) < 0.0) {
  2337 + map_drb(ipix, iphibar) = 0.0;
  2338 + }
  2339 + }
  2340 + }
  2341 +
  2342 + // Return
  2343 + return;
  2344 +}
  2345 +
  2346 +
  2347 +/***********************************************************************//**
2083 2348 * @brief Return weighted DRG map
2084 2349 *
2085 2350 * @return Weighted DRG map.
... ...