Michael Mayer / ctools

Commit 3992ea547e74edc4bc4b4ebe6a3de1cc31758fbf

Authored by Jürgen Knödlseder
1 parent fe29729a

Add ctbkgcube tool. 

The ctbkgcube tool generates a background model for cube-style analysis.
Inline Side-by-side
Showing 22 changed files with 500 additions and 234 deletions
.gitignore
  View file @3992ea5
... ... @@ -54,6 +54,7 @@ src/ctselect/ctselect
54 54 src/ctskymap/ctskymap
55 55 src/ctexpcube/ctexpcube
56 56 src/ctpsfcube/ctpsfcube
  57 +src/ctbkgcube/ctbkgcube
57 58 src/ctcubemask/ctcubemask
58 59 test/pfiles/
59 60 test/*.fits
... ...
ChangeLog
  View file @3992ea5
1   -2014-07-20 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2014-07-22 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.8.0 released
4 4 ========================
... ... @@ -6,7 +6,8 @@
6 6 Adapt to GammaLib 0.9.0 interface
7 7 ctbin now always merges event lists into a single counts cube
8 8 Add ctexpcube tool (generates exposure cube)
9   - Add ctpsfcube tool (generated PSF cube)
  9 + Add ctpsfcube tool (generates PSF cube)
  10 + Add ctbkgcube tool (generates background cube)
10 11 Add ctcubemask tool (masks event bins in event cube)
11 12 Support energy dispersion in ctobssim, ctlike and ctmodel
12 13 Correct ctobssim for thread safety
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -20 July 2014
  5 +22 July 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -97,7 +97,15 @@ definition file. The cube definition may be taken from a counts cube
97 97 or may be specified by user parameters.
98 98  
99 99  
100   -12. ctcubemask - CTA cube filtering
  100 +12. ctbkgcube - CTA background cube generation
  101 +----------------------------------------------
  102 +The first version of the ctbkgcube tool is included in the release.
  103 +ctbkgcube computes the background cube for an observation definition
  104 +file. The cube definition may be taken from a counts cube or may be
  105 +specified by user parameters.
  106 +
  107 +
  108 +13. ctcubemask - CTA cube filtering
101 109 -----------------------------------
102 110 The first version of the ctcubemask tool is included in the release.
103 111 ctcubemask allows filtering of bins in a counts cube by setting the
... ... @@ -105,12 +113,12 @@ bin content to -1. Filtering can be done by energy, for a circular
105 113 Region of Interest (ROI) or using ds9 exclusion regions.
106 114  
107 115  
108   -13. Examples
  116 +14. Examples
109 117 ------------
110 118 None
111 119  
112 120  
113   -14. Scripts
  121 +15. Scripts
114 122 -----------
115 123 Add logL and event statistics to cspull output file.
116 124  
... ... @@ -119,11 +127,11 @@ Convert all tabs into spaces.
119 127 The obsutils script has now been moved into the ctools module namespace.
120 128  
121 129  
122   -15. caldb
  130 +16. caldb
123 131 ---------
124 132 None
125 133  
126 134  
127   -16. models
  135 +17. models
128 136 ----------
129 137 None
... ...
pyext/ctbkgcube.i 0 → 100644
  View file @3992ea5
  1 +/***************************************************************************
  2 + * ctbkgcube - CTA background cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctbkgcube.i
  23 + * @brief CTA background cube tool definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "ctbkgcube.hpp"
  30 +#include "GTools.hpp"
  31 +%}
  32 +
  33 +/***********************************************************************//**
  34 + * @class ctbkgcube
  35 + *
  36 + * @brief CTA background cube tool
  37 + ***************************************************************************/
  38 +class ctbkgcube : public GApplication {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + ctbkgcube(void);
  43 + explicit ctbkgcube(const GObservations& obs);
  44 + ctbkgcube(int argc, char *argv[]);
  45 + ctbkgcube(const ctbkgcube& app);
  46 + virtual ~ctbkgcube(void);
  47 +
  48 + // Methods
  49 + void clear(void);
  50 + void execute(void);
  51 + void run(void);
  52 + void save(void);
  53 + const GObservations& obs(void) const;
  54 + void get_parameters(void);
  55 + void get_obs(void);
  56 + void get_ebounds(void);
  57 + void set_from_cntmap(const std::string& filename);
  58 + void fill_cube(GCTAObservation* obs);
  59 +};
  60 +
  61 +/***********************************************************************//**
  62 + * @brief CTA background cube tool Python extension
  63 + ***************************************************************************/
  64 +%extend ctbkgcube {
  65 + ctbkgcube copy() {
  66 + return (*self);
  67 + }
  68 +}
... ...
pyext/ctools/ctools.i
  View file @3992ea5
... ... @@ -53,4 +53,5 @@
53 53 %include "ctskymap.i"
54 54 %include "ctexpcube.i"
55 55 %include "ctpsfcube.i"
  56 +%include "ctbkgcube.i"
56 57 %include "ctcubemask.i"
... ...
src/ctbkgcube/ctbkgcube.cpp
  View file @3992ea5
... ... @@ -31,11 +31,11 @@
31 31 #include <cstdio>
32 32 #include "ctbkgcube.hpp"
33 33 #include "GTools.hpp"
34   -#include "GWcs.hpp"
35 34  
36 35 /* __ Method name definitions ____________________________________________ */
37   -#define G_GET_EBOUNDS "ctbkgcube::get_ebounds()"
38 36 #define G_GET_PARAMETERS "ctbkgcube::get_parameters()"
  37 +#define G_GET_OBS "ctbkgcube::get_obs()"
  38 +#define G_GET_EBOUNDS "ctbkgcube::get_ebounds()"
39 39 #define G_SET_FROM_CNTMAP "ctbkgcube::set_from_cntmap(std::string&)"
40 40  
41 41 /* __ Debug definitions __________________________________________________ */
... ... @@ -230,7 +230,7 @@ void ctbkgcube::execute(void)
230 230 *
231 231 * This method reads the task parameters from the parfile, sets up the
232 232 * observation container, loops over all CTA observations in the container
233   - * and generates an background cube from the CTA observations.
  233 + * and generates a background cube from the CTA observations.
234 234 ***************************************************************************/
235 235 void ctbkgcube::run(void)
236 236 {
... ... @@ -248,18 +248,6 @@ void ctbkgcube::run(void)
248 248 log << std::endl;
249 249 }
250 250  
251   - // Set energy dispersion flag for all CTA observations and save old
252   - // values in save_edisp vector
253   - std::vector<bool> save_edisp;
254   - save_edisp.assign(m_obs.size(), false);
255   - for (int i = 0; i < m_obs.size(); ++i) {
256   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
257   - if (obs != NULL) {
258   - save_edisp[i] = obs->response()->apply_edisp();
259   - obs->response()->apply_edisp(m_apply_edisp);
260   - }
261   - }
262   -
263 251 // Write observation(s) into logger
264 252 if (logTerse()) {
265 253 log << std::endl;
... ... @@ -275,83 +263,149 @@ void ctbkgcube::run(void)
275 263 // Write header
276 264 if (logTerse()) {
277 265 log << std::endl;
278   - log.header1("Generate background cube");
  266 + log.header1("Prepare model");
279 267 }
280 268  
281   - // If there are no models associated with the observations then
282   - // load now the model definition from the XML file
283   - if (m_obs.models().size() == 0) {
284   - m_bgmdl = (*this)["bgmdl"].filename();
285   - m_obs.models(GModels(m_bgmdl));
286   -
287   - // TODO: throw a message when no background definition file is given
  269 + // Copy models from observation container
  270 + m_bkgmdl = m_obs.models();
  271 +
  272 + // Remove all model that are not CTA background models from the
  273 + // container
  274 + int num = m_bkgmdl.size();
  275 + for (int i = num-1; i >= 0; --i) {
  276 +
  277 + // Flag removal
  278 + bool remove = true;
  279 +
  280 + // Do we have a CTA, HESS or VERITAS specific model?
  281 + if (m_bkgmdl[i]->is_valid("CTA", "") ||
  282 + m_bkgmdl[i]->is_valid("HESS", "") ||
  283 + m_bkgmdl[i]->is_valid("VERITAS", "")) {
  284 +
  285 + // Do we have a background model?
  286 + if (dynamic_cast<GModelData*>(m_bkgmdl[i]) != NULL) {
  287 +
  288 + // ... then keep model
  289 + remove = false;
  290 +
  291 + } // endif: had a background model
  292 + } // endif: had a CTA, HESS or VERITAS model
  293 +
  294 + // Log results
  295 + if (logTerse()) {
  296 + if (remove) {
  297 + log << gammalib::parformat("Remove model");
  298 + }
  299 + else {
  300 + log << gammalib::parformat("Keep model");
  301 + }
  302 + log << m_bkgmdl[i]->name();
  303 + log << " ";
  304 + log << m_bkgmdl[i]->type();
  305 + log << "(";
  306 + log << m_bkgmdl[i]->instruments();
  307 + log << ")";
  308 + log << std::endl;
  309 + }
  310 +
  311 + // Remove model if requested
  312 + if (remove) {
  313 + m_bkgmdl.remove(i);
  314 + }
  315 +
  316 + } // endfor: looped over all background models
  317 +
  318 + // Write header
  319 + if (logTerse()) {
  320 + log << std::endl;
  321 + log.header1("Generate background cube");
288 322 }
289 323  
  324 + // Initialise observation counter
  325 + int n_observations = 0;
  326 +
290 327 // Loop over all observations in the container
291 328 for (int i = 0; i < m_obs.size(); ++i) {
292 329  
293 330 // Get CTA observation
294 331 GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
295   -
296   - // Fill background cube
297   - fill_cube(obs, m_obs.models());
298   -
299   - }
300 332  
301   - // Restore energy dispersion flag for all CTA observations
302   - for (int i = 0; i < m_obs.size(); ++i) {
303   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  333 + // Continue only if observation is a CTA observation
304 334 if (obs != NULL) {
305   - obs->response()->apply_edisp(save_edisp[i]);
306   - }
  335 +
  336 + // Write header for observation
  337 + if (logTerse()) {
  338 + if (obs->name().length() > 1) {
  339 + log.header3("Observation "+obs->name());
  340 + }
  341 + else {
  342 + log.header3("Observation "+gammalib::str(i));
  343 + }
  344 + }
  345 +
  346 + // Increment number of observations
  347 + n_observations++;
  348 +
  349 + // Fill the cube
  350 + fill_cube(obs);
  351 +
  352 + } // endif: CTA observation found
  353 +
  354 + } // endfor: looped over observations
  355 +
  356 + // Log results
  357 + if (logTerse()) {
  358 + log << gammalib::parformat("Number of observations");
  359 + log << n_observations << std::endl;
307 360 }
308 361  
309 362 // Return
310 363 return;
311 364 }
312 365  
  366 +
313 367 /***********************************************************************//**
314 368 * @brief Generate background cube
315 369 *
316   - * @param[in] obs CTA observation pointer.
317   - * @param[in] models Model container.
318   - *
  370 + * @param[in] obs Pointer to CTA observation.
319 371 *
  372 + * Fills
320 373 ***************************************************************************/
321   -void ctbkgcube::fill_cube(GCTAObservation* obs, const GModels& models)
  374 +void ctbkgcube::fill_cube(GCTAObservation* obs)
322 375 {
323 376 // Continue only if observation pointer is valid
324 377 if (obs != NULL) {
325 378  
326   - obs->events(m_bkgcube);
327   - GCTAEventCube* cube =
328   - const_cast<GCTAEventCube*>(dynamic_cast<const GCTAEventCube*>(obs->events()));
329 379 // Initialise statistics
330 380 double sum = 0.0;
331   - GCTAEventBin* bin = (*cube)[10000];
332   - std::cout << bin->counts() << std::endl;
333   - // Loop over all events in counts map
334   - for (int i = 0; i < cube->size(); ++i) {
335   -
  381 +
  382 + // Set GTI of actual observations as the GTI of the event cube
  383 + m_bkgcube.gti(obs->events()->gti());
  384 +
  385 + // Loop over all bins in background cube
  386 + for (int i = 0; i < m_bkgcube.size(); ++i) {
  387 +
336 388 // Get event bin
337   - GCTAEventBin* bin = (*cube)[i];
  389 + GCTAEventBin* bin = m_bkgcube[i];
338 390  
339 391 // Compute model value for event bin
340 392 double model =
341   - models.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
  393 + m_bkgmdl.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
342 394 bin->size();
343   -
344   - // Store value
345   - bin->counts( bin->counts() + model );
346   -
347   - // Sum all events
348   - sum += model;
349 395  
  396 + // Add existing number of counts
  397 + model += bin->counts();
  398 +
  399 + // Store cumulated value
  400 + bin->counts(model);
  401 +
  402 + // Sum all events
  403 + sum += model;
350 404 }
351   -
  405 +
352 406 // Log results
353 407 if (logTerse()) {
354   - log << gammalib::parformat("Model events in cube");
  408 + log << gammalib::parformat("Background events in cube");
355 409 log << sum << std::endl;
356 410 }
357 411  
... ... @@ -361,9 +415,11 @@ void ctbkgcube::fill_cube(GCTAObservation* obs, const GModels&amp; models)
361 415 return;
362 416 }
363 417  
  418 +
364 419 /***********************************************************************//**
365 420 * @brief Save background cube
366 421 *
  422 + * Save the background cube into the file specified by the outfile parameter.
367 423 ***************************************************************************/
368 424 void ctbkgcube::save(void)
369 425 {
... ... @@ -373,11 +429,26 @@ void ctbkgcube::save(void)
373 429 log.header1("Save background cube");
374 430 }
375 431  
  432 + // Create energies container from energy boundaries
  433 + GEnergies energies;
  434 + for (int i = 0; i < m_ebounds.size(); ++i) {
  435 + energies.append(m_ebounds.elogmean(i));
  436 + }
  437 +
376 438 // Get output filename
377   - m_outfile = (*this)["outfile"].filename();
  439 + std::string filename = (*this)["outfile"].filename();
  440 +
  441 + // Create empty FITS file
  442 + GFits fits;
378 443  
379   - // Save background cube
380   - m_map.save(m_outfile, clobber());
  444 + // Write background cube
  445 + m_bkgcube.map().write(fits);
  446 +
  447 + // Write energies
  448 + energies.write(fits);
  449 +
  450 + // Save FITS file
  451 + fits.saveto(filename, clobber());
381 452  
382 453 // Return
383 454 return;
... ... @@ -387,32 +458,43 @@ void ctbkgcube::save(void)
387 458 /***********************************************************************//**
388 459 * @brief Get application parameters
389 460 *
  461 + * @exception GException::invalid_value
  462 + * No background model definition XML file specified.
  463 + *
390 464 * Get all task parameters from parameter file or (if required) by querying
391 465 * the user. The parameters are read in the correct order.
392   - *
393   - * @todo Setup background cube from counts map
394 466 ***************************************************************************/
395 467 void ctbkgcube::get_parameters(void)
396 468 {
397   - // If we do not have any observations in the container then get an
398   - // input file name or observation descriptor file
  469 + // If we do not have any observations in the container then load them
  470 + // using the "infile" parameter
399 471 if (m_obs.size() == 0) {
400 472 get_obs();
401 473 }
402 474  
403   - // Make sure that response is set
404   - set_response();
405   -
406   - // Read energy dispersion flag
407   - m_apply_edisp = (*this)["edisp"].boolean();
  475 + // If there are no models associated with the observations then load now
  476 + // the model definition from the XML file
  477 + if (m_obs.models().size() == 0) {
  478 + if ((*this)["bkgmdl"].is_undefined()) {
  479 + std::string msg = "No background model definition XML file "
  480 + "specified. Please set the \"bkgmdl\" parameter "
  481 + "to the background model definition XML "
  482 + "filename.";
  483 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  484 + }
  485 + std::string bkgmdl = (*this)["bkgmdl"].filename();
  486 + GModels models(bkgmdl);
  487 + m_obs.models(models);
  488 + }
408 489  
409 490 // If no counts map is specified then setup the background cube from
410 491 // the user parameters
411 492 std::string cntmap = (*this)["cntmap"].filename();
412   - if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) {
  493 + if ((gammalib::toupper(cntmap) == "NONE") ||
  494 + (gammalib::strip_whitespace(cntmap) == "")) {
413 495  
414 496 // Get user parameters for counts map definition
415   - std::string wcs = (*this)["proj"].string();
  497 + std::string proj = (*this)["proj"].string();
416 498 std::string coordsys = (*this)["coordsys"].string();
417 499 double xref = (*this)["xref"].real();
418 500 double yref = (*this)["yref"].real();
... ... @@ -423,21 +505,16 @@ void ctbkgcube::get_parameters(void)
423 505 // Get energy definition
424 506 get_ebounds();
425 507  
426   - // Define background cube
427   - m_map = GSkymap(wcs, coordsys, xref, yref,
428   - -binsz, binsz, nxpix, nypix,
429   - m_ebounds.size());
430   - GGti gti = m_obs[0]->events()->gti();
431   - m_bkgcube = GCTAEventCube(m_map, m_ebounds, gti);
432   -
433   - // Loop over all events in background cube and
434   - // set all pixels to 0
435   - for (int i = 0; i < m_bkgcube.size(); ++i)
436   - {
437   - // Get event bin
438   - GCTAEventBin* bin = m_bkgcube[i];
439   - bin->counts(0);
440   - }
  508 + // Define skymap for background cube
  509 + GSkymap map(proj, coordsys, xref, yref,
  510 + -binsz, binsz, nxpix, nypix,
  511 + m_ebounds.size());
  512 +
  513 + // Allocate background cube using a dummy GTI
  514 + GGti gti;
  515 + gti.append(GTime(0.0), GTime(1.0)); // Dummy GTI
  516 + m_bkgcube = GCTAEventCube(map, m_ebounds, gti);
  517 +
441 518 }
442 519  
443 520 // ... otherwise setup the background cube from the counts map
... ... @@ -450,7 +527,7 @@ void ctbkgcube::get_parameters(void)
450 527  
451 528 // Read output filename (if needed)
452 529 if (m_read_ahead) {
453   - m_outfile = (*this)["outfile"].filename();
  530 + std::string filename = (*this)["outfile"].filename();
454 531 }
455 532  
456 533 // Return
... ... @@ -461,6 +538,10 @@ void ctbkgcube::get_parameters(void)
461 538 /***********************************************************************//**
462 539 * @brief Get observation definition
463 540 *
  541 + * @exception GException::invalid_value
  542 + * No event list, counts cube or observation definition file
  543 + * specified.
  544 + *
464 545 * Get observation definition from the user parameters.
465 546 ***************************************************************************/
466 547 void ctbkgcube::get_obs(void)
... ... @@ -468,6 +549,15 @@ void ctbkgcube::get_obs(void)
468 549 // Get input filename
469 550 std::string filename = (*this)["infile"].filename();
470 551  
  552 + // Check whether infile parameter is undefined
  553 + if (gammalib::toupper(filename) == "NONE") {
  554 + std::string msg = "No event list, counts cube or observation definition "
  555 + "file specified. Please set the \"infile\" parameter "
  556 + "to either an event list filename, a counts cube "
  557 + "filename or an observation definition filename.";
  558 + throw GException::invalid_value(G_GET_OBS, msg);
  559 + }
  560 +
471 561 // Try first to open as FITS file
472 562 try {
473 563  
... ... @@ -479,7 +569,6 @@ void ctbkgcube::get_obs(void)
479 569  
480 570 // Append CTA observation to container
481 571 m_obs.append(obs);
482   -
483 572  
484 573 }
485 574  
... ... @@ -497,59 +586,11 @@ void ctbkgcube::get_obs(void)
497 586  
498 587  
499 588 /***********************************************************************//**
500   - * @brief Set observation response
501   - *
502   - * Set response for all observations that so have no response.
503   - ***************************************************************************/
504   -void ctbkgcube::set_response(void)
505   -{
506   - // Loop over all observations
507   - for (int i = 0; i < m_obs.size(); ++i) {
508   -
509   - // Is this observation a CTA observation?
510   - GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
511   -
512   - // Yes ...
513   - if (obs != NULL) {
514   -
515   - // Set response if we don't have one
516   - if (!obs->hasresponse()) {
517   -
518   - // Load response information
519   - std::string database = (*this)["caldb"].string();
520   - std::string irf = (*this)["irf"].string();
521   -
522   - // Set calibration database. If specified parameter is a
523   - // directory then use this as the pathname to the calibration
524   - // database. Otherwise interpret this as the instrument name,
525   - // the mission being "cta"
526   - GCaldb caldb;
527   - if (gammalib::dir_exists(database)) {
528   - caldb.rootdir(database);
529   - }
530   - else {
531   - caldb.open("cta", database);
532   - }
533   -
534   - // Set reponse
535   - obs->response(irf, caldb);
536   -
537   - } // endif: observation already has a response
538   -
539   - } // endif: observation was a CTA observation
540   -
541   - } // endfor: looped over all observations
542   -
543   - // Return
544   - return;
545   -}
546   -
547   -
548   -/***********************************************************************//**
549 589 * @brief Get energy boundaries from parameters
550 590 *
551 591 * @exception GException::invalid_value
552   - * Invalid extension name encountered.
  592 + * Energy bin definition file not defined or invalid extension
  593 + * name encountered.
553 594 *
554 595 * Get the energy boundaries from the user parameters.
555 596 ***************************************************************************/
... ... @@ -561,10 +602,18 @@ void ctbkgcube::get_ebounds(void)
561 602 // If we have the binning given by a file then try to get the boundaries
562 603 // from that file
563 604 if (ebinalg == "FILE") {
564   -
  605 +
565 606 // Get filename
566 607 std::string filename = (*this)["ebinfile"].filename();
567 608  
  609 + // Check whether infile parameter is undefined
  610 + if (gammalib::toupper(filename) == "NONE") {
  611 + std::string msg = "No energy bin definition file specified. "
  612 + "Please set the \"ebinfile\" parameter to a "
  613 + "file containing the energy bin definition.";
  614 + throw GException::invalid_value(G_GET_EBOUNDS, msg);
  615 + }
  616 +
568 617 // Open fits file to check which extension is given
569 618 GFits file(filename);
570 619  
... ... @@ -655,22 +704,14 @@ void ctbkgcube::set_from_cntmap(const std::string&amp; filename)
655 704  
656 705 // Get energy definition
657 706 m_ebounds = cube->ebounds();
658   -
659   -
660   - // Define background cube
661   - m_map = GSkymap(proj, coordsys, xref, yref,
662   - dx, dy, nx, ny,
663   - m_ebounds.size());
664   - GGti gti = m_obs[0]->events()->gti();
665   - m_bkgcube = GCTAEventCube(m_map, m_ebounds, gti);
666   - // Loop over all events in background cube and
667   - // set all pixels to 0
668   - for (int i = 0; i < m_bkgcube.size(); ++i)
669   - {
670   - // Get event bin
671   - GCTAEventBin* bin = m_bkgcube[i];
672   - bin->counts(0);
673   - }
  707 +
  708 + // Define skymap for background cube
  709 + GSkymap map(proj, coordsys, xref, yref,
  710 + dx, dy, nx, ny,
  711 + m_ebounds.size());
  712 +
  713 + // Allocate background cube
  714 + m_bkgcube = GCTAEventCube(map, m_ebounds, obs.events()->gti());
674 715  
675 716 } // endif: WCS projection was valid
676 717  
... ... @@ -706,14 +747,10 @@ void ctbkgcube::set_from_cntmap(const std::string&amp; filename)
706 747 void ctbkgcube::init_members(void)
707 748 {
708 749 // Initialise members
709   - m_outfile.clear();
710   - m_apply_edisp = false;
711   -
712   - // Initialise protected members
713 750 m_read_ahead = false;
714 751 m_obs.clear();
715 752 m_bkgcube.clear();
716   - m_map.clear();
  753 + m_bkgmdl.clear();
717 754 m_ebounds.clear();
718 755  
719 756 // Set logger properties
... ... @@ -731,15 +768,11 @@ void ctbkgcube::init_members(void)
731 768 ***************************************************************************/
732 769 void ctbkgcube::copy_members(const ctbkgcube& app)
733 770 {
734   - // Copy attributes
735   - m_outfile = app.m_outfile;
736   - m_apply_edisp = app.m_apply_edisp;
737   -
738   - // Copy protected members
  771 + // Copy members
739 772 m_read_ahead = app.m_read_ahead;
740 773 m_obs = app.m_obs;
741 774 m_bkgcube = app.m_bkgcube;
742   - m_map = app.m_map;
  775 + m_bkgmdl = app.m_bkgmdl;
743 776 m_ebounds = app.m_ebounds;
744 777  
745 778 // Return
... ...
src/ctbkgcube/ctbkgcube.hpp
  View file @3992ea5
1 1 /***************************************************************************
2 2 * ctbkgcube - CTA background cube tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,36 +55,42 @@ public:
55 55 ctbkgcube& operator=(const ctbkgcube& app);
56 56  
57 57 // Methods
58   - void clear(void);
59   - void execute(void);
60   - void run(void);
61   - void save(void);
62   - GObservations& obs(void) { return m_obs; }
63   - void get_parameters(void);
64   - void get_obs(void);
65   - void set_response(void);
66   - void get_ebounds(void);
67   - void set_from_cntmap(const std::string& filename);
68   - void fill_cube(GCTAObservation* obs, const GModels& models);
  58 + void clear(void);
  59 + void execute(void);
  60 + void run(void);
  61 + void save(void);
  62 + const GObservations& obs(void) const;
  63 + void get_parameters(void);
  64 + void get_obs(void);
  65 + void get_ebounds(void);
  66 + void set_from_cntmap(const std::string& filename);
  67 + void fill_cube(GCTAObservation* obs);
69 68  
70 69 protected:
71 70 // Protected methods
72   - void init_members(void);
73   - void copy_members(const ctbkgcube& app);
74   - void free_members(void);
75   -
76   - // User parameters
77   - std::string m_outfile; //!< Output background cube file
78   - std::string m_bgmdl; //!< background model
79   - bool m_apply_edisp; //!< Apply energy dispersion?
  71 + void init_members(void);
  72 + void copy_members(const ctbkgcube& app);
  73 + void free_members(void);
80 74  
81 75 // Protected members
82 76 bool m_read_ahead; //!< Read ahead parameters
83 77 GObservations m_obs; //!< Observation container
84 78 GCTAEventCube m_bkgcube; //!< Background cube
85   - GSkymap m_map; //!< Background map
  79 + GModels m_bkgmdl; //!< CTA background models
86 80 GEbounds m_ebounds; //!< Energy boundaries
87   -
88 81 };
89 82  
  83 +
  84 +/***********************************************************************//**
  85 + * @brief Return observation container
  86 + *
  87 + * @return Observation container.
  88 + ***************************************************************************/
  89 +inline
  90 +const GObservations& ctbkgcube::obs(void) const
  91 +{
  92 + // Return observation container
  93 + return m_obs;
  94 +}
  95 +
90 96 #endif /* CTBKGCUBE_HPP */
... ...
src/ctbkgcube/ctbkgcube.par
  View file @3992ea5
... ... @@ -4,26 +4,23 @@
4 4 # #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7   -# File: ctexpcube.par #
  7 +# File: ctbkgcube.par #
8 8 # Version: 00-01-00 #
9 9 # #
10   -# Author: Juergen Knoedlseder #
11   -# Jurgen.Knodlseder@irap.omp.eu #
12   -# IRAP #
  10 +# Author: Chia-Chun Lu #
  11 +# c.lu@mpi-hd.mpg.de #
  12 +# MPIK #
13 13 # #
14   -# Purpose: Parameter file for the ctexpcube application #
  14 +# Purpose: Parameter file for the ctbkgcube application #
15 15 # #
16 16 #############################################################################
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Observation definition file"
21   -bgmdl, f, a, "NONE",,, "Background definition file"
22   -cntmap, f, a, "NONE",,, "Counts map for background cube definition"
  20 +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file"
  21 +bkgmdl, f, a, "NONE",,, "Background model XML file"
  22 +cntmap, f, a, "NONE",,, "Counts cube for background cube definition"
23 23 outfile, f, a, "bkgcube.fits",,, "Output background cube file"
24   -caldb, s, a, "dummy",,, "Calibration database"
25   -irf, s, a, "cta_dummy_irf",,, "Instrument response function"
26   -edisp, b, h, no,,, "Apply energy dispersion?"
27 24 #
28 25 # Energy binning parameters
29 26 #==========================
... ... @@ -50,4 +47,4 @@ chatter, i, h, 2,0,4, &quot;Chattiness of output&quot;
50 47 clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
51 48 debug, b, h, no,,, "Debugging mode activated"
52 49 mode, s, h, "ql",,, "Mode of automatic parameters"
53   -logfile, f, h, "ctexpcube.log",,, "Log filename"
  50 +logfile, f, h, "ctbkgcube.log",,, "Log filename"
... ...
src/ctbkgcube/main.cpp
  View file @3992ea5
1 1 /***************************************************************************
2   - * ctbkgcube - CTA exposure cube tool *
  2 + * ctbkgcube - CTA background cube tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2014 by Chia-Chun Lu *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ...
src/ctexpcube/ctexpcube.par
  View file @3992ea5
... ... @@ -17,7 +17,7 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Observation definition file"
  20 +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file"
21 21 cntmap, f, a, "NONE",,, "Counts map for exposure cube definition"
22 22 outfile, f, a, "expcube.fits",,, "Output exposure cube file"
23 23 caldb, s, a, "dummy",,, "Calibration database"
... ...
src/ctpsfcube/ctpsfcube.par
  View file @3992ea5
... ... @@ -17,7 +17,7 @@
17 17 #
18 18 # General parameters
19 19 #===================
20   -infile, f, a, "NONE",,, "Observation definition file"
  20 +infile, f, a, "NONE",,, "Event list, counts cube or observation definition file"
21 21 cntmap, f, a, "NONE",,, "Counts map for exposure cube definition"
22 22 outfile, f, a, "psfcube.fits",,, "Output exposure cube file"
23 23 caldb, s, a, "dummy",,, "Calibration database"
... ...
test/test_ctbin.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# ========================= #
  26 +# Test class for ctbin tool #
  27 +# ========================= #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctbin tool.
... ...
test/test_ctbkgcube.py 0 → 100755
  View file @3992ea5
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the ctbkgcube tool.
  4 +#
  5 +# Copyright (C) 2014 Juergen Knoedlseder
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import ctools
  23 +
  24 +
  25 +# ============================= #
  26 +# Test class for ctbkgcube tool #
  27 +# ============================= #
  28 +class Test(gammalib.GPythonTestSuite):
  29 + """
  30 + Test class for ctbkgcube tool.
  31 + """
  32 + # Constructor
  33 + def __init__(self):
  34 + """
  35 + Constructor.
  36 + """
  37 + # Call base class constructor
  38 + gammalib.GPythonTestSuite.__init__(self)
  39 +
  40 + # Set members
  41 + self.events_name = "data/crab_events.fits"
  42 + self.bkg_model = "data/crab.xml"
  43 +
  44 + # Return
  45 + return
  46 +
  47 + # Set test functions
  48 + def set(self):
  49 + """
  50 + Set all test functions.
  51 + """
  52 + # Set test name
  53 + self.name("ctbkgcube")
  54 +
  55 + # Append tests
  56 + self.append(self.test_functional, "Test ctbkgcube functionality")
  57 +
  58 + # Return
  59 + return
  60 +
  61 + # Test ctbkgcube functionnality
  62 + def test_functional(self):
  63 + """
  64 + Test ctbkgcube functionnality.
  65 + """
  66 + # Set-up ctbkgcube
  67 + bkgcube = ctools.ctbkgcube()
  68 + bkgcube["infile"].filename(self.events_name)
  69 + bkgcube["bkgmdl"].filename(self.bkg_model)
  70 + bkgcube["cntmap"].filename("NONE")
  71 + bkgcube["outfile"].filename("psfcube.fits")
  72 + bkgcube["ebinalg"].string("LOG")
  73 + bkgcube["emin"].real(0.1)
  74 + bkgcube["emax"].real(100.0)
  75 + bkgcube["enumbins"].integer(20)
  76 + bkgcube["nxpix"].integer(10)
  77 + bkgcube["nypix"].integer(10)
  78 + bkgcube["binsz"].real(0.4)
  79 + bkgcube["coordsys"].string("CEL")
  80 + bkgcube["proj"].string("CAR")
  81 + bkgcube["xref"].real(83.63)
  82 + bkgcube["yref"].real(22.01)
  83 +
  84 + # Run tool
  85 + self.test_try("Run ctbkgcube")
  86 + try:
  87 + bkgcube.run()
  88 + self.test_try_success()
  89 + except:
  90 + self.test_try_failure("Exception occured in ctbkgcube.")
  91 +
  92 + # Save background cube
  93 + self.test_try("Save background cube")
  94 + try:
  95 + bkgcube.save()
  96 + self.test_try_success()
  97 + except:
  98 + self.test_try_failure("Exception occured in saving background cube.")
  99 +
  100 + # Return
  101 + return
... ...
test/test_ctcubemask.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# ============================== #
  26 +# Test class for ctcubemask tool #
  27 +# ============================== #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctcubemask tool.
... ...
test/test_ctexpcube.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# ============================= #
  26 +# Test class for ctexpcube tool #
  27 +# ============================= #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctexpcube tool.
... ...
test/test_ctlike.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# ========================== #
  26 +# Test class for ctlike tool #
  27 +# ========================== #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctlike tool.
... ...
test/test_ctmodel.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# =========================== #
  26 +# Test class for ctmodel tool #
  27 +# =========================== #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctmodel tool.
... ...
test/test_ctools.sh
  View file @3992ea5
... ... @@ -42,6 +42,7 @@ ctselect=../src/ctselect/ctselect
42 42 ctskymap=../src/ctskymap/ctskymap
43 43 ctexpcube=../src/ctexpcube/ctexpcube
44 44 ctpsfcube=../src/ctpsfcube/ctpsfcube
  45 +ctbkgcube=../src/ctbkgcube/ctbkgcube
45 46 ctcubemask=../src/ctcubemask/ctcubemask
46 47  
47 48  
... ... @@ -444,6 +445,51 @@ $ECHO &quot; ok&quot;
444 445  
445 446  
446 447 #
  448 +# Test ctbkgcube
  449 +# ==============
  450 +$ECHO -n "Test ctbkgcube: "
  451 +#
  452 +# Run 1
  453 +$ctbkgcube infile="data/crab_events.fits" \
  454 + bkgmdl="data/crab.xml" \
  455 + cntmap="NONE" \
  456 + outfile="bkgcube1.fits" \
  457 + emin=0.1 \
  458 + emax=100.0 \
  459 + enumbins=20 \
  460 + nxpix=10 \
  461 + nypix=10 \
  462 + binsz=0.4 \
  463 + coordsys="CEL" \
  464 + xref=83.63 \
  465 + yref=22.01 \
  466 + proj="CAR"
  467 +$ECHO -n "."
  468 +if [ -s "bkgcube1.fits" ]
  469 +then
  470 + $ECHO -n "."
  471 +else
  472 + $ECHO " bkgcube1.fits file is not valid"
  473 + exit 1
  474 +fi
  475 +#
  476 +# Run 2
  477 +$ctbkgcube infile="data/crab_events.fits" \
  478 + bkgmdl="data/crab.xml" \
  479 + cntmap="data/crab_cntmap.fits" \
  480 + outfile="bkgcube2.fits"
  481 +$ECHO -n "."
  482 +if [ -s "bkgcube2.fits" ]
  483 +then
  484 + $ECHO -n "."
  485 +else
  486 + $ECHO " bkgcube2.fits file is not valid"
  487 + exit 1
  488 +fi
  489 +$ECHO " ok"
  490 +
  491 +
  492 +#
447 493 # Test ctcubemask
448 494 # ===============
449 495 $ECHO -n "Test ctcubemask: "
... ...
test/test_ctpsfcube.py
  View file @3992ea5
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ============================ #
26   -# Test class for ctobssim tool #
27   -# ============================ #
  25 +# ============================= #
  26 +# Test class for ctpsfcube tool #
  27 +# ============================= #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctpsfcube tool.
... ...
test/test_ctselect.py
  View file @3992ea5
... ... @@ -23,7 +23,7 @@ import ctools
23 23  
24 24  
25 25 # ============================ #
26   -# Test class for ctobssim tool #
  26 +# Test class for ctselect tool #
27 27 # ============================ #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
... ...
test/test_ctskymap.py
  View file @3992ea5
... ... @@ -23,7 +23,7 @@ import ctools
23 23  
24 24  
25 25 # ============================ #
26   -# Test class for ctobssim tool #
  26 +# Test class for ctskymap tool #
27 27 # ============================ #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
... ...
test/test_python.py
  View file @3992ea5
... ... @@ -28,9 +28,10 @@ import test_ctbin
28 28 import test_ctlike
29 29 import test_ctmodel
30 30 import test_ctskymap
31   -import test_ctcubemask
32 31 import test_ctexpcube
33 32 import test_ctpsfcube
  33 +import test_ctbkgcube
  34 +import test_ctcubemask
34 35 import test_pipelines
35 36  
36 37  
... ... @@ -51,9 +52,10 @@ if __name__ == &#39;__main__&#39;:
51 52 suite_ctlike = test_ctlike.Test()
52 53 suite_ctmodel = test_ctmodel.Test()
53 54 suite_ctskymap = test_ctskymap.Test()
54   - suite_ctcubemask = test_ctcubemask.Test()
55 55 suite_ctexpcube = test_ctexpcube.Test()
56 56 suite_ctpsfcube = test_ctpsfcube.Test()
  57 + suite_ctbkgcube = test_ctbkgcube.Test()
  58 + suite_ctcubemask = test_ctcubemask.Test()
57 59 suite_pipelines = test_pipelines.Test()
58 60  
59 61 # Setup unit tests
... ... @@ -63,9 +65,10 @@ if __name__ == &#39;__main__&#39;:
63 65 suite_ctlike.set()
64 66 suite_ctmodel.set()
65 67 suite_ctskymap.set()
66   - suite_ctcubemask.set()
67 68 suite_ctexpcube.set()
68 69 suite_ctpsfcube.set()
  70 + suite_ctbkgcube.set()
  71 + suite_ctcubemask.set()
69 72 suite_pipelines.set()
70 73  
71 74 # Append tests to container
... ... @@ -75,9 +78,10 @@ if __name__ == &#39;__main__&#39;:
75 78 suites.append(suite_ctlike)
76 79 suites.append(suite_ctmodel)
77 80 suites.append(suite_ctskymap)
78   - suites.append(suite_ctcubemask)
79 81 suites.append(suite_ctexpcube)
80 82 suites.append(suite_ctpsfcube)
  83 + suites.append(suite_ctbkgcube)
  84 + suites.append(suite_ctcubemask)
81 85 suites.append(suite_pipelines)
82 86  
83 87 # Set PFILES environment variable
... ...