From e991511e54f5ee19e93596921bb0166b32c91af5 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 20 Jul 2014 03:18:12 +0200 Subject: [PATCH] Adapt to ctools coding conventions, clean-up and test code, and add a unit test. --- AUTHORS | 1 + ChangeLog | 3 ++- NEWS | 17 ++++++++++++----- configure.ac | 1 + pyext/ctcubemask.i | 8 ++++---- src/Makefile.am | 4 ++-- src/ctcubemask/ctcubemask.cpp | 241 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------- src/ctcubemask/ctcubemask.hpp | 8 ++++---- src/ctcubemask/ctcubemask.par | 14 +++++++------- test/data/exclusion.reg | 4 ++++ test/test_ctools.sh | 23 +++++++++++++++++++++++ 11 files changed, 170 insertions(+), 154 deletions(-) create mode 100644 test/data/exclusion.reg diff --git a/AUTHORS b/AUTHORS index 72430a1..7a3d3f6 100644 --- a/AUTHORS +++ b/AUTHORS @@ -2,3 +2,4 @@ Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) Michael Mayer (michael.mayer@desy.de) +Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de) diff --git a/ChangeLog b/ChangeLog index 73dd803..eea2455 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-07-18 Juergen Knoedlseder +2014-07-20 Juergen Knoedlseder * Version 0.8.0 released ======================== @@ -6,6 +6,7 @@ Adapt to GammaLib 0.9.0 interface ctbin now always merges event lists into a single counts cube Add ctexpcube tool + Add ctcubemask tool Support energy dispersion in ctobssim, ctlike and ctmodel Correct ctobssim for thread safety Enhance documentation diff --git a/NEWS b/NEWS index 4ccdd38..48573eb 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -18 July 2014 +20 July 2014 1. Introduction @@ -88,12 +88,19 @@ file. The cube definition may be taken from a counts map or may be specified by user parameters. -11. Examples +11. 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 ------------ None -12. Scripts +13. Scripts ----------- Add logL and event statistics to cspull output file. @@ -102,11 +109,11 @@ Convert all tabs into spaces. The obsutils script has now been moved into the ctools module namespace. -13. caldb +14. caldb --------- None -14. models +15. models ---------- None diff --git a/configure.ac b/configure.ac index 2f33138..9023bb7 100644 --- a/configure.ac +++ b/configure.ac @@ -883,6 +883,7 @@ AC_CONFIG_FILES([Makefile src/ctskymap/Makefile src/ctselect/Makefile src/ctexpcube/Makefile + src/ctcubemask/Makefile src/ctools-setup test/Makefile test/test_cscripts.sh diff --git a/pyext/ctcubemask.i b/pyext/ctcubemask.i index 5b4139a..a7582ed 100644 --- a/pyext/ctcubemask.i +++ b/pyext/ctcubemask.i @@ -1,5 +1,5 @@ /*************************************************************************** - * ctcubemask - CTA mask cube tool * + * ctcubemask.i - CTA mask cube tool * * ----------------------------------------------------------------------- * * copyright (C) 2014 Chia-Chun Lu * * ----------------------------------------------------------------------- * @@ -20,7 +20,7 @@ ***************************************************************************/ /** * @file ctcubemask.i - * @brief CTA data selection tool Python interface definition + * @brief CTA mask cube interface definition * @author Chia-Chun Lu */ %{ @@ -49,9 +49,9 @@ public: void execute(void); void run(void); void save(void); - GObservations& obs(void) { return m_obs; } + GObservations& obs(void); void get_parameters(void); - void apply_mask(GCTAObservation* obs, const std::string& filename); + void apply_mask(GCTAObservation* obs); }; diff --git a/src/Makefile.am b/src/Makefile.am index 73d13d6..7d789f6 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -12,7 +12,7 @@ SUBDIRS = ctbin \ ctselect \ ctskymap \ ctexpcube \ - ctcubemask + ctcubemask # Scripts to be distributed scripts = ctools-setup \ @@ -31,7 +31,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \ ctselect/libctselect.la \ ctskymap/libctskymap.la \ ctexpcube/libctexpcube.la \ - ctcubemask/libctcubemask.la + ctcubemask/libctcubemask.la # Distribute scripts EXTRA_DIST = $(scripts) diff --git a/src/ctcubemask/ctcubemask.cpp b/src/ctcubemask/ctcubemask.cpp index 2a88ce4..75587ad 100644 --- a/src/ctcubemask/ctcubemask.cpp +++ b/src/ctcubemask/ctcubemask.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctfilter - CTA cube filter tool * + * ctcubemask.cpp - CTA cube filter tool * * ----------------------------------------------------------------------- * * copyright (C) 2014 by Chia-Chun Lu * * ----------------------------------------------------------------------- * @@ -34,9 +34,8 @@ /* __ Method name definitions ____________________________________________ */ -#define G_RUN "ctcubemask::run()" -#define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*,"\ - "std::string&)" +#define G_RUN "ctcubemask::run()" +#define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*)" /* __ Debug definitions __________________________________________________ */ @@ -96,7 +95,7 @@ ctcubemask::ctcubemask(GObservations obs) : GApplication(CTCUBEMASK_NAME, CTCUBE * @param[in] argv Array of command line arguments. ***************************************************************************/ ctcubemask::ctcubemask(int argc, char *argv[]) : - GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv) + GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv) { // Initialise members init_members(); @@ -151,7 +150,7 @@ ctcubemask::~ctcubemask(void) * * @param[in] app Application. ***************************************************************************/ -ctcubemask& ctcubemask::operator= (const ctcubemask& app) +ctcubemask& ctcubemask::operator=(const ctcubemask& app) { // Execute only if object is not identical if (this != &app) { @@ -287,33 +286,9 @@ void ctcubemask::run(void) // Save event file name (for possible saving) m_infiles[i] = obs->eventfile(); - // Get temporary file name - std::string filename = std::tmpnam(NULL); - - // Save observation in temporary file - obs->save(filename, true); - - // Log saved FITS file - if (logExplicit()) { - GFits tmpfile(filename); - log << std::endl; - log.header1("FITS file content of temporary file"); - log << tmpfile << std::endl; - tmpfile.close(); - } - - // Check temporary file - std::string message = check_infile(filename); - if (message.length() > 0) { - throw GException::app_error(G_RUN, message); - } - // Load observation from temporary file, including event selection - apply_mask(obs, filename); + apply_mask(obs); - // Remove temporary file - std::remove(filename.c_str()); - } // endif: had a CTA observation } // endfor: looped over all observations @@ -409,7 +384,7 @@ void ctcubemask::get_parameters(void) try { // Load event list in CTA observation - obs.load_binned(m_infile); + obs.load(m_infile); // Append CTA observation to container m_obs.append(obs); @@ -455,115 +430,119 @@ void ctcubemask::get_parameters(void) /***********************************************************************//** - * @brief Select events + * @brief Apply mask to counts cube * * @param[in] obs CTA observation. - * @param[in] filename File name. * - * Select events from a FITS file by making use of the selection possibility - * of the cfitsio library on loading a file. A selection string is created - * from the specified criteria that is appended to the filename so that - * cfitsio will automatically filter the event data. This selection string - * is then applied when opening the FITS file. The opened FITS file is then - * saved into a temporary file which is the loaded into the actual CTA - * observation, overwriting the old CTA observation. The ROI, GTI and EBounds - * of the CTA event list are then set accordingly to the specified selection. - * Finally, the temporary file created during this process is removed. - * - * Good Time Intervals of the observation will be limited to the time - * interval [m_tmin, m_tmax]. If m_tmin=m_tmax=0, no time selection is - * performed. - * - * @todo Use INDEF instead of 0.0 for pointing as RA/DEC selection + * Masks out event bins by setting their content to -1. ***************************************************************************/ -void ctcubemask::apply_mask(GCTAObservation* obs, const std::string& filename) +void ctcubemask::apply_mask(GCTAObservation* obs) { - std::cout<< "start appply filter" << std::endl; - // Open FITS file - GFits file(filename); - GCTAEventCube* cube = (GCTAEventCube*) obs->events(); - GEbounds ebounds = cube->ebounds(); - GSkymap map = cube->map(); - - // Apply energy selection - int e_idx1 = 0; - int e_idx2 = 0; - int n_ebin = ebounds.size(); - - // Loop over ebounds to find the first energy band - for ( int i = 0 ; i < ebounds.size() ; i++){ - double emin = ebounds.emin(i).TeV(); - if ( emin > m_emin ){ - e_idx1 = i; - break; - } - } + // Get pointer to CTA event cube. Continue only if we have a cube + GCTAEventCube* cube = dynamic_cast(const_cast(obs->events())); + if (cube != NULL) { + + // Extract event cube and energy boundaries + GSkymap map = cube->map(); + const GEbounds& ebounds = cube->ebounds(); + + // Initialise energy selection + int n_ebin = ebounds.size(); + int e_idx1 = 0; + int e_idx2 = n_ebin; + + // Loop over ebounds to find the first valid energy band + for (int i = 0; i < n_ebin; ++i) { + double emin = ebounds.emin(i).TeV() + 1.0e-6; // Rounding tolerance + if (emin >= m_emin) { + e_idx1 = i; + break; + } + } - std::cout << "finish looping over ebounds to find 1st" << std::endl; - std::cout << "e_ix1 " << e_idx1 << std::endl; + // Loop over ebounds to find the last valid energy band + for (int i = n_ebin-1; i >= 0; i--) { + double emax = ebounds.emax(i).TeV() - 1.0e-6; // Rounding tolerance + if (emax <= m_emax) { + e_idx2 = i; + break; + } + } - // Loop over ebounds to find the last energy band - for ( int i = ebounds.size()-1 ; i < ebounds.size() ; i--){ - double emin = ebounds.emin(i).TeV(); - if ( emin < m_emax ){ - e_idx2 = i-1; - break; - } - } - std::cout << "finish looping over ebounds to find 1st" << std::endl; - std::cout << "e_ix2 " << e_idx2 << std::endl; - - // Set all pixels outside the desired energy bands negative - int npix = map.npix(); - for ( int i = 0 ; i < e_idx1 ; i++){ - for ( int pixel = 0 ; pixel < npix ; pixel++){ - map(pixel,i) = -1; - } - } - for ( int i = e_idx2 ; i < ebounds.size() ; i++){ - for ( int pixel = 0 ; pixel < npix ; pixel++){ - map(pixel,i) = -1; - } - } - std::cout << "finish energy band filter" << std::endl; - - - - // Set all pixels inside the desired energy bands - // but outside ROI or inside exlusion regions negative - GSkyRegions regs = GSkyRegions(m_regfile); - GSkyRegionCircle roi = GSkyRegionCircle(m_ra, m_dec, m_rad); - for ( int i = e_idx1 ; i <= e_idx2 ; i++){ - for ( int pixel = 0 ; pixel < npix ; pixel++){ - GSkyDir dir = map.inx2dir(pixel); - if ( roi.contains(dir) == false || - regs.contains(dir) == true){ - map(pixel,i) = -1; - } - } - } - std::cout << "finish roi selection" << std::endl; - - GCTAEventCube newcube( map, ebounds, obs->events()->gti()); - obs->events(newcube); + // Set all pixels outside the desired energy bands to -1.0 + int npix = map.npix(); + for (int i = 0; i < e_idx1; ++i) { + for (int pixel = 0; pixel < npix; ++pixel) { + map(pixel,i) = -1.0; + } + } + for (int i = e_idx2 + 1; i < n_ebin; ++i) { + for (int pixel = 0; pixel < npix; ++pixel) { + map(pixel,i) = -1.0; + } + } - // Log selected FITS file - if (logExplicit()) { - log << std::endl; - log.header1("FITS file content after filter"); - log << file << std::endl; - } + // Log selected energy band + if (logTerse()) { + log << gammalib::parformat("Selected energy band"); + log << ebounds.emin(e_idx1).TeV() << " - "; + log << ebounds.emax(e_idx2).TeV() << " TeV"; + log << std::endl; + } - // Check if we have an IMAGE HDU - if (!file.contains("IMAGE")) { - std::string message = "No \"IMAGE\" extension found in FITS file. "; - throw GException::app_error(G_APPLY_MASK, message); - } + // Set all pixels inside selected energy bands but outside ROI + // to -1.0 + GSkyRegionCircle roi(m_ra, m_dec, m_rad); + for (int i = e_idx1; i <= e_idx2; ++i) { + for (int pixel = 0; pixel < npix; ++pixel) { + GSkyDir dir = map.inx2dir(pixel); + if (!roi.contains(dir)) { + map(pixel,i) = -1.0; + } + } + } + + // Log selected energy band + if (logTerse()) { + log << gammalib::parformat("Selected ROI"); + log << "RA=" << m_ra << " deg, "; + log << "DEC=" << m_dec << " deg, "; + log << "Radius=" << m_rad << " deg"; + log << std::endl; + } + + // Set all pixels inside selected energy bands and inside exclusion + // regions to -1.0 + if (m_regfile != "NONE") { + GSkyRegions regions = GSkyRegions(m_regfile); + for (int i = e_idx1; i <= e_idx2; ++i) { + for (int pixel = 0; pixel < npix; ++pixel) { + GSkyDir dir = map.inx2dir(pixel); + if (regions.contains(dir)) { + map(pixel,i) = -1.0; + } + } + } + if (logTerse()) { + log << gammalib::parformat("Exclusion regions"); + log << std::endl; + log.indent(1); + log << regions.print(); + log << std::endl; + log.indent(0); + } + } + else { + if (logTerse()) { + log << gammalib::parformat("Exclusion regions"); + log << "None" << std::endl; + } + } + // Put back map into the event cube + cube->map(map); - // Get CTA event list pointer - GCTAEventList* list = - static_cast(const_cast(obs->events())); + } // endif: observation contained an event cube // Return return; @@ -798,7 +777,7 @@ void ctcubemask::save_xml(void) * nothing if the observation pointer is not valid. ***************************************************************************/ void ctcubemask::save_counts_map(const GCTAObservation* obs, - const std::string& outfile) const + const std::string& outfile) const { // Save only if observation is valid if (obs != NULL) { diff --git a/src/ctcubemask/ctcubemask.hpp b/src/ctcubemask/ctcubemask.hpp index 196c98d..2561263 100644 --- a/src/ctcubemask/ctcubemask.hpp +++ b/src/ctcubemask/ctcubemask.hpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctcubemask - CTA cube filter tool * + * ctcubemask.hpp - CTA cube filter tool * * ----------------------------------------------------------------------- * * copyright (C) 2014 by Chia-Chun Lu * * ----------------------------------------------------------------------- * @@ -35,7 +35,7 @@ /* __Definitions _________________________________________________________ */ #define CTCUBEMASK_NAME "ctcubemask" -#define CTCUBEMASK_VERSION "00-05-00" +#define CTCUBEMASK_VERSION "00-01-00" /***********************************************************************//** @@ -53,7 +53,7 @@ public: virtual ~ctcubemask(void); // Operators - ctcubemask& operator= (const ctcubemask& app); + ctcubemask& operator=(const ctcubemask& app); // Methods void clear(void); @@ -62,7 +62,7 @@ public: void save(void); GObservations& obs(void) { return m_obs; } void get_parameters(void); - void apply_mask(GCTAObservation* obs, const std::string& filename); + void apply_mask(GCTAObservation* obs); protected: // Protected methods diff --git a/src/ctcubemask/ctcubemask.par b/src/ctcubemask/ctcubemask.par index f9c8822..0452e24 100644 --- a/src/ctcubemask/ctcubemask.par +++ b/src/ctcubemask/ctcubemask.par @@ -4,23 +4,23 @@ # # # ------------------------------------------------------------------------- # # # -# File: ctfilter.par # -# Version: 00-04-00 # +# File: ctcubemask.par # +# Version: 00-01-00 # # # # Author: Chia-Chun Lu # # c.lu@mpi-hd.mpg.de # # IRAP # # # -# Purpose: Parameter file for the ctselect application # +# Purpose: Parameter file for the ctcubemask application # # # ############################################################################# # # File information #================= -infile, f, a, "test/data/crab_events.fits.gz",,,"Input event list or observation definition file" -regfile, f, a, "exlusion.reg",,,"exclusion region file in ds9 format" -outfile, f, a, "filtered_cube.fits",,,"Output event list or observation definition file" -prefix, s, h, "filtered_",,,"Prefix for counts maps in observation definition file" +infile, f, a, "cntmap.fits",,,"Input counts cube or observation definition file" +regfile, f, a, "NONE",,,"Exclusion region file in ds9 format" +outfile, f, a, "filtered_cube.fits",,,"Output event cube or observation definition file" +prefix, s, h, "filtered_",,,"Prefix for counts cube in observation definition file" # # Selection parameters #===================== diff --git a/test/data/exclusion.reg b/test/data/exclusion.reg new file mode 100644 index 0000000..80eaed4 --- /dev/null +++ b/test/data/exclusion.reg @@ -0,0 +1,4 @@ +# Region file format: DS9 version 4.1 +global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1 +fk5 +circle(83.63, 21.5, 0.4) diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 72bb3f8..a98ac73 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -41,6 +41,7 @@ ctobssim=../src/ctobssim/ctobssim ctselect=../src/ctselect/ctselect ctskymap=../src/ctskymap/ctskymap ctexpcube=../src/ctexpcube/ctexpcube +ctcubemask=../src/ctcubemask/ctcubemask # @@ -387,3 +388,25 @@ else $ECHO " expcube2.fits file is not valid" exit 1 fi + + +# +# Test ctcubemask +# =============== +$ECHO -n "Test ctcubemask: " +$ctcubemask infile="data/crab_events.fits" \ + regfile="data/exclusion.reg" \ + outfile="filtered_cube.fits" \ + ra=83.63 \ + dec=22.01 \ + rad=2.0 \ + emin=0.1 \ + emax=100.0 +$ECHO -n "." +if [ -s "filtered_cube.fits" ] +then + $ECHO -n "." +else + $ECHO " filtered_cube.fits file is not valid" + exit 1 +fi -- libgit2 0.21.4