diff --git a/inst/cta/include/GCTACubeBackground.hpp b/inst/cta/include/GCTACubeBackground.hpp index 65cca91..a7270b1 100644 --- a/inst/cta/include/GCTACubeBackground.hpp +++ b/inst/cta/include/GCTACubeBackground.hpp @@ -65,7 +65,7 @@ public: void clear(void); GCTACubeBackground* clone(void) const; std::string classname(void) const; - void fill(const GObservations& obs, GLog* log = 0); + void fill(const GObservations& obs, GLog* log = NULL); double integral(const double& logE) const; void read(const GFits& fits); void write(GFits& file) const; diff --git a/inst/cta/include/GCTACubeExposure.hpp b/inst/cta/include/GCTACubeExposure.hpp index 7e89728..ab36efa 100644 --- a/inst/cta/include/GCTACubeExposure.hpp +++ b/inst/cta/include/GCTACubeExposure.hpp @@ -74,7 +74,7 @@ public: GCTACubeExposure* clone(void) const; std::string classname(void) const; void set(const GCTAObservation& obs); - void fill(const GObservations& obs, GLog* log = 0); + void fill(const GObservations& obs, GLog* log = NULL); const GSkymap& cube(void) const; const GEbounds& ebounds(void) const; const GGti& gti(void) const; diff --git a/inst/cta/include/GCTACubePsf.hpp b/inst/cta/include/GCTACubePsf.hpp index 7d1b948..8fcc251 100644 --- a/inst/cta/include/GCTACubePsf.hpp +++ b/inst/cta/include/GCTACubePsf.hpp @@ -82,7 +82,7 @@ public: GCTACubePsf* clone(void) const; std::string classname(void) const; void set(const GCTAObservation& obs); - void fill(const GObservations& obs, GLog* log = 0); + void fill(const GObservations& obs, GLog* log = NULL); const GSkymap& map(void) const; const GEbounds& ebounds(void) const; const GNodeArray& deltas(void) const; diff --git a/inst/cta/include/GCTAObservation.hpp b/inst/cta/include/GCTAObservation.hpp index dd2ef4b..3cb7ea3 100644 --- a/inst/cta/include/GCTAObservation.hpp +++ b/inst/cta/include/GCTAObservation.hpp @@ -124,12 +124,12 @@ public: void deadc(const double& deadc); void eventfile(const std::string& filename); const std::string& eventfile(void) const; + const std::string& eventtype(void) const; void dispose_events(void); const double& lo_user_thres(void) const; const double& hi_user_thres(void) const; void n_tels(const int& tels); const int& n_tels(void) const; - const std::string& eventtype(void) const; protected: // Protected methods @@ -414,6 +414,18 @@ const std::string& GCTAObservation::eventfile(void) const /***********************************************************************//** + * @brief Return type of events + * + * @return Type of events. + ***************************************************************************/ +inline +const std::string& GCTAObservation::eventtype(void) const +{ + return m_eventtype; +} + + +/***********************************************************************//** * @brief Return user low energy threshold * * @return User low energy threshold. @@ -460,15 +472,4 @@ const int& GCTAObservation::n_tels(void) const return m_n_tels; } - -/***********************************************************************//** - * @brief Return type of events - * - * @return Type of events. - ***************************************************************************/ -inline -const std::string& GCTAObservation::eventtype(void) const -{ - return m_eventtype; -} #endif /* GCTAOBSERVATION_HPP */ diff --git a/inst/cta/pyext/GCTACubeBackground.i b/inst/cta/pyext/GCTACubeBackground.i index db43967..72e0a3b 100644 --- a/inst/cta/pyext/GCTACubeBackground.i +++ b/inst/cta/pyext/GCTACubeBackground.i @@ -52,7 +52,7 @@ public: void clear(void); GCTACubeBackground* clone(void) const; std::string classname(void) const; - void fill(const GObservations& obs, GLog* log=0); + void fill(const GObservations& obs, GLog* log = NULL); double integral(const double& logE) const; void read(const GFits& fits); void write(GFits& file) const; diff --git a/inst/cta/pyext/GCTACubeExposure.i b/inst/cta/pyext/GCTACubeExposure.i index 57f64cb..edbce61 100644 --- a/inst/cta/pyext/GCTACubeExposure.i +++ b/inst/cta/pyext/GCTACubeExposure.i @@ -65,7 +65,7 @@ public: GCTACubeExposure* clone(void) const; std::string classname(void) const; void set(const GCTAObservation& obs); - void fill(const GObservations& obs, GLog* log=0); + void fill(const GObservations& obs, GLog* log = NULL); const GSkymap& cube(void) const; const GEbounds& ebounds(void) const; const GGti& gti(void) const; diff --git a/inst/cta/pyext/GCTACubePsf.i b/inst/cta/pyext/GCTACubePsf.i index caf2455..de64a18 100644 --- a/inst/cta/pyext/GCTACubePsf.i +++ b/inst/cta/pyext/GCTACubePsf.i @@ -74,7 +74,7 @@ public: GCTACubePsf* clone(void) const; std::string classname(void) const; void set(const GCTAObservation& obs); - void fill(const GObservations& obs, GLog* log=0); + void fill(const GObservations& obs, GLog* log = NULL); const GSkymap& map(void) const; const GEbounds& ebounds(void) const; const GNodeArray& deltas(void) const; diff --git a/inst/cta/pyext/GCTAObservation.i b/inst/cta/pyext/GCTAObservation.i index 3a73485..cd4abef 100644 --- a/inst/cta/pyext/GCTAObservation.i +++ b/inst/cta/pyext/GCTAObservation.i @@ -96,12 +96,12 @@ public: void deadc(const double& deadc); void eventfile(const std::string& filename); const std::string& eventfile(void) const; + const std::string& eventtype(void) const; void dispose_events(void); const double& lo_user_thres(void) const; const double& hi_user_thres(void) const; void n_tels(const int& tels); const int& n_tels(void) const; - const std::string& eventtype(void) const; }; diff --git a/inst/cta/src/GCTACubeBackground.cpp b/inst/cta/src/GCTACubeBackground.cpp index 0fe30b8..b711428 100644 --- a/inst/cta/src/GCTACubeBackground.cpp +++ b/inst/cta/src/GCTACubeBackground.cpp @@ -264,6 +264,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const * @brief Fill background cube from observation container * * @param[in] obs Observation container. + * @param[in] log Pointer to logger (optional). * * @exception GException::invalid_value * No event list found in CTA observations. @@ -289,60 +290,63 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) // Get observation and continue only if it is a CTA observation const GCTAObservation* cta = dynamic_cast(obs[i]); - if (cta != NULL) { - - // Skip observation if we don't have an unbinned observation - if (cta->eventtype() != "EventList") { - - // Log that we skip the this observation - if (log != NULL) { - *log << "Warning: Skipping binned observation "; - *log << "\"" << cta->name() << "\"" <instrument()+" observation "; + *log << "\"" << obs[i]->name() << "\"" <eventtype() != "EventList") { + if (log != NULL) { + *log << "Warning: Skipping binned observation "; + *log << "\"" << cta->name() << "\"" <roi(); - - // Set GTI of actual observations as the GTI of the event cube - eventcube.gti(cta->gti()); + // Extract region of interest from CTA observation + GCTARoi roi = cta->roi(); - // Get observation livetime - double livetime = cta->livetime(); + // Set GTI of actual observations as the GTI of the event cube + eventcube.gti(cta->gti()); - // Loop over all bins in background cube - for (int i = 0; i < eventcube.size(); ++i) { + // Get observation livetime + double livetime = cta->livetime(); - // Get event bin - GCTAEventBin* bin = eventcube[i]; + // Loop over all bins in background cube + for (int i = 0; i < eventcube.size(); ++i) { - // Continue only if binned in contained in ROI - if (roi.contains(*bin)) { + // Get event bin + GCTAEventBin* bin = eventcube[i]; - // Compute model value for event bin. The model value is - // given in counts/MeV/s/sr. - double model = obs.models().eval(*bin, *cta); + // Continue only if binned in contained in ROI + if (roi.contains(*bin)) { - // Multiply by livetime to get the correct weighting for - // each observation. We divide by the total livetime later - // to get the background model in units of counts/MeV/s/sr. - model *= livetime; + // Compute model value for event bin. The model value is + // given in counts/MeV/s/sr. + double model = obs.models().eval(*bin, *cta); - // Add existing number of counts - model += bin->counts(); + // Multiply by livetime to get the correct weighting for + // each observation. We divide by the total livetime later + // to get the background model in units of counts/MeV/s/sr. + model *= livetime; - // Store cumulated value (units: counts/MeV/sr) - bin->counts(model); + // Add existing number of counts + model += bin->counts(); - } // endif: bin was contained in RoI + // Store cumulated value (units: counts/MeV/sr) + bin->counts(model); - } // endfor: looped over all bins + } // endif: bin was contained in RoI - // Accumulate livetime - total_livetime += livetime; + } // endfor: looped over all bins - } // endif: cta observation was vaild + // Accumulate livetime + total_livetime += livetime; } // endfor: looped over all observations diff --git a/inst/cta/src/GCTACubeExposure.cpp b/inst/cta/src/GCTACubeExposure.cpp index d2dfb4e..e1aa388 100644 --- a/inst/cta/src/GCTACubeExposure.cpp +++ b/inst/cta/src/GCTACubeExposure.cpp @@ -364,6 +364,7 @@ void GCTACubeExposure::set(const GCTAObservation& obs) * @brief Fill exposure cube from observation container * * @param[in] obs Observation container. + * @param[in] log Pointer to logger (optional). * * @exception GException::invalid_value * No event list found in CTA observations. @@ -384,66 +385,74 @@ void GCTACubeExposure::fill(const GObservations& obs, GLog* log) // Get observation and continue only if it is a CTA observation const GCTAObservation* cta = dynamic_cast(obs[i]); - if (cta != NULL) { - // Skip observation if we don't have an unbinned observation - if (cta->eventtype() != "EventList") { - - // Log that we skip the this observation - if (log != NULL) { - *log << "Warning: Skipping binned observation "; - *log << "\"" << cta->name() << "\"" <instrument()+" observation "; + *log << "\"" << obs[i]->name() << "\"" <roi(); + // Skip observation if we don't have an unbinned observation + if (cta->eventtype() != "EventList") { + if (log != NULL) { + *log << "Warning: Skipping binned observation "; + *log << "\"" << cta->name() << "\"" <(cta->response()); - const GSkyDir& pnt = cta->pointing().dir(); + // Extract region of interest from CTA observation + GCTARoi roi = cta->roi(); - // Continue only if response is valid - if (rsp != NULL) { + // Get references on CTA response and pointing direction + const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); + const GSkyDir& pnt = cta->pointing().dir(); + + // Skip observation if we don't have an unbinned observation + if (rsp == NULL) { + if (log != NULL) { + *log << "Warning: Observation \"" << cta->name(); + *log << "\" contains no IRF response." <aeff(theta, 0.0, 0.0, 0.0, logE) * - cta->livetime(); + // Get logE/TeV + double logE = m_ebounds.elogmean(iebin).log10TeV(); - } // endfor: looped over energy bins + // Add to exposure cube (effective area * livetime) + m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + cta->livetime(); - } // endif: pixel within RoI + } // endfor: looped over energy bins - } // endfor: looped over all pixels + } // endif: pixel within RoI - // Append GTIs and increment livetime - m_gti.extend(cta->gti()); - m_livetime += cta->livetime(); - - } // endif: response was valid - - } // endif: observation was a CTA observation + } // endfor: looped over all pixels + // Append GTIs and increment livetime + m_gti.extend(cta->gti()); + m_livetime += cta->livetime(); + } // endfor: looped over observations // Return diff --git a/inst/cta/src/GCTACubePsf.cpp b/inst/cta/src/GCTACubePsf.cpp index 1eb03e6..dcaa9b4 100644 --- a/inst/cta/src/GCTACubePsf.cpp +++ b/inst/cta/src/GCTACubePsf.cpp @@ -450,79 +450,87 @@ void GCTACubePsf::fill(const GObservations& obs, GLog* log) // Get observation and continue only if it is a CTA observation const GCTAObservation* cta = dynamic_cast(obs[i]); - if (cta != NULL) { - // Skip observation if we don't have an unbinned observation - if (cta->eventtype() != "EventList") { - - // Log that we skip the this observation - if (log != NULL) { - *log << "Warning: Skipping binned observation "; - *log << "\"" << cta->name() << "\"" << std::endl; - } - continue; + // Skip observation if it's not CTA + if (cta == NULL) { + if (log != NULL) { + *log << "Warning: Skipping "+obs[i]->instrument()+" observation "; + *log << "\"" << obs[i]->name() << "\"" <roi(); + // Skip observation if we don't have an unbinned observation + if (cta->eventtype() != "EventList") { + if (log != NULL) { + *log << "Warning: Skipping binned observation "; + *log << "\"" << cta->name() << "\"" <(cta->response()); - const GSkyDir& pnt = cta->pointing().dir(); + // Extract region of interest from CTA observation + GCTARoi roi = cta->roi(); - // Continue only if response is valid - if (rsp != NULL) { + // Get references on CTA response and pointing direction + const GCTAResponseIrf* rsp = dynamic_cast(cta->response()); + const GSkyDir& pnt = cta->pointing().dir(); + + // Skip observation if we don't have an unbinned observation + if (rsp == NULL) { + if (log != NULL) { + *log << "Warning: Observation \"" << cta->name(); + *log << "\" contains no IRF response." <aeff(theta, 0.0, 0.0, 0.0, logE) * - cta->livetime(); + // Loop over all energy bins + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - // Accumulate weights - exposure(pixel, iebin) += weight; + // Get logE/TeV + double logE = m_ebounds.elogmean(iebin).log10TeV(); - // Loop over delta values - for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { + // Compute exposure weight + double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) * + cta->livetime(); - // Compute delta in radians - double delta = m_deltas[idelta] * gammalib::deg2rad; + // Accumulate weights + exposure(pixel, iebin) += weight; - // Set map index - int imap = offset(idelta, iebin); + // Loop over delta values + for (int idelta = 0; idelta < m_deltas.size(); ++idelta) { - // Add on PSF cube - m_cube(pixel, imap) += - rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; + // Compute delta in radians + double delta = m_deltas[idelta] * gammalib::deg2rad; - } // endfor: looped over delta bins + // Set map index + int imap = offset(idelta, iebin); - } // endfor: looped over energy bins + // Add on PSF cube + m_cube(pixel, imap) += + rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight; - } // endif: pixel was within RoI + } // endfor: looped over delta bins - } // endfor: looped over all pixels + } // endfor: looped over energy bins - } // endif: response was valid + } // endif: pixel was within RoI - } // endif: observation was a CTA observation + } // endfor: looped over all pixels } // endfor: looped over observations