Sergio Hernández Cadena / gammalib

Commit 07fad2df2527d8037d2dde6cfe38f361b8d46e0f

Authored by Jürgen Knödlseder
1 parent 48b8040c

Multiply DRG pixels with solidangle before background computation

Inline Side-by-side
Showing 3 changed files with 128 additions and 59 deletions
inst/com/include/GCOMObservation.hpp
  View file @07fad2d
... ... @@ -140,29 +140,30 @@ public:
140 140  
141 141 protected:
142 142 // Protected methods
143   - void init_members(void);
144   - void copy_members(const GCOMObservation& obs);
145   - void free_members(void);
146   - void load_dre(const GFilename& drename);
147   - void load_drb(const GFilename& drbname);
148   - void load_drg(const GFilename& drgname);
149   - void load_drx(const GFilename& drxname);
150   - bool check_dri(const GCOMDri& map) const;
151   - void read_attributes(const GFitsHDU* hdu);
152   - void write_attributes(GFitsHDU* hdu) const;
153   - void compute_drb_phinor(const GCOMDri& drm);
154   - void compute_drb_bgdlixa(const GCOMDri& drm,
155   - const int& nrunav = 3,
156   - const int& navgr = 3,
157   - const int& nincl = 13,
158   - const int& nexcl = 0);
159   - void get_bgdlixa_phibar_indices(const int& iphibar,
160   - const int& nincl,
161   - const int& nexcl,
162   - int* isel1,
163   - int* iex1,
164   - int* iex2,
165   - int* isel2) const;
  143 + void init_members(void);
  144 + void copy_members(const GCOMObservation& obs);
  145 + void free_members(void);
  146 + void load_dre(const GFilename& drename);
  147 + void load_drb(const GFilename& drbname);
  148 + void load_drg(const GFilename& drgname);
  149 + void load_drx(const GFilename& drxname);
  150 + bool check_dri(const GCOMDri& map) const;
  151 + void read_attributes(const GFitsHDU* hdu);
  152 + void write_attributes(GFitsHDU* hdu) const;
  153 + void compute_drb_phinor(const GCOMDri& drm);
  154 + void compute_drb_bgdlixa(const GCOMDri& drm,
  155 + const int& nrunav = 3,
  156 + const int& navgr = 3,
  157 + const int& nincl = 13,
  158 + const int& nexcl = 0);
  159 + GSkyMap get_weighted_drg_map(void) const;
  160 + void get_bgdlixa_phibar_indices(const int& iphibar,
  161 + const int& nincl,
  162 + const int& nexcl,
  163 + int* isel1,
  164 + int* iex1,
  165 + int* iex2,
  166 + int* isel2) const;
166 167  
167 168 // Protected members
168 169 std::string m_instrument; //!< Instrument name
... ...
inst/com/src/GCOMObservation.cpp
  View file @07fad2d
... ... @@ -1214,7 +1214,7 @@ void GCOMObservation::compute_drb_phinor(const GCOMDri&amp; drm)
1214 1214 // Get DRI sky maps
1215 1215 const GSkyMap& map_dre = cube->dre().map();
1216 1216 const GSkyMap& map_drm = drm.map();
1217   - const GSkyMap& map_drg = m_drg.map();
  1217 + GSkyMap map_drg = get_weighted_drg_map();
1218 1218 GSkyMap& map_drb = const_cast<GSkyMap&>(m_drb.map());
1219 1219  
1220 1220 // Get data space dimensions
... ... @@ -1298,7 +1298,7 @@ void GCOMObservation::compute_drb_bgdlixa(const GCOMDri&amp; drm,
1298 1298 // Get references to relevant sky maps
1299 1299 const GSkyMap& map_dre = cube->dre().map();
1300 1300 const GSkyMap& map_drm = drm.map();
1301   - const GSkyMap& map_drg = m_drg.map();
  1301 + GSkyMap map_drg = get_weighted_drg_map();
1302 1302 GSkyMap& map_drb = const_cast<GSkyMap&>(m_drb.map());
1303 1303 GSkyMap& map_dri = const_cast<GSkyMap&>(dri.map());
1304 1304  
... ... @@ -1573,6 +1573,36 @@ void GCOMObservation::compute_drb_bgdlixa(const GCOMDri&amp; drm,
1573 1573  
1574 1574  
1575 1575 /***********************************************************************//**
  1576 + * @brief Return weighted DRG map
  1577 + *
  1578 + * @return Weighted DRG map.
  1579 + *
  1580 + * Returns a DRG as sky map where each pixel of the map is multiplied by the
  1581 + * solidangle of the pixel.
  1582 + ***************************************************************************/
  1583 +GSkyMap GCOMObservation::get_weighted_drg_map(void) const
  1584 +{
  1585 + // Initialise
  1586 + GSkyMap drg = m_drg.map();
  1587 +
  1588 + // Get data space dimensions
  1589 + int npix = drg.npix();
  1590 + int nphibar = drg.nmaps();
  1591 +
  1592 + // Weight map
  1593 + for (int ipix = 0; ipix < npix; ++ipix) {
  1594 + double omega = drg.solidangle(ipix);
  1595 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1596 + drg(ipix,iphibar) *= omega;
  1597 + }
  1598 + }
  1599 +
  1600 + // Return weighted DRG map
  1601 + return drg;
  1602 +}
  1603 +
  1604 +
  1605 +/***********************************************************************//**
1576 1606 * @brief Compute Phibar index range for BGDLIXA background method
1577 1607 *
1578 1608 * @param[in] iphibar Phibar layer index.
... ...
inst/com/src/GCOMResponse.cpp
  View file @07fad2d
... ... @@ -242,7 +242,7 @@ GCOMResponse* GCOMResponse::clone(void) const
242 242 * @param[in] event Observed event.
243 243 * @param[in] photon Incident photon.
244 244 * @param[in] obs Observation.
245   - * @return Instrument response function (cm^2 sr^-1)
  245 + * @return Instrument response function (\f$cm^2 sr^{-1}\f$)
246 246 *
247 247 * @exception GException::invalid_argument
248 248 * Observation is not a COMPTEL observation.
... ... @@ -250,22 +250,29 @@ GCOMResponse* GCOMResponse::clone(void) const
250 250 * @exception GException::invalid_value
251 251 * Response not initialised with a valid IAQ
252 252 *
253   - *
254 253 * Returns the instrument response function for a given observed photon
255 254 * direction as function of the assumed true photon direction. The result
256 255 * is given by
257   - * \f[IRF = \frac{IAQ \times DRG \times DRX}{ontime \times ewidth}\f]
  256 + *
  257 + * \f[
  258 + * {\tt IRF} = \frac{{\tt IAQ} \times {\tt DRG} \times {\tt DRX}}
  259 + * {T \times {\tt TOFCOR}}
  260 + * \f]
  261 + *
258 262 * where
259   - * \f$IRF\f$ is the instrument response function,
260   - * \f$IAQ\f$ is the COMPTEL response matrix (sr-1),
261   - * \f$DRG\f$ is the geometry factor (cm2),
262   - * \f$DRX\f$ is the exposure (s),
263   - * \f$ontime\f$ is the ontime (s), and
264   - * \f$ewidth\f$ is the energy width (MeV).
  263 + * - \f${\tt IRF}\f$ is the instrument response function (\f$cm^2 sr^{-1}\f$),
  264 + * - \f${\tt IAQ}\f$ is the COMPTEL response matrix (\f$sr^{-1}\f$),
  265 + * - \f${\tt DRG}\f$ is the geometry factor (probability),
  266 + * - \f${\tt DRX}\f$ is the exposure (\f$cm^2 s\f$),
  267 + * - \f$T\f$ is the ontime (\f$s\f$), and
  268 + * - \f${\tt TOFCOR}\f$ is a correction that accounts for the Time of Flight
  269 + * selection window.
265 270 *
266   - * The observed photon direction is spanned by the 3 values (Chi,Psi,Phibar).
267   - * (Chi,Psi) is the scatter direction of the event, given in sky coordinates.
268   - * Phibar is the Compton scatter angle, computed from the energy deposits.
  271 + * The observed photon direction is spanned by the three values \f$\Chi\f$,
  272 + * \f$\Psi\f$, and \f$\bar{\varphi})\f$. \f$\Chi\f$ and \f$\Psi\f$ is the
  273 + * scatter direction of the event, given in sky coordinates.
  274 + * \f$\bar{\varphi}\f$ is the Compton scatter angle, computed from the
  275 + * energy deposits in the two detector planes.
269 276 ***************************************************************************/
270 277 double GCOMResponse::irf(const GEvent& event,
271 278 const GPhoton& photon,
... ... @@ -360,14 +367,14 @@ double GCOMResponse::irf(const GEvent&amp; event,
360 367 // Continue only if IAQ is positive
361 368 if (iaq > 0.0) {
362 369  
363   - // Get DRG value (units: probability)
  370 + // Get DRG value (unit: probability)
364 371 double drg = observation->drg().map()(obsDir.dir(), iphibar);
365 372  
366   - // Get DRX value (units: cm^2 sec)
  373 + // Get DRX value (unit: cm^2 sec)
367 374 double drx = observation->drx().map()(srcDir);
368 375  
369   - // Get ontime
370   - double ontime = observation->ontime(); // sec
  376 + // Get ontime (unit: s)
  377 + double ontime = observation->ontime();
371 378  
372 379 // Get ToF correction
373 380 double tofcor = cube->dre().tof_correction();
... ... @@ -530,6 +537,29 @@ void GCOMResponse::load(const std::string&amp; rspname)
530 537 *
531 538 * Read the COMPTEL response from IAQ FITS file and convert the IAQ values
532 539 * into a probability per steradian.
  540 + *
  541 + * The IAQ values are divided by the solid angle of a Phigeo bin which is
  542 + * given by
  543 + *
  544 + * \f{eqnarray*}{
  545 + * \Omega & = & 2 \pi \left[
  546 + * \left(
  547 + * 1 - \cos \left( \varphi_{\rm geo} +
  548 + * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
  549 + * \right) -
  550 + * \left(
  551 + * 1 - \cos \left( \varphi_{\rm geo} -
  552 + * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
  553 + * \right) \right] \\
  554 + * &=& 2 \pi \left[
  555 + * \cos \left( \varphi_{\rm geo} -
  556 + * \frac{1}{2} \Delta \varphi_{\rm geo} \right) -
  557 + * \cos \left( \varphi_{\rm geo} +
  558 + * \frac{1}{2} \Delta \varphi_{\rm geo} \right)
  559 + * \right] \\
  560 + * &=& 4 \pi \sin \left( \varphi_{\rm geo} \right)
  561 + * \sin \left( \frac{1}{2} \Delta \varphi_{\rm geo} \right)
  562 + * \f}
533 563 ***************************************************************************/
534 564 void GCOMResponse::read(const GFitsImage& image)
535 565 {
... ... @@ -549,8 +579,10 @@ void GCOMResponse::read(const GFitsImage&amp; image)
549 579 m_phibar_bin_size = image.real("CDELT2");
550 580  
551 581 // Get axes minima (values of first bin)
552   - m_phigeo_min = m_phigeo_ref_value + (1.0-m_phigeo_ref_pixel) * m_phigeo_bin_size;
553   - m_phibar_min = m_phibar_ref_value + (1.0-m_phibar_ref_pixel) * m_phibar_bin_size;
  582 + m_phigeo_min = m_phigeo_ref_value + (1.0-m_phigeo_ref_pixel) *
  583 + m_phigeo_bin_size;
  584 + m_phibar_min = m_phibar_ref_value + (1.0-m_phibar_ref_pixel) *
  585 + m_phibar_bin_size;
554 586  
555 587 // Compute IAQ size. Continue only if size is positive
556 588 int size = m_phigeo_bins * m_phibar_bins;
... ... @@ -566,13 +598,16 @@ void GCOMResponse::read(const GFitsImage&amp; image)
566 598  
567 599 } // endif: size was positive
568 600  
  601 + // Precompute variable for conversion
  602 + double phigeo_min = m_phigeo_min * gammalib::deg2rad;
  603 + double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
  604 + double omega0 = gammalib::fourpi * std::sin(0.5 * phigeo_bin);
  605 +
569 606 // Convert IAQ matrix from probability per Phigeo bin into a
570 607 // probability per steradian
571   - double omega0 = gammalib::fourpi *
572   - std::sin(0.5 * m_phigeo_bin_size * gammalib::deg2rad);
573 608 for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
574   - double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
575   - double omega = omega0 * std::sin(phigeo * gammalib::deg2rad);
  609 + double phigeo = iphigeo * phigeo_bin + phigeo_min;
  610 + double omega = omega0 * std::sin(phigeo);
576 611 for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
577 612 m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
578 613 }
... ... @@ -774,7 +809,7 @@ void GCOMResponse::free_members(void)
774 809 * @param[in] obs Observation.
775 810 * @param[out] gradients Gradients matrix.
776 811 * @return Instrument response to point source for all events in
777   - * observation.
  812 + * observation (\f$cm^2\f$).
778 813 *
779 814 * @exception GException::invalid_argument
780 815 * Observation is not a COMPTEL observation.
... ... @@ -935,7 +970,7 @@ GVector GCOMResponse::irf_ptsrc(const GModelSky&amp; model,
935 970 * @param[in] obs Observation.
936 971 * @param[out] gradients Gradients matrix.
937 972 * @return Instrument response to radial source for all events in
938   - * observation.
  973 + * observation (\f$cm^2\f$).
939 974 *
940 975 * @todo Implement method.
941 976 ***************************************************************************/
... ... @@ -967,7 +1002,7 @@ GVector GCOMResponse::irf_radial(const GModelSky&amp; model,
967 1002 * @param[in] obs Observation.
968 1003 * @param[out] gradients Gradients matrix.
969 1004 * @return Instrument response to elliptical source for all events in
970   - * observation.
  1005 + * observation (\f$cm^2\f$).
971 1006 *
972 1007 * @todo Implement method.
973 1008 ***************************************************************************/
... ... @@ -999,7 +1034,7 @@ GVector GCOMResponse::irf_elliptical(const GModelSky&amp; model,
999 1034 * @param[in] obs Observation.
1000 1035 * @param[out] gradients Gradients matrix.
1001 1036 * @return Instrument response to diffuse source for all events in
1002   - * observation.
  1037 + * observation (\f$cm^2\f$).
1003 1038 *
1004 1039 * @exception GException::invalid_argument
1005 1040 * Observation is not a COMPTEL observation.
... ... @@ -1067,8 +1102,9 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1067 1102 GVector irfs(nevents);
1068 1103  
1069 1104 // Initialise some variables
1070   - double phigeo_min = m_phigeo_min * gammalib::deg2rad;
  1105 + double phigeo_min = m_phigeo_min * gammalib::deg2rad;
1071 1106 double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
  1107 + double omega0 = 2.0 * std::sin(0.5 * phigeo_bin);
1072 1108 const GSkyMap& drx = obs_ptr->drx().map();
1073 1109 const double* drg = obs_ptr->drg().map().pixels();
1074 1110  
... ... @@ -1092,18 +1128,20 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1092 1128 for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
1093 1129  
1094 1130 // Determine Phigeo value in radians
1095   - double phigeo = phigeo_min + iphigeo * phigeo_bin;
  1131 + double phigeo = phigeo_min + iphigeo * phigeo_bin;
  1132 + double sin_phigeo = std::sin(phigeo);
1096 1133  
1097 1134 // Determine number of azimuthal integration steps and step size
1098 1135 // in radians
1099   - int naz = int(gammalib::twopi * phigeo / phigeo_bin + 0.5);
  1136 + double length = gammalib::twopi * sin_phigeo;
  1137 + int naz = int(length / phigeo_bin + 0.5);
1100 1138 if (naz < 2) {
1101 1139 naz = 2;
1102 1140 }
1103 1141 double az_step = gammalib::twopi / double(naz);
1104 1142  
1105 1143 // Computes solid angle of integration bin multiplied by Jaccobian
1106   - double omega = phigeo_bin * az_step * std::sin(phigeo);
  1144 + double omega = omega0 * az_step * sin_phigeo;
1107 1145  
1108 1146 // Loop over azimuth angle
1109 1147 double az = 0.0;
... ... @@ -1121,7 +1159,7 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1121 1159 // Set photon
1122 1160 GPhoton photon(skyDir, bin->energy(), bin->time());
1123 1161  
1124   - // Get model sky intensity for photon
  1162 + // Get model sky intensity for photon (unit: sr^-1)
1125 1163 double intensity = model.spatial()->eval(photon);
1126 1164  
1127 1165 // Fall through if intensity is zero
... ... @@ -1129,10 +1167,10 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1129 1167 continue;
1130 1168 }
1131 1169  
1132   - // Multiply intensity by DRX value
  1170 + // Multiply intensity by DRX value (unit: cm^2 s sr^-1)
1133 1171 intensity *= drx(skyDir);
1134 1172  
1135   - // Multiply intensity by solid angle
  1173 + // Multiply intensity by solid angle (unit: cm^2 s)
1136 1174 intensity *= omega;
1137 1175  
1138 1176 // Loop over Phibar
... ... @@ -1152,7 +1190,7 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1152 1190 // Get DRI index
1153 1191 int idri = ipix + iphibar * npix;
1154 1192  
1155   - // Compute IRF value
  1193 + // Compute IRF value (unit: cm^2)
1156 1194 double irf = iaq * drg[idri] * iaq_norm * intensity;
1157 1195  
1158 1196 // Add IRF value if it is positive
... ...