Iftach Sadeh / ctools

Commit d33e30aabadddd8d73c13b7d0d519599442d5673

Authored by Jürgen Knödlseder
2 parents 61049399 e607c995

Merge branch 'faster-ctbin' into bugfix-1.7.1

Inline Side-by-side
Showing 2 changed files with 45 additions and 64 deletions
src/ctbin/ctbin.cpp
  View file @d33e30a
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& 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 @d33e30a
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
... ...