Olga Sergijenko / ctools

Commit b8e6842195111e67bf12aa401bdf47da50365e89

Authored by chiachun
1 parent 01c053e7

ctcubemask. add ctbin

Inline Side-by-side
Showing 6 changed files with 1096 additions and 2 deletions
pyext/ctcubemask.i 0 → 100644
  View file @b8e6842
  1 +/***************************************************************************
  2 + * ctcubemask - CTA mask cube tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 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 ctcubemask.i
  23 + * @brief CTA data selection tool Python interface definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "ctcubemask.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class ctcubemask
  35 + *
  36 + * @brief CTA data selection tool Python interface
  37 + ***************************************************************************/
  38 +class ctcubemask : public GApplication {
  39 +public:
  40 + // Constructors and destructors
  41 + ctcubemask(void);
  42 + explicit ctcubemask(GObservations obs);
  43 + ctcubemask(int argc, char *argv[]);
  44 + ctcubemask(const ctcubemask& app);
  45 + virtual ~ctcubemask(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + void execute(void);
  50 + void run(void);
  51 + void save(void);
  52 + GObservations& obs(void) { return m_obs; }
  53 + void get_parameters(void);
  54 + void apply_mask(GCTAObservation* obs, const std::string& filename);
  55 +};
  56 +
  57 +
  58 +/***********************************************************************//**
  59 + * @brief CTA data selection tool Python extension
  60 + ***************************************************************************/
  61 +%extend ctcubemask {
  62 + ctcubemask copy() {
  63 + return (*self);
  64 + }
  65 +}
... ...
src/Makefile.am
  View file @b8e6842
... ... @@ -11,7 +11,8 @@ SUBDIRS = ctbin \
11 11 ctobssim \
12 12 ctselect \
13 13 ctskymap \
14   - ctexpcube
  14 + ctexpcube \
  15 + ctcubemask
15 16  
16 17 # Scripts to be distributed
17 18 scripts = ctools-setup \
... ... @@ -29,7 +30,8 @@ libctools_la_LIBADD = ctbin/libctbin.la \
29 30 ctobssim/libctobssim.la \
30 31 ctselect/libctselect.la \
31 32 ctskymap/libctskymap.la \
32   - ctexpcube/libctexpcube.la
  33 + ctexpcube/libctexpcube.la \
  34 + ctcubemask/libctcubemask.la
33 35  
34 36 # Distribute scripts
35 37 EXTRA_DIST = $(scripts)
... ...
src/ctcubemask/ctcubemask.cpp 0 → 100644
  View file @b8e6842
  1 +/***************************************************************************
  2 + * ctfilter - CTA cube filter 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 ctcubemask.cpp
  23 + * @brief CTA mask cube implementation
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctcubemask.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_RUN "ctcubemask::run()"
  38 +#define G_APPLY_MASK "ctcubemask::apply_mask(GCTAObservation*,"\
  39 + "std::string&)"
  40 +
  41 +/* __ Debug definitions __________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +
  46 +/*==========================================================================
  47 + = =
  48 + = Constructors/destructors =
  49 + = =
  50 + ==========================================================================*/
  51 +
  52 +/***********************************************************************//**
  53 + * @brief Void constructor
  54 + ***************************************************************************/
  55 +ctcubemask::ctcubemask(void) : GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
  56 +{
  57 + // Initialise members
  58 + init_members();
  59 +
  60 + // Write header into logger
  61 + log_header();
  62 +
  63 + // Return
  64 + return;
  65 +}
  66 +
  67 +
  68 +/***********************************************************************//**
  69 + * @brief Observations constructor
  70 + *
  71 + * This constructor creates an instance of the class that is initialised from
  72 + * an observation container.
  73 + ***************************************************************************/
  74 +ctcubemask::ctcubemask(GObservations obs) : GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
  75 +{
  76 + // Initialise members
  77 + init_members();
  78 +
  79 + // Set observations
  80 + m_obs = obs;
  81 +
  82 + // Write header into logger
  83 + log_header();
  84 +
  85 + // Return
  86 + return;
  87 +}
  88 +
  89 +
  90 +
  91 +
  92 +/***********************************************************************//**
  93 + * @brief Command line constructor
  94 + *
  95 + * @param[in] argc Number of arguments in command line.
  96 + * @param[in] argv Array of command line arguments.
  97 + ***************************************************************************/
  98 +ctcubemask::ctcubemask(int argc, char *argv[]) :
  99 + GApplication(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv)
  100 +{
  101 + // Initialise members
  102 + init_members();
  103 +
  104 + // Write header into logger
  105 + log_header();
  106 +
  107 + // Return
  108 + return;
  109 +}
  110 +
  111 +
  112 +/***********************************************************************//**
  113 + * @brief Copy constructor
  114 + *
  115 + * @param[in] app Application.
  116 + ***************************************************************************/
  117 +ctcubemask::ctcubemask(const ctcubemask& app) : GApplication(app)
  118 +{
  119 + // Initialise members
  120 + init_members();
  121 +
  122 + // Copy members
  123 + copy_members(app);
  124 +
  125 + // Return
  126 + return;
  127 +}
  128 +
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Destructor
  132 + ***************************************************************************/
  133 +ctcubemask::~ctcubemask(void)
  134 +{
  135 + // Free members
  136 + free_members();
  137 +
  138 + // Return
  139 + return;
  140 +}
  141 +
  142 +
  143 +/*==========================================================================
  144 + = =
  145 + = Operators =
  146 + = =
  147 + ==========================================================================*/
  148 +
  149 +/***********************************************************************//**
  150 + * @brief Assignment operator
  151 + *
  152 + * @param[in] app Application.
  153 + ***************************************************************************/
  154 +ctcubemask& ctcubemask::operator= (const ctcubemask& app)
  155 +{
  156 + // Execute only if object is not identical
  157 + if (this != &app) {
  158 +
  159 + // Copy base class members
  160 + this->GApplication::operator=(app);
  161 +
  162 + // Free members
  163 + free_members();
  164 +
  165 + // Initialise members
  166 + init_members();
  167 +
  168 + // Copy members
  169 + copy_members(app);
  170 +
  171 + } // endif: object was not identical
  172 +
  173 + // Return this object
  174 + return *this;
  175 +}
  176 +
  177 +
  178 +/*==========================================================================
  179 + = =
  180 + = Public methods =
  181 + = =
  182 + ==========================================================================*/
  183 +
  184 +/***********************************************************************//**
  185 + * @brief Clear instance
  186 + ***************************************************************************/
  187 +void ctcubemask::clear(void)
  188 +{
  189 + // Free members
  190 + free_members();
  191 + this->GApplication::free_members();
  192 +
  193 + // Initialise members
  194 + this->GApplication::init_members();
  195 + init_members();
  196 +
  197 + // Return
  198 + return;
  199 +}
  200 +
  201 +
  202 +/***********************************************************************//**
  203 + * @brief Execute application
  204 + *
  205 + * This method performs the event selection and saves the result
  206 + ***************************************************************************/
  207 +void ctcubemask::execute(void)
  208 +{
  209 + // Signal that some parameters should be read ahead
  210 + m_read_ahead = true;
  211 +
  212 + // Perform event selection
  213 + run();
  214 +
  215 + // Save results
  216 + save();
  217 +
  218 + // Return
  219 + return;
  220 +}
  221 +
  222 +
  223 +
  224 +/***********************************************************************//**
  225 + * @brief mask data cube
  226 + *
  227 + * This method reads in the application parameters and loops over all
  228 + * observations that were found to apply a mask on the event cube.
  229 + ***************************************************************************/
  230 +void ctcubemask::run(void)
  231 +{
  232 + // Switch screen logging on in debug mode
  233 + if (logDebug()) {
  234 + log.cout(true);
  235 + }
  236 +
  237 + // Get parameters
  238 + get_parameters();
  239 +
  240 + // Write parameters into logger
  241 + if (logTerse()) {
  242 + log_parameters();
  243 + log << std::endl;
  244 + }
  245 +
  246 + // Write observation(s) into logger
  247 + if (logTerse()) {
  248 + log << std::endl;
  249 + log.header1("Observations for applying mask");
  250 + log << m_obs << std::endl;
  251 + }
  252 +
  253 + // Write header
  254 + if (logTerse()) {
  255 + log << std::endl;
  256 + log.header1("Apply mask");
  257 + }
  258 +
  259 + // Initialise counters
  260 + int n_observations = 0;
  261 +
  262 + // Loop over all observation in the container
  263 + for (int i = 0; i < m_obs.size(); ++i) {
  264 +
  265 + // Initialise event input and output filenames
  266 + m_infiles.push_back("");
  267 +
  268 + // Get CTA observation
  269 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  270 +
  271 + // Continue only if observation is a CTA observation
  272 + if (obs != NULL) {
  273 +
  274 + // Write header for observation
  275 + if (logTerse()) {
  276 + if (obs->name().length() > 1) {
  277 + log.header3("Observation "+obs->name());
  278 + }
  279 + else {
  280 + log.header3("Observation");
  281 + }
  282 + }
  283 +
  284 + // Increment counter
  285 + n_observations++;
  286 +
  287 + // Save event file name (for possible saving)
  288 + m_infiles[i] = obs->eventfile();
  289 +
  290 + // Get temporary file name
  291 + std::string filename = std::tmpnam(NULL);
  292 +
  293 + // Save observation in temporary file
  294 + obs->save(filename, true);
  295 +
  296 + // Log saved FITS file
  297 + if (logExplicit()) {
  298 + GFits tmpfile(filename);
  299 + log << std::endl;
  300 + log.header1("FITS file content of temporary file");
  301 + log << tmpfile << std::endl;
  302 + tmpfile.close();
  303 + }
  304 +
  305 + // Check temporary file
  306 + std::string message = check_infile(filename);
  307 + if (message.length() > 0) {
  308 + throw GException::app_error(G_RUN, message);
  309 + }
  310 +
  311 + // Load observation from temporary file, including event selection
  312 + apply_mask(obs, filename);
  313 +
  314 + // Remove temporary file
  315 + std::remove(filename.c_str());
  316 +
  317 + } // endif: had a CTA observation
  318 +
  319 + } // endfor: looped over all observations
  320 +
  321 + // If more than a single observation has been handled then make sure that
  322 + // an XML file will be used for storage
  323 + if (n_observations > 1) {
  324 + m_use_xml = true;
  325 + }
  326 +
  327 + // Write observation(s) into logger
  328 + if (logTerse()) {
  329 + log << std::endl;
  330 + log.header1("Observations after selection");
  331 + log << m_obs << std::endl;
  332 + }
  333 +
  334 + // Return
  335 + return;
  336 +}
  337 +
  338 +
  339 +/***********************************************************************//**
  340 + * @brief Save the selected event list(s)
  341 + *
  342 + * This method saves the selected event list(s) into FITS file(s). There are
  343 + * two modes, depending on the m_use_xml flag.
  344 + *
  345 + * If m_use_xml is true, all selected event list(s) will be saved into FITS
  346 + * files, where the output filenames are constructued from the input
  347 + * filenames by prepending the m_prefix string to name. Any path information
  348 + * will be stripped form the input name, hence event files will be written
  349 + * into the local working directory (unless some path information is present
  350 + * in the prefix). In addition, an XML file will be created that gathers
  351 + * the filename information for the selected event list(s). If an XML file
  352 + * was present on input, all metadata information will be copied from this
  353 + * input file.
  354 + *
  355 + * If m_use_xml is false, the selected event list will be saved into a FITS
  356 + * file.
  357 + ***************************************************************************/
  358 +void ctcubemask::save(void)
  359 +{
  360 + // Write header
  361 + if (logTerse()) {
  362 + log << std::endl;
  363 + if (m_obs.size() > 1) {
  364 + log.header1("Save observations");
  365 + }
  366 + else {
  367 + log.header1("Save observation");
  368 + }
  369 + }
  370 +
  371 + // Case A: Save event file(s) and XML metadata information
  372 + if (m_use_xml) {
  373 + save_xml();
  374 + }
  375 +
  376 + // Case B: Save event file as FITS file
  377 + else {
  378 + save_fits();
  379 + }
  380 +
  381 + // Return
  382 + return;
  383 +}
  384 +
  385 +
  386 +/***********************************************************************//**
  387 + * @brief Get application parameters
  388 + *
  389 + * Get all task parameters from parameter file or (if required) by querying
  390 + * the user. Times are assumed to be in the native CTA MJD format.
  391 + *
  392 + * This method also loads observations if no observations are yet allocated.
  393 + * Observations are either loaded from a single CTA even list, or from a
  394 + * XML file using the metadata information that is stored in that file.
  395 + ***************************************************************************/
  396 +void ctcubemask::get_parameters(void)
  397 +{
  398 + // If there are no observations in container then add a single CTA
  399 + // observation using the parameters from the parameter file
  400 + if (m_obs.size() == 0) {
  401 +
  402 + // Get CTA event list file name
  403 + m_infile = (*this)["infile"].filename();
  404 +
  405 + // Allocate CTA observation
  406 + GCTAObservation obs;
  407 +
  408 + // Try first to open as FITS file
  409 + try {
  410 +
  411 + // Load event list in CTA observation
  412 + obs.load_binned(m_infile);
  413 +
  414 + // Append CTA observation to container
  415 + m_obs.append(obs);
  416 +
  417 + // Signal that no XML file should be used for storage
  418 + m_use_xml = false;
  419 +
  420 + }
  421 +
  422 + // ... otherwise try to open as XML file
  423 + catch (GException::fits_open_error &e) {
  424 +
  425 + // Load observations from XML file
  426 + m_obs.load(m_infile);
  427 +
  428 + // Signal that XML file should be used for storage
  429 + m_use_xml = true;
  430 +
  431 + }
  432 +
  433 + } // endif: there was no observation in the container
  434 +
  435 + // Get parameters
  436 + m_regfile = (*this)["regfile"].filename();
  437 + m_usepnt = (*this)["usepnt"].boolean();
  438 + if (!m_usepnt) {
  439 + m_ra = (*this)["ra"].real();
  440 + m_dec = (*this)["dec"].real();
  441 + }
  442 + m_rad = (*this)["rad"].real();
  443 + m_emin = (*this)["emin"].real();
  444 + m_emax = (*this)["emax"].real();
  445 +
  446 + // Optionally read ahead parameters so that they get correctly
  447 + // dumped into the log file
  448 + if (m_read_ahead) {
  449 + m_outfile = (*this)["outfile"].filename();
  450 + }
  451 +
  452 + // Return
  453 + return;
  454 +}
  455 +
  456 +
  457 +/***********************************************************************//**
  458 + * @brief Select events
  459 + *
  460 + * @param[in] obs CTA observation.
  461 + * @param[in] filename File name.
  462 + *
  463 + * Select events from a FITS file by making use of the selection possibility
  464 + * of the cfitsio library on loading a file. A selection string is created
  465 + * from the specified criteria that is appended to the filename so that
  466 + * cfitsio will automatically filter the event data. This selection string
  467 + * is then applied when opening the FITS file. The opened FITS file is then
  468 + * saved into a temporary file which is the loaded into the actual CTA
  469 + * observation, overwriting the old CTA observation. The ROI, GTI and EBounds
  470 + * of the CTA event list are then set accordingly to the specified selection.
  471 + * Finally, the temporary file created during this process is removed.
  472 + *
  473 + * Good Time Intervals of the observation will be limited to the time
  474 + * interval [m_tmin, m_tmax]. If m_tmin=m_tmax=0, no time selection is
  475 + * performed.
  476 + *
  477 + * @todo Use INDEF instead of 0.0 for pointing as RA/DEC selection
  478 + ***************************************************************************/
  479 +void ctcubemask::apply_mask(GCTAObservation* obs, const std::string& filename)
  480 +{
  481 + std::cout<< "start appply filter" << std::endl;
  482 + // Open FITS file
  483 + GFits file(filename);
  484 + GCTAEventCube* cube = (GCTAEventCube*) obs->events();
  485 + GEbounds ebounds = cube->ebounds();
  486 + GSkymap map = cube->map();
  487 +
  488 + // Apply energy selection
  489 + int e_idx1 = 0;
  490 + int e_idx2 = 0;
  491 + int n_ebin = ebounds.size();
  492 +
  493 + // Loop over ebounds to find the first energy band
  494 + for ( int i = 0 ; i < ebounds.size() ; i++){
  495 + double emin = ebounds.emin(i).TeV();
  496 + if ( emin > m_emin ){
  497 + e_idx1 = i;
  498 + break;
  499 + }
  500 + }
  501 +
  502 + std::cout << "finish looping over ebounds to find 1st" << std::endl;
  503 + std::cout << "e_ix1 " << e_idx1 << std::endl;
  504 +
  505 + // Loop over ebounds to find the last energy band
  506 + for ( int i = ebounds.size()-1 ; i < ebounds.size() ; i--){
  507 + double emin = ebounds.emin(i).TeV();
  508 + if ( emin < m_emax ){
  509 + e_idx2 = i-1;
  510 + break;
  511 + }
  512 + }
  513 + std::cout << "finish looping over ebounds to find 1st" << std::endl;
  514 + std::cout << "e_ix2 " << e_idx2 << std::endl;
  515 +
  516 + // Set all pixels outside the desired energy bands negative
  517 + int npix = map.npix();
  518 + for ( int i = 0 ; i < e_idx1 ; i++){
  519 + for ( int pixel = 0 ; pixel < npix ; pixel++){
  520 + map(pixel,i) = -1;
  521 + }
  522 + }
  523 + for ( int i = e_idx2 ; i < ebounds.size() ; i++){
  524 + for ( int pixel = 0 ; pixel < npix ; pixel++){
  525 + map(pixel,i) = -1;
  526 + }
  527 + }
  528 + std::cout << "finish energy band filter" << std::endl;
  529 +
  530 +
  531 +
  532 + // Set all pixels inside the desired energy bands
  533 + // but outside ROI or inside exlusion regions negative
  534 + GSkyRegions regs = GSkyRegions(m_regfile);
  535 + GSkyRegionCircle roi = GSkyRegionCircle(m_ra, m_dec, m_rad);
  536 + for ( int i = e_idx1 ; i <= e_idx2 ; i++){
  537 + for ( int pixel = 0 ; pixel < npix ; pixel++){
  538 + GSkyDir dir = map.inx2dir(pixel);
  539 + if ( roi.contains(dir) == false ||
  540 + regs.contains(dir) == true){
  541 + map(pixel,i) = -1;
  542 + }
  543 + }
  544 + }
  545 + std::cout << "finish roi selection" << std::endl;
  546 +
  547 + GCTAEventCube newcube( map, ebounds, obs->events()->gti());
  548 + obs->events(newcube);
  549 +
  550 + // Log selected FITS file
  551 + if (logExplicit()) {
  552 + log << std::endl;
  553 + log.header1("FITS file content after filter");
  554 + log << file << std::endl;
  555 + }
  556 +
  557 + // Check if we have an IMAGE HDU
  558 + if (!file.contains("IMAGE")) {
  559 + std::string message = "No \"IMAGE\" extension found in FITS file. ";
  560 + throw GException::app_error(G_APPLY_MASK, message);
  561 + }
  562 +
  563 +
  564 + // Get CTA event list pointer
  565 + GCTAEventList* list =
  566 + static_cast<GCTAEventList*>(const_cast<GEvents*>(obs->events()));
  567 +
  568 + // Return
  569 + return;
  570 +}
  571 +
  572 +
  573 +/*==========================================================================
  574 + = =
  575 + = Private methods =
  576 + = =
  577 + ==========================================================================*/
  578 +
  579 +/***********************************************************************//**
  580 + * @brief Initialise class members
  581 + ***************************************************************************/
  582 +void ctcubemask::init_members(void)
  583 +{
  584 + // Initialise parameters
  585 + m_infile.clear();
  586 + m_regfile.clear();
  587 + m_outfile.clear();
  588 + m_prefix.clear();
  589 + m_usepnt = false;
  590 + m_ra = 0.0;
  591 + m_dec = 0.0;
  592 + m_rad = 0.0;
  593 + m_emin = 0.0;
  594 + m_emax = 0.0;
  595 +
  596 + // Initialise protected members
  597 + m_obs.clear();
  598 + m_infiles.clear();
  599 + m_use_xml = false;
  600 + m_read_ahead = false;
  601 +
  602 + // Set logger properties
  603 + log.date(true);
  604 +
  605 + // Return
  606 + return;
  607 +}
  608 +
  609 +
  610 +/***********************************************************************//**
  611 + * @brief Copy class members
  612 + *
  613 + * @param[in] app Application.
  614 + ***************************************************************************/
  615 +void ctcubemask::copy_members(const ctcubemask& app)
  616 +{
  617 + // Copy parameters
  618 + m_infile = app.m_infile;
  619 + m_regfile = app.m_regfile;
  620 + m_outfile = app.m_outfile;
  621 + m_prefix = app.m_prefix;
  622 + m_usepnt = app.m_usepnt;
  623 + m_ra = app.m_ra;
  624 + m_dec = app.m_dec;
  625 + m_rad = app.m_rad;
  626 + m_emin = app.m_emin;
  627 + m_emax = app.m_emax;
  628 +
  629 + // Copy protected members
  630 + m_obs = app.m_obs;
  631 + m_infiles = app.m_infiles;
  632 + m_use_xml = app.m_use_xml;
  633 + m_read_ahead = app.m_read_ahead;
  634 +
  635 + // Return
  636 + return;
  637 +}
  638 +
  639 +
  640 +/***********************************************************************//**
  641 + * @brief Delete class members
  642 + ***************************************************************************/
  643 +void ctcubemask::free_members(void)
  644 +{
  645 + // Write separator into logger
  646 + if (logTerse()) {
  647 + log << std::endl;
  648 + }
  649 +
  650 + // Return
  651 + return;
  652 +}
  653 +
  654 +
  655 +/***********************************************************************//**
  656 + * @brief Check input filename
  657 + *
  658 + * @param[in] filename File name.
  659 + *
  660 + * This method checks if the input FITS file is correct.
  661 + ***************************************************************************/
  662 +std::string ctcubemask::check_infile(const std::string& filename) const
  663 +{
  664 + // Initialise message string
  665 + std::string message = "";
  666 +
  667 + // Open FITS file
  668 + GFits file(filename);
  669 +
  670 + // Check for IMAGE HDU
  671 + GFitsImage* image = NULL;
  672 + try {
  673 +
  674 + // Get pointer to FITS image
  675 + image = file.image("IMAGE");
  676 + }
  677 + catch (GException::fits_hdu_not_found& e) {
  678 + message = "No \"IMAGE\" extension found in input file \""
  679 + + m_outfile + "\".";
  680 + }
  681 +
  682 + // Return
  683 + return message;
  684 +}
  685 +
  686 +
  687 +/***********************************************************************//**
  688 + * @brief Set output file name.
  689 + *
  690 + * @param[in] filename Input file name.
  691 + *
  692 + * This converts an input filename into an output filename by prepending a
  693 + * prefix to the input filename. Any path will be stripped from the input
  694 + * filename.
  695 + ***************************************************************************/
  696 +std::string ctcubemask::set_outfile_name(const std::string& filename) const
  697 +{
  698 + // Split input filename into path elements
  699 + std::vector<std::string> elements = gammalib::split(filename, "/");
  700 +
  701 + // The last path element is the filename
  702 + std::string outname = m_prefix + elements[elements.size()-1];
  703 +
  704 + // Return output filename
  705 + return outname;
  706 +}
  707 +
  708 +
  709 +/***********************************************************************//**
  710 + * @brief Save counts map in FITS format.
  711 + *
  712 + * Save the counts map as a FITS file. The filename of the FITS file is
  713 + * specified by the m_outfile member.
  714 + ***************************************************************************/
  715 +void ctcubemask::save_fits(void)
  716 +{
  717 + // Get output filename
  718 + m_outfile = (*this)["outfile"].filename();
  719 +
  720 + // Get CTA observation from observation container
  721 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
  722 +
  723 + // Save event list
  724 + save_counts_map(obs, m_outfile);
  725 +
  726 + // Return
  727 + return;
  728 +}
  729 +
  730 +
  731 +/***********************************************************************//**
  732 + * @brief Save counts map(s) in XML format.
  733 + *
  734 + * Save the counts map(s) into FITS files and write the file path information
  735 + * into a XML file. The filename of the XML file is specified by the
  736 + * m_outfile member, the filename(s) of the counts map(s) are built by
  737 + * prepending the prefix given by the m_prefix member to the input counts
  738 + * map(s) filenames. Any path present in the input filename will be stripped,
  739 + * i.e. the counts map(s) will be written in the local working directory
  740 + * (unless a path is specified in the m_prefix member).
  741 + ***************************************************************************/
  742 +void ctcubemask::save_xml(void)
  743 +{
  744 + // Get output filename and prefix
  745 + m_outfile = (*this)["outfile"].filename();
  746 + m_prefix = (*this)["prefix"].string();
  747 +
  748 + // Issue warning if output filename has no .xml suffix
  749 + std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
  750 + if (suffix != ".xml") {
  751 + log << "*** WARNING: Name of observation definition output file \""+
  752 + m_outfile+"\"" << std::endl;
  753 + log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
  754 + log << "*** WARNING: This is not an error, but might be misleading."
  755 + " It is recommended" << std::endl;
  756 + log << "*** WARNING: to use the suffix \".xml\" for observation"
  757 + " definition files." << std::endl;
  758 + }
  759 +
  760 + // Loop over all observation in the container
  761 + for (int i = 0; i < m_obs.size(); ++i) {
  762 +
  763 + // Get CTA observation
  764 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  765 +
  766 + // Handle only CTA observations
  767 + if (obs != NULL) {
  768 +
  769 + // Set event output file name
  770 + std::string outfile = set_outfile_name(m_infiles[i]);
  771 +
  772 + // Store output file name in observation
  773 + obs->eventfile(outfile);
  774 +
  775 + // Save event list
  776 + save_counts_map(obs, outfile);
  777 +
  778 + } // endif: observation was a CTA observations
  779 +
  780 + } // endfor: looped over observations
  781 +
  782 + // Save observations in XML file
  783 + m_obs.save(m_outfile);
  784 +
  785 + // Return
  786 + return;
  787 +}
  788 +
  789 +
  790 +
  791 +/***********************************************************************//**
  792 + * @brief Save a single counts map into a FITS file
  793 + *
  794 + * @param[in] obs Pointer to CTA observation.
  795 + * @param[in] outfile Output file name.
  796 + *
  797 + * This method saves a single counts map into a FITS file. The method does
  798 + * nothing if the observation pointer is not valid.
  799 + ***************************************************************************/
  800 +void ctcubemask::save_counts_map(const GCTAObservation* obs,
  801 + const std::string& outfile) const
  802 +{
  803 + // Save only if observation is valid
  804 + if (obs != NULL) {
  805 +
  806 + // Save observation into FITS file
  807 + obs->save(outfile, clobber());
  808 +
  809 + } // endif: observation was valid
  810 +
  811 + // Return
  812 + return;
  813 +}
  814 +
... ...
src/ctcubemask/ctcubemask.hpp 0 → 100644
  View file @b8e6842
  1 +/***************************************************************************
  2 + * ctcubemask - CTA cube filter 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 ctcubemask.hpp
  23 + * @brief CTA mask cube tool definition
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +#ifndef CTCUBEMASK_HPP
  28 +#define CTCUBEMASK_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GammaLib.hpp"
  34 +#include "GCTALib.hpp"
  35 +
  36 +/* __Definitions _________________________________________________________ */
  37 +#define CTCUBEMASK_NAME "ctcubemask"
  38 +#define CTCUBEMASK_VERSION "00-05-00"
  39 +
  40 +
  41 +/***********************************************************************//**
  42 + * @class ctcubemask
  43 + *
  44 + * @brief CTA mask cube tool interface defintion.
  45 + ***************************************************************************/
  46 +class ctcubemask : public GApplication {
  47 +public:
  48 + // Constructors and destructors
  49 + ctcubemask(void);
  50 + explicit ctcubemask(GObservations obs);
  51 + ctcubemask(int argc, char *argv[]);
  52 + ctcubemask(const ctcubemask& app);
  53 + virtual ~ctcubemask(void);
  54 +
  55 + // Operators
  56 + ctcubemask& operator= (const ctcubemask& app);
  57 +
  58 + // Methods
  59 + void clear(void);
  60 + void execute(void);
  61 + void run(void);
  62 + void save(void);
  63 + GObservations& obs(void) { return m_obs; }
  64 + void get_parameters(void);
  65 + void apply_mask(GCTAObservation* obs, const std::string& filename);
  66 +
  67 +protected:
  68 + // Protected methods
  69 + void init_members(void);
  70 + void copy_members(const ctcubemask& app);
  71 + void free_members(void);
  72 + std::string check_infile(const std::string& filename) const;
  73 + std::string set_outfile_name(const std::string& filename) const;
  74 + void save_fits(void);
  75 + void save_xml(void);
  76 + void save_counts_map(const GCTAObservation* obs,
  77 + const std::string& outfile) const;
  78 +
  79 + // User parameters
  80 + std::string m_infile; //!< Input event list or XML file
  81 + std::string m_regfile; //!< ds9 region file
  82 + std::string m_outfile; //!< Output event list or XML file
  83 + std::string m_prefix; //!< Prefix for multiple counts maps
  84 + bool m_usepnt; //!< Use pointing instead of RA/DEC parameters
  85 + double m_ra; //!< RA of ROI centre
  86 + double m_dec; //!< DEC of ROI centre
  87 + double m_rad; //!< ROI radius
  88 + double m_emin; //!< Lower energy
  89 + double m_emax; //!< Upper energy
  90 +
  91 + // Protected members
  92 + GObservations m_obs; //!< Observations container
  93 + std::vector<std::string> m_infiles; //!< Input event filenames
  94 + bool m_use_xml; //!< Use XML file instead of FITS file
  95 + bool m_read_ahead; //!< Read ahead parameters
  96 +};
  97 +
  98 +#endif /* CTCUBEMASK_HPP */
... ...
src/ctcubemask/ctcubemask.par 0 → 100644
  View file @b8e6842
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctfilter.par #
  8 +# Version: 00-04-00 #
  9 +# #
  10 +# Author: Chia-Chun Lu #
  11 +# c.lu@mpi-hd.mpg.de #
  12 +# IRAP #
  13 +# #
  14 +# Purpose: Parameter file for the ctselect application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# File information
  19 +#=================
  20 +infile, f, a, "test/data/crab_events.fits.gz",,,"Input event list or observation definition file"
  21 +regfile, f, a, "exlusion.reg",,,"exclusion region file in ds9 format"
  22 +outfile, f, a, "filtered_cube.fits",,,"Output event list or observation definition file"
  23 +prefix, s, h, "filtered_",,,"Prefix for counts maps in observation definition file"
  24 +#
  25 +# Selection parameters
  26 +#=====================
  27 +usepnt, b, h, no,,, "Use pointing instead of RA/DEC parameters?"
  28 +ra, r, a, 83.63, 0,360, "RA for ROI centre (degrees)"
  29 +dec, r, a, 22.01, -90, 90, "Dec for ROI centre (degrees)"
  30 +rad, r, a, 3.0, 0,180, "Radius of ROI (degrees)"
  31 +emin, r, a, 0.1, 0, , "Lower energy limit (TeV)"
  32 +emax, r, a, 100.0, 0, , "Upper energy limit (TeV)"
  33 +#
  34 +# Standard parameters
  35 +#====================
  36 +chatter, i, h, 2,0,4, "Chattiness of output"
  37 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  38 +debug, b, h, no,,, "Debugging mode activated"
  39 +mode, s, h, "ql",,, "Mode of automatic parameters"
... ...
src/ctcubemask/main.cpp 0 → 100644
  View file @b8e6842
  1 +/***************************************************************************
  2 + * ctcubemask - CTA mask cube tool main code *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 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 ctcubemask/main.cpp
  23 + * @brief CTA mask cube tool main code
  24 + * @author Chia-Chun Lu
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctcubemask.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point
  37 + *
  38 + * @param[in] argc Number of arguments
  39 + * @param[in] argv Arguments
  40 + ***************************************************************************/
  41 +int main (int argc, char *argv[])
  42 +{
  43 + // Initialise return code
  44 + int rc = 1;
  45 +
  46 + // Create instance of application
  47 + ctcubemask application(argc, argv);
  48 +
  49 + // Run application
  50 + try {
  51 + // Execute application
  52 + application.execute();
  53 +
  54 + // Signal success
  55 + rc = 0;
  56 + }
  57 + catch (std::exception &e) {
  58 +
  59 + // Extract error message
  60 + std::string message = e.what();
  61 + std::string signal = "*** ERROR encounterted in the execution of"
  62 + " ctcubemask. Run aborted ...";
  63 +
  64 + // Write error in logger
  65 + application.log << message << std::endl;
  66 + application.log << signal << std::endl;
  67 +
  68 + // Write error on standard output
  69 + std::cout << message << std::endl;
  70 + std::cout << signal << std::endl;
  71 +
  72 + } // endcatch: catched any application error
  73 +
  74 + // Return
  75 + return rc;
  76 +}
... ...