Pierrick Martin / ctools

Commit 433c9b7b999534f0e2645840de390ff4ffffde95

Authored by Jürgen Knödlseder
2 parents 02f23c95 e3263b75

Merge branch '2714-joint-binned' into integration 

# Conflicts:
#	ChangeLog
#	doc/source/admin/release_history/1.6.rst
Inline Side-by-side
Showing 12 changed files with 575 additions and 174 deletions
ChangeLog
  View file @433c9b7
... ... @@ -4,6 +4,7 @@
4 4 ========================
5 5  
6 6 Apply event bin weighting in ctmodel for stacked computation (#2721)
  7 + Enable ctbin to produce counts cubes for a joint binned analysis (#2714)
7 8 Add Jupyter notebooks for H.E.S.S. DR1 (#2712)
8 9 Add code generator (#2709)
9 10 Add csbkgmodel script (#2694)
... ...
... ... @@ -61,6 +61,14 @@ of the filename is now done by the GApplicationPar class.
61 61  
62 62 ctbin - CTA event binning
63 63 -------------------------
  64 +Added the option to produce an individual counts cube for each event list in
  65 +the input observation definition XML file. The option is enabled by setting
  66 +the new hidden parameter "stack" to "no". The cube() method now has an optional
  67 +index argument, providing access to the i'th counts cube that was produced.
  68 +The cubes() method specifies the number of counts cubes that were produced.
  69 +If no counts cube was produced, the cube() method now throws an exception
  70 +(#2714).
  71 +
64 72 ctbin now accouts for varying energy thresholds when computing the weight
65 73 array (#2685).
66 74  
... ...
cscripts/Makefile.am
  View file @433c9b7
... ... @@ -164,6 +164,7 @@ testscripts = $(srcdir)/tests/__init__.py \
164 164 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data
165 165 testdata = $(top_srcdir)/test/data/crab.xml \
166 166 $(top_srcdir)/test/data/obs_unbinned.xml \
  167 + $(top_srcdir)/test/data/obs_unbinned_two.xml \
167 168 $(top_srcdir)/test/data/obs_stacked.xml \
168 169 $(top_srcdir)/test/data/obs_stacked_two.xml \
169 170 $(top_srcdir)/test/data/crab_events.fits \
... ...
doc/source/admin/release_history/1.6.rst
  View file @433c9b7
... ... @@ -43,6 +43,9 @@ Bug fixes
43 43 Improvements
44 44 ------------
45 45  
  46 +* [`2714 <https://cta-redmine.irap.omp.eu/issues/2714>`_] -
  47 + Enable :ref:`ctbin` to produce counts cubes for a joint binned analysis.
  48 + Rename ``outcube`` to ``outobs`` parameter.
46 49 * [`2712 <https://cta-redmine.irap.omp.eu/issues/2712>`_] -
47 50 Add Jupyter notebooks for H.E.S.S. DR1
48 51 * [`2685 <https://cta-redmine.irap.omp.eu/issues/2685>`_] -
... ...
doc/source/users/reference_manual/ctbin.rst
  View file @433c9b7
... ... @@ -14,17 +14,27 @@ a 3-dimensional data cube spanned by Right Ascension or Galactic longitude,
14 14 Declination or Galactic latitude, and energy. The energy binning may be either
15 15 linear, logarithmic, or custom defined using an input file. The events are
16 16 either taken from a single event list file or from the event lists that are
17   -specified in an observation definition file. In case that multiple event
18   -lists are given in an observation definition file, the tool will loop over
19   -all event lists and fill all events into a single counts cube.
  17 +specified in an observation definition file.
20 18  
21   -:ref:`ctbin` generates a counts cube FITS file comprising three extensions. The
  19 +In case that multiple event lists are given in an observation definition file
  20 +there are two options. By default, the tool will loop over all event lists and
  21 +fill all events into a single counts cube. If ``stack=no`` is specified, the
  22 +tool will bin each event list into a separate counts cube.
  23 +
  24 +:ref:`ctbin` generates counts cube FITS file(s) comprising three extensions. The
22 25 primary extension contains a 3-dimensional image that contains the counts
23 26 cube values. The next extension named ``EBOUNDS`` contains a binary table
24 27 that defines the energy boundaries of the counts cube. The last extension
25 28 named ``GTI`` contains a binary table that defines the Good Time Intervals
26 29 of all event lists that have been filled into the counts cube.
27 30  
  31 +By default, the tool will write out a single, stacked, counts cube. If ``stack=no``
  32 +is specified, the tool will write all counts cubes into a separate FITS file,
  33 +prefixing the input event list file name with the string that is specified by
  34 +the hidden ``prefix`` parameter. In addition, the tool will create an observation
  35 +definition XML file that lists all produced counts cubes, and that can be used
  36 +as input for any further analysis step.
  37 +
28 38  
29 39 General parameters
30 40 ------------------
... ... @@ -32,8 +42,18 @@ General parameters
32 42 ``inobs [file]``
33 43 Input event list or observation definition XML file.
34 44  
35   -``outcube [file]``
36   - Output counts cube file.
  45 +``outobs [file]``
  46 + Output counts cube file or observation definition XML file.
  47 +
  48 +``(stack = yes) [boolean]``
  49 + Specify whether stacking of the events in a single counts cube is requested.
  50 + By default this hidden parameter is set to yes, and the tool produces a
  51 + single counts cube FITS file on output. If ``no`` is specified, the tool
  52 + will produce a counts cube for each event list that is found in the input
  53 + observation, and will also write an observation definition XML file.
  54 +
  55 +``(prefix = cntcube_) [string]``
  56 + Prefix for output counts cube in observation definition XML file.
37 57  
38 58 ``ebinalg <FILE|LIN|LOG> [string]``
39 59 Algorithm for defining energy bins. For ``FILE``, the energy bins are defined
... ...
doc/source/users/reference_manual/ctselect.rst
  View file @433c9b7
... ... @@ -52,7 +52,7 @@ General parameters
52 52 Output event list or observation definition XML file.
53 53  
54 54 ``(prefix = "selected_") [string]``
55   - Prefix for output event lists in observation definition file.
  55 + Prefix for output event lists in observation definition XML file.
56 56  
57 57 ``(usepnt = yes) [boolean]``
58 58 Use pointing instead of RA/DEC parameters?
... ...
pyext/ctbin.i
  View file @433c9b7
1 1 /***************************************************************************
2 2 * ctbin - Event binning tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -46,9 +46,10 @@ public:
46 46 // Methods
47 47 void clear(void);
48 48 void run(void);
  49 + int cubes(void) const;
  50 + const GCTAEventCube& cube(const int& index = 0) const;
49 51 void save(void);
50 52 void publish(const std::string& name = "");
51   - const GCTAEventCube& cube(void) const;
52 53 };
53 54  
54 55  
... ...
src/ctbin/ctbin.cpp
  View file @433c9b7
... ... @@ -38,6 +38,7 @@
38 38 #endif
39 39  
40 40 /* __ Method name definitions ____________________________________________ */
  41 +#define G_CUBE "ctbin::cube(int&)"
41 42 #define G_GET_PARAMETERS "ctbin::get_parameters()"
42 43 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)"
43 44 #define G_SET_WEIGHTS "ctbin::set_weights(GCTAObservation*)"
... ... @@ -234,6 +235,9 @@ void ctbin::run(void)
234 235 // Get task parameters
235 236 get_parameters();
236 237  
  238 + // Initialise sky direction cache for stacking
  239 + init_sky_dir_cache();
  240 +
237 241 // Write input observation container into logger
238 242 log_observations(NORMAL, m_obs, "Input observation");
239 243  
... ... @@ -253,64 +257,171 @@ void ctbin::run(void)
253 257 std::string msg = " Including unbinned "+obs->instrument()+
254 258 " observation";
255 259 log_string(NORMAL, msg);
  260 +
  261 + } // endfor: looped over all unbinned observations
  262 +
  263 + // Set number of relevant observations
  264 + int nobs = obs_list.size();
  265 +
  266 + // Initialise arrays
  267 + if (nobs > 0) {
  268 + m_counts = std::vector<GSkyMap>(nobs);
  269 + m_weights = std::vector<GSkyMap>(nobs);
  270 + m_cubes = std::vector<GCTAEventCube>(nobs);
  271 + }
  272 + else {
  273 + m_counts.clear();
  274 + m_weights.clear();
  275 + m_cubes.clear();
256 276 }
257 277  
258 278 // Write header into logger
259   - log_header1(TERSE, gammalib::number("Bin observation", m_obs.size()));
  279 + log_header1(TERSE, gammalib::number("Bin observation", nobs));
260 280  
261 281 // Loop over all unbinned CTA observations in the container
262 282 #pragma omp parallel for
263   - for (int i = 0; i < obs_list.size(); ++i) {
  283 + for (int i = 0; i < nobs; ++i) {
264 284  
265 285 // Get pointer to observation
266 286 GCTAObservation* obs = obs_list[i];
267 287  
268 288 // Fill the cube
269   - fill_cube(obs);
  289 + GSkyMap counts = fill_cube(obs);
270 290  
271 291 // Set the counts cube weights
272   - set_weights(obs);
  292 + GSkyMap weights = set_weights(obs);
273 293  
274 294 // Dispose events to free memory
275 295 obs->dispose_events();
276 296  
  297 + // Post process cubes
  298 + #pragma omp critical(ctbin_fill_cube)
  299 + {
  300 + // If stacking of cubes is requested then add up counts cube
  301 + // and set weights in the first slot of the m_counts and
  302 + // m_weights arrays
  303 + if (m_stack) {
  304 + if (i == 0) {
  305 + m_counts[0] = counts;
  306 + m_weights[0] = weights;
  307 + }
  308 + else {
  309 + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
  310 + for (int pixel = 0; pixel < counts.npix(); ++pixel) {
  311 + m_counts[0](pixel, iebin) += counts(pixel, iebin);
  312 + if (weights(pixel, iebin) == 1.0) {
  313 + m_weights[0](pixel, iebin) = 1.0;
  314 + }
  315 + }
  316 + }
  317 + }
  318 + }
  319 +
  320 + // ... otherwise store counts and weights
  321 + else {
  322 + m_counts[i] = counts;
  323 + m_weights[i] = weights;
  324 + }
  325 +
  326 + } // end: omp critical section
  327 +
277 328 } // endfor: looped over observations
278 329  
279   - // Compute total livetime in all observations
280   - double livetime = 0.0;
281   - for (int i = 0; i < obs_list.size(); ++i) {
282   - livetime += obs_list[i]->livetime();
283   - }
  330 + // Combine GTIs of all observations and compute total ontime and livetime
  331 + for (int i = 0; i < nobs; ++i) {
  332 +
  333 + // Get pointer to observation
  334 + GCTAObservation* obs = obs_list[i];
284 335  
285   - // Compute livetime fraction per energy bin so that varying energy
286   - // thresholds are correctly taken into account in the weighting
287   - // array. Computational speed is not optimum, we have to see later
288   - // whether this is an issue.
289   - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
290   - double livetime_ebin = 0.0;
291   - for (int i = 0; i < obs_list.size(); ++i) {
292   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs_list[i]->events());
293   - std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
294   - if (usage[iebin]) {
295   - livetime_ebin += obs_list[i]->livetime();
  336 + // Get Good Time Intervals for observation
  337 + GGti gti = obs->gti();
  338 +
  339 + // If resulting GTI is not filled then copy over the reference from
  340 + // the observation and use it for the resulting GTI
  341 + if (m_gti.is_empty()) {
  342 + m_gti.reference(gti.reference());
  343 + }
  344 +
  345 + // Append GTIs
  346 + m_gti.extend(gti);
  347 +
  348 + // Update ontime and livetime
  349 + m_ontime += obs->ontime();
  350 + m_livetime += obs->livetime();
  351 +
  352 + } // endfor: looped over all observations
  353 +
  354 + // Post process data for a stacked cube
  355 + if (m_stack && nobs > 0) {
  356 +
  357 + // Compute livetime fraction per energy bin so that varying energy
  358 + // thresholds are correctly taken into account in the weighting
  359 + // array. Computational speed is not optimum, we have to see later
  360 + // whether this is an issue.
  361 + for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
  362 + double livetime_ebin = 0.0;
  363 + for (int i = 0; i < obs_list.size(); ++i) {
  364 + std::vector<bool> usage =
  365 + cube_layer_usage(m_ebounds, obs_list[i]->ebounds());
  366 + if (usage[iebin]) {
  367 + livetime_ebin += obs_list[i]->livetime();
  368 + }
296 369 }
297   - }
298   - if (livetime > 0.0) {
299   - double fraction = livetime_ebin / livetime;
300   - for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
301   - if (m_weights(pixel, iebin) > 0.0) {
302   - m_weights(pixel, iebin) *= fraction;
  370 + if (m_livetime > 0.0) {
  371 + double fraction = livetime_ebin / m_livetime;
  372 + for (int pixel = 0; pixel < m_counts[0].npix(); ++pixel) {
  373 + if (m_weights[0](pixel, iebin) > 0.0) {
  374 + m_weights[0](pixel, iebin) *= fraction;
  375 + }
303 376 }
304 377 }
305 378 }
306   - }
307 379  
308   - // Build event cube (needs to come before obs_cube() since this method
309   - // relies on correct setting of m_cube)
310   - m_cube = GCTAEventCube(m_counts, m_weights, m_ebounds, m_gti);
  380 + // Set a single cube in the observation container
  381 + obs_cube_stacked();
  382 +
  383 + } // endif: stacked cube post processing
  384 +
  385 + // Post process data for non-stacked cubes
  386 + else {
  387 +
  388 + // Loop over all relevant observations
  389 + for (int i = 0; i < nobs; ++i) {
  390 +
  391 + // Get pointer to observation
  392 + GCTAObservation* obs = obs_list[i];
  393 +
  394 + // Build counts cube
  395 + m_cubes[i] = GCTAEventCube(m_counts[i], m_weights[i], m_ebounds,
  396 + obs->gti());
  397 +
  398 + // Attach event cube in the correct slot of the observation
  399 + // container
  400 + obs->events(m_cubes[i]);
  401 +
  402 + // Construct unique filename for possible saving of the counts
  403 + // cube into FITS files. Each observation has a unique pair of
  404 + // instrument() and id() attributes, hence we construct a lower
  405 + // case file name from these attributes, making sure that all
  406 + // whitespaces are converted into "_" characters. Note that the
  407 + // id() attribute is allowed to be blank, which can only happen
  408 + // for a single observation.
  409 + std::string inst = gammalib::tolower(obs->instrument());
  410 + std::string id = gammalib::tolower(obs->id());
  411 + inst = gammalib::replace_segment(inst, " ", "_");
  412 + id = gammalib::replace_segment(id, " ", "_");
  413 + std::string filename = m_prefix + inst;
  414 + if (!id.empty()) {
  415 + filename += "_" + id;
  416 + }
  417 + filename += ".fits";
311 418  
312   - // Set a single cube in the observation container
313   - obs_cube();
  419 + // Set filename of observation
  420 + obs->eventfile(filename);
  421 +
  422 + } // endfor: looped over relevant observations
  423 +
  424 + } // endelse: non-stacked cube post processing
314 425  
315 426 // Write resulting observation container into logger
316 427 log_observations(NORMAL, m_obs, "Binned observation");
... ... @@ -326,37 +437,106 @@ void ctbin::run(void)
326 437  
327 438  
328 439 /***********************************************************************//**
  440 + * @brief Return event cube at index
  441 + *
  442 + * @param[in] index Cube index (0,...,cubes()-1).
  443 + * @return Reference to event cube.
  444 + *
  445 + * @exception GException::out_of_range
  446 + * Event cube index out of range.
  447 + *
  448 + * Returns a reference to the event cube at the given index.
  449 + ***************************************************************************/
  450 +const GCTAEventCube& ctbin::cube(const int& index) const
  451 +{
  452 + // Compile option: raise an exception if index is out of range
  453 + if (index < 0 || index >= cubes()) {
  454 + throw GException::out_of_range(G_CUBE, "Cube index", index, cubes());
  455 + }
  456 +
  457 + // Return cube
  458 + return m_cubes[index];
  459 +}
  460 +
  461 +
  462 +/***********************************************************************//**
329 463 * @brief Save counts cube
330 464 *
331   - * Saves the counts cube into a FITS file.
  465 + * Saves the counts cube. If the counts cube was stacked single FITS file
  466 + * will be written. Otherwise each generated counts cube will be written
  467 + * into a FITS file, and an observation definition XML file will be written
  468 + * out that contains an observation for each counts cube.
332 469 ***************************************************************************/
333 470 void ctbin::save(void)
334 471 {
335 472 // Write header
336 473 log_header1(TERSE, "Save counts cube");
337 474  
338   - // Save only if filename is valid and if there is at least one
339   - // observation
340   - if ((*this)["outcube"].is_valid() && m_obs.size() > 0) {
  475 + // If counts cube was stacked then save first observation
  476 + if (m_stack) {
341 477  
342   - // Get counts cube filename
343   - m_outcube = (*this)["outcube"].filename();
  478 + // Save only if filename is valid and if there is at least one
  479 + // observation
  480 + if ((*this)["outobs"].is_valid() && m_obs.size() > 0) {
344 481  
345   - // Get CTA observation from observation container
346   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  482 + // Get counts cube filename
  483 + GFilename outobs = (*this)["outobs"].filename();
347 484  
348   - // Save only if observation is valid
349   - if (obs != NULL) {
  485 + // Get CTA observation from observation container
  486 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  487 +
  488 + // Save only if observation is valid
  489 + if (obs != NULL) {
  490 +
  491 + // Log counts cube file name
  492 + log_value(NORMAL, "Counts cube file", outobs.url());
350 493  
351   - // Log counts cube file name
352   - log_value(NORMAL, "Counts cube file", m_outcube.url());
  494 + // Save cube
  495 + obs->save(outobs, clobber());
353 496  
354   - // Save cube
355   - obs->save(m_outcube, clobber());
  497 + } // endif: observation was valid
  498 +
  499 + } // endif: outobs file was valid
  500 +
  501 + } // endif: counts cube was stacked
  502 +
  503 + // ... otherwise write out all relevant counts cubes and the observation
  504 + // definition XML file
  505 + else {
  506 +
  507 + // Loop over all observations
  508 + for (int i = 0; i < m_obs.size(); ++i) {
  509 +
  510 + // Get CTA observation from observation container
  511 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
356 512  
357   - } // endif: observation was valid
  513 + // Continue only if observation is valid
  514 + if (obs != NULL) {
358 515  
359   - } // endif: outcube file was valid
  516 + // Retrieve filename
  517 + std::string filename = obs->eventfile();
  518 +
  519 + // Save counts cube
  520 + obs->save(filename, clobber());
  521 +
  522 + } // endif: observation was valid
  523 +
  524 + } // endfor: looped over all observations
  525 +
  526 + // Write XML file
  527 + if ((*this)["outobs"].is_valid()) {
  528 +
  529 + // Get observation definition XML filename
  530 + std::string outobs = (*this)["outobs"].filename();
  531 +
  532 + // Issue warning if output filename has no .xml suffix
  533 + log_string(TERSE, warn_xml_suffix(outobs));
  534 +
  535 + // Save XML file
  536 + m_obs.save(outobs);
  537 + }
  538 +
  539 + } // endelse: counts cube were not stacked
360 540  
361 541 // Return
362 542 return;
... ... @@ -382,8 +562,10 @@ void ctbin::publish(const std::string&amp; name)
382 562 // Write counts cube name into logger
383 563 log_value(NORMAL, "Counts cube name", user_name);
384 564  
385   - // Publish counts cube
386   - m_counts.publish(user_name);
  565 + // Publish first counts cube if it exists
  566 + if (cubes() > 0) {
  567 + m_counts[0].publish(user_name);
  568 + }
387 569  
388 570 // Return
389 571 return;
... ... @@ -402,17 +584,18 @@ void ctbin::publish(const std::string&amp; name)
402 584 void ctbin::init_members(void)
403 585 {
404 586 // Initialise members
405   - m_outcube.clear();
406 587 m_usepnt = false;
  588 + m_stack = true;
  589 + m_prefix.clear();
407 590 m_publish = false;
408 591 m_chatter = static_cast<GChatter>(2);
409 592  
410 593 // Initialise protected members
  594 + m_cubes.clear();
411 595 m_counts.clear();
412 596 m_weights.clear();
413 597 m_ebounds.clear();
414 598 m_gti.clear();
415   - m_cube.clear();
416 599 m_ontime = 0.0;
417 600 m_livetime = 0.0;
418 601  
... ... @@ -432,17 +615,18 @@ void ctbin::init_members(void)
432 615 void ctbin::copy_members(const ctbin& app)
433 616 {
434 617 // Copy attributes
435   - m_outcube = app.m_outcube;
436 618 m_usepnt = app.m_usepnt;
  619 + m_stack = app.m_stack;
  620 + m_prefix = app.m_prefix;
437 621 m_publish = app.m_publish;
438 622 m_chatter = app.m_chatter;
439 623  
440 624 // Copy protected members
  625 + m_cubes = app.m_cubes;
441 626 m_counts = app.m_counts;
442 627 m_weights = app.m_weights;
443 628 m_ebounds = app.m_ebounds;
444 629 m_gti = app.m_gti;
445   - m_cube = app.m_cube;
446 630 m_ontime = app.m_ontime;
447 631 m_livetime = app.m_livetime;
448 632  
... ... @@ -479,32 +663,23 @@ void ctbin::get_parameters(void)
479 663 // information and do not accept counts cubes.
480 664 setup_observations(m_obs, false, true, false);
481 665  
482   - // Create an event cube based on task parameters
483   - GCTAEventCube cube = create_cube(m_obs);
484   -
485   - // Get the skymap from the cube and initialise all counts cube bins and
486   - // weights to zero
487   - m_counts = cube.counts();
488   - m_counts = 0.0;
489   - m_weights = m_counts;
  666 + // Create an event cube to query task parameters
  667 + GCTAEventCube cube = this->ctool::create_cube(m_obs);
490 668  
491 669 // Get energy boundaries
492 670 m_ebounds = cube.ebounds();
493 671  
  672 + // Get stack and prefix parameters
  673 + m_stack = (*this)["stack"].boolean();
  674 + m_prefix = (*this)["prefix"].string();
  675 +
494 676 // Get remaining parameters
495 677 m_publish = (*this)["publish"].boolean();
496 678 m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
497 679  
498 680 // If needed later, query output filename now
499 681 if (read_ahead()) {
500   - (*this)["outcube"].query();
501   - }
502   -
503   - // Cache sky directions in m_counts
504   - m_dirs.resize(m_counts.npix());
505   - for (int pix = 0; pix < m_counts.npix(); ++pix) {
506   - m_dirs[pix] = m_counts.inx2dir(pix);
507   - m_dirs[pix].radec(m_dirs[pix].ra(), m_dirs[pix].dec());
  682 + (*this)["outobs"].query();
508 683 }
509 684  
510 685 // Set number of OpenMP threads
... ... @@ -524,6 +699,95 @@ void ctbin::get_parameters(void)
524 699  
525 700  
526 701 /***********************************************************************//**
  702 + * @brief Initialise sky direction cache for cube stack
  703 + *
  704 + * Initialises a cache that holds the sky directions of each counts cube
  705 + * pixel for a stacked analysis. This will speed up the computations later.
  706 + *
  707 + * If no stacking is requested, the method does nothing.
  708 + ***************************************************************************/
  709 +void ctbin::init_sky_dir_cache(void)
  710 +{
  711 + // Continue only is stacking is requested
  712 + if (m_stack) {
  713 +
  714 + // Create a sky map from task parameters
  715 + GSkyMap map = this->ctool::create_map(m_obs);
  716 +
  717 + // Initialise sky direction
  718 + m_dirs.resize(map.npix());
  719 +
  720 + // Compute sky direction for each pixel
  721 + for (int pixel = 0; pixel < map.npix(); ++pixel) {
  722 + m_dirs[pixel] = map.inx2dir(pixel);
  723 + }
  724 +
  725 + } // endif: stacking was requested
  726 +
  727 + // Return
  728 + return;
  729 +}
  730 +
  731 +
  732 +/***********************************************************************//**
  733 + * @brief Create counts cube sky map for current observation
  734 + *
  735 + * @param[in] obs CTA observation.
  736 + *
  737 + * Creates a counts cube sky map for current observations.
  738 + ***************************************************************************/
  739 +GSkyMap ctbin::create_cube(const GCTAObservation* obs)
  740 +{
  741 + // Read coordinate system and projection
  742 + std::string coordsys = (*this)["coordsys"].string();
  743 + std::string proj = (*this)["proj"].string();
  744 +
  745 + // Read sky map reference
  746 + double xref = 0.0;
  747 + double yref = 0.0;
  748 + bool usepnt = (*this)["usepnt"].boolean();
  749 + if (!usepnt) {
  750 + xref = (*this)["xref"].real();
  751 + yref = (*this)["yref"].real();
  752 + }
  753 +
  754 + // Read sky map pixel size and number of bins
  755 + double binsz = (*this)["binsz"].real();
  756 + int nxpix = (*this)["nxpix"].integer();
  757 + int nypix = (*this)["nypix"].integer();
  758 +
  759 + // Read energy boundaries
  760 + GEbounds ebounds = create_ebounds();
  761 +
  762 + // If requested, get pointing from observations
  763 + if (usepnt) {
  764 +
  765 + // Get pointing direction for observation
  766 + GSkyDir pnt = obs->pointing().dir();
  767 +
  768 + // Set xref/yref based on the coordinate system
  769 + if (gammalib::toupper(coordsys) == "GAL") {
  770 + xref = pnt.l_deg();
  771 + yref = pnt.b_deg();
  772 + }
  773 + else {
  774 + xref = pnt.ra_deg();
  775 + yref = pnt.dec_deg();
  776 + }
  777 +
  778 + } // endif: got pointing from observation
  779 +
  780 + // Initialise counts cube
  781 + GSkyMap cube = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz,
  782 + nxpix, nypix, ebounds.size());
  783 +
  784 +
  785 + // Return cube
  786 + return cube;
  787 +}
  788 +
  789 +
  790 +/***********************************************************************//**
527 791 * @brief Fill events into counts cube
528 792 *
529 793 * @param[in] obs CTA observation.
... ... @@ -533,10 +797,10 @@ void ctbin::get_parameters(void)
533 797 *
534 798 * Fills the events from an event list into the counts cube.
535 799 ***************************************************************************/
536   -void ctbin::fill_cube(GCTAObservation* obs)
  800 +GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
537 801 {
538   - // Make sure that the observation holds a CTA event list. If this
539   - // is not the case then throw an exception.
  802 + // Make sure that the observation holds a CTA event list. If this is
  803 + // not the case then throw an exception.
540 804 const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
541 805 (obs->events());
542 806 if (events == NULL) {
... ... @@ -547,7 +811,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
547 811 }
548 812  
549 813 // Get the RoI
550   - const GCTARoi& roi = events->roi();
  814 + const GCTARoi& roi = obs->roi();
551 815  
552 816 // Check for RoI sanity
553 817 if (!roi.is_valid()) {
... ... @@ -559,7 +823,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
559 823 }
560 824  
561 825 // Get counts cube usage flags
562   - std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
  826 + std::vector<bool> usage = cube_layer_usage(m_ebounds, obs->ebounds());
563 827  
564 828 // Initialise binning statistics
565 829 int num_outside_roi = 0;
... ... @@ -568,6 +832,13 @@ void ctbin::fill_cube(GCTAObservation* obs)
568 832 int num_outside_ebds = 0;
569 833 int num_in_map = 0;
570 834  
  835 + // Create counts cube
  836 + GSkyMap counts = create_cube(obs);
  837 +
  838 + // Extract counts cube boundaries
  839 + double pixel_x_max = counts.nx() - 0.5;
  840 + double pixel_y_max = counts.ny() - 0.5;
  841 +
571 842 // Fill counts sky map
572 843 for (int i = 0; i < events->size(); ++i) {
573 844  
... ... @@ -587,7 +858,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
587 858 // Determine sky pixel
588 859 GSkyPixel pixel;
589 860 try {
590   - pixel = m_counts.dir2pix(dir);
  861 + pixel = counts.dir2pix(dir);
591 862 }
592 863 catch (std::exception &e) {
593 864 num_invalid_wcs++;
... ... @@ -596,8 +867,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
596 867  
597 868 // Skip event if corresponding counts cube pixel is outside the
598 869 // counts cube map range
599   - if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx()-0.5) ||
600   - pixel.y() < -0.5 || pixel.y() > (m_counts.ny()-0.5)) {
  870 + if (pixel.x() < -0.5 || pixel.x() > pixel_x_max ||
  871 + pixel.y() < -0.5 || pixel.y() > pixel_y_max) {
601 872 num_outside_map++;
602 873 continue;
603 874 }
... ... @@ -614,31 +885,13 @@ void ctbin::fill_cube(GCTAObservation* obs)
614 885 }
615 886  
616 887 // Fill event in skymap
617   - #pragma omp critical(ctbin_fill_cube)
618   - m_counts(pixel, iebin) += 1.0;
  888 + counts(pixel, iebin) += 1.0;
619 889  
620 890 // Increment number of maps
621 891 num_in_map++;
622 892  
623 893 } // endfor: looped over all events
624 894  
625   - // Update time information
626   - #pragma omp critical(ctbin_fill_cube)
627   - {
628   - // If GTI is empty then set its time reference from the observation.
629   - // From then on we keep that time reference
630   - if (m_gti.is_empty()) {
631   - m_gti.reference(events->gti().reference());
632   - }
633   -
634   - // Append GTIs
635   - m_gti.extend(events->gti());
636   -
637   - // Update ontime and livetime
638   - m_ontime += obs->ontime();
639   - m_livetime += obs->livetime();
640   - }
641   -
642 895 // Log filling results
643 896 #pragma omp critical(ctbin_fill_cube)
644 897 {
... ... @@ -652,7 +905,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
652 905 }
653 906  
654 907 // Return
655   - return;
  908 + return counts;
656 909 }
657 910  
658 911  
... ... @@ -667,21 +920,10 @@ void ctbin::fill_cube(GCTAObservation* obs)
667 920 * Sets the counts cube weights for all bins that are considered for the
668 921 * specific observation to unity.
669 922 ***************************************************************************/
670   -void ctbin::set_weights(GCTAObservation* obs)
  923 +GSkyMap ctbin::set_weights(const GCTAObservation* obs)
671 924 {
672   - // Make sure that the observation holds a CTA event list. If this
673   - // is not the case then throw an exception.
674   - const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
675   - (obs->events());
676   - if (events == NULL) {
677   - std::string msg = "CTA Observation does not contain an event "
678   - "list. An event list is needed to fill the "
679   - "counts cube.";
680   - throw GException::invalid_value(G_SET_WEIGHTS, msg);
681   - }
682   -
683 925 // Get the RoI
684   - const GCTARoi& roi = events->roi();
  926 + const GCTARoi& roi = obs->roi();
685 927  
686 928 // Check for RoI sanity
687 929 if (!roi.is_valid()) {
... ... @@ -692,15 +934,32 @@ void ctbin::set_weights(GCTAObservation* obs)
692 934 throw GException::invalid_value(G_SET_WEIGHTS, msg);
693 935 }
694 936  
  937 + // Get the RoI centre and radius in radians
  938 + GSkyDir roi_centre = roi.centre().dir();
  939 + double roi_radius = roi.radius() * gammalib::deg2rad;
  940 +
695 941 // Get counts cube usage flags
696   - std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
  942 + std::vector<bool> usage = cube_layer_usage(m_ebounds, obs->ebounds());
  943 +
  944 + // Create weights cube
  945 + GSkyMap weights = create_cube(obs);
  946 +
  947 + // Set sky direction cache availability flag
  948 + bool cache = (m_stack && m_dirs.size() == weights.npix());
697 949  
698 950 // Loop over all pixels in counts cube
699   - for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
  951 + for (int pixel = 0; pixel < weights.npix(); ++pixel) {
700 952  
701 953 // Skip pixel if it is outside the RoI
702   - if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) {
703   - continue;
  954 + if (cache) {
  955 + if (m_dirs[pixel].dist(roi_centre) > roi_radius) {
  956 + continue;
  957 + }
  958 + }
  959 + else {
  960 + if (roi_centre.dist(weights.inx2dir(pixel)) > roi_radius) {
  961 + continue;
  962 + }
704 963 }
705 964  
706 965 // Loop over all energy layers of counts cube
... ... @@ -712,15 +971,14 @@ void ctbin::set_weights(GCTAObservation* obs)
712 971 }
713 972  
714 973 // Signal that bin was filled
715   - #pragma omp critical(ctbin_set_weights)
716   - m_weights(pixel, iebin) = 1.0;
  974 + weights(pixel, iebin) = 1.0;
717 975  
718 976 } // endfor: looped over energy layers of counts cube
719 977  
720 978 } // endfor: looped over pixels of counts cube
721 979  
722 980 // Return
723   - return;
  981 + return weights;
724 982 }
725 983  
726 984  
... ... @@ -735,8 +993,11 @@ void ctbin::set_weights(GCTAObservation* obs)
735 993 * information in case that a single CTA observation is provided to support
736 994 * binned analysis.
737 995 ***************************************************************************/
738   -void ctbin::obs_cube(void)
  996 +void ctbin::obs_cube_stacked(void)
739 997 {
  998 + // Set event cube in first slot
  999 + m_cubes[0] = GCTAEventCube(m_counts[0], m_weights[0], m_ebounds, m_gti);
  1000 +
740 1001 // If we have only a single CTA observation in the container, then
741 1002 // keep that observation and just attach the event cube to it. Reset
742 1003 // the filename, otherwise we still will have the old event filename
... ... @@ -751,7 +1012,7 @@ void ctbin::obs_cube(void)
751 1012 if (obs->eventtype() == "EventList") {
752 1013  
753 1014 // Assign cube to the observation
754   - obs->events(this->cube());
  1015 + obs->events(m_cubes[0]);
755 1016  
756 1017 }
757 1018  
... ... @@ -761,7 +1022,8 @@ void ctbin::obs_cube(void)
761 1022 else {
762 1023  
763 1024 // Create empty counts cube
764   - GCTAEventCube cube(m_counts, m_weights, m_ebounds, obs->gti());
  1025 + GCTAEventCube cube(m_counts[0], m_weights[0], m_ebounds,
  1026 + obs->gti());
765 1027  
766 1028 // Assign empty cube to observation
767 1029 obs->events(cube);
... ... @@ -787,7 +1049,7 @@ void ctbin::obs_cube(void)
787 1049 GCTAObservation obs;
788 1050  
789 1051 // Attach event cube to CTA observation
790   - obs.events(this->cube());
  1052 + obs.events(m_cubes[0]);
791 1053  
792 1054 // Compute average pointing direction for all CTA event lists
793 1055 double ra = 0.0;
... ...
src/ctbin/ctbin.hpp
  View file @433c9b7
1 1 /***************************************************************************
2 2 * ctbin - Event binning tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -39,7 +39,7 @@
39 39 *
40 40 * @brief Event binning tool
41 41 *
42   - * This class bins event list(s) into a single counts cube. The class can
  42 + * This class bins event list(s) into counts cubes. The class can
43 43 * operate on predefined observation containers, on individual event list
44 44 * FITS files, and on observation definition XML files.
45 45 *
... ... @@ -47,6 +47,13 @@
47 47 * XML definition file, the class will merge these events into a single
48 48 * counts cube.
49 49 *
  50 + * If the hidden parameter stack=no is used, one counts cube is generated for
  51 + * each individual observation. These cubes are saved to a path that can be
  52 + * specified by the hidden parameter prefix, followed by the corresponding
  53 + * observation id. An observation definition XML file containing the paths
  54 + * to the newly generated counts cubes is written to path given in the
  55 + * parameter outcube.
  56 + *
50 57 * Results are stored in an observation container that can be written to disk
51 58 * in form of a single FITS file. On output, the observation container will
52 59 * have merged the input event lists into a single observation.
... ... @@ -74,51 +81,55 @@ public:
74 81 // Methods
75 82 void clear(void);
76 83 void run(void);
  84 + int cubes(void) const;
  85 + const GCTAEventCube& cube(const int& index = 0) const;
77 86 void save(void);
78 87 void publish(const std::string& name = "");
79   - const GCTAEventCube& cube(void) const;
80 88  
81 89 protected:
82 90 // Protected methods
83   - void init_members(void);
84   - void copy_members(const ctbin& app);
85   - void free_members(void);
86   - void get_parameters(void);
87   - void fill_cube(GCTAObservation* obs);
88   - void set_weights(GCTAObservation* obs);
89   - void obs_cube(void);
  91 + void init_members(void);
  92 + void copy_members(const ctbin& app);
  93 + void free_members(void);
  94 + void get_parameters(void);
  95 + void init_sky_dir_cache(void);
  96 + GSkyMap create_cube(const GCTAObservation* obs);
  97 + GSkyMap fill_cube(const GCTAObservation* obs);
  98 + GSkyMap set_weights(const GCTAObservation* obs);
  99 + void obs_cube_stacked(void);
90 100  
91 101 // User parameters
92   - GFilename m_outcube; //!< Output counts map file name
93   - bool m_usepnt; //!< Use pointing instead of xref/yref parameters
94   - bool m_publish; //!< Publish counts cube?
95   - GChatter m_chatter; //!< Chattiness
  102 + bool m_usepnt; //!< Use pointing instead of xref/yref parameters
  103 + bool m_stack; //!< Output one stacked cube or multiple cubes
  104 + std::string m_prefix; //!< Prefix for multiple cubes
  105 + bool m_publish; //!< Publish counts cube?
  106 + GChatter m_chatter; //!< Chattiness
96 107  
97 108 // Protected members
98   - GSkyMap m_counts; //!< Event cube counts
99   - GSkyMap m_weights; //!< Event cube weights
100   - GEbounds m_ebounds; //!< Energy boundaries
101   - GGti m_gti; //!< Good time intervals
102   - GCTAEventCube m_cube; //!< Events cube (for cube() method)
103   - double m_ontime; //!< Total ontime
104   - double m_livetime; //!< Total livetime
  109 + std::vector<GCTAEventCube> m_cubes; //!< Event cubes
  110 + std::vector<GSkyMap> m_counts; //!< List of event cube counts
  111 + std::vector<GSkyMap> m_weights; //!< List of event cube weights
  112 + GEbounds m_ebounds; //!< Energy boundaries
  113 + GGti m_gti; //!< Stacked Good time intervals
  114 + double m_ontime; //!< Total ontime
  115 + double m_livetime; //!< Total livetime
105 116  
106 117 // Cache members
107   - std::vector<GSkyDir> m_dirs; //!< Cached GSkyDir for each pixel in m_counts
  118 + std::vector<GSkyDir> m_dirs; //!< Cached GSkyDir for all pixels
108 119 };
109 120  
110 121  
111 122 /***********************************************************************//**
112   - * @brief Return event cube
  123 + * @brief Return number of event cubes
113 124 *
114   - * @return Reference to event cube
  125 + * @return Number of event cubes.
115 126 *
116   - * Returns a reference to the event cube.
  127 + * Returns number of event cubes.
117 128 ***************************************************************************/
118 129 inline
119   -const GCTAEventCube& ctbin::cube(void) const
  130 +int ctbin::cubes(void) const
120 131 {
121   - return m_cube;
  132 + return m_cubes.size();
122 133 }
123 134  
124 135 #endif /* CTBIN_HPP */
... ...
src/ctbin/ctbin.par
  View file @433c9b7
... ... @@ -16,8 +16,10 @@
16 16 #
17 17 # General parameters
18 18 #===================
19   -inobs, f, a, events.fits,,, "Input event list or observation definition XML file"
20   -outcube, f, a, cntcube.fits,,, "Output counts cube file"
  19 +inobs, f, a, events.fits,,, "Input event list or observation definition XML file"
  20 +outobs, f, a, cntcube.fits,,, "Output counts cube file or observation definition XML file"
  21 +stack, b, h, yes,,, "Stack events into one cube?"
  22 +prefix, s, h, cntcube_,,, "Prefix for counts cubes in observation definition XML file"
21 23  
22 24 #
23 25 # Energy binning parameters
... ...
test/data/obs_unbinned_two.xml 0 → 100644
  View file @433c9b7
  1 +<?xml version="1.0" encoding="UTF-8" standalone="no"?>
  2 +<observation_list title="observation list">
  3 + <observation name="" id="00001" instrument="CTA">
  4 + <parameter name="EventList" file="$(TEST_DATA)/crab_events.fits" />
  5 + <parameter name="EffectiveArea" file="$(TEST_DATA)/irf_file.fits" />
  6 + <parameter name="PointSpreadFunction" file="$(TEST_DATA)/irf_file.fits" />
  7 + <parameter name="EnergyDispersion" file="$(TEST_DATA)/irf_file.fits" />
  8 + <parameter name="Background" file="$(TEST_DATA)/irf_file.fits" />
  9 + </observation>
  10 + <observation name="" id="00002" instrument="CTA">
  11 + <parameter name="EventList" file="$(TEST_DATA)/crab_events.fits" />
  12 + <parameter name="EffectiveArea" file="$(TEST_DATA)/irf_file.fits" />
  13 + <parameter name="PointSpreadFunction" file="$(TEST_DATA)/irf_file.fits" />
  14 + <parameter name="EnergyDispersion" file="$(TEST_DATA)/irf_file.fits" />
  15 + <parameter name="Background" file="$(TEST_DATA)/irf_file.fits" />
  16 + </observation>
  17 +</observation_list>
... ...
test/test_ctbin.py
  View file @433c9b7
... ... @@ -43,6 +43,9 @@ class Test(test):
43 43 # Call base class constructor
44 44 test.__init__(self)
45 45  
  46 + # Set members
  47 + self._inobs_two = self._datadir + '/obs_unbinned_two.xml'
  48 +
46 49 # Return
47 50 return
48 51  
... ... @@ -71,7 +74,7 @@ class Test(test):
71 74  
72 75 # Setup ctbin command
73 76 cmd = ctbin+' inobs="'+self._events+'"'+ \
74   - ' outcube="cntmap_cmd1.fits"'+\
  77 + ' outobs="cntmap_cmd1.fits"'+\
75 78 ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \
76 79 ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \
77 80 ' xref=83.63 yref=22.01 proj="CAR"'+ \
... ... @@ -87,7 +90,7 @@ class Test(test):
87 90  
88 91 # Setup ctbin command
89 92 cmd = ctbin+' inobs="events_that_do_not_exist.fits"'+ \
90   - ' outcube="cntmap_cmd2.fits"'+\
  93 + ' outobs="cntmap_cmd2.fits"'+\
91 94 ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \
92 95 ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \
93 96 ' xref=83.63 yref=22.01 proj="CAR"'+ \
... ... @@ -99,7 +102,7 @@ class Test(test):
99 102  
100 103 # Setup ctbin command with different ebounds than input
101 104 cmd = ctbin+' inobs="'+self._events+'"'+ \
102   - ' outcube="cntmap_cmd3.fits"'+\
  105 + ' outobs="cntmap_cmd3.fits"'+\
103 106 ' emin=2.0 emax=20.0 enumbins=10 ebinalg="LOG"'+ \
104 107 ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \
105 108 ' xref=83.63 yref=22.01 proj="CAR"'+ \
... ... @@ -113,6 +116,40 @@ class Test(test):
113 116 evt = gammalib.GCTAEventCube('cntmap_cmd3.fits')
114 117 self._check_cube(evt, 115)
115 118  
  119 + # Check joint (unstacked) binning
  120 +
  121 + # Setup ctbin command
  122 + cmd = ctbin+' inobs="'+self._inobs_two+'"'+ \
  123 + ' outobs="obs_unbinned_two_binned.xml"'+\
  124 + ' emin=1.0 emax=100.0 enumbins=10 ebinalg="LOG"'+ \
  125 + ' nxpix=40 nypix=40 binsz=0.1 coordsys="CEL"'+ \
  126 + ' xref=83.63 yref=22.01 proj="CAR"'+ \
  127 + ' logfile="ctbin_cmd4.log" chatter=1 stack=no'+ \
  128 + ' prefix="cntmap_cmd4_"'
  129 +
  130 + # Check if execution was successful
  131 + self.test_value(self._execute(cmd), 0,
  132 + 'Check successful execution from command line with stack=no')
  133 +
  134 + # Load counts cubes and check content.
  135 + evt = gammalib.GCTAEventCube('cntmap_cmd4_00001.fits')
  136 + self._check_cube(evt, 245)
  137 +
  138 + # Load counts cubes and check content.
  139 + evt = gammalib.GCTAEventCube('cntmap_cmd4_00002.fits')
  140 + self._check_cube(evt, 245)
  141 +
  142 + # Check observation definition XML output file
  143 + obs = gammalib.GObservations('obs_unbinned_two_binned.xml')
  144 + self.test_value(obs.size(), 2, 'Check for 2 observations in XML file')
  145 +
  146 + # Check for content of observations file
  147 + self.test_value(obs[0].eventfile().file(), 'cntmap_cmd4_cta_00001.fits',
  148 + 'Check first counts cube file name')
  149 + self.test_value(obs[1].eventfile().file(), 'cntmap_cmd4_cta_00002.fits',
  150 + 'Check second counts cube file name')
  151 +
  152 +
116 153 # Check ctbin --help
117 154 self._check_help(ctbin)
118 155  
... ... @@ -131,12 +168,12 @@ class Test(test):
131 168 # cube
132 169 self.test_value(binning.obs().size(), 0,
133 170 'Check that empty ctbin has an empty observation container')
134   - self.test_value(binning.cube().size(), 0,
  171 + self.test_value(binning.cubes(), 0,
135 172 'Check that empty ctbin has an empty counts cube')
136 173  
137 174 # Check that saving does not nothing
138 175 binning['logfile'] = 'ctbin_py0.log'
139   - binning['outcube'] = 'ctbin_py0.fits'
  176 + binning['outobs'] = 'ctbin_py0.fits'
140 177 binning.logFileOpen()
141 178 binning.save()
142 179 self.test_assert(not os.path.isfile('ctbin_py0.fits'),
... ... @@ -150,7 +187,7 @@ class Test(test):
150 187  
151 188 # Now set ctbin parameters
152 189 binning['inobs'] = self._events
153   - binning['outcube'] = 'ctbin_py1.fits'
  190 + binning['outobs'] = 'ctbin_py1.fits'
154 191 binning['ebinalg'] = 'LOG'
155 192 binning['emin'] = 1.0
156 193 binning['emax'] = 100.0
... ... @@ -186,11 +223,13 @@ class Test(test):
186 223 cpy_bin.logFileOpen()
187 224 cpy_bin.run()
188 225  
189   - # Check content of observation container and counts cube. Now an empty
190   - # event cube is expected since on input the observation is binned, and
191   - # any binned observation will be skipped.
192   - self._check_observation(cpy_bin, 0)
193   - self._check_cube(cpy_bin.cube(), 0)
  226 + # Check content of observation container and number of counts cubes.
  227 + # There should be a single binned observation in the observation
  228 + # container, which is the one that was produced in the run before.
  229 + # Since ctbin finds no unbinned observation in the container, the
  230 + # number of cubes should be zero.
  231 + self._check_observation(cpy_bin, 245)
  232 + self.test_value(cpy_bin.cubes(), 0, 'Check that there are no counts cubes')
194 233  
195 234 # Save counts cube
196 235 binning.save()
... ... @@ -206,7 +245,7 @@ class Test(test):
206 245 # counts cube
207 246 self.test_value(binning.obs().size(), 0,
208 247 'Check that empty ctbin has an empty observation container')
209   - self.test_value(binning.cube().size(), 0,
  248 + self.test_value(binning.cubes(), 0,
210 249 'Check that empty ctbin has an empty counts cube')
211 250  
212 251 # Prepare observation container for stacking of events into a
... ... @@ -215,7 +254,7 @@ class Test(test):
215 254  
216 255 # Set-up ctbin using an observation container
217 256 binning = ctools.ctbin(obs)
218   - binning['outcube'] = 'ctbin_py3.fits'
  257 + binning['outobs'] = 'ctbin_py3.fits'
219 258 binning['ebinalg'] = 'LOG'
220 259 binning['emin'] = 1.0
221 260 binning['emax'] = 100.0
... ... @@ -249,7 +288,7 @@ class Test(test):
249 288 pars = {'inobs': self._events, 'ebinalg': 'LIN', 'emin': 1.0,
250 289 'emax': 100.0, 'enumbins': 10, 'nxpix': 40, 'nypix': 40,
251 290 'binsz': 0.1, 'coordsys': 'CEL', 'proj': 'CAR',
252   - 'xref': 83.63, 'yref': 22.01, 'outcube': 'ctbin_py4.fits',
  291 + 'xref': 83.63, 'yref': 22.01, 'outobs': 'ctbin_py4.fits',
253 292 'logfile': 'ctbin_py4.log', 'chatter': 2}
254 293 binning = ctools.ctbin()
255 294 binning.pardict(pars)
... ... @@ -260,6 +299,42 @@ class Test(test):
260 299 evt = gammalib.GCTAEventCube('ctbin_py4.fits')
261 300 self._check_cube(evt, 245)
262 301  
  302 + # Test unstacked version
  303 + pars = {'inobs': self._datadir+'/obs_unbinned_two.xml', 'ebinalg': 'LIN', 'emin': 1.0,
  304 + 'emax': 100.0, 'enumbins': 10, 'nxpix': 40, 'nypix': 40,
  305 + 'binsz': 0.1, 'coordsys': 'CEL', 'proj': 'CAR',
  306 + 'xref': 83.63, 'yref': 22.01, 'outobs': 'ctbin_py5.xml',
  307 + 'logfile': 'ctbin_py5.log', 'chatter': 2,
  308 + 'prefix': 'cntcube_py5_', 'stack' : False}
  309 + binning = ctools.ctbin()
  310 + binning.pardict(pars)
  311 +
  312 + binning.logFileOpen()
  313 + binning.run()
  314 +
  315 + # Check individual cubes
  316 + self._check_cube(binning.cube(0), 245)
  317 + self._check_cube(binning.cube(1), 245)
  318 +
  319 + # Save counts cube
  320 + binning.save()
  321 +
  322 + # Load observations
  323 + obs = gammalib.GObservations('ctbin_py5.xml')
  324 + self.test_value(obs.size(), 2, 'Check number of output observations')
  325 +
  326 + # Check counts cubes
  327 + evt = obs[0].events()
  328 + self._check_cube(evt, 245)
  329 + evt = obs[1].events()
  330 + self._check_cube(evt, 245)
  331 +
  332 + # Load counts cubes and check content
  333 + evt = gammalib.GCTAEventCube('cntcube_py5_00001.fits')
  334 + self._check_cube(evt, 245)
  335 + evt = gammalib.GCTAEventCube('cntcube_py5_00002.fits')
  336 + self._check_cube(evt, 245)
  337 +
263 338 # Return
264 339 return
265 340  
... ...