diff --git a/include/GResponse.hpp b/include/GResponse.hpp index 3adaae1..26713e7 100644 --- a/include/GResponse.hpp +++ b/include/GResponse.hpp @@ -24,9 +24,15 @@ #include #include "GLog.hpp" #include "GEvent.hpp" +#include "GSkyDir.hpp" #include "GEnergy.hpp" #include "GTime.hpp" #include "GModelSky.hpp" +#include "GModelPointSource.hpp" +#include "GModelExtendedSource.hpp" +#include "GModelDiffuseSource.hpp" +#include "GIntegrand.hpp" +#include "GMatrix.hpp" /* __ Forward declarations _______________________________________________ */ class GObservation; @@ -66,19 +72,125 @@ public: virtual GResponse* clone(void) const = 0; virtual bool hasedisp(void) const = 0; virtual bool hastdisp(void) const = 0; - virtual double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, const GObservation& obs) const = 0; - virtual double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, const GObservation& obs) const = 0; virtual std::string print(void) const = 0; + // Virtual methods + virtual double irf(const GEvent& event, + const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_ptsrc(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_extended(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_diffuse(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred(const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_ptsrc(const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_extended(const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_diffuse(const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + protected: // Protected methods void init_members(void); void copy_members(const GResponse& rsp); void free_members(void); + + // Npred theta integration kernel + class npred_kern_theta : public GIntegrand { + public: + npred_kern_theta(const GResponse* rsp, + const GModelRadial* radial, + const GEnergy* srcEng, + const GTime* srcTime, + const GObservation* obs, + const GMatrix* rot) : + m_rsp(rsp), + m_radial(radial), + m_srcEng(srcEng), + m_srcTime(srcTime), + m_obs(obs), + m_rot(rot) { return; } + double eval(double theta); + protected: + const GResponse* m_rsp; //!< Pointer to response + const GModelRadial* m_radial; //!< Pointer to radial spatial model + const GEnergy* m_srcEng; //!< Pointer to true photon energy + const GTime* m_srcTime; //!< Pointer to true photon arrival time + const GObservation* m_obs; //!< Pointer to observation + const GMatrix* m_rot; //!< Rotation matrix + }; + + // Npred phi integration kernel + class npred_kern_phi : public GIntegrand { + public: + npred_kern_phi(const GResponse* rsp, + const GEnergy* srcEng, + const GTime* srcTime, + const GObservation* obs, + const GMatrix* rot, + double theta, + double sin_theta) : + m_rsp(rsp), + m_srcEng(srcEng), + m_srcTime(srcTime), + m_obs(obs), + m_rot(rot), + m_theta(theta), + m_cos_theta(std::cos(theta)), + m_sin_theta(sin_theta) { return; } + double eval(double phi); + protected: + const GResponse* m_rsp; //!< Pointer to response + const GModelRadial* m_radial; //!< Pointer to radial spatial model + const GEnergy* m_srcEng; //!< Pointer to true photon energy + const GTime* m_srcTime; //!< Pointer to true photon arrival time + const GObservation* m_obs; //!< Pointer to observation + const GMatrix* m_rot; //!< Rotation matrix + double m_theta; //!< Offset angle (radians) + double m_cos_theta; //!< cosine of offset angle + double m_sin_theta; //!< Sine of offset angle + }; + }; #endif /* GRESPONSE_HPP */ diff --git a/inst/cta/include/GCTAResponse.hpp b/inst/cta/include/GCTAResponse.hpp index b7151d5..39a4020 100644 --- a/inst/cta/include/GCTAResponse.hpp +++ b/inst/cta/include/GCTAResponse.hpp @@ -12,7 +12,7 @@ ***************************************************************************/ /** * @file GCTAResponse.hpp - * @brief GCTAResponse class interface definition. + * @brief CTA instrument response function class interface definition * @author J. Knodlseder */ @@ -73,13 +73,35 @@ public: virtual GCTAResponse* clone(void) const; virtual bool hasedisp(void) const { return false; } virtual bool hastdisp(void) const { return false; } - virtual double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, const GObservation& obs) const; - virtual double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, const GObservation& obs) const; + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; virtual std::string print(void) const; + // Overload virtual base class methods + virtual double irf_extended(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_diffuse(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + // Other Methods GCTAEventAtom* mc(const double& area, const GPhoton& photon, const GPointing& pnt, GRan& ran) const; @@ -89,41 +111,6 @@ public: void eps(const double& eps) { m_eps=eps; } const double& eps(void) const { return m_eps; } - // Model type dependent response methods - double irf_ptsrc(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double irf_extended(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelExtendedSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double irf_diffuse(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelDiffuseSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_ptsrc(const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_extended(const GModelExtendedSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_diffuse(const GModelDiffuseSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - // Low-level response methods double aeff(const double& theta, const double& phi, diff --git a/inst/cta/pyext/GCTAResponse.i b/inst/cta/pyext/GCTAResponse.i index 98f6dbc..f0081c0 100644 --- a/inst/cta/pyext/GCTAResponse.i +++ b/inst/cta/pyext/GCTAResponse.i @@ -27,6 +27,7 @@ * @brief CTA instrument response function ***************************************************************************/ class GCTAResponse : public GResponse { + public: // Constructors and destructors GCTAResponse(void); @@ -39,11 +40,33 @@ public: virtual GCTAResponse* clone(void) const; virtual bool hasedisp(void) const; virtual bool hastdisp(void) const; - virtual double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, const GObservation& obs) const; - virtual double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, const GObservation& obs) const; + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + + // Overload virtual base class methods + virtual double irf_extended(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_diffuse(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; // Other Methods GCTAEventAtom* mc(const double& area, const GPhoton& photon, @@ -54,41 +77,6 @@ public: void eps(const double& eps); const double& eps(void) const; - // Model type dependent response methods - double irf_ptsrc(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double irf_extended(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelExtendedSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double irf_diffuse(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelDiffuseSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_ptsrc(const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_extended(const GModelExtendedSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - double npred_diffuse(const GModelDiffuseSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const; - // Low-level response methods double aeff(const double& theta, const double& phi, @@ -126,6 +114,7 @@ public: // Analytical PSF implementation double psf_dummy(const double& delta, const double& sigma) const; double psf_dummy_sigma(const double& srcLogEng) const; + double psf_dummy_max(const double& sigma) const; }; diff --git a/inst/cta/src/GCTAResponse.cpp b/inst/cta/src/GCTAResponse.cpp index 9941b78..4a85f06 100644 --- a/inst/cta/src/GCTAResponse.cpp +++ b/inst/cta/src/GCTAResponse.cpp @@ -40,20 +40,14 @@ /* __ Method name definitions ____________________________________________ */ #define G_CALDB "GCTAResponse::caldb(std::string&)" -#define G_IRF "GCTAResponse::irf(GEvent&, GModelSky&, GEnergy&, GTime&,"\ - " GObservation&)" -#define G_NPRED "GCTAResponse::npred(GModelSky&, GEnergy&, GTime&,"\ +#define G_IRF "GCTAResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\ + " GEnergy&, GTime&, GObservation&)" +#define G_NPRED "GCTAResponse::npred(GSkyDir&, GEnergy&, GTime&,"\ " GObservation&)" #define G_IRF_EXTENDED "GCTAResponse::irf_extended(GInstDir&, GEnergy&,"\ " GTime&, GModelExtendedSource&, GEnergy&, GTime&, GObservation&)" #define G_IRF_DIFFUSE "GCTAResponse::irf_diffuse(GCTAInstDir&, GEnergy&,"\ " GTime&, GModelDiffuseSource&, GEnergy&, GTime&, GCTAObservation&)" -#define G_NPRED_PTSRC "GCTAResponse::npred_ptsrc("\ - "GModelPointSource&, GEnergy&, GTime&, GCTAObservation&)" -#define G_NPRED_EXTENDED "GCTAResponse::npred_extended("\ - "GModelExtendedSource&, GEnergy&, GTime&, GCTAObservation&)" -#define G_NPRED_DIFFUSE "GCTAResponse::npred_diffuse("\ - "GModelDiffuseSource&, GEnergy&, GTime&, GCTAObservation&)" #define G_READ "GCTAResponse::read_performance_table(std::string&)" /* __ Macros _____________________________________________________________ */ @@ -213,53 +207,70 @@ GCTAResponse* GCTAResponse::clone(void) const /***********************************************************************//** * @brief Return value of instrument response function * - * @param[in] event Event. - * @param[in] model Source model. + * @param[in] obsDir Observed photon direction. + * @param[in] obsEng Observed energy of photon. + * @param[in] obsTime Observed photon arrival time. + * @param[in] srcDir True photon arrival direction. * @param[in] srcEng True energy of photon. * @param[in] srcTime True photon arrival time. * @param[in] obs Observation. * - * @exception GCTAException::bad_observation_type - * Observation is not a CTA observation. * @exception GCTAException::bad_instdir_type - * Instrument direction is not a CTA instrument direction. - * - * @todo Throw exception if model pointer is invalid. + * Instrument direction is not a valid CTA instrument direction. ***************************************************************************/ -double GCTAResponse::irf(const GEvent& event, - const GModelSky& model, +double GCTAResponse::irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime, const GObservation& obs) const { + // Get pointer on CTA instrument direction + const GCTAInstDir* dir = dynamic_cast(&obsDir); + if (dir == NULL) + throw GCTAException::bad_instdir_type(G_IRF); + + // Get pointing direction zenith angle and azimuth + double zenith = 0.0; + double azimuth = 0.0; + + // Get radial offset and polar angles in camera + double theta = 0.0; + double phi = 0.0; + + // Get log10(E/TeV) of true photon energy. + double srcLogEng = srcEng.log10TeV(); + + // Determine angular separation between true and measured photon + // direction in radians + double delta = dir->dist(srcDir); + + // Get maximum angular separation for which PSF is significant + double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng); + // Initialise IRF value double irf = 0.0; - // Get pointer on CTA observation - const GCTAObservation* ctaobs = dynamic_cast(&obs); - if (ctaobs == NULL) - throw GCTAException::bad_observation_type(G_IRF); + // Compute only if we're sufficiently close to PSF + if (delta <= delta_max) { - // Get pointer on CTA instrument direction - const GCTAInstDir* ctadir = dynamic_cast(&event.dir()); - if (ctadir == NULL) - throw GCTAException::bad_instdir_type(G_IRF); + // Get point source IRF + irf = aeff(theta, phi, zenith, azimuth, srcLogEng); + irf *= psf(delta, theta, phi, zenith, azimuth, srcLogEng); - // Get model pointers - const GModelPointSource* ptsrc = dynamic_cast(&model); - const GModelExtendedSource* extsrc = dynamic_cast(&model); - const GModelDiffuseSource* difsrc = dynamic_cast(&model); - - // Call model dependent method - if (ptsrc != NULL) - irf = irf_ptsrc(*ctadir, event.energy(), event.time(), - *ptsrc, srcEng, srcTime, *ctaobs); - else if (extsrc != NULL) - irf = irf_extended(*ctadir, event.energy(), event.time(), - *extsrc, srcEng, srcTime, *ctaobs); - else if (difsrc != NULL) - irf = irf_diffuse(*ctadir, event.energy(), event.time(), - *difsrc, srcEng, srcTime, *ctaobs); + // Multiply-in energy dispersion + if (hasedisp()) { + + // Get log10(E/TeV) of measured photon energy. + double obsLogEng = obsEng.log10TeV(); + + // Multiply-in energy dispersion + irf *= edisp(obsLogEng, theta, phi, zenith, azimuth, srcLogEng); + + } // endif: energy dispersion was available + + } // endif: we were sufficiently close to PSF // Return IRF value return irf; @@ -267,45 +278,62 @@ double GCTAResponse::irf(const GEvent& event, /***********************************************************************//** - * @brief Return spatial integral of instrument response function + * @brief Return spatial integral of point spread function * - * @param[in] model Source model. + * @param[in] srcDir True photon arrival direction. * @param[in] srcEng True energy of photon. * @param[in] srcTime True photon arrival time. * @param[in] obs Observation. * * @exception GCTAException::bad_observation_type * Observation is not a CTA observation. - * - * @todo Throw exception if model pointer is invalid. + * @exception GException::no_list + * Observation does not contain a valid CTA event list. ***************************************************************************/ -double GCTAResponse::npred(const GModelSky& model, +double GCTAResponse::npred(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime, const GObservation& obs) const { - // Initialise Npred value - double npred = 0.0; - // Get pointer on CTA observation const GCTAObservation* ctaobs = dynamic_cast(&obs); if (ctaobs == NULL) throw GCTAException::bad_observation_type(G_NPRED); - // Get model pointers - const GModelPointSource* ptsrc = dynamic_cast(&model); - const GModelExtendedSource* extsrc = dynamic_cast(&model); - const GModelDiffuseSource* difsrc = dynamic_cast(&model); + // Get pointer on CTA events list + const GCTAEventList* events = dynamic_cast(ctaobs->events()); + if (events == NULL) + throw GException::no_list(G_NPRED); - // Call model dependent method - if (ptsrc != NULL) - npred = npred_ptsrc(*ptsrc, srcEng, srcTime, *ctaobs); - else if (extsrc != NULL) - npred = npred_extended(*extsrc, srcEng, srcTime, *ctaobs); - else if (difsrc != NULL) - npred = npred_diffuse(*difsrc, srcEng, srcTime, *ctaobs); + // Get pointer on CTA pointing + const GCTAPointing *pnt = ctaobs->pointing(srcTime); - // Return response value + // Get pointing direction zenith angle and azimuth + // @todo Not yet implemented (and not yet needed) + double zenith = 0.0; + double azimuth = 0.0; + + // Get radial offset and polar angles in camera + // @todo Not yet implemented (and not yet needed) + double theta = 0.0; + double phi = 0.0; + + // Get log10(E/TeV) of true photon energy. + double srcLogEng = srcEng.log10TeV(); + + // Get IRF components + double npred = aeff(theta, phi, zenith, azimuth, srcLogEng); + npred *= npsf(srcDir, srcLogEng, srcTime, *pnt, events->roi()); + + // Multiply-in energy dispersion + if (hasedisp()) { + + // Multiply-in energy dispersion + npred *= nedisp(srcDir, srcEng, srcTime, *pnt, events->ebounds()); + + } // endif: had energy dispersion + + // Return Npred return npred; } @@ -470,81 +498,6 @@ std::string GCTAResponse::print(void) const ==========================================================================*/ /***********************************************************************//** - * @brief Return value of point source instrument response function - * - * @param[in] obsDir Observed photon direction. - * @param[in] obsEng Observed energy of photon. - * @param[in] obsTime Observed photon arrival time. - * @param[in] model Point source model. - * @param[in] srcEng True energy of photon. - * @param[in] srcTime True photon arrival time. - * @param[in] obs CTA Observation. - * - * @todo The telescope zenith and azimuth angles as well as the radial offset - * and polar angles in the camera are not yet computed correctly (as - * the actual IRF implement does not need these values). - ***************************************************************************/ -double GCTAResponse::irf_ptsrc(const GCTAInstDir& obsDir, - const GEnergy& obsEng, - const GTime& obsTime, - const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const -{ - // Get point source location - GSkyDir srcDir = model.dir(); - - // Get CTA pointing - const GCTAPointing* pnt = obs.pointing(srcTime); - - // Get pointing direction zenith angle and azimuth - double zenith = 0.0; - double azimuth = 0.0; - - // Get radial offset and polar angles in camera - double theta = 0.0; - double phi = 0.0; - - // Get log10(E/TeV) of true photon energy. - double srcLogEng = srcEng.log10TeV(); - - // Determine angular separation between true and measured photon - // direction in radians - double delta = obsDir.dist(srcDir); - - // Get maximum angular separation for which PSF is significant - double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng); - - // Initialise IRF value - double irf = 0.0; - - // Compute only if we're sufficiently close to PSF - if (delta <= delta_max) { - - // Get point source IRF - irf = aeff(theta, phi, zenith, azimuth, srcLogEng); - irf *= psf(delta, theta, phi, zenith, azimuth, srcLogEng); - - // Multiply-in energy dispersion - if (hasedisp()) { - - // Get log10(E/TeV) of measured photon energy. - double obsLogEng = obsEng.log10TeV(); - - // Multiply-in energy dispersion - irf *= edisp(obsLogEng, theta, phi, zenith, azimuth, srcLogEng); - - } // endif: energy dispersion was available - - } // endif: we were sufficiently close to PSF - - // Return IRF value - return irf; -} - - -/***********************************************************************//** * @brief Return value of extended source instrument response function * * @param[in] obsDir Observed photon direction. @@ -553,7 +506,12 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir& obsDir, * @param[in] model Extended source model. * @param[in] srcEng True energy of photon. * @param[in] srcTime True photon arrival time. - * @param[in] obs CTA Observation. + * @param[in] obs Observation. + * + * @exception GCTAException::bad_observation_type + * Observation is not a valid CTA observation + * @exception GCTAException::bad_instdir_type + * Instrument direction is not a valid CTA instrument direction * * @todo The telescope zenith and azimuth angles as well as the radial offset * and polar angles in the camera are not yet computed correctly (as @@ -562,16 +520,26 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir& obsDir, * @todo The theta limits could still improved as only the theta range is * relevant for which the source and the PSF circle overlap. ***************************************************************************/ -double GCTAResponse::irf_extended(const GCTAInstDir& obsDir, +double GCTAResponse::irf_extended(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime, const GModelExtendedSource& model, const GEnergy& srcEng, const GTime& srcTime, - const GCTAObservation& obs) const + const GObservation& obs) const { - // Get CTA pointing - const GCTAPointing* pnt = obs.pointing(srcTime); + // Get pointer on CTA observation + const GCTAObservation* ctaobs = dynamic_cast(&obs); + if (ctaobs == NULL) + throw GCTAException::bad_observation_type(G_IRF_EXTENDED); + + // Get pointer on CTA instrument direction + const GCTAInstDir* dir = dynamic_cast(&obsDir); + if (dir == NULL) + throw GCTAException::bad_instdir_type(G_IRF_EXTENDED); + + // Get pointer on CTA pointing + const GCTAPointing *pnt = ctaobs->pointing(srcTime); // Get pointing direction zenith angle and azimuth double zenith = 0.0; @@ -589,7 +557,7 @@ double GCTAResponse::irf_extended(const GCTAInstDir& obsDir, double sigma = psf_dummy_sigma(srcLogEng); // Determine camera direction of measured photon - GCTADir obsCam(obsDir, *pnt); + GCTADir obsCam(*dir, *pnt); // Determine camera direction of source model centre GCTADir srcCam(model.dir(), *pnt); @@ -667,13 +635,13 @@ double GCTAResponse::irf_extended(const GCTAInstDir& obsDir, * * @todo Implement method. ***************************************************************************/ -double GCTAResponse::irf_diffuse(const GCTAInstDir& obsDir, +double GCTAResponse::irf_diffuse(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime, const GModelDiffuseSource& model, const GEnergy& srcEng, const GTime& srcTime, - const GCTAObservation& obs) const + const GObservation& obs) const { // Feature not yet implemented throw GException::feature_not_implemented(G_IRF_DIFFUSE, @@ -684,113 +652,6 @@ double GCTAResponse::irf_diffuse(const GCTAInstDir& obsDir, } -/***********************************************************************//** - * @brief Return spatial integral of point source model - * - * @param[in] model Point source model. - * @param[in] srcEng True energy of photon. - * @param[in] srcTime True photon arrival time. - * @param[in] obs Observation. - * - * @exception GException::no_list - * Observation does not contain an event list. - ***************************************************************************/ -double GCTAResponse::npred_ptsrc(const GModelPointSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const -{ - // Get point source location - GSkyDir srcDir = model.dir(); - - // Get pointer on CTA events list - const GCTAEventList* ptr = dynamic_cast(obs.events()); - if (ptr == NULL) - throw GException::no_list(G_NPRED_PTSRC); - - // Get pointing direction zenith angle and azimuth - const GCTAPointing *pnt = obs.pointing(srcTime); - double zenith = 0.0; - double azimuth = 0.0; - - // Get radial offset and polar angles in camera - double theta = 0.0; - double phi = 0.0; - - // Get log10(E/TeV) of true photon energy. - double srcLogEng = srcEng.log10TeV(); - - // Get IRF components - double nirf = aeff(theta, phi, zenith, azimuth, srcLogEng); - nirf *= npsf(srcDir, srcLogEng, srcTime, *pnt, ptr->roi()); - - // Multiply-in energy dispersion - if (hasedisp()) { - - // Multiply-in energy dispersion - nirf *= nedisp(srcDir, srcEng, srcTime, *pnt, ptr->ebounds()); - - } // endif: had energy dispersion - - // Return integrated IRF value - return nirf; -} - - -/***********************************************************************//** - * @brief Return spatial integral of extended source model - * - * @param[in] model Extended source model. - * @param[in] srcEng True energy of photon. - * @param[in] srcTime True photon arrival time. - * @param[in] obs Observation. - * - * @exception GException::feature_not_implemented - * Method not yet implemented. - * - * @todo Implement method. - ***************************************************************************/ -double GCTAResponse::npred_extended(const GModelExtendedSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const -{ - // Feature not yet implemented - throw GException::feature_not_implemented(G_NPRED_EXTENDED, - "Method for extended source not yet implemented."); - - // Return integrated IRF value - return 0.0; -} - - -/***********************************************************************//** - * @brief Return spatial integral of diffuse source model - * - * @param[in] model Diffuse source model. - * @param[in] srcEng True energy of photon. - * @param[in] srcTime True photon arrival time. - * @param[in] obs Observation. - * - * @exception GException::feature_not_implemented - * Method not yet implemented. - * - * @todo Implement method. - ***************************************************************************/ -double GCTAResponse::npred_diffuse(const GModelDiffuseSource& model, - const GEnergy& srcEng, - const GTime& srcTime, - const GCTAObservation& obs) const -{ - // Feature not yet implemented - throw GException::feature_not_implemented(G_NPRED_DIFFUSE, - "Method for extended source not yet implemented."); - - // Return integrated IRF value - return 0.0; -} - - /*========================================================================== = = = Low-level CTA response methods = diff --git a/inst/lat/include/GLATResponse.hpp b/inst/lat/include/GLATResponse.hpp index 57ac373..8d0b96f 100644 --- a/inst/lat/include/GLATResponse.hpp +++ b/inst/lat/include/GLATResponse.hpp @@ -32,8 +32,6 @@ #include "GModel.hpp" #include "GObservation.hpp" #include "GResponse.hpp" -#include "GEnergy.hpp" -#include "GTime.hpp" /***********************************************************************//** @@ -52,18 +50,28 @@ public: // Operators GLATResponse& operator= (const GLATResponse & rsp); - // Implement pure virtual base class methods - void clear(void); - GLATResponse* clone(void) const; - bool hasedisp(void) const { return false; } - bool hastdisp(void) const { return false; } - double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, - const GObservation& obs) const; - double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, - const GObservation& obs) const; - std::string print(void) const; + // Implemented pure virtual methods + virtual void clear(void); + virtual GLATResponse* clone(void) const; + virtual bool hasedisp(void) const { return false; } + virtual bool hastdisp(void) const { return false; } + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual std::string print(void) const; + + // Implemented virtual methods + virtual double irf(const GEvent& event, const GModelSky& model, + const GEnergy& srcEng, const GTime& srcTime, + const GObservation& obs) const; // Other Methods void caldb(const std::string& caldb); diff --git a/inst/lat/pyext/GLATResponse.i b/inst/lat/pyext/GLATResponse.i index c528172..8cd73f8 100644 --- a/inst/lat/pyext/GLATResponse.i +++ b/inst/lat/pyext/GLATResponse.i @@ -1,5 +1,5 @@ /*************************************************************************** - * GLATResponse.i - Fermi LAT Response class * + * GLATResponse.i - Fermi-LAT response class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2011 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * @@ -12,7 +12,7 @@ ***************************************************************************/ /** * @file GLATResponse.i - * @brief Fermi LAT Response class python bindings + * @brief Fermi-LAT Response class Python interface definition * @author J. Knodlseder */ %{ @@ -24,26 +24,37 @@ /***********************************************************************//** * @class GLATResponse * - * @brief Fermi LAT Response class python bindings + * @brief Fermi-LAT response class ***************************************************************************/ class GLATResponse : public GResponse { + public: // Constructors and destructors GLATResponse(void); GLATResponse(const GLATResponse& rsp); virtual ~GLATResponse(void); - // Implement pure virtual base class methods - void clear(void); - GLATResponse* clone(void) const; - bool hasedisp(void) const; - bool hastdisp(void) const; - double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, - const GObservation& obs) const; - double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, - const GObservation& obs) const; + // Implemented pure virtual methods + virtual void clear(void); + virtual GLATResponse* clone(void) const; + virtual bool hasedisp(void) const; + virtual bool hastdisp(void) const; + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + + // Implemented virtual methods + virtual double irf(const GEvent& event, const GModelSky& model, + const GEnergy& srcEng, const GTime& srcTime, + const GObservation& obs) const; // Other Methods void caldb(const std::string& caldb); diff --git a/inst/lat/src/GLATResponse.cpp b/inst/lat/src/GLATResponse.cpp index a7ddfcf..dc95ce1 100644 --- a/inst/lat/src/GLATResponse.cpp +++ b/inst/lat/src/GLATResponse.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GLATResponse.cpp - Fermi LAT Response class * + * GLATResponse.cpp - Fermi-LAT response class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2011 by Jurgen Knodlseder * * ----------------------------------------------------------------------- * @@ -12,7 +12,7 @@ ***************************************************************************/ /** * @file GLATResponse.cpp - * @brief Fermi LAT Response class implementation. + * @brief Fermi-LAT response class implementation * @author J. Knodlseder */ @@ -23,6 +23,7 @@ #include // access() function #include // std::getenv() function #include +#include "GLATInstDir.hpp" #include "GLATResponse.hpp" #include "GLATObservation.hpp" #include "GLATEventAtom.hpp" @@ -37,6 +38,8 @@ /* __ Method name definitions ____________________________________________ */ #define G_CALDB "GLATResponse::caldb(std::string&)" #define G_LOAD "GLATResponse::load(std::string&)" +#define G_IRF "GLATResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\ + " GEnergy&, GTime&, GObservation&)" #define G_AEFF "GLATResponse::aeff(int&)" #define G_PSF "GLATResponse::psf(int&)" #define G_EDISP "GLATResponse::edisp(int&)" @@ -44,7 +47,7 @@ "GTime&, GObservation&)" #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\ "GTime&, GObservation&)" -#define G_NPRED "GLATResponse::npred(GModel&, GEnergy&, GTime&,"\ +#define G_NPRED "GLATResponse::npred(GSkyDir&, GEnergy&, GTime&,"\ " GObservation&)" /* __ Macros _____________________________________________________________ */ @@ -69,7 +72,7 @@ ***************************************************************************/ GLATResponse::GLATResponse(void) : GResponse() { - // Initialise class members for clean destruction + // Initialise members init_members(); // Return @@ -84,7 +87,7 @@ GLATResponse::GLATResponse(void) : GResponse() ***************************************************************************/ GLATResponse::GLATResponse(const GLATResponse& rsp) : GResponse(rsp) { - // Initialise class members for clean destruction + // Initialise members init_members(); // Copy members @@ -130,7 +133,7 @@ GLATResponse& GLATResponse::operator= (const GLATResponse& rsp) // Free members free_members(); - // Initialise private members for clean destruction + // Initialise members init_members(); // Copy members @@ -198,7 +201,7 @@ void GLATResponse::caldb(const std::string& caldb) // Check if calibration database directory is accessible if (access(caldb.c_str(), R_OK) != 0) throw GException::caldb_not_found(G_CALDB, caldb); - + // Store the path to the calibration database m_caldb = caldb; @@ -208,6 +211,78 @@ void GLATResponse::caldb(const std::string& caldb) /***********************************************************************//** + * @brief Return value of point source IRF + * + * @param[in] obsDir Measured photon arrival direction. + * @param[in] obsEng Measured photon energy. + * @param[in] obsTime Measured photon arrival time. + * @param[in] srcDir True photon arrival direction. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @exception GLATException::bad_instdir_type + * Instrument direction is not a valid LAT instrument direction. + * + * @todo The IRF value is not devided by ontime of the event, but it is + * already time integrated. + ***************************************************************************/ +double GLATResponse::irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Get LAT instrument direction + const GLATInstDir* dir = dynamic_cast(&obsDir); + if (dir == NULL) + throw GLATException::bad_instdir_type(G_IRF); + + // Search for mean PSF + int ipsf = -1; + for (int i = 0; i < m_ptsrc.size(); ++i) { + if (m_ptsrc[i]->dir() == srcDir) { + ipsf = i; + break; + } + } + + // If mean PSF has not been found then create it now + if (ipsf == -1) { + + // Allocate new mean PSF + GLATMeanPsf* psf = new GLATMeanPsf(srcDir, static_cast(obs)); + + // Set source name + std::string name = "SRC("+str(srcDir.ra_deg())+","+str(srcDir.dec_deg())+")"; + psf->name(name); + + // Push mean PSF on stack + ((GLATResponse*)this)->m_ptsrc.push_back(psf); + + // Set index of mean PSF + ipsf = m_ptsrc.size()-1; + + // Debug option: dump mean PSF + #if G_DUMP_MEAN_PSF + std::cout << "Added new mean PSF \""+name+"\"" << std::endl; + std::cout << *psf << std::endl; + #endif + + } // endif: created new mean PSF + + // Get IRF value + double offset = dir->dist_deg(srcDir); + double irf = (*m_ptsrc[ipsf])(offset, srcEng.log10MeV()); + + // Return IRF value + return irf; +} + + +/***********************************************************************//** * @brief Return value of model instrument response function * * @param[in] event Event. @@ -265,7 +340,7 @@ double GLATResponse::irf(const GLATEventAtom& event, const GModelSky& model, /***********************************************************************//** * @brief Return value of model IRF for event bin - * + * * @param[in] event Event bin. * @param[in] model Source model. * @param[in] srcEng True energy of photon. @@ -404,14 +479,14 @@ double GLATResponse::irf(const GLATEventBin& event, const GModelSky& model, /***********************************************************************//** * @brief Return integral of instrument response function. * - * @param[in] model Source model. + * @param[in] srcDir True photon arrival direction. * @param[in] srcEng True energy of photon. * @param[in] srcTime True photon arrival time. * @param[in] obs Observation. * * @todo Not yet implemented. ***************************************************************************/ -double GLATResponse::npred(const GModelSky& model, const GEnergy& srcEng, +double GLATResponse::npred(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime, const GObservation& obs) const { diff --git a/inst/mwl/include/GMWLResponse.hpp b/inst/mwl/include/GMWLResponse.hpp index a6a5977..9f46bd3 100644 --- a/inst/mwl/include/GMWLResponse.hpp +++ b/inst/mwl/include/GMWLResponse.hpp @@ -12,7 +12,7 @@ ***************************************************************************/ /** * @file GMWLResponse.hpp - * @brief GMWLResponse class interface definition. + * @brief Multi-wavelength response class interface definition * @author J. Knodlseder */ @@ -26,11 +26,14 @@ /***********************************************************************//** * @class GMWLResponse * - * @brief Interface for the Multi-wavelength response class + * @brief Multi-wavelength response class * - * The Multi-wavelength response class is needed to evaluate the response - * to a source model. Since the multi-wavelength instrument classes work - * directly in photon space, the response is by definition unity. + * This class implements a dummy response class for multi-wavelength + * observations. Since the multi-wavelength instrument classes handles data + * that a provided in photon space, no instrument response is in fact needed. + * The dummy response implemented by this class provides a simple diagonal + * response matrix that allows integration of multi-wavelength observations + * using the standard instrument specific interface. ***************************************************************************/ class GMWLResponse : public GResponse { @@ -41,22 +44,25 @@ public: virtual ~GMWLResponse(void); // Operators - GMWLResponse& operator= (const GMWLResponse & rsp); + virtual GMWLResponse& operator= (const GMWLResponse & rsp); - // Reponse function computation methods - double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, - const GObservation& obs) const { return 1.0; } - double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, - const GObservation& obs) const { return 1.0; } - - // Pure virtual base class methods - void clear(void); - GMWLResponse* clone(void) const; - bool hasedisp(void) const { return false; } - bool hastdisp(void) const { return false; } - std::string print(void) const; + // Implemented pure virtual methods + virtual void clear(void); + virtual GMWLResponse* clone(void) const; + virtual bool hasedisp(void) const { return false; } + virtual bool hastdisp(void) const { return false; } + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const { return 1.0; } + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const { return 1.0; } + virtual std::string print(void) const; protected: // Protected methods diff --git a/pyext/GResponse.i b/pyext/GResponse.i index 19aa8e7..a6b948b 100644 --- a/pyext/GResponse.i +++ b/pyext/GResponse.i @@ -12,7 +12,7 @@ ***************************************************************************/ /** * @file GResponse.i - * @brief Abstract response base class python bindings + * @brief Abstract response base class Python interface definition * @author J. Knodlseder */ %{ @@ -25,9 +25,10 @@ /***********************************************************************//** * @class GResponse * - * @brief Abstract response base class python bindings + * @brief Abstract response base class ***************************************************************************/ class GResponse { + public: // Constructors and destructors GResponse(void); @@ -35,16 +36,65 @@ public: virtual ~GResponse(void); // Pure virtual methods - virtual void clear(void) = 0; - virtual GResponse* clone(void) const = 0; - virtual bool hasedisp(void) const = 0; - virtual bool hastdisp(void) const = 0; - virtual double irf(const GEvent& event, const GModelSky& model, - const GEnergy& srcEng, const GTime& srcTime, - const GObservation& obs) const = 0; - virtual double npred(const GModelSky& model, const GEnergy& srcEng, - const GTime& srcTime, - const GObservation& obs) const = 0; + virtual void clear(void) = 0; + virtual GResponse* clone(void) const = 0; + virtual bool hasedisp(void) const = 0; + virtual bool hastdisp(void) const = 0; + virtual double irf(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const = 0; + virtual double npred(const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const = 0; + + // Virtual methods + virtual double irf(const GEvent& event, + const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_ptsrc(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_extended(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double irf_diffuse(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred(const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_ptsrc(const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_extended(const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; + virtual double npred_diffuse(const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const; }; diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index d368e60..b48f6e0 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -18,16 +18,29 @@ #include #include // access() function #include "GResponse.hpp" +#include "GObservation.hpp" +#include "GIntegral.hpp" +#include "GVector.hpp" +#include "GSkyDir.hpp" #include "GException.hpp" +#include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_IRF_EXTENDED "GResponse::irf_extended(GInstDir&, GEnergy&,"\ + " GTime&, GModelExtendedSource&, GEnergy&, GTime&, GObservation&)" +#define G_IRF_DIFFUSE "GResponse::irf_diffuse(GInstDir&, GEnergy&,"\ + " GTime&, GModelDiffuseSource&, GEnergy&, GTime&, GObservation&)" +#define G_NPRED_EXTENDED "GResponse::npred_extended(GModelExtendedSource&,"\ + " GEnergy&, GTime&, GObservation&)" +#define G_NPRED_DIFFUSE "GResponse::npred_diffuse(GModelDiffuseSource&,"\ + " GEnergy&, GTime&, GObservation&)" /* __ Macros _____________________________________________________________ */ /* __ Coding definitions _________________________________________________ */ /* __ Debug definitions __________________________________________________ */ - +#define G_DEBUG_NPRED_EXTENDED 0 //!< Debug npred_extended /*========================================================================== @@ -118,10 +131,286 @@ GResponse& GResponse::operator= (const GResponse& rsp) = = ==========================================================================*/ +/***********************************************************************//** + * @brief Return value of instrument response function + * + * @param[in] event Event. + * @param[in] model Source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @todo Throw exception if model pointer is invalid. + ***************************************************************************/ +double GResponse::irf(const GEvent& event, + const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Initialise IRF value + double irf = 0.0; + + // Get model pointers + const GModelPointSource* ptsrc = dynamic_cast(&model); + const GModelExtendedSource* extsrc = dynamic_cast(&model); + const GModelDiffuseSource* difsrc = dynamic_cast(&model); + + // Call model dependent method + if (ptsrc != NULL) + irf = irf_ptsrc(event.dir(), event.energy(), event.time(), + *ptsrc, srcEng, srcTime, obs); + else if (extsrc != NULL) + irf = irf_extended(event.dir(), event.energy(), event.time(), + *extsrc, srcEng, srcTime, obs); + else if (difsrc != NULL) + irf = irf_diffuse(event.dir(), event.energy(), event.time(), + *difsrc, srcEng, srcTime, obs); + + // Return IRF value + return irf; +} + + +/***********************************************************************//** + * @brief Return value of point source instrument response function + * + * @param[in] obsDir Observed photon direction. + * @param[in] obsEng Observed energy of photon. + * @param[in] obsTime Observed photon arrival time. + * @param[in] model Point source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + ***************************************************************************/ +double GResponse::irf_ptsrc(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Get point source location + GSkyDir srcDir = model.dir(); + + // Compute IRF + double irf = this->irf(obsDir, obsEng, obsTime, srcDir, srcEng, srcTime, obs); + + // Return IRF + return irf; +} + + +/***********************************************************************//** + * @brief Return value of extended source instrument response function + * + * @param[in] obsDir Observed photon direction. + * @param[in] obsEng Observed energy of photon. + * @param[in] obsTime Observed photon arrival time. + * @param[in] model Exteded source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @todo Generic method is not yet implemented. + ***************************************************************************/ +double GResponse::irf_extended(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Feature not yet implemented + throw GException::feature_not_implemented(G_IRF_EXTENDED, + "Extended IRF not yet implemented."); + + // Return IRF + return 0.0; +} + + +/***********************************************************************//** + * @brief Return value of diffuse source instrument response function + * + * @param[in] obsDir Observed photon direction. + * @param[in] obsEng Observed energy of photon. + * @param[in] obsTime Observed photon arrival time. + * @param[in] model Diffuse source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @todo Generic method is not yet implemented. + ***************************************************************************/ +double GResponse::irf_diffuse(const GInstDir& obsDir, + const GEnergy& obsEng, + const GTime& obsTime, + const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Feature not yet implemented + throw GException::feature_not_implemented(G_IRF_DIFFUSE, + "Diffuse IRF not yet implemented."); + + // Return IRF + return 0.0; +} + + +/***********************************************************************//** + * @brief Return spatial integral of instrument response function + * + * @param[in] model Source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @todo Throw exception if model pointer is invalid. + ***************************************************************************/ +double GResponse::npred(const GModelSky& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Initialise Npred value + double npred = 0.0; + + // Get model pointers + const GModelPointSource* ptsrc = dynamic_cast(&model); + const GModelExtendedSource* extsrc = dynamic_cast(&model); + const GModelDiffuseSource* difsrc = dynamic_cast(&model); + + // Call model dependent method + if (ptsrc != NULL) + npred = npred_ptsrc(*ptsrc, srcEng, srcTime, obs); + else if (extsrc != NULL) + npred = npred_extended(*extsrc, srcEng, srcTime, obs); + else if (difsrc != NULL) + npred = npred_diffuse(*difsrc, srcEng, srcTime, obs); + + // Return response value + return npred; +} + + +/***********************************************************************//** + * @brief Return spatial integral of point source model + * + * @param[in] model Point source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + ***************************************************************************/ +double GResponse::npred_ptsrc(const GModelPointSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Get point source location + GSkyDir srcDir = model.dir(); + + // Compute Npred + double npred = this->npred(srcDir, srcEng, srcTime, obs); + + // Return Npred + return npred; +} + + +/***********************************************************************//** + * @brief Return spatial integral of extended source model + * + * @param[in] model Extended source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + ***************************************************************************/ +double GResponse::npred_extended(const GModelExtendedSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Initialise Npred value + double npred = 0.0; + + // Compute rotation matrix to convert from native coordinates given by + // (theta,phi) into celestial coordinates. + GMatrix ry; + GMatrix rz; + GMatrix rot; + ry.eulery(model.radial()->dec() - 90.0); + rz.eulerz(-model.radial()->ra()); + rot = transpose(ry * rz); + + // Set offset angle integration range + double theta_min = 0.0; + double theta_max = model.radial()->theta_max(); + + // Perform offset angle integration if interval is valid + if (theta_max > theta_min) { + + // Setup integration kernel + GResponse::npred_kern_theta integrand(this, + model.radial(), + &srcEng, + &srcTime, + &obs, + &rot); + + // Integrate over theta + GIntegral integral(&integrand); + npred = integral.romb(theta_min, theta_max); + + // Compile option: Show integration results + #if G_DEBUG_NPRED_EXTENDED + std::cout << "npred_extended:"; + std::cout << " theta_min=" << theta_min; + std::cout << " theta_max=" << theta_max; + std::cout << " npred=" << npred << std::endl; + #endif + + } // endif: offset angle range was valid + + // Return Npred + return npred; +} + + +/***********************************************************************//** + * @brief Return spatial integral of diffuse source model + * + * @param[in] model Diffuse source model. + * @param[in] srcEng True energy of photon. + * @param[in] srcTime True photon arrival time. + * @param[in] obs Observation. + * + * @exception GException::feature_not_implemented + * Method not yet implemented. + * + * @todo Implement method. + ***************************************************************************/ +double GResponse::npred_diffuse(const GModelDiffuseSource& model, + const GEnergy& srcEng, + const GTime& srcTime, + const GObservation& obs) const +{ + // Feature not yet implemented + throw GException::feature_not_implemented(G_NPRED_DIFFUSE, + "Method for extended source not yet implemented."); + + // Return Npred + return 0.0; +} + /*========================================================================== = = - = Private methods = + = Protected methods = = = ==========================================================================*/ @@ -157,6 +446,67 @@ void GResponse::free_members(void) } +/***********************************************************************//** + * @brief Kernel for offset angle Npred integration + * + * @param[in] theta Radial model offset angle (radians). + ***************************************************************************/ +double GResponse::npred_kern_theta::eval(double theta) +{ + // Initialise Npred value + double npred = 0.0; + + // Get radial model value + double model = m_radial->eval(theta); + + // Compute sine of offset angle + double sin_theta = std::sin(theta); + + // Setup phi integration kernel + GResponse::npred_kern_phi integrand(m_rsp, + m_srcEng, + m_srcTime, + m_obs, + m_rot, + theta, + sin_theta); + + // Integrate over phi + GIntegral integral(&integrand); + npred = integral.romb(0.0, twopi) * sin_theta * model; + + // Return Npred + return npred; +} + + +/***********************************************************************//** + * @brief Kernel for azimuth angle Npred integration + * + * @param[in] phi Azimuth angle (radians). + ***************************************************************************/ +double GResponse::npred_kern_phi::eval(double phi) +{ + // Compute sky direction vector in native coordinates + double cos_phi = std::cos(phi); + double sin_phi = std::sin(phi); + GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta); + + // Rotate from native into celestial system + GVector cel = *m_rot * native; + + // Set sky direction + GSkyDir srcDir; + srcDir.celvector(cel); + + // Compute Npred for this sky direction + double npred = m_rsp->npred(srcDir, *m_srcEng, *m_srcTime, *m_obs); + + // Return Npred + return npred; +} + + /*========================================================================== = = = Friends =