diff --git a/inst/com/src/GCOMResponse.cpp b/inst/com/src/GCOMResponse.cpp index b23b838..342a7fd 100644 --- a/inst/com/src/GCOMResponse.cpp +++ b/inst/com/src/GCOMResponse.cpp @@ -73,16 +73,10 @@ " GObservation&, GMatrix*)" #define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\ "GMatrix*)" -#define G_IRF_EXTENDED "GCOMResponse::irf_extended(GModelSky&, "\ - "GObservation&, GSkyDir&, double&, GMatrix*)" /* __ Macros _____________________________________________________________ */ /* __ Coding definitions _________________________________________________ */ -#define G_USE_DRM_CUBE -#define G_USE_MODEL_INTEGRAL_IN_IRF_EXTENDED -#define G_USE_INTEGRAL_IN_IRF_EXTENDED -//#define G_APPROXIMATE_DRX_IN_IRF_EXTENDED /* __ Debug definitions __________________________________________________ */ @@ -991,90 +985,10 @@ GVector GCOMResponse::irf_radial(const GModelSky& model, const GObservation& obs, GMatrix* gradients) const { - // Get pointer to radial model - const GModelSpatialRadial* radial = - static_cast(model.spatial()); - - // Get references to model direction and maximum model radius. Reduce - // the maximum model radius by a small amount to avoid rounding - // problems at the boundary of a sharp edged model. - const GSkyDir& model_dir = radial->dir(); - const double& theta_max = radial->theta_max() - 1.0e-12; - - // Call support method - GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients); - - // Return IRF vector - return irfs; -} - - -/***********************************************************************//** - * @brief Return instrument response to elliptical source - * - * @param[in] model Sky model. - * @param[in] obs Observation. - * @param[out] gradients Gradients matrix. - * @return Instrument response to elliptical source for all events in - * observation (\f$cm^2\f$). - * - * Returns the instrument response to an elliptical source for all events in - * the observations. See GCOMResponse::irf_extended for more information. - ***************************************************************************/ -GVector GCOMResponse::irf_elliptical(const GModelSky& model, - const GObservation& obs, - GMatrix* gradients) const -{ - // Get pointer to elliptical model - const GModelSpatialElliptical* elliptical = - static_cast(model.spatial()); - - // Get references to model direction and maximum model radius - const GSkyDir& model_dir = elliptical->dir(); - const double& theta_max = elliptical->theta_max(); - - // Call support method - GVector irfs = irf_extended(model, obs, model_dir, theta_max, gradients); - - // Return IRF value - return irfs; -} - + // Set constants + static const int iter_rho = 5; + static const int iter_omega = 5; -/***********************************************************************//** - * @brief Return instrument response to diffuse source - * - * @param[in] model Sky model. - * @param[in] obs Observation. - * @param[out] gradients Gradients matrix. - * @return Instrument response to diffuse source for all events in - * observation (\f$cm^2\f$). - * - * @exception GException::invalid_argument - * Observation is not a COMPTEL observation. - * @exception GException::invalid_value - * Response not initialised with a valid IAQ - * - * Returns the instrument response to a diffuse source for all events in - * the observations. The diffuse source may be energy dependent. - * - * The computation is done by integrating the diffuse model for each pixel - * in Chi and Psi over a circular region centred on the Chi/Psi pixel with - * a radius equal to the maximum Phigeo value. The radial integration is - * done by looping over all Phigeo bins of the response. For each Phigeo - * value, the azimuthal integration is done by stepping with an angular step - * size that corresponds to the Phigeo step size (which typically is 1 deg). - * - * @p gradients is an optional matrix where the number of rows corresponds - * to the number of events in the observation and the number of columns - * corresponds to the number of spatial model parameters. Since for point - * sources no gradients are computed, the method does not alter the - * content of @p gradients. - ***************************************************************************/ -GVector GCOMResponse::irf_diffuse(const GModelSky& model, - const GObservation& obs, - GMatrix* gradients) const -{ // Extract COMPTEL observation const GCOMObservation* obs_ptr = dynamic_cast(&obs); if (obs_ptr == NULL) { @@ -1082,7 +996,7 @@ GVector GCOMResponse::irf_diffuse(const GModelSky& model, std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL " "observations. Please specify a COMPTEL observation " "as argument."; - throw GException::invalid_argument(G_IRF_DIFFUSE, msg); + throw GException::invalid_argument(G_IRF_RADIAL, msg); } // Extract COMPTEL event cube @@ -1091,22 +1005,32 @@ GVector GCOMResponse::irf_diffuse(const GModelSky& model, std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does " "not contain a COMPTEL event cube. Please specify " "a COMPTEL observation containing and event cube."; - throw GException::invalid_argument(G_IRF_DIFFUSE, msg); + throw GException::invalid_argument(G_IRF_RADIAL, msg); } // Throw an exception if COMPTEL response is not set or if if (m_iaq.empty()) { std::string msg = "COMPTEL response is empty. Please initialise the " "response with an \"IAQ\"."; - throw GException::invalid_value(G_IRF_DIFFUSE, msg); + throw GException::invalid_value(G_IRF_RADIAL, msg); } else if (m_phigeo_bin_size == 0.0) { std::string msg = "COMPTEL response has a zero Phigeo bin size. " "Please initialise the response with a valid " "\"IAQ\"."; - throw GException::invalid_value(G_IRF_DIFFUSE, msg); + throw GException::invalid_value(G_IRF_RADIAL, msg); } + // Get pointer to radial model + const GModelSpatialRadial* radial = + static_cast(model.spatial()); + + // Get references to model direction and maximum model radius. Reduce + // the maximum model radius by a small amount to avoid rounding + // problems at the boundary of a sharp edged model. + const GSkyDir& model_dir = radial->dir(); + const double& theta_max = radial->theta_max() - 1.0e-12; + // Get number of Chi/Psi pixels, Phibar layers and event bins int npix = cube->naxis(0) * cube->naxis(1); int nphibar = cube->naxis(2); @@ -1118,13 +1042,19 @@ GVector GCOMResponse::irf_diffuse(const GModelSky& model, // Initialise some variables double phigeo_min = m_phigeo_min * gammalib::deg2rad; double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad; - double omega0 = 2.0 * std::sin(0.5 * phigeo_bin); + double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins; const GSkyMap& drx = obs_ptr->drx().map(); const double* drg = obs_ptr->drg().map().pixels(); // Compute IAQ normalisation (1/s): DEADC / ONTIME (s) - double iaq_norm = obs_ptr->deadc() / - (obs_ptr->ontime() * cube->dre().tof_correction()); + double iaq_norm = obs_ptr->deadc() / (obs_ptr->ontime() * cube->dre().tof_correction()); + + // Setup rotation matrix for conversion towards sky coordinates + GMatrix ry; + GMatrix rz; + ry.eulery(model_dir.dec_deg() - 90.0); + rz.eulerz(-model_dir.ra_deg()); + GMatrix rot = (ry * rz).transpose(); // Loop over Chi and Psi for (int ipix = 0; ipix < npix; ++ipix) { @@ -1138,85 +1068,52 @@ GVector GCOMResponse::irf_diffuse(const GModelSky& model, // Get reference to Phigeo sky direction const GSkyDir& phigeoDir = obsDir.dir(); - // Loop over Phigeo - for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) { - - // Determine Phigeo value in radians - double phigeo = phigeo_min + iphigeo * phigeo_bin; - double sin_phigeo = std::sin(phigeo); - - // Determine number of azimuthal integration steps and step size - // in radians - double length = gammalib::twopi * sin_phigeo; - int naz = int(length / phigeo_bin + 0.5); - if (naz < 2) { - naz = 2; - } - double az_step = gammalib::twopi / double(naz); - - // Computes solid angle of integration bin multiplied by Jaccobian - double omega = omega0 * az_step * sin_phigeo; - - // Loop over azimuth angle - double az = 0.0; - for (int iaz = 0; iaz < naz; ++iaz, az += az_step) { - - // Get sky direction - GSkyDir skyDir = phigeoDir; - skyDir.rotate(az, phigeo); - - // Fall through if sky direction is not contained in DRX - if (!drx.contains(skyDir)) { - continue; - } - - // Set photon - GPhoton photon(skyDir, bin->energy(), bin->time()); - - // Get model sky intensity for photon (unit: sr^-1) - double intensity = model.spatial()->eval(photon); - - // Fall through if intensity is zero - if (intensity == 0.0) { - continue; - } - - // Multiply intensity by DRX value (unit: cm^2 s sr^-1) - intensity *= drx(skyDir); - - // Multiply intensity by solid angle (unit: cm^2 s) - intensity *= omega; - - // Loop over Phibar - for (int iphibar = 0; iphibar < nphibar; ++iphibar) { - - // Get IAQ index - int i = iphibar * m_phigeo_bins + iphigeo; + // Compute angle between model centre and Phigeo sky direction (radians) + // and position angle of model with respect to Phigeo sky direction (radians) + double zeta = phigeoDir.dist(model_dir); - // Get IAQ value - double iaq = m_iaq[i]; + // Set radial model zenith angle range + double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0; + double rho_max = zeta + phigeo_max; + if (rho_min > theta_max) { + rho_min = theta_max; + } + if (rho_max > theta_max) { + rho_max = theta_max; + } - // Fall through if IAQ is not positive - if (iaq <= 0.0) { - continue; - } + // Perform model zenith angle integration if interval is valid + if (rho_max > rho_min) { - // Get DRI index - int idri = ipix + iphibar * npix; + // Initialise IRF vector + GVector irf(nphibar); - // Compute IRF value (unit: cm^2) - double irf = iaq * drg[idri] * iaq_norm * intensity; + // Setup integration kernel + com_radial_kerns_rho integrands(m_iaq, + *radial, + irf, + bin, + rot, + drx, + phigeo_bin, + m_phigeo_bins, + nphibar, + iter_omega); - // Add IRF value if it is positive - if (irf > 0.0) { - irfs[idri] += irf; - } + // Setup integrator + GIntegrals integral(&integrands); + integral.fixed_iter(iter_rho); - } // endfor: looped over Phibar + // Integrate over Phigeo + irf = integral.romberg(rho_min, rho_max, iter_rho); - } // endfor: looped over azimuth angle around locus of IAQ + // Add IRF to result + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + int idri = ipix + iphibar * npix; + irfs[idri] = irf[iphibar] * drg[idri] * iaq_norm; + } - } // endfor: looped over Phigeo angles + } // endif: Phigeo angle interval was valid } // endfor: looped over Chi and Psi pixels @@ -1226,41 +1123,20 @@ GVector GCOMResponse::irf_diffuse(const GModelSky& model, /***********************************************************************//** - * @brief Return instrument response to extended source + * @brief Return instrument response to elliptical source * * @param[in] model Sky model. * @param[in] obs Observation. - * @param[in] model_dir Centre direction of extended model. - * @param[in] theta_max Maximum theta angle (radians). * @param[out] gradients Gradients matrix. - * @return Instrument response to radial source for all events in + * @return Instrument response to elliptical source for all events in * observation (\f$cm^2\f$). * - * Returns the instrument response to an extended source for all events in - * the observations. - * - * The computation is done by integrating the extended model for each pixel - * in Chi and Psi over a circular region centred on the Chi/Psi pixel with - * a radius equal to the maximum Phigeo value. The integration is done in - * a spherical system, spanned by Phigeo and azimuth angle. Only coordiates - * are considered that plausibly overlap with the extended model. The - * integration is done by numerical summation. The integration step is by - * default set to 1/10 of the model size, with the constraint that it never - * becomes larger than the IAQ Phigeo binsize. This should allow to capture - * all extended models with sufficient precision in an acceptable amount of - * computing time. - * - * @p gradients is an optional matrix where the number of rows corresponds - * to the number of events in the observation and the number of columns - * corresponds to the number of spatial model parameters. Since no gradients - * are computed in this method the content of @p gradients if not altered. + * Returns the instrument response to an elliptical source for all events in + * the observations. See GCOMResponse::irf_extended for more information. ***************************************************************************/ -#if defined(G_USE_MODEL_INTEGRAL_IN_IRF_EXTENDED) -GVector GCOMResponse::irf_extended(const GModelSky& model, - const GObservation& obs, - const GSkyDir& model_dir, - const double& theta_max, - GMatrix* gradients) const +GVector GCOMResponse::irf_elliptical(const GModelSky& model, + const GObservation& obs, + GMatrix* gradients) const { // Set constants static const int iter_rho = 5; @@ -1273,7 +1149,7 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL " "observations. Please specify a COMPTEL observation " "as argument."; - throw GException::invalid_argument(G_IRF_EXTENDED, msg); + throw GException::invalid_argument(G_IRF_ELLIPTICAL, msg); } // Extract COMPTEL event cube @@ -1282,22 +1158,30 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does " "not contain a COMPTEL event cube. Please specify " "a COMPTEL observation containing and event cube."; - throw GException::invalid_argument(G_IRF_EXTENDED, msg); + throw GException::invalid_argument(G_IRF_ELLIPTICAL, msg); } // Throw an exception if COMPTEL response is not set or if if (m_iaq.empty()) { std::string msg = "COMPTEL response is empty. Please initialise the " "response with an \"IAQ\"."; - throw GException::invalid_value(G_IRF_EXTENDED, msg); + throw GException::invalid_value(G_IRF_ELLIPTICAL, msg); } else if (m_phigeo_bin_size == 0.0) { std::string msg = "COMPTEL response has a zero Phigeo bin size. " "Please initialise the response with a valid " "\"IAQ\"."; - throw GException::invalid_value(G_IRF_EXTENDED, msg); + throw GException::invalid_value(G_IRF_ELLIPTICAL, msg); } + // Get pointer to elliptical model + const GModelSpatialElliptical* elliptical = + static_cast(model.spatial()); + + // Get references to model direction and maximum model radius + const GSkyDir& model_dir = elliptical->dir(); + const double& theta_max = elliptical->theta_max(); + // Get number of Chi/Psi pixels, Phibar layers and event bins int npix = cube->naxis(0) * cube->naxis(1); int nphibar = cube->naxis(2); @@ -1356,16 +1240,16 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, GVector irf(nphibar); // Setup integration kernel - com_extended_kerns_rho integrands(m_iaq, - model, - irf, - bin, - rot, - drx, - phigeo_bin, - m_phigeo_bins, - nphibar, - iter_omega); + com_elliptical_kerns_rho integrands(m_iaq, + model, + irf, + bin, + rot, + drx, + phigeo_bin, + m_phigeo_bins, + nphibar, + iter_omega); // Setup integrator GIntegrals integral(&integrands); @@ -1387,19 +1271,42 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, // Return IRF vector return irfs; } -#else -GVector GCOMResponse::irf_extended(const GModelSky& model, - const GObservation& obs, - const GSkyDir& model_dir, - const double& theta_max, - GMatrix* gradients) const -{ - // Set constants - //const double step_fraction = 0.1; - const double step_fraction = 0.05; - const int iter_phigeo = 5; - const int iter_phi = 5; + +/***********************************************************************//** + * @brief Return instrument response to diffuse source + * + * @param[in] model Sky model. + * @param[in] obs Observation. + * @param[out] gradients Gradients matrix. + * @return Instrument response to diffuse source for all events in + * observation (\f$cm^2\f$). + * + * @exception GException::invalid_argument + * Observation is not a COMPTEL observation. + * @exception GException::invalid_value + * Response not initialised with a valid IAQ + * + * Returns the instrument response to a diffuse source for all events in + * the observations. The diffuse source may be energy dependent. + * + * The computation is done by integrating the diffuse model for each pixel + * in Chi and Psi over a circular region centred on the Chi/Psi pixel with + * a radius equal to the maximum Phigeo value. The radial integration is + * done by looping over all Phigeo bins of the response. For each Phigeo + * value, the azimuthal integration is done by stepping with an angular step + * size that corresponds to the Phigeo step size (which typically is 1 deg). + * + * @p gradients is an optional matrix where the number of rows corresponds + * to the number of events in the observation and the number of columns + * corresponds to the number of spatial model parameters. Since for point + * sources no gradients are computed, the method does not alter the + * content of @p gradients. + ***************************************************************************/ +GVector GCOMResponse::irf_diffuse(const GModelSky& model, + const GObservation& obs, + GMatrix* gradients) const +{ // Extract COMPTEL observation const GCOMObservation* obs_ptr = dynamic_cast(&obs); if (obs_ptr == NULL) { @@ -1407,7 +1314,7 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL " "observations. Please specify a COMPTEL observation " "as argument."; - throw GException::invalid_argument(G_IRF_EXTENDED, msg); + throw GException::invalid_argument(G_IRF_DIFFUSE, msg); } // Extract COMPTEL event cube @@ -1416,20 +1323,20 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does " "not contain a COMPTEL event cube. Please specify " "a COMPTEL observation containing and event cube."; - throw GException::invalid_argument(G_IRF_EXTENDED, msg); + throw GException::invalid_argument(G_IRF_DIFFUSE, msg); } // Throw an exception if COMPTEL response is not set or if if (m_iaq.empty()) { std::string msg = "COMPTEL response is empty. Please initialise the " "response with an \"IAQ\"."; - throw GException::invalid_value(G_IRF_EXTENDED, msg); + throw GException::invalid_value(G_IRF_DIFFUSE, msg); } else if (m_phigeo_bin_size == 0.0) { std::string msg = "COMPTEL response has a zero Phigeo bin size. " "Please initialise the response with a valid " "\"IAQ\"."; - throw GException::invalid_value(G_IRF_EXTENDED, msg); + throw GException::invalid_value(G_IRF_DIFFUSE, msg); } // Get number of Chi/Psi pixels, Phibar layers and event bins @@ -1441,23 +1348,16 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, GVector irfs(nevents); // Initialise some variables - double phigeo_min = m_phigeo_min * gammalib::deg2rad; - double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad; - double phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins; - const GSkyMap* drx = &(obs_ptr->drx().map()); - const double* drg = obs_ptr->drg().map().pixels(); - double cos_theta_max = std::cos(theta_max); + double phigeo_min = m_phigeo_min * gammalib::deg2rad; + double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad; + double omega0 = 2.0 * std::sin(0.5 * phigeo_bin); + const GSkyMap& drx = obs_ptr->drx().map(); + const double* drg = obs_ptr->drg().map().pixels(); // Compute IAQ normalisation (1/s): DEADC / ONTIME (s) double iaq_norm = obs_ptr->deadc() / (obs_ptr->ontime() * cube->dre().tof_correction()); - // Code option: multiply DRX at model direction into IAQ normalisation - #if defined(G_APPROXIMATE_DRX_IN_IRF_EXTENDED) - iaq_norm *= (*drx)(model_dir); - drx = NULL; - #endif - // Loop over Chi and Psi for (int ipix = 0; ipix < npix; ++ipix) { @@ -1470,150 +1370,32 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, // Get reference to Phigeo sky direction const GSkyDir& phigeoDir = obsDir.dir(); - // Compute angle between model centre and Phigeo sky direction (radians) - // and position angle of model with respect to Phigeo sky direction (radians) - double zeta = phigeoDir.dist(model_dir); - double phi0 = phigeoDir.posang(model_dir); - - // Setup rotation matrix for faster conversion towards sky coordinates - GMatrix ry; - GMatrix rz; - ry.eulery(phigeoDir.dec_deg() - 90.0); - rz.eulerz(-phigeoDir.ra_deg()); - GMatrix rot = (ry * rz).transpose(); - - // Setup azimuth angle range computation - double sin_zeta = std::sin(zeta); - double cos_zeta = std::cos(zeta); - - // Setup Phigeo integration range - double phigeo_low = zeta - theta_max; - double phigeo_up = zeta + theta_max; - if (phigeo_low < phigeo_min) { - phigeo_low = phigeo_min; - } - if (phigeo_up > phigeo_max) { - phigeo_up = phigeo_max; - } - - // - // Code option A: Use integral - // =========================== - #if defined(G_USE_INTEGRAL_IN_IRF_EXTENDED) - // Perform Phigeo angle integration if interval is valid - if (phigeo_up > phigeo_low) { - - // Initialise IRF vector - GVector irf(nphibar); - - // Setup integration kernel - com_extended_kerns_phigeo integrands(m_iaq, - model, - irf, - bin->energy(), - bin->time(), - rot, - drx, - phigeo_bin, - m_phigeo_bins, - nphibar, - ipix, - npix, - iaq_norm, - zeta, - phi0, - theta_max, - iter_phi); - - // Setup integrator - GIntegrals integral(&integrands); - integral.fixed_iter(iter_phigeo); - - // Integrate over Phigeo - irf = integral.romberg(phigeo_low, phigeo_up, iter_phigeo); - - // Add IRF to result - for (int iphibar = 0; iphibar < nphibar; ++iphibar) { - int idri = ipix + iphibar * npix; - irfs[idri] += irf[iphibar] * drg[idri]; - } - - } // endif: Phigeo angle interval was valid - - // - // Code option B: Perform numerical integration - // ============================================ - #else - // Setup Phigeo integration step-size and number of integration steps. By - // default the integration step size is one tenth of the model diameter, but - // not larger than the Phigeo binning of the IAQ. We also make sure that there - // are at least two Phigeo integration steps. - double phigeo_length = phigeo_up - phigeo_low; - double phigeo_step = step_fraction * phigeo_length; - if (phigeo_step > phigeo_bin) { - phigeo_step = phigeo_bin; - } - int nphigeo = int(phigeo_length/phigeo_step + 0.5); - if (nphigeo < 2) { - nphigeo = 2; - } - phigeo_step = phigeo_length / double(nphigeo); - - // Compute solid-angle for selected Phigeo integration step size - double omega0 = 2.0 * std::sin(0.5 * phigeo_step); - - // Integrate over Phigeo - double phigeo = phigeo_low; - for (int iphigeo = 0; iphigeo <= nphigeo; ++iphigeo, phigeo += phigeo_step) { - - // Compute half length of the arc (in radians) from a circle with - // radius phigeo that intersects with the model, defined as a circle - // with maximum radius m_theta_max - double dphi = 0.5 * gammalib::roi_arclength(phigeo, - zeta, - cos_zeta, - sin_zeta, - theta_max, - cos_theta_max); - - // If arc length is not positive then skip this Phigeo value - if (dphi <= 0.0) { - continue; - } + // Loop over Phigeo + for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) { - // Compute sine and cosine of Phigeo + // Determine Phigeo value in radians + double phigeo = phigeo_min + iphigeo * phigeo_bin; double sin_phigeo = std::sin(phigeo); - double cos_phigeo = std::cos(phigeo); // Determine number of azimuthal integration steps and step size - // in radians. By default the integration step size is one tenth of - // the arclength, but larger than the Phigeo binning of the IAQ. We - // also make sure that there are at least two integration steps. - double azimuth_length = 2.0 * dphi * sin_phigeo; - double az_step = step_fraction * azimuth_length; - if (az_step > phigeo_bin) { - az_step = phigeo_bin; - } - int naz = int(azimuth_length/az_step + 0.5); + // in radians + double length = gammalib::twopi * sin_phigeo; + int naz = int(length / phigeo_bin + 0.5); if (naz < 2) { naz = 2; } - az_step = 2.0 * dphi / double(naz); + double az_step = gammalib::twopi / double(naz); // Computes solid angle of integration bin multiplied by Jaccobian double omega = omega0 * az_step * sin_phigeo; // Loop over azimuth angle - double az = phi0 - dphi; - for (int iaz = 0; iaz <= naz; ++iaz, az += az_step) { + double az = 0.0; + for (int iaz = 0; iaz < naz; ++iaz, az += az_step) { // Get sky direction - double sin_az = std::sin(az); - double cos_az = std::cos(az); - GVector native(-cos_az*sin_phigeo, sin_az*sin_phigeo, cos_phigeo); - GVector dir = rot * native; - GSkyDir skyDir; - skyDir.celvector(dir); + GSkyDir skyDir = phigeoDir; + skyDir.rotate(az, phigeo); // Fall through if sky direction is not contained in DRX if (!drx.contains(skyDir)) { @@ -1637,67 +1419,39 @@ GVector GCOMResponse::irf_extended(const GModelSky& model, // Multiply intensity by solid angle (unit: cm^2 s) intensity *= omega; - // Extract IAQ value by linear inter/extrapolation in Phigeo - double phirat = phigeo / phigeo_bin; // 0.5 at bin centre - int i_phigeo = int(phirat); // index into which Phigeo falls - double eps = phirat - i_phigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[ - double iaq = 0.0; - - // Continue only if Phigeo index is valid - if (i_phigeo < m_phigeo_bins) { - - // Loop over Phibar - for (int iphibar = 0; iphibar < nphibar; ++iphibar) { - - // Get IAQ index - int i = iphibar * m_phigeo_bins + i_phigeo; - - // Get interpolated IAQ value - if (eps < 0.0) { // interpolate towards left - if (i_phigeo > 0) { - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; - } - else { - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; - } - } - else { // interpolate towards right - if (i_phigeo < m_phigeo_bins-1) { - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; - } - else { - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; - } - } - - // Fall through if IAQ is not positive - if (iaq <= 0.0) { - continue; - } - - // Get DRI index - int idri = ipix + iphibar * npix; - - // Compute IRF value (unit: cm^2) - double irf = iaq * drg[idri] * iaq_norm * intensity; - - // Add IRF value if it is positive - if (irf > 0.0) { - irfs[idri] += irf; - } - - } // endfor: looped over Phibar - - } // endif: Phigeo bin was valid + // Loop over Phibar + for (int iphibar = 0; iphibar < nphibar; ++iphibar) { + + // Get IAQ index + int i = iphibar * m_phigeo_bins + iphigeo; + + // Get IAQ value + double iaq = m_iaq[i]; + + // Fall through if IAQ is not positive + if (iaq <= 0.0) { + continue; + } + + // Get DRI index + int idri = ipix + iphibar * npix; + + // Compute IRF value (unit: cm^2) + double irf = iaq * drg[idri] * iaq_norm * intensity; + + // Add IRF value if it is positive + if (irf > 0.0) { + irfs[idri] += irf; + } + + } // endfor: looped over Phibar } // endfor: looped over azimuth angle around locus of IAQ } // endfor: looped over Phigeo angles - #endif } // endfor: looped over Chi and Psi pixels // Return IRF vector return irfs; } -#endif diff --git a/inst/com/src/com_helpers_response_vector.cpp b/inst/com/src/com_helpers_response_vector.cpp index b3d76a3..6845aef 100644 --- a/inst/com/src/com_helpers_response_vector.cpp +++ b/inst/com/src/com_helpers_response_vector.cpp @@ -52,147 +52,6 @@ ==========================================================================*/ /***********************************************************************//** - * @brief Kernel for radial integration of extended model - * - * @param[in] phigeo Phigeo angle (radians). - * @return Azimuthally integrated extended model. - ***************************************************************************/ -GVector com_extended_kerns_phigeo::eval(const double& phigeo) -{ - // Initialise kernel values - m_irfs = 0.0; - - // Continue only if phigeo is positive - if (phigeo > 0.0) { - - // Compute half length of the arc (in radians) from a circle with - // radius phigeo that intersects with the model, defined as a - // circle with maximum radius m_theta_max - double dphi = 0.5 * gammalib::roi_arclength(phigeo, - m_zeta, - m_cos_zeta, - m_sin_zeta, - m_theta_max, - m_cos_theta_max); - - // Continue only if arc length is positive - if (dphi > 0.0) { - - // Compute sine and cosine of Phigeo - double sin_phigeo = std::sin(phigeo); - double cos_phigeo = std::cos(phigeo); - - // Compute phi integration range - double phi_min = m_phi0 - dphi; - double phi_max = m_phi0 + dphi; - - // Precompute interpolated IAQ vector for Phigeo angle - double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre - int iphigeo = int(phirat); // index into which Phigeo falls - double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[ - - // Continue only if Phigeo index is valid - if (iphigeo < m_phigeo_bins) { - - // Setup kernel for azimuthal integration - com_extended_kerns_phi integrand(m_model, - m_srcEng, - m_srcTime, - m_rot, - m_drx, - sin_phigeo, - cos_phigeo); - - // Setup integrator - GIntegral integral(&integrand); - integral.fixed_iter(m_iter); - - // Integrate over azimuth angle - double irf = integral.romberg(phi_min, phi_max, m_iter) * sin_phigeo; - - // Multiply result with IAQ for all Phibar layers - for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) { - - // Initialise IAQ - double iaq; - - // Get IAQ index - int i = iphibar * m_phigeo_bins + iphigeo; - - // Get interpolated IAQ value - if (eps < 0.0) { // interpolate towards left - if (iphigeo > 0) { - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; - } - else { - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; - } - } - else { // interpolate towards right - if (iphigeo < m_phigeo_bins-1) { - iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1]; - } - else { - iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1]; - } - } - - // Compute IRF - m_irfs[iphibar] = irf * iaq * m_iaq_norm; - - } // endfor: looped over Phibar - - } // endif: Phigeo bin was valid - - } // endif: arc length was positive - - } // endif: phigeo was positive - - // Return kernel values - return m_irfs; -} - - -/***********************************************************************//** - * @brief Kernel for azimuthal integration of extended model - * - * @param[in] phi Azimuth angle (radians). - * @return Vector of azimuthally integrated model. - ***************************************************************************/ -double com_extended_kerns_phi::eval(const double& phi) -{ - // Compute sine and cosine of azimuth angle - double sin_phi = std::sin(phi); - double cos_phi = std::cos(phi); - - // Get sky direction - GVector native(-cos_phi*m_sin_phigeo, sin_phi*m_sin_phigeo, m_cos_phigeo); - GVector dir = m_rot * native; - GSkyDir skyDir; - skyDir.celvector(dir); - - // Set photon - GPhoton photon(skyDir, m_srcEng, m_srcTime); - - // Get model sky intensity for photon (unit: sr^-1) - double intensity = m_model.spatial()->eval(photon); - - // Continue only if intensity is positive - if (intensity > 0.0) { - - // Multiply-in DRX value if pointer is valid - if (m_drx != NULL) { - intensity *= (*m_drx)(skyDir); - } - - } // endif: intensity was valid - - // Return kernel values - return intensity; -} - - -/***********************************************************************//** * @brief Kernel for radial integration of radial models * * @param[in] rho Rho angle (radians). @@ -320,12 +179,12 @@ GVector com_radial_kerns_omega::eval(const double& omega) return m_irfs; } /***********************************************************************//** - * @brief Kernel for radial integration of extended model + * @brief Kernel for radial integration of elliptical models * * @param[in] rho Rho angle (radians). - * @return Azimuthally integrated extended model. + * @return Azimuthally integrated elliptical model. ***************************************************************************/ -GVector com_extended_kerns_rho::eval(const double& rho) +GVector com_elliptical_kerns_rho::eval(const double& rho) { // Initialise kernel values m_irfs = 0.0; @@ -338,17 +197,17 @@ GVector com_extended_kerns_rho::eval(const double& rho) double cos_rho = std::cos(rho); // Setup azimuthal integration kernel - com_extended_kerns_omega integrands(m_iaq, - m_model, - m_irfs, - m_bin, - m_rot, - m_drx, - m_phigeo_bin_size, - m_phigeo_bins, - m_phibar_bins, - sin_rho, - cos_rho); + com_elliptical_kerns_omega integrands(m_iaq, + m_model, + m_irfs, + m_bin, + m_rot, + m_drx, + m_phigeo_bin_size, + m_phigeo_bins, + m_phibar_bins, + sin_rho, + cos_rho); // Setup integrator GIntegrals integral(&integrands); @@ -365,12 +224,12 @@ GVector com_extended_kerns_rho::eval(const double& rho) /***********************************************************************//** - * @brief Kernel for azimuthal integration of extended model + * @brief Kernel for azimuthal integration of elliptical models * * @param[in] omega Omega angle (radians). - * @return Kernel value for extended model. + * @return Kernel value for elliptical model. ***************************************************************************/ -GVector com_extended_kerns_omega::eval(const double& omega) +GVector com_elliptical_kerns_omega::eval(const double& omega) { // Initialise kernel values m_irfs = 0.0; diff --git a/inst/com/src/com_helpers_response_vector.hpp b/inst/com/src/com_helpers_response_vector.hpp index 2a2b7ce..b0118b7 100644 --- a/inst/com/src/com_helpers_response_vector.hpp +++ b/inst/com/src/com_helpers_response_vector.hpp @@ -44,116 +44,7 @@ /***********************************************************************//** - * @class com_extended_kerns_phigeo - * - * @brief Kernel for Phigeo angle integration of extended model - * - * This class provides the kernel for the Phigeo angle integration of - * extended sky models. - ***************************************************************************/ -class com_extended_kerns_phigeo : public GFunctions { -public: - com_extended_kerns_phigeo(const std::vector& iaq, - const GModelSky& model, - GVector& irfs, - const GEnergy& srcEng, - const GTime& srcTime, - const GMatrix& rot, - const GSkyMap* drx, - const double& phigeo_bin_size, - const int& phigeo_bins, - const int& phibar_bins, - const int& ipix, - const int& npix, - const double& iaq_norm, - const double& zeta, - const double& phi0, - const double& theta_max, - const int& iter) : - m_iaq(iaq), - m_model(model), - m_irfs(irfs), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_rot(rot), - m_drx(drx), - m_phigeo_bin_size(phigeo_bin_size), - m_phigeo_bins(phigeo_bins), - m_phibar_bins(phibar_bins), - m_ipix(ipix), - m_npix(npix), - m_iaq_norm(iaq_norm), - m_zeta(zeta), - m_sin_zeta(std::sin(zeta)), - m_cos_zeta(std::cos(zeta)), - m_phi0(phi0), - m_theta_max(theta_max), - m_cos_theta_max(std::cos(theta_max)), - m_iter(iter) { } - int size(void) const { return m_irfs.size(); } - GVector eval(const double& phigeo); -protected: - const std::vector& m_iaq; //!< IAQ vector - const GModelSky& m_model; //!< Sky model - GVector& m_irfs; //!< IRF vector to update - const GEnergy& m_srcEng; //!< Source energy - const GTime& m_srcTime; //!< Source time - const GMatrix& m_rot; //!< Rotation matrix - const GSkyMap* m_drx; //!< DRX - const double& m_phigeo_bin_size; //!< Phigeo bin size - const int& m_phigeo_bins; //!< Number of phigeo bins - const int& m_phibar_bins; //!< Number of phibar bins - const int& m_ipix; //!< DRI pixel index - const int& m_npix; //!< Number of DRI pixels - const double& m_iaq_norm; //!< IAQ normalisation - const double& m_zeta; //!< Zeta angle - double m_sin_zeta; //!< Sine of zeta angle - double m_cos_zeta; //!< Cosine of zeta angle - const double& m_phi0; //!< Model position angle - const double& m_theta_max; //!< Maximum model radius - double m_cos_theta_max; //!< Cosine of max. mod. rad. - const int& m_iter; //!< Number of phi iterations -}; - - -/***********************************************************************//** - * @class com_extended_kerns_phi - * - * @brief Kernel for azimuth angle integration of extended model - * - * This class provides the kernel for the azimuth angle integration of - * extended sky models. - ***************************************************************************/ -class com_extended_kerns_phi : public GFunction { -public: - com_extended_kerns_phi(const GModelSky& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GMatrix& rot, - const GSkyMap* drx, - const double& sin_phigeo, - const double& cos_phigeo) : - m_model(model), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_rot(rot), - m_drx(drx), - m_sin_phigeo(sin_phigeo), - m_cos_phigeo(cos_phigeo) { } - double eval(const double& phi); -protected: - const GModelSky& m_model; //!< Sky model - const GEnergy& m_srcEng; //!< Source energy - const GTime& m_srcTime; //!< Source time - const GMatrix& m_rot; //!< Rotation matrix - const GSkyMap* m_drx; //!< DRX - const double& m_sin_phigeo; //!< Sine of Phigeo - const double& m_cos_phigeo; //!< Cosine of Phigeo -}; - - -/***********************************************************************//** - * @class com_extended_kerns_rho + * @class com_radial_kerns_rho * * @brief Kernel for rho angle integration of radial models * @@ -199,12 +90,12 @@ protected: /***********************************************************************//** - * @class com_extended_kerns_omega + * @class com_radial_kerns_omega * - * @brief Kernel for azimuth angle integration of extended models + * @brief Kernel for azimuth angle integration of radial models * * This class provides the kernel for the azimuth angle integration of - * extended sky models. + * radial sky models. ***************************************************************************/ class com_radial_kerns_omega : public GFunctions { public: @@ -245,35 +136,35 @@ protected: /***********************************************************************//** - * @class com_radial_kerns_rho + * @class com_elliptical_kerns_rho * - * @brief Kernel for rho angle integration of extended models + * @brief Kernel for rho angle integration of elliptical models * - * This class provides the kernel for the rho angle integration of extended + * This class provides the kernel for the rho angle integration of elliptical * sky models. ***************************************************************************/ -class com_extended_kerns_rho : public GFunctions { +class com_elliptical_kerns_rho : public GFunctions { public: - com_extended_kerns_rho(const std::vector& iaq, - const GModelSky& model, - GVector& irfs, - const GCOMEventBin* bin, - const GMatrix& rot, - const GSkyMap& drx, - const double& phigeo_bin_size, - const int& phigeo_bins, - const int& phibar_bins, - const int& iter) : - m_iaq(iaq), - m_model(model), - m_irfs(irfs), - m_bin(bin), - m_rot(rot), - m_drx(drx), - m_phigeo_bin_size(phigeo_bin_size), - m_phigeo_bins(phigeo_bins), - m_phibar_bins(phibar_bins), - m_iter(iter) { } + com_elliptical_kerns_rho(const std::vector& iaq, + const GModelSky& model, + GVector& irfs, + const GCOMEventBin* bin, + const GMatrix& rot, + const GSkyMap& drx, + const double& phigeo_bin_size, + const int& phigeo_bins, + const int& phibar_bins, + const int& iter) : + m_iaq(iaq), + m_model(model), + m_irfs(irfs), + m_bin(bin), + m_rot(rot), + m_drx(drx), + m_phigeo_bin_size(phigeo_bin_size), + m_phigeo_bins(phigeo_bins), + m_phibar_bins(phibar_bins), + m_iter(iter) { } int size(void) const { return m_irfs.size(); } GVector eval(const double& phigeo); protected: @@ -291,37 +182,37 @@ protected: /***********************************************************************//** - * @class com_radial_kerns_omega + * @class com_elliptical_kerns_omega * - * @brief Kernel for azimuth angle integration of extended models + * @brief Kernel for azimuth angle integration of elliptical models * * This class provides the kernel for the azimuth angle integration of - * extended sky models. + * elliptical sky models. ***************************************************************************/ -class com_extended_kerns_omega : public GFunctions { +class com_elliptical_kerns_omega : public GFunctions { public: - com_extended_kerns_omega(const std::vector& iaq, - const GModelSky& model, - GVector& irfs, - const GCOMEventBin* bin, - const GMatrix& rot, - const GSkyMap& drx, - const double& phigeo_bin_size, - const int& phigeo_bins, - const int& phibar_bins, - const double& sin_rho, - const double& cos_rho) : - m_iaq(iaq), - m_model(model), - m_irfs(irfs), - m_bin(bin), - m_rot(rot), - m_drx(drx), - m_phigeo_bin_size(phigeo_bin_size), - m_phigeo_bins(phigeo_bins), - m_phibar_bins(phibar_bins), - m_sin_rho(sin_rho), - m_cos_rho(cos_rho) { } + com_elliptical_kerns_omega(const std::vector& iaq, + const GModelSky& model, + GVector& irfs, + const GCOMEventBin* bin, + const GMatrix& rot, + const GSkyMap& drx, + const double& phigeo_bin_size, + const int& phigeo_bins, + const int& phibar_bins, + const double& sin_rho, + const double& cos_rho) : + m_iaq(iaq), + m_model(model), + m_irfs(irfs), + m_bin(bin), + m_rot(rot), + m_drx(drx), + m_phigeo_bin_size(phigeo_bin_size), + m_phigeo_bins(phigeo_bins), + m_phibar_bins(phibar_bins), + m_sin_rho(sin_rho), + m_cos_rho(cos_rho) { } int size(void) const { return m_irfs.size(); } GVector eval(const double& omega); protected: