Luigi Tibaldo / gammalib

Commit 67f50dcfb69fd82183c0c5d321d762771647048a

Authored by Jürgen Knödlseder
1 parent 4dfdd7c5

Refactor COMPTEL response computation (#3889)

Inline Side-by-side
Showing 3 changed files with 246 additions and 742 deletions
inst/com/src/GCOMResponse.cpp
  View file @67f50dc
... ... @@ -73,16 +73,10 @@
73 73 " GObservation&, GMatrix*)"
74 74 #define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\
75 75 "GMatrix*)"
76   -#define G_IRF_EXTENDED "GCOMResponse::irf_extended(GModelSky&, "\
77   - "GObservation&, GSkyDir&, double&, GMatrix*)"
78 76  
79 77 /* __ Macros _____________________________________________________________ */
80 78  
81 79 /* __ Coding definitions _________________________________________________ */
82   -#define G_USE_DRM_CUBE
83   -#define G_USE_MODEL_INTEGRAL_IN_IRF_EXTENDED
84   -#define G_USE_INTEGRAL_IN_IRF_EXTENDED
85   -//#define G_APPROXIMATE_DRX_IN_IRF_EXTENDED
86 80  
87 81 /* __ Debug definitions __________________________________________________ */
88 82  
... ... @@ -991,90 +985,10 @@ GVector GCOMResponse::irf_radial(const GModelSky& model,
991 985 const GObservation& obs,
992 986 GMatrix* gradients) const
993 987 {
994   - // Get pointer to radial model
995   - const GModelSpatialRadial* radial =
996   - static_cast<const GModelSpatialRadial*>(model.spatial());
997   -
998   - // Get references to model direction and maximum model radius. Reduce
999   - // the maximum model radius by a small amount to avoid rounding
1000   - // problems at the boundary of a sharp edged model.
1001   - const GSkyDir& model_dir = radial->dir();
1002   - const double& theta_max = radial->theta_max() - 1.0e-12;
1003   -
1004   - // Call support method
1005   - GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients);
1006   -
1007   - // Return IRF vector
1008   - return irfs;
1009   -}
1010   -
1011   -
1012   -/***********************************************************************//**
1013   - * @brief Return instrument response to elliptical source
1014   - *
1015   - * @param[in] model Sky model.
1016   - * @param[in] obs Observation.
1017   - * @param[out] gradients Gradients matrix.
1018   - * @return Instrument response to elliptical source for all events in
1019   - * observation (\f$cm^2\f$).
1020   - *
1021   - * Returns the instrument response to an elliptical source for all events in
1022   - * the observations. See GCOMResponse::irf_extended for more information.
1023   - ***************************************************************************/
1024   -GVector GCOMResponse::irf_elliptical(const GModelSky& model,
1025   - const GObservation& obs,
1026   - GMatrix* gradients) const
1027   -{
1028   - // Get pointer to elliptical model
1029   - const GModelSpatialElliptical* elliptical =
1030   - static_cast<const GModelSpatialElliptical*>(model.spatial());
1031   -
1032   - // Get references to model direction and maximum model radius
1033   - const GSkyDir& model_dir = elliptical->dir();
1034   - const double& theta_max = elliptical->theta_max();
1035   -
1036   - // Call support method
1037   - GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients);
1038   -
1039   - // Return IRF value
1040   - return irfs;
1041   -}
1042   -
  988 + // Set constants
  989 + static const int iter_rho = 5;
  990 + static const int iter_omega = 5;
1043 991  
1044   -/***********************************************************************//**
1045   - * @brief Return instrument response to diffuse source
1046   - *
1047   - * @param[in] model Sky model.
1048   - * @param[in] obs Observation.
1049   - * @param[out] gradients Gradients matrix.
1050   - * @return Instrument response to diffuse source for all events in
1051   - * observation (\f$cm^2\f$).
1052   - *
1053   - * @exception GException::invalid_argument
1054   - * Observation is not a COMPTEL observation.
1055   - * @exception GException::invalid_value
1056   - * Response not initialised with a valid IAQ
1057   - *
1058   - * Returns the instrument response to a diffuse source for all events in
1059   - * the observations. The diffuse source may be energy dependent.
1060   - *
1061   - * The computation is done by integrating the diffuse model for each pixel
1062   - * in Chi and Psi over a circular region centred on the Chi/Psi pixel with
1063   - * a radius equal to the maximum Phigeo value. The radial integration is
1064   - * done by looping over all Phigeo bins of the response. For each Phigeo
1065   - * value, the azimuthal integration is done by stepping with an angular step
1066   - * size that corresponds to the Phigeo step size (which typically is 1 deg).
1067   - *
1068   - * @p gradients is an optional matrix where the number of rows corresponds
1069   - * to the number of events in the observation and the number of columns
1070   - * corresponds to the number of spatial model parameters. Since for point
1071   - * sources no gradients are computed, the method does not alter the
1072   - * content of @p gradients.
1073   - ***************************************************************************/
1074   -GVector GCOMResponse::irf_diffuse(const GModelSky& model,
1075   - const GObservation& obs,
1076   - GMatrix* gradients) const
1077   -{
1078 992 // Extract COMPTEL observation
1079 993 const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
1080 994 if (obs_ptr == NULL) {
... ... @@ -1082,7 +996,7 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1082 996 std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1083 997 "observations. Please specify a COMPTEL observation "
1084 998 "as argument.";
1085   - throw GException::invalid_argument(G_IRF_DIFFUSE, msg);
  999 + throw GException::invalid_argument(G_IRF_RADIAL, msg);
1086 1000 }
1087 1001  
1088 1002 // Extract COMPTEL event cube
... ... @@ -1091,22 +1005,32 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1091 1005 std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1092 1006 "not contain a COMPTEL event cube. Please specify "
1093 1007 "a COMPTEL observation containing and event cube.";
1094   - throw GException::invalid_argument(G_IRF_DIFFUSE, msg);
  1008 + throw GException::invalid_argument(G_IRF_RADIAL, msg);
1095 1009 }
1096 1010  
1097 1011 // Throw an exception if COMPTEL response is not set or if
1098 1012 if (m_iaq.empty()) {
1099 1013 std::string msg = "COMPTEL response is empty. Please initialise the "
1100 1014 "response with an \"IAQ\".";
1101   - throw GException::invalid_value(G_IRF_DIFFUSE, msg);
  1015 + throw GException::invalid_value(G_IRF_RADIAL, msg);
1102 1016 }
1103 1017 else if (m_phigeo_bin_size == 0.0) {
1104 1018 std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1105 1019 "Please initialise the response with a valid "
1106 1020 "\"IAQ\".";
1107   - throw GException::invalid_value(G_IRF_DIFFUSE, msg);
  1021 + throw GException::invalid_value(G_IRF_RADIAL, msg);
1108 1022 }
1109 1023  
  1024 + // Get pointer to radial model
  1025 + const GModelSpatialRadial* radial =
  1026 + static_cast<const GModelSpatialRadial*>(model.spatial());
  1027 +
  1028 + // Get references to model direction and maximum model radius. Reduce
  1029 + // the maximum model radius by a small amount to avoid rounding
  1030 + // problems at the boundary of a sharp edged model.
  1031 + const GSkyDir& model_dir = radial->dir();
  1032 + const double& theta_max = radial->theta_max() - 1.0e-12;
  1033 +
1110 1034 // Get number of Chi/Psi pixels, Phibar layers and event bins
1111 1035 int npix = cube->naxis(0) * cube->naxis(1);
1112 1036 int nphibar = cube->naxis(2);
... ... @@ -1118,13 +1042,19 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1118 1042 // Initialise some variables
1119 1043 double phigeo_min = m_phigeo_min * gammalib::deg2rad;
1120 1044 double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
1121   - double omega0 = 2.0 * std::sin(0.5 * phigeo_bin);
  1045 + double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
1122 1046 const GSkyMap& drx = obs_ptr->drx().map();
1123 1047 const double* drg = obs_ptr->drg().map().pixels();
1124 1048  
1125 1049 // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
1126   - double iaq_norm = obs_ptr->deadc() /
1127   - (obs_ptr->ontime() * cube->dre().tof_correction());
  1050 + double iaq_norm = obs_ptr->deadc() / (obs_ptr->ontime() * cube->dre().tof_correction());
  1051 +
  1052 + // Setup rotation matrix for conversion towards sky coordinates
  1053 + GMatrix ry;
  1054 + GMatrix rz;
  1055 + ry.eulery(model_dir.dec_deg() - 90.0);
  1056 + rz.eulerz(-model_dir.ra_deg());
  1057 + GMatrix rot = (ry * rz).transpose();
1128 1058  
1129 1059 // Loop over Chi and Psi
1130 1060 for (int ipix = 0; ipix < npix; ++ipix) {
... ... @@ -1138,85 +1068,52 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1138 1068 // Get reference to Phigeo sky direction
1139 1069 const GSkyDir& phigeoDir = obsDir.dir();
1140 1070  
1141   - // Loop over Phigeo
1142   - for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
1143   -
1144   - // Determine Phigeo value in radians
1145   - double phigeo = phigeo_min + iphigeo * phigeo_bin;
1146   - double sin_phigeo = std::sin(phigeo);
1147   -
1148   - // Determine number of azimuthal integration steps and step size
1149   - // in radians
1150   - double length = gammalib::twopi * sin_phigeo;
1151   - int naz = int(length / phigeo_bin + 0.5);
1152   - if (naz < 2) {
1153   - naz = 2;
1154   - }
1155   - double az_step = gammalib::twopi / double(naz);
1156   -
1157   - // Computes solid angle of integration bin multiplied by Jaccobian
1158   - double omega = omega0 * az_step * sin_phigeo;
1159   -
1160   - // Loop over azimuth angle
1161   - double az = 0.0;
1162   - for (int iaz = 0; iaz < naz; ++iaz, az += az_step) {
1163   -
1164   - // Get sky direction
1165   - GSkyDir skyDir = phigeoDir;
1166   - skyDir.rotate(az, phigeo);
1167   -
1168   - // Fall through if sky direction is not contained in DRX
1169   - if (!drx.contains(skyDir)) {
1170   - continue;
1171   - }
1172   -
1173   - // Set photon
1174   - GPhoton photon(skyDir, bin->energy(), bin->time());
1175   -
1176   - // Get model sky intensity for photon (unit: sr^-1)
1177   - double intensity = model.spatial()->eval(photon);
1178   -
1179   - // Fall through if intensity is zero
1180   - if (intensity == 0.0) {
1181   - continue;
1182   - }
1183   -
1184   - // Multiply intensity by DRX value (unit: cm^2 s sr^-1)
1185   - intensity *= drx(skyDir);
1186   -
1187   - // Multiply intensity by solid angle (unit: cm^2 s)
1188   - intensity *= omega;
1189   -
1190   - // Loop over Phibar
1191   - for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1192   -
1193   - // Get IAQ index
1194   - int i = iphibar * m_phigeo_bins + iphigeo;
  1071 + // Compute angle between model centre and Phigeo sky direction (radians)
  1072 + // and position angle of model with respect to Phigeo sky direction (radians)
  1073 + double zeta = phigeoDir.dist(model_dir);
1195 1074  
1196   - // Get IAQ value
1197   - double iaq = m_iaq[i];
  1075 + // Set radial model zenith angle range
  1076 + double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0;
  1077 + double rho_max = zeta + phigeo_max;
  1078 + if (rho_min > theta_max) {
  1079 + rho_min = theta_max;
  1080 + }
  1081 + if (rho_max > theta_max) {
  1082 + rho_max = theta_max;
  1083 + }
1198 1084  
1199   - // Fall through if IAQ is not positive
1200   - if (iaq <= 0.0) {
1201   - continue;
1202   - }
  1085 + // Perform model zenith angle integration if interval is valid
  1086 + if (rho_max > rho_min) {
1203 1087  
1204   - // Get DRI index
1205   - int idri = ipix + iphibar * npix;
  1088 + // Initialise IRF vector
  1089 + GVector irf(nphibar);
1206 1090  
1207   - // Compute IRF value (unit: cm^2)
1208   - double irf = iaq * drg[idri] * iaq_norm * intensity;
  1091 + // Setup integration kernel
  1092 + com_radial_kerns_rho integrands(m_iaq,
  1093 + *radial,
  1094 + irf,
  1095 + bin,
  1096 + rot,
  1097 + drx,
  1098 + phigeo_bin,
  1099 + m_phigeo_bins,
  1100 + nphibar,
  1101 + iter_omega);
1209 1102  
1210   - // Add IRF value if it is positive
1211   - if (irf > 0.0) {
1212   - irfs[idri] += irf;
1213   - }
  1103 + // Setup integrator
  1104 + GIntegrals integral(&integrands);
  1105 + integral.fixed_iter(iter_rho);
1214 1106  
1215   - } // endfor: looped over Phibar
  1107 + // Integrate over Phigeo
  1108 + irf = integral.romberg(rho_min, rho_max, iter_rho);
1216 1109  
1217   - } // endfor: looped over azimuth angle around locus of IAQ
  1110 + // Add IRF to result
  1111 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1112 + int idri = ipix + iphibar * npix;
  1113 + irfs[idri] = irf[iphibar] * drg[idri] * iaq_norm;
  1114 + }
1218 1115  
1219   - } // endfor: looped over Phigeo angles
  1116 + } // endif: Phigeo angle interval was valid
1220 1117  
1221 1118 } // endfor: looped over Chi and Psi pixels
1222 1119  
... ... @@ -1226,41 +1123,20 @@ GVector GCOMResponse::irf_diffuse(const GModelSky&amp; model,
1226 1123  
1227 1124  
1228 1125 /***********************************************************************//**
1229   - * @brief Return instrument response to extended source
  1126 + * @brief Return instrument response to elliptical source
1230 1127 *
1231 1128 * @param[in] model Sky model.
1232 1129 * @param[in] obs Observation.
1233   - * @param[in] model_dir Centre direction of extended model.
1234   - * @param[in] theta_max Maximum theta angle (radians).
1235 1130 * @param[out] gradients Gradients matrix.
1236   - * @return Instrument response to radial source for all events in
  1131 + * @return Instrument response to elliptical source for all events in
1237 1132 * observation (\f$cm^2\f$).
1238 1133 *
1239   - * Returns the instrument response to an extended source for all events in
1240   - * the observations.
1241   - *
1242   - * The computation is done by integrating the extended model for each pixel
1243   - * in Chi and Psi over a circular region centred on the Chi/Psi pixel with
1244   - * a radius equal to the maximum Phigeo value. The integration is done in
1245   - * a spherical system, spanned by Phigeo and azimuth angle. Only coordiates
1246   - * are considered that plausibly overlap with the extended model. The
1247   - * integration is done by numerical summation. The integration step is by
1248   - * default set to 1/10 of the model size, with the constraint that it never
1249   - * becomes larger than the IAQ Phigeo binsize. This should allow to capture
1250   - * all extended models with sufficient precision in an acceptable amount of
1251   - * computing time.
1252   - *
1253   - * @p gradients is an optional matrix where the number of rows corresponds
1254   - * to the number of events in the observation and the number of columns
1255   - * corresponds to the number of spatial model parameters. Since no gradients
1256   - * are computed in this method the content of @p gradients if not altered.
  1134 + * Returns the instrument response to an elliptical source for all events in
  1135 + * the observations. See GCOMResponse::irf_extended for more information.
1257 1136 ***************************************************************************/
1258   -#if defined(G_USE_MODEL_INTEGRAL_IN_IRF_EXTENDED)
1259   -GVector GCOMResponse::irf_extended(const GModelSky& model,
1260   - const GObservation& obs,
1261   - const GSkyDir& model_dir,
1262   - const double& theta_max,
1263   - GMatrix* gradients) const
  1137 +GVector GCOMResponse::irf_elliptical(const GModelSky& model,
  1138 + const GObservation& obs,
  1139 + GMatrix* gradients) const
1264 1140 {
1265 1141 // Set constants
1266 1142 static const int iter_rho = 5;
... ... @@ -1273,7 +1149,7 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1273 1149 std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1274 1150 "observations. Please specify a COMPTEL observation "
1275 1151 "as argument.";
1276   - throw GException::invalid_argument(G_IRF_EXTENDED, msg);
  1152 + throw GException::invalid_argument(G_IRF_ELLIPTICAL, msg);
1277 1153 }
1278 1154  
1279 1155 // Extract COMPTEL event cube
... ... @@ -1282,22 +1158,30 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1282 1158 std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1283 1159 "not contain a COMPTEL event cube. Please specify "
1284 1160 "a COMPTEL observation containing and event cube.";
1285   - throw GException::invalid_argument(G_IRF_EXTENDED, msg);
  1161 + throw GException::invalid_argument(G_IRF_ELLIPTICAL, msg);
1286 1162 }
1287 1163  
1288 1164 // Throw an exception if COMPTEL response is not set or if
1289 1165 if (m_iaq.empty()) {
1290 1166 std::string msg = "COMPTEL response is empty. Please initialise the "
1291 1167 "response with an \"IAQ\".";
1292   - throw GException::invalid_value(G_IRF_EXTENDED, msg);
  1168 + throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
1293 1169 }
1294 1170 else if (m_phigeo_bin_size == 0.0) {
1295 1171 std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1296 1172 "Please initialise the response with a valid "
1297 1173 "\"IAQ\".";
1298   - throw GException::invalid_value(G_IRF_EXTENDED, msg);
  1174 + throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
1299 1175 }
1300 1176  
  1177 + // Get pointer to elliptical model
  1178 + const GModelSpatialElliptical* elliptical =
  1179 + static_cast<const GModelSpatialElliptical*>(model.spatial());
  1180 +
  1181 + // Get references to model direction and maximum model radius
  1182 + const GSkyDir& model_dir = elliptical->dir();
  1183 + const double& theta_max = elliptical->theta_max();
  1184 +
1301 1185 // Get number of Chi/Psi pixels, Phibar layers and event bins
1302 1186 int npix = cube->naxis(0) * cube->naxis(1);
1303 1187 int nphibar = cube->naxis(2);
... ... @@ -1356,16 +1240,16 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1356 1240 GVector irf(nphibar);
1357 1241  
1358 1242 // Setup integration kernel
1359   - com_extended_kerns_rho integrands(m_iaq,
1360   - model,
1361   - irf,
1362   - bin,
1363   - rot,
1364   - drx,
1365   - phigeo_bin,
1366   - m_phigeo_bins,
1367   - nphibar,
1368   - iter_omega);
  1243 + com_elliptical_kerns_rho integrands(m_iaq,
  1244 + model,
  1245 + irf,
  1246 + bin,
  1247 + rot,
  1248 + drx,
  1249 + phigeo_bin,
  1250 + m_phigeo_bins,
  1251 + nphibar,
  1252 + iter_omega);
1369 1253  
1370 1254 // Setup integrator
1371 1255 GIntegrals integral(&integrands);
... ... @@ -1387,19 +1271,42 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1387 1271 // Return IRF vector
1388 1272 return irfs;
1389 1273 }
1390   -#else
1391   -GVector GCOMResponse::irf_extended(const GModelSky& model,
1392   - const GObservation& obs,
1393   - const GSkyDir& model_dir,
1394   - const double& theta_max,
1395   - GMatrix* gradients) const
1396   -{
1397   - // Set constants
1398   - //const double step_fraction = 0.1;
1399   - const double step_fraction = 0.05;
1400   - const int iter_phigeo = 5;
1401   - const int iter_phi = 5;
1402 1274  
  1275 +
  1276 +/***********************************************************************//**
  1277 + * @brief Return instrument response to diffuse source
  1278 + *
  1279 + * @param[in] model Sky model.
  1280 + * @param[in] obs Observation.
  1281 + * @param[out] gradients Gradients matrix.
  1282 + * @return Instrument response to diffuse source for all events in
  1283 + * observation (\f$cm^2\f$).
  1284 + *
  1285 + * @exception GException::invalid_argument
  1286 + * Observation is not a COMPTEL observation.
  1287 + * @exception GException::invalid_value
  1288 + * Response not initialised with a valid IAQ
  1289 + *
  1290 + * Returns the instrument response to a diffuse source for all events in
  1291 + * the observations. The diffuse source may be energy dependent.
  1292 + *
  1293 + * The computation is done by integrating the diffuse model for each pixel
  1294 + * in Chi and Psi over a circular region centred on the Chi/Psi pixel with
  1295 + * a radius equal to the maximum Phigeo value. The radial integration is
  1296 + * done by looping over all Phigeo bins of the response. For each Phigeo
  1297 + * value, the azimuthal integration is done by stepping with an angular step
  1298 + * size that corresponds to the Phigeo step size (which typically is 1 deg).
  1299 + *
  1300 + * @p gradients is an optional matrix where the number of rows corresponds
  1301 + * to the number of events in the observation and the number of columns
  1302 + * corresponds to the number of spatial model parameters. Since for point
  1303 + * sources no gradients are computed, the method does not alter the
  1304 + * content of @p gradients.
  1305 + ***************************************************************************/
  1306 +GVector GCOMResponse::irf_diffuse(const GModelSky& model,
  1307 + const GObservation& obs,
  1308 + GMatrix* gradients) const
  1309 +{
1403 1310 // Extract COMPTEL observation
1404 1311 const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
1405 1312 if (obs_ptr == NULL) {
... ... @@ -1407,7 +1314,7 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1407 1314 std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
1408 1315 "observations. Please specify a COMPTEL observation "
1409 1316 "as argument.";
1410   - throw GException::invalid_argument(G_IRF_EXTENDED, msg);
  1317 + throw GException::invalid_argument(G_IRF_DIFFUSE, msg);
1411 1318 }
1412 1319  
1413 1320 // Extract COMPTEL event cube
... ... @@ -1416,20 +1323,20 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1416 1323 std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
1417 1324 "not contain a COMPTEL event cube. Please specify "
1418 1325 "a COMPTEL observation containing and event cube.";
1419   - throw GException::invalid_argument(G_IRF_EXTENDED, msg);
  1326 + throw GException::invalid_argument(G_IRF_DIFFUSE, msg);
1420 1327 }
1421 1328  
1422 1329 // Throw an exception if COMPTEL response is not set or if
1423 1330 if (m_iaq.empty()) {
1424 1331 std::string msg = "COMPTEL response is empty. Please initialise the "
1425 1332 "response with an \"IAQ\".";
1426   - throw GException::invalid_value(G_IRF_EXTENDED, msg);
  1333 + throw GException::invalid_value(G_IRF_DIFFUSE, msg);
1427 1334 }
1428 1335 else if (m_phigeo_bin_size == 0.0) {
1429 1336 std::string msg = "COMPTEL response has a zero Phigeo bin size. "
1430 1337 "Please initialise the response with a valid "
1431 1338 "\"IAQ\".";
1432   - throw GException::invalid_value(G_IRF_EXTENDED, msg);
  1339 + throw GException::invalid_value(G_IRF_DIFFUSE, msg);
1433 1340 }
1434 1341  
1435 1342 // Get number of Chi/Psi pixels, Phibar layers and event bins
... ... @@ -1441,23 +1348,16 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1441 1348 GVector irfs(nevents);
1442 1349  
1443 1350 // Initialise some variables
1444   - double phigeo_min = m_phigeo_min * gammalib::deg2rad;
1445   - double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
1446   - double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
1447   - const GSkyMap* drx = &(obs_ptr->drx().map());
1448   - const double* drg = obs_ptr->drg().map().pixels();
1449   - double cos_theta_max = std::cos(theta_max);
  1351 + double phigeo_min = m_phigeo_min * gammalib::deg2rad;
  1352 + double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
  1353 + double omega0 = 2.0 * std::sin(0.5 * phigeo_bin);
  1354 + const GSkyMap& drx = obs_ptr->drx().map();
  1355 + const double* drg = obs_ptr->drg().map().pixels();
1450 1356  
1451 1357 // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
1452 1358 double iaq_norm = obs_ptr->deadc() /
1453 1359 (obs_ptr->ontime() * cube->dre().tof_correction());
1454 1360  
1455   - // Code option: multiply DRX at model direction into IAQ normalisation
1456   - #if defined(G_APPROXIMATE_DRX_IN_IRF_EXTENDED)
1457   - iaq_norm *= (*drx)(model_dir);
1458   - drx = NULL;
1459   - #endif
1460   -
1461 1361 // Loop over Chi and Psi
1462 1362 for (int ipix = 0; ipix < npix; ++ipix) {
1463 1363  
... ... @@ -1470,150 +1370,32 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1470 1370 // Get reference to Phigeo sky direction
1471 1371 const GSkyDir& phigeoDir = obsDir.dir();
1472 1372  
1473   - // Compute angle between model centre and Phigeo sky direction (radians)
1474   - // and position angle of model with respect to Phigeo sky direction (radians)
1475   - double zeta = phigeoDir.dist(model_dir);
1476   - double phi0 = phigeoDir.posang(model_dir);
1477   -
1478   - // Setup rotation matrix for faster conversion towards sky coordinates
1479   - GMatrix ry;
1480   - GMatrix rz;
1481   - ry.eulery(phigeoDir.dec_deg() - 90.0);
1482   - rz.eulerz(-phigeoDir.ra_deg());
1483   - GMatrix rot = (ry * rz).transpose();
1484   -
1485   - // Setup azimuth angle range computation
1486   - double sin_zeta = std::sin(zeta);
1487   - double cos_zeta = std::cos(zeta);
1488   -
1489   - // Setup Phigeo integration range
1490   - double phigeo_low = zeta - theta_max;
1491   - double phigeo_up = zeta + theta_max;
1492   - if (phigeo_low < phigeo_min) {
1493   - phigeo_low = phigeo_min;
1494   - }
1495   - if (phigeo_up > phigeo_max) {
1496   - phigeo_up = phigeo_max;
1497   - }
1498   -
1499   - //
1500   - // Code option A: Use integral
1501   - // ===========================
1502   - #if defined(G_USE_INTEGRAL_IN_IRF_EXTENDED)
1503   - // Perform Phigeo angle integration if interval is valid
1504   - if (phigeo_up > phigeo_low) {
1505   -
1506   - // Initialise IRF vector
1507   - GVector irf(nphibar);
1508   -
1509   - // Setup integration kernel
1510   - com_extended_kerns_phigeo integrands(m_iaq,
1511   - model,
1512   - irf,
1513   - bin->energy(),
1514   - bin->time(),
1515   - rot,
1516   - drx,
1517   - phigeo_bin,
1518   - m_phigeo_bins,
1519   - nphibar,
1520   - ipix,
1521   - npix,
1522   - iaq_norm,
1523   - zeta,
1524   - phi0,
1525   - theta_max,
1526   - iter_phi);
1527   -
1528   - // Setup integrator
1529   - GIntegrals integral(&integrands);
1530   - integral.fixed_iter(iter_phigeo);
1531   -
1532   - // Integrate over Phigeo
1533   - irf = integral.romberg(phigeo_low, phigeo_up, iter_phigeo);
1534   -
1535   - // Add IRF to result
1536   - for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1537   - int idri = ipix + iphibar * npix;
1538   - irfs[idri] += irf[iphibar] * drg[idri];
1539   - }
1540   -
1541   - } // endif: Phigeo angle interval was valid
1542   -
1543   - //
1544   - // Code option B: Perform numerical integration
1545   - // ============================================
1546   - #else
1547   - // Setup Phigeo integration step-size and number of integration steps. By
1548   - // default the integration step size is one tenth of the model diameter, but
1549   - // not larger than the Phigeo binning of the IAQ. We also make sure that there
1550   - // are at least two Phigeo integration steps.
1551   - double phigeo_length = phigeo_up - phigeo_low;
1552   - double phigeo_step = step_fraction * phigeo_length;
1553   - if (phigeo_step > phigeo_bin) {
1554   - phigeo_step = phigeo_bin;
1555   - }
1556   - int nphigeo = int(phigeo_length/phigeo_step + 0.5);
1557   - if (nphigeo < 2) {
1558   - nphigeo = 2;
1559   - }
1560   - phigeo_step = phigeo_length / double(nphigeo);
1561   -
1562   - // Compute solid-angle for selected Phigeo integration step size
1563   - double omega0 = 2.0 * std::sin(0.5 * phigeo_step);
1564   -
1565   - // Integrate over Phigeo
1566   - double phigeo = phigeo_low;
1567   - for (int iphigeo = 0; iphigeo <= nphigeo; ++iphigeo, phigeo += phigeo_step) {
1568   -
1569   - // Compute half length of the arc (in radians) from a circle with
1570   - // radius phigeo that intersects with the model, defined as a circle
1571   - // with maximum radius m_theta_max
1572   - double dphi = 0.5 * gammalib::roi_arclength(phigeo,
1573   - zeta,
1574   - cos_zeta,
1575   - sin_zeta,
1576   - theta_max,
1577   - cos_theta_max);
1578   -
1579   - // If arc length is not positive then skip this Phigeo value
1580   - if (dphi <= 0.0) {
1581   - continue;
1582   - }
  1373 + // Loop over Phigeo
  1374 + for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
1583 1375  
1584   - // Compute sine and cosine of Phigeo
  1376 + // Determine Phigeo value in radians
  1377 + double phigeo = phigeo_min + iphigeo * phigeo_bin;
1585 1378 double sin_phigeo = std::sin(phigeo);
1586   - double cos_phigeo = std::cos(phigeo);
1587 1379  
1588 1380 // Determine number of azimuthal integration steps and step size
1589   - // in radians. By default the integration step size is one tenth of
1590   - // the arclength, but larger than the Phigeo binning of the IAQ. We
1591   - // also make sure that there are at least two integration steps.
1592   - double azimuth_length = 2.0 * dphi * sin_phigeo;
1593   - double az_step = step_fraction * azimuth_length;
1594   - if (az_step > phigeo_bin) {
1595   - az_step = phigeo_bin;
1596   - }
1597   - int naz = int(azimuth_length/az_step + 0.5);
  1381 + // in radians
  1382 + double length = gammalib::twopi * sin_phigeo;
  1383 + int naz = int(length / phigeo_bin + 0.5);
1598 1384 if (naz < 2) {
1599 1385 naz = 2;
1600 1386 }
1601   - az_step = 2.0 * dphi / double(naz);
  1387 + double az_step = gammalib::twopi / double(naz);
1602 1388  
1603 1389 // Computes solid angle of integration bin multiplied by Jaccobian
1604 1390 double omega = omega0 * az_step * sin_phigeo;
1605 1391  
1606 1392 // Loop over azimuth angle
1607   - double az = phi0 - dphi;
1608   - for (int iaz = 0; iaz <= naz; ++iaz, az += az_step) {
  1393 + double az = 0.0;
  1394 + for (int iaz = 0; iaz < naz; ++iaz, az += az_step) {
1609 1395  
1610 1396 // Get sky direction
1611   - double sin_az = std::sin(az);
1612   - double cos_az = std::cos(az);
1613   - GVector native(-cos_az*sin_phigeo, sin_az*sin_phigeo, cos_phigeo);
1614   - GVector dir = rot * native;
1615   - GSkyDir skyDir;
1616   - skyDir.celvector(dir);
  1397 + GSkyDir skyDir = phigeoDir;
  1398 + skyDir.rotate(az, phigeo);
1617 1399  
1618 1400 // Fall through if sky direction is not contained in DRX
1619 1401 if (!drx.contains(skyDir)) {
... ... @@ -1637,67 +1419,39 @@ GVector GCOMResponse::irf_extended(const GModelSky&amp; model,
1637 1419 // Multiply intensity by solid angle (unit: cm^2 s)
1638 1420 intensity *= omega;
1639 1421  
1640   - // Extract IAQ value by linear inter/extrapolation in Phigeo
1641   - double phirat = phigeo / phigeo_bin; // 0.5 at bin centre
1642   - int i_phigeo = int(phirat); // index into which Phigeo falls
1643   - double eps = phirat - i_phigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
1644   - double iaq = 0.0;
1645   -
1646   - // Continue only if Phigeo index is valid
1647   - if (i_phigeo < m_phigeo_bins) {
1648   -
1649   - // Loop over Phibar
1650   - for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
1651   -
1652   - // Get IAQ index
1653   - int i = iphibar * m_phigeo_bins + i_phigeo;
1654   -
1655   - // Get interpolated IAQ value
1656   - if (eps < 0.0) { // interpolate towards left
1657   - if (i_phigeo > 0) {
1658   - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
1659   - }
1660   - else {
1661   - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
1662   - }
1663   - }
1664   - else { // interpolate towards right
1665   - if (i_phigeo < m_phigeo_bins-1) {
1666   - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
1667   - }
1668   - else {
1669   - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
1670   - }
1671   - }
1672   -
1673   - // Fall through if IAQ is not positive
1674   - if (iaq <= 0.0) {
1675   - continue;
1676   - }
1677   -
1678   - // Get DRI index
1679   - int idri = ipix + iphibar * npix;
1680   -
1681   - // Compute IRF value (unit: cm^2)
1682   - double irf = iaq * drg[idri] * iaq_norm * intensity;
1683   -
1684   - // Add IRF value if it is positive
1685   - if (irf > 0.0) {
1686   - irfs[idri] += irf;
1687   - }
1688   -
1689   - } // endfor: looped over Phibar
1690   -
1691   - } // endif: Phigeo bin was valid
  1422 + // Loop over Phibar
  1423 + for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
  1424 +
  1425 + // Get IAQ index
  1426 + int i = iphibar * m_phigeo_bins + iphigeo;
  1427 +
  1428 + // Get IAQ value
  1429 + double iaq = m_iaq[i];
  1430 +
  1431 + // Fall through if IAQ is not positive
  1432 + if (iaq <= 0.0) {
  1433 + continue;
  1434 + }
  1435 +
  1436 + // Get DRI index
  1437 + int idri = ipix + iphibar * npix;
  1438 +
  1439 + // Compute IRF value (unit: cm^2)
  1440 + double irf = iaq * drg[idri] * iaq_norm * intensity;
  1441 +
  1442 + // Add IRF value if it is positive
  1443 + if (irf > 0.0) {
  1444 + irfs[idri] += irf;
  1445 + }
  1446 +
  1447 + } // endfor: looped over Phibar
1692 1448  
1693 1449 } // endfor: looped over azimuth angle around locus of IAQ
1694 1450  
1695 1451 } // endfor: looped over Phigeo angles
1696   - #endif
1697 1452  
1698 1453 } // endfor: looped over Chi and Psi pixels
1699 1454  
1700 1455 // Return IRF vector
1701 1456 return irfs;
1702 1457 }
1703   -#endif
... ...
inst/com/src/com_helpers_response_vector.cpp
  View file @67f50dc
... ... @@ -52,147 +52,6 @@
52 52 ==========================================================================*/
53 53  
54 54 /***********************************************************************//**
55   - * @brief Kernel for radial integration of extended model
56   - *
57   - * @param[in] phigeo Phigeo angle (radians).
58   - * @return Azimuthally integrated extended model.
59   - ***************************************************************************/
60   -GVector com_extended_kerns_phigeo::eval(const double& phigeo)
61   -{
62   - // Initialise kernel values
63   - m_irfs = 0.0;
64   -
65   - // Continue only if phigeo is positive
66   - if (phigeo > 0.0) {
67   -
68   - // Compute half length of the arc (in radians) from a circle with
69   - // radius phigeo that intersects with the model, defined as a
70   - // circle with maximum radius m_theta_max
71   - double dphi = 0.5 * gammalib::roi_arclength(phigeo,
72   - m_zeta,
73   - m_cos_zeta,
74   - m_sin_zeta,
75   - m_theta_max,
76   - m_cos_theta_max);
77   -
78   - // Continue only if arc length is positive
79   - if (dphi > 0.0) {
80   -
81   - // Compute sine and cosine of Phigeo
82   - double sin_phigeo = std::sin(phigeo);
83   - double cos_phigeo = std::cos(phigeo);
84   -
85   - // Compute phi integration range
86   - double phi_min = m_phi0 - dphi;
87   - double phi_max = m_phi0 + dphi;
88   -
89   - // Precompute interpolated IAQ vector for Phigeo angle
90   - double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
91   - int iphigeo = int(phirat); // index into which Phigeo falls
92   - double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
93   -
94   - // Continue only if Phigeo index is valid
95   - if (iphigeo < m_phigeo_bins) {
96   -
97   - // Setup kernel for azimuthal integration
98   - com_extended_kerns_phi integrand(m_model,
99   - m_srcEng,
100   - m_srcTime,
101   - m_rot,
102   - m_drx,
103   - sin_phigeo,
104   - cos_phigeo);
105   -
106   - // Setup integrator
107   - GIntegral integral(&integrand);
108   - integral.fixed_iter(m_iter);
109   -
110   - // Integrate over azimuth angle
111   - double irf = integral.romberg(phi_min, phi_max, m_iter) * sin_phigeo;
112   -
113   - // Multiply result with IAQ for all Phibar layers
114   - for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
115   -
116   - // Initialise IAQ
117   - double iaq;
118   -
119   - // Get IAQ index
120   - int i = iphibar * m_phigeo_bins + iphigeo;
121   -
122   - // Get interpolated IAQ value
123   - if (eps < 0.0) { // interpolate towards left
124   - if (iphigeo > 0) {
125   - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
126   - }
127   - else {
128   - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
129   - }
130   - }
131   - else { // interpolate towards right
132   - if (iphigeo < m_phigeo_bins-1) {
133   - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
134   - }
135   - else {
136   - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
137   - }
138   - }
139   -
140   - // Compute IRF
141   - m_irfs[iphibar] = irf * iaq * m_iaq_norm;
142   -
143   - } // endfor: looped over Phibar
144   -
145   - } // endif: Phigeo bin was valid
146   -
147   - } // endif: arc length was positive
148   -
149   - } // endif: phigeo was positive
150   -
151   - // Return kernel values
152   - return m_irfs;
153   -}
154   -
155   -
156   -/***********************************************************************//**
157   - * @brief Kernel for azimuthal integration of extended model
158   - *
159   - * @param[in] phi Azimuth angle (radians).
160   - * @return Vector of azimuthally integrated model.
161   - ***************************************************************************/
162   -double com_extended_kerns_phi::eval(const double& phi)
163   -{
164   - // Compute sine and cosine of azimuth angle
165   - double sin_phi = std::sin(phi);
166   - double cos_phi = std::cos(phi);
167   -
168   - // Get sky direction
169   - GVector native(-cos_phi*m_sin_phigeo, sin_phi*m_sin_phigeo, m_cos_phigeo);
170   - GVector dir = m_rot * native;
171   - GSkyDir skyDir;
172   - skyDir.celvector(dir);
173   -
174   - // Set photon
175   - GPhoton photon(skyDir, m_srcEng, m_srcTime);
176   -
177   - // Get model sky intensity for photon (unit: sr^-1)
178   - double intensity = m_model.spatial()->eval(photon);
179   -
180   - // Continue only if intensity is positive
181   - if (intensity > 0.0) {
182   -
183   - // Multiply-in DRX value if pointer is valid
184   - if (m_drx != NULL) {
185   - intensity *= (*m_drx)(skyDir);
186   - }
187   -
188   - } // endif: intensity was valid
189   -
190   - // Return kernel values
191   - return intensity;
192   -}
193   -
194   -
195   -/***********************************************************************//**
196 55 * @brief Kernel for radial integration of radial models
197 56 *
198 57 * @param[in] rho Rho angle (radians).
... ... @@ -320,12 +179,12 @@ GVector com_radial_kerns_omega::eval(const double&amp; omega)
320 179 return m_irfs;
321 180 }
322 181 /***********************************************************************//**
323   - * @brief Kernel for radial integration of extended model
  182 + * @brief Kernel for radial integration of elliptical models
324 183 *
325 184 * @param[in] rho Rho angle (radians).
326   - * @return Azimuthally integrated extended model.
  185 + * @return Azimuthally integrated elliptical model.
327 186 ***************************************************************************/
328   -GVector com_extended_kerns_rho::eval(const double& rho)
  187 +GVector com_elliptical_kerns_rho::eval(const double& rho)
329 188 {
330 189 // Initialise kernel values
331 190 m_irfs = 0.0;
... ... @@ -338,17 +197,17 @@ GVector com_extended_kerns_rho::eval(const double&amp; rho)
338 197 double cos_rho = std::cos(rho);
339 198  
340 199 // Setup azimuthal integration kernel
341   - com_extended_kerns_omega integrands(m_iaq,
342   - m_model,
343   - m_irfs,
344   - m_bin,
345   - m_rot,
346   - m_drx,
347   - m_phigeo_bin_size,
348   - m_phigeo_bins,
349   - m_phibar_bins,
350   - sin_rho,
351   - cos_rho);
  200 + com_elliptical_kerns_omega integrands(m_iaq,
  201 + m_model,
  202 + m_irfs,
  203 + m_bin,
  204 + m_rot,
  205 + m_drx,
  206 + m_phigeo_bin_size,
  207 + m_phigeo_bins,
  208 + m_phibar_bins,
  209 + sin_rho,
  210 + cos_rho);
352 211  
353 212 // Setup integrator
354 213 GIntegrals integral(&integrands);
... ... @@ -365,12 +224,12 @@ GVector com_extended_kerns_rho::eval(const double&amp; rho)
365 224  
366 225  
367 226 /***********************************************************************//**
368   - * @brief Kernel for azimuthal integration of extended model
  227 + * @brief Kernel for azimuthal integration of elliptical models
369 228 *
370 229 * @param[in] omega Omega angle (radians).
371   - * @return Kernel value for extended model.
  230 + * @return Kernel value for elliptical model.
372 231 ***************************************************************************/
373   -GVector com_extended_kerns_omega::eval(const double& omega)
  232 +GVector com_elliptical_kerns_omega::eval(const double& omega)
374 233 {
375 234 // Initialise kernel values
376 235 m_irfs = 0.0;
... ...
inst/com/src/com_helpers_response_vector.hpp
  View file @67f50dc
... ... @@ -44,116 +44,7 @@
44 44  
45 45  
46 46 /***********************************************************************//**
47   - * @class com_extended_kerns_phigeo
48   - *
49   - * @brief Kernel for Phigeo angle integration of extended model
50   - *
51   - * This class provides the kernel for the Phigeo angle integration of
52   - * extended sky models.
53   - ***************************************************************************/
54   -class com_extended_kerns_phigeo : public GFunctions {
55   -public:
56   - com_extended_kerns_phigeo(const std::vector<double>& iaq,
57   - const GModelSky& model,
58   - GVector& irfs,
59   - const GEnergy& srcEng,
60   - const GTime& srcTime,
61   - const GMatrix& rot,
62   - const GSkyMap* drx,
63   - const double& phigeo_bin_size,
64   - const int& phigeo_bins,
65   - const int& phibar_bins,
66   - const int& ipix,
67   - const int& npix,
68   - const double& iaq_norm,
69   - const double& zeta,
70   - const double& phi0,
71   - const double& theta_max,
72   - const int& iter) :
73   - m_iaq(iaq),
74   - m_model(model),
75   - m_irfs(irfs),
76   - m_srcEng(srcEng),
77   - m_srcTime(srcTime),
78   - m_rot(rot),
79   - m_drx(drx),
80   - m_phigeo_bin_size(phigeo_bin_size),
81   - m_phigeo_bins(phigeo_bins),
82   - m_phibar_bins(phibar_bins),
83   - m_ipix(ipix),
84   - m_npix(npix),
85   - m_iaq_norm(iaq_norm),
86   - m_zeta(zeta),
87   - m_sin_zeta(std::sin(zeta)),
88   - m_cos_zeta(std::cos(zeta)),
89   - m_phi0(phi0),
90   - m_theta_max(theta_max),
91   - m_cos_theta_max(std::cos(theta_max)),
92   - m_iter(iter) { }
93   - int size(void) const { return m_irfs.size(); }
94   - GVector eval(const double& phigeo);
95   -protected:
96   - const std::vector<double>& m_iaq; //!< IAQ vector
97   - const GModelSky& m_model; //!< Sky model
98   - GVector& m_irfs; //!< IRF vector to update
99   - const GEnergy& m_srcEng; //!< Source energy
100   - const GTime& m_srcTime; //!< Source time
101   - const GMatrix& m_rot; //!< Rotation matrix
102   - const GSkyMap* m_drx; //!< DRX
103   - const double& m_phigeo_bin_size; //!< Phigeo bin size
104   - const int& m_phigeo_bins; //!< Number of phigeo bins
105   - const int& m_phibar_bins; //!< Number of phibar bins
106   - const int& m_ipix; //!< DRI pixel index
107   - const int& m_npix; //!< Number of DRI pixels
108   - const double& m_iaq_norm; //!< IAQ normalisation
109   - const double& m_zeta; //!< Zeta angle
110   - double m_sin_zeta; //!< Sine of zeta angle
111   - double m_cos_zeta; //!< Cosine of zeta angle
112   - const double& m_phi0; //!< Model position angle
113   - const double& m_theta_max; //!< Maximum model radius
114   - double m_cos_theta_max; //!< Cosine of max. mod. rad.
115   - const int& m_iter; //!< Number of phi iterations
116   -};
117   -
118   -
119   -/***********************************************************************//**
120   - * @class com_extended_kerns_phi
121   - *
122   - * @brief Kernel for azimuth angle integration of extended model
123   - *
124   - * This class provides the kernel for the azimuth angle integration of
125   - * extended sky models.
126   - ***************************************************************************/
127   -class com_extended_kerns_phi : public GFunction {
128   -public:
129   - com_extended_kerns_phi(const GModelSky& model,
130   - const GEnergy& srcEng,
131   - const GTime& srcTime,
132   - const GMatrix& rot,
133   - const GSkyMap* drx,
134   - const double& sin_phigeo,
135   - const double& cos_phigeo) :
136   - m_model(model),
137   - m_srcEng(srcEng),
138   - m_srcTime(srcTime),
139   - m_rot(rot),
140   - m_drx(drx),
141   - m_sin_phigeo(sin_phigeo),
142   - m_cos_phigeo(cos_phigeo) { }
143   - double eval(const double& phi);
144   -protected:
145   - const GModelSky& m_model; //!< Sky model
146   - const GEnergy& m_srcEng; //!< Source energy
147   - const GTime& m_srcTime; //!< Source time
148   - const GMatrix& m_rot; //!< Rotation matrix
149   - const GSkyMap* m_drx; //!< DRX
150   - const double& m_sin_phigeo; //!< Sine of Phigeo
151   - const double& m_cos_phigeo; //!< Cosine of Phigeo
152   -};
153   -
154   -
155   -/***********************************************************************//**
156   - * @class com_extended_kerns_rho
  47 + * @class com_radial_kerns_rho
157 48 *
158 49 * @brief Kernel for rho angle integration of radial models
159 50 *
... ... @@ -199,12 +90,12 @@ protected:
199 90  
200 91  
201 92 /***********************************************************************//**
202   - * @class com_extended_kerns_omega
  93 + * @class com_radial_kerns_omega
203 94 *
204   - * @brief Kernel for azimuth angle integration of extended models
  95 + * @brief Kernel for azimuth angle integration of radial models
205 96 *
206 97 * This class provides the kernel for the azimuth angle integration of
207   - * extended sky models.
  98 + * radial sky models.
208 99 ***************************************************************************/
209 100 class com_radial_kerns_omega : public GFunctions {
210 101 public:
... ... @@ -245,35 +136,35 @@ protected:
245 136  
246 137  
247 138 /***********************************************************************//**
248   - * @class com_radial_kerns_rho
  139 + * @class com_elliptical_kerns_rho
249 140 *
250   - * @brief Kernel for rho angle integration of extended models
  141 + * @brief Kernel for rho angle integration of elliptical models
251 142 *
252   - * This class provides the kernel for the rho angle integration of extended
  143 + * This class provides the kernel for the rho angle integration of elliptical
253 144 * sky models.
254 145 ***************************************************************************/
255   -class com_extended_kerns_rho : public GFunctions {
  146 +class com_elliptical_kerns_rho : public GFunctions {
256 147 public:
257   - com_extended_kerns_rho(const std::vector<double>& iaq,
258   - const GModelSky& model,
259   - GVector& irfs,
260   - const GCOMEventBin* bin,
261   - const GMatrix& rot,
262   - const GSkyMap& drx,
263   - const double& phigeo_bin_size,
264   - const int& phigeo_bins,
265   - const int& phibar_bins,
266   - const int& iter) :
267   - m_iaq(iaq),
268   - m_model(model),
269   - m_irfs(irfs),
270   - m_bin(bin),
271   - m_rot(rot),
272   - m_drx(drx),
273   - m_phigeo_bin_size(phigeo_bin_size),
274   - m_phigeo_bins(phigeo_bins),
275   - m_phibar_bins(phibar_bins),
276   - m_iter(iter) { }
  148 + com_elliptical_kerns_rho(const std::vector<double>& iaq,
  149 + const GModelSky& model,
  150 + GVector& irfs,
  151 + const GCOMEventBin* bin,
  152 + const GMatrix& rot,
  153 + const GSkyMap& drx,
  154 + const double& phigeo_bin_size,
  155 + const int& phigeo_bins,
  156 + const int& phibar_bins,
  157 + const int& iter) :
  158 + m_iaq(iaq),
  159 + m_model(model),
  160 + m_irfs(irfs),
  161 + m_bin(bin),
  162 + m_rot(rot),
  163 + m_drx(drx),
  164 + m_phigeo_bin_size(phigeo_bin_size),
  165 + m_phigeo_bins(phigeo_bins),
  166 + m_phibar_bins(phibar_bins),
  167 + m_iter(iter) { }
277 168 int size(void) const { return m_irfs.size(); }
278 169 GVector eval(const double& phigeo);
279 170 protected:
... ... @@ -291,37 +182,37 @@ protected:
291 182  
292 183  
293 184 /***********************************************************************//**
294   - * @class com_radial_kerns_omega
  185 + * @class com_elliptical_kerns_omega
295 186 *
296   - * @brief Kernel for azimuth angle integration of extended models
  187 + * @brief Kernel for azimuth angle integration of elliptical models
297 188 *
298 189 * This class provides the kernel for the azimuth angle integration of
299   - * extended sky models.
  190 + * elliptical sky models.
300 191 ***************************************************************************/
301   -class com_extended_kerns_omega : public GFunctions {
  192 +class com_elliptical_kerns_omega : public GFunctions {
302 193 public:
303   - com_extended_kerns_omega(const std::vector<double>& iaq,
304   - const GModelSky& model,
305   - GVector& irfs,
306   - const GCOMEventBin* bin,
307   - const GMatrix& rot,
308   - const GSkyMap& drx,
309   - const double& phigeo_bin_size,
310   - const int& phigeo_bins,
311   - const int& phibar_bins,
312   - const double& sin_rho,
313   - const double& cos_rho) :
314   - m_iaq(iaq),
315   - m_model(model),
316   - m_irfs(irfs),
317   - m_bin(bin),
318   - m_rot(rot),
319   - m_drx(drx),
320   - m_phigeo_bin_size(phigeo_bin_size),
321   - m_phigeo_bins(phigeo_bins),
322   - m_phibar_bins(phibar_bins),
323   - m_sin_rho(sin_rho),
324   - m_cos_rho(cos_rho) { }
  194 + com_elliptical_kerns_omega(const std::vector<double>& iaq,
  195 + const GModelSky& model,
  196 + GVector& irfs,
  197 + const GCOMEventBin* bin,
  198 + const GMatrix& rot,
  199 + const GSkyMap& drx,
  200 + const double& phigeo_bin_size,
  201 + const int& phigeo_bins,
  202 + const int& phibar_bins,
  203 + const double& sin_rho,
  204 + const double& cos_rho) :
  205 + m_iaq(iaq),
  206 + m_model(model),
  207 + m_irfs(irfs),
  208 + m_bin(bin),
  209 + m_rot(rot),
  210 + m_drx(drx),
  211 + m_phigeo_bin_size(phigeo_bin_size),
  212 + m_phigeo_bins(phigeo_bins),
  213 + m_phibar_bins(phibar_bins),
  214 + m_sin_rho(sin_rho),
  215 + m_cos_rho(cos_rho) { }
325 216 int size(void) const { return m_irfs.size(); }
326 217 GVector eval(const double& omega);
327 218 protected:
... ...