Francois Brun / ctools

Commit c2675d4fe5c7b3fb3edf5e758dafbcaefad62234

Authored by Jürgen Knödlseder
1 parent b2608b4d

Rename some of the base class methods and reduce further some code redundancies … 

…(see #1360 for details).
Inline Side-by-side
Showing 13 changed files with 409 additions and 690 deletions
src/ctbin/ctbin.cpp
  View file @c2675d4
... ... @@ -34,7 +34,7 @@
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
36 36 #define G_FILL_CUBE "ctbin::fill_cube(GCTAObservation*)"
37   -#define G_GET_PARAMETERS "ctbin::get_parameters()"
  37 +#define G_GET_PARAMETERS "ctbin::get_parameters()"
38 38  
39 39 /* __ Debug definitions __________________________________________________ */
40 40  
... ... @@ -386,9 +386,6 @@ void ctbin::free_members(void)
386 386 /***********************************************************************//**
387 387 * @brief Get application parameters
388 388 *
389   - * @exception GException::invalid_value
390   - * Parameter "inobs" is required for ctbin.
391   - *
392 389 * Get all task parameters from parameter file or (if required) by querying
393 390 * the user. Most parameters are only required if no observation exists so
394 391 * far in the observation container. In this case, a single CTA observation
... ... @@ -397,47 +394,22 @@ void ctbin::free_members(void)
397 394 ***************************************************************************/
398 395 void ctbin::get_parameters(void)
399 396 {
400   - // If there are no observations in container then load them via user parameters
  397 + // If there are no observations in container then load them via user
  398 + // parameters
401 399 if (m_obs.size() == 0) {
402 400  
403   - // Throw exception if no infile is given
404   - if ((*this)["inobs"].filename() == "NONE" ||
405   - (*this)["inobs"].filename() == "") {
406   - std::string msg = "A valid file needs to be specified for the "
407   - "\"inobs\" parameter, yet \""+
408   - (*this)["inobs"].filename()+"\" was given."
409   - " Specify a vaild observation definition or "
410   - "event list FITS file to proceed.";
411   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
412   - }
  401 + // Throw exception if no input observation file is given
  402 + require_inobs(G_GET_PARAMETERS);
413 403  
414   - // Build observation container without response (not needed)
  404 + // Get observation container without response (not needed)
415 405 m_obs = get_observations(false);
416 406  
417 407 } // endif: there was no observation in the container
418 408  
419   - // Get parameters
420   - m_usepnt = (*this)["usepnt"].boolean();
421   -
422   - // Initialise event cube
423   - GCTAEventCube cube;
424   -
425   - // Check if pointing should be used
426   - if (m_usepnt) {
427   -
428   - // Build cube from user parameters and pointing information
429   - cube = get_cube(m_obs);
  409 + // Create an event cube based on task parameters
  410 + GCTAEventCube cube = create_cube(m_obs);
430 411  
431   - } // endif: pointing used as map center
432   -
433   - else {
434   -
435   - // Build cube from user parameters
436   - cube = get_cube();
437   -
438   - } // endelse: m_usepnt was false
439   -
440   - // Get the skymap and initialise the pixels
  412 + // Get the skymap from the cube and initialise all pixels to zero
441 413 m_cube = cube.map();
442 414 m_cube = 0.0;
443 415  
... ... @@ -454,6 +426,7 @@ void ctbin::get_parameters(void)
454 426 return;
455 427 }
456 428  
  429 +
457 430 /***********************************************************************//**
458 431 * @brief Fill events into counts cube
459 432 *
... ...
src/ctbkgcube/ctbkgcube.cpp
  View file @c2675d4
... ... @@ -510,7 +510,6 @@ void ctbkgcube::free_members(void)
510 510 * @brief Get application parameters
511 511 *
512 512 * @exception GException::invalid_value
513   - * Parameter "inobs" is required for ctbkgcube.
514 513 * No background model definition XML file specified.
515 514 *
516 515 * Get all task parameters from parameter file or (if required) by querying
... ... @@ -518,20 +517,12 @@ void ctbkgcube::free_members(void)
518 517 ***************************************************************************/
519 518 void ctbkgcube::get_parameters(void)
520 519 {
521   - // If there are no observations in container then load them via user parameters
  520 + // If there are no observations in container then load them via user
  521 + // parameters
522 522 if (m_obs.size() == 0) {
523 523  
524   - // Throw exception if no infile is given, since this tool needs an observation
525   - // including events
526   - if ((*this)["inobs"].filename() == "NONE" ||
527   - (*this)["inobs"].filename() == "") {
528   - std::string msg = "A valid file needs to be specified for the "
529   - "\"inobs\" parameter, yet \""+
530   - (*this)["inobs"].filename()+"\" was given."
531   - " Specify a vaild observation definition or "
532   - "event list FITS file to proceed.";
533   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
534   - }
  524 + // Throw exception if no input observation file is given
  525 + require_inobs(G_GET_PARAMETERS);
535 526  
536 527 // Build observation container
537 528 m_obs = get_observations();
... ... @@ -561,31 +552,16 @@ void ctbkgcube::get_parameters(void)
561 552 if ((gammalib::toupper(incube) == "NONE") ||
562 553 (gammalib::strip_whitespace(incube) == "")) {
563 554  
564   - // Check for usepnt flag
565   - bool usepnt = (*this)["usepnt"].boolean();
566   -
567   - // Check if pointing should be used
568   - if (usepnt) {
569   -
570   - // build cube from user parameters and pointing information
571   - m_bkgcube = get_cube(m_obs);
572   -
573   - } // endif: pointing used as map center
574   -
575   - else {
576   -
577   - // Build cube from user parameters
578   - m_bkgcube = get_cube();
579   -
580   - } // endelse: m_usepnt was false
  555 + // Create an event cube based on task parameters
  556 + m_bkgcube = create_cube(m_obs);
581 557  
582 558 } // endif: filename was not valid
583 559  
584 560 // ... otherwise setup the background cube from the counts map
585 561 else {
586 562  
587   - // Set background cube from counts map
588   - m_bkgcube = set_from_cntmap(incube);
  563 + // Load event cube from filename
  564 + m_bkgcube.load(incube);
589 565  
590 566 } // endelse: cube was loaded from file
591 567  
... ...
src/ctbutterfly/ctbutterfly.cpp
  View file @c2675d4
... ... @@ -453,7 +453,6 @@ void ctbutterfly::free_members(void)
453 453 * @brief Get application parameters
454 454 *
455 455 * @exception GException::invalid_value
456   - * Parameter "inobs" is required for ctbutterfly.
457 456 * Test source not found or no RA/DEC parameters found for test
458 457 * source.
459 458 *
... ... @@ -465,19 +464,12 @@ void ctbutterfly::free_members(void)
465 464 ***************************************************************************/
466 465 void ctbutterfly::get_parameters(void)
467 466 {
468   - // If there are no observations in container then load them via user parameters
  467 + // If there are no observations in container then load them via user
  468 + // parameters
469 469 if (m_obs.size() == 0) {
470 470  
471   - // Throw exception if no infile is given
472   - if ((*this)["inobs"].filename() == "NONE" ||
473   - (*this)["inobs"].filename() == "") {
474   - std::string msg = "A valid file needs to be specified for the "
475   - "\"inobs\" parameter, yet \""+
476   - (*this)["inobs"].filename()+"\" was given."
477   - " Specify a vaild observation definition or "
478   - "event list FITS file to proceed.";
479   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
480   - }
  471 + // Throw exception if no input observation file is given
  472 + require_inobs(G_GET_PARAMETERS);
481 473  
482 474 // Build observation container
483 475 m_obs = get_observations();
... ... @@ -506,8 +498,8 @@ void ctbutterfly::get_parameters(void)
506 498 throw GException::invalid_value(G_GET_PARAMETERS, msg);
507 499 }
508 500  
509   - // Get energy binning information
510   - m_ebounds = get_ebounds();
  501 + // Create energy boundaries from user parameters
  502 + m_ebounds = create_ebounds();
511 503  
512 504 // Get matrix file name and load if possible
513 505 std::string matrixfilename = (*this)["matrix"].filename();
... ...
src/ctexpcube/ctexpcube.cpp
  View file @c2675d4
... ... @@ -350,13 +350,9 @@ void ctexpcube::free_members(void)
350 350  
351 351 /***********************************************************************//**
352 352 * @brief Get application parameters
353   - * @exception GException::invalid_value
354   - * Parameter "inobs" is required for ctexpcube.
355 353 *
356 354 * Get all task parameters from parameter file or (if required) by querying
357 355 * the user. The parameters are read in the correct order.
358   - *
359   - * @todo Setup exposure cube from counts map
360 356 ***************************************************************************/
361 357 void ctexpcube::get_parameters(void)
362 358 {
... ... @@ -364,16 +360,8 @@ void ctexpcube::get_parameters(void)
364 360 // parameters
365 361 if (m_obs.size() == 0) {
366 362  
367   - // Throw exception if no infile is given
368   - if ((*this)["inobs"].filename() == "NONE" ||
369   - (*this)["inobs"].filename() == "") {
370   - std::string msg = "A valid file needs to be specified for the "
371   - "\"inobs\" parameter, yet \""+
372   - (*this)["inobs"].filename()+"\" was given."
373   - " Specify a vaild observation definition or "
374   - "event list FITS file to proceed.";
375   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
376   - }
  363 + // Throw exception if no input observation file is given
  364 + require_inobs(G_GET_PARAMETERS);
377 365  
378 366 // Build observation container
379 367 m_obs = get_observations();
... ... @@ -388,27 +376,9 @@ void ctexpcube::get_parameters(void)
388 376  
389 377 // Check for filename validity
390 378 if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
391   -
392   - // Initialise event cube
393   - GCTAEventCube cube;
394   -
395   - // Check for usepnt flag
396   - bool usepnt = (*this)["usepnt"].boolean();
397   -
398   - // Check if pointing should be used
399   - if (usepnt) {
400   -
401   - // build cube from user parameters and pointing information
402   - cube = get_cube(m_obs);
403   -
404   - } // endif: pointing used as map center
405   -
406   - else {
407   -
408   - // Build cube from user parameters
409   - cube = get_cube();
410 379  
411   - } // endelse: m_usepnt was false
  380 + // Create an event cube based on task parameters
  381 + GCTAEventCube cube = create_cube(m_obs);
412 382  
413 383 // Define exposure cube
414 384 m_expcube = GCTAExposure(cube);
... ... @@ -418,8 +388,8 @@ void ctexpcube::get_parameters(void)
418 388 // ... otherwise setup the exposure cube from the counts map
419 389 else {
420 390  
421   - // Set exposure cube from counts map
422   - GCTAEventCube cube = set_from_cntmap(incube);
  391 + // Load event cube from filename
  392 + GCTAEventCube cube(incube);
423 393  
424 394 // Define exposure cube
425 395 m_expcube = GCTAExposure(cube);
... ...
src/ctlike/ctlike.cpp
  View file @c2675d4
... ... @@ -341,9 +341,6 @@ void ctlike::save(void)
341 341 /***********************************************************************//**
342 342 * @brief Get application parameters
343 343 *
344   - * @exception GException::invalid_value
345   - * Parameter "inobs" is required for ctlike.
346   - *
347 344 * Get all required task parameters from the parameter file or (if specified)
348 345 * by querying the user. Observation dependent parameters will only be read
349 346 * if the observation container is actually empty. Observation dependent
... ... @@ -360,21 +357,14 @@ void ctlike::save(void)
360 357 ***************************************************************************/
361 358 void ctlike::get_parameters(void)
362 359 {
363   - // If there are no observations in container then load them via user parameters
  360 + // If there are no observations in container then load them via user
  361 + // parameters
364 362 if (m_obs.size() == 0) {
365 363  
366   - // Throw exception if no infile is given
367   - if ((*this)["inobs"].filename() == "NONE" ||
368   - (*this)["inobs"].filename() == "") {
369   - std::string msg = "A valid file needs to be specified for the "
370   - "\"inobs\" parameter, yet \""+
371   - (*this)["inobs"].filename()+"\" was given."
372   - " Specify a vaild observation definition or "
373   - "event list FITS file to proceed.";
374   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
375   - }
  364 + // Throw exception if no input observation file is given
  365 + require_inobs(G_GET_PARAMETERS);
376 366  
377   - // Build observation container
  367 + // Get observation container
378 368 m_obs = get_observations();
379 369  
380 370 } // endif: there was no observation in the container
... ...
src/ctmodel/ctmodel.cpp
  View file @c2675d4
... ... @@ -33,6 +33,7 @@
33 33 #include "GTools.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
  36 +#define G_GET_PARAMETERS "ctmodel::get_parameters()"
36 37  
37 38 /* __ Debug definitions __________________________________________________ */
38 39  
... ... @@ -411,13 +412,11 @@ void ctmodel::free_members(void)
411 412 ***************************************************************************/
412 413 void ctmodel::get_parameters(void)
413 414 {
414   - // If there are no observations in container then load them via user parameters
415   - if (m_obs.size() == 0) {
416   -
417   - // Build observation container
418   - m_obs = get_observations();
419   -
420   - } // endif: there was no observation in the container
  415 + // If there are no observations in container then load them via user
  416 + // parameters
  417 + if (m_obs.size() == 0) {
  418 + m_obs = get_observations();
  419 + }
421 420  
422 421 // Read model definition file if required
423 422 if (m_obs.models().size() == 0) {
... ... @@ -441,34 +440,18 @@ void ctmodel::get_parameters(void)
441 440 // Read cube definition file
442 441 std::string incube = (*this)["incube"].filename();
443 442  
444   - // If no cube file has been specified then read all parameters that
445   - // are necessary to create the cube from scratch (see method ctool::get_cube())
446   - if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
447   -
448   - // Get parameters
449   - bool usepnt = (*this)["usepnt"].boolean();
450   -
451   - // Initialise event cube
452   - GCTAEventCube cube;
453   -
454   - // Check if pointing should be used
455   - if (usepnt) {
456   -
457   - // buld cube from user parameters and pointing information
458   - m_cube = get_cube(m_obs);
459   -
460   - } // endif: pointing used as map center
461   -
462   - else {
463   -
464   - // Build cube from user parameters
465   - m_cube = get_cube();
466   -
467   - } // endelse: m_usepnt was false
  443 + // If no cube file has been specified then create a cube from
  444 + // the task parameters ...
  445 + if ((incube == "NONE") ||
  446 + (gammalib::strip_whitespace(incube) == "")) {
  447 +
  448 + // Create cube from scratch
  449 + m_cube = create_cube(m_obs);
468 450  
469   - } // endif: incube file was not vaild
  451 + }
470 452  
471   - // ... else: build map information from file
  453 + // ... otherwise load the cube from file and reset all bins
  454 + // to zero
472 455 else {
473 456  
474 457 // Load cube from given file
... ...
src/ctobssim/ctobssim.cpp
  View file @c2675d4
... ... @@ -559,13 +559,11 @@ void ctobssim::free_members(void)
559 559 ***************************************************************************/
560 560 void ctobssim::get_parameters(void)
561 561 {
562   - // If there are no observations in container then load them via user parameters
  562 + // If there are no observations in container then load them via user
  563 + // parameters
563 564 if (m_obs.size() == 0) {
564   -
565   - // Build observation container
566   - m_obs = get_observations();
567   -
568   - } // endif: there was no observation in the container
  565 + m_obs = get_observations();
  566 + }
569 567  
570 568 // Read model definition file if required
571 569 if (m_obs.models().size() == 0) {
... ... @@ -581,7 +579,7 @@ void ctobssim::get_parameters(void)
581 579 // dumped into the log file
582 580 if (read_ahead()) {
583 581 m_outevents = (*this)["outevents"].filename();
584   - m_prefix = (*this)["prefix"].string();
  582 + m_prefix = (*this)["prefix"].string();
585 583 }
586 584  
587 585 // Initialise random number generators. We initialise here one random
... ...
src/ctpsfcube/ctpsfcube.cpp
  View file @c2675d4
... ... @@ -341,29 +341,17 @@ void ctpsfcube::free_members(void)
341 341 /***********************************************************************//**
342 342 * @brief Get application parameters
343 343 *
344   - * @exception GException::invalid_value
345   - * Parameter "inobs" is required for ctbin.
346   - *
347 344 * Get all task parameters from parameter file or (if required) by querying
348 345 * the user. The parameters are read in the correct order.
349   - *
350   - * @todo Setup PSF cube from counts map
351 346 ***************************************************************************/
352 347 void ctpsfcube::get_parameters(void)
353 348 {
354   - // If there are no observations in container then load them via user parameters
  349 + // If there are no observations in container then load them via user
  350 + // parameters
355 351 if (m_obs.size() == 0) {
356 352  
357   - // Throw exception if no infile is given
358   - if ((*this)["inobs"].filename() == "NONE" ||
359   - (*this)["inobs"].filename() == "") {
360   - std::string msg = "A valid file needs to be specified for the "
361   - "\"inobs\" parameter, yet \""+
362   - (*this)["inobs"].filename()+"\" was given."
363   - " Specify a vaild observation definition or "
364   - "event list FITS file to proceed.";
365   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
366   - }
  353 + // Throw exception if no input observation file is given
  354 + require_inobs(G_GET_PARAMETERS);
367 355  
368 356 // Build observation container
369 357 m_obs = get_observations();
... ... @@ -374,46 +362,28 @@ void ctpsfcube::get_parameters(void)
374 362 std::string incube = (*this)["incube"].filename();
375 363  
376 364 // Get additional binning parameters
377   - double dmax = (*this)["amax"].real();
378   - int ndbins = (*this)["anumbins"].integer();
  365 + double amax = (*this)["amax"].real();
  366 + int anumbins = (*this)["anumbins"].integer();
379 367  
380 368 // Check for filename validity
381 369 if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) {
382 370  
383   - // Get parameters
384   - bool usepnt = (*this)["usepnt"].boolean();
385   -
386   - // Initialise event cube
387   - GCTAEventCube cube;
388   -
389   - // Check if pointing should be used
390   - if (usepnt) {
391   -
392   - // Build cube from user parameters and pointing information
393   - cube = get_cube(m_obs);
394   -
395   - } // endif: pointing used as map center
396   -
397   - else {
398   -
399   - // Build cube from user parameters
400   - cube = get_cube();
401   -
402   - } // endelse: m_usepnt was false
  371 + // Create an event cube based on task parameters
  372 + GCTAEventCube cube = create_cube(m_obs);
403 373  
404   - // Define psf cube
405   - m_psfcube = GCTAMeanPsf(cube, dmax, ndbins);
  374 + // Define psf cube
  375 + m_psfcube = GCTAMeanPsf(cube, amax, anumbins);
406 376  
407 377 }
408 378  
409 379 // ... otherwise setup the exposure cube from the counts map
410 380 else {
411 381  
412   - // Set exposure cube from counts map
413   - GCTAEventCube cube = set_from_cntmap(incube);
  382 + // Load event cube from filename
  383 + GCTAEventCube cube(incube);
414 384  
415   - // Define psf cube
416   - m_psfcube = GCTAMeanPsf(cube, dmax, ndbins);
  385 + // Define psf cube
  386 + m_psfcube = GCTAMeanPsf(cube, amax, anumbins);
417 387  
418 388 } // endelse: cube loaded from file
419 389  
... ...
src/ctselect/ctselect.cpp
  View file @c2675d4
... ... @@ -456,10 +456,7 @@ void ctselect::free_members(void)
456 456 /***********************************************************************//**
457 457 * @brief Get application parameters
458 458 *
459   - * @exception GException::invalid_value
460   - * Parameter "inobs" is required for ctselect.
461   - *
462   - * Get all task parameters from parameter file or (if required) by querying
  459 + * Get all user parameters from parameter file or (if required) by querying
463 460 * the user. Times are assumed to be in the native CTA MJD format.
464 461 *
465 462 * This method also loads observations if no observations are yet allocated.
... ... @@ -468,21 +465,14 @@ void ctselect::free_members(void)
468 465 ***************************************************************************/
469 466 void ctselect::get_parameters(void)
470 467 {
471   - // If there are no observations in container then load them via user parameters
  468 + // If there are no observations in container then load them via user
  469 + // parameters
472 470 if (m_obs.size() == 0) {
473 471  
474   - // Throw exception if no infile is given
475   - if ((*this)["inobs"].filename() == "NONE" ||
476   - (*this)["inobs"].filename() == "") {
477   - std::string msg = "A valid file needs to be specified for the "
478   - "\"inobs\" parameter, yet \""+
479   - (*this)["inobs"].filename()+"\" was given."
480   - " Specify a vaild observation definition or "
481   - "event list FITS file to proceed.";
482   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
483   - }
  472 + // Throw exception if no input observation file is given
  473 + require_inobs(G_GET_PARAMETERS);
484 474  
485   - // Build observation container without response (not needed)
  475 + // Get observation container without response (not needed)
486 476 m_obs = get_observations(false);
487 477  
488 478 } // endif: there was no observation in the container
... ...
src/ctskymap/ctskymap.cpp
  View file @c2675d4
... ... @@ -319,9 +319,6 @@ void ctskymap::save(void)
319 319 /***********************************************************************//**
320 320 * @brief Get application parameters
321 321 *
322   - * @exception GException::invalid_value
323   - * Parameter "inobs" is required for ctskymap.
324   - *
325 322 * Get all task parameters from parameter file or (if required) by querying
326 323 * the user. Most parameters are only required if no observation exists so
327 324 * far in the observation container. In this case, a single CTA observation
... ... @@ -334,39 +331,16 @@ void ctskymap::get_parameters(void)
334 331 // parameters
335 332 if (m_obs.size() == 0) {
336 333  
337   - // Throw exception if no infile is given
338   - if ((*this)["inobs"].filename() == "NONE" ||
339   - (*this)["inobs"].filename() == "") {
340   - std::string msg = "A valid file needs to be specified for the "
341   - "\"inobs\" parameter, yet \""+
342   - (*this)["inobs"].filename()+"\" was given."
343   - " Specify a vaild observation definition or "
344   - "event list FITS file to proceed.";
345   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
346   - }
  334 + // Throw exception if no input observation file is given
  335 + require_inobs(G_GET_PARAMETERS);
347 336  
348   - // Build observation container without response (not needed)
  337 + // Get observation container without response (not needed)
349 338 m_obs = get_observations(false);
350 339  
351 340 } // endif: there was no observation in the container
352 341  
353   - // Get parameters
354   - bool usepnt = (*this)["usepnt"].boolean();
355   -
356   - // Check if pointing should be used
357   - if (usepnt) {
358   -
359   - // buld map from user parameters and pointing information
360   - m_skymap = get_map(m_obs);
361   -
362   - } // endif: pointing used as map center
363   -
364   - else {
365   -
366   - // Build cube from user parameters
367   - m_skymap = get_map();
368   -
369   - } // endelse: m_usepnt was false
  342 + // Create sky map based on task parameters
  343 + m_skymap = create_map(m_obs);
370 344  
371 345 // Get remaining parameters
372 346 m_emin = (*this)["emin"].real();
... ... @@ -381,6 +355,7 @@ void ctskymap::get_parameters(void)
381 355 return;
382 356 }
383 357  
  358 +
384 359 /***********************************************************************//**
385 360 * @brief Map events into a sky map
386 361 *
... ...
src/cttsmap/cttsmap.cpp
  View file @c2675d4
... ... @@ -470,9 +470,6 @@ void cttsmap::free_members(void)
470 470 /***********************************************************************//**
471 471 * @brief Get application parameters
472 472 *
473   - * @exception GException::invalid_value
474   - * Parameter "inobs" is required for cttsmap.
475   - *
476 473 * @exception GException::invalid_value
477 474 * Test source not found or no RA/DEC parameters found for test
478 475 * source.
... ... @@ -489,18 +486,10 @@ void cttsmap::get_parameters(void)
489 486 // parameters
490 487 if (m_obs.size() == 0) {
491 488  
492   - // Throw exception if no infile is given
493   - if ((*this)["inobs"].filename() == "NONE" ||
494   - (*this)["inobs"].filename() == "") {
495   - std::string msg = "A valid file needs to be specified for the "
496   - "\"inobs\" parameter, yet \""+
497   - (*this)["inobs"].filename()+"\" was given."
498   - " Specify a vaild observation definition or "
499   - "event list FITS file to proceed.";
500   - throw GException::invalid_value(G_GET_PARAMETERS, msg);
501   - }
  489 + // Throw exception if no input observation file is given
  490 + require_inobs(G_GET_PARAMETERS);
502 491  
503   - // Build observation container
  492 + // Get observation container
504 493 m_obs = get_observations();
505 494  
506 495 } // endif: there was no observation in the container
... ... @@ -547,23 +536,8 @@ void cttsmap::get_parameters(void)
547 536 (*m_testsource)["DEC"].fix();
548 537 }
549 538  
550   - // Get sky map binning from user parameters and check if pointing
551   - // should be used
552   - GSkymap map;
553   - bool usepnt = (*this)["usepnt"].boolean();
554   - if (usepnt) {
555   -
556   - // Build cube from user parameters and pointing information
557   - map = get_map(m_obs);
558   -
559   - } // endif: pointing used as map center
560   -
561   - else {
562   -
563   - // Build cube from user parameters
564   - map = get_map();
565   -
566   - } // endelse: m_usepnt was false
  539 + // Create sky map based on task parameters
  540 + GSkymap map = create_map(m_obs);
567 541  
568 542 // Initialise maps from user parameters
569 543 init_maps(map);
... ...
src/support/ctool.cpp
  View file @c2675d4
... ... @@ -47,10 +47,8 @@
47 47  
48 48  
49 49 /* __ Method name definitions ____________________________________________ */
50   -#define G_GET_EBOUNDS "ctool::get_ebounds()"
51   -#define G_SET_FROM_CNTMAP "ctool::set_from_cntmap(std::string&)"
52 50 #define G_GET_MEAN_POINTING "ctool::get_mean_pointing(GObservations&)"
53   -#define G_GET_CUBE "ctool::get_cube(GObservations&)"
  51 +#define G_CREATE_EBOUNDS "ctool::create_ebounds()"
54 52  
55 53 /* __ Debug definitions __________________________________________________ */
56 54  
... ... @@ -277,255 +275,6 @@ void ctool::free_members(void)
277 275  
278 276  
279 277 /***********************************************************************//**
280   - * @brief Build a CTA observation with an empty event list
281   - *
282   - * Creates an empty CTA observation from user parameters. An empty event list
283   - * including RoI, GTI and energy boundary information is attached to the
284   - * observation. The method also sets the pointing direction using the ra and
285   - * dec parameter, the ROI based on ra, dec and rad, a single GTI based on
286   - * tmin and tmax, and a single energy boundary based on emin and emax. The
287   - * method furthermore sets the ontime, livetime and deadtime correction
288   - * factor.
289   - *
290   - * The following parameters are read:
291   - * ra - Right Ascension of pointing and RoI centre (deg)
292   - * dec - Declination of pointing and RoI centre (deg)
293   - * rad - Radius of RoI (deg)
294   - * deadc - Deadtime correction factor
295   - * tmin - Start time
296   - * tmax - Stop time
297   - * emin - Minimum energy (TeV)
298   - * emax - Maximum energy (TeV)
299   - ***************************************************************************/
300   -GCTAObservation ctool::setup_obs(void)
301   -{
302   - // Get CTA observation parameters
303   - double ra = (*this)["ra"].real();
304   - double dec = (*this)["dec"].real();
305   - double rad = (*this)["rad"].real();
306   - double deadc = (*this)["deadc"].real();
307   - double t_min = (*this)["tmin"].real();
308   - double t_max = (*this)["tmax"].real();
309   - double e_min = (*this)["emin"].real();
310   - double e_max = (*this)["emax"].real();
311   -
312   - // Allocate CTA observation and empty event list
313   - GCTAObservation obs;
314   - GCTAEventList list;
315   -
316   - // Set pointing direction
317   - GCTAPointing pnt;
318   - GSkyDir skydir;
319   - skydir.radec_deg(ra, dec);
320   - pnt.dir(skydir);
321   -
322   - // Set ROI
323   - GCTAInstDir instdir(skydir);
324   - GCTARoi roi(instdir, rad);
325   -
326   - // Set GTI
327   - GGti gti(m_cta_ref);
328   - GTime tstart;
329   - GTime tstop;
330   - tstart.set(t_min, m_cta_ref);
331   - tstop.set(t_max, m_cta_ref);
332   - gti.append(tstart, tstop);
333   -
334   - // Set energy boundaries
335   - GEbounds ebounds;
336   - GEnergy emin;
337   - GEnergy emax;
338   - emin.TeV(e_min);
339   - emax.TeV(e_max);
340   - ebounds.append(emin, emax);
341   -
342   - // Set CTA event list attributes
343   - list.roi(roi);
344   - list.gti(gti);
345   - list.ebounds(ebounds);
346   -
347   - // Attach event list to CTA observation
348   - obs.events(list);
349   -
350   - // Set observation ontime, livetime and deadtime correction factor
351   - obs.pointing(pnt);
352   - obs.ontime(gti.ontime());
353   - obs.livetime(gti.ontime()*deadc);
354   - obs.deadc(deadc);
355   -
356   - // Return CTA observation
357   - return obs;
358   -
359   -}
360   -
361   -
362   -/***********************************************************************//**
363   - * @brief Build a skymap based on user parameters
364   - *
365   - * Creates an empty skymap from input user parameters.
366   - *
367   - * The following parameters are read:
368   - * proj - Coordinate projection
369   - * coordsys - Coordinate system
370   - * binsz - Bin size (deg)
371   - * nxpix - Number of pixels in Right Ascension or longitude
372   - * nypix - Number of pixels in Declination or latitude
373   - ***************************************************************************/
374   -GSkymap ctool::get_map()
375   -{
376   - // Read task parameters
377   - double xref = (*this)["xref"].real();
378   - double yref = (*this)["yref"].real();
379   - std::string proj = (*this)["proj"].string();
380   - std::string coordsys = (*this)["coordsys"].string();
381   - double binsz = (*this)["binsz"].real();
382   - int nxpix = (*this)["nxpix"].integer();
383   - int nypix = (*this)["nypix"].integer();
384   -
385   - // Initialise sky map
386   - GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz,
387   - nxpix, nypix, 1);
388   -
389   - // Return sky map
390   - return map;
391   -}
392   -
393   -
394   -/***********************************************************************//**
395   - * @brief Build a skymap based on user parameters
396   - *
397   - * @param[in] obs Observation container.
398   - *
399   - * Creates an empty skymap from input user parameters. The reference
400   - * coordinate is extracted from the mean pointing of the CTA observations
401   - * in the observation container (see get_mean_pointing()).
402   - *
403   - * The following parameters are read:
404   - * proj - Coordinate projection
405   - * coordsys - Coordinate system
406   - * binsz - Bin size (deg)
407   - * nxpix - Number of pixels in Right Ascension or longitude
408   - * nypix - Number of pixels in Declination or latitude
409   - ***************************************************************************/
410   -GSkymap ctool::get_map(const GObservations& obs)
411   -{
412   - // Read task parameters
413   - std::string proj = (*this)["proj"].string();
414   - std::string coordsys = (*this)["coordsys"].string();
415   - double binsz = (*this)["binsz"].real();
416   - int nxpix = (*this)["nxpix"].integer();
417   - int nypix = (*this)["nypix"].integer();
418   -
419   - // Initialise reference coordinates
420   - double xref = 0.0;
421   - double yref = 0.0;
422   -
423   - // Get pointing from observations
424   - GSkyDir pnt = get_mean_pointing(obs);
425   -
426   - // Use correct coordinate system
427   - if (gammalib::toupper(coordsys) == "GAL") {
428   - xref = pnt.l_deg();
429   - yref = pnt.b_deg();
430   - }
431   - else {
432   - xref = pnt.ra_deg();
433   - yref = pnt.dec_deg();
434   - }
435   -
436   - // Initialise sky map
437   - GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz,
438   - nxpix, nypix, 1);
439   -
440   - // Return sky map
441   - return map;
442   -}
443   -
444   -
445   -/***********************************************************************//**
446   - * @brief Get the energy boundaries
447   - *
448   - * @exception GException::invalid_value
449   - * No valid energy boundary extension found.
450   - *
451   - * Get the energy boundaries according to the user parameters. The method
452   - * supports loading of energy boundary information from the EBOUNDS or
453   - * ENERGYBINS extension, or setting energy boundaries using a linear
454   - * or logarithmical spacing.
455   - *
456   - * The following parameters are read:
457   - * ebinalg - Energy binning algorithm
458   - * emin - Minimum energy (if ebinalg != FILE)
459   - * emax - Maximum energy (if ebinalg != FILE)
460   - * enumbins - Number of energy bins (if ebinalg != FILE)
461   - ***************************************************************************/
462   -GEbounds ctool::get_ebounds(void)
463   -{
464   - // Allocate energy boundaries
465   - GEbounds ebounds;
466   -
467   - // Get energy binning algorithm
468   - std::string ebinalg = (*this)["ebinalg"].string();
469   -
470   - // If energy binning algorithm is of type "FILE" (case sensitive), then
471   - // read energy boundaries from FITS file ...
472   - if (ebinalg == "FILE") {
473   -
474   - // Get filename
475   - std::string ebinfile = (*this)["ebinfile"].filename();
476   -
477   - // Open energy boundary file using the EBOUNDS or ENERGYBINS
478   - // extension. Throw an exception if opening fails.
479   - GFits file(ebinfile);
480   - if (file.contains("EBOUNDS")) {
481   - file.close();
482   - ebounds.load(ebinfile,"EBOUNDS");
483   - }
484   - else if (file.contains("ENERGYBINS")) {
485   - file.close();
486   - ebounds.load(ebinfile,"ENERGYBINS");
487   - }
488   - else {
489   - file.close();
490   - std::string msg = "No extension with name \"EBOUNDS\" or"
491   - " \"ENERGYBINS\" found in FITS file"
492   - " \""+ebinfile+"\".\n"
493   - "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
494   - " is required if the parameter \"ebinalg\""
495   - " is set to \"FILE\".";
496   - throw GException::invalid_value(G_GET_EBOUNDS, msg);
497   - }
498   -
499   - } // endif: ebinalg was "FILE"
500   -
501   - // ... otherwise use a linear or a logarithmically-spaced energy binning
502   - else {
503   -
504   - // Get task parameters
505   - double emin = (*this)["emin"].real();
506   - double emax = (*this)["emax"].real();
507   - int enumbins = (*this)["enumbins"].integer();
508   -
509   - // Initialise log mode for ebinning
510   - bool log = true;
511   -
512   - // Check if algorithm is linear
513   - if (ebinalg == "LIN") {
514   - log = false;
515   - }
516   -
517   - // Setup energy bins
518   - ebounds = GEbounds(enumbins, GEnergy(emin, "TeV"),
519   - GEnergy(emax, "TeV"), log);
520   -
521   - } // endelse: ebinalg was not "FILE"
522   -
523   - // Return energy boundaries
524   - return ebounds;
525   -}
526   -
527   -
528   -/***********************************************************************//**
529 278 * @brief Get observation container
530 279 *
531 280 * @param[in] get_response Indicates whether response information should
... ... @@ -546,11 +295,10 @@ GObservations ctool::get_observations(const bool& get_response)
546 295  
547 296 // If no observation definition file has been specified then read all
548 297 // parameters that are necessary to create an observation from scratch
549   - // (see method setup_obs)
550 298 if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) {
551 299  
552 300 // Setup a new CTA observation
553   - GCTAObservation cta_obs = setup_obs();
  301 + GCTAObservation cta_obs = create_cta_obs();
554 302  
555 303 // Get response if required
556 304 if (get_response) {
... ... @@ -622,110 +370,152 @@ GObservations ctool::get_observations(const bool& get_response)
622 370  
623 371  
624 372 /***********************************************************************//**
625   - * @brief Derives mean pointing from CTA observations
  373 + * @brief Create energy boundaries from user parameters
626 374 *
627   - * @param[in] obs Observation container.
628   - * @return Pointing direction.
  375 + * @exception GException::invalid_value
  376 + * No valid energy boundary extension found.
629 377 *
630   - * Computes the mean pointing direction from all CTA observations in the
631   - * observation container.
  378 + * Get the energy boundaries according to the user parameters. The method
  379 + * supports loading of energy boundary information from the EBOUNDS or
  380 + * ENERGYBINS extension, or setting energy boundaries using a linear
  381 + * or logarithmical spacing.
632 382 *
633   - * @todo This method does not work properly for wrap arounds in Right
634   - * Ascension. It certainly will also lead to bad results near the
635   - * celestial pole.
  383 + * The following parameters are read:
  384 + * ebinalg - Energy binning algorithm
  385 + * emin - Minimum energy (if ebinalg != FILE)
  386 + * emax - Maximum energy (if ebinalg != FILE)
  387 + * enumbins - Number of energy bins (if ebinalg != FILE)
636 388 ***************************************************************************/
637   -GSkyDir ctool::get_mean_pointing(const GObservations& obs)
  389 +GEbounds ctool::create_ebounds(void)
638 390 {
639   - // Initialise values
640   - double xref = 0.0;
641   - double yref = 0.0;
642   - int n_pnt = 0;
  391 + // Allocate energy boundaries
  392 + GEbounds ebounds;
643 393  
644   - // Loop over all observations to compute mean pointings
645   - for (int i = 0; i < obs.size(); ++i) {
  394 + // Get energy binning algorithm
  395 + std::string ebinalg = (*this)["ebinalg"].string();
646 396  
647   - // If we have a CTA observation then extract the pointing direction
648   - const GCTAObservation* cta_obs = dynamic_cast<const GCTAObservation*>(obs[i]);
649   - if (cta_obs != NULL) {
  397 + // If energy binning algorithm is of type "FILE" (case sensitive), then
  398 + // read energy boundaries from FITS file ...
  399 + if (ebinalg == "FILE") {
650 400  
651   - // Get pointing
652   - const GCTAPointing& pnt = cta_obs->pointing();
  401 + // Get filename
  402 + std::string ebinfile = (*this)["ebinfile"].filename();
653 403  
654   - // Add to coordinates
655   - xref += pnt.dir().ra_deg();
656   - yref += pnt.dir().dec_deg();
657   - n_pnt += 1;
  404 + // Open energy boundary file using the EBOUNDS or ENERGYBINS
  405 + // extension. Throw an exception if opening fails.
  406 + GFits file(ebinfile);
  407 + if (file.contains("EBOUNDS")) {
  408 + file.close();
  409 + ebounds.load(ebinfile,"EBOUNDS");
  410 + }
  411 + else if (file.contains("ENERGYBINS")) {
  412 + file.close();
  413 + ebounds.load(ebinfile,"ENERGYBINS");
  414 + }
  415 + else {
  416 + file.close();
  417 + std::string msg = "No extension with name \"EBOUNDS\" or"
  418 + " \"ENERGYBINS\" found in FITS file"
  419 + " \""+ebinfile+"\".\n"
  420 + "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
  421 + " is required if the parameter \"ebinalg\""
  422 + " is set to \"FILE\".";
  423 + throw GException::invalid_value(G_CREATE_EBOUNDS, msg);
  424 + }
  425 +
  426 + } // endif: ebinalg was "FILE"
658 427  
659   - } // endif: cta observation was valid
  428 + // ... otherwise use a linear or a logarithmically-spaced energy binning
  429 + else {
660 430  
661   - } // endfor: loop over all observations
  431 + // Get task parameters
  432 + double emin = (*this)["emin"].real();
  433 + double emax = (*this)["emax"].real();
  434 + int enumbins = (*this)["enumbins"].integer();
662 435  
663   - // Signal if no pointing is found
664   - if (n_pnt < 1) {
665   - std::string msg = "No valid CTA observation has been found in "
666   - "observation list, hence no pointing information "
667   - "could be extracted. Use the \"usepnt=no\" "
668   - "option and specify the pointing explicitly.";
669   - throw GException::invalid_value(G_GET_MEAN_POINTING, msg);
670   - }
  436 + // Initialise log mode for ebinning
  437 + bool log = true;
671 438  
672   - // Calculate mean coordinates
673   - double ra = xref / double(n_pnt);
674   - double dec = yref / double(n_pnt);
  439 + // Check if algorithm is linear
  440 + if (ebinalg == "LIN") {
  441 + log = false;
  442 + }
675 443  
676   - // Initialise sky dir
677   - GSkyDir dir = GSkyDir();
678   - dir.radec_deg(ra,dec);
  444 + // Setup energy bins
  445 + ebounds = GEbounds(enumbins, GEnergy(emin, "TeV"),
  446 + GEnergy(emax, "TeV"), log);
679 447  
680   - // Return sky direction
681   - return dir;
  448 + } // endelse: ebinalg was not "FILE"
  449 +
  450 + // Return energy boundaries
  451 + return ebounds;
682 452 }
683 453  
684 454  
685 455 /***********************************************************************//**
686   - * @brief Build CTA event cube based on user parameters
  456 + * @brief Create a skymap from user parameters
  457 + *
  458 + * @param[in] obs Observation container.
687 459 *
688   - * Creates an empty CTA event cube from input user parameters. It appends a
689   - * dummy GTI to the event cube.
  460 + * Creates an empty skymap from input user parameters. The reference
  461 + * coordinate is extracted from the mean pointing of the CTA observations
  462 + * in the observation container (see get_mean_pointing()).
690 463 *
691 464 * The following parameters are read:
692   - * xref - Reference in Right Ascension or longitude (deg)
693   - * yref - Reference in Declination or latitude (deg)
  465 + * usepnt - Use pointing for reference coordinate?
  466 + * xref - Reference in Right Ascension or longitude (if usepnt=no)
  467 + * yref - Reference in Declination or latitude (if usepnt=no)
694 468 * proj - Coordinate projection
695 469 * coordsys - Coordinate system
696 470 * binsz - Bin size (deg)
697 471 * nxpix - Number of pixels in Right Ascension or longitude
698 472 * nypix - Number of pixels in Declination or latitude
699   - * ebinalg - Energy binning algorithm
700   - * emin - Minimum energy (if ebinalg != FILE)
701   - * emax - Maximum energy (if ebinalg != FILE)
702   - * enumbins - Number of energy bins (if ebinalg != FILE)
703 473 ***************************************************************************/
704   -GCTAEventCube ctool::get_cube(void)
  474 +GSkymap ctool::create_map(const GObservations& obs)
705 475 {
706   - // Get skymap
707   - GSkymap map = get_map();
  476 + // Read task parameters
  477 + double xref = 0.0;
  478 + double yref = 0.0;
  479 + bool usepnt = (*this)["usepnt"].boolean();
  480 + if (!usepnt) {
  481 + xref = (*this)["xref"].real();
  482 + yref = (*this)["yref"].real();
  483 + }
  484 + std::string proj = (*this)["proj"].string();
  485 + std::string coordsys = (*this)["coordsys"].string();
  486 + double binsz = (*this)["binsz"].real();
  487 + int nxpix = (*this)["nxpix"].integer();
  488 + int nypix = (*this)["nypix"].integer();
708 489  
709   - // Set energy boundaries
710   - GEbounds ebounds = get_ebounds();
  490 + // If requested, get pointing from observations
  491 + if (usepnt) {
  492 +
  493 + // Get mean pointing
  494 + GSkyDir pnt = get_mean_pointing(obs);
711 495  
712   - // Set dummy GTI that is needed for event cube creation
713   - GGti gti;
714   - gti.append(GTime(0.0), GTime(0.1234));
  496 + // Use correct coordinate system
  497 + if (gammalib::toupper(coordsys) == "GAL") {
  498 + xref = pnt.l_deg();
  499 + yref = pnt.b_deg();
  500 + }
  501 + else {
  502 + xref = pnt.ra_deg();
  503 + yref = pnt.dec_deg();
  504 + }
715 505  
716   - // Extend skymap to the requested number of maps
717   - map.nmaps(ebounds.size());
  506 + } // endif: got mean pointing as reference
718 507  
719   - // Build event cube from given information
720   - GCTAEventCube cube = GCTAEventCube(map, ebounds, gti);
  508 + // Initialise sky map
  509 + GSkymap map = GSkymap(proj, coordsys, xref, yref, -binsz, binsz,
  510 + nxpix, nypix, 1);
721 511  
722   - // Return cube
723   - return cube;
  512 + // Return sky map
  513 + return map;
724 514 }
725 515  
726 516  
727 517 /***********************************************************************//**
728   - * @brief Build CTA event cube based on user parameters
  518 + * @brief Create a CTA event cube from user parameters
729 519 *
730 520 * @param[in] obs Observation container.
731 521 *
... ... @@ -746,13 +536,13 @@ GCTAEventCube ctool::get_cube(void)
746 536 * emax - Maximum energy (if ebinalg != FILE)
747 537 * enumbins - Number of energy bins (if ebinalg != FILE)
748 538 ***************************************************************************/
749   -GCTAEventCube ctool::get_cube(const GObservations& obs)
  539 +GCTAEventCube ctool::create_cube(const GObservations& obs)
750 540 {
751 541 // Get skymap
752   - GSkymap map = get_map(obs);
  542 + GSkymap map = create_map(obs);
753 543  
754 544 // Set energy boundaries
755   - GEbounds ebounds = get_ebounds();
  545 + GEbounds ebounds = create_ebounds();
756 546  
757 547 // Set dummy GTI that is needed for event cube creation
758 548 GGti gti;
... ... @@ -770,82 +560,108 @@ GCTAEventCube ctool::get_cube(const GObservations&amp; obs)
770 560  
771 561  
772 562 /***********************************************************************//**
773   - * @brief Get current resident set size (physical memory use) in Bytes
  563 + * @brief Create a CTA observation from user parameters
774 564 *
775   - * @return Physical memory use in Bytes.
  565 + * Creates an empty CTA observation from user parameters. An empty event list
  566 + * including RoI, GTI and energy boundary information is attached to the
  567 + * observation. The method also sets the pointing direction using the ra and
  568 + * dec parameter, the ROI based on ra, dec and rad, a single GTI based on
  569 + * tmin and tmax, and a single energy boundary based on emin and emax. The
  570 + * method furthermore sets the ontime, livetime and deadtime correction
  571 + * factor.
  572 + *
  573 + * The following parameters are read:
  574 + * ra - Right Ascension of pointing and RoI centre (deg)
  575 + * dec - Declination of pointing and RoI centre (deg)
  576 + * rad - Radius of RoI (deg)
  577 + * deadc - Deadtime correction factor
  578 + * tmin - Start time
  579 + * tmax - Stop time
  580 + * emin - Minimum energy (TeV)
  581 + * emax - Maximum energy (TeV)
776 582 ***************************************************************************/
777   -size_t ctool::get_current_rss(void)
  583 +GCTAObservation ctool::create_cta_obs(void)
778 584 {
779   - // Initialize resident set size
780   - size_t rss = 0;
  585 + // Get CTA observation parameters
  586 + double ra = (*this)["ra"].real();
  587 + double dec = (*this)["dec"].real();
  588 + double rad = (*this)["rad"].real();
  589 + double deadc = (*this)["deadc"].real();
  590 + double t_min = (*this)["tmin"].real();
  591 + double t_max = (*this)["tmax"].real();
  592 + double e_min = (*this)["emin"].real();
  593 + double e_max = (*this)["emax"].real();
781 594  
782   - // Determine resident set size (architecture dependent)
783   - // OSX
784   - #if defined(__APPLE__) && defined(__MACH__)
785   - #ifdef MACH_TASK_BASIC_INFO
786   - struct mach_task_basic_info info;
787   - mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
788   - if (task_info(mach_task_self( ), MACH_TASK_BASIC_INFO,
789   - (task_info_t)&info, &infoCount) == KERN_SUCCESS) {
790   - rss = (size_t)info.resident_size;
791   - }
792   - #else
793   - struct task_basic_info info;
794   - mach_msg_type_number_t info_count = TASK_BASIC_INFO_COUNT;
795   - if (task_info(mach_task_self(), TASK_BASIC_INFO,
796   - (task_info_t)&info, &info_count) == KERN_SUCCESS) {
797   - rss = (size_t)info.resident_size;
798   - }
799   - #endif
800   - // Linux
801   - #elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
802   - FILE* fp = NULL;
803   - if ((fp = fopen( "/proc/self/statm", "r" )) != NULL) {
804   - if (fscanf( fp, "%*s%ld", &rss ) == 1) {
805   - rss *= (size_t)sysconf(_SC_PAGESIZE);
806   - }
807   - fclose(fp);
808   - }
809   - // AIX, BSD, Solaris, and Unknown OS
810   - #else
811   - rss = 0;
812   - #endif
  595 + // Allocate CTA observation and empty event list
  596 + GCTAObservation obs;
  597 + GCTAEventList list;
  598 +
  599 + // Set pointing direction
  600 + GCTAPointing pnt;
  601 + GSkyDir skydir;
  602 + skydir.radec_deg(ra, dec);
  603 + pnt.dir(skydir);
  604 +
  605 + // Set ROI
  606 + GCTAInstDir instdir(skydir);
  607 + GCTARoi roi(instdir, rad);
  608 +
  609 + // Set GTI
  610 + GGti gti(m_cta_ref);
  611 + GTime tstart;
  612 + GTime tstop;
  613 + tstart.set(t_min, m_cta_ref);
  614 + tstop.set(t_max, m_cta_ref);
  615 + gti.append(tstart, tstop);
  616 +
  617 + // Set energy boundaries
  618 + GEbounds ebounds;
  619 + GEnergy emin;
  620 + GEnergy emax;
  621 + emin.TeV(e_min);
  622 + emax.TeV(e_max);
  623 + ebounds.append(emin, emax);
  624 +
  625 + // Set CTA event list attributes
  626 + list.roi(roi);
  627 + list.gti(gti);
  628 + list.ebounds(ebounds);
  629 +
  630 + // Attach event list to CTA observation
  631 + obs.events(list);
  632 +
  633 + // Set observation ontime, livetime and deadtime correction factor
  634 + obs.pointing(pnt);
  635 + obs.ontime(gti.ontime());
  636 + obs.livetime(gti.ontime()*deadc);
  637 + obs.deadc(deadc);
  638 +
  639 + // Return CTA observation
  640 + return obs;
813 641  
814   - // Return resident set size
815   - return rss;
816 642 }
817 643  
818 644  
819 645 /***********************************************************************//**
820   - * @brief Get a CTA ebent cube from from file
821   - *
822   - * @param[in] filename FITS file of count cube.
  646 + * @brief Throws exception in inobs parameter is not valid
823 647 *
824   - * @exception GException::invalid_value
825   - * Counts map not of WCS type
826   - *
827   - * Gets an event cube from a FITS file. This function is used in many ctools
828   - * to obtain a skymap definition from an event cube.
829   - *
830   - * @todo Is it necessary to restrict this method to WCS coordinates?
  648 + * @param[in] method Method name.
831 649 ***************************************************************************/
832   -GCTAEventCube ctool::set_from_cntmap(const std::string filename)
  650 +void ctool::require_inobs(const std::string& method)
833 651 {
834   - // Allocate CTA observation
835   - GCTAEventCube cube;
836   -
837   - // Load counts map in CTA observation
838   - cube.load(filename);
839   -
840   - // Check for valid WCS coordinate system
841   - const GWcs* wcs = dynamic_cast<const GWcs*>(cube.map().projection());
842   - if (wcs == NULL) {
843   - std::string msg = "Counts map project is not of WCS type.";
844   - throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
  652 + // Throw exception if no infile is given
  653 + if ((*this)["inobs"].filename() == "NONE" ||
  654 + (*this)["inobs"].filename() == "") {
  655 + std::string msg = "A valid file needs to be specified for the "
  656 + "\"inobs\" parameter, yet \""+
  657 + (*this)["inobs"].filename()+"\" was given."
  658 + " Specify a valid observation definition or "
  659 + "FITS file to proceed.";
  660 + throw GException::invalid_value(method, msg);
845 661 }
846 662  
847   - // Return cube
848   - return cube;
  663 + // Return
  664 + return;
849 665 }
850 666  
851 667  
... ... @@ -921,3 +737,111 @@ void ctool::set_obs_response(GCTAObservation* obs)
921 737 // Return
922 738 return;
923 739 }
  740 +
  741 +
  742 +/***********************************************************************//**
  743 + * @brief Derives mean pointing from CTA observations
  744 + *
  745 + * @param[in] obs Observation container.
  746 + * @return Pointing direction.
  747 + *
  748 + * Computes the mean pointing direction from all CTA observations in the
  749 + * observation container.
  750 + *
  751 + * @todo This method does not work properly for wrap arounds in Right
  752 + * Ascension. It certainly will also lead to bad results near the
  753 + * celestial pole.
  754 + ***************************************************************************/
  755 +GSkyDir ctool::get_mean_pointing(const GObservations& obs)
  756 +{
  757 + // Initialise values
  758 + double xref = 0.0;
  759 + double yref = 0.0;
  760 + int n_pnt = 0;
  761 +
  762 + // Loop over all observations to compute mean pointings
  763 + for (int i = 0; i < obs.size(); ++i) {
  764 +
  765 + // If we have a CTA observation then extract the pointing direction
  766 + const GCTAObservation* cta_obs = dynamic_cast<const GCTAObservation*>(obs[i]);
  767 + if (cta_obs != NULL) {
  768 +
  769 + // Get pointing
  770 + const GCTAPointing& pnt = cta_obs->pointing();
  771 +
  772 + // Add to coordinates
  773 + xref += pnt.dir().ra_deg();
  774 + yref += pnt.dir().dec_deg();
  775 + n_pnt += 1;
  776 +
  777 + } // endif: cta observation was valid
  778 +
  779 + } // endfor: loop over all observations
  780 +
  781 + // Signal if no pointing is found
  782 + if (n_pnt < 1) {
  783 + std::string msg = "No valid CTA observation has been found in "
  784 + "observation list, hence no pointing information "
  785 + "could be extracted. Use the \"usepnt=no\" "
  786 + "option and specify the pointing explicitly.";
  787 + throw GException::invalid_value(G_GET_MEAN_POINTING, msg);
  788 + }
  789 +
  790 + // Calculate mean coordinates
  791 + double ra = xref / double(n_pnt);
  792 + double dec = yref / double(n_pnt);
  793 +
  794 + // Initialise sky dir
  795 + GSkyDir dir = GSkyDir();
  796 + dir.radec_deg(ra,dec);
  797 +
  798 + // Return sky direction
  799 + return dir;
  800 +}
  801 +
  802 +
  803 +/***********************************************************************//**
  804 + * @brief Get current resident set size (physical memory use) in Bytes
  805 + *
  806 + * @return Physical memory use in Bytes.
  807 + ***************************************************************************/
  808 +size_t ctool::get_current_rss(void)
  809 +{
  810 + // Initialize resident set size
  811 + size_t rss = 0;
  812 +
  813 + // Determine resident set size (architecture dependent)
  814 + // OSX
  815 + #if defined(__APPLE__) && defined(__MACH__)
  816 + #ifdef MACH_TASK_BASIC_INFO
  817 + struct mach_task_basic_info info;
  818 + mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
  819 + if (task_info(mach_task_self( ), MACH_TASK_BASIC_INFO,
  820 + (task_info_t)&info, &infoCount) == KERN_SUCCESS) {
  821 + rss = (size_t)info.resident_size;
  822 + }
  823 + #else
  824 + struct task_basic_info info;
  825 + mach_msg_type_number_t info_count = TASK_BASIC_INFO_COUNT;
  826 + if (task_info(mach_task_self(), TASK_BASIC_INFO,
  827 + (task_info_t)&info, &info_count) == KERN_SUCCESS) {
  828 + rss = (size_t)info.resident_size;
  829 + }
  830 + #endif
  831 + // Linux
  832 + #elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
  833 + FILE* fp = NULL;
  834 + if ((fp = fopen( "/proc/self/statm", "r" )) != NULL) {
  835 + if (fscanf( fp, "%*s%ld", &rss ) == 1) {
  836 + rss *= (size_t)sysconf(_SC_PAGESIZE);
  837 + }
  838 + fclose(fp);
  839 + }
  840 + // AIX, BSD, Solaris, and Unknown OS
  841 + #else
  842 + rss = 0;
  843 + #endif
  844 +
  845 + // Return resident set size
  846 + return rss;
  847 +}
... ...
src/support/ctool.hpp
  View file @c2675d4
... ... @@ -66,18 +66,22 @@ protected:
66 66 void copy_members(const ctool& app);
67 67 void free_members(void);
68 68 const bool& read_ahead(void) const;
69   - GCTAObservation setup_obs(void);
70   - GSkymap get_map(void);
71   - GSkymap get_map(const GObservations& obs);
72   - GEbounds get_ebounds(void);
73 69 GObservations get_observations(const bool& get_response = true);
74   - GCTAEventCube get_cube(void);
75   - GCTAEventCube get_cube(const GObservations& obs);
76   - GSkyDir get_mean_pointing(const GObservations& obs);
77   - size_t get_current_rss(void);
78   - GCTAEventCube set_from_cntmap(const std::string filename);
  70 +
  71 + // Protected methods that create objects from user parameters
  72 + GEbounds create_ebounds(void);
  73 + GSkymap create_map(const GObservations& obs);
  74 + GCTAEventCube create_cube(const GObservations& obs);
  75 + GCTAObservation create_cta_obs(void);
  76 +
  77 + // Protected methods that check user parameters
  78 + void require_inobs(const std::string& method);
  79 +
  80 + // Protected support methods
79 81 void set_response(GObservations& obs);
80 82 void set_obs_response(GCTAObservation* obs);
  83 + GSkyDir get_mean_pointing(const GObservations& obs);
  84 + size_t get_current_rss(void);
81 85  
82 86 // Protected members
83 87 bool m_read_ahead; //!< Read ahead parameters
... ...