Integrate ctpsfcube tool.

* Checked the code for coding conventions. * Renamed the angular seperation user parameters to amax and anumbins. "a" stands for angular separation, "d" for delta is less explicit. * Add ctpsfcube executable unit test

Integrate ctpsfcube tool.
* Checked the code for coding conventions. * Renamed the angular seperation user parameters to amax and anumbins. "a" stands for angular separation, "d" for delta is less explicit. * Add ctpsfcube executable unit test
Jürgen Knödlseder
1 parent 73b1c340
Showing 13 changed files with 116 additions and 54 deletions
.gitignore
ChangeLog
NEWS
pyext/ctcubemask.i
pyext/ctpsfcube.i
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
src/ctcubemask/ctcubemask.par
src/ctpsfcube/ctpsfcube.cpp
src/ctpsfcube/ctpsfcube.hpp
src/ctpsfcube/ctpsfcube.par
src/ctpsfcube/main.cpp
test/test_ctools.sh
@@ -53,6 +53,7 @@ src/ctobssim/ctobssim
 src/ctselect/ctselect
 src/ctskymap/ctskymap
 src/ctexpcube/ctexpcube
+src/ctpsfcube/ctpsfcube
 src/ctcubemask/ctcubemask
 test/pfiles/
 test/*.fits
@@ -5,8 +5,9 @@
 	Adapt to GammaLib 0.9.0 interface
         ctbin now always merges event lists into a single counts cube
-        Add ctexpcube tool
-        Add ctcubemask tool
+        Add ctexpcube tool (generates exposure cube)
+        Add ctpsfcube tool (generated PSF cube)
+        Add ctcubemask tool (masks event bins in event cube)
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
         Enhance documentation
@@ -80,27 +80,37 @@ the use for 2D response functions in FITS files.
 ctmodel now also supports the handling of energy dispersion. Energy
 dispersion can be switched on or off using the hidden edisp parameter.
-10. ctexpcube - CTA exposure cube
----------------------------------
+
+10. ctexpcube - CTA exposure cube generation
+--------------------------------------------
 The first version of the ctexpcube tool is included in the release.
 ctexpcube computes the exposure cube for an observation definition
-file. The cube definition may be taken from a counts map or may be
+file. The cube definition may be taken from a counts cube or may be
 specified by user parameters.
-11. ctcubemask - CTA cube filtering
+11. ctpsfcube - CTA point spread function cube generation
+---------------------------------------------------------
+The first version of the ctpsfcube tool is included in the release.
+ctpsfcube computes the point spread function cube for an observation
+definition file. The cube definition may be taken from a counts cube
+or may be specified by user parameters.
+
+
+12. ctcubemask - CTA cube filtering
 -----------------------------------
 The first version of the ctcubemask tool is included in the release.
 ctcubemask allows filtering of bins in a counts cube by setting the
 bin content to -1. Filtering can be done by energy, for a circular
 Region of Interest (ROI) or using ds9 exclusion regions.
-12. Examples
+
+13. Examples
 ------------
 None
-13. Scripts
+14. Scripts
 -----------
 Add logL and event statistics to cspull output file.
@@ -109,11 +119,11 @@ Convert all tabs into spaces.
 The obsutils script has now been moved into the ctools module namespace.
-14. caldb
+15. caldb
 ---------
 None
-15. models
+16. models
 ----------
 None
 /***************************************************************************
- *                   ctcubemask.i - CTA mask cube tool                     *
+ *                    ctcubemask - CTA mask cube tool                      *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 Chia-Chun Lu                                        *
  * ----------------------------------------------------------------------- *
@@ -39,7 +39,7 @@ class ctcubemask : public GApplication  {
 public:
     // Constructors and destructors
     ctcubemask(void);
-    explicit ctcubemask(GObservations obs);
+    explicit ctcubemask(const GObservations& obs);
     ctcubemask(int argc, char *argv[]);
     ctcubemask(const ctcubemask& app);
     virtual ~ctcubemask(void);
 /***************************************************************************
- *                   ctpsfcube - CTA PSF cube tool                         *
+ *                     ctpsfcube - CTA PSF cube tool                       *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -40,7 +40,7 @@ class ctpsfcube : public GApplication  {
 public:
     // Constructors and destructors
     ctpsfcube(void);
-    explicit ctpsfcube(GObservations obs);
+    explicit ctpsfcube(const GObservations& obs);
     ctpsfcube(int argc, char *argv[]);
     ctpsfcube(const ctpsfcube& app);
     virtual ~ctpsfcube(void);
 /***************************************************************************
- *                  ctcubemask.cpp - CTA cube filter tool                  *
+ *                    ctcubemask - CTA cube filter tool                    *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -70,7 +70,8 @@ ctcubemask::ctcubemask(void) : GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
  * This constructor creates an instance of the class that is initialised from
  * an observation container.
  ***************************************************************************/
-ctcubemask::ctcubemask(GObservations obs) : GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
+ctcubemask::ctcubemask(const GObservations& obs) :
+            GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
 {
     // Initialise members
     init_members();
@@ -302,7 +303,7 @@ void ctcubemask::run(void)
     // Write observation(s) into logger
     if (logTerse()) {
         log << std::endl;
-        log.header1("Observations after masked");
+        log.header1("Observations after event bin masking");
         log << m_obs << std::endl;
     }
 /***************************************************************************
- *                  ctcubemask.hpp - CTA cube filter tool                  *
+ *                    ctcubemask - CTA cube filter tool                    *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -47,7 +47,7 @@ class ctcubemask : public GApplication  {
 public:
     // Constructors and destructors
     ctcubemask(void);
-    explicit ctcubemask(GObservations obs);
+    explicit ctcubemask(const GObservations& obs);
     ctcubemask(int argc, char *argv[]);
     ctcubemask(const ctcubemask& app);
     virtual ~ctcubemask(void);
@@ -9,7 +9,7 @@
 #                                                                           #
 #  Author:    Chia-Chun Lu                                                  #
 #             c.lu@mpi-hd.mpg.de                                            #
-#             IRAP                                                          #
+#             MPIK                                                          #
 #                                                                           #
 #  Purpose:   Parameter file for the ctcubemask application                 #
 #                                                                           #
 /***************************************************************************
- *                   ctpsfcube - CTA PSF cube tool                         *
+ *                    ctpsfcube - CTA PSF cube tool                        *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -68,10 +68,12 @@ ctpsfcube::ctpsfcube(void) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
 /***********************************************************************//**
  * @brief Observations constructor
  *
+ * @param[in] obs Observation container.
+ *
  * This method creates an instance of the class by copying an existing
  * observations container.
  ***************************************************************************/
-ctpsfcube::ctpsfcube(GObservations obs) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
+ctpsfcube::ctpsfcube(const GObservations& obs) : GApplication(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
 {
     // Initialise members
     init_members();
@@ -293,7 +295,6 @@ void ctpsfcube::run(void)
 /***********************************************************************//**
  * @brief Save PSF cube
- *
  ***************************************************************************/
 void ctpsfcube::save(void)
 {
@@ -349,19 +350,20 @@ void ctpsfcube::get_parameters(void)
         double      binsz    = (*this)["binsz"].real();
         int         nxpix    = (*this)["nxpix"].integer();
         int         nypix    = (*this)["nypix"].integer();
-	double      dmax     = (*this)["dmax"].real();
-	int         ndbins   = (*this)["ndbins"].integer();
+        double      dmax     = (*this)["amax"].real();
+        int         ndbins   = (*this)["anumbins"].integer();
         // Get energy definition
         get_ebounds();
-        // Define exposure cube
+        // Define PSF cube
         m_psfcube = GCTAMeanPsf(wcs, coordsys, xref, yref,
-				-binsz, binsz, nxpix, nypix,
-				m_ebounds, dmax, ndbins);
-    }
+                                -binsz, binsz, nxpix, nypix,
+                                m_ebounds, dmax, ndbins);
-    // ... otherwise setup the exposure cube from the counts map
+    } // endif: PSF cube set from user parameters
+
+    // ... otherwise setup the PSF cube from the counts map
     else {
         // Set PSF cube from counts map
@@ -401,7 +403,6 @@ void ctpsfcube::get_obs(void)
         // Append CTA observation to container
         m_obs.append(obs);
-            
     }
     // ... otherwise try to open as XML file
@@ -420,7 +421,7 @@ void ctpsfcube::get_obs(void)
 /***********************************************************************//**
  * @brief Set observation response
  *
- * Set response for all observations that so have no response.
+ * Set response for all observations that have no response.
  ***************************************************************************/
 void ctpsfcube::set_response(void)
 {
@@ -537,12 +538,12 @@ void ctpsfcube::get_ebounds(void)
 /***********************************************************************//**
- * @brief Set exposure cube definition from counts map
+ * @brief Set PSF cube definition from counts map
  *
  * @exception GException::invalid_value
  *            Invalid counts map projection or invalid events encountered.
  *
- * Set exposure cube definition from counts map.
+ * Set PSF cube definition from counts map.
  ***************************************************************************/
 void ctpsfcube::set_from_cntmap(const std::string& filename)
 {
@@ -552,7 +553,7 @@ void ctpsfcube::set_from_cntmap(const std::string&amp; filename)
     // Load counts map in CTA observation
     obs.load(filename);
-    // Set exposure cube from counts map
+    // Set PSF cube from counts map
     const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
     // Continue only if cube is valid
@@ -573,16 +574,16 @@ void ctpsfcube::set_from_cntmap(const std::string&amp; filename)
             double      dy       = wcs->cdelt(1);
             int         nx       = cube->map().nx();
             int         ny       = cube->map().ny();
-	    double      dmax     = (*this)["dmax"].real();
-	    int         ndbins   = (*this)["ndbins"].integer();
+            double      amax     = (*this)["amax"].real();
+            int         anumbins = (*this)["anumbins"].integer();
             // Get energy definition
             m_ebounds = cube->ebounds();
-            // Define exposure cube
+            // Define PSF cube
             m_psfcube = GCTAMeanPsf(proj, coordsys, xref, yref,
-				    dx, dy, nx, ny,
-				    m_ebounds, dmax, ndbins);
+                                    dx, dy, nx, ny,
+                                    m_ebounds, amax, anumbins);
         } // endif: WCS projection was valid
 /***************************************************************************
- *                   ctpsfcube - CTA PSF cube tool                         *
+ *                     ctpsfcube - CTA PSF cube tool                       *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -46,7 +46,7 @@ class ctpsfcube : public GApplication  {
 public:
     // Constructors and destructors
     ctpsfcube(void);
-    explicit ctpsfcube(GObservations obs);
+    explicit ctpsfcube(const GObservations& obs);
     ctpsfcube(int argc, char *argv[]);
     ctpsfcube(const ctpsfcube& app);
     virtual ~ctpsfcube(void);
@@ -4,14 +4,14 @@
 #                                                                           #
 # ------------------------------------------------------------------------- #
 #                                                                           #
-#  File:      ctexpcube.par                                                 #
+#  File:      ctpsfcube.par                                                 #
 #  Version:   00-01-00                                                      #
 #                                                                           #
-#  Author:    Juergen Knoedlseder                                           #
-#             Jurgen.Knodlseder@irap.omp.eu                                 #
-#             IRAP                                                          #
+#  Author:    Chia-Chun Lu                                                  #
+#             c.lu@mpi-hd.mpg.de                                            #
+#             MPIK                                                          #
 #                                                                           #
-#  Purpose:   Parameter file for the ctexpcube application                  #
+#  Purpose:   Parameter file for the ctpsfcube application                  #
 #                                                                           #
 #############################################################################
 #
@@ -19,7 +19,7 @@
 #===================
 infile,   f, a, "NONE",,, "Observation definition file"
 cntmap,   f, a, "NONE",,, "Counts map for exposure cube definition"
-outfile,  f, a, "expcube.fits",,, "Output exposure cube file"
+outfile,  f, a, "psfcube.fits",,, "Output exposure cube file"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 edisp,    b, h, no,,, "Apply energy dispersion?"
@@ -43,10 +43,10 @@ yref,     r, a, 22.01,,, &quot;Second coordinate of image center in degrees (DEC or g
 axisrot,  r, h, 0.0,,, "Rotation angle of image axis, in degrees"
 proj,     s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
 #
-# Histogram binning parameters
-#=============================
-dmax,     r, a, 0.3,,, "Upper bound of delta"
-ndbins,   i, a, 60,,,  "Number of delta bins"
+# Angular separation binning parameters
+#======================================
+amax,     r, h, 0.3,,, "Upper bound of angular separation between true and measued photon direction (degrees)"
+anumbins, i, h, 60,,,  "Number of angular separation bins"
 #
 # Standard parameters
 #====================
@@ -54,4 +54,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
 debug,   b, h, no,,, "Debugging mode activated"
 mode,    s, h, "ql",,, "Mode of automatic parameters"
-logfile, f, h, "ctexpcube.log",,, "Log filename"
+logfile, f, h, "ctpsfcube.log",,, "Log filename"
 /***************************************************************************
- *                  ctpsfcube - CTA PSF cube tool                          *
+ *                    ctpsfcube - CTA PSF cube tool                        *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Chia-Chun Lu                                     *
  * ----------------------------------------------------------------------- *
@@ -41,6 +41,7 @@ ctobssim=../src/ctobssim/ctobssim
 ctselect=../src/ctselect/ctselect
 ctskymap=../src/ctskymap/ctskymap
 ctexpcube=../src/ctexpcube/ctexpcube
+ctpsfcube=../src/ctpsfcube/ctpsfcube
 ctcubemask=../src/ctcubemask/ctcubemask
@@ -392,6 +393,53 @@ $ECHO &quot; ok&quot;
 #
+# Test ctpsfcube
+# ==============
+$ECHO -n "Test ctpsfcube: "
+#
+# Run 1
+$ctpsfcube infile="data/crab_events.fits" \
+           cntmap="NONE" \
+           outfile="psfcube1.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=200 \
+           nypix=200 \
+           binsz=0.02 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR"
+$ECHO -n "."
+if [ -s "psfcube1.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " psfcube1.fits file is not valid"
+  exit 1
+fi
+#
+# Run 2
+$ctpsfcube infile="data/crab_events.fits" \
+           cntmap="data/crab_cntmap.fits" \
+           outfile="psfcube2.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "psfcube2.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " psfcube2.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
 # Test ctcubemask
 # ===============
 $ECHO -n "Test ctcubemask: "