diff --git a/.gitignore b/.gitignore index 2224049..ac6551d 100755 --- a/.gitignore +++ b/.gitignore @@ -54,6 +54,7 @@ pyext/build/ pyext/setup.py pyext/ctools/__init__.py pyext/ctools/tools.py +pyext/ctools/_tools.so pyext/ctools/tools_wrap.cpp cscripts/__init__.py src/ctbin/ctbin diff --git a/ChangeLog b/ChangeLog index d8b1d4f..c201d8a 100755 --- a/ChangeLog +++ b/ChangeLog @@ -4,6 +4,16 @@ ======================== +2020-08-20 + + * Bug fix version 1.7.1 released + ================================ + + Copy Python module into pyext/ctools folder (#3322) + Add units to reference manual of response cube tools (#3321) + Speed-up ctbin for stacked event binning (#3275) + + 2020-07-15 * Version 1.7.0 released diff --git a/configure.ac b/configure.ac index 3c20906..2700c3c 100755 --- a/configure.ac +++ b/configure.ac @@ -1,7 +1,7 @@ ############################################################################# # Configuration for ctools - tools for science analysis of CTA data # # ------------------------------------------------------------------------- # -# Copyright (C) 2011-2019 Juergen Knoedlseder # +# Copyright (C) 2011-2020 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # diff --git a/doc/source/admin/release_history/1.7.1.rst b/doc/source/admin/release_history/1.7.1.rst new file mode 100644 index 0000000..2b335ba --- /dev/null +++ b/doc/source/admin/release_history/1.7.1.rst @@ -0,0 +1,21 @@ +.. _1.7.1: + +ctools 1.7.1 release (20 August 2020) +===================================== + +Overview +-------- + +ctools 1.7.1 is a bug fix release for ctools release 1.7. + +The following bugs were fixed: + +* [`3322 `_] - + Copy Python module into pyext/ctools folder +* [`3321 `_] - + Add units to reference manual of response cube tools + +The following improvements were done: + +* [`3275 `_] - + Speed-up ctbin for stacked event binning diff --git a/doc/source/admin/release_history/index.rst b/doc/source/admin/release_history/index.rst index 9446570..ceaf817 100755 --- a/doc/source/admin/release_history/index.rst +++ b/doc/source/admin/release_history/index.rst @@ -6,6 +6,8 @@ Release History .. toctree:: :maxdepth: 1 + ctools 1.8 <1.8> + ctools 1.7.1 <1.7.1> ctools 1.7 <1.7> ctools 1.6.3 <1.6.3> ctools 1.6.2 <1.6.2> diff --git a/doc/source/users/reference_manual/ctbkgcube.rst b/doc/source/users/reference_manual/ctbkgcube.rst index 9ef5c46..5956ffe 100755 --- a/doc/source/users/reference_manual/ctbkgcube.rst +++ b/doc/source/users/reference_manual/ctbkgcube.rst @@ -24,8 +24,9 @@ the ``incube`` parameter. :ref:`ctbkgcube` generates a background cube FITS file comprising two extensions. The primary extension contains a 3-dimensional image that contains the -background cube values. The next extension named ``EBOUNDS`` contains a -binary table that defines the energy boundaries of the background cube. +background cube values. The values are in units of MeV**(-1) s**(-1) sr**(-1). +The next extension named ``EBOUNDS`` contains a binary table that defines the +energy boundaries of the background cube. :ref:`ctbkgcube` generates also an output model XML file that can serve as input for a maximum likelihood analysis. The output model XML file is a copy of diff --git a/doc/source/users/reference_manual/ctedispcube.rst b/doc/source/users/reference_manual/ctedispcube.rst index 38eb884..6959229 100755 --- a/doc/source/users/reference_manual/ctedispcube.rst +++ b/doc/source/users/reference_manual/ctedispcube.rst @@ -30,11 +30,11 @@ value of ``binsz=1.0``). :ref:`ctedispcube` generates an energy dispersion cube FITS file comprising three extensions. The primary extension contains a 4-dimensional image that contains -the energy disperison values. The next extension named ``ENERGIES`` contains -a binary table that defines the energies of the energy dispersion cube. The -last extension named ``MIGRAS`` contains a binary table that defines the -migration values of the cube which is the ratio between reconstructed and true -photon energy. +the energy disperison values. The values are in units of MeV**(-1). The next +extension named ``ENERGIES`` contains a binary table that defines the energies of +the energy dispersion cube. The last extension named ``MIGRAS`` contains a binary +table that defines the migration values of the cube which is the ratio between +reconstructed and true photon energy. General parameters diff --git a/doc/source/users/reference_manual/ctexpcube.rst b/doc/source/users/reference_manual/ctexpcube.rst index 28859da..0f01772 100755 --- a/doc/source/users/reference_manual/ctexpcube.rst +++ b/doc/source/users/reference_manual/ctexpcube.rst @@ -38,10 +38,11 @@ well sampled (in particular at low energies). :ref:`ctexpcube` generates an exposure cube FITS file comprising three extensions. The primary extension contains a 3-dimensional image that contains the -exposure values. The next extension named ``ENERGIES`` contains a binary table -that defines the energies of the exposure cube. The last extension named ``GTI`` -contains a binary table that provides the Good Time Intervals of all -observations that have been used for the computation of the exposure cube. +exposure values. The values are in units of cm**2 s. The next extension named +``ENERGIES`` contains a binary table that defines the energies of the exposure +cube. The last extension named ``GTI`` contains a binary table that provides the +Good Time Intervals of all observations that have been used for the computation +of the exposure cube. General parameters diff --git a/doc/source/users/reference_manual/ctpsfcube.rst b/doc/source/users/reference_manual/ctpsfcube.rst index 0e20da0..f6cad60 100755 --- a/doc/source/users/reference_manual/ctpsfcube.rst +++ b/doc/source/users/reference_manual/ctpsfcube.rst @@ -30,10 +30,11 @@ manageable (with a typical value of ``binsz=1.0``). :ref:`ctpsfcube` generates a point spread function cube FITS file comprising three extensions. The primary extension contains a 4-dimensional image that contains -the point spread function values. The next extension named ``ENERGIES`` contains -a binary table that defines the energies of the point spread function cube. -The last extension named ``DELTAS`` contains a binary table that defines the -offset angles between true and measured arrival direction of the photon. +the point spread function values. The values are in units of sr**(-1). The next +extension named ``ENERGIES`` contains a binary table that defines the energies of +the point spread function cube. The last extension named ``DELTAS`` contains a +binary table that defines the offset angles between true and measured arrival +direction of the photon. General parameters diff --git a/pyext/Makefile.am b/pyext/Makefile.am index b576b13..351ef2c 100755 --- a/pyext/Makefile.am +++ b/pyext/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile for ctools Python extension module # # ------------------------------------------------------------------------- # -# Copyright (C) 2011-2017 Juergen Knoedlseder # +# Copyright (C) 2011-2020 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -160,6 +160,11 @@ if IS_ELCAPITAN install_name_tool -add_rpath ${GAMMALIB}/lib $$f; \ done endif + -for f in build/ctools/ctools/_*.so; do \ + dst=`echo $$f | sed 's/build\/ctools\///g'`; \ + cp $$f $$dst; \ + echo "Copied "$$f" -> "$$dst; \ + done # Install the ctools Python module install-exec-local: diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index b703f04..a41108e 100755 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2018 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -89,7 +89,6 @@ ctbin::ctbin(const GObservations& obs) : } - /***********************************************************************//** * @brief Command line constructor * @@ -221,9 +220,9 @@ void ctbin::clear(void) * @brief Run the event binning tool * * Gets the user parameters and loops over all CTA observations in the - * observation container to bin the events into a single counts cube. All - * observations in the observation container that do not contain CTA event - * lists will be skipped. + * observation container to bin the events into either a single or a set of + * counts cube. All observations in the observation container that do not + * contain CTA event lists will be skipped. ***************************************************************************/ void ctbin::run(void) { @@ -281,9 +280,15 @@ void ctbin::run(void) // Initialise arrays if (nobs > 0) { - m_counts = std::vector(nobs); - m_weights = std::vector(nobs); - m_cubes = std::vector(nobs); + + // Determine number of cubes + int ncubes = (m_stack) ? 1 : nobs; + + // Allocate cubes + m_counts = std::vector(ncubes); + m_weights = std::vector(ncubes); + m_cubes = std::vector(ncubes); + } else { m_counts.clear(); @@ -316,46 +321,25 @@ void ctbin::run(void) // Get pointer to observation GCTAObservation* obs = obs_list[i]; - // Fill the cube - GSkyMap counts = fill_cube(obs); + // Determine counts cube slot + int islot = (m_stack) ? 0 : i; + + // If slot is empty then allocate cubes + #pragma omp critical(ctbin_run) + if (m_counts[islot].is_empty()) { + m_counts[islot] = create_cube(obs); + m_weights[islot] = create_cube(obs); + } + + // Fill counts cube + fill_cube(obs, m_counts[islot]); // Set the counts cube weights - GSkyMap weights = set_weights(obs); + set_weights(obs, m_weights[islot]); // Dispose events to free memory obs->dispose_events(); - // Post process cubes - #pragma omp critical(ctbin_fill_cube) - { - // If stacking of cubes is requested then add up counts cube - // and set weights in the first slot of the m_counts and - // m_weights arrays - if (m_stack) { - if (m_counts[0].is_empty()) { - m_counts[0] = counts; - m_weights[0] = weights; - } - else { - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - for (int pixel = 0; pixel < counts.npix(); ++pixel) { - m_counts[0](pixel, iebin) += counts(pixel, iebin); - if (weights(pixel, iebin) == 1.0) { - m_weights[0](pixel, iebin) = 1.0; - } - } - } - } - } - - // ... otherwise store counts and weights - else { - m_counts[i] = counts; - m_weights[i] = weights; - } - - } // end: omp critical section - } // endfor: looped over observations // Combine GTIs of all observations and compute total ontime and livetime @@ -369,16 +353,16 @@ void ctbin::run(void) // If resulting GTI is not filled then copy over the reference from // the observation and use it for the resulting GTI - if (m_gti.is_empty()) { + if (m_gti.is_empty()) { m_gti.reference(gti.reference()); - } + } - // Append GTIs - m_gti.extend(gti); + // Append GTIs + m_gti.extend(gti); - // Update ontime and livetime - m_ontime += obs->ontime(); - m_livetime += obs->livetime(); + // Update ontime and livetime + m_ontime += obs->ontime(); + m_livetime += obs->livetime(); } // endfor: looped over all observations @@ -842,13 +826,14 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs) * @brief Fill events into counts cube * * @param[in] obs CTA observation. + * @param[in,out] counts Counts cube * * @exception GException::invalid_value * No event list or valid RoI found in observation. * * Fills the events from an event list into the counts cube. ***************************************************************************/ -GSkyMap ctbin::fill_cube(const GCTAObservation* obs) +void ctbin::fill_cube(const GCTAObservation* obs, GSkyMap& counts) { // Make sure that the observation holds a CTA event list. If this is // not the case then throw an exception. @@ -883,9 +868,6 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) int num_outside_ebds = 0; int num_in_map = 0; - // Create counts cube - GSkyMap counts = create_cube(obs); - // Extract counts cube boundaries double pixel_x_max = counts.nx() - 0.5; double pixel_y_max = counts.ny() - 0.5; @@ -936,6 +918,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) } // Fill event in skymap + #pragma omp atomic counts(pixel, iebin) += 1.0; // Increment number of maps @@ -956,7 +939,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) } // Return - return counts; + return; } @@ -964,6 +947,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) * @brief Set counts cube weights for a given observation * * @param[in] obs CTA observation. + * @param[in,out] weights Counts cube weights. * * @exception GException::invalid_value * No event list or valid RoI found in observation. @@ -971,7 +955,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) * Sets the counts cube weights for all bins that are considered for the * specific observation to unity. ***************************************************************************/ -GSkyMap ctbin::set_weights(const GCTAObservation* obs) +void ctbin::set_weights(const GCTAObservation* obs, GSkyMap& weights) { // Get the RoI const GCTARoi& roi = obs->roi(); @@ -992,9 +976,6 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs) // Get counts cube usage flags std::vector usage = cube_layer_usage(m_ebounds, obs->ebounds()); - // Create weights cube - GSkyMap weights = create_cube(obs); - // Set sky direction cache availability flag bool cache = (m_stack && m_dirs.size() == weights.npix()); @@ -1029,7 +1010,7 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs) } // endfor: looped over pixels of counts cube // Return - return weights; + return; } diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index e60f01a..063c32e 100755 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2018 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -44,8 +44,8 @@ * FITS files, and on observation definition XML files. * * If multiple event lists are specified in the observation container or the - * XML definition file, the class will merge these events into a single - * counts cube. + * XML definition file, the class will upon request stack these events into a + * single counts cube. * * If the hidden parameter stack=no is used, one counts cube is generated for * each individual observation. These cubes are saved to a path that can be @@ -94,8 +94,8 @@ protected: void get_parameters(void); void init_sky_dir_cache(void); GSkyMap create_cube(const GCTAObservation* obs); - GSkyMap fill_cube(const GCTAObservation* obs); - GSkyMap set_weights(const GCTAObservation* obs); + void fill_cube(const GCTAObservation* obs, GSkyMap& counts); + void set_weights(const GCTAObservation* obs, GSkyMap& weights); void obs_cube_stacked(void); // User parameters