Make sure that simulated events are within ROI. Add also an ROI margin of 0.5 de…

…g so that also sources outside but at the edge of the ROI will be simulated.

Make sure that simulated events are within ROI. Add also an ROI margin of 0.5 de…
…g so that also sources outside but at the edge of the ROI will be simulated.
Jürgen Knödlseder
1 parent 6e39d7c7
Showing 5 changed files with 45 additions and 27 deletions
ChangeLog
NEWS
src/ctobssim/ctobssim.cpp
src/ctobssim/ctobssim.hpp
src/ctobssim/ctobssim.par
-2013-02-21  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2013-04-09  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.7.0 released
         ========================
         Adapt to GammaLib-00-08-00 interface
+        ctobssim now only simulates events within ROI
 2012-12-21  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-07-00
 Author(s) : Juergen Knoedlseder
-February 21th, 2013
+April 9th, 2013
 1. Introduction
@@ -27,7 +27,11 @@ None
 4. ctobssim - CTA observation simulation
 ----------------------------------------
-None
+The application now checks for each event whether it lies within the
+circular region of interested defined by the region centre (parameters
+ra and dec) and region radius (parameter rad). A margin of 0.5 has been
+added internally to the radius so that also sources at the edge but
+slightly outside of the ROI will be simulated correctly.
 5. ctselect - CTA event selection
@@ -36,6 +36,9 @@
 #define G_SIMULATE_SOURCE       "ctobssim::simulate_source(GCTAObservation*,"\
                                                                 " GPhotons&)"
+/* __ Constants __________________________________________________________ */
+const double g_roi_margin = 0.5;      //!< Simulation radius margin (degrees)
+
 /* __ Debug definitions __________________________________________________ */
 /* __ Coding definitions _________________________________________________ */
@@ -66,10 +69,12 @@ ctobssim::ctobssim(void) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION)
 /***********************************************************************//**
  * @brief Observations constructor
  *
- * This constructor creates an instance of the class that is initialised from
- * an observation container.
+ * @param[in] obs Observation container.
+ *
+ * Constructs application by initialising from an observation container.
  ***************************************************************************/
-ctobssim::ctobssim(GObservations obs) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION)
+ctobssim::ctobssim(GObservations obs) :
+          GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION)
 {
     // Initialise members
     init_members();
@@ -90,9 +95,11 @@ ctobssim::ctobssim(GObservations obs) : GApplication(CTOBSSIM_NAME, CTOBSSIM_VER
  *
  * @param[in] argc Number of arguments in command line.
  * @param[in] argv Array of command line arguments.
+ *
+ * Constructs application using command line arguments for parameter setting.
  ***************************************************************************/
 ctobssim::ctobssim(int argc, char *argv[]) : 
-                    GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
+          GApplication(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -146,6 +153,7 @@ ctobssim::~ctobssim(void)
  * @brief Assignment operator
  *
  * @param[in] app Application.
+ * @return Application.
  ***************************************************************************/
 ctobssim& ctobssim::operator= (const ctobssim& app)
 {
@@ -178,7 +186,7 @@ ctobssim&amp; ctobssim::operator= (const ctobssim&amp; app)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear instance
+ * @brief Clear application
  ***************************************************************************/
 void ctobssim::clear(void)
 {
@@ -598,11 +606,6 @@ void ctobssim::set_list(GCTAObservation* obs)
  *
  * This method does nothing if the observation pointer is NULL. It also
  * verifies if the observation has a valid pointing and response.
- *
- * @todo Add margins to the data selections so that we generate also photons
- *       that fall slightly outside the ROI and the energy interval so that
- *       they can be scatter within the ROI and energy interval by the
- *       PSF and energy dispersion.
  ***************************************************************************/
 void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models,
                                GRan& ran, GLog* wrklog)
@@ -631,9 +634,9 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
         // Get pointer on event list (circumvent const correctness)
         GCTAEventList* events = static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
-        // Extract ROI
+        // Extract simulation region.
         GSkyDir dir = events->roi().centre().dir();
-        double  rad = events->roi().radius();
+        double  rad = events->roi().radius() + g_roi_margin;
         // Dump simulation cone information
         if (logNormal()) {
@@ -649,7 +652,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
         int indent = 0;
         // Loop over all Good Time Intervals
-        for (int it = 0; it <  events->gti().size(); ++it) {
+        for (int it = 0; it < events->gti().size(); ++it) {
             // Extract time interval
             GTime tmin = events->gti().tstart(it);
@@ -744,10 +747,14 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
                                                                *obs,
                                                                ran);
                                 if (event != NULL) {
-                                    event->event_id(m_event_id);
-                                    events->append(*event);
+                                    // Use event only if it falls within ROI
+                                    if (event->dir().dist_deg(events->roi().centre().dir()) <=
+                                        events->roi().radius()) {
+                                        event->event_id(m_event_id);
+                                        events->append(*event);
+                                        m_event_id++;
+                                    }
                                     delete event;
-                                    m_event_id++;
                                 }
                             } // endfor: looped over events
@@ -886,12 +893,18 @@ void ctobssim::simulate_background(GCTAObservation* obs,
                     // Get event pointer
                     GCTAEventAtom* event = (*list)[k];
-                    // Set event identifier
-                    event->event_id(m_event_id);
-                    m_event_id++;
+                    // Use event only if it falls within ROI
+                    if (event->dir().dist_deg(events->roi().centre().dir()) <=
+                        events->roi().radius()) {
+
+                        // Set event identifier
+                        event->event_id(m_event_id);
+                        m_event_id++;
+
+                        // Append event
+                        events->append(*event);
-                    // Append event
-                    events->append(*event);
+                    } // endif: event was within ROI
                 } // endfor: looped over all events
@@ -57,7 +57,7 @@ public:
     virtual ~ctobssim(void);
     // Operators
-    ctobssim& operator= (const ctobssim& app);
+    ctobssim& operator=(const ctobssim& app);
     // Methods
     void           clear(void);
@@ -87,7 +87,7 @@ protected:
     int           m_seed;       //!< Random number generator seed 
     double        m_ra;         //!< RA of pointing direction
     double        m_dec;        //!< DEC of pointing direction
-    double        m_rad;        //!< FOV radius
+    double        m_rad;        //!< FOV radius (degrees)
     double        m_tmin;       //!< Start time (MET)
     double        m_tmax;       //!< Stop time (MET)
     double        m_emin;       //!< Lower energy (TeV)
@@ -5,7 +5,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctobssim.par                                                  #
-#  Version:   00-07-00                                                      #
+#  Version:   00-08-00                                                      #
 #                                                                           #
 #  Author:    Jurgen Knodlseder                                             #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #