Justine Devin / gammalib

Commit e98298d1dfedc6a87cf610df9ce3ddbcf4ff2dda

Authored by Jürgen Knödlseder
1 parent 1fc29539

Implement npred_elliptical method.

Inline Side-by-side
Showing 2 changed files with 221 additions and 16 deletions
include/GResponse.hpp
  View file @e98298d
... ... @@ -36,6 +36,7 @@
36 36 #include "GEnergy.hpp"
37 37 #include "GTime.hpp"
38 38 #include "GModelSpatialRadial.hpp"
  39 +#include "GModelSpatialElliptical.hpp"
39 40 #include "GFunction.hpp"
40 41 #include "GMatrix.hpp"
41 42  
... ... @@ -168,6 +169,64 @@ protected:
168 169 double m_sin_theta; //!< Sine of offset angle
169 170 };
170 171  
  172 + // Npred theta integration kernel for elliptical model
  173 + class npred_elliptical_kern_theta : public GFunction {
  174 + public:
  175 + npred_elliptical_kern_theta(const GResponse* rsp,
  176 + const GModelSpatialElliptical* spatial,
  177 + const GEnergy* srcEng,
  178 + const GTime* srcTime,
  179 + const GObservation* obs,
  180 + const GMatrix* rot) :
  181 + m_rsp(rsp),
  182 + m_spatial(spatial),
  183 + m_srcEng(srcEng),
  184 + m_srcTime(srcTime),
  185 + m_obs(obs),
  186 + m_rot(rot) { }
  187 + double eval(double theta);
  188 + protected:
  189 + const GResponse* m_rsp; //!< Pointer to response
  190 + const GModelSpatialElliptical* m_spatial; //!< Pointer to spatial model
  191 + const GEnergy* m_srcEng; //!< Pointer to true photon energy
  192 + const GTime* m_srcTime; //!< Pointer to true photon arrival time
  193 + const GObservation* m_obs; //!< Pointer to observation
  194 + const GMatrix* m_rot; //!< Pointer to rotation matrix
  195 + };
  196 +
  197 + // Npred phi integration kernel for elliptical model
  198 + class npred_elliptical_kern_phi : public GFunction {
  199 + public:
  200 + npred_elliptical_kern_phi(const GResponse* rsp,
  201 + const GModelSpatialElliptical* spatial,
  202 + const GEnergy* srcEng,
  203 + const GTime* srcTime,
  204 + const GObservation* obs,
  205 + const GMatrix* rot,
  206 + double theta,
  207 + double sin_theta) :
  208 + m_rsp(rsp),
  209 + m_spatial(spatial),
  210 + m_srcEng(srcEng),
  211 + m_srcTime(srcTime),
  212 + m_obs(obs),
  213 + m_rot(rot),
  214 + m_theta(theta),
  215 + m_cos_theta(std::cos(theta)),
  216 + m_sin_theta(sin_theta) { }
  217 + double eval(double phi);
  218 + protected:
  219 + const GResponse* m_rsp; //!< Pointer to response
  220 + const GModelSpatialElliptical* m_spatial; //!< Pointer to spatial model
  221 + const GEnergy* m_srcEng; //!< Pointer to true photon energy
  222 + const GTime* m_srcTime; //!< Pointer to true photon arrival time
  223 + const GObservation* m_obs; //!< Pointer to observation
  224 + const GMatrix* m_rot; //!< Pointer to rotation matrix
  225 + double m_theta; //!< Offset angle (radians)
  226 + double m_cos_theta; //!< cosine of offset angle
  227 + double m_sin_theta; //!< Sine of offset angle
  228 + };
  229 +
171 230 };
172 231  
173 232 #endif /* GRESPONSE_HPP */
... ...
src/obs/GResponse.cpp
  View file @e98298d
... ... @@ -390,7 +390,7 @@ double GResponse::npred_ptsrc(const GSource&amp; source,
390 390  
391 391  
392 392 /***********************************************************************//**
393   - * @brief Return ROI integral of radial source model
  393 + * @brief Return spatial integral of radial source model
394 394 *
395 395 * @param[in] source Source.
396 396 * @param[in] obs Observation.
... ... @@ -406,32 +406,32 @@ double GResponse::npred_radial(const GSource&amp; source,
406 406 // Initialise Npred value
407 407 double npred = 0.0;
408 408  
409   - // Get radial model
410   - const GModelSpatialRadial* src =
  409 + // Get radial spatial model
  410 + const GModelSpatialRadial* spatial =
411 411 dynamic_cast<const GModelSpatialRadial*>(source.model());
412 412  
413   - // Continue only if model is valid
414   - if (src != NULL) {
  413 + // Continue only if spatial model is valid
  414 + if (spatial != NULL) {
415 415  
416 416 // Compute rotation matrix to convert from native coordinates given
417 417 // by (theta,phi) into celestial coordinates.
418 418 GMatrix ry;
419 419 GMatrix rz;
420 420 GMatrix rot;
421   - ry.eulery(src->dec() - 90.0);
422   - rz.eulerz(-src->ra());
  421 + ry.eulery(spatial->dec() - 90.0);
  422 + rz.eulerz(-spatial->ra());
423 423 rot = transpose(ry * rz);
424 424  
425 425 // Set offset angle integration range
426 426 double theta_min = 0.0;
427   - double theta_max = src->theta_max();
  427 + double theta_max = spatial->theta_max();
428 428  
429 429 // Perform offset angle integration if interval is valid
430 430 if (theta_max > theta_min) {
431 431  
432 432 // Setup integration kernel
433 433 GResponse::npred_radial_kern_theta integrand(this,
434   - src,
  434 + spatial,
435 435 &(source.energy()),
436 436 &(source.time()),
437 437 &obs,
... ... @@ -475,19 +475,74 @@ double GResponse::npred_radial(const GSource&amp; source,
475 475 *
476 476 * @param[in] source Source.
477 477 * @param[in] obs Observation.
478   - *
479   - * @exception GException::feature_not_implemented
480   - * Method not yet implemented.
481 478 ***************************************************************************/
482 479 double GResponse::npred_elliptical(const GSource& source,
483 480 const GObservation& obs) const
484 481 {
485   - // Feature not yet implemented
486   - throw GException::feature_not_implemented(G_NPRED_ELLIPTICAL,
487   - "Npred computation not implemented for elliptical models.");
  482 + // Initialise Npred value
  483 + double npred = 0.0;
  484 +
  485 + // Get elliptical spatial model
  486 + const GModelSpatialElliptical* spatial =
  487 + dynamic_cast<const GModelSpatialElliptical*>(source.model());
  488 +
  489 + // Continue only if spatial model is valid
  490 + if (spatial != NULL) {
  491 +
  492 + // Compute rotation matrix to convert from native coordinates given
  493 + // by (theta,phi) into celestial coordinates.
  494 + GMatrix ry;
  495 + GMatrix rz;
  496 + GMatrix rot;
  497 + ry.eulery(spatial->dec() - 90.0);
  498 + rz.eulerz(-spatial->ra());
  499 + rot = transpose(ry * rz);
  500 +
  501 + // Set offset angle integration range
  502 + double theta_min = 0.0;
  503 + double theta_max = spatial->theta_max();
  504 +
  505 + // Perform offset angle integration if interval is valid
  506 + if (theta_max > theta_min) {
  507 +
  508 + // Setup integration kernel
  509 + GResponse::npred_elliptical_kern_theta integrand(this,
  510 + spatial,
  511 + &(source.energy()),
  512 + &(source.time()),
  513 + &obs,
  514 + &rot);
  515 +
  516 + // Integrate over theta
  517 + GIntegral integral(&integrand);
  518 + npred = integral.romb(theta_min, theta_max);
  519 +
  520 + // Compile option: Show integration results
  521 + #if defined(G_DEBUG_NPRED_ELLIPTICAL)
  522 + std::cout << "GResponse::npred_elliptical:";
  523 + std::cout << " theta_min=" << theta_min;
  524 + std::cout << " theta_max=" << theta_max;
  525 + std::cout << " npred=" << npred << std::endl;
  526 + #endif
  527 +
  528 + } // endif: offset angle range was valid
  529 +
  530 + // Debug: Check for NaN
  531 + #if defined(G_NAN_CHECK)
  532 + if (isnotanumber(npred) || isinfinite(npred)) {
  533 + std::cout << "*** ERROR: GResponse::npred_elliptical:";
  534 + std::cout << " NaN/Inf encountered";
  535 + std::cout << " (npred=" << npred;
  536 + std::cout << ", theta_min=" << theta_min;
  537 + std::cout << ", theta_max=" << theta_max;
  538 + std::cout << ")" << std::endl;
  539 + }
  540 + #endif
  541 +
  542 + } // endif: elliptical model was valid
488 543  
489 544 // Return Npred
490   - return 0.0;
  545 + return npred;
491 546 }
492 547  
493 548  
... ... @@ -638,3 +693,94 @@ double GResponse::npred_radial_kern_phi::eval(double phi)
638 693 // Return Npred
639 694 return npred;
640 695 }
  696 +
  697 +
  698 +/***********************************************************************//**
  699 + * @brief Kernel for offset angle Npred integration of elliptical model
  700 + *
  701 + * @param[in] theta Radial model offset angle (radians).
  702 + ***************************************************************************/
  703 +double GResponse::npred_elliptical_kern_theta::eval(double theta)
  704 +{
  705 + // Initialise Npred value
  706 + double npred = 0.0;
  707 +
  708 + // Compute sine of offset angle
  709 + double sin_theta = std::sin(theta);
  710 +
  711 + // Setup phi integration kernel
  712 + GResponse::npred_elliptical_kern_phi integrand(m_rsp,
  713 + m_spatial,
  714 + m_srcEng,
  715 + m_srcTime,
  716 + m_obs,
  717 + m_rot,
  718 + theta,
  719 + sin_theta);
  720 +
  721 + // Integrate over phi
  722 + GIntegral integral(&integrand);
  723 + npred = integral.romb(0.0, twopi) * sin_theta;
  724 +
  725 + // Debug: Check for NaN
  726 + #if defined(G_NAN_CHECK)
  727 + if (isnotanumber(npred) || isinfinite(npred)) {
  728 + std::cout << "*** ERROR: GResponse::npred_radial_kern_theta::eval";
  729 + std::cout << "(theta=" << theta << "):";
  730 + std::cout << " NaN/Inf encountered";
  731 + std::cout << " (npred=" << npred;
  732 + std::cout << ", sin_theta=" << sin_theta;
  733 + std::cout << ")" << std::endl;
  734 + }
  735 + #endif
  736 +
  737 + // Return Npred
  738 + return npred;
  739 +}
  740 +
  741 +
  742 +/***********************************************************************//**
  743 + * @brief Kernel for azimuth angle Npred integration of elliptical model
  744 + *
  745 + * @param[in] phi Azimuth angle (radians).
  746 + ***************************************************************************/
  747 +double GResponse::npred_elliptical_kern_phi::eval(double phi)
  748 +{
  749 + // Compute sky direction vector in native coordinates
  750 + double cos_phi = std::cos(phi);
  751 + double sin_phi = std::sin(phi);
  752 + GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta);
  753 +
  754 + // Rotate from native into celestial system
  755 + GVector cel = *m_rot * native;
  756 +
  757 + // Set sky direction
  758 + GSkyDir srcDir;
  759 + srcDir.celvector(cel);
  760 +
  761 + // Set photon
  762 + GPhoton photon(srcDir, *m_srcEng, *m_srcTime);
  763 +
  764 + // Get elliptical model value
  765 + double model = m_spatial->eval(m_theta, phi);
  766 +
  767 + // Compute Npred for this sky direction
  768 + double npred = m_rsp->npred(photon, *m_obs) * model;
  769 +
  770 + // Debug: Check for NaN
  771 + #if defined(G_NAN_CHECK)
  772 + if (isnotanumber(npred) || isinfinite(npred)) {
  773 + std::cout << "*** ERROR: GResponse::npred_radial_kern_phi::eval";
  774 + std::cout << "(phi=" << phi << "):";
  775 + std::cout << " NaN/Inf encountered";
  776 + std::cout << " (npred=" << npred;
  777 + std::cout << ", model=" << model;
  778 + std::cout << ", cos_phi=" << cos_phi;
  779 + std::cout << ", sin_phi=" << sin_phi;
  780 + std::cout << ")" << std::endl;
  781 + }
  782 + #endif
  783 +
  784 + // Return Npred
  785 + return npred;
  786 +}
... ...