Domenico Tiziani / gammalib

Commit 3b85d07c344dae38bcadaa1ae43d20d46574f0d4

Authored by knodlseder
1 parent 2b014284

add source region phi range selection

Inline Side-by-side
Showing 2 changed files with 168 additions and 71 deletions
inst/cta/include/GCTAResponse.hpp
  View file @3b85d07
... ... @@ -161,6 +161,7 @@ public:
161 161 // Analytical PSF implementation
162 162 double psf_dummy(const double& delta, const double& sigma) const;
163 163 double psf_dummy_sigma(const double& srcLogEng) const;
  164 + double psf_dummy_max(const double& sigma) const;
164 165  
165 166 private:
166 167 // Private methods
... ... @@ -178,6 +179,7 @@ private:
178 179 const GCTADir* obsCam,
179 180 const GCTADir* srcCam,
180 181 double delta_max,
  182 + double src_max,
181 183 double zenith,
182 184 double azimuth,
183 185 double srcLogEng,
... ... @@ -190,6 +192,8 @@ private:
190 192 m_src_cam(srcCam),
191 193 m_delta_max(delta_max),
192 194 m_cos_delta_max(std::cos(delta_max)),
  195 + m_src_max(src_max),
  196 + m_cos_src_max(std::cos(src_max)),
193 197 m_zenith(zenith),
194 198 m_azimuth(azimuth),
195 199 m_srcLogEng(srcLogEng),
... ... @@ -204,6 +208,8 @@ private:
204 208 const GCTADir* m_src_cam; //!< Pointer to source location in camera
205 209 double m_delta_max; //!< Maximum PSF angle delta_max
206 210 double m_cos_delta_max; //!< Cosine of delta_max
  211 + double m_src_max; //!< Maximum source radius
  212 + double m_cos_src_max; //!< Cosine of source radius
207 213 double m_zenith; //!< Telescope zenith
208 214 double m_azimuth; //!< Telescope azimuth
209 215 double m_srcLogEng; //!< True photon energy
... ...
inst/cta/src/GCTAResponse.cpp
  View file @3b85d07
... ... @@ -558,6 +558,9 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir&amp; obsDir,
558 558 * @todo The telescope zenith and azimuth angles as well as the radial offset
559 559 * and polar angles in the camera are not yet computed correctly (as
560 560 * the actual IRF implement does not need these values).
  561 + *
  562 + * @todo The theta limits could still improved as only the theta range is
  563 + * relevant for which the source and the PSF circle overlap.
561 564 ***************************************************************************/
562 565 double GCTAResponse::irf_extended(const GCTAInstDir& obsDir,
563 566 const GEnergy& obsEng,
... ... @@ -585,34 +588,54 @@ double GCTAResponse::irf_extended(const GCTAInstDir&amp; obsDir,
585 588 // Get PSF sigma
586 589 double sigma = psf_dummy_sigma(srcLogEng);
587 590  
588   - // Get maximum angular separation for which PSF is significant
589   - double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
590   -
591 591 // Determine camera direction of measured photon
592 592 GCTADir obsCam(obsDir, *pnt);
593 593  
594 594 // Determine camera direction of source model centre
595 595 GCTADir srcCam(model.dir(), *pnt);
596 596  
597   - // Setup integration kernel
598   - GCTAResponse::irf_kern_theta integrand(this,
599   - pnt,
600   - model.radial(),
601   - &obsCam,
602   - &srcCam,
603   - delta_max,
604   - zenith,
605   - azimuth,
606   - srcLogEng,
607   - obsLogEng,
608   - sigma);
609   -
610   - // Integrate over theta
611   - GIntegral integral(&integrand);
612   - integral.eps(m_eps);
613   - double theta_min = (obsCam.theta() > delta_max) ? obsCam.theta() - delta_max : 0.0;
614   - double theta_max = obsCam.theta() + delta_max;
615   - double irf = integral.romb(theta_min, theta_max);
  597 + // Get maximum angular separation for which PSF is significant (in radians).
  598 + // Set offset angle range that bounds PSF.
  599 + double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
  600 + double theta_psf_min = (obsCam.theta() > delta_max) ? obsCam.theta() - delta_max : 0.0;
  601 + double theta_psf_max = obsCam.theta() + delta_max;
  602 +
  603 + // Get maximum source radius in radians.
  604 + // Set offset angle range that bounds the source.
  605 + double src_max = model.radial()->theta_max();
  606 + double theta_src_min = (srcCam.theta() > src_max) ? srcCam.theta() - src_max : 0.0;
  607 + double theta_src_max = srcCam.theta() + src_max;
  608 +
  609 + // Set overlapping offset angle range for integration (in radians).
  610 + double theta_min = std::max(theta_psf_min, theta_src_min);
  611 + double theta_max = std::min(theta_psf_max, theta_src_max);
  612 +
  613 + // Initialise IRF value
  614 + double irf = 0.0;
  615 +
  616 + // Perform offset angle integration if interval is valid
  617 + if (theta_max > theta_min) {
  618 +
  619 + // Setup integration kernel
  620 + GCTAResponse::irf_kern_theta integrand(this,
  621 + pnt,
  622 + model.radial(),
  623 + &obsCam,
  624 + &srcCam,
  625 + delta_max,
  626 + src_max,
  627 + zenith,
  628 + azimuth,
  629 + srcLogEng,
  630 + obsLogEng,
  631 + sigma);
  632 +
  633 + // Integrate over theta
  634 + GIntegral integral(&integrand);
  635 + integral.eps(m_eps);
  636 + irf = integral.romb(theta_min, theta_max);
  637 +
  638 + }
616 639  
617 640 // Compile option: Show integration results
618 641 #if G_DEBUG_IRF_EXTENDED
... ... @@ -856,7 +879,7 @@ double GCTAResponse::psf_delta_max(const double&amp; theta,
856 879 double sigma = psf_dummy_sigma(srcLogEng);
857 880  
858 881 // Set maximum angular separation
859   - double delta_max = sigma * 5.0;
  882 + double delta_max = psf_dummy_max(sigma);
860 883  
861 884 // Return PSF
862 885 return delta_max;
... ... @@ -922,9 +945,7 @@ double GCTAResponse::npsf(const GSkyDir&amp; srcDir,
922 945 double radroi = roi.radius() * deg2rad;
923 946 double psf = roi.centre().dist(srcDir);
924 947 double sigma = psf_dummy_sigma(srcLogEng);
925   -
926   - // Get maximum PSF radius
927   - double rmax = 5.0*sigma;
  948 + double rmax = psf_dummy_max(sigma);
928 949  
929 950 // If PSF is sufficiently enclosed by ROI, skip the numerical integration
930 951 // and assume that the integral is 1.0
... ... @@ -1033,6 +1054,18 @@ double GCTAResponse::psf_dummy_sigma(const double&amp; srcLogEng) const
1033 1054 }
1034 1055  
1035 1056  
  1057 +/***********************************************************************//**
  1058 + * @brief Returns radius beyond which PSF is negligible (in radians)
  1059 + *
  1060 + * @param[in] sigma Width of point spread function (radians).
  1061 + ***************************************************************************/
  1062 +double GCTAResponse::psf_dummy_max(const double& sigma) const
  1063 +{
  1064 + // Return radius
  1065 + return (5.0*sigma);
  1066 +}
  1067 +
  1068 +
1036 1069 /*==========================================================================
1037 1070 = =
1038 1071 = Private methods =
... ... @@ -1042,8 +1075,20 @@ double GCTAResponse::psf_dummy_sigma(const double&amp; srcLogEng) const
1042 1075 /***********************************************************************//**
1043 1076 * @brief Initialise class members
1044 1077 *
1045   - * We set the relative integration precision to 1e-4 as test images have
1046   - * revealed smooth with this precision.
  1078 + * We set the relative integration precision to 1e-4 as test images are
  1079 + * pretty smooth with this precision and computations are reasonably fast.
  1080 + * We know that integration is not perfect, however, and the results are
  1081 + * still a little noisy. We could still improve a little more the current
  1082 + * theta integration scheme, but the mathematics have to be worked out
  1083 + * for this first.
  1084 + *
  1085 + * The following timing was obtained for a 0.4 deg disk model on a 100x100
  1086 + * grid with 10 energy layers on a 64 Bit machine (fermi):
  1087 + * m_eps = 1e-4 : user 0m12.894s
  1088 + * m_eps = 5e-5 : user 0m21.929s
  1089 + * m_eps = 1e-5 : user 1m53.607s
  1090 + * So the time increases about linearly with an increased precision
  1091 + * requirement.
1047 1092 ***************************************************************************/
1048 1093 void GCTAResponse::init_members(void)
1049 1094 {
... ... @@ -1108,7 +1153,7 @@ void GCTAResponse::free_members(void)
1108 1153 void GCTAResponse::read_performance_table(const std::string& filename)
1109 1154 {
1110 1155 // Allocate line buffer
1111   - const int n = 1000;
  1156 + const int n = 1000;
1112 1157 char line[n];
1113 1158  
1114 1159 // Open performance table readonly
... ... @@ -1165,53 +1210,99 @@ void GCTAResponse::read_performance_table(const std::string&amp; filename)
1165 1210 *
1166 1211 * This method provides the integration kernel for the radial offset angle
1167 1212 * IRF integration. The radial offset angle is given in the camera system.
  1213 + *
  1214 + * @todo I tried to cover the wrap around properly, yet I have not
  1215 + * demonstrated mathematically that it works in all cases (but I did
  1216 + * many tests to check that it works).
1168 1217 ***************************************************************************/
1169 1218 double GCTAResponse::irf_kern_theta::eval(double theta)
1170 1219 {
1171   - // Compute half interval length (in radians)
1172   - double delta_phi = 0.5 * cta_roi_arclength(theta,
1173   - m_obs_cam->theta(),
1174   - m_obs_cam->costheta(),
1175   - m_obs_cam->sintheta(),
1176   - m_delta_max, m_cos_delta_max);
1177   -
1178   - // Set phi interval (in radians)
1179   - double phi_min = m_obs_cam->phi() - delta_phi;
1180   - double phi_max = m_obs_cam->phi() + delta_phi;
1181   -
1182   - // Precompute cosine and sine terms for azimuthal integration
1183   - double cos_theta = std::cos(theta);
1184   - double sin_theta = std::sin(theta);
1185   - double cos_obs = cos_theta * m_obs_cam->costheta();
1186   - double sin_obs = sin_theta * m_obs_cam->sintheta();
1187   - double cos_src = cos_theta * m_src_cam->costheta();
1188   - double sin_src = sin_theta * m_src_cam->sintheta();
1189   -
1190   - // Setup integration kernel
1191   - GCTAResponse::irf_kern_phi integrand(m_rsp,
1192   - m_pnt,
1193   - m_radial,
1194   - m_obs_cam,
1195   - m_src_cam,
1196   - theta,
1197   - m_zenith,
1198   - m_azimuth,
1199   - m_srcLogEng,
1200   - m_obsLogEng,
1201   - m_sigma,
1202   - cos_obs,
1203   - sin_obs,
1204   - cos_src,
1205   - sin_src);
1206   -
1207   - // Integrate over phi
1208   - GIntegral integral(&integrand);
1209   - integral.eps(m_rsp->m_eps);
1210   - integral.silent(true);
1211   - double value = integral.romb(phi_min, phi_max) * sin_theta;
  1220 + // Compute half interval length for PSF circle (in radians)
  1221 + double dphi_psf = 0.5 * cta_roi_arclength(theta,
  1222 + m_obs_cam->theta(),
  1223 + m_obs_cam->costheta(),
  1224 + m_obs_cam->sintheta(),
  1225 + m_delta_max, m_cos_delta_max);
  1226 +
  1227 + // Compute half interval length for model circle (in radians)
  1228 + double dphi_src = 0.5 * cta_roi_arclength(theta,
  1229 + m_src_cam->theta(),
  1230 + m_src_cam->costheta(),
  1231 + m_src_cam->sintheta(),
  1232 + m_src_max, m_cos_src_max);
  1233 +
  1234 + // Set the phi phase to common half space
  1235 + double phi_psf = m_obs_cam->phi();
  1236 + double phi_src = m_src_cam->phi();
  1237 + if (phi_src-phi_psf > pi)
  1238 + phi_src -= twopi;
  1239 + if (phi_src-phi_psf < -pi)
  1240 + phi_src += twopi;
  1241 +
  1242 + // Set phi interval for PSF and source region (in radians)
  1243 + double phi_psf_min = phi_psf - dphi_psf;
  1244 + double phi_psf_max = phi_psf + dphi_psf;
  1245 + double phi_src_min = phi_src - dphi_src;
  1246 + double phi_src_max = phi_src + dphi_src;
  1247 +
  1248 + // Set common phi interval. If one of the intervals covers the entire
  1249 + // azimuth range then just apply the other interval
  1250 + double phi_min;
  1251 + double phi_max;
  1252 + if (dphi_src >= pi) {
  1253 + phi_min = phi_psf_min;
  1254 + phi_max = phi_psf_max;
  1255 + }
  1256 + else if (dphi_psf >= pi) {
  1257 + phi_min = phi_src_min;
  1258 + phi_max = phi_src_max;
  1259 + }
  1260 + else {
  1261 + phi_min = std::max(phi_psf_min, phi_src_min);
  1262 + phi_max = std::min(phi_psf_max, phi_src_max);
  1263 + }
  1264 +
  1265 + // Initialise result
  1266 + double irf = 0.0;
  1267 +
  1268 + // Perform phi integration (only if interval is valid)
  1269 + if (phi_min < phi_max) {
  1270 +
  1271 + // Precompute cosine and sine terms for azimuthal integration
  1272 + double cos_theta = std::cos(theta);
  1273 + double sin_theta = std::sin(theta);
  1274 + double cos_obs = cos_theta * m_obs_cam->costheta();
  1275 + double sin_obs = sin_theta * m_obs_cam->sintheta();
  1276 + double cos_src = cos_theta * m_src_cam->costheta();
  1277 + double sin_src = sin_theta * m_src_cam->sintheta();
  1278 +
  1279 + // Setup integration kernel
  1280 + GCTAResponse::irf_kern_phi integrand(m_rsp,
  1281 + m_pnt,
  1282 + m_radial,
  1283 + m_obs_cam,
  1284 + m_src_cam,
  1285 + theta,
  1286 + m_zenith,
  1287 + m_azimuth,
  1288 + m_srcLogEng,
  1289 + m_obsLogEng,
  1290 + m_sigma,
  1291 + cos_obs,
  1292 + sin_obs,
  1293 + cos_src,
  1294 + sin_src);
  1295 +
  1296 + // Integrate over phi
  1297 + GIntegral integral(&integrand);
  1298 + integral.eps(m_rsp->m_eps);
  1299 + integral.silent(true);
  1300 + irf = integral.romb(phi_min, phi_max) * sin_theta;
  1301 +
  1302 + } // endif: phi intervals was valid
1212 1303  
1213 1304 // Return result
1214   - return value;
  1305 + return irf;
1215 1306 }
1216 1307  
1217 1308  
... ...