Luigi Tibaldo / ctools

Commit baa7cd9c734504fcd1c5260070d4516e950ad5ef

Authored by Jürgen Knödlseder
2 parents b6fde3a7 b130ceb3

Merge branch 'bugfix-1.7.1' into integration 

# Conflicts:
#	ChangeLog
#	NEWS
#	README.md
Inline Side-by-side
Showing 12 changed files with 104 additions and 81 deletions
.gitignore
  View file @baa7cd9
... ... @@ -54,6 +54,7 @@ pyext/build/
54 54 pyext/setup.py
55 55 pyext/ctools/__init__.py
56 56 pyext/ctools/tools.py
  57 +pyext/ctools/_tools.so
57 58 pyext/ctools/tools_wrap.cpp
58 59 cscripts/__init__.py
59 60 src/ctbin/ctbin
... ...
ChangeLog
  View file @baa7cd9
... ... @@ -4,6 +4,16 @@
4 4 ========================
5 5  
6 6  
  7 +2020-08-20
  8 +
  9 + * Bug fix version 1.7.1 released
  10 + ================================
  11 +
  12 + Copy Python module into pyext/ctools folder (#3322)
  13 + Add units to reference manual of response cube tools (#3321)
  14 + Speed-up ctbin for stacked event binning (#3275)
  15 +
  16 +
7 17 2020-07-15
8 18  
9 19 * Version 1.7.0 released
... ...
configure.ac
  View file @baa7cd9
1 1 #############################################################################
2 2 # Configuration for ctools - tools for science analysis of CTA data #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2011-2019 Juergen Knoedlseder #
  4 +# Copyright (C) 2011-2020 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ...
doc/source/admin/release_history/1.7.1.rst 0 → 100644
  View file @baa7cd9
  1 +.. _1.7.1:
  2 +
  3 +ctools 1.7.1 release (20 August 2020)
  4 +=====================================
  5 +
  6 +Overview
  7 +--------
  8 +
  9 +ctools 1.7.1 is a bug fix release for ctools release 1.7.
  10 +
  11 +The following bugs were fixed:
  12 +
  13 +* [`3322 <https://cta-redmine.irap.omp.eu/issues/3322>`_] -
  14 + Copy Python module into pyext/ctools folder
  15 +* [`3321 <https://cta-redmine.irap.omp.eu/issues/3321>`_] -
  16 + Add units to reference manual of response cube tools
  17 +
  18 +The following improvements were done:
  19 +
  20 +* [`3275 <https://cta-redmine.irap.omp.eu/issues/3275>`_] -
  21 + Speed-up ctbin for stacked event binning
... ...
doc/source/admin/release_history/index.rst
  View file @baa7cd9
... ... @@ -6,6 +6,8 @@ Release History
6 6 .. toctree::
7 7 :maxdepth: 1
8 8  
  9 + ctools 1.8 <1.8>
  10 + ctools 1.7.1 <1.7.1>
9 11 ctools 1.7 <1.7>
10 12 ctools 1.6.3 <1.6.3>
11 13 ctools 1.6.2 <1.6.2>
... ...
doc/source/users/reference_manual/ctbkgcube.rst
  View file @baa7cd9
... ... @@ -24,8 +24,9 @@ the ``incube`` parameter.
24 24  
25 25 :ref:`ctbkgcube` generates a background cube FITS file comprising two extensions.
26 26 The primary extension contains a 3-dimensional image that contains the
27   -background cube values. The next extension named ``EBOUNDS`` contains a
28   -binary table that defines the energy boundaries of the background cube.
  27 +background cube values. The values are in units of MeV**(-1) s**(-1) sr**(-1).
  28 +The next extension named ``EBOUNDS`` contains a binary table that defines the
  29 +energy boundaries of the background cube.
29 30  
30 31 :ref:`ctbkgcube` generates also an output model XML file that can serve as input
31 32 for a maximum likelihood analysis. The output model XML file is a copy of
... ...
doc/source/users/reference_manual/ctedispcube.rst
  View file @baa7cd9
... ... @@ -30,11 +30,11 @@ value of ``binsz=1.0``).
30 30  
31 31 :ref:`ctedispcube` generates an energy dispersion cube FITS file comprising three
32 32 extensions. The primary extension contains a 4-dimensional image that contains
33   -the energy disperison values. The next extension named ``ENERGIES`` contains
34   -a binary table that defines the energies of the energy dispersion cube. The
35   -last extension named ``MIGRAS`` contains a binary table that defines the
36   -migration values of the cube which is the ratio between reconstructed and true
37   -photon energy.
  33 +the energy disperison values. The values are in units of MeV**(-1). The next
  34 +extension named ``ENERGIES`` contains a binary table that defines the energies of
  35 +the energy dispersion cube. The last extension named ``MIGRAS`` contains a binary
  36 +table that defines the migration values of the cube which is the ratio between
  37 +reconstructed and true photon energy.
38 38  
39 39  
40 40 General parameters
... ...
doc/source/users/reference_manual/ctexpcube.rst
  View file @baa7cd9
... ... @@ -38,10 +38,11 @@ well sampled (in particular at low energies).
38 38  
39 39 :ref:`ctexpcube` generates an exposure cube FITS file comprising three extensions.
40 40 The primary extension contains a 3-dimensional image that contains the
41   -exposure values. The next extension named ``ENERGIES`` contains a binary table
42   -that defines the energies of the exposure cube. The last extension named ``GTI``
43   -contains a binary table that provides the Good Time Intervals of all
44   -observations that have been used for the computation of the exposure cube.
  41 +exposure values. The values are in units of cm**2 s. The next extension named
  42 +``ENERGIES`` contains a binary table that defines the energies of the exposure
  43 +cube. The last extension named ``GTI`` contains a binary table that provides the
  44 +Good Time Intervals of all observations that have been used for the computation
  45 +of the exposure cube.
45 46  
46 47  
47 48 General parameters
... ...
doc/source/users/reference_manual/ctpsfcube.rst
  View file @baa7cd9
... ... @@ -30,10 +30,11 @@ manageable (with a typical value of ``binsz=1.0``).
30 30  
31 31 :ref:`ctpsfcube` generates a point spread function cube FITS file comprising three
32 32 extensions. The primary extension contains a 4-dimensional image that contains
33   -the point spread function values. The next extension named ``ENERGIES`` contains
34   -a binary table that defines the energies of the point spread function cube.
35   -The last extension named ``DELTAS`` contains a binary table that defines the
36   -offset angles between true and measured arrival direction of the photon.
  33 +the point spread function values. The values are in units of sr**(-1). The next
  34 +extension named ``ENERGIES`` contains a binary table that defines the energies of
  35 +the point spread function cube. The last extension named ``DELTAS`` contains a
  36 +binary table that defines the offset angles between true and measured arrival
  37 +direction of the photon.
37 38  
38 39  
39 40 General parameters
... ...
pyext/Makefile.am
  View file @baa7cd9
1 1 #############################################################################
2 2 # Makefile for ctools Python extension module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2011-2017 Juergen Knoedlseder #
  4 +# Copyright (C) 2011-2020 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -160,6 +160,11 @@ if IS_ELCAPITAN
160 160 install_name_tool -add_rpath ${GAMMALIB}/lib $$f; \
161 161 done
162 162 endif
  163 + -for f in build/ctools/ctools/_*.so; do \
  164 + dst=`echo $$f | sed 's/build\/ctools\///g'`; \
  165 + cp $$f $$dst; \
  166 + echo "Copied "$$f" -> "$$dst; \
  167 + done
163 168  
164 169 # Install the ctools Python module
165 170 install-exec-local:
... ...
src/ctbin/ctbin.cpp
  View file @baa7cd9
1 1 /***************************************************************************
2 2 * ctbin - Event binning tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -89,7 +89,6 @@ ctbin::ctbin(const GObservations&amp; obs) :
89 89 }
90 90  
91 91  
92   -
93 92 /***********************************************************************//**
94 93 * @brief Command line constructor
95 94 *
... ... @@ -221,9 +220,9 @@ void ctbin::clear(void)
221 220 * @brief Run the event binning tool
222 221 *
223 222 * Gets the user parameters and loops over all CTA observations in the
224   - * observation container to bin the events into a single counts cube. All
225   - * observations in the observation container that do not contain CTA event
226   - * lists will be skipped.
  223 + * observation container to bin the events into either a single or a set of
  224 + * counts cube. All observations in the observation container that do not
  225 + * contain CTA event lists will be skipped.
227 226 ***************************************************************************/
228 227 void ctbin::run(void)
229 228 {
... ... @@ -281,9 +280,15 @@ void ctbin::run(void)
281 280  
282 281 // Initialise arrays
283 282 if (nobs > 0) {
284   - m_counts = std::vector<GSkyMap>(nobs);
285   - m_weights = std::vector<GSkyMap>(nobs);
286   - m_cubes = std::vector<GCTAEventCube>(nobs);
  283 +
  284 + // Determine number of cubes
  285 + int ncubes = (m_stack) ? 1 : nobs;
  286 +
  287 + // Allocate cubes
  288 + m_counts = std::vector<GSkyMap>(ncubes);
  289 + m_weights = std::vector<GSkyMap>(ncubes);
  290 + m_cubes = std::vector<GCTAEventCube>(ncubes);
  291 +
287 292 }
288 293 else {
289 294 m_counts.clear();
... ... @@ -316,46 +321,25 @@ void ctbin::run(void)
316 321 // Get pointer to observation
317 322 GCTAObservation* obs = obs_list[i];
318 323  
319   - // Fill the cube
320   - GSkyMap counts = fill_cube(obs);
  324 + // Determine counts cube slot
  325 + int islot = (m_stack) ? 0 : i;
  326 +
  327 + // If slot is empty then allocate cubes
  328 + #pragma omp critical(ctbin_run)
  329 + if (m_counts[islot].is_empty()) {
  330 + m_counts[islot] = create_cube(obs);
  331 + m_weights[islot] = create_cube(obs);
  332 + }
  333 +
  334 + // Fill counts cube
  335 + fill_cube(obs, m_counts[islot]);
321 336  
322 337 // Set the counts cube weights
323   - GSkyMap weights = set_weights(obs);
  338 + set_weights(obs, m_weights[islot]);
324 339  
325 340 // Dispose events to free memory
326 341 obs->dispose_events();
327 342  
328   - // Post process cubes
329   - #pragma omp critical(ctbin_fill_cube)
330   - {
331   - // If stacking of cubes is requested then add up counts cube
332   - // and set weights in the first slot of the m_counts and
333   - // m_weights arrays
334   - if (m_stack) {
335   - if (m_counts[0].is_empty()) {
336   - m_counts[0] = counts;
337   - m_weights[0] = weights;
338   - }
339   - else {
340   - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
341   - for (int pixel = 0; pixel < counts.npix(); ++pixel) {
342   - m_counts[0](pixel, iebin) += counts(pixel, iebin);
343   - if (weights(pixel, iebin) == 1.0) {
344   - m_weights[0](pixel, iebin) = 1.0;
345   - }
346   - }
347   - }
348   - }
349   - }
350   -
351   - // ... otherwise store counts and weights
352   - else {
353   - m_counts[i] = counts;
354   - m_weights[i] = weights;
355   - }
356   -
357   - } // end: omp critical section
358   -
359 343 } // endfor: looped over observations
360 344  
361 345 // Combine GTIs of all observations and compute total ontime and livetime
... ... @@ -369,16 +353,16 @@ void ctbin::run(void)
369 353  
370 354 // If resulting GTI is not filled then copy over the reference from
371 355 // the observation and use it for the resulting GTI
372   - if (m_gti.is_empty()) {
  356 + if (m_gti.is_empty()) {
373 357 m_gti.reference(gti.reference());
374   - }
  358 + }
375 359  
376   - // Append GTIs
377   - m_gti.extend(gti);
  360 + // Append GTIs
  361 + m_gti.extend(gti);
378 362  
379   - // Update ontime and livetime
380   - m_ontime += obs->ontime();
381   - m_livetime += obs->livetime();
  363 + // Update ontime and livetime
  364 + m_ontime += obs->ontime();
  365 + m_livetime += obs->livetime();
382 366  
383 367 } // endfor: looped over all observations
384 368  
... ... @@ -842,13 +826,14 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs)
842 826 * @brief Fill events into counts cube
843 827 *
844 828 * @param[in] obs CTA observation.
  829 + * @param[in,out] counts Counts cube
845 830 *
846 831 * @exception GException::invalid_value
847 832 * No event list or valid RoI found in observation.
848 833 *
849 834 * Fills the events from an event list into the counts cube.
850 835 ***************************************************************************/
851   -GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
  836 +void ctbin::fill_cube(const GCTAObservation* obs, GSkyMap& counts)
852 837 {
853 838 // Make sure that the observation holds a CTA event list. If this is
854 839 // not the case then throw an exception.
... ... @@ -883,9 +868,6 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
883 868 int num_outside_ebds = 0;
884 869 int num_in_map = 0;
885 870  
886   - // Create counts cube
887   - GSkyMap counts = create_cube(obs);
888   -
889 871 // Extract counts cube boundaries
890 872 double pixel_x_max = counts.nx() - 0.5;
891 873 double pixel_y_max = counts.ny() - 0.5;
... ... @@ -936,6 +918,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
936 918 }
937 919  
938 920 // Fill event in skymap
  921 + #pragma omp atomic
939 922 counts(pixel, iebin) += 1.0;
940 923  
941 924 // Increment number of maps
... ... @@ -956,7 +939,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
956 939 }
957 940  
958 941 // Return
959   - return counts;
  942 + return;
960 943 }
961 944  
962 945  
... ... @@ -964,6 +947,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
964 947 * @brief Set counts cube weights for a given observation
965 948 *
966 949 * @param[in] obs CTA observation.
  950 + * @param[in,out] weights Counts cube weights.
967 951 *
968 952 * @exception GException::invalid_value
969 953 * No event list or valid RoI found in observation.
... ... @@ -971,7 +955,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
971 955 * Sets the counts cube weights for all bins that are considered for the
972 956 * specific observation to unity.
973 957 ***************************************************************************/
974   -GSkyMap ctbin::set_weights(const GCTAObservation* obs)
  958 +void ctbin::set_weights(const GCTAObservation* obs, GSkyMap& weights)
975 959 {
976 960 // Get the RoI
977 961 const GCTARoi& roi = obs->roi();
... ... @@ -992,9 +976,6 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs)
992 976 // Get counts cube usage flags
993 977 std::vector<bool> usage = cube_layer_usage(m_ebounds, obs->ebounds());
994 978  
995   - // Create weights cube
996   - GSkyMap weights = create_cube(obs);
997   -
998 979 // Set sky direction cache availability flag
999 980 bool cache = (m_stack && m_dirs.size() == weights.npix());
1000 981  
... ... @@ -1029,7 +1010,7 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs)
1029 1010 } // endfor: looped over pixels of counts cube
1030 1011  
1031 1012 // Return
1032   - return weights;
  1013 + return;
1033 1014 }
1034 1015  
1035 1016  
... ...
src/ctbin/ctbin.hpp
  View file @baa7cd9
1 1 /***************************************************************************
2 2 * ctbin - Event binning tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -44,8 +44,8 @@
44 44 * FITS files, and on observation definition XML files.
45 45 *
46 46 * If multiple event lists are specified in the observation container or the
47   - * XML definition file, the class will merge these events into a single
48   - * counts cube.
  47 + * XML definition file, the class will upon request stack these events into a
  48 + * single counts cube.
49 49 *
50 50 * If the hidden parameter stack=no is used, one counts cube is generated for
51 51 * each individual observation. These cubes are saved to a path that can be
... ... @@ -94,8 +94,8 @@ protected:
94 94 void get_parameters(void);
95 95 void init_sky_dir_cache(void);
96 96 GSkyMap create_cube(const GCTAObservation* obs);
97   - GSkyMap fill_cube(const GCTAObservation* obs);
98   - GSkyMap set_weights(const GCTAObservation* obs);
  97 + void fill_cube(const GCTAObservation* obs, GSkyMap& counts);
  98 + void set_weights(const GCTAObservation* obs, GSkyMap& weights);
99 99 void obs_cube_stacked(void);
100 100  
101 101 // User parameters
... ...