Optionally replace point sources by small Gaussians and add reference documentation.

Also some code speed up by considering only relevant pixels for radial and elliptical models.

Optionally replace point sources by small Gaussians and add reference documentation.
Also some code speed up by considering only relevant pixels for radial and elliptical models.
Jürgen Knödlseder
1 parent 71ceba5f
Showing 5 changed files with 294 additions and 66 deletions
doc/source/reference_manual/ctmapcube.rst
doc/source/reference_manual/index.rst
src/ctmapcube/ctmapcube.cpp
src/ctmapcube/ctmapcube.hpp
src/ctmapcube/ctmapcube.par
+.. _ctmapcube:
+
+ctmapcube
+=========
+
+Generate a map cube from an input model.
+
+
+Synopsis
+--------
+
+This tool generates a map cube from all sky model components contained in
+an input model. A map cube is a series of sky maps, each spanned by Right
+Ascension or Galactic longitude and Declination or Galactic latitude. Each
+map corresponds to a specific photon energy. The energy binning of the map
+cube may be either linear, logarithmic, or custom defined using an input
+file.
+
+ctmapcube generates a map cube FITS file comprising two extensions. The
+primary extension contains a 3-dimensional image containing the map cube
+values. A second extension named ``ENERGIES`` contains a binary table that
+defines the energy values of each map in the cube.
+
+
+General parameters
+------------------
+
+``inmodel [file]``
+    Input model XML file.
+
+``outcube [file]``
+    Output map cube file.
+
+``(ptsrcsig = 1.0) [real]``
+    Sigma of Gaussian to be used for point sources (in arcmin).
+
+``ebinalg <FILE|LIN|LOG> [string]``
+    Algorithm for defining energy bins.
+ 	 	 
+``emin [real]``
+    Lower energy value for first energy bin (in TeV).
+ 	 	 
+``emax [real]``
+    Upper energy value for last energy bin (in TeV).
+ 	 	 
+``enumbins [integer]``
+    Number of energy bins.
+ 	 	 
+``ebinfile [file]``
+    Name of the file containing the energy bin definition.
+ 	 	 
+``nxpix [integer]``
+    Number of cube bins in Right Ascension or Galactic longitude.
+ 	 	 
+``nypix [integer]``
+    Number of cube bins in Declination or Galactic latitude.
+ 	 	 
+``binsz [real]``
+    Image scale (in degrees/pixel).
+ 	 	  	 	 
+``coordsys <CEL|GAL> [string]``
+    Coordinate system (CEL - celestial, GAL - galactic).
+ 	 	 
+``xref [real]``
+    Right Ascension / Galactic longitude of cube centre (J2000, in degrees).
+ 	 	 
+``yref [real]``
+    Declination / Galactic latitude of cube centre (J2000, in degrees).
+ 	 	 
+``proj <AIT|AZP|CAR|MER|MOL|STG|TAN> [string]``
+    Projection method.
+ 	 	 
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+    Specifies whether the map cube should be published on VO Hub.
+
+``(chatter = 2) [integer]``
+    Verbosity of the executable:
+     chatter = 0: no information will be logged
+     
+     chatter = 1: only errors will be logged
+     
+     chatter = 2: errors and actions will be logged
+     
+     chatter = 3: report about the task execution
+     
+     chatter = 4: detailed report about the task execution
+ 	 	 
+``(clobber = yes) [boolean]``
+    Specifies whether an existing output counts cube file should be overwritten.
+ 	 	 
+``(debug = no) [boolean]``
+    Enables debug mode. In debug mode the executable will dump any log file output to the console.
+ 	 	 
+``(mode = ql) [string]``
+    Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = ctmapcube.log) [string]``
+    Name of log file.
+
+
+Related tools or scripts
+------------------------
+
+None
@@ -20,16 +20,17 @@ ctools
    ctbutterfly --- Compute butterfly <ctbutterfly>
    ctcubemask --- Filter counts cube <ctcubemask>
    ctedispcube --- Generates energy dispersion cube <ctedispcube>
+   cterror --- Calculates likelihood profile errors <cterror>
    ctexpcube --- Generates exposure cube <ctexpcube>
    ctlike --- Performs maximum likelihood fitting <ctlike>
+   ctmapcube --- Generates a map cube <ctmapcube>
    ctmodel --- Computes model counts cube <ctmodel>
-   ctobssim --- Simulate CTA observations <ctobssim>
+   ctobssim --- Simulate observations <ctobssim>
    ctpsfcube --- Generates point spread function cube <ctpsfcube>
    ctselect --- Selects event data <ctselect>
    ctskymap --- Generates sky map <ctskymap>
-   cttsmap --- Generates Test Statistics map <cttsmap>
+   cttsmap --- Generates Test Statistic map <cttsmap>
    ctulimit --- Calculates upper limit <ctulimit>
-   cterror --- Calculates likelihood profile errors <cterror>
 General cscripts
@@ -151,7 +151,7 @@ ctmapcube&amp; ctmapcube::operator=(const ctmapcube&amp; app)
  ==========================================================================*/
 /***********************************************************************//**
- * @brief Clear instance
+ * @brief Clear tool
  ***************************************************************************/
 void ctmapcube::clear(void)
 {
@@ -245,6 +245,8 @@ void ctmapcube::save(void)
  * @brief Publish map cube
  *
  * @param[in] name Map cube name.
+ *
+ * Publishes the map cube on a VO Hub.
  ***************************************************************************/
 void ctmapcube::publish(const std::string& name)
 {
@@ -287,7 +289,8 @@ void ctmapcube::init_members(void)
 {
     // Initialise members
     m_outcube.clear();
-    m_publish = false;
+    m_ptsrcsig = 1.0;
+    m_publish  = false;
     // Initialise protected members
     m_models.clear();
@@ -306,8 +309,9 @@ void ctmapcube::init_members(void)
 void ctmapcube::copy_members(const ctmapcube& app)
 {
     // Copy attributes
-    m_outcube = app.m_outcube;
-    m_publish = app.m_publish;
+    m_outcube  = app.m_outcube;
+    m_ptsrcsig = app.m_ptsrcsig;
+    m_publish  = app.m_publish;
     // Copy protected members
     m_models = app.m_models;
@@ -332,10 +336,21 @@ void ctmapcube::free_members(void)
  * @brief Get application parameters
  *
  * Get all task parameters from parameter file or (if required) by querying
- * the user. The parameters are read in the correct order.
+ * the user.
  ***************************************************************************/
 void ctmapcube::get_parameters(void)
 {
+    // Read model definition file if required
+    if (m_models.size() == 0) {
+
+        // Get model filename
+        GFilename inmodel = (*this)["inmodel"].filename();
+
+        // Load models from file
+        m_models.load(inmodel);
+
+    } // endif: there were no models
+
     // Get energy binning parameters
     GEnergies energies = create_energies();
@@ -348,19 +363,9 @@ void ctmapcube::get_parameters(void)
     int         nxpix    = (*this)["nxpix"].integer();
     int         nypix    = (*this)["nypix"].integer();
-    // Read model definition file if required
-    if (m_models.size() == 0) {
-
-        // Get model filename
-        GFilename inmodel = (*this)["inmodel"].filename();
-
-        // Load models from file
-        m_models.load(inmodel);
-
-    } // endif: there were no models
-
-    // Get remaining parameters
-    m_publish = (*this)["publish"].boolean();
+    // Get remaining hidden parameters
+    m_ptsrcsig = (*this)["ptsrcsig"].real();
+    m_publish  = (*this)["publish"].boolean();
     // Read optionally output cube filenames
     if (read_ahead()) {
@@ -386,18 +391,19 @@ void ctmapcube::get_parameters(void)
 /***********************************************************************//**
- * @brief Create energy vector from user parameters
+ * @brief Create map cube energies from user parameters
  *
- * Get the energy vector according to the user parameters. The method
- * supports loading of energy information from the ENERGIES extension (or
- * any other extension specified by the "ebinfile" parameters), or setting
- * energies using a linear or logarithmical spacing.
+ * Creates the map cube energies according to the user parameters. If the
+ * @p ebinfile parameter is "FILE" the method loads the energy information
+ * from a FITS file. The name of the FITS file is given by the @p ebinfile
+ * parameter. If no extension name is specified the method will use assume
+ * the name to be "ENERGIES".
  *
- * The following parameters are read:
- *      ebinalg - Energy binning algorithm
- *      emin - Minimum energy (if ebinalg != FILE)
- *      emax - Maximum energy (if ebinalg != FILE)
- *      enumbins - Number of energies (if ebinalg != FILE)
+ * If the @p ebinfile parameter is "LIN" the method will create a linearly
+ * spaced series of @p enumbins energies in the interval @p emin to @p emax.
+ * Otherwise, the @p ebinfile parameter is assumed to be "LOG" and the method
+ * will create a logarithmically spaced series of @p enumbins energies in the
+ * interval @p emin to @p emax.
  ***************************************************************************/
 GEnergies ctmapcube::create_energies(void)
 {
@@ -510,6 +516,8 @@ void ctmapcube::create_cube(void)
 /***********************************************************************//**
  * @brief Add one model to map cube
  *
+ * @param[in] model Pointer to point source sky model.
+ *
  * Adds one model to map cube.
  ***************************************************************************/
 void ctmapcube::add_model(GModelSky* model)
@@ -520,21 +528,37 @@ void ctmapcube::add_model(GModelSky* model)
     // Get pointer to sky map
     GSkyMap* map = const_cast<GSkyMap*>(&m_cube.cube());
-    // Loop over all energy bins
-    for (int iebin = 0; iebin < energies.size(); ++iebin) {
+    // Get the model's bounding circle
+    GSkyDir centre;
+    double  radius;
+    get_bounding_circle(model, &centre, &radius);
-        // Loop over all spatial pixels
-        for (int ipixel = 0; ipixel < map->npix(); ++ipixel) {
+    // Add some margin to circle radius
+    radius += 0.1;
+
+    // Loop over all map cube pixels
+    for (int ipixel = 0; ipixel < map->npix(); ++ipixel) {
+
+        // Get pixel sky direction
+        GSkyDir dir = map->inx2dir(ipixel);
+
+        // Skip pixel if it's outside the bounding circle
+        if (centre.dist_deg(dir) > radius) {
+            continue;
+        }
+
+        // Loop over all energy bins
+        for (int iebin = 0; iebin < energies.size(); ++iebin) {
             // Set photon
-            GPhoton photon(map->inx2dir(ipixel), energies[iebin], GTime());
+            GPhoton photon(dir, energies[iebin], GTime());
             // Add model value
             map->operator()(ipixel, iebin) += model->value(photon);
-        } // endfor: looped over all pixels
+        } // endfor: looped over all energies
-    } // endfor: looped over all energies
+    } // endfor: looped over all pixels
     // Return
     return;
@@ -544,51 +568,141 @@ void ctmapcube::add_model(GModelSky* model)
 /***********************************************************************//**
  * @brief Add one point source model to map cube
  *
+ * @param[in] model Pointer to point source sky model.
+ *
  * Adds one point source model to map cube.
  ***************************************************************************/
 void ctmapcube::add_ptsrc_model(GModelSky* model)
 {
-    // Get point source spatial component
-    GModelSpatialPointSource* ptsrc = static_cast<GModelSpatialPointSource*>(model->spatial());
+    // Get pointer to sky map
+    GSkyMap* map = const_cast<GSkyMap*>(&m_cube.cube());
-    // Get point source direction
-    GSkyDir dir = ptsrc->dir();
+    // If Gaussian sigma is positive then replace the point sources by a
+    // Gaussian with the user defined width
+    if (m_ptsrcsig > 0.0) {
-    // Get cube energies
-    GEnergies energies = m_cube.energies();
+        // Get the point source centre (radius is ignored)
+        GSkyDir centre;
+        double  radius;
+        get_bounding_circle(model, &centre, &radius);
-    // Get pointer to sky map
-    GSkyMap* map = const_cast<GSkyMap*>(&m_cube.cube());
+        // Set the Gaussian sigma in degrees (m_ptsrcsig is in arcmin)
+        double sigma = m_ptsrcsig/60.0;
-    // Get test sky map for flux computation
-    GSkyMap test_map(*map);
+        // Create a Gaussian spatial model
+        GModelSpatialRadialGauss spatial(centre, sigma);
-    // Continue only if map contains point source sky direction
-    if (map->contains(dir)) {
+        // Create a new model reusing the spectral and temporal components
+        // of the original model
+        GModelSky skymodel(spatial, *(model->spectral()), *(model->temporal()));
-        // Get map pixel
-        int ipixel = map->dir2inx(dir);
+        // Add model
+        add_model(&skymodel);
-        // Loop over all energy bins
-        for (int iebin = 0; iebin < energies.size(); ++iebin) {
+    } // endif: Gaussian sigma was positive
-            // Set photon using the point source sky direction
-            GPhoton photon(dir, energies[iebin], GTime());
+    // ... otherwise set a single map pixel
+    else {
-            // Compute model value
-            double value = model->value(photon);
+        // Get point source spatial component
+        GModelSpatialPointSource* ptsrc =
+                    static_cast<GModelSpatialPointSource*>(model->spatial());
-            // Compute flux normalization
-            test_map(ipixel, iebin) = value;
-            double flux = test_map.flux(ipixel, iebin);
-            double norm = (flux > 0.0) ? value / flux : 0.0;
+        // Get point source direction
+        GSkyDir dir = ptsrc->dir();
-            // Add model to sky map
-            map->operator()(ipixel, iebin) += norm * value;
+        // Get cube energies
+        GEnergies energies = m_cube.energies();
-        } // endfor: looped over all energies
+        // Get test sky map for flux computation
+        GSkyMap test_map(*map);
+
+        // Continue only if map contains point source sky direction
+        if (map->contains(dir)) {
+
+            // Get map pixel
+            int ipixel = map->dir2inx(dir);
-    } // endif: map contained sky direction
+            // Loop over all energy bins
+            for (int iebin = 0; iebin < energies.size(); ++iebin) {
+
+                // Set photon using the point source sky direction
+                GPhoton photon(dir, energies[iebin], GTime());
+
+                // Compute model value
+                double value = model->value(photon);
+
+                // Compute flux normalization
+                test_map(ipixel, iebin) = value;
+                double flux = test_map.flux(ipixel, iebin);
+                double norm = (flux > 0.0) ? value / flux : 0.0;
+
+                // Add model to sky map
+                map->operator()(ipixel, iebin) += norm * value;
+
+            } // endfor: looped over all energies
+
+        } // endif: map contained sky direction
+
+    } // endelse: Gaussian sigma was non positive
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get bounding circle for a sky model
+ *
+ * @param[in] model Pointer to sky model.
+ * @param[out] dir Centre of bounding circle.
+ * @param[out] radius Radius of bounding circle (degrees).
+ *
+ * Determine for any spatial model the bounding circle. The bounding circle
+ * is a circular region that fully enclosed the model.
+ *
+ * The method does not verify the validity of the pointer arguments.
+ ***************************************************************************/
+void ctmapcube::get_bounding_circle(GModelSky* model,
+                                    GSkyDir*   dir,
+                                    double*    radius) const
+{
+    // Interpret the spatial model component as point source, radial source
+    // or elliptical source
+    GModelSpatialPointSource* ptsrc =
+                 dynamic_cast<GModelSpatialPointSource*>(model->spatial());
+    GModelSpatialRadial* radial =
+                 dynamic_cast<GModelSpatialRadial*>(model->spatial());
+    GModelSpatialElliptical* elliptical =
+                 dynamic_cast<GModelSpatialElliptical*>(model->spatial());
+
+    // If the spatial component as a point source then extract the source
+    // position and set the radius of the bounding box to 0.1 degrees
+    if (ptsrc != NULL) {
+        *dir    = ptsrc->dir();
+        *radius = 0.1;
+    }
+
+    // ... otherwise, if the spatial component is a radial source then extract
+    // the source position and set the radius of the bounding box from the
+    // value returned by theta_max() in radians
+    else if (radial != NULL) {
+        *dir    = radial->dir();
+        *radius = radial->theta_max() * gammalib::rad2deg;
+    }
+
+    // ... otherwise, if the spatial component is an elliptical source then
+    // extract the source position and set the radius of the bounding box from
+    // the value returned by theta_max() in radians
+    else if (elliptical != NULL) {
+        *dir    = elliptical->dir();
+        *radius = elliptical->theta_max() * gammalib::rad2deg;
+    }
+
+    // ... otherwise use the entire sky
+    else {
+        *radius = 180.0;
+    }
     // Return
     return;
@@ -71,9 +71,13 @@ protected:
     void      create_cube(void);
     void      add_model(GModelSky* model);
     void      add_ptsrc_model(GModelSky* model);
+    void      get_bounding_circle(GModelSky* model,
+                                  GSkyDir*   dir,
+                                  double*    radius) const;
     // User parameters
     GFilename m_outcube;                //!< Output map cube filename
+    double    m_ptsrcsig;               //!< Point source sigma (arcmin)
     bool      m_publish;                //!< Publish map cube?
     // Protected members
@@ -19,6 +19,7 @@
 #===================
 inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
 outcube,  f, a, "mapcube.fits",,, "Output map cube file"
+ptsrcsig, r, h, 1.0,0,60, "Point source Gaussian sigma in arcmin"
 #
 # Energy binning parameters