From a276a1d31ef988c8a1bcd6603b4e81474fd7a8fa Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 3 Mar 2021 22:30:49 +0100 Subject: [PATCH] Implement COMPTEL background model computation (#3557) --- ChangeLog | 1 + NEWS | 8 ++++++-- doc/source/admin/release_history/2.0.rst | 2 ++ inst/com/include/GCOMObservation.hpp | 19 +++++++++++++++++++ inst/com/pyext/GCOMObservation.i | 6 ++++++ inst/com/src/GCOMObservation.cpp | 514 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- 6 files changed, 547 insertions(+), 3 deletions(-) diff --git a/ChangeLog b/ChangeLog index 5083a96..04ac6ad 100644 --- a/ChangeLog +++ b/ChangeLog @@ -5,6 +5,7 @@ Store DRB, DRG and DRX as GCOMDri in GCOMObservation (#3560) Implement COMPTEL response for diffuse models (#3558) + Implement COMPTEL background model computation (#3557) Fix reading of scaled FITS columns, remove GFitsTableCol::scale() method (#3553) Add gammalib::compare() function (#3552) Add DRI constructor to GCOMObservation class diff --git a/NEWS b/NEWS index c5b5adf..773cb96 100644 --- a/NEWS +++ b/NEWS @@ -77,6 +77,7 @@ The following methods have been added: - GCOMObservation::drgname() - GCOMObservation::drxname() - GCOMObservation::response(GCOMResponse&) +- GCOMObservation::compute_drb() - GCOMDri::tof_correction() - GCOMDri::num_superpackets() - GCOMDri::num_used_superpackets() @@ -280,8 +281,8 @@ spatial model gradients (#3203). Store DRB, DRG and DRX as instances of GCOMDri instead of GSkyMap in GCOMObservation (#3560). -Remove internal response caching since caching is handled by the base class response -computation methods (#2973). +Added GCOMObservation::compute_drb() method to compute background models for COMPTEL. +So far the PHINOR and BGDLIXA background methods were implemented (#3557). Added method to set COMPTEL response using GCOMResponse instance. @@ -301,6 +302,9 @@ Added methods to access DRI filenames in GCOMObservation. Implemented vectorised response computation in GCOMResponse::irf_ptsrc() method (#3501). +Remove internal response caching since caching is handled by the base class response +computation methods (#2973). + 22. INTEGRAL/SPI interface -------------------------- diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 8aa79cc..555409a 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -33,6 +33,8 @@ Improvements * [`3558 `_] - Implement COMPTEL response for diffuse models +* [`3557 `_] - + Implement COMPTEL background model computation * [`3552 `_] - Add ``gammalib::compare()`` function * Add DRI constructor to ``GCOMObservation`` class diff --git a/inst/com/include/GCOMObservation.hpp b/inst/com/include/GCOMObservation.hpp index f5600fc..23898f6 100644 --- a/inst/com/include/GCOMObservation.hpp +++ b/inst/com/include/GCOMObservation.hpp @@ -131,6 +131,12 @@ public: void drbname(const GFilename& drbname); void drgname(const GFilename& drgname); void drxname(const GFilename& drxname); + void compute_drb(const std::string& method, + const GCOMDri& drm, + const int& nrunav = 3, + const int& navgr = 3, + const int& nincl = 13, + const int& nexcl = 0); protected: // Protected methods @@ -144,6 +150,19 @@ protected: bool check_dri(const GCOMDri& map) const; void read_attributes(const GFitsHDU* hdu); void write_attributes(GFitsHDU* hdu) const; + void compute_drb_phinor(const GCOMDri& drm); + void compute_drb_bgdlixa(const GCOMDri& drm, + const int& nrunav = 3, + const int& navgr = 3, + const int& nincl = 13, + const int& nexcl = 0); + void get_bgdlixa_phibar_indices(const int& iphibar, + const int& nincl, + const int& nexcl, + int* isel1, + int* iex1, + int* iex2, + int* isel2) const; // Protected members std::string m_instrument; //!< Instrument name diff --git a/inst/com/pyext/GCOMObservation.i b/inst/com/pyext/GCOMObservation.i index 1cc3732..9fd94bc 100644 --- a/inst/com/pyext/GCOMObservation.i +++ b/inst/com/pyext/GCOMObservation.i @@ -99,6 +99,12 @@ public: void drbname(const GFilename& drbname); void drgname(const GFilename& drgname); void drxname(const GFilename& drxname); + void compute_drb(const std::string& method, + const GCOMDri& drm, + const int& nrunav = 3, + const int& navgr = 3, + const int& nincl = 13, + const int& nexcl = 0); }; diff --git a/inst/com/src/GCOMObservation.cpp b/inst/com/src/GCOMObservation.cpp index 20d87c9..b785c77 100644 --- a/inst/com/src/GCOMObservation.cpp +++ b/inst/com/src/GCOMObservation.cpp @@ -57,7 +57,9 @@ const GObservationRegistry g_obs_com_registry(&g_obs_com_seed); #define G_LOAD_DRB "GCOMObservation::load_drb(GFilename&)" #define G_LOAD_DRG "GCOMObservation::load_drg(GFilename&)" #define G_DRM "GCOMObservation::drm(GModels&)" -#define G_ADD_DRM "GCOMObservation::add_drm(GSkyModel&)" +#define G_COMPUTE_DRB "GCOMObservation::compute_drb(std::string&, GCOMDri&, "\ + "int&, int&, int&, int&" +#define G_COMPUTE_DRB_PHINOR "GCOMObservation::compute_drb_phinor(GCOMDri&)" /* __ Macros _____________________________________________________________ */ @@ -731,6 +733,45 @@ GCOMDri GCOMObservation::drm(const GModels& models) const /***********************************************************************//** + * @brief Compute DRB cube + * + * @param[in] method Background method (PHINOR or BGDLIXA). + * @param[in] drm DRM cube. + * @param[in] nrunav BGDLIXA: number of bins used for running average. + * @param[in] navgr BGDLIXA: number of bins used for averaging. + * @param[in] nincl BGDLIXA: number of Phibar layers to include. + * @param[in] nexcl BGDLIXA: number of Phibar layers to exclude. + * + * Computes a COMPTEL DRB cube using either the PHINOR or BGDLIXA method. + * See the protected methods compute_drb_phinor() and compute_drb_bgdlixa() + * for more information. + ***************************************************************************/ +void GCOMObservation::compute_drb(const std::string& method, + const GCOMDri& drm, + const int& nrunav, + const int& navgr, + const int& nincl, + const int& nexcl) +{ + // Branch to relevant method based on specified background method + if (method == "PHINOR") { + compute_drb_phinor(drm); + } + else if (method == "BGDLIXA") { + compute_drb_bgdlixa(drm, nrunav, navgr, nincl, nexcl); + } + else { + std::string msg = "Unknown background method \""+method+"\". " + "Specify either \"PHINOR\" or \"BGDLIXA\"."; + throw GException::invalid_argument(G_COMPUTE_DRB, msg); + } + + // Return + return; +} + + +/***********************************************************************//** * @brief Print observation information * * @param[in] chatter Chattiness. @@ -1135,3 +1176,474 @@ void GCOMObservation::write_attributes(GFitsHDU* hdu) const // Return return; } + + +/***********************************************************************//** + * @brief Compute DRB cube using PHINOR method + * + * @param[in] drm DRM cube. + * + * @exception GException::invalid_value + * Observation does not contain an event cube + * @exception GException::invalid_argument + * DRM cube is incompatible with DRE + ***************************************************************************/ +void GCOMObservation::compute_drb_phinor(const GCOMDri& drm) +{ + // Extract COMPTEL event cube + const GCOMEventCube* cube = dynamic_cast(this->events()); + if (cube == NULL) { + std::string msg = "Observation \""+this->name()+"\" ("+this->id()+") " + "does not contain a COMPTEL event cube. Please " + "specify a COMPTEL observation containing and event " + "cube."; + throw GException::invalid_value(G_COMPUTE_DRB_PHINOR, msg); + } + + // Check DRM + if (!check_dri(drm)) { + std::string msg = "Specified DRM cube is incompatible with DRE. Please " + "specify a DRM with a data-space definition that is " + "identical to that of the DRE."; + throw GException::invalid_argument(G_COMPUTE_DRB_PHINOR, msg); + } + + // Initialise DRB by cloning DRG + m_drb = m_drg; + + // Get DRI sky maps + const GSkyMap& map_dre = cube->dre().map(); + const GSkyMap& map_drm = drm.map(); + const GSkyMap& map_drg = m_drg.map(); + GSkyMap& map_drb = const_cast(m_drb.map()); + + // Get data space dimensions + int npix = m_drg.map().npix(); + int nphibar = m_drg.nphibar(); + + // Phibar normalise DRB + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + double sum_dre = 0.0; + double sum_drm = 0.0; + double sum_drg = 0.0; + for (int ipix = 0; ipix < npix; ++ipix) { + sum_dre += map_dre(ipix, iphibar); + sum_drm += map_drm(ipix, iphibar); + sum_drg += map_drg(ipix, iphibar); + } + if (sum_drg > 0.0) { + double norm = (sum_dre - sum_drm) / sum_drg; + for (int ipix = 0; ipix < npix; ++ipix) { + map_drb(ipix, iphibar) *= norm; + } + } + else { + for (int ipix = 0; ipix < npix; ++ipix) { + map_drb(ipix, iphibar) = 0.0; + } + } + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Compute DRB cube using BGDLIXA method + * + * @param[in] drm DRM cube. + * @param[in] nrunav Number of bins used for running average. + * @param[in] navgr Number of bins used for averaging. + * @param[in] nincl Number of Phibar layers to include. + * @param[in] nexcl Number of Phibar layers to exclude. + * + * @exception GException::invalid_value + * Observation does not contain an event cube + * @exception GException::invalid_argument + * DRM cube is incompatible with DRE + * + * Computes a DRB cube using the BGDLIXA method that is documented in Rob + * van Dijk's PhD thesis. The revelant equations from the thesis that are + * implemented here are 3.12, 3.12 and 3.14. + ***************************************************************************/ +void GCOMObservation::compute_drb_bgdlixa(const GCOMDri& drm, + const int& nrunav, + const int& navgr, + const int& nincl, + const int& nexcl) +{ + // Extract COMPTEL event cube + const GCOMEventCube* cube = dynamic_cast(this->events()); + if (cube == NULL) { + std::string msg = "Observation \""+this->name()+"\" ("+this->id()+") " + "does not contain a COMPTEL event cube. Please " + "specify a COMPTEL observation containing and event " + "cube."; + throw GException::invalid_value(G_COMPUTE_DRB_PHINOR, msg); + } + + // Check DRM + if (!check_dri(drm)) { + std::string msg = "Specified DRM cube is incompatible with DRE. Please " + "specify a DRM with a data-space definition that is " + "identical to that of the DRE."; + throw GException::invalid_argument(G_COMPUTE_DRB_PHINOR, msg); + } + + // Initialise DRB and scratch DRI by cloning DRG + m_drb = m_drg; + GCOMDri dri = m_drg; + + // Get references to relevant sky maps + const GSkyMap& map_dre = cube->dre().map(); + const GSkyMap& map_drm = drm.map(); + const GSkyMap& map_drg = m_drg.map(); + GSkyMap& map_drb = const_cast(m_drb.map()); + GSkyMap& map_dri = const_cast(dri.map()); + + // Get data space dimensions + int nchi = m_drg.nchi(); + int npsi = m_drg.npsi(); + int nphibar = m_drg.nphibar(); + int npix = nchi * npsi; + int nsize = npix * nphibar; + + // Precompute half running average lengths + int navgr2 = int(navgr/2); + int nexcl2 = int(nexcl/2); + int nincl2 = int(nincl/2); + + // Initialise DRB and scratch DRI + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + for (int ipix = 0; ipix < npix; ++ipix) { + map_drb(ipix, iphibar) = 0.0; + map_dri(ipix, iphibar) = 0.0; + } + } + + // Phibar normalise DRB (Equation 3.12 in Rob van Dijk's PhD thesis). + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + double sum_dre = 0.0; + double sum_drm = 0.0; + double sum_drg = 0.0; + for (int ipix = 0; ipix < npix; ++ipix) { + sum_dre += map_dre(ipix, iphibar); + sum_drm += map_drm(ipix, iphibar); + sum_drg += map_drg(ipix, iphibar); + } + if (sum_drg > 0.0) { + double norm = (sum_dre - sum_drm) / sum_drg; + for (int ipix = 0; ipix < npix; ++ipix) { + map_dri(ipix, iphibar) = map_drg(ipix, iphibar) * norm; + } + } + } + + // Do 3D running average of DRE and Phibar normalised DRG (Equation 3.13 + // in Rob van Dijk's PhD thesis). The result is a Phibar normalised DRG + // that is normalised over the 3D region to the DRE. The 3D region is all + // Phibar layers and [-nrunav,+nrunav] Chi/Psi pixels around the pixel of + // consideration. + if (nrunav >= 1) { + + // Loop over Chi pixels + for (int ichi = 0; ichi < nchi; ++ichi) { + + // Compute running average window in Chi + int kchi_min = ichi - nrunav; + int kchi_max = ichi + nrunav; + if (kchi_min < 0) { + kchi_min = 0; + } + if (kchi_max >= nchi) { + kchi_max = nchi - 1; + } + + // Loop over Psi pixels + for (int ipsi = 0; ipsi < npsi; ++ipsi) { + + // Compute running average window in Psi + int kpsi_min = ipsi - nrunav; + int kpsi_max = ipsi + nrunav; + if (kpsi_min < 0) { + kpsi_min = 0; + } + if (kpsi_max >= npsi) { + kpsi_max = npsi - 1; + } + + // Initialise sums + double sum_dre = 0.0; + double sum_drm = 0.0; + double sum_dri = 0.0; + + // Compute running average + for (int kchi = kchi_min; kchi <= kchi_max; ++kchi) { + for (int kpsi = kpsi_min; kpsi <= kpsi_max; ++kpsi) { + int kpix = kchi + kpsi * nchi; + for (int kphibar = 0; kphibar < nphibar; ++kphibar) { + if (map_drg(kpix, kphibar) != 0.0) { + sum_dre += map_dre(kpix, kphibar); + sum_drm += map_drm(kpix, kphibar); + sum_dri += map_dri(kpix, kphibar); + } + } + } + } + + // Renormalise scratch array + if (sum_dri != 0.0) { + int ipix = ichi + ipsi * nchi; + double norm = (sum_dre - sum_drm) / sum_dri; + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + map_dri(ipix, iphibar) *= norm; + } + } + + } // endfor: looped over Phi + + } // endfor: looped over Chi + + } // endif: running averaging was requested + + // First part of Equation (3.14) in Rob van Dijk's PhD thesis. + // (DRB = B^0C_L / sum B^0C_L). The DRB is pre-computed by adjusting the + // scratch DRI over nincl Phibar layers. + int ksel1 = 0; + int kex1 = 0; + int kex2 = 0; + int ksel2 = 0; + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + + // Get Phibar index range for sums + get_bgdlixa_phibar_indices(iphibar, nincl, nexcl, + &ksel1, &kex1, &kex2, &ksel2); + + // Loop over Chi/Psi pixels + for (int ipix = 0; ipix < npix; ++ipix) { + + // Initialise DRB value + map_drb(ipix, iphibar) = 0.0; + + // Continue only if scratch DRI is non-zero + if (map_dri(ipix, iphibar) != 0.0) { + + // Initialise sum + double sum_dri = 0.0; + + // Take sum over [ksel1, kex1] + if (kex1 >= 0) { + for (int kphibar = ksel1; kphibar <= kex1; ++kphibar) { + sum_dri += map_dri(ipix, kphibar); + } + } + + // Take sum over [kxe2, ksel2] + if (kex2 < nphibar) { + for (int kphibar = kex2; kphibar <= ksel2; ++kphibar) { + sum_dri += map_dri(ipix, kphibar); + } + } + + // If sum is not zero then set DRB + if (sum_dri != 0.0) { + map_drb(ipix, iphibar) = map_dri(ipix, iphibar) / sum_dri; + } + + } // endif: scratch DRI was non-zero + + } // endfor: looped over Chi/Psi pixels + + } // endfor: looped over Phibar + + // Second part of Equation (3.14) in Rob van Dijk's PhD thesis that + // corresponds to G * [E / G]^s that will be stored in scratch DRI. + // + // NOTES: + // - we replaced DRE by DRE-DRM in the nominator of the running + // average since and source component should of course be + // subtracted + // - we implemented an alternative running average computation where + // the ratio between the DRE and DRG sums is taken. This avoids + // the divergence of the ratio at the edge of the DRG. + for (int ichi = 0; ichi < nchi; ++ichi) { + + // Compute running average window in Chi + int kchi_min = ichi - navgr2; + int kchi_max = ichi + navgr2; + if (kchi_min < 0) { + kchi_min = 0; + } + if (kchi_max >= nchi) { + kchi_max = nchi - 1; + } + + // Loop over Psi pixels + for (int ipsi = 0; ipsi < npsi; ++ipsi) { + + // Compute running average window in Psi + int kpsi_min = ipsi - navgr2; + int kpsi_max = ipsi + navgr2; + if (kpsi_min < 0) { + kpsi_min = 0; + } + if (kpsi_max >= npsi) { + kpsi_max = npsi - 1; + } + + // Loop over Phibar layers + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + + // Initialise sums + double sum_dre = 0.0; + double sum_drm = 0.0; + double sum_drg = 0.0; + + // Compute running average + for (int kchi = kchi_min; kchi <= kchi_max; ++kchi) { + for (int kpsi = kpsi_min; kpsi <= kpsi_max; ++kpsi) { + int kpix = kchi + kpsi * nchi; + sum_dre += map_dre(kpix, iphibar); + sum_drm += map_drm(kpix, iphibar); + sum_drg += map_drg(kpix, iphibar); + } + } + + // Multiply average by DRG and store it in scratch DRI + int ipix = ichi + ipsi * nchi; + if (sum_drg > 0.0) { + double norm = (sum_dre - sum_drm) / sum_drg; + map_dri(ipix, iphibar) = map_drg(ipix, iphibar) * norm; + } + else { + map_dri(ipix, iphibar) = 0.0; + } + + } // endfor: looped over Phibar layers + + } // endfor: looped over Psi pixels + + } // endfor: looped over Chi pixels + + // Third part of Equation (3.14) in Rob van Dijk's PhD thesis that sums + // G * [E / G]^s over a subset of Phibar layers (the first parenthesis) + // and multplies with B^0C_L / sum B^0C_L + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + + // Get Phibar index range for sums + get_bgdlixa_phibar_indices(iphibar, nincl, nexcl, + &ksel1, &kex1, &kex2, &ksel2); + + // Loop over Chi/Psi pixels + for (int ipix = 0; ipix < npix; ++ipix) { + + // If DRB is not empty then correct it + if (map_drb(ipix, iphibar) != 0.0) { + + // Initialise DRI sum + double sum_dri = 0.0; + + // Take sum over [ksel1, kex1] + if (kex1 >= 0) { + for (int kphibar = ksel1; kphibar <= kex1; ++kphibar) { + sum_dri += map_dri(ipix, kphibar); + } + } + + // Take sum over [kxe2, ksel2] + if (kex2 < nphibar) { + for (int kphibar = kex2; kphibar <= ksel2; ++kphibar) { + sum_dri += map_dri(ipix, kphibar); + } + } + + // Correct DRB + map_drb(ipix, iphibar) *= sum_dri; + + } // endif: DRB was not empty + + } // endfor: looped over Chi/Psi pixels + + } // endfor: looped over Phibar + + // Return + return; +} + + +/***********************************************************************//** + * @brief Compute Phibar index range for BGDLIXA background method + * + * @param[in] iphibar Phibar layer index. + * @param[in] nincl Number of Phibar layers to include. + * @param[in] nexcl Number of Phibar layers to exclude. + * @param[out] isel1 Start index for first sum. + * @param[out] iex1 Stop index for first sum. + * @param[out] iex2 Start index for second sum. + * @param[out] isel2 Stop index for second sum. + * + * This method is a helper method for the BGDLIXA background method. It + * computes the Phibar index range for the third step of the background + * model computation. The third step corresponds to Equation (3.14) in Rob + * van Dijk's PhD thesis. + * + * The Phibar sum will be taken over the index ranges [isel1, iex1] and + * [ixe2, isel2]. + ***************************************************************************/ +void GCOMObservation::get_bgdlixa_phibar_indices(const int& iphibar, + const int& nincl, + const int& nexcl, + int* isel1, + int* iex1, + int* iex2, + int* isel2) const +{ + // Get number of Phibar layers + int nphibar = m_drg.nphibar(); + int imax = nphibar - 1; + + // Precompute half running average lengths + int nexcl2 = int(nexcl/2); + int nincl2 = int(nincl/2); + + // Compute Phibar index range without respecting boundaries + *iex1 = iphibar - nexcl2 - 1; + *iex2 = iphibar + nexcl2 + 1; + *isel1 = iphibar - nincl2; + *isel2 = iphibar + nincl2; + + // If no Phibar layers are to be excluded then make sure that the + // Phibar layer range is continuous + if (nexcl == 0) { + (*iex2)--; + } + + // Make sure that start index of first sum and stop index of second + // some are within the validity range + if (*isel1 < 0) { + *isel1 = 0; + } + if (*isel2 > imax) { + *isel2 = imax; + } + + // Make sure that the sums use always ncincl Phibar layers, also near + // the edges + if (*isel1 == 0) { + *isel2 = nincl - 1; + if (nincl > imax) { + *isel2 = imax; + } + + } + if (*isel2 == imax) { + *isel1 = nphibar - nincl; + if (*isel1 < 0) { + *isel1 = 0; + } + } + + // Return + return; +} -- libgit2 0.21.4